text:
import QXGraphDecompositions as qxg
import LightGraphs
import QXTns

export convert_to_graph, convert_to_line_graph
export contraction_scheme
export flow_cutter_contraction_plan, min_fill_contraction_plan
export netcon

# **************************************************************************************** #
#                             Converting tensor networks to graphs
# **************************************************************************************** #

"""
    convert_to_graph(tn::TensorNetwork)

Return a SimpleGraph representing the graph structure of the tensor network `tn`.
"""
function convert_to_graph(tn::TensorNetwork)
    g = LightGraphs.SimpleGraph(length(tn))
    tensor_map = Dict{Symbol, Int64}([k => i for (i, k) in enumerate(keys(tn))])
    # TODO: add structures to tensor network to avoid N^2 complexity here
    for bond in bonds(tn)
        tensors = tn[bond]
        if length(tensors) == 2
            LightGraphs.add_edge!(g, tensor_map[tensors[1]], tensor_map[tensors[2]])
        end
    end
    g
end

convert_to_graph(tnc::TensorNetworkCircuit) = convert_to_graph(tnc.tn)

"""
    convert_to_line_graph(tn::TensorNetwork; use_hyperedges::Bool=false)

Create a labeled graph representing the line graph of 'tn'.

If `use_hyperedges` is false, each edge in `tn` is assigned a symbol which labels a unique
vertex in the line graph. If it is true, the hyperedges of `tn` are assigned symbols instead.

A dictionary is also returned which maps these symbols to the corresponding edges which are
represented as tensor sets.

# Keywords
- `use_hyperedges::Bool=false`: set if hyperedges should be used to create the line graph.
"""
function convert_to_line_graph(tn::TensorNetwork; use_hyperedges::Bool=false)
    edges = use_hyperedges ? get_hyperedges(tn) : [tn[ind] for ind in bonds(tn)]

    # Create the line graph and dict mapping labels to hyperedges.
    g = qxg.LabeledGraph(length(edges))
    symbol_map = Dict{Symbol, Array{Symbol, 1}}(qxg.labels(g) .=> edges)

    # Connect the vertices of g if corresponding hyperedges overlap.
    for i in 1:length(edges)-1
        for j in i+1:length(edges)
            if !isempty(intersect(edges[i], edges[j]))
                qxg.add_edge!(g, i, j)
            end
        end
    end
    g, symbol_map
end

convert_to_line_graph(tnc::TensorNetworkCircuit; kwargs...) = convert_to_line_graph(tnc.tn; kwargs...)

# **************************************************************************************** #
#                              Contraction Planning
# **************************************************************************************** #

"""
    flow_cutter_contraction_plan(tn::TensorNetwork;
                                 time::Integer=60,
                                 seed::Integer=-1,
                                 hypergraph::Bool=false)

Use flow cutter to create a contraction plan for 'tn'.

If flow cutter does not find a tree decomposotion to turn into a contraction plan within
the specified amount of time then the min fill heuristic is used to find a contraction plan.

# Keywords
- `time::Integer=60`: The number of seconds to run the flow cutter for.
- `seed::Integer=-1`: Sets the seed used by flow cutter. If negative then flow cutter will
                      choose a seed.
- `hypergraph::Bool=false`: Sets if hyperedges in `tn` should be accounted for.
"""
function flow_cutter_contraction_plan(tn::TensorNetwork;
                                      time::Integer=60,
                                      seed::Integer=-1,
                                      hypergraph::Bool=false)
    # Convert tn to a line graph and pass it to flow cutter to find an tree decomposition.
    lg, symbol_map = convert_to_line_graph(tn, use_hyperedges=hypergraph)

    # Use flow cutter to try find a tree decomposition of the line graph.
    td = qxg.flow_cutter(lg, time; seed=seed)

    # If a tree decomposition was found, convert it into a vertex elimination order for lg,
    # otherwise use the min fill heuristic to find an elimination order.
    if haskey(td, :treewidth)
        order = Symbol.(qxg.order_from_tree_decomposition(td))
    else
        tw, order = qxg.min_fill(lg)
    end

    # Convert the elimination order to an array of Index structs in tn.
    order = [symbol_map[lg_vertex] for lg_vertex in order]

    # Convert the elimination order into a contraction plan.
    order_to_contraction_plan(order, tn)
end

flow_cutter_contraction_plan(tnc::TensorNetworkCircuit; kwargs...) = flow_cutter_contraction_plan(tnc.tn; kwargs...)

"""
    min_fill_contraction_plan(tn::TensorNetwork;
                              hypergraph::Bool=false)

Use the min fill heuristic to create a contraction plan for 'tn'.

# Keywords
- `hypergraph::Bool=false`: set if hyperedges in `tn` should be accounted for.
"""
function min_fill_contraction_plan(tn::TensorNetwork;
                                   hypergraph::Bool=false)
    # Convert tn to a line graph and pass it to flow cutter to find an tree decomposition.
    lg, symbol_map = convert_to_line_graph(tn, use_hyperedges=hypergraph)

    # Use flow cutter to try find a tree decomposition of the line graph.
    tw, order = qxg.min_fill(lg)

    # Convert the elimination order to an array of Index structs in tn.
    order = [symbol_map[lg_vertex] for lg_vertex in order]

    # Convert the elimination order into a contraction plan.
    order_to_contraction_plan(order, tn)
end

min_fill_contraction_plan(tnc::TensorNetworkCircuit; kwargs...) = min_fill_contraction_plan(tnc.tn; kwargs...)

"""
    min_fill_contraction_plan(tensors::OrderedDict{Symbol, Array{Index, 1}})

Use the min fill heuristic to create a contraction plan for the set of tensors described by
the OrderedDict `tensors`.

The keys of `tensors` are assumed to be ids/names of tensors and the values are arrays of
indices belonging to the corrseponding tensor.
"""
function min_fill_contraction_plan(tensors::OrderedDict{Symbol, Array{Index, 1}})
    # Create a graph for the tensors and convert it to to a line graph.
    tensor_ids = collect(keys(tensors))
    g = qxg.LabeledGraph(length(tensor_ids))
    for i = 1:length(tensor_ids)-1
        for j = i+1:length(tensor_ids)
            if !isempty(intersect(tensors[tensor_ids[i]], tensors[tensor_ids[j]]))
                qxg.add_edge!(g, i, j)
            end
        end
    end
    lg = qxg.line_graph(g)

    # Use the min fill heuristic to find a vertex elimination order for the line graph.
    tw, order = qxg.min_fill(lg)

    # Convert the elimination order to an array of tensor symbol pairs.
    order = [parse.(Int, split(String(lg_vertex), '_')) for lg_vertex in order]
    order = [[tensor_ids[edge[1]], tensor_ids[edge[2]]] for edge in order]

    # Convert the elimination order into a contraction plan.
    order_to_contraction_plan(order, tensors)
end

# **************************************************************************************** #
#                         Contraction with automatic slicing
# **************************************************************************************** #

"""
    expand_slice_bonds_to_hyperindices(tn::TensorNetwork, bonds_to_slice::Vector{<: Index})

Given a list of bonds to slice it is necessary to expand these to include bonds in the same hyper edge group. In
this function for each of the bonds to slice we identify their hyper indices if any and return an array of
groups of hyper edges to slice.
"""
function expand_slice_bonds_to_hyperindices(tn::QXTns.TensorNetwork, bonds_to_slice::Vector{<: Index})
    bond_groups = Array{Array{<:Index, 1}, 1}()
    if length(bonds_to_slice) > 0
        push!(bond_groups, find_connected_indices(tn, bonds_to_slice[1]))
        if length(bonds_to_slice) > 1
            for b in bonds_to_slice[2:end]
                if !any([b in g for g in bond_groups])
                    push!(bond_groups, find_connected_indices(tn, b))
                end
            end
        end
    end
    bond_groups
end

"""
    contraction_scheme(tn::TensorNetwork, num::Integer)

Return an array of 'num' indices to slice in the tensor network 'tn' and a contraction plan
for the remaining tensor network.

# Keywords
- `time::Integer=120`: number of seconds to run quickbb when looking for contraction plan.
- `qbb_order::Symbol=:min_fill`: the branching order to be used by quickbb (:random or :min_fill).
- `lb::Bool=false`: set if a lowerbound for the treewidth should be computed.
- `score_function::Symbol=:direct_treewidth`: function to maximise when selecting vertices
                                            to remove. (:degree, :direct_treewidth)
- `hypergraph::Bool=true`: set if hyperedges exist in `tn` and should be accounted for.
"""
function contraction_scheme(tn::TensorNetwork, num::Integer;
                            time::Integer=120,
                            seed::Integer=-1,
                            score_function::Symbol=:direct_treewidth,
                            hypergraph::Bool=true)

    # Create the line graph for the given tn.
    lg, symbol_map = convert_to_line_graph(tn; use_hyperedges=hypergraph)

    # Use flow cutter to try find a tree decomposition of the line graph.
    td = qxg.flow_cutter(lg, time; seed=seed)
    flow_cutter_metadata = OrderedDict(1:length(td[:comments]) .=> td[:comments])

    # If a tree decomposition was found, convert it into a vertex elimination order for lg,
    # otherwise use the min fill heuristic to find an elimination order.
    if haskey(td, :treewidth)
        order = Symbol.(qxg.order_from_tree_decomposition(td))
        tw = td[:treewidth]
        method_used = "flow cutter"
    else
        tw, order = qxg.min_fill(lg)
        method_used = "min fill heuristic"
    end

    # Create a dictionary for metadata regarding the contraction plan.
    contraction_metadata = OrderedDict{String, Any}()
    contraction_metadata["Method used"] = method_used
    contraction_metadata["Treewidth"] = tw
    contraction_metadata["Time allocated"] = time
    contraction_metadata["Seed used"] = seed
    contraction_metadata["Returned metadata"] = flow_cutter_metadata
    contraction_metadata["Hypergraph used"] = hypergraph
    contraction_metadata["Hyperedge contraction method"] = "Netcon where possible, min fill heuristic otherwise."

    # Use the greedy treewidth deletion algorithm to select indices in tn to slice.
    scheme = qxg.greedy_treewidth_deletion(lg, num;
                                           score_function=score_function,
                                           elim_order=order)
    sliced_lg, edges_to_slice, modified_orders, treewidths = scheme

    # Create a dictionary for metadata regarding sliced edges.
    slicing_metadata = OrderedDict{String, Any}()
    slicing_metadata["Method used"] = "greedy treewidth deletion"
    slicing_metadata["Edges sliced"] = num
    slicing_metadata["Score fucntion used"] = score_function
    slicing_metadata["Treewidths after slicing consecutive edges"] = treewidths

    # Convert the contraction plan to an array of Index structs before returning.
    contraction_order = !isempty(modified_orders) ? modified_orders[end] : order
    contraction_order = [symbol_map[edge_symbol] for edge_symbol in contraction_order]
    edges_to_slice = [symbol_map[edge_symbol] for edge_symbol in edges_to_slice]
    contraction_plan = order_to_contraction_plan(contraction_order, tn)

    # Convert edges_to_slice to an array of Index structs to slice.
    # TODO: I'm not confident there are no edge cases were this method for converting
    # edges to index arrays doesn't work. For example, when there are two consecutive
    # 2-qubit gates sharing both a hyperedge and a normal edge. Would be nice if it could
    # be made a bit cleaner also.
    indices_to_slice = Array{Array{<:Index, 1}, 1}(undef, length(edges_to_slice))
    for (i, edge) in enumerate(edges_to_slice)
        indices = Array{Index, 1}()
        for i in 1:length(edge)-1
            for j in i+1:length(edge)
                union!(indices, intersect(tn[edge[i]].indices, tn[edge[j]].indices))
            end
        end
        indices_to_slice[i] = indices
    end
    if length(indices_to_slice) > 0
        indices_to_slice = vcat(indices_to_slice...)
    else # in case there are no edges to slice ensure return array has correct type
        indices_to_slice = Index[]
    end
    index_groups_to_slice = expand_slice_bonds_to_hyperindices(tn, indices_to_slice)

    metadata = OrderedDict{String, Any}()
    metadata["Determination of contraction plan"] = contraction_metadata
    metadata["Slicing"] = slicing_metadata

    index_groups_to_slice, contraction_plan, metadata
end

contraction_scheme(tnc::TensorNetworkCircuit, args...; kwargs...) = contraction_scheme(tnc.tn, args...; kwargs...)

# **************************************************************************************** #
#                                   Netcon functions
# **************************************************************************************** #
import TensorOperations.optimaltree

"""
    netcon(tn::TensorNetwork}

Call the TensorOperations netcon implementation on the tensors in `tensor_map`.

`tensor_map` is assumed to be a dictionary mapping tensor ids to arrays of indices connected
to the corresponding tensor.
"""
function netcon(tn::TensorNetwork, tensor_map::Union{OrderedDict{Symbol, <:Array{<:Index, 1}}, Nothing}=nothing)
    if tensor_map === nothing
        tensor_map = OrderedDict(keys(tn) .=> QXTns.inds.(values(tn)))
    end
    # Generate the arguments for the netcon method from the given network.
    tensor_indices, index_dims, tensor_labels = create_netcon_input(tensor_map)

    # Run the netcon method to find the optimal contraction plan.
    contraction_tree, cost = optimaltree(tensor_indices, index_dims)

    # Replace integers in contraction tree with node labels before returning
    plan, final_tensor_id = tree_to_contraction_plan(tn, contraction_tree, tensor_labels)
    plan
end

netcon(tnc::TensorNetworkCircuit) = netcon(tnc.tn)

"""
    create_netcon_input(tensor_map::OrderedDict{Symbol, <:Array{<:Index, 1}})

Create arguments for netcon to find the optimal contraction plan for the tensors
contatined in `tensor_map`.
"""
function create_netcon_input(tensor_map::OrderedDict{Symbol, <:Array{<:Index, 1}})
    # Collect the indices and corresponding dimensions of the network into arrays.
    num_tensors = length(tensor_map)

    # Allocate space for indices, dims and labels of tensors to contract.
    tensor_indices = Array{Array{Symbol, 1}, 1}(undef, num_tensors)
    tensor_dims = Array{Array{Int, 1}, 1}(undef, num_tensors)
    tensor_labels = Array{Symbol, 1}(undef, num_tensors)

    index_ids = Dict{Index, Symbol}()
    for (i, (tensor_id, tensor_inds)) in enumerate(tensor_map)
        tensor_indices[i] = [get!(index_ids, ind, Symbol("ind_$(i)_$(j)"))
                            for (j, ind) in enumerate(tensor_inds)]
        tensor_dims[i] = [QXTns.dim(ind) for ind in tensor_inds]
        tensor_labels[i] = tensor_id
    end

    # Make a dictionary of index dimensions for netcon.
    labels = reduce(vcat, tensor_indices); dims = reduce(vcat, tensor_dims)
    index_dims = Dict{Symbol, Int}(labels .=> dims)

    # The netcon implementation only needs tensor_indices and index_dims as
    # input but it returns a nested array of integers which need to be replaced
    # by node labels for contraction. To this end, tensor_labels is also
    # returned.
    tensor_indices, index_dims, tensor_labels
end

# **************************************************************************************** #
#       Converting edge elimination orders & contraction trees into contraction plans
# **************************************************************************************** #

"""
    order_to_contraction_plan(elimination_order::Array{<:Array{Symbol, 1}, 1},
                              tn::Union{TensorNetwork, OrderedDict{Symbol, Array{Index, 1}}}
                              )::Array{NTuple{3, Symbol}, 1}

Convert the given edge elimination order into a contraction plan for `tn`.

`tn` can be a TensorNetwork or an OrderedDict describing a set of tensors, in which case
the keys of `tn` are assumed to be tensor ids/names and the values are arrays containing the
indices of the corresponding tensor.
"""
function order_to_contraction_plan(elimination_order::Array{Array{Symbol, 1}, 1},
                                   tn::Union{TensorNetwork, OrderedDict{Symbol, Array{Index, 1}}}
                                   )::Array{NTuple{3, Symbol}, 1}
    # An array to hold the contraction plan.
    plan = Array{NTuple{3, Symbol}, 1}()

    # A dictionary to keep track of which tensors are replaced by which intermediate tensors
    # at different stages of the contraction process. Initially, before any pairwise
    # contractions, none of the tensors are replaced by intermediates, so all tensor ids are
    # mapped to themselves.
    intermediate_tensor_map = Dict{Symbol, Symbol}(keys(tn) .=> keys(tn))

    # Convert each edge in the elimination order to a set of pairwise contractions and
    # append it to plan.
    for (i, edge) in enumerate(elimination_order)
        if length(edge) == 2
            # For edges consisting of just two tensors, find the intermediates they belong
            # to and add their pairwise contraction to the plan.
            A_id = _get_intermediate_tensor(intermediate_tensor_map, edge[1])
            B_id = _get_intermediate_tensor(intermediate_tensor_map, edge[2])
            if !(A_id == B_id)
                # id for the new intermediate created by contracting A_id and B_id.
                I_id = Symbol("I$i")
                append!(plan, [(A_id, B_id, I_id)])

                # Update intermediate_tensor_map with new intermediate.
                intermediate_tensor_map[A_id] = I_id
                intermediate_tensor_map[B_id] = I_id
                intermediate_tensor_map[I_id] = I_id
            end

        elseif length(edge) > 2
            # For edges with more than 2 tensors, collect all of the tensors and
            # intermediate tensors that belong to the edge and find a contraction plan for
            # them. Append the contraction this plan to plan.
            tensors_to_contract = OrderedDict{Symbol, Array{Index, 1}}()
            for t_id in edge
                I_id = _get_intermediate_tensor(intermediate_tensor_map, t_id)
                inds = typeof(tn) <: TensorNetwork ? QXTns.inds(tn[t_id]) : tn[t_id]
                tensors_to_contract[I_id] = symdiff(get(tensors_to_contract, I_id, []), inds)
            end
            if length(tensors_to_contract) > 1
                # Check if netcon can be used on the given set of tensors. If not, use
                # a fallback method to find a contraction plan.
                tensor_sizes = prod.([QXTns.dim.(inds) for inds in values(tensors_to_contract)])
                if length(tensors_to_contract) < 37 && any(tensor_sizes .< 2^62)
                    local_contraction_plan = netcon(tn, tensors_to_contract)
                else
                    local_contraction_plan = min_fill_contraction_plan(tensors_to_contract)
                end
                append!(plan, local_contraction_plan)

                # Update intermediate_tensor_map with new intermediates.
                for (A_id, B_id, I_id) in local_contraction_plan
                    intermediate_tensor_map[A_id] = I_id
                    intermediate_tensor_map[B_id] = I_id
                    intermediate_tensor_map[I_id] = I_id
                end
            end
        end
    end
    plan
end

"""
    _get_intermediate_tensor(intermediate_tensor_map::Dict{Symbol, Symbol}, t_id::Symbol)::Symbol

Return the id of the intermediate tensor of which the tensor `t_id` is a factor.

During a tensor network contraction, tensors and intermediate tensors are contracted
together and replaced by the result. The dictionary `intermediate_tensor_map` is
assumed to describe which intermediate tensor has replaced a given tensor during a
contraction of a tensor network.
"""
function _get_intermediate_tensor(intermediate_tensor_map::Dict{Symbol, Symbol}, t_id::Symbol)::Symbol
    while true
        I_id = intermediate_tensor_map[t_id]
        # If t_id = I_id then t_id has not been replaced by an intermediate tensor yet.
        (t_id == I_id) && return I_id
        t_id = I_id
    end
end

"""
    tree_to_contraction_plan(tn::TensorNetwork, tree::Array{Any, 1}, labels::Array{Symbol, 1})::Tuple{Array{NTuple{3, Symbol}, 1}, Symbol}

Return the contraction plan described by the netcon contraction tree stored in 'tree'. The
tensors labeled by integers in `tree` are labeled by the corresponding symbol in `labels` in
the returned contraction plan.
"""
function tree_to_contraction_plan(tn::TensorNetwork, tree::Array{Any, 1}, labels::Array{Symbol, 1})::Tuple{Array{NTuple{3, Symbol}, 1}, Symbol}
    plan_A, A_id = tree_to_contraction_plan(tn, tree[1], labels)
    plan_B, B_id = tree_to_contraction_plan(tn, tree[2], labels)
    C_id = next_tensor_id!(tn)
    return [plan_A; plan_B; NTuple{3, Symbol}[(A_id, B_id, C_id)]], C_id
end

function tree_to_contraction_plan(::TensorNetwork, tree::Int, labels::Array{Symbol, 1})::Tuple{Array{NTuple{3, Symbol}, 1}, Symbol}
    [], labels[tree]
end

# **************************************************************************************** #
#       Simple pre contraction planning operations
# **************************************************************************************** #

"""
    contract_low_rank_tenors!(tn::TensorNetwork, max_rank=1)

Contract in tensors which have rank that's lower than or equal to
provided max rank
"""
function contract_low_rank_tenors!(tn::TensorNetwork, max_rank=1)
    # function to find first tensor with rank less than or equal to max rank
    # and with tensors connected to all its indices
    find_t = () -> findfirst(tn) do x
        if length(size(x)) <= max_rank
            ns = unique(vcat([tn[x] for x in QXTns.inds(x)]...))
            if length(ns) == length(size(x)) + 1 return true end
        end
        false
    end
    t = find_t()
    plan = Vector{NTuple{3, Symbol}}()
    while t !== nothing
        ns = neighbours(tn, t)
        # sort in increasing order of rank
        ns = sort(ns, by=x -> length(size(tn[x])))
        o = contract_pair!(tn, t, ns[1], mock=true)
        push!(plan, (t, ns[1], o))
        t = find_t()
    end
    plan
end


text:
export mut_frequency_vec, time_sample, ensemble_evolution

function bitaccumulate!(vec::Vector,b::Int64)
    for i in 1:bit_length(b)
        if readbit(b, i) == 1
            vec[i] += 1
        end
    end
end

function mut_frequency_vec(vp)
    n_gene = length(vp.equilibrium_sampler)
    vec = zeros(Int64,n_gene);
    for virus in vp.pop
        bitaccumulate!(vec,virus.genotype)
    end
    return vec ./ length(vp.pop)
end

function time_sample(vp, time_points)
	reduce(hcat,
        begin
        run_viruses_stop(vp, time)
        mut_frequency_vec(vp)
        end
        for time in time_points)
end

function ensemble_evolution(θ, ab_profile, time_points; n_ensemble = 100, antibody = false, kwds...)
    reduce((x,y)->cat(x,y;dims = 3),
        time_sample(initialize_viral_population(θ, ab_profile; antibody, kwds...), time_points)
        for _ in 1:n_ensemble)
end

ext: jl

max_stars_repo_name: StatPhysBio/EscapeSimulator

max_issues_repo_name: StatPhysBio/EscapeSimulator

max_forks_repo_name: StatPhysBio/EscapeSimulator


text:
using Bridge, StaticArrays, Distributions
using Test, Statistics, Random, LinearAlgebra
using Bridge.Models

include("partialparam.jl")
struct IntegratedDiffusion <: ContinuousTimeProcess{ℝ{2}}
    γ::Float64
end

βu(t, x::Float64, P::IntegratedDiffusion) = - (x+sin(x)) + 1/2
Bridge.b(t::Float64, x, P::IntegratedDiffusion) = ℝ{2}(x[2], βu(t, x[2], P))
Bridge.σ(t, x, P::IntegratedDiffusion) = ℝ{2}(0.0, P.γ)

Bridge.constdiff(::IntegratedDiffusion) = true

struct IntegratedDiffusionAux <: ContinuousTimeProcess{ℝ{2}}
    γ::Float64
end

βu(t, x::Float64, P::IntegratedDiffusionAux) = -x + 1/2
Bridge.b(t::Float64, x, P::IntegratedDiffusionAux) = ℝ{2}(x[2], βu(t, x[2], P))
Bridge.σ(t, P::IntegratedDiffusionAux) =  ℝ{2}(0.0, P.γ)
Bridge.σ(t, x, P::IntegratedDiffusionAux) = Bridge.σ(t, P)

Bridge.B(t, P::IntegratedDiffusionAux) = @SMatrix [0.0 1.0; 0.0 -1.0]
Bridge.β(t, P::IntegratedDiffusionAux) = ℝ{2}(0, 1/2)
Bridge.a(t, P::IntegratedDiffusionAux) = @SMatrix [0.0 0.0; 0.0 P.γ^2]

Bridge.constdiff(::IntegratedDiffusionAux) = true

# Generate Data
Random.seed!(1)

P = IntegratedDiffusion(γ)
Pt = IntegratedDiffusionAux(γ)

W = sample(tt, Wiener())

X = solve(Euler(), x0, W, P)

# Solve Backward Recursion

S2 = typeof(L)
S = typeof(L*L')
T = typeof(diag(L*L'))

N = length(tt)
Lt = zeros(S2, N)
Mt = zeros(S, N)
μt = zeros(T, N)

Bridge.partialbridgeode!(Bridge.R3(), tt, L, Σ, Lt, Mt, μt, Pt)

j = 10

@test norm((μt[j+1] - μt[j])/dt - (-Lt[j+1]*Bridge.β(tt[j+1], Pt))) < 0.01
@test norm((inv(Mt[j+1]) - inv(Mt[j]))/dt - (-Lt[j+1]*Bridge.a(tt[j+1], Pt)*Lt[j+1]')) < 0.01

Po = Bridge.PartialBridge(tt, P, Pt, L, v, Σ)

@test Po.L == Lt

W = sample(tt, Wiener())

Xo = copy(X)
solve!(Euler(), Xo, x0, W, Po)

@show LP = log(pdf(Normal(μt[1][1] + (Lt[1]*x0)[1], (Mt[1][1])^(-0.5)), v[1]))

# Likelihood

ll = llikelihood(Bridge.LeftRule(), Xo, Po)

@testset "MCMC" begin

    # MCMC parameter

    iterations = 10000
    subsamples = 0:100:iterations
    ρ = 0.9

    # initalization
    sample!(W, Wiener())
    solve!(Euler(), X, x0, W, Po)
    ll = llikelihood(Bridge.LeftRule(), X, Po)

    acc = 0

    Wo = copy(W)
    W2 = copy(W)

    XX = Any[]
    if 0 in subsamples
        push!(XX, copy(X))
    end
    for iter in 1:iterations
        # Proposal
        sample!(W2, Wiener())
        Wo.yy .= ρ*W.yy + sqrt(1-ρ^2)*W2.yy

        solve!(Euler(), Xo, x0, Wo, Po)
        llo = llikelihood(Bridge.LeftRule(), Xo, Po)
        if log(rand()) <= llo - ll
            X.yy .= Xo.yy
            W.yy .= Wo.yy
            ll = llo
            acc += 1
        end
        if iter in subsamples
            push!(XX, copy(X))
        end
    end
    @test 1 < acc < iterations

end


text:
# Much of this code is lifted from LineSearches.jl
# I modified it to accept StaticArrays and not allocate

# Some of the optimization code is adapted from Optim.jl

@with_kw struct BackTracking{TF, TI}
    c_1::TF = 1e-4
    ρ_hi::TF = 0.5
    ρ_lo::TF = 0.1
    iterations::TI = 1_000
    maxstep::TF = Inf
end

abstract type BackTrackingOrder end
struct Order2 <: BackTrackingOrder end
struct Order3 <: BackTrackingOrder end
struct Order0 <: BackTrackingOrder end
ordernum(::Order2) = 2
ordernum(::Order3) = 3
struct StaticOptimizationResults{Tx, Th, Tf}
    initial_x::Tx
    minimizer::Tx
    minimum::Tf
    iterations::Int
    g_converged::Bool
    g_tol::Tf
    f_calls::Int
    g_calls::Int
    g::Tx
    h::Th
end
struct StaticBFGS end
struct StaticNewton end

setresult(x::AbstractVector) = DiffResults.HessianResult(x)
setresult(x::Number) = DiffResults.DiffResult(x, x, x)
initialh(x::StaticVector{P,T}) where {P,T} = SMatrix{P,P,T}(I)
initialh(x::AbstractVector) = [i == j ? 1. : 0. for i in 1:size(x, 1), j in 1:size(x, 1)]
initialh(x::Number) = one(x)

setgradient!(res, f, x::AbstractVector) = ForwardDiff.gradient!(res, f, x)
setgradient!(res, f, x::Number) = ForwardDiff.derivative!(res, f, x)
sethessian!(res, f, x::AbstractVector) = ForwardDiff.hessian!(res, f, x)
function sethessian!(res, f, x::Number)
    dx = ForwardDiff.Dual(ForwardDiff.Dual(x, one(typeof(x))), one(typeof(x)))
    out = f(dx)
    val = out.value.value
    grad = out.partials.values[1].value
    hess = out.partials.values[1].partials[1]
    DiffResults.DiffResult(val, grad, hess)
end

function getgradient(res::DiffResults.ImmutableDiffResult)
    DiffResults.gradient(res)
end
function getgradient(res::DiffResults.MutableDiffResult)
    copy(DiffResults.gradient(res))
end
function getgradient(res::DiffResults.ImmutableDiffResult{1,T,Tuple{T}}) where T <: Number
    DiffResults.derivative(res)
end

function getgradient(res::DiffResults.ImmutableDiffResult{1,T,Tuple{T, T}}) where T <: Number
    DiffResults.derivative(res)
end

function gethessian(res::DiffResults.MutableDiffResult)
    copy(DiffResults.hessian(res))
end
function gethessian(res::DiffResults.ImmutableDiffResult)
    DiffResults.hessian(res)
end

function soptimize(f, x::Union{StaticVector{P,T}, TN, AbstractVector}; bto::BackTrackingOrder = Order2(), tol = 1e-8,
    updating = true, maxiter = 200) where {P,T, TN <: Number}
    res = setresult(x)
    ls = BackTracking()
    order = ordernum(bto)
    xinit = copy(x)
    x_new = copy(x)
    hx = gethessian(res)
    hold = initialh(x)::typeof(hx)
    jold = copy(x); s = copy(x)
    jx = jold
    @unpack c_1, ρ_hi, ρ_lo, iterations = ls
    iterfinitemax = -log2(eps(eltype(x)))
    sqrttol = sqrt(eps(Float64))
    α_0 = 1.
    f_calls = 0
    g_calls = 0
    for n = 1:maxiter
        ## Compute
        res = sethessian!(res, f, x); f_calls +=1; g_calls +=1; # Obtain gradient
        ϕ_0 = DiffResults.value(res)
        ## Check convergence
        isfinite(ϕ_0) || return StaticOptimizationResults(xinit, NaN*x,
        NaN, n, false, tol, f_calls, g_calls, jx, hx)
        jx = getgradient(res)
        norm(jx, Inf) < tol && return StaticOptimizationResults(xinit, x,
        ϕ_0, n, true, tol, f_calls, g_calls, jx, hx)
        hx = gethessian(res)
        ## Compute search directions
        s = -hx\jx
        dϕ_0 = dot(jx, s)
        if dϕ_0 >= 0. # If bad, reset search direction
            hx = hold
            s = -jx
            dϕ_0 = dot(jx, s)
        end
        ## Perform line search

        # Count the total number of iterations
        iteration = 0
        ϕx_0, ϕx_1 = ϕ_0, ϕ_0
        α_1, α_2 = α_0, α_0
        ϕx_1 = f(x + α_1*s); f_calls +=1;

        # Hard-coded backtrack until we find a finite function value
        iterfinite = 0
        while !isfinite(ϕx_1) && iterfinite < iterfinitemax
            iterfinite += 1
            α_1 = α_2
            α_2 = α_1/2
            ϕx_1 = f(x + α_2*s); f_calls += 1;
        end

        # Backtrack until we satisfy sufficient decrease condition
        while ϕx_1 > ϕ_0 + c_1 * α_2 * dϕ_0
            # If this part is reached we did not accept the initial linesearch
            # guess, so we will need to update on the next iterations
            # Increment the number of steps we've had to perform
            iteration += 1

            # Ensure termination
            if iteration > iterations
                error("Linesearch failed to converge, reached maximum iterations $(iterations).",
                α_2)
            end

            # Shrink proposed step-size:
            if order == 2 || iteration == 1
                # backtracking via quadratic interpolation:
                # This interpolates the available data
                #    f(0), f'(0), f(α)
                # with a quadractic which is then minimised; this comes with a
                # guaranteed backtracking factor 0.5 * (1-c_1)^{-1} which is < 1
                # provided that c_1 < 1/2; the backtrack_condition at the beginning
                # of the function guarantees at least a backtracking factor ρ.
                α_tmp = - (dϕ_0 * α_2^2) / ( 2 * (ϕx_1 - ϕ_0 - dϕ_0*α_2) )
            else
                div = 1. / (α_1^2 * α_2^2 * (α_2 - α_1))
                a = (α_1^2*(ϕx_1 - ϕ_0 - dϕ_0*α_2) - α_2^2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div
                b = (-α_1^3*(ϕx_1 - ϕ_0 - dϕ_0*α_2) + α_2^3*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div

                if isapprox(a, zero(a), atol = eps(Float64))
                    α_tmp = dϕ_0 / (2*b)
                else
                    # discriminant
                    d = max(b^2 - 3*a*dϕ_0, 0.)
                    # quadratic equation root
                    α_tmp = (-b + sqrt(d)) / (3*a)
                end
            end
            α_1 = α_2

            α_tmp = NaNMath.min(α_tmp, α_2*ρ_hi) # avoid too small reductions
            α_2 = NaNMath.max(α_tmp, α_2*ρ_lo) # avoid too big reductions

            # Evaluate f(x) at proposed position
            ϕx_0, ϕx_1 = ϕx_1, f(x + α_2*s); f_calls += 1;
        end
        alpha, fpropose = α_2, ϕx_1

        s = alpha*s
        x = x + s # Update x
        jold = copy(jx)
    end
    return StaticOptimizationResults(xinit, NaN*x,
    NaN, maxiter, false, tol, f_calls, g_calls, jx, hx)
end

function Base.show(io::IO, r::StaticOptimizationResults)
    @printf io "Results of Static Optimization Algorithm\n"
    @printf io " * Initial guess: [%s]\n" join(r.initial_x, ",")
    @printf io " * Minimizer: [%s]\n" join(r.minimizer, ",")
    @printf io " * Minimum: [%s]\n" join(r.minimum, ",")
    @printf io " * ∇f(x): [%s]\n" join(r.g, ",")
    @printf io " * Hf(x) is size [%s]\n" join(size(r.h), ",")
    @printf io " * Number of iterations: [%s]\n" join(r.iterations, ",")
    @printf io " * Number of function calls: [%s]\n" join(r.f_calls, ",")
    @printf io " * Number of gradient calls: [%s]\n" join(r.g_calls, ",")
    @printf io " * Converged: [%s]\n" join(r.g_converged, ",")
    return
end

### Modified Newton algorithm for root-finding
# I made this up, I do not know if it has nice convergence properties
# It first computes the value and derivative and performs a Newton step
# Then it computes the value and the potentially unnecessary derivative at the Newton point
# If this an improvement (closer to 0), the Newton step is accepted and it performs another Newton step
# If this is worse, it backtracks if the candidate is not finite
# Then it performs a Halley step using the function value at the candidate point
# If this is an improvement, the Halley step is accepted
# If it's still not an improvement, simple backtracking is done until it is
function snewton(f, x::Number; maxiter = 200, tol = 1e-8)
    res = DiffResults.DiffResult(x, (x,))
    iterfinitemax = -log2(eps(eltype(x)))
    α_0 = 1.
    res = ForwardDiff.derivative!(res, f, x) # Obtain gradient
    ϕ_0 = DiffResults.value(res)
    abs(ϕ_0) < tol && return (x = x, fx = ϕ_0, isroot = true, iter = 0)
    isfinite(ϕ_0) || return (x = NaN*x, fx = NaN, isroot = false, iter = 0)

    needsupdate = false
    for n = 1:maxiter
        if needsupdate
            res = ForwardDiff.derivative!(res, f, x) # Obtain gradient
            needsupdate = false
        end
        ϕ_0 = DiffResults.value(res)
        abs(ϕ_0) < tol && return (x = x, fx = ϕ_0, isroot = true, iter = n)
        isfinite(ϕ_0) || return (x = NaN*x, fx = NaN, isroot = false, iter = n)
        jx = DiffResults.derivative(res)
        x2 = x - ϕ_0/jx

        # Count the total number of iterations
        iteration = 0
        ϕx_0, ϕx_1 = ϕ_0, ϕ_0
        α_1, α_2 = α_0, α_0
        res = ForwardDiff.derivative!(res, f, x2) # Obtain gradient
        ϕx_1 = DiffResults.value(res)
        abs(ϕx_1) < tol && return (x = x2, fx = ϕx_1, isroot = true, iter = n)

        # Hard-coded backtrack until we find a finite function value
        iterfinite = 0
        x2old = x2
        while !isfinite(ϕx_1) && iterfinite < iterfinitemax
            needsupdate = true
            iterfinite += 1
            α_1 = α_2
            α_2 = α_1/2
            x2 = (1 - α_2)*x + α_2*x2old # convex combination
            ϕx_1 = f(x2)
        end

        α_h = 1.
        # Backtrack until we satisfy sufficient decrease condition
        if abs(ϕx_1) > abs(ϕx_0) # Closer to 0?
            needsupdate = true
            # Increment the number of steps we've had to perform
            # Interpolate available data using quadratic
            hx = 2(ϕx_1 - ϕx_0 - jx*(x2 - x) ) / (x2 - x)^2 # Taylor approximation
            x2 = x - 2ϕx_0*jx / (2jx^2 - ϕx_0*hx) # Halley step
            ϕx_1 = f(x2)
            iteration = 0
            # If still not closer to 0, do simple backtracking until we are
            while abs(ϕx_1) > abs(ϕx_0)
                iteration += 1
                α_1 = α_h
                α_2 = α_1/2
                x2 = (1 - α_2)*x + α_2*x2 # convex combination
                ϕx_1 = f(x2)
                iteration > 30 && error("Failed to converge")
            end
        end
        x = x2
    end
    return (x = NaN*x, fx = NaN, isroot = false, iter = maxiter)
end

### Copied from Tamas Papp on Discourse

"""
bisection(f, a, b; fa = f(a), fb = f(b), ftol, wtol)

Bisection algorithm for finding the root ``f(x) ≈ 0`` within the initial bracket
    `[a,b]`.

    Returns a named tuple

    `(x = x, fx = f(x), isroot = ::Bool, iter = ::Int, ismaxiter = ::Bool)`.

    Terminates when either

    1. `abs(f(x)) < ftol` (`isroot = true`),
    2. the width of the bracket is `≤wtol` (`isroot = false`),
    3. `maxiter` number of iterations is reached. (`isroot = false, maxiter = true`).

    which are tested for in the above order. Therefore, care should be taken not to make `wtol` too large.

"""
function bisection(f, a::Real, b::Real;
                   ftol = √eps(), wtol = 0., maxiter = 100)
    a, b = float(a), float(b)
    fa, fb = f(a), f(b)
    fa * fb ≤ 0 || error("Not a bracket")
    (isfinite(a) && isfinite(b)) || error("Not finite")
    _bisection(f, a, b, fa, fb, ftol, wtol, maxiter)
end

function _bisection(f, a, b, fa, fb, ftol, wtol, maxiter)
    iter = 0
    abs(fa) < ftol && return (x = a, fx = fa, isroot = true, iter = iter, ismaxiter = false)
    abs(fb) < ftol && return (x = b, fx = fb, isroot = true, iter = iter, ismaxiter = false)
    while true
        iter += 1
        m = middle(a, b)
        fm = f(m)
        abs(fm) < ftol && return (x = m, fx = fm, isroot = true, iter = iter, ismaxiter = false)
        abs(b-a) ≤ wtol && return (x = m, fx = fm, isroot = false, iter = iter, ismaxiter = false)
        if fa * fm > 0
            a, fa = m, fm
        else
            b, fb = m, fm
        end
        iter == maxiter && return (x = m, fx = fm, isroot = false, iter = iter, ismaxiter = true)
    end
end

sroot(f, x::Number) = snewton(f, x)
sroot(f, x::Tuple{T, T}) where T <: Number = @inbounds bisection(f, x[1], x[2])
function sroot(f, x::AbstractVector;
    updating = true, tol = 1e-8, maxiter = 200)
    f2(s) = sum(x -> x^2, f(s))
    soptimize(f2, x, updating = updating, tol = tol, maxiter = maxiter )
end

#=
Algorithm:
Keep the actual hessian around
Restrict it to the active variables
Only then do the inversion (linear system)

Inverting the restricted hessian will be identical to just solving a lower dimensional problem

To keep dimensions constant, zero out the inactive part of the hessian, but fill the
diagonal with 1. This will be block diagonal and the inverse of the active region
will be the inverse of the hessian of the lower dimensional problem. The inverse
of the inactive region will just be 1 on the diagonal.
=#

function constrained_soptimize(f, x::Union{StaticVector{P,T}, TN, AbstractVector};
    bto::BackTrackingOrder = Order2(), tol = 1e-7, lower = -Inf*x,
    upper = Inf*x, maxiter = 200) where {P,T, TN <: Number}
    res = setresult(x)
    ls = BackTracking()
    order = ordernum(bto)
    xinit = copy(x)
    x_new = copy(x)
    hx = gethessian(res)
    hold = initialh(x)::typeof(hx)
    jold = copy(x); s = copy(x)
    xold = copy(x)
    jx = copy(jold)
    @unpack c_1, ρ_hi, ρ_lo, iterations = ls
    iterfinitemax = -log2(eps(eltype(x)))
    sqrttol = sqrt(eps(Float64))
    α_0 = 1.
    f_calls = 0
    g_calls = 0
    clamp.(x, lower, upper) == x || error("Initial guess not in the feasible region")
    for n = 1:maxiter
        ## Compute the gradient if needed
        res = sethessian!(res, f, x); f_calls +=1; g_calls +=1; # Obtain hessian
        ϕ_0 = DiffResults.value(res)

        ## Check convergence
        isfinite(ϕ_0) || return StaticOptimizationResults(xinit, NaN*x,
        NaN, n, false, tol, f_calls, g_calls, jx, hx)
        jx = getgradient(res)

        # Binding set 1
        # true if free

        s1l = .!( ((x .<= lower) .& (jx .>= 0)) .| ((x .>= upper) .& (jx .<= 0)) )
        hx = gethessian(res)

        function cdiag(x, c, i)
            if c
                return x
            else
                if i == 1.
                    return 1.
                else
                    return 0.
                end
            end
        end
        isbinding(i, j) = i & j
        # Binding set 2
        binding = isbinding.(s1l, s1l')
        #sizex = size(x)[1]
        Ix = initialh(x) # An identity matrix of correct type
        Hbar = cdiag.(hx, binding, Ix)

        Sbargrad = (Hbar\jx)
        s2l = .!( ((x .<= lower) .& (Sbargrad .> 0)) .| ((x .>= upper) .& (Sbargrad .< 0)) )

        sl = s1l .& s2l
        binding = isbinding.(sl, sl')
        Hhat = cdiag.(hx, binding, Ix)

        # set jxc to 0 if var in binding set
        jxc = jx .* sl
        norm(jxc, Inf) < tol && return StaticOptimizationResults(xinit, x,
        ϕ_0, n, true, tol, f_calls, g_calls, jx, hx)

        s = (-Hhat\jx) .* sl
        dϕ_0 = dot(jx, s)

        if dϕ_0 >= 0. # If bad, reset search direction
            hx = hold
            binding = isbinding.(s1l, s1l')
            Shat = hx .* binding
            s = -Shat*jx
            dϕ_0 = dot(jx, s)
        end
        ## Perform line search

        # Count the total number of iterations
        iteration = 0
        ϕx_0, ϕx_1 = ϕ_0, ϕ_0
        α_1, α_2 = α_0, α_0
        ϕx_1 = f(clamp.(x + α_1*s, lower, upper)); f_calls +=1;

        # Hard-coded backtrack until we find a finite function value
        iterfinite = 0
        while !isfinite(ϕx_1) && iterfinite < iterfinitemax
            iterfinite += 1
            α_1 = α_2
            α_2 = α_1/2
            ϕx_1 = f(clamp.(x + α_2*s, lower, upper)); f_calls += 1;
        end
        # Backtrack until we satisfy sufficient decrease condition
        while ϕx_1 > ϕ_0 + c_1 * α_2 * dϕ_0
            # If this part is reached we did not accept the initial linesearch
            # guess, so we will need to update on the next iterations
            # Increment the number of steps we've had to perform
            iteration += 1

            # Ensure termination
            if iteration > iterations
                error("Linesearch failed to converge, reached maximum iterations $(iterations).",
                α_2)
            end

            # Shrink proposed step-size:
            if order == 2 || iteration == 1
                # backtracking via quadratic interpolation:
                # This interpolates the available data
                #    f(0), f'(0), f(α)
                # with a quadractic which is then minimised; this comes with a
                # guaranteed backtracking factor 0.5 * (1-c_1)^{-1} which is < 1
                # provided that c_1 < 1/2; the backtrack_condition at the beginning
                # of the function guarantees at least a backtracking factor ρ.
                α_tmp = - (dϕ_0 * α_2^2) / ( 2 * (ϕx_1 - ϕ_0 - dϕ_0*α_2) )
            else
                div = 1. / (α_1^2 * α_2^2 * (α_2 - α_1))
                a = (α_1^2*(ϕx_1 - ϕ_0 - dϕ_0*α_2) - α_2^2*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div
                b = (-α_1^3*(ϕx_1 - ϕ_0 - dϕ_0*α_2) + α_2^3*(ϕx_0 - ϕ_0 - dϕ_0*α_1))*div

                #if norm(a) <= eps(Float64) + sqrttol*norm(a)
                if isapprox(a, zero(a))
                    α_tmp = dϕ_0 / (2*b)
                else
                    # discriminant
                    d = max(b^2 - 3*a*dϕ_0, 0.)
                    # quadratic equation root
                    α_tmp = (-b + sqrt(d)) / (3*a)
                end
            end
            α_1 = α_2

            α_tmp = NaNMath.min(α_tmp, α_2*ρ_hi) # avoid too small reductions
            α_2 = NaNMath.max(α_tmp, α_2*ρ_lo) # avoid too big reductions

            # Evaluate f(x) at proposed position
            ϕx_0, ϕx_1 = ϕx_1, f(clamp.(x + α_2*s, lower, upper)); f_calls += 1;
        end
        alpha, fpropose = α_2, ϕx_1

        s = alpha*s
        xnew = clamp.(x + s, lower, upper) # Update x
        s = xnew - x
        xold = x
        x = xnew
        jold = copy(jx)
    end
    return StaticOptimizationResults(xinit, x,
    NaN, maxiter, false, tol, f_calls, g_calls, jx, hx)
end

max_stars_repo_name: aaowens/StaticOptim

max_issues_repo_name: aaowens/StaticOptim

max_forks_repo_name: aaowens/StaticOptim


text:
function dim_to_indices(arr::IndexedDimsArray{T, N}, inds::Tuple{Vararg{Any, N}}) where {T, N}
    return map(unwrap_index, arr.indexes, inds)
end

unwrap_index(acc::AcceleratedVector, ind) = findall(in(ind), acc)
unwrap_index(acc::AcceleratedVector, ::Colon) = Colon()
unwrap_index(acc::AcceleratedVector, ind::BypassIndex) = ind.index

function dim_getindex(arr::IndexedDimsArray{T, N}, inds::Vararg{Any, N}) where {T, N}
    # return Base.getindex(parent(arr), dim_to_indices(arr, inds)...)
    return invoke(getindex, Tuple{AbstractArray, Vararg{Any, N}}, arr, dim_to_indices(arr, inds)...)
end

function dim_setindex!(arr::IndexedDimsArray{T, N}, val, inds::Vararg{Any, N}) where {T, N}
    return Base.setindex!(parent(arr), val, dim_to_indices(arr, inds)...)
end

function Base.getindex(arr::IndexedDimsArray{T, N}, inds::Vararg{Any, N}) where {T, N}
    return dim_getindex(arr, inds...)
end

function Base.setindex!(arr::IndexedDimsArray{T, N}, val, inds::Vararg{Any, N}) where {T, N}
    return dim_setindex!(arr, inds...)
end


text:
#=
test_xml_preprocessor:
- Julia version: 1.1.0
- Author: bramb
- Date: 2019-03-21
=#

using Test, HyperCollate

test_subst_wrapping(xml::String, expected::String) = @test add_subst(xml) == expected

macro debug_on()
    return :(ENV["JULIA_DEBUG"] = "all")
end

macro debug_off()
    return :(ENV["JULIA_DEBUG"] = "")
end

@testset "automatic subst wrapping" begin
    @testset "wrapping non-nested add/del in subst" begin
        @testset "hoe zoet..." begin
            xml = """
            <xml>Hoe zoet moet nochtans zijn dit werven om<add>trachten naar</add> een vrouw,
            de ongewisheid vóór de liefelijke toestemming!</xml>
            """

            expected = """
            <xml>Hoe zoet moet nochtans zijn dit <subst>werven om<add>trachten naar</add></subst> een vrouw,
            de ongewisheid vóór de liefelijke toestemming!</xml>"""

            test_subst_wrapping(xml,expected)
        end

        @testset "simple del/add" begin
            xml = "<x>bla1 bla2<add>bla3</add> bla4</x>"
            expected = "<x>bla1 <subst>bla2<add>bla3</add></subst> bla4</x>"
            test_subst_wrapping(xml,expected)
        end

        @testset "don't wrap solitary del/add" begin
            xml = "<x>bla1 bla2 something <add>bla3</add> bla4</x>"
            expected = "<x>bla1 bla2 something <add>bla3</add> bla4</x>"
            test_subst_wrapping(xml,expected)
        end

        @testset "ignore whitespace between del/add" begin
            xml = "<x>bla1 bla2\n <add>bla3</add> bla4</x>"
            expected = "<x>bla1 <subst>bla2<add>bla3</add></subst> bla4</x>"
            test_subst_wrapping(xml,expected)
        end

        @testset "wrap del/del/add combo" begin
            xml = "<x>something not thisor this<add>but this</add> something else</x>"
            expected = "<x>something <subst>not thisor this<add>but this</add></subst> something else</x>"
            test_subst_wrapping(xml,expected)
        end
    end

    @testset "wrapping nested add/del in subst" begin
        @testset "brulez 01r nested del" begin
#             @debug_on()
            xml = """
            <xml><s>...weinig van pas komen
            zoo, o.m. in de sexueele opvoeding van den troop<add place="supralinear" hand="#RB" rend="grey pencil" resp="#EB">n</add>o<add place="supralinear" rend="grey pencil" hand="#RB" resp="#EB">p</add>volger...</s></xml>
            """
            expected = """
<xml><s>...weinig van pas komen
zoo, o.m. in de sexueele opvoeding van den troo<subst>p<add place="supralinear" hand="#RB" rend="grey pencil" resp="#EB">n</add></subst>o<add place="supralinear" hand="#RB" rend="grey pencil" resp="#EB">p</add>volger...</s></xml>"""
            test_subst_wrapping(xml,expected)
#             println(add_subst(xml))
#             @debug_off()
        end

        @testset "brulez 01r nested del 2" begin
#             @debug_on()
            xml = """<s>weinig van pas komen
,<add>:</add>
zoo, o.m. in de sexueele opvoeding van den troop<add>n</add>o<add>p</add>volger...
</s>"""
            expected = """<s>weinig van pas komen
<subst>,<add>:</add></subst>
zoo, o.m. in de sexueele opvoeding van den troo<subst>p<add>n</add></subst>o<add>p</add>volger...</s>"""
            test_subst_wrapping(xml,expected)
#             @debug_off()
        end

        @testset "del/add nested in add" begin
            xml = "<x>bla1 bla2\n <add>bla3 something<add>word</add></add> bla4</x>"
            expected = "<x>bla1 <subst>bla2<add>bla3 <subst>something<add>word</add></subst></add></subst> bla4</x>"
            test_subst_wrapping(xml,expected)
        end

        @testset "del/add nested in solitary del" begin
            xml = "
De te streng doorgedreven rationalisatie van zijn prinsenjeugd had dit met hem<add>zich</add> meegebracht.
"
            expected = "
De te streng doorgedreven rationalisatie van zijn prinsenjeugd had dit met <subst>hem<add>zich</add></subst> meegebracht.
"
            test_subst_wrapping(xml,expected)
        end

        @testset "Old Miss: del in add" begin
#             @debug_on()
            xml = "<s>Old Miss Hare<add><unclear>Scovell</unclear></add><add>McGlone<add>McGlone</add></add> always sings at this hour.</s>"
            expected = "<s>Old Miss <subst>Hare<add><unclear>Scovell</unclear></add><add><subst>McGlone<add>McGlone</add></subst></add></subst> always sings at this hour.</s>"
            test_subst_wrapping(xml,expected)
#             @debug_off()
        end

    end

    # @testset "handling milestones" begin
    #     xml = "<x>word <milestone/> word</x>"
    #     test_subst_wrapping(xml,xml)
    # end
end


text:
function run!(::Experiment{:Hypothetical}, models)
    results = DataFrame(id = Int[], action = String[], A = Int[], B = Int[])
    for id in eachindex(models)
        m = models[id]
        trays = [Tray(i) for i in 1:4] # C, D, A, B
        for day in 1:9
            cid = floor(Int, ((day-1)) % 2 + 1)
            add!(m, trays[cid]) # compartment to spend the day from 10am onwards
            add!(m, MaintenanceDiet)
            wait!(m, 6.75u"hr")
            remove!(m, Any)
            if day == 9
                add!(m, trays[3:4])
                add!(m, Kibble, 10) # test
                wait!(m, 15u"minute")
                remove!(m, Any)
                push!(results, [id; "cache"; countcache.(trays[3:4])])
                break # end of experiment
            else
                add!(m, trays[3:4])
                add!(m, PowderedKibble, 10)
                wait!(m, 15u"minute") # getting used to seeing them in the evening
                remove!(m, Any)
                add!(m, MaintenanceDiet)
                wait!(m, 3u"hr")
            end
            remove!(m, Any)
            wait!(m, 12u"hr") # morning
            add!(m, trays[cid == 1 ? 3 : 4]) # compartment to wake up hungry
            wait!(m, 1u"hr")
            remove!(m, Any)
            wait!(m, 1u"hr")
            add!(m, MaintenanceDiet)
        end
    end
    results
end


text:
function mvnrnd(mu::Array{T,1},sigma::Array{T,2},cases::Int64) where T <: Real
      TM = chol(sigma)
      r = broadcast( + ,  TM * randn(size(TM,1),cases) , mu )
end

function sumsq(xx::Vector{T}) where T <: Real
      return sum( xx.^2 )
end
function sumsq(xx::Array{T,2}) where T <: Real
      return sum( xx.^2 , 1)
end

using Distributions.pdf
using Distributions.MvNormal

function mvnpdf(mu::Array{T,1},sigma::Array{T,2},x::Array{T,2}) where T <: Real
      return pdf(MvNormal(mu,sigma),x)
end
function mvnpdf(mu::Array{T,1},sigma::Array{T,2},x::Array{T,1}) where T <: Real
      return pdf(MvNormal(mu,sigma),x)
end

#function mvnpdf_kaputt(mu::Array{Float64,1},sigma::Array{Float64,2},x::Array{Float64,2})
#      p = size(x,2)
#      r = chol(sigma)
#      #y = (2*pi)^(-1/2) * exp.(-sumsq(  (r'\ broadcast(-,x,mu)')'  )/2) ./ prod(diag(r))
#      pdf = (2*pi)^(-p/2) * exp(-sumsq(( broadcast(-,x,mu) )/r,2)/2) / prod(diag(r))
#end

struct GaussParam
      mu::Vector{Float64}
      sigma::Array{Float64,2}
end

struct GaussMixture
      mix::Vector{GaussParam}
      weights::Vector{Float64}
end
function GaussMixture(mix::Vector{GaussParam})
      weights = ones(length(mix)) * (1./length(mix))
      return GaussMixture(mix,weights)
end

function draw_gaussmixture(gm::GaussMixture , n::Int64)
      d   = length(gm.mix[1].mu)
      x   = NaN * ones(d,n)
      xsp = sample_proportional( gm.weights , n )
      for zi=1:length(gm.weights)
            x[:,xsp.==zi] = mvnrnd( gm.mix[zi].mu,gm.mix[zi].sigma, sum(xsp.==zi) )
      end
      return x
end

function pdf_gaussmixture_nonnormalized(gm::GaussMixture , x::Array{T,2}) where T <: Real
      p = zeros(size(x,2))
      for zi=1:length(gm.mix)
            p += mvnpdf( gm.mix[zi].mu , gm.mix[zi].sigma , x )
      end
      return p
end
function pdf_gaussmixture_nonnormalized(gm::GaussMixture , x::Array{T,1}) where T <: Real
      return pdf_gaussmixture_nonnormalized(gm,x[:,:])[1]
end

function sample_proportional( lambdas , n )
#fSampleProportional( lambdas , n ), samples n times integers 1:length(lambdas), according
#to the weights provided in lambdas.
    #lambdas = lambdas[:]'
    lambdas = lambdas ./ sum(lambdas)
    csl     = cumsum(lambdas)
    s       = rand(n,1)
    li      = broadcast( >= , csl , s' )
    fmres   = findmax( li.!=0, 1 )
    li       = ind2sub( size(li) , fmres[2][:] )[1]
end

function unifsamp_box( bounds::Array{T,2} , n::Int64 ) where T <: Real
    return broadcast( + , bounds[:,1] , broadcast( * , (bounds[:,2]-bounds[:,1]) , rand(size(bounds,1),n) ) )
end

function get_bounding_box( lb::Vector{T} , ub::Vector{T} ) where T <: Real
      G = [eye(length(lb));-eye(length(lb))]
      h = [ub;-lb]
      return (G,h)
end

function xygrid(xlim,ylim,xres,yres)
      grid_x = [ xlim[1]+ (0.5+i)*((xlim[2]-xlim[1])/(xres+1)) for i=1:xres, j=1:yres ]
      grid_y = [ ylim[1]+ (0.5+j)*((ylim[2]-ylim[1])/(yres+1)) for i=1:xres, j=1:yres ]
      return (grid_x,grid_y)
end


text:
macro c(ret_type, func, arg_types, lib)
  local args_in = Any[ symbol(string('a',x)) for x in 1:length(arg_types.args) ]
  quote
    $(esc(func))($(args_in...)) = ccall( ($(string(func)), $(Expr(:quote, lib)) ), $ret_type, $arg_types, $(args_in...) )
  end
end

macro ctypedef(fake_t,real_t)
  quote
    typealias $fake_t $real_t
  end
end

type aes_key_st
  rd_key::Void
  rounds::Int32
end
@ctypedef AES_KEY aes_key_st
@ctypedef DES_cblock Void
@ctypedef const_DES_cblock Void
type DES_ks
  ks::Void
end
@ctypedef DES_key_schedule DES_ks
type MD5state_st
  A::Uint32
  B::Uint32
  C::Uint32
  D::Uint32
  Nl::Uint32
  Nh::Uint32
  data::Void
  num::Uint32
end
@ctypedef MD5_CTX MD5state_st
type MD4state_st
  A::Uint32
  B::Uint32
  C::Uint32
  D::Uint32
  Nl::Uint32
  Nh::Uint32
  data::Void
  num::Uint32
end
@ctypedef MD4_CTX MD4state_st
type mdc2_ctx_st
  num::Uint32
  data::Void
  h::DES_cblock
  hh::DES_cblock
  pad_type::Int32
end
@ctypedef MDC2_CTX mdc2_ctx_st
type SHAstate_st
  h0::Uint32
  h1::Uint32
  h2::Uint32
  h3::Uint32
  h4::Uint32
  Nl::Uint32
  Nh::Uint32
  data::Void
  num::Uint32
end
@ctypedef SHA_CTX SHAstate_st
type SHA256state_st
  h::Void
  Nl::Uint32
  Nh::Uint32
  data::Void
  num::Uint32
  md_len::Uint32
end
@ctypedef SHA256_CTX SHA256state_st
type SHA512state_st
  h::Void
  Nl::Uint64
  Nh::Uint64
  u::
  num::Uint32
  md_len::Uint32
end
@ctypedef SHA512_CTX SHA512state_st
type rand_meth_st
  seed::Ptr{Void}
  bytes::Ptr{Void}
  cleanup::Ptr{Void}
  add::Ptr{Void}
  pseudorand::Ptr{Void}
  status::Ptr{Void}
end


text:
@testset "Matrix" begin

    @testset "Matrix: Determinant" begin
        M1 = [1 0; 2 2]
        M2 = rand(3,3)
        M3 = rand(4,4)

        @test determinant(M1) == det(M1)
        @test round(determinant(M2),digits = 4) == round(det(M2),digits = 4)
        @test round(determinant(M3),digits = 4) == round(det(M3),digits = 4)
    end

    @testset "Matrix: LU Decompose" begin
        mat = [
            2  -1 -2;
            -4   6  3;
            -4  -2  8
        ]

        L = [
            1  0 0;
            -2  1 0;
            -2 -1 1
        ]
        U = [
            2 -1 -2;
            0  4 -1;
            0  0  3
        ]

        @test lu_decompose(mat) == (L,U)
    end

    @testset "Matrix: rotation matrix" begin
        theta = pi/6
        a = [1;0]
        b = [0;1]
        @test rotation_matrix(theta)*a == [cos(theta);sin(theta)]
        @test rotation_matrix(theta)*b == [-sin(theta);cos(theta)]
    end

    @testset "Matrix: Gauss-Jordan Elimination" begin
        M1 = [1 2 3; 4 5 6];
        M2 = [1 2 3; 4 8 12];
        M3 = [1 2; 4 8];
        M4 = [1 4 8  1  4
              4 5 6  8  11
              1 3 2  4  8
              4 5 67 23 0];
        M5 = Float64[1 4 8  1  4
                     4 5 6  8  11
                     1 3 2  4  8
                     4 5 67 23 0];
        M6 = [1 4 8
              4 0 6
              1 3 2];

        R3 =  [ 1.0  0.0  0.0  0.0  -1.10637
                0.0  1.0  0.0  0.0   1.89743
                0.0  0.0  1.0  0.0  -0.444913
                0.0  0.0  0.0  1.0   1.07598];

        @test gauss_jordan(M1) == Float64[1 0 -1; 0 1 2]
        @test_throws AssertionError gauss_jordan(M2)
        @test_throws AssertionError gauss_jordan(M3)
        @test isapprox(gauss_jordan(M4), R3, atol=1e-5)
        @test isapprox(gauss_jordan(M5), R3, atol=1e-5)
        @test isapprox(gauss_jordan(M6), Float64[1 0 0; 0 1 0; -0 -0  1], atol=1e-5)
    end

end

max_stars_repo_name: frankschmitt/Julia

max_issues_repo_name: frankschmitt/Julia

max_forks_repo_name: frankschmitt/Julia


text:
using Base.Test
import Utils

@test Utils.allsame([1,1]) == true
@test Utils.allsame([1,2,1]) == false


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/constraints.jl
# Defines constraints supported by optimizer and how to store them.
#############################################################################

# Sets used in general constraints
const INEQ_SETS = Union{LT, GT, ET}
const VAR_SETS = Union{LT, GT, ET, ZO, MOI.Integer}

##### Utilities for checking that JuMP model contains variables used in expression
function check_inbounds!(m::Optimizer, vi::VI)
    if !(1 <= vi.value <= m._input_problem._variable_count)
        error("Invalid variable index $vi. ($(m._input_problem._variable_count) variables in the model.)")
    end
    return nothing
end
check_inbounds!(m::Optimizer, f::SAF) = foreach(x -> check_inbounds!(m, x.variable), f.terms)
check_inbounds!(m::Optimizer, f::VECOFVAR) = foreach(x -> check_inbounds!(m, x), f.variables)
function check_inbounds!(m::Optimizer, f::SQF)
    foreach(x -> check_inbounds!(m, x.variable), f.affine_terms)
    for term in f.quadratic_terms
        check_inbounds!(m, term.variable_1)
        check_inbounds!(m, term.variable_2)
    end
    return nothing
end

MOI.supports_constraint(::Optimizer, ::Type{VI}, ::Type{S}) where S <: VAR_SETS = true
MOI.supports_constraint(::Optimizer,::Type{T},::Type{S}) where {T<:Union{SAF,SQF},S<:INEQ_SETS} = true

MOI.is_valid(m::Optimizer, v::VI) = (1 <= v.value <= m._input_problem._variable_count)
MOI.is_valid(m::Optimizer, c::CI{VI,S}) where S <: VAR_SETS = (1 <= c.value <= m._input_problem._variable_count)
MOI.is_valid(m::Optimizer, c::CI{F,S}) where {F<:Union{SAF,SQF}, S<:INEQ_SETS} = (1 <= c.value <= length(_constraints(m,F,S)))

MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{VI,S}) where S<:VAR_SETS = length(_constraints(m,VI,S))
MOI.get(m::Optimizer, ::MOI.NumberOfConstraints{F,S}) where {F<:Union{SAF,SQF},S<:INEQ_SETS} = length(_constraints(m,F,S))

MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{VI,S}) where S<:VAR_SETS = collect(keys(_constraints(m,VI,S)))
MOI.get(m::Optimizer, ::MOI.ListOfConstraintIndices{F,S}) where {F<:Union{SAF,SQF},S<:INEQ_SETS} = collect(keys(_constraints(m,F,S)))

MOI.add_variable(m::Optimizer) = VI(m._input_problem._variable_count += 1)

function MOI.add_constraint(m::Optimizer, f::F, s::S) where {F<:Union{SAF,SQF},S<:INEQ_SETS}
    check_inbounds!(m, f)
    ci = CI{F, S}(m._input_problem._constraint_count += 1)
    _constraints(m, F, S)[ci] = (f, s)
    return ci
end
function MOI.add_constraint(m::Optimizer, f::VI, s::S) where S<:VAR_SETS
    check_inbounds!(m, f)
    ci = CI{VI,S}(f.value)
    _constraints(m,VI,S)[ci] = (f, s)
    return ci
end

result_index_1_error(v::T) where T = throw(MOI.ResultIndexBoundsError{T}(v, 1))
function MOI.get(model::Optimizer, v::MOI.VariablePrimal, vi::MOI.VariableIndex)
    check_inbounds!(model, vi)
    (v.result_index != 1) && result_index_1_error(v)
    return model._global_optimizer._continuous_solution[vi.value]
end
function MOI.get(m::Optimizer{R,S,T}, v::MOI.ConstraintPrimal, c::CI{VI,<:Any}) where {R,S,T}
    (v.result_index != 1) && result_index_1_error(v)
    return MOI.get(m, MOI.VariablePrimal(), MOI.VariableIndex(c.value))
end
function MOI.get(m::Optimizer{R,Q,T}, v::MOI.ConstraintPrimal, c::CI{F,S}) where {R,Q,T,F,S}
    (v.result_index != 1) && result_index_1_error(v)
    return _constraint_primal(m._global_optimizer, F, S)[c]
end

MOI.get(m::Optimizer, ::MOI.ConstraintFunction, c::CI{F,S}) where {F,S} = _constraints(m,F,S)[c][1]
MOI.get(m::Optimizer, ::MOI.ConstraintSet, c::CI{F,S}) where {F,S} = _constraints(m,F,S)[c][2]

function MOI.empty!(m::Optimizer{R,S,T}) where {R,S,T}

    MOI.empty!(m.subsolver_block)
    MOI.empty!(m._global_optimizer)
    m._input_problem = InputProblem()
    m._working_problem = ParsedProblem()

    m._termination_status_code = MOI.OPTIMIZE_NOT_CALLED
    m._result_status_code = MOI.OTHER_RESULT_STATUS
    m._run_time = 0.0
    m._objective_value  = -Inf
    m._objective_bound  =  Inf
    m. _relative_gap     = Inf
    m._iteration_count  = 0
    m._node_count       = 0

    return nothing
end

function MOI.is_empty(m::Optimizer{R,S,T}) where {R,S,T}

    flag = true
    flag &= MOI.is_empty(m._global_optimizer)
    flag &= isempty(m._input_problem)
    flag &= isempty(m._working_problem)
    flag &= isempty(m.subsolver_block)
    flag &= m._termination_status_code == MOI.OPTIMIZE_NOT_CALLED
    flag &= m._result_status_code == MOI.OTHER_RESULT_STATUS

    # set constructor reset on empty! and  to zero in initial parse! in parse.jl
    flag &= iszero(m._run_time)
    flag &= iszero(m._iteration_count)
    flag &= iszero(m._node_count)
    flag &= m._objective_value == -Inf
    flag &= m._objective_bound ==  Inf
    flag &= m. _relative_gap   == Inf

    return flag
end

MOI.supports_incremental_interface(m::Optimizer) = true
MOI.copy_to(model::Optimizer, src::MOI.ModelLike) = MOIU.default_copy_to(model, src)

#####
##### Set & get attributes of model
#####
MOI.supports(::Optimizer, ::MOI.Silent) = true

function MOI.set(m::Optimizer, s::MOI.Silent, value)
    push!(m._optimizer_attributes_set, s)
    unique!(m._optimizer_attributes_set)
    if value
        m._parameters.verbosity = 0
        m._parameters.log_on = false
    else
        m._parameters.verbosity = 1
    end
    return
end

function MOI.set(m::Optimizer, s::MOI.TimeLimitSec, ::Nothing)
    push!(m._optimizer_attributes_set, s)
    unique!(m._optimizer_attributes_set)
    m._parameters.time_limit = Inf
end
function MOI.set(m::Optimizer, s::MOI.TimeLimitSec, v::Int)
    push!(m._optimizer_attributes_set, s)
    unique!(m._optimizer_attributes_set)
    m._parameters.time_limit = v
end
function MOI.set(m::Optimizer, s::MOI.TimeLimitSec, v::Float64)
    push!(m._optimizer_attributes_set, s)
    unique!(m._optimizer_attributes_set)
    m._parameters.time_limit = v
end

MOI.get(m::Optimizer, ::MOI.ListOfOptimizerAttributesSet) = m._optimizer_attributes_set

function MOI.get(m::Optimizer, ::MOI.ListOfConstraintTypesPresent)
    constraint_types = []
    for S in (ZO, MOI.Integer)
        if MOI.get(m, MOI.NumberOfConstraints{VI,S}()) > 0
            push!(constraint_types, (VI,S))
        end
    end
    for S in (LT, GT, ET), T in (VI, SAF, SQF)
        if MOI.get(m, MOI.NumberOfConstraints{T,S}()) > 0
            push!(constraint_types, (VI,S))
        end
    end
    return constraint_types
end

MOI.get(m::Optimizer, v::MOI.ObjectiveValue) = !isone(v.result_index) ? result_index_1_error(v) : m._objective_value
MOI.get(m::Optimizer, v::MOI.PrimalStatus) = !isone(v.result_index) ? MOI.NO_SOLUTION : m._result_status_code
MOI.get(m::Optimizer, ::MOI.DualStatus) = MOI.NO_SOLUTION
MOI.get(m::Optimizer, ::MOI.ObjectiveBound) = m._objective_bound
MOI.get(m::Optimizer, ::MOI.NumberOfVariables) = m._input_problem._variable_count
MOI.get(m::Optimizer, ::MOI.SolverName) = "EAGO: Easy Advanced Global Optimization"
MOI.get(m::Optimizer, ::MOI.SolverVersion) = "0.7.1"
MOI.get(m::Optimizer, ::MOI.TerminationStatus) = m._termination_status_code
MOI.get(m::Optimizer, ::MOI.SolveTimeSec) = m._run_time
MOI.get(m::Optimizer, ::MOI.NodeCount) = m._node_count
MOI.get(m::Optimizer, ::MOI.ResultCount) = (m._result_status_code === MOI.FEASIBLE_POINT) ? 1 : 0
MOI.get(m::Optimizer, ::MOI.TimeLimitSec) = m._parameters.time_limit
MOI.get(m::Optimizer, ::MOI.Silent) = m._parameters.verbosity == 0
MOI.get(m::Optimizer, ::MOI.ListOfVariableIndices) = [VI(i) for i = 1:m._input_problem._variable_count]

function MOI.get(m::Optimizer, ::MOI.RelativeGap)
    b = MOI.get(m, MOI.ObjectiveBound())
    v = MOI.get(m, MOI.ObjectiveValue())
    return relative_gap(b,v)
end

_to_sym(d) = error("EAGO only supports raw parameters with Symbol or String names.")
_to_sym(d::String) = Symbol(d)
_to_sym(d::Symbol) = d
function MOI.get(m::Optimizer, p::MOI.RawOptimizerAttribute)
    s = _to_sym(p.name)
    s in EAGO_PARAMETERS ? getfield(m._parameters, s) : getfield(m, s)
end

raw_param_name(p::MOI.RawOptimizerAttribute) = _to_sym(p.name)
raw_param_name(p) = nothing
function MOI.set(m::Optimizer, p::MOI.RawOptimizerAttribute, x)
    inds = findall(x -> raw_param_name(x) == raw_param_name(p), m._optimizer_attributes_set)
    deleteat!(m._optimizer_attributes_set, inds)
    push!(m._optimizer_attributes_set, p)
    s = _to_sym(p.name)
    if (s == :relaxed_optimizer) || (s == :upper_optimizer)
        setfield!(m, s, Incremental(x))
    else
        s in EAGO_PARAMETERS ? setfield!(m._parameters, s, x) : setfield!(m, s, x)
    end
end
MOI.get(m::Optimizer, p::MOI.RawStatusString) = string(m._global_optimizer._end_state)

#####
##### Support, set, and evaluate objective functions
#####
MOI.supports(::Optimizer, ::MOI.TimeLimitSec) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveSense) = true
MOI.supports(::Optimizer, ::MOI.ObjectiveFunction{F}) where {F <: Union{VI, SAF, SQF}} = true

function MOI.set(m::Optimizer, ::MOI.NLPBlock, nlp_data::MOI.NLPBlockData)
    if nlp_data.has_objective
        m._input_problem._objective = nothing
    end
    m._input_problem._nlp_data = nlp_data
end

function MOI.set(m::Optimizer, ::MOI.ObjectiveFunction{T}, f::T) where T <: Union{VI,SAF,SQF}
    check_inbounds!(m, f)
    m._input_problem._objective = f
end
MOI.get(m::Optimizer, ::MOI.ObjectiveFunction{T}) where T <: Union{VI,SAF,SQF} = m._input_problem._objective
MOI.get(m::Optimizer, ::MOI.ObjectiveFunctionType) = typeof(m._input_problem._objective)

MOI.set(m::Optimizer, ::MOI.ObjectiveSense, s::MOI.OptimizationSense) = m._input_problem._optimization_sense = s
MOI.get(m::Optimizer, ::MOI.ObjectiveSense) = m._input_problem._optimization_sense

max_stars_repo_name: MatthewStuber/EAGO

max_issues_repo_name: MatthewStuber/EAGO

max_forks_repo_name: MatthewStuber/EAGO


text:
## ColorMap in Config for Imager
using Colors,ColorSchemes,YAML

cs = [RGB(HSL(360*i/255,1,0.5)) for i in 0:255]
cs = ColorSchemes.turbo.colors

cm = Dict("ColorMap"=>map(c->round.(Int,255*[c.r,c.g,c.b]),cs))
YAML.write_file(joinpath(@__DIR__,"colormap.yaml"),cm)

max_stars_repo_name: Experica/Imager

max_issues_repo_name: Experica/Imager

max_forks_repo_name: Experica/Imager


text:
#### Build full Gram matrix (ReproducingKernel_0)

function _gram!(
        A::AbstractMatrix,
        nodes::AbstractVecOfSVecs,
        kernel::ReproducingKernel_0,
    )
    n₁ = length(nodes)
    @inbounds for j in 1:n₁
        for i in 1:j
            A[i,j] = _rk(kernel, nodes[i], nodes[j])
        end
    end
    return Hermitian(A, :U)
end

function _gram(
        nodes::AbstractVecOfSVecs,
        kernel::ReproducingKernel_0,
    )
    n₁ = length(nodes)
    T  = eltype(eltype(nodes))
    _gram!(zeros(T, n₁, n₁), nodes, kernel)
end

#### Incrementally add column to Gram matrix (ReproducingKernel_0)

function _gram!(
        A::AbstractMatrix,
        new_node::SVector,
        curr_nodes::AbstractVecOfSVecs,
        kernel::ReproducingKernel_0,
    )
    n₁ = length(curr_nodes)
    @inbounds for i in 1:n₁
        A[i, n₁+1] = _rk(kernel, curr_nodes[i], new_node)
    end
    @inbounds A[n₁+1, n₁+1] = _rk(kernel, new_node, new_node)
    return Hermitian(A, :U)
end

#### Build full Gram matrix (ReproducingKernel_1)

function _gram!(
        A::AbstractMatrix,
        nodes::AbstractVecOfSVecs{n},
        d_nodes::AbstractVecOfSVecs{n},
        d_dirs::AbstractVecOfSVecs{n},
        kernel::ReproducingKernel_1,
    ) where {n}
    n₁  = length(nodes)
    n₂  = length(d_nodes)
    A11 = A
    A12 = uview(A, 1    : n₁,    n₁+1 : n₁+n₂)
    A22 = uview(A, n₁+1 : n₁+n₂, n₁+1 : n₁+n₂)

    @inbounds for j in 1:n₁
        # Top-left block (n₁ × n₁)
        for i in 1:j
            A11[i,j] = _rk(kernel, nodes[i], nodes[j])
        end
    end

    ε² = kernel.ε^2
    @inbounds for j in 1:n₂
        # Top-right block (n₁ × n₂)
        for i in 1:n₁
            A12[i,j] = _∂rk_∂e(kernel, nodes[i], d_nodes[j], d_dirs[j])
        end

        # Bottom-right block (n₂ × n₂)
        for i in 1:j-1
            A22[i,j] = _∂²rk_∂²e(kernel, d_nodes[j], d_nodes[i], d_dirs[j], d_dirs[i])
        end
        A22[j,j] = ε²
    end

    return Hermitian(A, :U)
end

function _gram(
        nodes::AbstractVecOfSVecs{n},
        d_nodes::AbstractVecOfSVecs{n},
        d_dirs::AbstractVecOfSVecs{n},
        kernel::ReproducingKernel_1,
    ) where {n}
    n₁ = length(nodes)
    n₂ = length(d_nodes)
    T  = promote_type(eltype(eltype(nodes)), eltype(eltype(d_nodes)), eltype(eltype(d_dirs)))
    _gram!(zeros(T, n₁+n₂, n₁+n₂), nodes, d_nodes, d_dirs, kernel)
end

#### Incrementally add column to Gram matrix (ReproducingKernel_1)

function _gram!(
        A::AbstractMatrix,
        new_node::SVector{n},
        curr_nodes::AbstractVecOfSVecs{n},
        curr_d_nodes::AbstractVecOfSVecs{n},
        curr_d_dirs::AbstractVecOfSVecs{n},
        kernel::ReproducingKernel_1,
    ) where {n}
    n₁ = length(curr_nodes)
    n₂ = length(curr_d_nodes)
    @assert size(A) == (n₁+1+n₂, n₁+1+n₂)

    # Top-left block (n₁+1 × n₁+1), right column (n₁+1 terms)
    @inbounds for i in 1:n₁
        A[i, n₁+1] = _rk(kernel, new_node, curr_nodes[i])
    end
    @inbounds A[n₁+1, n₁+1] = _rk(kernel, new_node, new_node)

    # Top-right block (n₁+1 × n₂), bottom row (n₂ terms)
    @inbounds for j in 1:n₂
        A[n₁+1, n₁+1+j] = _∂rk_∂e(kernel, new_node, curr_d_nodes[j], curr_d_dirs[j])
    end

    return Hermitian(A, :U)
end

function _gram!(
        A::AbstractMatrix,
        d_node::SVector{n},
        d_dir::SVector{n},
        curr_nodes::AbstractVecOfSVecs{n},
        curr_d_nodes::AbstractVecOfSVecs{n},
        curr_d_dirs::AbstractVecOfSVecs{n},
        kernel::ReproducingKernel_1,
    ) where {n}
    n₁ = length(curr_nodes)
    n₂ = length(curr_d_nodes)
    @assert size(A) == (n₁+n₂+1, n₁+n₂+1)

    # Top-right block, (n₁ × n₂+1), right column (n₁ terms)
    @inbounds for i in 1:n₁
        A[i, n₁+n₂+1] = _∂rk_∂e(kernel, curr_nodes[i], d_node, d_dir)
    end

    # Bottom-right block (n₂+1 × n₂+1), right column (n₂+1 terms)
    ε² = kernel.ε^2
    @inbounds for i in 1:n₂
        A[n₁+i, n₁+n₂+1] = _∂²rk_∂²e(kernel, d_node, curr_d_nodes[i], d_dir, curr_d_dirs[i])
    end
    @inbounds A[n₁+n₂+1, n₁+n₂+1] = ε²

    return Hermitian(A, :U)
end

#### Elastic Cholesky

Base.@kwdef struct ElasticCholesky{T, AType <: AbstractMatrix{T}} <: LinearAlgebra.Factorization{T}
    maxcols::Int
    ncols::Base.RefValue{Int} = Ref(0)
    colperms::Vector{Int}     = zeros(Int, maxcols)
    A::AType                  = zeros(T, maxcols, maxcols)
    U::Matrix{T}              = zeros(T, maxcols, maxcols)
    U⁻ᵀb::Vector{T}           = zeros(T, maxcols)
end
ElasticCholesky{T}(maxcols::Int) where {T} = ElasticCholesky{T,Matrix{T}}(; maxcols = maxcols)
ElasticCholesky(A::AbstractMatrix{T}) where {T} = ElasticCholesky{T,typeof(A)}(; maxcols = size(A,2), A = A)

Base.eltype(::ElasticCholesky{T}) where {T} = T
Base.size(C::ElasticCholesky) = (C.ncols[], C.ncols[])
Base.parent(C::ElasticCholesky) = C.A
Base.empty!(C::ElasticCholesky) = (C.ncols[] = 0; C)
Base.show(io::IO, mime::MIME"text/plain", C::ElasticCholesky{T}) where {T} = (print(io, "ElasticCholesky{T}\nU factor:\n"); show(io, mime, UpperTriangular(C.U[C.colperms[1:C.ncols[]], C.colperms[1:C.ncols[]]])))

function LinearAlgebra.ldiv!(x::AbstractVector{T}, C::ElasticCholesky{T}, b::AbstractVector{T}, ::Val{permview} = Val(false)) where {T, permview}
    @unpack U, U⁻ᵀb, colperms, ncols = C
    J = uview(colperms, 1:ncols[])
    U = UpperTriangular(uview(U, J, J))
    U⁻ᵀb = uview(U⁻ᵀb, 1:ncols[])
    if permview
        x = uview(x, J)
        b = uview(b, J)
    end
    ldiv!(U⁻ᵀb, U', b)
    ldiv!(x, U, U⁻ᵀb)
    return x
end

function Base.insert!(C::ElasticCholesky{T}, j::Int, B::AbstractMatrix{T}) where {T}
    @unpack A, colperms, ncols = C
    @inbounds colperms[ncols[] + 1] = j
    rows = uview(colperms, 1 : ncols[] + 1)
    @inbounds for i in rows
        A[i,j] = B[i,j]
    end
    return C
end

"""
    LinearAlgebra.cholesky!(C::ElasticCholesky, v::AbstractVector{T}) where {T}

Update the Cholesky factorization `C` as if the column `v` (and by symmetry, the corresponding row `vᵀ`)
were inserted into the underlying matrix `A`. Specifically, let `L` be the lower-triangular cholesky factor
of `A` such that `A = LLᵀ`, and let `v = [d; γ]` such that the new matrix `A⁺` is given by

```
A⁺ = [A  d]
     [dᵀ γ].
```

Then, the corresponding updated cholesky factor `L⁺` of `⁺` is:

```
L⁺ = [L  e]
     [eᵀ α]
```

where `e = L⁻¹d`, `α = √τ`, and `τ = γ - e⋅e > 0`. If `τ ≤ 0`, then `A⁺` is not positive definite.

See:
    https://igorkohan.github.io/NormalHermiteSplines.jl/dev/Normal-Splines-Method/#Algorithms-for-updating-Cholesky-factorization
"""
function LinearAlgebra.cholesky!(
        C::ElasticCholesky{T},
        j::Int,
        v::AbstractVector{T},
        ::Val{fill_parent},
    ) where {T, fill_parent}
    @unpack maxcols, A, U, colperms, ncols = C
    @assert length(v) == ncols[] + 1 <= maxcols

    @inbounds if ncols[] == 0
        # Initialize first entry of `A`
        colperms[1] = j
        if fill_parent
            A[j,j] = v[1]
        end
        U[j,j] = sqrt(v[1])
        ncols[] = 1
    else
        # Fill `A` with new column
        colperms[ncols[] + 1] = j
        if fill_parent
            rows = uview(colperms, 1:ncols[] + 1)
            copyto!(uview(A, rows, j), v)
        end

        # Update `U` with new column
        J = uview(colperms, 1:ncols[])
        d = uview(A, J, j)
        γ = A[j,j]
        e = uview(U, J, j)
        Uᵀ = UpperTriangular(uview(U, J, J))'
        ldiv!(e, Uᵀ, d)
        τ = γ - e⋅e
        α = √max(τ, 0) # `τ` should be positive by construction
        U[j,j] = max(α, eps(T)) # if `α < ϵ` you have bigger problems...

        # Increment column counter
        ncols[] += 1
    end

    return C
end

# Update the `j`th column of the factorization `C.U`, assuming the corresponding column `j` of `C.A` has been filled
function LinearAlgebra.cholesky!(C::ElasticCholesky{T}, j::Int) where {T}
    @unpack maxcols, A, colperms, ncols = C
    @assert ncols[] + 1 <= maxcols

    @inbounds colperms[ncols[] + 1] = j
    rows = uview(colperms, 1:ncols[] + 1)
    v = uview(A, rows, j)
    cholesky!(C, j, v, Val(false))

    return C
end

# Update columns `J` of the factorization `C.U`, assuming the corresponding columns `J` of `C.A` have been filled
function LinearAlgebra.cholesky!(C::ElasticCholesky, J = axes(C.A, 2))
    for j in J
        cholesky!(C, j)
    end
    return C
end


text:
# TODO: optimize this
function Base.string(dt::DateTime)
    y,m,d = yearmonthday(days(dt))
    h,mi,s = hour(dt),minute(dt),second(dt)
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    mm = lpad(m,2,"0")
    dd = lpad(d,2,"0")
    hh = lpad(h,2,"0")
    mii = lpad(mi,2,"0")
    ss = lpad(s,2,"0")
    ms = millisecond(dt) == 0 ? "" : string(millisecond(dt)/1000.0)[2:end]
    return "$yy-$mm-$(dd)T$hh:$mii:$ss$(ms)"
end
Base.show(io::IO,x::DateTime) = print(io,string(x))
function Base.string(dt::Date)
    y,m,d = yearmonthday(value(dt))
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    mm = lpad(m,2,"0")
    dd = lpad(d,2,"0")
    return "$yy-$mm-$dd"
end
Base.show(io::IO,x::Date) = print(io,string(x))

### Parsing
const english = (UTF8String=>Int)["january"=>1,"february"=>2,"march"=>3,"april"=>4,
                 "may"=>5,"june"=>6,"july"=>7,"august"=>8,"september"=>9,
                 "october"=>10,"november"=>11,"december"=>12]
const abbrenglish = (UTF8String=>Int)["jan"=>1,"feb"=>2,"mar"=>3,"apr"=>4,
                     "may"=>5,"jun"=>6,"jul"=>7,"aug"=>8,"sep"=>9,
                     "oct"=>10,"nov"=>11,"dec"=>12]
const MONTHTOVALUE = (UTF8String=>Dict{UTF8String,Int})["english"=>english]
const MONTHTOVALUEABBR = (UTF8String=>Dict{UTF8String,Int})["english"=>abbrenglish]

# Date/DateTime Parsing
abstract Slot{P<:AbstractTime}

immutable DelimitedSlot{P<:AbstractTime} <: Slot{P}
    i::Int
    period::Type{P}
    width::Int
    option::Int
    locale::String
end

immutable FixedWidthSlot{P<:AbstractTime} <: Slot{P}
    i::Int
    period::Type{P}
    width::Int
    option::Int
    locale::String
end

immutable DateFormat
    slots::Array{Slot,1}
    begtran # optional transition from the start of a string to the 1st slot
    trans #trans[i] == how to transition FROM slots[i] TO slots[i+1]
end

immutable DayOfWeekSlot <: AbstractTime end

duplicates(slots) = any(map(x->count(y->x.period==y.period,slots),slots) .> 1)

function DateFormat(f::String,locale::String="english")
    slots = Slot[]
    trans = {}
    begtran = match(r"^.*?(?=[ymuUdHMSsEe])",f).match
    ss = split(f,r"^.*?(?=[ymuUdHMSsEe])")
    s = split(begtran == "" ? ss[1] : ss[2],r"[^ymuUdHMSsEe]+|(?<=([ymuUdHMSsEe])(?!\1))")
    for (i,k) in enumerate(s)
        k == "" && break
        tran = i >= endof(s) ? r"$" : match(Regex("(?<=$(s[i])).*(?=$(s[i+1]))"),f).match
        slot = tran == "" || tran == r"$" ? FixedWidthSlot : DelimitedSlot
        width = length(k)
        typ = 'E' in k ? DayOfWeekSlot : 'e' in k ? DayOfWeekSlot :
              'y' in k ? Year : 'm' in k ? Month :
              'u' in k ? Month : 'U' in k ? Month :
              'd' in k ? Day : 'H' in k ? Hour :
              'M' in k ? Minute : 'S' in k ? Second : Millisecond
        option = 'E' in k ? 2 : 'e' in k ? 1 :
                 'U' in k ? 2 : 'u' in k ? 1 : 0
        push!(slots,slot(i,typ,width,option,locale))
        push!(trans,tran)
    end
    duplicates(slots) && throw(ArgumentError("Two separate periods of the same type detected"))
    return DateFormat(slots,begtran,trans)
end

const SLOTERROR = ArgumentError("Non-digit character encountered")
slotparse(slot,x) = !ismatch(r"[^0-9\s]",x) ? slot.period(x) : throw(SLOTERROR)
function slotparse(slot::Slot{Month},x)
    if slot.option == 0
        ismatch(r"[^0-9\s]",x) ? throw(SLOTERROR) : return Month(x)
    elseif slot.option == 1
        return Month(MONTHTOVALUEABBR[slot.locale][lowercase(x)])
    else
        return Month(MONTHTOVALUE[slot.locale][lowercase(x)])
    end
end
slotparse(slot::Slot{Millisecond},x) = !ismatch(r"[^0-9\s]",x) ? slot.period(parsefloat("."*x)*1000.0) : throw(SLOTERROR)
slotparse(slot::Slot{DayOfWeekSlot},x) = nothing

function getslot(x,slot::DelimitedSlot,df,cursor)
    endind = first(search(x,df.trans[slot.i],cursor+1))
    if endind == 0 # we didn't find the next delimiter
        s = x[cursor:end]
        return (endof(x)+1, isdigit(s) ? slotparse(slot,s) : default(slot.period))
    end
    return endind+1, slotparse(slot,x[cursor:(endind-1)])
end
getslot(x,slot,df,cursor) = (cursor+slot.width, slotparse(slot,x[cursor:(cursor+slot.width-1)]))

function parse(x::String,df::DateFormat)
    x = strip(replace(x, r"#.*$", ""))
    x = replace(x,df.begtran,"")
    isempty(x) && throw(ArgumentError("Cannot parse empty format string"))
    (typeof(df.slots[1]) <: DelimitedSlot && first(search(x,df.trans[1])) == 0) && throw(ArgumentError("Delimiter mismatch. Couldn't find first delimiter, \"$(df.trans[1])\", in date string"))
    periods = Period[]
    cursor = 1
    for slot in df.slots
        cursor, pe = getslot(x,slot,df,cursor)
        pe != nothing && push!(periods,pe)
        cursor > endof(x) && break
    end
    return sort!(periods,rev=true,lt=periodisless)
end

slotformat(slot::Slot{Year},dt) = lpad(string(value(slot.period(dt))),slot.width,"0")[(end-slot.width+1):end]
slotformat(slot,dt) = lpad(string(value(slot.period(dt))),slot.width,"0")
function slotformat(slot::Slot{Month},dt)
    if slot.option == 0
        return lpad(month(dt),slot.width,"0")
    elseif slot.option == 1
        return VALUETOMONTHABBR[slot.locale][month(dt)]
    else
        return VALUETOMONTH[slot.locale][month(dt)]
    end
end
function slotformat(slot::Slot{DayOfWeekSlot},dt)
    if slot.option == 1
        return VALUETODAYOFWEEKABBR[slot.locale][dayofweek(dt)]
    else # == 2
        return VALUETODAYOFWEEK[slot.locale][dayofweek(dt)]
    end
end
slotformat(slot::Slot{Millisecond},dt) = rpad(string(millisecond(dt)/1000.0)[3:end], slot.width, "0")

function format(dt::TimeType,df::DateFormat)
    f = ""
    for slot in df.slots
        f *= slotformat(slot,dt)
        f *= typeof(df.trans[slot.i]) <: Regex ? "" : df.trans[slot.i]
    end
    return f
end

# UI
const ISODateTimeFormat = DateFormat("yyyy-mm-ddTHH:MM:SS.s")
const ISODateFormat = DateFormat("yyyy-mm-dd")
const RFC1123Format = DateFormat("e, dd u yyyy HH:MM:SS")

DateTime(dt::String,format::String;locale::String="english") = DateTime(dt,DateFormat(format,locale))
DateTime(dt::String,df::DateFormat=ISODateTimeFormat) = DateTime(parse(dt,df)...)

Date(dt::String,format::String;locale::String="english") = Date(dt,DateFormat(format,locale))
Date(dt::String,df::DateFormat=ISODateFormat) = Date(parse(dt,df)...)

format(dt::TimeType,f::String;locale::String="english") = format(dt,DateFormat(f,locale))

# vectorized
DateTime{T<:String}(y::AbstractArray{T},format::String;locale::String="english") = DateTime(y,DateFormat(format,locale))
function DateTime{T<:String}(y::AbstractArray{T},df::DateFormat=ISODateTimeFormat)
    return reshape(DateTime[DateTime(parse(y[i],df)...) for i in 1:length(y)], size(y))
end
Date{T<:String}(y::AbstractArray{T},format::String;locale::String="english") = Date(y,DateFormat(format,locale))
function Date{T<:String}(y::AbstractArray{T},df::DateFormat=ISODateFormat)
    return reshape(Date[Date(parse(y[i],df)...) for i in 1:length(y)], size(y))
end

format{T<:TimeType}(y::AbstractArray{T},format::String;locale::String="english") = Dates.format(y,DateFormat(format,locale))
function format(y::AbstractArray{Date},df::DateFormat=ISODateFormat)
    return reshape([Dates.format(y[i],df) for i in 1:length(y)], size(y))
end
function format(y::AbstractArray{DateTime},df::DateFormat=ISODateTimeFormat)
    return reshape([Dates.format(y[i],df) for i in 1:length(y)], size(y))
end


text:
using Pkg
Pkg.activate(normpath(joinpath(@__DIR__, "..")))

using ATools
using NReco

using ArgParse
using Configurations
using DataFrames
using Glob
using HDF5
using Unitful

import Unitful: ns, mm

function define_paths(config::NReco.CalConfig)
    in_path  = joinpath(config.input_dir, config.conf_dir)
    out_path = joinpath(config.plot_dir , config.conf_dir)
    if !isdir(out_path)
        mkpath(out_path)
    end
    return in_path, out_path
end

function calibrate_lors()
    s = ArgParseSettings()

    @add_arg_table s begin
        "--conf", "-c"
            help     = "Calibration configuration"
			arg_type = String
        "--trueout", "-t"
            help     = "Output true position LORs?"
            arg_type = Bool
            default  = false
    end
    confArgs     = parse_args(s)

    conf         = from_toml(NReco.CalConfig, confArgs["conf"])

    (path_in ,
     path_out)   = define_paths(conf)

    (nsim, rmin,
     rmax, ndf ) = read_evtpar(glob("evt*.h5", path_in))

    cal_func     = NReco.calibration_function(conf.cal_func, rmin, rmax)

    ndfq         = NReco.filter_energies(ndf, conf.qmin, conf.qmax)

    NReco.calculate_interaction_radius!(ndfq, cal_func, conf.cal_func.cal_std, rmax)

    # We need/want units here. More generalised use of units should be implemented.
    ## Get the true interaction radius (Why not saved?)
    transform!(ndfq, [:xt1, :yt1] => ByRow(rxy) => :rt1, [:xt2, :yt2] => ByRow(rxy) => :rt2)
    units_ndfq = ATools.set_units(ndfq)

    # Calculate the x, y, z of the interaction.
    xint1, yint1, zint1 = ATools.radial_correction(units_ndfq.xr1 / mm, units_ndfq.yr1 / mm,
                                                   units_ndfq.zr1 / mm, units_ndfq.r1x / mm)
    xint2, yint2, zint2 = ATools.radial_correction(units_ndfq.xr2 / mm, units_ndfq.yr2 / mm,
                                                   units_ndfq.zr2 / mm, units_ndfq.r2x / mm)
    # Calculate the interaction time.
    t1 = ATools.interaction_time(units_ndfq, :r1x, :ta1, rmax, conf.cal_func.nLXe)
    t2 = ATools.interaction_time(units_ndfq, :r2x, :ta2, rmax, conf.cal_func.nLXe)

    # Calculate LORs.
    dt   = uconvert.(ns, t2 - t1) ./ ns
    mLor = ATools.MlemLor.(dt, xint1, yint1, zint1, xint2, yint2, zint2,
        units_ndfq[!, :q1], units_ndfq[!, :q2], units_ndfq[!, :E1], units_ndfq[!, :E2])
    # Will want a dataset/table with metadata too, to be decided.
    mlor_filename = joinpath(path_out, conf.conf_dir[1:end-1] * "_mlor.h5")
    ATools.write_lors_hdf5(mlor_filename, mLor)
    if confArgs["trueout"]
        flight1   = ATools.time_of_flight(units_ndfq, [:xs, :ys, :zs], [:xt1, :yt1, :zt1])
        flight2   = ATools.time_of_flight(units_ndfq, [:xs, :ys, :zs], [:xt2, :yt2, :zt2])
        dt        = uconvert.(ns, flight2 - flight1) ./ ns
        true_mLor = ATools.MlemLor.(dt, ndfq.xt1, ndfq.yt1, ndfq.zt1,
            ndfq.xt2, ndfq.yt2, ndfq.zt2,
            units_ndfq[!, :q1], units_ndfq[!, :q2], units_ndfq[!, :E1], units_ndfq[!, :E2])
        mlor_filename = joinpath(path_out, conf.conf_dir[1:end-1] * "_Truemlor.h5")
        ATools.write_lors_hdf5(mlor_filename, true_mLor, "true_info")
    end
end
calibrate_lors()

max_stars_repo_name: andLaing/NReco

max_issues_repo_name: andLaing/NReco

max_forks_repo_name: andLaing/NReco


text:
using DataStructures: SortedDict, OrderedDict

abstract type AbstractBags end

Base.getindex(b::AbstractBags, i) = b.bags[i]
Base.setindex!(b::AbstractBags, v, i) = setindex!(b.bags, v, i)
Base.iterate(b::AbstractBags) = iterate(b.bags)
Base.iterate(b::AbstractBags, s) = iterate(b.bags, s)
Base.length(b::AbstractBags) = length(b.bags)

# one instance belongs to only one bag
# sorted from left to right
struct AlignedBags <: AbstractBags
    bags::Vector{UnitRange{Int}}
end

AlignedBags() = AlignedBags(Vector{UnitRange{Int}}())
AlignedBags(ks::UnitRange{Int}...) = AlignedBags(collect(ks))
function AlignedBags(k::Vector{T}) where {T<:Integer}
    b = AlignedBags()
    !isempty(k) || return b
    a, v = 1, k[1]
    s = Set{T}(v)
    for (i, x) in enumerate(k[2:end])
        if x != v
            !(x in s) || error("Scattered bags")
            push!(b.bags, a:i)
            v = x
            a = i+1
            push!(s, x)
        end
    end
    push!(b.bags, a:length(k))
    b
end

# one instance may belong to more bags
struct ScatteredBags <: AbstractBags
    bags::Vector{Vector{Int}}
end

ScatteredBags() = ScatteredBags(Vector{Vector{Int}}())
function ScatteredBags(k::Vector{T}) where {T<:Integer}
    !isempty(k) || return ScatteredBags()
    d = SortedDict{T, Vector{Int}}()
    for (i, x) in enumerate(k)
        if !(x in keys(d))
            d[x] = Int[]
        end
        push!(d[x], i)
    end
    ScatteredBags(collect(values(d)))
end

Zygote.@nograd function length2bags(ls::Vector{Int})
    ls = vcat([0], cumsum(ls))
    bags = map(i -> i[1]+1:i[2],zip(ls[1:end-1],ls[2:end]))
    bags = map(b -> isempty(b) ? (0:-1) : b,bags)
    AlignedBags(bags)
end

# """
# 		function bags(k::Vector)

# 		creates AlignedBags if it is possible to do so, otherwise creates ScatteredBags instance
# """
bags(b::AbstractBags) = b
bags(b::Vector{UnitRange{Int}}) = AlignedBags(b)
function bags(k::Vector{T}) where {T<:Integer}
    try
        return AlignedBags(k)
    catch ErrorException
        return ScatteredBags(k)
    end
end

"""
    function remapbag(b::Bags,idcs::Vector{Int})

    bags corresponding to indices with collected indices
"""
function remapbag(b::AlignedBags, idcs::VecOrRange)
    rb = AlignedBags(Vector{UnitRange{Int}}(undef, length(idcs)))
    offset = 1
    for (i,j) in enumerate(idcs)
        rb[i] = (b[j] == 0:-1) ? b[j] : b[j] .- b[j].start .+ offset
        offset += length(b[j])
    end
    rb, Array{Int}(vcat([collect(b[i]) for i in idcs]...))
end

function remapbag(b::ScatteredBags, idcs::VecOrRange)
    rb = ScatteredBags(Vector{Vector{Int}}(undef, length(idcs)))
    m = OrderedDict{Int, Int}((v => i for (i, v) in enumerate(unique(vcat(b.bags[idcs]...)))))
    for (i,j) in enumerate(idcs)
        rb[i] = [m[v] for v in b[j]]
    end
    rb, collect(keys(m))
end

adjustbags(bags::AlignedBags, mask::T) where {T<:Union{Vector{Bool}, BitArray{1}}} = length2bags(map(b -> sum(@view mask[b]), bags))

Base.vcat(b1::AbstractBags, b2::AbstractBags) = _catbags([b1, b2])
Base.vcat(bs::AbstractBags...) = _catbags(collect(bs))

function _catbags(bs::Vector{AlignedBags})
    nbs = AlignedBags()
    offset = 0
    for b in bs
        !isempty(b) || continue
        append!(nbs.bags, [bb .+ (isempty(bb) ? 0 : offset) for bb in b])
        offset += max(0, mapreduce(i -> isempty(i) ? 0 : maximum(i), max, b))
    end
    mask = length.(nbs.bags) .== 0
    if sum(mask) > 0
        nbs[mask] = fill(0:-1, sum(mask))
    end
    nbs
end

function _catbags(bs::Vector{ScatteredBags})
    nbs = ScatteredBags()
    offset = 0
    for b in bs
        !isempty(b) || continue
        append!(nbs.bags, [bb .+ offset for bb in b])
        offset += max(0, mapreduce(i -> isempty(i) ? 0 : maximum(i), max, b))
    end
    nbs
end


text:
"Fourier transform of 2D Wavelet"
function Ψ(w; ξ0, L, c)
    ϕ = atan(w[2], w[1])
    t = (norm(w) - ξ0) / ξ0
    if (-π / 2 < ϕ < π / 2) && (-1 < t < 1)
        c * exp(-abs2(t) / (1 - abs2(t))) * cos(ϕ)^(L - 1)
    else
        0.
    end
end

rotate(w, θ) = [cos(θ) -sin(θ); sin(θ) cos(θ)] * w
scale(w, j) = 2^j * w

"Fourier transform of 2D Wavelet with scaling and rotation"
Ψ(w, j, θ; ξ0, L, c) = Ψ(scale(rotate(w, θ), j); ξ0, L, c)

"Square matrix of wavelet in frequency domain"
wavelet_matrix(ws, j=0, θ=0.; ξ0, L, c) =
    [Ψ(SVector(w1, w2), j, θ; ξ0, L, c) for w1 in ws, w2 in ws]


text:
# extend `nothing` to passthrough `first`,`last`, and allied
Base.iterate(x::Nothing) = nothing
Base.firstindex(x::Nothing) = nothing
Base.lastindex(x::Nothing) = nothing
Base.getindex(collection::Nothing, key::Any) = nothing
Base.getindex(collection::Any, key::Nothing) = nothing
Base.first(x::Nothing) = nothing
Base.last(x::Nothing) = nothing


text:
#module CBuffers
#export AbstractCBuffer, CBuffer, CBuffer2D, open, close, isgood, write!, next!
import Base:open, close
using Base.Filesystem: touch, abspath

#TODO FIX!
#const lib = abspath("./libbuffer")
const lib = "$(homedir())/gits/JSort/src/buffer/libbuffer.so"
abstract type Buffer end
struct Buffer2D <: Buffer end

function makebuffer(path::AbstractString; mode::Symbol=:read)
    if mode == :read
        ptr = ccall((:makeBufferR, lib), Ptr{Buffer}, (Cstring, ), path)
    elseif mode == :write
        ptr = ccall((:makeBufferW, lib), Ptr{Buffer}, (Cstring, ), path)
    else
        error("Only :read and :write supported for CBuffer")
    end
end

function close(buffer::Ptr{Buffer})
    ccall((:close, lib), Cvoid, (Ptr{Buffer},), buffer)
end

function next!(buffer::Ptr{Buffer})::UInt32
    val = ccall((:next, lib), UInt32, (Ptr{Buffer},), buffer)
end

function isgood(buffer::Ptr{Buffer})
    good = ccall((:good, lib), Bool, (Ptr{Buffer},), buffer)
end

function isgood(buffer::Ptr{Buffer2D})
    good = ccall((:good, lib), Bool, (Ptr{Buffer2D},), buffer)
end

function makebuffer2d(path::AbstractString, mode::Symbol=:write)
    ptr = ccall((:makeBufferW2D, lib), Ptr{Buffer2D}, (Cstring, ), path)
end

function write!(buffer::Ptr{Buffer}, x::UInt32)
    ccall((:fill, lib), Cvoid, (Ptr{Buffer}, UInt32), buffer, x)
end
write!(buffer::Ptr{Buffer}, x) = write!(buffer, reinterpret(UInt32, x))

function write!(buffer::Ptr{Buffer2D}, mat::AbstractMatrix)
    for row in 1:size(mat, 1)
        #println("Writing [$(mat[row, 1]), $(mat[row, 2])]")
        ccall((:fill2D, lib), Cvoid, (Ptr{Buffer2D}, UInt32, UInt32),
              buffer, mat[row, 1], mat[row, 2])
    end
end

function write!(buffer::Ptr{Buffer2D}, x::UInt32, y::UInt32)
    ccall((:fill2D, lib), Cvoid, (Ptr{Buffer2D}, UInt32, UInt32),
          buffer, x, y)
end
write!(buffer::Ptr{Buffer2D}, x, y) = write!(buffer, reinterpret(UInt32, x),
                                             reinterpret(UInt32, y))

function close(buffer::Ptr{Buffer2D})
    ccall((:close2D, lib), Cvoid, (Ptr{Buffer2D},), buffer)
end

abstract type AbstractCBuffer end
struct CBuffer <: AbstractCBuffer
    ptr::Ptr{Buffer}
    mode::Symbol
end
struct CBuffer2D <: AbstractCBuffer
    ptr::Ptr{Buffer2D}
    mode::Symbol
end
CBuffer(path::AbstractString, mode::Symbol=:read) = open(path, CBuffer, mode=mode)
CBuffer2D(path::AbstractString) = open(path, CBuffer2D)

function open(path::AbstractString, T::Type{AbstractCBuffer}; mode::Symbol=:read)
    open(path, T, mode=Val{mode})
end

function open(path::AbstractString, T::Type{CBuffer}; mode=:read)
    if mode != :read && mode != :write
        error("Only reading xor writing supported for CBuffer")
    end
    buffer_ptr = makebuffer(path, mode=mode)
    return CBuffer(buffer_ptr, mode)
end

function open(path::AbstractString, T::Type{CBuffer2D}; mode::Symbol=:write)
    if mode != :write
        error("Only writing supported for CBuffer2D")
    end
    touch(path)  # The file is not created by the C-code
    buffer_ptr = makebuffer2d(path)
    return CBuffer2D(buffer_ptr, :write)
end

function close(buffer::AbstractCBuffer)
    close(buffer.ptr)
end

function next!(buffer::AbstractCBuffer)
    next!(buffer.ptr)
end

write!(buffer::CBuffer2D, m::AbstractMatrix) = write!(buffer.ptr, m)
write!(buffer::CBuffer2D, x, y) = write!(buffer.ptr, x, y)
write!(buffer::CBuffer, x) = write!(buffer.ptr, x)
isgood(buffer::AbstractCBuffer) = isgood(buffer.ptr)
isgood(buffer::CBuffer) = isgood(buffer.ptr)

function main()
    buffer = makebuffer("duck")
    fact = read("duck") |> x -> reinterpret(UInt32, x)
    i = 1
    test = similar(fact)
    @allocated while isgood(buffer)
        test[i] = next!(buffer)
        if test[i] != fact[i]
            println(i)
        end
        i += 1
    end
    close(buffer)
end

function main2()
    buffer = makebuffer2d("ede.bin")
    mat = UInt32[
        10 2
        11 3
        4 5
        12 6]
    write!(buffer, mat)
    close(buffer)
end

function main3()
    buffer = makebuffer2d("ede.bin")
    @time mat = [r+c for r in 1:10_00_000, c in 0:1]
    @time write!(buffer, mat)
    @time close(buffer)
end

#end

max_stars_repo_name: Caronthir/JSort

max_issues_repo_name: Caronthir/JSort

max_forks_repo_name: Caronthir/JSort


text:
# Bertrand Nash example
# Generates random equilibrium data
function bestResponsePrice(thetas, intercept, prices, pmax, FIRM_ME)
	if FIRM_ME == :1
		me = 1; other = 2
	else
		me = 2; other = 1
	end
	p = (intercept + thetas[other] * prices[other]) / 2 / abs(thetas[me])
	if p <= 0
		return 0
	elseif p >= pmax
		return pmax
	else
		return p
	end
end

#solve for equilibrium prices.  True model is linear with random effects
function solveNashPrices(thetas1, thetas2, intercept1, intercept2, pmax; MAX_ITERS = 20, TOL=1e-6, trace=false)
	prices = fill(convert(Float64, pmax), 2)
	for i = 1:MAX_ITERS
		prices_old = copy(prices)
		prices[1] = bestResponsePrice(thetas1, intercept1, prices, pmax, :1)
		prices[2] = bestResponsePrice(thetas2, intercept2, prices, pmax, :2)
		dist = norm(prices - prices_old)
		if trace
			println( prices', dist )
		end
		if dist < TOL
			return prices
		end
	end
	error("Max Iters reached")
end

solveNashPricesFit(thetas1, thetas2, intercept1, intercept2, GDP, pmax) =
	solveNashPrices(thetas1[1:2], thetas2[1:2], thetas1[3] * GDP + intercept1, thetas2[3] * GDP + intercept2, pmax)

function genErrors()
	eps1, eps2, eps3 = rand(Normal(5, .2), 3)
	eps3 += eps1 + eps2
	return eps1, eps2, eps3
end

max_stars_repo_name: nlperic/InverseVIs

max_issues_repo_name: nlperic/InverseVIs