text:
##########################################
## Extract class and product from a Symbol
##########################################

function get_class(output_var::Symbol)
    s = string(output_var)
    if occursin("pseudo", s)
        :pseudo
    elseif occursin("obs", s)
        :obs
    elseif occursin("state", s)
        :states
    elseif occursin("stdshock", s)
        :stdshocks
    elseif occursin("shock", s)
        :shocks
    else
        error("Invalid output_var: " * s)
    end
end

function get_product(output_var::Symbol)
    s = string(output_var)
    if occursin("bddhistforecast4q", s)
        :bddhistforecast4q
    elseif occursin("histforecast4q", s)
        :histforecast4q
    elseif occursin("bddhistforecast", s)
        :bddhistforecast
    elseif occursin("histforecast", s)
        :histforecast
    elseif occursin("hist4q", s)
        :hist4q
    elseif occursin("histut", s)
        :histut
    elseif occursin("hist", s)
        :hist
    elseif occursin("bddforecast4q", s)
        :bddforecast4q
    elseif occursin("forecast4q", s)
        :forecast4q
    elseif occursin("bddforecastut", s)
        :bddforecastut
    elseif occursin("forecastut", s)
        :forecastut
    elseif occursin("bddforecast", s)
        :bddforecast
    elseif occursin("forecast", s)
        :forecast
    elseif occursin("shockdec", s)
        :shockdec
    elseif occursin("dettrend", s)
        :dettrend
    elseif occursin("trend", s)
        :trend
    elseif occursin("irf", s)
        :irf
    elseif occursin("decomp", s)
        :decomp
    else
        error("Invalid output_var: " * s)
    end
end

function get_class_longname(class::Symbol)
    if class == :pseudo
        :pseudoobservable
    elseif class == :obs
        :observable
    elseif class == :states
        :state
    elseif class in [:shocks, :stdshocks]
        :shock
    end
end

#####################################
## Preparing data for transformations
#####################################

"""
```
resize_population_forecast(population_forecast::DataFrame, nperiods::Int,
    mnemonic::Symbol)
```

Extends or shrinks the population forecasts to be `nperiods` in
length. If `population_forecast` must be extended, the last value is
simply repeated as many times as necessary.
"""
function resize_population_forecast(population_forecast::DataFrame, nperiods::Int,
                                    mnemonic::Symbol)

    nperiods_act = size(population_forecast, 1)
    if nperiods_act >= nperiods
        # Keep first nperiods of population forecast
        return population_forecast[1:nperiods, :]

    else
       # Extend forecast by n_filler_periods
        n_filler_periods = nperiods - nperiods_act

        last_provided = population_forecast[end, :date]

        first_filler = iterate_quarters(last_provided, 1)
        last_filler  = iterate_quarters(last_provided, n_filler_periods)
        filler_dates = quarter_range(first_filler, last_filler)
        filler = DataFrame(date = filler_dates)
        filler[!,mnemonic] .= population_forecast[end, mnemonic]
        return vcat(population_forecast, filler)
    end
end

"""
```
parse_transform(t::Symbol)
```

Parse the module name out of a Symbol to recover the transform associated with
an observable or pseudoobservable. Returns a function.
"""
parse_transform(t::Symbol) = eval(Symbol(split(string(t),".")[end]))

"""
```
load_population_growth(m; verbose = :low)
```

Returns `DataFrame`s of growth rates for HP-filtered (unless
`hpfilter_population(m) == false`) population data and forecast.
"""
function load_population_growth(m::AbstractDSGEModel; verbose::Symbol = :low)

    data_verbose = verbose == :none ? :none : :low

    # Parse population mnemonic
    mnemonic = parse_population_mnemonic(m)[1]

    if isnull(mnemonic)
        return DataFrame(), DataFrame()

    else
        # Read in unfiltered series
        vint = data_vintage(m)
        data_file = inpath(m, "raw", "population_data_levels_$vint.csv")
        unfiltered_data = read_population_data(data_file; verbose = data_verbose)
        unfiltered_forecast = if use_population_forecast(m)
            forecast_file = inpath(m, "raw", "population_forecast_$vint.csv")
            read_population_forecast(forecast_file, get(mnemonic); verbose = data_verbose)
        else
            DataFrame()
        end
        if (isempty(unfiltered_forecast) && use_population_forecast(m))
            error("Must set use_population_forecast setting to false if no population_forecast_$vint.csv file is provided as data. This file should be manually added to the saveroot in the raw sub-directory if the user wants to use a population forecast.")
        end

        # HP filter if necessary
        data, forecast = transform_population_data(unfiltered_data, unfiltered_forecast, get(mnemonic);
                                                   use_hpfilter = hpfilter_population(m),
                                                   verbose = :none)
        if hpfilter_population(m)
            data_mnemonic = :dlfiltered_population_recorded
            forecast_mnemonic = :dlfiltered_population_forecast
        else
            data_mnemonic = :dlpopulation_recorded
            forecast_mnemonic = :dlpopulation_forecast
        end

        # Prepare output variables
        data = data[!, [:date, data_mnemonic]]
        rename!(data, data_mnemonic => :population_growth)

        if use_population_forecast(m)
            forecast = forecast[!, [:date, forecast_mnemonic]]
            rename!(forecast, forecast_mnemonic => :population_growth)
        else
            forecast = DataFrame()
        end

        return data, forecast
    end
end

"""
```
get_population_series(mnemonic, population_data, population_forecast,
    start_date, end_date)
```

Returns the population series between `start_date` and `end_date`. If `end_date`
is after the last period in `population_forecast`, the last forecasted
population is repeated until `end_date` using
`resize_population_forecast`. Returns a `Vector{Float64}`.
"""
function get_population_series(mnemonic::Symbol, population_data::DataFrame,
                               population_forecast::DataFrame, start_date::Dates.Date,
                               end_date::Dates.Date)

    last_historical_date = population_data[end, :date]

    # If no population forecast, use last period of population data instead
    if isempty(population_forecast)
        population_forecast = DataFrame(deepcopy(population_data[end, :]))
        last_forecast_date = iterate_quarters(last_historical_date, 1)
        population_forecast[1, :date] = last_forecast_date
    else
        last_forecast_date = population_forecast[end, :date]
    end

    # Calculate number of periods that are in the future and extend
    # population forecast by the right number of periods
    if end_date > last_forecast_date
        n_fcast_periods = subtract_quarters(end_date, last_historical_date)
        population_forecast = resize_population_forecast(population_forecast, n_fcast_periods,
                                                         mnemonic)
    end

    population_insample = if start_date <= population_data[end, :date]

        padding = if start_date < population_data[1, :date]
            # Start date is before population data; compute number of NaNs to prepend
            @warn "Start date $start_date is before population data begins: prepending NaNs"
            n_nans = subtract_quarters(population_data[1, :date], start_date)

            DataFrame(date = quarter_range(start_date, iterate_quarters(population_data[1,:date], -1)))
        else
            DataFrame()
        end

        unpadded_data = if population_data[1, :date] < end_date < population_data[end, :date]
            # Dates entirely in past
            population_data[start_date .<= population_data[!, :date] .<= end_date, :]
        else
            # Dates span past and forecast
            data  = population_data[start_date .<= population_data[!, :date], :]
            fcast = population_forecast[population_forecast[!, :date] .<= end_date, :]
            vcat(data, fcast)
        end

        padding, unpadded_data = reconcile_column_names(padding, unpadded_data)
        padded_data = vcat(padding, unpadded_data)
        # na2nan!(padded_data)
        padded_data

    elseif population_forecast[1, :date] <= start_date <= population_forecast[end, :date]
        # Dates entirely in forecast
        population_forecast[start_date .<= population_forecast[!, :date] .<= end_date, :]
    else
        # start_date comes after population_forecast[end, :date]
        error("Start date $start_date comes after population forecast ends")
    end

    return population_insample[!, mnemonic]
end

"""
```
get_mb_population_series(product, mnemonic, population_data, population_forecast, date_list)
```

Returns the appropriate population series for the `product`.
"""
function get_mb_population_series(product::Symbol, population_data::DataFrame,
                                  population_forecast::DataFrame,
                                  date_list::Vector{Dates.Date})

    if product == :irf
        # Return empty vector for IRFs, which don't correspond to real dates
        return Float64[]
    elseif isempty(population_data) && isempty(population_forecast)
        # Return empty vector if population data and forecasts not provided
        return Float64[]
    else
        start_date = if product in [:hist4q, :forecast4q, :bddforecast4q]
            iterate_quarters(date_list[1], -3)
        elseif product in [:histut, :hist, :forecastut, :forecast, :bddforecastut, :bddforecast,
                           :shockdec, :dettrend, :trend]
            date_list[1]
        else
            error("Invalid product: $product")
        end
        end_date = date_list[end]

        pop_growth = get_population_series(:population_growth, population_data, population_forecast,
                                           start_date, end_date)
        pop_growth[ismissing.(pop_growth)] .= NaN
        return convert(Vector{Float64}, pop_growth)
    end
end

##############
## Other utils
##############

"""
```
get_mb_metadata(m, input_type, cond_type, output_var, forecast_output_file = "";
    forecast_string = "")
```

Returns the `metadata` dictionary from `read_forecast_metadata`, as well as
`mb_metadata`, the dictionary that we will save to the means and bands file.
"""
function get_mb_metadata(m::AbstractDSGEModel, input_type::Symbol, cond_type::Symbol,
                         output_var::Symbol, forecast_output_file::String = "";
                         forecast_string::String = "")

    if isempty(forecast_output_file)
        forecast_output_file = get_meansbands_input_file(m, input_type, cond_type, output_var,
                                                          forecast_string = forecast_string)
    end

    class   = get_class(output_var)
    product = get_product(output_var)
    metadata = JLD2.jldopen(forecast_output_file, "r") do jld
        read_forecast_metadata(jld)
    end

    class_long       = get_class_longname(class)
    variable_indices = metadata[Symbol(class_long, "_indices")]
    date_indices     = product == :irf ? Dict{Date, Int}() : metadata[:date_indices]

    # Make sure date lists are valid. This is vacuously true for and IRFs, which
    # are not time-dependent and hence have empty `date_indices`.
    date_list          = collect(keys(date_indices))   # unsorted array of actual dates
    date_indices_order = collect(values(date_indices)) # unsorted array of date indices
    check_consistent_order(date_list, date_indices_order)

    mb_metadata = Dict{Symbol,Any}(
                   :para            => input_type,
                   :cond_type       => cond_type,
                   :product         => product,
                   :class           => class,
                   :indices         => sort(variable_indices, by = x -> variable_indices[x]),
                   :date_inds       => sort(date_indices, by = x -> date_indices[x]),
                   :forecast_string => forecast_string)

    if product in [:shockdec, :irf]
        mb_metadata[:shock_indices] = metadata[:shock_indices]
    end

    return mb_metadata
end

"""
```
get_y0_index(m::AbstractDSGEModel, product::Symbol)
```

Returns the index of the period *before* the start of saved results
for a particular forecast-step product. For regular forecasts, this is
the last period for which we have full-quarter data. For running 4q
forecasts, this corresponds to 4 periods before the first forecast
period. For shockdecs, this is one period before
`index_shockdec_start(m)`, and for histories, it is the last period of
the presample.
"""
function get_y0_index(m::AbstractDSGEModel, product::Symbol)
    if product in [:forecastut, :forecast, :bddforecastut, :bddforecast]
        return index_forecast_start(m) - 1
    elseif product in [:forecast4q, :bddforecast4q]
        # We subtract 4 because there is 1 transform that actually
        # needs us to go 4 periods. Later, we can use y0_index + 1
        # to index out the data we need for all the other forecasts.
        return index_forecast_start(m) - 4
    elseif product in [:shockdec, :dettrend, :trend]
        return n_presample_periods(m) + index_shockdec_start(m) - 1
    elseif product in [:hist, :histut, :hist4q]
        return index_mainsample_start(m) - 1
    elseif product == :irf
        return -1
    else
        error("get_y0_index not implemented for product = $product")
    end
end

function get_yt_index(m::AbstractDSGEModel, product::Symbol)
    if product in [:forecastut, :forecast, :bddforecastut, :bddforecast]
        return index_forecast_start(m) - 1
    elseif product in [:forecast4q, :bddforecast4q]
        return index_forecast_start(m) - 1
    elseif product in [:shockdec, :dettrend, :trend, :hist, :histut, :hist4q, :irf]
        return -1
        # Note: below, I don't immediately know what we want for yt (end) and I don't think I need to implmented
#=    elseif product in [:shockdec, :dettrend, :trend]
        return n_presample_periods(m) + index_shockdec_start(m) - 1
    elseif product in [:hist, :histut, :hist4q]
        return index_mainsample_start(m) - 1
    elseif product == :irf
        return -1=#
    else
        error("get_yt_index not implemented for product = $product")
    end
end

"""
```
check_consistent_order(l1, l2)
```

Checks to make sure that l1 and l2 are ordered consistently.
"""
function check_consistent_order(l1, l2)
    @assert length(l1) == length(l2)

    # cache old pairs in a Dict
    original_pairs = Dict{eltype(l1), eltype(l2)}()
    for (i,item) in enumerate(l1)
        original_pairs[item] = l2[i]
    end

    # sort
    l1_sorted = sort(l1)
    l2_sorted = sort(l2)

    # make sure each pair has same pair as before
    for i in 1:length(l1)
        @assert original_pairs[l1_sorted[i]] == l2_sorted[i] "Lists not consistently ordered"
    end

    return true
end


text:
function features(board::Board)
    black_stones = zeros(Int, board.num_rows, board.num_cols)
    white_stones = zeros(Int, board.num_rows, board.num_cols)
    for r in 1:board.num_rows
        for c in 1:board.num_cols
            if board[r, c].color === black
                black_stones[r, c] = 1
            elseif board[r, c].color === white
                white_stones[r, c] = 1
            end
        end
    end
    return Int[vec(black_stones); vec(white_stones)]
end

function features(game_state::GameState)
    feats = Int[]
    state = game_state
    for lag in 1:8
        append!(feats, features(state.board))
        if state.prev !== nothing
            state = state.prev
        end
    end
    push!(feats, Int(state.next_player == black))
    feats
end

ext: jl


text:
export RAC

using ..Frames: AbstractFrame, ICRF

struct RAC <: AbstractFrame end

function Rotation(::RAC, ::ICRF, r, v)
    normal = cross(r, v)
    normal = normal / norm(normal)
    tangential = v / norm(v)
    orthogonal = cross(v, normal)
    orthogonal = orthogonal / norm(orthogonal)
    m = hcat(orthogonal, tangential, normal)
    Rotation{RAC(), icrf}(m)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE
using Test, Tensors

let parameters = MemoryParameterState(E = 200.0e3,
                                      nu = 0.3,
                                      R0 = 100.0,
                                      Kn = 20.0,
                                      nn = 3.0,
                                      C1 = 10000.0,
                                      D1 = 100.0,
                                      C2 = 50000.0,
                                      D2 = 1000.0,
                                      Q0 = 100.0,
                                      QM = 500.0,
                                      mu = 100.0,
                                      b = 30.0,
                                      eta = 0.5,
                                      m = 0.5,
                                      pt = 0.0,
                                      xi = 0.3),
    mat = Memory(parameters=parameters),

    tostrain(tens::Symm2) = copy(tovoigt(tens; offdiagscale=2.0)),
    tostress(tens::Symm2) = copy(tovoigt(tens)),

    n_cycles = 30,
    points_per_cycle = 40,
    t = range(0.0; stop=Float64(n_cycles), length=n_cycles * points_per_cycle + 1),
    dtime = t[end] / (length(t) - 1),

    # We initialize these manually to automatically get the correct type.
    times = [copy(mat.drivers.time)],
    stresses = [tostress(mat.variables.stress)],
    strains = [tostrain(mat.drivers.strain)],
    plastic_strains = [tostrain(mat.variables.plastic_strain)],
    cumeqs = [copy(mat.variables.cumeq)],
    qs = [copy(mat.variables.q)],
    Rs = [copy(mat.variables.R)],
    zetas = [tostrain(mat.variables.zeta)]

    function snapshot!()
        push!(times, mat.drivers.time)
        push!(stresses, tostress(mat.variables.stress))
        push!(strains, tostrain(mat.drivers.strain))
        push!(plastic_strains, tostrain(mat.variables.plastic_strain))
        push!(cumeqs, copy(mat.variables.cumeq))
        push!(qs, copy(mat.variables.q))
        push!(Rs, copy(mat.variables.R))
        push!(zetas, tostrain(mat.variables.zeta))
    end

    # Amplitude 1
    ea = 0.003
    strains11 = ea * sin.(2*pi*t)
    for dstrain11 in diff(strains11)
        uniaxial_increment!(mat, dstrain11, dtime)
        update_material!(mat)
        snapshot!()
    end
    R1 = copy(Rs[end])

    # Amplitude 2
    ea = 0.005
    strains11 = ea * sin.(2*pi*t)
    for dstrain11 in diff(strains11)
        uniaxial_increment!(mat, dstrain11, dtime)
        update_material!(mat)
        snapshot!()
    end
    R2 = copy(Rs[end])

    # Amplitude 3
    ea = 0.007
    strains11 = ea * sin.(2*pi*t)
    for dstrain11 in diff(strains11)
        uniaxial_increment!(mat, dstrain11, dtime)
        update_material!(mat)
        snapshot!()
    end
    R3 = copy(Rs[end])

    # Amplitude 4 - evanescence
    ea = 0.003
    strains11 = ea * sin.(2*pi*t)
    for _ in 1:3
        for dstrain11 in diff(strains11)
            uniaxial_increment!(mat, dstrain11, dtime)
            update_material!(mat)
            snapshot!()
        end
    end

    R4 = copy(Rs[end])

    @test R2 > R1
    @test R3 > R2
    @test R4 < R3
    @test isapprox(R1, R4; atol=1.0)
end


text:
#include("MetabolicFluxAnalysis.jl")
using MetabolicFluxAnalysis
reload("MetabolicFluxAnalysis")

@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# write your own tests here
@test 1 == 1

println(pwd())
path_project     = splitdir( pwd() )
path_to_datasets = joinpath( path_project[1] , "datasets")
println(path_to_datasets)

if(false)

    @testset "MFA.Basics_01" begin include("basic_tests_01.jl") end
    @testset "MFA.Basics_02" begin include("basic_tests_02.jl") end
    @testset "Fluxi.Basics_01" begin include("fluxi_tests_01.jl") end
    @testset "MFA.Polyhedra_01" begin include("tests_polyhedra.jl") end
    @testset "MFA.Sampling_01" begin include("test_samp.jl") end
    @testset "MFA.ImportCobra" begin include("test_cobra_import.jl") end
    @testset "Fluxi.AntoniewiczWorkflow" begin include("antoniewicz_toy_network_test.jl") end
    @testset "Fluxi.WXTransform_01" begin include("test_wx_transform.jl") end
    @testset "Fluxi.CMFA_Net_01" begin include("cmfa_network_test.jl") end
    @testset "Fluxi.ToyNetwork_01" begin include("fluxi_tests_01.jl") end
end

#@testset "Fluxi.AntoniewiczWorkflow2" begin include("antoniewicz_toy_network_test_02.jl") end
#@testset "Fluxi.FWSimTests_01" begin include("fw_sim_tests_01.jl") end
#@testset "MFA.ImportJSON" begin include("test_json_io.jl") end
@testset "Fluxi.OCL_01" begin include("test_opencl_backend.jl") end

# PARAMETRIATION FAILS SOMETIMES.. (Parametrization of polytope..) should try to make this work reliably..
@testset "MFA.ImportFTBL" begin include("test_ftbl_import.jl") end

# @testset "Fluxi.OpenCL_01" begin include("kernel_compilation_01.jl") end # does not work, not needed anymore..

#include("basic_tests_01.jl")
#include("basic_tests_02.jl")
#include("fluxi_tests_01.jl")
#include("antoniewicz_toy_network_test.jl")


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

# This benchmark should demonstrate the issue with Time To SafePoint when
# running some threads that allocate and some threads in heavy math code
# that doesn't offer a safepoint.

function work(i,j)
    out = 1
    for x in 1:i
       for y in 1:j
        out = out * x * y
       end
    end
    out
end

mutable struct ListNode
  key::Int64
  next::ListNode
  ListNode() = new()
  ListNode(x)= new(x)
  ListNode(x,y) = new(x,y);
end

function list(n=128)
    start::ListNode = ListNode(1)
    current::ListNode = start
    for i = 2:(n*1024^2)
        current = ListNode(i,current)
    end
    return current.key
end

function allocate_in_background()
   list(32)
end

function testhelper()
   y = 1024
   if (Threads.threadid() == 1)
      x = 256*1024
      y = work(x,x)
   else
      allocate_in_background()
   end
   y
end

function test()
  n = Threads.nthreads() * 2
  Threads.@threads for i in 1:n
      testhelper()
   end
end

@gctime test()

max_stars_repo_name: vtjnash/GCBenchmarks

max_issues_repo_name: vtjnash/GCBenchmarks

max_forks_repo_name: vtjnash/GCBenchmarks


text:
function sanity_checks_basic( X, x, S, h, minx, maxx )
    jassert( 1 ≤ S ≤ length( X ), "polynomial order $S must be at least one and no greater than the sample size")
    jassert( 0.0 < h < Inf, "bad bandwidth $h")
    jassert( minx ≤ minimum( X ), "data points less than $minx" )
    jassert( maxx ≥ maximum( X ), "data points greater than $maxx" )
    jassert( minx ≤ minimum( x ), "evaluation points less than $minx" )
    jassert( maxx ≥ maximum( x ), "evaluation points greater than $maxx" )
    return nothing
end

function sanity_checks_g( g, S, z, zl, zr )
    for t ∈ eachindex( z ), j ∈ 1:S
        jassert( g( -z[t], j, z[t], zr[t] ) == 0.0, "function g used is not zero at lower boundary" )
        jassert( g( zr[t], j, z[t], zr[t] ) == 0.0, "function g used is not zero at upper boundary" )
    end
    return nothing
end


text:
module CauchyBorn

using JuLIP, StaticArrays, Calculus

import JuLIPMaterials.CLE: elastic_moduli

using JuLIPMaterials: ee, Vec3, Mat3

abstract type Wcb end

export Wcb

Wcb(at::AbstractAtoms) = Wcb(at, calculator(at))

"""
Construct using
```
W = Wcb(at)
W = Wcb(at, calc)
```

## KW Arguments

* `precompute=true`: This (possibly) pre-computes several partial derivatives
of W in the reference state. To prevent this precomputation use the kw-argument
`precompute=false`.
* `normalise = :volume`: this normalises W against the volume of the unit cell.
Alternatively, to normalise against the number of atoms use `normalise = :atoms`.
Or, simply pass in a scalar `normalise = v0`.

## Usage

`W` can be used to evaluate the following quantities:

### For a Simple Lattice (unit cell contains 1 atom)

* `W(F)` : with `F` a 3 x 3 matrix is the energy / undeformed volume;
* `grad(W, F)` : with `F` a 3 x 3 matrix is the first Piola-Kirchhoff stress at
`F`, i.e., the jacobian of `W(F)` w.r.t. `F`
* `div_grad(W,F)` : finite-difference implementation of div ∂W(F)
"""
function Wcb(at::AbstractAtoms, calc::AbstractCalculator;
             normalise=:volume, T = Float64, kwargs...)  # switch to AbstractAtoms{T}
   @assert length(at) == 1
   set_calculator!(at, calc)
   # compute the normalisation factor; volume or what?
   if normalise == :volume
      v0 = det(cell(at))
   elseif normalise == :atoms
      v0 = T(length(at))
   elseif normalise isa Number
      v0 = T(normalise)
   else
      error("Wcb: unrecognised kwarg `normalise = $normalise`")
   end
   # simple lattice case
   if length(at) == 1
      return Wcb1(at, v0; kwargs...)
   elseif length(at) == 2 && length(unique(chemical_symbols(at))) == 1
      return Wcb2(at, v0; kwargs...)
   else
      error("""`Wcb`: so far, only single species 1-lattice and 2-lattice
      have been implemented. If you need a more general case, please file
      an issue at https://github.com/cortner/JuLIPMaterials.jl""")
   end
end

set_rel_defm!(W::Wcb, F) = set_defm!(W.at, F * W.C0')

# ================ Simple Lattice Cauchy--Born  ==============

"""
`struct Wcb1` : Simple Lattice Cauchy--Born potential, see documentation
for `Wcb`
"""
struct Wcb1{TA, TC, T} <: Wcb
   calc::TC
   at::TA        # unit cell
   C0::Mat3{T}   # original cell matrix (cell(at))
   v0::T         # volume of original cell
end

Wcb1(at::AbstractAtoms, v0) = Wcb1(calculator(at), at, Mat3(cell(at)), v0)

(W::Wcb1)(args...) = evaluate(W, args...)

evaluate(W::Wcb1, F) = energy( W.calc, set_rel_defm!(W, F) ) / W.v0

grad(W::Wcb1, F) = (- virial( W.calc, set_rel_defm!(W, F) ) * inv(F)') / W.v0

div_grad(W::Wcb1, F, x::Vec3{T}; h = 1e-5) where T =
      sum( (grad(W, F(x+h*ee(i)))[:,i] - grad(W, F(x-h*ee(i)))[:,i]) / (2*h)
            for i = 1:3 )

# ============= Single Species 2-Lattice =================

"""
`struct Wcb2` : single-species 2-lattice Lattice Cauchy--Born potential,
see documentation for `Wcb`
"""
struct Wcb2{T, TA}
   at::TA        # unit cell
   C0::Mat3{T}   # original cell matrix (cell(at))
   v0::T         # volume of original cell
   p0::Vec3{T}
   p1::Vec3{T}
   dpdpW::Mat3{T}
   dpdpW_inv::Mat3{T}
end

function Wcb2(at::AbstractAtoms; precompute = true)
   @assert length(at) == 2 && length(unique(chemical_symbols(at))) == 1
   if precompute
      dpdpW = dpdpWcb2(at)
      dpdpW_inv = inv(dpdpW)
   else
      dpdpW = zero(Mat3{T})
      dpdpW_inv = zero(Mat3{T})
   end
   X = positions(at)
   return Wcb2(at, Mat3(cell(at)), det(cell(at)), X[1], X[2], dpdpW, dpdpW_inv)
end

set_F_and_p!(W::Wcb2, F, p) = set_positions!(set_rel_defm!(W, F), [W.p0, W.p0+p])

evaluate(W::Wcb2, F, p) = energy( set_F_and_p!(W, F, p) )

grad_p(W::Wcb2, F, p) = - forces( set_F_and_p!(W, F, p) )[2]

# grad_F(W::Wcb2, F, p) = 0  # TODO

function dpdpWcb(W::Wcb2{T}) where T
   J = Calculus.jacobian(p -> grad_p(W, one(SMatrix{3,3,T}), p))
   dpdpW = J(W.p1)
   return Mat3{T}(0.5 * (dpdpW + dpdpW'))
end

# # TODO: replace with get_shift
# function (W::WcbQuad)(F)
#     # TODO this is fishy - why is the initial position not reset?
#     p0 = positions(W.at)[2]
#     p1 = p0 - W.dpdpW_inv * DpWcb(F, p0, W.at, W.calc)
#     p2 = p1 - W.dpdpW_inv * DpWcb(F, p1, W.at, W.calc)
#     return p2
# end
#
#
# function DW_infp(W::WcbQuad, F)
#     p = W(F)
#     at = W.at
#     set_defm!(at, F)
#     X = positions(at)
#     X[2] = p
#     set_positions!(at, X)
#     return stress(W.calc, at)
# end
#
#
# function elastic_moduli(W::WcbQuad)
#    F0 = W.F0 |> Matrix
#    Ih = eye(3)
#    h = eps()^(1/3)
#    C = zeros(3,3,3,3)
#    for i = 1:3, a = 1:3
#       Ih[i,a] += h
#       Sp = DW_infp(W, Ih * F0)
#       Ih[i,a] -= 2*h
#       Sm = DW_infp(W, Ih * F0)
#       C[i, a, :, :] = (Sp - Sm) / (2*h)
#       Ih[i,a] += h
#    end
#    # symmetrise it - major symmetries C_{iajb} = C_{jbia}
#    for i = 1:3, a = 1:3, j=1:3, b=1:3
#       t = 0.5 * (C[i,a,j,b] + C[j,b,i,a])
#       C[i,a,j,b] = t
#       C[j,b,i,a] = t
#    end
#    # minor symmetries - C_{iajb} = C_{iabj}
#    for i = 1:3, a = 1:3, j=1:3, b=1:3
#       t = 0.5 * (C[i,a,j,b] + C[i,a,b,j])
#       C[i,a,j,b] = t
#       C[i,a,b,j] = t
#    end
#    return C
# end

# =============== GENERAL MULTI-LATTICE FUNCTIONALITY ================

# TODO
end


text:
module FilterUtils

export lowpass_butterworth

using DSP

function lowpass_butterworth(f0, dt, n::Int=4)
    filter = digitalfilter(Lowpass(f0;fs=1/dt), Butterworth(n))
    return data -> filtfilt(filter, data)
end

end # module


text:
using Test, LinearAlgebra, SparseArrays, Random
using QDLDL
rng = Random.MersenneTwister(2701)

@testset "inertia checks" begin

  m = 20
  n = 30

  A = random_psd(n) + Diagonal(rand(n))
  B = sprandn(m,n,0.2)
  C = -Diagonal(rand(m))
  M = [A B' ; B C]

  b = randn(m+n)
  F = qdldl(M)

  @test positive_inertia(F) == n

end

nothing


text:
const METIS_NOPTIONS = 40
## Return codes
const METIS_OK           = 1  # normal return
const METIS_ERROR_INPUT  = -2 # erroneous inputs and/or options
const METIS_ERROR_MEMORY = -3 # insufficient memory
const METIS_ERROR        = -4 # Other errors
## Operation type codes
const METIS_OP_PMETIS = 0
const METIS_OP_KMETIS = 1
const METIS_OP_OMETIS = 2
## (1-based) positions in options vector
const METIS_OPTION_PTYPE     = 1
const METIS_OPTION_OBJTYPE   = 2
const METIS_OPTION_CTYPE     = 3
const METIS_OPTION_IPTYPE    = 4
const METIS_OPTION_RTYPE     = 5
const METIS_OPTION_DBGLVL    = 6
const METIS_OPTION_NITER     = 7
const METIS_OPTION_NCUTS     = 8
const METIS_OPTION_SEED      = 9
const METIS_OPTION_NO2HOP    = 10
const METIS_OPTION_MINCONN   = 11
const METIS_OPTION_CONTIG    = 12
const METIS_OPTION_COMPRESS  = 13
const METIS_OPTION_CCORDER   = 14
const METIS_OPTION_PFACTOR   = 15
const METIS_OPTION_NSEPS     = 16
const METIS_OPTION_UFACTOR   = 17
const METIS_OPTION_NUMBERING = 18
const METIS_OPTION_HELP      = 19
const METIS_OPTION_TPWGTS    = 20
const METIS_OPTION_NCOMMON   = 21
const METIS_OPTION_NOOUTPUT  = 22
const METIS_OPTION_BALANCE   = 23
const METIS_OPTION_GTYPE     = 24
const METIS_OPTION_UBVEC     = 25
## Partitioning Schemes ## Warning that types need to be adapted
const METIS_PTYPE_RB   = 0
const METIS_PTYPE_KWAY = 1
## Graph types for meshes
const METIS_GTYPE_DUAL  = 0
const METIS_GTYPE_NODAL = 1
## Coarsening Schemes
const METIS_CTYPE_RM   = 0
const METIS_CTYPE_SHEM = 1
## Initial partitioning schemes
const METIS_IPTYPE_GROW    = 0
const METIS_IPTYPE_RANDOM  = 1
const METIS_IPTYPE_EDGE    = 2
const METIS_IPTYPE_NODE    = 3
const METIS_IPTYPE_METISRB = 4
## Refinement schemes
const METIS_RTYPE_FM        = 0
const METIS_RTYPE_GREEDY    = 1
const METIS_RTYPE_SEP2SIDED = 2
const METIS_RTYPE_SEP1SIDED = 3
## Debug levels (bit positions)
const METIS_DBG_INFO       = 1 # Shows various diagnostic messages
const METIS_DBG_TIME       = 2 # Perform timing analysis
const METIS_DBG_COARSEN    = 4 # Show the coarsening progress
const METIS_DBG_REFINE     = 8 # Show the refinement progress
const METIS_DBG_IPART      = 16 # Show info on initial partitioning
const METIS_DBG_MOVEINFO   = 32 # Show info on vertex moves during refinement
const METIS_DBG_SEPINFO    = 64 # Show info on vertex moves during sep refinement
const METIS_DBG_CONNINFO   = 128 # Show info on minimization of subdomain connectivity
const METIS_DBG_CONTIGINFO = 256 # Show info on elimination of connected components
const METIS_DBG_MEMORY     = 2048 # Show info related to wspace allocation
## Types of objectives
const METIS_OBJTYPE_CUT  = 0
const METIS_OBJTYPE_VOL  = 1
const METIS_OBJTYPE_NODE = 2


text:
try
    ts = time_ns()
    let
        test = ByteStreamTestHarness("write-end-pop", 15)

        execute(test, Write("cat"))

        execute(test, InputEnded(false))
        execute(test, BufferEmpty(false))
        execute(test, Eof(false))
        execute(test, BytesRead(0))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(12))
        execute(test, BufferSize(3))
        execute(test, Peek("cat"))

        execute(test, EndInput())

        execute(test, InputEnded(true))
        execute(test, BufferEmpty(false))
        execute(test, Eof(false))
        execute(test, BytesRead(0))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(12))
        execute(test, BufferSize(3))
        execute(test, Peek("cat"))

        execute(test, Pop(3))

        execute(test, InputEnded(true))
        execute(test, BufferEmpty(true))
        execute(test, Eof(true))
        execute(test, BytesRead(3))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))
    end
    let
        test = ByteStreamTestHarness("write-pop-end", 15)

        execute(test, Write("cat"))

        execute(test, InputEnded(false))
        execute(test, BufferEmpty(false))
        execute(test, Eof(false))
        execute(test, BytesRead(0))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(12))
        execute(test, BufferSize(3))
        execute(test, Peek("cat"))

        execute(test, Pop(3))

        execute(test, InputEnded(false))
        execute(test, BufferEmpty(true))
        execute(test, Eof(false))
        execute(test, BytesRead(3))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))

        execute(test, EndInput())

        execute(test, InputEnded(true))
        execute(test, BufferEmpty(true))
        execute(test, Eof(true))
        execute(test, BytesRead(3))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))
    end
    let
        test = ByteStreamTestHarness("write-pop2-end", 15)

        execute(test, Write("cat"))

        execute(test, InputEnded(false))
        execute(test, BufferEmpty(false))
        execute(test, Eof(false))
        execute(test, BytesRead(0))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(12))
        execute(test, BufferSize(3))
        execute(test, Peek("cat"))

        execute(test, Pop(1))

        execute(test, InputEnded(false))
        execute(test, BufferEmpty(false))
        execute(test, Eof(false))
        execute(test, BytesRead(1))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(13))
        execute(test, BufferSize(2))
        execute(test, Peek("at"))

        execute(test, Pop(2))

        execute(test, InputEnded(false))
        execute(test, BufferEmpty(true))
        execute(test, Eof(false))
        execute(test, BytesRead(3))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))

        execute(test, EndInput())

        execute(test, InputEnded(true))
        execute(test, BufferEmpty(true))
        execute(test, Eof(true))
        execute(test, BytesRead(3))
        execute(test, BytesWritten(3))
        execute(test, RemainingCapacity(15))
        execute(test, BufferSize(0))
    end
    tc = (Int(time_ns() - ts)) / 1e9
    @info "byte_stream_one_write passed......... Time: $tc"

catch e
    println(stderr, e)
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
abstract type TensorProductBasis <: Function end

"""
Constructs a Chebyshev basis of the form [T_{0}(x), T_{1}(x), ..., T_{n-1}(x)] where T_j(.) is the j-th Chebyshev polynomial of the first kind.
```julia
ChebyshevBasis(n)
```
Arguments:
- `n`: number of terms in the polynomial expansion.
"""
struct ChebyshevBasis <: TensorProductBasis
    n::Int
end

function (basis::ChebyshevBasis)(x)
    return map(j -> cos(j*acos(x)), 1:basis.n)
end

"""
Constructs a sine basis of the form [sin(x), sin(2*x), ..., sin(n*x)].
```julia
SinBasis(n)
```
Arguments:
- `n`: number of terms in the sine expansion.
"""
struct SinBasis <: TensorProductBasis
    n::Int
end

function (basis::SinBasis)(x)
    return map(j -> sin(j*x), 1:basis.n)
end

"""
Constructs a cosine basis of the form [cos(x), cos(2*x), ..., cos(n*x)].
```julia
CosBasis(n)
```
Arguments:
- `n`: number of terms in the cosine expansion.
"""
struct CosBasis <: TensorProductBasis
    n::Int
end

function (basis::CosBasis)(x)
    return map(j -> cos(j*x), 1:basis.n)
end

#auxiliary function
function fourier(i::Int, x::Real)
    return iseven(i) ? cos(i*x/2) : sin(i*x/2)
end

"""
Constructs a Fourier basis of the form F_j(x) = j is even ? cos((j÷2)*x) : sin((j÷2)*x) => [F_0(x), F_1(x), ..., F_n(x)].
```julia
FourierBasis(n)
```
Arguments:
- `n`: number of terms in the Fourier expansion.
"""
struct FourierBasis <: TensorProductBasis
    n::Int
end

function (basis::FourierBasis)(x)
    return map(j -> fourier(j,x), 1:basis.n)
end

#auxiliary function
##Source: https://github.com/ranocha/PolynomialBases.jl/blob/master/src/legendre.jl
function legendre_poly(x, p::Integer)
    a::typeof(x) = one(x)
    b::typeof(x) = x

    if p <= 0
        return a
    elseif p == 1
        return b
    end

    for j in 2:p
        a, b = b, ((2j-1)*x*b - (j-1)*a) / j
    end

    b
end

"""
Constructs a Legendre basis of the form [P_{0}(x), P_{1}(x), ..., P_{n-1}(x)] where P_j(.) is the j-th Legendre polynomial.
```julia
LegendreBasis(n)
```
Arguments:
- `n`: number of terms in the polynomial expansion.
"""
struct LegendreBasis <: TensorProductBasis
    n::Int
end

function (basis::LegendreBasis)(x)
    f = k -> legendre_poly(x,k-1)
    return map(f, 1:basis.n)
end

"""
Constructs a Polynomial basis of the form [1, x, ..., x^(n-1)].
```julia
PolynomialBasis(n)
```
Arguments:
- `n`: number of terms in the polynomial expansion.
"""
struct PolynomialBasis <: TensorProductBasis
    n::Int
end

function (basis::PolynomialBasis)(x)
    return [evalpoly(x, (I+zeros(basis.n,basis.n))[k,:]) for k in 1:basis.n]
end


text:
module base64
export NoPadding, StdPadding
export CorruptInput, InvalidAlphabet, InvalidPadCharacter
export Encoding, NewEncoding, WithPadding, Strict
export StdEncoding, URLEncoding, RawURLEncoding, RawStdEncoding
export EncodedLen, Encode!, EncodeToString
export DecodedLen, Decode!, DecodeString

const StdPadding = 0x3d
const NoPadding = 0x00

const std_codec = [UInt8(x) for x in ['A':'Z'; 'a':'z'; '0':'9'; '+'; '/']]
const url_codec = [UInt8(x) for x in ['A':'Z'; 'a':'z'; '0':'9'; '-'; '_']]

struct CorruptInput <: Exception end
struct InvalidAlphabet <: Exception end
struct InvalidPadCharacter <: Exception end

struct Encoding
	encode::Array{UInt8, 1}
	decodeMap::Array{UInt8, 1}
	padChar::UInt8
  strict::Bool

  Encoding(encode::Array{UInt8, 1}, decodeMap::Array{UInt8, 1}, padChar::UInt8, strict::Bool) = new(encode, decodeMap, padChar, strict)
  Encoding(;encode::Array{UInt8, 1}, decodeMap::Array{UInt8, 1} = fill(0xff, 256), padChar::UInt8 = StdPadding, strict::Bool=false) = new(encode, decodeMap, padChar, strict)
end

@inline function assemble_64(n1::UInt8, n2::UInt8, n3::UInt8, n4::UInt8, n5::UInt8, n6::UInt8, n7::UInt8, n8::UInt8) :: Union{Nothing, UInt64}
  n1|n2|n3|n4|n5|n6|n7|n8 != 0xff && return UInt64(n1)<<58 | UInt64(n2)<<52 | UInt64(n3)<<46 | UInt64(n4)<<40 | UInt64(n5)<<34 | UInt64(n6)<<28 | UInt64(n7)<<22 | UInt64(n8)<<16
end

@inline function assemble_32(n1::UInt8, n2::UInt8, n3::UInt8, n4::UInt8) :: Union{Nothing, UInt32}
  n1|n2|n3|n4 != 0xff && return UInt32(n1)<<26 | UInt32(n2)<<20 | UInt32(n3)<<14 | UInt32(n4)<<8
end

@inline function put_64!(dst::AbstractArray{UInt8, 1}, val::UInt64) :: Nothing
  length(dst) < 8 && throw(BoundsError("destinary array out of bounds"))
  dst[1] = val>>56 & 0xff
  dst[2] = val>>48 & 0xff
  dst[3] = val>>40 & 0xff
  dst[4] = val>>32 & 0xff
  dst[5] = val>>24 & 0xff
  dst[6] = val>>16 & 0xff
  dst[7] = val>>8 & 0xff
	dst[8] = val & 0xff
	nothing
end

@inline function put_32!(dst::AbstractArray{UInt8, 1}, val::UInt32) :: Nothing
  length(dst) < 4 && throw(BoundsError("destinary array out of bounds"))
  dst[1] = val>>24 & 0xff
  dst[2] = val>>16 & 0xff
  dst[3] = val>>8 & 0xff
	dst[4] = val & 0xff
	nothing
end

function NewEncoding(alphabet::Array{UInt8, 1}) :: Encoding
  length(alphabet) != 64 && throw(InvalidAlphabet())
  check = in(alphabet)
  (check(0x0d) || check(0x0a)) && throw(InvalidAlphabet())
  encoder = Encoding(encode=alphabet)
  for i in 1:64
    @inbounds encoder.decodeMap[alphabet[i]] = i & 0xff - 1
  end
  return encoder
end
NewEncoding(alphabet::String) = NewEncoding(Array{UInt8, 1}(alphabet))

function WithPadding(enc::Encoding, pad_char::UInt8) :: Encoding
  (pad_char == 0x0d || pad_char == 0x0a || pad_char > 0x3d || in(pad_char, enc.encode)) && throw(InvalidPadCharacter())
  Encoding(enc.encode, enc.decodeMap, pad_char, enc.strict)
end
WithPadding(enc::Encoding, pad_char::Char) = WithPadding(enc, UInt8(pad_char))

function Strict(enc::Encoding) :: Encoding
  Encoding(enc.encode, enc.decodeMap, enc.padChar, true)
end

const StdEncoding = NewEncoding(std_codec)
const URLEncoding = NewEncoding(url_codec)

const RawStdEncoding = WithPadding(StdEncoding, NoPadding)
const RawURLEncoding = WithPadding(URLEncoding, NoPadding)

@inline function EncodedLen(enc::Encoding, n::Int) :: Int
  if enc.padChar == NoPadding
    return convert(Int, floor((n*8 + 5) / 6))
  end
  convert(Int, floor((n+2) / 3) * 4)
end

function Encode!(enc::Encoding, dst::Array{UInt8, 1}, src::Array{UInt8, 1}) :: Nothing
  isempty(src) && return nothing
  di, si, n = 0, 0, convert(Int, floor(length(src) / 3) * 3)
  while si < n
    v = UInt(src[si+1])<<16 | UInt(src[si+2])<<8 | UInt(src[si+3])
    @inbounds dst[di+1] = enc.encode[v>>18&0x3f + 1]
    @inbounds dst[di+2] = enc.encode[v>>12&0x3f + 1]
    @inbounds dst[di+3] = enc.encode[v>>6&0x3f + 1]
    @inbounds dst[di+4] = enc.encode[v&0x3f + 1]

    si += 3
    di += 4
  end

  rem = length(src) - si
  rem == 0 && return nothing

  v = UInt(src[si+1]) << 16
  if rem == 2
    v |= UInt(src[si+2]) << 8
  end

  @inbounds dst[di+1] = enc.encode[v>>18&0x3f + 1]
  @inbounds dst[di+2] = enc.encode[v>>12&0x3f + 1]

  if rem == 2
    @inbounds dst[di+3] = enc.encode[v>>6&0x3f + 1]
		if enc.padChar != NoPadding
			@inbounds dst[di+4] = enc.padChar
    end
  elseif rem == 1
    if enc.padChar != NoPadding
      @inbounds dst[di+3] = enc.padChar
      @inbounds dst[di+4] = enc.padChar
    end
  end
  nothing
end

@inline function EncodeToString(enc::Encoding, src::Array{UInt8, 1}) :: String
  buf = fill(enc.padChar, EncodedLen(enc, length(src)))
  Encode!(enc, buf, src)
  String(buf)
end
EncodeToString(enc::Encoding, src::String) = EncodeToString(enc, Array{UInt8, 1}(src))

@inline function DecodedLen(enc::Encoding, n::Int) :: Int
  if enc.padChar == NoPadding
    return convert(Int, floor(n * 6 / 8))
  end
  convert(Int, floor(n / 4 * 3))
end

function decode_quantum!(enc::Encoding, dst::AbstractArray{UInt8, 1}, src::Array{UInt8, 1}, si::Int) :: Tuple{Int, Int}
	dbuf = zeros(4)
	dlen = 4

  j = 1
  while j < length(dbuf)
    if length(src) == si
      if j == 1
        return si, 0
      elseif j == 2 || enc.padChar != NoPadding
        throw(CorruptInput())
      end
      dlen = j - 1
      break
    end
    i = src[si + 1]
    si += 1
    o = enc.decodeMap[i]
    if o != 0xff
      dbuf[j] = o
      j += 1
      continue
    end
    if i == 0x0d || i == 0x0a
      j -= 1
      continue
    end
    if i != enc.padChar
      throw(CorruptInput())
    end
    if j == 1 || j == 2
      throw(CorruptInput())
    elseif j == 3
      while si < length(src) && (src[si+1] == 0x0d || src[si+1] == 0x0a)
        si += 1
      end
      (si == length(src) || src[si+1] != enc.padChar) && throw(CorruptInput())
      si += 1
    end

    while si < length(src) && (src[si+1] == 0x0d || src[si+1] == 0x0a)
      si += 1
    end

    (si < length(src)) && throw(CorruptInput())
    dlen = j - 1
    break
  end

  v = UInt(dbuf[1])<<18 | UInt(dbuf[2])<<12 | UInt(dbuf[3])<<6 | UInt(dbuf[4])
  dbuf[3], dbuf[2], dbuf[1] = v>>0 & 0xff, v>>8 & 0xff, v>>16 & 0xff
  if dlen == 4
    @inbounds dst[3] = dbuf[3]
    dbuf[3] = 0
  end
  if dlen == 4 || dlen == 3
    @inbounds dst[2] = dbuf[2]
    if enc.strict && dbuf[3] != 0
      throw(CorruptInput())
    end
    dbuf[2] = 0
  end
  if dlen == 4 || dlen == 3 || dlen == 2
    @inbounds dst[1] = dbuf[1]
    if enc.strict && (dbuf[2] != 0 || dbuf[3] != 0)
      throw(CorruptInput())
    end
  end
  return si, dlen-1
end

function Decode!(enc::Encoding, dst::Array{UInt8, 1}, src::Array{UInt8, 1}) :: Int
  isempty(src) && return nothing
  si, n = 0, 0
  while length(src) - si >= 8 && length(dst) - n >= 8
    dn = assemble_64(
      enc.decodeMap[src[si+1]],
			enc.decodeMap[src[si+2]],
			enc.decodeMap[src[si+3]],
			enc.decodeMap[src[si+4]],
			enc.decodeMap[src[si+5]],
			enc.decodeMap[src[si+6]],
			enc.decodeMap[src[si+7]],
			enc.decodeMap[src[si+8]],
    )
    if !isnothing(dn)
      put_64!(view(dst, n+1:length(dst)), dn)
      n += 6
      si += 8
    else
      si, ni = decode_quantum!(enc, view(dst, n+1:length(dst)), src, si)
      n += ni
    end
  end

  while length(src)-si >= 4 && length(dst)-n >= 4
    dn = assemble_32(
      enc.decodeMap[src[si+1]],
			enc.decodeMap[src[si+2]],
			enc.decodeMap[src[si+3]],
			enc.decodeMap[src[si+4]],
    )
    if !isnothing(dn)
      put_32!(view(dst, n+1:length(dst)), dn)
      n += 3
      si += 4
    else
      si, ni = decode_quantum!(enc, view(dst, n+1:length(dst)), src, si)
      n += ni
    end
  end

  if si < length(src)
    si, ni = decode_quantum!(enc, view(dst, n+1:length(dst)), src, si)
    n += ni
  end
  return n
end

@inline function DecodeString(enc::Encoding, src::Array{UInt8, 1}) :: String
  buf = fill(enc.padChar, DecodedLen(enc, length(src)))
  n = Decode!(enc, buf, src)
  String(buf[1:n])
end
DecodeString(enc::Encoding, src::String) = DecodeString(enc, Array{UInt8, 1}(src))

end


text:
import JSON3
using Sockets

function main(text)
  jobj = JSON3.read(text)
  coordinates = jobj["coordinates"]
  len = length(coordinates)
  x = y = z = 0

  for coord in coordinates
    x += coord["x"]
    y += coord["y"]
    z += coord["z"]
  end

  println(x / len)
  println(y / len)
  println(z / len)
end

function test(text)
  x = @timed main(text)
  println("Elapsed: $(x[2]), Allocated: $(x[3]), GC Time: $(x[4])")
end

text = open("/tmp/1.json") do file
  read(file, String)
end

for i in 1:4
  test(text)
end

function notify(msg)
  try
    socket = connect("localhost", 9001)
    write(socket, msg)
    close(socket)
  catch
    # standalone usage
  end
end

notify("Julia JSON3\t$(getpid())")

test(text)

notify("stop")

max_stars_repo_name: insigniaapple/benchmarks

max_issues_repo_name: insigniaapple/benchmarks

max_forks_repo_name: insigniaapple/benchmarks


text:
using DataStructures, Compat
using Compat: String
if VERSION >= v"0.5.0-dev+7720"
    using Base.Test
else
    using BaseTestNext
    const Test = BaseTestNext
end

@testset "CircularDeque" begin
    D = CircularDeque{Int}(5)
    @test eltype(D) == Int
    @test capacity(D) == 5
    @test length(D) == 0
    @test isempty(D)
    @test_throws BoundsError front(D)
    @test_throws BoundsError back(D)
    push!(D, 1)
    @test front(D) === back(D) === 1
    push!(D, 2)
    @test front(D) === 1
    @test back(D)  === 2
    @test length(D) == 2
    for i = 3:5
        push!(D, i)
    end
    @test_throws BoundsError push!(D, 6)
    @test shift!(D) === 1
    @test front(D) === 2
    @test back(D) === 5
    io = IOBuffer()
    print(io, Int)
    intstr = String(take!(io))
    print(io, D)
    @test String(take!(io)) == "CircularDeque{$intstr}([2,3,4,5])"
    push!(D, 6)
    @test front(D) === 2
    @test back(D) === 6
    @test pop!(D) === 6
    @test front(D) === 2
    @test back(D) === 5
    unshift!(D, 7)
    @test front(D) === 7
    @test back(D) === 5
    @test_throws BoundsError unshift!(D, 8)
    @test shift!(D) === 7
    @test shift!(D) === 2
    @test pop!(D) === 5
    @test shift!(D) === 3
    @test pop!(D) === 4
    @test_throws BoundsError pop!(D)
    @test_throws BoundsError shift!(D)
    @test isempty(D)
    push!(D, 10)
    @test !isempty(D)
    empty!(D)
    @test isempty(D)
    @test_throws BoundsError front(D)
    push!(D, 20)
    @test front(D) == back(D) == 20
    empty!(D)
    for i = 1:5
        push!(D, i)
    end
    @test shift!(D) == 1
    push!(D, 6)
    for i = 2:6
        @test back(D) === 6
        @test D[1] === i
        @test D[7-i] === 6
        @test shift!(D) === i
    end

    # Test that unshift! works on an empty deque, and that front/back give the right answer
    D = CircularDeque{Int}(5)
    unshift!(D, 30)
    @test front(D) == back(D) == 30
    empty!(D)
    unshift!(D, 40)
    @test front(D) == back(D) == 40

    # Test iteration over loop
    D = CircularDeque{Int}(5)
    for i in 1:5 push!(D, i) end
    @test collect([i for i in D]) == collect(1:5)
end

nothing


text:
int main() {
      double x ;
      x = 2 * 3.14 ;
      return 0 ;
}

max_stars_repo_name: jwilk/jtc

max_issues_repo_name: jwilk/jtc

max_forks_repo_name: jwilk/jtc


text:
using Pkg

Pkg.develop(PackageSpec(path=pwd()))
Pkg.instantiate()

using Documenter, BED

makedocs(
    format = Documenter.HTML(
        edit_branch = "develop"
    ),
    sitename = "BED.jl",
    modules = [BED],
    pages = [
        "Home" => "index.md",
        "API Reference" => "lib/public.md"
    ],
    authors = replace(join(Pkg.TOML.parsefile("Project.toml")["authors"], ", "), r" <.*?>" => "" ) * ", The BioJulia Organisation, and other contributors."

)
deploydocs(
    repo = "github.com/BioJulia/BED.jl.git",
    devbranch = "develop",
    push_preview = true
)


text:
module PLE

using Revise
using PyPlot
using StaticArrays
using LinearAlgebra
using StatsBase
using Distributions
using Seaborn
using LsqFit

include("MyStructs.jl")
include("helper.jl")
include("traverse.jl")
include("checkcollision.jl")
include("assignobstacles.jl")
include("gettrajectory.jl")
include("visualize.jl")
include("getparticle.jl")
include("postprocess.jl")

set_context("poster")
set_context("poster")

end


text:
module Utils

using CSV, DataFrames
using Dates
using Base.Unicode

export @exported_enum, qstrtodict, pretty, basiccensor, get_if_set
export addsysnotoname, diff_between, mult_dict!
export loadtoframe, age_in_years, isapprox, ≈

function addsysnotoname(names, sysno)::Array{Symbol,1}
   a = Array{Symbol,1}(undef, 0)
   for n in names
      push!(a, Symbol(String(n) * "_$sysno"))
   end
   a
    # Symbol.(String.( names ).*"_sys_$sysno")
end

function get_if_set(key, dict::Dict, default; operation=nothing)
   v = default
   if haskey(dict, key)
      v = dict[key]
      if operation !== nothing
         v = operation( v )
      end
   end
   v
end

import Base.isapprox

"""
true if the keys are the same and all the elements
   compare approx equal, else false
"""
function isapprox( d1::Dict, d2::Dict ) :: Bool
   k1 = Set( collect( keys(d1)))
   k2 = Set( collect( keys(d2)))
   if k1 != k2
      # println(" resturn ne $(k1) !== $(k2)")
      return false
   end
   for k in k1
      if ! (d1[k] ≈ d2[k])
         # println( "return ne $k $(d1[k]) ≈ $(d2[k])")
         return false
      end
   end
   return true
end

const ≈ = isapprox

"""
return a dict with all the elements in common where m2[k]-m1[k] is possible.
"""
function diff_between(m2::Dict, m1::Dict)::Dict
   out = Dict()
   keybs = intersect(keys(m1), keys(m2))
   for k in keybs
      v1 = m1[k]
      v2 = m2[k]
      try
         d = v2 - v1
         out[k] = d
      catch

      end # exception
   end # loop[]
   out
end

function mult_dict!( m :: Dict, n :: Number )
   for (k,v) in m
      try
         m[k] = v*n
      catch
         ;
      end
   end
end

"""
parse an html query string like "sex=male&joe=22&bill=21342&z=12.20"
into a dict. If the value looks like a number, it's parsed into either an Int64 or Float64
"""
function qstrtodict(query_string::AbstractString)::Dict{AbstractString,Any}
   d = Dict{AbstractString,Any}()
   strip(query_string)
   if (query_string == "")
      return d
   end
   as = split(query_string, "&")
   for a in as
      try
         aa = split(a, "=")
         k = aa[1]
         v = aa[2]
         try
            v = parse(Int64, v)
         catch
            try
               v = parse(Float64, v)
            catch
            end
         end
         d[k] = v
      catch

      end
   end
   d
end

"""
returns the string converted to a form suitable to be used as (e.g.) a Symbol,
with leading/trailing blanks removed, forced to lowercase, and with various
characters replaced with '_' (at most '_' in a run).
"""
function basiccensor(s::AbstractString)::AbstractString
   s = strip(lowercase(s))
   s = replace(s, r"[ \-,\t–]" => "_")
   s = replace(s, r"[=\:\)\('’‘]" => "")
   s = replace(s, r"[\";:\.\?\*”“]" => "")
   s = replace(s, r"_$" => "")
   s = replace(s, r"^_" => "")
   s = replace(s, r"^_" => "")
   s = replace(s, r"\/" => "_or_")
   s = replace(s, r"\&" => "_and_")
   s = replace(s, r"\+" => "_plus_")
   s = replace(s, r"_\$+$" => "")
   if occursin(r"^[\d].*", s)
      s = string("v_", s) # leading digit
   end
   s = replace(s, r"__" => "_")
   s = replace(s, r"__" => "_")
   s = replace(s, r"__" => "_") # fixme neater way?
   s = replace(s, r"^_" => "")
   s = replace(s, r"_$" => "")
   return s
end

"""
a_string_or_symbol_like_this => "A String Or Symbol Like This"
"""
function pretty(a)
   s = string(a)
   s = strip(lowercase(s))
   s = replace(s, r"[_]" => " ")
   Unicode.titlecase(s)
end

"""
 macro to define an enum and automatically
 add export statements for its elements
 see: https://discourse.julialang.org/t/export-enum/5396
"""
macro exported_enum(name, args...)
   esc(quote
      @enum($name, $(args...))
      export $name
      for a in $args
         local av = string(a)
         :(export $av)
      end
              # $([:(export $arg) for arg in args]...)
   end)
end

"""
load a file into a dataframe and force all the identifiers into
lower case
"""
function loadtoframe(filename::AbstractString)::DataFrame
    df = CSV.File(filename, delim = '\t') |> DataFrame
    lcnames = Symbol.(lowercase.(string.(names(df))))
    rename!(df, lcnames)
    df
end

"""
Age now for someone with this birthday.
If today is 2020-01-27 then:

   * dob(1958-01-26) = 62
   * dob(1958-01-27) = 62
   * dob(1958-01-28) = 61

And so on.
"""
function age_in_years(
   dob :: Dates.TimeType,
   to_date :: Dates.TimeType = Dates.now() ) :: Integer
   @assert dob <= to_date
   y_to = year(to_date)
   m_to = month(to_date)
   d_to = day(to_date)
   y_dob = year(dob)
   m_dob = month(dob)
   d_dob = day(dob)
   age = y_to - y_dob
   if m_dob > m_to # check if you've not yet had your birthday ..
      age -= 1
   elseif (m_dob == m_to ) && (d_dob > d_to )
      age -= 1
   end
   age
end

end # module


text:
@discord_object struct Attachment
    id::Snowflake
    filename::String
    size::Int
    url::String
    proxy_url::String
    height::Int
    width::Int
end


text:
function degree_based_matching(A::SparseMatrixCSC{T,Int},B::SparseMatrixCSC{T,Int}) where T

	m = size(A,1)
	n = size(B,1)

    noise_factor = .01
	d_A = A*ones(m) .+ noise_factor*rand(m)
	d_B = B*ones(n) .+ noise_factor*rand(n)

	p_A = sort(1:m,by=i->d_A[i]) #breaks ties with some variation
	p_B = sort(1:n,by=j->d_B[j])

	return p_A,p_B

end

function low_rank_matching(U::Array{Float64,2},V::Array{Float64,2})
    n,d1 = size(U)
    m,d2 = size(V)
    @assert d1 == d2

    matching_weights = zeros(d1)
    matchings = Array{Dict{Int,Int},1}(undef,d1)

    for i = 1:d1
        matchings[i],matching_weights[i] = rank_one_matching(U[:,i],V[:,i])
    end

    D = zeros(d1,d2)
    for j in 1:d2
        for i in 1:d1
            w = 0.0
            for (v_i,v_j) in matchings[j]
                edge_w = U[v_i,i]*V[v_j,i]
                if edge_w > 0
                    w += edge_w
                end
            end
            D[i,j] = matching_weights[i]/w
        end
    end

    d_j = maximum(D,dims= 1)
    opt_j = argmin(d_j).I[2]

    return matchings[opt_j], matching_weights[opt_j]

end

function rank_one_matching(u::Array{Float64,1},v::Array{Float64,1})

    #get maximum matching by rearrangement theorem
    u_perm = sortperm(u,rev=true)
    v_perm = sortperm(v,rev=true)

    matching_weight = 0
    Match_mapping = Dict{Int,Int}()
    for (i,j) in zip(u_perm,v_perm)

        #only match positive to positive and negative to negative
        if (u[i] > 0 && v[j] > 0 ) || (u[i] < 0 && v[j] < 0)
            Match_mapping[i] = j
            matching_weight += u[i]*v[j]
        end
    end
    return Match_mapping, matching_weight

end

function search_Krylov_space(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,U::Array{Float64,2},V::Array{Float64,2})

    best_score = -1
    best_i = -1
    best_j = -1
    best_matching = Dict{Int,Int}()

    Triangle_check = Dict{Array{Int,1},Int}()
    A_unique_nnz = length(A.values)
    B_unique_nnz = length(B.values)

    if A_unique_nnz > B_unique_nnz
        for i in 1:A_unique_nnz
            Triangle_check[A.indices[i,:]] = 1
        end
        Input_tensor = B
    else
        for i in 1:B_unique_nnz
            Triangle_check[B.indices[i,:]] = 1
        end
        Input_tensor = A
    end

    for i in 1:size(U,2)
       for j in 1:size(V,2)

            if A_unique_nnz > B_unique_nnz
                matched_tris,gaped_tris,matching = TAME_score(Triangle_check,Input_tensor,V[:,j],U[:,i])
            else
                matched_tris,gaped_tris,matching = TAME_score(Triangle_check,Input_tensor,U[:,i],V[:,j])
            end

            if matched_tris > best_score
                best_matching = matching
                best_score = matched_tris
                best_i = i
                best_j = j
            end
        end
    end
    return best_score, best_i, best_j, best_matching
end

#used when we don't want to recreate the triangle matching dictionary multiple times
function TAME_score(Triangle_Dict::Dict{Array{Int,1},Int},Input_tensor,
                    u::Array{Float64,1},v::Array{Float64,1})

    Match_mapping, _ = rank_one_matching(u,v)
    TAME_score(Triangle_Dict,Input_tensor,Match_mapping)

end

function TAME_score(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,U::Array{Float64,2},V::Array{Float64,2};return_timings=false)

    if return_timings
        (Match_mapping, _), matching_time = @timed low_rank_matching(U,V)
        (triangle_count, gaped_triangles,_), scoring_time = @timed TAME_score(A,B,Match_mapping)
        return triangle_count, gaped_triangles,Match_mapping, matching_time, scoring_time
    else
        Match_mapping,weight = low_rank_matching(U,V)
        TAME_score(A,B,Match_mapping)
    end

end

function TAME_score(A::ThirdOrderSymTensor, B::ThirdOrderSymTensor, X::SparseMatrixCSC{Float64,Int64};return_timings=false)

    if return_timings
        x ,bipartite_matching_time = @timed bipartite_matching(X)
        (triangle_count, gaped_triangles,matching), scoring_time = @timed TAME_score(A,B,Dict(i => j for (i,j) in enumerate(x.match)))
        return triangle_count, gaped_triangles, matching, bipartite_matching_time, scoring_time
    else
        x = bipartite_matching(X)
        return TAME_score(A,B,Dict(i => j for (i,j) in enumerate(x.match)))
    end

end

function TAME_score(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,X::Array{Float64,2};return_timings=false)

    if return_timings
        (_,_,matching,_) ,matching_time = @timed bipartite_matching_primal_dual(X)
        (triangle_count, gaped_triangles,inverted_matching), scoring_time = @timed TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
        return triangle_count, gaped_triangles,inverted_matching, matching_time, scoring_time, matching
    else
        _,_,matching,_ = bipartite_matching_primal_dual(X)
        return TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
    end
end

function TAME_score(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,x::Array{Float64,1};return_timings=false)

    X = reshape(x,A.n,B.n)
    if return_timings
        (_,_,matching,_) ,matching_time = @timed bipartite_matching_primal_dual(X)
        (triangle_count, gaped_triangles,inverted_matching), scoring_time = @timed TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
        return triangle_count, gaped_triangles,inverted_matching, matching_time, scoring_time
    else
        _,_,matching,_ = bipartite_matching_primal_dual(X)
        return TAME_score(A,B,Dict(j => i for (i,j) in enumerate(matching)))
    end

end

"""
   Match_mapping is expected to map V_A -> V_B
"""
function TAME_score(A::ThirdOrderSymTensor,B::ThirdOrderSymTensor,Match_mapping::Dict{Int,Int})

    match_len = length(Match_mapping)

    Triangle_check = Dict{Array{Int,1},Int}()
    gaped_triangles = 0
    triangle_count = 0

    if size(A.indices,1) > size(B.indices,1)

        for i in 1:size(A.indices,1)
            Triangle_check[A.indices[i,:]] = 1
        end

        #invert to map V_B indices to V_A
        Match_mapping = Dict(value => key for (key, value) in Match_mapping)

        for i in 1:size(B.indices,1)
            v_i,v_j,v_k = B.indices[i,:]

            matched_triangle =
              sort([get(Match_mapping,v_i,-1),get(Match_mapping,v_j,-1),get(Match_mapping,v_k,-1)])
    #        println(B.indices[i,:]," -> ",matched_triangle)
            match = get(Triangle_check,matched_triangle,0)
            if match == 1
                triangle_count += 1
            else
                gaped_triangles += 1
            end
        end

    else
        for i in 1:size(B.indices,1)
            Triangle_check[B.indices[i,:]] = 1
        end

        for i in 1:size(A.indices,1)
            v_i,v_j,v_k = A.indices[i,:]
            matched_triangle =
               sort([get(Match_mapping,v_i,-1), get(Match_mapping,v_j,-1),get(Match_mapping,v_k,-1)])

            match = get(Triangle_check,matched_triangle,0)
            if match == 1
                triangle_count += 1
            else
                gaped_triangles += 1
            end
        end
    end

   # triangles, triangle_weight = count_triangles(sub_A,sub_B)

    return triangle_count, gaped_triangles, Match_mapping #sub_A, sub_B

end

function TAME_score(Triangle_Dict::Dict{Array{Int,1},Int},Input_tensor::ThirdOrderSymTensor,
                    Match_mapping::Dict{Int,Int})

    triangle_count = 0
    gaped_triangles = 0

    for i in 1:length(Input_tensor.values)
        v_i,v_j,v_k = Input_tensor.indices[i,:]

        matched_triangle =
          sort([get(Match_mapping,v_i,-1),get(Match_mapping,v_j,-1),get(Match_mapping,v_k,-1)])

        match = get(Triangle_Dict,matched_triangle,0)
        if match == 1
            triangle_count += 1
        else
            gaped_triangles += 1
        end
    end
    return triangle_count, gaped_triangles, Match_mapping
end

"""------------------------------------------------------------------------------
  A function for solving the a maximum bipartite matching problem with a dense
  m x n input X.

  Inputs
  ------
  * X - (Matrix{Float64}):
    The weights of the input bipartite graph.
  * tol - (Float):
    The tolerance to determine whether or not the flow through an edge is at
    capacity.
  * 'normalize_weights' - (Bool):
    Whether or not to normalize the maximum edge weight to be 1.0, implictly
    alters the tolerance to which the problem is solved.

  Outputs
  -------
  * val - (Float):
    The value of the maximum weighted matching.
  * noute - (Int):
    The cardinality of the matching found.
  * match1 - (Array{Int,1}):
    The mapping linking left hand side nodes to the right hand side. The ith
    entry maps node i to a right hand side node.
  * match2 - (Array{Int,1}):
    The mapping linking right hand side nodes to the right hand side. The jth
    entry maps node j to a left hand side node. The last n entries in the array
    represent dummy nodes which are used to account non-matches which may arise
    from negative weights.
------------------------------------------------------------------------------"""
function bipartite_matching_primal_dual(X::Matrix{Float64};tol::Float64=1e-8,
                                       normalize_weights::Bool=false)
    #to get the access pattern right, we must match the right hand side to the left hand side.

	m,n = size(X)
	@assert m >= n  #error occurs when m < n

    # variables used for the primal-dual algorithm
    # normalize ai values # updated on 2-19-2019
    if normalize_weights
        X ./= maximum(abs.(X))
    end

	#initialize variables
    alpha=zeros(Float64,n)
    bt=zeros(Float64,m+n)#beta
    match1 = zeros(Int64,n)
    match2 = zeros(Int64,n+m)
    queue=zeros(Int64,n)
    t=zeros(Int64,m+n)
    tmod = zeros(Int64,m+n)
    ntmod=0

    # initialize the primal and dual variables

	for j = 1:n
		for i=1:m
			if X[i,j] > alpha[j]
			   alpha[j]=X[i,j]
			end
		end
    end

    # dual variables (bt) are initialized to 0 already
    # match1 and match2 are both 0, which indicates no matches

    j=1
    @inbounds while j<=n

        for i=1:ntmod
            t[tmod[i]]=0
        end
        ntmod=0
        # add i to the stack
        head=1
        tail=1
        queue[head]=j
        while head <= tail && match1[j]==0 #queue empty + i is unmatched

 		 	k=queue[head]
			#println("begining of queue loop")

            for i=1:m+1 #iterate over column k

				if i == m+1 #check the dummy node
					i = k+m
				end
				if i == k+m #dummy nodes don't have weight
					if 0.0 < alpha[k] + bt[i] - tol
						continue
					end
                elseif X[i,k] < alpha[k] + bt[i] - tol
					continue
                end # skip if tight

                if t[i]==0
                    tail=tail+1 #put the potential match in the queue
                    if tail <= m
                        queue[tail]=match2[i]
					end
                    t[i]=k  #try vertex k for vertex j
                    ntmod=ntmod+1
                    tmod[ntmod]=i
                    if match2[i]<1  #if i is unmatched
                        while i>0   #unfurl out to get an augmented path
                            match2[i]=t[i]
                            k=t[i]
                            temp=match1[k]
                            match1[k]=i
							i=temp
						end
                        break
                    end
                end
            end
            head=head+1

        end

        #if node j was unable to be matched, update flows and search for new augmenting path
		if match1[j] < 1
            theta=Inf
            for i=1:head-1
                t1=queue[i]
                for t2=1:m
                    if t[t2] == 0 && alpha[t1] + bt[t2] - X[t2,t1] < theta
                        theta = alpha[t1] + bt[t2] - X[t2,t1]
                    end
                end
				#check t1's dummy node
				if t[t1 + m] == 0 && alpha[t1] + bt[t1 + m] < theta
                        theta = alpha[t1] + bt[t1 + m]
				end
            end

            for i=1:head-1
                alpha[queue[i]] -= theta
            end
            for i=1:ntmod
                bt[tmod[i]] += theta
            end
            continue
        end

        j=j+1
    end

	#count
    val=0.0
    for j=1:n
        for i=1:m
            if i==match1[j]
                val=val+X[i,j]
            end
        end
    end

	#count how many are properly matched
    noute = 0
    for j=1:n
        if match1[j]<=m
            noute=noute+1
        end
	end

    return val,noute,match1, match2
end


text:
#=
BFSParentsMultiSource:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-10
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra

function insert(P, R, s, n)
    for i = 1:s
        for j = 1:n
            if(P[i, j]!=nothing && R[i, j]==nothing) #*
                R[i, j] = P[i, j]
            end
        end
    end

    return R

end

#A is the input matrix, S are the sources, n is the number of nodes, s is the number of sources, Output: Parents matrix
function ms_bfsp(A, S, n, s)

    index = GBMatrix{Int64}(s, n)
            for i = 1:n
                for j = 1:s
                    index[j, i] = i
                end
            end

    #frontier matrix -> F=S source nodes insrted in frontier
    F = GBMatrix{Int64}(n, n)
    F = S

    #parent matrix -> P[S]=S because the source is already visited
    P = GBMatrix{Int64}(n, n)
    P = S

    #result matrix -> R[S]=0 because the parents of the source are considered node 0
    R = GBMatrix{Int64}(s, n)
    LR = GBMatrix{Int64}(s, n)
    R = copy(S)
    for i = 1:n
        for j = 1:s
            if(R[j, i]!=nothing)
                R[j, i] = 0
            end
        end
    end

    print("STARTING ...\n\n")
    for i = 1:n-1
            F = mul(F, A, Semirings.MIN_FIRST, mask=P, desc=Descriptors.RC)
            P = F[:, :, mask=F, desc=Descriptors.S]
            F = index[:, :, mask=F, desc=Descriptors.S]
            R = insert(P, R, s, n)
            if(R == LR)
                break
            end
            print("\nFRONTIER UPDATED\n\n")
            print(F)
            LR = copy(R)
            print("\n\nITERATION FINISHED\n\n")

    end

    print("COMPLETED\n\n")
    print("RESULT\n\n")
    print(R)

    return R

end


text:
export VoiceRegion

"""
A region for a [`Guild`](@ref)'s voice server.
More details [here](https://discordapp.com/developers/docs/resources/voice#voice-region-object).
"""
struct VoiceRegion
    id::String
    name::String
    vip::Bool
    optimal::Bool
    deprecated::Bool
    custom::Bool
end
@boilerplate VoiceRegion :dict :docs :lower :merge


text:
using JuAFEM

# Create a quadrature rule
quad_rule = QuadratureRule{2, RefCube}(:legendre, 2)

# Create a function space
func_space = Lagrange{2, RefCube, 1}()

# Use these to create an FEValues object
cell_values = CellScalarValues(quad_rule, func_space)
fieldnames(cell_values) |> display
println()

# Coordinates of 1 cell in grid
x = Vec{2, Float64}[
  Vec{2}((0.0, 0.0)),
  Vec{2}((1.5, 0.0)),
  Vec{2}((2.0, 2.0)),
  Vec{2}((0.0, 1.0))
];

# Use reinit! to update cell_values on an element
reinit!(cell_values, x)

# Query cell_values shape function in quadrature point 3
shape_value(cell_values, 3, 1) |> display
println()

# DErivative of shape function 2 in quadrature point 1
shape_gradient(cell_values, 1, 2) |> display
println()

# We can also evaluate on the finite element basis
T = [0.0, 1.0, 2.0, 1.5]; # nodal values
function_value(cell_values, 3, T) |> display # value of T in 3rd quad point
println()
function_gradient(cell_values, 1, T) |> display # value of grad(T) in 1st quad point
println()

# Or for vector valued functions
u = Vec{2, Float64}[
  Vec{2}((0.0, 0.0)),
  Vec{2}((3.5, 2.0)),
  Vec{2}((2.0, 2.0)),
  Vec{2}((2.0, 1.0))
]; # nodal vectors

function_value(cell_values, 2, u) |> display # value of u in 2nd quad point
println()
function_symmetric_gradient(cell_values, 3, u)  |> display# sym(grad(u)) in 3rd quad point


text:
# # Example: Happy Eyeballs

using Sockets
using Julio
using Test

# An implementation of Happy Eyeballs (RFC 8305)
# <https://datatracker.ietf.org/doc/html/rfc8305>
#
# (TODO: verify)

function happy_eyeballs(host, port; delay = 0.3)
    addrs = getalladdrinfo(host)
    winner = Julio.Promise()
    decided = nothing
    allsockets = TCPSocket[]
    Julio.withtaskgroup() do tg
        issuccess = false
        closing = Threads.Atomic{Bool}(false)
        try
            failureevents = []
            for ip in addrs
                failed = Julio.Promise{Nothing}()
                push!(failureevents, failed)
                socket = TCPSocket()
                push!(allsockets, socket)
                Julio.spawn!(tg) do
                    inject_latencies(ip == addrs[end])                      #src
                    try
                        connect(socket, ip, port)
                    catch err
                        if err isa Base.IOError
                            failed[] = nothing
                            return
                        elseif closing[]
                            return
                        end
                        rethrow()
                    end
                    Julio.tryput!(winner, socket)
                end
                Julio.select(
                    (fetch, winner) => socket -> begin
                        decided = socket
                        true
                    end,
                    (fetch, failed) => _ -> false,
                    (sleep, delay) => _ -> false,
                ) && break
            end
            if decided === nothing
                for failed in failureevents
                    Julio.select(
                        (fetch, winner) => socket -> begin
                            decided = socket
                            true
                        end,
                        (fetch, failed) => _ -> false,
                    ) && break
                end
            end
            issuccess = true
        finally
            closing[] = true
            for socket in allsockets
                if !issuccess || socket !== decided
                    close(socket)
                end
            end
        end
    end
    for socket in allsockets                                                #src
        if socket !== decided                                               #src
            @test !isopen(socket)                                           #src
        end                                                                 #src
    end                                                                     #src
    return decided
end

# * [Nathaniel J. Smith - Trio: Async concurrency for mere mortals - PyCon 2018 - YouTube](https://www.youtube.com/watch?v=oLkfnc_UMcE)
# * [Two Approaches to Structured Concurrency - 250bpm](https://250bpm.com/blog:139/)

function test_happy_eyeballs()
    socket = happy_eyeballs("httpbin.org", 80)
    try
        @test socket isa TCPSocket
        @test isopen(socket)
    finally
        close(socket)
    end
end

# For reproducing "TCPSocket is not in initialization state" error:         #src
@isdefined(inject_latencies) || begin                                       #src
    inject_latencies(_) = nothing                                           #src
end                                                                         #src


text:
"""
    abstract type CombinatorialInstance

An instance of a combinatorial problem. An instance of this type only contains
the data required to solve a problem (e.g., the weights, the costs), and no
solver-specific data.
"""
abstract type CombinatorialInstance end

"""
    abstract type CombinatorialObjective

A specific objective function for a combinatorial problem. In general, the
objective is to either minimise or maximise the natural objective function,
but more specific objectives can be provided with other subtypes.

Objects of type `CombinatorialObjective` are supposed to be given as arguments
when building a `CombinatorialInstance`. In any case,
`objective(::CombinatorialInstance)` can be used to retrieve the objective of
an instance.
"""
abstract type CombinatorialObjective end

"""
    abstract type CombinatorialSolution

The solution returned by a solver when calling the `solve` function.
This type should be specified for each combinatorial instance to have
a unified solution type.

Two fields must be implemented: `instance` must be the instance that has been
solved, and `variables` must indicate the actual solution. `variables` can only
contain the discrete items that are chosen as part of the solution (i.e., the
elements of this object can be used as indices for the rewards of a
`CombinatorialInstance`). Solution objects can define more properties
if needed.
"""
abstract type CombinatorialSolution end

"""
    abstract type CombinatorialAlgorithm

One specific combinatorial algorithm to solve an instance. Given as a parameter
to `solve`, it specifies the way the instance should be solved.

As often as possible, generic objects are provided with usual names (e.g.,
a greedy algorithm can be proposed for many different problems: the actual
implementation will depend on the type of the combinatorial instance at hand).
"""
abstract type CombinatorialAlgorithm end

"""
    abstract type CombinatorialLinearFormulation

The base type for a linear formulation (LP, IP) for a given combinatorial
problem. The same problem might have several different formulations.
"""
abstract type CombinatorialLinearFormulation <: CombinatorialAlgorithm end

"""
    abstract type CombinatorialVariation

Variations of the base combinatorial set, e.g. to include other constraints.
A variation is an object that is built on top of a combinatorial instance
that contains all the required information for the supplementary constraints
"""
abstract type CombinatorialVariation <: CombinatorialInstance end

"""
    solve(i::CombinatorialInstance, algo::CombinatorialAlgorithm)

Solve the given combinatorial instance `i` using the provided algorithm
`algo` (potentially containing parameters for the solving process). The
returned value is a `CombinatorialSolution`.

If the requested algorithm is not available for this type of instance,
a `MethodError` is thrown.
"""
function solve end

"""
    objective(i::CombinatorialInstance)

Returns the objective associated with this instance.

By default, instances are supposed to provide an `objective` member, but this
function can be overridden otherwise (for instance, when only one objective
is supported).
"""
function objective(i::CombinatorialInstance)
    return i.objective
end

function objective(i::CombinatorialVariation)
    return objective(i.instance)
end

"""
    dimension(i::CombinatorialInstance)

Returns the dimension of this instance.

All instances should implement this method.
"""
function dimension end # i::CombinatorialInstance

function dimension(i::CombinatorialVariation)
    return dimension(i.instance)
end

"""
    formulation(i::CombinatorialInstance, f::CombinatorialLinearFormulation)

Returns the LP formulation `f` for this combinatorial instance. It returns
two arguments:

- the JuMP model
- the decision variables

For `CombinatorialVariation`, it also returns a third argument: the
supplementary constraint encoded by the `CombinatorialVariation` on top fo the
`CombinatorialInstance`.
"""
function formulation end # i::CombinatorialInstance, f::CombinatorialLinearFormulation; solver=nothing

"""
    reward(i::CombinatorialInstance, variable)

Returns the reward for the given `variable`.
"""
function reward(i::CombinatorialInstance, var)
    return i.rewards[var]
end

"""
    value(s::CombinatorialSolution)

Returns the value of a solution, i.e. its total reward. For a maximisation
problem, an infeasible solution has a `value` of `-Inf`.

    value(s::MinBudgetedSolution, budget::Int)

Returns the value of the solution corresponding to the given budget. This
function is likely to throw an error for `SingleMinBudgetedSolution`, but is
ensured not to for `MultipleMinBudgetedSolution`.
"""
function value(s::CombinatorialSolution)
    if -1 in s.variables || length(s.variables) == 0
        return -Inf
    end

    return sum(reward(s.instance, i) for i in s.variables)
end

"""
    make_solution(::CombinatorialInstance, ::Dict{K, Float64}) where {K}

Creates a solution object of the right type for the given combinatorial
instance. The dictionary gives the solution that should be included in the
returned object.

As only binary decision variables are supported, the caller of this function
may provide a dictionary such that *either* only the variables that are
included in the solution (i.e. `variable => 1.0` entries) *or* all the
variables, with associated values of 0 (not in the solution) or 1
(in the solution). For instance, if the solution contains only the edge
`1`-`2` out of the two possible edges (the other one being `1`-`3`), there are
two equivalent dictionaries that can be passed: either
`Dict(Edge(1, 2) => 1.0)` or `Dict(Edge(1, 2) => 1.0, Edge(1, 3) => 0.0)`.
"""
function make_solution end # i::CombinatorialInstance, vals::Dict{K, Float64}