text:
using CSFML.LibCSFML
using ModernGL
using Images, Colors

include("ModernGL_utils.jl")

# https://www.tomdalling.com/blog/modern-opengl/02-textures/

### Init

mode = sfVideoMode(700, 720, 32)
# settings = SF::ContextSettings.new(
#   depth: 24, stencil: 8, antialiasing: 4,
#   major: 3, minor: 0
# )
window = sfRenderWindow_create(mode, "SFML window", sfResize | sfClose, C_NULL)
@assert window != C_NULL
sfWindow_setVerticalSyncEnabled(window,true)
sfWindow_setActive(window, sfTrue)
println(createcontextinfo()) # used by ModernGL_utils;
glEnable(GL_DEPTH_TEST)
glDepthFunc(GL_LESS)
glClearColor(0.4 ,0.2, 0.3, 1.0)

### shaders

vsh = """
    $(get_glsl_version_string())
    //layout(location = 0) in vec3 v_position;
    //layout(location = 1) in vec3 v_color;

    in vec3 v_position;
    in vec2 v_texCoord;
    out vec2 fragTexCoord;

    //in vec3 v_color;
    //out vec3 initial_Color;
    void main() {
        fragTexCoord = v_texCoord;
        gl_Position = vec4(v_position, 1.0);
    }
    """

fsh = """
    $(get_glsl_version_string())
    // in vec3 initial_Color;

    in vec2 fragTexCoord;

    uniform sampler2D tex;
    out vec4 outColor;
    void main() {
        //outColor = vec4(0.5, 0.2, 0.3, 1.0);
        //outColor = vec4(frag_colour, 1.0 );
        outColor = texture(tex, fragTexCoord);
    }
    """

vertexShader = createShader(vsh, GL_VERTEX_SHADER)
fragmentShader = createShader(fsh, GL_FRAGMENT_SHADER)

function bindAttrib(shader_program)
    glBindAttribLocation(shader_program, 0, "v_position");
    glBindAttribLocation(shader_program, 1, "v_TextCoord");
end

program = createShaderProgram( bindAttrib, vertexShader, fragmentShader)
#glUseProgram(program)
#println("program : $program")
print("shader ---> ")
glCheckError()

### vertex, colors and buffers

vao = glGenVertexArray()
glBindVertexArray(vao)

points_vbo = glGenBuffer()
glBindBuffer(GL_ARRAY_BUFFER, points_vbo)

vertexData = GLfloat[
        #  X     Y     Z       U     V
         0.0, 0.8, 0.0,        0.5, 1.0,
        -0.8,-0.8, 0.0,        0.0, 0.0,
         0.8,-0.8, 0.0,        1.0, 0.0
]

glBufferData(GL_ARRAY_BUFFER, sizeof(vertexData), vertexData, GL_STATIC_DRAW)
println("points_vbo : $points_vbo")
glCheckError()

glEnableVertexAttribArray(0);
glVertexAttribPointer(0, 3, GL_FLOAT, GL_FALSE, 5*sizeof(GLfloat) , C_NULL)

glEnableVertexAttribArray(1)
glVertexAttribPointer(1, 2, GL_FLOAT, GL_TRUE, 5*sizeof(GLfloat) , C_NULL)

glBindBuffer(GL_ARRAY_BUFFER, 0)
glBindVertexArray(0)

### textures

#test = Array{Int32}(1)

# glGenTextures(n::GLsizei, textures::Ptr{GLuint})::Cvoid

minMagFiler = 0
wrapMode = 0

tex = GLint[0]
glGenTextures(1, tex)
glBindTexture(GL_TEXTURE_2D, tex[1]);
println("text ---> $tex")

## glTexParameteri(target::GLenum, pname::GLenum, param::GLint)::Cvoid

minMagFiler = GL_LINEAR
wrapMode = GL_CLAMP_TO_EDGE

glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, minMagFiler)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MAG_FILTER, minMagFiler)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_S, wrapMode)
glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_WRAP_T, wrapMode)

img = load("/julia/playgrounds/resources/hazard_256x256.png")
#img = rand(GLuint, 256, 256)
w, h = size(img)
rawTexData = []
for i in 1:w
    for j in 1:h
        push!(raw, RGBA(img[i,j]).r)
        push!(raw, RGBA(img[i,j]).g)
        push!(raw, RGBA(img[i,j]).b)
        #push!(raw, RGBA(img[i,j]).alpha)
    end
end

# glTexImage2D( target::GLenum,
                    # level::GLint,
                    # internalformat::GLint,
                    # width::GLsizei,
                    # height::GLsizei,
                    # border::GLint,
                    # format::GLenum,
                    # type_::GLenum,
                    # pixels::Ptr{Cvoid})
                    # (returns)::Cvoid

# glTexImage2D(GL_TEXTURE_2D,
#                  0,
#                  TextureFormatForBitmapFormat(bitmap.format()),
#                  (GLsizei)bitmap.width(),
#                  (GLsizei)bitmap.height(),
#                  0,
#                  TextureFormatForBitmapFormat(bitmap.format()),
#                  GL_UNSIGNED_BYTE,
#                  bitmap.pixelBuffer())

glTexImage2D(GL_TEXTURE_2D,
                0,
                GL_RGB,
                256,
                256,
                0,
                GL_RGB,
                GL_UNSIGNED_BYTE,
                rawTexData )
glCheckError()

glBindTexture(GL_TEXTURE_2D, 0)

###
event_ref = Ref{sfEvent}()
clock = sfClock_create()
# Let's express time in seconds.
function get_time(clock::Ptr{sfClock})
    sfTime_asSeconds(sfClock_getElapsedTime(clock))
end
function restart(clock::Ptr{sfClock})
     sfTime_asSeconds(sfClock_restart(clock))
end

target_frequency = 24.0
running = true
t0 = 9999999
while (running)

    while Bool(sfRenderWindow_isOpen(window)) && running
        frame_timestamp = get_time(clock)
        frame_time = frame_timestamp - t0
        global t0 = frame_timestamp
        actual_frequency = 1 / frame_time

        # process events
        while Bool(sfRenderWindow_pollEvent(window, event_ref))
            # close window : exit
            event_ref[].type == sfEvtClosed && sfRenderWindow_close(window)
            event_ref[].type == sfEvtResized && println("Trigger sfEvtResized.")
            event_ref[].type == sfEvtLostFocus && println("Trigger sfEvtLostFocus.")
            event_ref[].type == sfEvtGainedFocus && println("Trigger sfEvtGainedFocus.")
            event_ref[].type == sfEvtTextEntered && println("Trigger sfEvtTextEntered: $(event_ref[].text.unicode)")
            event_ref[].type == sfEvtKeyPressed && println("Trigger sfEvtKeyPressed: $(event_ref[].key.code)")
            event_ref[].type == sfEvtKeyReleased && println("Trigger sfEvtKeyReleased: $(event_ref[].key.code)")
            event_ref[].type == sfEvtMouseWheelMoved && println("Trigger sfEvtMouseWheelMoved: $(event_ref[].mouseWheel.x), $(event_ref[].mouseWheel.y)")
            event_ref[].type == sfEvtMouseWheelScrolled && println("Trigger sfEvtMouseWheelScrolled: $(event_ref[].mouseWheelScroll.wheel)")
            event_ref[].type == sfEvtMouseButtonPressed && println("Trigger sfEvtMouseButtonPressed: $(event_ref[].mouseButton.button)")
            event_ref[].type == sfEvtMouseButtonReleased && println("Trigger sfEvtMouseButtonReleased: $(event_ref[].mouseButton.x), $(event_ref[].mouseButton.y)")
            event_ref[].type == sfEvtMouseMoved && println("Trigger sfEvtMouseMoved: $(event_ref[].mouseMove.x), $(event_ref[].mouseMove.y)")

            if event_ref[].type == sfEvtClosed
                global running = false
                println("Render window closed.")
            elseif event_ref[].type == sfEvtResized
                glViewport(9, 9, event_ref[].size.width, event_ref[].size.height)
                println("Viewport updated.")
            elseif event_ref[].type == sfEvtKeyPressed
                if event_ref[].key.code == sfKeyQ
                    global running = false
                    println("Quiting application.")
                end
            end
        end
        time_event_processing = round(get_time(clock) - t0, digits=9)

        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)

        glUseProgram(program)

        glActiveTexture(GL_TEXTURE0)
        glBindTexture(GL_TEXTURE_2D, tex[1])
        # gProgram->setUniform("tex", 0);
        # void Program::setUniform(const GLchar* name, OGL_TYPE v0) \
        # glUniform1 ## TYPE_SUFFIX (uniform(name), v0);
        id = glGetUniformLocation(program,"tex")
        println("id ---> $id")
        glUniform1ui(id, 0)

        #glCheckError()

        glEnableVertexAttribArray(0)
        #glVertexAttribPointer(1, 2, GL_FLOAT, GL_TRUE,  5*sizeof(GLfloat), (const GLvoid*)(3 * sizeof(GLfloat)))

        glBindVertexArray(vao)
        glDrawArrays(GL_TRIANGLES, 0 , 3)

        glBindVertexArray(0)
        glUseProgram(0)

        time_render = round(get_time(clock) - time_event_processing - t0, digits=9)
        sfRenderWindow_display(window)
        time_to_sleep = 1 / target_frequency - time_render - time_event_processing

        if time_to_sleep > 0.005
            sleep(time_to_sleep-0.005)
        end
    end
end

sfRenderWindow_destroy(window)

max_stars_repo_path: GL_playground/tomdalling_openGL 0x02.jl

max_stars_repo_name: leroy-walton/playgrounds

max_issues_repo_path: GL_playground/tomdalling_openGL 0x02.jl

max_issues_repo_name: leroy-walton/playgrounds

max_forks_repo_path: GL_playground/tomdalling_openGL 0x02.jl

max_forks_repo_name: leroy-walton/playgrounds


text:
module BuiltInFunctions

    using ..Functions
    using ..HyperDimensional

    include("ConstantFunction.jl")
    include("VarFunction.jl")
    include("PowFunction.jl")
    include("ProductFunction.jl")
    include("SumFunction.jl")

    function _test()
        sum = ProductFunction{Float64, 2}(
                AbstractMathFun[
                    ConstFunction{Float64}(2.0), VarFunction{Float64}(1),
                    VarFunction{Float64}(2)])
        emit(sum, :TestSum)
    end

    function test()
        @assert (_test()(2.0, 3.0) == 12.0) "Incorrect Function Output!"
    end
end

ext: jl


text:
using MarketData
FactCheck.setstyle(:compact)
FactCheck.onlystats(true)

facts("In REPL plotting wrappers work for histogram") do

    h = uhist(cl);

    context("attributes are correct") do
        @fact h.title  --> "AAPL Close"
        @fact h.xlabel --> ""
        @fact h.ylabel --> ""
    end
end


text:
"""
    mutable struct AndersonUpdater <: AbstractFieldUpdater

    AndersonUpdater(; nhist = 50, iter_euler = 1, lam_euler = 0.01)

A field updating scheme implementing Anderson mixing, which obtains a new field configuration
based on a combination of previous configurations and residuals.

The number of configurations stored in history is set by the `nhist` argument.
A given number `iter_euler` of simple Euler update steps with size `lam_euler`
can be taken at the beginning of iteration to build of the necessary configuration histories.
"""
mutable struct AndersonUpdater <: AbstractFieldUpdater
    nhist      :: Int
    iter       :: Int
    iter_hist  :: Int
    iter_euler :: Int
    lam        :: Float64
    lam_euler  :: Float64
    setup      :: Bool

    res_hist   :: Vector{Dict{Int,FieldGrid{Float64}}} # Residual history for each field
    field_hist :: Vector{Dict{Int,FieldGrid{Float64}}} # Omega field history

    # Coefficient vectors, allocated at each step
    Umn        :: Matrix{Float64}
    Vm         :: Vector{Float64}
    Cn         :: Vector{Float64}

    # Step directions
    WW         :: FieldGrid{Float64}
    DD         :: FieldGrid{Float64}

    system     :: Option{FieldSystem}

    function AndersonUpdater(; nhist::Integer = 50, iter_euler::Integer = 1, lam_euler::Real = 0.01)
        # Must take at least 1 initial step before mixing
        if iter_euler <= 1; iter_euler = 1; end

        res_hist = [Dict() for i = 1:nhist+1]
        field_hist = [Dict() for i = 1:nhist+1]

        return new(nhist, 0, 0, iter_euler, 0.0, lam_euler, false, res_hist, field_hist,
            zeros(nhist, nhist), zeros(nhist), zeros(nhist), zeros(0,0,0), zeros(0,0,0),
             nothing
        )
    end
end

#==============================================================================#

Base.show(io::IO, updater::AndersonUpdater) = @printf(io, "AndersonUpdater(nhist = %d, iter = %d)", updater.nhist, updater.iter)

function setup!(updater::AndersonUpdater, sys::FieldSystem)
    updater.system = sys

    # If already bound to same system and setup, keep histories
    if updater.setup
        @warn "Using previous iterations in AndersonUpdater history."
    else
        if sys.ensemble != Canonical
            error("Anderson mixing can only be used for systems in the canonical ensemble.")
        end

        # Allocate step directions
        updater.WW = zeros(sys.dims)
        updater.DD = zeros(sys.dims)

        # Allocate history arrays
        for dict in updater.res_hist
            for (mid, res) in sys.residuals
                dict[mid] = zeros(sys.dims)
            end
        end
        for dict in updater.field_hist
            for (mid, omega) in sys.fields
                dict[mid] = zeros(sys.dims)
            end
        end

        # Store initial fields
        for (mid, omega) in sys.fields
            updater.field_hist[1][mid] .= omega
        end

        # Reset fields
        updater.iter = 0
        updater.iter_hist = 0
        updater.lam = 0.0
        updater.setup = true
    end

    return nothing
end

#==============================================================================#

function step!(updater::AndersonUpdater)
    @assert !isnothing(updater.system)
    sys = updater.system

    # Calculate lam and iter_hist
    updater.iter += 1
    if updater.iter < updater.nhist+1
        updater.iter_hist = updater.iter - 1
        updater.lam = 1.0 - 0.9^updater.iter
    else
        updater.iter_hist = updater.nhist
        updater.lam = 1.0
    end

    # Iteration numbers
    iter = updater.iter
    iter_hist = updater.iter_hist

    # Update density and residuals
    residuals!(sys)

    # Update the residual history
    for (mid, res) in sys.residuals
        updater.res_hist[iter_hist+1][mid] .= res
    end

    # First history iteration
    if iter <= updater.iter_euler
        # Update with simple Euler step
        for (mid, omega) in sys.fields
            @. omega += updater.lam_euler * sys.residuals[mid]
        end
    # Second iteration and beyond
    else
        # Reallocate appropriately sized coefficient arrays if needed
        if iter < updater.nhist+2
            updater.Umn = zeros(iter_hist, iter_hist)
            updater.Vm = zeros(iter_hist)
            updater.Cn = zeros(iter_hist)
        else
            fill!(updater.Umn, 0.0)
            fill!(updater.Vm, 0.0)
            fill!(updater.Cn, 0.0)
        end

        # Calculate coefficients of Umn and Vm arrays
        res_iter = updater.res_hist[iter_hist+1]
        field_iter = updater.field_hist[iter_hist+1]

        for mm = 1:iter_hist
            res_mm = updater.res_hist[iter_hist+1-mm]

            for nn = mm:iter_hist
                res_nn = updater.res_hist[iter_hist+1-nn]

                # Compute difference coefficients of Umn matrix
                #   elem = (curr - mm) * (curr - nn)
                for (mid, res) in res_iter
                    rmm = res_mm[mid]
                    rnn = res_nn[mid]

                    elem = 0.0
                    @simd for i in eachindex(res)
                        @inbounds elem += (res[i] - rmm[i]) * (res[i] - rnn[i])
                    end

                    updater.Umn[mm, nn] += elem
                end
                updater.Umn[nn, mm] = updater.Umn[mm, nn]
            end

            # Compute coefficients of Vm matrix
            #   elem = (curr - mm) * curr
            for (mid, res) in res_iter
                rmm = res_mm[mid]

                elem = 0.0
                @simd for i in eachindex(res)
                    @inbounds elem += (res[i] - rmm[i]) * res[i]
                end

                updater.Vm[mm] += elem
            end
        end

        # Inverting Umn matrix
        try
            inv!(lu!(updater.Umn))
            mul!(updater.Cn, updater.Umn, updater.Vm)
        catch err
            error("Error inverting Umn matrix on iteration $(iter): $(err)")
        end

        # Update the fields
        for (mid, omega) in sys.fields
            updater.WW .= omega
            updater.DD .= res_iter[mid]

            for nn = 1:iter_hist
                res_nn = updater.res_hist[iter_hist+1-nn]
                field_nn = updater.field_hist[iter_hist+1-nn]

                @. updater.WW += updater.Cn[nn] * (field_nn[mid] - field_iter[mid])
                @. updater.DD += updater.Cn[nn] * (res_nn[mid] - res_iter[mid])
            end

            @. omega = updater.WW + updater.lam * updater.DD
        end
    end

    # Update the history arrays
    if iter < updater.nhist+1
        # If not full, add to the array end
        for (mid, omega) in sys.fields
            updater.field_hist[iter+1][mid] .= omega
        end
    else
        # If full, shift all down one slot
        for nh = 1:updater.nhist
            for mid in keys(sys.fields)
                updater.res_hist[nh][mid] .= updater.res_hist[nh+1][mid]
                updater.field_hist[nh][mid] .= updater.field_hist[nh+1][mid]
            end
        end
        for (mid, omega) in sys.fields
            updater.field_hist[end][mid] .= omega
        end
    end

    return nothing
end


text:
# OBABO (PILE) integrator.

# Ceriotti, M., Parrinello, M., Markland, T. E., & Manolopoulos, D. E. (2010).
# Efficient stochastic thermostatting of path integral molecular dynamics. The
# Journal of Chemical Physics, 133(12), 124104.

struct OBABO{P,D,A,N} <: Integrator
    "Reciprocal temperature (mol / kJ)."
    beta::Float64
    nmt::NormalModeTransformer
    nms::NormalModeState
    oscprop::OscillatorPropagator
    forceapp::ForceApplicator
    thermostat::Maybe{Thermostat}
    "Instantaneous temperature for latest step (K)."
    temperature::Ref{Float64}
end

function OBABO(
        model::WaterModel, dt::Float64, beta::Float64, gamma0::Float64,
        fc::ForceContainer, s::State{P,D,A,N};
        extra_forces::Dict{String,<:Vector{<:AdHocForce}}=Dict{String,Vector{AdHocForce}}(),
        use_thermostat::Bool=true) where {P,D,A,N}
    Cmat = make_Cmat(P)
    omegas = [2P * sin(pi * k / P) / (hbar * beta) for k in 0:(P-1)] # 1 / ps

    nmt = NormalModeTransformer(Cmat)
    nms = NormalModeState(s)
    # A
    oscprop = OscillatorPropagator(dt, omegas)
    # B
    forceapp = ForceApplicator(model, extra_forces, 0.5 * dt, fc)
    # O
    if use_thermostat
        thermostat = Thermostat(beta, gamma0, 0.5 * dt, omegas)
    else
        thermostat = nothing
    end

    OBABO{P,D,A,N}(beta, nmt, nms, oscprop, forceapp, thermostat, NaN)
end

is_constrained(::Type{OBABO}) = false

function N_dof(::OBABO{P,D,A,N}) where {P,D,A,N}
    N * A * D * P
end

thermostat_enabled(int::OBABO) = !isnothing(int.thermostat)

function step!(s::State{P,D,A,N}, int::OBABO{P,D,A,N}) where {P,D,A,N}
    if thermostat_enabled(int)
        in_normal_modes(int.nms, s, int.nmt) do
            apply_thermostat!(int.nms, int.thermostat)
        end
    end

    apply_force!(s, int.forceapp; fresh_force=false)

    in_normal_modes(int.nms, s, int.nmt) do
        propagate_oscillators!(int.nms, int.oscprop)
    end

    apply_force!(s, int.forceapp)

    if thermostat_enabled(int)
        in_normal_modes(int.nms, s, int.nmt) do
            apply_thermostat!(int.nms, int.thermostat)
        end
    end
    int.temperature[] = temperature(MS_fict(P), N_dof(int), s)

    nothing
end

function energy(int::OBABO{P,D,A,N}, s::State{P,D,A,N}) where {P,D,A,N}
    result = Dict{String,Float64}()

    result["kinetic"] = kinetic_energy(MS_fict(P), s)
    result["spring"] = V_springs(int.beta, s)

    mergewith!(result, eval_pot(int.forceapp, s)) do
        error("Unexpected energy contribution.")
    end

    result
end

function get_force(int::OBABO)
    int.forceapp.current_force.initialized || error("Force not initialized.")

    int.forceapp.current_force.f
end

get_temperature(int::OBABO) = int.temperature[]

registered_integrators["OBABO"] = OBABO


text:
function ComplexEigenDecomposition(arg0::RealMatrix)
    return ComplexEigenDecomposition((RealMatrix,), arg0)
end

function get_d(obj::ComplexEigenDecomposition)
    return jcall(obj, "getD", FieldMatrix, ())
end

function get_determinant(obj::ComplexEigenDecomposition)
    return jcall(obj, "getDeterminant", jdouble, ())
end

function get_eigenvalues(obj::ComplexEigenDecomposition)
    return jcall(obj, "getEigenvalues", Vector{JComplex}, ())
end

function get_eigenvector(obj::ComplexEigenDecomposition, arg0::jint)
    return jcall(obj, "getEigenvector", FieldVector, (jint,), arg0)
end

function get_v(obj::ComplexEigenDecomposition)
    return jcall(obj, "getV", FieldMatrix, ())
end

function get_vt(obj::ComplexEigenDecomposition)
    return jcall(obj, "getVT", FieldMatrix, ())
end

function has_complex_eigenvalues(obj::ComplexEigenDecomposition)
    return jcall(obj, "hasComplexEigenvalues", jboolean, ())
end


text:
using SossMLJ
using Documenter

import Literate

include(joinpath(dirname(dirname(@__FILE__)), "test", "examples-list.jl"))

pages_before_examples = [
    "Home" => "index.md",
]
pages_examples = ["Examples" => ["$(example[1])" => "example-$(example[2]).md" for example in EXAMPLES]]
pages_after_examples = [
]
pages = vcat(
    pages_before_examples,
    pages_examples,
    pages_after_examples,
)

# Use Literate.jl to generate Markdown files for each of the examples
for example in EXAMPLES
    input_file = joinpath(EXAMPLESROOT, "example-$(example[2]).jl")
    Literate.markdown(input_file, DOCSOURCE)
end

makedocs(;
    modules=[SossMLJ],
    format=Documenter.HTML(),
    pages=pages,
    repo="https://github.com/cscherrer/SossMLJ.jl/blob/{commit}{path}#L{line}",
    sitename="SossMLJ.jl",
    authors="Chad Scherrer, Thibaut Lienart, Dilum Aluthge, Anthony Blaom, and contributors",
    strict=true,
)

deploydocs(;
    repo="github.com/cscherrer/SossMLJ.jl",
)


text:
using Flux, DiffEqFlux, OrdinaryDiffEq
using Zygote

u0 = Float32[2.; 0.]
datasize = 30
tspan = (0.0f0,1.5f0)

function trueODEfunc(du,u,p,t)
    true_A = [-0.1 2.0; -2.0 -0.1]
    du .= ((u.^3)'true_A)'
end
t = range(tspan[1],tspan[2],length=datasize)
prob = ODEProblem(trueODEfunc,u0,tspan)
ode_data = Array(solve(prob,Tsit5(),saveat=t))

dudt = Chain(x -> x.^3,
             Dense(2,50,tanh),
             Dense(50,2))
#= ps = Flux.params(dudt) =#

p = DiffEqFlux.destructure(dudt)
ps = Flux.params(p)
dudt_(u,p,t) = DiffEqFlux.restructure(dudt,p)(u)
prob = ODEProblem(dudt_,u0,tspan,p)

diffeq_adjoint(p,u0,prob,Tsit5())

forward(u->diffeq_adjoint(p,u,prob,Tsit5()),u0)

n_ode = x->neural_ode(dudt,x,tspan,Tsit5(),saveat=t,reltol=1e-7,abstol=1e-9)

pred = n_ode(u0) # Get the prediction using the correct initial condition

function predict_n_ode()
  n_ode(u0)
end
loss_n_ode() = sum(abs2,ode_data .- predict_n_ode())

data = Iterators.repeated((), 10)
opt = ADAM(0.1)
cb = function () #callback function to observe training
  display(loss_n_ode())
end

# Display the ODE with the initial parameter values.
cb()

Flux.train!(loss_n_ode, ps, data, opt, cb = cb)


text:
using DifferentialEquations

function rescale_jacobian!(dbi_dpj_n, db_dpj, dbi_dpi)
    for col in 1:size(dbi_dpj_n, 2)
        true_sum_of_non_diagonals = db_dpj[col] - dbi_dpi[col]
        actual_sum_of_non_diagonals = 0.0
        for row in 1:size(dbi_dpj_n, 1)
            if row != col
                actual_sum_of_non_diagonals += dbi_dpj_n[row, col]
            end
        end
        scale = true_sum_of_non_diagonals / actual_sum_of_non_diagonals
        for row in 1:size(dbi_dpj_n, 1)
            if row == col
                dbi_dpj_n[row, col] = dbi_dpi[col]
            else
                dbi_dpj_n[row, col] *= scale
            end
        end
    end
end

function integrated_explanation()
    cnt = 5
    d = rand(cnt)
    p = rand(cnt)

    b = total_burden(d, p)
    b_i = exact_burden(d, p)

    # Approximate values of the Jacobian.
    dbi_dpj_n = order_n_partial(d, p, cnt)
    # Approximate column sum of the Jacobian.
    db_dpj_n = sum(dbi_dpj_n, dims = 1)

    # The true column sum of the Jacobian.
    db_dpj = similar(p)
    partial_wrt_prevalence!(db_dpj, d, p)

    # The true diagonal of the Jacobian.
    dbi_dpi = b_i ./ p
    dbi_dpj_n
    rescale_jacobian!(dbi_dpj_n, db_dpj, dbi_dpi)
end

"""
t from 0 to 1.
p(t) = prevalence * t
dp = prevalence * dt
"""
struct BurdenPath
    d::Vector{Float64}
    prevalence::Vector{Float64}
    order::Int
end

"""
    (bp::BurdenPath)(db, b, params, t)

This is the derivative along the path from prevalence=0
at t=0 to full prevalence at t=1.
The arguments are the output derivatives, `db`, the input
total burden, `b` at time t, and parameters, which are unused.
"""
function (bp::BurdenPath)(db, b, params, t)
    if t > 0
        # Calulate an approximate Jacobian and then dot it with
        # the line integral direction, from (0,0,0...) to (p_1, p_2, p_3...).
        p = bp.prevalence * t
        # Approximate values of the Jacobian for first `db.order` terms.
        dbi_dpj_n = order_n_partial(bp.d, p, bp.order)

        # The true column sum of the Jacobian is easy-enough to calculate.
        db_dpj = similar(p)
        partial_wrt_prevalence!(db_dpj, bp.d, p)

        # The true diagonal of the Jacobian is very simple.
        dbi_dpi = b ./ p

        # We use the true values to fix the approximate values.
        rescale_jacobian!(dbi_dpj_n, db_dpj, dbi_dpi)
        db .= dbi_dpj_n * bp.prevalence
    else
        db .= bp.d .* bp.prevalence
    end
end

"""
    integrated_burden(weights, prevalences, order)

Calculates the burden for given weights and prevalences using integration.
The total sum of burden, by cause, is difficult to calculate because
individuals with many comorbidities contribute significantly to the total
burden. This means many of the combinatorial terms have to be included.
This function works around that problem by calculating the exact burden at
low prevalence and calculating the change in burden as prevalences increase.
In this way, the cross-terms, coming from comorbidities, contribute gradually
to the total burden, even if we don't explicitly calculate comorbidities for
individuals with many health outcomes.
"""
function integrated_burden(weights, prevalences, order)
    bp = BurdenPath(weights, prevalences, order)
    b0 = zeros(length(weights))
    tspan = (0.0, 1.0)
    prob = ODEProblem(bp, b0, tspan)
    # Tsit5(), BS3()
    sol = solve(prob, Tsit5(), save_everystep=false)
    sol[end]
end


text:
module DepthRenderer

using ModernGL
import GLFW
using Geometry: I4, TriangleMesh, CameraIntrinsics, SceneGraphNode, num_triangles

function scale_depth(x, near, far)
    far .* near ./ (far .- (far .- near) .* x)
end

function compute_projection_matrix(fx, fy, cx, cy, near, far, skew=0f0)
    proj = I4(Float32)
    proj[1, 1] = fx
    proj[2, 2] = fy
    proj[1, 2] = skew
    proj[1, 3] = -cx
    proj[2, 3] = -cy
    proj[3, 3] = near + far
    proj[3, 4] = near * far
    proj[4, 4] = 0.0f0
    proj[4, 3] = -1f0
    return proj
end

function compute_ortho_matrix(left, right, bottom, top, near, far)
    ortho = I4(Float32)
    ortho[1, 1] = 2f0 / (right-left)
    ortho[2, 2] = 2f0 / (top-bottom)
    ortho[3, 3] = - 2f0 / (far - near)
    ortho[1, 4] = - (right + left) / (right - left)
    ortho[2, 4] = - (top + bottom) / (top - bottom)
    ortho[3, 4] = - (far + near) / (far - near)
    return ortho
end

function perspective_matrix(width, height, fx, fy, cx, cy, near, far)
    # (height-cy) is used instead of cy because of the difference between
    # image coordinate systems between OpenCV and OpenGL. In the former,
    # the origin is at the top-left of the image while in the latter the
    # origin is at the bottom-left.
    proj_matrix = compute_projection_matrix(
            fx, fy, cx, (height-cy),
            near, far, 0.f0)
    ndc_matrix = compute_ortho_matrix(0, width, 0, height, near, far)
    ndc_matrix * proj_matrix
end

include("shaders.jl")

############
# Renderer #
############

mutable struct Renderer
    window::GLFW.Window
    cam::CameraIntrinsics
    compute_depth_shader::GLuint
    show_depth_shader::GLuint
    depth_texture::GLuint
    show_depth_vao::GLuint
    pos_attr::Int
    depth_image::Matrix{Float32}
    perspective_matrix::Matrix{Float32}
    mesh_vaos::Dict{TriangleMesh,GLuint}
end

function setup_for_show_in_window()

    # texture with which we'll show depth
    depth_texture = Ref(GLuint(0))
    glGenTextures(1, depth_texture)
    glBindTexture(GL_TEXTURE_2D, depth_texture[])
    glTexParameteri(GL_TEXTURE_2D, GL_TEXTURE_MIN_FILTER, GL_NEAREST)
    glBindTexture(GL_TEXTURE_2D, 0)

    # shader
    show_depth_shader, depth_pos_attr, depth_tex_attr = make_show_depth_shader()

    # vertices for showing depth image
    a = Float32[-1, -1, 0, 0]
    b = Float32[-1, 1, 0, 1]
    c = Float32[1, 1, 1, 1]
    d = Float32[1, -1, 1, 0]
    screen_vertices = hcat(a, b, c, a, c, d)
    show_depth_vao = Ref(GLuint(0))
    glGenVertexArrays(1, show_depth_vao)
    glBindVertexArray(show_depth_vao[])
    screen_vbo = Ref(GLuint(0))
    glGenBuffers(1, screen_vbo)
    glBindBuffer(GL_ARRAY_BUFFER, screen_vbo[])
    glBufferData(GL_ARRAY_BUFFER, sizeof(screen_vertices), Ref(screen_vertices, 1), GL_STATIC_DRAW)
    glVertexAttribPointer(
        depth_pos_attr, 2, GL_FLOAT, GL_FALSE,
        4 * sizeof(Float32), C_NULL)
    glEnableVertexAttribArray(depth_pos_attr)
    glVertexAttribPointer(
        depth_tex_attr, 2, GL_FLOAT, GL_FALSE,
        4 * sizeof(Float32), convert(Ptr{Cvoid}, 2 * sizeof(Float32)))
    glEnableVertexAttribArray(depth_tex_attr)
    glBindVertexArray(0)

    (depth_texture[]::GLuint, show_depth_vao[]::GLuint, show_depth_shader)
end

function Renderer(cam::CameraIntrinsics)

    # GLFW window
    window_hint = [
        (GLFW.SAMPLES,      0),
        (GLFW.DEPTH_BITS,   24),
        (GLFW.ALPHA_BITS,   8),
        (GLFW.RED_BITS,     8),
        (GLFW.GREEN_BITS,   8),
        (GLFW.BLUE_BITS,    8),
        (GLFW.STENCIL_BITS, 0),
        (GLFW.AUX_BUFFERS,  0),
        (GLFW.CONTEXT_VERSION_MAJOR, 3),
        (GLFW.CONTEXT_VERSION_MINOR, 3),
        (GLFW.OPENGL_PROFILE, GLFW.OPENGL_CORE_PROFILE),
        (GLFW.OPENGL_FORWARD_COMPAT, GL_TRUE),
    ]
    for (key, value) in window_hint
        GLFW.WindowHint(key, value)
    end
    window = GLFW.CreateWindow(cam.width, cam.height, "DepthRenderer")
    GLFW.MakeContextCurrent(window)

    # OpenGL setup for main rendering pipeline
    compute_depth_shader, pos_attr = make_compute_depth_shader()
    glEnable(GL_DEPTH_TEST)
    glViewport(0, 0, cam.width, cam.height)
    glClear(GL_DEPTH_BUFFER_BIT)

    # OpenGL setup for showing depth image in GLFW window
    (depth_texture, show_depth_vao, show_depth_shader) = setup_for_show_in_window()

    p = perspective_matrix(
        cam.width, cam.height,
        cam.fx, cam.fy,
        cam.cx, cam.cy,
        cam.near, cam.far)

    Renderer(
        window, cam, compute_depth_shader, show_depth_shader,
        depth_texture, show_depth_vao, pos_attr,
        Matrix{Float32}(undef, cam.width, cam.height), p,
        Dict{TriangleMesh,GLuint}())
end

function destroy!(r::Renderer)
    GLFW.DestroyWindow(r.window)
end

function get_depth_image!(renderer::Renderer; show_in_window=false)
    # depth buffer values in row major order
    data = Matrix{Float32}(undef, renderer.cam.width, renderer.cam.height)
    glReadPixels(0, 0, renderer.cam.width, renderer.cam.height, GL_DEPTH_COMPONENT, GL_FLOAT, Ref(data, 1))

    # compute actual depth in distance units from depth buffer values
    depth_image = scale_depth(data, renderer.cam.near, renderer.cam.far)

    if show_in_window
        glClearColor(0.0f0, 0.0f0, 0.0f0, 1.0f0)
        glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT)
        glActiveTexture(GL_TEXTURE0)
        glBindTexture(GL_TEXTURE_2D, renderer.depth_texture)
        glUseProgram(renderer.show_depth_shader)
        glUniform1i(glGetUniformLocation(renderer.show_depth_shader, "depth_texture"), 0)
        depth_normalized = depth_image ./ maximum(depth_image)
        glTexImage2D(
            GL_TEXTURE_2D, 0, GL_R32F, renderer.cam.width, renderer.cam.height,
            0, GL_RED, GL_FLOAT, Ref(depth_normalized, 1))
        glBindVertexArray(renderer.show_depth_vao)
        glDrawArrays(GL_TRIANGLES, 0, 6)
        glBindVertexArray(0)
        GLFW.SwapBuffers(renderer.window)
        GLFW.PollEvents()
    else
        glFlush()
    end
    glClear(GL_DEPTH_BUFFER_BIT)
    depth_image
end

function register_mesh!(renderer::Renderer, mesh::TriangleMesh)
    # vertices should be 3xN
    @assert size(mesh.vertices)[1] == 3

    # indices should be 3XM
    @assert size(mesh.indices)[1] == 3

    # TODO allow to change vertex data

    # create a vertex array object for the mesh and bind it
    vao = Ref(GLuint(0))
    glGenVertexArrays(1, vao)
    glBindVertexArray(vao[])

    # copy vertex data into an OpenGL buffer
    vbo = Ref(GLuint(0))
    glGenBuffers(1, vbo)
    glBindBuffer(GL_ARRAY_BUFFER, vbo[])
    glBufferData(GL_ARRAY_BUFFER, sizeof(mesh.vertices), Ref(mesh.vertices, 1), GL_STATIC_DRAW)

    # element buffer object for indices
    ebo = Ref(GLuint(0))
    glGenBuffers(1, ebo)
    glBindBuffer(GL_ELEMENT_ARRAY_BUFFER, ebo[])
    glBufferData(GL_ELEMENT_ARRAY_BUFFER, sizeof(mesh.indices), Ref(mesh.indices, 1), GL_STATIC_DRAW)

    # set vertex attribute pointers
    glVertexAttribPointer(renderer.pos_attr, 3, GL_FLOAT, GL_FALSE, 3 * sizeof(Float32), C_NULL)
    glEnableVertexAttribArray(renderer.pos_attr)

    # unbind it
    glBindVertexArray(0)

    renderer.mesh_vaos[mesh] = vao[]
    nothing
end

function draw!(renderer::Renderer, mesh::TriangleMesh, model::Matrix{Float32}, view::Matrix{Float32})
    if !haskey(renderer.mesh_vaos, mesh)
        error("mesh not registered")
    end
    vao = renderer.mesh_vaos[mesh]
    mvp = renderer.perspective_matrix * view * model
    glUseProgram(renderer.compute_depth_shader)
    glUniformMatrix4fv(0, 1, GL_FALSE, Ref(mvp, 1))
    glBindVertexArray(vao)
    glDrawElements(GL_TRIANGLES, num_triangles(mesh) * 3, GL_UNSIGNED_INT, C_NULL)
    glBindVertexArray(0)
end

######################
# draw a scene graph #
######################

function draw!(r::Renderer, node::SceneGraphNode, model::Matrix{Float32}, view::Matrix{Float32})
    model = model * node.transform
    if !isnothing(node.mesh)
        draw!(r, node.mesh, model, view)
    end
    for child in node.children
        draw!(r, child, model, view)
    end
end

export Renderer, register_mesh!, draw!, get_depth_image!, destroy!

###########
# File IO #
###########

import FileIO

function normalize_and_save_depth_image(depths::Matrix{Float64}, max_depth::Float64, fname::String)
    # clip and normalize depths
    depths = min.(depths, max_depth) ./ max_depth

    # save the file
    FileIO.save(fname, depths'[end:-1:1,:])
end

export normalize_and_save_depth_image

end # module Renderer

max_stars_repo_name: probcomp/DepthRenderer

max_issues_repo_name: probcomp/DepthRenderer

max_forks_repo_name: probcomp/DepthRenderer


text:
using LinearAlgebra
using DataStructures

const positions = map(line -> findall("#", line),
    readlines(length(ARGS) > 0 ? ARGS[1] : "example1.txt"))
const asteroids = Array{Vector,1}(undef,0)
for i in 1:length(positions)
    for j in 1:(length(positions[i]))
        push!(asteroids, Vector([positions[i][j][1] - 1, i - 1]))
    end
end

function direction(a, b)
    # Aha! Burned by a rounding error!
    v = [a[2] - b[2], b[1] - a[1]]
    normalized = v / norm(v)
    d = atand(normalized[2], normalized[1])
    # change to north-up azimuth increasing clockwise.
    #d = (450.0 - d) % 360.0
    d = (360 + d) % 360;
    return round(d, digits=3)
end

function directions(asteroids, a)
    return map(b -> direction(a, b), filter(b -> b != a, asteroids))
end

print("Day 10 Part 1: ");
v = findmax(map(a -> length(unique(directions(asteroids, a))), asteroids))
laser = asteroids[v[2]]
println("$(v[1]) are visible from $(laser)")

function print_field(a)
    xmin = findmin([x[1] for x in a])[1]
    xmax = findmax([x[1] for x in a])[1]
    ymin = findmin([y[1] for y in a])[1]
    ymax = findmax([y[1] for y in a])[1]
    for y in ymin:ymax
        for x in xmin:xmax
            print([x,y] in a ? '#' : '.')
        end
        println()
    end
end

#print_field(asteroids)

#for b in filter(x -> x != laser, asteroids)
#    println("angle from $(laser) to $(b) is $(direction(laser, b))")
#end

q = SortedDict()
for b in asteroids
    if b != laser
        d = direction(laser, b)
        if !(d in keys(q))
            q[d] = []
        end
        push!(q[d], b)
    end
end

# Well, that was ugly, but now we have a relation of angle -> [points].
# Next, we need to sort the points by their distance to the laser.
for target in values(q)
    sort!(target, by=(x -> norm(x - laser)))
end

const print_these = [1,2,3,10,20,50,100,199,200,201,299];
const th = Dict(1 => "st", 2 => "nd", 3 => "rd");

# Now spin the laser, remove the closest vertex each time, until we reach i == 200.
i = 0
while length(q) > 0
    for targets in keys(q)
        #println("   $(targets) => $(q[targets]): ");
        target = popfirst!(q[targets])
        if isempty(q[targets])
            delete!(q, targets)
        end
        global i += 1
        if i == 200
            global vaporize = target
        end
        if i in print_these
            println("The $(i)$(get(th, i, "th")) asteroid to be vaporized is at $(target[1]),$(target[2]).")
        end
    end
end

println("Day 10 Part 2: $(vaporize[1] * 100 + vaporize[2])")


text:
### Compute acceptance or rejection rate of MCChain

function acceptance(c::MCChain; lags::Range=1:size(c.samples, 1), reject::Bool=false)
  rlen = length(lags)
  @assert lags[end] <= size(c.samples, 1) "Range of acceptance rate not within post-burnin range of MCmc chain"

  if reject
    return (rlen-sum(c.diagnostics["accept"][lags]))*100/rlen
  else
    return sum(c.diagnostics["accept"][lags])*100/rlen
  end
end

### describe() provides summary statistics for MCChain objects

describe(c::MCChain) = describe(STDOUT, c)

function describe(io, c::MCChain)
  nsamples, npars = size(c.samples)
  for i in 1:npars

    col = c.samples[:, i]
    println(io, "Parameter $i")

    if sum(isnan(col)) != 0
      println(col, "Monte Carlo chain for parameter $i contains NaNs, which are not supported.")
      return
    end

    qs = quantile(col, [0, 1])
    varimse = mcvar_imse(col)
    variid = var(col)/nsamples
    mcerror = sqrt(varimse)
    ss = nsamples*variid./varimse
    act = varimse./variid

    statNames = ["Min", "Mean", "Max", "MC Error", "ESS", "AC Time"]
    statVals = [qs[1], mean(col), qs[2], mcerror, ss, act]
    for i = 1:6
        println(io, string(rpad(statNames[i], 10, " "), " ", string(statVals[i])))
    end

    println(io, )
  end
end


text:
# variables
""
function variable_transformer_voltage(pm::AbstractIVRModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    variable_transformer_voltage_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    variable_transformer_voltage_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    variable_transformer_voltage_excitation_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    variable_transformer_voltage_excitation_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end

""
function variable_transformer_current(pm::AbstractIVRModel; nw::Int=nw_id_default, bounded::Bool=true, report::Bool=true, kwargs...)
    variable_transformer_current_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    variable_transformer_current_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    variable_transformer_current_series_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    variable_transformer_current_series_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    variable_transformer_current_excitation_real(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    variable_transformer_current_excitation_imaginary(pm, nw=nw, bounded=bounded, report=report; kwargs...)

    expression_transformer_power(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    # expression_transformer_series_power(pm, nw=nw, bounded=bounded, report=report; kwargs...)
    expression_transformer_excitation_power(pm, nw=nw, bounded=bounded, report=report; kwargs...)
end

# constraints
""
function constraint_current_balance(pm::AbstractIVRModel, n::Int, i, bus_arcs, bus_arcs_dc, bus_arcs_xfmr, bus_gens, bus_loads, bus_gs, bus_bs)
    vr = var(pm, n, :vr, i)
    vi = var(pm, n, :vi, i)

    cr = var(pm, n, :cr)
    ci = var(pm, n, :ci)
    crt = var(pm, n, :crt)
    cit = var(pm, n, :cit)
    crdc = var(pm, n, :crdc)
    cidc = var(pm, n, :cidc)

    crg = var(pm, n, :crg)
    cig = var(pm, n, :cig)
    crd = var(pm, n, :crd)
    cid = var(pm, n, :cid)

    JuMP.@constraint(pm.model,  sum(cr[a] for a in bus_arcs)
                                + sum(crdc[d] for d in bus_arcs_dc)
                                + sum(crt[t] for t in bus_arcs_xfmr)
                                ==
                                sum(crg[g] for g in bus_gens)
                                - sum(crd[d] for d in bus_loads)
                                - sum(gs for gs in values(bus_gs))*vr + sum(bs for bs in values(bus_bs))*vi
                                )
    JuMP.@constraint(pm.model,  sum(ci[a] for a in bus_arcs)
                                + sum(cidc[d] for d in bus_arcs_dc)
                                + sum(cit[t] for t in bus_arcs_xfmr)
                                ==
                                sum(cig[g] for g in bus_gens)
                                - sum(cid[d] for d in bus_loads)
                                - sum(gs for gs in values(bus_gs))*vi - sum(bs for bs in values(bus_bs))*vr
                                )

end

""
function constraint_transformer_core_excitation(pm::AbstractIVRModel, n::Int, t)
    cert = var(pm, n, :cert, t)
    ceit = var(pm, n, :ceit, t)

    JuMP.@constraint(pm.model, cert == 0.0)
    JuMP.@constraint(pm.model, ceit == 0.0)
end

""
function constraint_transformer_core_excitation(pm::AbstractIVRModel, n::Int, t, int_a, int_b)
    cert = var(pm, n, :cert, t)
    ceit = var(pm, n, :ceit, t)

    voltage_harmonics_ntws = _PMs.ref(pm, n, :xfmr, t, "voltage_harmonics_ntws")

    et = reduce(vcat,[[var(pm, nw, :ert, t),var(pm, nw, :eit, t)]
                       for nw in voltage_harmonics_ntws])

    sym_exc_a = Symbol("exc_a_",n,"_",t)
    sym_exc_b = Symbol("exc_b_",n,"_",t)

    JuMP.register(pm.model, sym_exc_a, length(et), int_a; autodiff=true)
    JuMP.register(pm.model, sym_exc_b, length(et), int_b; autodiff=true)

    JuMP.add_NL_constraint(pm.model, :($(cert) == $(sym_exc_a)($(et...))))
    JuMP.add_NL_constraint(pm.model, :($(ceit) == $(sym_exc_b)($(et...))))
end

""
function constraint_transformer_core_voltage_drop(pm::AbstractIVRModel, n::Int, t, f_idx, xsc)
    ert = var(pm, n, :ert, t)
    eit = var(pm, n, :eit, t)

    vrt = var(pm, n, :vrt, f_idx)
    vit = var(pm, n, :vit, f_idx)

    csrt = var(pm, n, :csrt, f_idx)
    csit = var(pm, n, :csit, f_idx)

    JuMP.@constraint(pm.model, vrt == ert - xsc * csit)
    JuMP.@constraint(pm.model, vit == eit + xsc * csrt)
end

""
function constraint_transformer_core_voltage_balance(pm::AbstractIVRModel, n::Int, t, t_idx, tr, ti)
    ert = var(pm, n, :ert, t)
    eit = var(pm, n, :eit, t)

    vrt = var(pm, n, :vrt, t_idx)
    vit = var(pm, n, :vit, t_idx)

    JuMP.@constraint(pm.model, vrt == tr * ert - ti * eit)
    JuMP.@constraint(pm.model, vit == tr * eit + ti * ert)
end

""
function constraint_transformer_core_current_balance(pm::AbstractIVRModel, n::Int, t, f_idx, t_idx, tr, ti)
    cert = var(pm, n, :cert, t)
    ceit = var(pm, n, :ceit, t)

    csrt_fr = var(pm, n, :csrt, f_idx)
    csit_fr = var(pm, n, :csit, f_idx)

    csrt_to = var(pm, n, :csrt, t_idx)
    csit_to = var(pm, n, :csit, t_idx)

    JuMP.@constraint(pm.model, csrt_fr + tr * csrt_to + ti * csit_to
                                == cert
                    )
    JuMP.@constraint(pm.model, csit_fr + tr * csit_to - ti * csrt_to
                                == ceit
                    )
end

""
function constraint_transformer_winding_config(pm::AbstractIVRModel, n::Int, nh, i, idx, r, re, xe, gnd)
    vr = var(pm, n, :vr, i)
    vi = var(pm, n, :vi, i)

    vrt = var(pm, n, :vrt, idx)
    vit = var(pm, n, :vit, idx)

    crt = var(pm, n, :crt, idx)
    cit = var(pm, n, :cit, idx)

    # h ∈ 𝓗⁺ ⋃ 𝓗⁻
    if !is_zero_sequence(nh)
        JuMP.@constraint(pm.model, vrt == vr - r * crt)
        JuMP.@constraint(pm.model, vit == vi - r * cit)
    end

    # h ∈ 𝓗⁰, gnd == true -> cnf ∈ {Ye, Ze}
    if is_zero_sequence(nh) && gnd == 1
        JuMP.@constraint(pm.model, vrt == vr - (r + 3re) * crt + 3xe * cit)
        JuMP.@constraint(pm.model, vit == vi - (r + 3re) * cit - 3xe * crt)
    end

    # h ∈ 𝓗⁰, gnd == false -> cnf ∈ {D, Y, Z}
    if is_zero_sequence(nh) && gnd != 1
        JuMP.@constraint(pm.model, crt == 0)
        JuMP.@constraint(pm.model, cit == 0)
    end
end

""
function constraint_transformer_winding_current_balance(pm::AbstractIVRModel, n::Int, nh, idx, r, b_sh, g_sh, cnf)
    vrt = var(pm, n, :vrt, idx)
    vit = var(pm, n, :vit, idx)

    crt = var(pm, n, :crt, idx)
    cit = var(pm, n, :cit, idx)

    csrt = var(pm, n, :csrt, idx)
    csit = var(pm, n, :csit, idx)

    # h ∈ 𝓗⁺ ⋃ 𝓗⁻
    if !is_zero_sequence(nh)
        JuMP.@constraint(pm.model, crt == csrt - g_sh * vrt + b_sh * vit)
        JuMP.@constraint(pm.model, cit == csit - g_sh * vit - b_sh * vrt)
    end

    # h ∈ 𝓗⁰, cnf ∈ {Y(e), Z(e)}
    if is_zero_sequence(nh) && cnf in ['Y','Z']
        JuMP.@constraint(pm.model, crt == csrt - g_sh * vrt + b_sh * vit)
        JuMP.@constraint(pm.model, cit == csit - g_sh * vit - b_sh * vrt)
    end

    # h ∈ 𝓗⁰, cnf ∈ {D}
    if is_zero_sequence(nh) && cnf in ['D']
        JuMP.@constraint(pm.model, crt == csrt - g_sh * vrt + b_sh * vit - vrt / r)
        JuMP.@constraint(pm.model, crt == csit - g_sh * vit - b_sh * vrt - vit / r)
    end
end

""
function constraint_voltage_magnitude_rms(pm::AbstractIVRModel, i, vminrms, vmaxrms)
    w = [var(pm, nw, :w, i) for nw in _PMs.nw_ids(pm)]
    @assert vminrms>0
    @assert vmaxrms>vminrms

    JuMP.@constraint(pm.model, vminrms^2 <= sum(w)               )
    JuMP.@constraint(pm.model,              sum(w)  <= vmaxrms^2 )
end

""
function constraint_voltage_thd(pm::AbstractIVRModel, i, fundamental, thdmax)
    harmonics = Set(_PMs.nw_ids(pm))
    nonfundamentalharmonics = setdiff(harmonics, [fundamental])
    w = [var(pm, nw, :w, i) for nw in nonfundamentalharmonics]
    wfun = var(pm, fundamental, :w, i)

    JuMP.@constraint(pm.model, sum(w) <= thdmax^2*(wfun))
end

function constraint_voltage_harmonics_relative_magnitude(pm::AbstractIVRModel, n::Int, i, rm, fundamental)
    w  = var(pm, n, :w, i)
    wfun  = var(pm, fundamental, :w, i)

    JuMP.@constraint(pm.model, w <= rm*wfun)
end

function constraint_load_constant_power(pm::AbstractIVRModel, n::Int, i, bus, pd, qd)
    vr = var(pm, n, :vr, bus)
    vi = var(pm, n, :vi, bus)
    crd = var(pm, n, :crd, i)
    cid = var(pm, n, :cid, i)

    JuMP.@constraint(pm.model, pd == vr*crd  + vi*cid)
    JuMP.@constraint(pm.model, qd == vi*crd  - vr*cid)

end

function constraint_load_constant_current(pm::AbstractIVRModel, n::Int, i, bus, multiplier)
    crd = var(pm, n, :crd, i)
    cid = var(pm, n, :cid, i)
    fundamental = 1
    crd_fund = var(pm, fundamental, :crd, i)
    cid_fund = var(pm, fundamental, :cid, i)

    JuMP.@constraint(pm.model, crd == multiplier * crd_fund)
    JuMP.@constraint(pm.model, cid == multiplier * cid_fund)
end

function constraint_vm_auxiliary_variable(pm::AbstractIVRModel, n::Int, i)
    vr = var(pm, n, :vr, i)
    vi = var(pm, n, :vi, i)
    w  = var(pm, n, :w, i)
    #TODO choose whether this is relaxation or exact
    JuMP.@constraint(pm.model, w >= vr^2  + vi^2)
end

"reference bus angle constraint"
function constraint_ref_bus(pm::AbstractIVRModel, n::Int, i::Int)
    if n == 1 #fundamental frequency, fix reference angle
        JuMP.@constraint(pm.model, var(pm, n, :vi)[i] == 0.0)
        #TODO you should be able to set this more freely, but it helps a lot with stability, keeping it in for now.
        JuMP.@constraint(pm.model, var(pm, n, :vr)[i] == 1.0)
    else #fix harmonic voltage at reference bus to 0+j0
        JuMP.@constraint(pm.model, var(pm, n, :vi)[i] == 0.0)
        JuMP.@constraint(pm.model, var(pm, n, :vr)[i] == 0.0)
    end
end

function constraint_current_limit_rms(pm::AbstractIVRModel, f_idx, c_rating)
    (l, f_bus, t_bus) = f_idx
    t_idx = (l, t_bus, f_bus)

    crf =  [var(pm, n, :cr, f_idx) for n in _PMs.nw_ids(pm)]
    cif =  [var(pm, n, :ci, f_idx) for n in _PMs.nw_ids(pm)]

    crt =  [var(pm, n, :cr, t_idx) for n in _PMs.nw_ids(pm)]
    cit =  [var(pm, n, :ci, t_idx) for n in _PMs.nw_ids(pm)]

    JuMP.@constraint(pm.model, sum(crf^2 + cif^2) <= c_rating^2)
    JuMP.@constraint(pm.model, sum(crt^2 + cit^2) <= c_rating^2)
end

function constraint_active_filter(pm::AbstractIVRModel, i, fundamental)
    pgfun = var(pm, fundamental, :pg, i)
    pg =  [var(pm, n, :pg, i) for n in _PMs.nw_ids(pm)]

    JuMP.@NLconstraint(pm.model, pgfun == 0)
    JuMP.@NLconstraint(pm.model, sum(pg[n] for n in _PMs.nw_ids(pm)) == 0)
end

""
function objective_current_distortion_minimization(pm::AbstractIVRModel; gen_id=1, fundamental=1)
    harmonics = Set(_PMs.nw_ids(pm))
    nonfundamentalharmonics = setdiff(harmonics, [fundamental])

    crg = [var(pm, n, :crg, gen_id) for n in nonfundamentalharmonics]
    cig = [var(pm, n, :cig, gen_id) for n in nonfundamentalharmonics]

    pg = var(pm, fundamental, :pg, gen_id)

    JuMP.@NLconstraint(pm.model, pg <= 1.001*0.643386)

    #minimize magnitude of nonfundamental harmonics
    JuMP.@objective(pm.model, Min, sum(crg.^2 + cig.^2))
end

""
function objective_voltage_distortion_minimization(pm::AbstractIVRModel; bus_id=6, fundamental=1, gen_id=1)
    harmonics = Set(_PMs.nw_ids(pm))
    nonfundamentalharmonics = setdiff(harmonics, [fundamental])

    vr = [var(pm, n, :vr, bus_id) for n in nonfundamentalharmonics]
    vi = [var(pm, n, :vi, bus_id) for n in nonfundamentalharmonics]

    pg = var(pm, fundamental, :pg, gen_id)

    JuMP.@NLconstraint(pm.model, pg <= 1.001*0.643386)

    #minimize magnitude of nonfundamental harmonics
    JuMP.@objective(pm.model, Min, sum(vr.^2 + vi.^2))
end

#solution_processors=[sol_data_model!]
function sol_data_model!(pm::AbstractIVRModel, solution::Dict)
    _PMs.apply_pm!(_sol_data_model_ivr!, solution)
end

""
function _sol_data_model_ivr!(solution::Dict)
    if haskey(solution, "bus")
        for (i, bus) in solution["bus"]
            if haskey(bus, "vr") && haskey(bus, "vi")
                bus["vm"] = hypot(bus["vr"], bus["vi"])
                bus["va"] = atan(bus["vi"], bus["vr"])*180/pi
            end
        end
    end

    if haskey(solution, "branch")
        for (i, branch) in solution["branch"]
            if haskey(branch, "pf") && haskey(branch, "pt")
                branch["ploss"] = branch["pf"] + branch["pt"]
            end
            if haskey(branch, "qf") && haskey(branch, "qt")
                branch["qloss"] = branch["qf"] + branch["qt"]
            end
        end
    end

    if haskey(solution, "xfmr")
        for (i, xfmr) in solution["xfmr"]
            if haskey(xfmr, "pt_fr") && haskey(xfmr, "pt_to")
                xfmr["ptloss"] = xfmr["pt_fr"] + xfmr["pt_to"]
            end
            if haskey(xfmr, "qt_fr") && haskey(xfmr, "qt_to")
                xfmr["qtloss"] = xfmr["qt_fr"] + xfmr["qt_to"]
            end
        end
    end

    if haskey(solution, "gen")
        for (i, gen) in solution["gen"]
            if haskey(gen, "crg") && haskey(gen, "cig")
                gen["cm"] = hypot(gen["crg"], gen["cig"])
            end
        end
    end

    if haskey(solution, "load")
        for (i, load) in solution["load"]
            if haskey(load, "crd") && haskey(load, "cid")
                load["cm"] = hypot(load["crd"], load["cid"])
            end
        end
    end
end

function append_indicators!(result, hdata)
    solu = result["solution"]
    fundamental = "1"
    harmonics = Set(n for (n,nw) in solu["nw"])
    nonfundamentalharmonics = setdiff(harmonics, [fundamental])

    for (i,bus) in solu["nw"][fundamental]["bus"]
        vfun = solu["nw"][fundamental]["bus"][i]["vr"] + im*solu["nw"][fundamental]["bus"][i]["vi"]
        v   = [solu["nw"][n]["bus"][i]["vr"] + im*solu["nw"][n]["bus"][i]["vi"]  for n in harmonics]
        vnonfun = [solu["nw"][n]["bus"][i]["vr"] + im*solu["nw"][n]["bus"][i]["vi"]  for n in nonfundamentalharmonics]

        rms = sqrt(sum(abs.(v).^2))
        thd = sqrt(sum(abs.(vnonfun).^2))/abs(vfun)

        delete!(solu["nw"][fundamental]["bus"][i], "w")
        delete!(solu["nw"][fundamental]["bus"][i], "vr")
        delete!(solu["nw"][fundamental]["bus"][i], "vi")

        solu["nw"][fundamental]["bus"][i]["vrms"] = rms
        solu["nw"][fundamental]["bus"][i]["vthd"] = thd
        solu["nw"][fundamental]["bus"][i]["vmaxrms"] = hdata["nw"][fundamental]["bus"][i]["vmaxrms"]
        solu["nw"][fundamental]["bus"][i]["vminrms"] = hdata["nw"][fundamental]["bus"][i]["vminrms"]
        solu["nw"][fundamental]["bus"][i]["thdmax"] = hdata["nw"][fundamental]["bus"][i]["thdmax"]

    end

    for (i,gen) in solu["nw"][fundamental]["gen"]
        gencost = hdata["nw"][fundamental]["gen"][i]["cost"]
        pg = gen["pg"]
        if length(gencost) >0
            gen["totcost"] = sum(gencost[c]*(pg)^c for c in 1:length(gencost))
        end

        cfun = solu["nw"][fundamental]["gen"][i]["crg"] + im*solu["nw"][fundamental]["gen"][i]["cig"]
        call   = [solu["nw"][n]["gen"][i]["crg"] + im*solu["nw"][n]["gen"][i]["cig"]  for n in harmonics]
        cnonfun = [solu["nw"][n]["gen"][i]["crg"] + im*solu["nw"][n]["gen"][i]["cig"]  for n in nonfundamentalharmonics]

        rms = sqrt(sum(abs.(call).^2))
        thd = sqrt(sum(abs.(cnonfun).^2))/abs(cfun)

        solu["nw"][fundamental]["gen"][i]["crms"] = rms
        solu["nw"][fundamental]["gen"][i]["cthd"] = thd
        solu["nw"][fundamental]["gen"][i]["c_rating"] = hdata["nw"][fundamental]["gen"][i]["c_rating"]
        solu["nw"][fundamental]["gen"][i]["pmax"] = hdata["nw"][fundamental]["gen"][i]["pmax"]
    end

    if haskey(solu["nw"][fundamental], "load")
        for (i,load) in solu["nw"][fundamental]["load"]

            cfun = solu["nw"][fundamental]["load"][i]["crd"] + im*solu["nw"][fundamental]["load"][i]["cid"]
            call   = [solu["nw"][n]["load"][i]["crd"] + im*solu["nw"][n]["load"][i]["cid"]  for n in harmonics]
            cnonfun = [solu["nw"][n]["load"][i]["crd"] + im*solu["nw"][n]["load"][i]["cid"]  for n in nonfundamentalharmonics]

            rms = sqrt(sum(abs.(call).^2))
            thd = sqrt(sum(abs.(cnonfun).^2))/abs(cfun)

            solu["nw"][fundamental]["load"][i]["crms"] = rms
            solu["nw"][fundamental]["load"][i]["cthd"] = thd
        end
    end
end


text:
export dualize

# MOI dualize
function dualize(primal_model::MOI.ModelLike; dual_names::DualNames = DualNames("", ""))
    # Creates an empty dual problem
    dual_problem = DualProblem{Float64}()
    return dualize(primal_model, dual_problem, dual_names)
end

function dualize(primal_model::MOI.ModelLike, dual_problem::DualProblem{T};
                 dual_names::DualNames = DualNames("", "")) where T
    # Dualize with the optimizer already attached
    return dualize(primal_model, dual_problem, dual_names)
end

function dualize(primal_model::MOI.ModelLike, dual_problem::DualProblem{T}, dual_names::DualNames) where T
    # Throws an error if objective function cannot be dualized
    supported_objective(primal_model)

    # Query all constraint types of the model
    con_types = MOI.get(primal_model, MOI.ListOfConstraints())
    supported_constraints(con_types) # Throws an error if constraint cannot be dualized

    # Set the dual model objective sense
    set_dual_model_sense(dual_problem.dual_model, primal_model)

    # Get Primal Objective Coefficients
    primal_objective = get_primal_objective(primal_model)

    # Add variables to the dual model and their dual cone constraint.
    # Return a dictionary for dual variables with primal constraints
    dual_obj_affine_terms = add_dual_vars_in_dual_cones(dual_problem.dual_model, primal_model,
                                                        dual_problem.primal_dual_map,
                                                        dual_names, con_types)

    # Fill Dual Objective Coefficients Struct
    dual_objective = get_dual_objective(dual_problem.dual_model,
                                        dual_obj_affine_terms, primal_objective)

    # Add dual objective to the model
    set_dual_objective(dual_problem.dual_model, dual_objective)

    # Add dual equality constraint and get the link dictionary
    add_dual_equality_constraints(dual_problem.dual_model, primal_model,
                                  dual_problem.primal_dual_map, dual_names,
                                  primal_objective, con_types)

    return dual_problem
end

# JuMP dualize
function dualize(model::JuMP.Model; dual_names::DualNames = DualNames("", ""))
    # Create an empty JuMP model
    JuMP_model = JuMP.Model()

    if JuMP.mode(model) != JuMP.AUTOMATIC # Only works in AUTOMATIC mode
        error("Dualization does not support solvers in $(model.moi_backend.mode) mode")
    end
    # Dualize and attach to the model
    dualize(backend(model), DualProblem(backend(JuMP_model)); dual_names = dual_names)

    return JuMP_model
end

function dualize(model::JuMP.Model, factory::OptimizerFactory; dual_names::DualNames = DualNames("", ""))
    # Dualize the JuMP model
    dual_JuMP_model = dualize(model; dual_names = dual_names)
    # Set the optimizer
    JuMP.set_optimizer(dual_JuMP_model, factory)
    return dual_JuMP_model
end

# dualize docs
"""
    dualize(args...; kwargs...)

The `dualize` function works in three different ways. The user can provide:

* A `MathOptInterface.ModelLike`

The function will return a `DualProblem` struct that has the dualized model
and `PrimalDualMap{Float64}` for users to identify the links between primal and dual model.

* A `MathOptInterface.ModelLike` and a `DualProblem{T}`

* A `JuMP.Model`

The function will return a JuMP model with the dual representation of the problem.

* A `JuMP.Model` and an optimizer factory

The function will return a JuMP model with the dual representation of the problem with
the `OptimizerFactory` attached. The `OptimizerFactory` is the solver and its key arguments
that users provide in JuMP models, i.e. `with_optimizer(GLPK.Optimizer)`.

On each of these methods, the user can provide the keyword argument `dual_names`.
`dual_names` must be a `DualNames` struct. It allows users to set more intuitive names
for the dual variables and dual constraints created.
"""
function dualize end


text:
import Base64: base64decode

# from https://github.com/JuliaLang/julia/pull/36425
function detectwsl()
    Sys.islinux() &&
    isfile("/proc/sys/kernel/osrelease") &&
    occursin(r"Microsoft|WSL"i, read("/proc/sys/kernel/osrelease", String))
end

"""
    maybe_convert_path_to_wsl(path)

Return the WSL path if the system is using the Windows Subsystem for Linux (WSL) and return `path` otherwise.
WSL mounts the windows drive to /mnt/ and provides a utility tool to convert windows
paths into WSL paths. This function will try to use this tool to automagically
convert paths pasted from windows (with the right click -> copy as path functionality)
into paths Pluto can understand.

Example:
$(raw"C:\Users\pankg\OneDrive\Desktop\pluto\bakery_pnl_ready2.jl")
→
"/mnt/c/Users/pankg/OneDrive/Desktop/pluto/bakery_pnl_ready2.jl"

but "/mnt/c/Users/pankg/OneDrive/Desktop/pluto/bakery_pnl_ready2.jl" stays the same

"""
function maybe_convert_path_to_wsl(path)
    try
        if detectwsl()
            # wslpath utility prints path to stderr if it fails to convert
            # (it fails for WSL-valid paths)
            return readchomp(pipeline(`wslpath -u $(path)`; stderr=devnull))
        end
    catch e
        return path
    end
    return path
end

const adjectives = [
	"groundbreaking"
	"revolutionary"
	"important"
	"novel"
	"fun"
	"interesting"
	"fascinating"
	"exciting"
	"surprising"
	"remarkable"
	"wonderful"
	"stunning"
	"mini"
	"small"
	"tiny"
	"cute"
	"friendly"
	"wild"
]

const nouns = [
	"discovery"
	"experiment"
	"story"
	"journal"
	"notebook"
	"revelation"
	"computation"
	"creation"
	"analysis"
	"invention"
	"blueprint"
	"report"
	"science"
	"magic"
	"program"
	"notes"
	"lecture"
	"theory"
	"proof"
	"conjecture"
]

"""
Generate a filename like `"Cute discovery"`. Does not end with `.jl`.
"""
function cutename()
    titlecase(rand(adjectives)) * " " * rand(nouns)
end

function new_notebooks_directory()
    try
        path = joinpath(first(DEPOT_PATH), "pluto_notebooks")
        if !isdir(path)
            mkdir(path)
        end
        path
    catch
        homedir()
    end
end

"""
Standard Pluto file extensions, including `.jl` and `.pluto.jl`. Pluto can open files with any extension, but the default extensions are used when searching for notebooks, or when trying to create a nice filename for something else, like the backup file.
"""
const pluto_file_extensions = [
    ".pluto.jl",
    ".Pluto.jl",
	".nb.jl",
    ".jl",
    ".plutojl",
    ".pluto",
	".nbjl",
	".pljl",
	".pluto.jl.txt", # MacOS can create these .txt files sometimes
	".jl.txt",
]

endswith_pluto_file_extension(s) = any(endswith(s, e) for e in pluto_file_extensions)

function embedded_notebookfile(html_contents::AbstractString)::String
	if !occursin("</html>", html_contents)
		throw(ArgumentError("Pass the contents of a Pluto-exported HTML file as argument."))
	end

	m = match(r"pluto_notebookfile.*\"data\:.*base64\,(.*)\"", html_contents)
	if m === nothing
		throw(ArgumentError("Notebook does not have an embedded notebook file."))
	else
		String(base64decode(m.captures[1]))
	end
end

"""
Does the path end with a pluto file extension (like `.jl` or `.pluto.jl`) and does the first line say `### A Pluto.jl notebook ###`?
"""
is_pluto_notebook(path::String) = endswith_pluto_file_extension(path) && readline(path) == "### A Pluto.jl notebook ###"

function without_pluto_file_extension(s)
    for e in pluto_file_extensions
        if endswith(s, e)
            return s[1:prevind(s, ncodeunits(s), ncodeunits(e))]
        end
    end
    s
end

"""
Return `base` * `suffix` if the file does not exist yet.

If it does, return `base * sep * string(n) * suffix`, where `n` is the smallest natural number such that the file is new. (no 0 is not a natural number you snake)
"""
function numbered_until_new(base::AbstractString; sep::AbstractString=" ", suffix::AbstractString=".jl", create_file::Bool=true, skip_original::Bool=false)
    chosen = base * suffix
    n = 1
    while (n == 1 && skip_original) || isfile(chosen)
        chosen = base * sep * string(n) * suffix
        n += 1
	end
	if create_file
		touch(chosen)
	end
    chosen
end

backup_filename(path) = numbered_until_new(without_pluto_file_extension(path); sep=" backup ", suffix=".jl", create_file=false, skip_original=true)

"Like `cp` except we create the file manually (to fix permission issues). (It's not plagiarism if you use this function to copy homework.)"
function readwrite(from::AbstractString, to::AbstractString)
    write(to, read(from, String))
end

function tryexpanduser(path)
	try
		expanduser(path)
	catch ex
		path
	end
end

tamepath = abspath ∘ tryexpanduser

"Block until reading the file two times in a row gave the same result."
function wait_until_file_unchanged(filename::String, timeout::Real, last_contents::String="")::Nothing
	new_contents = try
        read(filename, String)
    catch
        ""
    end

    @info "Waiting for file to stabilize..."# last_contents new_contents

	if last_contents == new_contents
		# yayyy
        return
	else
        sleep(timeout)
		wait_until_file_unchanged(filename, timeout, new_contents)
	end
end

max_stars_repo_path: src/notebook/path helpers.jl

max_issues_repo_path: src/notebook/path helpers.jl

max_forks_repo_path: src/notebook/path helpers.jl


text:
include("activate.jl")

using LinearAlgebra
import Random

using MultivariatePolynomials
using DynamicPolynomials

using FiniteDifferences
import JSON
using Statistics

const N = 15

DynamicPolynomials.@polyvar X[1:15]

const max_degree = 4
coefficient_magnitude = 10
noise_magnitude = 1

const var_monomials = MultivariatePolynomials.monomials(X, 0:max_degree)

Random.seed!(42)
all_coeffs = map(var_monomials) do m
    d = MultivariatePolynomials.degree(m)
    return coefficient_magnitude * rand()
end
random_vector = rand(length(all_coeffs))
cutoff = quantile(random_vector, 0.95)
all_coeffs[findall(<(cutoff), random_vector)]  .= 0.0

const true_poly = dot(all_coeffs, var_monomials)

function evaluate_poly(coefficients)
    poly = dot(coefficients, var_monomials)
    return function p(x)
        return MultivariatePolynomials.subs(poly, Pair(X, x)).a[1]
    end
end

const training_data = map(1:500) do _
    x = 0.1 * randn(N)
    y = MultivariatePolynomials.subs(true_poly, Pair(X, x)) + noise_magnitude * randn()
    return (x, y.a[1])
end

const extended_training_data = map(training_data) do (x, y)
    x_ext = getproperty.(MultivariatePolynomials.subs.(var_monomials, X => x), :α)
    return (x_ext, y)
end

const test_data = map(1:1000) do _
    x = 0.4 * randn(N)
    y = MultivariatePolynomials.subs(true_poly, Pair(X, x)) + noise_magnitude * randn()
    return (x, y.a[1])
end

const extended_test_data = map(test_data) do (x, y)
    x_ext = getproperty.(MultivariatePolynomials.subs.(var_monomials, X => x), :α)
    return (x_ext, y)
end

function f(coefficients)
    return 0.5 / length(extended_training_data) * sum(extended_training_data) do (x, y)
        return (dot(coefficients, x) - y)^2
    end
end

function f_test(coefficients)
    return 0.5 / length(extended_test_data) * sum(extended_test_data) do (x, y)
        return (dot(coefficients, x) - y)^2
    end
end

function coefficient_errors(coeffs)
    return 0.5 * sum(eachindex(all_coeffs)) do idx
        return (all_coeffs[idx] - coeffs[idx])^2
    end
end

function grad!(storage, coefficients)
    storage .= 0
    for (x, y) in extended_training_data
        p_i = dot(coefficients, x) - y
        @. storage += x * p_i
    end
    storage ./= length(training_data)
    return nothing
end

function build_callback(trajectory_arr)
    return function callback(state)
        return push!(
            trajectory_arr,
            (Tuple(state)[1:5]..., f_test(state.x), coefficient_errors(state.x)),
        )
    end
end

#Check the gradient using finite differences just in case
gradient = similar(all_coeffs)

#Disable for now.
FrankWolfe.check_gradients(grad!, f, gradient)

max_iter = 100_000
random_initialization_vector = rand(length(all_coeffs))

#lmo = FrankWolfe.LpNormLMO{1}(100 * maximum(all_coeffs))

lmo = FrankWolfe.LpNormLMO{1}(0.95 * norm(all_coeffs, 1))

# L estimate
num_pairs = 10000
L_estimate = -Inf
gradient_aux = similar(gradient)
for i in 1:num_pairs
    global L_estimate
    x = compute_extreme_point(lmo, randn(size(all_coeffs)))
    y = compute_extreme_point(lmo, randn(size(all_coeffs)))
    grad!(gradient, x)
    grad!(gradient_aux, y)
    new_L = norm(gradient - gradient_aux) / norm(x - y)
    if new_L > L_estimate
        L_estimate = new_L
    end
end

# L1 projection
# inspired by https://github.com/MPF-Optimization-Laboratory/ProjSplx.jl
function projnorm1(x, τ)
    n = length(x)
    if norm(x, 1) ≤ τ
        return x
    end
    u = abs.(x)
    # simplex projection
    bget = false
    s_indices = sortperm(u, rev=true)
    tsum = zero(τ)

    @inbounds for i in 1:n-1
        tsum += u[s_indices[i]]
        tmax = (tsum - τ) / i
        if tmax ≥ u[s_indices[i+1]]
            bget = true
            break
        end
    end
    if !bget
        tmax = (tsum + u[s_indices[n]] - τ) / n
    end

    @inbounds for i in 1:n
        u[i] = max(u[i] - tmax, 0)
        u[i] *= sign(x[i])
    end
    return u
end

# gradient descent

xgd = FrankWolfe.compute_extreme_point(lmo, random_initialization_vector)
training_gd = Float64[]
test_gd = Float64[]
coeff_error = Float64[]
time_start = time_ns()
gd_times = Float64[]
for iter in 1:max_iter
    global xgd
    grad!(gradient, xgd)
    xgd = projnorm1(xgd - gradient / L_estimate, lmo.right_hand_side)
    push!(training_gd, f(xgd))
    push!(test_gd, f_test(xgd))
    push!(coeff_error, coefficient_errors(xgd))
    push!(gd_times, (time_ns() - time_start) * 1e-9)
end

@info "Gradient descent training loss $(f(xgd))"
@info "Gradient descent test loss $(f_test(xgd))"
@info "Coefficient error $(coefficient_errors(xgd))"

x00 = FrankWolfe.compute_extreme_point(lmo, random_initialization_vector)
x0 = deepcopy(x00)

# lazy AFW
trajectory_lafw = []
callback = build_callback(trajectory_lafw)
@time x_lafw, v, primal, dual_gap, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=max_iter,
    line_search=FrankWolfe.Adaptive(),
    print_iter=max_iter ÷ 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy=true,
    gradient=gradient,
    callback=callback,
    L=L_estimate,
);

@info "Lazy AFW training loss $(f(x_lafw))"
@info "Test loss $(f_test(x_lafw))"
@info "Coefficient error $(coefficient_errors(x_lafw))"

trajectory_bcg = []
callback = build_callback(trajectory_bcg)

x0 = deepcopy(x00)
@time x_bcg, v, primal, dual_gap, _ = FrankWolfe.blended_conditional_gradient(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=max_iter,
    line_search=FrankWolfe.Adaptive(),
    print_iter=max_iter ÷ 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    weight_purge_threshold=1e-10,
    callback=callback,
    L=L_estimate,
)

@info "BCG training loss $(f(x_bcg))"
@info "Test loss $(f_test(x_bcg))"
@info "Coefficient error $(coefficient_errors(x_bcg))"

x0 = deepcopy(x00)

#  compute reference solution using lazy AFW
trajectory_lafw_ref = []
callback = build_callback(trajectory_lafw_ref)
@time _, _, primal_ref, _, _ = FrankWolfe.away_frank_wolfe(
    f,
    grad!,
    lmo,
    x0,
    max_iteration=2 * max_iter,
    line_search=FrankWolfe.Adaptive(),
    print_iter=max_iter ÷ 10,
    emphasis=FrankWolfe.memory,
    verbose=true,
    lazy=true,
    gradient=gradient,
    callback=callback,
    L=L_estimate,
);

open(joinpath(@__DIR__, "polynomial_result.json"), "w") do f
    data = JSON.json((
        trajectory_arr_lafw=trajectory_lafw,
        trajectory_arr_bcg=trajectory_bcg,
        function_values_gd=training_gd,
        function_values_test_gd=test_gd,
        coefficient_error_gd=coeff_error,
        gd_times=gd_times,
        ref_primal_value=primal_ref,
    ))
    return write(f, data)
end

#Count missing\extra terms
print("\n Number of extra terms in GD: ", sum((all_coeffs .== 0) .* (xgd .!= 0)))
print("\n Number of missing terms in GD: ", sum((all_coeffs .!= 0) .* (xgd .== 0)))

print("\n Number of extra terms in BCG: ", sum((all_coeffs .== 0) .* (x_bcg .!= 0)))
print("\n Number of missing terms in BCG: ", sum((all_coeffs .!= 0) .* (x_bcg .== 0)))

print("\n Number of missing terms in Lazy AFW: ", sum((all_coeffs .== 0) .* (x_lafw .!= 0)))
print("\n Number of extra terms in Lazy AFW: ", sum((all_coeffs .!= 0) .* (x_lafw .== 0)))


text:
# The goal is to estimate the fractal dimension D of a curve by box counting
# algorith as N = eps^{-D}, where N is the number of boxes intersected
# by the curve at this level of refinement eps.

"""A hyper prism"""
struct BoundingBox{D, T}
    low::SVector{D, T}
    high::SVector{D, T}
end

"""Physical size of the bounding box; D vector of sizes"""
size{D, T}(box::BoundingBox{D, T}) = box.high - box.low

"""Length of bounding box in the i-th axis"""
size{D, T}(box::BoundingBox{D, T}, i::Int) = size(box)[i]

"""One refinement of a 2d box results in 4 boxes"""
function split{T}(box::BoundingBox{2, T})
    X, Y = box.low, box.high
    M = 0.5(X + Y)
    [BoundingBox(X, M),
     BoundingBox(SVector{2, T}(X[1], M[2]), SVector{2, T}(M[1], Y[2])),
     BoundingBox(SVector{2, T}(M[1], X[2]), SVector{2, T}(Y[1], M[2])),
     BoundingBox(M, Y)]
end

"""Bbox of segment"""
BoundingBox{D, T}(line::Segment{D, T}) = BoundingBox(SVector(Tuple(min(line.A[i], line.B[i]) for i in 1:D)),
                                                     SVector(Tuple(max(line.A[i], line.B[i]) for i in 1:D)))

"""Bbox of collection of segments"""
function BoundingBox{D, T}(segments::Vector{Segment{D, T}})

    low = SVector(Tuple(minimum(min(seg.A[i], seg.B[i]) for seg in segments) for i in 1:D))

    high = SVector(Tuple(maximum(max(seg.A[i], seg.B[i]) for seg in segments) for i in 1:D))

    return BoundingBox(low, high)
end

"""Bbox of a curve"""
function BoundingBox{D, T}(curve::Curve{D, T})
    segs = Vector{Segment{D, T}}(ne(curve))
    for (i, segment) in enumerate(iter_segments(curve))
        segs[i] = segment
    end
    BoundingBox(segs)
end

"""Does line collide with a box?"""
# NOTE: Liang-Barsky algorithm
function collides{D, T}(line::Segment{D, T}, box::BoundingBox{D, T})
    x, y = line.A
    dx = line.B[1] - x
    dy = line.B[2] - y

    p = [-dx, dx, -dy, dy]
    q = [x-box.low[1], box.high[1]-x, y-box.low[2], box.high[2]-y]

    u1 = -Inf
    u2 = Inf
    for i in 1:4
        (p[i] == 0 && q[i] < 0) && return false

        t = q[i]/p[i]

        if (p[i] < 0 && u1 < t)
            u1 = max(0, t)
        elseif (p[i] > 0 && u2 > t)
            u2 = min(1, t)
        end
    end

   (u1 > u2) && return false

    true
end

"""
Box counter holds after each refinement a pair of curve and box which
might intersect it. The curve is a collection if indices which refer
to segments of the unsplit curve.
"""
struct BoxCounter{D, T}
    segments::Vector{Segment{D, T}}
    seg_box::Deque{Pair{Vector{Int}, BoundingBox{D, T}}}
    stop_length::Float64
end

"""Box size at the current level"""
size{D, T}(bc::BoxCounter{D, T}) = size(first(bc.seg_box).second)

"""Construct BoxCounter for the curve"""
function BoxCounter{D, T}(curve::Curve{D, T})
    segs = Vector{Segment{D, T}}(ne(curve))
    for (i, segment) in enumerate(iter_segments(curve))
        segs[i] = segment
    end
    box = BoundingBox(segs)
    # Make the box slightly larger so that we don't hit corners
    X, Y = box.low, box.high
    center = (X+Y)/2.

    # Increaze the box a bit in an attempt/delay point isects
    # NOTE: as a result of this each run gives a slighly different
    # number - ideally avarage several runs
    shift1 = SVector(((Y - center).*(1. + rand(D)))...)
    shift2 = SVector(((Y - center).*(1. + rand(D)))...)

    box = BoundingBox(center - shift1, center + shift2)

    seg_box = Deque{Pair{Vector{Int}, BoundingBox{D, T}}}()
    push!(seg_box, Pair(collect(1:length(segs)), box))

    BoxCounter(segs, seg_box, min_segment(curve)/2)
end

# One box intersected. That box is bounding box
start{D, T}(bc::BoxCounter{D, T}) = (1, size(bc))

# Refine and count
function next{D, T}(bc::BoxCounter{D, T}, state)

    count = 0
    # Go over all the boxes on the current level
    for i in 1:length(bc.seg_box)

        index_segments, boxes = shift!(bc.seg_box)
        # Are children intersected
        for box in split(boxes)
            isected_segments = Vector{Int}()
            # Collect isects
            for index in index_segments
                collides(bc.segments[index], box) && push!(isected_segments, index)
            end
            # Increate count if there are and send for next round
            if !isempty(isected_segments)
                count += 1
                push!(bc.seg_box, Pair(isected_segments, box))
            end
        end
    end
    # We have halved the mesh
    (count, last(state)/2.)
end

# FIXME: It doesn't make sence to go much beyond the resolution of the curve
done{D, T}(c::BoxCounter{D, T}, state) = true #minimum(last(state)) < c.stop_length

"""Char size of the hyper box - best fit by cuve"""
volume{D, T}(p::SVector{D, T}) = (prod(p))^(1./D)

"""Remember volumes and counts as the refinement is done with box counter"""
mutable struct FractalDim{D, T}
    bc::BoxCounter{D, T}

    volumes::Vector{Float64}
    counts::Vector{Int}
    next_state::Tuple{Int, SVector{D, T}}
end

"""Estimate fd by lin fit of the stepped FractalDim"""
function estimate_fd{D, T}(fd::FractalDim{D, T}, slice=3)
    @assert !isempty(fd.volumes) && !isempty(fd.counts)

    if length(fd.counts) < slice
        x = log.(fd.volumes)
        y = log.(fd.counts)
    else
        x = log.(fd.volumes[end-slice:end])
        y = log.(fd.counts[end-slice:end])
    end
    -last(Polynomials.polyfit(x, y, 1).a)
end

"""Curve constructor"""
function FractalDim{D, T}(c::Curve{D, T})
    bc = BoxCounter(c)
    volumes = Vector{Float64}()
    counts = Vector{Int}()
    next_state = (1, size(bc))

    FractalDim(bc, volumes, counts, next_state)
end

"""Start stepping"""
function start{D, T}(fd::FractalDim{D, T})
    state = start(fd.bc)

    push!(fd.counts, first(state))
    push!(fd.volumes, volume(last(state)))
    fd.next_state = state

    state
end

"""Evolve"""
function next{D, T}(fd::FractalDim{D, T}, state)
    state = next(fd.bc, state)

    push!(fd.counts, first(state))
    push!(fd.volumes, volume(last(state)))
    fd.next_state = state

    state
end

"""Let box counter decide"""
done{D, T}(fd::FractalDim{D, T}, state) = done(fd.bc)

"""Evolve the started iterator nlevel times"""
function fractal_dim!{D, T}(fd::FractalDim{D, T}, nlevels::Int)
    state = fd.next_state
    for i in 1:nlevels
        state = next(fd, state)
    end
    (estimate_fd(fd), fd)
end

"""
Step the estimator for curve returning the counter which can be
further evolved if the precision is not okay
"""
fractal_dim{D, T}(c::Curve{D, T}, nlevels::Int) = fractal_dim!(FractalDim(c), nlevels)


text:
using BenchmarkTools
using Random

const SUITE = BenchmarkGroup()


text:
@testset "NDIndex" begin

x = NDIndex((1,2,3))
y = NDIndex((1,static(2),3))
z = NDIndex(static(3), static(3), static(3))

@testset "constructors" begin
    @test static(CartesianIndex(3, 3, 3)) === z == Base.setindex(Base.setindex(x, 3, 1), 3, 2)
    @test @inferred(ArrayInterface.Static.dynamic(z)) === CartesianIndex(3, 3, 3)
    @test @inferred(ArrayInterface.Static.known(z)) === (3, 3, 3)
    @test Tuple(@inferred(NDIndex{0}())) === ()
    @test @inferred(NDIndex{3}(1, static(2), 3)) === y
    @test @inferred(NDIndex{3}((1, static(2), 3))) === y
    @test @inferred(NDIndex{3}((1, static(2)))) === NDIndex(1, static(2), static(1))
    @test @inferred(NDIndex(x, y)) === NDIndex(1, 2, 3, 1, static(2), 3)
    @test @inferred(NDIndex{3,Tuple{Int,Int,Int}}((1,2, 3))) === x
end

@test @inferred(ArrayInterface.known_length(x)) === 3
@test @inferred(length(x)) === 3
@test @inferred(length(typeof(x))) === 3
@test @inferred(y[2]) === 2
@test @inferred(y[static(2)]) === static(2)

@test @inferred(-y) === NDIndex((-1,-static(2),-3))
@test @inferred(y + y) === NDIndex((2,static(4),6))
@test @inferred(y - y) === NDIndex((0,static(0),0))
@test @inferred(zero(x)) === NDIndex(static(0),static(0),static(0))
@test @inferred(oneunit(x)) === NDIndex(static(1),static(1),static(1))
@test @inferred(x * 3) === NDIndex((3,6,9))
@test @inferred(3 * x) === NDIndex((3,6,9))

@test @inferred(min(x, z)) === x
@test @inferred(max(x, z)) === NDIndex(3, 3, 3)
@test !@inferred(isless(y, x))
@test @inferred(isless(x, z))
@test @inferred(ArrayInterface.Static.lt(oneunit(z), z)) === static(true)

end


text:
"""
    AbstractEddyViscosityClosure <: AbstractTurbulenceClosure

Abstract supertype for turbulence closures that are defined by an isotropic viscosity
and isotropic diffusivities.
"""
abstract type AbstractEddyViscosityClosure{TD} <: AbstractIsotropicDiffusivity{TD} end

@inline viscosity(closure::AbstractEddyViscosityClosure, diffusivities) = diffusivities.νₑ
@inline diffusivity(::AbstractEddyViscosityClosure, diffusivities, ::Val{tracer_index}) where tracer_index = diffusivities.κₑ[tracer_index]


text:
#=
--- Day 6: Custom Customs ---

As your flight approaches the regional airport where you'll switch to a much larger plane, customs declaration forms are distributed to the passengers.

The form asks a series of 26 yes-or-no questions marked a through z. All you need to do is identify the questions for which anyone in your group answers "yes". Since your group is just you, this doesn't take very long.

However, the person sitting next to you seems to be experiencing a language barrier and asks if you can help. For each of the people in their group, you write down the questions for which they answer "yes", one per line. For example:

abcx
abcy
abcz

In this group, there are 6 questions to which anyone answered "yes": a, b, c, x, y, and z. (Duplicate answers to the same question don't count extra; each question counts at most once.)

Another group asks for your help, then another, and eventually you've collected answers from every group on the plane (your puzzle input). Each group's answers are separated by a blank line, and within each group, each person's answers are on a single line. For example:

abc

a
b
c

ab
ac

a
a
a
a

b

This list represents answers from five groups:

    The first group contains one person who answered "yes" to 3 questions: a, b, and c.
    The second group contains three people; combined, they answered "yes" to 3 questions: a, b, and c.
    The third group contains two people; combined, they answered "yes" to 3 questions: a, b, and c.
    The fourth group contains four people; combined, they answered "yes" to only 1 question, a.
    The last group contains one person who answered "yes" to only 1 question, b.

In this example, the sum of these counts is 3 + 3 + 3 + 1 + 1 = 11.

For each group, count the number of questions to which anyone answered "yes". What is the sum of those counts?
=#

function read_input()::Array{Array{String}}

    groups::Array{Array{String}} = []
    answers::Array{String} = []
    open("./input", "r") do f
        while ! eof(f)
            s::String = readline(f)
            if s == ""
                push!(groups,answers)
                answers = []
            else
                push!(answers,s)
            end
        end
    end

    push!(groups,answers)
    return groups
end

function count_each_group_yesses(groups::Array{Array{String}})::Int32

    count::Int32=0
    for group in groups
        answers::String = join(group)
        count += length(unique(answers))
    end
    return count
end

groups = read_input()
yesses_count = count_each_group_yesses(groups)
println("Total sum of different yesses for each group: ", yesses_count)

#=
--- Part Two ---

As you finish the last group's customs declaration, you notice that you misread one word in the instructions:

You don't need to identify the questions to which anyone answered "yes"; you need to identify the questions to which everyone answered "yes"!

Using the same example as above:

abc

a
b
c

ab
ac

a
a
a
a

b

This list represents answers from five groups:

    In the first group, everyone (all 1 person) answered "yes" to 3 questions: a, b, and c.
    In the second group, there is no question to which everyone answered "yes".
    In the third group, everyone answered yes to only 1 question, a. Since some people did not answer "yes" to b or c, they don't count.
    In the fourth group, everyone answered yes to only 1 question, a.
    In the fifth group, everyone (all 1 person) answered "yes" to 1 question, b.

In this example, the sum of these counts is 3 + 0 + 1 + 1 + 1 = 6.

For each group, count the number of questions to which everyone answered "yes". What is the sum of those counts?
=#

function count_each_group_yesses_intersect(groups::Array{Array{String}})::Int32

    count::Int32=0
    for group in groups
        count += length(intersect(group...))
    end
    return count
end

intersect_count = count_each_group_yesses_intersect(groups)
println("Total sum of common yesses for each group: ", intersect_count)


text:
seed!(1337)
x = randn(10000)

@testset "default params" begin
    p = @inferred histogram(x)
    @test_reference(
        "references/histogram/default.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    @test_reference(
        "references/histogram/default_nocolor.txt",
        @io2str(print(IOContext(::IO, :color=>false), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(x, closed = :left)
    @test_reference(
        "references/histogram/default.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    hist = fit(Histogram, x, closed = :left)
    p = @inferred histogram(hist)
    @test_reference(
        "references/histogram/default.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(x*100)
    @test_reference(
        "references/histogram/default_1e2.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(x*0.01)
    @test_reference(
        "references/histogram/default_1e-2.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(x, xscale = log10)
    @test_reference(
        "references/histogram/log10.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(x, xlabel = "custom label", xscale = log10)
    @test_reference(
        "references/histogram/log10_label.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram([0.1f0, 0.1f0, 0f0])
    @test_reference(
        "references/histogram/float32.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
end

@testset "hist params" begin
    hist = fit(Histogram, x, nbins = 5, closed = :right)
    p = @inferred histogram(hist)
    @test_reference(
        "references/histogram/hist_params.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(x, nbins = 5, closed = :right)
    @test_reference(
        "references/histogram/hist_params.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = nothing
    p = @test_logs (:warn, r"`bins`.+deprecated") @inferred histogram(x, bins = 5, closed = :right)
    @test_reference(
        "references/histogram/hist_params.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = nothing
    p = @test_logs (:warn, r"deprecated") @inferred histogram(x, 5, closed = :right)
    @test_reference(
        "references/histogram/hist_params.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
end

@testset "keyword arguments" begin
    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :blue,
        margin = 7,
        padding = 3,
    )
    @test_reference(
        "references/histogram/parameters1.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :blue,
        margin = 7,
        padding = 3,
        labels = false,
    )
    @test_reference(
        "references/histogram/parameters1_nolabels.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )

    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :yellow,
        border = :solid,
        symb = "=",
        width = 50
    )
    @test_reference(
        "references/histogram/parameters2.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
    # same but with Char as symb
    p = @inferred histogram(
        x,
        title = "My Histogram",
        xlabel = "Absolute Frequency",
        color = :yellow,
        border = :solid,
        symb = '=',
        width = 50
    )
    @test_reference(
        "references/histogram/parameters2.txt",
        @io2str(print(IOContext(::IO, :color=>true), p)),
        render = BeforeAfterFull()
    )
end


text:
@doc raw"""
```
l = segment_image(Algorithm::GrowCut, feature_vector::AbstractArray, seeds::AbstractArray{Int, 2}, t1::Int=9, t2::Int=9; max_iter::Int=1000, converge_at::Int=1, rule::Tuple{Function, Function}=(colour_diff,find_maxC))
```

Using a provided `feature_vector` (default rule uses RGB image) and an array of
seeds, `GrowCut` will use Cellular Automata to segment an image into different
regions. The algorithm takes two thresholds which can be used to control the
smoothness of the segmentation.

# Output

Returns a Array of Integers `l` that specifies the label assigned to each pixel.

# Details

The algorithm uses Cellular Automata to iteratively apply a rule to each pixel
which then competes with its neighbours to label the pixel. By default the
algorithm takes an RGB image for the `feature_vector` however a custom rule can
be defined for other `feature_vector` types. To improve performance after the
first iteration only the neighbourhood surrounding a pixel that changed in the
previous iteration is checked.

# Arguments

The function arguments are described in more detail below.

##  `feature_vector`

An `AbstractArray` containing the feature vector for each pixel. The default rule
assumes this is an RGB value for the pixel. Note: RGB images should be converted
to RGB{Float64} to prevent overflow.

##  `seeds`

An `Array{Int, 2}` specifying initial labels. Any pixel that is a seed should be
set to the desired integer label and all other positions should be 0. This can
be generated using `set_seed_pixels`

## `t1` and `t2`

`Int`s that specifiy the thresholds to be used for smoothing. `t1` is the maximum
number of "enemies" that a pixel can be surrounded by and still attack other
pixels. `t2` is the maximum number of "enemies" before the pixel is forced to
become the weakest enemy pixel. Note: these values can produce scenarios where
the algorithm is unable to converge and will instead reach its maximum number of
iterations.

## `maxiter`

An `Int` that specifies the maximum number of iterations before execution is
aborted. If left unspecified a default value of 1000 is used. This should never
be reached under normal conditions.

## `converge_at`

An `Int` that specifies the minimum number of changes across the image for
convergence. If left unspecified a default value of 1 is used.

## `rule`

A `Tuple{Function, Function}` that specifies the rule to be used for the
function. The first element is the rule that determines the strength of the
attacking pixel and the second element is used to calculate any extra
information needed by the rule. By default `colour_diff` and `find_maxC` are used. A
template for custom rules can be found using `help_rule_template` or generated into
a new file using `create_rule`.

# Example

Segment the "toucan" image using `GrowCut`

```julia
using Images, PictureSegmentation
img = testimage("toucan")

#convert to RGB{Float64}
img=RGB{Float64}.(img)

#create array of seeds in initialized to 0
seeds = zeros(Int,axes(img))

#use set_seed_pixels to set the inital seeds
set_seed_pixels(clicks, img, suppress_feedback=false)

#segment the image
l=segment_image(GrowCut(),img,clicks,9,9)

#create and display segmented image from labels
img2 = convert_labels(l)
img2
```

# Advanced

By default GrowCut uses the rule:

```math
\begin{aligned}
\textbf{for } \forall p &  \in  P                               &\\
                        & l_p^{t+1} =  l_p^{t}                  &\\
                        & \theta_p^{t+1} =  \theta_p^{t}        &\\
                        & \textbf{for } \forall q  \in  N(P)    & \\
                        &    \qquad    \textbf{if } g(\| RGB_p - RGB_q \|_2) \cdot \theta_q^t > \theta_p^t \\
                        &       \qquad \qquad l_p^{t+1} =  l_q^{t}  & \\
                        & \qquad  \qquad  \theta_p^{t+1} = g(\| RGB_p - RGB_q \|_2) \cdot \theta_q^t  \\
                         & \qquad     \textbf{end }  & \\
                        & \textbf{end }                         &\\
\textbf{end }\\
\end{aligned}
```
Where q is a pixel in the neighbourhood of p, l is the label of each pixel, θ is
the current strength of each pixel and RGB is the RGB vector of the pixel. g is
defined as ``g(x)=1-\dfrac{x}{\text{max}\|RGB\|}``.

An example of a custom rule for Lab colour is defined below:

```math
\begin{aligned}
\textbf{for } \forall p &  \in  P                               &\\
                        & l_p^{t+1} =  l_p^{t}                  &\\
                        & \theta_p^{t+1} =  \theta_p^{t}        &\\
                        & \textbf{for } \forall q  \in  N(P)    & \\
                        &    \qquad    \textbf{if } g(\| Lab_p - Lab_q \|_2) \cdot \theta_q^t > \theta_p^t \\
                        &       \qquad \qquad l_p^{t+1} =  l_q^{t}  & \\
                        & \qquad  \qquad  \theta_p^{t+1} = g(\| Lab_p - Lab_q \|_2) \cdot \theta_q^t  \\
                         & \qquad     \textbf{end }  & \\
                        & \textbf{end }                         &\\
\textbf{end }\\
\end{aligned}
```

where ``g(x)=1-\dfrac{x}{\text{max}\|Lab\|}``.
This would be implemented as per follows:

```julia

#calculate maxLab as needed by rule
function calculate_maxLab(LAB_color_vector::AbstractArray)
    maxLab = 0
    for i in CartesianIndices(feature_vector)
        current_Lab = sqrt(LAB_color_vector[i].l^2 + LAB_color_vector[i].a^2 + LAB_color_vector[i].b^2)
        if current_Lab > maxLab
            maxC = current_Lab
        end
    end
    return maxLab
end

#calculates the adjusted strength as per rule
function example_LAB_rule(θ::Array{Float64,1}, LAB_color_vector::AbstractArray, target_pixel::CartesianIndex, attacking_pixel::CartesianIndex, maxLab)

    LabTarget = LAB_color_vector[target_pixel]
    LabAttacker = LAB_color_vector[attacking_pixel]

    #calculate difference between LAB vectors
    LABdiff=Lab(LabTarget.l-LabAttacker.l,LabTarget.a-LabAttacker.a,LabTarget.b-LabAttacker.b)

    #find the absolute of the difference
    absLABdiff=sqrt(LABdiff.l^2+LABdiff.a^2+LABdiff.b^2)

    #calculate strength
    adjusted_strength = (1-(absLABdiff/maxLab))*θ[attacking_pixel]
    return adjusted_strength
end
```

The rule can then be passed to `Growcut` as per following:

```
l=segment_image(GrowCut(), img, clicks, rule=(example_LAB_rule,calculate_maxLab))
```

# Reference
1. V. Vezhnevets and V. Konouchine, "“GrowCut” - Interactive Multi-Label N-D Image Segmentation By Cellular Automata", in Graphicon, Novosibirsk Akademgorodok, Russia, 2005.
"""
function segment_image(Algorithm::GrowCut, feature_vector::AbstractArray, seeds::AbstractArray{Int, 2}, t1::Int=9, t2::Int=9; max_iter::Int=1000, converge_at::Int=1, rule::Tuple{Function, Function}=(colour_differance, find_maximum))
    # initialize values
    calculate_adjusted_strength = rule[1]
    extra_function = rule[2]
    extras = extra_function(feature_vector)
    labels = copy(seeds)
    E = zeros(Int, axes(labels))
    θ = set_strength(labels)
    labelsₜ₊₁ = copy(seeds)
    Eₜ₊₁ = zeros(Int, axes(labels))
    θₜ₊₁ = set_strength(labels)
    regions = maximum(seeds)
    iter = 0
    count = length(labels)
    changed = Array{Tuple{UnitRange{Int},UnitRange{Int}},1}(undef,1)
    changedₜ₊₁ = Array{Tuple{UnitRange{Int},UnitRange{Int}},1}(undef,0)
    rₘ, cₘ = size(labels)
    changed[1] = 1:rₘ, 1:cₘ

    # overflow warning
    if typeof(feature_vector) == Array{RGB,2} && typeof(feature_vector) != Array{RGB{Float64},2}
        @warn "feature_vector is RGB and not of type RGB{Float64}. This may cause overflow errors."
    end

    # iterate till convergence or maximum iterations reached
    while count > converge_at && iter < max_iter
        iter += 1
        count = 0
        copyto!(θₜ₊₁, θ)
        copyto!(labelsₜ₊₁, labels)
        copyto!(Eₜ₊₁, E)
        # loop through active pixel neighbourhoods (all pixels on first iteration)
        for j in CartesianIndices(changed)
            @inbounds for i in CartesianIndices(changed[j])
                r, c=i.I
                r₀ = max(r-1, 1)
                r₁ = min(r+1, rₘ)
                c₀ = max(c-1, 1)
                c₁ = min(c+1, cₘ)
                N₈ = r₀:r₁, c₀:c₁
                Eₜ₊₁[i] = count_enemy(labels,i,N₈)
                if E[i] < t2
                    θₜ₊₁, labelsₜ₊₁, count, changedₜ₊₁ = update_pixel(θ, θₜ₊₁, labels, labelsₜ₊₁, feature_vector, E, i, extras, t1, count, N₈, changedₜ₊₁, calculate_adjusted_strength)
                else
                    θₜ₊₁, labelsₜ₊₁, count, changedₜ₊₁ = occupy_pixel(θ, θₜ₊₁, labels, labelsₜ₊₁, feature_vector, i, extras,count, N₈, changedₜ₊₁, calculate_adjusted_strength)
                end
            end
        end

        # update values
        changed = changedₜ₊₁
        changedₜ₊₁ = empty(changedₜ₊₁)
        E = Eₜ₊₁
        θ = θₜ₊₁
        labels = labelsₜ₊₁
    end
    return labels
end