text:
# Example from Soare et al.

using JLD2;
using Distributed;
using Printf;
using IterTools;
using Distributions

@everywhere include("runit_bis.jl");
@everywhere include("../thresholds.jl");
include("../experiment_helpers.jl");

# setup

dist = Gaussian();
dim = 2
μ = zeros(dim);
µ[1] = 2.0;
# μ = [.9, .7, .5];

arms = Vector{Float64}[]
for k = 1:dim
    v = zeros(dim)
    v[k] = 1.0
    push!(arms, v)
end
ω = 0.1
v = zeros(dim);
v[1] = cos(ω);
v[2] = sin(ω);
push!(arms, v)

pep = LinearBestArm(dist, arms);

K = length(arms)
spanning_weights = ones(K)
spanning_weights[K] = 0.0
spanning_weights /= sum(spanning_weights)

# srs = [
#     ConvexGame(CTracking),
#     LearnerK(CTracking),
#     RoundRobin(),
#     FixedWeights(spanning_weights),
#     Static("G"),
#     Static("XY"),
#     XYAdaptive(),
#     LinGapE(),
# ]
srs = [SLT3C(), SLGapE(), LT3C("G"), LinGapE()]

# δs = (0.1, 0.01, 0.0001);
δs = (0.1,);
βs = GK16.(δs);

repeats = 100;
seed = 123;

# compute

@time data = map(  # TODO: replace by pmap (it is easier to debug with map)
    ((sr, i),) -> runit(seed + i, sr, μ, pep, βs, δs),
    Iterators.product(srs, 1:repeats),
);

dump_stats(pep, μ, δs, βs, srs, data, repeats);

# save

@save isempty(ARGS) ? "experiment_bai1.dat" : ARGS[1] dist μ pep srs data δs βs repeats seed

# visualise by loading viz_bai1.jl


text:
using HTTP
using Gumbo
using JSON

url = "https://docs.github.com/en/developers/webhooks-and-events/webhooks/webhook-events-and-payloads"
cache_dir = "cache"
cache_file = "$cache_dir/cache.html"

function read_cache()
    files = readdir(cache_dir)
    if isempty(files)
        return []
    end
    json = filter((f) -> endswith(f, ".json"), files)
    return Dict(map((f) -> (f[1:end-5], open(v->read(v, String), "$cache_dir/$f")), json))
end

function visitor(source::HTMLElement, predicate::Function, ans::Vector{HTMLElement})
    for e in source.children
        if !isa(e, HTMLElement)
            continue
        end
        if predicate(e)
            push!(ans, e)
        end
        visitor(e, predicate, ans)
    end
end

function extractor(source::HTMLElement)
    result = ""
    for e in source.children
        if !isa(e, HTMLText)
            result *= extractor(e)
            continue
        end
        result *= e.text
    end
    return result
end

function fetch()
    r = HTTP.request("GET", url)
    if r.status != 200
        return
    end

    body = String(r.body)
    open(cache_file, "w") do io
        write(io, body)
    end
    return body
end

function parse_html(body::String)
    html = parsehtml(body).root
    b = first(filter(e -> tag(e) == :body, html.children))
    raw_titles = HTMLElement[]
    raw_blocks = HTMLElement[]
    visitor(b,
            (e) -> tag(e) == :h2 &&
                getattr(e, "class", "") == "" &&
                getattr(e, "id", "") != "webhook-payload-object-common-properties",
            raw_titles)
    visitor(b,
            (e) -> tag(e) == :code &&
                getattr(e, "class", "") == "hljs language-json",
            raw_blocks)
    titles = map(extractor, raw_titles)
    blocks = map(extractor, raw_blocks)
    return Dict(zip(titles, blocks))
end

function main()
    if isempty(ARGS)
        error("key is required.")
    end

    payloads = read_cache()
    if isempty(payloads)
        if isfile(cache_file)
            body = open(f->read(f, String), cache_file)
        else
            body = fetch()
        end

        payloads = parse_html(body)
        for (t, b) in payloads
            open("$cache_dir/$t.json", "w") do io
                write(io, b)
            end
        end
    end

    key = first(ARGS)
    if !haskey(payloads, key)
        error("$key cannot be found.")
    end

    JSON.print(JSON.parse(payloads[key]), 4)
end

main()

ext: jl

max_stars_repo_name: CORDEA/gh-webhook-payload-gen

max_issues_repo_name: CORDEA/gh-webhook-payload-gen

max_forks_repo_name: CORDEA/gh-webhook-payload-gen


text:
sparse_options = (; conic_solver = sparse_hypatia)

insts = OrderedDict()

insts["test"] = [
    ((1, 2, 15, 50.0, true),),
    ((1, 2, 15, 50.0, false),),
    ((3, 1, 40, 0.0, true), sparse_options),
    ((3, 1, 40, 0.0, false), sparse_options),
    ((2, 2, 40, 100.0, true), sparse_options),
]

function polyregression_insts(use_nat::Bool, m_max::Int)
    return [((2, 2, m, 25.0, use_nat), sparse_options) for m in vcat(30, 30:5:m_max)]
end

insts["nat"] = polyregression_insts(true, 80)
insts["ext"] = polyregression_insts(false, 70)

return (PolyRegression, insts)


text:
module ComponentArrays

using ArrayInterface
# using ChainRulesCore
using LinearAlgebra
using Requires

const FlatIdx = Union{Int, CartesianIndex, AbstractArray}
const FlatOrColonIdx = Union{FlatIdx, Colon}

include("utils.jl")
include("lazyarray.jl")
include("axis.jl")
include("componentindex.jl")
include("componentarray.jl")
include("set_get.jl")
include("similar_convert_copy.jl")
include("broadcasting.jl")
include("math.jl")
include("show.jl")
include("plot_utils.jl")
# include(joinpath("if_required", "chainrulescore.jl"))

# If using for differential equations, the Array(sol) overload in RecursiveArrayTools will
# concatenate the ComponentVectors while preserving their ComponentArrayness
required(filename) = include(joinpath("if_required", filename))

function __init__()
    @require ConstructionBase="187b0558-2788-49d3-abe0-74a17ed4e7c9" required("constructionbase.jl")
    @require DiffEqBase="2b5f629d-d688-5b77-993f-72d75c75574e" begin
        @require RecursiveFactorization="f2c3362d-daeb-58d1-803e-2bc74f2840b4" required("diffeqbase.jl")
    end
    @require SciMLBase="0bca4576-84f4-4d90-8ffe-ffa030f20462" required("scimlbase.jl")
    @require RecursiveArrayTools="731186ca-8d62-57ce-b412-fbd966d074cd" required("recursivearraytools.jl")
    @require StaticArrays="90137ffa-7385-5640-81b9-e52037218182" required("staticarrays.jl")
    @require Zygote="e88e6eb3-aa80-5325-afca-941959d7151f" required("zygote.jl")
    # @require ChainRulesCore="d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4" required("chainrulescore.jl")
end

export AbstractAxis, Axis, PartitionedAxis, ShapedAxis, ViewAxis, FlatAxis

export ComponentArray, ComponentVector, ComponentMatrix

export getdata, getaxes, fastindices, labels, label2index, valkeys

end


text:
# test ro_view
#using ODLCommonTools
#using FactCheck

function sum_cols(A::AbstractMatrix{T}) where T

  val = zero(T)
  for i=1:size(A, 2)
    tmp = ro_sview(A, :, i)
    val += sum(tmp)
  end

  return val
end

function sumit(A::AbstractArray{T}) where T

  val = zero(T)
  for i=1:length(A)
    val += A[i]
  end

  return val
end

@testset "----- Testing ROView -----" begin

  a = rand(100, 1000)
  a2 = rand(100, 1000, 1)

  b = ROView(a)
  @test ( size(b) )== size(a)
  @test ( length(b) )== length(a)

  @test isapprox( norm(a - b), 0.0)

  if Base.JLOptions().can_inline == true
    # check that allocating a ROView does not allocate
    if !ODLCommonTools.safe_views
      @time sum_cols(a)
      bytes_alloc = @allocated sum_cols(a)
      @test  bytes_alloc  < 100

      # test ROView of sview
      atmp = sview(a2, :, :, 1)
      @time sum_cols(atmp)
      bytes_alloc = @allocated sum_cols(atmp)
      @test  bytes_alloc  < 100
    end

    # check that indexing a ROView does not allocate
    sumit(b)
    bytes_alloc = @allocated sumit(b)
    @test  bytes_alloc  < 100
  end

  # test reinterpert
  A = rand(9, 9)
  Av = unsafe_aview(A, 6:9, 1)

  Av2 = reinterpret(Complex128, Av, (2, 1))

  @test ( Av2[1] )== Complex128(Av[1], Av[2])
  @test ( Av2[2] )== Complex128(Av[3], Av[4])

  # test that composition works

  # ro_sview of an sview
  arr = rand(2, 3, 4)
  v1 = sview(arr, :, :, 1)
  v2 = sview(v1, :, 1)
  v3 = ro_sview(v1, :, 1)
  for i=1:length(v2)
    @test v2[i] == v3[i]
  end

  # test sview of ro_sview
  v1 = ro_sview(arr, :, :, 1)
  v2 = sview(v1, :, 1)
  v3 = arr[:, 1, 1]

  for i=1:length(v2)
    @test v2[i] == v3[i]
  end

end


text:
module BSONSerializer

using Mongoc: BSON, BSONObjectId

using Dates, InteractiveUtils

export @BSONSerializable

include("types.jl")
include("encoding.jl")
include("codegen.jl")
include("serialize.jl")

end # module


text:
module Examples_sep
using LinearAlgebra
using Random

srcDir = dirname(@__FILE__)*"\\"
include(srcDir *"Utils_states.jl")
using .Utils_states ; const us = Utils_states

export get_sep_examples,
       gen_ρ_RAND

## Separable States
"""
Separable state examples:
"""
function get_sep_example()
    ρ = Dict()

##  T1CD12 i-ii
    ρ["S1a"] = get_ρ_T1CD12(1)
    ρ["S1b"] = get_ρ_T1CD12(2)
    # ρ["T1CD12iii"]  = get_ρ_T1CD12(3)

## T2CD12 i,ii, v
    ρ["S2a"]   = get_ρ_T2CD12(1)
    ρ["S2b"]   = get_ρ_T2CD12(2)
    # ρ["T2CD12iii"]  = get_ρ_T2CD12(3)
    # ρ["T2CD12iv"]  = get_ρ_T2CD12(4)
    # ρ["T2CD12v"]    = get_ρ_T2CD12(5)
## Ex26CD12
    ρ["S3"]   = get_ρ_Ex26CD12()

## Ex3-8DNY20
    ρ["S4"] = get_ρ_DNY20()

## DTT00a
    ρ["S5a"] = get_ρ_DTT00a(3)
    ρ["S5b"] = get_ρ_DTT00b()
## Random batch
## ℝeal
    ρ["RRai"]   = gen_ρ_RAND((3,3), 5)
    ρ["RRaii"]  = gen_ρ_RAND((3,3), 10)
    ρ["RRaiii"] = gen_ρ_RAND((3,3), 15)
    ρ["RRaiv"]  = gen_ρ_RAND((3,3), 20)
    ρ["RRav"]   = gen_ρ_RAND((3,3), 25)

    ρ["RRbi"]   = gen_ρ_RAND((2,4), 5)
    ρ["RRbii"]  = gen_ρ_RAND((2,4), 10)
    ρ["RRbiii"] = gen_ρ_RAND((2,4), 15)
    ρ["RRbiv"]  = gen_ρ_RAND((2,4), 20)
    ρ["RRbv"]   = gen_ρ_RAND((2,4), 25)

    ρ["RRci"]   = gen_ρ_RAND((2,5), 5)
    ρ["RRcii"]  = gen_ρ_RAND((2,5), 10)
    ρ["RRciii"] = gen_ρ_RAND((2,5), 15)
    ρ["RRciv"]  = gen_ρ_RAND((2,5), 20)
    ρ["RRcv"]   = gen_ρ_RAND((2,5), 25)

    ρ["RRdv"]   = gen_ρ_RAND((4,4), 25)

## ℂomplex
    ρ["RCai"]   = gen_ρ_RAND((3,3), 5, true)
    ρ["RCaii"]  = gen_ρ_RAND((3,3), 10, true)
    ρ["RCaiii"] = gen_ρ_RAND((3,3), 15, true)
    ρ["RCaiv"]  = gen_ρ_RAND((3,3), 20, true)
    ρ["RCav"]   = gen_ρ_RAND((3,3), 25, true)

    ρ["RCbi"]   = gen_ρ_RAND((2,4), 5, true)
    ρ["RCbii"]  = gen_ρ_RAND((2,4), 10, true)
    ρ["RCbiii"] = gen_ρ_RAND((2,4), 15, true)
    ρ["RCbiv"]  = gen_ρ_RAND((2,4), 20, true)
    ρ["RCbv"]   = gen_ρ_RAND((2,4), 25, true)

    ρ["RCci"]   = gen_ρ_RAND((2,5), 5, true)
    ρ["RCcii"]  = gen_ρ_RAND((2,5), 10, true)
    ρ["RCciii"] = gen_ρ_RAND((2,5), 15, true)
    ρ["RCciv"]  = gen_ρ_RAND((2,5), 20, true)
    ρ["RCcv"]   = gen_ρ_RAND((2,5), 25, true)

    ρ["RCd"]   = gen_ρ_RAND((3,4), 25, true)
    ρ["RCe"]   = gen_ρ_RAND((2,6), 25, true)
 ## Real
 #    ρ["RANDa1"]  = gen_ρ_RAND(3, 5)
 #    ρ["RANDa2"]  = gen_ρ_RAND(3, 10)
 #    ρ["RANDa3"]  = gen_ρ_RAND(3, 15)
 #    ρ["RANDa4"]  = gen_ρ_RAND(3, 20)
 #    ρ["RANDa5"]  = gen_ρ_RAND(3, 25)
 #
 #    ρ["RANDb1"]  = gen_ρ_RAND(4, 5)
 #    ρ["RANDb2"]  = gen_ρ_RAND(4, 10)
 #    ρ["RANDb3"]  = gen_ρ_RAND(4, 15)
 #    ρ["RANDb4"]  = gen_ρ_RAND(4, 20)
 #    ρ["RANDb5"]  = gen_ρ_RAND(4, 25)
 # ## Complex
 #    ρ["RANDℂa1"]  = gen_ρ_RAND(3,  5,true)
 #    ρ["RANDℂa2"]  = gen_ρ_RAND(3, 10,true)
 #    ρ["RANDℂa3"]  = gen_ρ_RAND(3, 15,true)
 #    ρ["RANDℂa4"]  = gen_ρ_RAND(3, 20,true)
 #    ρ["RANDℂa5"]  = gen_ρ_RAND(3, 25,true)
 #
 #    ρ["RANDℂb1"]  = gen_ρ_RAND(4,  5,true)
 #    ρ["RANDℂb2"]  = gen_ρ_RAND(4, 10,true)
 #    ρ["RANDℂb3"]  = gen_ρ_RAND(4, 15,true)
 #    ρ["RANDℂb4"]  = gen_ρ_RAND(4, 20,true)
 #    ρ["RANDℂb5"]  = gen_ρ_RAND(4, 25,true)

    for key in keys(ρ)
        ρ[key] = Utils_states.maketraceone(ρ[key])
    end
    return ρ
end

## Examples

""" http://arxiv.org/abs/1210.0111v2"""
function get_ρ_T1CD12(i)
    ρ = Dict()
    # |00><00| + |11><11|
    ρ[1] =  psep(1,1,2) + psep(2,2,2)
    # |00><00| + |11><11| + (|0> + |1>)⊗(|0> + |1>)(<0| + <1|)⊗(<0| + <1|)
    ρ[2] = psep(1,1,2) + psep(2,2,2) + sq(ψ(2))
    return ρ[i]
end

""" http://arxiv.org/abs/1210.0111v2"""
function get_ρ_T2CD12(i)
    ρ = Dict()
    # |00><00| + |11><11| + |12><12|
    ρ[1] = psep(1,1,3) + psep(2,2,3) + psep(2,3,3)
    # |00><00| + |01><01| + |11><11| + |12><12|
    ρ[2] = psep(1,1,3) + psep(1,2,3) + psep(2,2,3) + psep(2,3,3)
    # |00><00| + |01><01| + |02><02| + |11><11| + |12><12|
    ρ[5] = psep(1,1,3) + psep(1,2,3) +  psep(1,3,3) + psep(2,2,3) + psep(2,3,3)
    return ρ[i]
end

""" http://arxiv.org/abs/1210.0111v2"""
function get_ρ_Ex26CD12()
    ρ = [4 0 0 0 0 0
         0 4 2 0 0 2
         0 2 2 1 -1 0
         0 0 1 2 1 -1
         0 0 -1 1 5 1
         0 2 0 -1 1 2]
    return ρ
end

"""https://arxiv.org/abs/2011.08132v1"""
function get_ρ_DNY20()
    H_flat = [i₁*j₁ + i₂*j₂ for i₁ in 1:3 for i₂ in 1:3 for j₁ in 1:3 for j₂ in 1:3 ]
    return  reshape(H_flat,9,9)
end

"""https://arxiv.org/abs/quant-ph/9904005v2 page 7 eq. (11)"""
function get_ρ_DTT00a(n)
    #@assert 0 ≤ f ≤ 1/n  # to be separable
    f=1/(n)
    Ψn = (1/sqrt(n))*sum([kron(us.eᵢ(n,i),us.eᵢ(n,i)) for i in 1:n])
    ρᵥᵥf = f* us.sq(Ψn) + (1-f)/(n^2)*I(n^2)
    ρf =  us.take_Pᵀ(ρᵥᵥf,1,(n,n))
end

"""https://arxiv.org/abs/quant-ph/9904005v2 page 9 eq. (13)"""
function get_ρ_DTT00b()
    N = 2/(sqrt(5 + sqrt(5)))
    h = 0.5*(sqrt(1 + sqrt(5)))
    Ñ = sqrt(2/sqrt(5))

    v(i) = N .* [cos(2*π*i/5),sin(2*π*i/5),h]
    w(j) = Ñ .* sqrt(cos(π/5))*[cos(2*π*j/5),sin(2*π*j/5),cos(4*π*j/5),sin(4*π*j/5)]

    (1/140)*(I(12) - sum([kron(us.sq(v(i)),us.sq(w(i))) for i in 1:4 ]))
end

"""generates a matrix of the form ∑ʳaᵀa⊗bᵀb, with a,b ∈ uniform random entry-wise.  d ∈ {2,3,4} , r ∈ [9]"""
function gen_ρ_RAND(d, r::Integer, isℂ = false)
    Random.seed!(343) ; d₁,d₂ = d ; ρ = zeros(d₁*d₂,d₁*d₂)
    if isℂ
        a = [randn(1,d₁) + im*randn(1,d₁) for i in 1:r ]
        b = [randn(1,d₂) + im*randn(1,d₂) for i in 1:r ]
    else
        a = [randn(1,d₁) for i in 1:r ]
        b = [randn(1,d₂) for i in 1:r ]
    end
    for i in 1:r
        A = transpose(a[i])*conj(a[i])
        B = transpose(b[i])*conj(b[i])
        ρ = kron(A,B) + ρ
    end
    return ρ
end

# """ http://arxiv.org/abs/quant-ph/9605038v2 """
# function get_ρ_HHH1(p  = 0.5)
#     # Random.seed!(343)
#     a  = 0.5; b = 0.8; # arbitary possitive numbers
#     ψ₁ = a*Utils_states.ψ(2,2,2) + b*Utils_states.ψ(1,1,2)
#     ψ₂ = a*Utils_states.ψ(2,1,2) + b*Utils_states.ψ(1,2,2)
#     ρ_temp = p*Utils_states.sq(ψ₁) + (1 - p)*Utils_states.sq(ψ₂)
#     return ρ_temp
# end

end

max_stars_repo_name: JAndriesJ/sep-rank

max_issues_repo_name: JAndriesJ/sep-rank

max_forks_repo_name: JAndriesJ/sep-rank


text:
using PlanktonIndividuals, Serialization

grid = RectilinearGrid(size = (1, 1, 1), x = (0,32), y = (0,32), z = (0,-32))

model = PlanktonModel(CPU(), grid;
                      N_species = 5,
                      N_individual = 1024,
                      max_individuals = 1024*10)

function tot_mass(nut, g)
    mass = zeros(g.Nx, g.Ny, g.Nz)
    for i in 1:g.Nx
        for j in 1:g.Ny
            for k in 1:g.Nz
                mass[i,j,k] = nut[i+g.Hx, j+g.Hy, k+g.Hz] * PlanktonIndividuals.Grids.volume(i+g.Hx, j+g.Hy, k+g.Hz, g)
            end
        end
    end
    return sum(mass)
end

TP = tot_mass(model.nutrients.PO4.data, grid) +
     tot_mass(model.nutrients.DOP.data, grid) +
     tot_mass(model.nutrients.POP.data, grid)
TP = TP + sum(model.individuals.phytos.sp1.data.Pq .+
              model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC) +
          sum(model.individuals.phytos.sp2.data.Pq .+
              model.individuals.phytos.sp2.data.Bm .* model.individuals.phytos.sp2.p.R_PC) +
          sum(model.individuals.phytos.sp3.data.Pq .+
              model.individuals.phytos.sp3.data.Bm .* model.individuals.phytos.sp3.p.R_PC) +
          sum(model.individuals.phytos.sp4.data.Pq .+
              model.individuals.phytos.sp4.data.Bm .* model.individuals.phytos.sp4.p.R_PC) +
          sum(model.individuals.phytos.sp5.data.Pq .+
              model.individuals.phytos.sp5.data.Bm .* model.individuals.phytos.sp5.p.R_PC)

sim = PlanktonSimulation(model, ΔT = 60, iterations = 10)

update!(sim)

TPt = tot_mass(model.nutrients.PO4.data, grid) +
      tot_mass(model.nutrients.DOP.data, grid) +
      tot_mass(model.nutrients.POP.data, grid)
TPt=TPt + sum(model.individuals.phytos.sp1.data.Pq .+
              model.individuals.phytos.sp1.data.Bm .* model.individuals.phytos.sp1.p.R_PC) +
          sum(model.individuals.phytos.sp2.data.Pq .+
              model.individuals.phytos.sp2.data.Bm .* model.individuals.phytos.sp2.p.R_PC) +
          sum(model.individuals.phytos.sp3.data.Pq .+
              model.individuals.phytos.sp3.data.Bm .* model.individuals.phytos.sp3.p.R_PC) +
          sum(model.individuals.phytos.sp4.data.Pq .+
              model.individuals.phytos.sp4.data.Bm .* model.individuals.phytos.sp4.p.R_PC) +
          sum(model.individuals.phytos.sp5.data.Pq .+
              model.individuals.phytos.sp5.data.Bm .* model.individuals.phytos.sp5.p.R_PC)

@testset "PlanktonIndividuals 0D tests:" begin
    @test isapprox(TP,TPt; atol=1e1)
end


text:
using Pkg
Pkg.activate(ENV["JULIA_DQMC"])

include(joinpath(ENV["JULIA_DQMC"], "src/dqmc_framework.jl"))

input = "input.in.xml" # choose DQMC input file

using Suppressor, DataFrames, CSV

function measure_mc_memory_usage()
    df = DataFrame(L=Int[], B=Int[], mem=Float64[])

    @suppress for L in [4,8,10,12,14]
        for beta in [5,10,20,40]
            p = Params()
            p.output_file = "live.out.h5.running"
            xml2parameters!(p, input)
            p.lattice_file = ""
            p.L = L
            p.beta = beta
            p.slices = beta * 10
            deduce_remaining_parameters(p)
            mc = DQMC(p)
            init!(mc)

            push!(df, [L, beta, memory_usage(mc)])
        end
    end

    return df
end

df = measure_mc_memory_usage()

CSV.write("mc_memory.csv", df)

using PyPlot, Polynomials
plt.style.use("publication_tex")
# plt.style.use("default")

function plot_memusage(df)
    g = groupby(df, :L);

    for (i, gl) in enumerate(g)
        x = gl[:B]
        y = gl[:mem]
        p = polyfit(x, y, 1)
        m = round(p.a[2], digits=1)
        b = round(p.a[1], digits=1)
        plot(x, y, "o-", color="C$i", label="L = $(gl[:L][1]) ($(m)x + $b)")
        # plot(1:40, p.(1:40), color="k")
    end
    handles, labels = gca().get_legend_handles_labels()
    legend(reverse(handles), reverse(labels), frameon=false, prop=Dict("size" => 15))
    xlim([3,42])
    ylim([-100, 1800])
    xlabel(L"inverse temperature $\beta$")
    ylabel("static memory usage in MB")
    tight_layout()
    savefig("mc_memory.pdf")
    nothing
end

plot_memusage(df)

max_stars_repo_name: crstnbr/dqmc

max_issues_repo_name: crstnbr/dqmc

max_forks_repo_name: crstnbr/dqmc


text:
### A Pluto.jl notebook ###
# v0.14.5

using Markdown
using InteractiveUtils

# ╔═╡ 2c554897-e486-418f-9cd4-a3bb3f00fc31
using Pkg, DrWatson

# ╔═╡ 19cc5ec6-f7cb-4736-98c3-498d91eb2414
begin
	@quickactivate "NumericalMethodsforEngineers"
	using BarycentricInterpolation
end

# ╔═╡ 53285b2c-3704-4adc-9bdd-9dea619b5140
md" ## BarycentricInterpolation.jl"

# ╔═╡ 74272021-5f45-4161-bea4-8386697de280
begin
	x = [1.0, 3.0, 6.0, 5.0]
	y = [1.0, 5.0, 10.0, 9.0]
	xi = [2.0, 4.5]
end

# ╔═╡ 587d934e-161e-403c-831b-a61f313f5ecc
indxs = sortperm(x)

# ╔═╡ 7ba7dd56-4ee9-4c5c-97a8-5f99a9ba507e
p = ArbitraryPolynomial(x[indxs]);

# ╔═╡ 880d1ae8-5a5a-463a-8e35-a73bdc40208a
interpolate(p, y[indxs], xi)

# ╔═╡ c2fcf437-bd86-4b24-8c77-29702cf5373d
md" ## End of BarycentricInterpolation.jl"

# ╔═╡ Cell order:
# ╟─53285b2c-3704-4adc-9bdd-9dea619b5140
# ╠═2c554897-e486-418f-9cd4-a3bb3f00fc31
# ╠═19cc5ec6-f7cb-4736-98c3-498d91eb2414
# ╠═74272021-5f45-4161-bea4-8386697de280
# ╠═587d934e-161e-403c-831b-a61f313f5ecc
# ╠═7ba7dd56-4ee9-4c5c-97a8-5f99a9ba507e
# ╠═880d1ae8-5a5a-463a-8e35-a73bdc40208a
# ╟─c2fcf437-bd86-4b24-8c77-29702cf5373d


text:
using Pkg
 for p in ("ArgParse", "Knet", "AutoGrad", "Gym")
    if !haskey(Pkg.installed(),p)
        Pkg.add(p)
        if p == "Gym"
            ENV["GYM_ENVS"] = "atari:algorithmic:box2d:classic_control"
            Pkg.build("Gym")
        end
    end
end

"""
julia actor_critic_discrete.jl

This example implements the online version of Actor-Critic algorithm.
It is a variant of REINFORCE algorithm. The TD-error is
used as the advantage function.
"""
module ACTOR_CRITIC

using Gym, ArgParse, Knet, AutoGrad, Statistics

function predict(w, x, name; nh=1)
    inp = x
    for i=1:nh
        inp = relu.(w[name*"w_$i"] * inp .+ w[name*"b_$i"])
    end
    out = w[name*"w_out"] * inp .+ w[name*"b_out"]
    return out
end

function sample_action(linear)
    linear = Array(linear)
    probs = vec(exp.(linear) ./ sum(exp.(linear); dims=1))
    c_probs = cumsum(probs)
    return findmax(c_probs .> rand())[2]
end

AutoGrad.@zerograd sample_action(linear)

function actor(w, ob; nh=1)
    linear = predict(w, ob, "actor";nh=nh)
    action = sample_action(linear)
    return action, linear
end

critic(w, ob;nh=1) = predict(w, ob, "critic"; nh=nh)

function loss(w, ob, t, env, o; output_of_step=nothing)
    action, linear = actor(w, ob; nh=length(o["actor"]))
    v_s = critic(w, ob; nh=length(o["critic"]))

    ob, reward, done, _ = step!(env, action-1)
    ob = convert(o["atype"], reshape(ob, size(ob, 1), 1))
    push!(output_of_step, (ob, reward, done))

    v_sp1 = done ? Float32(0.0) : critic(w, ob; nh=length(o["critic"]))

    reward = Float32(reward)

    δ = reward .+ o["gamma"] .* v_sp1 .- v_s

    critic_loss = sum(δ .* δ)# size of δ is (1,1)
    actor_loss = sum(-logp(linear; dims=1)[action] .*  AutoGrad.getval(δ)) / size(linear, 2)

    return actor_loss + critic_loss
end

lossgradient = grad(loss)

function play_episode!(w, opts, env, o)
    ob = reset!(env)
    total = 0
    ob = convert(o["atype"], reshape(ob, size(ob, 1), 1))
    for t=1:env.spec.max_episode_steps
        output = Any[]
        g = lossgradient(w, ob, t, env, o; output_of_step=output)
        update!(w, g, opts)

        #converted next state, reward and isTerminal
        ob, reward, done = output[1]

        total += reward
        o["render"] && render(env)

        done && break
    end
    return total
end

function init_weights(name, input, hiddens, output, atype)
    w = Dict()
    inp = input
    for i=1:length(hiddens)
        w[name*"w_$i"] = 0.01*randn(hiddens[i], inp)
        w[name*"b_$i"] = zeros(hiddens[i])
        inp = hiddens[i]
    end

    w[name*"w_out"] = 0.01*randn(output, hiddens[end])
    w[name*"b_out"] = zeros(output, 1)

    for k in keys(w)
        w[k] = convert(atype, w[k])
    end

    return w
end

function main(args=ARGS)
    s = ArgParseSettings()
    s.description="(c) Ozan Arkan Can, 2018. Demonstration of the online Actor-Critic algorithm."
    @add_arg_table s begin
        ("--env_id"; default="CartPole-v0"; help="environment name")
        ("--actor"; nargs='+'; arg_type=Int; default=[100]; help="number of hiddens for the actor")
        ("--critic"; nargs='+'; arg_type=Int; default=[100]; help="number of hiddens for the critic")
        ("--episodes"; arg_type=Int; default=20; help="number of episodes")
        ("--gamma"; arg_type=Float32; default=Float32(0.99); help="doscount factor")
        ("--lr"; arg_type=Float32; default=Float32(0.01); help="learning rate")
        ("--render"; help = "render the environment"; action = :store_true)
        ("--usegpu"; action=:store_true; help="use GPU or not")
    end

    Knet.seed!(12345)
    isa(args, AbstractString) && (args=split(args))
    if in("--help", args) || in("-h", args)
        ArgParse.show_help(s; exit_when_done=false)
        return
    end

    o = parse_args(args, s)
    o["atype"] = !o["usegpu"] ? Array{Float32} : KnetArray{Float32}

    env = GymEnv(o["env_id"])
    seed!(env, 12345)

    INPUT = env.observation_space.shape[1]
    OUTPUT = env.action_space.n

    w = init_weights("actor", INPUT, o["actor"], OUTPUT, o["atype"])
    merge!(w, init_weights("critic", INPUT, o["critic"], 1, o["atype"]))
    opts = Dict()
    for k in keys(w)
        opts[k] = Adam(lr=o["lr"])
    end

    for i=1:o["episodes"]
        total = play_episode!(w, opts, env, o)
        println("episode $i , total rewards: $total")
    end
end

PROGRAM_FILE=="actor_critic_discrete.jl" && main(ARGS)

end


text:
module MCParticles
using Distributions

using StatsFuns: logsumexp

export Particle, StaticParticle, DynamicParticle
export weights, location, logdensity
export pequal
export reweight
export mutate

include("helper-functions.jl")
include("particle-constructor.jl")
include("particle-equal.jl")
include("particle-array-accessor.jl")
include("particle-reweight.jl")
include("particle-mutate.jl")

end


text:
export ModulatedDevice
export solve, setup_Δϵᵣ!

mutable struct ModulatedDevice{D} <: AbstractDevice{D}
    grid::Grid{D}
    ϵᵣ::Array{Complex}
    Δϵᵣ::Array{Complex}
    src::Array{Complex}
    ω::AbstractVector{Float}
    Ω::Float
    nsidebands::Integer
    sharedpml::Bool
    modes::Array{Mode}
end

"    ModulatedDevice(grid::Grid, ω::AbstractVector{Float}, Ω::Float, nsidebands::Integer; sharedpml::Bool=true)"
function ModulatedDevice(grid::Grid, ω::AbstractVector{Float}, Ω::Float, nsidebands::Integer; sharedpml::Bool=true)
    ϵᵣ = ones(Complex, size(grid));
    Δϵᵣ = zeros(Complex, size(grid));
    src = zeros(Complex, size(grid));
    return ModulatedDevice(grid, ϵᵣ, Δϵᵣ, src, ω, Ω, nsidebands, sharedpml, Array{Mode}(undef, 0))
end

"    ModulatedDevice(grid::Grid, ω::Float, Ω::Float, nsidebands::Integer; sharedpml::Bool=true)"
ModulatedDevice(grid::Grid, ω::Float, Ω::Float, nsidebands::Integer; sharedpml::Bool=true) =
 ModulatedDevice(grid, [ω], Ω, nsidebands, sharedpml=sharedpml);

"    setup_Δϵᵣ!(d::ModulatedDevice, shapes::AbstractVector{<:Shape})"
setup_Δϵᵣ!(d::ModulatedDevice, shapes::AbstractVector{<:Shape}) = _compose_shapes!(d.Δϵᵣ, d.grid, shapes)

"    setup_Δϵᵣ!(d::ModulatedDevice, region, value)"
setup_Δϵᵣ!(d::ModulatedDevice, region, value) = _mask_values!(d.Δϵᵣ, d.grid, region, value)

"    solve(d::ModulatedDevice)"
function solve(d::ModulatedDevice)
    (ϵ₀, μ₀, c₀) = normalize_parameters(d);

    Nω = length(d.ω);
    nsidebands = d.nsidebands;
    nfrequencies = 2*nsidebands+1;
    n = -nsidebands:1:nsidebands;
    Ω = d.Ω;

    fields = Array{FieldTM}(undef, Nω, nfrequencies);

    for i in eachindex(d.ω)
        @info "Frequency: $i/$Nω"
        ω = d.ω[i];
        ωn = ω .+ Ω*n;

        length(d.modes) > 0 && ( d.src=zeros(Complex, size(d.grid)) ) # Only reset if we are using modes
        for mode in d.modes
            # TODO: handle the polarization here
            setup_mode!(d, TM, ω, mode.neff, mode.pt, mode.dir, mode.width);
        end

        Tϵ = sparse(Diagonal(ϵ₀*d.ϵᵣ[:])); #TODO: check reshape vs [:]
        TΔϵ = sparse(Diagonal(ϵ₀*d.Δϵᵣ[:])); #TODO: check reshape vs [:]

        # Construct derivates
        δxb = δ(x̂, Backward, d.grid);
        δxf = δ(x̂, Forward,  d.grid);
        δyb = δ(ŷ, Backward, d.grid);
        δyf = δ(ŷ, Forward,  d.grid);

        # Reshape Mz into a vector
        b0 = 1im*ω*d.src[:]; #TODO: check reshape vs [:]
        b = zeros(ComplexF64, length(d.grid)*nfrequencies, 1);
        b[(nsidebands*length(d.grid))+1:(nsidebands+1)*length(d.grid), 1] = b0;

        @info "Calculating: system matrix"
        @info "Sidebands (freqs): $nsidebands ($nfrequencies)"

        As = Array{SparseMatrixCSC}(undef, nfrequencies);
        Sxf = Array{SparseMatrixCSC}(undef, d.sharedpml ? 1 : nfrequencies);
        Sxb = Array{SparseMatrixCSC}(undef, d.sharedpml ? 1 : nfrequencies);
        Syf = Array{SparseMatrixCSC}(undef, d.sharedpml ? 1 : nfrequencies);
        Syb = Array{SparseMatrixCSC}(undef, d.sharedpml ? 1 : nfrequencies);
        if d.sharedpml
            @info "PML: shared"
            (Sxf[1], Sxb[1], Syf[1], Syb[1]) = S_create(d.grid, ω);
            A1 = Sxb[1]*δxb*1/μ₀*Sxf[1]*δxf + Syb[1]*δyb*1/μ₀*Syf[1]*δyf;
            for j = 1:nfrequencies
                As[j] = A1 + ωn[j]^2*Tϵ;
            end
        else
            @info "PML: per-frequency"
            for j = 1:nfrequencies
                (Sxf[j], Sxb[j], Syf[j], Syb[j]) = S_create(d.grid, ωn[j]);
                As[j] = Sxb[j]*δxb*1/μ₀*Sxf[j]*δxf + Syb[j]*δyb*1/μ₀*Syf[j]*δyf + ωn[j]^2*Tϵ;
            end
        end
        if nsidebands > 0
            @info "Calculating: coupling matrix"
            stencil_Cp = spdiagm( 1 => 0.5*ωn[1:end-1].^2);
            Cp = kron(stencil_Cp, conj(TΔϵ));
            stencil_Cm = spdiagm(-1 => 0.5*ωn[2:end].^2);
            Cm = kron(stencil_Cm, TΔϵ);

            @info "Calculating: total matrix"
            A = blockdiag(As...) + Cp + Cm;
        else
            A = As[1];
        end

        ez = dolinearsolve(A, b, CNSym);

        @info "Processing results"

        for j = 1:nfrequencies
            ezi = ez[(j-1)*length(d.grid)+1:j*length(d.grid)];
            hxi = -1/1im/ωn[j]/μ₀*Syf[d.sharedpml ? 1 : j]*δyf*ezi;
            hyi = 1/1im/ωn[j]/μ₀*Sxf[d.sharedpml ? 1 : j]*δxf*ezi;
            fields[i, j] = FieldTM(d.grid, ωn[j], ezi, hxi, hyi)
        end
    end

    return fields
end


text:
using Oceananigans.Solvers: solve_for_pressure!, solve_poisson_equation!, set_source_term!
using Oceananigans.Solvers: poisson_eigenvalues
using Oceananigans.Models.HydrostaticFreeSurfaceModels: _compute_w_from_continuity!
using Oceananigans.BoundaryConditions: regularize_field_boundary_conditions

function poisson_solver_instantiates(arch, grid, planner_flag)
    solver = FFTBasedPoissonSolver(arch, grid, planner_flag)
    return true  # Just making sure the FFTBasedPoissonSolver does not error/crash.
end

function random_divergent_source_term(arch, grid)
    default_bcs = FieldBoundaryConditions()
    u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
    v_bcs = regularize_field_boundary_conditions(default_bcs, grid, :v)
    w_bcs = regularize_field_boundary_conditions(default_bcs, grid, :w)

    Ru = CenterField(arch, grid, u_bcs)
    Rv = CenterField(arch, grid, v_bcs)
    Rw = CenterField(arch, grid, w_bcs)
    U = (u=Ru, v=Rv, w=Rw)

    Nx, Ny, Nz = size(grid)
    set!(Ru, rand(Nx, Ny, Nz))
    set!(Rv, rand(Nx, Ny, Nz))
    set!(Rw, rand(Nx, Ny, Nz))

    fill_halo_regions!(Ru, arch)
    fill_halo_regions!(Rv, arch)
    fill_halo_regions!(Rw, arch)

    # Compute the right hand side R = ∇⋅U
    ArrayType = array_type(arch)
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, U.u.data, U.v.data, U.w.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    return R, U
end

function random_divergence_free_source_term(arch, grid)
    default_bcs = FieldBoundaryConditions()
    u_bcs = regularize_field_boundary_conditions(default_bcs, grid, :u)
    v_bcs = regularize_field_boundary_conditions(default_bcs, grid, :v)
    w_bcs = regularize_field_boundary_conditions(default_bcs, grid, :w)

    # Random right hand side
    Ru = CenterField(arch, grid, u_bcs)
    Rv = CenterField(arch, grid, v_bcs)
    Rw = CenterField(arch, grid, w_bcs)
    U = (u=Ru, v=Rv, w=Rw)

    Nx, Ny, Nz = size(grid)
    set!(Ru, rand(Nx, Ny, Nz))
    set!(Rv, rand(Nx, Ny, Nz))
    set!(Rw, zeros(Nx, Ny, Nz))

    fill_halo_regions!(Ru, arch, nothing, nothing)
    fill_halo_regions!(Rv, arch, nothing, nothing)
    fill_halo_regions!(Rw, arch, nothing, nothing)

    event = launch!(arch, grid, :xy, _compute_w_from_continuity!, U, grid,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    fill_halo_regions!(Rw, arch, nothing, nothing)

    # Compute the right hand side R = ∇⋅U
    ArrayType = array_type(arch)
    R = zeros(Nx, Ny, Nz) |> ArrayType
    event = launch!(arch, grid, :xyz, divergence!, grid, Ru.data, Rv.data, Rw.data, R,
                    dependencies=Event(device(arch)))
    wait(device(arch), event)

    return R
end

#####
##### Regular rectilinear grid Poisson solver
#####

function divergence_free_poisson_solution(arch, grid, planner_flag=FFTW.MEASURE)
    ArrayType = array_type(arch)
    FT = eltype(grid)

    solver = FFTBasedPoissonSolver(arch, grid, planner_flag)
    R, U = random_divergent_source_term(arch, grid)

    p_bcs = FieldBoundaryConditions(grid, (Center, Center, Center))
    ϕ   = CenterField(arch, grid, p_bcs)  # "kinematic pressure"
    ∇²ϕ = CenterField(arch, grid, p_bcs)

    # Using Δt = 1 but it doesn't matter since velocities = 0.
    solve_for_pressure!(ϕ.data, solver, arch, grid, 1, datatuple(U))

    compute_∇²!(∇²ϕ, ϕ, arch, grid)

    return CUDA.@allowscalar interior(∇²ϕ) ≈ R
end

#####
##### Test that Poisson solver error converges as error ~ N⁻²
#####

ψ(::Type{Bounded}, n, x) = cos(n*x/2)
ψ(::Type{Periodic}, n, x) = cos(n*x)

k²(::Type{Bounded}, n) = (n/2)^2
k²(::Type{Periodic}, n) = n^2

function analytical_poisson_solver_test(arch, N, topo; FT=Float64, mode=1)
    grid = RegularRectilinearGrid(FT, topology=topo, size=(N, N, N), x=(0, 2π), y=(0, 2π), z=(0, 2π))
    solver = FFTBasedPoissonSolver(arch, grid)

    xC, yC, zC = nodes((Center, Center, Center), grid, reshape=true)

    TX, TY, TZ = topology(grid)
    Ψ(x, y, z) = ψ(TX, mode, x) * ψ(TY, mode, y) * ψ(TZ, mode, z)
    f(x, y, z) = -(k²(TX, mode) + k²(TY, mode) + k²(TZ, mode)) * Ψ(x, y, z)

    solver.storage .= convert(array_type(arch), f.(xC, yC, zC))

    solve_poisson_equation!(solver)

    ϕ = real(Array(solver.storage))

    L¹_error = mean(abs, ϕ - Ψ.(xC, yC, zC))

    return L¹_error
end

function poisson_solver_convergence(arch, topo, N¹, N²; FT=Float64, mode=1)
    error¹ = analytical_poisson_solver_test(arch, N¹, topo; FT, mode)
    error² = analytical_poisson_solver_test(arch, N², topo; FT, mode)

    rate = log(error¹ / error²) / log(N² / N¹)

    TX, TY, TZ = topo
    @info "Convergence of L¹-normed error, $(typeof(arch)), $FT, ($(N¹)³ -> $(N²)³), topology=($TX, $TY, $TZ): $rate"

    return isapprox(rate, 2, rtol=5e-3)
end

#####
##### Vertically stretched Poisson solver
#####

get_grid_size(TX, TY, TZ, Nx, Ny, Nz) = (Nx, Ny, Nz)
get_grid_size(::Type{Flat}, TY, TZ, Nx, Ny, Nz) = (Ny, Nz)
get_grid_size(TX, ::Type{Flat}, TZ, Nx, Ny, Nz) = (Nx, Nz)

get_xy_interval_kwargs(TX, TY, TZ) = (x=(0, 1), y=(0, 1))
get_xy_interval_kwargs(TX, ::Type{Flat}, TZ) = (x=(0, 1),)
get_xy_interval_kwargs(::Type{Flat}, TY, TZ) = (y=(0, 1),)

function vertically_stretched_poisson_solver_correct_answer(FT, arch, topo, Nx, Ny, zF)
    Nz = length(zF) - 1
    sz = get_grid_size(topo..., Nx, Ny, Nz)
    xy_intervals = get_xy_interval_kwargs(topo...)
    vs_grid = VerticallyStretchedRectilinearGrid(FT; architecture=arch, topology=topo, size=sz, z_faces=zF, xy_intervals...)
    solver = FourierTridiagonalPoissonSolver(arch, vs_grid)

    p_bcs = FieldBoundaryConditions(vs_grid, (Center, Center, Center))
    ϕ   = CenterField(arch, vs_grid, p_bcs)  # "kinematic pressure"
    ∇²ϕ = CenterField(arch, vs_grid, p_bcs)

    R = random_divergence_free_source_term(arch, vs_grid)

    set_source_term!(solver, R)
    solve_poisson_equation!(solver)

    # interior(ϕ) = solution(solver) or solution!(interior(ϕ), solver)
    CUDA.@allowscalar interior(ϕ) .= real.(solver.storage)
    compute_∇²!(∇²ϕ, ϕ, arch, vs_grid)

    return CUDA.@allowscalar interior(∇²ϕ) ≈ R
end

#####
##### Run pressure solver tests
#####

PB = (Periodic, Bounded)
topos = collect(Iterators.product(PB, PB, PB))[:]

two_dimensional_topologies = [(Flat,     Bounded,  Bounded),
                              (Bounded,  Flat,     Bounded),
                              (Bounded,  Bounded,  Flat),
                              (Flat,     Periodic, Bounded),
                              (Periodic, Flat,     Bounded),
                              (Periodic, Bounded,  Flat)]

@testset "Poisson solvers" begin
    @info "Testing Poisson solvers..."

    for arch in archs
        @testset "Poisson solver instantiation [$(typeof(arch))]" begin
            @info "  Testing Poisson solver instantiation [$(typeof(arch))]..."
            for FT in float_types

                grids_3d = [RegularRectilinearGrid(FT, size=(2, 2, 2), extent=(1, 1, 1)),
                            RegularRectilinearGrid(FT, size=(1, 2, 2), extent=(1, 1, 1)),
                            RegularRectilinearGrid(FT, size=(2, 1, 2), extent=(1, 1, 1)),
                            RegularRectilinearGrid(FT, size=(2, 2, 1), extent=(1, 1, 1))]

                grids_2d = [RegularRectilinearGrid(FT, size=(2, 2), extent=(1, 1), topology=topo)
                            for topo in two_dimensional_topologies]

                grids = []
                push!(grids, grids_3d..., grids_2d...)

                for grid in grids
                    @test poisson_solver_instantiates(arch, grid, FFTW.ESTIMATE)
                    @test poisson_solver_instantiates(arch, grid, FFTW.MEASURE)
                end
            end
        end

        @testset "Divergence-free solution [$(typeof(arch))]" begin
            @info "  Testing divergence-free solution [$(typeof(arch))]..."

            for topo in topos
                for N in [7, 16]

                    grids_3d = [RegularRectilinearGrid(topology=topo, size=(N, N, N), extent=(1, 1, 1)),
                                RegularRectilinearGrid(topology=topo, size=(1, N, N), extent=(1, 1, 1)),
                                RegularRectilinearGrid(topology=topo, size=(N, 1, N), extent=(1, 1, 1)),
                                RegularRectilinearGrid(topology=topo, size=(N, N, 1), extent=(1, 1, 1))]

                    grids_2d = [RegularRectilinearGrid(size=(N, N), extent=(1, 1), topology=topo)
                                for topo in two_dimensional_topologies]

                    grids = []
                    push!(grids, grids_3d..., grids_2d...)

                    for grid in grids
                        N == 7 && @info "    Testing $(topology(grid)) topology on square grids [$(typeof(arch))]..."
                        @test divergence_free_poisson_solution(arch, grid)
                    end
                end
            end

            Ns = [11, 16]
            for topo in topos
                @info "    Testing $topo topology on rectangular grids with even and prime sizes [$(typeof(arch))]..."
                for Nx in Ns, Ny in Ns, Nz in Ns
                    grid = RegularRectilinearGrid(topology=topo, size=(Nx, Ny, Nz), extent=(1, 1, 1))
                    @test divergence_free_poisson_solution(arch, grid)
                end
            end

            # Do a couple at Float32 (since its too expensive to repeat all tests...)
            Float32_grids = [RegularRectilinearGrid(Float32, topology=(Periodic, Bounded, Bounded), size=(16, 16, 16), extent=(1, 1, 1)),
                             RegularRectilinearGrid(Float32, topology=(Bounded, Bounded, Periodic), size=(7, 11, 13), extent=(1, 1, 1))]

            for grid in Float32_grids
                @test divergence_free_poisson_solution(arch, grid)
            end
        end

        @testset "Convergence to analytic solution [$(typeof(arch))]" begin
            @info "  Testing convergence to analytic solution [$(typeof(arch))]..."
            for topo in topos
                @test poisson_solver_convergence(arch, topo, 2^6, 2^7)
                @test poisson_solver_convergence(arch, topo, 67, 131, mode=2)
            end
        end
    end

    # Vertically stretched topologies to test.
    vs_topos = [
        (Periodic, Periodic, Bounded),
        (Periodic, Bounded,  Bounded),
        (Bounded,  Periodic, Bounded),
        (Bounded,  Bounded,  Bounded),
        (Flat,     Bounded,  Bounded),
        (Flat,     Periodic, Bounded),
        (Bounded,  Flat,     Bounded),
        (Periodic, Flat,     Bounded)
    ]

    for arch in archs, topo in vs_topos
        @testset "Vertically stretched Poisson solver [FACR, $(typeof(arch)), $topo]" begin
            @info "  Testing vertically stretched Poisson solver [FACR, $(typeof(arch)), $topo]..."

            @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8, 8, 1:8)
            @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 7, 7, 1:7)
            @test vertically_stretched_poisson_solver_correct_answer(Float32, arch, topo, 8, 8, 1:8)

            zF_even = [1, 2, 4, 7, 11, 16, 22, 29, 37]      # Nz = 8
            zF_odd  = [1, 2, 4, 7, 11, 16, 22, 29, 37, 51]  # Nz = 9

            for zF in [zF_even, zF_odd]
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8,  8, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 16, 8, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8, 16, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 8, 11, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 5,  8, zF)
                @test vertically_stretched_poisson_solver_correct_answer(Float64, arch, topo, 7, 13, zF)
            end
        end
    end
end


text:
using LinearAlgebra

import Gurobi
const GRBENV = Gurobi.Env()

include("LPBenchmark.jl")

function run_gurobi(fname::String, fsol::String="")

    grb = Gurobi.Optimizer(GRBENV)
    bridged = MOIB.full_bridge_optimizer(grb, Float64)

    # Read file
    load_problem!(bridged, fname)

    # MOI parameters
    MOI.set(bridged, MOI.Silent(), false)
    MOI.set(bridged, MOI.TimeLimitSec(), 10_000)
    MOI.set(bridged, MOI.RawParameter("Threads"), 1)
    MOI.set(bridged, MOI.RawParameter("Crossover"), 0)  # No crossover
    MOI.set(bridged, MOI.RawParameter("Method"), 2)  # barrier

    # Solve instance
    run(bridged)
    return nothing
end

if abspath(PROGRAM_FILE) == @__FILE__
    run_gurobi(ARGS[1])
end

max_stars_repo_name: mtanneau/LPBenchmarks

max_issues_repo_name: mtanneau/LPBenchmarks

max_forks_repo_name: mtanneau/LPBenchmarks


text:
abstract type Dopping
end

struct Intrinsic <: Dopping
end
Nd(D::Intrinsic,T,Ef,Ec,Ev) = 0.0/Dunit^3
Na(D::Intrinsic,T,Ef,Ec,Ev) = 0.0/Dunit^3
Nd(D::Intrinsic) = 0.0/Dunit^3
Na(D::Intrinsic) = 0.0/Dunit^3
Ed(D::Dopping) = 0.0*Eunit
struct NDoppingC <: Dopping
    Nd::Number
end
Nd(D::NDoppingC,T,Ef,Ec,Ev) = D.Nd
Na(D::NDoppingC,T,Ef,Ec,Ev) = 0/Dunit^3
Nd(D::NDoppingC) = D.Nd
Na(D::NDoppingC) = 0/Dunit^3
function NSilicon(N)
    S = Silicon()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,NDoppingC(N),S.χ,S.T0)
end
function NSiliconNU(N)
    S = SiliconNU()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,NDoppingC(N),S.χ,S.T0)
end

struct PDoppingC <: Dopping
    Na::Number
end
Nd(D::PDoppingC,T,Ef,Ec,Ev) = 0.0/Dunit^3
Na(D::PDoppingC,T,Ef,Ec,Ev) = D.Na
Nd(D::PDoppingC) = 0.0/Dunit^3
Na(D::PDoppingC) = D.Na
function PSilicon(N)
    S = Silicon()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,PDoppingC(N),S.χ,S.T0)
end
function PSiliconNU(N)
    S = SiliconNU()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,PDoppingC(N),S.χ,S.T0)
end
struct NDopping <: Dopping
    Nd::Number
    Ed::Number
end
Nd(D::NDopping,T,Ef,Ec,Ev) = D.Nd/(1+2*exp((Ef-(Ec-D.Ed))/kb/T))
Na(D::NDopping,T,Ef,Ec,Ev) = 0.0/Dunit^3
Nd(D::NDopping) = D.Nd
Na(D::NDopping) = 0.0/Dunit^3
Ed(D::NDopping) = D.Ed
function NSilicon(Nd,Ed)
    S = Silicon()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,NDopping(Nd,Ed),S.χ,S.T0)
end
function NSiliconNU(Nd,Ed)
    S = SiliconNU()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,NDopping(Nd,Ed),S.χ,S.T0)
end
struct PDopping <: Dopping
    Na::Number
    Ea::Number
end
Na(D::PDopping,T,Ef,Ec,Ev) = D.Na/(1+4*exp((Ev+D.Ea-Ef)/kb/T))
Nd(D::PDopping,T,Ef,Ec,Ev) = 0.0/Dunit^3
Na(D::PDopping) = D.Na
Nd(D::PDopping) = 0.0/Dunit^3
Ed(D::PDopping) = D.Ea
function PSilicon(Na,Ea)
    S = Silicon()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,PDopping(Na,Ea),S.χ,S.T0)
end
function PSiliconNU(Na,Ea)
    S = SiliconNU()
    SemiconductorMaterial(S.ϵₘ,S.Nv300,S.Nc300,S.Eg,PDopping(Na,Ea),S.χ,S.T0)
end


text:
using Base: LogicalIndex
@testset "event" begin

    for platform in testplatforms

        #Execute add agent without problems
        @test_nowarn begin
            m = @agent(
                0,

                [active,id1]::Identity,
                loc::Local,

                UpdateLocal =
                begin
                    addAgent(
                        active = 1,
                        id1 = 1,
                        loc = 0.5
                    )
                end
            )

            m = compile(m,platform=platform)
        end

        #Linear growth
        m = @agent(
            0,

            [active,id1]::Identity,
            loc::Local,

            UpdateLocal =
            begin
                if active == 1
                    addAgent(
                        active = 1,
                        id1 = 4,
                        loc = 67.
                    )
                    active = 0
                end
            end
        )

        m = compile(m,platform=platform,debug=false)
        # println(m.program)

        @test_nowarn begin

            com = Community(m,N=1)
            n = com.N

            com.id1 .= 4
            com.active .= 1
            com.loc .= 67.

            comt = m.evolve(com,dt=1,tMax=10,nMax=11)

            if comt[end].N != 11 error() end
            if !prod(comt[end].loc .== 67.) error() end
            if !prod(comt[end].id1 .== 4) error() end
            if !any(comt[end].id .== 1:11) error() end

        end

        #Execute Death without problems
        @test_nowarn begin
            m = @agent(
                0,

                [active,id1]::Identity,
                loc::Local,

                UpdateLocal = begin
                    if loc < t
                        removeAgent()
                    end
                end
            )

            m = compile(m,platform=platform)
        end

        m = @agent(
            1,

            [active,id1]::Identity,
            [loc]::Local,

            UpdateLocal = begin
                id1 = 3
                x = x
                if loc < t
                    removeAgent()
                end
            end,
        )

        m = compile(m,platform=platform)
        # println(m.program)

        @test_nowarn begin

            com = Community(m,N=100)
            n = com.N

            com.loc .= 0:99
            com.x .= 1:100

            comt = m.evolve(com,dt=1,tMax=105,nMax=100)

            if !any(comt.N[2:102] .== Array(100:-1:0)) error() end
            if !any(comt.N[102:end] .== 0) error() end
            for i in 1:100
                if !Bool(prod(comt[i].id .== comt[i].x)) error() end
            end

        end

        #Add and remove agents at the same time
        m = @agent(
            0,

            [active,id1]::Identity,
            [loc]::Local,

            UpdateLocal = begin
                if active == 1
                    addAgent(id1 = 4,
                            active = 1,
                            loc = 67
                            )
                    removeAgent()
                end
            end,
        )

        m = compile(m,platform=platform)
        # println(prettify(m.program))

        @test_nowarn begin

            com = Community(m,N=1)
            n = com.N

            com.id1 .= 4
            com.loc .= 67

            comt = m.evolve(com,dt=1,tMax=105,nMax=100)

            if any(comt.N .!= 1) error() end
            if any([comt.id1[i,i] for i in 1:105] .!= 4) error() end
            if any([comt.loc[i,i] for i in 1:105] .!= 67) error() end
            if any(comt[end].id .!= 106) error() end
        end

    end

end


text:
@testset "Triangulation: " begin
    P1 = SMatrix{3, 4, Float64}(
        0.999701, -0.0171452, 0.0174497,
        0.0174497, 0.999695, -0.0174524,
        -0.017145, 0.0177517, 0.999695,
        -500, -100, -100,
    )
    P2 = SMatrix{3, 4, Float64}(
        0.99969, 0.0177543, -0.0174497,
        -0.0174497, 0.999695, 0.0174524,
        0.0177543, -0.0171425, 0.999695,
        500, -100, -100,
    )
    K = SMatrix{3, 3, Float64}(
        7291.67, 0, 0,
        0, 7291.67, 0,
        639.5, 511.5, 1,
    )

    P1 = K * P1
    P2 = K * P2

    p1 = SVector{2, Float64}(146, 642.288)
    p2 = SVector{2, Float64}(1137.31, 385.201)

    t = triangulate(p1, p2, P1, P2)
    t /= t[4]
    @test isapprox(t, SVector{4, Float64}(0, 100, 10000, 1); atol=1e-1)
end

@testset "Triangulation: Horizontal setup" begin
    P1 = SMatrix{3, 4, Float64}(I)
    P2 = hcat(SMatrix{3, 3, Float64}(I), SVector{3, Float64}(-1000, 0, 0))
    K = SMatrix{3, 3, Float64}(
        7291.67, 0, 0,
        0, 7291.67, 0,
        639.5, 511.5, 1,
    )

    P1 = K * P1
    P2 = K * P2

    p1 = SVector{2, Float64}(1004.08, 511.5)
    p2 = SVector{2, Float64}(274.917, 511.5)

    t = triangulate(p1, p2, P1, P2)
    t /= t[4]
    @test isapprox(t, SVector{4, Float64}(500, 0, 10000, 1); atol=1e-1)
end


text:
VERSION >= v"0.4.0-dev+6521" && __precompile__(true)

module MarketModels

using TimeSeries, Distributions

# dev use only
using MarketData, TimeSeriesQueries, GaussianMixtures

export simulate
export rday, ryear

###### include ##################

include("sim.jl")
include("const.jl")

end


text:
module OutputALASmod

using OutputLSmod

type OutputALAS

 solved::Int64 #renvoi 0 si réussi
 solve_message::AbstractString #message de sortie
 xtab::Array #tableau des itérés
 dtab::Array #tableau des directions
 rhotab::Array{Float64,1}
 steptab::Array{Float64,1}

 realisabilite::Array{Float64,1}
 dualrealisabilite::Array{Float64,1}
 obj::Array{Float64,1}

 inner_output_linesearch::Any #en fait une liste de OutputLS

end

#Initialisation
function OutputALAS(x0::Vector,dj::Vector,feas::Float64,dualfeas::Float64,rho::Vector,ht::Float64)
 return OutputALAS(0,"Success",x0,dj,[norm(rho,Inf)],[],[feas],[dualfeas],[ht],[])
end

#Mise à jour
function Update(oa::OutputALAS,xk::Vector,rho::Vector,feas::Float64,dualfeas::Float64,dj::Vector,step::Float64,outputls::OutputLSmod.OutputLS,ht::Float64,n::Int64)

 oa.realisabilite=[oa.realisabilite;feas]
 oa.dualrealisabilite=[oa.dualrealisabilite;dualfeas]
 oa.obj=[oa.obj;ht]
 oa.rhotab=[oa.rhotab;norm(rho,Inf)]
 oa.steptab=[oa.steptab;step]
 oa.xtab=[oa.xtab xk]
 oa.dtab=[oa.dtab dj]
 oa.inner_output_linesearch=[oa.inner_output_linesearch;outputls]

Print(oa,n,2,j=length(oa.steptab))

 return oa
end

#Print
function Print(oa::OutputALAS,n::Int64,verbose::Int64;j::Int64=1)
 if verbose>1
  nb_comp=Int64(0.5*(length(oa.xtab)/length(oa.realisabilite)-n))
  while j <= length(oa.steptab)
   #println("k :",j," xjk=",oa.xtab[1:n,j+1]," sg=",oa.xtab[n+1:n+nb_comp,j+1]," sh=",oa.xtab[n+nb_comp+1:n+2*nb_comp,j+1],"\n\|c(xjk)\|=",oa.realisabilite[j+1]," JL(xjk)=",oa.dualrealisabilite[j+1]," rho= ",oa.rhotab[j+1]," \nalpha=",oa.steptab[j]," d=",oa.dtab[1:n+2*nb_comp,j+1])
   if n<=4
    if nb_comp>0
     #print_with_color(:blue, "	k: $j xjk=$(oa.xtab[1:n,j+1]) sg=$(oa.xtab[n+1:n+nb_comp,j+1]) sh=$(oa.xtab[n+nb_comp+1:n+2*nb_comp,j+1]) \n\|c(xjk)\|=$(oa.realisabilite[j+1]) JL(xjk)=$(oa.dualrealisabilite[j+1]) rho=$(oa.rhotab[j+1]) \nalpha=$(oa.steptab[j]) d=$(oa.dtab[1:n+2*nb_comp,j+1]) fpen=$(oa.obj[j]) \n")
     print_with_color(:blue, "	k: $j \|xjk\|=$(norm(oa.xtab[1:n,j+1],Inf)) \|c(xjk)\|=$(oa.realisabilite[j+1]) L'(xjk)=$(oa.dualrealisabilite[j+1]) rho=$(oa.rhotab[j+1]) \nalpha=$(oa.steptab[j]) \|d\|=$(norm(oa.dtab[1:n+2*nb_comp,j+1],Inf)) fpen=$(oa.obj[j]) \n")
    else
     print_with_color(:blue, "	k: $j xjk=$(oa.xtab[1:n,j+1]) \n\|c(xjk)\|=$(oa.realisabilite[j+1]) JL(xjk)=$(oa.dualrealisabilite[j+1]) rho=$(oa.rhotab[j+1]) \nalpha=$(oa.steptab[j]) d=$(oa.dtab[1:n+2*nb_comp,j+1]) fpen=$(oa.obj[j]) \n")
    end
   else
    if nb_comp>0
     #print_with_color(:blue, "	k: $j xjk=[$(oa.xtab[1,j+1])..$(oa.xtab[n,j+1])] sg=[$(oa.xtab[n+1,j+1])..$(oa.xtab[n+nb_comp,j+1])] sh=[$(oa.xtab[n+nb_comp+1,j+1])..$(oa.xtab[n+2*nb_comp,j+1])] \n\|c(xjk)\|=$(oa.realisabilite[j+1]) JL(xjk)=$(oa.dualrealisabilite[j+1]) rho=$(oa.rhotab[j+1]) \nalpha=$(oa.steptab[j]) d=(oa.dtab[1:n+2*nb_comp,j+1]) fpen=$(oa.obj[j]) \n")
     print_with_color(:blue, "	k: $j \|xjk\|=$(norm(oa.xtab[1:n,j+1],Inf)) \|c(xjk)\|=$(oa.realisabilite[j+1]) L'(xjk)=$(oa.dualrealisabilite[j+1]) rho=$(oa.rhotab[j+1]) \nalpha=$(oa.steptab[j]) \|d\|=$(norm(oa.dtab[1:n+2*nb_comp,j+1],Inf)) fpen=$(oa.obj[j]) \n")
    else
     print_with_color(:blue, "	k: $j xjk=[$(oa.xtab[1,j+1])..$(oa.xtab[n,j+1])] \n\|c(xjk)\|=$(oa.realisabilite[j+1]) JL(xjk)=$(oa.dualrealisabilite[j+1]) rho=$(oa.rhotab[j+1]) \nalpha=$(oa.steptab[j]) d=(oa.dtab[1:n+2*nb_comp,j+1]) fpen=$(oa.obj[j]) \n")
    end
   end
   OutputLSmod.Print(oa.inner_output_linesearch[j],n,verbose)
   j+=1
  end
  println("")
 end

 return true
end

#end of module
end


text:
"""
Package: ForecastPlots

    pacf(x,
         type = "step-real",
         lag = Integer(ceil(10*log10(length(x)))),
         alpha = (0.95,0.99);
         plot = true)

Compute the partial auto-correlation for an univariate series with an optional plot.

There are two versions; the "step" version estimates the auto-regressive parameters of an increasing model, the "real" version estimates the actual partial auto-correlation by eliminating the linear information provided by the lags. When using the default type "stepwise-real" both versions are calculated.

The distribution used to estimate the confidence intervals is an approximation of a Fisher Transformation via a Normal Distribution. There is a plot recipe for the returned object.

# Arguments
- `x`: Vector of data.
- `type` = Valid values are "stepwise", "real" and "stepwise-real".
- `lag`: Maximum number of lags.
- `alpha`: A tuple with two thresholds (t1,t2) with t1 <= t2 to plot confidence intervals. The default values are 0.95 and 0.99.
- `plot`: When true displays a plot with the results.

# Returns
A CF object

# Examples
```julia-repl
x = rand(100);
pacf(x);
pacf(x; plot=false)
20×2 Matrix{Float64}:
 -0.0475437   -0.101204
  0.0407006    0.058935
 -0.0747269    0.0366026
```
"""
function pacf(x::AbstractVector{<:Real};
              type::String = "stepwise-real",
              lag::Integer = Integer(ceil(10*log10(length(x)))),
              alpha::Tuple{Real,Real} = (0.95,0.99),
              plot::Bool = true,
              kw...)

    N = length(x)
    @assert type in  ["stepwise","real","stepwise-real"] "The options for `type` are `stepwise`, `real` and `stepwise-real`"
    @assert 1 <= lag <= N-4
    @assert length(alpha) == 2
    @assert 0.0 < alpha[1] < alpha[2] < 1.0

    # Confidence Intervals
    function fse(v::Real)::AbstractFloat
        1/sqrt(v-3)
    end

    a1 = alpha[1]
    a2 = alpha[2]
    z1 = quantile(Normal(), a1 + (1-a1)/2)
    ci1 = z1*fse(N-lag)
    z2 = quantile(Normal(), a2 + (1-a2)/2)
    ci2 = z2*fse(N-lag)
    ci = (ci1,ci2)

    call = "pacf(x"*
        "; type=\""*type*
        "\", lag="*string(lag)*
        ", alpha="*string(alpha)*")"

    if (type::String == "stepwise")
        pac = pacf_stepwise(x; lag, alpha)
        plot && display(cf_gr((cres=pac, N=N, type="pacf_stepwise", lag, alpha, ci, auto=true, kw...)))
        return pac
    end

    if (type::String == "real")
        pac = pacf_real(x; lag, alpha)
        plot && display(cf_gr((cres=pac, N=N, type="pacf_real", lag, alpha, ci, auto=true, kw...)))
        return pac
    end

    if type::String == "stepwise-real"
        pacs = pacf_stepwise(x; lag, alpha)
        pacr = pacf_real(x; lag, alpha)
        plot && display(cf_gr((cres=hcat(pacs, pacr), N=N, type="pacf_stepwise-real", lag, alpha, ci, auto=true, kw...)))
        return hcat(pacs, pacr)
    end

end

function pacf_lag(x::AbstractVector{<:Real};
                   lag::Integer = Integer(ceil(10*log10(length(x)))),
                   alpha::Tuple{Real,Real} = (0.95,0.99))

    R = vcat(cov(x,x),acf(x; type = "cov", lag, alpha, plot=false))
    p = length(R)-1

    A = zeros(p,p)
    for i in p:-1:1
        A[i,i:end] = R[1:p-i+1]
        A[i:end,i] = R[1:p-i+1]
    end
    b = R[2:end]

    lambda = 10^-36
    (A'*A+lambda*I(p))\(A'*b)

end

function pacf_stepwise(x::AbstractVector{<:Real};
                   lag::Integer = Integer(ceil(10*log10(length(x)))),
                   alpha::Tuple{Real,Real} = (0.95,0.99))

    #TODO optimize with Levinson-Durbin algorithm
    pac = []
    for i in 1:lag
        push!(pac, pacf_lag(x; lag = i)[end])
    end
    pac

end

function pacf_real(x::AbstractVector{<:Real};
             lag::Integer = Integer(ceil(10*log10(length(x)))),
             alpha::Tuple{Real,Real} = (0.95,0.99))

    # prevent singularities
    lambda = 10^-36

    N = length(x)
    A = zeros(N-lag,lag+1)
    for i in 1:(N-lag)
        A[i,:] = x[i:i+lag]
    end
    bt1 = A[:,1]
    A = A[:,2:end]
    A = hcat(A,repeat([1.0],N-lag))

    pacf_res = []
    for i in 2:lag+1
        A1i = drop(A, c=i-1)
        lsq = (A1i'*A1i+lambda*I(lag))\(A1i'*bt1)
        r1 = round.(bt1.-A1i*lsq, digits=10)

        bt0 = A[:,i-1]
        lsq = (A1i'*A1i+lambda*I(lag))\(A1i'*bt0)
        r0 = round.(bt0.-A1i*lsq, digits=10)
        push!(pacf_res,cor(r0,r1))

    end

    pacf_res
end

"""
Package: ForecastPlots

    drop(M;r,c)

Drop rows and columns from a matrix.
"""
function drop(M::AbstractMatrix; r=[], c=[])
    s = size(M)
    dr = collect(1:s[1])
    dc = collect(1:s[2])
    splice!(dr,r)
    splice!(dc,c)
    M[dr,dc]
end


text:
using OptimTestProblems
using OptimTestProblems.MultivariateProblems

import LinearAlgebra: norm

using Test

verbose = false

@testset "Bounded univariate problems" begin
    uvp = OptimTestProblems.UnivariateProblems.examples
    for (name, p) in uvp
        verbose && print_with_color(:green, "Problem: $name \n")
        for (miz, mia) in zip(p.minimizers, p.minima)
            @test p.f(miz) ≈ mia
        end
    end

end

@testset "Unconstrained multivariate problems" begin
    muvp = MultivariateProblems.UnconstrainedProblems.examples
    for (name, p) in muvp
        verbose && print_with_color(:green, "Problem: $name \n")
        soltest = !any(isnan, p.solutions)

        if startswith(name, "Penalty Function I")
            # The provided solutions are not exact
            tol = 1e-16
        else
            tol = 1e-32
        end

        f = objective(p)
        soltest && @test f(p.solutions) ≈ p.minimum

        gs = similar(p.initial_x)
        g! = gradient(p)
        g!(gs, p.solutions)
        soltest && @test norm(gs, Inf) < tol

        fg! = objective_gradient(p)
        fgs = similar(gs)
        g!(gs, p.initial_x)

        @test fg!(fgs, p.initial_x) ≈ f(p.initial_x)
        @test norm(fgs.-gs, Inf)  < eps(eltype(gs))
    end
end

@testset "Constrained multivariate problems" begin
    mcvp = MultivariateProblems.ConstrainedProblems.examples
    for (name, p) in mcvp
        verbose && print_with_color(:green, "Problem: $name \n")
        soltest = all(isfinite, p.solutions)

        f = objective(p)
        soltest && @test f(p.solutions) ≈ p.minimum

        if !isempty(p.constraintdata.lx)
            @test all(p.solutions .> p.constraintdata.lx)
        end
        if !isempty(p.constraintdata.ux)
            @test all(p.solutions .< p.constraintdata.ux)
        end

        lc = p.constraintdata.lc
        if !isempty(lc)
            c = fill(eltype(lc)(0.0), size(lc))
            p.constraintdata.c!(c, p.solutions)
            @test all(c .>= p.constraintdata.lc)
            @test all(c .<= p.constraintdata.uc)
        end

        # TODO: How do we test for optimality here?
        gs = similar(p.initial_x)
        g! = gradient(p)
        g!(gs, p.solutions)
        #soltest && @test norm(gs, Inf) < tol

        fg! = objective_gradient(p)
        fgs = similar(gs)
        g!(gs, p.initial_x)

        @test fg!(fgs, p.initial_x) ≈ f(p.initial_x)
        @test norm(fgs.-gs, Inf)  < eps(eltype(gs))
    end
end


text:
using Example
using Test

@testset "Example.jl" begin
    # Write your tests here.
    @test my_f(2, 1) == 5
    @test my_f(2, 3) == 7
end

@testset "Derivative test.jl" begin
    # Write your tests here.
    @test derivative_of_my_f(2, 1) == 2
end


text:
#Old way of organizing package, should be deprecated in the future
# Dielectrics
struct Sic <: OptProp end
struct Cbn <: OptProp end

function permittivity(material::Cbn,w) :: Complex{Float64}
    eps_fin = 4.46 + 0.0*im
    w_lo    = 2.451e14 # rad/s
    w_to    = 1.985e14 # rad/s
    gamma   = 9.934e11 # rad/s
    return eps_fin*(w^2-w_lo^2 + im*gamma*w)/(w^2-w_to^2 + im*gamma*w)
end

function permittivity(material::Sic,w) :: Complex{Float64}
    eps_fin = 6.7 + 0.0*im
    w_lo    = 1.827e14 # rad/s
    w_to    = 1.495e14 # rad/s
    gamma   = 8.971e11 # rad/s
    return eps_fin*(w^2-w_lo^2 + im*gamma*w)/(w^2-w_to^2 + im*gamma*w)
end


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

JLLWrappers.@generate_main_file_header("Dierckx")
JLLWrappers.@generate_main_file("Dierckx", UUID("cd4c43a9-7502-52ba-aa6d-59fb2a88580b"))
end  # module Dierckx_jll


text:
# Removing outliers using a StatisticalOutlierRemoval filter
# http://www.pointclouds.org/documentation/tutorials/statistical_outlier.php

using PCLCommon
using PCLIO
using PCLFilters
using PCLVisualization

pcd_file = Pkg.dir("PCLIO", "test", "data", "table_scene_lms400.pcd")
cloud = PointCloud{PointXYZ}(pcd_file)
cloud_filtered = PointCloud{PointXYZ}()

sor = StatisticalOutlierRemoval{PointXYZ}()
setInputCloud(sor, cloud)
setMeanK(sor, 50)
setStddevMulThresh(sor, 1.0)
@time filter(sor, cloud_filtered)

println("PointCloud before filtering: $(length(cloud)) data points")
println("PointCloud before filtering: $(length(cloud_filtered)) data points")

if isdefined(:vis) && vis
    viewer = PCLVisualizer("pcl visualizeer")
    red_handler = PointCloudColorHandlerCustom(cloud, 255, 0, 0)
    green_handler = PointCloudColorHandlerCustom(cloud, 0, 255, 0)
    addPointCloud(viewer, cloud, red_handler, id="cloud")
    addPointCloud(viewer, cloud_filtered, green_handler, id="cloud_filtered")
    spin(viewer)
    close(viewer)
end


text:
# Tue 13 Nov 2012 04:13:36 PM EST
#
# NeHe Tut 9 - Make some colored stars and play w/ alpha blending a bit more
#
# Q - quit
# T - turn "twinkle" on/off
# PageUp/Down - move camera closer/further away
# Up/Down - increase/decrease tilt about x-axis

# load necessary GLUT/OpenGL routines

using OpenGL
@OpenGL.version "1.0"
@OpenGL.load
using GLUT

# initialize variables

STAR_NUM = 50

type star
    r::Int
    g::Int
    b::Int
    dist::Float64
    angle::Float64
end

tempr = rand(1:256)
tempg = rand(1:256)
tempb = rand(1:256)

stars = [star(tempr,tempg,tempb,0.0,0.0)] # Julia doesn't like it when you try to initialize an empty array of
                                        # a composite type and try to fill it afterwards, so we start with a 1-element
                                        # vector and tack on values in a loop

for loop = 1:STAR_NUM-1
    tempr = rand(1:256)
    tempg = rand(1:256)
    tempb = rand(1:256)
    stars = push!(stars,star(tempr,tempg,tempb,loop/STAR_NUM*5.0,0.0))
end

global window

global tilt    = 90.0
global zoom    = -15.0
global spin    = 0.0

global twinkle = false

global tex     = Array(Uint32,1) # generating 1 texture

width          = 640
height         = 480

# load textures from images

function LoadGLTextures()
    global tex

    img, w, h = glimread(expanduser("~/.julia/GLUT/Examples/NeHe/tut9/Star.bmp"))

    glGenTextures(1,tex)
    glBindTexture(GL_TEXTURE_2D,tex[1])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, GL_LINEAR)
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_LINEAR)
    glTexImage2D(GL_TEXTURE_2D, 0, 3, w, h, 0, GL_RGB, GL_UNSIGNED_BYTE, img)
end

# function to init OpenGL context

function initGL(w::Integer,h::Integer)
    glViewport(0,0,w,h)
    LoadGLTextures()
    glClearColor(0.0, 0.0, 0.0, 0.0)
    glClearDepth(1.0)
    glShadeModel(GL_SMOOTH)

    # enable texture mapping and alpha blending
    glEnable(GL_TEXTURE_2D)
    glEnable(GL_BLEND)
    glblendfunc(GL_SRC_ALPHA, GL_ONE)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)
end

# prepare Julia equivalents of C callbacks that are typically used in GLUT code

function ReSizeGLScene(w::Int32,h::Int32)
    if h == 0
        h = 1
    end

    glViewport(0,0,w,h)

    glMatrixMode(GL_PROJECTION)
    glLoadIdentity()

    gluPerspective(45.0,w/h,0.1,100.0)

    glMatrixMode(GL_MODELVIEW)

    return nothing
end

_ReSizeGLScene = cfunction(ReSizeGLScene, Void, (Int32, Int32))

function DrawGLScene()
    global zoom
    global stars
    global tilt
    global twinkle
    global STAR_NUM
    global tex
    global spin

    glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
    glLoadIdentity()

    glBindTexture(GL_TEXTURE_2D,tex[1])

    for loop = 1:STAR_NUM

        glLoadIdentity()

        glTranslate(0.0, 0.0, zoom)

        glRotate(tilt,1.0,0.0,0.0)
        glRotate(stars[loop].angle, 0.0, 1.0, 0.0)

        glTranslate(stars[loop].dist, 0.0, 0.0)

        glRotate(-stars[loop].angle, 0.0, 1.0, 0.0)
        glRotate(-tilt,1.0,0.0,0.0)

        if twinkle
            glColor4ub(stars[STAR_NUM - loop + 1].r,stars[STAR_NUM - loop + 1].g,stars[STAR_NUM - loop + 1].b,255)

            glBegin(GL_QUADS)
                glTexCoord(0.0, 0.0)
                glVertex(-1.0, -1.0, 0.0)
                glTexCoord(1.0, 0.0)
                glVertex(1.0, -1.0, 0.0)
                glTexCoord(1.0, 1.0)
                glVertex(1.0, 1.0, 0.0)
                glTexCoord(0.0, 1.0)
                glVertex(-1.0, 1.0, 0.0)
            glEnd()
        end

        # main star

        glRotate(spin, 0.0, 0.0, 1.0)
        glColor4ub(stars[loop].r, stars[loop].g, stars[loop].b, 255)

        glBegin(GL_QUADS)
            glTexCoord(0.0, 0.0)
            glVertex(-1.0, -1.0, 0.0)
            glTexCoord(1.0, 0.0)
            glVertex(1.0, -1.0, 0.0)
            glTexCoord(1.0, 1.0)
            glVertex(1.0, 1.0, 0.0)
            glTexCoord(0.0, 1.0)
            glVertex(-1.0, 1.0, 0.0)
        glEnd()

        spin              +=0.01
        stars[loop].angle +=loop/STAR_NUM
        stars[loop].dist  -=0.01

        if stars[loop].dist < 0.0
            stars[loop].dist  +=5.0
            stars[loop].r     = rand(1:256)
            stars[loop].g     = rand(1:256)
            stars[loop].b     = rand(1:256)
        end

    end

    glutSwapBuffers()

    return nothing
end

_DrawGLScene = cfunction(DrawGLScene, Void, ())

function keyPressed(the_key::Char,x::Int32,y::Int32)
    global twinkle

    if the_key == int('q')
        glutDestroyWindow(window)
    elseif the_key == int('t')
        println("Twinkle was: $twinkle")
        twinkle = (twinkle ? false : true)
        println("Twinkle is now: $twinkle")
    end

    return nothing # keyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_keyPressed = cfunction(keyPressed, Void, (Char, Int32, Int32))

function specialKeyPressed(the_key::Int32,x::Int32,y::Int32)
    global zoom
    global tilt

    if the_key == GLUT_KEY_PAGE_UP
        zoom -= 0.02
    elseif the_key == GLUT_KEY_PAGE_DOWN
        zoom += 0.02
    elseif the_key == GLUT_KEY_UP
        tilt -= 0.5
    elseif the_key == GLUT_KEY_DOWN
        tilt += 0.5
    end

    return nothing # specialKeyPressed returns "void" in C. this is a workaround for Julia's "automatically return the value of the last expression in a function" behavior.
end

_specialKeyPressed = cfunction(specialKeyPressed, Void, (Int32, Int32, Int32))

# run GLUT routines

glutInit()
glutInitDisplayMode(GLUT_RGBA | GLUT_DOUBLE | GLUT_ALPHA | GLUT_DEPTH)
glutInitWindowSize(width, height)
glutInitWindowPosition(0, 0)

window = glutCreateWindow("NeHe Tut 9")

glutDisplayFunc(_DrawGLScene)
glutFullScreen()

glutIdleFunc(_DrawGLScene)
glutReshapeFunc(_ReSizeGLScene)
glutKeyboardFunc(_keyPressed)
glutSpecialFunc(_specialKeyPressed)

initGL(width, height)

glutMainLoop()


text:
############################################################
# Linked list with Tuple{S, Union{T, Nothing}}
############################################################

mutable struct Orbit4{T, S} <: AbstractOrbit{T,S}
    pt::T # starting point
    last::T
    elts::Dict{T, Tuple{S, Union{T, Nothing}}} # idx -> (content, nextptr)
end

function Orbit4(pt::T, val::S) where {T,S}
    return Orbit4(pt, pt, Dict{T, Tuple{S, Union{T,Nothing}}}(pt=>(val,nothing)))
end

Orbit4(pt::T) where {T} = Orbit4(pt, nothing)

@inline Base.in(x, orb::Orbit4) = haskey(orb.elts, x)

@inline function Base.push!(orb::Orbit4{T,S}, tup::Tuple{T,S}) where {T,S}
    new, val = tup
    orb.elts[orb.last] = (first(orb.elts[orb.last]), new)
    orb.last = new
    push!(orb.elts, new => (val, nothing))
    return orb
end

@inline Base.push!(orb::Orbit4{T,Nothing}, n::T) where T = push!(orb, (n,nothing))

@inline Base.getindex(orb::Orbit4, o) = first(orb.elts[o])

@inline islast(orb::Orbit4, s) = s == orb.last
@inline function next(orb::Orbit4{T}, s) where {T}
    return (last(orb.elts[s]), last(orb.elts[s]))::Tuple{T,T}
end

@inline Base.iterate(orb::Orbit4) = orb.pt, orb.pt
@inline Base.iterate(orb::Orbit4, s) = (islast(orb, s) ? nothing : next(orb, s))
@inline Base.length(orb::Orbit4) = length(orb.elts)

@inline function Base.:(==)(o1::Orbit4, o2::Orbit4)
    o1.pt == o2.pt && o1.last == o2.last || return false
    return o1.elts == o2.elts
end

@inline Base.first(orb::Orbit4) = orb.pt


text:
function assignterminals!(system::LegSystem,R::Vector{Float64},C::Vector{Float64},
    old=Dict("a"=>0),restart="no",assim="no",sample="no")

    # construct and define terminals
    termctr = 0;
    for i in 1:length(system.branches.ID)
        if isempty(system.branches.children[i])
            termctr+=1;
            if assim == "no" && sample == "no"
                system.branches.term[i] = LegHemodynamics.ArterialTerminal();
            end
            # compliance, lower compartments
            if assim == "no" && sample == "no"
                if restart == "no"
                    push!(system.branches.term[i].C,C[termctr]);
                elseif restart == "yes"
                    temp = old[i]["C"];
                    push!(system.branches.term[i].C,temp[1]);
                end
            elseif assim == "yes" || sample == "yes"
                system.branches.term[i].C[1] = C[termctr];
            end
            # resistance, lower compartments
            if assim == "no" && sample == "no"
                if restart == "no"
                    push!(system.branches.term[i].R,system.solverparams.rho*
                        system.branches.c0[i][end]/system.branches.A0[i][end]);
                    push!(system.branches.term[i].R,R[termctr]-system.branches.term[i].R[1]);
                    # println("Second resistance, artery $i: $(system.branches.term[i].R[end])")
                elseif restart == "yes"
                    temp = old[i]["R"];
                    push!(system.branches.term[i].R,temp[1]);
                    push!(system.branches.term[i].R,temp[2]);
                end
            elseif assim == "yes" || sample == "yes"
                system.branches.term[i].R[1] = system.solverparams.rho*
                    system.branches.c0[i][end]/system.branches.A0[i][end];
                system.branches.term[i].R[2] = R[termctr];
            end
        else
            if assim == "no" && sample == "no"
                system.branches.term[i] = LegHemodynamics.ArterialTerminal();
                if !isempty(system.branches.term[i].C)
                    system.branches.term[i].C[1] += NaN;
                else
                    push!(system.branches.term[i].C,NaN);
                end
                if !isempty(system.branches.term[i].R)
                    system.branches.term[i].R[1] += NaN;
                else
                    push!(system.branches.term[i].R,NaN);
                end
            end
        end
    end

end

max_stars_repo_name: dcanuto/leg-hemodynamics

max_issues_repo_name: dcanuto/leg-hemodynamics

max_forks_repo_name: dcanuto/leg-hemodynamics


text:
module Spandex

include("Intlist.jl")
include("SparseArray.jl")
include("DirectedGraph.jl")
include("SegmentTree.jl")
include("PriorityQueue.jl")

include("SparseMatrix.jl")
include("EliminationGraph.jl")
include("PermuteTable.jl")

include("SparseArithmetic.jl")
include("Normalization.jl")
include("CholeskySolver.jl")

end


text:
##---------------------CALLBACKS-------------------##

# Update the figure according to the date range that was picked
# Refresh the figure in the interval defined by INTERVAL. We will only do this
# if we are looking at live data, i.e. if the date range contains the last
# update.
callback!(app,
    Output("log-graph", "figure"),
    [
        Input("date-picker-range", "start_date"),
        Input("date-picker-range", "end_date"),
        Input("refresh-interval", "n_intervals"),
    ],
    State("log-graph", "figure") # we'll pull the data out of it
) do start_date_str, end_date_str, n_intervals, fig
    # parse dates as DateTime
    start_date = DateTime(start_date_str)
    end_date   = DateTime(end_date_str)

    # the `fig` variable comes is a`JSON3.Object{...}`:
    #     fig = {
    #      "data": [
    #                {
    #                      "x": [],
    #                      "y": [],
    #                   "name": "Temperature",
    #                   "line": {
    #                              "color": "OrangeRed"
    #                           }
    #                },
    #                {
    #                       "x": [],
    #                       "y": [],
    #                    "name": "Humidity",
    #                   "yaxis": "y2",
    #                    "line": {
    #                               "color": "SteelBlue"
    #                            }
    #                }
    #              ],
    #    "layout": {
    #                      "yaxis": {
    #                                      "title": {
    #                                                  "text": "Temperature (°C)",
    #                                                  "font": {
    #                                                             "color": "OrangeRed"
    #                                                          }
    #                                               },
    #                                   "tickfont": {
    #                                                  "color": "OrangeRed"
    #                                               },
    #                                      "range": [
    #                                                 -1,
    #                                                 4
    #                                               ],
    #                                  "autorange": true
    #                               },
    #                     "yaxis2": {
    #                                       "title": {
    #                                                   "text": "relative Humidity (%)",
    #                                                   "font": {
    #                                                              "color": "SteelBlue"
    #                                                           }
    #                                                },
    #                                    "tickfont": {
    #                                                   "color": "SteelBlue"
    #                                                },
    #                                  "overlaying": "y",
    #                                        "side": "right",
    #                                       "range": [
    #                                                  -1,
    #                                                  4
    #                                                ],
    #                                   "autorange": true
    #                               },
    #                 "uirevision": "no",
    #                      "xaxis": {
    #                                      "range": [
    #                                                 -1,
    #                                                 6
    #                                               ],
    #                                  "autorange": true
    #                               }
    #              }
    # }
    # This is immutable. We convert it to a mutable Dict with `copy`
    dfig = copy(fig) # this of type `Dict`

    # in the first go, load the data into the dataframe and put that data into
    # the figure
    if n_intervals == 0
        load_logdata!(df, LOGDATA_PATH)
        replace_figure_data!(dfig, df, (start_date, end_date))
    end

    ctx = Dash.callback_context()
    isempty(ctx.triggered) && return dfig
    if ctx.triggered[1].prop_id == "refresh-interval.n_intervals"
        # We have an auto update. So we'll refresh the figure. We'll only do
        # this if the selected end date is after in the future
        end_date < now() && return dfig

        load_logdata!(df, LOGDATA_PATH)
        update_figure!(dfig, df, end_date)
        return dfig
    elseif ctx.triggered[1].prop_id in ["date-picker-range.start_date", "date-picker-range.end_date"]
        replace_figure_data!(dfig, df, (start_date, end_date))
        return dfig
    else
        @warn "Callback was not handled"
        return dfig
    end
end

# Update the temperature stats
callback!(app,
    Output("temperature-stats", "children"),
    [
        Input("date-picker-range", "start_date"),
        Input("date-picker-range", "end_date"),
        Input("refresh-interval", "n_intervals"),
    ]
) do start_date_str, end_date_str, n_intervals
    # Don't do anything for the initial run (n_intervals == 0)
    n_intervals == 0 && return html_div()

    start_date = Dates.DateTime(start_date_str)
    end_date = Dates.DateTime(end_date_str)
    stats = statistics(df, :temperature, (start_date, end_date))
    generate_stats(stats, "°C")
end

# Update the temperature stats
callback!(app,
    Output("humidity-stats", "children"),
    [
        Input("date-picker-range", "start_date"),
        Input("date-picker-range", "end_date"),
        Input("refresh-interval", "n_intervals"),
    ]
) do start_date_str, end_date_str , n_intervals
    # Don't do anything for the initial run (n_intervals == 0)
    n_intervals == 0 || isempty(df[]) && return html_div()

    start_date = Dates.DateTime(start_date_str)
    end_date = Dates.DateTime(end_date_str)
    stats = statistics(df, :humidity, (start_date, end_date))
    generate_stats(stats, "%")
end


text:
# Similar to test for standard A*
using Graphs
using Test

# The "touring Romania" example from Russell and Norvig
g1_wedges = [
    (1, 20, 75.),
    (1, 16, 140.), #2, Arad -- Sibiu
    (1, 17, 118.),
    (20, 13, 71.),
    (13, 16, 151.),
    (16, 6, 99.),
    (16, 15, 80.), #7, Sibiu -- Rimnicu Vilcea
    (6, 2, 211.),
    (17, 10, 111.),
    (10, 11, 70.),
    (11, 4, 75.),
    (4, 3, 120.),
    (3, 15, 146.),
    (3, 14, 138.),
    (15, 14, 97.), #15, Rimnicu Vilcea -- Pitesti
    (14, 2, 101.), #16, Pitesti -- Bucharest
    (2, 7, 90.),
    (2, 18, 85.),
    (18, 8, 98.),
    (8, 5, 86.),
    (18, 19, 142.),
    (19, 9, 92.),
    (9, 12, 87.) ]

g1_heuristics = [
    366, 0, 160, 242, 161, 176, 77, 151, 226, 244, 241, 234, 380, 100, 193,
    253, 329, 80, 199, 374 ]

# Create a dictionary to map indices to neighbours and edge weights
nbrs_wt_dict = Dict{Int64, Dict{Int64, Float64}}()

for (u, v, wt) in g1_wedges
    if haskey(nbrs_wt_dict, u)
        nbrs_wt_dict[u][v] = wt
    else
        nbrs_wt_dict[u] = Dict(v => wt)
    end
end

# Setup edge weight function based on dictionary
edge_wt_fn(u::Int64, v::Int64) = nbrs_wt_dict[u][v]

struct TestVisitorImplicit <: AbstractDijkstraVisitor
    nbrs_wt_dict::Dict{Int64, Dict{Int64, Float64}}
    test_graph::SimpleVListGraph
    goal_vtx::Int64
end

function Graphs.include_vertex!(vis::TestVisitorImplicit, u::Int64, v::Int64,
                                d::Float64, nbrs::Vector{Int64})

    # Terminate if goal vertex reached
    if v == vis.goal_vtx
        return false
    end

    # Otherwise, insert all neighbours from dict
    # AND add to graph if not already there
    for n in collect(keys(vis.nbrs_wt_dict[v]))
        if ~(n in vis.test_graph.vertices)
            add_vertex!(test_graph, n)
        end

        # Neighbors will have INDICES of vertices, not vertices themselves
        push!(nbrs, vertex_index(vis.test_graph, n))
    end

    return nbrs
end

# Create initial graph with only start
start = 1
start_idx = 1
goal = 2
test_graph = SimpleVListGraph([start])

# Create visitor with goal vertex 2 and heuristic
vis = TestVisitorImplicit(nbrs_wt_dict, test_graph, goal)
heur(n) = convert(Float64, g1_heuristics[n])

a_star_sp_states = a_star_implicit_shortest_path!(test_graph, edge_wt_fn, start_idx, vis, heur)

# Now actually extract the shortest path indices
sp_idxs = shortest_path_indices(a_star_sp_states.parent_indices, test_graph, start_idx, vertex_index(test_graph, goal))
sp = [test_graph.vertices[v] for v in sp_idxs]

@test sp == [1, 16, 15, 14, 2]


text:
using ModelingToolkit

### Define a differential equation system

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
de = ODESystem(eqs, t, [x,y,z], [σ,ρ,β])
ode_f = ODEFunction(de)

### Use in DifferentialEquations.jl

using OrdinaryDiffEq
u₀ = ones(3)
tspan = (0.0,100.0)
p = [10.0,28.0,10/3]
prob = ODEProblem(ode_f,u₀,tspan,p)
sol = solve(prob,Tsit5())

using Plots
plot(sol,vars=(1,2,3))

generate_function(de)[1]

generate_function(de)[2]

jac = calculate_jacobian(de)

jac_expr = generate_jacobian(de)

ModelingToolkit.generate_factorized_W(de)[1]

@variables x y z
@parameters σ ρ β

# Define a nonlinear system
eqs = [0 ~ σ*(y-x),
       0 ~ x*(ρ-z)-y,
       0 ~ x*y - β*z]
ns = NonlinearSystem(eqs, [x,y,z], [σ,ρ,β])
nlsys_func = generate_function(ns)

nl_f = @eval eval(nlsys_func[2])
# Make a closure over the parameters for for NLsolve.jl
f2 = (du,u) -> nl_f(du,u,(10.0,26.0,2.33))

using NLsolve
nlsolve(f2,ones(3))

@derivatives D3'''~t
@derivatives D2''~t
@variables u(t), x(t)
eqs = [D3(u) ~ 2(D2(u)) + D(u) + D(x) + 1
       D2(x) ~ D(x) + 2]
de = ODESystem(eqs, t, [u,x], [])
de1 = ode_order_lowering(de)

de1.eqs

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t Dx'~x Dy'~y Dz'~z
eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       D(z) ~ x*y - β*z]
J = [Dx(eqs[1].rhs) Dy(eqs[1].rhs) Dz(eqs[1].rhs)
 Dx(eqs[2].rhs) Dy(eqs[2].rhs) Dz(eqs[2].rhs)
 Dx(eqs[3].rhs) Dy(eqs[3].rhs) Dz(eqs[3].rhs)]

J = expand_derivatives.(J)

using LinearAlgebra
luJ = lu(J,Val(false))

luJ.L

invJ = inv(luJ)

function lorenz(du,u,p,t)
 du[1] = p[1]*(u[2]-u[1])
 du[2] = u[1]*(p[2]-u[3]) - u[2]
 du[3] = u[1]*u[2] - p[3]*u[3]
end

u = [x,y,z]
du = similar(u)
p = [σ,ρ,β]
lorenz(du,u,p,t)
du

J = [Dx(du[1]) Dy(du[1]) Dz(du[1])
     Dx(du[2]) Dy(du[2]) Dz(du[2])
     Dx(du[3]) Dy(du[3]) Dz(du[3])]
J = expand_derivatives.(J)

using SparseArrays
function SparseArrays.SparseMatrixCSC(M::Matrix{T}) where {T<:ModelingToolkit.Expression}
    idxs = findall(!iszero, M)
    I = [i[1] for i in idxs]
    J = [i[2] for i in idxs]
    V = [M[i] for i in idxs]
    return SparseArrays.sparse(I, J, V, size(M)...)
end
sJ = SparseMatrixCSC(J)

@parameters σ(..)
eqs = [D(x) ~ σ(t-1)*(y-x),
       D(y) ~ x*(σ(t^2)-z)-y,
       D(z) ~ x*y - β*z]

@derivatives Dₜ'~t
Dₜ(x*(σ(t^2)-z)-y)
expand_derivatives(Dₜ(x*(σ(t^2)-z)-y))

_f(x) = 2x + x^2
_f(x)

f(x) = 2x + x^2
@register f(x)

f(x)

function ModelingToolkit.derivative(::typeof(f), args::NTuple{1,Any}, ::Val{1})
    2 + 2args[1]
end
expand_derivatives(Dx(f(x)))

using SciMLTutorials
SciMLTutorials.tutorial_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])

max_stars_repo_name: SciML/SciMLTutorialsOutput

max_issues_repo_name: SciML/SciMLTutorialsOutput

max_forks_repo_name: SciML/SciMLTutorialsOutput


text:
import Base.length

abstract type AbstractColoredNode end

struct RedNode <: AbstractColoredNode end; const R = RedNode()
struct BlackNode <: AbstractColoredNode end; const B = BlackNode()
struct Empty end; const E = Empty()
length(e::Empty) = 1

function balance(b::BlackNode, v::Vector, z, d)
    if v[1] == R
        if length(v[2]) == 4 && v[2][1] == R
            return [R, [B, v[2][2], v[2][3], v[2][4]], v[3], [B, v[4], z, d]]
        elseif length(v[4]) == 4 && v[4][1] == R
            return [R, [B, v[2], v[3], v[4][2]], v[4][3], [B, v[4][4], z, d]]
        end
    end
    [b, v, z, d]
end

function balance(b::BlackNode, a, x, v::Vector)
    if v[1] == R
        if length(v[2]) == 4 && v[2][1] == R
            return [R, [B, a, x, v[2][2]], v[2][3], [B, v[2][4], v[3], v[4]]]
        elseif length(v[4]) == 4 && v[4][1] == R
            return [R, [B, a, x, v[2]], v[3], [B, v[4][2], v[4][3], v[4][4]]]
        end
    end
    [b, a, x, v]
end

function balance(b::BlackNode, a::Vector, x, v::Vector)
    if v[1] == R
        if length(v[2]) == 4 && v[2][1] == R
            return [R, [B, a, x, v[2][2]], v[2][3], [B, v[2][4], v[3], v[4]]]
        elseif length(v[4]) == 4 && v[4][1] == R
            return [R, [B, a, x, v[2]], v[3], [B, v[4][2], v[4][3], v[4][4]]]
        end
    end
    [b, a, x, v]
end

balance(node, l, a, r) = [node, l, a, r]

function ins(v::Vector, x::Number)
    if length(v) == 4
        if x < v[3]
            return balance(v[1], ins(v[2], x), v[3], v[4])
        elseif x > v[3]
            return balance(v[1], v[2], v[3], ins(v[4], x))
        end
    end
    v
end

ins(t, a) = [R, E, a, E]

insert(v, a) = (t = ins(v, a); t[1] = B; t)

function testRB()
    t = E
    for i in rand(collect(1:20), 10)
        t = insert(t, i)
    end
    println(replace(string(t), r"lackNode\(\)|edNode\(\)|Any|mpty\(\)" => ""))
end

testRB()

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# process-based-bank-model.jl
# SimLynx.jl Process-Based Discrete-Event Simulation Example

using SimLynx
SimLynx.greet()

using Distributions: Exponential, Uniform
using Random

const N_TELLERS = 2
const MEAN_INTERARRIVAL_TIME = 4.0
const MIN_SERVICE_TIME = 2.0
const MAX_SERVICE_TIME = 10.0

tellers = nothing

"Process to generate n customers arriving into the system."
@process generator(n::Integer) begin
    dist = Exponential(MEAN_INTERARRIVAL_TIME)
    for i = 1:n
        @schedule now customer(i)
        # We don't want to wait after the last customer, which could extend the
        # simulation time past the last customer leaving.
        if i < n
            work(rand(dist))
        end
    end
end

"The ith customer into the system."
@process customer(i::Integer) begin
    dist = Uniform(MIN_SERVICE_TIME, MAX_SERVICE_TIME)
    @with_resource tellers begin
        work(rand(dist))
    end
end

"Run the simulation for n customers."
function run_simulation(n::Integer)
    @assert n > 0
    println("SimLynx.jl Process-Based Discrete-Event Simulation Example")
    println("$N_TELLERS Teller, Single Queue Bank Model")
    println("  Inter-arrival time = Exponential($MEAN_INTERARRIVAL_TIME)")
    println("  Service time = Uniform($MIN_SERVICE_TIME, $MAX_SERVICE_TIME)")
    println("  Number of customers = $n")
    @simulation begin
        global tellers = Resource(N_TELLERS, "tellers")
        @schedule at 0.0 generator(n)
        start_simulation()
        println("--- Simulation results after $(current_time()) minutes ---")
        print_stats(tellers.allocated, title="Teller Allocation Statistics")
        print_stats(tellers.queue_length, title="Teller Queue Length Statistics")
        plot_history(tellers.queue_length, title="Teller Queue Length History")
        print_stats(tellers.wait, title="Teller Queue Wait Time Statistics")
        plot_history(tellers.wait, title="Teller Queue Wait Time History")
    end
end

run_simulation(100)


text:
"""
    TypeE()

Lorem ipsum dolor sit amet.
"""
struct TypeE

end


text:
module Optimize

using LinearAlgebra
using ProximalOperators

abstract type AbstractOptimizer end;
abstract type AbstractSubspaceOptimizer end;
abstract type AbstractScalarizationMethod end;

include("./strridge.jl")
include("./admm.jl")
include("./sr3.jl")

#Nullspace for implicit sindy
include("./adm.jl")

export init, init!, fit!, set_threshold!, get_threshold
export STRRidge, ADMM, SR3
export ADM

include("./pareto.jl")
export ParetoCandidate
export point, parameter, iter, threshold

export WeightedSum, WeightedExponentialSum, GoalProgramming
export weights

export ParetoFront
export assert_dominance, conditional_add!, set_candidate!

end


text:
using Qecsim
using Documenter

DocMeta.setdocmeta!(Qecsim, :DocTestSetup, :(using Qecsim); recursive=true)

makedocs(;
    modules=[Qecsim],
    authors="David K. Tuckett, Terry Farrelly",
    repo="https://github.com/qecsim/Qecsim.jl/blob/{commit}{path}#{line}",
    sitename="Qecsim.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://qecsim.github.io/Qecsim.jl",
        assets=String[],
    ),
    pages=[
        "Overview" => "index.md",
        "API" => [
            "api/core.md",
            "api/models.md",
            "api/index.md",
        ],
    ],
)

deploydocs(;
    repo="github.com/qecsim/Qecsim.jl",
    devbranch = "main",
)