text:
@testset "Testing Placements" begin
    # Strategy: Use a bunch of different architectures and a linegraph from
    # Chessboard. Make sure that all get mapped to the minimum.

    # Try the 2D architectures.
    # Use only the "ChessboardColor()" coloring scheme to avoid mismatch between
    # white and black squares that may happen with HashColor
    architectures = (
        architecture(4, Rectangle2D(), ChessboardColor()),
        architecture(4, Hexagonal2D(), ChessboardColor()),
    )

    taskgraphs = (
        linegraph(16, AllGray()), linegraph(16, OddEven()), linegraph(16, Quarters())
    )

    move_generators = (SA.CachedMoveGenerator,)

    iterator = Iterators.product(architectures, taskgraphs, move_generators)

    for (A, T, movegen) in iterator
        m = Map(Chess(), A, T)

        sa = SA.SAStruct(m)
        move_generator = movegen(sa)
        @time SA.place!(sa; movegen = move_generator)
        # Number of links in the linegraph is 15.
        #
        # All maps should reach this.
        @test Mapper2.SA.map_cost(sa) == 15

        SA.record(m, sa)
        # Route and see that post-routing number of links is correct.
        route!(m)
        @test MapperCore.total_global_links(m) == 15
    end

    # Now test 3D placements
    architectures = (architecture(3, Rectangle3D(), ChessboardColor()),)

    taskgraphs = (
        linegraph(27, AllGray()), linegraph(27, OddEven()), linegraph(27, Quarters())
    )

    move_generators = (SA.CachedMoveGenerator,)

    iterator3d = Iterators.product(architectures, taskgraphs, move_generators)

    for (A, T, movegen) in iterator3d
        m = Map(Chess(), A, T)

        sa = SA.SAStruct(m)
        move_generator = movegen(sa)
        @time SA.place!(sa; movegen = move_generator)

        # When checking the number of links, give some leeway for the mapper
        # to not quite reach the minimum.
        #
        # Notable, the Quarters() taskgraph has trouble reaching a final
        # objective of 26. I think the minimum may in fact be 27.
        @test Mapper2.SA.map_cost(sa) <= 27

        SA.record(m, sa)
        # Route and see that post-routing number of links is correct.
        route!(m)
        @test MapperCore.total_global_links(m) <= 27
    end
end

@testset "Testing Fanout Placement" begin
    m = Example1.make_fanout()
    sa = Mapper2.SA.SAStruct(m)

    Mapper2.SA.place!(sa)
    # The example is small enough that it should achieve this result every time.
    @test Mapper2.SA.map_cost(sa) == 14
end

ext: jl


text:
function test_parameterized_hist(bin_type::VariableWidth, search_algorithm, parallelization)
    if run_1D_tests
        test_parameterized_hist_1D(bin_type, search_algorithm, parallelization)
    end

    if run_2D_tests
        test_parameterized_hist_2D(bin_type, search_algorithm, parallelization)
    end
end

function test_parameterized_hist_1D(bin_type::VariableWidth, search_algorithm, parallelization)
    @testset "regression 1 1x4 and 4x1" begin
        img1 = [0x00, 0x40, 0x80, 0xc0]

        h = create_fast_histogram(bin_type, search_algorithm, parallelization, [range(0x00; stop = 0xff, length = 5)])

        increment_bins!(h, img1)
        fh_counts = counts(h)
        sb_counts = StatsBase.fit(StatsBase.Histogram, vec(img1), 0x00:UInt8(64):0x10f).weights
        @test fh_counts == sb_counts

        # Also test the transpose of img1
        zero!(h)
        increment_bins!(h, transpose(img1))
        fh_counts = counts(h)
        sb_counts = StatsBase.fit(StatsBase.Histogram, vec(transpose(img1)), 0x00:UInt8(64):0x10f).weights
        @test fh_counts == sb_counts
    end

    @testset "random regressions (1D, Float32)" begin
        for i = 1:100
            img_size = rand(2:64, 2)
            img1 = rand(img_size...)

            h = create_fast_histogram(
                bin_type,
                search_algorithm,
                parallelization,
                [range(Float32(0.0); stop = Float32(1.0), length = 17)],
            )

            increment_bins!(h, img1)

            fh_counts = counts(h)
            sb_counts =
                StatsBase.fit(
                    StatsBase.Histogram,
                    vec(img1),
                    range(Float32(0.0); stop = Float32(1.0), length = 17),
                ).weights

            @test fh_counts == sb_counts

            if fh_counts != sb_counts
                break
            end
        end
    end

    @testset "random regressions (1D, Float32, variable width)" begin
        for i = 1:100
            img_size = rand(2:64, 2)
            img1 = rand(img_size...)

            edge = [0.0, 0.1, 0.5, 0.9, 1.0]
            h = create_fast_histogram(bin_type, search_algorithm, parallelization, [edge])

            increment_bins!(h, img1)

            fh_counts = counts(h)
            sb_counts = StatsBase.fit(StatsBase.Histogram, vec(img1), edge).weights

            @test fh_counts == sb_counts

            if fh_counts != sb_counts
                break
            end
        end
    end
end

function test_parameterized_hist_2D(bin_type::VariableWidth, search_algorithm, parallelization)
    @testset "regression 1 2x2" begin
        img1 = [
            0x00 0x40
            0x80 0xc0
        ]

        img2 = [
            0x00 0x00
            0x40 0x00
        ]

        edge = range(0x00; stop = 0xff, length = 5)
        h = create_fast_histogram(bin_type, search_algorithm, parallelization, [edge, edge])

        increment_bins!(h, img1, img2)

        fh_counts = counts(h)
        sb_counts =
            StatsBase.fit(
                StatsBase.Histogram,
                (vec(img1), vec(img2)),
                (0x00:UInt8(64):0x10f, 0x00:UInt8(64):0x10f),
            ).weights

        @test fh_counts == sb_counts
    end

    @testset "random regressions (2D, UInt8, square)" begin
        for i = 1:100
            img_size = rand(2:64, 2)
            img1 = rand(0x00:0xff, img_size...)
            img2 = rand(0x00:0xff, img_size...)

            edge = range(0x00; stop = 0xff, length = 17)
            h = create_fast_histogram(bin_type, search_algorithm, parallelization, [edge, edge])

            increment_bins!(h, img1, img2)

            fh_counts = counts(h)

            # StatsBase needs one more edge than the number of weights, but that doesn't fit evenly into
            # the 8-bit range [0x00, 0xff], so we need to extend the upper bound a bit to get it to do
            # the "logical equivalent" to what FastHistograms does.
            sb_counts =
                StatsBase.fit(
                    StatsBase.Histogram,
                    (vec(img1), vec(img2)),
                    (0x00:UInt8(16):0x10f, 0x00:UInt8(16):0x10f),
                ).weights

            @test fh_counts == sb_counts

            if fh_counts != sb_counts
                break
            end
        end
    end

    @testset "random regressions (2D, Float32, square)" begin
        for i = 1:100
            img_size = rand(2:64, 2)
            img1 = rand(img_size...)
            img2 = rand(img_size...)

            edge = range(Float32(0.0); stop = Float32(1.0), length = 17)
            h = create_fast_histogram(bin_type, search_algorithm, parallelization, [edge, edge])

            increment_bins!(h, img1, img2)

            fh_counts = counts(h)
            sb_counts =
                StatsBase.fit(
                    StatsBase.Histogram,
                    (vec(img1), vec(img2)),
                    (
                        range(Float32(0.0); stop = Float32(1.0), length = 17),
                        range(Float32(0.0); stop = Float32(1.0), length = 17),
                    ),
                ).weights

            @test fh_counts == sb_counts

            if fh_counts != sb_counts
                break
            end
        end
    end

    @testset "random regressions (2D, Float32, non-square)" begin
        for i = 1:100
            img_size = rand(2:64, 2)
            img1 = rand(img_size...)
            img2 = rand(img_size...)

            h = create_fast_histogram(
                bin_type,
                search_algorithm,
                parallelization,
                [
                    range(Float32(0.0); stop = Float32(1.0), length = 9),
                    range(Float32(0.0); stop = Float32(1.0), length = 17),
                ],
            )

            increment_bins!(h, img1, img2)

            fh_counts = counts(h)
            sb_counts =
                StatsBase.fit(
                    StatsBase.Histogram,
                    (vec(img1), vec(img2)),
                    (
                        range(Float32(0.0); stop = Float32(1.0), length = 9),
                        range(Float32(0.0); stop = Float32(1.0), length = 17),
                    ),
                ).weights

            @test fh_counts == sb_counts

            if fh_counts != sb_counts
                break
            end
        end
    end

    @testset "random regressions (2D, Float32, non-square, variable width)" begin
        for i = 1:100
            img_size = rand(2:64, 2)
            img1 = rand(img_size...)
            img2 = rand(img_size...)

            edge1 = [0.0, 0.1, 0.5, 0.9, 1.0]
            edge2 = [0.0, 0.5, 0.6, 1.0]
            h = create_fast_histogram(bin_type, search_algorithm, parallelization, [edge1, edge2])

            increment_bins!(h, img1, img2)

            fh_counts = counts(h)
            sb_counts = StatsBase.fit(StatsBase.Histogram, (vec(img1), vec(img2)), (edge1, edge2)).weights

            @test fh_counts == sb_counts

            if fh_counts != sb_counts
                break
            end
        end
    end
end


text:
# Unit tests for functions in src/algebra.jl file
using COSMO, Test, Random, LinearAlgebra, SparseArrays

rng = Random.MersenneTwister(9191)

@testset "Algebra" begin

  # clip()
  n = 50
  x = 10 * rand(rng, n) .- 5
  x[10] = -8
  x[20] = 8
  x_c = zeros(n)
  x_c2 = zeros(n)
  @. x_c = COSMO.clip(x, -2, 2)
  @test (maximum(x_c) <= 2. && minimum(x_c) >= -2.)
  @. x_c2 = COSMO.clip(x, -2, 2, -100, 100)
  @test (x_c2[10] == -100 && x_c2[20] == 100)

  # scaled_norm
  E = Diagonal([1, 2, 3.])
  v = [-1, 5, 4]
  @test COSMO.scaled_norm(E, v, 1) == 23
  @test COSMO.scaled_norm(E, v, 2) == norm([-1; 10; 12], 2)
  @test COSMO.scaled_norm(E, v, Inf) == 12
  @test_throws ArgumentError COSMO.scaled_norm(E, v, 3);

  # col_norms!()
  v = zeros(4)
  A = [1. 2 3 4; 2 -1 30 4.1]
  COSMO.col_norms!(v, A)
  @test v == [2; 2; 30; 4.1]

  A = sparse([1. 2 3 4; 2 -1 30 4.1])
  B = sparse([1 3 -2 4.2; -100 -1 -2 -100])
  COSMO.col_norms!(v, A, reset = true)
  COSMO.col_norms!(v, B, reset = false)
  @test v == [100; 3; 30; 100]

  # row_norms!()
  A = Matrix(A')
  B = Matrix(B')
  COSMO.row_norms!(v, A, reset = true)
  COSMO.row_norms!(v, B, reset = false)
  @test v == [100; 3; 30; 100]
  A = sparse(A)
  B = sparse(B)
  COSMO.row_norms!(v, A, reset = true)
  COSMO.row_norms!(v, B, reset = false)
  @test v == [100; 3; 30; 100]

  # lmul!()
  L = Diagonal(rand(rng, 5))
  x = rand(5)
  y = copy(x)
  COSMO.lmul!(L, x)
  @test norm(x - L * y, Inf ) <= 1e-10

  y = copy(x)
  COSMO.lmul!(I, x)
  @test y == x

  # lrmul!()
  L = Diagonal(rand(rng, 5))
  Mo = rand(rng, 5, 5)
  M = deepcopy(Mo)
  R = Diagonal(rand(rng, 5))
  COSMO.lrmul!(L, M, R)
  @test norm(M -  L * Mo * R, Inf) <= 1e-10

  M = deepcopy(Mo)
  COSMO.lrmul!(L, M, I)
  @test norm(M -  L * Mo, Inf) <= 1e-10

  M = deepcopy(Mo)
  COSMO.lrmul!(I, M, I)
  @test M == Mo

  # symmetrize!()
  A = rand(rng, 3, 3)
  B = deepcopy(A)
  COSMO.symmetrize_upper!(A)
  @test Symmetric(A) == (0.5 * (B + B'))
  @test_throws AssertionError COSMO.symmetrize_upper!(rand(2, 3));

  # is_pos_def()
  A = Symmetric(randn(rng, 4, 4))
  Apos = A + 4 * Matrix(1.0I, 4, 4)
  Aneg = A - 4 * Matrix(1.0I, 4, 4)
  @test COSMO.is_pos_def(Apos, 1e-6)
  @test !COSMO.is_pos_def(Aneg, 1e-6)
  @test COSMO.is_neg_def(Aneg, 1e-6)
  @test !COSMO.is_neg_def(Apos, 1e-6)

end
nothing


text:
using SciPy
using Test

@testset "optimize.linprog" begin
    A = [-3 1;
          1 2]
    b = [6, 4]
    c = [-1, 4]
    x0_bounds = (nothing, nothing)
    x1_bounds = (-3, nothing)
    res = optimize.linprog(c, A_ub=A, b_ub=b, bounds=[x0_bounds, x1_bounds])
    # sample
    # @test res == Dict{Any, Any}("fun" => -21.99999984082497,
    #                           "nit" => 6,
    #                           "status" => 0,
    #                           "message" => "Optimization terminated successfully.",
    #                           "success" => true,
    #                           "x" => [9.999999890483512, -2.9999999875853645],
    #                           "slack" => [38.9999996590359, 8.468721723886574e-8], "con" => Float64[])
    @test res["status"] == 0
end


text:
__precompile__()
"""
Specialized logging for testing and examples in WebSockets.

To avoid type piracy, defines _show for types where specialized show methods exist,
and falls back to 'show' where they don't.

When logging to a file, it's beneficial to keep a file open during the whole logging sesssion,
since opening and closing files is so slow it may affect the sequence of things.

This module dispatches on an ad-hoc type AbstractDevice, instead of putting
additional data on an IOContext method like 'show' does. When AbstracDevice points to
Base.NullDevice(), the input arguments are processed before the call is made, but
that is all.

In Julia 0.7, current logging functionality is replaced with macros, which
can be even faster. Macros can also retrieve the current function's name without using stacktraces.
With this module, each function defines id = "thisfunction".

Methods in this file have no external dependencies. Methods with dependencies are
loaded from separate files with @require. This adds to loading time.
"""
module logutils_ws
using Requires
import Base.text_colors
import Base.color_normal
import Base.text_colors
import Base.show
@require HttpServer include("log_httpserver.jl")
@require HTTP include("log_http.jl")
@require WebSockets include("log_ws.jl")
export clog
export clog_notime
export zlog
export zlog_notime
export logto
export loggingto
export directto_abstractdevice
export AbstractDevice
export DataDispatch
export zflush

const ColorDevices = Union{Base.TTY, IOBuffer, IOContext, Base.PipeEndpoint}
const BlackWDevices = Union{IOStream, IOBuffer} # and endpoint...
const LogDevices = Union{ColorDevices, BlackWDevices, Base.DevNullStream}
abstract type AbstractDevice{T} end
struct NullDevice <: AbstractDevice{NullDevice}
        s::Base.DevNullStream
    end
struct ColorDevice{S<:ColorDevices}<:AbstractDevice{ColorDevice}
        s::S
        end
struct BlackWDevice{S<:Union{IOStream, IOBuffer}}<:AbstractDevice{BlackWDevice}
        s::S
        end
mutable struct LogDevice
        s::AbstractDevice
        end
_devicecategory(::AbstractDevice{S}) where S = S
const CURDEVICE = LogDevice(NullDevice(Base.DevNullStream()))
struct DataDispatch
    data::Array{UInt8,1}
    contenttype::String
end
"""
Redirect coming zlog calls to a stream. Default is no logging.
clog calls will duplicate to STDOUT.
"""
logto(io::ColorDevices) =  CURDEVICE.s = ColorDevice(io)
logto(io::BlackWDevices) =  CURDEVICE.s = BlackWDevice(io)
logto(io::Base.DevNullStream) =  CURDEVICE.s = NullDevice(io)

"""
Returns the current logging stream
"""
loggingto() = CURDEVICE.s.s

"""
Log to (default) nothing, or streams as given in CURDEVICE.
First argument is expected to be function id.
"""
function zlog(vars::Vararg)
    _zlog(CURDEVICE.s, vars...)
    nothing
end
"""
Log to (default) NullDevice, or the device given in CURDEVICE.
Falls back to Base.showdefault when no special methods are defined.
"""
function zlog_notime(vars::Vararg)
    _zlog_notime(CURDEVICE.s, vars...)
    nothing
end
"""
Log to the given device, but also to STDOUT if that's not the given device.
"""
function clog(vars::Vararg)
    _zlog(CURDEVICE.s, vars...)
    _devicecategory(CURDEVICE.s) != ColorDevice && _zlog(ColorDevice(STDOUT), vars...)
    nothing
end
"""
Log to the given device, but also to STDOUT if that's not the given device.
"""
function clog_notime(vars::Vararg)
    _zlog_notime(CURDEVICE.s, vars...)
    _devicecategory(CURDEVICE.s) != ColorDevice  && _zlog_notime(ColorDevice(STDOUT), vars...)
    nothing
end
"""
Flushing file write buffers, only has an effect on streams.
We do not flush automatically after every log, because that
has a relatively dramatic effect on logging speeds.
"""
zflush() = isa(CURDEVICE.s.s, IOStream) && flush(CURDEVICE.s.s)

## Below are internal functions starting with _
_zlog(::NullDevice, ::Vararg) = nothing
_zlog_notime(::NullDevice, ::Vararg) = nothing
"First argument padded and emphasized.
End the original argument list with a :normal and a new line argument."
function _zlog(d::AbstractDevice, vars::Vararg{Any,N}) where N
    if N == 1
        _log(d, :normal, _tg(), "  ", :bold, :cyan, vars[1], :normal, "\n")
    else
        lid = 26
        pads = repeat(" ", max(0, lid - length(_string(vars[1])))) #rpad don't work with color codes
        _log(d, :normal, _tg(), "  ", :bold, :cyan, vars[1], :normal, pads, vars[2:end]..., :normal, "\n")
    end
    nothing
end
"zlog, but no time stamp and no different format first argument."
function _zlog_notime(d::AbstractDevice, vars::Vararg{Any,N}) where N
    # End the original argument list with a :normal and a new line argument.
    _log(d, :bold, :cyan, vars[1:end]..., :normal, "\n")
    nothing
end
"Write to blackwdevices, inapplicable formatting neglegted. No linefeed."
function _log(bwd::BlackWDevice, vars::Vararg)
    _show.(bwd, vars)
    nothing
end

"Write colordevices. Does not reset color codes after, no linefeed."
function _log(cd::ColorDevice, vars::Vararg)
    buf = ColorDevice(IOBuffer())
    _show.(buf, vars)
    write(cd.s, take!(buf.s))
    nothing
end
"Write directly to colordevice/ IOBuffers"
function _log(cd::ColorDevice{Base.AbstractIOBuffer{Array{UInt8,1}}}, vars::Vararg)
    _show.(cd, vars)
    nothing
end

"Return a string by emulating a bwdevice, where color codes are not output"
function _string(vars::Vararg)
    buf = BlackWDevice(IOBuffer())
    _show.(buf, vars)
    buf.s |> take! |> String
end

"_show takes just device and one other argument.
It falls back to normal show for nondefined _show methods."
_show(d::AbstractDevice, var) = show(d.s, var)
_show(d::AbstractDevice, s::AbstractString) = write(d.s, s);nothing

function _show(d::ColorDevice, sy::Symbol)
    co =  get(text_colors, sy, :NA)
    if co != :NA
        write(d.s, co)
    else
        write(d.s, sy)
    end
    nothing
end
function _show(d::ColorDevices, sy::Symbol)
    co =  get(text_colors, sy, :NA)
    if co != :NA
        write(d.s, co)
    else
        write(d.s, sy)
    end
    nothing
end

function _show(d::BlackWDevice, sy::Symbol)
    co =  get(text_colors, sy, :NA)
    if co == :NA
        write(d.s, ":", sy)
    end
    nothing
end

function _show(d::AbstractDevice, ex::Exception)
    _log(d, Base.warn_color(), typeof(ex), "\n")
    showerror(d.s, ex, [])
    _log(d, :normal, "\n")
end
function _show(d::AbstractDevice, err::ErrorException)
    _log(d, Base.error_color(), typeof(err), "\n")
    showerror(d.s, err, [])
    _log(d, :normal, "\n")
end
function _show(d::AbstractDevice, err::Base.UVError)
    _log(d, Base.error_color(), typeof(err), "\n")
    showerror(d.s, err, [])
    _log(d, :normal, "\n")
end

"Print dict, no heading, three pairs per line, truncate end to fit"
function _show(d::AbstractDevice, di::Dict)
    linelength = displaysize(STDOUT)[2]
    indent = 8
    npa = 3
    plen = div(linelength - indent, npa)
    pairs = collect(di)
    lpa = length(pairs)
    _log(d, :bold, :blue)
    for i = 1:npa:lpa
        write(d.s, " "^indent)
        write(d.s, _pairpad(pairs[i], plen))
        i+1 <= lpa   &&  write(d.s, _pairpad(pairs[i + 1], plen))
        i+2 <= lpa   &&  write(d.s, _pairpad(pairs[i + 2], plen))
        write(d.s, "\n")
    end
    _log(d, :normal)
    nothing
end
"Print array of pairs, no heading, three pairs per line, truncate end to fit"
function _show(d::AbstractDevice, pairs::Vector{Pair{SubString{String},SubString{String}}})
    linelength = 95
    indent = 8
    npa = 3
    plen = div(linelength - indent, npa)
    lpa = length(pairs)
    _log(d, :bold, :blue)
    for i = 1:npa:lpa
        write(d.s, " "^indent)
        write(d.s, _pairpad(pairs[i], plen))
        i+1 <= lpa   &&  write(d.s, _pairpad(pairs[i + 1], plen))
        i+2 <= lpa   &&  write(d.s, _pairpad(pairs[i + 2], plen))
        write(d.s, "\n")
    end
    _log(d, :normal)
    nothing
end
"Print dict, no heading, two pairs per line, truncate end to fit"
function _show(d::AbstractDevice, di::Dict{String, Function})
    linelength = 95
    indent = 8
    npa = 2
    plen = div(linelength - indent, npa)
    pairs = collect(di)
    lpa = length(pairs)
    _log(d, :bold, :blue)
    for i = 1:npa:lpa
        write(d.s, " "^indent)
        write(d.s, _pairpad(pairs[i], plen))
        i+1 <= lpa   &&  write(d.s, _pairpad(pairs[i + 1], plen))
        write(d.s, "\n")
    end
    _log(d, :normal)
    nothing
end

_pairpad(pa::Pair, plen::Int) = Base.cpad(_limlen(_string(pa), plen) , plen )
_string(pa::Pair) = _string(pa[1]) * " => " * _string(pa[2])
function _show(d::AbstractDevice, f::Function)
    mt = typeof(f).name.mt
    fnam = splitdir(string(mt.defs.func.file))[2]
    write(d.s, string(f) * " at " * fnam * ":"
        * string(mt.defs.func.line))
    nothing
end

"Type info not printed here as it is assumed the type is given by the context."
function _show(d::ColorDevice, stream::Base.LibuvStream)
    _log(d, "(", :bold, _uv_status(stream)..., :normal)
    nba = Base.nb_available(stream.buffer)
    nba > 0 && print(d.s, ", ", Base.nb_available(stream.buffer)," bytes waiting")
    print(d.s, ")")
    nothing
end
function _uv_status(x)
    s = x.status
    if x.handle == Base.C_NULL
        if s == Base.StatusClosed
            return :red, "✘" #"closed"
        elseif s == Base.StatusUninit
            return :red, "null"
        end
        return :red, "invalid status"
    elseif s == Base.StatusUninit
        return :yellow, "uninit"
    elseif s == Base.StatusInit
        return :yellow, "init"
    elseif s == Base.StatusConnecting
        return :yellow, "connecting"
    elseif s == Base.StatusOpen
        return :green, "✓"   # "open"
    elseif s == Base.StatusActive
        return :green, "active"
    elseif s == Base.StatusPaused
        return :red, "paused"
    elseif s == Base.StatusClosing
        return :red, "closing"
    elseif s == Base.StatusClosed
        return :red, "✘" #"closed"
    elseif s == Base.StatusEOF
        return :yellow, "eof"
    end
    return :red, "invalid status"
end

function _show(d::AbstractDevice, serv::Base.LibuvServer)
    _log(d, typeof(serv), "(", :bold, _uv_status(serv)..., :normal, ")")
    nothing
end

"Data as a truncated string"
function _show(d::AbstractDevice, datadispatch::DataDispatch)
    _showdata(d, datadispatch.data, datadispatch.contenttype)
    write(d.s, "\n")
    nothing
end

function _showdata(d::AbstractDevice, data::Array{UInt8,1}, contenttype::String)
    if ismatch(r"(text|script|html|xml|julia|java)", lowercase(contenttype))
        _log(d, :green, "\Data length: ", length(data), " ", :bold, :blue)
        s = data |> String |> _limlen
        write(d.s, replace(s, r"\s+", " "))
    else
        _log(d, :green, "\Data length: ", length(data), "  ", :blue)
        write(d.s,  data |> _limlen)
    end
    nothing
end

"Truncates for logging"
_limlen(data::AbstractString) = _limlen(data, 74)
function _limlen(data::AbstractString, linelength::Int)
    le = length(data)
   if le <  linelength
        return  normalize_string(string(data), stripcc = true)
    else
        adds = " … "
        addlen = length(adds)
        truncat = 2 * div(linelength, 3)
        tail = linelength - truncat - addlen - 1
        truncstring = String(data)[1:truncat] * adds * String(data)[end-tail:end]
        return normalize_string(truncstring, stripcc = true)
    end
end
function _limlen(data::Union{Vector{UInt8}, Vector{Float64}})
    le = length(data)
    maxlen = 12 # elements, not characters
	if le <  maxlen
		return  string(data)
    else
        adds =  " ..... "
        addlen = 2
        truncat = 2 * div(maxlen, 3)
        tail = maxlen - truncat - addlen - 1
        return string(data[1:truncat])[1:end-1] * adds * string(data[end-tail:end])[7:end]
	end
end

"Time group. show() converts to string only when necessary."
_tg() = Dates.Time(now())

"For use in show(io::IO, obj) methods. Hook into this logger's dispatch mechanism."
function directto_abstractdevice(io::IO, obj)
	if isa(io, ColorDevices)
		buf = ColorDevice(IOBuffer())
	else
		buf = BlackWDevice(IOBuffer())
	end
	_show(buf, obj)
    write(io, take!(buf.s))
    nothing
end

nothing
end # module


text:
using Poly
using StaticArrays
using LinearAlgebra
using BenchmarkTools
using SIMD
using Base.Threads

num_threads = Base.Threads.nthreads()
LinearAlgebra.BLAS.set_num_threads(num_threads)
thread = false
if num_threads != 1
    thread = true
end

function jacobi_naive(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    for t = 1:timesteps-1
        for i = 3:n-2
            b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        for ii = 3:n-2
            a[ii] = b[ii]
        end
    end
end

function jacobi_simple(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    @inbounds for t = 1:timesteps-1
        for i = 3:n-2
            b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        for ii = 3:n-2
            a[ii] = b[ii]
        end
    end
end

function jacobi_optimized(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    @simd for t = 1:timesteps-1
        @simd for i = 3:n-2
            @inbounds b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        @simd for ii = 3:n-2
            @inbounds a[ii] = b[ii]
        end
    end
end

function jacobi_optimized_thread(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    @threads for t = 1:timesteps-1
        @simd for i = 3:n-2
            @inbounds b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        @simd for ii = 3:n-2
            @inbounds a[ii] = b[ii]
        end
    end
end

function jacobi_poly(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    @poly_loop tile=0 for t = 1:timesteps-1
        for i = 3:n-2
            b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        for ii = 3:n-2
            a[ii] = b[ii]
        end
    end
end

function jacobi_poly_rt(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    @poly_loop tile=0 for t = 1:$timesteps-1
        for i = 3:$n-2
            b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        for ii = 3:$n-2
            a[ii] = b[ii]
        end
    end
end

function jacobi_poly_thread(a::Array{T,1}, b::Array{T,1}, timesteps, n) where {T}
    @poly_loop thread=true tile=0 for t = 1:timesteps-1
        for i = 3:n-2
            b[i] = (a[i - 1] + a[i] + a[i + 1])/3
        end
        for ii = 3:n-2
            a[ii] = b[ii]
        end
    end
end

n = 1024
A = rand(n)
B = zeros(n)
t = 512

basic = @benchmark jacobi_naive(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
simple = @benchmark jacobi_simple(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
if thread
    optimized = @benchmark jacobi_optimized_thread(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
    poly = @benchmark jacobi_poly_thread(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
    polyrt = nothing
else
    optimized = @benchmark jacobi_optimized(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
    poly = @benchmark jacobi_poly(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
    polyrt = @benchmark jacobi_poly_rt(in, B, $t, $n) setup=(in=copy(A); B=zeros(n))
end

@show allocs(basic)
@show allocs(simple)
@show allocs(optimized)
@show allocs(poly)
if !thread
    @show allocs(polyrt)
end
@show minimum(basic)
@show minimum(simple)
@show minimum(optimized)
@show minimum(poly)
if !thread
    @show minimum(polyrt)
end

rbasic = ratio(minimum(basic), minimum(poly))
rsimple = ratio(minimum(simple), minimum(poly))
roptimized = ratio(minimum(optimized), minimum(poly))

if !thread
    rbasicrt = ratio(minimum(basic), minimum(polyrt))
    rsimplert = ratio(minimum(simple), minimum(polyrt))
    roptimizedrt = ratio(minimum(optimized), minimum(polyrt))
end

@show rbasic
@show rsimple
@show roptimized

if !thread
    @show rbasicrt
    @show rsimplert
    @show roptimizedrt
end

jbasic = judge(minimum(poly), minimum(basic))
jsimple = judge(minimum(poly), minimum(simple))
joptimized = judge(minimum(poly), minimum(optimized))

if !thread
    jbasicrt = judge(minimum(polyrt), minimum(basic))
    jsimplert = judge(minimum(polyrt), minimum(simple))
    joptimizedrt = judge(minimum(polyrt), minimum(optimized))
end

@show jbasic
@show jsimple
@show joptimized

if !thread
    @show jbasicrt
    @show jsimplert
    @show joptimizedrt
end


text:
@doc raw"""
Tableau for variational Lobatto IIIA method with s stages

Uses Lobatto IIIA for the coefficients $a_{ij}$ and its symplectic conjugate, Lobatto IIIB, for the coefficients $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIA(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIA$(s)"), 2s-2, TableauLobattoIIIA(T,s), (-1)^(s+1), get_lobatto_nullvector(s))
end

TableauVPLobattoIIIA(s) = TableauVPLobattoIIIA(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIIB method with s stages

Uses Lobatto IIIB for the coefficients $a_{ij}$ and its symplectic conjugate, Lobatto IIIA, for the coefficients $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIB(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIB$(s)"), 2s-2, TableauLobattoIIIB(T,s), (-1)^(s+1), get_lobatto_nullvector(s))
end

TableauVPLobattoIIIB(s) = TableauVPLobattoIIIB(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIIC method with s stages

Uses Lobatto IIIC for the coefficients $a_{ij}$ and its symplectic conjugate, Lobatto IIIC̄, for the coefficients $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIC(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIC$(s)"), 2s-2, TableauLobattoIIIC(T,s), (-1)^(s+1))
end

TableauVPLobattoIIIC(s) = TableauVPLobattoIIIC(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIIC̄ method with s stages

Uses Lobatto IIIC̄ for the coefficients $a_{ij}$ and its symplectic conjugate, Lobatto IIIC, for the coefficients $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIC̄(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIC̄$(s)"), 2s-2, TableauLobattoIIIC̄(T,s), (-1)^(s+1))
end

TableauVPLobattoIIIC̄(s) = TableauVPLobattoIIIC̄(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIID method with s stages

Uses Lobatto IIID for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIID(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIID$(s)"), 2s-s, TableauLobattoIIID(T,s), (-1)^s)
end

TableauVPLobattoIIID(s) = TableauVPLobattoIIID(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIIE method with s stages

Uses Lobatto IIIE for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIE(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIE$(s)"), 2, TableauLobattoIIIE(T,s), (-1)^s)
end

TableauVPLobattoIIIE(s) = TableauVPLobattoIIIE(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIIF method with s stages

Uses Lobatto IIIF for the coefficients $a_{ij}$ and its symplectic conjugate for the coefficients $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIF(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIF$(s)"), 2, TableauLobattoIIIF(T,s), (-1)^s)
end

TableauVPLobattoIIIF(s) = TableauVPLobattoIIIF(Float64,s)

@doc raw"""
Tableau for variational Lobatto IIIG method with s stages

Uses Lobatto IIIG for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIG(::Type{T}, s) where {T}
    TableauVPRK(Symbol("LobattoIIIG$(s)"), 2, TableauLobattoIIIG(T,s), (-1)^s)
end

TableauVPLobattoIIIG(s) = TableauVPLobattoIIIG(Float64,s)

@doc raw"""
Tableau for variational Gauss-Legendre method with s stages

Uses Gauss coefficients for both $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPGLRK(::Type{T}, s) where {T}
    TableauVPRK(Symbol("GLRK", s), 2s, TableauGauss(T, s), (-1)^s)
end

TableauVPGLRK(s) = TableauVPGLRK(Float64,s)

@doc raw"""
Tableau for variational symmetric Runge-Kutta method with 3 stages

Uses SRK3 coefficients for both $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPSRK3(T=Float64)
    srk = TableauSRK3(T)
    TableauVPRK(:SRK3, srk.o, srk, (-1)^3)
end

@doc raw"""
Tableau for an implicit partitioned Lobatto IIIA Runge-Kutta method with s stages

Uses Lobatto IIIA for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIAIIIA(::Type{T}, s) where {T}
    lob = TableauLobattoIIIA(T,s)
    TableauVPRK(Symbol("LobattoIIIAIIIA$s"), 2s-2, lob, lob, (-1)^(s+1))
end

TableauVPLobattoIIIAIIIA(s) = TableauVPLobattoIIIAIIIA(Float64,s)

@doc raw"""
Tableau for an implicit partitioned Lobatto IIIB Runge-Kutta method with s stages

Uses Lobatto IIIB for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPLobattoIIIBIIIB(::Type{T}, s) where {T}
    lob = TableauLobattoIIIB(T,s)
    TableauVPRK(Symbol("LobattoIIIBIIIB$s"), 2s-2, lob, lob, (-1)^(s+1))
end

TableauVPLobattoIIIBIIIB(s) = TableauVPLobattoIIIBIIIB(Float64,s)

@doc raw"""
Tableau for an implicit partitioned Radau IIA Runge-Kutta method with s stages

Uses Radau IIA for both coefficients $a_{ij}$ and $\bar{a}_{ij}$.
"""
function TableauVPRadauIIAIIA(::Type{T}, s) where {T}
    rad = TableauRadauIIA(T,s)
    TableauVPRK(Symbol("RadauIIAIIA$s"), 2s-1, rad, rad, (-1)^(s+1))
end

TableauVPRadauIIAIIA(s) = TableauVPRadauIIAIIA(Float64,s)


text:
module HuggingFace

using ..Transformers

export get_or_download_hgf_file,
  get_or_download_hgf_config,
  get_or_download_hgf_weight,
  @hgf_str,
  load_config,
  load_model,
  load_model!,
  load_state_dict,
  load_state,
  save_config,
  save_model

include("./download.jl")
include("./weight.jl")
include("./configs/config.jl")
include("./models/models.jl")

"""
  `hgf"<model-name>:<item>"`

Get `item` from `model-name`. This will ensure the required
data are downloaded and registered. `item` can be "config",
"tokenizer", and model related like "model", or "formaskedlm", etc. Use [`get_model_type`](@ref) to see what
model/task are supported.
"""
macro hgf_str(name)
  :(load_hgf_pretrained($(esc(name))))
end

"""
  `get_model_type(::Val{model})`

See the list of supported model type of given model.
For example, use `get_mdoel_type(Val(:bert))` to
see all model/task that `bert` support.
"""
get_model_type

get_model_type(model, task) = error("Model $model doesn't support this kind of task: $task")

"""
  `load_hgf_pretrained(name)`

The underlying function of [`@hgf_str`](@ref).
"""
function load_hgf_pretrained(name)
  model_name, item = rsplit(name, ':'; limit=2)
  get_or_download_hgf_config(model_name)
  cfg = load_config(model_name)

  item == "config" &&
    return cfg

  item == "tokenizer" &&
    error("tokenizer not support yet.")

  model_type = cfg.model_type
  model_cons = get_model_type((Val ∘ Symbol ∘ lowercase)(model_type), (Val ∘ Symbol ∘ lowercase)(item))

  get_or_download_hgf_weight(model_name)
  model = load_model(model_cons, model_name; config=cfg)

  return model
end

end


text:
function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray::Nothing = nothing) where PT
    hull = PT[]
    if length(ps) == 0
        return hull
    end
    prev = sign_sense == 1 ? first(ps) : last(ps)
    cur = prev
    for j in (sign_sense == 1 ? (2:length(ps)) : ((length(ps)-1):-1:1))
        while prev != cur && counterclockwise(cur - prev, ps[j] - prev) >= 0
            cur = prev
            pop!(hull)
            if !isempty(hull)
                prev = last(hull)
            end
        end
        push!(hull, cur)
        prev = cur
        cur = ps[j]
    end
    push!(hull, cur)
    return hull
end

function getsemihull(ps::Vector{PT}, sign_sense, counterclockwise, yray) where PT
    hull = getsemihull(ps, sign_sense, counterclockwise)
    while length(hull) >= 2 && counterclockwise(hull[end] - hull[end - 1], yray) >= 0
        pop!(hull)
    end
    return hull
end

function _planar_hull(d::FullDim, points, lines, rays, counterclockwise, rotate)
    line = nothing
    lineleft = false
    lineright = false
    function checkleftright(r::Union{Ray, Line})
        cc = counterclockwise(r, line)
        if !isapproxzero(cc)
            if cc < 0 || islin(r)
                lineleft = true
            end
            if cc > 0 || islin(r)
                lineright = true
            end
        end
    end
    for l in lines
        if !isapproxzero(l)
            if line === nothing
                line = l
            else
                checkleftright(l)
            end
        end
    end
    xray = yray = nothing
    for r in rays
        isapproxzero(r) && continue
        if line === nothing
            if xray === nothing
                xray = yray = coord(r)
            else
                if Line(xray) ≈ linearize(r) && !(Ray(xray) ≈ r)
                    line = Line(xray)
                    checkleftright(Ray(yray))
                elseif Line(yray) ≈ linearize(r) && !(Ray(yray) ≈ r)
                    line = Line(yray)
                    checkleftright(Ray(xray))
                else
                    x_right = counterclockwise(r, xray) > 0
                    y_left = counterclockwise(r, yray) < 0
                    if x_right
                        if y_left
                            line = Line(xray)
                            lineleft = lineright = true
                        else
                            xray = coord(r)
                        end
                    else
                        if y_left
                            yray = coord(r)
                        end
                    end
                end
            end
        else
            checkleftright(r)
        end
    end
    if line === nothing
        if xray === nothing
            sweep_norm = rotate(basis(eltype(points), d, 1))
            if iszero(sum(abs, sweep_norm))
                sweep_norm = rotate(basis(eltype(points), d, 2))
            end
        else
            sweep_norm = rotate(xray)
        end
    else
        sweep_norm = rotate(coord(line))
    end
    sort!(points, by = x -> dot(x, sweep_norm))
    _points = eltype(points)[]
    _lines = eltype(lines)[]
    _rays = eltype(rays)[]
    if line === nothing
        append!(_points, getsemihull(points, 1, counterclockwise, yray))
        if yray === nothing
            append!(_points, getsemihull(points, -1, counterclockwise, yray)[2:end-1])
        else
            push!(_rays, Ray(xray))
            if !(Ray(xray) ≈ Ray(yray))
                push!(_rays, Ray(yray))
            end
        end
    else
        push!(_lines, line)
        push!(_points, first(points))
        if lineleft
            if lineright
                push!(_lines, Line(sweep_norm))
            else
                push!(_rays, Ray(sweep_norm))
            end
        elseif lineright
            push!(_rays, Ray(-sweep_norm))
        else
            if !(dot(first(points), sweep_norm) ≈ dot(last(points), sweep_norm))
                push!(_points, last(points))
            end
        end
    end
    return _points, _lines, _rays
end

counterclockwise(x, y) = x[1] * y[2] - x[2] * y[1]
rotate(x) = convert(typeof(x), StaticArrays.SVector(-x[2], x[1]))

function planar_hull(vr::VRepresentation)
    d = FullDim(vr)
    vrep(_planar_hull(FullDim(vr), collect(points(vr)), lines(vr), rays(vr), counterclockwise, rotate)...; d = d)
end


text:
using Blocks
using Compat

const datafile = joinpath(dirname(@__FILE__), "test.csv")
const nloops = 10

function testfn(f::Function, s::AbstractString, exp_res)
    println("\t$(s)...")
    ret = f()
    println("\t\tresult: $(ret)")
    @assert (ret == exp_res)
    println("\t\tresult: $(ret)")

    t = @elapsed for i in 1:nloops f(); end
    println("\t\ttime: $(t/nloops)")
end

#f_df_pmap() = pmapreduce(x->nrow(x), +, Block(File(datafile)) |> as_io |> as_recordio |> (x)->as_dataframe(x; header=false))
#f_df_map() = sum(map(x->nrow(x), Block(File(datafile)) |> as_io |> as_recordio |> (x)->as_dataframe(x; header=false)))
a_a1_pmap() = pmapreduce(x->sum(2*x), +, Block([1:1000;], 1, nworkers()))
a_a1_map() = sum(map(x->sum(2*x), Block([1:1000;], 1, nworkers())))
a_a2_pmap() = pmapreduce(x->sum(2*x), +, Block(reshape([1:1000;],10,10,10), [1,2]))
a_a2_map = ()->sum(map(x->sum(2*x), Block(reshape([1:1000;],10,10,10), [1,2])))
f_ios_pmap() = pmapreduce(x->length(x), +, Block(File(datafile)) |> as_io |> as_recordio |> as_lines)
f_ios_map() = sum(map(x->length(x), Block(File(datafile)) |> as_io |> as_recordio |> as_lines))
f_stream_pmap() = pmapreduce(x->length(x), +, @prepare Block(open(datafile), 1000) |> as_recordio |> as_lines)
f_stream_map() = sum(map(x->length(x), @prepare Block(open(datafile), 1000) |> as_recordio |> as_lines))

function procaff()
    b = Block(File(datafile)) |> as_io |> as_recordio |> as_lines
    # create random affinities
    wrkrids = workers()
    #nw = length(wrkrids)
    b.affinity = shuffle!(wrkrids) #map(x->randsample(wrkrids, min(3,nw)), 1:nw)
    pmapreduce(x->length(x), +, b)
end

function load_pkgs()
    println("loading packages...")
    @everywhere using Blocks
    @everywhere using Base.FS
end

function do_all_tests()
    println("running tests...")
    #testfn(f_df_pmap, "pmap file->dataframe", 100)
    #testfn(f_df_map, "map file->dataframe", 100)
    testfn(a_a1_pmap, "pmap array->array", 1001000)
    testfn(a_a1_map, "map array->array", 1001000)
    testfn(a_a2_pmap, "pmap array->array", 1001000)
    testfn(a_a2_map, "map array->array", 1001000)
    testfn(f_ios_pmap, "pmap file->iostream", 100)
    testfn(f_ios_map, "map file->iostream", 100)
    testfn(f_stream_pmap, "pmap stream->lines", 100)
    testfn(f_stream_map, "map stream->lines", 100)
    testfn(procaff, "pmap processor affinity", 100)
end

load_pkgs()
do_all_tests()

println("adding 3 more processors...")
addprocs(3)
println("\tnworkers: $(nworkers())")
load_pkgs()
do_all_tests()


text:
using PolyChaos, Test

deg = 5
op = genLaguerreOrthoPoly(deg,1.23)

n = 10
X = [ exp(3 + rand()), exp.(3 .+ rand(n)) ]

@testset "evaluate univariate basis at point(s)" begin
    for x in X
        @test evaluate(collect(0:deg),x,op.α,op.β) == evaluate(collect(0:deg),x,op)
        for d in 0:deg
            if typeof(x) <: Real
                @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op) == evaluate(d,[x],op.α,op.β)
            else
                @test evaluate(d,x,op.α,op.β) == evaluate(d,x,op)
            end
        end

        res = hcat(map(i->evaluate(i,x,op.α,op.β),collect(0:op.deg))...)
        if typeof(x) <: Real
            @test evaluate(x,op) == evaluate([x], op) == evaluate(collect(0:op.deg),[x],op) == evaluate(collect(0:op.deg),[x],op.α,op.β) == res
        else
            @test evaluate(x, op) == evaluate(collect(0:op.deg),x,op) == evaluate(collect(0:op.deg),x,op.α,op.β) == res
        end
    end
end

ops = [  genLaguerreOrthoPoly(deg,1.23),
        genHermiteOrthoPoly(deg,3.),
        LegendreOrthoPoly(deg),
        JacobiOrthoPoly(deg,4.3,10.)
     ]

mop = MultiOrthoPoly(ops,deg)

nop = length(ops)
X = [ exp.(3 .+ rand(nop)) , exp.(3 .+ rand(n,nop)) ]

@testset "evaluate multivariate basis at point(s)" begin
    for x in X
        @test evaluate(mop.ind,x,coeffs(mop)...) == evaluate(mop.ind,x,mop) == evaluate(x,mop)
        if typeof(x) <: Vector{<:Real}
            @test evaluate(mop.ind,x,coeffs(mop)...) == evaluate(mop.ind,x,mop) == evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...)
            @test evaluate(x,mop) == evaluate(mop.ind,reshape(x,1,length(x)),mop) == evaluate(mop.ind,reshape(x,1,length(x)),coeffs(mop)...)
        else
            @test evaluate(x,mop) == evaluate(mop.ind,x,mop) == evaluate(mop.ind,x,coeffs(mop)...)
        end
        for d in eachrow(mop.ind[1:3,:])
            @test evaluate(collect(d),x,coeffs(mop)...) == evaluate(collect(d),x,mop)
        end
    end
end


text:
"""
    Diagnostics

Accumulate mean fields and covariance statistics on the computational grid.
"""
module Diagnostics

export DiagnosticsGroup,
    setup_atmos_default_diagnostics,
    setup_atmos_core_diagnostics,
    setup_dump_state_and_aux_diagnostics

using CuArrays
using Dates
using FileIO
using JLD2
using KernelAbstractions
using MPI
using OrderedCollections
using Printf
using StaticArrays
import KernelAbstractions: CPU

using ..ConfigTypes
using ..DGMethods
using ..DGMethods:
    number_state_conservative,
    vars_state_conservative,
    number_state_auxiliary,
    vars_state_auxiliary,
    vars_state_gradient_flux,
    number_state_gradient_flux
using ..Mesh.Interpolation
using ..MPIStateArrays
using ..VariableTemplates
using ..Writers

using CLIMAParameters
using CLIMAParameters.Planet: planet_radius
struct EarthParameterSet <: AbstractEarthParameterSet end
const param_set = EarthParameterSet()

Base.@kwdef mutable struct Diagnostic_Settings
    mpicomm::MPI.Comm = MPI.COMM_WORLD
    dg::Union{Nothing, DGModel} = nothing
    Q::Union{Nothing, MPIStateArray} = nothing
    starttime::Union{Nothing, String} = nothing
    output_dir::Union{Nothing, String} = nothing
end
const Settings = Diagnostic_Settings()

"""
    init(mpicomm, dg, Q, starttime, output_dir)

Initialize the diagnostics collection module -- save the parameters into
`Settings`.
"""
function init(
    mpicomm::MPI.Comm,
    dg::DGModel,
    Q::MPIStateArray,
    starttime::String,
    output_dir::String,
)
    Settings.mpicomm = mpicomm
    Settings.dg = dg
    Settings.Q = Q
    Settings.starttime = starttime
    Settings.output_dir = output_dir
end

include("variables.jl")
include("helpers.jl")
include("atmos_common.jl")
include("thermo.jl")
include("groups.jl")

"""
    __init()__

Module initialization function. Currently, only fills in all currently
defined diagnostic variables.
"""
function __init__()
    setup_variables()
end

end # module Diagnostics

max_stars_repo_name: leios/CLIMA

max_issues_repo_name: leios/CLIMA

max_forks_repo_name: leios/CLIMA


text:
using SafeTestsets

@safetestset "read tests" begin include("read_tests.jl") end
@safetestset "write tests" begin include("write_tests.jl") end


text:
###############################################################################
#
# Initialize CanopyOpticals
#
###############################################################################
"""
    create_canopy_rads(FT, rt_dim::RTDimensions)

Create a [`CanopyRads`](@ref) struct, given
- `FT` Floating number type
- `rt_dim` [`RTDimensions`](@ref) type struct
"""
function create_canopy_rads(FT, rt_dim::RTDimensions)
    @unpack nAzi, nIncl, nLayer, nLevel, nWL, nWLF = rt_dim;

    return CanopyRads{FT}(nAzi   = nAzi  ,
                          nIncl  = nIncl ,
                          nLayer = nLayer,
                          nLevel = nLevel,
                          nWL    = nWL   ,
                          nWLF   = nWLF  )
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
using Distributed
@everywhere using Pkg
@everywhere Pkg.activate(@__DIR__)

@everywhere begin
    using RiceBPH
    using RiceBPH.PostProcess
    using ProgressMeter
    using CSV, DataFrames
    const testdir = "kiem-dinh"
    const resultdir = "results"
end

mkpath(testdir)

let cases = collect(Iterators.product([0.01, 0.05], [0.1, 0.05, 0.01, 0.005]))
    @showprogress pmap(cases) do (p0, a)
        df = PostProcess.batch_test_rice(resultdir, p0; alpha=a, clean=true)
        CSV.write(joinpath(testdir, "test-p0+$p0-alpha+$a.csv"), df)
    end
end

max_stars_repo_name: ndgnuh/RiceBPH

max_issues_repo_name: ndgnuh/RiceBPH

max_forks_repo_name: ndgnuh/RiceBPH


text:
using CryptoMiniSat
using Test

#using Aqua
#Aqua.test_all(CryptoMiniSat, deps_compat = true)

@testset "Low-level CryptoMiniSat" begin
    p = CryptoMiniSat.cms_new()
    CryptoMiniSat.cms_new_vars(p, 3)
    @test CryptoMiniSat.cms_nvars(p) == 3

    CryptoMiniSat.cms_add_clause(p, [Lit(0,false)])
    CryptoMiniSat.cms_add_clause(p, [Lit(1,true)])
    CryptoMiniSat.cms_add_clause(p, [Lit(0,true),Lit(1,false),Lit(2,false)])

    @test CryptoMiniSat.cms_solve(p) == CryptoMiniSat.L_TRUE

    @test CryptoMiniSat.cms_get_model(p) == [0x00,0x01,0x00]

    @test CryptoMiniSat.cms_solve_with_assumptions(p, [Lit(2,true)]) == CryptoMiniSat.L_FALSE
    @test CryptoMiniSat.cms_solve(p) == CryptoMiniSat.L_TRUE

    CryptoMiniSat.cms_add_clause(p, [Lit(2,true)])
    @test CryptoMiniSat.cms_solve(p) == CryptoMiniSat.L_FALSE

    CryptoMiniSat.cms_free(p)
end

@testset "Top-level PicoSAT-like interface" begin
    cnf = Any[[1, -5, 4], [-1, 5, 3, 4], [-3, -4]]
    @test CryptoMiniSat.solve(cnf) == [-1, -2, -3, -4, -5]
    @test sort(collect(CryptoMiniSat.itersolve(cnf))) == Any[[-1, -2, -3, -4, -5], [-1, -2, -3, 4, -5], [-1, -2, -3, 4, 5], [-1, -2, 3, -4, -5], [-1, 2, -3, -4, -5], [-1, 2, -3, 4, -5], [-1, 2, -3, 4, 5], [-1, 2, 3, -4, -5], [1, -2, -3, -4, 5], [1, -2, -3, 4, -5], [1, -2, -3, 4, 5], [1, -2, 3, -4, -5], [1, -2, 3, -4, 5], [1, 2, -3, -4, 5], [1, 2, -3, 4, -5], [1, 2, -3, 4, 5], [1, 2, 3, -4, -5], [1, 2, 3, -4, 5]]
end


text:
# io/z-test.jl

using Test: @testset

list = [
"caller_name.jl"
#"fld-read-deprecated.jl"
"fld-write-deprecated.jl"
"ir_dump.jl"
"prompt.jl"
"shows.jl"
]

for file in list
	@testset "$file" begin
		include(file)
	end
end


text:
export check_ctranspose, check_hermitian, check_positive_definite

"""
    check_ctranspose(op)
Cheap check that the operator and its conjugate transposed are related.
"""
function check_ctranspose(op::AbstractLinearOperator{T}) where {T <: Union{AbstractFloat, Complex}}
  (m, n) = size(op)
  x = rand(n)
  y = rand(m)
  yAx = dot(y, op * x)
  xAty = dot(x, op' * y)
  ε = eps(real(eltype(op)))
  return abs(yAx - conj(xAty)) < (abs(yAx) + ε) * ε^(1 / 3)
end

function check_ctranspose(op::AbstractLinearOperator{T}) where {T <: Integer}
  (m, n) = size(op)
  x = convert(Vector{T}, (floor.(10 * rand(n)))) .- 5
  y = convert(Vector{T}, (floor.(10 * rand(m)))) .- 5
  yAx = dot(y, op * x)
  xAty = dot(x, op' * y)
  return yAx == xAty
end

check_ctranspose(M::AbstractMatrix) = check_ctranspose(LinearOperator(M))

"""
    check_hermitian(op)
Cheap check that the operator is Hermitian.
"""
function check_hermitian(op::AbstractLinearOperator{T}) where {T <: Union{AbstractFloat, Complex}}
  m, n = size(op)
  m == n || throw(LinearOperatorException("shape mismatch"))
  v = rand(n)
  w = copy(op * v)  # copy necessary to guard against in-place operators
  s = dot(w, w)  # = (Av)'(Av) = v' A' A v.
  y = op * w
  t = dot(v, y)  # = v' A A v.
  ε = eps(real(eltype(op)))
  return abs(s - t) < (abs(s) + ε) * ε^(1 / 3)
end

function check_hermitian(op::AbstractLinearOperator{T}) where {T <: Integer}
  m, n = size(op)
  m == n || throw(LinearOperatorException("shape mismatch"))
  v = convert(Vector{T}, (floor.(10 * rand(n)))) .- 5
  w = copy(op * v)
  s = dot(w, w)  # = (Av)'(Av) = v' A' A v.
  y = op * w
  t = dot(v, y)  # = v' A A v.
  return s == t
end

check_hermitian(M::AbstractMatrix) = check_hermitian(LinearOperator(M))

"""
    check_positive_definite(op; semi=false)
Cheap check that the operator is positive (semi-)definite.
"""
function check_positive_definite(
  op::AbstractLinearOperator{T};
  semi = false,
) where {T <: Union{AbstractFloat, Complex}}
  m, n = size(op)
  m == n || throw(LinearOperatorException("shape mismatch"))
  v = rand(n)
  w = op * v
  vw = dot(v, w)
  ε = eps(real(eltype(op)))
  if imag(vw) > sqrt(ε) * abs(vw)
    return false
  end
  vw = real(vw)
  return semi ? (vw ≥ 0) : (vw > 0)
end

function check_positive_definite(op::AbstractLinearOperator{T}; semi = false) where {T <: Integer}
  m, n = size(op)
  m == n || throw(LinearOperatorException("shape mismatch"))
  v = convert(Vector{T}, (floor.(10 * rand(n)))) .- 5
  w = op * v
  vw = dot(v, w)
  return semi ? (vw ≥ 0) : (vw > 0)
end

check_positive_definite(M::AbstractMatrix; kwargs...) =
  check_positive_definite(LinearOperator(M); kwargs...)


text:
# Colormap parameters

# To avoid over-specialization, we use `Any` as the value type,
# but the actual type should be exactly as used here.
const colormaps_sequential = Dict{String, Any}(
    #single hue
    #name         hue    w     d     c     s     b     wcolor              dcolor
    "blues"   => (255.0, 0.3,  0.25, 0.88, 0.6,  0.75, RGB(1.0, 1.0, 0.0), RGB(0.0, 0.0, 1.0)),
    "greens"  => (120.0, 0.15, 0.18, 0.88, 0.55, 0.9,  RGB(1.0, 1.0, 0.0), RGB(0.0, 0.0, 1.0)),
    "grays"   => (  0.0, 0.0,  0.0,  1.0,  0.0,  0.75, RGB(1.0, 1.0, 0.0), RGB(0.0, 0.0, 1.0)),
    "oranges" => ( 20.0, 0.5,  0.4,  0.83, 0.95, 0.85, RGB(1.0, 1.0, 0.0), RGB(1.0, 0.0, 0.0)),
    "purples" => (265.0, 0.15, 0.2,  0.88, 0.5,  0.7,  RGB(1.0, 0.0, 1.0), RGB(1.0, 0.0, 0.0)),
    "reds"    => ( 12.0, 0.15, 0.25, 0.8,  0.85, 0.6,  RGB(1.0, 1.0, 0.0), RGB(0.3, 0.1, 0.1)),
)

const colormaps_diverging = Dict{String, Any}(
    #name          h1     h2     w     d1    d2    c     s     b     wcolor              dcolor1             dcolor2
    "rdbu"     => ( 12.0, 255.0, 0.2,  0.6,  0.0,  0.85, 0.6,  0.65, RGB(1.0, 1.0, 0.0), RGB(1.0, 0.0, 0.0), RGB(0.0, 0.0, 1.0)),
)


text:
using BaseBenchmarks
using BenchmarkTools
using Distributed

addprocs(1)

# re-tune the entire suite

BaseBenchmarks.loadall!(tune = false)
warmup(BaseBenchmarks.SUITE)
tune!(BaseBenchmarks.SUITE; verbose = true)
BenchmarkTools.save(BaseBenchmarks.PARAMS_PATH, params(BaseBenchmarks.SUITE))


text:
"""
    AbstractConnectStrategy

Base type for strategies for how to (re)connect vertices when doing structural
mutation.
"""
abstract type AbstractConnectStrategy end
struct ConnectAll <: AbstractConnectStrategy end
struct ConnectNone <: AbstractConnectStrategy end

"""
    AbstractAlignSizeStrategy

Base type for strategies for how to align size ([`nin`](@ref)/[`nout`](@ref)) when doing structural mutation.

Note that all strategies are not guaranteed to work in all cases.

Default strategies should however be selected based on case so that things always work out.
"""
abstract type AbstractAlignSizeStrategy end

"""
    NoSizeChange <: AbstractAlignSizeStrategy
    NoSizeChange()

Don't do any size change.
"""
struct NoSizeChange <: AbstractAlignSizeStrategy end

"""
    ChangeNinOfOutputs <: AbstractAlignSizeStrategy
    ChangeNinOfOutputs(Δoutsize)

Just sets [`nin`](@ref) of each output to the provided utility. Sometimes you just know the answer...
"""
struct ChangeNinOfOutputs <: AbstractAlignSizeStrategy
    Δoutsize::Int
end

"""
    FailAlignSizeNoOp <: AbstractAlignSizeStrategy
    FailAlignSizeNoOp()

Don't do any size change and return failure status.

Note that this means that graphs will most likely be left corrupted state if used as a fallback.
"""
struct FailAlignSizeNoOp <: AbstractAlignSizeStrategy end

"""
    SizeAlignFailError <: Exception

Sizes could not be aligned.
"""
struct SizeAlignFailError <: Exception
    msg::String
end
Base.showerror(io::IO, e::SizeAlignFailError) = print(io, "SizeAlignFailError: ", e.msg)

"""
    FailAlignSizeError <: AbstractAlignSizeStrategy
    FailAlignSizeError()

Throws `SizeAlignFailError`.
"""
struct FailAlignSizeError <: AbstractAlignSizeStrategy end

"""
    FailAlignSizeWarn <: AbstractAlignSizeStrategy
    FailAlignSizeWarn()
    FailAlignSizeWarn(;andthen, msgfun)

Logs warning and then proceeds with the next action.
"""
struct FailAlignSizeWarn{S, F} <: AbstractAlignSizeStrategy
    andthen::S
    msgfun::F
end
FailAlignSizeWarn(;andthen=FailAlignSizeNoOp(), msgfun=(vin,vout) -> "Could not align sizes of $(vin) and $(vout)!") = FailAlignSizeWarn(andthen, msgfun)

const doc_mapstrat = """
`mapstrat` can be used to wrap the [`AbstractΔSizeStrategy`](@ref) in a [`DecoratingJuMPΔSizeStrategy`](@ref). Main intended usecase is [`WithUtilityFun`](@ref).
"""

"""
    AlignSizeBoth <: AbstractAlignSizeStrategy
    AlignSizeBoth()
    AlignSizeBoth(;fallback, mapstrat)

Align sizes by changing both input and output.
Fallback to another strategy (default [`FailAlignSizeWarn`](@ref)) if size change is not possible.

$doc_mapstrat
"""
struct AlignSizeBoth{S, F} <: AbstractAlignSizeStrategy
    fallback::S
    mapstrat::F
end
AlignSizeBoth(;fallback=FailAlignSizeWarn(),mapstrat=identity) = AlignSizeBoth(fallback, mapstrat)

"""
    DecreaseBigger <: AbstractAlignSizeStrategy
    DecreaseBigger()
    DecreaseBigger(;fallback, mapstrat)

Try to align size by decreasing in the direction (in/out) which has the bigger size.
Fallback to another strategy (default [`AlignSizeBoth`](@ref)) if size change is not possible.

$doc_mapstrat
"""
struct DecreaseBigger{S, F} <: AbstractAlignSizeStrategy
    fallback::S
    mapstrat::F
end
DecreaseBigger(;fallback=AlignSizeBoth(),mapstrat=identity) = DecreaseBigger(fallback, mapstrat)

"""
    IncreaseSmaller <: AbstractAlignSizeStrategy
    IncreaseSmaller()
    IncreaseSmaller(;fallback, mapstrat)

Try to align size by increasing in the direction (in/out) which has the smaller size.
Fallback to another strategy (default [`DecreaseBigger`](@ref)) if size change is not possible.

$doc_mapstrat
"""
struct IncreaseSmaller{S,F} <: AbstractAlignSizeStrategy
    fallback::S
    mapstrat::F
end
IncreaseSmaller(;fallback=DecreaseBigger(),mapstrat=identity) = IncreaseSmaller(fallback, mapstrat)

"""
    CheckNoSizeCycle <: AbstractAlignSizeStrategy
    CheckNoSizeCycle()
    CheckNoSizeCycle(;ifok, ifnok)

Check if a size change in one direction causes a change in the other direction and execute
strategy `ifnok` (default [`FailAlignSizeWarn`](@ref)) if this is the case.
Motivation is that removing will result in the computation graph being in an invalid state
as one of the vertices must fulfill the impossible criterion
[`nout(v)`](@ref) `==` [`nout(v)`](@ref) `+ a` where `a > 0`.

If no such cycle is detected, then proceed to execute strategy `ifok` (default [`IncreaseSmaller`](@ref)).

Will execute strategy `ifok` if vertex shall not to be removed.
"""
struct CheckNoSizeCycle{S1,S2} <: AbstractAlignSizeStrategy
    ifok::S1
    ifnok::S2
end
CheckNoSizeCycle(;ifok=IncreaseSmaller(), ifnok=FailAlignSizeWarn(msgfun = (vin,vout) -> "Can not remove vertex $(vin)! Size cycle detected!")) = CheckNoSizeCycle(ifok, ifnok)

"""
    CheckCreateEdgeNoSizeCycle <: AbstractAlignSizeStrategy
    CheckCreateEdgeNoSizeCycle()
    CheckCreateEdgeNoSizeCycle(;ifok, ifnok)

Check if adding an edge creates the same type of size cycle that [`CheckNoSizeCycle`](@ref) checks for
and execute `ifnok` (default [`FailAlignSizeWarn`](@ref)) if this is the case.
Motivation is that removing will result in the computation graph being in an invalid state as one of
the vertices must fulfill the impossible criterion [`nout(v)`](@ref) `==` [`nout(v)`](@ref) `+ a` where `a > 0`.

If no such cycle is detected, then proceed to execute strategy `ifok` (default [`IncreaseSmaller`](@ref)).

Will check both at `prealignsizes` (i.e before edge is added) and at `postalignsizes` (i.e after edge is added).
"""
struct CheckCreateEdgeNoSizeCycle{S1, S2} <: AbstractAlignSizeStrategy
    ifok::S1
    ifnok::S2
end
CheckCreateEdgeNoSizeCycle(;ifok=IncreaseSmaller(), ifnok=FailAlignSizeWarn(msgfun = (vin,vout) -> "Can not add edge between $(vin) and $(vout)! Size cycle detected!")) = CheckCreateEdgeNoSizeCycle(ifok, ifnok)

"""
    CheckAligned <:AbstractAlignSizeStrategy
    CheckAligned()
    CheckAligned(ifnot)

Check if sizes are already aligned before making a change and return "go ahead" (`true`) if this is the case.
If not, proceed to execute another strategy (default [`CheckNoSizeCycle`](@ref)).
"""
struct CheckAligned{S} <:AbstractAlignSizeStrategy
    ifnot::S
end
CheckAligned() = CheckAligned(CheckNoSizeCycle())

"""
    PostAlign <: AbstractAlignSizeStrategy
    PostAlign()
    PostAlign(s::AbstractAlignSizeStrategy)
    PostAlign(s::AbstractAlignSizeStrategy, fallback)

Align sizes using a [`AbstractΔSizeStrategy`](@ref).

This is a post-align strategy, i.e it will be applied after a structural change has been made.
"""
struct PostAlign{S,F} <: AbstractAlignSizeStrategy
    sizestrat::S
    fallback::F
end
PostAlign() = PostAlign(DefaultJuMPΔSizeStrategy())
PostAlign(s::AbstractJuMPΔSizeStrategy; fallback = FailAlignSizeWarn()) = PostAlign(AlignNinToNout(s, ΔSizeFailNoOp()), fallback)
PostAlign(s::AlignNinToNout; fallback=FailAlignSizeWarn()) = PostAlign(s, fallback)
PostAlign(s; fallback=FailAlignSizeWarn()) = PostAlign(s, fallback)

"""
    RemoveStrategy
    RemoveStrategy()
    RemoveStrategy(rs::AbstractConnectStrategy)
    RemoveStrategy(as::AbstractAlignSizeStrategy)
    RemoveStrategy(rs::AbstractConnectStrategy, as::AbstractAlignSizeStrategy)

Strategy for removal of a vertex.

Consists of an [`AbstractConnectStrategy`](@ref) for how to treat inputs and outputs of
the removed vertex and an [`AbstractAlignSizeStrategy`](@ref) for how to align sizes of
inputs and outputs.
"""
struct RemoveStrategy{SC<:AbstractConnectStrategy, SA<:AbstractAlignSizeStrategy}
    reconnect::SC
    align::SA
end
RemoveStrategy() = RemoveStrategy(ConnectAll(), CheckAligned())
RemoveStrategy(rs::AbstractConnectStrategy) = RemoveStrategy(rs, CheckAligned())
RemoveStrategy(as::AbstractAlignSizeStrategy) = RemoveStrategy(ConnectAll(), as)

"""
    remove!(v::MutationVertex, strategy=RemoveStrategy())

Removes `v` from the graph by removing it from its [`inputs`](@ref) and [`outputs`](@ref).

It is possible to supply a strategy for how to 1) reconnect the inputs and outputs
of `v` and 2) align the input and output sizes of the `inputs` and `outputs` of `v`.

Default strategy is to first set `nin==nout` for `v` and then connect all its
[`inputs`](@ref) to all its [`outputs`](@ref).
"""
function remove!(v::MutationVertex, strategy=RemoveStrategy())
    prealignsizes(strategy.align, v, vx -> vx == v) || return false

    undo_fun = remove_with_undo!(v, strategy)

    # Note, success from prealignsizes strategies is implicit as we return if it fails
    success = postalignsizes(strategy.align, v)

    if !success
        undo_fun()
    end
    return success
end

remove_with_undo!(v) = remove_with_undo!(v, RemoveStrategy(ConnectNone()))
function remove_with_undo!(v, strategy)
    undo_fun = create_remove_undo(strategy.align, strategy.reconnect, v)

    remove!(v, inputs, outputs, strategy.reconnect)
    remove!(v, outputs, inputs, strategy.reconnect)
    return undo_fun
end

function create_remove_undo(sa, sc, v)
    undo_rm_inputs = create_remove_undo_one_direction(v, inputs, outputs)
    undo_rm_outputs = create_remove_undo_one_direction(v, outputs, inputs)
    return function()
        undo_rm_inputs()
        undo_rm_outputs()
    end
end

## Helper function to avoid code duplication. I don't expect this to be able to do
## anything useful unless f1=inputs and f2=outputs or vise versa.
function create_remove_undo_one_direction(v::AbstractVertex, f1, f2)
    prev = copy(f1(v))
    prevf1f2 = Dict(vf1 => copy(f2(vf1)) for vf1 in prev)
    return function()
        set_edges!(f1, v => prev)
        foreach(vf1f2 -> set_edges!(f2, vf1f2), prevf1f2)
    end
end

function set_edges!(f, (v, newins)::Pair)
    empty!(f(v))
    append!(f(v), newins)
end

## Helper function to avoid code duplication. I don't expect this to be able to do
## anything useful unless f1=inputs and f2=outputs or vise versa.
function remove!(v::AbstractVertex, f1::Function, f2::Function, s::AbstractConnectStrategy)
    for f1_v in f1(v)
        f1_f2_v = f2(f1_v)
        inds = findall(vx -> vx == v, f1_f2_v)
        deleteat!(f1_f2_v, inds)
        # Connect f2 of v to f2 of v1, e.g connect the inputs to v as inputs to v1 if f2 == inputs
        connect!(s, f1_f2_v, inds, f2(v))
    end
end

"""
    connect!(v, to, inds, items, ::AbstractConnectStrategy)

Does connection of `items` to `to` at position `inds` depending on the given strategy.
"""
#Not sure broadcasting insert! like this is supposed to work, but it does...
connect!(::ConnectAll, to, inds, items) = insert!.(Ref(to), inds, items)
function connect!(::ConnectNone, args...) end

tot_nin(v::AbstractVertex) = tot_nin(trait(v), v)
tot_nin(t::DecoratingTrait, v) = tot_nin(base(t), v)
tot_nin(::SizeInvariant, v) = length(unique(nin(v))) == 1 ? unique(nin(v))[] : nothing
tot_nin(::SizeTransparent, v) = sum(nin(v))
tot_nin(t::MutationTrait, v) = tot_nin(t, base(v))
tot_nin(::MutationTrait, v::CompVertex) = tot_nin(v.computation)
tot_nin(::MutationTrait, v::InputSizeVertex) = nothing
tot_nin(f) = nin(f)[1]

# Boilerplate
prealignsizes(s::AbstractAlignSizeStrategy, v, will_rm::Function) = prealignsizes(s, v, v, will_rm)
prealignsizes(s::AbstractAlignSizeStrategy, vin, vout, will_rm) = true

# Failure cases
prealignsizes(::FailAlignSizeError, vin, vout, will_rm) = throw(SizeAlignFailError("Could not align sizes of $(nameorrepr(vin)) and $(nameorrepr(vout))!"))
function prealignsizes(s::FailAlignSizeWarn, vin, vout, will_rm)
    @warn s.msgfun(vin, vout)
    return prealignsizes(s.andthen, vin, vout, will_rm)
end
prealignsizes(::FailAlignSizeNoOp, vin, vout, will_rm) = false

# Actual actions
function prealignsizes(s::CheckAligned, vin, vout, will_rm)
    nout(vin) == tot_nin(vout) && return true
    return prealignsizes(s.ifnot, vin, vout, will_rm)
end

function prealignsizes(s::CheckNoSizeCycle, vin, vout, will_rm)
    if will_rm(vout) && vin===vout
        isinsizecycle(vin) && return prealignsizes(CheckAligned(s.ifnok), vin, vout, will_rm)
    end
    return prealignsizes(s.ifok, vin, vout, will_rm)
end

function prealignsizes(s::CheckCreateEdgeNoSizeCycle, vin, vout, will_rm)
    if vout in all_in_Δsize_graph(vin, Output())
        # Previous alg created a MetaGraph of all_in_Δsize_graph and added an edge from vin to vout
        # but it didn't use that graph for anything. Should probably add a temporary edge before this....
        isinsizecycle(vin) && return prealignsizes(CheckAligned(s.ifnok), vin, vout, will_rm)
    end
    return prealignsizes(s.ifok, vin, vout, will_rm)
end

function prealignsizes(s::ChangeNinOfOutputs, vin, vout, will_rm)
    expected = nout(vin) + s.Δoutsize
    Δsize!(ΔNout{Exact}(vin, s.Δoutsize;fallback=ΔSizeFailNoOp()), all_in_Δsize_graph(vin, Output()))
    return nout(vin) == expected
end

function prealignsizes(s::Union{IncreaseSmaller, DecreaseBigger}, vin, vout, will_rm)
    Δinsize = nout(vin) - tot_nin(vout)
    Δoutsize = -Δinsize

    strat = proceedwith(s, Δinsize) ? Δninstrat(vout, Δinsize) : Δnoutstrat(vin, Δoutsize)
    success = prealignsizes(s.mapstrat(strat), vin, vout)

    if !success
        return prealignsizes(s.fallback, vin, vout, will_rm)
    end
    return true
end
proceedwith(::DecreaseBigger, Δ::Integer) = Δ <= 0
proceedwith(::IncreaseSmaller, Δ::Integer) = Δ >= 0

Δninstrat(v, Δ) = Δninstrat(trait(v), v, Δ)
Δninstrat(t::DecoratingTrait, v, Δ) = Δninstrat(base(t), v, Δ)
Δninstrat(::Immutable, v, Δ) = ΔSizeFailNoOp()
Δninstrat(::SizeTransparent, v, Δ) = ΔNoutExact(v => Δ; fallback=ΔSizeFailNoOp())
function Δninstrat(::MutationSizeTrait, v, Δ)
    Δs = Δnin2Δnout(v => Δ)
    Δs[v] = 0
    ΔNoutExact(Δs;fallback=ΔSizeFailNoOp())
end

Δnoutstrat(v, Δ) = Δnoutstrat(trait(v), v, Δ)
Δnoutstrat(t::DecoratingTrait, v, Δ) = Δnoutstrat(base(t), v, Δ)
Δnoutstrat(::Immutable, v, Δ) = ΔSizeFailNoOp()
Δnoutstrat(::SizeTransparent, v, Δ) = ΔNoutExact(v => Δ; fallback=ΔSizeFailNoOp())
Δnoutstrat(::MutationSizeTrait, v, Δ) = ΔNoutExact(v => Δ, (inputs(v) .=> 0)...; fallback=ΔSizeFailNoOp())

function prealignsizes(s::AlignSizeBoth, vin, vout, will_rm)

    strat = AlignNinToNoutVertices(vin, vout, 1:length(nin(vout)), AlignNinToNout(), ΔSizeFailNoOp())
    success = prealignsizes(s.mapstrat(strat), vin, vout)

    if !success
        return prealignsizes(s.fallback, vin, vout, will_rm)
    end
    return true
end

function prealignsizes(s::AbstractΔSizeStrategy, vin, vout)
    vin_all = all_in_Δsize_graph(vin, Output())
    vout_all = all_in_Δsize_graph(vout, Input())
    return Δsize!(s, union(vin_all, vout_all))
end

# Boilerplate
postalignsizes(s::AbstractAlignSizeStrategy, v) = postalignsizes(s, v, v, missing)
postalignsizes(s::AbstractAlignSizeStrategy, vin, vout, pos) = true

# Failure cases
postalignsizes(::FailAlignSizeNoOp, vin, vout, pos) = false
postalignsizes(::FailAlignSizeError, vin, vout, pos) = throw(SizeAlignFailError("Could not align sizes of $(nameorrepr(vin)) and $(nameorrepr(vout))!"))
function postalignsizes(s::FailAlignSizeWarn, vin, vout, pos)
    @warn s.msgfun(vin, vout)
    return postalignsizes(s.andthen, vin, vout, pos)
end

function postalignsizes(s::CheckCreateEdgeNoSizeCycle, vin, vout, pos)
    any(isinsizecycle, inputs(vout)) && return postalignsizes(s.ifnok, vin, vout, pos)
    return postalignsizes(s.ifok, vin, vout, pos)
end

# Ok, this one actually does something...
function postalignsizes(s::PostAlign, vin, vout, pos)
    vin_all = all_in_Δsize_graph(vin, Output())
    vout_all = all_in_Δsize_graph(vout, Input())

    success = Δsize!(s.sizestrat, union(vin_all, vout_all))
    if !success
        return postalignsizes(s.fallback, vin, vout, pos)
    end
    return success
end

"""
    insert!(vin::AbstractVertex, factory::Function, outselect::Function=identity)

Replace `vin` as input to all outputs of `vin` with vertex produced by `factory`.

The function `factory` takes `vin` as input.

Example:

```text
Before:

    vin
   / | \\
  /  |  \\
 v₁ ... vₙ

After:

    vin
     |
   vnew₁
     ⋮
   vnewₖ
   / | \\
  /  |  \\
 v₁ ... vₙ
```
Note that the connection `vin` -> `vnew₁` as well as all connections `vnewₚ -> vnewₚ₊₁` is done by `factory` in the above example.

The function `outselect` can be used to select a subset of outputs to replace (default all).

"""
function Base.insert!(vin::AbstractVertex, factory, outselect=identity)

    prevouts = copy(outselect(outputs(vin)))
    deleteat!(outputs(vin), findall(vx -> vx in prevouts, outputs(vin)))

    vnew = factory(vin)
    for vout in prevouts
        inds = findall(vx -> vx == vin, inputs(vout))
        inputs(vout)[inds] .= vnew
        push!(outputs(vnew), vout)
    end
    return true # TODO: Check sizes and try to align if not ok?
end

"""
    create_edge!(from::AbstractVertex, to::AbstractVertex; [pos], [strategy])

Create and edge from `from` to `to` at index `pos` in [`inputs(to)`](@ref).

Sizes will be adjusted based on given `strategy`.
"""
function create_edge!(from::AbstractVertex, to::AbstractVertex; pos = length(inputs(to))+1, strategy = default_create_edge_strat(to))

    prealignsizes(strategy, from, to, v -> false) || return false

    push!(outputs(from), to) # Order should never matter for outputs
    insert!(inputs(to), pos, from)

    # Note, success from prealignsizes strategies is implicit as we return if it fails
    success = postalignsizes(strategy, from, to, pos)

    # Revert what we did if failed
    if !success
        deleteat!(outputs(from), length(outputs(from)))
        deleteat!(inputs(to), pos)
    end
    return success
end

default_create_edge_strat(v::AbstractVertex) = default_create_edge_strat(trait(v),v)
default_create_edge_strat(t::DecoratingTrait,v) = default_create_edge_strat(base(t),v)
default_create_edge_strat(::SizeInvariant,v) = CheckCreateEdgeNoSizeCycle(ifok=AlignSizeBoth())
default_create_edge_strat(::SizeAbsorb,v) = NoSizeChange()
function default_create_edge_strat(::SizeStack,v)
    alignstrat = TruncateInIndsToValid(AlignNinToNout(DefaultJuMPΔSizeStrategy(), ΔSizeFailNoOp()))
    CheckCreateEdgeNoSizeCycle(ifok=PostAlign(alignstrat))
end

"""
    remove_edge!(from::AbstractVertex, to::AbstractVertex; [nr], [strategy])

Remove edge from `from` to `to`.

If there are multiple edges from `from` to `to` then `nr` can be used to distinguish which one shall be removed.

Sizes will be adjusted based on given `strategy`.
"""
function remove_edge!(from::AbstractVertex, to::AbstractVertex; nr = 1, strategy = default_remove_edge_strat(to))

    prealignsizes(strategy, from, to, v -> false) || return false

    in_ind  = findall(vx -> vx == from, inputs(to))[nr]
    out_ind = findall(vx -> vx == to, outputs(from))[1] #which output should not matter

    deleteat!(inputs(to), in_ind)
    deleteat!(outputs(from), out_ind)

    # Note, success from prealignsizes strategies is implicit as we return if it fails
    success = postalignsizes(strategy, from, to, in_ind)

    # Revert what we did if failed
    if !success
        insert!(inputs(to), in_ind, from)
        insert!(outputs(from), out_ind, to)
    end
    return success
end

default_remove_edge_strat(v::AbstractVertex) = default_remove_edge_strat(trait(v),v)
default_remove_edge_strat(t::DecoratingTrait,v) = default_remove_edge_strat(base(t),v)
default_remove_edge_strat(::SizeStack,v) = PostAlign()
default_remove_edge_strat(::SizeInvariant,v) = NoSizeChange()
default_remove_edge_strat(::SizeAbsorb,v) = NoSizeChange()


text:
# --------------------------------------------------------------------------
# ACE.jl and SHIPs.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

@testset "Basis Orthogonality"  begin

##

using Test
using ACE, JuLIP, JuLIP.Testing, QuadGK, LinearAlgebra
using ACE: TransformedJacobi, transform, transform_d, alloc_B, alloc_temp
using JuLIP: evaluate!

##

@info("Testing ortho-normality of Jacobi Polynomials")

α, β = 1 + rand(), 1 + rand()
N = 30
J = Jacobi(α, β, N)
Bj = alloc_B(J)
integrandJ = let J=J, B=Bj, α=α, β=β
   x ->  begin
            evaluate!(B, nothing, J, x)
            return B*B' * (1 - x)^α * (1+x)^β
         end
   end

Q = quadgk(integrandJ, -1, 1)[1]
println(@test round.(Q, digits=8) == Matrix(I, N+1, N+1))

##

@info("Testing ortho-normality of r-basis")

N = 20
rl, ru = 0.5, 3.0
fcut =  PolyCutoff2s(2, rl, ru)
trans = PolyTransform(2, 1.0)
P = TransformedJacobi(N, trans, fcut)
B = alloc_B(P)
tmp = alloc_temp(P)

integrand = let P = P, B = B, tmp = tmp
   r ->  begin
            evaluate!(B, tmp, P, r)
            return B * B' * abs(transform_d(P.trans, r))
         end
   end

Q = quadgk(integrand, rl, ru)[1]
println(@test round.(Q, digits=8) == Matrix(I, N+1, N+1))

##

@info("Testing (near-) orthonormality of r-basis via sampling")

# P = TransformedJacobi(...) from previous cell
Nsamples = 100_000

G = let
   G = zeros(length(P), length(P))
   for n = 1:Nsamples
      evaluate!(B, tmp, P, ACE.Utils.rand_radial(P))
      G += B * B'
   end
   G
end

println(@test cond(G) < 1.1)

##

@info("Testing (near-)orthonormality of Ylm-basis via sampling")

SH = ACE.SphericalHarmonics.SHBasis(5)

function gramian(SH::ACE.SphericalHarmonics.SHBasis, Nsamples=100_000)
   lenY = length(SH)
   G = zeros(ComplexF64, lenY, lenY)
   Y = alloc_B(SH)
   tmp = alloc_temp(SH)
   for n = 1:Nsamples
      evaluate!(Y, tmp, SH, ACE.Utils.rand_sphere())
      for i = 1:lenY, j = 1:lenY
         G[i,j] += Y[i] * Y[j]'
      end
   end
   return G / Nsamples
end

G = gramian(SH)
println(@test cond(G) < 1.1)

##

@info("Testing (near-)orthonormality of A-basis via sampling")

aceB = SHIPBasis( SparseSHIP(3, 5), trans, fcut )
function evalA(aceB, tmp, Rs)
   Zs = zeros(Int16, length(Rs))
   ACE.precompute_A!(tmp, aceB, Rs, Zs, 1)
   return tmp.A[1]
end

function A_gramian(aceB, Nsamples = 100_000)
   tmp = alloc_temp(aceB)
   lenA = length(tmp.A[1])
   G = zeros(ComplexF64, lenA, lenA)
   for n = 1:Nsamples
      R = ACE.rand_vec(aceB.J)
      A = evalA(aceB, tmp, [R])
      for i = 1:lenA, j = 1:lenA
         G[i,j] +=  A[i] * A[j]'
      end
   end
   return G
end

G = A_gramian(aceB)
println(@test cond(G) < 1.2)

end


text:
using StaticCompiler
using Test
using Libdl
using LinearAlgebra
using LoopVectorization
using ManualMemory
using Distributed
using StrideArraysCore
addprocs(1)
@everywhere using StaticCompiler, StrideArraysCore

remote_load_call(path, args...) = fetch(@spawnat 2 load_function(path)(args...))

@testset "Basics" begin

    simple_sum(x) = x + one(typeof(x))

    # This probably needs a macro
    for T ∈ (Int, Float64, Int32, Float32, Int16, Float16)
        _, path, = compile(simple_sum, (T,))
        @test remote_load_call(path, T(1)) == T(2)
    end
end

fib(n) = n <= 1 ? n : fib(n - 1) + fib(n - 2) # This needs to be defined globally due to https://github.com/JuliaLang/julia/issues/40990

@testset "Recursion" begin
    _, path = compile(fib, (Int,))
    @test remote_load_call(path, 10) == fib(10)

    # Trick to work around #40990
    _fib2(_fib2, n) = n <= 1 ? n : _fib2(_fib2, n-1) + _fib2(_fib2, n-2)
    fib2(n) = _fib2(_fib2, n)

    _, path = compile(fib2, (Int,))
    @test remote_load_call(path, 20) == fib(20)
    #@test compile(fib2, (Int,))[1](20) == fib(20)
end

@testset "Loops" begin
    function sum_first_N_int(N)
        s = 0
        for a in 1:N
            s += a
        end
        s
    end
    _, path = compile(sum_first_N_int, (Int,))
    @test remote_load_call(path, 10) == 55

    function sum_first_N_float64(N)
        s = Float64(0)
        for a in 1:N
            s += Float64(a)
        end
        s
    end
    _, path = compile(sum_first_N_float64, (Int,))
    @test remote_load_call(path, 10) == 55.

    function sum_first_N_int_inbounds(N)
        s = 0
        @inbounds for a in 1:N
            s += a
        end
        s
    end
    _, path = compile(sum_first_N_int_inbounds, (Int,))
    @test remote_load_call(path, 10) == 55

    function sum_first_N_float64_inbounds(N)
        s = Float64(0)
        @inbounds for a in 1:N
            s += Float64(a)
        end
        s
    end
    _, path = compile(sum_first_N_float64_inbounds, (Int,))
    @test remote_load_call(path, 10) == 55.
end

# Arrays with different input types Int32, Int64, Float32, Float64, Complex?
@testset "Arrays" begin
    function array_sum(n, A)
        s = zero(eltype(A))
        for i in 1:n
            s += A[i]
        end
        s
    end
    for T ∈ (Int, Complex{Float32}, Complex{Float64})
        _, path = compile(array_sum, (Int, Vector{T}))
        @test remote_load_call(path, 10, T.(1:10)) == T(55)
    end
end

@testset "Array allocations" begin
    function f(N)
        v = Vector{Float64}(undef, N)
        for i ∈ eachindex(v)
            v[i] = i*i
        end
        v
    end
    _, path = compile(f, (Int,))
    @test remote_load_call(path, 5) == [1.0, 4.0, 9.0, 16.0, 25.0]
end

# This is also a good test of loading and storing from the same object
@testset "Load & Store Same object" begin
    global const x = Ref(0)
    counter() = x[] += 1
    _, path = compile(counter, ())
    @spawnat 2 global counter = load_function(path)
    @test fetch(@spawnat 2 counter()) == 1
    @test fetch(@spawnat 2 counter()) == 2
end

# This is also a good test of loading and storing from the same object
counter = let x = Ref(0)
    () -> x[] += 1
end
@testset "Closures" begin
    #this currently segfaults during compilation
    @test_skip begin
        _, path = compile(counter, ())
        @spawnat 2 global counter_comp = load_function(path)
        @test fetch(@spawnat 2 counter_comp()) == 1
        @test fetch(@spawnat 2 counter_comp()) == 2
    end
end

@testset "Error handling" begin
    _, path = compile(sqrt, (Int,))
    tsk = @spawnat 2 begin
        try
            load_function(path)(-1)
        catch e;
            e
        end
    end
    @test fetch(tsk) isa DomainError
end

# Julia wants to treat Tuple (and other things like it) as plain bits, but LLVM wants to treat it as something with a pointer.
# We need to be careful to not send, nor receive an unwrapped Tuple to a compiled function.
# The interface made in `compile` should handle this fine.
@testset "Send and receive Tuple" begin
    foo(u::Tuple) = 2 .* reverse(u) .- 1

    _, path = compile(foo, (NTuple{3, Int},))
    @test remote_load_call(path, (1, 2, 3)) == (5, 3, 1)
end

# Just to call external libraries
@testset "BLAS" begin
    function mydot(a::Vector{Float64})
        N = length(a)
        BLAS.dot(N, a, 1, a, 1)
    end
    a = [1.0, 2.0]

    mydot_compiled, path = compile(mydot, (Vector{Float64},))
    # Works locally for me, but not on CI. Need some improvements to pointer relocation to be robust.
    @test_skip remote_load_call(path, a) == 5.0
    @test mydot_compiled(a) ≈ 5.0

    # This will need some more work apparently
    @test_skip begin
        _, path = compile((*), (Matrix{Float64}, Matrix{Float64}))
        A, B = rand(10, 11), rand(11, 12)
        @test remote_load_call(path, A, B) ≈ A * B
    end
end

@testset "Strings" begin
    function hello(name)
        "Hello, " * name * "!"
    end
    hello_compiled, path = compile(hello, (String,))
    @test remote_load_call(path, "world") == "Hello, world!"
end

@testset "Hello World" begin
    function hello(N)
        println("Hello World $N")
        N
    end
    # We'll need to be able to relocate a bunch of UV stuff for this, and deal with dynamic dispatch.
    @test_skip begin
        hello_compiled, path = compile(hello, (Int,))
        @test_skip remote_load_call(path, 1) == 1
    end
end

# This is a trick to get stack allocated arrays inside a function body (so long as they don't escape).
# This lets us have intermediate, mutable stack allocated arrays inside our
@testset "Alloca" begin
    function f(N)
        # this can hold at most 100 Int values, if you use it for more, you'll segfault
        buf = ManualMemory.MemoryBuffer{100, Int}(undef)
        GC.@preserve buf begin
            # wrap the first N values in a PtrArray
            arr = PtrArray(pointer(buf), (N,))
            arr .= 1 # mutate the array to be all 1s
            sum(arr) # compute the sum. It is very imporatant that no references to arr escape the function body
        end
    end
    _, path = compile(f, (Int,))
    @test remote_load_call(path, 20) == 20
end

# I can't beleive this works.
@testset "LoopVectorization" begin
    function mul!(C, A, B)
        # note: @tturbo does NOT work
        @turbo for n ∈ indices((C,B), 2), m ∈ indices((C,A), 1)
            Cmn = zero(eltype(C))
            for k ∈ indices((A,B), (2,1))
                Cmn += A[m,k] * B[k,n]
            end
            C[m,n] = Cmn
        end
    end

    C = Array{Float64}(undef, 10, 12)
    A = rand(10, 11)
    B = rand(11, 12)

    _, path = compile(mul!, (Matrix{Float64}, Matrix{Float64}, Matrix{Float64},))
    # remote_load_call(path, C, A, B) This won't work because @spawnat copies C
    C .= fetch(@spawnat 2 (load_function(path)(C, A, B); C))
    @test C ≈ A*B
end

@testset "Standalone Dylibs" begin
    # Test function
    # (already defined)
    # fib(n) = n <= 1 ? n : fib(n - 1) + fib(n - 2)

    #Compile dylib
    name = repr(fib)
    filepath = compile_shlib(fib, (Int,), "./", name)
    @test occursin("fib.$(Libdl.dlext)", filepath)

    # Open dylib
    ptr = Libdl.dlopen(filepath, Libdl.RTLD_LOCAL)
    fptr = Libdl.dlsym(ptr, "julia_$name")
    @test fptr != C_NULL
    @test ccall(fptr, Int, (Int,), 10) == 55
end

@testset "Standalone Executables" begin
    # Minimal test with no `llvmcall`
    @inline function foo()
        v = 0.0
        n = 1000
        for i=1:n
            v += sqrt(n)
        end
        return 0
    end

    filepath = compile_executable(foo, (), tempdir())

    r = run(`$filepath`);
    @test isa(r, Base.Process)
    @test r.exitcode == 0

    @static if VERSION>v"1.8.0-DEV" # The llvmcall here only works on 1.8+
        @inline function puts(s::Ptr{UInt8}) # Can't use Base.println because it allocates
            Base.llvmcall(("""
            ; External declaration of the puts function
            declare i32 @puts(i8* nocapture) nounwind

            define i32 @main(i8*) {
            entry:
               %call = call i32 (i8*) @puts(i8* %0)
               ret i32 0
            }
            """, "main"), Int32, Tuple{Ptr{UInt8}}, s)
        end

        @inline function print_args(argc::Int, argv::Ptr{Ptr{UInt8}})
            for i=1:argc
                # Get pointer
                p = unsafe_load(argv, i)
                # Print string at pointer location (which fortunately already exists isn't tracked by the GC)
                puts(p)
            end
            return 0
        end

        filepath = compile_executable(print_args, (Int, Ptr{Ptr{UInt8}}), tempdir())

        r = run(`$filepath Hello, world!`);
        @test isa(r, Base.Process)
        @test r.exitcode == 0
    end
end

# data structures, dictionaries, tuples, named tuples


text: # generic helper functions


text:
using Documenter, OhMyREPL

makedocs(
    sitename = "OhMyREPL",
    pages = Any[
        "Home" => "index.md",
        "Installation" => "installation.md",
        "Features" => Any[
            "features/syntax_highlighting.md",
            "features/bracket_highlighting.md",
            "features/bracket_complete.md",
            "features/prompt.md",
            "features/rainbow_brackets.md",
            "features/markdown_highlight.md",
            ],
        "Internals" => Any[
            "internals/passes.md"
            ]
    ]
)

deploydocs(
    repo = "github.com/KristofferC/OhMyREPL.jl.git",
)


text:
# detection (aka segmentation) - detecting sources on images

import WCS
import .SEP
using .Coordinates: match_coordinates
using .Model

# return the world coordinates of all objects in the catalog as a 2xN matrix
# where `result[:, i]` is the (latitude, longitude) of the i-th object.
function _worldcoords(catalog::SEP.Catalog, wcs::WCS.WCSTransform)
    pixcoords = Array{Float64}(2, length(catalog.x))
    for i in eachindex(catalog.x)
        pixcoords[1, i] = catalog.x[i]
        pixcoords[2, i] = catalog.y[i]
    end
    return WCS.pix_to_world(wcs, pixcoords)
end

# Get angle offset between +x axis and +Dec axis (North) from a WCS transform.
# This assumes there is no skew, meaning the x and y axes are perpindicular in
# world coordinates. It is also only based on the CD matrix.
function _x_vs_n_angle(wcs::WCS.WCSTransform)
    cd = wcs[:cd]::Matrix{Float64}
    sgn = sign(det(cd))
    n_vs_y_rot = atan2(sgn * cd[1, 2],  sgn * cd[1, 1])  # angle of N CCW
                                                         # from +y axis
    return -(n_vs_y_rot + pi/2)  # angle of +x CCW from N
end

"""
    detect_sources(images)

Detect sources in a set of (possibly overlapping) `Image`s and combine
duplicates. Returns a catalog and a 2-d array of `ImagePatch`.

"""
function detect_sources(images::Vector{<:Image})

    catalogs = SEP.Catalog[]
    for image in images

        calpixels = Model.calibrated_pixels(image)

        # Run background, just to get background rms.
        #
        # We're using sky subtracted (and calibrated) image data, but,
        # we still run a background analysis, just to determine the
        # rough image noise for the purposes of setting a
        # threshold. This could be slightly suboptimal in terms of
        # selecting faint sources: If the image noise varies
        # significantly across the image, the threshold might be too
        # high or too low in some places. However, we don't really
        # trust the variable background RMS without first masking
        # sources. (We could add this.)
        bkg = SEP.Background(calpixels; boxsize=(256, 256),
                             filtersize=(3, 3))
        noise = SEP.global_rms(bkg)
        push!(catalogs, SEP.extract(calpixels, 1.3; noise=noise))
    end

    # The image catalogs only have pixel positions. To match objects in
    # different images, we need to convert these pixel positions to world
    # coordinates in each image.
    catalogs_worldcoords = [_worldcoords(catalog, image.wcs)
                            for (catalog, image) in zip(catalogs, images)]

    # Initialize the "joined" catalog to all the coordinates in the first
    # image. `detections` tracks image index and object index of detections.
    joined_ra = length(images) > 0 ? catalogs_worldcoords[1][1, :] : Float64[]
    joined_dec = length(images) > 0 ? catalogs_worldcoords[1][2, :] : Float64[]
    detections = (length(images) > 0 ?
                  [[(1, j)] for j in 1:length(catalogs[1].x)] :
                  Vector{Tuple{Int, Int}}[])

    # Search for matches in remaining images
    for i in 2:length(images)
        ra = catalogs_worldcoords[i][1, :]
        dec = catalogs_worldcoords[i][2, :]
        idx, dist = match_coordinates(ra, dec, joined_ra, joined_dec)
        for j in eachindex(idx)
            # If there is an object in the joined catalog within 1 arcsec,
            # add the (image, object) index to detections for that object.
            # Otherwise, add the position and (image, object) index as a new
            # object in the joined catalog.
            if dist[j] < (1.0 / 3600.0)
                push!(detections[idx[j]], (i, j))
            else
                push!(joined_ra, ra[j])
                push!(joined_dec, dec[j])
                push!(detections, [(i, j)])
            end
        end
    end

    # Initialize joined output catalog and patches
    nobjects = length(joined_ra)
    result = Array{CatalogEntry}(nobjects)
    patches = Array{ImagePatch}(nobjects, length(images))

    # Precalculate some angles
    x_vs_n_angles = [_x_vs_n_angle(image.wcs) for image in images]

    # Loop over output catalog:
    # - Create catalog entry based on "best" detection.
    # - Create patches based on detection in each image.
    for i in 1:nobjects
        world_center = [joined_ra[i], joined_dec[i]]

        # which detection in each band is the "best" (has most pixels)
        # We use this to initialize flux in each band
        best = fill((0, 0), NUM_BANDS)
        npix = fill(0, NUM_BANDS)
        for (j, catidx) in detections[i]
            b = images[j].b
            np = catalogs[j].npix[catidx]
            if np > npix[b]
                best[b] = (j, catidx)
                npix[b] = np
            end
        end

        # set gal_fluxes (and star_fluxes) to best detection flux or 0 if not
        # detected.
        gal_fluxes = [j != 0 ? catalogs[j].flux[catidx] : 0.0
                      for (j, catidx) in best]
        star_fluxes = copy(gal_fluxes)

        # use the single best band for shape parameters
        j, catidx = best[indmax(npix)]
        gal_axis_ratio = catalogs[j].b[catidx] / catalogs[j].a[catidx]

        # SEP angle is CCW from +x axis and in [-pi/2, pi/2].
        # Add offset of x axis from N to make angle CCW from N
        gal_angle = catalogs[j].theta[catidx] + x_vs_n_angles[j]

        # A 2-d symmetric gaussian has CDF(r) =  1 - exp(-(r/sigma)^2/2)
        # The half-light radius is then r = sigma * sqrt(2 ln(2))
        sigma = sqrt(catalogs[j].a[catidx] * catalogs[j].b[catidx])
        gal_radius_px = sigma * sqrt(2.0 * log(2.0))

        result[i] = CatalogEntry(world_center,
                                 false,  # is_star
                                 star_fluxes,
                                 gal_fluxes,
                                 0.5,  # gal_frac_dev
                                 gal_axis_ratio,
                                 gal_angle,
                                 gal_radius_px)

        # create patches based on detections
        for (j, catidx) in detections[i]
            box = (Int(catalogs[j].xmin[catidx]):Int(catalogs[j].xmax[catidx]),
                   Int(catalogs[j].ymin[catidx]):Int(catalogs[j].ymax[catidx]))
            box = Model.dilate_box(box, 0.2)
            minbox = Model.box_around_point(images[j].wcs, world_center, 5.0)
            patches[i, j] = ImagePatch(images[j], Model.enclose_boxes(box, minbox))
        end

        # fill patches for images where there was no detection
        for j in 1:length(images)
            if !isassigned(patches, i, j)
                box = Model.box_around_point(images[j].wcs, world_center, 5.0)
                patches[i, j] = ImagePatch(images[j], box)
            end
        end
    end

    return result, patches
end


text:
using Statistics
using Formatting
using Crayons.Box
using NumericalDataManipulation.Common
using NumericalDataManipulation.CommonMath

struct DataGrid
    knots::Vector{Float64}
    Δmin::Float64
    Δmax::Float64
    Δmean::Float64
    Δlog10::Float64
    Δ𝐆::Float64
    Δknots::Vector{Float64}
    function DataGrid(knots::Vector{Float64})
        @assert(!isempty(knots), "Data grid cannot be empty.")
        Δknots, Δmin, Δmax, Δmean, Δlog10, Δ𝐆 = give_Δknots(knots)
        grid = new(knots, Δmin, Δmax, Δmean, Δlog10, Δ𝐆, Δknots)
    end
    function DataGrid(range::StepRange)
        @assert(issorted(range), "The range must be sorted in an ascending order.")
        knots = collect(range)
        Δknots, Δmin, Δmax, Δmean, Δlog10, Δ𝐆 = give_Δknots(knots)
        grid = new(knots, Δmin, Δmax, Δmean, Δlog10, Δ𝐆, Δknots)
    end
end
function Base.show(io::IO, grid::DataGrid)
    compact = get(io, :compact, false)

    knots = grid.knots; a = knots[1]; b = knots[end]; N = length(knots)
    Δmin = grid.Δmin
    Δmax = grid.Δmax
    Δmean = grid.Δmean
    Δlog10 = grid.Δlog10
    Δ𝐆 = grid.Δ𝐆
    Δknots = grid.Δknots

    if compact
        printfmt(io, "$(LIGHT_GRAY_FG("DataGrid([{1:.3e}, {2:.3e}]($N), Δlog10 = {3:.6e}, Δmin = {4:.6e}, Δmax = {5:.3e}, Δmean = {6:.5e}, Δ𝐆 = {7:.6e})"))",
            a, b, Δlog10, Δmin, Δmax, Δmean, Δ𝐆)
    else
        printfmt(io, "$(LIGHT_GRAY_FG("DataGrid([{1:.3e}, {2:.3e}]($N), Δlog10 = {3:.6e}, Δmin = {4:.6e}, Δmax = {5:.3e}, Δmean = {6:.5e}, Δ𝐆 = {7:.6e})"))",
            a, b, Δlog10, Δmin, Δmax, Δmean, Δ𝐆)
    end

end
function Base.show(io::IO, ::MIME"text/plain", grid::DataGrid)
    show(io, grid)
end

function give_Δknots(knots::Vector{Float64})
    @assert(length(knots) ≥ 3, "Grid must contain at least 3 knots.")
    @assert(issorted(knots), "The knots must be sorted in an ascending order.")
    Δknots = knots[2:end] - knots[1:end-1]
    Δmin = minimum(Δknots)
    Δmax = maximum(Δknots)
    Δmean = mean(Δknots)
    Δ𝐆 = mean_log10_golden_lower(Δmin, Δmax)
    Δlog10 = mean_log10(Δmin, Δmax)
    return Δknots, Δmin, Δmax, Δmean, Δlog10, Δ𝐆
end


text:
using Knet:RNN

"""Abstract type for RNN modules"""
abstract type AbstractRNN <: KnetModule end

for (RNNModule, rnnType) in zip([:RNNRelu, :RNNTanh, :LSTM, :GRU],
                                [:(:relu), :(:tanh), :(:lstm), :(:gru)])

    eval(:(begin
      type $(RNNModule) <: AbstractRNN
           w::Param
           cfg::RNN
           y::Bool
           hy::Bool
           cy::Bool
           _opt
      end

     rname = split(string($(RNNModule)), ".")[end]
"""
# Constructor
    $(rname)(input::Int, hidden::Int; y=true, o...)
        `y` determines whether or not to return output (state history)
        `o` might contain all keyword arguments to `Knet.rnnforw` and `Knet.rnninit`

# Fields
    `w::Param` stores trainable parameters
    `cfg::Knet.RNN` an object stores the RNN configuration
    `y::Bool` stores whether or not to return output (state history)
    `hy::Bool` stores whether or not to return next hidden state
    `cy::Bool` stores whether or not to return next cell state

# Forward execution
    r(ctx, x[, hy, cy]; o...)
where `o...` is runtime kwargs of `Knet.rnnforw`
"""
      function $(RNNModule)(input::Int, hidden::Int;
                            y=true, hy=false, cy=false,
                            o...)
           @assert (y || hy || cy) "You must return something"
           r, w = rnninit(input, hidden;
                          usegpu=false, rnnType=$(rnnType),o...)
           return $(RNNModule)(Param(w), r, y, hy, cy, o)
      end

      (rnn::$(RNNModule))(ctx, x, hx=nothing, cx=nothing; o...) =
           _forw(ctx, rnn, x, hx, cx; o...)

     end))
end

import Knet.rnnparams

"`rnnparams(r::AbstractRNN; o...)` gets parameter matrices and biases of module `r`"
rnnparams(r::AbstractRNN; o...) = rnnparams(r.cfg, aval(r.w); o...)

function _forw(ctx, rnn::AbstractRNN, x, hx, cx; o...)
    y, hy, cy, _ = rnnforw(rnn.cfg, val(ctx, rnn.w), x, hx, cx;
                           hy=rnn.hy, cy=rnn.cy, o...)
    # Post-process state
    gets = rnn.hy || rnn.cy
    next_state = (rnn.cy && cy !== nothing) ? (hy, cy) : hy
    return (gets && rnn.y) ? (y, next_state...) : gets ? next_state : y
end

function convert_params!(rnn::AbstractRNN, atype)
    r,w = rnninit(rnn.cfg.inputSize, rnn.cfg.hiddenSize;
                  usegpu=atype<:KnetArray,
                  rnn._opt...)
    rnn.cfg = r
    copy!(w, aval(rnn.w))
    setval!(rnn.w, w)
end


text:
# Lazily generate a C code so that it can be used as input to GCC.

function maybe_generate_code(size, file)
    if !isfile(file)
        command = `csmith --max-funcs $size`
        run(pipeline(command, stdout = file))
    end
end

function collect_csmith_statistics(fpath)
    found = false
    file = open(fpath, "r")

    for line in eachline(file)
        if line == "/************************ statistics *************************"
            found = true
        end
    end

    if !found
        exit("Provided file seems not to be a valid csmith file")
    end
end

# Convert the input vector to parameters that GCC can understeand, and
# launch it measuring the time taken.

function run_gcc(x)
	file_prefix = "processed/"
	compiler_prefix = "/tmp/gcc11_autopar/usr/local/bin/"

	#if !isdir(prefix)
	#    mkdir(prefix)
	#end

	file = file_prefix * x[1]
	lto_min_partition = x[2]
	run_parallel = x[3]

	file_splitted = split(x[1], ".")
	extension = file_splitted[length(file_splitted)]

	if (extension == "cc")
		compiler = compiler_prefix * "g++"
	else
		compiler = compiler_prefix * "gcc"
	end

	if (run_parallel == 1)
		parallel = ""
	else
		parallel = "-fparallel-jobs=" * string(run_parallel)
	end

	proc = "/tmp/gcc11_autopar/usr/local/bin/gcc"
	argv = ["--param=promote-statics=1", "-I/usr/include/csmith/", "-Wno-all", "--param=lto-min-partition=$lto_min_partition", "-c", "-o", "temp.o", "-O2", file]

	if parallel != ""
		append!(argv, [parallel])
	end

	command = `$proc $argv`

	t1 = time_ns()
	run(command)
	t2 = time_ns()

	rm("temp.o")

	return (t2 - t1)/1e9
end

#Quick test
#x = run_gcc(["0-0-10.cc", 1, false])
#println(x)

max_stars_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_issues_repo_name: giulianobelinassi/LTO-Timing-Analysis

max_forks_repo_name: giulianobelinassi/LTO-Timing-Analysis