text:
# using Revise, Plots
using Test
using BifurcationKit, SparseArrays, Parameters, LinearAlgebra, Setfield, ForwardDiff
const BK = BifurcationKit
const FD = ForwardDiff

norminf = x -> norm(x, Inf)
function Laplacian1D(N, lx, bc = :none)
	hx = 2lx/N
	Δ = spdiagm(0 => -2ones(N), 1 => ones(N-1), -1 => ones(N-1) )
	Δ[1,end]=1; Δ[end,1]=1
	D = spdiagm(1 => ones(N-1), -1 => -ones(N-1) )
	D[1,end]=-1; D[end,1]=1
	D = D / (2hx)
	Δ = Δ / hx^2
	return Δ, D
end

# add to f the nonlinearity
@views function NL!(f, u, p)
	@unpack r, μ, ν, c3, c5, γ = p
	n = div(length(u), 2)
	u1 = u[1:n]
	u2 = u[n+1:2n]

	ua = u1.^2 .+ u2.^2

	f1 = f[1:n]
	f2 = f[n+1:2n]

	f1 .+= @. r * u1 - ν * u2 - ua * (c3 * u1 - μ * u2) - c5 * ua^2 * u1 + γ
	f2 .+= @. r * u2 + ν * u1 - ua * (c3 * u2 + μ * u1) - c5 * ua^2 * u2

	return f
end

function Fcgl!(f, u, p, t = 0)
	mul!(f, p.Δ, u)
	NL!(f, u, p)
end

Fcgl(u, p, t = 0) = Fcgl!(similar(u), u, p)

# remark: I checked this against finite differences
@views function Jcgl(u, p)
	@unpack r, μ, ν, c3, c5, Δ = p

	n = div(length(u), 2)
	u1 = u[1:n]
	u2 = u[n+1:2n]

	ua = u1.^2 .+ u2.^2

	f1u = zero(u1)
	f2u = zero(u1)
	f1v = zero(u1)
	f2v = zero(u1)

	@. f1u =  r - 2 * u1 * (c3 * u1 - μ * u2) - c3 * ua - 4 * c5 * ua * u1^2 - c5 * ua^2
	@. f1v = -ν - 2 * u2 * (c3 * u1 - μ * u2)  + μ * ua - 4 * c5 * ua * u1 * u2
	@. f2u =  ν - 2 * u1 * (c3 * u2 + μ * u1)  - μ * ua - 4 * c5 * ua * u1 * u2
	@. f2v =  r - 2 * u2 * (c3 * u2 + μ * u1) - c3 * ua - 4 * c5 * ua * u2 ^2 - c5 * ua^2

	jacdiag = vcat(f1u, f2v)

	Δ + spdiagm(0 => jacdiag, n => f1v, -n => f2u)
end

jet = BK.getJet(Fcgl, Jcgl)
####################################################################################################
n = 50
	l = pi

	Δ, D = Laplacian1D(n, l, :Periodic)
	par_cgl = (r = 0.0, μ = 0.5, ν = 1.0, c3 = -1.0, c5 = 1.0, Δ = blockdiag(Δ, Δ), Db = blockdiag(D, D), γ = 0.0, δ = 1.0, N = 2n)
	sol0 = zeros(par_cgl.N)

_sol0 = zeros(2n)
_J0 = Jcgl(_sol0, par_cgl)
_J1 = FD.jacobian(z->Fcgl(z, par_cgl), _sol0) |> sparse
@test _J0 ≈ _J1

eigls = EigArpack(1.0, :LM)
	eigls = DefaultEig()
	# eigls = eig_MF_KrylovKit(tol = 1e-8, dim = 60, x₀ = rand(ComplexF64, Nx*Ny), verbose = 1)
	opt_newton = NewtonPar(tol = 1e-9, verbose = true, eigsolver = eigls, maxIter = 20)
	out, hist, flag = @time newton(jet[1], jet[2], sol0, par_cgl, opt_newton, normN = norminf)

opts_br = ContinuationPar(dsmin = 0.001, dsmax = 0.15, ds = 0.001, pMax = 2.5, detectBifurcation = 3, nev = 9, plotEveryStep = 50, newtonOptions = (@set opt_newton.verbose = false), maxSteps = 1060, nInversion = 8, maxBisectionSteps=20)
br, = continuation(jet[1], jet[2], vec(sol0), par_cgl, (@lens _.r), opts_br, verbosity = 0)

####################################################################################################
# we test the jacobian
_J0 = jet[2](sol0, par_cgl)
_J1 = FD.jacobian(z->jet[1](z, par_cgl), sol0) |> sparse
@test _J0 == _J1
####################################################################################################
function guessFromHopfO2(branch, ind_hopf, eigsolver, M, z1, z2 = 0.; phase = 0, k = 1.)
	specialpoint = branch.specialpoint[ind_hopf]
	@show specialpoint.ind_ev

	# parameter value at the Hopf point
	p_hopf = specialpoint.param

	# frequency at the Hopf point
	# @show branch.eig[specialpoint.idx].eigenvals[specialpoint.ind_ev:specialpoint.ind_ev+4]
	ωH = imag(branch.eig[specialpoint.idx].eigenvals[specialpoint.ind_ev]) |> abs

	# vec_hopf is the eigenvector for the eigenvalues iω
	vec_hopf1 = geteigenvector(eigsolver, br.eig[specialpoint.idx][2], specialpoint.ind_ev)
	vec_hopf1 ./=  norm(vec_hopf1)

	vec_hopf2 = geteigenvector(eigsolver, br.eig[specialpoint.idx][2], specialpoint.ind_ev - 2)
	vec_hopf2 ./=  norm(vec_hopf2)

	 orbitguess = [real.(specialpoint.x .+
	 			z1 .* vec_hopf1 .* exp(2pi * complex(0, 1) .* (ii/(M-1) - phase)) .+
				z2 .* vec_hopf2 .* exp(2pi * complex(0, 1) .* (ii/(M-1) - phase))) for ii=0:M-1]

	 X = LinRange(-pi, pi, n) |> collect; X = vcat(X, X)
	 # orbitguess = [real.(specialpoint.x .+ z1 .* exp.(complex(0, 1) .* (2pi * ii/(M-1) .- k .* X))) for ii=0:M-1]

	return p_hopf, 2pi/ωH, orbitguess, specialpoint.x, vec_hopf1, vec_hopf2
end
####################################################################################################
# we test TWProblem: travelling wave problem
# number of time slices in the periodic orbit
M = 50

# TW ansatz
r_hopf, Th, orbitguess2, hopfpt, eigvec = guessFromHopfO2(br, 2, opt_newton.eigsolver, M, 1. + 0.0im, 1+0.0im; k = 2.) #TW

uold = copy(orbitguess2[1][1:2n])
# plot(uold[1:end-1]; linewidth = 5)

# we create a TW problem
probTW = BK.TWProblem(Fcgl, Jcgl, par_cgl.Db, copy(uold))
probTW(vcat(uold,.1), par_cgl)

# we test the sparse formulation of the problem jacobian
_sol0 = rand(2n+1)
_J1 = FD.jacobian(z->probTW(z, par_cgl), _sol0) |> sparse
_J0 = probTW(Val(:JacFullSparse), _sol0, par_cgl)
@test _J1 ≈ _J0

# we test the matrix-free formulation of the problem jacobian
_sol0 = rand(2n+1)
_dsol0 = rand(2n+1)
_out1 = FD.derivative(t -> probTW(_sol0 .+ t .* _dsol0, par_cgl), 0)
_out0 = probTW(_sol0, par_cgl, _dsol0)
@test _out0 ≈ _out1

BK.VFtw(probTW, uold, (user=par_cgl, s=Tuple(0.,)))

# we test the ∂
BK.applyD(probTW, rand(2n))

# we test update section
BK.updateSection!(probTW, probTW.u₀)
####################################################################################################
# test newton method, not meant to converge
newton(probTW, vcat(uold,.1), par_cgl, NewtonPar(verbose = true, maxIter = 5), jacobian = :AutoDiff)
newton(probTW, vcat(uold,.1), par_cgl, NewtonPar(verbose = true, maxIter = 5), jacobian = :FullLU)
####################################################################################################
# test continuation method witth different Generalised eigensolvers
optn = NewtonPar(tol = 1e-8)
opt_cont_br = ContinuationPar(pMin = -1., pMax = 1., newtonOptions = optn, maxSteps = 3, detectBifurcation = 2)
continuation(probTW, vcat(uold,.1), par_cgl, (@lens _.r), opt_cont_br; jacobian = :FullLU, verbosity = 0)

@set! opt_cont_br.newtonOptions.eigsolver = BK.DefaultGEig(B = diagm(0=>vcat(ones(2n),0)))
continuation(probTW, vcat(uold,.1), par_cgl, (@lens _.r), opt_cont_br; jacobian = :FullLU, verbosity = 0)

BK.GEigArpack(nothing, :LR)
@set! opt_cont_br.newtonOptions.eigsolver = EigArpack(nev = 5, which = :LM, sigma = 0.2, v0 = rand(2n+1))
continuation(probTW, vcat(uold,.1), par_cgl, (@lens _.r), opt_cont_br; jacobian = :AutoDiff, verbosity = 0)

@set! opt_cont_br.newtonOptions.linsolver = GMRESIterativeSolvers(N = 2n+1)
@set! opt_cont_br.newtonOptions.eigsolver = EigArpack(nev = 4, ncv = 2n+1, tol = 1e-3, v0 = rand(2n+1))
@set! opt_cont_br.detectBifurcation = 0
continuation(probTW, vcat(uold,.1), par_cgl, (@lens _.r), opt_cont_br; jacobian = :MatrixFreeAD, verbosity = 0)


text:
@testset "Table" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 5, Int)
    @variable(m, -5 <= y <= 5, Int)
    @variable(m, -5 <= z <= 5, Int)

    tab = [
        1 2 3
        1 3 2
        2 1 3
        2 3 1
        3 1 2
        3 2 1
        4 4 4
    ]

    @constraint(m, [x, y, z] in CS.TableSet(tab))
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    # check if impossible values got removed
    constr_indices = constraint.indices
    for ind in constr_indices
        @test sort(CS.values(com.search_space[ind])) == 1:4
    end

    @test !CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 4])
    @test CS.is_constraint_solved(constraint, constraint.fct, constraint.set, [1, 2, 3])

    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        4,
    )
    @test CS.fix!(com, com.search_space[constr_indices[2]], 2)
    @test !CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        4,
    )

    # need to create a backtrack_vec to reverse pruning
    dummy_backtrack_obj = CS.BacktrackObj(com)
    dummy_backtrack_obj.step_nr = 1
    push!(com.backtrack_vec, dummy_backtrack_obj)
    # reverse previous fix
    CS.reverse_pruning!(com, 1)
    com.c_backtrack_idx = 1

    @test CS.still_feasible(
        com,
        constraint,
        constraint.fct,
        constraint.set,
        constr_indices[1],
        4,
    )

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 5, Int)
    @variable(m, -5 <= y <= 5, Int)
    @variable(m, -5 <= z <= 5, Int)

    tab = [
        1 2 3
        1 3 2
        2 1 3
        2 3 1
        3 1 2
        3 2 1
        4 4 4
    ]

    @constraint(m, [x, y, z] in CS.TableSet(tab))
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    constr_indices = constraint.indices

    @test CS.fix!(com, com.search_space[constr_indices[2]], 2)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices[[1, 3]]
        @test sort(CS.values(com.search_space[ind])) == [1, 3]
    end

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, -5 <= x <= 5, Int)
    @variable(m, -5 <= y <= 5, Int)
    @variable(m, -5 <= z <= 5, Int)

    tab = [
        1 2 3
        1 3 2
        2 1 3
        2 3 1
        3 1 2
        3 2 1
        4 4 4
    ]

    @constraint(m, [x, y, z] in CS.TableSet(tab))
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]
    constr_indices = constraint.indices

    @test CS.fix!(com, com.search_space[constr_indices[2]], 4)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    for ind in constr_indices
        @test CS.isfixed(com.search_space[ind])
        @test CS.value(com.search_space[ind]) == 4
    end
end

@testset "table is_constraint_violated test" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x <= 5, Int)
    @variable(m, 0 <= y <= 3, Int)
    @constraint(m, [x, y] in CS.TableSet([
        1 2
        3 4
        2 1
        2 3
    ]))
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 3; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x <= 5, Int)
    @variable(m, 0 <= y <= 3, Int)
    @constraint(m, [x, y] in CS.TableSet([
        1 2
        3 4
        2 1
        2 3
    ]))
    optimize!(m)
    com = CS.get_inner_model(m)

    constraint = com.constraints[1]

    variables = com.search_space
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
    @test !CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end

@testset "all modulo" begin
    function modulo(m, x, y, z)
        lbx = !(x isa Integer)  ? round(Int, JuMP.lower_bound(x)) : x
        ubx = !(x isa Integer) ? round(Int, JuMP.upper_bound(x)) : x
        lby = !(y isa Integer) ? round(Int, JuMP.lower_bound(y)) : y
        uby = !(y isa Integer) ? round(Int, JuMP.upper_bound(y)) : y

        table =  transpose(reduce(hcat,[ [i,j,i % j] for i in lbx:ubx, j in lby:uby if j != 0]))
        @constraint(m, [x, y, z] in CS.TableSet(table))
    end

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_solutions" => true))
    @variable(m, 1 <= x[1:2] <= 10, Int)
    modulo(m,x[1],2,x[2])
    optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xv = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test xv[1] % 2 == xv[2]
        push!(results, (xv[1], xv[2]))
    end

    found_nr = 0
    for i in 1:10, j in 1:10
        if i % 2 == j
            @test (i,j) in results
            found_nr += 1
        end
    end
    @test found_nr == nresults
end

@testset "vector affine" begin
    function modulo_affine(m, x, y, z)
        lbx = !(x isa Integer)  ? round(Int, JuMP.lower_bound(x)) : x
        ubx = !(x isa Integer) ? round(Int, JuMP.upper_bound(x)) : x
        lby = !(y isa Integer) ? round(Int, JuMP.lower_bound(y)) : y
        uby = !(y isa Integer) ? round(Int, JuMP.upper_bound(y)) : y

        table =  transpose(reduce(hcat,[ [i,j,i % j] for i in lbx:2*ubx, j in lby:uby if j != 0]))
        @constraint(m, [x+z, y, 2x-z] in CS.TableSet(table))
    end

    m = Model(optimizer_with_attributes(CS.Optimizer, "logging" => [], "all_solutions" => true))
    @variable(m, 1 <= x[1:2] <= 10, Int)
    modulo_affine(m,x[1],2,x[2])
    optimize!(m)

    @test JuMP.termination_status(m) == MOI.OPTIMAL
    nresults = JuMP.result_count(m)
    results = Set{Tuple{Int,Int}}()
    for i in 1:nresults
        xv = convert.(Int, round.(JuMP.value.(x; result=i)))
        @test (xv[1]+xv[2]) % 2 == 2*xv[1]-xv[2]
        push!(results, (xv[1], xv[2]))
    end

    found_nr = 0
    for i in 1:10, j in 1:10
        if (i+j) % 2 == 2*i-j
            @test (i,j) in results
            found_nr += 1
        end
    end
    @test found_nr == nresults
end

ext: jl


text:
module Day17

using ..ReTest, ..DocStringExtensions

import ..solve, ..PuzzleInfo, ..TEST_DATA_DIRECTORY

abstract type Day17PuzzleInfo <: PuzzleInfo end

"""
Day 17 part 1 puzzle information.

$(FIELDS)
"""
struct Part1PuzzleInfo <: Day17PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
Day 17 part 2 puzzle information.

$(FIELDS)
"""
struct Part2PuzzleInfo <: Day17PuzzleInfo
    """Path to input data file."""
    datapath::String
end

"""
    $(FUNCTIONNAME)(input)

Read target area limits from input stream `input`.
"""
function readtargetarealimits(input::IO)
    targetareastring = readline(input)
    @assert targetareastring[1:13] == "target area: "
    xlimitstring, ylimitstring = split(targetareastring[14:end], ", ")
    @assert xlimitstring[1:2] == "x="
    @assert ylimitstring[1:2] == "y="
    xlimits = [parse(Int, lim) for lim in split(xlimitstring[3:end], "..")]
    ylimits = [parse(Int, lim) for lim in split(ylimitstring[3:end], "..")]
    xlimits, ylimits
end

"""
    $(FUNCTIONNAME)(xlimits, ylimits, initialvelocity)

Whether probe launched with `initialvelocity` will enter target area with horizontal
limits `xlimits` and vertical limits `ylimits` in any step.
"""
function hitstargetarea(xlimits, ylimits, initialvelocity)
    pos = [0, 0]
    velocity = copy(initialvelocity)
    while pos[1] <= xlimits[2] && pos[2] >= ylimits[1]
        if xlimits[1] <= pos[1] <= xlimits[2] && ylimits[1] <= pos[2] <= ylimits[2]
            return true
        end
        pos += velocity
        velocity[1] = max(0, velocity[1] - 1)
        velocity[2] -= 1
    end
    false
end

function solve(info::Day17PuzzleInfo)
    xlimits, ylimits = open(readtargetarealimits, info.datapath)
    ypos(velocity, nstep) = nstep * velocity - (nstep * (nstep - 1) ÷ 2)
    if isa(info, Part1PuzzleInfo)
        if ylimits[1] > 0
            ypos(ylimits[1], ylimit[1])d
        elseif ylimits[1] < 0
            ypos(-ylimits[1] - 1, -ylimits[1] - 1)
        else
            error("Target area lower boundary allows infinite velocity")
        end
    else
        numvel = 0
        # Lower bound on x-velocity corresponds to just hitting left limit of target
        # area as x-velocity reaches 0: vx * (vx + 1) ÷ 2 >= xlimits[1] for integer vx
        # Upper bound on x-velocity correspondings to hitting right limit of target area
        # in one step
        for vx in ceil(Int, (sqrt(1 + 8 * xlimits[1]) - 1) / 2):xlimits[2]
            # Lower bound on y-velocity corresponds to hitting lower limit in one step
            # Upper bound on y-velocity corresponds to hitting lower limit in
            # (2 * vy + 1) steps
            for vy in ylimits[1]:(-ylimits[1] - 1)
                hitstargetarea(xlimits, ylimits, [vx, vy]) && (numvel += 1)
            end
        end
        numvel
    end
end

@testset begin
    @test solve(Part1PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day17.txt"))) == 45
    @test solve(Part2PuzzleInfo(joinpath(TEST_DATA_DIRECTORY, "day17.txt"))) == 112
end

end # module


text:
"""
    StateVecFuncs

  Provides a set of helper functions that operate on StateVec.
"""
module StateVecFuncs

using Requires
@init @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
  using .Plots
  export plot_state
end

using DelimitedFiles, WriteVTK, ..Grids, ..StateVecs
export surface_val, first_elem_above_surface_val
export domain_average!, distribute!, total_covariance!
export extrap!, assign_ghost!, integrate_ode!
export export_state, UseVTK, UseDat

"""
    surface_val(sv::StateVec, name::Symbol, grid::Grid, i_sd=1)

Get the value of variable `name` on the surface (`z_min`) by
interpolating between the first ghost point and the first element
above the surface
"""
function surface_val(sv::StateVec, name::Symbol, grid::Grid, i_sd=1)
  return (sv[name, grid.n_ghost, i_sd]+sv[name, 1+grid.n_ghost, i_sd])/2
end

"""
    first_elem_above_surface_val(sv::StateVec, name::Symbol, grid::Grid, i_sd=1)

Get the value of variable `name` on the first element above the
surface (`z_min+Δz/2`).
"""
function first_elem_above_surface_val(sv::StateVec, name::Symbol, grid::Grid, i_sd=1)
  return sv[name, 1+grid.n_ghost, i_sd]
end

"""
    domain_average!(dst::StateVec, src::StateVec, weight::StateVec
                    dst_idxs, src_idxs, weight_idx, grid::Grid)

Compute the domain average in state vector `dst`, given state vectors `src`
and `weight`, the grid `grid` and index iterators of source, destination and
weight names `src_idxs`, `dst_idxs`, and `weight_idx` respectively.

Formulaically, a domain-averaged variable ``⟨ϕ⟩`` is computed from

``⟨ϕ⟩ = Σ_i a_i \\overline{ϕ}_i``

Where variable ``\\overline{ϕ}_i`` represents ``ϕ`` decomposed across multiple
sub-domains, which are weighted by area fractions ``a_i``.

Note that `domain_average!` is the inverse function of `distribute!`.
"""
function domain_average!(dst::StateVec, src::StateVec, weight::StateVec,
                         dst_idxs, src_idxs, weight_idx, grid::Grid)
  for k in over_elems(grid)
    for (dst_idx, src_idx) in zip(dst_idxs, src_idxs)
      dst[dst_idx, k] = 0
      for i in over_sub_domains(src)
        dst[dst_idx, k] += src[src_idx, k, i]*weight[weight_idx, k, i]
      end
    end
  end
end

"""
    distribute!(dst::StateVec, src::StateVec, dst_idxs, src_idxs, grid::Grid)

Distributes values in the state vector `src`, to state vectors `dst` given
the grid `grid` and index iterators of source and destination names `src_idxs`
and `dst_idxs` respectively.

Formulaically, a domain-decomposed variable ``\\overline{ϕ}_i`` is computed from

``\\overline{ϕ}_i = ⟨ϕ⟩``

Where variable ``⟨ϕ⟩`` is the domain-averaged variable, computed
across multiple sub-domains.

Note that `distribute!` is the inverse function of `domain_average!`.
"""
function distribute!(dst::StateVec, src::StateVec, dst_idxs, src_idxs, grid::Grid)
  for k in over_elems(grid), i in over_sub_domains(dst)
    for (dst_idx, src_idx) in zip(dst_idxs, src_idxs)
      dst[dst_idx, k, i] = src[src_idx, k]
    end
  end
end

"""
    total_covariance!(dst::StateVec, src::StateVec, cv::StateVec, weights::StateVec,
                      dst_idxs, src_idxs, cv_idxs, weight_idx,
                      grid::Grid, decompose_ϕ_ψ::Function)

Computes the total covariance in state vector `dst`, given
 - `src` source state vector
 - `cv` state vector containing co-variances
 - `weights` state vector containing weights
 - `dst_idxs` indexes for destination state vector
 - `cv_idxs` indexes for state vector containing co-variances
 - `weight_idx` indexes for state vector containing weights
 - `grid` the grid
 - `decompose_ϕ_ψ` a function that receives the covariance index and
                   returns the indexes for each variable. For example:
                   `:ϕ_idx, :ψ_idx = decompose_ϕ_ψ(:cv_ϕ_ψ)`

Formulaically, a total covariance between variables ``ϕ`` and ``ψ`` is computed from

``⟨ϕ^*ψ^*⟩ = Σ_i a_i \\overline{ϕ_i'ψ_i'} + Σ_i Σ_j a_i a_j \\overline{ϕ}_i (\\overline{ψ}_i - \\overline{ψ}_j)``

Where variable ``\\overline{ϕ}_i`` represents ``ϕ`` decomposed across multiple
sub-domains, which are weighted by area fractions ``a_i``.
"""
function total_covariance!(dst::StateVec, src::StateVec, cv::StateVec, weights::StateVec,
                           dst_idxs, cv_idxs, weight_idx, grid::Grid, decompose_ϕ_ψ::Function)
  for k in over_elems(grid), (dst_idx, cv_idx) in zip(dst_idxs, cv_idxs)
    _ϕ, _ψ = decompose_ϕ_ψ(cv_idx)
    dst[dst_idx, k] = 0
    for i in over_sub_domains(weights)
      dst[dst_idx, k] += weights[weight_idx, k, i]*cv[cv_idx, k, i]
      for j in over_sub_domains(src)
        dst[dst_idx, k] += weights[weight_idx, k, i]*
                           weights[weight_idx, k, j]*
                           src[_ϕ, k, i]*(src[_ψ, k, i] - src[_ψ, k, j])
      end
    end
  end
end

"""
    extrap!(sv::StateVec, name::Symbol, grid::Grid, i_sd=1)

Extrapolate variable `name` to the first ghost point.
"""
function extrap!(sv::StateVec, name::Symbol, grid::Grid, i_sd=1)
  k = 1+grid.n_ghost
  sv[name, k-1, i_sd] = 2*sv[name, k, i_sd] - sv[name, k+1, i_sd]
  k = grid.n_elem - grid.n_ghost
  sv[name, k+1, i_sd] = 2*sv[name, k, i_sd] - sv[name, k-1, i_sd]
end

"""
    assign_ghost!(sv::StateVec, name::Symbol, val, grid::Grid, i_sd=1)

Assign value `val` to variable `name` for all ghost points.
"""
function assign_ghost!(sv::StateVec, name::Symbol, val, grid::Grid, i_sd=1)
  for k in over_elems_ghost(grid)
    sv[name, k, i_sd] = val
  end
end

"""
    integrate_ode!(sv::StateVec,
                   name::Symbol,
                   grid::Grid,
                   func::Function,
                   y_0::AbstractFloat,
                   args::NamedTuple,
                   i_sd::Int = 1)

Integrates the Ordinary Differential Equation

  ``\\frac{dy}{dz} = f(z)``

using Newton method, given the arguments

 - `sv` a state vector
 - `name` the name of the variable integrated
 - `grid` the grid that the variable lives on
 - `func` the function ``f(z)`` to be integrated, which accepts arguments `func(z, args)`
 - `y_0` the boundary condition ``y|_{z=0} = y_0``
 - `args` a `NamedTuple` of arguments passed to `func`
and optionally,
 - `i_sd` the i-th sub-domain of the variable (default is 1)
"""
function integrate_ode!(sv::StateVec,
                        name::Symbol,
                        grid::Grid,
                        func::Function,
                        y_0::AbstractFloat,
                        args::NamedTuple,
                        i_sd::Int = 1)
  k = first_elem_above_surface(grid)
  sv[name, k-1, i_sd] = 0
  for k in over_elems_real(grid)
    f = func(get_z(grid, k), args)
    sv[name, k, i_sd] = sv[name, k-1, i_sd] + grid.dz * f
  end
  for k in over_elems_real(grid)
    sv[name, k, i_sd] += y_0
  end
  assign_ghost!(sv, name, 0.0, grid, i_sd)
end

abstract type ExportType end
"""
    UseVTK

A singleton used to indicate to use a `.vtk`
file extension when exporting a `StateVec`.
"""
struct UseVTK <: ExportType end

"""
    UseDat

A singleton used to indicate to use a `.dat` (DAT)
file extension when exporting a `StateVec`. A DAT
file is a generic data file which may contain data
in binary or text format.
"""
struct UseDat <: ExportType end

"""
    export_state(sv::StateVec, grid::Grid, dir, filename, ::ExportType)

Export StateVec to a human-readable file `filename` in directory `dir`.
"""
function export_state(sv::StateVec, grid::Grid, dir, filename, ::ExportType = UseDat())
  domain = over_elems(grid)
  headers = [ length(over_sub_domains(sv, name))==1 ? string(name) : string(name)*"_"*string(i_sd)
             for name in sv.var_names for i_sd in over_sub_domains(sv, name)]
  n_vars = length(headers)
  n_elem = length(domain)
  data = reshape([sv[name, k, i_sd] for name in sv.var_names for
    i_sd in over_sub_domains(sv, name) for k in domain], n_elem, n_vars)
  z = grid.z[domain]
  data_all = hcat(z, data)
  open(string(dir, filename*"_vs_z.dat"),"w") do file
    write(file, join(headers, ", ")*"\n")
    writedlm(file, data_all)
  end
end

"""
    export_state(sv::StateVec, dir, filename, ::ExportType = UseDat())

Export StateVec to a human-readable file `filename` in directory `dir`.
The `z`-axis is _not_ included in the export.
"""
function export_state(sv::StateVec, dir, filename, ::ExportType = UseDat())
  open(string(dir, filename*".dat"),"w") do file
    print(file, sv)
  end
end

"""
    export_state(sv::StateVec, grid::Grid, dir, filename, ::UseVTK)

Export state vector `sv` to a compressed file `filename` in directory `dir`,
including the `z`-axis, given by the grid `grid`.
"""
function export_state(sv::StateVec, grid::Grid, dir, filename, ::UseVTK)
  domain = over_elems(grid)
  headers = [ length(over_sub_domains(sv, name))==1 ? string(name) : string(name)*"_"*string(i_sd)
             for name in sv.var_names for i_sd in over_sub_domains(sv, name)]
  n_vars = length(headers)
  n_elem = length(domain)
  data = reshape([sv[name, k, i_sd] for name in sv.var_names for
    i_sd in over_sub_domains(sv, name) for k in domain], n_elem, n_vars)
  z = grid.z[domain]
  fields = ((headers[i], data[:,i]) for i in 1:n_vars)
  points = [1, length(domain)]
  cells = Array{MeshCell{Array{Int,1}}, 1}(undef, n_elem)
  for k in domain
    cells[k] = MeshCell(VTKCellTypes.VTK_LINE, [k, k+1])
  end
  vtkfile = vtk_grid("$(filename)", z, cells; compress=false)
  for (name, v) in fields
    vtk_point_data(vtkfile, v, name)
  end
  outfiles = vtk_save(vtkfile)
end

@init @require Plots = "91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin

"""
    plot_state(sv::StateVec,
                    grid::Grid,
                    name_idx::Symbol = nothing,
                    directory::AbstractString,
                    filename::AbstractString,
                    i_sd = 1,
                    include_ghost = false,
                    xlims::Union{Nothing, Tuple{R, R}} = nothing,
                    ylims::Union{Nothing, Tuple{R, R}} = nothing
                    ) where R

Save the plot variable along the z-direction in `StateVec` given the
grid, `grid`, variable name `name_idx`, directory `directory`,
filename `filename`, sub-domain `i_sd`, and a `Bool`, `include_ghost`,
indicating to include include or exclude the ghost points.
"""
function plot_state(sv::StateVec,
                    grid::Grid,
                    directory::AbstractString,
                    filename::AbstractString,
                    name_idx::Symbol = nothing,
                    i_sd = 1,
                    include_ghost = false,
                    xlims::Union{Nothing, Tuple{R, R}} = nothing,
                    ylims::Union{Nothing, Tuple{R, R}} = nothing
                    ) where R
  if name_idx == nothing
    domain_range = include_ghost ? over_elems(grid) : over_elems_real(grid)
    for name_idx in sv.var_names
      x = [grid.z[k] for k in domain_range]
      y = [sv[name_idx, k, i_sd] for k in domain_range]
      plot(y, x)
    end
    plot!(title = "state vector vs z", xlabel = "state vector", ylabel = "z")
    if xlims != nothing; plot!(xlims = xlims); end
    if ylims != nothing; plot!(ylims = xlims); end
  else
    x_name = string(name_idx)
    domain_range = include_ghost ? over_elems(grid) : over_elems_real(grid)
    x = [grid.z[k] for k in domain_range]
    y = [sv[name_idx, k, i_sd] for k in domain_range]
    plot(y, x)
    plot!(title = x_name * " vs z", xlabel = x_name, ylabel = "z")
    if xlims != nothing; plot!(xlims = xlims); end
    if ylims != nothing; plot!(ylims = xlims); end
  end
  savefig(joinpath(directory, filename))
end

end

end

max_stars_repo_name: christophernhill-dev/CLIMA

max_issues_repo_name: christophernhill-dev/CLIMA

max_forks_repo_name: christophernhill-dev/CLIMA


text:
using DiffEqBase, Catalyst, Plots, OrdinaryDiffEq, DataFrames, CSVFiles, LinearAlgebra
using ReactionNetworkImporters

# parameters
doplot = false
networkname = "Repressilator"
tf = 40000
nsteps = 400

# BNG simulation data
datadir  = joinpath(@__DIR__,"../data/repressilator")
fname    = joinpath(datadir, "Repressilator.net")
gdatfile = joinpath(datadir, "Repressilator.gdat")
print("getting gdat file...")
gdatdf = DataFrame(load(File{format"CSV"}(gdatfile), header_exists=true, spacedelim=true))
println("done")

# load the BNG reaction network in DiffEqBio
prnbng = loadrxnetwork(BNGNetwork(), fname)
rn = prnbng.rn; u0 = Float64[]; p = Float64[]
boprob = ODEProblem(rn, u0, (0.,tf), p)

# BNG simulation data testing
bngsol = gdatdf[!,:pTetR]

# note solvers run _much_ faster the second time
bsol = solve(boprob, Tsit5(), abstol=1e-12, reltol=1e-12, saveat=tf/nsteps);
@unpack pTetR = rn

if doplot
    plotlyjs()
    p1 = plot(gdatdf[!,:time], gdatdf[!,:pTetR], label=:BNGPTetR)
    plot!(p1, bsol.t, bsol[pTetR], label=:DEBPTetR)
    display(p1)
end

@test all(bsol.t .== gdatdf[!,:time])
@test all( abs.(gdatdf[!,:pTetR] - bsol[pTetR]) .< 1e-6 .* abs.(gdatdf[!,:pTetR]))


text:
# Mu' = f
type SteadyStateProblem{uType,isinplace,F,MM} <: AbstractSteadyStateProblem{uType,isinplace}
  f::F
  u0::uType
  mass_matrix::MM
end

function SteadyStateProblem(f,u0; iip = isinplace(f,3),mass_matrix=I)
  SteadyStateProblem{typeof(u0),iip,typeof(f),typeof(mass_matrix)}(f,u0,mass_matrix)
end

SteadyStateProblem(prob::AbstractODEProblem) =
      SteadyStateProblem(prob.f,prob.u0,iip=isinplace(prob),mass_matrix=prob.mass_matrix)


text:
module SymmetricDominance

export ModelParameters, simulate, tmrca, nsegsites, distances, spectra

typealias WithoutNewAllele Val{0}
typealias WithNewAllele Val{1}

type GeneDB
    currentid::Int
    currentstate::Int
    id::Vector{Int}
    epoch::Vector{Int}
    parent::Vector{Int}
    state::Vector{Int}

    GeneDB() = new(0, 0, Int[], Int[], Int[], Int[])
end

function insert!(::Type{WithoutNewAllele}, gdb::GeneDB, epoch::Int, parentid::Int)
    gdb.currentid += 1
    push!(gdb.id, gdb.currentid)
    push!(gdb.epoch, epoch)
    push!(gdb.parent, parentid)
    push!(gdb.state, gdb.state[parentid])
    gdb.currentid
end

function insert!(::Type{WithoutNewAllele}, gdb::GeneDB, epoch::Int, parentid::Int, state::Int)
    gdb.currentid += 1
    push!(gdb.id, gdb.currentid)
    push!(gdb.epoch, epoch)
    push!(gdb.parent, parentid)
    push!(gdb.state, state)
    state > gdb.currentstate && (gdb.currentstate = state)
    gdb.currentid
end

function insert!(::Type{WithNewAllele}, gdb::GeneDB, epoch::Int, parentid::Int)
    gdb.currentid += 1
    push!(gdb.id, gdb.currentid)
    push!(gdb.epoch, epoch)
    push!(gdb.parent, parentid)
    gdb.currentstate += 1
    push!(gdb.state, gdb.currentstate)
    gdb.currentid
end

zerof(x) = 0
function select(gdb::GeneDB, geneid::Int, fields::Symbol...)
    if geneid == 0
        ntuple(length(fields), zerof)
    else
        tuple([getfield(gdb, field)[geneid] for field in fields]...)
    end
end

immutable ModelParameters
    popsize::Int
    numberofloci::Int
    heterozygousfitness::Float64
    homozygousfitness::Float64
    recombinationrates::Vector{Float64}
    mutationrates::Vector{Float64}
end

immutable Organism
    geneids::Matrix{Int}
    Organism(nloci::Int) = new(Array(Int, nloci, 2))
end

immutable Population
    data::Vector{Organism}
    Population(n::Int, nloci::Int) = new([Organism(nloci) for _ = 1:n])
end
Base.start(pop::Population) = start(pop.data)
Base.next(pop::Population, iter) = next(pop.data, iter)
Base.done(pop::Population, iter) = done(pop.data, iter)
Base.length(pop::Population) = length(pop.data)

nloci(pop::Population) = size(pop[1].geneids, 1)
nloci(org::Organism) = size(org.geneids, 1)

Base.getindex(org::Organism, locus::Int, chr::Int) = org.geneids[locus, chr]
Base.setindex!(org::Organism, val::Int, locus::Int, chr::Int) = org.geneids[locus, chr] = val
Base.getindex(pop::Population, orgid::Int) = pop.data[orgid]
Base.getindex(pop::Population, orgid::Int, locus::Int, chr::Int) = pop.data[orgid].geneids[locus, chr]
Base.setindex!(pop::Population, val::Int, orgid::Int, locus::Int, chr::Int) = pop[orgid][locus, chr] = val

isidbystate(gdb::GeneDB, g1::Int, g2::Int) = select(gdb, g1, :state) == select(gdb, g2, :state)

mutate!(gdb::GeneDB, epoch::Int, geneid::Int) = insert!(WithNewAllele, gdb, epoch, geneid)

function getparentids!(ps, n)
    rand!(ps, 1:n)
    while ps[1] == ps[2]
        ps[2] = rand(1:n)
    end
    nothing
end

function hascoalesced(gdb::GeneDB, pop::Population, locus::Int, lastcoal::Int)
    coalesced = true
    ancestor = getancestor(gdb, pop[1, locus, 1], lastcoal)
    for org in pop , chr = 1:2
        if ancestor != getancestor(gdb, org[locus, chr], lastcoal)
            coalesced = false
            break
        end
    end
    coalesced
end

function getancestor(gdb::GeneDB, geneid::Int, lastcoal::Int)
    epoch = select(gdb, geneid, :epoch)[1]
    while epoch > lastcoal
        geneid = select(gdb, geneid, :parent)[1]
        epoch = select(gdb, geneid, :epoch)[1]
    end
    geneid
end

function evolve!(gdb::GeneDB, parpop::Population, params::ModelParameters, t::Int, termon::Int)

    # unpacking parameters
    n = params.popsize
    heterofit = params.heterozygousfitness
    homofit = params.homozygousfitness
    # The last element is just there as a placeholder, and its value does not affect runs.
    recombs = [params.recombinationrates; 0.0]
    muts = params.mutationrates
    nloci = params.numberofloci

    # normalize mutation rates
    maxfit = max(heterofit, homofit)
    heterofit /= maxfit
    homofit /= maxfit

    mutarray = Array(Bool, nloci, 2) # boolean value for each gene if it will be mutated.
    ps = Array(Int, 2) # indices of parents of an offspring.
    parchrs = Array(Int, 2) # a gene from which chromosome is passed on to offspring.

    ncoals = [0 for _ = 1:nloci]
    lastcoals = [1 for _ = 1:nloci]

    chpop = Population(n, nloci)

    gen = 1 # current generation
    for gen = 1:t
        for i = 1:n # iterate over offspring
            while true
                getparentids!(ps, n)
                # determine if mutations occur.
                for locus = 1:nloci, chr = 1:2
                    mutarray[locus, chr] = rand() < muts[locus] ? true : false
                end
                # process the first locus, which is under selection. A offspring is homozygous only when it
                # inherits identical-by-state genes from both parents without mutation. Otherwise, the offspring
                # is heterozygous.
                rand!(parchrs, 1:2)
                if isidbystate(gdb, parpop[ps[1], 1, parchrs[1]], parpop[ps[2], 1, parchrs[2]]) &&
                    mutarray[1,1] == mutarray[1,2] == false

                    if rand() > homofit
                        continue
                    end

                else
                    if rand() > heterofit
                        continue
                    end
                end

                for par = 1:2,  locus = 1:nloci
                    if mutarray[locus, par]
                        chpop[i, locus, par] = mutate!(gdb, gen, parpop[ps[par], locus, parchrs[par]])
                    else
                        chpop[i, locus, par] = parpop[ps[par], locus, parchrs[par]]
                    end
                    parchrs[par] = rand() < recombs[locus] ? 3 - parchrs[par] : parchrs[par]
                end
                break
            end
        end
        parpop, chpop = chpop, parpop
        for locus in 1:nloci
            if hascoalesced(gdb, parpop, locus, lastcoals[locus])
                lastcoals[locus] = gen
                ncoals[locus] += 1
                reset!(gdb, parpop, gen, locus)
                write(STDERR, string("INFO: turn over ", ncoals[locus], " of locus ", locus, " on generation ", gen, "\n"))
            end
        end
        if termon == minimum(ncoals)
            break
        end
    end
    parpop, gen
end

function initialize!(pop::Population)
    # Initialize a parental population. Genes are distinct.
    gdb = GeneDB()
    for org in pop, locus = 1:nloci(pop), chr = 1:2
        org[locus, chr] = insert!(WithNewAllele, gdb, 1, 0)
    end
    gdb
end

function reset!(gdb::GeneDB, pop::Population, gen::Int, locus::Int)
    # Reset lineage info of genes so that all genes can be distinguished even for identical allelic state.
    for org in pop, chr = 1:2
        org[locus, chr] = insert!(WithoutNewAllele, gdb, gen, org[locus, chr])
    end
    nothing
end

function reinitialize!(pop::Population, oldgdb::GeneDB)
    gdb = GeneDB()
    for org in pop, locus = 1:nloci(pop), chr = 1:2
        org[locus, chr] = insert!(WithoutNewAllele, gdb, 1, 0, select(oldgdb, org[locus, chr], :state)[1])
    end
    gdb
end

function simulate(params::ModelParameters, burnin::Int, t::Int, turmon::Int)

    # This is a parental population, a population of offspring is created within evolve! function.
    pop, t = Population(params.popsize, params.numberofloci)

    # Burnin
    # Execute the exact-same sequence as main-loop of evolution and throws out lineage information afterwords.
    # This loop runs exacctly "burnin" generations regardless of the presence of coalescence.
    gdb = initialize!(pop) # all genes are distinct
    pop = evolve!(gdb, pop, params, burnin, -1)

    # Main loop of evolution
    # This loop terminates upon the first coalescence or after "t" generations.
    gdb = reinitialize!(pop, gdb)
    pop, t = evolve!(gdb, pop, params, t, turmon)
    pop, gdb, t
end

function toarray(gdb::GeneDB, pop::Population, field::Symbol)
    Int[select(gdb, org[locus, chr], field)[1] for org in pop, locus = 1:nloci(pop), chr = 1:2]
end

function counts(gdb::GeneDB, pop::Population)
    alleles = toarray(gdb, pop, :state)
    nl = nloci(pop)
    # allele count
    adata = [Dict{Int, Int}() for _ = 1:nl]
    # genotype count
    gdata = [Dict{NTuple{2, Int}, Int}() for _ = 1:nl]
    # haplotype count
    hdata = Dict{NTuple{nl, Int}, Int}()
    for org in 1:length(pop)
        for locus = 1:nl
            g = tuple(sort(vec(alleles[org, locus, :]))...)
            gdata[locus][g] = get(gdata[locus], g, 0) + 1
            for chr = 1:2
                a = alleles[org, locus, chr]
                adata[locus][a] = get(adata[locus], a, 0) + 1
            end
        end
        for chr = 1:2
            h = tuple(vec(alleles[org, :, chr])...)
            hdata[h] = get(hdata, h, 0) + 1
        end
    end
    adata, gdata, hdata
end

function spectra(gdb::GeneDB, pop::Population)
    adata, gdata, hdata = counts(gdb, pop)

    nl = nloci(pop)
    afs = [Dict{Int, Int}() for _ = 1:nl]
    gfs = [Dict{Int, Int}() for _ = 1:nl]
    hfs = Dict{Int, Int}()
    for locus = 1:nl
        for v in values(adata[locus])
            afs[locus][v] = get(afs[locus], v, 0) + 1
        end
    end
    for locus = 1:nl
        for v in values(gdata[locus])
            gfs[locus][v] = get(gfs[locus], v, 0) + 1
        end
    end
    for v in values(hdata)
        hfs[v] = get(hfs, v, 0) + 1
    end
    afs, gfs, hfs
end

function history(gdb::GeneDB, idx::Int, lastcoal::Int)
    val = [idx]
    epoch = select(gdb, idx, :epoch)[1]
    while epoch > lastcoal
        idx = select(gdb, idx, :parent)[1]
        epoch = select(gdb, idx, :epoch)[1]
        push!(val, idx)
    end
    reverse(val)
end

function distances(gdb::GeneDB, pop::Population, lastcoal::Int)
    alleles = toarray(gdb, pop, :id)
    states = Dict{Int, Int}()
    for org in pop, locus = 1:nloci(pop), chr = 1:2
        oid = org[locus, chr]
        states[oid] = select(gdb, oid, :state)[1]
    end

    dists = Array(Int, 4, 0)

    idx = 1

    for locus = 1:nloci(pop)
        lineages = sort(unique(alleles[:,locus,:]))
        dists = hcat(dists, Array(Int, 4, binomial(length(lineages), 2)))
        nl = length(lineages)
        h = Dict{Int, Vector{Int}}()
        for i = 1:nl
            l = lineages[i]
            h[l] = history(gdb, l, lastcoal)
        end
        for i = 1:(nl-1), j = (i+1):nl
            h1, h2 = h[lineages[i]], h[lineages[j]]
            ca = intersect(h1, h2)
            dists[1, idx] = locus
            dists[2, idx] = states[lineages[i]]
            dists[3, idx] = states[lineages[j]]
            if length(ca) > 0
                locca1 = findfirst(h1, ca[end])
                locca2 = findfirst(h2, ca[end])
                dists[4, idx] = length(h1) - locca1 + length(h2) - locca2
            else
                dists[4, idx] = -1
            end
            idx += 1
        end
    end
    dists'
end

function nsegsites(gdb::GeneDB, pop::Population, lastcoal::Int)
    # The number of segregating sites in the number of mutations in a tree up to the most recent common ancestor.
    lineages = toarray(gdb, pop, :id)

    nl = nloci(pop)
    nms = Array(Int, nl)
    for locus = 1:nl
        # First, check if every loci have MRCA. If not, number of mutations is not well-defined. In this case, simply return -1.
        if !hascoalesced(gdb, pop, locus, lastcoal)
            nms[locus] = -1
            continue
        end

        # Then, construct a set of mutations along a tree up until MRCA. This will include
        # MRCA. However, this should not be in the final number, as all samples have the same allele.
        ids = unique(lineages[:, locus, :])
        hists = Vector{Int}[history(gdb, id, lastcoal) for id in ids]
        commons = copy(hists[1])
        for h in hists
            commons = intersect(commons, h)
        end
        ca = commons[end]
        muts = Int[]
        for h in hists
            muts = union(muts, h[findfirst(h, ca):end])
        end
        nms[locus] = length(muts) - 1
    end
    nms
end

function tmrca(gdb::GeneDB, pop::Population)
    # Computes time until the most recent common ancestor.
    nl= nloci(pop)
    ts = Array(Int, nl)

    for locus = 1:nl
        path = IntSet(history(gdb, pop[1, locus, 1], 0))
        for org in pop, chr = 1:2
            intersect!(path, IntSet(history(gdb, org[locus, chr], 0)))
        end
        ts[locus] = select(gdb, maximum(path), :epoch)[1]
    end
    ts
end

end


text:
"""
	struct ExtractVector{T}
		item::Int
	end

represents an array of a fixed length, typically a feature vector of numbers of type T

```jloctest
julia> sc = ExtractVector(4)
julia> sc([2,3,1,4]).data
3×1 Array{Float32,2}:
 2.0
 3.0
 1.0
```
"""
struct ExtractVector{T} <: AbstractExtractor
	n::Int
end
ExtractVector(n::Int) = ExtractVector{FloatType}(n)

(s::ExtractVector{T})(::V) where {T,V<:Union{Missing, Nothing}} = ArrayNode(fill(missing, s.n, 1))
(s::ExtractVector{T})(::ExtractEmpty) where {T}= ArrayNode(Matrix{T}(undef, s.n, 0))
(s::ExtractVector)(v) = s(missing)
function (s::ExtractVector{T})(v::V) where {T,V<:Vector}
	isempty(v) && return s(missing)
	if length(v) > s.n
		@warn "array too long, truncating"
		x = reshape(T.(v[1:s.n]), :, 1)
		return(ArrayNode(x))
	elseif length(v) < s.n
		x = Matrix{Union{Missing, T}}(missing, s.n, 1)
		x[1:length(v)] .= v
		return(ArrayNode(x))
	else
		x = reshape(T.(v), :, 1)
		return(ArrayNode(x))
	end
end

Base.length(e::ExtractVector) = e.n
Base.hash(e::ExtractVector, h::UInt) = hash(e.n, h)
Base.:(==)(e1::ExtractVector, e2::ExtractVector) = e1.n === e2.n


text:
module test_loops

using LinearAlgebra
using SparseArrays
using SparseMatricesCSR
using ThreadedSparseCSR
using BenchmarkTools

struct ThreadBuff{IT,JT,VT}
    # mrange::UnitRange{Int64}
    I::IT
    J::JT
    V::VT
    result::Vector{Float64}
end

function coomul!(A, IB, JB, VB, C)
    nz = length(IB)
    @assert length(JB) == nz
    @assert length(VB) == nz
    A .= zero(eltype(A))
    @inbounds @simd for i in 1:nz
        ir = IB[i]
        jc = JB[i]
        A[ir] += VB[i] * C[jc]
    end
end

function parcoomul!(R, IB, JB, VB, tbs, U)
    function innercoomul!(A, IB, JB, VB, C)
        nz = length(IB)
        @assert length(JB) == nz
        @assert length(VB) == nz
        A .= zero(eltype(A))
        @inbounds @simd for i in 1:nz
            ir = IB[i]
            jc = JB[i]
            A[ir] += VB[i] * C[jc]
        end
    end
    tasks = []
    for th in 1:length(tbs)
        push!(tasks, Threads.@spawn begin
            innercoomul!(tbs[th].result, tbs[th].I, tbs[th].J, tbs[th].V, U)
        end)
    end
    Threads.wait(tasks[1])
    R .= tbs[1].result
    for th in 2:length(tbs)
        Threads.wait(tasks[th])
        R .+= tbs[th].result
    end
    R
end

function parcsrmul!(y::AbstractVector, A::SparseMatrixCSR, v::AbstractVector)
  A.n == size(v, 1) || throw(DimensionMismatch())
  A.m == size(y, 1) || throw(DimensionMismatch())
  fill!(y, zero(eltype(y)))
  o = getoffset(A)
  Threads.@threads for row = 1:size(y, 1)
    @inbounds for nz in nzrange(A,row)
      col = A.colval[nz]+o
      y[row] += A.nzval[nz]*v[col]
    end
  end
  return y
end

function _execute()
    N = 2000000
    K = sprand(N, N, 0.00004)
    K = K + K'
    # @show K
    # K = sparse([5, 5, 10, 4, 3, 6, 9, 10, 1, 2, 8, 4, 5, 10, 4, 2, 4, 8], [1, 2, 2, 3, 4, 4, 4, 4, 5, 5, 5, 6, 8, 8, 9, 10, 10, 10], [2.27568e-01, 1.10446e-01, 5.29148e-01, 3.24832e-01, 3.24832e-01, 6.16751e-03, 8.60792e-01, 1.01654e+00, 2.27568e-01, 1.10446e-01, 7.70209e-01, 6.16751e-03, 7.70209e-01, 9.17418e-01, 8.60792e-01, 5.29148e-01, 1.01654e+00, 9.17418e-01], 10, 10)
    U = rand(N)
    R = zeros(N)

    @show nth = Base.Threads.nthreads()
    nz = nnz(K)
    @show chunk = Int(floor(nz / nth))

    IK, JK, VK = findnz(K)

    threadbuffs = ThreadBuff[]
    for th in 1:nth
        mrange = th < nth ? (chunk*(th-1)+1:chunk*(th)+1-1) : (chunk*(th-1)+1:nz)
        push!(threadbuffs, ThreadBuff(view(IK, mrange), view(JK, mrange), view(VK, mrange), deepcopy(R)))
    end

    Kr = SparseMatricesCSR.sparsecsr(IK, JK, VK, N, N)

    oU = deepcopy(U)
    oU .= rand()

    U .= oU
    mul!(R, K, U)
    R1 = deepcopy(R)

    U .= oU
    coomul!(R, IK, JK, VK, U)
    @show norm(R - R1) / norm(R1)

    U .= oU
    parcoomul!(R, IK, JK, VK, threadbuffs, U)
    @show norm(R - R1) / norm(R1)

    U .= oU
    parcsrmul!(R, Kr, U)
    @show norm(R - R1) / norm(R1)

    @btime ThreadedSparseCSR.tmul!($R, $Kr, $U)
    @btime ThreadedSparseCSR.bmul!($R, $Kr, $U)
    @btime mul!($R, $K, $U)
    @btime $coomul!($R, $IK, $JK, $VK, $U)
    @btime $parcoomul!($R, $IK, $JK, $VK, $threadbuffs, $U)
    @btime SparseMatricesCSR.mul!($R, $Kr, $U)
    @btime $parcsrmul!($R, $Kr, $U)

    nothing
end

function allrun()
    _execute()
    return true
end # function allrun

end # module
nothing

 using .Main.test_loops
 test_loops.allrun()


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MOI.Test
const MOIU = MOI.Utilities
const MOIB = MOI.Bridges

import SDPA
const optimizer = SDPA.Optimizer()
MOI.set(optimizer, MOI.Silent(), true)

@testset "SolverName" begin
    @test MOI.get(optimizer, MOI.SolverName()) == "SDPA"
end

@testset "supports_default_copy_to" begin
    @test MOIU.supports_allocate_load(optimizer, false)
    @test !MOIU.supports_allocate_load(optimizer, true)
end

# UniversalFallback is needed for starting values, even if they are ignored by SDPA
const cache = MOIU.UniversalFallback(MOIU.Model{Float64}())
const cached = MOIU.CachingOptimizer(cache, optimizer)
const bridged = MOIB.full_bridge_optimizer(cached, Float64)
# test 1e-3 because of rsoc3 test, otherwise, 1e-5 is enough
const config = MOIT.TestConfig(atol=1e-3, rtol=1e-3)

@testset "Unit" begin
    MOIT.unittest(bridged, config, [
        # `NumberOfThreads` not supported.
        "number_threads",
        # `TimeLimitSec` not supported.
        "time_limit_sec",
        # SingleVariable objective of bridged variables, will be solved by objective bridges
        "solve_time", "raw_status_string", "solve_singlevariable_obj",
        # Quadratic functions are not supported
        "solve_qcp_edge_cases", "solve_qp_edge_cases",
        # Integer and ZeroOne sets are not supported
        "solve_integer_edge_cases", "solve_objbound_edge_cases",
        "solve_zero_one_with_bounds_1",
        "solve_zero_one_with_bounds_2",
        "solve_zero_one_with_bounds_3"])
end
@testset "Linear tests" begin
    # See explanation in `MOI/test/Bridges/lazy_bridge_optimizer.jl`.
    # This is to avoid `Variable.VectorizeBridge` which does not support
    # `ConstraintSet` modification.
    MOIB.remove_bridge(bridged, MOIB.Constraint.ScalarSlackBridge{Float64})
    MOIT.contlineartest(bridged, config, [
        # `MOI.UNKNOWN_RESULT_STATUS` instead of `MOI.INFEASIBILITY_CERTIFICATE`
        "linear8a",
        "linear12"
    ])
end
@testset "Conic tests" begin
    MOIT.contconictest(bridged, config, [
        "lin3", "soc3", "normone2", "norminf2",
        # Missing bridges
        "rootdets",
        # Does not support power and exponential cone
        "pow", "dualpow", "logdet", "exp", "dualexp"])
end


text:
# This script downloads and bundles CA certificates from Amazon Web Services (AWS).
# Right now only the RDS root certificate is included, but more could be added in the future.

# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder

name = "AmazonCACerts"
# Info here: https://docs.aws.amazon.com/AmazonRDS/latest/UserGuide/UsingWithRDS.SSL.html
certificate_version = "2019"
certificate_filename = "rds-ca-$certificate_version-root.pem"
tarball_version = VersionNumber(certificate_version)

# Collection of sources required to build AmazonCACerts
sources = [
    FileSource("https://s3.amazonaws.com/rds-downloads/$certificate_filename",
    "d464378fbb8b981d2b28a1deafffd0113554e6adfb34535134f411bf3c689e73",
    filename="$certificate_filename")
]

# Bash recipe for building across all platforms
script = """
cd \$WORKSPACE/srcdir/
mkdir -p \$prefix/share
cp $certificate_filename \$prefix/share/$certificate_filename
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [AnyPlatform()]

# The products that we will ensure are always built
products = [
    FileProduct("share/$certificate_filename", :rds_ca_root_certificate)
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, tarball_version, sources, script, platforms, products, dependencies)


text:
"""
    @sync ex

Run expression `ex` and synchronize the GPU afterwards.

See also: [`synchronize`](@ref).
"""
macro sync(code)
    quote
        local ret = $(esc(code))
        synchronize()
        ret
    end
end

function darwin_version()
    # extract the trailing `-darwinXXX` bit from the triple
    machine = Sys.MACHINE
    VersionNumber(machine[findfirst("darwin", machine)[end]+1:end])
end

const _macos_version = Ref{VersionNumber}()
function macos_version()
    if !isassigned(_macos_version)
        verstr = read(`sw_vers -productVersion`, String)
        _macos_version[] = parse(VersionNumber, verstr)
    end
    _macos_version[]
end

function versioninfo(io::IO=stdout)
    println(io, "macOS $(macos_version()), Darwin $(darwin_version())")
    println(io)

    println(io, "Toolchain:")
    println(io, "- Julia: $VERSION")
    println(io, "- LLVM: $(LLVM.version())")
    println(io)

    env = filter(var->startswith(var, "JULIA_METAL"), keys(ENV))
    if !isempty(env)
        println(io, "Environment:")
        for var in env
            println(io, "- $var: $(ENV[var])")
        end
        println(io)
    end

    devs = devices()
    if isempty(devs)
        println(io, "No Metal devices.")
    elseif length(devs) == 1
        println(io, "1 device:")
    else
        println(io, length(devs), " devices:")
    end
    for (i, dev) in enumerate(devs)
        println(io, "- $(dev.name) ($(Base.format_bytes(dev.currentAllocatedSize)) allocated)")
    end

    return
end


text:
@testset "ZeroMeanUnitVarianceScaling" begin
    random_datas = [rand(10), rand(5,5), rand(3,3,3)]

    for random_data in random_datas
        scaling = ZeroMeanUnitVarianceScaling(random_data)

        @test scaling.μ == mean(random_data)
        @test scaling.σ == std(random_data)

        @test mean(scaling.(random_data)) ≈ 0 atol=1e-10
        @test std(scaling.(random_data)) ≈ 1
        @test all(inv(scaling).(scaling.(random_data)) .≈ random_data)
    end
end

@testset "MinMaxScaling" begin
    random_datas = [rand(10), rand(5,5), rand(3,3,3)]

    ranges = ((a=0, b=10), (a=-5, b=-1), (a=1, b=2), (a=-1, b=20))

    for random_data in random_datas, range in ranges
        a = range.a
        b = range.b
        scaling = MinMaxScaling(random_data, a=a, b=b)

        @test scaling.a == a
        @test scaling.b == b
        @test scaling.data_min == minimum(random_data)
        @test scaling.data_max == maximum(random_data)

        @test minimum(scaling.(random_data)) ≈ a atol=1e-10
        @test maximum(scaling.(random_data)) ≈ b atol=1e-10
        @test all(inv(scaling).(scaling.(random_data)) .≈ random_data)
    end
end


text:
using ZXCalculus, LightGraphs
using Test
using ZXCalculus: Phase

struct TestZXDiagram{T, P} <: AbstractZXDiagram{T, P} end

test_zxd = TestZXDiagram{Int, Phase}();

@test_throws MethodError LightGraphs.nv(test_zxd)
@test_throws MethodError LightGraphs.ne(test_zxd)
@test_throws MethodError LightGraphs.degree(test_zxd, 1)
@test_throws MethodError LightGraphs.indegree(test_zxd, 1)
@test_throws MethodError LightGraphs.outdegree(test_zxd, 1)
@test_throws MethodError LightGraphs.neighbors(test_zxd, 1)
@test_throws MethodError LightGraphs.outneighbors(test_zxd, 1)
@test_throws MethodError LightGraphs.inneighbors(test_zxd, 1)
@test_throws MethodError LightGraphs.rem_edge!(test_zxd, 1, 2)
@test_throws MethodError LightGraphs.add_edge!(test_zxd, 1, 2, 1)
@test_throws MethodError LightGraphs.has_edge(test_zxd, 1, 2)

@test_throws MethodError print(test_zxd)
@test_throws MethodError Base.copy(test_zxd)

@test_throws MethodError ZXCalculus.nqubits(test_zxd)
@test_throws MethodError spiders(test_zxd)
@test_throws MethodError tcount(test_zxd)
@test_throws MethodError ZXCalculus.get_inputs(test_zxd)
@test_throws MethodError ZXCalculus.get_outputs(test_zxd)
@test_throws MethodError ZXCalculus.scalar(test_zxd)
@test_throws MethodError ZXCalculus.spider_sequence(test_zxd)
@test_throws MethodError ZXCalculus.round_phases!(test_zxd)

@test_throws MethodError spider_type(test_zxd, 1)
@test_throws MethodError phase(test_zxd, 1)
@test_throws MethodError rem_spider!(test_zxd, 1)
@test_throws MethodError rem_spiders!(test_zxd, [1, 2])
@test_throws MethodError ZXCalculus.qubit_loc(test_zxd, 1)
@test_throws MethodError ZXCalculus.column_loc(test_zxd, 1)
@test_throws MethodError ZXCalculus.add_global_phase!(test_zxd, 3)
@test_throws MethodError ZXCalculus.add_power!(test_zxd, 4)
@test_throws MethodError ZXCalculus.generate_layout!(test_zxd, [])

@test_throws MethodError ZXCalculus.set_phase!(test_zxd, 1, Phase(1//1))
@test_throws MethodError push_gate!(test_zxd, 1, Val(:X), Phase(1//2))
@test_throws MethodError pushfirst_gate!(test_zxd, 1, 2, Val(:CNOT))
@test_throws MethodError ZXCalculus.add_spider!(test_zxd, Phase(1//1))
@test_throws MethodError ZXCalculus.insert_spider!(test_zxd, Phase(1//2), [2, 3])


text:
@with_kw mutable struct DIRCOLSolverMTOptions{T} <: DirectSolverOptions{T}
    "NLP Solver to use. See MathOptInterface for available NLP solvers"
    nlp::MathOptInterface.AbstractOptimizer = Ipopt.Optimizer()

    "Options dictionary for the nlp solver"
    opts::Dict{Symbol,Any} = Dict{Symbol,Any}()

    "Print output to console"
    verbose::Bool = true

    "Minimum Time Options"
    R_min_time::T = 1.0

    h_max::T = 1.0
    h_min::T = 1.0e-3
end

"""
$(TYPEDEF)
Direct Collocation Solver.
Uses a commerical NLP solver to solve the Trajectory Optimization problem.
"""
mutable struct DIRCOLSolverMT{T,Q} <: DirectSolver{T}
    opts::DIRCOLSolverMTOptions{T}
    stats::Dict{Symbol,Any}
    Z::PrimalsMT{T}
    X_::VectorTrajectory{T}
    ∇F::PartedMatTrajectory{T}
    C::PartedVecTrajectory{T}
    ∇C::PartedMatTrajectory{T}
    fVal::VectorTrajectory{T}
    p::Vector{Int}
end

DIRCOLSolverMT(prob::Problem, opts::DIRCOLSolverMTOptions=DIRCOLSolverMTOptions{Float64}(),
    Z::PrimalsMT{T}=PrimalsMT(prob,true)) where {T,Q} = AbstractSolver(prob, opts, Z)

type(::PrimalsMT{T}) where T = T

function AbstractSolver(prob::Problem, opts::DIRCOLSolverMTOptions, Z::PrimalsMT{T}=PrimalsMT(prob, true)) where T
    n,m,N = size(prob)
    X_ = [zeros(T,n) for k = 1:N-1] # midpoints

    part_f = create_partition2(prob.model)
    constraints = prob.constraints
    p = num_stage_constraints(constraints)
    c_stage = [stage(constraints[k]) for k = 1:N-1]
    c_part = [create_partition(c_stage[k]) for k = 1:N-1]
    c_part2 = [create_partition2(c_stage[k],n,m) for k = 1:N-1]

    # Create Trajectories
    ∇F         = [PartedMatrix(zeros(T,n,length(prob.model)),part_f)         for k = 1:N]
    C          = [PartedVector(T,constraints[k],:stage)     for k = 1:N-1]
    ∇C         = [PartedMatrix(T,constraints[k],n,m,:stage) for k = 1:N-1]
    C          = [C...,  PartedVector(T,constraints[N],:terminal)]
    ∇C         = [∇C..., PartedMatrix(T,constraints[N],n,m,:terminal)]

    c_term = terminal(constraints[N])
    p_N = num_constraints(c_term)
    fVal = [zeros(T,n) for k = 1:N]
    p = num_constraints(prob)

    solver = DIRCOLSolverMT{T,HermiteSimpson}(opts, Dict{Symbol,Any}(), Z, X_, ∇F, C, ∇C, fVal, p)
    reset!(solver)
    return solver
end

function reset!(solver::DIRCOLSolverMT{T,Q}) where {T, Q<:QuadratureRule}
    state = Dict{Symbol,Any}(:iterations=>0, :c_max=>T[], :cost=>T[])
    solver.stats[:iterations] = 0
    solver.stats[:c_max] = T[]
    solver.stats[:cost] = T[]
end

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


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

JLLWrappers.@generate_main_file_header("ranger")
JLLWrappers.@generate_main_file("ranger", UUID("0a592095-bd93-57bb-89c0-a7c16a918b1d"))
end  # module ranger_jll


text:
using TimeSeries
using Dates
using PowerSystems

dates = collect(
    DateTime("1/1/2024  0:00:00", "d/m/y  H:M:S"):Hour(1):DateTime(
        "1/1/2024  23:00:00",
        "d/m/y  H:M:S",
    ),
)

nodes14() = [
    Bus(1, "Bus 1", "REF", 0.0, 1.06, (min = 0.94, max = 1.06), 69, nothing, nothing),
    Bus(
        2,
        "Bus 2",
        "PV",
        -0.08691739674931762,
        1.045,
        (min = 0.94, max = 1.06),
        69,
        nothing,
        nothing,
    ),
    Bus(
        3,
        "Bus 3",
        "PV",
        -0.22200588085367873,
        1.01,
        (min = 0.94, max = 1.06),
        69,
        nothing,
        nothing,
    ),
    Bus(
        4,
        "Bus 4",
        "PQ",
        -0.18029251173101424,
        1.019,
        (min = 0.94, max = 1.06),
        69,
        nothing,
        nothing,
    ),
    Bus(
        5,
        "Bus 5",
        "PQ",
        -0.15323990832510212,
        1.02,
        (min = 0.94, max = 1.06),
        69,
        nothing,
        nothing,
    ),
    Bus(
        6,
        "Bus 6",
        "PV",
        -0.24818581963359368,
        1.07,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        7,
        "Bus 7",
        "PQ",
        -0.23335052099164186,
        1.062,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        8,
        "Bus 8",
        "PV",
        -0.2331759880664424,
        1.09,
        (min = 0.94, max = 1.06),
        18,
        nothing,
        nothing,
    ),
    Bus(
        9,
        "Bus 9",
        "PQ",
        -0.2607521902479528,
        1.056,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        10,
        "Bus 10",
        "PQ",
        -0.26354471705114374,
        1.051,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        11,
        "Bus 11",
        "PQ",
        -0.2581341963699613,
        1.057,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        12,
        "Bus 12",
        "PQ",
        -0.2630211182755455,
        1.055,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        13,
        "Bus 13",
        "PQ",
        -0.2645919146023404,
        1.05,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
    Bus(
        14,
        "Bus 14",
        "PQ",
        -0.27995081201989047,
        1.036,
        (min = 0.94, max = 1.06),
        13.8,
        nothing,
        nothing,
    ),
]
branches14_dc(nodes14) = [
    Line(
        "Line1",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[1], to = nodes14[2]),
        0.01938,
        0.05917,
        (from = 0.0264, to = 0.0264),
        18.046,
        1.04,
    ),
    Line(
        "Line2",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[1], to = nodes14[5]),
        0.05403,
        0.22304,
        (from = 0.0246, to = 0.0246),
        4.896,
        1.04,
    ),
    HVDCLine(
        "DCLine3",
        true,
        0.0,
        Arc(from = nodes14[2], to = nodes14[3]),
        (min = -600.0, max = 600),
        (min = -600.0, max = 600),
        (min = -600.0, max = 600),
        (min = -600.0, max = 600),
        (l0 = 0.01, l1 = 0.001),
    ),
    HVDCLine(
        "DCLine4",
        true,
        0.0,
        Arc(from = nodes14[2], to = nodes14[4]),
        (min = -600.0, max = 600),
        (min = -600.0, max = 600),
        (min = -600.0, max = 600),
        (min = -600.0, max = 600),
        (l0 = 0.01, l1 = 0.001),
    ),
    #Line("Line3",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[3]),   0.04699, 0.19797, (from=0.0219, to=0.0219), 5.522, 1.04),
    #Line("Line4",  true, 0.0, 0.0, Arc(from=nodes14[2],to=nodes14[4]),   0.05811, 0.17632, (from=0.017,  to=0.017), 6.052, 1.04),
    Line(
        "Line5",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[2], to = nodes14[5]),
        0.05695,
        0.17388,
        (from = 0.0173, to = 0.0173),
        6.140,
        1.04,
    ),
    Line(
        "Line6",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[3], to = nodes14[4]),
        0.06701,
        0.17103,
        (from = 0.0064, to = 0.0064),
        6.116,
        1.04,
    ),
    Line(
        "Line7",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[4], to = nodes14[5]),
        0.01335,
        0.04211,
        (from = 0.0, to = 0.0),
        25.434,
        1.04,
    ),
    TapTransformer(
        "Trans3",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[4], to = nodes14[7]),
        0.0,
        0.20912,
        0.0,
        0.978,
        20.0,
    ),
    TapTransformer(
        "Trans1",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[4], to = nodes14[9]),
        0.0,
        0.55618,
        0.0,
        0.969,
        20.0,
    ),
    TapTransformer(
        "Trans2",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[5], to = nodes14[6]),
        0.0,
        0.25202,
        0.0,
        0.932,
        20.0,
    ),
    Line(
        "Line8",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[6], to = nodes14[11]),
        0.09498,
        0.19890,
        (from = 0.0, to = 0.0),
        5.373,
        1.04,
    ),
    Line(
        "Line9",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[6], to = nodes14[12]),
        0.12291,
        0.25581,
        (from = 0.0, to = 0.0),
        2.020,
        1.04,
    ),
    Line(
        "Line10",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[6], to = nodes14[13]),
        0.06615,
        0.13027,
        (from = 0.0, to = 0.0),
        4.458,
        1.04,
    ),
    Transformer2W(
        "Trans4",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[7], to = nodes14[8]),
        0.0,
        0.17615,
        0.0,
        20.0,
    ),
    Line(
        "Line16",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[7], to = nodes14[9]),
        0.0,
        0.11001,
        (from = 0.0, to = 0.0),
        12.444,
        1.04,
    ),
    Line(
        "Line11",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[9], to = nodes14[10]),
        0.03181,
        0.08450,
        (from = 0.0, to = 0.0),
        5.097,
        1.04,
    ),
    Line(
        "Line12",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[9], to = nodes14[14]),
        0.12711,
        0.27038,
        (from = 0.0, to = 0.0),
        3.959,
        1.04,
    ),
    Line(
        "Line13",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[10], to = nodes14[11]),
        0.08205,
        0.19207,
        (from = 0.0, to = 0.0),
        7.690,
        1.04,
    ),
    Line(
        "Line14",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[12], to = nodes14[13]),
        0.22092,
        0.19988,
        (from = 0.0, to = 0.0),
        6.378,
        1.04,
    ),
    Line(
        "Line15",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[13], to = nodes14[14]),
        0.17093,
        0.34802,
        (from = 0.0, to = 0.0),
        10.213,
        1.04,
    ),
]

branches14(nodes14) = [
    Line(
        "Line1",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[1], to = nodes14[2]),
        0.01938,
        0.05917,
        (from = 0.0264, to = 0.0264),
        18.046,
        1.04,
    ),
    Line(
        "Line2",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[1], to = nodes14[5]),
        0.05403,
        0.22304,
        (from = 0.0246, to = 0.0246),
        4.896,
        1.04,
    ),
    Line(
        "Line3",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[2], to = nodes14[3]),
        0.04699,
        0.19797,
        (from = 0.0219, to = 0.0219),
        5.522,
        1.04,
    ),
    Line(
        "Line4",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[2], to = nodes14[4]),
        0.05811,
        0.17632,
        (from = 0.017, to = 0.017),
        6.052,
        1.04,
    ),
    Line(
        "Line5",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[2], to = nodes14[5]),
        0.05695,
        0.17388,
        (from = 0.0173, to = 0.0173),
        6.140,
        1.04,
    ),
    Line(
        "Line6",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[3], to = nodes14[4]),
        0.06701,
        0.17103,
        (from = 0.0064, to = 0.0064),
        6.116,
        1.04,
    ),
    Line(
        "Line7",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[4], to = nodes14[5]),
        0.01335,
        0.04211,
        (from = 0.0, to = 0.0),
        25.434,
        1.04,
    ),
    TapTransformer(
        "Trans3",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[4], to = nodes14[7]),
        0.0,
        0.20912,
        0.0,
        0.978,
        20.0,
    ),
    TapTransformer(
        "Trans1",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[4], to = nodes14[9]),
        0.0,
        0.55618,
        0.0,
        0.969,
        20.0,
    ),
    TapTransformer(
        "Trans2",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[5], to = nodes14[6]),
        0.0,
        0.25202,
        0.0,
        0.932,
        20.0,
    ),
    Line(
        "Line8",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[6], to = nodes14[11]),
        0.09498,
        0.19890,
        (from = 0.0, to = 0.0),
        5.373,
        1.04,
    ),
    Line(
        "Line9",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[6], to = nodes14[12]),
        0.12291,
        0.25581,
        (from = 0.0, to = 0.0),
        2.020,
        1.04,
    ),
    Line(
        "Line10",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[6], to = nodes14[13]),
        0.06615,
        0.13027,
        (from = 0.0, to = 0.0),
        4.458,
        1.04,
    ),
    Transformer2W(
        "Trans4",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[7], to = nodes14[8]),
        0.0,
        0.17615,
        0.0,
        20.0,
    ),
    Line(
        "Line16",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[7], to = nodes14[9]),
        0.0,
        0.11001,
        (from = 0.0, to = 0.0),
        12.444,
        1.04,
    ),
    Line(
        "Line11",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[9], to = nodes14[10]),
        0.03181,
        0.08450,
        (from = 0.0, to = 0.0),
        5.097,
        1.04,
    ),
    Line(
        "Line12",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[9], to = nodes14[14]),
        0.12711,
        0.27038,
        (from = 0.0, to = 0.0),
        3.959,
        1.04,
    ),
    Line(
        "Line13",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[10], to = nodes14[11]),
        0.08205,
        0.19207,
        (from = 0.0, to = 0.0),
        7.690,
        1.04,
    ),
    Line(
        "Line14",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[12], to = nodes14[13]),
        0.22092,
        0.19988,
        (from = 0.0, to = 0.0),
        6.378,
        1.04,
    ),
    Line(
        "Line15",
        true,
        0.0,
        0.0,
        Arc(from = nodes14[13], to = nodes14[14]),
        0.17093,
        0.34802,
        (from = 0.0, to = 0.0),
        10.213,
        1.04,
    ),
]

thermal_generators14(nodes14) = [
    ThermalStandard(
        "Bus1",
        true,
        nodes14[1],
        2.0,
        -0.169,
        TechThermal(
            2.324,
            PrimeMovers.ST,
            ThermalFuels.COAL,
            (min = 0.0, max = 3.332),
            (min = 0.0, max = 0.1),
            nothing,
            nothing,
        ),
        ThreePartCost((430.292599, 2000.0), 0.0, 0.0, 0.0),
    ),
    ThermalStandard(
        "Bus2",
        true,
        nodes14[2],
        0.40,
        0.42,
        TechThermal(
            1.4,
            PrimeMovers.ST,
            ThermalFuels.COAL,
            (min = 0.0, max = 1.40),
            (min = -0.4, max = 0.5),
            nothing,
            nothing,
        ),
        ThreePartCost((2500.0, 2000.0), 0.0, 0.0, 0.0),
    ),
    ThermalStandard(
        "Bus3",
        true,
        nodes14[3],
        0.0,
        0.23,
        TechThermal(
            1.0,
            PrimeMovers.ST,
            ThermalFuels.COAL,
            (min = 0.0, max = 1.0),
            (min = 0.0, max = 0.4),
            nothing,
            nothing,
        ),
        ThreePartCost((100.0, 4000.0), 0.0, 0.0, 0.0),
    ),
    ThermalStandard(
        "Bus6",
        true,
        nodes14[6],
        0.0,
        0.12,
        TechThermal(
            1.0,
            PrimeMovers.ST,
            ThermalFuels.COAL,
            (min = 0.0, max = 1.0),
            (min = -0.06, max = 0.24),
            nothing,
            nothing,
        ),
        ThreePartCost((100.0, 4000.0), 0.0, 0.0, 0.0),
    ),
    ThermalStandard(
        "Bus8",
        true,
        nodes14[8],
        0.0,
        0.174,
        TechThermal(
            1.0,
            PrimeMovers.ST,
            ThermalFuels.COAL,
            (min = 0.0, max = 1.0),
            (min = -0.06, max = 0.24),
            nothing,
            nothing,
        ),
        ThreePartCost((100.0, 4000.0), 0.0, 0.0, 0.0),
    ),
]

loadz1_ts = [
    0.792729978
    0.723201574
    0.710952098
    0.677672816
    0.668249175
    0.67166919
    0.687608809
    0.711821241
    0.756320618
    0.7984057
    0.827836527
    0.840362459
    0.84511032
    0.834592803
    0.822949221
    0.816941743
    0.824079963
    0.905735139
    0.989967048
    1
    0.991227765
    0.960842114
    0.921465115
    0.837001437
]

loadz2_ts = [
    0.831093782
    0.689863228
    0.666058513
    0.627033103
    0.624901388
    0.62858924
    0.650734211
    0.683424321
    0.750876413
    0.828347191
    0.884248576
    0.888523615
    0.87752169
    0.847534405
    0.8227661
    0.803809323
    0.813282799
    0.907575962
    0.98679848
    1
    0.990489904
    0.952520972
    0.906611479
    0.824307054
]

loadz3_ts = [
    0.871297342
    0.670489749
    0.642812243
    0.630092987
    0.652991383
    0.671971681
    0.716278493
    0.770885833
    0.810075243
    0.85562361
    0.892440566
    0.910660449
    0.922135467
    0.898416969
    0.879816542
    0.896390855
    0.978598576
    0.96523761
    1
    0.969626503
    0.901212601
    0.81894251
    0.771004923
    0.717847996
]

loads14(nodes14) = [
    PowerLoad(
        "Bus2",
        true,
        nodes14[2],
        LoadModels.ConstantPower,
        0.217,
        0.127,
        0.217,
        0.127,
    ),
    PowerLoad("Bus3", true, nodes14[3], LoadModels.ConstantPower, 0.942, 0.19, 0.942, 0.19),
    PowerLoad(
        "Bus4",
        true,
        nodes14[4],
        LoadModels.ConstantPower,
        0.478,
        -0.039,
        0.478,
        -0.039,
    ),
    PowerLoad(
        "Bus5",
        true,
        nodes14[5],
        LoadModels.ConstantPower,
        0.076,
        0.016,
        0.076,
        0.016,
    ),
    PowerLoad(
        "Bus6",
        true,
        nodes14[6],
        LoadModels.ConstantPower,
        0.112,
        0.075,
        0.112,
        0.075,
    ),
    PowerLoad(
        "Bus9",
        true,
        nodes14[9],
        LoadModels.ConstantPower,
        0.295,
        0.166,
        0.295,
        0.166,
    ),
    PowerLoad(
        "Bus10",
        true,
        nodes14[10],
        LoadModels.ConstantPower,
        0.09,
        0.058,
        0.09,
        0.058,
    ),
    PowerLoad(
        "Bus11",
        true,
        nodes14[11],
        LoadModels.ConstantPower,
        0.035,
        0.018,
        0.035,
        0.018,
    ),
    PowerLoad(
        "Bus12",
        true,
        nodes14[12],
        LoadModels.ConstantPower,
        0.061,
        0.016,
        0.061,
        0.016,
    ),
    PowerLoad(
        "Bus13",
        true,
        nodes14[13],
        LoadModels.ConstantPower,
        0.135,
        0.058,
        0.135,
        0.058,
    ),
    PowerLoad(
        "Bus14",
        true,
        nodes14[14],
        LoadModels.ConstantPower,
        0.149,
        0.050,
        0.149,
        0.050,
    ),
]

timeseries_DA14 = [
    TimeArray(dates, loadz1_ts),
    TimeArray(dates, loadz1_ts),
    TimeArray(dates, loadz3_ts),
    TimeArray(dates, loadz1_ts),
    TimeArray(dates, loadz2_ts),
    TimeArray(dates, loadz3_ts),
    TimeArray(dates, loadz2_ts),
    TimeArray(dates, loadz2_ts),
    TimeArray(dates, loadz2_ts),
    TimeArray(dates, loadz2_ts),
    TimeArray(dates, loadz2_ts),
];

battery14(nodes14) = [GenericBattery(
    name = "Bat",
    primemover = PrimeMovers.BA,
    available = true,
    bus = nodes14[1],
    energy = 5.0,
    capacity = (min = 5.0, max = 100.0),
    rating = 70,
    activepower = 10.0,
    inputactivepowerlimits = (min = 0.0, max = 50.0),
    outputactivepowerlimits = (min = 0.0, max = 50.0),
    reactivepower = 0.0,
    reactivepowerlimits = (min = -50.0, max = 50.0),
    efficiency = (in = 0.80, out = 0.90),
)]


text:
using JuSwarm, Test, SafeTestsets
@time begin
@time @safetestset "Sphere Function" begin include("sphere_function_test.jl") end
@time @safetestset "Rosenbrock Function" begin include("rosenbrock_function_test.jl") end
@time @safetestset "Ackley Function" begin include("ackley_function_test.jl") end
@time @safetestset "Rastrigin Function" begin include("rastrigin_function_test.jl") end
end

max_stars_repo_name: GrantHecht/JuSwarm

max_forks_repo_name: GrantHecht/JuSwarm


text:
"""
Random graph model exhibiting low GCC and high ALCC. Creates a fixed number
of cliques each with a random number of nodes. Nodes can partake in multiple
cliques.

Based on "Random clique covers for graphs with local density and global sparsity"
by Sinead A. Williamson and Mauricio Tec. The original paper may be found
at https://arxiv.org/pdf/1704.03913.pdf.
"""

using MatrixNetworks
using SpecialFunctions # Used for gamma
using LinearAlgebra # Matrix transpose
using SparseArrays
using Distributions: Poisson

"""
`rcc`
=====================

Generate a random graph by creating random clique covers.

Arguments
---------
    - `a::Float64 (> 0)`: The expected number of vertices per clique
    - `s::Float64 (∈ [0,1))`: Controls how close the degree distribution is to a power law
    - `c::Float64 (> -s)`: Controls (but not equal to) the expected number of cliques each vertex is a part of.
                           Larger c = less overlap
    - `N::Int64 (=10)`: The total number of cliques

Preconditions
-------------
    - `a > 0`
    - `s ∈ [0,1)`
    - `c > -s`
    - `N >= 0`

Examples
--------
~~~~
random_clique_cover(10.,0.9,0.1;N=100) # 100-clique graph with ~10 nodes per clique with high degree skew and medium clique overlap
random_clique_cover(50.,0.1,5.) # 10-clique graph with ~50 nodes per clique with more uniform degree distribution and high clique overlap
~~~~

The random cliques are generated using the following process. Imagine a fixed
number of restaurant customers coming up to an (infinitely long) buffet one by one.
The first customer takes the first n dishes, where n is random. Each customer
thereafter
    1) samples each already-tried dish at random, with the probability
    of taking a dish decreasing with the number of times it has been tried
    2) selects the next k (also random) dishes that have not yet been tried.

This graph model works by substituting customers for cliques and dishes for nodes. Bon appétit!
"""
function rcc(a::Float64, s::Float64, c::Float64; N::Int64=10)
    @assert a > 0
    @assert s >= 0 && s < 1
    @assert c > -s
    @assert N >= 0

    P(x) = rand(Poisson(x)) # Sample from Poisson distribution with λ = x

    ### Approximate Γ(n+c+s−1)/Γ(n+c) to avoid overflow with large n ###
    function gamma_ratio_approx(n) # Permits large N, since it avoids using large numbers
        radical = ((ℯ^(1-s)) * ((n+c+s)^(s+0.5)) / ((n+c+1)^1.5)) ^ (1/(n+c))
        return (radical * (n+c+s)/(n+c+1)) ^ (n+c) # Derived from Stirling's approximation for gamma
    end

    agammr = a*gamma(c+1) / gamma(c+s) # Constant term for determining the number of new nodes to add to a clique

    cliques = [[i for i = 1:P(a)]] # Stores cliques
    clique_count = [1 for i = 1:size(cliques[1],1)] # Stores how many cliques each node is part of

    for n = 1:N-1
        push!(cliques, []) # New empty clique
        for k = 1:size(clique_count,1) # Decide which of existing nodes to place in new clique
            if rand(Float64) <= (clique_count[k]-s) / (n+c)
                push!(cliques[n+1], k)
                clique_count[k] += 1
            end
        end

        for j = 1:P(agammr*gamma_ratio_approx(n)) # Add new nodes to graph and clique
            push!(cliques[n+1], size(clique_count,1)+1)
            push!(clique_count, 1)
        end
    end

    Z = zeros(Int64, N, size(clique_count,1)) # Convert clique list to edge clique cover.
                                              # Z[i,j] = 1 iff clique i contains node j.
    for i = 1:N
        for j in cliques[i]
            Z[i,j] = 1
        end
    end

    mat = Z' * Z # Convert edge clique cover to corresponding adjacency matrix
    mat -= Diagonal(mat) # Simplify matrix
    mat = min.(mat, 1)

    return MatrixNetwork(sparse(mat))
end


text:
############################################################################################
# External libraries

############################################################################################
# Load sub-container
include("parameterinfo.jl")
include("modelwrapper.jl")
include("tagged.jl")
include("objective.jl")
#!NOTE: Remove Soss dependency from ModelWrappers because of heavy deps. Can make separate BaytesSoss later on.
#include("_soss.jl")
############################################################################################
# Export


text:
using RetroSignalModel
using RetroSignalModel: RtgMTK
using Documenter

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

deploydocs(;
    branch="gh-pages",
    devbranch="main",
    repo="github.com/ntumitolab/RetroSignalModel.jl.git"
)


text:
using Documenter, MyExample

makedocs(sitename="My Documentation")

max_stars_repo_name: filippigianluca/MyExample

max_issues_repo_name: filippigianluca/MyExample


text:
# MD5 sum of a file
md5sum(fn) = split(readstring(`md5sum $fn`))[1]

## tools.jl tests
#################
const todelete = []
function tempfn(ext="")
    out = tempname()*ext
    push!(todelete, out)
    out
end
function delfls()
    for f in todelete
        if isfile(f)
            rm(f)
        end
    end
end

## .agr readers:
g1 = BM._read_agr("testfiles/wiki.agr")
g1b = BM._read_agr("testfiles/wiki.bin")
@test g1.v == [-9999 -9999 5 2
                -9999 20 100 36
                3 8 35 10
                32 42 50 6
                88 75 27 9
                13 5 1 -9999]
@test g1.nc==4
@test g1.nr==6
@test g1.xll==0.0
@test g1.yll==-40.0
@test g1.dx==50.0
@test g1.NA==-9999
# write and read back:
tfn1 = tempfn(".agr")
write_agr(g1, tfn1)
g2 = BM._read_agr(tfn1, Float64)
@test g1==g2
tfn2 = tempfn(".agr")
write_agr(g2, tfn2)
@test md5sum(tfn1)==md5sum(tfn2)

# write and read back binary:
tfn3 = tempfn(".bin")
write_agr(g1, tfn3)
g2 = BM._read_agr(tfn3)
g1_ = BM._read_agr("testfiles/wiki.agr", Float32)
@test g1_==g2
@test md5sum(tfn3)==md5sum("testfiles/wiki.bin")

# larger file:
g3 = BM._read_agr("testfiles/t2.bin", Float32)
@test g3.nc==121
@test g3.nr==81
@test g3.xll==678750.0f0
@test g3.yll==140250.0f0
@test g3.dx==25.0f0
@test g3.NA==-9999.0f0

# larger file:
# @test_throws ErrorException BM._read_agr("testfiles/t2.bin", Float64)

tfn4 = tempfn(".bin")
write_agr(g3, tfn4)
@test md5sum("testfiles/t2.bin")==md5sum(tfn4)

# Gridded
gg = Gridded(g1)
@test BM.AGR(gg, NA_agr=-9999)==g1

## .xyn reader
po = read_xyn("testfiles/poly.xyn")
@test length(po)==1
@test size(po[1])==(2,24)

po = read_xyn("testfiles/multipoly.xyzn", hasz=true)
@test length(po)==4
@test size(po[1])==(3,4)
@test size(po[2])==(3,5)
@test size(po[3])==(3,4)
@test size(po[4])==(3,3)

# tidy up temp-files
delfls()

## IN poly
@test BM.leftorright(0.5,0.5, 1,0,1,1)==-1
@test BM.leftorright(1.5,.5, 1,0,1,1)==1
@test BM.leftorright(1,0.5, 1,0,1,1)==0

poc = BM.concat_poly(po)
@test !inpoly([0.0,0], poc[1])

poly = Float64[0 0
               0 1
               1 1
               1 0
               0 0]'
p1 = [0.5, 0.5]
p2 = [0.5, 0.99]
p22 = [0.5, 1] # on edge
p23 = [0.5, 0] # on edge
p24 = [0, 0]   # on corner
p25 = [0, .4]   # on edge
p3 = [0.5, 1.1]

@test inpoly(p1, poly)
@test inpoly(p2, poly)
@test inpoly(p22, poly)
@test inpoly(p23, poly)
@test inpoly(p24, poly)
@test inpoly(p25, poly)
@test !inpoly(p3, poly)

# clockwise poly
poly = Float64[0 0
               1 0
               1 1
               0 1
               0 0]'

@test inpoly(p1, poly)
@test inpoly(p2, poly)
@test inpoly(p22, poly)
@test inpoly(p23, poly)
@test inpoly(p24, poly)
@test inpoly(p25, poly)
@test !inpoly(p3, poly)

# cross-over poly
poly = Float64[0 0
               1 0
               0 1
               1 1
               0 0]'
if VERSION>=v"0.5-"
    eval(:(@test_broken inpoly(p1, poly) )) # should be true
end
@test inpoly(p2, poly)
@test inpoly(p22, poly)
@test inpoly(p23, poly)
@test inpoly(p24, poly)
@test !inpoly(p25, poly) # different
@test !inpoly(p3, poly)

# with interior region
poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside interior poly: i.e. labeled as outside
@test !inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               # in-interior
               0.4 0.4
               0.4 0.2
               0.2 0.2
               0.2 0.4
               0.4 0.4
               # interior
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside in-interior poly
@test inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior
               0.1 0.1
               0.1 0.6
               0.6 0.6
               # in-interior
               0.4 0.4
               0.2 0.4
               0.2 0.2
               0.4 0.2
               0.4 0.4
               # interior
               0.6 0.6
               0.6 0.1
               0.1 0.1
               # exterior
               0 0
               0 1
               1 1
               1 0
               0 0]'
# inside in-interior poly
@test inpoly([0.3,0.3], poly)
@test !inpoly([0.3,0.5], poly)

poly = Float64[0 0
               # interior #1
               0.1 0.1
               0.1 0.6
               0.4 0.6
               0.4 0.6
               0.4 0.1
               0.1 0.1
               0 0
               # interior #2
               0.6 0.4
               0.6 0.6
               0.8 0.6
               0.8 0.4
               0.6 0.4
               0 0
               # exterior
               0 1
               1 1
               1 0
               0 0]'
@test !inpoly([0.2,0.4], poly)
@test !inpoly([0.3,0.15], poly)
@test inpoly([0.5,0.4], poly)
@test inpoly([0.5,0.2], poly)
@test !inpoly([0.7,0.5], poly)

#####
# boxcar
#####
#nr,nc = 4,5
T = Float32
nr,nc = 20,31
#nr,nc = 600,500
window = 3
# orig = (1.0:nr)''*(1:nc)'
# weights = ones(Int,nr,nc)
srand(1)
orig = rand(T,nr,nc)
weights = rand(nr,nc)
weightsb = bitrand(nr,nc)
weightsbb = convert(Matrix{Bool}, weightsb)
# with weights
filt1 = BM.boxcar(orig, window, weights)
M = BM.boxcar_matrix(T, window, weights)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = reshape(M*reshape(orig, length(orig)), size(orig))
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test_approx_eq filt1[i] filt2[i]
end
# with weightsb
filt1 = BM.boxcar(orig, window, weightsb)
M = BM.boxcar_matrix(T, window, weightsb)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = reshape(M*reshape(orig, length(orig)), size(orig))
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test_approx_eq filt1[i] filt2[i]
end
# with weightsbb
filt1 = BM.boxcar(orig, window, weightsbb)
M = BM.boxcar_matrix(T, window, weightsbb)
@test size(M,1)==size(M,2)
@test size(M,1)==length(orig)
filt2 = reshape(M*reshape(orig, length(orig)), size(orig))
@test eltype(orig)==eltype(filt1)
@test eltype(orig)==eltype(filt2)
for i=eachindex(orig)
    @test_approx_eq filt1[i] filt2[i]
end

####
# piecewiselinear
####
f = BM.piecewiselinear([78.0], [-12.0])
@test f(-1)==-12
@test f(100)==-12

f = BM.piecewiselinear([0, 1, 2], [-1, 1, -2])
@test f(0)==-1
@test f(0.5)==0
@test f(1)==1
@test f(1.5)==-0.5
@test f(2)==-2


text:
function TableauSymplecticProjection(name, q::Tableau{T}, p::Tableau{T}, d=[]; R∞=1) where {T}

    @assert q.s == p.s

    o = min(q.o, p.o)

    a_q = q.a
    a_p = p.a

    α_q = zeros(T, q.s, 2)
    α_q[:,1] .= 0.5

    α_p = zeros(T, p.s, 2)
    α_p[:,1] .= 0.5

    a_q̃ = Array(transpose(hcat(zero(q.b), q.b)))
    a_p̃ = Array(transpose(hcat(zero(p.b), p.b)))

    α_q̃ = [[0.0  0.0]
           [0.5  R∞*0.5]]
    α_p̃ = [[0.0  0.0]
           [0.5  R∞*0.5]]

    b_q = q.b
    b_p = p.b
    β_q = [0.5, R∞*0.5]
    β_p = [0.5, R∞*0.5]

    c_q = q.c
    c_p = p.c
    c_λ = [0.0, 1.0]
    d_λ = [0.0, 1.0]

    if length(d) == 0
        return TableauVPARK(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ)
    else
        @assert length(d) == q.s == p.s

        return TableauVPARK(name, o,
                            a_q, a_p, α_q, α_p,
                            a_q̃, a_p̃, α_q̃, α_p̃,
                            b_q, b_p, β_q, β_p,
                            c_q, c_p, c_λ, d_λ,
                            d)
    end

end

"Tableau for Gauss-Lobatto IIIA-IIIB method with s stages and symplectic projection."
function TableauLobattoIIIAIIIBpSymplectic(s)
    TableauSymplecticProjection(Symbol("LobattoIIIAIIIB($s)pSymplectic"), TableauLobattoIIIA(s), TableauLobattoIIIB(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Lobatto IIIB-IIIA method with s stages and symplectic projection."
function TableauLobattoIIIBIIIApSymplectic(s)
    TableauSymplecticProjection(Symbol("LobattoIIIBIIIA($s)pSymplectic"), TableauLobattoIIIB(s), TableauLobattoIIIA(s), get_lobatto_nullvector(s); R∞=-1^(s+1))
end

"Tableau for Gauss-Legendre method with s stages and symplectic projection."
function TableauGausspSymplectic(s)
    glrk = TableauGauss(s)
    TableauSymplecticProjection(Symbol("GLRK($s)pSymplectic"), glrk, glrk; R∞=(-1)^s)
end


text:
const DESCRIPTORS = joinpath(RTS_GMLC_DIR, "user_descriptors.yaml")

function create_rts_system(forecast_resolution = Dates.Hour(1))
    data = PowerSystemTableData(RTS_GMLC_DIR, 100.0, DESCRIPTORS)
    return System(data; forecast_resolution = forecast_resolution)
end

"""Allows comparison of structs that were created from different parsers which causes them
to have different UUIDs."""
function compare_values_without_uuids(x::T, y::T)::Bool where {T <: PowerSystemType}
    match = true

    for (fieldname, fieldtype) in zip(fieldnames(T), fieldtypes(T))
        if fieldname == :internal
            continue
        end

        val1 = getfield(x, fieldname)
        val2 = getfield(y, fieldname)

        # Recurse if this is a PowerSystemType.
        if val1 isa PowerSystemType
            if !compare_values_without_uuids(val1, val2)
                match = false
            end
            continue
        end

        if val1 != val2
            @error "values do not match" fieldname repr(val1) repr(val2)
            match = false
        end
    end

    return match
end

"""Return the first component of type component_type that matches the name of other."""
function get_component_by_name(sys::System, component_type, other::Component)
    for component in get_components(component_type, sys)
        if get_name(component) == get_name(other)
            return component
        end
    end

    error("Did not find component $component")
end

"""Return the Branch in the system that matches another by case-insensitive arc
names."""
function get_branch(sys::System, other::Branch)
    for branch in get_components(Branch, sys)
        if lowercase(other.arc.from.name) == lowercase(branch.arc.from.name) &&
           lowercase(other.arc.to.name) == lowercase(branch.arc.to.name)
            return branch
        end
    end

    error("Did not find branch with buses $(other.arc.from.name) ", "$(other.arc.to.name)")
end


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

mutable struct CreateFirewallRuleWithAccountParameters <: SwaggerModel
    name::Any # spec type: Union{ Nothing, String } # spec name: name
    properties::Any # spec type: Union{ Nothing, CreateOrUpdateFirewallRuleProperties } # spec name: properties

    function CreateFirewallRuleWithAccountParameters(;name=nothing, properties=nothing)
        o = new()
        validate_property(CreateFirewallRuleWithAccountParameters, Symbol("name"), name)
        setfield!(o, Symbol("name"), name)
        validate_property(CreateFirewallRuleWithAccountParameters, Symbol("properties"), properties)
        setfield!(o, Symbol("properties"), properties)
        o
    end
end # type CreateFirewallRuleWithAccountParameters

const _property_map_CreateFirewallRuleWithAccountParameters = Dict{Symbol,Symbol}(Symbol("name")=>Symbol("name"), Symbol("properties")=>Symbol("properties"))
const _property_types_CreateFirewallRuleWithAccountParameters = Dict{Symbol,String}(Symbol("name")=>"String", Symbol("properties")=>"CreateOrUpdateFirewallRuleProperties")
Base.propertynames(::Type{ CreateFirewallRuleWithAccountParameters }) = collect(keys(_property_map_CreateFirewallRuleWithAccountParameters))
Swagger.property_type(::Type{ CreateFirewallRuleWithAccountParameters }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_CreateFirewallRuleWithAccountParameters[name]))}
Swagger.field_name(::Type{ CreateFirewallRuleWithAccountParameters }, property_name::Symbol) =  _property_map_CreateFirewallRuleWithAccountParameters[property_name]

function check_required(o::CreateFirewallRuleWithAccountParameters)
    (getproperty(o, Symbol("name")) === nothing) && (return false)
    (getproperty(o, Symbol("properties")) === nothing) && (return false)
    true
end

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


text:
module TestArithmetic

using Test
using CategoricalDistributions
using StableRNGs
rng = StableRNG(123)

L = ["yes", "no"]
d1 = UnivariateFinite(L, rand(rng, 2), pool=missing)
d2 = UnivariateFinite(L, rand(rng, 2), pool=missing)

@testset "arithmetic" begin

    # addition and subtraction:
    for op in [:+, :-]
        quote
            s = $op(d1, d2 )
            @test $op(pdf.(d1, L), pdf.(d2, L)) ≈ pdf.(s, L)
        end |> eval
    end

    # negative:
    d_neg = -d1
    @test pdf.(d_neg, L) == -pdf.(d1, L)

    # multiplication by scalar:
    d3 = d1*42
    @test pdf.(d3, L) ≈ pdf.(d1, L)*42
    d3 = 42*d1
    @test pdf.(d3, L) ≈ pdf.(d1, L)*42

    # division by scalar:
    d3 = d1/42
    @test pdf.(d3, L) ≈ pdf.(d1, L)/42

    # "probabilities" that aren't:
    d = UnivariateFinite(L, randn(rng, 2) + im*randn(2), pool=missing)
    @test pdf.(3*d -(4*d)/2, L) ≈ pdf.(d, L)
end

p = [0.1, 0.9]
P = vcat(fill(p', 10^5)...);
slow = fill(UnivariateFinite(L, p, pool=missing), 10^5);
fast = UnivariateFinite(L, P, pool=missing);
# @assert pdf(slow, L) == pdf(fast, L)

@testset "performant arithmetic for UnivariateFiniteArray" begin
    @test pdf(slow + slow, L) == pdf(fast + fast, L)
    t_slow = @elapsed @eval slow + slow
    t_fast = @elapsed @eval fast + fast
    @test t_slow/t_fast > 10

    @test pdf(slow - slow, L) == pdf(fast - fast, L)
    t_slow = @elapsed @eval slow - slow
    t_fast = @elapsed @eval fast - fast
    @test t_slow/t_fast > 10

    @test pdf(42*slow, L) == pdf(42*fast, L)
    @test pdf(slow*42, L) == pdf(fast*42, L)
    t_slow = @elapsed @eval 42*slow
    t_fast = @elapsed @eval 42*fast
    @test t_slow/t_fast > 10
    t_slow = @elapsed @eval slow*42
    t_fast = @elapsed @eval fast*42
    @test t_slow/t_fast > 10

    @test pdf(slow/42, L) == pdf(fast/42, L)
    t_slow = @elapsed @eval slow/42
    t_fast = @elapsed @eval fast/42
    @test t_slow/t_fast > 10
end

end # module

true


text:
function nlp_cvx_105_013(
    optimizer,
    objective_tol,
    primal_tol,
    dual_tol,
    termination_target = TERMINATION_TARGET_LOCAL,
    primal_target = PRIMAL_TARGET_LOCAL,
)
    model = Model(optimizer)

    @variable(model, x, start = 0.1)
    @variable(model, y)

    @objective(model, Min, -x + y)
    @NLconstraint(model, exp(x - 2.0) - 0.5 <= y)
    @NLconstraint(model, log(x) + 0.5 >= y)

    optimize!(model)

    check_status(model, FEASIBLE_PROBLEM, termination_target, primal_target)
    check_objective(model, -3 / 2, tol = objective_tol)
    return check_solution([x, y], [2, 1 / 2], tol = primal_tol)
end