text:
export append, append!

"""
    append!(data::PopData, data2::PopData)
Add the rows of `data2` to the end of `data`. This will add the samples present
in the second `PopData` object to the first `PopData` object (mutating it).
**Note** that this is a simple appending, and you risk corrupting your `PopData` if
the two `PopData` objects do not have identical loci.

**Example**
```
julia> cats = nancycats()
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 237
  Loci: 9
  Populations: 17
  Longitude: absent
  Latitude: absent

julia> purrfect_pairs = cross(cats, "N200", "N7", generation = "F1")
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 100
  Loci: 9
  Populations: 1
  Longitude: absent
  Latitude: absent

julia> append!(cats, purrfect_pairs);

julia> cats
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 337
  Loci: 9
  Populations: 18
  Longitude: absent
  Latitude: absent
```
"""
function Base.append!(data::PopData, data2::PopData)
    if "parents" ∉ names(data.meta) && "parents" ∈ names(data2.meta)
        len = length(data.meta.name)
        insertcols!(
            data.meta,
            :parents => Vector{Union{Missing, Tuple{String,String}}}(undef, len)
            )
    elseif "parents" ∉ names(data2.meta) && "parents" ∈ names(data.meta)
        len = length(data2.meta.name)
        insertcols!(
            data2.meta,
            :parents => Vector{Union{Missing, Tuple{String,String}}}(undef, len)
            )
    end

    append!(data.meta, data2.meta)

    data.loci.name = decompress(data.loci.name)
    append!(data.loci, data2.loci)
    return data
end

"""
    append(data::PopData, data2::PopData)
Add the rows of `data2` to the end of `data`. This will combine the samples present
in both `PopData` objects and return a new `PopData` object. **Note** that this is
a simple appending, and you risk corrupting your `PopData` if the two `PopData`
objects do not have identical loci.

**Example**
```
julia> cats = nancycats()
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 237
  Loci: 9
  Populations: 17
  Longitude: absent
  Latitude: absent

julia> purrfect_pairs = cross(cats, "N200", "N7", generation = "F1")
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 100
  Loci: 9
  Populations: 1
  Longitude: absent
  Latitude: absent

julia> merged_cats = append(cats, purrfect_pairs)
PopData Object
  Markers: Microsatellite
  Ploidy: 2
  Samples: 337
  Loci: 9
  Populations: 18
  Longitude: absent
  Latitude: absent
```
"""
function append(data::PopData, data2::PopData)
    tmp = PopData(copy(data.meta), copy(data.loci))
    tmp2 = data2
    if "parents" ∉ names(data.meta) && "parents" ∈ names(data2.meta)
        len = length(tmp.meta.name)
        insertcols!(
            tmp.meta,
            :parents => Vector{Union{Missing, Tuple{String,String}}}(undef, len)
            )
    elseif "parents" ∉ names(data2.meta) && "parents" ∈ names(data.meta)
        tmp2 = PopData(copy(data2.meta), copy(data2.loci))
        len = length(tmp2.meta.name)
        insertcols!(
            tmp2.meta,
            :parents => Vector{Union{Missing, Tuple{String,String}}}(undef, len)
            )
    end

    append!(tmp.meta, tmp2.meta)

    tmp.loci.name = decompress(tmp.loci.name)
    append!(tmp.loci, tmp2.loci)
    return tmp
end

function allele_pool(locus::T) where T <: GenoArray
  Tuple(Base.Iterators.flatten(skipmissing(locus)))
end

function Base.sort(x::NTuple{N,T}) where N where T <: Signed
  Tuple(sort(SVector(x)))
end

function allele_pool(data::PopData)
  # index dataframe by locus
  idx_df = groupby(data.loci, [:locus])
  # instantiate dict to store alleles
  allele_dict = Dict{String,NTuple}()
  # pull out loci names
  loc = getindex.(keys(idx_df), :locus)
  [allele_dict[i] = allele_pool(idx_df[(;locus = i)].genotype) for i in loc]
  return String.(loc), allele_dict
end

"""
```
simulate_sample(alleles::Dict{String,NTuple}, loc::Vector{String}; ploidy::Int)
```
Using a global allele pool given by a Dict{loci,alleles} and a list of loci (`loc`), simulate
an individual with a given `ploidy`. Returns a Vector of genotypes.

**Example**
```
julia> cats = nancycats() ;
julia> loc, alleles = allele_pool(cats) ;
julia> simulate_parent(alleles, loc, ploidy = 2)
9-element Array{Array{Int16,1},1}:
 [139, 129]
 [146, 146]
 [145, 141]
 [126, 126]
 [150, 148]
 [148, 140]
 [185, 199]
 [91, 113]
 [208, 208]
```
"""
function simulate_sample(alleles::Dict{String,NTuple}, loc::Vector{String}; ploidy::Int)
  map(i -> rand(alleles[i], ploidy) ,loc)
end


text:
@testset "Long Sequence" begin
    str = random_dna(50000)
    @test ReferenceSequence(str) == LongDNASeq(str)
    str = ("N"^300 * random_dna(1000))^5 * "N"^300
    @test ReferenceSequence(str) == LongDNASeq(str)
end

ext: jl


text:
function double_conv(in_channels, out_channels; activation=relu, bn_momentum=0.1f0)
    Chain(Conv((3,3), in_channels=>out_channels,  activation; pad=1),
          BatchNorm(out_channels; momentum=bn_momentum),
          Conv((3,3), out_channels=>out_channels, activation; pad=1),
          BatchNorm(out_channels; momentum=bn_momentum)
         )
end

struct UNetUpBlock{U,C}
    upsampling_op::U
    conv_op::C
end

Flux.@functor UNetUpBlock

function UNetUpBlock(in_ch_up, in_ch_concat, out_ch; activation=relu, bn_momentum=0.1f0)
    up_op = Upsample(:bilinear; scale=2)
    conv_op = double_conv(in_ch_up + in_ch_concat, out_ch; activation=activation, bn_momentum=bn_momentum)
    UNetUpBlock(up_op, conv_op)
end

function Base.show(io::IO, b::UNetUpBlock)
    println(io, "UpBlock")
    println(io, b.upsampling_op)
    println(io, b.conv_op)
end

function (b::UNetUpBlock)(up_input, concat_input; dims=3)
    up = cat(b.upsampling_op(up_input), concat_input; dims=dims)
    return b.conv_op(up)
end

"""
    UNet(in_channels::Integer=3, num_classes::Integer=1; init_channels::Integer=16, stages::Integer=4, final_activation=sigmoid)
    UNet(m::Classifier; num_classes=1, decoder_channels=(16,32,64,128,256,512,1024), final_activation=sigmoid)

There are two options to instantiate a [UNet](https://arxiv.org/pdf/1505.04597.pdf):

1. with a given number of input channels and output classes
2. with a classifier from Metalhead, either ResNet or VGG

In the first case you get a UNet with a simple backbone using double convs, where you can specify the number of input channels.
In the second case, you can only specify the number of classes.
"""
struct UNet{E,D,S}
    encoder::E
    decoder::D
    segmentation_head::S
end

Flux.@functor UNet

function UNet(in_channels::Integer=3, num_classes::Integer=1;
              init_channels::Integer=16,
              stages::Integer=4,
              final_activation=sigmoid)

    down_stage1 = double_conv(in_channels, init_channels)
    down_stages = ntuple(i -> Chain(MaxPool((2,2)),
                                    double_conv(init_channels*2^(i-1), init_channels*2^i)),
                        stages)

    encoder = Chain(down_stage1, down_stages...)
    decoder = ntuple(i -> UNetUpBlock(init_channels*2^(stages-i+1),
                                      init_channels*2^(stages-i),
                                      init_channels*2^(stages-i)),
                    stages)

    segmentation_head = Conv((1,1), init_channels=>num_classes, final_activation)

    UNet(encoder, decoder, segmentation_head)
end

function UNet(m::Classifier; num_classes=1, decoder_channels=(16,32,64,128,256,512,1024), final_activation=sigmoid)
    enc = encoder(m)
    enc_channels = encoder_channels(enc)
    decoder_channels = decoder_channels[1:length(enc_channels)-1]
    decoder_channels = (decoder_channels..., last(enc_channels))

    decoder = ntuple(i -> UNetUpBlock(decoder_channels[i+1], enc_channels[i], decoder_channels[i]), length(enc_channels)-1)[end:-1:1]

    segmentation_head = Conv((1,1), decoder_channels[1]=>num_classes, final_activation)
    UNet(enc, decoder, segmentation_head)
end

function decode(ops::Tuple, ft::Tuple)
    up = first(ops)(ft[end], ft[end-1])
    decode(Base.tail(ops), (ft[1:end-2]..., up))
end

decode(::Tuple{}, ft::NTuple{1, T}) where T = first(ft)

function (u::UNet)(input)
    encoder_features = Flux.activations(u.encoder, input)
    up = decode(u.decoder, encoder_features)
    return u.segmentation_head(up)
end

function Base.show(io::IO, u::UNet)
    println(io, "UNet:")
    print(io, "\n")
    println(io, "Encoder:")
    Flux._big_show(io, u.encoder)

    println(io, "\n")
    println(io, "Decoder:")
    for l in u.decoder
        println(io, l)
    end
    println(io, "Segmentation head:")
    println(io, u.segmentation_head)
end


text:
using BenchmarkTools
using FillArrays
using LinearAlgebra: triu, tril

# Subtract the overhead from benchmark times
BenchmarkTools.DEFAULT_PARAMETERS.overhead = BenchmarkTools.estimate_overhead()

const SUITE = BenchmarkGroup()

###
### Eye
###

g = addgroup!(SUITE, "Eye", [])

eye1float = Eye{Float64}(1)
eye10float = Eye{Float64}(10)
eye1000float = Eye{Float64}(1000)
eye10int = Eye{Int}(10)
eye1000int = Eye{Int}(1000)

r1 = addgroup!(g, "reduction", [])
r2 = addgroup!(g, "properties", ["properties"])
r3 = addgroup!(g, "any/all", [])
r4 = addgroup!(g, "iterate", [])
r5 = addgroup!(g, "identities", [])

dimstring(a) = string("n=", size(a, 1))
fulldimstring(a) = string("size=", size(a))
funop(fun, op) = string(fun,"(", op, ", a)")

for a in (eye10float, eye1000float, eye10int, eye1000int)
    for fun in (sum,)
        r1[string(fun), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($a)
    end
    for fun in (isone, iszero)
        r2[string(fun), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($a)
    end
    for (fun, op) in ((any, isone), (all, isone))
        r3[funop(fun, op), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($op, $a)
    end
    for fun in (collect,)
        r4[string(fun), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($a)
    end
    for fun in (permutedims, triu, tril, inv)
        r5[string(fun), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($a)
    end
end

for a in (eye1float,)
    for (fun, op) in ((any, isone), (all, isone))
        r3[funop(fun, op), string(eltype(a)), dimstring(a)] = @benchmarkable $fun($op, $a)
    end
end

###
### Zeros
###

g1 = addgroup!(SUITE, "Zeros", [])

zeros10 = Zeros(10)
zeros10x10 = Zeros(10, 10)

z1 = addgroup!(g1, "properties", ["properties"])

for a in (zeros10, zeros10x10)
    for fun in (iszero, )
        z1[string(fun), string(eltype(a)), fulldimstring(a)] = @benchmarkable $fun($a)
    end
end

# If a cache of tuned parameters already exists, use it, otherwise, tune and cache
# the benchmark parameters. Reusing cached parameters is faster and more reliable
# than re-tuning `SUITE` every time the file is included.
paramspath = joinpath(dirname(@__FILE__), "params.json")

if isfile(paramspath)
    loadparams!(SUITE, BenchmarkTools.load(paramspath)[1], :evals);
else
    tune!(SUITE)
    BenchmarkTools.save(paramspath, params(SUITE));
end

result = run(SUITE, verbose = true)


text:
eval(quote
using .Javis: Luxor, background, sethue, Video, Background, Object, render

export ground,
       luxify,
       nfolds

# default canvas background
function ground(args...)
    background("white") # canvas background
    sethue("black") # pen color
end

# convert a lazysets to a polygonal overapproximation with the given tolerance,
# and then to the function `poly` from Luxor.jl that receives array of points
function luxify(X::LazySet; tol=1e-3)

    # ε-close outer polygonal approximation of X
    P = overapproximate(X, HPolygon, tol)

    # convert to vertex representation
    V = convert(VPolygon, P)
    out = vertices_list(P)

    out_luxor = [Luxor.Point(Tuple(p)) for p in out]

    #return Luxor.poly(out_luxor, :stroke, close=true)
    return out_luxor
end

# useful function to split the array x into n equal parts (if possible)
# if the `rev` flag to return the partition in reverse order
# this function is useful to split in equal parts the animation time for an array of sets
function nfolds(x::AbstractVector, n::Int; rev=true)
    p = length(x)
    m = min(n, p)
    s = p / m
    out = [x[round(Int64, (i-1)*s)+1:min(p,round(Int64, i*s))] for i in 1:m]
    return rev ? reverse(out) : out
end

# animate an array of sets using Javis.jl
# TODO pass colors, eg. with default cols = Javis.distinguishable_colors(length(X))
function Javis.animate(X::AbstractVector{ST}; kwargs...) where {N, ST<:LazySet{N}}

    # size of the animation
    sz = get(kwargs, :size, (500, 500))

    # number of sets to render
    nsets = length(X)

    # maximum time of the animation
    tmax = get(kwargs, :tmax, 10*nsets)

    # number of time divisions (pass either an integer or a partition)
    ntdiv = get(kwargs, :ntdiv, nsets)
    tdiv = nfolds(1:tmax, ntdiv)

    # convert to Luxor.jl object
    Xlux = luxify.(X)

    # pick scaling factor
    α = get(kwargs, :scale, 50.0) # FIXME what is a good default? use radius of box_approximation(X[1]) ?

    # create empty video
    vid = Video(sz...)
    Background(1:tmax, ground)

    # add objects to the canvas
    for (dt, p) in zip(tdiv, Xlux)
        Object(dt, (args...) -> Javis.poly(α .* p, :stroke, close=true))
    end

    # render video
    render(vid)
end

end)


text:
module Allocs

using Base.StackTraces: StackTrace, StackFrame, lookup
using Base: InterpreterIP

# --- Raw results structs, originally defined in C ---

# The C jl_bt_element_t object contains either an IP pointer (size_t) or a void*.
const BTElement = Csize_t;

# matches jl_raw_backtrace_t on the C side
struct RawBacktrace
    data::Ptr{BTElement} # in C: *jl_bt_element_t
    size::Csize_t
end

# matches jl_raw_alloc_t on the C side
struct RawAlloc
    type::Ptr{Type}
    backtrace::RawBacktrace
    size::Csize_t
    task::Ptr{Cvoid}
    timestamp::UInt64
end

# matches jl_profile_allocs_raw_results_t on the C side
struct RawResults
    allocs::Ptr{RawAlloc}
    num_allocs::Csize_t
end

"""
    Profile.Allocs.@profile [sample_rate=0.0001] expr

Profile allocations that happen during `expr`, returning
both the result and and AllocResults struct.

A sample rate of 1.0 will record everything; 0.0 will record nothing.

```julia
julia> Profile.Allocs.@profile sample_rate=0.01 peakflops()
1.03733270279065e11

julia> results = Profile.Allocs.fetch()

julia> last(sort(results.allocs, by=x->x.size))
Profile.Allocs.Alloc(Vector{Any}, Base.StackTraces.StackFrame[_new_array_ at array.c:127, ...], 5576)
```

!!! note
    The current implementation of the Allocations Profiler does not
    capture types for all allocations. Allocations for which the profiler
    could not capture the type are represented as having type
    `Profile.Allocs.UnknownType`.

    You can read more about the missing types and the plan to improve this, here:
    https://github.com/JuliaLang/julia/issues/43688.

!!! compat "Julia 1.8"
    The allocation profiler was added in Julia 1.8.
"""
macro profile(opts, ex)
    _prof_expr(ex, opts)
end
macro profile(ex)
    _prof_expr(ex, :(sample_rate=0.0001))
end

function _prof_expr(expr, opts)
    quote
        $start(; $(esc(opts)))
        try
            $(esc(expr))
        finally
            $stop()
        end
    end
end

"""
    Profile.Allocs.start(sample_rate::Real)

Begin recording allocations with the given sample rate
A sample rate of 1.0 will record everything; 0.0 will record nothing.
"""
function start(; sample_rate::Real)
    ccall(:jl_start_alloc_profile, Cvoid, (Cdouble,), Float64(sample_rate))
end

"""
    Profile.Allocs.stop()

Stop recording allocations.
"""
function stop()
    ccall(:jl_stop_alloc_profile, Cvoid, ())
end

"""
    Profile.Allocs.clear()

Clear all previously profiled allocation information from memory.
"""
function clear()
    ccall(:jl_free_alloc_profile, Cvoid, ())
    return nothing
end

"""
    Profile.Allocs.fetch()

Retrieve the recorded allocations, and decode them into Julia
objects which can be analyzed.
"""
function fetch()
    raw_results = ccall(:jl_fetch_alloc_profile, RawResults, ())
    return decode(raw_results)
end

# decoded results

struct Alloc
    type::Any
    stacktrace::StackTrace
    size::Int
    task::Ptr{Cvoid}
    timestamp::UInt64
end

struct AllocResults
    allocs::Vector{Alloc}
end

# Without this, the Alloc's stacktrace prints for lines and lines and lines...
function Base.show(io::IO, a::Alloc)
    stacktrace_sample = length(a.stacktrace) >= 1 ? "$(a.stacktrace[1]), ..." : ""
    print(io, "$Alloc($(a.type), $StackFrame[$stacktrace_sample], $(a.size))")
end

const BacktraceCache = Dict{BTElement,Vector{StackFrame}}

# copied from julia_internal.h
const JL_BUFF_TAG = UInt(0x4eadc000)
const JL_GC_UNKNOWN_TYPE_TAG = UInt(0xdeadaa03)

struct CorruptType end
struct BufferType end
struct UnknownType end

function load_type(ptr::Ptr{Type})
    if UInt(ptr) < UInt(4096)
        return CorruptType
    elseif UInt(ptr) == JL_BUFF_TAG
        return BufferType
    elseif UInt(ptr) == JL_GC_UNKNOWN_TYPE_TAG
        return UnknownType
    end
    return unsafe_pointer_to_objref(ptr)
end

function decode_alloc(cache::BacktraceCache, raw_alloc::RawAlloc)::Alloc
    Alloc(
        load_type(raw_alloc.type),
        stacktrace_memoized(cache, load_backtrace(raw_alloc.backtrace)),
        UInt(raw_alloc.size),
        raw_alloc.task,
        raw_alloc.timestamp
    )
end

function decode(raw_results::RawResults)::AllocResults
    cache = BacktraceCache()
    allocs = [
        decode_alloc(cache, unsafe_load(raw_results.allocs, i))
        for i in 1:raw_results.num_allocs
    ]
    return AllocResults(allocs)
end

function load_backtrace(trace::RawBacktrace)::Vector{BTElement}
    out = Vector{BTElement}()
    for i in 1:trace.size
        push!(out, unsafe_load(trace.data, i))
    end

    return out
end

function stacktrace_memoized(
    cache::BacktraceCache,
    trace::Vector{BTElement},
    c_funcs::Bool=true
)::StackTrace
    stack = StackTrace()
    for ip in trace
        frames = get(cache, ip) do
            res = lookup(ip)
            cache[ip] = res
            return res
        end
        for frame in frames
            # Skip frames that come from C calls.
            if c_funcs || !frame.from_c
                push!(stack, frame)
            end
        end
    end
    return stack
end

# Precompile once for the package cache.
@assert precompile(start, ())
@assert precompile(stop, ())
@assert precompile(fetch, ())

end

max_stars_repo_name: GarnierAdrian/julia

max_issues_repo_name: GarnierAdrian/julia

max_forks_repo_name: GarnierAdrian/julia


text:
using DiffEqFlux, OrdinaryDiffEq, DiffEqSensitivity
using CUDA, Test, Zygote
CUDA.allowscalar(false)

H = CuArray(rand(Float32, 2, 2))
ann = FastChain(FastDense(1, 4, tanh))
p = initial_params(ann)

function func(x, p, t)
    ann([t],p)[1]*H*x
end

x0 = CuArray(rand(Float32, 2))
x1 = CuArray(rand(Float32, 2))

prob = ODEProblem(func, x0, (0.0f0, 1.0f0))

function evolve(p)
    solve(prob, Tsit5(), p=p, save_start=false,
          save_everystep=false, abstol=1e-4, reltol=1e-4,
          sensealg=QuadratureAdjoint(autojacvec=ZygoteVJP())).u[1]
end

function cost(p)
    x = evolve(p)
    c = sum(abs,x - x1)
    #println(c)
    c
end

grad = Zygote.gradient(cost,p)[1]
@test !iszero(grad[1])
@test iszero(grad[2:4])
@test !iszero(grad[5])
@test iszero(grad[6:end])


text:
# This is a Julia version of Solution of the optimal control problem
# based on code written by Andrea Walther. See:
# Walther, Andrea, and Narayanan, Sri Hari Krishna. Extending the Binomial Checkpointing
# Technique for Resilience. United States: N. p., 2016. https://www.osti.gov/biblio/1364654.

using Checkpointing
using Zygote

include("optcontrolfunc.jl")

function header()
        println("**************************************************************************")
        println("*              Solution of the optimal control problem                   *")
        println("*                                                                        *")
        println("*                     J(y) = y_2(1) -> min                               *")
        println("*           s.t.   dy_1/dt = 0.5*y_1(t) + u(t),            y_1(0)=1      *")
        println("*                  dy_2/dt = y_1(t)^2 + 0.5*u(t)^2         y_2(0)=0      *")
        println("*                                                                        *")
        println("*                  the adjoints equations fulfill                        *")
        println("*                                                                        *")
        println("*         dl_1/dt = -0.5*l_1(t) - 2*y_1(t)*l_2(t)          l_1(1)=0      *")
        println("*         dl_2/dt = 0                                      l_2(1)=1      *")
        println("*                                                                        *")
        println("*   with Revolve for Online and (Multi-Stage) Offline Checkpointing      *")
        println("*                                                                        *")
        println("**************************************************************************")

        println("**************************************************************************")
        println("*        The solution of the optimal control problem above is            *")
        println("*                                                                        *")
        println("*        y_1*(t) = (2*e^(3t)+e^3)/(e^(3t/2)*(2+e^3))                     *")
        println("*        y_2*(t) = (2*e^(3t)-e^(6-3t)-2+e^6)/((2+e^3)^2)                 *")
        println("*          u*(t) = (2*e^(3t)-e^3)/(e^(3t/2)*(2+e^3))                     *")
        println("*        l_1*(t) = (2*e^(3-t)-2*e^(2t))/(e^(t/2)*(2+e^3))                *")
        println("*        l_2*(t) = 1                                                     *")
        println("*                                                                        *")
        println("**************************************************************************")

        return
end

function muoptcontrol(scheme, steps)
    println( "\n STEPS    -> number of time steps to perform")
    println("SNAPS    -> number of checkpoints")
    println("INFO = 1 -> calculate only approximate solution")
    println("INFO = 2 -> calculate approximate solution + takeshots")
    println("INFO = 3 -> calculate approximate solution + all information ")
    println(" ENTER:   STEPS, SNAPS, INFO \n")

    # F   : output
    # F_H : input
    # L   : seed the output adjoint
    # L_H : set input adjoint to 0
    F = [1.0, 0.0]
    F_H = [0.0, 0.0]
    t = 0.0
    h = 1.0/steps
    model = Model(F, F_H, t, h)
    # Just make sure it's all zero.
    shadowmodel = Model([0.0,0.0], [0.0,0.0], 0.0, 0.0)

    function foo(model::Model)
        @checkpoint_struct scheme model shadowmodel for i in 1:steps
            model.F_H .= model.F
            advance(model)
            model.t += h
        end
        return model.F[2]
    end
    g = Zygote.gradient(foo,model)

    F = model.F
    L = [g[1].F[1], g[1].F[2]]

    F_opt = Array{Float64, 1}(undef, 2)
    L_opt = Array{Float64, 1}(undef, 2)
    opt_sol(F_opt,1.0)
    opt_lambda(L_opt,0.0)
    println("\n\n")
    println("y_1*(1)  = " , F_opt[1] , " y_2*(1)  = " , F_opt[2])
    println("y_1 (1)  = " , F[1] , "  y_2 (1)  = " , F[2] , " \n\n")
    println("l_1*(0)  = " , L_opt[1] , "  l_2*(0)  = " , L_opt[2])
    println("l_1 (0)  = " , L[1]     , "  sl_2 (0)  = " , L[2] , " ")
    return F, L, F_opt, L_opt
end

# main(10,3,3)


text:
using Plots;

bool(x::Int) = (x!=0); # helper function since Julia doesn't like to use ints as bools
int(x::Float64) = floor(Int, x); # converts floats to int8
int(xs::Array{Int8,1}) = sum([x<<(length(xs)-i) for (i,x) in enumerate(xs)])

plot(grid) = heatmap(1:size(grid,2), 1:size(grid,1),
                reverse(grid, dims=1), c = :greys, axis=([], false),
                legend=false, size=(1920, 1080))

function eca(N::Int, it::Int, rule::Int, k::Int=1)
    mask = 1
    world = zeros(Int8, it, N)  # Create grid
    world[1,int(N/2)] = 1;      # start with a single one in the middle

    # main update loop
    for i in 1:size(world, 1)-1
        for j in 1:size(world, 2)
            b = mod1((j-k+N), N)
            e = mod1((j+k+N), N)
            if b < e
                val = int(world[i,b:e])
            else
                val = int(vcat(world[i, b:N], world[i, 1:e]))
                if bool(rule & (mask<<val))
                end
            end
            if bool(rule & (mask<<val))
                world[i+1, j] = 1
            end
        end
    end
    return world;
end

function k_color(N::Int, it::Int, rule::Int, k::Int=3) # Creates
    world = zeros(Int8, it, N)  # Create grid
    world[1,div(N,2)] = 1;      # start with a single one in the middle
    rule_basek = reverse([parse(Int8, c) for c in string(rule, base=k)])

    # main update loop
    for i in 1:size(world, 1)-1
        for j in 1:size(world, 2)
            b = mod1((j-1+N), N)
            e = mod1((j+1+N), N)
            if b < e
                val = sum(world[i,b:e])
            else
                val = sum([world[i,b], world[i,j], world[i,e]])
            end
            world[i+1, j] = rule_basek[val+1]
        end
    end
    return world;
end

max_stars_repo_name: Stevellen/cellular-automata

max_issues_repo_name: Stevellen/cellular-automata

max_forks_repo_name: Stevellen/cellular-automata


text:
#  Copyright 2017, Iain Dunning, Joey Huchette, Miles Lubin, and contributors
#  This Source Code Form is subject to the terms of the Mozilla Public
#  License, v. 2.0. If a copy of the MPL was not distributed with this
#  file, You can obtain one at https://mozilla.org/MPL/2.0/.
#############################################################################
# JuMP
# An algebraic modeling language for Julia
# See https://github.com/jump-dev/JuMP.jl
#############################################################################

# This file contains an example bridge used for tests.

using JuMP

const MOIB = MOI.Bridges
const MOIBC = MOI.Bridges.Constraint

"""
    Nonnegative <: MOI.AbstractScalarSet

Scalar set of nonnegative numbers.
"""
struct Nonnegative <: MOI.AbstractScalarSet end

MOI.copy(set::Nonnegative) = set

"""
    NonnegativeBridge{T}

The `NonnegativeBridge` replaces a constraint `func`-in-`Nonnegative` into
`func`-in-`GreaterThan{T}`.
"""
struct NonnegativeBridge{T,F<:MOI.AbstractScalarFunction} <:
       MOIBC.AbstractBridge
    constraint_index::MOI.ConstraintIndex{F,MOI.GreaterThan{T}}
end

function MOIBC.bridge_constraint(
    ::Type{NonnegativeBridge{T,F}},
    model,
    f::F,
    s::Nonnegative,
) where {T,F}
    ci = MOIU.normalize_and_add_constraint(model, f, MOI.GreaterThan(zero(T)))
    return NonnegativeBridge{T,F}(ci)
end

function MOI.supports_constraint(
    ::Type{NonnegativeBridge{T}},
    ::Type{<:MOI.AbstractScalarFunction},
    ::Type{Nonnegative},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{<:NonnegativeBridge})
    return Tuple{DataType}[]
end
function MOIB.added_constraint_types(::Type{NonnegativeBridge{T,F}}) where {T,F}
    return [(F, MOI.GreaterThan{T})]
end
function MOIBC.concrete_bridge_type(
    ::Type{NonnegativeBridge{T}},
    F::Type{<:MOI.AbstractScalarFunction},
    ::Type{Nonnegative},
) where {T}
    # In the constructor, the function `f` of type `F` is passed to
    # `MOIU.normalize_and_add_constraint` which removes the constrant from `f` but
    # does not change its type so the type of the function in `MOI.GreaterThan`
    # will also be `F`.
    return NonnegativeBridge{T,F}
end

# Attributes, Bridge acting as an model
function MOI.get(
    ::NonnegativeBridge{T,F},
    ::MOI.NumberOfConstraints{F,MOI.GreaterThan{T}},
) where {T,F}
    return 1
end
function MOI.get(
    bridge::NonnegativeBridge{T,F},
    ::MOI.ListOfConstraintIndices{F,MOI.GreaterThan{T}},
) where {T,F}
    return [bridge.constraint_index]
end

# Indices
function MOI.delete(model::MOI.ModelLike, bridge::NonnegativeBridge)
    return MOI.delete(model, bridge.constraint_index)
end

# Attributes, Bridge acting as a constraint
function MOI.get(
    model::MOI.ModelLike,
    attr::Union{MOI.ConstraintPrimal,MOI.ConstraintDual},
    bridge::NonnegativeBridge,
)
    return MOI.get(model, attr, bridge.constraint_index)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
function errors_main()
    lim = 10

    uf = UnionFinder(lim)

    @test_throws BoundsError union!(uf, lim, 0)
    @test_throws BoundsError union!(uf, 0, lim)
    @test_throws BoundsError union!(uf, lim, lim + 1)
    @test_throws BoundsError union!(uf, lim + 1, lim)

    @test_throws ArgumentError union!(uf, [lim, lim, lim], [lim, lim])

    @test_throws BoundsError union!(uf, [lim], [0])
    @test_throws BoundsError union!(uf, [0], [lim])
    @test_throws BoundsError union!(uf, [lim], [lim + 1])
    @test_throws BoundsError union!(uf, [lim + 1], [lim])

    @test_throws BoundsError union!(uf, [(lim, 0)])
    @test_throws BoundsError union!(uf, [(0, lim)])
    @test_throws BoundsError union!(uf, [(lim, lim + 1)])
    @test_throws BoundsError union!(uf, [(lim + 1, lim)])

    @test_throws BoundsError find!(uf, lim + 1)
    @test_throws BoundsError find!(uf, 0)

    @test_throws BoundsError size!(uf, lim + 1)
    @test_throws BoundsError size!(uf, 0)

    cf = CompressedFinder(uf)

    @test_throws BoundsError find(cf, lim + 1)
    @test_throws BoundsError find(cf, 0)

    @test_throws ArgumentError UnionFinder(0)
end

errors_main()


text:
let stages = (POST_EPISODE_STAGE, PRE_EPISODE_STAGE)
    @testset "DoEveryNEpisode stage=$(stage),s2=$(s2)" for stage in stages, s2 in stages
        hook = DoEveryNEpisode((x...) -> true; stage)
        if stage === s2
            @test hook(stage, nothing, nothing)
        else
            @test hook(s2, nothing, nothing) === nothing
        end
    end
end


text:
## Exercise 3-2
## Write a function named rightjustify that takes a string named s as a parameter and prints the string with enough leading spaces so that the last letter of the string is in column 70 of the display.
println("Ans.")

function rightjustify(string::String)
    return ' ' ^ (70 - length(string)) * string
end

println(rightjustify("Monty"))

println("End.")


text:
export linearInter

function linearInter(T::Array{Float64},domain::Array{Float64},xin::Array{Float64};doDerivative::Bool=false)

	dim = round(Int,length(domain)/2)
	n   = (dim==1) ? length(T) : size(T)
	h   = sdiag(1./collect(n)) * vec(domain[2:2:end]-domain[1:2:end])
	np  = round(Int,length(xin)/dim)
	x   = reshape(copy(xin),np,dim)

	# map x from [h/2,domain-h/2] -> [1,m],
	for i=1:dim
		x[:,i] = (x[:,i]-domain[2*i-1])/h[i] + 0.5;
	end
	Tc = zeros(np)                                # initialize output
	Valid(j) = (0.<x[:,j]) .& (x[:,j].<(n[j]+1))      # determine indices of valid points

	if dim==1
		valid = find(Valid(1))
	elseif dim==2
		valid = find( Valid(1) .& Valid(2) );
	elseif dim==3
		valid = find( Valid(1) .& Valid(2) .& Valid(3) );
	end

	if isempty(valid)
		if !doDerivative
			return Tc
		else
			return Tc, spzeros(np,dim*np)
		end
	end

	pad = 1; TP = zeros(tuple(collect(n)+2*pad...));                 # pad data to reduce cases
	P     = round.(Int,floor.(x)); x = x-P;                        # split x into integer/remainder
	p(j)  = P[valid,j]
	xi(j) = x[valid,j]

	# increments for linearized ordering
	i1 = 1; i2 = size(T,1)+2*pad; i3 = (size(T,1)+2*pad)*(size(T,2)+2*pad);

	if dim==1
		TP[pad+(1:n)] = vec(T)
		p = pad + p(1);
	    Tc[valid] = TP[p].* (1-xi(1)) + TP[p+1].*xi(1);   # compute weighted sum

	    if doDerivative
			dT = zeros(np)
			dT[valid] = TP[p+1]-TP[p]
			dT = sdiag(dT[:,1]/h[1])
		end;
	elseif dim==2
	    TP[pad+(1:n[1]),pad+(1:n[2])] = T
	    TP = vec(TP)
	    p  = (pad + p(1)) + i2*(pad + p(2) - 1)

	    # compute Tc as weighted sum
	    Tc[valid] = (TP[p]   .* (1-xi(1))  + TP[p+i1]    .*xi(1)) .* (1-xi(2))  + (TP[p+i2] .* (1-xi(1)) + TP[p+i1+i2] .*xi(1)) .* xi(2)

		if doDerivative
			dT = zeros(np,dim)
			dT[valid,1] = (TP[p+i1]-TP[p]).*(1-xi(2)) + (TP[p+i1+i2]-TP[p+i2]).*xi(2)
			dT[valid,2] = (TP[p+i2]-TP[p]).*(1-xi(1)) + (TP[p+i1+i2]-TP[p+i1]).*xi(1)
			dT = spdiagm((dT[:,1]/h[1], dT[:,2]/h[2]),(0,np), np, 2*np)
			# dT = [sdiag(dT[:,1]/h[1])  sdiag(dT[:,2]/h[2])]
		end
	  elseif dim==3
	    TP[pad+(1:n[1]),pad+(1:n[2]),pad+(1:n[3])] = T
	    p  = (pad + p(1)) + i2*(pad + p(2) - 1) + i3*(pad + p(3) -1);
	    # compute Tc as weighted sum
	    Tc[valid] = ((TP[p].*(1-xi(1))+TP[p+i1].*xi(1)).*(1-xi(2)) +
		  (TP[p+i2].*(1-xi(1))+TP[p+i1+i2].*xi(1)).*(xi(2))).*(1-xi(3)) +
	      +((TP[p+i3].*(1-xi(1))+TP[p+i1+i3].*xi(1)).*(1-xi(2)) +
	      +(TP[p+i2+i3].*(1-xi(1))+TP[p+i1+i2+i3].*xi(1)).*(xi(2))).*(xi(3))

	    if doDerivative
			dT = zeros(np,dim)
			dT[valid,1] = ((TP[p+i1]-TP[p]).*(1-xi(2))+(TP[p+i1+i2]-TP[p+i2]).*xi(2)).*(1-xi(3)) +
	                      ((TP[p+i1+i3]-TP[p+i3]).*(1-xi(2))+(TP[p+i1+i2+i3]-TP[p+i2+i3]).*xi(2)).*(xi(3))
		    dT[valid,2] = ((TP[p+i2]-TP[p]).*(1-xi(1))+(TP[p+i1+i2]-TP[p+i1]).*xi(1)).*(1-xi(3)) +
	                      ((TP[p+i2+i3]-TP[p+i3]).*(1-xi(1))+(TP[p+i1+i2+i3]-TP[p+i1+i3]).*xi(1)).*(xi(3))
		    dT[valid,3] = ((TP[p+i3].*(1-xi(1))+TP[p+i1+i3].*xi(1)).*(1-xi(2)) +
	                      (TP[p+i2+i3].*(1-xi(1))+TP[p+i1+i2+i3].*xi(1)).*(xi(2))) -
	                      ((TP[p].*(1-xi(1))+TP[p+i1].*xi(1)).*(1-xi(2)) +
	                      (TP[p+i2].*(1-xi(1))+TP[p+i1+i2]*xi(1)).*(xi(2)));
			dT = [sdiag(dT[:,1]/h[1])  sdiag(dT[:,2]/h[2])  sdiag(dT[:,3]/h[3])]
		end
	end
	if doDerivative
		return Tc,dT
	else
		return Tc
	end
end


text:
"""
    hmul!(C::HMatrix,A::HMatrix,B::HMatrix,a,b,compressor)

Similar to `mul!` : compute `C <-- A*B*a + B*b`, where `A,B,C` are
hierarchical matrices and `compressor` is a function/functor used in the
intermediate stages of the multiplication to avoid growring the rank of
admissible blocks after addition is performed.
"""
function hmul!(C::HMatrix,A::HMatrix,B::HMatrix,a,b,compressor)
    b == true || rmul!(C,b)
    @timeit_debug "constructing plan" begin
        plan = plan_hmul(C,A,B,a,1)
    end
    @timeit_debug "executing plan" begin
        execute!(plan,compressor)
    end
    return C
end

# disable `mul` of hierarchial matrices
function mul!(C::HMatrix,A::HMatrix,B::HMatrix,a::Number,b::Number)
    msg = "use `hmul` to multiply hierarchical matrices"
    error(msg)
end

"""
    struct HMulNode{S,T} <: AbstractTree

Tree data structure representing the following computation:
```
    C <-- C + a * ∑ᵢ Aᵢ * Bᵢ
```
where `C = target(node)`, and `Aᵢ,Bᵢ` are pairs stored in `sources(node)`.

This structure is used to group the operations required when multiplying
hierarchical matrices so that they can later be executed in a way that minimizes
recompression of intermediate computations.
"""
mutable struct HMulNode{T} <: AbstractTree
    target::T
    children::Matrix{HMulNode{T}}
    sources::Vector{Tuple{T,T}}
    multiplier::Float64
end

function HMulNode(C::HMatrix,a::Number)
    T    = typeof(C)
    chdC = children(C)
    m,n  = size(chdC)
    HMulNode{T}(C,Matrix{HMulNode{T}}(undef,m,n),T[],a)
end

function build_HMulNode_structure(C::HMatrix,a)
    node = HMulNode(C,a)
    chdC = children(C)
    m,n  = size(chdC)
    for i in 1:m
        for j in 1:n
            child = build_HMulNode_structure(chdC[i,j],a)
            node.children[i,j] = child
        end
    end
    node
end

function plan_hmul(C::T,A::T,B::T,a,b) where {T<:HMatrix}
    @assert b == 1
    # root = HMulNode(C)
    root = build_HMulNode_structure(C,a)
    # recurse
    _build_hmul_tree!(root,A,B)
    return root
end

function _build_hmul_tree!(tree::HMulNode,A::HMatrix,B::HMatrix)
    C = tree.target
    if isleaf(A) || isleaf(B) || isleaf(C)
        push!(tree.sources,(A,B))
    else
        ni,nj = blocksize(C)
        _ ,nk = blocksize(A)
        A_children = children(A)
        B_children = children(B)
        C_children = children(C)
        for i=1:ni
            for j=1:nj
                child = tree.children[i,j]
                for k=1:nk
                    _build_hmul_tree!(child,A_children[i,k],B_children[k,j])
                end
            end
        end
    end
    return tree
end

function Base.show(io::IO,::MIME"text/plain",tree::HMulNode)
    print(io,"HMulNode with $(size(children(tree))) children and $(length(sources(tree))) pairs")
end

function Base.show(io::IO,tree::HMulNode)
    print(io,"HMulNode with $(size(children(tree))) children and $(length(sources(tree))) pairs")
end

function Base.show(io::IO,::MIME"text/plain",tree::Adjoint{<:Any,<:HMulNode})
    p = parent(tree)
    print(io,"adjoint HMulNode with $(size(children(p))) children and $(length(sources(p))) pairs")
end

function Base.show(io::IO,tree::Adjoint{<:Any,<:HMulNode})
    p = parent(tree)
    print(io,"adjoint HMulNode with $(size(children(p))) children and $(length(sources(p))) pairs")
end

# compress the operator L = C + ∑ a*Aᵢ*Bᵢ
function (paca::PartialACA)(plan::HMulNode)
    _aca_partial(plan,:,:,paca.atol,paca.rank,paca.rtol)
end

# getters
target(node::HMulNode) = node.target
sources(node::HMulNode) = node.sources
multiplier(node::HMulNode) = node.multiplier

# Trees interface
Trees.children(node::HMulNode) = node.children
Trees.children(node::HMulNode,idxs...) = node.children[idxs]
Trees.parent(node::HMulNode)   = node.parent
Trees.isleaf(node::HMulNode)   = isempty(children(node))
Trees.isroot(node::HMulNode)   = parent(node) === node

# AbstractMatrix interface
Base.size(node::HMulNode) = size(target(node))
Base.eltype(node::HMulNode) = eltype(target(node))

Base.getindex(node::HMulNode,::Colon,j::Int) = getcol(node,j)

function getcol(node::HMulNode,j)
    m,n = size(node)
    T   = eltype(node)
    col = zeros(T,m)
    getcol!(col,node,j)
    return col
end

function getcol!(col,node::HMulNode,j)
    a = multiplier(node)
    C = target(node)
    m,n = size(C)
    T  = eltype(C)
    ej = zeros(T,n)
    ej[j] = 1
    # compute j-th column of ∑ Aᵢ Bᵢ
    for (A,B) in sources(node)
        m,k = size(A)
        k,n = size(B)
        tmp = zeros(T,k)
        jg   = j + offset(B)[2] # global index on hierarchila matrix B
        getcol!(tmp,B,jg)
        _hgemv_recursive!(col,A,tmp,offset(A))
    end
    # multiply by a
    rmul!(col,a)
    # add the j-th column of C if C has data
    # jg  = j + offset(C)[2] # global index on hierarchila matrix B
    # cj  = getcol(C,jg)
    # axpy!(1,cj,col)
    if hasdata(C)
        d  = data(C)
        cj = getcol(d,j)
        axpy!(1,cj,col)
    end
    return col
end

adjoint(node::HMulNode) = Adjoint(node)
Base.size(adjnode::Adjoint{<:Any,<:HMulNode}) = reverse(size(adjnode.parent))
Trees.children(adjnode::Adjoint{<:Any,<:HMulNode}) = adjoint(children(adjnode.parent))

Base.getindex(adjnode::Adjoint{<:Any,<:HMulNode},::Colon,j::Int) = getcol(adjnode,j)

function getcol(adjnode::Adjoint{<:Any,<:HMulNode},j)
    m,n = size(adjnode)
    T   = eltype(adjnode)
    col = zeros(T,m)
    getcol!(col,adjnode,j)
    return col
end

function getcol!(col,adjnode::Adjoint{<:Any,<:HMulNode},j)
    node  = parent(adjnode)
    a     = multiplier(node)
    C     = target(node)
    T     = eltype(C)
    Ct    = adjoint(C)
    m,n   = size(Ct)
    ej    = zeros(T,n)
    ej[j] = 1
    # compute j-th column of ∑ adjoint(Bᵢ)*adjoint(Aᵢ)
    for (A,B) in sources(node)
        At,Bt = adjoint(A), adjoint(B)
        tmp = zeros(T,size(At,1))
        # _hgemv_recursive!(tmp,At,ej,offset(At))
        jg  = j + offset(At)[2] # global index on hierarchila matrix B
        getcol!(tmp,At,jg)
        _hgemv_recursive!(col,Bt,tmp,offset(Bt))
    end
    # multiply by a
    rmul!(col,conj(a))
    # add the j-th column of Ct if it has data
    # jg  = j + offset(Ct)[2] # global index on hierarchila matrix B
    # cj  = getcol(Ct,jg)
    # axpy!(1,cj,col)
    if hasdata(Ct)
        d  = data(Ct)
        cj = getcol(d,j)
        axpy!(1,cj,col)
    end
    return col
end

function execute!(node::HMulNode,compressor)
    execute_node!(node,compressor)
    C = target(node)
    flush_to_children!(C)
    @threads for chd in children(node)
        execute!(chd,compressor)
    end
    return node
end

# non-recursive execution
function execute_node!(node::HMulNode,compressor)
    C = target(node)
    isempty(sources(node)) && (return node)
    a = multiplier(node)
    if isleaf(C) && !isadmissible(C)
        d = data(C)
        for (A,B) in sources(node)
            # _mul_leaf!(C,A,B,a,compressor)
            _mul_dense!(d,A,B,a)
        end
    else
        R = compressor(node)
        setdata!(C,R)
    end
    return node
end

function _mul_dense!(C::Matrix,A,B,a)
    Adata = isleaf(A) ? A.data : A
    Bdata = isleaf(B) ? B.data : B
    if Adata isa HMatrix
        if Bdata isa Matrix
            _mul131!(C, Adata, Bdata, a)
        elseif Bdata isa RkMatrix
            _mul132!(C, Adata, Bdata, a)
        end
    elseif Adata isa Matrix
        if Bdata isa Matrix
            _mul111!(C, Adata, Bdata, a)
        elseif Bdata isa RkMatrix
            _mul112!(C, Adata, Bdata, a)
        elseif Bdata isa HMatrix
            _mul113!(C, Adata, Bdata, a)
        end
    elseif Adata isa RkMatrix
        if Bdata isa Matrix
            _mul121!(C, Adata, Bdata, a)
        elseif Bdata isa RkMatrix
            _mul122!(C, Adata, Bdata, a)
        elseif Bdata isa HMatrix
            _mul123!(C, Adata, Bdata, a)
        end
    end
end

"""
    flush_to_children!(H::HMatrix,compressor)

Transfer the blocks `data` to its children. At the end, set `H.data` to `nothing`.
"""
function flush_to_children!(H::HMatrix)
    T = eltype(H)
    isleaf(H)  && (return H)
    hasdata(H) || (return H)
    R::RkMatrix{T}   = data(H)
    _add_to_children!(H,R)
    setdata!(H,nothing)
    return H
end

function _add_to_children!(H,R::RkMatrix)
    shift = pivot(H) .- 1
    for block in children(H)
        irange   = rowrange(block) .- shift[1]
        jrange   = colrange(block) .- shift[2]
        bdata    = data(block)
        tmp      = RkMatrix(R.A[irange,:],R.B[jrange,:])
        if bdata === nothing
            setdata!(block,tmp)
        else
            axpy!(true,tmp,bdata)
        end
    end
end

"""
    flush_to_leaves!(H::HMatrix,compressor)

Transfer the blocks `data` to its leaves. At the end, set `H.data` to `nothing`.

# See also: [`flush_to_children!`](@ref)
"""
function flush_to_leaves!(H::HMatrix)
    T = eltype(H)
    isleaf(H)  && (return H)
    hasdata(H) || (return H)
    R::RkMatrix{T}   = data(H)
    _add_to_leaves!(H,R)
    setdata!(H,nothing)
    return H
end

function _add_to_leaves!(H,R::RkMatrix)
    shift = pivot(H) .- 1
    for block in Leaves(H)
        irange   = rowrange(block) .- shift[1]
        jrange   = colrange(block) .- shift[2]
        bdata    = data(block)
        tmp      = RkMatrix(R.A[irange,:],R.B[jrange,:])
        if bdata === nothing
            setdata!(block,tmp)
        else
            axpy!(true,tmp,bdata)
        end
    end
end

_mul111!(C::Union{Matrix,SubArray,Adjoint},A::Union{Matrix,SubArray,Adjoint},B::Union{Matrix,SubArray,Adjoint},a::Number) = mul!(C,A,B,a,true)

function _mul112!(C::Union{Matrix,SubArray,Adjoint}, M::Union{Matrix,SubArray,Adjoint}, R::RkMatrix, a::Number)
    buffer = M * R.A
    _mul111!(C, buffer, R.Bt, a)
    return C
end

function _mul113!(C::Union{Matrix,SubArray,Adjoint}, M::Union{Matrix,SubArray,Adjoint}, H::HMatrix, a::Number)
    T = eltype(C)
    if hasdata(H)
        mat = data(H)
        if mat isa Matrix
            _mul111!(C, M, mat, a)
        elseif mat isa RkMatrix
            _mul112!(C, M, mat, a)
        else
            error()
        end
    end
    for child in children(H)
        shift  = pivot(H) .- 1
        irange = rowrange(child) .- shift[1]
        jrange = colrange(child) .- shift[2]
        Cview  = @views C[:, jrange]
        Mview  = @views M[:, irange]
        _mul113!(Cview, Mview, child, a)
    end
    return C
end

function _mul121!(C::Union{Matrix,SubArray,Adjoint}, R::RkMatrix, M::Union{Matrix,SubArray,Adjoint}, a::Number)
    _mul111!(C, R.A, R.Bt * M, a)
end

function _mul122!(C::Union{Matrix,SubArray,Adjoint}, R::RkMatrix, S::RkMatrix, a::Number)
    if rank(R) < rank(S)
        _mul111!(C, R.A, (R.Bt * S.A) * S.Bt, a)
    else
        _mul111!(C, R.A * (R.Bt * S.A), S.Bt, a)
    end
    return C
end

function _mul123!(C::Union{Matrix,SubArray,Adjoint}, R::RkMatrix, H::HMatrix, a::Number)
    T = promote_type(eltype(R), eltype(H))
    tmp = zeros(T, size(R.Bt, 1), size(H, 2))
    _mul113!(tmp, R.Bt, H, 1)
    _mul111!(C, R.A, tmp, a)
    return C
end

function _mul131!(C::Union{Matrix,SubArray,Adjoint}, H::HMatrix, M::Union{Matrix,SubArray,Adjoint}, a::Number)
    if isleaf(H)
        mat = data(H)
        if mat isa Matrix
            _mul111!(C, mat, M, a)
        elseif mat isa RkMatrix
            _mul121!(C, mat, M, a)
        else
            error()
        end
    end
    for child in children(H)
        shift  = pivot(H) .- 1
        irange = rowrange(child) .- shift[1]
        jrange = colrange(child) .- shift[2]
        Cview  = view(C, irange, :)
        Mview  = view(M, jrange, :)
        _mul131!(Cview, child, Mview, a)
    end
    return C
end

function _mul132!(C::Union{Matrix,SubArray,Adjoint}, H::HMatrix, R::RkMatrix, a::Number)
    T = promote_type(eltype(H),eltype(R))
    buffer = zeros(T, size(H, 1), size(R.A, 2))
    _mul131!(buffer,H,R.A,1)
    _mul111!(C, buffer, R.Bt, a,)
    return C
end

############################################################################################
# Specializations on gemv:
# The routines below provide specialized version of mul!(C,A,B,a,b) when `A` and
# `B` are vectors
############################################################################################

# 1.2.1
function mul!(y::AbstractVector,R::RkMatrix,x::AbstractVector,a::Number,b::Number)
    tmp = R.Bt*x
    mul!(y,R.A,tmp,a,b)
end

# 1.2.1
function mul!(y::AbstractVector,adjR::Adjoint{<:Any,<:RkMatrix},x::AbstractVector,a::Number,b::Number)
    R = parent(adjR)
    tmp = R.At*x
    mul!(y,R.B,tmp,a,b)
end

# 1.3.1
"""
    mul!(y::AbstractVector,H::HMatrix,x::AbstractVector,a,b;global_index=true,threads=true)

Perform `y <-- H*x*a + y*b` in place.
"""
function mul!(y::AbstractVector,A::HMatrix,x::AbstractVector,a::Number=1,b::Number=0;
                            global_index=true,threads=true)
    # since the HMatrix represents A = Pr*H*Pc, where Pr and Pc are row and column
    # permutations, we need first to rewrite C <-- b*C + a*(Pc*H*Pb)*B as
    # C <-- Pr*(b*inv(Pr)*C + a*H*(Pc*B)). Following this rewrite, the
    # multiplication is performed by first defining B <-- Pc*B, and C <--
    # inv(Pr)*C, doing the multiplication with the permuted entries, and then
    # permuting the result  C <-- Pr*C at the end.
    ctree     = coltree(A)
    rtree     = rowtree(A)
    # permute input
    if global_index
        x         = x[ctree.loc2glob]
        y         = permute!(y,rtree.loc2glob)
        rmul!(x,a) # multiply in place since this is a new copy, so does not mutate exterior x
    else
        x = a*x # new copy of x
    end
    iszero(b) ? fill!(y,zero(eltype(y))) : rmul!(y,b)
    # offset in case A is not indexed starting at (1,1); e.g. A is not the root
    # of and HMatrix
    offset = pivot(A) .- 1
    if threads
        # TODO: test the various threaded implementations and chose one.
        # Currently there are two main choices:
        # 1. spawn a task per leaf, and let julia scheduler handle the tasks
        # 2. create a static partition of the leaves and try to estimate the
        #    cost, then spawn one task per block of the partition. In this case,
        #    test if the hilbert partition is really faster than col_partition
        #    or row_partition
        #    Right now the hilbert partition is chosen by default without proper
        #    testing.
        @timeit_debug "hilbert partition" begin
            nt        = Threads.nthreads()
            partition = hilbert_partitioning(A,nt,_cost_gemv)
        end
        @timeit_debug "threaded multiplication" begin
            _hgemv_static_partition!(y,x,partition,offset)
        end
        # _hgemv_threads!(y,A,x,offset)  # threaded implementation
    else
        @timeit_debug "serial multiplication" begin
            _hgemv_recursive!(y,A,x,offset) # serial implementation
        end
    end
    # permute output
    global_index && invpermute!(y,loc2glob(rtree))
    return y
end

"""
    _hgemv_recursive!(C,A,B,offset)

Internal function used to compute `C[I] <-- C[I] + A*B[J]` where `I =
rowrange(A) - offset[1]` and `J = rowrange(B) - offset[2]`.

The `offset` argument is used on the caller side to signal if the original
hierarchical matrix had a `pivot` other than `(1,1)`.
"""
function _hgemv_recursive!(C::AbstractVector,A::Union{HMatrix,Adjoint{<:Any,<:HMatrix}},B::AbstractVector,offset)
    T = eltype(A)
    if isleaf(A)
        irange = rowrange(A) .- offset[1]
        jrange = colrange(A) .- offset[2]
        d   = data(A)
        if T <: SMatrix
            # FIXME: there is bug with gemv and static arrays, so we convert
            # them to matrices of n × 1
            # (https://github.com/JuliaArrays/StaticArrays.jl/issues/966#issuecomment-943679214).
            mul!(view(C,irange,1:1),d,view(B,jrange,1:1),1,1)
        else
            # C and B are the "global" vectors handled by the caller, so a view
            # is needed.
            mul!(view(C,irange),d,view(B,jrange),1,1)
        end
    else
        for block in children(A)
            _hgemv_recursive!(C,block,B,offset)
        end
    end
    return C
end

function _hgemv_threads!(C::AbstractVector,A::HMatrix,B::AbstractVector,offset)
    nt        = Threads.nthreads()
    # make `nt` copies of C and run in parallel. The tree is partitioned up to a
    # given granularity to avoid creating too many small tasks.
    Cthreads  = [zero(C) for _ in 1:nt]
    blocks = filter_tree(A,true) do x
        (isleaf(x) || length(x)<1000*1000)
    end
    sort!(blocks;lt=(x,y)->length(x)<length(y),rev=true)
    n = length(blocks)
    @sync for i in 1:n
        block = blocks[i]
        Threads.@spawn begin
            id = Threads.threadid()
            _hgemv_recursive!(Cthreads[id],block,B,offset)
        end
    end
    # reduce
    for Ct in Cthreads
        axpy!(1,Ct,C)
    end
    return C
end

function _hgemv_static_partition!(C::AbstractVector,B::AbstractVector,partition,offset)
    # create a lock for the reduction step
    T = eltype(C)
    mutex = ReentrantLock()
    nt    = length(partition)
    times = zeros(nt)
    Threads.@threads for n in 1:nt
        id = Threads.threadid()
        times[id] =
        @elapsed begin
            leaves = partition[n]
            Cloc   = zero(C)
            for leaf in leaves
                irange = rowrange(leaf) .- offset[1]
                jrange = colrange(leaf) .- offset[2]
                data   = leaf.data
                if T <: SVector
                    mul!(view(Cloc,irange,1:1),data,view(B,jrange,1:1),1,1)
                else
                    mul!(view(Cloc,irange),data,view(B,jrange),1,1)
                end
            end
            # reduction
            lock(mutex) do
                axpy!(1,Cloc,C)
            end
        end
    end
    tmin,tmax = extrema(times)
    if tmax/tmin > 1.1
        @debug "gemv: ratio of tmax/tmin = $(tmax/tmin)"
    end
    return C
end

"""
    hilbert_partitioning(H::HMatrix,np,cost)

Partiotion the leaves of `H` into `np` sequences of approximate equal cost (as
determined by the `cost` function) while also trying to maximize the locality of
each partition.
"""
function hilbert_partitioning(H::HMatrix,np,cost)
    # the hilbert curve will be indexed from (0,0) × (N-1,N-1), so set N to be
    # the smallest power of two larger than max(m,n), where m,n = size(H)
    m,n = size(H)
    N   = max(m,n)
    N   = nextpow(2,N)
    # sort the leaves by their hilbert index
    leaves = Leaves(H) |> collect
    hilbert_indices = map(leaves) do leaf
        # use the center of the leaf as a cartesian index
        i,j = pivot(leaf) .- 1 .+ size(leaf) .÷ 2
        hilbert_cartesian_to_linear(N,i,j)
    end
    p = sortperm(hilbert_indices)
    permute!(leaves,p)
    # now compute a quasi-optimal partition of leaves based `cost_mv`
    cmax      = find_optimal_cost(leaves,np,cost,1)
    partition = build_sequence_partition(leaves,np,cost,cmax)
    return partition
end

# TODO: benchmark the different partitioning strategies for gemv. Is the hilber
# partition really faster than the simpler alternatives (row partition, col partition)?
function row_partitioning(H::HMatrix,np=Threads.nthreads())
    # sort the leaves by their row index
    leaves = filter(x -> isleaf(x),H)
    row_indices = map(leaves) do leaf
        # use the center of the leaf as a cartesian index
        i,j = pivot(leaf)
        return i
    end
    p = sortperm(row_indices)
    permute!(leaves,p)
    # now compute a quasi-optimal partition of leaves based `cost_mv`
    cmax = find_optimal_cost(leaves,np,cost_mv,1)
    partition = build_sequence_partition(leaves,np,cost_mv,cmax)
    return partition
end

function col_partitioning(H::HMatrix,np=Threads.nthreads())
    # sort the leaves by their row index
    leaves = filter(x -> isleaf(x),H)
    row_indices = map(leaves) do leaf
        # use the center of the leaf as a cartesian index
        i,j = pivot(leaf)
        return j
    end
    p = sortperm(row_indices)
    permute!(leaves,p)
    # now compute a quasi-optimal partition of leaves based `cost_mv`
    cmax = find_optimal_cost(leaves,np,cost_mv,1)
    partition = build_sequence_partition(leaves,np,cost_mv,cmax)
    return partition
end

function rmul!(R::RkMatrix, b::Number)
    m, n = size(R)
    if m > n
        rmul!(R.B, conj(b))
    else
        rmul!(R.A, b)
    end
    return R
end

function rmul!(H::HMatrix, b::Number)
    b == true && (return H) # short circuit. If inlined, rmul!(H,true) --> no-op
    if hasdata(H)
        rmul!(data(H), b)
    end
    for child in children(H)
        rmul!(child, b)
    end
    return H
end

"""
    _cost_gemv(A::Union{Matrix,SubArray,Adjoint})

A proxy for the computational cost of a matrix/vector product.
"""
function _cost_gemv(R::RkMatrix)
    rank(R)*sum(size(R))
end
function _cost_gemv(M::Matrix)
    length(M)
end
function _cost_gemv(H::HMatrix)
    acc = 0.0
    if isleaf(H)
        acc += _cost_gemv(data(H))
    else
        for c in children(H)
            acc += cost_gemv(c)
        end
    end
    return acc
end

max_stars_repo_name: WaveProp/HMatrices

max_issues_repo_name: WaveProp/HMatrices

max_forks_repo_name: WaveProp/HMatrices


text:
using Pkg
Pkg.add("CSV")
Pkg.add("SpecialFunctions")

using CSV
using LinearAlgebra
using SpecialFunctions

X = CSV.read("C:/Users/dhojo/Desktop/X.csv")
y = CSV.read("C:/Users/dhojo/Desktop/y.csv")
n = size(X)[1]

X = convert(Matrix,X)
y = convert(Matrix,y)

r_0 = 0.06
s_0 = 6
Q_0 = [r_0 0;0 s_0]
a_0 = 6
b_0 = 600^2

## the exact value of marginal likelihood

M = X' * X + Q_0
R = I(n) - X * inv(M) * X'
ML = π^(-n/2) * b_0^(a_0/2) * (gamma((n+a_0)/2) / gamma(a_0/2)) *
  (det(Q_0)^(1/2) / det(M)^(1/2)) * ((y'*R*y)[1,1]+b_0)^(-(n+a_0)/2)
log(ML)

logML = -(n/2)*log(pi) + (a_0/2)*log(b_0) + log( gamma((n+a_0)/2)/gamma(a_0/2) ) +
  (1/2)*log( det(Q_0)/det(M) ) - ((n+a_0)/2)*log((y'*R*y)[1,1]+b_0)
logML


text:
"""
DimensionalData.jl provides types and methods for indexing with named dimensions,
using named dimensions in Base and Statistics methods instead of Integer dims,
and selecting data from dimension values instead of using indices directly.

Dimensions are simply types that wrap values. They both store dimension values
and are used for dimension lookup or indices, ranges or dimension number.
`X`, Y`, `Z` and `Time` are the unexported defaults, add this line to use them:
```julia
using DimensionalData: X, Y, Z, Time
```

Selectors find indices in the dimension based on values `At`, `Near`, or `Between`
the passed in value(s).

These are some examples of valid syntax:

```julia
# Indexing with dim wrappers
a[X(1:10), Y(1:4)]
# Dim wrappers and a selector
a[X(1:10), Y<|At(25.7)]
# Unitful.jl selectors. Without dim wrappers selectors must be in the right order
a[Near(23"s"), Between(10.5u"m", 50.5u"m")]
# Dim type used instead of a dimension number
mean(a; dims=X)
```
"""
module DimensionalData

using RecipesBase, Statistics

using Base: tail, OneTo

export AbstractDimension, Dim

export Selector, Near, Between, At

export AbstractDimensionalArray, DimensionalArray

export dims, refdims, metadata, name, shortname,
       val, label, units, order, bounds, <|

include("interface.jl")
include("dimension.jl")
include("selector.jl")
include("array.jl")
include("methods.jl")
include("primitives.jl")
include("utils.jl")
include("plotrecipes.jl")

end


text:
using QOptCtrl

I = pauli_I
X = pauli_X
Y = pauli_Y
Z = pauli_Z

J = 1

prob = StateTransfer(
    J*kron(Z, Z), # drift Hamiltonian
    [kron(X, I), kron(Y, I), kron(I, X), kron(I, Y)], # control Hamiltonians
    10, # number of steps
    4, # transfer time
    kron(ket_0, ket_0), # initial state |00>
    kron(ket_0, ket_1), # target state |01>
    rand(10, 4) #initial guess
)

init_ctrl_arr,  final_ctrl_arr = OPTIMIZE(prob)
visualise_pulses(init_ctrl_arr, final_ctrl_arr)


text:
include("../common/layers.jl")

using StatsBase

mutable struct EmbeddingDot
    embed
    params
    grads
    cache
    function EmbeddingDot(W)
        self = new()
        self.embed = Embedding(W)
        self.params = self.embed.params
        self.grads = self.embed.grads
        self.cache = nothing
        return self
    end
end

function forward!(self::EmbeddingDot, h, idx)
    target_W = forward!(self.embed, idx)
    out = sum(target_W .* h, dims = 2)
    self.cache = (h, target_W)
    return out
end

function backward!(self::EmbeddingDot, dout)
    h, target_W = self.cache
    dout = reshape(dout, size(dout, 1), 1)

    dtarget_W = dout .* h
    backward!(self.embed, dtarget_W)
    dh = dout .* target_W
    return dh
end

mutable struct UnigramSampler
    sample_size
    vocab_size
    word_p
    function UnigramSampler(corpus, power, sample_size)
        self = new()
        self.sample_size = sample_size
        self.vocab_size = nothing
        self.word_p = nothing

        counts = Dict()
        for word_id in corpus
            if nothing == get(counts, word_id, nothing)
                get!(counts, word_id, 1)
            else
                counts[word_id] += 1
            end
        end

        vocab_size = length(counts)
        self.vocab_size = vocab_size

        self.word_p = zeros(vocab_size)
        for i = 1:vocab_size
            self.word_p[i] = counts[i]
        end

        self.word_p = self.word_p .^ power
        self.word_p ./= sum(self.word_p)

        return self
    end
end

function get_negative_sample(self::UnigramSampler, target)
    batch_size = size(target, 1)

    negative_sample = zeros(Int32, batch_size, self.sample_size)

    for i = 1:batch_size
        p = copy(self.word_p)
        target_idx = target[i]
        p[target_idx] = 0
        p ./= sum(p)
        negative_sample[i, :] = sample(self.vocab_size, ProbabilityWeights(p), self.sample_size, replace = false)
    end

    # TODO: GPUで処理する際の演算

    return negative_sample
end

mutable struct NegativeSamplingLoss <: AbstractLayer
    sample_size
    sampler
    loss_layers
    embed_dot_layers
    params
    grads
    function NegativeSamplingLoss(W, corpus; power = 0.75, sample_size = 5)
        self = new()
        self.sample_size = sample_size
        self.sampler = UnigramSampler(corpus, power, sample_size)
        self.loss_layers = [SigmoidWithLoss() for _ = 1:(sample_size+1)]
        self.embed_dot_layers = [EmbeddingDot(W) for _ = 1:(sample_size+1)]

        self.params, self.grads = [], []
        for layer in self.embed_dot_layers
            push!(self.params, layer.params)
            push!(self.grads, layer.grads)
        end
        return self
    end
end

function forward!(layer::NegativeSamplingLoss, h, target)
    batch_size = size(target, 1)
    negative_sample = get_negative_sample(layer.sampler, target)

    # 正例の順伝播
    score = forward!(layer.embed_dot_layers[1], target)
    correct_label = ones(Int32, batch_size)
    loss = forward!(layer.loss_layers[1], score, correct_label)

    # 負例の順伝播
    negative_label = zeros(Int32, batch_size)
    for i = 1:sample_size
        negative_target = negative_sample[:, i]
        score = forward!(layer.embed_dot_layers[1+i], h, negative_target)
        loss .+= forward!(layer.loss_layers[1+i], score, negative_label)
    end

    return loss
end

function backward!(layer::NegativeSamplingLoss; dout = 1)
    dh = 0
    for (l0, l1) in zip(layer.loss_layers, layer.embed_dot_layers)
        dscore = backward!(l0, dout)
        hf .+= backward!(l1, dscore)
    end

    return dh
end


text:
# Compressor Codec
# ================

struct ZstdCompressor <: TranscodingStreams.Codec
    cstream::CStream
    level::Int
end

function Base.show(io::IO, codec::ZstdCompressor)
    print(io, summary(codec), "(level=$(codec.level))")
end

# Same as the zstd command line tool (v1.2.0).
const DEFAULT_COMPRESSION_LEVEL = 3

"""
    ZstdCompressor(;level=$(DEFAULT_COMPRESSION_LEVEL))

Create a new zstd compression codec.

Arguments
---------
- `level`: compression level (1..$(MAX_CLEVEL))
"""
function ZstdCompressor(;level::Integer=DEFAULT_COMPRESSION_LEVEL)
    if !(1 ≤ level ≤ MAX_CLEVEL)
        throw(ArgumentError("level must be within 1..$(MAX_CLEVEL)"))
    end
    return ZstdCompressor(CStream(), level)
end

const ZstdCompressorStream{S} = TranscodingStream{ZstdCompressor,S} where S<:IO

"""
    ZstdCompressorStream(stream::IO; kwargs...)

Create a new zstd compression stream (see `ZstdCompressor` for `kwargs`).
"""
function ZstdCompressorStream(stream::IO; kwargs...)
    x, y = splitkwargs(kwargs, (:level,))
    return TranscodingStream(ZstdCompressor(;x...), stream; y...)
end

# Methods
# -------

function TranscodingStreams.initialize(codec::ZstdCompressor)
    code = initialize!(codec.cstream, codec.level)
    if iserror(code)
        zstderror(codec.cstream, code)
    end
    return
end

function TranscodingStreams.finalize(codec::ZstdCompressor)
    if codec.cstream.ptr != C_NULL
        code = free!(codec.cstream)
        if iserror(code)
            zstderror(codec.cstream, code)
        end
        codec.cstream.ptr = C_NULL
    end
    return
end

function TranscodingStreams.startproc(codec::ZstdCompressor, mode::Symbol, error::Error)
    code = reset!(codec.cstream, 0 #=unknown source size=#)
    if iserror(code)
        error[] = ErrorException("zstd error")
        return :error
    end
    return :ok
end

function TranscodingStreams.process(codec::ZstdCompressor, input::Memory, output::Memory, error::Error)
    cstream = codec.cstream
    cstream.ibuffer.src = input.ptr
    cstream.ibuffer.size = input.size
    cstream.ibuffer.pos = 0
    cstream.obuffer.dst = output.ptr
    cstream.obuffer.size = output.size
    cstream.obuffer.pos = 0
    if input.size == 0
        code = finish!(cstream)
    else
        code = compress!(cstream)
    end
    Δin = Int(cstream.ibuffer.pos)
    Δout = Int(cstream.obuffer.pos)
    if iserror(code)
        error[] = ErrorException("zstd error")
        return Δin, Δout, :error
    else
        return Δin, Δout, input.size == 0 && code == 0 ? :end : :ok
    end
end