text:
#    Copyright 2015 Raytheon BBN Technologies
#
#     Licensed under the Apache License, Version 2.0 (the "License");
#     you may not use this file except in compliance with the License.
#     You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#     Unless required by applicable law or agreed to in writing, software
#     distributed under the License is distributed on an "AS IS" BASIS,
#     WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#     See the License for the specific language governing permissions and
#     limitations under the License.

using SCS

import Convex
import SparseArrays: spzeros

"""
 ϕ represents the isomorphism between complex and real matrices.
 e.g., see [Invariant semidefinite programs](http://arxiv.org/abs/1007.2905)
 by Bachoc et al.  If two arguments are given, they are taken to be the real and
 imaginary parts of a complex matrix.
"""
function ϕ(r,i)
    if size(r) != size(i)
        error("ϕ requires both arguments to be of the same size.")
    end
    [r i; -i r]
end

function ϕ(c)
    r = real(c)
    i = imag(c)
    ϕ(r,i)
end

"""
Extract the real part of a complex matrix represented as a real matrix
"""
function ϕr(m)
    m[1:div(end,2),1:div(end,2)]
end

"""
Extract the imag part of a complex matrix represented as a real matrix.
"""
function ϕi(m)
    m[1:div(end,2),div(end,2)+1:end]
end

"""
Extract the real and imag parts of a complex matrix represented as a real matrix.
"""
function ϕinv(m)
   ϕr(m) + 1im*ϕi(m)
end

# NOTE: This function is unused
"""
Compute the trace of the real representation of a complex matrix.
"""
function retrϕ(m)
    Convex.tr(ϕr(m))
end

function ket(i,d)
    v = spzeros(Float64,d,1)
    v[i+1] = 1.0
    return v
end

function bra(i,d)
    return ket(i,d)'
end

"""
E_(id_dim, ρ_rim)

Generates linear map E such that E(ρ) → 1 ⊗ ρ

Note that the dual of this map is F such that F*vec(σ ⊗ ρ) → trace(σ) vec(ρ),
in other words, E is the dual of the partial trace.
"""
function E_(id_dim, ρ_dim)
    M = spzeros(Float64,id_dim^2*ρ_dim^2,ρ_dim^2)
    for m in 0:ρ_dim-1
        for n in 0:ρ_dim-1
            for k in 0:id_dim-1
                M += Base.kron(ket(k,id_dim),ket(m,ρ_dim),ket(k,id_dim),ket(n,ρ_dim))*Base.kron(bra(m,ρ_dim),bra(n,ρ_dim))
            end
        end
    end
    return M
end

"""
involution(m)

Permutes the elements of a matrix so that it transforms the column
major representation of a linear map L into a matrix C that is
positive iff L is completely positive, Hermitian iff L maps
(vectorized) Hermitian matrices to Hermitian matrices. In other words,
it corresponds to the Choi-Jamiolkoski isomorphism.
"""
function involution(m)
    dsq = size(m,1) # we assume the matrix is square
    d   = Int(sqrt(dsq))
    return reshape(permutedims(reshape(m,(d,d,d,d)),[2,4,1,3]),(dsq,dsq))
end

"""
dnormcptp(L1,L2)

Computes the diamond norm of a linear superoperator `L` (i.e., a
linear transformation of operators). The superoperator must be
represented in column major form. In other words, it must be given
by a matrix that, when multiplying a vectorized (column major)
operator, it should result in the vectorized (column major)
representation of the result of the transformation.
"""
function dnormcptp2 end

let # wat09b
    global dnormcptp2
    local prev_dy, F

    prev_dy = -1

    function dnormcptp2(L1,L2)
        J = involution(L1-L2)

        dx = size(J,1) |> sqrt |> x -> round(Int,x)
        dy = dx

        if prev_dy != dy
            F = E_(dy,dx)
            prev_dy = dy
        end

        Jr = real(J)
        Ji = imag(J)

        Wr = Convex.Variable(dy*dx, dy*dx)
        Wi = Convex.Variable(dy*dx, dy*dx)

        Mr = Convex.Variable(dy*dx, dy*dx)
        Mi = Convex.Variable(dy*dx, dy*dx)

        ρr = Convex.Variable(dx, dx)
        ρi = Convex.Variable(dx, dx)

        prob = Convex.maximize( Convex.tr( Jr*Wr + Ji*Wi ) )

        prob.constraints += Convex.tr(ρr) == 1
        prob.constraints += Convex.tr(ρi) == 0

        Mr = reshape(F*vec(ρr), dy*dx, dy*dx)
        Mi = reshape(F*vec(ρi), dy*dx, dy*dx)

        prob.constraints += Convex.isposdef( ϕ(ρr,ρi) )
        prob.constraints += Convex.isposdef( ϕ(Wr,Wi) )
        prob.constraints += Convex.isposdef( ϕ(Mr,Mi) - ϕ(Wr,Wi) )

        Convex.solve!(prob, SCSSolver(verbose=0, eps=1e-6, max_iters=5_000,
            acceleration_lookback=1, alpha=1.9))

        if prob.status != :Optimal
            #println("DNORM_CPTP warning.")
            #println("Input: $(L)")
            #println("Input's Choi spectrum: $(eigvals(liou2choi(L)))")
            @warn("Diamond norm calculation did not converge.")
        end

        return 2*prob.optval
    end
end

"""
dnormcptp(L1,L2)

Computes the diamond norm distance between two linear completely
positive and trace preserving superoperators `L1` and `L2` . The
superoperators must be represented in column major form.
"""
function dnormcptp end

let # wat09b
    global dnormcptp
    local prev_dy, F

    prev_dy = -1

    function dnormcptp(L1,L2)
        J = involution(L1-L2)

        dx = size(J,1) |> sqrt |> x -> round(Int,x)
        dy = dx

        if prev_dy != dy
            F = E_(dy,dx)'
            prev_dy = dy
        end

        Jr = real(J)
        Ji = imag(J)

        Zr = Convex.Variable(dy*dx, dy*dx)
        Zi = Convex.Variable(dy*dx, dy*dx)

        pZr = reshape(F*vec(Zr), dx, dx)
        pZi = reshape(F*vec(Zi), dx, dx)

        prob = Convex.minimize( Convex.opnorm( ϕ(pZr, pZi) ) )

        prob.constraints += Convex.isposdef( ϕ(Zr,Zi) )
        prob.constraints += Convex.isposdef( ϕ(Zr,Zi) - ϕ(Jr,Ji) )

        Convex.solve!(prob, SCSSolver(verbose=0, eps=1e-6, max_iters=5_000,
            acceleration_lookback=1, alpha=1.9))

        if prob.status != :Optimal
            #println("DNORM_CPTP warning.")
            #println("Input: $(L)")
            #println("Input's Choi spectrum: $(eigvals(liou2choi(L)))")
            @warn("Diamond norm calculation did not converge.")
        end

        return 2*prob.optval
    end
end

"""
dnorm2(L1,L2)

Computes the diamond norm of a linear superoperator `L` (i.e., a
linear transformation of operators). The superoperator must be
represented in column major form. In other words, it must be given
by a matrix that, when multiplying a vectorized (column major)
operator, it should result in the vectorized (column major)
representation of the result of the transformation.
"""
function dnorm2 end

let # wat13b
    global dnorm2
    local prev_dy, F

    prev_dy = -1

    function dnorm2(L1,L2)
        J = involution(L1-L2)

        dx = size(J,1) |> sqrt |> x -> round(Int,x)
        dy = dx

        if prev_dy != dy
            F = E_(dy,dx)'
            prev_dy = dy
        end

        Jr = real(J)
        Ji = imag(J)

        Y0r = Convex.Variable(dy*dx, dy*dx)
        Y0i = Convex.Variable(dy*dx, dy*dx)

        Y1r = Convex.Variable(dy*dx, dy*dx)
        Y1i = Convex.Variable(dy*dx, dy*dx)

        prob = Convex.minimize( Convex.opnorm(ϕ( ρ0r, ρ0i )) + Convex.opnorm(ϕ( ρ1r, ρ1i )))

        ρ0r = reshape(F*vec(Y0r), dx, dx)
        ρ0i = reshape(F*vec(Y0i), dx, dx)

        ρ1r = reshape(F*vec(Y1r), dx, dx)
        ρ1i = reshape(F*vec(Y1i), dx, dx)

        prob.constraints += Convex.isposdef( ϕ(Y0r,Y0i) )
        prob.constraints += Convex.isposdef( ϕ(Y1r,Y1i) )
        prob.constraints += Convex.isposdef( ϕ( [ Y0r -Jr ; -Jr' Y1r ], [ Y0i -Xi ; Xi' Y1i ] ) )

        Convex.solve!(prob, SCSSolver(verbose=0, eps=1e-6, max_iters=5_000,
            acceleration_lookback=1, alpha=1.9))

        if prob.status != :Optimal
            #println("DNORM_CPTP warning.")
            #println("Input: $(L)")
            #println("Input's Choi spectrum: $(eigvals(liou2choi(L)))")
            @warn("Diamond norm calculation did not converge.")
        end

        return 2*prob.optval
    end
end

"""
dnorm(L; solver=Convex.get_default_solver())

Computes the diamond norm of a linear superoperator `L` (i.e., a
linear transformation of operators). The superoperator must be
represented in column major form. In other words, it must be given
by a matrix that, when multiplying a vectorized (column major)
operator, it should result in the vectorized (column major)
representation of the result of the transformation.
"""
function dnorm end

let # wat13b
    global dnorm
    local prev_dx, M

    prev_dx = -1

    function dnorm(L)
        J = involution(L)

        dx = size(J,1) |> sqrt |> x -> round(Int,x)
        dy = dx

        if prev_dx != dx
            M = E_(dy,dx)
            prev_dx = dx
        end

        Jr = real(J)
        Ji = imag(J)

        Xr  = Convex.Variable(dy*dx, dy*dx)
        Xi  = Convex.Variable(dy*dx, dy*dx)
        ρ0r = Convex.Variable(dx, dx)
        ρ0i = Convex.Variable(dx, dx)
        ρ1r = Convex.Variable(dx, dx)
        ρ1i = Convex.Variable(dx, dx)

        prob = Convex.maximize( Convex.tr( Jr*Xr + Ji*Xi ) )

        prob.constraints += Convex.tr(ρ0r) == 1
        prob.constraints += Convex.tr(ρ0i) == 0
        prob.constraints += Convex.tr(ρ1r) == 1
        prob.constraints += Convex.tr(ρ1i) == 0

        Mρ0r = reshape(M * vec(ρ0r), dy*dx, dy*dx)
        Mρ0i = reshape(M * vec(ρ0i), dy*dx, dy*dx)
        Mρ1r = reshape(M * vec(ρ1r), dy*dx, dy*dx)
        Mρ1i = reshape(M * vec(ρ1i), dy*dx, dy*dx)

        prob.constraints += Convex.isposdef( ϕ(ρ0r,ρ0i) )

        prob.constraints += Convex.isposdef( ϕ(ρ1r,ρ1i) )

        prob.constraints += Convex.isposdef( ϕ( [ Mρ0r Xr ; Xr' Mρ1r ], [ Mρ0i Xi ; -Xi' Mρ1i ] ) )

        Convex.solve!(prob, SCSSolver(verbose=0, eps=1e-6, max_iters=5_000,
            acceleration_lookback=1, alpha=1.9))

        if prob.status != :Optimal
            #println("DNORM warning.")
            #println("Input: $(L)")
            #println("Input's Choi spectrum: $(eigvals(liou2choi(L)))")
            @warn("Diamond norm calculation did not converge.")
        end

        return prob.optval
    end
end


text:
using Omega
using Distributions
import Omega: solve, isconditioned
using Random: MersenneTwister

struct ϵType end
const ϵ = ϵType()

rng = MersenneTwister(0);
N = 3;
xs = rand(rng, N)
m = 0.8
c = 1.2
linear_model(x, m, c) = m * x + c
obs_model(x) = linear_model(x, m, c) + randn(rng) * 0.1;
ys = obs_model.(xs);

M = @~ Normal(0, 1);
C = @~ Normal(0, 1);

## Simpler model
x = 0.3
ϵ_ = ϵ ~ Normal(0, 1)
y = M .* x .+ C .+ ϵ_
y_ = 21.0

y_ = m * x + c + randn()

e3 = (y .== y_)
ypost = M |ᶜ e3

# y_ = mx + c + ϵ

solve(::typeof(ϵ_), ω) = @show y_ - M(ω) * x - C(ω)
solve(::Member{StdNormal{Float64},ϵType}, ω) = @assert false

test() = complete(ypost)
# ##

# Y_class(i, ω) = linear_model(xs[i], M(ω), C(ω)) + (ϵ ∘ i ~ Normal(0, 0.1))(ω);

# EID = ComposedFunction{Type{ϵ}, Int64}
# Y⃗ = Mv(1:N, Y_class)
# evidence = Y⃗ ==ₚ ys
# joint_posterior = @joint(M, C) |ᶜ evidence

# solve(::StdNormal, i, ω) = (ϵ ∘ i ~ Normal(0, 0.1))(ω) / 0.1 # -
# solve(::typeof(Y_class), i, ω) = Y⃗(ω)[i] # -
# solve(::Normal, i::EID, ω) = Y_class(i, ω) - M(ω) - C(ω) # -
# solve(::typeof(Y⃗), ω) = ys # -
# solve(::typeof(evidence), ω) = true #

# ## Transform into conditioend model
# transform(::StdNormal, i, ω) = (ϵ ∘ i ~ Normal(0, 0.1))(ω) / 0.1 # -
# solve(::Normal, i::EID, ω) = Y_class(i, ω) - M(ω) - C(ω) #

# ## Problem is that we're doing all this work, really we just needed to evaluate StdNormal
# #

# isconditioned(::StdNormal) = true
# isconditioned(x) = false
# isconditioned(::Normal) = true

# # Propagation model
# # In this model, we define propagation functions which define which values of varaiables in the model conditioend on evidence
# function propagate(::typeof(evidence), x)
#   if x
#     Y⃗ => ys
#   else
#     nothing
#   end
# end

ext: jl

max_stars_repo_name: zenna/expect

max_issues_repo_name: zenna/expect

max_forks_repo_name: zenna/expect


text:
using PowerModels

function post_pfls{T}(pm::GenericPowerModel{T})

    PowerModels.variable_complex_voltage(pm, bounded=true) # overloaded: phase angle variables for DC
    variable_active_generation(pm)
    variable_reactive_generation(pm) # does nothing for DC
    variable_load_shed(pm) # TODO: add a PR to PowerModels; also overloaded
    variable_conductance_shed(pm) # variable that has 0 and vmax^2 as bounds; for each bus
    PowerModels.variable_active_line_flow(pm) # overloaded: p(i,j) = -p(j,i)
    PowerModels.variable_reactive_line_flow(pm) # overloaded: does nothing for DC

    objective_min_load_shed(pm) # TODO: add PR to PowerModels

    PowerModels.constraint_theta_ref(pm) # overloaded: does nothing for SOC -
    PowerModels.constraint_complex_voltage(pm) # overloaded: |W_ij|² ≦ Wᵢ⋅Wⱼ - does not return reference

    for (i, bus) in pm.set.buses
        constraint_active_kcl_ls(pm, bus)
        constraint_reactive_kcl_ls(pm, bus)
        #PowerModels.constraint_active_kcl_shunt(pm, bus)
        #PowerModels.constraint_reactive_kcl_shunt(pm, bus)
    end

    for (i, branch) in pm.set.branches
        constraint_active_ohms(pm, branch)
        constraint_reactive_ohms(pm, branch)
        #PowerModels.constraint_active_ohms_yt(pm, branch)
        #PowerModels.constraint_reactive_ohms_yt(pm, branch)

        PowerModels.constraint_phase_angle_difference(pm, branch)
        PowerModels.constraint_thermal_limit_from(pm, branch)
        PowerModels.constraint_thermal_limit_to(pm, branch)
    end

end

function post_ots{T}(pm::GenericPowerModel{T})
    variable_line_indicator(pm)
    PowerModels.variable_complex_voltage_on_off(pm) # adds phase angle variables for DC
    variable_active_generation(pm)
    variable_reactive_generation(pm) # does nothing for DC
    variable_load_shed(pm)
    PowerModels.variable_active_line_flow(pm)
    PowerModels.variable_reactive_line_flow(pm) # does nothing for DC

    objective_min_load_shed(pm)

    PowerModels.constraint_theta_ref(pm) # overloaded: does nothing for SOC
    PowerModels.constraint_complex_voltage(pm) # overloaded: |W_ij|² ≦ Wᵢ⋅Wⱼ - does not return reference

    for (i, bus) in pm.set.buses
        constraint_active_kcl_ls(pm, bus)
        constraint_reactive_kcl_ls(pm, bus)
    end

    for (i, branch) in pm.set.branches
        PowerModels.constraint_active_ohms_yt_on_off(pm, branch)
        PowerModels.constraint_reactive_ohms_yt_on_off(pm, branch)

        PowerModels.constraint_phase_angle_difference_on_off(pm, branch)
        PowerModels.constraint_thermal_limit_from_on_off(pm, branch)
        PowerModels.constraint_thermal_limit_to_on_off(pm, branch)
    end

end

function variable_line_indicator{T}(pm::GenericPowerModel{T})
    @variable(pm.model, 0 <= line_z[l in pm.set.branch_indexes] <= 1, start = PowerModels.getstart(pm.set.branches, 1, "line_z_start", 1.0))
    return line_z
end

function variable_active_generation{T}(pm::GenericPowerModel{T}; bounded = true)
    if bounded
        @variable(pm.model, 0 <= pg[i in pm.set.gen_indexes] <= pm.set.gens[i]["pmax"], start = PowerModels.getstart(pm.set.gens, i, "pg_start"))
    else
        @variable(pm.model, pg[i in pm.set.gen_indexes], start = PowerModels.getstart(pm.set.gens, i, "pg_start"))
    end
    return pg
end

function variable_reactive_generation{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T}; bounded = true)
    for i in pm.set.gen_indexes
        if (pm.set.gens[i]["qmin"] > 0 && pm.set.gens[i]["qmax"] > 0)
            pm.set.gens[i]["qmin"] = 0
        end
        if (pm.set.gens[i]["qmin"] < 0 && pm.set.gens[i]["qmax"] < 0)
            pm.set.gens[i]["qmax"] = 0
        end
    end
    if bounded
        @variable(pm.model, pm.set.gens[i]["qmin"] <= qg[i in pm.set.gen_indexes] <= pm.set.gens[i]["qmax"], start = PowerModels.getstart(pm.set.gens, i, "qg_start"))
    else
        @variable(pm.model, qg[i in pm.set.gen_indexes], start = PowerModel.getstart(pm.set.gens, i, "qg_start"))
    end
    return qg
end

function variable_reactive_generation{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T})
    # do nothing for dc model - no reactive variables
end

function variable_reactive_generation{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T}; bounded = true)
    for i in pm.set.gen_indexes
        if (pm.set.gens[i]["qmin"] > 0 && pm.set.gens[i]["qmax"] > 0)
            pm.set.gens[i]["qmin"] = 0
        end
        if (pm.set.gens[i]["qmin"] < 0 && pm.set.gens[i]["qmax"] < 0)
            pm.set.gens[i]["qmax"] = 0
        end
    end
    if bounded
        @variable(pm.model, pm.set.gens[i]["qmin"] <= qg[i in pm.set.gen_indexes] <= pm.set.gens[i]["qmax"], start = PowerModels.getstart(pm.set.gens, i, "qg_start"))
    else
        @variable(pm.model, qg[i in pm.set.gen_indexes], start = PowerModel.getstart(pm.set.gens, i, "qg_start"))
    end
    return qg
end

function variable_load_shed{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T})
    ld_min = [i => 0.0 for i in pm.set.bus_indexes]
    ld_max = [i => 0.0 for i in pm.set.bus_indexes]
    for i in pm.set.bus_indexes
        if (pm.set.buses[i]["pd"] > 0)
            ld_max[i] = 1.0
        end
    end
    @variable(pm.model, ld_min[i] <= ld[i in pm.set.bus_indexes] <= ld_max[i], start = 0.0)
    return ld
end

function variable_load_shed{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T})
    ld_min = [i => 0.0 for i in pm.set.bus_indexes]
    ld_max = [i => 0.0 for i in pm.set.bus_indexes]
    for i in pm.set.bus_indexes
        if (pm.set.buses[i]["pd"] > 0 || pm.set.buses[i]["qd"] != 0)
            ld_max[i] = 1.0
        end
    end
    @variable(pm.model, ld_min[i] <= ld[i in pm.set.bus_indexes] <= ld_max[i], start = 0.0)
    return ld
end

function variable_load_shed{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T})
    ld_min = [i => 0.0 for i in pm.set.bus_indexes]
    ld_max = [i => 0.0 for i in pm.set.bus_indexes]
    for i in pm.set.bus_indexes
        if (pm.set.buses[i]["pd"] > 0 || pm.set.buses[i]["qd"] != 0)
            ld_max[i] = 1.0
        end
    end
    @variable(pm.model, ld_min[i] <= ld[i in pm.set.bus_indexes] <= ld_max[i], start = 0.0)
    return ld
end

function variable_conductance_shed{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T})
    @variable(pm.model, 0.0 <= ld_c[i in pm.set.bus_indexes] <= 1.0, start = 0.0)
    return ld_c
end

function variable_conductance_shed{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T})
    @variable(pm.model, 0.0 <= ld_c[i in pm.set.bus_indexes] <= pm.set.buses[i]["vmax"]^2, start = 0.0)
    return ld_c
end

function variable_conductance_shed{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T})
    @variable(pm.model, 0.0 <= ld_c[i in pm.set.bus_indexes] <= 1.0, start = 0.0)
    return ld_c
end

function objective_min_load_shed{T}(pm::GenericPowerModel{T})
    c_pd = [bp => 1 for bp in pm.set.bus_indexes]
    for (i, bus) in pm.set.buses
        if (bus["pd"] < 0)
            c_pd[i] = 0
        end
    end
    ld = getvariable(pm.model, :ld)
    @objective(pm.model, Min, sum{pm.set.buses[i]["pd"] * c_pd[i] * ld[i], i=pm.set.bus_indexes} )
end

function constraint_active_kcl_ls{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T}, bus)
    i = bus["index"]
    bus_branches = pm.set.bus_branches[i]
    bus_gens = pm.set.bus_gens[i]

    pg = getvariable(pm.model, :pg)
    p_expr = pm.model.ext[:p_expr]
    ld = getvariable(pm.model, :ld)
    ld_c = getvariable(pm.model, :ld_c)

    c = @constraint(pm.model, sum{p_expr[a], a in bus_branches} == sum{pg[g], g in bus_gens} - bus["pd"]*(1-ld[i]) - bus["gs"]*ld_c[i])
    #c = @constraint(pm.model, sum{p_expr[a], a in bus_branches} == sum{pg[g], g in bus_gens} - bus["pd"]*(1-ld[i]))
    return Set([c])
end

function constraint_reactive_kcl_ls{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T}, bus)
    # do nothing
    return Set()
end

function constraint_active_kcl_ls{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T}, bus)
    i = bus["index"]
    bus_branches = pm.set.bus_branches[i]
    bus_gens = pm.set.bus_gens[i]

    w = getvariable(pm.model, :w)
    p = getvariable(pm.model, :p)
    pg = getvariable(pm.model, :pg)
    ld = getvariable(pm.model, :ld)
    ld_c = getvariable(pm.model, :ld_c)

    c = @constraint(pm.model, sum{p[a], a in bus_branches} == sum{pg[g], g in bus_gens} - bus["pd"]*(1-ld[i]) - bus["gs"]*ld_c[i])
    #c = @constraint(pm.model, sum{p[a], a in bus_branches} == sum{pg[g], g in bus_gens} - bus["pd"]*(1-ld[i]))
    return Set([c])
end

function constraint_reactive_kcl_ls{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T}, bus)
    i = bus["index"]
    bus_branches = pm.set.bus_branches[i]
    bus_gens = pm.set.bus_gens[i]

    w = getvariable(pm.model, :w)
    q = getvariable(pm.model, :q)
    qg = getvariable(pm.model, :qg)
    ld = getvariable(pm.model, :ld)
    ld_c = getvariable(pm.model, :ld_c)

    c = @constraint(pm.model, sum{q[a], a in bus_branches} == sum{qg[g], g in bus_gens} - bus["qd"]*(1-ld[i]) + bus["bs"]*ld_c[i])
    #c = @constraint(pm.model, sum{q[a], a in bus_branches} == sum{qg[g], g in bus_gens} - bus["qd"]*(1-ld[i]))
    return Set([c])
end

function constraint_active_kcl_ls{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T}, bus)
    i = bus["index"]
    bus_branches = pm.set.bus_branches[i]
    bus_gens = pm.set.bus_gens[i]

    v = getvariable(pm.model, :v)
    p = getvariable(pm.model, :p)
    pg = getvariable(pm.model, :pg)
    ld = getvariable(pm.model, :ld)
    ld_c = getvariable(pm.model, :ld_c)

    c = @NLconstraint(pm.model, sum{p[a], a in bus_branches} == sum{pg[g], g in bus_gens} - bus["pd"]*(1-ld[i]) - bus["gs"]*ld_c[i]*v[i]^2)
    #c = @constraint(pm.model, sum{p[a], a in bus_branches} == sum{pg[g], g in bus_gens} - bus["pd"]*(1-ld[i]))
    return Set([c])
end

function constraint_reactive_kcl_ls{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T}, bus)
    i = bus["index"]
    bus_branches = pm.set.bus_branches[i]
    bus_gens = pm.set.bus_gens[i]

    v = getvariable(pm.model, :v)
    q = getvariable(pm.model, :q)
    qg = getvariable(pm.model, :qg)
    ld = getvariable(pm.model, :ld)
    ld_c = getvariable(pm.model, :ld_c)

    c = @NLconstraint(pm.model, sum{q[a], a in bus_branches} == sum{qg[g], g in bus_gens} - bus["qd"]*(1-ld[i]) + bus["bs"]*ld_c[i]*v[i]^2)
    #c = @constraint(pm.model, sum{q[a], a in bus_branches} == sum{qg[g], g in bus_gens} - bus["qd"]*(1-ld[i]))
    return Set([c])
end

function constraint_active_ohms{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T}, branch)
    i = branch["index"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    p_fr = getvariable(pm.model, :p)[f_idx]
    t_fr = getvariable(pm.model, :t)[f_bus]
    t_to = getvariable(pm.model, :t)[t_bus]

    b = branch["b"]

    c = @constraint(pm.model, p_fr == -b*(t_fr - t_to))
    return Set([c])
end

function constraint_reactive_ohms{T <: PowerModels.AbstractDCPForm}(pm::GenericPowerModel{T}, branch)
    # do nothing
    return Set()
end

function constraint_active_ohms{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T}, branch)
    i = branch["index"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    p_fr = getvariable(pm.model, :p)[f_idx]
    p_to = getvariable(pm.model, :p)[t_idx]
    w_fr = getvariable(pm.model, :w)[f_bus]
    w_to = getvariable(pm.model, :w)[t_bus]
    wr = getvariable(pm.model, :wr)[(f_bus, t_bus)]
    wi = getvariable(pm.model, :wi)[(f_bus, t_bus)]

    g = branch["g"]
    b = branch["b"]

    c1 = @constraint(pm.model, p_fr == g*w_fr + (-g*wr) + (-b*wi) )
    c2 = @constraint(pm.model, p_to == g*w_to + (-g*wr) + (b*wi) )
    return Set([c1, c2])
end

function constraint_reactive_ohms{T <: PowerModels.AbstractWRForm}(pm::GenericPowerModel{T}, branch)
    i = branch["index"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    q_fr = getvariable(pm.model, :q)[f_idx]
    q_to = getvariable(pm.model, :q)[t_idx]
    w_fr = getvariable(pm.model, :w)[f_bus]
    w_to = getvariable(pm.model, :w)[t_bus]
    wr = getvariable(pm.model, :wr)[(f_bus, t_bus)]
    wi = getvariable(pm.model, :wi)[(f_bus, t_bus)]

    g = branch["g"]
    b = branch["b"]

    c1 = @constraint(pm.model, q_fr == -b*w_fr - (-b*wr) + (-g*wi) )
    c2 = @constraint(pm.model, q_to == -b*w_to - (-b*wr) + (g*wi) )
    return Set([c1, c2])
end

function constraint_active_ohms{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T}, branch)
    i = branch["index"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    p_fr = getvariable(pm.model, :p)[f_idx]
    p_to = getvariable(pm.model, :p)[t_idx]
    v_fr = getvariable(pm.model, :v)[f_bus]
    v_to = getvariable(pm.model, :v)[t_bus]
    t_fr = getvariable(pm.model, :t)[f_bus]
    t_to = getvariable(pm.model, :t)[t_bus]

    g = branch["g"]
    b = branch["b"]

    c1 = @NLconstraint(pm.model, p_fr == g*v_fr^2 + -g*v_fr*v_to*cos(t_fr-t_to) + -b*v_fr*v_to*sin(t_fr-t_to) )
    c2 = @NLconstraint(pm.model, p_to == g*v_to^2 + -g*v_fr*v_to*cos(t_to-t_fr) + -b*v_fr*v_to*sin(t_to-t_fr) )
    return Set([c1, c2])
end

function constraint_reactive_ohms{T <: PowerModels.AbstractACPForm}(pm::GenericPowerModel{T}, branch)
    i = branch["index"]
    f_bus = branch["f_bus"]
    t_bus = branch["t_bus"]
    f_idx = (i, f_bus, t_bus)
    t_idx = (i, t_bus, f_bus)

    q_fr = getvariable(pm.model, :q)[f_idx]
    q_to = getvariable(pm.model, :q)[t_idx]
    v_fr = getvariable(pm.model, :v)[f_bus]
    v_to = getvariable(pm.model, :v)[t_bus]
    t_fr = getvariable(pm.model, :t)[f_bus]
    t_to = getvariable(pm.model, :t)[t_bus]

    g = branch["g"]
    b = branch["b"]

    c1 = @NLconstraint(pm.model, q_fr == -b*v_fr^2 + b*v_fr*v_to*cos(t_fr-t_to) + -g*v_fr*v_to*sin(t_fr-t_to) )
    c2 = @NLconstraint(pm.model, q_to == -b*v_to^2 + b*v_fr*v_to*cos(t_to-t_fr) + -g*v_fr*v_to*sin(t_to-t_fr) )
    return Set([c1, c2])
end

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
module IfOnTop

include("./../../constraints/geometric/constants.jl")
using JuMP
using .Constants

export if_on_top

function if_on_top(m)

    #if (flbz_front_left_bot_z[i] == flbz_front_left_bot_z[j] + hob_height_of_box[j] &&
    #    flbz_front_left_bot_z[i] - hob_height_of_box[j] >= 0 &&
    #    flbx_front_left_bot_x[j] > flbx_front_left_bot_x[i] - box_i_x_axis_cover &&
    #    flbx_front_left_bot_x[j] < flbx_front_left_bot_x[i] + box_i_x_axis_cover &&
    #    flby_front_left_bot_y[j] > flby_front_left_bot_y[i] - box_i_y_axis_cover &&
    #    flby_front_left_bot_y[j] < flby_front_left_bot_y[i] + box_i_y_axis_cover)

    # IF 1
    @variable(m, if11_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @variable(m, if12_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @variable(m, if1_on_top[i=1:length(boxes), j=1:length(boxes); i!=j] >= 0, Bin)

    #Hz1 >= x-b+e
    #H(1-z1) >= b-x-e
    @constraint(m, c1_if11_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if11_on_top][i, j] >=
        m[:flbz_front_left_bot_z][i] - (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) + 0.001
    )
    @constraint(m, c2_if11_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if11_on_top][i, j]) >=
        (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - m[:flbz_front_left_bot_z][i] - 0.001
    )

    #Hz2 >= b-x+e
    #H(1-z2) >= x-b-e
    @constraint(m, c1_if12_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * m[:if12_on_top][i, j] >=
        (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - m[:flbz_front_left_bot_z][i] + 0.001
    )
    @constraint(m, c2_if12_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        H * (1 - m[:if12_on_top][i, j]) >=
        m[:flbz_front_left_bot_z][i] - (m[:flbz_front_left_bot_z][j] + m[:hob_height_of_box][j]) - 0.001
    )

    #If if11, if12 = 1, then box i and j have z-bot/top value within interval -e <= z <= +e, hence, expression = 1
    @constraint(m, c_if1_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if11_on_top][i, j] + m[:if12_on_top][i, j] >=
        2 * m[:if1_on_top][i, j]
    )

    #Xj < Xi + Xcoveri
    @variable(m, if2_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if2_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * (1 - m[:if2_on_top][i, j]) >=
        m[:flbx_front_left_bot_x][j] + 0.001 - (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i])
    )

    @constraint(m, c2_if2_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * m[:if2_on_top][i, j] >=
        (m[:flbx_front_left_bot_x][i] + m[:box_x_axis_cover][i]) - m[:flbx_front_left_bot_x][j] - 0.001
    )

    # Xj + Xcoverj > Xi
    @variable(m, if3_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if3_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * (1 - m[:if3_on_top][i, j]) >=
        m[:flbx_front_left_bot_x][i] + 0.001 - (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j])
    )

    @constraint(m, c2_if3_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        L * m[:if3_on_top][i, j] >=
        (m[:flbx_front_left_bot_x][j] + m[:box_x_axis_cover][j]) - m[:flbx_front_left_bot_x][i] - 0.001
    )

    # Yj < Yi + Ycoveri
    @variable(m, if4_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if4_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * (1 - m[:if4_on_top][i, j]) >=
        m[:flby_front_left_bot_y][j] + 0.001 - (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i])
    )

    @constraint(m, c2_if4_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * m[:if4_on_top][i, j] >=
        (m[:flby_front_left_bot_y][i] + m[:box_y_axis_cover][i]) - m[:flby_front_left_bot_y][j] - 0.001
    )

    # Yj + Ycoverj > Yi
    @variable(m, if5_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c1_if5_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * (1 - m[:if5_on_top][i, j]) >=
        m[:flby_front_left_bot_y][i] + 0.001 - (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j])
    )

    @constraint(m, c2_if5_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        W * m[:if5_on_top][i, j] >=
        (m[:flby_front_left_bot_y][j] + m[:box_y_axis_cover][j]) - m[:flby_front_left_bot_y][i] - 0.001
    )

    @variable(m, if_on_top[i=1:length(boxes), j=1:length(boxes); i!=j], Bin)
    @constraint(m, c_if_on_top[i=1:length(boxes), j=1:length(boxes); i!=j],
        m[:if1_on_top][i, j] +
        m[:if2_on_top][i, j] +
        m[:if3_on_top][i, j] +
        m[:if4_on_top][i, j] +
        m[:if5_on_top][i, j] >=
        5 * m[:if_on_top][i, j]
    )

    return m
end

end


text:
# Common facilities

#typealias
const AttributeDict = Dict{String, Any}

#################################################
#
#  vertex types
#
################################################

struct KeyVertex{K}
    index::Int
    key::K
    KeyVertex(idx::Int, key::K) where {K} = new{K}(idx, key)
    KeyVertex{K}(idx::Int, key::K) where {K} = new{K}(idx, key)
end

# KeyVertex{K}(idx::Int, key::K) = KeyVertex{K}(idx, key)
make_vertex(g::AbstractGraph{V}, key) where {V<:KeyVertex} = V(num_vertices(g) + 1, key)
vertex_index(v::KeyVertex) = v.index

mutable struct ExVertex
    index::Int
    label::String
    attributes::AttributeDict

    ExVertex(i::Int, label::AbstractString) = new(i, label, AttributeDict())
end

make_vertex(g::AbstractGraph{ExVertex}, label::AbstractString) = ExVertex(num_vertices(g) + 1, String(label))
vertex_index(v::ExVertex) = v.index
attributes(v::ExVertex, g::AbstractGraph) = v.attributes

const ProvidedVertexType = Union{KeyVertex, ExVertex}

vertex_index(v::V, g::AbstractGraph{V}) where {V<:ProvidedVertexType} = vertex_index(v)

function vertex_index(v::V, g::AbstractGraph{V}) where {V}
    @graph_requires g vertex_list
    if applicable(vertex_index, v)
        return vertex_index(v)
    end
    return vertex_index(v, vertices(g)) # slow linear search
end

vertex_index(v, vs::AA) where {AA <: AbstractArray} = something(findfirst(isequal(v), vs), 0) # findfirst(vs, v)

#################################################
#
#  edge types
#
################################################

struct Edge{V}
    index::Int
    source::V
    target::V
    Edge(i::Int, s::V, t::V) where {V} = new{V}(i, s, t)
end
const IEdge = Edge{Int}

# Edge{V}(i::Int, s::V, t::V) = Edge{V}(i, s, t)
make_edge(g::AbstractGraph{V,E}, s::V, t::V) where {V,E<:Edge} = Edge(num_edges(g) + 1, s, t)

revedge(e::Edge{V}) where {V} = Edge(e.index, e.target, e.source)

edge_index(e::Edge) = e.index
source(e::Edge) = e.source
target(e::Edge) = e.target
source(e::Edge{V}, g::AbstractGraph{V}) where {V} = e.source
target(e::Edge{V}, g::AbstractGraph{V}) where {V} = e.target

mutable struct ExEdge{V}
    index::Int
    source::V
    target::V
    attributes::AttributeDict
    ExEdge(i::Int, s::V, t::V) where {V} = new{V}(i, s, t, AttributeDict())
    ExEdge(i::Int, s::V, t::V, attrs::AttributeDict) where {V} = new{V}(i, s, t, attrs)
    ExEdge{V}(i::Int, s::V, t::V) where {V} = new{V}(i, s, t, AttributeDict())
    ExEdge{V}(i::Int, s::V, t::V, attrs::AttributeDict) where {V} = new{V}(i, s, t, attrs)
end

==(e1::ExEdge{V}, e2::ExEdge{V}) where {V} = (e1.index == e2.index &&
                                       e1.source == e2.source &&
                                       e1.target == e2.target)

# ExEdge{V}(i::Int, s::V, t::V) = ExEdge{V}(i, s, t, AttributeDict())
# ExEdge{V}(i::Int, s::V, t::V, attrs::AttributeDict) = ExEdge{V}(i, s, t, attrs)
make_edge(g::AbstractGraph{V}, s::V, t::V) where {V} = ExEdge(num_edges(g) + 1, s, t)

revedge(e::ExEdge{V}) where {V} = ExEdge{V}(e.index, e.target, e.source, e.attributes)

edge_index(e::ExEdge) = e.index
source(e::ExEdge) = e.source
target(e::ExEdge) = e.target
source(e::ExEdge{V}, g::AbstractGraph{V}) where {V} = e.source
target(e::ExEdge{V}, g::AbstractGraph{V}) where {V} = e.target
attributes(e::ExEdge, g::AbstractGraph) = e.attributes

#################################################
#
#  iteration
#
################################################

# general reindexed vector

struct ReindexedVec{T, Vec<:AbstractVector, I<:AbstractVector{Int}}
    src::Vec
    inds::I
    ReindexedVec(a::AT, inds::I) where {T, AT <: AbstractVector{T}, I <: AbstractVector{Int}} =
        new{T,typeof(a),typeof(inds)}(a, inds)
    ReindexedVec{T}(a::AT, inds::I) where {T, AT <: AbstractVector{T}, I <: AbstractVector{Int}} =
        new{T,typeof(a),typeof(inds)}(a, inds)
end

# ReindexedVec{T}(a::AbstractVector{T}, inds::AbstractVector{Int}) =
#     ReindexedVec{T,typeof(a),typeof(inds)}(a, inds)

length(a::ReindexedVec) = length(a.inds)
isempty(a::ReindexedVec) = isempty(a.inds)
getindex(a::ReindexedVec, i::Integer) = a.src[a.inds[i]]

# start(a::ReindexedVec) = start(a.inds)
# done(a::ReindexedVec, s) = done(a.inds, s)
# next(a::ReindexedVec, s) = ((i, s) = next(a.inds); (a.src[i], s))

# iterating over targets

struct TargetIterator{G<:AbstractGraph,EList}
    g::G
    lst::EList
    TargetIterator(g::G, lst::EList) where {G<:AbstractGraph,EList} = new{G,EList}(g, lst)
    TargetIterator{G,EList}(g::G, lst::EList) where {G<:AbstractGraph,EList} = new{G,EList}(g, lst)
end

# TargetIterator{G<:AbstractGraph,EList}(g::G, lst::EList) =
#     TargetIterator{G,EList}(g, lst)

Base.length(a::TargetIterator) = length(a.lst)
isempty(a::TargetIterator) = isempty(a.lst)
getindex(a::TargetIterator, i::Integer) = target(a.lst[i], a.g)

# start(a::TargetIterator) = start(a.lst)
# done(a::TargetIterator, s) = done(a.lst, s)
# next(a::TargetIterator, s::Int) = ((e, s) = next(a.lst, s); (target(e, a.g), s)) # likely deprecated

# Base.length(iter::TargetIterator) = length(iter.lst)
# Base.eltype(iter::TargetIterator) = ??
function Base.iterate(it::TargetIterator, (el, i)=(0, 0))
	return i >= length(it) ? nothing : (target(it.lst[i+1], it.g), (target(it.lst[i+1], it.g), i + 1))
end

# iterating over sources

struct SourceIterator{G<:AbstractGraph,EList}
    g::G
    lst::EList
    SourceIterator(g::G, lst::EList) where {G<:AbstractGraph,EList} =
        new{G,EList}(g, lst)
    SourceIterator{G,EList}(g::G, lst::EList) where {G<:AbstractGraph,EList} =
        new{G,EList}(g, lst)
end

# SourceIterator{G<:AbstractGraph,EList}(g::G, lst::EList) =
#     SourceIterator{G,EList}(g, lst)

length(a::SourceIterator) = length(a.lst)
isempty(a::SourceIterator) = isempty(a.lst)
getindex(a::SourceIterator, i::Integer) = source(a.lst[i], a.g)

# start(a::SourceIterator) = start(a.lst)
# done(a::SourceIterator, s) = done(a.lst, s)
# next(a::SourceIterator, s::Int) = ((e, s) = next(a.lst, s); (source(e, a.g), s))

# Base.length(iter::SourceIterator) = length(iter.lst)
# Base.eltype(iter::SourceIterator) = ??
function Base.iterate(it::SourceIterator, (el, i)=(0, 0))
  return i >= length(it) ? nothing : (source(it.lst[i+1], it.g), (source(it.lst[i+1], it.g), i + 1))
end

#################################################
#
#  Edge Length Visitors
#
################################################

abstract type AbstractEdgePropertyInspector{T} end

# edge_property_requirement{T, V}(visitor::AbstractEdgePropertyInspector{T}, g::AbstractGraph{V}) = nothing

mutable struct ConstantEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T}
  value::T
end

edge_property(visitor::ConstantEdgePropertyInspector{T}, e, g) where {T} = visitor.value

mutable struct VectorEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T}
  values::Vector{T}
end

edge_property(visitor::VectorEdgePropertyInspector{T}, e, g::AbstractGraph{V}) where {T,V} = visitor.values[edge_index(e, g)]

edge_property_requirement(visitor::AbstractEdgePropertyInspector{T}, g::AbstractGraph{V}) where {T, V} = @graph_requires g edge_map

mutable struct AttributeEdgePropertyInspector{T} <: AbstractEdgePropertyInspector{T}
  attribute::String
end

function edge_property(visitor::AttributeEdgePropertyInspector{T},edge::ExEdge, g) where {T}
    convert(T,edge.attributes[visitor.attribute])
end
#################################################
#
#  convenient functions
#
################################################

isnz(x::Bool) = x
isnz(x::Number) = x != zero(x)

intrange(n::Integer) = 1:convert(Int,n)

multivecs(::Type{T}, n::Int) where {T} = [T[] for _ =1:n]

function collect_edges(graph::AbstractGraph{V,E}) where {V,E}
    local edge_list
    if implements_edge_list(graph)
        collect(edges(graph))

    elseif implements_vertex_list(graph) && implements_incidence_list(graph)
        edge_list = Array{E}(undef, 0)
        sizehint!(edge_list, num_edges(graph))

        for v in vertices(graph)
            for e in out_edges(v, graph)
                push!(edge_list, e)
            end
        end
        edge_list
    else
        throw(ArgumentError("graph must implement either edge_list or incidence_list."))
    end
end

struct WeightedEdge{E,W}
    edge::E
    weight::W
end

isless(a::WeightedEdge{E,W}, b::WeightedEdge{E,W}) where {E,W} = a.weight < b.weight

function collect_weighted_edges(graph::AbstractGraph{V,E}, weights::AbstractEdgePropertyInspector{W}) where {V,E,W}

    edge_property_requirement(weights, graph)

    wedges = Array{WeightedEdge{E,W}}(undef, 0)
    sizehint!(wedges, num_edges(graph))

    if implements_edge_list(graph)
        for e in edges(graph)
            w = edge_property(weights, e, graph)
            push!(wedges, WeightedEdge(e, w))
        end

    elseif implements_vertex_list(graph) && implements_incidence_list(graph)
        for v in vertices(graph)
            for e in out_edges(v, graph)
                w = edge_property(weights, e, graph)
                push!(wedges, WeightedEdge(e, w))
            end
        end
    else
        throw(ArgumentError("graph must implement either edge_list or incidence_list."))
    end

    return wedges
end

function collect_weighted_edges(graph::AbstractGraph{V,E}, weights::AbstractVector{W}) where {V,E,W}
    visitor::AbstractEdgePropertyInspector{D} = VectorEdgePropertyInspector(edge_dists)
    collect_weighted_edges(graph, visitor)
end


text:
function luhn(id::AbstractString)
    id = split(replace(id, ' ' => ""), "")
    length(id) < 2 && return false
    all(all(isdigit, s) for s in id) || return false

    numbers = Int[]
    for val in id
        push!(numbers, parse(Int, val))
    end

    for i in length(numbers) - 1:-2:1
        numbers[i] = (numbers[i] * 2) % 9
    end

    return sum(numbers) % 10 == 0
end


text:
using seisimu
using PyCall
@pyimport matplotlib.pyplot as plt
include("extrapcump.jl")
pvel = 5000;
rho = 2.8;
DZ = 4000;
HX = 4000;
HY = 100;
pkf = 10;
T = 10;
ext = 20;
iflag = 2;

# === cuda ===#
model1 = initmodel(pvel,nothing,rho,DZ,HX,HY,pkf,T);

CPU_mem = 14
GPU_mem = 2.8
num_stream = 2
threadim = [8,4,8]

sn = 1
loc = [2000 2000 50]
ot = [0.05]
mt = [1 1 1 0 0 0]
wavelet = Ricker(pkf,model1.medium.dt)
waveform = wavelet
sou = initsource(sn, loc, ot, mt, waveform, model1)

@time extrap!(model1, threadim, CPU_mem, GPU_mem; num_stream = num_stream)

#=== openmp ===#
model2 = initmodel(pvel,nothing,rho,DZ,HX,HY,pkf,T; ext = 20);

sn = 1
loc = [2000 2000 50]
ot = [0.05]
mt = [1 1 1 0 0 0]
wavelet = Ricker(pkf,model2.medium.dt)
waveform = wavelet
sou = initsource(sn, loc, ot, mt, waveform, model2)

for it in 1:model2.medium.nT
    addsou!(model2.wf,sou,it)
    extrap!(model2)
    if mod(it,300) == 0
      plt.show(plt.imshow(model2.wf.vx[:,:,23]))
  end
end


text:
##### Numeric simplification

"""
```julia
simplify(x; expand=false,
            threaded=false,
            thread_subtree_cutoff=100,
            rewriter=nothing)
```

Simplify an expression (`x`) by applying `rewriter` until there are no changes.
`expand=true` applies [`expand`](/api/#expand) in the beginning of each fixpoint iteration.
"""
function simplify(x;
                  expand=false,
                  polynorm=nothing,
                  threaded=false,
                  thread_subtree_cutoff=100,
                  rewriter=nothing)
    if polynorm !== nothing
        Base.depwarn("simplify(..; polynorm=$polynorm) is deprecated, use simplify(..; expand=$polynorm) instead",
                        :simplify)
    end

    f = if rewriter === nothing
        if threaded
            threaded_simplifier(thread_subtree_cutoff)
        elseif expand
            serial_expand_simplifier
        else
            serial_simplifier
        end
    else
        Fixpoint(rewriter)
    end

    PassThrough(f)(x)
end

Base.@deprecate simplify(x, ctx; kwargs...)  simplify(x; rewriter=ctx, kwargs...)

"""
    substitute(expr, dict)

substitute any subexpression that matches a key in `dict` with
the corresponding value.
"""
function substitute(expr, dict; fold=true)
    haskey(dict, expr) && return dict[expr]

    if istree(expr)
        if fold
            canfold=true
            args = map(arguments(expr)) do x
                x′ = substitute(x, dict; fold=fold)
                canfold = canfold && !(x′ isa Symbolic)
                x′
            end
            canfold && return operation(expr)(args...)
            args
        else
            args = map(x->substitute(x, dict), arguments(expr))
        end
        similarterm(expr, operation(expr), args, metadata=metadata(expr))
    else
        expr
    end
end

"""
    occursin(needle::Symbolic, haystack::Symbolic)

Determine whether the second argument contains the first argument. Note that
this function doesn't handle associativity, commutativity, or distributivity.
"""
Base.occursin(needle::Symbolic, haystack::Symbolic) = _occursin(needle, haystack)
Base.occursin(needle, haystack::Symbolic) = _occursin(needle, haystack)
Base.occursin(needle::Symbolic, haystack) = _occursin(needle, haystack)
function _occursin(needle, haystack)
    isequal(needle, haystack) && return true

    if istree(haystack)
        args = arguments(haystack)
        for arg in args
            occursin(needle, arg) && return true
        end
    end
    return false
end


text:
using MAT

@testset "testfunctions" begin
@testset "single-obj" begin
    function test_optim(test_f, args...)
        f = test_f(args...)
        y, x = EfficientGlobalOptimization.optimum(f)
        abs(f(x)-y) < 1e-6
    end
    function test_inbounds(test_f, args...)
        f = test_f(args...)
        y, x = EfficientGlobalOptimization.optimum(f)
        lb = EfficientGlobalOptimization.lowerbounds(f)
        ub = EfficientGlobalOptimization.upperbounds(f)
        all(lb .<=x .<= ub)
    end
    @test test_inbounds(EfficientGlobalOptimization.Branin_Forrester)
    @test test_optim(EfficientGlobalOptimization.Branin_Forrester)

    for i = 2:6
        @test test_optim(EfficientGlobalOptimization.Ackley, i)
        @test test_optim(EfficientGlobalOptimization.Rosenbrock, i)
        @test test_optim(EfficientGlobalOptimization.SumSphere, i)
        @test test_inbounds(EfficientGlobalOptimization.Ackley, i)
        @test test_inbounds(EfficientGlobalOptimization.Rosenbrock, i)
        @test test_inbounds(EfficientGlobalOptimization.SumSphere, i)
    end
end

@testset "mult-obj" begin
    x = range(0.0, 1.0; length=11)
    zdt1 = EfficientGlobalOptimization.ZDT(1,3)
    zdt2 = EfficientGlobalOptimization.ZDT(2,3)
    zdt3 = EfficientGlobalOptimization.ZDT(3,3)
    nx = length(x)
    y1 = zeros(nx,nx,nx,2)
    y2 = zeros(nx,nx,nx,2)
    y3 = zeros(nx,nx,nx,2)
    for k=1:nx
        for j=1:nx
            for i=1:nx
                xx = [x[i],x[j],x[k]]
                y1[i,j,k,:] = zdt1(xx)
                y2[i,j,k,:] = zdt2(xx)
                y3[i,j,k,:] = zdt3(xx)
            end
        end
    end
    f = matopen("ZDT.mat")
    y1t = read(f, "zdt1")
    y2t = read(f, "zdt2")
    y3t = read(f, "zdt3")
    @test y1 ≈ y1t
    @test y2 ≈ y2t
    @test y3 ≈ y3t
end

@testset "constrained" begin
    f = EfficientGlobalOptimization.PressureVesselDesign()
    y, x = EfficientGlobalOptimization.optimum(f)
    lb = lowerbounds(f)
    ub = upperbounds(f)
    @test all(lb .<=x .<= ub)
    @test isapprox(y, f(x), atol=1e-6)
    g = constrains(f)
    @test length(g) == 4
    z = map(gg->gg(x), g)
    @test maximum(z) < sqrt(eps())
end
end


text:
content = readlines("input.txt")
template = content[1]
rules = Dict{String, String}()

foreach(split.(content[3:end], " -> ")) do rule
    rules[rule[1]] = rule[2]
end

generatePairs(str) = [str[i]*str[i+1] for i in 1:length(str)-1]

function counterCount(collection)
    dict = Dict()
    for item in collection
        dict[item] = get(dict, item, 0) + 1
    end
    return dict
end

function polymerizationCountsSorted(template, nSteps)
    pairs = generatePairs(template)

    pairCounts = counterCount(pairs)

    for _ in 1:nSteps
        subCounter = Dict()
        for (key, val) in pairCounts
            if !haskey(rules, key) return pairCounts end

            insertionCh = rules[key]
            pair1, pair2 = key[1] * insertionCh, insertionCh * key[2]

            subCounter[pair1] = get(subCounter, pair1, 0) + val
            subCounter[pair2] = get(subCounter, pair2, 0) + val
        end
        pairCounts = subCounter
    end

    counts = Dict{Char, Int}()

    for (keyPair, occurrences) in pairCounts
        counts[keyPair[1]] = get(counts, keyPair[1], 0) + occurrences
    end
    counts[template[end]] = get(counts, template[end], 0) + 1

    return sort(collect(counts), by = x -> x[2])
end

diffMostLeastCommon(countsSorted) = countsSorted[end][2] - countsSorted[begin][2]

#! Part 1

template10Counts = polymerizationCountsSorted(template, 10)
diffMostLeastCommon(template10Counts)

#! Part 2

template40Counts = polymerizationCountsSorted(template, 40)
diffMostLeastCommon(template40Counts)

max_stars_repo_path: Day 14/main.jl

max_issues_repo_path: Day 14/main.jl

max_forks_repo_path: Day 14/main.jl


text:
using DifferentialEquations
using Plots

function SIR(du,u,p,t)
    s,i,r = u
    b,k = p
    du[1] = -b*s*i
    du[2] = b*s*i - k*i
    du[3] = k*i
end

#initian conditions
p = [0.5/7900000.0,0.33]
u0 = [7900000.0,10.0,0.0]
# u0 = [1.,1.27e-6,0.]
tspan = (0.0,140.0)

#solve
sir = ODEProblem(SIR,u0,tspan,p)
sol_sir = solve(sir)

#plot
plot(sol_sir)

max_stars_repo_name: Ved-Mahajan/JustForFun

max_issues_repo_name: Ved-Mahajan/JustForFun

max_forks_repo_name: Ved-Mahajan/JustForFun


text:
using DiffEqBase: AbstractODESolution

function shoot((ws, angle, w))
  t = TrebuchetState(;wind_speed=ws, release_angle=angle, weight=w)
  if w <= 0
      @warn "negative weight"
      return (t, sling_end(t).x)
  end
  simulate(t)
  (t, endDist(t))
end

shoot(x...) = shoot(x)

ft2m(x) = x*0.3048
lb2kg(x) = x*0.45359237

Lengths(::Val{:ft}, args...) = Lengths(ft2m.(args)...)
Masses(::Val{:lb}, args...) = Masses(lb2kg.(args)...)
Solution() = Solution([], [], [], [], [], [], [], -1, -1)

function TrebuchetState(;wind_speed=1.0, release_angle=deg2rad(45), weight=98.09)
    release_angle = asin(sin(release_angle))
    l = Lengths(Val{:ft}(), 5.0, 6.792, 1.75, 2.0, 6.833, 2.727, 0.1245)
    m = Masses(Val{:lb}(), promote(weight, 0.328, 10.65)...)
    c = Constants(promote(wind_speed, 1.0, 1.0, 9.80665, release_angle)...)
    t = TrebuchetState(l, m, c, 60.0)
    t.i = Inertias(promote(lb2kg(1.0) |> ft2m |> ft2m , t.i.ia)...)
    return t
end

Base.display(s::Solution) = Base.show(stdin, s)
Base.show(io::IO, ::MIME"text/plain", s::Solution) = Base.show(io, s)
Base.show(io::IO, s::Solution) = println(io, "Solution($(length(s.WeightCG)))")

approx_less(a, b, atol=1e-2) = isapprox(a, b, atol=atol) || a < b

Base.:+(a::Vec, b::Vec) = Vec(a.x + b.x, a.y + b.y)
Base.:-(a::Vec, b::Vec) = Vec(a.x - b.x, a.y - b.y)
Base.:/(a::Vec, b::Number) = Vec(a.x/b, a.y/b)
polar(r, θ) = Vec(r*cos(θ), r*sin(θ))

function ∠(a::Vec)
    θ = atan(a.y/a.x)
    a.x >= 0 && return θ
    a.y >= 0 && return π + θ
    -π + θ
end

function sling_end(t)
   l, a = t.l, t.a
   b, e = l.b, l.e
   aq, sq = a.aq, a.sq

   p = polar(b, π/2 + aq)
   q = polar(e, π/2 + aq + sq)
   p + q
end

function sling_velocity(t::TrebuchetState)
    l, a, aw_ = t.l, t.a, t.aw
    e, b = l.e, l.b
    aq, sq = a.aq, a.sq
    aw, sw = aw_.aw, aw_.sw
    sling_velocity(t, b, e, aq, sq, aw, sw)
end

sling_velocity(t::TrebuchetState, u::Array) = sling_velocity(t, t.l.b, t.l.e, u[1], u[3], u[4], u[6])
sling_velocity(t::TrebuchetState, b, e, aq, sq, aw, sw) =
  polar(-e*(sw + aw),  sq + aq) + polar(-b*aw, aq)

projectile_angle(t::TrebuchetState, u::Array) =
  ∠(sling_velocity(t, u))

endTime(t::TrebuchetState) = t.sol.Time[end]
endDist(t::TrebuchetState) = t.sol.Projectile[end][1]


text:
"""
    inout_map_property(𝑃::Predicate,
                       partition::AbstractVector{<:AbstractVector{Int}},
                       blocks::AbstractVector{Int},
                       n::Int
                      )

Map a property to the dimensions of analyzed blocks.

### Input

- `𝑃`         -- property
- `partition` -- block partition; elements are start and end indices of a block
- `blocks`    -- list of all output block indices in the partition
- `n`         -- total number of input dimensions

### Output

A new property of reduced dimension, if necessary.

### Notes

The return type depends on the ambient dimension (`n`), the block indices in the
partition and on the type of property:

- The property `𝑃` is returned unchanged whenever `n` matches the dimensions
  in `blocks`.
- If the property is a `HalfSpace` (resp. `HPolyhedron`), the function `project`
  from `LazySets` returns a new property with a `HalfSpace` (resp. `HPolyhedron`)
  in the reduced dimensions, according to `blocks`.
- Otherwise, the dimensional reduction is implemented via a (lazy) `LinearMap`.
"""
function inout_map_property(𝑃::Predicate,
                            partition::AbstractVector{<:AbstractVector{Int}},
                            blocks::AbstractVector{Int},
                            n::Int
                           )

    # create a sorted list of all dimensions in `blocks` => available variables for projection
    proj = vcat(partition[blocks]...)

    # no change in the dimension => return the original property
    length(proj) == n && return 𝑃

    return project(𝑃, proj)
end


text:
# # Dihedral symmetry of the Robinson form

# Based on version from SumOfSquares by Benoît Legat
# **Adapted from**: Example 5.4 of [GP04]
#
# [GP04] Gatermann, Karin and Parrilo, Pablo A.
# *Symmetry groups, semidefinite programs, and sums of squares*.
# Journal of Pure and Applied Algebra 192.1-3 (2004): 95-128.

using LinearAlgebra
using SparseArrays

using SymbolicWedderburn
using SymbolicWedderburn.StarAlgebras
using DynamicPolynomials
using SumOfSquares
using SCS

OPTIMIZER = optimizer_with_attributes(
    SCS.Optimizer,
    "acceleration_lookback" => 10,
    "max_iters" => 3_000,
    "alpha" => 1.2,
    "eps" => 1e-6,
    "linear_solver" => SCS.DirectSolver,
)

@polyvar x y
const robinson_form =
    x^6 + y^6 - x^4 * y^2 - y^4 * x^2 - x^4 - y^4 - x^2 - y^2 + 3x^2 * y^2 + 1

# The Robinson form is invariant under the following action of the Dihedral group on monomials:
# The action of each element of the groups is to map the variables `x, y` to:
#
# | id | rotation | reflection |
# |----|----------|------------|
# | 0  | x, y     | y, x       |
# | 1  | -y, x    | -x, y      |
# | 2  | -x, -y   | -y, -x     |
# | 3  | y, -x    | x, -y      |

# implementation of the dihedral group:
include(joinpath(@__DIR__, "dihedral.jl"))
# definitions of general actions on polynomials:
include(joinpath(@__DIR__, "action_polynomials.jl"))

struct DihedralAction <: OnMonomials end

SymbolicWedderburn.coeff_type(::DihedralAction) = Float64
function SymbolicWedderburn.action(
    ::DihedralAction,
    el::DihedralElement,
    mono::AbstractMonomial,
)
    x, y = variables(mono)
    if iseven(el.reflection + el.id)
        var_x, var_y = x, y
    else
        var_x, var_y = y, x
    end
    sign_x = 1 <= el.id <= 2 ? -1 : 1
    sign_y = 2 <= el.id ? -1 : 1
    return mono([x, y] => [sign_x * var_x, sign_y * var_y])
end

G = DihedralGroup(4)
for g in G
    @assert SymbolicWedderburn.action(DihedralAction(), g, robinson_form) ==
            robinson_form
end

# A simple Gramm-matrix like formulation of SOS problem:
# * a single large pdd constraint and
# * one linear constraint for each monomial in the monomial basis of f.
no_symmetry = let f = robinson_form
    stats = @timed let m = Model(OPTIMIZER)
        @variable m t
        @objective m Max t
        @constraint m f - t in SOSCone()
        m
    end
    m = stats.value
    optimize!(m)
    @info (m,) termination_status(m) objective_value(m) creation_time =
        (stats.time) solve_time(m)

    @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2)

    (
        status = termination_status(m),
        objective = objective_value(m),
        creation_t = stats.time,
        solve_t = solve_time(m),
    )
end

# Let's decompose basis into semi-simple summands and prepare the approximately
# reduced optimization problem. Since this is an isomorphism we reduce the
# number of variables from N² to n₁² + … + nₖ², where Σᵢ nᵢ = N. That is instead
# of one psd constraint of size N we have k constraints of sizes n₁,… nₖ.
# The number of linear constraints remains the same as before.
semisimple_dec = let f = robinson_form, G = DihedralGroup(4)
    basis = monomials([x, y], 0:(DynamicPolynomials.maxdegree(f)÷2))

    stats = @timed begin
        sa_basis, symmetry_adaptation_time, = @timed let
            ssimple_basis =
                SymbolicWedderburn.symmetry_adapted_basis(
                    Rational{Int},
                    G,
                    DihedralAction(),
                    basis,
                    semisimple = true,
                )
            sp = sparse.(ssimple_basis)
            droptol!.(sp, 1e-12)
        end

        m = SOSModel(OPTIMIZER)

        let basis = basis, R = sa_basis
            @variable m t
            @objective m Max t

            soses = @variable m [r in R] SOSPoly(
                FixedPolynomialBasis(r * basis),
            )
            @constraint m f - t == sum(soses)
        end
        m
    end
    m = stats.value
    optimize!(m)

    @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2)

    @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time =
        (stats.time) symmetry_adaptation_time
    (
        status = termination_status(m),
        objective = objective_value(m),
        creation_t = stats.time,
        solve_t = solve_time(m),
    )
end

# Now given that f must be constant on orbits of monomials under the action of G
# let's exploit that to reduce the number of cnstraints (one for each orbit).
# Here we keep the single large psd constraint of size N×N.
include(joinpath(@__DIR__, "util.jl"))

orbit_dec = let f = robinson_form, G = DihedralGroup(4), T = Float64
    basis_monoms = monomials([x, y], 0:DynamicPolynomials.maxdegree(f))

    invariant_vs, symmetry_adaptation_time, = @timed let G = G
        tblG = SymbolicWedderburn.Characters.CharacterTable(Rational{Int}, G)
        invariant_vectors(
            tblG,
            DihedralAction(),
            StarAlgebras.Basis{UInt32}(basis_monoms),
        )
    end

    M =
        let basis_full = StarAlgebras.Basis{UInt32}(basis_monoms),
            basis_half =
                monomials([x, y], 0:(DynamicPolynomials.maxdegree(f)÷2))

            [basis_full[x*y] for x in basis_half, y in basis_half]
        end

    stats = @timed let m = JuMP.Model(OPTIMIZER)
        JuMP.@variable m t
        JuMP.@objective m Max t
        P = JuMP.@variable(m, [1:size(M, 1), 1:size(M, 1)], PSD)

        # preallocating
        M_orb = similar(M, T)

        C = DynamicPolynomials.coefficients(f - t, basis_monoms)

        for iv in invariant_vs
            c = dot(C, iv)

            # invariant_constraint! is defined locally in util.jl
            # it essentially amounts to averaging constraint matrices defined by M
            # with weights provided by iv and storing the result in (dense) matrix M_orb
            M_orb = invariant_constraint!(M_orb, M, iv)

            #=
            if !iszero(M_orb)
                monoms = basis_monoms[unique!(M[map(!iszero, M_orb)])]
                @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c
                @show sparse(M_orb)
            else
                monoms = basis_monoms[unique!(M[SparseArrays.nonzeroinds(iv)])]
                @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c
                @info "M_orb is zero!"
            end
            =#

            JuMP.@constraint m dot(M_orb, P) == c
        end
        m
    end

    m = stats.value
    optimize!(m)

    @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2)

    @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time =
        (stats.time) symmetry_adaptation_time
    (
        status = termination_status(m),
        objective = objective_value(m),
        creation_t = stats.time,
        solve_t = solve_time(m),
    )
end

wedderburn_dec = let f = robinson_form, G = DihedralGroup(4), T = Float64
    basis_half = monomials([x, y], 0:(DynamicPolynomials.maxdegree(f)÷2))

    wedderburn, symmetry_adaptation_time, = @timed WedderburnDecomposition(
        T,
        G,
        DihedralAction(),
        monomials([x, y], 0:DynamicPolynomials.maxdegree(f)),
        basis_half,
    )

    M = [
        SymbolicWedderburn.basis(wedderburn)[x*y] for x in basis_half,
        y in basis_half
    ]

    stats = @timed let m = JuMP.Model(OPTIMIZER)
        JuMP.@variable m t
        JuMP.@objective m Max t
        psds = [
            JuMP.@variable(m, [1:d, 1:d] in PSDCone()) for
            d in size.(direct_summands(wedderburn), 1)
        ]

        # preallocating
        Mπs = zeros.(T, size.(psds))
        M_orb = similar(M, T)
        tmps = SymbolicWedderburn._tmps(wedderburn)

        C = DynamicPolynomials.coefficients(
            f - t,
            SymbolicWedderburn.basis(wedderburn),
        )

        for iv in invariant_vectors(wedderburn)
            c = dot(C, iv)

            # invariant_constraint! is defined locally in util.jl
            # it essentially amounts to averaging constraint matrices defined by M
            # with weights provided by iv and storing the result in (dense) matrix M_orb
            M_orb = invariant_constraint!(M_orb, M, iv)

            #=
            if !iszero(M_orb)
                monoms = SymbolicWedderburn.basis(wedderburn)[unique!(M[map(!iszero, M_orb)])]
                @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c
                @show sparse(M_orb)
            else
                monoms = SymbolicWedderburn.basis(wedderburn)[unique!(M[SparseArrays.nonzeroinds(iv)])]
                @info "Coefficient of the $monoms: $([DynamicPolynomials.coefficient(f-t, m) for m in monoms])" dot_C_v=c
                @info "M_orb is zero!"
            end
            =#

            Mπs = SymbolicWedderburn.diagonalize!(Mπs, M_orb, wedderburn, tmps)

            JuMP.@constraint m sum(
                dot(Mπ, Pπ) for (Mπ, Pπ) in zip(Mπs, psds) if !iszero(Mπ)
            ) == c
        end
        m
    end
    m = stats.value
    optimize!(m)

    @info value(m[:t])
    @assert isapprox(value(m[:t]), -3825 / 4096, rtol = 1e-2)

    @info (m,) termination_status(m) objective_value(m) solve_time(m) creation_time =
        (stats.time)
    (
        status = termination_status(m),
        objective = objective_value(m),
        creation_t = stats.time,
        solve_t = solve_time(m),
    )
end


text:
using HorizonSideRobots

function cross_border_cells_marking(r::Robot)
    back_track = []
    vertical_shift = 0
    horizontal_shift = 0
    while ((isborder(r, Sud)) && (isborder(r, West))) == false
        push!(back_track, show_path_after_moving(r, West))
        push!(back_track, show_path_after_moving(r, Sud))
    end

    n = length(back_track)
    for i in 1:n
        if (i % 2) == false
            vertical_shift += back_track[i]
        else
            horizontal_shift += back_track[i]
        end
    end

    for i in 1:horizontal_shift
        move!(r, Ost)
    end
    putmarker!(r)
    while isborder(r, Ost) == false
        move!(r, Ost)
    end
    for i in 1:vertical_shift
        move!(r, Nord)
    end
    putmarker!(r)
    while isborder(r, Nord) == false
        move!(r, Nord)
    end

    while isborder(r, West) == false
        move!(r, West)
    end

    for i in 1:horizontal_shift
        move!(r, Ost)
    end
    putmarker!(r)
    while isborder(r, West) == false
        move!(r, West)
    end
    while isborder(r, Sud) == false
        move!(r, Sud)
    end
    for i in 1:vertical_shift
        move!(r, Nord)
    end
    putmarker!(r)
    while isborder(r, Sud) == false
        move!(r, Sud)
    end

    while n > 0
        if n % 2 == 1
            side = Ost
        else
            side = Nord
        end
        move_n(r, side, back_track[n])
        n -= 1
    end
end

function move_n(r::Robot, side::HorizonSide, n::Int)
    for i in 1:n
        move!(r, side)
    end
end

function show_path_after_moving(r::Robot, side::HorizonSide)
    n = 0
    while isborder(r, side)==false
        move!(r, side)
        n = n + 1
    end
    return n
end


text:
module VTUFileHandler

using Base64, CodecZlib
using XMLParser

include(joinpath(".","VTUFileHandler","defs.jl"))

"""
    VTUHeader(::Type{T},input::Vector{UInt8}) where {T<:Union{UInt32,UInt64}}

Computes the VTUHeader based on the headertype and a Base64 decoded input data array.

# Constructor
- `::Type{T}`: headertype, either UInt32 or UInt64
- `input::Vector{UInt8}`: input data

# Fields
- `num_blocks::T` : number of blocks
- `blocksize::T` : size of blocks
- `last_blocksize::T` : size of last block (can be different)
- `compressed_blocksizes::T` : size of compressed blocks
"""
struct VTUHeader{T<:Union{UInt32,UInt64}}
	num_blocks::T
	blocksize::T
	last_blocksize::T
	compressed_blocksizes::Vector{T}
	function VTUHeader(::Type{T}) where {T<:Union{UInt32,UInt64}}
		return new{T}(zero(T),zero(T),zero(T),T[])
	end
	function VTUHeader(::Type{T},input::Vector{UInt8}) where {T<:Union{UInt32,UInt64}}
		dat = reinterpret(T,input)
		num_blocks = dat[1]
		blocksize = dat[2]
		last_blocksize = dat[3]
		compressed_blocksizes=Vector{T}(undef,num_blocks)
		for (j,i) = enumerate(4:length(dat))
			compressed_blocksizes[j] = dat[i]
		end
		return new{T}(num_blocks,blocksize,last_blocksize,compressed_blocksizes)
	end
	function VTUHeader(a::T,b::T,c::T,d::Vector{T}) where {T<:Union{UInt32,UInt64}}
		return new{T}(a,b,c,d)
	end
end

include(joinpath(".","VTUFileHandler","vtuheader_utils.jl"))
include(joinpath(".","VTUFileHandler","io.jl"))

"""
    VTUDataField(x::Vector{T}) where {T} = new{T}(x)

Container for VTU field data.

# Constructor
- `x::Vector{T}`: Vector with VTU field data

# Fields
- `dat::Vector{T}` : VTU field data
"""
struct VTUDataField{T}
	dat::Vector{T}
	VTUDataField(x::Vector{T}) where {T} = new{T}(x)
	VTUDataField(::Type{T}) where {T} = new{T}()
	VTUDataField(x::Base.ReinterpretArray{T, A, B, Vector{C}, D}) where {T,A,B,C,D} = new{T}(x)
end

"""
    VTUData(dataarrays::Vector{XMLElement},appendeddata::Vector{XMLElement},headertype::Union{Type{UInt32},Type{UInt64}},offsets::Vector{Int},compressed_dat::Bool)

Container for VTU data.

# Constructor
- `dataarrays::Vector{XMLElement}`: Vector with all XML elements with tag `DataArray` of VTU file
- `appendeddata::Vector{XMLElement}`: Vector with all XML elements with tag `AppendedData`
- `headertype::Union{Type{UInt32},Type{UInt64}}` : Type of VTU header
- `offsets::Vector{Int}` : Offset of each field data in the compressed appended data
- `compressed_dat::Bool` : True if data is compressed

# Fields
- `names::Vector{String}` : Name of each data field
- `header::Vector{VTUHeader}` : Vector with [`VTUHeader`](@ref)s.
- `data::Vector{VTUDataField}` : Vector with [`VTUDataField`](@ref)s
- `interp_data::Vector{VTUDataField{Float64}}` : Vector with [`VTUDataField`](@ref)s
- `idat::Vector{Int}` : Vector indexing the `data`-fields onto which the math operators should be applied

For a [`VTUDataField`](@ref) to appear in `data`, the appropriate keyword has to be added via [`add_uncompress_keywords`](@ref) or [`add_interpolation_keywords`](@ref).
For a [`VTUDataField`](@ref) to appear in `interp_data` the appropriate keyword has to be added via [`add_interpolation_keywords`](@ref).
For more information, see [`VTUKeyWords`](@ref).
"""
struct VTUData
	names::Vector{String}
	header::Vector{VTUHeader}
	data::Vector{VTUDataField}
	interp_data::Vector{VTUDataField{Float64}}
	idat::Vector{Int}
	VTUData(names,header,data,interp_data,idat) = new(names,header,data,interp_data,idat)
	function VTUData(dataarrays::Vector{XMLElement},appendeddata::Vector{XMLElement},headertype::Union{Type{UInt32},Type{UInt64}},offsets::Vector{Int},compressed_dat::Bool)
		names = Vector{String}()
		data = Vector{VTUDataField}()
		header = Vector{VTUHeader}()
		for (i,el) in enumerate(dataarrays)
			_type = replace(getAttribute(el,"type"),"\""=>"")
			type = eval(Meta.parse(_type))
			_format = getAttribute(el,"format")
			_name = getAttribute(el,"Name")
			if findfirst(x->x==_name,vtukeywords.interpolation_keywords) != nothing	|| findfirst(x->x==_name,vtukeywords.uncompress_keywords) != nothing
				if _format == "appended" || _format == "\"appended\""
					vtuheader = readappendeddata!(data,i,appendeddata,offsets,type,headertype,compressed_dat)
				else
					vtuheader = readdataarray!(data,el,type,headertype,compressed_dat)
				end
				push!(names,_name)
				push!(header,vtuheader)
			end
		end
		idat = findall(map(y->findfirst(x->x==y,vtukeywords.interpolation_keywords)!=nothing,names))
		interp_data = data[idat]
		return new(names,header,data,interp_data,idat)
	end
end

include(joinpath(".","VTUFileHandler","vtudata_utils.jl"))
include(joinpath(".","VTUFileHandler","vtudata_math.jl"))

function deletefieldata!(xmlroot)
	dataarrays = getElements(xmlroot,"DataArray")
	fielddata = getElements(xmlroot,"FieldData")[1]
	_offset = parse(Int,replace(getAttribute(dataarrays[length(fielddata.content)+1],"offset"),"\""=>""))
	empty!(fielddata.content)
	els = getElements(xmlroot,"DataArray")
	for el in els
		setAttribute(el,"offset",string(parse(Int,replace(getAttribute(el,"offset"),"\""=>""))-_offset))
	end
	appendeddata = getElements(xmlroot,"AppendedData")
	appendeddata[1].content[1] = String(deleteat!(collect(appendeddata[1].content[1]),2:_offset+1))
end

"""
    VTUFile(name::String)

Loads a VTU file.
Don't forget to set the proper fieldnames via `set_uncompress_keywords` and `set_interpolation_keywords`.

# Constructor
- `name::String`: path to vtu file

# Fields
- `name::String`: path to vtu file; destination for file writing
- `xmlroot::XMLElement`: VTU file in XML represantation
- `dataarrays::Vector{XMLElement}`: Vector with all XML elements with tag `DataArray` of VTU file
- `headertype::Union{Type{UInt32},Type{UInt64}}`: type of header
- `offsets::Vector{Int}`: Offset of each field data in the compressed appended data
- `data::VTUData`: Conatainer with  [`VTUData`](@ref)
- `compressed_dat::Bool` : True if data is compressed

# Example
```julia
set_uncompress_keywords("temperature","points")
set_interpolation_keywords("temperature")
vtufile = VTUFile("./path-to-vtu/example.vtu");
```
"""
mutable struct VTUFile
	name::String
	xmlroot::XMLElement
	dataarrays::Vector{XMLElement}
	appendeddata::Vector{XMLElement}
	headertype::Union{Type{UInt32},Type{UInt64}}
	offsets::Vector{Int}
	data::VTUData
	compressed_dat::Bool
	VTUFile(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data) = new(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data,true)
	VTUFile(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data,compr_dat) = new(name,xmlroot,dataarrays,appendeddata,headertype,offsets,data,compr_dat)
	function VTUFile(name::String)
		state = IOState(name);
		xmlroot = readXMLElement(state);
		if !isempty(getElements(xmlroot,"FieldData"))
			deletefieldata!(xmlroot)
		end
		dataarrays = getElements(xmlroot,"DataArray")
		appendeddata = getElements(xmlroot,"AppendedData")
		vtkfile = getElements(xmlroot,"VTKFile")
		@assert length(vtkfile) == 1
		compressed_dat = false
		if hasAttributekey(vtkfile[1],"compressor")
			compressed_dat = true
		end
		#headertype = eval(Symbol(getAttribute(first(vtkfile),"header_type")))
		attr = replace(getAttribute(first(vtkfile),"header_type"),"\""=>"")
		headertype = eval(Meta.parse(attr))
		offsets = Vector{Int}()
		for el in dataarrays
			if hasAttributekey(el,"offset")
				_offset = replace(getAttribute(el,"offset"),"\""=>"")
				offset = parse(Int,_offset)
				push!(offsets,offset)
			else
				push!(offsets,0)
			end
		end
		retval = new(name,xmlroot,dataarrays,appendeddata,headertype,offsets,VTUData(dataarrays,appendeddata,headertype,offsets,compressed_dat),compressed_dat)
		return retval
	end
end

include(joinpath(".","VTUFileHandler","utils.jl"))
include(joinpath(".","VTUFileHandler","vtufile_math.jl"))
include(joinpath(".","VTUFileHandler","globaltolocal.jl"))

export VTUHeader, VTUFile, set_uncompress_keywords, add_uncompress_keywords, set_interpolation_keywords, add_interpolation_keywords, add!, minus!, pow!, div!, min!, max!

end #module VTUHandler

max_stars_repo_name: baxmittens/VTUFileHandler

max_issues_repo_name: baxmittens/VTUFileHandler

max_forks_repo_name: baxmittens/VTUFileHandler


text:
@testset "spin_half.jl" begin

    @testset "nspins" begin
        for N in 0:10
            @test nspins(2^N) == N
            @test nspins(collect(1:2^N)) == N
            @test nspins(zeros(2^N, 2^N)) == N
        end
        @test_throws ArgumentError nspins(collect(1:3))
        @test_throws ArgumentError nspins(zeros(3,3))
        @test_throws ArgumentError nspins(zeros(2,4))
    end

    @testset "single_spin_op" begin
        @test single_spin_op(σz, 1, 2) ≈ Diagonal([1,1,-1,-1])
        @test single_spin_op(σz, 2, 2) ≈ Diagonal([1,-1,1,-1])
        @test single_spin_op(σz, 2, 3) ≈ kron(Diagonal([1,-1,1,-1]), 𝟙)

        @test_throws DomainError single_spin_op(σz, 2, 1) # outside range
        @test_throws DomainError single_spin_op(σz, -1, 3) # bogus site - not in range [1, N]
        @test_throws DomainError single_spin_op(σz, 0, 3) # bogus site - not in range [1, N]
        @test_throws DomainError correlator(𝟙, 1, 2, -3) # negative length

        @test_throws MethodError single_spin_op(σz, 1.5, 3) # bogus site - not integer
        @test_throws MethodError correlator(𝟙, 1, 2, 2.5) # bogus length
    end

    @testset "correlator" begin
        @test correlator(σx, 1, 2, 2) ≈ σx ⊗ σx
        @test correlator(σz, 1, 2, 3) ≈ σz ⊗ σz ⊗ 𝟙
        @test correlator(σz, 2, 1, 3) ≈ σz ⊗ σz ⊗ 𝟙 # symmetry

        @test_throws DomainError correlator(𝟙, 1, 1, 2) # operators on same site
        @test_throws DomainError correlator(𝟙, 1, 3, 2) # operator outside
        @test_throws DomainError correlator(𝟙, -1, 2, 3) # negative site
        @test_throws DomainError correlator(𝟙, 1, 2, -3) # negative length

        @test_throws MethodError correlator(𝟙, 1, 1.5, 3) # bogus site
        @test_throws MethodError correlator(𝟙, 1, 2, 2.5) # bogus length
    end

    @test op_list(σy, 2) ≈ [σy ⊗ 𝟙, 𝟙 ⊗ σy]
end


text:
"""
    cl2(x::Float64)::Float64

Returns the value of the Clausen function ``\\operatorname{Cl}_2(x)``
for a real angle ``x`` of type `Float64`.  The Clausen function is
defined as

```math
\\operatorname{Cl}_2(x) = -\\int_0^x \\log|2\\sin(t/2)| dt
```

Author: Alexander Voigt

License: MIT

# Example
```jldoctest; setup = :(using ClausenFunctions), output = false
julia> cl2(1.0)
1.0139591323607684
```
"""
function cl2(x::Float64)::Float64
    pi28 = pi*pi/8.0

    (x, sgn) = range_reduce_even(x)

    if x == zero(x) || x == pi
        return zero(x)
    end

    if x < 0.5*pi
        P = (1.3888888888888889e-02, -4.3286930203743071e-04,
             3.2779814789973427e-06, -3.6001540369575084e-09)
        Q = (1.0000000000000000e+00, -3.6166589746694121e-02,
             3.6015827281202639e-04, -8.3646182842184428e-07)
        y = x*x
        y2 = y*y
        p = P[1] + y * P[2] + y2 * (P[3] + y * P[4])
        q = Q[1] + y * Q[2] + y2 * (Q[3] + y * Q[4])
        sgn*x*(1.0 - log(x) + y*p/q)
    else
        P = (6.4005702446195512e-01, -2.0641655351338783e-01,
             2.4175305223497718e-02, -1.2355955287855728e-03,
             2.5649833551291124e-05, -1.4783829128773320e-07)
        Q = (1.0000000000000000e+00, -2.5299102015666356e-01,
             2.2148751048467057e-02, -7.8183920462457496e-04,
             9.5432542196310670e-06, -1.8184302880448247e-08)
        y = pi - x
        z = y*y - pi28
        z2 = z*z
        z4 = z2*z2
        p = P[1] + z * P[2] + z2 * (P[3] + z * P[4]) + z4 * (P[5] + z * P[6])
        q = Q[1] + z * Q[2] + z2 * (Q[3] + z * Q[4]) + z4 * (Q[5] + z * Q[6])
        sgn*y*p/q
    end
end


text:
#=

Log linear optimal growth model, with log utility, CD production and
multiplicative lognormal shock, so that

    y = f(k, z) = z k^alpha

with z ~ LN(mu, s).

=#

alpha = 0.4
beta = 0.96
mu = 0
s = 0.1

ab = alpha * beta
c1 = log(1 - ab) / (1 - beta)
c2 = (mu + alpha * log(ab)) / (1 - alpha)
c3 = 1 / (1 - beta)
c4 = 1 / (1 - ab)

# Utility
u(c) = log(c)

u_prime(c) = 1 / c

# Deterministic part of production function
f(k) = k^alpha

f_prime(k) = alpha * k^(alpha - 1)

# True optimal policy
c_star(y) = (1 - alpha * beta) * y

# True value function
v_star(y) = c1 + c2 * (c3 - c4) + c4 * log(y)


text:
using POMDPs
using POMDPModels
using POMDPModelTools
using StaticArrays
using JLD2, FileIO
using Random
using Logging
using Distributions
using Statistics
using GridInterpolations
using LocalFunctionApproximation
using LocalApproximationValueIteration
using CMSSPs

rng = MersenneTwister(1234)

flight_policy_name = "dreamr-uf-params111.jld2"

arguments = ["./paramsets/scale-1.toml","./paramsets/simtime-1.toml",
        "./paramsets/cost-1.toml", "5"]

scale_file = arguments[1]
simtime_file = arguments[2]
cost_file = arguments[3]
num_trials = parse(Int64, arguments[4])

# First, parse parameter files to get filenames and construct params object
params = parse_params(scale_file=scale_file, simtime_file=simtime_file, cost_file=cost_file)

uf_mdp = get_uf_mdp(MultiRotorUAVState, MultiRotorUAVAction, params)

flight_policy = load_localapproxvi_policy_from_jld2(flight_policy_name)

costs = Vector{Float64}(undef, 0)

for i = 1 : num_trials

    curr_state = generate_start_state(MultiRotorUAVState, params, rng)
    goal_state = get_state_at_rest(MultiRotorUAVState, Point(0.0, 0.0))

    tot_cost = 0.0

    @show value(flight_policy, get_relative_state(uf_mdp, curr_state, goal_state))

    while true

        # @show curr_state
        inf_tp_dist = TPDistribution([typemax(Int64)], [1.0])

        a = get_best_intramodal_action(flight_policy, 0, inf_tp_dist, curr_state, goal_state)
        # @show a
        # readline()
        (next_state, reward) = generate_sr(uf_mdp, curr_state, a.action, rng)

        curr_state = next_state

        tot_cost += -1.0*reward

        if isterminal(uf_mdp, get_relative_state(uf_mdp, curr_state, goal_state))
            @info "Reached!"
            break
        end
    end

    push!(costs, tot_cost)
end

@show costs


text:
# -- PRIVATE FUNCTIONS NOT EXPORTED -------------------------------------------------------------- #
function _reorder_logic_with_callback(symbol_array::Array{String,1}; callback::Union{Function,Nothing} = nothing)

    # if there is callback logic -> use it ...
    if (isnothing(callback) == false)
        return callback(symbol_array)
    end

    # otherwise - alphabetical
    return sort(symbol_array)
end
# ------------------------------------------------------------------------------------------------ #

"""
    reorder_molecular_symbol_array(symbol_array::Array{String,1};
          callback::Union{Function,Nothing} = nothing)::Array{String,1}

Reorder species symbol names alphabetically (default implementation) unless a callback function is passed, in which case sorting is done based on the user-defined routine.

Input arguments:
`symbol_array::Array{String,1}` - array holding symbol names for all molecular species in the network.
`callback::Union{Function,Nothing}` - user-defined function to reorder species symbol names (optional).

Output arguments:
`symbol_array::Array{String,1}` - sorted molecular species symbol array.

"""
function reorder_molecular_symbol_array(symbol_array::Array{String,1}; callback::Union{Function,Nothing} = nothing)::Array{String,1}
    return _reorder_logic_with_callback(symbol_array; callback = callback)
end

"""
    reorder_reaction_symbol_array(symbol_array::Array{String,1};
          callback::Union{Function,Nothing} = nothing)::Array{String,1}

Reorder reaction symbol names alphabetically (default implementation) unless a callback function is passed, in which case sorting is done based on the user-defined routine.

Input arguments:
`symbol_array::Array{String,1}` - array holding symbol names for all reactions in the network.
`callback::Union{Function,Nothing}` - user-defined function to reorder reaction names (optional).

Output arguments:
`symbol_array::Array{String,1}` - sorted reaction symbol array.

"""
function reorder_reaction_symbol_array(symbol_array::Array{String,1}; callback::Union{Function,Nothing} = nothing)::Array{String,1}
    return return _reorder_logic_with_callback(symbol_array; callback = callback)
end

max_stars_repo_name: varnerlab/CellFreeModelGenerationKit

max_issues_repo_name: varnerlab/CellFreeModelGenerationKit

max_forks_repo_name: varnerlab/CellFreeModelGenerationKit


text:
"""
    convergence_rates(problem, solution, method, dts; kwargs...)

Compute the errors rates of `method` on `problem` by comparing to `solution`
on the set of `dt` values in `dts`. Extra `kwargs` are passed to `solve`

`solution` should be a function with a method `solution(u0, p, t)`.
"""
function convergence_errors(prob, sol, method, dts; kwargs...)
  errs = map(dts) do dt
       # copy the problem so we don't mutate u0
      prob_copy = deepcopy(prob)
      u = solve(prob_copy, method; dt=dt, kwargs...)
      norm(u .- sol(prob.u0, prob.p, prob.tspan[end]))
  end
  return errs
end

"""
  convergence_order(problem, solution, method, dts; kwargs...)

Estimate the order of the rate of convergence of `method` on `problem` by comparing to
`solution` the set of `dt` values in `dts`.

`solution` should be a function with a method `solution(u0, p, t)`.
"""
function convergence_order(prob, sol, method, dts; kwargs...)
  errs = convergence_errors(prob, sol, method, dts; kwargs...)
  # find slope coefficient in log scale
  _,order_est = hcat(ones(length(dts)), log2.(dts)) \ log2.(errs)
  return order_est
end

"""
    DirectSolver

A linear solver which forms the full matrix of a linear operator and its LU factorization.
"""
struct DirectSolver end

DirectSolver(args...) = DirectSolver()

function (::DirectSolver)(x,A,b,matrix_updated; kwargs...)
  n = length(x)
  M = mapslices(y -> mul!(similar(y), A, y), Matrix{eltype(x)}(I,n,n), dims=1)
  x .= M \ b
end


text:
"""
    DD.set_parameters!(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
        θ°,
        var_p_names,
        var_p_aux_names,
        θ_local_names,
        θ_local_aux_names,
        θ_local_obs_ind,
        critical_change = is_critical_update(
            PP, θ_local_aux_names, θ_local_obs_ind
        ),
    )

Set parameters of Guided proposals `PP°` in an MCMC setting.

...
# Arguments
- `PP`: containers indicating how `PP°` should look like before `θ°` update
- `θ°`: a vector with new parameters to be potentially set inside `PP°`
- `var_p_names`: list of all variable parameter names in `PP`
- `var_p_aux_names`: list of all variable parameter names of auxiliary laws in
`PP`
- `θ_local_names`: list of pairs with relevant `θ°` entries for the target law
in a format (idx-of-param-in-θ°::Int64, param-name-in-law::Symbol)
- `θ_local_aux_names`: list of pairs with relevant `θ°` entries for the aux law
in a format (idx-of-param-in-θ°::Int64, param-name-in-law::Symbol)
- `θ_local_obs_ind`: list of pairs with relevant `θ°` entries for the
observations in a format (idx-of-param-in-θ°::Int64, idx-in-obs::Int64)
- `critical_change`: boolean for whether `θ°` update alters any critical
parameters prompting for recomputation of the guiding term
...
"""
function DD.set_parameters!(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
        θ°,
        var_p_names,
        var_p_aux_names,
        θ_local_names,
        θ_local_aux_names,
        θ_local_obs_ind,
        critical_change = is_critical_update(
            PP, θ_local_aux_names, θ_local_obs_ind
        ),
    )
    # Resetting the state from the previous update
    # --------------------------------------------
    # first, equalize observations
    equalize_obs_params!(PP, PP°) && (critical_change = true)

    # second, equalize parameters
    if !DD.same_entries(PP, PP°, var_p_names)
        equalize_law_params!(
            PP, PP°, var_p_names, var_p_aux_names
        ) && (critical_change = true)
    end

    # Setting the new parameters
    # --------------------------
    DD.set_parameters!(
        PP°, θ°, θ_local_names, θ_local_aux_names, θ_local_obs_ind
    )
    critical_change
end

"""
    equalize_obs_params!(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp}
    )

Go through a collection of `GuidProp` in PP° and make sure that all observation
parameters are the same. If not, equalize them and return a `critical_change`
flag indicating that `GuidProp` laws need to be recomputed.
"""
function equalize_obs_params!(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
    )
    critical_change = false
    for i in eachindex(PP, PP°)
        if PP[i].obs.θ != PP°[i].obs.θ
            critical_change = true
            PP°[i].obs.θ .= PP[i].obs.θ
        end
    end
    critical_change
end

"""
    DD.same_entries(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
        entries
    )
Go through the collections of guided proposals `PP` and `PP°` and compare if
they share the same values of fields listed in `entries`.
"""
function DD.same_entries(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
        entries
    )
    DD.same_entries(PP[1].P_target, PP°[1].P_target, entries)
end

"""
    equalize_law_params!(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
        var_p_names,
        var_p_aux_names,
    )

Go through two collections of `GuidProp` and compare their variable parameters.
If any of the parameters in `PP°` ends up being different than in `PP` then
equalize them and if any of them belongs to the auxiliary law then return a
`critical_change` flag.
"""
function equalize_law_params!(
        PP::AbstractArray{<:GuidProp},
        PP°::AbstractArray{<:GuidProp},
        var_p_names,
        var_p_aux_names,
    )
    critical_change = false
    for i in eachindex(PP, PP°)
        DD.set_parameters!(PP[i].P_target, PP°[i].P_target, var_p_names)
        ( # check for a critical change
            critical_change ||
            DD.same_entries(PP[i].P_aux, PP°[i].P_aux, var_p_aux_names[i]) ||
            (critical_change = true)
        )
        DD.set_parameters!(PP[i].P_aux, PP°[i].P_aux, var_p_aux_names[i])
    end
    return critical_change
end

"""
    DD.set_parameters!(
        PP::AbstractArray{<:GuidProp},
        θ°,
        θ_local_names,
        θ_local_aux_names,
        θ_local_obs_ind
    )

Go throught the collection of `GuidProp`s in `PP` and set the new parameters
`θ°` in relevant fields. `θ_local_names` lists which fields in the target
law need changing. `θ_local_aux_names` lists the same for the auxiliary laws
and `θ_local_obs_ind` for the observations.
"""
function DD.set_parameters!(
        PP::AbstractArray{<:GuidProp},
        θ°,
        θ_local_names,
        θ_local_aux_names,
        θ_local_obs_ind
    )
    for i in eachindex(PP)
        DD.set_parameters!(
            PP[i], θ°, θ_local_names, θ_local_aux_names[i], θ_local_obs_ind[i]
        )
    end
end

"""
    DD.set_parameters!(
        P::GuidProp,
        θ°,
        θ_local_names,
        θ_local_aux_names,
        θ_local_obs_ind
    )

Set the new parameters `θ°` in `P`. `θ_local_names` lists which fields in the
target law need changing. `θ_local_aux_names` lists the same for the auxiliary
law and `θ_local_obs_ind` for the observation.
"""
function DD.set_parameters!(
        P::GuidProp,
        θ°,
        θ_local_names,
        θ_local_aux_names,
        θ_local_obs_ind
    )
    DD.set_parameters!(P.P_target, θ°, θ_local_names)
    DD.set_parameters!(P.P_aux, θ°, θ_local_aux_names)
    OBS.set_parameters!(P.obs, θ°, θ_local_obs_ind)
    P.guiding_term_solver.integrator.p = P.P_aux # for safe measure
end

"""
    is_critical_update(
        PP::AbstractArray{<:GuidProp},
        θ_local_aux_names,
        θ_local_obs_ind
    )

Check if the update of parameters that updates the fields of the auxiliary law
listed in `θ_local_aux_names` and the observations listed in `θ_local_obs_ind`
is critical i.e. whether it prompts for recomputation of the guiding term.
"""
function is_critical_update(
        PP::AbstractArray{<:GuidProp},
        θ_local_aux_names,
        θ_local_obs_ind
    )
    for i in eachindex(PP)
        is_critical_update(
            PP[i], θ_local_aux_names[i], θ_local_obs_ind[i]
        ) && return true
    end
    false
end

function is_critical_update(P::GuidProp, θ_local_aux_names, θ_local_obs_ind)
    DD.has_any(P.P_aux, θ_local_aux_names) || length(θ_local_obs_ind) > 0
end

"""
    DD.set_parameters!(P::GuidProp, θ°::Dict)

Convenience parameter setter
"""
function DD.set_parameters!(P::GuidProp, θ°::Dict)
    DD.set_parameters!(P.P_target, θ°)
    DD.set_parameters!(P.P_aux, θ°)
    P.guiding_term_solver.integrator.p = P.P_aux # for safe measure
end

"""
    DD.set_parameters!(PP::AbstractArray{<:GuidProp}, θ°::Dict)

Convenience parameter setter
"""
function DD.set_parameters!(PP::AbstractArray{<:GuidProp}, θ°::Dict)
    for i in eachindex(PP)
        DD.set_parameters!(PP[i], θ°)
    end
end

function set_obs!(P::GuidProp, obs)
    P.obs.obs .= obs
    P.P_aux.vT .= obs
    P.guiding_term_solver.integrator.p = P.P_aux
end


text:
import Images: imresize, normedview
import ImageTransformations: imresize!
export GreyMeanWrapper, GreyChanWrapper, DownsizeWrapper, MultiFrameWrapper, state_img
using PyCall

# --------------------------------
# Greyscale Wrappers
# --------------------------------
mat_like_state(::Type{T}, env) where T = Matrix{T}(size(state(env))[1:2])

abstract type GreyscaleGymWrapper <: AbstractGymWrapper end

"""
`GreyMeanWrapper(env::AbstractGymEnv, ::Type{T} = UInt8, zero_centre::Bool = false)`

For each pixel in `state(env)`, assumed to be encoded as RGB UInt8s in [0,255].
Average the R, G, and B channels, and normalise so that
`state(GreyMeanWrapper(env))` is a Matrix with all values in the interval [0,1]

If `zero_centre` is true subtract 0.5 from each channel, so that all values lie
in the interval [-0.5, 0.5]. Accordingly, `zero_centre` is not supported for
types that don't support negative values, e.g. `UInt8` and Images.jl's
FixedPoint types.
"""
struct GreyMeanWrapper{T} <: GreyscaleGymWrapper
    env::AbstractGymEnv
    zero_centre::Bool
    img::Matrix{T} # initialise on creation, update each `step!`
end

function GreyMeanWrapper(env::AbstractGymEnv, ::Type{T} = UInt8,
  zero_centre::Bool = false) where T
    img = mat_like_state(T, env)
    rgb2greymean!(img, state(env), zero_centre)
    GreyMeanWrapper{T}(env, zero_centre, img)
end

# @generated function rgb2greymean!(wrap_img::Matrix{FPT}, state_img::Array{UT,3}, zero_centre::Bool) where
#   {UT <: Unsigned, FPT <: Normed}
#     @assert sizeof(UT) == sizeof(FPT)
#     # must divide each UInt8 by 3 rather than their sum, to avoid overflow
#     quote
#         wrap_img .= reinterpret.($FPT,
#                       @view(state_img[:,:,1]) .÷ 0x3 .+
#                       @view(state_img[:,:,2]) .÷ 0x3 .+
#                       @view(state_img[:,:,3]) .÷ 0x3
#                     )
#     end
# end

# fast conversion of PyArray{UInt8} to Array{N0f8}
# @benchmark unsafe_wrap(Array, convert(Ptr{N0f8}, pyimg.data), size(pyimg))

"""
Returned array's eltype == input array's eltype == UInt8
"""
function rgb2greymean!(wrap_img::Matrix{UInt8}, state_img::AbstractArray{UInt8,3},
    zero_centre::Bool)
    # must divide each UInt8 by 3 rather than their sum, to avoid overflow
    wrap_img .= (@view(state_img[:,:,1]) .÷ 0x3 .+
                 @view(state_img[:,:,2]) .÷ 0x3 .+
                 @view(state_img[:,:,3]) .÷ 0x3)
end

"""
Returned array's eltype (TR) != input array's eltype (TI)
and TR <: AbstractFloat and TI <: Integer (pixel vals ∈ [0,255]),
so we scale by 255.0 so output pixels are ∈ [0, 1.0]
"""
function rgb2greymean!(wrap_img::Matrix{TR}, state_img::AbstractArray{TI,3},
    zero_centre::Bool) where {TR <: AbstractFloat, TI <: Integer}
    centre_const = zero_centre ? TR(0.5) : TR(0.0)
    wrap_img .= (TR.(@view(state_img[:,:,1])) .+
                 TR.(@view(state_img[:,:,2])) .+
                 TR.(@view(state_img[:,:,3]))) ./ TR(255.0*3.0)
end

function Reinforce.step!(wrpenv::GreyMeanWrapper, args...)
    r, s = Reinforce.step!(wrpenv.env, args...)
    rgb2greymean!(wrpenv.img, state(wrpenv.env), wrpenv.zero_centre)
    r, state(wrpenv)
end

function Reinforce.state(wrpenv::GreyMeanWrapper)
    wrpenv.img
end

"""
`GreyChanWrapper=>(channel::Int, )`
Turn the rgb to a grey image by selecting just one Channel, and subtracting
128 from all pixels.

state(env) assumed to be encoded as RGB UInt8s in [0,255].
"""
struct GreyChanWrapper{T} <: GreyscaleGymWrapper
    env::AbstractGymEnv
    chan::Int
    zero_centre::Bool
    img::Matrix{T} # allocate on creation, update each `step!`
end

function GreyChanWrapper(env::AbstractGymEnv, chan = 1, zero_centre::Bool = false)
    img = mat_like_state(env)
    rgb2greychan!(img, state(env), chan, zero_centre)
    GreyChanWrapper(env, chan, zero_centre, img)
end

function rgb2greychan!(wrap_img, state_img, chan::Int, zero_centre::Bool)
    zero_factor = zero_centre ? 128.0 : 0.0
    wrap_img .= ((@view state_img[:, :, chan]) .- zero_factor) ./ 255.0
end

function Reinforce.step!(wrpenv::GreyChanWrapper, args...)
    r, s = Reinforce.step!(wrpenv.env, args...)
    rgb2greychan!(wrpenv.img, state(wrpenv.env), wrpenv.chan, wrpenv.zero_centre)
    r, state(wrpenv)
end

function Reinforce.state(wrpenv::GreyChanWrapper)
    wrpenv.img
end

# --------------------------------
# Multi-Frame Wrapper
# --------------------------------
"""
MultiFrameWrapper=>(num_frames::Int)
"""
mutable struct MultiFrameWrapper{T,N} <: GreyscaleGymWrapper
    env::AbstractGymEnv
    nframes::Int
    frame_idx::Int # final index into all-frames, will be incremented on each
        # step, like a sliding window, to avoid copying
    recent_frames::Array{T, N} # allocate on creation, update each step
        # size(recent_frames) == (size(state(wrpenv.env))..., wrpenv.nframes)
end

function init_with_repeated_frame(env, nframes)
    res = cat(ndims(state(env))+1, (state(env) for _ in 1:nframes)...)
    # @dbg ndims(state(env))+1 size(res) nframes size(state(env))
    res
end

"""
`MultiFrameWrapper{T,N}(env, nframes::Int=2)`

N is the rank of the recent_frames tensor, so is one more than the rank of
state(env), i.e. ndims(state(env))+1, e.g.
Atari RGB: state(env) is (m, n, c) so is 3 + 1 = 4
Atari Grey: state(env) is (m, n)   so is 2 + 1 = 3
"""
MultiFrameWrapper(env::AbstractGymEnv, nframes = 2) =
    MultiFrameWrapper(env, nframes, 1, init_with_repeated_frame(env, nframes))

"""
`Reinforce.state{N}(wrpenv::MultiFrameWrapper{T, N})`

Return a view of the last `wrpenv.nframes` frames. Returns `wrpenv.nframes` of
`wrpenv.recent_frames` starting from wrpenv.recent_frames[:, :,( :,)
wrpenv.frame_idx] wrapping around the last dim as needed. The returned array is
of size `(m, n, c * wrpenv.nframes)`, where (m, n, c) are the width, height, and
number of channels of the wrapped env (i.e. of `state(wrpenv.env)`)

assumes N ∈ [3,4], i.e. state(wrpenv.env) is 2D (Greyscale), or 3D (RGB)
"""
function Reinforce.state{T, N}(wrpenv::MultiFrameWrapper{T, N})
    # assert(N ∈ (3,4))
    newsize = ones(Int, 3)
    newsize[1:N-1] .= size(wrpenv.recent_frames)[1:N-1]
    newsize[3] *= wrpenv.nframes # stack all channels from all frames along dim 3
    # @dbg size(wrpenv.recent_frames) newsize wrpenv.nframes N wrpenv.frame_idx
    res = reshape(recent_frames(wrpenv), newsize...)
    # @dbg "state" size(res) typeof(res)
    res
end

"""
`state_img{N}(wrpenv::MultiFrameWrapper{T, N})`
Get an image representation of the state for rendering
returns state(env) with frames stacked side by side
"""
function state_img{T, N}(wrpenv::MultiFrameWrapper{T, N})
    newsize = ones(Int, 3)
    newsize[1:N-1] .= size(wrpenv.recent_frames)[1:N-1]
    newsize[2] *= wrpenv.nframes # stack all images side by side
    newsize[3] == 1 && pop!(newsize) # squeeze for Greyscale
    # @dbg size(wrpenv.recent_frames) newsize wrpenv.nframes N wrpenv.frame_idx
    res = reshape(recent_frames(wrpenv), newsize...)
    # @dbg "state_img" size(res) typeof(res)
    res
end

function recent_frames{T, N}(wrpenv::MultiFrameWrapper{T, N})
    recent_window_indices =
        wraprangeidxs(wrpenv.frame_idx,
                      mod1(wrpenv.frame_idx - 1, wrpenv.nframes),
                      wrpenv.nframes)
    # @dbg N recent_window_indices
    view(wrpenv.recent_frames, colons(N-1)..., recent_window_indices)
end

"""
`state_img(env::AbstractEnvironment)`
Get an image representation of the state
returns state(env) by default
"""
state_img(env::AbstractEnvironment) = state(env)
state_img(wrpenv::AbstractGymWrapper) = state_img(wrpenv.env)

"""
performs:
wrpenv.recent_frames[:, :, ..., :, wrpenv.frame_idx] .= state(wrpenv.env)
and increments `wrpenv.frame_idx`, wrapping it back to 1 if it's `> wrpenv.nframes`
"""
function set_next_frame!{T, N}(wrpenv::MultiFrameWrapper{T, N}, state_img)
    env_state_ndims = N-1
    setindex!(wrpenv.recent_frames, state_img,
              colons(env_state_ndims)..., wrpenv.frame_idx)
    wrpenv.frame_idx = mod1(wrpenv.frame_idx + 1, wrpenv.nframes)
end

function Reinforce.step!(wrpenv::MultiFrameWrapper, args...)
    r, s = Reinforce.step!(wrpenv.env, args...)
    set_next_frame!(wrpenv, state(wrpenv.env))
    r, state(wrpenv)
end

# --------------------------------
# Downsize Wrapper
# --------------------------------
"""
DownsizeWrapper=>(newsize::Tuple{Int,Int})

Return an wrapped env such that all channels of state(env) are scaled to
`newsize`. The input env's state is assumed to be of size `(h, w, c)` where w,h
are the width and height of the image, and c is the number of channels).  The
wrapped state(env), will thus have size `newsize`

Assumes state(env) stays constant size.
"""
struct DownsizeWrapper{T, N} <: GreyscaleGymWrapper
    env::AbstractGymEnv
    scaled_img::Array{T, N} # allocate on creation, update each step
end

"""
`DownsizeWrapper(env, scaler::Tuple{Float64}=(0.5, 0.5))`
"""
function DownsizeWrapper(env, scaler = 0.5)
    newsize = size(state(env)) |> collect
    newsize[1:2] .= round.(newsize[1:2].*scaler)
    DownsizeWrapper(env, Tuple(newsize))
end

"""
`DownsizeWrapper(env, newsize::Tuple{Int,Int})`
"""
function DownsizeWrapper(env, newsize::Tuple{Int,Int})
    DownsizeWrapper(env, imresize(state(env), newsize))
end

"""
`DownsizeWrapper(env, newsize::Int...)`
"""
DownsizeWrapper(env, newsize::Int...) = DownsizeWrapper(env, newsize)

"""
`Reinforce.state(wrpenv::DownsizeWrapper)`

see `DownsizeWrapper(env, scaler)`.
"""
function Reinforce.state(wrpenv::DownsizeWrapper)
    wrpenv.scaled_img
end

function Reinforce.step!(wrpenv::DownsizeWrapper, args...)
    r, s = Reinforce.step!(wrpenv.env, args...)
    imresize!(wrpenv.scaled_img, state(wrpenv.env))
    r, state(wrpenv)
end

downsized(wrapped_env::AbstractGymWrapper) = downsized(wrapped_env.env)
downsized(wrapped_env::DownsizeWrapper) = wrapped_env
downsized(env::GymEnv) = error("env $env doesn't have a DownsizeWrapper")