text:
module TestBootsampleParametric

using Bootstrap
using Bootstrap.Datasets
using FactCheck

using DataFrames
using StatsBase
using GLM

function test_bootsample(bs, ref, raw_data, n)

    show(IOBuffer(), bs)
    @fact issubtype(typeof(bs), ParametricBootstrapSample) --> true
    t0 = original(bs)
    @fact length(t0) --> length(ref)
    [@fact t --> roughly(r) for (t, r) in zip(t0, ref)]

    t1 = straps(bs)
    @fact length(t1) --> length(t0)
    [@fact length(t) --> n for t in t1]
    [@fact (minimum(t) <= tr && maximum(t) >= tr) --> true for (t, tr) in zip(t1, t0)]
    [@fact eltype(t) --> eltype(tr) for (t, tr) in zip(t1, t0)]

    @fact Bootstrap.data(bs) --> raw_data ## TODO: define scope

    @fact nrun(sampling(bs)) --> n
    @fact nrun(bs) --> n

    @fact length(bias(bs)) --> length(ref)
    [@fact (b[1] > -Inf && b[1] < Inf) --> true for b in bias(bs)]
    @fact length(se(bs)) --> length(ref)
    [@fact s --> greater_than(0) for s in se(bs)]

    [@fact original(bs, i) --> original(bs)[i]  for i in 1:nvar(bs)]
    [@fact straps(bs, i) --> straps(bs)[i]  for i in 1:nvar(bs)]
    [@fact bias(bs, i) --> bias(bs)[i]  for i in 1:nvar(bs)]
    [@fact se(bs, i) --> se(bs)[i]  for i in 1:nvar(bs)]

    m = model(bs)
    @fact issubtype(typeof(m), Bootstrap.Model) --> true

    return Void
end

function test_ci(bs)

    cim_all = (BasicConfInt(), PercentileConfInt(), NormalConfInt(), BCaConfInt())
    for cim in cim_all
        c = ci(bs, cim)
        [@fact (x[1] >= x[2]  && x[1] <= x[3]) --> true for x in c]
        [@fact x[1] --> roughly(t0) for (x, t0) in zip(c, original(bs))]
    end

    return Void
end

n = 100

## log cannot be specified in formula, see GLM issues
city2 = DataFrame(twenty = log10(city[:,:U]),
                  thirty = log10(city[:,:X]))

facts("Parametric distributions") do

    r = randn(20)
    ref = mean(r)

    context("Basic resampling: Normal distribution") do
        bs = bootstrap(r, mean, Model(Normal), BasicSampling(n))
        test_bootsample(bs, ref, r, n)
        test_ci(bs)
    end

    context("Balanced resampling: Normal distribution") do
        bs = bootstrap(r, mean, Model(Normal), BalancedSampling(n))
        test_bootsample(bs, ref, r, n)
        test_ci(bs)
    end

    ref = mean(fit(Exponential, aircondit))

    context("Basic resampling: Exponential distribution") do
        bs = bootstrap(aircondit, mean, Model(Exponential), BasicSampling(n))
        test_bootsample(bs, ref, aircondit, n)
        test_ci(bs)
    end

    context("Balanced resampling: Exponential distribution") do
        bs = bootstrap(aircondit, mean, Model(Exponential), BalancedSampling(n))
        test_bootsample(bs, ref, aircondit, n)
        test_ci(bs)
    end

end

facts("Linear regression models") do

    ref = coef(fit(LinearModel, @formula(thirty ~ twenty), city2))

    context("Residual resampling") do
        bs = bootstrap(city2, coef, Model(LinearModel, @formula(thirty ~ twenty)), ResidualSampling(n))
        test_bootsample(bs, ref, city2, n)
    end

    context("Wild resampling: Rademacher") do
        bs = bootstrap(city2, coef, Model(LinearModel, @formula(thirty ~ twenty)), WildSampling(n, rademacher))
        test_bootsample(bs, ref, city2, n)
        #@fact typeof(noise(sampling(bs))) --> Function
    end

    context("Wild resampling: Mammen") do
        bs = bootstrap(city2, coef, Model(LinearModel, @formula(thirty ~ twenty)), WildSampling(n, mammen))
        test_bootsample(bs, ref, city2, n)
        #@fact typeof(noise(sampling(bs))) --> Function
    end

end

facts("Generalized linear regression models") do

    ref = coef(fit(GeneralizedLinearModel, @formula(thirty ~ twenty), city2, Normal()))
    max_iter = 200
    conv_tol = 1e-3

    context("Residual resampling") do
        bs = bootstrap(city2, coef,
                       Model(GeneralizedLinearModel, @formula(thirty ~ twenty), Normal(), maxIter = max_iter, convTol = conv_tol),
                       ResidualSampling(n))
        test_bootsample(bs, ref, city2, n)
    end

    context("Residual resampling with link function") do
        bs = bootstrap(city2, coef,
                       Model(GeneralizedLinearModel, @formula(thirty ~ twenty), Normal(), IdentityLink(), maxIter = max_iter, convTol = conv_tol),
                       ResidualSampling(n))
        test_bootsample(bs, ref, city2, n)
    end

    context("Wild resampling: Rademacher") do
        bs = bootstrap(city2, coef,
                       Model(GeneralizedLinearModel, @formula(thirty ~ twenty), Normal(), maxIter = max_iter, convTol = conv_tol),
                       WildSampling(n, rademacher))
        test_bootsample(bs, ref, city2, n)
        @pending typeof(noise(sampling(bs))) --> Function
    end

    context("Wild resampling with link function: Mammen") do
        bs = bootstrap(city2, coef,
                       Model(GeneralizedLinearModel, @formula(thirty ~ twenty), Normal(), IdentityLink(), maxIter = max_iter, convTol = conv_tol),
                       WildSampling(n, mammen))
        test_bootsample(bs, ref, city2, n)
        @pending typeof(noise(sampling(bs))) --> Function
    end

end

end


text:
using JuMP
using Ipopt
using Polynomials, MorePolynomials

struct BoundaryConstraint # make a default setting if unbounded
    initialState::Array # make sure these are Nx1 size matricies
    finalState::Array # incorperate lower bound and ub
end

mutable struct TrajProblem
    stateVectorGuess::AbstractPolynomial
    controlVectorGuess::AbstractPolynomial
    stateVectorDiff::AbstractPolynomial
    diffMat::AbstractArray # differentiation matrix for state
    boundaryConstraints::BoundaryConstraint # add checks that these are dimensionally correct
    τ::AbstractVector
    numStates::Int   # can determine these from size of guess
    numControls::Int
    weights::AbstractVector
end

# work around to make current problem globally available
mutable struct CurrentProblem
    nullableProblem::Union{TrajProblem,Nothing}
end

const CURRENT_PROBLEM = CurrentProblem(nothing)
setCurrentProblem(problem::TrajProblem) = (CURRENT_PROBLEM.nullableProblem = problem)
getCurrentProblem() = CURRENT_PROBLEM.nullableProblem

function solve(problem::TrajProblem) # remove 2.5s to get this to work, and work out why it doesn't work with transformation
    setCurrentProblem(problem)
    model = Model(with_optimizer(Ipopt.Optimizer, max_cpu_time=60.0)) # remove cpu time at some point
    stateLen = length(problem.stateVectorGuess) # number of points on each state vector
    controlLen = length(problem.controlVectorGuess)
    numStates = problem.numStates # number of state vectors
    numControls = problem.numControls
    # assign state and control varibles
    @variable(model, x[1:numStates, 1:stateLen])
    @variable(model, u[1:numControls, 1:controlLen])
    # set state and control variables to problem guess
    [set_start_value(x[i,j], problem.stateVectorGuess[j][2]) for i = 1:numStates, j = 1:stateLen]
    [set_start_value(u[i,j], problem.controlVectorGuess[j][2]) for i = 1:numControls, j = 1:stateLen-1]

    #register custom defined functions
    register(model, :stateDerUpdate!, numStates * stateLen + 1, stateDerUpdate!, ΔstateDerUpdate!)

    # state boundary constraints
    # initial state
    [fix(x[i,1], problem.boundaryConstraints.initialState[i]) for i = 1:length(problem.boundaryConstraints.initialState)]   # adjust this for lb and ub
    # final state
    [fix(x[i,end], problem.boundaryConstraints.finalState[i]) for i = 1:length(problem.boundaryConstraints.finalState)]   # adjust this for lb and ub

    @NLexpression(model, dynamics[i=1:stateLen-1], 2*x[i] + 2*u[i]*sqrt(abs(x[i])))
    @NLexpression(model, objecFun[i=1:stateLen-1], x[i] + u[i]^2)

    for j = 1:stateLen - 1
        @NLconstraint(model, stateDerUpdate!(j,x...) - 2.5*dynamics[j] == 0)
    end
    @NLconstraint(model, x[1] + 2.5*sum(problem.weights[i] * dynamics[i] for i = 1:stateLen-1) - x[stateLen] == 0)

    @NLobjective(model, Min, 2.5*sum(problem.weights[i] * objecFun[i] for i = 1:stateLen-1))
    optimize!(model)
    return value.(x), value.(u)
end

function stateDerUpdate!(ζr, stateVector...) # make sure timestep vector matches length of state vector
    problem = getCurrentProblem()
    stateArray = zeros(typeof(problem.stateVectorGuess[1][1]), length(problem.stateVectorGuess))
    [stateArray[i] = stateVector[i] for i = 1:length(stateVector)]
    ζr = convert(Int, ζr)
    τ = problem.τ
    # assemble state and control vector from tuple inputs
    if stateArray == problem.stateVectorGuess.(τ)
        return problem.stateVectorDiff(τ[ζr])
    else
        problem.stateVectorGuess[:] = stateArray
        problem.stateVectorDiff[:] = problem.diffMat*stateArray
        return problem.stateVectorDiff(τ[ζr])
    end

end

function ΔstateDerUpdate!(g, ζr,stateVector...)
    problem = getCurrentProblem()
    stateArray = zeros(typeof(problem.stateVectorGuess[1][1]), length(problem.stateVectorGuess))
    [stateArray[i] = stateVector[i] for i = 1:length(stateVector)]
    ζr = convert(Int, ζr)
    τ = problem.τ
    g[1] = 0 # gradient of xr and xc is zero as they're just indexes and will not be changed
    g[2:end] = problem.diffMat[ζr,:]
    return g
end

####################################################################################################################################
### example code  ##################################################################################################################
####################################################################################################################################
numCollocationPoints = 15

lgrPoints = lgr_points(numCollocationPoints)
yploc = LagrangePoly([lgrPoints...,1],ones(numCollocationPoints+1))
uploc = LGRPoly(ones(numCollocationPoints))
boundaryConstraints = BoundaryConstraint([2],[1])
problem = TrajProblem(yploc, uploc, derivative(yploc),derivmatrix(yploc), boundaryConstraints,[lgrPoints...,1],1,1,lgr_weights(uploc))
x,u = solve(problem)

using Plots
plotly()
transform(x) = 2.5.*x .+ 2.5
p1 = plot(transform([lgrPoints...,1]),x[1:end],labels=["" ""], ylabel="x")
p2 = plot(transform(lgrPoints),u[1:end],labels=["" ""], ylabel="u", xlabel="t")
plot(p1,p2,layout = (2,1))

ext: jl

max_stars_repo_name: sprockmonty/JuliaFYP

max_issues_repo_name: sprockmonty/JuliaFYP

max_forks_repo_name: sprockmonty/JuliaFYP


text:
using HTTP
using ..Fed: Configuration, serialize_payload, deserialize_payload

"""
    Node(host, port, fit, config)

Orchestrate the training on local data of the clients participating in the
federated training. A `Node` is registerating itself to the server on startup,
and then wait for server requests for training on its local data.

The `fit` function is called "proxied function". It should define the training
loop and respect the following signature:

    fit(weights::Vector{Float32})::Vector{Float32}

Where `weights` is a flatten 1D array contianing the model parameters.
"""
struct Node
    host::String
    port::Int
    fit::Function
    config::Configuration
end

"""
    register_to_server(node::Node)

Register the `node` to the server, letting it knows that a node is available to
take part in the training.

!!! info

    **API Client**: this is an HTTP client request.
"""
function register_to_server(node::Node)::HTTP.Response
    endpoint = node.config.base.serverurl * node.config.base.register_node
    payload = "http://$(node.host):$(node.port)"
    response = HTTP.request("POST", endpoint, [], payload)
    return response
end

"""
    fit_service(node::Node, request::HTTP.Request)

Fit the received weights contained in the `request` with the node's local data,
and return the updated weights to the caller. This function is calling the
"proxied function" `fit` from the `Node`.

!!! info

    **API Endpoint**: this is an HTTP client endpoint.
"""
function fit_service(node::Node, request::HTTP.Request)::HTTP.Response
    weights = deserialize_payload(node.config.payload_serde, request.body, "http://127.0.0.1:9090")
    weights = node.fit(weights)
    payload = serialize_payload(node.config.payload_serde, weights)
    return HTTP.Response(200, payload)
end


text:
using GtkReactive, ImageCore, ImageView
using ImageView: sliceinds
using Test
import AxisArrays
using AxisArrays: Axis

@testset "CLim" begin
    cl = Signal(CLim(0.0, 0.8))
    push!(cl, CLim(0, 1))
    Reactive.run_till_now()
    @test value(cl) === CLim{Float64}(0, 1)

    cmin = RGB(0.2,0.4,0.1)
    cmax = RGB(1.0,0.8,0.95)
    smm = scaleminmax(RGB{N0f8}, cmin, cmax)
    @test @inferred(smm(cmin)) == RGB(0,0,0)
    @test @inferred(smm(cmax)) == RGB(1,1,1)
    @test @inferred(smm((cmin+1.00001*cmax)/2)) == RGB{N0f8}(0.5,0.5,0.5)
end

@testset "NaN and Inf" begin
    # Grayscale
    A = collect(reshape(range(-0.1, stop=1.1, length=25), (5, 5)))
    A[2,2] = NaN
    A[3,3] = -Inf
    A[4,4] = Inf

    target = Gray{N0f8}[
        0   0.0625  0.375   0.6875  1;
        0   0       0.4375  0.75    1;
        0   0.1875  0       0.8125  1;
        0   0.25    0.5625  1       1;
        0   0.3125  0.625   0.9375  1
    ]

    img = Signal(A)
    clim = Signal(CLim(0.1, 0.9))
    enabled, histsigs, imgcsig = ImageView.prep_contrast(img, clim)
    imgc = value(imgcsig)
    @test eltype(imgc) == Gray{N0f8}
    @test imgc == target

    img = Signal(Gray.(A))
    clim = Signal(CLim(0.1, 0.9))
    enabled, histsigs, imgcsig = ImageView.prep_contrast(img, clim)
    imgc = value(imgcsig)
    @test eltype(imgc) == Gray{N0f8}
    @test imgc == target

    # RGB
    Ar = collect(reshape(range(-0.1, stop=1.1, length=3*25), (3, 5, 5)))
    Arc = copy(Ar)
    Ar[1,2,2] = NaN
    Ar[2,3,3] = -Inf
    Ar[3,4,4] = Inf
    target = colorview(RGB, N0f8.((clamp.(Arc, 0.1, 0.9) .- 0.1)./0.8))
    c = target[2,2]
    target[2,2] = RGB{N0f8}(0, green(c), blue(c))
    c = target[3,3]
    target[3,3] = RGB{N0f8}(red(c), 0, blue(c))
    c = target[4,4]
    target[4,4] = RGB{N0f8}(red(c), green(c), 1)

    img = Signal(colorview(RGB, Ar))
    cmin, cmax = RGB(0.1,0.1,0.1), RGB(0.9, 0.9, 0.9)
    clim = Signal(CLim(cmin, cmax))
    enabled, histsigs, imgcsig = ImageView.prep_contrast(img, clim)
    imgc = value(imgcsig)
    @test eltype(imgc) == RGB{N0f8}
    @test imgc == target
end

@testset "sliceinds" begin
    A = rand(3,3,3)
    @test @inferred(sliceinds(A, (1:3, 1:3), Axis{3}(1))) === (1:3, 1:3, 1)
    @test @inferred(sliceinds(A, (1:3, 1:3), Axis{2}(1))) === (1:3, 1, 1:3)
    @test @inferred(sliceinds(A, (1:3, 1:3), Axis{1}(1))) === (1, 1:3, 1:3)
    @test_throws TypeError sliceinds(A, (1:3, 1:3), Axis{1}(1), Axis{2}(1))
    @test_throws BoundsError sliceinds(A, (1:3, 1:3), Axis{3}(1), Axis{2}(1))
    @test_throws BoundsError sliceinds(A, (1:3, 1:3), Axis{2}(1), Axis{2}(1))

    B = AxisArray(A, :X, :Y, :Z)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:Z}(1))) === (1:3, 1:3, 1)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:Y}(1))) === (1:3, 1, 1:3)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:X}(1))) === (1, 1:3, 1:3)
    @test_throws MethodError sliceinds(B, (1:3, 1:3), Axis{:X}(1), Axis{:Y}(1))
    @test_throws MethodError sliceinds(B, (1:3, 1:3), Axis{:Z}(1), Axis{:Y}(1))
    @test_throws MethodError sliceinds(B, (1:3, 1:3), Axis{:Y}(1), Axis{:Y}(1))

    A = rand(3,3,3,3)
    zr = ZoomRegion(axes(A)[1:2])
    @test_throws TypeError sliceinds(A, (1:3, 1:3), Axis{3}(1))
    @test @inferred(sliceinds(A, (1:3, 1:3), Axis{3}(1), Axis{4}(1))) === (1:3, 1:3, 1, 1)
    @test @inferred(sliceinds(A, (1:3, 1:3), Axis{2}(1), Axis{4}(1))) === (1:3, 1, 1:3, 1)
    @test @inferred(sliceinds(A, (1:3, 1:3), Axis{1}(1), Axis{3}(1))) === (1, 1:3, 1, 1:3)
    @test_throws BoundsError sliceinds(A, (1:3, 1:3), Axis{2}(1), Axis{2}(1))

    B = AxisArray(A, :W, :X, :Y, :Z)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:Y}(1), Axis{:Z}(1))) === (1:3, 1:3, 1, 1)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:Z}(1), Axis{:Y}(1))) === (1:3, 1:3, 1, 1)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:X}(1), Axis{:Z}(1))) === (1:3, 1, 1:3, 1)
    @test @inferred(sliceinds(B, (1:3, 1:3), Axis{:X}(1), Axis{:W}(1))) === (1, 1, 1:3, 1:3)
    @test_throws BoundsError sliceinds(B, (1:3, 1:3), Axis{:Y}(1), Axis{:Y}(1))
end

@testset "SliceData" begin
    A = reshape(1:9, 3, 3)
    zr, sd = roi(axes(A), (1,2))
    v = @inferred(slice2d(A, value(zr), sd))
    @test v == A
    zr, sd = roi(axes(A), (2,1))
    v = @inferred(slice2d(A, value(zr), sd))
    @test v == A'

    A = reshape(1:27, 3, 3, 3)
    for (slicedims, target) in (((1, 2), view(A, :, :, 1)),
                                ((1, 3), view(A, :, 1, :)),
                                ((2, 3), view(A, 1, :, :)),
                                ((2, 1), view(A, :, :, 1)'),
                                ((3, 1), view(A, :, 1, :)'),
                                ((3, 2), view(A, 1, :, :)'))
        zr, sd = roi(axes(A), slicedims)
        v = @inferred(slice2d(A, value(zr), sd))
        @test v == target
    end

    A = reshape(1:81, 3, 3, 3, 3)
    for (slicedims, target) in (((1, 2), view(A, :, :, 1, 1)),
                                ((1, 3), view(A, :, 1, :, 1)),
                                ((2, 3), view(A, 1, :, :, 1)),
                                ((2, 4), view(A, 1, :, 1, :)),
                                ((3, 4), view(A, 1, 1, :, :)),
                                ((1, 4), view(A, :, 1, 1, :)),
                                ((2, 1), view(A, :, :, 1, 1)'),
                                ((3, 1), view(A, :, 1, :, 1)'),
                                ((3, 2), view(A, 1, :, :, 1)'),
                                ((4, 2), view(A, 1, :, 1, :)'),
                                ((4, 3), view(A, 1, 1, :, :)'),
                                ((4, 1), view(A, :, 1, 1, :)'))
        zr, sd = roi(axes(A), slicedims)
        v = @inferred(slice2d(A, value(zr), sd))
        @test v == target
    end

    B = AxisArray(A, :x, :y, :z, :time)
    for (slicedims, target) in (((:x, :y), view(B, :, :, 1, 1)),
                                ((:x, :z), view(B, :, 1, :, 1)),
                                ((:y, :z), view(B, 1, :, :, 1)),
                                ((:y, :time), view(B, 1, :, 1, :)),
                                ((:z, :time), view(B, 1, 1, :, :)),
                                ((:x, :time), view(B, :, 1, 1, :)),
                                ((:y, :x), view(B, :, :, 1, 1)'),
                                ((:z, :x), view(B, :, 1, :, 1)'),
                                ((:z, :y), view(B, 1, :, :, 1)'),
                                ((:time, :y), view(B, 1, :, 1, :)'),
                                ((:time, :z), view(B, 1, 1, :, :)'),
                                ((:time, :x), view(B, :, 1, 1, :)'))
        zr, sd = roi(AxisArrays.axes(B), slicedims)
        v = @inferred(slice2d(B, value(zr), sd))
        @test v == target
    end
end

@testset "Canvas size" begin
    @test ImageView.canvas_size((1000,1000), (300, 200)) == (300, 200)
    @test ImageView.canvas_size((1000,1000), (8, 5)) == (160, 100)
    @test ImageView.canvas_size((1000,1000), (3000, 2000)) == (1000, 667)
end

nothing


text:
using Plots, Test

const PLOTS_DEFAULTS = Dict(:theme => :wong2, :fontfamily => :palantino)
Plots.__init__()

@testset "Loading theme" begin
    pl = plot(1:5)
    @test pl[1][1][:seriescolor] == RGBA(colorant"black")
    @test guidefont(pl[1][:xaxis]).family == "palantino"
end

empty!(PLOTS_DEFAULTS)
Plots.__init__()


text:
# Line plots

function PlotLine(label_id::String, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, y::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer=length(x), offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    ImPlot_PlotLine(label_id, x, y, count, offset, stride)
end

function PlotLine(label_id::String, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, y::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer=length(x), offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    ImPlot_PlotLine(label_id, Float64.(x), Float64.(y), count, offset, stride)
end

function PlotLine(label_id::String, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer=length(x), xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0, stride::Integer = sizeof(T)) where {T<:ImPlotData}
    ImPlot_PlotLine(label_id, x, count, xscale, x0, offset, stride)
end

function PlotLine(label_id::String, x::Union{AbstractArray{T},Ref{T},Ptr{T}}, count::Integer=length(x), xscale::Real = 1.0, x0::Real = 0.0, offset::Integer = 0, stride::Integer = sizeof(Float64)) where {T<:Real}
    ImPlot_PlotLine(label_id, Float64.(x), count, xscale, x0, offset, stride)
end

function PlotLine(label_id::String, x::AbstractArray{T}, y::AbstractArray{T};
                  count::Integer = min(length(x), length(y)), offset::Integer = 0,
                  stride::Integer = 1) where {T<:ImPlotData}
    ImPlot_PlotLine(label_id, x, y, count, offset, stride * sizeof(T))
end

function PlotLine(label_id::String, x::AbstractArray{T1}, y::AbstractArray{T2}; kwargs...) where {T1<:Real, T2<:Real}
    PlotLine(label_id, promote(x,y)...; kwargs...)
end

function PlotLine(label_id::String, y::AbstractArray{T}; count::Integer=length(y),
                  xscale::Real = 1.0, x0::Real = 0.0, offset::Integer=0,
                  stride::Integer=1) where {T<:ImPlotData}
    ImPlot_PlotLine(label_id, y, count, xscale, x0, offset, stride * sizeof(T))
end

function PlotLine(label_id::String, x::UnitRange{<:Integer}, y::AbstractArray{T}; xscale::Real = 1.0,
                  x0::Real = 0.0) where {T<:ImPlotData}

    count::Cint = length(x) <= length(y) ? length(x) : throw("Range out of bounds")
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride::Cint = sizeof(T)
    ImPlot_PlotLine(label_id, y, count,  xscale, x0, offset, stride)
end

function PlotLine(label_id::String, x::StepRange, y::AbstractArray{T}; xscale::Real = 1.0, x0::Real = 0.0) where {T<:ImPlotData}

    x.stop < 1 && throw("Range out of bounds")
    count::Cint = length(x) <= length(y) ? length(x) : throw("Range out of bounds")
    offset::Cint = x.start >= 1 ? x.start - 1 : throw("Range out of bounds")
    stride = Cint(x.step * sizeof(T))
    ImPlot_PlotLine(label_id, y, count, xscale, x0, offset, stride)
end

PlotLine(label_id::String, x::UnitRange{<:Integer}, y::AbstractArray{<:Real}; kwargs...) = PlotLine(label_id, x, Float64.(y); kwargs...)
PlotLine(label_id::String, x::StepRange, y::AbstractArray{<:Real}; kwargs...) = PlotLine(label_id, x, Float64.(y); kwargs...)

# xfield, yfield should be propertynames of eltype(structvec)
function PlotLine(
    label_id::String,
    structvec::Vector{T}, xfield::Symbol, yfield::Symbol;
    count::Integer = length(structvec), offset::Integer = 0,
    stride::Integer = 1
) where T

    Tx = fieldtype(T, xfield)
    Ty = fieldtype(T, yfield)
    x_offset = fieldoffset(T, Base.fieldindex(T, xfield))
    y_offset = fieldoffset(T, Base.fieldindex(T, yfield))
    x_ptr = (pointer(structvec, 1) + x_offset) |> Ptr{Tx}
    y_ptr = (pointer(structvec, 1) + y_offset) |> Ptr{Ty}

    if !T.mutable
        # this is somewhat illegal and is used only to pass a pointer through AbstractArray argument into ccall
        x = unsafe_wrap(Vector{Tx}, x_ptr, size(structvec); own = false)
        y = unsafe_wrap(Vector{Ty}, y_ptr, size(structvec); own = false)
        stride = stride * sizeof(T)
    else # two new vectors every 1/60 second...
        x = Vector{Tx}(undef, length(structvec))
        y = Vector{Ty}(undef, length(structvec))
        for (i, val) in enumerate(structvec)
            x[i] = getproperty(val, xfield)
            y[i] = getproperty(val, yfield)
        end
        if Tx !== Ty
            x, y = promote(x, y)
        end
        stride = stride * sizeof(eltype(x))
    end

    ImPlot_PlotLine(label_id, x, y, count, offset, stride)
end

PlotLineG(label_id::String, getter, data, count, offset = 0) =
ImPlot_PlotLineG(label_id, getter, data, count, offset)

line_functions = [
    (Cfloat, Cfloat) => ImPlot_PlotLine_FloatPtrFloatPtr
    (Cdouble, Cdouble) => ImPlot_PlotLine_doublePtrdoublePtr
    (Int32, Int32) => ImPlot_PlotLine_S32PtrS32Ptr
    (UInt64, UInt64) => ImPlot_PlotLine_U64PtrU64Ptr
    (UInt16, UInt16) => ImPlot_PlotLine_U16PtrU16Ptr
    (Int64, Int64) => ImPlot_PlotLine_S64PtrS64Ptr
    (UInt8, UInt8) => ImPlot_PlotLine_U8PtrU8Ptr
    (Int16, Int16) => ImPlot_PlotLine_S16PtrS16Ptr
    (UInt32, UInt32) => ImPlot_PlotLine_U32PtrU32Ptr
    (Int8, Int8) => ImPlot_PlotLine_S8PtrS8Ptr
]
for ((T1, T2), func) in line_functions
    @eval ImPlot_PlotLine(label_id::String, xs::AbstractArray{<:$T1}, ys::AbstractArray{<:$T2}, count, offset, stride) = $func(label_id, xs, ys, count, offset, stride)
end

line_functions_nox = [
    UInt16 => ImPlot_PlotLine_U16PtrInt
    Cfloat => ImPlot_PlotLine_FloatPtrInt
    Int64 => ImPlot_PlotLine_S64PtrInt
    UInt8 => ImPlot_PlotLine_U8PtrInt
    Int16 => ImPlot_PlotLine_S16PtrInt
    UInt32 => ImPlot_PlotLine_U32PtrInt
    Int8 => ImPlot_PlotLine_S8PtrInt
    Cdouble => ImPlot_PlotLine_doublePtrInt
    Int32 => ImPlot_PlotLine_S32PtrInt
    UInt64 => ImPlot_PlotLine_U64PtrInt
]
for (T1, func) in line_functions_nox
    @eval ImPlot_PlotLine(label_id::String, values::AbstractArray{<:$T1}, count, xscale=1.0, x0=0, offset=0, stride=$(sizeof(T1))) = $func(label_id, values, count, xscale, x0, offset, stride)
end


text:
# Build the parameters type and p₀
t = empty_parameter_table()    # initialize table of parameters
add_parameter!(t, :xgeo, 2.17u"mmol/m^3",
    variance_obs = ustrip(upreferred(0.1 * 2.17u"mmol/m^3"))^2,
    description = "Geological mean P concentration",
    LaTeX = "\\state^\\mathrm{geo}")
add_parameter!(t, :τgeo, 1.0u"Myr",
    description = "Geological restoring timescale",
    LaTeX = "\\tau_\\mathrm{geo}")
add_parameter!(t, :k, 10.0u"μmol/m^3",
    optimizable = true,
    description = "Half-saturation constant (Michaelis-Menten)",
    LaTeX = "k_\\vec{u}")
add_parameter!(t, :z₀, 80.0u"m",
    description = "Depth of the euphotic layer base",
    LaTeX = "z_0")
add_parameter!(t, :w₀, 1.0u"m/d",
    optimizable = true,
    description = "Sinking velocity at surface",
    LaTeX = "w_0")
add_parameter!(t, :w′, 1/4.4625u"d",
    optimizable = true,
    description = "Vertical gradient of sinking velocity",
    LaTeX = "w'")
add_parameter!(t, :κDOP, 1/0.25u"yr",
    optimizable = true,
    description = "Remineralization rate constant (DOP to DIP)",
    LaTeX = "\\kappa")
add_parameter!(t, :κPOP, 1/5.25u"d",
    optimizable = true,
    description = "Dissolution rate constant (POP to DOP)",
    LaTeX = "\\kappa")
add_parameter!(t, :σ, 0.3u"1",
    description = "Fraction of quick local uptake recycling",
    LaTeX = "\\sigma")
add_parameter!(t, :τ, 30.0u"d",
    optimizable = true,
    description = "Maximum uptake rate timescale",
    LaTeX = "\\tau_\\vec{u}")
add_parameter!(t, :foo, 1.0u"m")
delete_parameter!(t, :foo)
add_parameter!(t, :foo, 1.0u"m")
delete_parameter!(t, size(t, 1))
initialize_Parameters_type(t)   # Generate the parameter type
p₀ = Parameters()
m_all = length(fieldnames(typeof(p₀)))
m = length(p₀)

@testset "Parameters" begin
    @test t isa DataFrames.DataFrame
    @test_throws ErrorException initialize_Parameters_type(t)
    @test_throws ErrorException add_parameter!(t, :τ, 30.0u"d")
    @test_throws ErrorException delete_parameter!(t, :not_a_parameter)
    @test size(t) == (m_all, 9) # m_all is # of params
    @test size(p₀) == (m,)
    @test length(p₀) == m
    println(p₀)
    @test p₀ isa AIBECS.Parameters{Float64}
    using DualNumbers, HyperDualNumbers
    println(p₀ * im)
    @test p₀ * im isa AIBECS.Parameters{Complex{Float64}}
    println(p₀ * ε)
    @test p₀ * ε isa AIBECS.Parameters{Dual{Float64}}
    println(p₀ * ε₁)
    @test p₀ * ε₁ isa AIBECS.Parameters{Hyper{Float64}}
end


text:
cp = CPEG.CPEGWorkspace()
g = cp.params.gravity
# check gravity at surface
r = g.R*normalize(SA[1, 2, 4.0])

# check gravity at surface is around 3.7 m/s^2
a = CPEG.gravity(g, r)
@test norm(a) < 3.8
@test norm(a) > 3.6

# check J2 calculation by doing both parts seperately

# position vector
r_1 = g.R*1.1*normalize(SA[1, -2, 1.8])

# accel from the gravity function
a1 = CPEG.gravity(g, r_1)

# gravity due to central body
a1_central = -r_1*g.μ/(norm(r_1)^3)

# gravity due to J2
nr = norm(r_1)
x,y,z = r_1
a1_j2 = -1.5*g.J2*(g.μ/nr^2)*(g.R/nr)^2*SA[(1-5*(z/nr)^2)*x/nr,
                                           (1-5*(z/nr)^2)*y/nr,
                                           (3-5*(z/nr)^2)*z/nr]

# test the sum of these is equal to the output of the gravity function
@test a1 ≈ (a1_central + a1_j2) rtol = 1e-12


text:
"""
    get_mesh_surface(points, facets)

Compute surface areas, centers and normals for each facet.
"""
function get_mesh_surface(points, facets)

    # Number of facets
    nfacet = size(facets, 2)

    # Facets
    x = points[:, facets]

    # Facet centers
    centers = reshape(mean(x, dims = 2), 3, nfacet)

    # Facet normals
    normals = reshape(
        mapslices(x -> (x[:, 1] - x[:, 2]) × (x[:, 3] - x[:, 2]), x, dims = [1, 2]),
        3,
        nfacet,
    )

    # Facet areas
    areas = 1 / 2 * reshape(mapslices(norm, normals, dims = 1), nfacet)
    total_area = sum(areas)

    # Normalize normals
    normals = normals ./ 2areas'

    total_area, areas, centers, normals
end


text:
using MarineEcosystemsJuliaCon2021
using Test

@testset "MarineEcosystemsJuliaCon2021.jl" begin
    # Write your tests here.
end


text:
@adjoint function Base.reduce(::typeof(hcat), V::AbstractVector{<:AbstractVector})
    return reduce(hcat, V), dV -> (nothing, collect(eachcol(dV)))
end


text:
using Compose: x_measure, y_measure

function check_arguments(arg::Vector, len)
    if length(arg)>1
        @assert length(arg) == len
    else
        arg = fill(arg[1], len)
    end
    arg
end
check_arguments(arg, len) = fill(arg, len)

struct HLineGeometry <: Gadfly.GeometryElement
    color::Union{Vector, Color, (Nothing)}
    size::Union{Vector, Measure, (Nothing)}
    style::Union{Vector, Symbol, (Nothing)}
    tag::Symbol

    HLineGeometry(color, size, style, tag) =
            new(color === nothing ? nothing : Gadfly.parse_colorant(color), size, style, tag)
end
HLineGeometry(; color=nothing, size=nothing, style=nothing, tag=empty_tag) =
        HLineGeometry(color, size, style, tag)

"""
    Geom.hline[(; color=nothing, size=nothing, style=nothing)]

Draw horizontal lines across the plot canvas at each position in the `yintercept` aesthetic.
"""
const hline = HLineGeometry

element_aesthetics(::HLineGeometry) = [:yintercept]

# Generate a form for the hline geometry
function render(geom::HLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.hline", aes, :yintercept)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size
    style = geom.style === nothing ? theme.line_style[1] : geom.style

    color = check_arguments(color, length(aes.yintercept))
    size = check_arguments(size, length(aes.yintercept))
    style = check_arguments(style, length(aes.yintercept))

    style = map(Gadfly.get_stroke_vector, style)

    root = compose(context(), svgclass("xfixed"))
    for (idx,y) in enumerate(aes.yintercept)
        compose!(root, (context(),
            Compose.line([(0w, y_measure(y)), (1w, y_measure(y))], geom.tag),
            stroke(color[idx]),
            linewidth(size[idx]),
            strokedash(style[idx])))
    end
    root
end

struct VLineGeometry <: Gadfly.GeometryElement
    color::Union{Vector, Color, (Nothing)}
    size::Union{Vector, Measure, (Nothing)}
    style::Union{Vector, Symbol, (Nothing)}
    tag::Symbol

    VLineGeometry(color, size, style, tag) =
            new(color === nothing ? nothing : Gadfly.parse_colorant(color), size, style, tag)
end
VLineGeometry(; color=nothing, size=nothing, style=nothing, tag=empty_tag) =
        VLineGeometry(color, size, style, tag)

"""
    Geom.vline[(; color=nothing, size=nothing, style=nothing)]

Draw vertical lines across the plot canvas at each position in the `xintercept` aesthetic.
"""
const vline = VLineGeometry

element_aesthetics(::VLineGeometry) = [:xintercept]

# Generate a form for the vline geometry
function render(geom::VLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    Gadfly.assert_aesthetics_defined("Geom.vline", aes, :xintercept)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size
    style = geom.style === nothing ? theme.line_style[1] : geom.style

    color = check_arguments(color, length(aes.xintercept))
    size = check_arguments(size, length(aes.xintercept))
    style = check_arguments(style, length(aes.xintercept))

    style = map(Gadfly.get_stroke_vector, style)

    root = compose(context(), svgclass("yfixed"))
    for (idx,x) in enumerate(aes.xintercept)
        compose!(root, (context(),
                Compose.line([(x_measure(x), 0h), (x_measure(x), 1h)], geom.tag),
                stroke(color[idx]),
                linewidth(size[idx]),
                strokedash(style[idx])))
    end
    root
end

struct ABLineGeometry <: Gadfly.GeometryElement
    color::Union{Vector, Color, (Nothing)}
    size::Union{Vector, Measure, (Nothing)}
    style::Union{Vector, Symbol, (Nothing)}
    tag::Symbol

    ABLineGeometry(color, size, style, tag) =
            new(color === nothing ? nothing : Gadfly.parse_colorant(color), size, style, tag)
end
ABLineGeometry(; color=nothing, size=nothing, style=nothing, tag::Symbol=empty_tag) =
        ABLineGeometry(color, size, style, tag)

"""
    Geom.abline[(; color=nothing, size=nothing, style=nothing)]

For each corresponding pair of elements in the `intercept` and `slope` aesthetics,
draw the lines `T(y) = slope * T(x) + intercept` across the plot canvas, where `T(⋅)` defaults to the identity function.
If unspecified, `intercept` defaults to [0] and `slope` to [1].

This geometry also works with nonlinear `Scale` transformations of the `y` and/or `x` variable, with one caveat: for log transformations of the `x` variable, the `intercept` is the `y`-value at `x=1` rather than at `x=0`.
"""
const abline = ABLineGeometry

element_aesthetics(geom::ABLineGeometry) = [:intercept, :slope]

function render(geom::ABLineGeometry, theme::Gadfly.Theme, aes::Gadfly.Aesthetics)
    if aes.intercept == nothing && aes.slope == nothing
        aes.intercept = [0]
        aes.slope = [1]
    elseif aes.intercept == nothing
        aes.intercept = fill(0,length(aes.slope))
    elseif aes.slope == nothing
        aes.slope = fill(1,length(aes.intercept))
    end
    Gadfly.assert_aesthetics_equal_length("Geom.line", aes, element_aesthetics(geom)...)

    color = geom.color === nothing ? theme.default_color : geom.color
    size = geom.size === nothing ? theme.line_width : geom.size
    style = geom.style === nothing ? theme.line_style[1] : geom.style

    color = check_arguments(color, length(aes.intercept))
    size = check_arguments(size, length(aes.intercept))
    style = check_arguments(style, length(aes.intercept))

    style = map(Gadfly.get_stroke_vector, style)

    # it would've been nice to just say x0, x1 = realmin(Float64), realmax(Float64).
    # but SVG() and PDF() have silent overflow errors when plotting way outside the
    # context's bounding box.  so instead, use the extrema of the data, and hope that
    # the line extends to the edges of the graph.

    if typeof(aes.y) <: Array{Function}
        lowx, highx = aes.xmin[1], aes.xmax[1]
    else
        lowx, highx = extrema(aes.x)
    end

    # extending the line to width 3 times x-range of data should be enough
    rangex = highx - lowx
    x0 = lowx - rangex
    x1 = highx + rangex

    y0s = [x0 * m + b for (m,b) in zip(aes.slope, aes.intercept)]
    y1s = [x1 * m + b for (m,b) in zip(aes.slope, aes.intercept)]

    root = compose(context(), svgclass("geometry"))
    for (idx,(y0,y1)) in enumerate(zip(y0s,y1s))
        compose!(root, (context(),
                Compose.line([(x_measure(x0),y_measure(y0)),
                              (x_measure(x1),y_measure(y1))], geom.tag),
                stroke(color[idx]),
                linewidth(size[idx]),
                strokedash(style[idx])))
    end
    root
end


text:
preptime(layers) = 2 * layers
remaining_time(current_time) = 60 - current_time
total_working_time(layers, current_time) = preptime(layers) + current_time


text:
"""
Creates Ohms constraints for shiftable PSTs / OLTCs

```
p[f_idx] == -b*((t[f_bus] - t_shift[f_idx]) - (t[t_bus] - t_shift[t_idx]))
```
"""
function constraint_variable_transformer_y_from(pm::GenericPowerModel{T}, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, g_shunt, tap_min, tap_max) where T <: PowerModels.AbstractDCPForm
    p_fr = PowerModels.var(pm, n, cnd, :p, f_idx)
    va_fr = PowerModels.var(pm, n, cnd, :va, f_bus)
    va_to = PowerModels.var(pm, n, cnd, :va, t_bus)
    va_shift_fr = PowerModels.var(pm, n, cnd, :va_shift, f_idx)
    va_shift_to = PowerModels.var(pm, n, cnd, :va_shift, t_idx)

    @constraint(pm.model, p_fr ==  g_shunt/2 + (-b)*((va_fr - va_shift_fr) - (va_to - va_shift_to)))
    # omit reactive constraint
end

"Do nothing, this model is symmetric"
function constraint_variable_transformer_y_to(pm::GenericPowerModel{T}, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, g, b, c, g_shunt, tap_min, tap_max) where T <: PowerModels.AbstractDCPForm
end
"Do nothing, there are no voltage magnitude variables"
function constraint_link_voltage_magnitudes(pm::GenericPowerModel{T}, n::Int, cnd::Int, f_bus, t_bus, f_idx, t_idx, tap_fr, tap_to) where T <: PowerModels.AbstractDCPForm
end

""
function constraint_kcl_shunt_aggregated(pm::GenericPowerModel{T}, n::Int, cnd::Int, i, bus_arcs, bus_arcs_dc, bus_gs, bus_bs) where T <: PowerModels.AbstractDCPForm
    p = PowerModels.var(pm, n, cnd, :p)
    p_dc = PowerModels.var(pm, n, cnd, :p_dc)
    pnode = PowerModels.var(pm, n, cnd, :pnode, i)

    @constraint(pm.model, sum(p[a] for a in bus_arcs) + sum(p_dc[a_dc] for a_dc in bus_arcs_dc) == pnode - sum(gs for gs in values(bus_gs))*1.0^2)
    # omit reactive constraint
end

function constraint_flexible_gen(pm::GenericPowerModel{T}, i::Int;  nw::Int = pm.cnw, cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
        constraint_flexible_active_gen(pm, nw, cnd, i)
end

function constraint_flexible_load(pm::GenericPowerModel{T}, i::Int;  nw::Int = pm.cnw, cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
        constraint_flexible_active_load(pm, nw, cnd,  i)
end

function constraint_fixed_load(pm::GenericPowerModel{T}, i::Int;  nw::Int = pm.cnw, cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
        constraint_fixed_active_load(pm, nw, cnd, i)
end

function constraint_load_gen_aggregation_sheddable(pm::GenericPowerModel{T}, i::Int; nw::Int = pm.cnw, cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    constraint_active_load_gen_aggregation_sheddable(pm, nw, cnd, i)
end

function constraint_load_gen_aggregation(pm::GenericPowerModel{T}, i::Int;  nw::Int = pm.cnw, cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
        constraint_active_load_gen_aggregation(pm, nw, cnd, i)
end

function constraint_redispatch_power_gen(pm::GenericPowerModel{T}, i::Int;  nw::Int = pm.cnw, cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
        constraint_redispatch_active_power_gen(pm, nw, cnd, i)
end

function constraint_gen_contingency(pm::GenericPowerModel{T}, n::Int, i::Int; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    constraint_active_power_gen_contingency(pm, n, cnd, i)
end

function constraint_branch_contingency(pm::GenericPowerModel{T}, n::Int, i::Int; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    constraint_active_power_branch_contingency(pm, n, cnd, i)
end

function constraint_second_stage_redispatch_power_gen(pm::GenericPowerModel{T}, n::Int, i::Int, first_stage_network_id; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    constraint_second_stage_redispatch_active_power_gen(pm, n, cnd, i, first_stage_network_id)
end

function constraint_second_stage_redispatch_power_load(pm::GenericPowerModel{T}, n::Int, i::Int, first_stage_network_id; cnd::Int = pm.ccnd) where T <: PowerModels.AbstractDCPForm
    constraint_second_stage_redispatch_active_power_load(pm, n, cnd, i, first_stage_network_id)
end

function constraint_redispatch_active_power_gen(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, pref::AbstractFloat) where T <: PowerModels.AbstractDCPForm
    pg_delta = PowerModels.var(pm, n, cnd, :pg_delta, i)
    pg = PowerModels.var(pm, n, cnd, :pg, i)

    @constraint(pm.model, pg_delta >= pg - pref)
    @constraint(pm.model, pg_delta >= -(pg - pref))
end
# DO Nothing
function constraint_redispatch_reactive_power_gen(pm::GenericPowerModel{T}, n::Int,  cnd::Int, i::Int, qref::AbstractFloat) where T <: PowerModels.AbstractDCPForm
end

function constraint_second_stage_redispatch_active_power_gen(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, first_stage_network_id) where T <: PowerModels.AbstractDCPForm
    pg_delta = PowerModels.var(pm, n, cnd, :pg_delta, i)
    pg = PowerModels.var(pm, n, cnd, :pg, i)
    pg_first_stage = PowerModels.var(pm, first_stage_network_id, cnd, :pg, i)

    @constraint(pm.model, pg_delta >= pg - pg_first_stage)
    @constraint(pm.model, pg_delta >= -(pg - pg_first_stage))
end

# DO Nothing
function constraint_second_stage_redispatch_reactive_power_gen(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, first_stage_network_id) where T <: PowerModels.AbstractDCPForm
end

function constraint_redispatch_active_power_load(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, pref::AbstractFloat) where T <: PowerModels.AbstractDCPForm
    pl_delta = PowerModels.var(pm, n, cnd, :pl_delta, i)
    pl = PowerModels.var(pm, n, cnd, :pl, i)

    @constraint(pm.model, pl_delta >= pl - pref)
    @constraint(pm.model, pl_delta >= -(pl - pref))
end

# DO NOTING
function constraint_redispatch_reactive_power_load(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, qref::AbstractFloat) where T <: PowerModels.AbstractDCPForm
end

function constraint_second_stage_redispatch_active_power_load(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, first_stage_network_id)  where T <: PowerModels.AbstractDCPForm
    pl_delta = PowerModels.var(pm, n, cnd, :pl_delta, i)
    pl = PowerModels.var(pm, n, cnd, :pl, i)
    pl_first_stage = PowerModels.var(pm, first_stage_network_id, cnd, :pl, i)

    @constraint(pm.model, pl_delta >= pl - pl_first_stage)
    @constraint(pm.model, pl_delta >= -(pl - pl_first_stage))
end

# DO NOTHING
function constraint_second_stage_redispatch_reactive_power_load(pm::GenericPowerModel{T}, n::Int, cnd::Int, i::Int, first_stage_network_id)  where T <: PowerModels.AbstractDCPForm
end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
# Modelo de bomba
type Bomba
	Bomba()=begin
		PP2=outers.PP2
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"Steam tables"
			)),
			Entalpia(),
			Pot_sinal(Dict{Symbol,Any}(
				:Brief=>"Potencia do motor da bomba",
				:PosX=>0.5,
				:PosY=>1
			)),
			Potencia(Dict{Symbol,Any}(
				:Brief=>"Potencia injetada pela bomba"
			)),
			Eficiencia(Dict{Symbol,Any}(
				:Brief=>"Eficiencia da bomba"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{in}"
			)),
			Corrente(Dict{Symbol,Any}(
				:Symbol=>"_{out}",
				:PosX=>0,
				:PosY=>0.2
			)),
			VolumeEspecifico(),
			[
				:(v_esp = PP2.Vesp(Fin.P,Fin.T)),
				:(H_IS = PP2.propPS(Fout.P,Fin.S)),
				:((Fout.H - Fin.H) * EF_B = H_IS - Fin.H),
				:([Fout.S,Fout.T] = PP2.propPH(Fout.P,Fout.H)),
				:(POT_EF = POT_BMB * EF_B),
				:(POT_EF = Fin.F * v_esp * (Fout.P - Fin.P)),
				:(Fout.F = Fin.F),
			],
			[
				"","","","","","","",
			],
			[:PP2,],
			[:H_IS,:POT_BMB,:POT_EF,:EF_B,:Fin,:Fout,:v_esp,]
		)
	end
	PP2::DanaPlugin
	H_IS::Entalpia
	POT_BMB::Pot_sinal
	POT_EF::Potencia
	EF_B::Eficiencia
	Fin::Corrente
	Fout::Corrente
	v_esp::VolumeEspecifico
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Bomba
function setEquationFlow(in::Bomba)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	#	(Fout.H - Fin.H) * Fin.F = POT_EF; # Forma alternativa
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
end
function atributes(in::Bomba,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/bomba1"
	drive!(fields,_)
	return fields
end
Bomba(_::Dict{Symbol,Any})=begin
	newModel=Bomba()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
"Julog terms are variables, constants, or compound terms."
abstract type Term end

"Julog constants."
struct Const <: Term
    name::Any
end

"Julog variables."
struct Var <: Term
    name::Union{Symbol,UInt}
end

"Julog compound terms (e.g. )."
struct Compound <: Term
    name::Symbol
    args::Vector{Term}
end

"Julog clauses are definite Horn clauses of the form [head] <<= [body]."
struct Clause
    head::Term
    body::Vector{Term}
end

"Substitution mapping from variables to terms."
const Subst = Dict{Var,Term}

"Check if two terms are exactly equal."
Base.:(==)(t1::Term, t2::Term) = false
Base.:(==)(t1::Const, t2::Const) = t1.name == t2.name
Base.:(==)(t1::Var, t2::Var) = t1.name == t2.name
Base.:(==)(t1::Compound, t2::Compound) =
    (t1.name == t2.name && length(t1.args) == length(t2.args) &&
            all(a1 == a2 for (a1, a2) in zip(t1.args, t2.args)))

"Compute hash of Julog term from name and arguments."
Base.hash(t::Term, h::UInt) = error("Not implemented.")
Base.hash(t::Const, h::UInt) = hash(t.name, h)
Base.hash(t::Var, h::UInt) = hash(t.name, h)
Base.hash(t::Compound, h::UInt) = hash(t.name, hash(Tuple(t.args), h))

"Check if two clauses are exactly equal."
Base.:(==)(c1::Clause, c2::Clause) =
    (c1.head == c2.head && length(c1.body) == length(c2.body) &&
     all(t1 == t2 for (t1, t2) in zip(c1.body, c2.body)))

"Compute hash of Julog clause from head and body."
Base.hash(c::Clause, h::UInt) = hash(c.head, hash(Tuple(c.body), h))

"Convert Julog term to Horn clause."
Base.convert(::Type{Clause}, term::Term) = Clause(term, [])

"Convert Horn clause to Julog term."
function Base.convert(::Type{Term}, clause::Clause)
    if length(clause.body) == 0 return clause.head end
    return Compound(:(=>), Term[Compound(:and, copy(clause.body)), clause.head])
end

"Show Julog terms as they would be parsed."
function Base.show(io::IO, t::Term)
    print(io, t.name)
end

function Base.show(io::IO, t::Var)
    (t.name isa UInt) ? print(io, "#", t.name) : print(io, t.name)
end

function Base.show(io::IO, t::Compound)
    if t.name == :cons && length(t.args) == 2
        # Handle lists separately
        head, tail = t.args[1], t.args[2]
        if isa(tail, Var)
            print(io, "[", repr(head), " | ", repr(tail), "]")
        elseif isa(tail, Compound) && tail.name == :cend
            print(io, "[", repr(head), "]")
        else
            print(io, "[", repr(head), ", ", repr(tail)[2:end-1], "]")
        end
    elseif t.name == :cend && length(t.args) == 0
        print(io, "[]")
    else
        # Print compound term as "name(args...)"
        print(io, t.name, "(", join([repr(a) for a in t.args], ", ")..., ")")
    end
end

"Show Julog clauses as they would be parsed."
function Base.show(io::IO, c::Clause)
    if length(c.body) == 0
        print(io, c.head)
    else
        print(io, c.head, " <<= ", join([repr(t) for t in c.body], " & ")...)
    end
end

"Show variable substitions."
function Base.show(io::IO, subst::Subst)
    print(io, "{", join(["$k => $v" for (k, v) in subst], ", ")..., "}")
end

"Get arguments of a term."
get_args(term::Term) = Term[]
get_args(term::Compound) = term.args

Base.getproperty(t::Const, f::Symbol) = f == :args ? Term[] : getfield(t, f)
Base.getproperty(t::Var, f::Symbol) = f == :args ? Term[] : getfield(t, f)

Base.propertynames(t::Const) = (:name, :args)
Base.propertynames(t::Var) = (:name, :args)


text:
using CLIMA.Mesh.BrickMesh
using Test
using MPI

MPI.Initialized() || MPI.Init()

@testset "Linear Parition" begin
  @test BrickMesh.linearpartition(1,1,1) == 1:1
  @test BrickMesh.linearpartition(20,1,1) == 1:20
  @test BrickMesh.linearpartition(10,1,2) == 1:5
  @test BrickMesh.linearpartition(10,2,2) == 6:10
end

@testset "Hilbert Code" begin
  @test BrickMesh.hilbertcode([0,0], bits=1) == [0, 0]
  @test BrickMesh.hilbertcode([0,1], bits=1) == [0, 1]
  @test BrickMesh.hilbertcode([1,1], bits=1) == [1, 0]
  @test BrickMesh.hilbertcode([1,0], bits=1) == [1, 1]
  @test BrickMesh.hilbertcode([0,0], bits=2) == [0, 0]
  @test BrickMesh.hilbertcode([1,0], bits=2) == [0, 1]
  @test BrickMesh.hilbertcode([1,1], bits=2) == [0, 2]
  @test BrickMesh.hilbertcode([0,1], bits=2) == [0, 3]
  @test BrickMesh.hilbertcode([0,2], bits=2) == [1, 0]
  @test BrickMesh.hilbertcode([0,3], bits=2) == [1, 1]
  @test BrickMesh.hilbertcode([1,3], bits=2) == [1, 2]
  @test BrickMesh.hilbertcode([1,2], bits=2) == [1, 3]
  @test BrickMesh.hilbertcode([2,2], bits=2) == [2, 0]
  @test BrickMesh.hilbertcode([2,3], bits=2) == [2, 1]
  @test BrickMesh.hilbertcode([3,3], bits=2) == [2, 2]
  @test BrickMesh.hilbertcode([3,2], bits=2) == [2, 3]
  @test BrickMesh.hilbertcode([3,1], bits=2) == [3, 0]
  @test BrickMesh.hilbertcode([2,1], bits=2) == [3, 1]
  @test BrickMesh.hilbertcode([2,0], bits=2) == [3, 2]
  @test BrickMesh.hilbertcode([3,0], bits=2) == [3, 3]

  @test BrickMesh.hilbertcode(UInt64.([14,3,4])) == UInt64.([0x0,0x0,0xe25])
end

@testset "Mesh to Hilbert Code" begin
  let
    etc = Array{Float64}(undef, 2, 4, 6)
    etc[:, :, 1] = [2.0  3.0  2.0  3.0; 4.0  4.0  5.0  5.0]
    etc[:, :, 2] = [3.0  4.0  3.0  4.0; 4.0  4.0  5.0  5.0]
    etc[:, :, 3] = [4.0  5.0  4.0  5.0; 4.0  4.0  5.0  5.0]
    etc[:, :, 4] = [2.0  3.0  2.0  3.0; 5.0  5.0  6.0  6.0]
    etc[:, :, 5] = [3.0  4.0  3.0  4.0; 5.0  5.0  6.0  6.0]
    etc[:, :, 6] = [4.0  5.0  4.0  5.0; 5.0  5.0  6.0  6.0]

    code_exect = UInt64[0x0000000000000000 0x1555555555555555 0xffffffffffffffff 0x5555555555555555 0x6aaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa; 0x0000000000000000 0x5555555555555555 0xffffffffffffffff 0x5555555555555555 0xaaaaaaaaaaaaaaaa 0xaaaaaaaaaaaaaaaa]

    code = centroidtocode(MPI.COMM_SELF, etc)

    @test code == code_exect
  end

  let
    nelem = 1
    d = 2

    etc = Array{Float64}(undef, d, d^2, nelem)
    etc[:, :, 1] = [2.0  3.0  2.0  3.0; 4.0  4.0  5.0  5.0]
    code = centroidtocode(MPI.COMM_SELF, etc)

    @test code == zeros(eltype(code), d, nelem)
  end
end

@testset "Vertex Ordering" begin
  @test ((1,), 1) == BrickMesh.vertsortandorder(1)

  @test ((1,2), 1) == BrickMesh.vertsortandorder(1, 2)
  @test ((1,2), 2) == BrickMesh.vertsortandorder(2, 1)

  @test ((1,2,3), 1) == BrickMesh.vertsortandorder(1, 2, 3)
  @test ((1,2,3), 2) == BrickMesh.vertsortandorder(3, 1, 2)
  @test ((1,2,3), 3) == BrickMesh.vertsortandorder(2, 3, 1)
  @test ((1,2,3), 4) == BrickMesh.vertsortandorder(2, 1, 3)
  @test ((1,2,3), 5) == BrickMesh.vertsortandorder(3, 2, 1)
  @test ((1,2,3), 6) == BrickMesh.vertsortandorder(1, 3, 2)

  err = ErrorException("Problem finding vertex ordering (1, 1, 2) with flips" *
                       "\n          (false, true, false)")
  @test_throws err BrickMesh.vertsortandorder(2, 1, 1)

  @test ((1,2,3,4), 1) == BrickMesh.vertsortandorder(1, 2, 3, 4)
  @test ((1,2,3,4), 2) == BrickMesh.vertsortandorder(1, 3, 2, 4)
  @test ((1,2,3,4), 3) == BrickMesh.vertsortandorder(2, 1, 3, 4)
  @test ((1,2,3,4), 4) == BrickMesh.vertsortandorder(2, 4, 1, 3)
  @test ((1,2,3,4), 5) == BrickMesh.vertsortandorder(3, 1, 4, 2)
  @test ((1,2,3,4), 6) == BrickMesh.vertsortandorder(3, 4, 1, 2)
  @test ((1,2,3,4), 7) == BrickMesh.vertsortandorder(4, 2, 3, 1)
  @test ((1,2,3,4), 8) == BrickMesh.vertsortandorder(4, 3, 2, 1)

  err = ErrorException("Problem finding vertex ordering (1, 1, 3, 3)" *
                       "\n            with flips" *
                       " (false, true, false, false, true)")
  @test_throws err BrickMesh.vertsortandorder(1, 3, 3, 1)
end

@testset "Mesh" begin
  let
    (etv, etc, etb, fc) = brickmesh((4:7,), (false,))
    etv_expect = [1 2 3
                  2 3 4]
    etb_expect = [1 0 0
                  0 0 1]
    fc_expect = Array{Int64,1}[]

    @test etv == etv_expect
    @test etb == etb_expect
    @test fc == fc_expect
    @test etc[:,:,1] == [4 5]
    @test etc[:,:,2] == [5 6]
    @test etc[:,:,3] == [6 7]
  end

  let
    (etv, etc, etb, fc) = brickmesh((4:7,), (true,))
    etv_expect = [1 2 3
                  2 3 4]
    etb_expect = [0 0 0
                  0 0 0]
    fc_expect = Array{Int64,1}[[3, 2, 1]]

    @test etv == etv_expect
    @test etb == etb_expect
    @test fc == fc_expect
    @test etc[:,:,1] == [4 5]
    @test etc[:,:,2] == [5 6]
    @test etc[:,:,3] == [6 7]
  end

  let
    (etv, etc, etb, fc) = brickmesh((2:5,4:6), (false,true))

    etv_expect = [ 1  2  5  6
                   2  3  6  7
                   3  4  7  8
                   5  6  9 10
                   6  7 10 11
                   7  8 11 12]'
    etb_expect = [1 0 0 1 0 0
                  0 0 1 0 0 1
                  0 0 0 0 0 0
                  0 0 0 0 0 0]
    fc_expect = Array{Int64,1}[[4, 4, 1, 2],
                               [5, 4, 2, 3],
                               [6, 4, 3, 4]]

    @test etv == etv_expect
    @test etb == etb_expect
    @test fc == fc_expect
    @test etc[:,:,1] == [2 3 2 3
                         4 4 5 5]
    @test etc[:,:,5] == [3 4 3 4
                         5 5 6 6]
  end

  let
    (etv, etc, etb, fc) = brickmesh((-1:2:1,-1:2:1,-1:1:1), (true,true,true))
    etv_expect = [1   5
                  2   6
                  3   7
                  4   8
                  5   9
                  6  10
                  7  11
                  8  12]
    etb_expect = zeros(Int64, 6, 2)

    fc_expect = Array{Int64,1}[[1, 2, 1, 3, 5,  7],
                               [1, 4, 1, 2, 5,  6],
                               [2, 2, 5, 7, 9, 11],
                               [2, 4, 5, 6, 9, 10],
                               [2, 6, 1, 2, 3,  4]]

    @test etv == etv_expect
    @test etb == etb_expect
    @test fc == fc_expect

    @test etc[:,:,1] == [-1  1 -1  1 -1  1 -1  1
                         -1 -1  1  1 -1 -1  1  1
                         -1 -1 -1 -1  0  0  0  0]

    @test etc[:,:,2] == [-1  1 -1  1 -1  1 -1  1
                         -1 -1  1  1 -1 -1  1  1
                          0  0  0  0  1  1  1  1]
  end

  let
    (etv, etc, etb, fc) = brickmesh((-1:1,-1:1,-1:1),
                                    (false,false,false),
                                    boundary=((11,12),(13,14),(15,16)))

    @test etb == [11  0 11  0 11  0 11  0
                   0 12  0 12  0 12  0 12
                  13 13  0  0 13 13  0  0
                   0  0 14 14  0  0 14 14
                  15 15 15 15  0  0  0  0
                   0  0  0  0 16 16 16 16]
  end

  let
    x = (1:1000,)
    p = (false,)
    b = ((1,2),)

    (etv, etc, etb, fc) = brickmesh(x, p, boundary=b)

    n = 50
    (etv_parts, etc_parts, etb_parts, fc_parts) = brickmesh(x, p, boundary=b,
                                                            part=1,
                                                            numparts=n)
    for j=2:n
      (etv_j, etc_j, etb_j, fc_j) = brickmesh(x, p, boundary=b,
                                              part=j, numparts=n)
      etv_parts = cat(etv_parts, etv_j; dims=2)
      etc_parts = cat(etc_parts, etc_j; dims=3)
      etb_parts = cat(etb_parts, etb_j; dims=2)
    end

    @test etv == etv_parts
    @test etc == etc_parts
    @test etb == etb_parts
  end

  let
    x = (-1:2:10,-1:1:1,-4:1:1)
    p = (true,false,true)
    b = ((1,2),(3,4),(5,6))

    (etv, etc, etb, fc) = brickmesh(x, p, boundary=b)

    n = 50
    (etv_parts, etc_parts, etb_parts, fc_parts) = brickmesh(x, p, boundary=b,
                                                            part=1,
                                                            numparts=n)
    for j=2:n
      (etv_j, etc_j, etb_j, fc_j) = brickmesh(x, p, boundary=b,
                                              part=j, numparts=n)
      etv_parts = cat(etv_parts, etv_j; dims=2)
      etc_parts = cat(etc_parts, etc_j; dims=3)
      etb_parts = cat(etb_parts, etb_j; dims=2)
    end

    @test etv == etv_parts
    @test etc == etc_parts
    @test etb == etb_parts
  end
end

@testset "Connect" begin
  let
    comm = MPI.COMM_SELF

    mesh = connectmesh(comm, partition(comm, brickmesh((0:10,),
                                                       (true,))...)[1:4]...)

    nelem = 10

    @test mesh[:elemtocoord][:,:, 1] == [0  1]
    @test mesh[:elemtocoord][:,:, 2] == [1  2]
    @test mesh[:elemtocoord][:,:, 3] == [2  3]
    @test mesh[:elemtocoord][:,:, 4] == [3  4]
    @test mesh[:elemtocoord][:,:, 5] == [4  5]
    @test mesh[:elemtocoord][:,:, 6] == [5  6]
    @test mesh[:elemtocoord][:,:, 7] == [6  7]
    @test mesh[:elemtocoord][:,:, 8] == [7  8]
    @test mesh[:elemtocoord][:,:, 9] == [8  9]
    @test mesh[:elemtocoord][:,:,10] == [9 10]

    @test mesh[:elemtoelem] == [10 1 2 3 4 5 6 7  8 9
                                 2 3 4 5 6 7 8 9 10 1]

    @test mesh[:elemtoface] == [2 2 2 2 2 2 2 2 2 2
                                1 1 1 1 1 1 1 1 1 1]

    @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))
    @test mesh[:elemtobndy] == zeros(Int, size(mesh[:elemtoordr]))

    @test mesh[:elems] == 1:nelem
    @test mesh[:realelems] == 1:nelem
    @test mesh[:ghostelems] == nelem.+(1:0)

    @test length(mesh[:sendelems]) == 0

    @test mesh[:nabrtorank] == Int[]
    @test mesh[:nabrtorecv] == UnitRange{Int}[]
    @test mesh[:nabrtosend] == UnitRange{Int}[]
  end

  let
    comm = MPI.COMM_SELF
    mesh = connectmesh(comm, partition(comm, brickmesh((0:4,5:9),
                                                       (false,true))...)[1:4]...)

    nelem = 16

    @test mesh[:elemtocoord][:,:, 1] == [0 1 0 1; 5 5 6 6]
    @test mesh[:elemtocoord][:,:, 2] == [1 2 1 2; 5 5 6 6]
    @test mesh[:elemtocoord][:,:, 3] == [1 2 1 2; 6 6 7 7]
    @test mesh[:elemtocoord][:,:, 4] == [0 1 0 1; 6 6 7 7]
    @test mesh[:elemtocoord][:,:, 5] == [0 1 0 1; 7 7 8 8]
    @test mesh[:elemtocoord][:,:, 6] == [0 1 0 1; 8 8 9 9]
    @test mesh[:elemtocoord][:,:, 7] == [1 2 1 2; 8 8 9 9]
    @test mesh[:elemtocoord][:,:, 8] == [1 2 1 2; 7 7 8 8]
    @test mesh[:elemtocoord][:,:, 9] == [2 3 2 3; 7 7 8 8]
    @test mesh[:elemtocoord][:,:,10] == [2 3 2 3; 8 8 9 9]
    @test mesh[:elemtocoord][:,:,11] == [3 4 3 4; 8 8 9 9]
    @test mesh[:elemtocoord][:,:,12] == [3 4 3 4; 7 7 8 8]
    @test mesh[:elemtocoord][:,:,13] == [3 4 3 4; 6 6 7 7]
    @test mesh[:elemtocoord][:,:,14] == [2 3 2 3; 6 6 7 7]
    @test mesh[:elemtocoord][:,:,15] == [2 3 2 3; 5 5 6 6]
    @test mesh[:elemtocoord][:,:,16] == [3 4 3 4; 5 5 6 6]

    @test mesh[:elemtoelem] ==
      [1   1   4  4  5  6   6  5   8   7  10   9  14   3   2  15
       2  15  14  3  8  7  10  9  12  11  11  12  13  13  16  16
       6   7   2  1  4  5   8  3  14   9  12  13  16  15  10  11
       4   3   8  5  6  1   2  7  10  15  16  11  12   9  14  13]

    @test mesh[:elemtoface] ==
      [1  2  2  1  1  1  2  2  2  2  2  2  2  2  2  2
       1  1  1  1  1  1  1  1  1  1  2  2  2  1  1  2
       4  4  4  4  4  4  4  4  4  4  4  4  4  4  4  4
       3  3  3  3  3  3  3  3  3  3  3  3  3  3  3  3]

    @test mesh[:elemtoordr] == ones(Int, size(mesh[:elemtoordr]))

    @test mesh[:elems] == 1:nelem
    @test mesh[:realelems] == 1:nelem
    @test mesh[:ghostelems] == nelem.+(1:0)

    @test length(mesh[:sendelems]) == 0

    @test mesh[:nabrtorank] == Int[]
    @test mesh[:nabrtorecv] == UnitRange{Int}[]
    @test mesh[:nabrtosend] == UnitRange{Int}[]
  end
end

@testset "Mappings" begin
  let
    comm = MPI.COMM_SELF
    x = (0:4,)
    mesh = connectmesh(comm, partition(comm, brickmesh(x, (true,))...)[1:4]...)

    N = 3
    d = length(x)
    nelem = prod(length.(x).-1)
    nface = 2d
    Nfp = (N+1)^(d-1)

    vmapM, vmapP = mappings(N, mesh[:elemtoelem], mesh[:elemtoface],
                            mesh[:elemtoordr])

    @test vmapM == reshape([1,4,5,8,9,12,13,16], Nfp, nface, nelem)
    @test vmapP == reshape([16,5,4,9,8,13,12,1], Nfp, nface, nelem)
  end

  let
    comm = MPI.COMM_SELF
    x = (-1:1, 0:1)
    p = (false, true)
    mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

    N = 2
    d = length(x)
    nelem = prod(length.(x).-1)
    nface = 2d
    Nfp = (N+1)^(d-1)

    vmapM, vmapP = mappings(N, mesh[:elemtoelem], mesh[:elemtoface],
                            mesh[:elemtoordr])

    @test vmapM == reshape([ 1, 4, 7,  # f=1 e=1
                             3, 6, 9,  # f=2 e=1
                             1, 2, 3,  # f=3 e=1
                             7, 8, 9,  # f=4 e=1
                            10,13,16,  # f=1 e=2
                            12,15,18,  # f=2 e=2
                            10,11,12,  # f=3 e=2
                            16,17,18], # f=4 e=2
                           Nfp, nface, nelem)

    @test vmapP == reshape([ 1, 4, 7,  # f=1 e=1
                            10,13,16,  # f=1 e=2
                             7, 8, 9,  # f=4 e=1
                             1, 2, 3,  # f=3 e=1
                             3, 6, 9,  # f=2 e=1
                            12,15,18,  # f=2 e=2
                            16,17,18,  # f=4 e=2
                            10,11,12], # f=3 e=2
                           Nfp, nface, nelem)
  end

  let
    comm = MPI.COMM_SELF
    x = (0:1, 0:1, -1:1)
    p = (false, true, false)
    mesh = connectmesh(comm, partition(comm, brickmesh(x, p)...)[1:4]...)

    N = 2
    d = length(x)
    nelem = prod(length.(x).-1)
    nface = 2d
    Np = (N+1)^d
    Nfp = (N+1)^(d-1)

    vmapM, vmapP = mappings(N, mesh[:elemtoelem], mesh[:elemtoface],
                            mesh[:elemtoordr])

    fmask = [ 1  3  1  7 1 19
              4  6  2  8 2 20
              7  9  3  9 3 21
             10 12 10 16 4 22
             13 15 11 17 5 23
             16 18 12 18 6 24
             19 21 19 25 7 25
             22 24 20 26 8 26
             25 27 21 27 9 27]

    @test vmapM == reshape([fmask[:]; fmask[:].+Np],
                           Nfp, nface, nelem)

    @test vmapP == reshape([fmask[:,1];
                            fmask[:,2];
                            fmask[:,4];
                            fmask[:,3];
                            fmask[:,5];
                            fmask[:,5].+Np;
                            fmask[:,1].+Np;
                            fmask[:,2].+Np;
                            fmask[:,4].+Np;
                            fmask[:,3].+Np;
                            fmask[:,6];
                            fmask[:,6].+Np],
                           Nfp, nface, nelem)
  end
end

@testset "Get Partition" begin
  let
    Nelem = 150
    (so, ss, rs) = BrickMesh.getpartition(MPI.COMM_SELF, Nelem:-1:1)
    @test so == Nelem:-1:1
    @test ss == [1, Nelem+1]
    @test rs == [1, Nelem+1]
  end

  let
    Nelem = 111
    code = [ones(1,Nelem); collect(Nelem:-1:1)']
    (so, ss, rs) = BrickMesh.getpartition(MPI.COMM_SELF, Nelem:-1:1)
    @test so == Nelem:-1:1
    @test ss == [1, Nelem+1]
    @test rs == [1, Nelem+1]
  end
end

@testset "Partition" begin
  (etv, etc, etb, fc) = brickmesh((-1:2:1,-1:2:1,-2:1:2), (true,true,true))
  (netv, netc, netb, nfc) = partition(MPI.COMM_SELF, etv, etc, etb, fc)[1:4]
  @test etv == netv
  @test etc == netc
  @test etb == netb
  @test fc == nfc
end

@testset "Comm Mappings" begin
  let
    N = 1
    d = 2
    nface = 2d

    commelems = [1,2,5]
    commfaces = BitArray(undef, nface, length(commelems))
    commfaces .= false
    nabrtocomm = [1:2, 3:3]

    vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

    @test vmapC == Int[]
    @test nabrtovmapC == UnitRange{Int64}[1:0, 1:0]
  end

  let
    N = 1
    d = 2
    nface = 2d

    commelems = [1, 2, 5]
    commfaces = BitArray([false false false;
                          false true  false;
                          false true  false;
                          false false true ])
    nabrtocomm = [1:2, 3:3]

    vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

    @test vmapC == [5, 6, 8, 19, 20]
    @test nabrtovmapC == UnitRange{Int64}[1:3, 4:5]
  end

  let
    N = 2
    d = 2
    nface = 2d

    commelems = [2, 4, 5]
    commfaces = BitArray([true  true  false;
                          false false false;
                          false true  true;
                          false false true ])
    nabrtocomm = [1:1, 2:3]

    vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

    @test vmapC == [10, 13, 16, 28, 29, 30, 31, 34, 37, 38, 39, 43, 44, 45]
    @test nabrtovmapC == UnitRange{Int64}[1:3, 4:14]
  end

  let
    N = 2
    d = 3
    nface = 2d

    commelems = [3, 4, 7, 9]
    commfaces = BitArray([true  true true  false;
                          false true false false;
                          false true false false;
                          false true false false;
                          false true true  true ;
                          false true false true ])
    nabrtocomm = [1:1, 2:4]

    vmapC, nabrtovmapC = commmapping(N, commelems, commfaces, nabrtocomm)

    @test vmapC == [55, 58, 61, 64, 67, 70, 73, 76, 79, 82, 83, 84, 85, 86,
                    87, 88, 89, 90, 91, 92, 93, 94, 96, 97, 98, 99, 100, 101,
                    102, 103, 104, 105, 106, 107, 108, 163, 164, 165, 166,
                    167, 168, 169, 170, 171, 172, 175, 178, 181, 184, 187,
                    217, 218, 219, 220, 221, 222, 223, 224, 225, 235, 236,
                    237, 238, 239, 240, 241, 242, 243]
    @test nabrtovmapC == UnitRange{Int64}[1:9, 10:68]
  end

end

max_stars_repo_name: emassoud/CLIMA

max_issues_repo_name: emassoud/CLIMA

max_forks_repo_name: emassoud/CLIMA


text:
using statistical_rethinking
using Test

@testset "statistical_rethinking.jl" begin
    # Write your tests here.
end


text:
#=
Examples of parsing a JSON file. Requires the JSON package.
=#

import JSON

# Parse vehicles.json which contains an array of objects
v = JSON.parsefile("vehicles.json")
model = v[1]["model"]
year = v[1]["year"]
println("Vehicle 1 is $model $year")

# Parse fruits.json which contains several name:value pairs
f = JSON.parsefile("fruits.json")
nberries = f["berries"]
println("Number of berries is $nberries")

max_stars_repo_name: wigging/julia-computing

max_issues_repo_name: wigging/julia-computing

max_forks_repo_name: wigging/julia-computing


text:
@testset "utils" begin
    source = joinpath(pkgdir(GraphMLDatasets), "test", "data")
    zipfile = joinpath(source, "test.zip")
    tgzfile = joinpath(source, "test.tgz")
    target_file = joinpath(source, "random.file")

    GraphMLDatasets.unzip_zip(zipfile)
    @test isfile(target_file)
    rm(target_file, force=true)

    GraphMLDatasets.unzip_tgz(tgzfile, source)
    @test isfile(target_file)
    rm(target_file, force=true)
end


text:
import Base: show

"""
$(TYPEDEF)
NOTE: Not used yet.
A status message from the server indicating the progress of
a request, e.g. addVariable or solveSession.
"""
struct StatusMessage
    requestId::String
    action::String
    state::String
    timestamp::String
    client::Client

function show(io::IO, s::StatusMessage)
    print(io, """StatusMessage:
    client: $(self.client)
    requestId: $(s.requestId)
    action: $(self.action)
    state: $(self.state)
    timestamp: $(self.timestamp)
    """
    )
end


text:
# Unit tests for JUDI linear operators (without PDE solves)
# Philipp Witte (pwitte.slim@gmail.com)
# May 2018
#
# Mathias Louboutin, mlouboutin3@gatech.edu
# Updated July 2020

# Tests
function test_transpose(Op)
    @test isequal(size(Op), size(conj(Op)))
    @test isequal(reverse(size(Op)), size(transpose(Op)))
    @test isequal(reverse(size(Op)), size(adjoint(Op)))
    @test isequal(reverse(size(Op)), size(transpose(Op)))
end

sub_dim(v::Vector{<:Integer}, i::Integer) = v[i:i]
sub_dim(v::Vector{<:Integer}, i::AbstractRange) = v[i]
sub_dim(v::Integer, i) = v

check_dims(v::Integer, vo::Integer) = v == vo
check_dims(v::Vector{<:Integer}, vo::Vector{<:Integer}) = v == vo
check_dims(v::Integer, vo::Vector{<:Integer}) = (length(vo) == 1 && v == vo[1])
check_dims(v::Vector{<:Integer}, vo::Integer) = (length(v) == 1 && v[1] == vo)

function test_getindex(Op, nsrc)
    so = size(Op)

    Op_sub = Op[1]
    @test isequal(Op_sub.model, Op.model)
    # Check sizes. Same dimensions but subsampled source
    @test issetequal(keys(size(Op_sub)), keys(size(Op)))
    s = size(Op_sub)
    for i=1:2
        for (k, v) ∈ s[i]
            vo = k == :src ? 1 : so[i][k]
            @test check_dims(v, sub_dim(vo, 1))
        end
    end

    inds = 1:nsrc
    Op_sub = Op[inds]
    @test isequal(Op_sub.model, Op.model)
    @test isequal(keys(size(Op_sub)), keys(size(Op)))
    s = size(Op_sub)
    for i=1:2
        for (k, v) ∈ s[i]
            vo = k == :src ? nsrc : so[i][k]
            @test check_dims(v, sub_dim(vo, inds))
        end
    end
end

model = example_model()
########################################################## judiModeling ###############################################
@testset "judiModeling Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiModeling nsrc=$(nsrc)" begin
        F_forward = judiModeling(model; options=Options())
        F_adjoint = adjoint(F_forward)

        test_transpose(F_forward)
        @test isequal(size(F_forward)[1], time_space(model.n))
        @test isequal(size(F_forward)[2], time_space(model.n))
        @test isequal(size(F_adjoint)[1], time_space(model.n))
        @test isequal(size(F_adjoint)[2], time_space(model.n))
        @test issetequal(keys(size(F_forward)[1]), [:time, :x, :z])
        @test issetequal(keys(size(F_forward)[2]), [:time, :x, :z])
        # Time is uninitialized until first multiplication since there is no info
        # on propagation time
        @test issetequal(values(size(F_forward)[1]), [[0], model.n...])
        @test Int(size(F_forward)[1]) == 0
        # Update size for check
        nt = 10
        size(F_forward)[1][:time] = [nt for i=1:nsrc]
        @test Int(size(F_forward)[1]) == nsrc * nt * prod(model.n)

        @test adjoint(F_adjoint) == F_forward
        @test isequal(typeof(F_forward), judiModeling{Float32, :forward})
        @test isequal(typeof(F_adjoint), judiModeling{Float32, :adjoint})

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        if VERSION>v"1.2"
            a = randn(Float32, model.n...)
            F2 = F_forward(;m=a)
            @test isapprox(F2.model.m, a)
            F2 = F_forward(Model(model.n, model.d, model.o, a))
            @test isapprox(F2.model.m, a)
            @test F2.model.n == model.n
        end
    end
end

@testset "judiPointSourceModeling Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiPointSourceModeling nsrc=$(nsrc)" begin
        src_geometry = example_src_geometry(nsrc=nsrc)
        F_forward = judiModeling(model; options=Options()) * judiProjection(src_geometry)'
        F_adjoint = adjoint(F_forward)

        test_transpose(F_forward)
        @test isequal(size(F_forward)[1], time_space_src(nsrc, src_geometry.nt, model.n))
        @test isequal(size(F_forward)[2], rec_space(src_geometry))
        @test isequal(size(F_adjoint)[1], rec_space(src_geometry))
        @test isequal(size(F_adjoint)[2], time_space_src(nsrc, src_geometry.nt, model.n))
        @test issetequal(keys(size(F_forward)[1]), (:src, :time, :x, :z))
        @test issetequal(keys(size(F_forward)[2]), (:src, :time, :rec))
        # With the composition, everything should be initialized
        @test issetequal(values(size(F_forward)[1]), (nsrc, src_geometry.nt, model.n...))
        @test issetequal(values(size(F_forward)[2]), (nsrc, src_geometry.nt, src_geometry.nrec))
        @test Int(size(F_forward)[1]) == prod(model.n) * sum(src_geometry.nt)
        @test Int(size(F_forward)[2]) == sum(src_geometry.nrec .* src_geometry.nt)

        @test adjoint(F_adjoint) == F_forward
        @test isequal(typeof(F_forward), judiPointSourceModeling{Float32, :forward})
        @test isequal(typeof(F_adjoint), judiDataModeling{Float32, :adjoint})

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        if VERSION>v"1.2"
            a = randn(Float32, model.n...)
            F2 = F_forward(;m=a)
            @test isapprox(F2.model.m, a)
            F2 = F_forward(Model(model.n, model.d, model.o, a))
            @test isapprox(F2.model.m, a)
            @test F2.model.n == model.n
        end
    end
end

@testset "judiDataModeling Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiDataModeling nsrc=$(nsrc)" begin
        rec_geometry = example_rec_geometry(nsrc=nsrc)
        F_forward = judiProjection(rec_geometry)*judiModeling(model; options=Options())
        F_adjoint = adjoint(F_forward)

        test_transpose(F_forward)
        @test isequal(size(F_forward)[1], rec_space(rec_geometry))
        @test isequal(size(F_forward)[2], time_space_src(nsrc, rec_geometry.nt, model.n))
        @test isequal(size(F_adjoint)[1], time_space_src(nsrc, rec_geometry.nt, model.n))
        # With the composition, everything should be initialized
        @test issetequal(values(size(F_forward)[2]), (nsrc, rec_geometry.nt, model.n...))
        @test issetequal(values(size(F_forward)[1]), (nsrc, rec_geometry.nt, rec_geometry.nrec))
        @test Int(size(F_forward)[2]) == prod(model.n) * sum(rec_geometry.nt)
        @test Int(size(F_forward)[1]) == sum(rec_geometry.nrec .* rec_geometry.nt)

        @test adjoint(F_adjoint) == F_forward
        @test isequal(typeof(F_forward), judiDataModeling{Float32, :forward})
        @test isequal(typeof(F_adjoint), judiPointSourceModeling{Float32, :adjoint})

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        if VERSION>v"1.2"
            a = randn(Float32, model.n...)
            F2 = F_forward(;m=a)
            @test isapprox(F2.model.m, a)
            F2 = F_forward(Model(model.n, model.d, model.o, a))
            @test isapprox(F2.model.m, a)
            @test F2.model.n == model.n
        end
    end
end

@testset "judiDataSourceModeling Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiDataSourceModeling nsrc=$(nsrc)" begin
        rec_geometry = example_rec_geometry(nsrc=nsrc)
        src_geometry = example_src_geometry(nsrc=nsrc)
        F_forward = judiModeling(model, src_geometry, rec_geometry; options=Options())
        F_adjoint = adjoint(F_forward)

        test_transpose(F_forward)
        @test isequal(size(F_forward)[2], rec_space(src_geometry))
        @test isequal(size(F_forward)[1], rec_space(rec_geometry))
        @test isequal(size(F_adjoint)[2], rec_space(rec_geometry))
        @test isequal(size(F_adjoint)[1], rec_space(src_geometry))
        @test issetequal(keys(size(F_forward)[1]), (:src, :time, :rec))
        @test issetequal(keys(size(F_forward)[2]), (:src, :time, :rec))
        # With the composition, everything should be initialized
        @test issetequal(values(size(F_forward)[2]), (nsrc, src_geometry.nt, src_geometry.nrec))
        @test issetequal(values(size(F_forward)[1]), (nsrc, rec_geometry.nt, rec_geometry.nrec))
        @test Int(size(F_forward)[2]) == sum(src_geometry.nrec .* src_geometry.nt)
        @test Int(size(F_forward)[1]) == sum(rec_geometry.nrec .* rec_geometry.nt)

        @test adjoint(F_adjoint) == F_forward
        @test isequal(typeof(F_forward), judiDataSourceModeling{Float32, :forward})
        @test isequal(typeof(F_adjoint), judiDataSourceModeling{Float32, :adjoint})

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        # get index
        test_getindex(F_forward, nsrc)
        test_getindex(F_adjoint, nsrc)

        if VERSION>v"1.2"
            a = randn(Float32, model.n...)
            F2 = F_forward(;m=a)
            @test isapprox(F2.model.m, a)
            F2 = F_forward(Model(model.n, model.d, model.o, a))
            @test isapprox(F2.model.m, a)
            @test F2.model.n == model.n
        end
    end
end

######################################################## judiJacobian ##################################################
@testset "judiJacobian Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiJacobian nsrc=$(nsrc)" begin
        rec_geometry = example_rec_geometry(nsrc=nsrc)
        src_geometry = example_src_geometry(nsrc=nsrc)
        wavelet = randn(Float32, src_geometry.nt[1])
        PDE = judiModeling(model, src_geometry, rec_geometry; options=Options())
        q = judiVector(src_geometry, wavelet)
        J = judiJacobian(PDE, q)
        # Check sizes
        @test isequal(size(J)[1], rec_space(rec_geometry))
        @test isequal(size(J)[2], space(model.n))
        @test issetequal(keys(size(J)[2]), (:x, :z))
        @test Int(size(J)[2]) == prod(model.n)
        @test Int(size(J)[1]) == sum(rec_geometry.nrec .* rec_geometry.nt)

        @test isequal(typeof(J), judiJacobian{Float32, :born, typeof(PDE)})
        @test isequal(J.F.rInterpolation.geometry, rec_geometry)
        @test isequal(J.F.qInjection.geometry, src_geometry)
        @test isequal(J.q.geometry, src_geometry)
        @test all(isequal(J.q.data[i], wavelet) for i=1:nsrc)
        test_transpose(J)

        # get index
        J_sub = J[1]
        @test isequal(J_sub.model, J.model)
        @test isequal(J_sub.F, J.F[1])
        @test isequal(J_sub.q, q[1])

        inds = nsrc > 1 ? (1:nsrc) : 1
        J_sub = J[inds]
        @test isequal(J_sub.model, J.model)
        @test isequal(J_sub.F, J.F[inds])
        @test isequal(J_sub.q, q[inds])

        inds = nsrc > 1 ? (1:nsrc) : 1
        J_sub = J[inds]
        @test isequal(J_sub.model, J.model)
        @test isequal(size(J_sub), size(J))
    end
end

######################################################## judiJacobian ##################################################
@testset "judiJacobianExtended Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiJacobianExtended nsrc=$(nsrc)" begin
        rec_geometry = example_rec_geometry(nsrc=nsrc)
        wavelet = randn(Float32, rec_geometry.nt[1])
        # Setup operators
        Pr = judiProjection(rec_geometry)
        F = judiModeling(model)
        Pw = judiLRWF(rec_geometry.dt, wavelet)
        w = judiWeights(randn(Float32, model.n); nsrc=nsrc)
        J = judiJacobian(Pr*F*Pw', w)
        # Check sizes
        @test isequal(size(J)[1], rec_space(rec_geometry))
        @test isequal(size(J)[2], space(model.n))
        @test issetequal(keys(size(J)[2]), (:x, :z))
        @test Int(size(J)[2]) == prod(model.n)
        @test Int(size(J)[1]) == sum(rec_geometry.nrec .* rec_geometry.nt)

        rec_geometry = example_rec_geometry(nsrc=nsrc)
        wavelet = randn(Float32, rec_geometry.nt[1])
        # Setup operators
        Pr = judiProjection(rec_geometry)
        F = judiModeling(model)
        Pw = judiLRWF(nsrc, rec_geometry.dt[1], wavelet)
        w = judiWeights(randn(Float32, model.n); nsrc=nsrc)
        PDE = Pr*F*Pw'
        J = judiJacobian(PDE, w)

        @test isequal(typeof(J), judiJacobian{Float32, :born, typeof(PDE)})
        @test isequal(J.F.rInterpolation.geometry, rec_geometry)
        for i=1:nsrc
            @test isapprox(J.F.qInjection.wavelet[i], wavelet)
            @test isapprox(J.q.data[i], w.data[i])
        end
        @test isequal(size(J)[2], space(model.n))
        test_transpose(J)

        # get index
        J_sub = J[1]
        @test isequal(J_sub.model, J.model)
        @test isapprox(J_sub.q.weights[1], J.q.weights[1])

        inds = 1:nsrc
        J_sub = J[inds]
        @test isequal(J_sub.model, J.model)
        @test isapprox(J_sub.q.weights, J.q.weights[inds])

        # Test Pw alone
        P1 = subsample(Pw, 1)
        @test isapprox(P1.wavelet, Pw[1].wavelet)
        @test isapprox(conj(Pw).wavelet, Pw.wavelet)
        @test size(conj(Pw)) == size(Pw)
        @test isapprox(transpose(Pw).wavelet, Pw.wavelet)

    end
end

####################################################### judiProjection #################################################

@testset "judiProjection Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiProjection nsrc=$(nsrc)" begin
        rec_geometry = example_rec_geometry(nsrc=nsrc)
        P = judiProjection(rec_geometry)
        # Check sizes and that's it's unitialized since not combined with any propagator
        @test size(P)[2] == time_space_src(get_nsrc(rec_geometry), rec_geometry.nt, 3)
        @test size(P)[1] == rec_space(rec_geometry)
        # Size is zero since un-initialized
        @test Int(size(P)[2]) == 0
        @test Int(size(P)[1]) ==  sum(rec_geometry.nrec .* rec_geometry.nt)

        @test isequal(typeof(P), judiProjection{Float32})
        @test isequal(P.geometry, rec_geometry)

        Pr_sub = P[1]
        @test isequal(get_nsrc(Pr_sub.geometry), 1)
        @test get_nsrc(Pr_sub.geometry) == 1

        inds = nsrc > 1 ? (1:nsrc) : 1
        Pr_sub = P[inds]
        @test isequal(Pr_sub.geometry, rec_geometry[inds])
        @test get_nsrc(Pr_sub.geometry) == nsrc
    end
end

@testset "judiLRWF Unit Test with $(nsrc) sources" for nsrc=[1, 2]
    @timeit TIMEROUTPUT "judiLRWF nsrc=$(nsrc)" begin
        src_geometry = example_src_geometry(nsrc=nsrc)
        wavelet = randn(Float32, src_geometry.nt[1])

        P = judiLRWF(src_geometry.dt, wavelet)
        @test isequal(typeof(P), judiLRWF{Float32})
        # Check sizes and that's it's unitialized since not combined with any propagator
        @test size(P)[1] == space_src(nsrc)
        @test size(P)[2] == time_space_src(nsrc, src_geometry.nt)
        # Size is zero since un-initialized
        @test Int(size(P)[1]) == 0
        @test Int(size(P)[2]) == 0

        for i=1:nsrc
            @test isequal(P.wavelet[i], wavelet)
            @test isequal(P.dt[i], src_geometry.dt[i])
        end
        @test length(P.wavelet) == nsrc

        P_sub = P[1]
        @test isequal(P_sub.wavelet[1], wavelet)
        @test isequal(P_sub.dt[1], src_geometry.dt[1])
        @test length(P_sub.wavelet) == 1

        inds = nsrc > 1 ? (1:nsrc) : 1
        Pr_sub = P[inds]
        for i=1:nsrc
            @test isequal(Pr_sub.wavelet[i], wavelet)
            @test isequal(Pr_sub.dt[i], src_geometry.dt[i])
        end
        @test length(Pr_sub.wavelet) == nsrc
    end
end


text: 1 -2 -3

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
# scale factors are standard deviations
# square them for (co)variance
struct MixtureResults{P,T}
    probs::SVector{P,T}
    scale_factors::SVector{P,T}
end

struct MixtureWorkingData{P,T,S,L,O,TI,F}
    state::DifferentiableObjects.BFGSState{P,T,S,L}
    ls::DifferentiableObjects.BackTracking2{O,T,TI}
    obj::F
    initial_x::PaddedMatrices.MutableFixedSizeVector{P,T,S,S}
end

function MixtureWorkingData(mahal::MahalanobisDistances{T}, ::Val{P}, ::Val{BT} = Val(2)) where {T,P,BT}
    initial_x = PaddedMatrices.MutableFixedSizeVector{P,T}(undef)
    MixtureWorkingData(
        DifferentiableObjects.BFGSState(Val(P), T),
        DifferentiableObjects.BackTracking2{T}(Val(BT)),
        OnceDifferentiable(mahal, initial_x),
        initial_x
    )
end

struct Weight{P,T}
    ω::SVector{P,T}
end
mixture_weights(w::Weight) = exp.(w.ω)
function log_jac(w::Weight)
    sum(w.ω)
end

@inline Χ²₃_pdf(x) = sqrt(x) * exp(-x/2) #/ sqrt(2π)
@inline Χ₃_pdf(x) = (x² = x^2; x² * exp(-x²/2)) #/ sqrt(2π)

@generated function (d::MahalanobisDistances)(x::PaddedMatrices.AbstractMutableFixedSizeVector{TwoPp1,T}) where {T,TwoPp1}
    P, R = divrem(TwoPp1,2)
    R == 1 || throw("Only odd dimensions supported; probability is vector presented as 1 less.")
    quote
        @nexprs 4 i -> loop_i = zero(T)
        # $(Expr(:meta, :inline))
        π, target = log_jac_probs(Simplex{$(P+1),$T,$P}(SVector{$P,$T}( @ntuple $P p -> x[p]  )))
        untransformed_weights = Weight{$(P+1),$T}(SVector{$(P+1),$T}( @ntuple $(P+1) p -> x[$P+p] ) )
        target -= log_jac(untransformed_weights)
        w = @. exp( - untransformed_weights.ω)
        # @show π, w
        # loop = zero(T)
        # @show w, π
        N = length(d)
        N4, r = divrem(N, 4)
        # We accumulate 4 seperate sums to increase the accuracy of the summayion (pairwise summation algorithm)
        # If we naively did target +=, we occasionally get convergence failure due to numerical problems.
        # Accumulating the loop on a seperate loop variable seemed to solve the problem.
        # Extending that to pairwise summation, with some number (4, currently) seperate accumulators
        # should further increase accuracy.

        # If we had a random random order, I would accumulate all four in a single loop.
        # But we are likely to sort the Mahalanobis distances, because of use in the PercentileMatching algorthm
        # So, sorting gurantees that each of the loops should be a sum of relatively similarly valued "m"s.
        # This, again, ought to increase accuracy. (summing numbers of different magnitudes is less accurate)
        @nexprs 4 i -> @inbounds for n ∈ 1:N4
            loop_i += log(sum((@ntuple $(P+1) p -> Χ₃_pdf(d[n+(i-1)*N4]*w[p])*w[p]*π[p] )))
        end
        # Catch remainder of the loop
        @inbounds for n ∈ N-r+1:N
            loop_1 += log(sum((@ntuple $(P+1) p -> Χ₃_pdf(d[n]*w[p])*w[p]*π[p] )))
        end
        # we normalize by 100N, so that the magnitude scale with sample size
        # The reason for downscaling the magnitude at all is that we optimize this function
        # via the BFGS algorithm. This is a quasi-Newton algorithm.
        # These algorithms are iterative, repeatedly taking steps to find a local minimizer.
        # As a quasi-Newton method, it builds an approximation to the Hessian from changes
        # in the gradient between each of these iterations.
        # Hessians are used to adjust the gradient in determining the next step.
        # However, on the first iteration, the identity matrix is used to approximate the Hessian.
        # Thus, if the gradient is extreme because the magnitude of the function is extreme,
        # it is prone to taking wildly huge steps, landing in poorly conditioned regions of the
        # parameter space, and end up getting stuck and terminating before finding its way back out.
        # (Badly behaved [very non-quadratic] regions can also lead to building poor Hessian approximations,
        #  further slowing down any potential recovery.)
        # By downscaling, we temper the behavior, ensuring more conservative steps before a reasonable
        # Hessian approximation is computed.
        - (target + sum((@ntuple 4 loop))) / ( $(20TwoPp1) * N)
    end
end

@generated function initial_val(v::SVector{P}) where P
    :(vcat((@SVector zeros($(P-1))), v))
end
@generated function extract_values(v::Union{SVector{TwoPp1,T},PaddedMatrices.AbstractFixedSizeVector{TwoPp1,T}}) where {TwoPp1,T}
    P = TwoPp1 >> 1
    quote
        probs = probabilities(Simplex{$(P+1),$T,$P}(SVector{$P,$T}(@ntuple $P p -> v[p])))
        mixture_sf = mixture_weights(Weight{$(P+1),T}(SVector{$(P+1),$T}(@ntuple $(P+1) p -> v[$P+p])))
        MixtureResults(probs, mixture_sf)
    end
end

function mixture_fit(mwd::MixtureWorkingData{TwoPp1, T}, sat_number, prop_point) where {TwoPp1,T}
    for i ∈ 1:20
        opt, scale = DifferentiableObjects.optimize_scale!(mwd.state, mwd.obj, randn!(mwd.initial_x), mwd.ls, 10f0, 1f-5)
        # opt = soptimize(mahals, 0.1 * (@SVector randn(2P+1)))
        all(isfinite, mwd.state.x_old) || continue
        return extract_values(mwd.state.x_old)
    end
    @warn "Mixture failed 20 times. Mixture scale factors `1`. Propogation point: $prop_point, Satellite number: $sat_number."
    extract_values(@SVector zeros(T, TwoPp1))
end

function mixture_fit!(res::Vector{MixtureResults{P,T}}, mahals::Vector{<:MahalanobisDistances}) where {P,T}

    @inbounds for i ∈ eachindex(mahals)
        res[i] = mixture_fit(mahals[i], Val(P))
    end

end

function sample_distances!(rng::AbstractRNG, d::WeightedSamples, res::MixtureResults{P}) where P
    @fastmath for n ∈ 0:P:length(d)-P
        dist = sqrt(abs2(randn(rng)) + abs2(randn(rng)) + abs2(randn(rng)))
        @inbounds for p ∈ 1:P
            d[n + p] = (res.scale_factors[p]*dist, res.probs[p])
        end
    end
    d
end
sample_distances!(d::WeightedSamples, res::MixtureResults) = sample_distances!(Random.GLOBAL_RNG, d, res)

struct MixturePcs{P,T}
    probs::SVector{P,T}
    scale_factors::SVector{P,T}
end

function sample_Pc(rng::AbstractRNG, res::MixtureResults{P,T}, r1, v1, c1, r2, v2, c2, HBR::T2) where {P,T,T2}
    MixturePcs(
        T2.(res.probs),
        SVector{P,T2}(ntuple(p -> pc2dfoster_RIC(r1, v1, c1, r2, v2, c2 * abs2(res.scale_factors[p]), HBR), Val(P)))
    )
end
function calc_Pc(rng, res::AbstractMatrix{T1}, r1, v1, c1, r2, v2, c2, HBR::T2) where {T1,T2}
    T = promote_type(T1,T2)
    Pc = zero(T)
    # column 1 must be of probs
    # column 2 must be of scale factors
    for i ∈ 1:size(res,2)
        p = res[i,1]
        s = res[i,2]
        ### Perhaps do a safer check?
        ### Throw an error, or provide a warning?
        (isfinite(p) && isfinite(s)) || continue
        Pc = muladd(pc2dfoster_RIC(r1,v1,c1,r2,v2,c2 * s^2, HBR), p, Pc)
    end
    Pc
end

function Distributions.logpdf(res::MixtureResults{P,T}, x::SVector{3}) where {P,T}
    probs = res.probs
    scale_factors = res.scale_factors

    m = T(-0.5)*(x' * x)
    ms = zero(T)
    @inbounds @fastmath @simd for p ∈ 1:P
        s = scale_factors[p]
        ms += probs[p] * SLEEFwrap.exp( m / s^2 - T(3)*SLEEFwrap.log(s) )
    end
    # log(ms) - T(1.5log(2π))
    log(ms) - 1.5log(2π)
end


text:
# Save some bytes by using ints for switch variable
const UpdateObservable = "0"
const OnjsCallback = "1"
const EvalJavascript = "2"
const JavascriptError = "3"
const JavascriptWarning = "4"
const RegisterObservable = "5"
const JSDoneLoading = "8"
const FusedMessage = "9"

"""
    request_to_sessionid(request; throw = true)

Returns the session and browser id from request.
With throw = false, it can be used to check if a request
contains a valid session/browser id for a websocket connection.
Will return nothing if request is invalid!
"""
function request_to_sessionid(request; throw=true)
    if length(request.target) >= 42 # for /36session_id/4browser_id/
        session_browser = split(request.target, "/", keepempty=false)
        if length(session_browser) == 2
            sessionid, browserid = string.(session_browser)
            if length(sessionid) == 36 && length(browserid) == 4
                return sessionid, browserid
            end
        end
    end
    if throw
        error("Invalid sessionid: $(request.target)")
    else
        return nothing
    end
end

function html(body)
    return HTTP.Response(200, ["Content-Type" => "text/html", "charset" => "utf-8"], body=body)
end

function response_404(body="Not Found")
    return HTTP.Response(404, ["Content-Type" => "text/html", "charset" => "utf-8"], body=body)
end

function replace_url(match_str)
    key_regex = r"(/assetserver/[a-z0-9]+-.*?):([\d]+):[\d]+"
    m = match(key_regex, match_str)
    key = m[1]
    path = assetserver_to_localfile(string(key))
    return path * ":" * m[2]
end

const ASSET_URL_REGEX = r"http://.*/assetserver/([a-z0-9]+-.*?):([\d]+):[\d]+"

"""
    handle_ws_message(session::Session, message)

Handles the incoming websocket messages from the frontend
"""
function handle_ws_message(session::Session, message)
    isempty(message) && return
    data = MsgPack.unpack(message)
    typ = data["msg_type"]
    if typ == UpdateObservable
        registered, obs = session.observables[data["id"]]
        Base.invokelatest(update_nocycle!, obs, data["payload"])
    elseif typ == JavascriptError
        show(stderr, JSException(data))
    elseif typ == JavascriptWarning
        @warn "Error in Javascript: $(data["message"])\n)"
    elseif typ == JSDoneLoading
        if data["exception"] !== "null"
            exception = JSException(data)
            show(stderr, exception)
            session.init_error[] = exception
        end
        session.on_websocket_ready(session)
    else
        @error "Unrecognized message: $(typ) with type: $(typeof(type))"
    end
end

function handle_ws_error(e)
    if !(e isa WebSockets.WebSocketClosedError || e isa Base.IOError)
        @warn "error in websocket handler!"
        Base.showerror(stderr, e, Base.catch_backtrace())
    end
end

function handle_ws_connection(application::Server, session::Session, websocket::WebSocket)
    # We need two tries here
    # First, isopen may throw -.-...
    # Second, we need the finally to guruantee to delete + close the session
    # The inner try is of course to not break the loop on error
    try
        @debug("opening ws connection for session: $(session.id)")
        while isopen(websocket)
            try
                bytes = read(websocket)
                handle_ws_message(session, bytes)
            catch e
                handle_ws_error(e)
            end
        end
    catch e
        handle_ws_error(e)
    finally
        # This always needs to happen, which is why we need a try catch!
        @debug("Closing: $(session.id)")
        close(session)
        delete!(application.sessions, session.id)
    end
end

"""
    handles a new websocket connection to a session
"""
function websocket_handler(context, websocket::WebSocket)
    request = context.request; application = context.application
    sessionid_browserid = request_to_sessionid(request, throw=true)
    sessionid, browserid = sessionid_browserid
    # Look up the connection in our sessions
    if haskey(application.sessions, sessionid)
        session = application.sessions[sessionid]
        if isopen(session)
            # Would be nice to not error here - but I think this should never
            # Happen, and if it happens, we need to debug it!
            error("Session already has connection")
        end
        session.connection[] = websocket
        handle_ws_connection(application, session, websocket)
    else
        # This happens when an old session trys to reconnect to a new app
        # We somehow need to figure out better, how to recognize this
        @debug("Unregistered session id: $sessionid. Sessions: $(collect(keys(application.sessions)))")
    end
end