text:
import Base.*, Base./, Base.+, Base.-, Base.==, Base.^, Base.exp, Base.isapprox

@doc raw"""
    tf = TransferFunction([1, 2], [3, 5, 8])
    tf = TransferFunction([1, 2], [3, 5, 8], 3.0)

construct a transfer function representing a linear, time-invariant system.

# example
to construct the transfer function

```math
G(s) = \frac{4e^{-2.2s}}{2s+1}
```

in Julia:

```jldoctest
tf = TransferFunction([4], [2, 1], 2.2)
# output
    4.0
----------- e^(-2.2*s)
2.0*s + 1.0
```

# attributes
* `numerator::Polynomial{Float64, :s}`: the polynomial in the numerator of the transfer function
* `denominator::Polynomial{Float64, :s}`: the polynomial in the denominator of the transfer function
* `time_delay::Float64`: the associated time delay
"""
struct TransferFunction
    numerator::Polynomial{Float64, :s}
    denominator::Polynomial{Float64, :s}
    time_delay::Float64
end

TransferFunction(num::Union{Vector{Int}, Vector{Float64}},
                 den::Union{Vector{Int}, Vector{Float64}},
                 td::Union{Int, Float64}) =
    TransferFunction(Polynomial(reverse(convert(Vector{Float64}, num)), :s),
                     Polynomial(reverse(convert(Vector{Float64}, den)), :s),
                     convert(Float64, td))

TransferFunction(num::Union{Vector{Int}, Vector{Float64}},
                 den::Union{Vector{Int}, Vector{Float64}}) =
    TransferFunction(num, den, 0.0)

const s = TransferFunction([1, 0], [1])

function *(tf1::TransferFunction, tf2::TransferFunction)
    g = TransferFunction(tf1.numerator * tf2.numerator,
                         tf1.denominator * tf2.denominator,
                         tf1.time_delay + tf2.time_delay)
    return pole_zero_cancellation(g)
end

*(k::Number, tf::TransferFunction) = TransferFunction(k * tf.numerator, tf.denominator, tf.time_delay)
*(tf::TransferFunction, k::Number) = *(k, tf)
function ^(tf::TransferFunction, x::Int)
    # x< 0 case handled by MethodError
    if x == 0
        return TransferFunction([1], [1])
    else
        tfp = tf
        for i = 2:x
            tfp = tfp * tf
        end
        return tfp
    end
end

"""
    tf_time_delay = exp(- 3 * s)

Conveniently create a time delay by exp(- θ * s).

# Example
g = 1 / (s + 1) * exp(-2.0 * s) # introduce time delay of 2.0
"""
function exp(tf::TransferFunction)
    # if numerator is zero, then the value is 1.0 and no time delay is introduced
    if degree(tf.numerator) == -1
        return TransferFunction([1], [1])
    end
    if ! ((degree(tf.numerator) == 1) && (degree(tf.denominator) == 0) && (tf.time_delay == 0.0))
        error("exp only functions with exp(-3.0 * s) for example, to introduce time delays")
    end
    return TransferFunction([1.0], [1.0], -tf.numerator[1])
end

function +(tf1::TransferFunction, tf2::TransferFunction)
    if (tf1.time_delay != tf2.time_delay)
        error("""we cannot add two transfer functions with two different time delays.

        if you are trying to construct a closed-loop transfer function with a time delay,
        instead use `ClosedLoopTransferFunction`.""")
    end
    g = TransferFunction(tf1.numerator * tf2.denominator + tf1.denominator * tf2.numerator,
                            tf1.denominator * tf2.denominator,
                            tf1.time_delay)
    return pole_zero_cancellation(g)
end

+(tf::TransferFunction, x::Number) = tf + TransferFunction([x], [1.0])
+(x::Number, tf::TransferFunction) = +(tf, x)
-(tf::TransferFunction, x::Number) = +(tf, -x)
-(x::Number, tf::TransferFunction) = +(-1 * tf, x)
-(tf::TransferFunction) = -1.0 * tf
-(tf1::TransferFunction, tf2::TransferFunction) = +(tf1, -1*tf2)

function /(tf1::TransferFunction, tf2::TransferFunction)
    g = TransferFunction(tf1.numerator * tf2.denominator,
                        tf1.denominator * tf2.numerator,
                        tf1.time_delay - tf2.time_delay)
    return pole_zero_cancellation(g)
end
/(tf::TransferFunction, x::Number) = TransferFunction(tf.numerator / x, tf.denominator, tf.time_delay)
/(x::Number, tf::TransferFunction) = TransferFunction([x], [1.0], 0.0) / tf

# note: must have polynomails written in same form
==(tf1::TransferFunction, tf2::TransferFunction) = (tf1.numerator == tf2.numerator) && (tf1.denominator == tf2.denominator) && (tf1.time_delay == tf2.time_delay)

@doc raw"""
    tf = zpk_form(tf)

write transfer function `tf` in zeros, poles, k-factor form:

$$g(s)=k\dfrac{\Pi_j (s-z_j)}{\Pi_j (s-p_j)}$$
where $z_j$ is zero $j$, $p_j$ is pole $j$, and $k$ is a constant factor (not equal to the zero-frequency gain) that uniquely specifies the transfer function.

this is achieved by multiplying by 1.0 in a fancy way such that the highest power of $s$ in the denominator is associated with a coefficient of $1$.

# Example

```jldoctest
g = 8.0 / (2 * s^2 + 3 * s + 4)
g_zpk = zpk_form(g)
# output
         4.0
---------------------
1.0*s^2 + 1.5*s + 2.0
```
"""
function zpk_form(tf::TransferFunction)
    coeff_highest_power_denom = tf.denominator[end]

    return TransferFunction(
        tf.numerator / coeff_highest_power_denom,
        tf.denominator / coeff_highest_power_denom,
        tf.time_delay)
end

function isapprox(tf1::TransferFunction, tf2::TransferFunction; atol=0.0)
    # to directly compare numerators and denominators, put in zpk form
    tf1s = zpk_form(tf1)
    tf2s = zpk_form(tf2)
    return isapprox(tf1s.time_delay, tf2s.time_delay, atol=atol) && isapprox(
	    tf1s.numerator.coeffs, tf2s.numerator.coeffs, atol=atol) && isapprox(tf1s.denominator.coeffs, tf2s.denominator.coeffs, atol=atol)
end

_zeros(tf::TransferFunction) = roots(tf.numerator)

_poles(tf::TransferFunction) = roots(tf.denominator)

_k(tf::TransferFunction) = tf.numerator[end] / tf.denominator[end]

@doc raw"""
    K = zero_frequency_gain(tf)

compute the (signed) zero frequency gain of a transfer function $g(s)$, which is:

```math
K := \lim_{s\rightarrow 0} G(s)
```

the zero-frequency gain "represents the ratio of the steady state value of the output with respect to a step input" [source](http://www.cds.caltech.edu/~murray/books/AM05/pdf/am06-xferfcns_16Sep06.pdf)

# example

```jldoctest
g = 5 / (3 * s + 1)
K = zero_frequency_gain(g)
# output
5.0
```

# arguments
* `tf::TransferFunction`: the transfer function

# returns
* `K::Float64`: the zero-frequency gain of the transfer function
"""
function zero_frequency_gain(tf::TransferFunction)
    tf_simplified = pole_zero_cancellation(tf)
    return evaluate(tf_simplified, 0.0)
end

@doc raw"""
    # compute the zeros, poles, and k-factor of a transfer function
    z, p, k = zeros_poles_k(tf)
    # construct a transfer function from its zeros, poles, and k-factor
    tf = zeros_poles_k(z, p, k, time_delay=0.0)

the representation of a transfer function in this context is:

$$g(s)=k\dfrac{\Pi_j (s-z_j)}{\Pi_j (s-p_j)}$$
where $z_j$ is zero $j$, $p_j$ is pole $j$, and $k$ is a constant factor (not equal to the zero-frequency gain) that uniquely specifies the transfer function.

* the zeros are the zeros of the numerator of the transfer function.
* the poles are the zeros of the denominator of the transfer function.
"""
zeros_poles_k(tf::TransferFunction) = _zeros(tf), _poles(tf), _k(tf)

# docstring above. this overloaded function is for *constructing* tfs
function zeros_poles_k(zeros::Array, poles::Array, k::Union{Float64, Int};
                       time_delay::Union{Int64, Float64}=0)
    # construct numerator polynomial
    top = length(zeros) == 0 ? Polynomial(1.0, :s) : fromroots(zeros, var=:s)

    # construct denominator polynomial
    bottom = length(poles) == 0 ? Polynomial(1.0, :s) : fromroots(poles, var=:s)

    # owing to numerical errors, the coeff's could be imaginary with tiny imaginary parts
    #  let's check that the imaginary parts are tiny, then convert the coefficients to real.
    @assert maximum(abs.(imag.(top.coeffs))) < 1e-6
    @assert maximum(abs.(imag.(bottom.coeffs))) < 1e-6
    top = Polynomial(real.(top.coeffs), :s) # reconstruct a polynomial with real coeffs
    bottom = Polynomial(real.(bottom.coeffs), :s)

    return TransferFunction(k * top, bottom, time_delay)
end

@doc raw"""
    z, p, gain = zeros_poles_gain(tf)

Compute the zeros, poles, and zero-frequency gain of a transfer function.

* the zeros are the zeros of the numerator of the transfer function.
* the poles are the zeros of the denominator of the transfer function.
* the zero-frequency gain is the transfer function evaluated at $s=0$
"""
zeros_poles_gain(tf::TransferFunction) = _zeros(tf), _poles(tf), zero_frequency_gain(tf)

"""
    evaluate(tf, z)

evaluate a `TransferFunction`, `tf`, at a particular number `z` (could be complex).

# example
```jldoctest eval
tf = TransferFunction([1], [3, 1])
evaluate(tf, 1.0)
# output
0.25
```

```jldoctest eval
```
"""
function evaluate(tf::TransferFunction, z::Number)
    return tf.numerator(z) / tf.denominator(z) * exp(-tf.time_delay * z)
end

"""
    proper(tf)

Return `true` if transfer function `tf` is proper and `false` otherwise.
"""
proper(tf::TransferFunction) = degree(tf.numerator) <= degree(tf.denominator)

"""
    strictly_proper(tf)

Return `true` if transfer function `tf` is strictly proper and `false` otherwise.
"""
strictly_proper(tf::TransferFunction) = degree(tf.numerator) < degree(tf.denominator)

"""
    tf = pole_zero_cancellation(tf, verbose=false, digits=8)

find (pole, zero) pairs such that pole = zero and return a new transfer function with those pairs cancelled.
this is achieved by comparing the poles and zeros with `isapprox`, with poles and zeros rounded to
`digits` digits (also applies to reconstruction).

# arguments
* `tf::TransferFunction`: the transfer function
* `verbose::Bool=false`: print off which poles, zeros are cancelled.
* `digits::Int`: number of digits to round poles and zeros to, for (i) cancelling and (ii) reconstruction.

# example
```jldoctest
pole_zero_cancellation(s * (s - 1) / (s * (s + 1)))
# output
1.0*s - 1.0
-----------
1.0*s + 1.0
```
"""
function pole_zero_cancellation(tf::TransferFunction; verbose::Bool=false, digits::Int=6)
    # compute poles, zeros, and k-factor of the transfer function
    zs, ps, k = zeros_poles_k(tf)
    zs = round.(zs, digits=digits)
    ps = round.(ps, digits=digits)

    # store in these boolean arrays whether we will cancel them or not
    canceled_zeros = [false for i = 1:length(zs)]
    canceled_poles = [false for i = 1:length(ps)]

    # loop through poles and zeros. if they are equal, cancel!
    for (i_z, z) in enumerate(zs)
        for (i_p, p) in enumerate(ps)
            # the pole and zero are equal...
            # *and* if this pole has not already been canceled...
	    if isapprox(p, z) && (! canceled_poles[i_p])
                canceled_zeros[i_z] = true
                canceled_poles[i_p] = true
                break
            end
        end
    end
    @assert sum(canceled_zeros) == sum(canceled_poles)
    if verbose && sum(canceled_poles) != 0
        println("canceling the following poles and zeros: ", ps[canceled_poles])
    end

    if sum(canceled_zeros) == 0
        return tf
    else
        # reconstruct
        return zeros_poles_k(zs[.! canceled_zeros], ps[.! canceled_poles], k, time_delay=tf.time_delay)
    end
end

"""
    o = system_order(tf::TransferFunction)

return the order of the numerator and denominator of the transfer function `tf`.

use [`pole_zero_cancellation`](@ref) first if you wish to cancel poles and zeros that are equal before determining the order.

# returns
`o::Tuple{Int, Int}`: (order of numerator, order of denominator)

# examples
```jldoctest
g = 1 / (s + 1)
system_order(g)
# output
(0, 1)

g = (s + 1) / ((s + 2) * (s + 3))
system_order(g)
# output
(1, 2)
```
"""
function system_order(tf::TransferFunction)
    return (degree(tf.numerator), degree(tf.denominator))
end

max_stars_repo_name: SimonEnsemble/Controlz

max_issues_repo_name: SimonEnsemble/Controlz

max_forks_repo_name: SimonEnsemble/Controlz


text:
module TestOrderings

using UnsafeAtomics: UnsafeAtomics, Ordering
using Test

const ORDERINGS = [:unordered, :monotonic, :acquire, :release, :acq_rel, :seq_cst]

function check_ordering(name::Symbol)
    ord = getfield(UnsafeAtomics, name)
    @test ord isa Ordering
    @test string(ord) == sprint(print, ord) == string(name)
    @test repr(ord) == "UnsafeAtomics.$name"
end

function test_ordering()
    @testset for name in ORDERINGS
        check_ordering(name)
    end
end

function test_aliases()
    @test UnsafeAtomics.acq_rel === UnsafeAtomics.acquire_release
    @test UnsafeAtomics.seq_cst === UnsafeAtomics.sequentially_consistent
end

end  # module

ext: jl


text:
using PkgOptions
using PkgOptionsDemo
using Test

@testset "PkgOptions.jl" begin
    @test PkgOptionsDemo.pkgoptions isa Dict
end


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# based on Oleg Mazurov's Java Implementation
# transliterated by Hamza Yusuf Çakır

global const nchunks = 150

struct Fannkuch
    n::Int64
    chunksz::Int32
    ntasks::Int32
    maxflips::Vector{Int32}
    chksums::Vector{Int32}
    taskId::Threads.Atomic{Int}

    function Fannkuch(n)
        nfact = factorial(n)

        chunksz = (nfact + nchunks - 1) ÷ nchunks
        ntasks = (nfact + chunksz - 1) ÷ chunksz

        maxflips = Vector{Int32}(undef, ntasks)
        chksums = Vector{Int32}(undef, ntasks)

        taskId = Threads.Atomic{Int}(0)

        new(n, chunksz, ntasks, maxflips, chksums, taskId)
    end
end

struct Perm
    p::Vector{Int32}
    pp::Vector{Int32}
    count::Vector{Int32}

    function Perm(n)
        p = zeros(Int32, n)
        pp = zeros(Int32, n)
        count = zeros(Int32, n)

        new(p, pp, count)
    end
end

Base.@propagate_inbounds function first_permutation(perm::Perm, idx)
    p = perm.p
    pp = perm.pp

    for i = Int32(1):Int32(length(p))
        p[i] = i - 1
    end

    for i = Int32(length(p)):Int32(-1):Int32(2)
        ifact = factorial(i-1)
        d = idx ÷ ifact
        perm.count[i] = d
        idx = idx % ifact

        for j = Int32(1):i
            pp[j] = p[j]
        end

        for j = Int32(1):i
            p[j] = j+d <= i ? pp[j+d] : pp[j+d-i]
        end
    end
end

Base.@propagate_inbounds function next_permutation(perm::Perm)
    p = perm.p
    count = perm.count

    first = p[2]
    p[2]  = p[1]
    p[1]  = first

    i = 2
    while (count[i] += 1) > i - 1
        count[i] = 0
        i += 1

        next = p[1] = p[2]

        for j = 1:i-1
            p[j] = p[j+1]
        end

        p[i] = first
        first = next
    end
end

Base.@propagate_inbounds function count_flips(perm::Perm)
    p = perm.p
    pp = perm.pp

    flips = 1
    first = p[1]

    if p[first + 1] != 0
        pp .= p

        while true # do..while(pp[first+1] != 0)
            flips += 1

            lo = 1; hi = first - 1
            while lo < hi
                t = pp[lo+1]
                pp[lo+1] = pp[hi+1]
                pp[hi+1] = t
                lo += 1
                hi -= 1
            end

            t = pp[first+1]
            pp[first+1] = first
            first = t

            (pp[first+1] == 0) && break
        end
    end

    return flips
end

Base.@propagate_inbounds function run_task(f::Fannkuch, perm::Perm, task)
    idxmin = task * f.chunksz
    idxmax = min(factorial(f.n), idxmin + f.chunksz)

    first_permutation(perm, idxmin)

    maxflips = 1
    chksum = 0

    let
        i = idxmin
        while true
            if perm.p[1] != 0
                flips = count_flips(perm)
                maxflips = max(maxflips, flips)
                chksum += iseven(i) ? flips : -flips
            end

            i += 1
            if i == idxmax
                break
            end

            next_permutation(perm)
        end
    end

    f.maxflips[task+1] = maxflips
    f.chksums[task+1] = chksum
end

function runf(f::Fannkuch)
    perm = Perm(f.n)

    taskId = f.taskId # atomic
    while (task = Threads.atomic_add!(taskId, 1)) < f.ntasks
        @inbounds run_task(f, perm, task)
    end
end

function fannkuchredux(n::Int)
    f = Fannkuch(n)

    Threads.@threads for i = 1:Threads.nthreads()
        runf(f)
    end

    chk = sum(f.chksums)
    res = maximum(f.maxflips)

    return (chk, res)
end

n = parse(Int, ARGS[1])
chk, res = fannkuchredux(n)
println(chk, "\nPfannkuchen(", n, ") = ", res)


text:
###
### LongSequence specific specializations of src/biosequence/transformations.jl
###

"""
    resize!(seq, size, [force::Bool])

Resize a biological sequence `seq`, to a given `size`. Does not resize the underlying data
array unless the new size does not fit. If `force`, always resize underlying data array.
"""
function Base.resize!(seq::LongSequence{A}, size::Integer, force::Bool=false) where {A}
    if size < 0
        throw(ArgumentError("size must be non-negative"))
    else
        if force | (seq_data_len(A, size) > seq_data_len(A, length(seq)))
            resize!(seq.data, seq_data_len(A, size))
        end
        seq.len = size
        return seq
    end
end

function Base.filter!(f, seq::LongSequence)
    writeindex = bitindex(seq, 0)
    chunk = zero(UInt64)
    bps = bits_per_symbol(seq)
    for i in eachindex(seq)
        encoded_symbol = extract_encoded_element(bitindex(seq, i), encoded_data(seq))
        symbol = decode(Alphabet(seq), encoded_symbol)
        if f(symbol)
            writeindex += bps
            chunk |= encoded_symbol << offset(writeindex)
            # If the offset is now zero, the chunk is full
            if iszero(offset(writeindex + bps))
                seq.data[index(writeindex)] = chunk
                chunk = zero(UInt64)
            end
        end
    end
    # If it's not zero, we need to write the first chunk
    if !iszero(offset(writeindex + bps))
        seq.data[index(writeindex)] = chunk
    end
    resize!(seq, div((writeindex + bps).val, bps))
end

function Base.map!(f, seq::LongSequence)
    for i in 1:lastindex(seq)
        unsafe_setindex!(seq, f(inbounds_getindex(seq, i)), i)
    end
    return seq
end

"""
    reverse!(seq::LongSequence)

Reverse a biological sequence `seq` in place.
"""
Base.reverse!(seq::LongSequence{<:Alphabet}) = _reverse!(seq, BitsPerSymbol(seq))

"""
    reverse(seq::LongSequence)

Create reversed copy of a biological sequence.
"""
Base.reverse(seq::LongSequence{<:Alphabet}) = _reverse(seq, BitsPerSymbol(seq))

# Fast path for non-inplace reversion
@inline function _reverse(seq::LongSequence{A}, B::BT) where {A <: Alphabet,
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    cp = LongSequence{A}(unsigned(length(seq)))
    reverse_data_copy!(identity, cp.data, seq.data, seq_data_len(seq) % UInt, B)
    return zero_offset!(cp)
end

_reverse(seq::LongSequence{<:Alphabet}, ::BitsPerSymbol) = reverse!(copy(seq))

# Generic fallback
function _reverse!(seq::LongSequence{<:Alphabet}, ::BitsPerSymbol)
    i, j = 1, lastindex(seq)
    @inbounds while i < j
        seq[i], seq[j] = seq[j], seq[i]
        i += 1
        j -= 1
    end
    return seq
end

@inline function _reverse!(seq::LongSequence{<:Alphabet}, B::BT) where {
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    # We need to account for the fact that the seq may not use all its stored data
    reverse_data!(identity, seq.data, seq_data_len(seq) % UInt, B)
    return zero_offset!(seq)
end

# Reversion of chunk bits may have left-shifted data in chunks, this function right shifts
# all chunks by up to 63 bits.
# This is written so it SIMD parallelizes - careful with changes
@inline function zero_offset!(seq::LongSequence{A}) where A <: Alphabet
    offs = (64 - offset(bitindex(seq, length(seq)) + bits_per_symbol(A()))) % UInt
    zero_offset!(seq, offs)
end

@inline function zero_offset!(seq::LongSequence{A}, offs::UInt) where A <: Alphabet
    rshift = offs
    lshift = 64 - rshift
    len = length(seq.data)
    @inbounds if !iszero(lshift)
        this = seq.data[1]
        for i in 1:len-1
            next = seq.data[i+1]
            seq.data[i] = (this >>> (unsigned(rshift) & 63)) | (next << (unsigned(lshift) & 63))
            this = next
        end
        seq.data[len] >>>= (unsigned(rshift) & 63)
    end
    return seq
end

# Reverse chunks in data vector and each symbol within a chunk. Chunks may have nonzero
# offset after use, so use zero_offset!
@inline function reverse_data!(pred, data::Vector{UInt64}, len::UInt, B::BT) where {
    BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    @inbounds @simd ivdep for i in 1:len >>> 1
        data[i], data[len-i+1] = pred(reversebits(data[len-i+1], B)), pred(reversebits(data[i], B))
    end
    @inbounds if isodd(len)
        data[len >>> 1 + 1] = pred(reversebits(data[len >>> 1 + 1], B))
    end
end

@inline function reverse_data_copy!(pred, dst::Vector{UInt64}, src::Vector{UInt64}, len::UInt,
    B::BT) where {BT <: Union{BitsPerSymbol{2}, BitsPerSymbol{4}, BitsPerSymbol{8}}}
    @inbounds @simd for i in eachindex(dst)
        dst[i] = pred(reversebits(src[len - i + 1], B))
    end
end

"""
    complement!(seq)

Make a complement sequence of `seq` in place.
"""
function complement!(seq::LongSequence{A}) where {A<:NucleicAcidAlphabet}
    seqdata = seq.data
    @inbounds for i in eachindex(seqdata)
        seqdata[i] = complement_bitpar(seqdata[i], Alphabet(seq))
    end
    return seq
end

function complement!(s::LongSubSeq{A}) where {A <: NucleicAcidAlphabet}
    bps = bits_per_symbol(A())
    bi = firstbitindex(s)
    i = 1
    stop = lastbitindex(s) + bps
    @inbounds while (!iszero(offset(bi)) & (bi < stop))
        s[i] = complement(s[i])
        bi += bps
        i += 1
    end
    @inbounds for j in index(bi):index(stop)-1
        s.data[j] = complement_bitpar(s.data[j], Alphabet(s))
        bi += 64
        i += symbols_per_data_element(s)
    end
    @inbounds while bi < stop
        s[i] = complement(s[i])
        bi += bps
        i += 1
    end
    return s
end

function reverse_complement!(seq::LongSequence{<:NucleicAcidAlphabet})
    pred = x -> complement_bitpar(x, Alphabet(seq))
    reverse_data!(pred, seq.data, seq_data_len(seq) % UInt, BitsPerSymbol(seq))
    return zero_offset!(seq)
end

function reverse_complement(seq::LongSequence{<:NucleicAcidAlphabet})
    cp = typeof(seq)(unsigned(length(seq)))
    pred = x -> complement_bitpar(x, Alphabet(seq))
    reverse_data_copy!(pred, cp.data, seq.data, seq_data_len(seq) % UInt, BitsPerSymbol(seq))
    return zero_offset!(cp)
end

function Random.shuffle!(seq::LongSequence)
    # Fisher-Yates shuffle
    @inbounds for i in 1:lastindex(seq) - 1
        j = rand(i:lastindex(seq))
        seq[i], seq[j] = seq[j], seq[i]
    end
    return seq
end


text:
#-------------------------------------------------------------------------------
# File: functions.jl
# Description: This file contains functions that are used in ILS and VNS.
# Date: October 31, 2019
# Author: Jonathan Fontaine, Killian Fretaud, Rémi Garcia,
#         Boualem Lamraoui, Benoît Le Badezet, Benoit Loger
#-------------------------------------------------------------------------------

# More readable
include("move_insertion.jl")
include("move_exchange.jl")

"""
    sequence_insert!(sequence::Array{Int, 1}, index_remove::Int, index_insert::Int)

Delete the item at `index_remove` and insert at `index_insert` in the sequence.
"""
function sequence_insert!(sequence::Array{Int, 1}, index_remove::Int, index_insert::Int)
    temp = sequence[index_remove]
    if index_remove < index_insert
        for k in index_remove:(index_insert-1)
            sequence[k] = sequence[k+1]
        end
    else
        for k in index_remove:-1:(index_insert+1)
            sequence[k] = sequence[k-1]
        end
    end
    sequence[index_insert] = temp
end

"""

"""
function penalize_costs!(costs::Array{Int, 2}, init_position::Int,
                         solution::Solution, instance::Instance)
    sequence = copy(solution.sequence)
    b0 = instance.nb_late_prec_day+1
    sequence_insert!(sequence, init_position, b0)
    if !is_sequence_valid(sequence, solution.length, instance)
        costs[b0, :] .= instance.nb_cars
    end

    for position in b0:(solution.length-1)
        sequence[position], sequence[position+1] = sequence[position+1], sequence[position]
        if !is_sequence_valid(sequence, solution.length, instance)
            costs[position+1, :] .= instance.nb_cars
        end
    end

    return costs
end

"""
    greedy_add!(solution::Solution, instance::Instance, k::Int, objectives::Int)

Reinserts last car in the sequence of `solution`.
"""
function greedy_add!(solution::Solution, instance::Instance, position_car::Int,
                     objectives::Int, valid_sequence::Bool = false)
    @assert objectives <= 3
    @assert objectives >= 1
    b0 = instance.nb_late_prec_day + 1
    costs = cost_move_insertion(solution, position_car, instance, objectives)
    if valid_sequence
        penalize_costs!(costs, position_car, solution, instance)
    end
    delta = weighted_sum(costs[b0, :])
    best_delta = delta
    best_position = b0
    for position in (b0+1):solution.length
        delta = weighted_sum(costs[position, :])
        if delta < best_delta
            best_delta = delta
            best_position = position
        end
    end
    #Limit case:
    # No position allows a valid sequence
    if valid_sequence && best_delta == instance.nb_cars
        return greedy_add!(solution, instance, position_car-1, objectives, valid_sequence)
    end
    solution.length += 1
    move_insertion!(solution, position_car, best_position, instance)

    return solution
end

"""
    remove!(solution_init::Solution, instance::Instance,
            k::Int, crit::Array{Int, 1})

Removes `k` critical cars of the sequence of `solution_init`.
"""
function remove!(solution::Solution, instance::Instance,
                 k::Int, crit::Array{Int, 1})
    indices = sort(randperm(length(crit))[1:k])
    crit_sort = sort(crit, rev = true)
    for i in 1:k
        position = crit_sort[indices[i]]
        solution.length -= 1
        move_insertion!(solution, position, instance.nb_cars, instance)
    end

    return solution
end

"""
    find_critical_cars(solution::Solution, instance::Instance,
                       objectives::BitArray{1})

Returns the set cars involved in at least one HPRC violation and/or one LPRC
violation.
"""
function find_critical_cars(solution::Solution, instance::Instance,
                            objectives::BitArray{1})
    @assert length(objectives) == 2
    @assert true in objectives
    if objectives[1]
        first_option = 1
    else
        first_option = instance.nb_HPRC+1
    end
    if objectives[2]
        last_option = instance.nb_HPRC+instance.nb_LPRC
    else
        last_option = instance.nb_HPRC
    end

    critical_cars = Set{Int}()
    b0 = instance.nb_late_prec_day + 1
    for index_car in b0:solution.length
        for option in first_option:last_option
            if solution.M1[option, index_car] > instance.RC_p[option]
                index_car_lim = index_car + min(instance.RC_p[option], solution.length-index_car)
                for index_car_add in index_car:index_car_lim
                    if instance.RC_flag[solution.sequence[index_car_add], option]
                        push!(critical_cars, index_car_add)
                    end
                end
            end
        end
    end
    return collect(critical_cars)
end

find_critical_cars(solution::Solution, instance::Instance, objectives::Int) =
    find_critical_cars(solution::Solution, instance::Instance, [trues(objectives) ; falses(2-objectives)])

"""
    cost(solution::Solution, instance::Instance, objectives::BitArray{1})

Returns the (partial) set of objective values (without weights) of `solution`.
"""
function cost(solution::Solution, instance::Instance, objectives::BitArray{1})
    @assert length(objectives) == 3
    @assert true in objectives
    @assert solution.length == instance.nb_cars
    cost_on_objective = zeros(Int, 3)

    if objectives[1]
        for car in 1:solution.length
            for option in 1:instance.nb_HPRC
                cost_on_objective[1] += max(0 , solution.M1[option, car] - instance.RC_p[option])
            end
        end
    end

    if objectives[2]
        for car in 1:solution.length
            for option in (instance.nb_HPRC+1):(instance.nb_HPRC+instance.nb_LPRC)
                cost_on_objective[2] += max(0 , solution.M1[option, car] - instance.RC_p[option])
            end
        end
    end

    if objectives[3]
        if solution.colors === nothing
            for i in (instance.nb_late_prec_day+2):solution.length
                if instance.color_code[solution.sequence[i]] != instance.color_code[solution.sequence[i-1]]
                    cost_on_objective[3] += 1
                end
            end
        else
            cost_on_objective[3] = -1
            index_batch = instance.nb_late_prec_day + 1
            while index_batch <= solution.length
                index_batch += solution.colors[index_batch].width
                cost_on_objective[3] += 1
            end
        end
    end

    return cost_on_objective
end

"""
    cost(solution::Solution, instance::Instance, objectives::Int)

Returns the set of unweighted values for objectives 1 to `objectives` of `solution`.
"""
cost(solution::Solution, instance::Instance, objectives::Int) =
    cost(solution::Solution, instance::Instance, [trues(objectives) ; falses(3-objectives)])

cost(solution::Solution, instance::Instance) =
    cost(solution::Solution, instance::Instance, 3)

"""
    weighted_sum(cost_solution::Array{Int, 1})

Returns the partial weighted sum of the objective values.
"""
function weighted_sum(cost_solution::Array{Int, 1})
    return sum(cost_solution[i] * WEIGHTS_OBJECTIVE_FUNCTION[i] for i in 1:length(cost_solution))
end

weighted_sum(solution::Solution, instance::Instance) =
    weighted_sum(cost(solution, instance))

"""
    print_cost(solution::Solution, instance::Instance)

Prints the objective values of the solution.
"""
function print_cost(solution::Solution, instance::Instance)
    cost_solution = cost(solution, instance, 3)
    println("\tHPRC violations: ", cost_solution[1])
    println("\tLPRC violations: ", cost_solution[2])
    println("\tPCC  violations: ", cost_solution[3])
    println("\tObjective value is ", weighted_sum(cost_solution))
end

"""

"""
#TODO: why n, why sequence, why not solution ?
function is_sequence_valid(sequence::Array{Int, 1}, n::Int, instance::Instance)
    counter = 1
    for car_pos in (instance.nb_late_prec_day+2):n
        if instance.color_code[sequence[car_pos-1]] == instance.color_code[sequence[car_pos]]
            counter += 1
        else
            counter = 1
        end
        if counter > instance.nb_paint_limitation
            return false
        end
    end
    return true
end


text:
###########################################################
# Part 1: Helper functions

"""
    precompilefuncs(resid, RJ, ::Val{N}) where N

Add code that precompiles the given `resid` and `RJ` functions together
with the dual-number arythmetic required by ForwardDiff.

!!! warning
    Internal function. Do not call directly

# Implementation (for developers)
"""
function precompilefuncs(resid, RJ, ::Val{N}) where N
    ccall(:jl_generating_output, Cint, ()) == 1 || return nothing

    tag = ForwardDiff.Tag{resid,Float64}
    dual = ForwardDiff.Dual{tag,Float64,N}
    duals = Array{dual,1}
    cfg = ForwardDiff.GradientConfig{tag,Float64,N,duals}
    mdr = DiffResults.MutableDiffResult{1,Float64,Tuple{Array{Float64,1}}}

    precompile(resid, (Array{Float64,1},)) || error("precompile")
    precompile(resid, (duals,)) || error("precompile")
    precompile(RJ, (Array{Float64,1},)) || error("precompile")

    for pred in (ForwardDiff.UNARY_PREDICATES ∪ Symbol[:-, :+, :log, :exp])
        pred ∈ (:iseven, :isodd) || precompile(getfield(Base, pred), (Float64,)) || error("precompile")
        precompile(getfield(Base, pred), (dual,)) || error("precompile")
    end

    for pred in ForwardDiff.BINARY_PREDICATES ∪ Symbol[:+, :-, :*, :/, :^]
        precompile(getfield(Base, pred), (Float64, Float64)) || error("precompile")
        precompile(getfield(Base, pred), (dual, Float64)) || error("precompile")
        precompile(getfield(Base, pred), (Float64, dual)) || error("precompile")
        precompile(getfield(Base, pred), (dual, dual)) || error("precompile")
    end

    precompile(ForwardDiff.extract_gradient!, (Type{tag}, mdr, dual)) || error("precompile")
    precompile(ForwardDiff.vector_mode_gradient!, (mdr, typeof(resid), Array{Float64,1}, cfg)) || error("precompile")

    # precompile(Tuple{typeof(ForwardDiff.extract_gradient!), Type{tag}, mdr, dual}) || error("precompile")
    # precompile(Tuple{typeof(ForwardDiff.vector_mode_gradient!), mdr, resid, Array{Float64, 1}, cfg}) || error("precompile")

    return nothing
end

"""
    funcsyms(mod::Module)

Create a pair of identifiers that does not conflict with existing identifiers
in the given module.

!!! warning
    Internal function. Do not call directly.

### Implementation (for developers)
We need two identifiers `resid_N` and `RJ_N` where "N" is some integer number.
The first is going to be the name of the function that evaluates the equation and
the second is going to be the name of the function that evaluates both the equation
and its gradient.
"""
function funcsyms end

let funcsyms_state = 0
    global funcsyms_counter() = (funcsyms_state += 1)
end
function funcsyms(mod::Module)
    fn1, fn2 = mod.eval(quote
        let nms = names(@__MODULE__; all = true)
            num = $(@__MODULE__).funcsyms_counter()
            local fn1 = Symbol("resid_$num")
            local fn2 = Symbol("RJ_$num")
            while fn1 ∈ nms || fn2 ∈ nms
                num = $(@__MODULE__).funcsyms_counter()
                fn1 = Symbol("resid_$num")
                fn2 = Symbol("RJ_$num")
            end
            fn1, fn2
        end
    end)
end

"""
    makefuncs(expr, vsyms [, params_expr]; mod::Module)

Create two functions that evaluate the residual and its gradient for the given
expression.

!!! warning
    Internal function. Do not call directly.

### Arguments
- `expr`: the expression
- `vsyms`: a list of variables in the expression.
- `params_expr`: an expression to be included in the function code before evaluating
    the expression. This allows for assigning values of parameters in the expression,
    i.e. variables that are not arguments to the funcion.

### Return value
Return a quote block to be evaluated in the module where the model is being defined.
The quote block contains definitions of the residual function and a second function
that evaluates both the residual and its gradient.
"""
function makefuncs(expr, vsyms, params_expr = nothing; mod::Module)
    fn1, fn2 = funcsyms(mod)
    x = gensym("x")
    nargs = length(vsyms)
    return quote
        function $fn1($x::AbstractVector{T}) where T <: Real
            ($(vsyms...),) = $x
            $(params_expr)
            $expr
        end
        const $fn2 = EquationGradient($fn1, Val{$nargs}())
        $(@__MODULE__).precompilefuncs($fn1, $fn2, Val{$nargs}())
        ($fn1, $fn2)
    end
end

"""
    initfuncs(mod::Module)

Initialize the given module before creating functions that evaluate residuals
and thier gradients.

"""
function initfuncs(mod::Module)
    if :MyTag ∉ names(mod; all = true)
        mod.eval(quote
            struct MyTag end
            struct EquationGradient{DR,CFG} <: Function
                fn1::Function
                dr::DR
                cfg::CFG
            end
            EquationGradient(fn1::Function, ::Val{N}) where N = EquationGradient(fn1,
                $(@__MODULE__).DiffResults.DiffResult(zero(Float64), zeros(Float64, N)),
                $(@__MODULE__).ForwardDiff.GradientConfig(fn1, zeros(Float64, N), $(@__MODULE__).ForwardDiff.Chunk{N}(), MyTag))
            function (s::EquationGradient)(x::AbstractVector{Float64})
                $(@__MODULE__).ForwardDiff.gradient!(s.dr, s.fn1, x, s.cfg)
                return s.dr.value, s.dr.derivs[1]
            end
        end)
    end
    return nothing
end

###########################################################
# Part 2: Evaluation data for models and equations

"""
    AbstractModelEvaluationData

Base type for all model evaluation structures.
Specific derived types would specialize in different types of models.

### Implementaion (for developers)
Derived types must specialize two functions
  * [`eval_R!`](@ref) - evaluate the residual
  * [`eval_RJ`](@ref) - evaluate the residual and its Jacobian
"""
abstract type AbstractModelEvaluationData end

"""
    eval_R!(res::AbstractArray{Float64,1}, point::AbstractArray{Float64, 2}, ::MED) where MED <: AbstractModelEvaluationData

Evaluate the model residual at the given point using the given model evaluation structure.
The residual is stored in the provided vector.

### Implementation details (for developers)
When creating a new type of model evaluation data, you must define a
method of this function specialized to it.

The `point` argument will be a 2d array, with the number of rows equal to
`maxlag+maxlead+1` and the number of columns equal to the number of `variables+shocks+auxvars` of the model.
The `res` vector will have the same length as the number of equations + auxiliary equations.
Your implementation must not modify `point` and must update `res`.

See also: [`eval_RJ`](@ref)
"""
function eval_R! end
export eval_R!
eval_R!(res::AbstractVector{Float64}, point::AbstractMatrix{Float64}, ::AMED) where AMED <: AbstractModelEvaluationData = throw(NotImplementedError(AMED))

"""
    eval_RJ(point::AbstractArray{Float64, 2}, ::MED) where MED <: AbstractModelEvaluationData

Evaluate the model residual and its Jacobian at the given point using the given model evaluation structure.
Return a tuple, with the first element being the residual and the second element being the Jacobian.

### Implementation details (for developers)
When creating a new type of model evaluation data, you must define a
method of this function specialized to it.

The `point` argument will be a 2d array, with the number of rows equal to
`maxlag+maxlead+1` and the number of columns equal to the number of `variables+shocks+auxvars` of the model.
Your implementation must not modify `point` and must return the tuple of (residual, Jacobian) evaluated
at the given `point`. The Jacobian is expected to be `SparseMatrixCSC` (*this might change in the future*).

See also: [`eval_R!`](@ref)
"""
function eval_RJ end
export eval_RJ
eval_RJ(point::AbstractMatrix{Float64}, ::AMED) where AMED <: AbstractModelEvaluationData = throw(NotImplementedError(AMED))

# Model Evaluation Data that doesn't exist.

"""
    struct NoModelEvaluationData <: AbstractModelEvaluationData

Specific type that indicates that the model cannot be evaluated.
This is used as a placeholder while the model is being defined.
During initialization, the actual model evaluation data is created.
"""
struct NoModelEvaluationData <: AbstractModelEvaluationData end
global const NoMED = NoModelEvaluationData()
eval_R!(res::AbstractVector{Float64}, point::AbstractMatrix{Float64}, ::NoModelEvaluationData) = throw(ModelNotInitError())
eval_RJ(point::AbstractMatrix{Float64}, ::NoModelEvaluationData) = throw(ModelNotInitError())

# The standard Model Evaluation Data used in the general case.

"""
    ModelEvaluationData <: AbstractModelEvaluationData

The standard model evaluation data used in the general case and by default.
"""
struct ModelEvaluationData{E <: AbstractEquation,I} <: AbstractModelEvaluationData
    alleqns::Vector{E}
    allinds::Vector{I}
    "Placeholder for the Jacobian matrix"
    J::SparseMatrixCSC{Float64,Int64}
    "Placeholder for the residual vector"
    R::Vector{Float64}
    rowinds::Vector{Vector{Int64}}
end

"""
    ModelEvaluationData(model::AbstractModel)

Create the standard evaluation data structure for the given model.
"""
function ModelEvaluationData(model::AbstractModel)
    time0 = 1 + model.maxlag
    alleqns = model.alleqns
    neqns = length(alleqns)
    allinds = @timer [[CartesianIndex((time0 + ti, vi)) for (ti, vi) in eqn.vinds] for eqn in alleqns]
    ntimes = 1 + model.maxlag + model.maxlead
    nvars = length(model.allvars)
    LI = LinearIndices((ntimes, nvars))
    II = @timer reduce(vcat, (fill(Int64(i), length(eqn.vinds)) for (i, eqn) in enumerate(alleqns)))
    JJ = @timer [LI[inds] for inds in allinds]
    M = @timer SparseArrays.sparse(II, reduce(vcat, JJ), similar(II), neqns, ntimes * nvars)
    M.nzval .= @timer 1:length(II)
    rowinds = @timer [copy(M[i,LI[inds]].nzval) for (i, inds) in enumerate(JJ)]
    ModelEvaluationData(alleqns, allinds, similar(M, Float64), Vector{Float64}(undef, neqns), rowinds)
end

function eval_R!(res::AbstractVector{Float64}, point::AbstractMatrix{Float64}, med::ModelEvaluationData)
    # med === NoMED && throw(ModelNotInitError())
    for (i, eqn, inds) in zip(1:length(med.alleqns), med.alleqns, med.allinds)
        res[i] = eqn.eval_resid(point[inds])
    end
    return nothing
end

function eval_RJ(point::AbstractMatrix{Float64}, med::ModelEvaluationData)
    # med === NoMED && throw(ModelNotInitError())
    neqns = length(med.alleqns)
    res = med.R
    jac = med.J
    for (i, eqn, inds, ri) in zip(1:neqns, med.alleqns, med.allinds, med.rowinds)
        res[i], jac.nzval[ri] = eqn.eval_RJ(point[inds])
    end
    return res, jac
end


text:
using NamedDims

# Abbreviations for things which have both names & keys:

NdaKa{L,T,N} = NamedDimsArray{L,T,N,<:KeyedArray{T,N}}
KaNda{L,T,N} = KeyedArray{T,N,<:NamedDimsArray{L,T,N}}

NdaKaVoM{L,T} = Union{NamedDimsArray{L,T,1,<:KeyedArray}, NamedDimsArray{L,T,2,<:KeyedArray}}
# NdaKaV{L,T} = NamedDimsArray{L,T,1,<:KeyedArray{T,1}}

# NamedDims functionality:

NamedDims.dimnames(A::KaNda{L}) where {L} = L
NamedDims.dimnames(A::KaNda{L,T,N}, d::Int) where {L,T,N} = d <= N ? L[d] : :_

NamedDims.dim(A::KaNda{L}, name) where {L} = NamedDims.dim(L, name)

Base.axes(A::KaNda{L}, s::Symbol) where {L} = axes(A, NamedDims.dim(L,s))
Base.size(A::KaNda{L,T,N}, s::Symbol) where {T,N,L} = size(A, NamedDims.dim(L,s))

NamedDims.rename(A::KaNda, names...) = KeyedArray(rename(parent(A), names...), axiskeys(A))

# Extra complication to make wrappers commutative:

hasnames(A::KaNda) = true
hasnames(A::NamedDimsArray) = true
hasnames(A) = false

NamedDims.unname(A::KaNda) = KeyedArray(unname(parent(A)), axiskeys(A))
keyless(A::NdaKa{L}) where {L} = NamedDimsArray(parent(parent(A)), L)

axiskeys(A::NdaKa) = axiskeys(parent(A))
axiskeys(A::NdaKa, d::Int) = axiskeys(parent(A), d)

axiskeys(A::NdaKa{L}, s::Symbol) where {L} = axiskeys(parent(A), NamedDims.dim(L,s))
axiskeys(A::KaNda{L}, s::Symbol) where {L} = axiskeys(A, NamedDims.dim(L,s))

haskeys(A::NdaKa) = true
haskeys(A::KeyedArray) = true
haskeys(A) = false

keys_or_axes(A) = haskeys(A) ? axiskeys(A) : axes(A)
keys_or_axes(A, d) = haskeys(A) ? axiskeys(A, d) : axes(A, d)

# Double un-wrappers:

keyless_unname(A::NdaKa) = parent(parent(A))
keyless_unname(A::KaNda) = parent(parent(A))
keyless_unname(A::NamedDimsArray) = parent(A)
keyless_unname(A::KeyedArray) = parent(A)
keyless_unname(A) = A

# Re-constructors:

function KeyedArray(A::NdaKa, r2::Tuple)
    r3 = unify_keys(axiskeys(parent(A)), r2)
    KeyedArray(keyless(A), r3)
end

function NamedDims.NamedDimsArray(A::KaNda{L}, L2::Tuple) where {L}
    L3 = NamedDims.unify_names(L, L2)
    NamedDimsArray(NamedDims.unname(A), L3)
end

# getproperty: it's useful to say for `(i,t) in enumerate(A.time)` etc.
# This will make saying ".data" slow (by 30ns), fixed in NamedDims.jl#78

Base.propertynames(A::NdaKa{L}, private::Bool=false) where {L} =
    private ? (L..., fieldnames(typeof(A))...) : L
Base.propertynames(A::KaNda{L}, private::Bool=false) where {L} =
    private ? (L..., fieldnames(typeof(A))...) : L

Base.getproperty(A::NdaKa{L}, s::Symbol) where {L} =
    Base.sym_in(s, L) ? axiskeys(A, NamedDims.dim(L, s)) : getfield(A, s)
Base.getproperty(A::KaNda{L}, s::Symbol) where {L} =
    Base.sym_in(s, L) ? axiskeys(A, NamedDims.dim(L, s)) : getfield(A, s)
Base.getproperty(A::NamedDimsArray{L}, s::Symbol) where {L} =
    Base.sym_in(s, L) ? axes(A, NamedDims.dim(L, s)) : getfield(A, s) # 🏴☠?

# Keyword indexing of KeyedArray:

@inline @propagate_inbounds function Base.getindex(A::KeyedArray; kw...)
    hasnames(A) || error("must have names!")
    inds = NamedDims.order_named_inds(Val(dimnames(A)); kw...)
    getindex(A, inds...)
end
@inline @propagate_inbounds function Base.view(A::KeyedArray; kw...)
    hasnames(A) || error("must have names!")
    inds = NamedDims.order_named_inds(Val(dimnames(A)); kw...)
    view(A, inds...)
end

# Any NamedDimsArray + KeyedArray combination is callable:

@inline @propagate_inbounds (A::NdaKa)(args...) = getkey(A, args...)

@inline @propagate_inbounds (A::KaNda)(;kw...) = getkey(A; kw...)
@inline @propagate_inbounds (A::NdaKa)(;kw...) = getkey(A; kw...)

@inline @propagate_inbounds function getkey(A; kw...)
    list = dimnames(A)
    issubset(kw.itr, list) || error("some keywords not in list of names!")
    args = map(s -> Base.sym_in(s, kw.itr) ? getfield(kw.data, s) : Colon(), list)
    A(args...)
end

# Constructors, including pirate method (A; kw...)

_construc_doc = """
    KeyedArray(A; i=2:3, j=["a", "b"])
    NamedDimsArray(A; i=2:3, j=["a", "b"])

These constructors make `KeyedArray(NamedDimsArray(A, names), keys)`
or `NamedDimsArray(KeyedArray(A, keys), names)`, which should be equivalent.

These perform less sanity checking than `wrapdims(A; kw...)`.
"""
@doc _construc_doc
function KeyedArray(A::AbstractArray; kw...)
    L = keys(values(kw))
    length(L) == ndims(A) || throw(ArgumentError("number of names must match number of dimensions"))
    R = values(values(kw))
    map(x -> axes(x, 1), R) == axes(A) || throw(ArgumentError("axes of keys must match axes of array"))
    KeyedArray(NamedDimsArray(A, L), R)
end

@doc _construc_doc
function NamedDims.NamedDimsArray(A::AbstractArray; kw...)
    L = keys(values(kw))
    length(L) == ndims(A) || throw(ArgumentError("number of names must match number of dimensions"))
    R = values(values(kw))
    map(x -> axes(x, 1), R) == axes(A) || throw(ArgumentError("axes of keys must match axes of array"))
    NamedDimsArray(KeyedArray(A, R), L)
end

"""
    named_axiskeys(arr)::NamedTuple

Return the [`axiskeys`](@ref) along with their names.
If there are duplicate names or unnamed axes, an error is thrown.

```jldoctest
julia> using AxisKeys

julia> arr = KeyedArray(rand(1,3), x=[1], y=[2,3,4]);

julia> named_axiskeys(arr)
(x = [1], y = [2, 3, 4])
```
"""
function named_axiskeys(A::AbstractArray)
    NT = NamedTuple{dimnames(A)}
    NT(axiskeys(A))
end


text:
###############################################################################
#
# low-level in-place arithmetic operators
# should probably be added into Nemo upstream
#
#=
function set!(x::fmpz, y::Int)
    ccall((:fmpz_set_si, Nemo.libflint), Nothing, (Ref{fmpz}, Int), x, y)
    return x
end

function set!(x::fmpq, y::fmpq)
    ccall((:fmpz_set, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}), x, y)
    return x
end

function set!(x::fmpq, n::Int, d::Int)
        ccall((:fmpq_set_si, Nemo.libflint), Nothing, (Ref{fmpq}, Int, UInt), x, n, UInt(d))
    return x
end

function set!(x::fmpq, n::fmpz, d::fmpz)
    ccall((:fmpq_set_fmpz_frac, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpz}, Ref{fmpz}), x, n, d)
    return x
end

function Nemo.mul!(z::fmpq, x::fmpq, y::fmpz)
    ccall((:fmpq_mul_fmpz, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), z, x, y)
    return z
end

function Nemo.mul!(z::fmpq, x::fmpq, y::Int)
    ccall((:fmpq_mul_si, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Int), z, x, y)
    return z
end

function Nemo.sub!(z::fmpq, x::fmpq, y::Int)
    ccall((:fmpq_sub_si, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Int), z, x, y)
    return z
end

function Nemo.inv!(x::fmpq, y::fmpq)
    ccall((:fmpq_inv, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}), x, y)
    return x
end

function pow!(c::fmpq, a::fmpq, b::Int)
    iszero(a) && b < 0 && throw(DivideError())
    ccall((:fmpq_pow_si, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Int), c, a, b)
    return c
end

function div!(z::fmpq, x::fmpq, y::fmpq)
    iszero(y) && throw(DivideError())
    ccall((:fmpq_div, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), z, x, y)
    return z
end

muleq!(x::fmpq, y::fmpq) = mul!(x, x, y)
muleq!(x::fmpq, y::fmpz) = mul!(x, x, y)
muleq!(x::fmpq, y::Int) = mul!(x, x, y)
subeq!(x::fmpq, y::Int) = Nemo.sub!(x, x, y)
poweq!(x::fmpq, y::Int) = pow!(x, x, y)
diveq!(x::fmpq, y::fmpq) = div!(x, x, y)
Nemo.inv!(x::fmpq) = Nemo.inv!(x, x)
=#
"""
    pow_eq!(a, b)
Equivalent to a ^= b
"""
function pow_eq!(a::fmpq, b::Int)::Nothing
    #iszero(a) && b < 0 && throw(DivideError())
    ccall((:fmpq_pow_si, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Int), a, a, b)
end

"""
    mul_eq!(a, b)
Equivalent to a *= b
"""
function mul_eq!(a::fmpq, b::fmpq)::Nothing
    ccall((:fmpq_mul, Nemo.libflint), Nothing,(Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), a, a, b)
end
function mul_eq!(a::fmpq, b::Int)::Nothing
    ccall((:fmpq_mul_si, Nemo.libflint), Nothing,(Ref{fmpq}, Ref{fmpq}, Int), a, a, b)
end

"""
    add_eq!(a, b)
Equivalent to a += b
"""
function add_eq!(a::fmpq, b::fmpq)::Nothing
    ccall((:fmpq_add, Nemo.libflint), Nothing,(Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), a, a, b)
end
function add_eq!(a::Vector{fmpq}, b::Vector{fmpq})::Nothing
    for i in 1:length(a)
        ccall((:fmpq_add, Nemo.libflint), Nothing,(Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), a[i], a[i], b[i])
    end
end

"""
    div_eq!(a, b)
Equivalent to a /= b
"""
function div_eq!(a::fmpq, b::fmpq)::Nothing
    ccall((:fmpq_div, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), a, a, b)
end
function div_eq!(a::fmpq, b::Int)::Nothing
    ccall((:fmpq_div_fmpz, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Ref{fmpz}), a, a, fmpz(b))
end

"""
    eq!(a, b)
Equivalent to a = b
"""
function eq!(a::fmpq, b::fmpq)::Nothing
    ccall((:fmpq_set, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}), a, b)
end

"""
    sub!(a, b)
Equivalent to a -= b
"""
function sub!(a::fmpq, b::fmpq)::Nothing
    ccall((:fmpq_sub, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}, Ref{fmpq}), a, a, b)
end

"""
    neg!(a, b)
Equivalent to a *= -1
"""
function neg!(a::fmpq)::Nothing
    ccall((:fmpq_neg, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}), a, a)
end

"""
    inv!(a, b)
Equivalent to a /= 1
"""
function inv!(a::fmpq)::Nothing
    ccall((:fmpq_inv, Nemo.libflint), Nothing, (Ref{fmpq}, Ref{fmpq}), a, a)
end


text:
type G2SwaptionEngine{AffineModelType, Y <: YieldTermStructure} <: PricingEngine
  model::G2{AffineModelType, Y}
  range::Float64
  intervals::Int

  call{Y}(::Type{G2SwaptionEngine}, model::G2{AffineModelType, Y}, range::Float64, intervals::Int) = new{AffineModelType, Y}(model, range, intervals)
end

# methods #
function _calculate!(pe::G2SwaptionEngine, swaption::Swaption)
  swap = swaption.swap

  # overriding pricing engine
  _calculate!(DiscountingSwapEngine(pe.model.ts), swap)
  swap.lazyMixin.calculated = true

  correction = swap.spread * abs(floating_leg_BPS(swap) / fixed_leg_BPS(swap))
  fixedRate = swap.fixedRate - correction
  swaption.results.value = gen_swaption(pe.model, swaption, fixedRate, pe.range, pe.intervals)

  return swaption
end


text:
using GivensBoson,LinearAlgebra
begin
	export adj_mat,boundary_condition

    abstract type BoundaryType end
    abstract type Free <: BoundaryType end
    abstract type Periodic <: BoundaryType end
    abstract type Cylinder <: BoundaryType end
    struct BoundaryCondition{T<:BoundaryType} name::String end
    boundary_condition(name::String) = BoundaryCondition{eval(Meta.parse(name))}(name)

    """
    adj_mat is a adjacent matrix builder which will return a matrix which is simple to use in site-site coupling system.
    The result is a matrix A_{ij}
    A_{ij}=k means the j-th bond of site-i is point to site-k
    """
    function adj_mat(Lx::Int,Ly::Int,bc::BoundaryCondition{Periodic};D=2)
        N = Lx*Ly
        adj = zeros(Int,N,4) # up right down left
        n(i::Int,j::Int) = (i-1)*Ly+j
        function periodic_index(x::Int,Lx::Int)
            if x == 0
                return Lx
            elseif x < Lx+1
                return x
            else
                return x-Lx
            end
        end
        pix = periodic_index

        #
        for i = 1:Lx # row
            for j = 1:Ly # column
                id = n(i,j)
                adj[id,1] = n(i,pix(j-1,Ly))
                adj[id,2] = n(pix(i+1,Lx),j)
                adj[id,3] = n(i,pix(j+1,Ly))
                adj[id,4] = n(pix(i-1,Lx),j)
            end
        end

        if D==1
            adj[:,2].=0
            adj[:,4].=0
        end
        return adj
    end

    function adj_mat(Lx::Int,Ly::Int,bc::BoundaryCondition{Free};D=2)
        adj = adj_mat(Lx,Ly,boundary_condition("Periodic"))
        n(i::Int,j::Int) = (i-1)*Ly+j
        for i = 1:Lx
            for j =1:Ly
                id = n(i,j)
                if j == 1
                    adj[id,1] = 0
                end
                if j == Ly
                    adj[id,3] =0
                end

                if i == 1
                    adj[id,4] = 0
                end
                if i == Lx
                    adj[id,2] =0
                end
            end
        end
        return adj
    end

    """ Lx = 5 ; Ly = 2
    |-1 2 3 4 5-|
    |-6 7 8 9 10-|
    -------
    """
    function adj_mat(Lx::Int,Ly::Int,bc::BoundaryCondition{Cylinder};D=2)
        adj = adj_mat(Lx,Ly,boundary_condition("Periodic"))
        n(i::Int,j::Int) = (i-1)*Ly+j
        for i = 1:Lx
            for j =1:Ly
                id = n(i,j)
                if j == 1
                    adj[id,1] = 0
                end
                if j == Ly
                    adj[id,3] =0
                end
            end
        end
        return adj
    end

    function adj_mat(L::Int,bc::BoundaryCondition;D)
        if D == 2
            return adj_mat(L,L,bc,D=2)
        elseif D==1
            return adj_mat(L,1,bc,D=1)
        else
            error("Invalid D")
        end
    end
    # L = 6
    # bc = boundary_condition("Free")
    # bc = boundary_condition("Periodic")
    # @show adj_mat(L,bc)
end
function hamiltonian(adj::Matrix)
    N = size(adj)[1]
    h = zeros(2N,2N)
    for id = CartesianIndices(adj)
        if adj[id] != 0
            if id[1] > adj[id]
                h[id[1],id[1]] += 0.5
                h[adj[id],adj[id]] += 0.5
                h[id[1]+N,id[1]+N] += 0.5
                h[adj[id]+N,adj[id]+N] += 0.5

                h[id[1]+N,adj[id]] += 0.5
                h[adj[id]+N,id[1]] += 0.5
                h[id[1],adj[id]+N] += 0.5
                h[adj[id],id[1]+N] += 0.5
            end
        end
    end
    return h
end

function initialize_givens_eigen_solver(ih::Matrix;perturbation::Int = 1000,hamiltonian_type = "Symmetry")
    N = Int(size(ih)[1]/2)
    η = diagm(vcat([1.0 for i=1:N],[-1.0 for i=1:N]))
    # S,V = eigen(η*ih)
    r = [i/N/1000 for i=1:N]
    r = vcat(r,r)
    # th = V*diagm(S+sort(vcat(-r,r)))*V^(-1)  #add perturbation
    th = ih + diagm(r)

    S,G = eigen(η*th)
    for i = 1:2N  # normalization
        c = transpose(G[:,i])*η*G[:,i]
        G[:,i] .=  G[:,i]./sqrt(abs(c))
    end

    for i = 1:N  # d~d^{\dagger} construction
        G[:,i] .= G[:,i+N]  #Why this construnction raise too large error?
    end
    if hamiltonian_type == "Symmetry"
        G[1:N,N+1:2N].= G[N+1:2N,1:N]
        G[N+1:2N,N+1:2N].= G[1:N,1:N]
    else hamiltonian_type == "AntiSymmetry"
        G[1:N,N+1:2N].= -G[N+1:2N,1:N]
        G[N+1:2N,N+1:2N].= G[1:N,1:N]
    end
    # the previous steps would raise small error, we need to use given eigen solver to eliminate them.
    return transpose(G)*ih*G,G
end

begin
	Lx = 4
	Ly = 3
	bc_type = "Cylinder"
    bc = boundary_condition(bc_type)
    N = Lx*Ly
end;

η = diagm(vcat([1.0 for i=1:N],[-1.0 for i=1:N]))
h = hamiltonian(adj_mat(Lx,Ly,bc))
s,G = initialize_givens_eigen_solver(h,hamiltonian_type="Symmetry")
print("-------initialize_givens_eigen_solver--------------\n")
@show norm(transpose(G)*hamiltonian(adj_mat(Lx,Ly,bc))*G-s)
@show norm(transpose(G)^(-1)*s*G^(-1)-hamiltonian(adj_mat(Lx,Ly,bc)))
@show norm(transpose(G)*η*G-η)
print("----------------------\n")

s,G = given_eigen_solver(h,hamiltonian_type = "Symmetry")
vscodedisplay(s)
vscodedisplay(G)
print("-------givens_eigen_solver--------------\n")
@show norm(transpose(G)*hamiltonian(adj_mat(Lx,Ly,bc))*G-s)
@show norm(transpose(G)^(-1)*s*G^(-1)-hamiltonian(adj_mat(Lx,Ly,bc)))
@show norm(transpose(G)*η*G-η)
print("----------------------\n")


text:
function get_provider()
    # Note: we CANNOT do something like have the `default` value be `get(ENV, "JULIA_FFTW_PROVIDER", "fftw")` here.
    # This is because the only way the Julia knows that a default has changed is if the values on-disk change; so
    # if your "default" value can be changed from the outside, you quickly run into cache invalidation issues.
    # So the default here _must_ be a constant.
    default_provider = "fftw"

    # Load the preference
    provider = @load_preference("provider", default_provider)

    # Ensure the provider matches one of the ones we support
    if provider ∉ ("fftw", "mkl")
        @error("Invalid provider setting \"$(provider)\"; valid settings include [\"fftw\", \"mkl\"], defaulting to \"fftw\"")
        provider = default_provider
    end
    return provider
end

# Read in preferences, see if any users have requested a particular backend
const fftw_provider = get_provider()

# We'll initialize `libfftw3` here (in the conditionals below), and
# it will get overwritten again in `__init__()`.  This allows us to
# `ccall` at build time, and also be relocatable for PackageCompiler.
const libfftw3 = Ref{String}()
const libfftw3f = Ref{String}()

"""
    set_provider!(provider; export_prefs::Bool = false)

Convenience wrapper for setting the FFT provider.  Valid values include `"fftw"`, `"mkl"`.
Also supports `Preferences` sentinel values `nothing` and `missing`; see the docstring for
`Preferences.set_preferences!()` for more information on what these values mean.
"""
function set_provider!(provider; export_prefs::Bool = false)
    if provider !== nothing && provider !== missing && provider ∉ ("fftw", "mkl")
        throw(ArgumentError("Invalid provider value '$(provider)'"))
    end
    set_preferences!(@__MODULE__, "provider" => provider; export_prefs, force = true)
    if provider != fftw_provider
        # Re-fetch to get default values in the event that `nothing` or `missing` was passed in.
        provider = get_provider()
        @info("FFTW provider changed; restart Julia for this change to take effect", provider)
    end
end

# If we're using fftw_jll, load it in
@static if fftw_provider == "fftw"
    import FFTW_jll
    libfftw3[] = FFTW_jll.libfftw3_path
    libfftw3f[] = FFTW_jll.libfftw3f_path

    # callback function that FFTW uses to launch `num` parallel
    # tasks (FFTW/fftw3#175):
    function spawnloop(f::Ptr{Cvoid}, fdata::Ptr{Cvoid}, elsize::Csize_t, num::Cint, callback_data::Ptr{Cvoid})
        @sync for i = 0:num-1
            Threads.@spawn ccall(f, Ptr{Cvoid}, (Ptr{Cvoid},), fdata + elsize*i)
        end
    end

    # If FFTW was built with threads, then they must be initialized before any FFTW planning routine.
    #   -- This initializes FFTW's threads support (defaulting to 1 thread).
    #      If this isn't called before the FFTW planner is created, then
    #      FFTW's threads algorithms won't be registered or used at all.
    #      (Previously, we called fftw_cleanup, but this invalidated existing
    #       plans, causing Base Julia issue #19892.)
    function fftw_init_threads()
        stat =  ccall((:fftw_init_threads,   libfftw3[]), Int32, ())
        statf = ccall((:fftwf_init_threads, libfftw3f[]), Int32, ())
        if stat == 0 || statf == 0
            error("could not initialize FFTW threads")
        end

        if nthreads() > 1
            cspawnloop = @cfunction(spawnloop, Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}, Csize_t, Cint, Ptr{Cvoid}))
            ccall((:fftw_threads_set_callback,  libfftw3[]),  Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), cspawnloop, C_NULL)
            ccall((:fftwf_threads_set_callback, libfftw3f[]), Cvoid, (Ptr{Cvoid}, Ptr{Cvoid}), cspawnloop, C_NULL)
        end
    end
end

# If we're using MKL, load it in and set library paths appropriately.
@static if fftw_provider == "mkl"
    import MKL_jll
    libfftw3[] = MKL_jll.libmkl_rt_path
    libfftw3f[] = MKL_jll.libmkl_rt_path
    const _last_num_threads = Ref(Cint(1))
end


text:
using Knet: Knet, gpu, conv4, value

include(Knet.dir("data","imagenet.jl"))

mutable struct ResNet
    w
    ms
    avgimg
    atype
end

function ResNet(;name="imagenet-resnet-101-dag", atype=nothing)
    if atype == nothing
        atype = gpu() >= 0 ? KnetArray{Float32} : Array{Float32}
    end

    model = matconvnet(name)

    avgimg = model["meta"]["normalization"]["averageImage"]
    avgimg = convert(Array{Float32}, avgimg)

    w, ms = get_params(model["params"], atype)
    ResNet(w, ms, avgimg, atype)
end

function (r::ResNet)(path::String)
    img = imgdata(path, r.avgimg)
    x = convert(r.atype, img)
    resnet101(r.w, x, r.ms)
end

function resnet(path)
    atype = gpu()>=0 ? KnetArray{Float64} : Array{Float64}
    modelname = "imagenet-resnet-101-dag"
    model = matconvnet(modelname)
    avgimg = model["meta"]["normalization"]["averageImage"]
    avgimg = convert(Array{Float64}, avgimg)
    img = imgdata(path, avgimg)
    x = convert(atype, img)

    w, ms = get_params(model["params"], atype)
    return resnet101(w,x,ms)
end

function resnet101(w,x,ms; mode=1)
    # layer 1
    conv1 = reslayerx1(w[1:3],x,ms; padding=3, stride=2, mode=mode)
    pool1 = pool(conv1; window=3, stride=2)

    # layer 2,3,4,5
    r2 = reslayerx5(w[4:33], pool1, ms; strides=[1,1,1,1], mode=mode)
    r3 = reslayerx5(w[34:72], r2, ms; mode=mode)
    r4 = reslayerx5(w[73:282], r3, ms; mode=mode)
    r5 = reslayerx5(w[283:312], r4, ms; mode=mode)

    windowsize = size(r5, 1)
    # fully connected layer
    pool5  = pool(r5; stride=1, window=windowsize, mode=2)
    return mat(pool5)
    #fc1000 = w[313] * mat(pool5) .+ w[314]
end

function batchnorm(w, x, ms; mode=1, epsilon=1e-5)
    mu, sigma = nothing, nothing
    if mode == 0
        d = ndims(x) == 4 ? (1,2,4) : (2,)
        s = prod(size(x,d...))
        mu = sum(x,d) / s
        x0 = x .- mu
        x1 = x0 .* x0
        sigma = sqrt(epsilon + (sum(x1, d)) / s)
    elseif mode == 1
        mu = popfirst!(ms)
        sigma = popfirst!(ms)
    end

    # we need value in backpropagation
    push!(ms, value(mu), value(sigma))
    xhat = (x.-mu) ./ sigma
    return w[1] .* xhat .+ w[2]
end

function reslayerx0(w,x,ms; padding=0, stride=1, mode=1)
    b  = conv4(w[1],x; padding=padding, stride=stride)
    bx = batchnorm(w[2:3],b,ms; mode=mode)
end

function reslayerx1(w,x,ms; padding=0, stride=1, mode=1)
    relu.(reslayerx0(w,x,ms; padding=padding, stride=stride, mode=mode))
end

function reslayerx2(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    ba = reslayerx1(w[1:3],x,ms; padding=pads[1], stride=strides[1], mode=mode)
    bb = reslayerx1(w[4:6],ba,ms; padding=pads[2], stride=strides[2], mode=mode)
    bc = reslayerx0(w[7:9],bb,ms; padding=pads[3], stride=strides[3], mode=mode)
end

function reslayerx3(w,x,ms; pads=[0,0,1,0], strides=[2,2,1,1], mode=1) # 12
    a = reslayerx0(w[1:3],x,ms; stride=strides[1], padding=pads[1], mode=mode)
    b = reslayerx2(w[4:12],x,ms; strides=strides[2:4], pads=pads[2:4], mode=mode)
    relu.(a .+ b)
end

function reslayerx4(w,x,ms; pads=[0,1,0], strides=[1,1,1], mode=1)
    relu.(x .+ reslayerx2(w,x,ms; pads=pads, strides=strides, mode=mode))
end

function reslayerx5(w,x,ms; strides=[2,2,1,1], mode=1)
    x = reslayerx3(w[1:12],x,ms; strides=strides, mode=mode)
    for k = 13:9:length(w)
        x = reslayerx4(w[k:k+8],x,ms; mode=mode)
    end
    return x
end

function get_params(params, atype)
    len = length(params["value"])
    ws, ms = [], []
    for k = 1:len
        name = params["name"][k]
        value = convert(Array{Float32}, params["value"][k])

        if endswith(name, "moments")
            push!(ms, reshape(value[:,1], (1,1,size(value,1),1)))
            push!(ms, reshape(value[:,2], (1,1,size(value,1),1)))
        elseif startswith(name, "bn")
            push!(ws, reshape(value, (1,1,length(value),1)))
        elseif startswith(name, "fc") && endswith(name, "filter")
            push!(ws, transpose(reshape(value,(size(value,3),size(value,4)))))
        elseif startswith(name, "conv") && endswith(name, "bias")
            push!(ws, reshape(value, (1,1,length(value),1)))
        else
            push!(ws, value)
        end
    end
    map(wi->convert(atype, wi), ws),
    map(mi->convert(atype, mi), ms)
end

max_stars_repo_name: alperkesen/AutoML

max_issues_repo_name: alperkesen/AutoML


text:
if isdefined(Base, :Grisu)
    import Base.Grisu
else
    import Grisu
end

function tostring(x::AbstractFloat, n::Int=15)
    #
    # Based on Base.Grisu._show(io::IO, x::AbstractFloat, mode, n::Int, typed, compact)
    # See julia/base/grisu/grisu.jl
    #
    mode = Grisu.SHORTEST
    #const typed = false
    #const compact = false

    isnan(x) && return "NaN"
    if isinf(x)
        if signbit(x)
            return "-Inf"
        else
            return "Inf"
        end
    end

    (len,pt,neg), buffer = Grisu.grisu(x,mode,n), Grisu.DIGITS
    pdigits = pointer(buffer)

    #=
    if mode == PRECISION
        while len > 1 && buffer[len] == 0x30
            len -= 1
        end
    end
    =#

    io = IOBuffer()
    if neg
        write(io, "-")
    end

    #=
    exp_form = pt <= -4 || pt > 6
    exp_form = exp_form || (pt >= len && abs(mod(x + 0.05, 10^(pt - len)) - 0.05) > 0.05) # see issue #6608
    if exp_form # .00001 to 100000.
        # => #.#######e###
        unsafe_write(io, pdigits, 1)
        write(io, '.')
        if len > 1
            unsafe_write(io, pdigits+1, len-1)
        else
            write(io, '0')
        end
        write(io, (typed && isa(x,Float32)) ? 'f' : 'e')
        write(io, dec(pt-1))
        typed && isa(x,Float16) && write(io, ")")
        return
    elseif pt <= 0
    =#
    if pt <= 0
        # => 0.00########
        write(io, "0.")
        while pt < 0
            write(io, '0')
            pt += 1
        end
        unsafe_write(io, pdigits, len)
    elseif pt >= len
        # => ########00.0
        unsafe_write(io, pdigits, len)
        while pt > len
            write(io, '0')
            len += 1
        end
        write(io, ".0")
    else # => ####.####
        unsafe_write(io, pdigits, pt)
        write(io, '.')
        unsafe_write(io, pdigits+pt, len-pt)
    end
    #typed && !compact && isa(x,Float32) && write(io, "f0")
    #typed && isa(x,Float16) && write(io, ")")
    return String(take!(io))
end


text:
using Plots
using Distributions
using Random
gr(size=(300,300))

function ellipse!(μ,Σ;kw...)
    for i=1:3
        plot!([tuple((i*sqrt(Σ)*[cos(t),sin(t)]+μ)...) for t in range(0,stop=2π,length=100)];
              seriestype=:shape,linealpha=0.5,linecolor=:gray,fillalpha=[0.4,0.2,0.1][i],kw...)
    end
    current()
end

function mixtureplot(X₁,X₂,μ₀,Σ₀,μ₁,Σ₁,Π)
    scatter(X₁,X₂;fillalpha=0.5,markerstrokewidth=0.5,
        marker_z=Π,mc=ColorGradient([:lightblue,:orange]),
        colorbar=:false)
    ellipse!(μ₀,Σ₀,fillcolor=:lightblue)
    ellipse!(μ₁,Σ₁,fillcolor=:orange)
    plot!(;bg=:transparent,xlims=(-4.5,8.5),ylims=(-3.5,10.5),
        leg=false,ticks=:none,ratio=:equal)
end

function scalein(x)
    (x-1/2)^101/(0.5^101)
end

Random.seed!(123);
n = 100
α = 0.4
𝒩₀ = MvNormal([1,1],[2.0 1.0; 1.0 2.0])
𝒩₁ = MvNormal([3.0,7.0],[1.5 0; 0 0.5])
X₁ = zeros(n)
X₂ = zeros(n)
Z = zeros(Bool,n)
for i=1:n
    Z[i] = rand(Bernoulli(α))
    X₁[i],X₂[i] = Z[i] ? rand(𝒩₁) : rand(𝒩₀)
end

α = 0.6
μ₀ = [3.0,3.0]
μ₁ = [1.0,6.0]
Σ₀ = 1.0*Matrix(I, 2, 2)
Σ₁ = 1.0*Matrix(I, 2, 2)
mixtureplot(X₁,X₂,μ₀,Σ₀,μ₁,Σ₁)

Π = [α*pdf(MvNormal(μ₁,Σ₁),[x₁,x₂]) /
       ((1-α)*pdf(MvNormal(μ₀,Σ₀),[x₁,x₂]) +
        α*pdf(MvNormal(μ₁,Σ₁),[x₁,x₂])) for (x₁,x₂) in zip(X₁,X₂)];

α = sum(Π)/n
μ₀ = [(1 .- Π) ⋅ X₁, (1 .- Π) ⋅ X₂] / sum(1 .- Π)
μ₁ = [Π ⋅ X₁, Π ⋅ X₂] / sum(Π)
Σ₀ = Matrix(Hermitian(sum((1-π)*([x₁,x₂] - μ₀) * ([x₁,x₂] - μ₀)' for (x₁,x₂,π) in zip(X₁,X₂,Π))/sum(1 .- Π)))
Σ₁ = Matrix(Hermitian(sum(π*([x₁,x₂] - μ₁) * ([x₁,x₂] - μ₁)' for (x₁,x₂,π) in zip(X₁,X₂,Π))/sum(Π)))
Π = [α*pdf(MvNormal(μ₁,Σ₁),[x₁,x₂]) /
       ((1-α)*pdf(MvNormal(μ₀,Σ₀),[x₁,x₂]) +
        α*pdf(MvNormal(μ₁,Σ₁),[x₁,x₂])) for (x₁,x₂) in zip(X₁,X₂)];

mixtureplot(X₁,X₂,μ₀,Σ₀,μ₁,Σ₁,Π)

max_stars_repo_name: seanrattana/courses

max_issues_repo_name: seanrattana/courses

max_forks_repo_name: seanrattana/courses


text:
abstract Benchmarkable

immutable Benchmark <: Benchmarkable
    name :: String
    ntimes :: Integer
    run :: Function
end
Base.show(b::Benchmark) = let name = b.name
			  print("Benchmark $name")
			  end

immutable BenchmarkGroup <: Benchmarkable
    name :: String
    benchmarks :: Vector{Benchmark}
end
Base.show(b::BenchmarkGroup) = let name = b.name
			       print("Benchmark Group $name")
			       end
function bench(name::String,
	       ntimes::Integer,
	       action::Function)
    Benchmark(name, ntimes, action)
end

function bgroup(name::String,
		benchmarks::Vector{Benchmarkable})
    BenchMarkGroup(name, benchmarks)
end

typealias Sample Vector{Float64}

immutable Payload
    sample   :: Sample
    analysis :: SampleAnalysis
    outliers :: Outliers
end

immutable Result
    info :: String
    payload :: Payload
end


text:
"""
    newton(f, fprime, x0, e, n)

Compute the root of the given function using Newton's method.
Uses the point at which tangent line intersects x-axis to find the next approximate of the root of function 'f',
as such it requires initial approximation `x0` and `fprime`, the derivative of function `f`.
Newton's method has very fast convergence provided a sufficiantly accurate initial approximation is chosen.

    newton(f::Function, fprime::Function, x0::Float64, e::Float64, n::Integer)

### OUTPUT:
  * `p` - approximation of the root of the function `f` (Float64) or,
  * `-1` - if `n` the maximum number of iterations was exceeded.
"""
function newton(f::Function, fprime::Function, x0::Float64, e::Float64, n::Integer)

    Afx = []

    for i in 1:n

        yprime = fprime(x0)

        if abs(yprime) < e
            warn("First derivative is zero")
            return x0
        end

        y = f(x0)
        x1 = x0 - y / yprime
        push!(Afx, f(x1))

        if abs(x1 - x0) < e
            return x1
        end

        x0 = x1
    end

    return -1
end


text:
import Base.∈

export constrained_dimensions,
       tosimplehrep,
       remove_redundant_constraints,
       remove_redundant_constraints!,
       linear_map,
       chebyshev_center

# default LP solver for floating-point numbers
function default_lp_solver(N::Type{<:AbstractFloat})
    GLPKSolverLP(method=:Simplex)
end

# default LP solver for rational numbers
function default_lp_solver(N::Type{<:Rational})
    GLPKSolverLP(method=:Exact)
end

# fallback method
function default_polyhedra_backend(P, N)
    require(:Polyhedra; fun_name="default_polyhedra_backend")
    error("no default backend for numeric type $N")
end

"""
    ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N})::Bool where {N<:Real}

Check whether a given point is contained in a polyhedron.

### Input

- `x` -- point/vector
- `P` -- polyhedron

### Output

`true` iff ``x ∈ P``.

### Algorithm

This implementation checks if the point lies inside each defining half-space.
"""
function ∈(x::AbstractVector{N}, P::AbstractPolyhedron{N})::Bool where {N<:Real}
    @assert length(x) == dim(P) "a $(length(x))-dimensional point cannot be " *
        "an element of a $(dim(P))-dimensional set"

    for c in constraints_list(P)
        if dot(c.a, x) > c.b
            return false
        end
    end
    return true
end

"""
    constrained_dimensions(P::AbstractPolyhedron)::Vector{Int} where {N<:Real}

Return the indices in which a polyhedron is constrained.

### Input

- `P` -- polyhedron

### Output

A vector of ascending indices `i` such that the polyhedron is constrained in
dimension `i`.

### Examples

A 2D polyhedron with constraint ``x1 ≥ 0`` is constrained in dimension 1 only.
"""
function constrained_dimensions(P::AbstractPolyhedron)::Vector{Int}
    constraints = constraints_list(P)
    if isempty(constraints)
        return Int[]
    end
    zero_indices = zeros(Int, dim(P))
    for constraint in constraints
        for i in constrained_dimensions(constraint)
            zero_indices[i] = i
        end
    end
    return filter(x -> x != 0, zero_indices)
end

"""
    tosimplehrep(constraints::AbstractVector{LC})
        where {N<:Real, LC<:LinearConstraint{N}}

Return the simple H-representation ``Ax ≤ b`` from a list of linear constraints.

### Input

- `constraints` -- a list of linear constraints

### Output

The tuple `(A, b)` where `A` is the matrix of normal directions and `b` is the
vector of offsets.
"""
function tosimplehrep(constraints::AbstractVector{LC}
                     ) where {N<:Real, LC<:LinearConstraint{N}}
    n = length(constraints)
    if n == 0
        A = Matrix{N}(undef, 0, 0)
        b = Vector{N}(undef, 0)
        return (A, b)
    end
    A = zeros(N, n, dim(first(constraints)))
    b = zeros(N, n)
    @inbounds begin
        for (i, Pi) in enumerate(constraints)
            A[i, :] = Pi.a
            b[i] = Pi.b
        end
    end
    return (A, b)
end

"""
     remove_redundant_constraints!(
         constraints::AbstractVector{<:LinearConstraint{N}};
         [backend]=default_lp_solver(N))::Bool where {N<:Real}

Remove the redundant constraints of a given list of linear constraints; the list
is updated in-place.

### Input

- `constraints` -- list of constraints
- `backend`     -- (optional, default: `default_lp_solver(N)`) the backend used
                   to solve the linear program
### Output

`true` if the function was successful and the list of constraints `constraints`
is modified by removing the redundant constraints, and `false` only if the
constraints are infeasible.

### Notes

Note that the result may be `true` even if the constraints are infeasible.
For example, ``x ≤ 0 && x ≥ 1`` will return `true` without removing any
constraint.
To check if the constraints are infeasible, use
`isempty(HPolyhedron(constraints)`.

### Algorithm

If there are `m` constraints in `n` dimensions, this function checks one by one
if each of the `m` constraints is implied by the remaining ones.

To check if the `k`-th constraint is redundant, an LP is formulated using the
constraints that have not yet being removed.
If, at an intermediate step, it is detected that a subgroup of the constraints
is infeasible, this function returns `false`.
If the calculation finished successfully, this function returns `true`.

For details, see [Fukuda's Polyhedra
FAQ](https://www.cs.mcgill.ca/~fukuda/soft/polyfaq/node24.html).
"""
function remove_redundant_constraints!(
        constraints::AbstractVector{<:LinearConstraint{N}};
        backend=default_lp_solver(N))::Bool where {N<:Real}

    A, b = tosimplehrep(constraints)
    m, n = size(A)
    non_redundant_indices = 1:m

    i = 1 # counter over reduced constraints

    for j in 1:m    # loop over original constraints
        α = A[j, :]
        Ar = A[non_redundant_indices, :]
        br = b[non_redundant_indices]
        br[i] = b[j] + one(N)
        lp = linprog(-α, Ar, '<', br, -Inf, Inf, backend)
        if lp.status == :Infeasible
            # the polyhedron is empty
            return false
        elseif lp.status == :Optimal
            objval = -lp.objval
            if _leq(objval, b[j])
                # the constraint is redundant
                non_redundant_indices = setdiff(non_redundant_indices, j)
            else
                # the constraint is not redundant
                i = i+1
            end
        else
            error("LP is not optimal; the status of the LP is $(lp.status)")
        end
    end

    deleteat!(constraints, setdiff(1:m, non_redundant_indices))
    return true
end

"""
    remove_redundant_constraints(
        constraints::AbstractVector{<:LinearConstraint{N}};
        backend=default_lp_solver(N)) where {N<:Real}

Remove the redundant constraints of a given list of linear constraints.

### Input

- `constraints` -- list of constraints
- `backend`     -- (optional, default: `default_lp_solver(N)`) the backend used
                   to solve the linear program
### Output

The list of constraints with the redundant ones removed, or an empty set if the
constraints are infeasible.

### Algorithm

See `remove_redundant_constraints!(::AbstractVector{<:LinearConstraint})` for
details.
"""
function remove_redundant_constraints(
        constraints::AbstractVector{<:LinearConstraint{N}};
        backend=default_lp_solver(N)) where {N<:Real}
    constraints_copy = copy(constraints)
    if remove_redundant_constraints!(constraints_copy, backend=backend)
        return constraints_copy
    else  # the constraints are infeasible
        return EmptySet{N}()
    end
end

"""
    linear_map(M::AbstractMatrix{N},
               P::AbstractPolyhedron{N};
               [algorithm]::Union{String, Nothing}=nothing,
               [check_invertibility]::Bool=true,
               [cond_tol]::Number=DEFAULT_COND_TOL,
               [inverse]::Union{AbstractMatrix{N}, Nothing}=nothing
              ) where {N<:Real}

Concrete linear map of a polyhedral set.

### Input

- `M`         -- matrix
- `P`         -- polyhedral set
- `algorithm` -- (optional; default: `nothing`, see the Notes section below)
                 algorithm to be used; possible choices are:
  - `"vrep"`     -- apply the linear map to each vertex of `P` (note that this
                    only works for *polytopes*)
  - `"inverse"`  -- compute the matrix inverse and apply it to each constraint
                    of `P`
  - `"division"` -- divide each constraint of `P` by `M` from the left
- `check_invertibility` -- (optional, deault: `true`) check if the linear map is
                           invertible, in which case this function uses the
                           matrix inverse; if this flag is set to `false`, we
                           assume that the matrix is invertible; otherwise, if
                           the invertibility check fails, we fall back to the
                           `"vrep"` algorithm; this option is ignored if the
                           `"vrep"` algorithm is used
- `cond_tol`  -- (optional; default: `DEFAULT_COND_TOL`) tolerance of matrix
                 condition (used to check whether the matrix is invertible)
- `inverse`   -- (optional; default: `nothing`) matrix inverse (if known)

### Output

The type of the result is "as close as possible" to the the type of `P`.
Let `(m, n)` be the size of `M`, where `m ≠ n` is allowed for rectangular maps.

To fix the type of the output to something different than the default value,
consider post-processing the result of this function with a call to a suitable
`convert` method.

In particular, the output depends on the type of `P`, on `m`, and the algorithm
that was used:

- If the vertex-based approach was used:

    - If `P` is a `VPolygon` and `m = 2` then the output is a `VPolygon`.
    - If `P` is a `VPolytope` then the output is a `VPolytope`.
    - Otherwise, the output is an `Interval` if `m = 1`, a `VPolygon` if `m = 2`
      and a `VPolytope` in other cases.

- If the invertibility criterion was used:

    - The types of `HalfSpace`, `Hyperplane`, `Line` and `AbstractHPolygon` are
      preserved.
    - If `P` is an `AbstractPolytope`, then the output is an `Interval` if `m = 1`,
      an `HPolygon` if `m = 2` and an `HPolytope` in other cases.
    - Otherwise, the output is an `HPolyhedron`.

### Notes

If `algorithm` is `nothing` (the default), we determine the algorithm as
follows.
If `M` is sparse or non-square, and has full row rank (which we check), we use
`"division"`.
If `M` is dense and square and invertible (either assumed or checked; see option
`check_invertibility`), we use `"inverse"`.
Otherwise, we use `"vrep"`.

### Algorithm

This function implements two algorithms for the linear map:

- If the matrix ``M`` is invertible (which we check via a sufficient condition),
  then ``y = M x`` implies ``x = \\text{inv}(M) y`` and we transform the
  constraint system ``A x ≤ b`` to ``A \\text{inv}(M) y ≤ b``.
- Otherwise, we transform the polyhedron to vertex representation and apply the
  map to each vertex, returning a polyhedron in vertex representation.

Note that the vertex representation (second approach) is only available if the
polyhedron is bounded, which we check.

If the matrix is known to be invertible, the the option `check_invertibility`
can be used to skip the invertibility test.
Note that we only check for invertibility if the `algorithm` is either
unspecified or specified as `"inverse"`.
If the matrix inverse is even known, it can be specified with the option
`inverse`, in which case we ignore the other options and also the original
matrix `M`.

The algorithms `"division"` and `"inverse"` give control about whether the full
matrix inverse is computed or only the left division on the normal vectors is
used.
Note that this helps as a workaround when `M` is sparse (since the `inv`
function is not available for sparse matrices) or rectangular.
For sparse and invertible matrices, either use the algorithm `"division"` or
convert `M` to a dense matrix (as in `linear_map(Matrix(M), P)`).

Internally, this function operates on the level of the `AbstractPolyhedron`
interface, but the actual algorithm uses dispatch on the concrete type of `P`,
depending on the algorithm that is used:

- `_linear_map_vrep(M, P)` if the vertex approach is used
- `_linear_map_hrep(M, P, use_inv)` if the invertibility criterion is used,
  where `use_inv` is determined by the algorithm

New subtypes of the interface should define their own `_linear_map_vrep`
(resp. `_linear_map_hrep`) for special handling of the linear map; otherwise
the fallback implementation for `AbstractPolyhedron` is used.
"""
function linear_map(M::AbstractMatrix{N},
                    P::AbstractPolyhedron{N};
                    algorithm::Union{String, Nothing}=nothing,
                    check_invertibility::Bool=true,
                    cond_tol::Number=DEFAULT_COND_TOL,
                    inverse::Union{AbstractMatrix{N}, Nothing}=nothing
                   ) where {N<:Real}
    # if `inverse` is specified, we ignore other inputs and do not check them
    if inverse != nothing
        return _linear_map_hrep(M, P, true; inverse=inverse)
    end

    @assert dim(P) == size(M, 2) "a linear map of size $(size(M)) cannot be " *
        "applied to a set of dimension $(dim(P))"

    if algorithm == nothing
        if issparse(M) || !issquare(M)
            if hasfullrowrank(M)
                algorithm = "division"
            else
                algorithm = "vrep"
            end
        elseif !check_invertibility || isinvertible(M; cond_tol=cond_tol)
            algorithm = "inverse"
        else
            algorithm = "vrep"
        end
    elseif algorithm == "inverse"
        # check invertibility
        if check_invertibility && !isinvertible(M; cond_tol=cond_tol)
            throw(ArgumentError("algorithm \"inverse\" requires an " *
                                "invertible matrix"))
        end
    end

    if algorithm == "vrep"
        return _linear_map_vrep(M, P)
    else
        if algorithm == "inverse"
            use_inv = true
        elseif algorithm == "division"
            use_inv = false
        else
            throw(ArgumentError("got unknown algorithm \"$algorithm\"; " *
                "available choices: \"vrep\", \"inverse\", \"division\""))
        end
        return _linear_map_hrep(M, P, use_inv)
    end
end

# handle different numeric types
function linear_map(M::AbstractMatrix{NM},
                    P::AbstractPolyhedron{NP};
                    kwargs...) where {NM<:Real, NP<:Real}
    N = promote_type(NM, NP)
    if N != NP
        error("conversion between numeric types of polyhedra not implemented " *
            "yet (see #1181)")
    else
        return linear_map(N.(M), P; kwargs...)
    end
end

function _linear_map_vrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N}) where {N<:Real}
    if !isbounded(P)
        throw(ArgumentError("the linear map in vertex representation for an " *
            "unbounded set is not defined"))
    end
    require(:Polyhedra; fun_name="linear_map",
            explanation="of a $(typeof(P)) by a non-invertible matrix")
    # since P is bounded, we pass an HPolytope and then convert it to vertex representation
    P = tovrep(HPolytope(constraints_list(P), check_boundedness=false))
    return _linear_map_vrep(M, P)
end

function _linear_map_hrep(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                          use_inv::Bool;
                          inverse::Union{Nothing, AbstractMatrix{N}}=nothing
                         ) where {N<:Real}
    constraints = _linear_map_hrep_helper(M, P, use_inv, inverse=inverse)
    return HPolyhedron(constraints)
end

function _linear_map_hrep_helper(M::AbstractMatrix{N}, P::AbstractPolyhedron{N},
                                 use_inv::Bool;
                                 inverse::Union{Nothing, AbstractMatrix{N}}=nothing
                                ) where {N<:Real}
    constraints_P = constraints_list(P)
    constraints_MP = similar(constraints_P)
    if use_inv || inverse != nothing
        if inverse == nothing
            inverse = inv(M)
        end
        @inbounds for (i, c) in enumerate(constraints_P)
            constraints_MP[i] =
                LinearConstraint(vec(_At_mul_B(c.a, inverse)), c.b)
        end
    else
        # take left division for each constraint c, transpose(M) \ c.a
        @inbounds for (i, c) in enumerate(constraints_P)
            constraints_MP[i] = LinearConstraint(_At_ldiv_B(M, c.a), c.b)
        end
    end
    return constraints_MP
end

"""
    plot_recipe(P::AbstractPolyhedron{N}, [ε]::N=zero(N)) where {N<:Real}

Convert a (bounded) polyhedron to a pair `(x, y)` of points for plotting.

### Input

- `P` -- bounded polyhedron
- `ε` -- (optional, default: `0`) ignored, used for dispatch

### Output

A pair `(x, y)` of points that can be plotted.

### Algorithm

We first assert that `P` is bounded (i.e., that `P` is a polytope).

One-dimensional polytopes are converted to an `Interval`.
Three-dimensional or higher-dimensional polytopes are not supported.

For two-dimensional polytopes (i.e., polygons) we compute their set of vertices
using `vertices_list` and then plot the convex hull of these vertices.
"""
function plot_recipe(P::AbstractPolyhedron{N}, ε::N=zero(N)) where {N<:Real}
    @assert dim(P) <= 2 "cannot plot a $(dim(P))-dimensional $(typeof(P))"
    @assert isbounded(P) "cannot plot an unbounded $(typeof(P))"

    if dim(P) == 1
        return plot_recipe(convert(Interval, P), ε)
    else
        vlist = transpose(hcat(convex_hull(vertices_list(P))...))
        if isempty(vlist)
            @warn "received a polyhedron with no vertices during plotting"
            return plot_recipe(EmptySet{N}(), ε)
        end
        x, y = vlist[:, 1], vlist[:, 2]

        if length(x) > 1
            # add first vertex to "close" the polygon
            push!(x, x[1])
            push!(y, y[1])
        end

        return x, y
    end
end

"""
    chebyshev_center(P::AbstractPolyhedron{N};
                     [get_radius]::Bool=false,
                     [backend]=default_polyhedra_backend(P, N),
                     [solver]=JuMP.with_optimizer(GLPK.Optimizer)
                     ) where {N<:AbstractFloat}

Compute the [Chebyshev center](https://en.wikipedia.org/wiki/Chebyshev_center)
of a polytope.

### Input

- `P`          -- polytope
- `get_radius` -- (optional; default: `false`) option to additionally return the
                  radius of the largest ball enclosed by `P` around the
                  Chebyshev center
- `backend`    -- (optional; default: `default_polyhedra_backend(P, N)`) the
                  backend for polyhedral computations
- `solver`     -- (optional; default: `JuMP.with_optimizer(GLPK.Optimizer)`) the
                  LP solver passed to `Polyhedra`

### Output

If `get_radius` is `false`, the result is the Chebyshev center of `P`.
If `get_radius` is `true`, the result is the pair `(c, r)` where `c` is the
Chebyshev center of `P` and `r` is the radius of the largest ball with center
`c` enclosed by `P`.

### Notes

The Chebyshev center is the center of a largest Euclidean ball enclosed by `P`.
In general, the center of such a ball is not unique (but the radius is).
"""
function chebyshev_center(P::AbstractPolyhedron{N};
                          get_radius::Bool=false,
                          backend=default_polyhedra_backend(P, N),
                          solver=JuMP.with_optimizer(GLPK.Optimizer)
                         ) where {N<:AbstractFloat}
    require(:Polyhedra; fun_name="chebyshev_center")
    # convert to HPolyhedron to ensure `polyhedron` is applicable (see #1505)
    Q = polyhedron(convert(HPolyhedron, P); backend=backend)
    c, r = Polyhedra.chebyshevcenter(Q, solver)

    if get_radius
        return c, r
    end
    return c
end