text:
ENV["MPLBACKEND"]="Qt5Agg"
include("CreateMesh.jl")
using LoopVectorization
using PyPlot
pygui(true)

function tplot(p, t, v)
    # Plot triangular mesh with nodes `p` and triangles `t`
    clf()
    tris = convert(Array{Int64}, hcat(t...)')
    # fig1, ax1 = plt.subplots()
    display(tripcolor(first.(p), last.(p), tris .- 1, v,
              cmap="viridis", edgecolors="none", linewidth=0))
    axis("equal")
    ylim([0, 1])
    xlim([0, 1])
    title("Low res.")
    xlabel("x")
    ylabel("y")
    colorbar()
    show()
    return
end

function qplot(x, y, v)
    clf()
    display( pcolor(x, y, v) )
    colorbar()
    xlabel("x")
    ylabel("y")
    show()
end

function Tanalytic2!( mesh, xc , yc, T, a, b, c, d, alp, bet )
    if mesh.type == "triangle" nn_el = 3.0; nn_el_int=3; end
    if mesh.type == "quad"     nn_el = 4.0; nn_el_int=4; end
    # Evaluate T analytic on barycentres
    for j=1:mesh.nn_el # as fast as above and less cryptic
        @. xc = 0.0
        @. yc = 0.0
        @. xc = xc + 1.0/mesh.nn_el * mesh.xn[mesh.e2n[j,:]]
        @. yc = yc + 1.0/mesh.nn_el * mesh.yn[mesh.e2n[j,:]]
    end
    @. T = exp(alp*sin(a*xc + c*yc) + bet*cos(b*xc + d*yc));
    # xn = mesh.xn
    # yn = mesh.yn
    # e2n = mesh.e2n
    #  for i=1:mesh.nel
    #     xc[i] = 0;
    #     yc[i] = 0;
    #     @avx for j=1:3
    #         # xc[i] += 1.0/nn_el * mesh.xn[mesh.e2n[j,i]]
    #         # yc[i] += 1.0/nn_el * mesh.yn[mesh.e2n[j,i]]
    #         xc[i] += 1.0/nn_el * xn[e2n[j,i]]
    #         yc[i] += 1.0/nn_el * yn[e2n[j,i]]
    #     end
    # end
    # for i=1:mesh.nel
    #     T[i] = exp(alp.sin(a*xc[i] + c*yc[i]) + bet*cos(b*xc[i] + d*yc[i]));
    # end
return
end

function main()

    # Create sides of mesh
    xmin, xmax = 0, 1
    ymin, ymax = 0, 1
    nx, ny     = 20, 20
    quad = false

    if quad==true
        mesh = MakeQuadMesh( nx, ny, xmin, xmax, ymin, ymax )
    elseif quad==false
        mesh = MakeTriangleMesh( nx, ny, xmin, xmax, ymin, ymax )
    end

    println(mesh.nel)

    # Generate function to be visualised on mesh
    alp = 0.1; bet = 0.3; a = 5.1; b = 4.3; c = -6.2; d = 3.4;
    xc = zeros(mesh.nel)
    yc = zeros(mesh.nel)
    T  = zeros(mesh.nel)
    ncalls = 4
    # A) Loop version with @avx
    # @printf("Looped, %d times:\n", ncalls)
    for icall=1:ncalls
        @time Tanalytic2!(mesh, xc ,yc , T, a, b, c, d, alp, bet)
    end

    if mesh.type == "triangle"  tplot(mesh.xyviz, mesh.e2nviz,                    T ) end
    if mesh.type == "quad"      qplot(mesh.xcviz, mesh.ycviz,  reshape(T, (nx,ny) ) ) end

end

main()


text:
using DynamicHMCModels, MCMCChains, Random

#Random.seed!(1233)

ProjDir = @__DIR__
cd(ProjDir)

include(joinpath(@__DIR__, "LBA_functions.jl"))

Base.@kwdef struct LBAModel{T}
   data::T
   N::Int
   Nc::Int
end

function make_transformation(model::LBAModel)
    as((v=as(Array,asℝ₊,Nc), A=asℝ₊, k=asℝ₊, tau=asℝ₊))
end

N = 10
v = [1.0, 1.5]
Nc = length(v)
data=simulateLBA(;Nd=N,v=v,A=.8,k=.2,tau=.4)

#dist = LBA(ν=[1.0,1.5,2.0],A=.8,k=.2,τ=.4)
#data = rand(dist,N)

model = LBAModel(; data=data, N=N, Nc=Nc)

function (model::LBAModel)(θ)
    @unpack data=model
    @unpack v,A,k,tau=θ
    d=LBA(ν=v,A=A,k=k,τ=tau)
    minRT = minimum(x->x[2],data)
    logpdf(d,data)+sum(logpdf.(TruncatedNormal(0,3,0,Inf),v)) +
    logpdf(TruncatedNormal(.8,.4,0,Inf),A)+logpdf(TruncatedNormal(.2,.3,0,Inf),k)+
    logpdf(TruncatedNormal(.4,.1,0,minRT),tau)
end

d = [(c,r) for (c,r) in zip(data.choice,data.rt)]

p = LBAModel(d,N,Nc)
p((v=fill(.5,Nc),A=.8,k=.2,tau=.4))

# Use Flux for the gradient.
P = TransformedLogDensity(make_transformation(p), p)
∇P = ADgradient(:ForwardDiff, P)

# Sample from the posterior.
results = mcmc_with_warmup(Random.GLOBAL_RNG, ∇P, 1000;
  warmup_stages = default_warmup_stages(local_optimization=nothing),
  reporter = NoProgressReport()
  )
posterior = P.transformation.(results.chain)

@show DynamicHMC.Diagnostics.EBFMI(results.tree_statistics)
println()
@show DynamicHMC.Diagnostics.summarize_tree_statistics(results.tree_statistics)
println()

parameter_names = ["v[1]", "v[2]", "A", "k", "tau"]

# Create a3d
a3d = Array{Float64, 3}(undef, 1000, 5, 1);
for j in 1:1
  for i in 1:1000
    a3d[i, 1:2, j] = values(posterior[i].v)
    a3d[i, 3, j] = values(posterior[i].A)
    a3d[i, 4, j] = values(posterior[i].k)
    a3d[i, 5, j] = values(posterior[i].tau)
  end
end

chns = MCMCChains.Chains(a3d,
  vcat(parameter_names),
  Dict(
    :parameters => parameter_names
  )
)

describe(chns)

ext: jl


text:
using DiffEqFlux, OrdinaryDiffEq, Optim, Flux, Zygote, Test

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))

fastdudt2 = FastChain((x,p) -> x.^3,
             FastDense(2,50,tanh),
             FastDense(50,2))
fast_n_ode = NeuralODE(fastdudt2,tspan,Tsit5(),saveat=t)

function fast_predict_n_ode(p)
  fast_n_ode(u0,p)
end

function fast_loss_n_ode(p)
    pred = fast_predict_n_ode(p)
    loss = sum(abs2,ode_data .- pred)
    loss,pred
end

staticdudt2 = FastChain((x,p) -> x.^3,
                        StaticDense(2,50,tanh),
                        StaticDense(50,2))
static_n_ode = NeuralODE(staticdudt2,tspan,Tsit5(),saveat=t)

function static_predict_n_ode(p)
  static_n_ode(u0,p)
end

function static_loss_n_ode(p)
    pred = static_predict_n_ode(p)
    loss = sum(abs2,ode_data .- pred)
    loss,pred
end

dudt2 = Chain((x) -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
n_ode = NeuralODE(dudt2,tspan,Tsit5(),saveat=t)

function predict_n_ode(p)
  n_ode(u0,p)
end

function loss_n_ode(p)
    pred = predict_n_ode(p)
    loss = sum(abs2,ode_data .- pred)
    loss,pred
end

p = initial_params(fastdudt2)
_p,re = Flux.destructure(dudt2)
@test fastdudt2(ones(2),_p) ≈ dudt2(ones(2))
@test staticdudt2(ones(2),_p) ≈ dudt2(ones(2))
@test fast_loss_n_ode(p)[1] ≈ loss_n_ode(p)[1]
@test static_loss_n_ode(p)[1] ≈ loss_n_ode(p)[1]
@test Zygote.gradient((p)->fast_loss_n_ode(p)[1], p)[1] ≈ Zygote.gradient((p)->loss_n_ode(p)[1], p)[1]
@test Zygote.gradient((p)->static_loss_n_ode(p)[1], p)[1] ≈ Zygote.gradient((p)->loss_n_ode(p)[1], p)[1]

#=
using BenchmarkTools
@btime Zygote.gradient((p)->static_loss_n_ode(p)[1], p)
@btime Zygote.gradient((p)->fast_loss_n_ode(p)[1], p)
@btime Zygote.gradient((p)->loss_n_ode(p)[1], p)
=#


text:
import .AbstractPaulis: AbstractPaulis, AbstractPauli, is_pauli_y
import ._op_term_macro_helper
import LightGraphs
import .Utils: property_graph, kron_alt, isapprox_zero, triprod, pow_of_minus_one

####
#### Constructors
####

const PauliTerm = OpTerm{Pauli}
const DensePauliTerm = DenseOpTerm{<:AbstractPauli}
const PauliSum = OpSum{Pauli}

## This ruins printing PauliTerm as an alias if `APauliTerm` is in the export list
"""
    APauliTerm

A `UnionAll` type representing a term of any type `<:AbstractPauli`.
"""
const APauliTerm = OpTerm{T} where {T <:AbstractPauli}

"""
    APauliSum

A `UnionAll` type representing a sum of terms of factors of any type `<:AbstractPauli`.
"""
const APauliSum = OpSum{T} where {T <:AbstractPauli}

# There is a printing bug in Julia that makes this a bad option
#const PauliTerm = APauliTerm{Pauli}
"""
    @pauli_str(str::String) -> PauliTerm

Construct a `PauliTerm` from `str`.

# Examples
```jldoctest
julia> pauli"XXY"
3-factor PauliTerm{Vector{Pauli}, Complex{Int64}}:
XXY * (1 + 0im)

julia> pauli"3 * XXY"
3-factor PauliTerm{Vector{Pauli}, ComplexF64}:
XXY * (3.0 + 0.0im)

julia> pauli"3.1 + 2.0im * XXY"
3-factor PauliTerm{Vector{Pauli}, ComplexF64}:
XXY * (3.1 + 2.0im)
```
"""
macro pauli_str(str)
    return _op_term_macro_helper(PauliTerm, str)
end

"""
    PauliSumA(::Type{PauliT}, matrix::AbstractMatrix{<:Number}; threads=true)

Construct a Pauli decomposition of `matrix`, that is, a `PauliSumA` representing `matrix`.
If `thread` is `true`, use a multi-threaded algorithm for increased performance.
"""
function OpSum{PauliT}(matrix::AbstractMatrix{<:Number}; threads::Bool=true) where PauliT <: AbstractPauli
    if threads
        return pauli_sum_from_matrix_threaded(PauliT, matrix)
    else
        return pauli_sum_from_matrix_one_thread(PauliT, matrix)
    end
end

function pauli_sum_from_matrix_one_thread(::Type{PauliT}, matrix::AbstractMatrix{<:Number}) where PauliT
    nside = LinearAlgebra.checksquare(matrix)
    n_qubits = ILog2.checkispow2(nside)
    denom = 2^n_qubits  # == nside
    s = OpSum{PauliT}()
    for pauli in pauli_basis(PauliT, n_qubits)
        mp = SparseArrays.sparse(pauli)  # Much faster than dense
        coeff = LinearAlgebra.dot(mp, matrix)
        if ! isapprox_zero(coeff)
            push!(s, (op_string(pauli), coeff / denom))  # a bit faster than PauliTermA for small `matrix` (eg 2x2)
        end
    end
    return s
end

function pauli_sum_from_matrix_threaded(::Type{PauliT}, matrix::AbstractMatrix{<:Number}) where PauliT
    nside = LinearAlgebra.checksquare(matrix)
    n_qubits = ILog2.checkispow2(nside)
    denom = 2^n_qubits  # == nside
    ## Create a PauliSumA for each thread, for accumulation.
    sums = [OpSum{PauliT}() for i in 1:Threads.nthreads()]
    Threads.@threads for j in 0:(4^n_qubits - 1)
        pauli = OpTerm{PauliT}(j, n_qubits)
        mp = SparseArrays.sparse(pauli)  # Much faster than dense
        coeff = LinearAlgebra.dot(mp, matrix)
        if ! isapprox_zero(coeff)
            push!(sums[Threads.threadid()], (op_string(pauli), coeff / denom))
        end
    end
    for ind in 2:length(sums)  # Collate results from all threads.
        add!(sums[1], sums[ind])
    end
    return sums[1]
end

function OpTerm{T}(index::Integer, n_paulis::Integer, coeff=_DEFAULT_COEFF) where T <: AbstractPauli
    return OpTerm(AbstractPaulis.pauli_vector(T, index, n_paulis), coeff)
end

####
####
####

function Base.zero(pterm::OpTerm{T}) where T <: AbstractPauli
    new_string = similar(op_string(pterm))
    fill!(new_string, one(T))
    return strip_typeof(pterm)(new_string, zero(pterm.coeff))
end

####
#### Math
####

"""
    cis(p::PauliTerm)::PauliSum

Compute ``\\exp(i p)``
"""
function Base.cis(pt::OpTerm{<:AbstractPauli})
    return OpSum([op_string(one(pt)), copy(op_string(pt))], [cos(pt.coeff), im * sin(pt.coeff)], true)
end

"""
    exp(p::PauliTerm)::PauliSum

Compute ``\\exp(p)``
"""
function Base.exp(pt::OpTerm{<:AbstractPauli})
    coeff = im * pt.coeff
    return OpSum([op_string(one(pt)), copy(op_string(pt))], [cos(coeff), -im * sin(coeff)], true)
end

"""
    numeric_function(pt::PauliTerm, f)::PauliSum

Compute `f(pt)` by decomposing `f` into odd and even functions.
"""
function numeric_function(pt::OpTerm{<:AbstractPauli}, f)
    c = pt.coeff
    fe = (f(c) + f(-c)) / 2  # Even term
    fo = (f(c) - f(-c)) / 2  # Odd term
    if isapprox_zero(fe)
        strings = [copy(op_string(pt))]
        coeffs = [fo]
    elseif isapprox_zero(fo)
        strings = [op_string(one(pt))]
        coeffs = [fe]
    else
        strings = [op_string(one(pt)), copy(op_string(pt))]
        coeffs = [fe, fo]
    end
    # sorting would take 30x longer
    return OpSum(strings, coeffs; already_sorted=true)
end

# Julia 1.5 does not have cispi
for f in (:cos, :sin, :tan, :sqrt, :sind, :sinpi, :cospi, :sinh, :tanh,
          :acos, :asin, :atan, :sec, :csc, :cot, :log, :log2, :log10,
          :log1p)
    @eval Base.$f(pt::OpTerm{<:AbstractPauli}) = numeric_function(pt, $f)
end

####
#### Conversion
####

function _multiply_coefficient(coeff, matrix)
    if isbool_and_equal(coeff, one(eltype(matrix)))
        return matrix
    end
    # Tests pass with the stuff below commented out. It does not seem useful.
    # if coeff isa eltype(matrix)
    #     return LinearAlgebra.lmul!(coeff, copy(matrix))
    # end
    # coeff1 = convert(promote_type(typeof(coeff), eltype(matrix)), coeff)
    # return coeff1 .* matrix
    return coeff .* matrix
end

Base.Matrix(pt::APauliTerm) = Matrix(Float64, pt)

# FIXME: Not type stable.
function Base.Matrix(::Type{Z4Group0}, pt::APauliTerm)
    matrix = kron_alt((Z4Group0.(m) for m in Matrix.(op_string(pt)))...)
    return _multiply_coefficient(pt.coeff, matrix)
end

function Base.Matrix(::Type{Float64}, pt::APauliTerm)
    matrix = kron_alt(Matrix.(op_string(pt))...)
    return _multiply_coefficient(pt.coeff, matrix)
end

SparseArrays.sparse(pt::APauliTerm) = SparseArrays.sparse(Float64, pt)

function SparseArrays.sparse(::Type{Float64}, pt::APauliTerm)
    matrix = kron_alt(SparseArrays.sparse.(op_string(pt))...)
    return _multiply_coefficient(pt.coeff, matrix)
end

# FIXME: broken, should throw inexact error earlier rather than return wrong type
function SparseArrays.sparse(::Type{Z4Group0}, pt::APauliTerm)
    matrix = kron_alt((Z4Group0.(m) for m in SparseArrays.sparse.(op_string(pt)))...)
    return _multiply_coefficient(pt.coeff, matrix)
end

####
#### Compare / predicates
####

"""
    isunitary(pt::PauliTermA)

Return `true` if `pt` is a unitary operator.
"""
IsApprox.isunitary(pt::APauliTerm, approx_test::AbstractApprox=Equal()) = IsApprox.isunitary(pt.coeff, approx_test)

"""
    IsApprox.ishermitian(pt::PauliTermA)

Return `true` if `pt` is a Hermitian operator.
"""
IsApprox.ishermitian(pt::APauliTerm, approx_test::AbstractApprox=Equal()) =
    IsApprox.isreal(pt.coeff, approx_test)

## TODO: anticommutes. for FermiOp as well

@inline IsApprox.commutes(p1::APauliTerm, p2::APauliTerm) = commutes(op_string(p1), op_string(p2))

#using LoopVectorization
## Each non-commuting factor introduces a phase factor of -1.
#  This is slower   iseven(count(x -> !commutes(x[1], x[2]), zip(v1, v2)))
@inline function IsApprox.commutes(v1::AbstractArray{<:AbstractPauli}, v2::AbstractArray{<:AbstractPauli})
    length(v1) == length(v2) || throw(DimensionMismatch("v1 and v2 of differing lengths"))
    c = 0
    @inbounds for i in eachindex(v1)
        if !commutes(v1[i], v2[i])
            c += 1
        end
    end
    return iseven(c)
end

"""
    commutes(pauli_sum::APauliSum)

Return `true` if the terms in `pauli_sum` are mutually commuting.
"""
IsApprox.commutes(ps::APauliSum) = commutes(op_strings(ps))

## The following is simpler, but slightly less performant than `commutes` below
## triprod could be more efficient
function acommutes(a::AbstractVector{<:AbstractVector{<:AbstractPauli}})
   return all(x -> commutes(x...), triprod(a))
end

"""
    commutes(pauli_strings::AbstractVector{<:AbstractVector{<:AbstractPauli}})

Return `true` if the strings in `pauli_strings` are mutually commuting.
"""
function IsApprox.commutes(a::AbstractVector{<:AbstractVector{<:AbstractPauli}})
    for i in eachindex(a), j in i+1:lastindex(a)
        @inbounds commutes(a[i], a[j]) || return false
    end
    return true
end

####
#### Algebra
####

function Base.conj(pt::APauliTerm)
    return Base.transpose(pt, conj)
end

function Base.transpose(pt::APauliTerm, conj_func=identity)
    num_ys = count(is_pauli_y, op_string(pt))
    fac = pow_of_minus_one(num_ys)
    return strip_typeof(pt)(op_string(pt), conj_func(pt.coeff) * fac)
end

function Base.adjoint(pt::APauliTerm)
    return strip_typeof(pt)(op_string(pt), conj(pt.coeff))
end

Base.inv(p::APauliTerm) = strip_typeof(p)(op_string(p), inv(p.coeff))

function Base.:^(p::APauliTerm, n::Integer)
    new_coeff = p.coeff^n
    if iseven(n)
        return strip_typeof(p)(fill(Pauli(:I), length(p)), new_coeff)
    else
        return strip_typeof(p)(op_string(p), new_coeff)
    end
end

Base.kron(ps1::APauliTerm, ps2::APauliTerm) = strip_typeof(ps1)(vcat(op_string(ps1), op_string(ps2)), ps1.coeff * ps2.coeff)

Base.kron(paulis::AbstractPauli...) = PauliTerm([paulis...])

## TODO: PauliTerm is hardcoded. Should be OpTerm{T} where T is some AbstractPauli
# TODO: @code_warntype shows red here
function Base.kron(ps::Union{APauliTerm, AbstractPauli}...)
    if ps[1] isa AbstractPauli
        v = typeof(ps[1])[]
    else
        v = eltype(ps[1])[]
    end
    coeffs = []
    for x in ps
        if x isa APauliTerm
            push!(coeffs, x.coeff)
            append!(v, x)
        else
            push!(v, x)
        end
    end
    if isempty(coeffs)
        return PauliTerm(v)
    elseif length(coeffs) == 1
        return PauliTerm(v, only(coeffs))
    else
        return PauliTerm(v, reduce(*, coeffs))
    end
end

function LinearAlgebra.eigvals(pt::APauliTerm)
    if all(isone, op_string(pt)) # pt is propto identity
        return fill(float(pt.coeff), 2^length(pt)) # All eigenvalues are one
    end
    vals = Vector{promote_type(Float64, typeof(pt.coeff))}(undef, 2^length(pt))
    pos_eigval = pt.coeff
    neg_eigval = -pos_eigval
    if real(neg_eigval) > real(pos_eigval)
        (pos_eigval, neg_eigval) = (neg_eigval, pos_eigval)
    end
    half_length = 2^(length(pt)-1)
    @inbounds for i in 1:half_length
        vals[i] = neg_eigval  # follow the usual order
        vals[i+half_length] = pos_eigval
    end
    return vals
end

####
#### Other ?
####

"""
    pauli_basis(::Type{PauliT}, n_qubits; coeff=_DEFAULT_COEFF)

Return a `Generator` over all `PauliTerm`s of `n_qubits`.
"""
function pauli_basis(::Type{PauliT}, n_qubits; coeff=_DEFAULT_COEFF) where PauliT
    return (OpTerm{PauliT}(i, n_qubits, coeff) for i in 0:(4^n_qubits - 1))
end

pauli_basis(n_qubits; coeff=_DEFAULT_COEFF) = pauli_basis(PauliDefault, n_qubits; coeff=coeff)

"""
    group_paulis(ps::APauliSum; sort_degree=false, reps=1)
    group_paulis(strings::AbstractVector{T}; sort_degree=false, reps=1) where {T<:AbstractVector{<:AbstractPauli}}

Return a partition of `ps` or `strings` into groups of mutually commuting strings or terms.
The return type is `Vector{typeof(ps)}` or `Vector{typeof(strings)}`.
`sort_degree` and `reps` are passed to `LightGraphs.greedy_color`.

# Examples
```julia
julia> ps = rand_op_sum(Pauli, 5, 5; coeff_func=randn)
5x5 PauliSum{Vector{Vector{Pauli}}, Vector{Float64}}:
IXXZZ * 0.7287656446540858
IXZIZ * 0.6497331352636981
YIZZX * -0.7581101586820902
YZIXI * -1.9067637258086527
ZYIII * 1.3523639501623812

julia> group_paulis(ps)
4-element Vector{PauliSum{Vector{Vector{Pauli}}, Vector{Float64}}}:
 2x5 PauliSum{Vector{Vector{Pauli}}, Vector{Float64}}:
IXXZZ * 0.7287656446540858
YZIXI * -1.9067637258086527
 1x5 PauliSum{Vector{Vector{Pauli}}, Vector{Float64}}:
IXZIZ * 0.6497331352636981
 1x5 PauliSum{Vector{Vector{Pauli}}, Vector{Float64}}:
YIZZX * -0.7581101586820902
 1x5 PauliSum{Vector{Vector{Pauli}}, Vector{Float64}}:
ZYIII * 1.3523639501623812

julia> all(commutes, group_paulis(ps))
true

julia> group_paulis(ps.strings)
4-element Vector{Vector{Vector{Pauli}}}:
 [IXXZZ, YZIXI]
 [IXZIZ]
 [ZYIII]
 [YIZZX]
```
"""
group_paulis(ps::APauliSum; sort_degree=false, reps=1) = _group_paulis(ps, ps.strings, sort_degree, reps)

group_paulis(strings::AbstractVector{T}; sort_degree=false, reps=1) where {T<:AbstractVector{<:AbstractPauli}} =
    _group_paulis(strings, strings, sort_degree, reps)

"""
    _group_paulis(_sum, strings)

Partition indices of `strings` into groups corresponding to mutually commuting terms
in `strings`. Return a partition of `_sum`, which must have the same indices as `strings`
according to partitioned indices. Note that `group_paulis(strings, strings)` simply
partitions `strings`.
"""
function _group_paulis(_sum, strings, sort_degree, reps)
    graph = property_graph(strings, !commutes) # edges for non-commuting pairs
    coloring = LightGraphs.greedy_color(graph; sort_degree=sort_degree, reps=reps)
    groups = [empty(_sum) for i in 1:coloring.num_colors]
    for (i, c) in enumerate(coloring.colors)
        push!(groups[c], _sum[i])
    end
    return groups
end

"""
    ishermitian(ps::PauliSumA)

Return `true` if `ps` is a Hermitian operator.
"""
IsApprox.ishermitian(ps::APauliSum, approx_test::AbstractApprox=Equal()) =
    all(x -> IsApprox.isreal(x, approx_test), ps.coeffs)


text:
module TreeUtil
    using ..PyUtils
    using ..Core

    export flatten, unflatten

    const _tree_util = PyNULL()

    function tree_flatten(tree)
        vals, f = _tree_util.tree_flatten(tree)
        vals = vals isa AbstractArray ? tuple(vals...) : vals
        return vals, f
    end

    function _flatten end
    function _unflatten end

    jl_flatten(x) = _flatten(x)
    jl_unflatten(f, x) = _unflatten(f, x)

    const _fwdfuncs = Dict(
        :Partial=>:Partial,
        :all_leaves=>:all_leaves,
        :build_tree=>:build_tree,
        :register_pytree_node=>:register_pytree_node,
        :tree_all=>:tree_all,
        :tree_leaves=>:tree_leaves,
        :tree_map=>:tree_map,
        :tree_multimap=>:tree_multimap,
        :tree_reduce=>:tree_reduce,
        :tree_structure=>:tree_structure,
        :tree_transpose=>:tree_transpose,
        :tree_unflatten=>:tree_unflatten,
        :treedef_children=>:treedef_children,
        :treedef_is_leaf=>:treedef_is_leaf,
        :treedef_tuple=>:treedef_tuple)
    for (jl,py) in _fwdfuncs
        @eval begin
            $jl(args...; kwargs...) = _tree_util.$py(args...;kwargs...)
        end
    end

    const flatten = tree_flatten
    const unflatten = tree_unflatten

    function __init__()
        copy!(_tree_util, jax.tree_util)

        # This is an horrinble hack.
        # it will register as a pytree node any wrapper of julia
        # objects.
        # It woudl be nice to internally dispatch on traits when calling
        TreeUtil.register_pytree_node(
            pytypeof(PyObject(identity)),
            jl_flatten,
            jl_unflatten,
        )
    end
end


text:
module DebugRequest

import ..DEBUG_LEVEL
using ..IOExtras
import ..Streams: Stream

include("IODebug.jl")

export debuglayer

"""
    debuglayer(stream::Stream) -> HTTP.Response

Wrap the `IO` stream in an `IODebug` stream and print Message data.
"""
function debuglayer(handler)
    return function(stream::Stream; verbose::Int=0, kw...)
        # if debugging, wrap stream.stream in IODebug
        if verbose >= 3 || DEBUG_LEVEL[] >= 3
            stream = Stream(stream.message, IODebug(stream.stream))
        end
        return handler(stream; verbose=verbose, kw...)
    end
end

end # module DebugRequest


text:
using Colors
using Images
using Boltzmann
using MNIST
using ImageView
using Devectorize

function run_mnist()
    HiddenUnits = 100;
    Epochs = 20;

    X, y = testdata()  # test data is smaller, no need to downsample
    X = X ./ (maximum(X) - minimum(X))

    Nfeatures = size(X,1)
    Nsamples = size(X,2)

    # Mean subtraction
    mu = mean(X,2)
    X = broadcast(-,X,mu)

    # Scaling to unit variance
    Xvars = sqrt(sum(X.^2,1))
    X = sqrt(Nfeatures) * broadcast(/,X,Xvars)

    m = GRBM(28*28, HiddenUnits, (28,28); momentum=0.0, dataset=X)

    # Attempt without weight decay
    info("Running With Gaussian Visible Units")
    fit(m, X;n_iter=Epochs,lr=0.000005,persistent=false,weight_decay="l2",decay_magnitude=0.05)
    chart_weights(m.W, (28, 28); annotation="Gaussian Visible")
    chart_weights_distribution(m.W;filename="grbm_distribution.pdf",bincount=200)

    return m
end

run_mnist()

println("Press RETURN when ready")
readline(STDIN)


text:
module AstroTime

using EarthOrientation
using Reexport

import Dates
import MacroTools

export @timescale

include("TimeScales.jl")
include("Periods.jl")
include("AstroDates.jl")
include("Epochs.jl")

@reexport using .TimeScales
@reexport using .Periods
@reexport using .AstroDates
@reexport using .Epochs

import .Epochs: format

function __init__()
    Dates.CONVERSION_SPECIFIERS['t'] = TimeScale
    Dates.CONVERSION_SPECIFIERS['D'] = Epochs.DayOfYearToken
    Dates.CONVERSION_DEFAULTS[TimeScale] = TimeScales.NotATimeScale()
    Dates.CONVERSION_DEFAULTS[Epochs.DayOfYearToken] = Int64(0)

    Dates.CONVERSION_TRANSLATIONS[Epoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
        TimeScale,
    )

    Dates.CONVERSION_TRANSLATIONS[TAIEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[UTCEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[UT1Epoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TTEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TCGEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TCBEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )

    Dates.CONVERSION_TRANSLATIONS[TDBEpoch] = (
        Dates.Year,
        Dates.Month,
        Dates.Day,
        Epochs.DayOfYearToken,
        Dates.Hour,
        Dates.Minute,
        Dates.Second,
        Dates.Millisecond,
    )
end

"""
    @timescale scale epoch [args...] body

Define a new time scale, the corresponding `Epoch` type alias, and a function that calculates
the offset from TAI for the new time scale.

### Arguments ###

- `scale`: The name of the time scale
- `epoch`: The name of the `Epoch` parameter that is passed to the `tai_offset` function
- `args`: Optional additional parameters that are passed to the `tai_offset` function
- `body`: The body of the `tai_offset` function

# Example

```jldoctest
julia> @timescale GMT ep tai_offset(UTC, ep)

julia> GMT isa TimeScale
true

julia> GMTEpoch
Epoch{GMT,T} where T
```
"""
macro timescale(scale::Symbol, ep::Symbol, args...)
    scale_type = Symbol(scale, "Type")
    epoch_type = Symbol(scale, "Epoch")
    name = String(scale)

    arglist = [
        (nothing, scale_type, false, nothing),
        (ep, Epoch, false, nothing),
    ]
    body = args[end]
    append!(arglist, MacroTools.splitarg.(args[1:end-1]))
    args = map(x->MacroTools.combinearg(x...), arglist)
    def = Dict{Symbol, Any}(
        :name => :((AstroTime.Epochs).tai_offset),
        :args => map(x->MacroTools.combinearg(x...), arglist),
        :body => body,
        :kwargs => Any[],
        :whereparams => (),
    )
    func = esc(MacroTools.combinedef(def))

    MacroTools.@esc scale scale_type epoch_type

    MacroTools.@q begin
        struct $scale_type <: TimeScale end
        const $scale = $scale_type()
        const $epoch_type = Epoch{$scale}
        Base.show(io::IO, ::$scale_type) = print(io, $name)
        AstroTime.TimeScales.tryparse(::Val{Symbol($name)}) = $scale

        Dates.CONVERSION_TRANSLATIONS[$epoch_type] = (
            Dates.Year,
            Dates.Month,
            Dates.Day,
            Dates.Hour,
            Dates.Minute,
            Dates.Second,
            Dates.Millisecond,
        )
        Dates.default_format(::Type{$epoch_type}) = Dates.ISODateTimeFormat

        $func

        nothing
    end
end

function update()
    EarthOrientation.update()
    nothing
end

end # module


text:
using VortexDistributions, LaTeXStrings

include("../util.jl")

benchmarks = benchmark2D(200, 4, 12)

using Plots

naive = [mean(benchmarks[i][1].times) for i in 1:length(benchmarks)]
optim = [mean(benchmarks[i][2].times) for i in 1:length(benchmarks)]
optimInterp = [mean(benchmarks[i][3].times) for i in 1:length(benchmarks)]
n = [benchmarks[i][4] for i in 1:length(benchmarks)]

benchPlot = plot(log.(2, n), log.(naive), label="naive", legend=:bottomright)
plot!(log.(2, n), log.(optim), label="optimised")
plot!(log.(2, n), log.(optimInterp), label="optimisedInterpolated")
xlabel!(L"log_2(N)")
ylabel!(L"ln(10^{-9})s")
title!("Benchmark comparison of 2D Plaquette \n detection 200 vorticies")

savefig(benchPlot, "results/benchmarkPlotMultiVort200.png")


text:
using DataStructures

function floodfill(snakes::Vector{Snake}, width::Int, height::Int, player::Int)
    snek = snakes[player]
    opponent = (player % 2) + 1
    opponent_snake = snakes[opponent]
    fields = zeros(2)
    fields[player] = length(snek.body) - 1
    fields[opponent] = length(opponent_snake.body) - 1

    visited = Set(get_body_no_head(snek))
    for b in get_body_no_head(opponent_snake)
        push!(visited, b)
    end

    stack = Tuple{Position, Int}[]
    push!(stack, (get_head(snek), player))
    push!(stack, (get_head(opponent_snake), opponent))
    while !isempty(stack)
        pos, p = popfirst!(stack)
        if !(pos in visited)
            push!(visited, pos)
            fields[p] += 1
            for d in (up, down, left, right)
                next_pos = move_from(pos, d)
                if !(next_pos in visited) && check_bounds(next_pos, height, width)
                    push!(stack, (next_pos, p))
                end
            end
        end
    end

    return fields[player] - fields[opponent]
end

function maxN_floodfill(snakes::Vector{Snake}, width::Int, height::Int)

    fields = zeros(length(snakes))
    visited = Set()
    for snek in snakes
        push!(visited, get_body_no_head(snek))
    end

    stack = Tuple{Position, Int}[]
    for (p, snek) in enumerate(snakes)
        push!(stack, (get_head(snek), p))
    end

    while !isempty(stack)
        pos, p = popfirst!(stack)
        if !(pos in visited)
            push!(visited, pos)
            fields[p] += 1
            for d in (up, down, left, right)
                next_pos = move_from(pos, d)
                if !(next_pos in visited) && check_bounds(next_pos, height, width)
                    push!(stack, (next_pos, p))
                end
            end
        end
    end

    return fields
end

function a_star_search(start_field:: Position, search_field::Position, board::Boardstate)

    queue = PriorityQueue{Position, Int}()
    came_from = Dict{Position, Tuple{Position, Direction}}()
    cost_so_far = Dict{Position, Int}()
    cost_so_far[start_field] = 0
    queue[start_field] = 0

    while !isempty(queue)
        node = dequeue!(queue)
        if node == search_field
            break
        end
        for (dir, nb) in get_neighbors(node)
            s = Inf
            for snake in board.snakes
                if nb in snake.body
                    idx = findfirst(b -> b==nb, snake.body)
                    s = length(snake.body) - 1 - idx
                end
            end
            # next pos is / will be free, in bounds and no circles
            if (!is_occupied_by_snake(nb, board) || h(start_field, nb) > s) && check_bounds(nb, board.height, board.width) && get(came_from, nb, nothing) === nothing
                c = cost_so_far[node] + h(nb, search_field)
                cost_so_far[nb] = cost_so_far[node] + 1
                came_from[nb] = (node, dir)
                queue[nb] = c
            end
        end
    end

    prev = search_field
    direction = nothing
    path = Vector{Tuple{Union{Direction, Nothing}, Position}}()
    push!(path, (direction, prev))
    while true
        prev, direction = get(came_from, prev, (nothing, nothing))
        if prev === nothing
            return nothing, nothing
        end
        push!(path, (direction, prev))
        if prev == start_field
            break
        end
    end
    return cost_so_far[search_field], path
end

function h(x::Position, y::Position)
    dx = abs(x.x - y.x)
    dy = abs(x.y - y.y)
    return dx + dy
end

max_stars_repo_name: rhotertj/solide-schlange

max_issues_repo_name: rhotertj/solide-schlange

max_forks_repo_name: rhotertj/solide-schlange


text:
# SPDX-License-Identifier: MIT

"""
    abstract type AbstractModel end

Abstract type for all the models.
"""
abstract type AbstractModel end

"""
    kldiv(model, qθ, Tθ)

KL divergence between the variational posterior qθ and the prior
of `model`. `Tθ` is the expected sufficient statistics w.r.t. `qθ`.
"""
kldiv(model::AbstractModel, qθ, Tθ) =
    kldiv(qθ, model.pθ, Tθ)

"""
    loglikelihood(model::AbstractModel, X, uTθ)

Return the log-likelihood of X given `model`. `uTθ` are the unpacked
sufficient statistics of the parameters.
"""
loglikelihood

"""
    elbo(model::AbstractModel, X, qθ, Tθ = getstats(model, qθ))

Evidence Lower-BOund of the `X` given `model` and the variational posterior
`qθ`. `Tθ` are the vectorized sufficient statistics of the parameters.
"""
function elbo(model::AbstractModel, X, qθ, Tθ = getstats(model, qθ); klscale=1,
              llhscale=1)
    uTθ = unpack(model, qθ, Tθ)
    llh = sum(loglikelihood(model, X, uTθ))
    kl = kldiv(model, qθ, Tθ)
    llhscale*llh - klscale*kl
end

"""
    getstats(model, qθ)

Return the expected value of the sufficient statistics of the parameters
`Tθ` w.r.t. `qθ`.
"""
getstats(model::AbstractModel, qθ::ExponentialFamilyDistribution) = μ(qθ)

"""
    unpack(model, qθ, Tθ)

Split the parameters statistics `Tθ` into its main components.
"""
unpack(model::AbstractModel, qθ::ExponentialFamilyDistribution,
       Tθ::AbstractVector) = unpack(qθ, Tθ)

"""
    newposterior(model, qθ, ∇μ; lrate=1, clip=true)

Create a new posterior whose parameter are estimated by a natural gradient
step with learning rate `lrate`.
"""
function newposterior(::AbstractModel, qθ::ExponentialFamilyDistribution, ∇μ;
                      lrate=1, clip=true)
    c∇μ = clip ? ∇μ ./ norm(∇μ) : ∇μ
    ηq = η(qθ)
    ∇̃ξ = ForwardDiff.derivative(t -> ξ(qθ, ηq + t * c∇μ), 0)
    ξᵗ⁺¹ = ξ(qθ, ηq) + lrate*∇̃ξ
    typeof(qθ)(unpack(qθ, ξᵗ⁺¹)...)
end

max_stars_repo_name: BUTSpeechFIT/PPCA

max_issues_repo_name: BUTSpeechFIT/PPCA

max_forks_repo_name: BUTSpeechFIT/PPCA


text:
using StochasticDiffEq, Plots, DiffEqDevTools, DiffEqProblemLibrary
using ParameterizedFunctions, DiffEqMonteCarlo, Base.Test, OrdinaryDiffEq
using Plots; plotly()
using BenchmarkTools

################################################################################

## Lorenz

f = @ode_def_nohes LorenzSDE begin
  dx = σ*(y-x)
  dy = x*(ρ-z) - y
  dz = x*y - β*z
end σ ρ β

σ = (du,u,p,t) -> begin
  for i in 1:3
    du[i] = 10.0 #Additive
  end
end

p = (10.0,28.0,2.66)
prob_sde_lorenz = SDEProblem(f,σ,ones(3),(0.0,10.0),p)

monte_prob = MonteCarloProblem(prob_sde_lorenz)
N = 100; fails = zeros(5,4); times = zeros(5,4)
offset = 0
for i in 1:5
  sol1 =solve(monte_prob,SRA1();num_monte=N,abstol=10,reltol=1/2^(2),verbose=false)
  sol2 =solve(monte_prob,SRA3();num_monte=N,abstol=10,reltol=1/2^(offset+i),verbose=false)
  sol3 =solve(monte_prob,SOSRA();num_monte=N,abstol=10,reltol=1/2^(offset+i),verbose=false)
  sol4 =solve(monte_prob,SOSRA2();num_monte=N,abstol=10,reltol=1/2^(offset+i),verbose=false)
  fails[i,1] = N - sum([sol1[i].retcode == :Success for i in 1:N])
  fails[i,2] = N - sum([sol2[i].retcode == :Success for i in 1:N])
  fails[i,3] = N - sum([sol3[i].retcode == :Success for i in 1:N])
  fails[i,4] = N - sum([sol4[i].retcode == :Success for i in 1:N])
  times[i,1] = sol1.elapsedTime
  times[i,2] = sol2.elapsedTime
  times[i,3] = sol3.elapsedTime
  times[i,4] = sol4.elapsedTime
end


text:
module ForwardDiffBenchmarks

using ForwardDiff
using BenchmarkTools

include(joinpath(dirname(dirname(@__FILE__)), "test", "utils.jl"))

name(f) = last(split(string(f), '.'))

const SUITE = BenchmarkGroup()

xs = map(n -> rand(MersenneTwister(1), n), (10, 100, 1000, 10000))
chunk_sizes = 1:10

val = addgroup!(SUITE, "value")
deriv = addgroup!(SUITE, "derivative")
grad = addgroup!(SUITE, "gradient")
jac = addgroup!(SUITE, "jacobian")
hess = addgroup!(SUITE, "hessian")

for f in (NUMBER_TO_NUMBER_FUNCS..., NUMBER_TO_ARRAY_FUNCS...)
    x = 1
    val[name(f)] = @benchmarkable $(f)($x)
    deriv[name(f)] = @benchmarkable ForwardDiff.derivative($f, $x)
end

for f in VECTOR_TO_NUMBER_FUNCS
    fval = addgroup!(val, name(f))
    fgrad = addgroup!(grad, name(f))
    fhess = addgroup!(hess, name(f))
    for x in xs
        xlen = length(x)
        fval[xlen] = @benchmarkable $(f)($x)
        for c in chunk_sizes
            fgrad[xlen, c] = @benchmarkable ForwardDiff.gradient($f, $x, $(Chunk{c}()))
            fhess[xlen, c] = @benchmarkable ForwardDiff.hessian($f, $x, $(Chunk{c}()))
        end
    end
end

for (f!, f) in VECTOR_TO_VECTOR_FUNCS
    fval = addgroup!(val, name(f))
    fjac = addgroup!(jac, name(f))
    for x in xs
        xlen = length(x)
        fval[xlen] = @benchmarkable $(f)($x)
        for c in chunk_sizes
            fjac[xlen, c] = @benchmarkable ForwardDiff.jacobian($f, $x, $(Chunk{c}()))
        end
    end
end

function runafew(seconds = 1)
    result = BenchmarkGroup()

    deriv = addgroup!(result, "derivative")
    tune!(SUITE["derivative"], seconds = seconds) # fast enough to require tuning
    for (f, b) in SUITE["derivative"]
        deriv[f] = run(b; seconds = seconds)
    end

    grad = addgroup!(result, "gradient")
    for (f, group) in SUITE["gradient"]
        grad[f] = run(group[1000, 10]; seconds = seconds)
    end

    hess = addgroup!(result, "hessian")
    for (f, group) in SUITE["hessian"]
        hess[f] = run(group[100, 10]; seconds = seconds)
    end

    jac = addgroup!(result, "jacobian")
    for (f, group) in SUITE["jacobian"]
        jac[f] = run(group[1000, 10]; seconds = seconds)
    end

    return result
end

end # module


text:
using FromFile
@from "subchain/subfolder/subchain1.jl" import a
@from "subchain/subchain2.jl" import B.b, B.SubB.b2
@from "subchain/subchain3.jl" import C: c
@from "subchain/subchain3.jl" import C2.SubC: c2
@from "subchain/subchain4.jl" using D  # d
@from "subchain/subchain5.jl" using E: e


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

using ModelBlocks
using Test
using OrdinaryDiffEq
using Plots

println("Block Test starting");

variables = Variables([
    Variable("X", 1, (0, 100), "kg", ""),
    Variable("Y", 1, (0, 100), "kg", ""),
    Variable("Z", 1, (0, 100), "kg", ""),
]);

parameters = Variables([
    Variable("Xin", 1.5, (0, 100), "kg/s", ""),
    Variable("Yin", 1.0, (0, 100), "kg/s", ""),
    Variable("Zout", 1.0, (0, 100), "kg/s", ""),
    Variable("Xout", .2, (0, 100), "kg/s", ""),
    Variable("Kxy", 0.5, (0, 100), "?", ""),
    Variable("R", .1, (0, 100), "?", ""),
]);

reactions = [
    SimpleReaction("X uptake", [], ["X"], "Xin"),
    SimpleReaction("Y uptake", [], ["Y"], "Yin"),
    SimpleReaction("Z clearance", ["Z"], [], "Zout"),
    SimpleReaction("X clearance", ["X"], [], "Xout"),
    SimpleReaction("Transformation", ["X", "Y"], ["Z"], "Kxy"),
    GeneralRateReaction("Random stuff", ["X"], ["Y"], (t, v, p) -> v.X * cos(t) * p.R),
    GeneralReaction("Really random stuff", (dvdt, t, v, p, e) -> begin
        dvdt.X += 0.1;
        dvdt.Y -= v.Z * 0.1;
        e["Extra"] = 3;
    end),
    GeneralRateReaction("Random stuff", ["X"], ["Y"], (t, v, p, e) -> e["Extra"]),
];

block = Block(variables, parameters, reactions);
setTimeRange!(block, 0:100);

@time solution = runBlock(block);

outputs = Variables([
    Variable("vector", 0.0, (0, 100), "?", ""),
    Variable("scalar", 0.0, (0, 100), "?", ""),
]);

setOutputDefinition!(block, outputs, (variables, parameters, timeRange, solution, outputs) -> begin
    outputs.vector = [solution.X[1], solution.X[5]];
    outputs.scalar = solution.Y[end];

    return outputs;
end);

@time fit = fitParameters([block], [
    ("Xin", 0., 10.0),
    ("Yin", 0., 10.0),
], Dict(
    "vector" => ([1, 2], [1, 1]),
    "scalar" => (3, 1)
); MaxTime = 3);

println("Block Test finished");

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
using Distributions
export SimMcmcSampler
export SimMcmcInsertDelete, SimMcmcInsertDeleteGibbs, SimMcmcInsertDeleteSubpath, SimMcmcInsertDeleteProportional
export SimMcmcInsertDeleteLengthCentered

abstract type SimMcmcSampler end

# Edit Allocation Sampler
# -----------------------

struct SimMcmcInsertDelete <: SimMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of trans-dimensional move
    len_dist::DiscreteUnivariateDistribution
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values to insert in interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td_add::Vector{Int} # Storage of where to insert/delete
    ind_td_del::Vector{Int} # Storage of where to insert/delete
    function SimMcmcInsertDelete(
        ;
        K=100,
        ν_ed=2, ν_td=2, β=0.4, len_dist=TrGeometric(0.8,1,K),
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td_add = zeros(Int, ν_td)
        ind_td_del = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_ed, ν_td, β,
            len_dist, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td_add, ind_td_del
            )
    end
end

function Base.show(io::IO, sampler::SimMcmcInsertDelete)
    title = "MCMC Sampler for SIM Models"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Description: interaction insertion/deletion with edit allocation updates.")
    println(io, "Parameters:")
    num_of_pars = 5
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

# Edit allocation sampler with edits proportional to path size

struct SimMcmcInsertDeleteProportional <: SimMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of trans-dimensional move
    len_dist::DiscreteUnivariateDistribution
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values to insert in interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SimMcmcInsertDeleteProportional(
        ;
        K=100,
        ν_ed=2, ν_td=2, β=0.4, len_dist=TrGeometric(0.8,1,K),
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_ed, ν_td, β,
            len_dist, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td
            )
    end
end

function Base.show(io::IO, sampler::SimMcmcInsertDeleteProportional)
    title = "MCMC Sampler for SIM Models"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Description: interaction insertion/deletion with edit allocation updates proportional to path size.")
    println(io, "Parameters:")
    num_of_pars = 5
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

struct SimMcmcInsertDeleteLengthCentered <: SimMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of trans-dimensional move
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values to insert in interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SimMcmcInsertDeleteLengthCentered(
        ;
        K=100,
        ν_ed=2, ν_td=2, β=0.5,
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_ed, ν_td, β, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td
            )
    end
end

function Base.show(io::IO, sampler::SimMcmcInsertDeleteLengthCentered)
    title = "MCMC Sampler for SIM Models"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Description: interaction insertion/deletion with edit allocation updates proportional to path size.")
    println(io, "Parameters:")
    num_of_pars = 4
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

# Edit Allocation Sampler (for LSP ground dist)
# ---------------------------------------------

struct SimMcmcInsertDeleteSubpath <: SimMcmcSampler
    ν_ed::Int  # Maximum number of edit operations
    ν_td::Int  # Maximum change in outer dimension
    β::Real  # Probability of trans-dimensional move
    len_dist::DiscreteUnivariateDistribution
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    init::McmcInitialiser
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions of interactions
    ind_add::Vector{Int} # Storage for indexing of additions of interactions
    vals::Vector{Int} # Storage for values to insert in interactions
    ind_update::Vector{Int} # Storage of which values have been updated
    ind_td::Vector{Int} # Storage of where to insert/delete
    function SimMcmcInsertDeleteSubpath(
        ;
        K=100,
        ν_ed=2, ν_td=2, β=0.4, len_dist=TrGeometric(0.8,1,K),
        desired_samples=1000, lag=1, burn_in=0, init=InitMode()
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν_ed)
        ind_add = zeros(Int, ν_ed)
        vals = zeros(Int, ν_ed)
        ind_update = zeros(Int, ν_ed)
        ind_td = zeros(Int, ν_td)
        par_info = Dict()
        par_info[:ν_ed] = "(maximum number of edit operations)"
        par_info[:ν_td] = "(maximum increase/decrease in dimension)"
        par_info[:len_dist] = "(distribution to sample length of path insertions)"
        par_info[:β] = "(probability of update move)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν_ed, ν_td, β,
            len_dist, K,
            desired_samples, burn_in, lag, init,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals,
            ind_update, ind_td
            )
    end
end

function Base.show(io::IO, sampler::SimMcmcInsertDeleteSubpath)
    title = "MCMC Sampler for SIM Models"
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Description: interaction insertion/deletion with edit allocation updates preserving subpaths.")
    println(io, "Parameters:")
    num_of_pars = 5
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+4)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end

# Gibbs Scan Sampler
# ------------------

struct SimMcmcInsertDeleteGibbs<: SimMcmcSampler
    ν::Int   # Maximum number of edit ops
    path_dist::PathDistribution  # Distribution used to introduce new interactions
    β::Real  # Extra probability of Gibbs move
    K::Int # Max number of interactions (used to determined how many pointers to store interactions)
    desired_samples::Int  # Final three set default values for MCMC samplers
    burn_in::Int
    lag::Int
    par_info::Dict
    curr_pointers::InteractionSequence{Int} # Storage for prev value in MCMC
    prop_pointers::InteractionSequence{Int} # Storage for curr value in MCMC
    ind_del::Vector{Int} # Storage for indexing of deletions in Gibbs scan
    ind_add::Vector{Int} # Storage for indexing of additions in Gibbs scan
    vals::Vector{Int} # Storage for valuse to insert in Gibbs scan
    function SimMcmcInsertDeleteGibbs(
        path_dist::PathDistribution;
        K=100,
        ν=4, β=0.0,
        desired_samples=1000, lag=1, burn_in=0
        )
        curr_pointers = [Int[] for i in 1:K]
        prop_pointers = [Int[] for i in 1:K]
        ind_del = zeros(Int, ν)
        ind_add = zeros(Int, ν)
        vals = zeros(Int, ν)
        par_info = Dict()
        par_info[:ν] = "(maximum number of edit operations)"
        par_info[:path_dist] = "(path distribution for insertions)"
        par_info[:β] = "(extra probability of Gibbs scan)"
        par_info[:K] = "(maximum number of interactions, used to initialise storage)"

        new(
            ν, path_dist, β, K,
            desired_samples, burn_in, lag,
            par_info,
            curr_pointers, prop_pointers, ind_del, ind_add, vals
            )
    end
end

function Base.show(io::IO, sampler::SimMcmcInsertDeleteGibbs)
    title = "MCMC Sampler for SIM Models via Gibbs + Insert/Delete Moves."
    n = length(title)
    println(io, title)
    println(io, "-"^n)
    println(io, "Parameters:")
    num_of_pars = 4
    for par in fieldnames(typeof(sampler))[1:num_of_pars]
        println(io, par, " = $(getfield(sampler, par))  ", sampler.par_info[par])
    end
    println(io, "\nDefault output parameters:")
    for par in fieldnames(typeof(sampler))[(num_of_pars+1):(num_of_pars+3)]
        println(io, par, " = $(getfield(sampler, par))  ")
    end
end


text:
using CBindingGen
import P4est_jll

# Get include directory from jll package
const incdir = joinpath(dirname(dirname(P4est_jll.libp4est_path)), "include")

# Manually set header files to consider
const hdrs = ["p4est.h", "p8est.h", "p4est_extended.h", "p8est_extended.h"]

# Convert symbols in header
cvts = convert_headers(hdrs, args = ["-I", incdir]) do cursor
	header = CodeLocation(cursor).file
	name   = string(cursor)

	# only wrap the libp4est headers
	startswith(header, "$(incdir)/") || return false

  # Ignore macro hacks
  startswith(name, "sc_extern_c_hack_") && return false

	return true
end

# Write generated C bindings to file
const bindings_filename = joinpath(@__DIR__, "libp4est.jl")
open(bindings_filename, "w+") do io
	generate(io, P4est_jll.libp4est_path => cvts)
end


text:
using JSON
using StatsBase

include("./BayesianBlocks.jl")

function readJSON()
    open("./config_demos.json") do f
        r = JSON.parse(f)
        return r
    end
end

function tst_tecno()
    S = readJSON()
    x=readdlm(S["dat1"], '\n')
    x=x[:]
    println(typeof(x))
    edges_in = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0]
    println(x)
    norm_hist = fit(Histogram, x, edges_in)
    # edges_norm = norm_hist.edges
    # counts_norm = norm_hist.weights
    # edges_bayes = bayesian_blocks(edges_norm, x)
    # bayes_hist = fir(Histogram, x, edges_bayes)
    # counts_bayes = bayes_hist.weights
    # println(counts_norm)
    # println(counts_bayes)
end
tst_tecno()

max_stars_repo_name: av-maslov/Recurrency

max_issues_repo_name: av-maslov/Recurrency

max_forks_repo_name: av-maslov/Recurrency


text:
export DIM1, DIM2, DIM3, Bounds

SVectorF{N} = SVector{N,Float}

struct Dim{N}
    Dim{1}() = new()
    Dim{2}() = new()
    Dim{3}() = new()
end

const DIM1 = Dim{1}()
const DIM2 = Dim{2}()
const DIM3 = Dim{3}()

struct Bounds{N,T}
    lower::SVector{N,T}
    upper::SVector{N,T}

    function Bounds(
        lower::SVector{N,T},
        upper::SVector{N,T},
    ) where {N,T<:Number}
        @check all(upper > lower)
        new{N,T}(lower, upper)
    end

    function Bounds(lower::T, upper::T) where {T<:Number}
        @check upper > lower
        new{1,T}(SVector(lower), SVector(upper))
    end
end

lower(bounds::Bounds) = bounds.lower

upper(bounds::Bounds) = bounds.upper

"""
    extents(bounds::Bounds{N,T}) => SVector{N,T}

Compute the extents of the region within the given bounds.
"""
extents(bounds::Bounds) = bounds.upper - bounds.lower

"""
    extended(bounds::Bounds{N,T}, lengths::SVector{N,T}) => Bounds{N,T}

Compute new bounds extended symmetrically by the given lengths.
"""
function extended(bounds::Bounds{N,T}, lengths::SVector{N,T}) where {N,T}
    half_lengths = 0.5 * lengths
    Bounds(bounds.lower - half_lengths, bounds.upper + half_lengths)
end

max_stars_repo_name: lars-frogner/Whirl

max_issues_repo_name: lars-frogner/Whirl

max_forks_repo_name: lars-frogner/Whirl


text:
 function solve(pfem)

   #=
   casetable = createcasetable(elementdf,
      data[:wobrange], data[:inclinationrange],
      data[:computefunction], data[:pdir])

   pmap(runcase, casetable)
   =#
   fem, dis_df, fm_df = PtFEM.p44(pfem)
   (fem, dis_df, fm_df)
 end


text:
module Bathymetry

import DataFrames: DataFrame
import Geodesy: LLA, euclidean_distance
import LinearAlgebra: norm

const Soundings = Vector{NTuple{3,Float64}}

function read(dirname)
  data = NTuple{3,Float64}[]
  for filename ∈ filter(s -> endswith(s, ".txt"), readdir(dirname, join=true))
    open(filename, "r") do io
      for s ∈ eachline(io)
        m = match(r"^<SOUNDG .*>\( (.*) \)", s)
        if m !== nothing
          for x ∈ split(m[1], " ), ( ")
            try
              y = parse.(Float64, split(x, ","))
              push!(data, (y[1], y[2], y[3]))
            catch
              @warn "Bad data: $s"
            end
          end
        end
      end
    end
  end
  data
end

function _bathy(s::Soundings, p)
  r = [norm(p .- x[1:2]) for x ∈ s]
  ndx = sortperm(r)
  r = r[ndx]
  r[1] == 0.0 && return s[ndx[1]][3]
  i = findfirst(r1 -> r1 > 20 * r[1], r)
  j = i === nothing ? length(r) : i - 1
  w = 1 ./ (r[1:j].^2)
  d = [x[3] for x ∈ s[ndx[1:j]]]
  sum(d .* w) / sum(w)
end

bathy(s::Soundings, p::LLA) = _bathy(s, (p.lon, p.lat))

function bathy(s::Soundings, p1::LLA, p2::LLA; spacing=10.0)
  d = euclidean_distance(p1, p2)
  n = round(Int, d / spacing)
  x = LinRange(p1.lon, p2.lon, n)
  y = LinRange(p1.lat, p2.lat, n)
  DataFrame(range=LinRange(0, d, n), location=LLA.(y, x), depth=[_bathy(s, (x[i], y[i])) for i ∈ 1:n])
end

end # module


text:
"
Auxiliary structure that keeps some generic information for STGP experiments.

# Fields
- `_functionSet::Array{FunctionNode}`: array that contains all the non-terminal
    nodes specified in the node information file.
- `_terminalSet::Array{TerminalNode}`: array that contains all the terminal nodes
    specified in the node information file.
- `_maxTreeDepth::UInt16`: maximum tree depth permitted.
- `_variables::Array{String}`: array that contains the names of the variables of
    the problem.

See also: [`GPExperimentInfo`](@ref)
"
struct STGPInfo <: GPExperimentInfo
    _functionSet::Array{FunctionNode}
    _terminalSet::Array{TerminalNode}
    _maxTreeDepth::UInt16
    _variables::Array{String}
end # struct

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
abstract type AbstractRobot end
#функции от обычного робота

mutable struct CheatRobot <:AbstractRobot
    robot::Robot
    coords::NTuple{2,Integer}

    CheatRobot(robot::Robot, coords::NTuple{2,Integer}) = new( robot, (0,0) )

    function move!(r::CheatRobot, side::HorizonSide)

        if side == Nord
            r.coords[2] += 1
        end

        function try_move_var2!(robot, side)
            if isborder(r,side)
                steps = 0
                nextside = next_side_pr(side)
                while isborder(robot, side)
                    if !isborder(robot, nextside)
                        move!(robot, nextside)
                        steps += 1
                    else
                        for j in 1:steps
                            move!(robot, anti_side(nextside))
                        end
                        return false
                    end
                end

                move!(robot, side)
                while isborder(robot, anti_side(nextside))
                    move!(robot, side)
                end

                move_for!(robot, anti_side(nextside), steps)
                return true
            else
                move!(r,side)
            end
        end
        # -----------------------------------------------------------------------------------------------------
        function try_move!(robot, side)
            if isborder(r,side)
                flag = true
                steps = 0
                nextside = next_side(side)
                while isborder(robot, side)
                    if !isborder(robot, nextside)
                        move!(robot, nextside)
                        steps += 1
                    else
                        for j in 1:steps
                            move!(robot, anti_side(nextside))
                        end
                        flag = false
                        try_move_var2!(r,side)
                    end
                end
                if (flag)
                    move!(robot, side)
                    while isborder(robot, anti_side(nextside))
                        move!(robot, side)
                    end

                    move_for!(robot, anti_side(nextside), steps)
                    return true
                end
            else
                move!(r,side)
            end
        end
        # -----------------------------------------------------------------------------------------------------
    end

    function anti_side(side::HorizonSide)::HorizonSide
        return HorizonSide((Int(side) + 2) % 4)
    end
    # по часовой
    function next_side(side::HorizonSide)
        return anti_side(HorizonSide((Int(side) + 1) % 4))
    end
    # против часовой
    function next_side_pr(side::HorizonSide)
        return HorizonSide((Int(side) + 1) % 4)
    end
    # -----------------------------------------------------------------------------------------------------
    function move_for!(robot,side::HorizonSide,steps)
        for step in 1:steps
            move!(robot,side)
        end
    end
    # -----------------------------------------------------------------------------------------------------
    function try_move_var2!(robot, side)
        if isborder(r,side)
            steps = 0
            nextside = next_side_pr(side)
            while isborder(robot, side)
                if !isborder(robot, nextside)
                    move!(robot, nextside)
                    steps += 1
                else
                    for j in 1:steps
                        move!(robot, anti_side(nextside))
                    end
                    return false
                end
            end

            move!(robot, side)
            while isborder(robot, anti_side(nextside))
                move!(robot, side)
            end

            move_for!(robot, anti_side(nextside), steps)
            return true
        else
            move!(r,side)
        end
    end
    # -----------------------------------------------------------------------------------------------------
    function try_move!(robot, side)
        if isborder(r,side)
            flag = true
            steps = 0
            nextside = next_side(side)
            while isborder(robot, side)
                if !isborder(robot, nextside)
                    move!(robot, nextside)
                    steps += 1
                else
                    for j in 1:steps
                        move!(robot, anti_side(nextside))
                    end
                    flag = false
                    try_move_var2!(r,side)
                end
            end
            if (flag)
                move!(robot, side)
                while isborder(robot, anti_side(nextside))
                    move!(robot, side)
                end

                move_for!(robot, anti_side(nextside), steps)
                return true
            end
        else
            move!(r,side)
        end
    end
    # -----------------------------------------------------------------------------------------------------
    # moving_in_angle! передвинет робота в Юго-Западный угол поля и вернёт в переменную back_bath ОБРАТНЫЙ путь робота
    function moving_in_angle!(robot::Robot)
        back_path = []
        while (!(isborder(robot, West)) || !(isborder(robot, Sud)))
            for side in (Sud, West)
                steps = 0
                while (!(isborder(robot,side)))
                    move!(robot,side)
                    steps += 1
                end
                push!(back_path, steps)
            end
        end
        return reverse!(back_path)
    end
    # moving_back_to_start! вернёт робота в начальное положение через обратный путь робота
    function moving_back_to_start!(robot::Robot, back_path::Vector)
        for (index,value) in enumerate(back_path)
            for step in 1:value
                move!(robot, HorizonSide((((index + 1)% 2) + 3) % 4))
            end
        end
    end
    # -----------------------------------------------------------------------------------------------------

end

max_stars_repo_name: Igetsuga/AlgorithmicsHomework

max_issues_repo_name: Igetsuga/AlgorithmicsHomework

max_forks_repo_name: Igetsuga/AlgorithmicsHomework


text:
function _find_span(n, m, M, pos)
    pos > M && return nothing
    pos < m && return nothing
    stride = (M - m) / n
    centers = (m + 0.5stride):stride:(M-0.5stride)
    span_pos = last(findmin(abs.(pos .- centers)))
    return (stride, centers[span_pos], span_pos)
end

"""
    geotiff(::Type{LT}, tiff_file; left::T=-180.0, right::T=180.0, bottom::T=-90.0, top::T=90.0) where {LT <: SimpleSDMLayer, T <: Number}

The geotiff function reads a geotiff file, and returns it as a matrix of the
correct type. The optional arguments `left`, `right`, `bottom`, and `left` are
defining the bounding box to read from the file. This is particularly useful if
you want to get a small subset from large files.

The first argument is the type of the `SimpleSDMLayer` to be returned.
"""
function geotiff(
    ::Type{LT},
    ::Type{ST},
    tiff_file;
    left = nothing,
    right = nothing,
    bottom = nothing,
    top = nothing,
) where {LT<:SimpleSDMLayer,ST<:SimpleSDMSource}

    left = isnothing(left) ? minimum(longitudes(ST)) : left
    right = isnothing(right) ? maximum(longitudes(ST)) : right
    bottom = isnothing(bottom) ? minimum(latitudes(ST)) : bottom
    top = isnothing(top) ? maximum(latitudes(ST)) : top

    # We do a bunch of checking that the required bounding box is not out of bounds
    # for the range of latitudes and longitudes.
    @assert right > left
    @assert top > bottom

    # This next block is reading the geotiff file, but also making sure that we
    # clip the file correctly to avoid reading more than we need.
    ArchGDAL.read(tiff_file) do dataset

        # The data we need is pretty much always going to be stored in the first
        # band, so this is what we will get for now. Note that this is not
        # reading the data yet, just retrieving the metadata.
        band = ArchGDAL.getband(dataset, 1)

        # This next bit of information is crucial, as it will allow us to assign
        # a matrix of the correct size, but also to get the right latitudes and
        # longitudes.
        width = ArchGDAL.width(dataset)
        height = ArchGDAL.height(dataset)

        global lon_stride, lat_stride
        global left_pos, right_pos
        global bottom_pos, top_pos

        lon_stride, left_pos, min_width = _find_span(width, minimum(longitudes(ST)), maximum(longitudes(ST)), left)
        _, right_pos, max_width = _find_span(width, minimum(longitudes(ST)), maximum(longitudes(ST)), right)

        lat_stride, top_pos, max_height = _find_span(height, minimum(latitudes(ST)), maximum(latitudes(ST)), top)
        _, bottom_pos, min_height = _find_span(height, minimum(latitudes(ST)), maximum(latitudes(ST)), bottom)

        max_height, min_height = height .- (min_height, max_height) .+ 1

        # We are now ready to initialize a matrix of the correct type.
        pixel_type = ArchGDAL.pixeltype(band)
        buffer = Matrix{pixel_type}(undef, length(min_width:max_width), length(min_height:max_height))
        ArchGDAL.read!(dataset, buffer, 1, min_height:max_height, min_width:max_width)
    end

    buffer = convert(Matrix{Union{Nothing,eltype(buffer)}}, rotl90(buffer))
    replace!(buffer, minimum(buffer) => nothing)

    return LT(buffer, left_pos-0.5lon_stride, right_pos+0.5lon_stride, bottom_pos-0.5lat_stride, top_pos+0.5lat_stride)

end


text:
using VSL
import Base.Test: @test

brng_types = [
    VSL_BRNG_MCG31, VSL_BRNG_R250, VSL_BRNG_MRG32K3A, VSL_BRNG_MCG59, VSL_BRNG_WH, VSL_BRNG_SOBOL, VSL_BRNG_NIEDERR,
    VSL_BRNG_MT19937, VSL_BRNG_MT2203, VSL_BRNG_SFMT19937, VSL_BRNG_NONDETERM, VSL_BRNG_ARS5, VSL_BRNG_PHILOX4X32X10
]

macro test_rand(expr)
    values = @eval $expr
    if isa(values, Array)
        for value in values
            @printf "%.6f\n" value
        end
    else
        @printf "%.6f\n" values
    end
end


text:
# const gjt = GeoJSONTables
# const shp = Shapefile

struct read_gjt{T<:Array}
    properties::T
    geomtery::T
end
""" Accepts a geojson table and returns it's properties and geometry """
function read_gjt(geojson_table)
    properties = []
    geometry = []
    for i in geojson_table
        push!(properties, get_props(i))
        push!(geometry, get_geom(i))
    end
    return read_gjt(properties, geometry)
end

""" Getter method for fetching properties of a geojson table"""
function get_props(geojson_table_row)
    keys_prop = Tuple(propertynames(geojson_table_row));
    vals_prop= ()
    for ky in keys_prop
        vals_prop = (vals_prop..., getproperty(geojson_table_row, ky));
    end
    return NamedTuple{keys_prop}(vals_prop);
end

""" Getter method for fetching geometry of a geojson table"""
function get_geom(geojson_table_row)
    gj_geom = getfield(geojson_table_row, :json).geometry
    keys_geom = Tuple(propertynames(gj_geom))
    vals_geom = ()
    for k in keys_geom
        vals_geom = (vals_geom..., getproperty(gj_geom, k))
    end
    return NamedTuple{keys_geom}(vals_geom)
end

struct read_shp{T<:Array}
    properties::T
    geometry::T
end

""" Accepts a shapefile table and returns it's properties and geometry """
function read_shp(shp_table)
    properties = []
    geometry = []
    for i in shp_table
        push!(properties, get_props(i))
        push!(geometry, get_geom(i))
    end
    return read_shp(properties, geometry)
end

""" Getter method for fetching properties of a geojson table"""
function get_props(shp_table_row)
    shp_prop = getfield(shp_table_row, :record)
    return shp_prop
end

""" Getter method for fetching geometry of a shapefile row"""
function get_geom(shp_table_row)
    vals_geom = getfield(shp_table_row, :geometry) #this
    return vals_geom
end

# function gdf_from_geojson(geojson_source)
#     parsed_geojson = GeoJSON.read(geojson_source)
#     feat = parsed_geojson.features

#     tbl = GeoTable(feat)

#     return DataFrame(tbl)
# end

# # # function read(parsed_geojson)

# function gdf_from_shp(shp)
#     tbl = GeoTable_sp(shp)
#     return DataFrame(tbl)
# end

# # function GJT_read()


text:
using Gadfly

set_default_plot_size(6inch, 6inch)

z = rand(1:8,100)
plot(x=rand(100), y=rand(100), shape=z, Geom.point,
    Scale.shape_discrete(levels=sort(unique(z))))


text:
export DegreeDominatingSet

struct DegreeDominatingSet end

"""
    update_dominated!(degree_queue, v, dominated, in_dom_set)

Check if a vertex is already dominated.
If not, make it dominated and update `degree_queue` by decreasing
the priority of the vertices adjacent to `v` by 1.
"""
function update_dominated!(
    g::AbstractGraph{T},
    degree_queue::PriorityQueue,
    v::Integer,
    dominated::BitArray{1},
    in_dom_set::BitArray{1}
    ) where T <: Integer

    @inbounds if !dominated[v]
        dominated[v] = true
        if !in_dom_set[v]
            degree_queue[v] -= 1
        end
        @inbounds @simd for u in neighbors(g, v)
            if !in_dom_set[u]
                degree_queue[u] -= ifelse(in_dom_set[u], 0, 1)
            end
        end
    end
end

"""
    dominating_set(g, DegreeDominatingSet())

Obtain a [dominating set](https://en.wikipedia.org/wiki/Dominating_set) using a greedy heuristic.

### Implementation Notes
A vertex is said to be dominated if it is in the dominating set or adjacent to a vertex
in the dominating set.
Initialise the dominating set to an empty set and iteratively choose the vertex that would
dominate the most undominated vertices.

### Performance
Runtime: ``\\mathcal{O}((|V|+|E|)*log(|V|))``
Memory: ``\\mathcal{O}(|V|)``
Approximation Factor: `ln(maximum(degree(g)))+2`
"""
function dominating_set(
    g::AbstractGraph{T},
    alg::DegreeDominatingSet
    ) where T <: Integer

    nvg = nv(g)
    in_dom_set = falses(nvg)
    dominated = falses(nvg)
    degree_queue = PriorityQueue(Base.Order.Reverse, enumerate(degree(g) .+ 1))
    length_ds = 0

    while !isempty(degree_queue) && peek(degree_queue)[2] > 0
        v = dequeue!(degree_queue)
        in_dom_set[v] = true
        length_ds += 1

        update_dominated!(g, degree_queue, v, dominated, in_dom_set)
        for u in neighbors(g, v)
            update_dominated!(g, degree_queue, u, dominated, in_dom_set)
        end
    end

    return LightGraphs.findall!(in_dom_set, Vector{T}(undef, length_ds))
end


text:
using Documenter
using TimeSeriesEcon
using ModelBaseEcon
using StateSpaceEcon

# Scans the src/Tutorials path and populates the tutorials array used below
# include("scan_tutorials.jl")

makedocs(
    sitename="StateSpaceEcon",
    format=Documenter.HTML(prettyurls=get(ENV, "CI", nothing) == "true"),
    modules=[StateSpaceEcon, ModelBaseEcon, TimeSeriesEcon],
    doctest=false,
    pages=[
        "index.md",
        "Tutorials" => [
            "Tutorials/index.md",
            "Tutorials/README.md",
            "Tutorials/US_SW07/main.md",
            "FRB/US" => "Tutorials/FRB-US/main.md",
            "Tutorials/TimeSeriesEcon/main.md",
        ],
        "Reference" => [
            "TimeSeriesEcon" => "Reference/TimeSeriesEcon.md",
            "ModelBaseEcon" => "Reference/ModelBaseEcon.md",
            "StateSpaceEcon" => "Reference/StateSpaceEcon.md",
        ],
        "Design Papers" => [
            "DesignPapers/index.md",
            "DesignPapers/final_conditions.md",
            "DesignPapers/log_variables.md",
        ],
        "Index" => "indexpage.md"
    ]
)

if get(ENV, "CI", nothing) == "true"
    deploydocs(repo="github.com/bankofcanada/DocsEcon.jl.git")
end


text:
import POMDPModelTools: obs_weight
import POMDPs: observation

struct P <: POMDP{Nothing, Nothing, Nothing} end

POMDPs.observation(::P, ::Nothing) = Deterministic(nothing)
@test obs_weight(P(), nothing, nothing, nothing, nothing) == 1.0

#=
@test !@implemented obs_weight(::P, ::Nothing, ::Nothing, ::Nothing, ::Nothing)
@test !@implemented obs_weight(::P, ::Nothing, ::Nothing, ::Nothing)
@test !@implemented obs_weight(::P, ::Nothing, ::Nothing)

obs_weight(::P, ::Nothing, ::Nothing, ::Nothing) = 1.0

@test @implemented obs_weight(::P, ::Nothing, ::Nothing, ::Nothing)
@test @implemented obs_weight(::P, ::Nothing, ::Nothing, ::Nothing, ::Nothing)
@test !@implemented obs_weight(::P, ::Nothing, ::Nothing)
=#


text:
using Sockets

include("commit_comment.jl")
event_request = create_event()
event_json = JSON.parse(IOBuffer(HTTP.payload(event_request)))
event = GitHub.event_from_payload!("commit_comment", event_json)

@testset "WebhookEvent" begin
    @test event.repository.name == "BenchmarkTrackers.jl"
    @test event.sender.login == "jrevels"
end # testset

@testset "EventListener" begin
    @test !(GitHub.has_valid_secret(event_request, "wrong"))
    @test GitHub.has_valid_secret(event_request, "secret")
    @test !(GitHub.is_valid_event(event_request, ["wrong"]))
    @test GitHub.is_valid_event(event_request, ["commit_comment"])
    @test !(GitHub.from_valid_repo(event, ["JuliaWeb/GitHub.jl"]))
    @test GitHub.from_valid_repo(event, ["JuliaCI/BenchmarkTrackers.jl"])
    @test GitHub.handle_event_request(event_request, x -> true,
                                    secret = "secret",
                                    events = ["commit_comment"],
                                    repos = ["JuliaCI/BenchmarkTrackers.jl"])

    @test begin
        listener = EventListener(x -> true;
                                secret = "secret",
                                repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                                events = ["commit_comment"],
                                forwards = ["http://bob.com", HTTP.URI("http://jim.org")])
        r = listener.handle_request(HTTP.Request())
        r.status == 400
    end
end

@testset "CommentListener" begin
    result = GitHub.handle_comment((e, m) -> m, event, GitHub.AnonymousAuth(), r"`RunBenchmarks\(.*?\)`", false, false)
    @test result.match == "`RunBenchmarks(\"binary\", \"unary\")`"
    @test begin
        listener = CommentListener((x, y) -> true, r"trigger";
                                secret = "secret",
                                repos = [Repo("JuliaCI/BenchmarkTrackers.jl"), "JuliaWeb/GitHub.jl"],
                                forwards = ["http://bob.com", HTTP.URI("http://jim.org")],
                                check_collab = false)
        r = listener.listener.handle_request(HTTP.Request())
        r.status == 400
    end
end

@testset "HTTPClientServer" begin
    auth = GitHub.JWTAuth(4123, "not_a_real_key.pem")

    function test_handler(event::WebhookEvent, phrase::RegexMatch)
       return HTTP.Messages.Response((phrase.match == "RunBenchmarks") ? 200 : 500)
    end
    listener = GitHub.CommentListener(test_handler, r"RunBenchmarks"; check_collab=false, auth=auth, secret=nothing)
    host = IPv4("127.0.0.1")
    port, sock = Sockets.listenany(host, 8001)

    srvrtask = @async GitHub.run(listener, sock, host, Int(port))

    server_started = false
    while !server_started
        try
            close(connect(host, port))
            server_started = true
        catch
            yield()
        end
    end

    resp = HTTP.request("POST", "http://$host:$port", event_request.headers, event_request.body)
    @test resp.status == 200

    try
        close(sock)
        wait(srvrtask)
    catch
    end
end


text:
export solve_trajopt_jump!

mutable struct SCPParam_TrajOpt <: SCPParamSpecial
  mu0         # initial penalty coefficient
  s0          # initial trust region size
  c           # step acceptance parameter
  τ_plus      # trust region expansion factor
  τ_minus     # trust region shrinkage factor
  k           # penalty scaling factor
  ftol        # convergence threshold for merit
  xtol        # convergence threhsold for x
  ctol        # constraint satisfaction threshold
  max_penalty_iteration
  max_convex_iteration
  max_trust_iteration

  ρ_vec::Vector   # trust region ratios
  mu_vec::Vector  # penalty coefficient
  s_vec::Vector   # trust region size
  xtol_vec::Vector
  ftol_vec::Vector
  ctol_vec::Vector
end

function SCPParam_TrajOpt(mu0, s0, c, τ_plus, τ_minus, k, ftol, xtol, ctol,max_penalty_iteration,max_convex_iteration,max_trust_iteration)
	SCPParam_TrajOpt(mu0, s0, c, τ_plus, τ_minus, k, ftol, xtol, ctol,
    max_penalty_iteration,max_convex_iteration,max_trust_iteration,[0.],[mu0], [s0], [0.], [0.], [0.])
end

#########
# JuMP
#########
function solve_trajopt_jump!(SCPS::SCPSolution, SCPP::SCPProblem, solver="Ipopt", max_iter=125, force=false; kwarg...)
	# Solves a sequential convex programming problem
	# Inputs:
	#		SCPS - SCP solution structure, including an initial trajectory
	#		SCPP - SCP problem definition structure
	#   max_iter - Set to enforce a maximum number of iterations (0 = unlimited)
	#   force - Set true to force further refinement iterations despite convergence

	N = SCPP.N
	param, model = SCPP.param, SCPP.PD.model

	param.alg = SCPParam_TrajOpt(SCPP.PD.model)
  mu0,s0 = param.alg.mu0, param.alg.s0
  ρ_vec,mu_vec,s_vec = param.alg.ρ_vec, param.alg.mu_vec, param.alg.s_vec
  τ_plus,τ_minus = param.alg.τ_plus, param.alg.τ_minus
  c = param.alg.c
  ftol,xtol,ctol = param.alg.ftol, param.alg.xtol, param.alg.ctol
  ftol_vec, xtol_vec, ctol_vec = param.alg.ftol_vec, param.alg.xtol_vec, param.alg.ctol_vec

  max_penalty_iteration,max_convex_iteration,max_trust_iteration =
    param.alg.max_penalty_iteration, param.alg.max_convex_iteration, param.alg.max_trust_iteration

  iter_cap = SCPS.iterations + max_iter
	SCPV = SCPVariables{JuMP.VariableRef, Array{JuMP.VariableRef}}()
	SCPS.SCPC = SCPConstraints(SCPP)
	SCPC = SCPS.SCPC

  initialize_model_params!(SCPP, SCPS.traj)
	push!(SCPS.J_true, cost_true(SCPS.traj, SCPS.traj, SCPP))
	param.obstacle_toggle_distance = model.clearance + 1. # TODO(ambyld): Generalize this, workspace-dependent

  constraints_satisfied = false
  xtol_satisfied = false

  old_penalty_traj, old_convex_traj = SCPS.traj, SCPS.traj

  for penalty_iteration in 1:max_penalty_iteration
    constraints_satisfied && break
    copy!(old_penalty_traj,SCPS.traj)

    for convex_iteration in 1:max_convex_iteration
      copy!(old_convex_traj,SCPS.traj)

      constraints_satisfied && break
      if xtol_satisfied
        xtol_satisfied = false
        break
      end
      for trust_iteration in 1:max_trust_iteration
        time_start = time_ns()

        if solver == "Mosek"
          SCPS.solver_model = Model(with_optimizer(Mosek.Optimizer; kwarg...))
        elseif solver == "Gurobi"
          SCPS.solver_model = Model(with_optimizer(Gurobi.Optimizer; kwarg...))
        elseif solver == "Ipopt"
          SCPS.solver_model = Model(with_optimizer(Ipopt.Optimizer; kwarg...))
        else 	# Default solver
          SCPS.solver_model = Model(with_optimizer(Ipopt.Optimizer; kwarg...))
        end

        # Set up, solve problem
        update_model_params!(SCPP, SCPS.traj)
        add_variables_jump!(SCPS, SCPV, SCPP)
        add_constraints_trajopt_jump!(SCPS, SCPV, SCPC, SCPP)   # TODO(acauligi)
        add_objective_trajopt_jump!(SCPS, SCPV, SCPC, SCPP)     # TODO(acauligi)

        set_start_value.(SCPV.X, SCPS.traj.X)
        set_start_value.(SCPV.U, SCPS.traj.U)
        set_start_value(SCPV.Tf, SCPP.tf_guess)

        JuMP.optimize!(SCPS.solver_model)

		    push!(SCPS.solver_status, JuMP.termination_status(SCPS.solver_model))
        if SCPS.solver_status[end] ∉ (MOI.OPTIMAL, MOI.LOCALLY_SOLVED)
          warn("TrajOpt iteration failed to find an optimal solution")
          push!(SCPS.iter_elapsed_times, (time_ns() - time_start)/10^9)
        end

        # Recover solution
		    new_traj = Trajectory(JuMP.value.(SCPV.X), JuMP.value.(SCPV.U), JuMP.value.(SCPV.Tf))
        push!(xtol_vec, evaluate_xtol(new_traj, old_convex_traj, SCPP))
        push!(SCPS.convergence_measure, xtol_vec[end])
		    push!(SCPS.J_full, JuMP.objective_value(SCPS.solver_model))

		    SCPS.dual = get_dual_jump(SCPC, SCPP)

        # grow or shrink trust region
        push!(ρ_vec, trust_region_ratio_trajopt(new_traj, old_convex_traj, SCPP))
        if ρ_vec[end] > c
          push!(s_vec, τ_plus*s_vec[end])
        else
          push!(s_vec, τ_minus*s_vec[end])
        end

        copy!(SCPS.traj, new_traj)
        iter_elapsed_time = (time_ns() - time_start)/10^9
        push!(SCPS.J_true, cost_true(new_traj, SCPS.traj, SCPP))
        push!(SCPS.iter_elapsed_times, iter_elapsed_time)
        SCPS.total_time += iter_elapsed_time
        SCPS.iterations += 1
        if s_vec[end] < xtol
          xtol_satisfied = true
          break
        end
      end

      push!(ftol_vec, evaluate_ftol(SCPS.traj, old_convex_traj, SCPP))
      push!(xtol_vec, evaluate_xtol(SCPS.traj, old_convex_traj, SCPP))
      if ftol_vec[end] < ftol || xtol[end] < xtol
        constraints_satisfied = true
        break
      end
    end

    if evaluate_ctol(SCPS.traj, old_penalty_traj, SCPP, SCPC) < ctol
  	 constraints_satisfied = true
     SCPS.converged = true
     break
  	else
      # increase penalty coefficient
      push!(mu_vec, mu_vec[end]*SCPP.param.alg.k)
    end
  end
end

function add_constraints_trajopt_jump!(SCPS::SCPSolution, SCPV::SCPVariables, SCPC::SCPConstraints, SCPP::SCPProblem)
  s = SCPP.param.alg.s_vec[end]
	solver_model, traj_prev = SCPS.solver_model, SCPS.traj

	update_model_params!(SCPP, traj_prev)

  for cc_list in values(SCPC.state_trust_region_ineq)
    for cc in cc_list
  		if cc.dimtype == :scalar
  			cc.con_reference = @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, k, i...) - s <= 0)
  		elseif cc.dimtype == :array
  			cc.con_reference = @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, k, i...) - s .<= 0)
  		end
    end
	end

  for cc_list in values(merge(SCPC.convex_state_boundary_condition_eq, SCPC.nonconvex_state_boundary_condition_convexified_eq))
    for cc in cc_list
      if cc.dimtype == :scalar
        cc.con_reference = @constraint(solver_model, [k=cc.params.k_timestep,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, cc.params, i...) == 0)
      elseif cc.dimtype == :array
        cc.con_reference = @constraint(solver_model, [k=cc.params.k_timestep,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, cc.params, i...) .== 0)
      end
    end
  end

  for cc_list in values(merge(SCPC.convex_state_boundary_condition_ineq, SCPC.nonconvex_state_boundary_condition_convexified_ineq))
    for cc in cc_list
  		if cc.dimtype == :scalar
  			cc.con_reference = @constraint(solver_model, [k=cc.params.k_timestep,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, cc.params, i...) <= 0)
  		elseif cc.dimtype == :array
  			cc.con_reference = @constraint(solver_model, [k=cc.params.k_timestep,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, cc.params, i...) .<= 0)
  		end
    end
	end

  for cc_list in values(SCPC.convex_state_eq)
    for cc in cc_list
  		if cc.dimtype == :scalar
  			cc.con_reference = @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, k, i...) == 0)
  		elseif cc.dimtype == :array
  			cc.con_reference = @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)], cc.func(SCPV, traj_prev, SCPP, k, i...) .== 0)
  		end
    end
	end

  # TODO: Store constraint reference to this
	if !SCPP.param.fixed_final_time
		@constraint(solver_model, SCPV.Tf >= 0.1)
	end
end

function add_objective_trajopt_jump!(SCPS::SCPSolution, SCPV::SCPVariables, SCPC::SCPConstraints, SCPP::SCPProblem)
	solver_model, traj_prev = SCPS.solver_model, SCPS.traj
	robot, model = SCPP.PD.robot, SCPP.PD.model
	x_dim, u_dim, N = model.x_dim, model.u_dim, SCPP.N

  mu = SCPP.param.alg.mu_vec[end]

	cost_expr = cost_true_convexified(SCPV, traj_prev, SCPP)

  # Inequality constraints
  for cc_list in values(merge(SCPC.convex_state_ineq,SCPC.nonconvex_state_convexified_ineq,SCPC.convex_control_ineq))
    for cc in cc_list
  		if cc.dimtype == :scalar
  			cc.var_reference = @variable(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)])
  			cc.con_reference = @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)],
  				mu*cc.func(SCPV, traj_prev, SCPP, k, i...) <= cc.var_reference[k,i])
  			@constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)],
  				-cc.var_reference[k,i] <= 0.)
  			cost_expr += sum(cc.var_reference[k,i] for k in cc.ind_time, i in Iterators.product(cc.ind_other...))
  			set_start_value.(cc.var_reference, 0.)
  		elseif cc.dimtype == :array
  			throw(":array type inequality constraint penalty not implemented")      # TODO(acauligi)
  		end
    end
	end

  # Equality constraints
  for cc_list in values(SCPC.nonconvex_state_convexified_eq)
    for cc in cc_list
  		if cc.dimtype == :scalar
  			cc.var_reference = @variable(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)])
  			cc.con_reference = @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)],
  				mu*cc.func(SCPV, traj_prev, SCPP, k, i...) <= cc.var_reference[k,i])
        @constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)],
  				-mu*cc.func(SCPV, traj_prev, SCPP, k, i...) <= cc.var_reference[k,i])         # TODO(acauligi): Need to track this second constraint satisfaction somehow
  			@constraint(solver_model, [k=cc.ind_time,i=Iterators.product(cc.ind_other...)],
  				-cc.var_reference[k,i] <= 0.)
  			cost_expr += sum(cc.var_reference[k,i] for k in cc.ind_time, i in Iterators.product(cc.ind_other...))
  			set_start_value.(cc.var_reference, 0.)
  		elseif cc.dimtype == :array
  			throw(":array type equality constraint penalty not implemented")      # TODO(acauligi)
  		end
    end
	end

  # Dynamics constraints
  for cc_list in values(SCPC.dynamics)
    for cc in cc_list
      if cc.dimtype == :scalar
        throw(":scalar type dynamics constraint penalty not implemented")      # TODO(acauligi)
      elseif cc.dimtype == :array
        # TODO(ambyld): Redesign this
        for k in cc.ind_time
          var_reference1 = @variable(solver_model, [j=1:x_dim])
          var_reference2 = @variable(solver_model, [j=1:x_dim])
          @constraint(solver_model, mu*cc.func(SCPV, traj_prev, SCPP, k) - var_reference1 .<= 0)
          @constraint(solver_model, var_reference2 - mu*cc.func(SCPV, traj_prev, SCPP, k) .<= 0)
          cost_expr += sum(var_reference1[j] + var_reference2[j] for j in 1:x_dim)
          set_start_value.(var_reference1, 0.)
          set_start_value.(var_reference2, 0.)
        end
      end
    end
  end

	@objective(solver_model, Min, cost_expr)
end

function evaluate_xtol(traj::Trajectory, traj_prev::Trajectory, SCPP::SCPProblem)
  return convergence_metric(traj,traj_prev,SCPP)
end

function evaluate_ftol(traj::Trajectory, traj_prev::Trajectory, SCPP::SCPProblem)
  abs(cost_true(traj, traj, SCPP)-cost_true(traj_prev, traj_prev, SCPP))/abs(cost_true(traj, traj, SCPP))
end

function evaluate_ctol(traj::Trajectory, traj_prev::Trajectory, SCPP::SCPProblem, SCPC::SCPConstraints)
	JNum, JDen = 0, 0
  JtestNum, JtestDen = 0, 0

  fprev = SCPC.convex_state_ineq[1].func
  for cc_list in values(merge(SCPC.convex_state_ineq, SCPC.nonconvex_state_ineq, SCPC.convex_state_boundary_condition_ineq, SCPC.nonconvex_state_boundary_condition_ineq, SCPC.convex_state_eq, SCPC.dynamics, SCPC.nonconvex_state_eq))
    for cc in cc_list
      fnext = cc.func
      for k in cc.ind_time, i in Iterators.product(cc.ind_other...)
        if fnext == fprev	# Evaluate and update maximum value while considering same class of constraints
          JtestNum = max(JtestNum, norm(fnext(traj, traj, SCPP, k, i...) -fnext(traj_prev, traj_prev, SCPP, k, i...)))
          JtestDen = max(JtestDen, norm(fnext(traj, traj, SCPP, k, i...)))
        else		# Entering new class of constraints
          JNum += JtestNum	# Add max from previous class of constraints to cost
          JDen += JtestDen
          fprev = fnext	# Update function handle representing current class of constraints
          JtestNum = norm(fnext(traj, traj, SCPP, k, i...) - fnext(traj_prev, traj_prev, SCPP, k, i...)) 	# Evaluate value of first member of new class of constraints
          JtestDen = norm(fnext(traj, traj, SCPP, k, i...))
        end
      end
    end
  end
	return JNum/JDen
end


text:
module TestingScientificCameras

using Test
using ScientificCameras

#@testset "Group" begin
#    @testset "SubGroup1" begin
#        @test booleanexpression
#        @test booleanexpression
#    end
#    @testset "SubGroup2" begin
#        @test booleanexpression
#        @test booleanexpression
#    end
#end

end # module


text:
# Model and discretization
model = Cartpole()
n,m = size(model)
N = 11
tf = 5.
dt = tf/(N-1)

# Initial and Final conditions
x0 = @SVector zeros(n)
xf = @SVector [0, pi, 0, 0]

# Objective
Q = 1.0e-2*Diagonal(@SVector ones(n))
Qf = 100.0*Diagonal(@SVector ones(n))
R = 1.0e-1*Diagonal(@SVector ones(m))
obj = LQRObjective(Q,R,Qf,xf,N)

# Constraints
u_bnd = 3.0
conSet = ConstraintList(n,m,N)
bnd = BoundConstraint(n,m, u_min=-u_bnd, u_max=u_bnd)
goal = GoalConstraint(xf)
add_constraint!(conSet, bnd, 1:N-1)
add_constraint!(conSet, goal, N:N)

# Initial conditions
X0 = [@SVector fill(0.0,n) for k = 1:N]
u0 = @SVector fill(0.01,m)
U0 = [u0 for k = 1:N-1]
Z = Traj(X0,U0, fill(dt, N))

# Inner constructor
prob = Problem{RK3}(model, obj, conSet, x0, xf, Z, N, 0.0, tf)
@test TO.integration(prob) == RK3
@test prob.x0 == x0
@test prob.xf == xf
@test prob.constraints === conSet
add_constraint!(conSet, goal, N-1)
@test length(TO.get_constraints(prob)) == length(conSet)
@test TO.num_constraints(TO.get_constraints(prob)) === conSet.p
@test prob.obj === obj
@test prob.tf ≈ tf
@test prob.N == N
@test states(prob) ≈ X0
@test controls(prob) ≈ U0

# Alternate constructor
prob = Problem(model, obj, xf, tf, x0=x0, constraints=conSet, X0=X0, U0=U0)
@test TO.integration(prob) == RK3
@test prob.x0 == x0
@test prob.xf == xf
@test prob.constraints === conSet
add_constraint!(conSet, goal, N-1)
@test length(TO.get_constraints(prob)) == length(conSet)
@test TO.num_constraints(TO.get_constraints(prob)) === conSet.p
@test prob.obj === obj
@test prob.tf ≈ tf
@test prob.N == N
@test states(prob) ≈ X0
@test controls(prob) ≈ U0

# Change integration
prob = Problem(model, obj, xf, tf, x0=x0, constraints=conSet, integration=RK2)
@test TO.integration(prob) == RK2

# Test defaults
prob = Problem(model, obj, xf, tf)
@test prob.x0 == zero(x0)
@test prob.N == N
@test all(all.(isnan, states(prob)))
@test controls(prob) == [zeros(m) for k = 1:N-1]
@test isempty(prob.constraints)
@test TO.integration(prob) == RK3
@test RobotDynamics.get_times(prob) ≈ range(0, tf; step=dt)

# Set initial trajectories
initial_states!(prob, 2 .* X0)
@test states(prob) ≈ 2 .* X0
initial_controls!(prob, 2 .* U0)
@test controls(prob) ≈ 2 .* U0

initial_trajectory!(prob, Z)
@test controls(prob) ≈ U0
@test states(prob) ≈ X0

# Use 2D matrices
X0_mat = hcat(X0...)
U0_mat = hcat(U0...)
initial_states!(prob, 3*X0)
initial_controls!(prob, 4*U0)
@test states(prob) ≈ 3 .* X0
@test controls(prob) ≈ 4 .* U0

prob = Problem(model, obj, xf, tf, X0=X0_mat, U0=U0_mat)
@test states(prob) ≈ X0
@test controls(prob) ≈ U0

# Use variable time steps
dts = rand(N-1)
dts = dts / sum(dts) * tf
prob = Problem(model, obj, xf, tf, dt=dts)
times = RobotDynamics.get_times(prob)
@test times[end] ≈ tf
@test times[1] ≈ 0
@test times[2] ≈ dts[1]
@test diff(times) ≈ dts

# Test initial and final conditions
prob = Problem(model, obj, Vector(xf), tf, x0=Vector(x0))
@test prob.x0 ≈ x0
@test prob.xf ≈ xf
prob = Problem(model, obj, MVector(xf), tf, x0=MVector(x0))
@test prob.x0 ≈ x0
@test prob.xf ≈ xf
@test prob.x0 isa MVector
@test prob.xf isa MVector

x0_ = rand(n)
TO.set_initial_state!(prob, x0_)
@test prob.x0 ≈ x0_
@test_throws DimensionMismatch TO.set_initial_state!(prob, rand(2n))

## Change initial and goal states
prob = Problem(model, copy(obj), xf, tf, x0=x0, constraints=deepcopy(conSet))
x0_new = @SVector rand(n)
TO.set_initial_state!(prob, x0_new)
@test TO.get_initial_state(prob) == x0_new
TO.set_initial_state!(prob, Vector(2*x0_new))
@test TO.get_initial_state(prob) ≈ 2*x0_new

# goal state, changing both objective and terminal constraint
@test conSet[2].xf ≈ xf
xf_new = @SVector rand(n)
TO.set_goal_state!(prob, xf_new)
@test prob.xf ≈ xf_new
@test prob.obj[1].q ≈ -Q*xf_new
@test prob.constraints[2].xf ≈ xf_new

# make sure it doesn't modify the orignal since they're copied
@test obj[1].q ≈ -Q*xf
@test conSet[2].xf ≈ xf

# don't modify the terminal constraint
prob = Problem(model, copy(obj), xf, tf, x0=x0, constraints=copy(conSet))
TO.set_goal_state!(prob, xf_new, constraint=false)
@test prob.xf ≈ xf_new
@test prob.obj[1].q ≈ -Q*xf_new
@test prob.obj[end].q ≈ -Qf*xf_new
@test prob.constraints[2].xf ≈ xf

# don't modify the objective, and leave off constraints
prob = Problem(model, copy(obj), xf, tf, x0=x0)
TO.set_goal_state!(prob, xf_new, objective=false)
@test prob.xf ≈ xf_new
@test prob.obj[1].q ≈ -Q*xf
@test prob.obj[end].q ≈ -Qf*xf

# check that it modifies the orignal objective and constraint list if not copied
prob = Problem(model, obj, xf, tf, x0=x0, constraints=copy(conSet))
TO.set_goal_state!(prob, xf_new)
@test obj[1].q ≈ -Q*xf_new
@test obj[end].q ≈ -Qf*xf_new
@test conSet[2].xf ≈ xf_new


text:
module Schemata

export Schema, TableSchema, ColumnSchema, # types
       compare,                           # core function
       readschema                         # read schema from config file

include("handle_validvalues.jl")
include("types.jl")
include("compare_table_to_schema.jl")
include("readwrite.jl")

using .handle_validvalues
using .types
using .compare_table_to_schema
using .readwrite

end