text:
const SWISS_DICT_LOAD_FACTOR = 0.97
const _u8x16 = NTuple{16, VecElement{UInt8}}

"""
    SwissDict([itr])

`SwissDict{K,V}()` constructs a hash table with keys of type `K` and values of type `V`.
Keys are compared with [`isequal`](@ref) and hashed with [`hash`](@ref).
Given a single iterable argument, constructs a [`SwissDict`](@ref) whose key-value pairs
are taken from 2-tuples `(key,value)` generated by the argument.

# Examples
```jldoctest
julia> SwissDict([("A", 1), ("B", 2)])
SwissDict{String, Int64} with 2 entries:
  "A" => 1
  "B" => 2
```

Alternatively, a sequence of pair arguments may be passed.

```jldoctest
julia> SwissDict("A"=>1, "B"=>2)
SwissDict{String, Int64} with 2 entries:
  "A" => 1
  "B" => 2
```
"""
mutable struct SwissDict{K,V} <: AbstractDict{K,V}
    slots::Vector{_u8x16}
    keys::Vector{K}
    vals::Vector{V}
    nbfull::Int
    count::Int
    age::UInt
    idxfloor::Int  # an index <= the indices of all used slots

    function SwissDict{K,V}() where {K, V}
        new(fill(_expand16(0x00),1), Vector{K}(undef, 16), Vector{V}(undef, 16), 0, 0, 0, 1)
    end
    function SwissDict{K,V}(d::SwissDict{K,V}) where {K, V}
        new(copy(d.slots), copy(d.keys), copy(d.vals), d.nbfull, d.count, d.age,
            d.idxfloor)
    end
    function SwissDict{K, V}(slots, keys, vals, nbfull, count, age, idxfloor) where {K, V}
        new(slots, keys, vals, nbfull, count, age, idxfloor)
    end
end
function SwissDict{K,V}(kv) where {K, V}
    h = SwissDict{K,V}()
    for (k,v) in kv
        h[k] = v
    end
    return h
end
SwissDict{K,V}(p::Pair) where {K,V} = setindex!(SwissDict{K,V}(), p.second, p.first)
function SwissDict{K,V}(ps::Pair...) where {K, V}
    h = SwissDict{K,V}()
    sizehint!(h, length(ps))
    for p in ps
        h[p.first] = p.second
    end
    return h
end
SwissDict() = SwissDict{Any,Any}()
SwissDict(kv::Tuple{}) = SwissDict()
Base.copy(d::SwissDict) = SwissDict(d)
Base.empty(d::SwissDict, ::Type{K}, ::Type{V}) where {K, V} = SwissDict{K, V}()

SwissDict(ps::Pair{K,V}...) where {K,V} = SwissDict{K,V}(ps)
SwissDict(ps::Pair...)                  = SwissDict(ps)

function SwissDict(kv)
    try
        dict_with_eltype((K, V) -> SwissDict{K, V}, kv, eltype(kv))
    catch e
        if !isiterable(typeof(kv)) || !all(x->isa(x,Union{Tuple,Pair}),kv)
            throw(ArgumentError("SwissDict(kv): kv needs to be an iterator of tuples or pairs"))
        else
            rethrow(e)
        end
    end
end

# SIMD utilities
@inline _expand16(u::UInt8) = ntuple(i->VecElement(u), Val(16))
_blsr(i::UInt32)= i & (i-Int32(1))

@inline _vcmp_eq(u::_u8x16, v::_u8x16) = Core.Intrinsics.llvmcall(("""
%cmp = icmp eq <16 x i8> %0, %1
%cmp16 = bitcast <16 x i1> %cmp to i16
%res = zext i16 %cmp16 to i32
ret i32 %res
"""), UInt32, Tuple{_u8x16,_u8x16}, u, v)

@inline _vcmp_le(u::_u8x16, v::_u8x16) = Core.Intrinsics.llvmcall(("""
%cmp = icmp ule <16 x i8> %0, %1
%cmp16 = bitcast <16 x i1> %cmp to i16
%res = zext i16 %cmp16 to i32
ret i32 %res
"""), UInt32, Tuple{_u8x16,_u8x16}, u, v)

@inline function _prefetchr(p::Ptr)
    ccall("llvm.prefetch", llvmcall, Cvoid, (Ref{Int8}, Int32, Int32, Int32), Ptr{Int8}(p), 0, 3, 1)
end

@inline function _prefetchw(p::Ptr)
    ccall("llvm.prefetch", llvmcall, Cvoid, (Ref{Int8}, Int32, Int32, Int32), Ptr{Int8}(p), 1, 3, 1)
end

@inline function _hashtag(u::Unsigned)
    #extracts tag between 0x02 and 0xff from lower bits, rotates tag bits to front
    u = u % UInt
    tag = u % UInt8
    if UInt === UInt64
        hi = ((u>>8) | (u<<56)) % Int
    else
        hi = ((u>>8) | (u<<24)) % Int
    end
    tag = tag > 1 ? tag : tag+0x02
    return (hi, tag)
end

Base.@propagate_inbounds function _slotget(slots::Vector{_u8x16}, i::Int)
    @boundscheck 0 < i <= length(slots)*16 || throw(BoundsError(slots, 1 + (i-1)>>4))
    GC.@preserve slots begin
        return unsafe_load(convert(Ptr{UInt8}, pointer(slots)), i)
    end
end

Base.@propagate_inbounds function _slotset!(slots::Vector{_u8x16},  v::UInt8, i::Int)
   @boundscheck 0 < i <= length(slots)*16 || throw(BoundsError(slots, 1 + (i-1)>>4))
    GC.@preserve slots begin
        return unsafe_store!(convert(Ptr{UInt8}, pointer(slots)), v, i)
    end
end

@inline function _find_candidates(v::_u8x16, tag::UInt8)
    match = _vcmp_eq(v, _expand16(tag))
    return (match, v[16].value === 0x00)
end

@inline _find_free(v::_u8x16) = _vcmp_le(v, _expand16(UInt8(1)))

# Basic operations

# get the index where a key is stored, or -1 if not present
ht_keyindex(h::SwissDict, key) = ht_keyindex(h::SwissDict, key, _hashtag(hash(key))...)
function ht_keyindex(h::SwissDict, key, i0, tag)
    slots = h.slots
    keys = h.keys
    sz = length(slots)
    i = i0 & (sz-1)
    _prefetchr(pointer(h.keys, i*16+1))
    _prefetchr(pointer(h.vals, i*16+1))
    #Todo/discuss: _prefetchr(pointer(h.keys, i*16+9))?
    @inbounds while true
        msk = slots[i+1]
        cands, done = _find_candidates(msk, tag)
        while cands != 0
            off = trailing_zeros(cands)
            idx = i*16 + off + 1
            isequal(keys[idx], key) && return idx
            cands = _blsr(cands)
        end
        done && break
        i = (i+1) & (sz-1)
    end
    return -1
end

# get the index where a key is stored, or -pos if not present
# and the key would be inserted at pos
# This version is for use by setindex! and get!. It never rehashes.
ht_keyindex2!(h::SwissDict, key) = ht_keyindex2!(h, key, _hashtag(hash(key))...)
@inline function ht_keyindex2!(h::SwissDict, key, i0, tag)
    slots = h.slots
    keys = h.keys
    sz = length(slots)
    i = i0 & (sz-1)
    _prefetchw(pointer(h.keys, i*16+1))
    _prefetchw(pointer(h.vals, i*16+1))
    #Todo/discuss: _prefetchr(pointer(h.keys, i*16+9))?
    @inbounds while true
        msk = slots[i+1]
        cands, done = _find_candidates(msk, tag)
        while cands != 0
            off = trailing_zeros(cands)
            idx = i*16 + off + 1
            isequal(keys[idx], key) && return idx, tag
            cands = _blsr(cands)
        end
        done && break
        i = (i+1) & (sz-1)
    end
    i = i0 & (sz-1)
    @inbounds while true
        msk = slots[i+1]
        cands = _find_free(msk)
        if cands != 0
            off = trailing_zeros(cands)
            idx = i*16 + off + 1
            return -idx, tag
        end
        i = (i+1) & (sz-1)
    end
end

function _setindex!(h::SwissDict, v, key, index, tag)
    @inbounds h.keys[index] = key
    @inbounds h.vals[index] = v
    h.count += 1
    h.age += 1
    so =  _slotget(h.slots, index)
    h.nbfull += (iszero(index & 0x0f) & (so==0x00))
    _slotset!(h.slots, tag, index)
    if index < h.idxfloor
        h.idxfloor = index
    end
    maybe_rehash_grow!(h)
end

function _delete!(h::SwissDict{K,V}, index) where {K,V}
    # Caller is responsible for maybe shrinking the SwissDict after the deletion.
    isbitstype(K) || isbitsunion(K) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), h.keys, index-1)
    isbitstype(V) || isbitsunion(V) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), h.vals, index-1)
    isboundary = iszero(index & 0x0f) #boundaries: 16, 32, ...
    @inbounds _slotset!(h.slots, ifelse(isboundary, 0x01, 0x00), index)
    h.count -= 1
    h.age += 1
    maybe_rehash_shrink!(h)
end

# fast iteration over active slots.
function _iterslots(h::SwissDict, start::Int)
    i0 = ((start-1) & (length(h.keys)-1))>>4 + 1
    off = (start-1) & 0x0f
    @inbounds sl = _find_free(h.slots[i0>>4 + 1])
    sl = ((~sl & 0xffff)>>off) << off
    return _iterslots(h, (i0, sl))
end

function _iterslots(h::SwissDict, state)
    i, sl = state
    while iszero(sl)
        i += 1
        i <= length(h.slots) || return nothing
        @inbounds msk = h.slots[i]
        sl = _find_free(msk)
        sl = (~sl & 0xffff)
    end
    return ((i-1)*16 + trailing_zeros(sl) + 1, (i, _blsr(sl)))
end

# Dictionary resize logic:
# Guarantee 40% of buckets and 15% of entries free, and at least 25% of entries filled
# growth when > 85% entries full or > 60% buckets full, shrink when <25% entries full.
# >60% bucket full should be super rare outside of very bad hash collisions or
# super long-lived Dictionaries (expected 0.85^16 = 7% buckets full at 85% entries full).
# worst-case hysteresis: shrink at 25% vs grow at 30% if all hashes collide.
# expected hysteresis is 25% to 42.5%.
function maybe_rehash_grow!(h::SwissDict)
        sz = length(h.keys)
        if h.count > sz * SWISS_DICT_LOAD_FACTOR || (h.nbfull-1) * 10 > sz * 6
            rehash!(h, sz<<2)
        end
    end

function maybe_rehash_shrink!(h::SwissDict)
   sz = length(h.keys)
   if h.count*4 < sz && sz > 16
       rehash!(h, sz>>1)
   end
end

function Base.sizehint!(d::SwissDict, newsz)
    newsz = _tablesz(newsz*2)  # *2 for keys and values in same array
    oldsz = length(d.keys)
    # grow at least 25%
    if newsz < (oldsz*5)>>2
        return d
    end
    rehash!(d, newsz)
end

function rehash!(h::SwissDict{K,V}, newsz = length(h.keys)) where {K, V}
    olds = h.slots
    oldk = h.keys
    oldv = h.vals
    sz = length(oldk)
    newsz = _tablesz(newsz)
    (newsz*SWISS_DICT_LOAD_FACTOR) > h.count || (newsz <<= 1)
    h.age += 1
    h.idxfloor = 1
    if h.count == 0
        resize!(h.slots, newsz>>4)
        fill!(h.slots, _expand16(0x00))
        resize!(h.keys, newsz)
        resize!(h.vals, newsz)
        h.nbfull = 0
        return h
    end
    nssz = newsz>>4
    slots = fill(_expand16(0x00), nssz)
    keys = Vector{K}(undef, newsz)
    vals = Vector{V}(undef, newsz)
    age0 = h.age
    nbfull = 0
    is = _iterslots(h, 1)
    count = 0
    @inbounds while is !== nothing
        i, s = is
        k = oldk[i]
        v = oldv[i]
        i0, t = _hashtag(hash(k))
        i = i0 & (nssz-1)
        idx = 0
        while true
            msk = slots[i + 1]
            cands = _find_free(msk)
            if cands != 0
                off = trailing_zeros(cands)
                idx = i*16 + off + 1
                break
            end
            i = (i+1) & (nssz-1)
        end
        _slotset!(slots, t, idx)
        keys[idx] = k
        vals[idx] = v
        nbfull += iszero(idx & 0x0f)
        count += 1
        if h.age != age0
            return rehash!(h, newsz)
        end
        is = _iterslots(h, s)
    end
    h.slots = slots
    h.keys = keys
    h.vals = vals
    h.nbfull = nbfull
    @assert h.age == age0
    @assert h.count == count
    return h
end

Base.isempty(t::SwissDict) = (t.count == 0)
Base.length(t::SwissDict) = t.count

"""
    empty!(collection) -> collection

Remove all elements from a `collection`.

# Examples
```jldoctest
julia> A = SwissDict("a" => 1, "b" => 2)
SwissDict{String, Int64} with 2 entries:
  "a" => 1
  "b" => 2

julia> empty!(A);

julia> A
SwissDict{String, Int64}()
```
"""
function Base.empty!(h::SwissDict{K,V}) where {K, V}
    fill!(h.slots, _expand16(0x00))
    sz = length(h.keys)
    empty!(h.keys)
    empty!(h.vals)
    resize!(h.keys, sz)
    resize!(h.vals, sz)
    h.nbfull = 0
    h.count = 0
    h.age += 1
    h.idxfloor = 1
    return h
end

function Base.setindex!(h::SwissDict{K,V}, v0, key0) where {K, V}
    key = convert(K, key0)
    _setindex!(h, v0, key)
end

function _setindex!(h::SwissDict{K,V}, v0, key::K) where {K, V}
    v = convert(V, v0)
    index, tag = ht_keyindex2!(h, key)

    if index > 0
        h.age += 1
        @inbounds h.keys[index] = key
        @inbounds h.vals[index] = v
    else
        _setindex!(h, v, key, -index, tag)
    end

    return h
end

"""
    get!(collection, key, default)

Return the value stored for the given key, or if no mapping for the key is present, store
`key => default`, and return `default`.

# Examples
```jldoctest
julia> d = SwissDict("a"=>1, "b"=>2, "c"=>3);

julia> get!(d, "a", 5)
1

julia> get!(d, "d", 4)
4

julia> d
SwissDict{String, Int64} with 4 entries:
  "a" => 1
  "b" => 2
  "c" => 3
  "d" => 4
```
"""
Base.get!(h::SwissDict{K,V}, key0, default) where {K,V} = get!(()->default, h, key0)

"""
    get!(f::Function, collection, key)

Return the value stored for the given key, or if no mapping for the key is present, store
`key => f()`, and return `f()`.

This is intended to be called using `do` block syntax:
```julia
get!(dict, key) do
    # default value calculated here
    time()
end
```
"""
function Base.get!(default::Callable, h::SwissDict{K,V}, key0) where {K, V}
    key = convert(K, key0)
    return _get!(default, h, key)
end

function _get!(default::Callable, h::SwissDict{K,V}, key::K) where {K, V}
    index, tag = ht_keyindex2!(h, key)

    index > 0 && return @inbounds h.vals[index]

    age0 = h.age
    v = convert(V, default())
    if h.age != age0
        index, tag = ht_keyindex2!(h, key)
    end
    if index > 0
        h.age += 1
        @inbounds h.keys[index] = key
        @inbounds h.vals[index] = v
    else
        _setindex!(h, v, key, -index, tag)
    end
    return v
end

function Base.getindex(h::SwissDict{K,V}, key) where {K, V}
    index = ht_keyindex(h, key)
    @inbounds return (index < 0) ? throw(KeyError(key)) : h.vals[index]::V
end

"""
    get(collection, key, default)

Return the value stored for the given key, or the given default value if no mapping for the
key is present.

# Examples
```jldoctest
julia> d = SwissDict("a"=>1, "b"=>2);

julia> get(d, "a", 3)
1

julia> get(d, "c", 3)
3
```
"""
function Base.get(h::SwissDict{K,V}, key, default) where {K, V}
    index = ht_keyindex(h, key)
    @inbounds return (index < 0) ? default : h.vals[index]::V
end

"""
    get(f::Function, collection, key)

Return the value stored for the given key, or if no mapping for the key is present, return
`f()`.  Use [`get!`](@ref) to also store the default value in the dictionary.

This is intended to be called using `do` block syntax

```julia
get(dict, key) do
    # default value calculated here
    time()
end
```
"""
function Base.get(default::Callable, h::SwissDict{K,V}, key) where {K, V}
    index = ht_keyindex(h, key)
    @inbounds return (index < 0) ? default() : h.vals[index]::V
end

"""
    haskey(collection, key) -> Bool

Determine whether a collection has a mapping for a given `key`.

# Examples
```jldoctest
julia> D = SwissDict('a'=>2, 'b'=>3)
SwissDict{Char, Int64} with 2 entries:
  'a' => 2
  'b' => 3

julia> haskey(D, 'a')
true

julia> haskey(D, 'c')
false
```
"""
Base.haskey(h::SwissDict, key) = (ht_keyindex(h, key) > 0)
Base.in(key, v::KeySet{<:Any, <:SwissDict}) = (ht_keyindex(v.dict, key) > 0)

"""
    getkey(collection, key, default)

Return the key matching argument `key` if one exists in `collection`, otherwise return `default`.

# Examples
```jldoctest
julia> D = SwissDict('a'=>2, 'b'=>3)
SwissDict{Char, Int64} with 2 entries:
  'a' => 2
  'b' => 3

julia> getkey(D, 'a', 1)
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)

julia> getkey(D, 'd', 'a')
'a': ASCII/Unicode U+0061 (category Ll: Letter, lowercase)
```
"""
function Base.getkey(h::SwissDict{K,V}, key, default) where {K, V}
    index = ht_keyindex(h, key)
    @inbounds return (index<0) ? default : h.keys[index]::K
end

function _pop!(h::SwissDict, index)
    @inbounds val = h.vals[index]
    _delete!(h, index)
    maybe_rehash_shrink!(h)
    return val
end

"""
    pop!(collection, key[, default])

Delete and return the mapping for `key` if it exists in `collection`, otherwise return
`default`, or throw an error if `default` is not specified.

# Examples
```jldoctest
julia> d = SwissDict("a"=>1, "b"=>2, "c"=>3);

julia> pop!(d, "a")
1

julia> pop!(d, "d")
ERROR: KeyError: key "d" not found
Stacktrace:
 [1] pop!(h::SwissDict{String, Int64}, key::String)
   @ DataStructures ~/.julia/dev/DataStructures/src/swiss_dict.jl:599
 [2] top-level scope
   @ none:1

julia> pop!(d, "e", 4)
4
```
"""
function Base.pop!(h::SwissDict, key)
    index = ht_keyindex(h, key)
    return index > 0 ? _pop!(h, index) : throw(KeyError(key))
end

function Base.pop!(h::SwissDict, key, default)
    index = ht_keyindex(h, key)
    return index > 0 ? _pop!(h, index) : default
end

function Base.pop!(h::SwissDict)
    isempty(h) && throw(ArgumentError("SwissDict must be non-empty"))
    is = _iterslots(h, h.idxfloor)
    @assert is !== nothing
    idx, s = is
    @inbounds key = h.keys[idx]
    @inbounds val = h.vals[idx]
    _delete!(h, idx)
    h.idxfloor = idx
    return key => val
end

"""
    delete!(collection, key)

Delete the mapping for the given key in a collection, and return the collection.

# Examples
```jldoctest
julia> d = SwissDict("a"=>1, "b"=>2)
SwissDict{String, Int64} with 2 entries:
  "a" => 1
  "b" => 2

julia> delete!(d, "b")
SwissDict{String, Int64} with 1 entry:
  "a" => 1
```
"""
function Base.delete!(h::SwissDict, key)
    index = ht_keyindex(h, key)
    if index > 0
        _delete!(h, index)
    end
    maybe_rehash_shrink!(h)
    return h
end

Base.@propagate_inbounds function Base.iterate(h::SwissDict, state = h.idxfloor)
    is = _iterslots(h, state)
    is === nothing && return nothing
    i, s = is
    @inbounds p = h.keys[i] => h.vals[i]
    return (p, s)
end

Base.@propagate_inbounds function Base.iterate(v::Union{KeySet{<:Any, <:SwissDict}, Base.ValueIterator{<:SwissDict}}, state=v.dict.idxfloor)
    is = _iterslots(v.dict, state)
    is === nothing && return nothing
    i, s = is
    return (v isa KeySet ? v.dict.keys[i] : v.dict.vals[i], s)
end


text:
"""
TODO: Description of logging system.
"""
module logger

importall common
using state, core, utility
export Logger, LoggerDefinition

"""
The base type used by all loggers.
"""
abstract Logger

"""
The base type used by all logger definitions.
"""
abstract LoggerDefinition

"""
Creates a fresh, empty output directory at the specified location, before
preparing each of the given loggers for the upcoming run.
"""
function prepare!(ls::Vector{Logger}, output::AbstractString)
  if isdir(output)
    rm(output; recursive = true)
  end
  mkdir(output)
  for l in ls
    prepare!(l, output)
  end
end
prepare!(l::Logger, output::AbstractString) = l

"""
Safely closes a list of provided loggers.
"""
close!(ls::Vector{Logger}) = for l in ls; close!(l); end

"""
Safely closes a given logger.
"""
close!(l::Logger) = l

"""
Calls a list of attached loggers with the state of the search.
"""
call!(ls::Vector{Logger}, s::State) = for l in ls; call!(l, s); end

"""
Calls a single logger with the state of the search.
"""
call!(l::Logger, s::State) = log!(l, s)

# Load each of the loggers.
include("logger/best_individual.jl")
include("logger/population_dump.jl")
include("logger/fitness.jl")

end

ext: jl


text:
struct PolicyValueRewardTest <: TrainingTest
    idx::Integer
    state::AbstractEnvState
    image::AbstractImage
    actions::Vector{<:Integer}
    value::Float32
    action::Integer
    reward::Float32
end

mutable struct PolicyValueRewardTests <: TrainingTests
    tests::Vector{PolicyValueRewardTest}
    images::AbstractArray{Float32, 4}
    actual_value::AbstractArray{Float32, 2}
    actual_reward::AbstractArray{Float32, 2}
    period::Int
    warmup::Int
    function PolicyValueRewardTests(tests; period=1, warmup=0)
        new(
            tests,
            tensorize(get_image.(tests)),
            Flux.unsqueeze(get_value.(tests), 1),
            Flux.unsqueeze(get_reward.(tests), 1),
            period,
            warmup,
        )
    end
end

function gettestfn(tests::PolicyValueRewardTests)
    function fn(n::Network)
        initial_output = initialinferencefn(tests)(n)
        acc, fails = getpolicyacc(tests, initial_output)
        value_mse = getvaluemse(tests, initial_output)
        recurrent_output = recurrentinferencefn(tests, initial_output)(n)
        reward_mse = getrewardmse(tests, recurrent_output)
        acc, fails, value_mse, reward_mse
    end
end

struct PolicyValueRewardTestsSVG
    tests::PolicyValueRewardTests
end
svg(x::PolicyValueRewardTests) = PolicyValueRewardTestsSVG(x)

function Base.show(io::IO, m::MIME"image/svg+xml", x::PolicyValueRewardTestsSVG)
    num_tests = length(x.tests)
    range = 1:(min(10, num_tests))
    subset = 1:num_tests
    subset_indices = subset[range]
    svgbatch = [gridsvg(x.tests[i]) for i in subset_indices]
    Base.show(io, m, GridSVG(svgbatch, main_title="RVP Tests"))
end


text:
# calling a function in a shared library:
lang = ccall( (:getenv, "libc"),Ptr{Uint8}, (Ptr{Uint8},), "LANGUAGE")
bytestring(lang) #> "en_US"

# test existence of library:
find_library(["libc"]) #> "libc"

max_stars_repo_path: Module 1/Chapter09/Old/callc.jl

max_issues_repo_path: Module 1/Chapter09/Old/callc.jl

max_forks_repo_path: Module 1/Chapter09/Old/callc.jl


text:
import .AxisArrays
const AA = AxisArrays

axissymbol(ax::AA.Axis{s}) where {s} = s

axissymbols(aa::AA.AxisArray) = map(axissymbol, AA.axes(aa))
axisvalues(aa::AA.AxisArray) = AA.axisvalues(aa)
values(aa::AA.AxisArray) = parent(aa)

function from_namedtuple(::Type{T}, nt::NamedTuple) where {T<:AA.AxisArray}
    return AA.AxisArray(nt.values; nt.axes...)::T
end

roottype(::Type{T}) where {T<:AA.AxisArray} = AA.AxisArray


text:
#= Monetary type, and low-level operations =#
using FixedPointDecimals

# Abstract class for Monetary-like things
"""
The abstract type of objects representing a single value in one currency, or a
collection of values in a set of currencies. These objects should behave like
`Monetary` or `Basket` objects.
"""
abstract type AbstractMonetary end

"""
A simpler variant of `Monetary` that is expected to eventually be the default.
The currency represented is part of the type and not the object. The value is
internally represented as a quantity of some number type. The usual way to
construct a `Currency` directly, if needed, is:

    Currency{:USD}(FixedDecimal{Int, 2}(1))  # 1.00 USD
"""
struct Currency{C, T} <: AbstractMonetary
    val::T

    (::Type{Currency{C}})(x::Real) where C = new{C,typeof(x)}(x)
end

"""
A representation of a monetary value, denominated in some currency. The currency
used is part of the type and not the object. The value is internally represented
as a quantity of some integer type. The usual way to construct a `Monetary`
directly, if needed, is:

    Monetary(:USD)      # 1.00 USD
    Monetary(:USD, 325) # 3.25 USD

Be careful about the decimal point, as the `Monetary` constructor takes an
integer, representing the number of smallest denominations of the currency.
Typically, this constructor is not called directly. It is easier to use the
`@usingcurrencies` macro and the `100USD` form instead.

Although this type is flexible enough to support values internally represented
as any integer type, such as `BigInt`, it is recommended to use the built-in
`Int` type on your architecture unless you need a bigger type. Do not mix
different kinds of internal types. To use a different internal representation,
change the type of the second argument to `Monetary`:

    Monetary(:USD, BigInt(100))

In some applications, the minor denomination of a currency is not precise
enough. It is sometimes useful to override the number of decimal points stored.
For these applications, a third type parameter can be provided, indicating the
number of decimal points to keep after the major denomination:

    Monetary{:USD, BigInt, 4}(10000)            # 1.0000 USD
    Monetary(:USD, BigInt(10000); precision=4)  # 1.0000 USD
"""
Monetary{C, I, f} = Currency{C, FixedDecimal{I, f}}

# TODO: deprecate this constructor
(::Type{Monetary{C, I, f}})(x::Integer) where {C, I, f} =
    Currency{C}(reinterpret(FixedDecimal{I,f}, x))
(::Type{Monetary{C, I, f}})(x::Real) where {C, I, f} =
    Currency{C}(FixedDecimal{I,f}(x))

function Monetary(T::Symbol, x; precision=decimals(T))
    if precision == -1
        throw(ArgumentError("Must provide precision for currency $T."))
    else
        Monetary{T, typeof(x), precision}(x)
    end
end

function Monetary(T::Symbol; precision=decimals(T), storage=Int)
    if precision == -1
        throw(ArgumentError("Must provide precision for currency $T."))
    else
        majorunit(Monetary{T, storage, precision})
    end
end

"""
    filltype(typ) → typ

Fill in default type parameters to get a fully-specified concrete type from a
partially-specified one.
"""
filltype(::Type{Monetary{T}}) where {T} = Monetary{T, Int, decimals(T)}
filltype(::Type{Monetary{T,U}}) where {T,U} = Monetary{T, U, decimals(T)}


text:
lines = readlines("input.txt")
heightmap = parse.(Int, reduce(hcat, collect.(lines)))
idx = CartesianIndices(heightmap)

o = CartesianIndex(1,1)
V = CartesianIndex(-1, 0), CartesianIndex(1, 0), CartesianIndex(0, -1), CartesianIndex(0, 1)
padded = fill(10, size(heightmap) .+ 2)
padded[idx .+ o] = heightmap[idx]

minimum = reduce(.&, [padded[idx .+ o .+ v] .> heightmap for v = V])
println(sum(minimum .* (heightmap .+ 1)))


text:
# ------------ Call sites ------------ #

@inline function (ctx::GenerateContext)(::typeof(cond),
                                        addr::Address,
                                        c::Function,
                                        args::Tuple,
                                        a::Function,
                                        a_args::Tuple,
                                        b::Function,
                                        b_args::Tuple)
    visit!(ctx, addr)
    ss = get_subselection(ctx, addr => :C)
    ret, cl, w = generate(ss, c, args...)
    if ret
        ss = get_subselection(ctx, addr => :A)
        branch_ret, branch_cl, branch_w = generate(ss, a, a_args...)
        br_tr = BranchTrace(cl, branch_cl)
        add_call!(ctx, addr, ConditionalBranchCallSite(br_tr, get_score(branch_cl) + get_score(cl), c, args, ret, a, a_args, branch_ret))
    else
        ss = get_subselection(ctx, addr => :B)
        branch_ret, branch_cl, branch_w = generate(ss, b, b_args...)
        br_tr = BranchTrace(cl, branch_cl)
        add_call!(ctx, addr, ConditionalBranchCallSite(br_tr, get_score(branch_cl) + get_score(cl), c, args, ret, b, b_args, ret))
    end
end


text:
using GLMakie, JLD2, FFTW, Statistics
using FourierCore, FourierCore.Grid, FourierCore.Domain, Printf
# jld_name = "high_res_tracer_"
# jld_name = "high_order_timestep_spatial_tracer_"
# jld_name = "delete_me_later"
# jld_name = "high_order_timestep_tracer_"
jld_name = "tracer_"
# WARNING!!!! due to bug need to multiply all θ̅a by 5/4
Ω = S¹(4π)^2
N = 2^8 # number of gridpoints
κ = 2 / 2^8 # fixed diffusivity # in general should load from file
grid = FourierGrid(N, Ω, arraytype = Array)
nodes, wavenumbers = grid.nodes, grid.wavenumbers
# in the future load from file
x = nodes[1]
y = nodes[2]
kˣ = wavenumbers[1]
kʸ = wavenumbers[2]

if jld_name == "tracer_"
    bugfix = 5 / 4
    extrapolate = true # only works if zeroth mode is included
else
    bugfix = 1
    extrapolate = false # only works if zeroth mode is included
end

effective_diffusivity = Float64[]
index_choices = 2:40
for index_choice in index_choices
    println("on index ", index_choice)
    jlfile = jldopen(jld_name * string(index_choice) * ".jld2", "a+")
    ψ_save = jlfile["ψ_save"]
    θ_save = jlfile["θ_save"]
    if jld_name == "tracer_"
        κ = 2 / 2^8
        kˣ = wavenumbers[1]
    else
        κ = jlfile["κ"]
        xa = jlfile["xnodes"]
        ya = jlfile["ynodes"]
        kˣ = jlfile["kˣ_wavenumbers"]
    end
    θ̅a = mean(jlfile["θ̅a"], dims = 2)[:] * bugfix # messed up for tracer_ okay for others

    # calculate effective diffusivity
    # factor of 2 comes from the fact that the domain is 4π (as opposed to 2π)
    # will need to renormalize by length(θ̅a)/2 to take inverse transform
    if index_choice > 1
        𝒦ᵉᶠᶠ = ((length(θ̅a) / maximum(real.(fft(θ̅a)))) / 2) / (kˣ[index_choice]^2) - κ
    else
        𝒦ᵉᶠᶠ = abs(mean(θ̅a))
    end
    push!(effective_diffusivity, 𝒦ᵉᶠᶠ)
    close(jlfile)
end

fig = Figure()
ax = Axis(fig[1, 1]; xlabel = "wavenumber", ylabel = "𝒦(k)", title = "Effective Diffusivity")
scatter!(ax, Array(kˣ[index_choices]), effective_diffusivity)

ylims!(ax, (0.0, 0.4))
xlims!(ax, (-0.1, kˣ[index_choices[end]+1]))
display(fig)

for i in eachindex(effective_diffusivity)
    txtstring = @sprintf "%1.3f" effective_diffusivity[i]
    if i == 1
        print("[", txtstring, ", ")

    elseif i == length(effective_diffusivity)
        print(txtstring, "]")
    else
        print(txtstring, ", ")
    end
end
##
if extrapolate
    # only works if zeroth mode is included
    # fit model for effective diffusivity by picking out the pattern in the tail
    tail_index = length(effective_diffusivity) - 1 # 19 was good
    k₀ = kˣ[tail_index+1]
    a = 1 / effective_diffusivity[tail_index]
    c = (1 / effective_diffusivity[tail_index+1] - a) / (kˣ[tail_index+2]^2 - kˣ[tail_index+1]^2)
    # model_𝒦ᵉᶠᶠ(k) = 1 / (a + c * (k^2 - k₀^2))
    slope_index = index_choices[end] - 4
    logk = log.(kˣ[2:length(effective_diffusivity)+1])
    slope = (log.(effective_diffusivity[slope_index]) - log.(effective_diffusivity[slope_index-10])) / (logk[slope_index] - logk[slope_index-10])
    fit_function(x) = slope * (x - logk[slope_index]) + log(effective_diffusivity[slope_index])

    model_𝒦ᵉᶠᶠ(k) = abs(k) ≥ eps(100.0) ? exp(fit_function(log(abs(k)))) : 0.0 # just to handle zero well

    lines!(ax, Array(kˣ[index_choices[2:end]]), model_𝒦ᵉᶠᶠ.(Array(kˣ[index_choices[2:end]])), color = :red)

    ##
    hr_N = 1024 * 2
    hr_grid = FourierGrid(hr_N, Ω, arraytype = Array)
    hr_x = hr_grid.nodes[1]
    hr_kˣ = hr_grid.wavenumbers[1]

    fft_scaling = length(hr_x) / 2
    wavenumberspace = model_𝒦ᵉᶠᶠ.(hr_kˣ[:]) * fft_scaling # copy(x[:]) .* 0.0
    wavenumberspace[1] = effective_diffusivity[1] * fft_scaling
    wavenumberspace[div(length(hr_x), 2)+1] = 0.0 # need this because of redefinition of kˣ

    effective_diffusivity_mod = copy(effective_diffusivity[2:end])
    NKeff = length(effective_diffusivity_mod)
    wavenumberspace[2:1+NKeff] .= effective_diffusivity_mod * fft_scaling
    wavenumberspace[length(hr_x)-NKeff+1:end] .= reverse(effective_diffusivity_mod) * fft_scaling
    realspace = real.(ifft(wavenumberspace))

    # realspace = real.(ifft([0.0, effective_diffusivity..., reverse(effective_diffusivity)...]))
    index_shift = div(length(hr_x), 2) # div(length(realspace) - 1, 2) # div(length(x),2)
    realspace_shifted = circshift(realspace, index_shift) # .- minimum(realspace)
    println("zero'th mode is chosen to be ", minimum(realspace))
    fig_kernel = Figure()
    ax_kernel = Axis(fig_kernel[1, 1]; xlabel = "x-x'", ylabel = "𝒦(x-x')", title = "Effective Diffusivity")
    lines!(ax_kernel, hr_x[:], realspace_shifted)
end

##
#=
fig = Figure()
ax = Axis(fig[1, 1])
logk = log.(kˣ[2:length(effective_diffusivity)+1])
scatter!(logk, log.(effective_diffusivity))
slope_index = index_choices[end] - 4
slope = (log.(effective_diffusivity[slope_index]) - log.(effective_diffusivity[slope_index-10]))/(logk[slope_index] - logk[slope_index-10] )
fit_function(x) = slope*(x - logk[slope_index]) + log(effective_diffusivity[slope_index])
lines!(logk, fit_function.(logk), color = :red, linewidth = 3)

fig = Figure()
ax = Axis(fig[1, 1])
logk = log.(kˣ[2:length(effective_diffusivity)+1])
scatter!(kˣ[2:length(effective_diffusivity)+1], effective_diffusivity)
fit_function(x) = -(x - 3) - 4
lines!(kˣ[2:length(effective_diffusivity)+1], exp.(fit_function.(logk)), color = :red, linewidth = 3)
=#

max_stars_repo_name: sandreza/FourierCore

max_issues_repo_name: sandreza/FourierCore

max_forks_repo_name: sandreza/FourierCore


text:
function add_additional_state_provider(obj::Propagator, arg0::AdditionalStateProvider)
    return jcall(obj, "addAdditionalStateProvider", void, (AdditionalStateProvider,), arg0)
end

function add_event_detector(obj::Propagator, arg0::EventDetector)
    return jcall(obj, "addEventDetector", void, (EventDetector,), arg0)
end

function clear_events_detectors(obj::Propagator)
    return jcall(obj, "clearEventsDetectors", void, ())
end

function get_additional_state_providers(obj::Propagator)
    return jcall(obj, "getAdditionalStateProviders", List, ())
end

function get_attitude_provider(obj::Propagator)
    return jcall(obj, "getAttitudeProvider", AttitudeProvider, ())
end

function get_events_detectors(obj::Propagator)
    return jcall(obj, "getEventsDetectors", Collection, ())
end

function get_frame(obj::Propagator)
    return jcall(obj, "getFrame", Frame, ())
end

function get_generated_ephemeris(obj::Propagator)
    return jcall(obj, "getGeneratedEphemeris", BoundedPropagator, ())
end

function get_initial_state(obj::Propagator)
    return jcall(obj, "getInitialState", SpacecraftState, ())
end

function get_managed_additional_states(obj::Propagator)
    return jcall(obj, "getManagedAdditionalStates", Vector{JString}, ())
end

function get_max_date(obj::BoundedPropagator)
    return jcall(obj, "getMaxDate", AbsoluteDate, ())
end

function get_min_date(obj::BoundedPropagator)
    return jcall(obj, "getMinDate", AbsoluteDate, ())
end

function get_mode(obj::Propagator)
    return jcall(obj, "getMode", jint, ())
end

function get_pv_coordinates(obj::PVCoordinatesProvider, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function is_additional_state_managed(obj::Propagator, arg0::JString)
    return jcall(obj, "isAdditionalStateManaged", jboolean, (JString,), arg0)
end

function propagate(obj::Propagator, arg0::AbsoluteDate)
    return jcall(obj, "propagate", SpacecraftState, (AbsoluteDate,), arg0)
end

function reset_initial_state(obj::Propagator, arg0::SpacecraftState)
    return jcall(obj, "resetInitialState", void, (SpacecraftState,), arg0)
end

function set_attitude_provider(obj::Propagator, arg0::AttitudeProvider)
    return jcall(obj, "setAttitudeProvider", void, (AttitudeProvider,), arg0)
end

function set_ephemeris_mode(obj::Propagator)
    return jcall(obj, "setEphemerisMode", void, ())
end

function set_ephemeris_mode(obj::Propagator, arg0::OrekitStepHandler)
    return jcall(obj, "setEphemerisMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::Propagator, arg0::OrekitStepHandler)
    return jcall(obj, "setMasterMode", void, (OrekitStepHandler,), arg0)
end

function set_master_mode(obj::Propagator, arg0::jdouble, arg1::OrekitFixedStepHandler)
    return jcall(obj, "setMasterMode", void, (jdouble, OrekitFixedStepHandler), arg0, arg1)
end

function set_slave_mode(obj::Propagator)
    return jcall(obj, "setSlaveMode", void, ())
end


text:
# old and deprecated wrappers

## deprecated in CUDNN 8.0

@cenum cudnnConvolutionFwdPreference_t::UInt32 begin
    CUDNN_CONVOLUTION_FWD_NO_WORKSPACE = 0
    CUDNN_CONVOLUTION_FWD_PREFER_FASTEST = 1
    CUDNN_CONVOLUTION_FWD_SPECIFY_WORKSPACE_LIMIT = 2
end

@cenum cudnnConvolutionBwdFilterPreference_t::UInt32 begin
    CUDNN_CONVOLUTION_BWD_FILTER_NO_WORKSPACE = 0
    CUDNN_CONVOLUTION_BWD_FILTER_PREFER_FASTEST = 1
    CUDNN_CONVOLUTION_BWD_FILTER_SPECIFY_WORKSPACE_LIMIT = 2
end

@cenum cudnnConvolutionBwdDataPreference_t::UInt32 begin
    CUDNN_CONVOLUTION_BWD_DATA_NO_WORKSPACE = 0
    CUDNN_CONVOLUTION_BWD_DATA_PREFER_FASTEST = 1
    CUDNN_CONVOLUTION_BWD_DATA_SPECIFY_WORKSPACE_LIMIT = 2
end

@checked function cudnnGetConvolutionBackwardDataAlgorithm(handle, wDesc, dyDesc, convDesc,
                                                           dxDesc, preference,
                                                           memoryLimitInBytes, algo)
    initialize_api()
    @runtime_ccall((:cudnnGetConvolutionBackwardDataAlgorithm, libcudnn()), cudnnStatus_t,
                   (cudnnHandle_t, cudnnFilterDescriptor_t, cudnnTensorDescriptor_t,
                    cudnnConvolutionDescriptor_t, cudnnTensorDescriptor_t,
                    cudnnConvolutionBwdDataPreference_t, Csize_t,
                    Ptr{cudnnConvolutionBwdDataAlgo_t}),
                   handle, wDesc, dyDesc, convDesc, dxDesc, preference, memoryLimitInBytes,
                   algo)
end

@checked function cudnnSetRNNDescriptor(handle, rnnDesc, hiddenSize, numLayers,
                                        dropoutDesc, inputMode, direction, mode, algo,
                                        mathPrec)
    initialize_api()
    @runtime_ccall((:cudnnSetRNNDescriptor, libcudnn()), cudnnStatus_t,
                   (cudnnHandle_t, cudnnRNNDescriptor_t, Cint, Cint,
                    cudnnDropoutDescriptor_t, cudnnRNNInputMode_t, cudnnDirectionMode_t,
                    cudnnRNNMode_t, cudnnRNNAlgo_t, cudnnDataType_t),
                   handle, rnnDesc, hiddenSize, numLayers, dropoutDesc, inputMode,
                   direction, mode, algo, mathPrec)
end

@checked function cudnnGetRNNDescriptor(handle, rnnDesc, hiddenSize, numLayers,
                                        dropoutDesc, inputMode, direction, mode, algo,
                                        mathPrec)
    initialize_api()
    @runtime_ccall((:cudnnGetRNNDescriptor, libcudnn()), cudnnStatus_t,
                   (cudnnHandle_t, cudnnRNNDescriptor_t, Ptr{Cint}, Ptr{Cint},
                    Ptr{cudnnDropoutDescriptor_t}, Ptr{cudnnRNNInputMode_t},
                    Ptr{cudnnDirectionMode_t}, Ptr{cudnnRNNMode_t}, Ptr{cudnnRNNAlgo_t},
                    Ptr{cudnnDataType_t}),
                   handle, rnnDesc, hiddenSize, numLayers, dropoutDesc, inputMode,
                   direction, mode, algo, mathPrec)
end

@checked function cudnnSetRNNDescriptor_v5(rnnDesc, hiddenSize, numLayers, dropoutDesc,
                                           inputMode, direction, mode, mathPrec)
    initialize_api()
    @runtime_ccall((:cudnnSetRNNDescriptor_v5, libcudnn()), cudnnStatus_t,
                   (cudnnRNNDescriptor_t, Cint, Cint, cudnnDropoutDescriptor_t,
                    cudnnRNNInputMode_t, cudnnDirectionMode_t, cudnnRNNMode_t,
                    cudnnDataType_t),
                   rnnDesc, hiddenSize, numLayers, dropoutDesc, inputMode, direction, mode,
                   mathPrec)
end

@checked function cudnnGetConvolutionForwardAlgorithm(handle, xDesc, wDesc, convDesc,
                                                      yDesc, preference,
                                                      memoryLimitInBytes, algo)
    initialize_api()
    @runtime_ccall((:cudnnGetConvolutionForwardAlgorithm, libcudnn()), cudnnStatus_t,
                   (cudnnHandle_t, cudnnTensorDescriptor_t, cudnnFilterDescriptor_t,
                    cudnnConvolutionDescriptor_t, cudnnTensorDescriptor_t,
                    cudnnConvolutionFwdPreference_t, Csize_t,
                    Ptr{cudnnConvolutionFwdAlgo_t}),
                   handle, xDesc, wDesc, convDesc, yDesc, preference, memoryLimitInBytes,
                   algo)
end

@checked function cudnnGetConvolutionBackwardFilterAlgorithm(handle, xDesc, dyDesc,
                                                             convDesc, dwDesc, preference,
                                                             memoryLimitInBytes, algo)
    initialize_api()
    @runtime_ccall((:cudnnGetConvolutionBackwardFilterAlgorithm, libcudnn()), cudnnStatus_t,
                   (cudnnHandle_t, cudnnTensorDescriptor_t, cudnnTensorDescriptor_t,
                    cudnnConvolutionDescriptor_t, cudnnFilterDescriptor_t,
                    cudnnConvolutionBwdFilterPreference_t, Csize_t,
                    Ptr{cudnnConvolutionBwdFilterAlgo_t}),
                   handle, xDesc, dyDesc, convDesc, dwDesc, preference, memoryLimitInBytes,
                   algo)
end


text:
using LinearAlgebra
include("../various_mle_things/incremental.jl")

"""
Computes remainder and minimizer for lsq problem min_x |hat(A)x-hat(b)|_2^2 where
ith row of A and b are discarded, given QR factorization of A,
remainder r and minimizer xstar of the full problem
"""
function lsq_loocv(A, AQR, r, xstar, i)
    n = size(A, 1)
    p = size(A, 2)
    @assert n>=p #more rows than columns
    @assert i>=1 && n>=i #i must be valid row
    Q = AQR.Q[1:n, 1:p] #get reduced QR factor
    R = AQR.R
    At = A'
    r = reshape(r, length(r), 1) #reshape into col vector
    xstar = reshape(xstar, length(xstar), 1) #reshape into col vector
    remainder = repeat(r, 1, 1) #initial value (r) is important to algorithm
    minimizer = repeat(xstar, 1, 1) #initial value (xstar) is important to algorithm
    rmv = @views remainder
    minv = @views minimizer
    Atv = @views At
    leverage_score = norm(Q[i, :])^2 #norm of ith row of Q squared
    gamma = r[i] / (1-leverage_score)
    minv[:, 1] .-= R\(R'\ Atv[:, i]) * gamma
    rmv[:, 1] .+= Q * (Q'[:, i]) * gamma
    return remainder, minimizer #col vectors
end

"""
LOOCV for least squares problem min_x |Ax-b|_2^2. Computes all n remainder and minimizer vectors. Primarily used for testing.
"""
function lsq_loocv(A, AQR, r, xstar)
    n = size(A, 1)
    p = size(A, 2)
    @assert n>=p #more rows than columns
    Q = AQR.Q[1:n, 1:p] #get reduced QR factor
    R = AQR.R
    At = A'
    r = reshape(r, length(r), 1) #reshape into col vector
    xstar = reshape(xstar, length(xstar), 1) #reshape into col vector
    remainders = repeat(r, 1, n) #length n b/c we have n subproblems
    minimizers = repeat(xstar, 1, n)
    rmv = @views remainders
    minv = @views minimizers
    Atv = @views At
    for i = 1:n
        leverage_score = norm(Q[i, :])^2 #norm of ith row of Q squared
        gamma = r[i] / (1-leverage_score)
        minv[:, i] .-= R\(R'\ Atv[:, i]) * gamma
        rmv[:, i] .+= Q * (Q'[:, i]) * gamma
    end
    return remainders, minimizers
end

function unitvec(n, i)
    v = zeros(n, 1)
    v[i] = 1
    return v
end

"""
loocv for linear system Kc = y, generalized to case where c is a matrix
    INPUTS: R: upper triangular Cholesky factor of K
            c: solution to Kc = y
            i: row up for deletion
"""
function lin_sys_loocv(c, R, i)
    n = size(R)[1]
    #println(n)
    if size(c, 2)==1
        c = reshape(c, length(c), 1)
    end
    @assert i <= size(c, 1)
    ret = Array{Float64, 2}(undef, size(c, 1) - 1, size(c, 2))
    for j = 1:size(ret, 2)
        #println(unitvec(n, i) )
        #println(R')
        ve = R'\ unitvec(n, i)
        ri = c[i, j] / norm(ve)^2
        ci = c[:, j] - ri * (R\ve)
        ret[:, j] = ci[[1:i-1;i+1:end]] #dimension of return vector is n-1, ith entry will be zero
    end
    return ret
end

"""
Use incremental cholesky in computation
"""
function lin_sys_loocv_IC(c::Array{T} where T, IC::IncrementalCholesky{Float64}, i::Int64)
    R = get_chol(IC).U
    L = get_chol(IC).L
    n = size(IC)[1]
    if size(c, 2)==1
        c = reshape(c, length(c), 1)
    end
    @assert i <= size(c, 1)
    ret = Array{Float64, 2}(undef, size(c, 1) - 1, size(c, 2))
    for j = 1:size(ret, 2)
        #println(unitvec(n, i) )
        #println(R')
        #ve = R'\ unitvec(n, i)
        ve = L\unitvec(n, i)
        ri = c[i, j] / norm(ve)^2
        ci = c[:, j] - ri * (R\ve)
        ret[:, j] = ci[[1:i-1;i+1:end]] #dimension of return vector is n-1, ith entry will be zero
    end
    return ret
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
@testset "qstate_test" begin
    @testset "MPSState & ExactState consistency test" begin
    	    mps = MPSState(5)
    	    exact = toExactState(mps)
    	    circuit = randomQCircuit(5, 40)
    	    apply!(mps, circuit)
    	    apply!(exact, circuit)
            config = rand(0:1, 5)
            for i =1:10
                @test abs(probability(mps, config)-probability(exact, config)) < 1e-10
            end
    end

    @testset "MPSState measure test" begin
        @testset "measure single qubit test" begin
            zero_state = MPSState(1)
            one_state = MPSState(1)
            apply!(one_state, X, 1)
            zero_result = measure!(zero_state, 1, 100)
            one_result = measure!(one_state, 1, 100)
            @test sum(zero_result) == 0
            @test prod(one_result) ==1
        end
        @testset "measure! single qubit test" begin
            state = MPSState(1)
            apply!(state, H, 1)
            result = collapse!(state, 1)
            result2 = measure!(state, 1, 100)
            if result==0
            	@test sum(result2) == 0
            else
            	@test prod(result2) ==1
            end
    	end
    	@testset "measure single qubit in X basis" begin
    	    state = MPSState(1)
            apply!(state, H, 1)
            basis = add!(QGateBlock(), H, 1)
            result = measure!(state, basis, 1, 1024)
            @test sum(result) == 0
    	end
        @testset "measure multiple qubits test" begin
            @testset "Throw away probability" for binary in [true, false]
                state = MPSState(2)
                apply!(state, H, 1)
                apply!(state, CNOT, [1,2])
                if binary
                    result = measure!(state, [1,2], 20; binary = binary)
                    for i = 1:20
                        @test result[i,1] == result[i,2]
                    end
                else
                    result = measure!(state, [1,2], 20; binary = binary)
                    for i=1:20
                        @test result[i]==0 || result[i]==3
                    end
                end
            end
            @testset "Keep probability as a by-product" begin
                state = MPSState(2)
                apply!(state, H, 1)
                apply!(state, CNOT, [1,2])
                    result, prob = measure!(state, [1,2], 20; probability=true)
                for i = 1:20
                    @test result[i,1] == result[i,2]
                    @test prob[i] ≈ 0.5
                end
            end
        end
        @testset "measure! multiple qubits test" begin
            state = MPSState(2)
            apply!(state, H, 1)
            apply!(state, CNOT, [1,2])
            result = collapse!(state, [1,2])
            result2 = measure!(state, [1,2], 100)
            for i = 1:100
                @test result2[i,:] == result
            end
        end
        @testset "measure! with given rng" begin
            state = MPSState(5)
            for i =1: 5
                apply!(state, H, rand(1:5))
                apply!(state, Y, rand(1:5))
            end
            rng = Random.MersenneTwister(1234)
            result = measure!(state, [1:5;], 5; rng = rng)
            Random.seed!(rng, 1234)
            result2 = measure!(state, [1:5;], 5; rng = rng)
            @test result == result2
        end
    end

    @testset "probability() test" begin
         state_m = MPSState(5)
         state_e = ExactState(5)
         apply!(state_m, H, 3)
         apply!(state_e, H, 3)
         config = rand(0:1, 5)
         @test probability(state_m, config) ≈ probability(state_e, config)
    end
end


text:
"""
The [Multinomial distribution](http://en.wikipedia.org/wiki/Multinomial_distribution)
generalizes the *binomial distribution*. Consider n independent draws from a Categorical
distribution over a finite set of size k, and let ``X = (X_1, ..., X_k)`` where ``X_i``
represents the number of times the element ``i`` occurs, then the distribution of ``X``
is a multinomial distribution. Each sample of a multinomial distribution is a k-dimensional
integer vector that sums to n.

The probability mass function is given by

```math
f(x; n, p) = \\frac{n!}{x_1! \\cdots x_k!} \\prod_{i=1}^k p_i^{x_i},
\\quad x_1 + \\cdots + x_k = n
```

```julia
Multinomial(n, p)   # Multinomial distribution for n trials with probability vector p
Multinomial(n, k)   # Multinomial distribution for n trials with equal probabilities
                    # over 1:k
```
"""
struct Multinomial{T<:Real} <: DiscreteMultivariateDistribution
    n::Int
    p::Vector{T}

    function Multinomial{T}(n::Integer, p::Vector{T}) where T
        if n < 0
            throw(ArgumentError("n must be a nonnegative integer."))
        end
        if !isprobvec(p)
            throw(ArgumentError("p = $p is not a probability vector."))
        end
        new{T}(round(Int, n), p)
    end
    Multinomial{T}(n::Integer, p::Vector{T}, ::NoArgCheck) where {T} = new{T}(round(Int, n), p)
end
Multinomial(n::Integer, p::Vector{T}) where {T<:Real} = Multinomial{T}(n, p)
Multinomial(n::Integer, k::Integer) = Multinomial{Float64}(round(Int, n), fill(1.0 / k, k))

# Parameters

ncategories(d::Multinomial) = length(d.p)
length(d::Multinomial) = ncategories(d)
probs(d::Multinomial) = d.p
ntrials(d::Multinomial) = d.n

params(d::Multinomial) = (d.n, d.p)
@inline partype(d::Multinomial{T}) where {T<:Real} = T

### Conversions
convert(::Type{Multinomial{T}}, d::Multinomial) where {T<:Real} = Multinomial(d.n, Vector{T}(d.p))
convert(::Type{Multinomial{T}}, n, p::Vector) where {T<:Real} = Multinomial(n, Vector{T}(p))

# Statistics

mean(d::Multinomial) = d.n .* d.p

function var(d::Multinomial{T}) where T<:Real
    p = probs(d)
    k = length(p)
    n = ntrials(d)

    v = Vector{T}(k)
    for i = 1:k
        @inbounds p_i = p[i]
        v[i] = n * p_i * (1 - p_i)
    end
    v
end

function cov(d::Multinomial{T}) where T<:Real
    p = probs(d)
    k = length(p)
    n = ntrials(d)

    C = Matrix{T}(k, k)
    for j = 1:k
        pj = p[j]
        for i = 1:j-1
            @inbounds C[i,j] = - n * p[i] * pj
        end

        @inbounds C[j,j] = n * pj * (1-pj)
    end

    for j = 1:k-1
        for i = j+1:k
            @inbounds C[i,j] = C[j,i]
        end
    end
    C
end

function mgf(d::Multinomial{T}, t::AbstractVector) where T<:Real
    p = probs(d)
    n = ntrials(p)
    s = zero(T)
    for i in 1:length(p)
        s += p[i] * exp(t[i])
    end
    return s^n
end

function cf(d::Multinomial{T}, t::AbstractVector) where T<:Real
    p = probs(d)
    n = ntrials(d)
    s = zero(Complex{T})
    for i in 1:length(p)
        s += p[i] * exp(im * t[i])
    end
    return s^n
end

function entropy(d::Multinomial)
    n, p = params(d)
    s = -lgamma(n+1) + n*entropy(p)
    for pr in p
        b = Binomial(n, pr)
        for x in 0:n
            s += pdf(b, x) * lgamma(x+1)
        end
    end
    return s
end

# Evaluation

function insupport(d::Multinomial, x::AbstractVector{T}) where T<:Real
    k = length(d)
    length(x) == k || return false
    s = 0.0
    for i = 1:k
        @inbounds xi = x[i]
        if !(isinteger(xi) && xi >= 0)
            return false
        end
        s += xi
    end
    return s == ntrials(d)  # integer computation would not yield truncation errors
end

function _logpdf(d::Multinomial, x::AbstractVector{T}) where T<:Real
    p = probs(d)
    n = ntrials(d)
    S = eltype(p)
    R = promote_type(T, S)
    s = R(lgamma(n + 1))
    t = zero(T)
    for i = 1:length(p)
        @inbounds xi = x[i]
        @inbounds p_i = p[i]
        t += xi
        s -= R(lgamma(xi + 1))
        @inbounds s += xi * log(p_i)
    end
    return ifelse(t == n, s, -R(Inf))
end

# Sampling

_rand!(d::Multinomial, x::AbstractVector{T}) where {T<:Real} = multinom_rand!(ntrials(d), probs(d), x)

sampler(d::Multinomial) = MultinomialSampler(ntrials(d), probs(d))

## Fit model

struct MultinomialStats <: SufficientStats
    n::Int  # number of trials in each experiment
    scnts::Vector{Float64}  # sum of counts
    tw::Float64  # total sample weight

    MultinomialStats(n::Int, scnts::Vector{Float64}, tw::Real) = new(n, scnts, Float64(tw))
end

function suffstats(::Type{Multinomial}, x::Matrix{T}) where T<:Real
    K = size(x, 1)
    n::T = zero(T)
    scnts = zeros(K)

    for j = 1:size(x,2)
        nj = zero(T)
        for i = 1:K
            @inbounds xi = x[i,j]
            @inbounds scnts[i] += xi
            nj += xi
        end

        if j == 1
            n = nj
        elseif nj != n
            error("Each sample in X should sum to the same value.")
        end
    end
    MultinomialStats(n, scnts, size(x,2))
end

function suffstats(::Type{Multinomial}, x::Matrix{T}, w::Array{Float64}) where T<:Real
    length(w) == size(x, 2) || throw(DimensionMismatch("Inconsistent argument dimensions."))

    K = size(x, 1)
    n::T = zero(T)
    scnts = zeros(K)
    tw = 0.

    for j = 1:size(x,2)
        nj = zero(T)
        @inbounds wj = w[j]
        tw += wj
        for i = 1:K
            @inbounds xi = x[i,j]
            @inbounds scnts[i] += xi * wj
            nj += xi
        end

        if j == 1
            n = nj
        elseif nj != n
            error("Each sample in X should sum to the same value.")
        end
    end
    MultinomialStats(n, scnts, tw)
end

fit_mle(::Type{Multinomial}, ss::MultinomialStats) = Multinomial(ss.n, ss.scnts * inv(ss.tw * ss.n))

function fit_mle(::Type{Multinomial}, x::Matrix{T}) where T<:Real
    ss = suffstats(Multinomial, x)
    Multinomial(ss.n, multiply!(ss.scnts, inv(ss.tw * ss.n)))
end

function fit_mle(::Type{Multinomial}, x::Matrix{T}, w::Array{Float64}) where T<:Real
    ss = suffstats(Multinomial, x, w)
    Multinomial(ss.n, multiply!(ss.scnts, inv(ss.tw * ss.n)))
end


text:
# This file is part of RLinearAlgebra.jl

module ProceduralTestLSVRRandCyclic

using Test, RLinearAlgebra, Random

Random.seed!(1010)

@testset "LSVR Random Cyclic -- Procedural" begin

    # Verify appropriate super type
    @test supertype(LinSysVecRowRandCyclic) == LinSysVecRowSampler

    # Test whether ordering is appropriately sampled, and changes after exhaustion
    @test let
        A = rand(10, 3)
        b = rand(10)
        x = rand(3)

        cyc = LinSysVecRowRandCyclic()

        # Generate random ordering
        α, β = RLinearAlgebra.sample(cyc, A, b, x, 1)
        order = copy(cyc.order)

        flag = true
        for j = 2:10
            α, β = RLinearAlgebra.sample(cyc, A, b, x, j)
            # Order should remain fixed
            flag = flag & (order == cyc.order)
        end

        # Order should change
        α, β = RLinearAlgebra.sample(cyc, A, b, x, 11)
        flag = flag & (order != cyc.order)
        flag = flag & (α == A[cyc.order[1],:])
        flag = flag & (β == b[cyc.order[1]])

        flag
    end

end

end


text:
function ReadData(fname::AbstractString)
    lines = []

    xMax, xMin = 0, 1
    yMax, yMin = 0, 1

    for line in eachline(fname)
        push!(lines, map(y-> map(x -> parse(Int, x), split(y, ",")), split(line, " -> ")))

        x0, y0 = lines[length(lines)][1]
        x1, y1 = lines[length(lines)][2]

        if min(x0, x1) < xMin
            xMin = min(x0, x1)
        end
        if max(x0, x1) > xMax
            xMax = max(x0, x1)
        end
        if min(y0, y1) < yMin
            yMin = min(y0, y1)
        end
        if max(y0, y1) > yMax
            yMax = max(y0, y1)
        end
    end

    return lines, (xMin, xMax), (yMin, yMax)
end

function IsVertical(line::Array)
    return line[1][1] == line[2][1]
end

function IsHorizontal(line::Array)
    return line[1][2] == line[2][2]
end

function IsDiagonal(line::Array)
    return abs(line[1][1]-line[2][1]) == abs(line[1][2]-line[2][2])
end

function Part1(data)
    lines = data[1]
    xMin, xMax = data[2]
    yMin, yMax = data[3]

    grid = zeros(Int, yMax+1, xMax+1)

    for line in lines
        if IsVertical(line)
            y0, y1 = 0, 0
            if line[1][2] > line[2][2]
                y0, y1 = line[2][2], line[1][2]
            else
                y0, y1 = line[1][2], line[2][2]
            end

            grid[(y0+1):(y1+1), line[1][1]+1] += fill(1, y1-y0+1)
        end

        if IsHorizontal(line)
            x0, x1 = 0, 0
            if line[1][1] > line[2][1]
                x0, x1 = line[2][1], line[1][1]
            else
                x0, x1 = line[1][1], line[2][1]
            end

            grid[line[1][2]+1, (x0+1):(x1+1)] += fill(1, x1-x0+1)
        end
    end

    #for i in 1:(yMax+1)
    #    println(grid[i,:])
    #end

    return count(x -> x>=2, grid)
end

function Part2(data)
    lines = data[1]
    xMin, xMax = data[2]
    yMin, yMax = data[3]

    grid = zeros(Int, yMax+1, xMax+1)

    for line in lines
        if IsVertical(line)
            y0, y1 = 0, 0
            if line[1][2] > line[2][2]
                y0, y1 = line[2][2], line[1][2]
            else
                y0, y1 = line[1][2], line[2][2]
            end

            grid[(y0+1):(y1+1), line[1][1]+1] += fill(1, y1-y0+1)

        elseif IsHorizontal(line)
            x0, x1 = 0, 0
            if line[1][1] > line[2][1]
                x0, x1 = line[2][1], line[1][1]
            else
                x0, x1 = line[1][1], line[2][1]
            end

            grid[line[1][2]+1, (x0+1):(x1+1)] += fill(1, x1-x0+1)

        else
            x0, x1 = line[1][1], line[2][1]
            y0, y1 = line[1][2], line[2][2]

            dx = x1 > x0 ? 1 : -1
            dy = y1 > y0 ? 1 : -1

            for i = 0:(abs(x1-x0))
                grid[y0 + dy*i + 1, x0 + dx*i + 1] += 1
            end
        end

    end

    #for i in 1:yMax+1
    #    println(grid[i,:])
    #end

    return count(x -> x>=2, grid)
end

#data = ReadData("../data/test5.txt")
data = ReadData("../data/day5.txt")

#println(Part1(data))
println(Part2(data))

# 0,0 -> 8,8
# 8,0 -> 0,8
# NOT 6,4 -> 2,0


text:
"""
$(TYPEDEF)

Adjacency describing edges per grid cell
"""
# abstract type CellEdges  <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing cells per grid edge
"""
# abstract type EdgeCells  <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing nodes per grid edge
"""
# abstract type EdgeNodes <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing cells per boundary or interior face
"""
abstract type BFaceCells <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing outer normals to boundary faces
"""
abstract type BFaceNormals <: AbstractGridComponent end

"""
$(TYPEDEF)

Adjacency describing edges per boundary or interior face
"""
abstract type BFaceEdges <: AbstractGridAdjacency end

"""
$(TYPEDEF)

Adjacency describing nodes per boundary or interior edge
"""
# abstract type BEdgeNodes <: AbstractGridAdjacency end

"""
$(SIGNATURES)

Prepare edge adjacencies (celledges, edgecells, edgenodes)

Currently depends on ExtendableSparse, we may want to remove this
adjacency.
"""
function prepare_edges!(grid::ExtendableGrid)
    Ti=eltype(grid[CellNodes])
    cellnodes=grid[CellNodes]
    geom=grid[CellGeometries][1]
    # Create cell-node incidence matrix
    cellnode_adj=asparse(cellnodes)

    # Create node-node incidence matrix for neigboring
    # nodes.
    nodenode_adj=cellnode_adj*transpose(cellnode_adj)

    # To get unique edges, we set the lower triangular part
    # including the diagonal to 0
    for icol=1:length(nodenode_adj.colptr)-1
        for irow=nodenode_adj.colptr[icol]:nodenode_adj.colptr[icol+1]-1
            if nodenode_adj.rowval[irow]>=icol
                nodenode_adj.nzval[irow]=0
            end
        end
    end
    dropzeros!(nodenode_adj)

    # Now we know the number of edges and
    nedges=length(nodenode_adj.nzval)

    if dim_space(grid)==2
        # Let us do the Euler test (assuming no holes in the domain)
        v=num_nodes(grid)
        e=nedges
        f=num_cells(grid)+1
        @assert v-e+f==2
    end

    if dim_space(grid)==1
        @assert nedges==num_cells(grid)
    end

    # Calculate edge nodes and celledges
    edgenodes=zeros(Ti,2,nedges)
    celledges=zeros(Ti,num_edges(geom),num_cells(grid))
    cen=local_celledgenodes(geom)

    for icell=1:num_cells(grid)
        for iedge=1:num_edges(geom)
            n1=cellnodes[cen[1,iedge],icell]
            n2=cellnodes[cen[2,iedge],icell]

            # We need to look in nodenod_adj for upper triangular part entries
            # therefore, we need to swap accordingly before looking
	       if (n1<n2)
		      n0=n1
		      n1=n2
		      n2=n0;
	       end

            for irow=nodenode_adj.colptr[n1]:nodenode_adj.colptr[n1+1]-1
                if nodenode_adj.rowval[irow]==n2
                    # If the coresponding entry has been found, set its
                    # value. Note that this introduces a different edge orientation
                    # compared to the one found locally from cell data
                    celledges[iedge,icell]=irow
                    edgenodes[1,irow]=n1
                    edgenodes[2,irow]=n2
                end
            end
        end
    end

    # Create sparse incidence matrix for the cell-edge adjacency
    celledge_adj=asparse(celledges)

    # The edge cell matrix is the transpose
    edgecell_adj=SparseMatrixCSC(transpose(celledge_adj))

    # Get the adjaency array from the matrix
    edgecells=zeros(Ti,2,nedges) ## for 3D we need more here!
    for icol=1:length(edgecell_adj.colptr)-1
        ii=1
        for irow=edgecell_adj.colptr[icol]:edgecell_adj.colptr[icol+1]-1
            edgecells[ii,icol]=edgecell_adj.rowval[irow]
            ii+=1
        end
    end

    grid[EdgeCells]=edgecells
    grid[CellEdges]=celledges
    grid[EdgeNodes]=edgenodes
    true
end

#ExtendableGrids.instantiate(grid, ::Type{CellEdges})=prepare_edges!(grid) && grid[CellEdges]
#ExtendableGrids.instantiate(grid, ::Type{EdgeCells})=prepare_edges!(grid) && grid[EdgeCells]
#ExtendableGrids.instantiate(grid, ::Type{EdgeNodes})=prepare_edges!(grid) && grid[EdgeNodes]

function prepare_bfacecells!(grid)
    cn   = grid[CellNodes]
    bn   = grid[BFaceNodes]
    dim  = dim_space(grid)

    abn=asparse(bn)
    # The maximum number of nodes adjacent to bfaces may
    # be less than the number of nodes in the grid.
    # Therefore we add a the missing number of rows.
    if abn.m<num_nodes(grid)
        abn=SparseMatrixCSC(num_nodes(grid), abn.n,abn.colptr,abn.rowval,abn.nzval)
    end

    abc  = asparse(cn)'*abn
    abcx = dropzeros!(SparseMatrixCSC(abc.m,abc.n, abc.colptr, abc.rowval,
                                      map( i-> i==dim,  abc.nzval)))
    grid[BFaceCells] = VariableTargetAdjacency(abcx)
    true
end

function prepare_bedges!(grid)
    Ti           = eltype(grid[CellNodes])
    bgeom        = grid[BFaceGeometries][1]
    bfacenodes   = grid[BFaceNodes]

    # Create bface-node incidence matrix
    bfacenode_adj = asparse(bfacenodes)
    nodenode_adj = bfacenode_adj*transpose(bfacenode_adj)

    # To get unique edges, we set the lower triangular part
    # including the diagonal to 0
    for icol=1:length(nodenode_adj.colptr)-1
        for irow=nodenode_adj.colptr[icol]:nodenode_adj.colptr[icol+1]-1
            if nodenode_adj.rowval[irow]>=icol
                nodenode_adj.nzval[irow]=0
            end
        end
    end
    dropzeros!(nodenode_adj)

    # Now we know the number of bedges and
    nbedges=length(nodenode_adj.nzval)

    bedgenodes   = zeros(Ti, 2, nbedges)
    bfaceedges   = zeros(Ti, num_edges(bgeom), num_bfaces(grid))

    cen        = local_celledgenodes(bgeom)
    for ibface=1:num_bfaces(grid)
        for ibedge=1:num_edges(bgeom)
            n1 = bfacenodes[cen[1,ibedge], ibface]
            n2 = bfacenodes[cen[2,ibedge], ibface]

            # We need to look in nodenod_adj for upper triangular part entries
            # therefore, we need to swap accordingly before looking
           if (n1<n2)
              n0=n1
              n1=n2
              n2=n0;
           end

            for irow = nodenode_adj.colptr[n1]:nodenode_adj.colptr[n1+1]-1
                if nodenode_adj.rowval[irow]==n2
                    # If the coresponding entry has been found, set its
                    # value. Note that this introduces a different edge orientation
                    # compared to the one found locally from cell data
                    bfaceedges[ibedge, ibface] = irow
                    bedgenodes[1, irow]=n1
                    bedgenodes[2, irow]=n2
                end
            end
        end

    end
    grid[BEdgeNodes] = bedgenodes
    grid[BFaceEdges] = bfaceedges
    true
end

#ExtendableGrids.instantiate(grid, ::Type{BEdgeNodes})=prepare_bedges!(grid) && grid[BEdgeNodes]
ExtendableGrids.instantiate(grid, ::Type{BFaceEdges})=prepare_bedges!(grid) && grid[BFaceEdges]
ExtendableGrids.instantiate(grid, ::Type{BFaceCells})=prepare_bfacecells!(grid) && grid[BFaceCells]

normal!(normal,nodes,coord,::Type{Val{1}}) = normal[1]=1.0

function normal!(normal,nodes,coord,::Type{Val{2}})
    normal[1] = -(coord[2,nodes[1]]-coord[2,nodes[2]])
    normal[2] =   coord[1,nodes[1]]-coord[1,nodes[2]]
    d=norm(normal)
    normal[1]/=d
    normal[2]/=d
end

function normal!(normal,nodes,coord,::Type{Val{3}})
    ax=coord[1,nodes[1]]-coord[1,nodes[2]]
    ay=coord[2,nodes[1]]-coord[2,nodes[2]]
    az=coord[3,nodes[1]]-coord[3,nodes[2]]

    bx=coord[1,nodes[1]]-coord[1,nodes[3]]
    by=coord[2,nodes[1]]-coord[2,nodes[3]]
    bz=coord[3,nodes[1]]-coord[3,nodes[3]]

    normal[1]=(ay*bz-by*az)
    normal[2]=(az*bx-bz*ax)
    normal[3]=(ax*by-bx*ay)

    d=norm(normal)
    normal[1]/=d
    normal[2]/=d
    normal[3]/=d
end

function midpoint!(mid,nodes,coord)
    dim=size(coord,1)
    nn=size(nodes,1)
    for i=1:dim
        mid[i]=sum(coord[i,nodes])/nn
    end
end

function adjust!(normal,cmid,bmid)
    d=dot(normal,bmid-cmid)
    if d<0.0
        normal.*= -1
    end
end

function prepare_bfacenormals!(grid)
    bfnodes=grid[BFaceNodes]
    nbf=size(bfnodes,2)
    cellnodes=grid[CellNodes]
    dim=dim_space(grid)
    bfcells=grid[BFaceCells]
    bfnormals=zeros(dim,nbf)
    coord=grid[Coordinates]
    cmid=zeros(dim)
    bmid=zeros(dim)
    for ibf=1:nbf
        icell=bfcells[1,ibf]
        @views normal!(bfnormals[:,ibf],bfnodes[:,ibf],coord,Val{dim})
        @views midpoint!(cmid,cellnodes[:,icell],coord)
        @views midpoint!(bmid,bfnodes[:,ibf],coord)
        @views adjust!(bfnormals[:,ibf],cmid,bmid)
    end
    grid[BFaceNormals]=bfnormals
    true
end

ExtendableGrids.instantiate(grid, ::Type{BFaceNormals})=prepare_bfacenormals!(grid) && grid[BFaceNormals]


text:
export getFGMRESmem,FGMRES_relaxation

mutable struct FGMRESmem{T}
	v_prec              ::Array{T} ## memory for the result of the preconditioner
	Az					::Array{T} ## memory for the result of A*z.
	V    				::Array{T}
	Z					::Array{T}
end

function resetMem(mem::FGMRESmem{T}) where T
mem.v_prec[:].=0.0;
mem.Az[:].=0.0;
mem.V[:].=0.0;
mem.Z[:].=0.0;
end

function getFGMRESmem(n::Int64,T::Type,k::Int64,m::Int64=1)
## Effectively we do not use v_prec and Az.
v_prec = zeros(T,0);
Az = zeros(T,0);
if m==1
	Az = zeros(T,n);
	v_prec = zeros(T,n);
else
	Az = zeros(T,n,m);
	v_prec = zeros(T,n,m);
end
return FGMRESmem{T}(v_prec,Az,zeros(T,n*m,k),zeros(T,n*m,k));
end

function getEmptyFGMRESmem(T::Type)
return FGMRESmem{T}(zeros(T,0),zeros(T,0),zeros(T,0),zeros(T,0));
end
import Base.isempty
function isempty(mem::FGMRESmem{T}) where T
return size(mem.V,1)==0;
end

function FGMRES_relaxation(AT::SparseMatrixCSC,r0::Array{T},x0::Array{T},inner::Int64,prec::Function,TOL::Float64,
                 verbose::Bool,numCores::Int64, mem::FGMRESmem{T} =  getEmptyFGMRESmem(T)) where T
return FGMRES_relaxation(getAfun(AT,zeros(T,size(r0)),numCores),r0,x0,inner,prec,TOL,verbose,numCores,mem);
end

function FGMRES_relaxation(Afun::Function,r0::Array{T},x0::Array{T},inner::Int64,prec::Function,TOL::Float64,verbose::Bool,
                 numCores::Int64, mem::FGMRESmem{T} = getEmptyFGMRESmem(T)) where T

# For multiple right hand sides we apply a trick: we consider one huge block diagonal system for Ax_i = b_i

m = size(r0,2);
n = size(r0,1);
if isempty(mem)
	# println("WARNING: Allocating memory in GMRES: ",n=,"n",", m = ",m)
	mem = getFGMRESmem(n,T,inner,m);
else
    resetMem(mem);
	if size(mem.V,2)!= inner
		error("FGMRES_relaxation: size of Krylov subspace is different than inner");
	end
end

rnorm0 = norm(r0);

H = zeros(T,inner,inner);
xi = zeros(T,inner);
t = zeros(T,inner);

constOne = one(T);
constZero = zero(T);

Z = mem.Z;
AZ = mem.V;

w = zeros(eltype(r0),0);
rnorms = zeros(inner);

for j = 1:inner
	if j==1
		z = prec(r0);
	else
		z = prec(w);
	end

	Z[:,j] .= reshape(z,n*m); # no memory problem with this line....

	w = Afun(z); # w = A'*z;# w = SpMatMul(A,z,w,numCores);

	AZ[:,j] .= reshape(w,n*m);  # no memory problem with this line....

	BLAS.gemv!('C', constOne, AZ, reshape(w,n*m),constZero,t);

	xi[j] = dot(w,r0);

	H[:,j] = t;
	H[j,:] = t';
	H = 0.5*H + 0.5*H';
	t[:] = pinv(H)*xi;

	rnorms[j] = sqrt(abs(real(dot(t,H*t)-2.0*dot(t,xi)+rnorm0^2)));

    if verbose
		if j==1
			resprev = rnorm0;
		else
			resprev = rnorms[j-1];
		end
        println(string("Inner iter ",j,": FGMRES: residual norm: ",rnorms[j],", relres: ",rnorms[j]/rnorm0,", factor: ",rnorms[j]/resprev));
    end
    if rnorms[j] < TOL
        rnorms = rnorms[1:j];
        break;
    end
end

BLAS.gemv!('N', constOne, Z, t,constZero,reshape(w,n*m)); # w = Z*t  #This is the correction that corresponds to the residual.

addVectors(constOne,w,x0);

return x0,rnorms; # w is the correction vector;
end


text:
### A Pluto.jl notebook ###
# v0.14.8

using Markdown
using InteractiveUtils

# ╔═╡ 41201500-ad8c-451d-99a7-ab4942e23d33
begin
    import Pkg
    # activate a clean environment
    Pkg.activate(mktempdir())

    Pkg.add(name="Images")

	using Images
end

# ╔═╡ 11e70fe0-d987-11eb-245c-db81de0329fb
url = "https://i.imgur.com/VGPeJ6s.jpg"

# ╔═╡ 843a6c5c-f382-474b-addb-8941e6752770
download(url, "philip.jpg")

# ╔═╡ 44f16bc6-0ad2-4dfe-80c0-db2ea763befc
philip = load("philip.jpg")

# ╔═╡ e167ad86-4a57-4d0f-acb9-94e050613803
typeof(philip)

# ╔═╡ b8db912c-3478-44d9-8931-40a8cafabde1
size(philip)

# ╔═╡ 95275236-b5ae-41ce-9592-21c8dee273fa
size(philip)[1]

# ╔═╡ 0ef602dc-2b73-4e04-9dc1-2ad2cfda81cd
begin
	(h,w) = size(philip)
	head = philip[h÷4:h÷2, w÷4:w÷2]
	head
end

# ╔═╡ fe7d0ded-f2d6-47c6-a4f4-8e4e55630f90
[head head]

# ╔═╡ 95ed376a-a5e2-4f29-8672-6a2dcf6a7cb1
[
	head reverse(head, dims=2)
	reverse(head, dims=1) reverse(reverse(head, dims=1), dims=2)
]

# ╔═╡ 1836c020-ae36-41b0-b3b2-dd62715b99da
begin
	# This is a comment, copy will create new memory for new variable
	new_phil = copy(head)
	new_phil
end

# ╔═╡ e4747f83-17cd-48ab-89f1-429038d6435b
begin
	new_phil2 = copy(new_phil)
	new_phil2[100:200, 100:200] .= RGBX(0.8, 0.5, 0.1)
	new_phil2
end

# ╔═╡ c0ec7536-5d06-4e0d-9616-a9cfb8fb0d5b
function redify(color)
	return RGB(color.r, 0, 0)
end

# ╔═╡ d15f4003-8e95-4a88-ae6c-d6dcaae14b09
redify.(new_phil)

# ╔═╡ 90d93c97-f629-4c5d-9caf-0169b93dd656

# ╔═╡ 9682f121-e078-4438-9956-ef1b2b116fba
decimate(image, ratio=5) = image[1:ratio:end, 1:ratio:end]

# ╔═╡ 971d69b4-3c92-41a9-966e-18c75af352de
begin
	poor_phil = decimate(philip, 10)
	poor_phil
end

# ╔═╡ Cell order:
# ╠═41201500-ad8c-451d-99a7-ab4942e23d33
# ╠═11e70fe0-d987-11eb-245c-db81de0329fb
# ╠═843a6c5c-f382-474b-addb-8941e6752770
# ╠═44f16bc6-0ad2-4dfe-80c0-db2ea763befc
# ╠═e167ad86-4a57-4d0f-acb9-94e050613803
# ╠═b8db912c-3478-44d9-8931-40a8cafabde1
# ╠═95275236-b5ae-41ce-9592-21c8dee273fa
# ╠═0ef602dc-2b73-4e04-9dc1-2ad2cfda81cd
# ╠═fe7d0ded-f2d6-47c6-a4f4-8e4e55630f90
# ╠═95ed376a-a5e2-4f29-8672-6a2dcf6a7cb1
# ╠═1836c020-ae36-41b0-b3b2-dd62715b99da
# ╠═e4747f83-17cd-48ab-89f1-429038d6435b
# ╠═c0ec7536-5d06-4e0d-9616-a9cfb8fb0d5b
# ╠═d15f4003-8e95-4a88-ae6c-d6dcaae14b09
# ╠═971d69b4-3c92-41a9-966e-18c75af352de
# ╠═90d93c97-f629-4c5d-9caf-0169b93dd656
# ╠═9682f121-e078-4438-9956-ef1b2b116fba

max_stars_repo_name: codewithkaranjeswani/my-pluto-notebooks

max_issues_repo_name: codewithkaranjeswani/my-pluto-notebooks

max_forks_repo_name: codewithkaranjeswani/my-pluto-notebooks


text:
module TestCV
    using GaussianProcesses
    using GaussianProcesses: get_params, set_params!, update_mll!, get_value
    using Test, Random
    using Distributions: Normal, Uniform, logpdf, MvNormal
    import Calculus

    @testset "leave-one-out" begin
        Random.seed!(1)
        n,p = 10,1
        f_star(x::Real) = abs(x-5)*cos(2*x)
        σ_y = 0.8
        X_distr = Uniform(-2,2)
        ϵ_distr = Normal(0,σ_y)
        x = sort(rand(X_distr, n))
        Y = f_star.(x) .+ rand(ϵ_distr,n)
        k = SEIso(0.5, 0.8)
        logNoise = log(σ_y)
        m = MeanLin([1.0])
        gp = GPE(Matrix(x'), Y, m, k, logNoise)
        optimize!(gp; domean=false, kern=true, noise=true)

        μi,σi2 = GaussianProcesses.predict_LOO(gp)

        CV = 0.0
        @testset "predictions" begin
            for i in 1:n
                xV, yV = gp.x[:,[i]], gp.y[i]
                T = [j for j in 1:n if j!=i]
                xT, yT = gp.x[:,T], gp.y[T]
                gpT = GPE(xT, yT, gp.mean, gp.kernel, gp.logNoise)
                pred_i = predict_y(gpT, xV)
                @test pred_i[1][1] ≈ μi[i] atol=1e-5
                @test  pred_i[2][1] ≈ σi2[i] atol=1e-5
                CV += logpdf(Normal(μi[i], √σi2[i]), gp.y[i])
            end
        end
        @testset "CVmetric" begin
            @test CV ≈ GaussianProcesses.logp_LOO(gp)
        end
        @testset "gradient" begin
            target = function(θ)
                θprev = get_params(gp.kernel)
                set_params!(gp.kernel, θ)
                update_mll!(gp)
                CV = GaussianProcesses.logp_LOO(gp)
                set_params!(gp.kernel, θprev) # put it back
                return CV
            end
            grad_numerical = Calculus.gradient(target, get_params(k))
            update_mll!(gp)
            grad_analytical = GaussianProcesses.dlogpdθ_LOO(gp; noise=false, kern=true, domean=false)
            @test grad_numerical ≈ grad_analytical  atol=1e-6
        end
        @testset "logNoise gradient" begin
            target = function(θ)
                θprev = get_params(gp; noise=true, kern=false, domean=false)
                set_params!(gp, θ; noise=true, kern=false, domean=false)
                update_mll!(gp)
                CV = GaussianProcesses.logp_LOO(gp)
                set_params!(gp, θprev; noise=true, kern=false, domean=false) # put it back
                return CV
            end
            grad_numerical = Calculus.gradient(target, Float64[get_value(gp.logNoise)])
            update_mll!(gp)
            grad_analytical = GaussianProcesses.dlogpdθ_LOO(gp; noise=true, kern=false, domean=false)
            @test grad_numerical ≈ grad_analytical  atol=1e-6
        end
    end

    @testset "folds" begin
        Random.seed!(1)
        n,p = 20,1
        folds = [1:5, 6:14, 15:20]
        f_star(x::Real) = abs(x-5)*cos(2*x)
        σ_y = 0.8
        X_distr = Uniform(-2,2)
        ϵ_distr = Normal(0,σ_y)
        x = sort(rand(X_distr, n))
        Y = f_star.(x) .+ rand(ϵ_distr,n)
        k = SEIso(0.5, 0.8)
        logNoise = log(σ_y)
        m = MeanLin([1.0])
        gp = GPE(Matrix(x'), Y, m, k, logNoise)
        optimize!(gp; domean=false, kern=true, noise=true)

        μ,Σ = GaussianProcesses.predict_CVfold(gp, folds)

        CV = 0.0
        @testset "predictions" begin
            for (ifold,V) in enumerate(folds)
                xV, yV = gp.x[:,V], gp.y[V]
                T = setdiff(1:n, V)
                xT, yT = gp.x[:,T], gp.y[T]
                gpT = GPE(xT, yT, gp.mean, gp.kernel, gp.logNoise)
                pred_V = predict_y(gpT, xV; full_cov=true)
                @test pred_V[1] ≈ μ[ifold] atol=1e-5
                @test Matrix(pred_V[2]) ≈ Σ[ifold] atol=1e-5
                CV += logpdf(MvNormal(μ[ifold], pred_V[2]), gp.y[V])
            end
        end
        @testset "CVmetric" begin
            @test CV ≈ GaussianProcesses.logp_CVfold(gp, folds) atol=1e-5
        end
        @testset "gradient" begin
            target = function(θ)
                θprev = get_params(gp.kernel)
                set_params!(gp.kernel, θ)
                update_mll!(gp)
                CV = GaussianProcesses.logp_CVfold(gp, folds)
                set_params!(gp.kernel, θprev) # put it back
                return CV
            end
            grad_numerical = Calculus.gradient(target, get_params(k))
            update_mll!(gp)
            grad_analytical = GaussianProcesses.dlogpdθ_CVfold(gp, folds; noise=false, kern=true, domean=false)
            @test grad_numerical ≈ grad_analytical  atol=1e-6
        end
        @testset "logNoise gradient" begin
            target = function(θ)
                θprev = get_params(gp; noise=true, kern=false, domean=false)
                set_params!(gp, θ; noise=true, kern=false, domean=false)
                update_mll!(gp)
                CV = GaussianProcesses.logp_CVfold(gp, folds)
                set_params!(gp, θprev; noise=true, kern=false, domean=false) # put it back
                return CV
            end
            grad_numerical = Calculus.gradient(target, Float64[get_value(gp.logNoise)])
            update_mll!(gp)
            grad_analytical = GaussianProcesses.dlogpdθ_CVfold(gp, folds; noise=true, kern=false, domean=false)
            @test grad_numerical ≈ grad_analytical  atol=1e-6
        end
    end
end


text:
using SimpleDirectMediaLayer
const SDL2 = SimpleDirectMediaLayer
using Test

include("lib/SDL.jl")
include("lib/SDL_ttf.jl")
include("lib/SDL_mixer.jl")

# Integration tests
@testset "example1" begin
    include(joinpath(@__DIR__,"..","src","examples","example1.jl"))
end


text:
# The color plan below is better in terminal dark mode
# To use these print funcitons:
# using ABG: title, message, warning, item, done
"""
Print a title.
"""
function title(msg::AbstractString)
    printstyled('\n', msg, '\n', bold = true, color = :cyan)
    printstyled(repeat('=', length(msg) + 2), '\n', color=31)
end

"""
Print a subtitle.
"""
function subtitle(msg::AbstractString)
    printstyled('\n', msg, '\n', color = :cyan)
    printstyled(repeat('-', length(msg) + 2), '\n', color=31)
end

"""
Print a message.
"""
function message(msg::AbstractString)
    printstyled('\n', msg, '\n'; color = :light_magenta)
end

"""
Print a warning.
"""
function warning(msg::AbstractString)
    msg = "!!! Warning: " * msg
    printstyled('\n', msg, '\n'; color=229)
end

"""
Print an item.
"""
function item(it::AbstractString)
    printstyled("\n- $it\n"; color=74)
end

"""
Print a `Done` message.
"""
function done(msg::AbstractString = "Done")
    printstyled(" ... $msg\n"; color=40)
end

function test_my_styled_printing()
    for i in 1:256
        printstyled(' ', lpad("$i", 3, '0'), ' ', color=i)
        (i % 16 == 0) && println()
    end
    println()

    title("This is a title using bold, color cyan and 154")
    subtitle("This is a subtitle using color cyan and 154")
    item("Item 1 using color 74")
    item("Item 2")
    message("A description message using color :light_magenta")
    warning("This is a warning using color 229")
    println("An average messages using color default")
    done("Done.  Using color 40. OK")
end


text:
"""
    s = inpoly2(vert, node, [edge=Nothing, atol=0.0, rtol=eps(), outformat=InOnOut{1,0,-1}])

Check all points defined by `vert` are inside, outside or on bounds of
polygon defined by `node` and `edge`.

`vert` and `node` are matrices with the x-coordinates in the first
and the y-coordinates in the second column.

`edge` is an integer matrix of same size as `node`. The first column contains
the index of the starting node, the second column the index of the ending node.
The polygon needs to be closed. It may contain unconnected cycles.

A point is considered on boundary, if its Euclidian distance to any edge is less
than or equal `max(atol, rtol * span)`, where span is the maximum extension of the
polygon in either direction.

If a point is considered on bounds, no further check is made to determine it it is inside or
outside numerically.

Expected effort of the algorithm is `(N+M)*log(M)`, where `M` is the number of points
and `N` is the number of polygon edges.
"""
function inpoly2(vert, node, edge=zeros(Int); atol::T=0.0, rtol::T=NaN, iyperm=nothing, outformat=InOnBit) where T<:AbstractFloat

    rtol = !isnan(rtol) ? rtol : iszero(atol) ? eps(T)^0.85 : zero(T)
    poly = PolygonMesh(node, edge)
    points = PointsInbound(vert)
    npoints = length(points)

    vmin = minimum(points)
    vmax = maximum(points)
    pmin = minimum(poly)
    pmax = maximum(poly)

    lbar = sum(pmax - pmin)
    tol = max(abs(rtol * lbar), abs(atol))

    ac = areacount(poly)
    stat = ac > 1 ? falses(npoints,2,ac) : falses(npoints,2)
    # flip coordinates so expected efford is minimal

    dvert = vmax - vmin
    ix = dvert[1] < dvert[2] ? 1 : 2
    if isnothing(iyperm)
        iyperm = sortperm(points, 3 - ix)
    end

    inpoly2!(points, iyperm, poly, ix, tol, stat)

    convertout(outformat, InOnBit, stat)
end

"""
    inpoly2_mat(vert, node, edge, fTolx, ftoly, stats)

INPOLY2_MAT the local m-code version of the crossing-number
test. Loop over edges; do a binary-search for the first ve-
rtex that intersects with the edge y-range; do crossing-nu-
mber comparisons; break when the local y-range is exceeded.
"""
function inpoly2!(points, iyperm, poly, ix::Integer, veps::T, stat::S) where {N,T<:AbstractFloat,S<:AbstractArray{Bool,N}}

    nvrt = length(points)   # number of points to be checked
    nedg = edgecount(poly)  # number of edges of the polygon mesh
    vepsx = vepsy = veps
    iy = 3 - ix

    #----------------------------------- loop over polygon edges
    for epos = 1:nedg

        inod = edgeindex(poly, epos, 1)  # from
        jnod = edgeindex(poly, epos, 2)  # to
        # swap order of vertices
        if vertex(poly, inod, iy) > vertex(poly, jnod, iy)
            inod, jnod = jnod, inod
        end

        #------------------------------- calc. edge bounding-box
        xone = vertex(poly, inod, ix)
        yone = vertex(poly, inod, iy)
        xtwo = vertex(poly, jnod, ix)
        ytwo = vertex(poly, jnod, iy)

        xmin0 = min(xone, xtwo)
        xmax0 = max(xone, xtwo)
        xmin = xmin0 - vepsx
        xmax = xmax0 + vepsx
        ymin = yone - vepsy
        ymax = ytwo + vepsy

        ydel = ytwo - yone
        xdel = xtwo - xone
        xysq = xdel^2 + ydel^2
        feps = sqrt(xysq) * veps

        # find top points[:,iy] < ymin by binary search
        ilow = searchfirst(points, iy, iyperm, ymin)
        #------------------------------- calc. edge-intersection
        # loop over all points with y ∈ [ymin,ymax)
        for jpos = ilow:nvrt
            jorig = iyperm[jpos]
            ypos = vertex(points, jorig, iy)
            ypos > ymax && break
            xpos = vertex(points, jorig, ix)

            if xpos >= xmin
                if xpos <= xmax
                    #--------- inside extended bounding box of edge
                    mul1 = ydel * (xpos - xone)
                    mul2 = xdel * (ypos - yone)
                    if abs(mul2 - mul1) <= feps
                        #------- distance from line through edge less veps
                        mul3 = xdel * (2xpos-xone-xtwo) + ydel * (2ypos-yone-ytwo)
                        if abs(mul3) <= xysq ||
                           hypot(xpos- xone, ypos - yone) <= veps ||
                           hypot(xpos- xtwo, ypos - ytwo) <= veps
                            # ---- round boundaries around endpoints of edge
                            setonbounds!(poly, stat, jorig, epos)
                        end
                        if mul1 < mul2 && yone <= ypos < ytwo
                            #----- left of line && ypos exact to avoid multiple counting
                            flipio!(poly, stat, jorig, epos)
                        end
                    elseif mul1 < mul2 && yone <= ypos < ytwo
                        #----- left of line && ypos exact to avoid multiple counting
                        flipio!(poly, stat, jorig, epos)
                    end
                end
            else # xpos < xmin - left of bounding box
                if yone <= ypos <  ytwo
                    #----- ypos exact to avoid multiple counting
                    flipio!(poly, stat, jorig, epos)
                end
            end
        end
    end
    stat
end

"""
    search lowest iy coordinate >= ymin
"""
function searchfirst(points::PointsInbound, iy::Integer, iyperm, ymin)
    ilow = 0
    iupp = length(points) + 1
    @inbounds while ilow < iupp - 1
        imid = ilow + (iupp-ilow) >>> 0x01
        if vertex(points, iyperm[imid], iy) < ymin
            ilow = imid
        else
            iupp = imid
        end
    end
    iupp
end

"""
    flipio!(poly, stat, p, epos)

Flip boolean value of in-out-status of point `p` for all areas associated with edge `epos`.
"""
function flipio!(poly::PolygonMesh, stat::AbstractArray{Bool}, i::Integer, j::Integer)
    statop!((stat,j,a) -> begin stat[j,1,a] = !stat[j,1,a]; end, poly, stat, i, j)
end

"""
    setonbounds!(poly, stat, p, epos)

Set boolean value of on-bounds-status of point `p` for all areas associated with edge `epos`.
"""
function setonbounds!(poly::PolygonMesh, stat::AbstractArray{Bool}, i::Integer, j::Integer)
    statop!((stat, j,a) -> begin stat[j,2,a] = true; end, poly, stat, i, j)
end
function statop!(f!::Function, poly::PolygonMesh{A}, stat::AbstractArray{Bool,N}, p::Integer, ed::Integer) where {A,N}
    for k in 1:max(A, 1)
        area = N == 3 ? areaindex(poly, ed, k) : 1
        if area > 0
            f!(stat, p, area)
        end
    end
end


text:
module Constants

include("./read_json.jl")
using JSON
using .ReadJson

export setFileName, boxes, area_of_box, L, W, H, pb_ss_var, max_LW, lin_plus, lin_plus_2, lin_min, lin_min_2, guillotine_ds_var, masses, sum_of_masses, cm_x_interval_1, cm_x_interval_2, cm_y_interval_1, cm_y_interval_2, guillotine_fix_var, mass_density, max_allowed_weight, max_allowed_weight_constant, mult_boxes_ss_var, number_of_boxes, lb_number_of_boxes, mass_volume, resultfile, resultjson, filename

#JSON
filename = "SPECIFY FULL PATH OF EXPERIMENT FILE"
json = ReadJson.get_data_from_json(filename)
boxes = ReadJson.set_boxes(json)
masses = ReadJson.set_masses(json)
L = ReadJson.set_pallet_length(json)
W = ReadJson.set_pallet_width(json)
H = 2
resultfile = "SPECIFY FULL PATH OF TXT RESULT FILE"
resultjson = "SPECIFY FULL PATH OF JSON RESULT FILE"

#VARIABLES

# ====== TESTING START =====

# A box is made up of its [length, width, height]. The last box is the sentinel => Height = 0
#boxes = [[1, 1, 1], [1, 1, 1], [1, 1, 0]]
#boxes = [[1, 1, 1], [1, 1, 1], [1, 1, 1], [2, 2, 0]]
#boxes = [[3, 3, 1], [1, 1, 1], [3, 3, 0]]
#boxes = [[1, 1, 1], [2, 2, 2], [2, 2, 0]]
#boxes = [[10, 5, 7], [7, 8, 3], [3, 3, 3]]

# A box must have a mass (in kg)
#masses = [1, 1, 1, 0]
#masses = [3, 1, 0]

#Euro pallet dimensions: H must be somewhat larger than the total height of the boxes as some error can occur.
#L = 2.1
#W = 2.1
#H = 3.1

# ===== TESTING FINISHED =====

#partial base static stability variable. The one that decides how much static stability there should be.
pb_ss_var = 0.7
mult_boxes_ss_var = 0.05
number_of_boxes = 2

#guillotine_cut_variable. The one that decides how many percentages of H are allowed as a guillotine cut
guillotine_ds_var = 0.6
guillotine_fix_var = 0.05

#load balancing area. Can be split into two more variables to have separate for x and y.
load_balancing_lower = 0.4
load_balancing_upper = 0.6

#load bearing. Relationship between area and load bearing. 2 = 1*1 box able to carry 2 kg.
#CHOOSE THIS OR NEXT ONE IN SIMPLE LB CALC
load_bearing_factor = 1000
max_allowed_weight = [load_bearing_factor * (boxes[i][1] * boxes[i][2]) for i in 1:length(boxes)-1]
lb_number_of_boxes = 3

# In this calc, it doesn't depend on the size of box, but rather on the box walls and how many layers they consist of.
max_allowed_weight_constant = 30

#CONSTANTS

area_of_box = [boxes[i][1]*boxes[i][2] for i in 1:length(boxes)]

cm_x_interval_1 = load_balancing_lower * boxes[length(boxes)][1]
cm_x_interval_2 = load_balancing_upper * boxes[length(boxes)][1]
cm_y_interval_1 = load_balancing_lower * boxes[length(boxes)][2]
cm_y_interval_2 = load_balancing_upper * boxes[length(boxes)][2]

#linearisation constants
max_LW = max(L, W)
lin_plus = [i*(0.5*(L+W))/9 for i in 1:9]
lin_plus_2 = [lin_plus[i]^2 for i in 1:9]
lin_min = [
    -0.5*max_LW,
    -0.7*0.5*max_LW,
    -0.4*0.5*max_LW,
    -0.1*0.5*max_LW,
    0.1*0.5*max_LW,
    0.3*0.5*max_LW,
    0.4*0.5*max_LW,
    0.7*0.5*max_LW,
    0.5*max_LW
    ]
lin_min_2 = [lin_min[i]^2 for i in 1:9]

sum_of_masses = sum(masses)
mass_density = [masses[i]/(boxes[i][1] * boxes[i][2]) for i in 1:length(boxes)]
mass_volume = [(boxes[i][1]*10 * boxes[i][2]*10 * boxes[i][3]*10) for i in 1:length(boxes)]

end


text:
using Test
using Logging
# using Revise
# using Plasm
using LarSurf
# Logging.configure(level==Logging.Debug)

# include("../src/LarSurf.jl")
# include("../src/block.jl")

@testset "Map inner face ID to outer ID" begin
    data_size = [2,3,4]
    # block_size = [2,2,2]

#     fid = 15
#     face_cart, axis = LarSurf.grid_face_id_to_cartesian(data_size, fid)
#     @test face_cart .= [2, 1, 3]
#     @test axis = 1

#     fid = 37
#     @test face_cart .= [1, 1, 2]
#     @test axis = 2
#     fid = 40
#     @test face_cart .= [1, 1, 4]
#     @test axis = 2
#     fid = 60
#     @test face_cart .= [2, 2, 4]
#     @test axis = 2
#     fid = 85
#     @test face_cart .= [1, 3, 5]
#     @test axis = 3
    fid = 83
    face_cart, axis = LarSurf.grid_face_id_to_cartesian(data_size, fid)
    println(face_cart, " ", axis)

#     @test collect(faces) == [1, 13, 22]
    # [1, 13, 2, 14, 3, 15, 5, 17, 6, 18, 7, 19, 37, 45, 38, 46, 39, 47, 69, 72, 74, 77]
end

@testset "one face should have the same id from both neighbooring voxels" begin
    ids1 = LarSurf.get_face_ids_from_cube_in_grid([2,3,4], [2,2,4],true)
    ids2 = LarSurf.get_face_ids_from_cube_in_grid([2,3,4], [2,2,5],false)
    @test ids1[3] == ids2[3]
end


text:
struct pReLU{N,Aθ,Aγ,AΔ,Aη} <: AbstractLayer{N}
    θ::Aθ
    γ::Aγ
    Δ::AΔ
    η::Aη
    function pReLU(θ::AbstractArray, γ::AbstractArray, Δ::AbstractArray, η::AbstractArray)
        @assert size(θ) == size(γ) == size(Δ) == size(η)
        return new{ndims(θ), typeof(θ), typeof(γ), typeof(Δ), typeof(η)}(θ, γ, Δ, η)
    end
end

function pReLU(::Type{T}, n::Int...) where {T}
    θ = zeros(T, n...)
    γ = ones(T, n...)
    Δ = zeros(T, n...)
    η = zeros(T, n...)
    return pReLU(θ, γ, Δ, η)
end

pReLU(n::Int...) = pReLU(Float64, n...)

Base.repeat(l::pReLU, n::Int...) = pReLU(
    repeat(l.θ, n...), repeat(l.γ, n...), repeat(l.Δ, n...), repeat(l.η, n...)
)

energies(layer::pReLU, x::AbstractArray) = energies(dReLU(layer), x)
free_energies(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = free_energies(dReLU(layer), inputs)
sample_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = sample_from_inputs(dReLU(layer), inputs)
mean_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = mean_from_inputs(dReLU(layer), inputs)
var_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = var_from_inputs(dReLU(layer), inputs)
meanvar_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = meanvar_from_inputs(dReLU(layer), inputs)
mode_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = mode_from_inputs(dReLU(layer), inputs)
mean_abs_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = mean_abs_from_inputs(dReLU(layer), inputs)
std_from_inputs(layer::pReLU, inputs::Union{Real,AbstractArray} = 0) = sqrt.(var_from_inputs(layer, inputs))

function ∂free_energies(layer::pReLU, inputs::Union{Real,AbstractArray} = 0)
    drelu = dReLU(layer)

    lp = ReLU( drelu.θp, drelu.γp)
    ln = ReLU(-drelu.θn, drelu.γn)
    Fp = free_energies(lp,  inputs)
    Fn = free_energies(ln, -inputs)
    F = -logaddexp.(-Fp, -Fn)
    pp = exp.(F - Fp)
    pn = exp.(F - Fn)
    μp, νp = meanvar_from_inputs(lp,  inputs)
    μn, νn = meanvar_from_inputs(ln, -inputs)
    μ2p = @. (νp + μp^2) / 2
    μ2n = @. (νn + μn^2) / 2

    abs_γ = abs.(layer.γ)

    ∂θ = @. -(pp * μp - pn * μn)
    ∂γ = @. (pp * μ2p / (1 + layer.η) + pn * μ2n / (1 - layer.η)) .* sign.(layer.γ)
    ∂Δ = @. -pp * μp / (1 + layer.η) - pn * μn / (1 - layer.η)
    ∂η = @. (
        pp * (-abs_γ * μ2p + layer.Δ * μp) / (1 + layer.η)^2 +
        pn * ( abs_γ * μ2n - layer.Δ * μn) / (1 - layer.η)^2
    )
    return (θ = ∂θ, γ = ∂γ, Δ = ∂Δ, η = ∂η)
end

struct pReLUStats{A<:AbstractArray}
    x::A; xp1::A; xn1::A; xp2::A; xn2::A;
    function pReLUStats(layer::AbstractLayer, data::AbstractArray; wts = nothing)
        xp = max.(data, 0)
        xn = min.(data, 0)
        x = batchmean(layer, data; wts)
        xp1 = batchmean(layer, xp; wts)
        xn1 = batchmean(layer, xn; wts)
        xp2 = batchmean(layer, xp.^2; wts)
        xn2 = batchmean(layer, xn.^2; wts)
        return new{typeof(x)}(x, xp1, xn1, xp2, xn2)
    end
end

Base.size(stats::pReLUStats) = size(stats.x)

function ∂energy(layer::pReLU, stats::pReLUStats)
    @assert size(layer) == size(stats)
    ∂γ = @. (stats.xp2 / (1 + layer.η) + stats.xn2 / (1 - layer.η)) / 2
    ∂γ .*= sign.(layer.γ)
    ∂Δ = @. -(stats.xp1 / (1 + layer.η) - stats.xn1 / (1 - layer.η))
    abs_γ = abs.(layer.γ)
    ∂η = @. (
        (-abs_γ * stats.xp2 / 2 + layer.Δ * stats.xp1) / (1 + layer.η)^2 +
        ( abs_γ * stats.xn2 / 2 + layer.Δ * stats.xn1) / (1 - layer.η)^2
    )
    return (θ = -stats.x, γ = ∂γ, Δ = ∂Δ, η = ∂η)
end

suffstats(layer::pReLU, data::AbstractArray; wts = nothing) = pReLUStats(layer, data; wts)


text:
using ClangCompiler
using ClangCompiler.LLVM
using ClangCompiler.LLVM_full_jll

# source file
src = joinpath(@__DIR__, "HowToUseLLJIT.cpp")

# compilation flags
llvm_include_dir = joinpath(LLVM_full_jll.artifact_dir, "include") |> normpath
@assert isdir(llvm_include_dir)

julia_include_dir = joinpath(Sys.BINDIR, "..", "include", "julia") |> normpath
@assert isdir(julia_include_dir)

args = get_compiler_args()
push!(args, "-I$llvm_include_dir")
push!(args, "-I$julia_include_dir")

# create JIT, generate LLVM IR, and call function
jit = LLJIT(;tm=JITTargetMachine())
irgen = IRGenerator(src, args)
cc = CXCompiler(irgen, jit)
link_process_symbols(cc)
compile(cc)

addr = lookup(jit, "lljit_main")
@eval main() = ccall($(pointer(addr)), Cint, ())

main()

# clean up
dispose(cc)


text:
struct HmmScoreTable
        c_cmat::Matrix{Int}   #char level confusion matrix
        c_f1::Vector{Float64} #char level f1 score
        c_p::Vector{Float64}  #char level precision
        c_r::Vector{Float64}  #char level recall

        w_f1::Float64         #word level f1 score
        w_p::Float64          #word level precision
        w_r::Float64          #word level recall

        n::Int                #how many tables combined

end

HmmScoreTable(xs::Vector{Vector{String}}, ys::Vector{Vector{String}}) = mapreduce(p->HmmScoreTable(p...), +, zip(xs, ys))
"""
        x : Standard output
        y : KongYiji output
"""
function HmmScoreTable(x::Vector{<:AbstractString}, y::Vector{<:AbstractString})
        #### char level
        # confusion matrix
        c_cmat, c_f1, c_p, c_r = fill(0, (4, 4)), fill(0., 4), fill(0., 4), fill(0., 4)
        idx = Dict(:B=>1, :M=>2, :E=>3, :S=>4)
        nchr = mapreduce(length, +, x)
        xa = fill(:M, nchr)
        p = 1
        for w in x
                nw = length(w)
                xa[p] = :B
                xa[p + nw - 1] = :E
                if nw == 1 xa[p] = :S end
                p += nw
        end
        ya = fill(:M, nchr)
        p = 1
        for w in y
                nw = length(w)
                ya[p] = :B
                ya[p + nw - 1] = :E
                if nw == 1 ya[p] = :S end
                p += nw
        end
        for i = 1:nchr c_cmat[idx[ya[i]],idx[xa[i]]] += 1 end
        # precision, recall, f1-score
        c_f1, c_p, c_r = fill(0., 4), fill(0., 4), fill(0., 4)
        for i = 1:4
                num, den = c_cmat[i,i], sum(c_cmat[i,:])
                c_p[i] = num == den ? 1. : (0. + num) / den
                num, den = c_cmat[i,i], sum(c_cmat[:,i])
                c_r[i] = num == den ? 1. : (0. + num) / den
                c_f1[i] = f1(c_p[i], c_r[i])
        end

        ##### word level
        r = length(eachmatch(x, y)) - 2
        w_p = r / length(y)
        w_r = r / length(x)
        w_f1 = f1(w_p, w_r)
        return HmmScoreTable(c_cmat, c_f1, c_p, c_r, w_f1, w_p, w_r, 1)
end

function eachmatch(x::Vector{<:AbstractString}, y::Vector{<:AbstractString})
        ix, iy, px, py, nx, ny = 1, 1, 1, 1, length(x), length(y)
        ret = Tuple{Int, Int}[]
        push!(ret, (0,0))
        #todo sizehint?
        while ix <= nx && iy <= ny
                if px == py
                        if x[ix] == y[iy] push!(ret, (ix,iy)) end
                        nxi, nyi = length(x[ix]), length(y[iy])
                        if nxi == nyi px += nxi; py += nyi; ix += 1; iy += 1
                        elseif nxi < nyi px += nxi; ix += 1
                        else py += nyi; iy += 1
                        end
                elseif px < py
                        while px < py && ix <= nx px += length(x[ix]); ix += 1 end
                else
                        while py < px && iy <= ny py += length(y[iy]); iy += 1 end
                end
        end
        push!(ret, (nx + 1,ny + 1))
        return ret
end

+(a::HmmScoreTable, b::HmmScoreTable) = HmmScoreTable(map(x->getfield(a, x) + getfield(b, x), fieldnames(HmmScoreTable))...)

function show(io::IO, tb::HmmScoreTable)
        n = size(tb.c_cmat, 1)
        fm(d) = string(trunc(d * 100, sigdigits=3), "%")
        fm1(d) = trunc(d, sigdigits=5)
        char_f1, char_p, char_r = map(x->x / tb.n, [tb.c_f1, tb.c_p, tb.c_r])
        word_f1, word_p, word_r = map(x->x / tb.n, [tb.w_f1, tb.w_p, tb.w_r])
        mat = Array{Any}(missing, (n + 4, n + 1))
        mat[1:n,1:n] .= tb.c_cmat
        mat[n + 1,1:n] .= sum(tb.c_cmat, dims=1)[1,:]
        mat[1:n,n + 1] .= sum(tb.c_cmat, dims=2)[:,1]
        mat[n + 1,n + 1] = sum(tb.c_cmat)
        mat[n + 2,1:n] .= fm.(char_p); mat[n + 2,n + 1] = ""
        mat[n + 3,1:n] .= fm.(char_r); mat[n + 3,n + 1] = ""
        mat[n + 4,1:n] .= fm1.(char_f1); mat[n + 4,n + 1] = ""
        char_avg_f1, char_avg_p, char_avg_r = map(x->sum(x)/length(x), [char_f1, char_p, char_r])

        utb = UselessTable(mat; cnames=(:B, :M, :E, :S, ""),
                                rnames=(:B, :M, :E, :S, "", :Precision, :Recall, :F1),
                                topleft="O\\S",
                                foots=["Char level avg.F1: $(fm1(char_avg_f1)), avg.precison: $(fm(char_avg_p)), avg.recall: $(fm(char_avg_r))",
                                      "Word level F1: $(fm1(word_f1)), precison: $(fm(word_p)), recall: $(fm(word_r))",
                                ],
                                heads=["KongYiji(1) HMM Score Table $(tb.n) combined"],
        )
        show(io, utb)
end

f1(x, y) = 2 / (1. / x + 1. / y)