text:
using NLopt
using Base.Test

count = 0 # keep track of # function evaluations

function myfunc(x::Vector, grad::Vector)
    if length(grad) > 0
        grad[1] = 0
        grad[2] = 0.5/sqrt(x[2])
    end

    global count
    count::Int += 1
    println("f_$count($x)")

    sqrt(x[2])
end

function myconstraint(x::Vector, grad::Vector, a, b)
    if length(grad) > 0
        grad[1] = 3a * (a*x[1] + b)^2
        grad[2] = -1
    end
    (a*x[1] + b)^3 - x[2]
end

opt = Opt(:LD_MMA, 2)
lower_bounds!(opt, [-Inf, 0.])
xtol_rel!(opt,1e-4)

min_objective!(opt, myfunc)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,2,0), 1e-8)
inequality_constraint!(opt, (x,g) -> myconstraint(x,g,-1,1), 1e-8)

(minf,minx,ret) = optimize(opt, [1.234, 5.678])
println("got $minf at $minx after $count iterations (returned $ret)")

@test_approx_eq_eps minx[1] 1/3 1e-5
@test_approx_eq_eps minx[2] 8/27 1e-5
@test_approx_eq_eps minf sqrt(8/27) 1e-5
@test ret == :XTOL_REACHED


text:
# Copyright 2017 Stephen L. Smith and Frank Imeson
# Modified copyright 2020 A. Kudryavtsev and M. Khachay
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

#####################################################
#########  GTSP Utilities ###########################

""" tour type that stores the order array and the length of the tour
"""
mutable struct Tour
    tour::Array{Int64, 1}
    cost::Int64
end

""" return the vertex before tour[i] on tour """
@inline function prev_tour(tour, i)
    i != 1 && return tour[i - 1]
    return tour[length(tour)]
end

#####################################################
#########  PCGTSP Utilities ###########################

mutable struct Constraints
    ascendants::Array{Int64, 1}
    descendants::Array{Int64, 1}
end

function get_relative_ordering(set_idx::Int64, set_orderings::Array{Int64, 2})
    constraints = Constraints(Int64[], Int64[])

    # Ascendants.
    path = Set{Int64}()
    push!(path, set_idx)
    while !isempty(path)
        x = pop!(path)
        for (idx, ordering) in enumerate(set_orderings[x, :])
            if ordering == 1 && !in(idx, constraints.ascendants)
                push!(path, idx)
                push!(constraints.ascendants, idx)
            end
        end
    end
    # unique!(constraints.ascendants)

    # Descendants.
    path = Set{Int64}()
    push!(path, set_idx)
    while !isempty(path)
        x = pop!(path)
        for (idx, ordering) in enumerate(set_orderings[x, :])
            if ordering == -1 && !in(idx, constraints.descendants)
                push!(path, idx)
                push!(constraints.descendants, idx)
            end
        end
    end
    # unique!(constraints.descendants)

    return constraints
end

function calc_order_constraints(sets::Array{Any, 1}, set_orderings::Array{Int64, 2})
    sets_length = length(sets)
    order_constraints = Array{Constraints, 1}(undef, sets_length)
    Threads.@threads for set_idx in 1:sets_length
        order_constraints[set_idx] = get_relative_ordering(set_idx, set_orderings)
    end

    return order_constraints
end

@inline function calc_bounds(
    tour::Array{Int64, 1},
    set_idx::Int64,
    order_constraints::Array{Constraints, 1},
    member::Array{Int64, 1},
)
    min_insert_idx = 1
    max_insert_idx = length(tour)
    constraints = order_constraints[set_idx]
    for i in 1:length(tour)
        if in(member[tour[i]], constraints.ascendants)
            min_insert_idx = i + 1
        elseif in(member[tour[i]], constraints.descendants) &&
               max_insert_idx == length(tour)
            max_insert_idx = i
        end
    end

    if max_insert_idx < min_insert_idx
        max_insert_idx = min_insert_idx
    end

    return min_insert_idx, max_insert_idx
end

######################################################
#############  Randomizing tour before insertions ####

""" some insertions break tie by taking first minimizer -- this
randomization helps avoid getting stuck choosing same minimizer """
function pivot_tour!(RNGs::Array{MersenneTwister, 1}, tour::Array{Int64, 1})
    pivot = rand(RNGs[Threads.threadid()], 1:length(tour))
    tour = [tour[pivot:end]; tour[1:(pivot - 1)]]
end

function randomize_sets!(
    RNGs::Array{MersenneTwister, 1},
    sets::Array{Any, 1},
    sets_to_insert::Array{Int64, 1},
)
    for i in sets_to_insert
        shuffle!(RNGs[Threads.threadid()], sets[i])
    end
end

function findmember(num_vertices::Int64, sets::Array{Any, 1})
    """  create an array containing the set number for each vertex """
    member = zeros(Int64, num_vertices)
    num_verts = 0
    for i in 1:length(sets)
        set = sets[i]
        num_verts += length(set)
        for vertex in set
            if member[vertex] != 0
                error("vertex ", vertex, " belongs to more than one set")
            else
                member[vertex] = i
            end
        end
    end
    return member
end

struct Distsv
    set_vert::Array{Int64, 2}
    vert_set::Array{Int64, 2}
    min_sv::Array{Int64, 2}
end

function set_vertex_dist(dist::Array{Int64, 2}, num_sets::Int, member::Array{Int64, 1})
    """
 Computes the minimum distance between each set and each vertex
 Also compute the minimum distance from a set to a vertex, ignoring direction
 This is used in insertion to choose the next set.
 """
    numv = size(dist, 1)
    dist_set_vert = typemax(Int64) * ones(Int64, num_sets, numv)
    mindist = typemax(Int64) * ones(Int64, num_sets, numv)
    dist_vert_set = typemax(Int64) * ones(Int64, numv, num_sets)

    for i in 1:numv
        for j in 1:numv
            set = member[j]
            if dist[j, i] < dist_set_vert[set, i]
                dist_set_vert[set, i] = dist[j, i]
            end
            if dist[j, i] < mindist[set, i]  # dist from set containing j to vertex i
                mindist[set, i] = dist[j, i]
            end
            set = member[i]
            if dist[j, i] < dist_vert_set[j, set]  # dist from j to set containing i
                dist_vert_set[j, set] = dist[j, i]
            end
            if dist[j, i] < mindist[set, j] # record as distance from set containing i to j
                mindist[set, j] = dist[j, i]
            end
        end
    end
    return Distsv(dist_set_vert, dist_vert_set, mindist)
end

function set_vertex_distance(dist::Array{Int64, 2}, sets::Array{Any, 1})
    """
 Computes the minimum distance between each set and each vertex
 """
    numv = size(dist, 1)
    nums = length(sets)
    dist_set_vert = typemax(Int64) * ones(Int64, nums, numv)
    # dist_vert_set = typemax(Int64) * ones(Int64, numv, nums)
    for i in 1:nums
        for j in 1:numv
            for k in sets[i]
                newdist = min(dist[k, j], dist[j, k])
                dist_set_vert[i, j] > newdist && (dist_set_vert[i, j] = newdist)
            end
        end
    end
    return dist_set_vert
end

""" Find the set with the smallest number of vertices """
function min_set(sets::Array{Any, 1})
    min_size = length(sets[1])
    min_index = 1
    for i in 2:length(sets)
        set_size = length(sets[i])
        if set_size < min_size
            min_size = set_size
            min_index = i
        end
    end
    return min_index
end

############################################################
############ Trial acceptance criteria #####################

"""
decide whether or not to accept a trial based on simulated annealing criteria
"""
function accepttrial(
    RNGs::Array{MersenneTwister, 1},
    trial_cost::Int64,
    current_cost::Int64,
    temperature::Float64,
)
    if trial_cost <= current_cost
        accept_prob = 2.0
    else
        accept_prob = exp((current_cost - trial_cost) / temperature)
    end
    return (rand(RNGs[Threads.threadid()]) < accept_prob ? true : false)
end

"""
decide whether or not to accept a trial based on simple probability
"""
function accepttrial_noparam(
    RNGs::Array{MersenneTwister, 1},
    trial_cost::Int64,
    current_cost::Int64,
    prob_accept::Float64,
)
    if trial_cost <= current_cost
        return true
    end
    return (rand(RNGs[Threads.threadid()]) < prob_accept ? true : false)
end

###################################################
################ Tour checks ######################

"""  Compute the length of a tour  """
@inline function tour_cost(tour::Array{Int64, 1}, dist::Array{Int64, 2})
    tour_length = max(dist[tour[end], tour[1]], 0)
    for i in 1:(length(tour) - 1)
        tour_length += dist[tour[i], tour[i + 1]]
    end
    return tour_length
end

"""
Checks if a tour is feasible in that it visits each set exactly once.
"""
function tour_feasibility(
    tour::Array{Int64, 1},
    membership::Array{Int64, 1},
    num_sets::Int64,
)
    length(tour) != num_sets && return false

    set_test = falses(num_sets)
    for v in tour
        set_v = membership[v]
        if set_test[set_v]
            return false  # a set is visited twice in the tour
        end
        set_test[set_v] = true
    end
    for visited_set in set_test
        !visited_set && return false
    end
    return true
end

#####################################################
#############  Incremental Shuffle ##################

@inline function incremental_shuffle!(
    RNGs::Array{MersenneTwister, 1},
    a::AbstractVector,
    i::Int,
)
    j = i + floor(Int, rand(RNGs[Threads.threadid()]) * (length(a) + 1 - i))
    a[j], a[i] = a[i], a[j]
    return a[i]
end

""" rand_select for randomize over all minimizers """
@inline function rand_select(RNGs::Array{MersenneTwister, 1}, a::Array{Int64, 1}, val::Int)
    inds = Int[]
    for i in 1:length(a)
        a[i] == val && (push!(inds, i))
    end
    return rand(RNGs[Threads.threadid()], inds)
end

ext: jl

max_stars_repo_name: AndreiKud/PCGLNS

max_issues_repo_name: AndreiKud/PCGLNS

max_forks_repo_name: AndreiKud/PCGLNS


text:
"""
	generate_base_expressions(model::Model)
Generate fast base methods using Symbolics symbolic computing tools.
"""
function generate_base_expressions(model::Model;
	M_analytical = true,
	C_analytical = true,
	mapping = (a, x) -> I,
	nv = model.nq)

	nq = model.nq
	nu = model.nu
	nw = model.nw
	nc = model.nc

	# Declare variables
	@variables q[1:nq]
	@variables q̇[1:nv]

	# Lagrangian
	L = lagrangian(model, q, q̇)
	L = Symbolics.simplify.(L)

	if !(M_analytical && C_analytical)
		ddL = Symbolics.hessian(L, [q; q̇], simplify=true)
	end

	# Mass Matrix
	if M_analytical
		M = M_func(model, q)
		M = reshape(M, (nv, nv))
		M = Symbolics.simplify.(M)
	else
		# ddL = Symbolics.hessian(L, [q; q̇], simplify=true)
		M = ddL[nq .+ (1:nv), nq .+ (1:nv)]
	end

	# Bias terms
	if C_analytical
		C = C_func(model, q, q̇)
		C = simplify.(C)
	else
		dLq = mapping(model, q)' * Symbolics.gradient(L, q, simplify=true) # including mapping for orientation (e.g., attitude Jacobian)
		# dLq̇ = Symbolics.gradient(L, q̇, simplify=true)
		# ddL = Symbolics.hessian(L, [q; q̇], simplify=true)
		ddLq̇q = ddL[nq .+ (1:nv), 1:nq] * mapping(model, q)

		# Coriolis and Centrifugal forces Jacobians
		C = ddLq̇q * q̇ - dLq
		C = Symbolics.simplify.(C)
	end

	# Control input Jacobian
	B = B_func(model, q)
	B = reshape(B, (nu, nv))
	B = Symbolics.simplify.(B)

	# Disturbance input Jacobian
	A = A_func(model, q)
	A = reshape(A, (nw, nv))
	A = Symbolics.simplify.(A)

	# Kinematics
	k = kinematics(model, q)
	k = simplify.(k)

	# Build function
	expr = Dict{Symbol, Expr}()
	expr[:L]    = build_function([L], q, q̇, checkbounds=true)[1] # need to transpose to get a line vector
	expr[:M]    = build_function(M, q, checkbounds=true)[1]
	expr[:B]    = build_function(B, q, checkbounds=true)[1]
	expr[:A]    = build_function(A, q, checkbounds=true)[1]
	expr[:C]    = build_function(C, q, q̇, checkbounds=true)[1]
	expr[:k]    = build_function(k, q, checkbounds=true)[1]

	return expr
end

"""
	instantiate_base!(model,
		path::AbstractString="model/base.jld2")
Loads the base expressions from the `path`, evaluates them to generate functions,
stores them into the model.
"""
function instantiate_base!(model::Model, path::AbstractString="model/base.jld2")
	instantiate_base!(model.base, load_expressions(path))
	return nothing
end

"""
	instantiate_base!(model,
		path::AbstractString="model/base.jld2")
Evaluates the base expressions to generate functions, stores them into the model.
"""
function instantiate_base!(fct::BaseMethods, expr::Dict{Symbol,Expr})
	fct.L    = eval(expr[:L])
	fct.M    = eval(expr[:M])
	fct.B    = eval(expr[:B])
	fct.A    = eval(expr[:A])
	fct.C    = eval(expr[:C])
	fct.k    = eval(expr[:k])
	return nothing
end

"""
	generate_dynamics_expressions(model::Model)
Generate fast dynamics methods using Symbolics symbolic computing tools.
"""
function generate_dynamics_expressions(model::Model; derivs = false, nv = model.nq)
	nq = model.nq
	nu = model.nu
	nw = model.nw
	nc = model.nc
	# ncf = nc * dim(env)

	# Declare variables
	@variables q0[1:nq]
	@variables q1[1:nq]
	@variables u1[1:nu]
	@variables w1[1:nw]
	@variables Λ1[1:nv]
	@variables q2[1:nq]
	@variables h[1:1]

	# Expressions
	expr = Dict{Symbol, Expr}()

	# Dynamics
	d = dynamics(model, h, q0, q1, u1, w1, Λ1, q2)
	d = Symbolics.simplify.(d)

	# Functions
	expr[:d] = build_function(d, h, q0, q1, u1, w1, Λ1, q2, checkbounds=true)[1]

	if derivs
		∂q2 = Symbolics.jacobian(d, q2, simplify = true)
		dΛ1 = Symbolics.jacobian(d, Λ1, simplify = true)
		dθ =  Symbolics.jacobian(d, [q0; q1; u1; w1; h])#, simplify = true)

		expr[:∂q2] = build_function(∂q2, h, q0, q1, u1, w1, Λ1, q2, checkbounds=true)[1]
		expr[:dΛ1] = build_function(dΛ1, h, q0, q1, u1, w1, Λ1, q2, checkbounds=true)[1]
		expr[:dθ]  = build_function(dθ,  h, q0, q1, u1, w1, Λ1, q2, checkbounds=true)[1]
	end

	return expr
end

"""
	instantiate_dynamics!(model,
		path::AbstractString="dynamics.jld2")
Loads the dynamics expressoins from the `path`, evaluates them to generate functions,
stores them into the model.
"""
function instantiate_dynamics!(model::Model, path::AbstractString="model/dynamics.jld2";
	derivs = false)
	expr = load_expressions(path)
	instantiate_dynamics!(model.dyn, expr, derivs = derivs)
	return nothing
end

"""
	instantiate_dynamics!(model,
		path::AbstractString="dynamics.jld2")
Evaluates the dynamics expressions to generate functions, stores them into the dedicated struture.
"""
function instantiate_dynamics!(fct::DynamicsMethods, expr::Dict{Symbol,Expr};
	derivs = false)
	fct.d = eval(expr[:d])

	if derivs
		fct.dθ = eval(expr[:dθ])
		fct.dΛ1 = eval(expr[:dΛ1])
		fct.∂q2 = eval(expr[:∂q2])
	end

	return nothing
end


text:
# Autogenerated wrapper script for pandoc_jll for aarch64-linux-gnu
export pandoc

JLLWrappers.@generate_wrapper_header("pandoc")
JLLWrappers.@declare_executable_product(pandoc)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        pandoc,
        "bin/pandoc",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
mutable struct Simulation{T <: SimulationModel}
    status::BUILD_STATUS
    problem::Union{Nothing, SciMLBase.DEProblem}
    sys::PSY.System
    perturbations::Vector{<:Perturbation}
    x0_init::Vector{Float64}
    initialized::Bool
    tstops::Vector{Float64}
    callbacks::DiffEqBase.CallbackSet
    solution::Union{Nothing, SciMLBase.AbstractODESolution}
    simulation_folder::String
    simulation_inputs::SimulationInputs
    console_level::Base.CoreLogging.LogLevel
    file_level::Base.CoreLogging.LogLevel
    multimachine::Bool
end

get_system(sim::Simulation) = sim.sys
get_simulation_inputs(sim::Simulation) = sim.simulation_inputs

function Simulation(
    ::Type{T},
    sys::PSY.System;
    tspan,
    initial_conditions,
    perturbations = Vector{Perturbation}(),
    simulation_folder::String = "",
    console_level = Logging.Warn,
    file_level = Logging.Debug,
) where {T <: SimulationModel}
    PSY.set_units_base_system!(sys, "DEVICE_BASE")
    simulation_inputs = SimulationInputs(sys, tspan)

    if isempty(initial_conditions)
        @debug "Initial Condition set to flat start"
        bus_count = get_bus_count(simulation_inputs)
        var_count = get_variable_count(simulation_inputs)
        initial_conditions = zeros(var_count)
        initial_conditions[1:bus_count] .= 1.0
    end

    return Simulation{T}(
        BUILD_INCOMPLETE,
        nothing,
        sys,
        perturbations,
        initial_conditions,
        false,
        Vector{Float64}(),
        DiffEqBase.CallbackSet(),
        nothing,
        simulation_folder,
        simulation_inputs,
        console_level,
        file_level,
        false,
    )
end

function Simulation!(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
) where {T <: SimulationModel}
    return Simulation!(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
end

function Simulation(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbation::Perturbation;
    kwargs...,
) where {T <: SimulationModel}
    return Simulation(T, system, simulation_folder, tspan, [perturbation]; kwargs...)
end

"""
Builds the simulation object and conducts the indexing process. The initial conditions are stored in the system.

# Accepted Key Words
- `initialize_simulation::Bool : Runs the initialization routine. If false, simulation runs based on the operation point stored in System`
- `system_to_file::Bool`: Serializes the initialized system
- `console_level::Logging`: Sets the level of logging output to the console. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
- `file_level::Logging`: Sets the level of logging output to file. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
"""
function Simulation!(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
) where {T <: SimulationModel}
    check_folder(simulation_folder)
    # Instantiates the Simulation object
    sim = Simulation(
        T,
        system;
        tspan = tspan,
        initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Debug),
    )
    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "initialized_system.json"))
    end
    return sim
end

"""
Initializes the simulations and builds the indexing. The input system is not modified during the initialization

# Accepted Key Words
- `initialize_simulation::Bool : Runs the initialization routine. If false, simulation runs based on the operation point stored in System`
- `system_to_file::Bool`: Serializes the original input system
- `console_level::Logging`: Sets the level of logging output to the console. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debugg
- `file_level::Logging`: Sets the level of logging output to file. Can be set to Logging.Error, Logging.Warn, Logging.Info or Logging.Debug
"""
function Simulation(
    ::Type{T},
    system::PSY.System,
    simulation_folder::String,
    tspan::NTuple{2, Float64},
    perturbations::Vector{<:Perturbation} = Vector{Perturbation}();
    kwargs...,
) where {T <: SimulationModel}
    check_folder(simulation_folder)
    simulation_system = deepcopy(system)
    # Instantiates the Simulation object
    sim = Simulation(
        T,
        simulation_system;
        tspan = tspan,
        initial_conditions = get(kwargs, :initial_conditions, Vector{Float64}()),
        simulation_folder = simulation_folder,
        perturbations = perturbations,
        console_level = get(kwargs, :console_level, Logging.Warn),
        file_level = get(kwargs, :file_level, Logging.Debug),
    )
    build!(sim; kwargs...)
    if get(kwargs, :system_to_file, false)
        PSY.to_json(system, joinpath(simulation_folder, "input_system.json"))
    end
    return sim
end

function reset!(sim::Simulation{T}) where {T <: SimulationModel}
    @info "Rebuilding the simulation after reset"
    sim.simulation_inputs = SimulationInputs(get_system(sim), sim.simulation_inputs.tspan)
    build!(sim)
    @info "Simulation reset to status $(sim.status)"
    return
end

function configure_logging(sim::Simulation, file_mode)
    return IS.configure_logging(
        console = true,
        console_stream = stderr,
        console_level = sim.console_level,
        file = true,
        filename = joinpath(sim.simulation_folder, SIMULATION_LOG_FILENAME),
        file_level = sim.file_level,
        file_mode = file_mode,
        tracker = nothing,
        set_global = false,
    )
end

function _build_inputs!(sim::Simulation{T}) where {T <: SimulationModel}
    simulation_system = get_system(sim)
    sim.status = BUILD_INCOMPLETE
    build!(sim.simulation_inputs, T, simulation_system)
    @debug "Simulation Inputs Created"
    return
end

function _build_perturbations!(sim::Simulation)
    @info "Attaching Perturbations"
    if isempty(sim.perturbations)
        @debug "The simulation has no perturbations"
        return DiffEqBase.CallbackSet(), [0.0]
    end
    system = get_system(sim)
    perturbations = sim.perturbations
    perturbations_count = length(perturbations)
    callback_vector = Vector{DiffEqBase.DiscreteCallback}(undef, perturbations_count)
    tstops = Vector{Float64}(undef, perturbations_count)
    for (ix, pert) in enumerate(perturbations)
        @debug pert
        condition = (x, t, integrator) -> t in [pert.time]
        affect = get_affect(system, pert)
        callback_vector[ix] = DiffEqBase.DiscreteCallback(condition, affect)
        tstops[ix] = pert.time
    end
    sim.tstops = tstops
    sim.callbacks = DiffEqBase.CallbackSet((), tuple(callback_vector...))
    return
end

function _initialize_simulation!(sim::Simulation; kwargs...)
    if get(kwargs, :initialize_simulation, true)
        @info("Initializing Simulation States")
        sim.initialized = calculate_initial_conditions!(sim)
    else
        sim.initialized = true
        sim.status = SIMULATION_INITIALIZED
    end
    return
end

function _build!(sim::Simulation{ImplicitModel}; kwargs...)
    check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
    sim.status = BUILD_INCOMPLETE
    _build_inputs!(sim)
    _initialize_simulation!(sim; kwargs...)
    _build_perturbations!(sim)
    if sim.status != BUILD_FAILED
        try
            simulation_inputs = get_simulation_inputs(sim)
            simulation_pre_step!(simulation_inputs, Float64)
            var_count = get_variable_count(simulation_inputs)
            dx0 = zeros(Float64, var_count)
            sim.problem = SciMLBase.DAEProblem(
                system_implicit!,
                dx0,
                sim.x0_init,
                get_tspan(sim.simulation_inputs),
                simulation_inputs,
                differential_vars = get_DAE_vector(simulation_inputs);
                kwargs...,
            )
            sim.multimachine = (get_global_vars(simulation_inputs)[:ω_sys_index] != 0)
            sim.status = BUILT
            @info "Simulations status = $(sim.status)"
        catch e
            bt = catch_backtrace()
            @error "DiffEq DAEProblem failed to build" exception = e, bt
            sim.status = BUILD_FAILED
        end
    else
        @error "The simulation couldn't be initialized correctly. Simulations status = $(sim.status)"
    end
    return
end

function _build!(sim::Simulation{MassMatrixModel}; kwargs...)
    check_kwargs(kwargs, SIMULATION_ACCEPTED_KWARGS, "Simulation")
    sim.status = BUILD_INCOMPLETE
    _build_inputs!(sim)
    _initialize_simulation!(sim; kwargs...)
    _build_perturbations!(sim)

    if sim.status != BUILD_FAILED
        try
            simulation_inputs = get_simulation_inputs(sim)
            # TODO: Not use Real here. It has a bad performance hit.
            simulation_pre_step!(simulation_inputs, Real)
            sim.problem = SciMLBase.ODEProblem(
                SciMLBase.ODEFunction(
                    system_mass_matrix!,
                    mass_matrix = get_mass_matrix(simulation_inputs),
                ),
                sim.x0_init,
                get_tspan(sim.simulation_inputs),
                simulation_inputs;
                kwargs...,
            )
            sim.multimachine = (get_global_vars(simulation_inputs)[:ω_sys_index] != 0)
            sim.status = BUILT
            @info "Simulations status = $(sim.status)"
        catch e
            bt = catch_backtrace()
            @error "DiffEq ODEProblem failed to build" exception = e, bt
            sim.status = BUILD_FAILED
        end
    else
        @error "The simulation couldn't be initialized correctly. Simulations status = $(sim.status)"
    end
    return
end

function build!(sim; kwargs...)
    logger = configure_logging(sim, "w")
    try
        Logging.with_logger(logger) do
            _build!(sim; kwargs...)
        end
    finally
        close(logger)
    end
    return
end

function simulation_pre_step!(
    sim::Simulation,
    reset_simulation::Bool,
    ::Type{T},
) where {T <: Real}
    reset_simulation = sim.status == CONVERTED_FOR_SMALL_SIGNAL || reset_simulation
    if sim.status == BUILD_FAILED
        error(
            "The Simulation status is $(sim.status). Can not continue, correct your inputs and build the simulation again.",
        )
    elseif sim.status != BUILT && !reset_simulation
        error(
            "The Simulation status is $(sim.status). Use keyword argument reset_simulation = true",
        )
    end

    reset_simulation && reset!(sim)
    simulation_pre_step!(sim.simulation_inputs, T)
    return
end

function execute!(sim::Simulation{ImplicitModel}, solver; kwargs...)
    @debug "status before execute" sim.status
    simulation_pre_step!(sim, get(kwargs, :reset_simulation, false), Float64)
    sim.status = SIMULATION_STARTED
    sim.solution = SciMLBase.solve(
        sim.problem,
        solver;
        callback = sim.callbacks,
        tstops = sim.tstops,
        kwargs...,
    )
    if sim.solution.retcode == :Success
        sim.status = SIMULATION_FINALIZED
    else
        sim.status = SIMULATION_FAILED
    end
    return
end

function execute!(sim::Simulation{MassMatrixModel}, solver; kwargs...)
    @debug "status before execute" sim.status
    simulation_pre_step!(sim, get(kwargs, :reset_simulation, false), Real)
    sim.status = SIMULATION_STARTED
    sim.solution = SciMLBase.solve(
        sim.problem,
        solver;
        callback = sim.callbacks,
        tstops = sim.tstops,
        kwargs...,
    )
    if sim.solution.retcode == :Success
        sim.status = SIMULATION_FINALIZED
    else
        sim.status = SIMULATION_FAILED
    end
    return
end


text:
mutable struct gpuCTM <: TopicModel
	"GPU accelerated correlated topic model."

	K::Int
	M::Int
	V::Int
	N::Vector{Int}
	C::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	mu::Vector{Float32}
	sigma::Symmetric{Float32}
	invsigma::Symmetric{Float32}
	beta::Matrix{Float32}
	lambda::VectorList{Float32}
	lambda_dist::Vector{Float32}
	vsq::VectorList{Float32}
	logzeta::Vector{Float32}
	phi::MatrixList{Float32}
	elbo::Float32
	device::cl.Device
	context::cl.Context
	queue::cl.CmdQueue
	mu_kernel::cl.Kernel
	sigma_kernel::cl.Kernel
	beta_kernel::cl.Kernel
	beta_norm_kernel::cl.Kernel
	lambda_kernel::cl.Kernel
	vsq_kernel::cl.Kernel
	logzeta_kernel::cl.Kernel
	phi_kernel::cl.Kernel
	phi_norm_kernel::cl.Kernel
	C_buffer::cl.Buffer{Int}
	N_cumsum_buffer::cl.Buffer{Int}
	J_cumsum_buffer::cl.Buffer{Int}
	terms_buffer::cl.Buffer{Int}
	terms_sortperm_buffer::cl.Buffer{Int}
	counts_buffer::cl.Buffer{Int}
	mu_buffer::cl.Buffer{Float32}
	sigma_buffer::cl.Buffer{Float32}
	invsigma_buffer::cl.Buffer{Float32}
	beta_buffer::cl.Buffer{Float32}
	lambda_buffer::cl.Buffer{Float32}
	lambda_dist_buffer::cl.Buffer{Float32}
	vsq_buffer::cl.Buffer{Float32}
	logzeta_buffer::cl.Buffer{Float32}
	phi_buffer::cl.Buffer{Float32}
	phi_count_buffer::cl.Buffer{Float32}

	function gpuCTM(corp::Corpus, K::Integer)
		check_corp(corp)
		K > 0 || throw(ArgumentError("Number of topics must be a positive integer."))

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]

		mu = zeros(Float32, K)
		sigma = Symmetric(Matrix{Float32}(I, K, K))
		invsigma = copy(sigma)
		beta = rand(Dirichlet(V, 1.0f0), K)'
		lambda = [zeros(Float32, K) for _ in 1:M]
		lambda_dist = zeros(Float32, M)
		vsq = [ones(Float32, K) for _ in 1:M]
		logzeta = fill(0.5f0, M)
		phi = [fill(Float32(1/K), K, N[d]) for d in 1:M]
		elbo = 0f0

		device, context, queue = cl.create_compute_context()

		mu_program = cl.Program(context, source=CTM_MU_c) |> cl.build!
		sigma_program = cl.Program(context, source=CTM_SIGMA_c) |> cl.build!
		beta_program = cl.Program(context, source=CTM_BETA_c) |> cl.build!
		beta_norm_program = cl.Program(context, source=CTM_BETA_NORM_c) |> cl.build!
		lambda_program = cl.Program(context, source=CTM_LAMBDA_c) |> cl.build!
		vsq_program = cl.Program(context, source=CTM_VSQ_c) |> cl.build!
		logzeta_program = cl.Program(context, source=CTM_LOGZETA_c) |> cl.build!
		phi_program = cl.Program(context, source=CTM_PHI_c) |> cl.build!
		phi_norm_program = cl.Program(context, source=CTM_PHI_NORM_c) |> cl.build!

		mu_kernel = cl.Kernel(mu_program, "update_mu")
		sigma_kernel = cl.Kernel(sigma_program, "update_sigma")
		beta_kernel = cl.Kernel(beta_program, "update_beta")
		beta_norm_kernel = cl.Kernel(beta_norm_program, "normalize_beta")
		lambda_kernel = cl.Kernel(lambda_program, "update_lambda")
		vsq_kernel = cl.Kernel(vsq_program, "update_vsq")
		logzeta_kernel = cl.Kernel(logzeta_program, "update_logzeta")
		phi_kernel = cl.Kernel(phi_program, "update_phi")
		phi_norm_kernel = cl.Kernel(phi_norm_program, "normalize_phi")

		model = new(K, M, V, N, C, copy(corp), topics, mu, sigma, invsigma, beta, lambda, lambda_dist, vsq, logzeta, phi, elbo, device, context, queue, mu_kernel, sigma_kernel, beta_kernel, beta_norm_kernel, lambda_kernel, vsq_kernel, logzeta_kernel, phi_kernel, phi_norm_kernel)
		return model
	end
end

function Elogpeta(model::gpuCTM, d::Int)
	"Compute E_q[log(P(eta))]."

	x = 0.5 * (logdet(model.invsigma) - model.K * log(2pi) - dot(diag(model.invsigma), model.vsq[d]) - dot(model.lambda[d] - model.mu, model.invsigma * (model.lambda[d] - model.mu)))
	return x
end

function Elogpz(model::gpuCTM, d::Int)
	"Compute E_q[log(P(z))]."

	counts = model.corp[d].counts
	x = dot(model.phi[d]' * model.lambda[d], counts) - model.C[d] * (sum(exp.(model.lambda[d] + 0.5 * model.vsq[d] .- model.logzeta[d])) + model.logzeta[d] - 1)
	return x
end

function Elogpw(model::gpuCTM, d::Int)
	"Compute E_q[log(P(w))]."

	terms, counts = model.corp[d].terms, model.corp[d].counts
	x = sum(model.phi[d] .* log.(@boink model.beta[:,terms]) * counts)
	return x
end

function Elogqeta(model::gpuCTM, d::Int)
	"Compute E_q[log(q(eta))]."

	x = -entropy(MvNormal(model.lambda[d], diagm(model.vsq[d])))
	return x
end

function Elogqz(model::gpuCTM, d::Int)
	"Compute E_q[log(q(z))]."

	counts = model.corp[d].counts
	x = -sum([c * entropy(Categorical(model.phi[d][:,n])) for (n, c) in enumerate(counts)])
	return x
end

function update_elbo!(model::gpuCTM)
	"Update the evidence lower bound."

	model.elbo = 0
	for d in 1:model.M
		model.elbo += Elogpeta(model, d) + Elogpz(model, d) + Elogpw(model, d) - Elogqeta(model, d) - Elogqz(model, d)
	end

	return model.elbo
end

const CTM_MU_c =
"""
kernel void
update_mu(	long K,
			long M,
			const global float *lambda,
			global float *mu)

			{
			long i = get_global_id(0);

			float acc = 0.0f;

			for (long d=0; d<M; d++)
				acc += lambda[K * d + i];

			mu[i] = acc / M;
			}
			"""

function update_mu!(model::gpuCTM)
	"Update mu."
	"Analytic."

	model.queue(model.mu_kernel, model.K, nothing, model.K, model.M, model.lambda_buffer, model.mu_buffer)
end

const CTM_SIGMA_c =
"""
kernel void
update_sigma( 	long K,
		        long M,
		        const global float *mu,
		        const global float *lambda,
		        const global float *vsq,
		        global float *sigma)

		        {
		        long i = get_global_id(0);
		        long j = get_global_id(1);

		        float acc = 0.0f;

		        if (i == j)
		            for (long d=0; d<M; d++)
		              acc += vsq[K * d + i] + (lambda[K * d + i] - mu[i]) * (lambda[K * d + j] - mu[j]);

		        if (i != j)
		            for (long d=0; d<M; d++)
		              acc += (lambda[K * d + i] - mu[i]) * (lambda[K * d + j] - mu[j]);

		        sigma[K * i + j] = acc / M;
		        }
		        """

function update_sigma!(model::gpuCTM)
	"Update sigma."
	"Analytic"

	model.queue(model.sigma_kernel, (model.K, model.K), nothing, model.K, model.M, model.mu_buffer, model.lambda_buffer, model.vsq_buffer, model.sigma_buffer)

	@host model.sigma_buffer
	model.invsigma = inv(model.sigma)
	@buffer model.invsigma
end

const CTM_BETA_c =
"""
kernel void
update_beta(long K,
			const global long *J_cumsum,
			const global long *terms_sortperm,
			const global long *counts,
			const global float *phi,
			global float *beta)

			{
			long i = get_global_id(0);
			long j = get_global_id(1);

			float acc = 0.0f;

			for (long w=J_cumsum[j]; w<J_cumsum[j+1]; w++)
				acc += counts[terms_sortperm[w]] * phi[K * terms_sortperm[w] + i];

			beta[K * j + i] = acc;
			}
			"""

const CTM_BETA_NORM_c =
"""
kernel void
normalize_beta(	long K,
				long V,
				global float *beta)

				{
				long i = get_global_id(0);

				float normalizer = 0.0f;

				for (long j=0; j<V; j++)
					normalizer += beta[K * j + i];

				for (long j=0; j<V; j++)
					beta[K * j + i] /= normalizer;
				}
				"""

function update_beta!(model::gpuCTM)
	"Update beta."
	"Analytic."

	model.queue(model.beta_kernel, (model.K, model.V), nothing, model.K, model.J_cumsum_buffer, model.terms_sortperm_buffer, model.counts_buffer, model.phi_buffer, model.beta_buffer)
	model.queue(model.beta_norm_kernel, model.K, nothing, model.K, model.V, model.beta_buffer)
end

const CTM_LAMBDA_c =
"""
$(LINSOLVE_c)

kernel void
update_lambda(	long niter,
				float ntol,
				long K,
				global float *phi_count,
				const global long *C,
				const global long *N_cumsum,
				const global long *counts,
				const global float *mu,
				const global float *invsigma,
				const global float *vsq,
				const global float *logzeta,
				const global float *phi,
				global float *lambda,
				global float *lambda_dist,
				local float *lambda_old,
				local float *lambda_grad,
				local float *lambda_hess)

				{
				long d = get_group_id(0);
				long z = get_local_id(0);

				float acc;

				lambda_old[z] = lambda[K * d + z];

				phi_count[K * d + z] = 0.0f;
				for (long n=N_cumsum[d]; n<N_cumsum[d+1]; n++)
					phi_count[K * d + z] += phi[K * n + z] * counts[n];

				barrier(CLK_LOCAL_MEM_FENCE);

				for (long _=0; _<niter; _++)
				{
					acc = 0.0f;
					for (long j=0; j<K; j++)
						acc += invsigma[K * j + z] * (mu[j] - lambda[K * d + j]);

					lambda_grad[z] = acc + phi_count[K * d + z] - C[d] * exp(lambda[K * d + z] + 0.5f * vsq[K * d + z] - logzeta[d]);

					for (long j=0; j<K; j++)
						lambda_hess[K * j + z] = -invsigma[K * j + z];

					lambda_hess[K * z + z] -= C[d] * exp(lambda[K * d + z] + 0.5f * vsq[K * d + z] - logzeta[d]);

					barrier(CLK_LOCAL_MEM_FENCE);

					acc = 0.0f;
					for (long i=0; i<K; i++)
						acc += lambda_grad[i] * lambda_grad[i];

					barrier(CLK_LOCAL_MEM_FENCE);

					linsolve(K, z, lambda_hess, lambda_grad);

					barrier(CLK_LOCAL_MEM_FENCE);

					lambda[K * d + z] -= lambda_grad[z];

					barrier(CLK_LOCAL_MEM_FENCE);

					if (sqrt(acc) < ntol)
						break;
				}

				if (z == 0)
				{
					acc = 0.0f;
					for (long i=0; i<K; i++)
						acc += pow(lambda[K * d + i] - lambda_old[i], 2);

					lambda_dist[d] = sqrt(acc);
				}
				}
				"""

function update_lambda!(model::gpuCTM, niter::Int, ntol::Float32)
	"Update lambda."
	"Newton's method."

	model.queue(model.lambda_kernel, model.M * model.K, model.K, niter, ntol, model.K, model.phi_count_buffer, model.C_buffer, model.N_cumsum_buffer, model.counts_buffer, model.mu_buffer, model.invsigma_buffer, model.vsq_buffer, model.logzeta_buffer, model.phi_buffer, model.lambda_buffer, model.lambda_dist_buffer)
	@host model.lambda_dist_buffer
end

const CTM_VSQ_c =
"""
kernel void
update_vsq(	long niter,
			float ntol,
			long K,
			const global long *C,
			const global float *invsigma,
			const global float *lambda,
			const global float *logzeta,
			global float *vsq)

			{
			long d = get_global_id(0);
			long i = get_global_id(1);

			float vsq_grad;
			float vsq_invhess;
			float p;

			for (long _=0; _<niter; _++)
			{
				float rho = 1.0f;

				vsq_grad = -0.5f * (invsigma[K * i + i] + C[d] * exp(lambda[K * d + i] + 0.5f * vsq[K * d + i] - logzeta[d]) - 1 / vsq[K * d + i]);
				vsq_invhess = -1 / (0.25f * C[d] * exp(lambda[K * d + i] + 0.5f * vsq[K * d + i] - logzeta[d]) + 0.5f / (vsq[K * d + i] * vsq[K * d + i]));
				p = vsq_grad * vsq_invhess;

				while (vsq[K * d + i] - rho * p <= 0)
					rho *= 0.5f;

				vsq[K * d + i] -= rho * p;

				if (rho * fabs(vsq_grad) < ntol)
					break;
			}

			vsq[K * d + i] += $(EPSILON32);
			}
			"""

function update_vsq!(model::gpuCTM, niter::Int, ntol::Float32)
	"Update vsq."
	"Newton's method with back-tracking line search."

	model.queue(model.vsq_kernel, (model.M, model.K), nothing, niter, ntol, model.K, model.C_buffer, model.invsigma_buffer, model.lambda_buffer, model.logzeta_buffer, model.vsq_buffer)
end

const CTM_LOGZETA_c =
"""
kernel void
update_logzeta(	long K,
				const global float *lambda,
				const global float *vsq,
				global float *logzeta)

				{
				long d = get_global_id(0);

				float maxval = 0.0f;

				for (long i=0; i<K; i++)
				{
					float x = lambda[K * d + i] + 0.5f * vsq[K * d + i];
					if (x > maxval)
						maxval = x;
				}

				float acc = 0.0f;

				for (long i=0; i<K; i++)
					acc += exp(lambda[K * d + i] + 0.5f * vsq[K * d + i] - maxval);

				logzeta[d] = maxval + log(acc);
				}
				"""

function update_logzeta!(model::gpuCTM)
	"Update logzeta."
	"Analytic."

	model.queue(model.logzeta_kernel, model.M, nothing, model.K, model.lambda_buffer, model.vsq_buffer, model.logzeta_buffer)
end

const CTM_PHI_c =
"""
kernel void
update_phi(	long K,
			const global long *N_cumsum,
			const global long *terms,
			const global float *beta,
			const global float *lambda,
			global float *phi)

			{
			long i = get_global_id(0);
			long d = get_global_id(1);

			for (long n=N_cumsum[d]; n<N_cumsum[d+1]; n++)
				phi[K * n + i] = log(beta[K * terms[n] + i]) + lambda[K * d + i] + $(EPSILON32);
			}
			"""

const CTM_PHI_NORM_c =
"""
kernel void
normalize_phi(	long K,
				global float *phi)

				{
				long dn = get_global_id(0);

				float maxval = 0.0f;

				for (long i=0; i<K; i++)
					if (phi[K * dn + i] > maxval)
						maxval = phi[K * dn + i];

				float normalizer = 0.0f;

				for (long i=0; i<K; i++)
				{
					phi[K * dn + i] -= maxval;
					normalizer += exp(phi[K * dn + i]);
				}

				for (long i=0; i<K; i++)
					phi[K * dn + i] = exp(phi[K * dn + i]) / normalizer;
				}
				"""

function update_phi!(model::gpuCTM)
	"Update phi."
	"Analytic."

	model.queue(model.phi_kernel, (model.K, model.M), nothing, model.K, model.N_cumsum_buffer, model.terms_buffer, model.beta_buffer, model.lambda_buffer, model.phi_buffer)
	model.queue(model.phi_norm_kernel, sum(model.N), nothing, model.K, model.phi_buffer)
end

function train!(model::gpuCTM; iter::Integer=150, tol::Real=1.0, niter::Integer=1000, ntol::Real=1/model.K^2, viter::Integer=10, vtol::Real=1/model.K^2, checkelbo::Real=1, printelbo::Bool=true)
	"Coordinate ascent optimization procedure for GPU accelerated correlated topic model variational Bayes algorithm."

	ntol = Float32(ntol)
	check_model(model)
	all([tol, ntol, vtol] .>= 0)										|| throw(ArgumentError("Tolerance parameters must be nonnegative."))
	all([iter, niter, viter] .>= 0)										|| throw(ArgumentError("Iteration parameters must be nonnegative."))
	(isa(checkelbo, Integer) & (checkelbo > 0)) | (checkelbo == Inf) 	|| throw(ArgumentError("checkelbo parameter must be a positive integer or Inf."))
	all([isempty(doc) for doc in model.corp]) ? (iter = 0) : update_buffer!(model)
	(checkelbo <= iter) && update_elbo!(model)

	for k in 1:iter
		for v in 1:viter
			update_phi!(model)
			update_logzeta!(model)
			update_vsq!(model, niter, ntol)
			update_lambda!(model, niter, ntol)

			if median(model.lambda_dist) < vtol
				break
			end
		end
		update_beta!(model)
		update_sigma!(model)
		update_mu!(model)

		if check_elbo!(model, checkelbo, printelbo, k, tol)
			break
		end
	end

	(iter > 0) && update_host!(model)
	model.topics = [reverse(sortperm(vec(model.beta[i,:]))) for i in 1:model.K]
	nothing
end


text:
using Pkg
using Test

@testset "Dubins_package" begin
    @test run(`julia src/Dubins/test.jl`).exitcode == 0
    snoop_file = "src/Dubins/test-snoop.txt"
    @test isfile(snoop_file)
    lines = readlines(open(snoop_file))
    count1 = count(x -> occursin("Dubins", x), lines)
    @test count1 > 1
    @test run(`julia -J src/Dubins/sys.so src/Dubins/check.jl`).exitcode == 0
    snoop_file = "src/Dubins/check-snoop.txt"
    lines = readlines(open(snoop_file))
    count2 = count(x -> occursin("Dubins", x), lines)
    @test count2 < count1
    rm("src/Dubins/sys.so")
    rm("src/Dubins/sys.a")
    rm("src/Dubins/test-snoop.txt")
    rm("src/Dubins/check-snoop.txt")
end


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

"""
The `AbstractString` type is the supertype of all string implementations in
Julia. Strings are encodings of sequences of [Unicode](https://unicode.org/)
code points as represented by the `AbstractChar` type. Julia makes a few assumptions
about strings:

* Strings are encoded in terms of fixed-size "code units"
  * Code units can be extracted with `codeunit(s, i)`
  * The first code unit has index `1`
  * The last code unit has index `ncodeunits(s)`
  * Any index `i` such that `1 ≤ i ≤ ncodeunits(s)` is in bounds
* String indexing is done in terms of these code units:
  * Characters are extracted by `s[i]` with a valid string index `i`
  * Each `AbstractChar` in a string is encoded by one or more code units
  * Only the index of the first code unit of an `AbstractChar` is a valid index
  * The encoding of an `AbstractChar` is independent of what precedes or follows it
  * String encodings are [self-synchronizing] – i.e. `isvalid(s, i)` is O(1)

[self-synchronizing]: https://en.wikipedia.org/wiki/Self-synchronizing_code

Some string functions that extract code units, characters or substrings from
strings error if you pass them out-of-bounds or invalid string indices. This
includes `codeunit(s, i)` and `s[i]`. Functions that do string
index arithmetic take a more relaxed approach to indexing and give you the
closest valid string index when in-bounds, or when out-of-bounds, behave as if
there were an infinite number of characters padding each side of the string.
Usually these imaginary padding characters have code unit length `1` but string
types may choose different "imaginary" character sizes as makes sense for their
implementations (e.g. substrings may pass index arithmetic through to the
underlying string they provide a view into). Relaxed indexing functions include
those intended for index arithmetic: `thisind`, `nextind` and `prevind`. This
model allows index arithmetic to work with out-of- bounds indices as
intermediate values so long as one never uses them to retrieve a character,
which often helps avoid needing to code around edge cases.

See also: [`codeunit`](@ref), [`ncodeunits`](@ref), [`thisind`](@ref),
[`nextind`](@ref), [`prevind`](@ref)
"""
AbstractString

## required string functions ##

"""
    ncodeunits(s::AbstractString) -> Int

Return the number of code units in a string. Indices that are in bounds to
access this string must satisfy `1 ≤ i ≤ ncodeunits(s)`. Not all such indices
are valid – they may not be the start of a character, but they will return a
code unit value when calling `codeunit(s,i)`.

See also: [`codeunit`](@ref), [`checkbounds`](@ref), [`sizeof`](@ref),
[`length`](@ref), [`lastindex`](@ref)
"""
ncodeunits(s::AbstractString)

"""
    codeunit(s::AbstractString) -> Type{<:Union{UInt8, UInt16, UInt32}}

Return the code unit type of the given string object. For ASCII, Latin-1, or
UTF-8 encoded strings, this would be `UInt8`; for UCS-2 and UTF-16 it would be
`UInt16`; for UTF-32 it would be `UInt32`. The unit code type need not be
limited to these three types, but it's hard to think of widely used string
encodings that don't use one of these units. `codeunit(s)` is the same as
`typeof(codeunit(s,1))` when `s` is a non-empty string.

See also: [`ncodeunits`](@ref)
"""
codeunit(s::AbstractString)

"""
    codeunit(s::AbstractString, i::Integer) -> Union{UInt8, UInt16, UInt32}

Return the code unit value in the string `s` at index `i`. Note that

    codeunit(s, i) :: codeunit(s)

I.e. the value returned by `codeunit(s, i)` is of the type returned by
`codeunit(s)`.

See also: [`ncodeunits`](@ref), [`checkbounds`](@ref)
"""
@propagate_inbounds codeunit(s::AbstractString, i::Integer) = typeof(i) === Int ?
    throw(MethodError(codeunit, (s, i))) : codeunit(s, Int(i))

"""
    isvalid(s::AbstractString, i::Integer) -> Bool

Predicate indicating whether the given index is the start of the encoding of a
character in `s` or not. If `isvalid(s, i)` is true then `s[i]` will return the
character whose encoding starts at that index, if it's false, then `s[i]` will
raise an invalid index error or a bounds error depending on if `i` is in bounds.
In order for `isvalid(s, i)` to be an O(1) function, the encoding of `s` must be
[self-synchronizing](https://en.wikipedia.org/wiki/Self-synchronizing_code) this
is a basic assumption of Julia's generic string support.

See also: [`getindex`](@ref), [`iterate`](@ref), [`thisind`](@ref),
[`nextind`](@ref), [`prevind`](@ref), [`length`](@ref)

# Examples

```jldoctest
julia> str = "αβγdef";

julia> isvalid(str, 1)
true

julia> str[1]
'α': Unicode U+03B1 (category Ll: Letter, lowercase)

julia> isvalid(str, 2)
false

julia> str[2]
ERROR: StringIndexError("αβγdef", 2)
Stacktrace:
[...]
```
"""
@propagate_inbounds isvalid(s::AbstractString, i::Integer) = typeof(i) === Int ?
    throw(MethodError(isvalid, (s, i))) : isvalid(s, Int(i))

"""
    iterate(s::AbstractString, i::Integer) -> Union{Tuple{<:AbstractChar, Int}, Nothing}

Return a tuple of the character in `s` at index `i` with the index of the start
of the following character in `s`. This is the key method that allows strings to
be iterated, yielding a sequences of characters. If `i` is out of bounds in `s`
then a bounds error is raised. The `iterate` function, as part of the iteration
protocol may assume that `i` is the start of a character in `s`.

See also: [`getindex`](@ref), [`checkbounds`](@ref)
"""
@propagate_inbounds iterate(s::AbstractString, i::Integer) = typeof(i) === Int ?
    throw(MethodError(iterate, (s, i))) : iterate(s, Int(i))

## basic generic definitions ##

eltype(::Type{<:AbstractString}) = Char # some string types may use another AbstractChar

"""
    sizeof(str::AbstractString)

Size, in bytes, of the string `str`. Equal to the number of code units in `str` multiplied by
the size, in bytes, of one code unit in `str`.

# Examples
```jldoctest
julia> sizeof("")
0

julia> sizeof("∀")
3
```
"""
sizeof(s::AbstractString) = ncodeunits(s) * sizeof(codeunit(s))
firstindex(s::AbstractString) = 1
lastindex(s::AbstractString) = thisind(s, ncodeunits(s))

function getindex(s::AbstractString, i::Integer)
    @boundscheck checkbounds(s, i)
    @inbounds return isvalid(s, i) ? iterate(s, i)[1] : string_index_err(s, i)
end

getindex(s::AbstractString, i::Colon) = s
# TODO: handle other ranges with stride ±1 specially?
# TODO: add more @propagate_inbounds annotations?
getindex(s::AbstractString, v::AbstractVector{<:Integer}) =
    sprint(io->(for i in v; write(io, s[i]) end), sizehint=length(v))
getindex(s::AbstractString, v::AbstractVector{Bool}) =
    throw(ArgumentError("logical indexing not supported for strings"))

function get(s::AbstractString, i::Integer, default)
# TODO: use ternary once @inbounds is expression-like
    if checkbounds(Bool, s, i)
        @inbounds return s[i]
    else
        return default
    end
end

## bounds checking ##

checkbounds(::Type{Bool}, s::AbstractString, i::Integer) =
    1 ≤ i ≤ ncodeunits(s)
checkbounds(::Type{Bool}, s::AbstractString, r::AbstractRange{<:Integer}) =
    isempty(r) || (1 ≤ minimum(r) && maximum(r) ≤ ncodeunits(s))
checkbounds(::Type{Bool}, s::AbstractString, I::AbstractArray{<:Real}) =
    all(i -> checkbounds(Bool, s, i), I)
checkbounds(::Type{Bool}, s::AbstractString, I::AbstractArray{<:Integer}) =
    all(i -> checkbounds(Bool, s, i), I)
checkbounds(s::AbstractString, I::Union{Integer,AbstractArray}) =
    checkbounds(Bool, s, I) ? nothing : throw(BoundsError(s, I))

## construction, conversion, promotion ##

string() = ""
string(s::AbstractString) = s

(::Type{Vector{UInt8}})(s::AbstractString) = unsafe_wrap(Vector{UInt8}, String(s))
(::Type{Array{UInt8}})(s::AbstractString) = unsafe_wrap(Vector{UInt8}, String(s))
(::Type{Vector{T}})(s::AbstractString) where {T<:AbstractChar} = collect(T, s)

Symbol(s::AbstractString) = Symbol(String(s))
Symbol(x...) = Symbol(string(x...))

convert(::Type{T}, s::T) where {T<:AbstractString} = s
convert(::Type{T}, s::AbstractString) where {T<:AbstractString} = T(s)

## string & character concatenation ##

"""
    *(s::Union{AbstractString, AbstractChar}, t::Union{AbstractString, AbstractChar}...) -> AbstractString

Concatenate strings and/or characters, producing a [`String`](@ref). This is equivalent
to calling the [`string`](@ref) function on the arguments. Concatenation of built-in
string types always produces a value of type `String` but other string types may choose
to return a string of a different type as appropriate.

# Examples
```jldoctest
julia> "Hello " * "world"
"Hello world"

julia> 'j' * "ulia"
"julia"
```
"""
(*)(s1::Union{AbstractChar, AbstractString}, ss::Union{AbstractChar, AbstractString}...) = string(s1, ss...)

one(::Union{T,Type{T}}) where {T<:AbstractString} = convert(T, "")

## generic string comparison ##

"""
    cmp(a::AbstractString, b::AbstractString) -> Int

Compare two strings. Return `0` if both strings have the same length and the character
at each index is the same in both strings. Return `-1` if `a` is a prefix of `b`, or if
`a` comes before `b` in alphabetical order. Return `1` if `b` is a prefix of `a`, or if
`b` comes before `a` in alphabetical order (technically, lexicographical order by Unicode
code points).

# Examples
```jldoctest
julia> cmp("abc", "abc")
0

julia> cmp("ab", "abc")
-1

julia> cmp("abc", "ab")
1

julia> cmp("ab", "ac")
-1

julia> cmp("ac", "ab")
1

julia> cmp("α", "a")
1

julia> cmp("b", "β")
-1
```
"""
function cmp(a::AbstractString, b::AbstractString)
    a === b && return 0
    a, b = Iterators.Stateful(a), Iterators.Stateful(b)
    for (c, d) in zip(a, b)
        c ≠ d && return ifelse(c < d, -1, 1)
    end
    isempty(a) && return ifelse(isempty(b), 0, -1)
    return 1
end

"""
    ==(a::AbstractString, b::AbstractString) -> Bool

Test whether two strings are equal character by character (technically, Unicode
code point by code point).

# Examples
```jldoctest
julia> "abc" == "abc"
true

julia> "abc" == "αβγ"
false
```
"""
==(a::AbstractString, b::AbstractString) = cmp(a, b) == 0

"""
    isless(a::AbstractString, b::AbstractString) -> Bool

Test whether string `a` comes before string `b` in alphabetical order
(technically, in lexicographical order by Unicode code points).

# Examples
```jldoctest
julia> isless("a", "b")
true

julia> isless("β", "α")
false

julia> isless("a", "a")
false
```
"""
isless(a::AbstractString, b::AbstractString) = cmp(a, b) < 0

# faster comparisons for symbols

cmp(a::Symbol, b::Symbol) = Int(sign(ccall(:strcmp, Int32, (Cstring, Cstring), a, b)))

isless(a::Symbol, b::Symbol) = cmp(a, b) < 0

## character index arithmetic ##

"""
    length(s::AbstractString) -> Int
    length(s::AbstractString, i::Integer, j::Integer) -> Int

The number of characters in string `s` from indices `i` through `j`. This is
computed as the number of code unit indices from `i` to `j` which are valid
character indices. With only a single string argument, this computes the
number of characters in the entire string. With `i` and `j` arguments it
computes the number of indices between `i` and `j` inclusive that are valid
indices in the string `s`. In addition to in-bounds values, `i` may take the
out-of-bounds value `ncodeunits(s) + 1` and `j` may take the out-of-bounds
value `0`.

See also: [`isvalid`](@ref), [`ncodeunits`](@ref), [`lastindex`](@ref),
[`thisind`](@ref), [`nextind`](@ref), [`prevind`](@ref)

# Examples
```jldoctest
julia> length("jμΛIα")
5
```
"""
length(s::AbstractString) = @inbounds return length(s, 1, ncodeunits(s))

function length(s::AbstractString, i::Int, j::Int)
    @boundscheck begin
        0 < i ≤ ncodeunits(s)+1 || throw(BoundsError(s, i))
        0 ≤ j < ncodeunits(s)+1 || throw(BoundsError(s, j))
    end
    n = 0
    for k = i:j
        @inbounds n += isvalid(s, k)
    end
    return n
end

@propagate_inbounds length(s::AbstractString, i::Integer, j::Integer) =
    length(s, Int(i), Int(j))

"""
    thisind(s::AbstractString, i::Integer) -> Int

If `i` is in bounds in `s` return the index of the start of the character whose
encoding code unit `i` is part of. In other words, if `i` is the start of a
character, return `i`; if `i` is not the start of a character, rewind until the
start of a character and return that index. If `i` is equal to 0 or `ncodeunits(s)+1`
return `i`. In all other cases throw `BoundsError`.

# Examples
```jldoctest
julia> thisind("α", 0)
0

julia> thisind("α", 1)
1

julia> thisind("α", 2)
1

julia> thisind("α", 3)
3

julia> thisind("α", 4)
ERROR: BoundsError: attempt to access String
  at index [4]
[...]

julia> thisind("α", -1)
ERROR: BoundsError: attempt to access String
  at index [-1]
[...]
```
"""
thisind(s::AbstractString, i::Integer) = thisind(s, Int(i))

function thisind(s::AbstractString, i::Int)
    z = ncodeunits(s) + 1
    i == z && return i
    @boundscheck 0 ≤ i ≤ z || throw(BoundsError(s, i))
    @inbounds while 1 < i && !isvalid(s, i)
        i -= 1
    end
    return i
end

"""
    prevind(str::AbstractString, i::Integer, n::Integer=1) -> Int

* Case `n == 1`

  If `i` is in bounds in `s` return the index of the start of the character whose
  encoding starts before index `i`. In other words, if `i` is the start of a
  character, return the start of the previous character; if `i` is not the start
  of a character, rewind until the start of a character and return that index.
  If `i` is equal to `1` return `0`.
  If `i` is equal to `ncodeunits(str)+1` return `lastindex(str)`.
  Otherwise throw `BoundsError`.

* Case `n > 1`

  Behaves like applying `n` times `prevind` for `n==1`. The only difference
  is that if `n` is so large that applying `prevind` would reach `0` then each remaining
  iteration decreases the returned value by `1`.
  This means that in this case `prevind` can return a negative value.

* Case `n == 0`

  Return `i` only if `i` is a valid index in `str` or is equal to `ncodeunits(str)+1`.
  Otherwise `StringIndexError` or `BoundsError` is thrown.

# Examples
```jldoctest
julia> prevind("α", 3)
1

julia> prevind("α", 1)
0

julia> prevind("α", 0)
ERROR: BoundsError: attempt to access String
  at index [0]
[...]

julia> prevind("α", 2, 2)
0

julia> prevind("α", 2, 3)
-1
```
"""
prevind(s::AbstractString, i::Integer, n::Integer) = prevind(s, Int(i), Int(n))
prevind(s::AbstractString, i::Integer)             = prevind(s, Int(i))
prevind(s::AbstractString, i::Int)                 = prevind(s, i, 1)

function prevind(s::AbstractString, i::Int, n::Int)
    n < 0 && throw(ArgumentError("n cannot be negative: $n"))
    z = ncodeunits(s) + 1
    @boundscheck 0 < i ≤ z || throw(BoundsError(s, i))
    n == 0 && return thisind(s, i) == i ? i : string_index_err(s, i)
    while n > 0 && 1 < i
        @inbounds n -= isvalid(s, i -= 1)
    end
    return i - n
end

"""
    nextind(str::AbstractString, i::Integer, n::Integer=1) -> Int

* Case `n == 1`

  If `i` is in bounds in `s` return the index of the start of the character whose
  encoding starts after index `i`. In other words, if `i` is the start of a
  character, return the start of the next character; if `i` is not the start
  of a character, move forward until the start of a character and return that index.
  If `i` is equal to `0` return `1`.
  If `i` is in bounds but greater or equal to `lastindex(str)` return `ncodeunits(str)+1`.
  Otherwise throw `BoundsError`.

* Case `n > 1`

  Behaves like applying `n` times `nextind` for `n==1`. The only difference
  is that if `n` is so large that applying `nextind` would reach `ncodeunits(str)+1` then
  each remaining iteration increases the returned value by `1`. This means that in this
  case `nextind` can return a value greater than `ncodeunits(str)+1`.

* Case `n == 0`

  Return `i` only if `i` is a valid index in `s` or is equal to `0`.
  Otherwise `StringIndexError` or `BoundsError` is thrown.

# Examples
```jldoctest
julia> nextind("α", 0)
1

julia> nextind("α", 1)
3

julia> nextind("α", 3)
ERROR: BoundsError: attempt to access String
  at index [3]
[...]

julia> nextind("α", 0, 2)
3

julia> nextind("α", 1, 2)
4
```
"""
nextind(s::AbstractString, i::Integer, n::Integer) = nextind(s, Int(i), Int(n))
nextind(s::AbstractString, i::Integer)             = nextind(s, Int(i))
nextind(s::AbstractString, i::Int)                 = nextind(s, i, 1)

function nextind(s::AbstractString, i::Int, n::Int)
    n < 0 && throw(ArgumentError("n cannot be negative: $n"))
    z = ncodeunits(s)
    @boundscheck 0 ≤ i ≤ z || throw(BoundsError(s, i))
    n == 0 && return thisind(s, i) == i ? i : string_index_err(s, i)
    while n > 0 && i < z
        @inbounds n -= isvalid(s, i += 1)
    end
    return i + n
end

## string index iteration type ##

struct EachStringIndex{T<:AbstractString}
    s::T
end
keys(s::AbstractString) = EachStringIndex(s)

length(e::EachStringIndex) = length(e.s)
first(::EachStringIndex) = 1
last(e::EachStringIndex) = lastindex(e.s)
iterate(e::EachStringIndex, state=firstindex(e.s)) = state > ncodeunits(e.s) ? nothing : (state, nextind(e.s, state))
eltype(::Type{<:EachStringIndex}) = Int

"""
    isascii(c::Union{AbstractChar,AbstractString}) -> Bool

Test whether a character belongs to the ASCII character set, or whether this is true for
all elements of a string.

# Examples
```jldoctest
julia> isascii('a')
true

julia> isascii('α')
false

julia> isascii("abc")
true

julia> isascii("αβγ")
false
```
"""
isascii(c::Char) = bswap(reinterpret(UInt32, c)) < 0x80
isascii(s::AbstractString) = all(isascii, s)
isascii(c::AbstractChar) = UInt32(c) < 0x80

## string map, filter, has ##

function map(f, s::AbstractString)
    out = IOBuffer(sizehint=sizeof(s))
    for c in s
        c′ = f(c)
        isa(c′, AbstractChar) || throw(ArgumentError(
            "map(f, s::AbstractString) requires f to return AbstractChar; try map(f, collect(s)) or a comprehension instead"))
        write(out, c′::AbstractChar)
    end
    String(take!(out))
end

function filter(f, s::AbstractString)
    out = IOBuffer(sizehint=sizeof(s))
    for c in s
        f(c) && write(out, c)
    end
    String(take!(out))
end

## string first and last ##

"""
    first(s::AbstractString, n::Integer)

Get a string consisting of the first `n` characters of `s`.

```jldoctest
julia> first("∀ϵ≠0: ϵ²>0", 0)
""

julia> first("∀ϵ≠0: ϵ²>0", 1)
"∀"

julia> first("∀ϵ≠0: ϵ²>0", 3)
"∀ϵ≠"
```
"""
first(s::AbstractString, n::Integer) = @inbounds s[1:min(end, nextind(s, 0, n))]

"""
    last(s::AbstractString, n::Integer)

Get a string consisting of the last `n` characters of `s`.

```jldoctest
julia> last("∀ϵ≠0: ϵ²>0", 0)
""

julia> last("∀ϵ≠0: ϵ²>0", 1)
"0"

julia> last("∀ϵ≠0: ϵ²>0", 3)
"²>0"
```
"""
last(s::AbstractString, n::Integer) = @inbounds s[max(1, prevind(s, ncodeunits(s)+1, n)):end]

"""
    reverseind(v, i)

Given an index `i` in [`reverse(v)`](@ref), return the corresponding index in
`v` so that `v[reverseind(v,i)] == reverse(v)[i]`. (This can be nontrivial in
cases where `v` contains non-ASCII characters.)

# Examples
```jldoctest
julia> r = reverse("Julia")
"ailuJ"

julia> for i in 1:length(r)
           print(r[reverseind("Julia", i)])
       end
Julia
```
"""
reverseind(s::AbstractString, i::Integer) = thisind(s, ncodeunits(s)-i+1)

"""
    repeat(s::AbstractString, r::Integer)

Repeat a string `r` times. This can be written as `s^r`.

See also: [`^`](@ref)

# Examples
```jldoctest
julia> repeat("ha", 3)
"hahaha"
```
"""
repeat(s::AbstractString, r::Integer) = repeat(String(s), r)

"""
    ^(s::Union{AbstractString,AbstractChar}, n::Integer)

Repeat a string or character `n` times. This can also be written as `repeat(s, n)`.

See also: [`repeat`](@ref)

# Examples
```jldoctest
julia> "Test "^3
"Test Test Test "
```
"""
(^)(s::Union{AbstractString,AbstractChar}, r::Integer) = repeat(s, r)

# reverse-order iteration for strings and indices thereof
iterate(r::Iterators.Reverse{<:AbstractString}, i=lastindex(r.itr)) = i < firstindex(r.itr) ? nothing : (r.itr[i], prevind(r.itr, i))
iterate(r::Iterators.Reverse{<:EachStringIndex}, i=lastindex(r.itr.s)) = i < firstindex(r.itr.s) ? nothing : (i, prevind(r.itr.s, i))

## code unit access ##

"""
    CodeUnits(s::AbstractString)

Wrap a string (without copying) in an immutable vector-like object that accesses the code units
of the string's representation.
"""
struct CodeUnits{T,S<:AbstractString} <: DenseVector{T}
    s::S
    CodeUnits(s::S) where {S<:AbstractString} = new{codeunit(s),S}(s)
end

length(s::CodeUnits) = ncodeunits(s.s)
sizeof(s::CodeUnits{T}) where {T} = ncodeunits(s.s) * sizeof(T)
size(s::CodeUnits) = (length(s),)
elsize(s::CodeUnits{T}) where {T} = sizeof(T)
@propagate_inbounds getindex(s::CodeUnits, i::Int) = codeunit(s.s, i)
IndexStyle(::Type{<:CodeUnits}) = IndexLinear()
iterate(s::CodeUnits, i=1) = (@_propagate_inbounds_meta; i == length(s)+1 ? nothing : (s[i], i+1))

write(io::IO, s::CodeUnits) = write(io, s.s)

unsafe_convert(::Type{Ptr{T}},    s::CodeUnits{T}) where {T} = unsafe_convert(Ptr{T}, s.s)
unsafe_convert(::Type{Ptr{Int8}}, s::CodeUnits{UInt8}) = unsafe_convert(Ptr{Int8}, s.s)

"""
    codeunits(s::AbstractString)

Obtain a vector-like object containing the code units of a string.
Returns a `CodeUnits` wrapper by default, but `codeunits` may optionally be defined
for new string types if necessary.
"""
codeunits(s::AbstractString) = CodeUnits(s)

max_stars_repo_name: taqtiqa-mark/julia

max_issues_repo_name: taqtiqa-mark/julia

max_forks_repo_name: taqtiqa-mark/julia


text:
export PoissonSparseLU,solve!,mul!
using SparseArrays
using LinearAlgebra

struct PoissonSparseLU
    spLxy::SparseMatrixCSC{Float64,Int64}
    fspLxy::Any

    function PoissonSparseLU(sp)
        tlu=@elapsed fsp=lu(sp)
        println("LU factorization time: ",tlu)
        new(sp,fsp),tlu
    end
end

function solve!(Xxy,Bxy,sp::PoissonSparseLU)
    ldiv!(view(Xxy,1:length(Xxy)),sp.fspLxy,view(Bxy,1:length(Bxy)))
    Xxy
end

function mul!(Bxy,sp::PoissonSparseLU,Xxy)
    LinearAlgebra.mul!(view(Bxy,1:length(Xxy)),sp.spLxy,view(Xxy,1:length(Xxy)))
end


text:
import Base.ht_keyindex

# Cf. https://github.com/JuliaStats/StatsBase.jl/issues/135
immutable UnitWeights <: AbstractVector{Int}
end
Base.getindex(w::UnitWeights, ::Integer...) = 1
Base.getindex(w::UnitWeights, ::AbstractVector) = w

# About the type inference limitation which prompts this workaround, see
# https://github.com/JuliaLang/julia/issues/10880
Base.@pure eltypes(T) = Tuple{map(eltype, T.parameters)...}
Base.@pure vectypes(T) = Tuple{map(U -> Vector{U}, T.parameters)...}

# Internal function needed for now so that n is inferred
function _freqtable{T<:Real}(x::Tuple,
                             skipmissing::Bool = false,
                             weights::AbstractVector{T} = UnitWeights(),
                             subset::Union{Void, AbstractVector{Int}, AbstractVector{Bool}} = nothing)
    if !isa(subset, Void)
        x = map(y -> y[subset], x)
        weights = weights[subset]
    end

    n = length(x)
    l = map(length, x)
    vtypes = eltypes(typeof(x))

    for i in 1:n
        if l[1] != l[i]
            error("arguments are not of the same length: $l")
        end
    end

    if !isa(weights, UnitWeights) && length(weights) != l[1]
        error("'weights' (length $(length(weights))) must be of the same length as vectors (length $(l[1]))")
    end

    d = Dict{vtypes, eltype(weights)}()

    for (i, el) in enumerate(zip(x...))
        index = ht_keyindex(d, el)

        if index > 0
            @inbounds d.vals[index] += weights[i]
        else
            @inbounds d[el] = weights[i]
        end
    end

    if skipmissing
        filter!((k, v) -> !any(ismissing, k), d)
    end

    keyvec = collect(keys(d))

    dimnames = Vector{Vector}(n)
    for i in 1:n
        s = Set{vtypes.parameters[i]}()
        for j in 1:length(keyvec)
            push!(s, keyvec[j][i])
        end

        # convert() is needed for Union{T, Missing}, which currently gives a Vector{Any}
        # which breaks inference of the return type
        dimnames[i] = convert(Vector{vtypes.parameters[i]}, unique(s))
        try
            sort!(dimnames[i])
        catch err
            err isa MethodError || rethrow(err)
        end
    end

    a = zeros(eltype(weights), map(length, dimnames)...)::Array{eltype(weights), n}
    na = NamedArray(a, tuple(dimnames...)::vectypes(vtypes), ntuple(i -> "Dim$i", n))

    for (k, v) in d
        na[k...] = v
    end

    na
end

freqtable{T<:Real}(x::AbstractVector...;
                   skipmissing::Bool = false,
                   weights::AbstractVector{T} = UnitWeights(),
                   subset::Union{Void, AbstractVector{Int}, AbstractVector{Bool}} = nothing) =
    _freqtable(x, skipmissing, weights, subset)

# Internal function needed for now so that n is inferred
function _freqtable{n}(x::NTuple{n, AbstractCategoricalVector}, skipmissing::Bool = false)
    len = map(length, x)
    miss = map(v -> eltype(v) >: Missing, x)
    lev = map(v -> eltype(v) >: Missing && !skipmissing ? [levels(v); missing] : levels(v), x)
    dims = map(length, lev)
    # First entry is for missing values (only correct and used if present)
    ord = map((v, d) -> Int[d; CategoricalArrays.order(v.pool)], x, dims)

	for i in 1:n
	    if len[1] != len[i]
	        error(string("arguments are not of the same length: ", tuple(len...)))
	    end
	end

    sizes = cumprod([dims...])
    a = zeros(Int, dims)
    missingpossible = any(miss)

    @inbounds for i in 1:len[1]
        ref = x[1].refs[i]
        el = ord[1][ref + 1]
        anymiss = missingpossible & (ref <= 0)

        for j in 2:n
            ref = x[j].refs[i]
            anymiss |= missingpossible & (ref <= 0)
            el += (ord[j][ref + 1] - 1) * sizes[j - 1]
        end

        if !(missingpossible && skipmissing && anymiss)
            a[el] += 1
        end
    end

    NamedArray(a, lev, ntuple(i -> "Dim$i", n))
end

freqtable(x::AbstractCategoricalVector...; skipmissing::Bool = false) = _freqtable(x, skipmissing)

function freqtable(d::AbstractDataFrame, x::Symbol...; args...)
    a = freqtable([d[y] for y in x]...; args...)
    setdimnames!(a, x)
    a
end


text:
function make_grids_for_chunklist_timeseries(clt::ACLT; oversample_factor::Real = 1 ) where { ACLT<:AbstractChunkListTimeseries }
    map(i->RvSpectMLBase.make_grid_for_chunk(clt,i,oversample_factor=oversample_factor), 1:num_chunks(clt))
end

# Wrapper code to deal with weird data structures
function pack_chunk_list_timeseries_to_matrix(timeseries::ACLT, chunk_grids::Union{AR,AAV}; alg::Symbol = :TemporalGP,
    oversample_factor::Real = 1, smooth_factor::Real=1, verbose::Bool = false ) where {
        ACLT<:AbstractChunkListTimeseries, RT<:AbstractRange, AR<:AbstractArray{RT,1}, T<:Real, AV<:AbstractVector{T}, AAV<:AbstractArray{AV,1} }
    @assert alg == :Linear || alg == :GP  || alg == :Sinc || alg == :TemporalGP # TODO: Eventually move to traits-based system?
    num_obs = length(timeseries)
    num_λ = sum(length.(chunk_grids))
    flux_matrix = Array{Float64,2}(undef,num_λ,num_obs)
    var_matrix = Array{Float64,2}(undef,num_λ,num_obs)
    λ_vec = Array{Float64,1}(undef,num_λ)
    chunk_map = Array{UnitRange{Int64}}(undef,length(chunk_grids))
    mean_flux = zeros(num_λ)
    mean_var = zeros(num_λ)
    dfluxdlnλ = zeros(num_λ)
    d2fluxdlnλ2 = zeros(num_λ)
    if (alg == :GP) && (maximum(length.(chunk_grids))>1024)
        @error "To use a GP with more than 1024 points in a chunk, use a more efficient GP package, e.g., :TemporalGP."
    end
    if alg == :Sinc    # Setup workspace for Sync.  TODO: Put into functor
        filter_size=23
        kaiserB=13
        sinc_filter = SincInterpolation.create_filter_curve(filter_size*21; filter_size=filter_size, kaiserB=kaiserB)
    end

    for t in 1:num_obs
       idx_start = 0
       for c in 1:length(chunk_grids)
           idx = (idx_start+1):(idx_start+length(chunk_grids[c]))
           if verbose && !issorted(timeseries.chunk_list[t].data[c].λ)
               flush(stdout);  println("t= ",t, " c= ",c," idx= ", idx, " size(flux)= ",size(flux_matrix), " issorted(λ)=", issorted(timeseries.chunk_list[t].data[c].λ))
           end
           if t == 1
               λ_vec[idx] .= chunk_grids[c]
               chunk_map[c] = idx
           end
           mean_flux_in_chunk = mean(timeseries.chunk_list[t].data[c].flux)
           if mean_flux_in_chunk > 0
               if alg == :Linear
                   LinearInterpolation.interp_chunk_to_grid_linear!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c])
               elseif alg == :Sinc
                   SincInterpolation.interp_chunk_to_grid_sinc!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c], Filter=sinc_filter)
               elseif alg == :GP
                   GPInterpolation.interp_chunk_to_grid_gp!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c])
               elseif alg == :TemporalGP
                   RvSpectML.TemporalGPInterpolation.interp_chunk_to_grid_gp_temporal!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c] , use_logx=true, use_logy=false, smooth_factor=smooth_factor)
               end
               mean_flux[idx] .+= flux_matrix[idx,t]
               mean_var[idx]  .+= var_matrix[idx,t]
               dfluxdlnλ[idx] .+= calc_dfluxdlnlambda(view(flux_matrix,idx,t),view(λ_vec,idx))
               d2fluxdlnλ2[idx] .+= calc_d2fluxdlnlambda2(view(flux_matrix,idx,t),view(λ_vec,idx))

               flux_matrix[idx,t] ./= mean_flux_in_chunk
               var_matrix[idx,t]  ./= mean_flux_in_chunk^2
           else  # no flux in chunk?!?
               @warn "No flux in chunk " * string(c) * " at time " * string(t) * "."
           end
           idx_start += length(chunk_grids[c])
       end # for c
   end # for t

    # TODO Think if can propagate uncertainties better
    var_matrix .*= sqrt(oversample_factor)

    #flush(stdout)
    #println("# mean flux pre-normalize = ", mean(mean_flux))
    idx_start = 0
    for c in 1:length(chunk_grids)
        idx = (idx_start+1):(idx_start+length(chunk_grids[c]))
        mean_flux_in_chunk = mean(mean_flux[idx])
        #println("Normalizing chunk ", c, " (",idx,") by ", mean_flux_in_chunk)
        #dmeanfluxdlnλ[idx] .= calc_dfluxdlnlambda(vec(sum(view(flux_matrix,idx,:),dims=2)),vec(sum(view(var_matrix,idx,:), dims=2)) ) ./mean_flux_in_chunk
        mean_flux[idx] ./= mean_flux_in_chunk
        mean_var[idx] ./= mean_flux_in_chunk^2
        #dfluxdlnλ[idx] ./= mean_flux_in_chunk
        #d2fluxdlnλ2[idx] ./= mean_flux_in_chunk
        dfluxdlnλ[idx] .= calc_dfluxdlnlambda(mean_flux[idx],view(λ_vec,idx))
        d2fluxdlnλ2[idx] .= calc_d2fluxdlnlambda2(mean_flux[idx],view(λ_vec,idx))

        #dmeanfluxdlnλ[idx] .= calc_dfluxdlnlambda(vec(sum(view(flux_matrix,idx,:),dims=2)),view(λ_vec,idx))  ./mean_flux_in_chunk
        idx_start += length(chunk_grids[c])
    end
    #println("# mean flux post-normalize = ", mean(mean_flux))
    #flush(stdout)
    #mean_flux /= num_obs
    mean_var .*= sqrt(oversample_factor)
    #dmeanfluxdlnλ ./= num_obs
    return ( matrix=SpectralTimeSeriesCommonWavelengths(λ_vec,flux_matrix,var_matrix,chunk_map, TheoreticalInstrument1D() ), mean_flux=mean_flux, mean_var=mean_var, deriv=dfluxdlnλ, deriv2=d2fluxdlnλ2 )
end

# Wrapper code to deal with weird data structures
function pack_shifted_chunk_list_timeseries_to_matrix(timeseries::ACLT, chunk_grids::Union{AR,AAV}; alg::Symbol = :TemporalGP,
    oversample_factor::Real = 1, smooth_factor::Real=1, remove_rv_est::Bool = true, verbose::Bool = false ) where {
        ACLT<:AbstractChunkListTimeseries, RT<:AbstractRange, AR<:AbstractArray{RT,1}, T<:Real, AV<:AbstractVector{T}, AAV<:AbstractArray{AV,1} }
    @assert alg == :Linear || alg == :GP  || alg == :Sinc || alg == :TemporalGP # TODO: Eventually move to traits-based system?
    num_obs = length(timeseries)
    num_λ = sum(length.(chunk_grids))
    flux_matrix = Array{Float64,2}(undef,num_λ,num_obs)
    var_matrix = Array{Float64,2}(undef,num_λ,num_obs)
    λ_vec = Array{Float64,1}(undef,num_λ)
    chunk_map = Array{UnitRange{Int64}}(undef,length(chunk_grids))
    mean_flux = zeros(num_λ)
    mean_var = zeros(num_λ)
    dfluxdlnλ = zeros(num_λ)
    d2fluxdlnλ2 = zeros(num_λ)
    if remove_rv_est
        @assert haskey(first(timeseries.metadata),:rv_est)
    end
    if (alg == :GP) && (maximum(length.(chunk_grids))>1024)
        @error "To use a GP with more than 1024 points in a chunk, use a more efficient GP package, e.g., :TemporalGP."
    end
    if alg == :Sinc    # Setup workspace for Sync.  TODO: Put into functor
        filter_size=23
        kaiserB=13
        sinc_filter = SincInterpolation.create_filter_curve(filter_size*21; filter_size=filter_size, kaiserB=kaiserB)
    end
    for t in 1:num_obs
       idx_start = 0
       boost_factor = remove_rv_est ? calc_doppler_factor(timeseries.metadata[t][:rv_est]) : 1.0
       for c in 1:length(chunk_grids)
           idx = (idx_start+1):(idx_start+length(chunk_grids[c]))
           if verbose
               flush(stdout);  println("t= ",t, " c= ",c," idx= ", idx, " size(flux)= ",size(flux_matrix))
           end
           if t == 1
               λ_vec[idx] .= chunk_grids[c]
               chunk_map[c] = idx
           end
           mean_flux_in_chunk = mean(timeseries.chunk_list[t].data[c].flux)
           if mean_flux_in_chunk > 0
               #println(" mean flux in chunk ", c, " at time ", t, " = ", mean_flux_in_chunk)
               if alg == :Linear
                   LinearInterpolation.interp_chunk_to_shifted_grid_linear!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c], boost_factor )
                   view(flux_matrix,idx,t) .= calc_running_mean_weighted(view(flux_matrix,idx,t))
                   view(var_matrix,idx,t) .= calc_running_mean_weighted(view(var_matrix,idx,t))
               elseif alg == :Sinc
                   SincInterpolation.interp_chunk_to_shifted_grid_sinc!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c], boost_factor, Filter=sinc_filter )
                   view(flux_matrix,idx,t) .= calc_running_mean_weighted(view(flux_matrix,idx,t))
                   view(var_matrix,idx,t) .= calc_running_mean_weighted(view(var_matrix,idx,t))
               elseif alg == :GP
                   GPInterpolation.interp_chunk_to_shifted_grid_gp!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c], boost_factor)
               elseif alg == :TemporalGP
                   #TemporalGPsInterpolation.
                   interp_chunk_to_shifted_grid_gp_temporal!(view(flux_matrix,idx,t), view(var_matrix,idx,t), timeseries.chunk_list[t].data[c], chunk_grids[c], boost_factor, use_logx=true, use_logy=false, smooth_factor=smooth_factor)
               end
               if t == 1
                  λ_vec[idx] .= chunk_grids[c]
                  chunk_map[c] = idx
               end
            else  # no flux in chunk?!?
                @warn "No flux in chunk " * string(c) * " at time " * string(t) * "."
            end

           mean_flux[idx] .+= view(flux_matrix,idx,t)
           mean_var[idx]  .+= view(var_matrix,idx,t)
           #dfluxdlnλ[idx] .+= calc_dfluxdlnlambda(view(flux_matrix,idx,t),view(λ_vec,idx))
           #d2fluxdlnλ2[idx] .+= calc_d2fluxdlnlambda2(view(flux_matrix,idx,t),view(λ_vec,idx))

           flux_matrix[idx,t] ./= mean_flux_in_chunk
           var_matrix[idx,t]  ./= mean_flux_in_chunk^2

           idx_start += length(chunk_grids[c])
       end
    end

    # TODO Think if can propagate uncertainties better
    var_matrix .*= sqrt(oversample_factor)

    #flush(stdout)
    #println("# mean flux pre-normalize = ", mean(mean_flux))
    idx_start = 0
    for c in 1:length(chunk_grids)
        idx = (idx_start+1):(idx_start+length(chunk_grids[c]))
        mean_flux_in_chunk = mean(mean_flux[idx])
        #println("Normalizing chunk ", c, " (",idx,") by ", mean_flux_in_chunk)
        #dmeanfluxdlnλ[idx] .= calc_dfluxdlnlambda(vec(sum(view(flux_matrix,idx,:),dims=2)),vec(sum(view(var_matrix,idx,:), dims=2)) ) ./mean_flux_in_chunk
        #sum_in_chunk = vec(sum(view(flux_matrix,idx,:),dims=2))
        mean_flux[idx] ./= mean_flux_in_chunk
        mean_var[idx] ./= mean_flux_in_chunk^2
        #dfluxdlnλ[idx] ./= mean_flux_in_chunk
        #d2fluxdlnλ2[idx] ./= mean_flux_in_chunk
        if length(idx) >= 2
            dfluxdlnλ[idx] .= calc_dfluxdlnlambda(mean_flux[idx],view(λ_vec,idx))
        else
            dfluxdlnλ[idx] .= 0.0
        end
        if length(idx) >= 3
            d2fluxdlnλ2[idx] .= calc_d2fluxdlnlambda2(mean_flux[idx],view(λ_vec,idx))
        else
            d2fluxdlnλ2[idx] .= 0.0
        end
        idx_start += length(chunk_grids[c])
    end

    #println("# mean flux post-normalize = ", mean(mean_flux))
    #flush(stdout)
    #mean_flux /= num_obs
    mean_var .*= sqrt(oversample_factor)
    #dmeanfluxdlnλ ./= num_obs
    return ( matrix=SpectralTimeSeriesCommonWavelengths(λ_vec,flux_matrix,var_matrix,chunk_map, TheoreticalInstrument1D() ), mean_flux=mean_flux, mean_var=mean_var, deriv=dfluxdlnλ, deriv2=d2fluxdlnλ2 )
end

"""   repack_flux_vector_to_chunk_matrix(λ, flux, var, chunk_map, λc; alg )
Warning:  This doesn't work yet
"""
function repack_flux_vector_to_chunk_matrix(λ::AA1, flux::AA2, var::AA3, chunk_map::CMT, λc::AA4; alg::Symbol = :TemporalGP,
    oversample_factor::Real = 1, smooth_factor::Real=1, verbose::Bool = false ) where {
        T1<:Real, AA1<:AbstractVector{T1}, T2<:Real, AA2<:AbstractVector{T2}, T3<:Real, AA3<:AbstractVector{T3}, T4<:Real, AA4<:AbstractVector{T4},
        CMT<:AbstractArray{UnitRange{Int64}} } #RT<:AbstractRange, AR<:AbstractArray{RT,1}, T<:Real, AV<:AbstractVector{T}, AAV<:AbstractArray{AV,1} }
    @assert alg == :Linear || alg == :GP  || alg == :Sinc || alg == :TemporalGP # TODO: Eventually move to traits-based system?
	@warn "This doesn't work yet"
    num_chunks = length(chunk_map)
    num_λ_per_chunk = minimum(length.(chunk_map))
    flux_matrix = zeros(num_λ_per_chunk,num_chunks)
    var_matrix = zeros(num_λ_per_chunk,num_chunks)
    if (alg == :GP) && (num_λ_per_chunk>1024)
        @error "To use a GP with more than 1024 points in a chunk, use a more efficient GP package, e.g., :TemporalGP."
    end
    if alg == :Sinc    # Setup workspace for Sync.  TODO: Put into functor
        filter_size=23
        kaiserB=13
        sinc_filter = SincInterpolation.create_filter_curve(filter_size*21; filter_size=filter_size, kaiserB=kaiserB)
    end

   idx_start = 0
   for c in 1:length(chunk_map)
       idx_in = (idx_start+1):(idx_start+length(chunk_map[c]))
       idx_out_begin = floor(Int,(idx_in[1]+idx_in[end])//2)-floor(Int,num_λ_per_chunk//2)
       idx_out_end = idx_out_begin+num_λ_per_chunk-1
       idx_out = idx_out_begin:idx_out_end
       Δλ_in_chunk = λ[idx_out].-λc[c]
       mean_flux_in_chunk = mean(flux[idx_out])
       chunk_of_spectrum = ChunkOfSpectrum{eltype(λ),eltype(flux),eltype(var)}(view(λ,idx_out), view(flux,idx_out), view(var,idx_out))
       Δv = 400  # TODO: Figure out good choice for Δv
       Δλ = Δv /RvSpectML.speed_of_light_mps*λc[c]
       if mean_flux_in_chunk > 0
           if alg == :Linear
               LinearInterpolation.interp_chunk_to_grid_linear!(view(flux_matrix,:,c), view(var_matrix,:,c), chunk_of_spectrum, chunk_grid)
               view(flux_matrix,idx,:,c) .= calc_running_mean_weighted(view(flux_matrix,idx,:,c))
               view(var_matrix,idx,:,c) .= calc_running_mean_weighted(view(var_matrix,idx,:,c))
           elseif alg == :Sinc
               chunk_grid = range(λc[c]-(num_λ_per_chunk-1)//2*Δλ, stop=λc[c]+(num_λ_per_chunk-1)//2*Δλ, length=num_λ_per_chunk )
               SincInterpolation.interp_chunk_to_grid_sinc!(view(flux_matrix,:,c), view(var_matrix,:,c), chunk_of_spectrum, chunk_grid, Filter=sinc_filter)
               view(flux_matrix,idx,:,c) .= calc_running_mean_weighted(view(flux_matrix,idx,:,c))
               view(var_matrix,idx,:,c) .= calc_running_mean_weighted(view(var_matrix,idx,:,c))
           elseif alg == :GP
               GPInterpolation.interp_chunk_to_grid_gp!(view(flux_matrix,:,c), view(var_matrix,:,c), chunk_of_spectrum, chunk_grid)
           elseif alg == :TemporalGP
               TemporalGPsInterpolation.interp_chunk_to_grid_gp_temporal!(view(flux_matrix,:,c), view(var_matrix,:,c), chunk_of_spectrum, chunk_grid , use_logx=true, use_logy=false, smooth_factor=smooth_factor)
           end
       else  # no flux in chunk?!?
           @warn "No flux in chunk " * string(c) * " at time " * string(t) * "."
       end
       idx_start += length(chunk_map[c])
   end # for c
   return flux_matrix
    #flush(stdout)
    #println("# mean flux pre-normalize = ", mean(mean_flux))
    #=
    idx_start = 0
    for c in 1:length(chunk_grids)
        idx = (idx_start+1):(idx_start+length(chunk_grids[c]))
        mean_flux_in_chunk = mean(mean_flux[idx])
        #println("Normalizing chunk ", c, " (",idx,") by ", mean_flux_in_chunk)
        #dmeanfluxdlnλ[idx] .= calc_dfluxdlnlambda(vec(sum(view(flux_matrix,idx,:),dims=2)),vec(sum(view(var_matrix,idx,:), dims=2)) ) ./mean_flux_in_chunk
        mean_flux[idx] ./= mean_flux_in_chunk
        mean_var[idx] ./= mean_flux_in_chunk^2
        #dfluxdlnλ[idx] ./= mean_flux_in_chunk
        #d2fluxdlnλ2[idx] ./= mean_flux_in_chunk
        dfluxdlnλ[idx] .= calc_dfluxdlnlambda(mean_flux[idx],view(λ_vec,idx))
        d2fluxdlnλ2[idx] .= calc_d2fluxdlnlambda2(mean_flux[idx],view(λ_vec,idx))

        #dmeanfluxdlnλ[idx] .= calc_dfluxdlnlambda(vec(sum(view(flux_matrix,idx,:),dims=2)),view(λ_vec,idx))  ./mean_flux_in_chunk
        idx_start += length(chunk_grids[c])
    end
    =#
    #println("# mean flux post-normalize = ", mean(mean_flux))
    #flush(stdout)
    #mean_flux /= num_obs
    #mean_var .*= sqrt(oversample_factor)
    #dmeanfluxdlnλ ./= num_obs
    #return flux ( matrix=SpectralTimeSeriesCommonWavelengths(λ_vec,flux_matrix,var_matrix,chunk_map, TheoreticalInstrument1D() ), mean_flux=mean_flux, mean_var=mean_var, deriv=dfluxdlnλ, deriv2=d2fluxdlnλ2 )
end


text:
using PowerPlots
using Test

using PowerModels
# import Ipopt

using Memento
using Memento.TestUtils

prev_level = getlevel(getlogger(PowerModels))
PowerModels.logger_config!("error") # silence PowerModels logger while testing
data = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")

@testset "PowerPlots.jl" begin

    @testset "VegaLite Backend" begin
        ## Test vegalite plotting
        case = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")

        # get the logger of PowerPlots
        logger = getlogger(PowerPlots)
        start_level = getlevel(logger)
        PowerPlots.logger_config!("info") # show everything

        # case5.m should not cause any warning messages besides this one
        function is_unexpected_warning(output::String)
            output != nothing && output != """Data column "ComponentType" does not exist for DC line"""
        end

        @test_nolog(logger, "warn", is_unexpected_warning, powerplot(case)) # vanilla function call should not cause any unexpected error messages
        @test_nolog(logger, "error", r".+", powerplot(case)) # vanilla function call should not cause any error messages

        # Color attribute tests
        @test_nolog(logger, "error", r".+", powerplot(case; bus_color=:red, branch_color=["blue", :green])) # valid function call should not cause any error messages
        @test_nolog(logger, "warn", is_unexpected_warning, powerplot(case; bus_color=:red, branch_color=["blue", :green])) # test valid colors
        @test_warn(logger, r"Ignoring unexpected attribute (.*)$", powerplot(case; fake_attr=nothing)) # test ignoring of invalid attributes
        @test_warn(logger, r"Color value for (.*) should be given as symbol or string$", powerplot(case; bus_color=0)) # test invalid color type warning
        @test_warn(logger, r"Invalid color (.*) given for (.*)$", powerplot(case; bus_color=:goosegray)) # test invalid CSS color warning
        @test_warn(logger, r"Invalid color (.*) given for (.*)$", powerplot(case; bus_color=[:red, 0])) # test array with invalid CSS colors

        # Numeric attribute tests
        @test_nolog(logger, "error", r".+", powerplot(case; width=100, height="100")) # valid function call should not cause any error messages
        @test_nolog(logger, "warn", is_unexpected_warning, powerplot(case; width=100, height="100")) # test valid numeric attribute
        @test_warn(logger, r"Invalid number (.*) given for (.*)$", powerplot(case; width="abcd")) # test invalid numeric string given for numeric attribute
        @test_warn(logger, r"Value for (.*) should be given as a number or numeric String$", powerplot(case; width=:zero)) # test invalid datatype given for numeric attribute

        # Data label tests
        @test_nolog(logger, "error", r".+", powerplot(case; bus_data=:ComponentType, gen_data="ComponentType",
            bus_data_type=:ordinal, gen_data_type="nominal")) # valid function call should not cause any error messages
        @test_nolog(logger, "warn", is_unexpected_warning, powerplot(case; bus_data=:ComponentType,
            gen_data="ComponentType", bus_data_type=:ordinal, gen_data_type="nominal")) # test valid function call
        @test_warn(logger, r"Value for (.*) should be given as a String or Symbol$", powerplot(case; bus_data=0)) # test invalid datatype passed into data label
        @test_warn(logger, r"Data column :blah does not exist for (.*)$", powerplot(case; bus_data=:blah)) # test invalid data column
        @test_warn(logger, r"Data type :blah not a valid VegaLite data type$", powerplot(case; bus_data_type=:blah)) # test invalid data type

        PowerPlots.logger_config!(start_level) # restore logger to initial level
    end

    @testset "PowerPlotsLogger" begin
        # get the logger of PowerPlots
        logger = getlogger(PowerPlots)

        start_level = getlevel(logger) # get initial log level
        # test logger configuration functions
        PowerPlots.logger_config!("info")
        @test getlevel(logger) == "info"
        PowerPlots.logger_config!("warn")
        @test getlevel(logger) == "warn"
        PowerPlots.silence()
        @test getlevel(logger) == "error"

        PowerPlots.logger_config!(start_level) # restore logger to initial level
    end

    @testset "PowerModelsDataFrame" begin
        case = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")
        df = PowerModelsDataFrame(data)
        @test typeof(df) <:PowerModelsDataFrame
        @test size(df.branch) == (7,20)

        data_mn = PowerModels.replicate(data,2)
        df_mn = PowerModelsDataFrame(data_mn)
        @test typeof(df_mn) <:PowerModelsDataFrame # these functions didn't error
        @test size(df_mn.branch) == (14,21)

        comp_df = comp_dict_to_dataframe(data["branch"])
        @test size(comp_df) == (7,19)
    end

    @testset "PowerModelsGraph and Layouts" begin
        case = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")
        PMG = PowerModelsGraph(case)
        positions = layout_graph_kamada_kawai!(PMG)
        @test size(positions) == (10,)
        @test typeof(positions) == Vector{Vector{Float64}}
    end

    @testset "Experimental" begin
        using PowerPlots.Experimental
        case = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")
        plot1 = powerplot(case)
        PowerPlots.Experimental.add_zoom!(plot1)
        @test true # what do I test here?

        case = PowerModels.parse_file("$(joinpath(dirname(pathof(PowerModels)), ".."))/test/data/matpower/case5.m")
        plot1 = powerplot(case)
        PowerPlots.Experimental.cartesian2geo!(plot1)
        @test true # what do I test here?

    end

end

PowerModels.logger_config!(prev_level); # reset PowerModels logger to previous level


text:
# Test Metrics from Distances.jl
m, n = 3, 5
sz_img = (3, 3)
sz_img_3 = (3, 3, 3)

dist_list = [SqEuclidean(),
            Euclidean(),
            Cityblock(),
            Minkowski(2),
            Minkowski(2.5),
            Hamming(),
            TotalVariation(),
            ZNCC(),
            SumAbsoluteDifference(),
            SumSquaredDifference(),
            MeanAbsoluteError(),
            MeanSquaredError(),
            RootMeanSquaredError(),
            ]

for dist in dist_list
    @testset "$(typeof(dist))" begin
        @info "test: $(typeof(dist))"
        # Gray image
        type_list = generate_test_types([Bool, Float32, N0f8], [Gray])
        A = [0.0 0.0 0.0; 1.0 1.0 1.0; 0.0 0.0 0.0]
        B = [1.0 1.0 1.0; 0.0 0.0 0.0; 1.0 1.0 1.0]
        for T in type_list
            @testset "$T" begin
                test_colwise(dist, n, sz_img, T)
                test_pairwise(dist, m, n, sz_img, T)
                test_ndarray(dist, sz_img_3, T)

                a = A .|> T
                b = B .|> T
                test_numeric(dist, a, b, T)
                test_numeric(dist, channelview(a), channelview(b), T)
            end
        end
        test_cross_type(dist, A, B, type_list)

        # RGB image
        type_list = generate_test_types([Float32, N0f8], [RGB])
        A = [RGB(0.0, 0.0, 0.0) RGB(0.0, 1.0, 0.0) RGB(0.0, 1.0, 1.0)
            RGB(0.0, 0.0, 1.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 1.0, 0.0)
            RGB(1.0, 1.0, 1.0) RGB(1.0, 0.0, 1.0) RGB(0.0, 0.0, 0.0)]
        B = [RGB(0.0, 0.0, 0.0) RGB(0.0, 0.0, 1.0) RGB(1.0, 1.0, 1.0)
            RGB(0.0, 1.0, 0.0) RGB(1.0, 0.0, 0.0) RGB(1.0, 0.0, 1.0)
            RGB(0.0, 1.0, 1.0) RGB(1.0, 1.0, 0.0) RGB(0.0, 0.0, 0.0)]
        for T in type_list
            @testset "$T" begin
                test_colwise(dist, n, sz_img, T)
                test_pairwise(dist, m, n, sz_img, T)
                test_ndarray(dist, sz_img_3, T)

                a = A .|> T
                b = B .|> T

                test_numeric(dist, a, b, T)
            end
        end
        test_cross_type(dist, A, B, type_list)
    end
end


text:
# check if vector is sorted and elements are unique
function _checksu(v::AbstractVector)
    isu = true
    for i in 2:length(v)
        if v[i] <= v[i-1]
            isu = false
            break
        end
    end
    isu
end

# a vector whose elements are sorted and unique
struct SUVector{T} <: AbstractVector{T}
    # elements are stored in a standard vector
    data::Vector{T}
    # constructor
    function SUVector(data::AbstractVector{T}) where T
        # check sorted and unique properties
        @assert _checksu(data)
        # NOTE: data is copied
        new{T}(copy(data))
    end
end

# empty constructor
SUVector{T}() where T = SUVector(T[])

# formatted printing
show(io::IO, v::SUVector) = show(io::IO, v.data)

# shape of the vector
size(v::SUVector) = size(v.data)

# indexed reading
getindex(v::SUVector, i::Int) = getindex(v.data, i)

# indexed writing
function setindex!(v::SUVector, vi, i::Int)
    # check to make sure order and uniqueness of elements is preserved
    l = length(v)
    if l > 1
        errmsg = "This doesn't preserve order or uniqueness of elements."
        if i == 1
            @assert vi < v[i+1] errmsg
        elseif i == l
            @assert v[i-1] < vi errmsg
        else
            @assert v[i-1] < vi < v[i+1] errmsg
        end
    end
    setindex!(v.data, vi, i)
end

# pushing an element onto the end of a sorted unique vector
function push!(v::SUVector{T}, el::T) where T
    if !isempty(v)
        @assert el > v[end] "Must be larger than current last element."
    end
    push!(v.data, el)
end

# appending elements onto the end of a sorted unique vector
function append!(v::SUVector{T}, els::AbstractVector{T}) where T
    if !isempty(v)
        errmsg = "Doesn't preserve order or uniqueness of elements."
        @assert _checksu(els) && els[1] > v[end] errmsg
    end
    append!(v.data, els)
end

# intersection of two sorted unique vectors
function intersect(v1::SUVector{T}, v2::SUVector{T}) where T
    s1 = length(v1)
    s2 = length(v2)
    # empty intersection to be added to
    int = SUVector{T}()
    # step through vectors simultaneously
    i1 = 1
    i2 = 1
    while i1 <= s1 && i2 <= s2
        if v1[i1] < v2[i2]
            # not an intersection, move towards one
            i1 += 1
        elseif v2[i2] < v1[i1]
            # not an intersection, move towards one
            i2 += 1
        else
            # an intersection
            push!(int, v1[i1])
            i1 += 1
            i2 += 1
        end
    end
    int
end

# symmetric difference of two sorted unique vectors
function symdiff(v1::SUVector{T}, v2::SUVector{T}) where T
    s1 = length(v1)
    s2 = length(v2)
    # empty symmetric difference to be added to
    sdiff = SUVector{T}()
    if s1 == 0
        # v1 is empty, so return all of v2
        append!(sdiff, v2)
    elseif s2 == 0
        # v2 is empty, so return all of v1
        append!(sdiff, v1)
    else
        # step through vectors simultaneously
        i1 = 1
        i2 = 1
        while i1 <= s1 && i2 <= s2
            if v1[i1] < v2[i2]
                # add to sdiff and move towards and intersection
                push!(sdiff, v1[i1])
                i1 += 1
            elseif v2[i2] < v1[i1]
                # add to sdiff and move towards and intersection
                push!(sdiff, v2[i2])
                i2 += 1
            else
                # intersection
                i1 += 1
                i2 += 1
            end
        end
        # once one vector is finished, add the rest of the other vector
        if i1 > s1 && i2 <= s2
            append!(sdiff, view(v2, i2:s2))
        elseif i2 > s2 && i1 <= s1
            append!(sdiff, view(v1, i1:s1))
        end
    end
    sdiff
end

# union of two sorted unique vectors
function union(v1::SUVector{T}, v2::SUVector{T}) where T
    s1 = length(v1)
    s2 = length(v2)
    # empty union to be added to
    uni = SUVector{T}()
    if s1 == 0
        # v1 is empty, return all of v2
        append!(uni, v2)
    elseif s2 == 0
        # v2 is empty, return all of v1
        append!(uni, v1)
    else
        # step through vectors simultaneously
        i1 = 1
        i2 = 1
        while i1 <= s1 && i2 <= s2
            if v1[i1] < v2[i2]
                # add to union and move towards an intersection
                push!(uni, v1[i1])
                i1 += 1
            elseif v2[i2] < v1[i1]
                # add to union and move towards and intersection
                push!(uni, v2[i2])
                i2 += 1
            else
                # intersection, only add one but increment both
                push!(uni, v1[i1])
                i1 += 1
                i2 += 1
            end
        end
        # when finished one of the vectors, add the rest of the other
        if i1 > s1 && i2 <= s2
            append!(uni, view(v2, i2:s2))
        elseif i2 > s2 && i1 <= s1
            append!(uni, view(v1, i1:s1))
        end
    end
    uni
end

max_stars_repo_name: aormorningstar/SBRG

max_issues_repo_name: aormorningstar/SBRG

max_forks_repo_name: aormorningstar/SBRG


text:
__precompile__(true)
module Setfield
using MacroTools
using MacroTools: isstructdef, splitstructdef, postwalk

include("lens.jl")
include("sugar.jl")
include("functionlenses.jl")
include("settable.jl")
include("experimental.jl")

# To correctly dispatch to `show(::IO, ::CustomLens)` when it is defined by a
# user, we avoid defining the generic `show(::IO, ::Lens)`.  This way, we can
# safely call `show` inside `ComposedLens` without worrying about the
# `StackOverflowError` that can be easily triggered in the previous approach.
# See also:
# * https://github.com/jw3126/Setfield.jl/pull/86
# * https://github.com/jw3126/Setfield.jl/pull/88
for n in names(Setfield, all=true)
    T = getproperty(Setfield, n)
    if T isa Type && T <: Lens && (T === ComposedLens || has_atlens_support(T))
        @eval Base.show(io::IO, l::$T) = _show(io, nothing, l)
    end
end

end