text:
using Test
import ArchGDAL as AG

@testset "test_feature.jl" begin
    AG.read("data/point.geojson") do dataset
        layer = AG.getlayer(dataset, 0)
        AG.getfeature(layer, 0) do f1
            AG.getfeature(layer, 2) do f2
                @test sprint(print, f1) == """
                Feature
                  (index 0) geom => POINT
                  (index 0) FID => 2.0
                  (index 1) pointname => point-a
                """
                AG.getgeom(f1) do g1
                    @test sprint(print, g1) == "Geometry: POINT (100 0)"
                end
                fid1 = AG.getfid(f1)
                @test fid1 == 0

                AG.getgeom(f2, "fake name") do g
                    @test sprint(print, g) == "NULL Geometry"
                end

                @test sprint(print, f2) == """
                Feature
                  (index 0) geom => POINT
                  (index 0) FID => 0.0
                  (index 1) pointname => a
                """
                AG.getgeom(f2, "") do g2
                    @test sprint(print, g2) == "Geometry: POINT (100 0)"
                end
                fid2 = AG.getfid(f2)
                @test fid2 == 2

                AG.clone(f1) do f3
                    @test AG.equals(AG.getgeom(f1), AG.getgeom(f3)) == true
                end
                AG.setfid!(f1, fid2)
                AG.setfid!(f2, fid1)
                @test AG.getfid(f1) == 2
                @test AG.getfid(f2) == 0

                @test AG.findgeomindex(f1, "geom") == -1
                @test AG.findgeomindex(f1, "") == 0
                @test AG.findgeomindex(f2, "geom") == -1
                @test AG.findgeomindex(f2, "") == 0
                @test AG.gettype(AG.getgeomdefn(f1, 0)) == AG.wkbPoint
                @test AG.gettype(AG.getgeomdefn(f2, 0)) == AG.wkbPoint
            end
        end

        AG.getfeature(layer, 0) do f
            @test AG.toWKT(AG.getgeom(f, 0)) == "POINT (100 0)"
            AG.setgeom!(f, 0, AG.createpoint(0, 100))
            @test AG.toWKT(AG.getgeom(f, 0)) == "POINT (0 100)"
            AG.createpolygon([(0.0, 100.0), (100.0, 0.0)]) do poly
                return AG.setgeom!(f, 0, poly)
            end
            @test AG.toWKT(AG.getgeom(f, 0)) == "POLYGON ((0 100,100 0))"

            AG.setstylestring!(f, "@Name")
            @test AG.getstylestring(f) == "@Name"
            AG.setstylestring!(f, "NewName")
            @test AG.getstylestring(f) == "NewName"

            AG.createstyletable() do st
                AG.addstyle!(st, "name", "style")
                AG.setstyletable!(f, st)
                @test AG.findstylestring(AG.getstyletable(f), "name") == "style"
            end

            AG.setnativedata!(f, "nativedata1")
            @test AG.getnativedata(f) == "nativedata1"
            AG.setnativedata!(f, "nativedata2")
            @test AG.getnativedata(f) == "nativedata2"

            AG.setmediatype!(f, "mediatype1")
            @test AG.getmediatype(f) == "mediatype1"
            AG.setmediatype!(f, "mediatype2")
            @test AG.getmediatype(f) == "mediatype2"

            @test AG.validate(f, AG.F_VAL_NULL, false) == true
            @test AG.validate(f, AG.F_VAL_GEOM_TYPE, false) == false
            @test AG.validate(f, AG.F_VAL_WIDTH, false) == true
            @test AG.validate(f, AG.F_VAL_ALLOW_NULL_WHEN_DEFAULT, false) ==
                  true
            @test AG.validate(f, AG.F_VAL_ALLOW_DIFFERENT_GEOM_DIM, false) ==
                  true

            @testset "Missing and Null Semantics" begin
                @test isnothing(AG.getdefault(f, 1))
                AG.setdefault!(AG.getfielddefn(f, 1), "default value")
                @test AG.getdefault(f, 1) == "default value"

                @test AG.isfieldsetandnotnull(f, 1)
                @test AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1)
                @test AG.getfield(f, 1) == "point-a"

                AG.unsetfield!(f, 1)
                @test !AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1) # carried over from earlier
                @test isnothing(AG.getfield(f, 1))

                # unset & notnull: missing
                AG.fillunsetwithdefault!(f)
                # nothing has changed
                @test isnothing(AG.getfield(f, 1))
                # because it is a nullable field
                @test AG.isnullable(AG.getfielddefn(f, 1))
                # even though it is not a null value
                @test !AG.isfieldnull(f, 1)
                # the field is still not set
                @test !AG.isfieldset(f, 1)

                # set & notnull: value
                AG.fillunsetwithdefault!(f, notnull = false)
                # now the field is set to the default
                @test AG.getfield(f, 1) == AG.getdefault(f, 1)
                @test !AG.isfieldnull(f, 1) # still as expected
                @test AG.isfieldset(f, 1) # the field is now set

                # set the field to be notnullable
                AG.setnullable!(AG.getfielddefn(f, 1), false)
                # now if we unset the field
                AG.unsetfield!(f, 1)
                @test !AG.isfieldnull(f, 1)
                @test !AG.isfieldset(f, 1)
                @test isnothing(AG.getfield(f, 1))
                # and we fill unset with default again
                AG.fillunsetwithdefault!(f)
                # the field is set to the default
                @test AG.getfield(f, 1) == AG.getdefault(f, 1)

                # set & null: missing
                @test !AG.isfieldnull(f, 1)
                @test AG.isfieldset(f, 1)
                AG.setfieldnull!(f, 1)
                @test AG.isfieldnull(f, 1)
                @test AG.isfieldset(f, 1)
                @test ismissing(AG.getfield(f, 1))

                # unset & null: N/A (but nothing otherwise)
                AG.unsetfield!(f, 1)
                # Observe that OGRUnset and OGRNull are mutually exclusive
                @test !AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1) # notice the field is notnull

                # setting the field for a notnullable column
                AG.setnullable!(AG.getfielddefn(f, 1), false)
                AG.setfield!(f, 1, "value")
                @test AG.getfield(f, 1) == "value"
                @test AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1)
                AG.setfield!(f, 1, missing)
                @test AG.getfield(f, 1) == AG.getdefault(f, 1)
                @test AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1)
                AG.setfield!(f, 1, nothing)
                @test isnothing(AG.getfield(f, 1))
                @test !AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1)

                # setting the field for a nullable column
                AG.setnullable!(AG.getfielddefn(f, 1), true)
                AG.setfield!(f, 1, "value")
                @test AG.getfield(f, 1) == "value"
                @test AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1)
                AG.setfield!(f, 1, missing)
                @test ismissing(AG.getfield(f, 1))
                @test AG.isfieldset(f, 1)
                @test AG.isfieldnull(f, 1) # different from that of notnullable
                AG.setfield!(f, 1, nothing)
                @test isnothing(AG.getfield(f, 1))
                @test !AG.isfieldset(f, 1)
                @test !AG.isfieldnull(f, 1)
            end
        end
    end

    @testset "In-Memory Driver" begin
        AG.create(AG.getdriver("MEMORY")) do output
            layer = AG.createlayer(dataset = output, geom = AG.wkbPolygon)
            AG.createfielddefn("int64field", AG.OFTInteger64) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("float64field", AG.OFTReal) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("intlistfield", AG.OFTIntegerList) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn(
                "int64listfield",
                AG.OFTInteger64List,
            ) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("float64listfield", AG.OFTRealList) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("stringlistfield", AG.OFTStringList) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("binaryfield", AG.OFTBinary) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("datetimefield", AG.OFTDateTime) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("booleanfield", AG.OFTInteger) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("int16field", AG.OFTInteger) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("int32field", AG.OFTInteger) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfielddefn("float32field", AG.OFTReal) do fielddefn
                return AG.addfielddefn!(layer, fielddefn)
            end
            AG.createfeature(layer) do feature
                AG.setfield!(feature, 0, Int64(1))
                AG.setfield!(feature, 1, Float64(1.0))
                AG.setfield!(feature, 2, Int32[1, 2])
                AG.setfield!(feature, 3, Int64[1, 2])
                AG.setfield!(feature, 4, Float64[1.0, 2.0])
                AG.setfield!(feature, 5, ["1", "2.0"])
                AG.setfield!(feature, 6, UInt8[1, 2, 3, 4])
                AG.setfield!(feature, 7, Dates.DateTime(2016, 9, 25, 21, 17, 0))
                AG.setfield!(feature, 8, true)
                AG.setfield!(feature, 9, Int16(1))
                AG.setfield!(feature, 10, Int32(1))
                AG.setfield!(feature, 11, Float32(1.0))
                for i in 1:AG.nfield(feature)
                    @test !AG.isfieldnull(feature, i - 1)
                    @test AG.isfieldsetandnotnull(feature, i - 1)
                end
                @test sprint(print, AG.getgeom(feature)) == "NULL Geometry"
                AG.getgeom(feature) do geom
                    @test sprint(print, geom) == "NULL Geometry"
                end
                @test sprint(print, AG.getgeom(feature, 0)) == "NULL Geometry"
                AG.getgeom(feature, 0) do geom
                    @test sprint(print, geom) == "NULL Geometry"
                end

                AG.addfeature(layer) do newfeature
                    AG.setfrom!(newfeature, feature)
                    @test AG.getfield(newfeature, 0) == 1
                    for i in 1:AG.nfield(newfeature)
                        @test !AG.isfieldnull(newfeature, i - 1)
                        @test AG.isfieldsetandnotnull(newfeature, i - 1)
                    end
                    @test AG.getfield(newfeature, 1) ≈ 1.0
                    @test AG.getfield(newfeature, 2) == Int32[1, 2]
                    @test AG.getfield(newfeature, 3) == Int64[1, 2]
                    @test AG.getfield(newfeature, 4) ≈ Float64[1.0, 2.0]
                    @test AG.getfield(newfeature, 5) == String["1", "2.0"]
                    @test AG.getfield(newfeature, 6) == UInt8[1, 2, 3, 4]
                    @test AG.getfield(newfeature, 7) ==
                          Dates.DateTime(2016, 9, 25, 21, 17, 0)
                    @test AG.getfield(newfeature, 8) == true
                    @test AG.getfield(newfeature, 9) == 1
                    @test AG.getfield(newfeature, 10) == 1
                    @test AG.getfield(newfeature, 11) == 1.0

                    AG.createfeature(layer) do lastfeature
                        AG.setfrom!(lastfeature, feature)
                        AG.setfield!(lastfeature, 0, 45)
                        AG.setfield!(lastfeature, 1, 18.2)
                        AG.setfield!(lastfeature, 5, ["foo", "bar"])
                        @test AG.getfield(lastfeature, 0) == 45
                        @test AG.getfield(lastfeature, 1) ≈ 18.2
                        @test AG.getfield(lastfeature, 2) == Int32[1, 2]
                        @test AG.getfield(lastfeature, 3) == Int64[1, 2]
                        @test AG.getfield(lastfeature, 4) ≈ Float64[1.0, 2.0]
                        @test AG.getfield(lastfeature, 5) ==
                              String["foo", "bar"]
                        @test AG.getfield(lastfeature, 6) == UInt8[1, 2, 3, 4]
                        @test AG.getfield(lastfeature, 7) ==
                              Dates.DateTime(2016, 9, 25, 21, 17, 0)

                        @test AG.getfield(newfeature, 0) == 1
                        @test AG.getfield(newfeature, 1) ≈ 1.0
                        @test AG.getfield(newfeature, 5) == String["1", "2.0"]
                        AG.setfrom!(
                            newfeature,
                            lastfeature,
                            collect(Cint, 0:AG.nfield(newfeature)),
                        )
                        @test AG.getfield(newfeature, 0) == 45
                        @test AG.getfield(newfeature, 1) ≈ 18.2
                        @test AG.getfield(newfeature, 5) == String["foo", "bar"]

                        @test AG.isfieldsetandnotnull(newfeature, 5)
                        AG.setfieldnull!(newfeature, 5)
                        @test !AG.isfieldsetandnotnull(newfeature, 5)
                        @test AG.isfieldset(newfeature, 5)
                        @test AG.isfieldnull(newfeature, 5)
                        @test ismissing(AG.getfield(newfeature, 5))
                    end
                    @test AG.nfeature(layer) == 1
                end
                @test AG.nfeature(layer) == 2
            end
            @test AG.nfeature(layer) == 3
        end
    end
end


text:
"""
    execute(command[; client::Client=get_global_client()])

Sends a RESP compliant command to the Redis host and returns the result. The command is either an
array of command keywords, or a single command string. Defaults to using the globally set Client.

# Examples
```julia-repl
julia> execute(["SET", "key", "value"])
"OK"
julia> execute("GET key")
"value"
```
"""
function execute(command::AbstractArray, client::Client=get_global_client())
    if client.is_subscribed
        throw(RedisError("SUBERROR", "Cannot execute commands while a subscription is open in the same Client instance"))
    end

    @lock client.lock begin
        flush!(client)
        write(client.socket, resp(command))
        msg = recv(client.socket)

        if msg isa Exception
            throw(msg)
        end

        return msg
    end
end
function execute(command::AbstractString, client::Client=get_global_client())
    execute(split_on_whitespace(command), client)
end

"""
    execute_without_recv(command[; client::Client=get_global_client()])

Sends a RESP compliant command to the Redis host without reading the returned result.
"""
function execute_without_recv(command::AbstractArray, client::Client=get_global_client())
    @lock client.lock begin
        flush!(client)
        write(client.socket, resp(command))
        return
    end
end
function execute_without_recv(command::AbstractString, client::Client=get_global_client())
    execute_without_recv(split_on_whitespace(command), client)
end

"""
    execute(command, pipe::Pipeline)

Add a RESP compliant command to a pipeline client, optionally filter out MULTI transaction responses
before the EXEC call, e.g. "QUEUED".

# Examples
```julia-repl
julia> pipe = Pipeline();

julia> execute(["SET", "key", "value"]; client=pipe);

julia> execute(["GET", "key"]; client=pipe);

julia> execute(pipe)
2-element Array{String,1}:
 "OK"
 "value"
```
"""
function execute(command::AbstractArray, pipe::Pipeline)
    add!(pipe, command)
    return
end
function execute(command::AbstractString, pipe::Pipeline)
    execute(split_on_whitespace(command), pipe)
end

"""
    execute(pipe::Pipeline[; filter_multi_exec=true])

Execute commands batched in a pipeline client, optionally filter out MULTI transaction responses
before the EXEC call, e.g. "QUEUED".

# Examples
```julia-repl
julia> pipe = Pipeline();

julia> set("key", "value"; client=pipe);

julia> get("key"; client=pipe);

julia> multi(; client=pipe);

julia> get("key"; client=pipe);

julia> get("key"; client=pipe);

julia> exec(; client=pipe);

julia> execute(pipe)
2-element Array{String,1}:
 "OK"
 "value"
 ["value", "value"]  # Only the response from final exec() call is returned
```
"""
function execute(pipe::Pipeline)
    if pipe.client.is_subscribed
        throw(RedisError("SUBERROR", "Cannot execute Pipeline while a subscription is open in the same Client instance"))
    end

    @lock pipe.client.lock begin
        try
            flush!(pipe.client)
            write(pipe.client.socket, join(pipe.resp))
            messages = [recv(pipe.client.socket) for _ in 1:length(pipe.resp)]

            if pipe.filter_multi_exec
                return messages[pipe.multi_exec_bitmask]
            end

            return messages
        finally
            flush!(pipe)
        end
    end
end

ext: jl


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

precompile(Tuple{Type{Base.Generator{I, F} where F where I}, typeof(Pkg.TOML.table2dict), Array{Any, 1}})
precompile(Tuple{Type{Base.IOContext{IO_t} where IO_t<:IO}, Base.GenericIOBuffer{Array{UInt8, 1}}, Base.TTY})
precompile(Tuple{Type{Base.SHA1}, Base.SubString{String}})
precompile(Tuple{Type{Base.SHA1}, String})
precompile(Tuple{Type{NamedTuple{(:env,), Tuple{Pkg.Types.EnvCache}}}, Tuple{Pkg.Types.EnvCache}})
precompile(Tuple{Type{Pkg.Display.DiffEntry}, Base.UUID, String, Nothing, Pkg.Display.VerInfo})
precompile(Tuple{Type{Pkg.Display.VerInfo}, Base.SHA1, Nothing, Base.VersionNumber, Bool, Nothing})
precompile(Tuple{Type{Pkg.Display.VerInfo}, Base.SHA1, Nothing, Base.VersionNumber, Bool, Pkg.Types.GitRepo})
precompile(Tuple{Type{Pkg.Types.GitRepo}, String, String})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), NamedTuple{(:env,), Tuple{Pkg.Types.EnvCache}}, Type{Pkg.Types.Context}})
precompile(Tuple{typeof(Base.:(==)), Base.SHA1, Nothing})
precompile(Tuple{typeof(Base.:(==)), Base.UUID, Base.UUID})
precompile(Tuple{typeof(Base.:(==)), Nothing, Base.UUID})
precompile(Tuple{typeof(Base.:(==)), Nothing, Nothing})
precompile(Tuple{typeof(Base.:(==)), Pkg.REPLMode.CommandKind, Pkg.REPLMode.CommandKind})
precompile(Tuple{typeof(Base.:(==)), String, Nothing})
precompile(Tuple{typeof(Base.append!), Array{String, 1}, Array{Base.SubString{String}, 1}})
precompile(Tuple{typeof(Base.close), LibGit2.GitRepo})
precompile(Tuple{typeof(Base.collect), Base.Generator{Array{Any, 1}, typeof(Pkg.TOML.table2dict)}})
precompile(Tuple{typeof(Base.deepcopy_internal), Array{Base.Dict{String, Any}, 1}, Base.IdDict{Any, Any}})
precompile(Tuple{typeof(Base.deepcopy_internal), Base.Dict{Any, Any}, Base.IdDict{Any, Any}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol, Symbol, Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol, Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{typeof(Base.diff_names), Tuple{Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{typeof(Base.eachmatch), Base.Regex, String})
precompile(Tuple{typeof(Base.empty), Array{Union{Nothing, Base.SubString{String}}, 1}, Type{Base.SubString{String}}})
precompile(Tuple{typeof(Base.getindex), Array{Base.Dict{String, Any}, 1}, Int64})
precompile(Tuple{typeof(Base.getproperty), Pkg.REPLMode.Command, Symbol})
precompile(Tuple{typeof(Base.getproperty), Pkg.TOML.Table, Symbol})
precompile(Tuple{typeof(Base.haskey), Base.Dict{Base.UUID, Base.Dict{K, V} where V where K}, Base.UUID})
precompile(Tuple{typeof(Base.haskey), Pkg.TOML.Table, String})
precompile(Tuple{typeof(Base.isempty), Array{Base.Dict{String, Any}, 1}})
precompile(Tuple{typeof(Base.iterate), Array{Base.Dict{String, Any}, 1}, Int64})
precompile(Tuple{typeof(Base.iterate), Array{Base.Dict{String, Any}, 1}})
precompile(Tuple{typeof(Base.length), Array{Base.Dict{String, Any}, 1}})
precompile(Tuple{typeof(Base.merge), Base.Dict{String, Any}, Base.Dict{String, String}})
precompile(Tuple{typeof(Base.merge_names), Tuple{Symbol}, Tuple{Symbol}})
precompile(Tuple{typeof(Base.merge_types), Tuple{Symbol, Symbol}, Type{NamedTuple{(:callbacks,), Tuple{Base.Dict{Symbol, Tuple{Ptr{Nothing}, Any}}}}}, Type{NamedTuple{(:credentials,), Tuple{LibGit2.CachedCredentials}}}})
precompile(Tuple{typeof(Base.push!), Array{Base.SubString{String}, 1}, Base.SubString{String}})
precompile(Tuple{typeof(Base.read), Base.IOStream, Type{String}})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Base.UUID, Base.Dict{K, V} where V where K}, Base.Dict{String, Any}, Base.UUID})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Array{Any, 1}, String})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Array{Base.Dict{String, Any}, 1}, String})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Base.Dict{Any, Any}, String})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{String, Any}, Pkg.TOML.Table, String})
precompile(Tuple{typeof(Base.similar), Array{Base.Dict{String, Any}, 1}})
precompile(Tuple{typeof(Base.string), String, Base.UUID})
precompile(Tuple{typeof(Base.string), String, Module, String, String, String, Int64})
precompile(Tuple{typeof(Base.sym_in), Symbol, Tuple{Symbol}})
precompile(Tuple{typeof(Base.write), Base.TTY, Array{UInt8, 1}})
precompile(Tuple{typeof(Core.Compiler.vect), Tuple{Nothing, Int64}})
precompile(Tuple{typeof(Core.Compiler.vect), Tuple{Pkg.REPLMode.OptionKind, Int64}})
precompile(Tuple{typeof(Pkg.Display.filter_manifest), Pkg.Display.InProject{Base.Dict{String, Any}}, Base.Dict{String, Any}})
precompile(Tuple{typeof(Pkg.Display.in_project), Base.Dict{String, Any}})
precompile(Tuple{typeof(Pkg.Display.manifest_diff), Pkg.Types.Context, Base.Dict{String, Any}, Base.Dict{String, Any}})
precompile(Tuple{typeof(Pkg.Display.name_ver_info), Base.Dict{String, Any}})
precompile(Tuple{typeof(Pkg.Display.not_in_project), Base.Dict{String, Any}})
precompile(Tuple{typeof(Pkg.Display.print_diff), Base.TTY, Pkg.Types.Context, Array{Pkg.Display.DiffEntry, 1}, Bool})
precompile(Tuple{typeof(Pkg.GitTools.showprogress), Base.TTY, Pkg.GitTools.MiniProgressBar})
precompile(Tuple{typeof(Pkg.GitTools.transfer_progress), Ptr{LibGit2.TransferProgress}, Base.Dict{Symbol, Any}})
precompile(Tuple{typeof(Pkg.REPLMode.do_status!), Pkg.Types.Context, Array{Union{Pkg.Types.VersionRange, String, Pkg.REPLMode.Command, Pkg.REPLMode.Option, Pkg.REPLMode.Rev}, 1}})
precompile(Tuple{typeof(Pkg.REPLMode.promptf)})
precompile(Tuple{typeof(Pkg.TOML.SOME), Array{String, 1}})
precompile(Tuple{typeof(Pkg.TOML.SOME), Pkg.TOML.Table})
precompile(Tuple{typeof(Pkg.TOML.insertpair), Pkg.TOML.Parser{Base.IOStream}, Pkg.TOML.Table, String, Array{String, 1}, Int64})
precompile(Tuple{typeof(Pkg.TOML.insertpair), Pkg.TOML.Parser{Base.IOStream}, Pkg.TOML.Table, String, String, Int64})
precompile(Tuple{typeof(Pkg.Types.printpkgstyle), Base.TTY, Symbol, String})
precompile(Tuple{typeof(REPL.LineEdit.complete_line), Pkg.REPLMode.PkgCompletionProvider, REPL.LineEdit.PromptState})
precompile(Tuple{typeof(REPL.LineEdit.write_prompt), REPL.Terminals.TerminalBuffer, typeof(Pkg.REPLMode.promptf)})

@assert precompile(Tuple{typeof(Pkg.REPLMode.do_help!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}, REPL.LineEditREPL})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_status!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_add_or_develop!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}, Pkg.REPLMode.CommandKind})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_rm!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_up!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_pin!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_free!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_test!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})
@assert precompile(Tuple{typeof(Pkg.REPLMode.do_build!), Pkg.Types.Context, Vector{Pkg.REPLMode.Token}})

max_stars_repo_name: perrutquist/julia

max_issues_repo_name: perrutquist/julia

max_forks_repo_name: perrutquist/julia


text:
module POMDPFiles

using Reexport
using POMDPs
using POMDPModelTools
using Printf
@reexport using POMDPXFiles # for POMDPAlphas

import POMDPs: action, value

export
	POMDPFile,

	read_alpha

include("read.jl")
include("write.jl")

end # module


text:
export BPoint, cartesian, cross, angle_vector, project, angle_scalar, ×, unit, norm_angles, dim_count

struct BPoint
    ρ::Float64
    angles::Array{Float64, 1}

    function BPoint(length, angles::T...) where T <: Number
        new(Float64(length), Float64[Float64(angle) for angle in angles])
    end

    BPoint(length, angles::Array{Float64, 1}) = new(Float64(length), angles)
    BPoint(n::Number) = new(Float64(n), [0])

    "Convert Cartesian to Spherical"
    function BPoint(coordinates::Vector{Float64})
        ρ::Float64 = 0
        for coord in coordinates
            ρ += coord ^ 2
        end

        partialp = ρ
        ρ = sqrt(ρ)

        addPI = coordinates[end] < 0
        deleteat!(coordinates, length(coordinates))

        for i in 1:length(coordinates)
            value = acos(coordinates[i] / sqrt(partialp))
            partialp -= coordinates[i] ^ 2
            coordinates[i] = value
        end

        if addPI
            coordinates[end] = 2 * pi - coordinates[end]
        end

        BPoint(ρ, coordinates)
    end

    Base.length(p::BPoint) = p.ρ

    Base.:+(p::BPoint, n::Number)::BPoint = BPoint(p.ρ + n, p.angles)
    Base.:-(p::BPoint, n::Number)::BPoint = BPoint(p.ρ - n, p.angles)
    Base.:*(p::BPoint, n::Number)::BPoint = BPoint(p.ρ * n, p.angles)
    Base.:/(p::BPoint, n::Number)::BPoint = BPoint(p.ρ / n, p.angles)
    Base.:^(p::BPoint, pow::Number)::BPoint = BPoint(p.ρ ^ pow, p.angles)

    Base.:+(n::Number, p::BPoint)::BPoint = BPoint(p.ρ + n, p.angles)
    Base.:-(n::Number, p::BPoint)::BPoint = BPoint(n - p.ρ, p.angles)
    Base.:*(n::Number, p::BPoint)::BPoint = BPoint(p.ρ * n, p.angles)
    Base.:/(n::Number, p::BPoint)::BPoint = BPoint(n / p.ρ, p.angles)
    Base.:^(n::Number, pow::BPoint)::BPoint = BPoint(pow ^ p.ρ, p.angles)

    function Base.:+(p1::BPoint, p2::BPoint)::BPoint
        angles = norm_angles(p1, p2)
        return BPoint(p1.ρ + p2.ρ, angles)
    end

    function Base.:*(p1::BPoint, p2::BPoint)::BPoint
        angles = norm_angles(p1, p2)
        return BPoint(p1.ρ * p2.ρ, angles)
    end

    function Base.:/(p1::BPoint, p2::BPoint)::BPoint
        angles = norm_angles(p1, p2)
        return BPoint(p1.ρ / p2.ρ, angles)
    end

    function Base.:-(p1::BPoint, p2::BPoint)::BPoint
        angles = norm_angles(p1, p2)
        return BPoint(p1.ρ - p2.ρ, angles)
    end

end

dim_count(p)::Int64 = length(p.angles)

"Project p1 onto p2"
project(p1::BPoint, p2::BPoint)::BPoint = BPoint(p1.ρ * p2.ρ * cos(angle_scalar(p1, p2)), p2.angles)

"Convert to Cartesian"
function cartesian(p)::Array{Float64, 1}
     coordinates = zeros(Float64, dim_count(p) + 1)
     coordinates[1] = p.ρ * cos(p.angles[1])

     sin_product = p.ρ * sin(p.angles[1])
     for i in 2:(length(coordinates) - 1)
         coordinates[i] = sin_product * cos(p.angles[i])
         sin_product *= sin_product * sin(p.angles[i])
     end

     coordinates[end] = sin_product
     coordinates
end

"Angle Number Between Numbers"
function angle_vector(p1::BPoint, p2::BPoint)::Array{Float64, 1}
    check(p1, p2)
    angles = zeros(Float64, dim_count(p1))
    for i in 1:dim_count(p1)
        angles[i] = p1.angles[i] - p2.angles[i]
    end
    angles
end

angle_scalar(p1::BPoint, p2::BPoint)::Float64 = sum(angle_vector(p1, p2))

"Unit Number"
unit(p::BPoint)::BPoint = BPoint(1, p.angles)
unit(angles::Float64...)::BPoint = BPoint(1, collect(angles))
unit(angles::AbstractArray{Float64})::BPoint = BPoint(1, angles)

function check(p1::BPoint, p2::BPoint)
    if dim_count(p1) != dim_count(p2)
        throw(DimensionMismatch("Dimension Mismatch between $p1 and $p2"))
    end
end

"Normalized Angles"
function norm_angles(p1::BPoint, p2::BPoint)::Array{Float64, 1}
    check(p1, p2)
    angles = zeros(Float64, dim_count(p1))
    for i in 1:dim_count(p1)
        angles[i] = (p1.ρ * p1.angles[i] + p2.ρ * p2.angles[i]) / (p1.ρ + p2.ρ)
    end
    angles
end

"Cross Product"
function ×(p1::BPoint, p2::BPoint)::BPoint
    check(p1, p2)
    angle_add = angle_scalar(p1, p2) > pi ? true : false
    anglev = angle_vector(p1, p2)
    for i in 1:dim_count(p1)
        anglev[i] = angle_add ? anglev[i] - (pi/2) : anglev[i] + (pi/2)
    end

    #TODO Change to be area of parallelpid
    return BPoint(p1.ρ * p2.ρ, anglev)
end

cross(p1::BPoint, p2::BPoint)::BPoint = p1 × p2


text:
using Test, FilePaths, Bonsai
using HTTP: Response
using FilePathsBase: /
using JSON3

@testset "Path" begin
	file = joinpath(Path(@__DIR__), "data/a.json")

	# writing
	res = Response()
	Bonsai.write(res, file)
	k, v = res.headers[1]
	@test k == "content-type"
	@test v == "application/json"
	@test JSON3.read(res.body, Dict) isa Dict
end

max_stars_repo_name: onetonfoot/Tree

max_issues_repo_name: onetonfoot/Tree

max_forks_repo_name: onetonfoot/Tree


text:
function GRASP_single_thread(M1::Number, M2::Number, weights::Array{Float64, 1}, tasks_by_w::Dict{Int, LTask}, prec::Dict{Int, Array}, bj::Array{Int, 1}, CTS::Constants)
    #it_results=Array{Tuple{Number, Number, String, String}, 1}()
    LS = init_ls(CTS)
    QC = init_qc(CTS)
    TIME = init_timer(CTS)
    crit = sample(["number", "dist"], Weights(weights))

    #it_time = @elapsed while LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
    while LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
        cranes_status = "idle"
        while cranes_status != "All cranes are busy" && cranes_status != "Next time period" && cranes_status != "LS is completed"
            cranes_status = randomizedConstructionHeuristic("dist", crit, tasks_by_w, prec, bj, LS, TIME, QC, CTS)
        end
        next_time_period(TIME, QC, CTS)
    end

    if check_solution(prec, LS, CTS) == true
        makespan = total_makespan(LS, CTS)
        best_LS = deepcopy(LS)
        best_makespan = makespan
        if makespan < M1
            #perform local_search_1
            gamma = rand(35:5:95)/100
            for ind in ["number", "minimal"]
                new_LS = remove_tasks(1-gamma, LS, CTS)
                check, TIME, QC = get_current_state(new_LS, CTS)
                while new_LS.tasks_left > 0 && TIME.period < TIME.horizon_plan
                    cranes_status = "idle"
                    while cranes_status != "All cranes are busy" && cranes_status != "Next time period" && cranes_status != "LS is completed"
                        cranes_status = deterministicConstructionHeuristic(ind, crit, tasks_by_w, prec, bj, new_LS, TIME, QC, CTS)
                    end
                    next_time_period(TIME, QC, CTS);
                end

                if check_solution(prec, new_LS, CTS) == true
                    new_makespan = total_makespan(new_LS, CTS)
                    if new_makespan < makespan
                        best_LS = deepcopy(new_LS)
                        best_makespan = new_makespan
                    end
                else
                    break
                end
            end
        end

        if best_makespan < M2
            #perform local_search_2
            r = remove_single_moves(prec, best_makespan, best_LS, CTS)
            remove_useless_travel(r)
            return((total_makespan(r, CTS), r, crit, makespan))
        else
            return((best_makespan, best_LS, crit, makespan))
        end
    else
        return((CTS.H, 0, crit, CTS.H))
    end
end

function GRASP_multi_thread(beta::Number, prec::Dict{Int, Array}, task_times::Array{Int, 2}, bj::Array{Int,1}, CTS::Constants)
    tasks_by_w = OperationalShippingProblem(beta, task_times, bj, CTS)
    worst_makespan = 0
    best_makespan = CTS.H
    best_LS = init_ls(CTS)
    iter = 0
    no_imp_it = 0
    #reactive weights
    weights = [0.5, 0.5]
    sort_crit = ["number", "dist"]
    crit_results = Dict{String, Tuple{Number, Int}}()
    for i in sort_crit
        crit_results[i] = (0,0)
    end
    # statistical filter
    M1 = CTS.H
    M2 = CTS.H
    sol_ini = Array{Int, 1}()
    sol_local = Array{Int, 1}()

    while iter < 15000 && no_imp_it < 1000
        #update indicator probabilities
        if iter%500 == 0
            sum_inv = 0
            for (key, value) in crit_results
                sum_inv += ((value[1]-worst_makespan)/(best_makespan-worst_makespan))^2
            end
            if sum_inv != 0
                for i = 1:2
                    weights[i] = ((crit_results[sort_crit[i]][1]-worst_makespan)/(best_makespan-worst_makespan))^2/sum_inv
                end
            end
        end

        s1 = @spawnat 2 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        s2 = @spawnat 3 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        s3 = @spawnat 4 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        s4 = @spawnat 5 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        s5 = @spawnat 6 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        s6 = @spawnat 7 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        s7 = @spawnat 8 GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        r = GRASP_single_thread(M1, M2, weights, tasks_by_w, prec, bj, CTS)
        iter += 8
        wait(s1)
        crit_results[s1[3]] = ((crit_results[s1[3]][1]*crit_results[s1[3]][end]+s1[1])/(crit_results[s1[3]][end]+1), crit_results[s1[3]][end]+1)
        if s1[1] < best_makespan
            best_makespan = s1[1]
            best_LS = deepcopy(s1[2])
            no_imp_it = 0
        elseif s1[1] > worst_makespan
            worst_makespan = s1[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        wait(s2)
        crit_results[s2[3]] = ((crit_results[s2[3]][1]*crit_results[s2[3]][end]+s2[1])/(crit_results[s2[3]][end]+1), crit_results[s2[3]][end]+1)
        if s2[1] < best_makespan
            best_makespan = s2[1]
            best_LS = deepcopy(s2[2])
            no_imp_it = 0
        elseif s2[1] > worst_makespan
            worst_makespan = s2[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        wait(s3)
        crit_results[s3[3]] = ((crit_results[s3[3]][1]*crit_results[s3[3]][end]+s3[1])/(crit_results[s3[3]][end]+1), crit_results[s3[3]][end]+1)
        if s3[1] < best_makespan
            best_makespan = s3[1]
            best_LS = deepcopy(s3[2])
            no_imp_it = 0
        elseif s3[1] > worst_makespan
            worst_makespan = s3[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        wait(s4)
        crit_results[s4[3]] = ((crit_results[s4[3]][1]*crit_results[s4[3]][end]+s4[1])/(crit_results[s4[3]][end]+1), crit_results[s4[3]][end]+1)
        if s4[1] < best_makespan
            best_makespan = s4[1]
            best_LS = deepcopy(s4[2])
            no_imp_it = 0
        elseif s4[1] > worst_makespan
            worst_makespan = s4[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        wait(s5)
        crit_results[s5[3]] = ((crit_results[s5[3]][1]*crit_results[s5[3]][end]+s5[1])/(crit_results[s5[3]][end]+1), crit_results[s5[3]][end]+1)
        if s5[1] < best_makespan
            best_makespan = s5[1]
            best_LS = deepcopy(s5[2])
            no_imp_it = 0
        elseif s5[1] > worst_makespan
            worst_makespan = s5[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        wait(s6)
        crit_results[s6[3]] = ((crit_results[s6[3]][1]*crit_results[s6[3]][end]+s6[1])/(crit_results[s6[3]][end]+1), crit_results[s6[3]][end]+1)
        if s6[1] < best_makespan
            best_makespan = s6[1]
            best_LS = deepcopy(s6[2])
            no_imp_it = 0
        elseif s6[1] > worst_makespan
            worst_makespan = s6[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        wait(s7)
        crit_results[s7[3]] = ((crit_results[s7[3]][1]*crit_results[s7[3]][end]+s7[1])/(crit_results[s7[3]][end]+1), crit_results[s7[3]][end]+1)
        if s7[1] < best_makespan
            best_makespan = s7[1]
            best_LS = deepcopy(s7[2])
            no_imp_it = 0
        elseif s7[1] > worst_makespan
            worst_makespan = s7[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end
        crit_results[r[3]] = ((crit_results[r[3]][1]*crit_results[r[3]][end]+r[1])/(crit_results[r[3]][end]+1), crit_results[r[3]][end]+1)
        if r[1] < best_makespan
            best_makespan = r[1]
            best_LS = deepcopy(r[2])
            no_imp_it = 0
        elseif r[1] > worst_makespan
            worst_makespan = r[1]
            no_imp_it += 1
        else
            no_imp_it += 1
        end

        #update statistical filter
        if iter < 400
            append!(sol_local, [s1[1], s2[1], s3[1], s4[1], s5[1], s6[1], s7[1], r[1]])
            append!(sol_ini, [s1[4], s2[4], s3[4], s4[4], s5[4], s6[4], s7[4], r[4]])
        elseif iter == 400
            sol_ratio = Array{Float64,1}()
            for i = 1:length(sol_ini)
                push!(sol_ratio, sol_ini[i]/sol_local[i])
            end
            M1 = median(sol_ratio)*best_makespan
            M2 = nquantile(sol_ratio, 4)[2]*best_makespan
        end
    end

    return(best_makespan, best_LS)
end

max_stars_repo_licenses: Zlib MIT

max_issues_repo_licenses: Zlib MIT

max_forks_repo_licenses: Zlib MIT


text:
#///////////////////////////////////////
#// File Name: ileqg.jl
#// Author: Haruki Nishimura (hnishimura@stanford.edu)
#// Date Created: 2020/10/28
#// Description: Iterative Linear-Exponential-Quadratic-Gaussian algorithm
#///////////////////////////////////////

using LinearAlgebra
using ForwardDiff
using Printf
using Random
using Distributions

"""
Simulate noiseless dynamics x_{t+1} = f(x_t, u_t) with u_array from x_0.
"""
function simulate_dynamics(problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                           x_0::Vector{Float64}, u_array::Vector{Vector{Float64}};
                           f_returns_jacobian=false)
    @assert problem.N == length(u_array)
    x_array = Vector{Vector{Float64}}(undef, problem.N + 1)
    x_array[1] = copy(x_0);
    if f_returns_jacobian
        A_array = Vector{Matrix{Float64}}(undef, problem.N)
        B_array = Vector{Matrix{Float64}}(undef, problem.N)
        for ii = 1 : problem.N
            x_array[ii + 1], A_array[ii], B_array[ii] =
                problem.f(x_array[ii], u_array[ii], true)
        end
        return x_array, A_array, B_array
    else
        for ii = 1 : problem.N
            x_array[ii + 1] = problem.f(x_array[ii], u_array[ii])
        end
        return x_array
    end
end;

"""
Simulate noisy dynamics x_{t+1} = f(x_t, u_t) + w_t with w_t ~ N(0, W(t)) and u_array from x_0.
"""
function simulate_dynamics(problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                           x_0::Vector{Float64}, u_array::Vector{Vector{Float64}},
                           rng::AbstractRNG)
    @assert problem.N == length(u_array)
    x_array = Vector{Vector{Float64}}(undef, problem.N + 1)
    x_array[1] = x_0;
    for ii = 1 : problem.N
        d = MvNormal(zeros(length(x_0)), problem.W(ii - 1));
        x_array[ii + 1] = problem.f(x_array[ii], u_array[ii]) + rand(rng, d)
    end
    return x_array
end;

"""
Simulate noiseless dynamics x_{t+1} = f(x_t, u_t) with x_array and state feedback policy
π(x_t) = L_t*(x_t - ̄x_t) + l_t.
"""
function simulate_dynamics(problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                           x_array::Vector{Vector{Float64}},
                           l_array::Vector{Vector{Float64}},
                           L_array::Vector{Matrix{Float64}};
                           f_returns_jacobian=false)
    @assert problem.N == length(l_array) && problem.N == length(L_array)
    u_array = Vector{Vector{Float64}}(undef, problem.N)
    x_array_new = Vector{Vector{Float64}}(undef, problem.N + 1)
    x_array_new[1] = copy(x_array[1]);
    if f_returns_jacobian
        A_array = Vector{Matrix{Float64}}(undef, problem.N)
        B_array = Vector{Matrix{Float64}}(undef, problem.N)
        for ii = 1 : problem.N
            u_array[ii] = l_array[ii] + L_array[ii]*(x_array_new[ii] - x_array[ii])
            x_array_new[ii + 1], A_array[ii], B_array[ii] =
                problem.f(x_array_new[ii], u_array[ii], true)
        end
        return x_array_new, u_array, A_array, B_array
    else
        for ii = 1 : problem.N
            u_array[ii] = l_array[ii] + L_array[ii]*(x_array_new[ii] - x_array[ii])
            x_array_new[ii + 1] = problem.f(x_array_new[ii], u_array[ii])
        end
        return x_array_new, u_array
    end
end

"""
Simulate noisy dynamics x_{t+1} = f(x_t, u_t) + w_t with w_t ~ N(0, W(t)), x_array and state feedback
policy π(x_t) = L_t*(x_t - ̄x_t) + l_t.
"""
function simulate_dynamics(problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                           x_array::Vector{Vector{Float64}},
                           l_array::Vector{Vector{Float64}},
                           L_array::Vector{Matrix{Float64}},
                           rng::AbstractRNG)
    @assert problem.N == length(l_array) && problem.N == length(L_array)
    u_array = Vector{Vector{Float64}}(undef, problem.N)
    x_array_new = Vector{Vector{Float64}}(undef, problem.N + 1)
    x_array_new[1] = copy(x_array[1]);
    for ii = 1 : problem.N
        d = MvNormal(zeros(length(x_array_new[ii])), problem.W(ii - 1));
        u_array[ii] = l_array[ii] + L_array[ii]*(x_array_new[ii] - x_array[ii])
        x_array_new[ii + 1] = problem.f(x_array_new[ii], u_array[ii]) + rand(rng, d)
    end
    return x_array_new, u_array
end

"""
Integrate cost given state and control trajectories
"""
function integrate_cost(problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                        x_array::Vector{Vector{Float64}}, u_array::Vector{Vector{Float64}})
    @assert problem.N == length(u_array) && problem.N + 1 == length(x_array)
    cost = 0.0;
    for ii = 1 : problem.N
        cost += problem.c(ii - 1, x_array[ii], u_array[ii])
    end
    cost += problem.h(x_array[end])
    return cost
end;

"""
    ILEQGSolver(problem::FiniteHorizonRiskSensitiveOptimalControlProblems, kwargs...)

iLQG and iLEQG Solver for `problem`.

# Optional Keyword Arguments
- `μ_min::Float64` -- minimum value for Hessian regularization parameter `μ` (> 0).
  Default: `1e-6`.
- `Δ_0::Float64` -- minimum multiplicative modification factor (> 0) for `μ`.
  Default: `2.0`.
- `λ::Float64` -- multiplicative modification factor in (0, 1) for line search
  step size `ϵ`. Default: `0.5`.
- `d::Float64` -- convergence error norm threshold (> 0). If the maximum l2
  norm of the change in nominal control over the horizon is less than `d`, the
  solver is considered to be converged. Default: `1e-2`.
- `iter_max::Int64` -- maximum iteration number. Default: 100.
- `ϵ_init::Float64` -- initial step size in (`ϵ_min`, 1] to start the
  backtracking line search with. If `adaptive_ϵ_init` is `true`, then this
  value is overridden by the solver's adaptive initialization functionality
  after the first iLEQG iteration. If `adaptive_ϵ_init` is `false`, the
  specified value of `ϵ_init` is used across all the iterations as the initial
  step size. Default:`1.0`.
- `adaptive_ϵ_init::Bool` -- if `true`, `ϵ_init` is adaptively changed based on
  the last step size `ϵ` of the previous iLEQG iteration. Default: `false`.
   - If the first line search iterate `ϵ_init_prev` in the previous iLEQG
     iteration is successful, then `ϵ_init` for the next iLEQG iteration is set
     to `ϵ_init = ϵ_init_prev / λ` so that the initial line search step increases.
   - Otherwise `ϵ_init = ϵ_last` where `ϵ_last` is the line search step accepted
     in the previous iLEQG iteration.
- `ϵ_min::Float64` -- minimum value of step size `ϵ` to terminate the line
  search. When `ϵ_min` is reached, the last candidate nominal trajectory is accepted
  regardless of the Armijo condition and the current iLEQG iteration is
  finished. Default: `1e-6`.
- `f_returns_jacobian::Bool` -- if `true`, Jacobian matrices of the dynamics function
  are user-provided. This can reduce computation time since automatic
  differentiation is not used. Default: `false`.
"""
mutable struct ILEQGSolver
    μ_min::Float64   # Minimum damping parameter for regularization
    μ::Float64       # Damping parameter for regularization
    Δ_0::Float64     # Minimum Modification factor for μ
    Δ::Float64       # Modification factor for μ
    λ::Float64       # Multiplicative factor for line search step parameter
    d::Float64       # Convergence error norm thresholds
    iter_max::Int64  # Maximum iteration
    ϵ_init_auto::Bool # Automatic initialization of ϵ_init from the previous iLEQG iteration.
    ϵ_init::Float64  # Initial step size for backtracking line search
    ϵ_min::Float64   # Minimum step size for backtracking line search
    f_returns_jacobian::Bool

    x_array::Vector{Vector{Float64}} # Nominal state trajectory
    l_array::Vector{Vector{Float64}} # Linear control schedule
    L_array::Vector{Matrix{Float64}} # Feedback gain schedule

    A_array::Union{Nothing, Vector{Matrix{Float64}}} # Sequence of Jacobians [dx_1/dx_0, ..., dx_{T}/dx_{T-1}]
    B_array::Union{Nothing, Vector{Matrix{Float64}}} # Sequence of Jacobians [dx_1/du_0, ..., dx_{T}/du_0]

    value_current::Float64 # current value (cost-to-go)
    iter_current::Int64 # current iteration number
    d_current::Float64  # current error norm
    ϵ_history::Array{Tuple{Float64, Float64}} # Array of (ϵ, -value_improvement)
    ϵ_init_init::Float64 # initial value of ϵ_init
end;

function ILEQGSolver(problem::FiniteHorizonRiskSensitiveOptimalControlProblem;
                     μ_min=1e-6, Δ_0=2.0, λ=0.5, d=1e-2, iter_max=100,
                     ϵ_init=1.0, adaptive_ϵ_init=false,
                     ϵ_min=1e-6, f_returns_jacobian=false)
    @assert 0 < λ < 1 "λ has to be in (0, 1)"
    @assert d > 0 "d > 0 is necessary"
    @assert μ_min > 0 "μ_min > 0 is necessary"
    @assert Δ_0 > 0 "Δ_0 > 0 is necessary"
    @assert 0 < ϵ_init <= 1 "ϵ_init has to be in (0, 1]"
    @assert ϵ_init > ϵ_min "ϵ_init > ϵ_min is necessary"
    @assert 0 < ϵ_min < 1 "ϵ_min has to be in (0, 1)"

    x_array = Vector{Vector{Float64}}(undef, problem.N + 1)
    l_array = Vector{Vector{Float64}}(undef, problem.N)
    L_array = Vector{Matrix{Float64}}(undef, problem.N)
    return ILEQGSolver(μ_min, μ_min, Δ_0, Δ_0, λ, d, iter_max, adaptive_ϵ_init, ϵ_init, ϵ_min, f_returns_jacobian,
                       x_array, l_array, L_array, nothing, nothing, Inf, 0, Inf, Tuple{Float64, Float64}[], ϵ_init)
end;

"""
Initialize iLEQG solver
"""
function initialize!(ileqg::ILEQGSolver, problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                     x_0::Vector{Float64}, u_array::Vector{Vector{Float64}}, θ::Float64)
    ileqg.μ, ileqg.Δ = 0.0, ileqg.Δ_0;
    ileqg.d_current = Inf;
    ileqg.iter_current = 0;
    ileqg.ϵ_init = ileqg.ϵ_init_init
    ileqg.ϵ_history = Tuple{Float64, Float64}[];
    if ileqg.f_returns_jacobian
        ileqg.x_array, ileqg.A_array, ileqg.B_array =
            simulate_dynamics(problem, x_0, u_array, f_returns_jacobian=true)
    else
        ileqg.x_array = simulate_dynamics(problem, x_0, u_array, f_returns_jacobian=false);
        ileqg.A_array, ileqg.B_array = nothing, nothing
    end
    ileqg.l_array = copy(u_array);
    x_dim, u_dim = length(ileqg.x_array[1]), length(ileqg.l_array[1])
    for ii = 1 : problem.N
        ileqg.L_array[ii] = zeros(u_dim, x_dim)
    end
    approx_result = approximate_model(problem, ileqg.l_array, ileqg.x_array);
    dp_result = solve_approximate_dp(approx_result, ileqg.L_array, θ=θ, μ=ileqg.μ);
    ileqg.value_current = dp_result.s_array[1];
end

"""
Model Approximation
"""
struct ApproximationResult
    q_array::Vector{Float64}
    q_vec_array::Vector{Vector{Float64}}
    Q_array::Vector{Matrix{Float64}}
    r_array::Vector{Vector{Float64}}
    R_array::Vector{Matrix{Float64}}
    P_array::Vector{Matrix{Float64}}
    A_array::Vector{Matrix{Float64}}
    B_array::Vector{Matrix{Float64}}
    W_array::Vector{Matrix{Float64}}
end

"""
Linearize the dynamics and quadratize the costs around the nominal trajectory.
"""
function approximate_model(problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                           u_array::Vector{Vector{Float64}},
                           x_array::Vector{Vector{Float64}},
                           A_array_input::Union{Nothing, Vector{Matrix{Float64}}}=nothing,
                           B_array_input::Union{Nothing, Vector{Matrix{Float64}}}=nothing)
    f, c, h, W = problem.f, problem.c, problem.h, problem.W;
    # linearization & quadratization functions
    fx(x, u) = ForwardDiff.jacobian(a -> f(a, u), x);
    fu(x, u) = ForwardDiff.jacobian(b -> f(x, b), u);
    cx(k, x, u) = ForwardDiff.gradient(a -> c(k, a, u), x);
    cu(k, x, u) = ForwardDiff.gradient(b -> c(k, x, b), u);
    cux(k, x, u) = ForwardDiff.jacobian(a -> ForwardDiff.gradient(b -> c(k, a, b), u), x);
    cxx(k, x, u) = Symmetric(ForwardDiff.hessian(a -> c(k, a, u), x));
    cuu(k, x, u) = Symmetric(ForwardDiff.hessian(b -> c(k, x, b), u));
    hx(x) = ForwardDiff.gradient(a -> h(a), x);
    hxx(x) = Symmetric(ForwardDiff.hessian(a -> h(a), x));

    # compute approximations
    @assert problem.N == length(u_array)
    N = problem.N;
    if !isnothing(A_array_input)
        @assert length(A_array_input) == N
    end
    if !isnothing(B_array_input)
        @assert length(B_array_input) == N
    end
    q_array = Vector{Float64}(undef, N + 1)
    q_vec_array = Vector{Vector{Float64}}(undef, N + 1)
    Q_array = Vector{Matrix{Float64}}(undef, N + 1)
    r_array = Vector{Vector{Float64}}(undef, N)
    R_array = Vector{Matrix{Float64}}(undef, N)
    P_array = Vector{Matrix{Float64}}(undef, N)
    A_array = Vector{Matrix{Float64}}(undef, N)
    B_array = Vector{Matrix{Float64}}(undef, N)
    W_array = Vector{Matrix{Float64}}(undef, N)
    # Threads.@threads for ii = 1 : N
    for ii = 1 : N
        @inbounds x, u = x_array[ii], u_array[ii];
        @inbounds q_array[ii] = c(ii - 1, x, u);
        @inbounds q_vec_array[ii] = cx(ii - 1, x, u);
        @inbounds Q_array[ii] = cxx(ii - 1, x, u);
        @inbounds r_array[ii] = cu(ii - 1, x, u);
        @inbounds R_array[ii] = cuu(ii - 1, x, u);
        @inbounds P_array[ii] = cux(ii - 1, x, u);
        if isnothing(A_array_input)
            @inbounds A_array[ii] = fx(x, u);
        else
            @inbounds A_array[ii] = A_array_input[ii] # copy(A_array_input[ii])
        end
        if isnothing(B_array_input)
            @inbounds B_array[ii] = fu(x, u);
        else
            @inbounds B_array[ii] = B_array_input[ii] # copy(B_array_input[ii])
        end
        @inbounds W_array[ii] = W(ii - 1);
    end
    @inbounds q_array[N + 1] = h(x_array[N + 1]);
    @inbounds q_vec_array[N + 1] = hx(x_array[N + 1]);
    @inbounds Q_array[N + 1] = hxx(x_array[N + 1]);

    approx_result = ApproximationResult(q_array, q_vec_array, Q_array,
                                        r_array, R_array, P_array,
                                        A_array, B_array, W_array)
    return approx_result
end;

"""
Dynamic Programming Result
"""
struct DynamicProgrammingResult
    s_array::Vector{Float64}
    s_vec_array::Vector{Vector{Float64}}
    S_array::Vector{Matrix{Float64}}
    g_array::Vector{Vector{Float64}}
    G_array::Vector{Matrix{Float64}}
    H_array::Vector{Matrix{Float64}}
end

"""
Compute approximate dynamic programming solution using Riccati(-like) recursiong with
risk-sensitivity θ, while optimizing the policy parameters L_array and dl_array.
"""
function solve_approximate_dp!(ileqg::ILEQGSolver, approx_result::ApproximationResult,
                               verbose::Bool; θ::Float64)
    N = length(approx_result.W_array);
    s_array = Vector{Float64}(undef, N + 1)
    s_vec_array = Vector{Vector{Float64}}(undef, N + 1)
    S_array = Vector{Matrix{Float64}}(undef, N + 1)
    g_array = Vector{Vector{Float64}}(undef, N)
    G_array = Vector{Matrix{Float64}}(undef, N)
    H_array = Vector{Matrix{Float64}}(undef, N)

    # terminal condition
    s_array[end] = approx_result.q_array[end];
    s_vec_array[end] = approx_result.q_vec_array[end];
    S_array[end] = Symmetric(approx_result.Q_array[end]);

    # riccati-like recursion, computation of optimal gains with regularization
    dl_array_new = Vector{Vector{Float64}}(undef, N);
    all_hessians_psd = false
    while !all_hessians_psd
        for ii in Iterators.reverse(1 : N)
            @inbounds q, q_vec, Q = approx_result.q_array[ii], approx_result.q_vec_array[ii], approx_result.Q_array[ii];
            @inbounds r, R, P = approx_result.r_array[ii], approx_result.R_array[ii], approx_result.P_array[ii];
            @inbounds A, B = approx_result.A_array[ii], approx_result.B_array[ii];
            @inbounds W = approx_result.W_array[ii];
            @inbounds M = Symmetric(inv(W) - θ.*S_array[ii + 1]);
            @assert isposdef(M) "M = $(M): (inv(W_{ii}) - θ*S_{ii + 1}) is not PSD at ii = $(ii)"
            @inbounds D = Matrix(1.0I, size(W)) + θ.*S_array[ii + 1]/M
            @inbounds g_array[ii] = r + B'*D*s_vec_array[ii + 1];
            @inbounds G_array[ii] = P + B'*(D*S_array[ii + 1])*A;
            @inbounds H_array[ii] = R + B'*(D*S_array[ii + 1])*B + ileqg.μ*Matrix(1.0I, size(R));
            @inbounds H_array[ii] = Symmetric(H_array[ii])
            if !isposdef(H_array[ii])
                increase_μ_and_Δ!(ileqg)
                if verbose
                    println("------Hessian not PSD at ii = $(ii). Increasing μ to $(ileqg.μ) and Δ to $(ileqg.Δ)")
                end
                break;
            end
            @inbounds L = -H_array[ii]\G_array[ii];
            @inbounds ileqg.L_array[ii] = L;
            @inbounds dl = -H_array[ii]\g_array[ii];
            @inbounds dl_array_new[ii] = dl;
            @inbounds s_array[ii] = q + s_array[ii + 1] + 0.5*dl'*H_array[ii]*dl + dl'*g_array[ii];
            if θ == 0.0 # iLQG
                @inbounds s_array[ii] += 0.5*tr(W*S_array[ii + 1])
            else
                @inbounds s_array[ii] += θ/2*s_vec_array[ii + 1]'/M*s_vec_array[ii + 1] - 1/(2*θ)*logdet(W*M);
            end
            @inbounds s_vec_array[ii] = q_vec + A'*D*s_vec_array[ii + 1] + L'*H_array[ii]*dl + L'*g_array[ii] + G_array[ii]'*dl;
            @inbounds S_array[ii] = Q + A'*D*S_array[ii + 1]*A + L'*H_array[ii]*L + L'*G_array[ii] + G_array[ii]'*L;
            @inbounds S_array[ii] = Symmetric(S_array[ii])
            if ii == 1
                all_hessians_psd = true
            end
        end
        if all_hessians_psd
            if verbose
                println("------Approximate dynamic programming solved.")
            end
        end
    end
    dp_result = DynamicProgrammingResult(s_array, s_vec_array, S_array,
                                         g_array, G_array, H_array);

    return dp_result, dl_array_new;
end

"""
Compute approximate dynamic programming solution using Riccati(-like) recursion under the policy
π(x_t) = L_t*(x_t - ̄x_t) + l_t + dl_t, with risk-sensitivity θ.
"""
function solve_approximate_dp(approx_result::ApproximationResult,
                              L_array::Vector{Matrix{Float64}},
                              dl_array::Union{Nothing, Vector{Vector{Float64}}}=nothing;
                              θ::Float64, μ::Float64)
    N = length(approx_result.W_array);
    @assert N == length(L_array)
    if !isnothing(dl_array)
        @assert N == length(dl_array)
    end
    s_array = Vector{Float64}(undef, N + 1)
    s_vec_array = Vector{Vector{Float64}}(undef, N + 1)
    S_array = Vector{Matrix{Float64}}(undef, N + 1)
    g_array = Vector{Vector{Float64}}(undef, N)
    G_array = Vector{Matrix{Float64}}(undef, N)
    H_array = Vector{Matrix{Float64}}(undef, N)

    # terminal condition
    s_array[end] = approx_result.q_array[end];
    s_vec_array[end] = approx_result.q_vec_array[end];
    S_array[end] = Symmetric(approx_result.Q_array[end]);

    # riccati-like recursion
    for ii in Iterators.reverse(1 : N)
        @inbounds q, q_vec, Q = approx_result.q_array[ii], approx_result.q_vec_array[ii], approx_result.Q_array[ii];
        @inbounds r, R, P = approx_result.r_array[ii], approx_result.R_array[ii], approx_result.P_array[ii];
        @inbounds A, B = approx_result.A_array[ii], approx_result.B_array[ii];
        @inbounds W = approx_result.W_array[ii];
        @inbounds M = Symmetric(inv(W) - θ.*S_array[ii + 1]);
        @assert isposdef(M) "M = $(M): (inv(W_{ii}) - θ*S_{ii + 1}) is not PSD at ii = $(ii)"
        @inbounds D = Matrix(1.0I, size(W)) + θ.*S_array[ii + 1]/M
        @inbounds g_array[ii] = r + B'*D*s_vec_array[ii + 1];
        @inbounds G_array[ii] = P + B'*(D*S_array[ii + 1])*A;
        @inbounds H_array[ii] = R + B'*(D*S_array[ii + 1])*B + μ*Matrix(1.0I, size(R));
        @inbounds H_array[ii] = Symmetric(H_array[ii])
        @inbounds L = L_array[ii];
        if isnothing(dl_array)
            dl = zeros(size(L, 1))
        else
            @inbounds dl = dl_array[ii];
        end
        @inbounds s_array[ii] = q + s_array[ii + 1] + 0.5*dl'*H_array[ii]*dl + dl'*g_array[ii];
        if θ == 0.0 # iLQG
            @inbounds s_array[ii] += 0.5*tr(W*S_array[ii + 1])
        else
            @inbounds s_array[ii] += θ/2*s_vec_array[ii + 1]'/M*s_vec_array[ii + 1] - 1/(2*θ)*logdet(W*M);
        end
        @inbounds s_vec_array[ii] = q_vec + A'*D*s_vec_array[ii + 1] + L'*H_array[ii]*dl + L'*g_array[ii] + G_array[ii]'*dl;
        @inbounds S_array[ii] = Q + A'*D*S_array[ii + 1]*A + L'*H_array[ii]*L + L'*G_array[ii] + G_array[ii]'*L;
        @inbounds S_array[ii] = Symmetric(S_array[ii])
    end
    dp_result = DynamicProgrammingResult(s_array, s_vec_array, S_array,
                                         g_array, G_array, H_array)
    return dp_result
end;

"""
Increase μ and Δ of iLEQG solver
"""
function increase_μ_and_Δ!(ileqg::ILEQGSolver)
    ileqg.Δ = max(ileqg.Δ_0, ileqg.Δ*ileqg.Δ_0);
    ileqg.μ = max(ileqg.μ_min, ileqg.μ*ileqg.Δ);
end;

"""
Decrease μ and Δ of iLEQG solver
"""
function decrease_μ_and_Δ!(ileqg::ILEQGSolver)
    ileqg.Δ = min(1/ileqg.Δ_0, ileqg.Δ/ileqg.Δ_0);
    new_μ_cand = ileqg.μ*ileqg.Δ;
    if new_μ_cand >= ileqg.μ_min;
        ileqg.μ = new_μ_cand;
    else
        ileqg.μ = 0.0;
    end
end;

"""
Perform Backtracking Line Search
"""
function line_search!(ileqg::ILEQGSolver, problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                      dl_array_new::Vector{Vector{Float64}}, θ::Float64, verbose=true)

    expected_cost_current = ileqg.value_current;
    if verbose
        cost_str = @sprintf "%3.3f" expected_cost_current
        println("----Current expected cost: $(cost_str)")
    end
    ϵ = ileqg.ϵ_init
    count = 0;
    while true
        count += 1;
        if verbose
            println("----Performing line search with ϵ == $(ϵ)")
        end
        l_array_new = ileqg.l_array .+ ϵ.*dl_array_new;
        if ileqg.f_returns_jacobian
            x_array_new, u_array_new, A_array_new, B_array_new =
                simulate_dynamics(problem, ileqg.x_array, l_array_new, ileqg.L_array,
                                  f_returns_jacobian=true)
            approx_result_new = approximate_model(problem, u_array_new, x_array_new,
                                                  A_array_new, B_array_new);
        else
            x_array_new, u_array_new = simulate_dynamics(problem, ileqg.x_array,
                                                         l_array_new, ileqg.L_array,
                                                         f_returns_jacobian=false)
            approx_result_new = approximate_model(problem, u_array_new, x_array_new);
        end
        dp_result_new =
        try
            solve_approximate_dp(approx_result_new, ileqg.L_array,
                                 θ=θ, μ=ileqg.μ)
        catch e
            nothing;
        end
        if isnothing(dp_result_new)
            ϵ *= ileqg.λ;
            if verbose
                println("----Approximate DP not PosDef. Re-doing with ϵ == $(ϵ)");
            end
            continue;
        end
        expected_cost_new = dp_result_new.s_array[1];
        push!(ileqg.ϵ_history, (ϵ, expected_cost_new - expected_cost_current));
        if expected_cost_new ≈ expected_cost_current || expected_cost_new < expected_cost_current
            ileqg.d_current = maximum(norm.(ileqg.l_array .- u_array_new));
            if verbose
                cost_str = @sprintf "%3.3f" expected_cost_new;
                err_norm_str = @sprintf "%3.3f" ileqg.d_current
                println("------Accepted. New expected cost: $(cost_str). d == $(err_norm_str)")
            end
            ileqg.value_current = expected_cost_new;
            ileqg.x_array = x_array_new;
            ileqg.l_array = u_array_new;
            if ileqg.f_returns_jacobian
                ileqg.A_array = approx_result_new.A_array
                ileqg.B_array = approx_result_new.B_array
            else
                ileqg.A_array = nothing
                ileqg.B_array = nothing
            end
            break;
        else
            ϵ *= ileqg.λ;
            if ϵ < ileqg.ϵ_min
                ileqg.d_current = maximum(norm.(ileqg.l_array .- u_array_new));
                if verbose
                    cost_str = @sprintf "%3.3f" expected_cost_new;
                    err_norm_str = @sprintf "%3.3f" ileqg.d_current
                    println("------New expected cost: $(cost_str). Terminating as ϵ_min reached. d == $(err_norm_str)")
                end
                ileqg.value_current = expected_cost_new;
                ileqg.x_array = x_array_new;
                ileqg.l_array = u_array_new;
                if ileqg.f_returns_jacobian
                    ileqg.A_array = approx_result_new.A_array
                    ileqg.B_array = approx_result_new.B_array
                else
                    ileqg.A_array = nothing
                    ileqg.B_array = nothing
                end
                break;
            elseif verbose
                cost_str = @sprintf "%3.3f" expected_cost_new;
                println("------New expected cost: $(cost_str). Decreasing ϵ to $(ϵ)")
            end
        end
    end
    if ileqg.ϵ_init_auto
        if count == 1
            ileqg.ϵ_init = min(ileqg.ϵ_init_init, ϵ/ileqg.λ)
        else
            while ϵ < ileqg.ϵ_min
                ϵ = ϵ/ileqg.λ;
            end
            ileqg.ϵ_init = ϵ
        end
    end
end

"""
single iteration of iLEQG
"""
function step!(ileqg::ILEQGSolver, problem::FiniteHorizonRiskSensitiveOptimalControlProblem, θ::Float64, verbose=true)
    ileqg.iter_current += 1;
    if verbose
        println("--ILEQG iteration $(ileqg.iter_current)")
        println("----Approximating model around current l_array and x_array")
    end
    approx_result = approximate_model(problem, ileqg.l_array, ileqg.x_array,
                                      ileqg.A_array, ileqg.B_array)

    if verbose
        println("----Solving approximate dynamic programming")
    end
    ~, dl_array =
        solve_approximate_dp!(ileqg, approx_result, verbose, θ=θ);
    line_search!(ileqg, problem, dl_array, θ, verbose);
end

"""
    solve!(ileqg::ILEQGSolver, problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
    x_0::Vector{Float64}, u_array::Vector{Vector{Float64}}; θ::Float64, verbose=true)

Given `problem`, and `ileqg` solver, solve iLQG (if `θ == 0`) or iLEQG (if `θ > 0`)
with current state `x_0` and nominal control schedule `u_array = [u_0, ..., u_{N-1}]`.

# Return Values (Ordered)
- `x_array::Vector{Vector{Float64}}` -- nominal state trajectory `[x_0,...,x_N]`.
- `l_array::Vector{Vector{Float64}}` -- nominal control schedule `[l_0,...,l_{N-1}]`.
- `L_array::Vector{Matrix{Float64}}` -- feedback gain schedule `[L_0,...,L_{N-1}]`.
- `value::Float64` -- optimal cost-to-go (i.e. value) found by the solver.
- `ϵ_history::Vector{Float64}` -- history of line search step sizes used during
  the iLEQG iteration. Mainly for debugging purposes.

# Notes
- Returns a time-varying affine state-feedback policy `π_k` of the form
  `π_k(x) = L_k(x - x_k) + l_k`.
"""
function solve!(ileqg::ILEQGSolver,
                problem::FiniteHorizonRiskSensitiveOptimalControlProblem,
                x_0::Vector{Float64}, u_array::Vector{Vector{Float64}};
                θ::Float64, verbose=true)
    initialize!(ileqg, problem, x_0, u_array, θ);
    while true
        step!(ileqg, problem, θ, verbose);
        if ileqg.d > ileqg.d_current && ileqg.μ <= ileqg.μ_min
            if verbose
                err_norm_str = @sprintf "%3.3f" ileqg.d_current
                println("ILEQG Converged. d == $(err_norm_str)")
            end
            break;
        elseif ileqg.iter_current == ileqg.iter_max
            if verbose
                println("Maximum iteration number reached.")
            end
            break;
        end
    end
    x_array, l_array, L_array = copy(ileqg.x_array), copy(ileqg.l_array), copy(ileqg.L_array)
    ϵ_history = copy(ileqg.ϵ_history)
    value = ileqg.value_current;
    return x_array, l_array, L_array, value, ϵ_history
end


text:
################### stepwise regression with replacement (SRR) #################
# one-lookahead greedy two-step algorithm, based on stepwise regression algorithms
function srr(A::AbstractMatOrFac, b::AbstractVector, k::Int, δ::Real = 1e-12,
                                            x = spzeros(eltype(A), size(A, 2));
                                            maxiter = 4k, initialization::Int = 1,
                                            l::Int = 1)
    P = StepwiseRegression(A, b, x.nzind)
    # initialize support with k atoms using
    if initialization == 1 # oblivious algorithm
        oblivious_acquisition!(P, x, k-nnz(x)) # initialize with k largest inner products
    elseif initialization == 2 # forward regression
        for _ in 1:k
            update!(P, x)
        end
    elseif initialization == 3 # random indices
        random_acquisition!(P, x, k-nnz(x)) # initialize with k largest inner products
    end
    resnorm = norm(residual!(P, x))
    for i in 1:maxiter
        oldnorm = resnorm
        for _ in 1:l
            forward_step!(P, x, 0, 0) || break # add unless we found a solution
        end
        while nnz(x) > k # delete until we have an active set of size k
            backward_step!(P, x, Inf, Inf)
        end
        resnorm = norm(residual!(P, x))
        if resnorm ≤ δ || oldnorm ≤ resnorm # until convergence or stationarity
            break
        end
    end
    return x
end

@inline function argmaxinner!(P::StepwiseRegression, k::Int)
    Ar = P.A' * P.r
    @. Ar = abs(Ar)
    partialsortperm(Ar, 1:k, rev = true)
end

################################## Subspace Pursuit ############################
struct SubspacePursuit{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T}} <: AbstractMatchingPursuit{T}
    A::AT
    b::B
    k::Int # maximum number of non-zeros

    # temporary storage
    r::B # residual
    Ar::B # inner products between measurement matrix and residual
    Ai::AT # space for A[:, x.nzind] and its qr factorization
end
const SP = SubspacePursuit

function SP(A::AbstractMatOrFac, b::AbstractVector, k::Integer)
    2k > length(b) && error("2k = $(2k) > $(length(b)) = length(b) is invalid for Subspace Pursuit")
    n, m = size(A)
    T = eltype(A)
    r, Ar = copy(b), zeros(T, m)
    Ai = zeros(T, (n, 2k))
    SP(A, b, k, r, Ar, Ai)
end

# adds indices of k atoms with largest inner products with residual
# solves least-squares problem on expanded support
# could use threshold on P.Ar for adaptive stopping
# this is oblivious algorithm applied to residual
function sp_acquisition!(P::SP, x, k::Int = P.k)
    residual!(P, x)
    i = argmaxinner!(P, k)
    @. x[i] = NaN
    solve!(P, x)
end
# IDEA: could pre-allocate nz arrays to be of length 2K
# IDEA: could add efficient qr updating
function update!(P::SP, x::AbstractVector = spzeros(size(P.A, 2)), ε::Real = 0.)
    nnz(x) == P.k || throw("nnz(x) = $(nnz(x)) ≠ $(P.k) = k")
    sp_acquisition!(P, x)
    i = partialsortperm(abs.(x.nzval), 1:nnz(x)-P.k) # find the smallest atoms
    sort!(i) # for deleteat! to work, indices need to be sorted
    deleteat!(x.nzind, i)
    deleteat!(x.nzval, i)
    solve!(P, x) # optimize all active atoms
end

# calculates k-sparse approximation to Ax = b via subspace pursuit
# could also stop if indices are same between iterations
function sp(A::AbstractMatOrFac, b::AbstractVector, k::Int, δ::Real = 1e-12; maxiter = 16k)
    P = SP(A, b, k)
    x = spzeros(size(A, 2))
    sp_acquisition!(P, x, P.k)
    resnorm = norm(residual!(P, x))
    for i in 1:maxiter
        oldnorm = resnorm
        update!(P, x)
        resnorm = norm(residual!(P, x))
        if resnorm ≤ δ || oldnorm ≤ resnorm
            break
        end
    end
    return x
end

# ordinary least squares solve
@inline function solve!(P::SP, x::AbstractVector, b::AbstractVector = P.b)
    F = factorize!(P, x)
    ldiv!(x.nzval, F, b)     # optimize all active atoms
end

############################# OMP with replacement #############################
struct OMPR{T, AT<:AbstractMatOrFac{T}, B<:AbstractVector{T}, FT} <: AbstractMatchingPursuit{T}
    A::AT
    b::B
    k::Int
    # l # number of atoms to replace

    # temporary storage
    r::B # residual
    Ar::B # inner products between measurement matrix and residual
    QA::AT # (k, m)
    AiQR::FT # updatable QR factorization of Ai
end

# similar to RMP, or OMP?
function OMPR(A::AbstractMatOrFac, b::AbstractVector, k::Int)
    n, m = size(A)
    T = eltype(A)
    r, Ar = zeros(T, n), zeros(T, m)
    AiQR = UpdatableQR(T, n, n) # initializing empty qr factorization with maximum rank n
    QA = zeros(T, (n, m))
    OMPR(A, b, k, r, Ar, QA, AiQR)
end

# η is stepsize
function update!(P::OMPR, x::AbstractVector, η::Real = 1.)
    nnz(x) == P.k || throw("nnz(x) = $(nnz(x)) ≠ $(P.k) = k")
    residual!(P, x)
    copy!(P.Ar, x)
    mul!(P.Ar, P.A', P.r, η, 1)
    begin # argmax of abs(P.Ar) over i not in x.nzind     # argmax( [abs(ar) for ar in P.Ar]
        m = 0.
        i = 0
        for (j, Arj) in enumerate(P.Ar)
            if j ∉ x.nzind
                fArj = abs(Arj)
                if fArj > m
                    m = fArj
                    i = j
                end
            end
        end
    end

    if i == 0
        return x
    end

    # add non-zero index to active set
    x[i] = NaN
    qr_i = findfirst(==(i), x.nzind) # index in qr where new atom should be added

    # least-squares solve for active atoms (this would be two-stage algorithm)
    # ldiv!(x.nzval, P.AiQR, P.b)

    # gradient descent for active atoms
    @. x.nzval = P.Ar[x.nzind]

    j = argmin(abs, x.nzval) # delete index
    deleteat!(x.nzind, j) # delete the value from array
    deleteat!(x.nzval, j)

    if qr_i ≠ j # update qr factorization using Givens rotations
        a = @view P.A[:, i]
        add_column!(P.AiQR, a, qr_i)
        remove_column!(P.AiQR, j)
    end

    # least-squares solve for active atoms
    ldiv!!(x.nzval, P.AiQR, P.b, P.r)
    return x
end

# k is desired sparsity level
# l is cardinality of maximum replacement per iteration, l = k corresponds to sp
function ompr(A::AbstractMatOrFac, b::AbstractVector, k::Int, δ::Real,
                                x = spzeros(size(A, 2)); maxiter = size(A, 1))
    P = OMPR(A, b, k)
    if nnz(x) < P.k # make sure support set is of size k
        @. x = 0
        dropzeros!(x)
        oblivious_acquisition!(P, x, P.k)
    end
    resnorm = norm(residual!(P, x))
    for i in 1:maxiter
        oldnorm = resnorm
        update!(P, x)
        resnorm = norm(residual!(P, x))
        if resnorm ≤ δ || oldnorm ≤ resnorm
            break
        end
    end
    return x
end


text:
"""
This file is for Gym acrobot environment https://github.com/openai/gym/tree/master/gym/envs/classic_control.
These are based on https://github.com/WilsonWangTHU/mbbl
"""

using LyceumBase
using Shapes
using Distributions
using Base: @propagate_inbounds

struct Acrobot{D,SP,OP,AP,RP} <: AbstractEnvironment
    state::D
    obs::D
    action::D
    time::D
    nextobs::D
    ctrlrange::D

    sp::SP
    op::OP
    ap::AP
    rp::RP
    function Acrobot{T}() where T<:AbstractFloat
        nobs = 6
        nstate = 4
        naction = 1

        ctrlrange = zeros(T, 2)
        ctrlrange[1] = -1.;  ctrlrange[2] = 1.

        sp = VectorShape(T, nstate)
        op = VectorShape(T, nobs)
        ap = VectorShape(T, naction)
        rp = ScalarShape(T)
        s = zeros(T, nstate)
        o = zeros(T, nobs)
        a = zeros(T, naction)
        env = new{typeof(s),
        typeof(sp), typeof(op), typeof(ap), typeof(rp)}(s,o,a,zeros(T,1),
                                                        zeros(T, nobs),
                                                        ctrlrange,
                                                        sp,op,ap,rp)
        reset!(env)
    end
end

function LyceumBase.tconstruct(::Type{Acrobot}, n::Integer)
    Tuple(Acrobot{Float64}() for _ in 1:n)
end
Acrobot() = first(tconstruct(Acrobot, 1))

@inline LyceumBase.statespace(env::Acrobot) = env.sp
@inline LyceumBase.obsspace(env::Acrobot) = env.op
@inline LyceumBase.actionspace(env::Acrobot) = env.ap
@inline LyceumBase.rewardspace(env::Acrobot) = env.rp

@inline function LyceumBase.getstate!(state, env::Acrobot)
    #state .= env.state
    state[1] = env.state[1] # unroll for speed
    state[2] = env.state[2]
    state[3] = env.state[3]
    state[4] = env.state[4]
end
@inline function LyceumBase.setstate!(env::Acrobot, state)
    #env.state .= state
    env.state[1] = state[1] # unroll for speed
    env.state[2] = state[2]
    env.state[3] = state[3]
    env.state[4] = state[4]
end

@inline LyceumBase.getaction!(action, env::Acrobot) = action .= env.action
@inline LyceumBase.setaction!(env::Acrobot, action) = env.action .= action

@propagate_inbounds function LyceumBase.isdone(state, action, obs, env::Acrobot)
    false
end
@propagate_inbounds function LyceumBase.getobs!(obs, env::Acrobot)
    s = env.state
    o = env.obs
    o[2], o[1] = sincos(s[1])
    o[4], o[3] = sincos(s[2])
    o[5] = s[3]
    o[6] = s[4]

    obs .= env.obs
end
@propagate_inbounds function LyceumBase.getreward(s, a, o, env::Acrobot)
    h1 = o[1]  # Height of first arm
    h2 = o[1] * o[3] - o[2] * o[4]  # Height of second arm
    return - (h1 + h2)

end
@propagate_inbounds function LyceumBase.geteval(s, a, o, env::Acrobot)
    h1 = o[1]  # Height of first arm
    h2 = o[1] * o[3] - o[2] * o[4]  # Height of second arm
    return - (h1 + h2)

end
@propagate_inbounds function LyceumBase.reset!(env::Acrobot)
    randreset!(env)
end
@propagate_inbounds function LyceumBase.randreset!(env::Acrobot)
    env.state .= rand(Uniform(-0.1, 0.1), 4)
    getobs!(env.obs, env)
    env.time[1] = 0.0
    env
end

@inline function _wrappi(v)
    diff = 2pi
    while v >  π; v -= diff; end;
    while v < -π; v += diff; end;
    v
end

@propagate_inbounds function LyceumBase.step!(env::Acrobot)

    state, act, ctrlrange = env.state, env.action, env.ctrlrange
    """
        corresponds to the one in Benchmark paper
        if action[0] < -.33:
            action = 0
        elif action[0] < .33:
            action = 1
        else:
            action = 2
        Also in Gym env
        AVAIL_TORQUE = [-1., 0., +1] , torque = self.AVAIL_TORQUE[a]
    """
    act_d = 0.0
    if act[1] < -0.33; act_d = -1.0; elseif act[1] < 0.33; act_d = 0.0; else; act_d = 1.0; end;

    # BOOK default
    ns = MArray( _rk4(_dsdt, state, act_d, SA_F64[0, 0.2]) )

    #wrap function/ maybe there is more efficient way
    ns[1] = _wrappi(ns[1])
    ns[2] = _wrappi(ns[2])
    ns[3] = clamp(ns[3], -4 * π, 4 * π)
    ns[4] = clamp(ns[4], -9 * π, 9 * π)

    setstate!(env, ns)
    env.time[1] += timestep(env)
    env
end

function _dsdt(s::AbstractArray, a)#, t::AbstractFloat)
        m1 = 1.; m2 = 1.
        l1 = 1.; lc1 = 0.5; lc2 = 0.5
        I1 = 1.; I2 = 1.
        g = 9.8
        theta1, theta2, dtheta1, dtheta2 = s
        sintheta2, costheta2 = sincos(theta2)

        d1       = m1 * lc1^2 + m2 * (l1^2 + lc2^2 + 2 * l1 * lc2 * costheta2) + I1 + I2
        d2       = m2 * (lc2^2 + l1 * lc2 * costheta2) + I2
        phi2     = m2 * lc2 * g * cos(theta1 + theta2 - π / 2.)
        phi1     = - (m2 * l1 * lc2 * dtheta2^2 * sintheta2) - (2 * m2 * l1 * lc2 * dtheta2 * dtheta1 * sintheta2)
                    + ((m1 * lc1 + m2 * l1) * g * cos(theta1 - π / 2) + phi2)
        ddtheta2 = (a + d2 / d1 * phi1 - m2 * l1 * lc2 * dtheta1^2 * sintheta2 - phi2) / (m2 * lc2^2 + I2 - d2^2 / d1)
        ddtheta1 = -(d2 * ddtheta2 + phi1) / d1
        return SVector{4, eltype(s)}(dtheta1, dtheta2, ddtheta1, ddtheta2)
end

function _rk4(derivs::Function, y0::AbstractArray, a, t::AbstractArray)
    #y0 in our case should always be length 4, since we include the action
    yout = SVector{4, eltype(y0)}(y0)
    for i in 1:(length(t)-1)
        thist = t[i]
        dt = t[i + 1] - thist
        dt2 = dt / 2.0
        y0 = yout
        k1 = derivs(y0, a)#, thist)
        k2 = derivs(y0 + dt2 * k1, a)#, thist + dt2)
        k3 = derivs(y0 + dt2 * k2, a)#, thist + dt2)
        k4 = derivs(y0 + dt * k3, a)#, thist + dt)
        yout = y0 + dt / 6.0 * (k1 + (2 * k2) + (2 * k3) + k4)
    end
    return yout # don't store intermediate computations, just return final
end
@inline LyceumBase.timestep(env::Acrobot) = 1.0
@inline Base.time(env::Acrobot) = env.time[1]


text:
@polly function kernel_bicg(A,s,q,p,r)
    n,m = size(A)

    for i = 1:m
        s[i] = zero(eltype(s))
    end

    for i = 1:n
        q[i] = zero(eltype(q));
        for j = 1:m
            s[j] = s[j] + r[i] * A[i,j]
            q[i] = q[i] + A[i,j] * p[j]
        end
    end
end

let
    m = 1900
    n = 2100

    A = zeros(Float32, n, m)
    s = zeros(Float32, m)
    q = zeros(Float32, n)
    p = zeros(Float32, m)
    r = zeros(Float32, n)

    for i = 1:m
        p[i] = (i % m) / m
    end
    for i = 1:n
        r[i] = (i % n) / n
        for j = 1:m
            A[i,j] = (i*(j+1) % n)/n
        end
    end

    SUITE["bicg"] = @benchmarkable kernel_bicg(A, s, q, p, r) setup = (A = copy($A); s = copy($s); q = copy($q); p = copy($p); r = copy($r))
end


text:
using DataGraphs
using Graphs
using MultiDimDictionaries
using NamedGraphs

g = NamedDimGraph(grid((4,)), ["A", "B", "C", "D"])
dg = NamedDimDataGraph{String,Symbol}(g)

@show has_vertex(dg, "A")
@show has_vertex(dg, "D")
@show !has_vertex(dg, 0)
@show !has_vertex(dg, 5)

@show has_edge(dg, "A", "B")
@show has_edge(dg, "A" => "B")
@show !has_edge(dg, "A", "C")
@show !has_edge(dg, "A" => "C")

@show !isassigned(dg, "A")
@show !isassigned(dg, "B")
@show !isassigned(dg, "C")
@show !isassigned(dg, "D")

@show !isassigned(dg, NamedDimEdge("A", "B"))
@show !isassigned(dg, "A" => "B")
@show !isassigned(dg, NamedDimEdge("A" => "B"))
@show !isassigned(dg, "A" => "C")

dg["A"] = "V1"
dg["B"] = "V2"
dg["C"] = "V3"
dg["D"] = "V4"

# Error: does not have vertex
# dg[1] = "X"

@show isassigned(dg, "A")
@show dg["A"] == "V1"
@show dg["B"] == "V2"
@show dg["C"] == "V3"
@show dg["D"] == "V4"

dg["A" => "B"] = :E12
dg["B" => "C"] = :E23
dg[NamedDimEdge("C", "D")] = :E34
@show isassigned(dg, NamedDimEdge("B", "C"))
@show isassigned(dg, "C" => "D")
@show dg[NamedDimEdge("A", "B")] == :E12
@show dg["B" => "C"] == :E23
@show dg["C" => "D"] == :E34


text:
"""
    SimClock(initial_value)

A simulated clock (for simulation or testing purpose)
"""
mutable struct SimClock <: AbstractClock
    value
end

"""
    now(clock::RealTimeClock)

Return current instant from a simulated clock (for simulation or testing purpose)
"""
function now(clock::SimClock)
    clock.value
end

"""
    set(clock::SimClock, new_value)

Set time (with `new_value` of a simulated clock
"""
function set(clock::SimClock, value)
    if value > clock.value
        clock.value = value
    else
        throw(ClockException("time must be ascending"))
    end
end


text:
"""
    attachmap(f,phi)
"""
function attachmap(f,phi)
  k = MapGrad()
  apply_kernel_to_field(k,f,phi)
end

"""
    attachmap(f::AbstractArray,phi::AbstractArray)
"""
function attachmap(f::AbstractArray,phi::AbstractArray)
  k = AddMap()
  apply(k,f,phi)
end

struct MapGrad <: Kernel end

@inline apply_kernel!(cache,k::MapGrad,fx,phix) = fx

function apply_kernel_gradient(k::MapGrad,f,phi)
  g = field_gradient(f)
  jac = field_gradient(phi)
  k = PhysGrad()
  apply_kernel_to_field(k,g,jac)
end

struct PhysGrad <: Kernel end

function kernel_cache(k::PhysGrad,a,b)
  _attachmap_checks(a,b)
  Ta = eltype(a)
  Tb = eltype(b)
  T = return_type(*,return_type(inv,Tb),Ta)
  r = zeros(T,size(a))
  CachedArray(r)
end

function _attachmap_checks(a,b)
  @assert ndims(a) == 2 "attachmap: map can only be attached to bases of fields."
  ia,ja = size(a)
  ib = length(b)
  @assert ia == ib "attachmap: basis and jacobian size mismatch."
end

@inline function apply_kernel!(r,k::PhysGrad,a,b)
  _attachmap_checks(a,b)
  s = size(a)
  setsize!(r,s)
  np, ni = s
  for p in 1:np
    @inbounds jacinv = inv(b[p])
    for i in 1:ni
      @inbounds r[p,i] = jacinv * a[p,i]
    end
  end
  r
end

struct AddMap <: Kernel end

@inline function apply_kernel!(cache,k::AddMap,fi,phii)
  attachmap(fi,phii)
end

function apply_gradient(k::AddMap,f,phi)
  g = gradient(f)
  jac = gradient(phi)
  k = PhysGrad()
  apply_to_field_array(k,g,jac)
end

function kernel_evaluate(k::AddMap,x,f,phi)
  evaluate(f,x)
end


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/Materials.jl/blob/master/LICENSE

using Test, Tensors

let nu = 0.3,
    yield_strength=100.0,
    parameters = PerfectPlasticParameterState(youngs_modulus=200.0e3,
                                              poissons_ratio=nu,
                                              yield_stress=yield_strength),
    epsilon=1e-3,
    mat,  # scope the name to this level; actual definition follows later
    tostrain(vec) = fromvoigt(Symm2, vec; offdiagscale=2.0),
    tostress(vec) = fromvoigt(Symm2, vec),
    uniaxial_stress(sigma) = tostress([sigma, 0, 0, 0, 0, 0])
    let dtime=0.25
        # elastic straining
        dstrain_dtime = tostrain(epsilon*[1.0, -nu, -nu, 0.0, 0.0, 0.0])
        ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_dtime*dtime)
        mat = PerfectPlastic(parameters=parameters, ddrivers=ddrivers)
        integrate_material!(mat)
        update_material!(mat)
        @test isapprox(mat.variables.stress, uniaxial_stress(yield_strength / 2))

        mat.ddrivers = ddrivers
        integrate_material!(mat)
        update_material!(mat)
        # We should now be at the yield surface.
        @test isapprox(mat.variables.stress, uniaxial_stress(yield_strength))
        @test isapprox(mat.variables.cumeq, 0.0; atol=1.0e-12)

        # plastic straining
        # von Mises material, plastically incompressible, so plastic nu=0.5.
        dstrain_dtime = tostrain(epsilon*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0])
        ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_dtime*dtime)
        mat.ddrivers = ddrivers
        integrate_material!(mat)
        update_material!(mat)
        @test isapprox(mat.variables.stress, uniaxial_stress(yield_strength); atol=1.0e-12)
        @test isapprox(mat.variables.cumeq, dtime*epsilon)

        # return to elastic state
        dstrain_dtime = tostrain(-epsilon*[1.0, -nu, -nu, 0.0, 0.0, 0.0])
        ddrivers = PerfectPlasticDriverState(time=dtime, strain=dstrain_dtime*dtime)
        mat.ddrivers = ddrivers
        integrate_material!(mat)
        update_material!(mat)
        @test isapprox(mat.variables.stress, uniaxial_stress(yield_strength / 2); atol=1.0e-12)
    end

    let dtime=1.0
        # loading in reverse direction to plastic state
        # The 0.75 term: one 0.25 cancels the current elastic stress state,
        # and the remaining 0.5 reaches the yield surface.
        # The 0.25 term: plastic strain.
        dstrain_dtime = (-0.75*tostrain(epsilon*[1.0, -nu, -nu, 0.0, 0.0, 0.0])
                         -0.25*tostrain(epsilon*[1.0, -0.5, -0.5, 0.0, 0.0, 0.0]))
        ddrivers = PerfectPlasticDriverState(time=1.0, strain=dstrain_dtime*dtime)
        mat.ddrivers = ddrivers
        integrate_material!(mat)
        integrate_material!(mat)
        update_material!(mat)
        @test isapprox(mat.variables.stress, uniaxial_stress(-yield_strength))
    end
end


text:
"""
   AbstractBlockBandedLayout

isa a `MemoryLayout` that indicates that the array implements the block-banded
interface.
"""
abstract type AbstractBlockBandedLayout <: AbstractBlockLayout end

"""
   AbstractBandedBlockBandedLayout

isa a `MemoryLayout` that indicates that the array implements the banded-block-banded
interface.
"""
abstract type AbstractBandedBlockBandedLayout <: AbstractBlockBandedLayout end

struct BandedBlockBandedLayout <: AbstractBandedBlockBandedLayout end
struct BlockBandedLayout <: AbstractBlockBandedLayout end

struct BandedBlockBandedColumns{LAY} <: AbstractBandedBlockBandedLayout end
struct BandedBlockBandedRows{LAY} <: AbstractBandedBlockBandedLayout end
struct BlockBandedColumns{LAY} <: AbstractBlockBandedLayout end
struct BlockBandedRows{LAY} <: AbstractBlockBandedLayout end

const BandedBlockBandedColumnMajor = BandedBlockBandedColumns{ColumnMajor}
const BandedBlockBandedRowMajor = BandedBlockBandedRows{ColumnMajor}
const BlockBandedColumnMajor = BlockBandedColumns{ColumnMajor}
const BlockBandedRowMajor = BlockBandedRows{ColumnMajor}

transposelayout(::BandedBlockBandedColumns{Lay}) where Lay = BandedBlockBandedRows{Lay}()
transposelayout(::BandedBlockBandedRows{Lay}) where Lay = BandedBlockBandedColumns{Lay}()
transposelayout(::BlockBandedColumns{Lay}) where Lay = BlockBandedRows{Lay}()
transposelayout(::BlockBandedRows{Lay}) where Lay = BlockBandedColumns{Lay}()

conjlayout(::Type{T}, ::BandedBlockBandedColumns{Lay}) where {T<:Complex,Lay} = BandedBlockBandedColumns{typeof(conjlayout(T,Lay))}()
conjlayout(::Type{T}, ::BandedBlockBandedRows{Lay}) where {T<:Complex,Lay} = BandedBlockBandedRows{typeof(conjlayout(T,Lay))}()
conjlayout(::Type{T}, ::BlockBandedColumns{Lay}) where {T<:Complex,Lay} = BlockBandedColumns{typeof(conjlayout(T,Lay))}()
conjlayout(::Type{T}, ::BlockBandedRows{Lay}) where {T<:Complex,Lay} = BlockBandedRows{typeof(conjlayout(T,Lay))}()

# AbstractBandedMatrix must implement

# A BlockBandedMatrix is a BlockMatrix, but is not a BandedMatrix
abstract type AbstractBlockBandedMatrix{T} <: AbstractBlockMatrix{T} end
MemoryLayout(::Type{<:AbstractBlockBandedMatrix}) = BlockBandedLayout()

"""
    blockbandwidths(A)

Returns a tuple containing the upper and lower blockbandwidth of `A`.
"""
blockbandwidths(A::AbstractVecOrMat) = blockbandwidths(MemoryLayout(A), axes(A), A)
blockbandwidths(_, _, A) = (blocksize(A,1)-1 , blocksize(A,2)-1)

"""
    blockbandwidth(A,i)

Returns the lower blockbandwidth (`i==1`) or the upper blockbandwidth (`i==2`).
"""
blockbandwidth(A, k::Integer) = blockbandwidths(A)[k]

"""
    bandrange(A)

Returns the range `-blockbandwidth(A,1):blockbandwidth(A,2)`.
"""
blockbandrange(A::AbstractMatrix) = -blockbandwidth(A,1):blockbandwidth(A,2)

# start/stop indices of the i-th column/row, bounded by actual matrix size
@inline blockbanded_blockcolstart(A, i::Integer) = Block(max(i-colblockbandwidth(A,2)[i], 1))
@inline blockbanded_blockcolstop(A, i::Integer) = Block(max(min(i+colblockbandwidth(A,1)[i], blocksize(A, 1)), 0))
@inline blockbanded_blockrowstart(A, i::Integer) = Block(max(i-rowblockbandwidth(A,1)[i], 1))
@inline blockbanded_blockrowstop(A, i::Integer) = Block(max(min(i+rowblockbandwidth(A,2)[i], blocksize(A, 2)), 0))
for Func in (:blockbanded_blockcolstart, :blockbanded_blockcolstop, :blockbanded_blockrowstart, :blockbanded_blockrowstop)
    @eval $Func(A, i::Block{1}) = $Func(A, Int(i))
end

blockbanded_blockcolstart(A, i::BlockRange) = blockbanded_blockcolstart(A, minimum(i))
blockbanded_blockrowstart(A, i::BlockRange) = blockbanded_blockrowstart(A, minimum(i))
blockbanded_blockcolstop(A, i::BlockRange) = blockbanded_blockcolstop(A, maximum(i))
blockbanded_blockrowstop(A, i::BlockRange) = blockbanded_blockrowstop(A, maximum(i))

const AllBlockBandedLayout{UPLO,UNIT} = Union{AbstractBlockBandedLayout,TriangularLayout{UPLO,UNIT,<:AbstractBlockBandedLayout}}

@inline blockcolsupport(::AllBlockBandedLayout, A, i) = isempty(i) ? (Block(1):Block(0)) : blockbanded_blockcolstart(A,i):blockbanded_blockcolstop(A,i)
@inline blockrowsupport(::AllBlockBandedLayout, A, i) = isempty(i) ? (Block(1):Block(0)) : blockbanded_blockrowstart(A,i):blockbanded_blockrowstop(A,i)

# length of i-the column/row
@inline blockcollength(A, i) = max(Int(blockcolstop(A, i)) - Int(blockcolstart(A, i)) + 1, 0)
@inline blockrowlength(A, i) = max(Int(blockrowstop(A, i)) - Int(blockrowstart(A, i)) + 1, 0)

# this gives the block bandwidth in each block column/row
@inline colblockbandwidths(A::AbstractMatrix) = Fill.(blockbandwidths(A), blocksize(A,2))
@inline rowblockbandwidths(A::AbstractMatrix) = Fill.(blockbandwidths(A), blocksize(A,1))

@inline colblockbandwidth(bs, i::Int) = colblockbandwidths(bs)[i]
@inline rowblockbandwidth(bs, i::Int) = rowblockbandwidths(bs)[i]

"""
    isblockbanded(A)

returns true if a matrix implements the block banded interface.
"""
isblockbanded(::AbstractBlockBandedMatrix) = true
isblockbanded(_) = false

# override bandwidth(A,k) for each AbstractBandedMatrix
# override inbands_getindex(A,k,j)

# return id of first empty diagonal intersected along row k
function _firstblockdiagrow(A::AbstractMatrix, k::Int)
    a, b = blockrowstart(A, k), blockrowstop(A, k)
    c = a == 1 ? b+1 : a-1
    c-k
end

# return id of first empty diagonal intersected along column j
function _firstblockdiagcol(A::AbstractMatrix, j::Int)
    a, b = blockcolstart(A, j), blockcolstop(A, j)
    r = a == 1 ? b+1 : a-1
    j-r
end

######################################
#  RaggedMatrix interface
######################################

@inline function blockbanded_colstart(A, i::Integer)
    bs = blockcolstart(A,findblock(axes(A,2),i))
    if isempty(axes(A,1))
        1
    elseif Int(bs) ≤ blocksize(A, 2)
        first(axes(A,1)[bs])
    else
        size(A,1)+1
    end
end

@inline function blockbanded_colstop(A, i::Integer)
    CS = blockcolstop(A,findblock(axes(A,2),i))
    CS == Block(0) && return 0
    last(axes(A,1)[CS])
end

@inline blockbanded_rowstart(A, i::Integer) =
    first(axes(A,2)[blockrowstart(A,findblock(axes(A,1),i))])

@inline function blockbanded_rowstop(A, i::Integer)
    CS = blockrowstop(A,findblock(axes(A,1),i))
    CS == Block(0) && return 0
    last(axes(A,2)[CS])
end

@inline blockbanded_colsupport(A, j::Integer) = blockbanded_colstart(A, j):blockbanded_colstop(A, j)
@inline blockbanded_rowsupport(A, j::Integer) = blockbanded_rowstart(A, j):blockbanded_rowstop(A, j)

@inline blockbanded_rowsupport(A, j) = isempty(j) ? (1:0) : blockbanded_rowstart(A,minimum(j)):blockbanded_rowstop(A,maximum(j))
@inline blockbanded_colsupport(A, j) = isempty(j) ? (1:0) : blockbanded_colstart(A,minimum(j)):blockbanded_colstop(A,maximum(j))

@inline colsupport(::AbstractBlockBandedLayout, A, j) = blockbanded_colsupport(A, j)
@inline rowsupport(::AbstractBlockBandedLayout, A, j) = blockbanded_rowsupport(A, j)

# default implementation loops over all indices, including zeros
function fill!(A::AbstractBlockBandedMatrix, val::Any)
  iszero(val) || throw(BandError(A))
  fill!(A.data, val)
  A
end

function ==(A::AbstractBlockBandedMatrix, B::AbstractBlockBandedMatrix)
    axes(A) == axes(B) || return false
    blockisequal(axes(A), axes(B)) || return Base.invoke(==, Tuple{AbstractArray,AbstractArray}, A, B)
    l,u = max.(blockbandwidths(A), blockbandwidths(B))
    N = blocksize(A,1)
    for J = blockaxes(A,2), K = max(Block(1),J-u):min(J+l,Block(N))
        view(A, K, J) == view(B, K, J) || return false
    end
    return true
end