text:
module InducedLSVMapDefinition
using ValidatedNumerics, TaylorSeries
using ..DynamicDefinition, ..Contractors, ..BasisDefinition

import ..C2BasisDefinition: C2Basis, dual_val, dual_der

export ApproxInducedLSV, preim, nbranches, plottable

"""
This class introduces an ApproximatedLSV
the induced map for the Liverani-Saussol-Vaienti maps
on the interval I = [0.5, 1].
The interval I is then mapped to [0,1]
"""

struct ApproxInducedLSV <: MarkovDynamic
	TnListPlottable::Array{Function, 1}
	nbranches::Integer
	domains::Array{Interval, 1}
	α::Real
	T::Function
end

function derleft(D::ApproxInducedLSV, x)
	@assert 0<=x<=0.5
	α = D.α
	return 1+2^α*(α+1)*x^α # check that α must be an interv)al?
end

function derderleft(D::ApproxInducedLSV, x)
	@assert 0<=x<=0.5
	α = D.α
	return 2^α*(α+1)*α*x^(α-1) # check that α must be an interval?
end

# coordinate change that maps [0.5, 1] to [0, 1]
@inline CoordinateChange(x) = 2*x-1
# inverse coordinate change
@inline InvCoordinateChange(x) = x/2+0.5

function ShootingLSV(n, y, α, rigstep = 10; T = Float64)
	x = [Interval{T}(0.5, 1); Interval{T}(0, 0.5)*ones(Interval{T}, n-1)]
	f(x) = 0<=x<=0.5 ? x*(1+(2*x)^α) : 2x-1
	fprime(x) = 0<=x<=0.5 ? 1+(α+1)*(2*x)^α : 2
	return ShootingMethod(f, fprime, n, x, y, rigstep)
end

function GetDomains(α , branches; T = Float64)
    domains=Interval{T}[]
    left = Interval{T}(0.5)
    for i in branches:-1:2
    	right = ShootingLSV(i, 0.5, α)[1]
		push!(domains,union(left,right))
    	left = right
    #for i=2:branches
    #    push!(domains, interval(ShootingLSV(i, 0.5, α)[1].lo, ShootingLSV(i-1, 0.5, α)[1].hi))
    end
    push!(domains, union(left, Interval{T}(0.75)))
    push!(domains, union(Interval{T}(0.75), Interval{T}(1)))

    return domains
end

function DefineTnPlottable(α, branches)
    Tleft(x)=x*(1+(2*x)^α)

    Tnlist = Array{Function, 1}(undef, branches+1)
    Tnlist[branches+1] = x->2*x-1

    for i=branches:-1:2
        Tnlist[i] = x-> Tleft(Tnlist[i+1](x))
    end

    right = mid(ShootingLSV(branches, 0.5, α)[1])

    Tnlist[1] = x-> 0.5*1/(right-0.5)*(x-0.5)+0.5
    return Tnlist
end

"""
This constructor builds the induced LSV map on [0.5, 1],
truncated with k branches
"""
function ApproxInducedLSV(α::Real, k::Int)
	nbranches = k+1
	domains = GetDomains(α, k)
	TnListPlottable = DefineTnPlottable(α, k)
	return ApproxInducedLSV(TnListPlottable, nbranches, domains, α, x-> 0<=x<=0.5 ? x*(1+(2*x)^α) : 2x-1)
end

domains(D::ApproxInducedLSV)=CoordinateChange.(D.domains)

DynamicDefinition.nbranches(S::ApproxInducedLSV)=S.nbranches

DynamicDefinition.is_full_branch(S::ApproxInducedLSV) = true

function _T(x, domains, TnList)
	x = InvCoordinateChange(x)
	for (i, I) in enumerate(domains)
		if x in I
			return CoordinateChange(TnList[i](x))
		end
	end
	return 0
end

function _T(x::Taylor1, domains, TnList)
	x = InvCoordinateChange(x)
	for (i, I) in enumerate(domains)
		if x[0] in I
			#@info x, I
			return CoordinateChange(TnList[i](x))
		end
	end
	return 0
end

#function _T(x::Taylor1{Interval}, domains, TnList)
#	x = InvCoordinateChange(x)
#	for (i, I) in enumerate(domains)
#		if x[0] ⊆ I
			#@info x, I
#			return CoordinateChange(TnList[i](x))
#		end
#	end
#	return 0
#end

DynamicDefinition.plottable(D::ApproxInducedLSV, x) = _T(x, D.domains, D.TnListPlottable)

function DynamicDefinition.preim(D::ApproxInducedLSV, k, y, ϵ)
	@assert 1 <= k <= D.nbranches
	_y = InvCoordinateChange(y)
	if k == 1
		right = ShootingLSV(D.nbranches-1, 0.5, D.α)[1]
		_x = (2*_y-1)*(right-0.5)+0.5
		return CoordinateChange(_x)
	elseif k == nbranches
		_x = (_y+1)/2
		return CoordinateChange(_x)
	else
		_x = ShootingLSV(D.nbranches-k+1, _y, D.α)[1]
		return CoordinateChange(_x)
	end
end

#returns the preimage of a point in a branch with the der and the second derivative
function preimwithder_derder(D::ApproxInducedLSV, k, y, ϵ) #::NTuple{Interval, 3}
	@assert 1 <= k <= D.nbranches

	y = Interval(y) #hack, please check
	_y = InvCoordinateChange(y)

	if k == 1 # the manufactured branch
		right = ShootingLSV(D.nbranches-1, 0.5, D.α)[1]
		_x = (2*_y-1)*(right-0.5)+0.5
		return CoordinateChange(_x), 0.5/(right-0.5), 0

	elseif k == D.nbranches # the linear branch
		_x = (_y+1)/2
		return CoordinateChange(_x), 2, 0
	else
		orbit_x = ShootingLSV(D.nbranches-k+1, _y, D.α)

		der = 2
		derder = 0
		for i in 2:length(orbit_x)
			dx = derleft(D, orbit_x[i])
			ddx = derderleft(D, orbit_x[i])
			derder = ddx*der^2+dx*derder
			der*=dx
		end

		return CoordinateChange(orbit_x[1]), der, derder/2
		# this /2 follows from the coordinate change
		# i.e., we are looking at the map F(x) = ψ(f(ϕ(x))) where
		# ϕ maps linearly [0,1]->[0.5, 1] and ψ is the
		# inverse map
		# by a direct computation we get  F''(x) = ψ'(f(ϕ(x)))f''(ϕ(x))(ϕ'(x))^2
	end
end

"""
Return (in an iterator) the pairs (i, (x, |T'(x)|)) where x is a preimage of p[i], which
describe the "dual" L* evaluation(p[i])
"""
function Base.iterate(S::DualComposedWithDynamic{T, ApproxInducedLSV}, state = (1, 1)) where T<:C2Basis
	i, k = state

	if i == length(S.basis)+1
			return nothing
	end

	n = length(S.basis.p)

	if i <= n
		x, der, derder = preimwithder_derder(S.dynamic, k, S.basis.p[i], S.ϵ)
		ret = x, (f, fprime) -> dual_val(f, fprime, x, der, derder)
	else
		x, der, derder = preimwithder_derder(S.dynamic, k, S.basis.p[i-n], S.ϵ)
		ret = x, (f, fprime) -> dual_der(f, fprime, x, der, derder)
	end
	if k == nbranches(S.dynamic)
		return ((i, ret), (i+1, 1))
	else
		return ((i, ret), (i, k+1))
	end
end

function iterate_LSV(x, i, α)
	@assert i>0
	x = 2*x-1
	for j in 2:i
		x = x*(1+(2*x)^α)
	end
	return x
end

end

using TaylorSeries: Taylor1

function bound_gamma_norm_gprime(α ,small_k, big_k, γ; T = Float64)
	#domains = reverse(InducedLSVMapDefinition.GetDomains(α, big_k+1; T= T))
	sum = 0.0
	for i in small_k:big_k
		right = InducedLSVMapDefinition.ShootingLSV(i-1, 0.5, α)[1]
		@info right
		#@info domains[i].lo, domains[i].hi
		f(x) = InducedLSVMapDefinition.iterate_LSV(x, i, α)
		fprime(x) = f(Taylor1([x, 1], 1))[1]
		fsecond(x) = 2*f(Taylor1([x, 1], 2))[2]
		tempered_distorsion(x)= abs((f(x)^γ)*1/fprime(x))
		dom = Interval(0.5, right.hi)
		tol = diam(dom)*2^(-10)
		dist = maximise(tempered_distorsion, dom, tol= tol)
		@info dist[1]
		sum+= dist[1]
	end
	return sum
end

function bound_inf_norm_gsecond_over_gprime(α ,small_k, big_k, γ; T = Float64)
	#domains = reverse(InducedLSVMapDefinition.GetDomains(α, big_k+1; T= T))
	glob_dist = 0.0
	right = InducedLSVMapDefinition.ShootingLSV(small_k-1, 0.5, α)[1]
	for i in small_k:big_k
		left = InducedLSVMapDefinition.ShootingLSV(i, 0.5, α)[1]
		@info i, left, right
		#@info domains[i].lo, domains[i].hi
		f(x) = InducedLSVMapDefinition.iterate_LSV(x, i, α)
		fprime(x) = f(Taylor1([x, 1], 1))[1]
		fsecond(x) = 2*f(Taylor1([x, 1], 2))[2]
		tempered_distorsion(x)= abs(fsecond(x)/(fprime(x))^2)
		dom = Interval(left.lo, right.hi)
		tol = diam(dom)*2^(-10)
		dist = maximise(tempered_distorsion, dom, tol= tol)
		glob_dist = max(dist[1].hi, glob_dist)
		@info dist[1]
		right = left
	end
	return glob_dist
end

max_stars_repo_name: orkolorko/ResponseManneville

max_issues_repo_name: orkolorko/ResponseManneville

max_forks_repo_name: orkolorko/ResponseManneville


text:
requiredTypes = [ "AbstractLevel" ]

for reqType in requiredTypes
    if !isdefined( Symbol( uppercase( string( reqType[ 1 ] ) ) * reqType[ 2:end ] ) )
        include(  reqType * ".jl" )
    end  # if !isdefined( Symbol( ...
end

export Rank

type Rank <: AbstractLevel

    # Name of the Rank ex: LT, CPN
    Name::String

    # Code of the rank ex: LT => 250, CPN =>220
    RankCode::Int

    # Maximum and minimum stay in this Rank
    MinStay::Int
    MaxStay::Int

    # A list of allowed Assignment
    NextLevel::Array{String}

    Rank() = new("", 0, 0, 0, Array{String}())
    Rank(nm::String,
         code::Int;
         min::Int = 0,
         max::Int = 0,
         next::Array{String} = Array{String}()) = new(nm, code, min, max, next)

end

ext: jl

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
using CImGui

include("../Renderer.jl")
using .Renderer

let counter = 0
global function counter_ui()
    flag = CImGui.ImGuiWindowFlags_NoTitleBar |
           CImGui.ImGuiWindowFlags_NoResize |
           CImGui.ImGuiWindowFlags_AlwaysAutoResize |
           CImGui.ImGuiWindowFlags_NoSavedSettings |
           CImGui.ImGuiWindowFlags_NoFocusOnAppearing |
           CImGui.ImGuiWindowFlags_NoNav
    CImGui.Begin("Counter", Ref(true), flag)
        spacing = CImGui.GetStyle().ItemInnerSpacing.x
        CImGui.PushButtonRepeat(true)
        CImGui.ArrowButton("##left", CImGui.ImGuiDir_Left) && (counter-=1;)
        CImGui.SameLine(0.0, spacing)
        CImGui.ArrowButton("##right", CImGui.ImGuiDir_Right) && (counter+=1;)
        CImGui.PopButtonRepeat()
        CImGui.SameLine()
        CImGui.Text("$counter")
    CImGui.End()
end
end

Renderer.render(counter_ui, width=230, height=50, title="App: Counter-Vanilla")


text:
# Julia wrapper for header: /usr/local/include/libavcodec/videotoolbox.h
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    av_videotoolbox_alloc_context,
    av_videotoolbox_default_init,
    av_videotoolbox_default_init2,
    av_videotoolbox_default_free

function av_videotoolbox_alloc_context()
    ccall((:av_videotoolbox_alloc_context,libavcodec),Ptr{AVVideotoolboxContext},())
end

function av_videotoolbox_default_init(avctx)
    ccall((:av_videotoolbox_default_init,libavcodec),Cint,(Ptr{AVCodecContext},),avctx)
end

function av_videotoolbox_default_init2(avctx,vtctx)
    ccall((:av_videotoolbox_default_init2,libavcodec),Cint,(Ptr{AVCodecContext},Ptr{AVVideotoolboxContext}),avctx,vtctx)
end

function av_videotoolbox_default_free(avctx)
    ccall((:av_videotoolbox_default_free,libavcodec),Void,(Ptr{AVCodecContext},),avctx)
end

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
"""
    PhysicalMedium{T<:AbstractFloat,Dim,FieldDim}

An abstract type used to represent the physical medium, the dimension of
the field, and the number of spatial dimensions.
"""
abstract type PhysicalMedium{T<:AbstractFloat,Dim,FieldDim} end

"Extract the dimension of the field of this physical property"
field_dimension(p::PhysicalMedium{T,Dim,FieldDim}) where {T,Dim,FieldDim} = FieldDim

"Extract the dimension of the field of this type of physical property"
field_dimension(p::Type{P}) where {Dim,FieldDim,T,P<:PhysicalMedium{T,Dim,FieldDim}} = FieldDim

"Extract the dimension of the space that this physical property lives in"
spatial_dimension(p::PhysicalMedium{T,Dim,FieldDim}) where {Dim,FieldDim,T} = Dim

"Extract the dimension of the space that this type of physical property lives in"
spatial_dimension(p::Type{P}) where {Dim,FieldDim,T,P<:PhysicalMedium{T,Dim,FieldDim}} = Dim

"""
A basis for regular functions, that is, smooth functions. A series expansion in this basis should converge to any regular function within a ball.
"""
function regular_basis_function(medium::P, ω::T) where {P<:PhysicalMedium,T}
    error("No regular basis function implemented for this physics type.")
end

"""
Basis of outgoing wave. A series expansion in this basis should converge to any scattered field outside of a ball which contains the scatterer.
"""
function outgoing_basis_function(medium::P, ω::T) where {P<:PhysicalMedium,T}
    error("No outgoing basis function implmented for this physics type.")
end

"""
the field inside an AbstractParticle a some given point x.
"""
internal_field

"""
A tuples of vectors of the field close to the boundary of the shape. The field is calculated from sim::FrequencySimulation, but the PhysicalMedium inside and outside of the shape are assumed to be given by inside_medium and outside_medium.
"""
boundary_data

# estimate_regular_basisorder(medium::P, ka) where P<:PhysicalMedium = estimate_regular_basisorder(P, ka)

"""
    estimate_regular_basis_order(medium::P, ω::Number, radius::Number; tol = 1e-6)
"""
function estimate_regular_basisorder(medium::P, ω::Number, radius::Number; tol = 1e-6) where P<:PhysicalMedium

    @error "This is not complete"

    k = ω / real(medium.c)
    vs = regular_basis_function(medium, ω)

    # A large initial guess
    L = Int(round(4 * abs(k*radius)))

    xs = radius .* rand(spatial_dimension(medium),10)

    l = nothing
    while isnothing(l)
        meanvs = [
            mean(norm(vs(l, xs[:,i])) for i in axes(xs,2))
        for l = 1:L]

        meanvs = mean(abs.(vs(L, xs[:,i])) for i in axes(xs,2))
        normvs = [
            norm(meanvs[basisorder_to_basislength(P,i-1):basisorder_to_basislength(P,i)])
        for i = 1:L]
        l = findfirst(normvs .< tol)
        L = L + Int(round(abs(k * radius / 2.0))) + 1
    end

    return l
end


text:
using Test, FixedEffects

method_s = [:lsmr, :lsmr_threads, :lsmr_cores]

p1 = repeat(1:5, inner = 2)
p2 = repeat(1:5, outer = 2)
x = [ 0.5548445405298847 , 0.9444014472663531 , 0.0510866660400604 , 0.9415750229576445 , 0.697755708534771  , 0.9664962514198971 , 0.12752269572311858, 0.4633531422366297 , 0.03341608526498096, 0.1647934493047556]
fes = [FixedEffect(p1), FixedEffect(p2)]
r_ols =  [-0.2015993617092453,  0.2015993617092464, -0.2015993617092463,  0.2015993617092462, -0.2015993617092465,  0.2015993617092467, -0.2015993617092465,  0.2015993617092470, -0.2015993617092468,  0.20159936170924628]

c_lsmr,_,_ = solve_coefficients!(copy(x), deepcopy(fes), method=:lsmr)

for method in method_s
	println(method)
	(c, iter, conv) = solve_coefficients!(copy(x), deepcopy(fes), method = method)
	@test c ≈ c_lsmr
	(r, iter, conv) = solve_residuals!(copy(x), deepcopy(fes), method = method)
	@test r ≈ r_ols
	(r, iter, conv) = solve_residuals!([x x x x x], deepcopy(fes), method = method)
	@test r ≈ [r_ols r_ols r_ols r_ols r_ols]
end

if isdefined(FixedEffects, :FixedEffectSolverLSMRGPU)
	push!(method_s, :lsmr_gpu)
end
for method in method_s
	if method != :qr
		println("$method Float32")
		(c, iter, conv) = solve_coefficients!(copy(x), deepcopy(fes), method=method, double_precision = false)
		@test c ≈ c_lsmr rtol = 1e-3
		(r, iter, conv) = solve_residuals!(copy(x),deepcopy(fes), method=method, double_precision = false)
		@test Float32.(r) ≈ Float32.(r_ols)
		(r, iter, conv) = solve_residuals!([x x x x x], deepcopy(fes), method=method, double_precision = false)
		@test r ≈ [r_ols r_ols r_ols r_ols r_ols] rtol = 1e-3
	end
end

fe = FixedEffect([1, 2])
@test_throws "FixedEffects must have the same length as y" ỹ = solve_residuals!(ones(100), [fe])


text:
export rule

@rule typeof(*)(:in, Marginalisation) (m_out::PointMass, m_A::PointMass) = PointMass(mean(m_in1) * mean(m_in2))

@rule typeof(*)(:in, Marginalisation) (m_out::F, m_A::PointMass, meta::AbstractCorrection) where { F <: NormalDistributionsFamily } = begin
    A = mean(m_A)
    W = correction!(meta, A' * precision(m_out) * A)
    return convert(promote_variate_type(F, NormalWeightedMeanPrecision), A' * weightedmean(m_out), W)
end


text:
using ArgParse
using CairoMakie
using NCDatasets
using Statistics
const NC = NCDatasets

function quick_plots(ncfile, out_path)
    NC.Dataset(ncfile) do ds
        group_names = keys(ds.group)
        for group_name in group_names
            pName1, pName2 = split(group_name, ".")
            group = ds.group[group_name]  # fetch group data
            xaxis_values = group[pName1][:]  # physical values along x-axis
            yaxis_values = group[pName2][:]  # physical values along y-axis
            nxvals = length(xaxis_values)
            nyvals = length(yaxis_values)

            case_names = group["case"][:]
            n_cases = length(case_names)
            h_res = 600
            fig = Figure(resolution = (h_res * (1 + n_cases), h_res))
            gl_list = [fig[1, i] = GridLayout() for i in 1:n_cases]
            text_size = 23.0
            for (case_ind, case_name) in enumerate(case_names)
                gl = gl_list[case_ind]

                ax1 = Axis(gl[1, 1], title = case_name, titlesize = text_size)
                ax2 = Axis(gl[2, 1])

                z_data = group["loss_data"][:, :, case_ind, :]
                z_mean = nanmean(z_data, dims = 3, drop = true)
                z_std = nanstd(z_data, dims = 3, drop = true)

                if nxvals == 1
                    scatter!(ax1, vec(z_mean), yaxis_values, markersize = 15)
                    scatter!(ax2, vec(z_std), yaxis_values, markersize = 15)
                    linkxaxes!(ax1, ax2)
                elseif nyvals == 1
                    scatter!(ax1, xaxis_values, vec(z_mean), markersize = 15)
                    scatter!(ax2, xaxis_values, vec(z_std), markersize = 15)
                    linkyaxes!(ax1, ax2)
                else
                    hm1 = heatmap!(ax1, xaxis_values, yaxis_values, z_mean)
                    hidexdecorations!(ax1, ticks = false)
                    Colorbar(gl[1:2, 2], hm1)
                    hm2 = heatmap!(ax2, xaxis_values, yaxis_values, z_std, colorrange = extrema(z_mean))
                end

            end  # end cases

            # Sum of all cases
            gl = fig[1, n_cases + 1]

            ax1 = Axis(gl[1, 1], title = "Cases average", titlesize = text_size)
            ax2 = Axis(gl[2, 1])

            z_data = group["loss_data"][:, :, :, :]
            z_case_mean = nanmean(z_data, dims = 3, drop = true)
            z_mean = nanmean(z_case_mean, dims = 3, drop = true)
            z_std = nanstd(z_case_mean, dims = 3, drop = true)

            if nxvals == 1
                scatter!(ax1, vec(z_mean), yaxis_values, markersize = 15)
                scatter!(ax2, vec(z_std), yaxis_values, markersize = 15)
                linkxaxes!(ax1, ax2)
                # axis labels
                xtitle = "Loss ($pName1 = $(only(xaxis_values)))"
                ytitle = pName2
            elseif nyvals == 1
                scatter!(ax1, xaxis_values, vec(z_mean), markersize = 15)
                scatter!(ax2, xaxis_values, vec(z_std), markersize = 15)
                linkyaxes!(ax1, ax2)
                # axis labels
                xtitle = pName1
                ytitle = "Loss ($pName2 = $(only(yaxis_values)))"
            else
                hm1 = heatmap!(ax1, xaxis_values, yaxis_values, z_mean)
                hidexdecorations!(ax1, ticks = false)
                Colorbar(gl[1:2, 2], hm1)
                hm2 = heatmap!(ax2, xaxis_values, yaxis_values, z_std, colorrange = extrema(z_mean))
                # axis labels
                xtitle = pName1
                ytitle = pName2
            end

            Label(fig[2, :], text = xtitle, textsize = text_size)  # x-axis
            Label(fig[:, 0], text = ytitle, textsize = text_size, rotation = pi / 2)  # y-axis
            save(joinpath(out_path, "$(group_name).png"), fig)
        end  # end groups
    end  # end NC.Dataset
end  # end quick_plots()

# https://discourse.julialang.org/t/nanmean-for-3d-array/62569/4
# -->
nansum((s, n), x) = isnan(x) ? (s, n) : (s + x, n + 1)
division((s, n)) = s / n
nanmean(a; dims = :, drop = false) = begin
    res = division.(reduce(nansum, a, init = (zero(eltype(a)), 0), dims = dims))
    drop ? dropdims(res, dims = dims) : res
end

nanvar(a; dims = :, drop = false) = nanmean((a .- nanmean(a, dims = dims)) .^ 2, dims = dims, drop = drop)
nanstd(a; dims = :, drop = false) = sqrt.(nanvar(a, dims = dims, drop = drop))
# <--

function parse_commandline_lm()

    s = ArgParseSettings(; description = "Run data path input")

    @add_arg_table s begin
        "--sim_dir"
        help = "Grid search simulations directory"
        arg_type = String
    end

    return ArgParse.parse_args(s)  # parse_args(ARGS, s)
end

if abspath(PROGRAM_FILE) == @__FILE__
    args = parse_commandline_lm()
    sim_dir = args["sim_dir"]
    out_path = joinpath(sim_dir, "figures")
    mkpath(out_path)
    ncfile = joinpath(sim_dir, "loss.nc")
    quick_plots(ncfile, out_path)
end


text:
# Silhouette

# this function returns r of size (k, n), such that
# r[i, j] is the sum of distances of all points from cluster i to point j
#
function sil_aggregate_dists(k::Int, a::AbstractVector{Int}, dists::AbstractMatrix{T}) where T<:Real
    n = length(a)
    S = typeof((one(T)+one(T))/2)
    r = zeros(S, k, n)
    @inbounds for j = 1:n
        for i = 1:j-1
            r[a[i],j] += dists[i,j]
        end
        for i = j+1:n
            r[a[i],j] += dists[i,j]
        end
    end
    return r
end

"""
    silhouettes(assignments::AbstractVector, [counts,] dists)
    silhouettes(clustering::ClusteringResult, dists)

Compute *silhouette* values for individual points w.r.t. given clustering.

Returns the ``n``-length vector of silhouette values for each individual point.

# Arguments
 - `assignments::AbstractVector{Int}`: the vector of point assignments
   (cluster indices)
 - `counts::AbstractVector{Int}`: the optional vector of cluster sizes (how many
   points assigned to each cluster; should match `assignments`)
 - `clustering::ClusteringResult`: the output of some clustering method
 - `dists::AbstractMatrix`: ``n×n`` matrix of pairwise distances between
   the points
"""
function silhouettes(assignments::AbstractVector{<:Integer},
                     counts::AbstractVector{<:Integer},
                     dists::AbstractMatrix{T}) where T<:Real

    n = length(assignments)
    k = length(counts)
    k >= 2 || throw(ArgumentError("silhouettes() not defined for the degenerated clustering with a single cluster."))
    for j = 1:n
        (1 <= assignments[j] <= k) || throw(ArgumentError("Bad assignments[$j]=$(assignments[j]): should be in 1:$k range."))
    end
    sum(counts) == n || throw(ArgumentError("Mismatch between assignments ($n) and counts (sum(counts)=$(sum(counts)))."))
    size(dists) == (n, n) || throw(DimensionMismatch("The size of a distance matrix ($(size(dists))) doesn't match the length of assignment vector ($n)."))

    # compute average distance from each cluster to each point --> r
    r = sil_aggregate_dists(k, assignments, dists)
    # from sum to average
    @inbounds for j = 1:n
        for i = 1:k
            c = counts[i]
            if i == assignments[j]
                c -= 1
            end
            if c == 0
                r[i,j] = 0
            else
                r[i,j] /= c
            end
        end
    end

    # compute a and b
    # a: average distance w.r.t. the assigned cluster
    # b: the minimum average distance w.r.t. other cluster
    a = similar(r, n)
    b = similar(r, n)

    for j = 1:n
        l = assignments[j]
        a[j] = r[l, j]

        v = typemax(eltype(b))
        for i = 1:k
            @inbounds rij = r[i,j]
            if (i != l) && (rij < v)
                v = rij
            end
        end
        b[j] = v
    end

    # compute silhouette score
    sil = a   # reuse the memory of a for sil
    for j = 1:n
        if counts[assignments[j]] == 1
            sil[j] = 0
        else
            #If both a[i] and b[i] are equal to 0 or Inf, silhouettes is defined as 0
            @inbounds sil[j] = a[j] < b[j] ? 1 - a[j]/b[j] :
                               a[j] > b[j] ? b[j]/a[j] - 1 :
                               zero(eltype(r))
        end
    end
    return sil
end

silhouettes(R::ClusteringResult, dists::AbstractMatrix) =
    silhouettes(assignments(R), counts(R), dists)

function silhouettes(assignments::AbstractVector{<:Integer}, dists::AbstractMatrix)
    counts = fill(0, maximum(assignments))
    for a in assignments
        counts[a] += 1
    end
    silhouettes(assignments, counts, dists)
end


text:
mutable struct FunctionCallingAffect{funcFunc, funcatType, funcatCacheType}
    func::funcFunc
    funcat::funcatType
    funcat_cache::funcatCacheType
    func_everystep::Bool
    func_start::Bool
    funciter::Int
end

function (affect!::FunctionCallingAffect)(integrator,force_func = false)
    # see OrdinaryDiffEq.jl -> integrator_utils.jl, function funcvalues!
    while !isempty(affect!.funcat) && integrator.tdir*top(affect!.funcat) <= integrator.tdir*integrator.t # Perform funcat
        affect!.funciter += 1
        curt = pop!(affect!.funcat) # current time
        if curt != integrator.t # If <t, interpolate
            if typeof(integrator) <: ODEIntegrator
                # Expand lazy dense for interpolation
                DiffEqBase.addsteps!(integrator)
            end
            if typeof(integrator.u) <: Union{Number,SArray}
                curu = integrator(curt)
            else
                curu = first(get_tmp_cache(integrator))
                integrator(curu,curt) # inplace since func_func allocates
            end
            affect!.func(curu, curt, integrator)
        else # ==t, just func
            affect!.func(integrator.u, integrator.t, integrator)
        end
    end
    if affect!.func_everystep || force_func
        affect!.funciter += 1
        affect!.func(integrator.u, integrator.t, integrator)
    end
    u_modified!(integrator, false)
end

function functioncalling_initialize(cb, u, t, integrator)
    if cb.affect!.funciter != 0
        if integrator.tdir > 0
            cb.affect!.funcat = BinaryMinHeap(cb.affect!.funcat_cache)
        else
            cb.affect!.funcat = BinaryMaxHeap(cb.affect!.funcat_cache)
        end
        cb.affect!.funciter = 0
    end
    cb.affect!.func_start && cb.affect!(integrator)
end

"""
    FunctionCallingCallback(func;
                    funcat=Vector{Float64}(),
                    func_everystep=isempty(funcat),
                    func_start = true,
                    tdir=1)

A `DiscreteCallback` applied after every step to call `func(u,t,integrator)`
If `func_everystep`, every step of the integrator gives a `func` call.
If `funcat` is specified, the function is called at the given times, using
interpolation if necessary.
If the time `tdir` direction is not positive, i.e. `tspan[1] > tspan[2]`,
`tdir = -1` has to be specified.
"""
function FunctionCallingCallback(func;
                        funcat=Vector{Float64}(),
                        func_everystep=isempty(funcat),
                        func_start = true,
                        tdir=1)
    # funcat conversions, see OrdinaryDiffEq.jl -> integrators/type.jl
    funcat_vec = collect(funcat)
    if tdir > 0
        funcat_internal = BinaryMinHeap(funcat_vec)
    else
        funcat_internal = BinaryMaxHeap(funcat_vec)
    end
    affect! = FunctionCallingAffect(func, funcat_internal,
                                    funcat_vec, func_everystep, func_start, 0)
    condtion = (u, t, integrator) -> true
    DiscreteCallback(condtion, affect!;
                     initialize = functioncalling_initialize,
                     save_positions=(false,false))
end

export FunctionCallingCallback


text:
@testset "airy" begin
    @test_throws AmosException airyai(200im)
    @test_throws AmosException airybi(200)

    for T in [Float16, Float32, Float64,Complex{Float16}, Complex{Float32},Complex{Float64}]
        @test airyai(T(1.8)) ≈ 0.0470362168668458052247
        @test airyaiprime(T(1.8)) ≈ -0.0685247801186109345638
        @test airybi(T(1.8)) ≈ 2.595869356743906290060
        @test airybiprime(T(1.8)) ≈ 2.98554005084659907283
    end
    for T in [Complex{Float16}, Complex{Float32}, Complex{Float64}]
        z = convert(T,1.8 + 1.0im)
        @test airyaix(z) ≈ airyai(z) * exp(2/3 * z * sqrt(z))
        @test airyaiprimex(z) ≈ airyaiprime(z) * exp(2/3 * z * sqrt(z))
        @test airybix(z) ≈ airybi(z) * exp(-abs(real(2/3 * z * sqrt(z))))
        @test airybiprimex(z) ≈ airybiprime(z) * exp(-abs(real(2/3 * z * sqrt(z))))
    end
    @test_throws MethodError airyai(complex(big(1.0)))

    for x = -3:3
        @test airyai(x) ≈ airyai(complex(x))
        @test airyaiprime(x) ≈ airyaiprime(complex(x))
        @test airybi(x) ≈ airybi(complex(x))
        @test airybiprime(x) ≈ airybiprime(complex(x))
        if x >= 0
            @test airyaix(x) ≈ airyaix(complex(x))
            @test airyaiprimex(x) ≈ airyaiprimex(complex(x))
        else
            @test_throws DomainError airyaix(x)
            @test_throws DomainError airyaiprimex(x)
        end
        @test airybix(x) ≈ airybix(complex(x))
        @test airybiprimex(x) ≈ airybiprimex(complex(x))
    end
end

@testset "bessel functions" begin
    bessel_funcs = [(bessely0, bessely1, bessely), (besselj0, besselj1, besselj)]
    @testset "$z, $o" for (z, o, f) in bessel_funcs
        @test z(Float32(2.0)) ≈ z(Float64(2.0))
        @test o(Float32(2.0)) ≈ o(Float64(2.0))
        @test z(Float16(2.0)) ≈ z(Float64(2.0))
        @test o(Float16(2.0)) ≈ o(Float64(2.0))
        @test z(2) ≈ z(2.0)
        @test o(2) ≈ o(2.0)
        @test z(2.0 + im) ≈ f(0, 2.0 + im)
        @test o(2.0 + im) ≈ f(1, 2.0 + im)
    end
    @testset "besselj error throwing" begin
        @test_throws MethodError besselj(1.2,big(1.0))
        @test_throws MethodError besselj(1,complex(big(1.0)))
        @test_throws MethodError besseljx(1,big(1.0))
        @test_throws MethodError besseljx(1,complex(big(1.0)))
    end
    @testset "besselh" begin
        true_h133 = 0.30906272225525164362 - 0.53854161610503161800im
        @test besselh(3,1,3) ≈ true_h133
        @test besselh(-3,1,3) ≈ -true_h133
        @test besselh(Float32(3),1,Float32(3)) ≈ true_h133
        @test besselh(Float16(3),1,Float16(3)) ≈ true_h133
        @test besselh(3,2,3) ≈ conj(true_h133)
        @test besselh(-3,2,3) ≈ -conj(true_h133)
        @test besselh(1,1,0) ≈ complex(0,-Inf)
        @test besselh(1,2,0) ≈ complex(0,Inf)
        @testset "Error throwing" begin
            @test_throws AmosException besselh(1,5,0)
            @test_throws MethodError besselh(1,big(1.0))
            @test_throws MethodError besselh(1,complex(big(1.0)))
            @test_throws MethodError besselhx(1,big(1.0))
            @test_throws MethodError besselhx(1,complex(big(1.0)))
        end
    end
    @testset "besseli" begin
        true_i33 = 0.95975362949600785698
        @test besseli(3,3) ≈ true_i33
        @test besseli(-3,3) ≈ true_i33
        @test besseli(3,-3) ≈ -true_i33
        @test besseli(-3,-3) ≈ -true_i33
        @test besseli(Float32(-3),Complex{Float32}(-3,0)) ≈ -true_i33
        @test besseli(Float16(-3),Complex{Float16}(-3,0)) ≈ -true_i33
        true_im3p1_3 = 0.84371226532586351965
        @test besseli(-3.1,3) ≈ true_im3p1_3
        for i in [-5 -3 -1 1 3 5]
            @test besseli(i,0) == 0.0
            @test besseli(i,Float32(0)) == 0
            @test besseli(i,Complex{Float32}(0)) == 0
            @test besseli(i,Float16(0)) == 0
            @test besseli(i,Complex{Float16}(0)) == 0
        end
        @testset "Error throwing" begin
            @test_throws AmosException besseli(1,1000)
            @test_throws DomainError besseli(0.4,-1.0)
            @test_throws MethodError besseli(1,big(1.0))
            @test_throws MethodError besseli(1,complex(big(1.0)))
            @test_throws MethodError besselix(1,big(1.0))
            @test_throws MethodError besselix(1,complex(big(1.0)))
        end
    end
    @testset "besselj" begin
        @test besselj(0,0) == 1
        for i in [-5 -3 -1 1 3 5]
            @test besselj(i,0) == 0
            @test besselj(i,Float32(0)) == 0
            @test besselj(i,Complex{Float32}(0)) == 0.0
        end

        j33 = besselj(3,3.)
        @test besselj(3,3) == j33
        @test besselj(-3,-3) == j33
        @test besselj(-3,3) == -j33
        @test besselj(3,-3) == -j33
        @test besselj(3,3f0) ≈ j33
        @test besselj(3,complex(3.)) ≈ j33
        @test besselj(3,complex(3f0)) ≈ j33
        @test besselj(3,complex(3)) ≈ j33

        j43 = besselj(4,3.)
        @test besselj(4,3) == j43
        @test besselj(-4,-3) == j43
        @test besselj(-4,3) == j43
        @test besselj(4,-3) == j43
        @test besselj(4,3f0) ≈ j43
        @test besselj(4,complex(3.)) ≈ j43
        @test besselj(4,complex(3f0)) ≈ j43
        @test besselj(4,complex(3)) ≈ j43

        @test j33 ≈ 0.30906272225525164362
        @test j43 ≈ 0.13203418392461221033
        @test besselj(0.1, complex(-0.4)) ≈ 0.820421842809028916 + 0.266571215948350899im
        @test besselj(3.2, 1.3+0.6im) ≈ 0.01135309305831220201 + 0.03927719044393515275im
        @test besselj(1, 3im) ≈ 3.953370217402609396im
        @test besselj(1.0,3im) ≈ besselj(1,3im)

        true_jm3p1_3 = -0.45024252862270713882
        @test besselj(-3.1,3) ≈ true_jm3p1_3
        @test besselj(Float16(-3.1),Float16(3)) ≈ true_jm3p1_3

        @testset "Error throwing" begin
            @test_throws DomainError    besselj(0.1, -0.4)
            @test_throws AmosException besselj(20,1000im)
            @test_throws MethodError besselj(big(1.0),3im)
        end
    end

    @testset "besselk" begin
        true_k33 = 0.12217037575718356792
        @test besselk(3,3) ≈ true_k33
        @test besselk(Float32(3),Float32(3)) ≈ true_k33
        @test besselk(Float16(3),Float16(3)) ≈ true_k33
        @test besselk(-3,3) ≈ true_k33
        true_k3m3 = -0.1221703757571835679 - 3.0151549516807985776im
        @test besselk(3,complex(-3)) ≈ true_k3m3
        @test besselk(-3,complex(-3)) ≈ true_k3m3
        # issue #6564
        @test besselk(1.0,0.0) == Inf
        @testset "Error throwing" begin
            @test_throws AmosException besselk(200,0.01)
            @test_throws DomainError besselk(3,-3)
            @test_throws MethodError besselk(1,big(1.0))
            @test_throws MethodError besselk(1,complex(big(1.0)))
            @test_throws MethodError besselkx(1,big(1.0))
            @test_throws MethodError besselkx(1,complex(big(1.0)))
        end
    end

    @testset "bessely" begin
        y33 = bessely(3,3.)
        @test bessely(3,3) == y33
        @test bessely(3.,3.) == y33
        @test bessely(3,Float32(3.)) ≈ y33
        @test bessely(-3,3) ≈ -y33
        @test y33 ≈ -0.53854161610503161800
        @test bessely(3,complex(-3)) ≈ 0.53854161610503161800 - 0.61812544451050328724im
        @testset "Error throwing" begin
            @test_throws AmosException bessely(200.5,0.1)
            @test_throws DomainError bessely(3,-3)
            @test_throws DomainError bessely(0.4,-1.0)
            @test_throws DomainError bessely(0.4,Float32(-1.0))
            @test_throws DomainError bessely(1,Float32(-1.0))
            @test_throws DomainError bessely(0.4,BigFloat(-1.0))
            @test_throws DomainError bessely(1,BigFloat(-1.0))
            @test_throws DomainError bessely(Cint(3),Float32(-3.))
            @test_throws DomainError bessely(Cint(3),Float64(-3.))

            @test_throws MethodError bessely(1.2,big(1.0))
            @test_throws MethodError bessely(1,complex(big(1.0)))
            @test_throws MethodError besselyx(1,big(1.0))
            @test_throws MethodError besselyx(1,complex(big(1.0)))
        end
    end

    @testset "sphericalbesselj" begin
        @test sphericalbesselj(1, 1)      ≈ 0.3011686789397568
        @test sphericalbesselj(10, 5.5)   ≈ 0.0009369210263385842
        @test sphericalbesselj(1.25, 5.5) ≈ -0.1123558799930763
        @test sphericalbesselj(1.25, -5.5+0im) ≈ 0.079447604649286 + 0.079447604649286im

        @test sphericalbesselj(0, 0.01) ≈ 0.999983333416666
        @test sphericalbesselj(0, 0)    == 1.0
        @test sphericalbesselj(1, 0)    == 0.0
        @test sphericalbesselj(1, 0.01) ≈ 0.003333300000119047

        @test_throws DomainError sphericalbesselj(1.25, -5.5)
    end

    @testset "sphericalbessely" begin
        @test sphericalbessely(1, 1)      ≈ -1.381773290676036
        @test sphericalbessely(10, 5.5)   ≈ -10.89087037026398
        @test sphericalbessely(1.25, 5.5) ≈ 0.148322390312228
        @test sphericalbessely(1.25, -5.5+0im) ≈ -0.054015441306998 - 0.104879767991574im

        @test sphericalbessely(0, 1e-5) ≈ -99999.9999950000000
        @test sphericalbessely(1, 1e-5) ≈ -1e10

        @test_throws DomainError sphericalbessely(1.25, -5.5)
        @test_throws AmosException sphericalbessely(1, 0)
    end

    @testset "besselhx" begin
        for elty in [Complex{Float16},Complex{Float32},Complex{Float64}]
            z = convert(elty, 1.0 + 1.9im)
            @test besselhx(1.0, 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
            @test besselhx(Float32(1.0), 1, z) ≈ convert(elty,-0.5949634147786144 - 0.18451272807835967im)
        end
        @testset "Error throwing" begin
            @test_throws MethodError besselh(1,1,big(1.0))
            @test_throws MethodError besselh(1,1,complex(big(1.0)))
            @test_throws MethodError besselhx(1,1,big(1.0))
            @test_throws MethodError besselhx(1,1,complex(big(1.0)))
        end
    end
    @testset "scaled bessel[ijky] and hankelh[12]" begin
        for x in (1.0, 0.0, -1.0), y in (1.0, 0.0, -1.0), nu in (1.0, 0.0, -1.0)
            z = Complex{Float64}(x + y * im)
            z == zero(z) || @test hankelh1x(nu, z) ≈ hankelh1(nu, z) * exp(-z * im)
            z == zero(z) || @test hankelh2x(nu, z) ≈ hankelh2(nu, z) * exp(z * im)
            (nu < 0 && z == zero(z)) || @test besselix(nu, z) ≈ besseli(nu, z) * exp(-abs(real(z)))
            (nu < 0 && z == zero(z)) || @test besseljx(nu, z) ≈ besselj(nu, z) * exp(-abs(imag(z)))
            z == zero(z) || @test besselkx(nu, z) ≈ besselk(nu, z) * exp(z)
            z == zero(z) || @test besselyx(nu, z) ≈ bessely(nu, z) * exp(-abs(imag(z)))
        end
        @test besselkx(1, 0) == Inf
        for i = [-5 -3 -1 1 3 5]
            @test besseljx(i,0) == 0
            @test besselix(i,0) == 0
            @test besseljx(i,Float32(0)) == 0
            @test besselix(i,Float32(0)) == 0
            @test besseljx(i,Complex{Float32}(0)) == 0
            @test besselix(i,Complex{Float32}(0)) == 0
            @test besseljx(i,Float16(0)) == 0
            @test besselix(i,Float16(0)) == 0
            @test besseljx(i,Complex{Float16}(0)) == 0
            @test besselix(i,Complex{Float16}(0)) == 0
        end
        @testset "Error throwing" begin
            @test_throws AmosException hankelh1x(1, 0)
            @test_throws AmosException hankelh2x(1, 0)
            @test_throws AmosException besselix(-1.01, 0)
            @test_throws AmosException besseljx(-1.01, 0)
            @test_throws AmosException besselyx(1, 0)
            @test_throws DomainError besselix(0.4,-1.0)
            @test_throws DomainError besseljx(0.4, -1.0)
            @test_throws DomainError besselkx(0.4,-1.0)
            @test_throws DomainError besselyx(0.4,-1.0)
        end
    end
    @testset "issue #6653" begin
        @testset "$f" for f in (besselj,bessely,besseli,besselk,hankelh1,hankelh2)
            @test f(0,1) ≈ f(0,Complex{Float64}(1))
            @test f(0,1) ≈ f(0,Complex{Float32}(1))
            @test f(0,1) ≈ f(0,Complex{Float16}(1))
        end
    end
end


text:
@testset "Test Reserves from Thermal Dispatch" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_uc = build_system("c_sys5_uc"; add_reserves = true)
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, model_template, c_sys5_uc; use_parameters = p)
        moi_tests(op_problem, p, 648, 0, 120, 216, 72, false)
    end
end

@testset "Test Reserves from Thermal Standard UC" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalBasicUnitCommitment),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :UpReserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_uc = build_system("c_sys5_uc"; add_reserves = true)
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, model_template, c_sys5_uc; use_parameters = p)
        moi_tests(op_problem, p, 1008, 0, 240, 216, 192, true)
    end
end

@testset "Test Upwards Reserves from Renewable Dispatch" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(RenewableDispatch, RenewableFullDispatch),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_re = build_system("c_sys5_re"; add_reserves = true)
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, model_template, c_sys5_re; use_parameters = p)
        moi_tests(op_problem, p, 360, 0, 72, 48, 72, false)
    end
end

@testset "Test Reserves from Storage" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
        :Storage => DeviceModel(GenericBattery, BookKeeping),
        # Added here to test it doesn't add reserve variables
        :Ren => DeviceModel(RenewableDispatch, FixedOutput),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_bat = build_system("c_sys5_bat"; add_reserves = true)
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, model_template, c_sys5_bat; use_parameters = p)
        moi_tests(op_problem, p, 408, 0, 192, 264, 96, false)
    end
end

@testset "Test Reserves from Hydro" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(HydroEnergyReservoir, HydroDispatchRunOfRiver),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_hyd = build_system("c_sys5_hyd"; add_reserves = true)
    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, model_template, c_sys5_hyd; use_parameters = p)
        moi_tests(op_problem, p, 240, 0, 24, 96, 72, false)
    end
end

@testset "Test Reserves from with slack variables" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_uc = build_system("c_sys5_uc"; add_reserves = true)
    for p in [true, false]
        op_problem = OperationsProblem(
            TestOpProblem,
            model_template,
            c_sys5_uc;
            use_parameters = p,
            services_slack_variables = true,
            balance_slack_variables = true,
        )
        moi_tests(op_problem, p, 504, 0, 120, 192, 24, false)
    end
end

@testset "Test AGC" begin
    c_sys5_reg = build_system("c_sys5_reg")
    # End of the system creation code.
    devices = Dict(
        :Loads => DeviceModel(PowerLoad, StaticPowerLoad),
        :Regulation_thermal =>
            DeviceModel(RegulationDevice{ThermalStandard}, DeviceLimitedRegulation),
    )
    services = Dict(:AGC => ServiceModel(AGC, PIDSmoothACE))

    @test_throws ArgumentError template_agc_reserve_deployment(devices = devices)

    template_agc = template_agc_reserve_deployment()
    agc_problem = OperationsProblem(AGCReserveDeployment, template_agc, c_sys5_reg)
    # These values might change as the AGC model is refined
    moi_tests(agc_problem, false, 720, 0, 480, 0, 384, false)
end

@testset "Test GroupReserve from Thermal Dispatch" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
        :GroupReserve => ServiceModel(StaticReserveGroup{ReserveDown}, GroupReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_uc = build_system("c_sys5_uc"; add_reserves = true)
    services = get_components(Service, c_sys5_uc)
    contributing_services = Vector{Service}()
    for service in services
        push!(contributing_services, service)
    end
    groupservice = StaticReserveGroup{ReserveDown}(;
        name = "init",
        available = true,
        requirement = 0.0,
        ext = Dict{String, Any}(),
    )
    add_service!(c_sys5_uc, groupservice, contributing_services)

    for p in [true, false]
        op_problem =
            OperationsProblem(TestOpProblem, model_template, c_sys5_uc; use_parameters = p)
        moi_tests(op_problem, p, 648, 0, 120, 240, 72, false)
    end
end

@testset "Test GroupReserve Errors" begin
    devices = Dict{Symbol, DeviceModel}(
        :Generators => DeviceModel(ThermalStandard, ThermalDispatch),
        :Loads => DeviceModel(PowerLoad, PSI.StaticPowerLoad),
    )
    branches = Dict{Symbol, DeviceModel}()
    services_template = Dict{Symbol, PSI.ServiceModel}(
        :Reserve => ServiceModel(VariableReserve{ReserveUp}, RangeReserve),
        :DownReserve => ServiceModel(VariableReserve{ReserveDown}, RangeReserve),
        :ORDC => ServiceModel(ReserveDemandCurve{ReserveUp}, StepwiseCostReserve),
        :GroupReserve => ServiceModel(StaticReserveGroup{ReserveDown}, GroupReserve),
    )
    model_template = OperationsProblemTemplate(
        CopperPlatePowerModel,
        devices,
        branches,
        services_template,
    )
    c_sys5_uc = build_system("c_sys5_uc"; add_reserves = true)
    services = get_components(Service, c_sys5_uc)
    contributing_services = Vector{Service}()
    for service in services
        push!(contributing_services, service)
    end
    groupservice = StaticReserveGroup{ReserveDown}(;
        name = "init",
        available = true,
        requirement = 0.0,
        ext = Dict{String, Any}(),
    )
    add_service!(c_sys5_uc, groupservice, contributing_services)

    off_service = VariableReserve{ReserveUp}("Reserveoff", true, 0.6, 10)
    push!(groupservice.contributing_services, off_service)

    @test_logs(
        (:error, r"is not stored"),
        match_mode = :any,
        @test_throws InfrastructureSystems.InvalidValue op_problem = OperationsProblem(
            TestOpProblem,
            model_template,
            c_sys5_uc;
            use_parameters = false,
        )
    )

end


text:
function findnearest(p, data)
    maxindex = size(data,2);
    distvector = zeros(maxindex,1);
    for i = 1:maxindex
        distvector[i] = norm(data[:,i]- p);
    end
#println(size(vec(distvector)))
    index = sortperm(vec(distvector));
    ni1 = index[1];
    ni2 = index[2];
    n1 = data[:,ni1];
    n2 = data[:,ni2];
   return n1, n2, ni1, ni2, distvector
end


text:
__precompile__()

module Postprocessing

using Gencode, Preprocessing

export fetchsolution, vis, exasim, producecode

include("createcgcells.jl");
include("createcggrid.jl");
include("getcelltype.jl");
include("getsolution.jl");
include("getmeansolution.jl");
include("fetchsolution.jl");
include("fetchresidual.jl");
include("vtuwrite.jl");
include("pvdwrite.jl");
include("vis.jl");
include("exasim.jl");
include("producecode.jl");
include("generatecode.jl");

end


text:
import Pkg
Pkg.activate(".")
using Revise
using SimpleSDMLayers # dev version
using Statistics

temperature = worldclim(1)
precipitation = worldclim(12)
wc_vars = worldclim(1:19)

# Extrema overload
minimum(temperature)
maximum(temperature)
extrema(temperature)

# Broadcast type
broadcast(abs, temperature)
maximum(abs(temperature))

# Max-min
@time min(temperature, precipitation)
@time max(temperature, precipitation)
@time test = reduce(max, wc_vars)
wc_vars_resp = convert.(SimpleSDMResponse, wc_vars)
@time test = reduce(max, wc_vars_resp)

minifloat = SimpleSDMPredictor([-1 0 1 nothing])
miniint = SimpleSDMResponse([-1.0 0.0 1.0 nothing])
mini32 = SimpleSDMPredictor([Float32.([-1.0 0.0 1.0])... nothing])
min(miniint, minifloat)
@time reduce(min, [miniint, minifloat, mini32])

# Mean-std
mean(wc_vars)
std(wc_vars)
# Attempt 1
mininothing = SimpleSDMPredictor([nothing 1.0 nothing nothing])
coordsnothing = map(x -> findall(isnothing, x.grid), [minifloat, miniint, mini32, mininothing])
uniquenothing = unique(coordsnothing)
unique(reduce(vcat, uniquenothing)) |> sort
# Need InvertedIndex...
# Attempt 2
coordsnotnothing = map(x -> findall(!isnothing, x.grid), [minifloat, miniint, mini32, mininothing])
uniquenot = unique(coordsnotnothing)
@time reduce(intersect, uniquenot)
@time intersect(uniquenot...)

# Math operations
temperature + precipitation
temperature - precipitation
temperature * precipitation
temperature + precipitation
reduce(+, wc_vars)


text:
include("plot_measurements.jl")
include("plot_observables.jl")
include("plot_ranking.jl")
include("plot_intervals.jl")


text:
using DataValues
using Dates

a, b = DataValue("14:00:00"), DataValue("15:15:00")

Dates.value(DateTime(a, "HH:MM:SS") - DateTime(b, "HH:MM:SS"))

c = DateTime(a, "HH:MM:SS") - DateTime(b, "HH:MM:SS")


text:
"""
    function get_PTDF_thermal_constprice(args..., kwargs...)

Identify generators which cost function is a constant price, locate their buses and return
the corresponding PTDF matrix. Skips Slack generators, Generators which Active Power is
binding and Generators with non zero slope. Note: For a quadratic cost function of the form:
αPg²+ βPg + γ, a constant price generator will have the parameter α == 0.

# Arguments
- `sys::System`:                            Power system in p.u. (from PowerSystems.jl)
- `res::PowerSimulations.
    OperationsProblemResults`:              Results of the solved OPF for the system
                                             (from PowerSimulations.jl)
- `PTDF_matrix::Any`:                       PTDF matrix as a direct output from
                                            PTDF_matrix = PTDF(sys), see PowerSystems.jl

# Keywords
- `dual_gen_tol::Float64 = 1e-1`:           Tolerance to identify any binding generators

# Throws
- `ERROR`:                                  PTDF_matrix has no field data. The argument must
                                            be a direct output from PTDF_matrix = PTDF(sys)
                                            see PowerSystems.jl
"""
function get_PTDF_thermal_constprice(
    sys::System,
    res::PowerSimulations.OperationsProblemResults,
    PTDF_matrix ::Array;
    dual_gen_tol::Float64 = 1e-1
    )

    # Define elements
    gen_loc = 0
    gen_constprice = 0 #number of generators with constant price (counter starter)
    gens_thermal_constantprice_busnumber = Array{Int64}(undef, 0)
    nl,nb = size(PTDF_matrix)
    PTDF_constantprice = Array{Float64}(undef,(nl, 0))

    # Get all thermal Generators
    gens_thermal = get_components(ThermalStandard, sys)

    # Get dual variables from the OPF results
    dualPgub = get_duals(res)[:P_ub__ThermalStandard__RangeConstraint]
    dualPglb = get_duals(res)[:P_lb__ThermalStandard__RangeConstraint]

    # Identify Generators and Locate their Buses
    for gen_thermal in gens_thermal
        gen_loc = gen_loc +1
        if get_bustype(gen_thermal.bus) ≠ BusTypes.REF  #skips slack
            bindub = !isapprox(dualPgub[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            bindlb = !isapprox(dualPglb[1,gen_loc], 0.0; atol = dual_gen_tol) #is binding?
            if !bindub && !bindlb #skips binding Pg gens Dual = 0
                #Get gen costs parameters Gencost = αPg²+ βPg + γ
                α, β = get_cost(get_variable(get_operation_cost(gen_thermal)))
                γ = get_fixed(get_operation_cost(gen_thermal))
                if α == 0 && (β ≠ 0 || γ ≠ 0)#skips non zero slope
                    gen_constprice = gen_constprice +1
                    resize!(gens_thermal_constantprice_busnumber,gen_constprice)
                    gens_thermal_constantprice_busnumber[gen_constprice] =
                     gen_thermal.bus.number
                end
            end
        end
    end

    # Parse the PTDF matrix
    if gens_thermal_constantprice_busnumber == Any[]
        PTDF_constantprice = Array{Float64}(undef,(nl,0))
    else
        sort!(gens_thermal_constantprice_busnumber)
        PTDF_constantprice = PTDF_matrix[:,gens_thermal_constantprice_busnumber]
    end

    return gens_thermal_constantprice_busnumber, PTDF_constantprice
end

max_stars_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_issues_repo_name: invenia/MaximumGeneratorProfitTCRDD

max_forks_repo_name: invenia/MaximumGeneratorProfitTCRDD


text:
# © Przemyslaw Szufel, 2018, run under Julia 1.0.0

#Usage
#julia -p 5 schelling2.jl 2x2 35X35 slurm=false plot_field=false plot_ascii=false max_step=100

using Random, Printf, Statistics
import Future

const cpu_grid = [2,2]
plot_field = false
slurm = false
plot_ascii = false
const field_size = [35,35]
max_step = 100
for arg in ARGS
    if startswith(arg,"max_step=")
        max_step = parse(Int,arg[length("max_step=")+1:end])
    elseif occursin("x", arg)
        ix = findfirst("x",arg)[1]
        cpu_grid[1]=parse(Int,arg[1:(ix-1)])
        cpu_grid[2]=parse(Int,arg[(ix+1):end])
    elseif occursin("X",arg)
        ix = findfirst("X",arg)[1]
        field_size[1]=parse(Int,arg[1:(ix-1)])
        field_size[2]=parse(Int,arg[(ix+1):end])
    elseif startswith(arg,"slurm=")
        slurm = parse(Bool,arg[length("slurm=")+1:end])
    elseif startswith(arg,"plot_field=")
        plot_field = parse(Bool,arg[length("plot_field=")+1:end])
    elseif startswith(arg,"plot_ascii=")
        plot_ascii = parse(Bool,arg[length("plot_ascii=")+1:end])
    else
        println("Sample usage\njulia  -p 5 schelling2.jl 2x2 35X35 slurm=false plot_field=false plot_ascii=false max_step=100")
        exit(1)
    end
end

const GRID_SIZE = (cpu_grid[1],cpu_grid[2])
plot_field=true
plot_ascii=true
#println(plot_field, plot_ascii, GRID_SIZE)

if slurm
    using ClusterManagers
    addprocs_slurm(GRID_SIZE[1]*GRID_SIZE[2]+1,job_name="schelling2", account="GC71-37", time="01:00:00", exename="/lustre/tetyda/home/pszufe/julia0.6.2_v002/usr/bin/julia")
    cd("/lustre/tetyda/home/pszufe")
else
    #aws version requires machinefile to be present, e.g.
    #julia --machinefile machinefile schelling2.jl .....
    #or use -p parameter
    #julia -p 17 schelling2.jl ............
    min_procs = GRID_SIZE[1]*GRID_SIZE[2]+1
    if nprocs() <= min_procs
        println("For this config at least ",min_procs, " workers are required")
        exit(1)
    end
end

if plot_field
    include("schelling2_plot.jl");
end

const p2p = true

const mywks = workers()[1:(GRID_SIZE[1]*GRID_SIZE[2]+1)]

@everywhere using Random, Printf, Statistics
@everywhere include("schelling2_everywhere.jl");
const modelC = ModelConfig(GRID_SIZE, field_size[1],field_size[2],[0.1,0.45,0.45],5,mywks[1:length(mywks)-1],mywks[length(mywks)])

@everywhere using ParallelDataTransfer
sendto(mywks, modelC = modelC)
sendto(modelC.master_id, p2p = p2p)

const n_rngs = length(modelC.slaves)+1
const rngs = accumulate(Future.randjump, [0; [big(10)^20 for i in 1:n_rngs-1]], init = MersenneTwister(0))

sendto(modelC.master_id, rng=rngs[length(rngs)])

const futures = Dict{Int,Distributed.Future}()

const worker_locations =  Dict{Int,Tuple{Int,Int}}()
const worker_grid = zeros(Int,GRID_SIZE)

function set_grid()
    ix=0
    for w in modelC.slaves
        i =  Int(floor(ix / GRID_SIZE[2])) +1
        j =  ix - ( (i-1)*GRID_SIZE[2] ) +1
        ix += 1
        #print("Creating worker ",w," ix=",ix," i=",i," j=",j,"\n")
        #creates a global variable at each worker
        sendto(w, node = NodeData(i, j, zeros(Int8,(0,0)),Vector{Tuple{Int,Int}}(),rngs[ix]))
        futures[w] = @spawnat w init_node()
        worker_locations[w] = (i,j)
        worker_grid[i,j] = w
    end
end
set_grid()

print("Waiting for ",length(modelC.slaves)," slaves\n")
flush(stdout)
sendto(modelC.master_id, worker_grid=worker_grid)
sendto(modelC.master_id, worker_locations=worker_locations)
for w in modelC.slaves
    init_ok = fetch(futures[w])
    sendto(w, worker_grid = worker_grid)
end
println("Completed creating workers max slave id ",maximum(modelC.slaves)," p2p=",p2p)
flush(stdout)

include("schelling2_master.jl")
fetch(@spawnat modelC.master_id include("schelling2_master.jl"))

function get_aggField(modelC::ModelConfig, GRID_SIZE::Tuple{Int,Int}, worker_grid::Matrix{Int},node_colors::Bool)::Matrix{Int8}
    dim_i=modelC.dim_i
    dim_j=modelC.dim_j
    aggField = zeros(Int8,modelC.dim_i*GRID_SIZE[1],modelC.dim_j*GRID_SIZE[2])
    for gi=1:GRID_SIZE[1], gj=1:GRID_SIZE[2]
        begin
            field = remotecall_fetch(get_field_interior,worker_grid[gi,gj])
            for i=1:modelC.dim_i, j=1:modelC.dim_j
                ii=dim_i*(gi-1)+i
                jj=dim_j*(gj-1)+j
                aggField[ii,jj]=field[i,j]
                if node_colors && aggField[ii,jj]==0
                    aggField[ii,jj]=length(modelC.probs_w)+GRID_SIZE[2]^(gi-1)+(gj-1)
                end
            end
        end
    end
    return aggField
end

function main_loop(modelC :: ModelConfig,worker_grid::Matrix{Int},worker_locations::Dict{Int,Tuple{Int,Int}},plotField::Bool,plotAscii::Bool,maxStep::Int)
    img = nothing
    for stepNo in 1:maxStep
        res = remotecall_fetch(do_item, modelC.master_id )
        println("[",stepNo,"] ",res)
        flush(stdout)

        if plotField
            aggField = get_aggField(modelC,modelC.GRID_SIZE,worker_grid,true)
            img = plot_schelling(aggField,stepNo,img)
        end
    end

    if plotAscii
        aggField = get_aggField(modelC,(min(GRID_SIZE[1],2),min(GRID_SIZE[2],2)),worker_grid,false)
        for i in 1:size(aggField)[1]
            for j=1:size(aggField)[2]
                if aggField[i,j]==0
                    print(".")
                else
                    print(aggField[i,j])
                end
            end
            print("\n")
        end
    end

end

main_loop(modelC,worker_grid,worker_locations,plot_field,plot_ascii,max_step)
if plot_field
    println("Press <ENTER>")
    readline()
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
using LazyArrays, FillArrays, ArrayLayouts, StaticArrays, Base64, Test
import LazyArrays: PaddedLayout, LayoutVector, MemoryLayout, paddeddata, ApplyLayout, sub_materialize, CachedVector
import Base: setindex

# padded block arrays have padded data that is also padded. This is to test this
struct PaddedPadded <: LayoutVector{Int} end

MemoryLayout(::Type{PaddedPadded}) = PaddedLayout{UnknownLayout}()
Base.size(::PaddedPadded) = (10,)
Base.getindex(::PaddedPadded, k::Int) = k ≤ 5 ? 1 : 0
paddeddata(a::PaddedPadded) = a

@testset "Padded" begin
    @testset "PaddedLayout" begin
        A = Vcat([1,2,3], Zeros(7))
        B = Vcat([1,2], Zeros(8))

        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.args[1] isa Vector
        @test C.args[2] isa Zeros{Float64}
        @test C == Vector(A) + Vector(B)

        @test colsupport(A, 1) == 1:3

        B = Vcat([1,2], Ones(8))

        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.args[1] isa Vector{Float64}
        @test C.args[2] isa Ones{Float64}
        @test C == Vector(A) + Vector(B)

        B = Vcat([1,2], randn(8))

        C = @inferred(A+B)
        @test C isa Vcat
        @test C == Vector(A) + Vector(B)

        B = Vcat(SVector(1,2), Ones(8))
        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.args[1] isa Vector{Float64}
        @test C.args[2] isa Ones{Float64}
        @test C == Vector(A) + Vector(B)

        A = Vcat(SVector(3,4), Zeros(8))
        B = Vcat(SVector(1,2), Ones(8))
        C = @inferred(A+B)
        @test C isa Vcat{Float64,1}
        @test C.args[1] isa SVector{2,Int}
        @test C.args[2] isa Ones{Float64}
        @test C == Vector(A) + Vector(B)

        @testset "multiple scalar" begin
            # We only do 1 or 2 for now, this should be redesigned later
            A = Vcat(1, Zeros(8))
            @test MemoryLayout(A) isa PaddedLayout{ScalarLayout}
            @test paddeddata(A) == 1
            B = Vcat(1, 2, Zeros(8))
            @test paddeddata(B) == [1,2]
            @test MemoryLayout(B) isa PaddedLayout{ApplyLayout{typeof(vcat)}}
            C = Vcat(1, cache(Zeros(8)));
            @test paddeddata(C) == [1]
            @test MemoryLayout(C) isa PaddedLayout{ApplyLayout{typeof(vcat)}}
            D = Vcat(1, 2, cache(Zeros(8)));
            @test paddeddata(D) == [1,2]
            @test MemoryLayout(D) isa PaddedLayout{ApplyLayout{typeof(vcat)}}
        end

        @testset "PaddedPadded" begin
            @test colsupport(PaddedPadded()) ≡ Base.OneTo(10)
            @test stringmime("text/plain", PaddedPadded()) == "10-element PaddedPadded:\n 1\n 1\n 1\n 1\n 1\n 0\n 0\n 0\n 0\n 0"
            @test dot(PaddedPadded(), PaddedPadded()) == 5
            @test dot(PaddedPadded(), 1:10) == dot(1:10, PaddedPadded()) == 15

            @test PaddedPadded()[1:7] isa Vcat

            @test norm(PaddedPadded()) ≈ sqrt(5)
        end

        @testset "Matrix padded" begin
            a = Vcat(randn(2,3), Zeros(3,3))
            b = Vcat(randn(3,3), Zeros(2,3))
            @test a + b isa Vcat
            @test a + b == Matrix(a) + Matrix(b)
        end
    end

    @testset "copyto!" begin
        a = Vcat(1, Zeros(10));
        c = cache(Zeros(11));
        @test MemoryLayout(typeof(a)) isa PaddedLayout
        @test MemoryLayout(typeof(c)) isa PaddedLayout{DenseColumnMajor}
        @test copyto!(c, a) ≡ c;
        @test c.datasize[1] == 1
        @test c == a

        a = Vcat(1:3, Zeros(10))
        c = cache(Zeros(13));
        @test MemoryLayout(typeof(a)) isa PaddedLayout
        @test MemoryLayout(typeof(c)) isa PaddedLayout{DenseColumnMajor}
        @test copyto!(c, a) ≡ c;
        @test c.datasize[1] == 3
        @test c == a

        @test dot(a,a) ≡ dot(a,c) ≡ dot(c,a) ≡ dot(c,c) ≡ 14.0

        a = Vcat(1:3, Zeros(5))
        c = cache(Zeros(13));
        @test copyto!(c, a) ≡ c;
        @test c.datasize[1] == 3
        @test c[1:8] == a

        a = cache(Zeros(13)); b = cache(Zeros(15));
        @test a ≠ b
        b = cache(Zeros(13));
        a[3] = 2; b[3] = 2; b[5]=0;
        @test a == b
    end

    @testset "vcat and padded" begin
        x,y = Vcat([1,2,3],Zeros(5)), Vcat(5, 1:7)
        @test MemoryLayout(x) isa PaddedLayout
        @test MemoryLayout(y) isa ApplyLayout{typeof(vcat)}
        @test x .+ y == y .+ x == Vector(x) .+ Vector(y)
        @test x .+ y isa Vcat
        @test y .+ x isa Vcat

        c = cache(Zeros(8));
        @test c + x == x + c
        @test c + y == y + c
    end

    @testset "vcat and Zeros" begin
        x,y = Vcat([1,2,3],Zeros(5)), Vcat(5, 1:7)
        @test x .+ Zeros(8) == Zeros(8) .+ x == x
        @test y .+ Zeros(8) == Zeros(8) .+ y == y
        @test x .* Zeros(8) ≡ Zeros(8) .* x ≡ Zeros(8)
        @test y .* Zeros(8) ≡ Zeros(8) .* y ≡ Zeros(8)
        @test x .\ Zeros(8) ≡ Zeros(8) ./ x ≡ Zeros(8)
        @test y .\ Zeros(8) ≡ Zeros(8) ./ y ≡ Zeros(8)
    end

    @testset "vcat and BroadcastArray" begin
        x,y,z = Vcat([1,2,3],Zeros(5)), Vcat(5, 1:7), BroadcastArray(exp,1:8)
        @test x .+ z == z .+ x == Array(x) .+ Array(z)
        @test y .+ z == z .+ y == Array(y) .+ Array(z)

        @test x .+ z isa CachedVector
        @test y .+ z isa BroadcastVector
    end
    @testset "subpadded" begin
        A = Hcat(1:10, Zeros(10,10))
        V = view(A,3:5,:)
        @test MemoryLayout(V) isa PaddedLayout
        @test A[parentindices(V)...] == copy(V) == Array(A)[parentindices(V)...]
        V = view(A,3:5,1:4)
        @test MemoryLayout(V) isa PaddedLayout
        @test @inferred(paddeddata(V)) == reshape(3:5,3,1)

        v = view(A,2,1:5)
        @test MemoryLayout(v) isa PaddedLayout
        @test paddeddata(v) == [2]
        @test A[2,1:5] == copy(v) == sub_materialize(v)

        @testset "Padded subarrays" begin
            a = Vcat([1,2,3],[4,5,6])
            @test sub_materialize(view(a,2:6)) == a[2:6]
            a = Vcat([1,2,3], Zeros(10))
            c = cache(Zeros(10)); c[1:3] = 1:3;
            v = view(a,2:4)
            w = view(c,2:4);
            @test MemoryLayout(typeof(a)) isa PaddedLayout{DenseColumnMajor}
            @test MemoryLayout(v) isa PaddedLayout{DenseColumnMajor}
            @test sub_materialize(v) == a[2:4] == sub_materialize(w)
            @test sub_materialize(v) isa Vcat
            @test sub_materialize(w) isa Vcat
            A = Vcat(Eye(2), Zeros(10,2))
            V = view(A, 1:5, 1:2)
            @test sub_materialize(V) == A[1:5,1:2]
        end
    end
    @testset "hcat" begin
        z = cache(Zeros(1,3));
        z[1] = 3;
        H = Hcat([1 2], z);
        @test MemoryLayout(H) isa PaddedLayout
        @test paddeddata(H) == [1 2 3]
        H = Hcat(1, 3, Zeros(1,3))
        @test MemoryLayout(H) isa PaddedLayout
        @test paddeddata(H) == [1 3]
    end
    @testset "padded broadcast" begin
        @testset "vector" begin
            a = Vcat([1], Zeros(3))
            b = Vcat([1,2],Zeros(2))
            v = Vcat(1, 1:3)
            n = Vcat(1, Zeros{Int}(3))
            c = cache(Zeros(4));
            @test @inferred(a + b) isa Vcat
            @test @inferred(a .* b) isa Vcat
            @test paddeddata(a .* b) == [1]
            @test a + v == v + a
            @test a .* v == v .* a

            @test n .+ n ≡ Vcat(2,Zeros{Int}(3))
            @test n .+ v ≡ n .+ v
            @test n .+ v ≡ v .+ n ≡ Vcat(2,1:3)
            @test n .+ b == b .+ n
            @test n .+ a == a .+ n
            @test n .+ c == c .+ n
            @test c.datasize == (0,)

            @test @inferred(n .* n) ≡ @inferred(n .* n) ≡ n
            @test @inferred(n .* a) ≡ @inferred(a .* n) ≡ Vcat(1,Zeros{Float64}(3))
            @test @inferred(n .* v) ≡ @inferred(v .* n) ≡ Vcat(1,Zeros{Int}(3))
            @test n .* c == c .* n

            @test view(a, 1:4) .* v == v .* view(a, 1:4)
        end
        @testset "matrix" begin
            a = Vcat([1 2], Zeros(3,2))
            b = Vcat([1 2; 3 4], Zeros(2,2))
            @test MemoryLayout(a + a) isa PaddedLayout
            @test a + a isa Vcat
            @test a + a == 2a
            @test a + b == b + a

            c = PaddedArray(randn(2,2), 4, 2)
            d = PaddedArray(randn(1,2), 4, 2)
            e = PaddedArray(randn(3,1), 4, 2)
            @test c + c == 2c
            @test c + a == a + c
            @test c + b == b + c
            @test c + d == d + c
            @test d + a == a + d
            @test d + b == b + d
            @test e + a == a + e
            @test e + b == b + e
            @test e + c == c + e
            @test e + d == d + e
        end
    end

    @testset "hvcat" begin
        P = ApplyArray(hvcat, 2, randn(4,5), Zeros(4,6), Zeros(6,5), Zeros(6,6))
        @test eltype(P) == Float64
        @test MemoryLayout(P) isa PaddedLayout
        @test P == hvcat(P.args...)
        @test @inferred(colsupport(P,3)) == Base.OneTo(4)
        @test @inferred(colsupport(P,7)) == Base.OneTo(0)
        @test @inferred(rowsupport(P,3)) == Base.OneTo(5)
        @test @inferred(rowsupport(P,7)) == Base.OneTo(0)
        @test copyto!(similar(P), P) == P

        @test P[3,:] isa Vcat
        @test P[3,1:11] == P[3,:]
        @test P[:,3] isa Vcat
        @test P[1:10,3] == P[:,3]
        @test P[6,:] isa Vcat
        @test P[1:10,6] == P[:,6]
        @test P[:,6] isa Vcat
        @test P[6,1:11] == P[6,:]
    end
    @testset "setindex" begin
        a = ApplyArray(setindex, 1:6, 5, 2)
        @test a == [1; 5; 3:6]
        @test_throws BoundsError a[7]
        a = ApplyArray(setindex, 1:6, [9,8,7], 1:3)
        @test a == [9; 8; 7; 4:6]
        @test_throws BoundsError a[7]

        a = ApplyArray(setindex, Zeros(5,5), 2, 2, 3)
        @test a[2,3] === 2.0
        @test a == setindex!(zeros(5,5),2,2,3)

        a = ApplyArray(setindex, Zeros(5,5), [4,5], 2:3, 3)
        @test a == setindex!(zeros(5,5),[4,5], 2:3, 3)

        a = ApplyArray(setindex, Zeros(5,5), [1 2 3; 4 5 6], 2:3, 3:5)
        @test a == setindex!(zeros(5,5),[1 2 3; 4 5 6], 2:3, 3:5)

        a = ApplyArray(setindex, Zeros(5), [1,2], Base.OneTo(2))
        @test MemoryLayout(a) isa PaddedLayout{DenseColumnMajor}
        @test paddeddata(a) == 1:2

        a = ApplyArray(setindex, Zeros(5,5), [1 2 3; 4 5 6], Base.OneTo(2), Base.OneTo(3))
        @test a == setindex!(zeros(5,5),[1 2 3; 4 5 6], 1:2, 1:3)
        @test MemoryLayout(a) isa PaddedLayout{DenseColumnMajor}
        @test paddeddata(a) == [1 2 3; 4 5 6]

        # need to add bounds checking
        @test_broken ApplyArray(setindex, Zeros(5,5), [1 2; 4 5], 2:3, 3:5)

        @test PaddedArray(1, 3) == PaddedVector(1,3) == [1; zeros(2)]
        @test PaddedArray(1, 3, 3) == PaddedMatrix(1, 3, 3) == [1 zeros(1,2); zeros(2,3)]
    end

    @testset "adjtrans" begin
        a = Vcat(1, Zeros(3))
        @test MemoryLayout(a') isa DualLayout{<:PaddedLayout}
        @test MemoryLayout(transpose(a)) isa DualLayout{<:PaddedLayout}
        @test paddeddata(a') ≡ 1
        b = Vcat(SVector(1,2), Zeros(3))
        @test paddeddata(b') ≡ SVector(1,2)'
        @test paddeddata(transpose(b)) ≡ transpose(SVector(1,2))
    end

    @testset "norm" begin
        a = Vcat(1, Zeros(3))
        c = cache(Zeros(4)); c[1] = 1
        @test norm(a) ≡ LinearAlgebra.normInf(c) ≡ LinearAlgebra.norm2(c) ≡ LinearAlgebra.norm1(c) ≡ LinearAlgebra.normp(c,2) ≡ 1.0
    end
end


text:
using StaticArrays
using FastClosures
"""
gellmann_matrices(;type=Float64)
Return the Gell-Mann matrices.
"""
function gellmann_matrices(;type=Float64, struct_const=false)
	λ1 = SMatrix{3, 3, Complex{type}}([0 1 0; 1 0 0; 0 0 0])
	λ2 = SMatrix{3, 3, Complex{type}}([0 -1im 0; 1im 0 0; 0 0 0])
	λ3 = SMatrix{3, 3, Complex{type}}([1 0 0; 0 -1 0; 0 0 0])
	λ4 = SMatrix{3, 3, Complex{type}}([0 0 1; 0 0 0; 1 0 0])
	λ5 = SMatrix{3, 3, Complex{type}}([0 0 -1im; 0 0 0;1im 0 0])
	λ6 = SMatrix{3, 3, Complex{type}}([0 0 0; 0 0 1; 0 1 0])
	λ7 = SMatrix{3, 3, Complex{type}}([0 0 0; 0 0 -1im; 0 1im 0])
	λ8 = SMatrix{3, 3, Complex{type}}(
		 [1/sqrt(type(3)) 0 0; 0 1/sqrt(type(3)) 0; 0 0 -2/sqrt(type(3))])
	λ9 = SMatrix{3, 3, Complex{type}}([1 0 0; 0 1 0; 0 0 1])
	λ = (λ1, λ2, λ3, λ4, λ5, λ6, λ7, λ8, λ9)

	if struct_const == true
		index = ((1, 2, 3), (1, 4, 7), (1, 6, 5), (2, 4, 6), (2, 5, 7), (3, 4, 5),
				 (3, 7, 6), (4, 5, 8), (6, 7, 8))

		f = zeros(type, 9)
		f[1] = type(1)
		f[2: 7] .= type(1)/2
		f[8: 9] .= sqrt(type(3))/2
		f = Tuple(f)

		return λ, index, f
	else
		return λ
	end
end

function vector2matrix(v, λ)
	M = v[1]*λ[1]
	for i = 2: 9
		@inbounds M = M + v[i]*λ[i]
	end
	return M
end
"""
	fidelity(U, U_ideal, dims)
Fidelity of the evolution operator U with respect to the ideal evolution operator U_ideal.
The Hilbert space is truncated in the `ndims` first dimensions.
"""
function fidelity(U, U_ideal, dims)
	U_ideal_dag = adjoint(U_ideal)
	M = U_ideal_dag[1: dims, 1: dims]*U[1: dims, 1: dims]
	M_dag = adjoint(M)
	fidelity = (tr(M*M_dag) + abs(tr(M))^2)/(dims*(dims + 1))
	return real(fidelity)
end
"""
	gell_mann_decomp(M)
Decomposition of a generic 3x3 Hermitian matrix M in the Gell-Mann basis.
"""
function gell_mann_decomp(M)
	c = zeros(9)
	c[9] = real(M[1,1] + M[2,2] + M[3,3])/3
	c[8] = (-real(M[3,3]) + c[9])*sqrt(3.0)/2
	c[1], c[2] = real(M[2,1]), imag(M[2,1])
	c[3] = real(M[1,1]) - c[9] - c[8]/sqrt(3.0)
	c[4], c[5] = real(M[3,1]), imag(M[3,1])
	c[6], c[7] = real(M[3,2]), imag(M[3,2])
	return c
end
"""
	H0(t, α, η, f; cr_terms::Bool = false)
Hamiltonian of an isolated 3 level system in the rotating frame.
# Arguments
- `t::Real`: time.
- `α::Real`: anarmonicity.
- `η::Real`: parameter that characterizes the strength of the 1<->2 transition relative
to the 0<->1 transition.
- `f::Function`: pulse envelope as a sunction of t.
- `cr_terms::Bool = false`: If true, counter rotating terms are taken in account.
"""
function H0(t, α, η, f; cr_terms::Bool=false)
	type = typeof(t)
	ft = f(t)/2
	if cr_terms == false
		H0t = @SVector [ft,
						0,
						0,
						ft*η,
						0,
						0,
						0,
						-α/sqrt(type(3)),
						α/3]
	else
		cos2t = cos(2*t)
		sin2t = sin(2*t)
		H0t = @SVector [ft*(1 + cos2t),
						-ft*sin2t,
						0,
						ft*(1 + cos2t)*η,
						ft*sin2t*η,
						0,
						0,
						-α/sqrt(type(3)),
						α/3]
	end
	return H0t
end
#Differential equation function of the evolution operator.
function dU!(dU, U, args, t)
	Omega = args[1]
	λ = args[2]
	H = args[3]
	Omega_t = args[4]

	Omega_t = Omega(t)
	#Hamiltonian
	H .= zero(Complex{typeof(t)})
	for i = 1: 9
		H .= H .+ Omega_t[i]*λ[i]
	end
	H .= -1im*H
	mul!(dU, H, U)
end
function dU(U, args, t)
	Ω = args[1]
	λ = args[2]

	Ωt = Ω(t)
	H = vector2matrix(Ωt, λ)
	dU = -1im*H*U
	return dU
end

function solveU2(H, tf, λ; type::DataType = Float64)
	U0 = Matrix{Complex{type}}(I, 3, 3)
	tspan = (type(0), tf)
	args = (H, λ, zeros(Complex{type}, 3, 3), zeros(9))
	prob = DE.ODEProblem(dU!, U0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-11, abstol=1e-14)
	if sol.retcode == :Success
		return sol
	else
		error("Could not obtain the evolution operator. Numerical integration was
			  unsuccessful.")
	end
end
function solveU(H, tf, λ; type::DataType = Float64)
	U0 = SMatrix{3,3}(Matrix{Complex{type}}(I, 3, 3))
	tspan = (zero(tf), tf)
	args = (H, λ)
	prob = DE.ODEProblem(dU, U0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-11, abstol=1e-14)
	if sol.retcode == :Success
		return sol
	else
		error("Could not obtain the evolution operator. Numerical integration was
			  unsuccessful.")
	end
end
"""
	comm_ijk!(commAB, A, B, i, j, k, f_ijk)
.
"""
function comm_ijk!(commAB, A, B, i, j, k, f_ijk)
	commAB[k] = commAB[k] + (A[i]*B[j] - A[j]*B[i])*2*f_ijk
	commAB[j] = commAB[j] + (A[k]*B[i] - A[i]*B[k])*2*f_ijk
	commAB[i] = commAB[i] + (A[j]*B[k] - A[k]*B[j])*2*f_ijk
	return nothing
end
"""
	comm!(commAB, A, B)
Compute the commutator of A and B. The result is stored in commAB.
"""
# function comm!(commAB, A, B)
# 	type = eltype(A)
# 	commAB .= zero(type)
# 	w1, w2 = type(1)/2, sqrt(type(3))/2
# 	comm_ijk!(commAB, A, B, 1, 2, 3, 1)
# 	comm_ijk!(commAB, A, B, 1, 4, 7, w1)
# 	comm_ijk!(commAB, A, B, 1, 6, 5, w1)
# 	comm_ijk!(commAB, A, B, 2, 4, 6, w1)
# 	comm_ijk!(commAB, A, B, 2, 5, 7, w1)
# 	comm_ijk!(commAB, A, B, 3, 4, 5, w1)
# 	comm_ijk!(commAB, A, B, 3, 7, 6, w1)
# 	comm_ijk!(commAB, A, B, 4, 5, 8, w2)
# 	comm_ijk!(commAB, A, B, 6, 7, 8, w2)
# 	return nothing
# end
const index_comm = ((1, 2, 3), (1, 4, 7), (1, 6, 5),
					(2, 4, 6), (2, 5, 7), (3, 4, 5),
					(3, 7, 6), (4, 5, 8), (6, 7, 8))

type = Float64
f = zeros(type, 9)
f[1] = 2*type(1)
f[2: 7] .= type(1)
f[8: 9] .= sqrt(type(3))
const struct_const = Tuple(f)
function comm!(commAB, A, B)
	commAB .= zero(eltype(A))
	for n = 1: 9
		i, j, k = index_comm[n]
		@inbounds commAB[k] += (A[i]*B[j] - A[j]*B[i])*struct_const[n]
		@inbounds commAB[j] += (A[k]*B[i] - A[i]*B[k])*struct_const[n]
		@inbounds commAB[i] += (A[j]*B[k] - A[k]*B[j])*struct_const[n]
	end
	return nothing
end
"""
	V_I(t, α, η, f, c; cr_terms::Bool = false)
Hamiltonian in the interaction picture with respect to H0.
# Arguments
- `t::Real`: time.
- `α::Real`: anarmonicity.
- `η::Real`: parameter that characterizes the strength of the 1<->2 transition relative
to the 0<->1 transition.
- `f::Function`: pulse envelope as a sunction of t.
- `c::Real`: integral of the pulse envelope as a function of t.
- `cr_terms::Bool = false`: If true, counter rotating terms are taken in account.
"""
function V_I(t, α, η, f, c; cr_terms::Bool = false)
	type = typeof(t)
	ft = f(t)/2		#note the division by 2!!!
	ct = c(t)
	cosc, sinc = cos(ct), sin(ct)
	cosα, sinα = cos(α*t), sin(α*t)
	ηf = η*ft

	if cr_terms == true
		cos_dp2, sin_dp2 = cos((α + 2)*t), sin((α + 2)*t)
		V_It = @SVector [ft*cos(2*t),
						 -ft*sin(2*t)*cos(2*ct),
						 ft*sin(2*t)*sin(2*ct),
						 ηf*cosc*(cosα + cos_dp2),
						 ηf*cosc*(sinα + sin_dp2),
						 ηf*sinc*(sinα + sin_dp2),
						 -ηf*sinc*(cosα + cos_dp2),
						 0,
						 0]
	else
		V_It = @SVector [0,
						 0,
						 0,
						 ηf*cosc*cosα,
						 ηf*cosc*sinα,
						 ηf*sinc*sinα,
						 -ηf*sinc*cosα,
						 0,
						 0]
	end
	return V_It
end
"""
	U(system_params, λ; ideal_dynamics=false)
Evolution operator.
# Arguments
- `system_params::Dictionary`: dictionary with the parameters of the system.
- `λ::Array{Array{<:Complex, 2}}`: Array with the Gell-Mann
"""
function U(system_params, λ; ideal_dynamics::Bool = false)
	tf = system_params["tf"]
	α = system_params["alpha"]
	η = system_params["eta"]
	pulse = system_params["pulse"]

	if ideal_dynamics == true
		η = zero(η)
	end
	H = @closure t-> H0(t, α, η, pulse)
	sol = solveU2(H, tf, λ)
	return sol(tf)
end
#Drag correction of the leakage.
function correction_drag(t, α, η, pulse, pulse_derivative; high_order::Bool = false)
	f = pulse(t)
	f2 = f*f
	df = pulse_derivative(t)
	η2 = η*η

	if high_order == true
		η4 = η2*η2
		f3, f4 = f2*f, f2*f2
		f5 = f4*f
		Δ = (η2 - 4)*f2/(4*α) - (η4 - 7*η2 + 12)*f4/(16*α^3)
		gx = (η2 - 4)*f3/(8*α^2) - (13*η4 - 76*η2 + 112)*f5/(128*α^4)
		gy = -df/α + 33*(η2 - 2)*f2*df/(24*α^3)
	else
		Δ = (η2 - 4)*f2/(4*α)
		gx = zero(t)
		gy = -df/α
	end
	return Δ, gx, gy
end
#Correction Hamiltonian for drag.
function Wdrag(t, α, η, pulse, pulse_derivative; high_order::Bool = false)
	type = typeof(t)
	Δt, gx, gy = correction_drag(t, α, η, pulse, pulse_derivative; high_order=high_order)
	cost, sint = cos(t), sin(t)
	gt = gx*cost + gy*sint
	ηg = η*gt
	Δprime = Δt/2

	Wt = @SVector [gt*cost,
				   -gt*sint,
				   Δprime,
				   ηg*cost,
				   ηg*sint,
				   0,
				   0,
				   -(sqrt(type(3))*Δprime),
				   Δprime]
	return Wt
end

function Udrag(system_params, λ; high_order::Bool = false)
	tf = system_params["tf"]
	α = system_params["alpha"]
	η = system_params["eta"]
	pulse = system_params["pulse"]
	pulse_derivative = system_params["pulse_derivative"]

	H_drag = @closure (
		t-> H0(t, α, η, pulse)
		+ Wdrag(t, α, η, pulse, pulse_derivative, high_order= high_order))

	sol = solveU(H_drag, tf, λ)
	return sol(tf)
end
#Lindbladian
function L(ρ, c)
	c_dagger = adjoint(c)
	n = c_dagger*c
	return c*ρ*c_dagger - (n*ρ + ρ*n)/2
end

function dρ(ρ, args, t)
	H = args[1]
	c1 = args[2]
	c2 = args[3]
	λ = args[4]

	Ht = vector2matrix(H(t), λ)
	return 1im*(ρ*Ht - Ht*ρ) + L(ρ, c1) + L(ρ, c2)
end

function solveρ(system_params, λ, c1, c2, ρ0; type::DataType = Float64)
	tf = system_params["tf"]
	α = system_params["alpha"]
	η = system_params["eta"]
	pulse = system_params["pulse"]

	H = @closure t-> H0(t, α, η, pulse)

	tspan = (zero(tf), tf)
	args = (H, c1, c2, λ)
	prob = DE.ODEProblem(dρ, ρ0, tspan, args)
	sol = DE.solve(prob, DE.Vern9(), saveat=[tf], reltol=1e-10, abstol=1e-14)
	if sol.retcode == :Success
		return sol
	else
		error("Could not solve the master equation.")
	end
end

function ρbloch(a)
	u, v, w = a
	ρ11, ρ12 = (1 + w)/2, (u - 1im*v)/2
	ρ21, ρ22 = (u + 1im*v)/2, (1 - w)/2
	ρ = @SMatrix [ρ11 ρ12 0; ρ21 ρ22 0; 0 0 0]
	return ρ
end

function average_fidelity(system_params, λ, c1, c2)
	tf = system_params["tf"]
	Uideal = U(system_params, λ; ideal_dynamics = true)

	dict = Dict("+x"=> (1.0, 0.0, 0.0), "-x"=> (-1.0, 0.0, 0.0),
				"+y"=> (0.0, 1.0, 0.0), "-y"=> (0.0, -1.0, 0.0),
				"+z"=> (0.0, 0.0, 1.0), "-z"=> (0.0, 0.0, -1.0))
	F = 0.0
	for tag in ("+x", "-x", "+y", "-y", "+z", "-z")
		a = dict[tag]
		ρ0 = ρbloch(a)
		ρf = solveρ(system_params, λ, c1, c2, ρ0)(tf)
		ρfideal = Uideal*ρ0*adjoint(Uideal)
		F += tr(ρfideal*ρf)/6
	end
	if imag(F) > 10.0^(-10.0)
		error("Complex fidelity!")
	else
		return real(F)
	end
end

max_stars_repo_name: thalesfr/MagnusCorrection

max_issues_repo_name: thalesfr/MagnusCorrection

max_forks_repo_name: thalesfr/MagnusCorrection


text:
localized = true
default_locale = "la"

include("loremla.jl")

    word()=word_list[rand(1:length(word_list))]

    words(nb=3)= [ (word_list[rand(1:length(word_list))]) for i=1:nb ]

    sentence( nb_words=6, variable_nb_words=true)= ( nb_words <= 0 && (return " ");
        variable_nb_words && (nb_words = rand(1:nb_words));
        palabras = words(nb_words);
        palabras[1] = ucfirst(palabras[1]);
        sal = join(palabras," ");
        sal*= ".";
        sal)

    sentences( nb=3)= [sentence() for i=0:nb]

    paragraph( nb_sentences=3, variable_nb_sentences=true) =
       (nb_sentences <= 0 && (return " ");
        variable_nb_sentences && (nb_sentences = rand(1:nb_sentences));
        sal = join(sentences(nb_sentences)," ");
        sal*= ".";
        sal)

    paragraphs( nb=3)= [paragraph() for i=1:nb]

    text( max_nb_chars=200)=
      (
        texto = [];
        max_nb_chars < 5 && (return "text() can only generate text of at least 5 characters");
        if max_nb_chars < 25
                size = 0
                while size < max_nb_chars
                    palabra = ( size>0? " " : "") * word();
                    push!(texto, palabra)
                    size += length(palabra)
                end
            texto[1] = ucfirst(texto[1]);
            texto[end] *= "."

        elseif max_nb_chars < 100
                size = 0
                while size < max_nb_chars
                    sentencia = (size>0? " " : "") * sentence()
                    push!(texto, sentencia)
                    size += length(sentencia)
                end
        else
                size = 0
                while size < max_nb_chars
                    parrafo = (size>0? " " : "") * paragraph()
                    push!(texto, parrafo)
                    size += length(parrafo)
                end
        end ;

        join(texto,"") )


text:
"""
Original Kings law (1910) represents the relation between voltage and velocity
in a hotwire anemometer.
The law is given by:

`E^2 = A + B*U^0.5`

This function estimates `A` and `B`.
"""
linear_king_fit(E, U) = linear_fit(sqrt.(U), E .* E)

"Type that represents a Linear (original) King's law"
struct LinearKingFit <: LeastSquares
    coefs::Array{Float64,1}
    LinearKingFit(coefs) = new(copy(coefs))
end
LinearKingFit(E,U) = LinearKingFit(linear_king_fit(E, U))

(f::LinearKingFit)(E) = ( (E.*E .- f.coefs[1]) ./ f.coefs[2] ) .^ 2

"Equation that computes the error of the modified King's law "
kingfun(x, a) = a[1] + a[2] * x[2] ^ a[3] - x[1]*x[1]

"""
Uses nonlinear least squares to fit the modified King's law:

`E^2 = A + B*U^n`

The Original (linear) King's law is used to estimate `A` and `B` when `n=1/2`.
This initial value is used as an initial guess for fitting the nonlinear modified King's law
using the function `nonlinear_fit`.
"""
function king_fit(E, U, eps=1e-8, maxiter=200)

    f = linear_king_fit(E, U)
    a0 = [f[1], f[2], 0.5]

    coefs, converged, niter = nonlinear_fit(hcat(E,U), kingfun, a0, eps, maxiter)
    return coefs
end

"Type that represents the modified King's law "
struct KingFit <: LeastSquares
    coefs::Array{Float64,1}
    KingFit(coefs) = new(copy(coefs))
end
KingFit(E,U, eps=1e-8, maxiter=200) = KingFit(king_fit(E, U, eps, maxiter))

(f::KingFit)(E) = ( (E*E - f.coefs[1]) / f.coefs[2]) ^ (1 / f.coefs[3])


text:
module ActiveSetMethods

using LinearAlgebra
using SparseArrays
using Printf
using Logging

import MathOptInterface

const MOI = MathOptInterface
const MOIU = MathOptInterface.Utilities

include("status.jl")
include("parameters.jl")
include("model.jl")

include("algorithms.jl")

include("MOI_wrapper.jl")

end

max_stars_repo_name: exanauts/ActiveSetMethods

max_issues_repo_name: exanauts/ActiveSetMethods

max_forks_repo_name: exanauts/ActiveSetMethods


text:
using Dispersal, Test, Unitful, DimensionalData
using Dispersal: layer, cyclic_index, precalc_timeindex, timestep
using DynamicGrids: Extent
using Unitful: d

suitlayer1 = [0.1 0.2;
              0.3 0.4]

@testset "layer returns the layer value" begin
    @test layer(suitlayer1, (1, 1), 34325) == 0.1
    @test layer(suitlayer1, (2, 2), 7685) == 0.4
end

@testset "sequence cycling" begin
    @test cyclic_index(-4, 2) == 2
    @test cyclic_index(-3, 2) == 1
    @test cyclic_index(-2, 2) == 2
    @test cyclic_index(-1, 2) == 1
    @test cyclic_index(0, 2) == 2
    @test cyclic_index(1, 2) == 1
    @test cyclic_index(2, 2) == 2
    @test cyclic_index(3, 2) == 1
    @test cyclic_index(4, 2) == 2
    @test cyclic_index(20, 10) == 10
    @test cyclic_index(21, 10) == 1
    @test cyclic_index(27, 10) == 7
end

a = cat([0.1 0.2; 0.3 0.4], [1.5 1.6; 1.7 1.8]; dims=3)
dimz = X(1:2), Y(1:2), Ti(0d:10d:10d)
suitseq = DimensionalArray(a, dimz)
init = zero(suitseq)

@testset "sequence of layers" begin
    rule = ExactExponentialGrowth(; intrinsicrate=1.0)
    ruleset = Ruleset(rule; timestep=1d)
    data = DynamicGrids.SimData(Extent(init=init, aux=(suit=suitseq,), tspan=1:1), ruleset)

    @testset "layers sequences are interpolated over timespans" begin
        ind = ((1, 1), (1, 2), (2, 1), (2, 2))
        times = (0d, 20d), (16d, 14d), (19d, 11d), (5d, 45d)#, (15d, 55d)

        for (t1, t2) in times
            interp1 = precalc_timeindex(suitseq, rule, data, t1)
            interp2 = precalc_timeindex(suitseq, rule, data, t2)
            @test interp1 == interp2
        end
    end
    @testset "layers returns first frame values at 0.5 through the timespan" begin
        t = 5d
        interp = precalc_timeindex(suitseq, rule, data, t)
        @test layer(suitseq, (1, 1), interp) == 0.1
        @test layer(suitseq, (1, 2), interp) == 0.2
        @test layer(suitseq, (2, 1), interp) == 0.3
        @test layer(suitseq, (2, 2), interp) == 0.4
    end

    @testset "layers returns second frame values at 1.5 times through the timespan" begin
        t = 15d
        interp = precalc_timeindex(suitseq, rule, data, t)
        @test layer(suitseq, (1, 1), interp) == 1.5
        @test layer(suitseq, (1, 2), interp) == 1.6
        @test layer(suitseq, (2, 1), interp) == 1.7
        @test layer(suitseq, (2, 2), interp) == 1.8
    end

    @testset "layers returns the second frame one timestep in" begin
        t = 10d
        interp = precalc_timeindex(suitseq, rule, data, t)
        @test layer(suitseq, (1, 1), interp) ≈ 1.5
        @test layer(suitseq, (1, 2), interp) ≈ 1.6
        @test layer(suitseq, (2, 1), interp) ≈ 1.7
        @test layer(suitseq, (2, 2), interp) ≈ 1.8
    end

end

@testset "start and stop times" begin
    A = rand(4, 5, 10)
    @test Dispersal.timestep(A) == 1
    @test Dispersal.starttime(A) == 1
    @test Dispersal.stoptime(A) == 10
    A = DimensionalArray(rand(4, 5, 10), (X, Y, Ti(20:10:110)))
    @test Dispersal.timestep(A) == 10
    @test Dispersal.starttime(A) == 20
    @test Dispersal.stoptime(A) == 110
end

@testset "LayerCopy" begin
    init = [0 0]
    l = cat([1 2], [3 4]; dims=3)
    @test LayerCopy(:l, 1) === LayerCopy(layerkey=:l)
    @test LayerCopy{:a,:b}(layerkey=:l, timeindex=2) === LayerCopy{:a,:b}(:l, 2)
    ruleset = Ruleset(LayerCopy(layerkey=:l))
    output = ArrayOutput(init; tspan=1:3, aux=(l=l,))
    sim!(output, ruleset)
    @test output == [[0 0], [3 4], [1 2]]
end

@testset "Layer on intrinsic growth rate" begin

    popSizeInit = [ 1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

    intrinsicRate = cat([ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0],
                        [ 2.0 2.0 2.0;
                        2.0 2.0 2.0;
                        2.0 2.0 2.0],
                        [ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0]; dims=3)

    popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(intrinsicRate=intrinsicRate,));
    growthRule = Ruleset(DiscreteGrowthMap(layerkey=:intrinsicRate));
    sim!(popSizeGrids, growthRule);
    @test popSizeGrids[1] == intrinsicRate[:,:,1] .* popSizeInit
    @test popSizeGrids[2] == intrinsicRate[:,:,2] .* popSizeGrids[1]
    @test popSizeGrids[3] == intrinsicRate[:,:,3] .* popSizeGrids[2]
    @test popSizeGrids[4] == intrinsicRate[:,:,1] .* popSizeGrids[3]
    @test popSizeGrids[5] == intrinsicRate[:,:,2] .* popSizeGrids[4]
    @test popSizeGrids[6] == intrinsicRate[:,:,3] .* popSizeGrids[5]

    popSizeInit = [ 1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

    intrinsicRate = cat([ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0],
                        [ 2.0 2.0 2.0;
                        2.0 2.0 2.0;
                        2.0 2.0 2.0],
                        [ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0]; dims=3)

    popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(intrinsicRate=intrinsicRate,));
    growthRule = Ruleset(ExactExponentialGrowthMap(layerkey=:intrinsicRate));
    sim!(popSizeGrids, growthRule);
    @test popSizeGrids[1] == popSizeInit
    @test popSizeGrids[2] == popSizeGrids[1].*exp.(intrinsicRate[:,:,2])
    @test popSizeGrids[3] == popSizeGrids[2].*exp.(intrinsicRate[:,:,3])
    @test popSizeGrids[4] == popSizeGrids[3].*exp.(intrinsicRate[:,:,1])
    @test popSizeGrids[5] == popSizeGrids[4].*exp.(intrinsicRate[:,:,2])
    @test popSizeGrids[6] == popSizeGrids[5].*exp.(intrinsicRate[:,:,3])

end

@testset "Test double layer on intrinsic growth rate and carying capacity" begin
    popSizeInit = [ 1.0 4.0 7.0;
                    2.0 5.0 8.0;
                    3.0 6.0 9.0]

    intrinsicRate = cat([ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0],
                        [ 2.0 2.0 2.0;
                        2.0 2.0 2.0;
                        2.0 2.0 2.0],
                        [ 1.0 1.0 1.0;
                        1.0 1.0 1.0;
                        1.0 1.0 1.0]; dims=3)

    carryingCapacity = cat([ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0],
                            [ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0],
                            [ 10.0 10.0 10.0;
                            10.0 10.0 10.0;
                            10.0 10.0 10.0]; dims=3)

    popParameter = cat(intrinsicRate, carryingCapacity; dims = 4)

    popSizeGrids = ArrayOutput(popSizeInit; tspan=1:6, aux=(popParameter=popParameter,));
    growthRule = Ruleset(ExactLogisticGrowthMap2(layerkey=:popParameter));
    sim!(popSizeGrids, growthRule);
end