text:
module Interfaces

export isinterfacetype, @interface, @implements, InterfaceImplementationError

struct InterfaceImplementationError <: Exception
    msg::String
end

isinterfacetype(::Type{T}) where {T} = false
isinterfacetype(::Type{Type{T}}) where {T} = isinterfacetype(T)

function interface end

struct Interface
    name
    exprs::Vector{Any}
end

function Base.show(io::IO, ::MIME"text/plain", x::Interface)
    print(io, "Interface: $(x.name)")
    foreach(ex -> showexpr(io, ex), x.exprs)
    return
end

unwrapblock(ex::Expr) = ex.head == :block ? ex.args[1] : ex

function showexpr(io::IO, arg, indent=1)
    if arg.head == :call
        print(io, "\n" * "  "^indent * "* method definition: `$arg`")
    elseif arg.head == :(::)
        print(io, "\n" * "  "^indent * "* method definition with required return type: `$arg`")
    elseif arg.head == :<:
        print(io, "\n" * "  "^indent * "* subtyping: `T <: $(arg.args[2])`")
    elseif arg.head == :if
        print(io, "\n" * "  "^indent * "* conditional requirements:")
        print(io, "\n" * "  "^(indent + 1) * "* if $(unwrapblock(arg.args[1]))")
        showexpr(io, arg.args[2], indent + 2)
        while length(arg.args) > 2
            if arg.args[3].head == :elseif
                arg = arg.args[3]
                print(io, "\n" * "  "^(indent + 1) * "* elseif $(unwrapblock(arg.args[1]))")
                showexpr(io, arg.args[2], indent + 2)
            else
                # else block
                print(io, "\n" * "  "^(indent + 1) * "* else")
                showexpr(io, arg.args[3], indent + 2)
                break
            end
        end
    elseif arg.head == :||
        print(io, "\n" * "  "^indent * "* one of the following required:")
        while true
            showexpr(io, arg.args[1], indent + 1)
            arg.args[2].head == :|| || break
            arg = arg.args[2]
        end
        showexpr(io, arg.args[2], indent + 1)
    elseif arg.head == :block
        for x in arg.args
            showexpr(io, x, indent)
        end
    elseif arg.head == :macrocall && arg.args[1] == Symbol("@optional")
        print(io, "\n" * "  "^indent * "* optional interface requirement:")
        showexpr(io, arg.args[3], indent + 1)
    end
    return
end

function implements end

implements(::Type{Type{T}}, ::Type{Type{IT}}, mods::Vector{Module}=[parentmodule(T)]; debug::Bool=false) where {T, IT} =
    implements(T, IT, mods; debug=debug)

function implemented(f, args, mods)
    impls = Base.methods(f, args, mods)
    return length(impls) > 0
end

"""
Replace all arguments `T` with `sym` (default `:T`) in the given `expr::Expr`.

Used to swap `::SomeInterface` -> `::T`, because `@interface` uses syntax like:

    @interface SomeInterface begin
        foo(::SomeInterface)
    end

to express "`T` implements the interface `SomeInterface` if there is a method `foo(::T)`",
i.e. we need there to be a `foo` method that accepts the type `T` not the type `SomeInterface`.
"""
function recursiveswapsymbols!(IT, expr, sym=:T)
    for i = 1:length(expr.args)
        arg = expr.args[i]
        if arg == IT
            expr.args[i] = sym
        elseif arg isa Expr
            recursiveswapsymbols!(IT, arg, sym)
        end
    end
    return
end

"""
Convert arguments to just their type constraints.
Used as part of converting method signature syntax to a something we can check for with
`methods(f, args)`.

    # Given `foo(x::Int, ::Float64, z)`
    args = [:(1::Int), :(::Float64), :z]
    convertargs(:SomeInterface, :foo, args) == :(Tuple{Int, Float64, Any})
"""
function convertargs(IT, nm, args)
    isempty(args) && throw(ArgumentError("invalid `$IT` interface method with zero arguments: `$nm()`"))
    for i = 1:length(args)
        arg = args[i]
        if arg isa Symbol
            args[i] = :Any
        elseif arg.head === :(::) && length(arg.args) == 1
            recursiveswapsymbols!(IT, arg)
            args[i] = arg.args[1]
        elseif arg.head === :(::) && length(arg.args) == 2
            recursiveswapsymbols!(IT, arg)
            args[i] = arg.args[2]
        else
            throw(ArgumentError("invalid `$IT` interface method argument for method `$nm`: `$arg`"))
        end
    end
    return Expr(:curly, :Tuple, args...)
end

unconvertargs(::Type{T}) where {T <: Tuple} = Any[Expr(:(::), fieldtype(T, i)) for i = 1:fieldcount(T)]

"""
Extract the (quoted) function name and argument types from a method signature.

    methodparts(:SomeInterface, :(foo(x::Int, ::Float64, z))) == (:foo, :(Tuple{Int, Float64, Any}))
"""
function methodparts(IT, x::Expr)
    @assert x.head === :call
    methodname = x.args[1]
    args = convertargs(IT, methodname, x.args[2:end])
    return methodname, args
end

"""
Construct the method signature from a function name and argument types.

    methodsig(:foo, Tuple{Int64, Float64, Any}) == :(foo(::Int64, ::Float64, ::Any))
"""
methodsig(f, args) = Expr(:call, f, unconvertargs(args)...)

function requiredmethod(IT, nm, args, shouldthrow)
    return quote
        check = Interfaces.implemented($nm, $args, mods)
        check || ($shouldthrow && Interfaces.missingmethod(loglevel, $IT, $nm, $args, mods))
        check
    end
end

function requiredreturn(IT, nm, args, shouldthrow, RT_sym, __RT__)
    return quote
        check = $(requiredmethod(IT, nm, args, shouldthrow))
        if check
            $RT_sym = Interfaces.returntype($nm, $args)
            # @show $RT_sym, $nm, $args, Interfaces.isinterfacetype($__RT__)
            check &= Interfaces.isinterfacetype($__RT__) ?  Interfaces.implements($RT_sym, $__RT__) : $RT_sym <: $__RT__
            check || ($shouldthrow && Interfaces.invalidreturntype(loglevel, $IT, $nm, $args, $RT_sym, $__RT__))
        end
        check
    end
end

@noinline function missingmethod(loglevel, IT, f, args, mods)
    log_or_throw(loglevel, "missing `$IT` interface method definition: `$(Interfaces.methodsig(f, args))`, in module(s): `$mods`")
end

@noinline function invalidreturntype(loglevel, IT, f, args, RT1, RT2)
    log_or_throw(loglevel, "invalid return type for `$IT` interface method definition: `$(methodsig(f, args))`; inferred $RT1, required $RT2")
end

@noinline function subtypingrequired(loglevel, IT, T)
    log_or_throw(loglevel, "interface `$IT` requires implementing types to subtype, like: `struct $T <: $IT`")
end

@noinline function atleastonerequired(loglevel, IT, expr)
    log_or_throw(loglevel, "for `$IT` interface, one of the following method definitions is required: `$expr`")
end

function log_or_throw(loglevel, msg)
    if loglevel === :required
        log_required(msg)
    elseif loglevel === :optional
        log_optional(msg)
    else
        throw(InterfaceImplementationError(msg))
    end
end

function log_required(msg)
    printstyled("[ InterfaceImplementationError: "; color=Base.error_color())
    println(msg)
end

function log_optional(msg)
    printstyled("[ @optional definition missing: "; color=Base.warn_color())
    println(msg)
end

# `shouldthrow` indicates whether or not it is an error for a particular expression not
# to be satisfied. We set this to false for the individual expressions in `ex1 || ex2`.
function toimplements!(IT, arg::Expr, shouldthrow::Bool=true)
    if arg.head == :call
        # required method definition
        nm, args = methodparts(IT, arg)
        return requiredmethod(IT, nm, args, shouldthrow)
    elseif arg.head == :(::)
        # required method definition and required return type
        nm, args = methodparts(IT, arg.args[1])
        annotation = arg.args[2]
        if !isa(annotation, Symbol) && annotation.head == :where
            # `::(T where T)` or `::(T where T<:Foo)`
            sym = annotation.args[1]
            __RT__ = annotation.args[2] isa Symbol ? Any : annotation.args[2].args[2]
        else # `::Foo` or `::Union{Foo,Bar}` or `::Type{Foo}`
            sym = gensym()
            __RT__ = annotation
        end
        return requiredreturn(IT, nm, args, shouldthrow, sym, __RT__)
    elseif arg.head == :<:
        RT = arg.args[2]
        if RT == IT
            return quote
                check = T <: $RT
                check || Interfaces.subtypingrequired(loglevel, $RT, T)
                check
            end
        else
            return quote
                check = Interfaces.isinterfacetype($RT) ? Interfaces.implements(T, $RT, mods) : T <: $RT
                check || Interfaces.subtypingrequired(loglevel, $RT, T)
                check
            end
        end
    elseif arg.head == :if
        # conditional requirement
        origarg = arg
        recursiveswapsymbols!(IT, arg.args[1])
        arg.args[2] = toimplements!(IT, arg.args[2])
        while length(arg.args) > 2
            if arg.args[3].head == :elseif
                arg = arg.args[3]
                recursiveswapsymbols!(IT, arg.args[1])
                arg.args[2] = toimplements!(IT, arg.args[2])
            else
                # else block
                arg.args[3] = toimplements!(IT, arg.args[3])
                break
            end
        end
        return origarg
    elseif arg.head == :||
        # one of many required
        argcopy = copy(arg)
        origarg = arg
        while true
            arg.args[1] = toimplements!(IT, arg.args[1], false)
            arg.args[2].head == :|| || break
            arg = arg.args[2]
        end
        arg.args[2] = toimplements!(IT, arg.args[2], false)
        return quote
            check = $origarg
            check || Interfaces.atleastonerequired(loglevel, $IT, $(Meta.quot(argcopy)))
            check
        end
    elseif arg.head == :block
        # not expected at top-level of @interface block
        # but can be block of if-else or || expressions
        map!(x -> toimplements!(IT, x, shouldthrow), arg.args, arg.args)
        return arg
    elseif arg.head == :macrocall && arg.args[1] == Symbol("@optional")
        return quote
            if debug
                oldlvl = loglevel
                loglevel = :optional
                check = $(toimplements!(IT, arg.args[3], shouldthrow))
                loglevel = oldlvl
                check
            else
                true
            end
        end
    else
        io = IOBuffer()
        dump(io, arg)
        ex = String(take!(io))
        throw(ArgumentError("unsupported expression in @interface block for `$IT`: `$ex`"))
    end
end

macro interface(IT, alias_or_block, maybe_block=nothing)
    @assert IT isa Symbol || IT.head == :.
    if alias_or_block isa Symbol
        alias = alias_or_block
        block = maybe_block
    else
        alias = IT
        block = alias_or_block
    end
    @assert block isa Expr && block.head == :block
    Base.remove_linenums!(block)
    if IT !== alias
        # Swap the alias for the InterfaceType itself, so later we can generically swap the
        # InterfaceType for an arbitrary symbol e.g. :T, regardless of if an alias was used.
        recursiveswapsymbols!(alias, block, IT)
    end
    iface = Interface(IT, deepcopy(block.args))
    filter!(x -> !(x isa String), block.args)
    toimplements!(IT, block)
    return esc(quote
        Interfaces.isinterfacetype(::Type{$IT}) = true
        Interfaces.interface(::Type{$IT}) = $iface
        function Interfaces.implements(::Type{T}, ::Type{$IT}, mods::Vector{Module}=[parentmodule(T)]; debug::Bool=false) where {T}
            loglevel = debug ? :required : :none
            $block
        end
    end)
end

macro implements(T, IT, debug=false)
    x = debug isa Bool ? debug :
        (debug isa Expr && debug.head == :(=)) ? debug.args[2] :
        throw(ArgumentError("unsupported argument in @implements: `$debug`"))
    return esc(quote
        @assert Interfaces.implements($T, $IT; debug=$x)
        Interfaces.implements(::Type{$T}, ::Type{$IT}) = true
    end)
end

@noinline returntype(@nospecialize(f), @nospecialize(args)) = Base.return_types(f, args)[1]

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

struct PrivateEndpointConnectionsApi <: SwaggerApi
    client::Swagger.Client
end

"""

Deletes the specified private endpoint connection associated with the storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: privateEndpointConnectionName::String (required)
Return: Nothing
"""
function _swaggerinternal_privateEndpointConnectionsDelete(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsDelete", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsDelete", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "privateEndpointConnectionsDelete", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "privateEndpointConnectionsDelete", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "privateEndpointConnectionsDelete", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "privateEndpointConnectionsDelete", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "DELETE", Nothing, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/privateEndpointConnections/{privateEndpointConnectionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "privateEndpointConnectionName", privateEndpointConnectionName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function privateEndpointConnectionsDelete(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsDelete(_api, resourceGroupName, accountName, api_version, subscriptionId, privateEndpointConnectionName; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function privateEndpointConnectionsDelete(_api::PrivateEndpointConnectionsApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsDelete(_api, resourceGroupName, accountName, api_version, subscriptionId, privateEndpointConnectionName; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Gets the specified private endpoint connection associated with the storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: privateEndpointConnectionName::String (required)
Return: PrivateEndpointConnection
"""
function _swaggerinternal_privateEndpointConnectionsGet(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsGet", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsGet", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "privateEndpointConnectionsGet", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "privateEndpointConnectionsGet", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "privateEndpointConnectionsGet", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "privateEndpointConnectionsGet", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", PrivateEndpointConnection, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/privateEndpointConnections/{privateEndpointConnectionName}", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "privateEndpointConnectionName", privateEndpointConnectionName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function privateEndpointConnectionsGet(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsGet(_api, resourceGroupName, accountName, api_version, subscriptionId, privateEndpointConnectionName; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function privateEndpointConnectionsGet(_api::PrivateEndpointConnectionsApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsGet(_api, resourceGroupName, accountName, api_version, subscriptionId, privateEndpointConnectionName; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

List all the private endpoint connections associated with the storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Return: PrivateEndpointConnectionListResult
"""
function _swaggerinternal_privateEndpointConnectionsList(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsList", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsList", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "privateEndpointConnectionsList", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "privateEndpointConnectionsList", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "privateEndpointConnectionsList", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "privateEndpointConnectionsList", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "GET", PrivateEndpointConnectionListResult, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/privateEndpointConnections", ["azure_auth"])
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function privateEndpointConnectionsList(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsList(_api, resourceGroupName, accountName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function privateEndpointConnectionsList(_api::PrivateEndpointConnectionsApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsList(_api, resourceGroupName, accountName, api_version, subscriptionId; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

"""

Update the state of specified private endpoint connection associated with the storage account.
Param: resourceGroupName::String (required)
Param: accountName::String (required)
Param: api_version::String (required)
Param: subscriptionId::String (required)
Param: privateEndpointConnectionName::String (required)
Param: properties::PrivateEndpointConnection (required)
Return: PrivateEndpointConnection
"""
function _swaggerinternal_privateEndpointConnectionsPut(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String, properties; _mediaType=nothing)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsPut", :maxLength, resourceGroupName, 90)
    Swagger.validate_param("resourceGroupName", "privateEndpointConnectionsPut", :minLength, resourceGroupName, 1)

    Swagger.validate_param("accountName", "privateEndpointConnectionsPut", :maxLength, accountName, 24)
    Swagger.validate_param("accountName", "privateEndpointConnectionsPut", :minLength, accountName, 3)

    Swagger.validate_param("api_version", "privateEndpointConnectionsPut", :minLength, api_version, 1)

    Swagger.validate_param("subscriptionId", "privateEndpointConnectionsPut", :minLength, subscriptionId, 1)

    _ctx = Swagger.Ctx(_api.client, "PUT", PrivateEndpointConnection, "/subscriptions/{subscriptionId}/resourceGroups/{resourceGroupName}/providers/Microsoft.Storage/storageAccounts/{accountName}/privateEndpointConnections/{privateEndpointConnectionName}", ["azure_auth"], properties)
    Swagger.set_param(_ctx.path, "resourceGroupName", resourceGroupName)  # type String
    Swagger.set_param(_ctx.path, "accountName", accountName)  # type String
    Swagger.set_param(_ctx.path, "subscriptionId", subscriptionId)  # type String
    Swagger.set_param(_ctx.path, "privateEndpointConnectionName", privateEndpointConnectionName)  # type String
    Swagger.set_param(_ctx.query, "api-version", api_version)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/json"] : [_mediaType])
    return _ctx
end

function privateEndpointConnectionsPut(_api::PrivateEndpointConnectionsApi, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String, properties; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsPut(_api, resourceGroupName, accountName, api_version, subscriptionId, privateEndpointConnectionName, properties; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function privateEndpointConnectionsPut(_api::PrivateEndpointConnectionsApi, response_stream::Channel, resourceGroupName::String, accountName::String, api_version::String, subscriptionId::String, privateEndpointConnectionName::String, properties; _mediaType=nothing)
    _ctx = _swaggerinternal_privateEndpointConnectionsPut(_api, resourceGroupName, accountName, api_version, subscriptionId, privateEndpointConnectionName, properties; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export privateEndpointConnectionsDelete, privateEndpointConnectionsGet, privateEndpointConnectionsList, privateEndpointConnectionsPut

ext: jl


text:
using Mass2

srand(hash("mockpluses.jl"))
pg = TupacLikeTwoExponentialPulseGenerator(Int, Normal(1000,0.5))
getcleanpulse(pg,1000)
gettriggeredpulse!(pg,10)

function compute_whitenoise_filter(pulse, selection_good)
	filter = mean(pulse[selection_good])
	filter -= mean(filter) # 0 average
	normalization = (maximum(filter)-minimum(filter))./dot(filter, filter)
	filter*normalization
end
function filter1lag(pulse, filter)
	out = Array(Float64, length(pulse))
	for i=1:length(pulse)
		out[i] = dot(pulse[i], filter)
	end
	out
end
cleanpulses = [getcleanpulse(pg,1000) for i = 1:10000]
filter = compute_whitenoise_filter(cleanpulses, trues(length(cleanpulses)))
cleanpulses_filt_values = filter1lag(cleanpulses, filter)
@test abs(mean(cleanpulses_filt_values)-1000)<1
@test abs(std(cleanpulses_filt_values)-.655)<.1


text:
export JSONRPCEndpoint, send_notification, send_request, send_success_response, send_error_response

include("core.jl")
include("typed.jl")
include("interface_def.jl")


text:
using Documenter
using jInv.Mesh
using jInv.Utils
using jInv.ForwardShare
using jInv.LinearSolvers
using jInv.InverseSolve

makedocs()

deploydocs(
    deps   = Deps.pip("mkdocs", "python-markdown-math"),
    repo = "github.com/JuliaInv/jInv.jl.git",
    julia  = "0.5"
)


text:
cd("INSERT_PATH")

using Distributions
# using Polynomials
using PolynomialRoots
using DataFrames

include("toolBox.jl")
include("polyDNRR.jl")

#########################################################
# Parameter Declaration
#########################################################

ss = 270

## used data sets
# w1, lam1 = 0.6, 1.8e04 # f_{2,1}, (ss200), dataseeds: 23,24,21,32 , etas = 3.5,4.9,6.3,7.5 (sampseeds=94,95,93,96)
# w1, lam1 = 0.6, 3e04 # f_{2,1}, (ss200), dataseed: 44 , eta = 2.7 (sampseed=99)

w1, lam1 = 0.6, 1e05 # f_{2,1}, (ss200), dataseeds: 23,24,21,32 , etas = 3.5,4.9,6.3,7.5 (sampseeds=94,95,93,96)

lam2 = 1e-02 * lam1

# map parameters and initial condition
theta, x0 = [0.05, 2.55, 0., -0.99], 1.
# theta, x0 = [1., 0., -1.71], 0.1

# corresponding deterministic orbit
xdet = zeros(ss)
xdet[1] = polyMap(theta, x0)
for i in 2:ss
  xdet[i] = polyMap(theta, xdet[i - 1])
end

# orbit contaminated with dynamical noise
# first argument is  type of noise: 0 is parametric, 1-3 nonparametric
dataSeed = 2
data = genData(3, ss, theta, x0, w1, lam1, lam2, dataSeed)
eta = sqrt(w1*1/lam1+(1-w1)*1/lam2)/std(data)

#########################################################
# GSB Reconstruction - Denoising
#########################################################

# Dynamical Noise Only

degree, k, gibbsIter, burnIn, dela, delb, wdist, gEps1, gEps2, thLow, thUp, zLow, zUp, papr, pbpr, thin, samplerSeed =
  5, ss, 150000, 50000, 1e05, 1e-03, 1e-03, 1e-03, 1e-03, -10.0, 10.0, -10.0, 10.0, 1., 1., 10, 999;

thetas, x0s, noise, clusters, ps, deltas, ys = 0, 0, 0, 0, 0, 0, 0;
gc()

filename = "/Results"
savelocation = string(pwd(), filename, "/seed$samplerSeed")
mkpath(savelocation)
writedlm(string(savelocation,"/xdet.txt"), xdet)

@time thetas = polyDNRR(data, degree, k, gibbsIter, burnIn, dela, delb, wdist, gEps1, gEps2, thLow, thUp, zLow, zUp, papr, pbpr, thin, samplerSeed);

thetas = 0;
gc()

max_stars_repo_name: KosKal/DNRR

max_issues_repo_name: KosKal/DNRR

max_forks_repo_name: KosKal/DNRR


text:
module TestAqua

import Aqua
import ConcurrentCollections

test() = Aqua.test_all(ConcurrentCollections; unbound_args = false)

end  # module


text:
#
# ```@setup signalsim
# # the example command line input for this script
# ARGS = ["example.toml",  "143", "143", "TT", "10", "--plot"]
# ```
#
# # [Signal Sims (signalsim.jl)](@id signalsim)
# This script runs signal-only simulations, and then computes spectra from those simulations
# after applying and then correcting for the mode-coupling of the mask. The point-source holes
# in the Planck likelihood masks are insufficiently apodized, which violates an assumption
# in the estimation of analytic covariance matrices. We wish to correct for these, so we
# generate a bunch of simulations and compare the analytic and sample covariance.
#
# ```@raw html
# <pre class="shell">
# <code class="language-shell hljs">$ julia signalsim.jl example.toml 143 143 TT 10</code></pre>
# ```
#
# ## Configuration
# We set up various information we'll need later.
#
configfile, freq1, freq2, spec, nsims = ARGS

using TOML
using Plots
using Healpix
using JLD2, UUIDs  # for saving sim arrays
include("util.jl")

config = TOML.parsefile(configfile)
nside = config["general"]["nside"]
run_name = config["general"]["name"]
spectrapath = joinpath(config["scratch"], "rawspectra")
lmax = nside2lmax(nside)
lmax_planck = min(2508, lmax)
splits = "1", "2"  # planck never uses any other splits

# Next, we prepare the input signal used to generate these simulations.

signal11, th = signal_and_theory(freq1, freq1, config)
signal12, th = signal_and_theory(freq1, freq2, config)
signal22, th = signal_and_theory(freq2, freq2, config)
signal = Dict(("1", "1") => signal11, ("1", "2") => signal12, ("2", "2") => signal22)

# We'll store it in a nice array for access later, when simulating.

𝐂 = zeros(2, 2, lmax+1)
X, Y = spec
inds = 1:(lmax_planck+1)
𝐂[1,1,inds] .= signal[splits[1], splits[1]][X * X][inds]
𝐂[1,2,inds] .= signal[splits[1], splits[2]][X * Y][inds]
𝐂[2,1,inds] .= signal[splits[1], splits[2]][X * Y][inds]
𝐂[2,2,inds] .= signal[splits[2], splits[2]][Y * Y][inds];

# Next, we pre-allocate some structures that we'll re-use between simulation iterations.

m1 = PolarizedHealpixMap{Float64, RingOrder}(nside)
m2 = PolarizedHealpixMap{Float64, RingOrder}(nside)
a1 = [Alm(lmax, lmax) for i in 1:3]
a2 = [Alm(lmax, lmax) for i in 1:3]

X, Y = Symbol(spec[1]), Symbol(spec[2])
run_name = config["general"]["name"]
masktype1 = (X == :T) ? "T" : "P"
masktype2 = (Y == :T) ? "T" : "P"
mapid1 = "P$(freq1)hm$(splits[1])"
mapid2 = "P$(freq2)hm$(splits[2])"

# We'll apply the likelihood masks to the simulations.

maskfile1 = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid1)_mask$(masktype1).fits")
maskfile2 = joinpath(config["scratch"], "masks", "$(run_name)_$(mapid2)_mask$(masktype2).fits")
mask1 = readMapFromFITS(maskfile1, 1, Float64)
mask2 = readMapFromFITS(maskfile2, 1, Float64)

# Next, we generate the mode-coupling matrix.

if spec == "EE"
    @time M = mcm(:EE_BB, map2alm(mask1), map2alm(mask2))
else
    @time M = mcm(Symbol(X,Y), map2alm(mask1), map2alm(mask2))
end

# We write a function that does one iteration of this signal-only simulation and then
# estimates spectra from it.

## map T,E,B => 1,2,3
channelindex(X) = findfirst(first(X), "TEB")

function sim_iteration(𝐂, m1, m2, a1, a2, M, spec::String)
    ## get indices of the spectrum
    c₁, c₂ = channelindex(spec[1]), channelindex(spec[2])

    ## zero out alms
    for i in 1:3
        fill!(a1[i].alm, 0.0)
        fill!(a2[i].alm, 0.0)
    end

    ## synthesize polarized spectrum into m1
    synalm!(𝐂, [a1[c₁], a2[c₂]])
    alm2map!(a1, m1)
    alm2map!(a2, m2)

    ## same signal, but different masks
    mask!(m1, mask1, mask1)
    mask!(m2, mask2, mask2)

    ## subtract monopole if TT
    if spec[1] == 'T'
        monopole, dipole = fitdipole(m1.i * mask1)
        subtract_monopole_dipole!(m1.i, monopole, dipole)
    end
    if spec[2] == 'T'
        monopole, dipole = fitdipole(m2.i * mask2)
        subtract_monopole_dipole!(m2.i, monopole, dipole)
    end

    ## apply pixel weights and then map2alm
    Healpix.applyFullWeights!(m1)
    Healpix.applyFullWeights!(m2)
    map2alm!(m1, a1; niter=0)
    map2alm!(m2, a2; niter=0)

    if spec == "EE"
        pCl_EE = SpectralVector(alm2cl(a1[c₁], a2[c₂]))
        pCl_BB = SpectralVector(zeros(length(pCl_EE)))
        @spectra Cl_EE, Cl_BB = M \ [pCl_EE; pCl_BB]
        return Cl_EE
    end

    ## otherwise easy mode coupling
    pCl_XY = SpectralVector(alm2cl(a1[c₁], a2[c₂]))
    return M \ pCl_XY
end

# Let's take a look at a simulated spectrum. Note the deviation at ``\ell > 2n_{\mathrm{side}}``.
# This is due to pixelization, and is unimportant. For example, in the full Planck analysis we would
# expect to start seeing these problems at ``\ell \sim 2n_{\mathrm{side}} = 4096``.

if "--plot" in ARGS
    @time simTT = sim_iteration(𝐂, m1, m2, a1, a2, M, "TT")
    plot(simTT .* eachindex(simTT).^2, label="sim")
    plot!(signal12["TT"][1:(lmax_planck+1)] .* eachindex(0:lmax_planck).^2,
        xlim=(0,lmax_planck), label="input")
end

# This script generates many simulations. We'll run this for `nsims` iterations and
# then save the spectra to disk.

simpath = joinpath(config["scratch"], "signalsims", "$(freq1)_$(freq2)_$(spec)")
mkpath(simpath)

for sim_index in 1:parse(Int,nsims)
    @time cl = sim_iteration(𝐂, m1, m2, a1, a2, M, spec)
    @save "$(simpath)/$(uuid4()).jld2" cl=cl
end

max_forks_repo_name: sgiardie/PSpipe


text:
empty_message_str = "Empty message passed to encode!"
big_message_str = "Can only safely encode 254 bytes at a time"

"""
    encode(ard::Arduino, payload::Vector{UInt8})
Zero Allocation encoding of a message block
"""
function encode(ard::Arduino, payload::AbstractVector{UInt8})
    length(payload) == 0 && error("empty_message_str")
    length(payload) > 254 && error("big_message_str")

    n = length(payload)
    ard.msg_out_length = n + 2

    ind = 0x01
    acc = 0x01
    for x in Iterators.reverse(payload)
        if iszero(x)
            ard.msg_out_buffer[ind] = acc
            acc = 0x00
        else
            ard.msg_out_buffer[ind] = x
        end

        ind += 0x01
        acc += 0x01
    end
    ard.msg_out_buffer[ard.msg_out_length-1] = acc

    # Reverse the msg_buffer
    reverse!(ard.msg_out_buffer, 1, ard.msg_out_length-1)
    # Add on end flag to message
    ard.msg_out_buffer[ard.msg_out_length] = 0x00
    # Return a view into the msg buffer of just critical part of the buffer
    return view(ard.msg_out_buffer, 1:ard.msg_out_length)
end

"""
    decode(msg)
Zero allocation decoding of message block.
"""
function decode(ard::Arduino, msg::AbstractVector{UInt8})
    length(msg) == 0 && error("empty_message_str")
    length(msg) > 256 && error(big_message_str)

    push_ind = 1
    ard.msg_in_buffer

    pop_ind = 1
    n = msg[pop_ind]
    pop_ind += 1

    c = 0

    b = msg[pop_ind]
    pop_ind += 1

    while b ≠ 0x00 && pop_ind <= 256 && push_ind <= 256
        c += 1
        if c < n
            ard.msg_in_buffer[push_ind] = b
            push_ind += 1
        else
            ard.msg_in_buffer[push_ind] = 0x00
            push_ind += 1

            n = b
            c = 0
        end

        b = msg[pop_ind]
        pop_ind += 1
    end

    ard.msg_in_length = push_ind - 1

    return view(ard.msg_in_buffer, 1:ard.msg_in_length)
end


text:
# Point3D, FreeVector3D. The difference is that a FreeVector3D is only rotated when its frame is changed,
# whereas a Point3D is also translated
for VectorType in (:FreeVector3D, :Point3D)
    @eval begin
        # TODO: consider storing as a homogeneous vector
        struct $VectorType{V<:AbstractVector}
            v::V
            frame::CartesianFrame3D

            @inline function $VectorType(frame::CartesianFrame3D, v::V) where {V}
                @boundscheck length(v) == 3 || throw(DimensionMismatch())
                new{V}(v, frame)
            end
        end

        $VectorType(::Type{T}, frame::CartesianFrame3D) where {T} = $VectorType(frame, zero(SVector{3, T}))
        $VectorType(frame::CartesianFrame3D, x::Number, y::Number, z::Number) = $VectorType(frame, SVector(x, y, z))

        Base.convert(::Type{$VectorType{V}}, p::$VectorType{V}) where {V} = p
        Base.convert(::Type{$VectorType{V}}, p::$VectorType) where {V} = $VectorType(p.frame, convert(V, p.v))

        Base.zero(p::$VectorType) = $VectorType(p.frame, zero(p.v))

        Base.:/(p::$VectorType, s::Number) = $VectorType(p.frame, p.v / s)
        Base.:*(p::$VectorType, s::Number) = $VectorType(p.frame, p.v * s)
        Base.:*(s::Number, p::$VectorType) = $VectorType(p.frame, s * p.v)
        Base.:-(p::$VectorType) = $VectorType(p.frame, -p.v)

        Random.rand(::Type{$VectorType}, ::Type{T}, frame::CartesianFrame3D) where {T} = $VectorType(frame, rand(SVector{3, T}))
        Base.show(io::IO, p::$VectorType) = print(io, "$($(string(VectorType))) in \"$(string(p.frame))\": $(p.v)")
        Base.isapprox(x::$VectorType, y::$VectorType; atol::Real = 1e-12) = x.frame == y.frame && isapprox(x.v, y.v; atol = atol)

        """
        $(SIGNATURES)

        Return `x` transformed to `CartesianFrame3D` `t.from`.
        """
        transform(x::$VectorType, t::Transform3D) = t * x

        Base.eltype(::Type{$VectorType{V}}) where {V} = eltype(V)
    end
end

"""
$(TYPEDEF)

A `Point3D` represents a position in a given coordinate system.

A `Point3D` is a [bound vector](https://en.wikipedia.org/wiki/Euclidean_vector#Overview).
Applying a `Transform3D` to a `Point3D` both rotates and translates the
`Point3D`.
"""
Point3D

"""
$(TYPEDEF)

A `FreeVector3D` represents a [free vector](https://en.wikipedia.org/wiki/Euclidean_vector#Overview).

Examples of free vectors include displacements and velocities of points.

Applying a `Transform3D` to a `FreeVector3D` only rotates the `FreeVector3D`.
"""
FreeVector3D

# Point3D-specific
Base.:-(p1::Point3D, p2::Point3D) = begin @framecheck(p1.frame, p2.frame); FreeVector3D(p1.frame, p1.v - p2.v) end
Base.:*(t::Transform3D, point::Point3D) = begin @framecheck(t.from, point.frame); Point3D(t.to, rotation(t) * point.v + translation(t)) end
Base.:\(t::Transform3D, point::Point3D) = begin @framecheck point.frame t.to; Point3D(t.from, rotation(t) \ (point.v - translation(t))) end

# FreeVector3D-specific
FreeVector3D(p::Point3D) = FreeVector3D(p.frame, p.v)
Base.:-(v1::FreeVector3D, v2::FreeVector3D) = begin @framecheck(v1.frame, v2.frame); FreeVector3D(v1.frame, v1.v - v2.v) end
LinearAlgebra.cross(v1::FreeVector3D, v2::FreeVector3D) = begin @framecheck(v1.frame, v2.frame); FreeVector3D(v1.frame, v1.v × v2.v) end
LinearAlgebra.dot(v1::FreeVector3D, v2::FreeVector3D) = begin @framecheck(v1.frame, v2.frame); v1.v ⋅ v2.v end
Base.:*(t::Transform3D, vector::FreeVector3D) = begin @framecheck(t.from, vector.frame); FreeVector3D(t.to, rotation(t) * vector.v) end
Base.:\(t::Transform3D, point::FreeVector3D) = begin @framecheck point.frame t.to; FreeVector3D(t.from, rotation(t) \ point.v) end
LinearAlgebra.norm(v::FreeVector3D) = norm(v.v)
LinearAlgebra.normalize(v::FreeVector3D, p = 2) = FreeVector3D(v.frame, normalize(v.v, p))

# Mixed Point3D and FreeVector3D
Base.:+(p1::FreeVector3D, p2::FreeVector3D) = begin @framecheck(p1.frame, p2.frame); FreeVector3D(p1.frame, p1.v + p2.v) end
Base.:+(p::Point3D, v::FreeVector3D) = begin @framecheck(p.frame, v.frame); Point3D(p.frame, p.v + v.v) end
Base.:+(v::FreeVector3D, p::Point3D) = p + v
Base.:-(p::Point3D, v::FreeVector3D) = begin @framecheck(p.frame, v.frame); Point3D(p.frame, p.v - v.v) end
LinearAlgebra.cross(p::Point3D, v::FreeVector3D) = begin @framecheck(p.frame, v.frame); FreeVector3D(p.frame, p.v × v.v) end


text:
push!(LOAD_PATH, "Modules")
using Colors
using IntervalSets
using Revise
using JuliRay
Base.@irrational ° 0.0174532925199432957692369076848861271344 (big(pi)/big(180))

N=60
for i in 1
    𝒑(u) = [u...,u[1]^2+u[2]^2]
    D1 = -1..1
    D2 = -1..1
    mesh = (m,m)
    obj = rgbColor(ParametricSurface(𝒑, D1, D2, mesh=mesh, smooth=false), RGB(0.1,1,0.1))
    lng = 30°
    lat = 20°
    lookat = [0,0,0]
    zoom = 0.3
    render(obj,camera=LngLatCamera(lng=lng,lat=lat,pers=0.2,zoom=zoom,lookat=lookat),name="Paraboloid",index=i)
end

max_stars_repo_name: hyrodium/JuliRay

max_issues_repo_name: hyrodium/JuliRay

max_forks_repo_name: hyrodium/JuliRay


text:
# PSO

# Import Packages

using Interpolations
using PyPlot

pygui(true)

# Form Camber

x = y = 0: 10: 100
z = rand(length(x), length(y)) .* 10
# fig1 = surf(x, y, z)
xq = yq = 0: 100
interp_cubic = CubicSplineInterpolation((x, y), z)
zq = interp_cubic(xq, yq)
fig1 = surf(xq, yq, zq)

# Create Pops

N = 100
iter = 1000
ω = 0.8
c1 = 1
c2 = 1
poplimit = [x[1], x[end]]
vlimit = [-3, 3]

pop = zeros(N, 3)
v = zeros(N, 2)
localBest = zeros(N, 3)
globalBest = zeros(N, 3)

for ii = 1: N

    global pop, v, localBest

    pop[ii, 1: 2] = rand(1, 2) * 100
    pop[ii, 3] = interp_cubic(pop[ii, 1], pop[ii, 2])
    localBest[ii, :] = pop[ii, :]
    v[ii, :] = rand(1, 2) * 2 .* vlimit[1] .+ vlimit[2]

end

scatter3D(pop[:, 1], pop[:, 2], pop[:, 3])

for jj = 1: 2

    for kk = 1: N

        global pop, v, localBest, globalBest

        # Update Pops

        popnew = zeros(1, 3)
        popnew[1: 2] = pop[kk, 1: 2] .+ v[kk, :]

        if (poplimit[1] ≤ popnew[1, 1] ≤ poplimit[2]) && (poplimit[1] ≤ popnew[1, 2] ≤ poplimit[2])

            popnew[3] = interp_cubic(popnew[1], popnew[2])
            pop[kk, :] = popnew[1, :]

            if pop[kk, 3] ≥ localBest[kk, 3]

                localBest[kk, :] = pop[kk, :]

                if localBest[kk, 3] ≥ globalBest[3]

                    globalBest = localBest[kk, :]

                end

            end

        else

            popnew[3] = 0
            pop[kk, :] = popnew[1, :]

        end

        # Update v

        vnew = zeros(1, 2)
        vrand = rand(2, 2) .* 2 .- 1
        vchange = [c1 * (globalBest[1] - pop[kk, 1]) c2 * (localBest[kk, 1] - pop[kk, 1]); c1 * (globalBest[2] - pop[kk, 2]) c2 * (localBest[kk, 2] - pop[kk, 2])]
        vnew = ω * v[kk, :] .+ vrand .* vchange

        if (vlimit[1] ≤ vnew[1, 1] ≤ vlimit[2]) && (vlimit[1] ≤ vnew[1, 2] ≤ vlimit[2])

            v[kk, :] = vnew[1, :]

        end

    end

end

scatter3D(globalBest[1], globalBest[2], globalBest[3])

println("max is at $(globalBest)")

max_stars_repo_name: mizu-bai/Try-Julia

max_issues_repo_name: mizu-bai/Data-Science-Demos

max_forks_repo_name: mizu-bai/Data-Science-Demos


text:
# # Wright-Fisher model of evolution

# This is one of the simplest models of population genetics that demonstrates the
# use of [`sample!`](@ref).
# We implement a simple case of the model where we study haploids (cells with a single set
# of chromosomes) while for simplicity, focus only on one locus (a specific gene).
# In this example we will be dealing with a population of constant size.

# ## A neutral model

# * Imagine a population of `n` haploid individuals.
# * At each generation, `n` offsprings replace the parents.
# * Each offspring chooses a parent at random and inherits its genetic material.

using Agents
n = 100
nothing # hide

# Let's define an agent. The genetic value of an agent is a number (`trait` field).
mutable struct Haploid <: AbstractAgent
    id::Int
    trait::Float64
end

# And make a model without any spatial structure:
m = ABM(Haploid)

# Create `n` random individuals:
for i in 1:n
    add_agent!(m, rand())
end

# To create a new generation, we can use the `sample!` function. It chooses
# random individuals with replacement from the current individuals and updates
# the model. For example:
sample!(m, nagents(m))
nothing # hide

# The model can be run for many generations and we can collect the average trait
# value of the population. To do this we will use a model-step function (see [`step!`](@ref))
# that utilizes [`sample!`](@ref):

modelstep_neutral!(m) = sample!(m, nagents(m))
nothing # hide

# We can now run the model and collect data. We use `dummystep` for the agent-step
# function (as the agents perform no actions).
using Statistics: mean

data, _ = run!(m, dummystep, modelstep_neutral!, 20; adata = [(:trait, mean)])
data

# As expected, the average value of the "trait" remains around 0.5.

# ## A model with selection

# We can sample individuals according to their trait values, supposing that their
# fitness is correlated with their trait values.

m = ABM(Haploid)
for i in 1:100
    add_agent!(m, rand())
end

modelstep_selection!(m::ABM) = sample!(m, nagents(m), :trait)

data, _ = run!(m, dummystep, modelstep_selection!, 20; adata = [(:trait, mean)])
data

# Here we see that as time progresses, the trait becomes closer and closer to 1,
# which is expected - since agents with higher traits have higher probability of being
# sampled for the next "generation".


text:
using Documenter, TypedTables

makedocs(;
    modules=[TypedTables],
    format=Documenter.HTML(),
    pages=[
        "Home" => "index.md",
        "Quick start tutorial" => "man/tutorial.md",
        "Table types" => [
            "Table" => "man/table.md",
            "FlexTable " => "man/flextable.md",
            "DictTable " => "man/dicttable.md",
        ],
        "Data storage" => [
            "Data representation" => "man/data.md",
            "Input and output" => "man/io.md",
        ],
        "Basic data manipulation" => [
            "Mapping data" => "man/map.md",
            "Finding data" => "man/filter.md",
            "Reducing data" => "man/reduce.md",
        ],
        "Grouping and joining data" => [
            "Grouping data" => "man/group.md",
            "Joining data" => "man/join.md",
        ],
        "Accelerations" => [
            "Acceleration indices" => "man/acceleratedarrays.md",
        ],
        "API reference" => "man/reference.md"
    ],
    repo="https://github.com/JuliaData/TypedTables.jl/blob/{commit}{path}#L{line}",
    sitename="TypedTables.jl",
)

deploydocs(;
    repo="github.com/JuliaData/TypedTables.jl",
    devbranch = "main"
)


text:
function CartesianOrbit(arg0::Orbit)
    return CartesianOrbit((Orbit,), arg0)
end

function CartesianOrbit(arg0::PVCoordinates, arg1::Frame, arg2::AbsoluteDate, arg3::jdouble)
    return CartesianOrbit((PVCoordinates, Frame, AbsoluteDate, jdouble), arg0, arg1, arg2, arg3)
end

function CartesianOrbit(arg0::TimeStampedPVCoordinates, arg1::Frame, arg2::jdouble)
    return CartesianOrbit((TimeStampedPVCoordinates, Frame, jdouble), arg0, arg1, arg2)
end

function add_kepler_contribution(obj::CartesianOrbit, arg0::PositionAngle, arg1::jdouble, arg2::Vector{jdouble})
    return jcall(obj, "addKeplerContribution", void, (PositionAngle, jdouble, Vector{jdouble}), arg0, arg1, arg2)
end

function get_a(obj::CartesianOrbit)
    return jcall(obj, "getA", jdouble, ())
end

function get_a_dot(obj::CartesianOrbit)
    return jcall(obj, "getADot", jdouble, ())
end

function get_date(obj::Orbit)
    return jcall(obj, "getDate", AbsoluteDate, ())
end

function get_e(obj::CartesianOrbit)
    return jcall(obj, "getE", jdouble, ())
end

function get_e_dot(obj::CartesianOrbit)
    return jcall(obj, "getEDot", jdouble, ())
end

function get_equinoctial_ex(obj::CartesianOrbit)
    return jcall(obj, "getEquinoctialEx", jdouble, ())
end

function get_equinoctial_ex_dot(obj::CartesianOrbit)
    return jcall(obj, "getEquinoctialExDot", jdouble, ())
end

function get_equinoctial_ey(obj::CartesianOrbit)
    return jcall(obj, "getEquinoctialEy", jdouble, ())
end

function get_equinoctial_ey_dot(obj::CartesianOrbit)
    return jcall(obj, "getEquinoctialEyDot", jdouble, ())
end

function get_frame(obj::Orbit)
    return jcall(obj, "getFrame", Frame, ())
end

function get_hx(obj::CartesianOrbit)
    return jcall(obj, "getHx", jdouble, ())
end

function get_hx_dot(obj::CartesianOrbit)
    return jcall(obj, "getHxDot", jdouble, ())
end

function get_hy(obj::CartesianOrbit)
    return jcall(obj, "getHy", jdouble, ())
end

function get_hy_dot(obj::CartesianOrbit)
    return jcall(obj, "getHyDot", jdouble, ())
end

function get_i(obj::CartesianOrbit)
    return jcall(obj, "getI", jdouble, ())
end

function get_i_dot(obj::CartesianOrbit)
    return jcall(obj, "getIDot", jdouble, ())
end

function get_jacobian_wrt_cartesian(obj::Orbit, arg0::PositionAngle, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getJacobianWrtCartesian", void, (PositionAngle, Vector{Vector{jdouble}}), arg0, arg1)
end

function get_jacobian_wrt_parameters(obj::Orbit, arg0::PositionAngle, arg1::Vector{Vector{jdouble}})
    return jcall(obj, "getJacobianWrtParameters", void, (PositionAngle, Vector{Vector{jdouble}}), arg0, arg1)
end

function get_keplerian_mean_motion(obj::Orbit)
    return jcall(obj, "getKeplerianMeanMotion", jdouble, ())
end

function get_keplerian_period(obj::Orbit)
    return jcall(obj, "getKeplerianPeriod", jdouble, ())
end

function get_le(obj::CartesianOrbit)
    return jcall(obj, "getLE", jdouble, ())
end

function get_le_dot(obj::CartesianOrbit)
    return jcall(obj, "getLEDot", jdouble, ())
end

function get_lm(obj::CartesianOrbit)
    return jcall(obj, "getLM", jdouble, ())
end

function get_lm_dot(obj::CartesianOrbit)
    return jcall(obj, "getLMDot", jdouble, ())
end

function get_lv(obj::CartesianOrbit)
    return jcall(obj, "getLv", jdouble, ())
end

function get_lv_dot(obj::CartesianOrbit)
    return jcall(obj, "getLvDot", jdouble, ())
end

function get_mu(obj::Orbit)
    return jcall(obj, "getMu", jdouble, ())
end

function get_pv_coordinates(obj::Orbit)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, ())
end

function get_pv_coordinates(obj::Orbit, arg0::AbsoluteDate, arg1::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (AbsoluteDate, Frame), arg0, arg1)
end

function get_pv_coordinates(obj::Orbit, arg0::Frame)
    return jcall(obj, "getPVCoordinates", TimeStampedPVCoordinates, (Frame,), arg0)
end

function get_type(obj::CartesianOrbit)
    return jcall(obj, "getType", OrbitType, ())
end

function has_derivatives(obj::CartesianOrbit)
    return jcall(obj, "hasDerivatives", jboolean, ())
end

function interpolate(obj::CartesianOrbit, arg0::AbsoluteDate, arg1::Stream)
    return jcall(obj, "interpolate", CartesianOrbit, (AbsoluteDate, Stream), arg0, arg1)
end

function interpolate(obj::TimeInterpolable, arg0::AbsoluteDate, arg1::Collection)
    return jcall(obj, "interpolate", TimeInterpolable, (AbsoluteDate, Collection), arg0, arg1)
end

function shifted_by(obj::CartesianOrbit, arg0::jdouble)
    return jcall(obj, "shiftedBy", CartesianOrbit, (jdouble,), arg0)
end

function to_string(obj::CartesianOrbit)
    return jcall(obj, "toString", JString, ())
end


text:
using JuMP
using ParameterJuMP
using Polyhedra
mutable struct BalancedRealPolytope{T, VT <: AbstractVector{T}, D<:Polyhedra.FullDim} <: Polyhedra.VPolytope{T}
    d::D
    points::Vector{VT}
    optimizer_constructor
    model::Union{Nothing, MOI.ModelLike}
    sum_con::Union{Nothing, MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}}
    z::Union{Nothing, Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}}}
    t_0::Union{Nothing, MOI.VariableIndex}
    function BalancedRealPolytope{T, VT, D}(
        d::Polyhedra.FullDim, points::Polyhedra.PointIt,
        optimizer_constructor=nothing) where {T, VT, D}
        new{T, VT, D}(Polyhedra.FullDim_convert(D, d), Polyhedra.lazy_collect(points), optimizer_constructor, nothing, nothing, nothing, nothing)
    end
end
function BalancedRealPolytope(d::Polyhedra.FullDim, points::Polyhedra.PointIt, args...)
    return BalancedRealPolytope{Polyhedra.coefficient_type(eltype(points)), eltype(points),
                                typeof(d)}(d, points, args...)
end

Polyhedra.FullDim(v::BalancedRealPolytope) = v.d
Polyhedra.vvectortype(::Type{<:BalancedRealPolytope{T, AT}}) where {T, AT} = AT
function Polyhedra.similar_type(PT::Type{<:BalancedRealPolytope}, d::Polyhedra.FullDim, ::Type{T}) where {T}
    BalancedRealPolytope{T, Polyhedra.similar_type(Polyhedra.vvectortype(PT), d, T), typeof(d)}
end

Polyhedra.vreptype(::Type{BalancedRealPolytope{T, AT, D}}) where {T, AT, D} = Polyhedra.Hull{T, AT, D}
Polyhedra.fulltype(::Type{BalancedRealPolytope{T, AT, D}}) where {T, AT, D} = Polyhedra.Hull{T, AT, D}

Base.length(idxs::Polyhedra.PointIndices{T, <:BalancedRealPolytope{T}}) where {T} = 2length(idxs.rep.points)
Base.isempty(idxs::Polyhedra.PointIndices{T, <:BalancedRealPolytope{T}}) where {T} = isempty(idxs.rep.points)
function Polyhedra.startindex(idxs::Polyhedra.PointIndices{T, <:BalancedRealPolytope{T}}) where {T}
    if isempty(idxs.rep.points)
        return nothing
    else
        return eltype(idxs)(1)
    end
end
function Base.get(rep::BalancedRealPolytope{T}, idx::Polyhedra.PointIndex{T}) where {T}
    return sign(idx.value) * rep.points[abs(idx.value)]
end
function Polyhedra.nextindex(rep::BalancedRealPolytope{T}, idx::Polyhedra.PointIndex{T}) where {T}
    if idx.value < 0
        if -idx.value >= length(rep.points)
            return nothing
        else
            return typeof(idx)(idx.value - 1)
        end
    else
        if idx.value >= length(rep.points)
            return typeof(idx)(-1)
        else
            return typeof(idx)(idx.value + 1)
        end
    end
end

function _abs_constraint(model::MOI.ModelLike, t::MOI.VariableIndex, q::MOI.VariableIndex, T::Type)
    MOI.add_constraint(model, MOI.Utilities.operate(-, T, q, t), MOI.GreaterThan(zero(T)))
    MOI.add_constraint(model, MOI.Utilities.operate(+, T, q, t), MOI.GreaterThan(zero(T)))
end

function _build_model(brp::BalancedRealPolytope{T}) where T
    @assert brp.model === nothing
    n = length(brp.points)
    brp.model = _instantiate(brp.optimizer_constructor, T)
    t = MOI.add_variables(brp.model, n)
    q = MOI.add_variables(brp.model, n)
    # |t| ≤ q
    for i in 1:n
        _abs_constraint(brp.model, t[i], q[i], T)
    end
    # sum |t| ≤ sum q ≤ 1
    brp.sum_con = MOI.add_constraint(brp.model, _sum(q, T), MOI.LessThan(one(T)))
    return _convex_combination(t, brp.points, brp.d)
end

function _parametrized_model(brp::BalancedRealPolytope, v)
    hull = _build_model(brp)
    brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.EqualTo(v[i])) for i in 1:brp.d]
end
function _ratio_model(brp::BalancedRealPolytope{T}, v) where T
    hull = _build_model(brp)
    brp.t_0 = MOI.add_variable(brp.model)
    for i in 1:brp.d
        push!(hull[i].terms, MOI.ScalarAffineTerm(-v[i], brp.t_0))
    end
    brp.z = [MOI.add_constraint(brp.model, hull[i], MOI.EqualTo(zero(T))) for i in 1:brp.d]
end
function _update_model(p::BalancedRealPolytope{T}, v) where T
    t = MOI.add_variable(p.model)
    q = MOI.add_variable(p.model)
    _abs_constraint(p.model, t, q, T)
    MOI.modify(p.model, p.sum_con, MOI.ScalarCoefficientChange(q, one(T)))
    for i in 1:p.d
        MOI.modify(p.model, p.z[i], MOI.ScalarCoefficientChange(t, v[i]))
    end
end

function _fix(brp::BalancedRealPolytope, v)
    for i in 1:brp.d
        if brp.t_0 === nothing
            MOI.set(brp.model, MOI.ConstraintSet(), brp.z[i], MOI.EqualTo(v[i]))
        else
            MOI.modify(brp.model, brp.z[i], MOI.ScalarCoefficientChange(brp.t_0, -v[i]))
        end
    end
end


text:
"""
    QuasiSteady{Order} <: NoStateModel

2D quasi-steady aerodynamic model with parameters ``a, b, a_0, \\alpha_0, c_{d0}, c_{m0}``.
"""
struct QuasiSteady{Order} <: NoStateModel end

"""
    Steady <: NoStateModel

2D steady aerodynamic model with parameters ``a, b, a_0, \\alpha_0, c_{d0}, c_{m0}``.
"""
const Steady = QuasiSteady{0}

# --- Constructors --- #

"""
    Steady()

Initialize an object of type [`Steady`](@ref) which represents a 2D
steady aerodynamic model.
"""
Steady()

"""
    QuasiSteady()

Initialize an object of type [`QuasiSteady`](@ref) which represents a 2D
quasi-steady aerodynamic model.
"""
QuasiSteady() = QuasiSteady{2}()

# --- Traits --- #

number_of_parameters(::Type{<:QuasiSteady}) = 6

# --- Convenience Methods --- #

function set_parameters!(p, model::QuasiSteady; a, b, a0, alpha0, cd0, cm0)

    p[1] = a
    p[2] = b
    p[3] = a0
    p[4] = alpha0
    p[5] = cd0
    p[6] = cm0

    return p
end

function separate_parameters(model::QuasiSteady, p)

    return (a=p[1], b=p[2], a0=p[3], alpha0=p[4], cd0=p[5], cm0=p[6])
end

# --- Internal Methods --- #

# steady state loads
function quasisteady0_loads(a, b, ρ, a0, α0, cd0, cm0, u, v)
    # normal force at reference point
    N = a0*ρ*b*u*(v - u*α0)
    # axial force at reference point
    A = -a0*ρ*b*(v - u*α0)^2 + ρ*b*u^2*cd0
    # moment at reference point
    M = 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N

    return SVector(N, A, M)
end

function quasisteady0_loads_a(a, b, ρ, a0, α0, u, v)

    N = a0*ρ*b*u*(v - u*α0)

    N_a = 0

    A_a = 0

    M_a = b*N

    return SVector(N_a, A_a, M_a)
end

function quasisteady0_loads_b(a, b, ρ, a0, α0, cd0, cm0, u, v)

    N = a0*ρ*b*u*(v - u*α0)

    N_b = a0*ρ*u*(v - u*α0)

    A_b = -a0*ρ*(v - u*α0)^2 + ρ*u^2*cd0

    M_b = 4*ρ*b*u^2*cm0 + (1/2 + a)*N + (b/2 + a*b)*N_b

    return SVector(N_b, A_b, M_b)
end

function quasisteady0_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, u, v)
    N_ρ = a0*b*u*(v - u*α0)
    A_ρ = -a0*b*(v - u*α0)^2 + b*u^2*cd0
    M_ρ = 2*b^2*u^2*cm0 + (b/2 + a*b)*N_ρ
    return SVector(N_ρ, A_ρ, M_ρ)
end

function quasisteady0_loads_a0(a, b, ρ, α0, u, v)
    # normal force at reference point
    N_a0 = ρ*b*u*(v - u*α0)
    A_a0 = -ρ*b*(v - u*α0)^2
    M_a0 = (b/2 + a*b)*N_a0

    return SVector(N_a0, A_a0, M_a0)
end

function quasisteady0_loads_α0(a, b, ρ, a0, α0, u, v)
    N_α0 = -a0*ρ*b*u^2
    A_α0 = 2*a0*b*u*(v - u*α0)
    M_α0 = (b/2 + a*b)*N_α0

    return SVector(N_α0, A_α0, M_α0)
end

function quasisteady0_loads_cd0(b, ρ, u)
    N_cd0 = 0
    A_cd0 = ρ*b*u^2
    M_cd0 = 0
    return SVector(N_cd0, A_cd0, M_cd0)
end

function quasisteady0_loads_cm0(b, ρ, u)
    N_cm0 = 0
    A_cm0 = 0
    M_cm0 = 2*ρ*b^2*u^2
    return SVector(N_cm0, A_cm0, M_cm0)
end

function quasisteady0_loads_u(a, b, ρ, a0, α0, cd0, cm0, u, v)
    N_u = a0*ρ*b*v - 2*a0*ρ*b*u*α0
    A_u = 2*a0*α0*ρ*b*(v - u*α0)^2 + 2*ρ*b*u*cd0
    M_u = 4*ρ*b^2*u*cm0 + (b/2 + a*b)*N_u
    return SVector(N_u, A_u, M_u)
end

function quasisteady0_loads_v(a, b, ρ, a0, α0, u, v)
    N_v = a0*ρ*b*u
    A_v = -2*a0*ρ*b*(v - u*α0)
    M_v = (b/2 + a*b)*N_v
    return SVector(N_v, A_v, M_v)
end

# circulatory loads + added mass effects, neglecting accelerations
function quasisteady1_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # normal force at reference point
    N = tmp1*(v + d*ω - u*α0) + tmp2*u/b*ω
    # axial force at the reference point
    A = -a0*ρ*b*(v + d*ω - u*α0)^2 + ρ*b*u^2*cd0
    # moment at reference point
    M = -tmp2*u*ω + 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N

    return SVector(N, A, M)
end

function quasisteady1_loads_a(a, b, ρ, a0, α0, u, v, ω)

    tmp1 = a0*ρ*u*b
    tmp2 = pi*ρ*b^3

    N = tmp1*(v + (b/2 - a*b)*ω - u*α0) + tmp2*u/b*ω
    N_a = -tmp1*b*ω
    A_a = 2*a0*ρ*b^2*(v + (b/2 - a*b)*ω - u*α0)*ω
    M_a = b*N + (b/2 + a*b)*N_a

    return SVector(N_a, A_a, M_a)
end

function quasisteady1_loads_b(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)

    tmp1 = a0*ρ*u*b
    tmp1_b = a0*ρ*u

    tmp2 = pi*ρ*b^3
    tmp2_b = 3*pi*ρ*b^2

    d = b/2 - a*b
    d_b = 1/2 - a

    N = tmp1*(v + d*ω - u*α0) + tmp2*u/b*ω

    N_b = tmp1_b*(v + d*ω - u*α0) + tmp1*d_b*ω + tmp2_b*u/b*ω - tmp2*u/b^2*ω
    A_b = -2*a0*ρ*b*(v + d*ω - u*α0)*d_b - a0*ρ*(v + d*ω - u*α0)^2  + ρ*u^2*cd0
    M_b = -tmp2_b*u*ω + 4*ρ*b*u^2*cm0 + (1/2 + a)*N + (b/2 + a*b)*N_b

    return SVector(N_b, A_b, M_b)
end

function quasisteady1_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)

    tmp1_ρ = a0*u*b

    tmp2_ρ = pi*b^3

    d = b/2 - a*b

    N_ρ = tmp1_ρ*(v + d*ω - u*α0) + tmp2_ρ*u/b*ω
    A_ρ = -a0*b*(v + d*ω - u*α0)^2 + b*u^2*cd0
    M_ρ = -tmp2_ρ*u*ω + 2*b^2*u^2*cm0 + (b/2 + a*b)*N_ρ

    return SVector(N_ρ, A_ρ, M_ρ)
end

function quasisteady1_loads_a0(a, b, ρ, a0, α0, u, v, ω)

    tmp1_a0 = ρ*u*b

    tmp2 = pi*ρ*b^3

    d = b/2 - a*b

    N_a0 = tmp1_a0*(v + d*ω - u*α0)
    A_a0 = -ρ*b*(v + d*ω - u*α0)^2
    M_a0 = (b/2 + a*b)*N_a0

    return SVector(N_a0, A_a0, M_a0)
end

function quasisteady1_loads_α0(a, b, ρ, a0, α0, u, v, ω)
    tmp1 = a0*ρ*u*b
    N_α0 = -tmp1*u
    A_α0 = 2*a0*ρ*b*u*(v + (b/2 - a*b)*ω - u*α0)
    M_α0 = (b/2 + a*b)*N_α0
    return SVector(N_α0, A_α0, M_α0)
end

function quasisteady1_loads_cd0(b, ρ, u)
    N_cd0 = 0
    A_cd0 = ρ*b*u^2
    M_cd0 = 0
    return SVector(N_cd0, A_cd0, M_cd0)
end

function quasisteady1_loads_cm0(b, ρ, u)
    N_cm0 = 0
    A_cm0 = 0
    M_cm0 = 2*ρ*b^2*u^2
    return SVector(N_cm0, A_cm0, M_cm0)
end

function quasisteady1_loads_u(a, b, ρ, a0, α0, cd0, cm0, u, v, ω)
    tmp1 = a0*ρ*u*b
    tmp1_u = a0*ρ*b
    tmp2 = pi*ρ*b^3
    d = b/2 - a*b
    N_u = tmp1_u*(v + d*ω - u*α0) - tmp1*α0 + tmp2/b*ω
    A_u = 2*a0*ρ*b*α0*(v + d*ω - u*α0) + 2*ρ*b*u*cd0
    M_u = -tmp2*ω + (b/2 + a*b)*N_u + 4*ρ*b^2*u*cm0
    return SVector(N_u, A_u, M_u)
end

function quasisteady1_loads_v(a, b, ρ, a0, α0, u, v, ω)
    d = b/2 - a*b
    N_v = a0*ρ*u*b
    A_v = -2*a0*ρ*b*(v + d*ω - u*α0)
    M_v = (b/2 + a*b)*N_v
    return SVector(N_v, A_v, M_v)
end

function quasisteady1_loads_ω(a, b, ρ, a0, u)
    tmp1 = a0*ρ*u*b
    tmp2 = pi*ρ*b^3
    d = b/2 - a*b
    N_ω = tmp1*d + tmp2*u/b
    A_ω = 2*a0*ρ*a*b^2*(v + d*ω - u*α0)
    M_ω = -tmp2*u + (b/2 + a*b)*N_ω
    return SVector(N_ω, A_ω, M_ω)
end

# quasi-steady loads + added mass effects
function quasisteady2_loads(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)
    # circulatory load factor
    tmp1 = a0*ρ*u*b
    # non-circulatory load factor
    tmp2 = pi*ρ*b^3
    # constant based on geometry
    d = b/2 - a*b
    # normal force at reference point
    N = tmp1*(v + d*ω - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    # axial force at reference point
    A = -a0*ρ*b*(v + d*ω - u*α0)^2 + ρ*b*u^2*cd0
    # moment at reference point
    M = -tmp2*(vdot/2 + u*ω + (b/8 - a*b/2)*ωdot) + 2*ρ*b^2*u^2*cm0 + (b/2 + a*b)*N

    return SVector(N, A, M)
end

function quasisteady2_loads_a(a, b, ρ, a0, α0, u, v, ω, vdot, ωdot)

    tmp1 = a0*ρ*u*b
    tmp2 = pi*ρ*b^3

    N = tmp1*(v + (b/2 - a*b)*ω - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)
    N_a = -tmp1*b*ω - tmp2*ωdot
    A_a = 2*a0*ρ*b^2*(v + (b/2 - a*b)*ω - u*α0)*ω
    M_a = tmp2*b/2*ωdot + b*N + (b/2 + a*b)*N_a

    return SVector(N_a, A_a, M_a)
end

function quasisteady2_loads_b(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)

    tmp1 = a0*ρ*u*b
    tmp1_b = a0*ρ*u

    tmp2 = pi*ρ*b^3
    tmp2_b = 3*pi*ρ*b^2

    d = b/2 - a*b
    d_b = 1/2 - a

    N = tmp1*(v + d*ω - u*α0) + tmp2*(vdot/b + u/b*ω - a*ωdot)

    N_b = tmp1_b*(v + d*ω - u*α0) + tmp1*d_b*ω +
        tmp2_b*(vdot/b + u/b*ω - a*ωdot) + tmp2*(-vdot/b^2 - u/b^2*ω)
    A_b = -2*a0*ρ*b*(v + d*ω - u*α0)*d_b - a0*ρ*(v + d*ω - u*α0)^2 + ρ*u^2*cd0
    M_b = -tmp2_b*(vdot/2 + u*ω + b*(1/8 - a/2)*ωdot) - tmp2*(1/8 - a/2)*ωdot +
        4*ρ*b*u^2*cm0 + (1/2 + a)*N + (b/2 + a*b)*N_b

    return SVector(N_b, A_b, M_b)
end

function quasisteady2_loads_ρ(a, b, ρ, a0, α0, cd0, cm0, u, v, ω, vdot, ωdot)

    tmp1_ρ = a0*u*b

    tmp2_ρ = pi*b^3

    d = b/2 - a*b

    N_ρ = tmp1_ρ*(v + d*ω - u*α0) + tmp2_ρ*(vdot/b + u/b*ω - a*ωdot)
    A_ρ = -a0*b*(v + d*ω - u*α0)^2 + b*u^2*cd0
    M_ρ = -tmp2_ρ*(vdot/2 + u*ω + (b/8 - a*b/2)*ωdot) + 2*b^2*u^2*cm0 + (b/2 + a*b)*N_ρ

    return SVector(N_ρ, A_ρ, M_ρ)
end

quasisteady2_loads_a0(args...) = quasisteady1_loads_a0(args...)

quasisteady2_loads_α0(args...) = quasisteady1_loads_α0(args...)

quasisteady2_loads_cd0(args...) = quasisteady1_loads_cd0(args...)

quasisteady2_loads_cm0(args...) = quasisteady1_loads_cm0(args...)

quasisteady2_loads_u(args...) = quasisteady1_loads_u(args...)

quasisteady2_loads_v(args...) = quasisteady1_loads_v(args...)

quasisteady2_loads_ω(args...) = quasisteady1_loads_ω(args...)

function quasisteady2_loads_vdot(a, b, ρ)
    tmp = pi*ρ*b^3
    N_vdot = tmp/b
    A_vdot = 0
    M_vdot = -tmp/2 + (b/2 + a*b)*N_vdot
    return SVector(N_vdot, A_vdot, M_vdot)
end

function quasisteady2_loads_ωdot(a, b, ρ)
    tmp = pi*ρ*b^3
    N_ωdot = -tmp*a
    A_ωdot = 0
    M_ωdot = -tmp*(b/8 - a*b/2) + (b/2 + a*b)*N_ωdot
    return SVector(N_ωdot, A_ωdot, M_ωdot)
end


text:
using Turing, ContinuousBenchmarks

include(joinpath(ContinuousBenchmarks.STAN_DATA_DIR, "lda-stan.data.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "lda-stan.model.jl"))
#include(joinpath(ContinuousBenchmarks.BENCH_DIR, "lda-stan.run.jl"))
include(joinpath(ContinuousBenchmarks.STAN_MODELS_DIR, "lda.model.jl"))

# setchunksize(100)
setadbackend(:reverse_diff)
# setadbackend(:forward_diff)

# setadsafe(false)

turnprogress(false)

for (modelc, modeln) in zip([
    # "ldamodel_vec",
    "ldamodel"
    ], [
      # "LDA-vec",
      "LDA"
      ])
    data = map(ind -> ldastandata[1][ind], ["K", "V", "M", "N", "w", "doc", "beta", "alpha"])
    @tbenchmark(HMC(2, 0.005, 10), eval(modelc), (data...))
    bench_res = @tbenchmark(HMC(3000, 0.005, 10), eval(modelc), (data...))
    bench_res[4].names = ["phi[$k]" for k in 1:ldastandata[1]["K"]]
    logd = build_logd(modeln, bench_res...)
    #logd["stan"] = lda_stan_d
    #logd["time_stan"] = lda_time
    print_log(logd)
    send_log(logd)
end


text:
module AutomatonTests

using Test
using WeaklyHard
using WeaklyHard: nbits, shift, shift_rowhit, _childexists, consword

##########################
#### Common Variables ####
##########################

k_max           = 20
n_constraints   = 6
n_tests         = 1_000
types           = [Int8, Int16, Int32, Int64, BigInt]

#####################
#### Shift tests ####
#####################

function test_shift_correctness(w::T, c) where {T <: Integer}
    newWord = (w << 1) | T(c)
    if c == M
        return newWord == shift(w, c)
    else
        return newWord - T(2^(nbits(newWord)-1)) == shift(w, c)
    end # if
end # function

function test_shift_rowhit_correctness(w::T, c, n::T) where {T <: Integer}
    newWord = (w << 1) | T(c)
    if c == M
        return newWord == shift_rowhit(w, c, n)
    else
        return ((newWord & n == n) ? n : newWord) == shift_rowhit(w, c, n)
    end # if
end # function

@testset "Shift Tests" begin
    chars = [H, M]
    for T in types
        for _ in 1:n_tests
            @test test_shift_correctness(rand(T(1):T(8*sizeof(T)-1)), rand(chars))
        end # for
    end # for
end # testset

@testset "Shift_RowHit Tests" begin
    chars = [H, M]
    for T in types
        n = consword( rand( T(1) : max( T(sizeof(T) - 1), T(1) ) ) )
        for _ in 1:n_tests
            @test test_shift_rowhit_correctness(rand(T(1):T(8*sizeof(T)-1)), rand(chars), n)
        end # for
    end # for
end # testset

####################
### Vertex tests ###
####################

function test_vertex_correctness(Lambda::Set{T}) where {T <: Constraint}
    if length(Lambda) == 1
        L = pop!(Lambda)
    else
        L = dominant_set(Lambda)
    end
    automaton = build_automaton(L)

    seq = random_sequence(automaton, 5000)
    c = is_satisfied(L, seq)
    if !c
        return false
    end

    for n in vertices(automaton)
        if (!is_satisfied(L, Int64(n.w))
                || (_childexists(n, :miss) && !is_satisfied(L, Int64(n.miss.w)))
                || !is_satisfied(L, Int64(n.hit.w)))
            return false
        end
    end

    true
end # function

@testset "Vertices Tests" begin
    constraint_names = ["RowMissConstraint", "RowHitConstraint", "AnyMissConstraint", "AnyHitConstraint"]

    tested = Set{Set{Constraint}}()
    for _ in 1:n_tests

        Lambda = Set{Constraint}()
        nbr_c  = rand(1:n_constraints)
        for _ in 1:nbr_c
            k = rand(1:k_max)
            x = rand(1:max(k-1, 1))
            lambda_type = rand(constraint_names)
            if lambda_type == "RowHitConstraint"
                x = x #floor(Int, x/2)
            end
            push!(Lambda, getfield(WeaklyHard, Symbol(lambda_type))(x, k))
        end # for

        L = dominant_set(Lambda)

        if !(L in tested)
            push!(tested, L)

            c = test_vertex_correctness(Lambda)
            @test c
            if !c

                @show "###############"
                @show "Lambda:"
                for l in Lambda
                    @show l
                end
                @show "---------------"
                @show "Lambda_star:"
                Lambda_star = dominant_set(Lambda)
                for l in Lambda_star
                    @show l
                end
                @show "###############"

                return
            end
        end
    end # for

    println(" ")
    println("Executed tests on $(length(tested)) different automata")
end # testset

############################
##### Transitions tests ####
############################

function test_transition_correctness(Lambda::Set{T}) where {T <: Constraint}

    if length(Lambda) == 1
        L = pop!(Lambda)
    else
        L = Lambda
    end

    automaton = build_automaton(L)
    trans     = transitions(automaton)

    for v in vertices(automaton)
        v1 = v.hit
        if !((v.w, v1.w, H) in trans)
            return false
        end

        if _childexists(v, :miss)
            v2 = v.miss
            if !((v.w, v2.w, M) in trans)
                return false
            end
        end
    end # for

    return true
end # function

@testset "Transitions Tests" begin
    constraint_names = ["RowMissConstraint", "RowHitConstraint", "AnyMissConstraint", "AnyHitConstraint"]

    tested = Set{Set{Constraint}}()
    for _ in 1:n_tests

        Lambda = Set{Constraint}()
        nbr_c  = rand(1:n_constraints)
        for _ in 1:nbr_c
            k = rand(1:k_max)
            x = rand(1:max(k-1, 1))
            lambda_type = rand(constraint_names)
            if lambda_type == "RowHitConstraint"
                x = x #floor(Int, x/2)
            end
            push!(Lambda, getfield(WeaklyHard, Symbol(lambda_type))(x, k))
        end # for

        L = dominant_set(Lambda)

        if !(L in tested)
            push!(tested, L)

            c = test_transition_correctness(Lambda)
            @test c
            if !c

                @show "###############"
                @show "Lambda:"
                @show Lambda
                @show "---------------"
                @show "Lambda_star:"
                Lambda_star = dominant_set(Lambda)
                @show Lambda_star
                @show "###############"

                return
            end
        end
    end # for

    println(" ")
    println("Executed tests on $(length(tested)) different automata")
end # testset

@testset "Set-construction Tests" begin
    constraint_names = ["RowMissConstraint", "RowHitConstraint", "AnyMissConstraint", "AnyHitConstraint"]
    n_constraints    = 3
    k_max            = 15
    n_tests          = 1_000

    for _ in 1:n_tests
        Lambda = Set{Constraint}()
        n = rand(1:n_constraints)
        for _ in 1:n
            k = rand(1:k_max)
            x = rand(0:k)
            lambda_type = rand(constraint_names)
            if lambda_type == "RowHitConstraint"
                x = floor(Int, x/2)
            end
            push!(Lambda, getfield(WeaklyHard, Symbol(lambda_type))(x, k))
        end # for

        Lambda_star = dominant_set(Lambda)

        G = build_automaton(Lambda_star)
        minimize_automaton!(G)
        all_seq = all_sequences(G, k_max+1) # +1 to account for RowMissConstraints
        for seq in all_seq
            if !is_satisfied(Lambda_star, seq)
                @show Lambda
                @show Lambda_star
                @show seq
                @show bitstring(seq)
            end
            @test is_satisfied(Lambda_star, seq)
        end # for
    end # for
end # testset

end # module


text:
include("helper.jl")

include("chebyshevtransform.jl")

include("lyap.jl")
include("constrained_lyap.jl")

include("bary.jl")
include("clenshaw.jl")
include("ultraspherical.jl")
include("fourier.jl")

include("hesseneigs.jl")

include("LowRankMatrix.jl")
include("AlmostBandedMatrix.jl")
include("BandedBlockBandedMatrix.jl")
include("RaggedMatrix.jl")


text:
function FieldTimeSpanMap(arg0::Object, arg1::Field)
    return FieldTimeSpanMap((Object, Field), arg0, arg1)
end

function add_valid_after(obj::FieldTimeSpanMap, arg0::Object, arg1::FieldAbsoluteDate)
    return jcall(obj, "addValidAfter", void, (Object, FieldAbsoluteDate), arg0, arg1)
end

function add_valid_before(obj::FieldTimeSpanMap, arg0::Object, arg1::FieldAbsoluteDate)
    return jcall(obj, "addValidBefore", void, (Object, FieldAbsoluteDate), arg0, arg1)
end

function get(obj::FieldTimeSpanMap, arg0::FieldAbsoluteDate)
    return jcall(obj, "get", Object, (FieldAbsoluteDate,), arg0)
end

function get_transitions(obj::FieldTimeSpanMap)
    return jcall(obj, "getTransitions", SortedSet, ())
end


text:
nbc = 3
p = ParamMPCC(nbc)

@test p.precmpcc == 1e-3
#@test p.prec_oracle == (r,s,t,prec) -> max(min(r,s,t),prec)
#@test p.rho_restart == (r,s,t,rho) -> rho
@test p.paramin == sqrt(eps(Float64))
@test p.initrst() == (1.0,0.5,0.5)
@test p.updaterst(1.0,1.0,1.0) == (0.1,0.1,0.05)
#@test p.solve_sub_pb == x -> x
@test p.tb(1.0,1.0,1.0) == 1.0
@test p.ite_max_alas == 100
@test p.ite_max_viol == 30
@test p.rho_init == ones(nbc)
@test p.rho_update == 1.5
@test p.rho_max == 4000.0
@test p.goal_viol == 0.75
@test p.ite_max_armijo == 100
@test p.tau_armijo == 0.1
@test p.armijo_update == 0.9
@test p.ite_max_wolfe == 10
@test p.tau_wolfe == 0.6
@test p.wolfe_update == 5.0
@test p.verbose == 0
#@test p.uncmin == x -> x
#@test p.penalty == x -> x
#@test p.direction == x -> x
#@test p.linesearch == x -> x


text:
using Plots
include("scf_callbacks.jl")

function default_n_bands(model)
    min_n_bands = div(model.n_electrons, filled_occupation(model))
    n_extra = model.temperature == 0 ? 0 : max(4, ceil(Int, 0.2 * min_n_bands))
    min_n_bands + n_extra
end

"""
Obtain new density ρ by diagonalizing `ham`.
"""
function next_density(ham::Hamiltonian;
                      n_bands=default_n_bands(ham.basis.model),
                      ψ=nothing, n_ep_extra=3,
                      eigensolver=lobpcg_hyper, kwargs...)
    if ψ !== nothing
        @assert length(ψ) == length(ham.basis.kpoints)
        for ik in 1:length(ham.basis.kpoints)
            @assert size(ψ[ik], 2) == n_bands + n_ep_extra
        end
    end

    # Diagonalize
    eigres = diagonalize_all_kblocks(eigensolver, ham, n_bands + n_ep_extra; guess=ψ,
                                     n_conv_check=n_bands, kwargs...)
    eigres.converged || (@warn "Eigensolver not converged" iterations=eigres.iterations)

    # Update density from new ψ
    occupation, εF = find_occupation(ham.basis, eigres.λ)
    ρnew = compute_density(ham.basis, eigres.X, occupation)

    (ψ=eigres.X, eigenvalues=eigres.λ, occupation=occupation, εF=εF, ρ=ρnew,
     diagonalization=eigres)
end

"""
Solve the Kohn-Sham equations with a SCF algorithm, starting at ρ.
"""
@timing function self_consistent_field(basis::PlaneWaveBasis;
                                       n_bands=default_n_bands(basis.model),
                                       ρ=guess_density(basis),
                                       ψ=nothing,
                                       tol=1e-6,
                                       maxiter=100,
                                       solver=scf_nlsolve_solver(),
                                       eigensolver=lobpcg_hyper,
                                       n_ep_extra=3,
                                       determine_diagtol=ScfDiagtol(),
                                       mixing=SimpleMixing(),
                                       is_converged=ScfConvergenceEnergy(tol),
                                       callback=ScfDefaultCallback(),
                                       compute_consistent_energies=true,
                                       enforce_symmetry=false,
                                      )
    T = eltype(basis)
    model = basis.model

    # All these variables will get updated by fixpoint_map
    if ψ !== nothing
        @assert length(ψ) == length(basis.kpoints)
        for ik in 1:length(basis.kpoints)
            @assert size(ψ[ik], 2) == n_bands + n_ep_extra
        end
    end
    occupation = nothing
    eigenvalues = nothing
    ρout = ρ
    εF = nothing
    n_iter = 0
    energies = nothing
    ham = nothing
    info = (n_iter=0, ρin=ρ)   # Populate info with initial values

    # We do density mixing in the real representation
    # TODO support other mixing types
    function fixpoint_map(x)
        n_iter += 1
        ρin = from_real(basis, x)

        # Build next Hamiltonian, diagonalize it, get ρout
        if n_iter == 1 # first iteration
            _, ham = energy_hamiltonian(basis, nothing, nothing;
                                        ρ=ρin, eigenvalues=nothing, εF=nothing)
        else
            # Note that ρin is not the density of ψ, and the eigenvalues
            # are not the self-consistent ones, which makes this energy non-variational
            energies, ham = energy_hamiltonian(basis, ψ, occupation;
                                               ρ=ρin, eigenvalues=eigenvalues, εF=εF)
        end

        # Diagonalize `ham` to get the new state
        nextstate = next_density(ham; n_bands=n_bands, ψ=ψ, eigensolver=eigensolver,
                                 miniter=1, tol=determine_diagtol(info),
                                 n_ep_extra=n_ep_extra)
        ψ, eigenvalues, occupation, εF, ρout = nextstate
        enforce_symmetry && (ρout = DFTK.symmetrize(ρout))

        # Compute the energy of the new state
        if compute_consistent_energies
            energies, _ = energy_hamiltonian(basis, ψ, occupation;
                                             ρ=ρout, eigenvalues=eigenvalues, εF=εF)
        end

        # Update info with results gathered so far
        info = (ham=ham, basis=basis, energies=energies, converged=false, ρin=ρin, ρout=ρout,
                eigenvalues=eigenvalues, occupation=occupation, εF=εF, n_iter=n_iter, ψ=ψ,
                diagonalization=nextstate.diagonalization, stage=:iterate)

        # Apply mixing and pass it the full info as kwargs
        ρnext = mix(mixing, basis, ρin, ρout; info...)
        enforce_symmetry && (ρnext = DFTK.symmetrize(ρnext))
        info = merge(info, (ρnext=ρnext, ))

        callback(info)
        is_converged(info) && return x

        ρnext.real
    end

    # Tolerance and maxiter are only dummy here: Convergence is flagged by is_converged
    # inside the fixpoint_map. Also we do not use the return value of fpres but rather the
    # one that got updated by fixpoint_map
    fpres = solver(fixpoint_map, ρout.real, maxiter; tol=eps(T))

    # We do not use the return value of fpres but rather the one that got updated by fixpoint_map
    # ψ is consistent with ρout, so we return that. We also perform
    # a last energy computation to return a correct variational energy
    energies, ham = energy_hamiltonian(basis, ψ, occupation;
                                       ρ=ρout, eigenvalues=eigenvalues, εF=εF)

    # Callback is run one last time with final state to allow callback to clean up
    info = (ham=ham, basis=basis, energies=energies, converged=fpres.converged,
            ρ=ρout, eigenvalues=eigenvalues, occupation=occupation, εF=εF,
            n_iter=n_iter, ψ=ψ, diagonalization=info.diagonalization, stage=:finalize)
    callback(info)
    info
end


text:
using TableWidgets, Interact, StatsPlots, JuliaDB, Blink

function mypipeline(t)
    (t isa Observables.AbstractObservable) || (t = Observable{Any}(t))

    filter = addfilter(t)
    editor = dataeditor(map(DataFrame, filter))
    plotter = dataviewer(editor)

    components = OrderedDict{Symbol, Any}(
        :filter => filter,
        :editor => editor,
        :plotter => plotter
    )
    wdg = Widget(
        components,
        layout = x -> tabulator(components)
    )
end

function mypipeline()
    wdg = filepicker()
    widget(mypipeline∘loadtable, wdg, init = wdg) # initialize the widget as a filepicker, when the filepicker gets used, replace with the output of `mypipeline` called with the loaded table
end

##

w = Window(async = false)
body!(w, mypipeline())


text:
using MAT
using Convex
using ECOS
using Plots
using LinearAlgebra
using LineSearches
using ProgressBars
using WoodburyMatrices
using Printf: @sprintf

basis = matread("../Curvelet_Basis_Construction/G_32_32.mat")
G = basis["G_mat"]
s = basis["scales"]
G = G.*(2 .^ (-2 .*s))
Gn = maximum(sqrt.(sum(G.^2, dims=2)))
G ./= Gn

function d_oed(X::Matrix{T}, B, ϵ, λ; maxiter=2500, ls=MoreThuente(), g_rtol::T = sqrt(eps(T)), g_atol::T = eps(T)) where T <: Number
    #initial vector
    n = size(X, 1)
    μ = ones(n, 1)
    μp = Variable(n, 1)
    gvec = similar(μ)
    s = similar(μ)
    #setup linesearch functions
    #linesearch - this needs to make sure that any interior functions actually have function barriers around global variables

    function Ml(μ, X, ϵ)
        println(Diagonal(μ))
        M = ϵ*I + X'*Diagonal(μ)*X
        Mc = cholesky(Hermitian(M))
        l = μ.*diag(X*(Mc\X'))
        return (M, Mc, l)
    end

    Ml(μ) = Ml(μ, X, ϵ)

    function f(μ, λ)
        M, Mc, l = Ml(μ)
        return -logdet(M) + λ*maximum(l)
    end

    f(μ) = f(μ, λ)

    function gi!(gvec, i, Mc, μ, u, lmind, X, λ, T)
        ld_term = @views -X[i,:]'*(Mc\X[i,:])
        tmp_A = SymWoodbury(Mc, X[i,:], -μ[i])
        tmp_b = @views tmp_A \ X[i,:]
        tmp_1 = (u'*tmp_b)^2
        tmp_d = @views X[i,:]'*tmp_b
        tmp_den = (1+μ[i]*tmp_d)^2
        lev_term_1 =  -μ[lmind]*((1+μ[i]*tmp_d)*tmp_1 - μ[i]*tmp_1*tmp_d) / tmp_den
        lev_term_2 = (i==lmind) ? u'*(Mc\u) : zero(T)
        gvec[i] = ld_term + λ*(lev_term_1 + lev_term_2)
    end

    gi!(gvec, i, Mc, μ, u, lmind) = gi!(gvec, i, Mc, μ, u, lmind, X, λ, T)

    function g!(gvec, μ)
        M, Mc, l = Ml(μ)
        lm, lmind = findmax(l)
        lmind=lmind[1]
        u = @views X[lmind, :]
        for i = 1:n
            gi!(gvec, i, Mc, μ, u, lmind)
        end
    end

    function fg!(gvec, μ)
        M, Mc, l = Ml(μ)
        lm, lmind = findmax(l)
        lmind=lmind[1]
        u = @views X[lmind, :]
        for i = 1:n
            gi!(gvec, i, Mc, μ, u, lmind)
        end
        return -logdet(M) + λ*maximum(l)
    end

    ϕ(α) = f(μ.+α.*s)

    function dϕ(α)
        g!(gvec, μ .+ α.*s)
        return dot(gvec, s)
    end

    function ϕdϕ(α)
        phi = fg!(gvec, μ .+ α.*s)
        dphi = dot(gvec, s)
        return (phi, dphi)
    end

    #get initial gradient and iterate
    g!(gvec, μ)
    fμ = f(μ)
    gnorm = norm(gvec)
    gtol = max(g_rtol*gnorm, g_atol)
    iter = ProgressBar(1:maxiter)
    set_description(iter, string(@sprintf("Starting Up")))

    for i in iter
        if gnorm < gtol
            break
        else
            #set step direction
            s .= -gvec
            dϕ_0 = dot(s, gvec)
            #linesearch gradient
            α, fμ = ls(ϕ, dϕ, ϕdϕ, 1.0, fμ, dϕ_0)
            # α = 0.001
            # fμ = f(μ)
            @. μ = μ + α*s
            #project into constrained space - decrease budget as we go along

            projection = minimize(sumsquares(μp-μ), [0<=μp, μp<=1, sum(μp)<=max(n-i, B)])
            solve!(projection, () -> ECOS.Optimizer(verbose=false))
            μ .= μp.value
            #update gradient
            g!(gvec, μ)
            gnorm = norm(gvec)

            set_description(iter, string(@sprintf("Objective: %.5f", fμ)))
        end
    end

    return μ, round.(μ)
end

resuround, res = d_oed(G, 23, 1.0, 10.0, maxiter=10)

max_stars_repo_name: jbmuir/DAS-Reconstruction

max_issues_repo_name: jbmuir/DAS-Reconstruction

max_forks_repo_name: jbmuir/DAS-Reconstruction


text:
using Compat.Test

@testset "simpler implementations" begin

    # bisection
    xrt = Roots.bisection(sin, 3.0, 4.0)
    @test isapprox(xrt, pi)

    xrt = Roots.bisection(sin, 3.0, 4.0, xatol=1e-3)
    @test abs(sin(xrt)) >= 1e-7  # not to close

    xrt = Roots.bisection(sin, big(3.0), big(4.0))
    @test isapprox(xrt, pi)

    # secant_method
    fpoly(x) = x^5 - x - 1
    xrt = Roots.secant_method(fpoly, 1.0)
    @test abs(fpoly(xrt)) <= 1e-15

    xrt = Roots.secant_method(fpoly, (1, 2))
    @test abs(fpoly(xrt)) <= 1e-14

    # dfree
    fpoly(x) = x^5 - x - 1
    xrt = Roots.dfree(fpoly, 1.0)
    @test abs(fpoly(xrt)) <= 1e-14

end


text:
module Optimize

using LinearAlgebra
using Statistics

using ProgressMeter
using DocStringExtensions

abstract type AbstractOptimizerHistory end;

abstract type AbstractProximalOperator end;

abstract type AbstractOptimizer{T} end;
abstract type AbstractSubspaceOptimizer{T} <: AbstractOptimizer{T} end;

# Overload the norm
LinearAlgebra.norm(x, p, λ) = norm(x[abs.(x) .> λ], p)

# Dispatch on vector
function (opt::AbstractOptimizer{T} where T)(
        X::AbstractArray, A::AbstractArray, Y::AbstractVector{V} where V,
        args...; kwargs...
    )
    Y = reshape(Y, 1, length(Y))
    return opt(X, A, Y, args...; kwargs...)
end

function (opt::AbstractOptimizer{T} where T)(
        X::AbstractArray, A::AbstractArray, Y::Adjoint{V, AbstractVector{V}} where V,
        args...; kwargs...
    )
    Y = reshape(Y, 1, length(Y))
    return opt(X, A, Y, args...; kwargs...)
end

"""
$(SIGNATURES)

Set the threshold(s) of an optimizer to (a) specific value(s).
"""
function set_threshold!(opt::AbstractOptimizer{T}, threshold::T) where T <: AbstractVector
    @assert all(threshold .> zero(T)) "Threshold must be positive definite"
    opt.λ .= threshold
end

function set_threshold!(opt::AbstractOptimizer{T}, threshold::T) where T <: Number
    @assert all(threshold .> zero(T)) "Threshold must be positive definite"
    opt.λ = threshold
end

"""
$(SIGNATURES)

Get the threshold(s) of an optimizer.
"""
get_threshold(opt::AbstractOptimizer) = opt.λ

"""
$(SIGNATURES)

Initialize the optimizer with the least square solution for explicit and `zeros` for implicit optimization.
"""
init(o::AbstractOptimizer, A::AbstractArray, Y::AbstractArray) = A \ Y

init(o::AbstractSubspaceOptimizer, A::AbstractArray, Y::AbstractArray) = zeros(eltype(A), size(A, 2), size(Y, 2))

"""
$(SIGNATURES)

Initialize the optimizer with the least square solution for explicit and `zeros` for implicit optimization in place.
"""
init!(X::AbstractArray, o::AbstractOptimizer, A::AbstractArray, Y::AbstractArray) =  ldiv!(X, qr(A, Val(true)), Y)

include("./utils.jl")

include("./proximals.jl")
export SoftThreshold, HardThreshold,ClippedAbsoluteDeviation

# Remove the trace right now
#include("./history.jl")

include("./stlsq.jl")
include("./admm.jl")
include("./sr3.jl")

#Nullspace for implicit sindy
include("./adm.jl")
include("./implicit.jl")

# Init the progressmeters
# For a general optimizer
default_prg_msg(o::AbstractOptimizer) = summary(o)

function init_progress(opt::AbstractOptimizer, X, A, Y, maxiters, start)
    Progress(
        maxiters, default_prg_msg(opt), start
    )
end

include("./sparseregression.jl")
export sparse_regression!
export init, init!, set_threshold!, get_threshold
export STLSQ, ADMM, SR3
export ImplicitOptimizer, ADM

end


text:
#math_coordinates.jl:
#-------------------------------------------------------------------------------

#==Types
===============================================================================#

#Tag data as being part of a given coordinate system:
struct CoordSystem{ID}; end
const DeviceCoord = CoordSystem{:dev}
const AxisCoord = CoordSystem{:axis}
const NormCoord = CoordSystem{:norm} #Normalized (ex: relative to axis-delimited viewport)
const DataCoord = CoordSystem{:data}

#Never actually used:
#... but could simplify names of mapping functions/mitigate incorrect use.
struct TypedCoord{CT<:CoordSystem}
	v::DReal
end
coord(s::Symbol, v::DReal) = TypedCoord{CoordSystem{s}}(v)

#Annotation coordinates can match data or be normalized to plot bounds (0 -> 1):
#NOTE: not used
const AnnotationCoord = Union{TypedCoord{NormCoord}, TypedCoord{DataCoord}}

#Dispatchable scale:
abstract type AxisScale end #Note: "Scale" too generic a name

struct LinScale{T} <: AxisScale
	tgtmajor::DReal #Targeted number of major grid lines
	tgtminor::Int   #Targeted number of minor grid lines
	(t::Type{LinScale{_T}})(a1, a2) where _T = error("$t not supported")
	(::Type{LinScale{1}})(a1, a2) = new{1}(a1, a2) #Default linear
	(::Type{LinScale{:dB10}})(a1, a2) = new{:dB10}(a1, a2)
	(::Type{LinScale{:dB20}})(a1, a2) = new{:dB20}(a1, a2)
end
LinScale(t=1; tgtmajor=3.5, tgtminor=4) = LinScale{t}(tgtmajor, tgtminor)

struct LogScale{T} <: AxisScale
	#TODO: Configure how many major/minor lines are desired?
	(t::Type{LogScale{_T}})() where _T = error("$t not supported")
	(::Type{LogScale{:e}})() = new{:e}()
	(::Type{LogScale{2}})() = new{2}()
	(::Type{LogScale{10}})() = new{10}()
end
LogScale(t) = LogScale{t}()

function AxisScale(t::Symbol; kwargs...)
	if :ln == t; return LogScale(:e); end
	if :log2 == t; return LogScale(2); end
	if :log10 == t; return LogScale(10); end
	if :log == t; return LogScale(10); end #Alias
	if :lin == t; return LinScale(1; kwargs...); end
	return LinScale(t; kwargs...)
end
AxisScale() = AxisScale(:lin)

#Nonlinear transform type definitions
#-------------------------------------------------------------------------------
abstract type NLTransform{NDIMS} end #Nonlinear transform

#InputXfrm: Maps raw input data to axis coordinates:

#InputXfrm1DSpec: Generates specialized (efficient) code:
struct InputXfrm1DSpec{T} <: NLTransform{1}
end
InputXfrm1DSpec(::LinScale{1}) = InputXfrm1DSpec{:lin}()
InputXfrm1DSpec(::LinScale{:dB10}) = InputXfrm1DSpec{:dB10}()
InputXfrm1DSpec(::LinScale{:dB20}) = InputXfrm1DSpec{:dB20}()
InputXfrm1DSpec(::LogScale{:e}) = InputXfrm1DSpec{:ln}()
InputXfrm1DSpec(::LogScale{2}) = InputXfrm1DSpec{:log2}()
InputXfrm1DSpec(::LogScale{10}) = InputXfrm1DSpec{:log10}()

struct InputXfrm1D <: NLTransform{1} #Generates non-specialized code
	spec::InputXfrm1DSpec
end
InputXfrm1D(s::AxisScale) = InputXfrm1D(InputXfrm1DSpec(s))
struct InputXfrm2D <: NLTransform{2} #Generates non-specialized code
	x::InputXfrm1DSpec #x-transform
	y::InputXfrm1DSpec #y-transform
end
InputXfrm2D(x::InputXfrm1D, y::InputXfrm1D) = InputXfrm2D(x.spec, y.spec)
InputXfrm2D(xs::AxisScale, ys::AxisScale) = InputXfrm2D(InputXfrm1DSpec(xs), InputXfrm1DSpec(ys))

#Offsettable 2D position:
#TODO: Should reloffset exist as is?
#      Alt: v.x/y is *either* "read" or "relative" (use relx::Bool, rely::Bool)
mutable struct Pos2DOffset
	v::Point2D #Position ("Read"able coordinates) - set NaN to use offsets only
	reloffset::Vector2D #Relative offset (Normalized to [0,1] graph bounds)
	offset::Vector2D #Absolute offset (device units)
end

#==Mapping/interpolation functions
===============================================================================#
#Map input data to axis coordinates:
map2axis(::Type{T}, ::InputXfrm1DSpec{:lin}) where T<:Number = (x::T)->DReal(x)
map2axis(::Type{T}, ::InputXfrm1DSpec{:dB10}) where T<:Number = (x::T)->(5*log10(DReal(abs2(x))))
map2axis(::Type{T}, ::InputXfrm1DSpec{:dB20}) where T<:Number = (x::T)->(10*log10(DReal(abs2(x))))
map2axis(::Type{T}, ::InputXfrm1DSpec{:ln}) where T<:Number =
	(x::T)->(x<0 ? DNaN : log(DReal(x)))
map2axis(::Type{T}, ::InputXfrm1DSpec{:log2}) where T<:Number =
	(x::T)->(x<0 ? DNaN : log2(DReal(x)))
map2axis(::Type{T}, ::InputXfrm1DSpec{:log10}) where T<:Number =
	(x::T)->(x<0 ? DNaN : log10(DReal(x)))
#TODO: find a way to show negative values for log10?

#Map axis coord --> user-readable coord (Typically reverse of map2axis):
axis2read(::Type{T}, ::InputXfrm1DSpec) where T<:Number = (x::T)->DReal(x) #NOTE: lin/dB values remain as-is - for readability
axis2read(::Type{T}, ::InputXfrm1DSpec{:ln}) where T<:Number = (x::T)->(exp(x))
axis2read(::Type{T}, ::InputXfrm1DSpec{:log2}) where T<:Number = (x::T)->(2.0^x)
axis2read(::Type{T}, ::InputXfrm1DSpec{:log10}) where T<:Number = (x::T)->(10.0^x)

axis2read(v::T, s::InputXfrm1DSpec) where T<:Number = axis2read(T,s)(v) #One-off conversion

#Map user-readable coord --> axis coord (Typically same as map2axis):
read2axis(::Type{T}, ixf::InputXfrm1DSpec) where T<:Number = map2axis(T, ixf)
#Exception: "Readable" coordinates match axis coordinates when in dB:
read2axis(::Type{T}, ::InputXfrm1DSpec{:dB10}) where T<:Number = (x::DReal)->x
read2axis(::Type{T}, ::InputXfrm1DSpec{:dB20}) where T<:Number = (x::DReal)->x

read2axis(v::T, s::InputXfrm1DSpec) where T<:Number = read2axis(T,s)(v) #One-off conversion

#Map Point2D:
#-------------------------------------------------------------------------------
#TODO: define one-off conversion for map2axis to simplify following?:
map2axis(pt::Point2D, xixf::InputXfrm1DSpec, yixf::InputXfrm1DSpec) =
	Point2D(map2axis(DReal, xixf)(pt.x), map2axis(DReal, yixf)(pt.y))
map2axis(pt::Point2D, ixf::InputXfrm2D) = map2axis(pt, ixf.x, ixf.y)
axis2read(pt::Point2D, xixf::InputXfrm1DSpec, yixf::InputXfrm1DSpec) =
	Point2D(axis2read(pt.x, xixf), axis2read(pt.y, yixf))
axis2read(pt::Point2D, ixf::InputXfrm2D) = axis2read(pt, ixf.x, ixf.y)
read2axis(pt::Point2D, xixf::InputXfrm1DSpec, yixf::InputXfrm1DSpec) =
	Point2D(read2axis(pt.x, xixf), read2axis(pt.y, yixf))
read2axis(pt::Point2D, ixf::InputXfrm2D) = read2axis(pt, ixf.x, ixf.y)

#Map entire data vector:
#-------------------------------------------------------------------------------
function map2axis(d::Vector{T}, xf::InputXfrm1DSpec) where T<:Number
	#Apparently, passing functions as arguments is not efficient in Julia.
	#-->Specializing on InputXfrm1DSpec, hoping to improve efficiency on dmap:
	dmap = map2axis(T, xf)

	result = Array{DReal}(undef, length(d))
	for i in 1:length(d)
		result[i] = dmap(d[i])
	end
	return result
end
map2axis(d::Vector{T}, ::InputXfrm1DSpec{:lin}) where T<:Number = d #Optimization: Linear scale does not need re-scaling

#Extents mapping functions, depending on axis type:
#-------------------------------------------------------------------------------
#Shortcut - hardwire type used in extents:
_extread2axis(t::InputXfrm1DSpec) = read2axis(DReal, t)
_extaxis2read(t::InputXfrm1DSpec) = axis2read(DReal, t)

function read2axis(ext::PExtents1D, xf::InputXfrm1D)
	emap = _extread2axis(xf.spec)
	return PExtents1D(emap(ext.min), emap(ext.max))
end
function read2axis(ext::PExtents2D, xf::InputXfrm2D)
	xmap = _extread2axis(xf.x)
	ymap = _extread2axis(xf.y)
	return PExtents2D(
		xmap(ext.xmin), xmap(ext.xmax),
		ymap(ext.ymin), ymap(ext.ymax)
	)
end

function axis2read(ext::PExtents1D, xf::InputXfrm1D)
	emap = _extaxis2read(xf.spec)
	return PExtents1D(emap(ext.min), emap(ext.max))
end
function axis2read(ext::PExtents2D, xf::InputXfrm2D)
	xmap = _extaxis2read(xf.x)
	ymap = _extaxis2read(xf.y)
	return PExtents2D(
		xmap(ext.xmin), xmap(ext.xmax),
		ymap(ext.ymin), ymap(ext.ymax)
	)
end

#Pos2DOffset mapping functions:
#-------------------------------------------------------------------------------
function map2dev(pos::Pos2DOffset, xf::Transform2D, ixf::InputXfrm2D, graphbb::BoundingBox)
	pt = read2axis(pos.v, ixf)
	pt = map2dev(xf, pt)
	x = pt.x; y = pt.y
	if isnan(x); x = graphbb.xmin; end
	if isnan(y); y = graphbb.ymax; end
	x += pos.reloffset.x * width(graphbb) + pos.offset.x
	y -= pos.reloffset.y * height(graphbb) + pos.offset.y
	return Point2D(x, y)
end

#Last line


text:
"""
    mid_point(M, p, q, x)
    mid_point!(M, y, p, q, x)

Compute the mid point between `p` and `q`. If there is more than one mid point
of (not neccessarily minimizing) geodesics (e.g. on the sphere), the one nearest
to `x` is returned (in place of `y`).
"""
mid_point(M::AbstractManifold, p, q, ::Any) = mid_point(M, p, q)
mid_point!(M::AbstractManifold, y, p, q, ::Any) = mid_point!(M, y, p, q)

function mid_point(M::Circle, p, q, x)
    if distance(M, p, q) ≈ π
        X = 0.5 * log(M, p, q)
        Y = log(M, p, x)
        return exp(M, p, (sign(X) == sign(Y) ? 1 : -1) * X)
    end
    return mid_point(M, p, q)
end

function mid_point(M::Sphere, p, q, x)
    if isapprox(M, p, -q)
        X = log(M, p, x) / distance(M, p, x) * π
    else
        X = log(M, p, q)
    end
    return exp(M, p, 0.5 * X)
end
function mid_point!(M::Sphere, y, p, q, x)
    if isapprox(M, p, -q)
        X = log(M, p, x) / distance(M, p, x) * π
    else
        X = log(M, p, q)
    end
    y .= exp(M, p, 0.5 * X)
    return y
end

function prox_TV2(::Euclidean, λ, pointTuple::Tuple{T,T,T}, p::Int=1) where {T}
    w = @SVector [1.0, -2.0, 1.0]
    x = SVector(pointTuple)
    if p == 1 # Example 3.2 in Bergmann, Laus, Steidl, Weinmann, 2014.
        m = min.(Ref(λ), abs.(x .* w) / (dot(w, w)))
        s = sign.(sum(x .* w))
        return x .- m .* s .* w
    elseif p == 2 # Theorem 3.6 ibd.
        t = λ * sum(x .* w) / (1 + λ * dot(w, w))
        return x .- t .* w
    else
        throw(
            ErrorException(
                "Proximal Map of TV2(Euclidean,λ,pT,p) not implemented for p=$(p) (requires p=1 or 2)",
            ),
        )
    end
end
@doc raw"""
    reflect(M, f, x)
    reflect!(M, q, f, x)

reflect the point `x` from the manifold `M` at the point `f(x)` of the
function ``f: \mathcal M → \mathcal M``, i.e.,

````math
    \operatorname{refl}_f(x) = \operatorname{refl}_{f(x)}(x),
````

Compute the result in `q`.

see also [`reflect`](@ref reflect(M::AbstractManifold, p, x))`(M,p,x)`.
"""

reflect(M::AbstractManifold, pr::Function, x) = reflect(M, pr(x), x)
reflect!(M::AbstractManifold, q, pr::Function, x) = reflect!(M, q, pr(x), x)

@doc raw"""
    reflect(M, p, x)
    reflect!(M, q, p, x)
reflect the point `x` from the manifold `M` at point `p`, i.e.

````math
    \operatorname{refl}_p(x) = \exp_p(-\log_p x).
````

where exp and log denote the exponential and logarithmic map on `M`.
This can also be done in place of `q`.
"""
reflect(M::AbstractManifold, p, x) = exp(M, p, -log(M, p, x))
reflect!(M::AbstractManifold, q, p, x) = exp!(M, q, p, -log(M, p, x))