text:
@testset "Shape" begin
    @testset "2D Box" begin
        rectangle = Box([[0.0,0.0], [2.0, 3.0]])
        o = origin(rectangle)
        @test volume(rectangle) ≈ 6.0
        @test name(rectangle) == "Rectangle"

        @test o ∈ rectangle
        @test (o+[0.0,3.1]) ∉ rectangle
        @test iscongruent(rectangle, congruent(rectangle, SVector(3.0,4.0)))
        @test rectangle ≅ congruent(rectangle, SVector(3.0,4.0)) # test infix

        smaller_rectangle = Box(o, [1.0, 2.0])
        @test smaller_rectangle ⊆ rectangle
        @test rectangle == bounding_box(rectangle)
        @test rectangle ⊆ Sphere(o, outer_radius(rectangle)) # a 2D sphere is a circle

        # Test different ways to construct Rectangle produce same results
        @test Box((0.0,0.0), (1.0, 2.0)) == Box([0.0,0.0], [1.0, 2.0]) # Tuple or Vector
        @test Box((0.0,0.0), (1.0, 2.0)) == Box([1.0, 2.0]) # Assumes origin at zero

        @testset "Boundary functions" begin
            x, y = boundary_functions(rectangle)
            @test x.(0:0.1:1) ≈ [0.0, 0.8, 1.6, 2.0, 2.0, 2.0, 1.2, 0.4, 0.0, 0.0, 0.0]
            @test y.(0:0.1:1) ≈ [0.0, 0.0, 0.0, 0.6, 1.8, 3.0, 3.0, 3.0, 2.4, 1.2, 0.0]
            @test_throws(DomainError,x(-eps(Float64)))
            @test_throws(DomainError,x(1.0 + eps(Float64)))
            @test_throws(DomainError,y(-eps(Float64)))
            @test_throws(DomainError,y(1.0 + eps(Float64)))
        end
    end

    @testset "Circle" begin
        radius = 2.0
        o = [6.7,8.9]
        circle = Sphere(o, radius)
        circle_bounding_box = bounding_box(circle)
        @test volume(circle) / volume(circle_bounding_box) ≈ 0.7853981633974483
        @test name(circle) == "Circle"

        @test outer_radius(circle) == radius
        @test o ∈ circle
        @test (o+[0.0,radius+0.1]) ∉ circle
        @test iscongruent(circle, congruent(circle,SVector(3.0,4.0)))
        @test circle ≅ congruent(circle,SVector(3.0,4.0))
        @test volume(circle) ≈ 12.566370614359172

        smaller_circle = Sphere(o, radius-2)
        @test smaller_circle ⊆ circle
        @test circle ⊆ circle_bounding_box

        # Test different ways to construct Circle produce same results
        @test Sphere((0.0,0.0), 1.0) == Sphere([0.0,0.0], 1.0) # Tuple or Vector
        @test Sphere((0.0,0.0), 1.0) == Sphere(2, 1.0) # Assumes origin at zero

        @testset "Boundary functions" begin
            x, y = boundary_functions(Sphere([-1.0,2.0],3.0))
            @test x.(0:0.1:1) ≈
            [2.0, 1.4270509831248424, -0.07294901687515765, -1.927050983124842, -3.427050983124842, -4.0, -3.4270509831248424, -1.9270509831248428, -0.07294901687515831, 1.427050983124842, 2.0]
            @test y.(0:0.1:1) ≈
            [2.0, 3.763355756877419, 4.853169548885461, 4.853169548885461, 3.7633557568774196, 2.0000000000000004, 0.2366442431225808, -0.8531695488854605, -0.8531695488854609, 0.23664424312257992, 1.9999999999999993]
            @test_throws(DomainError,x(-eps(Float64)))
            @test_throws(DomainError,x(1.0 + eps(Float64)))
            @test_throws(DomainError,y(-eps(Float64)))
            @test_throws(DomainError,y(1.0 + eps(Float64)))
        end
    end

    @testset "Time of flight" begin
        time_of_flight = TimeOfFlightPlaneWaveToPoint([-10.0,0.0], 40.0; minimum_x = 0.0)
        time_of_flight_bounding_box = bounding_box(time_of_flight)
        ratio = volume(time_of_flight) / volume(time_of_flight_bounding_box)
        # Geometric arguments dictate that the ratio must be between 0.5 and 1.0
        @test ratio > 0.5
        @test ratio < 1.0
        @test name(time_of_flight) == "Time of flight from planar source to the focal point"
    end

    @testset "Time of flight from point" begin
        time_of_flight = TimeOfFlightPointWaveToPoint([-10.0,0.0],40.0)
        time_of_flight_bounding_box = bounding_box(time_of_flight)
        ratio = volume(time_of_flight) / volume(time_of_flight_bounding_box)
        # Geometric arguments dictate that the ratio must be between 0.5 and 1.0
        @test ratio > 0.5
        @test ratio < 1.0
        @test name(time_of_flight) == "Time of flight from point source"

        @testset "Boundary functions" begin
            x, y = boundary_functions(TimeOfFlightPointWaveToPoint([-1.0,0.0],3.0))
            @test x.(0:0.1:1) ≈ [0.0,1.2182846930812632,1.9095426112571139,1.9095426112571139,1.2182846930812632,0.0,0.0,0.0,0.0,0.0,0.0]
            @test y.(0:0.1:1) ≈ [-2.82842712474619,-2.019706171808505,-0.7311373286046446,0.7311373286046446,2.0197061718085054,2.82842712474619,1.6970562748477145,0.5656854249492386,-0.5656854249492386,-1.6970562748477145,-2.8284271247461903]
            @test_throws(DomainError,x(-eps(Float64)))
            @test_throws(DomainError,x(1.0 + eps(Float64)))
            @test_throws(DomainError,y(-eps(Float64)))
            @test_throws(DomainError,y(1.0 + eps(Float64)))
        end
    end

    @testset "Sphere" begin
        radius = 4.2
        o = SVector(3.4, -2.1, 2.0)
        sphere = Sphere(o, radius)
        @test name(sphere) == "Sphere"
        @test outer_radius(sphere) == radius
        @test o ∈ sphere
        @test (o+[0.0,radius+0.1,0.0]) ∉ sphere
        @test iscongruent(sphere, congruent(sphere,SVector(3.0,4.0,1.0)))
        @test volume(sphere) ≈ 310.33908869221415

        smaller_sphere = Sphere(o, radius-2)
        @test smaller_sphere ⊆ sphere

        # Test different ways to construct Rectangle produce same results
        @test Sphere((0.0,0.0,0.0), 1.0) == Sphere([0.0,0.0,0.0], 1.0)
    end

    # Causing Segmentation fault on Julia 1.1
    @testset "Plot Shapes" begin
        # Just try each to see if we have any errors (yes thats a very low bar)
        rectangle = Box([[0.0,0.0],[2.0,3.0]])
        plot(rectangle)

        circle = Sphere([-1.0,2.0],2.0)
        plot!(circle)

        # timeofflight = TimeOfFlightPlaneWaveToPoint([-1.0,0.0],3.0)
        # plot!(timeofflight)

        timeofflightfrompoint = TimeOfFlightPointWaveToPoint([-1.0,0.0],3.0)
        plot!(timeofflightfrompoint)

        @test true
    end

    @testset "Infinite volume shapes" begin
        normal = rand(3);
        origin = rand(3);
        width = 1.0 + rand();
        p = Plate(normal, width, origin)
        h = Halfspace(normal, origin)

        shapes = [p,h];

        @test isempty(findall( (volume.(shapes) .== Inf) .== 0))

        xs = [rand(3) for i = 1:10];
        [[x ∈ s for x in xs] for s in shapes]

    end
end


text:
function stress_cell_mass_matrix(basis, quad, detjac)
    dim = dimension(basis)
    sdim = number_of_symmetric_degrees_of_freedom(dim)
    return mass_matrix(basis, quad, detjac, sdim)
end

function stress_cell_rhs_operator(basis, quad, stiffness, jacobian)
    dim = dimension(basis)
    sdim = number_of_symmetric_degrees_of_freedom(dim)
    nf = number_of_basis_functions(basis)
    op = zeros(sdim * nf, dim * nf)
    vectosymmconverter = vector_to_symmetric_matrix_converter()
    detjac = prod(jacobian)

    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jacobian)
        NK = sum([make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim])
        vals = basis(p)
        NI = interpolation_matrix(vals, sdim)

        op .+= NI' * stiffness * NK * detjac * w
    end
    return op
end

function stress_cell_rhs(basis, quad, stiffness, celldisp, jacobian)

    op = stress_cell_rhs_operator(basis, quad, stiffness, jacobian)
    return op * celldisp
end

ext: jl


text:
mutable struct Record
    id::Integer
    recordtype::String
    name::String
    data::String
    priority::Union{Nothing, Integer}
    port::Union{Nothing, Integer}
    ttl::Integer
    weight::Union{Nothing, Integer}
    flags::Union{Nothing, Integer}
    tag::Union{Nothing, String}

    function Record(data::Dict{String})
        new(
            data["id"],
            data["type"],
            data["name"],
            data["data"],
            data["priority"],
            data["port"],
            data["ttl"],
            data["weight"],
            data["flags"],
            data["tag"],
        )
    end
end

function show(io::IO, r::Record)
    print(io, "Record ($(r.recordtype), $(r.name))")
end

# List all domain records
function getalldomainrecords!(client::AbstractClient, domainname::String)
    uri = joinpath(ENDPOINT, "domains", domainname, "records")
    getalldata!(client, uri, Record)
end

function getalldomainrecords!(client::AbstractClient, domain::Domain)
    getalldomainrecords!(client, domain.name)
end

# Retrieve an existing domain record
function getdomainrecord!(client::AbstractClient, domainname::String, record_id::Integer)
    uri = joinpath(ENDPOINT, "domains", domainname, "records", "$record_id")
    Record(getdata!(client, uri))
end

function getdomainrecord!(client::AbstractClient, domain::Domain, record_id::Integer)
    getdomainrecord!(client, domain.name, record_id)
end

# Create a new domain record
function createdomainrecord!(client::AbstractClient, domainname::String,
                             recordtype::String; kwargs...)
    postbody = Dict{String, Any}("type" => recordtype)
    merge!(postbody, Dict{String, Any}([String(k[1]) => k[2] for k in kwargs]))

    if recordtype in ("A", "AAAA", "CAA", "CNAME", "MX", "TXT", "SRV", "NS")
        if !(recordtype in ("MX", "NS") || haskey(postbody, "name"))
            error("'name' is a required argument for type '$(recordtype)'")
        end
        if !haskey(postbody, "data")
            error("'data' is a required argument for type '$(recordtype)'")
        end
    end
    if recordtype in ("SRV", "MX") && !haskey(postbody, "priority")
        error("'priority' is a required argument for type '$(recordtype)'")
    end
    if recordtype == "SRV"
        if !haskey("postbody", "port")
            error("'port' is a required argument for type '$(recordtype)'")
        end
        if !haskey(postbody, "weight")
            error("'weight' is a required argument for type '$(recordtype)'")
        end
    end
    if recordtype == "CAA"
        if !haskey(postbody, "flags")
            error("'flags' is a required argument for type '$(recordtype)'")
        end
        if !haskey(postbody, "tag")
            error("'tag' is a required argument for type '$(recordtype)'")
        end
    end

    uri = joinpath(ENDPOINT, "domains", domainname, "records")
    Record(postdata!(client, uri, postbody))
end

function createdomainrecord!(client::AbstractClient, domain::Domain,
                             recordtype::String; kwargs...)
    createdomainrecord!(client, domain.name, recordtype; kwargs...)
end

# Update a Domain Record
function updatedomainrecord!(client::AbstractClient, domainname::String,
                             record_id::Integer; kwargs...)
    putbody = Dict{String, Any}([String(k[1]) => k[2] for k in kwargs])

    uri = joinpath(ENDPOINT, "domains", domainname, "records", "$(record_id)")
    Record(putdata!(client, uri, putbody))
end

function updatedomainrecord!(client::AbstractClient, domain::Domain, record_id::Integer; kwargs...)
    updatedomainrecord!(client, domain.name, record_id; kwargs...)
end

function updatedomainrecord!(client::AbstractClient, domainname::String, record::Record; kwargs...)
    data = updatedomainrecord!(client, domainname, record.id; kwargs...)
    for k in kwargs
        if k[1] in fieldnames(Domain)
            setfield!(record, Symbol(k[1]), k[2])
        end
    end
    data
end

function updatedomainrecord!(client::AbstractClient, domain::Domain, record::Record; kwargs...)
    updatedomainrecord!(client, domain.name, record; kwargs...)
end

# Delete a Domain Record
function deletedomainrecord!(client::AbstractClient, domainname::String, record_id::Integer)
    uri = joinpath(ENDPOINT, "domains", domainname, "records", "$(record_id)")
    deletedata!(client, uri)
end

function deletedomainrecord!(client::AbstractClient, domainname::String, record::Record)
    deletedomainrecord!(client, domainname, record.id)
end


text:
using Revise
	using PseudoArcLengthContinuation, LinearAlgebra, Plots, SparseArrays, Setfield, Parameters
	const PALC = PseudoArcLengthContinuation

f1(u, v) = u * u * v
norminf = x -> norm(x, Inf)

function plotsol(x; kwargs...)
	N = div(length(x), 2)
	plot!(x[1:N], label="u"; kwargs...)
	plot!(x[N+1:2N], label="v"; kwargs...)
end

function Fbru!(f, x, p)
	@unpack α, β, D1, D2, l = p
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h
	c1 = D1 / l^2 / h2
	c2 = D2 / l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	# Dirichlet boundary conditions
	f[1]   = c1 * (α	  - 2u[1] + u[2] ) + α - (β + 1) * u[1] + f1(u[1], v[1])
	f[end] = c2 * (v[n-1] - 2v[n] + β / α)			 + β * u[n] - f1(u[n], v[n])

	f[n]   = c1 * (u[n-1] - 2u[n] +  α   ) + α - (β + 1) * u[n] + f1(u[n], v[n])
	f[n+1] = c2 * (β / α  - 2v[1] + v[2])			 + β * u[1] - f1(u[1], v[1])

	for i=2:n-1
		  f[i] = c1 * (u[i-1] - 2u[i] + u[i+1]) + α - (β + 1) * u[i] + f1(u[i], v[i])
		f[n+i] = c2 * (v[i-1] - 2v[i] + v[i+1])			  + β * u[i] - f1(u[i], v[i])
	end
	return f
end

function Fbru(x, p)
	f = similar(x)
	Fbru!(f, x, p)
	f
end

function Jbru_sp(x, p)
	@unpack α, β, D1, D2, l = p
	# compute the Jacobian using a sparse representation
	n = div(length(x), 2)
	h = 1.0 / n; h2 = h*h

	c1 = D1 / p.l^2 / h2
	c2 = D2 / p.l^2 / h2

	u = @view x[1:n]
	v = @view x[n+1:2n]

	diag   = zeros(eltype(x), 2n)
	diagp1 = zeros(eltype(x), 2n-1)
	diagm1 = zeros(eltype(x), 2n-1)

	diagpn = zeros(eltype(x), n)
	diagmn = zeros(eltype(x), n)

	@. diagmn = β - 2 * u * v
	@. diagm1[1:n-1] = c1
	@. diagm1[n+1:end] = c2

	@. diag[1:n]    = -2c1 - (β + 1) + 2 * u * v
	@. diag[n+1:2n] = -2c2 - u * u

	@. diagp1[1:n-1]   = c1
	@. diagp1[n+1:end] = c2

	@. diagpn = u * u
	J = spdiagm(0 => diag, 1 => diagp1, -1 => diagm1, n => diagpn, -n => diagmn)
	return J
end

function finalise_solution(z, tau, step, contResult)
	n = div(length(z.u), 2)
	printstyled(color=:red, "--> Solution constant = ", norm(diff(z.u[1:n])), " - ", norm(diff(z.u[n+1:2n])), "\n")
	return true
end

n = 500
####################################################################################################
# test for the Jacobian expression
# using ForwardDiff
# sol0 = rand(2n)
# J0 = ForwardDiff.jacobian(x-> Fbru(x, par_bru), sol0) |> sparse
# J1 = Jbru_sp(sol0, par_bru)
# J0 - J1
####################################################################################################
# different parameters to define the Brusselator model and guess for the stationary solution
par_bru = (α = 2., β = 5.45, D1 = 0.008, D2 = 0.004, l = 0.3)
	sol0 = vcat(par_bru.α * ones(n), par_bru.β/par_bru.α * ones(n))

# # parameters for an isola of stationary solutions
# par_bru = (α = 2., β = 4.6, D1 = 0.0016, D2 = 0.008, l = 0.061)
# 	xspace = LinRange(0, par_bru.l, n)
# 	sol0 = vcat(		par_bru.α .+ 2 .* sin.(pi*xspace/par_bru.l),
# 			par_bru.β/par_bru.α .- 0.5 .* sin.(pi*xspace/par_bru.l))

# eigls = EigArpack(1.1, :LM)
# 	opt_newton = PALC.NewtonPar(eigsolver = eigls)
# 	out, hist, flag = @time PALC.newton(
# 		x ->    Fbru(x, par_bru),
# 		x -> Jbru_sp(x, par_bru),
# 		sol0, opt_newton, normN = norminf)
#
# 	plot();plotsol(out);plotsol(sol0, label = "sol0",line=:dash)
####################################################################################################
eigls = EigArpack(1.1, :LM)
opts_br_eq = ContinuationPar(dsmin = 0.001, dsmax = 0.01, ds = 0.001, pMax = 1.9, detectBifurcation = 2, nev = 21, plotEveryNsteps = 50, newtonOptions = NewtonPar(eigsolver = eigls, tol = 1e-9), maxSteps = 1060)

	br, _ = @time continuation(
		(x, p) ->    Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		sol0, par_bru.l,
		opts_br_eq, verbosity = 0,
		plot = true,
		plotSolution = (x; kwargs...) -> (plotsol(x; label="", kwargs... )),
		printSolution = (x, p) -> x[div(n,2)], normC = norminf)
#################################################################################################### Continuation of the Hopf Point using Jacobian expression
ind_hopf = 2
	# hopfpt = PALC.HopfPoint(br, ind_hopf)
	optnew = opts_br_eq.newtonOptions
	hopfpoint, _, flag = @time newtonHopf(
		(x, p) ->    Fbru(x, @set par_bru.l = p),
		(x, p) -> Jbru_sp(x, @set par_bru.l = p),
		br, ind_hopf,
		(@set optnew.verbose=true), normN = norminf)
	flag && printstyled(color=:red, "--> We found a Hopf Point at l = ", hopfpoint.p[1], ", ω = ", hopfpoint.p[2], ", from l = ", br.bifpoint[ind_hopf].param, "\n")

# br_hopf, u1_hopf = @time PALC.continuationHopf(
# 	(x, p, β) ->  Fbru(x, setproperties(par_bru, (l=p, β=β))),
# 	(x, p, β) -> Jbru_sp(x, setproperties(par_bru, (l=p, β=β))),
# 	br, ind_hopf, par_bru.β,
# 	ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, newtonOptions = opt_newton), verbosity = 2, normC = norminf)
# PALC.plotBranch(br_hopf, xlabel="beta", ylabel = "l", label="")
#
# # test with analytical Hessian but with dummy expression ;)
# d2Fbru(x, p, dx1, dx2) = dx1 .* dx2
#
# hopfpoint, hist, flag = @time PALC.newtonHopf(
# 	(x, p) ->  Fbru(x, @set par_bru.l = p),
# 	(x, p) -> Jbru_sp(x, @set par_bru.l = p),
# 	(x, p) -> transpose(Jbru_sp(x, @set par_bru.l = p)),
# 	d2Fbru,
# 	br, ind_hopf,
# 	(@set opt_newton.verbose = true), normN = norminf)
#
# br_hopf, u1_hopf = @time PALC.continuationHopf(
# 	(x, p, β) ->  Fbru(x, setproperties(par_bru, (l=p, β=β))),
# 	(x, p, β) -> Jbru_sp(x, setproperties(par_bru, (l=p, β=β))),
# 	(x, p, β) -> transpose(Jbru_sp(x, setproperties(par_bru, (l=p, β=β)))),
# 	β -> d2Fbru,
# 	br, ind_hopf, par_bru.β,
# 	ContinuationPar(dsmin = 0.001, dsmax = 0.05, ds= 0.01, pMax = 6.5, pMin = 0.0, newtonOptions = opt_newton), verbosity = 2, normC = norminf)

####################################################################################################Continuation of Periodic Orbit
# number of time slices
M = 51
l_hopf, Th, orbitguess2, hopfpt, vec_hopf = PALC.guessFromHopf(br, ind_hopf, opts_br_eq.newtonOptions.eigsolver, M, 2.7; phase = 0.25)
#
# orbitguess_f2 = orbitguess2[1]; for ii=2:M; global orbitguess_f2 = hcat(orbitguess_f2, orbitguess2[ii]);end
orbitguess_f2 = reduce(vcat, orbitguess2)
orbitguess_f = vcat(vec(orbitguess_f2), Th) |> vec

poTrap = p -> PeriodicOrbitTrapProblem(
			x ->	Fbru(x, @set par_bru.l = p),
			x -> Jbru_sp(x, @set par_bru.l = p),
			real.(vec_hopf),
			hopfpt.u,
			M)

poTrap(l_hopf + 0.01)(orbitguess_f) |> plot
PALC.plotPeriodicPOTrap(orbitguess_f, n, M; ratio = 2)

using ForwardDiff
d1Fbru(x, p, dx) = ForwardDiff.derivative(t -> Fbru(x .+ t .* dx, p), 0.)

ls0 = GMRESIterativeSolvers(N = 2n, tol = 1e-9)#, Pl = lu(I + par_cgl.Δ))
poTrapMF = p -> PeriodicOrbitTrapProblem(
			x ->	Fbru(x, @set par_bru.l = p),
			x -> (dx -> d1Fbru(x, (@set par_bru.l = p), dx)),
			real.(vec_hopf),
			hopfpt.u,
			M, ls0)

deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [zero(orbitguess_f)])
# deflationOp = DeflationOperator(2.0, (x,y) -> dot(x[1:end-1], y[1:end-1]),1.0, [outpo_f])
####################################################################################################
opt_po = PALC.NewtonPar(tol = 1e-9, verbose = true, maxIter = 20)
	outpo_f, _, flag = @time PALC.newton(poTrap(l_hopf + 0.01),
			orbitguess_f,
			# ig,
			# (br_po.sol[11].u),
			(@set opt_po.maxIter = 250),
			# opt_po,
			# deflationOp,
			# :FullLU;
			normN = norminf,
			callback = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", PALC.amplitude(x, n, M; ratio = 2));true)
			)
	flag && printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", PALC.amplitude(outpo_f, n, M; ratio = 2),"\n")
	PALC.plotPeriodicPOTrap(outpo_f, n, M; ratio = 2)

PALC.plotPeriodicPOTrap(orbitguess_f, n, M; ratio = 2)
opt_po = @set opt_po.eigsolver = EigKrylovKit(tol = 1e-5, x₀ = rand(2n), verbose = 2)
opt_po = @set opt_po.eigsolver = DefaultEig()
# opt_po = @set opt_po.eigsolver = EigArpack(; tol = 1e-5, v0 = rand(2n))
opts_po_cont = ContinuationPar(dsmin = 0.001, dsmax = 0.03, ds= 0.01, pMax = 3.0, maxSteps = 100, newtonOptions = opt_po, nev = 5, precisionStability = 1e-6, detectBifurcation = 1)
	br_po, _ , _= @time PALC.continuationPOTrap(poTrap,
			outpo_f, l_hopf + 0.01,
			opts_po_cont;
			verbosity = 2,	plot = true,
			# callbackN = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", PALC.amplitude(x, n, M));true),
			plotSolution = (x;kwargs...) -> heatmap!(reshape(x[1:end-1], 2*n, M)'; ylabel="time", color=:viridis, kwargs...),
			normC = norminf)

####################################################################################################
using IncompleteLU
Jpo = poTrap(l_hopf + 0.01)(Val(:JacFullSparse), orbitguess_f)
Jpo = poTrap(l_hopf + 0.01)(Val(:JacCyclicSparse), orbitguess_f)

Precilu = @time ilu(Jpo, τ = 0.005)
ls = GMRESIterativeSolvers(verbose = false, tol = 1e-4, N = size(Jpo,1), restart = 30, maxiter = 150, log=true, Pl = Precilu)
	@time ls(Jpo, rand(ls.N))

opt_po = PALC.NewtonPar(tol = 1e-10, verbose = true, maxIter = 20)
	outpo_f, _, flag = @time PALC.newton(poTrap(l_hopf + 0.01),
			orbitguess_f,
			(@set opt_po.linsolver = ls), :BorderedMatrixFree;
			normN = norminf,
			# callback = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", amplitude(x), " T = ", x[end], ", T0 = ",orbitguess_f[end]);true)
			)
	printstyled(color=:red, "--> T = ", outpo_f[end], ", amplitude = ", PALC.amplitude(outpo_f, n, M; ratio = 2),"\n")
	PALC.plotPeriodicPOTrap(outpo_f, n, M; ratio = 2)

opts_po_cont = ContinuationPar(dsmin = 0.0001, dsmax = 0.05, ds= 0.01, pMax = 2.2, maxSteps = 3000, newtonOptions = (@set opt_po.linsolver = ls))
	br_pok, _ , _= @time PALC.continuationPOTrap(poTrap,
			outpo_f, l_hopf + 0.01,
			opts_po_cont , :FullMatrixFree;
			verbosity = 2,
			plot = true,
			# callbackN = (x, f, J, res, iteration, options; kwargs...) -> (println("--> amplitude = ", PALC.amplitude(x, n, M));true),
			plotSolution = (x;kwargs...) -> heatmap!(reshape(x[1:end-1], 2*n, M)'; ylabel="time", color=:viridis, kwargs...), normC = norminf)


text:
# # Fast-Diffusion Problem
#
# This example demonstrates the use of 'NonLinearDiffusion' operator to solve time-dependent non-linear diffusion PDEs with coefficient having dependence on the unknown function.
# Here we consider a fast diffusion problem with Dirichlet BCs on unit interval:
# ∂ₜu = ∂ₓ(k*∂ₓu)
# k = 1/u²
# u(x=0,t) = exp(-t)
# u(x=1,t) = 1/(1.0 + exp(2t))
# u(x, t=0) = u₀(x)
using Test
using DiffEqOperators, OrdinaryDiffEq

@testset "1D Nonlinear fast diffusion equation" begin
    # The analytical solution for this is given by :

    u_analytic(x, t) = 1 / sqrt(x^2 + exp(2*t))

    #
    # Reproducing it numerically
    #

    nknots = 100
    h = 1.0/(nknots+1)
    knots = range(h, step=h, length=nknots)
    n = 1                                   # Outer differential order
    m = 1                                   # Inner differential order
    approx_ord = 2

    u0 = u_analytic.(knots,0.0)
    du = similar(u0)
    t0 = 0.0
    t1 = 1.0

    function f(du,u,p,t)
        bc = DirichletBC(exp(-t),(1.0 + exp(2*t))^(-0.5))
        l = bc*u
        k = l.^(-2)                        # Diffusion Coefficient
        nonlinear_diffusion!(du,n,m,approx_ord,k,l,h,nknots)
    end

    prob = ODEProblem(f, u0, (t0, t1))
    alg = KenCarp4()
    sol = solve(prob,alg)

    for t in t0:0.1:t1
        @test u_analytic.(knots, t) ≈ sol(t) rtol=1e-3
    end

end


text:
using RchannelExperiments
using Base.Test
const R=RchannelExperiments
const RI=RchannelExperiments.RcImages
# write your own tests here

x = linspace(-5pi,5pi,156);
vec = sin.(x).*cos.(x/5)+0.005x.^2;
peaks, ips, proms, doms = RI.findpeaks(vec);
@test peaks ≈ [1.96984, 0.497992, 0.709033, 0.964469, 0.362738, 1.24738] atol=1e-5
@test maximum(peaks)==maximum(vec)
@test all(doms.==[148, 17, 27, 45, 10, 121])
@test proms ≈ [2.9112, 0.369982, 0.694124, 1.47311, 0.10339, 2.18874] atol=1e-5

x = linspace(-5pi,5pi,156);
vec = sin.(x).*cos.(x/5)+0.005x.^2;
peaks, ips, proms, doms = RI.findpeaks(vec);
peaks_l, ips_l, proms_l, doms_l = RI.findpeaks_onesided(vec, :left);
peaks_r, ips_r, proms_r, doms_r = RI.findpeaks_onesided(vec, :right);
@test peaks ≈ peaks_r
@test peaks ≈ peaks_l
rs_d = [1,2,3,4,5]
ls_d = [6]
rs_p = [1,2,4,5]
ls_p = [3,6]
@test proms[rs_p] ≈ proms_r[rs_p]
@test proms[ls_p] ≈ proms_l[ls_p]
@test doms[rs_d] ≈ doms_r[rs_d]
@test doms[ls_d] ≈ doms_l[ls_d]


text:
# This file contains the ingredients to create a PackageManager for BinDeps
import BinDeps: Binaries, can_use, package_available, bindir, libdir,
                generate_steps, LibraryDependency, provider, provides
import Base: show

mutable struct BP <: Binaries
    url::String
    hash::String
    prefix::Prefix
end

show(io::IO, p::BP) = write(io, "BinaryProvider for $(p.url)")

# We are cross-platform baby, and we never say no to a party
can_use(::Type{BP}) = true
package_available(p::BP) = true
libdir(p::BP, dep) = @static if Compat.Sys.iswindows()
    joinpath(p.prefix, "bin")
else
    joinpath(p.prefix, "lib")
end

# We provide our own overload of provides() for BP
macro BP_provides(url, hash, dep, opts...)
    return quote
        prefix = Prefix(joinpath(dirname(@__FILE__), "usr"))
        activate(prefix)
        return provides(BP, ($url, $hash, prefix), $(esc(dep)), $(opts...))
    end
end
provider(::Type{BP}, data; opts...) = BP(data...)

function generate_steps(dep::LibraryDependency, p::BP, opts)
    () -> begin
        install(p.url, p.hash; prefix=p.prefix, verbose=true)
    end
end


text:
nq = 7
nv = 6
nx = nq + nv

r = 0.5
c1 = @SVector [r, r, r]
c2 = @SVector [r, r, -r]
c3 = @SVector [r, -r, r]
c4 = @SVector [r, -r, -r]
c5 = @SVector [-r, r, r]
c6 = @SVector [-r, r, -r]
c7 = @SVector [-r, -r, r]
c8 = @SVector [-r, -r, -r]

corner_offset = [c1, c2, c3, c4, c5, c6, c7, c8]

function attitude_jacobian(q)
	s = q[1]
	v = q[2:4]

	[-transpose(v);
	 s * I + skew(v)]
end

function G_func(q)
	return [I zeros(eltype(q), 3, 3);
		zeros(eltype(q), 4, 3) attitude_jacobian(q)]
end

function quaternion_rotation_matrix(q)
	r, i, j, k  = q

	r11 = 1.0 - 2.0 * (j^2.0 + k^2.0)
	r12 = 2.0 * (i * j - k * r)
	r13 = 2.0 * (i * k + j * r)

	r21 = 2.0 * (i * j + k * r)
	r22 = 1.0 - 2.0 * (i^2.0 + k^2.0)
	r23 = 2.0 * (j * k - i * r)

	r31 = 2.0 * (i * k - j * r)
	r32 = 2.0 * (j * k + i * r)
	r33 = 1.0 - 2.0 * (i^2.0 + j^2.0)

	SMatrix{3,3}([r11 r12 r13;
	              r21 r22 r23;
				  r31 r32 r33])
end

function kinematics(q)
    p = q[1:3]
    quat = q[4:7]

    R = quaternion_rotation_matrix(quat)

	p1 = p + R * c1
	p2 = p + R * c2
	p3 = p + R * c3
	p4 = p + R * c4
	p5 = p + R * c5
	p6 = p + R * c6
	p7 = p + R * c7
	p8 = p + R * c8

    SVector{24}([p1;
	              p2;
				  p3;
				  p4;
				  p5;
				  p6;
				  p7;
				  p8])
end

ϕ(x1)

function ϕ(q)
	k = kinematics(q)
	idx = collect([3, 6, 9, 12, 15, 18, 21, 24])

	return k[idx]
end

P(q) = ForwardDiff.jacobian(kinematics, q) * G_func(q)
N(q) = ForwardDiff.jacobian(ϕ, q) * G_func(q)

mass = 1.0
inertia = 1.0
gravity = 9.81
mass_matrix = Diagonal(mass * ones(3))
inertia_matrix = Diagonal(inertia * ones(3))
mi_matrix = cat(mass_matrix, inertia_matrix, dims=(1,2))
function bias_vector(q, v)
	ω = v[4:6]
	[0.0; 0.0; gravity * mass; cross(ω, inertia_matrix * ω)]
end

h = 0.01

softplus(x) = log.(1.0 .+ exp.(x))
softminus(x) = x .- softplus.(x)

function L_multiply(q)
	s = q[1]
	v = q[2:4]

	SMatrix{4,4}([s -transpose(v);
	              v s * I + skew(v)])
end

softminus(1.0)
N(ones(nq))
function dynamics(x, u, t)
	xq = x[1:nq]
	pos = xq[1:3]
	quat = xq[4:7]

	xv = x[nq .+ (1:nv)]
	vel = xv[1:3]
	omg = xv[4:6]

	τ = [0.0; 0.0; 0.0; 0.0; 0.0; 0.0]

	# k_impact = 10.0
	# λ = [0.0; -k_impact * softminus(ϕ(q))]
	k_impact = 2.35
	b_impact = 1.0
	λ = zeros(24)
	impact_idx = collect([3, 6, 9, 12, 15, 18, 21, 24])
	@show λ[impact_idx] = -k_impact * softminus.(ϕ(xq)) + max.(0.0, - b_impact .* N(xq) * xv)
	# λ = [0.0; -k_impact * min(0.0, ϕ(q))]
	# λ = [0.0; -k_impact * softminus(ϕ(q)) + max(0.0, - b_impact * v[2])]

	xv⁺ = xv + h * (mi_matrix \ (τ - bias_vector(xq, xv) + transpose(P(xq)) * λ))
	pos⁺ = pos + h * xv⁺[1:3]
	quat⁺ = 0.5 * h * L_multiply(quat) * [sqrt((0.5 / h)^2.0 - xv⁺[4:6]' * xv⁺[4:6]); xv⁺[4:6]]

	return [pos⁺; quat⁺; xv⁺]
end
sqrt(2.0)
pos1 = [0.0; 0.0; 1.0]
quat1 = [0.5 * sqrt(2.0); 0.5 * sqrt(2.0); 0.0; 0.0]
vel1 = [0.0; 0.0; 0.0]
omg1 = [0.0; 0.0; 0.0]

x1 = [pos1; quat1; vel1; omg1]

x_hist = [x1]
T = 500

for t = 1:T-1
	push!(x_hist, dynamics(x_hist[end], nothing, t))
	# println(ϕ(x_hist[end][1:nq]))
end
visualize!(vis, x_hist, Δt = h)

x_hist[end]
# plot(hcat(x_hist...)[1:2, :]', )

# include(joinpath(pwd(), "models/visualize.jl"))
# vis = Visualizer()
# render(vis)
# visualize!(vis, x_hist, Δt = h)
#
# function visualize!(vis, q;
#         Δt = 0.1)
#
# 	default_background!(vis)
#
#     setobject!(vis["box"], GeometryBasics.Rect(Vec(-1.0 * r,
# 		-1.0 * r,
# 		-1.0 * r),
# 		Vec(2.0 * r, 2.0 * r, 2.0 * r)),
# 		MeshPhongMaterial(color = RGBA(0.0, 0.0, 0.0, 1.0)))
#
#     for i = 1:8
#         setobject!(vis["corner$i"], GeometryBasics.Sphere(Point3f0(0),
#             convert(Float32, 0.05)),
#             MeshPhongMaterial(color = RGBA(1.0, 165.0 / 255.0, 0.0, 1.0)))
#     end
#
#     anim = MeshCat.Animation(convert(Int, floor(1.0 / Δt)))
#
#     for t = 1:length(q)
#         MeshCat.atframe(anim, t) do
#             settransform!(vis["box"],
# 				compose(Translation(q[t][1:3]...), LinearMap(UnitQuaternion(q[t][4:7]...))))
#
#             for i = 1:8
#                 settransform!(vis["corner$i"],
#                     Translation((q[t][1:3] + UnitQuaternion(q[t][4:7]...) * (corner_offset[i]))...))
#             end
#         end
#     end
#     # settransform!(vis["/Cameras/default"], compose(Translation(-1, -1, 0),LinearMap(RotZ(pi/2))))
#     MeshCat.setanimation!(vis, anim)
# end

max_issues_repo_name: thowell/DirectMotionPlanning


text:
import ZeroRPC

module TestServer

function hello(name)
  ret = @sprintf "Hello, %s" name
  return ret
end

function add(a, b)
  return a+b
end

export hello, add

end # module

server = ZeroRPC.Server(TestServer)
ZeroRPC.bind(server, "tcp://127.0.0.1:4242")
ZeroRPC.run(server)


text:
using ApproxFunBase, Test
import ApproxFunBase: PointSpace, HeavisideSpace, PiecewiseSegment, dimension, Vec, checkpoints

@testset "Spaces" begin
    @testset "PointSpace" begin
        @test eltype(domain(PointSpace([0,0.1,1])) ) == Float64

        f=Fun(x->(x-0.1),PointSpace([0,0.1,1]))
        @test roots(f) == [0.1]

        a=Fun(exp,space(f))
        @test f/a == Fun(x->(x-0.1)*exp(-x),space(f))

        f = Fun(space(f),[1.,2.,3.])
    end

    @testset "Derivative operator for HeavisideSpace" begin
        H = HeavisideSpace([-1.0,0.0,1.0])
        @test Fun(H, [1.0])(1.0) == 0.0
        @test Fun(H, [0.0,1.0])(1.0) == 1.0

        H=HeavisideSpace([1,2,3])
        D=Derivative(H)
        @test domain(D)==PiecewiseSegment([1,2,3])
        @test D[1,1]==-1
        @test D[1,2]==1

        H=HeavisideSpace([1,2,3,Inf])
        D=Derivative(H)
        @test domain(D)==PiecewiseSegment([1,2,3,Inf])
        @test D[1,1]==-1
        @test D[2,2]==-1
        @test D[1,2]==1

        S = HeavisideSpace([-1.0,0.0,1.0])
        @test Derivative(S) === Derivative(S,1)

        a = HeavisideSpace(0:0.25:1)
        @test dimension(a) == 4
        @test @inferred(points(a)) == 0.125:0.25:0.875
    end

    @testset "DiracDelta integration and differentiation" begin
        δ = DiracDelta()
        h = integrate(δ)
        @test domain(h) == PiecewiseSegment([0,Inf])
        @test h(-2) == 0
        @test h(2) == 1

        δ = 0.3DiracDelta(0.1) + 3DiracDelta(2.3)
        h = integrate(δ)
        @test domain(h) == PiecewiseSegment([0.1,2.3,Inf])
        @test h(-2) == 0
        @test h(2) == 0.3
        @test h(3) == 3.3

        δ = (0.3+1im)DiracDelta(0.1) + 3DiracDelta(2.3)
        h = integrate(δ)
        @test domain(h) == PiecewiseSegment([0.1,2.3,Inf])
        @test h(-2) == 0
        @test h(2) == 0.3+1im
        @test h(3) == 3.3+1im
    end

    @testset "Multivariate" begin
        a = HeavisideSpace(0:0.25:1)
        @test @inferred(dimension(a^2)) == dimension(a)^2
        @test @inferred(domain(a^2)) == domain(a)^2
        @test @inferred(points(a^2)) == vec(Vec.(points(a), points(a)'))
        @test  @inferred(checkpoints(a^2)) == vec(Vec.(checkpoints(a)', checkpoints(a)))

        aa2 = TensorSpace(a , a^2)
        @test dimension(aa2) == dimension(a)^3
        @test @inferred(domain(aa2)) == domain(a)^3
        @test @inferred(points(aa2)) == vec(Vec.(points(a), points(a)', reshape(points(a), 1,1,4)))
        @test  @inferred(checkpoints(aa2)) == vec(Vec.(reshape(checkpoints(a), 1,1,length(checkpoints(a))), checkpoints(a)', checkpoints(a)))

        @test dimension(a^3) == dimension(a)^3
        @test @inferred(domain(a^3)) == domain(a)^3
        @test_broken @inferred(points(a^3)) == vec(Vec.(points(a), points(a)', reshape(points(a), 1,1,4)))
    end
end


text:
export TracyWidom

"""
Tracy-Widom distribution

The probability distribution of the normalized largest eigenvalue of a random
Hermitian matrix.

The cdf of Tracy-Widom is given by

``
F_2 (s) = lim_{n→∞} Pr(√2 n^{1/6} (λₙ - √(2n) ≤ s)
``

References:

1. doi:10.1016/0370-2693(93)91114-3
2. doi:10.1007/BF02100489

Numerical routines adapted from Alan Edelman's course notes for MIT 18.338,
Random Matrix Theory, 2016.
"""
immutable TracyWidom <: ContinuousUnivariateDistribution end

"""
Probability density function of the Tracy-Widom distribution

Computes the Tracy-Widom distribution by directly solving the
Painlevé II equation using the ode23 numerical integrator

# Arguments
* `d::TracyWidom` or `Type{TracyWidom}`: an instance of `TracyWidom` or the type itself
* `t::Real`: The point at which to evaluate the pdf
* `t0::Real = -8.0`: The point at which to start integrating
"""
function pdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -8.0))
    t≤t0 && return 0.0
    t≥5  && return 0.0

    ts, y = _solve_painleve_ii(t0, t)

    ΔF2=exp(-y[end-1][3]) - exp(-y[end][3]) # the cumulative distribution
    f2=ΔF2/(ts[end]-ts[end-1])              # the density at t
end
pdf(d::Type{TracyWidom}, t::Real) = pdf(d(), t)

"""
Cumulative density function of the Tracy-Widom distribution

Computes the Tracy-Widom distribution by directly solving the
Painlevé II equation using the ode23 numerical integrator

See `pdf(::TracyWidom)` for a description of the arguments.
"""
function cdf{S<:Real}(d::TracyWidom, t::S, t0::S = convert(S, -8.0))
    t≤t0 && return 0.0
    t≥5  && return 1.0

    ts, y = _solve_painleve_ii(t0, t)

    F2=exp(-y[end][3])
end
cdf(d::Type{TracyWidom}, t::Real) = cdf(d(), t)

# An internal function which sets up the Painleve II differential equation and
# runs it through the ode23 numerical integrator
function _solve_painleve_ii{S<:Real}(t0::S, t::S)
    deq(t, y) = [y[2], t*y[1]+2y[1]^3, y[4], y[1]^2]
    a0 = airy(t0)
    T = typeof(a0)
    y0=T[a0, airy(1, t0), 0, airy(t0)^2]    # Initial conditions
    (ts, y)=ode23(deq, y0, [t0, t])         # Solve the ODE
end

"""
Samples the largest eigenvalue of the n × n GUE matrix
"""
function rand(d::TracyWidom, n::Int)
    n > 1 || error("Cannot construct $n × $n matrix")
    if n < 100
        k = n
    else #Exploit the fact that the eigenvector is concentrated in the top left corner
	k = round(Int, n-10*n^(1/3)-1)
    end
    a=randn(n-k+1)
    b=[χ(i) for i=(n-1):-1:k]
    v=eigmax(SymTridiagonal(a, b))
end
rand(d::Type{TracyWidom}, t::Integer) = rand(d(), t)


text:
############################################
# Default implementation `DefaultConfig` accepts nearly
# every `AbstractConfig` default, except that a
# database is used.

struct DefaultConfig <: AbstractConfig end
global const default_config = DefaultConfig()
use_db( ::DefaultConfig ) = ArrayDB

############################################
# `AlgorithmConfig` is a struct with fields defining the method outputs.
@with_kw struct AlgorithmConfig{ D <: Union{Symbol,AbstractDescentConfig} } <: AbstractConfig @deftype Float64

    eps_crit = _eps_crit( default_config )
    gamma_crit = _gamma_crit(default_config)
    max_critical_loops :: Int = max_critical_loops(default_config)

    db :: Union{Nothing,Type{<:AbstractDB}} = ArrayDB

    count_nonlinear_iterations :: Bool = count_nonlinear_iterations( default_config )

    delta_0 :: NumOrVec64 = get_delta_0(default_config)
    delta_max :: NumOrVec64 = get_delta_max(default_config)

    max_evals :: Int = max_evals( default_config )
    max_iter :: Int = max_iter( default_config )

    # relative stopping
    # stop if ||Δf|| ≤ ε ||f||
    f_tol_rel :: NumOrVec64 = f_tol_rel( default_config )
    # stop if ||Δx|| ≤ ε ||x||
    x_tol_rel:: NumOrVec64 = x_tol_rel(default_config)

    # absolute stopping
    f_tol_abs:: NumOrVec64  = f_tol_abs(default_config)
    x_tol_abs :: NumOrVec64  = x_tol_abs(default_config)

    # stop if ω ≤ omega_tol_rel && Δ .≤ Δ_tol_rel
    omega_tol_rel = omega_tol_rel(default_config)
    delta_tol_rel :: NumOrVec64 = delta_tol_rel(default_config)

    # stop if ω <= omega_tol_abs
    omega_tol_abs = omega_tol_abs(default_config)

    # stop if Δ .<= Δ_tol_abs
    delta_tol_abs :: NumOrVec64 = delta_tol_abs(default_config)

    descent_method :: D = descent_method( default_config )

    strict_acceptance_test :: Bool = strict_acceptance_test(default_config)
    nu_success = _nu_success( default_config )
    nu_accept = _nu_accept( default_config )

    mu = _mu( default_config )
    beta = _beta( default_config )

    radius_update_method :: Symbol = radius_update_method(default_config)
    gamma_grow = _gamma_grow(default_config)
    gamma_shrink = _gamma_shrink(default_config)
    gamma_shrink_much = _gamma_shrink_much(default_config)

    @assert descent_method isa AbstractDescentConfig || ( descent_method isa Symbol && descent_method ∈ [:steepest_descent, :ps, :pascoletti_serafini, :ds, :directed_search] ) "`descent_method` must be one of `:steepest_descent, :ps, :pascoletti_serafini, :ds, :directed_search`."
end

for fn in fieldnames(AlgorithmConfig)
    if fn ∉ [ :use_db, ]
        @eval $fn( ac :: AlgorithmConfig ) = getfield( ac, Symbol($fn) )
    end
end

use_db( ac :: AlgorithmConfig ) = ac.db

#####################################################
# outer constructors for the lazy user and backwards compatibility
EmptyConfig() = default_config
AlgoConfig(args...; kwargs...) = AlgorithmConfig(args...; kwargs...)

export EmptyConfig, AlgoConfig


text:
import YAML

struct ValidationInfo
    field_descriptor::Dict
    struct_name::AbstractString
    ist_struct::InfrastructureSystemsType
    field_type::Any
    limits::Union{NamedTuple{(:min, :max)}, NamedTuple{(:min, :max, :zero)}}
end

function read_validation_descriptor(filename::AbstractString)
    if occursin(r"(\.yaml)|(\.yml)"i, filename)
        data = open(filename) do file
            return YAML.load(file)
        end
    elseif occursin(r"(\.json)"i, filename)
        data = open(filename) do file
            return JSON3.read(file, Dict)
        end
    else
        error("Filename is not a YAML or JSON file.")
    end

    if data isa Dict
        if haskey(data, "auto_generated_structs")
            descriptors = data["auto_generated_structs"]
        else
            descriptors = []
        end
        if haskey(data, "struct_validation_descriptors")
            for descr in data["struct_validation_descriptors"]
                push!(descriptors, descr)
            end
        end
    else
        throw(DataFormatError("{filename} has invalid format"))
    end

    return descriptors
end

# Get validation info for one struct.
function get_config_descriptor(config::Vector, name::AbstractString)
    for item in config
        if item["struct_name"] == name
            return item
        end
    end

    @warn(
        "struct $name does not exist in validation configuration file, validation skipped"
    )
    return nothing
end

# Get validation info for one field of one struct.
function get_field_descriptor(struct_descriptor::Dict, fieldname::AbstractString)
    for field in struct_descriptor["fields"]
        if field["name"] == fieldname
            return field
        end
    end

    throw(
        DataFormatError(
            "field $fieldname does not exist in $(struct_descriptor["struct_name"]) validation config",
        ),
    )
end

function validate_fields(
    components::Components,
    ist_struct::T,
) where {T <: InfrastructureSystemsType}
    type_name = strip_parametric_type(strip_module_name(repr(T)))
    struct_descriptor = get_config_descriptor(components.validation_descriptors, type_name)
    isnothing(struct_descriptor) && return true
    is_valid = true

    for (field_name, fieldtype) in zip(fieldnames(T), fieldtypes(T))
        field_value = getfield(ist_struct, field_name)
        if isnothing(field_value)  # Many structs are of type Union{Nothing, xxx}.

        elseif fieldtype <: Union{Nothing, InfrastructureSystemsType} &&
               !(fieldtype <: InfrastructureSystemsType)
            # Recurse. Components are validated separately and do not need to
            # be validated twice.
            if !validate_fields(components, getfield(ist_struct, field_name))
                is_valid = false
            end
        else
            field_descriptor = get_field_descriptor(struct_descriptor, string(field_name))
            if !haskey(field_descriptor, "valid_range")
                continue
            end
            valid_range = field_descriptor["valid_range"]
            limits = get_limits(valid_range, ist_struct)
            valid_info = ValidationInfo(
                field_descriptor,
                struct_descriptor["struct_name"],
                ist_struct,
                fieldtype,
                limits,
            )
            if !validate_range(valid_range, valid_info, field_value)
                is_valid = false
            end
        end
    end
    return is_valid
end

function get_limits(valid_range::String, ist_struct::InfrastructureSystemsType)
    # Gets min and max values from activepowerlimits for activepower, etc.
    function recur(d, a, i = 1)
        if i <= length(a)
            d = getfield(d, Symbol(a[i]))
            recur(d, a, i + 1)
        else
            return d
        end
    end

    valid_range, ist_struct
    vr = recur(ist_struct, split(valid_range, "."))

    if isnothing(vr)
        limits = (min = nothing, max = nothing)
    else
        limits = get_limits(vr, ist_struct)
    end

    return limits
end

function get_limits(valid_range::Dict, unused::InfrastructureSystemsType)
    # Gets min and max value defined for a field,
    # e.g. "valid_range": {"min":-1.571, "max":1.571}.
    return (min = valid_range["min"], max = valid_range["max"])
end

function get_limits(
    valid_range::Union{NamedTuple{(:min, :max)}, NamedTuple{(:max, :min)}},
    unused::InfrastructureSystemsType,
)
    # Gets min and max value defined for a field,
    # e.g. "valid_range": {"min":-1.571, "max":1.571}.
    return (min = valid_range.min, max = valid_range.max)
end

function validate_range(::String, valid_info::ValidationInfo, field_value)
    # Validates activepower against activepowerlimits, etc.
    is_valid = true
    if !isnothing(valid_info.limits)
        is_valid = check_limits_impl(valid_info, field_value)
    end

    return is_valid
end

function validate_range(
    ::Union{
        Dict,
        NamedTuple{(:min, :max)},
        NamedTuple{(:max, :min)},
        NamedTuple{(:min, :max, :zero)},
    },
    valid_info::ValidationInfo,
    field_value,
)
    return check_limits(valid_info, field_value)
end

function check_limits(valid_info::ValidationInfo, field_value::Real)
    # Validates numbers.
    return check_limits_impl(valid_info, field_value)
end

function check_limits(valid_info::ValidationInfo, field_value::NamedTuple)
    # Validates up/down, min/max, from/to named tuples.
    @assert_op length(field_value) == 2
    result1 = check_limits_impl(valid_info, field_value[1])
    result2 = check_limits_impl(valid_info, field_value[2])
    return result1 && result2
end

function check_limits_impl(valid_info::ValidationInfo, field_value::Real)
    is_valid = true
    action_function = get_validation_action(valid_info.field_descriptor)
    if (
        (
            !isnothing(valid_info.limits.min) &&
            field_value + LIM_TOL < valid_info.limits.min
        ) || (
            !isnothing(valid_info.limits.max) &&
            field_value - LIM_TOL > valid_info.limits.max
        )
    ) && !(haskey(valid_info.limits, :zero) && field_value == 0.0)
        is_valid = action_function(valid_info, field_value)
    end
    return is_valid
end

function get_validation_action(field_descriptor::Dict)
    action = get(field_descriptor, "validation_action", "error")
    if action == "warn"
        action_function = validation_warning
    elseif action == "error"
        action_function = validation_error
    else
        error("Invalid validation action $action")
    end
    return action_function
end

function validation_warning(valid_info::ValidationInfo, field_value)
    valid_range = valid_info.field_descriptor["valid_range"]
    field_name = valid_info.field_descriptor["name"]
    # field_value is not displayed because it may conflict with PowerSystems unit
    # conversion in Base.show. Leaving in the call stack in case we need it later.
    @warn "Invalid range" valid_info.struct_name field_name valid_range valid_info.ist_struct
    return true
end

function validation_error(valid_info::ValidationInfo, field_value)
    valid_range = valid_info.field_descriptor["valid_range"]
    field_name = valid_info.field_descriptor["name"]
    # field_value is not displayed because it may conflict with PowerSystems unit
    # conversion in Base.show. Leaving in the call stack in case we need it later.
    @error "Invalid range" valid_info.struct_name field_name valid_range valid_info.ist_struct
    return false
end

"""
Validates a struct using only information within the struct.
"""
function validate_struct(ist::InfrastructureSystemsType)
    return true
end


text:
# # Estimating π using Monte-Carlo method
#
# ## Idea
#
# Let's compute an approximation of π using the [Monte Carlo
# method](https://en.wikipedia.org/wiki/Monte_Carlo_method).  The idea
# is to draw points from the uniform distribution on a unit square and
# count the ratio of points that are inside the quadrant of the unit
# circle.  Since the probably ``p`` for a point to be inside the
# quadrant is the area of the quadrant (``= π / 4``) divided by the
# area of the unit square (``= 1``), we can estimate ``π`` by
# estimating the probability ``p = π / 4``:

xs = rand(10_000)
ys = rand(10_000)
p = count(xs.^2 .+ ys.^2 .< 1) / length(xs)
4 * p

# > ![Illustration of monte-Carlo method for computing pi](https://upload.wikimedia.org/wikipedia/commons/8/84/Pi_30K.gif)
# >
# > --- [File:Pi 30K.gif - Wikipedia](https://en.wikipedia.org/wiki/File:Pi_30K.gif)
# > by
# > [nicoguaro](https://commons.wikimedia.org/wiki/User:Nicoguaro)
# > is licensed under
# > [CC BY 3.0](https://creativecommons.org/licenses/by/3.0/deed.en).

# ## RNG on GPU
#
# We try to do this computation on a GPU using FoldsCUDA.jl:

using CUDA
using FLoops
using FoldsCUDA

# As of writing, `CUDA.CURAND` does not provide the API usable inside
# the loop body (i.e., the device API).  However, we can use
# pure-Julia pseudo number generator quite easily.  In particular, we
# use a [Counter-based random number generator
# (CBRNG)](https://en.wikipedia.org/wiki/Counter-based_random_number_generator_(CBRNG))
# provided by [Random123.jl](https://github.com/sunoru/Random123.jl)
# ([documentation](https://sunoru.github.io/RandomNumbers.jl/stable/man/random123/)).

using Random123

# In this example, we use
# [`Random123.Philox2x`](https://sunoru.github.io/RandomNumbers.jl/stable/lib/random123/#Random123.Philox2x).
# This RNG gives us two `UInt64`s for each counter which wraps around
# at `typemax(UInt64)`:

rng_a = Philox2x(0)
rng_b = Philox2x(0)
set_counter!(rng_b, typemax(UInt64))
rand(rng_b, UInt64, 2)
@assert rng_a == rng_b

# ## Using counter-based RNG
#
# Let's create a helper function that divides
# `UInt64(0):typemax(UInt64)` into `n` equal intervals:

function counters(n)
    stride = typemax(UInt64) ÷ n
    return UInt64(0):stride:typemax(UInt64)-stride
end
nothing  # hide
#-

# This lets us use "independent" RNG for each `ctr`-th iteration:

function monte_carlo_pi(n, m = 10_000, ex = has_cuda_gpu() ? CUDAEx() : ThreadedEx())
    @floop ex for ctr in counters(n)
        rng = set_counter!(Philox2x(0), ctr)
        nhits = 0
        for _ in 1:m
            x = rand(rng)
            y = rand(rng)
            nhits += x^2 + y^2 < 1
        end
        @reduce(tot = 0 + nhits)
    end
    return 4 * tot / (n * m)
end
nothing  # hide
#-

πₐₚₚᵣₒₓ = monte_carlo_pi(2^12)


text:
#=
   Likelihhod is chi^2 comparing observations of the radial velocity of the star to a model,
   where the model computes the velocity of the star as the linear superposition of the Keplerian orbit
   induced by each planet, i.e., neglecting mutual planet-planet interactions
   Priors based on SAMSI reference priors from Ford & Gregory 2006
=#

module RvModelKeplerian

include("rv_model_param.jl")       # Code to convert parameter vector to named parameters
include("utils_internal.jl")       # Manually add code to support autodiff of mod2pi
include("rv_model_keplerian.jl")   # Code to compute RV model for Keplerian orbits
include("stat_model.jl")           # Code to compute target density, plogtarget

end # module RvModelKeplerian


text:
########################################
#            Side constants            #
########################################

const TOP = Top()
const LEFT = Left()
const RIGHT = Right()
const BOTTOM = Bottom()

################################################################################
#                              Utility functions                               #
################################################################################

to2tuple(x1, x2) = (x1, x2)

"""
Recursively convert every element in the collection to a Point2f0.
"""
rec_point(x::Vector{<: AbstractFloat}) = Point2f0(x)
rec_point(x) = rec_point.(x)

"""
Recurses through the given structure until it reaches a Point2f0, then projects it.
"""
rec_project(source::Projection, dest::Projection, array) = rec_project.(source, dest, array)
rec_project(source::Projection, dest::Projection, point::Point2f0) = transform(source, dest, point)

########################################
#          Grid tessellation           #
########################################

"""
    grid_triangle_faces(lats, lons)

Returns a series of triangle indices from naive triangulation.
"""
function grid_triangle_faces(lons, lats)
    faces = Array{Int, 2}(undef, (length(lons)-1)*(length(lats)-1)*2, 3)

    xmax = length(lons)

    i = 1
    # TODO optimize this with broadcast
    for lon in eachindex(lons)[1:end-1]

        for lat in eachindex(lats)[1:end-1]

            cpos = lon + (lat-1)*xmax

            faces[i, :] = [cpos, cpos+1, cpos+xmax+1]

            faces[i+1, :] = [cpos, cpos+xmax, cpos+xmax+1]

            i += 2

        end
    end

    return faces

end

"""
    gridpoints(xs, ys)

Returns a Vector of Points of a grid formed by xs and ys.
"""
gridpoints(xs, ys) = vec([Point2f0(x, y) for y in ys, x in xs])

"""
    triangulated_grid(xs, ys) -> points, faces

Takes in two ranges, and returns a triangulated regular grid based on those ranges.
"""
triangulated_grid(xs, ys) = (gridpoints(xs, ys), grid_triangle_faces(xs, ys))

date_regex(dirname, ext) = Regex("$(dirname)_(\\d{4})-(\\d{2}).$(uppercase(ext))")
imflip(img) = reverse(vec(transpose(reverse(img; dims=2))))

########################################
#       Nested array flattening        #
########################################

"""
    to_nansep_vec([f::Function,] data::Vector{Vector{T}}) where T

Flattens the given Vector of Vectors into a single Vector,
while inserting NaN separations between each individual sub-
vector.

If `f` is given, then `f` will be executed on each subvector
before it is merged into the main vector.  If it is not, the data
will remain unchanged.
"""
function to_nansep_vec(f::Function, data::AbstractVector{T}, outtyp = Point2f0) where T

    length_of_data = sum(length.(data)) + length(data)

    lvec = Vector{outtyp}(undef, length_of_data)

    pos = 1

    for (i, datum) in enumerate(data)
        lvec[pos:(pos+length(datum)-1)] .= f(datum)
        pos += length(datum)
        lvec[pos] = outtyp(NaN)
        pos += 1
    end

    return lvec

end

to_nansep_vec(data::Vector{Vector{T}}, outtyp = Point2f0) where T = to_nansep_vec(identity, data, outtyp)

########################################
#     Matrix-based grid generation     #
########################################

"""
    xygrid(lons::Vector, lats::Vector; projection = LatLon()) -> (x::Matrix, y::Matrix)

Generates a grid of x and y coordinates in lat-lon space from the given ranges.

Optionally, also projects them into the specified CRS before returning them.

!!! note
    Materializing the matrices is necessary for C-interoperability,
    since passing the matrix to C implies that it can't
    be stored lazily and generated on the fly.
"""
function xygrid(lons::AbstractVector{<: Number}, lats::AbstractVector{<: Number}; projection = LonLat(), source = LonLat())

    lats = Float64.(lats)
    lons = Float64.(lons)

    xs = [lon for lon in lons, lat in lats]  # xs / longitudes
    ys = [lat for lon in lons, lat in lats]  # ys / latitudes

    # this will internally reinterpret the matrices as vectors,
    # and transform them in place in C.  The matrices will
    # then hold the modified values.
    # Proj4.transform!(source, projection, vec(xs), vec(ys))

    return xs, ys

end

"""
    wrap_lon(x)

Wraps the input longitude (or Vector of longitudes)
into the domain `(-180, 180)`.
"""
wrap_lon(x) = -180 .+ (360 .+ ((x.+180) .% 360)) .% 360

# function Base.:(==)(a::Projection, b::Projection)
#     return Bool(
#         Proj4.proj_is_equivalent_to(
#             a.rep,
#             b.rep,
#             Proj4.PJ_COMP_EQUIVALENT
#         )
#     )
# end


text:
#=

  P-median problem in ConstraintSolver.jl

  Model and data from the OPL Manual, which describes the problem:
  """
  The P-Median problem is a well known problem in Operations Research.
  The problem can be stated very simply, like this: given a set of customers
  with known amounts of demand, a set of candidate locations for warehouses,
  and the distance between each pair of customer-warehouse, choose P
  warehouses to open that minimize the demand-weighted distance of serving
  all customers from those P warehouses.
  """

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

using ConstraintSolver, JuMP
using Cbc, GLPK, Ipopt
const CS = ConstraintSolver
include("constraints_utils.jl")

function p_median(problem, print_solutions=true,all_solutions=true,timeout=6)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)
    ipopt_optimizer = optimizer_with_attributes(Ipopt.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> all_solutions,
                                                            # "all_optimal_solutions"=>all_solutions,
                                                            "logging"=>[],

                                                            "traverse_strategy"=>:BFS,
                                                            # "traverse_strategy"=>:DFS,
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf,

                                                            # https://wikunia.github.io/ConstraintSolver.jl/stable/options/#branch_strategy-(:Auto)
                                                            "branch_strategy" => :IMPS, # default
                                                            # "branch_strategy" => :ABS, # Activity Based Search
                                                            # "activity.decay" => 0.999, # default 0.999
                                                            # "activity.max_probes" => 10, # default, 10
                                                            # "activity.max_confidence_deviation" => 20, # default 20

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>timeout,

                                                            # "backtrack" => false, # default true
                                                            # "backtrack_sorting" => false, # default true

                                                            # "lp_optimizer" => cbc_optimizer,
                                                            # "lp_optimizer" => glpk_optimizer,
                                                            # "lp_optimizer" => ipopt_optimizer,
                                        ))
    demand = problem[:demand]
    distance = problem[:distance]
    p = problem[:p]
    num_customers,num_warehouses = size(distance)

    println("p:$p")
    println("demand:$demand")
    println("distance:$distance")
    println("num_warehouses:$num_warehouses num_customers:$num_customers")

    @variable(model, open_warehouse[1:num_warehouses], Bin)
    @variable(model, ship_to_customer[1:num_customers,1:num_warehouses], Bin)
    @variable(model, 0 <= z <= sum(demand)*num_customers, Int)

    @constraint(model, z == sum(demand[c]*distance[c,w]*ship_to_customer[c,w]
                            for c in 1:num_customers, w in 1:num_warehouses))

    for c in 1:num_customers
        @constraint(model, sum([ship_to_customer[c,w] for w in 1:num_warehouses]) == 1)
    end

    @constraint(model, p == sum(open_warehouse))

    for c in 1:num_customers, w in 1:num_warehouses
        @constraint(model, ship_to_customer[c,w] <= open_warehouse[w])
    end

    @objective(model,Min,z)

    # Solve the problem
    optimize!(model)

    status = JuMP.termination_status(model)
    # println("status:$status")
    num_sols = 0
    if status == MOI.OPTIMAL
        num_sols = MOI.get(model, MOI.ResultCount())
        println("num_sols:$num_sols\n")
        if print_solutions
            for sol in 1:num_sols
                println("solution #$sol")
                open_warehouse_val = convert.(Integer,JuMP.value.(open_warehouse; result=sol))
                ship_to_customer_val = convert.(Integer,JuMP.value.(ship_to_customer; result=sol))
                z_val = convert.(Integer,JuMP.value.(z; result=sol))
                println("open_warehouse:$open_warehouse_val")
                println("ship_to_customer:$ship_to_customer_val")
                println("z:$z_val")
                println()

            end
        end
    else
        println("status:$status")
    end

    return status, num_sols
end

problem = Dict(
    :demand => [100,80,80,70],
    :distance => resize_matrix([[2, 10, 50],
                                [2, 10, 52],
                                [50, 60,  3],
                                [40, 60,  1]]),
    :p => 2,
)
@time p_median(problem)

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
using Test
using ChartParsers
using ChartParsers: ActiveArc, PassiveArc, rule, start, stop, constituents, Chart, combine, head, AbstractGrammar

@testset "arc in chart detection" begin
    chart = Chart{SimpleRule, Symbol}(2)
    a1 = PassiveArc(Arc(0, 1, SimpleRule(:NP => Symbol[]), [], 1))
    a2 = ActiveArc(Arc(0, 0, SimpleRule(:S => [:NP, :VP]), [], 1))
    a3 = combine(a2, a1)

    @test a1 ∉ chart
    @test a2 ∉ chart
    @test a3 ∉ chart
    push!(chart, a1)
    @test a1 ∈ chart
    @test a2 ∉ chart
    @test a3 ∉ chart
    push!(chart, a2)
    @test a1 ∈ chart
    @test a2 ∈ chart
    @test a3 ∉ chart
    push!(chart, a3)
    @test a1 ∈ chart
    @test a2 ∈ chart
    @test a3 ∈ chart
end

@testset "Example from nlp-with-prolog" begin
    # example taken from http://cs.union.edu/~striegnk/courses/nlp-with-prolog/html/node71.html#l9.sec.bottomup

    grammar = SimpleGrammar([
            :S => [:NP, :VP, :PP],
            :S => [:NP, :VP],
            :NP => [:PN],
            :VP => [:IV],
            :PP => [:P, :NP],
        ], Dict(
            "mia" => [:PN],
            "danced" => [:IV]
        ), :S)

    tokens = split("mia danced")

    @testset "BottomUp" begin
        parser = ChartParser(tokens, grammar, BottomUp())
        parses = collect(parser)
        @test length(parses) == 5

        complete_parses = collect(Iterators.filter(is_complete(parser), parser))
        @test length(complete_parses) == 1
        p = first(complete_parses)
        @test head(p) == :S
        @test rule(p).rule == (:S => [:NP, :VP])
        @test rule(constituents(p)[1]).rule == (:NP => [:PN])
        @test rule(constituents(p)[2]).rule == (:VP => [:IV])
    end

    @testset "TopDown" begin
        parser = ChartParser(tokens, grammar, TopDown())
        parses = collect(parser)

        # The top down parser currently doesn't yield the terminal productions,
        # but the bottom up parser does. Should we change that?
        @test length(parses) == 3

        complete_parses = collect(Iterators.filter(is_complete(parser), parser))
        @test length(complete_parses) == 1
        p = first(complete_parses)
        @test head(p) == :S
        @test rule(p).rule == (:S => [:NP, :VP])
        @test rule(constituents(p)[1]).rule == (:NP => [:PN])
        @test rule(constituents(p)[2]).rule == (:VP => [:IV])
    end
end

@testset "longer example" begin
    tokens = split("mary sat on the table yesterday")
    grammar = SimpleGrammar([
            :S => [:NP, :VP],
            :NP => [:PN],
            :VP => [:V, :NP],
            :VP => [:V, :PP],
            :VP => [:VP, :AV],
            :PP => [:P, :NP],
            :NP => [:D, :N]
        ], Dict(
        "mary" => [:PN],
        "sat" => [:V],
        "on" => [:P],
        "the" => [:D],
        "table" => [:N],
        "yesterday" => [:AV]
        ), :S)
    parser = ChartParser(tokens, grammar)
    complete_parses = @inferred collect(Iterators.filter(is_complete(parser), parser))
    @test length(complete_parses) == 1

    p = first(complete_parses)
    @test head(p) == :S
    @test rule(p).rule == (:S => [:NP, :VP])
    @test rule(constituents(p)[1]).rule == (:NP => [:PN])
    @test rule(constituents(p)[2]).rule == (:VP => [:VP, :AV])
    @test rule(constituents(constituents(p)[2])[1]).rule == (:VP => [:V, :PP])
end

abstract type GrammaticalSymbol end
struct Token <: GrammaticalSymbol end
struct S <: GrammaticalSymbol end
struct VP <: GrammaticalSymbol end
struct NP <: GrammaticalSymbol end

struct TypedRule <: AbstractRule{GrammaticalSymbol}
    rule::Pair{GrammaticalSymbol, NTuple{N, GrammaticalSymbol} where N}
end
ChartParsers.rhs(t::TypedRule) = last(t.rule)
ChartParsers.lhs(t::TypedRule) = first(t.rule)
ChartParsers.score(t::TypedRule) = 1.0

struct TypedGrammar <: AbstractGrammar{TypedRule}
    productions::Vector{TypedRule}
    labels::Dict{String, Vector{GrammaticalSymbol}}
end

ChartParsers.productions(g::TypedGrammar) = g.productions

function ChartParsers.terminal_productions(g::TypedGrammar, tokens)
    R = ChartParsers.rule_type(g)
    result = Arc{R}[]
    for (i, token) in enumerate(tokens)
        for label in get(g.labels, token, GrammaticalSymbol[])
            push!(result, Arc{R}(i - 1, i, TypedRule(label => (Token(),)), [token], 1))
        end
    end
    result
end

ChartParsers.start_symbol(g::TypedGrammar) = S()

@testset "typed rules" begin
    tokens = split("mia danced")
    grammar = TypedGrammar(TypedRule.([
            S() => (VP(), NP()),
            S() => (NP(), VP())
        ]),
        Dict("mia" => [NP()], "danced" => [VP()]))
    parser = ChartParser(tokens, grammar)
    complete_parses = @inferred collect(Iterators.filter(is_complete(parser), parser))
    @test length(complete_parses) == 1
    @test rule(first(complete_parses)).rule == (S() => (NP(), VP()))
end

@testset "Weighted terminal productions" begin
    grammar = SimpleWeightedGrammar([
        (:A => [:B, :C], 1),
        (:B => [:E], 1),
        (:C => [:F], 0.5),
        (:C => [:G], 0.5)],
        Dict("a" => [:D => 0.6, :E => 0.4],
             "b" => [:F => 0.9, :G => 0.1]),
        :A)
    tokens = split("a b")
    parser = ChartParser(tokens, grammar, BottomUp())
    parses = collect(parser)
    @test [rule(p).rule for p in parses] == [
        :F => Symbol[],
        :D => Symbol[],
        :C => [:F],
        :E => Symbol[],
        :B => [:E],
        :A => [:B, :C],
        :G => Symbol[],
        :C => [:G],
        :A => [:B, :C]
    ]

    p = first(filter(is_complete(parser), parses))
    @test [rule(c).rule for c in constituents(p)] == [:B => [:E], :C => [:F]]
end


text:
module NLreg

    using DataArrays, StatsBase, DataFrames
    using Distributions: FDist, ccdf
    using Base.LinAlg.QRCompactWY

    export     # Partially linear models
        AsympOff,     # Asymptotic regression expressed as an offset
        AsympOrig,    # Asymptotic regression through origin
        AsympReg,     # Asymptotic regression model
        Biexp,        # Bi-exponential model
        BolusSD1,     # 1-compartment, single bolus dose model
        Chwirut,      # Model used in the NIST nonlinear regression examples
        Gompertz,     # Gompertz growth model
        Logis3P,      # 3-parameter logistic
        Logis4P,      # 4-parameter logistic
        MicMen,       # Michaelis-Menten model
              # types for models functions and model fits
        CompositeModF,   # Nonlinear regression model function with parameter transformation
        CompositePLModF, # Partially linear regression model function with parameter trans
        NLregModF,       # Nonlinear regression model function
        NonlinearLS,     # Nonlinear least squares model
        PLinearLS,       # Partially linear regression model
        PLregModF,       # Partially linear regression model function
        NLMM,            # Nonlinear mixed-effects model
        SimpleNLMM,      # Simple population nonlinear mixed-effects model
              # parameter transformations
        Dtrans,          # diagonal parameter transformation
        IdTr,            # scalar identity transformation
        LogTr,           # scalar logarithmic transformation
        Ptrans,          # abstract parameter transformation
        Strans,          # scalar parameter transformation
              # Full, nonlinear regression models
        Logsd1,       # 1-compartment, single bolus dose model with logged parameters

        formula,      # define the nonlinear model formula
        incr!,        # increment the spherical random effects
        initpars,     # create initial values for the parameters
        lowerbd,      # lower bounds on NLMM parameters
        npars,        # total number of parameters in an NLregMod
        pnames,       # names of parameters in a model
        pnls!,        # penalized nonlinear least squares fit
        prss!,        # penalized rss for b = lambda * (u + fac*delu)
        setpars!,     # set new parameter values (beta + theta) in an NLMM
        theta,        # extract covariance parameters
        theta!,       # set covariance parameters
        updtMM!,      # update the model matrix in a PLregMod
        updtL!,       # update L and solve for delu
        updtmu!       # update mu and tgrad

    typealias FP FloatingPoint

    include("nonlinreg.jl")
    include("plreg.jl")
    include("models.jl")
#    include("nlmm.jl")
#    include("simplenlmm.jl")
    include("ptrans.jl")
    include("compositemodel.jl")

end # module


text:
using VideoIO
using Images
using ImageDraw

io = VideoIO.openvideo("video.mp4") # for testing purposes

function posindex(x)
	return floor(Int, x%600),ceil(Int, x/600)
end

function firstlast(imgdash)
	dia1pos = 0
	dia2pos = 0
	changed = true
	for (id,i) in enumerate(imgdash)
		if i == 1 && changed == true
			dia1pos = posindex(id)
			changed = false
		elseif i == 1
			dia2pos = posindex(id)
		else
			continue
		end
	end
	return dia1pos, dia2pos
end

function centerradius(dia1, dia2)
	center = (floor(Int, (dia1[1]+dia2[1])/2 ),floor(Int, (dia1[2]+dia2[2])/2))
	radius = floor(Int,sqrt((dia1[1]-dia2[1])^2 + (dia1[2]-dia2[2])^2) /2)
	return center, radius
end
centers = []
function framestackbuilder(i)
	img = convert(Array{Float64}, Gray.(i)) .> 0.5
	first, last = firstlast(img)
	center, radius = centerradius(first, last)
	append!(centers, center)
	draw!(RGB.(img), Polygon(RectanglePoints(Point(center[2]-radius-5, center[1]-radius-5), Point(center[2]+radius+5, center[1]+radius+5))), RGB(0,1,0))
	if length(centers) < 7
		a =  1
	else
		a =  length(centers)-6
	end
	for i in centers[a:end]
		draw!(img, LineSegment(Point(10, 10), Point(100, 100)), RGB{N0f8}(1))
	end
	return img
end
# framestack = map(x->framestackbuilder(x), io)
encoder_options = (crf=23, preset="medium")
framerate=24
img1  = load("./lotofdata/screenshot_1.jpeg")
open_video_out("video1.mp4", img1, framerate=framerate, encoder_options=encoder_options) do writer
    for frame in io
		img = convert(Array{Float64}, Gray.(frame)) .> 0.5
		first, last = firstlast(img)
		center, radius = centerradius(first, last)
		append!(centers, center)
		img = draw(RGB.(img), Polygon(RectanglePoints(Point(center[2]-radius-5, center[1]-radius-5), Point(center[2]+radius+5, center[1]+radius+5))), RGB(0,1,0))
		if length(centers) < 7
			a =  1
		else
			a =  length(centers)-6
		end
		for i in centers[a:end]
			img = draw(img, LineSegment(Point(i[1], i[1]), Point(i[2]-100, i[2]-100)), RGB{N0f8}(1))
		end
        write(writer, RGB.(img))
    end
end

max_stars_repo_name: ashwani-rathee/Sigil

max_issues_repo_name: ashwani-rathee/Sigil

max_forks_repo_name: ashwani-rathee/Sigil


text:
struct ObjectiveData{C,X,U,XX,UU,UX}
    costs::C
    gradient_state::Vector{X}
    gradient_action::Vector{U}
    hessian_state_state::Vector{XX}
    hessian_action_action::Vector{UU}
    hessian_action_state::Vector{UX}
end

function objective_data(dynamics::Vector{Dynamics{T}}, costs) where T
	gradient_state = [[zeros(d.num_state) for d in dynamics]...,
        zeros(dynamics[end].num_next_state)]
    gradient_action = [zeros(d.num_action) for d in dynamics]
    hessian_state_state = [[zeros(d.num_state, d.num_state) for d in dynamics]...,
        zeros(dynamics[end].num_next_state, dynamics[end].num_next_state)]
    hessian_action_action = [zeros(d.num_action, d.num_action) for d in dynamics]
    hessian_action_state = [zeros(d.num_action, d.num_state) for d in dynamics]
    ObjectiveData(costs, gradient_state, gradient_action, hessian_state_state, hessian_action_action, hessian_action_state)
end

function reset!(data::ObjectiveData)
    H = length(data.gradient_state)
    for t = 1:H
        fill!(data.gradient_state[t], 0.0)
        fill!(data.hessian_state_state[t], 0.0)
        t == H && continue
        fill!(data.gradient_action[t], 0.0)
        fill!(data.hessian_action_action[t], 0.0)
        fill!(data.hessian_action_state[t], 0.0)
    end
end


text:
#= Possible components for the gui
    - window displaying the simulation
    - ⏵/⏸ button for the simulation
    - slider for variables (temperature, adhesiveness, etc)
    - display for much of MH steps
    - count for number of cells
=#

function CellGUI(CPM::CellPotts{2})

    #--------Simulation Screen--------
    #Start by creating a blank figure
    fig = Figure(resolution = (1200, 1200), backgroundcolor = RGBf0(0.98, 0.98, 0.98))

    #Create observables that will change as the simulation progresses
    timestep = Node(1) #will increase by one every step
    maxLength = 1000 #Window history for plot
    energies = Node(fill(CPM.energy, maxLength)) #vector of past CPM emergies

    #The energies plot updates overtime and needs dynamic axes limits
    lim = @lift begin
        xmin = max(0, $timestep - length($energies))
        xmax = max(length($energies), $timestep)
        ymin, ymax =  extrema($energies)

        (xmin, xmax, ymin-1, ymax+1)
    end

    #Name the axes on the figure
    axSim = fig[1, 1] = Axis(fig, title = "Simulation")
    axEnergy = fig[1, 2] = Axis(fig, title = "Energy", limits = lim, height = 300, tellheight = false, valign = :top)
    colsize!(fig.layout, 1, Relative(2/3)) #Make cells heatmap plot relatively larger

    #The first plot will show a heatmap of the cell simulation
    # heatmap_node is an array that updates when timestep updates
    heatmap_node = @lift begin
        currentTime = $timestep
        MHStep!(CPM)
        CPM.visual
    end

    #Create the heatmap
    heatmap!(axSim,
            heatmap_node,
            show_axis = false,
            colormap = :Purples) #:Greys_3
    tightlimits!.(axSim)
    hidedecorations!.(axSim) #removes axis numbers

    #Similiarly, overlay a heatmap showing the articulation points
    artic_node = @lift begin
        currentTime = $timestep
        reshape(CPM.graph.isArticulation, CPM.M.graphDimension)
    end

    #Give transparency to the articulation points
    cm = [RGBA(0,0,0,0.0), RGBA(1,0,0,0.5)]
    heatmap!(axSim,
            artic_node,
            show_axis = false,
            colormap = cm)

    #Finally, add the edge borders to the cells
    edgeConnectors = Edge2Grid(CPM.M.graphDimension)
    (m,n) = CPM.M.graphDimension

    #Generate all of the edge Connections by putting a point on each cell corner
    horizontal = [Point2f0(x, y) => Point2f0(x+1, y) for x in 0.5:m-0.5, y in 0.5:m+0.5]
    vertical = [Point2f0(x, y) => Point2f0(x, y+1) for y in 0.5:n-0.5, x in 0.5:n+0.5]
    points = vcat(horizontal[:],vertical[:])

    #Determine the transparency of the linesegments
    gridflip = rotl90(CPM.visual) #https://github.com/JuliaPlots/Makie.jl/issues/205

    #Cell borders are outlined in black
    black = RGBA{Float64}(0.0,0.0,0.0,1.0);
    clear = RGBA{Float64}(0.0,0.0,0.0,0.0);

    #Loop through all the grid connected and assign the correct color
    currentEdgeColors = [gridflip[edges[1]]==gridflip[edges[2]] ? clear : black for edges in edgeConnectors];

    #For each time update, recolor all of the edges
    lineColors_node = @lift begin
        currentTime = $timestep

        gridflip = rotl90(CPM.visual)

        for (i,edges) in enumerate(edgeConnectors)
            currentEdgeColors[i] = gridflip[edges[1]]==gridflip[edges[2]] ? clear : black
        end

        currentEdgeColors
    end

    #Plot all of the line segments onto the simulation
    linesegments!(
        axSim,
        points,
        color = lineColors_node,
        linewidth = 2
    )

    #On the 2nd axis, plot the energies
    xvals = @lift( $lim[1]+1:$lim[2] )
    lines!(axEnergy, xvals, energies)

    #--------Buttons--------
    #Currently 2 buttons: a play/pause and a stop button
    #Place the buttons below the simulation and align to the left
    fig[2,1] = buttongrid = GridLayout(tellwidth = false, halign = :left)
    buttonsLabels = ["▶","■","Divide!","Kill"]

    #Loop through and assign button labels, width, and color
    buttongrid[1,1:length(buttonsLabels)] = [Button(
        fig,
        label = lab,
        width = 70,
        buttoncolor = :grey) for lab in buttonsLabels]

    #Set the buttons to individual variables
    playpause, stop, cellDivideButton, cellDeathButton = contents(buttongrid)

    #If the play/pause button is clicked, change the label
    on(playpause.clicks) do clicks
        if isodd(clicks)
            playpause.label = "𝅛𝅛"
        else
            playpause.label = "▶"
        end
    end

    #partition a random cell when button is clicked
    on(cellDivideButton.clicks) do clicks
        CellDivision!(CPM,rand(1:maximum(CPM.cell.ids)))
    end

    #Choose a random cell to kill
    on(cellDeathButton.clicks) do clicks
        CellDeath!(CPM,rand(1:maximum(CPM.cell.ids)))
    end

    display(fig)

    runsim = true
    stop.clicks[] = 0 #for stop button
    while runsim

        #Is the pause button pushed?
        if playpause.label[] == "▶"
            timestep[] += 1
            appendEnergy!(energies,CPM)
            notify(timestep)
            notify(energies)
        end

        #Has the stop button been pushed?
        if  stop.clicks[] == 1
            runsim = false
            GLMakie.destroy!(GLMakie.global_gl_screen()) #close the window
        end

        sleep(eps())
    end
end

function appendEnergy!(energies,CPM)
    popfirst!(energies[])
    push!(energies[], CPM.energy)
end

#This is very ugly and maybe one day I'll make it better
#This function takes in the grid dimension and returns pairs of all adjacent squares (wraps around)
# ╔═══╤═══╗
# ║ 1 │ 3 ║
# ╠═══╪═══╣
# ║ 2 │ 4 ║
# ╚═══╧═══╝
# Edge2Grid(2) = [[2, 1], [4, 3], [1, 2], [3, 4], [2, 1], [4, 3], [2, 4], [1, 3], [4, 2], [3, 1], [2, 4], [1, 3]]

function Edge2Grid(dim)
    gridIndices = LinearIndices(dim)

    x1 = reverse(reshape(gridIndices,dim),dims=1)'[:]
    x2 = circshift(x1,dim[2])

    y1 = reverse(reshape(reverse(gridIndices),dim),dims=2)[:]
    y2 = circshift(y1,dim[1])

    append!(x1,x1[1:dim[1]])
    append!(x2,x2[1:dim[1]])
    append!(y1,y1[1:dim[1]])
    append!(y2,y2[1:dim[1]])

    return [[id1,id2] for (id1,id2) in zip([x1;y1],[x2;y2])]
end

#3D method for GUI
function CellGUI(CPM::CellPotts{3})

    timestep = Node(1)

    # SSAO attributes are per scene (will need to play with these), too slow for animation?
    scene = Scene()

    GLMakie.enable_SSAO[] = true
    scene[:SSAO][:radius][] = 5.0
    scene[:SSAO][:blur][] = 3
    scene[:SSAO][:bias][] = 0.025

    #General voxel
    voxel = Rect3D(Point3f0(-0.5), Vec3f0(1))

    lim = FRect3D((0,0,0), CPM.M.graphDimension)

    mesh_node = @lift begin
        currentTime = $timestep
        MHStep!(CPM)
         #Positions for the voxels
         #use the cell indices to color the cell (could also use cell type)
        [Point3f0(idx.I...) for idx in CartesianIndices(CPM.visual) if CPM.visual[idx] ≠ 0]
    end

    color_node = @lift begin
        currentTime = $timestep
        colors = filter(!isequal(0), CPM.visual)
    end

    meshscatter!(scene,
        mesh_node,
        marker=voxel,
        markersize=1,
        color=color_node,
        colormap=:RdYlBu_11, #see https://juliagraphics.github.io/ColorSchemes.jl/stable/basics/
        limits = lim,
        ssao=true)

    display(scene)

        runsim = true
        while runsim
            timestep[] += 1
            notify(timestep)
            sleep(eps())
        end
end


text:
# getFullsizeModelData: getting the data of full-size model which is entered by users
function getFullsizeModelData(m::JuMP.Model)
  # parameter
  # input:
  # m: the model entered by users using JuMP
  # output:
  # fullsizeObjSense: the objective sense of full-size model (Max or Min)
  # fullsizeObjCoeff: the coefficients of objective function of full-size model
  # fullsizeConstLB: the lower bound of each constraint of full-size model, which is a vector
  # fullsizeConstUB: the upper bound of each constraint of full-size model, which is a vector
  # fullsizeConstMatrix: the constraint matrix of full-size model
  # fullsizeVarLB: the lower bound of each varible of full-size model, which is a vector
  # fullsizeVarUB: the upper bound of each varible of full-size model, which is a vector
  # fullsizeNumVar: number of variables of full-size model
  # fullsizeNumConst: number of constraints of full-size model

  fullsizeObjSense = m.objSense
  fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB = JuMP.prepProblemBounds(m)
  fullsizeConstMatrix = JuMP.prepConstrMatrix(m)
  fullsizeVarLB = m.colLower
  fullsizeVarUB = m.colUpper

  fullsizeNumVar = m.numCols
  fullsizeNumConst = length(fullsizeConstLB)

  return fullsizeObjSense, fullsizeObjCoeff, fullsizeConstLB, fullsizeConstUB, fullsizeConstMatrix, fullsizeVarLB, fullsizeVarUB, fullsizeNumVar, fullsizeNumConst
end

# identifying formulation of model, which is envelopment form or multiplier form
function identifyModelForm(fullsizeObjCoeff::Vector{Float64})
  # parameters
  # input
  # fullsizeObjCoeff: the coefficients of objective function of full-size model
  # output
  # modelForm: formulation of model, which is envelopment form or multiplier form
  modelForm = :Envelopment
  objSum = sum(fullsizeObjCoeff)

  # if summatiom of coefficients of objective function is 1 or -1, the model is envelopment form
  if objSum == 1 || objSum == -1
    modelForm = :Envelopment
  else
    modelForm = :Multiplier
  end

  return modelForm
end

# getting the right hand side of model by its constraints lower and upper bound
function getRHS(constLB::Vector{Float64}, constUB::Vector{Float64})
  # parameters
  # input
  # constLB: the lower bound of each constraint of  model, which is a vector
  # constUB: the upper bound of each constraint of  model, which is a vector
  # output
  # rhs: right hand side of model, which is a vector

  numConst = length(constLB)
  rhs = Array(Float64, numConst)

  # if the lower bound of constraint is -Inf, the constraint is "<=" its rhs and its rhs is its upper bound
  lb = find(x -> x != -Inf, constLB)
  rhs[lb] = constLB[lb]
  # if the upper bound of constraint is Inf, the constraint is ">=" its rhs and its rhs is its lower bound
  ub = find(x -> x != Inf, constUB)
  rhs[ub] = constUB[ub]

  return rhs
end

# transforming model from mutiplier form to envelopment form
function modelTransform(fullsizeObjSense::Symbol, fullsizeObjCoeff::Vector{Float64}, fullsizeConstLB::Vector{Float64}, fullsizeConstUB::Vector{Float64},
                        fullsizeConstMatrix::SparseMatrixCSC, fullsizeVarLB::Vector{Float64}, fullsizeVarUB::Vector{Float64}, fullsizeNumVar::Int, fullsizeNumConstr::Int)
  # parameter
  # input
  # data of model which is multiplier form
  # output
  # data of model which is transformed to envelopment form

  # transforming coefficients of objective function
  newObjCoeff = getRHS(fullsizeConstLB, fullsizeConstUB)
  # transforming constraint matrix
  newMatrix = transpose(fullsizeConstMatrix)
  newNumConst, newNumVar = size(newMatrix)

  # initialization
  newVarLB = fill(NaN, newNumVar)
  newVarUB = fill(NaN, newNumVar)
  newConstLB = fill(NaN, newNumConst)
  newConstUB = fill(NaN, newNumConst)

  # transforming objective sense
  newObjSense = :Min

  # the rule of tranforming is different by objective sense
  if fullsizeObjSense == :Max
    newObjSense = :Min

    # transforming lower bound and upper bound of variables
    # finding the indices of constraints whose senses are ">="
    constSenseMore = find(x -> x == Inf, fullsizeConstUB)
    newVarLB[constSenseMore] = -Inf
    newVarUB[constSenseMore] = 0

    # finding the indices of constraints whose senses are "<="
    constSenseLess = find(x -> x == -Inf, fullsizeConstLB)
    newVarLB[constSenseLess] = 0
    newVarUB[constSenseLess] = Inf

    # finding the indices of constraints whose senses are "=="
    constSenseEqual = find(x -> isnan(x), newVarLB)
    newVarLB[constSenseEqual] = -Inf
    newVarUB[constSenseEqual] = Inf

    # transforming lower bound and upper bound of variables
    # finding the indices of variables which are ">= 0"
    varSenseMore = find(x -> x == 0, fullsizeVarLB)
    newConstLB[varSenseMore] = fullsizeObjCoeff[varSenseMore]
    newConstUB[varSenseMore] = Inf

    # finding the indices of variables which are "<= 0"
    varSenseLess = find(x -> x == 0, fullsizeVarUB)
    newConstLB[varSenseLess] = -Inf
    newConstUB[varSenseLess] = fullsizeObjCoeff[varSenseLess]

    # finding the indices of variables which are unbounded
    vrsUnbounded = find(x -> isnan(x), newConstLB)
    newConstLB[vrsUnbounded] = fullsizeObjCoeff[vrsUnbounded]
    newConstUB[vrsUnbounded] = fullsizeObjCoeff[vrsUnbounded]

  else
    newObjSense = :Max

    # transforming lower bound and upper bound of variables
    # finding the indices of constraints whose senses are ">="
    constSenseMore = find(x -> x == Inf, fullsizeConstUB)
    newVarLB[constSenseMore] = 0
    newVarUB[constSenseMore] = Inf

    # finding the indices of constraints whose senses are "<="
    constSenseLess = find(x -> x == -Inf, fullsizeConstLB)
    newVarLB[constSenseLess] = -Inf
    newVarUB[constSenseLess] = 0

    # finding the indices of constraints whose senses are "=="
    constSenseEqual = find(x -> isnan(x), newVarLB)
    newVarLB[constSenseEqual] = -Inf
    newVarUB[constSenseEqual] = Inf

    # transforming lower bound and upper bound of variables
    # finding the indices of variables which are ">= 0"
    varSenseMore = find(x -> x == 0, fullsizeVarLB)
    newConstLB[varSenseMore] = -Inf
    newConstUB[varSenseMore] = fullsizeObjCoeff[varSenseMore]

    # finding the indices of variables which are "<= 0"
    varSenseLess = find(x -> x == 0, fullsizeVarUB)
    newConstLB[varSenseLess] = fullsizeObjCoeff[varSenseLess]
    newConstUB[varSenseLess] = Inf

    # finding the indices of variables which are unbounded
    vrsUnbounded = find(x -> isnan(x), newConstLB)
    newConstLB[vrsUnbounded] = fullsizeObjCoeff[vrsUnbounded]
    newConstUB[vrsUnbounded] = fullsizeObjCoeff[vrsUnbounded]
  end

  return newObjSense, newObjCoeff, newConstLB, newConstUB, newMatrix, newVarLB, newVarUB, newNumVar, newNumConst
end

# finding theta, which is the index w.r.t variable vector
# parameters
# input
# fullsizeObjCoeff: the coefficients of objective function of full-size model
# output: the index of theta w.r.t variable vector
findThetaIdx(fullsizeObjCoeff::Vector{Float64}) = find(x -> x != 0, fullsizeObjCoeff)[1]

# getting the coefficients of constraints of target DMU
function getTargetDMU(fullsizeConstLB::Vector{Float64}, fullsizeConstUB::Vector{Float64}, fullsizeConstMatrix::SparseMatrixCSC, thetaIdx::Int)
  # parameters
  # input
  # fullsizeConstLB: the lower bound of each constraint of full-size model, which is a vector
  # fullsizeConstUB: the upper bound of each constraint of full-size model, which is a vector
  # fullsizeConstMatrix: the constraint matrix of full-size model
  # thetaIdx: the index of theta w.r.t variable vector
  # output
  # targetDMU: the coefficients of constraints of target DMU, which is a vector

  # the coefficients of output constraints of target DMU is equal to the rhs of output constraints
  rhs = getRHS(fullsizeConstLB, fullsizeConstUB)
  # the coefficients of input constraints of target DMU is equal to the coefficients of theta
  coeffOfTheta = full(fullsizeConstMatrix[:, thetaIdx])

  targetDMU = rhs-coeffOfTheta
  return targetDMU[:, 1]
end

# finding the index of target DMU w.r.t variable vector by finding the vector from constraint matrix which is same as target DMU
function findTargetDMUIdx(constMatrix::SparseMatrixCSC, targetDMU::Vector{Float64})
  # parameters
  # input
  # constMatrix: the constraint matrix of model
  # targetDMU: the coefficients of constraints of target DMU, which is a vector
  # output
  # targetDMUIdx: the index of target DMU w.r.t variable vector

  dimension, scale = size(constMatrix)
  targetDMUIdx = 0
  numSame = scale
  colmn = 1
  check = collect(1:scale)
  idxOfSame = collect(1:scale)

  while colmn <= dimension && numSame != 0
    idxOfSame = find(x -> x == targetDMU[colmn], constMatrix[colmn, check])
    check = check[idxOfSame]
    numSame = length(idxOfSame)

    colmn = colmn + 1
  end

  if numSame != 0
    targetDMUIdx = check[1]
  else
    error("couldn't find target DMU in constraint matrix")
  end

  return targetDMUIdx
end

# generating the coefficients of objective function, lower and upper bound of varibles of small model by set, targetDMUIdx and thetaIdx
function generateSmallV(v::Vector{Float64}, lpUB::Int, set::Vector{Int}, targetDMUIdx::Int, thetaIdx::Int)
  # parameters
  # input
  # v: vector of full-size model, which may be coefficients of objective function, lower or upper bound of varibles
  # lpUB: the limit size of LP
  # set: the selected DMUs, which are the indices w.r.t the variable vector
  # targetDMUIdx: the index of target DMU w.r.t variable vector
  # thetaIdx: the index of theta w.r.t variable vector
  # output
  # smallV: vector of small model, which may be coefficients of objective function, lower or upper bound of varibles

  smallV = Array(Float64, lpUB)
  smallV[1:lpUB-2] = v[set]
  smallV[lpUB-1] = v[targetDMUIdx]
  smallV[lpUB] = v[thetaIdx]

  return smallV
end


text:
@testset "$TEST $G" begin

# TODO one commented out test here intermittenlty fails on Travis
g = G(4)
add_edge!(g, 1,2)
add_edge!(g, 1,3)
add_edge!(g, 1,4)
add_edge!(g, 2,1)
add_edge!(g, 2,3)
add_edge!(g, 2,4)
add_edge!(g, 3,1)
add_edge!(g, 3,2)
add_edge!(g, 3,4)
add_edge!(g, 4,1)
add_edge!(g, 4,2)
add_edge!(g, 4,3)

# Creating custom adjacency matrix
distmx = zeros(4,4)

# Populating custom adjacency matrix
distmx[1,2] = 1
distmx[1,3] = 5
distmx[1,4] = 6
distmx[2,1] = 1
distmx[2,3] = 4
distmx[2,4] = 10
distmx[3,1] = 5
distmx[3,2] = 4
distmx[3,4] = 3
distmx[4,1] = 6
distmx[4,2] = 10
distmx[4,3] = 3

# Testing Kruskal's algorithm
mst = minimum_spanning_tree(g, distmx)
vec_mst = Vector{E}()
push!(vec_mst, E(1, 2))
push!(vec_mst, E(3, 4))
push!(vec_mst, E(2, 3))

@test mst == vec_mst

#second test
g = G(8)
add_edge!(g, 2,8)
add_edge!(g, 1,3)
add_edge!(g, 6,8)
add_edge!(g, 1,8)
add_edge!(g, 3,4)
add_edge!(g, 2,4)
add_edge!(g, 2,6)
add_edge!(g, 3,8)
add_edge!(g, 6,5)
add_edge!(g, 2,3)
add_edge!(g, 5,7)
add_edge!(g, 1,7)
add_edge!(g, 4,7)
add_edge!(g, 7,3)
add_edge!(g, 1,5)
add_edge!(g, 5,8)

distmx_sec = zeros(8, 8)
distmx_sec[2,8] = 0.19
distmx_sec[8,2] = 0.19
distmx_sec[1,3] = 0.26
distmx_sec[3,1] = 0.26
distmx_sec[6,8] = 0.28
distmx_sec[8,2] = 0.28
distmx_sec[1,8] = 0.16
distmx_sec[8,1] = 0.16
distmx_sec[3,4] = 0.17
distmx_sec[4,3] = 0.17
distmx_sec[2,4] = 0.29
distmx_sec[4,2] = 0.29
distmx_sec[2,6] = 0.32
distmx_sec[6,2] = 0.32
distmx_sec[3,8] = 0.34
distmx_sec[8,3] = 0.34
distmx_sec[5,6] = 0.35
distmx_sec[6,5] = 0.35
distmx_sec[2,3] = 0.36
distmx_sec[3,2] = 0.36
distmx_sec[5,7] = 0.93
distmx_sec[7,5] = 0.93
distmx_sec[1,7] = 0.58
distmx_sec[7,1] = 0.58
distmx_sec[4,7] = 0.52
distmx_sec[7,4] = 0.52
distmx_sec[3,7] = 0.40
distmx_sec[7,3] = 0.40
distmx_sec[1,5] = 0.38
distmx_sec[5,1] = 0.38
distmx_sec[5,8] = 0.37
distmx_sec[8,5] = 0.37

mst2 = minimum_spanning_tree(g, distmx_sec)
vec2 = Vector{E}()
push!(vec2, E(1, 8))
push!(vec2, E(3, 4))
push!(vec2, E(2, 8))
push!(vec2, E(1, 3))
push!(vec2, E(6, 8))
push!(vec2, E(5, 6))
push!(vec2, E(3, 7))

@test mst2 == vec2

n = 10; m = 20
g = G(n, m)
g = blockdiag(g, g)
d = blockdiag(sparse(rand(n, n)), sparse(rand(n, n)))
d = d + d'
mst = minimum_spanning_tree(g, d)
#@test_skip length(mst) == 2n - 2 #TODO
for e in mst
    @test has_edge(g, e)
    u , v  = src(e), dst(e)
    u , v = u < v ? (u,v) : (v,u)
    @test  1 <= u < v <= 10 || 11 <= u < v <= 20
end

end # testset


text:
#=
  All different in Turing.jl

  Well, I have to try, don't I. :-)

  With enforce_alldiff true (counting number unique values)

  Distributions of variable num_unique (num:0)
  6.00000 =>   10000  (1.000000)

  With enforce_alldiff false:
  Distributions of variable num_unique (num:0)
  4.00000 =>    5029  (0.502900)
  3.00000 =>    2340  (0.234000)
  5.00000 =>    2295  (0.229500)
  2.00000 =>     197  (0.019700)
  6.00000 =>     138  (0.013800)
  1.00000 =>       1  (0.000100)

  Cf ~/blog/all_different.blog
     ~/webppl/all_different.wppl

=#

using Turing
include("jl_utils.jl")

@model function all_different(n=6,base=10,enforce_alldiff=true)
    # digits = tzeros(n)
    # for i in 1:n
    #     digits[i] ~ DiscreteUniform(0,base-1)
    # end
    # This is much faster: 6.5s vs 33.9s!
    digits ~ filldist(DiscreteUniform(0,base-1), n)

    if enforce_alldiff
        for i in 1:n, j in i+1:n
            true ~ Dirac(digits[i] != digits[j])
        end
    end

    # num_duplicates ~ Dirac(sum([
    #                              sum([digits[i] == d for i in 1:n]) > 1
    #                              for d in 0:base-1
    #                             ]))

    # Faster
    num_unique ~ Dirac(length(unique(digits)))

    # This helps
    if enforce_alldiff
        # true ~ Dirac(num_duplicates == 0)
        true ~ Dirac(num_unique == n)
    end
    # println("digits: ", digits, " num_unique: ", num_unique)

end

n = 6
base = 6
model = all_different(n,base,true)

num_chns = 4
# chns = sample(model, Prior(), 1000)
# chns = sample(model, MH(), 10_000)
# chns = sample(model, PG(15), 10_000)
# chns = sample(model, IS(), 10_000)
chns = sample(model, SMC(), 10_000)
# chns = sample(model, SMC(), MCMCThreads(), 10_000, num_chns)
# chns = sample(model, SGLD(), 10_000)

display(chns)
# display(plot(chns))

# show_var_dist_pct(chns, :num_duplicates)
show_var_dist_pct(chns, :num_unique)


text:
export satellite_plan_milp

function satellite_plan_milp(problem::SatPlanningProblem; sat_id::Integer=1,
    contact_reward::Real=0, timeout::Real=900,
    enforce_first::Bool=false)

    # Initialize MILP problem
    milp = Model(with_optimizer(Gurobi.Optimizer, Presolve=0, TimeLimit=timeout, OutputFlag=0))

    # Initialize Variables
    @variable(milp, x[keys(problem.lt_collects)], Bin)

    # Add Objective
    @objective(milp, Max, sum(opp.reward*x[id] for (id, opp) in problem.lt_collects))

    # Enforce constraint that first opportunity in time must be taken to
    # ensure consistency with other planning methods
    if enforce_first == true
        @constraint(milp, x[1] == 1)
    end

    # Define non-repeat constraints if allowed
    if problem.solve_allow_repeats == false
        # Group opportunities by image
        for req in problem.requests
            # Get all collects for that request
            collects = collect(filter(x -> x.location.id == req.id, problem.opportunities))

            if length(collects) > 1
                @constraint(milp, sum(x[col.id] for col in collects) <= 1)
            end
        end
    end

    for (sopp_id, start_opp) in problem.lt_collects
        for (eopp_id, end_opp) in problem.lt_collects

            if (sopp_id == eopp_id) || (start_opp.t_start >= end_opp.t_end)
               # Conditions to skip insertion
            else
                # Check to see if each transition constraint is valid
                for constraint in problem.constraints
                    if constraint(start_opp, end_opp) == false
                        @constraint(milp, x[sopp_id] + x[eopp_id] <= 1)
                        break # Exit inner loop early
                    end
                end
            end
        end
    end

    # Solve problem
    optimize!(milp)

    # Extract plan and reward from MILP
    plan = Opportunity[]
    reward = 0.0

    for idx in keys(problem.lt_collects)
        if JuMP.value(x[idx]) != 0.0
            push!(plan, problem.lt_collects[idx])
        end
    end

    # Sort plan in ascending order to be viable
    sort!(plan, by = x -> x.t_start)

    reward = objective_value(milp)

    return plan, reward
end