text:
import Pkg; Pkg.activate(@__DIR__); Pkg.instantiate();
#Write Altro folder instead of @__DIR__
using TrajectoryOptimization
using Plots, LinearAlgebra

#number individual systems
agent_number=10

#defining product system dynamics
function dynamics!(ẋ::AbstractVector{T},x::AbstractVector{T},u::AbstractVector{T}) where T

	for k=1:10
		ẋ[3*k-2] = u[2*k-1]*cos(x[3*k]);
		ẋ[3*k-1] = u[2*k-1]*sin(x[3*k]);
	    ẋ[3*k] = u[2*k];
	end
end

#Total number of states
n=3*agent_number
#Total number of control inputs
m=2*agent_number

#defining descrete model
model = Model(dynamics!,n,m)
model_d = rk4(model)
n = model.n;
m = model.m;

T = Float64;

#inital point for product system
x0 = [0.0;2.0;0.0;
0.0;4.0;0.0;
0.0;6.0;0.0;
0.0;8.0;0.0;
2.0;0.0;pi/2;
4.0;0.0;pi/2;
6.0;0.0;pi/2;
8.0;0.0;pi/2;
0.0;0.0;pi/4;
0.0;10.0;-pi/4]

#goal point for product system
xf = [10.0;2.0;0.0;
10.0;4.0;0.0;
10.0;6.0;0.0;
10.0;8.0;0.0;
2.0;10.0;pi/2;
4.0;10.0;pi/2;
6.0;10.0;pi/2;
8.0;10.0;pi/2;
10.0;10.0;pi/4;
10.0;0.0;-pi/4]

#number of points in trajectory(number of time samples)
N = 104
#sampling time
dt = 0.1
#total time
tf = (N-1)*dt

#LQR Cost Matrices (same as default)
Q = 0.01*Diagonal(I,n)
Qf = 1000.0*Diagonal(I,n)
R = 0.01*Diagonal(I,m)
obj = LQRObjective(Q,R,Qf,xf,N)

#creating upperbound and lower bound for control input with bound constraint
constraints = Constraints(N)
u_min = Array{Float64}(undef,2*agent_number)
u_max = Array{Float64}(undef,2*agent_number)
u_c=2;
for k=1:2*agent_number
	u_min[k]=-u_c;
	u_max[k]=u_c;
end
bnd = BoundConstraint(n,m,u_min=u_min,u_max=u_max)
for k = 1:N-1
     constraints[k] += bnd
 end

 #constraint for obstacle
 #this function is defining penalty for collision constraint
function My_rectangle(x)
	M=10000
	ϵ=1.46
		for i=1:10
        	temp=max((x[3*i-2]-5)^2,(x[3*i-1]-5)^2)
			M=min(temp,M)
    	end
    return -(M-ϵ^2);
end

#constraint for obstacle
#this function is defining penalty for collision constraint
function My_constraint(x)
	ϵ=0.64
	M=1000;
	for i=1:9 #agent number -1
		for j=i+1:10  #agent number
		temp=max((x[3*i-2]-x[3*j-2])^2 , (x[3*i-1]-x[3*j-1])^2);
		M=min(temp,M);
		end
	end
	return -(M - ϵ^2)
end

function My_obs(c,x,u)
    c[1]=My_constraint(x);
    c[2]=My_rectangle(x);
	return nothing
end

obs= Constraint{Inequality}(My_obs,n,m,2,:obs)
for k = 1:N-1
     constraints[k] += obs
 end

#goal constraint
goal = goal_constraint(xf);
constraints[N] += goal

prob = Problem(model_d,obj,constraints=constraints,x0=x0,N=N,dt=dt,xf=xf)

# options for solver
max_con_viol = 1.0e-8
verbose=false

opts_ilqr = iLQRSolverOptions{T}(verbose=verbose,
    iterations=300)

opts_al = AugmentedLagrangianSolverOptions{T}(verbose=verbose,
    opts_uncon=opts_ilqr,
    iterations=40,
    cost_tolerance=1.0e-5,
    cost_tolerance_intermediate=1.0e-4,
    constraint_tolerance=max_con_viol,
    penalty_scaling=10.,
    penalty_initial=1.)

opts_pn = ProjectedNewtonSolverOptions{T}(verbose=verbose,
    feasibility_tolerance=max_con_viol,
    solve_type=:feasible)

opts_altro = ALTROSolverOptions{T}(verbose=verbose,
    opts_al=opts_al,
    R_inf=1.0e-8,
    resolve_feasible_problem=false,
    opts_pn=opts_pn,
    projected_newton=true,
    projected_newton_tolerance=1.0e-4);

	solver = AugmentedLagrangianSolver(prob)

	#here we call the solver to solve
    @time begin
    a=solve!(prob, opts_altro); # solve with ALTRO
    end

#parsing the solution to differenet arrays

	x=Array{Array{Float64,1}}(undef,agent_number)
	z=Array{Array{Float64,1}}(undef,agent_number)
	for i=1:agent_number
    	x[i] = [prob.X[k][3*i-2] for k = 1:N]
    	z[i] = [prob.X[k][3*i-1] for k = 1:N]
	end
#plotting
	a=plot()

	rectangle(w, h, x, y) = Shape(x .+ [0,w,w,0], y .+ [0,0,h,h])
	plot!(rectangle(2,2,4,4), opacity=.5)
	plot!(rectangle(2.8,2.8,3.6,3.6), opacity=.3)

	plot!(x[1],z[1],xlabel="x",ylabel="y",label="1",legend=:topleft,color=:blue,width=2,ratio=:equal,title="Multi-Agent Example Trajectories")
	plot!(x[2],z[2],xlabel="x",ylabel="y",label="2",legend=:topleft,width=2,ratio=:equal)
	plot!(x[3],z[3],xlabel="x",ylabel="y",label="3",legend=:topleft,width=2,ratio=:equal)
	plot!(x[4],z[4],xlabel="x",ylabel="y",label="4",legend=:topleft,width=2,ratio=:equal)
	plot!(x[5],z[5],xlabel="x",ylabel="y",label="5",legend=:topleft,width=2,ratio=:equal)
	plot!(x[6],z[6],xlabel="x",ylabel="y",label="6",legend=:topleft,width=2,ratio=:equal)
	plot!(x[7],z[7],xlabel="x",ylabel="y",label="7",legend=:topleft,width=2,ratio=:equal)
	plot!(x[8],z[8],xlabel="x",ylabel="y",label="8",legend=:topleft,width=2,ratio=:equal)
	plot!(x[9],z[9],xlabel="x",ylabel="y",label="9",legend=:topleft,width=2,ratio=:equal)
	plot!(x[10],z[10],xlabel="x",ylabel="y",label="10",legend=:topleft,width=2,ratio=:equal)

	plot!((x[1][1],z[1][1]),marker=:circle,color=:red,label="")
    plot!((x[1][end],z[1][end]),marker=:utriangle,color=:green,label="")

	plot!((x[2][1],z[2][1]),marker=:circle,color=:red,label="")
    plot!((x[2][end],z[2][end]),marker=:utriangle,color=:green,label="")

	plot!((x[3][1],z[3][1]),marker=:circle,color=:red,label="")
    plot!((x[3][end],z[3][end]),marker=:utriangle,color=:green,label="")

	plot!((x[4][1],z[4][1]),marker=:circle,color=:red,label="")
    plot!((x[4][end],z[4][end]),marker=:utriangle,color=:green,label="")

	plot!((x[5][1],z[5][1]),marker=:circle,color=:red,label="")
    plot!((x[5][end],z[5][end]),marker=:utriangle,color=:green,label="")

	plot!((x[6][1],z[6][1]),marker=:circle,color=:red,label="")
    plot!((x[6][end],z[6][end]),marker=:utriangle,color=:green,label="")

	plot!((x[7][1],z[7][1]),marker=:circle,color=:red,label="")
    plot!((x[7][end],z[7][end]),marker=:utriangle,color=:green,label="")

	plot!((x[8][1],z[8][1]),marker=:circle,color=:red,label="")
    plot!((x[8][end],z[8][end]),marker=:utriangle,color=:green,label="")

	plot!((x[9][1],z[9][1]),marker=:circle,color=:red,label="")
    plot!((x[9][end],z[9][end]),marker=:utriangle,color=:green,label="")

	plot!((x[10][1],z[10][1]),marker=:circle,color=:red,label="")
    plot!((x[10][end],z[10][end]),marker=:utriangle,color=:green,label="")

	io = open("nom_tr.txt", "w")
	for i in prob.U
		for k=1:m-1
    		print(io,i[k])
    		print(io," ")
		end
    	println(io,i[m])
	end
	close(io)

max_stars_repo_name: MehrdadZareian/GMARA

max_issues_repo_name: MehrdadZareian/GMARA

max_forks_repo_name: MehrdadZareian/GMARA


text:
"""
Package: Smoothers

    sma(x, n, center=false)

Smooth a vector of data using a simple moving average

# Arguments
- `x`: Vector of data.
- `n`: Size of the moving average.
- `center`: returns a vector of the same size of x with missing values on the tails.

# Returns
Vector of moving average values

# Examples
```julia-repl
julia> sma(1:5,3)
3-element Vector{Float64}:
 2.0
 3.0
 4.0

julia> sma(1:5,3,true)
5-element Vector{Union{Missing, Float64}}:
  missing
 2.0
 3.0
 4.0
  missing
```
"""
@inline function sma(x::AbstractVector{T}, n::Integer, center::Bool=false) where T<:Real

    n == 1 && return x
    N = length(x)
    @assert 1 <= n <= N

    P = Base.promote_op(/, T, T)
    res = Vector{P}(undef, N-n+1)

    # initial moving average value
    res[1] = ma = sum(x[1:n])/n
    for i in 1:N-n
        @inbounds res[1+i] = ma += (x[n+i] - x[i]) / P(n)
    end

    center ? vcat(repeat([missing], n÷2),res,repeat([missing], n-n÷2-1)) : res

end

ext: jl


text:
module BinomialGPU

using CUDA
using Random

# user-level API
include("rand_binomial.jl")
export rand_binomial!

# CUDA kernels
include("kernels.jl")

end


text:
using Revise
using Plots
using Optim
using ExponentialFamilies
using Combinatorics
using IterTools
using RCall
using Random
using LaTeXStrings

support_grid = -3:0.2:3

EfronTwoTowers = MixtureModel( [Uniform(-1.7, -0.7),
                                Uniform(0.7, 2.7)],
                                [1/3, 2/3])

dem_poly = DiscreteExponentialFamilyModel(support_grid, x->[x x^2 x^3])

tmp_probs = probs(dem_poly, [1.0; -1.0; 0.0])

plot(support_grid,tmp_probs)

tau = -3:0.05:3
@rput tau
R"library(splines)"
R"""Q1 <- scale(splines::ns(tau, 5), center = TRUE, scale = FALSE)
    Q1 <- apply(Q1, 2, function(w) w / sqrt(sum(w * w)))
 """
@rget Q1

dem_spline = DiscreteExponentialFamilyModel(tau, Q1)
d1 = DiscreteExponentialFamily(dem_poly, [1.0; -1.0; 0.0])

d2 = DiscreteExponentialFamily(dem_spline, [1.0; -1.0; 0.0; 1.0; 1.0])

plot(support(d2), probs(d2), seriestype=:sticks)
Zs = rand(d1,1000)

function deconvolved_empirical_moment(Zs, r)
    u_stats = Vector{Float64}(undef, size(Zs,1))
    for (i, row) in enumerate(eachrow(Zs))
        u_stats[i] = mean([prod(comb) for comb in combinations(row,r)])
    end
    mean(u_stats)
end

deconvolved_empirical_moments(Zs, rs) =  deconvolved_empirical_moment.(Ref(Zs), rs)

function empirical_moment(Zs, r::Integer)
    mean(Zs.^r)
end

empirical_moments(Zs, rs) =  empirical_moment.(Ref(Zs), rs)

function deconvolve_moments(dem::DiscreteExponentialFamilyModel, Zs, R; deconv=true)

    if deconv
        m_hat = deconvolved_empirical_moments(Zs,1:R)
    else
        m_hat = empirical_moments(Zs,1:R)
    end

    function my_loss(α)
        moments_α = [moment(dem, α, r) for r in 1:R]
        sum( (moments_α .- m_hat).^2)
    end

    α_init = ones(eltype(dem), size(dem.Q, 2))

    optim_tst = optimize(my_loss, α_init)
    min1 = Optim.minimizer(optim_tst)
    DiscreteExponentialFamily(dem, min1)
end

Random.seed!(1)
n=50_000
μs=rand(MarronWandSkewedBimodal, n)
εs = rand(Normal(0,1),n, 10).*sqrt(10)
Zs = εs .+ μs
Zs_mean = vec(mean(Zs; dims=2))

tmp = deconvolve_moments(dem_spline, Zs, 5)

@rput Zs_mean

R"""
 library(deconvolveR)
 deconv_res <- deconv(tau, Zs_mean, family='Normal')
 a_g_hat <- deconv_res$mle
 deconv_res$mle
 """

@rget a_g_hat
efron_g = DiscreteExponentialFamily(dem_spline, a_g_hat)

pgfplots()
plot_1 = plot(support(tmp), [probs(tmp)./0.05 probs(efron_g)./0.05 pdf.(MarronWandSkewedBimodal,
                                    support(tmp))],
                                    color=["purple" "orange" "black"],
                                    width=[2 2 0.7],
                                    alpha=[1 1 0.7],
                                    linestyle=[:dot :dash :solid],
                                    ylab="Density",
                                    label=["Deconvolved Moments" "Oracle MLE" "Truth"],
                                    legend=:topleft,
                                    xlab=L"\mu")

scale_lapl = rand(Uniform(0,sqrt(5)), n)
εs_lapl = rand(Laplace(0,1),n, 10).*scale_lapl
Zs_lapl = εs_lapl .+ μs
Zs_lapl_mean = vec(mean(Zs; dims=2))

tmp_lapl = deconvolve_moments(dem_spline, Zs_lapl, 5)

@rput Zs_lapl_mean

R"""
 library(deconvolveR)
 deconv_res2 <- deconv(tau, Zs_lapl_mean, family='Normal')
 a_g_hat2 <- deconv_res2$mle
 """

@rget a_g_hat2
efron_g_lapl = DiscreteExponentialFamily(dem_spline, a_g_hat2)

pgfplots()
plot_2 = plot(support(tmp_lapl), [probs(tmp_lapl)./0.05 pdf.(MarronWandSkewedBimodal,
                                    support(tmp_lapl))],
                                    color=["purple"  "black"],
                                    width=[2  0.7],
                                    alpha=[1  0.7],
                                    linestyle=[:dot :solid],
                                    ylab="Density",
                                    label=["Deconvolved Moments" "Truth"],
                                    legend=:topleft,
                                    xlab=L"\mu")

upscale = 0.6#8 #8x upscaling in resolution
default(size=(2000*upscale,400*upscale))
pl = plot(plot_1, plot_2, title=["A" "B"],
            fglegend = :transparent,
            grid = false,
            size=(1500*upscale,600*upscale))
savefig(pl,"deconv_plots2.pdf")
moment_kernel.(cmbs)
moment_kernel(xs) = prod(xs)

fit(dem::DiscreteExponentialFamilyModel, R::Integer, Zs)

R = 3
m_hat = empirical_moments(Zs, 1:3)

my_loss([1.0;-1.0;0.0])

moment(dem_poly, 2, α_init)
moment(d1 )

α_init = ones(eltype(d1), size(d1.Q, 2))

min2 = Optim.minimizer(optim_tst)

moment.(d1, 1:5)

probs(d1.dbn)

moment(dem, [1.0,-1.0,1.0])


text:
foo(x::Integer) = "an integer" #Declare type of function argument
foo(x::ASCIIString) = "a string"
function bar(a, b)
x::Int64 = 0 #Declare type of local variable
y = a+b #Type of variable will be inferred
return y
end

max_stars_repo_path: Module 2/Chapter03/testfunc.jl

max_issues_repo_path: Module 2/Chapter03/testfunc.jl

max_forks_repo_path: Module 2/Chapter03/testfunc.jl


text:
using HDF5
using LaTeXStrings
using Plasmons
using Plots

pgfplotsx()

function data_for_figure_1(n::Int = 32)
    x, y, z = h5open("input_square_sheet_$(n)x$(n).h5", "r") do io
        read(io["x"]), read(io["y"]), read(io["z"])
    end
    χ, ε, ωs, qs = h5open("result_square_sheet_$(n)x$(n).h5", "r") do io
        ωs = map(d->real(read(attributes(d)["ħω"])), io["/ε"])
        number_qs = n + 1
        qs = collect(0:(π / (number_qs - 1)):π)
        direction = (1.0, 0.0, 0.0)
        # Each row corresponds to an ħω, each column -- to a q
        χ = dispersion((read(d) for d in io["/χ"]), map(q-> q .* direction, qs), x, y, z)
        ε = dispersion((read(d) for d in io["/ε"]), map(q-> q .* direction, qs), x, y, z)
        χ, ε, ωs, qs
    end
    h5open("dispersion_square_sheet_$(n)x$(n).h5", "w") do io
        io["ε"] = ε
        io["χ"] = χ
        io["ħω"] = ωs
        io["q"] = qs
    end
    nothing
end

function figure_1()
    theme(
        :solarized_light;
        # fg = RGB(88 / 255, 110 / 255, 117 / 255),
        fg = RGB(7 / 255, 54 / 255, 66 / 255),
        fgtext = RGB(7 / 255, 54 / 255, 66 / 255),
        fgguide = RGB(7 / 255, 54 / 255, 66 / 255),
        fglegend = RGB(7 / 255, 54 / 255, 66 / 255),
    )
    χ, ε, ωs, qs = h5open("dispersion_square_sheet_32x32.h5", "r") do io
        read(io["χ"]), read(io["ε"]), read(io["ħω"]), read(io["q"])
    end
    # Compute loss function
    g₁ = heatmap(
        qs,
        ωs,
        transpose(@. -imag(1 / ε)),
        ylims = (0, 15),
        xlabel = raw"$q$, $1/a$",
        ylabel = raw"$\hbar\omega$, eV",
        title = L"$-\mathrm{Im}\left[1 / \varepsilon(q, \omega) \right]$",
    )
    g₂ = heatmap(
        qs,
        ωs,
        transpose(@. -imag(χ)),
        ylims = (0, 15),
        xlabel = nothing,
        ylabel = nothing,
        title = L"$-\mathrm{Im}\left[\chi(q, \omega)\right]$",
    )
    savefig(plot(g₁, g₂, size = (700, 350)), "dispersion_square_sheet_32x32.png")
    nothing
end


text:
export gammatrick!, gammatrick

#
# EVALUATION + DIFFERENTATION
#
function evaluate(H::AbstractPolynomialHomotopy{T}, x::Vector{S}, t::Number) where {T, S}
    evaluate!(zeros(H.target, promote_type(T, S)), H, x, t)
end

function evaluate(H::AbstractPolynomialHomotopy{T}, x::Vector{S}, t::Number, cfg::PolynomialHomotopyConfig, precomputed=false) where {T, S}
    evaluate!(zeros(H.target, promote_type(T, S)), H, x, t, cfg, precomputed)
end

function jacobian(H::AbstractPolynomialHomotopy{T}, x::AbstractVector, t, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    u = similar(jacobian_target(cfg))
    jacobian!(u, H, x, t, cfg, precomputed)
    u
end

function dt(H::AbstractPolynomialHomotopy{T}, x::AbstractVector, t, cfg::PolynomialHomotopyConfig, precomputed=false) where {T<:Number}
    u = similar(value_start(cfg))
    dt!(u, H, x, t, cfg, precomputed)
    u
end

#
# EQUALITY
#
function ==(H1::T, H2::T) where {T<:AbstractPolynomialHomotopy}
    H1.start == H2.start && H1.target == H2.target
end
function Base.isequal(H1::T, H2::T) where {T<:AbstractPolynomialHomotopy}
    Base.isequal(H1.start, H2.start) && Base.isequal(H1.target, H2.target)
end

function Base.show(io::IO, H::AbstractPolynomialHomotopy)
    println(io, typeof(H), " with $(length(H.start)) polynomials in $(nvariables(H)) variables.")
end

"""
    gammatrick!(H::AbstractPolynomialHomotopy{Complex} [, seed::Int]])

Scale the coefficients of the start system of `H` with a random complex number picked
uniformly from the (complex) unit circle. Use this to make the paths ``z(t)`` generic.

    gammatrick!(H::AbstractPolynomialHomotopy{Complex}, γ::Complex)

You can also pass a scaling factor directly.
"""
function gammatrick!(H::AbstractPolynomialHomotopy{T}, γ::Union{AbstractFloat, Complex}) where {T<:Complex}
    FP.scale_coefficients!.(H.start, convert(T, γ))
    γ
end
gammatrick!(H::AbstractPolynomialHomotopy{T}) where T = gammatrick!(H, randomgamma(T))
gammatrick!(H::AbstractPolynomialHomotopy{T}, seed::Int) where T = gammatrick!(H, randomgamma(T, seed))

"""
    gammatrick(H::AbstractPolynomialHomotopy{Complex} , γ::Number)

Scale the coefficients of the start system of `H` with `γ`.

    gammatrick(H::AbstractPolynomialHomotopy{Complex})

A a random complex number `γ` is picked uniformly from the (complex) unit circle and then
scale the coefficients of the start system of `H` with `γ`.
This returns the new `H` and `γ`.
"""
function gammatrick(H::AbstractPolynomialHomotopy)
    K = deepcopy(H)
    γ = gammatrick!(K)
    K, γ
end
function gammatrick(H::AbstractPolynomialHomotopy, γ::Number)
    K = deepcopy(H)
    gammatrick!(K, γ)
    K
end

randomgamma(::Type{T}) where {T<:Complex} = convert(T, exp(im * (rand() * 2π - π)))
function randomgamma(::Type{T}, seed::Int) where {T<:Complex}
    srand(seed)
    θ = rand() * 2π - π
    convert(T, exp(im * θ))
end

function homogenize(F::Vector{<:FP.Polynomial}, G::Vector{<:FP.Polynomial})
    F_ishomogenous = all(FP.ishomogenous, F)
    G_ishomogenous = all(FP.ishomogenous, G)

    if F_ishomogenous && G_ishomogenous
        F1 = map(f -> FP.homogenize(f, respect_homogenous=true), F)
        G1 = map(g -> FP.homogenize(g, respect_homogenous=true), G)
        (F1, G1)
    else
        F1 = map(f -> FP.homogenize(f, respect_homogenous=false), F)
        G1 = map(g -> FP.homogenize(g, respect_homogenous=false), G)
        (F1, G1)
    end
end

function homogenize(H::AbstractPolynomialHomotopy)
    G, F = homogenize(H.start, H.target)
    typeof(H)(G, F)
end
function dehomogenize(H::AbstractPolynomialHomotopy)
    typeof(H)(FP.dehomogenize.(H.start), FP.dehomogenize.(H.target))
end

function ishomogenized(H::AbstractPolynomialHomotopy)
    all(FP.ishomogenized, H.start) && all(FP.ishomogenized, H.target)
end
function ishomogenous(H::AbstractPolynomialHomotopy)
    all(FP.ishomogenous, H.start) && all(FP.ishomogenous, H.target)
end

nvariables(H::AbstractPolynomialHomotopy) = FP.nvariables(H.start[1])
Base.length(H::AbstractPolynomialHomotopy) = length(H.start)


text:
@testset "Unit tests" begin
    include("unit/auto-generated-source-files.jl")
    include("unit/fhir-to-julia.jl")
    include("unit/other-fhir-versions.jl")
    include("unit/requests.jl")
end


text:
"""
    SegmentedMean{V <: AbstractVector{<:Number}}} <: AbstractAggregation

[`AbstractAggregation`](@ref) implementing segmented mean aggregation:

``
f(\\{x_1, \\ldots, x_k\\}) = \\frac{1}{k} \\sum_{i = 1}^{k} x_i
``

Stores a vector of parameters `ψ` that are filled into the resulting matrix in case an empty bag is encountered.

See also: [`AbstractAggregation`](@ref), [`AggregationStack`](@ref), [`mean_aggregation`](@ref),
    [`SegmentedMax`](@ref), [`SegmentedSum`](@ref), [`SegmentedPNorm`](@ref), [`SegmentedLSE`](@ref).
"""
struct SegmentedMean{V <: AbstractVector{<:Number}} <: AbstractAggregation
    ψ::V
end

Flux.@functor SegmentedMean

SegmentedMean(T::Type, d::Int) = SegmentedMean(zeros(T, d))
SegmentedMean(d::Int) = SegmentedMean(Float32, d)

Flux.@forward SegmentedMean.ψ Base.getindex, Base.length, Base.size, Base.firstindex, Base.lastindex,
        Base.first, Base.last, Base.iterate, Base.eltype

Base.vcat(as::SegmentedMean...) = reduce(vcat, as |> collect)
function Base.reduce(::typeof(vcat), as::Vector{<:SegmentedMean})
    SegmentedMean(reduce(vcat, [a.ψ for a in as]))
end

function (a::SegmentedMean)(x::Maybe{AbstractMatrix{T}}, bags::AbstractBags,
                               w::Optional{AbstractVecOrMat{T}}=nothing) where T
    _check_agg(a, x)
    segmented_mean_forw(x, a.ψ, bags, w)
end

segmented_mean_forw(::Missing, ψ::AbstractVector, bags::AbstractBags, w) = repeat(ψ, 1, length(bags))
function segmented_mean_forw(x::AbstractMatrix, ψ::AbstractVector, bags::AbstractBags, w::Optional{AbstractVecOrMat})
    t = promote_type(eltype(x), eltype(ψ))
    y = zeros(t, size(x, 1), length(bags))
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                y[i, bi] = ψ[i]
            end
        else
            for j in b
                for i in 1:size(x, 1)
                    y[i, bi] += _weight(w, i, j, t) * x[i, j]
                end
            end
            @views y[:, bi] ./= _bagnorm(w, b)
        end
    end
    y
end

function segmented_mean_back(Δ, y, x, ψ, bags, w)
    dx = zero(x)
    dψ = zero(ψ)
    dw = isnothing(w) ? ZeroTangent() : zero(w)
    @inbounds for (bi, b) in enumerate(bags)
        if isempty(b)
            for i in eachindex(ψ)
                dψ[i] += Δ[i, bi]
            end
        else
            ws = _bagnorm(w, b)
            for j in b
                for i in 1:size(x, 1)
                    dx[i, j] += _weight(w, i, j, eltype(x)) * Δ[i, bi] / _weightsum(ws, i)
                    ∇dw_segmented_mean!(dw, Δ, x, y, w, ws, i, j, bi)
                end
            end
        end
    end
    dx, dψ, NoTangent(), dw
end

function segmented_mean_back(Δ, y, x::Missing, ψ, bags, w)
    dψ = zero(ψ)
    @inbounds for (bi, b) in enumerate(bags)
        for i in eachindex(ψ)
            dψ[i] += Δ[i, bi]
        end
    end
    ZeroTangent(), dψ, NoTangent(), ZeroTangent()
end

∇dw_segmented_mean!(dw::ZeroTangent, Δ, x, y, w::Nothing, ws, i, j, bi) = nothing
function ∇dw_segmented_mean!(dw::AbstractVector, Δ, x, y, w::AbstractVector, ws, i, j, bi)
    dw[j] += Δ[i, bi] * (x[i, j] - y[i, bi]) / ws
end
function ∇dw_segmented_mean!(dw::AbstractMatrix, Δ, x, y, w::AbstractMatrix, ws, i, j, bi)
    dw[i, j] += Δ[i, bi] * (x[i, j] - y[i, bi]) / ws[i]
end

function ChainRulesCore.rrule(::typeof(segmented_mean_forw), args...)
    y = segmented_mean_forw(args...)
    grad = Δ -> (NoTangent(), segmented_mean_back(Δ, y, args...)...)
    y, grad
end


text:
using Distributions, MCMC

# Using Optim, M0 and M1 would be ~ 5.01, acceptance rate ~ 60%
#using Optim

using Base.Test
srand(1)

L = 1

# Target distribution shape, unscaled
g(x) = exp(-dot(x, x)/2)
log_g(x) = -dot(x, x)/2

# Candidate (starting) density, use Normal(0, 2)
g0(x) = 1/(sqrt(2*pi)*2)*exp(-dot(x, x)/8)
log_g0(x) = -dot(x, x)/8-log(2 * sqrt(2*pi))

# Compute M0 such that candidate dominates target
h0(x) = -g(x)/g0(x)

#res0 = optimize(h0, ones(L))
#M0 = -res0.f_minimum

M0 = 10.0
log_M0 = log(M0)

# Candidate (starting) Distribution, use Normal(0, 2)
g1(x) = exp(sum(logpdf(Normal(0, 2), x)))
log_g1(x) = sum(logpdf(Normal(0, 2), x))

# Compute M1 such that candidate dominates target
h1(x) = -g(x)/sum(g1(x))
#res1 = optimize(h1, ones(L))
#M1 = -res1.f_minimum
M1 = 10.0
log_M1 = log(M1)

mcmodel = model(log_g, init=ones(L))

mcchain = run(mcmodel, [ARS(log_g0, log_M0), ARS(log_g1, log_M1)], SerialMC(1000:1:10000))

println()
println("M: ", [M0 M1])
println()

println("Acceptance rate: ", acceptance(mcchain[1]))
describe(mcchain[1])

println("Acceptance rate: ", acceptance(mcchain[2]))
describe(mcchain[2])

theta = hcat(mcchain[1].samples, mcchain[2].samples);
for i in sample(1:length(theta), 200)
  @test_approx_eq h0(theta[i])/M0 h1(theta[i])/M1
end


text:
##### temperature function
@inline function tempFunc(temp, p)
    TFunc = exp(-p.Ea/(8.3145*(temp+273.15)))*(1.0-exp(temp+273.15 - p.T⁺))
    return max(0.0, TFunc)
end

##### calculate photosynthesis rate (mmolC/individual/second)
@inline function calc_PS(par, temp, Chl, Bm, Sz, p)
    αI  = par * p.α * p.Φ
    PCm = p.PCmax * Sz^p.PC_b * tempFunc(temp, p)
    PS  = PCm * (1.0 - exp(-αI / max(1.0e-10, PCm) * Chl / max(1.0e-10, Bm))) * Bm
    return PS
end

@kernel function calc_inorganic_uptake_kernel!(plank, nuts, p)
    i = @index(Global)
    @inbounds plank.PS[i] = calc_PS(nuts.par[i], nuts.T[i], plank.Chl[i], plank.Bm[i], plank.Sz[i], p) * plank.ac[i]
end
function calc_inorganic_uptake!(plank, nuts, p, arch::Architecture)
    kernel! = calc_inorganic_uptake_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, nuts, p)
    wait(device(arch), event)
    return nothing
end

##### calculate respiration (mmolC/individual/second)
@kernel function calc_respir_kernel!(plank, T, p)
    i = @index(Global)
    @inbounds plank.resp[i] = p.respir_a * plank.Sz[i]^p.respir_b * plank.Bm[i] * tempFunc(T[i], p) * plank.ac[i]
end
function calc_respir!(plank, T, p, arch)
    kernel! = calc_respir_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, T, p)
    wait(device(arch), event)
    return nothing
end

##### update C quotas
@kernel function update_quotas_kernel!(plank, ΔT)
    i = @index(Global)
    @inbounds plank.Bm[i]  += (plank.PS[i] - plank.resp[i]) * ΔT
    @inbounds plank.age[i] += ΔT / 3600.0 * plank.ac[i]
end
function update_quotas!(plank, ΔT, arch)
    kernel! = update_quotas_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, ΔT)
    wait(device(arch), event)
    return nothing
end

##### update cell size
@kernel function update_cellsize_kernel!(plank, p)
    i = @index(Global)
    @inbounds plank.Sz[i]  = plank.Bm[i] / (p.Cquota * p.Nsuper)
    @inbounds plank.Chl[i] = plank.Bm[i] * p.Chl2C
end
function update_cellsize!(plank, p, arch)
    kernel! = update_cellsize_kernel!(device(arch), 256, (size(plank.ac,1)))
    event = kernel!(plank, p)
    wait(device(arch), event)
    return nothing
end


text:
"""file with functions used in wind farm optimization, including constraints
created May 4, 2020
author: PJ Stanley and Jared Thomas
contributors: Nicholas F. Baker and Wesley Holt
"""

"""
    turbine_spacing(turbine_x,turbine_y)

Calculate the distance between turbines in a wind farm. There is an infinite gradient
of this function if two points are exactly the same. This can be avoided by returning the
square of the turbine spacing rather than the actual distance, but it makes the gradients
scale much more poorly. Because it is very vanishinly rare to have turbines exactly in the
same location, this function leaves the square root in the calculations.

# Arguments
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
"""
function turbine_spacing(turbine_x, turbine_y)
    nturbines = length(turbine_x)
    spacing_vec =
        zeros(typeof(turbine_x[1]), Int((nturbines) * (nturbines - 1) / 2))
    k = 1
    for i in 1:nturbines
        for j in i+1:nturbines
            spacing_vec[k] = sqrt(
                (turbine_x[j] - turbine_x[i])^2 +
                (turbine_y[j] - turbine_y[i])^2,
            )
            k += 1
        end
    end
    return spacing_vec
end

"""
    circle_boundary(center,radius,turbine_x,turbine_y)

calculate the distance from each turbine to a circular boundary. Negative means the
turbine is inside the boundary

# Arguments
- `center::Float`: circular boundary center [x,y]
- `radius::Float`: circulat boundary radius
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
"""
function circle_boundary(center, radius, turbine_x, turbine_y)
    nturbines = length(turbine_x)
    boundary_vec = zeros(typeof(turbine_x[1]), nturbines)
    for i in 1:nturbines
        boundary_vec[i] =
            (center[1] - turbine_x[i])^2 + (center[2] - turbine_y[i])^2 -
            radius^2
    end
    return boundary_vec
end

"""
    convex_boundary(boundary_vertices,boundary_normals,turbine_x,turbine_y)

calculate the distance from each turbine to a possibly non-circular, but convex boundary. Negative means the
turbine is inside the boundary

# Arguments
- `boundary_vertices::Array{Float,2}`: vertices of the convex hull CCW in order s.t.
        boundaryVertices[i] -> first point of face for unit_normals[i]
- `boundary_normals::Array{Float,2}`: unit normal vector for each boundary face
        CCW where boundaryVertices[i] is the first point of the corresponding face
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
"""
function convex_boundary(boundary_vertices, boundary_normals, turbine_x, turbine_y)
    nturbines = length(turbine_x)
        nVertices = size(boundary_vertices)[1]
        # initialize array to hold distances from each point to each face
        face_distance = zeros(typeof(turbine_x[1]),(nturbines, nVertices))

        # loop through points and find distance to each face
        for i = 1:nturbines
                # determine if point is inside or outside of each face, and distance from each face
                for j = 1:nVertices
                        # define the vector from the point of interest to the first point of the face
                        pa = [boundary_vertices[j, 1]-turbine_x[i], boundary_vertices[j, 2]-turbine_y[i]]
                        # find perpendicular distance from point to current surface (vector projection)
                        d_vec = sum(pa .* boundary_normals[j,:]) .* boundary_normals[j,:]
                        # calculate the sign of perpendicular distance from point to current face (- is inside, + is outside)
                        face_distance[i, j] = -sum(d_vec .* boundary_normals[j,:])
                end
        end

        return vcat(face_distance...)
end

"""
    splined_boundary(turbine_x, turbine_y, bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies)

calculate the distance from each turbine to a closed boundary made up of zero or
more reflex angles (concavities). Boundary will have three or four user-selected
"corners", such that the "sides" between corners (that will be splined) are
injective functions (meaning that for every x-coord, there exists only one
corresponding y-coord). Returns four values for every turbine, corresponding to
the distance from the turb to the upper, lower, left, and right splined "sides".
A negative return value means the turb is inside the boundary for that "side".
Returns a single array of {Float64} of length {length(turbine_x) * 4}. Note that
all boundary coordinates must be in the first quadrant of the Cartesian
coordinate system (+x and +y values only)

# Arguments
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
- `bndry_x_clsd::Array{Float}`: x locations of boundary vertices, CCW in
        order s.t. boundaryVertices[0] is the NE corner of boundary, and with
        the first vertex duplicated at the end for completeness of calcs.
- `bndry_y_clsd::Array{Float}`: y locations of boundary vertices, CCW in
        order s.t. boundaryVertices[0] is the NE corner of boundary, and with
        the first vertex duplicated at the end for completeness of calcs.
- `bndry_corner_indcies::Array{Float}`: An array of 3 or 4 indicies in the
        bndry_x/y_clsd arrays that correspond to the three four "corners" used
        between splined "sides"
"""
function splined_boundary(turbine_x, turbine_y, bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies)
    """ Returns if the passed turbines are within the passed closed boundary """
    num_turbs = Int8(length(turbine_x))
    num_sides = Int8(length(bndry_corner_indcies)-1)
    x_min_indx = 3                      # Default to work w/ squared boundaries
    if num_sides == 3                   # If we only have 3 corners
        x_min_indx = 2                  # denote that it's a triangle boundary
    end

    # Check to make sure our points are in
    num_cons = Int(num_turbs * 4)
    bndry_cons = zeros(typeof(turbine_x[1]), num_cons)   # 4 values (2 x and 2 y) for each turb

    x_max = bndry_x_clsd[bndry_corner_indcies[1]]           # Our maximum x-value
    x_min = bndry_x_clsd[bndry_corner_indcies[x_min_indx]]  # Our min x-value

    # For every turbine
    for cntr in 1:num_turbs
        place = (cntr-1)*4
        #- Calc x-vals
        bndry_cons[place+1] = (x_max - turbine_x[cntr])   # Positive good, neg bad
        bndry_cons[place+2] = (turbine_x[cntr] - x_min)   # pos good, neg bad

        #- Calc y-vals
        y_max,y_min = getUpDwnYvals(turbine_x[cntr], bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies)
        bndry_cons[place+3] = (y_max - turbine_y[cntr])
        bndry_cons[place+4] = (turbine_y[cntr] - y_min)
    end

    return -bndry_cons # Invert values so Negative is inside boundary
end

"""
    splined_boundary_discreet_regions(turbine_x, turbine_y, bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies, turbs_per_region)

Uses FLOWFarm's splined_boundary() function to calculate the turbine-boundary
constraints for one or more discreet regions, with pre-allocated turbines
for each region. Returns four values for every turbine, corresponding to the
distance from each turb to the upper, lower, left, and right splined "sides" for
the region to which it was allocated. A negative return value means the turb is
outside the "side" of boundary for which it has been allocated. Returns a single
array of {Float64} of length {length(turbine_x) * 4}. Note that all boundary
coordinates must be in the first quadrant of the Cartesian coordinate system
(+x and +y values only)

# Arguments
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
- `bndry_x_clsd::Array{Float}`: x locations of boundary vertices, CCW in
        order s.t. boundaryVertices[0] is the NE corner of boundary, and with
        the first vertex duplicated at the end for completeness of calcs.
- `bndry_y_clsd::Array{Float}`: y locations of boundary vertices, CCW in
        order s.t. boundaryVertices[0] is the NE corner of boundary, and with
        the first vertex duplicated at the end for completeness of calcs.
- `bndry_corner_indcies::Array{Float}`: An array of 3 or 4 indicies in the
        bndry_x/y_clsd arrays that correspond to the three four "corners" used
        between splined "sides"
- 'turbs_per_region::Array{Int}`: An array of length equivalent to the number of
        discrete boundary regions, with each element denoting howmany turbines
        are apportioned to the corresponding region. sum(turbs_per_region) must
        be equivalent to the total number of turbines in the windfarm
"""
function splined_boundary_discreet_regions(turbine_x, turbine_y, bndry_x_clsd, bndry_y_clsd, num_bndry_verts, bndry_corner_indcies, turbs_per_region)
    """ Goes through numerous discrete splined-boundary regions and returns if the apportioned turbines are within their region """
    num_regions = length(turbs_per_region)
    bndry_constraints = zeros(typeof(turbine_x[1]), sum(turbs_per_region)*4)#[ Float64[] for i in 1:num_regions ]  # To hold cnstrnts
    #-- Loop through and do all regions --#
    prev_turb_index = 1
    bndry_vert_index = 1
    for cntr in 1:num_regions
        next_turb_index = ((turbs_per_region[cntr]-1) + prev_turb_index)  # Next index for our Turbines
        region_turbine_x = turbine_x[prev_turb_index:next_turb_index]   # Simplified list of turbines preallocated to this region
        region_turbine_y = turbine_y[prev_turb_index:next_turb_index]
        cnstrnts_index_strt = ((prev_turb_index-1)*4)+1
        cnstrnts_index_end = ((next_turb_index-1)*4)+4
        bndry_constraints[cnstrnts_index_strt:cnstrnts_index_end] = splined_boundary(region_turbine_x, region_turbine_y, bndry_x_clsd[cntr], bndry_y_clsd[cntr], bndry_corner_indcies[bndry_vert_index:(bndry_vert_index+(num_bndry_verts[cntr]-1))])
        bndry_vert_index += num_bndry_verts[cntr]
        prev_turb_index += turbs_per_region[cntr]
    end

    # Make it a long 1D array for SNOPT
    return bndry_constraints
end

"""
    getUpDwnYvals(turbine_x, bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies)

Supplements FLOWFarm's splined_boundary() function by calculating (for a given x
location) the maximum and minimum y-value permitted to remain "inside" the
boundary. If turbine_x is located left of the boundary's leftmost vertex or
right of the boundary's rightmost vertex, it return's that corresponding
vertex's y-value as the max and min, as default. Returns two values, the minimum
and maximum interior y-values withing a boundary for the given turbine_x value.
Note that all boundary coordinates must be in the first quadrant of the
Cartesian coordinate system (+x and +y values only)

# Arguments
- `turbine_x::Array{Float}`: x-value of the turbine being examined
- `bndry_x_clsd::Array{Float}`: x locations of boundary vertices, CCW in
        order s.t. boundaryVertices[0] is the NE corner of boundary, and with
        the first vertex duplicated at the end for completeness of calcs.
- `bndry_y_clsd::Array{Float}`: y locations of boundary vertices, CCW in
        order s.t. boundaryVertices[0] is the NE corner of boundary, and with
        the first vertex duplicated at the end for completeness of calcs.
- `bndry_corner_indcies::Array{Float}`: An array of 3 or 4 indicies in the
        bndry_x/y_clsd arrays that correspond to the three four "corners" used
        between splined "sides"
"""
function getUpDwnYvals(turbine_x, bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies)
    # Given that there are 4 "sides" (with 3&4 below, 1&2 above),
    # returns the splined y-vals the given x-coord falls between
    bool_triangle = false                # Default to work w/ squared boundaries
    if length(bndry_corner_indcies) == 4 # If we only have 3 corners (4 closed)
        bool_triangle = true             # denote that it's a triangle boundary
    end

    #-- Upper edge calculation
    if (bool_triangle) # If we're dealing with a 3-sided boundary
        # If it's out of bounds to the left of the upper boundary
        if (turbine_x <= bndry_x_clsd[bndry_corner_indcies[2]])
            # Give it the y-value of our leftmost point
            y_max = bndry_y_clsd[bndry_corner_indcies[2]]
        # If it's right of the leftmost corner, but left of the rightmost corner
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[1]])
            y_max = linear(bndry_x_clsd[bndry_corner_indcies[2]:-1:bndry_corner_indcies[1]],
                            bndry_y_clsd[bndry_corner_indcies[2]:-1:bndry_corner_indcies[1]], turbine_x)  # Make it the left upper spline #############
        else # Otherwise it's right of our top right corner
            # Give it the y-value of our rightmost corner
            y_max = bndry_y_clsd[bndry_corner_indcies[1]]
        end
    else # If we're dealing with a 4-sided boundary
        # If left of the bottom left corner
        if (turbine_x <= bndry_x_clsd[bndry_corner_indcies[3]])
            # Give it the y-value of our bottom left corner
            y_max = bndry_y_clsd[bndry_corner_indcies[3]]
        # If it's to the left of the right upper spline
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[2]])
            y_max = linear(bndry_x_clsd[bndry_corner_indcies[3]:-1:bndry_corner_indcies[2]],
                            bndry_y_clsd[bndry_corner_indcies[3]:-1:bndry_corner_indcies[2]], turbine_x)
        # If it's to the left of the rightmost point
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[1]])
            y_max = linear(bndry_x_clsd[bndry_corner_indcies[2]:-1:bndry_corner_indcies[1]],
                            bndry_y_clsd[bndry_corner_indcies[2]:-1:bndry_corner_indcies[1]], turbine_x)
        else # Otherwise, if it's out of bounds to the right
            # Give it the y-value of our rightmost point
            y_max = bndry_y_clsd[bndry_corner_indcies[1]]
        end
    end
    #-- Lower edge calculation.
    if (bool_triangle)
        # If it's to the left of the leftmost lower spline
        if (turbine_x <= bndry_x_clsd[bndry_corner_indcies[2]])
            y_min = bndry_y_clsd[bndry_corner_indcies[2]] # leftmost y-value is corner #2
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[3]])
            # Spline from bottom left to bottom middle
            y_min = linear(bndry_x_clsd[bndry_corner_indcies[2]:bndry_corner_indcies[3]],
                            bndry_y_clsd[bndry_corner_indcies[2]:bndry_corner_indcies[3]], turbine_x)
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[1]])
            # Spline from bottom middle to bottom right
            y_min = linear(bndry_x_clsd[bndry_corner_indcies[3]:bndry_corner_indcies[4]],
                            bndry_y_clsd[bndry_corner_indcies[3]:bndry_corner_indcies[4]], turbine_x)
        else
            # Give it the y-value of our rightmost point
            y_min = bndry_y_clsd[bndry_corner_indcies[1]]
        end
    else # if it's a square
        # If it's to the left of the leftmost lower spline
        if (turbine_x <= bndry_x_clsd[bndry_corner_indcies[3]])
            y_min = bndry_y_clsd[bndry_corner_indcies[3]] # leftmost y-value is corner #3
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[4]])
            # Spline from bottom left to bottom middle
            y_min = linear(bndry_x_clsd[bndry_corner_indcies[3]:bndry_corner_indcies[4]],
                            bndry_y_clsd[bndry_corner_indcies[3]:bndry_corner_indcies[4]], turbine_x)
        elseif (turbine_x < bndry_x_clsd[bndry_corner_indcies[5]])
            # Spline from bottom middle to bottom right
            y_min = linear(bndry_x_clsd[bndry_corner_indcies[4]:bndry_corner_indcies[5]],
                            bndry_y_clsd[bndry_corner_indcies[4]:bndry_corner_indcies[5]], turbine_x)
        else
            # Give it the y-value of our rightmost point
            y_min = bndry_y_clsd[bndry_corner_indcies[1]]
        end
    end

    return y_max, y_min
end

"""
    ray_trace_boundary(boundary_vertices,boundary_normals,turbine_x,turbine_y;discrete=false,smooth_max_factor=20.0)

Calculate the distance from each turbine to the nearest point on the boundary.
Negative means the turbine is inside the boundary.

# Arguments
- `boundary_vertices::Array{Float,2}`: vertices of the boundary CCW in order s.t.
        boundaryVertices[i] -> first point of face for unit_normals[i]
- `boundary_normals::Array{Float,2}`: unit normal vector for each boundary face
        CCW where boundaryVertices[i] is the first point of the corresponding face
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
- `discrete::Bool`: true if there are multiple discrete regions, false if there is only one region
- `smooth_max_factor::Float`: factor to control smoothness for smoothmax. The lower the number, the smoother the constraint function
"""
function ray_trace_boundary(boundary_vertices,boundary_normals,turbine_x,turbine_y;discrete=false,smooth_max_factor=100.0)

    # single region
    if discrete == false

        # number of turbines and boundary vertices
        nturbines = length(turbine_x)
        nvertices = size(boundary_vertices)[1]

        # initialize constraint output values
        c = zeros(typeof(turbine_x[1]),(nturbines))

        # initialize array to hold distances from each turbine to closest boundary face
        turbine_to_face_distance = zeros(typeof(turbine_x[1]),(nvertices))

        # add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
        boundary_vertices = [boundary_vertices; boundary_vertices[1,1] boundary_vertices[1,2]]

        # iterate through each turbine location
        for i = 1:nturbines

            # initialize intersection counter
            intersection_counter = 0

            # get vector from turbine to the first vertex in first face
            turbine_to_first_facepoint = boundary_vertices[1, :] - [turbine_x[i]; turbine_y[i]]

            # iterate through each boundary
            for j = 1:nvertices

                # check if y-coordinate of turbine is less than at least one y-coordinate of the two boundary vertices
                if !(boundary_vertices[j, 2] < turbine_y[i] && boundary_vertices[j+1, 2] < turbine_y[i])        # (this might not be necessary)

                    # check if x-coordinate of turbine is between the x-coordinates of the two boundary vertices
                    if boundary_vertices[j, 1] < turbine_x[i] < boundary_vertices[j+1, 1] || boundary_vertices[j, 1] > turbine_x[i] > boundary_vertices[j+1, 1]

                        # check to see if the turbine is below the boundary
                        if turbine_y[i] < (boundary_vertices[j+1, 2] - boundary_vertices[j, 2]) / (boundary_vertices[j+1, 1] - boundary_vertices[j, 1]) * (turbine_x[i] - boundary_vertices[j, 1]) + boundary_vertices[j, 2]

                            # the vertical ray intersects the boundary
                            intersection_counter += 1

                        end

                    end

                end

                # define the vector from the turbine to the second point of the face
                turbine_to_second_facepoint = boundary_vertices[j+1, :] - [turbine_x[i]; turbine_y[i]]

                # find perpendicular distance from turbine to current face (vector projection)
                boundary_vector = boundary_vertices[j+1, :] - boundary_vertices[j, :]

                # check if perpendicular distance is the shortest
                if sum(boundary_vector .* -turbine_to_first_facepoint) > 0 && sum(boundary_vector .* turbine_to_second_facepoint) > 0

                    # perpendicular distance from turbine to face
                    turbine_to_face_distance[j] = abs(sum(turbine_to_first_facepoint .* boundary_normals[j,:]))

                # check if distance to first facepoint is shortest
                elseif sum(boundary_vector .* -turbine_to_first_facepoint) <= 0

                    # distance from turbine to first facepoint
                    turbine_to_face_distance[j] = sqrt(sum(turbine_to_first_facepoint.^2))

                # distance to second facepoint is shortest
                else

                    # distance from turbine to second facepoint
                    turbine_to_face_distance[j] = sqrt(sum(turbine_to_second_facepoint.^2))

                end

                # reset for next face iteration
                turbine_to_first_facepoint = turbine_to_second_facepoint        # (for efficiency, so we don't have to recalculate for the same vertex twice)

            end

            # magnitude of the constraint value
            c[i] = -ff.smooth_max(-turbine_to_face_distance, s=smooth_max_factor)

            # sign of the constraint value (- is inside, + is outside)
            if mod(intersection_counter, 2) == 1
                c[i] = -c[i]
            end

        end

        return c

    # multiple discrete regions
    else

        # number of turbines
        nturbines = length(turbine_x)

        # number of regions
        nregions = length(boundary_vertices)

        # initialize constraint output values
        c = zeros(typeof(turbine_x[1]),(nturbines))

        # initialize turbine status vector
        status = zeros(Int64, nturbines)

        # iterate through each turbine location
        for i = 1:nturbines

            # iterate through each region
            for k = 1:nregions

                # number of boundary vertices
                nvertices = size(boundary_vertices[k])[1]

                # add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
                boundary_vertices_closed = [boundary_vertices[k]; boundary_vertices[k][1,1] boundary_vertices[k][1,2]]

                # initialize intersection counter
                intersection_counter = 0

                # get vector from turbine to the first vertex in first face
                turbine_to_first_facepoint = boundary_vertices_closed[1, :] - [turbine_x[i]; turbine_y[i]]

                # iterate through each boundary
                for j = 1:nvertices

                    # check if y-coordinate of turbine is less than at least one y-coordinate of the two boundary vertices
                    if !(boundary_vertices_closed[j, 2] < turbine_y[i] && boundary_vertices_closed[j+1, 2] < turbine_y[i])        # (this might not be necessary)

                        # check if x-coordinate of turbine is between the x-coordinates of the two boundary vertices
                        if boundary_vertices_closed[j, 1] < turbine_x[i] < boundary_vertices_closed[j+1, 1] || boundary_vertices_closed[j, 1] > turbine_x[i] > boundary_vertices_closed[j+1, 1]

                            # check to see if the turbine is below the boundary
                            if turbine_y[i] < (boundary_vertices_closed[j+1, 2] - boundary_vertices_closed[j, 2]) / (boundary_vertices_closed[j+1, 1] - boundary_vertices_closed[j, 1]) * (turbine_x[i] - boundary_vertices_closed[j, 1]) + boundary_vertices_closed[j, 2]

                                # the vertical ray intersects the boundary
                                intersection_counter += 1

                            end

                        end

                    end

                end

                # check if the turbine is in the current region
                if mod(intersection_counter, 2) == 1

                    # initialize array to hold distances from each turbine to closest boundary face
                    turbine_to_face_distance = zeros(typeof(turbine_x[1]),(nvertices))

                    # iterate through each boundary
                    for j = 1:nvertices

                        # define the vector from the turbine to the second point of the face
                        turbine_to_second_facepoint = boundary_vertices_closed[j+1, :] - [turbine_x[i]; turbine_y[i]]

                        # find perpendicular distance from turbine to current face (vector projection)
                        boundary_vector = boundary_vertices_closed[j+1, :] - boundary_vertices_closed[j, :]

                        # check if perpendicular distance is the shortest
                        if sum(boundary_vector .* -turbine_to_first_facepoint) > 0 && sum(boundary_vector .* turbine_to_second_facepoint) > 0

                            # perpendicular distance from turbine to face
                            turbine_to_face_distance[j] = abs(sum(turbine_to_first_facepoint .* boundary_normals[k][j,:]))

                        # check if distance to first facepoint is shortest
                        elseif sum(boundary_vector .* -turbine_to_first_facepoint) <= 0

                            # distance from turbine to first facepoint
                            turbine_to_face_distance[j] = sqrt(sum(turbine_to_first_facepoint.^2))

                        # distance to second facepoint is shortest
                        else

                            # distance from turbine to second facepoint
                            turbine_to_face_distance[j] = sqrt(sum(turbine_to_second_facepoint.^2))

                        end

                        # reset for next face iteration
                        turbine_to_first_facepoint = turbine_to_second_facepoint        # (for efficiency, so we don't have to recalculate for the same vertex twice)

                    end

                    # magnitude of the constraint value
                    c[i] = ff.smooth_max(-turbine_to_face_distance, s=smooth_max_factor)
                    status[i] = 1

                end

            end

            # check if the turbine is in none of the regions
            if status[i] == 0

                # initialize array to hold distances from each turbine to closest boundary face
                turbine_to_face_distance = zeros(typeof(turbine_x[1]),0)

                # iterate through each region
                for k = 1:nregions

                    # number of boundary vertices
                    nvertices = size(boundary_vertices[k])[1]

                    # add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
                    boundary_vertices_closed = [boundary_vertices[k]; boundary_vertices[k][1,1] boundary_vertices[k][1,2]]

                    # get vector from turbine to the first vertex in first face
                    turbine_to_first_facepoint = boundary_vertices_closed[1, :] - [turbine_x[i]; turbine_y[i]]

                    # iterate through each boundary
                    for j = 1:nvertices

                        # define the vector from the turbine to the second point of the face
                        turbine_to_second_facepoint = boundary_vertices_closed[j+1, :] - [turbine_x[i]; turbine_y[i]]

                        # find perpendicular distance from turbine to current face (vector projection)
                        boundary_vector = boundary_vertices_closed[j+1, :] - boundary_vertices_closed[j, :]

                        # check if perpendicular distance is the shortest
                        if sum(boundary_vector .* -turbine_to_first_facepoint) > 0 && sum(boundary_vector .* turbine_to_second_facepoint) > 0

                            # perpendicular distance from turbine to face
                            push!(turbine_to_face_distance, abs(sum(turbine_to_first_facepoint .* boundary_normals[k][j,:])))

                        # check if distance to first facepoint is shortest
                        elseif sum(boundary_vector .* -turbine_to_first_facepoint) <= 0

                            # distance from turbine to first facepoint
                            push!(turbine_to_face_distance, sqrt(sum(turbine_to_first_facepoint.^2)))

                        # distance to second facepoint is shortest
                        else

                            # distance from turbine to second facepoint
                            push!(turbine_to_face_distance, sqrt(sum(turbine_to_second_facepoint.^2)))

                        end

                        # reset for next face iteration
                        turbine_to_first_facepoint = turbine_to_second_facepoint        # (for efficiency, so we don't have to recalculate for the same vertex twice)

                    end

                end

                # magnitude of the constraint value
                c[i] = -ff.smooth_max(-turbine_to_face_distance, s=smooth_max_factor)
                status[i] = 1

            end

        end

        return c

    end

end

"""
    winding_number_boundary(boundary_vertices,boundary_normals,turbine_x,turbine_y;discrete=false,smooth_max_factor=20.0)

Calculate the distance from each turbine to the nearest point on the boundary.
Negative means the turbine is inside the boundary.

# Arguments
- `boundary_vertices::Array{Float,2}`: vertices of the boundary CCW in order s.t.
        boundaryVertices[i] -> first point of face for unit_normals[i]
- `boundary_normals::Array{Float,2}`: unit normal vector for each boundary face
        CCW where boundaryVertices[i] is the first point of the corresponding face
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
- `discrete::Bool`: true if there are multiple discrete regions, false if there is only one region
- `smooth_max_factor::Float`: factor to control smoothness for smoothmax. The lower the number, the smoother the constraint function
"""
function winding_number_boundary_1(boundary_vertices,boundary_normals,turbine_x,turbine_y;discrete=false,smooth_max_factor=100.0)

    # number of turbines and boundary vertices
    nturbines = length(turbine_x)
    nvertices = size(boundary_vertices)[1]

    # initialize constraint output values
    c = zeros(typeof(turbine_x[1]),(nturbines))

    # initialize array to hold distances from each turbine to closest boundary face
    turbine_to_face_distance = zeros(typeof(turbine_x[1]),(nvertices))

    # add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
    boundary_vertices = [boundary_vertices; boundary_vertices[1,1] boundary_vertices[1,2]]
    boundary_lengths = zeros(nvertices)
    for i = 1:nvertices
        boundary_lengths[i] = sqrt(sum(boundary_vertices[i+1,:] - boundary_vertices[i,:]).^2)
    end

    # iterate through each turbine
    for i = 1:nturbines

        # initialize angle sum
        angle_sum = 0

        # get vector and distance from turbine to the first vertex in first face
        turbine_to_first_facepoint = boundary_vertices[1, :] - [turbine_x[i]; turbine_y[i]]
        turbine_to_first_facepoint_length = sqrt(sum(turbine_to_first_facepoint.^2))

        # iterate through each boundary
        for j = 1:nvertices

            # define the vector and distance from the turbine to the second point of the face
            turbine_to_second_facepoint = boundary_vertices[j+1, :] - [turbine_x[i]; turbine_y[i]]
            turbine_to_second_facepoint_length = sqrt(sum(turbine_to_second_facepoint.^2))

            # find angle
            theta = acos(sum(turbine_to_first_facepoint.*turbine_to_second_facepoint)/(turbine_to_first_facepoint_length*turbine_to_second_facepoint_length))

            # determine sign of angle using cross product
            clockwise = turbine_to_first_facepoint[1]*turbine_to_second_facepoint[2] - turbine_to_first_facepoint[2]*turbine_to_second_facepoint[1] < 0

            # add angle to total angle sum
            if clockwise
                angle_sum -= theta
            else
                angle_sum += theta
            end

            # find perpendicular distance from turbine to face (angle projection)
            perpendicular_distance = turbine_to_first_facepoint_length*turbine_to_second_facepoint_length*sin(theta)/boundary_lengths[j]

            # find the minimum distance to the face
            turbine_to_face_distance[j] = min(turbine_to_first_facepoint_length, turbine_to_second_facepoint_length, perpendicular_distance)

            # reset for next face iteration
            turbine_to_first_facepoint = turbine_to_second_facepoint        # (for efficiency, so we don't have to recalculate for the same vertex twice)
            turbine_to_first_facepoint_length = turbine_to_second_facepoint_length

        end

        # magnitude of the constraint value
        c[i] = -ff.smooth_max(-turbine_to_face_distance, s=smooth_max_factor)
        # c[i] = maximum(turbine_to_face_distance)

        # determine winding number
        winding_number = Int(round(angle_sum/(2*pi)))

        if winding_number != 0
            c[i] = -c[i]
        end

    end

    return c

end

"""
    winding_number_boundary(boundary_vertices,boundary_normals,turbine_x,turbine_y;discrete=false,smooth_max_factor=20.0)

Calculate the distance from each turbine to the nearest point on the boundary.
Negative means the turbine is inside the boundary.

# Arguments
- `boundary_vertices::Array{Float,2}`: vertices of the boundary CCW in order s.t.
        boundaryVertices[i] -> first point of face for unit_normals[i]
- `boundary_normals::Array{Float,2}`: unit normal vector for each boundary face
        CCW where boundaryVertices[i] is the first point of the corresponding face
- `turbine_x::Array{Float}`: turbine x locations
- `turbine_y::Array{Float}`: turbine y locations
- `discrete::Bool`: true if there are multiple discrete regions, false if there is only one region
- `smooth_max_factor::Float`: factor to control smoothness for smoothmax. The lower the number, the smoother the constraint function
"""
function winding_number_boundary_2(boundary_vertices,boundary_normals,turbine_x,turbine_y;discrete=false,smooth_max_factor=100.0)

    # number of turbines and boundary vertices
    nturbines = length(turbine_x)
    nvertices = size(boundary_vertices)[1]

    # initialize constraint output values
    c = zeros(typeof(turbine_x[1]),(nturbines))

    # initialize array to hold distances from each turbine to closest boundary face
    turbine_to_face_distance = zeros(typeof(turbine_x[1]),(nvertices))

    # add the first boundary vertex again to the end of the boundary vertices vector (to form a closed loop)
    boundary_vertices = [boundary_vertices; boundary_vertices[1,1] boundary_vertices[1,2]]

    # iterate through each turbine
    for i = 1:nturbines

        # initialize winding number
        winding_number = 0

        # get vector from turbine to the first vertex in first face
        turbine_to_first_facepoint = boundary_vertices[1, :] - [turbine_x[i]; turbine_y[i]]

        # iterate through each boundary
        for j = 1:nvertices

            # check if x-coordinate of turbine is greater than the x-coordinate of the first boundary vertex
            if boundary_vertices[j, 1] <= turbine_x[i]

                # check if x coordinate of turbine is less than x-coordinate of second boundary vertex
                if turbine_x[i] < boundary_vertices[j+1, 1]

                    # check for intersection of ray and left-to-right boundary
                    if (boundary_vertices[j+1, 1] - boundary_vertices[j, 1]) * (turbine_y[i] - boundary_vertices[j, 2]) - (turbine_x[i] - boundary_vertices[j, 1]) * (boundary_vertices[j+1, 2] - boundary_vertices[j, 2]) > 0
                        winding_number += 1
                    end
                end

            # x-coordinate of turbine is less than the x-coordinate of the first boundary vertex
            else

                # check if x coordinate of turbine is greater than x-coordinate of second boundary vertex
                if turbine_x[i] > boundary_vertices[j+1, 1]

                    # check for intersection of ray and right-to-left boundary
                    if (boundary_vertices[j+1, 1] - boundary_vertices[j, 1]) * (turbine_y[i] - boundary_vertices[j, 2]) - (turbine_x[i] - boundary_vertices[j, 1]) * (boundary_vertices[j+1, 2] - boundary_vertices[j, 2]) < 0
                        winding_number -= 1
                    end
                end
            end

            # define the vector from the turbine to the second point of the face
            turbine_to_second_facepoint = boundary_vertices[j+1, :] - [turbine_x[i]; turbine_y[i]]

            # find perpendicular distance from turbine to current face (vector projection)
            boundary_vector = boundary_vertices[j+1, :] - boundary_vertices[j, :]

            # check if perpendicular distance is the shortest
            if sum(boundary_vector .* -turbine_to_first_facepoint) > 0 && sum(boundary_vector .* turbine_to_second_facepoint) > 0

                # perpendicular distance from turbine to face
                turbine_to_face_distance[j] = abs(sum(turbine_to_first_facepoint .* boundary_normals[j,:]))

            # check if distance to first facepoint is shortest
            elseif sum(boundary_vector .* -turbine_to_first_facepoint) <= 0

                # distance from turbine to first facepoint
                turbine_to_face_distance[j] = sqrt(sum(turbine_to_first_facepoint.^2))

            # distance to second facepoint is shortest
            else

                # distance from turbine to second facepoint
                turbine_to_face_distance[j] = sqrt(sum(turbine_to_second_facepoint.^2))

            end

            # reset for next face iteration
            turbine_to_first_facepoint = turbine_to_second_facepoint        # (for efficiency, so we don't have to recalculate for the same vertex twice)

        end

        # magnitude of the constraint value
        c[i] = -ff.smooth_max(-turbine_to_face_distance, s=smooth_max_factor)

        # sign of the constraint value (- is inside, + is outside)
        if winding_number != 0
            c[i] = -c[i]
        end

    end

    return c

end

"""
    VR_bounary_startup(bndry_x_clsd, bndry_y_clsd, start_dist, turb_min_spacing, num_turbs)

Determines if the requested number of turbines can be placed along the closed
boundary with spacing and corner constraints. If the requested <num_turbs> is
too many, places as many turbines as possible along the boundary, and returns
the number of turbines not placed. NOTE: A shortcoming is that the
smallest-angled corner limits the spacing of all turbines. in the worst case,
a very thin boundary area would prevent any more than one turbine being placed
on the boundary, though more would be optimal. Future work would check to make
sure this corner (and the length of its adjacent sides) don't actually require
limiting the minimum distance between turbines.

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
- `start_dist::Float64`: the distance (positive or negative) along the boundary
        from the first boundary point where the turbines will begin to be placed
- `turb_spacing::Float64`: the fixed distance along the boundary's edge between
        adjacent turbines
- 'num_turbs::Float64`: the number of turbines to be placed around the boundary.
        Note that this function assumes VR_bounary_startup() has already been
        run so that the user won't attempt to place too many turbines.
- 'bndry_seg_length::Array{Int}`: an array of the lengths between adjacent
        boundary verticies, corresponding to how they appear in bndry_x and _y
"""
function VR_boundary_startup(bndry_x_clsd, bndry_y_clsd, start_dist, turb_min_spacing, num_turbs)
    #- Get arc-length -#
    bndry_seg_length = getPerimeterLength(bndry_x_clsd,bndry_y_clsd)
    bndry_tot_length = sum(bndry_seg_length)

    #- Determine how much distance between turbs for worst case (shallow corner) -#
    limiting_angle = calcSmallestAngle(bndry_x_clsd, bndry_y_clsd)  # The smallest angle in our boundary
    turb_corner_spacing = turb_min_spacing/sind(limiting_angle/2.0) # Min spacing neeed for the smallest angle
    if (turb_corner_spacing > turb_min_spacing)   # If the corners are tighter than our minimum distance
        turb_min_spacing = turb_corner_spacing    # That dictates our new minimum distance
    end

    #- Determine how many turbines we can actually place with corner constraint -#
    max_num_turbs = floor(Int, bndry_tot_length/turb_min_spacing) # Get max number of turbines that can fit the perimeter
    num_leftover_turbs = 0
    if (max_num_turbs < num_turbs)      # If the max is less than the number allocated
        num_leftover_turbs = floor(Int, num_turbs - max_num_turbs) # Note how many turbines we aren't placing
        num_turbs = floor(Int, max_num_turbs)       # Place as many as will fit
    end

    #- Evenly space the turbines we'll place -#
    turb_min_spacing = bndry_tot_length / num_turbs
    #- Make sure our start point is within the boundary length bounds
    start_dist = mod(start_dist, bndry_tot_length)
    #- Place the correct number of turbs along the boundary <start_dist> away from the first vertex -#
    turbine_x, turbine_y = VR_boundary(bndry_x_clsd, bndry_y_clsd, start_dist, turb_min_spacing, num_turbs)

    # Return the x- and y- coordinates of every turbine, and how many weren't placed
    return turbine_x,turbine_y, num_leftover_turbs
end

"""
    VR_boundary(bndry_x_clsd, bndry_y_clsd, start_dist, turb_spacing, num_turbs, bndry_seg_length)

Uses the Boundary portion of Boundary-Grid variable reduction method
place turbines along a closed wind farm boundary and perturb their location with
one (1) variable <start_dist>.  NOTE: Use of this function assumes prior use of
VR_bounary_startup(), which ensures the number of turbines placed on the
boundary doesn't violate any minimum spacing rules eiter along the boundary or
around corners.

# Arguments
- `bndry_x::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
- `start_dist::Float64`: the distance (positive or negative) along the boundary
        from the first boundary point where the turbines will begin to be placed
- `turb_spacing::Float64`: the fixed distance along the boundary's edge between
        adjacent turbines
- 'num_turbs::Float64`: the number of turbines to be placed around the boundary.
        Note that this function assumes VR_bounary_startup() has already been
        run so that the user won't attempt to place too many turbines.
"""
function VR_boundary(bndry_x_clsd, bndry_y_clsd, start_dist, turb_spacing, num_turbs)
    # Initialize necessary variables
    bndry_seg_length = getPerimeterLength(bndry_x_clsd,bndry_y_clsd)
    num_segs = length(bndry_x_clsd)-1
    # Initialize turbine locations
    turbine_x = zeros(num_turbs)
    turbine_y = zeros(num_turbs)

    #- Figure out where the starting point should be -#
    # "leg" is distance until next placed turbine
    # "seg" is distance between boundary verticies
    leg_remaining = start_dist
    curr_seg = 1
    for i in 1:num_segs   # Looping through the segments till we get there
        curr_seg = i
        if (bndry_seg_length[i] < leg_remaining)  # If this segment length is less than the start length
            leg_remaining -= bndry_seg_length[i]  # Clip the start length and move to the next one
        else                                        # Otherwise the start point is on this segment
            percent_to_start = leg_remaining / bndry_seg_length[i] # How far along our segment we are
            # Translate how far along the segment we are to actual x- y-coords
            turbine_x[1] = bndry_x_clsd[i] + ((bndry_x_clsd[i+1] - bndry_x_clsd[i]) * percent_to_start)
            turbine_y[1] = bndry_y_clsd[i] + ((bndry_y_clsd[i+1] - bndry_y_clsd[i]) * percent_to_start)
            break
        end
    end
    # Get how much distance is left on this leg after the starting point
    percent_left_of_segment = 1 - abs((turbine_x[1] - bndry_x_clsd[curr_seg]) / (bndry_x_clsd[curr_seg+1] - bndry_x_clsd[curr_seg]))
    seg_remaining = bndry_seg_length[curr_seg] * percent_left_of_segment
    leg_remaining = turb_spacing

    #- Place the rest of the turbines -#
    for i in 2:num_turbs    # For every turbine we have to place
        if(seg_remaining > leg_remaining)   #- If there's space on this leg to place the next turbine
            percent_to_place = (bndry_seg_length[curr_seg] - seg_remaining + turb_spacing) / bndry_seg_length[curr_seg]
            seg_remaining -= turb_spacing # Take out the distance we used
        else                                #- If there isn't enough space on this boundary segment to place the next turbine
            while(seg_remaining < leg_remaining)           # Keep checking boundary segments until we reach the next placement
                leg_remaining -= seg_remaining             # Subtract that much till our next placement
                curr_seg = mod(curr_seg, (num_segs))+1   # Increment what segment we're onspot
                seg_remaining = bndry_seg_length[curr_seg] # Reset how much segment length we have left
            end
            percent_to_place = leg_remaining / bndry_seg_length[curr_seg]
            seg_remaining = bndry_seg_length[curr_seg] - leg_remaining
        end
        #- Place the turbines the appropriate distance from the segment start point -#
        turbine_x[i] = bndry_x_clsd[curr_seg] + ((bndry_x_clsd[curr_seg+1] - bndry_x_clsd[curr_seg]) * percent_to_place)
        turbine_y[i] = bndry_y_clsd[curr_seg] + ((bndry_y_clsd[curr_seg+1] - bndry_y_clsd[curr_seg]) * percent_to_place)
        leg_remaining = turb_spacing  # Reset how much length till the next turbine is placed
    end

    return turbine_x, turbine_y
end

"""
    iea37cs4BndryVRIntPM(bndry_x_clsd, bndry_y_clsd, bndry_corner_indcies, turbine_x, turbine_y, turb_diam, turb_min_space, num_turbs_to_place)

Uses the Variable reduction method for placing boundary turbines, and the
Partition Method (from splined_boundary()) for random interior points,
maintaining proper spacing from all previously placed turbines.

# Arguments
- `bndry_x_clsd::Array{Float,1}` : 1-D array of x-coordinates for the vertices
        around a singlar closed boundary
- `bndry_y_clsd::Array{Float,1}` : 1-D array of y-coordinates for the vertices
        around a singlar closed boundary
- `bndry_corner_indcies::Float64`: The indicies within <bndry_x_clsd> and
        <bndry_y_clsd> which denote the "corners" adjacent turbines
- 'turb_min_space::Float64`: For proximity knowledge, the minimum spacing
        required between any two turbines
- 'num_bndry_turbs::Float64`: The number of turbines desired to be placed along
        the boundary. If too many are selected (due to spacing condtraints), the
        remaining will be placed in the interior
- 'num_tot_turbs::Float64`: The number of total turbines to be placed both on
        the boundary and in the interior
"""
function iea37cs4BndryVRIntPM(bndry_x_clsd, bndry_y_clsd, bndry_corner_indicies, turb_min_space, num_bndry_turbs, num_tot_turbs)
    #-- Place all the boundary turbines we can --#
    bndry_tot_len = sum(getPerimeterLength(bndry_x_clsd,bndry_y_clsd))
    #- Make a random starting point along the boundary -#
    start_dist = rand(Float64) * bndry_tot_len
    #- Place the boundary turbines -#
    turbine_x_bndry, turbine_y_bndry, num_leftover_turbs = VR_boundary_startup(bndry_x_clsd, bndry_y_clsd, start_dist, turb_min_space, num_bndry_turbs)
    #- Determine how many will be placed in the interior -#
    num_bndry_turbs = num_bndry_turbs - num_leftover_turbs
    num_interior_turbs = num_tot_turbs - num_bndry_turbs

    # Initialize full list of turbine locations
    turbine_x = zeros(num_tot_turbs)
    turbine_y = zeros(num_tot_turbs)
    #- Fill in the ones ew've already placed along the boundary
    turbine_x[1:num_bndry_turbs] = turbine_x_bndry
    turbine_y[1:num_bndry_turbs] = turbine_y_bndry

    #-- Initialize interior space --#w
    num_sides = length(bndry_corner_indicies)-1
    #- Get the x-values -#
    x_min_indx = 3                      # Default to work w/ squared boundaries
    if num_sides == 3                   # If we only have 3 corners
        x_min_indx = 2                  # denote that it's a triangle boundary
    end
    x_max = bndry_x_clsd[bndry_corner_indicies[1]]           # Our maximum x-value
    x_min = bndry_x_clsd[bndry_corner_indicies[x_min_indx]]  # Our min x-value
    turbine_x[num_bndry_turbs+1:end] = (x_max - x_min)*rand(Float64, num_interior_turbs) .+ x_min # Get random x-values for the interior turbines

    #-- Get the y-values --#
    #- Determine the upper and lower splines to use for the given x -#
    # Fake for-loop here for proximity checking.
    i = num_bndry_turbs + 1   # Start after our last boundary turbine
    while (i <= num_tot_turbs)
        # Get the max and min y-value for this x-value
        y_min, y_max = getUpDwnYvals(turbine_x[i], bndry_x_clsd, bndry_y_clsd, bndry_corner_indicies)
        # Generate a random y-value within these limits
        turbine_y[i] = (y_max - y_min)*rand(Float64) + y_min # Get a random number in our bounds
        #- Check it doesn't conflict with aleady placed turbines -#
        for j in 1:(i-1) # Check the new ones we've place so far
            # If this turbine has a proximity conflict
            if (coordDist(turbine_x[i], turbine_y[i], turbine_x[j], turbine_y[j]) < turb_min_space)
                turbine_x[i] = (x_max - x_min)*rand(Float64) + x_min # Give it a new x-val
                i = i-1 # Redo the y-val too
                break # Stop checking for conflicts and redo the y-value
            end
        end
        i += 1
    end

    return turbine_x, turbine_y, num_interior_turbs
end


text:
"""
parse_tuple(ps, ret)

`ret` is followed by a comma so tries to parse the rest of the
tuple.
"""
function parse_tuple end

@static if VERSION > v"1.1-"
    function parse_tuple(ps::ParseState, @nospecialize(ret))
        op = PUNCTUATION(next(ps))
        if (isassignment(ps.nt) && ps.nt.kind != Tokens.APPROX)
            ret = EXPR{TupleH}(Any[ret, op])
        elseif closer(ps)
            ret = ErrorToken(EXPR{TupleH}(Any[ret, op]))
        else
            nextarg = @closer ps tuple parse_expression(ps)
            ret = EXPR{TupleH}(Any[ret, op, nextarg])
        end
        return ret
    end

    function parse_tuple(ps::ParseState, ret::EXPR{TupleH})
        op = PUNCTUATION(next(ps))
        if (isassignment(ps.nt) && ps.nt.kind != Tokens.APPROX)
            push!(ret, op)
        elseif closer(ps)
            push!(ret, op)
            ret = ErrorToken(ret)
        else
            nextarg = @closer ps tuple parse_expression(ps)
            if !(is_lparen(first(ret.args)))
                push!(ret, op)
                push!(ret, nextarg)
            else
                ret = EXPR{TupleH}(Any[ret, op, nextarg])
            end
        end
        return ret
    end
else
    function parse_tuple(ps::ParseState, @nospecialize(ret))
        op = PUNCTUATION(next(ps))
        if closer(ps) || (isassignment(ps.nt) && ps.nt.kind != Tokens.APPROX)
            ret = EXPR{TupleH}(Any[ret, op])
        else
            nextarg = @closer ps tuple parse_expression(ps)
            ret = EXPR{TupleH}(Any[ret, op, nextarg])
        end
        return ret
    end

    function parse_tuple(ps::ParseState, ret::EXPR{TupleH})
        op = PUNCTUATION(next(ps))
        if closer(ps) || (isassignment(ps.nt) && ps.nt.kind != Tokens.APPROX)
            push!(ret, op)
        else
            nextarg = @closer ps tuple parse_expression(ps)
            if !(is_lparen(first(ret.args)))
                push!(ret, op)
                push!(ret, nextarg)
            else
                ret = EXPR{TupleH}(Any[ret, op, nextarg])
            end
        end
        return ret
    end
end

"""
    parse_array(ps)
Having hit '[' return either:
+ A vect
+ A vcat
+ A comprehension
+ An array (vcat of hcats)
"""
function parse_array(ps::ParseState)
    args = Any[PUNCTUATION(ps)]

    if ps.nt.kind == Tokens.RSQUARE
        accept_rsquare(ps, args)
        ret = EXPR{Vect}(args)
    else
        first_arg = @nocloser ps newline @closesquare ps  @closer ps insquare @closer ps ws @closer ps wsop @closer ps comma parse_expression(ps)

        if ps.nt.kind == Tokens.RSQUARE
            if first_arg isa EXPR{Generator} || first_arg isa EXPR{Flatten}
                accept_rsquare(ps, args)

                if first_arg.args[1] isa BinaryOpCall && is_pairarrow(first_arg.args[1].op)
                    return EXPR{DictComprehension}(Any[args[1], first_arg, INSTANCE(ps)])
                else
                    return EXPR{Comprehension}(Any[args[1], first_arg, INSTANCE(ps)])
                end
            elseif ps.ws.kind == SemiColonWS
                push!(args, first_arg)
                accept_rsquare(ps, args)
                return EXPR{Vcat}(args)
            else
                push!(args, first_arg)
                accept_rsquare(ps, args)
                ret = EXPR{Vect}(args)
            end
        elseif ps.nt.kind == Tokens.COMMA
            etype = Vect
            push!(args, first_arg)
            accept_comma(ps, args)
            @closesquare ps parse_comma_sep(ps, args, false)
            accept_rsquare(ps, args)
            return EXPR{etype}(args)
        elseif ps.ws.kind == NewLineWS
            ret = EXPR{Vcat}(args)
            push!(ret, first_arg)
            while ps.nt.kind != Tokens.RSQUARE
                if ps.nt.kind == Tokens.ENDMARKER
                    break
                end
                a = @closesquare ps  parse_expression(ps)
                push!(ret, a)
            end
            accept_rsquare(ps, ret)
            update_span!(ret)
            return ret
        elseif ps.ws.kind == WS || ps.ws.kind == SemiColonWS
            first_row = EXPR{Hcat}(Any[first_arg])
            while ps.nt.kind != Tokens.RSQUARE && ps.ws.kind != NewLineWS && ps.ws.kind != SemiColonWS
                if ps.nt.kind == Tokens.ENDMARKER
                    break
                end
                a = @closesquare ps @closer ps ws @closer ps wsop parse_expression(ps)
                push!(first_row, a)
            end
            if ps.nt.kind == Tokens.RSQUARE && ps.ws.kind != SemiColonWS
                if length(first_row.args) == 1
                    first_row = EXPR{Vcat}(first_row.args)
                end
                push!(first_row, INSTANCE(next(ps)))
                pushfirst!(first_row, args[1])
                update_span!(first_row)
                return first_row
            else
                if length(first_row.args) == 1
                    first_row = first_row.args[1]
                else
                    first_row = EXPR{Row}(first_row.args)
                end
                ret = EXPR{Vcat}(Any[args[1], first_row])
                while ps.nt.kind != Tokens.RSQUARE
                    if ps.nt.kind == Tokens.ENDMARKER
                        break
                    end
                    first_arg = @closesquare ps @closer ps ws @closer ps wsop parse_expression(ps)
                    push!(ret, EXPR{Row}(Any[first_arg]))
                    while ps.nt.kind != Tokens.RSQUARE && ps.ws.kind != NewLineWS && ps.ws.kind != SemiColonWS
                        if ps.nt.kind == Tokens.ENDMARKER
                            break
                        end
                        a = @closesquare ps @closer ps ws @closer ps wsop parse_expression(ps)
                        push!(last(ret.args), a)
                    end
                    # if only one entry dont use :row
                    if length(last(ret.args).args) == 1
                        ret.args[end] = ret.args[end].args[1]
                    end
                    update_span!(ret)
                end
                accept_rsquare(ps, ret)
                update_span!(ret)
                return ret
            end
        else
            ret = EXPR{Vect}(args)
            push!(ret, first_arg)
            push!(ret, accept_rsquare(ps))
        end
    end
    return ret
end

"""
    parse_ref(ps, ret)

Handles cases where an expression - `ret` - is followed by
`[`. Parses the following bracketed expression and modifies it's
`.head` appropriately.
"""
function parse_ref(ps::ParseState, ret)
    next(ps)
    ref = parse_array(ps)
    return _parse_ref(ret, ref)
end

function _parse_ref(ret, ref::EXPR{Vect})
    args = Any[ret]
    for a in ref.args
        push!(args, a)
    end
    return EXPR{Ref}(args)
end

function _parse_ref(ret, ref::EXPR{Hcat})
    args = Any[ret]
    for a in ref.args
        push!(args, a)
    end
    return EXPR{TypedHcat}(args)
end

function _parse_ref(ret, ref::EXPR{Vcat})
    args = Any[ret]
    for a in ref.args
        push!(args, a)
    end
    return EXPR{TypedVcat}(args)
end

function _parse_ref(ret, ref)
    args = Any[ret]
    for a in ref.args
        push!(args, a)
    end
    return EXPR{TypedComprehension}(args)
end

"""
parse_curly(ps, ret)

Parses the juxtaposition of `ret` with an opening brace. Parses a comma
seperated list.
"""
function parse_curly(ps::ParseState, ret)
    args = Any[ret, PUNCTUATION(next(ps))]
    parse_comma_sep(ps, args, true)
    accept_rbrace(ps, args)
    return EXPR{Curly}(args)
end

function parse_braces(ps::ParseState)
    args = Any[PUNCTUATION(ps)]
    parse_comma_sep(ps, args, true)
    accept_rbrace(ps, args)
    return EXPR{Braces}(args)
end


text:
module UIVariables

#------------------------------------------------------------------------------
# EXports
#------------------------------------------------------------------------------
export  BasicValidation,
        Variable,
        on,
        emit,
        addGridColumn!,
        insertGridColumn!,
        deleteGridColumn!,
        replaceGridColumn!,
        gridOption!,
        gridColumnOption!,
        forEachGridColumn!,
        clearGridColumns!,
        gridDefaultColDef!,
        gridDefaultColDef!

using ..Glimmer
import Tables
import JSON, JSONTables
import Dates
using OrderedCollections
using Parameters
import UUIDs

function typedict(x::Union{AbstractUIVariable, AbstractUIValidation})
    res = OrderedDict()

    for fn in filter(x->!startswith(string(x), "_"), fieldnames(typeof(x)))
        key = fn
        val = getfield(x, fn)
        if (val === nothing)
            continue
        end
        if (isstructtype(typeof(val)) && !(val isa String))
            if (val isa Dict || val isa OrderedDict)
                # @info "DICT" key, val
                res[key] = val;
            else
                # @show key
                res[key] = typedict(val)
            end
        else
            res[key] = val;
        end
    end
    # if (isstructtype(typeof(x)))
    # Dict(fn=>getfield(x, fn) for fn ∈ filter(x->!startswith(string(x), "_"), fieldnames(typeof(x))))

    return res;
end

# ----------------------------------------------------------------------
# Validation
# ----------------------------------------------------------------------
@with_kw mutable struct BasicValidation <: AbstractUIValidation
    type::String = "any"
    min::Union{Nothing, Float64} = nothing
    max::Union{Nothing, Float64} = nothing
    valid_values::Union{Nothing, Vector{Any}} = nothing
    invalid_values::Union{Nothing, Vector{Any}} = nothing
end

function BasicValidation(min::Any, max::Any)
    # _, min, max = promote(0.0, min, max)
    return BasicValidation(type="number", min=min, max=max)
end

# ----------------------------------------------------------------------
# Variables
# ----------------------------------------------------------------------
@with_kw mutable struct Variable <: AbstractUIVariable
    name::String = "defaultVariableName"
    type::String = "any"
    value::Union{Nothing, Any} = nothing;
    validation::Union{Nothing, BasicValidation} = nothing

    _app::Union{Nothing, AbstractUIApp} = nothing

    _events_handlers::Dict{Symbol, Vector{Function}} = Dict{Symbol, Vector{Function}}()
    _meta_data::Dict{Symbol, Any} = Dict{Symbol, Any}()
end

"""
    rawValue(var::Variable) -> string

returns the internal value of the variable. in most cases this value would be identical to using var[], but in
some cases, such as `aggrid` type variable, this would return the internal structure of the grid definition, allowing the user
to manually manipulate it for better control on the grid behevior.
"""
function rawValue(var::Variable)
    return var.value;
end

# function on(func::Function, var::AbstractUIVariable)
#     # var._on_change_func = func
#     on(func, var, :valueChanged)
# end
"""
    on(func::Function, var::AbstractUIVariable, event::Symbol = :valueChanged)

Adds function `func` as listener to the variable. Whenever variables's value is set via var[] = val or one of
the defined UI Controls, func is called with val. Most variables have just one type of events associated with them: `valueChanged`.
But some, such as 'aggrid' have additional events that you can bind listeners to.
"""
function on(func::Function, var::AbstractUIVariable, event::Symbol = :valueChanged)
    if (!haskey(var._events_handlers, event))
        var._events_handlers[event] = Vector{Function}()
    end
    push!(var._events_handlers[event], func)
end

"""
    on(func::Function, var::AbstractUIVariable, event::Symbol = :valueChanged)

Emits the specific event on a variable, calling all the listeners.
"""
function emit(var::AbstractUIVariable, event::Symbol, val::Any = var.value)
    if (haskey(var._events_handlers, event))
        funcs = var._events_handlers[event]
        for func in funcs
            func(val)
        end
    end
end

function Base.getindex(var::Variable)
    if (var.type == "integer" || var.type == "int")
        if (isa(var.value, String))
            try
                return Int(floor(parse(Float64, var.value)))
            catch
                return 0.0
            end
        elseif (isa(var.value, Int64))
            return var.value
        elseif (isa(var.value, Float64))
            return Int(floor(var.value))
        else
            return var.value
        end
    elseif (var.type == "float64")
        if (isa(var.value, String))
            try
                return parse(Float64, var.value)
            catch
                return 0.0
            end
        elseif (isa(var.value, Int))
            return convert(Float64, var.value)
        elseif (isa(var.value, Float64))
            return var.value
        else
            return var.value
        end
    elseif (var.type == "aggrid")
        jtable = JSONTables.jsontable(JSON.json(var.value["rowData"]))
        return var._meta_data[:ag_grid_materializer](jtable)
    else
        return var.value
    end
end

function Base.setindex!(var::Variable, val::Any)
    # @info "setting variable $(var.name) to $(val) - app: $(var._app)"

    var.value = val;

    # special case of agGrid
    if (var.type == "aggrid")
        if !(val isa String)
            # var.value = GridUtils.table2agGrid(val)
            # var._meta_data[:ag_grid_variable] = val;
            # var._meta_data[:ag_grid_materializer] = Tables.materializer(val);
            # var._meta_data[:ag_grid_schema] = Tables.schema(val);
            if (val isa Dict || val isa OrderedDict)
                if (haskey(val, "_metaData"))
                    var._meta_data[:ag_grid_materializer] = val["_metaData"]["ag_grid_materializer"]
                    var._meta_data[:ag_grid_schema] = val["_metaData"]["ag_grid_schema"]
                    delete!(var.value, "_metaData")
                end
            end
        end
    end

    Glimmer.updateVariable!(var._app, var)
end

function Base.setindex!(var::Variable, val::Any, type::String)
    if (var.type == "aggrid")
        if (type == "options")
            var.value["options"] = val;
        elseif (type == "readonly")
            for colDef in var.value["columnDefs"]
                colDef["editable"] = !val
            end
        else
            @warn "unrecognized property of aggrid [$type]"
        end
        return;
    end

    if (type == "png")
        try
            val = Glimmer.base64png(val)
        catch e
            @error "Can't convert to SVG [$e]"
        end
    elseif (type == "svg")
        try
            val = Glimmer.base64svg(val)
        catch e
            @error "Can't convert to SVG [$e]"
        end
    end
    var.value = val;
    Glimmer.updateVariable!(var._app, var)
end

#---------------------------------------------------------
# Tables and AG-Grid Utility Functions
#---------------------------------------------------------

"""
    clearGridColumns!(var::Variable)

Clear all existing definitions of columns for this `var` grid.
"""
function clearGridColumns!(var::Variable)
    rawValue(var)["columnDefs"] = []
end

"""
    addGridColumn!(var::Variable, colDef::Dict)

Add a new grid column definition as the last column.
"""
function addGridColumn!(var::Variable, colDef::Dict)
    push!(rawValue(var)["columnDefs"], colDef)
end

"""
    insertGridColumn!(var::Variable, index::Int64, colDef::Dict)

Add a new grid column definition as the last column.
"""
function insertGridColumn!(var::Variable, index::Int64, colDef::Dict)
    insert!(rawValue(var)["columnDefs"], index, colDef)
end

"""
deleteGridColumn!(var::Variable, index::Int64)

Delete a grid column at a specific index.
"""
function deleteGridColumn!(var::Variable, index::Int64)
    deleteat!(rawValue(var)["columnDefs"], index:index)
end

"""
    deleteGridColumn!(var::Variable, name::String)

Delete a grid column at a specific index.
"""
function deleteGridColumn!(var::Variable, name::String)
    colDefs = rawValue(var)["columnDefs"]
    deleteat!(colDefs, findall(x->x["headerName"]==name || x["field"]==name,colDefs))
end

"""
    replaceGridColumn!(var::Variable, index::Int64, colDef::Dict)

Replace an existing grid column definition with a new one.
"""
function replaceGridColumn!(var::Variable, index::Int64, colDef::Dict)
    rawValue(var)["columnDefs"][index] = colDef
end

"""
    gridDefaultColDef!(var::Variable, colDef::Dict)

Set the default column definition for this grid variable.
"""
function gridDefaultColDef!(var::Variable, colDef::Dict)
    rawValue(var)["defaultColDef"] = colDef
end

"""
    gridOption!(var::Variable, option::String, value)

Set a grid options.
"""
function gridOption!(var::Variable, option::String, value)
    if (!haskey(rawValue(var), "options"))
        rawValue(var)["options"] = Dict();
    end
    rawValue(var)["options"][option] = value;
end

"""
    gridColumnOption!(var::Variable, colName::String, option::String, value)

Set a grid column's option.
"""
function gridColumnOption!(var::Variable, colName::String, option::String, value)
    found = false;
    for col in rawValue(var)["columnDefs"]
        if (haskey(col, "field") && col["field"] == colName)
            col[option] = value;
            found = true;
        end
    end
    if (!found)
        @warn "Can't find grid column with name [$colName]"
    end
end

"""
    forEachGridColumn!(var::Variable, option::String, value)

Set a column options for each exiting column definition in this grid variable.
"""
function forEachGridColumn!(var::Variable, option::String, value)
    for col in rawValue(var)["columnDefs"]
        col[option] = value
    end
end

end # module UIVariables


text:
#!/usr/bin/env julia

module Drivers

using ...ReconstructionEvaluations
using ..Synaptor
using ...ReconstructionEvaluations.io

const RE = ReconstructionEvaluations

"""

    run_synaptor_cfg( cfg, output_prefix, dist_thr, res,
                      keep_seg=false, keep_mapping=false )

Runs Synaptor postprocessing using a passed configuration filename. Also performs
duplicate removal within the passed distance threshold (dist_thr), assuming a voxel
resolution (res).

keep_seg and keep_mapping flags determine whether to keep these products
of the Synaptor postprocessing
"""
function run_synaptor_cfg( cfg, output_prefix, dist_thr, res,
  keep_seg=false, keep_mapping=false )

  @assert isfile(cfg)

  Synaptor.run_cfg_file(cfg)

  if !keep_seg Synaptor.rm_seg_file(output_prefix) end
  if !keep_mapping Synaptor.rm_mapping_file(output_prefix) end

  edge_filename = "$(output_prefix)_edges.csv"
  cons_edge_filename = "$(output_prefix)_edges_cons.csv"

  #The string interpolation is a bit funky, so the resolution
  # vector should be written as a string
  Synaptor.rm_duplicates( edge_filename, cons_edge_filename, dist_thr, string(res) )

  cons_edge_filename
end

"""

    h5_file_overlap_matrix( seg1fname, seg2fname, chunk_shape, verb=true )

Finds the overlap matrix for two h5 segmentations as a sparse matrix
"""
function h5_file_om( seg1fname, seg2fname, chunk_shape, verb=true )

  seg1 = io.read_h5(seg1fname, false)
  seg2 = io.read_h5(seg2fname, false)

  impl_om = RE.overlap_in_chunks(seg1, seg2, chunk_shape, verb)

  RE.om_from_impl(impl_om)
end

"""
Returns the number of splits and mergers for each segment in the rows of the
overlap matrix
"""
function splits_and_mergers( om )

  rs, cs = findn(om)

  splits = Dict{eltype(rs),Vector{eltype(cs)}}();
  mergers = Dict{eltype(cs),Vector{eltype(rs)}}();

  for i in eachindex(rs)
    r,c = rs[i], cs[i]

    if !haskey(splits,r) splits[r] = [c]
    else push!(splits[r],c)
    end

    if !haskey(mergers,c) mergers[c] = [r]
    else push!(mergers[c],r)
    end

  end

  splits  = filter( (k,v) -> length(v) > 1, splits )
  mergers = filter( (k,v) -> length(v) > 1, mergers )

  splits, mergers
end

end#module end


text:
using Polyhedra
using SumOfSquares

"""
    struct PolySet{T, B, U} <: AbstractSet{U}
        degree::Int
        p::GramMatrix{T, B, U}
    end

Set ``\\{\\, x \\mid p(x) \\le 1 \\,\\}`` where `p` is a homogeneous polynomial
of degree `degree`.
"""
struct PolySet{T, B, U} <: AbstractSet{U}
    degree::Int
    p::GramMatrix{T, B, U}
end

function space_variables(set::PolySet)
    return MP.variables(set.p)
end

"""
    struct ConvexPolySet{T, B, U} <: AbstractSet{U}
        degree::Int
        p::GramMatrix{T, B, U}
        convexity_proof::Union{Nothing, SumOfSquares.SymMatrix{T}} # may be nothing after applying LinearMap
    end

Set ``\\{\\, x \\mid p(x) \\le 1 \\,\\}`` where `p` is a homogeneous polynomial
of degree `degree`.
"""
struct ConvexPolySet{T, B, U} <: AbstractSet{U}
    degree::Int
    p::GramMatrix{T, B, U}
    convexity_proof::Union{Nothing, SumOfSquares.SymMatrix{T}} # may be nothing after applying LinearMap
end
function ConvexPolySet(
    degree::Int,
    p::GramMatrix{T, B, U},
    convexity_proof::SumOfSquares.SymMatrix{T}) where {T, B, U}
    return ConvexPolySet{T, B, U}(degree, p, convexity_proof)
end
function ConvexPolySet(
    degree::Int,
    p::GramMatrix{S},
    convexity_proof::SumOfSquares.SymMatrix{T}) where {S, T}
    U = promote_type(S, T)
    _convert(mat) = SumOfSquares.SymMatrix(convert(Vector{U}, mat.Q), mat.n)
    return ConvexPolySet(
        degree, GramMatrix(_convert(p.Q), p.basis), _convert(convexity_proof))
end

function space_variables(set::ConvexPolySet)
    return MP.variables(set.p)
end
function dimension(set::ConvexPolySet)
    return MP.nvariables(set.p)
end
function gauge1(set::ConvexPolySet)
    return set.p
end
function zero_eliminate(set::ConvexPolySet, I)
    vars = space_variables(set)[I]
    K = findall(mono -> all(var -> iszero(MP.degree(mono, var)), vars),
                set.p.basis.monomials)
    Q = SumOfSquares.square_getindex(set.p.Q, K)
    monos = set.p.basis.monomials[K]
    J = setdiff(1:dimension(set), I)
    monos = DynamicPolynomials.MonomialVector(
        monos.vars[J],
        Vector{Int}[z[J] for z in monos.Z]
    )
    p = SumOfSquares.GramMatrix(Q, MB.MonomialBasis(monos))
    return ConvexPolySet(set.degree, p, nothing)
end

convexity_proof(set::ConvexPolySet) = set.convexity_proof

function scaling_function(set::Union{PolySet, ConvexPolySet})
    @assert dimension(set) == 2
    # We convert the GramMatrix to a polynomial to avoid having to do the
    # conversion for every substitution.
    p = MP.polynomial(set.p)
    vars = MP.variables(p)
    @assert length(vars) == 2
    vx, vy = vars
    return (x, y) -> p(vx => x, vy => y)^(1 / set.degree)
end

function _print_gauge_function(set::ConvexPolySet; digits=6)
    print(" ")
    p = MP.polynomial(set.p)
    if digits !== nothing
        p = round(p, digits=digits)
    end
    println(p)
end

"""
    struct ConvexPolynomialSet{T, U} <: AbstractSet{U}
        degree::Int
        q::GramMatrix{T, MonoBasis, u}
        z::SpaceVariable
        x::Vector{SpaceVariable}
    end

Set ``\\{\\, (z, x) \\mid p(z, x) \\le 0 \\,\\}`` or
``H \\{\\, (z, x) \\mid q(z, x) \\le z^{\\texttt{degree}} \\,\\}`` where `p` and
`q` are homogeneous polynomials of degree `degree` and `H` is a householder
matrix.
"""
struct ConvexPolynomialSet{T, B, U} <: AbstractSet{U}
    degree::Int
    q::GramMatrix{T, B, U}
    z::SpaceVariable
    x::Vector{SpaceVariable}
end
perspective_gauge0(set) = set.q - set.z^set.degree
perspective_variable(set::ConvexPolynomialSet) = set.z
space_variables(set::ConvexPolynomialSet) = set.x
function gauge1(set::ConvexPolynomialSet{T}) where T
    return subs(set.q, perspective_variable(set) => one(T))
end