text:
# # Lagrange Polynomial Models

#src We generate a documentation page from this source file.
#src Let's tell the user:
#md # This file is automatically generated from source code.
#md # For usage examples refer to [Summary & Quick Examples](@ref lagrange_summary).

# ## Intro and Prerequisites

# Polyoniaml interpolation models are a common choice for surrogate modeling.
# In our setting we want to construct models for ``n``-variate objectives and
# use polynomials of degree 1 or 2.
# We hence need a basis for the space ``Π_n^d`` of polynomials.
# Given a point set that is suited for interpolation (a *poised* set) we can
# use the Lagarange basis ``\{l_i\}`` with ``l_i(x_j) = δ_{i,j}`` to easily
# find the coefficients for vector valued models.

# We use `DynamicPolynomials` for polynomial arithmetic, `NLopt` to optimize
# polynomials and some more packages:
using DynamicPolynomials
import NLopt
import Combinatorics

# ## Surrogate Interface Implementations

# The model itself is defined only by its vector of Lagrange basis polynomials
# and the coefficients.

@with_kw struct LagrangeModel{
        B <: AbstractArray{<:AbstractPolynomialLike},
        G <: AbstractArray{<:AbstractArray{<:AbstractPolynomialLike}},
        V <: AbstractVector{<:AbstractVector{<:AbstractFloat} } } <: SurrogateModel
    basis :: B
    grads :: G
    coeff :: V
    fully_linear :: Bool = false
end

fully_linear( lm :: LagrangeModel ) = lm.fully_linear

# There is a multitude of configuration parameters, most of which will
# be explained later:

@with_kw mutable struct LagrangeConfig <: SurrogateConfig

    "Degree of the surrogate model polynomials."
    degree :: Int = 2

    "Enlargement parameter to consider more points for inclusion."
    θ_enlarge :: Real = 2

    "Quality parameter in Λ-Poisedness Algorithm."
    LAMBDA :: Real = 1.5

    "Whether or not the interpolation sets must be Λ-poised (and the models fully linear)."
    allow_not_linear :: Bool = false

    "Whether or not to try to construct a new interpolation set in each iteration."
    optimized_sampling :: Bool = true

    ## if optimized_sampling = false, shall we try to use saved sites?
    save_path :: String = ""
    io_lock :: Union{Nothing, Threads.ReentrantLock} = nothing

    algo1_max_evals :: Int = -1
    algo2_max_evals :: Int = -1

    algo1_solver :: Symbol = :LN_BOBYQA
    algo2_solver :: Symbol = :LN_BOBYQA

    max_evals :: Int64 = typemax(Int64);

    @assert 1 <= degree <= 2 "Only linear and quadratic models are supported."
    @assert LAMBDA > 1 "`LAMBDA` must be > 1."
    @assert let algo_str = string( algo1_solver );
        length( algo_str ) > 2 && string(algo_str[2]) == "N"
    end "`algo1_solver` must be a derivative-free NLopt algorithm."
    @assert let algo_str = string( algo2_solver );
        length( algo_str ) > 2 && string(algo_str[2]) == "N"
    end "`algo2_solver` must be a derivative-free NLopt algorithm."
end

# Overwrite `lock` and `unlock` so we can use `nothing` as a "lock":
function Base.lock(::Nothing) end
function Base.unlock(::Nothing) end

# The required method implementations are straightforward.
# Note, thate we allow the models to be combined to vector functions if they
# share the same configuration to avoid redundant efforts whilst constructing models.

max_evals( cfg :: LagrangeConfig ) :: Int = cfg.max_evals
combinable( cfg :: LagrangeConfig ) :: Bool = true

# We also need to introduce our own implementation for `isequal` and `hash` for
# `LagrangeConfig`s to be combinable, see [the docs too](https://docs.julialang.org/en/v1/base/base/).
function Base.hash( cfg :: LagrangeConfig, h :: UInt )
	return hash( getfield.( cfg, Tuple( fn for fn ∈ fieldnames(LagrangeConfig) ) ), h )
end
function Base.isequal( cfg1 :: LagrangeConfig, cfg2 :: LagrangeConfig )
	all( isequal( getfield(cfg1, fn), getfield(cfg2, fn) ) for fn in fieldnames( LagrangeConfig) )
end

# The `LagrangeMeta` simply holds the database indices of the results
# we want to interpolate at.
# We also store the output indices of the model for convenience and carry
# polynomials that act on ``[0,1]^n``.

@with_kw struct LagrangeMeta{
        CB <: Union{Nothing, Vector{<:AbstractPolynomialLike}},
        LB <: Union{Nothing, Vector{<:AbstractPolynomialLike}},
        P <: Union{Nothing, AbstractVector{<:AbstractVector{<:Real}}}
    } <: SurrogateMeta
    interpolation_indices :: Vector{Int} = []
    out_indices :: Vector{Int} = []
    canonical_basis :: CB = nothing
    lagrange_basis :: LB = nothing
    stamp_points :: P = nothing     ## used only if `unoptimized_sampling == false`
    fully_linear :: Bool = false
end

get_saveable_type( T :: LagrangeMeta ) = LagrangeMeta{Nothing,Nothing}
get_saveable( meta :: LagrangeMeta ) = LagrangeMeta(;
    interpolation_indices = meta.interpolation_indices, out_indices = meta.output_indices )

export LagrangeConfig, LagrangeMeta, LagrangeModel

# ## Construction

# ### A Bit of Theory
# The canonical basis is obtained by calculating the non-negative
# integral solutions to the euqation
# ```math
# x_1 + … + x_n \le d.
# ```
# These solutions can be found using the `Combinatorics`
# package via `multiexponents(n,d)` (`d` must be successively increased).

function non_negative_ineq_solutions(deg, n_vars)
	Iterators.flatten( ( collect( Combinatorics.multiexponents( n_vars, d )) for d = 0 : deg ) )
end

function get_poly_basis( deg, n_vars)
	exponents = non_negative_ineq_solutions(deg, n_vars )
	polys = let
		@polyvar x[1:n_vars]
		[ prod(x.^e) for e in exponents ]
	end
	return polys
end

# We are going to use the canonical basis to
# determine a **poised set** of points.
# This does in fact work with any polynomial basis for ``Π_n^d``. \
# In the process of doing so, we also modify (a copy of?)
# the basis so that it becomes the Lagrange basis for the returned point set.
#
# The Larange basis is formed by normalizing and orthogonalizing with respect to the point set:

function orthogonalize_polys( poly_arr, x, i )
	## normalize i-th polynomial with respect to `x`
	p_i = poly_arr[i] / poly_arr[i](x)

	## orthogonalize
	return [ j != i ? poly_arr[j] - ( poly_arr[j](x) * p_i ) : p_i for j = eachindex( poly_arr ) ]
end

# We use Algorithm 6.2 and Algorithm 6.3 from the book
# "Introduction to Derivative-Free Optimization"
# by Conn et. al. \
# Algorithm 6.2 makes the set poised (suited for interpolation) and
# returns the corresponding Lagrange basis.
# Algorithm 6.3 takes the poised set and the Lagrange basis and
# tries to make it ``Λ``-poised. ``Λ`` must be greater 1 and
# a smaller value makes the set more suited for good models.

"""
    get_poised_set( basis, points; solver = :LN_BOBYQA, max_solver_evals = -1 )

Compute a point set suited for polynomial interpolation.

Input:
* `basis`: A vector of polynomials constituting a basis for the polynomial space.
* `points`: (optional) A set of candidate points to be tried for inclusion into the poised set.
* `solver`: NLopt solver to use. Should be derivative-free.
* `max_solver_evals`: Maximum number of evaluations in each optimization run.

Return:
* `poised_points :: Vector{T}` where `T` is either a `Vector{F}` or an `SVector{n_vars, F}` and `F` is the precision of the points in `points`, but at least `Float32`.
* `lagrange_basis :: Vector{<:AbstractPolynomialLike}`: The Lagrange basis corresponding to `poised_points`.
* `point_indices`: An array indicating which points from `points` are also in `poised_points`. A positive entry corresponds to the index of a poised point in `points`. If a poised point is new, then the entry is `-1`.
"""
function get_poised_set( basis, points :: AbstractArray{T} = Vector{MIN_PRECISION}[];
		solver = :LN_BOBYQA, max_solver_evals = -1 ) where {
		T <: AbstractArray{<:Real}
	}

	p = length(basis)
	@assert p > 0 "`basis` must not be an empty array."

    @logmsg loglevel3 "Trying to find a poised set with $(p) points."

	vars = variables( basis[end] )
	n_vars = length(vars)
	@assert n_vars > 0 "The number of variables must be positive."

	if max_solver_evals < 0
		max_solver_evals = 2000 * n_vars
	end

	F = promote_type( eltype( T ), MIN_PRECISION )
	#P_type = n_vars > 100 ? Vector{Vector{F}} : Vector{SVector{n_vars, F}}
    P_type = Vector{Vector{F}}
	ZERO_TOL = min(eps(F) * 100, eps(Float16) * 10)

	## indicates which points from points have been accepted
	point_indices = fill(-1, p)
	not_accepted_indices = collect( eachindex( points ) )
	## return array of points that form a poised set
	poised_points = P_type(undef, p)

	new_basis = basis
	for i = 1 : p
		_points = points[not_accepted_indices]

		## find the point that maximizes the i-th polynomial
		## if the polynomial is constant, then the first remaining point is used (j = 1)
		l_max, j = if isempty(_points)
			0.0, 0
		else
			findmax( abs.( [ new_basis[i]( x ) for x in _points ] ) )
		end

		if l_max > ZERO_TOL
			## accept the `j`-th point from `_points`
			poised_points[i] = _points[j]
			### indicate what the actual point index was
			point_indices[i] = not_accepted_indices[j]
			### delete from further consideration
			deleteat!(not_accepted_indices, j)
		else
			## no point was suitable to add to the set
			## trying to find the maximizer for a | l_i(x) |
			opt = NLopt.Opt( solver, n_vars )
			opt.lower_bounds = zeros(F, n_vars )
            opt.upper_bounds = ones(F, n_vars )
            opt.maxeval = max_solver_evals
            opt.xtol_rel = 1e-3
            opt.max_objective = (x,g) -> abs( new_basis[i](x) )

            ## try to find a good starting point
			x₀_tmp = [ rand(F, n_vars) for i = 1 : 50 * n_vars ]
            x₀ = x₀_tmp[argmax( abs.(new_basis[i].(x₀_tmp)) ) ]

			_, ξ, ret = NLopt.optimize(opt, x₀)

			poised_points[i] = ξ
		end

		new_basis = orthogonalize_polys( new_basis, poised_points[i], i )
	end

	return poised_points, new_basis, point_indices
end

"""
    make_set_lambda_poised( basis, points;
        LAMBDA = 1.5, solver = :LN_BOBYQA, max_solver_evals = -1, max_loops = -1, skip_indices = [1,] )

Make the output of `get_poised_set` even better suited for interpolation.

Input:
* `basis`: A vector of polynomials constituting a Lagrange basis for the polynomial space.
* `points`: The vector of points belonging to the Lagrange basis.
* `LAMBDA :: Real > 1`: Determines the quality of the interpolation.
* `solver`: NLopt solver to use. Should be derivative-free.
* `max_solver_evals`: Maximum number of evaluations in each optimization run.
* `max_loops`: Maximum number of loops that try to make the set Λ-poised.
* `skip_indices`: Inidices of points to discard last.

Return:
* `poised_points :: Vector{T}` where `T` is either a `Vector{F}` or an `SVector{n_vars, F}` and `F` is the precision of the points in `points`, but at least `Float32`.
* `lagrange_basis :: Vector{<:AbstractPolynomialLike}`: The Lagrange basis corresponding to `poised_points`.
* `point_indices`: An array indicating which points from `points` are also in `poised_points`. A positive entry corresponds to the index of a poised point in `points`. If a poised point is new, then the entry is `-1`.
"""
function make_set_lambda_poised( basis, points :: AbstractArray{T};
		LAMBDA :: Real = 1.5, solver = :LN_BOBYQA, max_solver_evals = -1,
		max_loops = -1, skip_indices = [1,] ) where {
		T <: AbstractArray{<:Real}
	}

	@assert length(basis) == length(points) "Polynomial array `basis` and point array `points` must have the same length."
	if length(points) > 0
		n_vars = length(points[1])
		@assert n_vars > 0 "The number of variables must be positive."

		F = promote_type( eltype( T ), MIN_PRECISION )
		#P_type = n_vars > 100 ? Vector{Vector{F}} : Vector{SVector{n_vars, F}}
        P_type = Vector{Vector{F}}

		if max_loops < 0
			max_loops = length(basis) * 100
		end

		if max_solver_evals < 0
			max_solver_evals = 2000 * n_vars
		end

       	@logmsg loglevel3 "Trying $(max_loops) times to make a set poised with Λ = $(LAMBDA)."

		new_basis = basis
		new_points = P_type(points)
		point_indices = collect(eachindex(new_points))

		for k = 1 : max_loops
            iₖ = -1
		    xₖ = points[1]
			for (i, polyᵢ) in enumerate(new_basis)
				opt = NLopt.Opt( solver, n_vars )
				opt.lower_bounds = zeros(F, n_vars)
				opt.upper_bounds = ones(F, n_vars)
				opt.maxeval = max_solver_evals
				opt.xtol_rel = 1e-3
				opt.max_objective = (x,g) -> abs( polyᵢ( x ) )

                x₀ = points[i]

				abs_lᵢ, xᵢ, _ = NLopt.optimize(opt, x₀)

				if abs_lᵢ > LAMBDA
					iₖ = i
					xₖ = xᵢ
					if iₖ ∉ skip_indices
						## i is not prioritized we can brake here
						break
					end#if
				end#if
			end#for

			if iₖ > 0
                @logmsg loglevel4 "Discarding point $(iₖ)."
				## perform a point swap
				new_points[iₖ] = xₖ
				point_indices[iₖ] = -1
				## adapt coefficients of lagrange basis
				new_basis = orthogonalize_polys( new_basis, xₖ, iₖ )
			else
				## we are done, the set is lambda poised
				break
			end#if
		end#for

		return new_points, new_basis, point_indices
	else
		return points, basis, collect(eachindex(points))
	end

end

# And a convenient function that combines both steps:

function get_lambda_poised_set( basis, points; solver1 = :LN_BOBYQA, solver2 = :LN_BOBYQA, max_solver_evals1 = -1, max_solver_evals2 = -1, LAMBDA = 1.5, max_lambda_loops = -1 )
	lagrange_points, lagrange_basis, lagrange_indices = get_poised_set(
		basis, points; solver = solver1, max_solver_evals = max_solver_evals1 )
	lambda_points, lambda_basis, lambda_indices = make_set_lambda_poised(
		lagrange_basis, lagrange_points; LAMBDA, max_loops = max_lambda_loops,
		solver = solver2, max_solver_evals = max_solver_evals2 )
	combined_indices = [ i < 0 ? i : lagrange_indices[j] for (j,i) in enumerate( lambda_indices ) ]
	return lambda_points, lambda_basis, combined_indices
end

# We actually only try to find points suitable points in the hypercube ``[0,1]^n``.
# The points can be (un)scaled with the usual methods.
# But for `Polynomial`s we can actually use substition to make evaluation more effective.

"Return vector of polynomials that unscales variables from [0,1]^n to [lb,ub]."
function get_unscaling_poly( vars, lb, ub )
    ## we don't have to check for Inf here because of finite trust region
    w = ub .- lb
    return vars .* w .+ lb
end

"Return vector of polynomials that scales variables from [lb, ub] to [0,1]^n."
function get_scaling_poly( vars, lb, ub )
    w = ub .- lb
    return ( vars .- lb ) ./ w
end

# ### Method Implementations

# We will use the functions from above in the `prepare_XXX` routines:\
# The initial `prepare_init_model` function should return a meta object that can be used
# to build an initial surrogate model.
# We delegate the work to `prepare_update_model`.
function prepare_init_model( cfg :: LagrangeConfig, objf :: AbstractObjective, mop :: AbstractMOP,
	id :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
	ensure_fully_linear = true, kwargs...)

    n_vars = num_vars( mop )

	meta = LagrangeMeta(;
        canonical_basis = get_poly_basis( cfg.degree, n_vars ),
        out_indices = output_indices(objf, mop)
    )
	return prepare_update_model(nothing, objf, meta, mop, id, db, ac; ensure_fully_linear, kwargs... )
end

# Usually, `prepare_update_model` would only accept a model as its first argument.
# Because of the trick from above, we actually allow `nothing`, too.

"""
    _consume_points(data_base, poised_points, poised_indices, candidate_indices)

Helper to return array of database indices for `poised_points` and
`poised_indices`. Add result to database if index is -1.
`candidate_indices` are the database indices of the points from the trust region.
"""
function _consume_points( db, poised_points, poised_indices, candidate_indices, lb, ub, F )
    interpolation_indices = Int[]
    for (i,ind) in enumerate(poised_indices)
        if ind < 0
            ## we need an additional new site
            new_db_id = new_result!(db, _unscale(poised_points[i], lb, ub), F[] )
            push!(interpolation_indices, new_db_id)
        else
            ## we could recycle a candidate point
            push!(interpolation_indices, candidate_indices[ind])
        end
    end
    return interpolation_indices
end

function _scale_poly_basis( poised_basis, lb, ub )
    ## we modify the basis so that the input is scaled to [0,1]^n with respect to
    ## the enlarged trust region bounds, because the poisedness algos sought points there
    poly_vars = variables( poised_basis[1] )
    scaling_poly = get_scaling_poly( poly_vars, lb, ub )

    zero_pol = sum( 0 .* poly_vars ) # TODO remove once https://github.com/JuliaAlgebra/DynamicPolynomials.jl/issues/92 is fixed

    return [ subs(p, poly_vars => scaling_poly) + zero_pol for p in poised_basis ]
end

function prepare_update_model( mod :: Union{Nothing, LagrangeModel}, objf :: AbstractObjective,
    meta :: LagrangeMeta,  mop :: AbstractMOP, iter_data :: AbstractIterData,
    db :: AbstractDB, algo_config :: AbstractConfig;
    ensure_fully_linear = true, kwargs... )

    x = get_x( iter_data )
    fx = get_fx( iter_data )
    F = eltype(fx)
    x_index = get_x_index( iter_data )
    n_vars = length(x)
    Δ = get_delta( iter_data )

    cfg = model_cfg(objf)
    lb, ub = local_bounds(mop, x, Δ * cfg.θ_enlarge )

    if cfg.optimized_sampling
        ## Find points in current trust region …
        candidate_indices = [x_index; results_in_box_indices( db, lb, ub, [x_index,] )]
        ## … and scale them to [0,1]^n
        candidate_points = [_scale(ξ, lb, ub) for ξ in get_site.(db, candidate_indices)]

        ## Get a poised set and lagrange basis
        poised_points, poised_basis, poised_indices = get_poised_set(
            meta.canonical_basis, candidate_points;
            solver = cfg.algo1_solver, max_solver_evals = cfg.algo1_max_evals
        )

        fully_linear = false
        ## Make set even better
        if ensure_fully_linear || !cfg.allow_not_linear
            ### We would like to keep x if possible
            skip_indices = let l = findfirst( i -> i == 1, poised_indices );
                isnothing(l) ? [] : [l,]
            end

            poised_points, poised_basis, indices_2 = make_set_lambda_poised(
                poised_basis, poised_points;
                LAMBDA = cfg.LAMBDA, solver = cfg.algo2_solver,
                max_solver_evals = cfg.algo2_max_evals, skip_indices
            )
            poised_indices = [ i < 0 ? i : poised_indices[j] for (j,i) in enumerate( indices_2 ) ]
            fully_linear = true
        end

        interpolation_indices = _consume_points( db, poised_points, poised_indices, candidate_indices, lb, ub, F)
        scaled_basis = _scale_poly_basis( poised_basis, lb, ub )

        return LagrangeMeta(;
            interpolation_indices,
            out_indices = meta.out_indices,
            canonical_basis = meta.canonical_basis,
            lagrange_basis = scaled_basis,
            fully_linear
        )

    else
        ## unoptimized sampling: we only look for a good point set once
        ## in the very first iteration and store the basis and the points
        ## in the meta data which is then passed through in subsequent iterations
        lpoints, lbasis = if isnothing(meta.lagrange_basis)
            candidate_points = [ fill(.5, n_vars) ]
            lpoints, lbasis, _ = get_lambda_poised_set(
                meta.canonical_basis, candidate_points;
                solver1 = cfg.algo1_solver, solver2 = cfg.algo2_solver,
                max_solver_evals1 = cfg.algo1_max_evals, max_solver_evals2 = cfg.algo2_max_evals,
                LAMBDA = cfg.LAMBDA )

            lpoints, _scale_poly_basis( lbasis, lb, ub )
        else
            meta.stamp_points, meta.lagrange_basis
        end

        candidate_indices = [x_index,]
        lindices = fill(-1, length(lpoints))

        ## check if x (scaled to [0,1] wrt trust region bounds) is center of `lpoints`
        #src TODO does using `≈` make problems for small trust region radii? `==` always fails
        x_s = _scale(x, lb, ub)
        x_in_points_index = findfirst(χ -> χ ≈ x_s, lpoints )
        if !isnothing(x_in_points_index)
             candidate_indices[ x_in_points_index ] = 1
        end

        interpolation_indices = _consume_points( db, lpoints, lindices, candidate_indices, lb, ub, F)

        return LagrangeMeta(;
            interpolation_indices,
            out_indices = meta.out_indices,
            lagrange_basis = lbasis,
            stamp_points = lpoints,
            fully_linear = true
        )
    end
end#function

# The improvement preparation enforces a Λ-poised set:
function prepare_improve_model( mod :: Union{Nothing, LagrangeModel}, objf :: AbstractObjective, meta :: LagrangeMeta,
    mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, algo_config :: AbstractConfig;
    kwargs... )
    return prepare_update_model( mod, objf, meta, mop, iter_data, db, algo_config; ensure_fully_linear = true, kwargs...)
end

# Now, in the 2-phase construction process, first all `prepare_` functions are called for all surrogate models.
# Then, the unevaluated results are evaluated and we can proceed with the model building.
# As before, `_init_model` simply delegates work to `update_model`. \
# Not much is left to do, only to retrieve the correct values from the database to use as
# coefficients.
# We also store the gradient (vector of polynomials) for each basis polynomial.

function _init_model( cfg :: LagrangeConfig, objf :: AbstractObjective, mop :: AbstractMOP,
	iter_data :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig, meta :: LagrangeMeta; kwargs... )
	return update_model( nothing, objf, meta, mop, iter_data, db, ac; kwargs... )
end

function update_model( mod::Union{Nothing,LagrangeModel}, objf:: AbstractObjective,
    meta :: LagrangeMeta, mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
	kwargs... )

    coeff = [ c[ meta.out_indices ] for c in get_value.(db, meta.interpolation_indices) ]

    return LagrangeModel(;
        coeff, fully_linear = meta.fully_linear,
        basis = copy(meta.lagrange_basis),
        ## NOTE I don't know why I need to copy here
        ## but if i don't copy then testing fails:
        ## the meta data does hold a valid Lagrange basis but the model does not !?
        grads = [ differentiate( p, variables(p) ) for p in meta.lagrange_basis ]
    ), meta
end

function improve_model( mod::Union{Nothing,LagrangeModel}, objf:: AbstractObjective,
    meta :: LagrangeMeta, mop :: AbstractMOP, iter_data :: AbstractIterData, db :: AbstractDB, ac :: AbstractConfig;
	kwargs... )
    return update_model( mod, objf, meta, mop, iter_data, db, algo_config; kwargs...)
end

# ## Evaluation
# The evaluation of some output is
# ```math
# \sum_{i=1}^p c_i l_i( x ),
# ```
# where ``p = \dim Π_n^d``.

function _eval_poly_vec( poly_vec, x )
    [ p(x) for p in poly_vec ]
end

function eval_models( lm :: LagrangeModel, x̂ :: Vec, ℓ :: Int)
    return sum( c[ℓ] * p(x̂) for (c,p) in zip( lm.coeff, lm.basis ) )
end

function eval_models( lm :: LagrangeModel, x̂ :: Vec )
    return sum( c * p(x̂) for (c,p) in zip( lm.coeff, lm.basis ) )
end

function get_gradient( lm :: LagrangeModel, x̂ :: Vec, ℓ :: Int )
    sum( c[ℓ] * _eval_poly_vec(p,x̂) for (c,p) in zip( lm.coeff, lm.grads ) )
end

function get_jacobian( lm :: LagrangeModel, x̂ :: Vec )
    grad_evals = [ _eval_poly_vec(p,x̂) for p in lm.grads ]
    no_out = length(lm.coeff[1])
    return transpose( hcat( (sum( c[ℓ] * g for (c,g) in zip( lm.coeff, grad_evals) ) for ℓ = 1 : no_out)... ) )
end

# ## [Summary & Quick Examples](@id lagrange_summary)

# 1. To use the default configuration for a scalar objective `f` do
#    ```julia
#    add_objective!(mop, f, LagrangeConfig())
#    ```
# 2. For a vector valued objective do
#    ```julia
#    add_vector_objective!(mop, f, LagrangeConfig(); n_out = 2)
#    ```
# 3. If you want a linear polyonmial only:
#    ```julia
#    add_objective!(mop, f, LagrangeConfig(;degree=1))
#    ```
# 4. By default, a new interpolation set is built in every iteration.
#    To use a "stamp" instead, turn of optimized sampling:
#    ```julia
#    add_objective!(mop, f, LagrangeConfig(;optimized_sampling=true))
#    ```

# ### Complete usage example
# ```julia
# using Morbit
# Morbit.print_all_logs()
# mop = MixedMOP(3)
#
# F = x -> [ sum( ( x .- 1 ).^2 ); sum( ( x .+ 1 ).^2 ) ]
#
# add_vector_objective!( mop, F, LagrangeConfig() )
#
# x_fin, f_fin, _ = optimize( mop, [-π, ℯ, 0])
# ```


text:
using JLD
using Plots
using SignalAnalysis
using Plots

filename = string("../JLD/PopulationScatter_adexp.jld")#, py"target_num_spikes")#,py"specimen_id)

if isfile(filename)
    opt_vec = load(filename, "opt_vec")
    extremum_param = load(filename,"extremum_param")
    vmgtt = load(filename, "vmgtt")
    vmgtv = load(filename, "vmgtv")
end

println("data loaded")

function vecplot(v)
    y = hcat(v...)'
    x = 1:length(v)
    #t = [i/1000.0 for i in x]
    plot(x, y, leg = :none,
    xaxis=("t", extrema(x)),
    yaxis=(string("V_{m}"), extrema(y)))
end
p1 = vecplot(opt_vec)
p1|>display
p2 = vecplot(vmgtv[:])
p2|>display
s_a = signal(opt_vec, length(opt_vec)/last(vmgtt))
s_b = signal(vmgtv, length(vmgtt)/last(vmgtt))
p = plot(s_a)
p2 = plot!(p, s_b)
display(plot(s_a))
display(plot(s_b))
savefig(p2, "aligned_VM.png")

ext: jl


text:
module CodecBzip2

export
    Bzip2Compressor,
    Bzip2CompressorStream,
    Bzip2Decompressor,
    Bzip2DecompressorStream

import TranscodingStreams:
    TranscodingStreams,
    TranscodingStream,
    Memory,
    Error,
    initialize,
    finalize,
    splitkwargs

import Bzip2_jll:
    libbzip2

include("libbzip2.jl")
include("compression.jl")
include("decompression.jl")

end # module


text:
using TheDataMustFlow
using MacroTools


text:
# Autogenerated wrapper script for RadeonProRender_jll for x86_64-w64-mingw32-cxx03
export libNorthstar64, libProRenderGLTF, libRadeonProRender64, libRprLoadStore64, libTahoe64

using CompilerSupportLibraries_jll
JLLWrappers.@generate_wrapper_header("RadeonProRender")
JLLWrappers.@declare_library_product(libNorthstar64, "Northstar64.dll")
JLLWrappers.@declare_library_product(libProRenderGLTF, "ProRenderGLTF.dll")
JLLWrappers.@declare_library_product(libRadeonProRender64, "RadeonProRender64.dll")
JLLWrappers.@declare_library_product(libRprLoadStore64, "RprLoadStore64.dll")
JLLWrappers.@declare_library_product(libTahoe64, "Tahoe64.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll)
    JLLWrappers.@init_library_product(
        libNorthstar64,
        "bin\\Northstar64.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libProRenderGLTF,
        "bin\\ProRenderGLTF.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libRadeonProRender64,
        "bin\\RadeonProRender64.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libRprLoadStore64,
        "bin\\RprLoadStore64.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libTahoe64,
        "bin\\Tahoe64.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
#export ARCSpectral, ARCSpectral_abs, ARCLDLt, ARCqKOp
export ALL_solvers

# Valid combinations
#

ALL_solvers = Function[]

include("ARCSpectral.jl")      #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ARCSpectral_abs.jl") #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRSpectral.jl") #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRSpectral_abs.jl") #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ARCLDLt.jl")         #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ARCLDLt_abs.jl")      #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRLDLt.jl")             #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRLDLt_abs.jl")         #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ARCqKOp.jl")			   #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ARCqKsparse.jl")          #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ARCqKdense.jl")          #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRKOp.jl")              #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRKsparse.jl")          #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("TRKdense.jl")           #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ST_TROp.jl")            #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ST_TRsparse.jl")        #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ST_TRdense.jl")        #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ST_ARCOp.jl")            #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ST_ARCsparse.jl")        #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

include("ST_ARCdense.jl")         #works in 0.7/1.1
push!(ALL_solvers,eval(fname))

## Needs further work to understand the HSL methods
## Will update them in the future
# include("ARCMA97.jl")
# push!(ALL_solvers,eval(fname))
#
# include("ARCMA97_abs.jl")
# push!(ALL_solvers,eval(fname))
#
# include("TRMA97.jl")
# push!(ALL_solvers,eval(fname))
#
# include("TRMA97_abs.jl")
# push!(ALL_solvers,eval(fname))
#
#
#
#
# include("ARCMA57.jl")
# push!(ALL_solvers,eval(fname))
#
# include("ARCMA57_abs.jl")
# push!(ALL_solvers,eval(fname))
#
# include("TRMA57.jl")
# push!(ALL_solvers,eval(fname))
#
# include("TRMA57_abs.jl")
# push!(ALL_solvers,eval(fname))


text:
if basename(pwd()) == "aoc"
    cd("2017/7")
end

function mermaid(filename::AbstractString)
    open(replace(filename, ".txt" => ".mmd"), "w") do f
        write(f, "graph TD\r\n")
        for line in eachline(filename; keep = true)
            if occursin(r"->", line)
                write(f, replace(line, ", " => " & ", "->" => "-->", r"\(\d+\)" => ""))
            end
        end
    end
end

struct Disc
    value::Int
    towers::Vector
end
Disc(value) = Disc(value, Disc[])

function loaddata(filename::AbstractString)
    lines = readlines(filename)
    discs = Dict{String,Disc}()
    for line in lines
        name, num = match(r"^(\S+) \((\d+)\)", line)
        push!(discs, name => Disc(parse(Int, num)))
    end
    for line in lines
        if occursin(" -> ", line)
            name, rhs = split(line, r" \(\d+\) -> ")
            othernames = split(rhs, ", ")
            push!(discs[name].towers, (discs[n] for n in othernames)...)
        end
    end
    discs
end

function findroot(discs)
    only(setdiff(values(discs), (tower for disc in values(discs) for tower in disc.towers)))
end

part1(filename::AbstractString) = part1(loaddata(filename))
function part1(discs)
    root = findroot(discs)
    for (k, v) in pairs(discs)
        if v == root
            return k
        end
    end
end

part2(filename::AbstractString) = part2(loaddata(filename))
function part2(discs)
    disc = findunbalanced(findroot(discs))
    weights = weight.(disc.towers)
    if weights[1] == weights[2]
        weights[2]
    else
        weights[3]
    end
end

function isbalanced(discs)
    weights = weight.(discs.towers)
    weights[1] == weights[2] == weights[3]
end

function findunbalanced(disc)
    if !isbalanced(disc)
        disc
    elseif !isisempty(disc.towers)
        something((findunbalanced(d) for d in disc.towers)..., Some(nothing))
    end
end

function weight(disc)
    disc.value + sum(weight(tower) for tower in disc.towers; init = 0)
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
include("../src/sir.jl")
using PyPlot
K = 50
Np = 2500
σ_o = 0.1
σ_d = 0.0
Δ = 5
Nth = 500
obsfun = obs1

x_trj, w_trj, y, x_true = assimilate(K, Np,
							σ_o, σ_d,
							Δ, Nth,
							obsfun)
x_trj = permutedims(x_trj,[3,2,1])
trj_len = K*Δ
y = permutedims(y,[2,1])
w_trj = permutedims(w_trj,[2,1])
x_true = permutedims(x_true,[2,1])

function plot_trj()
	#n_t = 10
	#trj_len = size(x_true)[1] - n_t + 1
	time_arr = LinRange(0, trj_len -1, trj_len)
	time_arr = dt*time_arr

	for i = 1:d
		fig, ax = subplots()
		ax.grid(true)
		ax.plot(time_arr, x_trj[:,:,i],"b",lw=0.9)
		ax.plot(time_arr, x_trj[:,end,i],"b",lw=0.9,label="\$ x^{($i)}\$ orbit of particle")
		if i==d
		    ax.plot(time_arr[1:Δ:end], y[:,1],"k",lw=3.0,label="Observed \$ x^{($i)}\$")
		end
		ax.legend(fontsize=28)
		ax.xaxis.set_tick_params(labelsize=28)
		ax.yaxis.set_tick_params(labelsize=28)
		ax.set_xlabel("Time",fontsize=28)
		ax.set_title("# particles = $(Np), partial observation", fontsize=28)
	end
end
function plot_rmse()
    filter_mean = reshape(sum(x_trj,dims=2)/Np, trj_len, d)
    norm2(x,dim) = sum(x -> x^2, x, dims=dim)
    rmse = norm2(x_true .- filter_mean, 2)
	rmse .= sqrt.(rmse)/d
	time_arr = LinRange(0, trj_len -1, trj_len)
	time_arr = dt*time_arr
	fig, ax = subplots()
	ax.grid(true)
	ax.plot(time_arr, rmse,"r",lw=3.0,label="RMSE")
	ax.legend(fontsize=28)
	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
	ax.set_xlabel("Time",fontsize=28)
	ax.set_title("# particles = $(Np), partial observation", fontsize=28)
end
function plot_weight_stats()
	fig, ax = subplots()
	ax.grid(true)
	ax.plot(w_trj,"b",lw=0.2)
	ax.plot(w_trj[:,end],"b",lw=1.0,label="particle weight orbit")
	ax.legend(fontsize=28)
	ax.xaxis.set_tick_params(labelsize=28)
	ax.yaxis.set_tick_params(labelsize=28)
end


text:
module CloneCmd

using TOML
using GitHub
using ..Options
using ..IonBase: read_github_auth
using ..SearchCmd: search_exact_package, fetch_repo_from_url
using Comonicon.Tools: prompt, cmd_error

function clone_url(url::String)
    if endswith(url, "jl.git")
        _clone(url, basename(url)[1:end-7])
    else
        _clone(url, basename(url))
    end
end

function clone_package(package::String, ion::Options.Ion=Options.read())
    info = search_exact_package(package, "", ion)
    isnothing(info) && cmd_error("cannot find $package in registries")
    uuid, reg, pkginfo = info
    pkg_toml = TOML.parsefile(joinpath(reg.path, pkginfo["path"], "Package.toml"))
    _clone(pkg_toml["repo"], pkg_toml["name"])
end

function _clone(url::String, to::String)
    username = readchomp(`git config user.name`)
    auth = GitHub.authenticate(read_github_auth())
    rp = fetch_repo_from_url(url; auth=auth)

    local has_access
    try
        has_access = iscollaborator(rp, username; auth=auth)
    catch e
        has_access = false
    end

    if has_access
        git_clone(url, to)
    elseif prompt("do not have access to $url, fork?")
        @info "fork upstream repo: $rp"
        owned_repo = fork_repo(rp, auth)
        git_clone(owned_repo.clone_url.uri, to)
        @info "setting upstream to $url"
        set_upstream(url, to)
    end
end

function fork_repo(repo, auth)
    return create_fork(repo; auth=auth)
end

function set_upstream(url::String, to::String)
    cd(joinpath(pwd(), to)) do
        run(`git remote add upstream $url`)
        run(`git fetch upstream`)
        run(`git branch --set-upstream-to=upstream/master`)
    end
end

function git_clone(url, to)
    run(`git clone $url $to`)
end

# Copied from IsURL.jl
# Source: https://github.com/sindresorhus/is-absolute-url (MIT license)
const windowsregex = r"^[a-zA-Z]:[\\]"
const urlregex = r"^[a-zA-Z][a-zA-Z\d+\-.]*:"

"""
    isurl(str)
Checks if the given string is an absolute URL.
# Examples
```julia-repl
julia> isurl("https://julialang.org")
true
julia> isurl("mailto:someone@example.com")
true
julia> isurl("/foo/bar")
false
```
"""
function isurl(str::AbstractString)
    return !occursin(windowsregex, str) && occursin(urlregex, str)
end

function isgithub(url::AbstractString)
    if startswith(url, "http") && startswith(split(url, "//")[2], "github.com")
        return true
    else
        return false
    end
end

function clone(url_or_package::String)
    if isurl(url_or_package)
        clone_url(url_or_package)
    else
        clone_package(url_or_package)
    end
    return
end

end

"""
clone a package or url and setup the local repo. If the current local
git user do not have push access to remote github repo, it will fork
this repo then clone the repo and set the input url as upstream.

# Arguments

- `url_or_package`: name of the package or url.
"""
@cast function clone(url_or_package::String)
    CloneCmd.clone(url_or_package)
end


text:
# Utils
function module_dir()
    return joinpath(@__DIR__, "..")
end

# Activate package
using Pkg
Pkg.activate(module_dir())

using LinearAlgebra
using LyceumBase
using LyceumMuJoCo
using LyceumMuJoCoViz
using MuJoCo
using Plots
using JLD2
using Random

include("methods.jl")

################################################################################
# Demo
################################################################################

timestep=0.01
jm, jd, mjsim = mj_model("dzhanibekov.xml", timestep=timestep)
tsim = 1.0
N = Int(floor(tsim/timestep))

# jd.qpos
jd.qvel .= [0,0,0,20π,0,0.01]

traj = zeros(20, 0)
E = zeros(N)
Plin = zeros(N)
Pang = zeros(N)
tcompute = 0.0
for i = 1:N
    jd.ctrl .= 0.0
    tcompute += @elapsed mj_step(jm, jd);
    traj = hcat(traj, deepcopy(getstate(mjsim)))
    E[i] = energy(jm, jd)
	Plin[i] = norm(linear_momentum(jm, jd))
	# Pang[i] = angular_momentum(jm, jd)
end
tsim / tcompute
plot(E)
plot(Plin)
# plot(Pang)

LyceumMuJoCoViz.visualize(mjsim, trajectories=traj)

jd.qpos
jd.qvel

timestep=0.01
jm, jd, mjsim = mj_model("dzhanibekov.xml", timestep=timestep)
jd.qvel .= [0,0,0,20π,0,0]
# jd.qvel .= [1,0,0,0.0,0,0]

mj_step(jm, jd)
MuJoCo.mj_subtreeVel(jm, jd)
jd.subtree_linvel
jd.subtree_angmom

jm.body_mass
mjsim.d.xipos

function linear_momentum(jm, jd)
	MuJoCo.mj_subtreeVel(jm, jd)
    p = Vector(sum(jd.subtree_linvel, dims=2)[:,1])
    return p
end

#
# ################################################################################
# ################################################################################
#
# timestep=0.01
# jm, jd, mjsim = mj_model("half_cheetah.xml", timestep=timestep)
# jd.qpos[2] += 2.0
# jd.qvel .= [0,0,0, 0,0,0, 0,0,5]
# # mj_step(jm, jd)
# MuJoCo.mj_subtreeVel(jm, jd)
# jd.subtree_linvel
# jd.subtree_angmom[:,1:4]
# jd.subtree_angmom[:,5:8]
# traj = zeros(28, 0)
#
#
#
# #
# # jm, jd, mjsim = mj_model("astronaut.xml", timestep=timestep)
# # tsim = 0.10
# # N = Int(floor(tsim/timestep))
#
#
# # jd.qpos
# # jd.qvel .= [1; zeros(26)]
#
# # traj = zeros(83, 0)
# E = zeros(N)
# Plin = zeros(N)
# Pang = zeros(N)
# tcompute = 0.0
# for i = 1:N
#     jd.ctrl .= 0.0
#     tcompute += @elapsed mj_step(jm, jd);
#     traj = hcat(traj, deepcopy(getstate(mjsim)))
#     E[i] = energy(jm, jd)
# 	Plin[i] = norm(linear_momentum(jm, jd))
# 	# Pang[i] = angular_momentum(jm, jd)
# end
# tsim / tcompute
# plot(E)
# plot(Plin)
# # plot(Pang)
#
# LyceumMuJoCoViz.visualize(mjsim, trajectories=traj)
#
#
#
# MuJoCo.mj_subtreeVel(jm, jd)
# jd.subtree_linvel
# jd.subtree_angmom
#
#


text:
@testset "rounding" begin
    st_johns = TimeZone("America/St_Johns")     # UTC-3:30 (or UTC-2:30)
    eucla = TimeZone("Australia/Eucla")         # UTC+8:45
    colombo = TimeZone("Asia/Colombo")          # See note below

    #=
    On 1996-05-25 at 00:00, the Asia/Colombo time zone in Sri Lanka moved from Indian
    Standard Time (UTC+5:30) to Lanka Time (UTC+6:30). On 1996-10-26 at 00:30, Lanka
    Time was revised from UTC+6:30 to UTC+6:00, marking a -00:30 transition. Transitions
    like these are doubly unusual (compared to the more common DST transitions) as it is both
    an half-hour transition and a transition that lands at midnight (causing
    1996-10-26T00:00 to be ambiguous; midnights are rarely ambiguous). In 2006,
    Asia/Colombo returned to Indian Standard Time, causing another -00:30 transition
    from 00:30 to 00:00.
    =#

    periods = (Year, Month, Week, Day, Hour, Minute, Second)

    @testset "basic" begin
        dt = DateTime(2016, 2, 5, 13, 10, 20, 500)
        for tz in [winnipeg, st_johns, eucla, colombo]
            lzdt = LaxZonedDateTime(dt, tz)

            @test floor(lzdt, Year) == LaxZonedDateTime(DateTime(2016), tz)
            @test floor(lzdt, Month) == LaxZonedDateTime(DateTime(2016, 2), tz)
            @test floor(lzdt, Week) == LaxZonedDateTime(DateTime(2016, 2), tz)
            @test floor(lzdt, Day) == LaxZonedDateTime(DateTime(2016, 2, 5), tz)
            @test floor(lzdt, Hour) == LaxZonedDateTime(DateTime(2016, 2, 5, 13), tz)
            @test floor(lzdt, Minute) == LaxZonedDateTime(DateTime(2016, 2, 5, 13, 10), tz)
            @test floor(lzdt, Second) == LaxZonedDateTime(DateTime(2016, 2, 5, 13, 10, 20), tz)

            @test ceil(lzdt, Year) == LaxZonedDateTime(DateTime(2017), tz)
            @test ceil(lzdt, Month) == LaxZonedDateTime(DateTime(2016, 3), tz)
            @test ceil(lzdt, Week) == LaxZonedDateTime(DateTime(2016, 2, 8), tz)
            @test ceil(lzdt, Day) == LaxZonedDateTime(DateTime(2016, 2, 6), tz)
            @test ceil(lzdt, Hour) == LaxZonedDateTime(DateTime(2016, 2, 5, 14), tz)
            @test ceil(lzdt, Minute) == LaxZonedDateTime(DateTime(2016, 2, 5, 13, 11), tz)
            @test ceil(lzdt, Second) == LaxZonedDateTime(DateTime(2016, 2, 5, 13, 10, 21), tz)

            @test round(lzdt, Year) == LaxZonedDateTime(DateTime(2016), tz)
            @test round(lzdt, Month) == LaxZonedDateTime(DateTime(2016, 2), tz)
            @test round(lzdt, Week) == LaxZonedDateTime(DateTime(2016, 2, 8), tz)
            @test round(lzdt, Day) == LaxZonedDateTime(DateTime(2016, 2, 6), tz)
            @test round(lzdt, Hour) == LaxZonedDateTime(DateTime(2016, 2, 5, 13), tz)
            @test round(lzdt, Minute) == LaxZonedDateTime(DateTime(2016, 2, 5, 13, 10), tz)
            @test round(lzdt, Second) == LaxZonedDateTime(DateTime(2016, 2, 5, 13, 10, 21), tz)
        end
    end

    @testset "transition" begin
        lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 1, 45), winnipeg)
        @test floor(lzdt, Day) == LaxZonedDateTime(DateTime(2016, 3, 13), winnipeg)
        @test floor(lzdt, Hour) == LaxZonedDateTime(DateTime(2016, 3, 13, 1), winnipeg)
        @test ceil(lzdt, Day) == LaxZonedDateTime(DateTime(2016, 3, 14), winnipeg)
        @test ceil(lzdt, Hour) == LaxZonedDateTime(DateTime(2016, 3, 13, 3), winnipeg)
        @test round(lzdt, Day) == LaxZonedDateTime(DateTime(2016, 3, 13), winnipeg)
        @test round(lzdt, Hour) == LaxZonedDateTime(DateTime(2016, 3, 13, 3), winnipeg)

        lzdt = LaxZonedDateTime(DateTime(1996, 10, 26, 0, 45), colombo)
        @test floor(lzdt, Hour) == LaxZonedDateTime(ZonedDateTime(1996, 10, 26, colombo, 2))
        @test isambiguous(floor(lzdt, Day))

        lzdt = LaxZonedDateTime(DateTime(1996, 10, 25, 23, 45), colombo)
        @test ceil(lzdt, Hour) == LaxZonedDateTime(ZonedDateTime(1996, 10, 26, colombo, 1))
        @test isambiguous(ceil(lzdt, Day))

        # Using `round` when either `ceil` or `floor` would result in an ambiguous or
        # missing date results in something unrepresentable (because `round` can't
        # figure out whether to use the result from `floor` or `ceil`).
        lzdt = LaxZonedDateTime(DateTime(1996, 10, 25, 23, 45), colombo)
        @test round(lzdt, Hour) == LaxZonedDateTime(ZonedDateTime(1996, 10, 26, colombo, 1))
        @test !isrepresentable(round(lzdt, Day))
    end

    @testset "non-existent" begin
        dne_lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 2, 0, 0, 1), winnipeg)
        @test !isvalid(dne_lzdt)
        @test !isambiguous(dne_lzdt)
        @test isnonexistent(dne_lzdt)
        @test isrepresentable(dne_lzdt)

        @testset "$(nameof(P))" for P in periods, f in (floor, ceil, round)
            result = f(dne_lzdt, P)

            if P <: DatePeriod
                @test isrepresentable(result)
            else
                @test !isrepresentable(result)
            end
        end
    end

    @testset "ambiguous" begin
        amb_lzdt = LaxZonedDateTime(DateTime(2015, 11, 1, 1, 0, 0, 1), winnipeg)
        @test !isvalid(amb_lzdt)
        @test isambiguous(amb_lzdt)
        @test !isnonexistent(amb_lzdt)
        @test isrepresentable(amb_lzdt)

        @testset "$(nameof(P))" for P in periods, f in (floor, ceil, round)
            result = f(amb_lzdt, P)

            if P <: DatePeriod
                @test isrepresentable(result)
            else
                @test !isrepresentable(result)
            end
        end
    end

    @testset "non-existent rounded" begin
        dne_lzdt = LaxZonedDateTime(DateTime(2016, 3, 13, 2), winnipeg)
        @test isnonexistent(dne_lzdt)

        @test floor(dne_lzdt, Hour) == dne_lzdt
        @test ceil(dne_lzdt, Hour) == dne_lzdt
        @test round(dne_lzdt, Hour) == dne_lzdt
    end

    @testset "ambiguous rounded" begin
        amb_lzdt = LaxZonedDateTime(DateTime(2015, 11, 1, 1), winnipeg)
        @test isambiguous(amb_lzdt)

        @test floor(amb_lzdt, Hour) == amb_lzdt
        @test ceil(amb_lzdt, Hour) == amb_lzdt
        @test round(amb_lzdt, Hour) == amb_lzdt
    end
end


text:
# Get the cumulative k-path length
function cumdists(kvs::AVec{<:AVec{<:Real}})
    N     = length(kvs)
    D     = length(first(kvs))
    dists = Vector{Float64}(undef, N)

    dists[1] = 0.0
    @inbounds for i in 2:N
        s = 0.0
        for j in 1:D # use a loop to avoid allocations
            s += abs2(kvs[i][j] - kvs[i-1][j])
        end
        dists[i] = dists[i-1] + sqrt(s)
    end
    return dists
end

# ---------------------------------------------------------------------------------------- #
# INTERPOLATIONS OF ::KPath

"""
$(TYPEDEF)
$(TYPEDFIELDS)
"""
struct KPathInterpolant{D} <: AbstractPath{SVector{D, Float64}}
    kpaths  :: Vector{Vector{SVector{D, Float64}}}
    labels  :: Vector{Dict{Int, Symbol}}
    basis   :: SVector{D, SVector{D, Float64}}
    setting :: Ref{BasisEnum}
end

size(kpi::KPathInterpolant) = (sum(length, kpi.kpaths),)
Base.@propagate_inbounds function composite_index(kpi::KPathInterpolant, i::Int)
    # index into the `i`th point in the "flattened" `kpi.kpaths`
    @boundscheck i < 0 && throw(BoundsError(kpi, i))
    j = 1
    i′ = i′′ = i
    while (i′′ -= length(kpi.kpaths[j])) > 0
        i′ = i′′
        j += 1
        @boundscheck j > length(kpi.kpaths) && throw(BoundsError(kpi, i))
    end
    @boundscheck i′ > length(kpi.kpaths[j]) && throw(BoundsError(kpi, i))
    return j, i′
end
Base.@propagate_inbounds function getindex(kpi::KPathInterpolant, i::Int)
    j, i′ = composite_index(kpi, i)
    return @inbounds kpi.kpaths[j][i′]
end
IndexStyle(::Type{<:KPathInterpolant}) = IndexLinear()
Base.@propagate_inbounds function setindex!(kpi::KPathInterpolant{D},
                                            val::SVector{D, Float64}, i::Int) where D
    j, i′ = composite_index(kpi, i)
    @inbounds kpi.kpaths[j][i′] = val
    return kpi
end

"""
    cartesianize!(kpi::KPathInterpolant)

Transform an interpolated **k**-path `kpi` in a lattice basis to a Cartesian basis with
(primitive) reciprocal lattice vectors `basis(kpi)`.
Modifies `kpi` in-place.
"""
function cartesianize!(kpi::KPathInterpolant)
    setting(kpi) === CARTESIAN && return kp
    for i in eachindex(kpi)
        @inbounds kpi[i] = cartesianize(kpi[i], kpi.basis)
    end
    set_setting!(kpi, CARTESIAN)
    return kpi
end

"""
    latticize!(kpi::KPathInterpolant{D})

Transform an interpolated **k**-path `kpi` in a Cartesian basis to a lattice basis with
(primitive) reciprocal lattice vectors `basis(kpi)`.
Modifies `kpi` in-place.
"""
function latticize!(kpi::KPathInterpolant)
    setting(kpi) === LATTICE && return kp
    basismatrix = hcat(kpi.basis...)
    for i in eachindex(kpi)
        @inbounds kpi[i] = latticize(kpi[i], basismatrix)
    end
    set_setting!(kpi, LATTICE)
    return kpi
end

"""
    interpolate(kp::KPath, length::Integer)
    interpolate(kp::KPath; length::Integer, density::Real) --> KPathInterpolant

Return an interpolant of `kp` with approximately `length` points distributed approximately
equidistantly across the full **k**-path (equidistance is measured in a Cartesian metric).

Note that the interpolant may contain fewer or more points than `length` (typically fewer).
`length` can also be provided as a keyword argument.

As an alternative to specifying the desired total number of interpolate points via `length`,
a desired density per unit (reciprocal) length can be specified via the keyword argument
`density`.

See also [`interpolate(::AbstractVector{::AbstractVector{<:Real}}, ::Integer)`](@ref).
"""
function interpolate(kp::KPath{D}, N::Integer) where D
    kpᶜ = setting(kp) === CARTESIAN ? kp : cartesianize(kp)
    distss = map(paths(kpᶜ)) do path
        map(1:length(path)-1) do i
            norm(points(kpᶜ)[path[i]] - points(kpᶜ)[path[i+1]])
        end
    end
    totaldist = sum(dists->sum(dists), distss)

    kipaths = [Vector{SVector{D, Float64}}() for _ in 1:length(paths(kp))]
    labels = [Dict{Int, Symbol}() for _ in 1:length(paths(kp))]
    for (j, (path, dists)) in enumerate(zip(paths(kp), distss))
        push!(labels[j], 1 => first(path))
        for i in 1:length(path)-1
            Nᵢ = max(convert(Int, div(N*dists[i], totaldist, RoundUp)), 2)
            append!(kipaths[j], range(points(kp)[path[i]], points(kp)[path[i+1]];
                                      length=Nᵢ)[1:end-1])
            push!(labels[j], length(kipaths[j])+1 => path[i+1])
        end
        push!(kipaths[j], points(kp)[last(path)])
    end

    return KPathInterpolant(kipaths, labels, basis(kp), Ref(setting(kp)))
end

function interpolate(kp::KPath;
            length::Union{Integer, Nothing}=nothing,
            density::Union{Real, Nothing}=nothing)
    if length === nothing && density === nothing
        throw(ArgumentError("must set either `length` or `density` keyword arguments explicitly or use positional argument for `length`"))
    elseif length !== nothing && density !== nothing
        throw(ArgumentError("cannot define `length` and `density` simultaneously"))
    elseif length !== nothing && density === nothing
        return interpolate(kp, length)
    elseif length === nothing && density !== nothing
        return interpolate_via_density(kp, density)
    end

    error("unreachable reached; please file an issue")
end

function interpolate_via_density(kp::KPath{D}, density::Real) where D
    kpᶜ = setting(kp) === CARTESIAN ? kp : cartesianize(kp)

    kipaths = [Vector{SVector{D, Float64}}() for _ in 1:length(paths(kp))]
    labels = [Dict{Int, Symbol}() for _ in 1:length(paths(kp))]
    for (j, path) in enumerate(paths(kpᶜ)) # over connected paths
        push!(labels[j], 1 => first(path))
        for i in 1:length(path)-1 # over linear segments
            klabᵢ   = path[i]
            klabᵢ₊₁ = path[i+1]

            # compute required number of points for current segment
            kᶜᵢ   = points(kpᶜ)[klabᵢ]   # start
            kᶜᵢ₊₁ = points(kpᶜ)[klabᵢ₊₁] # end
            distᵢ = norm(kᶜᵢ₊₁ - kᶜᵢ)
            Nᵢ    = max(ceil(Int, distᵢ*density), 2)

            # compute interpolation of current segment
            kᵢ    = points(kp)[klabᵢ]
            kᵢ₊₁  = points(kp)[klabᵢ₊₁]
            ks    = range(kᵢ, kᵢ₊₁, length=Nᵢ)[1:end]

            append!(kipaths[j], ks)
            push!(labels[j], length(kipaths[j]) => klabᵢ₊₁)
        end
    end

    return KPathInterpolant(kipaths, labels, basis(kp), Ref(setting(kp)))
end

"""
    splice(kp::KPath, N::Integer) --> KPathInterpolant

Return an interpolant of `kp` with `N` points inserted into each **k**-path segment of `kp`.

See also [`splice(::AbstractVector{::AbstractVector{<:Real}}, ::Integer)`](@ref).
"""
function splice(kp::KPath{D}, N::Integer) where D
    kipaths = [Vector{SVector{D, Float64}}() for _ in 1:length(paths(kp))]
    labels = [Dict{Int, Symbol}() for _ in 1:length(paths(kp))]
    for (j, path) in enumerate(paths(kp))
        push!(labels[j], 1 => first(path))
        for i in 1:length(path)-1
            append!(kipaths[j], (@view range(points(kp)[path[i]],
                                             points(kp)[path[i+1]],
                                             length=N+2)[1:end-1]))
            push!(labels[j], length(kipaths[j])+1 => path[i+1])
        end
        push!(kipaths[j], points(kp)[last(path)])
    end

    return KPathInterpolant(kipaths, labels, kp.basis, Ref(setting(kp)))
end

# ---------------------------------------------------------------------------------------- #
# INTERPOLATIONS OF ::AVec{<:AVec{<:Real}}

# TODO: Merge the implementations of `interpolate` and `splice` for `AVec{<:AVec{<:Real}}`
#       vs. `KPath`.

"""
    interpolate(kvs::AbstractVector{<:AbstractVector{<:Real}}, N::Integer)
        --> Vector{<:Vector{<:Real}}

Return an interpolated **k**-path between discrete **k**-points in `kvs`, with approximately
`N` interpolation points in total (typically fewer).

Note that, in general, it is not possible to do this so that all interpolated **k**-points
are equidistant; samples are however *exactly* equidistant across each linear segment
defined by points in `kvs` and *approximately* equidistant across all segments.

See also [`interpolate(::KPath, ::Integer)`](@ref) and [`splice`](@ref).
"""
function interpolate(kvs::AVec{<:AVec{<:Real}}, N::Integer)
    Nkpairs = length(kvs)-1
    elT     = typeof(first(first(kvs))/1)
    dists   = Vector{elT}(undef, Nkpairs)
    @inbounds for i in 1:Nkpairs
        dists[i] = norm(kvs[i] .- kvs[i+1])
    end
    totaldist  = sum(dists)
    N_per_dist = N/totaldist

    # TODO: preallocate
    kvpath = [float.(kvs[1])]
    @inbounds for i in 1:Nkpairs
        # try to maintain an even distribution of k-points along path
        Nᵢ = round(Int64, dists[i]*N_per_dist, RoundUp) # points in current segment
        Nᵢ = max(Nᵢ, 2) # must be at least two points
        new_kvs = range(kvs[i],kvs[i+1],length=Nᵢ)
        append!(kvpath, (@view new_kvs[2:end]))         # append `new_kvs` to `kvpath`
    end
    return kvpath
end

"""
    splice(kvs::AbstractVector{<:AbstractVector{<:Real}}, N::Integer)
        --> Vector{<:Vector{<:Real}}

Return an interpolated **k**-path between the discrete **k**-points in `kvs`, with `N`
interpolation points inserted in each segment defined by pairs of adjacent **k**-points.

See also [`splice(::KPath, ::Integer)`](@ref) and [`interpolate`](@ref).
"""
function splice(kvs::AVec{<:AVec{<:Real}}, N::Integer)
    D       = length(first(kvs))
    Nkpairs = length(kvs)-1
    N⁺²     = N+2

    elT    = typeof(first(first(kvs))/1)
    kvpath = [Vector{elT}(undef, D) for _ in Base.OneTo(Nkpairs+1 + Nkpairs*N)]
    kvpath[1] = kvs[1]
    start, stop = 2, N⁺²
    @inbounds for i in 1:Nkpairs
        new_kvs = range(kvs[i],kvs[i+1],length=N⁺²)
        @views kvpath[start:stop] .= new_kvs[2:end] # insert `new_kvs` in `kvpath`
        start  = stop+1
        stop  += N⁺²-1
    end
    return kvpath
end


text:
#using FactCheck #depracated, so removed
using QSWalk
using Test
using SparseArrays
using LinearAlgebra

include("utils.jl")
include("dirac.jl")
include("operator.jl")
include("demoralization.jl")
include("evolution.jl")


text:
using Documenter

deploydocs(
    repo = "github.com/tkf/PyBase.jl.git",
    julia = "1.0",
    target = "build",
    deps = nothing,
    make = nothing,
)


text:
module TupleBenchmarks

include(joinpath(dirname(@__FILE__), "..", "utils", "RandUtils.jl"))

using .RandUtils
using BenchmarkTools

const SUITE = BenchmarkGroup()

###############
# issue #5274 #
###############

struct TupleWrapper{N, T}
    data::NTuple{N, T}
end

Base.eltype(::TupleWrapper{N,T}) where {N,T} = T
Base.length(::TupleWrapper{N,T}) where {N,T} = N

function get_index(n::NTuple, i::Int)
    @inbounds v = n[i]
    return v
end

function get_index(n::TupleWrapper, i::Int)
    @inbounds v = n.data[i]
    return v
end

function sum_tuple(n::Union{NTuple{N, T}, TupleWrapper{N, T}}) where {N, T}
    s = zero(T)
    for i in 1:N
        s += get_index(n, i)
    end
    return s
end

const TUPLE_SUM_SIZES = (3, 8, 30, 60)
const TUPLE_SUM_TYPES = (Float32, Float64)

g = addgroup!(SUITE, "index", ["sum"])

for s in TUPLE_SUM_SIZES, T in TUPLE_SUM_TYPES
    tup = tuple(samerand(T, s)...)
    tupwrap = TupleWrapper(tup)
    g["sumelt", "NTuple", s, T] = @benchmarkable sum_tuple($tup) time_tolerance=0.40
    g["sumelt", "TupleWrapper", s, T] = @benchmarkable sum_tuple($tupwrap) time_tolerance=0.40
end

#####################
# Fixed Size Arrays #
#####################

# Short fixed size array implementation

abstract type FixedArray{T, N} <: AbstractArray{T, N} end

Base.IndexStyle(::Type{<: FixedArray}) = IndexLinear()
Base.getindex(fsa::FixedArray, i::Int) = fsa.data[i]

struct FixedVector{L, T} <: FixedArray{T, 1}
    data::NTuple{L, T}
end

Base.size(::FixedVector{L}) where {L} = (L,)
Base.size(::Type{FixedVector{L, T}}) where {L, T} = (L,)
Base.length(::FixedVector{L}) where {L} = L

struct FixedMatrix{R, C, T, RC} <: FixedArray{T, 2}
    data::NTuple{RC, T}
end

Base.size(::FixedMatrix{R, C}) where {R, C} = (R, C)
Base.size(::Type{FixedMatrix{R, C, T, RC}}) where {R, C, T, RC} = (R, C)
Base.length(::FixedMatrix{R, C, T, RC}) where {R, C, T, RC} = RC

# Reductions

@inline function perf_reduce(op, a::FixedArray)
    if length(a) == 1
        return a[1]
    else
        s = op(a[1], a[2])
        for j = 3:length(a)
            s = op(s, a[j])
        end
        return s
    end
end

perf_minimum(a::FixedArray) = perf_reduce(min, a)

@inline function perf_reduce(op, v0, a::FixedArray)
    if length(a) == 0
        return v0
    else
        s = v0
        @inbounds @simd for j = 1:length(a)
            s = op(s, a[j])
        end
        return s
    end
end

perf_sum(v::FixedArray{T}) where {T} = perf_reduce(+, zero(T), v)

@inline function perf_mapreduce(f, op, v0, a1::FixedArray)
    if length(a1) == 0
        return v0
    else
        s = op(v0, f(a1[1]))
        for j = 2:length(a1)
            s = op(s, f(a1[j]))
        end
        return s
    end
end

perf_sumabs2(a::FixedArray{T}) where {T} = perf_mapreduce(abs2, +, zero(T), a)

# Linear Algebra

@generated function perf_matvec(A::FixedMatrix{R, C, T}, b::FixedVector{C, T}) where {R, C, T}
    sA = size(A)
    indA = LinearIndices(sA)
    exprs = Expr(:tuple, [reduce((ex1,ex2) -> :(+($ex1,$ex2)),
                [:(A[$(indA[k, j])]*b[$j]) for j = 1:sA[2]]) for k = 1:sA[1]]...)
    return quote
        @inbounds return FixedVector{R, T}($exprs)
    end
end

@generated function perf_matmat(A::FixedMatrix{R1, C, T}, B::FixedMatrix{C, R2, T}) where {R1, R2, C, T}
    sA, sB = size(A), size(B)
    indA, indB = LinearIndices(sA), LinearIndices(sB)
    exprs =  Expr(:tuple, [reduce((ex1,ex2) -> :(+($ex1,$ex2)),
                [:(A[$(indA[k1, j])] * B[$(indB[j, k2])]) for j = 1:sA[2]]) for k1 = 1:sA[1], k2 = 1:sB[2]]...)
    result_type = FixedMatrix{R1, R2, T, (R1 * R2)}
    return quote
        @inbounds return $result_type($exprs)
    end
end

# Benchmarks #
##############
v2, v4, v8, v16 = [FixedVector((rand(i)...,)) for i in (2, 4, 8, 16)]
m2x2, m4x4, m8x8, m16x16 = [FixedMatrix{i,i, Float64, i*i}((rand(i*i)...,)) for i in (2, 4, 8, 16)]

# Reductions
g = addgroup!(SUITE, "reduction", ["tuple"])

for mv in (v2, v4, v8, v16, m2x2, m4x4, m8x8, m16x16)
    g["sum", size(mv)] = @benchmarkable perf_sum($mv)
    g["sumabs", size(mv)] = @benchmarkable perf_sumabs2($mv)
    g["minimum", size(mv)] = @benchmarkable perf_minimum($mv)
end

# Linear algebra

g = addgroup!(SUITE, "linear algebra", ["tuple"])

for (m, v) in zip((m2x2, m4x4, m8x8, m16x16), (v2, v4, v8, v16 ))
    g["matvec", size(m), size(v)] = @benchmarkable perf_matvec($m, $v)
    g["matmat", size(m), size(m)] = @benchmarkable perf_matmat($m, $m)
end

function _add( a::NTuple{4,Float32}, b::NTuple{4,Float32} )
    (a[1]+b[1],a[2]+b[2],a[3]+b[3],a[4]+b[4])
end

function _mul( a::NTuple{4,Float32}, b::NTuple{4,Float32} )
   (a[1]*b[1],a[2]*b[2],a[3]*b[3],a[4]*b[4])
end

function _madd( a::NTuple{4,Float32}, b::NTuple{4,Float32}, c::NTuple{4,Float32} )
   _add(_mul(a,b),c)
end

perf_tuple_11899(t) = _madd(t, t, t)
g = addgroup!(SUITE, "misc", ["tuple"])
t = (Float32(1.0), Float32(2.0), Float32(3.0), Float32(4.0))
g["11899"] = @benchmarkable perf_tuple_11899($t)

end # module


text:
const e  = 2.71828182845904523536
const pi = 3.14159265358979323846

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Arb"
version = v"2.18.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/fredrik-johansson/arb.git", "137f3b5e13e31163634e91f8547aaae7c1ad8c32")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
cd arb/

if [[ ${target} == *musl* ]]; then
   export CFLAGS=-D_GNU_SOURCE=1
elif [[ ${target} == *mingw* ]]; then
   # /lib is hardcoded in many places
   sed -i -e "s#/lib\>#/$(basename ${libdir})#g" configure
   # MSYS_NT-6.3 is not detected as MINGW
   extraflags=--build=MINGW${nbits}
fi

./configure --prefix=$prefix --disable-static --enable-shared --with-gmp=$prefix --with-mpfr=$prefix --with-flint=$prefix ${extraflags}
make -j${nproc}
make install LIBDIR=$(basename ${libdir})
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line

platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct("libarb", :libarb)
]

# Dependencies that must be installed before this package can be built
dependencies = [
    Dependency(PackageSpec(name="FLINT_jll", uuid="e134572f-a0d5-539d-bddf-3cad8db41a82"))
    Dependency(PackageSpec(name="GMP_jll", uuid="781609d7-10c4-51f6-84f2-b8444358ff6d"))
    Dependency(PackageSpec(name="MPFR_jll", uuid="3a97d323-0669-5f0c-9066-3539efd106a3"))
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: jbrea/Yggdrasil

max_issues_repo_name: jbrea/Yggdrasil

max_forks_repo_name: jbrea/Yggdrasil


text:
#========================
    NSE scripts
========================#
"""
    Element

A script or table elemental block
"""
Base.@kwdef mutable struct Element <:Leaf
    key::String     = ""
    value::String   = ""
end
function StructTypes.construct(::Type{T}, value::AbstractString) where T<:Element
    return Element(value = string(value))
end
function StructTypes.construct(::Type{Vector{T}}, value::AbstractString) where T<:Element
    return Vector{T}([Element(value = string(value))])
end

"""
    Table

Set of Elements and subtables related to a script's execution

All fields can be empty
"""
Base.@kwdef mutable struct Table <:Marsh
   key::String               = ""
   tables::Vector{Table}     = Vector{Table}()
   elements::Vector{Element} = Vector{Element}()
end
StructTypes.names(::Type{Table}) = ((:elements, :elem), (:tables, :table))

"""
    Script
"""
Base.@kwdef mutable struct Script <:Marsh
    id::String                  = ""
    output::String              = ""
    elements::Vector{Element}   = Vector{Element}()
    tables::Vector{Table}       = Vector{Table}()
end
StructTypes.names(::Type{Script}) = ((:tables, :table), (:elements, :elem))
# //TODO Tackle some inconsistencies
function StructTypes.construct(::Type{Vector{T}}, xml::XMLDict.XMLDictElement) where T<:Script
    if haskey(xml, "script")
        xml = xml["script"]
    end
    if !(xml isa Vector) xml = [xml] end
    return Vector{T}([StructTypes.construct(Script, script) for script in xml])
end

#= end NSE scripts =#