text:
import Base: setindex!, sizehint!, empty!, isempty, length, copy, empty,
             getindex, getkey, haskey, iterate, @propagate_inbounds,
             pop!, delete!, get, get!, isbitstype, in, hashindex, isbitsunion,
             isiterable, dict_with_eltype, KeySet, Callable, _tablesz, filter!

# the load factor after which the dictionary `rehash` happens
const ROBIN_DICT_LOAD_FACTOR = 0.70

"""
    RobinDict([itr])

`RobinDict{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 [`RobinDict`](@ref) whose key-value pairs
are taken from 2-tuples `(key,value)` generated by the argument.

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

Alternatively, a sequence of pair arguments may be passed.

```jldoctest
julia> RobinDict("A"=>1, "B"=>2)
RobinDict{String,Int64} with 2 entries:
  "B" => 2
  "A" => 1
```
"""
mutable struct RobinDict{K,V} <: AbstractDict{K,V}
    hashes::Vector{UInt32}
    keys::Array{K,1}
    vals::Array{V,1}
    count::Int
    idxfloor::Int
end

function RobinDict{K, V}() where {K, V}
    n = 16
    RobinDict{K, V}(zeros(UInt32, n), Vector{K}(undef, n), Vector{V}(undef, n), 0, 0)
end

function RobinDict{K, V}(d::RobinDict{K, V}) where {K, V}
    RobinDict{K, V}(copy(d.hashes), copy(d.keys), copy(d.vals), d.count, d.idxfloor)
end

function RobinDict{K,V}(kv) where V where K
    h = RobinDict{K,V}()
    for (k,v) in kv
        h[k] = v
    end
    return h
end
RobinDict{K,V}(p::Pair) where {K,V} = setindex!(RobinDict{K,V}(), p.second, p.first)
function RobinDict{K,V}(ps::Pair...) where V where K
    h = RobinDict{K,V}()
    sizehint!(h, length(ps))
    for p in ps
        h[p.first] = p.second
    end
    return h
end

RobinDict() = RobinDict{Any,Any}()
RobinDict(kv::Tuple{}) = RobinDict()
copy(d::RobinDict) = RobinDict(d)
empty(d::RobinDict, ::Type{K}, ::Type{V}) where {K, V} = RobinDict{K, V}()

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

RobinDict(d::AbstractDict{K, V}) where {K, V} = RobinDict{K, V}(d)

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

hash_key(key) = (hash(key)%UInt32) | 0x80000000
desired_index(hash, sz) = (hash & (sz - 1)) + 1

function calculate_distance(h::RobinDict{K, V}, index) where {K, V}
    @assert isslotfilled(h, index)
    sz = length(h.keys)
    @inbounds index_init = desired_index(h.hashes[index], sz)
    return (index - index_init + sz) & (sz - 1)
end

# insert algorithm
function rh_insert!(h::RobinDict{K, V}, key::K, val::V) where {K, V}
    # table full
    @assert h.count != length(h.keys)

    ckey, cval, chash = key, val, hash_key(key)
    sz = length(h.keys)
    index_init = desired_index(chash, sz)

    index_curr = index_init
    probe_distance = 0
    probe_current = 0
    @inbounds while true
        if (isslotempty(h, index_curr)) || (isslotfilled(h, index_curr) && isequal(h.keys[index_curr], ckey))
            break
        end
        probe_distance = calculate_distance(h, index_curr)

        if probe_current > probe_distance
            h.vals[index_curr], cval = cval, h.vals[index_curr]
            h.keys[index_curr], ckey = ckey, h.keys[index_curr]
            h.hashes[index_curr], chash = chash, h.hashes[index_curr]
            probe_current = probe_distance
        end
        probe_current += 1
        index_curr = (index_curr & (sz - 1)) + 1
    end

    @inbounds if isslotfilled(h, index_curr) && isequal(h.keys[index_curr], ckey)
        h.vals[index_curr] = cval
        return index_curr
    end

    @inbounds if isslotempty(h, index_curr)
        h.count += 1
    end

    @inbounds h.vals[index_curr] = cval
    @inbounds h.keys[index_curr] = ckey
    @inbounds h.hashes[index_curr] = chash

    @assert probe_current >= 0

    if h.idxfloor == 0
        h.idxfloor = index_curr
    else
        h.idxfloor = min(h.idxfloor, index_curr)
    end
    return index_curr
end

function rh_insert_for_rehash!(h_new::RobinDict{K, V}, key::K, val::V, hash::UInt32) where {K, V}
    # table full
    @assert h_new.count != length(h_new.keys)

    ckey, cval, chash = key, val, hash
    sz = length(h_new.keys)
    index_init = desired_index(chash, sz)

    index_curr = index_init
    probe_distance = 0
    probe_current = 0
    @inbounds while true
        if (isslotempty(h_new, index_curr))
            break
        end
        probe_distance = calculate_distance(h_new, index_curr)

        if probe_current > probe_distance
            h_new.vals[index_curr], cval = cval, h_new.vals[index_curr]
            h_new.keys[index_curr], ckey = ckey, h_new.keys[index_curr]
            h_new.hashes[index_curr], chash = chash, h_new.hashes[index_curr]
            probe_current = probe_distance
        end
        probe_current += 1
        index_curr = (index_curr & (sz - 1)) + 1
    end

    @inbounds if isslotempty(h_new, index_curr)
        h_new.count += 1
    end

    @inbounds h_new.vals[index_curr] = cval
    @inbounds h_new.keys[index_curr] = ckey
    @inbounds h_new.hashes[index_curr] = chash

    @assert probe_current >= 0

    if h_new.idxfloor == 0
        h_new.idxfloor = index_curr
    else
        h_new.idxfloor = min(h_new.idxfloor, index_curr)
    end
    return index_curr
end

#rehash! algorithm
function rehash!(h::RobinDict{K,V}, newsz = length(h.keys)) where {K, V}
    oldk = h.keys
    oldv = h.vals
    oldh = h.hashes
    sz = length(oldk)
    newsz = _tablesz(newsz)
    if h.count == 0
        resize!(h.keys, sz)
        resize!(h.vals, sz)
        resize!(h.hashes, newsz)
        fill!(h.hashes, 0)
        h.count = 0
        h.idxfloor = 0
        return h
    end

    h.keys = Vector{K}(undef, newsz)
    h.vals = Vector{V}(undef, newsz)
    h.hashes = zeros(UInt32,newsz)
    h.count = 0
    h.idxfloor = 0

    for i = 1:sz
        @inbounds if oldh[i] != 0
            k = oldk[i]
            v = oldv[i]
            rh_insert_for_rehash!(h, k, v, oldh[i])
        end
    end
    return h
end

function sizehint!(d::RobinDict, 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

@propagate_inbounds isslotfilled(h::RobinDict, index) = (h.hashes[index] != 0)
@propagate_inbounds isslotempty(h::RobinDict, index) = (h.hashes[index] == 0)

function setindex!(h::RobinDict{K,V}, v0, key0) where {K, V}
    key = convert(K, key0)
    isequal(key, key0) || throw(ArgumentError("$key0 is not a valid key for type $K"))
    _setindex!(h, key, v0)
end

function _setindex!(h::RobinDict{K,V}, key::K, v0) where {K, V}
    v = convert(V, v0)
    sz = length(h.keys)
    (h.count > ROBIN_DICT_LOAD_FACTOR * sz) && rehash!(h, sz<<2)
    index = rh_insert!(h, key, v)
    @assert index > 0
    h
end

isempty(d::RobinDict) = (d.count == 0)
length(d::RobinDict) = d.count

"""
    empty!(collection) -> collection

Remove all elements from a `collection`.

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

julia> empty!(A);

julia> A
RobinDict{String,Int64} with 0 entries
```
"""
function empty!(h::RobinDict{K,V}) where {K, V}
    sz = length(h.keys)
    empty!(h.hashes)
    empty!(h.keys)
    empty!(h.vals)
    resize!(h.keys, sz)
    resize!(h.vals, sz)
    resize!(h.hashes, sz)
    h.count = 0
    h.idxfloor = 0
    return h
end

function rh_search(h::RobinDict{K, V}, key::K) where {K, V}
    sz = length(h.keys)
    chash = hash_key(key)
    index = desired_index(chash, sz)
    cdibs = 0
    @inbounds while true
        if isslotempty(h, index)
            return -1
        elseif cdibs > calculate_distance(h, index)
            return -1
        elseif h.hashes[index] == chash && (h.keys[index] === key || isequal(h.keys[index], key))
            return index
        end
        index = (index & (sz - 1)) + 1
    end
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 = RobinDict("a"=>1, "b"=>2, "c"=>3);

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

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

julia> d
RobinDict{String,Int64} with 4 entries:
  "c" => 3
  "b" => 2
  "a" => 1
  "d" => 4
```
"""
get!(collection, key, default)

get!(h::RobinDict{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
```
"""
get!(f::Function, collection, key)

function get!(default::Callable, h::RobinDict{K,V}, key0::K) where {K, V}
    key = convert(K, key0)
    return _get!(default, h, key)
end

function _get!(default::Callable, h::RobinDict{K,V}, key::K) where V where K
    index = rh_search(h, key)

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

    v = convert(V, default())
    rh_insert!(h, key, v)
    return v
end

function getindex(h::RobinDict{K, V}, key0) where {K, V}
    key = convert(K, key0)
    index = rh_search(h, key)
    @inbounds return (index < 0) ? throw(KeyError(key)) : h.vals[index]
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 = RobinDict("a"=>1, "b"=>2);

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

julia> get(d, "c", 3)
3
```
"""
get(collection, key, default)

function get(h::RobinDict{K,V}, key0, default) where {K, V}
    key = convert(K, key0)
    index = rh_search(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
```
"""
get(::Function, collection, key)

function get(default::Callable, h::RobinDict{K,V}, key0) where {K, V}
    key = convert(K, key0)
    index = rh_search(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 = RobinDict('a'=>2, 'b'=>3)
RobinDict{Char,Int64} with 2 entries:
  'a' => 2
  'b' => 3

julia> haskey(D, 'a')
true

julia> haskey(D, 'c')
false
```
"""
haskey(h::RobinDict, key) = (rh_search(h, key) > 0)
in(key, v::KeySet{<:Any, <:RobinDict}) = (rh_search(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 = RobinDict('a'=>2, 'b'=>3)
RobinDict{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 getkey(h::RobinDict{K,V}, key0, default) where {K, V}
    key = convert(K, key0)
    index = rh_search(h, key)
    @inbounds return (index < 0) ? default : h.keys[index]::K
end

# backward shift deletion by not keeping any tombstones
function rh_delete!(h::RobinDict{K, V}, index) where {K, V}
    @assert index > 0

    # this assumes that there is a key/value present in the dictionary at index
    index0 = index
    sz = length(h.keys)
    @inbounds while true
        index0 = (index0 & (sz - 1)) + 1
        if isslotempty(h, index0) || calculate_distance(h, index0) == 0
            break
        end
    end
    #index0 represents the position before which we have to shift backwards

    # the backwards shifting algorithm
    curr = index
    next = (index & (sz - 1)) + 1

    @inbounds while next != index0
        h.vals[curr] = h.vals[next]
        h.keys[curr] = h.keys[next]
        h.hashes[curr] = h.hashes[next]
        curr = next
        next = (next & (sz-1)) + 1
    end

    #curr is at the last position, reset back to normal
    isbitstype(K) || isbitsunion(K) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), h.keys, curr-1)
    isbitstype(V) || isbitsunion(V) || ccall(:jl_arrayunset, Cvoid, (Any, UInt), h.vals, curr-1)
    @inbounds h.hashes[curr] = 0x0

    h.count -= 1
    # this is necessary because key at idxfloor might get deleted
    h.idxfloor = get_next_filled(h, h.idxfloor)
    return h
end

function _pop!(h::RobinDict, index)
    @inbounds val = h.vals[index]
    rh_delete!(h, index)
    return val
end

function pop!(h::RobinDict{K, V}, key0) where {K, V}
    key = convert(K, key0)
    index = rh_search(h, key)
    return index > 0 ? _pop!(h, index) : throw(KeyError(key))
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 = RobinDict("a"=>1, "b"=>2, "c"=>3);

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

julia> pop!(d, "d")
ERROR: KeyError: key "d" not found
Stacktrace:
[...]

julia> pop!(d, "e", 4)
4
```
"""
pop!(collection, key, default)

function pop!(h::RobinDict{K, V}, key0, default) where {K, V}
    key = convert(K, key0)
    index = rh_search(h, key)
    return index > 0 ? _pop!(h, index) : default
end

function pop!(h::RobinDict)
    isempty(h) && throw(ArgumentError("dict must be non-empty"))
    idx = h.idxfloor
    @inbounds key = h.keys[idx]
    @inbounds val = h.vals[idx]
    rh_delete!(h, idx)
    key => val
end

"""
    delete!(collection, key)

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

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

julia> delete!(d, "b")
RobinDict{String,Int64} with 1 entry:
  "a" => 1
```
"""
function delete!(h::RobinDict{K, V}, key0) where {K, V}
    key = convert(K, key0)
    index = rh_search(h, key)
    if index > 0
        rh_delete!(h, index)
    end
    return h
end

function get_idxfloor(h::RobinDict)
    @inbounds for i = 1:length(h.keys)
        if isslotfilled(h, i)
            return i
        end
    end
    return 0
end

function get_next_filled(h::RobinDict, i)
    L = length(h.keys)
    (1 <= i <= L) || return 0
    for j = i:L
        @inbounds if isslotfilled(h, j)
            return  j
        end
    end
    return 0
end

@propagate_inbounds _iterate(t::RobinDict{K,V}, i) where {K,V} = i == 0 ? nothing : (Pair{K,V}(t.keys[i],t.vals[i]), i == typemax(Int) ? 0 : get_next_filled(t, i+1))
@propagate_inbounds function iterate(t::RobinDict)
    _iterate(t, t.idxfloor)
end
@propagate_inbounds iterate(t::RobinDict, i) = _iterate(t, get_next_filled(t, i))

filter!(f, d::RobinDict) = Base.filter_in_one_pass!(f, d)

ext: jl


text:
using SavitzkyGolay

t = LinRange(-4, 4, 500)
y1 = collect(0:20)
y2 = exp.(-t.^2) .+ 0.05 .* (1.0 .+ randn(length(t)))

# Plain call
y1_ = savitzky_golay(y1, 11, 2)
y2_ = savitzky_golay(y2, 21, 4)

# With constructor
f1 = SGolay(11, 2)
y1c_ = f1(y1)

f2 = SGolay(21, 4)
y2c_ = f2(y2)


text:
#=
SimpleMVO
=========
    Mean-Variance Optimization Model

Author: Azamat Berdyshev, Shen Wang
Date: 01/30/2016
=#

type SimpleMVO{R<:Real, S<:AbstractString} <: AbstractModel
    sense::Symbol
    vars::Vector{Expr}
    _objective::Expr
    _default_constraints::Vector{Expr}
    constraints::Dict{Symbol,Expr}
    assets::AssetsCollection{R,S}
    objVal::R
    weights::Vector{R}
    status::Symbol

    # Inner constructor
    function SimpleMVO(assets::AssetsCollection{R, S},
                        r_min::R,
                        constraints::Dict{Symbol,Expr}=Dict{Symbol,Expr}(),
                        short_sale::Bool=false)

        n = length(assets)
        Σ = getCovariance(assets)
        μ = getReturns(assets)

        # if no short sale => add corresponding constraint
        if short_sale
            vars = [:(w[1:$n])]
        else
            vars = [:(w[1:$n] >= 0)]
        end

        _objective = :(dot(w,$Σ*w))

        _default_constraints = [:(dot($μ,w) ≥ $r_min),
                                :(dot(ones($n),w) == 1)]

        new(:Min,
            vars,
            _objective,
            _default_constraints,
            constraints,
            assets,
            NaN,
            fill(NaN,n),
            :Unsolved)
    end
end

# Outer constructor
SimpleMVO{R<:Real, S<:AbstractString}(
            assets::AssetsCollection{R, S},
            r_min::R,
            constraints=Dict{Symbol,Expr}()::Dict{Symbol,Expr};
            short_sale=false::Bool) = SimpleMVO{R, S}(
                                                assets,
                                                r_min,
                                                constraints,
                                                short_sale)


text:
module Run
using FromFile: @from

@from "./ArtefactModelling.jl" import ArtefactModelling

# ArtefactModelling.geo_test()
# ArtefactModelling.vis_test()
# ArtefactModelling.csg_test()
ArtefactModelling.prism_test()

end # module


text:
using DelimitedFiles

# Lanternfish pack structure
mutable struct Lanternfish
  days:: Int
  n :: Int
  groups :: Vector{Int}
end

function Lanternfish(fish::Vector{Int})
  n = 0
  groups = zeros(Int, 9)
  map(x -> groups[x+1] += 1 , fish)
  n = sum(groups)
  return Lanternfish(0,n,groups)
end

# One-day reproduction cycle.
function nextDay!(lf::Lanternfish)
  permute!(lf.groups,[2,3,4,5,6,7,8,9,1])
  lf.groups[7] += lf.groups[9]
  lf.n += lf.groups[9]
  lf.days += 1
end

# Advance reproduction cycle several days.
function nextDays!(lf::Lanternfish,nDays::Int)
  for i in 1:nDays
    nextDay!(lf)
  end
end

# Read input
function readInput() :: Vector{Int}
  data = readdlm("./inputs/day6.txt", ',', Int)
  return reshape(data,length(data))
end

println("Problem 1: ")
fish = readInput()
lf = Lanternfish(fish)
nextDays!(lf,80)
println("Number of days: ", lf.days)
println("Number of fish: ", lf.n)
println("Final state: ", lf.groups)

println("Problem 2: ")
fish = readInput()
lf = Lanternfish(fish)
nextDays!(lf,256)
println("Number of days: ", lf.days)
println("Number of fish: ", lf.n)
println("Final state: ", lf.groups)


text:
####################################
# Test de operadores de cruce
#   BLX
####################################
using  Test
include("../../algoritmos-geneticos/operadores_cruce/BLX.jl")

@testset "BLX " begin
    n = 10
    c1 = rand(n)
    c2 = rand(n)
    alpha = 0.3
    h1, h2 = BLX(c1,c2, alpha)
    @test length(h1) == n
    @test length(h2) == n
end


text:
using HeavyChain

chain1 = create_chain()
chain2 = create_chain()

dt = 1/1000
t = 0:dt:20
u = sin.(t)

all_xL_1 = Array{Float64}(length(t))
all_xL_2 = Array{Float64}(length(t))

simulate!(chain1, 0.0, dt)
@time for i = 1:length(t)
    simulate!(chain1, u[i], dt)
    all_xL_1[i] = getendpoint(chain1)
end

simulate2!(chain2, 0.0, dt)
@time for i = 1:length(t)
    simulate2!(chain2, u[i], dt)
    all_xL_2[i] = getendpoint(chain2)
end

using Plots
plot(t, all_xL_1, lab="C")
plot!(t, all_xL_2, lab="julia")

max_stars_repo_name: NikonPic/REINFORCEMENT-JULIA

max_issues_repo_name: NikonPic/REINFORCEMENT-JULIA

max_forks_repo_name: NikonPic/REINFORCEMENT-JULIA


text:
module TestCore

using ContextualMacros
using Test

ContextualMacros.@def @f
ContextualMacros.@def @g

@testset begin
    ex = ContextualMacros.expandwith(
        @__MODULE__,
        quote
            @f 1 2 3
            @g 4 5
        end;
        f = ctx -> QuoteNode((:f, ctx.args)),
        g = ctx -> QuoteNode((:g, ctx.args)),
    )
    nodes = filter(x -> x isa QuoteNode, ex.args)
    @test nodes == [QuoteNode((:f, (1, 2, 3))), QuoteNode((:g, (4, 5)))]
end

end  # module


text:
# Utilities for visualization
w_pw = 256
h_pw = 144

w_lr = 768
h_lr = 432

w_HD = 1920
h_HD = 1080

w_4k = 3840
h_4k = 2160

@doc raw"""
    get_coords(fractal::FractalData)

Get the coordinates for the selected `FractalData` object.
"""
function get_coords(fractal::FractalData)
    return fractal.xmin, fractal.xmax, fractal.ymin, fractal.ymax
end

@doc raw"""
    set_coords(fractal::FractalData, xmin, xmax, ymin, ymax)

Set the coordinates for the desired FractalData object.
"""
function set_coords(fractal::FractalData, xmin, xmax, ymin, ymax)
    fractal.xmin, fractal.xmax, fractal.ymin, fractal.ymax = xmin,
        xmax,
        ymin,
        ymax
    return
end

# functions to move on the fractal
@doc raw"""
    move_center!(fractal::FractalData, nStepsX::Int, nStepsY::Int)

Move the center of a FractalData objoct of n steps in the `x` or `y` direction.
`nStepsX=1` equals to a 1% movement to the right, while `nStepsY=1` equals to a
1% movement upside.
"""
function move_center!(fractal::FractalData, nStepsX::Int, nStepsY::Int)

    if nStepsX != 0
        xc = (fractal.xmax + fractal.xmin) / 2
        width = fractal.xmax - fractal.xmin
        dx = width / 100

        xc += nStepsX * dx
        fractal.xmin = xc - width / 2
        fractal.xmax = xc + width / 2
    end
    if nStepsY != 0
        yc = (fractal.ymax + fractal.ymin) / 2
        height = fractal.ymax - fractal.ymin
        dy = height / 100

        yc += nStepsY * dy
        fractal.ymin = yc - height / 2
        fractal.ymax = yc + height / 2
    end
    return preview_fractal(fractal) #preview changes
end

@doc raw"""
    move_up!(fractal::FractalData, nSteps::Int = 1)

Move nSteps up. `nSteps = 1` equals to a 1% movement up.
"""
function move_up!(fractal::FractalData, nSteps::Int = 1)
    move_center!(fractal::FractalData, 0, nSteps)
end

@doc raw"""
    move_down!(fractal::FractalData, nSteps::Int = 1)

Move nSteps down. `nSteps = 1` equals to a 1% movement down.
"""
function move_down!(fractal::FractalData, nSteps::Int = 1)
    move_center!(fractal::FractalData, 0, -nSteps)
end

@doc raw"""
    move_left!(fractal::FractalData, nSteps::Int = 1)

Move nSteps left. `nSteps = 1` equals to a 1% movement left.
"""
function move_left!(fractal::FractalData, nSteps::Int = 1)
    move_center!(fractal::FractalData, -nSteps, 0)
end

@doc raw"""
    move_right!(fractal::FractalData, nSteps::Int = 1)

Move nSteps right. `nSteps = 1` equals to a 1% movement right.
"""
function move_right!(fractal::FractalData, nSteps::Int = 1)
    move_center!(fractal::FractalData, nSteps, 0)
end

@doc raw"""
    zoom!(fractal::FractalData, zoom_factor::Real)

Zoom at the current position by a factor `zoom_factor`.
"""
function zoom!(fractal::FractalData, zoom_factor::Real)
    xc = (fractal.xmax + fractal.xmin) / 2
    width = fractal.xmax - fractal.xmin
    width /= zoom_factor

    fractal.xmin = xc - width / 2
    fractal.xmax = xc + width / 2

    yc = (fractal.ymax + fractal.ymin) / 2
    height = fractal.ymax - fractal.ymin
    height /= zoom_factor

    fractal.ymin = yc - height / 2
    fractal.ymax = yc + height / 2

    if typeof(fractal.xmax) != BigFloat &&
       2 + ceil(Int, -log10(width / fractal.width)) > 16

        @info "Upgrading fractal to high resolution computations"

        fractal.xmin = BigFloat(fractal.xmin)
        fractal.xmax = BigFloat(fractal.xmax)
        fractal.ymin = BigFloat(fractal.ymin)
        fractal.ymax = BigFloat(fractal.ymax)
    end

    return preview_fractal(fractal) #preview changes
end

# custom colorbars, generate them using https://cssgradient.io/

@doc raw"""
    pumpkin(nRepeat::Int = 1)

Pumpkin colormap. `nRepeat` specifies the number of times the colormap loops over.
"""
function pumpkin(nRepeat::Int = 1)
    return ColorGradient(repeat([:black, :red, :orange], nRepeat))
end

@doc raw"""
    fire_and_ice(nRepeat::Int = 1)

"Fire and Ice" colormap. `nRepeat` specifies the number of times the colormap loops over.
"""
function fire_and_ice(nRepeat::Int = 1)
    cmap = ColorGradient(vcat(cgrad(:inferno).colors, cgrad(:ice).colors))
    return ColorGradient(repeat(cmap.colors, nRepeat))
end

@doc raw"""
    deep_space(nRepeat::Int = 1)

"Deep space" colormap. `nRepeat` specifies the number of times the colormap loops over.
"""
function deep_space(nRepeat::Int = 1)
    return ColorGradient(repeat(["#000000", "#193e7c", "#dde2ff"], nRepeat))
end

@doc raw"""
    alien_space(nRepeat::Int = 1)

"Alien space" colormap. `nRepeat` specifies the number of times the colormap loops over.
"""
function alien_space(nRepeat::Int = 1)
    return ColorGradient(repeat(
        ["#1a072a", "#ff3e24", "#ffa805", "#7b00ff"],
        nRepeat,
    ))
end

@doc raw"""
    funky_rainbow(nRepeat::Int = 1)

"Funky rainbow" colormap. `nRepeat` specifies the number of times the colormap loops over.
"""
function funky_rainbow(nRepeat::Int = 1)
    return ColorGradient(repeat(
        [
         "#833ab4",
         "#fdcd1d",
         "#45d0fc",
         "#fc45ae",
         "#747473",
         "#454cfc",
         "#fc4545",
         "#d27625",
         "#00ff10",
         "#45a8fc",
         "#833ab4",
        ],
        nRepeat,
    ))
end

@doc raw"""
    cycle_cmap(cmap, nRepeat::Int = 1)

Function to make a colormap cyclic.

# Arguments
- `cmap` may be either a `Symbol` or a `ColorGradient`.
- `nRepeat` specifies the number of times the colormap loops over.
"""
function cycle_cmap(cmap::Symbol, nRepeat::Int = 1)
    return ColorGradient(repeat(cgrad(cmap).colors, nRepeat))
end

function cycle_cmap(cmap::ColorGradient, nRepeat::Int = 1)
    return ColorGradient(repeat(cmap.colors, nRepeat))
end

# utility functions to convert from bits to precision digits
"""
    bits_to_digits(bits::Int)

Converts bits to the number of digits of accuracy.
"""
function bits_to_digits(bits::Int)
    return floor(bits * log10(2))
end

"""
    digits_to_bits(digits::Int)

Converts a number of figures to the corresponding number of bits required to achieve such accuracy.
"""
function digits_to_bits(digits::Int)
    return ceil(Int, digits / log10(2))
end

# utilities to display or preview the fractal
@doc raw"""
    display_fractal(fractal::FractalData; scale = :none, filename = :none, offset=0)

    display_fractal(
        fractal::Matrix;
        colormap = cycle_cmap(:inferno, 3),
        background_color = :white,
        scale = :linear,
        filename = :none,
        offset = 0
    )

Function to display a fractal, can either use a matrix or a `FractalData` object.

#Arguments
- `fractal`: can either be a `Matrix` or a `FractalData` object.
- `colormap`: the colormap which needs to be used, for `FractalData` uses the one stored insde the object.
- `background_color`: the color of the backgroud, also applies to transparent zones, such as when `Inf` values are not plotted. For `FractalData` uses the one stored insde the object.
- `scale`: a scale function which can be applied to the fractal image before plotting. If `:linear` is selected no scale function is applied. If `:none` is selected, uses the default scale function stored inside the `FractalData` object.
- `filename`: if speccified, save the resulting image at that location.
- `offset`: value to be summed to the image before plotting it.

# Examples
```jldoctest
julia> cmap1 = Mandelbrot.cycle_cmap(:inferno, 5)
       xmin1 = -1.744453831814658538530
       xmax1 = -1.744449945239591698236
       ymin1 = 0.022017835126305555133
       ymax1 = 0.022020017997233506531

julia> fractal1_data = FractalData(xmin1, xmax1, ymin1, ymax1, width = Mandelbrot.w_4k, height = Mandelbrot.h_4k, colormap = cmap1, maxIter = 1500, scale_function = log)

julia> computeMandelbrot!(fractal1_data)

julia> display_fractal(fractal1_data, filename = "mandelbrot1.png")
```
"""
function display_fractal(
    fractal::Matrix;
    colormap = cycle_cmap(:inferno, 3),
    background_color = :white,
    scale = :linear,
    filename = :none,
    offset = 0
)
    img = deepcopy(fractal) .+ offset
    if scale == :log
        img = log.(img) # normalize image to have nicer colors
    elseif scale == :exp
        img = exp.(img)
    elseif typeof(scale) <: Function
        img = scale.(img)
    end

    plot = heatmap(
        img,
        colorbar = :none,
        color = colormap,
        axis = false,
        size = (size(img)[2], size(img)[1]),
        grid = false,
        framestyle = :none,
        background_color = background_color
    )

    if filename != :none
        savefig(filename)
    end

    return plot
end

# # version using the structure
function display_fractal(fractal::FractalData; scale = :none, filename = :none, offset=0)
    if scale == :none
        display_fractal(
            fractal.fractal,
            colormap = fractal.colormap,
            background_color = fractal.background_color,
            scale = fractal.scale_function,
            filename = filename,
            offset = offset,
        )
    else
        display_fractal(
            fractal.fractal,
            colormap = fractal.colormap,
            background_color = fractal.background_color,
            scale = scale,
            filename = filename,
            offset = offset
        )
    end
end

@doc raw"""
    preview_fractal( fractal_data::FractalData; scale = :linear, use_GPU = false)

Function used to preview the fractal contained in a `FractalData` object before computing it.

See also: [`display_fractal`](@ref)
"""
function preview_fractal(
    fractal_data::FractalData;
    scale = :linear,
    use_GPU = false,
)
    if use_GPU
        pixels = computeMandelbrot_GPU(
            fractal_data.xmin,
            fractal_data.xmax,
            fractal_data.ymin,
            fractal_data.ymax,
            w_pw,
            h_pw,
            fractal_data.maxIter,
            true,
        )
    else
        pixels = computeMandelbrot(
            fractal_data.xmin,
            fractal_data.xmax,
            fractal_data.ymin,
            fractal_data.ymax,
            w_pw,
            h_pw,
            fractal_data.maxIter,
            true,
        )
    end

    return display_fractal(
        pixels,
        colormap = fractal_data.colormap,
        scale = scale,
        filename = :none,
    )
end

# function create_animation_old(
#     coords_stop::NTuple{4,AbstractFloat},
#     coords_start = :auto;
#     width = w_lr,
#     height = h_lr,
#     maxIter = :auto,
#     colormap = cycle_cmap(:inferno, 2),
#     scale = :linear,
#     filename = "mandelbrot_$(rand(UInt))$(rand(UInt)).gif",
#     n_frames = 100,
#     fps = 10,
# )
#     p = Progress(n_frames)
#     update!(p, 0)
#     jj = 0
#
#     xmin_f, xmax_f, ymin_f, ymax_f = coords_stop
#     xc_f = (xmax_f + xmin_f) / 2
#     yc_f = (ymax_f + ymin_f) / 2
#
#     if coords_start == :auto
#         xmin_i = xc_f - 1.5
#         xmax_i = xc_f + 1.5
#         ymin_i = yc_f - 1.2
#         ymax_i = yc_f + 1.2
#     else
#         xmin_i, xmax_i, ymin_i, ymax_i = coords_start
#     end
#
#     xc_i = (xmax_i + xmin_i) / 2
#     yc_i = (ymax_i + ymin_i) / 2
#
#
#     dxc = xc_f - xc_i
#     dyc = yc_f - yc_i
#
#     dxc_arr = range(0, dxc, length = n_frames)
#     dyc_arr = range(0, dyc, length = n_frames)
#
#     dW_i = abs(xmax_i - xmin_i)
#     dW_f = abs(xmax_f - xmin_f)
#     dH_i = abs(ymax_i - ymin_i)
#     dH_f = abs(ymax_f - ymin_f)
#
#     dW_arr = 10 .^ range(log10(dW_i), log10(dW_f), length = n_frames)
#     dH_arr = 10 .^ range(log10(dH_i), log10(dH_f), length = n_frames)
#
#     xc_arr = dxc_arr .+ xc_i
#     yc_arr = dyc_arr .+ yc_i
#
#     xmin_arr = xc_arr .- dW_arr / 2
#     xmax_arr = xc_arr .+ dW_arr / 2
#     ymin_arr = yc_arr .- dH_arr / 2
#     ymax_arr = yc_arr .+ dH_arr / 2
#
#     if maxIter == :auto
#         nIters1 = 50
#         nIters2 = 50
#         nIters3 = 50
#         anim = @animate for i = 1:n_frames
#             nIters_new = 50 +
#                          round(Int, log10(4 / abs(xmax_arr[i] - xmin_arr[i]))^4)
#
#             nIters1 = nIters2
#             nIters2 = nIters3
#             nIters3 = max(nIters3, nIters_new)
#             nIters = ceil(Int, (nIters1 + nIters2 + nIters3) / 3)
#             img = computeMandelbrot(
#                 xmin_arr[i],
#                 xmax_arr[i],
#                 ymin_arr[i],
#                 ymax_arr[i],
#                 width,
#                 height,
#                 nIters,
#                 false,
#             )
#
#             display_fractal(
#                 img,
#                 colormap = colormap,
#                 scale = scale,
#                 filename = :none,
#             )
#
#             jj += 1
#             update!(p, jj)
#         end
#     else
#         anim = @animate for i = 1:n_frames
#             img = computeMandelbrot(
#                 xmin_arr[i],
#                 xmax_arr[i],
#                 ymin_arr[i],
#                 ymax_arr[i],
#                 width,
#                 height,
#                 maxIter,
#                 false,
#             )
#
#             display_fractal(
#                 img,
#                 colormap = colormap,
#                 scale = scale,
#                 filename = :none,
#             )
#
#             jj += 1
#             update!(p, jj)
#         end
#     end
#     gif(anim, filename, fps = fps)
#
# end

# v2

# function create_animation_old_v2(
#     coords_stop::NTuple{4,AbstractFloat},
#     coords_start = :auto;
#     width = w_lr,
#     height = h_lr,
#     maxIter = :auto,
#     colormap = cycle_cmap(:inferno, 2),
#     scale = :linear,
#     filename = "mandelbrot_$(rand(UInt))$(rand(UInt)).gif",
#     n_frames = 100,
#     fps = 10,
# )
#     p = Progress(n_frames)
#     update!(p, 0)
#     jj = 0
#
#     xmin_f, xmax_f, ymin_f, ymax_f = coords_stop
#     xc_f = (xmax_f + xmin_f) / 2
#     yc_f = (ymax_f + ymin_f) / 2
#
#     if coords_start == :auto
#         xmin_i = xc_f - 1.5
#         xmax_i = xc_f + 1.5
#         ymin_i = yc_f - 1.2
#         ymax_i = yc_f + 1.2
#     else
#         xmin_i, xmax_i, ymin_i, ymax_i = coords_start
#     end
#
#     xc_i = (xmax_i + xmin_i) / 2
#     yc_i = (ymax_i + ymin_i) / 2
#
#
#     dxc = xc_f - xc_i
#     dyc = yc_f - yc_i
#
#     dxc_arr = range(0, dxc, length = n_frames)
#     dyc_arr = range(0, dyc, length = n_frames)
#
#     dW_i = abs(xmax_i - xmin_i)
#     dW_f = abs(xmax_f - xmin_f)
#     dH_i = abs(ymax_i - ymin_i)
#     dH_f = abs(ymax_f - ymin_f)
#
#     dW_arr = 10 .^ range(log10(dW_i), log10(dW_f), length = n_frames)
#     dH_arr = 10 .^ range(log10(dH_i), log10(dH_f), length = n_frames)
#
#     xc_arr = dxc_arr .+ xc_i
#     yc_arr = dyc_arr .+ yc_i
#
#     xmin_arr = xc_arr .- dW_arr / 2
#     xmax_arr = xc_arr .+ dW_arr / 2
#     ymin_arr = yc_arr .- dH_arr / 2
#     ymax_arr = yc_arr .+ dH_arr / 2
#
#     images = zeros(n_frames, height, width)
#
#     if maxIter == :auto
#         nIters = 50
#         maxVal = 0
#         anim = @animate for i = 1:n_frames
#             nIters_new = 50 +
#                          round(Int, log10(4 / abs(xmax_arr[i] - xmin_arr[i]))^4)
#
#             nIters = max(nIters, nIters_new)
#             img = computeMandelbrot(
#                 xmin_arr[i],
#                 xmax_arr[i],
#                 ymin_arr[i],
#                 ymax_arr[i],
#                 width,
#                 height,
#                 nIters,
#                 false,
#             )
#
#             maxVal = max(maxVal, maximum(img))
#             img[1, 1] = maxVal
#             img[1, 2] = 0.0
#
#             display_fractal(
#                 img,
#                 colormap = colormap,
#                 scale = scale,
#                 filename = :none,
#             )
#
#             jj += 1
#             update!(p, jj)
#         end
#     else
#         anim = @animate for i = 1:n_frames
#             img = computeMandelbrot(
#                 xmin_arr[i],
#                 xmax_arr[i],
#                 ymin_arr[i],
#                 ymax_arr[i],
#                 width,
#                 height,
#                 maxIter,
#                 false,
#             )
#
#             maxVal = max(maxVal, maximum(img))
#             img[1, 1] = maxVal
#             img[1, 2] = 0.0
#
#             display_fractal(
#                 img,
#                 colormap = colormap,
#                 scale = scale,
#                 filename = :none,
#             )
#
#             jj += 1
#             update!(p, jj)
#         end
#     end
#     gif(anim, filename, fps = fps)
#
# end

#v3

@doc raw"""
    create_animation(
        coords_stop::NTuple{4,AbstractFloat},
        coords_start = :auto;
        width = w_lr,
        height = h_lr,
        maxIter = :auto,
        offset = 0,
        colormap = :inferno,
        cycle_colormap = true,
        scale = :linear,
        filename = "default",
        n_frames = 100,
        fps = 10,
    )

Computes an animated gif from `coords_start` to `coords_stop`.

# Arguments
- coords_stop: a 4-tuple containing (xmin, xmax, ymin, ymax) of the ending point.
- coords_start: a 4-tuple containing (xmin, xmax, ymin, ymax) of the starting point. If `:auto`, start from a zoomed out Mandelbrot fractal
- width: width in pixels
- height: height in pixels
- maxIter: maximum number of iterations, if `:auto` tries to automatically determine it for each frame.
- offset: value to be added to the fractal matrix before `scale` is applied. It is better not to use it with tracendental scale functions.
- colormap: a colormap, either a valid `Symbol` or a `ColorGradient`
- cycle_colormap: choose whether to automatically cycle the colormap.
- scale: scale funciton to be applied to the fractal matrix before it is plotted.
- filename: if set, save the animation at the specified position, else save it at the current working directory with a random name in the format `mandelbrot_*.gif`.
- n_frames: number of frames to be computed.
- fps: frames per secod.

# Examples
```jldoctest
julia> xmin = 0.307567454839614329536
julia> xmax = 0.307567454903142214608
julia> ymin = 0.023304267108419154581
julia> ymax = 0.023304267156089095573

julia> Mandelbrot.create_animation((xmin, xmax, ymin, ymax), n_frames = 500, scale = log10, colormap = Mandelbrot.fire_and_ice())
```

"""
function create_animation(
    coords_stop::NTuple{4,AbstractFloat},
    coords_start = :auto;
    width = w_lr,
    height = h_lr,
    maxIter = :auto,
    offset = 0,
    colormap = :inferno,
    cycle_colormap = true,
    scale = :linear,
    filename = "default",
    n_frames = 100,
    fps = 10,
)
    if filename == "default"
        filename = "mandelbrot_$(rand(UInt))$(rand(UInt)).gif"
    end

    p = Progress(n_frames)
    update!(p, 0)
    jj = 0

    xmin_f, xmax_f, ymin_f, ymax_f = coords_stop
    xc_f = (xmax_f + xmin_f) / 2
    yc_f = (ymax_f + ymin_f) / 2

    if coords_start == :auto
        xmin_i = xc_f - 1.5
        xmax_i = xc_f + 1.5
        ymin_i = yc_f - 1.2
        ymax_i = yc_f + 1.2
    else
        xmin_i, xmax_i, ymin_i, ymax_i = coords_start
    end

    xc_i = (xmax_i + xmin_i) / 2
    yc_i = (ymax_i + ymin_i) / 2

    dxc = xc_f - xc_i
    dyc = yc_f - yc_i

    dxc_arr = range(0, dxc, length = n_frames)
    dyc_arr = range(0, dyc, length = n_frames)

    dW_i = abs(xmax_i - xmin_i)
    dW_f = abs(xmax_f - xmin_f)
    dH_i = abs(ymax_i - ymin_i)
    dH_f = abs(ymax_f - ymin_f)

    dW_arr = 10 .^ range(log10(dW_i), log10(dW_f), length = n_frames)
    dH_arr = 10 .^ range(log10(dH_i), log10(dH_f), length = n_frames)

    xc_arr = dxc_arr .+ xc_i
    yc_arr = dyc_arr .+ yc_i

    xmin_arr = xc_arr .- dW_arr / 2
    xmax_arr = xc_arr .+ dW_arr / 2
    ymin_arr = yc_arr .- dH_arr / 2
    ymax_arr = yc_arr .+ dH_arr / 2

    #compute order of magnitude of the zoom
    if cycle_colormap
        magnitude = dW_i / dW_f
        nRepetitions = ceil(Int, log10(magnitude))
        cmap = cycle_cmap(colormap, nRepetitions)
    else
        cmap = colormap
    end

    images = zeros(height, width, n_frames)

    if maxIter == :auto
        nIters = 50
        @info "Computing frames"
        for i = 1:n_frames
            nIters_new = 50 +
                         round(Int, log10(4 / abs(xmax_arr[i] - xmin_arr[i]))^4)

            nIters = max(nIters, nIters_new)
            images[:, :, i] .= computeMandelbrot(
                xmin_arr[i],
                xmax_arr[i],
                ymin_arr[i],
                ymax_arr[i],
                width,
                height,
                nIters,
                false,
            )

            jj += 1
            update!(p, jj)
        end
        maxAll = maximum(images)
        minAll = minimum(images)

        @info "Creating animation"
        update!(p, 0)
        jj = 0

        anim = @animate for i = 1:n_frames
            img = images[:, :, i]
            img[1, 1] = maxAll
            img[1, 2] = minAll

            img .+= 10 #remove zeros

            display_fractal(
                img,
                colormap = cmap,
                scale = scale,
                filename = :none,
                offset = offset
            )
            jj += 1
            update!(p, jj)
        end

    else
        @info "Computing frames"
        nIters = maxIter
        for i = 1:n_frames
            images[:, :, i] .= computeMandelbrot(
                xmin_arr[i],
                xmax_arr[i],
                ymin_arr[i],
                ymax_arr[i],
                width,
                height,
                nIters,
                false,
            )

            jj += 1
            update!(p, jj)
        end
        maxAll = maximum(images)
        minAll = minimum(images)

        @info "Creating animation"
        update!(p, 0)
        jj = 0

        anim = @animate for i = 1:n_frames
            img = images[:, :, i]
            img[1, 1] = maxAll
            img[1, 2] = minAll

            img .+= 10 #remove zeros

            display_fractal(
                img,
                colormap = cmap,
                scale = scale,
                filename = :none,
                offset = offset
            )
            jj += 1
            update!(p, jj)
        end
    end
    gif(anim, filename, fps = fps)

end


text: @time @testset "Reachability.Discretization" begin include("unit_discretization.jl") end


text:
module PropertyDestructuring

using ExprTools

export @destructure

function _is_property_assignee(ex)
    Meta.isexpr(ex, :tuple) && (Meta.isexpr(ex.args[1], :parameters) || any(_is_property_assignee, ex.args))
end

function _is_destructure_property_assignment(ex)
    Meta.isexpr(ex, :(=)) && return _is_property_assignee(ex.args[1])
    if Meta.isexpr(ex, [:local, :global, :const])
        Meta.isexpr(ex.args[1], :(=)) && return _is_property_assignee(ex.args[1].args[1])
    end
    false
end

function _is_destructure_property_for(ex)
    Meta.isexpr(ex, :for) && (
        _is_destructure_property_assignment(ex.args[1]) ||
        Meta.isexpr(ex.args[1], :block) && any(_is_destructure_property_assignment, ex.args[1].args)
    )
end

function _destructure_property_for_multi(ex)
    blk = _destructure(ex.args[2])
    variables = Dict{Symbol, Vector{Any}}()
    in_clauses = Expr(:block)
    for assignment_src in ex.args[1].args
        new_assignee = _deconstruct_assignees!(assignment_src.args[1], variables)
        push!(in_clauses.args, :($new_assignee = $(assignment_src.args[2])))
    end
    assignments = Expr(:block, (
        _to_assignment_ex(sym, var, :local) for sym in keys(variables) for var in variables[sym]
    )...)
    Expr(:for, in_clauses, quote
        $assignments
        $blk
    end)
end

function _destructure_property_for_complex(ex)
    blk = _destructure(ex.args[2])
    new_assignee, variables = _deconstruct_assignees(ex.args[1].args[1])
    src = _destructure(ex.args[1].args[2])
    assignments = Expr(:block, (
        _to_assignment_ex(sym, var, :local) for sym in keys(variables) for var in variables[sym]
    )...)
    :(for $new_assignee in $src
        $assignments
        $blk
    end)
end

function _destructure_property_for(ex)
    # @assert Meta.isexpr(ex, :for) && _is_destructure_property_assignment(ex.args[1])
    Meta.isexpr(ex.args[1], :block) && return _destructure_property_for_multi(ex)
    Meta.isexpr(ex.args[1].args[1], :tuple) && Meta.isexpr(ex.args[1].args[1].args[1], :parameters) ||
        return _destructure_property_for_complex(ex)
    variables = ex.args[1].args[1].args[1].args
    any(Meta.isexpr(v, :(::)) for v in variables) && return _destructure_property_for_complex(ex)
    blk = _destructure(ex.args[2])
    lh = Expr(:tuple, variables...)
    _el = gensym(:el)
    gvars = Expr(:tuple, (:($_el.$var) for var in variables)...)
    src0 = _destructure(ex.args[1].args[2])
    src = :(($gvars for $_el in $(src0)))
    Expr(:for, :($lh = $src), blk)
end

_deconstruct_assignees!(ex, _variables) = ex
function _deconstruct_assignees!(ex::Expr, variables)
    if Meta.isexpr(ex, :tuple) && Meta.isexpr(ex.args[1], :parameters)
        new_symbol = gensym(string(ex))
        variables[new_symbol] = ex.args[1].args
        return new_symbol
    end
    Expr(ex.head, (_deconstruct_assignees!(arg, variables) for arg in ex.args)...)
end

function _deconstruct_assignees(ex)
    variables = Dict{Symbol, Vector{Any}}()
    new_ex = _deconstruct_assignees!(ex, variables)
    (new_ex, variables)
end

_to_assignment_ex(src, var) = _to_assignment_ex(src, var, nothing)
_to_assignment_ex(src, var, context) = isnothing(context) ?
    :($(var) = getproperty($src, $(QuoteNode(var)))) :
    Expr(context, :($(var) = getproperty($src, $(QuoteNode(var)))))
function _to_assignment_ex(src, ex::Expr, context)
    if !Meta.isexpr(ex, :(::), 2)
        # FALLBACK, maybe cause Error
        return if context === nothing
            :($(ex) = getproperty($src, $(QuoteNode(ex))))
        else
            Expr(context, :($(ex) = getproperty($src, $(QuoteNode(ex)))))
        end
    end
    var, _type = ex.args
    if context === nothing
        :($(var)::$(_type) = getproperty($src, $(QuoteNode(var))))
    else
        Expr(context, :($(var)::$(_type) = getproperty($src, $(QuoteNode(var)))))
    end
end

function _destructure_property_assignment_complex(ex, context)
    # @assert _is_destructure_property_assignment(ex)
    new_assignee, variables = _deconstruct_assignees(ex.args[1])
    src = _destructure(ex.args[2])
    assignments = Expr(:block, (_to_assignment_ex(sym, var, context) for sym in keys(variables) for var in variables[sym])...)
    quote
        $new_assignee = $src
        $assignments
        $src
    end
end

function _destructure_property_assignment(ex, context=nothing)
    # @assert _is_destructure_property_assignment(ex)
    Meta.isexpr(ex, [:local, :global, :const]) && return _destructure_property_assignment(ex.args[1], ex.head)
    Meta.isexpr(ex.args[1], :tuple) && Meta.isexpr(ex.args[1].args[1], :parameters) ||
        return _destructure_property_assignment_complex(ex, context)
    variables = ex.args[1].args[1].args
    src = _destructure(ex.args[2])
    assignments = Expr(:block, (_to_assignment_ex(src, var, context) for var in variables)...)
    quote
        $assignments
        $src
    end
end

# original(inspired by): https://github.com/FluxML/MacroTools.jl/blob/65c55530b63918daac5f041b144c50d4e34e7984/src/utils.jl#L438-L452
splitarg(name::Symbol) = (name, :Any, false, nothing)
function splitarg(arg_expr)
    if Meta.isexpr(arg_expr, :kw)
        # デフォルト値が設定されている場合
        (splitarg(arg_expr.args[1])[1:3]..., arg_expr.args[2])
    elseif Meta.isexpr(arg_expr, :(...))
        # `...` が指定されている場合
        (splitarg(arg_expr.args[1])[1:2]..., true, nothing)
    elseif Meta.isexpr(arg_expr, :(::))
        # 型アノテーションが存在する場合
        _name = length(arg_expr.args) > 1 ? arg_expr.args[1] : nothing
        _type = arg_expr.args[end]
        (_name, _type, false, nothing)
    else
        (arg_expr, :Any, false, nothing)
    end
end

# original(inspired by): https://github.com/FluxML/MacroTools.jl/blob/65c55530b63918daac5f041b144c50d4e34e7984/src/utils.jl#L414-L418
function combinearg(_name, _type, _isslurp, _default)
    ex = isnothing(_name) ? :(::$_type) : :($_name::$_type)
    if _isslurp
        ex = Expr(:(...), ex)
    end
    if !isnothing(_default)
        ex = Expr(:kw, ex, _default)
    end
    ex
end

function _destructure_function_definition(ex)
    # @assert ExprTools.isdef(ex)
    def_dict = ExprTools.splitdef(ex)
    body = def_dict[:body]
    # args = ExprTools.splitarg.(def_dict[:args])
    # if !any(_is_property_assignee, (arg[1] for arg in args))
    assignments = Expr[]
    if haskey(def_dict, :args)
        for (i, arg) in enumerate(def_dict[:args])
            (_name, _type, _slurp, _default) = splitarg(arg)
            if _is_property_assignee(_name)
                newvar = gensym()
                def_dict[:args][i] = combinearg(newvar, _type, _slurp, _default)
                variables = _name.args[1].args
                append!(assignments, [:($(var) = getproperty($newvar, $(QuoteNode(var)))) for var in variables])
            end
        end
    end
    body = Expr(body.head, _destructure.(body.args)...)
    if !isempty(assignments)
        body = Expr(:block, assignments..., body)
    end
    def_dict[:body] = body
    ExprTools.combinedef(def_dict)
end

isdef(ex) = false
isdef(ex::Expr) = !isnothing(ExprTools.splitdef(ex; throw=false))

_destructure(ex) = ex

function _destructure(ex::Expr)
    _is_destructure_property_for(ex) && return _destructure_property_for(ex)
    _is_destructure_property_assignment(ex) && return _destructure_property_assignment(ex)
    isdef(ex) && return _destructure_function_definition(ex)
    Expr(ex.head, _destructure.(ex.args)...)
end

macro destructure(ex::Expr)
    VERSION < v"1.7.0-DEV.364" || return esc(ex)
    esc(_destructure(ex))
end

end


text:
function balance_str(model::MetNet, ider; digits = 50)
    meti = metindex(model, ider)
    rxns = model.rxns[met_rxns(model, ider)] |> sort
    b_str = []
    for rxn in rxns
        s = round(S(model, meti, rxn), digits = digits)
        push!(b_str, "($s)$rxn")
    end
    return "$(model.mets[meti]): " * join(b_str, " + ") * " == " *
                "$(round(model.b[meti], digits = digits))"
end

function balance_str(model::MetNet, state::AbstractMetState, ider; digits = 50)
    meti = metindex(model, ider)
    rxns = model.rxns[met_rxns(model, ider)] |> sort
    b_str = []
    b = 0.0
    for rxn in rxns
        s = round(S(model, meti, rxn), digits = digits)
        f = round(av(model, state, rxn), digits = digits)
        sf = round(s*f, digits = digits)
        b += s*f
        push!(b_str, "($s*$f=$sf)$(rxn)")
    end
    return "$(model.mets[meti]): " * join(b_str, " + ") * " [tot: $(round(b, digits = digits))]" * " == " *
                "$(round(model.b[meti], digits = digits))"
end
balance_str(model, state::AbstractMetState, iders::Vector; digits = 50) =
    [balance_str(model, state, ider; digits = digits) for ider in iders]


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

using Markdown
using InteractiveUtils

# ╔═╡ 7818dfb7-1574-4989-826d-824f8b1397d3
begin
    import Pkg
    Pkg.activate(Base.current_project())

	using QuantumOptics
end

# ╔═╡ de09dfbf-d4ae-400e-b1ae-d140c02d5857
begin
	using PyPlot
	ion()
end

# ╔═╡ b960f032-f0d2-11eb-018b-09493fde5a0e
md"## Intrinsic decoherence dynamics in the three-coupled harmonic oscillators interaction"

# ╔═╡ e989ff96-be84-4a14-8688-1a488aba28db
begin
	const dim = 30

	fbasis = FockBasis(dim)
end

# ╔═╡ 34d5d453-f526-4353-bff9-1212b1649804
begin
	a = destroy(fbasis)
	n = number(fbasis)
end

# ╔═╡ 3b135ff8-38a2-43ab-b5bc-21bf3ba58d2d
function aT(pos)
	if pos == 1
		return tensor(a, identityoperator(fbasis), identityoperator(fbasis))
	elseif pos == 2
		return tensor(identityoperator(fbasis), a, identityoperator(fbasis))
	elseif pos == 3
		return tensor(identityoperator(fbasis), identityoperator(fbasis), a)
	end
end

# ╔═╡ 6e4d9b8a-c987-4d23-a80e-dfbddedcb3e2
function nT(pos)
	if pos == 1
		return tensor(n, identityoperator(fbasis), identityoperator(fbasis))
	elseif pos == 2
		return tensor(identityoperator(fbasis), n, identityoperator(fbasis))
	elseif pos == 3
		return tensor(identityoperator(fbasis), identityoperator(fbasis), n)
	end
end

# ╔═╡ bd5de8ee-7b51-46e3-98c1-95acb088ab6b
Hfree(ω) = ω * (nT(1) + nT(2) + nT(3))

# ╔═╡ 0a4f0325-d2c1-4cb5-8eef-a66eee7e2473
Hl(λ) = λ * (dagger(aT(1)) * aT(2) + dagger(aT(2)) * aT(1))

# ╔═╡ 379f0235-dc11-4a2f-a182-95782b4b5c17
Hg(g) = g * (aT(3) * (dagger(aT(1)) + dagger(aT(2))) + dagger(aT(3)) * (aT(1) + aT(2)))

# ╔═╡ 32dc8f18-b5c4-42d1-a86c-658f2f1270a1
H(ω, λ, g) = Hfree(ω) + Hl(λ) + Hg(g)

# ╔═╡ cd4e925f-52c1-4a2c-93c1-48fc4d48aa16
ψ₀(α) = tensor(coherentstate(fbasis, α), fockstate(fbasis, 0), fockstate(fbasis, 0))

# ╔═╡ fd1236b5-2b44-4a85-bb3b-787ad2bc3dab
tlist = [0:0.1:50;]

# ╔═╡ f01cf97c-a18d-448e-ac57-c923cea6dc73
tout, ψₜ = timeevolution.schroedinger(tlist, ψ₀(4), H(4, 0.5, 0.1))

# ╔═╡ 5d81d922-a86e-4937-a7e8-7caba47a0bdb
begin
	expn1 = expect(nT(1), ψₜ)
	expn2 = expect(nT(2), ψₜ)
	expn3 = expect(nT(3), ψₜ)
end

# ╔═╡ 63519656-c468-467e-a37c-3eadc730425a
begin
	plot(tout, expn1)
	plot(tout, expn2)
	plot(tout, expn3)
end

# ╔═╡ Cell order:
# ╟─b960f032-f0d2-11eb-018b-09493fde5a0e
# ╠═7818dfb7-1574-4989-826d-824f8b1397d3
# ╠═e989ff96-be84-4a14-8688-1a488aba28db
# ╠═34d5d453-f526-4353-bff9-1212b1649804
# ╠═3b135ff8-38a2-43ab-b5bc-21bf3ba58d2d
# ╠═6e4d9b8a-c987-4d23-a80e-dfbddedcb3e2
# ╠═bd5de8ee-7b51-46e3-98c1-95acb088ab6b
# ╠═0a4f0325-d2c1-4cb5-8eef-a66eee7e2473
# ╠═379f0235-dc11-4a2f-a182-95782b4b5c17
# ╠═32dc8f18-b5c4-42d1-a86c-658f2f1270a1
# ╠═cd4e925f-52c1-4a2c-93c1-48fc4d48aa16
# ╠═fd1236b5-2b44-4a85-bb3b-787ad2bc3dab
# ╠═f01cf97c-a18d-448e-ac57-c923cea6dc73
# ╠═5d81d922-a86e-4937-a7e8-7caba47a0bdb
# ╠═de09dfbf-d4ae-400e-b1ae-d140c02d5857
# ╠═63519656-c468-467e-a37c-3eadc730425a

max_stars_repo_name: rurz/IntrinsicDecoherence

max_issues_repo_name: rurz/IntrinsicDecoherence

max_forks_repo_name: rurz/IntrinsicDecoherence


text:
"""
`PrincipalAxes` is an abstract type to be used
as the argument of a function.
"""
abstract PrincipalAxes


text:
using Test
using Classes
using Suppressor
using MacroTools: striplines

@test superclass(Class) === nothing

@test isclass(AbstractClass) == false       # not concrete
@test isclass(Int) == false                 # not <: AbstractClass

@class Foo <: Class begin
   foo::Int

   Foo() = Foo(0)

   # Although Foo is immutable, subclasses might not be,
   # so it's still useful to define this method.
   function Foo(self::absclass(Foo))
        self.foo = 0
    end
end

@test classof(AbstractFoo) == Foo
@test classof(Foo) == Foo

@test superclass(Foo) == Class
@test_throws Exception superclass(AbstractFoo)

function clean_str(s)
    s = replace(s, r"\n" => " ")
    s = replace(s, r"\s\s+" => " ")
end

@class mutable Bar <: Foo begin
    bar::Int

    # Mutable classes can use this pattern
    function Bar(self::Union{Nothing, absclass(Bar)}=nothing)
        self = (self === nothing ? new() : self)
        superclass(Bar)(self)
        Bar(self, 0)
    end
end

@class mutable Baz <: Bar begin
   baz::Int

   function Baz(self::Union{Nothing, absclass(Baz)}=nothing)
        self = (self === nothing ? new() : self)
        superclass(Baz)(self)
        Baz(self, 0)
    end
end

# emitted by class Foo
# @method foo(obj::Foo) = obj.foo

@method function sum(obj::Bar)
    return obj.foo + obj.bar
end

@test Set(superclasses(Baz)) == Set([Foo, Bar, Class])
@test Set(subclasses(Foo)) == Set(Any[Bar, Baz])

x = Foo(1)
y = Bar(10, 11)
z = Baz(100, 101, 102)

@test fieldnames(Foo) == (:foo,)
@test fieldnames(Bar) == (:foo, :bar)
@test fieldnames(Baz) == (:foo, :bar, :baz)

@method foo(x::Foo) = x.foo

@test foo(x) == 1
@test foo(y) == 10
@test foo(z) == 100

@test sum(y) == 21

@method get_bar(x::Bar) = x.bar

@test get_bar(y) == 11
@test get_bar(z) == 101

@test_throws Exception get_bar(x)

# Mutable
@method set_foo!(x::Foo, value) = (x.foo = value)
set_foo!(z, 1000)
@test foo(z) == 1000

# Immutable
@test_throws Exception x.foo = 1000

# test that where clause is amended properly
@method zzz(obj::Foo, bar::T) where {T} = T

@test zzz(x, :x) == Symbol
@test zzz(y, 10.6) == Float64

# Test other @class structural errors
@test_throws(LoadError, eval(Meta.parse("@class X2 x y")))
@test_throws(LoadError, eval(Meta.parse("@class (:junk,)")))

# Test that classof fails if an abstract class has multiple concrete classes
@class Blink <: Baz
struct RenegadeStruct <: AbstractBlink end

@test_throws Exception classof(AbstractBlink)

# Test that parameterized type is handled properly
@class TupleHolder{NT <: NamedTuple} begin
    nt::NT
end

nt = (foo=1, bar=2)
NT = typeof(nt)
th = TupleHolder{NT}(nt)

@test typeof(th).parameters[1] == NT
@test th.nt.foo == 1
@test th.nt.bar == 2

# Test updating using instance of parent class
bar = Bar(1, 2)
baz = Baz(100, 101, 102)

upd = Baz(bar, 555)
@test upd.foo == 1 && upd.bar == 2 && upd.baz == 555

# ...and with parameterized types
@class mutable SubTupleHolder{NT <: NamedTuple} <: Baz begin
    nt::NT
end

sub = SubTupleHolder{NT}(z, nt)
@test sub.nt.foo == 1 && sub.nt.bar == 2 && sub.foo == 1000 && sub.bar == 101 && sub.baz == 102

xyz = SubTupleHolder(sub, 10, 20, 30, (foo=111, bar=222))
@test xyz.nt.foo == 111 && xyz.nt.bar == 222 && xyz.foo == 10

# Test method structure)
# "First argument of method whatever must be explicitly typed"
@test_throws(LoadError, eval(Meta.parse("@method whatever(i) = i")))

expr = quote
    abstract type AbstractX <: AbstractClass end
    struct X{} <: AbstractX
        function X()
            #= /Users/rjp/.julia/dev/Classes/src/Classes.jl:136 =#
            new()
        end
        function X(_self::T) where T <: AbstractX
            _self
        end
    end
    Classes.superclass(::Type{X}) = begin
            Class
        end
    X
end

expected1 = striplines(expr)
emitted1  = striplines(Classes._defclass(:X, Class, false, nothing, []))

@test string(expected1) == string(emitted1)

expr = quote
    abstract type AbstractX <: AbstractClass end
    mutable struct X{NT <: NamedTuple} <: AbstractX
        i::Int
        j::Int
        function X{NT}(i::Int, j::Int) where NT <: NamedTuple
            new{NT}(i, j)
        end
        function X(_self::T, i::Int, j::Int) where {T <: AbstractX, NT <: NamedTuple}
            _self.i = i
            _self.j = j
            _self
        end
    end
    Classes.superclass(::Type{X}) = begin
            Class
        end
    X
end

expected2 = striplines(expr)
emitted2  = striplines(Classes._defclass(:X, Class, true, [:(NT <: NamedTuple)], [:(i::Int), :(j::Int)]))

@test string(expected2) == string(emitted2)


text:
# This README was generated directly from
# [this source file](https://github.com/fredrikekre/Literate.jl/blob/master/examples/README.jl)
# running these commands from the package root of Literate.jl:


text:
module Sim
include("./name2idx/parameters.jl")
include("./name2idx/species.jl")
include("./set_model.jl")
include("./observable.jl")

using .C
using .V

using Sundials
using SteadyStateDiffEq

# Options for ODE solver
const ABSTOL = 1e-9
const RELTOL = 1e-9
const DTMIN = 1e-8

const normalization = true
#=
if true, simulation results in each observable
are divided by their maximum values
=#
const dt = 1.0
t = collect(0.0:dt:5400.0)  # 0, 1, 2, ..., 5400 [sec.]

const conditions = ["EGF", "HRG"]

simulations = Array{Float64,3}(
    undef, length(observables), length(t), length(conditions)
)

function solveode(
        f::Function,
        u0::Vector{Float64},
        t::Vector{Float64},
        p::Vector{Float64})::Union{ODESolution{}, Nothing}
    local sol::ODESolution{}, is_successful::Bool
    try
        prob = ODEProblem(f,u0,(t[1],t[end]),p)
        sol = solve(
            prob,CVODE_BDF(),
            abstol=ABSTOL,reltol=RELTOL,dtmin=DTMIN,saveat=dt,verbose=false
        )
        is_successful = ifelse(sol.retcode === :Success, true, false)
    catch
        is_successful = false
    finally
        if !is_successful
            GC.gc()
        end
    end
    return is_successful ? sol : nothing
end

function get_steady_state(
        f::Function,
        u0::Vector{Float64},
        p::Vector{Float64})::Vector{Float64}
    local sol::SteadyStateSolution{}, is_successful::Bool
    try
        prob = ODEProblem(diffeq,u0,(0.0,Inf),p)
        prob = SteadyStateProblem(prob)
        sol = solve(
            prob,
            DynamicSS(
                CVODE_BDF();abstol=ABSTOL,reltol=RELTOL
            ),
            dt=dt,dtmin=DTMIN,verbose=false
        )
        is_successful = ifelse(sol.retcode === :Success, true, false)
    catch
        is_successful = false
    finally
        if !is_successful
            GC.gc()
        end
    end
    return is_successful ? sol.u : []
end

function simulate!(p::Vector{Float64}, u0::Vector{Float64})::Union{Bool, Nothing}
    # get steady state
    p[C.Ligand] = p[C.no_ligand]
    u0 = get_steady_state(diffeq,u0,p)
    if isempty(u0)
        return false
    end
    # add ligand
    for (i,condition) in enumerate(conditions)
        if condition == "EGF"
            p[C.Ligand] = p[C.EGF]
        elseif condition == "HRG"
            p[C.Ligand] = p[C.HRG]
        end
        sol = solveode(diffeq,u0,t,p)
        if sol === nothing
            return false
        else
            @inbounds @simd for j in eachindex(t)
                simulations[observables_index("Phosphorylated_MEKc"),j,i] = (
                    sol.u[j][V.ppMEKc]
                )
                simulations[observables_index("Phosphorylated_ERKc"),j,i] = (
                    sol.u[j][V.pERKc] + sol.u[j][V.ppERKc]
                )
                simulations[observables_index("Phosphorylated_RSKw"),j,i] = (
                    sol.u[j][V.pRSKc] + sol.u[j][V.pRSKn]*(p[C.Vn]/p[C.Vc])
                )
                simulations[observables_index("Phosphorylated_CREBw"),j,i] = (
                    sol.u[j][V.pCREBn]*(p[C.Vn]/p[C.Vc])
                )
                simulations[observables_index("dusp_mRNA"),j,i] = (
                    sol.u[j][V.duspmRNAc]
                )
                simulations[observables_index("cfos_mRNA"),j,i] = (
                    sol.u[j][V.cfosmRNAc]
                )
                simulations[observables_index("cFos_Protein"),j,i] = (
                    (sol.u[j][V.pcFOSn] + sol.u[j][V.cFOSn])*(p[C.Vn]/p[C.Vc])
                    + sol.u[j][V.cFOSc] + sol.u[j][V.pcFOSc]
                )
                simulations[observables_index("Phosphorylated_cFos"),j,i] = (
                    sol.u[j][V.pcFOSn]*(p[C.Vn]/p[C.Vc]) + sol.u[j][V.pcFOSc]
                )
            end
        end
    end
end
end # module

max_stars_repo_name: himoto/paramEstim

max_issues_repo_name: himoto/paramEstim

max_forks_repo_name: himoto/paramEstim


text:
#the statmech Vumps
#it made sense to seperate both vumpses as
# - leading_boundary primarily works on MPSMultiline
# - they search for different eigenvalues
# - ham vumps should use Lanczos, this has to use arnoldi
# - this vumps updates entire collumns (state[:,i]); incompatible with InfiniteMPS
# - (a)c-prime takes a different number of arguments
# - it's very litle duplicate code, but together it'd be a bit more convoluted (primarily because of the indexing way)

"
    leading_boundary(state,opp,alg,envs=environments(state,ham))

    approximate the leading eigenvector for opp
"
function leading_boundary(state::InfiniteMPS,H,alg,envs=environments(state,H))
    (st,pr,de) = leading_boundary(convert(MPSMultiline,state),Multiline([H]),alg,envs)
    return convert(InfiniteMPS,st),pr,de
end

function leading_boundary(state::MPSMultiline, H,alg::Vumps,envs = environments(state,H))

    galerkin  = 1+alg.tol_galerkin
    iter       = 1

    temp_ACs = similar.(state.AC);
    temp_Cs = similar.(state.CR);

    while true

        eigalg = Arnoldi(tol=alg.tol_galerkin/10)

        @sync for col in 1:size(state,2)

            @Threads.spawn begin
                (vals_ac,vecs_ac) = eigsolve(RecursiveVec($state.AC[:,col]), 1, :LM, eigalg) do x
                    y = similar(x.vecs);

                    @sync for i in 1:length(x)
                        @Threads.spawn y[mod1(i+1,end)] = ac_prime(x[i],i,$col,$state,$envs)
                    end

                    RecursiveVec(y)
                end

                $temp_ACs[:,col] = vecs_ac[1][:]
            end

            @Threads.spawn begin
                (vals_c,vecs_c) = eigsolve(RecursiveVec($state.CR[:,col]), 1, :LM, eigalg) do x
                    y = similar(x.vecs);

                    @sync for i in 1:length(x)
                        @Threads.spawn y[mod1(i+1,end)] = c_prime(x[i],i,$col,$state,$envs)
                    end

                    RecursiveVec(y)
                end
                $temp_Cs[:,col] = vecs_c[1][:]
            end
        end

        for row in 1:size(state,1),col in 1:size(state,2)
            QAc,_ = leftorth!(temp_ACs[row,col], alg=TensorKit.QRpos())
            Qc,_  = leftorth!(temp_Cs[row,col], alg=TensorKit.QRpos())
            temp_ACs[row,col] = QAc*adjoint(Qc)
        end

        state = MPSMultiline(temp_ACs,state.CR[:,end]; tol = alg.tol_gauge, maxiter = alg.orthmaxiter)
        recalculate!(envs,state);

        (state,envs) = alg.finalize(iter,state,H,envs) :: Tuple{typeof(state),typeof(envs)};

        galerkin = calc_galerkin(state, envs)
        alg.verbose && @info "vumps @iteration $(iter) galerkin = $(galerkin)"

        if (galerkin <= alg.tol_galerkin) || iter>=alg.maxiter
            iter>=alg.maxiter && @warn "vumps didn't converge $(galerkin)"
            return state, envs, galerkin
        end

        iter += 1
    end
end


text:
@export function getshaderiv(shader::Cuint, name::GLenum)
   param::Cint = 0
   @c glGetShaderiv(shader, name, &param)
   return param
end

function getshaderiv(shader::Cuint, name::Symbol)
   @assert name in keys(ShaderParamMap) "Unrecognized shader parameter: $name"
   param = ShaderParamMap[name]
   return getshaderiv(shader, param)
end

@export shadertype(shader::Cuint)  = getshaderiv(shader, :type)
@export iscompiled(shader::Cuint)   = getshaderiv(shader, :compile_status) == GL_TRUE

@export function shaderinfo(shader::Cuint)
   len::Cint = 0
   maxlen = getshaderiv(shader, :info_log_length)
   message = StringArray(maxlen)
   @c glGetShaderInfoLog(shader, maxlen, &len, message)
   return String(message[1:len])
end

@export function getprogramiv(program::Cuint, name::GLenum)
   param::Cint = 0
   @c glGetProgramiv(program, name, &param)
   return param
end

function getprogramiv(program::Cuint, name::Symbol)
   @assert name in keys(ProgramParamMap) "Unrecognized program parameter: $name"
   param = ProgramParamMap[name]
   return getprogramiv(program, param)
end

@export islinked(program::Cuint) = getprogramiv(program, :link_status)     == GL_TRUE
@export isvalid(program::Cuint)  = getprogramiv(program, :validate_status) == GL_TRUE

@export function programinfo(program::Cuint)
   len::Cint = 0
   maxlen = getprogramiv(program, :info_log_length)
   message = StringArray(maxlen)
   @c glGetProgramInfoLog(program, maxlen, &len, message)
   return String(message[1:len])
end

@export function active_uniform(program::Cuint, index::Integer)
   buff_size = getprogramiv(program, :active_uniform_max_length)
   len::Cint = 0
   size::Cint = 0
   type::Cuint = 0
   _name = StringArray(buff_size)

   @c glGetActiveUniform(program, index, buff_size, &len, &size, &type, _name)

   name = String(_name[1:len])
   location = glGetUniformLocation(program, name)

   return name, location, type, size
end

@export function active_uniforms(program::Cuint)
   count = getprogramiv(program, :active_uniforms)
   return map(0:count-1) do index
      active_uniform(program, index)
   end
end

@export function active_attribute(program::Cuint, index::Integer)
   buff_size = getprogramiv(program, :active_attribute_max_length)
   len::Cint = 0
   size::Cint = 0
   type::Cuint = 0
   name = StringArray(buff_size)

   @c glGetActiveAttrib(program, index, buff_size, &len, &size, &type, name)

   return String(name[1:len]), size, type
end

@export function active_attributes(program::Cuint)
   count = getprogramiv(program, :active_uniforms)
   return map(0:count-1) do index
      active_attribute(program, index)
   end
end

@export function create_shader(type::Symbol, source::AbstractString; strict::Bool = true)
   @assert type in keys(ShaderTypeMap) "Unrecognized shader type '$type'"
   shader = glCreateShader(ShaderTypeMap[type])

   @c glShaderSource(shader, 1, Ptr{Cchar}[pointer(source)], C_NULL)
   glCompileShader(shader)

   if !iscompiled(shader) && strict
      info = shaderinfo(shader)
      glDeleteShader(shader)
      error("Failed to compile shader: $info")
   end

   return shader
end

function create_shader(filename::AbstractString, type::Symbol = :auto; strict::Bool = true)
   ext = extof(filename)
   source = read(filename, String)
   if ext != "" && type == :auto
      stype = Symbol(ext[2:end])

      @assert stype in keys(ShaderTypeMap) "
         Param `type` is set to auto, but could not automatically detect
         the type of the shader from the extension `$ext`.

         You can manually provide the shader type using the `type` parameter.
      "

      return create_shader(stype, source; strict)
   else
      return create_shader(type, source; strict)
   end
end

@export function create_program(shaders::Cuint...; strict::Bool = true)
   program = glCreateProgram()

   if length(shaders) > 0
      foreach(shaders) do shader
         glAttachShader(program, shader)
      end

      glLinkProgram(program)
      glValidateProgram(program)

      if strict && !islinked(program)
         info = programinfo(program)
         deleteshaders!(shaders...)
         glDeleteProgram(program)
         error("Failed to link program: $info")
      end

      if strict && !isvalid(program)
         info = programinfo(program)
         deleteshaders!(shaders...)
         glDeleteProgram(program)
         error("Failed to validate program: $info")
      end

      foreach(shaders) do shader
         glDetachShader(program, shader)
         glDeleteShader(shader)
      end
   end

   return program
end

@export bindprogram!(program::Cuint) = glUseProgram(program)
@export deleteshaders!(shaders::Cuint...) = foreach(glDeleteShader, shaders)
@export deleteprograms!(programs::Cuint...) = foreach(glDeleteProgram, programs)

_suffix(::Type{<:Signed}) = "i"
_suffix(::Type{<:Unsigned}) = "ui"
_suffix(::Type{<:AbstractFloat}) = "f"

@export function uniform(program::Cuint, name::StringLike)
   return glGetUniformLocation(program, string(name))
end

@export function uniforms(program::Cuint, names::StringLike...)
   return map(names) do name
      uniform(program, name)
   end
end

@export macro uniform(expr)
   @assert Meta.isexpr(expr, :(=)) "Wrong usage of @uniform macro"

   left, value = expr.args
   @assert Meta.isexpr(left, :(.)) "Wrong usage of @uniform macro"

   program, name = left.args

   :( uniform!(uniform($program, $name), $value) ) |> esc
end

@export macro uniforms(program, block)
   exprs = map(block.args) do expr
      if Meta.isexpr(expr, :(=))
         name, value = expr.args
         return :(@uniform $(program).$(name) = $(value))
      end
      return expr
   end

   quote
      bindprogram!($(program))
      $(exprs...)
   end |> esc
end

@export function uniform!(location::Cint, vec::SVector{N, T}) where {T <: Real, N}
   @assert location > -1  "Attemping to set value of uniform with negative location: $location"
   @assert 4 >= N >= 1 "Expected uniform vector to have 1-4 dimensions, got $N"
   fname = Symbol("glUniform$N" * _suffix(T))
   getfield(ModernGL, fname)(location, vec...)
end

uniform!(location::Cint, values::Real...) = uniform!(location, SVector(values...))
uniform!(location::Cint, vec::Vector{<:Real}) = uniform!(location, SVector(vec...))

function uniform!(location::Cint, matrix::SMatrix{N, M}; transpose = false) where {N, M}
   @assert location > -1 "Attemping to set value of uniform with negative location: $location"
   @assert 4 >= N >= 1 "Expected uniform matrice to have 1-4 columns, got $N"
   @assert 4 >= M >= 1 "Expected uniform matrice to have 1-4 rows, got $M"

   suffix = N == M ? "$(N)fv" : "$(N)x$(M)fv"
   fname = Symbol("glUniformMatrix$suffix")
   getfield(ModernGL, fname)(location, 1, gbool(transpose), matrix)
end

max_stars_repo_name: jorge-brito/Glutils

max_issues_repo_name: jorge-brito/Glutils

max_forks_repo_name: jorge-brito/Glutils


text:
using AlgRemez
using Test

@testset "AlgRemez.jl" begin
    # Write your tests here.
    y = 1
    z = 2
    n = 5
    lambda_low = 0.0004
    lambda_high = 64
    println("The order is $n")
    coeff_plus,coeff_minus = AlgRemez.calc_coefficients(y,z,n,lambda_low,lambda_high)
    println(coeff_plus.α)
    println(coeff_plus.β)

    funcplus = AlgRemez.fittedfunction(coeff_plus)
    println("1: sqrt(x)")
    println("Approximated value at x=0.1: ",funcplus(0.1))
    println("Exact value at x=0.1: ", sqrt(0.1))

    funcminus = AlgRemez.fittedfunction(coeff_minus)
    println("2: 1/sqrt(x)")
    println("Approximated value at x=0.1: ",funcminus(0.1))
    println("Exact value at x=0.1: ", 1/sqrt(0.1))

    n = 15
    println("The order is $n")
    coeff_plus,coeff_minus = AlgRemez.calc_coefficients(y,z,n,lambda_low,lambda_high)
    funcplus = AlgRemez.fittedfunction(coeff_plus)
    println("1: sqrt(x)")
    println("Approximated value at x=0.1: ",funcplus(0.1))
    println("Exact value at x=0.1: ", sqrt(0.1))

    funcminus = AlgRemez.fittedfunction(coeff_minus)
    println("2: 1/sqrt(x)")
    println("Approximated value at x=0.1: ",funcminus(0.1))
    println("Exact value at x=0.1: ", 1/sqrt(0.1))

    y = 1
    z = 4
    n = 5
    lambda_low = 0.0004
    lambda_high = 64
    println("The order is $n")
    coeff_plus,coeff_minus = AlgRemez.calc_coefficients(y,z,n,lambda_low,lambda_high)
    println(coeff_plus.α)
    println(coeff_plus.β)

    funcplus = AlgRemez.fittedfunction(coeff_plus)
    println("3: x^($y/$z)")
    println("Approximated value at x=0.1: ",funcplus(0.1))
    println("Exact value at x=0.1: ", 0.1^(y/z))

    funcminus = AlgRemez.fittedfunction(coeff_minus)
    println("4: x^(-$y/$z)")
    println("Approximated value at x=0.1: ",funcminus(0.1))
    println("Exact value at x=0.1: ", 0.1^(-y/z))

    n = 15
    println("The order is $n")
    coeff_plus,coeff_minus = AlgRemez.calc_coefficients(y,z,n,lambda_low,lambda_high)
    funcplus = AlgRemez.fittedfunction(coeff_plus)
    println("3: x^($y/$z)")
    println("Approximated value at x=0.1: ",funcplus(0.1))
    println("Exact value at x=0.1: ", 0.1^(y/z))

    funcminus = AlgRemez.fittedfunction(coeff_minus)
    println("4: x^(-$y/$z)")
    println("Approximated value at x=0.1: ",funcminus(0.1))
    println("Exact value at x=0.1: ", 0.1^(-y/z))
end


text:
using Gadfly

function lossDemo()
    qs = [.2, 1, 2]
    loss(x, q) = abs(x).^q

    x = -2:0.01:2

    function lossPlot(q)
        plot(x=x, y=loss(x, q), Geom.line, Theme(default_color=colorant"black", line_width=1px))
    end

    [lossPlot(q) for q in qs]
end

function hingeLossDemo(ε)
    hinge(x) = (abs(x)-ε).*(abs(x).>=ε)

    x = -2:0.01:2
    plot(x=x, y=hinge(x), Geom.line, Theme(default_color=colorant"black", line_width=1px))
end


text:
# Tests for Generalized DF DEA Models
@testset "GeneralizedDFDEAModel" begin

    ## Test Generalized DF DEA Model with Zofío and Prieto (2006) data
    X = [5 3; 2 4; 4 2; 4 8; 7 9]
    Y = [7 4; 10 8; 8 10; 5 4; 3 6]

    # alpha = 0.5 CRS equals Input Oriented CRS
    deaio = dea(X, Y, orient = :Input, rts = :CRS)
    deagdf05crs = deagdf(X, Y, 0.5, rts = :CRS)
    @test efficiency(deaio) ≈ efficiency(deagdf05crs) atol = 1e-7

    @test nobs(deagdf05crs) == 5
    @test ninputs(deagdf05crs) == 2
    @test noutputs(deagdf05crs) == 2

    # alphpa = 0.5 VRS
    deagdf05vrs = deagdf(X, Y, 0.5, rts = :VRS)

    @test nobs(deagdf05vrs) == 5
    @test ninputs(deagdf05vrs) == 2
    @test noutputs(deagdf05vrs) == 2
    @test efficiency(deagdf05vrs) ≈ [0.682;
                               1;
                               1;
                               0.250;
                               0.360] atol = 1e-3

    # Test no slacks
    deagdfnoslack = deagdf(X, Y, 0.5, rts = :VRS, slack = false)
    @test efficiency(deagdfnoslack) == efficiency(deagdf05vrs)
    @test isempty(slacks(deagdfnoslack, :X)) == 1
    @test isempty(slacks(deagdfnoslack, :Y)) == 1

     ## Test if one-by-one DEA using evaluation and reference sets match initial results
     deagdf05crs_ref_eff = zeros(size(X, 1))

     deagdf05vrs_ref_eff = zeros(size(X, 1))

     Xref = X[:,:]
     Yref = Y[:,:]

     for i = 1:size(X, 1)
         Xeval = X[i:i,:]
         Xeval = Xeval[:,:]
         Yeval = Y[i:i,:]
         Yeval = Yeval[:,:]

         deagdf05crs_ref_eff[i] = efficiency(deagdf(Xeval, Yeval, 0.5, rts = :CRS, Xref = Xref, Yref = Yref))[1]

         deagdf05vrs_ref_eff[i] = efficiency(deagdf(Xeval, Yeval, 0.5, rts = :VRS, Xref = Xref, Yref = Yref))[1]
     end

     @test deagdf05crs_ref_eff ≈ efficiency(deagdf05crs)
     @test deagdf05vrs_ref_eff ≈ efficiency(deagdf05vrs)

    # Print
    show(IOBuffer(), deagdf05crs)
    show(IOBuffer(), deagdfnoslack)

    # Test errors
    @test_throws ErrorException deagdf([1; 2 ; 3], [4 ; 5], 0.5) #  Different number of observations
    @test_throws ErrorException deagdf([1; 2], [4 ; 5], 0.5, Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws ErrorException deagdf([1 1; 2 2], [4 4; 5 5], 0.5, Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws ErrorException deagdf([1 1; 2 2], [4 4; 5 5], 0.5, Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ErrorException deagdf([1; 2; 3], [4; 5; 6], 0.5, rts = :Error) # Invalid returns to scale

end


text:
@inline function rk4(f, x, τ)
    τ½ = τ / 2

    k = f(x)
    dx = @. k / 6

    k = f(@. x + τ½ * k)
    dx = @. dx + k / 3

    k = f(@. x + τ½ * k)
    dx = @. dx + k / 3

    k = f(@. x + τ * k)
    dx = @. dx + k / 6

    return @. x + τ * dx
end

function double_gyre(x, y, t; A = 0.25, ϵ = 0.25, ω = 2π)
    f(x, t) = ϵ * sin.(ω * t) .* x .^ 2 .+ (1 .- 2ϵ * (sin.(ω * t))) .* x
    df(x, t) = 2 * ϵ * sin.(ω * t) .* x .+ (1 .- 2ϵ * (sin.(ω * t)))

    return -π * A * sin.(π * f(x, t)) .* cos.(π * y),
    π * A * cos.(π * f(x, t)) .* sin.(π * y) .* df(x, t),
    1
end

function double_gyre_map(
    x;
    steps = 40,
    h = 0.25,
    return_intermediate = false,
    A = 0.25,
    ϵ = 0.25,
    ω = 2π,
    invert_time = false,
)
    factor = invert_time ? -1 : 1
    f = λ -> factor .* double_gyre(λ[1], λ[2], λ[3], A = A, ϵ = ϵ, ω = ω)
    if return_intermediate
        values = zeros(Float64, (steps, 3))
    end
    for i = 1:steps
        x = rk4(f, x, h)
        if return_intermediate
            values[i, :] = x
        end
    end
    if return_intermediate
        return values
    else
        return x
    end
end


text:
using Documenter, CUDAnative

makedocs(
    modules = [CUDAnative],
    format = :html,
    sitename = "CUDAnative.jl",
    pages = [
        "Home"    => "index.md",
        "Manual"  => [
            "man/usage.md",
            "man/troubleshooting.md",
            "man/performance.md",
            "man/hacking.md"
        ],
        "Library" => [
            "lib/compilation.md",
            "lib/reflection.md",
            "lib/profiling.md",
            "Device Code" => [
                "lib/device/intrinsics.md",
                "lib/device/array.md",
                "lib/device/libdevice.md"
            ]
        ]
    ],
    doctest = true
)

deploydocs(
    repo = "github.com/JuliaGPU/CUDAnative.jl.git",
    julia = "",
    osname = "",
    # no need to build anything here, re-use output of `makedocs`
    target = "build",
    deps = nothing,
    make = nothing
)


text:
module PkgTest

using BenchmarkTools, Random, Statistics, Test, OTSM
using COSMO, SCS #, MosekTools
using LinearAlgebra

@testset "cat" begin
A = Matrix{Matrix{Int}}(undef, 2, 2)
A[1, 1] = [1 1 1; 1 1 1]
A[1, 2] = [2 2; 2 2]
A[2, 1] = [3 3 3]
A[2, 2] = [3 3]
@test cat(A) == [1 1 1 2 2; 1 1 1 2 2; 3 3 3 3 3]
end

@testset "otsm_pba (port wine, MAXDIFF)" begin
A, maxdiff_optim, _ = portwine_data()
Smaxdiff = [A[i]'A[j] for i in 1:4, j in 1:4]
OTSM.verify_input_data!(Smaxdiff, 3, true) # set S[i, i] to 0
for r in 1:3, init_fun in [init_eye, init_sb, init_tb]
    @info "rank = $r, init = $init_fun"
    @time Ô_ba, ts_ba, = otsm_pba(Smaxdiff, r;
    verbose = true, O = init_fun(Smaxdiff, r))
    @show ts_ba
    @test ts_ba ≈ maxdiff_optim[r]
    @show test_optimality(Ô_ba, Smaxdiff, method = :wzl)[[1, end]]
    @show test_optimality(Ô_ba, Smaxdiff, method = :lww)[[1, end]]
end
# bm = @benchmark tsm_ba(S, 1)
# display(bm); println()
end

@testset "otsm_pba (port wine, MAXBET)" begin
A, _, maxbet_optim = portwine_data()
Smaxbet = [A[i]'A[j] for i in 1:4, j in 1:4]
for r in 1:3, init_fun in [init_eye, init_sb, init_tb, init_lww1]
    @info "rank = $r, init = $init_fun"
    @time Ô_ba, ts_ba, = otsm_pba(Smaxbet, r;
    verbose = true, O = init_fun(Smaxbet, r))
    @show ts_ba
    if r == 3 && init_fun ∈ [init_eye, init_lww1]
        # r=3, init=init_eye or init_lww1 leads to an inferior local maximum
        @test ts_ba ≈ 818.7063749651282
    else
        @test ts_ba ≈ maxbet_optim[r]
    end
    # check global optimality certificate
    @show test_optimality(Ô_ba, Smaxbet, method = :wzl)[[1, end]]
    @show test_optimality(Ô_ba, Smaxbet, method = :lww)[[1, end]]
end
# bm = @benchmark tsm_ba(S, 1)
# display(bm); println()
end

@testset "SDP relaxation (port wine, MAXDIFF)" begin
A, maxdiff_optim, _ = portwine_data()
Smaxdiff = [A[i]'A[j] for i in 1:4, j in 1:4]
OTSM.verify_input_data!(Smaxdiff, 3, true) # set S[i, i] to 0
for r in 1:3
    for solver in [
        #Mosek.Optimizer(LOG=0),
        COSMO.Optimizer(max_iter=5000, verbose=false), # takes long at r = 1
        SCS.Optimizer(verbose=0)
        ]
        @info "rank = $r, solver = $solver"
        @time Ô_sdp, ts_sdp, = otsm_sdp(Smaxdiff, r, solver = solver)
        @show ts_sdp
        @test ts_sdp ≈ maxdiff_optim[r]
    end
end
end

@testset "Procrustes (random data)" begin
Random.seed!(123)
n, d, m = 10, 3, 50
A, S, Ā, A_manopt = generate_procrustes_data(m, n, d)
@info "Proximal block ascent algorithm:"
for init_fun in [init_eye, init_sb, init_tb]
    @info "init = $init_fun"
    @time Ô_ba, ts_ba, = otsm_pba(S, d; verbose = true, O = init_fun(S, d))
    @show ts_ba
    @show test_optimality(Ô_ba, S, method = :wzl)[[1, end]]
    @show test_optimality(Ô_ba, S, method = :lww)[[1, end]]
end
@info "SDP relaxation (failed to find global solution):"
for solver in [
    #Mosek.Optimizer(LOG=0),
    COSMO.Optimizer(max_iter=5000, verbose=false),
    SCS.Optimizer(verbose=0)
    ]
    @info "solver = $solver"
    @time Ô_sdp, ts_sdp, = otsm_sdp(S, d, solver = solver)
    @show ts_sdp
end
end

@testset "MaxBet (Liu-Wang-Wang Example 5.1)" begin
A = [4.3299 2.3230 -1.3711 -0.0084 -0.7414;
2.3230 3.1181 1.0959 0.1285 0.0727;
-1.3711 1.0959 6.4920 -1.9883 -0.1878;
-0.0084 0.1285 -1.9883 2.4591 1.8463;
-0.7414 0.0727 -0.1878 1.8463 5.8875]
ns    = [2, 3]
m     = length(ns)
r     = 1
start = [1; cumsum(ns)[1:end-1] .+ 1]
stop  = cumsum(ns)
S     = [A[start[i]:stop[i], start[j]:stop[j]] for i in 1:m, j in 1:m]
@info "Proximal block ascent algorithm:"
for init_fun in [init_eye, init_sb, init_tb, init_lww1]
    @info "init = $init_fun"
    @time Ô_ba, ts_ba, = otsm_pba(S, r; verbose = true, O = init_fun(S, r))
    @show test_optimality(Ô_ba, S, method = :wzl)[[1, end]]
    @show test_optimality(Ô_ba, S, method = :lww)[[1, end]]
    if init_fun == init_eye
        # init_eye leads to a inferior local maximum
        @test abs(ts_ba - 14.10123) < 1e-3
    elseif init_fun == init_sb
        # the test matrix by Shapiro-Botha (init_sb) is singular with largest eigenvalue
        # 0 (with multiplicity 3). Different machines give different init_sb. Some of
        # them lead to inferior local maximum.
        @test abs(ts_ba - 14.10123) < 1e-3 || abs(ts_ba - 14.73022) < 1e-3
    else
        @test abs(ts_ba - 14.73022) < 1e-3
    end
end
end

@testset "MaxBet (Liu-Wang-Wang Example 5.2)" begin
A = Float64.([45 -20 5 6 16 3;
-20 77 -20 -25 -8 -21;
5 -20 74 47 18 -32;
6 -25 47 54 7 -11;
16 -8 18 7 21 -7;
3 -21 -32 -11 -7 70])
ns    = [2, 2, 2]
m     = length(ns)
r     = 1
start = [1; cumsum(ns)[1:end-1] .+ 1]
stop  = cumsum(ns)
S     = [A[start[i]:stop[i], start[j]:stop[j]] for i in 1:m, j in 1:m]
@info "Proximal block ascent algorithm:"
for init_fun in [init_eye, init_sb, init_tb, init_lww1]
    @info "init = $init_fun"
    @time Ô_ba, ts_ba, = otsm_pba(S, r; verbose = true, O = init_fun(S, r))
    @test abs(ts_ba - 378.9624) < 1e-3
    @show test_optimality(Ô_ba, S, method = :wzl)[[1, end]]
    @show test_optimality(Ô_ba, S, method = :lww)[[1, end]]
end
end

@testset "Ten Berge 1977 p270" begin
m    = 3
d    = 3
A    = Vector{Matrix{Float64}}(undef, d)
A[1] = Matrix([Matrix(I, d, d); Matrix(I, d, d); zeros(d, d)])
A[2] = Matrix([-Matrix(I, d, d); zeros(d, d); Matrix(I, d, d)])
A[3] = Matrix([zeros(d, d); Matrix(I, d, d); Matrix(I, d, d)])
Stb  = [A[i]'A[j] for i in 1:m, j in 1:m]
OTSM.verify_input_data!(Stb, 1, true) # make it MAXDIFF
for r in 1:3, init in [init_eye, init_tb, init_sb]
    Ô_ba, ts_ba, obj_ba, history = otsm_pba(Stb, r; verbose = false, log = true, O = init(Stb, r))
    println("r = $r, init = $init, obj = $ts_ba, niters = $(length(history[:tracesum]))")
    @show test_optimality(Ô_ba, Stb)[1]
end
end

@testset "Procrustes (restrict to SO(d))" begin
# under this setting, the solution from `otsm_pba(soconstr=false)` has some Oi
# with determinant -1;
# but `otsm_pba(soconstr=true)` will enforce the solution in SO(d)
Random.seed!(123)
n, d, m = 10, 3, 50
A, S, Ā, A_manopt = generate_procrustes_data(m, n, d, 1.0)
@info "Proximal block ascent algorithm:"
for init_fun in [init_eye, init_sb, init_tb]
    @info "init = $init_fun"
    @time Ô_ba, ts_ba, = otsm_pba(S, d; verbose = true, O = init_fun(S, d))
    @test any(det.(Ô_ba) .< 0)
    @time Ô_ba, ts_ba, = otsm_pba(S, d; verbose = true, O = init_fun(S, d), soconstr=true)
    # @show det.(Ô_ba)
    @test all(det.(Ô_ba) .≥ 0)
    @show test_optimality(Ô_ba, S, method = :wzl)[[1, end]]
    @show test_optimality(Ô_ba, S, method = :lww)[[1, end]]
end
end

end # PkgTest module


text:
import Base.(<)
#=
== : get a new Assoc where all of the elements of input Assoc mataches the given Element.
=#
function <(A::Assoc, E::Union{AbstractString,Number})
    tarIndex = searchsortedfirst(A.val,E)
    rowkey, colkey, valkey = findnz(A.A)
    mapping = find( x-> x < tarIndex, valkey)
    rowkey , colkey = unique(rowkey[mapping]), unique(colkey[mapping])
    sort!(rowkey)
    sort!(colkey)
    return A[rowkey,colkey]
end

<(E::Union{AbstractString,Number},A::Assoc) = (A > E)

########################################################
# D4M: Dynamic Distributed Dimensional Data Model
# Architect: Dr. Jeremy Kepner (kepner@ll.mit.edu)
# Software Engineer: Alexander Chen (alexc89@mit.edu)
########################################################


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule libsqm_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("libsqm")
JLLWrappers.@generate_main_file("libsqm", UUID("12e50552-4c0f-5870-af95-b9f374f51e8d"))
end  # module libsqm_jll


text:
const SRange{T} = Union{OneToSRange{T},UnitSRange{T},StepSRange{T},LinSRange{T},StepSRangeLen{T}}
const MRange{T} = Union{OneToMRange{T},UnitMRange{T},StepMRange{T},LinMRange{T},StepMRangeLen{T}}
const UnionRange{T} = Union{SRange{T},MRange{T}}
const IRange{T} = Union{OneTo{T},UnitRange{T},StepRange{T},LinRange{T}, StepRangeLen{T}}

# TODO: this should be in ArrayInterface
ArrayInterface.can_setindex(::Type{X}) where {X<:AbstractRange} = false

# TODO: this should be in ArrayInterface
ArrayInterface.ismutable(::Type{X}) where {X<:AbstractRange} = false
ArrayInterface.ismutable(::Type{X}) where {X<:MRange} = true

"""
    is_static(x) -> Bool

Returns `true` if `x` is static.
"""
is_static(::T) where {T} = is_static(T)
is_static(::Type{T}) where {T} = false
is_static(::Type{T}) where {T<:SRange} = true

"""
    Mutability
"""
abstract type Mutability end

struct MutableTrait <: Mutability end
const Mutable = MutableTrait()

struct ImmutableTrait <: Mutability end
const Immutable = ImmutableTrait()

struct StaticTrait <: Mutability end
const Static = StaticTrait()

Mutability(::T) where {T} = Mutability(T)
Mutability(::Type{T}) where {T<:AbstractArray} = Immutable
Mutability(::Type{T}) where {T<:Array} = Mutable
Mutability(::Type{T}) where {T<:MRange} = Mutable
Mutability(::Type{T}) where {T<:StaticArray} = Static
Mutability(::Type{T}) where {T<:SRange} = Static

"""
    as_mutable(x)

If `x` is mutable then returns `x`, otherwise returns a comparable but mutable
type to `x`.
"""
as_mutable(x::OneToMRange) = x
as_mutable(x::Union{OneTo,OneToSRange}) = OneToMRange(last(x))

as_mutable(x::UnitMRange) = x
as_mutable(x::Union{UnitRange,UnitSRange}) = UnitMRange(first(x), last(x))

as_mutable(x::StepMRange) = x
as_mutable(x::Union{StepRange,StepSRange}) = StepMRange(first(x), step(x), last(x))

as_mutable(x::LinMRange) = x
as_mutable(x::Union{LinRange,LinSRange}) = LinMRange(first(x), last(x), length(x))

as_mutable(x::StepMRangeLen) = x
as_mutable(x::Union{StepRangeLen,StepSRangeLen}) = StepMRangeLen(first(x), step(x), length(x), x.offset)

"""
    as_immutable(x)

If `x` is immutable then returns `x`, otherwise returns a comparable but immutable
type to `x`.
"""
as_immutable(x::OneTo) = x
as_immutable(x::OneToRange) = OneTo(last(x))

as_immutable(x::UnitRange) = x
as_immutable(x::StaticUnitRange) = UnitRange(first(x), last(x))

as_immutable(x::StepRange) = x
as_immutable(x::AbstractStepRange) = StepRange(first(x), step(x), last(x))

as_immutable(x::LinRange) = x
as_immutable(x::AbstractLinRange) = LinRange(first(x), last(x), length(x))

as_immutable(x::StepRangeLen) = x
as_immutable(x::AbstractStepRangeLen) = StepRangeLen(first(x), step(x), length(x), x.offset)

"""
    as_static(x)

If `x` is static then returns `x`, otherwise returns a comparable but static
type to `x`.
"""
as_static(x::OneToSRange) = x
as_static(x::Union{OneTo,OneToMRange}) = OneToSRange(last(x))

as_static(x::UnitSRange) = x
as_static(x::Union{UnitRange,UnitMRange}) = UnitSRange(first(x), last(x))

as_static(x::StepSRange) = x
as_static(x::Union{StepRange,StepMRange}) = StepSRange(first(x), step(x), last(x))

as_static(x::LinSRange) = x
as_static(x::Union{LinRange,LinMRange}) = LinSRange(first(x), last(x), length(x))

as_static(x::StepSRangeLen) = x
as_static(x::Union{StepRangeLen,StepMRangeLen}) = StepSRangeLen(first(x), step(x), length(x), x.offset)

"""
    is_dynamic(x) -> Bool

Returns true if the size of `x` is dynamic/can change.
"""
is_dynamic(::T) where {T} = is_dynamic(T)
is_dynamic(::Type{T}) where {T<:AbstractArray} = false
is_dynamic(::Type{T}) where {T<:MRange} = true
is_dynamic(::Type{T}) where {T<:Array} = true


text:
using Base.Test
import Roots.fzero

# test problems from Table 1 of paper referenced in fzero.jl
# 1
@test_approx_eq fzero(x -> sin(x) - x/2, pi/2, pi) 1.89549426703398094714
@test_approx_eq Roots.a42a(x -> sin(x) - x/2, pi/2, pi) 1.89549426703398094714

# 2
function test2(x)
    -2*sum([(2i - 5)^2/(x - i^2)^3 for i = 1:20])
end
@test_approx_eq fzero(test2, 1 + 1e-9, 4 - 1e-9) 3.0229153472730568
@test_approx_eq_eps fzero(test2, 100 + 1e-9, 121 - 1e-9) 110.02653274766949 1e-9

# 3
function test3(x, a, b)
    a*x*exp(b*x)
end
@test_approx_eq fzero(x -> test3(x, -40, -1), -9, 31) 0.0
@test_approx_eq fzero(x -> test3(x, -100, -2), -9, 31) 0.0
@test_approx_eq fzero(x -> test3(x, -200, -3), -9, 31) 0.0

# 4
function test4(x, n, a)
    x^n - a
end
@test_approx_eq fzero(x -> test4(x, 4, 0.2), 0, 5) 0.668740304976422
@test_approx_eq fzero(x -> test4(x, 12, 1), 0, 5) 1.0
@test_approx_eq fzero(x -> test4(x, 8, 1), -0.95, 4.05) 1.0
@test_approx_eq fzero(x -> test4(x, 12, 1), -0.95, 4.05) 1.0

# 5
@test_approx_eq fzero(x -> sin(x) - 0.5, 0, 1.5) 0.5235987755982989

# 6
function test6(x, n)
    2x*exp(-n) - 2*exp(-n*x) + 1
end
@test_approx_eq fzero(x -> test6(x, 1), 0, 1) 0.42247770964123665
@test_approx_eq fzero(x -> test6(x, 20), 0, 1) 0.03465735902085385
@test_approx_eq fzero(x -> test6(x, 100), 0, 1) 0.006931471805599453

# 7
function test7(x, n)
    (1 + (1 - n)^2)*x - (1 - n*x)^2
end
@test_approx_eq fzero(x -> test7(x, 5), 0, 1) 0.0384025518406219
@test_approx_eq fzero(x -> test7(x, 10), 0, 1) 0.0099000099980005
@test_approx_eq fzero(x -> test7(x, 20), 0, 1) 0.0024937500390620117

# 8
function test8(x, n)
    x^2 - (1 - x)^n
end
@test_approx_eq fzero(x -> test8(x, 2), 0, 1) 0.5
@test_approx_eq fzero(x -> test8(x, 5), 0, 1) 0.345954815848242
@test_approx_eq fzero(x -> test8(x, 10), 0, 1) 0.24512233375330725
@test_approx_eq fzero(x -> test8(x, 20), 0, 1) 0.16492095727644096

# 9
function test9(x, n)
    (1 + (1 - n)^4)*x - (1 - n*x)^4
end
@test_approx_eq fzero(x -> test9(x, 1), 0, 1) 0.2755080409994844
@test_approx_eq fzero(x -> test9(x, 2), 0, 1) 0.1377540204997422
@test_approx_eq fzero(x -> test9(x, 20), 0, 1) 7.668595122185337e-6

# 10
function test10(x, n)
    exp(-n*x)*(x - 1) + x^n
end
@test_approx_eq fzero(x -> test10(x, 1), 0, 1) 0.401058137541547
@test_approx_eq fzero(x -> test10(x, 5), 0, 1) 0.5161535187579336
@test_approx_eq fzero(x -> test10(x, 20), 0, 1) 0.5527046666784878

# 11
function test11(x, n)
    (n*x - 1)/((n - 1)*x)
end
@test_approx_eq fzero(x -> test11(x, 2), 0.01, 1) 0.5
@test_approx_eq fzero(x -> test11(x, 15), 0.01, 1) 0.0666666666666666
@test_approx_eq fzero(x -> test11(x, 20), 0.01, 1) 0.05

# 12
function test12(x, n)
    x^(1/n) - n^(1/n)
end
@test_approx_eq fzero(x -> test12(x, 2), 1, 100) 2
@test_approx_eq fzero(x -> test12(x, 6), 1, 100) 6
@test_approx_eq fzero(x -> test12(x, 33), 1, 100) 33

# 13
function test13(x)
    if x == 0
        return 0
    else
        return x/exp(x^2)
    end
end
@test_approx_eq fzero(test13, -1, 4) 0

# 14
function test14(x, n)
    if x >= 0
        return n/20*(x/1.5 + sin(x) - 1)
    else
        return -n/20
    end
end
@test_approx_eq fzero(x -> test14(x, 1), -1e4, pi/2) 0.6238065189616123
@test_approx_eq fzero(x -> test14(x, 40), -1e4, pi/2) 0.6238065189616123

# 15
function test15(x, n)
    if x > 2e-3/(1+ n)
        return e - 1.859
    elseif x < 0
        return -0.859
    else
        return exp((n + 1)*x/2*1e3) - 1.859
    end
end
@test_approx_eq fzero(x -> test15(x, 20), -1e4, 1e4) 0.00005905130559421971
@test_approx_eq fzero(x -> test15(x, 40), -1e4, 1e4) 0.000030245790670210097
@test_approx_eq fzero(x -> test15(x, 100), -1e4, 1e4) 0.000012277994232461523
@test_approx_eq fzero(x -> test15(x, 1000), -1e4, 1e4) 1.2388385788997142e-6

## various tests
## issue #29, basically f(a) or f(b) so big we get NaN for guess from bracket.
f = x ->  x + exp(x)
@test_approx_eq fzero(f, 0, [-1e6, 1e6]) -0.5671432904097838

## test infinite range
@test_approx_eq fzero(x -> x, [-Inf, Inf]) 0.0

## issue with fzeros and roots near 'b'
@test 0 <  maximum(fzeros(x -> sin(x) - 1/1000*x, 0, pi)) < pi

## test secant_method_no_bracket
@test_throws Roots.Dithering Roots.secant_method_no_bracket(x -> (cos(x) * cosh(x))^2, 4.5, 5.0)
@test_throws Roots.StateConverged Roots.secant_method_no_bracket(x -> (cos(x) * cosh(x))^2, 4.5, 4.8, ftol=cbrt(eps()))


text:
export flattenednames

flattenednames(nt::Type{NTuple{N, T}}; prefix = "") where {N, T} =
    Iterators.flatten([
        flattenednames(T; prefix = "$(prefix)[$i]") for i in 1:N
    ]) |> collect
flattenednames(::Type{NamedTuple{(), Tuple{}}}; prefix = "") = ()
flattenednames(::Type{T}; prefix = "") where {T <: Real} = (prefix,)
flattenednames(::Type{T}; prefix = "") where {T <: SArray} =
    ntuple(i -> "$prefix[$i]", length(T))
function flattenednames(::Type{T}; prefix = "") where {T <: SHermitianCompact}
    N = size(T, 1)
    [["$prefix[$i,$j]" for i in j:N] for j in 1:N] |>
    Iterators.flatten |>
    collect
end
function flattenednames(::Type{T}; prefix = "") where {T <: NamedTuple}
    map(1:fieldcount(T)) do i
        Ti = fieldtype(T, i)
        name = fieldname(T, i)
        flattenednames(
            Ti,
            prefix = prefix == "" ? string(name) : string(prefix, '.', name),
        )
    end |>
    Iterators.flatten |>
    collect
end