text:
# This is an experiment in adding indexing to vectors. The idea is
# that if a DataFrame column is indexed, the following will have fast
# lookups:
#     df[:(2 .< col1 .< 7), :]
#     df[:( datecol .>= "2011-01-01" ), "col3"]
#     df[:( col .== "red" ), :]
#
# Keeping indexing with columns has advantages:
#     - Column ordering is less likely to be messed up inadvertantly.
#     - Multiple columns can be indexed; there's no main key grouping.
#     - More flexible logic combinations are possible.
#     - You can mix and match keyed and non-keyed comparisons.
#     - The DataFrame structure and indexing are not disturbed.
#     - It's probably less coding and maintenance.
# The disadvantages are:
#     - It's yet another vector type.
#     - It's somewhat slower than having one overall set of DataFrame
#       keys in a sorted DataFrame. It should still be pretty fast (no
#       speed checks, yet).
#     - You can't do data.table/pandas shortcuts like df["A"] for
#       df[:( keycol .== "A" ), :]. (But, df["A"] is less descriptive if
#       you don't know what df's keys are.)

# An IndexedVector is a pointer to the original column along with an
# indexing vector equal to `order(orig)`. A comparison operation like
# `idv .> 3` returns an Indexer type. The Indexer type includes a
# pointer to the IndexedVector along with a vector of Range1's.
# DataVector's and DataFrame's can be indexed with Indexers. It's fast
# because you're using a slice of the already indexed vector.

# Indexer's can be combined with `|` and `&`. In the case where the
# IndexedVector is the same, the Indexer is reduced or expanded as
# appropriate. This includes: `0 .< idv .< 3` or `idv .== 1 | idv .==
# 4`. If Indexers have different IndexedVectors (like `idv1 .== 1 |
# idv2 .== 1`), then the result is converted to a BitVector.

# We handle NA's by excluding them from the index. `order` puts NA's
# at the front. We use the following function to exclude those. This
# does make the indexing a little trickier as the length of the index
# can be less than the length of the DataArray.

indexorder(x) = order(x)
function indexorder{T}(v::AbstractDataVector{T})
    Nna = sum(isna(v))
    order(v)[Nna + 1 : end]
end

type IndexedVector{T} <: AbstractVector{T}
    x::AbstractVector{T}
    idx::Vector{Int}
end
IndexedVector(x::AbstractVector) = IndexedVector(x, indexorder(x))

ref{T,I<:Real}(v::IndexedVector{T},i::AbstractVector{I}) = IndexedVector(v.x[i])
ref{T}(v::IndexedVector{T},i::Real) = v.x[i]
ref{T}(v::IndexedVector{T}, i) = IndexedVector(v.x[i])
assign(v::IndexedVector, i, val::Real) = IndexedVector(assign(v.x, i, val))
assign(v::IndexedVector, i, val) = IndexedVector(assign(v.x, i, val))

# to make assign in a DataFrame work:
upgrade_vector{T}(v::IndexedVector{T}) = v
function insert_single_column!{T}(df::DataFrame,
                               dv::IndexedVector{T},
                               col_ind::ColumnIndex)
    dv_n, df_n = length(dv), nrow(df)
    if df_n != 0
        if dv_n != df_n
            #dv = repeat(dv, df_n)
            error("New columns must have the same length as old columns")
        end
    end
    if has(df.colindex, col_ind)
        j = df.colindex[col_ind]
        df.columns[j] = dv
    else
        if typeof(col_ind) <: String || typeof(col_ind) <: Symbol
            push!(df.colindex, col_ind)
            push!(df.columns, dv)
        else
            if isnextcol(df, col_ind)
                push!(df.colindex, nextcolname(df))
                push!(df.columns, dv)
            else
                println("Column does not exist: $col_ind")
                error("Cannot assign to non-existent column")
            end
        end
    end
    return dv
end

order(x::IndexedVector) = x.idx
sort(x::IndexedVector) = x.x[x.idx]   # Return regular array?
## sort(x::IndexedVector) = IndexedVector(x.x[x.idx], 1:length(x.x))   # or keep this an IndexedVector, like this?

type Indexer
    r::Vector{Range1}
    iv::IndexedVector
end

function intersect(v1::Vector{Range1}, v2::Vector{Range1})
    # Assumes that the Range1's in each vector are sorted and don't overlap.
    # (Actually, it doesn't, but it should to get more speed.)
    res = Range1[]
    # This does more work than it needs to.
    for r1 in v1
        for r2 in v2
            isoverlap = !((r1[1] > r2[end]) | (r2[1] > r1[end]))
            if isoverlap
                push!(res, max(r1[1], r2[1]):min(r1[end], r2[end]))
            end
        end
    end
    res
end

function union(v1::Vector{Range1}, v2::Vector{Range1})
    # Assumes that the Range1's in each vector are sorted and don't overlap.
    # (Actually, it doesn't, but it should to get more speed.)
    # TODO Check for zero length.
    compare(v1, v2) = (v1[end][end] > v2[end][end] ? v2 : v1,
                       v1[end][end] > v2[end][end] ? v1 : v2)
    res = Range1[]
    v1 = copy(v1) # Destructively operate on these
    v2 = copy(v2)
    while length(v1) > 0 && length(v2) > 0
        # right part, working right to left
        (left, right) = compare(v1, v2)
        while length(right) > 0 && right[end][1] > left[end][end]
            push!(res, pop!(right))
        end
        if length(right) == 0 break; end
        # overlap
        r_end = max(right[end][end], left[end][end])
        overlap = false
        while length(left) > 0 && length(right) > 0 && right[end][end] > left[end][1]
            r_start = min(right[end][1], left[end][1])
            overlap = true
            (left, right) = compare(v1, v2)
            if right[end][1] >= r_start
                pop!(right)
            end
        end
        if overlap
            if length(right) == 0 && length(left) > 0 && r_start <= left[end][1]
                r_start = left[end][1]
                pop!(left)
            end
            push!(res, r_start : r_end)
        end
    end
    # Rest of v1 or v2 (no overlaps here)
    while length(v1) > 0
        push!(res, pop!(v1))
    end
    while length(v2) > 0
        push!(res, pop!(v2))
    end
    reverse(res)
end

function (!)(x::Indexer)   # Negate the Indexer
    res = Range1[1 : x.r[1][1] - 1]
    for i in 1:length(x.r) - 1
        push!(res, x.r[i][end] + 1 : x.r[i + 1][1] - 1)
    end
    push!(res, x.r[end][end] + 1 : length(x.iv.idx))
    Indexer(res, x.iv)
end

function (&)(x1::Indexer, x2::Indexer)
    if is(x1.iv, x2.iv)
        Indexer(intersect(x1.r, x2.r), x1.iv)
    else
        bool(x1) & bool(x2)
    end
end
(&)(x1::BitVector, x2::Indexer) = x1 & bool(x2)
(&)(x1::Indexer, x2::BitVector) = x2 & bool(x1)

function (|)(x1::Indexer, x2::Indexer)
    if is(x1.iv, x2.iv)
        Indexer(union(x1.r, x2.r), x1.iv)
    else
        bool(x1) | bool(x2)
    end
end
(|)(x1::BitVector, x2::Indexer) = x1 | bool(x2)
(|)(x1::Indexer, x2::BitVector) = x2 | bool(x1)

function bool(ix::Indexer)
    res = falses(length(ix.iv.idx))
    for i in ix.iv.idx[[ix.r...]]
        res[i] = true
    end
    res
end

# `ref` -- each Range1 of the Indexer (i.r...) is applied to the indexing vector (i.iv.idx)

ref(x::IndexedVector, i::Indexer) = x[i.iv.idx[[i.r...]]]
ref(x::AbstractVector, i::Indexer) = x[i.iv.idx[[i.r...]]]
ref(x::AbstractDataVector, i::Indexer) = x[i.iv.idx[[i.r...]]]

# df[MultiRowIndex, SingleColumnIndex] => (Sub)?AbstractDataVector
function ref(df::DataFrame, row_inds::Indexer, col_ind::ColumnIndex)
    selected_column = df.colindex[col_ind]
    return df.columns[selected_column][row_inds.iv.idx[[row_inds.r...]]]
end

# df[MultiRowIndex, MultiColumnIndex] => (Sub)?DataFrame
function ref{T <: ColumnIndex}(df::DataFrame, row_inds::Indexer, col_inds::AbstractVector{T})
    selected_columns = df.colindex[col_inds]
    new_columns = {dv[row_inds.iv.idx[[row_inds.r...]]] for dv in df.columns[selected_columns]}
    return DataFrame(new_columns, Index(df.colindex.names[selected_columns]))
end

typealias ComparisonTypes Union(Number, String)   # defined mainly to avoid warnings

# element-wise (in)equality operators
# these may need range checks
# Should these results be sorted? Could be a counting sort.
.=={T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : search_sorted_last(a.x, v, a.idx)], a)
.=={T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : search_sorted_last(a.x, v, a.idx)], a)
.>={T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : length(a.idx)], a)
.<={T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[1 : search_sorted_last(a.x, v, a.idx)], a)
.>={T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[1 : search_sorted_last(a.x, v, a.idx)], a)
.<={T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[search_sorted_first(a.x, v, a.idx) : length(a.idx)], a)
.>{T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[search_sorted_first_gt(a.x, v, a.idx) : length(a.idx)], a)
.<{T<:ComparisonTypes}(a::IndexedVector{T}, v::T) = Indexer(Range1[1 : search_sorted_last_lt(a.x, v, a.idx)], a)
.<{T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[search_sorted_first_gt(a.x, v, a.idx) : length(a.idx)], a)
.>{T<:ComparisonTypes}(v::T, a::IndexedVector{T}) = Indexer(Range1[1 : search_sorted_last_lt(a.x, v, a.idx)], a)

function search_sorted_first_gt{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    res = search_sorted_last(a, x, idx)
    if res == 0 return 1 end
    if res == length(a) && a[idx[res]] != x return(length(a)+1) end
    a[idx[res]] == x ? res + 1 : res
end
function search_sorted_last_lt{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    res = search_sorted_first(a, x, idx)
    if res > length(idx) return length(idx) end
    if res == 1 && a[idx[res]] != x return(0) end
    a[idx[res]] == x ? res - 1 : res
end

function in{T}(a::IndexedVector{T}, y::Vector{T})
    res = a .== y[1]
    for i in 2:length(y)
        res = res | (a .== y[i])
    end
    res
end

between{T}(a::IndexedVector{T}, v1::T, v2::T, ) = Indexer(Range1[search_sorted_first(a.x, v1, a.idx) : search_sorted_last(a.x, v2, a.idx)], a)

size(a::IndexedVector) = size(a.x)
length(a::IndexedVector) = length(a.x)
ndims(a::IndexedVector) = 1
eltype{T}(a::IndexedVector{T}) = T

## print(io, a::IndexedVector) = print(io, a.x)
show(io::IO, a::IndexedVector) = show(io, a.x)
repl_show(io::IO, a::IndexedVector) = repl_show(io, a.x)

function search_sorted_last{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    ## Index of the last value of vector a that is less than or equal to x.
    ## Returns 0 if x is less than all values of a.
    ## idx is an indexing vector equal in length to a that sorts a
    ## @assert length(a) == length(idx)
    lo = 0
    hi = length(idx) + 1
    while lo < hi-1
        i = (lo+hi)>>>1
        if isless(x,a[idx[i]])
            hi = i
        else
            lo = i
        end
    end
    lo
end

function search_sorted_first{I<:Integer}(a::AbstractVector, x, idx::AbstractVector{I})
    ## Index of the first value of vector a that is greater than or equal to x.
    ## Returns length(a) + 1 if x is greater than all values in a.
    ## idx is an indexing vector equal in length to a that sorts a
    ## @assert length(a) == length(idx)
    lo = 0
    hi = length(idx) + 1
    while lo < hi-1
        i = (lo+hi)>>>1
        if isless(a[idx[i]],x)
            lo = i
        else
            hi = i
        end
    end
    hi
end

# the following is needed for show(df)
maxShowLength(v::IndexedVector) = length(v) > 0 ? max([length(_string(x)) for x = v.x]) : 0


text:
# configurations
# --------------

"""
Configurations for abstract interpretation performed by JET.
These configurations will be active for all the entries.

---
- `strict_condition_check::Bool = false` \\
  Enables strict condition check.
  JET reports an error if a condition expression type is "non-boolean". In a case when
    the condition type is `Union`, JET will report if either of union split case is
    non-boolean type, but this can lead to lots of false positive error reports when the
    code is not well-typed, because Julia `Base` defines generic functions that are commonly
    used at a conditional context but also may return "non-boolean" values, e.g.:
    - `!(::Function) -> Function`
    - `!(::Missing) -> Missing`
    - `==(::Missing, ::Any) -> Missing`
    - `==(::Any, ::Missing) -> Missing`
    and thus loosely-typed conditional expression often becomes e.g. `Union{Bool, Missing}`,
    and consequently JET will report it as "non-boolean" type
    (NOTE: in Julia `Missing` is certainly not valid conditional type).
  If this configuration is set to `false`, JET enables an heuristic to avoid those false
    positive error reports and won't report an error if a condition expression type is
    `Union` and either of its union split case is `Function` or `Missing`.

  The effect of this configuration can be described with the following examples:

  * with `strict_condition_check::Bool = false` (default)
    ```julia
    julia> test_f() = Dict('a' => 1, :b => 2) # ::Dict{Any,Int}
    test_f (generic function with 1 method)

    julia> @report_call test_f()
    No errors !
    Dict{Any, Int64}
    ```

  * with `strict_condition_check::Bool = true`
    ```julia
    julia> test_f() = Dict('a' => 1, :b => 2) # ::Dict{Any,Int}
    test_f (generic function with 1 method)

    julia> @report_call strict_condition_check = true test_f()
    ═════ 1 possible error found ═════
    ┌ @ REPL[2]:1 Main.Dict(Main.=>('a', 1), Main.=>(:b, 2))
    │┌ @ dict.jl:125 Base.Dict(ps)
    ││┌ @ dict.jl:129 Base.dict_with_eltype(#308, kv, Base.eltype(kv))
    │││┌ @ abstractdict.jl:539 Base.grow_to!(Base.dict_with_eltype(DT_apply, _5), kv)
    ││││┌ @ dict.jl:145 Base.grow_to!(dest2, itr, st)
    │││││┌ @ dict.jl:159 Base.setindex!(new, v, k)
    ││││││┌ @ dict.jl:383 Base.ht_keyindex2!(h, key)
    │││││││┌ @ dict.jl:328 goto %35 if not Base.isequal(key, Base.getindex(keys, index))
    ││││││││ for 1 of union split cases, non-boolean (Missing) used in boolean context: goto %35 if not Base.isequal(key::Symbol, Base.getindex(keys::Vector{Any}, index::Int64)::Any)::Union{Missing, Bool}
    │││││││└───────────────
    Dict{Any, Int64}
    ```

---
---

You can also configure any of the keyword parameters that `$InferenceParams` or `$OptimizationParams` can take,
e.g. `max_methods::Int = 3`, `union_splitting::Int = 4`.
Listed below are selections of those parameters that can have a potent influence on JET analysis.

---
- `ipo_constant_propagation::Bool = true` \\
  Enables constant propagation in abstract interpretation.
  It is _**highly**_ recommended that you keep this configuration `true` to get reasonable analysis result,
  because constant propagation can cut off lots of false positive errorenous code paths and
  thus produce more accurate and useful analysis results.
---
- `aggressive_constant_propagation::Bool = true` \\
  If `true`, JET will try to do constant propagation more "aggressively".
  It can lead to more accurate analysis as explained above, but also it may incur a performance cost.
  JET by default enables this configuration to get more accurate analysis result.
---
- `unoptimize_throw_blocks::Bool = false` \\
  Turn this on to skip analysis on code blocks that will eventually lead to a `throw` call.
  This configuration improves the analysis performance, but it's better to be turned off
  to get a "proper" analysis result, just because there may be other errors even in those "throw blocks".
---
"""
struct JETAnalysisParams
    strict_condition_check::Bool
    @jetconfigurable function JETAnalysisParams(; strict_condition_check::Bool = false,
                                                  )
        return new(strict_condition_check,
                   )
    end
end

# define these functions just to make them able to accept other JET configrations
JETInferenceParams(; ipo_constant_propagation::Bool        = true,
                     aggressive_constant_propagation::Bool = false,
                     unoptimize_throw_blocks::Bool         = true,
                     max_methods::Int                      = 3,
                     union_splitting::Int                  = 4,
                     apply_union_enum::Int                 = 8,
                     tupletype_depth::Int                  = 3,
                     tuple_splat::Int                      = 32,
                     __jetconfigs...) =
    return InferenceParams(; ipo_constant_propagation,
                             aggressive_constant_propagation,
                             unoptimize_throw_blocks,
                             max_methods,
                             union_splitting,
                             apply_union_enum,
                             tupletype_depth,
                             tuple_splat,
                             )
JETOptimizationParams(; inlining::Bool                = inlining_enabled(),
                        inline_cost_threshold::Int    = 100,
                        inline_nonleaf_penalty::Int   = 1000,
                        inline_tupleret_bonus::Int    = 250,
                        inline_error_path_cost::Int   = 20,
                        max_methods::Int              = 3,
                        tuple_splat::Int              = 32,
                        union_splitting::Int          = 4,
                        unoptimize_throw_blocks::Bool = true,
                        __jetconfigs...) =
    return OptimizationParams(; inlining,
                                inline_cost_threshold,
                                inline_nonleaf_penalty,
                                inline_tupleret_bonus,
                                inline_error_path_cost,
                                max_methods,
                                tuple_splat,
                                union_splitting,
                                unoptimize_throw_blocks,
                                )
# # assert here that they create same objects as the original constructors
@assert JETInferenceParams() == InferenceParams()
@assert JETOptimizationParams() == OptimizationParams()

const LOGGER_LEVEL_KEY     = :JET_LOGGER_LEVEL
const INFO_LOGGER_LEVEL    = 0
const DEBUG_LOGGER_LEVEL   = 1
const LOGGER_LEVELS        = Dict(INFO_LOGGER_LEVEL  => :info,
                                  DEBUG_LOGGER_LEVEL => :debug,
                                  )
const DEFAULT_LOGGER_LEVEL = INFO_LOGGER_LEVEL
const LOGGER_LEVELS_DESC   = let
    map(collect(LOGGER_LEVELS)) do (level, desc)
        if level == DEFAULT_LOGGER_LEVEL
            "`$level` (\"$desc\" level, default)"
        else
            "`$level` (\"$desc\" level)"
        end
    end |> (x->join(x,", "))
end
get_logger_level(io::IO)   = get(io, LOGGER_LEVEL_KEY, DEFAULT_LOGGER_LEVEL)::Int

"""
Logging configurations for JET analysis.

---
- `toplevel_logger::Union{Nothing,IO} = nothing` \\
  If `IO` object is given, it will track JET's toplevel analysis.
  Logging level can be specified with `$(repr(LOGGER_LEVEL_KEY))` `IO` property.
  Currently supported logging levels are either of $(LOGGER_LEVELS_DESC).

  Examples:
  * logs into `stdout`
    ```julia
    julia> report_file(filename; toplevel_logger = stdout)
    ```
  * logs into `io::IOBuffer` with "debug" logger level
    ```julia
    julia> report_file(filename; toplevel_logger = IOContext(io, $(repr(LOGGER_LEVEL_KEY)) => $DEBUG_LOGGER_LEVEL));
    ```
---
- `inference_logger::Union{Nothing,IO} = nothing` \\
  If `IO` object is given, it will track JET's abstract interpretation routine.
  Logging level can be specified with `$(repr(LOGGER_LEVEL_KEY))` `IO` property.
  Currently supported logging levels are either of $(LOGGER_LEVELS_DESC).

  Examples:
  * logs into `stdout`
    ```julia
    report_call(f, args...; inference_logger = stdout)
    ```
  * logs into `io::IOBuffer` with "debug" logger level
    ```julia
    julia> report_call(f, args...; inference_logger = IOContext(io, $(repr(LOGGER_LEVEL_KEY)) => $DEBUG_LOGGER_LEVEL))
    ```
---
!!! tip
    Of course you can specify both `toplevel_logger` and `inference_logger` at the same time like below:
    ```julia
    julia> report_and_watch_file(filename;
                                 toplevel_logger = IOContext(logger_io, $(repr(LOGGER_LEVEL_KEY)) => $DEBUG_LOGGER_LEVEL),
                                 inference_logger = inference_io)
    ```
"""
struct JETLogger
    # logger to track toplevel interpretation
    toplevel_logger::Union{Nothing,IO}
    # logger to track abstract interpretation routine
    inference_logger::Union{Nothing,IO}
    @jetconfigurable function JETLogger(; toplevel_logger::Union{Nothing,IO} = nothing,
                                          inference_logger::Union{Nothing,IO} = nothing,
                                          )
        if isa(toplevel_logger, IO)
            @assert get_logger_level(toplevel_logger) in keys(LOGGER_LEVELS) "toplevel_logger's $LOGGER_LEVEL_KEY should be either of $LOGGER_LEVELS_DESC"
        end
        if isa(inference_logger, IO)
            @assert get_logger_level(inference_logger) in keys(LOGGER_LEVELS) "inference_logger's $LOGGER_LEVEL_KEY should be either of $LOGGER_LEVELS_DESC"
        end
        new(toplevel_logger, inference_logger)
    end
end

# XXX we need to consider world range ?
struct AnalysisResult
    linfo::MethodInstance
    argtypes::Vector{Any}
    overridden_by_const::CC.BitVector
    cache::Vector{InferenceErrorReportCache}
end

mutable struct AnalyzerState
    ### native ###

    native::NativeInterpreter

    ### JET.jl specific ###

    ## general ##

    # additional configurations for abstract interpretation
    analysis_params::JETAnalysisParams

    # key for the JET's global cache
    cache_key::UInt

    # reports found so far
    reports::Vector{InferenceErrorReport}

    # stashes `UncaughtExceptionReport`s that are not caught so far
    uncaught_exceptions::Vector{UncaughtExceptionReport}

    # keeps locations where `SeriousExceptionReport` are reported, in order to avoid
    # duplicated `UncaughtExceptionReport`s on same `throw` calls
    throw_locs::Vector{LineInfoNode}

    # keeps reports that should be updated when returning back the parent frame (i.e. the next time we get back to inter-procedural context)
    to_be_updated::Set{InferenceErrorReport}

    # keeps track of the current inference frame (needed for report cache reconstruction)
    current_frame::Union{Nothing,InferenceState}

    # whether to enable global report cache restoration at the next call of `CC.get(wvc::WorldView{JETGlobalCache}, mi::MethodInstance, default)`
    cache_enabled::Bool

    # local report cache for constant analysis
    cache::Vector{AnalysisResult}

    ## virtual toplevel execution ##

    # will be used in toplevel analysis (skip inference on actually interpreted statements)
    concretized::BitVector

    # virtual toplevel module
    toplevelmod::Module

    # slots to represent toplevel global variables
    global_slots::Dict{Int,Symbol}

    ## debug ##

    logger::JETLogger

    # records depth of call stack
    depth::Int
end

"""
    mutable struct AnalyzerState
        ...
    end

The mutable object that holds various states that are consumed by all [`AbstractAnalyzer`](@ref)s.

---

    AnalyzerState(analyzer::AbstractAnalyzer) -> AnalyzerState

If `NewAnalyzer` implements the `AbstractAnalyzer` interface, `NewAnalyzer` should implement
this `AnalyzerState(analyzer::NewAnalyzer) -> AnalyzerState` interface.

A new `AnalyzerState` is supposed to be constructed using [JET configurations](@ref) passed
as keyword arguments `jetconfigs` of the [`NewAnalyzer(; jetconfigs...)`](@ref AbstractAnalyzer)
constructor, and the constructed `AnalyzerState` is usually kept within `NewAnalyzer` itself:
```julia
function NewAnalyzer(; jetconfigs...)
    ...
    state = AnalyzerState(; jetconfigs...)
    return NewAnalyzer(..., state)
end
AnalyzerState(analyzer::NewAnalyzer) = analyzer.state
```
"""
function AnalyzerState(analyzer::Analyzer) where Analyzer<:AbstractAnalyzer
    error("""
    missing `$AbstractAnalyzer` API:
    `$Analyzer` is required to implement the `$AnalyzerState(analyzer::$Analyzer) -> $AnalyzerState` interface.
    See the documentation of `$AbstractAnalyzer` and `$AnalyzerState`.
    """)
end

function AbstractAnalyzer(analyzer::Analyzer, state::AnalyzerState) where Analyzer<:AbstractAnalyzer
    error("""
    missing `$AbstractAnalyzer` API:
    `$Analyzer` is required to implement the `$AbstractAnalyzer(analyzer::$Analyzer, state::$AnalyzerState) -> $Analyzer` interface.
    See the documentation of `$AbstractAnalyzer`.
    """)
end

# define getter/setter methods
for fld in fieldnames(AnalyzerState)
    fld === :cache_key && continue # will be defined later (in order to take in `ReportPass` hash)
    getter = Symbol("get_", fld)
    setter = Symbol("set_", fld, '!')
    @eval (@__MODULE__) @inline $getter(analyzer::AbstractAnalyzer)    = getfield(AnalyzerState(analyzer), $(QuoteNode(fld)))
    @eval (@__MODULE__) @inline $setter(analyzer::AbstractAnalyzer, v) = setfield!(AnalyzerState(analyzer), $(QuoteNode(fld)), v)
end

# constructor for fresh analysis
@jetconfigurable function AnalyzerState(world::UInt     = get_world_counter();
                                        analysis_params = nothing,
                                        current_frame   = nothing,
                                        cache           = AnalysisResult[],
                                        inf_params      = nothing,
                                        opt_params      = nothing,
                                        concretized     = _CONCRETIZED,
                                        toplevelmod     = _TOPLEVELMOD,
                                        global_slots    = _GLOBAL_SLOTS,
                                        logger          = nothing,
                                        depth           = 0,
                                        jetconfigs...)
    isnothing(analysis_params) && (analysis_params = JETAnalysisParams(; jetconfigs...))
    isnothing(inf_params)      && (inf_params      = JETInferenceParams(; jetconfigs...))
    isnothing(opt_params)      && (opt_params      = JETOptimizationParams(; jetconfigs...))
    isnothing(logger)          && (logger          = JETLogger(; jetconfigs...))

    return AnalyzerState(NativeInterpreter(world; inf_params, opt_params),
                         analysis_params,
                         get_cache_key_per_config(analysis_params, inf_params),
                         InferenceErrorReport[],
                         UncaughtExceptionReport[],
                         LineInfoNode[],
                         Set{InferenceErrorReport}(),
                         current_frame,
                         false,
                         cache,
                         concretized,
                         toplevelmod,
                         global_slots,
                         logger,
                         depth,
                         )
end

# dummies for non-toplevel analysis
module __toplevelmod__ end
const _CONCRETIZED  = BitVector()
const _TOPLEVELMOD  = __toplevelmod__
const _GLOBAL_SLOTS = Dict{Int,Symbol}()

# constructor for sequential toplevel JET analysis
function AbstractAnalyzer(analyzer::T, concretized, toplevelmod) where {T<:AbstractAnalyzer}
    newstate = AnalyzerState(# update world age to take in newly added methods defined
                             # in a previous toplevel interpretation performed by `ConcreteInterpreter`
                             get_world_counter();
                             analysis_params = JETAnalysisParams(analyzer),
                             inf_params      = InferenceParams(analyzer),
                             opt_params      = OptimizationParams(analyzer),
                             concretized     = concretized, # or construct partial `CodeInfo` from remaining abstract statements ?
                             toplevelmod     = toplevelmod,
                             logger          = JETLogger(analyzer),
                             )
    newanalyzer = AbstractAnalyzer(analyzer, newstate)
    maybe_initialize_caches!(newanalyzer)
    return newanalyzer
end

# constructor to do additional JET analysis in the middle of parent (non-toplevel) interpretation
function AbstractAnalyzer(analyzer::T) where {T<:AbstractAnalyzer}
    newstate = AnalyzerState(get_world_counter(analyzer);
                             analysis_params = JETAnalysisParams(analyzer),
                             current_frame   = get_current_frame(analyzer),
                             cache           = get_cache(analyzer),
                             inf_params      = InferenceParams(analyzer),
                             opt_params      = OptimizationParams(analyzer),
                             logger          = JETLogger(analyzer),
                             depth           = get_depth(analyzer),
                             )
    newanalyzer = AbstractAnalyzer(analyzer, newstate)
    maybe_initialize_caches!(newanalyzer)
    return newanalyzer
end

function Base.show(io::IO, analyzer::AbstractAnalyzer)
    rn = length(get_reports(analyzer))
    en = length(get_uncaught_exceptions(analyzer))
    print(io, "AbstractAnalyzer with $(rn) reports and $(en) uncaught exceptions")
    frame = get_current_frame(analyzer)
    if !isnothing(frame)
        print(io, " at ")
        show(io, frame.linfo)
    end
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", analyzer::AbstractAnalyzer) =
    return analyzer
Base.show(io::IO, ::MIME"application/prs.juno.inline", analyzer::NativeInterpreter) =
    return analyzer

# XXX: this should be upstreamed
function Base.show(io::IO, frame::InferenceState)
    print(io, "InfernceState for ")
    show(io, frame.linfo)
    print(io, " at pc ", frame.currpc, '/', length(frame.src.code))
end
Base.show(io::IO, ::MIME"application/prs.juno.inline", frame::InferenceState) =
    return frame

# AbstractInterpreter API
# -----------------------

CC.InferenceParams(analyzer::AbstractAnalyzer)    = InferenceParams(get_native(analyzer))
CC.OptimizationParams(analyzer::AbstractAnalyzer) = OptimizationParams(get_native(analyzer))
CC.get_world_counter(analyzer::AbstractAnalyzer)  = get_world_counter(get_native(analyzer))

# JET only works for runtime inference
CC.lock_mi_inference(::AbstractAnalyzer, ::MethodInstance) = nothing
CC.unlock_mi_inference(::AbstractAnalyzer, ::MethodInstance) = nothing

CC.add_remark!(analyzer::AbstractAnalyzer, sv, s) = ReportPass(analyzer)(NativeRemark, analyzer, sv, s)

CC.may_optimize(analyzer::AbstractAnalyzer)      = true
CC.may_compress(analyzer::AbstractAnalyzer)      = false
CC.may_discard_trees(analyzer::AbstractAnalyzer) = false

# AbstractAnalyzer
# ----------------

JETAnalysisParams(analyzer::AbstractAnalyzer) = get_analysis_params(analyzer)

# maybe we want to strip off `@nospecialize`s below ?

@inline add_new_report!(report::InferenceErrorReport, analyzer::AbstractAnalyzer) =
    return stash_report!(report, analyzer)
@inline restore_cached_report!(cache::InferenceErrorReportCache, analyzer::AbstractAnalyzer) =
    return stash_report!(restore_cached_report(cache), analyzer)
@inline function stash_report!(report::InferenceErrorReport, analyzer::AbstractAnalyzer)
    if isa(report, UncaughtExceptionReport)
        push!(get_uncaught_exceptions(analyzer), report)
    else
        push!(get_reports(analyzer), report)
    end
    return report
end

function get_cache_key_per_config(analysis_params::JETAnalysisParams, inf_params::InferenceParams)
    h = @static UInt === UInt64 ? 0xa49bd446c0a5d90e : 0xe45361ac
    h = hash(analysis_params, h)
    h = hash(inf_params, h)
    return h
end

function get_cache_key(analyzer::AbstractAnalyzer)
    h = AnalyzerState(analyzer).cache_key
    h = hash(ReportPass(analyzer), h)
    return h
end

function maybe_initialize_caches!(analyzer::AbstractAnalyzer)
    cache_key = get_cache_key(analyzer)
    haskey(JET_REPORT_CACHE, cache_key) || (JET_REPORT_CACHE[cache_key] = IdDict())
    haskey(JET_CODE_CACHE, cache_key)   || (JET_CODE_CACHE[cache_key]   = IdDict())
end

# check if we're in a toplevel module
@inline istoplevel(analyzer::AbstractAnalyzer, sv::InferenceState)    = istoplevel(analyzer, sv.linfo)
@inline istoplevel(::AbstractAnalyzer, ::OptimizationState)           = false # optimization never happen for top-level code
@inline istoplevel(analyzer::AbstractAnalyzer, linfo::MethodInstance) = get_toplevelmod(analyzer) === linfo.def

is_global_slot(analyzer::AbstractAnalyzer, slot::Int)   = slot in keys(get_global_slots(analyzer))
is_global_slot(analyzer::AbstractAnalyzer, slot::Slot)  = is_global_slot(analyzer, slot_id(slot))
is_global_slot(analyzer::AbstractAnalyzer, sym::Symbol) = sym in values(get_global_slots(analyzer))

JETLogger(analyzer::AbstractAnalyzer) = get_logger(analyzer)

@inline with_toplevel_logger(@nospecialize(f), analyzer::AbstractAnalyzer, @nospecialize(filter = ≥(DEFAULT_LOGGER_LEVEL)); kwargs...) =
    with_logger(f, JETLogger(analyzer).toplevel_logger, filter, "toplevel"; kwargs...)

@inline with_inference_logger(@nospecialize(f), analyzer::AbstractAnalyzer, @nospecialize(filter = ≥(DEFAULT_LOGGER_LEVEL)); kwargs...) =
    with_logger(f, JETLogger(analyzer).inference_logger, filter, "inference"; kwargs...)

@inline function with_logger(
    @nospecialize(f), io::Union{Nothing,IO}, @nospecialize(filter), logger_name;
    @nospecialize(pre = identity))
    isnothing(io) && return
    level = get_logger_level(io)
    filter(level) || return
    pre(io)
    print(io, "[$logger_name-$(LOGGER_LEVELS[level])] ")
    f(io)
end

"""
    NativeRemark <: InferenceErrorReport

This special `InferenceErrorReport` wraps remarks by `NativeInterpreter`.
"remarks" are information that Julia's native compiler emits about how its type inference goes,
and those remarks are less interesting in term of "error checking", so currently any of JET's
pre-defined report passes doesn't make any use of `NativeRemark`.
"""
@reportdef struct NativeRemark <: InferenceErrorReport
    s::String
end
get_msg(::Type{NativeRemark}, analyzer::AbstractAnalyzer, sv, s) = s

# ignores `NativeRemark`s by default
(::SoundBasicPass)(::Type{NativeRemark}, analyzer::AbstractAnalyzer, sv::InferenceState, s) = return

# JETAnalyzer
# -----------

"""
    struct JETAnalyzer <: AbstractAnalyzer

JET.jl's default error analyzer.
"""
struct JETAnalyzer{RP<:ReportPass} <: AbstractAnalyzer
    report_pass::RP
    state::AnalyzerState
end

# AbstractAnalyzer API requirements

# NOTE `@aggressive_constprop` here makes sure `mode` to be propagated as constant
@aggressive_constprop @jetconfigurable function JETAnalyzer(;
    report_pass::Union{Nothing,T} = nothing,
    mode::Symbol                  = :basic,
    # default `InferenceParams` tuning
    aggressive_constant_propagation::Bool = true,
    unoptimize_throw_blocks::Bool         = false,
    # default `OptimizationParams` tuning
    inlining::Bool = false,
    jetconfigs...) where {T<:ReportPass}
    if isnothing(report_pass)
        # if `report_pass` isn't passed explicitly, here we configure it according to `mode`
        report_pass = mode === :basic ? BasicPass() :
                      mode === :sound ? SoundPass() :
                      throw(ArgumentError("`mode` configuration should be either of `:basic` or `:sound`"))
    end
    # NOTE we always disable inlining, because:
    # - our current strategy to find undefined local variables and uncaught `throw` calls assumes un-inlined frames
    # - the cost for inlining isn't necessary for JETAnalyzer
    @assert !inlining "inlining should be disabled for JETAnalyzer"
    state = AnalyzerState(; aggressive_constant_propagation,
                            unoptimize_throw_blocks,
                            inlining,
                            jetconfigs...)
    return JETAnalyzer(report_pass, state)
end
AnalyzerState(analyzer::JETAnalyzer) =
    return analyzer.state
AbstractAnalyzer(analyzer::JETAnalyzer, state::AnalyzerState) =
    return JETAnalyzer(ReportPass(analyzer), state)
ReportPass(analyzer::JETAnalyzer) =
    return analyzer.report_pass

ext: jl


text:
"
https://arxiv.org/pdf/cond-mat/0203500.pdf
"
function dynamicaldmrg(A::Union{MPSComoving,FiniteMPS},z,ham::MPOHamiltonian;init=copy(A),solvtol=Defaults.tol,tol=solvtol*length(A)*2,maxiter=Defaults.maxiter,verbose=Defaults.verbose)
    w=real(z);eta=imag(z)

    envs1 = environments(init,ham) #environments for h
    (ham2,envs2) = squaredenvs(init,ham,envs1) #environments for h^2
    mixedenvs = environments(init,A); #environments for <init | A>

    delta=2*tol

    numit = 0
    while delta>tol && numit<maxiter
        numit+=1
        delta=0

        for i in [1:(length(A)-1);length(A):-1:2]

            #the alternative is using gradient descent, which is at least sure to converge...
            @tensor tos[-1 -2;-3]:=leftenv(mixedenvs,i,init)[-1,1]*A.AC[i][1,-2,2]*rightenv(mixedenvs,i,init)[2,-3]

            (res,convhist) = @closure linsolve(-eta*tos,init.AC[i],GMRES(tol=solvtol)) do x
                y=(eta*eta+w*w)*x
                y-=2*w*ac_prime(x,i,init,envs1)
                y+=ac_prime(x,i,init,envs2)
            end

            delta = max(delta,norm(res-init.AC[i]))
            init.AC[i] = res

            convhist.converged == 0 && @info "($(i)) failed to converge $(convhist.normres)"
        end

        verbose && @info "ddmrg sweep delta : $(delta)"
    end

    a = @tensor leftenv(mixedenvs,1,init)[-1,1]*A.AC[1][1,-2,2]*rightenv(mixedenvs,1,init)[2,-3]*conj(init.AC[1][-1,-2,-3])
    a = a';

    cb = leftenv(envs1,1,A);
    for i in 1:length(A)
        cb = transfer_left(cb,ham[i],init.AL[i],A.AL[i]);
    end

    b = 0*a
    for i in 1:length(cb)
        b+=@tensor cb[i][1,2,3]*A.CR[end][3,4]*rightenv(envs1,length(A),A)[i][4,2,5]*conj(init.CR[end][1,5]);
    end

    v = b/eta-w/eta*a+1im*a
    return v,init
end

function squaredenvs(state::Union{MPSComoving,FiniteMPS},ham::MPOHamiltonian,envs=environments(state,ham))
    nham=conj(ham,transpo=true)*ham

    #to construct the squared caches we will first initialize environments
    #then make all data invalid so it will be recalculated
    #then initialize the right caches at the edge
    ncocache=environments(state,nham)

    #make sure the dependencies are incorrect, so data will be recalculated
    for i in 1:length(state)
        poison!(ncocache,i)
    end

    #impose the correct boundary conditions (important for comoving mps, should do nothing for finite mps)
    indmap = LinearIndices((ham.odim,ham.odim));

    nleft=leftenv(ncocache,1,state)
    nright=rightenv(ncocache,length(state),state)

    for i in 1:ham.odim
        for j in 1:ham.odim
            @tensor temp[-1 -2 -3;-4]:=leftenv(envs,1,state)[j][1,-3,-4]*conj(leftenv(envs,1,state)[i][1,-2,-1])
            copy!(nleft[indmap[i,j]].data,temp.data)

            @tensor temp[-1 -2 -3;-4]:=rightenv(envs,length(state),state)[j][-1,-2,1]*conj(rightenv(envs,length(state),state)[i][-4,-3,1])
            copy!(nright[indmap[i,j]].data,temp.data)
        end
    end

    return nham,ncocache
end


text:
function generate_dfpair(n, str_cols = 2, float_cols = 5)
    DfPair([
        DataFrame(union(
            Dict("scol_$i" => [randstring(rand(4:20)) for _ = 1:n] for i = 1:str_cols),
            Dict("fcol_$i" => rand(n) for i = 1:float_cols),
        )) for i = 1:2
    ]...)
end

function load_result(dir::String)
    res_dir = joinpath(dir, "result")
    Dict(split(resf, ".")[1] => ParquetFiles.load(joinpath(res_dir, resf)) |> DataFrame |> df -> hcat(df.i1, df.i2)
    for resf in readdir(res_dir))
end


text:
# # Sugarscape
# **Growing Artificial Societies**

# ![](sugar.gif)

# (Descriptions below are from [this page](http://jasss.soc.surrey.ac.uk/12/1/6/appendixB/EpsteinAxtell1996.html))

# ---

# "Growing Artificial Societies" (Epstein & Axtell 1996) is a reference book for scientists interested in agent-based modelling and computer simulation. It represents one of the most paradigmatic and fascinating examples of the so-called generative approach to social science (Epstein 1999). In their book, Epstein & Axtell (1996) present a computational model where a heterogeneous population of autonomous agents compete for renewable resources that are unequally distributed over a 2-dimensional environment. Agents in the model are autonomous in that they are not governed by any central authority and they are heterogeneous in that they differ in their genetic attributes and their initial environmental endowments (e.g. their initial location and wealth). The model grows in complexity through the different chapters of the book as the agents are given the ability to engage in new activities such as sex, cultural exchange, trade, combat, disease transmission, etc. The core of Sugarscape has provided the basis for various extensions to study e.g. norm formation through cultural diffusion (Flentge et al. 2001) and the emergence of communication and cooperation in artificial societies (Buzing et al. 2005). Here we analyse the model described in the second chapter of Epstein & Axtell's (1996) book within the Markov chain framework.

# ## Model structure

# The first model that Epstein & Axtell (1996) present comprises a finite population of agents who live in an environment. The environment is represented by a two-dimensional grid which contains sugar in some of its cells, hence the name Sugarscape. Agents' role in this first model consists in wandering around the Sugarscape harvesting the greatest amount of sugar they can find.

# ### Environment

# The environment is a 50×50 grid that wraps around forming a torus. Grid positions have both a sugar level and a sugar capacity c. A cell's sugar level is the number of units of sugar in the cell (potentially none), and its sugar capacity c is the maximum value the sugar level can take on that cell. Sugar capacity is fixed for each individual cell and may be different for different cells. The spatial distribution of sugar capacities depicts a sugar topography consisting of two peaks (with sugar capacity c = 4) separated by a valley, and surrounded by a desert region of sugarless cells (see Figure 1) - note, however, that the grid wraps around in both directions–.

# The Sugarscape obbeys the following rule:

# Sugarscape growback rule G$\alpha$:
#     At each position, sugar grows back at a rate of $\alpha$ units per time-step up to the cell's capacity c.

# ### Agents

# Every agent is endowed with individual (life-long) characteristics that condition her skills and capacities to survive in the Sugarscape. These individual attributes are:

# * A vision _v_, which is the maximum number of positions the agent can see in each of the four principal lattice directions: north, south, east and west.
# * A metabolic rate _m_, which represents the units of sugar the agent burns per time-step.
# * A maximum age _max-age_, which is the maximum number of time-steps the agent can live.

# Agents also have the capacity to accumulate sugar wealth _w_. An agent's sugar wealth is incremented at the end of each time-step by the sugar collected and decremented by the agent's metabolic rate. __Two agents are not allowed to occupy the same position in the grid.__

# The agents' behaviour is determined by the following two rules:

# #### Agent movement rule _M_:

# Consider the set of unoccupied positions within your vision (including the one you are standing on), identify the one(s) with the greatest amount of sugar, select the nearest one (randomly if there is more than one), move there and collect all the sugar in it. At this point, the agent's accumulated sugar wealth is incremented by the sugar collected and decremented by the agent's metabolic rate _m_. If at this moment the agent's sugar wealth is not greater than zero, then the agent dies.

# #### Agent replacement rule _R_:

# Whenever an agent dies it is replaced by a new agent of age 0 placed on a randomly chosen unoccupied position, having random attributes _v_, _m_ and _max-age_, and random initial wealth w0. All random numbers are drawn from uniform distributions with ranges specified in Table 1 below.

# ### Scheduling of events

# Scheduling is determined by the order in which the different rules _G_, _M_ and _R_ are fired in the model. Environmental rule _G_ comes first, followed by agent rule _M_ (which is executed by all agents in random order) and finally agent rule _R_ is executed (again, by all agents in random order).

# ### Parameterisation

# Our analysis corresponds to a model used by Epstein & Axtell (1996, pg. 33) to study the emergent wealth distribution in the agent population. This model is parameterised as indicated in Table 1 below (where U[a,b] denotes a uniform distribution with range [a,b]).

# Initially, each position of the Sugarscape contains a sugar level equal to its sugar capacity c, and the 250 agents are created at a random unoccupied initial location and with random attributes (using the uniform distributions indicated in Table 1).

# __Table 1__

# | Parameter                                | Value        |
# |------------------------------------------|--------------|
# | Lattice length L                         | 50           |
# | Number of sugar peaks                    | 2            |
# | Growth rate $\alpha$                     | 1            |
# | Number of agents N                       | 250          |
# | Agents' initial wealth w0 distribution   | U[5,25]      |
# | Agents' metabolic rate m distribution    | U[1,4]       |
# | Agents' vision v distribution            | U[1,6]       |
# | Agents' maximum age max-age distribution | U[60,100]    |

using Agents
using AgentsPlots
using Random

mutable struct SugarSeeker <: AbstractAgent
    id::Int
    pos::Tuple{Int,Int}
    vision::Int
    metabolic_rate::Int
    age::Int
    max_age::Int
    wealth::Int
end

# Functions `distances` and `sugar_caps` produce a matrix for the distribution of sugar capacities."

function distances(pos, sugar_peaks, max_sugar)
    all_dists = Array{Int,1}(undef, length(sugar_peaks))
    for (ind, peak) in enumerate(sugar_peaks)
        d = round(Int, sqrt(sum((pos .- peak) .^ 2)))
        all_dists[ind] = d
    end
    return minimum(all_dists)
end

function sugar_caps(dims, sugar_peaks, max_sugar, dia = 4)
    sugar_capacities = zeros(Int, dims)
    for i in 1:dims[1], j in 1:dims[2]
        sugar_capacities[i, j] = distances((i, j), sugar_peaks, max_sugar)
    end
    for i in 1:dims[1]
        for j in 1:dims[2]
            sugar_capacities[i, j] = max(0, max_sugar - (sugar_capacities[i, j] ÷ dia))
        end
    end
    return sugar_capacities
end

"Start a sugarscape simulation"
function sugarscape(;
    dims = (50, 50),
    sugar_peaks = ((10, 40), (40, 10)),
    growth_rate = 1,
    N = 250,
    w0_dist = (5, 25),
    metabolic_rate_dist = (1, 4),
    vision_dist = (1, 6),
    max_age_dist = (60, 100),
    max_sugar = 4,
)
    sugar_capacities = sugar_caps(dims, sugar_peaks, max_sugar, 6)
    sugar_values = deepcopy(sugar_capacities)
    space = GridSpace(dims)
    properties = Dict(
        :growth_rate => growth_rate,
        :N => N,
        :w0_dist => w0_dist,
        :metabolic_rate_dist => metabolic_rate_dist,
        :vision_dist => vision_dist,
        :max_age_dist => max_age_dist,
        :sugar_values => sugar_values,
        :sugar_capacities => sugar_capacities,
    )
    model = AgentBasedModel(
        SugarSeeker,
        space,
        scheduler = random_activation,
        properties = properties,
    )
    for ag in 1:N
        add_agent_single!(
            model,
            rand(vision_dist[1]:vision_dist[2]),
            rand(metabolic_rate_dist[1]:metabolic_rate_dist[2]),
            0,
            rand(max_age_dist[1]:max_age_dist[2]),
            rand(w0_dist[1]:w0_dist[2]),
        )
    end
    return model
end

model = sugarscape()

# Fig. 1: Spatial distribution of sugar capacities in the Sugarscape. Cells are coloured according to their sugar capacity.

heatmap(model.sugar_capacities)

#

function env!(model)
    ## At each position, sugar grows back at a rate of $\alpha$ units per time-step up to the cell's capacity c.
    togrow = findall(
        x -> model.sugar_values[x] < model.sugar_capacities[x],
        1:length(positions(model)),
    )
    model.sugar_values[togrow] .+= model.growth_rate
end

function movement!(agent, model)
    newsite = agent.pos
    ## find all unoccupied position within vision
    neighbors = nearby_positions(agent.pos, model, agent.vision)
    empty = collect(empty_positions(model))
    if length(empty) > 0
        ## identify the one(s) with greatest amount of sugar
        available_sugar = (model.sugar_values[x,y] for (x, y) in empty)
        maxsugar = maximum(available_sugar)
        if maxsugar > 0
            sugary_sites_inds = findall(x -> x == maxsugar, collect(available_sugar))
            sugary_sites = empty[sugary_sites_inds]
            ## select the nearest one (randomly if more than one)
            for dia in 1:(agent.vision)
                np = nearby_positions(agent.pos, model, dia)
                suitable = intersect(np, sugary_sites)
                if length(suitable) > 0
                    newsite = rand(suitable)
                    break
                end
            end
            ## move there and collect all the sugar in it
            newsite != agent.pos && move_agent!(agent, newsite, model)
        end
    end
    ## update wealth (collected - consumed)
    agent.wealth += (model.sugar_values[newsite...] - agent.metabolic_rate)
    model.sugar_values[newsite...] = 0
    ## age
    agent.age += 1
end

function replacement!(agent, model)
    ## If the agent's sugar wealth become zero or less, it dies
    if agent.wealth <= 0 || agent.age >= agent.max_age
        kill_agent!(agent, model)
        ## Whenever an agent dies, a young one is added to a random pos.
        ## New agent has random attributes
        add_agent_single!(
            model,
            rand(model.vision_dist[1]:model.vision_dist[2]),
            rand(model.metabolic_rate_dist[1]:model.metabolic_rate_dist[2]),
            0,
            rand(model.max_age_dist[1]:model.max_age_dist[2]),
            rand(model.w0_dist[1]:model.w0_dist[2]),
        )
    end
end

function agent_step!(agent, model)
    movement!(agent, model)
    replacement!(agent, model)
end

# The following animation shows the emergent unequal distribution of agents on resourceful areas.

anim = @animate for i in 1:50
    step!(model, agent_step!, env!, 1)
    p1 = heatmap(model.sugar_values)
    p2 = plotabm(model, as = 3, am = :square, ac = :blue)
    title!(p1, "Sugar levels")
    title!(p2, "Agents\n Step $i")
    p = plot(p1, p2)
end
gif(anim, "sugar.gif", fps = 8)

# ### Distribution of wealth across individuals

model2 = sugarscape()
adata, _ = run!(model2, agent_step!, env!, 20, adata = [:wealth])

anim2 = @animate for i in 0:20
    histogram(
        adata[adata.step .== i, :wealth],
        legend = false,
        color = :black,
        nbins = 15,
        title = "step $i",
    )
end
nothing # hide

# We see that the distribution of wealth shifts from a more or less uniform distribution to a skewed distribution.

gif(anim2, fps = 3)

# ## References

# BUZING P, Eiben A & Schut M (2005) Emerging communication and cooperation in evolving agent societies. Journal of Artificial Societies and Social Simulation 8(1)2. http://jasss.soc.surrey.ac.uk/8/1/2.html.

# EPSTEIN J M (1999) Agent-Based Computational Models And Generative Social Science. Complexity 4(5), pp. 41-60.

# EPSTEIN J M & Axtell R L (1996) Growing Artificial Societies: Social Science from the Bottom Up. The MIT Press.

# FLENTGE F, Polani D & Uthmann T (2001) Modelling the emergence of possession norms using memes. Journal of Artificial Societies and Social Simulation 4(4)3. http://jasss.soc.surrey.ac.uk/4/4/3.html.


text:
concept_dist_different(x) = abs(x[1] - x[2]) ≠ abs(x[3] - x[4])

"""
    dist_different
Local constraint ensuring that `concept(dist_different, x) = |x[1] - x[2]| ≠ |x[3] - x[4]|)`.
"""
const dist_different = Constraint(
    args_length = 4,
    concept = concept_dist_different,
    error = make_error(:dist_different),
)


text:
using AdaptiveParticleMCMC
# 'Particle' and 'ParticleScratch' data types for SequentialMonteCarlo0
mutable struct MyParticle
    s::Float64
    MyParticle() = new(0.0)
end
mutable struct MyParam
    μ::Float64
    MyParam() = new(0.0)
end
# SequentialMonteCarlo M!, lM, lG:
M!(x, rng, k, x_prev, par) = (x.s = par.μ + randn(rng); nothing)
lM(k, x_prev, x, par) = -.5(x.s-par.μ)^2
lG(k, x, scratch) = -.5*x.s^2
# (Trivial) prior & parameter update function
test_prior(theta) = 1
set_param!(par, theta) = (par.μ = theta[1]; nothing)
# Test run: N particles, n iterations, time series length T
N=16; T=10; n=10000
# SequentialMonteCarlo data types
state = SMCState(T, N, MyParticle, MyParam, set_param!, lG, M!, lM)
# Run the algorithms
out_pmmh = adaptive_pmmh([0.0], test_prior, state, n)
out_pg = adaptive_pg([0.0], test_prior, state, n)
out_pg_aswam = adaptive_pg([0.0], test_prior, state, n; algorithm=:aswam)


text:
# Convert between different parameterizations.

module Transform

using Compat

# TODO: don't import Model; transformations should operate on
# generic ParamSets
using ..SensitiveFloats
import ..Log

import ..Model: NUM_SOURCE_TYPES, ParamSet, NUM_BANDS, NUM_COLOR_COMPONENTS

export DataTransform, ParamBounds, ParamBox, SimplexBox,
       get_mp_transform, enforce_bounds!,
       VariationalParams, FreeVariationalParams

# A vector of variational parameters.  The outer index is
# of celestial objects, and the inner index is over individual
# parameters for that object (referenced using ParamIndex).
const VariationalParams{T <: Number} = Vector{Vector{T}}
const FreeVariationalParams{T <: Number} = Vector{Vector{T}}

#####################################################################################
# this is essentially a compatibility layer since Model has gotten rid of this code
# it's messy, but this doesn't really matter too much, since Transforms.jl is going
# the way of the dinosaur

struct CanonicalParams <: ParamSet
    pos::Vector{Int}
    gal_frac_dev::Int
    gal_axis_ratio::Int
    gal_angle::Int
    gal_radius_px::Int
    flux_loc::Vector{Int}
    flux_scale::Vector{Int}
    color_mean::Matrix{Int}
    color_var::Matrix{Int}
    a::Matrix{Int}
    k::Matrix{Int}
    CanonicalParams() =
        new([1, 2], 3, 4, 5, 6,
            collect(7:(7+NUM_SOURCE_TYPES-1)),  # flux_loc
            collect((7+NUM_SOURCE_TYPES):(7+2NUM_SOURCE_TYPES-1)), # flux_scale
            reshape((7+2NUM_SOURCE_TYPES):(7+2NUM_SOURCE_TYPES+(NUM_BANDS-1)*NUM_SOURCE_TYPES-1), (NUM_BANDS-1, NUM_SOURCE_TYPES)),  # color_mean
            reshape((7+2NUM_SOURCE_TYPES+(NUM_BANDS-1)*NUM_SOURCE_TYPES):(7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES-1), (NUM_BANDS-1, NUM_SOURCE_TYPES)),  # color_var
            reshape((7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES):(7+3NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES-1), (NUM_SOURCE_TYPES, 1)),  # a
            reshape((7+3NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES):(7+3NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES+NUM_COLOR_COMPONENTS*NUM_SOURCE_TYPES-1), (NUM_COLOR_COMPONENTS, NUM_SOURCE_TYPES))) # k
end
const ids = CanonicalParams()
Base.length(::Type{CanonicalParams}) = 6 + 3*NUM_SOURCE_TYPES + 2*(NUM_BANDS-1)*NUM_SOURCE_TYPES + NUM_COLOR_COMPONENTS*NUM_SOURCE_TYPES

struct UnconstrainedParams <: ParamSet
    pos::Vector{Int}
    gal_frac_dev::Int
    gal_axis_ratio::Int
    gal_angle::Int
    gal_radius_px::Int
    flux_loc::Vector{Int}
    flux_scale::Vector{Int}
    color_mean::Matrix{Int}
    color_var::Matrix{Int}
    a::Matrix{Int}
    k::Matrix{Int}
    UnconstrainedParams() =
        new([1, 2], 3, 4, 5, 6,
            collect(7:(7+NUM_SOURCE_TYPES-1)),  # flux_loc
            collect((7+NUM_SOURCE_TYPES):(7+2NUM_SOURCE_TYPES-1)), # flux_scale
            reshape((7+2NUM_SOURCE_TYPES):(7+2NUM_SOURCE_TYPES+(NUM_BANDS-1)*NUM_SOURCE_TYPES-1), (NUM_BANDS-1, NUM_SOURCE_TYPES)),  # color_mean
            reshape((7+2NUM_SOURCE_TYPES+(NUM_BANDS-1)*NUM_SOURCE_TYPES):(7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES-1), (NUM_BANDS-1, NUM_SOURCE_TYPES)),  # color_var
            reshape((7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES):
                    (7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES+(NUM_SOURCE_TYPES-1)-1), (NUM_SOURCE_TYPES - 1, 1)),  # a
            reshape((7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES+(NUM_SOURCE_TYPES-1)):
                    (7+2NUM_SOURCE_TYPES+2*(NUM_BANDS-1)*NUM_SOURCE_TYPES+(NUM_SOURCE_TYPES-1)+(NUM_COLOR_COMPONENTS-1)*NUM_SOURCE_TYPES-1), (NUM_COLOR_COMPONENTS-1, NUM_SOURCE_TYPES))) # k
end
const ids_free = UnconstrainedParams()
Base.length(::Type{UnconstrainedParams}) =  6 + 2*NUM_SOURCE_TYPES + 2*(NUM_BANDS-1)*NUM_SOURCE_TYPES + (NUM_COLOR_COMPONENTS-1)*NUM_SOURCE_TYPES + NUM_SOURCE_TYPES-1

################################
# Elementary functions.

"""
Unconstrain x in the unit interval to lie in R.
"""
function inv_logit(x::Number)
    @assert(x >= 0)
    @assert(x <= 1)
    -log(1.0 / x - 1)
end

"""
Convert x in R to lie in the unit interval.
"""
function logit(x::Number)
    1.0 / (1.0 + exp(-x))
end

"""
Convert an (n - 1)-vector of real numbers to an n-vector on the simplex, where
the last entry implicitly has the untransformed value 1.
"""
function constrain_to_simplex(x::Vector{T}) where {T<:Number}
    m = maximum(x)
    if m == Inf
        # If more than 1 entry in x is Inf, it may be because the
        # the last entry in z is 0. Here we set all those entries to the
        # same value, though that may not be strictly correct.
        z   = T[ x_entry .== Inf ? one(T) : zero(T) for x_entry in x]
        z ./= sum(z)
        push!(z, 0)
        return z
    else
        z = exp.(x .- m)
        z_sum = sum(z) + exp(-m)
        z ./= z_sum
        push!(z, inv(z_sum)*exp(-m))
        return z
    end
end

"""
Convert an n-vector on the simplex to an (n - 1)-vector in R^{n -1}.  Entries
are expressed relative to the last element.
"""
function unconstrain_simplex(z::Vector{T}) where {T<:Number}
    n = length(z)
    T[ log(z[i]) - log(z[n]) for i = 1:(n - 1)]
end

################################
# The transforms for Celeste.

struct ParamBox
    lb::Float64  # lower bound
    ub::Float64  # upper bound
    scale::Float64

    function ParamBox(lb, ub, scale)
        @assert lb > -Inf # Not supported
        @assert scale > 0.0
        @assert lb < ub
        new(lb, ub, scale)
    end
end

struct SimplexBox
    lb::Float64  # lower bound
    scale::Float64
    n::Int

    function SimplexBox(lb, scale, n)
        @assert n >= 2
        @assert 0.0 <= lb < 1 / n
        new(lb, scale, n)
    end
end

# The vector of transform parameters for a Symbol.
const ParamBounds = Dict{Symbol, Union{Vector{ParamBox}, Vector{SimplexBox}}}

###############################################
# Functions for a "free transform".

function unbox_parameter(param::Number, pb::ParamBox)
    positive_constraint = (pb.ub == Inf)

    # exp and the logit functions handle infinities correctly, so
    # parameters can equal the bounds.
    @assert(pb.lb .<= param .<= pb.ub,
                    string("unbox_parameter: param outside bounds: ",
                                 "$param ($(pb.lb), $(pb.ub))"))

    if positive_constraint
        return log(param - pb.lb) * pb.scale
    else
        param_bounded = (param - pb.lb) / (pb.ub - pb.lb)
        return inv_logit(param_bounded) * pb.scale
    end
end

function box_parameter(free_param::Number, pb::ParamBox)
    positive_constraint = (pb.ub == Inf)
    if positive_constraint
        return exp(free_param / pb.scale) + pb.lb
    else
        return logit(free_param / pb.scale) * (pb.ub - pb.lb) + pb.lb
    end
end

"""
Convert an unconstrained (n-1)-vector to a simplicial n-vector, z, such that
  - sum(z) = 1
  - z >= sb.lb
See notes for a derivation and reasoning.
"""
function simplexify_parameter(free_param::Vector{T},
                              sb::SimplexBox) where {T<:Number}
    @assert length(free_param) == (sb.n - 1)

    # Broadcasting doesn't work with DualNumbers and Floats. :(
    # z_sim is on an unconstrained simplex.
    z_sim = constrain_to_simplex(T[ p / sb.scale for p in free_param ])
    param = T[ (1 - sb.n * sb.lb) * p + sb.lb for p in z_sim ]

    param
end

"""
Invert the transformation simplexify_parameter() by converting an n-vector
on a simplex to R^{n - 1}.
"""
function unsimplexify_parameter(param::Vector{T},
                                sb::SimplexBox) where {T<:Number}
    @assert length(param) == sb.n
    @assert all(param .>= sb.lb)
    @assert(abs(sum(param) - 1) < 1e-14, abs(sum(param) - 1))

    # z_sim is on an unconstrained simplex.
    # Broadcasting doesn't work with DualNumbers and Floats. :(
    z_sim = T[ (p - sb.lb) / (1 - sb.n * sb.lb) for p in param ]
    free_param = T[ p * sb.scale for p in unconstrain_simplex(z_sim) ]

    free_param
end

##################
# Derivatives

"""
Return derivatives of an unscaled transform from free parameters to a simplex.

Args:
  - z_sim: A vector in a simplex (NB: the function returns derivatives of the
           function f(unconstrained) = simplex, but it the answer is expressed
           in terms of the output of the function.)

Returns:
  - jacobian: n by (n -1) matrix of derivatives of the simplex output (in rows)
              wrt the free parameters (in columns)
  - hessian_vec: An n-length vector of the hessian of each simplex output
                 parameter with respect to the (n-1) free input parameters.
"""
function simplex_derivatives(z_sim::Vector{T}) where {T<:Number}
    n = length(z_sim)
    hessian_vec = Vector{Array{T}}(n)

    for i = 1:n
        hessian_vec[i] = Matrix{T}(n - 1, n - 1)
        for j=1:(n - 1), k=1:(n - 1)
            if j != k
                if (j == i)
                    hessian_vec[i][j, k] = -z_sim[i] * z_sim[k] * (1 - 2 * z_sim[i])
                elseif (k == i)
                    hessian_vec[i][j, k] = -z_sim[i] * z_sim[j] * (1 - 2 * z_sim[i])
                else
                    hessian_vec[i][j, k] = 2 * z_sim[i] * z_sim[j] * z_sim[k]
                end
            else # j == k
                if i == j # All equal
                    hessian_vec[i][j, k] = z_sim[i] * (1 - z_sim[j]) * (1 - 2 * z_sim[k])
                else # j == k, but both are different from i
                    hessian_vec[i][j, k] = - z_sim[i] * z_sim[j] * (1 - 2 * z_sim[j])
                end
            end
        end
    end

    jacobian =
        T[ z_sim[i] * (i == j) - z_sim[i] * z_sim[j] for i=1:n, j=1:(n - 1) ]

    jacobian, hessian_vec
end

"""
Return the derivative and hessian of a simplex transform given the constrained
parameters.

Args:
  - param: The constrained parameter (NB: the derivatives are expressed
               as a function of the constrained parameterd despite being
                     the derivative of the function unconstrained -> constrained)
    - sb: A box simplex constraint
"""
function box_simplex_derivatives(param::Vector{T},
                                 sb::SimplexBox) where {T<:Number}

    @assert length(param) == sb.n

    # z_sim is on an unconstrained simplex.
    # Broadcasting doesn't work with DualNumbers and Floats. :(
    z_sim = T[ (p - sb.lb) / (1 - sb.n * sb.lb) for p in param ]

    jacobian, hessian_vec = simplex_derivatives(z_sim)
    for i in 1:sb.n
        hessian_vec[i] *= (sb.scale ^ 2) * (1 - sb.n * sb.lb)
    end
    jacobian *= sb.scale * (1 - sb.n * sb.lb)

    jacobian, hessian_vec
end

"""
Return the derivative and hessian of a box transform given the constrained
parameters.

Args:
  - param: The constrained parameter (NB: the derivatives are expressed
               as a function of the constrained parameterd despite being
                     the derivative of the function unconstrained -> constrained)
    - pb: A box constraint
"""
function box_derivatives(param::Number, pb::ParamBox)
    if pb.ub == Inf
        centered_param = param - pb.lb
        return pb.scale * centered_param, pb.scale ^ 2 * centered_param
    else
        param_range = pb.ub - pb.lb
        centered_param = (param - pb.lb) / param_range
        derivative = param_range * centered_param * (1 - centered_param) / pb.scale
        return derivative, derivative * (1 - 2 * centered_param) / pb.scale
    end
end

"""
A datatype containing derivatives of the transform from free to constrained
parameters.

Members:
  dparam_dfree: The Jacobian of the transformation
                contrainted_param = f(free param)
  d2param_dfree2: A vector of hessians.  Each element is the Hessian of one
                  component of the aforementioned f()
"""
struct TransformDerivatives{T<:Number}
    dparam_dfree::Matrix{T}
    d2param_dfree2::Vector{Matrix{T}}
    Sa::Int

    # TODO: use sparse matrices?
    function TransformDerivatives{T}(Sa::Int) where {T<:Number}
        dparam_dfree =
            zeros(T, Sa * length(CanonicalParams), Sa * length(UnconstrainedParams))
        d2param_dfree2 = Vector{Matrix{T}}(Sa * length(CanonicalParams))
        for i in 1:(Sa * length(CanonicalParams))
            d2param_dfree2[i] =
                zeros(T, Sa * length(UnconstrainedParams), Sa * length(UnconstrainedParams))
        end
        new(dparam_dfree, d2param_dfree2, Sa)
    end
end

"""
Populate a TransformDerivatives object in place.

Args:
  - bounds: A vector containing one ParamBounds for each active source in the
            same order as active_sources.
  - transform_derivatives: TransformDerivatives to be populated.

Returns:
  Update transform_derivatives in place.
"""
function get_transform_derivatives!(
        vp::VariationalParams,
        active_sources::Vector{Int},
        bounds::Vector{ParamBounds},
        transform_derivatives::TransformDerivatives)

    @assert transform_derivatives.Sa == length(active_sources)

    for param in fieldnames(ids), sa = 1:length(active_sources)
        s = active_sources[sa]
        constraint_vec = bounds[sa][param]

        if isa(constraint_vec[1], ParamBox) # It is a box constraint
            @assert(length(constraint_vec) ==
                    length(getfield(ids_free, param)) ==
                    length(getfield(ids, param)))

            # Get each components' derivatives one by one.
            for ind = 1:length(constraint_vec)
                @assert isa(constraint_vec[ind], ParamBox)
                vp_ind = getfield(ids, param)[ind]
                vp_free_ind = getfield(ids_free, param)[ind]

                jac, hess = box_derivatives(vp[s][vp_ind], constraint_vec[ind])

                vp_sf_ind = length(CanonicalParams) * (sa - 1) + vp_ind
                vp_free_sf_ind = length(UnconstrainedParams) * (sa - 1) + vp_free_ind

                transform_derivatives.dparam_dfree[vp_sf_ind, vp_free_sf_ind] = jac
                transform_derivatives.d2param_dfree2[vp_sf_ind][vp_free_sf_ind, vp_free_sf_ind] = hess
            end
        else # It is a simplex constraint

            # If a param is not a box constraint, it must have all simplex constraints.
            @assert all([ isa(constraint, SimplexBox)  for constraint in constraint_vec])

            param_size = size(getfield(ids, param))
            @assert length(constraint_vec) == param_size[2]
            for col=1:(param_size[2])
                vp_free_ind = getfield(ids_free, param)[:, col]
                vp_ind = getfield(ids, param)[:, col]
                vp_sf_ind = length(CanonicalParams) * (sa - 1) + vp_ind
                vp_free_sf_ind = length(UnconstrainedParams) * (sa - 1) + vp_free_ind

                jac, hess = Transform.box_simplex_derivatives(
                    vp[s][vp_ind], constraint_vec[col])

                transform_derivatives.dparam_dfree[vp_sf_ind, vp_free_sf_ind] = jac
                for row in 1:(param_size[1])
                    transform_derivatives.d2param_dfree2[
                         vp_sf_ind[row]][vp_free_sf_ind, vp_free_sf_ind] = hess[row]
                end
            end
        end
    end
end

function get_transform_derivatives(vp::VariationalParams,
                                   active_sources::Vector{Int},
                                   bounds::Vector{ParamBounds})
  transform_derivatives =
    TransformDerivatives{Float64}(length(active_sources))
  get_transform_derivatives!(vp, active_sources, bounds, transform_derivatives)
  transform_derivatives
end

######################
# Functions to take actual parameter vectors.

"""
Convert between variational parameter vectors and unconstrained parameters.

Args:
  - vp: The vector of contrainted variationl parameters
  - vp_free: The vector of uncontrainted parameters for optimization
  - bounds: ParamBounds describing the transformation
  - to_unconstrained: If true, converts vp to vp_free.  If false, converts
                      vp_free to vp.

Returns:
  If to_unconstrained is true, updates vp_free in place.
  If to_unconstrained is false, updates vp in place.
"""
function perform_transform!(vp::Vector{T},
                            vp_free::Vector{T},
                            bounds::ParamBounds,
                            to_unconstrained::Bool) where {T<:Number}
    for (param, constraint_vec) in bounds
        is_box = isa(bounds[param], Array{ParamBox})
        if is_box
            # Apply a box constraint to each parameter.
            @assert(length(getfield(ids, param)) == length(getfield(ids_free, param)) ==
                length(bounds[param]))
            for ind in 1:length(getfield(ids, param))
                constraint = constraint_vec[ind]
                free_ind = getfield(ids_free, param)[ind]
                vp_ind = getfield(ids, param)[ind]
                to_unconstrained ?
                    vp_free[free_ind] = unbox_parameter(vp[vp_ind], constraint):
                    vp[vp_ind] = box_parameter(vp_free[free_ind], constraint)
            end
        else
            # Apply a simplex constraint to each parameter.
            @assert isa(bounds[param], Array{SimplexBox})
            param_size = size(getfield(ids, param))

            # Each column is a simplex and should have its own
            # simplicial constraints.
            @assert length(bounds[param]) == param_size[2]
            for col in 1:(param_size[2])
                free_ind = getfield(ids_free, param)[:, col]
                vp_ind = getfield(ids, param)[:, col]
                constraint = constraint_vec[col]
                to_unconstrained ?
                    vp_free[free_ind] = unsimplexify_parameter(vp[vp_ind], constraint):
                    vp[vp_ind] = simplexify_parameter(vp_free[free_ind], constraint)
            end
        end
    end
end

function free_to_vp!(vp_free::Vector{T},
                     vp::Vector{T},
                     bounds::ParamBounds) where {T<:Number}
    perform_transform!(vp, vp_free, bounds, false)
end

function vp_to_free!(vp::Vector{T},
                     vp_free::Vector{T},
                     bounds::ParamBounds) where {T<:Number}
    perform_transform!(vp, vp_free, bounds, true)
end

#####################
# Conversion to and from variational parameter vectors and arrays.

"""
Transform VariationalParams to an array.

Args:
  - vp = variational parameters
  - omitted_ids = ids in ParamIndex

There is probably no use for this function, since you'll only be passing
trasformations to the optimizer, but I'll include it for completeness.
"""
function free_vp_to_array(vp::FreeVariationalParams{T},
                          omitted_ids::Vector{Int}) where {T<:Number}

    left_ids = setdiff(1:length(UnconstrainedParams), omitted_ids)
    new_P = length(left_ids)
    S = length(vp)
    x_new = zeros(T, new_P, S)

    for p1 in 1:length(left_ids), s=1:S
        p0 = left_ids[p1]
        x_new[p1, s] = vp[s][p0]
    end

    x_new
end

"""
Transform a parameter vector to variational parameters in place.

Args:
 - xs: A (param x sources) matrix created from free variational parameters.
 - vp_free: Free variational parameters.  Only the ids in kept_ids
            will be updated.
 - kept_ids: Ids to update (from ids_free)

Returns:
 - Update vp_free in place.
"""
function array_to_free_vp!(xs::Array{T},
                           vp_free::FreeVariationalParams{T},
                           kept_ids::Vector{Int}) where {T}
    for s in 1:length(vp_free), p1 in 1:length(kept_ids)
        vp_free[s][kept_ids[p1]] = xs[p1, s]
    end
end

#########################
# Define the exported variables.

"""
Functions to move between a single source's variational parameters and a
transformation of the data for optimization.

to_vp: A function that takes transformed parameters and returns
             variational parameters
from_vp: A function that takes variational parameters and returned
                 transformed parameters
to_vp!: A function that takes (transformed paramters, variational parameters)
                and updates the variational parameters in place
from_vp!: A function that takes (variational paramters, transformed parameters)
                    and updates the transformed parameters in place
...
transform_sensitive_float: A function that takes (sensitive float, model
    parameters) where the sensitive float contains partial derivatives with
    respect to the variational parameters and returns a sensitive float with total
    derivatives with respect to the transformed parameters.
bounds: The bounds for each parameter and each object in ElboArgs.
active_sources: The sources that are being optimized.    Only these sources'
    parameters are transformed into the parameter vector.
"""
struct DataTransform
    bounds::Vector{ParamBounds}
    active_sources::Vector{Int}
    S::Int
    # TODO: Maybe this should be initialized with ElboArgs with optional
    # custom bounds. Or maybe it should be part of ElboArgs with one transform
    # per celestial object rather than a single object containing an array of
    # transforms.
    function DataTransform(bounds::Vector{ParamBounds},
                           active_sources=collect(1:length(bounds)),
                           S=length(bounds))
        @assert length(bounds) == length(active_sources)
        @assert maximum(active_sources) <= S
        new(bounds, active_sources, S)
    end
end

function from_vp!(dt::DataTransform,
                  vp::VariationalParams{T},
                  vp_free::VariationalParams{T}) where {T}
    active_sources = dt.active_sources
    bounds = dt.bounds
    for i in eachindex(active_sources)
        vp_to_free!(vp[active_sources[i]], vp_free[i], bounds[i])
    end
    return nothing
end

function to_vp!(dt::DataTransform,
                vp_free::FreeVariationalParams{T},
                vp::VariationalParams{T}) where {T<:Number}
    active_sources = dt.active_sources
    bounds = dt.bounds
    for i in eachindex(active_sources)
        free_to_vp!(vp_free[i], vp[active_sources[i]], bounds[i])
    end
    return nothing
end

function from_vp(dt::DataTransform,
                 vp::VariationalParams{T}) where {T}
    vp_free = [zeros(T, length(ids_free)) for _ in 1:length(dt.active_sources)]
    from_vp!(dt, vp, vp_free)
    vp_free
end

function to_vp(dt::DataTransform, vp_free::FreeVariationalParams{T}) where {T}
    n_active_sources = length(dt.active_sources)
    @assert(n_active_sources == dt.S,
            string("to_vp is not supported when active_sources is a ",
                   "strict subset of all sources."))
    vp = [zeros(T, length(CanonicalParams)) for _ in 1:n_active_sources]
    to_vp!(dt, vp_free, vp)
    vp
end

function vp_to_array(dt::DataTransform,
                     vp::VariationalParams{T},
                     omitted_ids::Vector{Int}) where {T}
    vp_trans = from_vp(dt, vp)
    free_vp_to_array(vp_trans, omitted_ids)
end

function array_to_vp!(dt::DataTransform,
                      xs::Array{T},
                      vp::VariationalParams{T},
                      kept_ids::Vector{Int}) where {T}
    # This needs to update vp in place so that variables in omitted_ids
    # stay at their original values.
    vp_trans = from_vp(dt, vp)
    array_to_free_vp!(xs, vp_trans, kept_ids)
    to_vp!(dt, vp_trans, vp)
end

# Given a sensitive float with derivatives with respect to all the
# constrained parameters, calculate derivatives with respect to
# the unconstrained parameters.
#
# Note that all the other functions in DeterministicVI calculated derivatives with
# respect to the constrained parameterization.
function transform_sensitive_float!(dt::DataTransform,
                                    sf_free::SensitiveFloat,
                                    sf::SensitiveFloat,
                                    vp::VariationalParams{T},
                                    active_sources::Vector{Int}) where {T}

    if isnan(sf.v[])
       error("sf has NaN value:", sf_free.v[])
    end

    if any(isnan.(sf.d))
       error("sf has NaN derivatives:", sf_free.d[])
    end

    if any(isnan.(sf.h))
       error("sf has NaN hessian:", sf_free.h)
    end

    n_active_sources = length(active_sources)
    @assert size(sf.d) == (length(CanonicalParams), n_active_sources)

    transform_derivatives = get_transform_derivatives(vp, active_sources, dt.bounds)

    d, h = sf.d, sf.h
    free_d, free_h = sf_free.d, sf_free.h
    dparam_dfree = transform_derivatives.dparam_dfree
    d2param_dfree2 = transform_derivatives.d2param_dfree2

    sf_free.v[] = sf.v[]

    copy!(free_d, reshape(dparam_dfree' * d, length(UnconstrainedParams), n_active_sources))

    At_mul_B!(free_h, dparam_dfree, h * dparam_dfree)

    for i in eachindex(d)
        X, n = d2param_dfree2[i], d[i]
        for j in eachindex(X)
            free_h[j] += X[j] * n
        end
    end

    symmetrize!(free_h)

    if isnan(sf_free.v[])
        error("sf_free has NaN value:", sf_free.v[])
    end

    if any(isnan.(sf_free.d))
        error("sf_free has NaN derivatives:", sf_free.d[])
    end

    if any(isnan.(sf_free.h))
        error("sf_free has NaN hessian:", sf_free.h)
    end

    return sf_free
end

# Ensure exact symmetry, which is necessary for
# some numerical linear algebra routines.
function symmetrize!(A, c = 0.5)
    for i in 1:size(A, 1), j in 1:i
        A[i, j] = A[j, i] = c * (A[i, j] + A[j, i])
    end
    return A
end

function get_mp_transform(vp::VariationalParams{T},
                          active_sources::Vector{Int};
                          loc_scale=1.0,
                          loc_width=1e-4) where {T<:Number}
    bounds = Vector{ParamBounds}(length(active_sources))

    # Note that, for numerical reasons, the bounds must be on the scale
    # of reasonably meaningful changes.
    for si in 1:length(active_sources)
        s = active_sources[si]
        bounds[si] = ParamBounds()
        bounds[si][:pos] = Vector{ParamBox}(2)
        pos = vp[s][ids.pos]
        for axis in 1:2
            bounds[si][:pos][axis] =
                ParamBox(pos[axis] - loc_width, pos[axis] + loc_width, loc_scale)
        end
        bounds[si][:flux_loc] = Vector{ParamBox}(NUM_SOURCE_TYPES)
        bounds[si][:flux_scale] = Vector{ParamBox}(NUM_SOURCE_TYPES)
        for i in 1:NUM_SOURCE_TYPES
            bounds[si][:flux_loc][i] = ParamBox(-1.0, 10., 1.0)
            bounds[si][:flux_scale][i] = ParamBox(1e-4, 0.1, 1.0)
        end
        bounds[si][:color_mean] = Vector{ParamBox}(4 * NUM_SOURCE_TYPES)
        bounds[si][:color_var] = Vector{ParamBox}(4 * NUM_SOURCE_TYPES)
        for ind in 1:length(ids.color_mean)
            bounds[si][:color_mean][ind] = ParamBox(-10., 10., 1.0)
            bounds[si][:color_var][ind] = ParamBox(1e-4, 1., 1.0)
        end
        bounds[si][:gal_frac_dev] = ParamBox[ ParamBox(1e-2, 1 - 1e-2, 1.0) ]
        bounds[si][:gal_axis_ratio] = ParamBox[ ParamBox(1e-2, 1 - 1e-2, 1.0) ]
        bounds[si][:gal_angle] = ParamBox[ ParamBox(-10.0, 10.0, 1.0) ]
        bounds[si][:gal_radius_px] = ParamBox[ ParamBox(0.1, 70., 1.0) ]

        bounds[si][:a] = Vector{SimplexBox}(1)
        bounds[si][:a][1] = SimplexBox(0.005, 1.0, 2)

        bounds[si][:k] = Vector{SimplexBox}(NUM_SOURCE_TYPES)
        for i in 1:NUM_SOURCE_TYPES
            bounds[si][:k][i] = SimplexBox(0.01 / NUM_COLOR_COMPONENTS, 1.0, NUM_COLOR_COMPONENTS)
        end
    end

    DataTransform(bounds, active_sources, length(vp))
end

"""
Put the variational parameters within the bounds of the transform.
Returns:
  Updates vp in place.
"""
function enforce_bounds!(vp::VariationalParams{T},
                         active_sources::Vector{Int},
                         transform::DataTransform) where {T<:Number}
    for sa in eachindex(transform.active_sources), (param, constraint_vec) in transform.bounds[sa]
        s = active_sources[sa]
        is_box = isa(constraint_vec, Array{ParamBox})
        if is_box
            # Box parameters.
            for ind in 1:length(getfield(ids, param))
                constraint = constraint_vec[ind]
                if !(constraint.lb <=
                         vp[s][getfield(ids, param)[ind]] <=
                         constraint.ub)
                    # Don't set the value to exactly the lower bound to avoid Inf
                    diff = constraint.ub - constraint.lb
                    epsilon = diff == Inf ? 1e-12: diff * 1e-12
                    vp[s][getfield(ids, param)[ind]] =
                        min(vp[s][getfield(ids, param)[ind]], constraint.ub - epsilon)
                    vp[s][getfield(ids, param)[ind]] =
                        max(vp[s][getfield(ids, param)[ind]], constraint.lb + epsilon)
                end
            end
        else
            # Simplex parameters.
            param_size = size(getfield(ids, param))
            for col in 1:param_size[2]
                constraint = constraint_vec[col]
                param_sum = zero(T)
                for row in 1:param_size[1]
                    if !(constraint.lb <= vp[s][getfield(ids, param)[row, col]] <= 1.0)
                        # Don't set the value to exactly the lower bound to avoid Inf
                        epsilon = (1.0 - constraint.lb) * 1e-12
                        vp[s][getfield(ids, param)[row, col]] =
                            min(vp[s][getfield(ids, param)[row, col]], 1.0 - epsilon)
                        vp[s][getfield(ids, param)[row, col]] =
                            max(vp[s][getfield(ids, param)[row, col]],
                                    constraint.lb + epsilon)
                    end
                    param_sum += vp[s][getfield(ids, param)[row, col]]
                end
                if param_sum != 1.0
                    # Normalize in a way that maintains the lower bounds
                    rescale =
                      (1 - constraint.n * constraint.lb) /
                      (param_sum - constraint.n * constraint.lb)
                    vp[s][getfield(ids, param)[:, col]] =
                      constraint.lb +
                      rescale * (vp[s][getfield(ids, param)[:, col]] -
                                 constraint.lb)
                end
            end
        end
    end
end

end


text:
let
    pts = [rand(Point{2,Float64}) for i in 1:10]
    rect = GeometryBasics.Rect(Vec2(0.0,0.0),Vec2(1.0,1.0))
    controller = voronoi_controller(pts,rect)
    pts = pts .+ get_vels(controller)
    update_controller!(controller,pts)
end


text:
#### Types

abstract type Language end
struct Fortran end; export Fortran
struct Python end; export Python
struct Cython end; export Cython
struct Julia end; export Julia


text:
# CUDA's complex types are defined in terms of vector types (float2, double2),
# but those seem compatible with Julia's complex numbers, so use those.
const cuFloatComplex = Complex{Float32}
const cuDoubleComplex = Complex{Float64}

# aliases
const cuComplex = cuFloatComplex


text:
#####
##### Outer functions for filling halo regions for open boundary conditions.
#####

# TODO: support true open boundary conditions.
# For this we need to have separate functions for each of the six boundaries,
# and need to unroll a loop over the boundary normal direction.
# The syntax for `getbc` is also different for OpenBoundaryCondition than for others,
# because the boundary-normal index can vary (and array boundary conditions need to be
# 3D in general).

@kernel function set_east_west_u_velocity!(u, i_boundary, bc, grid, args...)
    j, k = @index(Global, NTuple)
    @inbounds u[i_boundary, j, k] = getbc(bc, j, k, grid, args...)
end

@kernel function set_north_south_v_velocity!(v, j_boundary, bc, grid, args...)
    i, k = @index(Global, NTuple)
    @inbounds v[i, j_boundary, k] = getbc(bc, i, k, grid, args...)
end

@kernel function set_top_bottom_w_velocity!(w, k_boundary, bc, grid, args...)
    i, j = @index(Global, NTuple)
    @inbounds w[i, j, k_boundary] = getbc(bc, i, j, grid, args...)
end

  fill_west_halo!(u, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :yz, set_east_west_u_velocity!,   u,           1, bc, grid, args...; dependencies=dep, kwargs...)
  fill_east_halo!(u, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :yz, set_east_west_u_velocity!,   u, grid.Nx + 1, bc, grid, args...; dependencies=dep, kwargs...)
 fill_south_halo!(v, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xz, set_north_south_v_velocity!, v,           1, bc, grid, args...; dependencies=dep, kwargs...)
 fill_north_halo!(v, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xz, set_north_south_v_velocity!, v, grid.Ny + 1, bc, grid, args...; dependencies=dep, kwargs...)
fill_bottom_halo!(w, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xy, set_top_bottom_w_velocity!,  w,           1, bc, grid, args...; dependencies=dep, kwargs...)
   fill_top_halo!(w, bc::OBC, arch, dep, grid, args...; kwargs...) = launch!(arch, grid, :xy, set_top_bottom_w_velocity!,  w, grid.Nz + 1, bc, grid, args...; dependencies=dep, kwargs...)


text:
#
# Static Hedging
#

function initial_state(ctx::CompilerContext{M}) where {M<:StaticHedgingModel}
    return get_pricing_date(ctx)
end

function get_riskfree_curve_symbol(m::StaticHedgingModel, kind::String, rf::CashRiskFactor) :: Symbol
    @assert haskey(m.riskfree_curve_map, kind) "Pricing Model doesn't know about a riskfree curve of type $kind. These are the available types of riskfree curves: $(collect(keys(m.riskfree_curve_map)))."
    riskfree_curve_map = m.riskfree_curve_map[kind]
    sym = risk_factor_symbol(rf)
    @assert haskey(riskfree_curve_map, sym) "Risk factor $sym not available in this pricing model. These are the available items: $(collect(keys(riskfree_curve_map)))."
    return riskfree_curve_map[sym]
end

function is_riskfree_discountfactor(model::StaticHedgingModel, kind::String, rf::DiscountFactor) :: Bool
    @assert haskey(model.riskfree_curve_map, kind) "Pricing Model doesn't know about a riskfree curve of type $kind. These are the available types of riskfree curves: $(collect(keys(model.riskfree_curve_map)))."
    riskfree_curve_map = model.riskfree_curve_map[kind]
    sym = market_data_symbol(rf)

    for curve_sym in values(riskfree_curve_map)
        if curve_sym == sym
            return true
        end
    end

    return false
end

# in StaticHedging, all risk factors are observed at the pricing_date
function lower!(ctx::CompilerContext{M}, rf::RiskFactor, at::Date) :: OptimizingIR.ImmutableValue where {M<:StaticHedgingModel}
    assert_at_initial_state(ctx, at)
    return add_input_riskfactor!(ctx, rf)
end

function lower!(ctx::CompilerContext{StaticHedgingModel{C}}, c::Unit{SpotCurrency{C}}, at::Date) :: OptimizingIR.ImmutableValue where {C<:Currencies.Currency}
    assert_at_initial_state(ctx, at)
    # same currency
    return lower!(ctx, Konst(1.0), at)
end

function lower!(ctx::CompilerContext{StaticHedgingModel{C1}}, c::Unit{SpotCurrency{C2}}, at::Date) :: OptimizingIR.ImmutableValue where {C1<:Currencies.Currency, C2<:Currencies.Currency}
    assert_at_initial_state(ctx, at)
    return lower!(ctx, c.o, at)
end

function lower!(ctx::CompilerContext{M}, c::Unit{Stock}, at::Date) where {M<:StaticHedgingModel}
    assert_at_initial_state(ctx, at)
    return lower!(ctx, c.o, at)
end

function lower!(ctx::CompilerContext{M}, o::ObservableAt, at::Date) :: OptimizingIR.ImmutableValue where {M<:StaticHedgingModel}
    return lower!(ctx, o.o, resolve_compile_time_value(ctx, o.at, at))
end

# Present Value for a Unit
function lower!(ctx::CompilerContext{M}, c::WhenAt{U}, at::Date) :: OptimizingIR.ImmutableValue where {U<:Unit, M<:StaticHedgingModel}
    assert_at_initial_state(ctx, at)
    pricing_date = get_pricing_date(ctx)

    unit = c.c
    underlying = unit.o
    mat = c.maturity

    if pricing_date > mat
        # ignore old cashflows
        return lower!(ctx, Worthless(), at)
    elseif pricing_date == mat
        # returns the underlying asset value if at maturity
        return lower!(ctx, c.c, at)
    end

    riskfree_discount_factor = DiscountFactor(get_riskfree_curve_symbol(ctx.model, ctx.attr["riskfree_curves"], underlying), mat)
    carry_type = get_carry_type(ctx.attr)

    if (
        carry_type == CARRY_TYPE_NONE # there's no carry curve
        || (carry_type == CARRY_TYPE_CURVE && !has_carry_curve(ctx.attr, underlying) ) # there's a carry curve, but not for the current underlying
        || unit.carryless # the unit contract was declared caryless
        )

        # caryless cashflows are discounted using only the riskfree curve
        discount_factor_process = lower!(ctx, riskfree_discount_factor, at)

    elseif carry_type == CARRY_TYPE_CURVE
        carry_discount_factor_observable = DiscountFactor(get_carry_curve_symbol(ctx.attr, underlying), mat)
        riskfree_discount_factor_process = lower!(ctx, riskfree_discount_factor, at)
        carry_discount_factor_process = lower!(ctx, carry_discount_factor_observable, at)
        discount_factor_process = add_instruction!(ctx, OP_BIN_SCALAR_MUL, riskfree_discount_factor_process, carry_discount_factor_process)

    #elseif carry_type == "curve_lookup"
    #    error("carry_type $carry_type not implemented")
        #contract_id = get(attr, :extern_contract_id, missing)
        #counterpart_id = get(attr, :extern_counterpart_id, missing)

    else
        error("Attribute carry_type='$carry_type' not supported.")
    end

    unit_process = lower!(ctx, unit, at)

    # unit_process * discount_factor_process
    return add_instruction!(ctx, OP_BIN_SCALAR_MUL, unit_process, discount_factor_process)
end

# TODO: is this a general rule, or only applies to StaticHedgingModel / FutureValueModel ?
function lower!(ctx::CompilerContext{M}, c::WhenAt{Scale}, state) :: OptimizingIR.ImmutableValue where {M<:StaticHedgingModel}
    assert_at_initial_state(ctx, state)
    pricing_date = get_pricing_date(ctx)

    mat = c.maturity
    scale_contract = c.c
    scale_observable = scale_contract.s
    scale_inner_contract = scale_contract.c

    if pricing_date > mat
        # ignore old cashflows
        return lower!(ctx, Worthless(), state)
    elseif pricing_date == mat
        # returns the underlying asset value if at maturity
        return lower!(ctx, c.c, state)
    end

    # obtains the observable value at the acquisition date
    arg1 = lower!(ctx, scale_observable, mat) # TODO: this should be a problem if we wanted to apply this as a general rule
    arg2 = lower!(ctx, WhenAt(mat, scale_inner_contract), state)

    # arg1 * arg2
    return add_instruction!(ctx, OP_BIN_SCALAR_MUL, arg1, arg2)
end


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

name = "cif_api"
version = v"0.4.2"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/COMCIFS/cif_api/archive/refs/tags/v0.4.2.tar.gz", "803fa1d0525bb51407754fa63b9439ba350178f45372103e84773ed4871b3924"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/cif_api-*

update_configure_scripts

if [[ ${target} == *mingw* ]]; then
    #remove win32 if branch to prevent file not found errors on make install?
    atomic_patch -p1 ${WORKSPACE}/srcdir/patches/mingw-remove-install-hook.patch
    autoreconf -vi
fi

#CPP flags needed to help *-musl-* builds configure scripts find sqlite3.h, unsure exactly why?
export CPPFLAGS="-I${includedir}"

./configure \
    --prefix=${prefix} \
    --libdir=${libdir} \
    --includedir=${includedir} \
    --build=${MACHTYPE} \
    --host=${target} \
    --with-docs=no \

make -j${nproc}
make install
"""

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

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

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency("ICU_jll"; compat="69.1")
    Dependency(PackageSpec(name="SQLite_jll", uuid="76ed43ae-9a5d-5a62-8c75-30186b810ce8"))
]

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

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# compatibility constructors for RANSAC.jl

function toimplicit(fitted::FittedSphere)
    return ImplicitSphere(fitted.center, fitted.radius)
end

function toimplicit(fitted::FittedPlane)
    return ImplicitPlane(fitted.point, fitted.normal)
end

function toimplicit(fitted::FittedCylinder)
    return ImplicitCylinder(fitted.axis, fitted.center, fitted.radius)
end

function toimplicit(fitted::FittedCone)
    return ImplicitCone(fitted.apex, fitted.axis, fitted.opang)
end

function toimplicit(fitted::ExtractedTranslational)
    cf = fitted.coordframe
    segs = fitted.contour
    c = fitted.center
    outw = fitted.outwards
    flipn = fitted.flipnormal
    return ImplicitTranslational(cf, segs, c, outw, flipn)
end

function scored2implicit(pcr::PointCloud, scored::ScoredShape)
    scored.candidate.shape isa FittedPlane && return [toimplicit(scored.candidate.shape)]
    result = AbstractImplicitSurface[]
    append!(result, [toimplicit(scored.candidate.shape)])
    inp = @view pcr.vertices[scored.inpoints]
    _, plusplanes = findOBB(inp)
    append!(result, plusplanes)
    return result
end

"""
    ransacresult2implicit(pcr::PointCloud, scored::Array{ScoredShape,1}, param)

Process ransac result shapes to implicit surfaces.
"""
function ransacresult2implicit(pcr, scored, params)
    implshapes = AbstractImplicitSurface[]
    for i in eachindex(scored)
        append!(implshapes, scored2implicit(pcr, scored[i]))
    end
    toremove = falses(size(implshapes,1))
    for i in eachindex(implshapes)
        for j in i+1:lastindex(implshapes)
            issame(implshapes[i], implshapes[j], params) || continue
            toremove[j] = true
        end
    end
    deleteat!(implshapes, toremove)
    return implshapes, toremove
end

function scored2implicit_ordered(pcr::PointCloud, scored::ScoredShape)
    inp = @view pcr.vertices[scored.inpoints]
    _, plusplanes = findOBB(inp)
    return plusplanes
end

function ransacresult2implicit_ordered(pcr, scored, params)
    implshapes = AbstractImplicitSurface[]
    append!(implshapes, [toimplicit(s.candidate.shape) for s in scored])
    for i in eachindex(scored)
        scored[i].candidate.shape isa FittedPlane && continue
        append!(implshapes, scored2implicit_ordered(pcr, scored[i]))
    end
    toremove = falses(size(implshapes,1))
    for i in eachindex(implshapes)
        for j in i+1:lastindex(implshapes)
            issame(implshapes[i], implshapes[j], params) || continue
            toremove[j] = true
        end
    end
    deleteat!(implshapes, toremove)
    return implshapes, toremove
end

function pairthem(a, b, p)
    @assert size(a) == size(b)
    pairing = []
    for i in 1:size(a,1)
        for j in 1:size(b,1)
            if a[i] == b[j]
                push!(pairing, [i, j])
            elseif CSGBuilding.issame(a[i], b[j], p)
                push!(pairing, [i, j])
            end
        end
    end
    na = CSGBuilding._name.(a)
    anames = [na[i]*string(pairing[i][1])*" " for i in eachindex(na)]
    bnames = [na[i]*string(pairing[i][2])*" " for i in eachindex(na)]

    syms = [string(gensym()) for _ in eachindex(na)]
    # 'name a' -> 'sym'
    asym = [anames[i]=>syms[i] for i in eachindex(anames)]
    # 'sym' -> 'name b'
    symb = [syms[i]=>bnames[i] for i in eachindex(anames)]
    pairing, asym, symb
end