text:
##plottting
using Statistics, Scallops, RayTraceEllipsoids
using Makie, Colors, GeometryTypes
using AbstractPlotting: textslider
x2h(x, a, b, positive) = Scallops.signit(positive)*b/a*sqrt(a^2 - x^2)
x2y(x, a, b, c, positive) = c + x2h(x, a, b, positive)
function Δx2y(x, m1, m2)
    abs(x2y(x, m1.rxy, m1.rz, m1.cz, m1.positive) - x2y(x, m2.rxy, m2.rz, m2.cz, m2.positive))
end
m1ltm2(x, m1, m2) = x2y(x, m1.rxy, m1.rz, m1.cz, m1.positive) < x2y(x, m2.rxy, m2.rz, m2.cz, m2.positive)
function getintersection(m1, m2)
    r = min(m1.rxy, m2.rxy)
    xs = r/2:r
    i = findfirst(x -> m1ltm2(x, m1, m2), xs)
    xs[i - 1]
end
y2h(y, c) = y - c
y2x(y, a, b, c, positive) = Scallops.signit(positive)*a*sqrt(1 - ((y - c)/b)^2)
function populateθ0!(layers)
    lensl = layers[:lens]
    x = getintersection(lensl.dis_membrane, lensl.pro_membrane)
    m = lensl.dis_membrane
    y = x2y(x, m.rxy, m.rz, m.cz, m.positive)
    for l in layers
        m = l.dis_membrane
        m.θ0 = atan(y2h(y, m.cz), y2x(y, m.rxy, m.rz, m.cz, m.positive))
    end
    layers
end
function θ2coordinate(θ, a, b, c)
    sinθ, cosθ = sincos(θ)
    r = a*b/sqrt((b*cosθ)^2 + (a*sinθ)^2)
    r*cosθ, c + r*sinθ
end
function ellipse_coordinates!(points, θ0, a, b, c, n, keepcorners)
    if θ0 < zero(θ0)
        if keepcorners
            θs = range(π - θ0, stop=2π + θ0, length = n)
        else
            θs = range(π - θ0, stop=2π + θ0, length = n + 1)[2:end]
        end
    else
        if keepcorners
            θs = range(θ0, stop=π - θ0, length = n)
        else
            θs = range(θ0, stop=π - θ0, length = n + 1)[1:end-1]
        end
    end
    for i in 1:n
        points[i] = θ2coordinate(θs[i], a, b, c)
    end
end
function getpoly(layer, n)
    points = Array{Point2f0}(undef, 2n)
    m1 = layer.dis_membrane
    m2 = layer.pro_membrane
    keepcorners = sign(m1.θ0) == sign(m2.θ0)
    ellipse_coordinates!(view(points, 1:n), m1.θ0, m1.rxy, m1.rz, m1.cz, n, keepcorners)
    ellipse_coordinates!(view(points, 2n:-1:n+1), m2.θ0, m2.rxy, m2.rz, m2.cz, n, keepcorners)
    points
end
function morph2layers2(morphz)
    layers = morph2layers(morphz)
    populateθ0!(layers)
end
function getvertices(layers, k, n)
    l = layers[k]
    getpoly(l, n)
end
layers2color(layers, k) = Gray(1 - (layers[k].tissue.ri - 1.32)/0.18)
layers2color(layers, k, m, M) = Gray(1 - (layers[k].tissue.ri - m)/(M - m))
struct ColoredPoint
    coord::Point2f0
    color::Gray{Float64}
end
ColoredPoint(r::DeadRay) = ColoredPoint(Point2f0(NaN, NaN), 1.)
ColoredPoint(r::AbstractRay) = ColoredPoint(Point2f0(r.orig[1], r.orig[3]), 1 - r.int)
function fillpoints!(points, l, b, ous)
    i = 1
    r = l(EmptyRay, b)
    points[i] = ColoredPoint(r)
    for ou in ous
        r = trace!(ou, r)
        i += 1
        points[i] = ColoredPoint(r)
    end
end
function getrays(ous, l; nrays = 21, ϵ = 1e-4)
    npoints = length(ous) + 1
    rays = [[ColoredPoint((NaN, NaN), Gray(0)) for j in 1:npoints] for i in 1:nrays]
    for (ray, b) in zip(rays, range(ϵ, stop=1 - ϵ, length=nrays))
        fillpoints!(ray, l, b, ous)
    end
    rays
end
mm2μm(x) = x*1000
function getranges(min_aperture, max_aperture, aperture_n, min_morph, max_morph, morph_n, min_distance, max_distance, distance_n, min_θ, max_θ, θ_n)
    apertures = range(min_aperture, stop = max_aperture, length = aperture_n)
    morphzs = collect(range(min_morph, stop = max_morph, length = morph_n - 1))
    push!(morphzs, 1)
    sort!(morphzs)
    distances = exp10.(range(log10(min_distance), stop = log10(max_distance), length = distance_n)) #in mm
    θs = range(min_θ, stop = max_θ, length = θ_n)
    (aperture = apertures, morphz = morphzs, distance = distances, θ = θs)
end
function getobservables(apertures, morphzs, distances, θs)
    aperture_slider, aperture = textslider(apertures, "Aperture (μm)", start=mean(apertures))
    morphz_slider, morphz = textslider(morphzs, "Morph", start=1)
    distance_slider, distance = textslider(distances, "Distance (mm)")
    θ_slider, θ = textslider(θs, "Angle (°)", start=0)
    (aperture = aperture_slider, morphz = morphz_slider, distance = distance_slider, θ = θ_slider), (aperture = aperture, morphz = morphz, distance = distance, θ = θ)
end
function ray_trace_applet(; min_aperture = 200, max_aperture = 400, aperture_n = max_aperture - min_aperture + 1, min_morph = 0.8, max_morph = 1.48, morph_n = 100, min_distance = 0.5, max_distance = 1000, distance_n = 100, min_θ = -38, max_θ = 38, θ_n = max_θ - min_θ + 1, limits = FRect(-400,-400,800,800), nvertices = 11, photoreceptor_radius = 5, absorption_coefficient = (dis_ret = 0.00667, pro_ret = 0.00667))
    ranges = getranges(min_aperture, max_aperture, aperture_n, min_morph, max_morph, morph_n, min_distance, max_distance, distance_n, min_θ, max_θ, θ_n)
    sliders, observables = getobservables(ranges...)
    all = lift(createobj, observables.morphz, photoreceptor_radius, absorption_coefficient)
    sc = Scene(scale_plot=false, limits=limits)
    vertices = lift(all) do all
        populateθ0!(all.layers)
        todo = [k for k in keys(all.layers) if k ∉ (:water, :mirror)]
        map(todo) do k
            getvertices(all.layers, k, nvertices)
        end
    end
    colors = lift(all) do all
        todo = [k for k in keys(all.layers) if k ∉ (:water, :mirror)]
        m, M = extrema(all.layers[k].tissue.ri for k in todo)
        Δ = M - m
        m -= 0.1Δ
        M += 0.1Δ
        if Δ == 0
            m = copy(M)
            M += 1
        end
        map(todo) do k
            layers2color(all.layers, k, m, M)
        end
    end
    poly!(sc, vertices, color = colors)
    rays = lift(observables.distance, observables.aperture, all, observables.θ) do distance, aperture, all, θ
        l = Light(mm2μm(distance), aperture, all.rc, deg2rad(θ), all.nodalz)
        getrays(all.ous, l)
    end
    segments = lift(rays) do rays
        [ray[i].coord => ray[i+1].coord for ray in rays for i in 1:length(ray) - 1]
    end
    colors = lift(rays) do rays
        [mean((ray[i].color, ray[i+1].color)) for ray in rays for i in 1:length(ray) - 1]
    end
    linesegments!(sc, segments, color=colors)
    scatter!(sc, lift(a -> [Point2f0(0, a.nodalz)], all), markersize=10, color=:red)
    photoreceptors = lift(all) do all
        map([Scallops.RETINAS...]) do retina
            z1, z2 = Scallops.getzs(all.layers[retina], photoreceptor_radius)
            Rectangle{Float32}(-photoreceptor_radius, z1, 2photoreceptor_radius, z2 - z1)
        end
    end
    poly!(sc, photoreceptors, color=:green)
    axis = sc[Axis];
    axis[:names, :axisnames] = ("X (μm)", "Y (μm)");
    vbox(hbox(sliders...), sc)
end


text:
# To be used instead of _PM.variable_ne_branch_indicator() - supports deduplication of variables
function variable_ne_branch_indicator(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    first_n = first_id(pm, nw, :hour, :scenario)
    if nw == first_n
        if !relax
            z_branch_ne = _PM.var(pm, nw)[:branch_ne] = JuMP.@variable(pm.model,
                [l in _PM.ids(pm, nw, :ne_branch)], base_name="$(nw)_branch_ne",
                binary = true,
                start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_branch, l), "branch_tnep_start", 1.0)
            )
        else
            z_branch_ne = _PM.var(pm, nw)[:branch_ne] = JuMP.@variable(pm.model,
                [l in _PM.ids(pm, nw, :ne_branch)], base_name="$(nw)_branch_ne",
                lower_bound = 0.0,
                upper_bound = 1.0,
                start = _PM.comp_start_value(_PM.ref(pm, nw, :ne_branch, l), "branch_tnep_start", 1.0)
            )
        end
    else
        z_branch_ne = _PM.var(pm, nw)[:branch_ne] = _PM.var(pm, first_n)[:branch_ne]
    end
    report && _PM.sol_component_value(pm, nw, :ne_branch, :built, _PM.ids(pm, nw, :ne_branch), z_branch_ne)
end

function variable_ne_branch_investment(pm::_PM.AbstractPowerModel; nw::Int=_PM.nw_id_default, relax::Bool=false, report::Bool=true)
    first_n = first_id(pm, nw, :hour, :scenario)
    if nw == first_n
        if !relax
            investment = _PM.var(pm, nw)[:branch_ne_investment] = JuMP.@variable(pm.model,
                [l in _PM.ids(pm, nw, :ne_branch)], base_name="$(nw)_branch_ne_investment",
                binary = true,
                start = 0
            )
        else
            investment = _PM.var(pm, nw)[:branch_ne_investment] = JuMP.@variable(pm.model,
                [l in _PM.ids(pm, nw, :ne_branch)], base_name="$(nw)_branch_ne_investment",
                lower_bound = 0.0,
                upper_bound = 1.0,
                start = 0
            )
        end
    else
        investment = _PM.var(pm, nw)[:branch_ne_investment] = _PM.var(pm, first_n)[:branch_ne_investment]
    end
    report && _PM.sol_component_value(pm, nw, :ne_branch, :investment, _PM.ids(pm, nw, :ne_branch), investment)
end

ext: jl


text:
using Gen
import Random
using FunctionalCollections

#####################
# uncollapsed model #
#####################

struct Params
    prob_outlier::Float64
    inlier_std::Float64
    outlier_std::Float64
    slope::Float64
    intercept::Float64
end

@compiled @gen function datum(x::Float64, params::Params)
    is_outlier::Bool = @addr(bernoulli(params.prob_outlier), :z)
    std::Float64 = is_outlier ? params.inlier_std : params.outlier_std
    y::Float64 = @addr(normal(x * params.slope + params.intercept, std), :y)
    return y
end

data = plate(datum)

function compute_data_change(inlier_std_change, outlier_std_change, slope_change, intercept_change)
    if all([c !== nothing && (c == NoChange() || !c[1]) for c in [
            inlier_std_change, outlier_std_change, slope_change, intercept_change]])
        NoChange()
    else
        nothing
    end
end

@compiled @gen function model(xs::Vector{Float64})
    n::Int = length(xs)
    inlier_std::Float64 = @addr(gamma(1, 1), :inlier_std)
    outlier_std::Float64 = @addr(gamma(1, 1), :outlier_std)
    slope::Float64 = @addr(normal(0, 2), :slope)
    intercept::Float64 = @addr(normal(0, 2), :intercept)
    params::Params = Params(0.5, inlier_std, outlier_std, slope, intercept)
    inlier_std_change::Union{Tuple{Bool,Float64},Nothing} = @change(:inlier_std)
    outlier_std_change::Union{Tuple{Bool,Float64},Nothing} = @change(:outlier_std)
    slope_change::Union{Tuple{Bool,Float64},Nothing} = @change(:slope)
    intercept_change::Union{Tuple{Bool,Float64},Nothing} = @change(:intercept)
    change::Union{NoChange,Nothing} = compute_data_change(
        inlier_std_change, outlier_std_change, slope_change, intercept_change)
    ys::PersistentVector{Float64} = @addr(data(xs, fill(params, n)), :data, change)
    return ys
end

###################
# collapsed model #
###################

import Distributions
import Gen: random, logpdf, get_static_argument_types
struct TwoNormals <: Distribution{Float64} end
const two_normals = TwoNormals()
function logpdf(::TwoNormals, x, mu, sigma1, sigma2)
    if sigma1 < 0 || sigma2 < 0
        return -Inf
    end
    l1 = Distributions.logpdf(Distributions.Normal(mu, sigma1), x) + log(0.5)
    l2 = Distributions.logpdf(Distributions.Normal(mu, sigma2), x) + log(0.5)
    m = max(l1, l2)
    m + log(exp(l1 - m) + exp(l2 - m))
end
function random(::TwoNormals, mu, sigma1, sigma2)
    mu + (rand() < 0.5 ? sigma1 : sigma2) * randn()
end
get_static_argument_types(::TwoNormals) = [Float64, Float64, Float64]

data_collapsed = plate(two_normals)

@compiled @gen function model_collapsed(xs::Vector{Float64})
    n::Int = length(xs)
    inlier_std::Float64 = @addr(gamma(1, 1), :inlier_std)
    outlier_std::Float64 = @addr(gamma(1, 1), :outlier_std)
    slope::Float64 = @addr(normal(0, 2), :slope)
    intercept::Float64 = @addr(normal(0, 2), :intercept)
    means::Vector{Float64} = broadcast(+, slope * xs, intercept)
    ys::PersistentVector{Float64} = @addr(
        data_collapsed(means, fill(inlier_std, n), fill(outlier_std, n)), :data)
    return ys
end

#######################
# inference operators #
#######################

@compiled @gen function slope_proposal(prev)
    slope::Float64 = get_assignment(prev)[:slope]
    @addr(normal(slope, 0.5), :slope)
end

@compiled @gen function intercept_proposal(prev)
    intercept::Float64 = get_assignment(prev)[:intercept]
    @addr(normal(intercept, 0.5), :intercept)
end

@compiled @gen function inlier_std_proposal(prev)
    inlier_std::Float64 = get_assignment(prev)[:inlier_std]
    @addr(normal(inlier_std, 0.5), :inlier_std)
end

@compiled @gen function outlier_std_proposal(prev)
    outlier_std::Float64 = get_assignment(prev)[:outlier_std]
    @addr(normal(outlier_std, 0.5), :outlier_std)
end

@compiled @gen function flip_z(prob::Float64)
    @addr(bernoulli(prob), :z)
end

data_proposal = at_dynamic(flip_z, Int)

function logsumexp(arr)
    min_arr = maximum(arr)
    min_arr + log(sum(exp.(arr .- min_arr)))
end

@compiled @gen function is_outlier_proposal(prev, i::Int)
	prev_z::Bool = get_assignment(prev)[:data => i => :z]
    @addr(data_proposal(i, (prev_z ? 0.0 : 1.0,)), :data)
end

@compiled @gen function observe_datum(y::Float64)
    @addr(dirac(y), :y)
end

observe_data = plate(observe_datum)

@compiled @gen function observer(ys::Vector{Float64})
    @addr(observe_data(ys), :data)
end

@gen function observer_collapsed(ys::Vector{Float64})
    for (i, y) in enumerate(ys)
        @addr(dirac(y), :data => i)
    end
end

Gen.load_generated_functions()

#####################
# generate data set #
#####################

Random.seed!(1)

prob_outlier = 0.5
true_inlier_noise = 0.5
true_outlier_noise = 5.0
true_slope = -1
true_intercept = 2
xs = collect(range(-5, stop=5, length=1000))
ys = Float64[]
for (i, x) in enumerate(xs)
    if rand() < prob_outlier
        y = true_slope * x + true_intercept + randn() * true_inlier_noise
    else
        y = true_slope * x + true_intercept + randn() * true_outlier_noise
    end
    push!(ys, y)
end

######################
# inference programs #
######################

function do_inference_collapsed(n)
    observations = get_assignment(simulate(observer_collapsed, (ys,)))

    # initial trace
    (trace, weight) = generate(model_collapsed, (xs,), observations)

    for i=1:n

        # steps on the parameters
        for j=1:5
            trace = mh(model_collapsed, slope_proposal, (), trace)
            trace = mh(model_collapsed, intercept_proposal, (), trace)
            trace = mh(model_collapsed, inlier_std_proposal, (), trace)
            trace = mh(model_collapsed, outlier_std_proposal, (), trace)
        end
		assignment = get_assignment(trace)
		println((assignment[:inlier_std], assignment[:outlier_std], assignment[:slope], assignment[:intercept]))
    end

	assignment = get_assignment(trace)
    return (assignment[:inlier_std], assignment[:outlier_std], assignment[:slope], assignment[:intercept])
end

function do_inference(n)
    observations = get_assignment(simulate(observer, (ys,)))

    # initial trace
    (trace, weight) = generate(model, (xs,), observations)

    for i=1:n

        # steps on the parameters
        for j=1:5
            trace = mh(model, slope_proposal, (), trace)
            trace = mh(model, intercept_proposal, (), trace)
            trace = mh(model, inlier_std_proposal, (), trace)
            trace = mh(model, outlier_std_proposal, (), trace)
        end

        # step on the outliers
        for j=1:length(xs)
            trace = mh(model, is_outlier_proposal, (j,), trace)
        end
		assignment = get_assignment(trace)
		println((assignment[:inlier_std], assignment[:outlier_std], assignment[:slope], assignment[:intercept]))
    end

    assignment = get_assignment(trace)
    return (assignment[:inlier_std], assignment[:outlier_std], assignment[:slope], assignment[:intercept])
end

#################
# run inference #
#################

using Test

(inlier_std, outlier_std, slope, intercept) = do_inference(100)
max_std = max(inlier_std, outlier_std)
min_std = min(inlier_std, outlier_std)
@test isapprox(min_std, 0.5, atol=1e-1)
@test isapprox(max_std, 5.0, atol=1e-0)
@test isapprox(slope, -1, atol=1e-1)
@test isapprox(intercept, 2, atol=2e-1)

(inlier_std, outlier_std, slope, intercept) = do_inference_collapsed(100)
max_std = max(inlier_std, outlier_std)
min_std = min(inlier_std, outlier_std)
@test isapprox(min_std, 0.5, atol=1e-1)
@test isapprox(max_std, 5.0, atol=1e-0)
@test isapprox(slope, -1, atol=1e-1)
@test isapprox(intercept, 2, atol=2e-1)

max_stars_repo_name: BoltzmannBrain/Gen

max_issues_repo_name: BoltzmannBrain/Gen

max_forks_repo_name: BoltzmannBrain/Gen


text:
@testset "Simulation Sequence" begin
    feedforward_chronologies = Dict(
        ("UC" => "HAUC") => Synchronize(periods = 24),
        ("HAUC" => "ED") => RecedingHorizon(),
        ("ED" => "AGC") => RecedingHorizon(),
    )
    ini_cond_chronology = InterStageChronology()
    order = Dict(1 => "DAUC", 2 => "HAUC", 3 => "ED", 4 => "AGC")
    intervals = Dict(
        "DAUC" => (Hour(24), Consecutive()),
        "HAUC" => (Hour(1), RecedingHorizon()),
        "ED" => (Minute(5), RecedingHorizon()),
        "AGC" => (Minute(1), Consecutive()),
    )
    horizons = Dict("DAUC" => 48, "HAUC" => 24, "ED" => 12, "AGC" => 6)

    test_sequence = SimulationSequence(
        order = order,
        feedforward_chronologies = feedforward_chronologies,
        step_resolution = Hour(24),
        horizons = horizons,
        intervals = intervals,
        ini_cond_chronology = ini_cond_chronology,
    )

    @test length(findall(x -> x == 4, test_sequence.execution_order)) == 24 * 60
    @test length(findall(x -> x == 3, test_sequence.execution_order)) == 24 * 12
    @test length(findall(x -> x == 2, test_sequence.execution_order)) == 24
    @test length(findall(x -> x == 1, test_sequence.execution_order)) == 1

    bad_order = Dict(1 => "DAUC", 5 => "HAUC", 3 => "ED", 4 => "AGC")
    @test_throws IS.InvalidValue SimulationSequence(
        order = bad_order,
        feedforward_chronologies = feedforward_chronologies,
        step_resolution = Hour(24),
        horizons = horizons,
        intervals = intervals,
        ini_cond_chronology = ini_cond_chronology,
    )

    @test_throws ArgumentError SimulationSequence(
        step_resolution = Hour(24),
        order = Dict(1 => "UC", 2 => "ED"),
        horizons = Dict("UC" => 24, "ED" => 12),
        intervals = Dict(
            "UC" => (Hour(24), Consecutive()),
            "ED" => (Hour(1), Consecutive()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_from_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        cache = Dict(("ED",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
        ini_cond_chronology = InterStageChronology(),
    )

    test_sequence = SimulationSequence(
        order = Dict(1 => "DAUC"),
        step_resolution = Hour(24),
        horizons = Dict("DAUC" => 24),
        intervals = Dict("DAUC" => (Hour(24), Consecutive())),
        ini_cond_chronology = InterStageChronology(),
    )

    @test isa(test_sequence.ini_cond_chronology, IntraStageChronology)
    @test test_sequence.execution_order == [1]

end

@testset "Test if Horizon and interval result in a discontinuous simulation" begin
    @test_throws IS.ConflictingInputsError SimulationSequence(
        step_resolution = Hour(24),
        order = Dict(1 => "UC", 2 => "ED"),
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        horizons = Dict("UC" => 24, "ED" => 12),
        intervals = Dict(
            "UC" => (Hour(2), Consecutive()),
            "ED" => (Hour(3), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_from_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        cache = Dict(("ED",) => TimeStatusChange(PSY.ThermalStandard, PSI.ON)),
        ini_cond_chronology = InterStageChronology(),
    )
end

@testset "Test if interval is shorter than resolution" begin
    @test_throws IS.ConflictingInputsError sequence = SimulationSequence(
        step_resolution = Hour(24),
        order = Dict(1 => "UC", 2 => "ED"),
        feedforward_chronologies = Dict(("UC" => "ED") => Synchronize(periods = 24)),
        horizons = Dict("UC" => 24, "ED" => 12),
        intervals = Dict(
            "UC" => (Minute(5), RecedingHorizon()),
            "ED" => (Minute(1), RecedingHorizon()),
        ),
        feedforward = Dict(
            ("ED", :devices, :Generators) => SemiContinuousFF(
                binary_from_stage = PSI.ON,
                affected_variables = [PSI.ACTIVE_POWER],
            ),
        ),
        ini_cond_chronology = InterStageChronology(),
    )
end


text:
module TestT3FFExplicit
end # module


text:
using Core: SSAValue
using Cassette

function isinteresting end
function alternative end

# Must return 8 exprs
function rewrite_call(ctx, stmt, extraslot, i)
    exprs = Any[]
    cond = stmt.args[1]        # already an SSAValue

    # order of expressions we pre-assign
    # SSA for each expression here
    #
    # input:
    #    f(x)
    #
    # After this pass:
    #    isinteresting(ctx, f, x) ? alternative(f, x) : f(x)
    #
    # After overdub_pass!:
    #    isinteresting(ctx, f, x) ? alternative(f, x) : overdub(ctx, f, x)

    ISINTERESTING,
    GOTOIFNOTINTERESTING,
    ALTERNATIVE,
    ALTERNATIVE_IN_SLOT,
    GOTO_LAST,
    ORIGINAL,
    ORIGINAL_IN_SLOT,
    LAST = i:i+20

    # ISINTERESTING
    push!(exprs, Expr(:call,
                      Expr(:nooverdub, isinteresting),
                      Expr(:contextslot),
                      stmt.args...))

    # GOTOIFNOTINTERESTING
    push!(exprs, Expr(:gotoifnot, SSAValue(ISINTERESTING), ORIGINAL))

    # ALTERNATIVE
    push!(exprs, Expr(:call,
                      Expr(:nooverdub, alternative),
                      Expr(:contextslot),
                      stmt.args...))

    # ALTERNATIVE_IN_SLOT
    push!(exprs, Expr(:(=), extraslot, SSAValue(ALTERNATIVE)))

    # GOTO_LAST
    push!(exprs, Core.GotoNode(LAST))

    # ORIGINAL
    push!(exprs, stmt)

    # ORIGINAL_IN_SLOT
    push!(exprs, Expr(:(=), extraslot, SSAValue(ORIGINAL)))

    # LAST
    push!(exprs, extraslot)

    exprs
end

function newslot!(ir)
    extraslot = gensym("tmp")
    push!(ir.slotnames, extraslot)
    push!(ir.slotflags, 0x00)
    Core.SlotNumber(length(ir.slotnames))
end

function rewrite_ir(ctx, ref)
    # turn
    #   f(x...)
    # into
    #   %i = interesting(ctx, f, x...)
    #   gotoifnot %i #g
    #   myslot = alternative(ctx, f, x...)
    #   goto #g+1
    #   #g: myslot = f(x...)

    ir = ref.code_info
    ir = copy(ir)

    Cassette.insert_statements!(ir.code, ir.codelocs,
        (stmt, i) -> Base.Meta.isexpr(stmt, :call) ? 8 : nothing,
        (stmt, i) -> (s = newslot!(ir); rewrite_call(ctx, stmt, s, i)))

    Cassette.insert_statements!(ir.code, ir.codelocs,
                                (stmt, i) -> Base.Meta.isexpr(stmt, :(=)) && stmt.args[1] isa Core.SlotNumber && stmt.args[2] isa Expr && stmt.args[2].head == :call ? 8 : nothing,
                                (stmt, i) -> rewrite_call(ctx, stmt.args[2], stmt.args[1], i))

    ir.ssavaluetypes = length(ir.code)

    # Core.Compiler.validate_code(ir)
    return ir
end

const CustomDispatchPass = Cassette.@pass rewrite_ir

## Test
#=

using Cassette

Cassette.@context AnyRational

c = AnyRational(pass=CPASS.CustomDispatchPass)

@inline function CPASS.isinteresting(ctx::AnyRational, f, args...)
    return any(x->x isa Rational, args)
end

@inline function CPASS.alternative(ctx::AnyRational, f, args...)
    Cassette.overdub(ctx, f, args...)
end

=#


text:
#= Instructions:

- Pkg.add("PkgBenchmark.jl")
- using PkgBenchmark

- results = benchmarkpkg("IntervalArithmetic")
- showall(results)

- results = judge("IntervalArithmetic", "v0.9.1")  # compare current version to that tag
- showall(results)
=#

using IntervalArithmetic

@benchgroup "Constructors" begin
    @bench "Interval" Interval(1, 2)
end

@benchgroup "@interval" begin

    @bench "@interval" @interval(0.1)
    @bench "pi" @interval(pi)
    @bench "expression" @interval(sin(0.1) + cos(0.2))
end

@benchgroup "Arithmetic" begin

    a = Interval(1, 2)
    b = Interval(3, 4)

    for op in (+, -, *, /)
        @bench string(op) $(op)($a, $b)
    end
end

@benchgroup "Elementary functions" begin
    for op in (exp, log, sin, tan)
        @bench string(op) $(op)($a)
    end
end

@benchgroup "Sum" begin

    sum1(N) = sum(Interval(i, i+1) for i in 1:N)
    sum2(N) = (one = Interval(1.0); sum(one / (i^2) for i in 1:N) )

    @bench "Sum1" sum1(1000)
    @bench "Sum2" sum2(1000)
end


text:
function dare(A,B,Q,R,S=0)
	#
	# solves discrete-time algebraic Riccati equation
	# X = A'XA - (A'XB+S)(R+B'XB)^{-1}(B'XA+S') + Q
	#
	# R^{-1} is not computed
	#
	# javier.cara@upm.es, 2016-02
	#

	# Dimensions
	n,m = size(B)

	if S== 0
		S = zeros(n,m)
	end

	L = [A zeros(n,n) B;Q -eye(n) S;S' zeros(m,n) R]
	M = [eye(n) zeros(n,n) zeros(n,m); zeros(n,n) -A' zeros(n,m);zeros(m,n) -B' zeros(m,m)]

	# Compute the eigenvalue decomposition
	(d,v) = eig(L,M)

	# Sort the eigenvalues
	ew = abs(d)
	pos = sortperm(ew)

	# Compute X
	v1 = v[1:n,pos[1:n]]
	v2 = v[n+1:2*n,pos[1:n]]
	X = real(v2/v1)

	return X

end

######
function dare_test()
    #=
    BibText
    @TECHREPORT{Benner95acollection,
    author = {Peter Benner and Alan J. Laub and Volker Mehrmann},
    title = {A Collection of Benchmark Examples for the Numerical Solution of
    Algebraic Riccati Equations II: Discrete-Time Case},
    institution = {FAK. F. MATHEMATIK, TU CHEMNITZ--ZWICKAU},
    year = {1995}
    }

    nov-2016
    =#

	# ------------------------------------------------------------------
	# ejm1 collection of benchmark examples
	A1 = [4 3;-4.5 -3.5]
	B1 = [1 -1]'
	R1 = 1
	Q1 = [9 6;6 4]
	X1 = (1+sqrt(5))/2*[9 6;6 4]

	X11 = dare(A1,B1,Q1,R1)

	# --------------------------------------------------------------------
	# ejm computed with mathematica
	Am = [1 -1 1;0 1 1;0 0 1]
	Bm = [1 0;1 0;0 1]
	Rm = [10 0;0 0.1]
	Qm = [10 0 0;0 1 0;0 0 0.1]
	Xm = [42.2835 -68.5247 -3.94783;-68.5247 154.043 16.0017;-3.94783 16.0017 8.33197]

	Xmm = dare(Am,Bm,Qm,Rm)

	#####

	return X1,X11,Xm,Xmm

end


text:
append!(empty!(LOAD_PATH), Base.DEFAULT_LOAD_PATH)
using Pkg

################################################################################
# Generate notebooks
# ################################################################################
exampledir = joinpath(@__DIR__, "..", "examples")
Pkg.activate(exampledir)
Pkg.instantiate()
include(joinpath(exampledir, "generate_notebooks.jl"))

################################################################################
# Build simulation environments
################################################################################
pkgdir = joinpath(@__DIR__, "..")
Pkg.activate(pkgdir)

using JLD2
using Symbolics
using LinearAlgebra
using Scratch

include(joinpath("../src/utils.jl"))
include(joinpath("../src/simulator/model.jl"))
include(joinpath("../src/simulator/dimensions.jl"))
include(joinpath("../src/simulator/residual.jl"))
include(joinpath("../src/simulator/codegen.jl"))

include(joinpath("../src/robots/integrator.jl"))
include(joinpath("../src/robots/codegen.jl"))

include(joinpath("../src/solver/cones.jl"))

include(joinpath(@__DIR__, "..", "src/robots/hopper/model.jl"))  # hopper special case
include(joinpath(@__DIR__, "..", "src/robots/hopper/simulator.jl"))

path_robots = @get_scratch!("robots")

robots = [:hopper, :biped, :quadruped]

for robot in robots
    # robot model
    include(joinpath(@__DIR__, "..", "src/robots", String(robot), "model.jl"))

    # expr path
    path_expr = joinpath(path_robots, String(robot) * ".jld2")

    # kinematics
    contact_kinematics = eval(Symbol(String(robot) * "_contact_kinematics"))
    contact_kinematics_jacobians = eval(Symbol(String(robot) * "_contact_kinematics_jacobians"))

    # codegen
    mass_matrix, dynamics_bias = codegen_dynamics(eval(robot))
    r_model, rz_model, rθ_model = codegen_residual(eval(robot), mass_matrix, dynamics_bias, contact_kinematics, contact_kinematics_jacobians,
        save=true, path=path_expr)
end


text:
@testset "TopologicalRelation" begin
  @testset "FullTopology" begin
    elems = connect.([(1,2,3),(4,3,2)])
    t = FullTopology(elems)
    ∂ = Boundary{2,0}(t)
    @test ∂(1) == [1,2,3]
    @test ∂(2) == [4,3,2]
  end

  @testset "GridTopology" begin
    # 3 segments
    t = GridTopology(3)
    ∂ = Boundary{1,0}(t)
    @test ∂(1) == [1,2]
    @test ∂(2) == [2,3]
    @test ∂(3) == [3,4]

    # 2x3 quadrangles
    t = GridTopology(2, 3)
    ∂ = Boundary{2,0}(t)
    @test ∂(1) == [1,2,5,4]
    @test ∂(2) == [2,3,6,5]
    @test ∂(3) == [4,5,8,7]
    @test ∂(4) == [5,6,9,8]
    @test ∂(5) == [7,8,11,10]
    @test ∂(6) == [8,9,12,11]
    ∂ = Boundary{1,0}(t)
    @test ∂(1) == [1,2]
    @test ∂(2) == [2,3]
    @test ∂(3) == [4,5]
    @test ∂(4) == [5,6]
    @test ∂(5) == [7,8]
    @test ∂(6) == [8,9]
    @test ∂(7) == [10,11]
    @test ∂(8) == [11,12]
    @test ∂(9) == [1,4]
    @test ∂(10) == [2,5]
    @test ∂(11) == [3,6]
    @test ∂(12) == [4,7]
    @test ∂(13) == [5,8]
    @test ∂(14) == [6,9]
    @test ∂(15) == [7,10]
    @test ∂(16) == [8,11]
    @test ∂(17) == [9,12]

    # 2x3x2 hexahedrons
    t = GridTopology(2, 3, 2)
    ∂ = Boundary{3,0}(t)
    @test ∂(1) == [1,2,5,4,13,14,17,16]
    @test ∂(2) == [2,3,6,5,14,15,18,17]
    @test ∂(3) == [4,5,8,7,16,17,20,19]
    @test ∂(12) == [20,21,24,23,32,33,36,35]

    # invalid relations
    t = GridTopology(2, 3)
    @test_throws AssertionError Boundary{3,0}(t)
    @test_throws AssertionError Coboundary{0,3}(t)
    @test_throws AssertionError Adjacency{3}(t)
    @test_throws AssertionError Boundary{0,2}(t)
    @test_throws AssertionError Coboundary{2,0}(t)
  end

  @testset "HalfEdgeTopology" begin
    # 2 triangles
    elems = connect.([(1,2,3),(4,3,2)])
    t = HalfEdgeTopology(elems)
    ∂ = Boundary{2,0}(t)
    @test ∂(1) == [2,3,1]
    @test ∂(2) == [3,2,4]
    ∂ = Boundary{2,1}(t)
    @test ∂(1) == [1,3,2]
    @test ∂(2) == [1,4,5]
    ∂ = Boundary{1,0}(t)
    @test ∂(1) == [3,2]
    @test ∂(2) == [1,2]
    @test ∂(3) == [3,1]
    @test ∂(4) == [2,4]
    @test ∂(5) == [4,3]
    𝒞 = Coboundary{0,1}(t)
    @test 𝒞(1) == [2,3]
    @test 𝒞(2) == [4,1,2]
    @test 𝒞(3) == [3,1,5]
    @test 𝒞(4) == [5,4]
    𝒞 = Coboundary{0,2}(t)
    @test 𝒞(1) == [1]
    @test 𝒞(2) == [2,1]
    @test 𝒞(3) == [1,2]
    @test 𝒞(4) == [2]
    𝒞 = Coboundary{1,2}(t)
    @test 𝒞(1) == [2,1]
    @test 𝒞(2) == [1]
    @test 𝒞(3) == [1]
    @test 𝒞(4) == [2]
    @test 𝒞(5) == [2]
    𝒜 = Adjacency{0}(t)
    @test 𝒜(1) == [2,3]
    @test 𝒜(2) == [4,3,1]
    @test 𝒜(3) == [1,2,4]
    @test 𝒜(4) == [3,2]

    # 2 triangles + 2 quadrangles
    elems = connect.([(1,2,6,5),(2,4,6),(4,3,5,6),(1,5,3)])
    t = HalfEdgeTopology(elems)
    ∂ = Boundary{2,0}(t)
    @test ∂(1) == [1,2,6,5]
    @test ∂(2) == [6,2,4]
    @test ∂(3) == [6,4,3,5]
    @test ∂(4) == [3,1,5]
    ∂ = Boundary{2,1}(t)
    @test ∂(1) == [1,3,5,6]
    @test ∂(2) == [3,9,4]
    @test ∂(3) == [4,7,8,5]
    @test ∂(4) == [2,6,8]
    ∂ = Boundary{1,0}(t)
    @test ∂(1) == [1,2]
    @test ∂(2) == [3,1]
    @test ∂(3) == [6,2]
    @test ∂(4) == [4,6]
    @test ∂(5) == [5,6]
    @test ∂(6) == [1,5]
    @test ∂(7) == [4,3]
    @test ∂(8) == [3,5]
    @test ∂(9) == [2,4]
    𝒞 = Coboundary{0,1}(t)
    @test 𝒞(1) == [1,6,2]
    @test 𝒞(2) == [9,3,1]
    @test 𝒞(3) == [2,8,7]
    @test 𝒞(4) == [7,4,9]
    @test 𝒞(5) == [5,8,6]
    @test 𝒞(6) == [3,4,5]
    𝒞 = Coboundary{0,2}(t)
    @test 𝒞(1) == [1,4]
    @test 𝒞(2) == [2,1]
    @test 𝒞(3) == [4,3]
    @test 𝒞(4) == [3,2]
    @test 𝒞(5) == [3,4,1]
    @test 𝒞(6) == [2,3,1]
    𝒞 = Coboundary{1,2}(t)
    @test 𝒞(1) == [1]
    @test 𝒞(2) == [4]
    @test 𝒞(3) == [2,1]
    @test 𝒞(4) == [2,3]
    @test 𝒞(5) == [3,1]
    @test 𝒞(6) == [4,1]
    @test 𝒞(7) == [3]
    @test 𝒞(8) == [3,4]
    @test 𝒞(9) == [2]
    𝒜 = Adjacency{0}(t)
    @test 𝒜(1) == [2,5,3]
    @test 𝒜(2) == [4,6,1]
    @test 𝒜(3) == [1,5,4]
    @test 𝒜(4) == [3,6,2]
    @test 𝒜(5) == [6,3,1]
    @test 𝒜(6) == [2,4,5]

    # invalid relations
    elems = connect.([(1,2,3),(4,3,2)])
    t = HalfEdgeTopology(elems)
    @test_throws AssertionError Boundary{3,0}(t)
    @test_throws AssertionError Coboundary{0,3}(t)
    @test_throws AssertionError Adjacency{3}(t)
    @test_throws AssertionError Boundary{0,2}(t)
    @test_throws AssertionError Coboundary{2,0}(t)
  end
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ b36fb690-3447-11eb-3629-edbddae4ebd1
input = readlines("./input.txt");

# ╔═╡ 3d1fca00-3449-11eb-3e0e-5784380ebdbc
struct PasswordAndPolicy
	low :: Int
	high :: Int
	policy_letter :: Char
	password :: String
end

# ╔═╡ 4ab2f908-3449-11eb-299d-53d5a0a3044e
function parseInput(line)
	matches = match(r"(\d+)-(\d+)\s+(\w):\s+(.+)",line)
	matches == nothing && return
	PasswordAndPolicy(
		parse(Int,matches.captures[1]),
		parse(Int,matches.captures[2]),
		matches.captures[3][1], # capture only first char
		matches.captures[4]
	)
end

# ╔═╡ 473a9452-3449-11eb-1dd1-25d770220737
passwordsAndPolicies = parseInput.(input)

# ╔═╡ 6e67d06a-344b-11eb-15c6-ed3cfcfd8a83
"""
	Checks if password is valid based on Part 1 rules
"""
function Part1Valid(p)
	count(==(p.policy_letter),p.password) ∈ UnitRange(p.low,p.high)
end

# ╔═╡ 959afa6a-3827-11eb-1fc9-5d4a250374cc
"""
	Evaluates Part 1
"""
function Part1(passwordsAndPolicies)
	count(
			p -> Part1Valid(p),
			passwordsAndPolicies
		)
end

# ╔═╡ b4d6dd9a-3827-11eb-18af-753996dfdc4d
Part1(passwordsAndPolicies)

# ╔═╡ 30635df6-3828-11eb-3ea0-ad5512977175
"""
	Checks if password is valid based on Part 2 rules
"""
function Part2Valid(p)
	(p.password[p.low] == p.policy_letter) ⊻ (p.password[p.high] == p.policy_letter)
end

# ╔═╡ 604e4e30-3828-11eb-3660-0196265f0647
"""
	Evaluates Part 2
"""
function Part2(passwordsAndPolicies)
	count(
			p -> Part2Valid(p),
			passwordsAndPolicies
		)
end

# ╔═╡ 6a0e2614-3828-11eb-3445-0bfb392e9a26
Part2(passwordsAndPolicies)

# ╔═╡ Cell order:
# ╠═b36fb690-3447-11eb-3629-edbddae4ebd1
# ╠═3d1fca00-3449-11eb-3e0e-5784380ebdbc
# ╠═4ab2f908-3449-11eb-299d-53d5a0a3044e
# ╠═473a9452-3449-11eb-1dd1-25d770220737
# ╠═6e67d06a-344b-11eb-15c6-ed3cfcfd8a83
# ╠═959afa6a-3827-11eb-1fc9-5d4a250374cc
# ╠═b4d6dd9a-3827-11eb-18af-753996dfdc4d
# ╠═30635df6-3828-11eb-3ea0-ad5512977175
# ╠═604e4e30-3828-11eb-3660-0196265f0647
# ╠═6a0e2614-3828-11eb-3445-0bfb392e9a26


text:
using ParamLevelSet
using jInv.Mesh
using SparseArrays
using LinearAlgebra

function heavi(x, type="global", eps=0.1, thr=0.0) # [h,d] =

    x .= x - thr;

    if type == "global"
        h = 0.5 * (1 + 2/pi*atan(pi* x /eps))
        d = 1.0 / (eps*((x .^ 2 * pi^2) / eps^2 + 1))

    elseif type == "compact"
        # compact
        # h = 0*x;
        # d = 0*x;
        # id = find((x < epsi) & (x > -epsi));
        # h(id) = 0.5*(1 + x(id)/epsi + 1/pi*sin(pi*x(id)/epsi));
        # h(x >= epsi) = 1;
        # h(x <=-epsi) = 0;
        # d(id) = 0.5*(1/epsi)*(1 + cos(pi*x(id)/epsi))
    end

    return h, d
end

function recon2d_slices_pals!(u::Array{T, 3}, A, p, niter::Int=50, nbasis::Int=4) where {T <: AbstractFloat}

    global Jt, J0, mesh, m, Hessian, grad, r

    mesh, m = init_pals(size(u), nbasis)
    sigmaH = getDefaultHeaviside();
    phat = similar(p)
    r = similar(p)
    At = A'
    stepsize = 0.01
    nslice = size(p,2)

    HW = size(u,1)*size(u,2)
    szpslice = size(p,1)*size(p,3)
    J = spzeros(prod(size(p)), nbasis*5) # A * J0

    for i=1:niter
        # u0, JBuilder = eval_f_grad(mesh, m; computeJacobian=1, sigma=sigmaH)
        u0, JBuilder = MeshFreeParamLevelSetModelFunc(mesh, m; computeJacobian = 1)

        J0 = getSparseMatrix(JBuilder)

        # construct J
        Threads.@threads for ss=1:nslice
            J[szpslice*(ss-1)+1:szpslice*ss, :] .= A * J0[HW*(ss-1)+1:HW*ss, :]

            u_slice = view(u, :, :, ss)
            Av = A*vec(u_slice)

            # phat_slice = reshape(Av, size(p,1), size(p,3))
            p_vec_ = reshape(Av, size(p,1), size(p,3))
            @views r[:,ss,:] = p_vec_ - p[:,ss,:]
        end

        Jt = J'
        Hessian = Jt*J
        Hinv = pinv(Matrix(Hessian))
        grad = Jt * vec(r)
        m .-= stepsize * Hinv * grad

        fres = sum(r .^ 2) / length(r)
        @show fres, m
        # println("residual: $fres, $m")
    end

    # u0, JBuilder = eval_f_grad(mesh,m;computeJacobian=1, sigma=sigmaH)
    u0, JBuilder = MeshFreeParamLevelSetModelFunc(mesh, m; computeJacobian = 1)
    u .= reshape(u0, size(u))
    return u, m
end

max_stars_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_issues_repo_name: JuliaTomo/DynamicTomographyWithFlow

max_forks_repo_name: JuliaTomo/DynamicTomographyWithFlow


text:
## This script makes the example graphs we'll use in subsequent studies

##
using Plots
using GraphRecipes
using NearestNeighbors
using Random
using SparseArrays
using LinearAlgebra
##
""" This function is weird because it only gets half of the knn, but it makes
a useful picture, so I'm keeping this bug. DO NOT USE FOR ANTYHING REAL. """
function gnk(n,k)
  xy = rand(2,n)
  T = KDTree(xy)
  idxs = knn(T, xy, k)[1]
  # form the edges for sparse
  ei = Int[]
  ej = Int[]
  for i=1:n
    for j=idxs[i]
      if i > j
        push!(ei,i)
        push!(ej,j)
      end
    end
  end
  return xy, ei, ej
end
#Random.seed!(3)
#xy, ei, ej = gnk(25,5)
Random.seed!(1)
xy, ei, ej = gnk(15,4)
gr(size=(300,300))
graphplot(ei, ej, x =xy[1,:], y=xy[2,:],
  markercolor=:black, markerstrokecolor=:white,
  markersize=4, linecolor=1, linealpha=0.8, linewidth=0.7,
  axis_buffer=0.02, background=nothing)

##
for i=1:size(xy,2)
  annotate!(xy[1,i],xy[2,i],"$i")
end
gui()
## We are going to remove edge 14,15 for simplicity
n = size(xy,2)
A = sparse(ei,ej,1,n,n)
A = A + A' # make it symmetric
A[14,15] = 0
A[15,14] = 0
A[6,7] = 1
A[7,6] = 1
dropzeros!(A)
## Remove vertex 4 as it's disconnected
verts = collect(1:15)
deleteat!(verts,4)
A = A[verts,verts]
xy = xy[:,verts]
##
ei,ej = findnz(triu(A))[1:2]
graphplot(ei, ej, x =xy[1,:], y=xy[2,:],
  markercolor=:black, markerstrokecolor=:white,
  markersize=4, linecolor=1, linealpha=0.8, linewidth=0.7,
  axis_buffer=0.02, background=nothing)
##
using DelimitedFiles
writedlm("undir-graph-14.xy", xy')
writedlm("undir-graph-14.edges", [ei ej])
## Save a few figures
pyplot(size=(300,300))
p = graphplot(ei, ej, x =xy[1,:], y=xy[2,:],
  markercolor=:white, markerstrokecolor=:black,
  markersize=8, linecolor=3, linealpha=1.0, linewidth=1.0,
  axis_buffer=0.02, background=:transparent,dpi=300)
savefig("undir-graph-14.pdf")

## Save a few figures
pyplot(size=(300,300))
p = graphplot(ei, ej, x =xy[1,:], y=xy[2,:],
  markercolor=:white, markerstrokecolor=:black,
  markersize=14, linecolor=3, linealpha=1.0, linewidth=1.0,
  axis_buffer=0.02, background=:transparent,dpi=300)
for i=1:size(A,1)
  annotate!(xy[1,i],xy[2,i],text("$i",9,:black))
end
savefig("undir-graph-14-labels.pdf")


text:
"""Transform anchor based upon a rotation of a surface of size w x h."""
function  transform_anchor(ax, ay, w, h, angle)

    θ = -deg2rad(angle)

    sinθ = sin(θ)
    cosθ = cos(θ)

    # Dims of the transformed rect
    tw = abs(w * cosθ) + abs(h * sinθ)
    th = abs(w * sinθ) + abs(h * cosθ)

    # Offset of the anchor from the center
    cax = ax - w * 0.5
    cay = ay - h * 0.5

    # Rotated offset of the anchor from the center
    rax = cax * costheta - cay * sintheta
    ray = cax * sintheta + cay * costheta

    return (
        tw * 0.5 + rax,
        th * 0.5 + ray
    )
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "libcleri"
version = v"0.12.1"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/transceptor-technology/libcleri.git", "a8a1adc9dcf4e889a4d44c17acb20d74d0c6cfe5")
]

# Bash recipe for building across all platforms
script = raw"""
export CFLAGS=$(pcre2-config --cflags)
export LDFLAGS=$(pcre2-config --libs8)
cd $WORKSPACE/srcdir
cd libcleri/
make -C Release FN="libcleri.${dlext}"
make -C Release FN="libcleri.${dlext}" install INSTALL_PATH=${prefix}
if [[ "${target}" == *-apple-* ]]; then
    install_name_tool -id @rpath/libcleri.dylib.0 ${prefix}/lib/libcleri.dylib
fi
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms(; experimental=true)

# The products that we will ensure are always built
products = [
    LibraryProduct("libcleri", :libcleri)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="PCRE2_jll", uuid="efcefdf7-47ab-520b-bdef-62a2eaa19f15"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6")

max_stars_repo_name: jameskermode/Yggdrasil

max_issues_repo_name: jameskermode/Yggdrasil

max_forks_repo_name: jameskermode/Yggdrasil


text:
@reexport module Foos

export Foo, phrase, EnglishFoo, GermanFoo

abstract type Foo end

phrase(::Foo) = error("abstract method")

struct EnglishFoo <: Foo end

phrase(::EnglishFoo) = "Hi everybody! I am "

struct GermanFoo <: Foo end

phrase(::GermanFoo) = "Guten Tag! Ich bin "

end # module


text:
function cal_flux_residual(loc_globaldata, loc_ghost_holder, loc_globaldata_mutable, configData)
	phi_i = zeros(Float64,4)
	phi_k = zeros(Float64,4)
	G_i = zeros(Float64,4)
    G_k = zeros(Float64,4)
	result = zeros(Float64,4)
	qtilde_i = zeros(Float64,4)
	qtilde_k = zeros(Float64,4)
	dist_length = length(loc_globaldata)

	# updateLocalGhost(loc_ghost_holder, globaldata)
	# temp = CuArray(loc_ghost_holder)
	for idx in 1:dist_length
		if loc_globaldata[idx].flag_1 == 0
			wallindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, loc_globaldata_mutable, configData, idx, phi_i, phi_k)
		elseif loc_globaldata[idx].flag_1 == 2
			outerindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, loc_globaldata_mutable, configData, idx, phi_i, phi_k)
		elseif loc_globaldata[idx].flag_1 == 1
			interiorindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, loc_globaldata_mutable, configData, idx, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k)
		end
	end
	# Array(temp)
	return nothing
end

function wallindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, loc_globaldata_mutable, configData, itm, phi_i, phi_k)
		Gxp = wall_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k)
		Gxn = wall_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k)
		Gyn = wall_dGy_neg(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k)
		# @. loc_globaldata[itm].flux_res = (Gxp + Gxn + Gyn) * 2
		@. loc_globaldata_mutable[5:8, itm] = (Gxp + Gxn + Gyn) * 2
		# if itm == 3
		# 	println(IOContext(stdout, :compact => false), Gxp)
		# 	println(IOContext(stdout, :compact => false), Gxp + Gxn)
		# 	println(IOContext(stdout, :compact => false), Gxp + Gxn + Gyn)
		# end
	return nothing
end

function outerindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, loc_globaldata_mutable, configData, itm, phi_i, phi_k)
	Gxp = outer_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k)
	Gxn = outer_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k)
	Gyp = outer_dGy_pos(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k)
	# @. loc_globaldata[itm].flux_res = Gxp + Gxn + Gyp
	@. loc_globaldata_mutable[5:8, itm] = Gxp + Gxn + Gyp
	return nothing
end

function interiorindices_flux_residual(loc_globaldata, loc_ghost_holder, dist_length, loc_globaldata_mutable, configData, itm, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k)
	# for itm in interiorindices
		Gxp = interior_dGx_pos(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k)
		Gxn = interior_dGx_neg(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k)
		Gyp = interior_dGy_pos(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k)
		Gyn = interior_dGy_neg(loc_globaldata, loc_ghost_holder, dist_length, itm, configData, phi_i, phi_k, G_i, G_k, result, qtilde_i, qtilde_k)
		# if itm == 1
		# 	println("=======")
		# 	println(IOContext(stdout, :compact => false), Gxp)
		# 	println(IOContext(stdout, :compact => false), Gxn + Gxp)
		# 	println(IOContext(stdout, :compact => false), Gxn + Gxp + Gyp)
		# 	println(IOContext(stdout, :compact => false), Gxn + Gxp + Gyp + Gyn)
		# 	println()
		# end
		# @. loc_globaldata[itm].flux_res = Gxp + Gxn + Gyp + Gyn
		@. loc_globaldata_mutable[5:8, itm] = Gxp + Gxn + Gyp + Gyn
	# end
	return nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 1683dfda-924d-48df-a512-3ae220545877
using Plots,SpectralVis

# ╔═╡ 15d67460-cc53-11eb-1d73-89e8d0bfb651
#using PackageCompiler;create_sysimage(:Plots, sysimage_path="sys_plots.so", precompile_execution_file="precompile_plots.jl")

# ╔═╡ 27c29b29-875a-48cc-9c3a-601429da7dc8
begin
	env = SPECENV
	html"""<style>main {max-width: 900px;}</style>"""
end

# ╔═╡ 576dfcef-d2de-4e4c-a8ad-4756730559f8
illum = D65_illuminant()

# ╔═╡ c291aa08-6965-463f-b5ff-5956a386361f
begin
	cfund = conefund(env.cmf)
	cfdef = shift_cone(cfund, :l, -5.0)
	#cfdef = deactivate_cone(cfund, :m)
	newcmf = cmfunc(cfdef)
	oldcmf = env.cmf
	oldr = matrixR(conefund(oldcmf))
	newr = matrixR(conefund(newcmf))
end

# ╔═╡ ecf55fa6-4aec-4fab-be55-1df606e1111b
#begin
#	plot(oldcmf.λ, [oldcmf.x, oldcmf.y, oldcmf.z], color = :black)
#	plot!(newcmf.λ, [newcmf.x, newcmf.y, newcmf.z], color = :green)
#end

# ╔═╡ 03c6c696-1ffe-444a-bd36-dd6684c2ef75
begin
	lev = 8
end

# ╔═╡ 7d169a7b-3554-46ad-a97e-9a7337a05cae
function deficiency(deftype::Symbol, c::Spectrum, illum::LSpec = D65_illuminant() , shift::Real = 0.0)
    cflms = conefund(SPECENV.cmf)
    if deftype == :protanopia
        #Acvd.l = zeros(length(Acvd.l))
		cvd = deactivate_cone(cflms, :l)
    elseif deftype == :deuteranopia
        cvd = deactivate_cone(cflms, :m)
    elseif deftype == :tritanopia
        cvd = deactivate_cone(cflms, :s)
    elseif deftype == :protanomaly
        cvd = shift_cone(cflms, :l, shift)
    elseif deftype == :deuteranomaly
        cvd = shift_cone(cflms, :m, shift)
    elseif deftype == :tritanomaly
        cvd = shift_cone(cflms, :s, shift)
    else
        DomainError(deftype, "Color anomaly type does not exist.")
    end
	#Acvd = matrixA(cvd)
	#col = adapt_spec(c, env, env.ex, :cubic, :natural)
	#Clms = fundamental_metamer(D65_illuminant(), col, conefund(env.cmf))
	#cf = conefund(env.cmf)
	#Rlms = matrix_R(cf)
	#t = typeof(cf)
	#Rcvd = matrix_R(cvd)
	cvd
end

# ╔═╡ b1822d4e-6e83-4044-87da-ba2d578e5a7c
begin
	defi = :protanopia
	col = adapt_spec(agfait(113), env, :linear, :cubic, :natural)
	cvd = deficiency(defi, col, D65_illuminant())
	acvd = matrixA(cvd)
	cf = conefund(env.cmf)
	clms = fundamental_metamer(D65_illuminant(), col, cf)
	rlms = matrixR(cf)
	t = typeof(cf)
	rcvd = matrixR(cvd, defi)
end

# ╔═╡ 0b3183e7-a866-4632-b79e-3022fdfeff0b
begin
	#p1 = contour(oldr, ratio = 1.0, levels = lev, fill = 0, zlimits = (-0.01,0.025))
	p2 = wireframe(collect(1:441), collect(1:441),rcvd, ratio = 1.0)
	plot(p2, size = (900,600))
end

# ╔═╡ 36e9c976-ac66-4330-9ce1-051281980eaf
begin
	match = env.cmf
	Acvd = conefund(match)
	typeof(Acvd)(Acvd.λ, zeros(length(Acvd.λ)), Acvd.m, Acvd.s)

end

# ╔═╡ Cell order:
# ╠═15d67460-cc53-11eb-1d73-89e8d0bfb651
# ╠═1683dfda-924d-48df-a512-3ae220545877
# ╠═27c29b29-875a-48cc-9c3a-601429da7dc8
# ╠═576dfcef-d2de-4e4c-a8ad-4756730559f8
# ╠═c291aa08-6965-463f-b5ff-5956a386361f
# ╠═ecf55fa6-4aec-4fab-be55-1df606e1111b
# ╠═03c6c696-1ffe-444a-bd36-dd6684c2ef75
# ╠═0b3183e7-a866-4632-b79e-3022fdfeff0b
# ╟─7d169a7b-3554-46ad-a97e-9a7337a05cae
# ╠═b1822d4e-6e83-4044-87da-ba2d578e5a7c
# ╠═36e9c976-ac66-4330-9ce1-051281980eaf


text:
using DiffusionProfile
using Test

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


text:
Lar = LinearAlgebraicRepresentation

""" Module for integration of polynomials over 3D volumes and surfaces """
function M(alpha::Int, beta::Int)::Float64
    a = 0
    for l=0:(alpha + 1)
        a += binomial(alpha+1,l) * (-1)^l/(l+beta+1)
    end
    return a/(alpha + 1)
end

""" The main integration routine """
function TT(
tau::Array{Float64,2},
alpha::Int, beta::Int, gamma::Int,
signedInt::Bool=false)
	vo,va,vb = tau[:,1],tau[:,2],tau[:,3]
	a = va - vo
	b = vb - vo
	s1 = 0.0
	for h=0:alpha
		for k=0:beta
			for m=0:gamma
				s2 = 0.0
				for i=0:h
					s3 = 0.0
					for j=0:k
						s4 = 0.0
						for l=0:m
							s4 += binomial(m,l) * a[3]^(m-l) * b[3]^l * M(
								h+k+m-i-j-l, i+j+l )
						end
						s3 += binomial(k,j) * a[2]^(k-j) * b[2]^j * s4
					end
					s2 += binomial(h,i) * a[1]^(h-i) * b[1]^i * s3;
				end
				s1 += binomial(alpha,h) * binomial(beta,k) * binomial(gamma,m) *
						vo[1]^(alpha-h) * vo[2]^(beta-k) * vo[3]^(gamma-m) * s2
			end
		end
	end
	c = cross(a,b)
	if signedInt == true
		return s1 * norm(c) * sign(c[3])
	else
		return s1 * norm(c)
	end
end

"""
	II(P::Lar.LAR, alpha::Int, beta::Int, gamma::Int, signedInt=false)

Basic integration function on 2D plane.

# Example  unit 3D triangle
```julia
julia> V = [0.0 1.0 0.0; 0.0 0.0 1.0; 0.0 0.0 0.0]
3×3 Array{Float64,2}:
 0.0  1.0  0.0
 0.0  0.0  1.0
 0.0  0.0  0.0

julia> FV = [[1,2,3]]
1-element Array{Array{Int64,1},1}:
 [1, 2, 3]

julia> P = V,FV
([0.0 1.0 0.0; 0.0 0.0 1.0; 0.0 0.0 0.0], Array{Int64,1}[[1, 2, 3]])

julia> Lar.II(P, 0,0,0)
0.5
```
"""
function II(
P::LAR,
alpha::Int, beta::Int, gamma::Int,
signedInt=false)::Float64
    w = 0
    V, FV = P
    if typeof(FV) == Array{Int64,2}
    	FV = [FV[:,k] for k=1:size(FV,2)]
    end
    for i=1:length(FV)
        tau = hcat([V[:,v] for v in FV[i]]...)
        if size(tau,2) == 3
        	term = TT(tau, alpha, beta, gamma, signedInt)
        	if signedInt
        		w += term
        	else
        		w += abs(term)
        	end
        elseif size(tau,2) > 3
        	println("ERROR: FV[$(i)] is not a triangle")
        else
        	println("ERROR: FV[$(i)] is degenerate")
        end
    end
    return w
end

"""
	III(P::Lar.LAR, alpha::Int, beta::Int, gamma::Int)::Float64

Basic integration function on 3D space.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]
3×4 Array{Float64,2}:
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]]
4-element Array{Array{Int64,1},1}:
 [1, 2, 4]
 [1, 3, 2]
 [4, 3, 1]
 [2, 3, 4]

julia> P = V,FV
([0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0],
Array{Int64,1}[[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]])

julia> Lar.III(P, 0,0,0)
0.16666666666666674
```
"""
function III(P::LAR, alpha::Int, beta::Int, gamma::Int)::Float64
    w = 0
    V, FV = P
    for i=1:length(FV)
        tau = hcat([V[:,v] for v in FV[i]]...)
        vo,va,vb = tau[:,1],tau[:,2],tau[:,3]
        a = va - vo
        b = vb - vo
        c = cross(a,b)
        w += c[1]/norm(c) * TT(tau, alpha+1, beta, gamma)
    end
    return w/(alpha + 1)
end

"""
	surface(P::Lar.LAR, signedInt::Bool=false)::Float64

`surface` integral on polyhedron `P`.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]
3×4 Array{Float64,2}:
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]]
4-element Array{Array{Int64,1},1}:
 [1, 2, 4]
 [1, 3, 2]
 [4, 3, 1]
 [2, 3, 4]

julia> P = V,FV
([0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0],
Array{Int64,1}[[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]])

julia> Lar.volume(P)
0.16666666666666674
```
"""
function surface(P::Lar.LAR, signedInt::Bool=false)::Float64
    return II(P, 0, 0, 0, signedInt)
end

"""
	volume(P::Lar.LAR)::Float64

`volume` integral on polyhedron `P`.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0]
3×4 Array{Float64,2}:
 0.0  1.0  0.0  0.0
 0.0  0.0  1.0  0.0
 0.0  0.0  0.0  1.0

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]]
4-element Array{Array{Int64,1},1}:
 [1, 2, 4]
 [1, 3, 2]
 [4, 3, 1]
 [2, 3, 4]

julia> P = V,FV
([0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0],
Array{Int64,1}[[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]])

julia> Lar.volume(P)
0.16666666666666674
```
"""
function volume(P::LAR)::Float64
    return III(P, 0, 0, 0)
end

#v,(vv,ev,fv,cv) = p.cuboidGrid((1,1,1),true)
#V = hcat([Array{Float64,1}(v[k,:]) for k=1:size(v,1)]...)
#FV = hcat([Array{Int64,1}(fv[k,:]+1) for k=1:size(fv,1)]...)
#EV = hcat([Array{Int64,1}(ev[k,:]+1) for k=1:size(ev,1)]...)
#model1 = Any[V,FV,EV]
#P = V,[FV[:,k] for k=1:size(FV,2)]
#surface(P,false)
#
#
#""" Surface integration """
#function surfIntegration(larModel)
#    V,FV,EV = model
#    FE = crossRelation(FV,EV)
#    if typeof(FV) == Array{Int64,2}
#    	FV = [FV[:,k] for k=1:size(FV,2)]
#    end
#    if typeof(V) == Array{Int64,2}
#    	if size(V,1) == 2
#    		V = vcat(V,zeros(1,size(V,2)))
#    	end
#    end
#    cochain = []
#    triangles = []
#    faceVertPairs = []
#	for face=1:length(FE)
#		push!(faceVertPairs, hcat([EV[:,e] for e in FE[face]]...))
#		row = [faceVertPairs[face][1] for k=1:length(FE[face])]
#		push!(triangles, vcat(faceVertPairs[face],row'))
#        P = V,triangles[face]
#        area = surface(P,false)
#        push!(cochain,area)
#    end
#    return cochain
#end

#    TODO: after having implemented ∂_3/∂_2
#def signedIntSurfIntegration(model,signedInt=False)
#    V,FV,EV = model
#    V = [v+[0.0] if len(v)==2 else v for v in V]
#    cochain = []
#    for face in FV
#        triangles = AA(C(AL)(face[0]))(TRANS([face[1-1],face[2]]))
#        P = V,triangles
#        area = Surface(P,signedInt)
#        cochain += [area]
#    return cochain
#

"""
	firstMoment(P::Lar.LAR)::Array{Float64,1}

First moments as terms of the Euler tensor. Remember that the integration algorithm is a boundary integration. Hence the model must be a boundary model. In this case, a 2-complex of triangles.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];

julia> P = V,FV;

julia> Lar.firstMoment(P)
3-element Array{Float64,1}:
 0.0416667
 0.0416667
 0.0416667
```
"""
function firstMoment(P::LAR)::Array{Float64,1}
    out = zeros(3)
    out[1] = III(P, 1, 0, 0)
    out[2] = III(P, 0, 1, 0)
    out[3] = III(P, 0, 0, 1)
    return out
end

"""
	secondMoment(P::Lar.LAR)::Array{Float64,1}

Second moments as terms of the Euler tensor.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];

julia> P = V,FV;

julia> Lar.secondMoment(P)
3-element Array{Float64,1}:
 0.0166667
 0.0166667
 0.0166667
```
"""
function secondMoment(P::LAR)::Array{Float64,1}
    out = zeros(3)
    out[1] = III(P, 2, 0, 0)
    out[2] = III(P, 0, 2, 0)
    out[3] = III(P, 0, 0, 2)
    return out
end

"""
	inertiaProduct(P::Lar.LAR)::Array{Float64,1}

Inertia products as terms of the Euler tensor.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];

julia> P = V,FV;

julia> Lar.inertiaProduct(P)
3-element Array{Float64,1}:
 0.00833333
 0.00833333
 0.00833333
```
"""
function inertiaProduct(P::LAR)::Array{Float64,1}
    out = zeros(3)
    out[1] = III(P, 0, 1, 1)
    out[2] = III(P, 1, 0, 1)
    out[3] = III(P, 1, 1, 0)
    return out
end

"""
	centroid(P::Lar.LAR)::Array{Float64,1}

Barycenter or `centroid` of polyhedron `P`.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];

julia> P = V,FV;

julia> Lar.centroid(P)
3-element Array{Float64,1}:
 0.25
 0.25
 0.25
```
"""
function centroid(P::LAR)::Array{Float64,1}
	return firstMoment(P)./volume(P)
end

"""
	inertiaMoment(P::Lar.LAR)::Array{Float64,1}

Inertia moments  of polyhedron `P`.

# Example # unit 3D tetrahedron
```julia
julia> V = [0.0 1.0 0.0 0.0; 0.0 0.0 1.0 0.0; 0.0 0.0 0.0 1.0];

julia> FV = [[1, 2, 4], [1, 3, 2], [4, 3, 1], [2, 3, 4]];

julia> P = V,FV;

julia> Lar.inertiaMoment(P)
3-element Array{Float64,1}:
 0.0333333
 0.0333333
 0.0333333
```
"""
function inertiaMoment(P::LAR)::Array{Float64,1}
    out = zeros(3)
    result = secondMoment(P)
    out[1] = result[2] + result[3]
    out[2] = result[3] + result[1]
    out[3] = result[1] + result[2]
    return out
end

function chainAreas(V::Array{Float64,2},EV::Array{Int64,2},chains::Array{Int64,2})
	FE = [chains[:,f] for f=1:size(chains,2)]
	return chainAreas(V,EV,FE)
end

""" Implementation using integr.jl """
function chainAreas(V::Array{Float64,2}, EV::Array{Int64,2},
				chains::Array{Array{Int64,1},1})
	if size(V,1) == 2
		V = vcat(V,zeros(1,size(V,2)))
	end
	pivots = [EV[:,abs(chain[1])][1] for chain in chains]
	out = zeros(length(pivots))
	for k=1:length(chains)
		area = 0
		triangles = [[] for h=1:length(chains[k])]
		for h=1:length(chains[k])
			edge = chains[k][h]
			v1,v2 = EV[:,abs(edge)]
			if sign(edge) == -1
				v1,v2=v2,v1
			end
			triangles[h] = Int[pivots[k],v1,v2]
		end
		P = V,hcat(triangles...)
		out[k] = Surface(P,true)
	end
	return out
end


text:
using Test
using JustSayIt

functions = [:email]

@testset "$(basename(@__FILE__))" begin
    @testset "Module Email available" begin
        @test @isdefined Email
    end
    @testset "function $f available" for f in functions
        @test f in names(Email; all=true)
    end;
end;


text:
# Find the smallest order of magnitude that is larger than xspan This is a
# little opaque because I want to avoid assuming the log function is defined
# over typeof(xspan)
function bounding_order_of_magnitude(xspan)
    DT = typeof(xspan)
    one_dt = one(DT)

    a = 1
    step = 1
    while xspan < 10.0^a * one_dt
        a -= step
    end

    b = 1
    step = 1
    while xspan > 10.0^b * one_dt
        b += step
    end

    while a + 1 < b
        c = div(a + b, 2)
        if xspan < 10.0^c * one_dt
            b = c
        else
            a = c
        end
    end

    return b
end

# Empty catchall
optimize_ticks() = {}

# Find some reasonable values for tick marks.
#
# This is basically Wilkinson's ad-hoc scoring method that tries to balance
# tight fit around the data, optimal number of ticks, and simple numbers.
#
# Args:
#   x_min: minimum value occuring in the data.
#   x_max: maximum value occuring in the data.
#
# Returns:
#   A Float64 vector containing tick marks.
#
function optimize_ticks{T}(x_min::T, x_max::T; extend_ticks::Bool=false)
    if x_min == x_max
        return [x_min], x_min - one(T), x_min + one(T)
    end

    # tick intervals and scores
    # TODO: these should perhaps be part of the theme
    const Q = {(1,1), (5, 0.9), (2, 0.7), (25, 0.5), (3, 0.2)}
    const n = length(Q)

    # number of ticks
    const k_min   = 2
    const k_max   = 10
    const k_ideal = 5

    # generalizing "order of magnitude"
    xspan = x_max - x_min
    z = bounding_order_of_magnitude(xspan)
    one_t = one(T)

    high_score = -Inf
    z_best = 0.0
    k_best = 0.0
    r_best = 0.0
    q_best = 0.0

    while k_max * 10.0^(z+1) * one_t > xspan
        for k in k_min:k_max
            for (q, qscore) in Q
                span = (k - 1) * q * 10.0^z * one_t
                if span < xspan
                    continue
                end

                r = ceil((x_max - span) / (q*10.0^z * one_t))
                while r*q*10.0^z * one_t <= x_min
                    has_zero = r <= 0 && abs(r) < k

                    # simplicity
                    s = has_zero ? 1.0 : 0.0

                    # granularity
                    g = 0 < k < 2k_ideal ? 1 - abs(k - k_ideal) / k_ideal : 0.0

                    # coverage
                    c = 1.5 * xspan/span

                    score = (1/4)g + (1/6)s + (1/3)c + (1/4)qscore

                    # strict limits on coverage
                    if span >= 2.0*xspan || span < xspan
                        score -= 1000
                    end

                    if score > high_score
                        (q_best, r_best, k_best, z_best) = (q, r, k, z)
                        high_score = score
                    end
                    r += 1
                end
            end
        end
        z -= 1
    end

    span = q_best * 10.0^z_best * one_t
    if extend_ticks
        S = Array(typeof(1.0 * one_t), int(3 * k_best))
        for i in 0:(3*k_best - 1)
            S[i+1] = (r_best + i - k_best) * span
        end
        viewmin, viewmax = S[k_best + 1], S[2 * k_best]
    else
        S = Array(typeof(1.0 * one_t), int(k_best))
        for i in 0:(k_best - 1)
            S[i+1] = (r_best + i) * span
        end
        viewmin, viewmax = S[1], S[end]
    end

    S, viewmin, viewmax
end

function optimize_ticks(x_min::Date, x_max::Date; extend_ticks::Bool=false)
    # This can be pretty simple. We are choosing ticks on one of three
    # scales: years, months, days.
    if year(x_max) - year(x_min) <= 1
        if year(x_max) == year(x_min) && month(x_max) - month(x_min) <= 1
            ticks = Date[]
            if x_max - x_min > days(7)
                # This will probably need to be smarter
                push!(ticks, x_min)
                while true
                    next_month = date(year(ticks[end]), month(ticks[end])) + month(1)
                    while ticks[end] + week(1) < next_month - days(2)
                        push!(ticks, ticks[end] + week(1))
                    end
                    push!(ticks, next_month)
                    if next_month >= x_max
                        break
                    end
                end
            else
                push!(ticks, x_min)
                while ticks[end] < x_max
                    push!(ticks, ticks[end] + day(1))
                end
            end

            viewmin, viewmax = ticks[1], ticks[end]
            ticks, viewmin, viewmax
        else
            ticks = Date[]
            push!(ticks, date(year(x_min), month(x_min)))
            while ticks[end] < x_max
                push!(ticks, ticks[end] + month(1))
            end
            viewmin, viewmax = ticks[1], ticks[end]

            ticks, x_min, x_max
        end
    else
        ticks, viewmin, viewmax =
            optimize_ticks(year(x_min), year(x_max), extend_ticks=extend_ticks)
        Date[date(y) for y in ticks], date(viewmin), date(viewmax)
    end
end


text:
################################################################################
#                   CORE DISPATCHVARIABLEREF METHOD EXTENSIONS
################################################################################
# Extend dispatch_variable_ref
function dispatch_variable_ref(model::InfiniteModel,
                               index::PointVariableIndex
                               )::PointVariableRef
    return PointVariableRef(model, index)
end

# Extend _add_data_object
function _add_data_object(model::InfiniteModel,
                          object::VariableData{<:PointVariable}
                          )::PointVariableIndex
    return MOIUC.add_item(model.point_vars, object)
end

# Extend _data_dictionary (type based)
function _data_dictionary(model::InfiniteModel, ::Type{PointVariable}
    )::MOIUC.CleverDict{PointVariableIndex, VariableData{PointVariable{GeneralVariableRef}}}
    return model.point_vars
end

# Extend _data_dictionary (reference based)
function _data_dictionary(vref::PointVariableRef
    )::MOIUC.CleverDict{PointVariableIndex, VariableData{PointVariable{GeneralVariableRef}}}
    return JuMP.owner_model(vref).point_vars
end

# Extend _data_object
function _data_object(vref::PointVariableRef
    )::VariableData{PointVariable{GeneralVariableRef}}
    object = get(_data_dictionary(vref), JuMP.index(vref), nothing)
    object === nothing && error("Invalid point variable reference, cannot find " *
                        "corresponding variable in the model. This is likely " *
                        "caused by using the reference of a deleted variable.")
    return object
end

# Extend _core_variable_object
function _core_variable_object(vref::PointVariableRef
    )::PointVariable{GeneralVariableRef}
    return _data_object(vref).variable
end

################################################################################
#                          DEFINTION HELPER METHODS
################################################################################
# Ensure parameter values match shape of parameter reference tuple stored in the
# infinite variable reference
function _check_tuple_shape(_error::Function,
                            ivref::Union{InfiniteVariableRef, DerivativeRef},
                            values::VectorTuple)::Nothing
    prefs = raw_parameter_refs(ivref)
    if !same_structure(prefs, values)
        _error("The dimensions and array formatting of the infinite parameter " *
               "values must match those of the parameter references for the " *
               "infinite variable/derivative.")
    end
    return
end

## Dispatch methods for checking the supports of parameters in point variable
# IndependentParameterRefs
function _check_element_support(_error::Function, prefs::Vector{IndependentParameterRef},
                                param_values::Vector{Float64},
                                counter::Int)::Int
    for pref in prefs
        if !supports_in_set(param_values[counter], infinite_set(pref))
            _error("Parameter values violate parameter bounds.")
        end
        counter += 1
    end
    return counter
end

# DependentParameterRefs
function _check_element_support(_error::Function, prefs::Vector{DependentParameterRef},
                                param_values::Vector{Float64},
                                counter::Int)::Int
    len = length(prefs)
    supp = reshape(param_values[counter:counter+len-1], len, 1)
    if !supports_in_set(supp, infinite_set(prefs))
        _error("Parameter values violate parameter bounds.")
    end
    return counter += len
end

# Used to ensure values don't violate parameter bounds
function _check_tuple_values(_error::Function,
                             ivref::Union{InfiniteVariableRef, DerivativeRef},
                             param_values::Vector{Float64})::Nothing
    raw_prefs = raw_parameter_refs(ivref)
    counter = 1
    for i in 1:size(raw_prefs, 1)
        prefs = map(e -> dispatch_variable_ref(e), raw_prefs[i, :])
        counter = _check_element_support(_error, prefs, param_values, counter)
    end
    return
end

# Update point variable info to consider the infinite variable
function _update_point_info(info::JuMP.VariableInfo,
                            ivref::Union{InfiniteVariableRef, DerivativeRef},
                            point::Vector{Float64})::JuMP.VariableInfo
    if JuMP.has_lower_bound(ivref) && !info.has_fix && !info.has_lb
        info = JuMP.VariableInfo(true, JuMP.lower_bound(ivref),
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 info.binary, info.integer)
    end
    if JuMP.has_upper_bound(ivref) && !info.has_fix && !info.has_ub
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 true, JuMP.upper_bound(ivref),
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 info.binary, info.integer)
    end
    if JuMP.is_fixed(ivref) && !info.has_fix  && !info.has_lb  && !info.has_ub
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 true, JuMP.fix_value(ivref),
                                 info.has_start, info.start,
                                 info.binary, info.integer)
    end
    if start_value_function(ivref) !== nothing && !info.has_start
        if _is_vector_start(ivref)
            start = start_value_function(ivref)(point)
        else
            prefs = raw_parameter_refs(ivref)
            start = start_value_function(ivref)(Tuple(point, prefs)...)
        end
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 true, Float64(start),
                                 info.binary, info.integer)
    end
    if JuMP.is_binary(ivref) && !info.integer
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 true, info.integer)
    end
    if JuMP.is_integer(ivref) && !info.binary
        info = JuMP.VariableInfo(info.has_lb, info.lower_bound,
                                 info.has_ub, info.upper_bound,
                                 info.has_fix, info.fixed_value,
                                 info.has_start, info.start,
                                 info.binary, true)
    end
    return info
end

# Define _make_variable (used by JuMP.build_variable)
function _make_variable(_error::Function, info::JuMP.VariableInfo, ::Type{Point};
                        infinite_variable_ref::Union{GeneralVariableRef,
                                                     Nothing} = nothing,
                        parameter_values::Union{Real, AbstractArray{<:Real},
                                                Tuple, Nothing} = nothing,
                        extra_kw_args...)::PointVariable{GeneralVariableRef}
    # check for unneeded keywords
    for (kwarg, _) in extra_kw_args
        _error("Keyword argument $kwarg is not for use with point variables.")
    end
    # ensure the needed arguments are given
    if parameter_values === nothing || infinite_variable_ref === nothing
        _error("Must specify the infinite variable/derivative and the values of its " *
               "infinite parameters")
    end
    # format tuple as VectorTuple
    pvalues = VectorTuple{Float64}(parameter_values)
    # check information and prepare format
    dispatch_ivref = dispatch_variable_ref(infinite_variable_ref)
    if !(dispatch_ivref isa Union{InfiniteVariableRef, DerivativeRef})
     _error("Expected an infinite variable/derivative reference dependency, but got a " *
            "variable reference of type $(typeof(dispatch_ivref)).")
    end
    _check_tuple_shape(_error, dispatch_ivref, pvalues)
    _check_tuple_values(_error, dispatch_ivref, pvalues.values)
    info = _update_point_info(info, dispatch_ivref, pvalues.values)
    # enforce parameter significant digits on the values
    prefs = parameter_list(dispatch_ivref)
    for i in eachindex(pvalues)
        pvalues[i] = round(pvalues[i], sigdigits = significant_digits(prefs[i]))
    end
    # make variable and return
    return PointVariable(_make_float_info(info), infinite_variable_ref,
                         pvalues.values)
end

## Dispatch methods for updating the supports of parameters in point variable
# IndependentParameterRefs
function _add_point_support(prefs::Vector{IndependentParameterRef},
                            param_values::Vector{Float64},
                            counter::Int)::Int
    for pref in prefs
        add_supports(pref, param_values[counter], check = false,
                     label = UserDefined)
        counter += 1
    end
    return counter
end

# DependentParameterRefs
function _add_point_support(prefs::Vector{DependentParameterRef},
                            param_values::Vector{Float64},
                            counter::Int)::Int
    len = length(prefs)
    supp = reshape(param_values[counter:counter+len-1], len, 1)
    add_supports(prefs, supp, check = false, label = UserDefined)
    return counter += len
end

# Used to add point variable support to parameter supports if necessary
function _update_param_supports(ivref::Union{InfiniteVariableRef, DerivativeRef},
                                param_values::Vector{Float64})::Nothing
    raw_prefs = raw_parameter_refs(ivref)
    model = JuMP.owner_model(ivref)
    counter = 1
    for i in 1:size(raw_prefs, 1)
        prefs = dispatch_variable_ref.(raw_prefs[i, :])
        counter = _add_point_support(prefs, param_values, counter)
    end
    return
end

# Used to update mapping infinite_to_points
function _update_infinite_point_mapping(pvref::PointVariableRef,
    ivref::Union{InfiniteVariableRef, DerivativeRef}
    )::Nothing
    push!(_point_variable_dependencies(ivref), JuMP.index(pvref))
    return
end

# Define _check_and_make_variable_ref (used by JuMP.add_variable)
function _check_and_make_variable_ref(model::InfiniteModel,
                                      v::PointVariable,
                                      name::String;
                                      add_support = true)::PointVariableRef
    ivref = dispatch_variable_ref(v.infinite_variable_ref)
    JuMP.check_belongs_to_model(ivref, model)
    data_object = VariableData(v, name)
    vindex = _add_data_object(model, data_object)
    vref = PointVariableRef(model, vindex)
    if add_support
        _update_param_supports(ivref, v.parameter_values)
    end
    _update_infinite_point_mapping(vref, ivref)
    return vref
end

################################################################################
#                         PARAMETER VALUE METHODS
################################################################################

"""
    infinite_variable_ref(vref::PointVariableRef)::GeneralVariableRef

Return the `InfiniteVariableRef` associated with the point variable `vref`.

**Example**
```julia-repl
julia> @infinite_variable(model, T(t))
T(t)

julia> vref = @point_variable(model, T(0))
T(0)

julia> infinite_variable_ref(vref)
T(t)
```
"""
function infinite_variable_ref(vref::PointVariableRef)::GeneralVariableRef
    return _core_variable_object(vref).infinite_variable_ref
end

"""
    raw_parameter_values(vref::PointVariableRef)::Vector{Float64}

Return the raw support point values associated with the point variable `vref`.
"""
function raw_parameter_values(vref::PointVariableRef)::Vector{Float64}
    return _core_variable_object(vref).parameter_values
end

"""
    parameter_values(vref::PointVariableRef)::Tuple

Return the support point associated with the point variable `vref`.

**Example**
```julia-repl
julia> @infinite_variable(model, T(t))
T(t)

julia> vref = @point_variable(model, T(0))
T(0)

julia> parameter_values(vref)
(0,)
```
"""
function parameter_values(vref::PointVariableRef)::Tuple
    prefs = raw_parameter_refs(infinite_variable_ref(vref))
    return Tuple(VectorTuple(raw_parameter_values(vref), prefs.ranges,
                             prefs.indices))
end

# Internal function used to change the parameter value tuple of a point variable
function _update_variable_param_values(vref::PointVariableRef,
                                       pref_vals::Vector{<:Real})::Nothing
    info = _variable_info(vref)
    ivref = infinite_variable_ref(vref)
    new_var = PointVariable(info, ivref, pref_vals)
    _set_core_variable_object(vref, new_var)
    return
end

################################################################################
#                           VARIABLE INFO METHODS
################################################################################
# Set info for point variables
function _update_variable_info(vref::PointVariableRef,
                               info::JuMP.VariableInfo)::Nothing
    ivref = infinite_variable_ref(vref)
    param_values = raw_parameter_values(vref)
    new_var = PointVariable(info, ivref, param_values)
    _set_core_variable_object(vref, new_var)
    return
end

################################################################################
#                                 DELETION
################################################################################
# Extend _delete_variable_dependencies (for use with JuMP.delete)
function _delete_variable_dependencies(vref::PointVariableRef)::Nothing
    # remove variable info constraints associated with vref
    _delete_info_constraints(vref)
    # remove the infinite variable dependency
    ivref = infinite_variable_ref(vref)
    filter!(e -> e != JuMP.index(vref), _point_variable_dependencies(ivref))
    return
end


text:
nom = ARGS[1]

f = open(nom)

lines = readlines(f)

data = []

for line in lines
    liste = split(line)
    for elt in liste
        push!(data, parse(Float64,elt))
    end
end

#print(data)

close(f)

nbClients = round(Int64, data[1])
nbServices = round(Int64, data[2])

couts = []
capa = []

demandes = []

# matrice client*services
liaisons = rand(nbClients, nbServices)

#récupération des coûts de liaison
indice = 3
for i = 1:nbClients
    for j = 1:nbServices
        liaisons[i,j] = data[indice]
        indice += 1
    end
end

#récupération des demandes
for i = 1:nbClients
    push!(demandes, data[indice])
    indice += 1
end

#récupération des coûts d'ouverture
for i = 1:nbServices
    push!(couts, data[indice])
    indice += 1
end

#récupération des capacités des services
for i = 1:nbServices
    push!(capa, data[indice])
    indice += 1
end

#enregistrement dans un fichier
nom2 = split(nom, ".")[1]
nom2 = nom2*".dat"
sortie = open(nom2, "w")

print(sortie, nbServices)
print(sortie, " ")
println(sortie, nbClients)

# enregistrement capacité/coûts pour les services
for i = 1:nbServices
    print(sortie, capa[i])
    print(sortie, " ")
    println(sortie, couts[i])
end

# enregistrement des demandes
for i = 1:nbClients
    print(sortie, demandes[i])
    if i < nbClients
        print(sortie, " ")
    else
        print(sortie, "\n")
    end
end

# enregistrement des coûts de liaison
for i = 1:nbServices
    for j = 1:nbClients
        print(sortie, liaisons[j,i])
        if j < nbClients
            print(sortie, " ")
        else
            print(sortie, "\n")
        end
    end
end

close(sortie)

max_stars_repo_name: smbct/IP-FLP

max_issues_repo_name: smbct/IP-FLP

max_forks_repo_name: smbct/IP-FLP


text:
using ModelingToolkit

@variables t
sts = @variables x1(t) x2(t) x3(t) x4(t)
params = @parameters u1(t) u2(t) u3(t) u4(t)
D = Differential(t)
eqs = [
    x1 + x2 + u1 ~ 0
    x1 + x2 + x3 + u2 ~ 0
    x1 + D(x3) + x4 + u3 ~ 0
    2*D(D(x1)) + D(D(x2)) + D(D(x3)) + D(x4) + u4 ~ 0
]
@named sys = ODESystem(eqs, t)

let pss = partial_state_selection(sys)
    @test length(equations(pss)) == 1
    @test length(states(pss)) == 2
    @test length(equations(ode_order_lowering(pss))) == 2
end

@parameters t σ ρ β
@variables x(t) y(t) z(t) a(t) u(t) F(t)
D = Differential(t)

eqs = [
       D(x) ~ σ*(y-x)
       D(y) ~ x*(ρ-z)-y + β
       0 ~ z - x + y
       0 ~ a + z
       u ~ z + a
    ]

lorenz1 = ODESystem(eqs,t,name=:lorenz1)
let al1 = alias_elimination(lorenz1)
    let lss = partial_state_selection(al1)
        @test length(equations(lss)) == 2
    end
end

max_stars_repo_name: paulrschrater/ModelingToolsJulia

max_issues_repo_name: paulrschrater/ModelingToolsJulia

max_forks_repo_name: paulrschrater/ModelingToolsJulia


text:
# EXPERIMENTAL - Threaded version 2. Try fewer threads, one per column of the return

import Base.Threads.@spawn

corkendallthreads_v2(x::Union{RealVector,RealOrMissingVector}, y::Union{RealVector,RealOrMissingVector}) = corkendall(float(copy(x)), float(copy(y)))# threads not used in this case

function corkendallthreads_v2(X::Union{RealMatrix,RealOrMissingMatrix}, y::Union{RealVector,RealOrMissingVector})
    n = size(X, 2)
    C = ones(float(eltype(X)), n)

    tasks = Array{Task,1}(undef, n)

    permy = sortperm(y)
    for i = 1:n
        tasks[i] = @spawn ck!(float(copy(y)), float(X[:,i]), permy)
    end

    for i = 1:n
        C[i] = fetch(tasks[i])
    end

    return C
end

function corkendallthreads_v2(x::Union{RealVector,RealOrMissingVector}, Y::Union{RealMatrix,RealOrMissingMatrix})
    n = size(Y, 2)
    C = ones(float(eltype(Y)), 1, n)

    tasks = Array{Task,1}(undef, n)

    permx = sortperm(x)
    for i = 1:n
        tasks[i] = @spawn ck!(float(copy(x)), float(Y[:,i]), permx)
    end

    for i = 1:n
        C[1,i] = fetch(tasks[i])
    end

    return C
end

function corkendallthreads_v2(X::Union{RealMatrix,RealOrMissingMatrix})
    n = size(X, 2)
    C = ones(float(eltype(X)), n, n)# avoids dependency on LinearAlgebra

    tasks = Array{Task,2}(undef, n, n)

    for j = 2:n
        permx = sortperm(X[:,j])
        for i = 1:j - 1
            tasks[j,i] = @spawn ck!(X[:,j], X[:,i], permx)
        end
    end

    for j = 2:n
        for i = 1:j - 1
            C[j,i] = fetch(tasks[j,i])
            C[i,j] = C[j,i]
        end
    end

    return C
end

function corkendallthreads_v2(X::Union{RealMatrix,RealOrMissingMatrix}, Y::Union{RealMatrix,RealOrMissingMatrix})
    nr = size(X, 2)
    nc = size(Y, 2)
    C = zeros(float(eltype(X)), nr, nc)
    tasks = Array{Task,1}(undef, nc)
    for j = 1:nc
        tasks[j] = @spawn corkendall(X, Y[:,j])
    end

    for j = 1:nc
        C[:,j] = fetch(tasks[j])
    end

    return C
end


text:
module TeraSort

using Elly
using StatsBase

export teragen, genrandfile
export terasort, sortrandfile
export terasample, samplerandfile
export cleanup, startcluster, stopcluster

include("util.jl")
include("gen.jl")
include("sort.jl")
include("merge.jl")
include("sample.jl")

end # module


text:
using Test
using MathOptInterface
const MOI = MathOptInterface

@testset "Functions" begin
    x = MOI.VariableIndex(1)
    y = MOI.VariableIndex(2)
    z = MOI.VariableIndex(3)
    @testset "Broadcast" begin
        xf = MOI.SingleVariable(x)
        yf = MOI.SingleVariable(y)
        zf = MOI.SingleVariable(z)
        function sum_indices(sv1::MOI.SingleVariable, sv2::MOI.SingleVariable)
            return sv1.variable.value + sv2.variable.value
        end
        @test sum_indices.(xf, [yf, zf]) == [3, 4]
    end
    @testset "Copy" begin
        @testset "VectorOfVariables" begin
            f = MOI.VectorOfVariables([x, y])
            f_copy = copy(f)
            f_copy.variables[2] = z
            @test f.variables[2] == y
        end
    end
end

@testset "Base.convert" begin
    @testset "SingleVariable" begin
        model = MOI.Utilities.Model{Float64}()
        x = MOI.add_variable(model)
        f = MOI.SingleVariable(x)
        f_vov = convert(MOI.VectorOfVariables, f)
        @test f_vov ≈ MOI.VectorOfVariables([x])
        f_vaf = convert(MOI.VectorAffineFunction{Float64}, f)
        @test f_vaf ≈ MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            [0.0],
        )
        f_vqf = convert(MOI.VectorQuadraticFunction{Float64}, f)
        @test f_vqf ≈ MOI.VectorQuadraticFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(1.0, x))],
            MOI.VectorQuadraticTerm{Float64}[],
            [0.0],
        )
    end
    @testset "ScalarAffineFunction" begin
        model = MOI.Utilities.Model{Float64}()
        x = MOI.add_variable(model)
        f = MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(2.0, x)], 1.0)
        @test_throws(MethodError, convert(MOI.VectorOfVariables, f))
        f_vaf = convert(MOI.VectorAffineFunction{Float64}, f)
        @test f_vaf ≈ MOI.VectorAffineFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
            [1.0],
        )
        f_vqf = convert(MOI.VectorQuadraticFunction{Float64}, f)
        @test f_vqf ≈ MOI.VectorQuadraticFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
            MOI.VectorQuadraticTerm{Float64}[],
            [1.0],
        )
    end
    @testset "ScalarQuadraticFunction" begin
        model = MOI.Utilities.Model{Float64}()
        x = MOI.add_variable(model)
        f = MOI.ScalarQuadraticFunction(
            [MOI.ScalarAffineTerm(2.0, x)],
            [MOI.ScalarQuadraticTerm(3.0, x, x)],
            1.0,
        )
        @test_throws(MethodError, convert(MOI.VectorOfVariables, f))
        @test_throws(MethodError, convert(MOI.VectorAffineFunction{Float64}, f))
        f_vqf = convert(MOI.VectorQuadraticFunction{Float64}, f)
        @test f_vqf ≈ MOI.VectorQuadraticFunction(
            [MOI.VectorAffineTerm(1, MOI.ScalarAffineTerm(2.0, x))],
            MOI.VectorQuadraticTerm{Float64}[MOI.VectorQuadraticTerm(
                1,
                MOI.ScalarQuadraticTerm(3.0, x, x),
            )],
            [1.0],
        )
    end
end