text:
using ImageFiltering
using ImageTransformations
using ImageCore

"""
    resizekern(source_size::Tuple{Int,Int}, dest_size::Tuple{Int,Int})

Create an image resize blur kernel, for use before reducing image size to avoid aliasing.
"""
function resizekern(source_size, dest_size)
    # Blur before resizing to prevent aliasing (kernel size dependent on both source and target image size)
    σ = 0.75 .* (source_size[1:2] ./ dest_size[1:2])
    if first(σ) < 1
        return ImageFiltering.KernelFactors.gaussian(σ)
    else
        return ImageFiltering.KernelFactors.IIRGaussian(σ)
    end
end

"""
    sizethatfits(original_size::Tuple{Int,Int}, target_shape::Tuple{Int,Int})

Takes an original image size, and fits it to a target shape for image resizing. Maintains aspect ratio.
"""
function sizethatfits(original_size, target_shape)
    channels = length(target_shape) == 2 ? 1 : target_shape[3]
    if original_size[1] > original_size[2]
        target_img_size = (
            target_shape[1],
            floor(Int, target_shape[2] * (original_size[2] / original_size[1])),
            channels
        )
    else
        target_img_size = (
        floor(Int, target_shape[1] * (original_size[1] / original_size[2])),
        target_shape[2],
        channels
        )
    end
    return target_img_size
end

"""
    resizePadImage(img::Array{T}, model::yolo) where {T<:ImageCore.Colorant}
    resizePadImage(img::Array{T}, target_img_size::Tuple{Int,Int}, kern) where {T<:.Color}
    resizePadImage(img::Array{T}, target_img::Array{U}, kern) where {T<:Color, U<:Float32}

Loads and prepares (resizes + pads) an image to fit within a given shape.
Returns the image and the padding.
"""
function resizePadImage(img::Array{T}, model::U) where {T<:ImageCore.Colorant, U<:Model}
    modelInputSize = getModelInputSize(model)
    if size(img) == modelInputSize[1:2]
        #TODO: Generalize for number of channels (i.e. some models trained on 1 channel)
        return permutedims(Float32.(channelview(img)), [3,2,1])[:,:,1:3]
    end
    img_size = size(img)
    target_img_size = sizethatfits(img_size, modelInputSize)
    kern = resizekern(img_size, target_img_size)
    return resizePadImage(img, modelInputSize, kern)
end
function resizePadImage(img::Array{T}, target_img_size::Tuple, kern) where {T<:ImageCore.Colorant}
    target_img = zeros(Float32, target_img_size[1:3])
    return resizePadImage(img, target_img, kern)
end
function resizePadImage(img::Array{T}, target_img::Array{U}, kern) where {T<:ImageCore.Colorant, U<:Float32}
    target_img_size = sizethatfits(size(img), size(target_img))
    # Determine top and left padding
    vpad_top = floor(Int, (size(target_img,1) - target_img_size[1]) / 2)
    hpad_left = floor(Int, (size(target_img,2) - target_img_size[2]) / 2)

    # Determine bottom and right padding accounting for rounding of top and left (to ensure accuate result image size if source has odd dimensions)
    vpad_bottom = size(target_img,1) - (vpad_top + target_img_size[1])
    hpad_right = size(target_img,2) - (hpad_left + target_img_size[2])

    padding = [hpad_left, vpad_top, hpad_right, vpad_bottom]

    vindex = (vpad_bottom+1):(size(target_img, 1) - vpad_top)
    hindex = (hpad_left+1):(size(target_img, 2) - hpad_right)

    if size(img,1) > size(target_img, 1) #Apply blur first if reducing size, to avoid aliasing
        imgr = ImageTransformations.imresize(ImageFiltering.imfilter(img, kern, NA()), target_img_size[1:2])
        target_img[vindex, hindex, :] .= permutedims(Float32.(channelview(imgr)), [3,2,1])[:,:,1:size(target_img,3)]
    elseif size(img,1) < size(target_img, 1)
        imgr = ImageTransformations.imresize(img, target_img_size[1:2])
        target_img[vindex, hindex, :] .= permutedims(Float32.(channelview(imgr)), [3,2,1])[:,:,1:size(target_img,3)]
    else
        target_img[vindex, hindex, :] .= permutedims(Float32.(channelview(img)), [3,2,1])[:,:,1:size(target_img,3)]
    end

    return target_img
end


text:
### Day 01 - Problem 1
## No Time for a Taxicab
## Author: Thanasis Georgiou

# Puzzle input
input = "R4, R3, R5, L3, L5, R2, L2, R5, L2, R5, R5, R5, R1, R3, L2, L2, L1, R5, L3, R1, L2, R1, L3, L5, L1, R3, L4, R2, R4, L3, L1, R4, L4, R3, L5, L3, R188, R4, L1, R48, L5, R4, R71, R3, L2, R188, L3, R2, L3, R3, L5, L1, R1, L2, L4, L2, R5, L3, R3, R3, R4, L3, L4, R5, L4, L4, R3, R4, L4, R1, L3, L1, L1, R4, R1, L4, R1, L1, L3, R2, L2, R2, L1, R5, R3, R4, L5, R2, R5, L5, R1, R2, L1, L3, R3, R1, R3, L4, R4, L4, L1, R1, L2, L2, L4, R1, L3, R4, L2, R3, L1, L5, R4, R5, R2, R5, R1, R5, R1, R3, L3, L2, L2, L5, R2, L2, R5, R5, L2, R3, L5, R5, L2, R4, R2, L1, R3, L5, R3, R2, R5, L1, R3, L2, R2, R1"

# Unpack complex number
function unpack(num)
    return real(num), imag(num)
end

# Current position
pos = 0 + 0im
dir = 0 + 1im

# History
history = zeros(Int8, 2000, 2000) # Large enough array (TM)

# HQ Location
hqLoc = 0 + 0im
hqFound = false

# Parse the input
steps = split(input, ", ")

# For each command
for step in steps
    # Find direction
    dir = begin
        command = step[1]
        if command == 'L'
            dir * 1im
        else
            dir * (-1im)
        end
    end

    # Move blocks
    blockDelta = parse(Int64, step[2:end])
    for i=1:1:blockDelta
        pos += dir

        # Store history
        x, y = unpack(pos)
        history[1000 - x, 1000 - y] += 1

        # If we have been here before, this is the HQ
        if hqFound == false && history[1000 - x, 1000 - y] > 1
            hqLoc = pos
            hqFound = true
        end
    end
end

# Calculate distance
x, y = unpack(pos)
#println("Position: $(real(pos)), $(imag(pos))")
println("Position: $x, $y")
println("Distance to position: $(abs(x) + abs(y))")

hqX, hqY = unpack(hqLoc)
println("HQ Location: $hqX, $hqY")
println("HQ Distance: $(abs(hqX) + abs(hqY))")

ext: jl


text:
include("euler/euler.jl")
using .Numbers: digits_sum
using .Primes: is_prime

function compute(n::Int)::Int
    harshad_numbers = collect(1:9)
    result = 0
    for _ ∈ 1:n - 2
        new_harshad_numbers = Vector{Int}()
        for i ∈ 0:9
            for harshad ∈ harshad_numbers
                if (harshad * 10 + i) % (digits_sum(harshad) + i) == 0
                    new_harshad = harshad * 10 + i
                    push!(new_harshad_numbers, new_harshad)
                    if is_prime(new_harshad ÷ digits_sum(new_harshad)) for j ∈ [1, 3, 7, 9] if is_prime(new_harshad * 10 + j) result += new_harshad * 10 + j end end end
                end
            end
        end
        harshad_numbers = new_harshad_numbers
    end
    return result
end

max_stars_repo_name: Dynortice/Project-Euler

max_issues_repo_name: Dynortice/Project-Euler

max_forks_repo_name: Dynortice/Project-Euler


text:
@info "Running functionality tests..."
@testset "Distributed Stochastic Programs" begin
    @testset "Distributed Sanity Check: $name" for (sp,res,name) in problems
        optimize!(sp)
        sp_nondist = copy(sp, procs = [1])
        add_scenarios!(sp_nondist, scenarios(sp))
        optimize!(sp_nondist,solver=GLPKSolverLP())
        @test scenariotype(sp) == scenariotype(sp_nondist)
        @test isapprox(stage_probability(sp), stage_probability(sp_nondist))
        @test nscenarios(sp) == nscenarios(sp_nondist)
        @test nscenarios(sp) == length(scenarios(sp))
        @test nsubproblems(sp) == nsubproblems(sp_nondist)
        @test isapprox(optimal_decision(sp), optimal_decision(sp_nondist))
        @test isapprox(optimal_value(sp), optimal_value(sp_nondist))
    end
    @testset "Distributed SP Constructs: $name" for (sp,res,name) in problems
        tol = 1e-2
        @test optimize!(sp) == :Optimal
        @test isapprox(optimal_decision(sp), res.x̄, rtol = tol)
        @test isapprox(optimal_value(sp), res.VRP, rtol = tol)
        @test isapprox(EWS(sp), res.EWS, rtol = tol)
        @test isapprox(EVPI(sp), res.EVPI, rtol = tol)
        @test isapprox(VSS(sp), res.VSS, rtol = tol)
        @test isapprox(EV(sp), res.EV, rtol = tol)
        @test isapprox(EEV(sp), res.EEV, rtol = tol)
    end
    @testset "Distributed Inequalities: $name" for (sp,res,name) in problems
        @test EWS(sp) <= VRP(sp)
        @test VRP(sp) <= EEV(sp)
        @test VSS(sp) >= 0
        @test EVPI(sp) >= 0
        @test VSS(sp) <= EEV(sp)-EV(sp)
        @test EVPI(sp) <= EEV(sp)-EV(sp)
    end
    @testset "Distributed Copying: $name" for (sp,res,name) in problems
        tol = 1e-2
        sp_copy = copy(sp)
        add_scenarios!(sp_copy, scenarios(sp))
        @test nscenarios(sp_copy) == nscenarios(sp)
        generate!(sp_copy)
        @test nsubproblems(sp_copy) == nsubproblems(sp)
        @test optimize!(sp_copy) == :Optimal
        optimize!(sp)
        @test isapprox(optimal_decision(sp_copy), optimal_decision(sp), rtol = tol)
        @test isapprox(optimal_value(sp_copy), optimal_value(sp), rtol = tol)
        @test isapprox(EWS(sp_copy), EWS(sp), rtol = tol)
        @test isapprox(EVPI(sp_copy), EVPI(sp), rtol = tol)
        @test isapprox(VSS(sp_copy), VSS(sp), rtol = tol)
        @test isapprox(EV(sp_copy), EV(sp), rtol = tol)
        @test isapprox(EEV(sp_copy), EEV(sp), rtol = tol)
    end
    @testset "Distributed Sampling" begin
        sampled_sp = sample(simple_model, sampler, 100, solver=GLPKSolverLP())
        @test nscenarios(sampled_sp) == 100
        @test nsubproblems(sampled_sp) == 100
        @test isapprox(stage_probability(sampled_sp), 1.0)
        sample!(sampled_sp, sampler, 100)
        @test nscenarios(sampled_sp) == 200
        @test nsubproblems(sampled_sp) == 200
        @test isapprox(stage_probability(sampled_sp), 1.0)
    end
    @testset "Distributed confidence intervals" begin
        glpk = GLPKSolverLP()
        try
            CI = confidence_interval(simple_model, sampler, solver = glpk, N = 200, log = false)
            @test lower(CI) <= upper(CI)
        catch end
        sol = optimize!(simple_model, sampler, solver = glpk, confidence = 0.95, tol = 1e-1, log = false)
        @test lower(confidence_interval(sol)) <= upper(confidence_interval(sol))
    end
end


text:
const ValBool = Union{Val{true}, Val{false}}

function foreach_thread(f, xs, rest...)
    tasks = Task[]
    i = 0
    for args in zip(xs, rest...)
        task = @task f(args...)
        task.sticky = false
        tid = i % Threads.nthreads()
        ccall(:jl_set_task_tid, Cvoid, (Any, Cint), task, tid)
        schedule(task)
        push!(tasks, task)
        i += 1
    end
    foreach(wait, tasks)
end


text:
# ----------------------------------- #
# Methods for saving figures to files #
# ----------------------------------- #

# TODO: add width and height and figure out how to convert from measures to the
#       pixels that will be expected in the SVG
function savefig_cairosvg(p::ElectronPlot, fn::AbstractString; dpi::Real=96)
    bas, ext = split(fn, ".")
    if !(ext in ["pdf", "png", "ps"])
        error("Only `pdf`, `png` and `ps` output supported")
    end
    # make sure plot window is active
    display(p)

    # write svg to tempfile
    temp = tempname()
    open(temp, "w") do f
        write(f, svg_data(p, ext))
    end

    # hand off to cairosvg for conversion
    run(`cairosvg $temp -d $dpi -o $fn`)

    # remove temp file
    rm(temp)

    # return plot
    p
end

# an alternative way to save plots -- no shelling out, but output less pretty

"""
`savefig(p::Plot, fn::AbstractString, js::Symbol)`

## Arguments

- `p::Plot`: Plotly Plot
- `fn::AbstractString`: Filename with extension (html, pdf, png)
- `js::Symbol`: One of the following:
    - `:local` - reference the javascript from PlotlyJS installation
    - `:remote` - reference the javascript from plotly CDN
    - `:embed` - embed the javascript in output (add's 1.7MB to size)
"""
function savefig_imagemagick(p::ElectronPlot, fn::AbstractString; js::Symbol=:local
                             #   sz::Tuple{Int,Int}=(8,6),
                             #   dpi::Int=300
                             )

    # Extract file type
    suf = split(fn, ".")[end]

    # if html we don't need a plot window
    if suf == "html"
        open(fn, "w") do f
            writemime(f, MIME"text/html"(), p, js)
        end
        return p
    end

    # for all the rest we need an active plot window
    show(p)

    # we can export svg directly
    if suf == "svg"
        open(fn, "w") do f
            write(f, svg_data(p))
        end
        return p
    end

    # now for the rest we need ImageMagick
    @eval import ImageMagick

    # construct a magic wand and read the image data from png
    wand = ImageMagick.MagickWand()
    # readimage(wand, _img_data(p, "svg"))
    ImageMagick.readimage(wand, base64decode(png_data(p)))
    ImageMagick.resetiterator(wand)

    # # set units to inches
    # status = ccall((:MagickSetImageUnits, ImageMagick.libwand), Cint,
    #       (Ptr{Void}, Cint), wand.ptr, 1)
    # status == 0 && error(wand)
    #
    # # calculate number of rows/cols
    # width, height = sz[1]*dpi, sz[2]*dpi
    #
    # # set resolution
    # status = ccall((:MagickSetImageResolution, ImageMagick.libwand), Cint,
    # (Ptr{Void}, Cdouble, Cdouble), wand.ptr, Cdouble(dpi), Cdouble(dpi))
    # status == 0 && error(wand)
    #
    # # set number of columns and rows
    # status = ccall((:MagickAdaptiveResizeImage, ImageMagick.libwand), Cint,
    #       (Ptr{Void}, Csize_t, Csize_t), wand.ptr, Csize_t(width), Csize_t(height))
    # status == 0 && error(wand)

    # finally write the image out
    ImageMagick.writeimage(wand, fn)

    p
end

function savefig(p::ElectronPlot, fn::AbstractString; js::Symbol=:local)
    suf = split(fn, ".")[end]

    # if html we don't need a plot window
    if suf == "html"
        open(fn, "w") do f
            writemime(f, MIME"text/html"(), p, js)
        end
        return p
    end

    # for all the rest we need raw svg data. to get that we'd have to display
    # the plot
    raw_svg = svg_data(p)

    # we can export svg directly
    if suf == "svg"
        open(fn, "w") do f
            write(f, raw_svg)
        end
        return p
    end

    # now we need to use librsvg/Cairo to finish
    try
        @eval import Rsvg
        @eval import Cairo
    catch e
        if isa(e, ArgumentError)
            msg = string("You need to install the Rsvg package use this",
                         " routine for file type $suf\n",
                         "Try insalling with `Pkg.add(\"Rsvg\")`")
            error(msg)
        else
            rethrow(e)
        end
    end

    if suf == "pdf"
        r = Rsvg.handle_new_from_data(raw_svg)
        cs = Cairo.CairoPDFSurface(fn, size(p.plot)...)
        ctx = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(ctx, r)
        Cairo.show_page(ctx)
        Cairo.finish(cs)
    elseif suf == "png"
        r = Rsvg.handle_new_from_data(raw_svg)
        cs = Cairo.CairoImageSurface(size(p.plot)...,Cairo.FORMAT_ARGB32)
        ctx = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(ctx, r)
        Cairo.write_to_png(cs, fn)
    elseif suf == "eps"
        r = Rsvg.handle_new_from_data(raw_svg)
        cs = Cairo.CairoEPSSurface(fn, size(p.plot)...)
        ctx = Cairo.CairoContext(cs)
        Rsvg.handle_render_cairo(ctx, r)
        Cairo.show_page(ctx)
        Cairo.finish(cs)
    else
        error("Only html, svg, png, pdf, eps output supported")
    end

    p
end

function png_data(p::ElectronPlot)
    raw = _img_data(p, "png")
    raw[length("data:image/png;base64,")+1:end]
end

function jpeg_data(p::ElectronPlot)
    raw = _img_data(p, "jpeg")
    raw[length("data:image/jpeg;base64,")+1:end]
end

function webp_data(p::ElectronPlot)
    raw = _img_data(p, "webp")
    raw[length("data:image/webp;base64,")+1:end]
end

# the above methods only work with an electron display. For now, we will
# implement them for all other displays simply by loading the plot into an
# electron display and calling the above methods.
for f in [:savefig_cairosvg, :savefig_imagemagick, :savefig,
          :png_data, :jpeg_data, :webp_data]
    @eval function $(f)(p::SyncPlot, args...; kwargs...)
        $(f)(ElectronPlot(p), args...; kwargs...)
        p
    end
end

const _mimeformats =  Dict("application/eps"         => "eps",
                           "image/eps"               => "eps",
                           "application/pdf"         => "pdf",
                           "image/png"               => "png",
                           "image/jpeg"              => "jpeg",
                           "application/postscript"  => "ps",
                           # "image/svg+xml"           => "svg"
)

# TODO: replace ElectronPlot with SyncPlot once I figure out how to get
#       img_data from within IJulia
for (mime, fmt) in _mimeformats
    @eval @compat function Base.show(io::IO, ::MIME{Symbol($mime)}, p::ElectronPlot)
        @eval import ImageMagick

        # construct a magic wand and read the image data from png
        wand = ImageMagick.MagickWand()
        ImageMagick.readimage(wand, base64decode(png_data(p)))
        ImageMagick.setimageformat(wand, $fmt)
        ImageMagick.writeimage(wand, io)

    end
end

for func in [:png_data, :jpeg_data, :wepb_data, :svg_data,
             :_img_data, :savefig, :savefig_cairosvg, :savefig_imagemagick]
    @eval function $(func)(::Plot, args...; kwargs...)
        msg = string("$($func) not available without a frontend. ",
                     "Try calling `$($func)(plot(p))` instead")
        error(msg)
    end
end


text:
using EarlyStopping, Dates, Test, InteractiveUtils

include("smoke.jl")

@testset "EarlyStopping" begin
    @testset "smoke" begin
        @testset "$C" for C in subtypes(StoppingCriterion)
            @eval @test_criteria $(C())
        end
    end

    @testset "criteria.jl" begin
        include("criteria.jl")
    end

    @testset "disjunction.jl" begin
        include("disjunction.jl")
    end

    @testset "object_oriented_api.jl" begin
        include("object_oriented_api.jl")
    end

    # to complete code coverage:
    @testset "api.jl" begin
        include("api.jl")
    end
end


text:
## FIXME: random OpSum{FermiOp} sometimes gets zeros and the size of the sum is wrong
@testset "OpSum" begin
    for _type = (Pauli, PauliI, FermiOp)
        m = rand(_type, 5, 10)
        @test size(OpSum(m)) == (5, 10)
    end
    for _type = (Pauli, PauliI, FermiOp)
        @test length(OpSum{_type}()) == 0
    end
    for _type in (Pauli, FermiOp)
        _sum = rand_op_sum(_type, 5, 5)
        @test 3 * _sum == _sum * 3
    end
end

@testset "plus, minus, zero, one" begin
    for op in (PauliSum(["XXZ", "YYX"]), FermiSum(["NNE", "++-"]))
        @test length(op - op) == 0
        @test iszero(op - op)
        @test zero(op) == (op - op)
        @test iszero(0 * op)
        @test isone(one(op))
        @test !isone(op)
        @test !iszero(op)
    end
    for term in (PauliTerm("XXX"), FermiTerm("---"))
        @test iszero(0 * term)
        @test iszero(term - term)
        @test iszero(zero(term))
        @test isone(one(term))
        @test !isone(term)
        @test !iszero(term)
        # TODO test that isone is not the expensive fallback method
    end
    ft = FermiTerm("---")
    @test iszero(ft * ft)
    pt = PauliTerm("XYZ")
    @test isone(pt * pt)
end

@testset "indexing" begin
    for _sum = (rand_op_sum(Pauli, 10, 10), rand_op_sum(FermiOp, 10, 10))
        @test _sum[2][3:6] == _sum[2, 3:6]
        @test _sum[2, 3:6] isa Vector
        @test length(_sum[2:5]) == 4
    end
end

@testset "OpTerm" begin
    for term in (PauliTerm("XYZ"), FermiTerm("--++"))
        @test reverse(term) != term
        @test reverse(reverse(term)) == term
        term2 = reverse(term)
        reverse!(term2)
        @test term2 == term
    end
    @test OpTerm(Paulis.X, Paulis.Y, Paulis.Z)  == OpTerm{Pauli}("XYZ")
    @test OpTerm(FermiOps.Raise, FermiOps.Lower, FermiOps.NumberOp) == OpTerm{FermiOp}("+-N")

    pt = PauliTerm("XXIIYY")
    @test weight(pt) == 4
    @test length(pt) == 6
    ft = FermiTerm("++II--")
    @test weight(ft) == 4
    @test length(ft) == 6

    @test OpTerm{Pauli}([Paulis.X, Paulis.Y]) == OpTerm([Paulis.X, Paulis.Y])
    (X, Y, Z) = (Paulis.X, Paulis.Y, Paulis.Z)
    @test PauliTerm(@SVector [X, Y, Z, Y ,Z]) ==  PauliTerm([X, Y, Z, Y ,Z])
end


text:
"""
    abstract type AbstractShape{S, T, N, L} end

The supertype for the various concrete shapes defined by `Shapes`. The,/StaticArrays.jl).
The `S` parameter is a `Tuple`-type specifying the dimensions, or size, of the
`AbstractShape`- such as `Tuple{3,4,5}` for a 3×4×5-sized array. The `T` parameter
specifies the underlying data type of the shape (e.g. the element type for an array
shape). The `L` parameter is the length of the array and is always equal to `prod(S)`.
Constructors may drop the `L` and `T` parameters if they are inferrable from the input
(e.g. `L` is always inferrable from `S`).
"""
abstract type AbstractShape{S,T,N,L} end

const AbstractScalarShape{T} = AbstractShape{Tuple{},T}
const AbstractVectorShape{S,T} = AbstractShape{Tuple{S},T}
const AbstractMatrixShape{S1,S2,T} = AbstractShape{Tuple{S1,S2},T}

Size(::Type{SH}) where {SH<:AbstractShape{S}} where {S<:Tuple} = Size(S)
Size(shape::SH) where {SH<:AbstractShape} = Size(SH)

Length(::Type{SH}) where {SH<:AbstractShape} = Length(Size(SH))
Length(shape::SH) where {SH<:AbstractShape} = Length(typeof(shape))

Base.eltype(::Type{SH}) where {S<:Tuple,T,SH<:AbstractShape{S,T}} = T

@pure Base.length(::Type{SH}) where {SH<:AbstractShape} = Int(Length(SH))
@inline Base.length(shape::AbstractShape) = length(typeof(shape))

@pure Base.size(::Type{SH}) where {SH<:AbstractShape} = Tuple(Size(SH))
@inline function Base.size(T::Type{<:AbstractShape}, d::Int)
    S = size(T)
    d > length(S) ? 1 : S[d]
end
@inline Base.size(shape::AbstractShape) = size(typeof(shape))
@inline Base.size(shape::AbstractShape, d::Int) = size(typeof(shape), d)

@pure Base.ndims(::Type{SH}) where {S<:Tuple,T,N,SH<:AbstractShape{S,T,N}} = N
Base.ndims(shape::AbstractShape) = ndims(typeof(shape))

Base.axes(shape::AbstractShape) = _axes(Size(shape))
Base.axes(::Type{SH}) where {SH<:AbstractShape} = _axes(Size(SH))
@pure _axes(::Size{sizes}) where {sizes} = map(SOneTo, sizes)

@inline concrete_eltype(::Type{SH}) where {SH<:AbstractShape} = default_datatype(eltype(SH))
@inline concrete_eltype(shape::AbstractShape) = concrete_eltype(typeof(shape))

"""
    checkaxes(shape::AbstractShape, x)

Throws an error if `axes(shape) != axes(x)``.
"""
@inline function checkaxes(shape::AbstractShape, x)
    if !checkaxes(Bool, shape, x)
        throw(DimensionMismatch("axes(x) must be equal to axes(shape)"))
    end
end

"""
    checkaxes(Bool, shape::AbstractShape, x)

Return `true` if `axes(shape) == axes(x)``.
"""
@inline function checkaxes(::Type{Bool}, shape::AbstractShape, x::AbstractArray)
    axes(x) == axes(shape)
end
@inline checkaxes(::Type{Bool}, ::AbstractScalarShape, ::Number) = true

"""
    checksize(shape::AbstractShape, x)

Throws an error if `size(shape) != size(x)`.
"""
@inline function checksize(shape::AbstractShape, x)
    if !checksize(Bool, shape, x)
        throw(DimensionMismatch("size(x) must be equal to size(shape)"))
    end
end

"""
    checksize(Bool, shape::AbstractShape, x)

Return `true` if `size(shape) == size(x)`.
"""
@inline function checksize(::Type{Bool}, shape::AbstractShape, x::AbstractArray)
    size(shape) == size(x)
end
@inline checksize(::Type{Bool}, ::AbstractScalarShape, ::Number) = true

struct Shape{S,T,N,L} <: AbstractShape{S,T,N,L}
    function Shape{S,T,N,L}() where {S<:Tuple,T,N,L}
        check_shape_params(S, T, Val{N}, Val{L})
        new{S,T,N,L}()
    end
end

@generated function (::Type{Shape{S,T,N}})() where {S<:Tuple,T,N}
    return quote
        @_inline_meta
        Shape{S,T,N,$(tuple_prod(S))}()
    end
end

@generated function (::Type{Shape{S,T}})() where {S<:Tuple,T}
    return quote
        @_inline_meta
        Shape{S,T,$(tuple_length(S)),$(tuple_prod(S))}()
    end
end

Shape(T::Type) = ScalarShape{T}()
Shape(T::Type, dim::Integer) = VectorShape{convert(Int, dim), T}()
Shape(T::Type, d1::Integer, d2::Integer) = MatrixShape{convert(Int, d1), convert(Int, d2), T}()
function Shape(T::Type, d1::Integer, d2::Integer, d3::Integer)
    d1 = convert(Int, d1)
    d2 = convert(Int, d2)
    d3 = convert(Int, d3)
    Shape{Tuple{d1, d2, d3}, T}()
end
function Shape(T::Type, dims::Vararg{Int, N}) where N
    dims = convert(NTuple{N, Int}, dims)
    Shape{tuple_to_size(dims), T, N}()
end

const ScalarShape{T} = Shape{Tuple{},T}
const VectorShape{S,T} = Shape{Tuple{S},T}
const MatrixShape{S1,S2,T} = Shape{Tuple{S1,S2},T}

ScalarShape(T::Type) = Shape(T)
VectorShape(T::Type, dim::Integer) = Shape(T, dim)
MatrixShape(T::Type, d1::Integer, d2::Integer) = Shape(T, d1, d2)

Shape(S::Size, T::Type) = Shape{tuple_to_size(Tuple(S)),T}()
Shape(::Type{A}) where {A<:AbstractArray} = Shape(Size(A), eltype(A))
function Shape(A::AbstractArray)
    size = Tuple(Size(A))
    check_dynamicsize(size)
    Shape{tuple_to_size(size),eltype(A)}()
end

check_dynamicsize(::Dims) = nothing
function check_dynamicsize(::Tuple)
    throw(ArgumentError(
        """Cannot create Shape from dynamically sized Array.
        Try Shape(s::Size, T) or Shape(T, dims...) instead."""
    ))
end

Adapt.adapt_storage(T::DataType, sh::Shape{S}) where {S <: Tuple} = Shape{S,T}()

# TODO unneccesary type params
struct MultiShape{S,T,N,L,namedtuple} <: AbstractShape{S,T,N,L}
    function MultiShape{S,T,N,L,namedtuple}() where {S,T,N,L,namedtuple}
        check_multishape_params(S, T, Val{N}, Val{L}, Val{namedtuple})
        new{S,T,N,L,namedtuple}()
    end
end

@generated function (::Type{MultiShape{namedtuple}})() where {namedtuple}
    check_multishape_namedtuple(Val{namedtuple})
    ntlength = namedtuple_length(namedtuple)
    ntsize = Tuple{ntlength}
    nteltype = promote_namedtuple_eltype(namedtuple)
    return quote
        @_inline_meta
        MultiShape{$ntsize,$nteltype,1,$ntlength,$namedtuple}()
    end
end

MultiShape(nt::NamedTuple) = MultiShape{nt}()
MultiShape(; shapes...) = MultiShape(values(shapes))
MultiShape(ms::MultiShape; shapes...) = MultiShape(merge(NamedTuple(ms), shapes))
MultiShape(shapes::Pair{Symbol, <:AbstractShape}...) = MultiShape(; shapes...)

@pure get(::Type{SH}) where {SH<:MultiShape} = SH.parameters[end]
get(ms::MultiShape) = get(typeof(ms))

NamedTuple(::SH) where {SH<:MultiShape} = get(SH)

Base.getproperty(ms::MultiShape, name::Symbol) = getfield(NamedTuple(ms), name)
Base.propertynames(ms::MultiShape) = propertynames(NamedTuple(ms))
Base.keys(ms::MultiShape) = keys(NamedTuple(ms))
Base.values(ms::MultiShape) = values(NamedTuple(ms))
Base.merge(ms1::MultiShape, ms2::MultiShape) = MultiShape(merge(NamedTuple(ms1), NamedTuple(ms2)))

@pure function namedtuple_length(nt::NamedTuple{<:Any,<:Tuple{Vararg{AbstractShape}}})
    sum(shape -> Int(Length(shape)), values(nt))
end

@generated function check_multishape_namedtuple(::Type{Val{namedtuple}}) where {namedtuple}
    if !(namedtuple isa NamedTuple)
        return :(throw(ArgumentError("MultiShape namedtuple parameter must be a namedtuple. Got: $namedtuple")))
    end

    if !all(x -> x isa AbstractShape, values(namedtuple))
        types = map(typeof, values(namedtuple))
        msg = "NamedTuple values must be AbstractShape. Got: $types"
        return :(throw(ArgumentError($msg)))
    end
end

@generated function check_multishape_params(
    ::Type{S},
    ::Type{T},
    ::Type{Val{N}},
    ::Type{Val{L}},
    ::Type{Val{namedtuple}},
) where {S,T,N,L,namedtuple}
    check_multishape_namedtuple(Val{namedtuple})
    check_shape_params(S, T, Val{N}, Val{L})

    ntlength = namedtuple_length(namedtuple)
    nteltype = promote_namedtuple_eltype(namedtuple)
    if !(nteltype <: T)
        msg = "Type mismatch in MultiShape parameters. Got type $T but promoted type of shapes was $nteltype"
        return :(throw(ArgumentError($msg)))
    end
    if L != ntlength
        msg = "Size mismatch in MultiShape parameters. Got length $L but sum of Shape lengths was $ntlength"
        return :(throw(ArgumentError($msg)))
    end
end

@pure Adapt.adapt_storage(T::DataType, ::MS) where {MS <: MultiShape} = MultiShape(map(s->Adapt.adapt(T, s), get(MS)))

@pure shapeindex(ms::MultiShape, name::Symbol) = Base.fieldindex(typeof(get(ms)), name)

@pure Base.firstindex(ms::MultiShape) = firstindex(get(ms))
@pure Base.lastindex(ms::MultiShape) = lastindex(get(ms))

@pure Base.getindex(ms::MultiShape, i::Int) = getindex(get(ms), i)

@pure function getoffset(ms::MultiShape, name::Symbol)
    idx = shapeindex(ms, name)
    idx == firstindex(ms) ? 0 : sum(i -> length(ms[i]), 1:(idx-1))
end


text:
abstract type Copula{d} <: Distributions.ContinuousMultivariateDistribution end
Base.length(::Copula{d}) where d = d

# The potential functions to code:
# Distributions._logpdf
# Distributions.cdf
# Distributions.fit(::Type{CT},u) where CT<:Mycopula
# Distributions._rand!
# Base.rand
# Base.eltype
# τ, τ⁻¹
# Base.eltype


text:
using SymPy
using Base.Test

## ODEs
x, a = Sym("x, a")
F = SymFunction("F")
ex = diff(F(x), x) - a*F(x)
ex1 = dsolve(ex)
ex2 = rhs(ex1) |> subs(Sym(:C1), 1) |> subs(a, 2)
@assert ex2 == exp(2x)

t, = @syms t
X, Y = map(SymFunction, ("X", "Y"))
eq = [Eq(diff(X(t),t), 12*t*X(t) + 8*Y(t)), Eq(diff(Y(t),t), 21*X(t) + 7*t*Y(t))]
dsolve(eq)

## version 0.4+ allow use of u'(x) in lieu of diff(u(x), x) and `ivpsolve`
u = SymFunction("u")
a,x, y, y0, y1 = symbols("a, x, y, y0, y1")

dsolve(u'(x) - a*u(x), x, (u, 0, 1))
dsolve(u'(x) - a*u(x), x, (u, 0, y1))
dsolve(u'(x) - a*u(x), x, (u, y0, y1))
dsolve(x*u'(x) + x*u(x) + 1, x, (u, 1, 1))
dsolve((u'(x))^2 - a*u(x), x, (u, 0, 1))
dsolve(u''(x) - a * u(x), x, (u, 0, 1), (u', 0, 0))

F, G, K = symbols("F, G, K", cls = symfunction)
eqn = F(x)*u'(y)*y + G(x)*u(y) + K(x)
dsolve(eqn, y, (u, 1, 0))

## dsolve eqn has two answers, but we want to eliminate based on initial condition
dsolve(u'(x) - (u(x)-1)*u(x)*(u(x)+1), x, (u, 0, 1//2))


text:
let
    function f_1(x::Vector)
        (x[1] - 5.0)^4
    end

    function g!_1(x::Vector, storage::Vector)
        storage[1] = 4.0 * (x[1] - 5.0)^3
    end

    function h!_1(x::Vector, storage::Matrix)
        storage[1, 1] = 12.0 * (x[1] - 5.0)^2
    end

    d = TwiceDifferentiableFunction(f_1, g!_1, h!_1)

    # Need to specify autodiff!
    @test_throws ErrorException Optim.optimize(DifferentiableFunction(f_1, g!_1), [0.0], Newton())
    Optim.optimize(DifferentiableFunction(f_1, g!_1), [0.0], Newton(), Optim.Options(autodiff = true))

    results = Optim.optimize(d, [0.0], Newton())
    @test_throws ErrorException Optim.x_trace(results)
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [5.0]) < 0.01

    eta = 0.9

    function f_2(x::Vector)
      (1.0 / 2.0) * (x[1]^2 + eta * x[2]^2)
    end

    function g!_2(x::Vector, storage::Vector)
      storage[1] = x[1]
      storage[2] = eta * x[2]
    end

    function h!_2(x::Vector, storage::Matrix)
      storage[1, 1] = 1.0
      storage[1, 2] = 0.0
      storage[2, 1] = 0.0
      storage[2, 2] = eta
    end

    d = TwiceDifferentiableFunction(f_2, g!_2, h!_2)
    results = Optim.optimize(d, [127.0, 921.0], Newton())
    @test_throws ErrorException Optim.x_trace(results)
    @assert Optim.g_converged(results)
    @assert norm(Optim.minimizer(results) - [0.0, 0.0]) < 0.01

    # Test Optim.newton for all twice differentiable functions in Optim.UnconstrainedProblems.examples
    for (name, prob) in Optim.UnconstrainedProblems.examples
    	if prob.istwicedifferentiable
    		ddf = TwiceDifferentiableFunction(prob.f, prob.g!,prob.h!)
    		res = Optim.optimize(ddf, prob.initial_x, Newton())
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
    	end
    end

    let
        prob=Optim.UnconstrainedProblems.examples["Himmelblau"]
        ddf = TwiceDifferentiableFunction(prob.f, prob.g!, prob.h!)
        res = optimize(ddf, [0., 0.], Newton())
        @assert norm(Optim.minimizer(res) - prob.solutions) < 1e-10
    end

    for (name, prob) in Optim.UnconstrainedProblems.examples
    	if prob.istwicedifferentiable
    		ddf = DifferentiableFunction(prob.f, prob.g!)
    		res = Optim.optimize(ddf, prob.initial_x, Newton(), Optim.Options(autodiff = true))
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
    		res = Optim.optimize(ddf.f, prob.initial_x, Newton(), Optim.Options(autodiff = true))
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
            res = Optim.optimize(ddf.f, ddf.g!, prob.initial_x, Newton(), Optim.Options(autodiff = true))
    		@assert norm(Optim.minimizer(res) - prob.solutions) < 1e-2
    	end
    end
end


text:
# gray-code implementation

export grayencode, graydecode, hamming_weight, nextgray

"gray-encode an integer."
function grayencode(n::Int) :: Int
    return xor(n, n >> 1)
end

"decode a gray-sequence integer"
function graydecode(n::Int) :: Int
    r = n
    while (n >>= 1) != 0
        r = xor(r, n)
    end
    return r
end

"return the number of ones in the base-2 representation of n."
function hamming_weight(n::Int) :: Int
    weight = 0
    i = one(n)
    while i > 0
        if i & n != 0
            weight += 1
        end
        i <<= 1
    end
    return weight
end

"return the first number after g in the gray sequence with hamming weight hw"
function nextgray(g::Int, hw::Int) :: Int
    wl = sizeof(g) * 8 # size in bits
    n = graydecode(g)
    for i in n+1:(1 << (wl-2))
        g = grayencode(i)
        if hamming_weight(g) == hw
            return g
        end
    end
    error("no number after $g of hamming weight $hw")
end

max_stars_repo_name: severinson/RaptorCodes

max_issues_repo_name: severinson/RaptorCodes

max_forks_repo_name: severinson/RaptorCodes


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

using MLJ
import DataFrames
import Statistics
using PrettyPrinting
using StableRNGs

rng = StableRNG(512)
Xraw = rand(rng, 300, 3)
y = exp.(Xraw[:,1] - Xraw[:,2] - 2Xraw[:,3] + 0.1*rand(rng, 300))
X = DataFrames.DataFrame(Xraw)

train, test = partition(eachindex(y), 0.7);

@load KNNRegressor
knn_model = KNNRegressor(K=10)

knn = machine(knn_model, X, y)

fit!(knn, rows=train)
ŷ = predict(knn, X[test, :]) # or use rows=test
rms(ŷ, y[test])

evaluate!(knn, resampling=Holdout(fraction_train=0.7, rng=StableRNG(666)),
          measure=rms) |> pprint

ensemble_model = EnsembleModel(atom=knn_model, n=20);

ensemble = machine(ensemble_model, X, y)
estimates = evaluate!(ensemble, resampling=CV())
estimates |> pprint

@show estimates.measurement[1]
@show mean(estimates.per_fold[1])

params(ensemble_model) |> pprint

B_range = range(ensemble_model, :bagging_fraction,
                lower=0.5, upper=1.0)
K_range = range(ensemble_model, :(atom.K),
                lower=1, upper=20);

tm = TunedModel(model=ensemble_model,
                tuning=Grid(resolution=10), # 10x10 grid
                resampling=Holdout(fraction_train=0.8, rng=StableRNG(42)),
                ranges=[B_range, K_range])

tuned_ensemble = machine(tm, X, y)
fit!(tuned_ensemble, rows=train);

best_ensemble = fitted_params(tuned_ensemble).best_model
@show best_ensemble.atom.K
@show best_ensemble.bagging_fraction

r = report(tuned_ensemble);

using PyPlot

figure(figsize=(8,6))

res = r.plotting
vals_b = res.parameter_values[:, 1]
vals_k = res.parameter_values[:, 2]

tricontourf(vals_b, vals_k, res.measurements)
xticks(0.5:0.1:1, fontsize=12)
xlabel("Bagging fraction", fontsize=14)
yticks([1, 5, 10, 15, 20], fontsize=12)
ylabel("Number of neighbors - K", fontsize=14)

ŷ = predict(tuned_ensemble, rows=test)
rms(ŷ, y[test])

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# ---
# title: 104. Maximum Depth of Binary Tree
# id: problem104
# author: zhwang
# date: 2022-02-20
# difficulty: Easy
# categories: Tree, Depth-first Search, Recursion
# link: <https://leetcode.com/problems/maximum-depth-of-binary-tree/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, return _its maximum depth_.
#
# A binary tree's **maximum depth**  is the number of nodes along the longest
# path from the root node down to the farthest leaf node.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/11/26/tmp-tree.jpg)
#
#
#
#     Input: root = [3,9,20,null,null,15,7]
#     Output: 3
#
#
# **Example 2:**
#
#
#
#     Input: root = [1,null,2]
#     Output: 2
#
#
# **Example 3:**
#
#
#
#     Input: root = []
#     Output: 0
#
#
# **Example 4:**
#
#
#
#     Input: root = [0]
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the tree is in the range `[0, 104]`.
#   * `-100 <= Node.val <= 100`
#
#
## @lc code=start
using LeetCode

maximum_depth_of_btree(::Nothing) = 0
maximum_depth_of_btree(root::TreeNode) = maximum(maximum_depth_of_btree, [root.left, root.right]) + 1
## @lc code=end

## @lc test=start
@testset "104.maximum-depth-of-binary-tree.jl" begin
    @test maximum_depth_of_btree(TreeNode{Int}([3,9,20,nothing,nothing,15,7])) == 3
    @test maximum_depth_of_btree(TreeNode{Int}([1,nothing,2])) == 2
    @test maximum_depth_of_btree(nothing) == 0
    @test maximum_depth_of_btree(TreeNode{Int}([5])) == 1
end
## @lc test=end

## @lc info=start
# link: [solution to 104](https://leetcode-cn.com/problems/maximum-depth-of-binary-tree/solution/104-er-cha-shu-de-zui-da-shen-du-bfs-dfs-2fbk/)
## @lc info=end


text:
function projection(d::Integer, q::Integer, elType=typeof(1.0))
    # Proj(deriv) = kron(diagm(0 => ones(elType, d)), [i==(deriv+1) ? 1 : 0 for i in 1:q+1]')

    # Slightly faster version of the above:
    D = d * (q + 1)
    Proj(deriv) = begin
        P = zeros(elType, d, D)
        @simd ivdep for i in deriv*d+1:D+1:d*D
            @inbounds P[i] = 1
        end
        return P
    end
    return Proj
end


text:
"""
    $(FUNCTIONNAME)(pol_type, obs_geom::ObsGeometry, τRayl, τAer, quad_points::QuadPoints, max_m, aerosol_optics, greek_rayleigh, τ_abs, brdf, architecture::AbstractArchitecture)

Perform Radiative Transfer calculations using given parameters

"""
function rt_run(pol_type::AbstractPolarizationType,   # Polarization type (IQUV)
                obs_geom::ObsGeometry,          # Solar Zenith, Viewing Zenith, Viewing Azimuthal
                τRayl,                          # Rayleigh optical depth
                τAer,                           # Aerosol optical depth and single-scattering albedo
                quad_points::QuadPoints,        # Quadrature points and weights
                max_m,                          # Max Fourier terms
                aerosol_optics,                 # AerosolOptics (greek_coefs, ω̃, k, fᵗ)
                greek_rayleigh::GreekCoefs,     # Greek coefficients of Rayleigh Phase Function
                τ_abs,                          # nSpec x Nz matrix of absorption
                brdf,                           # BRDF surface type
                architecture::AbstractArchitecture) # Whether to use CPU / GPU

    @unpack obs_alt, sza, vza, vaz = obs_geom   # Observational geometry properties
    @unpack qp_μ, wt_μ, qp_μN, wt_μN, iμ₀Nstart,μ₀, iμ₀,Nquad = quad_points # All quadrature points
    FT = eltype(sza)                    # Get the float-type to use
    Nz = length(τRayl)                  # Number of vertical slices
    nSpec = size(τ_abs, 1)              # Number of spectral points
    arr_type = array_type(architecture) # Type of array to use
    SFI = true                          # SFI flag
    NquadN = Nquad * pol_type.n         # Nquad (multiplied by Stokes n)
    dims = (NquadN,NquadN)              # nxn dims
    nAer  = length(aerosol_optics)      # Number of aerosols

    # Need to check this a bit better in the future!
    FT_dual = typeof(τAer[1])

    # Output variables: Reflected and transmitted solar irradiation at TOA and BOA respectively # Might need Dual later!!
    R = zeros(FT_dual, length(vza), pol_type.n, nSpec)
    T = zeros(FT_dual, length(vza), pol_type.n, nSpec)
    R_SFI = zeros(FT_dual, length(vza), pol_type.n, nSpec)
    T_SFI = zeros(FT_dual, length(vza), pol_type.n, nSpec)

    # Notify user of processing parameters
    @info "Processing on:" architecture
    @info "With FT:"  FT
    @info "Source Function Integration " SFI
    @info "Dimensions " NquadN,NquadN,nSpec
    println("Creating arrays...")

    # Create arrays
    @timeit "Creating layers" added_layer         = make_added_layer(FT_dual, arr_type, dims, nSpec)
    @timeit "Creating layers" added_layer_surface = make_added_layer(FT_dual, arr_type, dims, nSpec)
    @timeit "Creating layers" composite_layer     = make_composite_layer(FT_dual, arr_type, dims, nSpec)
    @timeit "Creating arrays" Aer𝐙⁺⁺ = arr_type(zeros(FT_dual, (dims[1], dims[2], nAer)))
    @timeit "Creating arrays" Aer𝐙⁻⁺ = similar(Aer𝐙⁺⁺)
    @timeit "Creating arrays" I_static = Diagonal(arr_type(Diagonal{FT}(ones(dims[1]))));

    println("Finished initializing arrays")

    # Loop over fourier moments
    for m = 0:max_m - 1

        println("Fourier Moment: ", m)

        # Azimuthal weighting
        weight = m == 0 ? FT(0.5) : FT(1.0)

        # Compute Z-moments of the Rayleigh phase matrix
        # For m>=3, Rayleigh matrices will be 0, can catch with if statement if wanted
        @timeit "Z moments" Rayl𝐙⁺⁺, Rayl𝐙⁻⁺ = Scattering.compute_Z_moments(pol_type, Array(qp_μ), greek_rayleigh, m, arr_type = arr_type);

        # Just for now (will change this later):
        # iBand = 1

        # Need to make sure arrays are 0:
        # TBD here

        # Compute aerosol Z-matrices for all aerosols
        for i = 1:nAer
            @timeit "Z moments"  Aer𝐙⁺⁺[:,:,i], Aer𝐙⁻⁺[:,:,i] = Scattering.compute_Z_moments(pol_type, Array(qp_μ), aerosol_optics[i].greek_coefs, m, arr_type = arr_type)
        end

        @timeit "Creating arrays" τ_sum_old = arr_type(zeros(FT, nSpec)) #Suniti: declaring τ_sum to be of length nSpec
        # @timeit "Creating arrays" τ_λ   = arr_type(zeros(FT, nSpec))

        # Loop over all layers and pre-compute all properties before performing core RT
        @timeit "Computing Layer Properties" computed_atmosphere_properties = construct_all_atm_layers(FT, nSpec, Nz, NquadN, τRayl, τAer, aerosol_optics, Rayl𝐙⁺⁺, Rayl𝐙⁻⁺, Aer𝐙⁺⁺, Aer𝐙⁻⁺, τ_abs, arr_type, qp_μ, μ₀, m)

        # Loop over vertical layers:
        @showprogress 1 "Looping over layers ..." for iz = 1:Nz  # Count from TOA to BOA

            # Construct the atmospheric layer
            # From Rayleigh and aerosol τ, ϖ, compute overall layer τ, ϖ
            computed_layer_properties = get_layer_properties(computed_atmosphere_properties, iz, arr_type)

            # Perform Core RT (doubling/elemental/interaction)
            rt_kernel!(pol_type, SFI, added_layer, composite_layer, computed_layer_properties, m, quad_points, I_static, architecture, qp_μN, iz)
        end

        # Create surface matrices:
        # arr_type(computed_atmosphere_properties.τ_sum_all[:,end])
        create_surface_layer!(brdf, added_layer, SFI, m, pol_type, quad_points, arr_type(computed_atmosphere_properties.τ_sum_all[:,end]), architecture);

        # One last interaction with surface:
        @timeit "interaction" interaction!(computed_atmosphere_properties.scattering_interfaces_all[end], SFI, composite_layer, added_layer, I_static)

        # Postprocess and weight according to vza
        postprocessing_vza!(iμ₀, pol_type, composite_layer, vza, qp_μ, m, vaz, μ₀, weight, nSpec, SFI, R, R_SFI, T, T_SFI)
    end

    # Show timing statistics
    print_timer()
    reset_timer!()

    #return R, T, R_SFI, T_SFI;
    if SFI
        return R_SFI;
    else
        return R;
    end

end

"""
    $(FUNCTIONNAME)(model::vSmartMOM_Model)

Perform Radiative Transfer calculations using parameters passed in through the
vSmartMOM_Model struct

"""
function rt_run(model::vSmartMOM_Model)

    return rt_run(model.params.polarization_type,
                  model.obs_geom::ObsGeometry,
                  model.τRayl,
                  model.τAer,
                  model.quad_points,
                  model.params.max_m,
                  model.aerosol_optics,
                  model.greek_rayleigh,
                  model.τ_abs,
                  model.brdf,
                  model.params.architecture)
end


text:
using FactCheck

facts("First facts group") do
  @fact 1 --> 1 "First group first test"
  @fact 2 --> 1 "First group second failing test"
end

facts("Second facts group") do
  @fact 1 --> 1 "Second group first test"
  @fact 2 --> 2 #Nameless fact
end

facts() do #Nameless group
  @fact 1 --> 1 #Nameless fact
  @fact 2 --> 1 #Failing nameless fact in nameless group
end
facts() do
end


text:
# advent of code 2020 day 6

function load_data(path:: String):: AbstractArray
    file = open(path)
    data_array = readlines(file)
    close(file)
    return data_array
end

function parse_group_one(answers_list)
    answers = ""
    for answer in answers_list
        answers *= answer
    end
    answers_array = [_c for _c in answers]
    answers_array = unique(answers_array)
    return length(answers_array)
end

function parse_group_two(answers_list)
    num_answers = 0
    for _c in 'a':'z'
        hits = 0
        for answers in answers_list
            hits += _c in answers
        end

        if length(answers_list) === hits
            num_answers += 1
        end
    end
    return num_answers
end

function process_array(data_array, foo)
    data_groups = []
    let answers = []
        for data in data_array
            if length(data) == 0
                push!(data_groups, foo(answers))
                answers = []
            else
                push!(answers, data)
            end
        end
        push!(data_groups, foo(answers))
    end
    return data_groups
end

data_array = load_data("data/day6.txt")
result_array = process_array(data_array, parse_group_one)
println("Result is $(sum(result_array)), with rule 1")
result_array = process_array(data_array, parse_group_two)
println("Result is $(sum(result_array)), with rule 2")


text:
using KitBase, Plots, Flux
using KitBase.JLD2

cd(@__DIR__)
include("../common.jl")

begin
    set = Setup(case = "layer", space = "1d1f3v", maxTime = 0.2, boundary = ["fix", "fix"], cfl = 0.5)
    ps = PSpace1D(-0.5, 0.5, 500, 1)
    vs = VSpace3D(-6.0, 6.0, 28, -6.0, 6.0, 64, -6.0, 6.0, 28)
    gas = Gas(Kn = 5e-3, K = 0.0)
    fw = function(x)
        prim = zeros(5)
        if x <= 0
            prim .= [1.0, 0.0, 1.0, 0.0, 1.0]
        else
            prim .= [1.0, 0.0, -1.0, 0.0, 2.0]
        end

        return prim_conserve(prim, ks.gas.γ)
    end
    ib = IB1F(fw, vs, gas)
    ks = SolverSet(set, ps, vs, gas, ib)
end

function get_sol(ks, ctrkt, ctrns, ctrad)
    solkt = zeros(ks.ps.nx, 5)
    solns = zero(solkt)
    solad = zero(solkt)
    for i in axes(solkt, 1)
        solkt[i, :] .= ctrkt[i].prim
        solns[i, :] .= ctrns[i].prim
        solad[i, :] .= ctrad[i].prim
        solkt[i, end] = 1 / solkt[i, end]
        solns[i, end] = 1 / solns[i, end]
        solad[i, end] = 1 / solad[i, end]
    end

    fc1 = (ctrkt[end÷2].f + ctrkt[end÷2+1].f) ./ 2
    hc1 = reduce_distribution(fc1, ks.vs.weights[:, 1, :], 2)
    fc3 = (ctrad[end÷2].f + ctrad[end÷2+1].f) ./ 2
    hc3 = reduce_distribution(fc3, ks.vs.weights[:, 1, :], 2)

    prim2 = (ctrns[end÷2].prim + ctrns[end÷2+1].prim) ./ 2
    sw2 = (ctrns[end÷2+1].w - ctrns[end÷2-1].w) ./ ks.ps.dx[end÷2] / 2
    τ = vhs_collision_time(prim2, ks.gas.μᵣ, ks.gas.ω)
    fc2 = chapman_enskog(ks.vs.u, ks.vs.v, ks.vs.w, prim2, sw2, zero(sw2), zero(sw2), ks.gas.K, τ)
    hc2 = reduce_distribution(fc2, ks.vs.weights[:, 1, :], 2)

    return solkt, solns, solad, hc1, hc2, hc3
end

###
# t = τ
###

begin
    @load "data/solkt_t.jld2" ctr
    ctrkt = deepcopy(ctr)
    @load "data/solns_t.jld2" ctr
    ctrns = deepcopy(ctr)
    @load "data/solkt_t.jld2" ctr
    ctrad = deepcopy(ctr)
end
solkt, solns, solad, hc1, hc2, hc3 = get_sol(ks, ctrkt, ctrns, ctrad)

begin
    begin
        plot(ks.ps.x[211:290], solkt[211:290, 1], lw=1.5, label="Kinetic", xlabel="x", ylabel="ρ")
        plot!(ks.ps.x[211:290], solns[211:290, 1], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[211:290], solad[211:290, 1], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_n_t.pdf")

    begin
        plot(ks.ps.x[211:290], solkt[211:290, 2], lw=1.5, label="Kinetic", xlabel="x", ylabel="U")
        plot!(ks.ps.x[211:290], solns[211:290, 2], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[211:290], solad[211:290, 2], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_u_t.pdf")

    begin
        plot(ks.ps.x[211:290], solkt[211:290, 3], lw=1.5, label="Kinetic", xlabel="x", ylabel="V")
        plot!(ks.ps.x[211:290], solns[211:290, 3], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[211:290], solad[211:290, 3], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_v_t.pdf")

    begin
        plot(ks.ps.x[211:290], solkt[211:290, 5], lw=1.5, label="Kinetic", xlabel="x", ylabel="T")
        plot!(ks.ps.x[211:290], solns[211:290, 5], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[211:290], solad[211:290, 5], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_t_t.pdf")
end

begin
    plot(ks.vs.v[1, :, 1], hc1, lw=1.5, label="Kinetic", xlabel="v", ylabel="f")
    plot!(ks.vs.v[1, :, 1], hc2, lw=1.5, line=:dash, label="NS")
    scatter!(ks.vs.v[1, :, 1], hc3, alpha=0.6, label="Adaptive")
end
savefig("figure/layer_f_t.pdf")

rg_ref = zeros(ks.ps.nx)
@inbounds Threads.@threads for i = 1:ks.ps.nx
    swx = (ctrkt[i+1].w - ctrkt[i-1].w) / (1e-6 + ks.ps.x[i+1] - ks.ps.x[i-1])
    rg_ref[i] = judge_regime(ks, ctrkt[i].f, ctrkt[i].prim, swx)
end

rg_kngll = zeros(ks.ps.nx)
@inbounds Threads.@threads for i = 1:ks.ps.nx
    swx = (ctrkt[i+1].w - ctrkt[i-1].w) / (1e-6 + ks.ps.x[i+1] - ks.ps.x[i-1])

    L = abs(ctrkt[i].w[1] / swx[1])
    ℓ = (1 / ctrkt[i].prim[end])^ks.gas.ω / ctrkt[i].prim[1] * sqrt(ctrkt[i].prim[end]) * ks.gas.Kn
    rg_kngll[i] = ifelse(ℓ / L > 0.05, 1, 0)
end

scatter(ks.ps.x[1:ks.ps.nx], rg_ref, alpha=0.7, label="NN", xlabel="x", ylabel="regime")
scatter!(ks.ps.x[1:ks.ps.nx], rg_kngll, alpha=0.7, label="KnGLL")
plot!(ks.ps.x[1:ks.ps.nx], rg_ref, lw=1.5, line=:dot, color=:gray27, label="True")
savefig("figure/layer_regime_t.pdf")

###
# t = 10τ
###

begin
    @load "data/solkt_10t.jld2" ctr
    ctrkt = deepcopy(ctr)
    @load "data/solns_10t.jld2" ctr
    ctrns = deepcopy(ctr)
    @load "data/solkt_10t.jld2" ctr
    ctrad = deepcopy(ctr)
end
solkt, solns, solad, hc1, hc2, hc3 = get_sol(ks, ctrkt, ctrns, ctrad)

begin
    begin
        plot(ks.ps.x[201:300], solkt[201:300, 1], lw=1.5, label="Kinetic", xlabel="x", ylabel="ρ")
        plot!(ks.ps.x[201:300], solns[201:300, 1], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[201:300], solad[201:300, 1], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_n_10t.pdf")

    begin
        plot(ks.ps.x[201:300], solkt[201:300, 2], lw=1.5, label="Kinetic", xlabel="x", ylabel="U")
        plot!(ks.ps.x[201:300], solns[201:300, 2], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[201:300], solad[201:300, 2], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_u_10t.pdf")

    begin
        plot(ks.ps.x[201:300], solkt[201:300, 3], lw=1.5, label="Kinetic", xlabel="x", ylabel="V")
        plot!(ks.ps.x[201:300], solns[201:300, 3], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[201:300], solad[201:300, 3], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_v_10t.pdf")

    begin
        plot(ks.ps.x[201:300], solkt[201:300, 5], lw=1.5, label="Kinetic", xlabel="x", ylabel="T")
        plot!(ks.ps.x[201:300], solns[201:300, 5], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[201:300], solad[201:300, 5], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_t_10t.pdf")
end

begin
    plot(ks.vs.v[1, :, 1], hc1, lw=1.5, label="Kinetic", xlabel="v", ylabel="f")
    plot!(ks.vs.v[1, :, 1], hc2, lw=1.5, line=:dash, label="NS")
    scatter!(ks.vs.v[1, :, 1], hc3, alpha=0.6, label="Adaptive")
end
savefig("figure/layer_f_10t.pdf")

rg_ref = zeros(ks.ps.nx)
@inbounds Threads.@threads for i = 1:ks.ps.nx
    swx = (ctrkt[i+1].w - ctrkt[i-1].w) / (1e-6 + ks.ps.x[i+1] - ks.ps.x[i-1])
    rg_ref[i] = judge_regime(ks, ctrkt[i].f, ctrkt[i].prim, swx)
end

rg_kngll = zeros(ks.ps.nx)
@inbounds Threads.@threads for i = 1:ks.ps.nx
    swx = (ctrkt[i+1].w - ctrkt[i-1].w) / (1e-6 + ks.ps.x[i+1] - ks.ps.x[i-1])

    L = abs(ctrkt[i].w[1] / swx[1])
    ℓ = (1 / ctrkt[i].prim[end])^ks.gas.ω / ctrkt[i].prim[1] * sqrt(ctrkt[i].prim[end]) * ks.gas.Kn
    rg_kngll[i] = ifelse(ℓ / L > 0.05, 1, 0)
end

begin
    scatter(ks.ps.x[1:ks.ps.nx], rg_ref, alpha=0.7, label="NN", xlabel="x", ylabel="regime")
    scatter!(ks.ps.x[1:ks.ps.nx], rg_kngll, alpha=0.7, label="KnGLL")
    plot!(ks.ps.x[1:ks.ps.nx], rg_ref, lw=1.5, line=:dot, color=:gray27, label="True")
end
savefig("figure/layer_regime_10t.pdf")

###
# t = 50τ
###

begin
    @load "data/solkt_50t.jld2" ctr
    ctrkt = deepcopy(ctr)
    @load "data/solns_50t.jld2" ctr
    ctrns = deepcopy(ctr)
    @load "data/solkt_50t.jld2" ctr
    ctrad = deepcopy(ctr)
end
solkt, solns, solad, hc1, hc2, hc3 = get_sol(ks, ctrkt, ctrns, ctrad)

begin
    begin
        plot(ks.ps.x[1:ks.ps.nx], solkt[:, 1], lw=1.5, label="Kinetic", xlabel="x", ylabel="ρ")
        plot!(ks.ps.x[1:ks.ps.nx], solns[:, 1], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[1:ks.ps.nx], solad[:, 1], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_n_50t.pdf")

    begin
        plot(ks.ps.x[1:ks.ps.nx], solkt[:, 2], lw=1.5, label="Kinetic", xlabel="x", ylabel="U")
        plot!(ks.ps.x[1:ks.ps.nx], solns[:, 2], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[1:ks.ps.nx], solad[:, 2], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_u_50t.pdf")

    begin
        plot(ks.ps.x[1:ks.ps.nx], solkt[:, 3], lw=1.5, label="Kinetic", xlabel="x", ylabel="V")
        plot!(ks.ps.x[1:ks.ps.nx], solns[:, 3], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[1:ks.ps.nx], solad[:, 3], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_v_50t.pdf")

    begin
        plot(ks.ps.x[1:ks.ps.nx], solkt[:, 5], lw=1.5, label="Kinetic", xlabel="x", ylabel="T")
        plot!(ks.ps.x[1:ks.ps.nx], solns[:, 5], lw=1.5, line=:dash, label="NS")
        scatter!(ks.ps.x[1:ks.ps.nx], solad[:, 5], alpha=0.6, label="Adaptive")
    end
    savefig("figure/layer_t_50t.pdf")
end

begin
    plot(ks.vs.v[1, :, 1], hc1, lw=1.5, label="Kinetic", xlabel="v", ylabel="f")
    plot!(ks.vs.v[1, :, 1], hc2, lw=1.5, line=:dash, label="NS")
    scatter!(ks.vs.v[1, :, 1], hc3, alpha=0.6, label="Adaptive")
end
savefig("figure/layer_f_50t.pdf")

rg_ref = zeros(ks.ps.nx)
@inbounds Threads.@threads for i = 1:ks.ps.nx
    swx = (ctrkt[i+1].w - ctrkt[i-1].w) / (1e-6 + ks.ps.x[i+1] - ks.ps.x[i-1])
    rg_ref[i] = judge_regime(ks, ctrkt[i].f, ctrkt[i].prim, swx)
end

rg_kngll = zeros(ks.ps.nx)
@inbounds Threads.@threads for i = 1:ks.ps.nx
    swx = (ctrkt[i+1].w - ctrkt[i-1].w) / (1e-6 + ks.ps.x[i+1] - ks.ps.x[i-1])

    L = abs(ctrkt[i].w[1] / swx[1])
    ℓ = (1 / ctrkt[i].prim[end])^ks.gas.ω / ctrkt[i].prim[1] * sqrt(ctrkt[i].prim[end]) * ks.gas.Kn
    rg_kngll[i] = ifelse(ℓ / L > 0.05, 1, 0)
end

begin
    scatter(ks.ps.x[1:ks.ps.nx], rg_ref, alpha=0.7, label="NN", xlabel="x", ylabel="regime")
    scatter!(ks.ps.x[1:ks.ps.nx], rg_kngll, alpha=0.7, label="KnGLL")
    plot!(ks.ps.x[1:ks.ps.nx], rg_ref, lw=1.5, line=:dot, color=:gray27, label="True")
end
savefig("figure/layer_regime_50t.pdf")

max_stars_repo_name: CSMMLab/Flowmachine

max_issues_repo_name: CSMMLab/Flowmachine

max_forks_repo_name: CSMMLab/Flowmachine


text:
using Test

using MathOptInterface
const MOI = MathOptInterface
const MOIT = MathOptInterface.Test
const MOIU = MathOptInterface.Utilities
const MOIB = MathOptInterface.Bridges

include("../utilities.jl")

mock = MOIU.MockOptimizer(MOIU.Model{Float64}())
config = MOIT.TestConfig()

bridged_mock = MOIB.Variable.Zeros{Float64}(mock)

x, cx = MOI.add_constrained_variable(bridged_mock, MOI.GreaterThan(0.0))
MOI.set(bridged_mock, MOI.VariableName(), x, "x")
MOI.set(bridged_mock, MOI.ConstraintName(), cx, "cx")
yz, cyz = MOI.add_constrained_variables(bridged_mock, MOI.Zeros(2))
MOI.set(bridged_mock, MOI.VariableName(), yz, ["y", "z"])
MOI.set(bridged_mock, MOI.ConstraintName(), cyz, "cyz")
y, z = yz
fx = MOI.SingleVariable(x)
fy = MOI.SingleVariable(y)
fz = MOI.SingleVariable(z)

@testset "SingleVariable objective" begin
    MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MIN_SENSE)
    MOI.set(bridged_mock, MOI.ObjectiveFunction{typeof(fx)}(), fx)
    @test MOI.get(bridged_mock, MOI.ObjectiveFunction{typeof(fx)}()) == fx
end

# Test before adding affine constraints are affine expressions cannot be
# unbridged when `Variable.ZerosBridge` is used.
@testset "Test bridged model" begin
    s = """
    variables: x, y, z
    cx: x >= 0.0
    cyz: [y, z] in MathOptInterface.Zeros(2)
    minobjective: x
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, s)
    MOIU.test_models_equal(bridged_mock, model, ["x", "y", "z"], ["cx", "cyz"])
end

c1, c2 = MOI.add_constraints(
    bridged_mock, [1.0fy + 1.0fz, 1.0fx + 1.0fy + 1.0fz],
    [MOI.EqualTo(0.0), MOI.GreaterThan(1.0)]
)
MOI.set(bridged_mock, MOI.ConstraintName(), c1, "con1")
MOI.set(bridged_mock, MOI.ConstraintName(), c2, "con2")
MOI.set(bridged_mock, MOI.ObjectiveSense(), MOI.MIN_SENSE)
obj = 1.0fx - 1.0fy - 1.0fz
MOI.set(bridged_mock, MOI.ObjectiveFunction{typeof(obj)}(), obj)

@test MOIB.Variable.unbridged_map(MOIB.bridge(bridged_mock, y), y, MOIB.Variable.IndexInVector(1)) === nothing
@test MOIB.Variable.unbridged_map(MOIB.bridge(bridged_mock, z), z, MOIB.Variable.IndexInVector(2)) === nothing

err = ErrorException(
    "Cannot delete constraint index of bridged constrained variables. Delete" *
    " the scalar variable or the vector of variables instead."
)
@test_throws err MOI.delete(bridged_mock, cyz)

err = ErrorException(
    "Cannot add two `VectorOfVariables`-in-`MathOptInterface.Zeros` on the" *
    " same first variable MathOptInterface.VariableIndex(-1)."
)
@test_throws err MOI.add_constraint(bridged_mock, MOI.VectorOfVariables(yz), MOI.Zeros(2))

err = ErrorException(
    "Cannot `VectorOfVariables`-in-`MathOptInterface.Zeros` for" *
    " which some variables are bridged but not the first one" *
    " `MathOptInterface.VariableIndex(12345679)`."
)
@test_throws err MOI.add_constraint(bridged_mock, MOI.VectorOfVariables([x, y]), MOI.Zeros(2))

err = ErrorException(
    "Cannot unbridge function because some variables are bridged by" *
    " variable bridges that do not support reverse mapping, e.g.," *
    " `ZerosBridge`."
)
@test_throws err MOI.get(bridged_mock, MOI.ObjectiveFunction{typeof(obj)}())
# With `c1`, the function does not contain any variable so it tests that it
# also throws an error even if it never calls `variable_unbridged_function`.
@test_throws err MOI.get(bridged_mock, MOI.ConstraintFunction(), c1)
@test_throws err MOI.get(bridged_mock, MOI.ConstraintFunction(), c2)

err = ArgumentError(
    "Variable bridge of type `MathOptInterface.Bridges.Variable.ZerosBridge{Float64}`" *
    " does not support accessing the attribute `MathOptInterface.Test.UnknownVariableAttribute()`."
)
@test_throws err MOI.get(bridged_mock, MOIT.UnknownVariableAttribute(), y)

@testset "Results" begin
    MOIU.set_mock_optimize!(mock,
        (mock::MOIU.MockOptimizer) -> MOIU.mock_optimize!(
            mock, [1.0],
            (MOI.ScalarAffineFunction{Float64}, MOI.EqualTo{Float64}) => 0.0,
            (MOI.ScalarAffineFunction{Float64}, MOI.GreaterThan{Float64}) => 1.0)
    )
    MOI.optimize!(bridged_mock)
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), x) == 1.0
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), y) == 0.0
    @test MOI.get(bridged_mock, MOI.VariablePrimal(), z) == 0.0

    @test MOI.get(bridged_mock, MOI.ConstraintPrimal(), cyz) == zeros(2)

    @test MOI.get(bridged_mock, MOI.ConstraintDual(), cx) == 0.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c1) == 0.0
    @test MOI.get(bridged_mock, MOI.ConstraintDual(), c2) == 1.0

    err = ArgumentError(
        "Bridge of type `MathOptInterface.Bridges.Variable.ZerosBridge{Float64}`" *
        " does not support accessing the attribute" *
        " `MathOptInterface.ConstraintDual(1)`."
    )
    @test_throws err MOI.get(bridged_mock, MOI.ConstraintDual(), cyz)
end

@testset "Query" begin
    @test MOI.get(bridged_mock, MOI.ConstraintFunction(), cyz).variables == yz
    @test MOI.get(mock, MOI.NumberOfVariables()) == 1
    @test MOI.get(mock, MOI.ListOfVariableIndices()) == [x]
    @test MOI.get(bridged_mock, MOI.NumberOfVariables()) == 3
    @test MOI.get(bridged_mock, MOI.ListOfVariableIndices()) == [x, y, z]
    @test MOI.get(mock, MOI.NumberOfConstraints{MOI.VectorOfVariables, MOI.Zeros}()) == 0
    @test MOI.get(bridged_mock, MOI.NumberOfConstraints{MOI.VectorOfVariables, MOI.Zeros}()) == 1
    @test MOI.get(bridged_mock, MOI.ListOfConstraintIndices{MOI.VectorOfVariables, MOI.Zeros}()) == [cyz]
end

@testset "Test mock model" begin
    s = """
    variables: x
    cx: x >= 0.0
    con1: x + 0.0 == 0.0
    con2: x + 0.0 >= 1.0
    minobjective: x
    """
    model = MOIU.Model{Float64}()
    MOIU.loadfromstring!(model, s)
    MOIU.test_models_equal(mock, model, ["x"], ["cx", "con1", "con2"])
end

@testset "Delete" begin
    test_delete_bridged_variables(bridged_mock, yz, MOI.Zeros, 3, (
        (MOI.SingleVariable, MOI.GreaterThan{Float64}, 1),
    ))
    @test MOI.is_valid(bridged_mock, x)
    @test !MOI.is_valid(bridged_mock, y)
    @test !MOI.is_valid(bridged_mock, z)
end


text:
"""
    continuous model
"""
abstract type Continuous <: Integration end

"""
    explicit midpoint integrator
"""
struct Midpoint <: Integration end

function fd(model::Model{Midpoint, FixedTime}, x, u, w, h, t)
    x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
end

function fd(model::Model{Midpoint, FreeTime}, x, u, w, h, t)
    h = u[end]
    x + h * f(model, x + 0.5 * h * f(model, x, u, w), u, w)
end

"""
    implicit midpoint integrator
"""
function fd(model::Model{Midpoint, FixedTime}, x⁺, x, u, w, h, t)
    x⁺ - (x + h * f(model, 0.5 * (x + x⁺), u, w))
end

function fd(model::Model{Midpoint, FreeTime}, x⁺, x, u, w, h, t)
    h = u[end]
    x⁺ - (x + h * f(model, 0.5 * (x + x⁺), u, w))
end

"""
    explicit RK3 integrator
"""
abstract type RK3 <: Integration end

function fd(model::Model{RK3, FixedTime}, z, u, w, h, t)
    # k1 = k2 = k3 = zero(z)
    k1 = h * f(model, z, u, w)
    k2 = h * f(model, z + 0.5 * k1, u, w)
    k3 = h * f(model, z - k1 + 2.0 * k2, u, w)
    z + (k1 + 4.0 * k2 + k3) / 6.0
end

abstract type RK4 <: Integration end

function fd(model::Model{RK4, FixedTime}, z, u, w, h, t)
	k1 = f(model, z, u, w) * h
	k2 = f(model, z + 0.5 * k1, u, w) * h
	k3 = f(model, z + 0.5 * k2, u, w) * h
	k4 = f(model, z + k3, u, w) * h
	z + (k1 + 2.0 * k2 + 2.0 * k3 + k4) / 6.0
end

"""
    discrete explicit dynamics
"""
abstract type Discrete <: Integration end


text:
@doc raw"""
Constructs constraint energy target data, and variable
# Constraints
`` varenergy[name, end]  >= paramenergytarget[name, end] ``
# LaTeX
`` x^{energy}_t  >= x^{energy}_{target} \text{ for } t = end ``
# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* time_series_data::Vector{DeviceTimeSeriesConstraintInfo} : Target reservoir storage forecast information
* cons_name::Symbol : energy target constraint name
* var_name::Symbol : the name of the Energy  variable
"""
function energy_target!(
    optimization_container::OptimizationContainer,
    target_data::Vector{T},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
) where {T <: DeviceTimeSeriesConstraintInfo}
    time_steps = model_time_steps(optimization_container)
    name_index = [get_component_name(d) for d in target_data]
    varenergy = get_variable(optimization_container, var_names[1])
    varslack_up = get_variable(optimization_container, var_names[2])
    varslack_dn = get_variable(optimization_container, var_names[3])

    target_constraint =
        add_cons_container!(optimization_container, cons_name, name_index, time_steps)

    for data in target_data, t in time_steps
        name = get_component_name(data)
        target_constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            varenergy[name, t] + varslack_up[name, t] + varslack_dn[name, t] ==
            data.multiplier * data.timeseries[t]
        )
    end

    return
end

@doc raw"""
Constructs constraint energy target data, and variable
# Constraints
`` varenergy[name, end]  >= paramenergytarget[name, end] ``
# LaTeX
`` x^{energy}_t  >= x^{energy}_{target} \text{ for } t = end ``
# Arguments
* optimization_container::OptimizationContainer : the optimization_container model built in PowerSimulations
* time_series_data::Vector{DeviceTimeSeriesConstraintInfo} : Target reservoir storage forecast information
* cons_names::Symbol : name of the constraint
* var_names::Symbol : the name of the energy variable
* param_reference::UpdateRef : UpdateRef to access the target parameter
"""
function energy_target_param!(
    optimization_container::OptimizationContainer,
    target_data::Vector{DeviceTimeSeriesConstraintInfo},
    cons_name::Symbol,
    var_names::Tuple{Symbol, Symbol, Symbol},
    param_reference::UpdateRef,
)
    time_steps = model_time_steps(optimization_container)
    name_index = [get_component_name(d) for d in target_data]
    varenergy = get_variable(optimization_container, var_names[1])
    varslack_up = get_variable(optimization_container, var_names[2])
    varslack_dn = get_variable(optimization_container, var_names[3])

    container_target = add_param_container!(
        optimization_container,
        param_reference,
        name_index,
        time_steps,
    )
    param_target = get_parameter_array(container_target)
    multiplier_target = get_multiplier_array(container_target)
    target_constraint =
        add_cons_container!(optimization_container, cons_name, name_index, time_steps)

    for d in target_data, t in time_steps
        name = get_component_name(d)
        param_target[name, t] =
            PJ.add_parameter(optimization_container.JuMPmodel, d.timeseries[t])
        multiplier_target[name, t] = d.multiplier
        target_constraint[name, t] = JuMP.@constraint(
            optimization_container.JuMPmodel,
            varenergy[name, t] + varslack_up[name, t] + varslack_dn[name, t] ==
            multiplier_target[name, t] * param_target[name, t]
        )
    end

    return
end


text:
using LinearAlgebra, SparseArrays

function simple_matrix(T::Type, nrow::Int, ncol::Int; symmetric::Bool = false)
  symmetric && @assert nrow == ncol
  if nrow == ncol == 1
    return ones(T, nrow, ncol)
  end
  U, _ = qr(rand(T, nrow, nrow))
  local V
  if symmetric
    V = U
  else
    V, _ = qr(rand(T, ncol, ncol))
  end
  S = T[(1 + (i - 1) / (nrow - 1)) * (i == j) for i = 1:nrow, j = 1:ncol]
  return U * S * V'
end

function simple_sparse_matrix(T::Type, nrow::Int, ncol::Int; symmetric::Bool = false)
  m = min(nrow, ncol)
  p = floor(Int, sqrt(m))
  q = div(m, p)
  M = kron(Matrix(I, q, q), simple_matrix(T, p, p, symmetric = symmetric))
  pq = p * q
  R = if m > pq
    simple_matrix(T, nrow - pq, ncol - pq, symmetric = symmetric)
  else
    spzeros(nrow - pq, ncol - pq)
  end
  A = [sparse(M) spzeros(pq, ncol - pq); spzeros(nrow - pq, pq) R]
end

simple_vector(T::Type, n::Int) = T[-(-one(T))^i for i = 1:n]


text:
abstract type AbstractLifeRule <: AbstractTDRule end

struct Life{T,A,C} <: AbstractLifeRule
    born::T
    survive::A
    radius::C
end

"""
Life(life_description; radius=1)

Returns a ```Life``` object given a tuple of tuples that follows the Golly notation ((b), (s)), where b stands for birth and s for survival.
These values indicates the number of neighbouring cells needed to birth a new one in the following generation, or to make the current alive
one survive.
"""
function Life(life_description; radius=1)
    born, survive = life_description[1], life_description[2]
    Life(born, survive, radius)
end

function (life::Life)(starting_array)
    return nextgen = life_evolution(starting_array, life.born, life.survive, life.radius)
end

function virtual_expansion(starting_array, radius)

    height, width = size(starting_array)
    nh, nw = height-radius+1, width-radius+1
    left = vcat(starting_array[nh:end,nw:end], starting_array[:,nw:end], starting_array[1:radius, nw:end])
    right = vcat(starting_array[nh:end, 1:radius], starting_array[:,1:radius], starting_array[1:radius,1:radius])
    middle = vcat(starting_array[nh:end,:], starting_array, starting_array[1:radius,:])

    hcat(left, middle, right)
end

function life_application(state, born, survive)

    past_value = state[size(state, 1)÷2+1, size(state, 2)÷2+1] #save past cell value
    state[size(state, 1)÷2+1, size(state, 2)÷2+1] = 0 #past cell value set to zero
    alive = sum(state) #the sum is the number of cell alive in the neighborhood since the central value is equal to zero

    if past_value == 1 && alive in survive
        return 1
    elseif past_value == 0 && alive in born
        return 1
    else
        return 0
    end
end

function life_evolution(starting_array, born, survive, radius)

    height, width = size(starting_array)
    output = zeros(typeof(starting_array[2]), height, width)
    virtual_output = virtual_expansion(starting_array, radius)

    for i = 1:height-radius+1, j = 1:width-radius+1
        output[i, j] = life_application(virtual_output[i:i+radius+1, j:j+radius+1], born, survive)
    end
    output
end


text:
abstract type HostMap{K,V} end
abstract type AbstractHashMap{K,V} <: HostMap{K,V} end
abstract type AbstractArrayMap{K,V} <: HostMap{K,V} end
abstract type AbstractBufferMap{K,V} <: HostMap{K,V} end

function hostmap(map::API.Map; K, V)
    fd = API.fd(map)
    @assert fd > 0 "Map \"$(API.name(map))\" is not yet allocated"
    def = API.definition(map)
    map_type = API.BPFMapType(def.type)
    jltype = maptype_to_jltype(Val(map_type))
    jltype{K,V}(fd)
end

struct HashMap{K,V} <: AbstractHashMap{K,V}
    fd::Cint
end
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_HASH}) = HashMap
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_STACK_TRACE}) = HashMap
struct ArrayMap{K,V} <: AbstractArrayMap{K,V}
    fd::Cint
end
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_ARRAY}) = ArrayMap

# perf buffers

const perf_buffer = Ptr{Cvoid}
struct PerfBuffer{K,V} <: AbstractBufferMap{K,V}
    fd::Cint
    buf::perf_buffer
end
PerfBuffer{K,V}(fd) where {K,V} = PerfBuffer{K,V}(fd, 8)
maptype_to_jltype(::Val{API.BPF_MAP_TYPE_PERF_EVENT_ARRAY}) = PerfBuffer
Base.@kwdef struct perf_buffer_opts
    sample_cb::Ptr{Cvoid} = C_NULL
    lost_cb::Ptr{Cvoid} = C_NULL
    ctx::Ptr{Cvoid} = C_NULL
end

function PerfBuffer{K,V}(map_fd::Cint, page_cnt::Integer, opts::perf_buffer_opts) where {K,V}
    opts_ref = ZeroInitRef(perf_buffer_opts, opts)
    buf = GC.@preserve opts_ref begin
        opts_ptr = Base.unsafe_convert(Ptr{perf_buffer_opts}, opts_ref)
        ccall((:perf_buffer__new, API.libbpf), perf_buffer, (Cint,Csize_t,Ptr{perf_buffer_opts}), map_fd, Csize_t(page_cnt), opts_ptr)
        #ccall((:perf_buffer__new, API.libbpf), perf_buffer, (Cint,Csize_t,Ptr{Cvoid},Ptr{Cvoid},Ptr{Cvoid},Ptr{perf_buffer_opts}), map_fd, Csize_t(page_cnt), opts.sample_cb, opts.lost_cb, opts.ctx, opts_ptr)
    end
    @assert Int(buf) > 0
    PerfBuffer{K,V}(map_fd, buf)
end
function PerfBuffer{K,V}(map_fd::Cint, page_cnt::Integer; sample_cb=C_NULL, lost_cb=C_NULL, ctx=C_NULL) where {K,V}
    opts = perf_buffer_opts(sample_cb, lost_cb, ctx)
    PerfBuffer{K,V}(map_fd, page_cnt, opts)
end
function poll(buf::PerfBuffer, timeout_ms::Integer)
    ccall((:perf_buffer__poll, API.libbpf), Cint, (perf_buffer,Cint), buf.buf, Cint(timeout_ms))
end

struct map_create_attr
    map_type::UInt32
    key_size::UInt32
    value_size::UInt32
    max_entries::UInt32
end

function HostMap(; map_type, key_type, value_type, max_entries)
    attr = Ref(map_create_attr(UInt32(map_type),
                               sizeof(key_type),
                               sizeof(value_type),
                               max_entries))
    fd = bpf(API.BPF_MAP_CREATE, attr)
    fd > 0 || Base.systemerror(fd)
    jltype = maptype_to_jltype(Val(map_type))
    jltype{key_type,value_type}(fd)
end

struct map_access_elem_attr
    map_fd::UInt32
    key::UInt64
    value_or_next::UInt64
    flags::UInt64
end

memset!(ptr::Ptr{T}) where T =
    ccall(:memset, Cvoid,
          (Ptr{T}, UInt8, UInt64),
          ptr, UInt8(0), sizeof(T))
"Creates a Ref{T} that's been zero-initialized before being stored. Necessary
to ensure that struct padding bytes are zeroed."
function ZeroInitRef(T, val; set=true)
    ref = Ref{T}()
    memset!(Base.unsafe_convert(Ptr{T}, ref))
    if set
        ref[] = val
    end
    ref
end
ZeroInitRef(val::T) where T = ZeroInitRef(T, val)
ZeroInitRef(::Type{T}) where T = ZeroInitRef(T, nothing; set=false)

function Base.getindex(map::AbstractHashMap{K,V}, idx) where {K,V}
    key = ZeroInitRef(K, idx)
    value = Ref{V}()
    key_ptr = Base.unsafe_convert(Ptr{K}, key)
    value_ptr = Base.unsafe_convert(Ptr{V}, value)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    value_ptr,
                                    0))
    ret = GC.@preserve key value begin
        bpf(API.BPF_MAP_LOOKUP_ELEM, attr)
    end
    ret >= 0 || Base.systemerror(ret)
    value[]
end
function Base.getindex(map::AbstractArrayMap{K,V}, idx) where {K,V}
    key = ZeroInitRef(K, idx-1)
    value = Ref{V}()
    key_ptr = Base.unsafe_convert(Ptr{K}, key)
    value_ptr = Base.unsafe_convert(Ptr{V}, value)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    value_ptr,
                                    0))
    ret = GC.@preserve key value begin
        bpf(API.BPF_MAP_LOOKUP_ELEM, attr)
    end
    ret >= 0 || Base.systemerror(ret)
    value[]
end

function Base.setindex!(map::AbstractHashMap{K,V}, value::U, idx) where {K,V,U}
    key_ref = ZeroInitRef(convert(K,idx))
    value_ref = ZeroInitRef(convert(V,value))
    key_ptr = Base.unsafe_convert(Ptr{K}, key_ref)
    value_ptr = Base.unsafe_convert(Ptr{V}, value_ref)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    value_ptr,
                                    0))
    ret = GC.@preserve key_ref value_ref begin
        bpf(API.BPF_MAP_UPDATE_ELEM, attr)
    end
    ret >= 0 || Base.systemerror(ret)
    value
end
function Base.setindex!(map::AbstractArrayMap{K,V}, value::U, idx) where {K,V,U}
    key_ref = ZeroInitRef(convert(K,idx-1))
    value_ref = ZeroInitRef(convert(V,value))
    key_ptr = Base.unsafe_convert(Ptr{K}, key_ref)
    value_ptr = Base.unsafe_convert(Ptr{V}, value_ref)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    value_ptr,
                                    0))
    ret = GC.@preserve key_ref value_ref begin
        bpf(API.BPF_MAP_UPDATE_ELEM, attr)
    end
    ret >= 0 || Base.systemerror(ret)
    value
end

function Base.delete!(map::AbstractHashMap{K,V}, idx) where {K,V}
    key = ZeroInitRef(K, idx)
    key_ptr = Base.unsafe_convert(Ptr{K}, key)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    0,
                                    0))
    GC.@preserve key begin
        @assert bpf(API.BPF_MAP_DELETE_ELEM, attr) >= 0
    end
    map
end

function Base.haskey(map::HostMap{K,V}, idx) where {K,V}
    key = ZeroInitRef(K, idx)
    value = Ref{V}()
    key_ptr = Base.unsafe_convert(Ptr{K}, key)
    value_ptr = Base.unsafe_convert(Ptr{V}, value)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    value_ptr,
                                    0))
    GC.@preserve key value begin
        bpf(API.BPF_MAP_LOOKUP_ELEM, attr) == 0
    end
end

function nextkey(map::HostMap{K,V}, idx) where {K,V}
    key = ZeroInitRef(K, idx)
    nkey = Ref{K}()
    key_ptr = Base.unsafe_convert(Ptr{K}, key)
    nkey_ptr = Base.unsafe_convert(Ptr{K}, nkey)
    attr = Ref(map_access_elem_attr(map.fd,
                                    key_ptr,
                                    nkey_ptr,
                                    0))
    ret = GC.@preserve key nkey begin
        bpf(API.BPF_MAP_GET_NEXT_KEY, attr)
    end
    if (ret == -1) && (Libc.errno() == Libc.ENOENT)
        return nothing
    end
    ret >= 0 || Base.systemerror(ret)
    nkey[]
end

struct HostMapKeySet{K,V,H<:HostMap{K,V}}
    map::H
end
function Base.length(hk::HostMapKeySet)
    ctr = 0
    for key in hk
        ctr += 1
    end
    ctr
end
Base.keys(map::H) where {K,V,H<:HostMap{K,V}} = HostMapKeySet{K,V,H}(map)
Base.IteratorSize(::Type{<:HostMapKeySet}) = Base.SizeUnknown()
Base.eltype(::Type{HostMapKeySet{K,V,H}}) where {K,V,H} = K
function Base.iterate(hmks::HostMapKeySet{K,V,H}) where {K,V,H}
    fakekey_ref = ZeroInitRef(K)
    realkey = if haskey(hmks.map, fakekey_ref[])
        fakekey_ref[]
    else
        nextkey(hmks.map, fakekey_ref[])
    end
    realkey === nothing && return nothing # empty
    return realkey, (realkey, realkey)
end
function Base.iterate(hmks::HostMapKeySet{K,V,H}, (lastkey, initial)) where {K,V,H}
    lastkey === nothing && return nothing
    nkey = nextkey(hmks.map, lastkey)
    if nkey === nothing || nkey == lastkey
        return nothing
    end
    return nkey, (nkey == initial ? nothing : nkey, initial)
end
Base.length(map::HostMap) = length(keys(map))


text:
# Nest
#
# If the line exceeds the print width it will be nested.
#
# This is done by replacing `PLACEHOLDER` nodes with `NEWLINE`
# nodes and updating the PTree's indent.
#
# `extra_width` provides additional width to consider from
# the top-level node.
#
# Example:
#
# arg1 op arg2
#
# the length of " op" will be considered when nesting arg1
#
# TODO: How to best format undoing a nest?

function walk(f, x::PTree, s::State)
    f(x, s)
    is_leaf(x) && return
    for n in x.nodes
        if n.typ === NEWLINE
            s.line_offset = x.indent
        else
            walk(f, n, s)
        end
    end
end

function reset_line_offset!(x::PTree, s::State)
    !is_leaf(x) && return
    s.line_offset += length(x)
end

function add_indent!(x, s, indent)
    indent == 0 && return
    lo = s.line_offset
    f = (x::PTree, s::State) -> x.indent += indent
    walk(f, x, s)
    s.line_offset = lo
end

function nest!(nodes::Vector{PTree}, s::State, indent; extra_width = 0)
    for n in nodes
        if n.typ === NEWLINE
            s.line_offset = indent
        else
            nest!(n, s, extra_width = extra_width)
        end
    end
end

function nest!(x::PTree, s::State; extra_width = 0)
    if is_leaf(x)
        s.line_offset += length(x)
        return
    end

    if x.typ === CSTParser.Import
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Export
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Using
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.ImportAll
        n_import!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.WhereOpCall
        n_wherecall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.ConditionalOpCall
        n_condcall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.BinaryOpCall
        n_binarycall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Curly
        n_curly!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Call
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.MacroCall
        n_macrocall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Ref
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.ChainOpCall
        n_chainopcall!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Comparison
        n_comparison!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.TupleH
        n_tuple!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Vect
        n_vect!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Vcat
        n_vect!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.TypedVcat
        n_call!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.StringH
        n_stringh!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Parameters
        n_params!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.Braces
        n_braces!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.InvisBrackets
        n_invisbrackets!(x, s, extra_width = extra_width)
    elseif x.typ === CSTParser.UnaryOpCall && x.nodes[2].typ === CSTParser.OPERATOR
        nest!(x.nodes[1], s, extra_width = extra_width + length(x.nodes[2]))
        nest!(x.nodes[2], s)
    else
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE && x.nodes[i+1].typ === CSTParser.Block
                s.line_offset = x.nodes[i+1].indent
            elseif n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === NOTCODE && x.nodes[i+1].typ === CSTParser.Block
                s.line_offset = x.nodes[i+1].indent
            elseif is_leaf(n)
                s.line_offset += length(n)
            else
                nest!(n, s, extra_width = extra_width)
            end
        end
    end
end

function n_invisbrackets!(x, s; extra_width = 0)
    for (i, n) in enumerate(x.nodes)
        if n.typ === NEWLINE && x.nodes[i+1].typ === CSTParser.Block
            s.line_offset = x.nodes[i+1].indent
        elseif n.typ === NEWLINE
            s.line_offset = x.indent
        elseif n.typ === NOTCODE && x.nodes[i+1].typ === CSTParser.Block
            s.line_offset = x.nodes[i+1].indent
        elseif is_leaf(n)
            s.line_offset += length(n)
        elseif i == length(x.nodes) - 1
            nest!(n, s, extra_width = extra_width + 1)
        else
            nest!(n, s, extra_width = extra_width)
        end
    end
end

# Import,Using,Export,ImportAll
function n_import!(x, s; extra_width = 0)
    line_width = s.line_offset + length(x) + extra_width
    idx = findfirst(n -> n.typ === PLACEHOLDER, x.nodes)
    if idx !== nothing && (line_width > s.margin || x.force_nest)
        sidx = findfirst(n -> is_colon(n), x.nodes)
        if sidx === nothing
            sidx = 2
        else
            # for the WHITESPACE node after the colon
            sidx += 1
        end
        x.indent = s.line_offset + sum(length.(x.nodes[1:sidx]))
        s.line_offset = x.indent

        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            else
                nest!(n, s)
            end
        end
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_tuple!(x, s; extra_width = 0)
    line_width = s.line_offset + length(x) + extra_width
    idx = findlast(n -> n.typ === PLACEHOLDER, x.nodes)
    opener = is_opener(x.nodes[1])
    # @info "ENTERING" x.typ s.line_offset length(x) extra_width
    if idx !== nothing && (line_width > s.margin || x.force_nest)
        # @debug "ENTERING" x.indent s.line_offset x.typ
        if opener
            x.nodes[end].indent = x.indent
        end
        line_offset = s.line_offset

        x.indent += s.indent_size
        if x.indent - s.line_offset > 1
            x.indent = s.line_offset
            if opener
                x.indent += 1
                x.nodes[end].indent = s.line_offset
            end
        end

        # @debug "DURING" x.indent s.line_offset x.typ
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA
                x.nodes[i].val = ","
                x.nodes[i].len = 1
                nest!(x.nodes[i], s)
            elseif opener && (i == 1 || i == length(x.nodes))
                nest!(n, s)
            else
                diff = x.indent - x.nodes[i].indent
                add_indent!(n, s, diff)
                nest!(n, s, extra_width = 1)
            end
        end

        if opener
            s.line_offset = x.nodes[end].indent + 1
        end
        # @info "EXITING" x.typ s.line_offset x.indent x.nodes[end].indent
    else
        opener && (extra_width += 1)
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_stringh!(x, s; extra_width = 0)
    # The indent of StringH is set to the the offset
    # of when the quote is first encountered in the source file.

    # This difference notes if there is a change due to nesting.
    diff = x.indent - s.line_offset

    # The new indent for the string is index of when a character in
    # the multiline string is FIRST encountered in the source file - the above difference
    x.indent = max(x.nodes[1].indent - diff, 0)
    nest!(x.nodes, s, x.indent, extra_width = extra_width)
end

n_braces!(x, s; extra_width = 0) = n_tuple!(x, s, extra_width = extra_width)
n_vect!(x, s; extra_width = 0) = n_tuple!(x, s, extra_width = extra_width)
n_params!(x, s; extra_width = 0) = n_tuple!(x, s, extra_width = extra_width)

function n_call!(x, s; extra_width = 0)
    line_width = s.line_offset + length(x) + extra_width
    idx = findlast(n -> n.typ === PLACEHOLDER, x.nodes)
    # @info "ENTERING" x.typ s.line_offset length(x) extra_width
    if idx !== nothing && (line_width > s.margin || x.force_nest)
        x.nodes[end].indent = x.indent
        line_offset = s.line_offset

        caller_len = length(x.nodes[1])

        x.indent += s.indent_size
        # @debug "ENTERING" x.indent s.line_offset x.typ
        if x.indent - s.line_offset > caller_len + 1
            x.indent = s.line_offset + caller_len + 1
            x.nodes[end].indent = s.line_offset
        end

        # @debug "DURING" x.indent s.line_offset x.typ
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif n.typ === PLACEHOLDER
                x.nodes[i] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA
                x.nodes[i].val = ","
                x.nodes[i].len = 1
                nest!(x.nodes[i], s)
            elseif i == 1 || i == length(x.nodes)
                n.typ === CSTParser.Parameters && (n.force_nest = true)
                nest!(n, s, extra_width = 1)
            else
                diff = x.indent - x.nodes[i].indent
                add_indent!(n, s, diff)
                n.typ === CSTParser.Parameters && (n.force_nest = true)
                nest!(n, s, extra_width = 1)
            end
        end

        s.line_offset = x.nodes[end].indent + 1
    else
        is_closer(x.nodes[end]) && (extra_width += 1)
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

n_curly!(x, s; extra_width = 0) = n_call!(x, s, extra_width = extra_width)
n_macrocall!(x, s; extra_width = 0) = n_call!(x, s, extra_width = extra_width)

# "A where B"
function n_wherecall!(x, s; extra_width = 0)
    line_width = s.line_offset + length(x) + extra_width
    # @debug "" s.line_offset x.typ line_width extra_width length(x)
    if line_width > s.margin || x.force_nest
        line_offset = s.line_offset
        # after "A where "
        idx = findfirst(n -> n.typ === PLACEHOLDER, x.nodes)
        Blen = sum(length.(x.nodes[idx+1:end]))

        # A, +7 is the length of " where "
        ew = Blen + 7 + extra_width
        nest!(x.nodes[1], s, extra_width = ew)

        for (i, n) in enumerate(x.nodes[2:idx-1])
            nest!(n, s)
        end

        # "B"

        has_braces = is_closer(x.nodes[end])
        if has_braces
            x.nodes[end].indent = x.indent
        end

        over = (s.line_offset + Blen + extra_width > s.margin) || x.force_nest
        # line_offset = s.line_offset
        x.indent += s.indent_size

        last_typ::Union{CSTParser.Head,Nothing} = nothing
        for (i, n) in enumerate(x.nodes[idx+1:end])
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif is_opener(n)
                if x.indent - s.line_offset > 1
                    x.indent = s.line_offset + 1
                    x.nodes[end].indent = s.line_offset
                end
                nest!(n, s)
            elseif n.typ === PLACEHOLDER && over
                x.nodes[i+idx] = Newline()
                s.line_offset = x.indent
            elseif n.typ === TRAILINGCOMMA && over
                x.nodes[i+idx].val = ","
                x.nodes[i+idx].len = 1
                nest!(x.nodes[i+idx], s)
            elseif has_braces
                nest!(n, s, extra_width = 1 + extra_width)
            else
                nest!(n, s, extra_width = extra_width)
            end
            (!is_leaf(n) || n.typ === CSTParser.IDENTIFIER) && (last_typ = n.typ)
        end

        # @info "" s.line_offset x.typ extra_width has_braces

        # Properly reset line offset in the case the last
        # argument is an IDENTIFIER.
        # if over && has_braces && last_typ === CSTParser.IDENTIFIER && x.nodes[end-2].typ === TRAILINGCOMMA
        if over &&
           has_braces &&
           last_typ === CSTParser.IDENTIFIER && x.nodes[end-2].typ === TRAILINGCOMMA
            s.line_offset = x.nodes[end].indent + 1
        end

        # @info "" s.line_offset x.typ extra_width has_braces
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

function n_chainopcall!(x, s; extra_width = 0)
    line_width = s.line_offset + length(x) + extra_width
    # @info "ENTERING" x.typ extra_width s.line_offset length(x)
    if line_width > s.margin || x.force_nest
        line_offset = s.line_offset
        x.indent = s.line_offset
        phs = reverse(findall(n -> n.typ === PLACEHOLDER, x.nodes))
        for (i, idx) in enumerate(phs)
            if i == 1
                x.nodes[idx] = Newline()
            else
                nidx = phs[i-1]
                l1 = sum(length.(x.nodes[1:idx-1]))
                l2 = sum(length.(x.nodes[idx:nidx-1]))
                width = line_offset + l1 + l2
                if x.force_nest || width > s.margin
                    x.nodes[idx] = Newline()
                end
            end
        end

        for (i, n) in enumerate(x.nodes)
            if i == length(x.nodes)
                nest!(n, s, extra_width = extra_width)
            elseif x.nodes[i+1].typ === WHITESPACE
                ew = length(x.nodes[i+1]) + length(x.nodes[i+2])
                nest!(n, s, extra_width = ew)
            elseif n.typ === NEWLINE
                s.line_offset = x.indent
            else
                nest!(n, s)
            end
        end

        s.line_offset = line_offset
        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE && !is_comment(x.nodes[i+1]) && !is_comment(x.nodes[i-1])
                # +1 for newline to whitespace conversion
                width = s.line_offset + 1
                if i == length(x.nodes) - 1
                    width += sum(length.(x.nodes[i+1:end])) + extra_width
                else
                    width += sum(length.(x.nodes[i+1:i+3]))
                end
                # @debug "" s.line_offset l  s.margin
                if width <= s.margin
                    x.nodes[i] = Whitespace(1)
                else
                    s.line_offset = x.indent
                end
            end
            s.line_offset += length(x.nodes[i])
        end

        s.line_offset = line_offset
        walk(reset_line_offset!, x, s)
    else
        nest!(x.nodes, s, x.indent, extra_width = extra_width)
    end
end

n_comparison!(x, s; extra_width = 0) = n_chainopcall!(x, s, extra_width = extra_width)
n_condcall!(x, s; extra_width = 0) = n_chainopcall!(x, s, extra_width = extra_width)

# arg1 op arg2
#
# nest in order of
#
# arg1 op
# arg2
function n_binarycall!(x, s; extra_width = 0)
    # If there's no placeholder the binary call is not nestable
    idxs = findall(n -> n.typ === PLACEHOLDER, x.nodes)
    line_width = s.line_offset + length(x) + extra_width
    # @info "ENTERING" x.typ extra_width s.line_offset length(x) idxs
    if length(idxs) == 2 && (line_width > s.margin || x.force_nest)
        line_offset = s.line_offset
        i1 = idxs[1]
        i2 = idxs[2]
        x.nodes[i1] = Newline()

        has_eq = CSTParser.defines_function(x.ref[]) || nest_assignment(x.ref[])
        if has_eq
            s.line_offset = x.indent + s.indent_size
            x.nodes[i2] = Whitespace(s.indent_size)
            add_indent!(x.nodes[end], s, s.indent_size)
        else
            x.indent = s.line_offset
        end

        # arg2
        nest!(x.nodes[end], s, extra_width = extra_width)

        # arg1 op
        s.line_offset = line_offset
        # extra_width for " op"
        ew = length(x.nodes[2]) + length(x.nodes[3])
        # @debug "DURING" ew s.line_offset x.typ
        nest!(x.nodes[1], s, extra_width = ew)
        for n in x.nodes[2:i1]
            nest!(n, s)
        end

        # @debug "BEFORE RESET" x.indent s.line_offset x.typ extra_width x.nodes[i2-2] length(x.nodes[i2+1])
        # Undo nest if possible, +1 for whitespace
        if !x.force_nest
            line_width = s.line_offset + length(x.nodes[i2+1]) + extra_width + 1
            if line_width <= s.margin
                x.nodes[i1] = Whitespace(1)
                if has_eq
                    x.nodes[i2] = Placeholder(0)
                    add_indent!(x.nodes[end], s, -s.indent_size)
                end
            end
        end

        walk(reset_line_offset!, x, s)
    else
        # Handles the case of a function def defined
        # as "foo(a)::R where {A,B} = body".
        #
        # In this case instead of it being parsed as:
        #
        # CSTParser.BinaryOpCall
        #  - CSTParser.WhereOpCall
        #  - OP
        #  - arg2
        #
        # It's parsed as:
        #
        # CSTParser.BinaryOpCall
        #  - CSTParser.BinaryOpCall
        #   - arg1
        #   - OP
        #   - CSTParser.WhereOpCall
        #    - R
        #    - ...
        #  - OP
        #  - arg2
        #
        # The result being extra width is trickier to deal with.
        idx = findfirst(n -> n.typ === CSTParser.WhereOpCall, x.nodes)
        return_width = 0
        if idx !== nothing && idx > 1
            return_width = length(x.nodes[idx].nodes[1]) + length(x.nodes[2])
        elseif idx === nothing
            return_width, _ = length_to_newline(x, 2)
        end

        # @info "" return_width
        # @debug "" s.line_offset return_width length(x.nodes[1])

        for (i, n) in enumerate(x.nodes)
            if n.typ === NEWLINE
                s.line_offset = x.indent
            elseif i == 1
                nest!(n, s, extra_width = return_width)
            elseif i == idx
                nest!(n, s, extra_width = extra_width)
            else
                nest!(n, s, extra_width = extra_width)
            end
        end
    end
end