text:
# Start from simple LP
# Solve it
# Copy and solve again
# Chg coeff, solve, change back solve
# del constr and solve
# del var and solve

# a simple LP example
#
#   maximize x + y
#
#   s.t. 2 x + 1 y <= 4
#        1 x + 2 y <= 4
#        x >= 0, y >= 0
#
#   solution: x = 1.3333333, y = 1.3333333, objv = 2.66666666
using Xpress, Compat.Test
using Compat.SparseArrays

@testset "Basics 4" begin

    model = Xpress.Model("lp_04", :maximize)

    add_cvars!(model, [1., 1.], [0., 0.], Inf)

    add_constrs!(model, Cint[1, 3], Cint[1, 2, 1, 2],
        [2., 1., 1., 2.], '<', [4., 4.])

    @test Xpress.get_obj(model) == [1, 1]
    @test Xpress.get_constrmatrix(model) == sparse([2. 1.; 1. 2.])
    @test Xpress.get_rhs(model) == [4, 4]
    @test Xpress.get_lb(model) == [0, 0]
    @test Xpress.get_ub(model) == [1e20, 1e20]
    @test Xpress.get_sense(model) == [:<=,:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [1.33333, 1.33333]; atol = 1e-3)

    @test isapprox(get_objval(model), 2.666666; atol = 1e-3)

    # PART 2:
    # copy and solve

    model2 = copy(model)

    @test Xpress.get_obj(model2) == [1, 1]
    @test Xpress.get_constrmatrix(model2) == sparse([2. 1.; 1. 2.])
    @test Xpress.get_rhs(model2) == [4, 4]
    @test Xpress.get_lb(model2) == [0, 0]
    @test Xpress.get_ub(model2) == [1e20, 1e20]
    @test Xpress.get_sense(model2) == [:<=,:<=]

    optimize(model2)

    ans = get_optiminfo(model2)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model2), [1.33333, 1.33333]; atol = 1e-3)

    @test isapprox(get_objval(model2), 2.666666; atol = 1e-3)

    # PART 3:
    # change coeff and solve

    #   maximize x + y
    #
    #   s.t. 2 x + 2 y <= 4
    #        1 x + 2 y <= 4
    #        x >= 0, y >= 0
    #
    #   solution: x = 0, y = 2, objv = 2

    chg_coeffs!(model, [1], [2],  [2.])

    @test Xpress.get_obj(model) == [1, 1]
    @test Xpress.get_constrmatrix(model) == sparse([2. 2.; 1. 2.])
    @test Xpress.get_rhs(model) == [4, 4]
    @test Xpress.get_lb(model) == [0, 0]
    @test Xpress.get_ub(model) == [1e20, 1e20]
    @test Xpress.get_sense(model) == [:<=,:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [0, 2]; atol = 1e-3)

    @test isapprox(get_objval(model), 2; atol = 1e-3)

    # PART 4:
    # change coeff and solve

    #   maximize x + y
    #
    #   s.t. 1 x + 2 y <= 4
    #        x >= 0, y >= 0
    #
    #   solution: x = 4, y = 0, objv = 4

    del_constrs!(model, [1])

    @test Xpress.get_obj(model) == [1, 1]
    @test Xpress.get_constrmatrix(model) == sparse([1. 2.])
    @test Xpress.get_rhs(model) == [4]
    @test Xpress.get_lb(model) == [0, 0]
    @test Xpress.get_ub(model) == [1e20, 1e20]
    @test Xpress.get_sense(model) == [:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [4, 0]; atol = 1e-3)

    @test isapprox(get_objval(model), 4; atol = 1e-3)

    # PART 5:
    # del var and solve

    #   maximize y
    #
    #   s.t.  2 y <= 4
    #           y >= 0
    #
    #   solution: y = 2, objv = 2

    del_vars!(model, [1])

    @test Xpress.get_obj(model) == [1]
    @test Xpress.get_constrmatrix(model) == sparse([2;]')
    @test Xpress.get_rhs(model) == [4]
    @test Xpress.get_lb(model) == [0]
    @test Xpress.get_ub(model) == [1e20]
    @test Xpress.get_sense(model) == [:<=]

    optimize(model)

    ans = get_optiminfo(model)
    @test ans.status_lp == :optimal

    @test isapprox(get_solution(model), [2]; atol = 1e-3)

    @test isapprox(get_objval(model), 2; atol = 1e-3)

end


text:
__precompile__()

module ClassicNNLS

export ClassicNNLSSolver

"""
Solves argmin ||X*d-Y||^2 subject to d(i) >= 0.

The original version of this code was developed by
Charles L. Lawson and Richard J. Hanson at Jet Propulsion Laboratory
1973 JUN 12, and published in the book
"SOLVING LEAST SQUARES PROBLEMS", Prentice-HalL, 1974.
Revised FEB 1995 to accompany reprinting of the book by SIAM.
"""
function ClassicNNLSSolver(X,Y,MAX_ITER=1000,tol=1e-8)

        N = size(X,2);

        # initialization
        global d = zeros(Float32,N,1);
        global w = transpose(X)*(Y-X*d);
        global R = ones(Bool,N,1);
        global P = zeros(Bool,N,1);
        global mw = maximum(w);

        # Main loop
        for iter in 1:MAX_ITER
            A = findmax(w);
            mw = A[1];
            mIX = A[2];

            # Look for the worst initial solution and include it in the active set
            P[mIX] = true;
            R[mIX] = false;

            # Solve in the active set
            s = zeros(Float32,N,1);
            SEL = P.>0;

            s[SEL] = (X[:,findall(SEL .== true)])\Y;
            # Non-negativity constraint
            neg_subset = s .< 0;
            global d;
            global s;
            global P;
            global R;
            while (sum(neg_subset) > 0)
                neg_subset = (s .< 0);
                alfa = (d[neg_subset]/(d[neg_subset]-s[neg_subset]));
                alfac = findmin(alfa);
                alfac = alfac[1];
                od = d;
                d = d + alfac*(s-d);
                # update the active set
                ToBeUpdated = findall(BitArray(d .== 0.0) .& BitArray(P .== true));

                L2Update = length(ToBeUpdated);
                P[ToBeUpdated] = repeat([false],inner=(1,1),outer=(L2Update,1));
                R[ToBeUpdated] = repeat([true],inner=(1,1),outer=(L2Update,1));

                s = zeros(Float32,N,1);
                SEL = P.>0;
                s[SEL] = (X[:,findall(SEL .== true)])\Y;
                neg_subset = s .< 0;
            end
            d = s;
            w[:] = X'*(Y-X*d);

            if (sum(R) == 0 || mw < tol)
                break
            end
        end

        return d;
end
end #module

ext: jl


text:
#########################
### Forest fire model ###
#########################
using Agents
using Random

mutable struct Tree{T<:Integer} <: AbstractAgent
  id::T
  pos::Tuple{T, T}
  status::Bool  # true is green and false is burning
end

mutable struct Forest{T<:AbstractSpace, Y<:AbstractVector, Z<:AbstractFloat} <: AbstractModel
  space::T
  agents::Y
  scheduler::Function
  f::Z  # probability that a tree will ignite
  d::Z  # forest density
  p::Z  # probability that a tree will grow in an empty space
end

mutable struct MyGrid{T<:Integer, Y<:AbstractVector} <: AbstractSpace
  dimensions::Tuple{T, T}
  space::SimpleGraph
  agent_positions::Y  # an array of arrays for each grid node
end

# we can put the model initiation in a function
function model_initiation(;f, d, p, griddims, seed)
  Random.seed!(seed)
  # initialize the model
  # we start the model without creating the agents first
  agent_positions = [Int64[] for i in 1:gridsize(griddims)]
  mygrid = MyGrid(griddims, grid(griddims, false, true), agent_positions)
  forest = Forest(mygrid, Array{Tree}(undef, 0), random_activation, f, d, p)

  # create and add trees to each node with probability d, which determines the density of the forest
  for node in 1:gridsize(forest.space.dimensions)
    pp = rand()
    if pp <= forest.d
      tree = Tree(node, (1,1), true)
      add_agent!(tree, node, forest)
    end
  end
  return forest
end

function forest_step!(forest)
  shuffled_nodes = Random.shuffle(1:gridsize(forest.space.dimensions))
  for node in shuffled_nodes  # randomly go through the cells and
    if length(forest.space.agent_positions[node]) == 0  # the cell is empty, maybe a tree grows here?
      p = rand()
      if p <= forest.p
        treeid = forest.agents[end].id +1
        tree = Tree(treeid, (1,1), true)
        add_agent!(tree, node, forest)
      end
    else
      treeid = forest.space.agent_positions[node][1]  # id of the tree on this cell
      tree = id_to_agent(treeid, forest)  # the tree on this cell
      if tree.status == false  # if it is has been burning, remove it.
        kill_agent!(tree, forest)
      else
        f = rand()
        if f <= forest.f  # the tree ignites on fire
          tree.status = false
        else  # if any neighbor is on fire, set this tree on fire too
          neighbor_cells = node_neighbors(tree, forest)
          for cell in neighbor_cells
            treeid = get_node_contents(cell, forest)
            if length(treeid) != 0  # the cell is not empty
              treen = id_to_agent(treeid[1], forest)
              if treen.status == false
                tree.status = false
                break
              end
            end
          end
        end
      end
    end
  end
end

forest = model_initiation(f=0.05, d=0.8, p=0.01, griddims=(20, 20), seed=2)
agent_properties = [:status, :pos]
steps_to_collect_data = collect(1:10)

# aggregators = [length, count]
# data = step!(dummystep, forest_step!, forest, 10, agent_properties, aggregators, steps_to_collect_data)
data = step!(dummystep, forest_step!, forest, 10, agent_properties, steps_to_collect_data)

# 9. explore data visually
visualize_data(data)

# or plot trees on a grid
for i in 1:10
  visualize_2D_agent_distribution(data, forest, Symbol("pos_$i"), types=Symbol("status_$i"), savename="step_$i", cc=Dict(true=>"green", false=>"red"))
end

# 10. Running batch
agent_properties = [:status, :pos]
data = batchrunner(dummystep, forest_step!, forest, 10, agent_properties, steps_to_collect_data, 10)
# Create a column with the mean and std of the :status_count columns from differen steps.
columnnames = vcat([:status_count], [Symbol("status_count_$i") for i in 1:9])
using StatsBase
# combine_columns!(data, columnnames, [StatsBase.mean, StatsBase.std])

# optionally write the results to file
write_to_file(df=data, filename="forest_model.csv")


text:
using Gen
using ZigZagBoomerang
using SparseArrays
using LinearAlgebra
using Distributions

"""
    genpdmp(Flow, θ,
        trace, selection::Selection;
        check=false, observations=EmptyChoiceMap(), c=1e-5)

c is a rejection sampler constant. Modifies θ.
Works best with PDMP samplers with high refreshment rate and
Crank-Nicholson parameter ρ close to one.

```julia
Γ = sparse(Diagonal(1.0, n, n)
μ = zeros(n)
θ = randn(n)
Flow = BouncyParticle(Γ, μ, 2.0, ρ = 0.995)

trace, _ = genpdmp(Flow, θ, trace, selection)
```

"""
function genpdmp(Flow, θ,
        trace, selection::Selection;
        check=false, observations=EmptyChoiceMap(), c=1e1, λ0=Flow.λref)

    args = get_args(trace)
    retval_grad = accepts_output_grad(get_gen_fn(trace)) ? zero(get_retval(trace)) : nothing
    argdiffs = map((_) -> NoChange(), args)

    (_, values_trie, gradient_trie) = choice_gradients(trace, selection, retval_grad)
    values = to_array(values_trie, Float64)
    gradient = to_array(gradient_trie, Float64)

    x = values

    ZigZagBoomerang.refresh!(θ, Flow)
    t = 0.0
    ∇ϕx = copy(θ)
    acc = num = 0
    function ∇ϕ!(y, x, values_trie, trace, args, argdiffs, selection, retval_grad)
        values = x

        values_trie[] = from_array(values_trie[], values)
        @assert keys(values_trie[].leaf_nodes)[1] == :slope
        trace[]  = update(trace[], args, argdiffs, values_trie[])[1]
        (_, _, gradient_trie) = choice_gradients(trace[], selection, retval_grad)
        gradient = to_array(gradient_trie, Float64)

        @. y = -gradient
        y
    end

    Ξ = ZigZagBoomerang.Trace(t, x, θ, Flow)
    τref = T = ZigZagBoomerang.waiting_time_ref(Flow)
    a, b = ZigZagBoomerang.ab(x, θ, c, Flow)
    t′ = t + poisson_time(a, b, rand())
    while t < T
        t, x, θ, (acc, num), c, a, b, t′, τref = ZigZagBoomerang.pdmp_inner!(Ξ, ∇ϕ!, ∇ϕx, t, x, θ, c, a, b, t′, τref, (acc, num), Flow,
            Ref(values_trie), Ref(trace), args, argdiffs, selection, retval_grad; adapt=false)
    end

    values = x
    trace = trace
    values_trie = from_array(values_trie, values)
    (trace, _, _) = update(trace, args, argdiffs, values_trie)

    check && check_observations(get_choices(trace), observations)

    return trace, true
end


text:
export moment, exp, log, series, dual, sparse_pol, sparse_decompose;

#----------------------------------------------------------------------
function Base.exp(A::Matrix{U}, pt::Vector{T}) where {U,T}
    N = size(A)
    R = fill(zero(eltype(pt)), N[1], N[2])
    for i  in 1:N[1]
        for j in 1:N[2]
            R[i,j] = pt[j]^A[i,j]
        end
    end
    R
end

#----------------------------------------------------------------------
function Base.log(A::Array{U}, pt::Vector{T}) where {U,T}
    N = size(A)
    R = fill(zero(eltype(pt)), N[1], N[2])
    for i  in 1:N[1]
        for j in 1:N[2]
            R[i,j] = real(log(A[i,j])/log(pt[i]))
        end
    end
    R
end

#------------------------------------------------------------------------
"""
```
dual(p::Polynomial) -> Series{C,M}
```
Compute the series associated to the polynomial p, replacing
the variables xi by their dual variables dxi. C is the type of coefficients
of the polynomial p and M its type of monomials.
"""
function dual(p::Polynomial{true,C}) where C
    s = Series{C, DynamicPolynomials.Monomial{true}}()
    for t in p
               s[t.x] = t.α
    end
    return s
end

function dual(t::Term{true,T}) where T
    s = Series{T, Monomial{true}}()
    s[t.x] = t.α
    return s
end

#----------------------------------------------------------------------
"""
```
moment(w::Vector{C}, P::Matrix{C}) -> Vector{Int64} -> C
```
Compute the moment function ``α -> ∑_{i} ω_{i} P_{i}^α``
associated to the sequence P of r points of dimension n, which is a matrix
of size r*n and the weights w.
"""
moment(w, P) = function(α)
  res = 0
  for j in 1:length(w)
      m = 1
      for i in 1:length(α)
        m *= P[i,j]^α[i]
      end
      res+=m*w[j];
  end
  res
end

#------------------------------------------------------------------------
"""
```
series(f::Function, L::Vector{M}) -> Series{C,M}
```
Compute the generating series ``\\sum_{x^{α} \\in L} f(α) z^α``
for a function  ``f: \\mathbb{N}^n \\rightarrow C`` and a sequence L of monomials.
"""
function series(f::Function, L::Vector)
   res = series(f(L[1].z), L[1])
   for m in L
         res[m] = f(m.z)
   end
   res
end

#----------------------------------------------------------------------
"""
```
moment(p::Polynomial, zeta::Vector{C}) -> Vector{Int64} -> C
```
Compute the moment function ``α \\rightarrow p(ζ^α)``.
"""
moment(p::Polynomial, zeta::Vector) =  function(V::Vector{Int})
    U = [zeta[i]^V[i] for i in 1:length(V)]
    p(U)
end

#----------------------------------------------------------------------
"""
```
series(w:: Vector{C}, P::Matrix{C}, L::Vector{M}) -> Series{C,M}
```
Compute the series of the moment sequence ``∑_{i} ω_{i} P_{i}^α`` for ``α \\in L``.
"""
function series(w:: Vector{C}, P, L::Vector{M}) where {C,M}
   series(moment(w,P), L)
end

#----------------------------------------------------------------------
"""
```
series(w:: Vector{C}, P::AbstractMatrix, X, d::Int64) -> Series{C,M}
```
Compute the series of the moment sequence ``∑_i ω_{i} P_{i}^α`` for ``|α| \\leq d``.
"""
function series(w::Vector{C}, P::AbstractMatrix, X, d::Int64) where C
    h = moment(w,P)
    L = monomials(X,seq(0:d))
   series(h,L)
end

#----------------------------------------------------------------------
"""
```
series(p::Polynomial, zeta, X, d::Int64) -> Series
```
Compute the series of moments ``p(ζ^α)`` for ``|α| \\leq d``.
"""
function series(p::Polynomial, zeta, X, d::Int64)
    h = moment(p,zeta)
    L = monomials(X,seq(0:d))
    series(h,L)
end

#----------------------------------------------------------------------
"""
```
sparse_pol(w, E, X) -> Polynomial{true,C}
```
Compute the polynomial ``∑ ωᵢ X^E[i,:]`` with coefficients ``ωᵢ`` and monomial exponents ``E[i,:]``.
"""
function sparse_pol(w::Vector, E::Matrix{Int64}, X)
    sum(w[j]*prod(X[i]^E[j,i] for i in 1:size(E,2)) for j in 1:length(w))
end

#------------------------------------------------------------------------
function sparse_decompose(f, zeta, X, d:: Int64)
    sigma = series(moment(f,zeta), monomials(X,seq(0:d)))
    w, Xi = decompose(sigma)
    w, log(Xi,zeta)
end


text:
export rnnConnectGraph

function rnnConnectGraph(L,T)
  connectGraph=zeros(Int64,L*T,L*T)
  for i=1:L*T
    if i%L==1
      connectGraph[i,i+1]=1
    elseif i%L!=0
      connectGraph[i,i+1]=1
      if i+L<L*T
        connectGraph[i,i+L]=1
      end
    else connectGraph[i,i-L+1]=1
    end
  end
  connectGraph
end


text:
using ComoniconOptions
using Documenter

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

makedocs(;
    modules=[ComoniconOptions],
    authors="Roger-Luo <rogerluo.rl18@gmail.com> and contributors",
    repo="https://github.com/Roger-luo/ComoniconOptions.jl/blob/{commit}{path}#{line}",
    sitename="ComoniconOptions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://Roger-luo.github.io/ComoniconOptions.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/Roger-luo/ComoniconOptions.jl",
)


text:
@testset "Solving ED Models" begin
    #5 Bus Test
    ED = PSI.EconomicDispatch(sys5b_uc, PSI.CopperPlatePowerModel; optimizer = GLPK_optimizer);
    res_5 = solve_op_model!(ED)
    @test isapprox(res_5.total_cost[:ED], 2400, atol = 1000)
    #14 Bus Test
    ED = PSI.EconomicDispatch(sys14, PSI.CopperPlatePowerModel; optimizer = OSQP_optimizer);
    res_14 = solve_op_model!(ED)
    @test isapprox(res_14.total_cost[:ED], 1000, atol = 100)
    #RTS Test
    ED = PSI.EconomicDispatch(sys_rts, PSI.CopperPlatePowerModel; optimizer = GLPK_optimizer);
    res_14 = solve_op_model!(ED)
end

@testset "Solving ED with PTDF Models" begin
    # 5 - Bus Test
    PTDF5,  = PowerSystems.buildptdf(branches5, nodes5)
    ED = PSI.EconomicDispatch(sys5b_uc, PSI.StandardPTDFForm; PTDF = PTDF5, optimizer = GLPK_optimizer);
    res_5 = solve_op_model!(ED)
    @test isapprox(res_5.total_cost[:ED], 3400, atol = 1000)

    PTDF5,  = PowerSystems.buildptdf(branches5, nodes5)
    ED = PSI.EconomicDispatch(sys5b_uc, PSI.StandardPTDFForm; PTDF = PTDF5, optimizer = GLPK_optimizer, parameters = false);
    res_5 = solve_op_model!(ED)
    @test isapprox(res_5.total_cost[:ED], 3400, atol = 1000)

    # 14 Bus Test
    PTDF14,  = PSY.buildptdf(sys14.branches, sys14.buses)
    ED = PSI.EconomicDispatch(sys14, PSI.StandardPTDFForm; PTDF = PTDF14, optimizer = OSQP_optimizer);
    res_14 = solve_op_model!(ED)
    @test isapprox(res_14.total_cost[:ED], 1000, atol = 100)

    PTDF14,  = PSY.buildptdf(sys14.branches, sys14.buses)
    ED = PSI.EconomicDispatch(sys14, PSI.StandardPTDFForm; PTDF = PTDF14, optimizer = OSQP_optimizer, parameters = false);
    res_14 = solve_op_model!(ED)
    @test isapprox(res_14.total_cost[:ED], 1000, atol = 100)
end

@testset "testing AngleDC-OPF 5-bus" begin
    # 5 - Bus Test
    ED = PSI.EconomicDispatch(sys5b_uc,PM.DCPlosslessForm; optimizer = GLPK_optimizer);
    res_5= solve_op_model!(ED)
    @test isapprox(res_5.total_cost[:ED], 3400, atol = 1000)

    ED = PSI.EconomicDispatch(sys5b_uc,PM.DCPlosslessForm; optimizer = GLPK_optimizer, parameters = false);
    res_5= solve_op_model!(ED)
    @test isapprox(res_5.total_cost[:ED], 3400, atol = 1000)

    # 14 - Bus Test
    ED = PSI.EconomicDispatch(sys14,PM.DCPlosslessForm; optimizer = ipopt_optimizer);
    res_14 = solve_op_model!(ED)
    @test isapprox(res_14.total_cost[:ED], 1000, atol = 100)

    ED = PSI.EconomicDispatch(sys14,PM.DCPlosslessForm; optimizer = ipopt_optimizer, parameters = false);
    res_14 = solve_op_model!(ED)
    @test isapprox(res_14.total_cost[:ED], 1000, atol = 100)
end

@testset "Test the equivelance of 5bus PM and PSI models" begin
     file = joinpath(dirname(dirname(pathof(PowerModels))),"test/data/matpower/case5.m")
     #ps5 = PowerSystems.parsestandardfiles(file) # this fails because PSY changes to mixed_units

     # adjusted standard file parsing without mixed units
     ps5 = PowerSystems.parse_file(file)
     #make_mixed_units(ps5)
     ps5 = PowerSystems.pm2ps_dict(ps5)
     # a single period simulation should be equivelant to a PM simulation
     for (k,l) in ps5["load"]
         l["scalingfactor"] = l["scalingfactor"][1]
     end
     bus5, gen5, stor5, branch5,load5,lz5,shunts5,service5 = ps_dict2ps_struct(ps5)
     sys5 = PSY.System(bus5,gen5,load5,branch5,stor5,100.0)
     ED5 = PSI.EconomicDispatch(sys5, PSI.PM.DCPlosslessForm; optimizer = ipopt_optimizer)
     res_5 = solve_op_model!(ED5)

     pm5 = PowerModels.parse_file(file)
     PM5 = build_generic_model(pm5,DCPPowerModel,PowerModels.post_opf)
     res_PM5 = solve_generic_model(PM5,ipopt_optimizer)

    @test isapprox(res_5.total_cost[:ED], res_PM5["objective"],atol = 1)
 end

  @testset "Test the equivelance of 14bus PM and PSI models" begin
     file = joinpath(dirname(dirname(pathof(PowerModels))),"test/data/matpower/case14.m")
     #ps14 = PowerSystems.parsestandardfiles(file) # this fails because PSY changes to mixed_units

      # adjusted stndard file parsing without mixed units
     ps14 = PowerSystems.parse_file(file)
     #make_mixed_units(ps14)
     ps14 = PowerSystems.pm2ps_dict(ps14)

      # a single period simulation should be equivelant to a PM simulation
     for (k,l) in ps14["load"]
         l["scalingfactor"] = l["scalingfactor"][1]
     end
     bus14, gen14, stor14, branch14,load14,lz14,shunts14,service14 = ps_dict2ps_struct(ps14)
     sys14 = PSY.System(bus14,gen14,load14,branch14,stor14,100.0)
     ED14 = PSI.EconomicDispatch(sys14, PSI.PM.DCPlosslessForm; optimizer = ipopt_optimizer)
     res_14 = solve_op_model!(ED14)

      pm14 = PowerModels.parse_file(file)
     PM14 = build_generic_model(pm14,DCPPowerModel,PowerModels.post_opf)
     res_PM14 = solve_generic_model(PM14,ipopt_optimizer)

      @test isapprox(res_14.total_cost[:ED], res_PM14["objective"], atol = 1)
 end

#=
@test try
    @info "testing ACP-OPF 5-bus"
    Net = PM.StandardACPForm
    m = Model(ipopt_optimizer);
    netinjection = PSI.instantiate_network(Net, sys5);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalDispatch, Net, sys5);
    PSI.construct_network!(m, [(device=Line, formulation=PSI.PiLine)], netinjection, Net, sys5)
    JuMP.@objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
true finally end

@test try
    @info "testing ACP- QCWForm 5-bus"
    Net = PM.QCWRForm
    m = Model(ipopt_optimizer);
    netinjection = PSI.instantiate_network(Net, sys5);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalDispatch, Net, sys5);
    PSI.construct_network!(m, [(device=Line, formulation=PSI.PiLine)], netinjection, Net, sys5)
    JuMP.@objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 3400, atol = 1000)
true finally end

@test try
    @info "testing AngleDC-OPF 14-bus"
    Net = PM.DCPlosslessForm
    m = Model(ipopt_optimizer);
    netinjection = PSI.instantiate_network(Net, sys14);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalDispatch, Net, sys14);
    PSI.construct_network!(m, [(device=Line, formulation=PSI.PiLine)], netinjection, Net, sys14)
    JuMP.@objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
finally end

@test try
    @info "testing ACP-OPF 14-bus"
    Net = PM.StandardACPForm
    m = Model(ipopt_optimizer);
    netinjection = PSI.instantiate_network(Net, sys14);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalDispatch, Net, sys14);
    PSI.construct_network!(m, [(device=Line, formulation=PSI.PiLine)], netinjection, Net, sys14)
    JuMP.@objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
true finally end

@test try
    @info "testing ACP-QCWForm 14-bus"
    Net = PM.QCWRForm
    m = Model(ipopt_optimizer);
    netinjection = PSI.instantiate_network(Net, sys14);
    PSI.construct_device!(m, netinjection, ThermalGen, PSI.ThermalDispatch, Net, sys14);
    PSI.construct_network!(m, [(device=Line, formulation=PSI.PiLine)], netinjection, Net, sys14)
    JuMP.@objective(m, Min, m.obj_dict[:objective_function])
    JuMP.optimize!(m)
    isapprox(JuMP.objective_value(m), 1200, atol = 1000)
true finally end
=#


text:
import CartesianGrids: dot
import ConstrainedSystems: init

_norm_sq(u) = dot(u,u)
_norm_sq(u::ConstrainedSystems.ArrayPartition) = sum(_norm_sq,u.x)
state_norm(u,t) = sqrt(_norm_sq(u))

# ensure that time marching makes use of
#init(prob;alg=LiskaIFHERK(),kwargs...) = init(prob,alg;internal_norm=state_norm,kwargs...)

function init(u0,tspan,sys::HeatConduction;alg=ConstrainedSystems.LiskaIFHERK(),kwargs...)
    prob = ODEProblem(sys.f,u0,tspan,sys)
    return init(prob, alg,dt=timestep(sys),internal_norm=state_norm,kwargs...)
end


text:
import Pkg; Pkg.activate(@__DIR__)

using ClusterTrees
using StaticArrays
using CompScienceMeshes
using LinearAlgebra

const P = SVector{3,Float64}

mbs = 1.0 / 2^(5-1)
points = [
    point(1,1,0)-point(mbs, mbs, 0) / 2,
    point(-1,-1,0)+point(mbs, mbs, 0) / 2]

ct, sz = point(0,0,0), 1.0
tree = ClusterTrees.LevelledTrees.LevelledTree(ct, sz, Int[])

for (i,pt) in pairs(points)
    dest = (smallest_box_size=mbs, target_point=pt)
    state = ClusterTrees.LevelledTrees.rootstate(tree, dest)
    ClusterTrees.update!(tree, state, i, dest) do tree, node, i
        push!(data(tree, node).values, i)
    end
end

ClusterTrees.print_tree(tree, root(tree))
block_tree = ClusterTrees.BlockTrees.BlockTree(tree,tree)

function updatestate(block_tree, chd, par_state)
    d = ClusterTrees.data(block_tree, chd)
    test_sector  = d[1].sector
    trial_sector = d[2].sector
    cs1 = ClusterTrees.LevelledTrees.center_size(test_sector, par_state[1][1], par_state[1][2])
    cs2 = ClusterTrees.LevelledTrees.center_size(trial_sector, par_state[2][1], par_state[2][2])
    return (cs1,cs2)
end

function listnearfarinteractions(block_tree, block, state, adm, nears, fars, level)
    adm(block_tree, block, state) && (push!(fars[level], block); return)
    !ClusterTrees.haschildren(block_tree, block) && (push!(nears,block); return)
    for chd ∈ children(block_tree, block)
        chd_state = updatestate(block_tree, chd, state)
        listnearfarinteractions(block_tree, chd, chd_state, adm, nears, fars, level+1)
    end
end

function isfar(block_tree, block, state)
    η = 1.1
    test_state, trial_state = state
    test_center, test_size = test_state
    trial_center, trial_size = trial_state
    center_dist = norm(test_center-trial_center)
    @show center_dist, test_size, trial_size
    if (center_dist - sqrt(3)*(test_size+trial_size)) / (test_size+trial_size) > η
        return true
    else
        return false
    end
end

nears = []
num_levels = length(tree.levels)
fars = [[] for l in 1:num_levels]
root_state = ((tree.center,tree.halfsize),(tree.center,tree.halfsize))
root_level = 1
listnearfarinteractions(block_tree, root(block_tree),
    root_state, isfar,nears, fars, root_level)


text:
using Statistics
include("asvar.jl")

function switchingtime(a::Float64,b::Float64,u::Float64=rand())
# generate switching time for rate of the form max(0, a + b s) + c
# under the assumptions that b > 0, c > 0
  if (b > 0)
    if (a < 0)
      return -a/b + switchingtime(0.0, b, u);
    else # a >= 0
      return -a/b + sqrt(a^2/b^2 - 2 * log(u)/b);
    end
  elseif (b == 0) # degenerate case
    if (a < 0)
      return Inf;
    else # a >= 0
      return -log(u)/a;
    end
  else # b <= 0
    if (a <= 0)
      return Inf;
    else # a > 0
      y = -log(u); t1=-a/b;
      if (y >= a * t1 + b *t1^2/2)
        return Inf;
      else
        return -a/b - sqrt(a^2/b^2 + 2 * y /b);
      end
    end
  end
end

function switchingtime_quadratic(a::Float64, b::Float64, c::Float64,u::Float64=rand())
  dummy::Float64 = (-b^3 + 6*a*b*c + 2*(-6*c^2*log(u) + c * sqrt(-3*a^2*b^2 + 16*a^3*c + 6*b^3*log(u) - 36*a*b*c*log(u) + 36*c^2*(log(u))^2)))^(1/3)
  return 1/(2c)*(-b + (b^2-4*a*c)/dummy + dummy)
end

function EstimateMomentAlongTrajectory(t_skeleton::Vector{Float64}, x_skeleton::Vector{Vector{Float64}}, v_skeleton::Vector{Vector{Float64}}; p = 1, ellipticQ::Bool = false, x_ref::Vector{Float64} = Vector{Float64}(undef,0))
  # assume every component follows dynamics x_t = x_ref + (x_0-x_ref)*cos(t) + v_0*sin(t)

  dim = length(x_skeleton[1])

  if ellipticQ && length(x_ref) == 0
    println("WARNING : EstimateMomentAlongTrajectory : No valid reference point provided.")
    x_ref = zeros(dim)
  end
  T = (t_skeleton[end]-t_skeleton[1])
  estimate = 0.0
  for i=2:length(t_skeleton)
    v = v_skeleton[i-1]
    x = x_skeleton[i-1]
    t0 = t_skeleton[i-1]
    t1 = t_skeleton[i]
    Δt = t1 - t0
    if ellipticQ
      if p==1
        estimate = estimate .+ 1/T * (x_ref * Δt + (x - x_ref) * sin(Δt) + v * (1-cos(Δt)))
      elseif p==2
        y = x - x_ref
        estimate = estimate .+ 1/(2*T) * (v .*(4 * x_ref + y) + Δt * (v.^2 + 2 * x_ref.^2 + y.^2) - (4 * x_ref + y * cos(Δt) + v * sin(Δt)).*(v*cos(Δt) - y *sin(Δt)))
      else
        error("Higher order moments (p >2) along Elliptic trajectories currently not implemented.")
      end
    else # assume linear trajectories
      estimate =  estimate .+ 1 ./((p+1)*v_skeleton[i-1] * T) .* ( x_skeleton[i].^(p+1) - x_skeleton[i-1].^(p+1))
    end
  end
  return estimate

end

function EllipticDynamics(t::Real, y0::Vector{Float64}, w0::Vector{Float64})
    # simulate dy/dt = w, dw/dt = -y

    y_new = y0 * cos(t) + w0 * sin(t);
    w_new = -y0 * sin(t) + w0 * cos(t);

    return (y_new, w_new);

end

function SkeletonIntoBatches(t_skeleton::Vector{Float64}, x_skeleton::Vector{Vector{Float64}}, v_skeleton::Vector{Vector{Float64}}, n_batches::Int, ellipticQ::Bool = false, x_ref::Vector{Float64} = Vector{Float64}(undef,0))

  t_skeletons::Vector{Vector{Float64}} = Vector{Vector{Float64}}(undef,0)
  x_skeletons::Vector{Vector{Vector{Float64}}} = Vector{Vector{Vector{Float64}}}(undef,0)
  v_skeletons = similar(x_skeletons)

  if n_batches < 1
    error("n_batches < 1")
  end

  if (ellipticQ && length(x_ref) == 0)
    x_ref = zeros(length(x_skeleton[1]))
  end
  T_max = t_skeleton[end]
  j_max = length(t_skeleton)
  T_batch = T_max / n_batches
  t_current = t_skeleton[1]
  x_current = x_skeleton[1]
  v_current = v_skeleton[1]
  j = 2
  for i = 1:n_batches
    t_skeleton_batch = Vector{Float64}(undef,0)
    x_skeleton_batch = Vector{Vector{Float64}}(undef,0)
    v_skeleton_batch = similar(x_skeleton_batch)
    # add initial position to skeleton
    push!(t_skeleton_batch, t_current)
    push!(x_skeleton_batch, x_current)
    push!(v_skeleton_batch, v_current)

    while (t_skeleton[j] < i*T_batch && j < j_max)
      t_current = t_skeleton[j]
      x_current = x_skeleton[j]
      v_current = v_skeleton[j]
      push!(t_skeleton_batch, t_current)
      push!(x_skeleton_batch, x_current)
      push!(v_skeleton_batch, v_current)
      j = j + 1
    end
    # j = min (j : t_skeleton[j] > i * T_batch)

    # add final position of this batch
    push!(t_skeleton_batch, i*T_batch)
    Δt = (i*T_batch - t_current)
    if  (ellipticQ)
      (y,v_current) = EllipticDynamics(Δt,x_current-x_ref, v_current)
      x_current = x_ref + y
    else
      x_current = x_current + Δt * v_current
    end
    t_current = i*T_batch
    push!(x_skeleton_batch, x_current)
    push!(v_skeleton_batch, v_current)

    # add skeleton to list of skeletons
    push!(t_skeletons, t_skeleton_batch)
    push!(x_skeletons, x_skeleton_batch)
    push!(v_skeletons, v_skeleton_batch)
  end
  return (t_skeletons, x_skeletons, v_skeletons)

end

function BatchMeansPDMP(t_skeleton::Vector{Float64}, x_skeleton::Vector{Vector{Float64}}, v_skeleton::Vector{Vector{Float64}}; n_batches::Int = 100, ellipticQ::Bool = false, x_ref::Vector{Float64} = Vector{Float64}(undef,0))

  dimension = length(x_skeleton[1])
  if (n_batches <= 0)
    error("A valid number of batches is required.")
  end
  T_batch = (t_skeleton[end]-t_skeleton[1])/n_batches
  (t_skeletons, x_skeletons, v_skeletons) = SkeletonIntoBatches(t_skeleton, x_skeleton, v_skeleton, n_batches, ellipticQ, x_ref)
  Z = Array{Float64,2}(undef,dimension,n_batches) # batch means
  for i=1:n_batches
    Z[:,i] = EstimateMomentAlongTrajectory(t_skeletons[i],x_skeletons[i],v_skeletons[i],ellipticQ=ellipticQ,x_ref=x_ref,p=1)
  end
  vec(var(Z, dims=2)*T_batch)
end

function ess_pdmp_components(t_skeleton::Vector{Float64}, x_skeleton::Vector{Vector{Float64}}, v_skeleton::Vector{Vector{Float64}}; n_batches::Int = 50, ellipticQ::Bool = false, x_ref::Vector{Float64} = Vector{Float64}(undef,0))

  asymptotic_variances = BatchMeansPDMP(t_skeleton, x_skeleton, v_skeleton, n_batches = n_batches, ellipticQ = ellipticQ, x_ref = x_ref)
  if any(asymptotic_variances .< 0)
    println("WARNING : ess_pdmp : at least one component of the asymptotic variance is negative.")
    println("asymptotic_variances = ", asymptotic_variances)
  end
  variances = EstimateMomentAlongTrajectory(t_skeleton,x_skeleton,v_skeleton, ellipticQ = ellipticQ, x_ref = x_ref, p=2) .- (EstimateMomentAlongTrajectory(t_skeleton,x_skeleton,v_skeleton, ellipticQ = ellipticQ, x_ref = x_ref).^2)
  if any(variances .< 0)
    println("WARNING : ess_pdmp : at least one component of the variance is negative.")
    println("variances = ", variances)
  end
  ess = variances./asymptotic_variances*(t_skeleton[end]-t_skeleton[1])
end

function ess_pdmp_squared_radius(t_skeleton::Vector{Float64}, x_skeleton::Vector{Vector{Float64}}, v_skeleton::Vector{Vector{Float64}}; n_batches::Int = 50, ellipticQ::Bool = false, x_ref::Vector{Float64} = Vector{Float64}(undef,0), true_mean::Float64 = NaN, true_variance::Float64 = NaN)

  n_samples = n_batches * 1000
  samples = ExtractSamples(t_skeleton, x_skeleton, v_skeleton, n_samples, ellipticQ=ellipticQ, x_ref = x_ref);
  X = sum(samples.^2, dims = 1)
  if isnan(true_mean)
    mean = Vector{Float64}(undef,0)
  else
    mean = [true_mean]
  end
  if isnan(true_variance)
    variance = Vector{Float64}(undef,0)
  else
    variance = [true_variance]
  end
  ESS(X,n_batches = n_batches, true_mean = mean, true_variance = variance)[1]
end

function plot_pdmp(t_skeleton, x_skeleton, v_skeleton; ellipticQ::Bool = false, x_ref::Vector{Float64} = Vector{Float64}(undef,0), name::String = "", n_samples::Int = 0)

  if (n_samples == 0)
    n_samples = Int(round(t_skeleton[end] - t_skeleton[1]))
  end
  samples = ExtractSamples(t_skeleton, x_skeleton, v_skeleton, n_samples, ellipticQ=ellipticQ, x_ref = x_ref);
  plot_ref = plot(samples[1,:], samples[2,:],legend = false, ratio=:equal);
  savefig(plot_ref, string(name,"-trace.pdf"))
  # plot_ref = bar(autocor(samples'))
  # if (length(name) > 0)
  #     savefig(plot_ref, string(name,"-acf.pdf"))
  # end
end

function ExtractSamples(t_skeleton::Vector{Float64}, x_skeleton::Vector{Vector{Float64}}, v_skeleton::Vector{Vector{Float64}}, n_samples::Int; x_ref::Vector{Float64} = Vector{Float64}(undef,0), ellipticQ::Bool=false)

    T = t_skeleton[end]
    Δt = T/(n_samples-1)
    dim = size(x_skeleton[1])[1]
    samples = Matrix{Float64}(undef, dim, n_samples)
    t = 0.0
    if (ellipticQ && length(x_ref) == 0)
        println("WARNING : ExtractSamples : No valid reference point provided.")
        x_ref = zeros(dim)
    end
    x = x_skeleton[1]
    v = v_skeleton[1]
    samples[:,1] = x
    skeleton_index = 2
    counter = 1
    for i in range(2, stop=n_samples)
        t_max = (i-1) * Δt
        while (counter + 1 < length(t_skeleton) && t_skeleton[counter + 1] < t_max)
            counter = counter + 1
            x = x_skeleton[counter]
            v = v_skeleton[counter]
            t = t_skeleton[counter]
        end
        if (ellipticQ)
            (y,v) = EllipticDynamics(t_max - t, x-x_ref, v)
            x = y + x_ref
        else
            x = x + v * (t_max - t)
        end
        t = t_max
        samples[:,i] = x
    end

    return samples
end

max_stars_repo_name: SebaGraz/ZZDiffusionBridge

max_issues_repo_name: SebaGraz/ZZDiffusionBridge

max_forks_repo_name: SebaGraz/ZZDiffusionBridge


text:
function FractionFormat()
    return FractionFormat(())
end

function FractionFormat(arg0::NumberFormat)
    return FractionFormat((NumberFormat,), arg0)
end

function FractionFormat(arg0::NumberFormat, arg1::NumberFormat)
    return FractionFormat((NumberFormat, NumberFormat), arg0, arg1)
end

function format(obj::FractionFormat, arg0::Fraction, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (Fraction, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format(obj::FractionFormat, arg0::Object, arg1::StringBuffer, arg2::FieldPosition)
    return jcall(obj, "format", StringBuffer, (Object, StringBuffer, FieldPosition), arg0, arg1, arg2)
end

function format_fraction(::Type{FractionFormat}, arg0::Fraction)
    return jcall(FractionFormat, "formatFraction", JString, (Fraction,), arg0)
end

function get_available_locales(::Type{FractionFormat})
    return jcall(FractionFormat, "getAvailableLocales", Vector{Locale}, ())
end

function get_improper_instance(::Type{FractionFormat})
    return jcall(FractionFormat, "getImproperInstance", FractionFormat, ())
end

function get_improper_instance(::Type{FractionFormat}, arg0::Locale)
    return jcall(FractionFormat, "getImproperInstance", FractionFormat, (Locale,), arg0)
end

function get_proper_instance(::Type{FractionFormat})
    return jcall(FractionFormat, "getProperInstance", FractionFormat, ())
end

function get_proper_instance(::Type{FractionFormat}, arg0::Locale)
    return jcall(FractionFormat, "getProperInstance", FractionFormat, (Locale,), arg0)
end

function parse(obj::FractionFormat, arg0::JString)
    return jcall(obj, "parse", Fraction, (JString,), arg0)
end

function parse(obj::FractionFormat, arg0::JString, arg1::ParsePosition)
    return jcall(obj, "parse", Fraction, (JString, ParsePosition), arg0, arg1)
end


text:
using NPZ
using LinearAlgebra
using Combinatorics
using ProgressBars
using Printf: @sprintf
using OnlineStats

G1 = npzread("Kernels/G1.npy")
G2 = npzread("Kernels/G2.npy")
G3 = npzread("Kernels/G3.npy")

function d_oed_exhaust(X1::Matrix{T}, X2::Matrix{T}, X3::Matrix{T}, B, ϵ, λ) where T <: Number
    X = [X3 zeros(T, size(X3));
         zeros(T, size(X3)) X3;
         X1 X2]
    # Xn = maximum(sqrt.(sum(X.^2, dims=2)))
    Xn = sqrt.(sum(X.^2, dims=2)) #normalize all rows to be 1?
    X ./= Xn

    n = size(X3, 1)
    xdas = ones(T, size(X1,1))

    function Ml(μ, X, ϵ)
        M = ϵ*I + X'*(Diagonal(μ)*X)
        Mc = cholesky(Hermitian(M))
        l = μ.^2 .*diag(X*(Mc\X'))
        return (M, Mc, l)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, λ)
        M, Mc, l = Ml(μ)
        return -logdet(M) + λ*maximum(l)
    end

    f(μ) = f(μ, λ)

    wrap_f(x) = f(vcat(x, x, xdas))

    combs = combinations(1:n,B)

    f_vals = zeros(length(combs))
    x = zeros(n, length(combs))

    iter = ProgressBar(enumerate(combs))
    stat = Extrema()
    for (i, comb) in iter
        x[comb,i] .= 1
        f_vals[i] = wrap_f(x[:,i])
        fit!(stat, f_vals[i])
        set_description(iter, "$(values(stat))")
    end

    return (f_vals, x)
end

f_vals, x = d_oed_exhaust(G1, G2, G3, 4, 1e-1, 10.0)

npzwrite("res2.npy", x[:,argmin(f_vals)])
# ;

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
using RayTracingWeekend

using BenchmarkTools
using Images
using LinearAlgebra
using InteractiveUtils
using StaticArrays
using Test

"Modified from BenchmarkTools @ballocated. Returns number of heap allocations, instead of bytes."
macro ballocs(args...)
    return esc(quote
        $BenchmarkTools.allocs($BenchmarkTools.minimum($BenchmarkTools.@benchmark $(args...)))
    end)
end

"Test passes if `args...` causes no heap allocations."
macro test_no_allocs(args...)
    return esc(quote
        # limit number of samples to 5... since the variance is extremely low
        @test (@ballocs $(args...) samples=3) == 0
    end)
end

@testset "Dep_Images" begin # dependencies on Images.jl
    img = rand(4, 3)
    img_rgb = rand(RGB, 4, 4)

    ex1 = [RGB(1,0,0) RGB(0,1,0) RGB(0,0,1);
        RGB(1,1,0) RGB(1,1,1) RGB(0,0,0)]

    ex2 = zeros(RGB, 2, 3)
    ex2[1,1] = RGB(1,0,0)
    ex2

    # The "Hello World" of graphics
    @fastmath function gradient(nx::Int, ny::Int)
        img = zeros(RGB, ny, nx)
        @inbounds for j in 1:nx
            for i in 1:ny # Julia is column-major, i.e. iterate 1 column at a time
                x = j; y = (ny-i);
                r = x/nx
                g = y/ny
                b = 0.2
                img[i,j] = RGB(r,g,b)
            end
        end
        img
    end

    gradient(20,10)
end

""" Temporary function to shoot rays through each pixel. Later replaced by `render`

    Args:
        scene: a function that takes a ray, returns the color of any object it hit
"""
function shoot_rays1(nx::Int, ny::Int, scene, ::Type{T}) where T
    lower_left_corner = SA[-2,-1,-1]
    horizontal = SA[4,0,0]
    vertical = SA[T(0),T(2),T(0)]
    origin = SA[T(0),T(0),T(0)]

    img = zeros(RGB{T}, ny, nx)
    Threads.@threads for j in 1:nx
        @inbounds for i in 1:ny # Julia is column-major, i.e. iterate 1 column at a time
            u = T(j/nx)
            v = T((ny-i)/ny) # Y-up!
            ray = Ray(origin, normalize(lower_left_corner + u*horizontal + v*vertical))
            #r = x/nx
            #g = y/ny
            #b = 0.2
            img[i,j] = rgb(scene(ray))
        end
    end
    img
end

# @fastmath speeds up a lot!
@inline @fastmath function hit_sphere1(center::Vec3{T}, radius::T, r::Ray{T}) where T
    oc = r.origin - center
    #a = r.dir ⋅ r.dir # unnecessary since `r.dir` is normalized
    a = 1
    b = 2oc ⋅ r.dir
    c = (oc ⋅ oc) - radius*radius
    discriminant = b*b - 4a*c
    discriminant > 0
end

@inline function sphere_scene1(r::Ray{T}) where T
    if hit_sphere1(SA[T(0), T(0), T(-1)], T(0.5), r) # sphere of radius 0.5 centered at z=-1
        return SA[T(1),T(0),T(0)] # red
    else
        skycolor(r)
    end
end

@inline function hit_sphere2(center::Vec3{T}, radius::T, r::Ray{T}) where T
    oc = r.origin - center
    #a = r.dir ⋅ r.dir # unnecessary since `r.dir` is normalized
    a = 1
    b = 2oc ⋅ r.dir
    c = (oc ⋅ oc) - radius^2
    discriminant = b*b - 4a*c
    if discriminant < 0
        return -1
    else
        return (-b - sqrt(discriminant)) / 2a
    end
end

@inline function sphere_scene2(r::Ray{T}) where T
    sphere_center = SA[T(0),T(0),T(-1)]
    t = hit_sphere2(sphere_center, T(0.5), r) # sphere of radius 0.5 centered at z=-1
    if t > T(0)
        n⃗ = normalize(point(r, t) - sphere_center) # normal vector. typed n\vec
        return T(0.5)*n⃗ + SA[T(0.5),T(0.5),T(0.5)] # remap normal to rgb
    else
        skycolor(r)
    end
end

@testset "RayTracingWeekend.jl" begin
    t_col = SA[0.4, 0.5, 0.1] # test color

    # Verify that no heap allocations occur for low-level functions
    @test_no_allocs squared_length($t_col) # 1.162 ns

    @test !near_zero(t_col)
    @test_no_allocs near_zero($t_col) # 1.382 ns

    #@show @btime rgb($t_col) # 1.382 ns
    @test rgb(t_col) == RGB{Float64}(0.4,0.5,0.1)
    @test_no_allocs rgb($t_col)

    _origin = SA[0.0,0.0,0.0]
    _v3_minusY = SA[0.0,-1.0,0.0]
    _t_ray1 = Ray(_origin, _v3_minusY)

    # Float32: 6.161 ns (0 allocations: 0 bytes)
    # Float64: 6.900 ns (0 allocations: 0 bytes)
    @test_no_allocs Ray($_origin, $_v3_minusY)
    @test_no_allocs point($_t_ray1, 0.5) # 1.412 ns
    @test_no_allocs skycolor($_t_ray1) # 1.402 ns
    @test_no_allocs rgb(skycolor($_t_ray1)) # 1.412 ns
    @test_no_allocs trand()
    @test_no_allocs random_between(50.0, 100.0) # 2.695 ns
    @test_no_allocs random_vec3(-1.0,1.0)
    @test_no_allocs random_vec2(-1.0f0,1.0f0)

    # REMEMBER: the times are somewhat random! e.g. use best timing of 5!
    @test_no_allocs random_vec3_in_sphere(Float64) # 19.716 ns
    @test_no_allocs random_vec3_on_sphere(Float64)
    @test_no_allocs random_vec2_in_disk(Float64)

    @benchmark random_vec2_in_disk(Float64) samples=3

    # varies depending on number of threads, maybe mem alignment?
    @test (@ballocs shoot_rays1(20, 10, $skycolor, $Float32) samples=3) < 85 # typically 81-83

    shoot_rays1(40, 20, sphere_scene1, Float32)

    @test (@ballocs shoot_rays1(40, 20, sphere_scene1, Float32) samples=3) < 85

    # # claforte: this got significantly worse after Vec3{T} was integrated...
    # # WAS:
    # # Float32: 262.286 μs (2 allocations: 234.45 KiB)
    # # NOW:
    # # Float32: 290.470 μs (2 allocations: 234.45 KiB)
    # # Float64: 351.164 μs (2 allocations: 468.83 KiB)
    # # Above was all using 1 single thread. With 16 threads:
    # #   64.612 μs (83 allocations: 241.75 KiB)
    # # Eliminate unnecessary r.dir ⋅ r.dir:
    # #   61.586 μs (82 allocations: 241.72 KiB)
    # #print("@btime main(200,100,sphere_scene2, Float32):")
    shoot_rays1(200,100,sphere_scene2, Float32)

    @assert RayTracingWeekend.reflect(SA[0.6,-0.8,0.0], SA[0.0,1.0,0.0]) == SA[0.6,0.8,0.0]

    default_camera(SA[0f0,0f0,0f0]) # Float32 camera

    get_ray(default_camera(SA[0f0,0f0,0f0]), 0.0f0, 0.0f0)

    @test_no_allocs get_ray(default_camera(SA[0f0,0f0,0f0]), 0.0f0, 0.0f0)

    ELEM_TYPE = Float64 # Float32/Float64

    t_default_cam = default_camera(SA{ELEM_TYPE}[0,0,0])

    # Using @paulmelis' style of hit(): @inbounds for i in eachindex(hittables) and Union{HitRecord, Nothing}
    #   951.447 μs (65574 allocations: 5.12 MiB)
    render(scene_2_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16) # 16 samples

    # # Using @paulmelis' style of hit(): @inbounds for i in eachindex(hittables) and Union{HitRecord, Nothing}
    # #   101.161 μs (4314 allocations: 459.41 KiB)
    # render(scene_2_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 1) # 1 sample

    # #render(scene_4_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

    # # unchanged angle
    # @assert refract((@SVector[0.6,-0.8,0]), (@SVector[0.0,1.0,0.0]), 1.0) == @SVector[0.6,-0.8,0.0]

    # # wider angle
    # t_refract_widerθ = refract(@SVector[0.6,-0.8,0.0], @SVector[0.0,1.0,0.0], 2.0)
    # @assert isapprox(t_refract_widerθ, @SVector[0.87519,-0.483779,0.0]; atol=1e-3)

    # # narrower angle
    # t_refract_narrowerθ = refract(@SVector[0.6,-0.8,0.0], @SVector[0.0,1.0,0.0], 0.5)
    # @assert isapprox(t_refract_narrowerθ, @SVector[0.3,-0.953939,0.0]; atol=1e-3)

    # #render(scene_diel_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)
    # #render(scene_diel_spheres(), default_camera(), 320, 32)

    # # Hollow Glass sphere using a negative radius
    # #render(scene_diel_spheres(-0.5; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

    # #render(scene_diel_spheres(; elem_type=ELEM_TYPE), default_camera((SA{ELEM_TYPE}[-2,2,1]), (SA{ELEM_TYPE}[0,0,-1]),
    # #																 (SA{ELEM_TYPE}[0,1,0]), ELEM_TYPE(20)), 96, 16)

    # #render(scene_blue_red_spheres(; elem_type=ELEM_TYPE), t_default_cam, 96, 16)

    # t_cam1 = default_camera([13,2,3], [0,0,0], [0,1,0], 20, 16/9, 0.1, 10.0; elem_type=ELEM_TYPE)

    # # took ~20s (really rough timing) in REPL, before optimization
    # # after optimization:
    # #   880.997 ms (48801164 allocations: 847.10 MiB)
    # # after switching to Float32 + reducing allocations using rand_vec3!():
    # #    79.574 ms (  542467 allocations:  14.93 MiB)
    # # after optimizing skycolor, rand*, probably more stuff I forgot...
    # #    38.242 ms (  235752 allocations:   6.37 MiB)
    # # after removing all remaining Color, Vec3, replacing them with @SVector[]...
    # #    26.790 ms (   91486 allocations: 2.28 MiB)
    # # Using convert(Float32, ...) instead of MyFloat(...):
    # #    25.856 ms (   70650 allocations: 1.96 MiB)
    # # Don't specify return value of Option{HitRecord} in hit()
    # # Don't initialize unnecessary elements in HitRecord():
    # #    38.961 ms (   70681 allocations: 1.96 MiB) # WORSE, probably because we're writing a lot more to stack...?
    # # Replace MyFloat by Float32:
    # #    36.726 ms (13889 allocations: 724.09 KiB)
    # # @inline lots of stuff:
    # #    14.690 ms (13652 allocations: 712.98 KiB)
    # # rand(Float32) to avoid Float64s:
    # #    14.659 ms (13670 allocations: 713.84 KiB)
    # # Re-measured:
    # #    14.069 ms (13677 allocations: 714.22 KiB)
    # # After parameterized Vec3{T}:
    # # Float32: 14.422 ms (26376 allocations: 896.52 KiB) # claforte: why are there 2X the allocations as previously?
    # # Float64: 12.772 ms (12868 allocations: 1.08 MiB) (10% speed-up! Thanks @woclass!)
    # # rand() using MersenneTwister _rng w/ Float64:
    # #    14.092 ms (12769 allocations: 1.07 MiB) (SLOWER?!)
    # # rand() using Xoroshiro128Plus _rng w/ Float64:
    # #    12.581 ms (12943 allocations: 1.08 MiB)
    # # Above was all using 1 single thread. With 16 threads:
    # #     1.789 ms (12926 allocations: 1.08 MiB) (WOW!)
    # # Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
    # # Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
    # #     2.168 ms (13791 allocations: 1.15 MiB)
    # # Using  bunch of @inbounds, @simd in low-level functions
    # #     2.076 ms (13861 allocations: 1.15 MiB)
    # # Lots of optimizations, up to `Using non-mutable HitRecord, Union{HitRecordMissing}, ismissing():`
    # #     2.042 ms (14825 allocations: 1.23 MiB)
    # #@btime render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 96, 1)

    # # took 5020s in Pluto.jl, before optimizations!
    # # after lots of optimizations, up to switching to Float32 + reducing allocations using rand_vec3!():
    # #   10.862032 seconds (70.09 M allocations: 1.858 GiB, 1.53% gc time)
    # # after optimizing skycolor, rand*, probably more stuff I forgot...
    # #    4.926 s (25,694,163 allocations: 660.06 MiB)
    # # after removing all remaining Color, Vec3, replacing them with @SVector[]...
    # #    3.541 s (9074843 allocations: 195.22 MiB)
    # # Using convert(Float32, ...) instead of MyFloat(...):
    # #    3.222 s (2055106 allocations: 88.53 MiB)
    # # Don't specify return value of Option{HitRecord} in hit()
    # # Don't initialize unnecessary elements in HitRecord():
    # #    5.016 s (2056817 allocations: 88.61 MiB) #  WORSE, probably because we're writing a lot more to stack...?
    # # Replace MyFloat by Float32:
    # #    5.185 s (1819234 allocations: 83.55 MiB) # Increase of 1% is probably noise
    # # Remove normalize() in reflect() (that function assumes the inputs are normalized)
    # #    5.040 s (1832052 allocations: 84.13 MiB)
    # # @inline lots of stuff:
    # #    2.110 s (1823044 allocations: 83.72 MiB)
    # # rand(Float32) to avoid Float64s:
    # #    2.063 s (1777985 allocations: 81.66 MiB)
    # # @woclass's "Use alias instead of new struct", i.e. `const HittableList = Vector{Hittable}`:
    # #    1.934 s (1796954 allocations: 82.53 MiB)
    # # @woclass's Vec3{T} with T=Float64: (7.8% speed-up!)
    # #    1.800 s (1711061 allocations: 131.03 MiB)
    # # rand() using Xoroshiro128Plus _rng w/ Float64:
    # #    1.808 s (1690104 allocations: 129.43 MiB) (i.e. rand is not a bottleneck?)
    # # Above was all using 1 single thread. With 16 threads:
    # #  265.331 ms (1645402 allocations: 126.02 MiB) (WOW!)
    # # Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
    # # Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
    # #  308.217 ms (1830162 allocations: 140.12 MiB)
    # # Using @inbounds, @simd in low-level functions:
    # #  302.952 ms (1892513 allocations: 144.88 MiB)
    # # Convert Camera and every Material structs to non-mutable:
    # #  301.042 ms (1849711 allocations: 141.61 MiB)  (i.e. - unchanged)
    # # Adapt @Christ_Foster's Base.getproperty w/ @inline @inbounds:
    # #  292.603 ms (1856398 allocations: 142.12 MiB) (ran multiple times, seems like real, 3-5% speed-up)
    # # Eliminate the off-by-half-a-pixel offset:
    # #  286.873 ms (1811412 allocations: 138.69 MiB) (ran multiple times, seems like ~2.5% speed-up)
    # # Fixed, per-thread RNGs with fixed seeds
    # #  286.575 ms (1884433 allocations: 144.26 MiB) (i.e. maybe a tiny bit faster considering this fixed seed has more allocations?)
    # # Make HitRecord non-mutable:
    # #   29.733 s (937962909 allocations: 69.88 GiB) (WTF!)
    # # Lots of optimizations, up to `Using non-mutable HitRecord, Union{HitRecordMissing}, ismissing():`
    # #  306.011 ms (1884433 allocations: 144.26 MiB) (Still slower... Hum)
    # # Using @paulmelis' style of hit(): @inbounds for i in eachindex(hittables) and Union{HitRecord, Nothing}
    # #  304.877 ms (1884433 allocations: 144.26 MiB)
    # # Extract the scene creation from the render() call:
    # #  300.344 ms (1883484 allocations: 144.21 MiB)
    # print("render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 200, 32):")
    # reseed!()
    # _scene_random_spheres = scene_random_spheres(; elem_type=ELEM_TYPE)
    # @btime render($_scene_random_spheres, $t_cam1, 200, 32)

    # # After some optimization, took ~5.6 hours:
    # #   20171.646846 seconds (94.73 G allocations: 2.496 TiB, 1.06% gc time)
    # # ... however the image looked weird... too blurry
    # # After removing all remaining Color, Vec3, replacing them with @SVector[]... took ~3.7 hours:
    # #   13326.770907 seconds (29.82 G allocations: 714.941 GiB, 0.36% gc time)
    # # Using convert(Float32, ...) instead of MyFloat(...):
    # # Don't specify return value of Option{HitRecord} in hit()
    # # Don't initialize unnecessary elements in HitRecord():
    # # Took ~4.1 hours:
    # #   14723.339976 seconds (5.45 G allocations: 243.044 GiB, 0.11% gc time) # WORSE, probably because we're writing a lot more to stack...?
    # # Replace MyFloat by Float32:
    # # Lots of other optimizations including @inline lots of stuff:
    # #    6018.101653 seconds (5.39 G allocations: 241.139 GiB, 0.21% gc time) (1.67 hour)
    # # @woclass's rand(Float32) to avoid Float64s: (expected to provide 2.2% speed up)
    # # @woclass's "Use alias instead of new struct", i.e. `const HittableList = Vector{Hittable}`
    # # @woclass's Vec3{T} with T=Float64: (7.8% speed-up!):
    # #    5268.175362 seconds (4.79 G allocations: 357.005 GiB, 0.47% gc time) (1.46 hours)
    # # Above was all using 1 single thread. With 16 threads: (~20 minutes)
    # #    1210.363539 seconds (4.94 G allocations: 368.435 GiB, 10.08% gc time)
    # # Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
    # # Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
    # #    1298.522674 seconds (5.43 G allocations: 404.519 GiB, 10.18% gc time)
    # # Using @inbounds, @simd in low-level functions:
    # #    1314.510565 seconds (5.53 G allocations: 411.753 GiB, 10.21% gc time) # NOTE: difference due to randomness?
    # # Adapt @Christ_Foster's Base.getproperty w/ @inline @inbounds: (expect 3-5% speed-up)
    # # Eliminate the off-by-half-a-pixel offset: (expect ~2.5% speed-up)
    # # Fixed, per-thread RNGs with fixed seeds (expecting no noticeable change in speed)
    # #  Using 16 threads: (21m22s)
    # #    1282.437499 seconds (5.53 G allocations: 411.742 GiB, 10.08% gc time) (i.e. 2.5% speed-up... currently GC- and memory-bound?)
    # #  Using 14 threads: (21m45s)
    # #    1305.767627 seconds (5.53 G allocations: 411.741 GiB, 9.97% gc time)
    # #print("@time render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 1920, 1000):")
    # #@time render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 1920, 1000)

    # t_cam2 = default_camera([3,3,2], [0,0,-1], [0,1,0], 20, 16/9, 2.0, norm([3,3,2]-[0,0,-1]);
    #                         elem_type=ELEM_TYPE)

    # # Before optimization:
    # #  5.993 s  (193097930 allocations: 11.92 GiB)
    # # after disabling: `Base.getproperty(vec::SVector{3}, sym::Symbol)`
    # #  1.001 s  ( 17406437 allocations: 425.87 MiB)
    # # after forcing Ray and point() to use Float64 instead of AbstractFloat:
    # #  397.905 ms (6269207 allocations: 201.30 MiB)
    # # after forcing use of Float32 instead of Float64:
    # #  487.680 ms (7128113 allocations: 196.89 MiB) # More allocations... something is causing them...
    # # after optimizing rand_vec3!/rand_vec2! to minimize allocations:
    # #  423.468 ms (6075725 allocations: 158.92 MiB)
    # # after optimizing skycolor, rand*, probably more stuff I forgot...
    # #  217.853 ms (2942272 allocations: 74.82 MiB)
    # # after removing all remaining Color, Vec3, replacing them with @SVector[]...
    # #   56.778 ms (1009344 allocations: 20.56 MiB)
    # # Using convert(Float32, ...) instead of MyFloat(...):
    # #   23.870 ms (210890 allocations: 8.37 MiB)
    # # Replace MyFloat by Float32:
    # #   22.390 ms (153918 allocations: 7.11 MiB)
    # # Various other changes, e.g. remove unnecessary normalize
    # #   20.241 ms (153792 allocations: 7.10 MiB)
    # # @inline lots of stuff:
    # #   18.065 ms (153849 allocations: 7.10 MiB)
    # # rand(Float32) to avoid Float64s:
    # #   16.035 ms (153777 allocations: 7.10 MiB)
    # # After @woclass's Vec3{T} with T=Float64:
    # #   16.822 ms (153591 allocations: 11.84 MiB)
    # # rand() using Xoroshiro128Plus _rng w/ Float64:
    # #   13.469 ms (153487 allocations: 11.83 MiB)
    # # Above was all using 1 single thread. With 16 threads:
    # #    6.537 ms (153599 allocations: 11.84 MiB)
    # # Above was all using max bounces=4, since this looked fine to me (except the negatively scaled sphere).
    # # Switching to max bounces=16 to match C++ version decreased performance by 7.2%:
    # #    6.766 ms (161000 allocations: 12.40 MiB)
    # # @inbounds and @simd in low-level functions
    # #    6.519 ms (160609 allocations: 12.37 MiB)
    # #render(scene_diel_spheres(; elem_type=ELEM_TYPE), t_cam2, 96, 16)

    # using Profile
    # render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 16, 1)
    # Profile.clear_malloc_data()
    # render(scene_random_spheres(; elem_type=ELEM_TYPE), t_cam1, 17, 13)

end # @testset


text:
using DensestCommonSubgraph
using FactCheck

# write your own tests here
facts("Testing DensestCommonSubgraph") do
  context("BinaryToGray conversion") do
    include("test_bitstr.jl")
    n = 5
    test_binary_to_gray(n)
  end

  context("DensityFunctions") do
    n = 1000

    p = 0.4
    k = 200
    G = erdos_renyi(n, p)
    for i=1:10
      ns = rand((n/10):n)
      S = randperm(n)[1:ns]
      d = subgraph_density(G, S)
      r = subgraph_relative_density(G, S)
      @fact d => roughly((ns-1)/2 * p; atol=1) "subgraph_density: $d expected: $((ns-1)/2 * p)"
      @fact r => roughly(p; atol=0.1) "relative_density: $r expected: $p"
    end

    H = planted_clique(n, p,  k)
    A = SparseMatrixCSC{Bool,Int64}[]
    push!(A, H)
    push!(A, G)
    for i=1:10
      ns = rand(1:n)
      S = randperm(n)[1:ns]
      dg = subgraph_density(G, S)
      dh = subgraph_density(H, S)
      da = common_subgraph_density(A, S)
      # println("dg: $dg dh: $dh dcs: $da")
      @fact da => roughly(min(dh, dg))
    end

  end

#   context("Faster DCS LP") do
#     n= 30
#     p= 1/n
#     for i=1:10
#       G  = get_dcs_test_instance(n, p);
#       @time p2 = dcs_lp(G);
#       @time p3 = DensestCommonSubgraph.faster_dcs_lp(G);
#       @time p1 = DensestCommonSubgraph.old_dcs_lp(G);
#       println("OLD LP: $(p1.optval), CURRENT: $(p2.optval) FASTER: $(p3.optval)")
#       @fact p1.optval => roughly(p2.optval, atol=1e-2) "Error in fast LP!"
#     end
#   end

  end

end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct PSSSimple <: PSS
        K_ω::Float64
        K_p::Float64
        ext::Dict{String, Any}
        states::Vector{Symbol}
        n_states::Int
        internal::InfrastructureSystemsInternal
    end

Parameters of a PSS that returns a proportional droop voltage to add to the reference for the AVR

# Arguments
- `K_ω::Float64`: Proportional gain for frequency, validation range: `(0, nothing)`
- `K_p::Float64`: Proportional gain for active power, validation range: `(0, nothing)`
- `ext::Dict{String, Any}`
- `states::Vector{Symbol}`
- `n_states::Int`: PSSSimple has no states
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct PSSSimple <: PSS
    "Proportional gain for frequency"
    K_ω::Float64
    "Proportional gain for active power"
    K_p::Float64
    ext::Dict{String, Any}
    states::Vector{Symbol}
    "PSSSimple has no states"
    n_states::Int
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function PSSSimple(K_ω, K_p, ext=Dict{String, Any}(), )
    PSSSimple(K_ω, K_p, ext, Vector{Symbol}(), 0, InfrastructureSystemsInternal(), )
end

function PSSSimple(; K_ω, K_p, ext=Dict{String, Any}(), states=Vector{Symbol}(), n_states=0, internal=InfrastructureSystemsInternal(), )
    PSSSimple(K_ω, K_p, ext, states, n_states, internal, )
end

# Constructor for demo purposes; non-functional.
function PSSSimple(::Nothing)
    PSSSimple(;
        K_ω=0,
        K_p=0,
        ext=Dict{String, Any}(),
    )
end

"""Get [`PSSSimple`](@ref) `K_ω`."""
get_K_ω(value::PSSSimple) = value.K_ω
"""Get [`PSSSimple`](@ref) `K_p`."""
get_K_p(value::PSSSimple) = value.K_p
"""Get [`PSSSimple`](@ref) `ext`."""
get_ext(value::PSSSimple) = value.ext
"""Get [`PSSSimple`](@ref) `states`."""
get_states(value::PSSSimple) = value.states
"""Get [`PSSSimple`](@ref) `n_states`."""
get_n_states(value::PSSSimple) = value.n_states
"""Get [`PSSSimple`](@ref) `internal`."""
get_internal(value::PSSSimple) = value.internal

"""Set [`PSSSimple`](@ref) `K_ω`."""
set_K_ω!(value::PSSSimple, val) = value.K_ω = val
"""Set [`PSSSimple`](@ref) `K_p`."""
set_K_p!(value::PSSSimple, val) = value.K_p = val
"""Set [`PSSSimple`](@ref) `ext`."""
set_ext!(value::PSSSimple, val) = value.ext = val


text:
# This is an approximate model of an acceleration controlled differential drive.
# The angular velocities and the orientation are updated using their closed form
# expression, the position along x and y axes are updated with the midway constant
# angular velocities. The error in this heuristic increases with time and
# with higher curvature.

L = 0.13    # m : distance between the wheels
r = 0.033   # m : wheel radius

# Creates a new state based upon the previous one assuming that the
# specified control that was executed continuously from the beginning
# of the time slot for the specified time.
# Units: x : m; y : m; θ : rad; ω : rad/s; α : rad/s²; t : s
function new_state(state, control, elapsed_time)
    x, y, θ, ω_left, ω_right = state
    α_left, α_right = control

    ω_left⁺ = ω_left + Δω(α_left, elapsed_time)
    ω_right⁺ = ω_right + Δω(α_right, elapsed_time)
    θ⁺ = θ + Δθ(ω_left, α_left, ω_right, α_right, elapsed_time)
    x⁺ = x + approximateΔx(θ, (ω_left + ω_left⁺)/2, (ω_right + ω_right⁺)/2, elapsed_time)
    y⁺ = y + approximateΔy(θ, (ω_left + ω_left⁺)/2, (ω_right + ω_right⁺)/2, elapsed_time)

    return [x⁺ y⁺ θ⁺ ω_left⁺ ω_right⁺]
end

# Change in the x position assuming constant angular velocities.
# Units: θ : rad; ω : rad/s; t : s
function approximateΔx(θ_init, ω_left, ω_right, elapsed_time)
    if ω_left != ω_right
        L/2.0 * (ω_left + ω_right)/(ω_right - ω_left) * (sin(θ_init + r/L * (ω_right - ω_left) * elapsed_time) - sin(θ_init))
    else
        r/2.0 * (ω_left + ω_right) * cos(θ_init) * elapsed_time
    end
end

# Change in the y position assuming constant angular velocities.
# Units: θ : rad; ω : rad/s; t : s
function approximateΔy(θ_init, ω_left, ω_right, elapsed_time)
    if ω_left != ω_right
        L/2.0 * (ω_left + ω_right)/(ω_right - ω_left) * (cos(θ_init) - cos(θ_init + r/L * (ω_right - ω_left) * elapsed_time))
    else
        r/2.0 * (ω_left + ω_right) * sin(θ_init) * elapsed_time
    end
end

# Change in the orientation after the specified time.
# Units: ω : rad/s; α : rad/s²; t : s
Δθ(ω_left_init, α_left, ω_right_init, α_right, elapsed_time) = r/L * (ω_right_init - ω_left_init) * elapsed_time + r/L * (α_right - α_left) * elapsed_time^2 / 2.0

# Change in the angular velocity over time.
# Units: α : rad/s²; t : s
Δω(α, elapsed_time) = α * elapsed_time

# Creates a new state given with the specified x, y and orientation
# components and the left and right wheel angular velocities.
# The position components should be in identical coordinate systems,
# e.g. world.
# Units: x : m; y : m; θ : rad; ω : rad/s
initial_state(x_world, y_world, θ_world, ω_left, ω_right) = [x_world y_world θ_world ω_left ω_right]

max_stars_repo_name: grgomrton/robot-model

max_issues_repo_name: grgomrton/robot-model

max_forks_repo_name: grgomrton/robot-model


text:
using ClimaTimeSteppers, Test

include("problems.jl")

@testset "saving default" begin
  sol = solve(deepcopy(const_prob), SSPRK33ShuOsher(); dt=1/32)
  ts = [const_prob.tspan[end]]
  @test sol.t == ts
  @test sol.u ≈ [const_sol(const_prob.u0, const_prob.p, t) for t in sol.t]
end
@testset "saveat isa Number" begin
  sol = solve(deepcopy(const_prob), SSPRK33ShuOsher(); dt=1/32, saveat=0.1)
  ts = 0:0.1:1
  @test all(abs.(sol.t .- ts) .< 1/32)
  @test all(isinteger, sol.t .* 32)
  @test sol.u ≈ [const_sol(const_prob.u0, const_prob.p, t) for t in sol.t]
end
@testset "saveat isa Range" begin
  ts = 0:0.1:1
  sol = solve(deepcopy(const_prob), SSPRK33ShuOsher(); dt=1/32, saveat=0:0.1:1)
  @test all(abs.(sol.t .- ts) .< 1/32)
  @test all(isinteger, sol.t .* 32)
  @test sol.u ≈ [const_sol(const_prob.u0, const_prob.p, t) for t in sol.t]
end
@testset "save_everystep" begin
  sol = solve(deepcopy(const_prob), SSPRK33ShuOsher(); dt=1/32, save_everystep=true)
  ts = 0:1/32:1
  @test sol.t == ts
  @test sol.u ≈ [const_sol(const_prob.u0, const_prob.p, t) for t in sol.t]
end
@testset "saveat too many" begin
  sol = solve(deepcopy(const_prob), SSPRK33ShuOsher(); dt=1/32, saveat=0:0.01:1)
  ts = 0:1/32:1
  @test sol.t == ts
  @test sol.u ≈ [const_sol(const_prob.u0, const_prob.p, t) for t in sol.t]
end


text:
# add one unit to the current limit
# SI+
export psp_incrementilimit

const INCREMENT_ILIMIT = [0x53, 0x49, 0x2B, 0x0D]

"""
    psp_incrementilimit(io)

Add one unit to the current limit.
"""
function psp_incrementilimit(io_psp::IO)
  write(io_psp, INCREMENT_ILIMIT)
  return nothing
end


text:
@testset "observable" begin

nsamples=10^6
ob = Observable()
v = randn(nsamples)
for x in v
    ob &= x
end

@test isapprox(err(ob),  std(v) / √(nsamples-1), atol=1e-9)
@test isapprox(mean(ob), sum_kbn(v) / nsamples, atol=1e-14)
@test isapprox(mean(ob), 0, atol=5*err(ob))

ob = Observable()
for i=1:nsamples
    ob &= norm(2 .* rand(2) .- 1) < 1 ? 1 : 0
end
ob *= 4
@test isapprox(mean(ob), π, atol = 5err(ob))

@test ob == ob

ob2 = copy(ob)
ob2 &= NaN
@test ob2 == ob
ob2 &= Inf
@test ob2 == ob

end #testset

max_stars_repo_name: JuliaTagBot/LittleScienceTools

max_issues_repo_name: JuliaTagBot/LittleScienceTools

max_forks_repo_name: JuliaTagBot/LittleScienceTools


text:
using FiniteDiff
using ComponentArrays: ComponentVector
using NQCDynamics: DynamicsMethods, DynamicsUtils, nbeads, Estimators

get_blank(sim::Simulation) = randn(size(sim))
get_blank(sim::RingPolymerSimulation) = RingPolymerArray(randn(size(sim)))

function test_motion!(sim)
    f(x) = DynamicsUtils.classical_hamiltonian(sim, x)

    v = get_blank(sim)
    r = get_blank(sim)
    u = ComponentVector(v=v, r=r)
    du = zero(u)

    grad = FiniteDiff.finite_difference_gradient(f, u)

    DynamicsMethods.motion!(du, u, sim, 0.0)

    # Rdot = dH/dP = dH/dV / mass
    @test DynamicsUtils.get_positions(du) ≈ DynamicsUtils.get_velocities(grad) ./ sim.atoms.masses' rtol=1e-3
    # Vdot = Pdot / mass = -dH/dR / mass
    @test DynamicsUtils.get_velocities(du) ≈ -DynamicsUtils.get_positions(grad) ./ sim.atoms.masses' rtol=1e-3
end

function test_velocity!(sim)
    r = get_blank(sim)
    v = get_blank(sim)
    du = zero(v)

    DynamicsUtils.velocity!(du, v, r, sim, 0.0)

    @test du ≈ v
end

function test_acceleration!(sim)
    f(x) = DynamicsUtils.classical_potential_energy(sim, x)

    r = get_blank(sim)
    v = get_blank(sim)
    dv = zero(v)

    grad = FiniteDiff.finite_difference_gradient(f, r)

    DynamicsMethods.ClassicalMethods.acceleration!(dv, v, r, sim, 0.0)

    @test dv ≈ -grad ./ sim.atoms.masses'
end

function test_acceleration!(sim::RingPolymerSimulation)
    f(x) = DynamicsUtils.classical_potential_energy(sim, x)

    r = get_blank(sim)
    v = get_blank(sim)
    dv = zero(v)

    grad = FiniteDiff.finite_difference_gradient(f, r)

    DynamicsMethods.ClassicalMethods.ring_polymer_acceleration!(dv, v, r, sim, 0.0)

    @test dv ≈ -grad ./ sim.atoms.masses'
end


text:
using ROOTDataFrames, DataFrames
using Base.Test

N = 500000
m = 25

df = DataFrame()

for i=1:m
    df[symbol("col$i")] = Float64[j for j=1:N]
end

for i=m+1:m+10
    df[symbol("dyncol$i")] = Vector{Float64}[rand(10) for j=1:N]
end
#
# for i=m+1:m+10
#     df[symbol("dyncol$i")] = Vector{Float64}[rand(10) for j=1:N]
# end

println("writing DataFrame")
writetree("test.root", df)


text:
# https://www.rp-photonics.com/brewsters_angle.html

# Load modules
using Plots, LaTeXStrings
pyplot()
using ThinFilmsTools

function main()
    # Define beam
    λ = [1064.]
    θ = LinRange(0.01,90,901)
    beam = PlaneWave(λ,θ)
    # Define layers: notice that lambda is outside the range
    layers = [
        LayerTMMO(RIdb.dummy(beam.λ,1.0,0.0)),
        LayerTMMO(RIdb.dummy(beam.λ,1.5,0.0), d=10.),
        LayerTMMO(RIdb.dummy(beam.λ,1.5,0.0)),
    ]
    # call main script
    sol = tmm_optics(beam, layers)
    return sol
end

sol = main()

# plot the R, T and A spectra
plot(SpectrumAngle1D(),
    sol.beam.θ,
    [sol.Spectra.Rp[1,:], sol.Spectra.Rs[1,:], sol.Spectra.Rs[1,:]./sol.Spectra.Rp[1,:]./1000.],
    label=["p-wave" "s-wave" "(Rs/Rp)/1000"],
    line=([:solid :dash :dashdot]),
    xlims=(sol.beam.θ[1], sol.beam.θ[end]),
    yaxis=("Reflectance", (0.,0.2)),
)
gui()

# plot the refractive index profile
plot(RIprofile(), sol)
gui()


text:
##############################################################################
#
# Test problem
#
#############################################################################

function rosenbrock(x)
  n = 6

  # Initializations
  f = 0

  evenIdx = 2:2:n
  oddIdx = 1:2:(n - 1)

  f1 = x[evenIdx] .- x[oddIdx] .^ 2
  f2 = 1 .- x[oddIdx]

  # Function
  f = sum(f1 .^ 2 .+ f2 .^ 2)

  return f
end


text:
function parseInput(raw_input)
    lines = split(raw_input, '\n')

    firstLine = split(lines[1], ' ')

    numberOfVertices = parse(Int, firstLine[1])
    numberOfEdges = parse(Int, firstLine[2])

    edges = Array{Union{Edge, Nothing}}(nothing, numberOfEdges)

    @assert numberOfEdges + 2 <= length(lines)

    for i in 1:numberOfEdges
        edgeArray = split(lines[i+1], ' ')
        edge = Edge(parse(Int, edgeArray[1])+1, parse(Int, edgeArray[2])+1)
        edges[i] = edge
    end

    return InputData(edges, numberOfEdges, numberOfVertices)
end

function solutionString(solution::OutputData)
    toReturn = ""
    toReturn *= string(solution.numberOfColors) * " " * string(solution.optimality ? 1 : 0) * "\n"
    for color in solution.edgeColors
        toReturn *= string(color) * " "
    end
    toReturn *= "\n"
end

function printSolution(solution::OutputData)
    println(solution.numberOfColors, " ", solution.optimality ? 1 : 0)
    for color in solution.edgeColors
        print(color, " ")
    end
    println()
end

function writeSolution(solution::OutputData, filename::String)
    open(filename, "w") do openedFile
        write(openedFile, solutionString(solution))
    end
    return
end

function getInputData(filename)
    if filename == ""
        throw(ArgumentError("You must specify a data file"))
    end

    input = InputData([], 0, 0)

    open(filename, "r") do openedFile
        input = parseInput(read(openedFile, String))
    end

    return input
end

max_stars_repo_name: corail-research/SeaPearlZoo

max_issues_repo_name: corail-research/SeaPearlZoo

max_forks_repo_name: corail-research/SeaPearlZoo


text:
using Test, Random, Query, StructArrays, Tables, StaticArrays
using GeometryBasics
using LinearAlgebra

@testset "GeometryBasics" begin
    @testset "embedding metadata" begin
        @testset "Meshes" begin

            @testset "per vertex attributes" begin
                points = rand(Point{3, Float64}, 8)
                tfaces = TetrahedronFace{Int}[(1, 2, 3, 4), (5, 6, 7, 8)]
                normals = rand(SVector{3, Float64}, 8)
                stress = LinRange(0, 1, 8)
                mesh = Mesh(meta(points, normals = normals, stress = stress), tfaces)

                @test hasproperty(coordinates(mesh), :stress)
                @test hasproperty(coordinates(mesh), :normals)
                @test coordinates(mesh).stress === stress
                @test coordinates(mesh).normals === normals
                @test coordinates(mesh).normals === normals
                @test GeometryBasics.faces(mesh) === tfaces
                @test propertynames(coordinates(mesh)) == (:position, :normals, :stress)

            end

            @testset "per face attributes" begin

                # Construct a cube out of Quads
                points = Point{3, Float64}[
                    (0.0, 0.0, 0.0), (2.0, 0.0, 0.0),
                    (2.0, 2.0, 0.0), (0.0, 2.0, 0.0),
                    (0.0, 0.0, 12.0), (2.0, 0.0, 12.0),
                    (2.0, 2.0, 12.0), (0.0, 2.0, 12.0)
                ]

                facets = QuadFace{Cint}[
                    1:4,
                    5:8,
                    [1,5,6,2],
                    [2,6,7,3],
                    [3, 7, 8, 4],
                    [4, 8, 5, 1]
                ]

                markers = Cint[-1, -2, 0, 0, 0, 0]
                # attach some additional information to our faces!
                mesh = Mesh(points, meta(facets, markers = markers))
                @test hasproperty(GeometryBasics.faces(mesh), :markers)
                # test with === to assert we're not doing any copies
                @test GeometryBasics.faces(mesh).markers === markers
                @test coordinates(mesh) === points
                @test metafree(GeometryBasics.faces(mesh)) === facets

            end

        end
        @testset "polygon with metadata" begin
            polys = [Polygon(rand(Point{2, Float32}, 20)) for i in 1:10]
            pnames = [randstring(4) for i in 1:10]
            numbers = LinRange(0.0, 1.0, 10)
            bin = rand(Bool, 10)
            # create just an array
            plain = meta(polys, name = pnames, value = numbers, category = bin)
            # create a MultiPolygon with the right type & meta information!
            multipoly = MultiPolygon(polys; name = pnames, value = numbers, category = bin)
            for x in (plain, multipoly)
                for (mp, p, n, num, b) in zip(x, polys, pnames, numbers, bin)
                    @test mp.polygon == p
                    @test mp.name == n
                    @test mp.value == num
                    @test mp.category == b
                end

                filtered = @from i in x begin
                    @where i.value < 0.7
                    @select i
                    @collect
                end
                @test length(filtered) == 7
            end
        end

        @testset "point with metadata" begin
            p = Point(1.1, 2.2)
            @test p isa AbstractVector{Float64}
            pm = GeometryBasics.PointMeta(1.1, 2.2; a=1, b=2)
            @test meta(pm) === (a=1, b=2)
            @test metafree(pm) === p
        end
    end

    @testset "view" begin
        @testset "TupleView" begin
            x = [1, 2, 3, 4, 5, 6]
            y = TupleView{2, 1}(x)
            @test y == [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6)]

            y = TupleView{2}(x)
            @test y == [(1, 2), (3, 4), (5, 6)]

            y = TupleView{2, 3}(x)
            @test y == [(1, 2), (4, 5)]

            y = TupleView{3, 1}(x)
            @test y == [(1, 2, 3), (2, 3, 4), (3, 4, 5), (4, 5, 6)]

            y = TupleView{2, 1}(x, connect = true)
            @test y == [(1, 2), (2, 3), (3, 4), (4, 5), (5, 6), (6, 1)]

        end

        @testset "connected views" begin
            numbers = [1, 2, 3, 4, 5, 6]
            x = connect(numbers, Point{2})

            @test x == Point[(1, 2), (3, 4), (5, 6)]

            line = connect(x, Line, 1)
            @test line == [Line(Point(1, 2), Point(3, 4)), Line(Point(3, 4), Point(5, 6))]

            triangles = connect(x, Triangle)
            @test triangles == [Triangle(Point(1, 2), Point(3, 4), Point(5, 6))]
            x = connect([1, 2, 3, 4, 5, 6, 7, 8], Point{2})
            tetrahedra = connect(x, NSimplex{4})
            @test tetrahedra == [Tetrahedron(x[1], x[2], x[3], x[4])]

            @testset "matrix non-copy point views" begin
                # point in row
                points = [1 2; 1 4; 66 77]
                comparison = [Point(1, 2), Point(1, 4), Point(66, 77)]
                @test connect(points, Point{2}) == comparison
                # point in column
                points = [1 1 66; 2 4 77]
                # huh, reinterpret array doesn't seem to like `==`
                @test all(((a,b),)-> a==b, zip(connect(points, Point{2}), comparison))
            end
        end

        @testset "face views" begin
            numbers = [1, 2, 3, 4, 5, 6]
            points = connect(numbers, Point{2})
            faces = connect([1, 2, 3], TriangleFace)
            triangles = connect(points, faces)
            @test triangles == [Triangle(Point(1, 2), Point(3, 4), Point(5, 6))]
            x = Point{3}(1.0)
            triangles = connect([x], [TriangleFace(1, 1, 1)])
            @test triangles == [Triangle(x, x, x)]
            points = connect([1, 2, 3, 4, 5, 6, 7, 8], Point{2})
            faces = connect([1, 2, 3, 4], SimplexFace{4})
            triangles = connect(points, faces)
            @test triangles == [Tetrahedron(points...)]
        end

        @testset "reinterpret" begin
            numbers = collect(reshape(1:6, 2, 3))
            points = reinterpret(Point{2, Int}, numbers)
            @test points[1] === Point(1, 2)
            @test points[2] === Point(3, 4)
            numbers[4] = 0
            @test points[2] === Point(3, 0)
        end
    end

    @testset "constructors" begin
        @testset "LineFace" begin

            points = connect([1, 2, 3, 4, 5, 6], Point{2})
            linestring = LineString(points)
            @test linestring == [Line(points[1], points[2]), Line(points[2], points[3])]

            points = rand(Point{2, Float64}, 4)
            linestring = LineString(points, 2)
            @test linestring == [Line(points[1], points[2]), Line(points[3], points[4])]

            linestring = LineString([points[1] => points[2], points[2] => points[3]])
            @test linestring == [Line(points[1], points[2]), Line(points[2], points[3])]

            faces = [1, 2, 3]
            linestring = LineString(points, faces)
            @test linestring == LineString([points[1] => points[2], points[2] => points[3]])
            a, b, c, d = Point(1, 2), Point(3, 4), Point(5, 6), Point(7, 8)
            points = [a, b, c, d]; faces = [1, 2, 3, 4]
            linestring = LineString(points, faces, 2)
            @test linestring == LineString([a => b, c => d])

            faces = [LineFace(1, 2), LineFace(3, 4)]
            linestring = LineString(points, faces)
            @test linestring == LineString([a => b, c => d])
        end

        @testset "Polygon" begin

            points = connect([1, 2, 3, 4, 5, 6], Point{2})
            polygon = Polygon(points)
            @test polygon == Polygon(LineString(points))

            points = rand(Point{2, Float64}, 4)
            linestring = LineString(points, 2)
            @test Polygon(points, 2) == Polygon(linestring)

            faces = [1, 2, 3]
            polygon = Polygon(points, faces)
            @test polygon == Polygon(LineString(points, faces))

            a, b, c, d = Point(1, 2), Point(3, 4), Point(5, 6), Point(7, 8)
            points = [a, b, c, d]; faces = [1, 2, 3, 4]
            polygon = Polygon(points, faces, 2)
            @test polygon == Polygon(LineString(points, faces, 2))

            faces = [LineFace(1, 2), LineFace(3, 4)]
            polygon = Polygon(points, faces)
            @test polygon == Polygon(LineString(points, faces))
        end

        @testset "Mesh" begin

            numbers = [1, 2, 3, 4, 5, 6]
            points = connect(numbers, Point{2})

            mesh = Mesh(points, [1,2,3])
            @test mesh == [Triangle(points...)]

            x = Point{3}(1.0)
            mesh = Mesh([x], [TriangleFace(1, 1, 1)])
            @test mesh == [Triangle(x, x, x)]

            points = connect([1, 2, 3, 4, 5, 6, 7, 8], Point{2})
            faces = connect([1, 2, 3, 4], SimplexFace{4})
            mesh = Mesh(points, faces)
            @test mesh == [Tetrahedron(points...)]

            points = rand(Point3f0, 8)
            tfaces = [GLTriangleFace(1, 2, 3), GLTriangleFace(5, 6, 7)]
            normals = rand(Vec3f0, 8)
            uv = rand(Vec2f0, 8)
            mesh = Mesh(points, tfaces)
            meshuv = Mesh(meta(points; uv=uv), tfaces)
            meshuvnormal = Mesh(meta(points; normals=normals, uv=uv), tfaces)

            @test mesh isa GLPlainMesh
            @test meshuv isa GLUVMesh3D
            @test meshuvnormal isa GLNormalUVMesh3D

        end

        @testset "Multi geometries" begin
            # coordinates from https://en.wikipedia.org/wiki/Well-known_text_representation_of_geometry#Geometric_objects
            points = Point{2, Int}[(10, 40), (40, 30), (20, 20), (30, 10)]
            multipoint = MultiPoint(points)
            @test size(multipoint) === size(points)
            @test multipoint[3] === points[3]

            linestring1 = LineString(Point{2, Int}[(10, 10), (20, 20), (10, 40)])
            linestring2 = LineString(Point{2, Int}[(40, 40), (30, 30), (40, 20), (30, 10)])
            multilinestring = MultiLineString([linestring1, linestring2])
            @test size(multilinestring) === (2,)
            @test multilinestring[1] === linestring1
            @test multilinestring[2] === linestring2

            polygon11 = Polygon(Point{2, Int}[(30, 20), (45, 40), (10, 40), (30, 20)])
            polygon12 = Polygon(Point{2, Int}[(15, 5), (40, 10), (10, 20), (5, 10), (15, 5)])
            multipolygon1 = MultiPolygon([polygon11, polygon12])
            @test size(multipolygon1) === (2,)
            @test multipolygon1[1] === polygon11
            @test multipolygon1[2] === polygon12

            polygon21 = Polygon(Point{2, Int}[(40, 40), (20, 45), (45, 30), (40, 40)])
            polygon22 = Polygon(LineString(Point{2, Int}[(20, 35), (10, 30), (10, 10), (30, 5), (45, 20), (20, 35)]),
                [LineString(Point{2, Int}[(30, 20), (20, 15), (20, 25), (30, 20)])])
            multipolygon2 = MultiPolygon([polygon21, polygon22])
            @test size(multipolygon2) === (2,)
            @test multipolygon2[1] === polygon21
            @test multipolygon2[2] === polygon22
        end

    end
end

@testset "decompose/triangulation" begin
    primitive = Sphere(Point3f0(0), 1)
    mesh = triangle_mesh(primitive)
    @test decompose(Point, mesh) isa Vector{Point3f0}
    @test decompose(Point, primitive) isa Vector{Point3f0}
    primitive = Rect2D(0, 0, 1, 1)
    mesh = triangle_mesh(primitive)

    @test decompose(Point, mesh) isa Vector{Point2f0}
    @test decompose(Point, primitive) isa Vector{Point2f0}

    primitive = Rect3D(0, 0, 0, 1, 1, 1)
    triangle_mesh(primitive)

    primitive = Sphere(Point3f0(0), 1)
    normal_mesh(primitive)
    primitive = Rect2D(0, 0, 1, 1)
    normal_mesh(primitive)
    primitive = Rect3D(0, 0, 0, 1, 1, 1)
    normal_mesh(primitive)

    points = decompose(Point2f0, Circle(Point2f0(0), 1))
    triangle_mesh(points)
    @test true # yay no errors so far!

end

@testset "Tests from GeometryTypes" begin
    include("geometrytypes.jl")
end

@testset "modifying meta" begin
    xx = rand(10)
    points = rand(Point3f0, 10)
    m = GeometryBasics.Mesh(meta(points, xx=xx), GLTriangleFace[(1,2,3), (3,4,5)])
    color = rand(10)
    m = GeometryBasics.pointmeta(m; color=color)

    @test hasproperty(m, :xx)
    @test hasproperty(m, :color)

    @test m.xx === xx
    @test m.color === color

    m, colpopt = GeometryBasics.pop_pointmeta(m, :color)
    m, xxpopt = GeometryBasics.pop_pointmeta(m, :xx)

    @test propertynames(m) == (:position,)
    @test colpopt === color
    @test xxpopt === xx

    @testset "creating meta" begin
        x = Point3f0[(1,3,4)]
        # no meta gets added, so should stay the same
        @test meta(x) === x
        @test meta(x, value=[1]).position === x
    end
    pos = Point2f0[(10, 2)]
    m = Mesh(meta(pos, uv=[Vec2f0(1, 1)]), [GLTriangleFace(1, 1, 1)])
    @test m.position === pos
end


text:
using StaticArrays
@testset "Spectral radius" begin
    @test_throws ArgumentError ρ(zeros(0, 0))
    # There is a specialized function for StaticArrays
    @test ρ(@SMatrix [1 2; 3 4]) ≈ ρ([1 2; 3 4])
    @test ρ(-@SMatrix [1 2; 3 4]) ≈ ρ(-[1 2; 3 4])
    @test ρ(@SMatrix [1 1; 0 1]) ≈ ρ([1 1; 0 1])
    # Matrix A1 from [Section 4, PEDJ16]
    A = @SMatrix [0.94  0.56; -0.35  0.73]
    @test ρ(A) ≈ ρ(Matrix(A))
    @test ρ(-A) ≈ ρ(Matrix(-A))
end


text:
using Test

withenv("BLINK_PORT"=>4321) do
    p = readlines(`$(joinpath(Sys.BINDIR, "julia")) -E 'using Blink; println(Blink.port[])'`)
    @test p[1] == "4321"
end


text:
__precompile__()

module TrussMorph

# dependencies
using JSON
using Makie
using LinearAlgebra

# abstract types

# exported APIs
export Truss, MaterialProperties, SectionProperties
export parse_truss_json, parse_load_json, draw_truss!, draw_load!, draw_deformed!
export compute_morph_path

export assemble_global_stiffness_matrix, assemble_load_vector, dof_permutation, get_weight_calculation_fn

include("common_types.jl")
include("json_parser.jl")
include("draw_utils.jl")
include("stiffness.jl")
include("morph.jl")
include("utils.jl")

end # module


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

using Markdown
using InteractiveUtils

# ╔═╡ 2ea7d655-35ed-4687-9873-e658528232a8
include("..\\..\\jl\\bin2int.jl")

# ╔═╡ d7b19f60-6348-11ec-11b8-abe5c0a6cee7
md"""
# Day 3: Binary Diagnostic

**Reference:** [Day 3 (Advent of Code 2021)](https://adventofcode.com/2021/day/3)

**Source:** [ljk233/AdventOfCode2021/days/03/d3.jl](https://github.com/ljk233/AdventOfCode2021/blob/main/days/03/d3.jl)
"""

# ╔═╡ 14d4296b-4c17-47b5-80f2-1dcf89690f86
md"""
## Imports
"""

# ╔═╡ d9e09af7-42a2-4981-ab2f-21e7308c45e4
md"""
## Functions
"""

# ╔═╡ 8ceda75d-d741-4c1a-a884-dba74255270a
"""
	balanceof(lines: String[]) -> Integer
"""
function balanceof(lines, pos)
	b = 0
	for line ∈ lines
		line[pos] == '1' ? b += 1 : b -= 1
	end
	b ≥ 0 ? 1 : -1
end

# ╔═╡ aeba97ec-fd96-4ef4-99b0-726665015217
md"""
## Part 1

Use the binary numbers in your diagnostic report to calculate the gamma rate and epsilon rate, then multiply them together.

**What is the power consumption of the submarine?**
"""

# ╔═╡ 7f03ad5e-bf73-4b16-958b-592ba4b6d923
"""
	main1(f: String) -> Integer
"""
function main1(f)
	lines = readlines(f)
	gamma, epsilon = "", ""
	for pos ∈ 1:length(first(lines))
		if balanceof(lines, pos) == 1
			gamma *= "1"
			epsilon *= "0"
		else
			gamma *= "0"
			epsilon *= "1"
		end
	end
	return bin2int(gamma) * bin2int(epsilon)
end

# ╔═╡ 8cc55754-6815-4c08-8538-b6f4c74d5c28
md"""
**Test #1.**
Expected result = 198.
"""

# ╔═╡ d0d5ecd6-fa56-4fca-a293-fb0a682e4a7b
main1("test.in")

# ╔═╡ 63a48804-1a55-48c6-832f-74f52ac1d583
md"""
**Solution.**
Confirmed result = 3320834.
"""

# ╔═╡ 6c58e558-21fe-4221-b324-1efd90c7911a
main1("data.in")

# ╔═╡ 0e681b8a-82bc-425c-a374-e36d360559d2
md"""
## Part 2

Use the binary numbers in your diagnostic report to calculate the oxygen generator rating and CO2 scrubber rating, then multiply them together.

**What is the life support rating of the submarine?**
"""

# ╔═╡ 5012e220-d2e2-436c-9ce4-dfe34f365df5
"""
	rating(lines: String[], iscommon: Bool) -> String
"""
function rating(lines, iscommon)::String
	for pos ∈ 1:length(lines[1])
		m = balanceof(lines, pos)
		if !iscommon
			m *= -1
		end
		m == 1 ? target = '1' : target = '0'
		lines = filter(line -> line[pos] == target, lines)
		length(lines) == 1 ? (return lines[1]) : nothing
	end
end

# ╔═╡ 774f7d4b-9d47-41eb-a6d7-e9d0fc579fbb
"""
	main2(f: String) -> Integer
"""
function main2(f)
	lines::Vector{String} = readlines(f)
	o2 = rating(lines, true) |> bin2int
	co2 = rating(lines, false) |> bin2int
	return o2*co2
end

# ╔═╡ 2e90f5d0-e5ee-4abb-932b-5638afb3a146
md"""
**Test #2.**
Expected result = 230.
"""

# ╔═╡ b5db6a00-eb1d-4cf0-bfad-d04b0d7d9968
main2("test.in")

# ╔═╡ 4d9b89d5-7f30-4031-a619-57fdb5936dc5
md"""
**Solution.**
Confirmed result = 4481199.
"""

# ╔═╡ 17e05631-6ee0-49f0-b273-c151812220af
main2("data.in")

# ╔═╡ 00000000-0000-0000-0000-000000000001
PLUTO_PROJECT_TOML_CONTENTS = """
[deps]
"""

# ╔═╡ 00000000-0000-0000-0000-000000000002
PLUTO_MANIFEST_TOML_CONTENTS = """
# This file is machine-generated - editing it directly is not advised

julia_version = "1.7.0"
manifest_format = "2.0"

[deps]
"""

# ╔═╡ Cell order:
# ╟─d7b19f60-6348-11ec-11b8-abe5c0a6cee7
# ╟─14d4296b-4c17-47b5-80f2-1dcf89690f86
# ╠═2ea7d655-35ed-4687-9873-e658528232a8
# ╟─d9e09af7-42a2-4981-ab2f-21e7308c45e4
# ╠═8ceda75d-d741-4c1a-a884-dba74255270a
# ╟─aeba97ec-fd96-4ef4-99b0-726665015217
# ╠═7f03ad5e-bf73-4b16-958b-592ba4b6d923
# ╟─8cc55754-6815-4c08-8538-b6f4c74d5c28
# ╠═d0d5ecd6-fa56-4fca-a293-fb0a682e4a7b
# ╟─63a48804-1a55-48c6-832f-74f52ac1d583
# ╠═6c58e558-21fe-4221-b324-1efd90c7911a
# ╟─0e681b8a-82bc-425c-a374-e36d360559d2
# ╠═5012e220-d2e2-436c-9ce4-dfe34f365df5
# ╠═774f7d4b-9d47-41eb-a6d7-e9d0fc579fbb
# ╟─2e90f5d0-e5ee-4abb-932b-5638afb3a146
# ╠═b5db6a00-eb1d-4cf0-bfad-d04b0d7d9968
# ╟─4d9b89d5-7f30-4031-a619-57fdb5936dc5
# ╠═17e05631-6ee0-49f0-b273-c151812220af
# ╟─00000000-0000-0000-0000-000000000001
# ╟─00000000-0000-0000-0000-000000000002


text:
# You know it's getting real when this is uncommented
# __precompile__()

module BinaryBuilder
using Compat, Compat.Libdl, Compat.LibGit2, Compat.Random
using Reexport
using ObjectFile
using Nullables
using GitHub

@reexport using BinaryProvider

include("compat.jl")
include("Auditor.jl")
include("Runner.jl")
include("RootFS.jl")
include("UserNSRunner.jl")
include("QemuRunner.jl")
include("Dependency.jl")
include("AutoBuild.jl")
include("Wizard.jl")

# This is a global github authentication token that is set the first time
# we authenticate and then reused
const _github_auth = Ref{GitHub.Authorization}()

function github_auth()
    if !isassigned(_github_auth)
        # If the user is feeding us a GITHUB_TOKEN token, use it!
        if length(get(ENV, "GITHUB_TOKEN", "")) == 40
            _github_auth[] = GitHub.authenticate(ENV["GITHUB_TOKEN"])
        else
            _github_auth[] = GitHub.AnonymousAuth()
        end
    end
    return _github_auth[]
end

# These globals store important information such as where we're downloading
# the rootfs to, and where we're unpacking it.  These constants are initialized
# by `__init__()` to allow for environment variable overrides from the user.
downloads_cache = ""
rootfs_cache = ""
shards_cache = ""
qemu_cache = ""
automatic_apple = false
use_squashfs = false
allow_ecryptfs = false
use_ccache = false
ccache_override = ""
sandbox_override = ""

function __init__()
    global downloads_cache, rootfs_cache, shards_cache, runner_override
    global qemu_cache, use_squashfs, automatic_apple, allow_ecryptfs
    global use_ccache, ccache_override, sandbox_override

    # If the user has overridden our rootfs tar location, reflect that here:
    def_dl_cache = joinpath(dirname(@__FILE__), "..", "deps", "downloads")
    downloads_cache = get(ENV, "BINARYBUILDER_DOWNLOADS_CACHE", def_dl_cache)
    downloads_cache = abspath(downloads_cache)

    # If the user has overridden our rootfs unpack location, reflect that here:
    def_rootfs_cache = joinpath(dirname(@__FILE__),  "..", "deps", "root")
    rootfs_cache = get(ENV, "BINARYBUILDER_ROOTFS_DIR", def_rootfs_cache)
    rootfs_cache = abspath(rootfs_cache)

    # If the user has overridden our shards unpack location, reflect that here:
    def_shards_cache = joinpath(dirname(@__FILE__),  "..", "deps", "shards")
    shards_cache = get(ENV, "BINARYBUILDER_SHARDS_DIR", def_shards_cache)
    shards_cache = abspath(shards_cache)

    # If the user has overridden our qemu unpack location, reflect that here:
    def_qemu_cache = joinpath(dirname(@__FILE__), "..", "deps", "qemu")
    qemu_cache = get(ENV, "BINARYBUILDER_QEMU_DIR", def_qemu_cache)
    qemu_cache = abspath(qemu_cache)

    # If the user has asked for squashfs mounting instead of tarball mounting,
    # use that here.  Note that on Travis, we default to using squashfs, unless
    # BINARYBUILDER_USE_SQUASHFS is set to "false", which overrides this
    # default. If we are not on Travis, we default to using tarballs and not
    # squashfs images as using them requires `sudo` access.
    if get(ENV, "BINARYBUILDER_USE_SQUASHFS", "") == "false"
        use_squashfs = false
    elseif get(ENV, "BINARYBUILDER_USE_SQUASHFS", "") == "true"
        use_squashfs = true
    else
        # If it hasn't been specified, but we're on Travis, default to "on"
        if get(ENV, "TRAVIS", "") == "true"
            use_squashfs = true
        end
    end

    # If the user has signalled that they really want us to automatically
    # accept apple EULAs, do that.
    if get(ENV, "BINARYBUILDER_AUTOMATIC_APPLE", "") == "true"
        automatic_apple = true
    end

    # If the user has overridden our runner selection algorithms, honor that
    runner_override = lowercase(get(ENV, "BINARYBUILDER_RUNNER", ""))
    if !(runner_override in ["", "userns", "qemu", "privileged"])
        Compat.@warn("Invalid runner value $runner_override, ignoring...")
        runner_override = ""
    end

    # If the user has signified that they want to allow mounting of ecryptfs
    # paths, then let them do so at their own peril.
    if get(ENV, "BINARYBUILDER_ALLOW_ECRYPTFS", "") == "true"
        allow_ecryptfs = true
    end

    # If the user has enabled `ccache` support, use it!
    if get(ENV, "BINARYBUILDER_USE_CCACHE", "false") == "true"
        use_ccache = true
    end

    # If the user has asked to use a particular `sandbox` executable, use it!
    sandbox_override = get(ENV, "BINARYBUILDER_SANDBOX_PATH", "")
    if !isempty(sandbox_override) && !isfile(sandbox_override)
        Compat.@warn("Invalid sandbox path $sandbox_override, ignoring...")
        sandbox_override = ""
    end

    # If the user asked for the ccache cache to be placed somewhere specific
    ccache_override = get(ENV, "BINARYBUILDER_CCACHE_PATH", "")
    if !isempty(ccache_override) && !isdir(ccache_override)
        Compat.@warn("Invalid ccache directory $ccache_override, ignoring...")
        ccache_override = ""
    end
end

"""
    versioninfo()

Helper function to print out some debugging information
"""
function versioninfo()
    Compat.@info("Julia versioninfo(): ")
    Base.versioninfo()

    # Get BinaryBuilder.jl's git sha
    repo = LibGit2.GitRepo(Pkg.dir("BinaryBuilder"))
    gitsha = hex(LibGit2.GitHash(LibGit2.GitCommit(repo, "HEAD")))
    Compat.@info("BinaryBuilder.jl version: $(gitsha)")
    @static if Compat.Sys.isunix()
        Compat.@info("Kernel version: $(readchomp(`uname -r`))")
    end

    # Dump if some important directories are encrypted:
    @static if Compat.Sys.islinux()
        print_enc(n, path) = begin
            is_encrypted, mountpoint = is_ecryptfs(path)
            if is_encrypted
                Compat.@info("$n is encrypted on mountpoint $mountpoint")
            else
                Compat.@info("$n is NOT encrypted on mountpoint $mountpoint")
            end
        end

        print_enc("pkg dir", dirname(@__FILE__))
        print_enc("rootfs dir", rootfs_dir())
        print_enc("shards dir", shards_dir())
    end

    # Dump any relevant environment variables:
    Compat.@info("Relevant envionment variables:")
    env_var_suffixes = [
        "AUTOMATIC_APPLE",
        "USE_SQUASHFS",
        "DOWNLOADS_CACHE",
        "ROOTFS_DIR",
        "SHARDS_DIR",
        "QEMU_DIR",
        "RUNNER",
        "ALLOW_ECRYPTFS",
        "USE_CCACHE",
    ]
    for e in env_var_suffixes
        envvar = "BINARYBUILDER_$(e)"
        if haskey(ENV, envvar)
            Compat.@info("  $(envvar): \"$(ENV[envvar])\"")
        end
    end

    # Print out the preferred runner stuff here:
    Compat.@info("Preferred runner: $(preferred_runner())")

    # Try to run 'echo julia' in Linux x86_64 environment
    Compat.@info("Trying to run `echo hello julia` within a Linux x86_64 environment...")

    runner = preferred_runner()(
        pwd();
        cwd="/workspace/",
        platform=Linux(:x86_64),
        verbose=true
    )
    run_interactive(runner, `/bin/bash -c "echo hello julia"`)

    # If we use ccache, dump the ccache stats
    if use_ccache
        Compat.@info("ccache stats:")
        runner = preferred_runner()(
            pwd();
            cwd="/workspace/",
            platform=Linux(:x86_64),
        )
        run_interactive(runner, `/usr/local/bin/ccache -s`)
    end
end

end # module


text:
module MyPackage

export Shape, getPosition, setPosition!
export Box
export Circle, getSize, setSize!
export Composite, numberOfShapes, addShape!

greet() = print("Hello World!")

abstract type Shape end
function getPosition(s::Shape)
    return s.position
end

function setPosition!(s::Shape, newPosition::Tuple{Float64, Float64})
    return s.position = newPosition
end

mutable struct Box <: Shape
    position::Tuple{Float64, Float64}
    eventCallbacks::Vector
    height::Float64
    width::Float64
end
Box() = Box((0.0, 0.0), [], 5, 5)

mutable struct Circle <: Shape
    position::Tuple{Float64, Float64}
    eventCallbacks::Vector
    size::Float64
end
Circle() = Circle((0.0, 0.0), [], 5)
getSize(c::Circle) = c.size
setSize!(c::Circle, aSize::Int64) = setSize!(c::Circle, Float64(aSize))
function setSize!(c::Circle, aSize::Float64)
    c.size = aSize
end

mutable struct Composite <: Shape
    position::Tuple{Float64, Float64}
    eventCallbacks::Vector
    shapes::Vector{Shape}
end
Composite() = Composite((0.0, 0.0), [], [])
function addShape!(c::Composite, childShape::Shape)
    push!(c.shapes, childShape)
end
numberOfShapes(c::Composite) = length(c.shapes)
end # module