text:
function get_time_series(
    psi_container::PSIContainer,
    device::PSY.Device,
    forecast_label::String,
)
    initial_time = model_initial_time(psi_container)
    @debug initial_time
    use_forecast_data = model_uses_forecasts(psi_container)
    time_steps = model_time_steps(psi_container)
    @debug "device $(PSY.get_name(device)) $(PSY.has_forecasts(device)) forecast"
    if use_forecast_data & PSY.has_forecasts(device)
        forecast = PSY.get_forecast(
            PSY.Deterministic,
            device,
            initial_time,
            forecast_label,
            length(time_steps),
        )
        return ts_vector = TS.values(PSY.get_data(forecast))
    else
        return ts_vector = ones(time_steps[end])
    end
end

function get_time_series(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{<:PSY.ElectricLoad},
    model::DeviceModel,
    get_constraint_values::Function,
)
    initial_time = model_initial_time(psi_container)
    @debug initial_time
    use_forecast_data = model_uses_forecasts(psi_container)
    time_steps = model_time_steps(psi_container)
    device_total = length(devices)

    constraint_data = Vector{DeviceRange}()
    ts_data_active = Vector{DeviceTimeSeries}()
    ts_data_reactive = Vector{DeviceTimeSeries}()

    for device in devices
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)

        if use_forecast_data
            active_power = PSY.get_maxactivepower(device)
            reactive_power = PSY.get_maxreactivepower(device)
            forecast = PSY.get_forecast(
                PSY.Deterministic,
                device,
                initial_time,
                "get_maxactivepower",
                length(time_steps),
            )
            ts_vector = TS.values(PSY.get_data(forecast))
        else
            active_power = PSY.get_activepower(device)
            reactive_power = PSY.get_reactivepower(device)
            ts_vector = ones(time_steps[end])
        end
        range_data = DeviceRange(name, get_constraint_values(device))
        add_device_services!(range_data, device, model)
        push!(constraint_data, range_data)
        push!(
            ts_data_active,
            DeviceTimeSeries(name, bus_number, active_power, ts_vector, range_data),
        )
        push!(
            ts_data_reactive,
            DeviceTimeSeries(name, bus_number, reactive_power, ts_vector, range_data),
        )

    end

    return ts_data_active, ts_data_reactive, constraint_data

end

function get_time_series(
    psi_container::PSIContainer,
    devices::IS.FlattenIteratorWrapper{<:PSY.HydroGen},
    model::DeviceModel,
    get_constraint_values::Function,
)
    initial_time = model_initial_time(psi_container)
    @debug initial_time
    use_forecast_data = model_uses_forecasts(psi_container)
    parameters = model_has_parameters(psi_container)
    time_steps = model_time_steps(psi_container)

    constraint_data = Vector{DeviceRange}()
    active_timeseries = Vector{DeviceTimeSeries}()

    for device in devices
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        if use_forecast_data
            active_power = PSY.get_rating(device)
            ts_vector = TS.values(PSY.get_data(PSY.get_forecast(
                PSY.Deterministic,
                device,
                initial_time,
                "get_rating",
                length(time_steps),
            )))
        else
            active_power = PSY.get_activepower(device)
            ts_vector = ones(time_steps[end])
        end
        range_data = DeviceRange(name, get_constraint_values(device))
        add_device_services!(range_data, device, model)
        push!(constraint_data, range_data)
        push!(
            active_timeseries,
            DeviceTimeSeries(name, bus_number, active_power, ts_vector, range_data),
        )
    end
    return active_timeseries, constraint_data
end


text:
using DrWatson
using DataFrames
using CSV
using BSON
using ValueHistories

"""
	function resize_images(images, isize::Int, channels=1)

resizes whole batch of images.
"""
function resize_images(images, isize::Int, channels=1) """ zjisiti jestli to takhle vlastně jde"""
	N = size(images, 4)
	channels_list = []
	for ch=1:channels
		push!(channels_list, reshape(imresize(images[:,:,ch,:], isize, isize, N), (isize, isize, 1, N)))
	end
	return cat(channels_list..., dims=3)
end

function scale_to_interval(X; range=[0,1])
	X_std = (X .- minimum(X))./(maximum(X)-minimum(X))
	return X_std .* (range[2]-range[1]) .+ range[1]
end

"""
	function preprocess_images(data, parameters)

Preprocess image data for ganomaly.
"""
function preprocess_images(data, parameters; range=[-1,1], denominator=16)
	(X_train, y_train), (X_val, y_val), (X_test, y_test) = data
	isize = maximum([size(X_train,1),size(X_train,2)]) # there is already fixed isize in parameters
	in_ch = parameters.in_ch

	residue = isize % denominator
	if residue != 0
		isize = isize + denominator - residue
		X_train = scale_to_interval(resize_images(X_train, isize, in_ch), range=range)
		X_val = scale_to_interval(resize_images(X_val, isize, in_ch), range=range)
		X_test = scale_to_interval(resize_images(X_test, isize, in_ch), range=range)
	end

	return (X_train, y_train), (X_val, y_val), (X_test, y_test)
end

function GANomalyHistory()
	history = Dict(
		"generator_loss" => Array{Float32}([]),
		"adversarial_loss" => Array{Float32}([]),
		"contextual_loss" => Array{Float32}([]),
		"encoder/latent_loss" => Array{Float32}([]),
		"discriminator_loss" => Array{Float32}([]),
		"val_generator_loss" => Array{Float32}([]),
		"val_discriminator_loss" => Array{Float32}([])
		)
	return history
end

"""
	function update_history(history, gl, dl, vgl=nothing, vdl=nothing)

do logging losses into history
"""
function update_history(history, gl, dl)
	push!(history["generator_loss"], gl[1])
	push!(history["adversarial_loss"], gl[2])
	push!(history["contextual_loss"], gl[3])
	push!(history["encoder/latent_loss"], gl[4])
	push!(history["discriminator_loss"], dl)
	return history
end

function update_val_history(history, vgl, vdl)
	push!(history["val_generator_loss"], vgl)
	push!(history["val_discriminator_loss"], vdl)
	return history
end

"""
	function prepare_dataloaders(data; batch_size::Int=64, iters::Int=10000)

Extracts normal data from validation dataset and returns training MLDataPattern.RandomBatches and
validation Flux.Data.DataLoader.
"""
function prepare_dataloaders(data; batch_size::Int=64, iters::Int=10000)
	train_loader = MLDataPattern.RandomBatches(data[1][1], size=batch_size, count=iters)
	# for cheching convergence I need to drop anomal samples from validation data
	val_data_ind = findall(x->x==0, data[2][2])
	if length(size(data[2][1]))==4
		val_data = data[2][1][:,:,:,val_data_ind]
	else
		val_data = data[2][1][:,val_data_ind]
	end
	val_loader = Flux.Data.DataLoader(val_data, batchsize=batch_size)
	return train_loader, val_loader
end

"""
	function prepare_dataloaders(data, params)

Extracts normal data from validation dataset and returns training MLDataPattern.RandomBatches and
validation Flux.Data.DataLoader.
"""
function prepare_dataloaders(data, params)
	return prepare_dataloaders(data, batch_size=params.batch_size, iters=params.iters)
end

"""
	clip_weights!(ps::Flux.Zygote.Params,c::Real)
	clip_weights!(ps::Flux.Zygote.Params,low::Real,high::Real)

Clips weights so they lie in the interval (-c,c)/(low,high).
"""
function clip_weights!(ps::Flux.Zygote.Params,low::Real,high::Real)
	@assert low <= high
	for p in ps
		T = eltype(p)
		p .= max.(p, T(low))
		p .= min.(p, T(high))
	end
end
clip_weights!(ps::Flux.Zygote.Params,c::Real) = clip_weights!(ps,-abs(c),abs(c))

"""
	function check_scaling(idims, scalings)

Function checks if input is compatible with scaling. Works with our nn_builders: conv_dcoder, conv_encoder.
"""
function check_scaling(idims, scalings)
    try
        ho = idims[1]/(reduce(*, scalings)) # height before reshaping
        wo = idims[2]/(reduce(*, scalings)) # width before reshaping
        # this ensures size compatibility, ho/wo are Ints
        (ho == floor(Int, ho)) ? ho = floor(Int, ho) : error("your input size and scaling is not compatible")
        (wo == floor(Int, wo)) ? wo = floor(Int, wo) : error("your input size and scaling is not compatible")
        return true
    catch e
        println(e)
        return false
    end
end

ext: jl


text:
using BinaryBuilder

name = "SuiteSparse32"
version = v"5.4.0"

# Collection of sources required to build SuiteSparse
sources = [
    ArchiveSource("https://github.com/DrTimothyAldenDavis/SuiteSparse/archive/v$(version).tar.gz",
                  "d9d62d539410d66550d0b795503a556830831f50087723cb191a030525eda770"),
]

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

# Disable OpenMP as it will probably interfere with blas threads and Julia threads
FLAGS+=(INSTALL="${prefix}" INSTALL_LIB="${libdir}" INSTALL_INCLUDE="${prefix}/include" CFOPENMP=)

FLAGS+=(LDFLAGS="${LDFLAGS} -L${libdir}")
if [[ ${target} == *mingw32* ]]; then
    FLAGS+=(UNAME=Windows)
else
    FLAGS+=(UNAME="$(uname)")
fi

mkdir -p ${prefix}/include
make -j${nproc} -C SuiteSparse_config "${FLAGS[@]}" config
cp SuiteSparse_config/SuiteSparse_config.h ${prefix}/include
for proj in SuiteSparse_config AMD BTF COLAMD KLU; do
    make -j${nproc} -C $proj "${FLAGS[@]}" static CFOPENMP="$CFOPENMP"
    [[ -d ${proj}/Include ]] && cp ${proj}/Include/*.h ${prefix}/include
done

# Move the static libraries into place
mkdir -p $libdir
pat="*.a"
if [[ ${target} == *mingw32* ]]; then
    pat="*.lib"
fi
find . -name "$pat" -exec mv {} ${libdir} \;

install_license LICENSE.txt
"""

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

# The products that we will ensure are always built
products = [
    FileProduct(["lib/libsuitesparseconfig.a", "bin/libsuitesparseconfig.lib"],  :libsuitesparseconfig_a),
    FileProduct(["lib/libamd.a", "bin/libamd.lib"],                              :libamd_a),
    FileProduct(["lib/libbtf.a", "bin/libbtf.lib"],                              :libbtf_a),
    FileProduct(["lib/libcolamd.a", "bin/libcolamd.lib"],                        :libcolamd_a),
    FileProduct(["lib/libklu.a", "bin/libklu.lib"],                              :libklu_a),
]

# Dependencies that must be installed before this package can be built
dependencies = [
]

# Build the tarballs
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: blenessy/Yggdrasil

max_issues_repo_name: blenessy/Yggdrasil

max_forks_repo_name: blenessy/Yggdrasil


text:
# MIT License
#
# Copyright (c) 2018 Martin Biel
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to deal
# in the Software without restriction, including without limitation the rights
# to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
# copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
# OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
# SOFTWARE.

struct RawTime
    name::String
    stages::Dict{Period, Int}
    col_delims::Vector{Col}
    row_delims::Vector{Row}
end

function parse_tim(filename::String)
    # Initialize auxiliary variables
    name       = "SLP"
    mode       = :NAME
    stage      = 1
    stages     = Dict{Period, Int}()
    col_delims = Vector{Col}()
    row_delims = Vector{Row}()
    # Parse the file
    open(filename) do io
        firstline = split(readline(io))
        if Symbol(firstline[1]) == :TIME
            name = join(firstline[2:end], " ")
        else
            throw(ArgumentError("`TIME` field is expected on the first line."))
        end
        for line in eachline(io)
            if mode == END
                # Parse finished
                break
            end
            words = split(line)
            first_word = Symbol(words[1])
            if first_word in TIM_MODES
                mode = first_word
                continue
            end
            if mode == PERIODS
                push!(col_delims, Symbol(words[1]))
                push!(row_delims, Symbol(words[2]))
                stages[Symbol(words[3])] = stage
                stage += 1
            else
                throw(ArgumentError("$(mode) is not a valid tim file mode."))
            end
        end
    end
    # Return tim data
    return RawTime(name, stages, col_delims, row_delims)
end

num_stages(tim::RawTime) = length(tim.stages)


text:
using CalibrationErrors, Distributions, StatsBase
using CalibrationErrors: perform, binid

using Random

Random.seed!(1234)

@testset "Binning indices" begin
    @testset "One dimension" begin
        @test_throws ArgumentError binid(UniformBinning(10), -0.5)
        @test binid(UniformBinning(10), 0) == 1
        @test binid(UniformBinning(10), 0.1) == 1
        @test binid(UniformBinning(10), 1) == 10
        @test_throws ArgumentError binid(UniformBinning(10), 1.5)
    end

    @testset "Two dimensions" begin
        @test_throws ArgumentError binid(UniformBinning(10), [-0.5, 0.5])
        @test binid(UniformBinning(10), [0, 0]) == 1
        @test binid(UniformBinning(10), [0.1, 0]) == 1
        @test binid(UniformBinning(10), [1, 1]) == 100
        @test_throws ArgumentError binid(UniformBinning(10), [1.5, 0.5])
    end
end

@testset "Simple example ($nclasses classes)" for nclasses in (2, 10, 100)
    # sample predictions and labels
    nsamples = 1_000
    predictions = rand(Dirichlet(nclasses, 1), nsamples)
    labels = rand(1:nclasses, nsamples)

    # bin data in bins of uniform width
    bins = perform(UniformBinning(10), predictions, labels)

    # check all bins
    for (idx, bin) in bins
        idxs = Int[]
        for i in axes(predictions, 2)
            idx == prod(max(1, ceil(10 * p)) for p in view(predictions, :, i)) &&
                push!(idxs, i)
        end

        @test bin.nsamples[] == length(idxs)
        @test bin.prediction_sum ≈ vec(sum(predictions[:, idxs], dims=2))
        @test bin.label_counts ≈ counts(labels[idxs], 1:nclasses)
    end
end


text:
module NumNN

    include("includes.jl")

end # module NumNN


text:
## io.jl  Some functions for reading/writing GMMs.
## (c) 2013--2014 David A. van Leeuwen

## This code is for exchange with our octave / matlab based system

## save a single GMM
function FileIO.save(filename::AbstractString, name::AbstractString, gmm::GMM)
    jldopen(filename, "w") do file
        # addrequire(file, GaussianMixtures)
        write(file, name, gmm)
    end
end
## save multiple GMMs
function FileIO.save(filename::AbstractString, name::AbstractString, gmms::Array{GMM})
    jldopen(filename, "w") do file
        # addrequire(file, GaussianMixtures)
        write(file, name, gmms)
    end
end


text:
module POLOSolvers

using MacroTools
using MacroTools: @q
using POLO
using POLO: AbstractLoss, AbstractTermination, AbstractLogger, ExecutionPolicy
using POLO.Execution: Serial, ParameterServerOptions, Master, Worker, Scheduler

macro define_polo_algorithm(algname, execution, paramdefs, ctypes)
    params = map(paramdef -> begin @capture(paramdef, param_Symbol::type_ = val_); param end, paramdefs)
    quote
        function $(esc(execution))(x₀::AbstractVector,loss::AbstractLoss,termination::AbstractTermination,logger::AbstractLogger; $(paramdefs...))
            xbegin = pointer(x₀, 1)
            xend = pointer(x₀, length(x₀) + 1)
            loss_c = @cfunction(POLO.loss_wrapper, Cdouble,
                                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
            termination_c = @cfunction(POLO.termination_wrapper, Cint,
                                       (Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
            log_c = @cfunction(POLO.log_wrapper, Nothing,
                               (Cint, Cdouble, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
            fval = ccall(POLO.$(execution), Cdouble,
                         (Ptr{Cdouble}, Ptr{Cdouble},
                          Ptr{Nothing}, Any,
                          Ptr{Nothing}, Any,
                          Ptr{Nothing}, Any,
                          $(ctypes.args...)),
                         xbegin, xend,
                         loss_c, loss,
                         termination_c, termination,
                         log_c, logger,
                         $(params...))
            return fval
        end
    end
end
macro define_inconsistent_polo_algorithm(algname, execution, paramdefs, ctypes)
    params = map(paramdef -> begin @capture(paramdef, param_Symbol::type_ = val_); param end, paramdefs)
    quote
        function $(esc(execution))(x₀::AbstractVector,loss::AbstractLoss; $(paramdefs...))
            xbegin = pointer(x₀, 1)
            xend = pointer(x₀, length(x₀) + 1)
            loss_c = @cfunction(POLO.loss_wrapper, Cdouble,
                                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
             fval = ccall(POLO.$(execution), Cdouble,
                         (Ptr{Cdouble}, Ptr{Cdouble},
                          Ptr{Nothing}, Any,
                          $(ctypes.args...)),
                         xbegin, xend,
                         loss_c, loss,
                         $(params...))
            return fval
        end
    end
end
macro define_polo_paramserver_algorithm(algname, execution, paramdefs, ctypes)
    params = map(paramdef -> begin @capture(paramdef, param_Symbol::type_ = val_); param end, paramdefs)
    @q begin
        function $(esc(execution))(x₀::AbstractVector,loss::AbstractLoss,poptions::Execution.ParameterServerOptions; $(paramdefs...))
            xbegin = pointer(x₀, 1)
            xend = pointer(x₀, length(x₀) + 1)
            loss_c = @cfunction(POLO.loss_wrapper, Cdouble,
                                (Ptr{Cdouble}, Ptr{Cdouble}, Ptr{Nothing}))
            fval = ccall(POLO.$(execution), Cdouble,
                         (Ptr{Cdouble}, Ptr{Cdouble},
                          Ptr{Nothing}, Any,
                          Ptr{Nothing},
                          $(ctypes.args...)),
                         xbegin, xend,
                         loss_c, loss,
                         options,
                         $(params...))
            return fval
        end
    end
end
macro define_polo_algorithms(alg)
    algsplit = splitdef(alg)
    algname = algsplit[:name]
    serial = Symbol(algname,:_s)
    consistent = Symbol(algname,:_mtc)
    inconsistent = Symbol(algname,:_mti)
    master = Symbol(algname,:_psm)
    worker = Symbol(algname,:_psw)
    scheduler = Symbol(algname,:_pss)
    structname = Symbol(titlecase(string(algname)))
    @q begin
        @define_polo_algorithm($algname,$serial,$(algsplit[:kwargs]),$(prettify(algsplit[:body])))
        @define_polo_algorithm($algname,$consistent,$(algsplit[:kwargs]),$(prettify(algsplit[:body])))
        @define_inconsistent_polo_algorithm($algname,$inconsistent,$(algsplit[:kwargs]),$(prettify(algsplit[:body])))
        @define_polo_paramserver_algorithm($algname,$master,$(algsplit[:kwargs]),$(prettify(algsplit[:body])))
        @define_polo_paramserver_algorithm($algname,$worker,$(algsplit[:kwargs]),$(prettify(algsplit[:body])))
        @define_polo_paramserver_algorithm($algname,$scheduler,$(algsplit[:kwargs]),$(prettify(algsplit[:body])))

        struct $(esc(structname)){Execution <: ExecutionPolicy} <: POLOAlgorithm
            x::Vector{Float64}

            function $(esc(structname))(execution::Execution) where Execution <: ExecutionPolicy
                return new{Execution}(Vector{Float64}())
            end
        end
        # Serial
        $(esc(structname))() = $(esc(structname))(Execution.Serial())
        function (alg::$(esc(structname)){Execution})(x₀::AbstractVector,loss::AbstractLoss,termination::AbstractTermination,logger::AbstractLogger; kw...) where Execution <: Serial
            resize!(alg.x,length(x₀))
            alg.x .= x₀
            return $serial(alg.x, loss, termination, logger; kw...)
        end
        function (alg::$(esc(structname)){Execution})(x₀::AbstractVector,loss::AbstractLoss) where Execution <: Serial
            termination = POLO.Utility.MaxIteration(100)
            return alg(alg.x,loss,termination,Utility.ProgressLogger.Value(termination))
        end
        function (alg::$(esc(structname)){Execution})(x₀::AbstractVector,loss::AbstractLoss,termination::AbstractTermination) where Execution <: Serial
            return alg(alg.x,loss,termination,Utility.ProgressLogger.Gradient(termination))
        end
        # Parameter server
        function (::$(esc(structname)){Execution})(x₀::AbstractVector,loss::AbstractLoss,poptions::ParameterServerOptions; kw...) where Execution <: Master
            return $master(x₀,loss,poptions; kw...)
        end
        function (::$(esc(structname)){Execution})(x₀::AbstractVector,loss::AbstractLoss,poptions::ParameterServerOptions; kw...) where Execution <: Worker
            return $worker(x₀,loss,poptions; kw...)
        end
        function (::$(esc(structname)){Execution})(x₀::AbstractVector,loss::AbstractLoss,poptions::ParameterServerOptions; kw...) where Execution <: Scheduler
            return $scheduler(x₀,loss,poptions; kw...)
        end
    end
end

abstract type POLOAlgorithm end
@define_polo_algorithms gradient(; γ::Float64 = 1.) = (Cdouble,)
@define_polo_algorithms momentum(; μ::Float64 = 0.9, ϵ::Float64 = 1e-1, γ::Float64 = 1.) = (Cdouble,Cdouble,Cdouble)
@define_polo_algorithms nesterov(; μ::Float64 = 0.9, ϵ::Float64 = 1e-1, γ::Float64 = 1.) = (Cdouble,Cdouble,Cdouble)
@define_polo_algorithms adagrad(; γ::Float64 = 1., ϵ::Float64 = 1e-6) = (Cdouble,Cdouble)
@define_polo_algorithms adadelta(; γ::Float64 = 1., ρ::Float64 = 0.95, ϵ::Float64 = 1e-6) = (Cdouble,Cdouble,Cdouble)
@define_polo_algorithms adam(; μ::Float64 = 0.99, ϵ::Float64 = 1e-2, γ::Float64 = 1., ρ::Float64 = 0.999, ϵ_rms::Float64 = 1e-8) = (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble)
@define_polo_algorithms nadam(; μ::Float64 = 0.99, ϵ::Float64 = 1e-2, γ::Float64 = 1., ρ::Float64 = 0.999, ϵ_rms::Float64 = 1e-8) = (Cdouble,Cdouble,Cdouble,Cdouble,Cdouble)

end


text:
# @testset "Escape Solve" begin
##
res = load(joinpath(@__DIR__, "escape_solve.jld2"))
prob,opts = Problems.DubinsCar(:escape)
solver = ALTROSolver(prob, opts, infeasible=true, R_inf=0.1)

##
ilqr = Altro.get_ilqr(solver)
let solver = ilqr
    Altro.errstate_jacobians!(solver.model, solver.G, solver.Z)
    Altro.dynamics_expansion!(solver)
    Altro.error_expansion!(solver.model, solver.D, solver.G)
    Altro.cost_expansion!(solver.obj.obj, solver.Efull, solver.Z)
    Altro.error_expansion!(solver.model, solver.Eerr, solver.Efull, solver.G, solver.Z)
end

@test Matrix.(ilqr.G)[1:prob.N] ≈ fill(Matrix(I,3,3), prob.N)
@test [D.∇f for D in ilqr.D] ≈ res["D"]
@test [E.xx for E in ilqr.Eerr] ≈ res["Q0"]
@test [E.uu for E in ilqr.Eerr] ≈ res["R0"]
@test [E.x for E in ilqr.Eerr] ≈ res["q0"]
@test [E.u for E in ilqr.Eerr] ≈ res["r0"]

conset = get_constraints(solver.solver_al)
let solver = ilqr
    Altro.evaluate_constraints!(conset, solver.Z)
    Altro.constraint_jacobians!(conset, solver.Z)
    Altro.algrad!(conset)
    Altro.alhess!(conset)
    # Altro.cost_expansion!(solver.obj.obj, solver.Efull, solver.Z)
    # Altro.error_expansion!(solver.model, solver.Eerr, solver.Efull, solver.G, solver.Z)
    # TO.cost_expansion!(solver.quad_obj, solver.obj, solver.Z, init=true, rezero=true)
    # TO.error_expansion!(solver.E, solver.quad_obj, solver.model, solver.Z, solver.G)
end
vals = [Vector.(cv.vals) for cv in conset]
jacs = [Matrix.(cv.jac) for cv in conset]
grad = [Vector.(cv.grad) for cv in conset]
hess = [Matrix.(cv.hess) for cv in conset]

@test vals[1] ≈ res["vals"][1]  # circle constraints
@test vals[2] ≈ res["vals"][2]  # bound constraint
@test vals[3] ≈ res["vals"][3]  # goal constraint
@test vals[4] ≈ res["vals"][4]  # infeasible constraint
@test [jac[:,1:3] for jac in jacs[1]] ≈ res["jacs"][1]    # circle constraint
@test jacs[2] ≈ res["jacs"][2]                            # bound constraint
@test [jac[:,1:3] for jac in jacs[3]] ≈ res["jacs"][3]    # goal constraint
@test [jac for jac in jacs[4]] ≈ res["jacs"][4]           # infeasible constraint
@test norm([jac[:,1:2] for jac in jacs[4]]) ≈ 0.0

@test [grad[1:3] for grad in grad[1]] ≈ res["grad"][1]    # circle constraint
@test grad[2] ≈ res["grad"][2]                            # bound constraint
@test [grad[1:3] for grad in grad[3]] ≈ res["grad"][3]    # goal constraint
@test [grad[4:end] for grad in grad[4]] ≈ res["grad"][4]  # infeasible constraint

@test [hess[1:3,1:3] for hess in hess[1]] ≈ res["hess"][1]      # circle constraint
@test hess[2] ≈ res["hess"][2]                                  # bound constraint
@test [hess[1:3,1:3] for hess in hess[3]] ≈ res["hess"][3]      # goal constraint
@test [hess[4:end,4:end] for hess in hess[4]] ≈ res["hess"][4]  # infeasible constraint

let solver = ilqr
    Altro.cost_expansion!(solver.obj, solver.Efull, solver.Z)
    Altro.error_expansion!(solver.model, solver.Eerr, solver.Efull, solver.G, solver.Z)
end
@test [E.xx for E in ilqr.Eerr] ≈ res["Q"]
@test [E.uu for E in ilqr.Eerr] ≈ res["R"]
@test [E.x for E in ilqr.Eerr] ≈ res["q"]
@test [E.u for E in ilqr.Eerr] ≈ res["r"]

ΔV = Altro.backwardpass!(ilqr)
@test ΔV ≈ res["dV"]

J = cost(ilqr)
J_new = Altro.forwardpass!(ilqr, J)
@test J_new ≈ res["J_new"]
# end


text:
using Documenter, ReadDatastores

makedocs(
    format = Documenter.HTML(),
    sitename = "ReadDatastores.jl",
    pages = [
        "Home"         => "index.md",
        "Datastore types" => "read-datastores.md",
        "Building & loading datastores" => "build-datastores.md",
        "Indexing & Iteration" => "indexing.md",
        "Additional Methods" => "additional-methods.md"
    ],
    authors = "Ben J. Ward, The BioJulia Organisation and other contributors."
)

deploydocs(
    repo = "github.com/BioJulia/ReadDatastores.jl.git",
    deps = nothing,
    make = nothing,
    push_preview = true
)


text:
module MySQLLoader

# https://github.com/JuliaDatabases/MySQL.jl
using MySQL # MySQL.jl v0.7.0
using Octo.Repo: ExecuteResult
using Octo.Backends: UnsupportedError

const current = Dict{Symbol, Any}(
    :conn => nothing,
)

current_conn() = current[:conn]

# db_connect
function db_connect(; kwargs...)
    args = (:hostname, :username, :password)
    (hostname, username, password) = getindex.(Ref(kwargs), args)
    options = filter(kv -> !(kv.first in args), kwargs)
    conn = MySQL.connect(hostname, username, password; options...)
    current[:conn] = conn
end

# db_disconnect
function db_disconnect()
    conn = current_conn()
    MySQL.disconnect(conn)
    current[:conn] = nothing
end

# query
function query(sql::String)
    conn = current_conn()
    table = MySQL.Query(conn, sql)
    collect(table)
end

function query(prepared::String, vals::Vector) # throw UnsupportedError
    throw(UnsupportedError("needs to be implemented"))
end

# execute
function execute(sql::String)::ExecuteResult
    conn = current_conn()
    MySQL.execute!(conn, sql)
    ExecuteResult()
end

function execute(prepared::String, vals::Vector)::ExecuteResult
    conn = current_conn()
    stmt = MySQL.Stmt(conn, prepared)
    MySQL.execute!(stmt, vals)
    ExecuteResult()
end

function execute(prepared::String, nts::Vector{<:NamedTuple})::ExecuteResult
    conn = current_conn()
    stmt = MySQL.Stmt(conn, prepared)
    MySQL.execute!(stmt, nts...)
    ExecuteResult()
end

end # module Octo.Backends.MySQLLoader


text:
"""
$SIGNATURES

Viewing a nested named tuple as a tree and return, as a tuple, the values
at the leaves, in the order they appear in the original tuple.

```julia-repl
julia> t = (X = (x = 1, y = 2), Y = 3)
julia> flat_values(t)
(1, 2, 3)
```
"""
function flat_values(params::NamedTuple)
    values = []
    for k in keys(params)
        value = getproperty(params, k)
        if value isa NamedTuple
            append!(values, flat_values(value))
        else
            push!(values, value)
        end
    end
    return Tuple(values)
end

#-#
## FOR ENCODING AND DECODING MODEL METADATA

# This method needs revisiting when new scitypes are introduced:
function encode_dic(s)
    if s isa Symbol
        return string(":", s)
    elseif s isa AbstractString
        return string(s)
    else # we have some more complicated object
        prestring = string("`", s, "`")
        # hack for objects with gensyms in their string representation:
        str = replace(prestring, '#'=>'_')
        return str
    end
end

encode_dic(v::Vector) = encode_dic.(v)
function encode_dic(d::AbstractDict)
    ret = LittleDict{}()
    for (k, v) in d
        ret[encode_dic(k)] = encode_dic(v)
    end
    return ret
end

function decode_dic(s::String)
    if !isempty(s)
        if  s[1] == ':'
            return Symbol(s[2:end])
        elseif s[1] == '`'
            return eval(Meta.parse(s[2:end-1]))
        else
            return s
        end
    else
        return ""
    end
end

decode_dic(v::Vector) = decode_dic.(v)
function decode_dic(d::AbstractDict)
    ret = LittleDict()
    for (k, v) in d
        ret[decode_dic(k)] = decode_dic(v)
    end
    return ret
end

# the inverse of a multivalued dictionary is a multivalued
# dictionary:
function inverse(d::LittleDict{S,Set{T}}) where {S,T}
    dinv = LittleDict{T,Set{S}}()
    for key in keys(d)
        for val in d[key]
            if val in keys(dinv)
                push!(dinv[val], key)
            else
                dinv[val] = Set([key,])
            end
        end
    end
    return dinv
end
function inverse(d::Dict{S,Set{T}}) where {S,T}
    dinv = Dict{T,Set{S}}()
    for key in keys(d)
        for val in d[key]
            if val in keys(dinv)
                push!(dinv[val], key)
            else
                dinv[val] = Set([key,])
            end
        end
    end
    return dinv
end
#-#

## SOME GENERAL PURPOSE MACROS

macro colon(p)
    Expr(:quote, p)
end

"""
$SIGNATURES

Take a dictionary of type `{T, <:Real}` and return the keys sorted according to
the value associated with them.

```julia-repl
julia> d = Dict("abc"=>5, "def"=>7, "ghi"=>2)
julia> keys_ordered_by_values(d)
3-element Array{String,1}:
 "ghi"
 "abc"
 "def"
```
"""
keys_ordered_by_values(d::Dict{T,<:Real} where T) = sort(collect(keys(d)), by=k->d[k])

"""
## `@curve`

The code,

    @curve var range code

evaluates `code`, replacing appearances of `var` therein with each
value in `range`. The range and corresponding evaluations are returned
as a tuple of arrays. For example,

    @curve  x 1:3 (x^2 + 1)

evaluates to

    ([1,2,3], [2, 5, 10])

This is convenient for plotting functions using, eg, the `Plots` package:

    plot(@curve x 1:3 (x^2 + 1))

A macro `@pcurve` parallelizes the same behaviour.  A two-variable
implementation is also available, operating as in the following
example:

    julia> @curve x [1,2,3] y [7,8] (x + y)
    ([1,2,3],[7 8],[8.0 9.0; 9.0 10.0; 10.0 11.0])

    julia> ans[3]
    3×2 Array{Float64,2}:
      8.0   9.0
      9.0  10.0
     10.0  11.0

N.B. The second range is returned as a *row* vector for consistency
with the output matrix. This is also helpful when plotting, as in:

    julia> u1, u2, A = @curve x range(0, stop=1, length=100) α [1,2,3] x^α
    julia> u2 = map(u2) do α "α = "*string(α) end
    julia> plot(u1, A, label=u2)

which generates three superimposed plots - of the functions x, x^2 and x^3 - each
labels with the exponents α = 1, 2, 3 in the legend.

"""
macro curve(var1, range, code)
    quote
        output = []
        N = length($(esc(range)))
        for i in eachindex($(esc(range)))
            local $(esc(var1)) = $(esc(range))[i]
            print((@colon $var1), "=", $(esc(var1)), "                   \r")
            flush(stdout)
            # print(i,"\r"); flush(stdout)
            push!(output, $(esc(code)))
        end
        collect($(esc(range))), [x for x in output]
    end
end

macro curve(var1, range1, var2, range2, code)
    quote
        output = Array{Float64}(undef, length($(esc(range1))), length($(esc(range2))))
        for i1 in eachindex($(esc(range1)))
            local $(esc(var1)) = $(esc(range1))[i1]
            for i2 in eachindex($(esc(range2)))
                local $(esc(var2)) = $(esc(range2))[i2]
                # @dbg $(esc(var1)) $(esc(var2))
                print((@colon $var1), "=", $(esc(var1)), " ")
                print((@colon $var2), "=", $(esc(var2)), "                    \r")
                flush(stdout)
                output[i1,i2] = $(esc(code))
            end
        end
        collect($(esc(range1))), permutedims(collect($(esc(range2)))), output
    end
end

macro pcurve(var1, range, code)
    quote
        N = length($(esc(range)))
        pairs = @distributed vcat for i in eachindex($(esc(range)))
            local $(esc(var1)) = $(esc(range))[i]
#            print((@colon $var1), "=", $(esc(var1)), "                    \r")
#            flush(stdout)
#            print(i,"\r"); flush(stdout)
            [( $(esc(range))[i], $(esc(code)) )]
        end
        sort!(pairs, by=first)
        collect(map(first,pairs)), collect(map(last, pairs))
    end
end

## FOR PRETTY PRINTING COLUMN TABLES

function pretty(X; showtypes=true, alignment=:l, kwargs...)
    names = schema(X).names |> collect
    if showtypes
        types = schema(X).types |> collect
        scitypes = schema(X).scitypes |> collect
        header = hcat(names, types, scitypes) |> permutedims
    else
        header  = names
    end
    show_color = MLJBase.SHOW_COLOR
    color_off()
    try
        PrettyTables.pretty_table(MLJBase.matrix(X),
                                  header; alignment=alignment, kwargs...)
    catch
        println("Trouble displaying table.")
    end
    show_color ? color_on() : color_off()
    return nothing
end

## RECURSIVE VERSIONS OF getproperty and setproperty!

# applying the following to `:(a.b.c)` returns `(:(a.b), :c)`
function reduce_nested_field(ex)
    ex.head == :. || throw(ArgumentError)
    tail = ex.args[2]
    tail isa QuoteNode || throw(ArgumentError)
    field = tail.value
    field isa Symbol || throw(ArgmentError)
    subex = ex.args[1]
    return (subex, field)
end

"""
    recursive_getproperty(object, nested_name::Expr)

Call getproperty recursively on `object` to extract the value of some
nested property, as in the following example:

    julia> object = (X = (x = 1, y = 2), Y = 3)
    julia> recursive_getproperty(object, :(X.y))
    2

"""
recursive_getproperty(obj, property::Symbol) = getproperty(obj, property)
function recursive_getproperty(obj, ex::Expr)
    subex, field = reduce_nested_field(ex)
    return recursive_getproperty(recursive_getproperty(obj, subex), field)
end

"""
    recursively_setproperty!(object, nested_name::Expr, value)

Set a nested property of an `object` to `value`, as in the following example:

```
julia> mutable struct Foo
           X
           Y
       end

julia> mutable struct Bar
           x
           y
       end

julia> object = Foo(Bar(1, 2), 3)
Foo(Bar(1, 2), 3)

julia> recursively_setproperty!(object, :(X.y), 42)
42

julia> object
Foo(Bar(1, 42), 3)
```

"""
recursive_setproperty!(obj, property::Symbol, value) =
    setproperty!(obj, property, value)
function recursive_setproperty!(obj, ex::Expr, value)
    subex, field = reduce_nested_field(ex)
    last_obj = recursive_getproperty(obj, subex)
    return recursive_setproperty!(last_obj, field, value)
end


text:
module AntisymmetricMatrices
using AutoHashEquals
#### More specialized matrix types ####

import Base.LinAlg: char_uplo, checksquare, HessenbergQ, @commutative
export isantisymmetric, Antisymmetric, AntisymTridiagonal, AntisymBlockdiagonal, pfaffian, pfaffiansign

function isantisymmetric(A::AbstractMatrix)
    n = checksquare(A)
    for i in 1:n, j in i:n
        if A[i,j]!=-ctranspose(A[j,i])
            return false
        end
    end
    return true
end

## (complex) antisymmetric tridiagonal matrices

"""
    AntisymTridiagonal(ev)

Construct a antisymmetric tridiagonal matrix from the first super-diagonal.
The result is of type `AntisymTridiagonal` and provides efficient specialized
eigensolvers, but may be converted into a regular matrix with
[`convert(Array, _)`](@ref) (or `Array(_)` for short).

# Example

```jldoctest
julia> ev = [7; 8; 9]
3-element Array{Int64,1}:
 7
 8
 9

julia> AntisymTridiagonal(ev)
4×4 AntisymTridiagonal{Int64}:
  ⋅   7   ⋅   ⋅
 -7   ⋅   8   ⋅
  ⋅  -8   ⋅   9
  ⋅   ⋅  -9   ⋅
```
"""
struct AntisymTridiagonal{T} <: AbstractMatrix{T}
    ev::Vector{T} # superdiagonal
    AntisymTridiagonal{T}(ev::Vector{T}) where T = new{T}(ev)
end

"""
    AntisymBlockdiagonal(λ)

Canonical form of an antisymmetric matrix as a vector of 2x2 antisymmetric blocks on the diagonal
followed by a number of 1x1 antisymmetric blocks (aka zeros on the diagonal).
The result is of type `AntisymBlockdiagonal` and provides efficient specialized
eigensolvers, but may be converted into a regular matrix with
[`convert(Array, _)`](@ref) (or `Array(_)` for short).

# Example

```jldoctest
julia> λ = [7; 8]
2-element Array{Int64,1}:
 7
 8

julia> AntisymBlockdiagonal(λ, 1)
5×5 AntisymBlockdiagonal{Int64}:
  ⋅   7   ⋅   ⋅   ⋅
 -7   ⋅   ⋅   ⋅   ⋅
  ⋅   ⋅   ⋅   8   ⋅
  ⋅   ⋅  -8   ⋅   ⋅
  ⋅   ⋅   ⋅   ⋅   ⋅
```
"""
struct AntisymBlockdiagonal{T} <: AbstractMatrix{T}
    λ::Vector{T}
    n::Int
    AntisymBlockdiagonal{T}(λ::Vector{T}, n::Int) where T = new{T}(λ, n)
end

AntisymTridiagonal(ev::Vector{T}) where T = AntisymTridiagonal{T}(ev)
AntisymBlockdiagonal(λ::Vector{T}, n::Int=0) where T = AntisymBlockdiagonal{T}(λ, n)

#######################################################################################################
# making the print function show the structure

function Base.replace_in_print_matrix(A::AntisymTridiagonal, i::Integer, j::Integer, s::AbstractString)
    i==j-1||i==j+1 ? s : Base.replace_with_centered_mark(s)
end

function Base.replace_in_print_matrix(A::AntisymBlockdiagonal, i::Integer, j::Integer, s::AbstractString)
    (i+1==j&&iseven(j))||(i==j+1&&iseven(i)) ? s : Base.replace_with_centered_mark(s)
end

function Base.show(io::IO, M::AntisymTridiagonal)
    println(io, summary(M), ":")
    print(io, "superdiag:")
    Base.print_matrix(io, (M.ev)')
end

function Base.show(io::IO, M::AntisymBlockdiagonal)
    println(io, summary(M), "with $(M.n) zero blocks and block values:")
    Base.print_matrix(io, (M.λ)')
end

#############################################################################################
# Creation and conversion

function AntisymTridiagonal(A::AbstractMatrix)
    if diag(A,1) == -1*diag(A,-1)
        AntisymTridiagonal(diag(A,1))
    else
        throw(ArgumentError("matrix is not antisymmetric; cannot convert to AntisymTridiagonal"))
    end
end

Base.convert(::Type{AntisymTridiagonal{T}}, S::AntisymTridiagonal) where {T} =
    AntisymTridiagonal(convert(Vector{T}, S.ev))
Base.convert(::Type{AbstractMatrix{T}}, S::AntisymTridiagonal) where {T} =
    AntisymTridiagonal(convert(Vector{T}, S.ev))
function Base.convert(::Type{Matrix{T}}, M::AntisymTridiagonal{T}) where T
    n = size(M, 1)
    Mf = zeros(T, n, n)
    @inbounds begin
        @simd for i = 1:n-1
            Mf[i+1,i] = -1*M.ev[i]
            Mf[i,i+1] = M.ev[i]
        end
    end
    return Mf
end
Base.convert(::Type{Matrix}, M::AntisymTridiagonal{T}) where {T} = convert(Matrix{T}, M)
Base.convert(::Type{Array}, M::AntisymTridiagonal) = convert(Matrix, M)
Base.full(M::AntisymTridiagonal) = convert(Array, M)

Base.convert(::Type{AntisymBlockdiagonal{T}}, S::AntisymBlockdiagonal) where {T} =
    AntisymBlockdiagonal(convert(Vector{T}, S.λ), S.n)
Base.convert(::Type{AbstractMatrix{T}}, S::AntisymBlockdiagonal) where {T} =
    AntisymBlockdiagonal(convert(Vector{T}, S.λ), S.n)
function Base.convert(::Type{Matrix{T}}, M::AntisymBlockdiagonal{T}) where T
    m = length(M.λ)
    n = M.n
    N = 2m+n
    Mf = zeros(T, N, N)
    for i = 1:m
        Mf[2*i,2*i-1] = -1*M.λ[i]
        Mf[2*i-1,2*i] = M.λ[i]
    end
    return Mf
end
Base.convert(::Type{Matrix}, M::AntisymBlockdiagonal{T}) where {T} = convert(Matrix{T}, M)
Base.convert(::Type{Array}, M::AntisymBlockdiagonal) = convert(Matrix, M)
Base.full(M::AntisymBlockdiagonal) = convert(Array, M)

Base.convert(::Type{Tridiagonal}, M::AntisymTridiagonal{T}) where T = Tridiagonal(-M.ev, zeros(T, size(M,1)), M.ev)

###################################################################################################
# Methods: size, similar, conj, copy, real, imag, transpose, ctranspose, diag, +, -, *, /, ==

Base.size(A::AntisymTridiagonal) = (length(A.ev)+1, length(A.ev)+1)
function Base.size(A::AntisymTridiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d<=2
        return length(A.ev)+1
    else
        return 1
    end
end

Base.size(A::AntisymBlockdiagonal) = (2*length(A.λ)+A.n, 2*length(A.λ)+A.n)
function Base.size(A::AntisymBlockdiagonal, d::Integer)
    if d < 1
        throw(ArgumentError("dimension must be ≥ 1, got $d"))
    elseif d<=2
        return 2*length(A.λ)+A.n
    else
        return 1
    end
end

Base.similar(S::AntisymTridiagonal, ::Type{T}) where {T} = AntisymTridiagonal{T}(similar(S.ev, T))
Base.similar(S::AntisymBlockdiagonal, ::Type{T}) where {T} = AntisymTridiagonal{T}(similar(S.λ, T), S.n)

for func in (:(Base.conj), :(Base.copy), :(Base.real), :(Base.imag))
    @eval ($func)(M::AntisymTridiagonal) = AntisymTridiagonal(($func)(M.ev))
    @eval ($func)(M::AntisymBlockdiagonal) = AntisymBlockdiagonal(($func)(M.λ), M.n)
end

Base.transpose(M::AntisymTridiagonal) = -1*M
Base.ctranspose(M::AntisymTridiagonal) = conj(-1*M)

Base.transpose(M::AntisymBlockdiagonal) = -1*M
Base.ctranspose(M::AntisymBlockdiagonal) = conj(-1*M)

function Base.diag(M::AntisymTridiagonal{T}, n::Integer=0) where T
    absn = abs(n)
    if n == 1
        return M.ev
    elseif n==-1
        return -1*M.ev
    elseif absn<size(M,1)
        return zeros(T,size(M,1)-absn)
    else
        throw(ArgumentError("$n-th diagonal of a $(size(M)) matrix doesn't exist!"))
    end
end

function Base.diag(M::AntisymBlockdiagonal{T}, n::Integer=0) where T
    absn = abs(n)
    if n == 1
        l = size(M.λ)
        d = zeros(T, 2*l+M.n-1)
        d[1:2:2l-1] = M.λ
        return d
    elseif n==-1
        l = size(M.λ)
        d = zeros(T, 2*l+M.n-1)
        d[1:2:2l-1] = -1*M.λ
        return d
    elseif absn<size(M,1)
        return zeros(T,size(M,1)-absn)
    else
        throw(ArgumentError("$n-th diagonal of a $(size(M)) matrix doesn't exist!"))
    end
end

Base.:+(A::AntisymTridiagonal, B::AntisymTridiagonal) = AntisymTridiagonal(A.ev+B.ev)
Base.:-(A::AntisymTridiagonal, B::AntisymTridiagonal) = AntisymTridiagonal(A.ev-B.ev)
@commutative Base.:*(A::AntisymTridiagonal, B::Number) = AntisymTridiagonal(A.ev*B)
Base.:/(A::AntisymTridiagonal, B::Number) = AntisymTridiagonal(A.ev/B)
Base.:(==)(A::AntisymTridiagonal, B::AntisymTridiagonal) = (A.ev==B.ev)

function Base.:+(A::AntisymBlockdiagonal, B::AntisymBlockdiagonal)
    if !(size(A, 1) == size(B, 1))
        throw(DimensionMismatch("A has first dimension $(size(S,1)), B has $(size(B,1)); must match"))
    end
    n = min(A.n, B.n)
    if n==A.n
        λ = copy(A.λ)
        λ[1:length(B.λ)] += B.λ
    else
        λ = copy(B.λ)
        λ[1:length(A.λ)] += A.λ
    end
    AntisymBlockdiagonal(λ, n)
end

Base.:-(A::AntisymBlockdiagonal, B::AntisymBlockdiagonal) = A+(-1*B)

@commutative Base.:*(A::AntisymBlockdiagonal, B::Number) = AntisymBlockdiagonal(A.λ*B, A.n)
Base.:/(A::AntisymBlockdiagonal, B::Number) = AntisymBlockdiagonal(A.λ/B, A.n)
Base.:(==)(A::AntisymBlockdiagonal, B::AntisymBlockdiagonal) = (A.λ==B.λ && A.n==B.n)

const tridiagtypes = Union{Diagonal, Bidiagonal, Tridiagonal, SymTridiagonal}
@commutative Base.:+(A::AntisymTridiagonal, B::tridiagtypes) = convert(Tridiagonal, A) + convert(Tridiagonal, B)
Base.:-(A::AntisymTridiagonal, B::tridiagtypes) = convert(Tridiagonal, A) - convert(Tridiagonal, B)
Base.:-(A::tridiagtypes, B::AntisymTridiagonal) = convert(Tridiagonal, A) - convert(Tridiagonal, B)

###################################################################################################
# Implementation of specialized matrix multiplication with A_mul_B

function Base.A_mul_B!(C::StridedVecOrMat, S::AntisymTridiagonal, B::StridedVecOrMat)
    m, n = size(B, 1), size(B, 2)
    if !(m == size(S, 1) == size(C, 1))
        throw(DimensionMismatch("A has first dimension $(size(S,1)), B has $(size(B,1)), C has $(size(C,1)) but all must match"))
    end
    if n != size(C, 2)
        throw(DimensionMismatch("second dimension of B, $n, doesn't match second dimension of C, $(size(C,2))"))
    end

    β = S.ev
    @inbounds begin
        @simd for j = 1:n
            C[1, j] = B[2, j]*β[1]
            @simd for i = 2:m - 1
                C[i, j] = B[i+1, j]*β[i] - B[i-1, j]*β[i-1]
            end
            C[m, j] = -B[m-1, j]*β[m-1]
        end
    end

    return C
end

###################################################################################################
# Implementation of specialized inv, det, and pfaffian

#Implements the inverse using the recurrence relation between principal minors
# special cased from inv_usmani in julia/base/linalg/tridiag.jl
#Reference:
#    R. Usmani, "Inversion of a tridiagonal Jacobi matrix",
#    Linear Algebra and its Applications 212-213 (1994), pp.413-414
#    doi:10.1016/0024-3795(94)90414-6

function inv_usmani(c::Vector{T}) where T
    n = length(c)+1
    θ = zeros(T, n+1) #principal minors of A
    θ[1] = 1
    n>=1 && (θ[2] = zero(T))
    for i=2:n
        θ[i+1] = (c[i-1]^2)*θ[i-1]
    end
    φ = zeros(T, n+1)
    φ[n+1] = 1
    n>=1 && (φ[n] = zero(T))
    for i=n-1:-1:1
        φ[i] = (c[i]^2)*φ[i+2]
    end
    α = Matrix{T}(n, n)
    for i=1:n, j=1:n
        sign = (i+j)%2==0 ? (+) : (-)
        if i<j
            α[i,j]=(sign)(prod(c[i:j-1]))*θ[i]*φ[j+1]/θ[n+1]
        elseif i==j
            α[i,i]=                       θ[i]*φ[i+1]/θ[n+1]
        else #i>j
            α[i,j]=(sign)(prod(-1*c[j:i-1]))*θ[j]*φ[i+1]/θ[n+1]
        end
    end
    α
end
Base.inv(A::AntisymTridiagonal) = inv_usmani(A.ev)

minv(q::Number) = -1./q

Base.inv(A::AntisymBlockdiagonal) = (A.n==0 || throw(error("Can't invert due to zero columns. Try pinv")); AntisymBlockdiagonal(minv.(A.λ)))
Base.pinv(A::AntisymBlockdiagonal) = AntisymBlockdiagonal(minv.(A.λ), A.n)

function pfaff(c::Vector{T}) where T
    n = length(c)+1
    if n == 0
        return one(T)
    elseif isodd(n)
        return zero(T)
    end
    return prod(c[1:2:n-1])
end

function pfaffsign(c::Vector{T}) where T
    n = length(c)+1
    if n == 0
        return one(T)
    elseif isodd(n)
        return zero(T)
    end
    return prod(sign.(c[1:2:n-1]))
end

pfaffian(A::AntisymTridiagonal) = pfaff(A.ev)
Base.det(A::AntisymTridiagonal) = pfaffian(A)^2
pfaffiansign(A::AntisymTridiagonal) = pfaffsign(A.ev)

pfaffian(A::AntisymBlockdiagonal{T}) where T = (A.n==0 ? prod(A.λ) : zero(T))
Base.det(A::AntisymBlockdiagonal) = pfaffian(A)^2
pfaffiansign(A::AntisymBlockdiagonal{T}) where T = (A.n==0 ? prod(sign.(A.λ)) : zero(T))

function Base.getindex(A::AntisymTridiagonal{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= size(A,1) && 1 <= j <= size(A,2))
        throw(BoundsError(A, (i,j)))
    end
    if i == j + 1
        return -A.ev[j]
    elseif i + 1 == j
        return A.ev[i]
    else
        return zero(T)
    end
end

function Base.setindex!(A::AntisymTridiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    throw(ArgumentError("cannot set entry ($i, $j)"))
    return x
end

function Base.getindex(A::AntisymBlockdiagonal{T}, i::Integer, j::Integer) where T
    if !(1 <= i <= size(A,1) && 1 <= j <= size(A,2))
        throw(BoundsError(A, (i,j)))
    end
    if i == j + 1 && iseven(i) && i <= 2*length(A.λ)
        return -A.λ[i>>1]
    elseif i + 1 == j && iseven(j) && j <= 2*length(A.λ)
        return A.λ[j>>1]
    else
        return zero(T)
    end
end

function Base.setindex!(A::AntisymBlockdiagonal, x, i::Integer, j::Integer)
    @boundscheck checkbounds(A, i, j)
    throw(ArgumentError("cannot set entry ($i, $j)"))
    return x
end

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

@auto_hash_equals struct AntisymTridiagonalSchur{T,S<:AbstractMatrix{T}} <: Factorization{T}
    W::S
    values::Vector{T}
    Vt::S
    AntisymTridiagonalSchur{T,S}(W::AbstractMatrix{T}, values::Vector{T}, Vt::AbstractMatrix{T}) where {T,S} = new(W, values, Vt)
end
AntisymTridiagonalSchur(W::AbstractMatrix{T}, values::Vector{T}, Vt::AbstractMatrix{T}) where T = AntisymTridiagonalSchur{T, typeof(W)}(W, values, Vt)

function Base.schurfact(A::AntisymTridiagonal{T}) where T<:Real
    if iseven(size(A, 1))
        C = Bidiagonal(A.ev[1:2:end], -A.ev[2:2:end], false)
    else
        C = Bidiagonal(vcat(A.ev[1:2:end], T[zero(T)]), -A.ev[2:2:end], false)
    end
    U, s, Vt = svd(C)
    if isodd(size(A, 1))
        @assert s[end]==0
        s = s[1:end-1]
        Vt = Vt[1:end-1, 1:end-1]
    end
    return AntisymTridiagonalSchur(U, s, Vt)
end

function unpackZ(F::AntisymTridiagonalSchur{T}) where T
    # n = length(F.values)
    # N = 2n
    W = F.W
    Vt = F.Vt
    m = size(W, 1)
    n = size(Vt, 1)
    N = m+n
    if isodd(N)
        @assert m==n+1
    else
        @assert m==n
    end
    Z = zeros(T, N, N)
    for i in 1:m, j in 1:m
        Z[2i-1, 2j-1] = W[i, j]
    end
    for i in 1:n, j in 1:n
        Z[2i, 2j] = Vt[i, j]
    end
    return Z
end

unpack(F::AntisymTridiagonalSchur) = F.values, unpackZ(F)

function Base.schur(A::AntisymTridiagonal{T}) where T<:Real
    F = schurfact(A)
    values, vecs = unpack(F)
    A0 = AntisymBlockdiagonal(values, size(A, 1)%2)
    evals = Complex{T}[sgn*im*v for v in values for sgn in [1, -1]]
    return A0, vecs, evals
end

##########################################################################################################################
# Next: add generic antisymmetric matrix schur factorization using first a reduction to tridiagonal

struct Antisymmetric{T, S<:AbstractMatrix} <: AbstractMatrix{T}
    data::S
    uplo::Char
end

Antisymmetric(A::AbstractMatrix, uplo::Char='U') = (checksquare(A); Antisymmetric{eltype(A),typeof(A)}(A, uplo))
Antisymmetric(A::AbstractMatrix, uplo::Symbol) = Antisymmetric(A, char_uplo(uplo))

Base.size(A::Antisymmetric, d) = size(A.data, d)
Base.size(A::Antisymmetric) = size(A.data)

@inline function Base.getindex(A::Antisymmetric{T}, i::Integer, j::Integer) where T
    @boundscheck checkbounds(A, i, j)
    @inbounds r = (i==j) ? zero(T) : (A.uplo == 'U') == (i < j) ? A.data[i, j] : -A.data[j, i]
    r
end

function Base.setindex!(A::Antisymmetric, v, i::Integer, j::Integer)
    throw(ArgumentError("Cannot set an index in a antisymmetric matrix"))
end

function Base.replace_in_print_matrix(A::Antisymmetric, i::Integer, j::Integer, s::AbstractString)
    i!=j ? s : Base.replace_with_centered_mark(s)
end

Base.similar(A::Antisymmetric, ::Type{T}) where {T} = Antisymmetric(similar(A.data, T))

Base.convert(::Type{Matrix}, A::Antisymmetric) = copyatri!(convert(Matrix, copy(A.data)), A.uplo)

function copyatri!(A::AbstractMatrix{T}, uplo::Char) where T
    n = checksquare(A)
    for i =1:n
        A[i, i] = zero(T)
    end
    if uplo == 'U'
        for i = 1:(n-1), j = (i+1):n
            A[j,i] = -A[i,j]
        end
    elseif uplo == 'L'
        for i = 1:(n-1), j = (i+1):n
            A[i,j] = -A[j,i]
        end
    else
        throw(ArgumentError("uplo argument must be 'U' (upper) or 'L' (lower), got $uplo"))
    end
    A
end

Base.convert(::Type{Array}, A::Antisymmetric) = convert(Matrix, A)
Base.full(A::Antisymmetric) = convert(Array, A)
Base.parent(A::Antisymmetric) = A.data
Base.copy(A::Antisymmetric{T,S}) where {T,S} = (B = copy(A.data); Antisymmetric{T,typeof(B)}(B,A.uplo))

function Base.copy!(dest::Antisymmetric, src::Antisymmetric)
    if src.uplo == dest.uplo
        copy!(dest.data, src.data)
    else
        transpose!(dest.data, -1*src.data)
    end
    return dest
end

isantisymmetric(A::Antisymmetric) = true
Base.transpose(A::Antisymmetric) = -A
Base.ctranspose(A::Antisymmetric{<:Real}) = -A
Base.trace(A::Antisymmetric) = zero(eltype(A))

Base.:-(A::Antisymmetric{Tv,S}) where {Tv,S<:AbstractMatrix} = Antisymmetric{Tv,S}(-A.data, A.uplo)
Base.:*(A::Antisymmetric, B::StridedMatrix) = full(A)*B
Base.:*(A::StridedMatrix, B::Antisymmetric) = A*full(B)

Base.:*(A::Antisymmetric, x::Number) = Antisymmetric(A.data*x, A.uplo)
Base.:/(A::Antisymmetric, x::Number) = Antisymmetric(A.data/x, A.uplo)

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

@auto_hash_equals struct AntisymmetricSchur{T,S<:AbstractMatrix{T}} <: Factorization{T}
    Q::HessenbergQ{T, S}
    F::AntisymTridiagonalSchur{T, S}
    AntisymmetricSchur{T,S}(Q::HessenbergQ{T, S}, F::AntisymTridiagonalSchur{T, S}) where {T,S} = new{T, S}(Q, F)
end
AntisymmetricSchur(Q::HessenbergQ{T, S}, F::AntisymTridiagonalSchur{T, S}) where {T, S} = AntisymmetricSchur{T, S}(Q, F)

function Base.schurfact(A::Antisymmetric{T}) where T<:Real
    F = hessfact!(copy(Array(A))) # converts A into a tridiagonal form using a unitary rotation F[:Q]
    Atri = AntisymTridiagonal(diag(F.factors, 1))
    G = schurfact(Atri) # diagonalizes the antisymmetric tridiagonal matrix
    return AntisymmetricSchur(F[:Q], G)
end

function Base.schur(A::Antisymmetric{T}) where T<:Real
    F = schurfact(A)
    Z = unpackZ(F.F)

    A0 = AntisymBlockdiagonal(F.F.values, size(A, 1)%2)

    evals = Complex{T}[sgn*im*v for v in F.F.values for sgn in [1, -1]]
    if isodd(size(A, 1))
        evals = vcat(evals, zero(T))
    end
    Q = F.Q
    return  A0, Q*Z, evals
end

unpack(F::AntisymmetricSchur) = (F.F.values, F.Q*unpackZ(F.F))

function pfaffian(A::Antisymmetric{T}) where T<:Real
    n = size(A, 1)
    if n == 0
        return one(T)
    elseif isodd(n)
        return zero(T)
    end
    F=hessfact!(copy(Array(A)))
    Atri = AntisymTridiagonal(diag(F.factors, 1))
    return det(F[:Q])*pfaffian(Atri)
end

function pfaffiansign(A::Antisymmetric{T}) where T<:Real
    n = size(A, 1)
    if n == 0
        return one(T)
    elseif isodd(n)
        return zero(T)
    end
    F=hessfact!(Array(A))
    Atri = AntisymTridiagonal(diag(F.factors, 1))
    return sign(det(F[:Q]))*pfaffiansign(Atri)
end

end # module


text:
using Distributions, Statistics, LinearAlgebra
using DataStructures: CircularBuffer

struct AFS <: AbstractOptimizer
    population_size::Int
    δ_init::Real
    λ::Real
    θ::Real
    m::Int
    η::Real
    step::Real
    mutation::Function

    AFS(; populationSize::Int = 100,
        δ_init::Real = 1.0,
        λ::Real = 1e-3,
        θ::Real = 1.0,
        m::Int = 5,
        η::Real = 1e-3,
        step = 1.0,
        mutation::Function = uniform(1.0),
        ) =
        new(populationSize, δ_init, λ, θ, m, η, step, mutation)
end

population_size(method::AFS) = method.population_size
default_options(method::AFS) = (iterations = 1000, abstol = 1e-32)
summary(m::AFS) = "($(m.δ_init)/$(m.λ)/$(m.θ)/$(m.η)/$(m.m)/$(m.step))"
show(io::IO, m::AFS) = print(io, summary(m))

mutable struct AFSState{T, IT} <: AbstractOptimizerState
    n::Int
    δ::Real
    board::Vector
    fitness::Vector{T}
    fittest::IT
    best_fitnesses::CircularBuffer{T}
    i::Int
end

value(s::AFSState) = minimum(s.fitness)
minimizer(s::AFSState) = s.fittest

function get_norms(population)
    return norm.(population)
end

function initial_state(method::AFS, options, objfun, population)
    T = typeof(value(objfun))
    n = length(first(population))
    fitness = fill(maxintfloat(T), method.population_size)
    return AFSState(n, copy(method.δ_init), deepcopy(population), fitness, population[1], CircularBuffer{T}(method.m), 0)
end

function get_centroid(School::Vector{Vector{T}}) where T<:AbstractFloat
    l = length(School)
    sum(School) ./ l
end

function scan_visual_field(i::Int, field::T, School::Vector{Vector{IT}}) where {T<:AbstractFloat, IT<:AbstractFloat}
    V = @view School[setdiff(1:length(School), i)]
    findall([norm(School[i] - V[j]) for j = 1:length(V)] .> field)
end

function random_move!(S::Vector{Vector{T}}, board::Vector{Vector{T}}, n::Int, field::IT, fish::Int, f::Function) where {T<:AbstractFloat, IT<:AbstractFloat}
    board[fish] = S[fish] + f(zeros(n)) .* field
end

function search!(S::Vector{Vector{T}}, board::Vector{Vector{T}}, n::Int, fish::Int, V::Vector, step::IT, f::Function) where {T<:AbstractFloat, IT<:AbstractFloat}
    j = rand(1:length(V))
    d = S[V[j]] - S[fish]
    board[fish] = S[fish] + step * f(zeros(n)) .* d ./ abs.(d)
end

function swarm!(S::Vector{Vector{T}}, board::Vector{Vector{T}}, n::Int, fish::Int, centroid::Vector{T}, step::IT, f::Function) where {T<:AbstractFloat, IT<:AbstractFloat}
    d = centroid - S[fish]
    board[fish] = S[fish] + step * f(zeros(n)) .* d ./ abs.(d)
end

function chase!(S::Vector{Vector{T}}, board::Vector{Vector{T}}, n::Int, fish::Int, neighbor::Int, step::IT, f::Function) where {T<:AbstractFloat, IT<:AbstractFloat}
    d = S[neighbor] - S[fish]
    board[fish] = S[fish] + step * f(zeros(n)) .* d ./ abs.(d)
end

function leap!(S::Vector{Vector{T}}, fish::Int, n::Int, max_norm::IT, f::Function) where {T<:AbstractFloat, IT<:AbstractFloat}
    S[fish] = S[fish] + f(zeros(n)) .* (max_norm * 2 .- S[fish])
end

function eval_fitness(constraints, objfun, x)
    value(objfun, x) + penalty(constraints, x)
end

function update_state!(objfun, constraints, state, population::AbstractVector{IT}, method::AFS, itr) where {IT}
    state.i += 1
    max_norm = first(findmax(get_norms(population)))
    v = state.δ * max_norm
    for fish = 1:method.population_size
        V = scan_visual_field(fish, v, population)
        if isempty(V)
            random_move!(population, state.board, state.n, v, fish, method.mutation)
        elseif length(V) / method.population_size > method.θ
            search!(population, state.board, state.n, fish, V, method.step, method.mutation)
        else
            c = get_centroid(population[V])
            state.fitness[fish] = eval_fitness(constraints, objfun, population[fish])
            if eval_fitness(constraints, objfun, c) < state.fitness[fish]
                swarm!(population, state.board, state.n, fish, c, method.step, method.mutation)
            else
                search!(population, state.board, state.n, fish, V, method.step, method.mutation)
            end
            best_neighbor = findmin(state.fitness[V])
            if first(best_neighbor) < state.fitness[fish]
                chase!(population, state.board, state.n, fish, V[last(best_neighbor)], method.step, method.mutation)
            else
                search!(population, state.board, state.n, fish, V, method.step, method.mutation)
            end
        end
    end

    for fish = 1:method.population_size
        o = apply!(constraints, state.board[fish])
        y_fit = value(objfun, o) + penalty(constraints, o)
        if y_fit < state.fitness[fish]
            population[fish] = o
            state.fitness[fish] = y_fit
        end
    end

    minfit = findmin(state.fitness)
    state.fittest = population[last(minfit)]
    pushfirst!(state.best_fitnesses, first(minfit))

    if state.i >= method.m
        stagnation = last(state.best_fitnesses) - first(state.best_fitnesses)
        if stagnation < method.η
            j = rand(1:method.population_size)
            leap!(population, j, state.n, max_norm, method.mutation)
            apply!(constraints, population[j])
        end
        if state.best_fitnesses[end] <= state.best_fitnesses[end-1]
            state.δ *= exp(-method.λ)
        end
    end

    return false

end


text:
using Base: Broadcast.DefaultArrayStyle, AbstractUnitRange

_broadcasted(style, f, x, r::AbstractRange) = invoke(Base.Broadcast.broadcasted, Tuple{typeof(style),typeof(f),Any,typeof(r)}, style, f, x, r)
_broadcasted(style, f, r::AbstractRange, x) = invoke(Base.Broadcast.broadcasted, Tuple{typeof(style),typeof(f),typeof(r),Any}, style, f, r, x)

for f in (+, -, *, /, \), R in (AbstractRange, AbstractUnitRange, StepRangeLen)
    @eval begin
        Base.Broadcast.broadcasted(style::DefaultArrayStyle{1}, ::typeof($f), x::Measurement, r::$R) = _broadcasted(style, $f, x, r)
        Base.Broadcast.broadcasted(style::DefaultArrayStyle{1}, ::typeof($f), r::$R, x::Measurement) = _broadcasted(style, $f, r, x)
    end
end


text:
using SparseTransforms
using Test

function test_spright()
    println("Testing input signal constructor...")
    n, b = 4, 2
    signal = TestSignal(4, [4, 6, 10, 15])
    signal = TestSignal(4, [4, 6, 10, 15], Float64.([2, 4, 1, 1]))
    @test signal.n == n

    println("Testing query methods...")
    b = get_b(signal)
    Ms = get_Ms(n, b)
    D = get_D(n)
    whts_and_inds = map(M -> compute_delayed_subtransform(signal, M, D, fwht), Ms)

    println("Testing reconstruction...")
    bin_cardinalities = map(M -> bin_cardinality(signal, M, D), Ms)

    println("Testing full transform...")
    # methods_list = [:simple, :random, :mle, :none]
    methods_list = [:simple, :nso, :nso, :none]
    spright_wht, utilization = spright(signal, methods_list; verbose=true, report=true)
    println("Used $(utilization / 2^n) of all time samples")

    println("Checking transform result...")
    true_wht = Dict(l => s for (l, s) in zip(signal.locs, signal.strengths))
    @test length(signal.locs) == length(spright_wht)
    for (i, x) in spright_wht
        @test true_wht[i] ≈ x
    end
    println("Transform is correct!")
end


text:
function Escape.loo(result::HLAModelResult)
    d = stan_input(result)
    N = d["N"]
    pw = Vector{Loo.PointwiseLoo}(undef, N)
    p = extract(stanfit(result))

    progress = Progress(N)
    @threads for i in 1:N
        ll = pointwise_loglikelihoods(result.sf, p, i)
        pw[i] = Loo.pointwise_loo(ll)
        next!(progress)
    end

    draws = result.sf.iter * result.sf.chains

    return Loo.LooResult(pw, (draws, length(pw)))
end

function pointwise_loglikelihoods(sf::StanInterface.Stanfit, p::Dict, i::Int)
    y = sf.data["y"][i]
    iter = sf.iter
    chains = sf.chains

    ll = Vector{Vector{Float64}}(undef, chains)
    for idx in eachindex(sf.result)
        theta = theta_i(sf, p, i)
        ll[idx] = y == 1 ? log.(theta) : log.(1 .- theta)
    end

    return ll
end

function theta_i(sf::StanInterface.Stanfit, p::Dict, i::Int)
    stan_input = sf.data

    rs = stan_input["rs"]
    R = stan_input["R"]
    D = stan_input["D"]
    X = stan_input["X"]
    idx = stan_input["idx"]
    phy = stan_input["phy"]

    r = rs[i]

    nu = zeros(length(p["beta_hla.1.1"]))

    for d in 1:D
        beta_hla = p["beta_hla.$r.$d"]

        nu += X[idx[i], d] .* beta_hla
    end

    b_phy = "b_phy.$r" in keys(p) ? p["b_phy.$r"] : zeros(length(p["beta_hla.1.1"]))
    theta = @. logistic(p["b0_hla.$r"] + (b_phy * logit(phy[r, idx[i]])) + nu)

    return theta
end


text:
module Orbits

using StaticArrays

export SimpleOrbit, KeplerianOrbit

abstract type AbstractOrbit end

Base.broadcastable(o::AbstractOrbit) = Ref(o)

"""
    relative_position(::AbstractOrbit, t)
"""
relative_position(::AbstractOrbit, t)

include("simple.jl")
include("kepler.jl")

end # module