text:
@testset "Thunk" begin
    MutateThunkException = ChainRulesCore.MutateThunkException

    @test @thunk(3) isa Thunk

    @testset "==" begin
        @test @thunk(3.2) == InplaceableThunk(x -> x + 3.2, @thunk(3.2))
        @test @thunk(3.2) == 3.2
        @test 3.2 == InplaceableThunk(x -> x + 3.2, @thunk(3.2))
    end

    @testset "iterate" begin
        a = [1.0, 2.0, 3.0]
        t = @thunk(a)
        for (i, j) in zip(a, t)
            @test i == j
        end

        @test nothing === iterate(@thunk ()) == iterate(())
    end

    @testset "show" begin
        rep = repr(Thunk(rand))
        @test occursin(r"Thunk\(.*rand.*\)", rep)
    end

    @testset "convert" begin
        @test convert(Thunk, ZeroTangent()) isa Thunk
    end

    @testset "unthunk" begin
        @test unthunk(@thunk(3)) == 3
        @test unthunk(@thunk(@thunk(3))) isa Thunk
    end

    @testset "erroring thunks should include the source in the backtrack" begin
        expected_line = (@__LINE__) + 2  # for testing it is at right palce
        try
            x = @thunk(error())
            unthunk(x)
        catch err
            err isa ErrorException || rethrow()
            st = stacktrace(catch_backtrace())
            # Should be 2nd last line, as last line will be the `error` function
            stackframe = st[2]
            @test stackframe.line == expected_line
            @test stackframe.file == Symbol(@__FILE__)
        end
    end

    @testset "Linear operators" begin
        x_real = [2.0 4.0; 8.0 5.0]
        x_complex = [(2.0+im) 4.0; 8.0 (5.0+4im)]
        @testset "$(typeof(x))" for x in (x_real, x_complex)
            x_thunked = @thunk(1.0 * x)
            @test unthunk(x_thunked') == x'
            @test unthunk(transpose(x_thunked)) == transpose(x)
        end
    end

    @testset "Broadcast" begin
        @testset "Array" begin
            was_unthunked = 0
            array_thunk = @thunk begin
                was_unthunked += 1
                [1.0, 2.0, 3.0]
            end

            was_unthunked = 0
            @test array_thunk .+ fill(10, 3) .+ fill(10, 3) == [21.0, 22.0, 23.0]
            @test was_unthunked == 1

            was_unthunked = 0
            @test array_thunk .+ 10.0 .+ 10.0 == [21.0, 22.0, 23.0]
            @test was_unthunked == 1
        end

        @testset "Scalar" begin
            was_unthunked = 0
            scalar_thunk = @thunk begin
                was_unthunked += 1
                sqrt(4.0)
            end

            was_unthunked = 0
            @test scalar_thunk .+ fill(10, 3) .+ fill(10, 3) == [22.0, 22.0, 22.0]
            @test was_unthunked == 1

            was_unthunked = 0
            @test scalar_thunk .+ 10.0 .+ 10.0 == 22.0
            @test was_unthunked == 1
        end
    end

    @testset "basic math" begin
        @test 1 == -@thunk(-1)
        @test 1 == @thunk(2) - 1
        @test 1 == 2 - @thunk(1)
        @test 1.0 == @thunk(1) / 1.0
        @test 1.0 == 1.0 / @thunk(1)

        @test 1 == real(@thunk(1 + 1im))
        @test 1 == imag(@thunk(1 + 1im))
        @test 1 + 1im == Complex(@thunk(1 + 1im))
        @test 1 + 1im == Complex(@thunk(1), @thunk(1))
    end

    @testset "Base functions" begin
        v = [1, 2, 3]
        t = @thunk(v)

        m = rand(3, 3)
        tm = @thunk(m)

        if VERSION >= v"1.2"
            @test 3 == mapreduce(_ -> 1, +, t)
            @test 3 == mapreduce((_, _) -> 1, +, v, t)
        end
        @test 10 == sum(@thunk([1 2; 3 4]))
        @test [4 6] == sum!([1 1], @thunk([1 2; 3 4]))

        @test fill(3.2, 3) == fill(@thunk(3.2), 3)
        @test v == vec(t)
        @test [1 2 3] == reshape(t, 1, 3)
        @test 1 == getindex(t, 1)
        @test_throws MutateThunkException setindex!(t, 0.0, 1)
        @test [4; 5; 6] == selectdim([1 2 3; 4 5 6], 1, 2)

        @test reverse(t) == reverse(v)
        @test reverse(t, 2) == reverse(v, 2)
        @test reverse(tm; dims=2) == reverse(m; dims=2)
    end

    @testset "LinearAlgebra" begin
        v = [1.0, 2.0, 3.0]
        tv = @thunk(v)
        a = [1.0 2.0; 3.0 4.0]
        t = @thunk(a)
        @test Array(a) == Array(t)
        @test Matrix(a) == Matrix(t)
        @test Diagonal(a) == Diagonal(t)
        @test LowerTriangular(a) == LowerTriangular(t)
        @test UpperTriangular(a) == UpperTriangular(t)
        @test Symmetric(a) == Symmetric(t)
        @test Hermitian(a) == Hermitian(t)

        if VERSION >= v"1.2"
            @test diagm(0 => v) == diagm(0 => tv)
            @test diagm(3, 4, 0 => v) == diagm(3, 4, 0 => tv)
            # Check against accidential type piracy
            # https://github.com/JuliaDiff/ChainRulesCore.jl/issues/472
            @test Base.which(diagm, Tuple{}()).module != ChainRulesCore
            @test Base.which(diagm, Tuple{Int,Int}).module != ChainRulesCore
        end
        @test tril(a) == tril(t)
        @test tril(a, 1) == tril(t, 1)
        @test triu(a) == triu(t)
        @test triu(a, 1) == triu(t, 1)
        @test tr(a) == tr(t)
        @test cross(v, v) == cross(v, tv)
        @test cross(v, v) == cross(tv, v)
        @test cross(v, v) == cross(tv, tv)
        @test dot(v, v) == dot(v, tv)
        @test dot(v, v) == dot(tv, v)
        @test dot(v, v) == dot(tv, tv)

        if VERSION >= v"1.2"
            @test_throws MutateThunkException ldiv!(2.0, deepcopy(t)) ==
                                              ldiv!(2.0, deepcopy(a))
            @test_throws MutateThunkException rdiv!(deepcopy(t), 2.0) ==
                                              rdiv!(deepcopy(a), 2.0)
        end

        @test mul!(deepcopy(a), a, a) == mul!(deepcopy(a), t, a)

        res = mul!(deepcopy(a), a, a, true, true)
        @test_throws MutateThunkException mul!(deepcopy(t), a, a, true, true)
        @test_throws MutateThunkException mul!(deepcopy(t), t, a, true, true)
        @test_throws MutateThunkException mul!(deepcopy(t), a, t, true, true)
        @test_throws MutateThunkException mul!(deepcopy(t), t, t, true, true)
        @test res == mul!(deepcopy(a), t, a, true, true)
        @test res == mul!(deepcopy(a), a, t, true, true)
        @test res == mul!(deepcopy(a), t, t, true, true)

        m = rand(3, 3)
        @test ger!(1.0, v, v, deepcopy(m)) == ger!(1.0, tv, v, deepcopy(m))
        @test ger!(1.0, v, v, deepcopy(m)) == ger!(1.0, v, tv, deepcopy(m))
        @test gemv!('C', 1.0, m, v, 1.0, deepcopy(v)) ==
              gemv!('C', 1.0, m, tv, 1.0, deepcopy(v))
        @test gemv('N', 1.0, m, v) == gemv('N', 1.0, m, tv)

        @test scal!(2, 2.0, v, 1) == scal!(2, @thunk(2.0), v, 1)
        @test_throws MutateThunkException LAPACK.trsyl!('C', 'C', m, m, @thunk(m))
    end
end


text:
function nthPower(x, power)
  if (power < 0)
    throw(DomainError(power, "argument must be nonnegative"))
  end

  result = 1
  for i = 1:power
    result *= x
  end
  result
end

function calcTerm(coeffient, term::Tuple{Int, Int})
  coeffient * nthPower(term[1], term[2])
end

function calcPolynomial(tuples...)
  result = 0
  for (tuple) in tuples
    result += calcTerm(tuple[1], (tuple[2], tuple[3]))
  end
  result
end

ext: jl


text:
import Symata: mxpra, MxprArgs

let mx = mxpra(:List, Any[1,2,3])
    @test  mhead(mx) == :List
    @test  margs(mx) == Any[1,2,3]
    @test  typeof(margs(mx)) == MxprArgs
end


text:
using Dates
using Setfield: @set

using TeaFiles.Body: ItemIOBlock, _first_position_ge_time, read_columns, read_items

# TODO These tests could be refactored a bit

function _test_impl_simple(
    times::Vector{Int64},
    func::Function,
    time::Int64,
    expected_index::Int
)
    buf = IOBuffer()
    written_bytes = write(buf, times)
    @test written_bytes == position(buf)
    @test position(buf) == length(times) * sizeof(Int64)

    # The only field is the time.
    block = ItemIOBlock{Int64}(buf, 0, buf.size, sizeof(Int64), 0)

    pos = func(block, time)
    expected_position = (expected_index - 1) * sizeof(Int64)
    @test pos == expected_position
end

struct Moo
    a::UInt8
    time::Int64
end

function _test_impl_padded(
    times::Vector{Int64},
    func::Function,
    time::Int64,
    expected_index::Int
)
    item_size = Int32(sizeof(Moo))
    @test item_size == 16  # Expecting :a to be padded
    time_field_offset = Int32(fieldoffset(Moo, 2))
    @test time_field_offset == 8

    items = [Moo('a', time) for time in times]

    buf = IOBuffer()
    written_bytes = write(buf, items)
    @test written_bytes == position(buf)
    @test position(buf) == length(times) * sizeof(Moo)

    # Use the field offset
    block = ItemIOBlock{Int64}(buf, 0, buf.size, sizeof(Moo), time_field_offset)

    pos = func(block, time)
    expected_position = (expected_index - 1) * sizeof(Moo)
    @test pos == expected_position
end

"""
Test block search with the example metadata.
"""
function _test_impl_example(
    times::Vector{Int64},
    func::Function,
    time::Int64,
    expected_index::Int
)
    metadata = _get_example_metadata()

    buf = IOBuffer()
    metadata_size = write(buf, metadata)
    items = [Tick(time, 10.0, 1) for time in times]
    written_bytes = write(buf, items)
    @test written_bytes == length(times) * sizeof(Tick)

    # The time field offset should be zero
    @test fieldoffset(Tick, 1) == 0

    # item_start is the same as metadata_size.
    block = ItemIOBlock{Int64}(buf, metadata_size, buf.size, sizeof(Tick), 0)

    pos = func(block, time)
    expected_position = metadata_size + (expected_index - 1) * sizeof(Tick)
    @test pos == expected_position
end

function _test_block_search(
    times::Vector{Int64},
    func::Function,
    time::Int64,
    expected_index::Int
)
    @testset "simple" begin _test_impl_simple(times, func, time, expected_index) end
    @testset "padded" begin _test_impl_padded(times, func, time, expected_index) end
    @testset "example" begin _test_impl_example(times, func, time, expected_index) end
end

@testset "read" begin
    @testset "_first_position_ge_time" begin
        @testset "empty" begin
            function _test(time, expected_index)
                times = Int64[]
                return _test_block_search(
                    times, _first_position_ge_time, time, expected_index
                )
            end
            # Regardless of the time we search for, we should find the start of the
            # item block.
            _test(0, 1)
            _test(1, 1)
            _test(2, 1)
        end

        @testset "nonempty" begin
            function _test(time, expected_index)
                times = [10, 10, 20, 30, 50, 50]
                return _test_block_search(
                    times, _first_position_ge_time, time, expected_index
                )
            end
            _test(0, 1)  # This time is before any data we have
            _test(10, 1)
            _test(20, 3)
            _test(30, 4)
            _test(40, 5)
            _test(50, 5)
            _test(60, 7)  # Strictly after all times, so point to end of block
            _test(70, 7)
        end
    end

    @testset "read_items" begin
        # Write the example to the buffer.
        metadata = _get_example_metadata()

        expected_times = [10, 20, 30, 40, 50]
        expected_prices = [12.0, 22.0, 32.0, 42.0, 52.0]
        expected_volumes = [1, 2, 3, 4, 5]

        expected_items = [
            Tick(time, price, volume)
            for (time, price, volume) in zip(
                expected_times, expected_prices, expected_volumes
            )
        ]

        buf = IOBuffer()
        write(buf, metadata)
        write(buf, expected_items)

        @test read_items(Tick, buf, metadata) == expected_items
        @test read_items(Tick, buf, metadata; lower=20) == expected_items[2:end]
        @test read_items(Tick, buf, metadata; lower=40) == expected_items[4:end]
        @test read_items(Tick, buf, metadata; lower=10, upper=40) == expected_items[1:3]
        @test read_items(Tick, buf, metadata; lower=40, upper=41) == expected_items[4:4]
        @test read_items(Tick, buf, metadata; lower=41, upper=42) == Tick[]

        # Verify that an incompatible item type results in an error.
        struct IncompatibleTick
            time::Int64
            price::Float32  # In `Tick`, this is a Float64
            volume::Int64
        end
        @test_throws ArgumentError read_items(IncompatibleTick, buf, metadata)
    end

    @testset "read_items_julia_datetime" begin
        @testset "compatible" begin
            metadata = _get_example_datetime_metadata()

            items = [DateTimeTick(DateTime(2000, 1, 1), 12.0, 1)]
            buf = IOBuffer()
            write(buf, metadata)
            write(buf, items)
            @test read_items(DateTimeTick, buf, metadata) == items

            # We should also be able to write any old integers in a Tick type, then read it
            # back as a DateTime, so long as the epoch is set correctly.
            buf = IOBuffer()
            write(buf, metadata)
            write(buf, [Tick(0, 12.0, 1)])
            @test read_items(DateTimeTick, buf, metadata) == [
                DateTimeTick(DateTime(0, 12, 31), 12.0, 1)
            ]
        end

        @testset "incompatible epoch" begin
            metadata = _get_example_metadata()
            buf = IOBuffer()
            write(buf, metadata)
            write(buf, [Tick(10, 12.0, 1)])

            # This should throw an error, because even though the `Tick` class stores times
            # as an Int64, which is binary compatible with a Julia DateTime, and each tick
            # is a millisecond, it has a different epoch.
            @test_throws ArgumentError read_items(DateTimeTick, buf, metadata)
        end

        @testset "incompatible ticks per day" begin
            metadata = _get_example_datetime_metadata()
            time_section = @set get_section(metadata, TimeSection).ticks_per_day = 1
            _replace_section!(metadata, time_section)

            buf = IOBuffer()
            write(buf, metadata)
            write(buf, [Tick(10, 12.0, 1)])

            @test_throws ArgumentError read_items(DateTimeTick, buf, metadata)
        end
    end

    @testset "read_items_namedtuple" begin
        @testset "compatible" begin
            metadata = _get_example_datetime_metadata()

            items = [DateTimeTick(DateTime(2000, 1, 1), 12.0, 1)]
            buf = IOBuffer()
            write(buf, metadata)
            write(buf, items)
            @test read_items(buf, metadata) == _structs_to_namedtuples(items)
        end
    end

    @testset "read_columns" begin
        # Write the example to the buffer.
        metadata = _get_example_metadata()

        expected_times = [10, 20, 30, 40, 50]
        expected_prices = [12.0, 22.0, 32.0, 42.0, 52.0]
        expected_volumes = [1, 2, 3, 4, 5]

        expected_items = [
            Tick(time, price, volume)
            for (time, price, volume) in zip(
                expected_times, expected_prices, expected_volumes
            )
        ]

        buf = IOBuffer()
        write(buf, metadata)
        write(buf, expected_items)

        times, prices, volumes = read_columns(buf, metadata)
        @test times == expected_times
        @test prices == expected_prices
        @test volumes == expected_volumes

        @test only(read_columns(buf, metadata; columns=[1])) == expected_times
        @test only(read_columns(buf, metadata; columns=[2])) == expected_prices
        @test only(read_columns(buf, metadata; columns=[3])) == expected_volumes

        @test (
            read_columns(buf, metadata; columns=[1, 2])
            == [expected_times, expected_prices]
        )
        @test (
            read_columns(buf, metadata; columns=[1, 3])
            == [expected_times, expected_volumes]
        )
        @test (
            read_columns(buf, metadata; columns=[2, 3])
            == [expected_prices, expected_volumes]
        )
    end
end


text:
normalise_ext(f) = endswith(f, ".bson.zstd" ) ? f : f * ".bson.zstd"
function bsave(fname, dict)
    open(normalise_ext(fname), "w") do f
        s = ZstdCompressorStream(f)
        bson(s, dict)
        close(s)
    end
end
function bload(fname, mod = @__MODULE__)
    open(normalise_ext(fname)) do f
        s = ZstdDecompressorStream(f)
        res = BSON.load(s, mod)
        close(s)
        res
    end
end


text:
"""
    struct SubMesh{N,T} <: AbstractMesh{N,T}

Create a view of a `parent` mesh over a given `domain`.

A submesh implements the interface for `AbstractMesh`; therefore you can iterate
over elements of the submesh just like you would with a mesh.
"""
struct SubMesh{N,T} <: AbstractMesh{N,T}
    parent::GenericMesh{N,T}
    domain::Domain
    dom2elt::Dict{DataType,Vector{Int}}
    function SubMesh{N,T}(mesh::GenericMesh,Ω::Domain) where {N,T}
        idxs = dom2elt(mesh,Ω)
        return new{N,T}(mesh,Ω,idxs)
    end
end
SubMesh(m::GenericMesh{N,T},args...;kwargs...) where {N,T} = SubMesh{N,T}(m,args...;kwargs...)

Base.view(m::GenericMesh,Ω::Domain)           = SubMesh(m,Ω)
Base.view(m::GenericMesh,ent::AbstractEntity) = SubMesh(m,Domain(ent))

# TODO: write tests for this
Base.view(m::SubMesh,Ω::Domain)           = view(mesh(m),intersect(Ω,domain(m)))
Base.view(m::SubMesh,ent::AbstractEntity) = view(m,Domain(ent))

Base.parent(m::SubMesh) = m.parent
domain(m::SubMesh) = m.domain

# ElementIterator for submesh
function Base.size(iter::ElementIterator{<:AbstractElement,<:SubMesh})
    E          = eltype(iter)
    submesh    = mesh(iter)
    idxs       = dom2elt(submesh,E)
    return (length(idxs),)
end

function Base.getindex(iter::ElementIterator{<:AbstractElement,<:SubMesh},i::Int)
    E      = eltype(iter)
    submsh = mesh(iter) # a SubMesh
    p_msh  = parent(submsh) # parent mesh
    iglob  = dom2elt(submsh,E)[i] # global index of element in parent mesh
    iter   = ElementIterator(p_msh,E) # iterator over parent mesh
    return iter[iglob]
end

function Base.iterate(iter::ElementIterator{<:AbstractElement,<:SubMesh}, state=1)
    state > length(iter) && (return nothing)
    iter[state], state + 1
end

"""
    dom2elt(m::SubMesh,[E])

A dictionary with keys being the element types of `m`, and values being the
element indices in the parent mesh. If a type `E` is given, return the values
associated with that key.
"""
dom2elt(m::SubMesh) = m.dom2elt
dom2elt(m::SubMesh,E::Type{<:AbstractElement}) = m.dom2elt[E]

function Base.keys(submesh::SubMesh)
    Ω, M = submesh.domain, submesh.parent
    ee = DataType[]
    for ent in entities(Ω)
        dict = M.ent2tags[ent]
        append!(ee, keys(dict))
    end
    return unique!(ee)
end

max_stars_repo_name: WaveProp/WavePropBase

max_issues_repo_name: WaveProp/WavePropBase

max_forks_repo_name: WaveProp/WavePropBase


text:
good = KNNClassifier(K=2)
bad = KNNClassifier(K=10)
ugly = ConstantClassifier()

r = [good, bad, ugly]

rng = StableRNGs.StableRNG(123)
X, y = make_blobs(rng=rng)

@testset_accelerated "integration" accel begin

    # evaluate the three models separately:
    resampling = Holdout(fraction_train=0.6)
    scores = map(r) do model
        e = evaluate(model, X, y, resampling=resampling, measure=LogLoss())
        return e.measurement[1]
    end
    @test scores == sort(scores)

    # find the scores using `Explicit` tuning:
    tmodel = TunedModel(models=r,
                        resampling=resampling,
                        measure=LogLoss(),
                        acceleration=accel)
    mach = machine(tmodel, X, y)
    fit!(mach, verbosity=0)

    history = report(mach).history
    _models = first.(history)
    _scores = map(history) do entry
        entry.measurement[1]
    end
    ms = sort(zip(_models, _scores) |> collect, by=last)

    # check ordering is best to worst again:
    @test first.(ms) == r

    # check scores are the same:
    @test last.(ms) ≈ scores

    # fail with ArgumentError when model types are wrong (e.g., are not instantiated).
    # this used to throw a very confusing MethodError.
    dcc = DeterministicConstantClassifier
    @test_throws ArgumentError TunedModel(; models=[dcc, dcc])
end

true

max_stars_repo_name: alan-turing-institute/MLJTuning

max_issues_repo_name: alan-turing-institute/MLJTuning

max_forks_repo_name: alan-turing-institute/MLJTuning


text:
# It is suggested that test code for MLJ.jl include files be placed in
# a file of the same name under "test/" (and included below) and that
# this test code be wrapped in a module. Any new module name will do -
# eg, `module TestDatasets` for code testing `datasets.jl`.

using MLJBase, Test

@testset "scientific trait" begin
    @test include("scientific_trait.jl")
end

@testset "equality" begin
  @test include("equality.jl")
end

@testset "static transforms" begin
  @test include("static.jl")
end

@testset "utilities" begin
  @test include("utilities.jl")
end

@testset "parameters" begin
  @test include("parameters.jl")
end

@testset "distributions" begin
  @test include("distributions.jl")
end

@testset "data" begin
  @test include("data.jl")
end

@testset "info" begin
  @test include("info.jl")
end

@testset "datasets" begin
  @test include("datasets.jl")
end

@testset "tasks" begin
  @test include("tasks.jl")
end

@testset "measures" begin
  @test include("measures.jl")
end

@testset "interface for LossFunctions" begin
  @test include("loss_functions_interface.jl")
end

@testset "@mlj_model" begin
  @test include("mlj_model_macro.jl")
end

@testset "metadatautils" begin
  @test include("metadata_utilities.jl")
end


text:
# Copyright 2019 Tobias Frilling
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

abstract type ABCanalysis end

"""
    ABCanalysis4Curve(curve::ABCcurve) <: ABCanalysis

Calculate an ABC analysis for the given curve.

# Fields
- `pareto`: Nearest point to a theoretically ideal Break Even point.
- `break_even`: Point on the curve where the gain (``dABC``) is 1.
- `demark_AB` Point on the curve at which most of the yield is already obtained; the smaller of the Pareto and Break Even points. Demarkation between *A* and *B*.
- `submarginal`: Point on the curve after which the gain can be considered trivial. Demarkation between *B* and *C*.
- `curve`: The given curve. Used for plotting.
"""
struct ABCanalysis4Curve <: ABCanalysis
    pareto     ::Tuple{Float64, Float64}
    break_even ::Tuple{Float64, Float64}
    demark_AB  ::Tuple{Float64, Float64}
    submarginal::Tuple{Float64, Float64}
    curve::ABCcurve

    function ABCanalysis4Curve(curve::ABCcurve)
        # euclidean distance (without sqrt)
        dist(p1,p2) = (p1 .- p2).^2 |> sum
        zipped_curve = zip(curve.effort, curve.yield)
        # detemines the index of the closest point to p
        index_min_dist_point(p) = map(x -> dist(p,x), zipped_curve) |> argmin

        abc_ideal = (0,1)
        pareto_index = index_min_dist_point(abc_ideal)

        gradients = map(i->Interpolations.gradient(curve.interpolation, i)[1], curve.effort)
        break_even_index =  abs.(gradients .- 1) |> argmin

        demark_AB_index = min(pareto_index, break_even_index)

        break_even_ideal = (curve.effort[max(pareto_index, break_even_index)], 1)
        submarginal_index = index_min_dist_point(break_even_ideal)

        (effort, yield) = (curve.effort, curve.yield)
        point(indx) = (effort[indx], yield[indx])

        new(point(pareto_index),
            point(break_even_index),
            point(demark_AB_index),
            point(submarginal_index),
            curve)
    end
end

"""
    ABCanalysis(curve::ABCcurve)

Convenience function for `ABCanalysis4Curve(curve)`
"""
ABCanalysis(curve::ABCcurve) = ABCanalysis4Curve(curve)

"""
    ABCanalysis4Data(data::Vector{<:Real}) <: ABCanalysis

Calculate an ABC analysis for the given data.

# Fields
Extends `ABCanalysis4Curve` with:
- `a_indices`: Indices of elements in the *A* set.
- `b_indices`: Indices of elements in the *B* set.
- `c_indices`: Indices of elements in the *C* set.
"""
struct ABCanalysis4Data <: ABCanalysis
    pareto     ::Tuple{Float64, Float64}
    break_even ::Tuple{Float64, Float64}
    demark_AB  ::Tuple{Float64, Float64}
    submarginal::Tuple{Float64, Float64}

    a_indices::Vector{Int64}
    b_indices::Vector{Int64}
    c_indices::Vector{Int64}

    curve::ABCcurve

    function ABCanalysis4Data(data::Vector{<:Real})
        curve = ABCcurve(data)
        ana = ABCanalysis4Curve(curve)
        perm = sortperm(data, rev=true)

        # have to subtract 1 because we added (0,0) to the curve data
        point2index(p) = findfirst(isequal(p[1]), curve.effort) - 1
        threshold_a = point2index(ana.demark_AB)
        threshold_b = point2index(ana.submarginal)

        new(ana.pareto,
            ana.break_even,
            ana.demark_AB,
            ana.submarginal,
            perm[1:threshold_a],
            perm[threshold_a+1:threshold_b],
            perm[threshold_b+1:end],
            curve)
    end
end

"""
    ABCanalysis(data::Vector{<:Real})

Convenience function for `ABCanalysis4Data(data)`
"""
ABCanalysis(data::Vector{<:Real}) = ABCanalysis4Data(data)

# text representation for REPL or Jupyter
function show_points(io::IO, analysis::ABCanalysis)
    println(io, "Pareto:          ", analysis.pareto)
    println(io, "Break Even:      ", analysis.break_even)
    println(io, "Demarkation A|B: ", analysis.demark_AB)
    print(io,   "Submarginal:     ", analysis.submarginal)
end

function Base.show(io::IO, ::MIME"text/plain", analysis::ABCanalysis)
    show_points(io, analysis)
end

function Base.show(io::IO, ::MIME"text/plain", analysis::ABCanalysis4Data)
    show_points(io, analysis)
    println(io)
    println(io, "A: $(analysis.a_indices |> length) elements")
    println(io, "B: $(analysis.b_indices |> length) elements")
    print(io, "C: $(analysis.c_indices |> length) elements")
end

# plotting
@recipe function f(ana::ABCanalysis; comparison=true, annotate=true)
    xlims --> (0,1)
    ylims --> (0,1)
    legend --> :bottomright
    ratio --> 1
    markersize --> 5

    if annotate
        # plot annotations
        fontsize = 8
        # these are the 'n = ???' in the plot
        a_size = findfirst(isequal(ana.demark_AB[1]), ana.curve.effort)
        b_size = findfirst(isequal(ana.submarginal[1]), ana.curve.effort) - a_size
        c_size = length(ana.curve.effort) - a_size - b_size

        # position for A and B is centred between set boundaries
        offset_y = 0.1
        a_xpos = ana.demark_AB[1] / 2
        b_xpos = ana.demark_AB[1] + (ana.submarginal[1] - ana.demark_AB[1])/2
        c_xpos = ana.submarginal[1] + 0.1

        annotations := [
                        # a_size - 1 because we added (0,0) in Curve
                        (a_xpos, offset_y, text("A\nn = $(a_size-1)", fontsize))
                        (b_xpos, offset_y, text("B\nn = $(b_size)", fontsize))
                        (c_xpos, offset_y, text("C\nn = $(c_size)", fontsize))
                       ]
    end

    @series begin
        comparison := comparison
        ana.curve
    end

    @series begin
        # set boundaries
        label := ""
        linecolor := :red
        seriestype := :path
        markershape := :none
        # path for the boundary lines. NaN starts a new line
        path = [
                0                  ana.demark_AB[2]
                ana.demark_AB[1]   ana.demark_AB[2]
                ana.demark_AB[1]   0
                NaN                NaN
                0                  ana.submarginal[2]
                ana.submarginal[1] ana.submarginal[2]
                ana.submarginal[1] 0
               ]
        path[:,1], path[:,2]
    end

    @series begin
        # marker for the analysis points
        label := ""
        seriestype := :scatter
        markercolors --> [:red, :green, :blue]
        [ana.break_even, ana.pareto, ana.submarginal]
    end
end


text:
"""
Genetic Algorithm
Authors: Carmine Spagnuolo and Giuseppe D'Ambrosio

"""
tstartglobal=time()
using Distributed
using eqjl
using pool

# uncomment if you need to install the packages

# import Pkg
# Pkg.add("Random")
using Random
# Pkg.add("Distributions")
using Distributions
# Pkg.add("JSON")
using JSON

# object function executed from workers
# the string should contain the path to the folder containing the function
# and the name of the module, separated by ;
# ex: path/to/folder;moduleName
const objfunc = string(pwd(), "/julia;objfunc")

# objective function
function objfuncF(pops)
	x = pops[1]
    y = pops[2]
    sin(4 * x) + sin(4 * y) + -2 * x + x^2 - 2 * y + y^2
end

# tournament selection
function selection(pop, groupSize)
	N = length(pop)
	selection = fill(0,N)
	nFitness = length(pop)

	for i in 1:N
		# first select a random contender
		contender = unique(rand(1:nFitness, groupSize))
		while length(contender) < groupSize
			contender = unique(vcat(contender, rand(1:nFitness, groupSize - length(contender))))
		end

		# find the best
		winner = first(contender)
		winnerFitness = pop[winner]
		for idx = 2:groupSize
			c = contender[idx]
			if winnerFitness < pop[c]
				winner = c
				winnerFitness = pop[c]
			end
		end

		selection[i] = winner
	end
	return selection
end

# uniform (binomial) crossover function
# crossover two parents to create two children
# cross_pb must be between [0,1]
function crossover(p1, p2, cross_pb)
	# children are copies of parents by default
	l = length(p1)
    c1 = copy(p1)
    c2 = copy(p2)
	j = rand(1:l)
	# check for recombination
	for i in (((1:l).+j.-2).%l).+1
		if rand() <= cross_pb
			val = p1[i]
    		p1[i] = p2[i]
    		p2[i] = val
		end
	end
	return c1, c2
end

# mutation operator
function mutation(recombinant, sigma)
	d = Distributions.Normal(0, sigma)
	recombinant .+= rand(d)
	return recombinant
end

# genetic algorithm search of the one max optimization problem
function genetic_algorithm()
	println("algorithm.jl starting genetic_algorithm()")
	global res

    # Getting settings filename
    eqjl.OUT_put("Params")
    algo_params_file = eqjl.IN_get()

	# Reading settings from file
    params = JSON.parsefile(algo_params_file)
    Random.seed!(params["seed"])
    num_iter = params["num_iter"]
    num_pop = params["num_pop"]
	sigma = params["sigma"]
	mate_pb = params["mate_pb"]
    mutate_pb = params["mutate_pb"]
    crossover_pb = params["crossover_pb"]
	e = params["e"]

	tmp_dir = joinpath(ENV["TURBINE_OUTPUT"], "tmp")

	# JLMap Pool
    global _pool = Pool(tmp_dir, "workers")
    @eval @everywhere _pool = _pool

	# Creating initial population of random Float
	individuals = ()->Float64[rand(0.0:0.01:5.0), rand(0.0:0.01:5.0)]
    pop = [individuals() for _ in (1:num_pop)]

	# Enumerate generations
	for gen in (1:num_iter)
		# Evaluate all candidates in the population using JLMap
		fitness = mapWorker(_pool, objfunc, pop)
		minFitness, fitidx = findmin(fitness)
		offspring = similar(pop)
		selected = selection(fitness, floor(Int, num_pop/2))

		# perform matingstate fitness
		offidx = Random.randperm(num_pop)
		eliteSize = isa(e, Int) ? e : round(Int, e * num_pop)
		offspringSize = num_pop - eliteSize

		for i in 1:2:offspringSize
			ij = (i == offspringSize) ? i-1 : i+1
			if rand() < crossover_pb
				offspring[i], offspring[ij] = crossover(pop[selected[offidx[i]]], pop[selected[offidx[ij]]], crossover_pb)
			else
				offspring[i], offspring[ij] = pop[selected[i]], pop[selected[ij]]
			end
		end

		fitidxs = sortperm(fitness)

		for j in 1:eliteSize
			subs = offspringSize+j
			offspring[subs] = copy(pop[fitidxs[j]])
		end

		for k in 1:offspringSize
			if rand() < mutate_pb
				mutation(offspring[k], sigma)
			end
		end

		for ik in 1:num_pop
			pop[ik] = offspring[ik]
		end

		fitness = mapWorker(_pool, objfunc, offspring)
		minfit, fitidx = findmin(fitness)
		best = pop[fitidx]
		best_eval = fitness[fitidx]
		res = string("f(", best, ") = ", best_eval)
	end
	tend=time()
	res = string(res, ", total elapsed: ", tend-tstartglobal)
    eqjl.OUT_put("DONE")
    println("###DONE")
    eqjl.OUT_put(res)
end

@async genetic_algorithm()

max_stars_repo_name: emews/EQ-JL

max_issues_repo_name: emews/EQ-JL

max_forks_repo_name: emews/EQ-JL


text:
using .JSort
using Statistics
using Distributions: Poisson, rand
import PyPlot; const plt = PyPlot
import PyCall
using ImageFiltering: imfilter, imfilter!, Kernel
using Base.Filesystem: splitext

function featurealign2d(path::String, pattern::Regex=r"edeb(?<b>\d)f(?<f>\d).bin", referenceid=1
                        ;weight::Weight=ThresholdSum(.3), metric::Function=L2, savepath::String="",
                        neighbourwidth=9, neighbourthreshold=80, plot=false)
    # Set savepath to `path` if it is not specified.
    # Save figures in a .figures subfolder
    if length(savepath) == 0
        savepath = path
    end
    savepath = joinpath(savepath, "figures")
    println("Saving to $savepath")
    mkpath(savepath)
    # Sort all back detectors corresponding to common front detectors
    frontdetectors = Dict(front => Dict{Int,String}() for front in 1:8)
    # Find all files
    anymatch = false
    for fname in readdir(path)
        m = match(pattern, fname)
        if m !== nothing
            f, b = parse.(Int, (m[:f], m[:b]))
            frontdetectors[f][b] = fname
            anymatch = true
        end
    end

    if !anymatch
        throw("Found no matching files using $pattern in $path")
    end

    model = [(x, y, w) -> leastsquares(x, y;order= 1)]
    name = ["Linear Regression"]

    lm = pyimport("sklearn.linear_model")
    loss = ["absolute_loss", "squared_loss"]
    models = [(x, y, w) -> begin
                  X = reshape(x, :, 1)
                  lr.fit(X, y, w)
                  return [lr.estimator_.intercept_; lr.estimator_.coef_]
              end
              for lr in [lm.RANSACRegressor(loss=l) for l in loss]]
    names = ["RANSAC "*l for l in loss]
	models = [model..., models...]
	names = [name..., names...]

    joinp(x) = joinpath(path, x)

    coefficientsx = [Float64[] for f in 1:8, b in 1:8]
    coefficientsy = [Float64[] for f in 1:8, b in 1:8]

    for (front, backdetectors) in frontdetectors
        reference_fname = backdetectors[referenceid]
        reference_df = ParticleDetector(reference_fname |> joinp)
        #edges = histogram_edges(reference_df[!, :e], reference_df[!, :Δe], nbins=1000)
        edges = histogram_edges(reference_df.e, reference_df.Δe, nbins=1000)
        reference_raw = JSort.histogram(reference_df.e, reference_df.Δe, edges...)
        # fig, ax = plt.subplots()
        # ax.matshow(transpose(reference_raw) .|> log10)
        # ax.invert_yaxis()
        # ax.set_title("Reference #$referenceid - Raw")
        reference = neighbourreduce(reference_raw, width=neighbourwidth, threshold=neighbourthreshold)
        # fig, ax = plt.subplots()
        # ax.matshow(transpose(reference) .|> log10)
        # ax.set_title("Reference #$referenceid - Reduced")
        if plot
        # ax.invert_yaxis()
          total = copy(reference_raw)
          alltotal = (e=copy(reference_df.e), Δe=copy(reference_df.Δe))
        end
        for (back, fname) in backdetectors
            if back == referenceid
                coefficientsx[front, back] = [0.0, 1.0]
                coefficientsy[front, back] = [0.0, 1.0]
                continue
            end
            target_df = ParticleDetector(fname |> joinp)
            target_raw = JSort.histogram(target_df.e, target_df.Δe, edges...)
            # fig, ax = plt.subplots()
            # ax.matshow(transpose(target_raw) .|> log10)
            # ax.set_title("Target b#$back - Raw")
            # ax.invert_yaxis()
            target = neighbourreduce(target_raw, width=neighbourwidth, threshold=80)

            points = featurealign2d_weighted(reference, target, numregions=15,
                                             limit=1000, searchwidth=40, weight=weight,
                                             metric=metric)
            #@show points
            rpx, tpx, rpy, tpy, weights = points
            # annotate_plot_2d(target, tpx, tpy, weights)
            # plt.show()
            # return
            # Convert from index basis to energy basis
            rex, tex = edges[1][rpx], edges[1][tpx]
            rey, tey = edges[2][rpy], edges[2][tpy]

			optimal = ()
			for (name, model) in zip(names, models)

                #coeffx, Σx = leastsquares(tex, rex, weights; order=1)
                #coeffy, Σy = leastsquares(tey, rey, weights; order=1);
				coeffx = model(tex, rex, weights)
				coeffy = model(tey, rey, weights)

				px = Poly(coeffx)
				py = Poly(coeffy)

                e =  px.(target_df.e)
                de = py.(target_df.Δe)

                calibrated = JSort.histogram(e, de, edges...);

                err_un = compare_2d(reference, target, metric=metric)
                err_ca = compare_2d(reference, calibrated, metric=metric)
				ratio = err_ca/err_un
				println("Error of $name: $ratio")

				if length(optimal) == 0
                    optimal =  (name, err_ca, err_un, coeffx, coeffy, calibrated)
				elseif optimal[2] > err_ca
					optimal = (name, err_ca, err_un, coeffx, coeffy, calibrated)
                end
			end
            println("Optimal is $(optimal[1]) with $(optimal[2])")
            name, err_ca, err_un, coeffx, coeffy, calibrated = optimal
            coefficientsx[front, back] = coeffx
            coefficientsy[front, back] = coeffy
            if plot
                total .+= calibrated
                append!(alltotal.e, Poly(coeffx).(target_df.e))
                append!(alltotal.Δe, Poly(coeffy).(target_df.Δe))

              fig, ax = plt.subplots(ncols=2, nrows=2, sharex=true, sharey=true,
                                    figsize=(20, 20))
              # ax[1].matshow(transpose(reference) .|> log10)
              # ax[1].scatter(rpx, rpy, marker="x", c="r", s=0.5)
              # ax[1].invert_yaxis()
              annotate_plot_2d(reference_raw, rpx, rpy, ax=ax[1])
              ax[1].set_title("f$front b$referenceid")
              # ax[2].matshow(transpose(target) .|> log10)
              # ax[2].scatter(tpx, tpy, marker="x", c="r", s=0.5)
              annotate_plot_2d(target_raw, tpx, tpy, weights, ax=ax[2])
              ax[2].set_title("f$front b$back")
              ax[2].invert_yaxis()
              ax[3].matshow(transpose(reference_raw) .|> log10, cmap="Reds")
              ax[3].matshow(transpose(target_raw) .|> log10, cmap="Blues")
              ax[3].annotate(string(err_un), xy=(1,0), xycoords="axes fraction",
                            horizontalalignment="right", verticalalignment="bottom")
              ax[3].set_title("Original")
              ax[3].invert_yaxis()
              ax[4].matshow(transpose(reference_raw) .|> log10, cmap="Reds")
              ax[4].matshow(transpose(calibrated) .|> log10, cmap="Blues")
              ax[4].set_title("Calibrated - "*name)
              ax[4].invert_yaxis()
              ax[4].annotate(string(err_ca), xy=(1,0), xycoords="axes fraction",
                            horizontalalignment="right", verticalalignment="bottom")
              fig.savefig(joinpath(savepath, "edef$(front)b$(back).png"), bbox_inches="tight", dpi=196)
              plt.close();

              fig, ax = plt.subplots(ncols=2, nrows=2, figsize=(10, 10))
              plotresiduals(rpx, tpx, coeffx, ax=ax[1:2])
              plotresiduals(rpy, tpy, coeffy, ax=ax[3:4])
              fig.savefig(joinpath(savepath, "residuef$(front)b$(back).png"), bbox_inches="tight", dpi=196)
              plt.close();
            end
        end
        if plot
          fig, ax = plt.subplots()
          ax.matshow(total' .|> log10)
          ax.invert_yaxis()
          fig.savefig(joinpath(savepath, "edef$front.png"), bbox_inches="tight", dpi=196)
          plt.close();
            fig, ax = plt.subplots()
            t = histogram(alltotal.e, alltotal.Δe, edges...)
            ax.matshow(t' .|> log10)
            ax.invert_yaxis()
            fig.savefig(joinpath(savepath, "tedef$front.png"), bbox_inches="tight", dpi=196)
            plt.close();
            p = ParticleDetector(alltotal.e, alltotal.Δe, 1, 1)
            write(joinpath(savepath, "edef$front.bin"), p)
        end
    end
    coefficientsx, coefficientsy
end

function featurealign2d(reference, target; numregions=5, width=nothing,
                        searchwidth=50, smoother=nothing, limit=50)
    if !isnothing(smoother)
        throw("Smoothing not implemented")
    end
    reference = copy(reference)
    target = copy(target)
    subarrays = splitregions2d(reference, numregions, width=width, limit=limit)
    imfilter!(target, target, Kernel.gaussian(3))
    target .= target ./ sum(target)
    imfilter!(reference, reference, Kernel.gaussian(3))
    reference .= reference ./ sum(reference)

    targetpoints_x    = Int64[]
    targetpoints_y    = Int64[]
    referencepoints_x = Int64[]
    referencepoints_y = Int64[]
    i = 0
    for array in subarrays
        x, y = featureoverlap2d(array, target, searchwidth)
        rx, ry = array.indices
        rm = (rx, ry) .|> middle

        push!(targetpoints_x, x)
        push!(targetpoints_y, y)
        push!(referencepoints_x, rm[1])
        push!(referencepoints_y, rm[2])
    end
    referencepoints_x, targetpoints_x, referencepoints_y, targetpoints_y
end

function featurealign2d_weighted(reference, target; numregions=5, width=nothing,
                                 searchwidth=50, smoother=nothing, limit=50,
                                 metric::Function=L2,
                                 weight::Weight=ThresholdArea(.2))
    if !isnothing(smoother)
        throw("Smoothing not implemented")
    end
    reference = copy(reference)
    target = copy(target)
    subarrays = splitregions2d(reference, numregions, width=width, limit=limit)
    imfilter!(target, target, Kernel.gaussian(3))
    target .= target ./ sum(target)
    imfilter!(reference, reference, Kernel.gaussian(3))
    reference .= reference ./ sum(reference)

    targetpoints_x    = Int64[]
    targetpoints_y    = Int64[]
    referencepoints_x = Int64[]
    referencepoints_y = Int64[]
    weights           = Float64[]
    i = 0
    for array in subarrays
        x, y, w = featureoverlap2d_weighted(array, target, searchwidth, weight, metric)
        rx, ry = array.indices
        rm = (rx, ry) .|> middle

        push!(targetpoints_x, x)
        push!(targetpoints_y, y)
        push!(referencepoints_x, rm[1])
        push!(referencepoints_y, rm[2])
        push!(weights, w)
    end
    referencepoints_x, targetpoints_x, referencepoints_y, targetpoints_y, weights
end

function featurealign2d_ensemble(reference, target; numregions=5, width=nothing,
                        searchwidth=50, smoother=nothing, limit=50)
    if !isnothing(smoother)
        throw("Smoothing not implemented")
    end
    reference = copy(reference)
    target = copy(target)
    subarrays = splitregions2d(reference, numregions, width=width, limit=limit)
    target = imfilter(target, Kernel.gaussian(3))
    target .= target ./ sum(target)
    N = sum(reference)

    targetpoints_x    = Int64[]
    targetpoints_y    = Int64[]
    referencepoints_x = Int64[]
    referencepoints_y = Int64[]
    i = 0
    for array in subarrays
        #tx, ty = featureoverlap2d(array, target, searchwidth)
        featureoverlap2d_ensemble(array, target, searchwidth, N)
        # println(ty, tx)
        # rx, ry = array.indices
        # rm = (rx, ry) .|> middle
        # weight by difference in region total sum?
        if (i += 1) > 1
            return
        end

        # push!(targetpoints_x, tm[1])
        # push!(targetpoints_y, tm[2])
        # push!(referencepoints_x, rm[1])
        # push!(referencepoints_y, rm[2])
    end
    referencepoints_x, targetpoints_x, referencepoints_y, targetpoints_y
end

function featureoverlap2d(feature, target, searchwidth, metric::Function=L2)
    window, _ = shiftsinrange(feature, target, searchwidth)
    distance = slidingwindow(feature, target, metric, searchwidth)
    min = argmin(distance)
    fig, ax = plt.subplots()
    ax.matshow(distance .|> log10)
    ax.contour(distance)
    ax.scatter(min[2]-1, min[1]-1, marker="x", c="r")
    # fig, ax = plt.subplots(ncols=2)
    # ax[1].matshow(feature)
    # X, Y = window[1] .+ min[1], window[2] .+ min[2]
    # ax[2].matshow(target[X, Y])
    # Convert from surface coordinates to feature indices
    offset = first.(window) .+ (length.(window) .÷ 2)
    return min[1] + offset[1], min[2] + offset[2]
end

function featureoverlap2d_weighted(feature, target, searchwidth, weight::Weight, metric::Function=L2)
    window, _ = shiftsinrange(feature, target, searchwidth)
    distance = slidingwindow(feature, target, metric, searchwidth)
    min = argmin(distance)
    # fig, ax = plt.subplots()
    # ax.matshow(distance .|> log10)
    # ax.contour(distance)
    # ax.scatter(min[2]-1, min[1]-1, marker="x", c="r")
    # fig, ax = plt.subplots()
    # d = distance .<= 1.2.*distance[min]
    # ax.matshow(d)
    w = weight(distance)
    # fig, ax = plt.subplots(ncols=2)
    # ax[1].matshow(feature)
    # X, Y = window[1] .+ min[1], window[2] .+ min[2]
    # ax[2].matshow(target[X, Y])
    # Convert from surface coordinates to feature indices
    offset = first.(window) .+ (length.(window) .÷ 2)
    return min[1] + offset[1], min[2] + offset[2], w
end

function featureoverlap2d_ensemble(feature, target, searchwidth, N::Number, members=1, metric=χ²)
    minima = []
    perturbed = similar(feature)
    #for member in ProgressBar(1:members)
    window, _ = shiftsinrange(feature, target, searchwidth)
    metric = L2
    for member in 1:members
        # Perturb the target to add some randomness
        perturbate!(perturbed, feature)
        perturbed = imfilter(perturbed, Kernel.gaussian(3))
        perturbed ./= N
        distance = slidingwindow(perturbed, target, feature.indices, metric, searchwidth)
        min = argmin(distance)
        fig, ax = plt.subplots()
        ax.matshow(distance)
        ax.scatter(min[2], min[1])
        fig, ax = plt.subplots(ncols=2)
        ax[1].matshow(perturbed)
        X, Y = window[1] .+ min[1], window[2] .+ min[2]
        ax[2].matshow(target[X, Y])
        push!(minima, min)
    end
    println(minima[1])
    # Convert from surface coordinates to feature indices
    return minima
    # Transform list of ranges to matrix of midpoints
    # mid = [m .|> middle for m in minima]
    # mid = transpose(hcat(mid...))
    # fig, ax = subplots()
    # ax.scatter(mid[:, 1], mid[:, 2], color="r", alpha=0.5)
    # midpoint = round.(Int, median(mid, dims=1))
    # ax.scatter(midpoint[1], midpoint[2], color="g")
    # x, y = feature.indices
    # wx = (x[end] - x[1])÷2
    # wy = (y[end] - y[1])÷2
    # mid = midpoint
    # [(mid[1]-wx):(mid[1]+wx), (mid[2]-wy):(mid[2]+wy)]
end

function automap(mat; width=50, searchwidth=100, metric=L2, weight=WellSum(4),
                 limit=50)
    score = similar(mat)
    I, J = size(mat)
    w = width÷2 |> Int
    sw = searchwidth÷2 |> Int
    template = @view mat[1:2, 1:2]
    searcharea = @view mat[1:2, 1:2]
    scores = Matrix(undef, searchwidth-width, searchwidth-width)
    @inbounds for i in 1:I, j in 1:J
        # Skip points `width` from the edge
        if !(searchwidth ≤ j ≤ (J-searchwidth)) || !(searchwidth ≤ i ≤ (I-searchwidth))
            score[i, j] = zero(eltype(mat))
            continue
        end
        template = @view mat[i-w:i+w-1, j-w:j+w-1]
        if sum(template) < limit
            score[i, j] = zero(eltype(mat))
            continue
        end
        searcharea = @view mat[i-sw:i+sw-1, j-sw:j+sw-1]
        slidingwindow!(scores, searcharea, template, metric)
        score[i, j] = weight(scores)
    end
    score
end

function squarefrommidpoint(midpoint, width)
    squarefrommidpoint(midpoint, width...)
end

function squarefrommidpoint(midpoint, xwidth, ywidth)
    wx = xwidth ÷ 2
    wy = ywidth ÷ 2
    [(midpoint[1]-wx):(midpoint[1]+wx), (midpoint[2]-wy):(midpoint[2]+wy)]
end

"""
    slidingwindow()

Compares `reference` to sliding windows of equal size of `target`
across a square of width `searchwidth` using the metric `metric`.
The search window is centered at the center of reference, under the
assumption that the reference is a view of a larger matrix.
Returns the resulting potential surface.
"""
function slidingwindow(reference::SubArray, target::AbstractMatrix, metric::Function, searchwidth::Int)
    return slidingwindow(reference, target, reference.indices, metric, searchwidth)
end

function slidingwindow(reference::AbstractMatrix, target::AbstractMatrix,
                       indices, metric::Function, searchwidth::Int)
    window, shifts = shiftsinrange(reference, target, indices, searchwidth)
    distance = Matrix(undef, shifts[1], shifts[2])
    for i in 1:shifts[1], j in 1:shifts[2]
        # The .- 1 accounts for zero shift
        distance[i, j] = metric(reference, target[window[1] .+ i .- 1,
                                                  window[2] .+ j .- 1])
    end
    return distance
end

function slidingwindow!(score::AbstractMatrix, reference::AbstractMatrix, template::AbstractMatrix, metric::Function)
    # Slide the template across the entire reference
    iw, jw = size(reference) .- size(template)
    w1, w2 = size(template)
    for i in 1:iw
        @simd for j in 1:jw
            score[i, j] = metric(template, @view reference[i:w1+i-1, j:w2+j-1])
        end
    end
end

function shiftsinrange(template::SubArray, searchregion::AbstractMatrix, searchwidth::Int)
    shiftsinrange(template, searchregion, template.indices, searchwidth)
end

function shiftsinrange(template, searchregion::AbstractMatrix, indices, searchwidth::Int)
    @assert searchwidth < size(template, 1) "$searchwidth > $(size(template, 1))"
    @assert searchwidth < size(template, 2) "$searchwidth > $(size(template, 2))"
    # Handle boundary issues by clipping the search range to within [1, size target]
    start = clip.(first.(indices) .- searchwidth, size(searchregion))
    stop  = clip.(last.(indices)  .+ searchwidth, size(searchregion))
    # Construct the indices defining the initial window in the upper left
    # corner of the search box
    window = range.(start, start .+ size(template) .- 1, step=1)
    # The number of shifts of the window possible within the search range
    shifts = stop .- start .- size(template)
    return window, shifts
end

# function featureoverlap_perturbed(feature, target, searchwidth, indices)
    # x, y = indices
    # widthx, widthy = size(feature)
    # startx  = clip(first(x) - searchwidth, size(target, 1))
    # stopx   = clip(last(x)  + searchwidth, size(target, 1))
    # starty  = clip(first(y) - searchwidth, size(target, 2))
    # stopy   = clip(last(y)  + searchwidth, size(target, 2))
    # windowx = startx:(startx + widthx - 1)
    # windowy = starty:(starty + widthy - 1)
    # lengthx, lengthy = (stopx - startx - widthx,
                        # stopy - starty - widthy)
    # overlap = zeros((lengthx, lengthy))
    # for j in 1:lengthy, i in 1:lengthx
        # shiftx = i - 1  # Account for zero shift
        # shifty = j - 1
        # diff = (feature .- target[windowx .+ shiftx,
                                  # windowy .+ shifty]).^2 ./ feature
        # overlap[i, j] = sum(diff)
    # end
    # #return overlap

    # # Find the best overlap
    # indices = argmin(overlap)
    # x = range(startx + indices[1] + 1, length=widthx)
    # y = range(starty + indices[2] + 1, length=widthy)
    # [x, y]
# end

function perturbate(raw::AbstractArray)
    perturbed = similar(raw)
    for i in eachindex(raw)
        mu = raw[i]# < 1 ? 1 : mat[i]
        perturbed[i] = rand(Poisson(mu))
    end
    perturbed
end

function perturbate!(perturbed::AbstractArray, raw::AbstractArray)
    for i in eachindex(raw)
        mu = raw[i]# < 1 ? 1 : mat[i]
        #perturbed[i] = mu < 2 ? 0 : rand(Poisson(mu))
        perturbed[i] = rand(Poisson(mu))
    end
end

function arraybox(array, startx::T, stopx::T, starty::T, stopy::T) where T <: Integer
    subarray = view(array, startx:stopx, starty:stopy)
end

function splitregions2d(array, numberofregions; width=nothing, limit=50)
    edges = boxbounds(array, numberofregions, width=width)
    arrays = []
    for edge in edges
        subarray = arraybox(array, edge)
        if sum(subarray) > limit
            push!(arrays, subarray)
        end
    end
    arrays
end

function arraybox(array, edges)
    arraybox(array, edges[1]..., edges[2]...)
end

function boxbounds(array, numberofregions; width=nothing)
    sizex, sizey = size(array)
    stepx, stepy = floor.(Int, size(array)./numberofregions)
    if isnothing(width)
        width = stepx
    end
    edges = []
    for x in 1:numberofregions
        startx = (x-1)stepx + 1
        stopx = x*stepx + 1
        lowerx, upperx = boxlim(startx, stopx, width) .|> x -> clip(x, sizex)
        for y in 1:numberofregions
            starty = (y-1)stepy + 1
            stopy = y*stepy + 1
            lowery, uppery = boxlim(starty, stopy, width) .|> y -> clip(y, sizey)
            push!(edges, ((lowerx, upperx), (lowery, uppery)))
        end
    end
    edges
end

middle(start, stop) = (stop-start)/2 + start |> x -> round(Int64, x)
function boxlim(start, stop, width)
    lower = middle(start, stop) - width/2
    upper = middle(start, stop) + width/2
    (lower, upper) .|> x -> floor(Int, x)
end

function downsample(mat::AbstractMatrix{T}, blocksize=2) where T
    xlen, ylen = size(mat)
    newx, newy = xlen÷blocksize, ylen÷blocksize
    @assert(xlen%blocksize == 0, "Dimension 1 is not divisible by blocksize, $xlen")
    @assert(ylen%blocksize == 0, "Dimension 2 is not divisible by blocksize, $ylen")
    new = Matrix{Float64}(undef, newx, newy)
    s = 0
    for i in 1:blocksize:xlen, j in 1:blocksize:ylen
        s = 0
        for i_ in i:(i+blocksize-1), j_ in j:(j+blocksize-1)
            s += mat[i_, j_]
        end
        new[i÷blocksize + 1, j÷blocksize + 1] = 1/blocksize*s
    end
    return new
end

function neighbourreduce(mat::AbstractMatrix{T}; width=3, threshold=10) where T
    @assert isodd(width)
    xlen, ylen = size(mat)
    new = Matrix{Float64}(undef, ylen, xlen)
    s = 0
    w = width÷2
    for i in 1:xlen, j in 1:ylen
        s = 0
        for i_ in (i-w):(i+w), j_ in (j-w):(j+w)
            if 0 < i_ < xlen && 0 < j_ < ylen
                s += mat[i_, j_]
            end
        end
        if s ≥ threshold
            new[i, j] = mat[i, j]
        else
            new[i, j] = 0.0
        end
    end
    return new
end

function annotate_plot_2d(mat, x, y, w=nothing; ax=nothing)
    if ax ≡ nothing
        fig, ax = plt.subplots()
    end

    ax.matshow(transpose(mat) .|> log10)
    ax.scatter(x, y, marker="x", c="r")
    if w !== nothing
        for (i, w_) in enumerate(w)
            ax.annotate(string(w_), (x[i], y[i]))
        end
    end
    ax.invert_yaxis()
end

function compare_2d(x, y;kernel=Kernel.gaussian(3), metric::Function=L2)
    x = imfilter(x, kernel)
    y = imfilter(y, kernel)
    # Normalize y to x
    # y .= y.*sum(x)/sum(y)
    # Normalize each
    x ./= sum(x)
    y ./= sum(y)
    return metric(x, y)
    #return (x .- y).^2
end

function compare_methods(reference, target; neighbourwidth=9, neighbourthreshold=8)
    edges = histogram_edges(reference[!, :e], reference[!, :Δe], nbins=1000)
    refmat_raw = JSort.histogram(reference[!, :e], reference[!, :Δe], edges...)
    tarmat_raw = JSort.histogram(target[!, :e], target[!, :Δe], edges...)

    Reds = alphamap("Reds")
    Blues = alphamap("Blues")

    fig, ax = plt.subplots()
    ax.set_title("Original")
    ax.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                  cmap=Reds, alpha=0.7)
    ax.pcolormesh(tarmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                  cmap=Blues, alpha=0.7)
    ax.set_ylim([0, 600])
    ax.set_xlim([0, 800])

    refmat = neighbourreduce(refmat_raw, width=neighbourwidth, threshold=neighbourthreshold)
    tarmat = neighbourreduce(tarmat_raw, width=neighbourwidth, threshold=neighbourthreshold)
    points = featurealign2d_weighted(refmat, tarmat, numregions=15, width=nothing,
                                 limit=1500, searchwidth=40, weight=ThresholdSum(.2))
    rpx, tpx, rpy, tpy, weights = points
    annotate_plot_2d(tarmat, tpx, tpy, weights)
    edgesx, edgesy = edges
    rpx, tpx = edgesx[rpx], edgesx[tpx]
    rpy, tpy = edgesy[rpy], edgesy[tpy]

    lm = pyimport("sklearn.linear_model")

    models = [leastsquares,
              (x, y) -> begin
                lr = lm.TheilSenRegressor(fit_intercept=true)
                X = reshape(x, :, 1)
                lr.fit(X, y)
                return [lr.intercept_; lr.coef_]
              end,
              (x, y) -> begin
                  lr = lm.RANSACRegressor()
                  X = reshape(x, :, 1)
                  lr.fit(X, y)
                  return [lr.estimator_.intercept_; lr.estimator_.coef_]
              end]
    names = ["OLS", "TheilSen", "RANSAC"]
    metric = (x, y) -> splitmetric(x, y, metric=L2)

    for (model, name) in zip(models, names)
        coeffx = model(tpx, rpx)
        coeffy = model(tpy, rpy)

        px = Poly(coeffx)
        py = Poly(coeffy)

        e  = px.(target[!, :e])
        de = py.(target[!, :Δe])

        calibmat = JSort.histogram(e, de, edges...);
        err1 = compare_2d(refmat, calibmat, metric=metric)
        err2 = compare_2d(refmat, tarmat, metric=metric)
        println("$name - Error ratio: $(err1/err2)")

        fig, ax = plt.subplots(ncols=2, nrows=2)
        plotresiduals(tpx, rpx, coeffx, ax=ax[1:2])
        plotresiduals(tpy, rpy, coeffy, ax=ax[3:4])
        fig.suptitle(name)

        fig, ax = plt.subplots(figsize=(15, 15))
        ax.set_title("Aligned")
        ax.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Reds, alpha=0.7)
        ax.pcolormesh(calibmat', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Blues, alpha=0.7)
        axins = ax.inset_axes([0.5, 0.5, 0.47, 0.47])
        axins.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Reds, alpha=0.7)
        axins.pcolormesh(calibmat', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Blues, alpha=0.7)
        axins.set_xlim(550, 800)
        axins.set_ylim(80, 220)
        ax.indicate_inset_zoom(axins)
        axins = ax.inset_axes([0.0, 0.5, 0.47, 0.47])
        axins.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Reds, alpha=0.7)
        axins.pcolormesh(calibmat', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Blues, alpha=0.7)
        axins.set_xlim(0, 200)
        axins.set_ylim(400, 600)
        ax.indicate_inset_zoom(axins)
    end

    loss = ["absolute_loss", "squared_loss"]
    models = [(x, y, w) -> begin
                  X = reshape(x, :, 1)
                  lr.fit(X, y, w)
                  return [lr.estimator_.intercept_; lr.estimator_.coef_]
              end
                for lr in
              [lm.RANSACRegressor(loss=l) for l in loss]]
    names = ["RANSAC "*l for l in loss]

    # Weighted
    for (name, model) in zip(names, models)
        coeffx = model(tpx, rpx, weights)
        coeffy = model(tpy, rpy, weights)

        px = Poly(coeffx)
        py = Poly(coeffy)

        e  = px.(target[!, :e])
        de = py.(target[!, :Δe])

        calibmat = JSort.histogram(e, de, edges...);
        err1 = compare_2d(refmat, calibmat, metric=metric)
        err2 = compare_2d(refmat, tarmat, metric=metric)
        println("$name - Error ratio: $(err1/err2)")

        fig, ax = plt.subplots(ncols=2, nrows=2)
        plotresiduals(tpx, rpx, coeffx, ax=ax[1:2])
        plotresiduals(tpy, rpy, coeffy, ax=ax[3:4])
        fig.suptitle(name)

        fig, ax = plt.subplots(figsize=(15, 15))
        ax.set_title("Aligned")
        ax.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Reds, alpha=0.7)
        ax.pcolormesh(calibmat', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Blues, alpha=0.7)
        axins = ax.inset_axes([0.5, 0.5, 0.47, 0.47])
        axins.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Reds, alpha=0.7)
        axins.pcolormesh(calibmat', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Blues, alpha=0.7)
        axins.set_xlim(550, 800)
        axins.set_ylim(80, 220)
        ax.indicate_inset_zoom(axins)
        axins = ax.inset_axes([0.0, 0.5, 0.47, 0.47])
        axins.pcolormesh(refmat_raw', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Reds, alpha=0.7)
        axins.pcolormesh(calibmat', norm=PyPlot.matplotlib.colors.LogNorm(),
                      cmap=Blues, alpha=0.7)
        axins.set_xlim(0, 200)
        axins.set_ylim(400, 600)
        ax.indicate_inset_zoom(axins)
    end
    println("bump")
end

function alphamap(name::String)
PyCall.py"""
import numpy as np
import matplotlib.pylab as pl
from matplotlib.colors import ListedColormap

def alpha_map(name):
    # Choose colormap
    cmap = pl.cm.get_cmap(name)

    # Get the colormap colors
    my_cmap = cmap(np.arange(cmap.N))

    # Set alpha
    my_cmap[:,-1] = np.linspace(0, 1, cmap.N)

    # Create new colormap
    my_cmap = ListedColormap(my_cmap)
    return my_cmap
"""
return PyCall.py"alpha_map"(name)
end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
##############################################################################
##
## Hcat specialization
##
##############################################################################

# hcat! for 2 arguments, only a vector or a data set is allowed

# Modify Dataset
function hcat!(ds1::Dataset, ds2::AbstractDataset;
               makeunique::Bool=false)
    u = add_names(index(ds1), index(ds2), makeunique=makeunique)
    for i in 1:length(u)
        ds1[!, u[i]] = ds2[:, i]
    end
    _modified(_attributes(ds1))
    return ds1
end

# definition required to avoid hcat! ambiguity

# Modify Dataset
hcat!(ds1::Dataset, ds2::Dataset;
      makeunique::Bool=false) =
    invoke(hcat!, Tuple{Dataset, AbstractDataset}, ds1, ds2,
           makeunique=makeunique)::Dataset

# Modify Dataset
# hcat!(ds::Dataset, x::AbstractVector; makeunique::Bool=false, copycols::Bool=true) =
#     hcat!(ds, Dataset(AbstractVector[x], [:x1], copycols=copycols),
#           makeunique=makeunique, copycols=copycols)

# Modify Dataset
# hcat!(x::AbstractVector, ds::Dataset; makeunique::Bool=false, copycols::Bool=true) =
    # hcat!(Dataset(AbstractVector[x], [:x1], copycols=copycols), ds,
          # makeunique=makeunique, copycols=copycols)

# Modify Dataset
hcat!(x, ds::Dataset; makeunique::Bool=false) =
    throw(ArgumentError("x must be AbstractDataset"))

# Modify Dataset
hcat!(ds::Dataset, x; makeunique::Bool=false) =
    throw(ArgumentError("x must be AbstractVector or AbstractDataset"))

# hcat! for 1-n arguments

# Modify Dataset
hcat!(ds::Dataset; makeunique::Bool=false) = ds

# Modify Dataset
# hcat!(a::Dataset, b, c...; makeunique::Bool=false, copycols::Bool=true) =
#     hcat!(hcat!(a, b, makeunique=makeunique, copycols=copycols),
#           c..., makeunique=makeunique, copycols=copycols)

# hcat

# # Create Dataset
# Base.hcat(ds::Dataset, x; makeunique::Bool=false, copycols::Bool=true) =
#     hcat!(copy(ds, copycols=copycols), x,
#           makeunique=makeunique, copycols=copycols)

# Create Dataset
Base.hcat(ds1::Dataset, ds2::AbstractDataset;
          makeunique::Bool=false) =
    hcat!(copy(ds1), ds2,
          makeunique=makeunique)

# Create Dataset
Base.hcat(ds1::Dataset, ds2::AbstractDataset, dsn::AbstractDataset...;
          makeunique::Bool=false) =
    hcat!(hcat(ds1, ds2, makeunique=makeunique), dsn...,
          makeunique=makeunique)


text:
# Plot tools using Gadfly library for Julia
# Berkeley 2015 - Maxime Legrand

#Pkg.add("Gadfly")
#Pkg.add("DataFrames")
#Pkg.add("Color")

#using Gadfly

"""Matrix plot
Input : M matrix, zeros boolean
Output : plotting of M with colored squares, zeros can be hidden by setting zeros=false"""
function matplot(M;zeros=true, title="",
    col1=Color.RGB(234/255,234/255,246/255), col2=Color.RGB(105/255,165/255,255/255))
    #col1=color("lavender"), col2=color("blue"))
    I,J = size(M)
    if !zeros
        is, js, values = findnz(M)
    else
        is, js, values = repeat([1:I], outer=[J]), repeat([1:J], inner=[I]), vec(M)
    end
    df = DataFrames.DataFrame(i=is-0.5, j=js-0.5, value=values)
    plot(df, x="j", y="i", color="value",
        Coord.cartesian(yflip=true, fixed=true),
        Scale.ContinuousColorScale(Scale.lab_gradient(col1, col2)),
        Scale.x_continuous, Scale.y_continuous,
        Guide.xticks(ticks=[0,J]), Guide.yticks(ticks=[0,I]),
        Guide.title(title),
        Theme(grid_line_width=0mm),
        Geom.rectbin)
end


text:
export
		get_matrix,
		get_outcome_names,
		set_definition,
		set_number_of_outcomes,
		is_nodetype_id
export DSL_DECISION, DSL_CHANCE, DSL_DETERMINISTIC, DSL_UTILITY, DSL_DISCRETE, DSL_CASTLOGIC,
	   DSL_DEMORGANLOGIC, DSL_NOISYMAXLOGIC, DSL_NOISYADDERLOGIC, DSL_PARENTSCONTIN,
	   DSL_SCC, DSL_DCHILDHPARENT, DSL_CONTINUOUS, DSL_TRUTHTABLE, DSL_CPT, DSL_NOISY_MAX,
	   DSL_NOISY_ADDER, DSL_CAST, DSL_DEMORGAN, DSL_LIST, DSL_TABLE, DSL_MAU, DSL_EQUATION,
	   DSL_EQUATION_SCC, DSL_DCHILD_HPARENT, DSL_DISTRIBUTION, DSL_HEQUATION, DSL_NO_DEFINITION

# -------------

const DSL_DECISION         = unsafe_load(cglobal((:EDSL_DECISION,        LIB_SMILE), Int32))
const DSL_CHANCE           = unsafe_load(cglobal((:EDSL_CHANCE,          LIB_SMILE), Int32))
const DSL_DETERMINISTIC    = unsafe_load(cglobal((:EDSL_DETERMINISTIC,   LIB_SMILE), Int32))
const DSL_UTILITY          = unsafe_load(cglobal((:EDSL_UTILITY,         LIB_SMILE), Int32))
const DSL_DISCRETE         = unsafe_load(cglobal((:EDSL_DISCRETE,        LIB_SMILE), Int32))
const DSL_CASTLOGIC        = unsafe_load(cglobal((:EDSL_CASTLOGIC,       LIB_SMILE), Int32))
const DSL_DEMORGANLOGIC    = unsafe_load(cglobal((:EDSL_DEMORGANLOGIC,   LIB_SMILE), Int32))
const DSL_NOISYMAXLOGIC    = unsafe_load(cglobal((:EDSL_NOISYMAXLOGIC,   LIB_SMILE), Int32))
const DSL_NOISYADDERLOGIC  = unsafe_load(cglobal((:EDSL_NOISYADDERLOGIC, LIB_SMILE), Int32))
const DSL_PARENTSCONTIN    = unsafe_load(cglobal((:EDSL_PARENTSCONTIN,   LIB_SMILE), Int32))
const DSL_SCC              = unsafe_load(cglobal((:EDSL_SCC,             LIB_SMILE), Int32))
const DSL_DCHILDHPARENT	   = unsafe_load(cglobal((:EDSL_DCHILDHPARENT,   LIB_SMILE), Int32))
const DSL_CONTINUOUS	   = unsafe_load(cglobal((:EDSL_CONTINUOUS,      LIB_SMILE), Int32))

const DSL_TRUTHTABLE       = unsafe_load(cglobal((:EDSL_TRUTHTABLE,      LIB_SMILE), Int32))
const DSL_CPT              = unsafe_load(cglobal((:EDSL_CPT,             LIB_SMILE), Int32))
const DSL_NOISY_MAX        = unsafe_load(cglobal((:EDSL_NOISY_MAX,       LIB_SMILE), Int32))
const DSL_NOISY_ADDER      = unsafe_load(cglobal((:EDSL_NOISY_ADDER,     LIB_SMILE), Int32))
const DSL_CAST             = unsafe_load(cglobal((:EDSL_CAST,            LIB_SMILE), Int32))
const DSL_DEMORGAN         = unsafe_load(cglobal((:EDSL_DEMORGAN,        LIB_SMILE), Int32))
const DSL_LIST             = unsafe_load(cglobal((:EDSL_LIST,            LIB_SMILE), Int32))
const DSL_TABLE            = unsafe_load(cglobal((:EDSL_TABLE,           LIB_SMILE), Int32))
const DSL_MAU              = unsafe_load(cglobal((:EDSL_MAU,             LIB_SMILE), Int32))
const DSL_EQUATION         = unsafe_load(cglobal((:EDSL_EQUATION,        LIB_SMILE), Int32))
const DSL_EQUATION_SCC     = unsafe_load(cglobal((:EDSL_EQUATION_SCC,    LIB_SMILE), Int32))
const DSL_DCHILD_HPARENT   = unsafe_load(cglobal((:EDSL_DCHILD_HPARENT,  LIB_SMILE), Int32))
const DSL_DISTRIBUTION     = unsafe_load(cglobal((:EDSL_DISTRIBUTION,    LIB_SMILE), Int32))
const DSL_HEQUATION        = unsafe_load(cglobal((:EDSL_HEQUATION,       LIB_SMILE), Int32))
const DSL_NO_DEFINITION    = unsafe_load(cglobal((:EDSL_NO_DEFINITION,   LIB_SMILE), Int32))

# -------------

function get_matrix( nodedef::NodeDefinition )
	ptr = ccall( (:nodedef_GetMatrix, LIB_SMILE), Ptr{Void}, (Ptr{Void},), nodedef.ptr )
	return DMatrix(ptr)
end

function get_outcome_names(nodedef::NodeDefinition)
	ptr = ccall( (:nodedef_GetOutcomeNames, LIB_SMILE), Ptr{Void}, (Ptr{Void},), nodedef.ptr )
	return IdArray(ptr)
end

get_type(nodedef::NodeDefinition) =
	ccall( (:nodedef_GetType, LIB_SMILE), Int32, (Ptr{Void},), nodedef.ptr)

set_definition( nodedef::NodeDefinition, dblarr::DoubleArray ) =
	ccall( (:nodedef_SetDefinition, LIB_SMILE), Int32, (Ptr{Void},Ptr{Void}), nodedef.ptr, dblarr.ptr )

set_number_of_outcomes( nodedef::NodeDefinition, idarr::IdArray ) =
	ccall( (:nodedef_SetNumberOfOutcomes_Names, LIB_SMILE), Int32, (Ptr{Void},Ptr{Void}), nodedef.ptr, idarr.ptr )
set_number_of_outcomes( nodedef::NodeDefinition, outcomeNumber::Integer ) =
	ccall( (:nodedef_SetNumberOfOutcomes_Int, LIB_SMILE), Int32, (Ptr{Void},Cint), nodedef.ptr, outcomeNumber )
get_number_of_outcomes(nodedef::NodeDefinition) =
	ccall( (:nodedef_GetNumberOfOutcomes, LIB_SMILE), Int32, (Ptr{Void},), nodedef.ptr)

function is_nodetype_id( typeID::Integer )
	return typeID == DSL_CPT  || typeID == DSL_TRUTHTABLE ||
	       typeID == DSL_LIST || typeID == DSL_NOISY_MAX  ||
	       typeID == DSL_MAU  || typeID == DSL_TABLE
end


text:
function Assembletet(IK,JK,VK,ke,i,j,k,l)

    DOF = [3*i-2 3*i-1 3*i 3*j-2 3*j-1 3*j 3*k-2 3*k-1 3*k 3*l-2 3*l-1 3*l]
    for n1 = 1:12
        for n2 = 1:12
            #KE[DOF[n1],DOF[n2]] = KE[DOF[n1],DOF[n2]] + ke[n1,n2]
            IK = append!(IK,DOF[n1])
            JK = append!(JK,DOF[n2])
            VK = append!(VK,ke[n1,n2])
        end
    end

    return IK, JK, VK
end

max_stars_repo_path: src/single-core juSFEM/assembletet.jl

max_stars_repo_name: Elliothuo/juSFEM

max_issues_repo_path: src/single-core juSFEM/assembletet.jl

max_issues_repo_name: Elliothuo/juSFEM

max_forks_repo_path: src/single-core juSFEM/assembletet.jl

max_forks_repo_name: Elliothuo/juSFEM


text:
using FCBClassification
using Base.Test

@test false


text:
"""
    @mutate_mtg!(node, args...,kwargs...)

Mutate the mtg nodes in place.

# Arguments

- `mtg`: the mtg to mutate
- `args...`: The computations to apply to the nodes (see examples)
- `kwargs...`: Optional keyword arguments for traversing and filtering (see details)

# Details

As for [`descendants`](@ref) and [`ancestors`](@ref), kwargs can be any filter from:

- `scale = nothing`: The scale to filter-in (i.e. to keep). Usually a Tuple-alike of integers.
- `symbol = nothing`: The symbol to filter-in. Usually a Tuple-alike of Strings.
- `link = nothing`: The link with the previous node to filter-in. Usually a Tuple-alike of Char.
- `all::Bool = true`: Return all filtered-in nodes (`true`), or stop at the first node that
is filtered out (`false`).
- `filter_fun = nothing`: Any filtering function taking a node as input, e.g. [`isleaf`](@ref).
- `traversal`: The type of tree traversal. By default it is using `AbstractTrees.PreOrderDFS`.

# Examples

```julia
# Importing an mtg from the package:
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)

# Compute a new attribute with the scales and add 2 to its values:
@mutate_mtg!(mtg, scaling = node.scales .+ 2)

# Compute several new attributes, some based on others:
@mutate_mtg!(mtg, x = length(node.name), y = node.x + 2, z = sum(node.y))

# We can also use it without parenthesis:

@mutate_mtg! mtg x = length(node.name)
```
"""
macro mutate_mtg!(mtg, args...)
    arguments = (args...,)
    # Get the value of the filters if any:
    flt, kwargs, args = MultiScaleTreeGraph.parse_macro_args(arguments)

    expr = quote
        # Check the filters consistency with the mtg:
        check_filters($(mtg), scale = $(flt.scale), symbol = $(flt.symbol), link = $(flt.link))
        # Traverse the mtg:
        traversed_mtg = $(kwargs[:traversal])($(mtg))
        for i00000000 in traversed_mtg
            # NB: using i00000000 to avoid naming clash with variables when using rewrite_expr!(i,...)
            if is_filtered(i00000000, $(flt.scale), $(flt.symbol), $(flt.link), $(flt.filter_fun))
                @mutate_node!(i00000000, $(args...))
            elseif !$(kwargs[:all])
                # In this case (all == false) we stop as soon as we reached the
                # first filtered-out value. The default behavior is defined in parse_macro_args
                # and is equal to true.
                break
            end
        end
    end
    esc(expr)
end

"""
    @mutate_node!(node, args...)

Mutate a single node in place.

# Arguments

- `node`: the node to mutate
- `args...`: The computations to apply to the node (see examples)

# See also

[`@mutate_mtg!`](@ref) to mutate all nodes of an mtg.

# Examples

```julia
# Importing an mtg from the package:
file = joinpath(dirname(dirname(pathof(MultiScaleTreeGraph))),"test","files","simple_plant.mtg")
mtg = read_mtg(file)

# Compute a new attribute with the scales and add 2 to its values:
@mutate_node!(mtg, scaling = node.scales .+ 2)

# The computation is only applied to the root node. To apply it to all nodes,
# see @mutate_mtg!

# Compute several new attributes, some based on others:
@mutate_node!(mtg, x = length(node.name), y = node.x + 2, z = sum(node.y))

# We can also use it without parenthesis:

@mutate_node! mtg x = length(node.name)
```
"""
macro mutate_node!(node, args...)
    arguments = (args...,)
    rewrite_expr!(:($node), arguments)
    expr = quote $(arguments...); nothing end
    esc(expr)
end

"""
    rewrite_expr!(arguments)

Re-write the call to the variables of a node in an expression to match their location: leave
it as it is if the variable is a node field, or add `attributes` after the node if it is
an attribute.

# Examples

```
test = :(x = node.name)
MultiScaleTreeGraph.rewrite_expr!(:mtg,test)
test
# :(mtg.attributes[:x] = mtg.name)

test = :(x = node.foo)
MultiScaleTreeGraph.rewrite_expr!(:mtg,test)
test
# :(mtg.attributes[:x] = mtg.attributes[:foo])

test = :(x = node.MTG.symbol)
MultiScaleTreeGraph.rewrite_expr!(:mtg,test)
test
# :(mtg.attributes[:x] = mtg.MTG.symbol)
```
"""
function rewrite_expr!(node_name, arguments::Expr)

    # For the Left-Hand Side (LHS) when using the form `x = RHS`. The LHS is treated
    # differently here because `x` must be treated as a field or an attribute of the node,
    # and never as an object from another module.
    # On the contrary the RHS can use variables from elsewhere, e.g. a constant defined in
    # the REPL
    if isa(arguments, Expr) && arguments.head == :(=) && isa(arguments.args[1], Symbol)
        arguments.args[1] = :($(node_name).attributes[$(QuoteNode(arguments.args[1]))])
        # if !(Symbol(replace(arg,"node."=>"")) in fieldnames(Node))
        # x.args[1] = :(node.attributes)
    end

    # For the RHS:
    for x in arguments.args
        arg = string(x)
        if isa(x, Expr) &&
            (x.head == :. || x.head == :ref) &&
            occursin(r"^node", arg) &&
            !occursin(string(node_name), arg)
            # x here is defined either as node.variable or node[:variable], we must replace
            # by node_name[:variable]
            if !(Symbol(replace(arg, "node." => "")) in fieldnames(Node))
                if any(match.(Regex.("\\." .* string.(fieldnames(Node))), arg) .!= nothing)
                    # If the expression contains node attributes, only replace the node name
                    # because they are reserved keywords, e.g.: node.MTG.scale becomes
                    # node_name.MTG.scale
                    x.args[1].args[1] = :($(node_name))
                else
                    x.args[1] = :($(node_name))
                    x.head = :ref
                end
            else
                x.args[1] = :($(node_name))
            end
        elseif isa(x, Expr) && x.head == :call && occursin("node", arg)
            # Call to a function, and we pass node as argument
            for i in 1:length(x.args)
                # The node is given as is to the function, e.g. fn(node):
                x.args[i] == :node ? x.args[i] = :($(node_name)) : nothing
            end
            rewrite_expr!(node_name, x)
        else
            rewrite_expr!(node_name, x)
        end
    end
end

function rewrite_expr!(node_name, arguments)
    nothing
end

function rewrite_expr!(node_name, arguments::Tuple)
    for x in arguments
        rewrite_expr!(node_name, x)
    end
end


text:
abstract type AbstractPeriod end

function checkInstantaneous(x::T) where {T <: TimePeriod}
    @assert x ≥ zero(T) "period start cannot be negative"
    x
end

checkInstantaneous(x) = x
struct Instantaneous{T} <: AbstractPeriod
    data::T
    Instantaneous{T}(x) where {T} = new(checkInstantaneous(x))
end
Instantaneous(x::T) where {T} = Instantaneous{T}(x)

Instantaneous() = Instantaneous(Millisecond(0))

function checkProlonged(x::I)  where {T <: TimePeriod, I <: AbstractInterval{T}}
    @assert leftendpoint(x) ≥ zero(T) "period start cannot be negative"
    @assert width(x) > zero(T) "period duration must be larger than zero"
    x
end

checkProlonged(x) = x

struct Prolonged{T} <: AbstractPeriod
    data::T
    Prolonged{T}(x) where {T} = new(checkProlonged(x))
end
Prolonged(x::T) where {T} = Prolonged{T}(x)

Prolonged() = Prolonged(Millisecond(0)..Millisecond(1))

AbstractPeriod(start::T) where {T <: TimePeriod} = Instantaneous(start)
AbstractPeriod(x::I) where {T <: TimePeriod, I <: AbstractInterval{T}} = Prolonged(x)
AbstractPeriod(start::T1, stop::T2) where {T1 <: TimePeriod, T2 <: TimePeriod} = Prolonged(start..stop)
AbstractPeriod(start::T, _::Nothing) where {T <: TimePeriod} = Instantaneous(start)

start(x::Instantaneous) = x.data
start(x::Prolonged) = leftendpoint(x.data)
duration(x::Instantaneous{T}) where {T <: TimePeriod} = zero(T)
duration(x::Prolonged) = width(x.data)

struct Temporal{V <: AbstractTimeLine, T <: AbstractPeriod}
    video::V
    time::T
    comment::String

    function Temporal{V, T}(video, time, comment) where {V <: AbstractTimeLine, T <: AbstractPeriod}
        @assert stop(time) ≤ duration(video) "POIs are outside the video's timeline (stop: $(stop(time)); duration: $(duration(video)); file: $(files(video))"
        new(video, time, comment)
    end

end
Temporal(video::V, time::T, comment::String) where {V <: AbstractTimeLine, T <: AbstractPeriod} = Temporal{V, T}(video, time, comment)

Temporal() = Temporal(WholeVideo(), Instantaneous(), "")


text:
abstract type Strategy end

"""
    getmove(board, s::Strategy)
"""
function getmove end

"""
    playgame(g::Game, s1::Strategy, s2::Strategy)
    playgame(board, s1::Strategy, s2::Strategy)
"""

#TODO: decide output
function playgame(board, s::Vararg{Strategy})

    while !isempty(getmoves(board))
        p = playerturn(board)
        m = getmove(board, s[p])
        board = play(board, m)
    end

    return board#TODO
end

playgame(g::Game, s::Vararg{Strategy}) = playgame(initialboard(g), s...)

"""
    name(s::Strategy)
"""
function name end

"""
    getvalues(strategy)
"""
function getvalues end

"""
    setvalue!(strategy, name, value)
"""
function setvalue! end

"""
    getvarsnames(strategy)
"""
function getvarsnames end


text:
# !!! this script needs Julia v 0.4 or greater !!!
# constructing:
d = Date(2014,9,1) #=> 2014-09-01
dt = DateTime(2014,9,1,12,30,59,1) #=> 2014-09-01T12:30:59.001
# accessors:
Dates.year(d)
Dates.month(d)
Dates.week(d)
Dates.day(d)
# functions:
Dates.isleapyear(d)
Dates.dayofyear(d)
Dates.monthname(d)
Dates.daysinmonth(d)

max_stars_repo_path: Module 1/Chapter02/Old/dates.jl

max_issues_repo_path: Module 1/Chapter02/Old/dates.jl

max_forks_repo_path: Module 1/Chapter02/Old/dates.jl


text:
using CompatHelper
using Documenter

DocMeta.setdocmeta!(CompatHelper, :DocTestSetup, :(using CompatHelper); recursive=true)

makedocs(;
    modules=[CompatHelper],
    authors="Dilum Aluthge, Brown Center for Biomedical Informatics, and contributors",
    repo="https://github.com/JuliaRegistries/CompatHelper.jl/blob/{commit}{path}#{line}",
    sitename="CompatHelper.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://JuliaRegistries.github.io/CompatHelper.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Configuration Options" => "options.md",
        "Acknowledgements" => "acknowledgements.md",
        "Other Environments" => "other-environments.md",
        "TroubleShooting" => "troubleshooting.md",
    ],
    strict=true,
)

deploydocs(; repo="github.com/JuliaRegistries/CompatHelper.jl")


text:
# Autogenerated wrapper script for libcgal2_julia_jll for x86_64-apple-darwin
export libcgal_julia_exact, libcgal_julia_inexact

using CGAL2_jll
using libcxxwrap_julia_jll
JLLWrappers.@generate_wrapper_header("libcgal2_julia")
JLLWrappers.@declare_library_product(libcgal_julia_exact, "@rpath/libcgal_julia_exact.0.dylib")
JLLWrappers.@declare_library_product(libcgal_julia_inexact, "@rpath/libcgal_julia_inexact.0.dylib")
function __init__()
    JLLWrappers.@generate_init_header(CGAL2_jll, libcxxwrap_julia_jll)
    JLLWrappers.@init_library_product(
        libcgal_julia_exact,
        "lib/libcgal_julia_exact.0.18.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libcgal_julia_inexact,
        "lib/libcgal_julia_inexact.0.18.0.dylib",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test
using StaticArrays
using ClimateMachine.VariableTemplates

@testset "Test complex models" begin

    include("complex_models.jl")

    FT = Float32

    # test getproperty
    m = ScalarModel()
    st = state(m, FT)
    vs = varsize(st)
    a_global = collect(1:vs)
    v = Vars{st}(a_global)
    @test v.x == FT(1)

    Nv = 4
    m = VectorModel{Nv}()
    st = state(m, FT)
    vs = varsize(st)
    a_global = collect(1:vs)
    v = Vars{st}(a_global)
    @test v.x == SVector{Nv, FT}(1:Nv)

    N = 3
    M = 6
    m = MatrixModel{N, M}()
    st = state(m, FT)
    vs = varsize(st)
    a_global = collect(1:vs)
    v = Vars{st}(a_global)
    @test v.x == SHermitianCompact{N, FT, M}(collect(1:(1 + M - 1)))

    Nv = 3
    N = 3
    M = 6
    m = CompositModel(Nv, N, M)
    st = state(m, FT)
    vs = varsize(st)
    a_global = collect(1:vs)
    v = Vars{st}(a_global)

    scalar_model = v.scalar_model
    @test v.scalar_model.x == FT(1)

    vector_model = v.vector_model
    @test v.vector_model.x == SVector{Nv, FT}([2, 3, 4])

    matrix_model = v.matrix_model
    @test v.matrix_model.x ==
          SHermitianCompact{N, FT, M}(collect(5:(5 + M - 1)))

    Nv = 3
    N = 5
    m = NTupleContainingModel(N, Nv)
    st = state(m, FT)
    vs = varsize(st)
    a_global = collect(1:vs)
    v = Vars{st}(a_global)

    offset = (Nv + 1) * N
    @test v.vector_model.x == SVector{Nv, FT}(collect(1:Nv) .+ offset)
    @test v.scalar_model.x == FT(1 + Nv) + offset

    unval(::Val{i}) where {i} = i
    @unroll_map(N) do i
        @test m.ntuple_model[i] isa NTupleModel
        @test m.ntuple_model[i].scalar_model isa ScalarModel
        @test v.ntuple_model[i].scalar_model.x ==
              FT(unval(i)) + (Nv) * (unval(i) - 1)
        @test v.vector_model.x == SVector{Nv, FT}(1:Nv) .+ offset
        @test v.scalar_model.x == FT(Nv + 1) + offset
    end

    # test flattenednames
    fn = flattenednames(st)
    j = 1
    for i in 1:N
        @test fn[j] === "ntuple_model[$i].scalar_model.x"
        j += 1
        for k in 1:Nv
            @test fn[j] === "ntuple_model[$i].vector_model.x[$k]"
            j += 1
        end
    end
    for k in 1:Nv
        @test fn[j] === "vector_model.x[$k]"
        j += 1
    end
    @test fn[j] === "scalar_model.x"

    # test flattened_tup_chain
    ftc = flattened_tup_chain(st)
    j = 1
    for i in 1:N
        @test ftc[j] === (:ntuple_model, i, :scalar_model, :x)
        j += 1
        @test ftc[j] === (:ntuple_model, i, :vector_model, :x)
        j += 1
    end
    @test ftc[j] === (:vector_model, :x)
    j += 1
    @test ftc[j] === (:scalar_model, :x)

    # test varsindex
    ntuple(N) do i
        i_val = Val(i)
        i_sm = varsindex(st, :ntuple_model, i_val, :scalar_model, :x)
        i_vm = varsindex(st, :ntuple_model, i_val, :vector_model, :x)
        nt_offset = (Nv + 1) - 1

        i_sm_correct = (i + nt_offset * (i - 1)):(i + nt_offset * (i - 1))
        @test i_sm == i_sm_correct

        offset = 1
        i_start = i + nt_offset * (i - 1) + offset
        i_vm_correct = (i_start):(i_start + Nv - 1)
        @test i_vm == i_vm_correct
    end

    convert_to_val(sym) = sym
    convert_to_val(i::Int) = Val(i)
    # test that getproperty matches varsindex
    ntuple(N) do i
        i_ϕ = varsindex(st, convert_to_val.(ftc[i])...)
        ϕ = getproperty(v, ftc[i])
        @test all(parent(v)[i_ϕ] .≈ ϕ)
    end

    # test getproperty with tup-chain
    @unroll_map(N) do i
        @test v.scalar_model.x == getproperty(v, (:scalar_model, :x))
        @test v.vector_model.x == getproperty(v, (:vector_model, :x))
        @test v.ntuple_model[i] == getproperty(v, (:ntuple_model, unval(i)))
        @test v.ntuple_model[i].scalar_model ==
              getproperty(v, (:ntuple_model, unval(i), :scalar_model))
        @test v.ntuple_model[i].scalar_model.x ==
              getproperty(v, (:ntuple_model, unval(i), :scalar_model, :x))
    end

end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
#==============================================================================
    Poincare plot example
        event based simulation of chaotic spiking neural network
        of quadratic integrate-and-fire neurons

    tested:     PASS (and it is fucking sweet)
    author:     discourse.julialang.org/u/RainerEngelken
 =============================================================================#

using   PyPlot
using   StaticArrays

function poincare()
    n,ϕ,𝚽=SMatrix{3,3}(0 .<[0 0 0;1 0 1;0 1 0]),randn(3),[]     # define adjacency matrix, initialize network
    for s=1:8^7                                                 # number of spikes in calculation
        m,j=findmax(ϕ)                                          # find next spiking neuron j
        ϕ.+=π/2-m                                               # evolve phases till next spike time
        ϕ[n[:,j]],ϕ[j]=atan.(tan.(ϕ[n[:,j]]).-1),-π/2           # update postsynaptic neurons, reset spiking neuron
        j==1 && append!(𝚽,ϕ[2:3])                               # save neuron 2 & 3 whenever neuron 1 spikes
    end
    plot(𝚽[1:2:end],𝚽[2:2:end],".k",ms=.1);axis("off")
end

max_stars_repo_name: mahiki/julia-is-great

max_issues_repo_name: mahiki/julia-is-great

max_forks_repo_name: mahiki/julia-is-great


text:
import PoseComposition: Pose, IDENTITY_POSE, IDENTITY_ORN, interp, ⊗, ⊘
import Rotations: QuatRotation, RotZYX
import StaticArrays: SVector, @SVector
import Test: @test, @testset

@testset "Pose group operations" begin
  p1 = Pose(1.0, 1.1, 1.2, QuatRotation(√(0.1), √(0.2), √(0.3), √(0.4)))
  p2 = Pose(2.0, 2.1, 2.2, QuatRotation(√(0.2), √(0.2), √(0.3), √(0.3)))
  p3 = Pose(3.0, 3.1, 3.2, QuatRotation(√(0.3), √(0.3), √(0.3), √(0.1)))

  function testPosesApprox(a, b; kwargs...)
    @test isapprox(a.pos, b.pos; kwargs...)
    @test isapprox(a.orientation, b.orientation; kwargs...)
  end

  testPosesApprox((p1 * p2) * p3, p1 * (p2 * p3))
  testPosesApprox((p1 * p2) / p2, p1)
  testPosesApprox((p1 / p2) * p2, p1)
  testPosesApprox(p1 * (p1 \ p2), p2)
  testPosesApprox(inv(p1) * p1, IDENTITY_POSE; atol=1e-14)
  testPosesApprox(p1 * inv(p1), IDENTITY_POSE; atol=1e-14)
end

@testset "Commutative diagram: a * b.pos == (a * b).pos" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  bpos = @SVector([7, 8, 9])
  b = Pose(bpos, RotZYX(1.0, 1.1, 1.2))
  @test a * bpos ≈ (a * b).pos
end

@testset "Associativity: (a * b) * cpos == a * (b * cpos)" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  b = Pose([7, 8, 9], RotZYX(1.0, 1.1, 1.2))
  cpos = @SVector([13, 14, 15])
  @test (a * b) * cpos ≈ a * (b * cpos)
end

@testset "Left division: a::Pose \\ bpos::StaticVector == inv(a) * bpos" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  bpos = @SVector([7, 8, 9])
  @test a \ bpos ≈ inv(a) * bpos
end

@testset "Generic vector (non-`StaticVector`) wrappers, sanity check" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  bpos = @SVector([7, 8, 9])
  bpos_ = [7, 8, 9]
  @test a * bpos_ ≈ a * bpos
  @test a \ bpos_ ≈ a \ bpos
end

@testset "Equivalence of `Pose * Matrix` with `Pose * StaticVector` on each column" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  bpoints = reshape(7:21, (3, 5))
  @test a * bpoints ≈ reduce(hcat,
      a * SVector{3}(bpoint) for bpoint in eachcol(bpoints))
end

@testset "Edge case: Pose * (3-by-0 Matrix)" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  bpoints = zeros(3, 0)
  @test a * bpoints == zeros(3, 0)
end

@testset "Left division: a::Pose \\ bpoints::Matrix == inv(a) * bpoints" begin
  a = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  bpoints = reshape(7:21, (3, 5))
  @test a \ bpoints ≈ inv(a) * bpoints
end

@testset "`interp` endpoint conditions" begin
  a = Pose([10, 11, 12], IDENTITY_ORN)
  b = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  @test interp(a, b, 0) ≈ a
  @test interp(a, b, 1) ≈ b
end

@testset "`t ↦ interp(b, t)` is a homomorphism from (ℝ, +) to (poses, ⊗)" begin
  b = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  @test interp(b, 0.5) ⊗ interp(b, 0.7) ≈ interp(b, 1.2)
end

@testset "`interp(b, t) == b ⊗ ... ⊗ b` (`t` times) when `t` is an integer" begin
  b = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  @test interp(b, -1) ≈ IDENTITY_POSE ⊘ b
  @test interp(b, 2) ≈ b ⊗ b
  @test interp(b, 0.5) ⊗ interp(b, 0.5) ≈ b
end

# Regression test for https://github.com/probcomp/PoseComposition.jl/issues/11
@testset "Can stringify, print and `show` with no crash" begin
  pose = Pose([1, 2, 3], RotZYX(0.4, 0.5, 0.6))
  string(pose)
  print(devnull, pose)
  show(devnull, "text/plain", pose)
end


text:
# ---
# title: 1110. Delete Nodes And Return Forest
# id: problem1110
# author: Indigo
# date: 2022-03-19
# difficulty: Medium
# categories: Tree, Depth-first Search
# link: <https://leetcode.com/problems/delete-nodes-and-return-forest/description/>
# hidden: true
# ---
#
# Given the `root` of a binary tree, each node in the tree has a distinct value.
#
# After deleting all nodes with a value in `to_delete`, we are left with a
# forest (a disjoint union of trees).
#
# Return the roots of the trees in the remaining forest.  You may return the
# result in any order.
#
#
#
# **Example 1:**
#
# **![](https://assets.leetcode.com/uploads/2019/07/01/screen-
# shot-2019-07-01-at-53836-pm.png)**
#
#
#
#     Input: root = [1,2,3,4,5,6,7], to_delete = [3,5]
#     Output: [[1,2,null,4],[6],[7]]
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in the given tree is at most `1000`.
#   * Each node has a distinct value between `1` and `1000`.
#   * `to_delete.length <= 1000`
#   * `to_delete` contains distinct values between `1` and `1000`.
#
#
## @lc code=start
using LeetCode

function del_nodes1110(node::TreeNode{Int}, to_delete::Vector{Int})
    to_deletes = Set(to_delete)
    res = TreeNode{Int}[]
    queue = [(node, true)]

    while !isempty(queue)
        root, flg = popfirst!(queue)
        if root.val ∈ to_deletes
            !isnothing(root.left) && push!(queue, (root.left))
            !isnothing(root.right) && push!(queue, (root.right))
            continue
        end
        for child in (:left, :right)
            isnothing(getproperty(root, child)) && continue
            if getproperty(root, child).val ∈ to_deletes
                !isnothing(getproperty(root, child).left) &&
                    push!(queue, (getproperty(root, child).left, true))
                !isnothing(getproperty(root, child).right) &&
                    push!(queue, (getproperty(root, child).right, true))
                setproperty!(root, child, nothing)
            else
                push!(queue, (getproperty(root, child), false))
            end
        end
        flg && push!(res, root)
    end
    return res
end
## @lc code=end