text:
"""
    VF2

An empty concrete type used to dispatch to [`vf2`](@ref) isomorphism functions.
"""
struct VF2 <: IsomorphismAlgorithm end

"""
    VF2State{G, T}

Structure that is internally used by vf2
"""
struct VF2State{G, T}
    g1::G
    g2::G
    core_1::Vector{T}
    core_2::Vector{T}
    in_1::Vector{T}
    in_2::Vector{T}
    out_1::Vector{T}
    out_2::Vector{T}

    function VF2State(g1::G, g2::G) where {G <: AbstractSimpleGraph{T}} where {T <: Integer}
        n1 = nv(g1)
        n2 = nv(g2)
        core_1 = zeros(T, n1)
        core_2 = zeros(T, n2)
        in_1 = zeros(T, n1)
        in_2 = zeros(T, n2)
        out_1 = zeros(T, n1)
        out_2 = zeros(T, n2)

        return new{G, T}(g1, g2, core_1, core_2, in_1, in_2, out_1, out_2)
    end
end

"""
    vf2(callback, g1, g2, problemtype; vertex_relation=nothing, edge_relation=nothing)

Iterate over all isomorphism between the graphs `g1` (or subgraphs thereof) and `g2`.
The problem that is solved depends on the value of `problemtype`:
- IsomorphismProblem(): Only isomorphisms between the whole graph `g1` and `g2` are considered.
- SubGraphIsomorphismProblem(): All isomorphism between subgraphs of `g1` and `g2` are considered.
- InducedSubGraphIsomorphismProblem(): All isomorphism between vertex induced subgraphs of `g1` and `g2` are considered.

Upon finding an isomorphism, the function `callback` is called with a vector `vmap` as an argument.
`vmap` is a vector where `vmap[v] == u` means that vertex `v` in `g2` is mapped to vertex `u` in `g1`.
If the algorithm should look for another isomorphism, then this function should return `true`.

### Optional Arguments
- `vertex_relation`: A binary function that takes a vertex from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched vertices.
- `edge_relation`: A binary function that takes an edge from `g1` and one from `g2`. An
    isomorphism only exists if this function returns `true` for all matched edges.

### References
Luigi P. Cordella, Pasquale Foggia, Carlo Sansone, Mario Vento
“A (Sub)Graph Isomorphism Algorithm for Matching Large Graphs”
"""
function vf2(callback::Function, g1::G, g2::G, problemtype::GraphMorphismProblem;
             vertex_relation::Union{Nothing, Function}=nothing,
             edge_relation::Union{Nothing, Function}=nothing) where {G <: AbstractSimpleGraph}
    if nv(g1) < nv(g2) || (problemtype == IsomorphismProblem() && nv(g1) != nv(g2))
        return
    end

    start_state = VF2State(g1, g2)
    start_depth = 1
    vf2match!(start_state, start_depth, callback, problemtype, vertex_relation, edge_relation)
    return
end

"""
    vf2check_feasibility(u, v, state, problemtype, vertex_relation, edge_relation)

Check whether two vertices of G₁ and G₂ can be matched. Used by [`vf2match!`](@ref).
"""
function vf2check_feasibility(u, v, state::VF2State, problemtype,
                              vertex_relation::Union{Nothing, Function},
                              edge_relation::Union{Nothing, Function})
    @inline function vf2rule_pred(u, v, state::VF2State, problemtype)
        if problemtype != SubGraphIsomorphismProblem()
            @inbounds for u2 in inneighbors(state.g1, u)
                if state.core_1[u2] != 0
                    found = false
                    # TODO can probably be replaced with has_edge for better performance
                    for v2 in inneighbors(state.g2, v)
                        if state.core_1[u2] == v2
                            found = true
                            break
                        end
                    end
                    found || return false
                end
            end
        end
        @inbounds for v2 in inneighbors(state.g2, v)
            if state.core_2[v2] != 0
                found = false
                for u2 in inneighbors(state.g1, u)
                    if state.core_2[v2] == u2
                        found = true
                        break
                    end
                end
                found || return false
            end
        end
        return true
    end

    @inline function vf2rule_succ(u, v, state::VF2State, problemtype)
        if problemtype != SubGraphIsomorphismProblem()
            @inbounds for u2 in outneighbors(state.g1, u)
                if state.core_1[u2] != 0
                    found = false
                    for v2 in outneighbors(state.g2, v)
                        if state.core_1[u2] == v2
                            found = true
                            break
                        end
                    end
                    found || return false
                end
            end
        end
        found = false
        @inbounds for v2 in outneighbors(state.g2, v)
            if state.core_2[v2] != 0
                found = false
                for u2 in outneighbors(state.g1, u)
                    if state.core_2[v2] == u2
                        found = true
                        break
                    end
                end
                found || return false
            end
        end
        return true
    end

    @inline function vf2rule_in(u, v, state::VF2State, problemtype)
        count1 = 0
        count2 = 0
        @inbounds for u2 in outneighbors(state.g1, u)
            if state.in_1[u2] != 0 && state.core_1[u2] == 0
                count1 += 1
            end
        end
        @inbounds for v2 in outneighbors(state.g2, v)
            if state.in_2[v2] != 0 && state.core_2[v2] == 0
                count2 += 1
            end
        end
        if problemtype == IsomorphismProblem()
            count1 == count2 || return false
        else
            count1 >= count2 || return false
        end
        count1 = 0
        count2 = 0
        @inbounds for u2 in inneighbors(state.g1, u)
            if state.in_1[u2] != 0 && state.core_1[u2] == 0
                count1 += 1
            end
        end
        @inbounds for v2 in inneighbors(state.g2, v)
            if state.in_2[v2] != 0 && state.core_2[v2] == 0
                count2 += 1
            end
        end
        problemtype == IsomorphismProblem() && return count1 == count2

        return count1 >= count2
    end

    @inline function vf2rule_out(u, v, state::VF2State, problemtype)
        count1 = 0
        count2 = 0
        @inbounds for u2 in outneighbors(state.g1, u)
            if state.out_1[u2] != 0 && state.core_1[u2] == 0
                count1 += 1
            end
        end
        @inbounds for v2 in outneighbors(state.g2, v)
            if state.out_2[v2] != 0 && state.core_2[v2] == 0
                count2 += 1
            end
        end
        if problemtype == IsomorphismProblem()
            count1 == count2 || return false
        else
            count1 >= count2 || return false
        end

        count1 = 0
        count2 = 0
        @inbounds for u2 in inneighbors(state.g1, u)
            if state.out_1[u2] != 0 && state.core_1[u2] == 0
                count1 += 1
            end
        end
        @inbounds for v2 in inneighbors(state.g2, v)
            if state.out_2[v2] != 0 && state.core_2[v2] == 0
                count2 += 1
            end
        end
        problemtype == IsomorphismProblem() && return count1 == count2

        return count1 >= count2
    end

    @inline function vf2rule_new(u, v, state::VF2State, problemtype)
        problemtype == SubGraphIsomorphismProblem() && return true
        count1 = 0
        count2 = 0
        @inbounds for u2 in inneighbors(state.g1, u)
            if state.in_1[u2] == 0 && state.out_1[u2] == 0
                count1 += 1
            end
        end
        @inbounds for v2 in inneighbors(state.g2, v)
            if state.in_2[v2] == 0 && state.out_2[v2] == 0
                count2 += 1
            end
        end
        if problemtype == IsomorphismProblem()
            count1 == count2 || return false
        else
            count1 >= count2 || return false
        end
        count1 = 0
        count2 = 0
        @inbounds for u2 in outneighbors(state.g1, u)
            if state.in_1[u2] == 0 && state.out_1[u2] == 0
                count1 += 1
            end
        end
        @inbounds for v2 in outneighbors(state.g2, v)
            if state.in_2[v2] == 0 && state.out_2[v2] == 0
                count2 += 1
            end
        end
        problemtype == IsomorphismProblem() && return count1 == count2

        return count1 >= count2
    end

    @inline function vf2rule_self_loops(u, v, state, problemtype)
        u_selflooped = has_edge(state.g1, u, u)
        v_selflooped = has_edge(state.g2, v, v)

        if problemtype == SubGraphIsomorphismProblem()
            return u_selflooped || !v_selflooped
        end
        return u_selflooped == v_selflooped
    end

    syntactic_feasability = vf2rule_pred(u, v, state, problemtype) &&
                            vf2rule_succ(u, v, state, problemtype) &&
                            vf2rule_in(u, v, state, problemtype)   &&
                            vf2rule_out(u, v, state, problemtype)  &&
                            vf2rule_new(u, v, state, problemtype)  &&
                            vf2rule_self_loops(u, v, state, problemtype)
    syntactic_feasability || return false

    if vertex_relation != nothing
        vertex_relation(u, v) || return false
    end
    if edge_relation != nothing
        E1 = edgetype(state.g1)
        E2 = edgetype(state.g2)
        for u2 in outneighbors(state.g1, u)
            state.core_1[u2] == 0 && continue
            v2 = state.core_1[u2]
            edge_relation(E1(u, u2), E2(v, v2)) || return false
        end
        for u2 in inneighbors(state.g1, u)
            state.core_1[u2] == 0 && continue
            v2 = state.core_1[u2]
            edge_relation(E1(u2, u), E2(v2, v)) || return false
        end
    end
    return true
end

"""
    vf2update_state!(state, u, v, depth)

Update state before recursing. Helper function for [`vf2match!`](@ref).
"""
function vf2update_state!(state::VF2State, u, v, depth)
@inbounds begin
     state.core_1[u] = v
     state.core_2[v] = u
     for w in outneighbors(state.g1, u)
         if state.out_1[w] == 0
             state.out_1[w] = depth
         end
     end
     for w in inneighbors(state.g1, u)
         if state.in_1[w] == 0
             state.in_1[w] = depth
         end
     end
     for w in outneighbors(state.g2, v)
         if state.out_2[w] == 0
             state.out_2[w] = depth
         end
     end
     for w in inneighbors(state.g2, v)
         if state.in_2[w] == 0
             state.in_2[w] = depth
         end
     end
end
end

"""
    vf2reset_state!(state, u, v, depth)

Reset state after returning from recursion. Helper function for [`vf2match!`](@ref).
"""
function vf2reset_state!(state::VF2State, u, v, depth)
@inbounds begin
    state.core_1[u] = 0
    state.core_2[v] = 0
    for w in outneighbors(state.g1, u)
        if state.out_1[w] == depth
            state.out_1[w] = 0
        end
    end
    for w in inneighbors(state.g1, u)
        if state.in_1[w] == depth
            state.in_1[w] = 0
        end
    end
    for w in outneighbors(state.g2, v)
        if state.out_2[w] == depth
            state.out_2[w] = 0
        end
    end
    for w in inneighbors(state.g2, v)
        if state.in_2[w] == depth
            state.in_2[w] = 0
        end
    end
end
end

"""
    vf2match!(state, depth, callback, problemtype, vertex_relation, edge_relation)

Perform isomorphic subgraph matching. Called by [`vf2`](@ref).
"""
function vf2match!(state, depth, callback::Function, problemtype::GraphMorphismProblem,
                   vertex_relation, edge_relation)
    n1 = Int(nv(state.g1))
    n2 = Int(nv(state.g2))
    # if all vertices of G₂ are matched we call the callback function. If the
    # algorithm should look for another isomorphism then callback has to return true
    if depth > n2
        keepgoing = callback(state.core_2)
        return keepgoing
    end
    # First we try if there is a pair of unmatched vertices u∈G₁ v∈G₂ that are connected
    # by an edge going out of the set M(s) of already matched vertices
    found_pair = false
    v = 0
    @inbounds for j = 1:n2
        if state.out_2[j] != 0 && state.core_2[j] == 0
            v = j
            break
        end
    end
    if v != 0
        @inbounds for u = 1:n1
            if state.out_1[u] != 0 && state.core_1[u] == 0
                found_pair = true
                if vf2check_feasibility(u, v, state, problemtype, vertex_relation, edge_relation)
                    vf2update_state!(state, u, v, depth)
                    keepgoing = vf2match!(state, depth + 1, callback, problemtype,
                                          vertex_relation, edge_relation)
                    keepgoing || return false
                    vf2reset_state!(state, u, v, depth)
                end
            end
        end
    end
    found_pair && return true
    # If that is not the case we try if there is a pair of unmatched vertices u∈G₁ v∈G₂ that
    # are connected  by an edge coming in from the set M(s) of already matched vertices
    v = 0
    @inbounds for j = 1:n2
        if state.in_2[j] != 0 && state.core_2[j] == 0
            v = j
            break
        end
    end
    if v != 0
        @inbounds for u = 1:n1
            if state.in_1[u] != 0 && state.core_1[u] == 0
                found_pair = true
                if vf2check_feasibility(u, v, state, problemtype, vertex_relation, edge_relation)
                    vf2update_state!(state, u, v, depth)
                    keepgoing = vf2match!(state, depth + 1, callback, problemtype,
                                          vertex_relation, edge_relation)
                    keepgoing || return false
                    vf2reset_state!(state, u, v, depth)
                end
            end
        end
    end
    found_pair && return true
    # If this is also not the case, we try all pairs of vertices u∈G₁ v∈G₂ that are not
    # yet matched
    v = 0
    @inbounds for j = 1:n2
        if state.core_2[j] == 0
            v = j
            break
        end
    end
    if v != 0
        @inbounds for u = 1:n1
            if state.core_1[u] == 0
                if vf2check_feasibility(u, v, state, problemtype, vertex_relation, edge_relation)
                    vf2update_state!(state, u, v, depth)
                    keepgoing = vf2match!(state, depth + 1, callback, problemtype,
                                          vertex_relation, edge_relation)
                    keepgoing || return false
                    vf2reset_state!(state, u, v, depth)
                end
            end
        end
    end
    return true
end

function has_induced_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                                 vertex_relation::Union{Nothing, Function}=nothing,
                                 edge_relation::Union{Nothing, Function}=nothing)::Bool
        result = false
        callback(vmap) = (result = true; return false)
        vf2(callback, g1, g2, InducedSubGraphIsomorphismProblem();
                       vertex_relation=vertex_relation, edge_relation=edge_relation)
        return result
end

function has_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                                 vertex_relation::Union{Nothing, Function}=nothing,
                                 edge_relation::Union{Nothing, Function}=nothing,)::Bool
    result = false
    callback(vmap) = (result = true; return false)
    vf2(callback, g1, g2, SubGraphIsomorphismProblem();
        vertex_relation=vertex_relation, edge_relation=edge_relation)
    return result
end

function has_isomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                         vertex_relation::Union{Nothing, Function}=nothing,
                         edge_relation::Union{Nothing, Function}=nothing)::Bool

    !could_have_isomorph(g1, g2) && return false

    result = false
    callback(vmap) = (result = true; return false)
    vf2(callback, g1, g2, IsomorphismProblem(),
        vertex_relation=vertex_relation,
        edge_relation=edge_relation)
    return result
end

function count_induced_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                                   vertex_relation::Union{Nothing, Function}=nothing,
                                   edge_relation::Union{Nothing, Function}=nothing)::Int
    result = 0
    callback(vmap) = (result += 1; return true)
    vf2(callback, g1, g2, InducedSubGraphIsomorphismProblem(),
        vertex_relation=vertex_relation,
        edge_relation=edge_relation)
    return result
end

function count_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                                vertex_relation::Union{Nothing, Function}=nothing,
                                edge_relation::Union{Nothing, Function}=nothing)::Int
    result = 0
    callback(vmap) = (result += 1; return true)
    vf2(callback, g1, g2, SubGraphIsomorphismProblem(), vertex_relation=vertex_relation, edge_relation=edge_relation)
    return result
end

function count_isomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                        vertex_relation::Union{Nothing, Function}=nothing,
                        edge_relation::Union{Nothing, Function}=nothing)::Int
    !could_have_isomorph(g1, g2) && return 0
    result = 0
    callback(vmap) = (result += 1; return true)
    vf2(callback, g1, g2, IsomorphismProblem(), vertex_relation=vertex_relation, edge_relation=edge_relation)
    return result
end

function all_induced_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                                 vertex_relation::Union{Nothing, Function}=nothing,
                                 edge_relation::Union{Nothing, Function}=nothing)::Channel{Vector{Tuple{eltype(g1),eltype(g2)}}}
        make_callback(c) = vmap -> (put!(c, collect(zip(vmap,1:length(vmap)))), return true)
    T = Vector{Tuple{eltype(g1), eltype(g2)}}
    ch::Channel{T} = Channel(ctype=T) do c
        vf2(make_callback(c), g1, g2, InducedSubGraphIsomorphismProblem(),
                       vertex_relation=vertex_relation,
                       edge_relation=edge_relation)
    end
end

function all_subgraphisomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                              vertex_relation::Union{Nothing, Function}=nothing,
                              edge_relation::Union{Nothing, Function}=nothing)::Channel{Vector{Tuple{eltype(g1), eltype(g2)}}}

    make_callback(c) = vmap -> (put!(c, collect(zip(vmap,1:length(vmap)))), return true)
    T = Vector{Tuple{eltype(g1), eltype(g2)}}
    ch::Channel{T} = Channel(ctype=T) do c
        vf2(make_callback(c), g1, g2, SubGraphIsomorphismProblem(),
            vertex_relation=vertex_relation,
            edge_relation=edge_relation)
    end
    return ch
end

function all_isomorph(g1::AbstractGraph, g2::AbstractGraph, alg::VF2;
                 vertex_relation::Union{Nothing, Function}=nothing,
                 edge_relation::Union{Nothing, Function}=nothing)::Channel{Vector{Tuple{eltype(g1),eltype(g2)}}}
    T = Vector{Tuple{eltype(g1), eltype(g2)}}
    !could_have_isomorph(g1, g2) && return Channel(_ -> return, ctype=T)
    make_callback(c) = vmap -> (put!(c, collect(zip(vmap,1:length(vmap)))), return true)
    ch::Channel{T} = Channel(ctype=T) do c
        vf2(make_callback(c), g1, g2, IsomorphismProblem(),
            vertex_relation=vertex_relation,
            edge_relation=edge_relation)
    end
    return ch
end


text:
using ForwardDiff

"""
    black_ivol(F, K, β)

The inverse black function: returns the value `s` such that

    black(K,s) == β

where `F > 0`, K > 0` and `max(F-K, 0) ≤ β ≤ F`.
"""
function black_ivol(F, K, β)
    F > zero(F) && K > zero(K) || throw(DomainError())
    if !(β < F)
        β == F ? (return Inf) : throw(DomainError())
    end

    if K ≥ F
        if !(zero(β) < β)
            β == zero(β) ? (return 0.0) : throw(DomainError())
        end
        _black_ivol(K/F, β/F)
    else
        L = F-K
        if !(L < β)
            β == L ? (return 0.0) : throw(DomainError())
        end
        _black_ivol(F/K, (β-L)/K)
    end
end

"""
    _black_ivol(K, β)

The kernel of the inverse black function: returns the value `s` such that

    black(K,s) == β

where `K ≥ 1` and `0 ≤ β ≤ 1`.
"""
function _black_ivol(K, β)
    maxiters = 4

    logK = log(K)
    x  = -logK
    s₂ = sqrt(-2*x) # point of inflection
    β₂ = 0.5 - K*Φ(-s₂) # cancelation of some terms
    ν₂ = oftype(β₂,invsqrt2π)
    if β < β₂
        s₁ = s₂ - β₂/ν₂
        h₁ = x/s₁;  t₁ = s₁/2
        β₁ = _black(h₁,t₁,K)
        if β < β₁
            # 1. Rational interpolation of transformation
            # (note scaling vs Jaeckel (2013))
            #   f(β)  = -2π*x/(3√3) * Φ(z)^3
            #            where z = x/(√3*s(β))
            #         -> 0   as β -> 0
            #   f′(β) = 2π*z² * ϕ(z)² / √K * exp(z² + s²/8)    # (eq. 4.32)/√K
            #         -> 1/√K as β -> 0
            #   f′′(β) =                                       # (eq. 4.33)/K
            #
            β₀ = 0.0; f₀ = 0.0;  d₀ = 1.0/sqrt(K)

            z₁ = x/(sqrt(3)*s₁)
            Φz₁ = Φ(z₁)
            f₁ = (-2π/(sqrt(27))) * x * Φz₁^3
            w₁ = exp(z₁^2 + s₁^2/8)
            d₁ = 2π * z₁^2 * Φz₁^2 / √K * w₁
            dd₁ = (π/6)*(z₁^2/s₁^3) * Φz₁ * w₁^2 / K *
                ((-8*sqrt(3))*s₁*x + (3*s₁^2*(s₁^2-8) - 8*x^2)*Φz₁/ϕ(z₁))

            Δ = (f₁ - f₀)/(β₁ - β₀)
            r = ((β₁ - β₀)/2*dd₁ + (d₁-d₀))/(d₁ - Δ)
            f = rational_cubic_interpolation(β, β₀,β₁, f₀,f₁, d₀,d₁, r)
            z = invΦ(cbrt(sqrt(27)*f/(-2π*x)))
            s = x/(sqrt(3)*z)

            # Use 3rd order Householder's method with objective function
            #     g(s)   = 1/log(b(s)) - 1/log(β)
            #     g′(s)  = -b′ / (b  * log(b)^2)
            #     g′′(s) = g′ * ( b′′/b′ - b′/b * (1+2/log(b)))
            #     g′′′(s) =  g′ *( b′′′/b′ +   2(b'/b)²·(1+3/ln(b)·(1+1/ln(b)))  -  3(b''/b)·(1+2/ln(b)))

            logβ = log(β)
            for iter = 1:maxiters
                h = x/s; t = s/2
                b = _black(h,t,K)
                logb = log(b)
                ν = _black_vega(h,t)
                bH2 = (h+t)*(h-t)/s
                bH3 = bH2^2 - 3*(h/s)^2 - 1/4

                g′ = -ν / (b*logb^2)
                N  = b/ν*(logβ-logb)*logb/logβ
                H2 = bH2 - ((ν/b)*(1+2/logb))
                H3 = bH3 + (ν/b)*(2*(ν/b)*(1+3/logb*(1+1/log(b))) - 3*bH2*(1+2/logb))

                ds = max(-0.5*s , householder_update(N, H2, H3))
                s += ds
                if abs(ds) < eps(typeof(s))*s
                    break
                end
            end
            return s
        else
            # 2. Rational interpolation between s₁ and s₂
            ν₁ = _black_vega(h₁,t₁)
            # choose r to obtain concave interpolant
            # (this is dominant over 2nd deriv zero at centre)
            Δ = (s₂-s₁)/(β₂-β₁)
            m,M = minmax(1/ν₁ - Δ, Δ - 1/ν₂)
            r = 1 + (1/ν₁-Δ)/(Δ-1/ν₂)
            s = rational_cubic_interpolation(β, β₁,β₂, s₁,s₂, 1/ν₁,1/ν₂, r)
            s₋ = s₁; s₊ = s₂
            @goto middle
        end
    else
        s₃ = s₂ + (1-β₂)/ν₂
        h₃ = x/s₃;  t₃ = s₃/2
        β₃ = _black(h₃,t₃,K)
        if β < β₃
            # 3. Direct rational interpolation between s₂ and s₃
            ν₃ = _black_vega(h₃,t₃)
            # choose r to obtain convex interpolant
            # (this is dominant over 2nd deriv zero at centre)
            Δ = (s₃-s₂)/(β₃-β₂)
            m,M = minmax(1/ν₃ - Δ, Δ - 1/ν₂)
            r = 1 + M/m
            s = rational_cubic_interpolation(β, β₂,β₃, s₂,s₃, 1/ν₂,1/ν₃, r)
            s₋ = s₂; s₊ = s₃

            @label middle # regions 2 & 3

            # Objective function
            #     g(s) = b(s) - β

            for iter = 1:maxiters
                h = x/s; t= s/2
                b = _black(h, t, K)
                ν = _black_vega(h, t)
                N = (β - b)/ν                # (β - b) / vega
                H2 = (h+t)*(h-t)/s           # vomma / vega
                H3 = H2^2 - 3*(h/s)^2 - 1/4  # ultima / vega

                ds = max(-0.5*s , householder_update(N, H2, H3))
                s += ds
                if abs(ds) < eps(typeof(s))
                    break
                end
            end
            return s

        else
            # 4. Rational interpolation of transformation:
            #  fᵤ(β)   = Φ(-s(β)/2)
            #           -> 0      as β -> 1
            #  fᵤ'(β)  = -1/2 * exp(1/2 * (x²/s² + x))
            #          = -1/(2*√K) * exp(h²/2)
            #           -> -1/(2*√K)   as β -> 1
            #  fᵤ''(β) = sqrt(π/2) * x² / s³ * exp(x² / s² + x + s²/8)
            #          = sqrt(π/2) * h²/(s*K) * exp(h² + s²/8)
            #
            # (note scaling 1/√K, 1/K vs Jaeckel (2013))

            β₄ = 1.0
            f₄ = 0.0
            d₄ = -1/(2*sqrt(K))

            f₃ = Φ(-t₃)
            d₃ = d₄ * exp(h₃^2/2)
            dd₃ = sqrt(π/2) * h₃^2 / (s₃*K) * exp(h₃^2 + t₃^2/2)

            Δ = (f₄ - f₃)/(β₄ - β₃)
            r = ((β₄-β₃)/2*dd₃ + (d₄-d₃))/(Δ - d₃)
            f = rational_cubic_interpolation(β, β₃,β₄, f₃,f₄, d₃,d₄, r)
            s = -2*invΦ(f)

            # b′ = _black_vega(h,t)           # vega
            # b′′ = b′ * (h+t)*(h-t)/s  # vomma
            # b′′′ = b′ * (4h^2 - (h+t)(h-t)*(1-(h+t)(h-t)))   # ultima

            # Use 3rd order Householder's method with objective function
            #     g(s)   = log(1-β) - log(1-b(s)) = -log(cb/cβ)
            #     g′(s)  = b′ / (1-b)
            #     g′′(s) = g′ * ( b′′/b′ + g′)
            #     g′′′(s) =  g′ *( b′′′/b′ +  g′ * (2g′ + 3b′′/b′))

            cβ = 1-β

            for iter = 1:maxiters
                h = x/s; t = s/2
                cb = _cblack(h,t,K)
                ν = _black_vega(h,t)
                bH2 = (h+t)*(h-t)/s
                bH3 = bH2^2 - 3*(h/s)^2 - 1/4

                g′ = ν / cb
                N  = log(cb/cβ) / g′
                H2 = bH2 + g′
                H3 = bH3 + g′ * (2*g′ + 3*bH2)

                ds = max(-0.5*s , householder_update(N, H2, H3))
                s += ds

                if abs(ds) < eps(typeof(s))
                    break
                end
            end
            return s

        end
    end
end

"""
    householder_update(N, [H2, [H3]])

The update incrememnt for `n`th order Householder's method (where `n` is the number of arguments) for finding the root of

    f(x) == 0

* `n == 1` is Newton's method
* `n == 2` is Halley's method

Arguments:

* `N`:  the Newton update `-f(x) / f′(x)`
* `H2`: the scaled 2nd order derivative `f′′(x) / f′(x)`
* `H3`: the scaled 3rd order derivative `f′′′(x) / f′(x)`
"""
householder_update(N) = N
householder_update(N, H2) = N/(1+H2*N/2)
householder_update(N, H2, H3) = N*(1+H2*N/2)/(1+N*(H2 + H3*N/6))

"""
    rational_cubic_interpolation(x, x₋, x₊, y₋, y₊, d₋, d₊, r)

Compute a shape preserving rational interpolation.

Arguments

* `x` the value at which to interpolate
* `x₋`,`x₊`: the lower and upper bounds of the interval
* `y₋`,`y₊`: the value of the function at the lower and upper bounds
* `d₋`,`d₊`: the value of the derivative at the lower and upper bounds
* `r`: the control parameter, which must be greater than -1.
  - If `r == 3`, then it is the usual Hermite cubic polynomial
  - As `r → ∞`, it approaches linear interpolation.

If `Δ = (y₊ - y₋)/(x₋ - x₊) > 0`, and `d₋ >= 0`, `d₊ >= 0`, then the interpolant is guaranteed monotonic increasing if

    r >= (d₋ + d₊)/Δ

[DG85, eq. 3.8].

If  `d₋ < Δ < d₊`, then the interpolant is convex if and only if

    r >= 1 + max(d₊ - Δ, Δ - d₋) / min(d₊ - Δ, Δ - d₋)

[DG85, eq. 3.18].

Note: for the purposes of computation, the bounds are actually interchangable (i.e. `x₋` can be the upper bound).

References

* R. Delbourgo and J. A. Gregory (1985) "Shape Preserving Piecewise Rational Interpolation", SIAM J. Sci. and Stat. Comput., 6(4), 967–976. doi:[10.1137/0906065](http://dx.doi.org/10.1137/0906065)

"""
function rational_cubic_interpolation(x, x₋, x₊, y₋, y₊, d₋, d₊, r)
    h = x₊ - x₋
    if h == 0
        return (y₋ + y₊)/2
    end

    # r should be greater than -1. We do not use  assert(r > -1)  here in order to allow values such as NaN to be propagated as they should.
    s = (x - x₋) / h
    t = (x₊ - x) / h

    if r < 2/eps(typeof(r))^2
        ((y₊*s + (r*y₊ - h*d₊)*t)*s^2 + ((r*y₋ + h*d₋)*s + y₋*t)*t^2) /
          (1 + (r - 3) * s * t)
    else
        # Linear interpolation without over-or underflow.
        y₊*s + y₋*t
    end
end

ext: jl


text:
test_gage(GeometricCenter,
    [1. 0. 1.;
     0. 1. 0.;],
    [0.5 0.5 0.5],
    "test 2")


text:
export fit,
       #fitA,
       #fitB,
       predict,
       FreeLSStateTomo,
       LSStateTomo,
       MLStateTomo

function build_state_predictor(obs::Vector{Matrix{T}}) where T
    return reduce(vcat,[vec(o)' for o in obs])
end

"""
`FreeLSStateTomo`

Free (unconstrained) least-squares state tomography algorithm. It is
constructed from a vector of observables corresponding to
measurements that are performed on the state being reconstructed.
"""
struct FreeLSStateTomo
    inputdim::Int
    outputdim::Int
    pred::Matrix{ComplexF64}
    function FreeLSStateTomo(obs::Vector)
        @assert all([LinearAlgebra.ishermitian(o) for o in obs]) "Observables must be Hermitian"
        pred = build_state_predictor(obs)
        outputdim = size(pred,1)
        inputdim = size(pred,2)
        new(inputdim,outputdim,pred)
    end
end

"""
`predict`

Predict outcomes of a tomography experiment for a given state (density matrix).
"""
function predict(method::FreeLSStateTomo, state)
    method.pred*vec(state)
end

"""
`fit(method, ...)`

Reconstruct a state from observations (i.e., perform state tomography). Exactly which fitting
routine is used depends on the type of the `method`. The possible types are

   + FreeLSStateTomo : unconstrained least-squares state tomography
   + LSStateTomo : least-squares state tomography constrained to physical states
   + MLStateTomo : maximum-likelihood, maximum-entropy state tomography constrained to physical states

"""
function fit(method::FreeLSStateTomo,
             means::Vector{Float64})
    if length(means) != method.outputdim
        error("The number of expected means does not match the required number of experiments")
    end
    d = round(Int, method.inputdim |> sqrt)
    reg = method.pred\means
    return reshape(reg,d,d), LinearAlgebra.norm(method.pred*reg-means,2)/length(means), :Optimal
end

function fit(method::FreeLSStateTomo,
             means::Vector{Float64},
             vars::Vector{Float64})
    if length(means) != method.outputdim
        error("The number of expected means does not match the required number of experiments")
    end
    d = round(Int, method.inputdim |> sqrt)
    if any(vars .< 0)
        error("Variances must be positive for generalized least squares.")
    end

    reg = (LinearAlgebra.Diagonal(1 ./ sqrt(vars)) * method.pred) \ (LinearAlgebra.Diagonal(1 ./ sqrt(vars)) * means)
    return reshape(reg,d,d),
           sqrt(dot(method.pred*reg-means,LinearAlgebra.Diagonal(vars)\(method.pred*reg-means)))/length(means),
           :Optimal
end

"""
`LSStateTomo`

Type corresponding to constrained least-squares state tomography
algorithm. It is constructed from a vector of observables
corresponding to measurements that are performed on the state being
reconstructed.

The outcome of said measurements for a given density matrix ρ can be
computed using the `predict` function , while the density matrix can
be estimated from observations (including variances) by using the
`fit` function.

"""
struct LSStateTomo
    inputdim::Int
    outputdim::Int
    realpred::Matrix{Float64}
    function LSStateTomo(obs::Vector)
        @assert all([LinearAlgebra.ishermitian(o) for o in obs]) "Observables must be Hermitian"
        pred = build_state_predictor(obs)
        outputdim = size(pred,1)
        inputdim = size(pred,2)
        realpred = [real(pred) imag(pred)];
        new(inputdim,outputdim,realpred)
    end
end

function predict(method::LSStateTomo, state)
    (method.realpred[:,1:round(Int,end/2)]+1im*method.realpred[:,round(Int,end/2)+1:end])*vec(state)
end

function fit(method::LSStateTomo,
             means::Vector{Float64},
             vars::Vector{Float64};
             #solver = MosekSolver(LOG=0))
             solver = SCS.SCSSolver(verbose=0, max_iters=10_000, eps = 1e-8))

    if length(means) != length(vars) || method.outputdim != length(means)
        error("Size of observations and/or predictons do not match.")
    end
    dsq = method.inputdim
    d = round(Int,sqrt(dsq))

    # We assume that the predictions are always real-valued
    # and we need to do the complex->real translation manually since
    # Convex.jl does not support complex numbers yet
    ivars = 1 ./ sqrt.(vars)

    ρr = Convex.Variable(d,d)
    ρi = Convex.Variable(d,d)

    constraints = Convex.tr(ρr) == 1
    constraints += Convex.tr(ρi) == 0
    constraints += Convex.isposdef([ρr ρi; -ρi ρr])

    # TODO: use quad_form instead of vecnorm? Have 1/vars are diagonal quadratic form
    problem = Convex.minimize( LinearAlgebra.norm( (means - method.realpred*[vec(ρr); vec(ρi)]) .* ivars, 2), constraints )

    Convex.solve!(problem, solver)

    return (ρr.value - 1im*ρi.value), problem.optval^2, problem.status
end

"""
`MLStateTomo`

Maximum-likelihood maximum-entropy quantum state tomography algorithm. It is
constructed from a collection of observables corresponding to
measurements that are performed on the state being reconstructed, as
well as a hedging factor β. If β=0, no hedging is applied. If β > 0
either a log determinant or log minimum eigenvalue penalty is applied.

**AT THE MOMENT HEDGING IS ONLY APPLIED IN fitA, AND IS NOT CURRENTLY WORKING**
"""
struct MLStateTomo
    effects::Vector{Matrix{ComplexF64}}
    dim::Int64
    β::Float64
    # TODO: perhaps for ML it is better to have the observables pecified as
    #       a vector of vectors of observables? Each vector of observable is a POVM
    #       that should add up to the identity
    function MLStateTomo(v::Vector,β=0.0)
        for e in v
            if !LinearAlgebra.ishermitian(e) || ! QuantumInfo.ispossemidef(e) || real(LinearAlgebra.tr(e))>1
                error("MLStateTomo state tomography is parameterized by POVM effects only.")
            end
        end
        sv = sum(v)
        if ! LinearAlgebra.isdiag(sv) || !isapprox(maximum(abs.(LinearAlgebra.diag(sv)))-minimum(abs.(LinearAlgebra.diag(sv))),0.0)
            error("POVM effects must add up to the identity.")
        end
        if !all([size(e,1)==size(e,2) for e in v]) || !all([size(v[1],1)==size(e,1) for e in v])
                error("All effects must be square matrices, and they must have have the same dimension.")
        end
        if β < 0
            error("Hedging penalty must be positive.")
        end
        new(v,size(v[1],1),β)
    end
end

function predict(method::MLStateTomo,ρ)
    Float64[ real(LinearAlgebra.tr(ρ*e)) for e in method.effects]
end

"""
`fitA(method::MLStateTomo, freq, solver=SCSSolver)`

This is a state tomography fitting routine using convex optimization. It's use is currently
discouraged simply because it is much slower than the iterative solver, and often does not
conver to a solution that meets the optimality criteria.
"""
function fitA(method::MLStateTomo,
             freq::Vector;
             #solver = MosekSolver(LOG=0))
             solver = SCSSolver(verbose=0, max_iters=100_000, eps = 1e-8))

    if length(method.effects) != length(freq)
        error("Vector of counts and vector of effects must have same length, but length(counts) == $(length(counts)) != $(length(method.effects))")
    end

    d = method.dim

    # ρr = Convex.Variable(d,d)
    ρr = Semidefinite(d)
    ρi = Convex.Variable(d,d)

    ϕ(m) = [real(m) -imag(m); imag(m) real(m)];
    ϕ(r,i) = [r i; -i r]

    obj = freq[1] * log(Convex.tr([ρr ρi; -ρi ρr]*ϕ(method.effects[1])))
    for i=2:length(method.effects)
        obj += freq[i] * log(Convex.tr([ρr ρi; -ρi ρr]*ϕ(method.effects[i])))
    end
    # add hedging regularization
    if method.β != 0.0
        obj += method.β*log(lambdamin([ρr ρi; -ρi ρr]))
    end

    constraints = Convex.tr(ρr) == 1
    constraints += Convex.tr(ρi) == 0
    constraints += Convex.isposdef([ρr ρi; -ρi ρr])

    problem = Convex.maximize(obj, constraints)

    Convex.solve!(problem, solver)

    return ρr.value-im*ρi.value, problem.optval, problem.status
end

"""
`fitB(method::MLStateTomo, freq, solver=SCSSolver)`

This is a state tomography fitting routine using convex optimization. It's use is currently
discouraged simply because it is much slower than the iterative solver, and often does not
conver to a solution that meets the optimality criteria, much like fitB.
"""
function fitB(method::MLStateTomo,
             freq::Vector;
             #solver = MosekSolver(LOG=0))
             solver = SCSSolver(verbose=0, max_iters=100_000, eps = 1e-8))

    if length(method.effects) != length(freq)
        error("Vector of counts and vector of effects must have same length, but length(counts) == $(length(counts)) != $(length(method.effects))")
    end

    d = method.dim

    ρr = Convex.Semidefinite(d)
    ρi = Convex.Variable(d,d)
    p  = Convex.Variable(length(method.effects),Positive())

    ϕ(m) = [real(m) -imag(m); imag(m) real(m)];
    ϕ(r,i) = [r i; -i r]
    obj = freq[1] * log(p[1])
    for i=2:length(method.effects)
        obj += freq[i] * log(p[i])
    end
    #obj += method.β!=0.0 ? method.β*logdet([ρr ρi; -ρi ρr]) : 0.0
    #obj += method.β!=0.0 ? method.β*log(lambdamin([ρr ρi; -ρi ρr])) : 0.0

    constraints = Convex.tr(ρr) == 1
    constraints += Convex.tr(ρi) == 0
    for i=1:length(method.effects)
        constraints += p[i] == Convex.tr([ρr ρi; -ρi ρr]*ϕ(method.effects[i]))/2
    end
    constraints += Convex.isposdef([ρr ρi; -ρi ρr])

    problem = Convex.maximize(obj, constraints)

    Convex.solve!(problem, solver)

    return ρr.value-im*ρi.value, problem.optval, problem.status
end

function R(ρ,E,obs)
    pr = Float64[real(LinearAlgebra.tr(ρ*Ei)) for Ei in E]
    R = obs[1]/pr[1]*E[1]
    for i in 2:length(E)
        R += obs[i]/pr[i]*E[i]
    end
    R
end

function LL(ρ,E,obs)
    pr = Float64[real(LinearAlgebra.tr(ρ*Ei)) for Ei in E]
    sum(obs.*log.(pr))
end

function fit(method::MLStateTomo,
             freq::Vector;
             δ=.005, # a more natural parameterization of "dilution"
             λ=0.0,  # entropy penalty
             tol=1e-9,
             maxiter=100_000,
             ρ0 = QuantumInfo.eye(ComplexF64,method.dim))
    ϵ=1/δ
    iter = 1
    ρk = copy(ρ0)
    ρktemp = similar(ρk)
    ρkm = similar(ρ0) #Array(ComplexF64,method.dim,method.dim)
    status = :Optimal
    while true
        copyto!(ρkm,ρk)
        if iter >= maxiter
            status = :MaxIter
            break
        end
        # diluted iterative scheme for ML (from PRA 75 042108 2007)
        Tk = λ > 0.0 ?
            (1 .+ ϵ * (R(ρk,method.effects,freq) .- QuantumInfo.eye(method.dim) .-
                       λ * (log(ρk) .- LinearAlgebra.tr(ρk*log(ρk)))))/(1+ϵ) :
                       (1 .+ ϵ * R(ρk,method.effects,freq)) / (1 + ϵ)
        ρktemp .= ρk * (1 .+ ϵ * Tk)/( 1 + ϵ)
        ρk .= (1 .+ ϵ * Tk) / (1 + ϵ) * ρktemp
        # LinearAlgebra.A_mul_B!(ρktemp, ρk, (1 .+ ϵ * Tk)/( 1 + ϵ)) # ρk = ρk * (1+ϵRk)/(1+ϵ)
        # LinearAlgebra.A_mul_B!(ρk, (1+ϵ*Tk) / (1+ϵ), ρktemp) # ρk = (1+ϵRk) * ρk/(1+ϵ)
        QuantumInfo.trnormalize!(ρk)
        if LinearAlgebra.norm(ρk-ρkm)/method.dim^2 < tol
            status = :Optimal
            break
        end
        iter += 1
    end
    return ρk, LL(ρk,method.effects,freq), status
end


text:
using Pkg
Pkg.activate(".")
using DeutschlandDaten
using Dates
using DataFrames
using DataFramesMeta
using StatsBase

deathdat = DeutschlandDaten.deaths_data()
popdat  = DeutschlandDaten.population()
popdatbl = DeutschlandDaten.population_bundesländer()

splitticks(x) = ( [ e[1] for e in x], [ e[2] for e in x])

mordat = [
    (
        alter = a,
        geschlecht=g,
        jahr=j,
        monat = monat,
        N=convert(Int, round(population(popdat; alter=a, geschlecht=g, jahr=j))),
        D=deaths(deathdat, alter=a, geschlecht=g, jahr=j, monat=monat)
    )
        for a in DeutschlandDaten.Altersgruppen_Monate
            for g in ["Männlich", "Weiblich"]
                for j in 2000:2021
                    for monat in 1:12
                        ] |> DataFrame;

mordat[:,:PD] = 10000 * mordat.D ./ mordat.N
mordat[:,:col] = collect(zip(mordat.jahr, mordat.geschlecht))
mordat[:,:row] = mordat.monat

function gliding_mean(x, window=12)
    y = Vector{Union{Missing,Float64}}(missing, length(x))
    for i in 1:length(x)
        y[i] = mean(skipmissing(x[max(1,i-div(window, 2)):min(end,i+window-div(window,2))]))
    end
    y
end
gliding_mean(mordat.PD)

@transform!(groupby(mordat, [:alter, :geschlecht]),
            :PD_12months = gliding_mean(:PD),
            :date = Date.(:jahr, :monat, 1)
            )

using StatsPlots
using Plots.PlotMeasures

savefig(plotyears(mordat), "images/Sterberaten.svg")

savefig(plotyears(filter(x-> x.date>=Date(2016,1,1), mordat)), "images/Sterberaten_2016_2021.svg")
for alter=DeutschlandDaten.Altersgruppen_Monate
    d = filter(x->x.alter==alter, mordat)
    us = unstack(d[:,[:row,:col,:PD]], :col, :PD)
    p = heatmap((Matrix(us[:,2:end]));
                title = "DE: Monatliche Sterberate, $alter alt, von 10'000",
                yticks = 1:12,
                xticks = (vcat(1:5:22, 23:5:44), repeat(collect(string.(2000:5:2021)),2)),
                xrotation = -90,
                c = :reds )
    plot!([22.5]; seriestype=:vline, color=:black, legend=nothing, width=2)
    annotate!([(12,10,"Männlich"), (23+12,10,"Weiblich")])
    savefig(p, "images/Sterberate_Monat_$(alter).svg")
end

mordat = [
    (
        alter = a,
        geschlecht=g,
        jahr=j,
        kw = kw,
        # date = Date(j, monat, 1),
        # N=population(alter=a, geschlecht=g, jahr=j),
        N=convert(Int, round(population(popdat; alter=a, geschlecht=g, jahr=j))),
        D=deaths(deathdat, alter=a, geschlecht=g, jahr=j, kw=kw)
    )
        for a in DeutschlandDaten.Altersgruppen_KW
            for g in ["Männlich", "Weiblich"]
                for j in 2000:2021
                    for kw in 1:52
                        ] |> DataFrame;

mordat[:,:PD] = 10000 * mordat.D ./ mordat.N
mordat[:,:col] = collect(zip(mordat.jahr, mordat.geschlecht))
mordat[:,:row] = mordat.kw

alter=DeutschlandDaten.Altersgruppen_KW|>first
for alter=DeutschlandDaten.Altersgruppen_KW
    d = filter(x->x.alter==alter, mordat)
    us = unstack(d[:,[:row,:col,:PD]], :col, :PD)
    p = heatmap((Matrix(us[:,2:end]));
                title = "Sterberate Deutschland im Alter $alter von 10'000",
                yticks = 1:5:52,
                xticks = (vcat(1:5:22, 23:5:44), repeat(collect(string.(2000:5:2021)),2)),
                xrotation = -90,
                c = :reds )
    plot!([22.5]; seriestype=:vline, color=:black, legend=nothing, width=2)
    annotate!([(12,45,"Männlich"), (23+12,45,"Weiblich")])
    savefig(p, "images/Sterberate_$(alter).svg")
end

mordat = ([
    (
        bundesland = bl,
        alter = a,
        geschlecht=g,
        jahr=j,
        monat = m,
        date = Date(j, m, 28),
        N = convert(Int, round(population(popdatbl, deathdat; bundesland = bl, alter=a, geschlecht=g, jahr=j))),
        D = deaths(deathdat; bundesland = bl, alter=a, geschlecht=g, jahr=j, monat=m)
    )
    for bl in DeutschlandDaten.Bundesländer
        for a in DeutschlandDaten.Altersgruppen_Bundesländer
            for g in ["Männlich", "Weiblich"]
                for j in 2000:2021
                    for m in 1:12
                        ]) |> DataFrame;

@transform!(mordat,
            :PD = 10000 .* :D ./ :N,
            :col = collect(zip(:jahr, :geschlecht)),
            :row = collect(zip(:bundesland, :monat)))

function partifact(x, fact=10)
    sp = sortperm(x)
    r = Vector{Int}[[sp[1]]]
    for i in sp[2:end]
        if x[i] > x[r[end][1]] * fact
            push!(r,[i])
        else
            push!(r[end],i)
        end
    end
    r
end

alter = last(DeutschlandDaten.Altersgruppen_Bundesländer)
for (alter) = (DeutschlandDaten.Altersgruppen_Bundesländer)
    mean_PD = combine(groupby(filter(x->x.alter ==  alter, mordat), :bundesland),
                      :PD => x->mean(skipmissing(x)))
    for (i, bli) = enumerate(partifact(mean_PD[:,2],2))
        bli = partifact(mean_PD[:,2],2)[end]
        bls = mean_PD.bundesland[bli]
        d = filter(x->x.alter ==  alter && x.bundesland in bls, mordat)
        us = unstack(d[:,[:row,:col,:PD]], :col, :PD)
        p = heatmap((Matrix(us[:,2:end]));
                    title = "Sterberate Bundesländer, $alter alt, von 10'000",
                    yticks = splitticks([ (i, x[1]) for (i,x) in enumerate(us[:,1]) if x[2]==6 ]),
                    xticks = (vcat(1:5:22, 23:5:44), repeat(collect(string.(2000:5:2021)),2)),
                    xrotation = -90,
                    c = :reds
                    )
        plot!([22.5]; seriestype=:vline, color=:black, legend=nothing, width=2)
        plot!((0:12:(size(us)[1])) .+ 0.5; seriestype=:hline, color=:black, legend=nothing, width=2)
        annotate!([(12,size(us)[1]-5,"Männlich"), (23+12,size(us)[1]-5,"Weiblich")])
        savefig(p, "images/Sterberate_$(alter)_$i.svg")
        println("$(alter)_$i: ", bls)
    end
end

max_stars_repo_name: gkappler/DeutschlandDaten

max_issues_repo_name: gkappler/DeutschlandDaten

max_forks_repo_name: gkappler/DeutschlandDaten


text:
@testset "optimize" begin
    # Tests for PR #302
    results = optimize(cos, 0, 2pi);
    @test norm(Optim.minimizer(results) - pi) < 0.01
    results = optimize(cos, 0.0, 2pi);
    @test norm(Optim.minimizer(results) - pi) < 0.01
    results = optimize(cos, 0, 2pi, Brent());
    @test norm(Optim.minimizer(results) - pi) < 0.01
    results = optimize(cos, 0.0, 2pi, Brent());
    @test norm(Optim.minimizer(results) - pi) < 0.01
    results = optimize(cos, 0, 2pi, method = Brent());
    @test norm(Optim.minimizer(results) - pi) < 0.01
    results = optimize(cos, 0.0, 2pi, method = Brent());
    @test norm(Optim.minimizer(results) - pi) < 0.01
end


text:
using MemoryMutate

# runs on my machine: archlinux x86_64

include("testutils.jl")

#= we have
- Julep: setfield! for mutable references to immutables
  https://github.com/JuliaLang/julia/issues/17115
  - "we propose making it possible to have setfield! modify fields inside of immutable objects that are wrapped in mutable objects"
  - "To support this proposal, the setfield! function will get a multi-arg form, with the following behaviors:"
    "setfield!(x, a, b, c, value) mutates the right most mutable object to change the value of its fields to be"
    "equivalent to copying the immutable objects and updating the referenced field."
  - "tl;dr The syntax:"
    "x.a.b.c = 3"
    "would now be valid, as long as at least one of the referenced fields is mutable."
- WIP: Make mutating immutables easier
  https://github.com/JuliaLang/julia/pull/21912
  - proposes an @ operator as in `x@a = 2` for an immutable `x`
  - "The ways this works is that under the hood, it creates a new immutable object with the specified field modified and then assigns it back to the appropriate place."
    "Syntax wise, everything to the left of the @ is what's being assigned to, everything to the right of the @ is what is to be modified."
- Make immutable mutable again - or make arrays of structures stack-allocated
  https://discourse.julialang.org/t/make-immutable-mutable-again-or-make-arrays-of-structures-stack-allocated/20800
  - colleced a few related links
=#

################################################################################

        struct G1; x :: Float32;          end # struct G2 { const float x; };
        struct F1; x :: Float32; g :: G1; end # struct F2 { const float x; const G2 g; };
        struct E1; x :: Float32; f :: F1; end # struct E2 { const float x; const F2 g; };
        struct D1; x :: Float32; e :: E1; end # struct D2 { const float x; const E2 g; };
        struct C1; x :: Float32; d :: D1; end # struct C2 { const float x; const D2 g; };
        struct B1; x :: Float32; c :: C1; end # struct B2 { const float x; const C2 g; };
mutable struct A1; x :: Float32; b :: B1; end # struct A2 {       float x;       B2 g; };

g1 = G1(0.0f0)
f1 = F1(0.0f0,g1)
e1 = E1(0.0f0,f1)
d1 = D1(0.0f0,e1)
c1 = C1(0.0f0,d1)
b1 = B1(0.0f0,c1)
a1 = A1(0.0f0,b1)

test101(a::A1,v::Float32) = @mem a.b.c.d.e.f.g.x = v
test102(a::A1,v::G1)      = @mem a.b.c.d.e.f.g = v
test103(a::A1,v::Float32) = @mem a.b.c.d.e.f.x = v
test104(a::A1,v::F1)      = @mem a.b.c.d.e.f = v
test105(a::A1,v::Float32) = @mem a.b.c.d.e.x = v
test106(a::A1,v::E1)      = @mem a.b.c.d.e = v
test107(a::A1,v::Float32) = @mem a.b.c.d.x = v
test108(a::A1,v::D1)      = @mem a.b.c.d = v
test109(a::A1,v::Float32) = @mem a.b.c.x = v
test110(a::A1,v::C1)      = @mem a.b.c = v
test111(a::A1,v::Float32) = @mem a.b.x = v
test112(a::A1,v::B1)      = @mem a.b = v
test113(a::A1,v::Float32) = @mem a.x = v
test114(a::A1,v::A1)      = @mem a = v

v  = 0f0
v += 1f0; @assert a1.b.c.d.e.f.g.x != v;        test101(a1,v);        @assert a1.b.c.d.e.f.g.x == v
v += 1f0; @assert a1.b.c.d.e.f.g   != G1(v);    test102(a1,G1(v));    @assert a1.b.c.d.e.f.g   == G1(v)
v += 1f0; @assert a1.b.c.d.e.f.x   != v;        test103(a1,v);        @assert a1.b.c.d.e.f.x   == v
v += 1f0; @assert a1.b.c.d.e.f     != F1(v,g1); test104(a1,F1(v,g1)); @assert a1.b.c.d.e.f     == F1(v,g1)
v += 1f0; @assert a1.b.c.d.e.x     != v;        test105(a1,v);        @assert a1.b.c.d.e.x     == v
v += 1f0; @assert a1.b.c.d.e       != E1(v,f1); test106(a1,E1(v,f1)); @assert a1.b.c.d.e       == E1(v,f1)
v += 1f0; @assert a1.b.c.d.x       != v;        test107(a1,v);        @assert a1.b.c.d.x       == v
v += 1f0; @assert a1.b.c.d         != D1(v,e1); test108(a1,D1(v,e1)); @assert a1.b.c.d         == D1(v,e1)
v += 1f0; @assert a1.b.c.x         != v;        test109(a1,v);        @assert a1.b.c.x         == v
v += 1f0; @assert a1.b.c           != C1(v,d1); test110(a1,C1(v,d1)); @assert a1.b.c           == C1(v,d1)
v += 1f0; @assert a1.b.x           != v;        test111(a1,v);        @assert a1.b.x           == v
v += 1f0; @assert a1.b             != B1(v,c1); test112(a1,B1(v,c1)); @assert a1.b             == B1(v,c1)
v += 1f0; @assert a1.x             != v;        test113(a1,v);        @assert a1.x             == v
v += 1f0; @assert a1               != A1(v,b1); test114(a1,A1(v,b1)); # @assert a1               == A1(v,b1) # call by "sharing"

code_native(io,test101,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test102,(A1,G1))     ; display_asm_stat_io(io) # (total = 4, movs = 3, mov = 2, vmov = 1)
code_native(io,test103,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test104,(A1,F1))     ; display_asm_stat_io(io) # (total = 5, movs = 4, mov = 2, vmov = 2)
code_native(io,test105,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test106,(A1,E1))     ; display_asm_stat_io(io) # (total = 11, movs = 9, mov = 9, vmov = 0)
code_native(io,test107,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test108,(A1,D1))     ; display_asm_stat_io(io) # (total = 6, movs = 5, mov = 1, vmov = 4)
code_native(io,test109,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test110,(A1,C1))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test111,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test112,(A1,B1))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test113,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test114,(A1,A1))     ; display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

################################################################################

a1 = A1(0.0f0,b1)

test101b(a::A1,v::Float32) = @mem getfield(a.b.c,:d).e.f.g.x = v
test102b(a::A1,v::G1)      = @mem getfield(a.b.c,:d).e.f.g = v
test103b(a::A1,v::Float32) = @mem getfield(a.b.c,:d).e.f.x = v
test104b(a::A1,v::F1)      = @mem getfield(a.b.c,:d).e.f = v
test105b(a::A1,v::Float32) = @mem getfield(a.b.c,:d).e.x = v
test106b(a::A1,v::E1)      = @mem getfield(a.b.c,:d).e = v
test107b(a::A1,v::Float32) = @mem getfield(a.b.c,:d).x = v
test108b(a::A1,v::D1)      = @mem getfield(a.b.c,:d) = v
test109b(a::A1,v::Float32) = @mem a.b.c.x = v
test110b(a::A1,v::C1)      = @mem a.b.c = v
test111b(a::A1,v::Float32) = @mem a.b.x = v
test112b(a::A1,v::B1)      = @mem a.b = v
test113b(a::A1,v::Float32) = @mem a.x = v
test114b(a::A1,v::A1)      = @mem a = v

v  = 0f0
v += 1f0; @assert a1.b.c.d.e.f.g.x != v;        test101b(a1,v);        @assert a1.b.c.d.e.f.g.x == v
v += 1f0; @assert a1.b.c.d.e.f.g   != G1(v);    test102b(a1,G1(v));    @assert a1.b.c.d.e.f.g   == G1(v)
v += 1f0; @assert a1.b.c.d.e.f.x   != v;        test103b(a1,v);        @assert a1.b.c.d.e.f.x   == v
v += 1f0; @assert a1.b.c.d.e.f     != F1(v,g1); test104b(a1,F1(v,g1)); @assert a1.b.c.d.e.f     == F1(v,g1)
v += 1f0; @assert a1.b.c.d.e.x     != v;        test105b(a1,v);        @assert a1.b.c.d.e.x     == v
v += 1f0; @assert a1.b.c.d.e       != E1(v,f1); test106b(a1,E1(v,f1)); @assert a1.b.c.d.e       == E1(v,f1)
v += 1f0; @assert a1.b.c.d.x       != v;        test107b(a1,v);        @assert a1.b.c.d.x       == v
v += 1f0; @assert a1.b.c.d         != D1(v,e1); test108b(a1,D1(v,e1)); @assert a1.b.c.d         == D1(v,e1)
v += 1f0; @assert a1.b.c.x         != v;        test109b(a1,v);        @assert a1.b.c.x         == v
v += 1f0; @assert a1.b.c           != C1(v,d1); test110b(a1,C1(v,d1)); @assert a1.b.c           == C1(v,d1)
v += 1f0; @assert a1.b.x           != v;        test111b(a1,v);        @assert a1.b.x           == v
v += 1f0; @assert a1.b             != B1(v,c1); test112b(a1,B1(v,c1)); @assert a1.b             == B1(v,c1)
v += 1f0; @assert a1.x             != v;        test113b(a1,v);        @assert a1.x             == v
v += 1f0; @assert a1               != A1(v,b1); test114b(a1,A1(v,b1)); # @assert a1               == A1(v,b1) # call by "sharing"

code_native(io,test101b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test102b,(A1,G1))     ; display_asm_stat_io(io) # (total = 4, movs = 3, mov = 2, vmov = 1)
code_native(io,test103b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test104b,(A1,F1))     ; display_asm_stat_io(io) # (total = 5, movs = 4, mov = 2, vmov = 2)
code_native(io,test105b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test106b,(A1,E1))     ; display_asm_stat_io(io) # (total = 11, movs = 9, mov = 9, vmov = 0)
code_native(io,test107b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test108b,(A1,D1))     ; display_asm_stat_io(io) # (total = 6, movs = 5, mov = 1, vmov = 4)
code_native(io,test109b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test110b,(A1,C1))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test111b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test112b,(A1,B1))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test113b,(A1,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test114b,(A1,A1))     ; display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

################################################################################

        struct G2; x :: Float32;          end # struct G2 { const float x; };
        struct F2; x :: Float32; g :: G2; end # struct F2 { const float x; const G2  g; };
        struct E2; x :: Float32; f :: F2; end # struct E2 { const float x; const F2  g; };
mutable struct D2; x :: Float32; e :: E2; end # struct D2 {       float x;       E2  g; };
        struct C2; x :: Float32; d :: D2; end # struct C2 { const float x; const D2& g; };
        struct B2; x :: Float32; c :: C2; end # struct B2 { const float x; const C2& g; };
mutable struct A2; x :: Float32; b :: B2; end # struct A2 {       float x;       B2& g; };

g2 = G2(6.0f0)
f2 = F2(5.0f0,g2)
e2 = E2(4.0f0,f2)
d2 = D2(3.0f0,e2)
c2 = C2(2.0f0,d2)
b2 = B2(1.0f0,c2)
a2 = A2(0.0f0,b2)

test201(a::A2,v::Float32) = @mem  a.b.c.d.e.f.g.x = v
test202(a::A2,v::G2)      = @mem  a.b.c.d.e.f.g = v
test203(a::A2,v::Float32) = @mem  a.b.c.d.e.f.x = v
test204(a::A2,v::F2)      = @mem  a.b.c.d.e.f = v
test205(a::A2,v::Float32) = @mem  a.b.c.d.e.x = v
test206(a::A2,v::E2)      = @mem  a.b.c.d.e = v
test207(a::A2,v::Float32) = @mem  a.b.c.d.x = v
test208(a::A2,v::D2)      = @yolo a.b.c.d = v
test209(a::A2,v::Float32) = @yolo a.b.c.x = v
test210(a::A2,v::C2)      = @yolo a.b.c = v
test211(a::A2,v::Float32) = @yolo a.b.x = v
test212(a::A2,v::B2)      = @mem  a.b = v
test213(a::A2,v::Float32) = @mem  a.x = v
test214(a::A2,v::A2)      = @mem  a = v

            v  = 0f0
            v += 1f0; @assert a2.b.c.d.e.f.g.x != v;        test201(a2,v);        @assert a2.b.c.d.e.f.g.x == v
            v += 1f0; @assert a2.b.c.d.e.f.g   != G2(v);    test202(a2,G2(v));    @assert a2.b.c.d.e.f.g   == G2(v)
            v += 1f0; @assert a2.b.c.d.e.f.x   != v;        test203(a2,v);        @assert a2.b.c.d.e.f.x   == v
            v += 1f0; @assert a2.b.c.d.e.f     != F2(v,g2); test204(a2,F2(v,g2)); @assert a2.b.c.d.e.f     == F2(v,g2)
            v += 1f0; @assert a2.b.c.d.e.x     != v;        test205(a2,v);        @assert a2.b.c.d.e.x     == v
            v += 1f0; @assert a2.b.c.d.e       != E2(v,f2); test206(a2,E2(v,f2)); @assert a2.b.c.d.e       == E2(v,f2)
            v += 1f0; @assert a2.b.c.d.x       != v;        test207(a2,v);        @assert a2.b.c.d.x       == v
try; global v += 1f0; @assert a2.b.c.d         != D2(v,e2); test208(a2,D2(v,e2)); @assert a2.b.c.d.x       == D2(v,e2).x # mutables with equal field values are not equal
                                                                                @assert a2.b.c.d.e       == D2(v,e2).e; catch e; print(e); end
try; global v += 1f0; @assert a2.b.c.x         != v;        test209(a2,v);        @assert a2.b.c.x         == v         ; catch e; print(e); end
try; global v += 1f0; @assert a2.b.c           != C2(v,d2); test210(a2,C2(v,d2)); @assert a2.b.c           == C2(v,d2)  ; catch e; print(e); end
try; global v += 1f0; @assert a2.b.x           != v;        test211(a2,v);        @assert a2.b.x           == v         ; catch e; print(e); end
            v += 1f0; @assert a2.b             != B2(v,c2); test212(a2,B2(v,c2)); @assert a2.b             == B2(v,c2)
            v += 1f0; @assert a2.x             != v;        test213(a2,v);        @assert a2.x             == v
            v += 1f0; @assert a2               != A2(v,b2); test214(a2,A2(v,b2)); # @assert a1               == A2(v,b2) # call by "sharing"

code_native(io,test201,(A2,Float32)); display_asm_stat_io(io) # (total = 5, movs = 4, mov = 3, vmov = 1)
code_native(io,test202,(A2,G2))     ; display_asm_stat_io(io) # (total = 7, movs = 6, mov = 5, vmov = 1)
code_native(io,test203,(A2,Float32)); display_asm_stat_io(io) # (total = 5, movs = 4, mov = 3, vmov = 1)
code_native(io,test204,(A2,F2))     ; display_asm_stat_io(io) # (total = 9, movs = 7, mov = 5, vmov = 2)
code_native(io,test205,(A2,Float32)); display_asm_stat_io(io) # (total = 5, movs = 4, mov = 3, vmov = 1)
code_native(io,test206,(A2,E2))     ; display_asm_stat_io(io) # (total = 14, movs = 12, mov = 12, vmov = 0)
code_native(io,test207,(A2,Float32)); display_asm_stat_io(io) # (total = 5, movs = 4, mov = 3, vmov = 1)
code_native(io,test208,(A2,D2))     ; display_asm_stat_io(io) # (total = 7, movs = 6, mov = 6, vmov = 0)
code_native(io,test209,(A2,Float32)); display_asm_stat_io(io) # (total = 28, movs = 18, mov = 15, vmov = 3)
code_native(io,test210,(A2,C2))     ; display_asm_stat_io(io) # (total = 18, movs = 15, mov = 14, vmov = 1)
code_native(io,test211,(A2,Float32)); display_asm_stat_io(io) # (total = 26, movs = 17, mov = 14, vmov = 3)
code_native(io,test212,(A2,B2))     ; display_asm_stat_io(io) # (total = 24, movs = 8, mov = 8, vmov = 0)
code_native(io,test213,(A2,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test214,(A2,A2))     ; display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

################################################################################

        struct G3; x :: Float32;          end # struct G3 { const float x; };
        struct F3; x :: Float32; g :: G3; end # struct F3 { const float x; const G3 g; };
        struct E3; x :: Float32; f :: F3; end # struct E3 { const float x; const F3 g; };
        struct D3; x :: Float32; e :: E3; end # struct D3 { const float x; const E3 g; };
        struct C3; x :: Float32; d :: D3; end # struct C3 { const float x; const D3 g; };
        struct B3; x :: Float32; c :: C3; end # struct B3 { const float x; const C3 g; };
        struct A3; x :: Float32; b :: B3; end # struct A3 { const float x; const B3 g; };

g3 = G3(6.0f0)
f3 = F3(5.0f0,g3)
e3 = E3(4.0f0,f3)
d3 = D3(3.0f0,e3)
c3 = C3(2.0f0,d3)
b3 = B3(1.0f0,c3)
a3 = A3(0.0f0,b3)

test301(a::Base.RefValue{A3},v::Float32) = @mem a[].b.c.d.e.f.g.x = v
test302(a::Base.RefValue{A3},v::G3)      = @mem a[].b.c.d.e.f.g = v
test303(a::Base.RefValue{A3},v::Float32) = @mem a[].b.c.d.e.f.x = v
test304(a::Base.RefValue{A3},v::F3)      = @mem a[].b.c.d.e.f = v
test305(a::Base.RefValue{A3},v::Float32) = @mem a[].b.c.d.e.x = v
test306(a::Base.RefValue{A3},v::E3)      = @mem a[].b.c.d.e = v
test307(a::Base.RefValue{A3},v::Float32) = @mem a[].b.c.d.x = v
test308(a::Base.RefValue{A3},v::D3)      = @mem a[].b.c.d = v
test309(a::Base.RefValue{A3},v::Float32) = @mem a[].b.c.x = v
test310(a::Base.RefValue{A3},v::C3)      = @mem a[].b.c = v
test311(a::Base.RefValue{A3},v::Float32) = @mem a[].b.x = v
test312(a::Base.RefValue{A3},v::B3)      = @mem a[].b = v
test313(a::Base.RefValue{A3},v::Float32) = @mem a[].x = v
test314(a::Base.RefValue{A3},v::A3)      = @mem a[] = v

a3r = Ref(a3)

v = 0f0
v += 1f0; @assert a3r[].b.c.d.e.f.g.x != v;        test301(a3r,v);        @assert a3r[].b.c.d.e.f.g.x == v
v += 1f0; @assert a3r[].b.c.d.e.f.g   != G3(v);    test302(a3r,G3(v));    @assert a3r[].b.c.d.e.f.g   == G3(v)
v += 1f0; @assert a3r[].b.c.d.e.f.x   != v;        test303(a3r,v);        @assert a3r[].b.c.d.e.f.x   == v
v += 1f0; @assert a3r[].b.c.d.e.f     != F3(v,g3); test304(a3r,F3(v,g3)); @assert a3r[].b.c.d.e.f     == F3(v,g3)
v += 1f0; @assert a3r[].b.c.d.e.x     != v;        test305(a3r,v);        @assert a3r[].b.c.d.e.x     == v
v += 1f0; @assert a3r[].b.c.d.e       != E3(v,f3); test306(a3r,E3(v,f3)); @assert a3r[].b.c.d.e       == E3(v,f3)
v += 1f0; @assert a3r[].b.c.d.x       != v;        test307(a3r,v);        @assert a3r[].b.c.d.x       == v
v += 1f0; @assert a3r[].b.c.d         != D3(v,e3); test308(a3r,D3(v,e3)); @assert a3r[].b.c.d         == D3(v,e3)
v += 1f0; @assert a3r[].b.c.x         != v;        test309(a3r,v);        @assert a3r[].b.c.x         == v
v += 1f0; @assert a3r[].b.c           != C3(v,d3); test310(a3r,C3(v,d3)); @assert a3r[].b.c           == C3(v,d3)
v += 1f0; @assert a3r[].b.x           != v;        test311(a3r,v);        @assert a3r[].b.x           == v
v += 1f0; @assert a3r[].b             != B3(v,c3); test312(a3r,B3(v,c3)); @assert a3r[].b             == B3(v,c3)
v += 1f0; @assert a3r[].x             != v;        test313(a3r,v);        @assert a3r[].x             == v
v += 1f0; @assert a3r[]               != A3(v,b3); test314(a3r,A3(v,b3)); @assert a3r[]               == A3(v,b3)

code_native(io,test301,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test302,(Base.RefValue{A3},G3))     ; display_asm_stat_io(io) # (total = 4, movs = 3, mov = 2, vmov = 1)
code_native(io,test303,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test304,(Base.RefValue{A3},F3))     ; display_asm_stat_io(io) # (total = 5, movs = 4, mov = 2, vmov = 2)
code_native(io,test305,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test306,(Base.RefValue{A3},E3))     ; display_asm_stat_io(io) # (total = 11, movs = 9, mov = 9, vmov = 0)
code_native(io,test307,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test308,(Base.RefValue{A3},D3))     ; display_asm_stat_io(io) # (total = 6, movs = 5, mov = 1, vmov = 4)
code_native(io,test309,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test310,(Base.RefValue{A3},C3))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test311,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test312,(Base.RefValue{A3},B3))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test313,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test314,(Base.RefValue{A3},A3))     ; display_asm_stat_io(io) # (total = 11, movs = 9, mov = 1, vmov = 8)

################################################################################

test301b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].b.c.d.e.f.g.x = v
test302b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::G3)      = @mem a[][][].b.c.d.e.f.g = v
test303b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].b.c.d.e.f.x = v
test304b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::F3)      = @mem a[][][].b.c.d.e.f = v
test305b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].b.c.d.e.x = v
test306b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::E3)      = @mem a[][][].b.c.d.e = v
test307b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].b.c.d.x = v
test308b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::D3)      = @mem a[][][].b.c.d = v
test309b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].b.c.x = v
test310b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::C3)      = @mem a[][][].b.c = v
test311b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].b.x = v
test312b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::B3)      = @mem a[][][].b = v
test313b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::Float32) = @mem a[][][].x = v
test314b(a::Base.RefValue{Base.RefValue{Base.RefValue{A3}}},v::A3)      = @mem a[][][] = v

a3rrr = Ref(Ref(Ref(a3)))

# @assert false

 v = 0f0
v += 1f0; @assert a3rrr[][][].b.c.d.e.f.g.x != v;        test301b(a3rrr,v);        @assert a3rrr[][][].b.c.d.e.f.g.x == v
v += 1f0; @assert a3rrr[][][].b.c.d.e.f.g   != G3(v);    test302b(a3rrr,G3(v));    @assert a3rrr[][][].b.c.d.e.f.g   == G3(v)
v += 1f0; @assert a3rrr[][][].b.c.d.e.f.x   != v;        test303b(a3rrr,v);        @assert a3rrr[][][].b.c.d.e.f.x   == v
v += 1f0; @assert a3rrr[][][].b.c.d.e.f     != F3(v,g3); test304b(a3rrr,F3(v,g3)); @assert a3rrr[][][].b.c.d.e.f     == F3(v,g3)
v += 1f0; @assert a3rrr[][][].b.c.d.e.x     != v;        test305b(a3rrr,v);        @assert a3rrr[][][].b.c.d.e.x     == v
v += 1f0; @assert a3rrr[][][].b.c.d.e       != E3(v,f3); test306b(a3rrr,E3(v,f3)); @assert a3rrr[][][].b.c.d.e       == E3(v,f3)
v += 1f0; @assert a3rrr[][][].b.c.d.x       != v;        test307b(a3rrr,v);        @assert a3rrr[][][].b.c.d.x       == v
v += 1f0; @assert a3rrr[][][].b.c.d         != D3(v,e3); test308b(a3rrr,D3(v,e3)); @assert a3rrr[][][].b.c.d         == D3(v,e3)
v += 1f0; @assert a3rrr[][][].b.c.x         != v;        test309b(a3rrr,v);        @assert a3rrr[][][].b.c.x         == v
v += 1f0; @assert a3rrr[][][].b.c           != C3(v,d3); test310b(a3rrr,C3(v,d3)); @assert a3rrr[][][].b.c           == C3(v,d3)
v += 1f0; @assert a3rrr[][][].b.x           != v;        test311b(a3rrr,v);        @assert a3rrr[][][].b.x           == v
v += 1f0; @assert a3rrr[][][].b             != B3(v,c3); test312b(a3rrr,B3(v,c3)); @assert a3rrr[][][].b             == B3(v,c3)
v += 1f0; @assert a3rrr[][][].x             != v;        test313b(a3rrr,v);        @assert a3rrr[][][].x             == v
v += 1f0; @assert a3rrr[][][]               != A3(v,b3); test314b(a3rrr,A3(v,b3)); @assert a3rrr[][][]               == A3(v,b3)

code_native(io,test301,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test302,(Base.RefValue{A3},G3))     ; display_asm_stat_io(io) # (total = 4, movs = 3, mov = 2, vmov = 1)
code_native(io,test303,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test304,(Base.RefValue{A3},F3))     ; display_asm_stat_io(io) # (total = 5, movs = 4, mov = 2, vmov = 2)
code_native(io,test305,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test306,(Base.RefValue{A3},E3))     ; display_asm_stat_io(io) # (total = 11, movs = 9, mov = 9, vmov = 0)
code_native(io,test307,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test308,(Base.RefValue{A3},D3))     ; display_asm_stat_io(io) # (total = 6, movs = 5, mov = 1, vmov = 4)
code_native(io,test309,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test310,(Base.RefValue{A3},C3))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test311,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test312,(Base.RefValue{A3},B3))     ; display_asm_stat_io(io) # (total = 10, movs = 9, mov = 5, vmov = 4)
code_native(io,test313,(Base.RefValue{A3},Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test314,(Base.RefValue{A3},A3))     ; display_asm_stat_io(io) # (total = 11, movs = 9, mov = 1, vmov = 8)

################################################################################

        struct G4; x :: Float32;                         end # struct G4 { const float x; };
        struct F4; x :: Float32; g :: G4;                end # struct F4 { const float x; const G4  g; };
        struct E4; x :: Float32; f :: F4;                end # struct E4 { const float x; const F4  g; };
        struct D4; x :: Float32; e :: Base.RefValue{E4}; end # struct D4 { const float x; const E4& g; };
        struct C4; x :: Float32; d :: D4;                end # struct C4 { const float x; const D4& g; };
        struct B4; x :: Float32; c :: C4;                end # struct B4 { const float x; const C4& g; };
mutable struct A4; x :: Float32; b :: B4;                end # struct A4 {       float x;       B4& g; };

g4 = G4(0.0f0)
f4 = F4(0.0f0,g4)
e4 = E4(0.0f0,f4)
d4 = D4(0.0f0,Ref(e4))
c4 = C4(0.0f0,d4)
b4 = B4(0.0f0,c4)
a4 = A4(0.0f0,b4)

test401(a::A4,v::Float32) = @mem  a.b.c.d.e[].f.g.x = v
test402(a::A4,v::G4)      = @mem  a.b.c.d.e[].f.g = v
test403(a::A4,v::Float32) = @mem  a.b.c.d.e[].f.x = v
test404(a::A4,v::F4)      = @mem  a.b.c.d.e[].f = v
test405(a::A4,v::Float32) = @mem  a.b.c.d.e[].x = v
test406(a::A4,v::E4)      = @mem  a.b.c.d.e[] = v
test407(a::A4,v::Float32) = @yolo a.b.c.d.x = v
test408(a::A4,v::D4)      = @yolo a.b.c.d = v
test409(a::A4,v::Float32) = @yolo a.b.c.x = v
test410(a::A4,v::C4)      = @yolo a.b.c = v
test411(a::A4,v::Float32) = @yolo a.b.x = v
test412(a::A4,v::B4)      = @mem  a.b = v
test413(a::A4,v::Float32) = @mem  a.x = v
test414(a::A4,v::A4)      = @mem  a = v

            v  = 0f0
            v += 1f0; @assert a4.b.c.d.e[].f.g.x != v;             test401(a4,v);             @assert a4.b.c.d.e[].f.g.x == v
            v += 1f0; @assert a4.b.c.d.e[].f.g   != G4(v);         test402(a4,G4(v));         @assert a4.b.c.d.e[].f.g   == G4(v)
            v += 1f0; @assert a4.b.c.d.e[].f.x   != v;             test403(a4,v);             @assert a4.b.c.d.e[].f.x   == v
            v += 1f0; @assert a4.b.c.d.e[].f     != F4(v,g4);      test404(a4,F4(v,g4));      @assert a4.b.c.d.e[].f     == F4(v,g4)
            v += 1f0; @assert a4.b.c.d.e[].x     != v;             test405(a4,v);             @assert a4.b.c.d.e[].x     == v
            v += 1f0; @assert a4.b.c.d.e[]       != E4(v,f4);      test406(a4,E4(v,f4));      @assert a4.b.c.d.e[]       == E4(v,f4)
try; global v += 1f0; @assert a4.b.c.d.x         != v;             test407(a4,v);             @assert a4.b.c.d.x         == v            ; catch e; print(e); end
try; global v += 1f0; @assert a4.b.c.d           != D4(v,Ref(e4)); test408(a4,D4(v,Ref(e4))); @assert a4.b.c.d.e[]       == e4           # references to equal immutables are not equal
                                                                                @assert a4.b.c.d.x         == v            ; catch e; print(e); end
try; global v += 1f0; @assert a4.b.c.x           != v;             test409(a4,v);             @assert a4.b.c.x           == v            ; catch e; print(e); end
try; global v += 1f0; @assert a4.b.c             != C4(v,d4);      test410(a4,C4(v,d4));      @assert a4.b.c             == C4(v,d4)     ; catch e; print(e); end
try; global v += 1f0; @assert a4.b.x             != v;             test411(a4,v);             @assert a4.b.x             == v            ; catch e; print(e); end
            v += 1f0; @assert a4.b               != B4(v,c4);      test412(a4,B4(v,c4));      @assert a4.b               == B4(v,c4)
            v += 1f0; @assert a4.x               != v;             test413(a4,v);             @assert a4.x               == v
            v += 1f0; @assert a4                 != A4(v,b4);      test414(a4,A4(v,b4));      # @assert a4               == A4(v,b4) # call by "sharing"

code_native(io,test401,(A4,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 4, vmov = 1)
code_native(io,test402,(A4,G4))     ; display_asm_stat_io(io) # (total = 8, movs = 7, mov = 6, vmov = 1)
code_native(io,test403,(A4,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 4, vmov = 1)
code_native(io,test404,(A4,F4))     ; display_asm_stat_io(io) # (total = 9, movs = 8, mov = 6, vmov = 2)
code_native(io,test405,(A4,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 4, vmov = 1)
code_native(io,test406,(A4,E4))     ; display_asm_stat_io(io) # (total = 15, movs = 13, mov = 13, vmov = 0)
code_native(io,test407,(A4,Float32)); display_asm_stat_io(io) # (total = 29, movs = 19, mov = 16, vmov = 3)
code_native(io,test408,(A4,D4))     ; display_asm_stat_io(io) # (total = 19, movs = 16, mov = 15, vmov = 1)
code_native(io,test409,(A4,Float32)); display_asm_stat_io(io) # (total = 28, movs = 18, mov = 15, vmov = 3)
code_native(io,test410,(A4,C4))     ; display_asm_stat_io(io) # (total = 18, movs = 15, mov = 14, vmov = 1)
code_native(io,test411,(A4,Float32)); display_asm_stat_io(io) # (total = 26, movs = 17, mov = 14, vmov = 3)
code_native(io,test412,(A4,B4))     ; display_asm_stat_io(io) # (total = 24, movs = 8, mov = 8, vmov = 0)
code_native(io,test413,(A4,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test414,(A4,A4))     ; display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

################################################################################

g4 = G4(0.0f0)
f4 = F4(0.0f0,g4)
e4 = E4(0.0f0,f4)
d4 = D4(0.0f0,Ref(e4))
c4 = C4(0.0f0,d4)
b4 = B4(0.0f0,c4)
a4 = A4(0.0f0,b4)

test401b(a::A4,v::Float32) = @mem  a.b.c.d.e.x.f.g.x = v
test402b(a::A4,v::G4)      = @mem  a.b.c.d.e.x.f.g = v
test403b(a::A4,v::Float32) = @mem  a.b.c.d.e.x.f.x = v
test404b(a::A4,v::F4)      = @mem  a.b.c.d.e.x.f = v
test405b(a::A4,v::Float32) = @mem  a.b.c.d.e.x.x = v
test406b(a::A4,v::E4)      = @mem  a.b.c.d.e.x = v
test407b(a::A4,v::Float32) = @yolo a.b.c.d.x = v
test408b(a::A4,v::D4)      = @yolo a.b.c.d = v
test409b(a::A4,v::Float32) = @yolo a.b.c.x = v
test410b(a::A4,v::C4)      = @yolo a.b.c = v
test411b(a::A4,v::Float32) = @yolo a.b.x = v
test412b(a::A4,v::B4)      = @mem  a.b = v
test413b(a::A4,v::Float32) = @mem  a.x = v
test414b(a::A4,v::A4)      = @mem  a = v

            v  = 0f0
            v += 1f0; @assert a4.b.c.d.e.x.f.g.x != v;             test401b(a4,v);             @assert a4.b.c.d.e.x.f.g.x == v
            v += 1f0; @assert a4.b.c.d.e.x.f.g   != G4(v);         test402b(a4,G4(v));         @assert a4.b.c.d.e.x.f.g   == G4(v)
            v += 1f0; @assert a4.b.c.d.e.x.f.x   != v;             test403b(a4,v);             @assert a4.b.c.d.e.x.f.x   == v
            v += 1f0; @assert a4.b.c.d.e.x.f     != F4(v,g4);      test404b(a4,F4(v,g4));      @assert a4.b.c.d.e.x.f     == F4(v,g4)
            v += 1f0; @assert a4.b.c.d.e.x.x     != v;             test405b(a4,v);             @assert a4.b.c.d.e.x.x     == v
            v += 1f0; @assert a4.b.c.d.e.x       != E4(v,f4);      test406b(a4,E4(v,f4));      @assert a4.b.c.d.e.x       == E4(v,f4)
try; global v += 1f0; @assert a4.b.c.d.x         != v;             test407b(a4,v);             @assert a4.b.c.d.x         == v            ;  catch e; print(e); end
try; global v += 1f0; @assert a4.b.c.d           != D4(v,Ref(e4)); test408b(a4,D4(v,Ref(e4))); @assert a4.b.c.d.e[]       == e4           # references to equal immutables are not equal
                                                                                @assert a4.b.c.d.x         == v            ;  catch e; print(e); end
try; global v += 1f0; @assert a4.b.c.x           != v;             test409b(a4,v);             @assert a4.b.c.x           == v            ;  catch e; print(e); end
try; global v += 1f0; @assert a4.b.c             != C4(v,d4);      test410b(a4,C4(v,d4));      @assert a4.b.c             == C4(v,d4)     ;  catch e; print(e); end
try; global v += 1f0; @assert a4.b.x             != v;             test411b(a4,v);             @assert a4.b.x             == v            ;  catch e; print(e); end
            v += 1f0; @assert a4.b               != B4(v,c4);      test412b(a4,B4(v,c4));      @assert a4.b               == B4(v,c4)
            v += 1f0; @assert a4.x               != v;             test413b(a4,v);             @assert a4.x               == v
            v += 1f0; @assert a4                 != A4(v,b4);      test414b(a4,A4(v,b4));      # @assert a4               == A4(v,b4) # call by "sharing"

code_native(io,test401b,(A4,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 4, vmov = 1)
code_native(io,test402b,(A4,G4))     ; display_asm_stat_io(io) # (total = 8, movs = 7, mov = 6, vmov = 1)
code_native(io,test403b,(A4,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 4, vmov = 1)
code_native(io,test404b,(A4,F4))     ; display_asm_stat_io(io) # (total = 9, movs = 8, mov = 6, vmov = 2)
code_native(io,test405b,(A4,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 4, vmov = 1)
code_native(io,test406b,(A4,E4))     ; display_asm_stat_io(io) # (total = 15, movs = 13, mov = 13, vmov = 0)
code_native(io,test407b,(A4,Float32)); display_asm_stat_io(io) # (total = 29, movs = 19, mov = 16, vmov = 3)
code_native(io,test408b,(A4,D4))     ; display_asm_stat_io(io) # (total = 19, movs = 16, mov = 15, vmov = 1)
code_native(io,test409b,(A4,Float32)); display_asm_stat_io(io) # (total = 28, movs = 18, mov = 15, vmov = 3)
code_native(io,test410b,(A4,C4))     ; display_asm_stat_io(io) # (total = 18, movs = 15, mov = 14, vmov = 1)
code_native(io,test411b,(A4,Float32)); display_asm_stat_io(io) # (total = 26, movs = 17, mov = 14, vmov = 3)
code_native(io,test412b,(A4,B4))     ; display_asm_stat_io(io) # (total = 24, movs = 8, mov = 8, vmov = 0)
code_native(io,test413b,(A4,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test414b,(A4,A4))     ; display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

################################################################################

# template for the next try, currently the same as for case 4

#         struct G5; x :: Float32;                         end # struct G5 { const float x; };
#         struct F5; x :: Float32; g :: G5;                end # struct F5 { const float x; const G5  g; };
#         struct E5; x :: Float32; f :: F5;                end # struct E5 { const float x; const F5  g; };
# mutable struct D5; x :: Float32; e :: Base.RefValue{E5}; end # struct D5 {       float x;       E5& g; };
#         struct C5; x :: Float32; d :: D5;                end # struct C5 { const float x; const D5& g; };
#         struct B5; x :: Float32; c :: C5;                end # struct B5 { const float x; const C5& g; };
# mutable struct A5; x :: Float32; b :: B5;                end # struct A5 {       float x;       B5& g; };
#
# g5 = G5(0.0f0)
# f5 = F5(0.0f0,g5)
# e5 = E5(0.0f0,f5)
# d5 = D5(0.0f0,Ref(e5))
# c5 = C5(0.0f0,d5)
# b5 = B5(0.0f0,c5)
# a5 = A5(0.0f0,b5)
#
# test501(a::A5,v::Float32) = @mem a.b.c.d.e[].f.g.x = v
# test502(a::A5,v::G5)      = @mem a.b.c.d.e[].f.g = v
# test503(a::A5,v::Float32) = @mem a.b.c.d.e[].f.x = v
# test504(a::A5,v::F5)      = @mem a.b.c.d.e[].f = v
# test505(a::A5,v::Float32) = @mem a.b.c.d.e[].x = v
# test506(a::A5,v::E5)      = @mem a.b.c.d.e[] = v
# test507(a::A5,v::Float32) = @mem a.b.c.d.x = v
# test508(a::A5,v::D5)      = @mem a.b.c.d = v
# test509(a::A5,v::Float32) = @mem a.b.c.x = v
# test510(a::A5,v::C5)      = @mem a.b.c = v
# test511(a::A5,v::Float32) = @mem a.b.x = v
# test512(a::A5,v::B5)      = @mem a.b = v
# test513(a::A5,v::Float32) = @mem a.x = v
# test514(a::A5,v::A5)      = @mem a = v
#
#             v  = 0f0
#             v += 1f0; @assert a5.b.c.d.e[].f.g.x != v;             test501(a5,v);             @assert a5.b.c.d.e[].f.g.x == v
#             v += 1f0; @assert a5.b.c.d.e[].f.g   != G5(v);         test502(a5,G5(v));         @assert a5.b.c.d.e[].f.g   == G5(v)
#             v += 1f0; @assert a5.b.c.d.e[].f.x   != v;             test503(a5,v);             @assert a5.b.c.d.e[].f.x   == v
#             v += 1f0; @assert a5.b.c.d.e[].f     != F5(v,g5);      test504(a5,F5(v,g5));      @assert a5.b.c.d.e[].f     == F5(v,g5)
#             v += 1f0; @assert a5.b.c.d.e[].x     != v;             test505(a5,v);             @assert a5.b.c.d.e[].x     == v
#             v += 1f0; @assert a5.b.c.d.e[]       != E5(v,f5);      test506(a5,E5(v,f5));      @assert a5.b.c.d.e[]       == E5(v,f5)
#             v += 1f0; @assert a5.b.c.d.x         != v;             test507(a5,v);             @assert a5.b.c.d.x         == v
# try; global v += 1f0; @assert a5.b.c.d           != D5(v,Ref(e5)); test508(a5,D5(v,Ref(e5))); @assert a5.b.c.d           == D5(v,Ref(e5)); catch e; print(e); end
# try; global v += 1f0; @assert a5.b.c.x           != v;             test509(a5,v);             @assert a5.b.c.x           == v            ; catch e; print(e); end
# try; global v += 1f0; @assert a5.b.c             != C5(v,d5);      test510(a5,C5(v,d5));      @assert a5.b.c             == C5(v,d5)     ; catch e; print(e); end
# try; global v += 1f0; @assert a5.b.x             != v;             test511(a5,v);             @assert a5.b.x             == v            ; catch e; print(e); end
#             v += 1f0; @assert a5.b               != B5(v,c5);      test512(a5,B5(v,c5));      @assert a5.b               == B5(v,c5)
#             v += 1f0; @assert a5.x               != v;             test513(a5,v);             @assert a5.x               == v
#             v += 1f0; @assert a5                 != A5(v,b5);      test514(a5,A5(v,b5));      # @assert a5               == A5(v,b5) # call by "sharing"
#
# code_native(io,test501,(A5,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 5, vmov = 1)
# code_native(io,test502,(A5,G5))     ; display_asm_stat_io(io) # (total = 8, movs = 7, mov = 6, vmov = 1)
# code_native(io,test503,(A5,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 5, vmov = 1)
# code_native(io,test504,(A5,F5))     ; display_asm_stat_io(io) # (total = 9, movs = 8, mov = 6, vmov = 2)
# code_native(io,test505,(A5,Float32)); display_asm_stat_io(io) # (total = 6, movs = 5, mov = 5, vmov = 1)
# code_native(io,test506,(A5,E5))     ; display_asm_stat_io(io) # (total = 15, movs = 13, mov = 13, vmov = 0)
# # code_native(io,test507,(A5,Float32)); display_asm_stat_io(io)
# # code_native(io,test508,(A5,D5))     ; display_asm_stat_io(io)
# # code_native(io,test509,(A5,Float32)); display_asm_stat_io(io)
# # code_native(io,test510,(A5,C5))     ; display_asm_stat_io(io)
# code_native(io,test511,(A5,Float32)); display_asm_stat_io(io) # (total = 25, movs = 16, mov = 13, vmov = 3)
# code_native(io,test512,(A5,B5))     ; display_asm_stat_io(io) # (total = 25, movs = 8, mov = 8, vmov = 0)
# code_native(io,test513,(A5,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
# code_native(io,test514,(A5,A5))     ; display_asm_stat_io(io) # (total = 3, movs = 2, mov = 2, vmov = 0)

################################################################################

# @macroexpand @mem a.b.c.d.e.f.g.x = v
# a,v = (a1,v)
#
# @macroexpand @mem a.b.c.d = v
# a,v = (a2,D2(v,e2))
#
# @macroexpand @mem a.b.c.d.e.f.g.x = v
# a,v = (a2,v)
#
# @macroexpand @mem a[].b.c.d.e.f.g.x = v
# a,v = (a3r,v)
#
# @macroexpand @mem a[] = v
# a,v = (a3r,A3(v,b3))
#
# @macroexpand @mem a.b = v
# a,v = (a2,B2(v,c2))
#
# @macroexpand @mem a[][][].b.c.d.e.f.g.x = v
# a,v = (a3rrr,v)
#
# @macroexpand @mem a.b.c.d.e[].f.g.x = v
# a,v = (a4,v)
#
# @macroexpand @mem a.b.c.d = v
# a,v = (a2,D2(v,e2))
#
# @macroexpand @mem a.b.c.d = v
# a,v = (a4,D4(v,Ref(e4)))
#
# @macroexpand @mem a[].b.c.d.e.f.g.x = v
# a,v = (a3r,v)

using StaticArrays

m = SArray{Tuple{4,3,2},Int64,3,4*3*2}(1:4*3*2)
r = Ref(m)

@assert r[][1,1,2] == 13
@mem r[][1,1,2] = 77
@assert r[][1,1,2] == 77
@assert r[][1,1,2] == unsafe_load(@ptr r[][1,1,2])

m = Ref(SVector{8,UInt8}(reverse([0xde, 0xad, 0xbe, 0xef, 0xba, 0xad, 0xf0, 0x0d])))
unsafe_load(@ptr m[][8]) # 0xde
unsafe_load(@ptr m[][7]) # 0xad
unsafe_load(@typedptr UInt32 m[][3]) # 0xbeefbaad

@assert pointer_from_objref(r) == @voidptr r[][1,1,1]
@voidptr r[][1,1,1]
@ptr r[][1,1,2]

const B6 = SArray{Tuple{4,3,2},Float32,3,4*3*2}
const A6 = Base.RefValue{B6}
a = Ref(B6(1:4*3*2))

test601(a :: A6, v :: Float32) = @mem a[][1,1,2] = v

v = 77f0
@assert a[][1,1,2] != v; test601(a,v); @assert a[][1,1,2] == v;

code_native(io,test601,(A6,Float32)); display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)

const B6b = SArray{Tuple{4,3,2},Float32,3,4*3*2}
mutable struct A6b
   b :: B6b
end

b = B6b(1:4*3*2)
a = A6b(b)

test601b(a :: A6b, v :: Float32) = @mem a.b[1,1,2] = v
test602b(a :: A6b, x :: Int64, y :: Int64, z :: Int64,  v :: Float32) = @mem a.b[x,y,z] = v

v = 77f0
x, y, z = 2, 1, 2
@assert a.b[1,1,2] != v; test601b(a,v); @assert a.b[1,1,2] == v;
@assert a.b[x,y,z] != v; test602b(a,x,y,z,v); @assert a.b[x,y,z] == v;

code_native(io,test601b,(A6b,Float32));                   display_asm_stat_io(io) # (total = 2, movs = 1, mov = 0, vmov = 1)
code_native(io,test602b,(A6b,Int64,Int64,Int64,Float32)); display_asm_stat_io(io) # (total = 5, movs = 1, mov = 0, vmov = 1)
# leaq    (%rcx,%rcx,2), %rax
# leaq    (%rsi,%rdx,4), %rcx
# leaq    (%rcx,%rax,4), %rax
# vmovss  %xmm0, -68(%rdi,%rax,4)

# TODO: a test of a non-isbitstype SArray; since it allocates an intermediate array and needs double dereferencing
#   same for MArray


text:
function getdim_error(x, condition::Symbol)
    throw(ArgumentError("Specified name ($(repr(condition))) does not match any dimension names for $(repr(x)))"))
end

function getdim_error(x, condition::Function)
    throw(ArgumentError("Method $(Symbol(condition)) is not true for any dimensions of $(repr(x)))"))
end

# I have to abuse @pure here to get type inference to propagate to the axes methods
# So `condition` also has to be pure, which shouldn't be hard because it should basically
# just be comparing symbols
macro defdim(name, condition)

    dim_noerror_name = Symbol(:dim_noerror_, name)

    name_dim = Symbol(name, :dim)
    name_dim_doc = """
        $name_dim(x) -> Int

    Returns the dimension corresponding to the $name.
    """

    nname = Symbol(:n, name)
    nname_doc = """
        $nname(x) -> Int

    Returns the size along the dimension corresponding to the $name.
    """

    has_name_dim = Symbol(:has_, name, :dim)
    has_name_dim_doc = """
        $has_name_dim(x) -> Bool

    Returns `true` if `x` has a dimension corresponding to $name.
    """

    name_axis = Symbol(name, :_axis)
    name_axis_doc = """
        $name_axis(x)

    Returns the axis corresponding to the $name dimension.
    """

    name_indices = Symbol(name, :_indices)
    name_indices_doc = """
        $name_indices(x)

    Returns the indices corresponding to the $name axis
    """

    name_keys = Symbol(name, :_keys)
    name_keys_doc = """
        $name_keys(x)

    Returns the keys corresponding to the $name axis
    """

    name_type = Symbol(name, :_axis_type)
    name_type_doc = """
        $name_type(x)

    Returns the key type corresponding to the $name axis.
    """

    name_format = Symbol(name, :_format)
    name_format_doc = """
        $name_format(x, default)

    Returns the appropriate `DataFormat` for $name.
    """

    to_name_format = Symbol(:to_, name, :_format)
    to_name_format_doc = """
        $to_name_format(dst, src [, default=getdim_error])

    Returns $name data from `src` in a format that is compatible with `dst`.
    `default` refers to a default dimension if data from `dst` or `src` requires
    data in an array format but doesn't specify dimension that corresponds to $name.
    """

    esc(quote
        Base.@pure function $dim_noerror_name(x::Tuple{Vararg{Symbol,N}}) where {N}
            for i in Base.OneTo(N)
                $condition(getfield(x, i)) && return i
            end
            return 0
        end

        @doc $name_dim_doc
        @inline function $name_dim(x)
            d = $dim_noerror_name(dimnames(x))
            if d === 0
                NamedIndicesMeta.getdim_error(x, $condition)
            else
                return d
            end
        end

        @doc $nname_doc
        $nname(x) = Base.size(x, $name_dim(x))

        @doc $has_name_dim_doc
        $has_name_dim(x) = !($dim_noerror_name(dimnames(x)) === 0)

        @doc $name_axis_doc
        $name_axis(x) = axes(x, $name_dim(x))

        @doc $name_keys_doc
        $name_keys(x) = keys($name_axis(x))

        @doc $name_indices_doc
        $name_indices(x) = values($name_axis(x))

        @doc $name_type_doc
        $name_type(x) = keytype($name_axis(x))

        @doc $name_format_doc
        function $name_format(x::AbstractArray)
            return NamedIndicesMeta.get_format(x, $condition)
        end

        @doc $to_name_format
        function $to_name_format(dst, src)
            return NamedIndicesMeta.to_format(dst, src, $condition)
        end

        nothing
    end)
end


text:
# Julia wrapper for header: /Users/jgoldfar/Documents/misc/julia/Silo/gen/../deps/usr/include/pmpio.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

function PMPIO_Init(numFiles::Cint,ioMode::PMPIO_iomode_t,mpiComm::Cint,mpiTag::Cint,createCb::PMPIO_CreateFileCallBack,openCb::PMPIO_OpenFileCallBack,closeCb::PMPIO_CloseFileCallBack,userData)
    ccall((:PMPIO_Init,libsilo),Ptr{PMPIO_baton_t},(Cint,PMPIO_iomode_t,Cint,Cint,PMPIO_CreateFileCallBack,PMPIO_OpenFileCallBack,PMPIO_CloseFileCallBack,Ptr{Nothing}),numFiles,ioMode,mpiComm,mpiTag,createCb,openCb,closeCb,userData)
end

function PMPIO_Finish(bat)
    ccall((:PMPIO_Finish,libsilo),Nothing,(Ptr{PMPIO_baton_t},),bat)
end

function PMPIO_WaitForBaton(Bat,fname,nsname)
    ccall((:PMPIO_WaitForBaton,libsilo),Ptr{Nothing},(Ptr{PMPIO_baton_t},Cstring,Cstring),Bat,fname,nsname)
end

function PMPIO_HandOffBaton(Bat,file)
    ccall((:PMPIO_HandOffBaton,libsilo),Nothing,(Ptr{PMPIO_baton_t},Ptr{Nothing}),Bat,file)
end

function PMPIO_GroupRank(Bat,rankInComm::Cint)
    ccall((:PMPIO_GroupRank,libsilo),Cint,(Ptr{PMPIO_baton_t},Cint),Bat,rankInComm)
end

function PMPIO_RankInGroup(Bat,rankInComm::Cint)
    ccall((:PMPIO_RankInGroup,libsilo),Cint,(Ptr{PMPIO_baton_t},Cint),Bat,rankInComm)
end


text:
using Plots, Statistics, MLDatasets, LinearAlgebra, Clustering, Printf, Dates, Random, XLSX
path1=realpath(dirname(@__FILE__)*"/../..")
include(string(path1,"/src/ANNBN.jl"))

path_project="~\\finance\\"
file1="GE.xlsx"
sheet1="GE"
range1="A1:G1260"
all_vars = XLSX.readdata(string(path_project,file1),sheet1,range1)

yy_train=float.(all_vars[2:end,5])
plot(yy_train)
xx_train=[yy_train[5:end-1] yy_train[4:end-2] yy_train[3:end-3] yy_train[2:end-4] yy_train[1:end-5]]

yy_train=yy_train[6:end]
aa=xx_train\yy_train
pred=xx_train*aa

plot(yy_train[1:100])
plot!(pred[1:100])

mean(abs.(yy_train-pred))
mean(abs.(yy_train[1:end-1]-yy_train[2:end]))

i_train=size(xx_train,1)
vars=size(xx_train,2)
neurons=Int64(floor(i_train/(vars+1)))
# neurons=50
inds_all,n_per_part=ann_by_parts.___clustering(neurons,xx_train,200)
cc1=100.0
a_all,a_layer1,layer1=ann_by_parts.train_layer_1_rbf(neurons,vars,i_train,n_per_part,
            inds_all,xx_train,yy_train,cc1)
predl1=ann_by_parts.predict_new_rbf(a_all,a_layer1,xx_train,i_train,neurons,
            n_per_part,inds_all,xx_train,vars,cc1)
# predl1_test=ann_by_parts.predict_new_rbf(a_all,a_layer1,xx_test,i_test,neurons,
            # n_per_part,inds_all,xx_train,vars,cc1)
maetr=mean(abs.(yy_train-predl1))
# scatter(yy_train,predl1,label=@sprintf("fann %1.4f",maetr))
# maete=mean(abs.(yy_test-predl1_test))

plot!(predl1[1:100])

nof_folds=20 # you may change these
i_fold=Int64(floor(0.8i_train))
include("ann_by_parts.jl")
maetrs,a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_fold_all,neurons_all=ann_by_parts.fit_nfolds!(xx_train,yy_train,nof_folds,vars,i_train,i_fold,cc1);
predl1_all=ann_by_parts.predict_nfolds(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_train,i_train,vars,cc1,xx_fold_all,nof_folds,maetrs,neurons_all)
predl1_train=sum(predl1_all)./sum(1 ./maetrs)
maetr=mean(abs.(yy_train-predl1_train))
# predl1_all=ann_by_parts.predict_nfolds(a_all_all,a_layer1_all,n_per_part_all,inds_all_all,xx_test,i_test,vars,cc1,xx_fold_all,nof_folds,maetrs,neurons_all)
# predl1_test=sum(predl1_all)./sum(1 ./maetrs)
# maete=mean(abs.(yy_test-predl1_test))

plot!(predl1_train[1:100])


text:
import Pkg

Pkg.update()

root_directory = dirname(@__DIR__)

amdgpu = Pkg.PackageSpec(path = root_directory)
Pkg.develop(amdgpu)
Pkg.build()
Pkg.precompile()

import AMDGPU


text:
"""
functions to controll stroke width , setting which texture is currently active and actions undoing
"""
module OtherKeyboardActions
using ModernGL, ..DisplayWords, ..StructsManag, Setfield, ..PrepareWindow,   ..DataStructs , Rocket, GLFW,Dictionaries,  ..ForDisplayStructs, ..TextureManag,  ..OpenGLDisplayUtils,  ..Uniforms, Match, Parameters,DataTypesBasic

"""
in case we want to  get new number set for manual modifications
    toBeSavedForBack - just marks weather we wat to save the info how to undo latest action
    - false if we invoke it from undoing
"""
function processKeysInfo(numbb::Identity{Int64}
                        ,actor::SyncActor{Any, ActorWithOpenGlObjects}
                        ,keyInfo::KeyboardStruct
                        ,toBeSavedForBack::Bool = true) where T

    valueForMasToSett = valueForMasToSetStruct(value = numbb.value)
    old = actor.actor.valueForMasToSet.value
    actor.actor.valueForMasToSet =valueForMasToSett

    textureList= actor.actor.textureToModifyVec
    # in case we increase number it should not be outside of the possible values
    if(!isempty(textureList))
        textureList[1].minAndMaxValue[2]= max(textureList[1].minAndMaxValue[2],numbb )
    end#if

    updateImagesDisplayed(actor.actor.currentlyDispDat
    , actor.actor.mainForDisplayObjects
    , actor.actor.textDispObj
    , actor.actor.calcDimsStruct ,valueForMasToSett )
# for undoing action
if(toBeSavedForBack)
    addToforUndoVector(actor, ()-> processKeysInfo( Option(old),actor, keyInfo,false ))
end

end#processKeysInfo

"""
In order to enable undoing last action we just invoke last function from list
"""
function processKeysInfoUndo(numbb::Identity{Bool},actor::SyncActor{Any, ActorWithOpenGlObjects},keyInfo::KeyboardStruct ) where T
    if(!isempty(actor.actor.forUndoVector))
     pop!(actor.actor.forUndoVector)()
    end
end#processKeysInfo

"""
when tab plus will be pressed it will increase stroke width
when tab minus will be pressed it will increase stroke width
"""
function processKeysInfo(annot::Identity{AnnotationStruct}
    ,actor::SyncActor{Any, ActorWithOpenGlObjects}
    ,keyInfo::KeyboardStruct
    ,toBeSavedForBack::Bool = true) where T

    textureList = actor.actor.textureToModifyVec
    if (!isempty(textureList))
        texture= textureList[1]
        oldsWidth = texture.strokeWidth
        texture.strokeWidth= oldsWidth+=annot.value.strokeWidthChange
        # for undoing action
        if(toBeSavedForBack)
        addToforUndoVector(actor, ()-> processKeysInfo( Option(AnnotationStruct(oldsWidth)),actor, keyInfo,false ))
        end#if
    end#if

end#processKeysInfo

end#OtherKeyboardActions