text:
const DEFAULT_BOOLEAN_STRINGS = ["T", "F", "t", "f", "TRUE", "FALSE", "true", "false"]
const DEFAULT_TRUE_STRINGS = ["T", "t", "TRUE", "true"]
const DEFAULT_FALSE_STRINGS = ["F", "f", "FALSE", "false"]

const DEFAULT_QUOTATION_CHARACTER = '"'
const DEFAULT_SEPARATOR = ','

const DEFAULT_MISSINGNESS_INDICATORS = ["", "NA", "#NA", "N/A", "#N/A", "NULL", "."]

function parse_bool(x::String)
  if contains(DEFAULT_TRUE_STRINGS, x)
    return true
  elseif contains(DEFAULT_FALSE_STRINGS, x)
    return false
  else
    error("Could not parse bool")
  end
end

##############################################################################
#
# Low-level text parsing
#
##############################################################################

let extract_cache = memio(500, false)
    global extract_string
    function extract_string(this, left::Int, right::Int, omitlist::Set)
        if right - left > length(extract_cache.ios)
            extract_cache_size = right - left
            extract_cache = memio(extract_cache_size,false)
        end
        seek(extract_cache, 0) # necessary?
        if length(this) >= 1
            while isvalid(this, right) && right > left && this[right] == ' '
                right -= 1
            end
            i = left
            while i <= right
                lasti = i
                ch, i = next(this, i)
                if !has(omitlist, lasti)
                    print(extract_cache, ch)
                end
            end
            return takebuf_string(extract_cache)
        else
            return ""
        end
    end
end

const NULLSET = Set()
extract_string(this, left::Int, right::Int) = extract_string(this, left, right, NULLSET)

const STATE_EXPECTING_VALUE = 0
const STATE_IN_BARE = 1
const STATE_IN_QUOTED = 2
const STATE_POSSIBLE_EOQUOTED = 3
const STATE_EXPECTING_SEP = 4

# Read one line of delimited text
# This is complex because delimited text can contain EOL inside quoted fields
function read_separated_line(io,
                             separator::Char,
                             quotation_character::Char)
    # Indexes into the current line for the current item
    left = 0
    right = 0

    # Was using RopeString for efficient appends, but rare case and makes
    # UTF-8 processing harder
    this = Base.chomp!(readline(io))

    # Short-circuit on the empty line
    if this == ""
      return Array(UTF8String, 0)
    end

    # 5-state machine. See list of possible states above
    state = STATE_EXPECTING_VALUE

    # Index of characters to remove
    omitlist = Set()

    # Where are we
    i = start(this)
    eol = false

    # Will eventually return a Vector of strings
    num_elems = 0
    ret = Array(ByteString, 0)

    # off we go! use manual loops because this can grow
    while true
        eol = done(this, i)
        if !eol
            this_i = i
            this_char, i = next(this, i)
        end
        if state == STATE_EXPECTING_VALUE
            if eol
                num_elems += 1
                push!(ret, "")
                break
            elseif this_char == ' '
                continue
            elseif this_char == separator
                num_elems += 1
                push!(ret, "")
            elseif this_char == quotation_character
                left = this_i + 1
                state = STATE_IN_QUOTED
            else
                left = this_i
                state = STATE_IN_BARE
            end
        elseif state == STATE_IN_BARE
            if eol
                right = this_i
                num_elems += 1
                push!(ret, extract_string(this, left, right))
                break
            elseif this_char == separator
                right = this_i - 1
                num_elems += 1
                push!(ret, extract_string(this, left, right))
                state = STATE_EXPECTING_VALUE
            else
                continue
            end
        elseif state == STATE_IN_QUOTED
            if eol
                this = string(this, "\n", Base.chomp!(readline(io)))
            elseif this_char == quotation_character
                state = STATE_POSSIBLE_EOQUOTED
            else
                continue
            end
        elseif state == STATE_POSSIBLE_EOQUOTED
            if eol
                right = this_i - 1
                num_elems += 1
                push!(ret, extract_string(this, left, right, omitlist))
                break
            elseif this_char == quotation_character
                add!(omitlist, this_i)
                state = STATE_IN_QUOTED
            elseif this_char == separator
                right = this_i - 2
                num_elems += 1
                push!(ret, extract_string(this, left, right, omitlist))
                empty!(omitlist)
                state = STATE_EXPECTING_VALUE
            elseif this_char == ' '
                right = this_i - 2
                num_elems += 1
                push!(ret, extract_string(this, left, right, omitlist))
                empty!(omitlist)
                state = STATE_EXPECTING_SEP
            else
                error("unexpected character after a quote")
            end
        elseif state == STATE_EXPECTING_SEP
            if eol
                break
            elseif this_char == ' '
                continue
            elseif this_char == separator
                state = STATE_EXPECTING_VALUE
            else
                error("expecting a separator but got something else")
            end
        end
    end
    ret
end

# Read data line-by-line
function read_separated_text(io::IO,
                             nrows::Int,
                             separator::Char,
                             quotation_character::Char)
    # Read one line to determine the number of columns
    i = 1
    sp = read_separated_line(io, separator, quotation_character)
    ncols = length(sp)

    # If the line is blank, return a 0x0 array to signify this
    if ncols == 0
        return Array(UTF8String, 0, 0)
    end

    # Otherwise, allocate an array to store all of the text we'll read
    text_data = Array(UTF8String, nrows, ncols)
    text_data[i, :] = sp

    # Loop until we've read nrows of text or run out of text
    while i < nrows
        sp = read_separated_line(io, separator, quotation_character)
        if length(sp) == ncols
            i += 1
            text_data[i, :] = sp
        else
            break
        end
    end

    # Return as much text as we read
    return text_data[1:i, :]
end

##############################################################################
#
# Inferential steps
#
##############################################################################

function determine_separator{T <: String}(filename::T)
  if ismatch(r"csv$", filename)
    return ','
  elseif ismatch(r"tsv$", filename)
    return '\t'
  elseif ismatch(r"wsv$", filename)
    return ' '
  else
    error("Unable to determine separator used in $filename")
  end
end

function determine_nrows{T <: String}(filename::T, header::Bool)
  total_lines = countlines(filename)
  if header
    return total_lines - 1
  else
    return total_lines
  end
end

function determine_column_names(io::IO,
                                separator::Char,
                                quotation_character::Char,
                                header::Bool)
  seek(io, 0)
  fields = read_separated_line(io, separator, quotation_character)

  if length(fields) == 0
    error("Failed to determine column names from an empty data source")
  end

  column_names = header ? fields : generate_column_names(length(fields))
  seek(io, 0)
  return column_names
end

function convert_to_dataframe{R <: String,
                              S <: String,
                              T <: String}(text_data::Matrix{R},
                                           missingness_indicators::Vector{S},
                                           column_names::Vector{T})
  # Keep a record of number of rows and columns
  nrows, ncols = size(text_data, 1), length(column_names)

  # Short-circuit if the text data is empty
  if nrows == 0
    column_types = {Any for i in 1:ncols}
    return DataFrame(column_types, column_names, 0)
  end

  # Store the columns as a set of DataVector's inside an Array of Any's
  columns = Array(Any, ncols)

  # Convert each column of text into a DataVector of the
  # appropriate type
  dtime = 0.0
  for j in 1:ncols
    is_missing = BitVector(nrows)
    for i in 1:nrows
      value_missing = contains(missingness_indicators, text_data[i, j])
      if value_missing
        text_data[i, j] = utf8("0")
        is_missing[i] = true
      else
        is_missing[i] = false
      end
    end
    values = Array(Int64, nrows)
    try
      for i in 1:nrows
        values[i] = parseint(text_data[i, j])
      end
    catch
      try
        values = Array(Float64, nrows)
        for i in 1:nrows
          values[i] = parsefloat(text_data[i, j])
        end
      catch
        try
          values = Array(Bool, nrows)
          for i in 1:nrows
            values[i] = parse_bool(text_data[i, j])
          end
        catch
          values = text_data[:, j]
        end
      end
    end
    columns[j] = DataArray(values, is_missing)
  end

  # Prepare the DataFrame we'll return
  df = DataFrame(columns, column_names)
  return df
end

##############################################################################
#
# Text input
#
##############################################################################

# Read at most N lines from an IO object
# Then return a minibatch of at most N rows as a DataFrame
# Add column_types, force_types option
function read_minibatch{R <: String,
                        S <: String}(io::IO,
                                     separator::Char,
                                     quotation_character::Char,
                                     missingness_indicators::Vector{R},
                                     column_names::Vector{S},
                                     minibatch_size::Int)
  # Represent data as an array of strings before type conversion
  text_data = read_separated_text(io, minibatch_size, separator, quotation_character)

  # Convert text data to a DataFrame
  return convert_to_dataframe(text_data, missingness_indicators, column_names)
end

# Read an entire data set into a DataFrame from an IO
# TODO: Do only IO-pass through the data
function read_table{R <: String,
                    S <: String}(io::IO,
                                 separator::Char,
                                 quotation_character::Char,
                                 missingness_indicators::Vector{R},
                                 header::Bool,
                                 column_names::Vector{S},
                                 nrows::Int)
  # Return to start of stream
  seek(io, 0)

  # Read first line to remove header in advance
  if header
    readline(io)
  end

  # Represent data as an array of strings before type conversion
  text_data = read_separated_text(io, nrows, separator, quotation_character)

  # Short-circuit if data set is empty except for a header line
  if size(text_data, 1) == 0
    column_types = {Any for i in 1:length(column_names)}
    return DataFrame(column_types, column_names, 0)
  else
    # Convert text data to a DataFrame
    df = convert_to_dataframe(text_data, missingness_indicators, column_names)
    return df
  end
end

function read_table{T <: String}(filename::T)
  # Do inference for missing configuration settings
  separator = determine_separator(filename)
  quotation_character = DEFAULT_QUOTATION_CHARACTER
  missingness_indicators = DEFAULT_MISSINGNESS_INDICATORS
  header = true
  nrows = determine_nrows(filename, header)
  io = open(filename, "r")
  column_names = determine_column_names(io, separator, quotation_character, header)
  df = read_table(io,
                  separator,
                  quotation_character,
                  missingness_indicators,
                  header,
                  column_names,
                  nrows)
  close(io)
  return df
end

##############################################################################
#
# Text output
#
##############################################################################

# Quotation rules
function in_quotes(val::String, quotation_character::Char)
  string(quotation_character, val, quotation_character)
end
function in_quotes(val::Real, quotation_character::Char)
  string(val)
end
function in_quotes(val::Any, quotation_character::Char)
  string(quotation_character, string(val), quotation_character)
end

# TODO: write_table should do more to react to the type of each column
# Need to increase precision of string representation of Float64's
function print_table(io::IO,
                     df::DataFrame,
                     separator::Char,
                     quotation_character::Char,
                     header::Bool)
  n, p = nrow(df), ncol(df)
  if header
    column_names = colnames(df)
    for j in 1:p
      if j < p
        print(io, in_quotes(column_names[j], quotation_character))
        print(io, separator)
      else
        println(io, in_quotes(column_names[j], quotation_character))
      end
    end
  end
  for i in 1:n
    for j in 1:p
      if j < p
        print(io, in_quotes(df[i, j], quotation_character))
        print(io, separator)
      else
        println(io, in_quotes(df[i, j], quotation_character))
      end
    end
  end
end

function print_table(io::IO,
                     df::DataFrame,
                     separator::Char,
                     quotation_character::Char)
  print_table(io, df, separator, quotation_character, true)
end

function print_table(df::DataFrame, separator::Char, quotation_character::Char)
  print_table(OUTPUT_STREAM, df, separator, quotation_character, true)
end

function print_table(df::DataFrame)
    print_table(OUTPUT_STREAM,
                df,
                DEFAULT_SEPARATOR,
                DEFAULT_QUOTATION_CHARACTER,
                true)
end

function write_table(filename::String,
                     df::DataFrame,
                     separator::Char,
                     quotation_character::Char)
  io = open(filename, "w")
  print_table(io, df, separator, quotation_character)
  close(io)
end

# Infer configuration settings from filename
function write_table(filename::String, df::DataFrame)
  separator = determine_separator(filename)
  quotation_character = DEFAULT_QUOTATION_CHARACTER
  write_table(filename, df, separator, quotation_character)
end

##############################################################################
#
# Binary serialization
#
##############################################################################

# Wrappers for serialization
function save(filename, d)
    f = open(filename, "w")
    serialize(f, d)
    close(f)
end

function load_df(filename)
    f = open(filename)
    dd = deserialize(f)()
    close(f)
    return dd
end

# end


text:
"""
    GuidRecursions{TL,TM⁺,TM, Tμ, TH, TLt0, TMt⁺0, Tμt0}

GuidRecursions defines a struct that contains all info required for computing the guiding term and
likelihood (including ptilde term) for a single shape

## Arguments
Suppose t is the specified (fixed) time grid. Then the elements of the struct are:

- `Lt`:     matrices L
- `Mt⁺``:    matrices M⁺ (inverses of M)
- `M`:      matrices M
- `μ``:      vectors μ
- `Ht`:     matrices H, where H = L' M L
- `Lt0`:    L(0) (so obtained from L(0+) after gpupdate step incorporating observation xobs0)
- `Mt⁺0`:   M⁺(0) (so obtained from M⁺(0+) after gpupdate step incorporating observation xobs0)
- `μt0`:    μ(0) (so obtained μ(0+) after gpupdate step incorporating observation xobs0)
"""
struct GuidRecursions{TL,TM⁺,TM, Tμ, TH, TLt0, TMt⁺0, Tμt0}
    Lt::Vector{TL}          # Lt on grid tt
    Mt⁺::Vector{TM⁺}        # Mt⁺ on grid tt
    Mt::Vector{TM}          # Mt on grid tt
    μt::Vector{Tμ}          # μt on grid tt
    Ht::Vector{TH}          # Ht on grid tt
    Lt0::TLt0               # Lt at time 0, after gpupdate step incorporating observation xobs0
    Mt⁺0::TMt⁺0             # inv(Mt) at time 0, after gpupdate step incorporating observation xobs0
    μt0::Tμt0               # μt at time 0, after gpupdate step incorporating observation xobs0

    function GuidRecursions(t,obsinfo)
        μT = zeros(PointF,obsinfo.n)
        Pnt = eltype(obsinfo.ΣT)
        Lt =  [copy(obsinfo.LT) for _ in t]
        Mt⁺ = [copy(obsinfo.ΣT) for _ in t]
        Mt = map(X -> InverseCholesky(lchol(X)),Mt⁺)
        μt = [copy(μT) for _ in t]
        H = obsinfo.LT' * (obsinfo.ΣT * obsinfo.LT )
        Ht = [copy(H) for _ in t]
        Lt0 = copy([obsinfo.L0; obsinfo.LT])

        m = size(obsinfo.Σ0)[1]
        n = size(obsinfo.ΣT)[2]
        if m==0
            Mt⁺0 = copy(obsinfo.ΣT)
        else
            Mt⁺0 = [copy(obsinfo.Σ0) zeros(Pnt,m,n); zeros(Pnt,n,m) copy(obsinfo.ΣT)]
        end
        μt0 = [0*obsinfo.xobs0; copy(μT)]
        new{eltype(Lt), eltype(Mt⁺), eltype(Mt),eltype(μt),eltype(Ht), typeof(Lt0), typeof(Mt⁺0), typeof(μt0)}(Lt, Mt⁺,Mt, μt, Ht, Lt0, Mt⁺0, μt0)
    end
end

"""
    gp_update!(Lt0₊, Mt⁺0₊::Array{Pnt,2}, μt0₊, (L0, Σ0, xobs0), Lt0, Mt⁺0, μt0) where Pnt

Guided proposal update for newly incoming observation at time zero.
Information on new observations at time zero is `(L0, Σ0, xobs0)`
Values just after time zero, `(Lt0₊, Mt⁺0₊, μt0₊)` are updated to time zero, the result being written into `(Lt0, Mt⁺0, μt0)`
"""
function gp_update!(Lt0₊, Mt⁺0₊::Array{Pnt,2}, μt0₊, (L0, Σ0, xobs0), Lt0, Mt⁺0, μt0) where Pnt
    Lt0 .= [L0; Lt0₊]
    m = size(Σ0)[1]
    n = size(Mt⁺0₊)[2]
    if m==0
        Mt⁺0 .= Mt⁺0₊
    else
        Mt⁺0 .= [Σ0 zeros(Pnt,m,n); zeros(Pnt,n,m) Mt⁺0₊]
    end
    μt0 .= [0*xobs0; μt0₊]
    nothing
end

struct Lm  end

"""
    guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, obsinfo; implicit=true, lowrank=false)

Solve backwards recursions in L, M, μ parametrisation on grid t

## Arguments
- `t`: time grid
- `(Lt, Mt⁺, μt)`: containers to write the solutions into
- `Paux`: auxiliary process
- `obsinfo`: of type ObsInfo containing information on the observations
- `implicit`: if true an implicit Euler backwards scheme is used (else explicit forward)

Case `lowrank=true` still gives an error: fixme!
"""
function guidingbackwards!(::Lm, t, (Lt, Mt⁺, μt), Paux, obsinfo, xT; implicit=true, lowrank=false) #FIXME: add lowrank
    Mt⁺[end] .= obsinfo.ΣT
    Lt[end] .= obsinfo.LT
    μt[end] .= μT = zeros(PointF,obsinfo.n)

    B̃ = Bridge.B(0, Paux, xT)          # does not depend on time
    β̃ = vec(Bridge.β(0,Paux, xT))       # does not depend on time
    σ̃T = Matrix(σ̃(0, Paux, xT))
    dt = t[2] - t[1]
    oldtemp = (0.5*dt) * Bridge.outer(Lt[end] * σ̃T)
    if lowrank  # TBA lowrank on σ̃T, and write into σ̃T
        error("not implemented")
    end
    for i in length(t)-1:-1:1
        dt = t[i+1]-t[i]
        if implicit
            Lt[i] .= Lt[i+1]/lu(I - dt* B̃, Val(false)) # should we use pivoting?
        else
            Lt[i] .=  Lt[i+1] * (I + B̃ * dt)
        end
        temp = (0.5dt) * Bridge.outer( (Lt[i]) * σ̃T)
        Mt⁺[i] .= Mt⁺[i+1] + oldtemp + temp
        oldtemp = temp
        μt[i] .= μt[i+1] + 0.5 * (Lt[i] + Lt[i+1]) * β̃ * dt  # trapezoid rule
    end
    (Lt[1], Mt⁺[1], μt[1])
end

"""
    update_guidrec!(Q, obsinfo)

Q::GuidedProposal
obsinfo::ObsInfo

Computes backwards recursion for (L,M⁺,μ), including gp-update step at time 0.
Next, the `guidrec` field or `Q` is updated.
Returns `Q`.
"""
function update_guidrec!(Q, obsinfo)
    for k in 1:obsinfo.nshapes  # for all shapes
        # solve backward recursions;
        Lt0₊, Mt⁺0₊, μt0₊ =  guidingbackwards!(Lm(), Q.tt, (Q.guidrec[k].Lt, Q.guidrec[k].Mt⁺,Q.guidrec[k].μt), Q.aux[k], obsinfo, auxiliary(Q,k).xT)
        # perform gpupdate step at time zero
        gp_update!(Lt0₊, Mt⁺0₊, μt0₊, (obsinfo.L0, obsinfo.Σ0, obsinfo.xobs0),Q.guidrec[k].Lt0, Q.guidrec[k].Mt⁺0, Q.guidrec[k].μt0)
        # compute Cholesky decomposition of Mt at each time on the grid, need to symmetrize gr.Mt⁺; else AHS  gives numerical roundoff errors when mT \neq 0
        S = map(X -> 0.5*(X+X'), Q.guidrec[k].Mt⁺)
        for i in eachindex(S)
            Q.guidrec[k].Mt[i] = InverseCholesky(lchol(S[i]))
            Q.guidrec[k].Ht[i] .= Q.guidrec[k].Lt[i]' * (Q.guidrec[k].Mt[i] * Q.guidrec[k].Lt[i] )
        end
    end
    Q
end

ext: jl


text:
module Oceananigans

export
    # Helper variables and macros for determining if machine is CUDA-enabled.
    HAVE_CUDA,
    @hascuda,

    Architecture,
    CPU,
    GPU,
    device,

    # Planetary Constants
    ConstantsCollection,
    PlanetaryConstants,
    Earth,
    EarthStationary,

    # Grids
    Grid,
    RegularCartesianGrid,

    # Fields
    Field,
    FaceField,
    CellField,
    FaceFieldX,
    FaceFieldY,
    FaceFieldZ,
    EdgeField,
    set!,

    # FieldSets (collections of related fields)
    FieldSet,
    VelocityFields,
    TracerFields,
    PressureFields,
    SourceTerms,
    StepperTemporaryFields,

    # Forcing functions
    Forcing,

    # Equation of state
    LinearEquationOfState,
    δρ,

    # Boundary conditions
    BoundaryConditions,
    BoundaryCondition,
    Default,
    Flux,
    Gradient,

    # Time stepping
    time_step!,
    time_step_kernel!,

    # Poisson solver
    PoissonSolver,
    PoissonSolverGPU,
    init_poisson_solver,
    solve_poisson_3d_ppn,
    solve_poisson_3d_ppn!,
    solve_poisson_3d_ppn_planned!,
    solve_poisson_3d_ppn_gpu!,
    solve_poisson_3d_ppn_gpu_planned!,

    # Model helper structs, e.g. configuration, clock, etc.
    ModelConfiguration,
    Clock,
    Model,

    # Model output writers
    OutputWriter,
    Checkpointer,
    restore_from_checkpoint,
    BinaryOutputWriter,
    NetCDFOutputWriter,
    write_output,
    read_output,

    # Model diagnostics
    Diagnostic,
    run_diagnostic,
    FieldSummary,
    NaNChecker,
    Nusselt_wT,
    Nusselt_Chi,

    # Package utilities
    prettytime

# Standard library modules
using
    Statistics,
    Printf

# Third-party modules
using
    FFTW,
    JLD,
    NetCDF,
    StaticArrays

import
    GPUifyLoops

const HAVE_CUDA = try
    using CUDAdrv, CUDAnative, CuArrays
    true
catch
    false
end

macro hascuda(ex)
    return HAVE_CUDA ? :($(esc(ex))) : :(nothing)
end

abstract type Architecture end
struct CPU <: Architecture end
struct GPU <: Architecture end

device(::CPU) = GPUifyLoops.CPU()
device(::GPU) = GPUifyLoops.CUDA()

@hascuda begin
    println("CUDA-enabled GPU(s) detected:")
    for (gpu, dev) in enumerate(CUDAnative.devices())
        println(dev)
    end
end

# @hascuda CuArrays.allowscalar(false)

abstract type Metadata end
abstract type ConstantsCollection end
abstract type EquationOfState end
abstract type Grid end
abstract type Field end
abstract type FaceField <: Field end
abstract type FieldSet end
abstract type OutputWriter end
abstract type Diagnostic end
abstract type AbstractPoissonSolver end

include("model_configuration.jl")
include("clock.jl")
include("planetary_constants.jl")
include("grids.jl")
include("fields.jl")
include("fieldsets.jl")
include("forcing.jl")

include("operators/operators.jl")

include("boundary_conditions.jl")
include("equation_of_state.jl")
include("poisson_solvers.jl")
include("models.jl")
include("time_steppers.jl")

include("output_writers.jl")
include("diagnostics.jl")

include("utils.jl")

end # module


text:
using PerformFootball
using Test

@testset "Lectura de Datos Básicos" begin
    path = "../src/Data/Info_Basica/Orange Lions V20-2 vs Orange Lions V20-1/The data report of Orange Lions V20-120200122.xlsx"
    @test leer_BasicInfo(path).match_id == "38986AF3-4EC0-4EE0-AD8F-7569BB343ADC"
    @test size(leer_BasicInfo(path).basic_info) == (8, 2)
    @test leer_BasicInfo(path).team_id == "4967"
    @test size(leer_BasicInfo(path).team_info) == (5, 53)
end

@testset "Lectura de Datos de Entrenamiento" begin
    path = "../src/Data/Info_Entrenamiento/Orange Lions V20-1/58D7EBEB71C74AAB897D20493B15D912.xlsx"
    @test leer_TrainInfo(path).train_id == "58D7EBEB71C74AAB897D20493B15D912"
    @test leer_TrainInfo(path).team_id == "4967"
    @test length(leer_TrainInfo(path).train_physical) == 4
    @test length(leer_TrainInfo(path).train_bpm) == 4
end

max_stars_repo_name: vntr-CC/PerformFoorball

max_issues_repo_name: vntr-CC/PerformFoorball

max_forks_repo_name: vntr-CC/PerformFoorball


text:
export NFileCache

"""
    NFileCache <: FileCache

Represents a number-of-files-constrained file cache. Configurable with maximum number of
files in the cache and how to choose which elements to discard.

    NFileCache(root, nfiles, discard_func; predicate=x->true)

Construct a new file cache. Arguments:
 - `root`: the root directory of the cache.
 - `nfiles`: number of files to retain.
 - `discard_func`: function that return the eviction order, see e.g. `DiscardLRU`
   and `DiscardLFU`.
 - `predicate`: function that determines whether a file should be tracked/included
   in the cache. Currently only used for existing files when setting up the file cache.
   The default is to track every existing file.

!!! warn
    This is a potentially destructive operation since the file cache may delete
    files in the `root` directory to fit the given constraints.

Example usage:

    # Create cache that keeps maximum 10 files and releases LRU objects
    fc = NFileCache(root, 10, DiscardLRU())

    # Add a new file:
    path = add!(fc, key)
    open(path, "w") do io
        write(path, filedata)
    end

    # Check to see if that file is there (and also increment its usage counters)
    hit!(fc, key)

    # Delete that file
    delete!(fc, key)
"""
mutable struct NFileCache <: FileCache
    # storage root
    root::String

    # Storage leaves
    entries::Dict{String,CacheEntry}
    # Not directly used by this cache but keeping this field lets us reuse much
    # of the code path between caches
    total_size::UInt64

    # maximum number of entires to keep
    max_entries::Int

    # Sorter function, see SizeConstrainedFileCache
    discard_ordering::Function

    function NFileCache(root::AbstractString, max_entries::Integer, discard_ordering::Function;
                        predicate::Function=x->true)
        fc = new(
            string(root),
            Dict{String,CacheEntry}(),
            UInt64(0),
            max_entries,
            discard_ordering,
        )

        # Always rebuild the cache so that it represents the correct data on-disk
        rebuild!(fc; predicate=predicate)

        # After rebuilding, immediately run a shrink check:
        n = length(fc.entries)
        if n > max_entries
            shrink!(fc, n - max_entries)
        end

        return fc
    end
end

"""
    add!(fc::NFileCache, key::AbstractString, new_size=0)

Reserves space for a new file within the file cache.
If the number of files exceed the maximum number of entries allowed in the cache
the cache delete a file as governed by the discard ordering.
"""
function add!(fc::NFileCache, key::AbstractString, new_size=0)
    # Remove this key from the cache if it already exists.
    delete!(fc, key)

    if length(fc.entries) == fc.max_entries
        shrink!(fc, 1)
    end

    fc.entries[key] = CacheEntry(new_size, time(), 1)
    fc.total_size += new_size

    # Return the filepath so that users can use this in `open()` or whatever
    path = filepath(fc, key)
    mkpath(dirname(path))
    return path
end

"""
    shrink!(fc::NFileCache, files_to_remove)

Shrinks a file cache by the amount given by discarding files within the file cache,
stopping only when the target has been met or the list of keys to discard is exhausted.
Returns the new total size of the file cache.
"""
function shrink!(fc::NFileCache, files_to_remove::Integer)
    # Sort the keys by the user criterion
    keys_to_discard = fc.discard_ordering(fc)
    # Delete elements until we are under our target size
    for _ in 1:files_to_remove
        delete!(fc, popfirst!(keys_to_discard))
    end
    return fc.total_size
end


text:
function _plot_fspec_grid(fs::Dict{Symbol,Dict{Int64,JAXTAM.FFTData}},
        obs_row)
    example = _recursive_first(fs)

    plt_1 = plot(fs; norm=:rms,   save=false, save_csv=true)
    plt_2 = plot(fs; norm=:leahy, save=false)

    plt_3 = plot(fs, norm=:leahy, save=false, freq_lims=(0, 1),     rebin=(:linear, 1),
        logx=false, logy=false, title_append=" - 0 to 1 Hz")
    plt_4 = plot(fs, norm=:leahy, save=false, freq_lims=(1, :end),  rebin=(:linear, 1),
        logx=false, logy=false, title_append=" - 1 to :end Hz")
    plt_4 = plot(fs, norm=:leahy, save=false, freq_lims=(50, :end), rebin=(:linear, 1),
        logx=false, logy=false, title_append=" - 50 to :end Hz")

    Plots.plot(plt_1, plt_2, plt_3, plt_4, layout=grid(4,1), size=(1140,600*4))

    savefig(example, obs_row, example.e_range; plot_name="grid")
end

"""
    report(::Mission, ::DataFrameRow; e_range::Tuple{Float64,Float64}, overwrite::Bool, nuke::Bool, update_masterpage::Bool)

Generates a report for the default energy range

Creates plots for:
    * Lightcurve (+ grouped lightcurves)
    * Periodigram (+ grouped periodograms)
    * Power Spectra
        * :rms, full range, log-rebinned, log-log plot
        * :leahy, full range, log-rebinned, log-log plot
        * :leahy, 0 to 1 Hz, no rebin, linear-linear plot
        * :leahy, 1 to end Hz, no rebin, linear-linear plot
        * :leahy, 50 to end Hz, no regin, linear-linear plot
    * Spectrogram
    * Pulsation search plot

Produces HTML report page

Updates the homepage
"""
function report(mission, obs_row; e_range=_mission_good_e_range(mission), overwrite=false, nuke=false, update_masterpage=true)
    path_jaxtam = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM")
    path_web    = abspath(_obs_path_local(mission, obs_row; kind=:web), "JAXTAM")

    if nuke
        GC.gc() # Required due to Feather.jl loading files lazily, meaning they can't be removed from disk
                # until garbace collection runs and un-lazily-loads them
        ispath(path_jaxtam) ? rm(path_jaxtam, recursive=true) : false
        ispath(path_web)    ? rm(path_web,    recursive=true) : false
    end

    if ismissing(_log_query(mission, obs_row, "meta", :downloaded; surpress_warn=true)) || !_log_query(mission, obs_row, "meta", :downloaded)
        try
            download(mission, obs_row)
        catch err
            if typeof(err) == JAXTAMError
                _log_add(mission, obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
                @warn err
                return nothing
            else
                rethrow(err)
            end
        end
    end

    if !ismissing(JAXTAM._log_query(mission, obs_row, "errors", :read_cl; surpress_warn=true))
        @warn "Error logged at :read_cl stage, no files to be analysed, skipping report gen"
        return nothing
    end

    images = _log_query(mission, obs_row, "images", e_range)

    img_count_groupless = ismissing(images) ? 0 : size(filter(x->ismissing(x[:group]), images), 1)
    # Expect five 'groupless' plots: lightcurve, periodogram, powerspectra, spectrogram, pulsations
    if img_count_groupless < 5 || overwrite
        try
            lc = JAXTAM.lcurve(mission, obs_row, 2.0^0; e_range=e_range)
            JAXTAM.plot(lc; save=true); JAXTAM.plot_groups(lc; save=true, size_in=(1140,400/2))
            pg = JAXTAM.pgram(lc); JAXTAM.plot(pg; save=true);
            pg = JAXTAM.pgram(lc; per_group=true); JAXTAM.plot_groups(pg; save=true, size_in=(1140,600/2));
            lc = nothing; pg = nothing; GC.gc()
        catch err
            if typeof(err) == JAXTAMError
                _log_add(mission, obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
                @warn err
            else
                rethrow(err)
            end
        end

        try
            lc = JAXTAM.lcurve(mission, obs_row, 2.0^-13; e_range=e_range)
            gtis = JAXTAM.gtis(mission, obs_row, 2.0^-13; lcurve_data=lc, e_range=e_range); lc = 0

            fs = JAXTAM.fspec(mission, obs_row, 2.0^-13, 128; gtis_data=gtis, e_range=e_range)
            @info "Plotting fspec grid";    JAXTAM._plot_fspec_grid(fs, obs_row)
            @info "Plotting fspec groups";  JAXTAM.plot_groups(fs; save=true, size_in=(1140,600/2))
            @info "Plotting sgram";         JAXTAM.plot_sgram(fs;  save=true, size_in=(1140,600))
            @info "Plotting pulses";        JAXTAM.plot_pulses_candle(fs; save=true, size_in=(1140,600/2))
            @info "Plotting pulses groups"; JAXTAM.plot_pulses_candle_groups(fs; save=true, size_in=(1140,600/2))
            fs = 0; GC.gc()
        catch err
            if typeof(err) == JAXTAMError
                _log_add(mission, obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
                @warn err
            else
                rethrow(err)
            end
        end

        # Disable second 64 s power spectra plots
        # fs = JAXTAM.fspec(mission_name, obs_row, 2.0^-13, 64)
        # JAXTAM._plot_fspec_grid(fs, obs_row, mission_name, 2.0^-13, "fspec/64.0/", "fspec.png")
        # JAXTAM.plot_groups(fs; save=true, size_in=(1140,600/2))
        # JAXTAM.plot_sgram(fs;  save=true, size_in=(1140,600/2))
        # JAXTAM.plot_pulses_candle(fs; save=true, size_in=(1140,600/2))
        # fs = 0; GC.gc()

        # _call_all_espec(mission, obs_row)
    end

    sp = try
        _webgen_subpage(mission, obs_row; e_range=e_range)
    catch err
        if typeof(err) == JAXTAMError
            _log_add(mission, obs_row, Dict{String,Any}("errors"=>Dict(err.step=>err)))
            @warn err
        else
            rethrow(err)
        end
        return nothing
    end

    if update_masterpage
        webgen_mission(mission)
    end

    return sp
end

function report(mission::Mission, obsid::String; e_range=_mission_good_e_range(mission), overwrite=false, nuke=false, update_masterpage=true)
    obs_row = master_query(mission, :obsid, obsid)

    return report(mission, obs_row; e_range=e_range, overwrite=overwrite, nuke=nuke, update_masterpage=update_masterpage)
end

"""
    report_all(::Mission, ::DataFrameRow; e_ranges=[(0.2,12.0), (2.0,10.0), (0.2,2.0)], overwrite::Bool, nuke::Bool, update_masterpage::Bool)

Calls `report` with three default energy ranges
"""
function report_all(mission::Mission, obs_row::DataFrames.DataFrameRow; e_ranges=[(0.2,12.0), (2.0,10.0), (0.2,2.0)], overwrite=false, nuke=false, update_masterpage=true)
    if nuke
        path_jaxtam = abspath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM")
        path_web    = abspath(_obs_path_local(mission, obs_row; kind=:web), "JAXTAM")

        GC.gc() # Required due to Feather.jl loading files lazily, meaning they can't be removed from disk
                # until garbace collection runs and un-lazily-loads them
        ispath(path_jaxtam) ? rm(path_jaxtam, recursive=true) : false
        ispath(path_web)    ? rm(path_web,    recursive=true) : false
    end

    for e_range in e_ranges
        report(mission, obs_row; e_range=e_range, overwrite=overwrite, update_masterpage=update_masterpage)
    end

    # Run it again to re-generate report pages so that they all have each others links in them
    for e_range in e_ranges
        println("Report path:\n\t$(report(mission, obs_row; e_range=e_range, overwrite=false))")
    end
end

function report_all(mission::Mission, obsid::String; e_ranges=[(0.2,12.0), (2.0,10.0), (0.2,2.0)], overwrite=false, nuke=false, update_masterpage=true)
    obs_row = _master_query(mission, :obsid, obsid)

    return report_all(mission, obs_row; e_ranges=e_ranges, overwrite=overwrite, nuke=nuke, update_masterpage=update_masterpage)
end


text:
using PersistenceDiagrams
using PersistenceDiagrams: Barcode, InfinityLine, ZeroPersistenceLine
using PersistenceDiagrams: dim_str, clamp_death, clamp_persistence, limits, set_default!

using Compat
using Test

using RecipesBase
using RecipesBase: apply_recipe

# Hack to avoid having to import Plots.
RecipesBase.is_key_supported(::Symbol) = true

"""
Idea: apply recipe and check the number of series on plots.
Not a perfect way to test, but at least it makes sure all points are plotted and checks that
there are no errors.
"""
series(args...; kwargs...) = apply_recipe(Dict{Symbol,Any}(kwargs), args...)
n_series(args...; kwargs...) = length(series(args...; kwargs...))

@testset "Helpers" begin
    @test dim_str(PersistenceDiagram([(1, 1)]; dim=0)) == "₀"
    @test dim_str(PersistenceDiagram([(1, 1)]; dim=1990)) == "₁₉₉₀"
    @test dim_str(PersistenceDiagram([(1, 1)])) == "ₓ"

    int1 = PersistenceInterval(3, Inf)
    int2 = PersistenceInterval(1, 2)
    @test clamp_death(int1, 5) == 5
    @test clamp_persistence(int1, 5) == 5
    @test clamp_death(int2, 5) == 2
    @test clamp_persistence(int2, 5) == 1

    diag1 = PersistenceDiagram([(1, 2), (2, 3), (2, Inf)])
    diag2 = PersistenceDiagram([(0, 2), (2, 3), (2, 5)])
    diag3 = PersistenceDiagram([(0, 2), (2, 3), (2, Inf)]; threshold=6)
    diag4 = PersistenceDiagram([(0, 2), (2, 3), (2, 5)]; threshold=7)

    @test limits((diag1,)) == (0, 3 * 1.25, 3 * 1.25)
    @test limits((diag2,)) == (0, 5, 5 * 1.25)
    @test limits((diag1, diag3)) == (0, 6, 6)
    @test limits((diag1, diag4)) == (0, 7, 7)
    @test limits((diag2, diag4)) == (0, 5, 7)
    @test limits((diag2, diag3)) == (0, 6, 6)

    d = set_default!(Dict(:a => 1), :a, 2)
    @test d[:a] == 1
    d = set_default!(Dict(:a => 1), :b, 2)
    @test d[:b] == 2
end

@testset "Recipes" begin
    @testset "InfinityLine" begin
        @test isequal(only(series(InfinityLine, InfinityLine(true))).args, ([NaN],))
        @test only(series(InfinityLine, InfinityLine(true); infinity=5)).args == ([5],)
        @test only(series(InfinityLine, InfinityLine(true); infinity=5)).plotattributes[:seriestype] ==
            :vline
        @test only(series(InfinityLine, InfinityLine(false); infinity=5)).plotattributes[:seriestype] ==
            :hline
    end

    @testset "ZeroPersistenceLine" begin
        @test only(series(ZeroPersistenceLine, ZeroPersistenceLine())).args == (identity,)
        @test only(
            series(ZeroPersistenceLine, ZeroPersistenceLine(); persistence=true)
        ).args == ([0],)
    end

    @testset "Diagram recipe" begin
        diag = PersistenceDiagram([(1, 2), (2, 3), (3, Inf)])
        @test only(series(typeof(diag), diag; letter=:x)).args == ([1, 2, 3],)
        @test only(series(typeof(diag), diag; letter=:y)).args == ([2, 3, Inf],)
        @test only(series(typeof(diag), diag; letter=:y, infinity=4)).args == ([2, 3, 4],)
    end

    @testset "Persistencediagram" begin
        @test only(series(Val{:persistencediagram}, 1:4, 1:4, nothing)).plotattributes[:markerstrokecolor] ==
            :auto
        @test only(series(Val{:persistencediagram}, 1:4, 1:4, nothing)).plotattributes[:seriestype] ==
            :scatter
    end

    @testset "Diagram plot" begin
        diag1 = PersistenceDiagram([(3, Inf), (1, 2), (3, 4)])
        diag2 = PersistenceDiagram([
            PersistenceInterval(1, 2; a=:a, b=:b, c=[1, 2, 3]),
            PersistenceInterval(3, 4; a=:b, b=:c, c=[1, 0]),
        ])
        diag3 = PersistenceDiagram(PersistenceInterval[])

        @test n_series((diag1,)) == 1 + 1 + 1
        @test n_series((diag1,); infinity=5) == 1 + 1 + 1
        @test n_series((diag1,); persistence=true) == 1 + 1 + 1
        @test n_series((diag1, diag2)) == 1 + 1 + 2
        @test n_series((diag1, diag2); infinity=5) == 1 + 1 + 2
        @test n_series((diag1, diag2); persistence=true) == 1 + 1 + 2
        @test n_series((diag3,); persistence=true) == 1 + 1 + 1
        @test n_series((diag3, diag3); persistence=true) == 1 + 1 + 2
    end

    @testset "Matching plot" begin
        diag1 = PersistenceDiagram([(3, 4), (1, 2), (3, 4)])
        diag2 = PersistenceDiagram([(1, 2), (3, 4)])
        match = matching(Bottleneck(), diag1, diag2)

        @test n_series(match) == 1 + 1 + 1 + 2
    end

    @testset "Barcode plot" begin
        diag1 = PersistenceDiagram([(3, Inf), (1, 2), (3, 4)])
        diag2 = PersistenceDiagram([
            PersistenceInterval(1, 2; a=:a, b=:b, c=[1, 2, 3]),
            PersistenceInterval(3, 4; a=:b, b=:c, c=[1, 0]),
        ])
        diag3 = PersistenceDiagram(PersistenceInterval[])

        @test n_series(Barcode(diag1)) == 1 + 1
        @test n_series(Barcode((diag1,)); infinity=5) == 1 + 1
        @test n_series(Barcode((diag1, diag2))) == 1 + 2
        @test n_series(Barcode([diag1, diag2]); infinity=5) == 1 + 2
        @test n_series(Barcode(typeof(diag1)[])) == 1
        @test n_series(Barcode(diag3)) == 2
        @test n_series(Barcode(diag3, diag3)) == 3
    end
end


text:
"""
Various utility functions for using across models, controllers and views.
"""
module Flash
using DocStringExtensionsMock

import Genie

export flash, flash_has_message

"""
$TYPEDSIGNATURES
"""
function init()
  Genie.Sessions.init()
end

"""
$TYPEDSIGNATURES

Returns the `flash` dict object associated with the current HTTP request.
"""
function flash()
  get(Genie.Requests.payload(), Genie.PARAMS_FLASH_KEY, "")
end

"""
$TYPEDSIGNATURES

Stores `value` onto the flash.
"""
function flash(value::Any) :: Nothing
  Genie.Sessions.set!(Genie.Sessions.session(Genie.Requests.payload()), Genie.PARAMS_FLASH_KEY, value)
  Genie.Requests.payload()[Genie.PARAMS_FLASH_KEY] = value

  nothing
end

"""
$TYPEDSIGNATURES

Checks if there's any value on the flash storage
"""
function flash_has_message() :: Bool
  ! isempty(flash())
end

end


text:
using DistributedFactorGraphs
using IncrementalInference

dfg = CloudGraphsDFG{SolverParams}("localhost", 7474, "neo4j", "test",
                            "testUser", "testRobot", "testSession",
                            "description of test session",
                            solverParams=SolverParams())

# Nuke the user
clearUser!!(dfg)

# Add some nodes.
v1 = addVariable!(dfg, :a, ContinuousScalar, labels = [:POSE])
v2 = addVariable!(dfg, :b, ContinuousScalar, labels = [:POSE])
v3 = addVariable!(dfg, :c, ContinuousScalar, labels = [:LANDMARK])
addFactor!(dfg, [:a], Prior(Normal(0,1)))
f1 = addFactor!(dfg, [:a; :b], LinearConditional(Normal(50.0,2.0)) )
f2 = addFactor!(dfg, [:b; :c], LinearConditional(Normal(50.0,2.0)) )

# Pull and solve this graph
dfgLocal = LightDFG{SolverParams}(solverParams=SolverParams())
DistributedFactorGraphs.getSubgraph(dfg, union(ls(dfg), lsf(dfg)), true, dfgLocal)

# Solve it
tree, smtasks = solveTree!(dfgLocal)


text:
export hmm_viterbi

import ..HMMModel

using mlpack._Internal.io

import mlpack_jll
const hmm_viterbiLibrary = mlpack_jll.libmlpack_julia_hmm_viterbi

# Call the C binding of the mlpack hmm_viterbi binding.
function hmm_viterbi_mlpackMain()
  success = ccall((:hmm_viterbi, hmm_viterbiLibrary), Bool, ())
  if !success
    # Throw an exception---false means there was a C++ exception.
    throw(ErrorException("mlpack binding error; see output"))
  end
end

" Internal module to hold utility functions. "
module hmm_viterbi_internal
  import ..hmm_viterbiLibrary

import ...HMMModel

# Get the value of a model pointer parameter of type HMMModel.
function IOGetParamHMMModel(paramName::String)::HMMModel
  HMMModel(ccall((:IO_GetParamHMMModelPtr, hmm_viterbiLibrary), Ptr{Nothing}, (Cstring,), paramName))
end

# Set the value of a model pointer parameter of type HMMModel.
function IOSetParamHMMModel(paramName::String, model::HMMModel)
  ccall((:IO_SetParamHMMModelPtr, hmm_viterbiLibrary), Nothing, (Cstring, Ptr{Nothing}), paramName, model.ptr)
end

# Serialize a model to the given stream.
function serializeHMMModel(stream::IO, model::HMMModel)
  buf_len = UInt[0]
  buf_ptr = ccall((:SerializeHMMModelPtr, hmm_viterbiLibrary), Ptr{UInt8}, (Ptr{Nothing}, Ptr{UInt}), model.ptr, Base.pointer(buf_len))
  buf = Base.unsafe_wrap(Vector{UInt8}, buf_ptr, buf_len[1]; own=true)
  write(stream, buf)
end
# Deserialize a model from the given stream.
function deserializeHMMModel(stream::IO)::HMMModel
  buffer = read(stream)
  HMMModel(ccall((:DeserializeHMMModelPtr, hmm_viterbiLibrary), Ptr{Nothing}, (Ptr{UInt8}, UInt), Base.pointer(buffer), length(buffer)))
end
end # module

"""
    hmm_viterbi(input, input_model; [verbose])

This utility takes an already-trained HMM, specified as `input_model`, and
evaluates the most probable hidden state sequence of a given sequence of
observations (specified as '`input`, using the Viterbi algorithm.  The computed
state sequence may be saved using the `output` output parameter.

For example, to predict the state sequence of the observations `obs` using the
HMM `hmm`, storing the predicted state sequence to `states`, the following
command could be used:

```julia
julia> using CSV
julia> obs = CSV.read("obs.csv")
julia> states = hmm_viterbi(obs, hmm)
```

# Arguments

 - `input::Array{Float64, 2}`: Matrix containing observations,
 - `input_model::HMMModel`: Trained HMM to use.
 - `verbose::Bool`: Display informational messages and the full list of
      parameters and timers at the end of execution.  Default value `false`.

# Return values

 - `output::Array{Int, 2}`: File to save predicted state sequence to.

"""
function hmm_viterbi(input,
                     input_model::HMMModel;
                     verbose::Union{Bool, Missing} = missing,
                     points_are_rows::Bool = true)
  # Force the symbols to load.
  ccall((:loadSymbols, hmm_viterbiLibrary), Nothing, ());

  IORestoreSettings("Hidden Markov Model (HMM) Viterbi State Prediction")

  # Process each input argument before calling mlpackMain().
  IOSetParamMat("input", input, points_are_rows)
  hmm_viterbi_internal.IOSetParamHMMModel("input_model", convert(HMMModel, input_model))
  if verbose !== nothing && verbose === true
    IOEnableVerbose()
  else
    IODisableVerbose()
  end

  IOSetPassed("output")
  # Call the program.
  hmm_viterbi_mlpackMain()

  return IOGetParamUMat("output", points_are_rows)
end


text:
using MH, BenchmarkTools

burn_N = 100_000
SEED = 1111
for N in [100_000, 1_000_000, 10_000_000, 100_000_000, 1_000_000_000]
    println("===================")
    println("N = $N")
    @benchmark mh_serial_naive(Float32(0.0), N, burn_N, ff)
    @benchmark mh_serial_optimized(Float32(0.0), N, burn_N, ff, SEED)

    @benchmark mh_threaded_optimized(Float32(0.0), N, burn_N, ff, SEED)
    @benchmark mh_threaded_naive(Float32(0.0), N, burn_N, ff, SEED)

end


text:
# # SIR model for the spread of COVID-19
# Adapted from the Agents.jl repository:
# https://github.com/JuliaDynamics/Agents.jl/blob/master/examples/sir.jl

# This example illustrates how to use `GraphSpace` and how to model agents on an graph
# (network) where the transition probabilities between each node (position) is not constant.
# ## SIR model

# A SIR model tracks the ratio of Susceptible, Infected, and Recovered individuals within a population.
# Here we add one more category of individuals: those who are infected, but do not know it.
# Transmission rate for infected and diagnosed individuals is lower than infected and undetected.
# We also allow a fraction of recovered individuals to catch the disease again, meaning
# that recovering the disease does not bring full immunity.

# ## Model parameters
# Here are the model parameters, some of which have default values.
# * `Ns`: a vector of population sizes per city. The amount of cities is just `C=length(Ns)`.
# * `β_und`: a vector for transmission probabilities β of the infected but undetected per city.
#   Transmission probability is how many susceptible are infected per day by an infected individual.
#   If social distancing is practiced, this number increases.
# * `β_det`: an array for transmission probabilities β of the infected and detected per city.
#   If hospitals are full, this number increases.
# * `infection_period = 30`: how many days before a person dies or recovers.
# * `detection_time = 14`: how many days before an infected person is detected.
# * `death_rate = 0.02`: the probability that the individual will die after the `infection_period`.
# * `reinfection_probability = 0.05`: The probability that a recovered person can get infected again.
# * `migration_rates`: A matrix of migration probability per individual per day from one city to another.
# * `Is = [zeros(C-1)..., 1]`: An array for initial number of infected but undetected people per city.
#   This starts as only one infected individual in the last city.

# Notice that `Ns, β, Is` all need to have the same length, as they are numbers for each
# city. We've tried to add values to the infection parameters similar to the ones you would hear
# on the news about COVID-19.

# The good thing with Agent based models is that you could easily extend the model we
# implement here to also include age as an additional property of each agent.
# This makes ABMs flexible and suitable for research of virus spreading.

# ## Making the model in Agents.jl
# We start by defining the `PoorSoul` agent type and the ABM
# cd(@__DIR__) #src
using Agents, Random, DataFrames, LightGraphs
using Distributions: Poisson, DiscreteNonParametric
using DrWatson: @dict
using Plots
# gr() # hide

mutable struct PoorSoul <: AbstractAgent
    id::Int
    pos::Int
    days_infected::Int  # number of days since is infected
    status::Symbol  # 1: S, 2: I, 3:R
end

function model_initiation(;
    Ns,
    migration_rates,
    β_und,
    β_det,
    infection_period = 14,
    reinfection_probability = 0.00,
    detection_time = 7,
    death_rate = 0.00,
    Is = [zeros(Int, length(Ns) - 1)..., 1],
    seed = 0,
)

    # Random.seed!(seed)
    @assert length(Ns) ==
    length(Is) ==
    length(β_und) ==
    length(β_det) ==
    size(migration_rates, 1) "length of Ns, Is, and B, and number of rows/columns in migration_rates should be the same "
    @assert size(migration_rates, 1) == size(migration_rates, 2) "migration_rates rates should be a square matrix"

    C = length(Ns)
    ## normalize migration_rates
    migration_rates_sum = sum(migration_rates, dims = 2)
    for c in 1:C
        migration_rates[c, :] ./= migration_rates_sum[c]
    end

    properties = @dict(
        Ns,
        Is,
        β_und,
        β_det,
        β_det,
        migration_rates,
        infection_period,
        reinfection_probability,
        detection_time,
        C,
        death_rate
    )
    space = GraphSpace(complete_digraph(C))
    model = ABM(PoorSoul, space; properties = properties)

    ## Add initial individuals
    for city in 1:C, n in 1:Ns[city]
        ind = add_agent!(city, model, 0, :S) # Susceptible
    end
    ## add infected individuals
    for city in 1:C
        inds = get_node_contents(city, model)
        for n in 1:Is[city]
            agent = model[inds[n]]
            agent.status = :I # Infected
            agent.days_infected = 1
        end
    end
    return model
end
nothing # hide

function agent_step!(agent, model)
    migrate!(agent, model)
    transmit!(agent, model)
    update!(agent, model)
    recover!(agent, model)
end

function migrate!(agent, model)
    pid = agent.pos
    d = DiscreteNonParametric(1:(model.C), model.migration_rates[pid, :])
    m = rand(d)
    if m ≠ pid
        move_agent!(agent, m, model)
    end
end

function transmit!(agent, model)
    rate = if agent.status != :I
        0.0
    elseif agent.days_infected < model.detection_time
        model.β_und[agent.pos]
    else
        model.β_det[agent.pos]
    end

    d = Poisson(rate)
    n = rand(d)
    n == 0 && return

    for contactID in get_node_contents(agent.pos, model)
        contact = model[contactID]
        if contact.status == :S
            contact.status = :I
            n -= 1
            n == 0 && return
        end
    end
end

update!(agent, model) = agent.status == :I && (agent.days_infected += 1)

function recover!(agent, model)
    if agent.days_infected >= model.infection_period
        agent.status = :R
        agent.days_infected = 0
    end
end


text:
const JobId = String
const ValueId = String
const MessageId = String

generated_value_ids = Set()
num_value_ids_issued = 0

function generate_value_id()
    global generated_value_ids
    # v = randstring(8)
    global num_value_ids_issued
    num_value_ids_issued += 1
    v = string(num_value_ids_issued)
    push!(generated_value_ids, v)
    return v
end

generated_message_ids = Set()
num_message_ids_issued = 0

function generate_message_id()
    global generated_message_ids
    # v = randstring(8)
    global num_message_ids_issued
    num_message_ids_issued += 1
    v = string(num_message_ids_issued)
    push!(generated_message_ids, v)
    return v
end

num_bang_values_issued = 0

function generate_bang_value()
    global num_bang_values_issued
    num_bang_values_issued += 1
    v = string(num_bang_values_issued)
    v
end

function get_num_bang_values_issued()
    global num_bang_values_issued
    num_bang_values_issued
end

function set_num_bang_values_issued(new_num_bang_values_issued)
    global num_bang_values_issued
    num_bang_values_issued = new_num_bang_values_issued
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
using SensorFeatureTracking
using Base: Test

@testset "blockMatchingFlow" begin
##
	im1 = zeros(Int32,41,61)
	im1[16:26,16:26] = 100
	im1[16:26,36:46] = 200

	reffeats = map(CartesianIndex{2}, [(16, 36), (26, 36), (16, 46), (26, 46)])

	im2 = zeros(Int32,41,61)
	im2[14:24,14:24] = 70
	im2[18:28,38:48] = 200

	feats = getApproxBestShiTomasi(im1,nfeatures=4, stepguess=0.8)
	@test feats == reffeats

	feats = getApproxBestShiTomasi(im1,nfeatures=8, stepguess=0.8)
	reffeats_im2 = map(CartesianIndex{2},[(14, 14),(24, 14),(14, 24),(24, 24),(18, 38),(28, 38),(18, 48),(28, 48)])

	trackerssd = BlockTracker(deepcopy(feats), search_size = 6)
	block_tracker!(trackerssd, im1, im2)
	@test reffeats_im2 == trackerssd.features

	trackersad = BlockTracker(deepcopy(feats), search_size = 6, matchFunction = compute_sad)
	block_tracker!(trackersad, im1, im2)
	@test reffeats_im2 == trackersad.features

	trackerncc = BlockTracker(deepcopy(feats), search_size = 6, matchFunction = compute_ncc)
	block_tracker!(trackerncc, im1, im2)
	@test reffeats_im2 == trackersad.features
##
end

##

max_stars_repo_name: dehann/SensorFeatureTracking

max_issues_repo_name: dehann/SensorFeatureTracking

max_forks_repo_name: dehann/SensorFeatureTracking


text:
using Revise
using LinearAlgebra
using ControlSystems
using Plots
using Printf
#pgfplots()
plotly()
# constant definitions
ρ_∞ = 1.112
θ_0 = -2
V_∞ = 28
m = 318
S = 6.11
l_t = 4.63
I_yy = 432
S_t = 1.14
l_w = 0.3
C_Lαw = 5.55
C_Lαt = 4.3
C_Lδt = 0.45
b = 15
T_0 = 0
e = 0.8
g = 9.81

function getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy, S_t, l_w, C_Lαw, C_Lαt, C_Lδt, b, T_0, e, g)
    C_Lα = C_Lαw + S_t * C_Lαt / S
    q_∞ = 0.5 * ρ_∞ * V_∞^2

    # state matrix
    A = [
          [
          2 * ( m * g * sind(θ_0) -T_0) / (m * V_∞)
          g * cosd(θ_0) * (1 - 2 * S * C_Lα / (π * b^2 * e) ) / V_∞
          0
          -g * cosd(θ_0)
          ]';

          [
          -2 * g * cosd(θ_0) / V_∞
          (m * g * sind(θ_0) - T_0) / (m * V_∞) - ρ_∞ * V_∞ * S * C_Lα / (2 * m)
          V_∞ - ρ_∞ * V_∞ * S_t * l_t * C_Lαt / (2 * m)
          -g * sind(θ_0)
          ]';

          [
          0
          ρ_∞ * V_∞ * (S * l_w * C_Lαw - S_t * l_t * C_Lαt) / (2 * I_yy)
          -ρ_∞ * V_∞ * S_t * l_t^2 * C_Lαt / (2 * I_yy)
          0
          ]';

          [
          0
          0
          1
          0
          ]'
        ]

    # input matrix
     B = [
            0
            -q_∞ * S_t * C_Lδt / m
            -q_∞ * S_t * l_t * C_Lδt / I_yy
            0
         ]

    # sensory matrix
    C = Matrix(I,4,4)

    # direct term
    D = 0
    return A,B,C,D
end

# Question 1
A,B,C,D = getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy, S_t, l_w, C_Lαw, C_Lαt, C_Lδt, b, T_0, e, g)

# eigenvalue plots
λ = eigvals(A)
scatter(λ, xlabel="Re(λ)", ylabel="Im(λ)", legend = false)

# Question 2
# state space
sys = ss(A,B,C,D)
stepplot(sys, 200, legend = false)
plot!(ylabel = ["V_x (m/s)" "V_z (m/s)" "ω_y (deg/s)" "θ (deg)"], xlabel = ["" "" "" "Time (s)"], title = ["Step Response" "" "" ""], size = [600,600])

phugω = imag(λ[4]) # phugiod frequency
@printf("The phugoid frequency is %f and the time period is %f", phugω, 2 * π / phugω)

# Question 3
λ_C_Lαw = [
           eigvals(getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy, S_t, l_w, C_Lαw*0.8, C_Lαt, C_Lδt, b, T_0, e, g)[1])';
           eigvals(getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy, S_t, l_w, C_Lαw*1.2, C_Lαt, C_Lδt, b, T_0, e, g)[1])'
          ]'

λ_C_Lαt = [
           eigvals(getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy, S_t, l_w, C_Lαw, C_Lαt*0.8, C_Lδt, b, T_0, e, g)[1])';
           eigvals(getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy, S_t, l_w, C_Lαw, C_Lαt*1.2, C_Lδt, b, T_0, e, g)[1])'
          ]'

λ_I_yy = [
           eigvals(getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy*0.8, S_t, l_w, C_Lαw, C_Lαt, C_Lδt, b, T_0, e, g)[1])';
           eigvals(getABCD(ρ_∞, θ_0, V_∞, m, S, l_t, I_yy*1.2, S_t, l_w, C_Lαw, C_Lαt, C_Lδt, b, T_0, e, g)[1])'
          ]'

λ_V_∞ = [
           eigvals(getABCD(ρ_∞, θ_0, V_∞*0.8, m, S, l_t, I_yy, S_t, l_w, C_Lαw, C_Lαt, C_Lδt, b, T_0, e, g)[1])';
           eigvals(getABCD(ρ_∞, θ_0, V_∞*1.2, m, S, l_t, I_yy, S_t, l_w, C_Lαw, C_Lαt, C_Lδt, b, T_0, e, g)[1])'
          ]'

# pole plots
function q3polePlot(λ,λ_ref)
    scatter([λ[:,1][1:2] λ[:,1][3:4]], xlabel="Re(λ)", ylabel="Im(λ)",markercolor = :red, layout = 2, title=["SPPO" "PHUGOID"], labels = ["80%" ""])
    scatter!([λ[:,2][1:2] λ[:,2][3:4]], xlabel="Re(λ)", ylabel="Im(λ)",markercolor = :green, layout = 2, title=["SPPO" "PHUGOID"], labels = ["120%" ""])
    scatter!([λ_ref[1:2] λ_ref[3:4] ], xlabel="Re(λ)", ylabel="Im(λ)",markercolor = :blue, layout = 2, title=["SPPO" "PHUGOID"], labels = ["Unperturbed" ""])
    scatter!(size = [800,400])
end

q3polePlot(λ_C_Lαw, λ)
q3polePlot(λ_C_Lαt, λ)
q3polePlot(λ_I_yy, λ)
q3polePlot(λ_V_∞, λ)

# freq and damp plots

function q3fdPlot(λ, λ_ref, mcolor, mname)
    λ = [λ λ_ref]
    sppo_freq = [sqrt(real(λ[:,i][1])^2 + imag(λ[:,i][1])^2) for i = 1:3]
    phug_freq = [sqrt(real(λ[:,i][3])^2 + imag(λ[:,i][3])^2) for i = 1:3]
    sppo_damp = [abs(real(λ[:,i][1])) / sppo_freq[i] for i = 1:3]
    phug_damp = [abs(real(λ[:,i][3])) / phug_freq[i] for i = 1:3]

    scatter!([sppo_freq[1]/sppo_freq[3] phug_freq[1]/phug_freq[3]], [sppo_damp[1]/sppo_damp[3] phug_damp[1]/phug_damp[3]], xlabel="Natural frequency (rad/s)", ylabel="Damping ratio",markercolor = mcolor, layout = 2, title=["SPPO" "PHUGOID"], markershape = :cross, labels = ["80% "*mname ""])
    scatter!([sppo_freq[2]/sppo_freq[3] phug_freq[2]/phug_freq[3]], [sppo_damp[2]/sppo_damp[3] phug_damp[2]/phug_damp[3]], xlabel="Natural frequency (rad/s)", ylabel="Damping ratio",markercolor = mcolor, layout = 2, title=["SPPO" "PHUGOID"], markershape = :xcross, labels = ["120% "*mname ""])
end

scatter([1 1], [1 1], xlabel="Natural frequency (rad/s)", ylabel="Damping ratio",markercolor = :black, layout = 2, title=["SPPO" "PHUGOID"], labels = ["Unperturbed" ""], size = [800,400])
q3fdPlot(λ_C_Lαw, λ, :blue, "C_Lαw")
q3fdPlot(λ_C_Lαt, λ, :red, "C_Lαt")
q3fdPlot(λ_I_yy, λ, :green, "I_yy")
q3fdPlot(λ_V_∞, λ, :pink, "V_∞")

# Question 4
l_β = (l_w + l_t) / 2 # position of spring
σ = 5:0.1:50
k_β = 0.5 * ρ_∞ * V_∞^2 * S_t * l_t * σ

A_s = -0.5 * ρ_∞ * V_∞ ^2 * S_t * l_β * C_Lαt
A_rs = ρ_∞ * V_∞^2 * S_t * C_Lαt / (2 * m) * [0,1,l_t * m / I_yy,0]
A_sr = 0.5 * ρ_∞ * V_∞ * S_t * l_β * C_Lαt * [0 1 l_t 0]

A_k = [A + A_rs * A_sr / (k_β[i] - A_s) for i = 1:length(k_β)] # state matrix with k spring
λ_k = map((x) -> eigvals(x), A_k)
kLength = length(λ_k)
# get modes

sppo_k = [λ_k[i][2] for i = 1:kLength]
phug_k = [λ_k[i][4] for i = 1:kLength]

sppo_freq = sqrt(real(λ[1])^2 + imag(λ[1])^2)
phug_freq = sqrt(real(λ[3])^2 + imag(λ[3])^2)
sppo_damp = abs(real(λ[1])) / sppo_freq
phug_damp = abs(real(λ[3])) / phug_freq

sppo_freq_k = [sqrt(real(λ_k[i][1])^2 + imag(λ_k[i][1])^2) for i = 1:kLength]
phug_freq_k = [sqrt(real(λ_k[i][3])^2 + imag(λ_k[i][3])^2) for i = 1:kLength]
sppo_damp_k = [abs(real(λ_k[i][1])) / sppo_freq_k[i] for i = 1:kLength]
phug_damp_k = [abs(real(λ_k[i][3])) / phug_freq_k[i] for i = 1:kLength]

# plot poles
scatter([sppo_k phug_k], xlabel="Re(λ)", ylabel="Im(λ)", layout = 2, legend = false, title=["SPPO" "PHUGOID"], size = [800,400], markercolor = :red)
# frequency and damping plots
scatter([sppo_freq_k/sppo_freq phug_freq_k/phug_freq], [sppo_damp_k/sppo_damp phug_damp_k/phug_damp], xlabel="Natural frequency (rad/s)", ylabel="Damping ratio",markercolor = :red, layout = 2, title=["SPPO" "PHUGOID"], legend = false, size = [800,400])

# Question 5a
Q = C

# starting value

R=10^( -3.2 + 0.2 * 1 )
K = lqr(sys, Q, R)
P = feedback(sys,K)
λ_closed = pole(P)
scatter([λ_closed[1:2] λ_closed[3:4]], xlabel="Re(λ)", ylabel="Im(λ)", layout = 2, markercolor = :red, labels = ["closed loop poles" ""], size = [800,400], title = ["SPPO" "PHUGOID"])
for i = 2:25
    R=10^( -3.2 + 0.2 * i )
    K = lqr(sys, Q, R)
    P = feedback(sys,K)
    λ_closed = pole(P)
    scatter!([λ_closed[1:2] λ_closed[3:4]], xlabel="Re(λ)", ylabel="Im(λ)", layout = 2, markercolor = :red, labels = ["" ""])
end
scatter!([λ[1:2] λ[3:4]], xlabel="Re(λ)", ylabel="Im(λ)", layout = 2, markercolor = :blue, labels = ["open loop poles" ""], markershape = :cross)

# Question 5b
R = 1
Q = Diagonal([5,0,0,5])
K = lqr(sys, Q, R)
P = feedback(sys,K)
stepplot(P,200, legend = false)
plot!(ylabel = ["V_x (m/s)" "V_z (m/s)" "ω_y (deg/s)" "θ (deg)"], xlabel = ["" "" "" "Time (s)"], title = ["Step Response" "" "" ""], size = [600,600])
#png("Files/MechanicsOfFlightCoursework/q5bstep.png")

max_stars_repo_name: sprockmonty/MechanicsOfFlightCoursework

max_issues_repo_name: sprockmonty/MechanicsOfFlightCoursework

max_forks_repo_name: sprockmonty/MechanicsOfFlightCoursework


text:
import FileIO:save
using Distributed
using DataFrames, DataFramesMeta, PyPlot
using Dates, Statistics

# Run gridsearch for the Gaussian task (2.3.1)
# (i.e. find optimal parameter values for a gaussian learner)
function run_gridsearch_gaussian(;learnertype = GLeakyIntegrator,
                            learnerparametervalues = collect(eps():0.03:0.9999999),
                            foldername = "gridsearch_gaussian",
                            filename = "leaky",
                            nseeds = 3,
                            # nsteps = 10^6,
                            nsteps = 10^5,
                            # changeprobabilities = [0.001, 0.0001], # With 10^6 steps
                            changeprobabilities = [0.1, 0.05, 0.01, 0.005], # With 10^5 steps
                            sigmas = [5., 2., 1., 0.5, 0.1], # Stays always same, except for robustness
                            Npar = 20)
    #rand(1:typemax(UInt64)-1, nseeds)
    magicseeds = [0x5e8b1e3e9f0e7628, 0xf911cd0e60e4ea6a, 0x7098506127432a75]
    specs = SpecsGaussian(; mu_0 = 0., sigma_0 = 1.,
                    sigmas = sigmas,
                    changeprobabilities = changeprobabilities,
                    nseeds = nseeds, seedenvpolicy = magicseeds,
                    seedlearner = magicseeds,
                    nsteps = nsteps)
    gridsearch(specs, learnertype, learnerparametervalues, foldername, filename, Npar = Npar)
end

function gridsearch(specs::SpecsGaussian, learnertype, learnerparametervalues,
                        foldername::String, filename::String; Npar = 20)

    println("Starting: " * Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS"))
    @show learnertype
    settings = collect(Iterators.product(learnerparametervalues,
                        specs.sigmas, specs.changeprobabilities))
    df = @distributed vcat for setting in settings
        # @show setting
        meanmerrors, stdmerrors, seerrors = grid_core(specs, learnertype,
                                                        setting, Npar = Npar)
        DataFrame(param = setting[1],
                sigma = setting[2],
                changeprobability = setting[3],
                meanmerrors = meanmerrors,
                stdmerrors = stdmerrors,
                seerrors = seerrors)
    end
    println("Finishing: " * Dates.format(Dates.now(), "yyyy-mm-dd_HH-MM-SS"))
    savepath = makehomesavepath(foldername); mkdir(savepath)
    save(joinpath(savepath, foldername * "_" * string(filename)*".jld2"),
        Dict("df" => df,  "learnertype" => string(learnertype),
        "specs" => specs,
        "seedenvpolicy" => specs.seedenvpolicy,
        "seedlearner" => specs.seedlearner))
end

function grid_core(specs::SpecsGaussian, learnertype, setting; Npar = 20)
    merrors = []
    for i in 1:length(specs.seedenvpolicy)
        parameterset = getgridparamset(learnertype, setting, Npar = Npar)
        # @show parameterset
        callback = run!(specs.envtype, learnertype, parameterset,
                        sigma = setting[2],
                        changeprobability = setting[3],
                        nsteps = specs.nsteps,
                        iseedenvpolicy = specs.seedenvpolicy[i],
                        iseedlearner = specs.seedlearner[i]
                        )
        merrortmp = computesquareddiff(callback.true_history[2:end],
                                    callback.estimated_history[2:end])
        push!(merrors, merrortmp)
    end
    merrors = collect(skipmissing(reduce(vcat, merrors)))
    meanmerrors = mean(merrors)
    stdmerrors = std(merrors)
    seerrors = stdmerrors / sqrt(length(merrors))
    meanmerrors, stdmerrors, seerrors
end


text:
"""
    NeuronIndices

Treat vertices as having parameters which represents neurons when formulating the size change problem.

This means that individual indices will be aligned so that the function to the largest extent possible is the same after resizing.
"""
struct NeuronIndices end

Δsizetype(vs::AbstractVector{<:AbstractVertex}) = partialsort!(Δsizetype.(vs), 1; by=Δsizetypeprio, rev=true)
Δsizetype(v::AbstractVertex, seen=Set()) = v in seen ? nothing : Δsizetype(trait(v), v, push!(seen, v))
Δsizetype(t::DecoratingTrait, v, seen) = Δsizetype(base(t), v, seen)
Δsizetype(t::MutationTrait, v, seen) = Δsizetype(t, base(v), seen)
Δsizetype(::MutationTrait, ::InputVertex, seen) = nothing
Δsizetype(::MutationTrait, v::CompVertex, seen) = Δsizetype(v.computation)
Δsizetype(f) = NeuronIndices()

# TODO: Label elemwise and concatenation to allow them to return some low prio Δsizetype? Now things like BatchNorm will not have an opinion on its Δsizetype
# Perhaps better idea: pass down trait to Δsizetype(f) and default to low prio Δsizetype for SizeTransparent? Requires that one always looks at all vertices, but that seems to be the requirement anyways
Δsizetype(::SizeTransparent, v, seen) = partialsort!(vcat(map(vi -> Δsizetype(vi, seen), inputs(v)), map(vo -> Δsizetype(vo, seen), outputs(v))), 1; by=Δsizetypeprio, rev=true)

Δsizetypeprio(::Nothing) = 0
Δsizetypeprio(::ScalarSize) = 1
Δsizetypeprio(::NeuronIndices) = 2

"""
    defaultutility(v::AbstractVertex)

Default function used to calculate utility of output neurons.

Implement either `defaultutility(f)` or `defaultutility(t, f)` where `f` is the computation performed
by `CompVertex` and `t` is [`trait(v)`](@ref) to set the default for `f`.

# Examples

```jldoctest
julia> using NaiveNASlib, Statistics

julia> struct Affine{T}
        W::Matrix{T}
       end;

julia> NaiveNASlib.defaultutility(l::Affine) = mean(abs, l.W; dims=2);

julia> NaiveNASlib.defaultutility(Affine(ones(2,3)))
2×1 Matrix{Float64}:
 1.0
 1.0
```
"""
defaultutility(v::AbstractVertex) = defaultutility(trait(v), v)
defaultutility(t, v::AbstractVertex) = defaultutility(t, base(v))
defaultutility(t, ::InputVertex) = 1
defaultutility(t, v::CompVertex) = defaultutility(t, v.computation)
defaultutility(t, f) = resolve_utility(t, defaultutility(f))

struct NoDefaultUtility end

defaultutility(f) = NoDefaultUtility()

resolve_utility(t::DecoratingTrait, u::NoDefaultUtility) = resolve_utility(base(t),u)
resolve_utility(::MutationTrait, ::NoDefaultUtility) = 1
resolve_utility(::SizeTransparent, ::NoDefaultUtility) = 0 # rationale is that utility of other vertices will be tied to this anyways
resolve_utility(t, val) = val

Δsize!(::Nothing, s::AbstractΔSizeStrategy, vs::AbstractVector{<:AbstractVertex}) = false

# Main entry point for computing and applying size changes with NeuronIndices

Δsize!(case::NeuronIndices, s::AbstractΔSizeStrategy, vs::AbstractVector{<:AbstractVertex}) = Δsize!(defaultutility, case,s , vs)
function Δsize!(utilityfun, case::NeuronIndices, s::AbstractΔSizeStrategy, vs::AbstractVector{<:AbstractVertex})
    success, ins, outs = solve_outputs_selection(s, vs, utilityfun)
    if success
        applyΔsize!(case, s, ins, outs, vs)
    end
    return success
end

function Δsize!(utilityfun, case::NeuronIndices, s::TruncateInIndsToValid, vs::AbstractVector{<:AbstractVertex})
    success, ins, outs = solve_outputs_selection(s.strategy, vs, utilityfun)
    if success
        for (vv, ininds) in ins
            for innr in eachindex(ininds)
                ininds[innr] = aligntomax(nin(vv)[innr], ininds[innr])
            end
            newins, newouts = align_outs_to_ins(vv, ins[vv], outs[vv])
            ins[vv] = newins
            outs[vv] = newouts
        end
        applyΔsize!(case, s, ins, outs, vs)
    end
    return success
end

aligntomax(maxval, ::Missing) = missing
function aligntomax(maxval, arr)
    arr = copy(arr)
    maxind = argmax(arr)
    while arr[maxind] > maxval
        newval = maxval
        while newval in arr && newval > -1
            newval -= 1
        end
        arr[maxind] = newval == 0 ? -1 : newval
        maxind = argmax(arr)
    end
    return arr
end

align_outs_to_ins(v, ins, outs) = align_outs_to_ins(trait(v), v, ins, outs)
align_outs_to_ins(t::DecoratingTrait, v, ins, outs) = align_outs_to_ins(base(t), v, ins, outs)
align_outs_to_ins(t, v, ins, outs) = ins,outs
function align_outs_to_ins(::SizeInvariant, v, ins::AbstractVector, outs)
    isempty(ins) && return ins, outs
    inds = ins[1]
    newins = repeat([inds], length(ins))
    newouts = ismissing(outs) ? outs : inds
    return newins, newouts
end

# Main entry point for applying size changes with NeuronIndices after we have compted them
# Peel off all DecoratingJuMPΔSizeStrategies until we find either an AbstractAfterΔSizeStrategy (which is needed in after_Δnin etc)
# or we hit the bottom (in the form of an AbstractΔSizeStrategy).
function applyΔsize!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, ins::AbstractDict, outs::AbstractDict, vs::AbstractVector{<:AbstractVertex})
    applyΔsize!(case, base(s), ins, outs, vs)
end
function applyΔsize!(case::NeuronIndices, s::AbstractAfterΔSizeStrategy, ins::AbstractDict, outs::AbstractDict, vs::AbstractVector{<:AbstractVertex})
    _applyΔsize!(case, s, ins, outs, vs)
end
function applyΔsize!(case::NeuronIndices, s::AbstractΔSizeStrategy, ins::AbstractDict, outs::AbstractDict, vs::AbstractVector{<:AbstractVertex})
    _applyΔsize!(case, s, ins, outs, vs)
end

function _applyΔsize!(case::NeuronIndices, s, ins::AbstractDict, outs::AbstractDict, vs::AbstractVector{<:AbstractVertex})

    Δnins = map(vs) do vi
        vins = ins[vi]
        ismissing(vins) && return false
        mapreduce(|, zip(vins, nin(vi)); init=false) do (vin, ninorg)
            ismissing(vin) && return false
            return length(vin) != ninorg
        end
    end
    Δnouts = [length(outs[vi]) != nout(vi) for vi in vs]

    for vi in vs
        applyΔsize!(case, s, vi, ins[vi], outs[vi])
    end

    for (i, vi) in enumerate(vs)
        after_Δnin(s, vi, ins[vi], Δnins[i])
        after_Δnout(s, vi, outs[vi], Δnouts[i])
    end
end

# Main entry point for applying size changes to a single vertex with NeuronIndices

function applyΔsize!(case::NeuronIndices, s::WithKwargs, v::AbstractVertex, ins::AbstractVector, outs::AbstractVector; kws...)
    applyΔsize!(case, base(s), v, ins, outs; kws..., s.kwargs...)
end
function applyΔsize!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, v::AbstractVertex, ins::AbstractVector, outs::AbstractVector; kws...)
     applyΔsize!(case, base(s), v, ins, outs; kws...)
end
function applyΔsize!(case::NeuronIndices, s::AbstractΔSizeStrategy, v::AbstractVertex, ins::AbstractVector, outs::AbstractVector; kws...)
    applyΔsize!(case, v, ins, outs; kws...)
end

function applyΔsize!(::NeuronIndices, v::InputSizeVertex, ins::AbstractVector, outs::AbstractVector; kwargs...)
    if !all(isempty, skipmissing(ins))
        throw(ArgumentError("Try to change input neurons of InputVertex $(name(v)) to $ins"))
    end
    if outs != 1:nout(v)
        throw(ArgumentError("Try to change output neurons of InputVertex $(name(v)) to $outs"))
    end
end

function applyΔsize!(case::NeuronIndices, v::AbstractVertex, ins::AbstractVector, outs::AbstractVector; kwargs...)
    applyΔsize!(case, base(v), ins, outs; kwargs...)
end
function applyΔsize!(case::NeuronIndices, v::CompVertex, ins::AbstractVector, outs::AbstractVector; kwargs...)
    # I don't love how we go back to Δsize! from here, but I don't want to export applyΔsize! and I can't think of another name right now...
    Δsize!(case, v.computation, ins, outs; kwargs...)
end
Δsize!(::NeuronIndices, f, ins::AbstractVector, outs::AbstractVector; kwargs...) = Δsize!(f, ins, outs; kwargs...)
"""
    Δsize!(f::F, ins::AbstractVector, outs::AbstractVector; kwargs...)

Apply the changes to `f` so that input neurons in `ins` and output neurons in `outs` are selected and/or inserted.

Argument `outs` is a vector of indices to select/insert while `ins` has one vector of indices per input vertex.

Shall be implemented for any type `F` which holds parameters for which the shape shall be modified by NaiveNASlib.

Tip: the function [`parselect`](@ref) can be used to change parameter arrays according to `ins` and `outs`.

Tip: `kwargs` can be passed using [`WithKwargs`](@ref).
"""
function Δsize!(f, ins::AbstractVector, outs::AbstractVector; kwargs...) end

"""
    parselect(pars::AbstractArray{T,N}, elements_per_dim...; newfun = (T, dim, size...) -> 0) where {T, N}

Return a new array of same type as `pars` which has a subset of the elements of `pars` as well as potentially
new elements.

Which elements to select/insert is determined by `elements_per_dim` which is a `Pair{Int, Vector{Int}}` mapping
dimension (first memeber) to which elements to select/insert in that dimension (second memeber).

For a each `dim=>elems` pair, the following holds: `selectdim(output, dim, i) == selectdim(pars, dim, elems[i])`
if `elems[i]` is positive and `selectdim(output, dim, i) .== newfun(T, dim, size)[j]` if `elems[i]` is the `j:th`
negative value and `size` is `sum(elems .< 0)`.

# Examples
```jldoctest
julia> using NaiveNASlib, NaiveNASlib.Extend

julia> pars = reshape(1:3*5, 3,5)
3×5 reshape(::UnitRange{Int64}, 3, 5) with eltype Int64:
 1  4  7  10  13
 2  5  8  11  14
 3  6  9  12  15

julia> NaiveNASlib.parselect(pars, 1 => [-1, 1,3,-1,2], 2=>[3, -1, 2, 1]; newfun = (T, d, s...) -> fill(-T(d), s))
5×4 Matrix{Int64}:
 -1  -2  -1  -1
  7  -2   4   1
  9  -2   6   3
 -1  -2  -1  -1
  8  -2   5   2
```
"""
function parselect(pars::AbstractArray{T,N}, elements_per_dim...; newfun = (T, dim, size...) -> 0) where {T, N}
    psize = collect(size(pars))
    assign = repeat(Any[Colon()], N)
    access = repeat(Any[Colon()], N)

    elements_per_dim = filter(!ismissing ∘ last, elements_per_dim)

    for (dim, elements) in elements_per_dim
        psize[dim] = length(elements)
    end
    newpars = similar(pars, psize...)

    for (dim, elements) in elements_per_dim
        indskeep = filter(ind -> ind > 0, elements)
        indsmap = elements .> 0
        newmap = .!indsmap

        assign[dim] = findall(indsmap)
        access[dim] = indskeep
        tsize = copy(psize)
        tsize[dim] = sum(newmap)
        selectdim(newpars, dim, newmap) .= newfun(T, dim, tsize...)
    end

    newpars[assign...] = pars[access...]
    return newpars
end
parselect(::Missing, args...;kwargs...) = missing

function solve_outputs_selection(s::LogΔSizeExec, vertices::AbstractVector{<:AbstractVertex}, utilityfun)
    @logmsg s.level s.msgfun(vertices[1])
    return solve_outputs_selection(base(s), vertices, utilityfun)
end

solve_outputs_selection(::ΔSizeFailNoOp, vs::AbstractVector{<:AbstractVertex}, ::Any) = false, Dict(v => [1:n for n in nin(v)] for v in vs), Dict(v => 1:nout(v) for v in vs)
solve_outputs_selection(s::ThrowΔSizeFailError, vertices::AbstractVector{<:AbstractVertex}, ::Any) = throw(ΔSizeFailError(s.msgfun(vertices)))
solve_outputs_selection(s::WithUtilityFun, vs::AbstractVector{<:AbstractVertex}, ::Any) = solve_outputs_selection(s.strategy, vs, s.utilityfun)
"""
    solve_outputs_selection(s::AbstractΔSizeStrategy, vertices::AbstractVector{<:AbstractVertex}, utilityfun)

Returns a tuple `(success, nindict, noutdict)` where `nindict[vi]` are new input neuron indices and `noutdict[vi]` are new output neuron indices for each vertex `vi` in `vertices`.

The function generally tries to maximize `sum(utilityfun(vi) .* selected[vi]) ∀ vi in vertices` where `selected[vi]` is all elements in `noutdict[vi]` larger than 0 (negative values
in `noutdict` indicates a new output shall be inserted at that position). This however is up to the implementation of the `AbstractΔSizeStrategy s`.

Since selection of outputs is not guaranteed to work in all cases, a flag `success` is also returned. If `success` is `false` then applying the new indices may (and probably will) fail.
"""
function solve_outputs_selection(s::AbstractJuMPΔSizeStrategy, vs::AbstractVector{<:AbstractVertex}, orgutilityfun)
    utilityfun = remaputilityfun(orgutilityfun, s, vs)
    model = selectmodel(NeuronIndices(), s, vs, utilityfun)
    case = NeuronIndices()

    # TODO: Cache and rescale utilityfun?

    # The binary variables `outselectvars` tells us which existing output indices to select
    # The integer variables `outinsertvars` tells us where in the result we shall insert -1 where -1 means "create a new output (e.g. a neuron)
    # Thus, the result will consist of all selected indices with possibly interlaced -1s
    outselectvars, outinsertvars, noutdict = createselectvars!(case, s, vs, (;model, utilityfun))

    objexpr = sizeobjective!(case, s, vs, (;model, outselectvars, outinsertvars, noutdict, utilityfun))
    for v in vs
        data = (;model, outselectvars, outinsertvars, noutdict, objexpr, utilityfun)
        vertexconstraints!(case, v, s, data)
        objexpr = selectobjective!(case, s, v, data)
    end

    @objective(model, Max, objexpr)

    JuMP.optimize!(model)

    !accept(case, s, model) && return   let fbstrat = fallback(s)
                                            solve_outputs_selection(fbstrat, add_participants!(fbstrat, copy(vs)), orgutilityfun)
                                        end

    return true, extract_ininds_and_outinds(s, outselectvars, outinsertvars)...
end

remaputilityfun(utilityfun, s::DecoratingJuMPΔSizeStrategy, vs, args...) = remaputilityfun(utilityfun, base(s), vs, args...)
function remaputilityfun(utilityfun, ::AbstractJuMPΔSizeStrategy, vs, floorval=1e-2, ceilval=1e5)
    cache = Dict{eltype(vs), Any}(v => utilityfun(v) for v in vs)

    # For any problem which one thinks can be solved by writing code one will always generate new problems
    # Learn to keep the number of moving parts in a codebase to a minimum or spend an eternity writing
    # code like this
    !any(vs) do v
        issizemutable(v) || return false
        any(>(0), cache[v])
    end && return v -> cache[v]

    # We would prefer the range of the utility values to stay well within typical solver tolerances
    # Not sure how to do this correctly. Some papers point towards scaling the cost/constraint matrix,
    # but it does not seem to be availble from JuMP

    # What this ad-hoc-y algorithm tries to do is to check if the magnitudes are within bounds [floorval, ceilval]
    # If they are not, it tries to adjust them to be so by uniform scaling.
    # The scaling is limited so that it never pushes something out of bounds.
    # If magnitudes outside both bounds exist (i.e both too large and too small values exist) it does not do anything.
    # TODO: Use quantiles to not let outliers prevent scaling? Some optimization criterion for maximizing the fraction of values within bounds?

    # Don't let 0s affect the result or else all those SizeTransparent vertices is going to push us to max scaling all the time.
    minval, maxval = extrema(abs, Iterators.filter(!=(0), Iterators.flatten(val for (v, val) in cache if issizemutable(v))))
    lb = max(minval, floorval)
    ub = min(maxval, ceilval)
    maxscaledown = lb / floorval
    scaledown = min(maxscaledown, maxval / ub)
    maxscaleup = ceilval / ub
    scaleup = min(maxscaleup, lb / minval)

    scal = scaledown > scaleup ? 1 / scaledown : scaleup

    for (v, val) in cache
        if issizemutable(v)
            cache[v] = val .* scal
        end
    end
    return v -> cache[v]
end

createselectvars!(case, s::DecoratingJuMPΔSizeStrategy, vs, data) = createselectvars!(case ,base(s), vs, data)
function createselectvars!(case, s::UpperInsertBound{<:Integer}, vs, data)
    isempty(vs) && return createselectvars!(case, base(s), vs, data)
    createselectvars!(case, UpperInsertBound(vs, s.vbounds; strategy=base(s), fallback=fallback(s)), vs, data)
end
function createselectvars!(case, s::UpperInsertBound, vs, data)
    outselectvars, outinsertvars, noutdict = createselectvars!(case, base(s), vs, data)
    for (v, ub) in s.vbounds
        JuMP.set_upper_bound.(outinsertvars[v], ub)
    end
    return outselectvars, outinsertvars, noutdict
end
function createselectvars!(case, ::AbstractJuMPΔSizeStrategy, vs, data)
    model = data.model
    utilityfun = data.utilityfun
    # We use the lenght of utilityfun as the number of outputs to select as we will use it to add utility to each variable,
    # so their lenghts must not mismatch
    # They might mismatch when there are parameters attached to some size transparent vertex (e.g BatchNorm or even
    # concatenation with some computed activation utility metric) and we are in the process of aligning sizes after
    # adding/removing an edge.
    # Forcing the user to provide a length of utilityfun matching nout(v) in this case is not really fair as it is we who
    # caused the mismatch.
    # TODO: The indicies alignment below will not be ideal in this case. Perhaps it can be mitigated with some strategy
    # providing information on how the edges have changed so that we can better align what hasn't changed.
    selectsizes = map(vs) do v
        outvals = utilityfun(v)
        return outvals isa Number ? nout(v) : length(outvals)
    end

    outselectvars = Dict(v => @variable(model, [1:nvals], Bin) for (v,nvals) in zip(vs, selectsizes))
    outinsertvars = Dict(v => @variable(model, [1:nout(v)], Int, lower_bound=0) for v in vs)

    # This is the same variable name as used when adjusting size only. It allows us to delegate alot of size changing strategies to ScalarSize.
    # Drawback is that it will yield a fair bit of back and forth so maybe it is too clever for its own good
    noutdict = Dict(v => @expression(model, sum(outselectvars[v]) + sum(outinsertvars[v])) for v in vs)

    return outselectvars, outinsertvars, noutdict
end

selectmodel(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, vs, utilityfun) = selectmodel(case, base(s), vs, utilityfun)
selectmodel(::NeuronIndices, ::AbstractJuMPΔSizeStrategy, vs, utilityfun) = JuMP.Model(JuMP.optimizer_with_attributes(Cbc.Optimizer, "loglevel"=>0, "seconds" => 20))
function selectmodel(case::NeuronIndices, s::TimeLimitΔSizeStrategy, vs, utilityfun)
    model = selectmodel(case, base(s), vs, utilityfun)
    JuMP.set_time_limit_sec(model, s.limit)
    return model
end

function vertexconstraints!(case::NeuronIndices, ::Immutable, v, s::AbstractJuMPΔSizeStrategy, data)
     JuMP.set_lower_bound.(data.outselectvars[v], 1)
     JuMP.set_upper_bound.(data.outinsertvars[v], 0)
     for (vi, outsize) in zip(inputs(v), nin(v))
        @constraint(data.model, data.noutdict[vi] == outsize)
     end
     # This is needed e.g. to apply (impossible) size constraints incase someone wants to change the input size of an output vertex
     vertexconstraints!(case ,s, v, data)
end

# ScalarSize constraints happen to be good here, but we must ensure that vertexconstraints!(::NeuronIndices,...) gets called instead of vertexconstraints!(::ScalarSize,...)
# For AlignNinToNout we could just have broken out the call to vertexconstraints!, but for AlignNinToNoutVertices we need to call it in the middle of the function
vertexconstraints!(::NeuronIndices, v::AbstractVertex, s::AlignNinToNout, data) = vertexconstraints!(ScalarSize(), v, s, data, NeuronIndices())
vertexconstraints!(::NeuronIndices, v::AbstractVertex, s::AlignNinToNoutVertices, data) = vertexconstraints!(ScalarSize(), v::AbstractVertex, s::AlignNinToNoutVertices, data, NeuronIndices())

function vertexconstraints!(case::NeuronIndices, s::AbstractJuMPΔSizeStrategy, v, data)
    insertconstraints!(case, s, v, data)
    sizeconstraint!(case, s, v, data)
    compconstraint!(case, s, v, (data..., vertex=v))
    inoutconstraint!(case, s, v, data)
end

insertconstraints!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, v, data) = insertconstraints!(case, base(s), v, data)
insertconstraints!(::NeuronIndices, ::AbstractJuMPΔSizeStrategy, v, data) = noinsertgaps!(data.model, data.outselectvars[v], data.outinsertvars[v])

"""
    noinsertgaps!(model, select, insert, maxinsert=length(outsel) * 10)

Add constraints so that `insert` does not create undefined gaps in the result of the neuron selection.

Assume `select` is a set of binary variables where `select[i] = 1` means select the output neuron at position `i`
and `insert[i] = N` means insert `N` new output neurons at the position after `i`.

An example of an undefined gap is if `select = [1, 1, 0]` and `insert = [0, 0, 1]` because this results in the
instruction to use existing output neurons `1 and 2` and then insert a new neuron at position `4`.
In this example position `3` is an undefined gap as one should neither put an existing neuron there nor shall
one insert new neurons. Running this method constrains `model` so that this solution is infeasible.
"""
function noinsertgaps!(model, select, insert, maxinsert=max(length(select) * 10, 200)) # TODO: get maxinsert from strategy instead?
    JuMP.set_upper_bound.(filter(!JuMP.has_upper_bound, insert), maxinsert)

    # This is redundant if we don't allow any inserts and the constraints here add significantly to the solution time
    all(var -> JuMP.upper_bound(var) == 0, insert) && return

    # See  https://discourse.julialang.org/t/help-with-constraints-to-select-and-or-insert-columns-to-a-matrix/63654
    insert_nogap = @variable(model, [1:length(insert)], Bin)

    @constraint(model, sum(insert) <= maxinsert)

    # insert[i] == 0 if insert_nogap[i] == 1
    @constraint(model, [i=1:length(insert)], insert[i] <= 2maxinsert * insert_nogap[i])
    # Monotonicity of insert_nogap, i.e insert_nogap[i] can only be 1 if insert_nogap[i+1] is 1
    @constraint(model, [i=2:length(insert)], insert_nogap[i] <= insert_nogap[i-1])
    # Force insert_nogap to have at least as many ones as the number of not selected neurons
    @constraint(model, length(insert) - sum(select) <= sum(insert_nogap))
end

sizeconstraint!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, v, data) = sizeconstraint!(case, base(s), v, data)
sizeconstraint!(::NeuronIndices, s::AbstractJuMPΔSizeStrategy, v, data) = sizeconstraint!(ScalarSize(), DefaultJuMPΔSizeStrategy(), v, data)
sizeconstraint!(case::NeuronIndices, s::ΔNoutMix, v, data) = sizeconstraint!(case, s.exact, v, data)
sizeconstraint!(::NeuronIndices, s::ΔNout{Exact}, v, data) = sizeconstraint!(ScalarSize(), s, v, data)

inoutconstraint!(case, s::DecoratingJuMPΔSizeStrategy, v, data) = inoutconstraint!(case, base(s), v, data)
inoutconstraint!(case, s, v, data) = inoutconstraint!(case, s, trait(v), v, data)
inoutconstraint!(case::NeuronIndices, s, t::DecoratingTrait, v, data) = inoutconstraint!(case, s, base(t), v, data)

function inoutconstraint!(::NeuronIndices, s, ::MutationTrait, v, data) end
function inoutconstraint!(case::NeuronIndices, s, t::SizeTransparent, v, data)
    onemismatch  = inoutconstraint!(case, s, t, v, data.model, data.outselectvars)
    onemismatch |= inoutconstraint!(case, s, t, v, data.model, data.outinsertvars)

    # We need to ensure that total sizes align if there was a mismatch
    if onemismatch
        ninconstraint!(case, s, v, data)
    end
end

function inoutconstraint!(::NeuronIndices, s, t::SizeStack, v, model, vardict::Dict)
    offs = 1
    var = vardict[v]
    onemismatch = false
    for (i, vi) in enumerate(inputs(v))
        var_i = vardict[vi]
        # Sizes mismatch when vertex/edge was removed (or edge added)
        if length(var_i) == nin(v)[i] && offs+length(var_i)-1 <= length(var)
            @constraint(model, var_i .== var[offs:offs+length(var_i)-1])
        else
            onemismatch = true
        end
        offs += length(var_i)
    end
    # Length var != nout means that we need to enforce the SizeStack constraint that sum(nin.(inputs(v))) == nout(v)
    # as this indicates that there is some parameter array which needs to be aligned with the new size
    return onemismatch || length(var) != nout(v)
end

function inoutconstraint!(::NeuronIndices, s, ::SizeInvariant, v, model, vardict::Dict)
    var = vardict[v]
    onemismatch = false
    for vi in inputs(v)
        # Sizes mismatch when vertex/edge was removed (or edge added)
        var_i = vardict[vi]
        if length(var_i) == length(var)
            @constraint(model, var_i .== var)
        else
            onemismatch = true
        end
    end
    return onemismatch
end

# Minus sign because we maximize objective function in case NeuronIndices
sizeobjective!(case::NeuronIndices, s::AbstractJuMPΔSizeStrategy, vertices, data) = -objective!(case, s, vertices, data)

objective!(case::NeuronIndices, s::ΔNoutMix, vertices, data) = objective!(case, s.relax, vertices, data)
objective!(case::NeuronIndices, s::ΔNout{Relaxed}, vertices, data) = noutrelax!(case, s.Δs, vertices, data)

objective!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, vertices, data) = objective!(case, base(s), vertices, data)
function objective!(::NeuronIndices, s::AbstractJuMPΔSizeStrategy, vertices, data)
    # Or else the problem becomes infeasible
    isempty(vertices) && return @expression(data.model, 0)

    scale = map(v -> maximum(abs, data.utilityfun(v)), vertices)
    noutvars = map(v -> data.noutdict[v], vertices)
    sizetargets = map(v -> sizetarget(data.utilityfun, v), vertices)
    # L1 norm prevents change in vertices which does not need to change.
    # Max norm tries to spread out the change so no single vertex takes most of the change.
    return norm!(SumNorm(0.1 => L1NormLinear(), 0.8 => MaxNormLinear()), data.model, @expression(data.model, objective[i=1:length(noutvars)], scale[i] * (noutvars[i] - sizetargets[i])), sizetargets)
end

function sizetarget(utilityfun, v)
    util = utilityfun(v)
    # Negative scalar utility is probably never useful, but lets make it consistent with an all negative array
    util isa Number && return util < 0 ? 1 : nout(v)
    return max(1, nout(v) - count(<(0), util))
end

function selectobjective!(case::NeuronIndices, s::AbstractJuMPΔSizeStrategy, v, data)
    utility = utilityobjective!(case, s, v, data)
    @objective(data.model, Max, utility)
    # Also makes us prefer to not insert
    insertlast = insertlastobjective!(case, s, v, data)
    return @expression(data.model, data.objexpr + utility + insertlast)
end
utilityobjective!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, v, data) = utilityobjective!(case, base(s), v, data)
utilityobjective!(::NeuronIndices, ::AbstractJuMPΔSizeStrategy, v, data) = @expression(data.model, sum(data.utilityfun(v) .* data.outselectvars[v]))

insertlastobjective!(case::NeuronIndices, s::DecoratingJuMPΔSizeStrategy, v, data) = insertlastobjective!(case::NeuronIndices, base(s), v, data)
function insertlastobjective!(::NeuronIndices, s, v, data)
    insvars = data.outinsertvars[v]

    isempty(insvars) && return @expression(data.model, 0)

    preferend = collect(length(insvars) : -1 : 1)
    scale = minimum(abs, data.utilityfun(v)) / sum(preferend)
    return @expression(data.model, -0.1scale*sum(insvars .* preferend))
end

extract_ininds_and_outinds(s::DecoratingJuMPΔSizeStrategy, outselectvars::Dict, outinsertvars::Dict) = extract_ininds_and_outinds(base(s), outselectvars, outinsertvars)
function extract_ininds_and_outinds(s, outselectvars::Dict, outinsertvars::Dict)
    outinds = Dict([v => extract_inds(s, outselectvars[v], outinsertvars[v]) for v in keys(outselectvars)])
    ininds = Dict([v => getall(outinds, inputs(v)) for v in keys(outselectvars)])
    return ininds, outinds
end

function extract_inds(::AbstractJuMPΔSizeStrategy, selectvars::T, insertvars::T) where T <: AbstractVector{JuMP.VariableRef}
    # insertvar is N at indices where a N new output neurons shall be added
    insert = round.(Int, JuMP.value.(insertvars))
    selected = findall(xi -> xi > 0, JuMP.value.(selectvars))
    return join_extracted_inds(selected, insert)
end

function join_extracted_inds(selected, insert)
    result = Vector{Int}(undef, length(selected) + sum(insert))

    # Expected result of algorithm: initialize result = selected, then insert insert[i] -1 elements after each index i in result

    i = 1 # result index
    j = 1 # selected index
    k = 1 # insert index
    while i <= length(result)
        if j <= i && j <= length(selected)
            result[i] = selected[j]
            j += 1
            i += 1
        end
        if k <= length(insert) && insert[k] > 0
            start = i
            stop = i+insert[k]-1
            result[start:stop] .= -1
            i += insert[k]
        end
        k += 1
    end
    return result
end