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

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

function _swaggerinternal_adminUsersAssign(_api::AdminUsersApi, team_id::String, user_id::String; is_restricted=nothing, is_ultra_restricted=nothing, channel_ids=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.assign", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_param(_ctx.form, "is_restricted", is_restricted)  # type Bool
    Swagger.set_param(_ctx.form, "is_ultra_restricted", is_ultra_restricted)  # type Bool
    Swagger.set_param(_ctx.form, "channel_ids", channel_ids)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Add an Enterprise user to a workspace.
Param: team_id::String (required)
Param: user_id::String (required)
Param: is_restricted::Bool
Param: is_ultra_restricted::Bool
Param: channel_ids::String
Return: Dict{String, Any}
"""
function adminUsersAssign(_api::AdminUsersApi, team_id::String, user_id::String; is_restricted=nothing, is_ultra_restricted=nothing, channel_ids=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersAssign(_api, team_id, user_id; is_restricted=is_restricted, is_ultra_restricted=is_ultra_restricted, channel_ids=channel_ids, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersAssign(_api::AdminUsersApi, response_stream::Channel, team_id::String, user_id::String; is_restricted=nothing, is_ultra_restricted=nothing, channel_ids=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersAssign(_api, team_id, user_id; is_restricted=is_restricted, is_ultra_restricted=is_ultra_restricted, channel_ids=channel_ids, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersInvite(_api::AdminUsersApi, team_id::String, email::String, channel_ids::String; custom_message=nothing, real_name=nothing, resend=nothing, is_restricted=nothing, is_ultra_restricted=nothing, guest_expiration_ts=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.invite", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "email", email)  # type String
    Swagger.set_param(_ctx.form, "channel_ids", channel_ids)  # type String
    Swagger.set_param(_ctx.form, "custom_message", custom_message)  # type String
    Swagger.set_param(_ctx.form, "real_name", real_name)  # type String
    Swagger.set_param(_ctx.form, "resend", resend)  # type Bool
    Swagger.set_param(_ctx.form, "is_restricted", is_restricted)  # type Bool
    Swagger.set_param(_ctx.form, "is_ultra_restricted", is_ultra_restricted)  # type Bool
    Swagger.set_param(_ctx.form, "guest_expiration_ts", guest_expiration_ts)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Invite a user to a workspace.
Param: team_id::String (required)
Param: email::String (required)
Param: channel_ids::String (required)
Param: custom_message::String
Param: real_name::String
Param: resend::Bool
Param: is_restricted::Bool
Param: is_ultra_restricted::Bool
Param: guest_expiration_ts::String
Return: Dict{String, Any}
"""
function adminUsersInvite(_api::AdminUsersApi, team_id::String, email::String, channel_ids::String; custom_message=nothing, real_name=nothing, resend=nothing, is_restricted=nothing, is_ultra_restricted=nothing, guest_expiration_ts=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersInvite(_api, team_id, email, channel_ids; custom_message=custom_message, real_name=real_name, resend=resend, is_restricted=is_restricted, is_ultra_restricted=is_ultra_restricted, guest_expiration_ts=guest_expiration_ts, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersInvite(_api::AdminUsersApi, response_stream::Channel, team_id::String, email::String, channel_ids::String; custom_message=nothing, real_name=nothing, resend=nothing, is_restricted=nothing, is_ultra_restricted=nothing, guest_expiration_ts=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersInvite(_api, team_id, email, channel_ids; custom_message=custom_message, real_name=real_name, resend=resend, is_restricted=is_restricted, is_ultra_restricted=is_ultra_restricted, guest_expiration_ts=guest_expiration_ts, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersList(_api::AdminUsersApi, team_id::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "GET", Dict{String, Any}, "/admin.users.list", ["slackAuth"])
    Swagger.set_param(_ctx.query, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.query, "cursor", cursor)  # type String
    Swagger.set_param(_ctx.query, "limit", limit)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

List users on a workspace
Param: team_id::String (required)
Param: cursor::String
Param: limit::Int32
Return: Dict{String, Any}
"""
function adminUsersList(_api::AdminUsersApi, team_id::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersList(_api, team_id; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersList(_api::AdminUsersApi, response_stream::Channel, team_id::String; cursor=nothing, limit=nothing, _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersList(_api, team_id; cursor=cursor, limit=limit, _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersRemove(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.remove", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Remove a user from a workspace.
Param: team_id::String (required)
Param: user_id::String (required)
Return: Dict{String, Any}
"""
function adminUsersRemove(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersRemove(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersRemove(_api::AdminUsersApi, response_stream::Channel, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersRemove(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersSetAdmin(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.setAdmin", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set an existing guest, regular user, or owner to be an admin user.
Param: team_id::String (required)
Param: user_id::String (required)
Return: Dict{String, Any}
"""
function adminUsersSetAdmin(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetAdmin(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersSetAdmin(_api::AdminUsersApi, response_stream::Channel, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetAdmin(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersSetExpiration(_api::AdminUsersApi, team_id::String, user_id::String, expiration_ts::Int32; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.setExpiration", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_param(_ctx.form, "expiration_ts", expiration_ts)  # type Int32
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set an expiration for a guest user
Param: team_id::String (required)
Param: user_id::String (required)
Param: expiration_ts::Int32 (required)
Return: Dict{String, Any}
"""
function adminUsersSetExpiration(_api::AdminUsersApi, team_id::String, user_id::String, expiration_ts::Int32; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetExpiration(_api, team_id, user_id, expiration_ts; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersSetExpiration(_api::AdminUsersApi, response_stream::Channel, team_id::String, user_id::String, expiration_ts::Int32; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetExpiration(_api, team_id, user_id, expiration_ts; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersSetOwner(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.setOwner", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set an existing guest, regular user, or admin user to be a workspace owner.
Param: team_id::String (required)
Param: user_id::String (required)
Return: Dict{String, Any}
"""
function adminUsersSetOwner(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetOwner(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersSetOwner(_api::AdminUsersApi, response_stream::Channel, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetOwner(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

function _swaggerinternal_adminUsersSetRegular(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = Swagger.Ctx(_api.client, "POST", Dict{String, Any}, "/admin.users.setRegular", ["slackAuth"])
    Swagger.set_param(_ctx.form, "team_id", team_id)  # type String
    Swagger.set_param(_ctx.form, "user_id", user_id)  # type String
    Swagger.set_header_accept(_ctx, ["application/json"])
    Swagger.set_header_content_type(_ctx, (_mediaType === nothing) ? ["application/x-www-form-urlencoded", "application/json"] : [_mediaType])
    return _ctx
end

"""

Set an existing guest user, admin user, or owner to be a regular user.
Param: team_id::String (required)
Param: user_id::String (required)
Return: Dict{String, Any}
"""
function adminUsersSetRegular(_api::AdminUsersApi, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetRegular(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx)
end

function adminUsersSetRegular(_api::AdminUsersApi, response_stream::Channel, team_id::String, user_id::String; _mediaType=nothing)
    _ctx = _swaggerinternal_adminUsersSetRegular(_api, team_id, user_id; _mediaType=_mediaType)
    Swagger.exec(_ctx, response_stream)
end

export adminUsersAssign, adminUsersInvite, adminUsersList, adminUsersRemove, adminUsersSetAdmin, adminUsersSetExpiration, adminUsersSetOwner, adminUsersSetRegular


text:
# readFastq_test.jl
@testset "readFastq test" begin
    record = readFastq("test_data/test.fastq")
    @test record.ID[1] == "testseq1"
    @test record.ID[2] == "testseq2"
    @test record.ID[3] == "testseq3"
    @test record.sequence[1] == "GGACTACCGGGGTATCTAATCCTGTTTGCTCCCCACGCTTTCGTACCTCAGTGTCAGTATTAGTCCAGGAAGCCGCCTTCGCCACTGGTGTTCCTCCGGATATCTACGCATTTCACCGCTACACCCGGAATTCCACTTCCCTCTACTATACTCTAGCCAGACAGTATCGAATGCAGTTCCCAGGTTGAGCCCAGGGCTTTCACATCCGACTAATCAAGCCACCTACGTACGCTTTACGCCCAGTAATTCCGATTAACGCTTGCACCCTCCGTATTACCGCGGCGGCTGGCACCTGATCTC"
    @test record.sequence[2] == "GATGGTAGGTGCCAGCAGCCGCGGTAATACGAAAGGTGCGAGCGTTAATCGGAATTACTGGGCGTAAAGCGCGCGTAGGTGGTGTGTTAAGTCGGATGTGAAAGCCCTGGGCTCAACCTGGGAACTGCATCCGATACTGGCACGCTAGAGTGCAGTAGAGGGAGGTGGAATTTCCGGTGTAGCGGTGAAATGCGTAGAGATCGGAAGGAACACCAGTGGCGAAGGCGGCCTCCTGGACTGACACTGACACTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCCGGTAGTCC"
    @test record.sequence[3] == "GATGCTTGGTGCCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGATTTACTGGGCGTAAAGCGCACGTAGGCGGACTATTAAGTCAGGGGTGAAATCCCGGGGCTCAACCCCGGAACTGCCTTTGATACTGGTAGTCTTGAGTTCGAGAGAGGTGAGTGGAATTCCGAGTGTAGAGGTGAAATTCGTAGATATTCGGAGGAACACCAGTGGCGAAGGCGGCTCACTGGCTCGATACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGATACCCCCGTAGTCC"
    @test record.quality[1] == "CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGEGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGFGGGFGGFGBGFGGGFFFC;>FFFFEGFGFGFFGFFF"
    @test record.quality[2] == "CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGGGGGFEGGGGGGFEGGFGFGFFGAFGFGGGFFGFGFF"
    @test record.quality[3] == "CCCCCGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGGFGGGGEGGGFGGGBCDGGGGG>DGDGFFGFGGGFFGG?FFFFFCFC@G?FB>FB"
end

ext: jl


text:
"""
Authors: Various
Compiled by: Eric Ford (see each function for credits)
Created: August 2020
Contact: https://github.com/eford/
"""

# util/physics.jl included separatesly from RvSpectMLBase since may be used early
#include("util/physics.jl")

include("find_line.jl")
# intended for internal use, so functions not exported

include("linelist.jl")
export assign_lines_to_orders
export add_line_boundaries_to_line_list, expand_line_boundaries_in_line_list

include("chunks.jl")
export make_orders_into_chunks, make_grid_for_chunk,  extract_chunk_list_timeseries_for_order
export make_chunk_list_timeseries_from_λ_ranges
export make_chunk_list_around_lines, make_chunk_list_timeseries_around_lines, make_chunk_list_tellric_free, make_chunk_list_timeseries_tellric_free
export make_order_list_timeseries
export find_orders_with_line, find_orders_in_range
export findall_line, find_line_best
export filter_bad_chunks, find_pixels_for_line_in_chunk, find_pixels_for_line_in_chunklist
export calc_normalization, calc_normalization_var_weighted
export set_rv_est!

include("spectra.jl")
export apply_doppler_boost!, normalize_spectrum!, normalize_spectra!
export get_λ_range, get_inst
export get_order_info
export make_vec_metadata_from_spectral_timeseries, discard_large_metadata
#export calc_snr

include("cs_algs.jl")
# Don't export to avoid risk of conflicts
# export searchsortednearest

include("files.jl")
export read_data_paths, code_to_include_param_jl
#export make_manifest

include("pipeline.jl")
using .Pipeline
export PipelinePlan
export make_plot,  save_plot, save_data, need_to, has_cache, read_cache, set_cache!  # Query pipeline
export need_to!, dont_need_to!, reset_all_needs!             # Write to pipeline
export make_plot!, dont_make_plot!, make_all_plots!,  make_no_plots!

max_stars_repo_name: RvSpectML/RvSpectMLTemplate

max_issues_repo_name: RvSpectML/RvSpectMLTemplate

max_forks_repo_name: RvSpectML/RvSpectMLTemplate


text:
#! /usr/bin/julia

# Rosetta Code, Knapsack problem/Continuous

immutable KPCSupply{S<:String, T<:Real}
    item::S
    weight::T
    value::T
    uvalue::T
end
Base.isless(a::KPCSupply, b::KPCSupply) = a.uvalue<b.uvalue

function KPCSupply{S<:String, T<:Real}(item::S, weight::T, value::T)
    KPCSupply(item, weight, value, value/weight)
end

function solve{S<:String, T<:Real}(store::Array{KPCSupply{S,T},1},
                                   capacity::T)
    ksack = KPCSupply{S,T}[]
    kweight = zero(T)
    for s in sort(store, rev=true)
        if kweight + s.weight <= capacity
            kweight += s.weight
            push!(ksack, s)
        else
            w = capacity-kweight
            v = w*s.uvalue
            push!(ksack, KPCSupply(s.item, w, v, s.uvalue))
            break
        end
    end
    return ksack
end

store = [KPCSupply("beef", 38//10, 36//1),
         KPCSupply("pork", 54//10, 43//1),
         KPCSupply("ham", 36//10, 90//1),
         KPCSupply("greaves", 24//10, 45//1),
         KPCSupply("flitch", 4//1, 30//1),
         KPCSupply("brawn", 25//10, 56//1),
         KPCSupply("welt", 37//10, 67//1),
         KPCSupply("salami", 3//1, 95//1),
         KPCSupply("sausage", 59//10, 98//1)]

sack = solve(store, 15//1)

println("The store contains:")
println("        Item Weight  Unit Price")
for s in store
    println(@sprintf("%12s %4.1f    %6.2f",
                     s.item, float(s.weight), float(s.uvalue)))
end

println()
println("The thief should take:")
println("        Item Weight  Value")
for s in sack
    println(@sprintf("%12s %4.1f    %6.2f",
                     s.item, float(s.weight), float(s.value)))
end
w = mapreduce(x->x.weight, +, sack)
v = mapreduce(x->x.value, +, sack)
println(@sprintf("%12s %4.1f    %6.2f",
                 "Total", float(w), float(v)))

max_stars_repo_name: MichaeLeroy/rosetta-code

max_issues_repo_name: MichaeLeroy/rosetta-code

max_forks_repo_name: MichaeLeroy/rosetta-code


text:
module Tricks

using Base: rewrap_unionall, unwrap_unionall, uncompressed_ast
using Base: CodeInfo

export static_hasmethod, static_methods, static_method_count, compat_hasmethod,
        static_fieldnames, static_fieldcount, static_fieldtypes

# This is used to create the CodeInfo returned by static_hasmethod.
_hasmethod_false(@nospecialize(f), @nospecialize(t)) = false
_hasmethod_true(@nospecialize(f), @nospecialize(t)) = true

"""
    static_hasmethod(f, type_tuple::Type{<:Tuple)

Like `hasmethod` but runs at compile-time (and does not accept a worldage argument).
"""
@generated function static_hasmethod(@nospecialize(f), @nospecialize(t::Type{T}),) where {T<:Tuple}
    # The signature type:
    method_insts = _method_instances(f, T)

    ftype = Tuple{f, fieldtypes(T)...}
    covering_method_insts = [mi for mi in method_insts if ftype <: mi.def.sig]

    method_doesnot_exist = isempty(covering_method_insts)
    ret_func = method_doesnot_exist ? _hasmethod_false : _hasmethod_true
    ci_orig = uncompressed_ast(typeof(ret_func).name.mt.defs.func)
    ci = ccall(:jl_copy_code_info, Ref{CodeInfo}, (Any,), ci_orig)

    # Now we add the edges so if a method is defined this recompiles
    if method_doesnot_exist
        # No method so attach to method table
        mt = f.name.mt
        typ = rewrap_unionall(Tuple{f, unwrap_unionall(T).parameters...}, T)
        ci.edges = Core.Compiler.vect(mt, typ)
    else  # method exists, attach edges to all instances
        ci.edges = covering_method_insts
    end
    return ci
end

function create_codeinfo_with_returnvalue(argnames, spnames, sp, value)
    expr = Expr(:lambda,
        argnames,
        Expr(Symbol("scope-block"),
            Expr(:block,
                Expr(:return, value),
            )
        )
    )
    if spnames !== nothing
        expr = Expr(Symbol("with-static-parameters"), expr, spnames...)
    end
    ci = ccall(:jl_expand, Any, (Any, Any), expr, @__MODULE__)
    return ci
end

"""
    static_methods(f, [type_tuple::Type{<:Tuple])
    static_methods(@nospecialize(f)) = _static_methods(Main, f, Tuple{Vararg{Any}})

Like `methods` but runs at compile-time (and does not accept a worldage argument).
"""
static_methods(@nospecialize(f)) = static_methods(f, Tuple{Vararg{Any}})
@generated function static_methods(@nospecialize(f) , @nospecialize(_T::Type{T})) where {T <: Tuple}
    list_of_methods = _methods(f, T)
    ci = create_codeinfo_with_returnvalue([Symbol("#self#"), :f, :_T], [:T], (:T,), :($list_of_methods))

    # Now we add the edges so if a method is defined this recompiles
    ci.edges = _method_table_all_edges_all_methods(f, T)
    return ci
end

function _method_table_all_edges_all_methods(f, T)
    mt = f.name.mt

    # We add an edge to the MethodTable itself so that when any new methods
    # are defined, it recompiles the function.
    mt_edges = Core.Compiler.vect(mt, Tuple{Vararg{Any}})

    # We want to add an edge to _every existing method instance_, so that
    # the deletion of any one of them will trigger recompilation of the function.
    world = typemax(UInt)
    method_insts = _method_instances(f, T)
    covering_method_insts = method_insts

    return vcat(mt_edges, covering_method_insts)
end

"""
    static_method_count(f, [type_tuple::Type{<:Tuple])
    static_method_count(@nospecialize(f)) = _static_methods(Main, f, Tuple{Vararg{Any}})
Returns `length(methods(f, tt))` but runs at compile-time (and does not accept a worldage argument).
"""
static_method_count(@nospecialize(f)) = static_method_count(f, Tuple{Vararg{Any}})
@generated function static_method_count(@nospecialize(f) , @nospecialize(_T::Type{T})) where {T <: Tuple}
    method_count = length(_methods(f, T))
    ci = create_codeinfo_with_returnvalue([Symbol("#self#"), :f, :_T], [:T], (:T,), :($method_count))

    # Now we add the edges so if a method is defined this recompiles
    ci.edges = _method_table_all_edges_all_methods(f, T)
    return ci
end

@static if VERSION < v"1.3"
    const compat_hasmethod = hasmethod
else
    const compat_hasmethod = static_hasmethod
end

Base.@pure static_fieldnames(t::Type) = Base.fieldnames(t)
Base.@pure static_fieldtypes(t::Type) = Base.fieldtypes(t)
Base.@pure static_fieldcount(t::Type) = Base.fieldcount(t)

# The below methods are copied and adapted from Julia Base:
# - https://github.com/JuliaLang/julia/blob/4931faa34a8a1c98b39fb52ed4eb277729120128/base/reflection.jl#L952-L966
# - https://github.com/JuliaLang/julia/blob/4931faa34a8a1c98b39fb52ed4eb277729120128/base/reflection.jl#L893-L896
# - https://github.com/JuliaLang/julia/blob/4931faa34a8a1c98b39fb52ed4eb277729120128/base/reflection.jl#L1047-L1055
# Like Base.methods, but accepts f as a _type_ instead of an instance.
function _methods(@nospecialize(f_type), @nospecialize(t_type),
                 mod::Union{Tuple{Module},AbstractArray{Module},Nothing}=nothing)
    tt = _combine_signature_type(f_type, t_type)
    lim, world = -1, typemax(UInt)
    mft = Core.Compiler._methods_by_ftype(tt, lim, world)
    ms = Base.Method[_get_method(m) for m in mft if (mod === nothing || m.method.module ∈ mod)]
    return Base.MethodList(ms, f_type.name.mt)
end
# Like Core.Compiler.method_instances, but accepts f as a _type_ instead of an instance.
function _method_instances(@nospecialize(f_type), @nospecialize(t_type))
    tt = _combine_signature_type(f_type, t_type)
    lim, world = -1, typemax(UInt)
    mft = Core.Compiler._methods_by_ftype(tt, lim, world)
    return Core.MethodInstance[_specialize_method(match) for match in mft]
end
# Like Base.signature_type, but starts with a type for f_type already.
function _combine_signature_type(@nospecialize(f_type::Type), @nospecialize(args::Type))
    u = unwrap_unionall(args)
    return rewrap_unionall(Tuple{f_type, u.parameters...}, args)
end
# MethodMatch is only defined in v1.6+, so the values returned from _methods_by_ftype need
# a bit of massaging here.
if VERSION < v"1.6"
    # _methods_by_ftype returns a triple
    _get_method((mtypes, msp, method)) = method
    # Core.Compiler.specialize_method(::MethodMatch) is only defined on v1.6+
    function _specialize_method(method_data)
        mtypes, msp, m = method_data
        instance = ccall(:jl_specializations_get_linfo, Ref{Core.MethodInstance}, (Any, Any, Any), m, mtypes, msp)
        return instance
    end
else
    # _methods_by_ftype returns a MethodMatch
    _get_method(method_match) = method_match.method
    _specialize_method = Core.Compiler.specialize_method
end

end  # module


text:
function golomb(n, L, ::Val{:raw})
    m = model(; kind=:golomb)

    # Add variables
    d = domain(0:L)
    foreach(_ -> variable!(m, d), 1:n)

    # Extract error function from usual_constraint
    e1 = (x; param=nothing, dom_size=n) -> error_f(
        usual_constraints[:all_different])(x; param, dom_size
    )
    e2 = (x; param=nothing, dom_size=n) -> error_f(
        usual_constraints[:all_equal_param])(x; param, dom_size
    )
    e3 = (x; param=nothing, dom_size=n) -> error_f(
        usual_constraints[:dist_different])(x; param, dom_size
    )

    # # Add constraints
    constraint!(m, e1, 1:n)
    constraint!(m, x -> e2(x; param=0), 1:1)
    for i in 1:(n - 1), j in (i + 1):n, k in i:(n - 1), l in (k + 1):n
        (i, j) < (k, l) || continue
        constraint!(m, e3, [i, j, k, l])
    end

    # Add objective
    objective!(m, o_dist_extrema)

    return m
end

function golomb(n, L, ::Val{:JuMP})
    m = JuMP.Model(CBLS.Optimizer)

    @variable(m, 0 ≤ X[1:n] ≤ L, Int)

    @constraint(m, X in AllDifferent()) # different marks
    @constraint(m, X in Ordered()) # for output convenience, keep them ordered

    # No two pairs have the same length
    for i in 1:(n - 1), j in (i + 1):n, k in i:(n - 1), l in (k + 1):n
        (i, j) < (k, l) || continue
        @constraint(m, [X[i], X[j], X[k], X[l]] in DistDifferent())
    end

    # Add objective
    @objective(m, Min, ScalarFunction(maximum))

    return m, X
end

"""
    golomb(n, L=n²)

Model the Golomb problem of `n` marks on the ruler `0:L`. The `modeler` argument accepts :raw, and :JuMP (default), which refer respectively to the solver internal model, the MathOptInterface model, and the JuMP model.
"""
golomb(n, L=n^2; modeler = :JuMP) = golomb(n, L, Val(modeler))


text:
if isdefined(Base, :datatype_name)
    datatype_name(x) = Base.datatype_name(x)
else
    datatype_name(t) = t.name.name
end

if VERSION < v"0.6.0-dev"
    include("fast-dates.jl")

    include("Str.jl")
else
    const Str = String
    import Base.Dates: CONVERSION_SPECIFIERS, TimeType, DatePart, tryparsenext, character_codes, genvar, CONVERSION_TRANSLATIONS, CONVERSION_DEFAULTS, _directives, DateLocale

    include("date-tryparse-internal.jl")
    const ISODateFormat = Base.Dates.ISODateFormat
    const ISODateTimeFormat = Base.Dates.ISODateTimeFormat
    const RFC1123Format = Base.Dates.RFC1123Format
end


text:
using Dates
# using GitCommand
using GitHub
using JSON
using Pkg
using Printf
using RegistryCI
using Test
using TimeZones

const AutoMerge = RegistryCI.AutoMerge

const path = joinpath(DEPOT_PATH[1], "registries", "General")
RegistryCI.test(path)


text:
# ---
# title: 48. Rotate Image
# id: problem48
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Array
# link: <https://leetcode.com/problems/rotate-image/description/>
# hidden: true
# ---
#
# You are given an _n_ x _n_ 2D `matrix` representing an image, rotate the image
# by 90 degrees (clockwise).
#
# You have to rotate the image [**in-place**](https://en.wikipedia.org/wiki/In-
# place_algorithm), which means you have to modify the input 2D matrix directly.
# **DO NOT** allocate another 2D matrix and do the rotation.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/mat1.jpg)
#
#
#
#     Input: matrix = [[1,2,3],[4,5,6],[7,8,9]]
#     Output: [[7,4,1],[8,5,2],[9,6,3]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2020/08/28/mat2.jpg)
#
#
#
#     Input: matrix = [[5,1,9,11],[2,4,8,10],[13,3,6,7],[15,14,12,16]]
#     Output: [[15,13,2,5],[14,3,4,1],[12,6,8,9],[16,7,10,11]]
#
#
# **Example 3:**
#
#
#
#     Input: matrix = [[1]]
#     Output: [[1]]
#
#
# **Example 4:**
#
#
#
#     Input: matrix = [[1,2],[3,4]]
#     Output: [[3,1],[4,2]]
#
#
#
#
# **Constraints:**
#
#   * `matrix.length == n`
#   * `matrix[i].length == n`
#   * `1 <= n <= 20`
#   * `-1000 <= matrix[i][j] <= 1000`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Test
import GDAL
import ArchGDAL; const AG = ArchGDAL

AG.registerdrivers() do
    @testset "Raster Tutorial" begin
        AG.read("data/utmsmall.tif") do dataset
            driver = AG.getdriver("GTiff")
            @test AG.shortname(driver) == "GTiff"
            @test AG.longname(driver) == "GeoTIFF"
            @test AG.width(dataset) == 100
            @test AG.height(dataset) == 100
            @test AG.nraster(dataset) == 1

            nad27_prefix = "PROJCS[\"NAD27 / UTM zone 11N\",GEOGCS[\"NAD27\",DATUM[\"North_American_Datum_1927\","
            @test startswith(AG.getproj(dataset), nad27_prefix) == true
            @test AG.getgeotransform(dataset) ≈
                    [440720.0,60.0,0.0,3.75132e6,0.0,-60.0]

            band = AG.getband(dataset, 1)
            @test AG.getblocksize(band) ≈ [100, 81]
            @test AG.getdatatype(band) == UInt8
            @test AG.getname(AG.getcolorinterp(band)) == "Gray"

            @test AG.getminimum(band) ≈ 0.0
            @test AG.getmaximum(band) ≈ 255.0

            @test AG.noverview(band) == 0

            # Reading Raster Data
            @test AG.width(band) == 100
            data = map(Float32, AG.read(dataset, 1))
            @test data[:,1] ≈ Float32[107.0f0,123.0f0,132.0f0,115.0f0,
                    132.0f0,132.0f0,140.0f0,132.0f0,132.0f0,132.0f0,107.0f0,
                    132.0f0,107.0f0,132.0f0,132.0f0,107.0f0,123.0f0,115.0f0,
                    156.0f0,148.0f0,107.0f0,132.0f0,107.0f0,115.0f0,99.0f0,
                    123.0f0,99.0f0,74.0f0,115.0f0,82.0f0,115.0f0,115.0f0,
                    107.0f0,123.0f0,123.0f0,99.0f0,123.0f0,123.0f0,115.0f0,
                    115.0f0,107.0f0,90.0f0,99.0f0,107.0f0,107.0f0,99.0f0,
                    123.0f0,107.0f0,140.0f0,123.0f0,123.0f0,115.0f0,99.0f0,
                    132.0f0,123.0f0,115.0f0,115.0f0,123.0f0,132.0f0,115.0f0,
                    123.0f0,132.0f0,214.0f0,156.0f0,165.0f0,148.0f0,115.0f0,
                    148.0f0,156.0f0,148.0f0,140.0f0,165.0f0,156.0f0,197.0f0,
                    156.0f0,197.0f0,140.0f0,173.0f0,156.0f0,165.0f0,148.0f0,
                    156.0f0,206.0f0,214.0f0,181.0f0,206.0f0,173.0f0,222.0f0,
                    206.0f0,255.0f0,214.0f0,173.0f0,214.0f0,255.0f0,214.0f0,
                    247.0f0,255.0f0,230.0f0,206.0f0,197.0f0]

            @test AG.metadatadomainlist(dataset) == ["IMAGE_STRUCTURE", "", "DERIVED_SUBDATASETS"]
            @test AG.metadata(dataset) == ["AREA_OR_POINT=Area"]
            @test AG.metadata(dataset, domain = "IMAGE_STRUCTURE") == ["INTERLEAVE=BAND"]
            @test AG.metadata(dataset, domain = "") == ["AREA_OR_POINT=Area"]
            @test AG.metadata(dataset, domain = "DERIVED_SUBDATASETS") == [
                "DERIVED_SUBDATASET_1_NAME=DERIVED_SUBDATASET:LOGAMPLITUDE:data/utmsmall.tif",
                "DERIVED_SUBDATASET_1_DESC=log10 of amplitude of input bands from data/utmsmall.tif"
            ]
        end

        # Techniques for Creating Files
        #@test GDAL.getmetadataitem(driver, "DCAP_CREATE", "") == "YES"
        #@test GDAL.getmetadataitem(driver, "DCAP_CREATECOPY", "") == "YES"

        AG.read("data/utmsmall.tif") do ds_src
            AG.write(ds_src, "tmp/utmsmall.tif")
        end
        rm("tmp/utmsmall.tif")
    end

    @testset "Vector Tutorial" begin
        AG.read("data/point.geojson") do dataset
            @test AG.nlayer(dataset) == 1
            layer = AG.getlayer(dataset, 0)
            @test (AG.getname(layer) in ["point", "OGRGeoJSON"]) == true
            # layerbyname = AG.getlayer(dataset, "point")
            # @test layerbyname.ptr == layer.ptr
            AG.resetreading!(layer)

            featuredefn = AG.getlayerdefn(layer)
            @test AG.nfield(featuredefn) == 2
            fielddefn = AG.getfielddefn(featuredefn, 0)
            @test AG.gettype(fielddefn) == GDAL.OFTReal
            fielddefn = AG.getfielddefn(featuredefn, 1)
            @test AG.gettype(fielddefn) == GDAL.OFTString

            AG.nextfeature(layer) do feature
                @test AG.asdouble(feature, 0) ≈ 2.0
                @test AG.asstring(feature, 1) == "point-a"
            end
            AG.nextfeature(layer) do feature # second feature
                @test AG.asdouble(feature, 0) ≈ 3.0
                @test AG.asstring(feature, 1) == "point-b"

                geometry = AG.getgeom(feature)
                @test AG.getgeomname(geometry) == "POINT"
                @test AG.getgeomtype(geometry) == GDAL.wkbPoint
                @test AG.nfield(featuredefn) == 2
                @test AG.getx(geometry, 0) ≈ 100.2785
                @test AG.gety(geometry, 0) ≈ 0.0893
                @test AG.getpoint(geometry, 0) == (100.2785,0.0893,0.0)
            end
        end

        @testset "Approach 1" begin
            pointshapefile = "tmp/point_out"
            AG.create("$pointshapefile.shp", "ESRI Shapefile") do dataset
                layer = AG.createlayer(dataset, "point_out", geom=GDAL.wkbPoint)
                AG.createfielddefn("Name", GDAL.OFTString) do fielddefn
                    AG.setwidth!(fielddefn, 32)
                    AG.createfield!(layer, fielddefn, true)
                end
                featuredefn = AG.getlayerdefn(layer)
                @test AG.getname(featuredefn) == "point_out"
                @test AG.nfeature(layer) == 0
                AG.createfeature(featuredefn) do feature
                    AG.setfield!(feature, AG.getfieldindex(feature, "Name"), "myname")
                    AG.setgeomdirectly!(feature, AG.unsafe_createpoint(100.123, 0.123))
                    AG.createfeature!(layer, feature)
                end
                @test AG.nfeature(layer) == 1
            end

            rm("$pointshapefile.dbf")
            rm("$pointshapefile.shp")
            rm("$pointshapefile.shx")
        end

        @testset "Approach 2" begin
            pointshapefile = "tmp/point_out"
            AG.create("$pointshapefile.shp", "ESRI Shapefile") do dataset
                layer = AG.createlayer(dataset, "point_out", geom=GDAL.wkbPoint)
                AG.createfielddefn("Name", GDAL.OFTString) do fielddefn
                    AG.setwidth!(fielddefn, 32)
                    AG.createfield!(layer, fielddefn, true)
                end
                featuredefn = AG.getlayerdefn(layer)
                @test AG.getname(featuredefn) == "point_out"
                @test AG.nfeature(layer) == 0
                AG.createfeature(layer) do feature
                    AG.setfield!(feature, AG.getfieldindex(feature, "Name"), "myname")
                    AG.setgeomdirectly!(feature, AG.unsafe_createpoint(100.123, 0.123))
                end
                @test AG.nfeature(layer) == 1
            end

            rm("$pointshapefile.dbf")
            rm("$pointshapefile.shp")
            rm("$pointshapefile.shx")
        end
    end
end


text:
using SolitonBEC
using Printf
using Plots
using Elliptic

### SIMULATE GENERATION OF SOLITON IN Li-7 CONDENSATE

# --------- Numerics ---------
num = Numerics(
  50e-3, #T
  5e-5, #dt
  10e-6 * E(2 * pi, 0.9), #S
  10e-9 * E(2 * pi, 0.9), #ds
)

# --------- Simulation ---------
khaykovich_gpe = Simulation(
  "GPE",
  "barrier",
)

#--------- Potential ---------
ellipse = Potential(
  "ellipse", # type
  5e-6, #width
  15e-6, #position
  5e-6, #energy
  0.9, #ϵ
  10e-6 #a
)

barrier = Potential(
  "barrier", #type
  5e-6, #width
  15e-6, #position
  5e-6, #energy
  0.9, #ϵ
  10e-6 #a
)

# --------- Apparata ---------
std_apparatus = Apparatus(
  6.941 * 1.660539e-27, #m (conversion AMU -> kg)
  -0.21e-9, # as
  2 * pi * 710, # ω_perp
  4e3, #N
  1.77e-11, #γ
  2 * pi * 4,# ω_z
)

# --------- InitialStates ---------
InitialState1 = InitialState(
  "sech", #type
  1.7e-6, # width
  1.5e6 #v0
)

## Configurations
configs = [
  (num, khaykovich_gpe, barrier, std_apparatus, InitialState1),
  #(num, khaykovich_gpe, barrier, std_apparatus, InitialState1),
]
mem_limit = 15000000000 #byte
cnt = 1

angle = LinRange(0, 2 * pi, 100)

for (num, sim, pot, app, state) in configs

  global cnt
  @printf("\n------Running simulation # %3i \n", cnt)
  estimate = mem_estimate(num)

  if estimate < mem_limit
    @printf("Estimated memory consumption: %4.1f MiB\n", estimate / (1024^2))
    @time run_dynamics(num, sim, app, pot, state)
  else
    @printf("Estimated memory consumption (%4.1f MiB) exceed maximum!\n", estimate / (1024^2))
  end

  cnt += 1

end

max_stars_repo_name: lorenzifrancesco/soliton-BEC

max_issues_repo_name: lorenzifrancesco/soliton-BEC

max_forks_repo_name: lorenzifrancesco/soliton-BEC


text:
#----------------- Get the A matrix
function getDivGrad(n1,n2,n3)

        # the Divergence
        D1 = kron(speye(n3),kron(speye(n2),ddx(n1)))
        D2 = kron(speye(n3),kron(ddx(n2),speye(n1)))
        D3 = kron(ddx(n3),kron(speye(n2),speye(n1)))
        # DIV from faces to cell-centers
        Div = [D1 D2 D3]

        return Div*Div';
end
#----------------- 1D finite difference on staggered grid
function ddx(n)
# generate 1D derivatives
    return d = spdiagm((-ones(n),ones(n)),[0,1],n,n+1)

end


text:
module MultiFESpacesTests

using Test
using Gridap

model = CartesianDiscreteModel(domain=(0.0,1.0,0.0,1.0), partition=(4,4))

order = 1
diritag = "boundary"
fespace = H1ConformingFESpace(Float64,model,order,diritag)

ufun1(x) = x[1] + x[2]
U1 = TrialFESpace(fespace,ufun1)

ufun2(x) = x[1] + x[2]
U2 = TrialFESpace(fespace,ufun2)

U = [U1,U2]

@test string(U) == "Vector of FESpace objects with 2 fields"

U = MultiFESpace([U1,U2])

@test string(U) == "MultiFESpace object with 2 fields"

@test length(U) == 2
@test U[1] === U1
@test U[2] === U2

V1, state = iterate(U)
@test V1 === U1

V2, state = iterate(U,state)
@test V2 === U2

@test num_free_dofs(U) == num_free_dofs(U1) + num_free_dofs(U2)

a(v,u) = varinner(v,u)

bfun(x) = x[2]
b(v) = varinner(v,CellField(trian,bfun))

trian = Triangulation(model)
quad = CellQuadrature(trian,degree=2)

v1 = CellBasis(V1)
v2 = CellBasis(V2)

u1 = CellBasis(U1)
u2 = CellBasis(U2)

mat11 = integrate(a(v1,u1),trian,quad)
mat12 = integrate(a(v1,u2),trian,quad)
mat22 = integrate(a(v2,u2),trian,quad)

mat = MultiCellMatrix([mat11,mat12,mat22],[(1,1),(1,2),(2,2)])

vec1 = integrate(b(v1),trian,quad)
vec2 = integrate(b(v2),trian,quad)

vec = MultiCellVector([vec1,vec2],[(1,),(2,)])

cellids = IdentityCellNumber(Int,length(vec1))
_vec = apply_constraints(U,vec,cellids)
@test _vec.blocks[1] === vec1
@test _vec.blocks[2] === vec2

_mat = apply_constraints_rows(U,mat,cellids)
@test _mat.blocks[1] === mat11
@test _mat.blocks[2] === mat12
@test _mat.blocks[3] === mat22

_mat = apply_constraints_cols(U,mat,cellids)
@test _mat.blocks[1] === mat11
@test _mat.blocks[2] === mat12
@test _mat.blocks[3] === mat22

dofs = celldofids(U)
@test dofs.blocks[1] === celldofids(U1)
@test dofs.blocks[2] === celldofids(U2)

n = num_free_dofs(U)
n1 = num_free_dofs(U[1])
n2 = num_free_dofs(U[2])
v = rand(n)

f1 = restrict_to_field(U,v,1)
@test f1 == v[1:n1]

f2 = restrict_to_field(U,v,2)
@test f2 == v[1+n1:n]

end # module MultiFESpacesTests


text:
module TestComponents

using Mimi
using Test

import Mimi:
    compdefs, compdef, compkeys, has_comp, first_period,
    last_period, compmodule, compname, compinstance, dim_keys, dim_values,
    set_first_last!

my_model = Model()

# Try running model with no components
@test length(compdefs(my_model)) == 0
@test length(my_model) == 0
@test_throws ErrorException run(my_model)

# Now add several components to the module
@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    """
    Test docstring.
    """
    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

@defcomp testcomp2 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

@defcomp testcomp3 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])
    cbox = Variable(index=[time, 5])    # anonymous dimension

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

# Testing that you cannot add a component without a time dimension
@test_throws ErrorException add_comp!(my_model, testcomp1)

# Start building up the model
set_dimension!(my_model, :time, 2015:5:2110)
add_comp!(my_model, testcomp1)

# Testing that you cannot add two components of the same name
@test_throws ErrorException add_comp!(my_model, testcomp1)

# Testing to catch adding component twice
@test_throws ErrorException add_comp!(my_model, testcomp1)

# Testing to catch if before or after does not exist
@test_throws ErrorException add_comp!(my_model, testcomp2, before=:testcomp3)

@test_throws ErrorException add_comp!(my_model, testcomp2, after=:testcomp3)

# Add more components to model
add_comp!(my_model, testcomp2)
add_comp!(my_model, testcomp3)

# Check addition of anonymous dimension
dvals = dim_values(my_model.md, Symbol(5))
dkeys = dim_keys(my_model.md, Symbol(5))
@test dvals == dkeys

comps = collect(compdefs(my_model))

# Test compdefs, compdef, compkeys, etc.
@test comps == collect(compdefs(my_model.md))
@test length(comps) == 3
@test compdef(my_model, :testcomp3).comp_id == comps[3].comp_id
@test_throws KeyError compdef(my_model, :testcomp4) #this component does not exist
@test [compkeys(my_model.md)...] == [:testcomp1, :testcomp2, :testcomp3]
@test has_comp(my_model.md, :testcomp1) == true
@test has_comp(my_model.md, :testcomp4) == false

@test compmodule(testcomp3) == Main.TestComponents
@test compname(testcomp3) == :testcomp3

@test length(my_model) == 3
add_comp!(my_model, testcomp3, :testcomp3_v2)
@test length(my_model) == 4

#------------------------------------------------------------------------------
#   Tests for component run periods when resetting the model's time dimension
#------------------------------------------------------------------------------

@defcomp testcomp1 begin
    var1 = Variable(index=[time])
    par1 = Parameter(index=[time])

    function run_timestep(p, v, d, t)
        v.var1[t] = p.par1[t]
    end
end

# 1. Test resetting the time dimension without explicit first/last values

comp_def = testcomp1
@test comp_def.first === nothing   # original component definition's first and last values are unset
@test comp_def.last === nothing

m = Model()
set_dimension!(m, :time, 2001:2005)
add_comp!(m, testcomp1, :C) # Don't set the first and last values here
comp_def = compdef(m.md, :C)      # Get the component definition in the model
@test comp_def.first === 2001
@test comp_def.last === 2005

update_param!(m, :C, :par1, zeros(5))
Mimi.build!(m)              # Build the model
ci = compinstance(m, :C)    # Get the component instance
@test ci.first == 2001 && ci.last == 2005      # no change

set_dimension!(m, :time, 2000:2020) # Reset the time dimension
comp_def = compdef(m.md, :C)        # Get the component definition in the model
@test comp_def.first === 2001 && comp_def.last === 2005 # no change

update_param!(m, :C, :par1, zeros(21))
Mimi.build!(m)               # Build the model
ci = compinstance(m, :C)    # Get the component instance
@test ci.first == 2001 && ci.last == 2005 # no change

set_first_last!(m, :C, first = 2000, last = 2020)
comp_def = compdef(m.md, :C)        # Get the component definition in the model
@test comp_def.first == 2000 && comp_def.last == 2020 # change!

# 2. Test resetting the time dimension with explicit first/last values

m = Model()
set_dimension!(m, :time, 2000:2100)
add_comp!(m, testcomp1, :C; first = 2010, last = 2090)

comp_def = compdef(m.md, :C)      # Get the component definition in the model
@test comp_def.first == 2010 && comp_def.last == 2090

set_dimension!(m, :time, 1950:2090)
update_param!(m, :C, :par1, zeros(141))
Mimi.build!(m)               # Build the model

ci = compinstance(m, :C) # Get the component instance
@test ci.first == 2010 && ci.last == 2090 # The component instance's first and last values are the same as in the comp def

set_dimension!(m, :time, 1940:2200) # Reset the time dimension
update_param!(m, :C, :par1, zeros(261)) # Have to reset the parameter to have the same width as the model time dimension

comp_def = compdef(m.md, :C)      # Get the component definition in the model
@test comp_def.first == 2010      # First and last values should still be the same
@test comp_def.last == 2090

Mimi.build!(m)               # Build the model
ci = compinstance(m, :C) # Get the component instance
@test ci.first == 2010      # The component instance's first and last values are the same as the comp def
@test ci.last == 2090

end #module


text:
export make_kovesi_c8

kovesi_c8_data = Tuple{UInt8, UInt8, UInt8}[
    (203, 125, 111),
    (203, 126, 111),
    (203, 128, 110),
    (204, 130, 110),
    (204, 133, 110),
    (204, 135, 110),
    (204, 137, 109),
    (204, 139, 109),
    (204, 142, 109),
    (204, 144, 109),
    (204, 146, 108),
    (204, 148, 108),
    (204, 151, 108),
    (204, 153, 107),
    (204, 155, 107),
    (204, 157, 107),
    (204, 160, 107),
    (204, 162, 106),
    (203, 164, 106),
    (203, 166, 106),
    (203, 168, 105),
    (203, 170, 105),
    (203, 173, 104),
    (203, 175, 104),
    (203, 177, 104),
    (202, 179, 103),
    (202, 181, 103),
    (202, 183, 103),
    (202, 186, 102),
    (201, 188, 102),
    (201, 190, 101),
    (201, 192, 101),
    (201, 194, 100),
    (200, 196, 100),
    (200, 198, 99),
    (200, 200, 99),
    (199, 203, 99),
    (199, 205, 98),
    (199, 207, 98),
    (198, 209, 97),
    (198, 211, 96),
    (197, 213, 96),
    (197, 215, 95),
    (196, 217, 95),
    (196, 220, 94),
    (195, 222, 94),
    (195, 224, 93),
    (194, 226, 92),
    (194, 228, 92),
    (193, 230, 91),
    (192, 232, 90),
    (192, 234, 90),
    (191, 236, 89),
    (191, 239, 88),
    (190, 241, 88),
    (189, 243, 87),
    (188, 245, 86),
    (188, 247, 85),
    (187, 249, 85),
    (186, 251, 84),
    (185, 253, 83),
    (184, 255, 83),
    (183, 255, 83),
    (182, 255, 83),
    (181, 255, 84),
    (180, 255, 85),
    (178, 255, 86),
    (177, 255, 88),
    (175, 254, 90),
    (174, 252, 92),
    (172, 251, 95),
    (171, 249, 97),
    (169, 248, 99),
    (167, 246, 101),
    (166, 244, 103),
    (164, 243, 105),
    (162, 241, 107),
    (161, 239, 108),
    (159, 238, 110),
    (157, 236, 112),
    (155, 234, 114),
    (154, 233, 116),
    (152, 231, 117),
    (150, 229, 119),
    (148, 228, 121),
    (146, 226, 122),
    (145, 225, 124),
    (143, 223, 126),
    (141, 221, 127),
    (139, 220, 129),
    (137, 218, 130),
    (135, 217, 132),
    (133, 215, 133),
    (131, 213, 135),
    (129, 212, 136),
    (127, 210, 138),
    (125, 209, 139),
    (123, 207, 141),
    (120, 205, 142),
    (118, 204, 143),
    (116, 202, 145),
    (114, 201, 146),
    (111, 199, 147),
    (109, 197, 149),
    (106, 196, 150),
    (104, 194, 151),
    (101, 193, 153),
    (99, 191, 154),
    (96, 190, 155),
    (93, 188, 157),
    (90, 187, 158),
    (87, 185, 159),
    (84, 184, 160),
    (81, 182, 162),
    (77, 180, 163),
    (74, 179, 164),
    (70, 177, 165),
    (66, 176, 167),
    (61, 174, 168),
    (57, 173, 169),
    (51, 171, 170),
    (46, 170, 172),
    (39, 168, 173),
    (31, 167, 174),
    (22, 165, 175),
    (8, 164, 176),
    (0, 162, 177),
    (0, 161, 177),
    (0, 159, 178),
    (0, 157, 178),
    (0, 155, 179),
    (0, 153, 179),
    (0, 152, 179),
    (0, 150, 179),
    (0, 148, 179),
    (1, 146, 179),
    (8, 144, 179),
    (14, 142, 179),
    (18, 140, 179),
    (21, 138, 179),
    (24, 137, 179),
    (27, 135, 179),
    (29, 133, 179),
    (31, 131, 179),
    (33, 129, 180),
    (35, 127, 180),
    (36, 125, 180),
    (37, 123, 180),
    (39, 121, 180),
    (40, 120, 180),
    (41, 118, 180),
    (42, 116, 180),
    (43, 114, 180),
    (43, 112, 180),
    (44, 110, 180),
    (45, 108, 180),
    (45, 106, 180),
    (46, 105, 180),
    (46, 103, 180),
    (47, 101, 180),
    (47, 99, 180),
    (47, 97, 180),
    (47, 95, 180),
    (48, 93, 180),
    (48, 91, 180),
    (48, 89, 179),
    (48, 87, 179),
    (48, 85, 179),
    (48, 83, 179),
    (48, 82, 179),
    (47, 80, 179),
    (47, 78, 179),
    (47, 76, 179),
    (47, 74, 179),
    (46, 72, 179),
    (46, 70, 179),
    (46, 68, 179),
    (45, 66, 179),
    (45, 63, 179),
    (44, 61, 179),
    (43, 59, 179),
    (43, 57, 179),
    (42, 55, 179),
    (41, 53, 179),
    (40, 51, 179),
    (40, 48, 179),
    (39, 46, 179),
    (38, 44, 178),
    (38, 42, 178),
    (38, 40, 178),
    (39, 38, 178),
    (42, 37, 177),
    (45, 36, 177),
    (50, 36, 176),
    (54, 37, 176),
    (59, 37, 175),
    (64, 38, 174),
    (69, 40, 173),
    (73, 41, 172),
    (77, 42, 171),
    (81, 43, 170),
    (85, 44, 169),
    (88, 46, 168),
    (92, 47, 167),
    (95, 48, 166),
    (98, 50, 165),
    (101, 51, 164),
    (104, 52, 163),
    (107, 54, 162),
    (110, 55, 161),
    (113, 56, 160),
    (115, 58, 159),
    (118, 59, 158),
    (120, 60, 157),
    (123, 62, 156),
    (125, 63, 155),
    (128, 65, 154),
    (130, 66, 153),
    (133, 68, 152),
    (135, 69, 151),
    (137, 70, 150),
    (139, 72, 149),
    (142, 73, 148),
    (144, 75, 147),
    (146, 76, 146),
    (148, 78, 145),
    (150, 79, 144),
    (152, 81, 143),
    (154, 82, 142),
    (156, 83, 141),
    (158, 85, 140),
    (160, 86, 138),
    (162, 88, 137),
    (164, 89, 136),
    (166, 91, 135),
    (168, 92, 134),
    (170, 94, 133),
    (172, 95, 132),
    (174, 97, 131),
    (175, 98, 130),
    (177, 100, 129),
    (179, 102, 127),
    (181, 103, 126),
    (183, 105, 125),
    (185, 106, 124),
    (186, 108, 123),
    (188, 109, 122),
    (190, 111, 121),
    (192, 112, 119),
    (193, 114, 118),
    (195, 115, 117),
    (197, 117, 116),
    (198, 118, 115),
    (200, 120, 114),
    (201, 121, 113),
    (202, 123, 112),
]

make_kovesi_c8() = matplotlib.colors.ListedColormap(
        [[(d ./ 255.0)..., 1] for d in kovesi_c8_data])


text:
using Documenter, AOPE
makedocs(sitename="Documentação")


text:
"""
    JuMP.termination_status(model::InfiniteModel)

Extend [`termination_status`](@ref JuMP.termination_status(::JuMP.Model)) to
return the `MOI.TerminationStatus` in accordance with the optimizer model.

**Example**
```julia-repl
julia> termination_status(model)
LOCALLY_SOLVED::TerminationStatusCode = 4
```
"""
function JuMP.termination_status(model::InfiniteModel)
    return JuMP.termination_status(optimizer_model(model))
end

"""
    JuMP.raw_status(model::InfiniteModel)

Extend [`raw_status`](@ref JuMP.raw_status(::JuMP.Model)) to return the status
reported by the solver in accordance with the optimizer model.

**Example**
```julia-repl
julia> raw_status(model) # Ipopt
"Solve_Succeeded"
```
"""
function JuMP.raw_status(model::InfiniteModel)
    return JuMP.raw_status(optimizer_model(model))
end

"""
     JuMP.primal_status(model::InfiniteModel; [result::Int = 1])

Extend [`primal_status`](@ref JuMP.primal_status(::JuMP.Model)) to return the
`MOI.PrimalStatus` reported in accordance with the optimizer model and the
result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> primal_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
```
"""
function  JuMP.primal_status(model::InfiniteModel; result::Int = 1)
    return JuMP.primal_status(optimizer_model(model); result = result)
end

"""
    JuMP.dual_status(model::InfiniteModel; [result::Int = 1])

Extend [`dual_status`](@ref JuMP.dual_status(::JuMP.Model)) to return the
`MOI.DualStatus` reported in accordance with the optimizer model and the
result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> dual_status(model)
FEASIBLE_POINT::ResultStatusCode = 1
```
"""
function JuMP.dual_status(model::InfiniteModel; result::Int = 1)
    return JuMP.dual_status(optimizer_model(model); result = result)
end

"""
    JuMP.solve_time(model::InfiniteModel)

Extend [`solve_time`](@ref JuMP.solve_time(::JuMP.Model)) to return the
time used by the solver to terminate reported in accordance with the optimizer
model. This will error if not supported by the solver.

**Example**
```julia-repl
julia> solve_time(model)
0.004999876022338867
```
"""
function JuMP.solve_time(model::InfiniteModel)
    return JuMP.solve_time(optimizer_model(model))
end

"""
    JuMP.has_values(model::InfiniteModel; [result::Int = 1])::Bool

Extend [`has_values`](@ref JuMP.has_values(::JuMP.Model)) to return a `Bool`
whether variable values are available in accordance with the optimizer model and
the result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> has_values(model)
true
```
"""
function JuMP.has_values(model::InfiniteModel; result::Int = 1)::Bool
    return JuMP.primal_status(model; result = result) != MOI.NO_SOLUTION
end

"""
    JuMP.has_duals(model::InfiniteModel; [result::Int = 1])::Bool

Extend [`has_duals`](@ref JuMP.has_duals(::JuMP.Model)) to return a `Bool`
whether constraint duals are available in accordance with the optimizer model and
the result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> has_duals(model)
true
```
"""
function JuMP.has_duals(model::InfiniteModel; result::Int = 1)::Bool
    return JuMP.dual_status(model; result = result) != MOI.NO_SOLUTION
end

"""
    JuMP.objective_bound(model::InfiniteModel)::Float64

Extend [`objective_bound`](@ref JuMP.objective_bound(::JuMP.Model)) to return
the objective bound in accordance with the optimizer model.

**Example**
```julia-repl
julia> objective_bound(model)
42.0
```
"""
function JuMP.objective_bound(model::InfiniteModel)::Float64
    return JuMP.objective_bound(optimizer_model(model))
end

"""
    JuMP.objective_value(model::InfiniteModel; [result::Int = 1])::Float64

Extend [`objective_value`](@ref JuMP.objective_value(::JuMP.Model)) to return
the objective value in accordance with the optimizer model and the
result index `result` of the most recent solution obtained.

**Example**
```julia-repl
julia> objective_value(model)
42.0
```
"""
function JuMP.objective_value(model::InfiniteModel; result::Int = 1)::Float64
    return JuMP.objective_value(optimizer_model(model); result = result)
end

"""
    JuMP.dual_objective_value(model::InfiniteModel; [result::Int = 1])::Float64

Extend [`dual_objective_value`](@ref JuMP.dual_objective_value(::JuMP.Model)) to
return the dual objective value in accordance with the optimizer model and the
result index `result` of the most recent solution obtained. Errors
if the solver does not support this.

**Example**
```julia-repl
julia> dual_objective_value(model)
42.0
```
"""
function JuMP.dual_objective_value(model::InfiniteModel; result::Int = 1)::Float64
    return JuMP.dual_objective_value(optimizer_model(model); result = result)
end

"""
    map_value(ref, key::Val{ext_key_name}, result::Int)

Map the value(s) of `ref` to its counterpart in the optimizer model type that is
distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default. Here `result` is the result index used in `value`.
"""
function map_value end

# Default method that depends on optimizer_model_variable --> making extensions easier
function map_value(vref::GeneralVariableRef, key, result::Int)
    opt_vref = optimizer_model_variable(vref, key)
    if opt_vref isa AbstractArray
        return JuMP.value.(opt_vref; result = result)
    else
        return JuMP.value(opt_vref; result = result)
    end
end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_value(cref::GeneralConstraintRef, key, result::Int)
    opt_cref = optimizer_model_constraint(cref)
    if opt_cref isa AbstractArray
        return JuMP.value.(opt_cref; result = result)
    else
        return JuMP.value(opt_cref; result = result)
    end
end

"""
    JuMP.value(vref::InfOptVariableRef; [result::Int = 1])

Extend [`JuMP.value`](@ref JuMP.value(::JuMP.VariableRef)) to return the value(s)
of `vref` in accordance with its reformulation variable(s) stored in the optimizer
model and the result index `result` of the most recent solution obtained. Use
[`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel)) to check
if a result exists before asking for values. For extensions, this only works if
[`optimizer_model_variable`](@ref) has been extended correctly and/or
[`map_value`](@ref) has been extended for variables.

**Example**
```julia-repl
julia> value(z)
42.0
```
"""
function JuMP.value(vref::InfOptVariableRef; result::Int = 1)
    return map_value(vref, Val(optimizer_model_key(JuMP.owner_model(vref))), result)
end

"""
    JuMP.value(cref::GeneralConstraintRef; [result::Int = 1])

Extend [`JuMP.value`](@ref JuMP.value(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the value(s) of `cref` in accordance with its reformulation constraint(s)
stored in the optimizer model and the result index `result` of the most recent
solution obtained. Use [`JuMP.has_values`](@ref JuMP.has_values(::InfiniteModel))
to check if a result exists before asking for values. For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_value`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> value(c1)
4-element Array{Float64,1}:
 -0.0
 20.9
 20.9
 20.9
```
"""
function JuMP.value(cref::GeneralConstraintRef; result::Int = 1)
    return map_value(cref, Val(optimizer_model_key(JuMP.owner_model(cref))), result)
end

"""
    map_optimizer_index(ref, key::Val{ext_key_name})

Map the `MathOptInterface` index(es) of `ref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default.
"""
function map_optimizer_index end

# Default method that depends on optimizer_model_variable --> making extensions easier
function map_optimizer_index(vref::GeneralVariableRef, key)
    opt_vref = optimizer_model_variable(vref, key)
    if opt_vref isa AbstractArray
        return JuMP.optimizer_index.(opt_vref)
    else
        return JuMP.optimizer_index(opt_vref)
    end
end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_optimizer_index(cref::GeneralConstraintRef, key)
    opt_cref = optimizer_model_constraint(cref)
    if opt_cref isa AbstractArray
        return JuMP.optimizer_index.(opt_cref)
    else
        return JuMP.optimizer_index(opt_cref)
    end
end

"""
    JuMP.optimizer_index(vref::InfOptVariableRef)

Extend [`JuMP.optimizer_index`](@ref JuMP.optimizer_index(::JuMP.VariableRef)) to
return the `MathOptInterface` index(es) of `vref` in accordance with its
reformulation variable(s) stored in the optimizer model. For extensions, this
only works if [`optimizer_model_variable`](@ref) has been extended correctly
and/or [`map_optimizer_index`](@ref) has been extended for variables.

**Example**
```julia-repl
julia> optimizer_index(x)
4-element Array{MathOptInterface.VariableIndex,1}:
 MathOptInterface.VariableIndex(2)
 MathOptInterface.VariableIndex(3)
 MathOptInterface.VariableIndex(4)
 MathOptInterface.VariableIndex(5)
```
"""
function JuMP.optimizer_index(vref::InfOptVariableRef)
    return map_optimizer_index(vref, Val(optimizer_model_key(JuMP.owner_model(vref))))
end

"""
    JuMP.optimizer_index(cref::GeneralConstraintRef)

Extend [`JuMP.optimizer_index`](@ref JuMP.optimizer_index(::JuMP.ConstraintRef{JuMP.Model}))
to return the `MathOptInterface` index(es) of `cref` in accordance with its
reformulation constraints(s) stored in the optimizer model. For extensions, this
only works if [`optimizer_model_constraint`](@ref) has been extended correctly
and/or [`map_optimizer_index`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> optimizer_index(c1)
4-element Array{MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}},1}:
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(1)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(2)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(3)
 MathOptInterface.ConstraintIndex{MathOptInterface.ScalarAffineFunction{Float64},MathOptInterface.GreaterThan{Float64}}(4)
```
"""
function JuMP.optimizer_index(cref::GeneralConstraintRef)
    return map_optimizer_index(cref, Val(optimizer_model_key(JuMP.owner_model(cref))))
end

"""
    map_dual(cref::GeneralConstraintRef, key::Val{ext_key_name}, result::Int)

Map the dual(s) of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default. Here `result` is the result index that is used in `dual`.
"""
function map_dual end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_dual(cref::GeneralConstraintRef, key, result::Int)
    opt_cref = optimizer_model_constraint(cref)
    if opt_cref isa AbstractArray
        return JuMP.dual.(opt_cref; result = result)
    else
        return JuMP.dual(opt_cref; result = result)
    end
end

"""
    JuMP.dual(cref::GeneralConstraintRef; [result::Int = 1])

Extend [`JuMP.dual`](@ref JuMP.dual(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the dual(s) of `cref` in accordance with its reformulation constraint(s)
stored in the optimizer model and the result index `result` of the most recent
solution obtained. Use [`JuMP.has_duals`](@ref JuMP.has_duals(::InfiniteModel))
to check if a result exists before asking for duals. For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_dual`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> dual(c1)
4-element Array{Float64,1}:
 -42.0
 -42.0
 32.3
 0.0
```
"""
function JuMP.dual(cref::GeneralConstraintRef; result::Int = 1)
    return map_dual(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                    result)
end

# Error redriect dor variable call
function JuMP.dual(vref::GeneralVariableRef; result::Int = 1)
    return JuMP.dual(JuMP.VariableRef(JuMP.Model(), MOI.VariableIndex(1)))
end

"""
    map_shadow_price(cref::GeneralConstraintRef, key::Val{ext_key_name})

Map the shadow price(s) of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `ref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable` and `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` and `optimizer_model_constraint` are used to make
these mappings by default.
"""
function map_shadow_price end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_shadow_price(cref::GeneralConstraintRef, key)
    opt_cref = optimizer_model_constraint(cref)
    if opt_cref isa AbstractArray
        return JuMP.shadow_price.(opt_cref)
    else
        return JuMP.shadow_price(opt_cref)
    end
end

"""
    JuMP.shadow_price(cref::GeneralConstraintRef)

Extend [`JuMP.shadow_price`](@ref JuMP.shadow_price(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the shadow price(s) of `cref` in accordance with its reformulation constraint(s)
stored in the optimizer model. Use [`JuMP.has_duals`](@ref JuMP.has_duals(::InfiniteModel))
to check if a result exists before asking for duals. For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_shadow_price`](@ref) has been extended for constraints.

**Example**
```julia-repl
julia> shadow_price(c1)
4-element Array{Float64,1}:
 42.0
 42.0
 -32.3
 -0.0
```
"""
function JuMP.shadow_price(cref::GeneralConstraintRef)
    return map_shadow_price(cref, Val(optimizer_model_key(JuMP.owner_model(cref))))
end

"""
    map_lp_rhs_perturbation_range(cref::GeneralConstraintRef,
                                  key::Val{ext_key_name}, toler::Float64)

Map the RHS perturbation range of `cref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `cref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_constraint`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_constraint` is used to make these mappings by default. Here
`toler` corresponds to the `feasibility_tolerance` used by `lp_rhs_perturbation_range`.
"""
function map_lp_rhs_perturbation_range end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_lp_rhs_perturbation_range(cref::GeneralConstraintRef, key,
                                       toler::Float64)
    opt_cref = optimizer_model_constraint(cref)
    if opt_cref isa AbstractArray
        return JuMP.lp_rhs_perturbation_range.(opt_cref;
                                               feasibility_tolerance = toler)
    else
        return JuMP.lp_rhs_perturbation_range(opt_cref;
                                              feasibility_tolerance = toler)
    end
end

"""
    JuMP.lp_rhs_perturbation_range(cref::GeneralConstraintRef;
                                   [feasibility_tolerance::Float64 = 1e-8])

Extend [`JuMP.lp_rhs_perturbation_range`](@ref JuMP.lp_rhs_perturbation_range(::JuMP.ConstraintRef{JuMP.Model, <:JuMP._MOICON}))
to return the range(s) of the RHS of `cref` for which the shadow price(s) are valid
in accordance with its reformulation constraint(s)
stored in the optimizer model. For extensions, this only
works if [`optimizer_model_constraint`](@ref) has been extended correctly and/or
[`map_lp_rhs_perturbation_range`](@ref) has been implemented.

**Example**
```julia-repl
julia> lp_rhs_perturbation_range(c1)
4-element Array{Tuple{Float64,Float64},1}:
 (-42.0, Inf)
 (-Inf, 42.0)
 (-Inf, 42.0)
 (-Inf, 42.0)
```
"""
function JuMP.lp_rhs_perturbation_range(cref::GeneralConstraintRef;
                                        feasibility_tolerance::Float64 = 1e-8)
    return map_lp_rhs_perturbation_range(cref, Val(optimizer_model_key(JuMP.owner_model(cref))),
                                         feasibility_tolerance)
end

"""
    map_lp_objective_perturbation_range(vref::InfOptVariableRef,
                                        key::Val{ext_key_name}, toler::Float64)

Map the reduced cost range(s) of `vref` to its counterpart in the optimizer
model type that is distininguished by its extension key `key` as type `Val{ext_key_name}`.
Here `vref` need refer to methods for both variable references and constraint
references. This only needs to be defined for reformulation extensions that cannot
readily extend `optimizer_model_variable`.
Such as is the case with reformuations that do not have a direct mapping between
variables and/or constraints in the original infinite form. Otherwise,
`optimizer_model_variable` is used to make these mappings by default. Here
`toler` corresponds to the `optimality_tolerance` used by
`lp_objective_perturbation_range`.
"""
function map_lp_objective_perturbation_range end

# Default method that depends on optimizer_model_constraint --> making extensions easier
function map_lp_objective_perturbation_range(vref::GeneralVariableRef, key,
                                             toler::Float64)
    opt_vref = optimizer_model_variable(vref)
    if opt_vref isa AbstractArray
        return JuMP.lp_objective_perturbation_range.(opt_vref;
                                                     optimality_tolerance = toler)
    else
        return JuMP.lp_objective_perturbation_range(opt_vref;
                                                    optimality_tolerance = toler)
    end
end

"""
    JuMP.lp_objective_perturbation_range(vref::InfOptVariableRef;
                                         [optimality_tolerance::Float64 = 1e-8])

Extend [`JuMP.lp_objective_perturbation_range`](@ref JuMP.lp_objective_perturbation_range(::JuMP.VariableRef))
to return the range(s) that the reduced cost(s) of `vref` remain valid
in accordance with its reformulation variables(s)
stored in the optimizer model. For extensions, this only
works if [`optimizer_model_variable`](@ref) has been extended correctly and/or
if [`map_lp_objective_perturbation_range`](@ref) has been implemented.

**Example**
```julia-repl
julia> lp_objective_perturbation_range(z)
(-2.0, Inf)
```
"""
function JuMP.lp_objective_perturbation_range(vref::GeneralVariableRef;
                                              optimality_tolerance::Float64 = 1e-8)
    return map_lp_objective_perturbation_range(vref, Val(optimizer_model_key(JuMP.owner_model(vref))),
                                               optimality_tolerance)
end


text:
using SimpleSolvers
using Test
using GeometricIntegrators
using GeometricIntegrators.Integrators.VPRK
using GeometricIntegrators.Utils
using GeometricProblems.LotkaVolterra2dSymmetric
using GeometricProblems.LotkaVolterra2dSymmetric: Δt, nt, reference_solution

SimpleSolvers.set_config(:nls_atol, 8eps())
SimpleSolvers.set_config(:nls_rtol, 2eps())

@testset "$(rpad("Lotka-Volterra 2D with symmetric Lagrangian",80))" begin
    ode  = lotka_volterra_2d_ode()
    iode = lotka_volterra_2d_iode()
    idae = lotka_volterra_2d_idae()

    int = Integrator(ode, TableauGauss(2), Δt)
    sol = integrate(ode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 4E-4

    int = IntegratorVPRKpMidpoint(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    # println(relative_maximum_error(sol.q, reference_solution))
    @test relative_maximum_error(sol.q, reference_solution) < 2E-3

    int = IntegratorVPRKpSymmetric(iode, TableauVPGLRK(2), Δt)
    sol = integrate(iode, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 5E-4

    int = Integrator(idae, TableauVSPARKGLRKpMidpoint(2), Δt)
    sol = integrate(idae, int, nt)
    # println(relative_maximum_error(sol.q, reference_solution))
    @test relative_maximum_error(sol.q, reference_solution) < 2E-3

    int = Integrator(idae, TableauVSPARKGLRKpSymmetric(2), Δt)
    sol = integrate(idae, int, nt)
    @test relative_maximum_error(sol.q, reference_solution) < 5E-4
end


text:
# Sparsity Patterns for the Hessian Matrix of rod
function sparse_append!(I, J, rows::AbstractVector{T}, cols::AbstractVector{T}) where T
    oner = ones(T, size(rows))
    onec = ones(T, size(cols))
    append!(I, kron(onec, rows))
    append!(J, kron(cols, oner))
end

function hessian_sparsity(xind, yind, zind, θind) where Ti

    T = eltype(xind)
    nv = length(xind)
    ne = nv-1

    vind = hcat(xind, yind, zind)
    tind = reshape(θind, length(θind), 1)
    I = Vector{T}()
    J = Vector{T}()

    # x - x interactions
    for i = 1:nv, j = max(i-2, 1):min(i+2, nv)
        sparse_append!(I, J, vind[i,:], vind[j,:])
    end

    # x - θ interactions
    for i = 1:nv, j = max(i-2, 1):min(i+1, ne)
        sparse_append!(I, J, vind[i,:], tind[j,:])
    end

    # θ - x interactions
    for i = 1:ne, j = max(i-1, 1):min(i+2, nv)
        sparse_append!(I, J, tind[i,:], vind[j,:])
    end

    # θ - θ interactions
    for i = 1:ne, j = max(i-1, 1):min(i+1, ne)
        sparse_append!(I, J, tind[i,:], tind[j,:])
    end

    V = ones(T, length(I))
    sparse(I, J ,V)
end

function banded_hessian_sparsity(ns)
    xind = 1 .+ 4*(0:ns+1)
    yind = xind .+ 1
    zind = xind .+ 2
    θind = 4 .+ 4*(0:ns)
    hessian_sparsity(xind, yind, zind, θind)
end

function block_hessian_sparsity(ns)
    xind = 1 .+ 3*(0:ns+1)
    yind = xind .+ 1
    zind = xind .+ 2
    θind = 3*(ns+2) .+ (1:ns+1)
    hessian_sparsity(xind, yind, zind, θind)
end

function coordinate_hessian_sparsity(ns)
    nv = ns + 2
    xind = 1:nv
    yind = xind .+ nv
    zind = yind .+ nv
    θind = 3*nv .+ (1:ns+1)
    hessian_sparsity(xind, yind, zind, θind)
end

# using Plots
# spy(banded_hessian_sparsity(ns))
# spy(block_hessian_sparsity(ns))
# spy(coordinate_hessian_sparsity(ns))

# ------------------------------------------------------------------------------
# Convert rod_delta to and from a vector representation
# ------------------------------------------------------------------------------
function banded_storage!(q, dr::rod_delta)
    N = length(dr.Δθ) - 1
    q[1:4:4*N+5] .= @view dr.Δx[:,1]
    q[2:4:4*N+6] .= @view dr.Δx[:,2]
    q[3:4:4*N+7] .= @view dr.Δx[:,3]
    q[4:4:4*N+4] .= @view dr.Δθ[:]
end

function banded_storage!(dr::rod_delta, q)
    N = length(dr.Δθ) - 1
    dr.Δx[:,1] .= @view q[1:4:4*N+5]
    dr.Δx[:,2] .= @view q[2:4:4*N+6]
    dr.Δx[:,3] .= @view q[3:4:4*N+7]
    dr.Δθ[:] .= @view q[4:4:4*N+4]
end

function block_storage!(q, dr::rod_delta)
    N = length(dr.Δθ) - 1
    q[1:3:3*N+4] .= @view dr.Δx[:,1]
    q[2:3:3*N+5] .= @view dr.Δx[:,2]
    q[3:3:3*N+6] .= @view dr.Δx[:,3]
    q[3*N+7:4*N+7] .= @view dr.Δθ[:]
end

function block_storage!(dr::rod_delta, q)
    N = length(dr.Δθ) - 1
    dr.Δx[:,1] .= @view q[1:3:3*N+4]
    dr.Δx[:,2] .= @view q[2:3:3*N+5]
    dr.Δx[:,3] .= @view q[3:3:3*N+6]
    dr.Δθ[:] .= @view q[3*N+7:4*N+7]
end

function coordinate_storage!(q, dr::rod_delta)
    N = length(dr.Δθ) - 1
    q[1:3*(N+2)] .= @view dr.Δx[:]
    q[3*N+7:4*N+7] .= view(dr.Δθ, :, 1)
end

function coordinate_storage!(dr::rod_delta, q)
    N = length(dr.Δθ) - 1
    dr.Δx[:,1] .= @view q[1:(N+2)]
    dr.Δx[:,2] .= @view q[N+3:2*(N+2)]
    dr.Δx[:,3] .= @view q[2*N+5:3*(N+2)]
    dr.Δθ .= @view q[3*N+7:4*N+7]
end

# ------------------------------------------------------------------------------
# Allocating versions of the above
# ------------------------------------------------------------------------------
for _storage = (:banded_storage, :block_storage, :coordinate_storage)
    _storage! = Symbol(_storage, "!")
    @eval begin
        function $(_storage)(dr::rod_delta{T}) where T
            N = length(dr.Δθ) - 1
            q = Vector{T}(undef, 4*N+7)
            $(_storage!)(q, dr)
            return q
        end

        function $(_storage)(q::Vector{T}) where T
            N = div(length(q) - 7, 4)
            dr = allocate_delta(T, N)
            $(_storage!)(dr, q)
            return dr
        end
    end
end

# Check on storage orders
# using BenchmarkTools
# ns = 500
# q = Float64.(1:4*ns+7)
# dr = DER.allocate_delta(Float64, ns)
# @btime DER.banded_storage!($dr, $q)
# @btime DER.banded_storage!($q, $dr)
# @btime DER.block_storage!($dr, $q)
# @btime DER.block_storage!($q, $dr)
# @btime DER.coordinate_storage!($dr, $q)
# @btime DER.coordinate_storage!($q, $dr)


text:
#Tests for chunking
cfn1 = tempname2()
nccreate(cfn1, "v1", "d1", 10, "d2", 20, chunksize = (5, 2))
ncwrite(rand(10, 20), cfn1, "v1")
v = NetCDF.open(cfn1, "v1")
@test v.chunksize == (2, 5)
@test NetCDF.eachchunk(v) == NetCDF.GridChunks((10,20),(5,2))


text:
#                        .-'''-.                               _..._
#                       '   _    \         _______          .-'_..._''.
#  /|                 /   /` '.   \        \  ___ `'.     .' .'      '.\
#  ||                .   |     \  '         ' |--.\  \   / .'
#  ||        .-,.--. |   '      |  '        | |    \  ' . '                                     .|
#  ||  __    |  .-. |\    \     / /  __     | |     |  '| |                 __                .' |_
#  ||/'__ '. | |  | | `.   ` ..' /.:--.'.   | |     |  || |              .:--.'.         _  .'     |
#  |:/`  '. '| |  | |    '-...-'`/ |   \ |  | |     ' .'. '             / |   \ |      .' |'--.  .-'
#  ||     | || |  '-             `" __ | |  | |___.' /'  \ '.          .`" __ | |     .   | / |  |
#  ||\    / '| |                  .'.''| | /_______.'/    '. `._____.-'/ .'.''| |   .'.'| |// |  |
#  |/\'..' / | |                 / /   | |_\_______|/       `-.______ / / /   | |_.'.'.-'  /  |  '.'
#  '  `'-'`  |_|                 \ \._,\ '/                          `  \ \._,\ '/.'   \_.'   |   /
#                                 `--'  `"                               `--'  `"             `'-'

using Base.Broadcast
using Base.Broadcast: Broadcasted, DefaultArrayStyle, instantiate

# Structural utilities

using Base: tail

tcat(x) = x
tcat(x, y, z...) = tcat((x..., y...), z...)

broadcast_args(x) = (x,)
broadcast_args(bc::Broadcasted) = tcat(map(broadcast_args, bc.args)...)

accum_sum(xs, dims = :) = reduce(accum, xs, dims = dims)

trim(x, Δ) = reshape(Δ, ntuple(i -> size(Δ, i), Val(ndims(x))))

unbroadcast(x::AbstractArray, Δ) =
  size(x) == size(Δ) ? Δ :
  length(x) == length(Δ) ? trim(x, Δ) :
    trim(x, accum_sum(Δ, dims = ntuple(i -> size(x, i) == 1 ? i : ndims(Δ)+1, Val(ndims(Δ)))))

unbroadcast(x::Union{Number,Ref}, Δ) = accum_sum(Δ)

# Trivial Special Cases
# TODO fix this up and use it

Jtrivial(f, a...) = nothing
Jtrivial(::typeof(+), a...) = a
Jtrivial(::typeof(-), a, b) = (a..., .-b...)

trivia(_) = (1,)
function trivia(bc::Broadcasted)
  t = map(trivia, bc.args)
  any(t -> t === nothing, t) && return
  Jtrivial(bc.f, t...)
end

Joutput(f, a...) = nothing
Joutput(::typeof(exp), x) = map(t -> y -> y*t, x)

function Jbroadcast(bc::Broadcasted)
  t = map(trivia, bc.args)
  any(t -> t === nothing, t) && return
  Joutput(bc.f, t...)
end

@inline function unbroadcast_t(x, y, ȳ, j::J) where J
  trim(x, j.(y).*ȳ)
end

@inline function unbroadcast_t(x::Number, y, ȳ, j::J) where J
  x̄ = zero(float(x))
  @simd for I in eachindex(y)
    @inbounds x̄ += j(y[I])*ȳ[I]
  end
  return x̄
end

function ∇broadcast_t(bc::Broadcasted, J)
  y = copy(instantiate(bc))
  back(ȳ) = map(unbroadcast_t, broadcast_args(bc), map(_ -> y, J), map(_ -> ȳ, J), J)
  return y, back
end

# Reverse Mode

# TODO: forward context appropriately
# multi-output map for better performance
function ∇broadcast_r(bc::Broadcasted)
  bc′, unflatten = _forward(Broadcast.flatten, bc)
  len = Val(length(bc′.args)+1)
  y∂b = broadcast(_forward, bc′.f, bc′.args...)
  y = map(x -> x[1], y∂b)
  ∂b = map(x -> x[2], y∂b)
  y, function (ȳ)
    dxs_zip = map((∂b, ȳ) -> ∂b(ȳ), ∂b, ȳ)
    dxs = ntuple(i -> map(x -> x[i], dxs_zip), len)
    (f = accum_sum(dxs[1]),
     args = map(unbroadcast, bc′.args, Base.tail(dxs)),
     axes = nothing) |> unflatten |> Base.tail
  end
end

function ∇broadcast_r(bc::Broadcasted{<:DefaultArrayStyle{0}})
  bc′, unflatten = _forward(Broadcast.flatten, bc)
  len = Val(length(bc′.args)+1)
  y, ∂b = broadcast(_forward, bc′.f, bc′.args...)
  y, function (ȳ)
    dxs = ∂b(ȳ)
    (f = dxs[1],
     args = Base.tail(dxs),
     axes = nothing) |> unflatten |> Base.tail
  end
end

∇broadcast(bc::Broadcasted, ::Nothing) = ∇broadcast_r(bc)
∇broadcast(bc::Broadcasted, J) = ∇broadcast_t(bc, J)
∇broadcast(bc::Broadcasted) = ∇broadcast(bc, Jbroadcast(bc))

@adjoint Broadcast.materialize(bc::Broadcasted{<:DefaultArrayStyle}) = ∇broadcast_r(bc)


text:
#!/usr/bin/env julia
include("soscan.jl")
selc7017 = sosread("data/Statistikdatabasen_2019-02-21 18_39_26.csv")
caplot(selc7017, :Kvinnor, 45)
caplot(selc7017, :Män, 45)
caplot(selc7017, :Kvinnor, 65)
caplot(selc7017, :Män, 65)
grid(1)
legend(loc = 3)
xlabel("Kohort")
ylabel("log(incidens)")
title("Incidens lungcancer Sverige")
savefig("../../images/SvLc7017CohAge.svg")

max_stars_repo_name: klpn/static-dust

max_issues_repo_name: klpn/static-dust

max_forks_repo_name: klpn/static-dust


text:
using SeisPlot, SeisAcoustic

# velocity and density model
nz = 101; nx = 301
vel = 3000 * ones(nz, nx);  # m/s
vel[51:end,:] .= 3500;  # m/s
rho = 2000 * ones(nz, nx);  # kg/m^3
rho[25:end,:] .= 2300;  # m/s

# top boundary condition
free_surface = true   #(pml or free_surface)

# grid size
h = 10

# time step size and maximum modelling length
dt = 0.001; tmax = 2.0;  # use second as unit

# organize these parameters into a structure
fparams = FdParams(rho, vel, free_surface, h; data_format=Float64);
tparams = TdParams(rho, vel, free_surface, h, h, dt, tmax; data_format=Float64);

# recordings
# initialize recordings
irx = collect(1 : 2 : fparams.nx);
irz = 2 * ones(length(irx));
rec_f = Recordings(irz, irx, dt, tmax, fparams);
rec_t = Recordings(irz, irx, tparams);

# # single source
# src = Source(2, 150, dt, fparams; amp=100000);
#
# # generate recordings
# get_recordings!(rec_f, src, fparams, print_flag=true, consistent_flag=true);
# multi_step_forward!(rec_t, src, tparams);

# multiple sources
isx = collect(10:50:nx-10)
isz = 3 * ones(Int64, length(isx))
srcs = get_multi_sources(isz, isx, dt, fparams; amp=100000);
get_recordings!(rec_f, srcs, fparams, print_flag=true, consistent_flag=true);
multi_step_forward!(rec_t, srcs, tparams);

SeisPlotTX(rec_f.p);
SeisPlotTX(rec_t.p);
figure(); plot(rec_f.p[:,40]); plot(rec_t.p[:,40]);

#===============================================================================
#            simultaneouse source simulation in frequency domain
#===============================================================================

# forward modeling in frequency domain
pre_f = get_wavefield_FDFD(src, fparams; tmax=tmax, print_flag=true);

# compare against the time-domain modelling result
src = Source(2, 150, fparams; dt=0.001, amp=100000);
path  = joinpath(homedir(), "Desktop/boundary.rsb");
multi_step_forward!(path, src, tparams; save_flag="pressure");
(hdr, pre_t) = read_RSdata(path);

#plotting
it = 400
figure(figsize=(9,3)); imshow(pre_f[:,:,it*2], cmap="seismic");
figure(figsize=(9,3)); imshow(pre_t[:,:,it], cmap="seismic");

iz = 51; ix = 151;
figure(); plot(pre_f[iz,ix,:]); plot(pre_t[iz,ix,:])

# get recordigns in frequency domain
irx = collect(1:2:fparams.nx);
irz = 2 * ones(length(irx));
rec = Recordings(irz, irx, fparams);

# get Helmholtz operator
H = get_helmholtz_LU(fparams, 20.0);
b = zeros(Complex{Float64}, fparams.Nz * fparams.Nx);
idx = isx + fparams


text:
include(joinpath(pwd(), "src/models/biped.jl"))

include(joinpath(pwd(), "src/objectives/velocity.jl"))
include(joinpath(pwd(), "src/objectives/acceleration.jl"))

include(joinpath(pwd(), "src/constraints/contact.jl"))
include(joinpath(pwd(), "src/constraints/pinned_foot.jl"))

# Visualize
include(joinpath(pwd(), "src/models/visualize.jl"))
vis = Visualizer()
open(vis)

θ = pi / 10.0
# q1 = initial_configuration(model, θ)
q1, qT = loop_configurations(model, θ)
visualize!(vis, model, [q1])

# Horizon
T = 11

# Time step
tf = 1.0
h = tf / (T-1)

# Bounds

# ul <= u <= uu
_uu = Inf * ones(model.m)
_uu[model.idx_u] .= 10.0
_ul = zeros(model.m)
_ul[model.idx_u] .= -10.0
ul, uu = control_bounds(model, T, _ul, _uu)

xl, xu = state_bounds(model, T, x1 = [q1; q1], xT = [Inf * ones(model.nq); qT])

# Objective
q_ref = linear_interp(q1, qT, T)
X0 = configuration_to_state(q_ref)

obj_penalty = PenaltyObjective(1.0e5, model.m)

Qq = Diagonal(ones(model.nq))
Q = cat(0.5 * Qq, 0.5 * Qq, dims = (1, 2))
QT = cat(0.5 * Qq, 10.0 * Diagonal(ones(model.nq)), dims = (1, 2))
R = Diagonal([1.0e-1 * ones(model.nu)..., zeros(model.m - model.nu)...])

obj_tracking = quadratic_tracking_objective(
    [t < T ? Q : QT for t = 1:T],
    [R for t = 1:T-1],
    [X0[t] for t = 1:T],
    [zeros(model.m) for t = 1:T]
    )
obj_velocity = velocity_objective(
    [Diagonal(1.0 * ones(model.nq)) for t = 1:T-1],
    model.nq)
obj_acceleration = acceleration_objective(
    [Diagonal(1.0e-1 * ones(model.nq)) for t = 1:T-1],
    model.nq)
obj = MultiObjective([obj_tracking, obj_penalty, obj_velocity, obj_acceleration])

# Constraints
con_contact = contact_constraints(model, T)
con_pinned = pinned_foot_constraint(model, q1, T)
con = multiple_constraints([con_contact, con_pinned])

# Problem
prob = trajectory_optimization_problem(model,
               obj,
               T,
               h = h,
               xl = xl,
               xu = xu,
               ul = ul,
               uu = uu,
               con = con
               )

# trajectory initialization
q_ref = linear_interp(q1, qT, T)
X0 = configuration_to_state(q_ref) # linear interpolation on state
U0 = [0.001 * rand(model.m) for t = 1:T-1] # random controls

# Pack trajectories into vector
Z0 = pack(X0, U0, prob)

#NOTE: may need to run examples multiple times to get good trajectories
# Solve nominal problem

@time Z̄ = solve(prob, copy(Z0), tol = 1.0e-3, c_tol = 1.0e-3, mapl = 5)

check_slack(Z̄, prob)

X̄, Ū = unpack(Z̄, prob)

# Visualize

visualize!(vis, model, state_to_configuration(X̄), Δt = h)

plot(hcat(Ū...)[1:model.nu,:]', linetype = :steppost)

max_stars_repo_name: zacmanchester/DirectMotionPlanning

max_issues_repo_name: zacmanchester/DirectMotionPlanning

max_forks_repo_name: zacmanchester/DirectMotionPlanning


text:
using Printf

function update_info(info::vis.InfoBox, agent::Agent, world::World, t::Float64,
                     cv::Array{Float64, 1}=[NaN, NaN])
  (x, y, sx, sy, dx, u) = diagnostic(agent, world, t)

  dp = dx[3]
  th = atan(sy, sx)
  th_car = atan(dp, agent.x[2]) + atan(sy, sx)

  vis.update_vector_thr!(info.ds_vec, info.scaling * x,
                         info.scaling * y, th,
                         0.2 * agent.x[2] / 50.0)
  vis.update_vector_thr!(info.u1_vec, info.scaling * x,
                         info.scaling * y, th,
                         0.2 * u[1] / 10.0)
  vis.update_vector_thr!(info.u2_vec, info.scaling * x,
                         info.scaling * y, th + pi / 2,
                         0.2 * u[2] / 0.1)
  vis.update_vector_xy!(info.c_vec, info.scaling * x,
                        info.scaling * y, info.scaling * cv[1],
                        info.scaling * cv[2])

  offset = 0.0
  vis.update_text!(info.id_text, @sprintf("Agent %d", agent.id), info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.s_text, @sprintf("s= %3.1e", agent.x[1]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.ds_text, @sprintf("ds= %3.1e", agent.x[2]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.p_text, @sprintf("p= %3.1e", agent.x[3]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.u1_text, @sprintf("u1= %+3.1e", u[1]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.u2_text, @sprintf("u2= %+3.1e", u[2]),
                   info.x, info.y - offset, 1.0)
end

function update_info_from_tape(info::vis.InfoBox,
                               agent_x::AbstractArray{Float64, 1},
                               agent_id::Int, diag_data::Tuple,
                               cv::Array{Float64, 1}=[NaN, NaN])
  (x, y, sx, sy, dx, u) = diag_data

  dp = dx[3]
  th = atan(sy, sx)
  th_car = atan(dp, agent_x[2]) + atan(sy, sx)

  vis.update_vector_thr!(info.ds_vec, info.scaling * x,
                         info.scaling * y, th,
                         0.2 * agent_x[2] / 50.0)
  vis.update_vector_thr!(info.u1_vec, info.scaling * x,
                         info.scaling * y, th,
                         0.2 * u[1] / 10.0)
  vis.update_vector_thr!(info.u2_vec, info.scaling * x,
                         info.scaling * y, th + pi / 2,
                         0.2 * u[2] / 0.1)
  vis.update_vector_xy!(info.c_vec, info.scaling * x,
                        info.scaling * y, info.scaling * cv[1],
                        info.scaling * cv[2])

  offset = 0.0
  vis.update_text!(info.id_text, @sprintf("Agent %d", agent_id),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.s_text, @sprintf("s= %3.1e", agent_x[1]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.ds_text, @sprintf("ds= %3.1e", agent_x[2]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.p_text, @sprintf("p= %3.1e", agent_x[3]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.u1_text, @sprintf("u1= %+3.1e", u[1]),
                   info.x, info.y - offset, 1.0)
  (xmin, ymin, xmax, ymax) = vis.bounding_box(info.id_text)
  offset += ymax - ymin
  vis.update_text!(info.u2_text, @sprintf("u2= %+3.1e", u[2]),
                   info.x, info.y - offset, 1.0)
end


text:
using Documenter

#    format = Documenter.HTML(prettyurls = get(ENV, "CI", nothing) == "true"),

makedocs(
    sitename = "CatmullRom.jl",
    authors = "Jeffrey Sarnoff",
    pages = Any[
        "Overview" => "index.md",
        "Perspecive" => "perspective.md",
        "Three Functions" => "threefunctions.md",
        "Examples" => "Examples.md",
        "Points along a path" => "pointsalongapath.md",
        "A few hints" => "hints.md",
        "References" => "references.md"
    ]
)

deploydocs(
    repo = "github.com/JeffreySarnoff/CatmullRom.jl.git",
    target = "build"
)


text:
map_unit_to_segment(xi::Tensors.Vec{1,T}, x::Vector{Tensors.Vec{2,T}}) where {T} = x[1]*(1-xi[]) + x[2]*xi[]

δ(i, j) = i == j

function polygon_area(verts)
    N = size(verts,1)
    return 0.5*abs(sum(verts[j][1]*verts[mod1(j+1,N)][2]-verts[mod1(j+1,N)][1]*verts[j][2] for j ∈ 1:N))
end

""" Compute volume of a simplex spanned by vertices `verts` """
function simplex_area(verts::Vector{Tensors.Vec{N, T}}) where {N,T}
    # Volume of reference simplex element is 1/n!
    ref_verts = reference_coordinates(Simplex{N})
    A, b = get_affine_map(ref_verts, verts)
    F = svd(A)
    return *(F.S...)/factorial(N)
end

""" Map points points from reference simplex to simplex K """
function mapPointsFromReference(::Type{Simplex{N}}, K::Vector{Tensors.Vec{N, T}}, points::Vector{Tensors.Vec{N, T}}) where {N,T}
    ref_verts = reference_coordinates(Simplex{N})
    A, b = get_affine_map(ref_verts, K)
    return [Tensors.Vec{N}(A*x+Tensors.Vec{N}(b)) for x in points]
end

"""
get_affine_map(x, ξ)
Get (A,b) such that ξ = A * x + b is the affine
mapping from the simplex with vertices x to the simplex of vertices ξ.
"""
# function get_affine_map(x, ξ)
#     @assert length(x[1]) == length(ξ[1]) "dimension mismatch"
#     n = length(x[1])
#     T = eltype(x[1])
#
#     B = zeros(T, n*(n+1), n*(n+1))
#     L = zeros(T, n*(n+1))
#
#     for i in 1:length(x)
#         for j in 1:n
#             row = (i-1) * n + j
#             B[row, n*(j-1)+1:n*j] = x[i]
#             L[row] = ξ[i][j]
#             B[row, n * n + j] = one(T)
#         end
#     end
#     X = B\L
#     return reshape(X[1:n*n], (n, n)), X[n*n+1:end]
# end
function get_affine_map(x, ξ)
    @assert length(x[1]) == length(ξ[1]) "dimension mismatch"
    n = length(x[1])
    T = eltype(x[1])
    A = zeros(T,n,n)
    b = zeros(T,n)
    K = zeros(n+1,n+1)
    c = zeros(n+1)
    for i in 1:n
        K[:,n+1] .= 1
        for j in 1:size(x,1)
            K[j,1:n] .= x[j]
            c[j] = ξ[j][i]
        end
        a = K\c
        A[i,:] .= a[1:n]
        b[i] = a[n+1]
    end
    return A,b
end

""" Tesselate convex 2D shape with `vertices` using Delaunay """
function get_2Dtesselation(vertices)
    #Map points to min_coord <= x <= max_coord
    maxc = 0.0; minc = 0.0
    for x in vertices
        lmin,lmax = minmax(x[1],x[2])
        minc = min(lmin,minc); maxc=max(maxc,lmax)
    end
    width = VoronoiDelaunay.max_coord - VoronoiDelaunay.min_coord
    normalize(x) = (x-minc)/maxc*width+VoronoiDelaunay.min_coord

    tess = VoronoiDelaunay.DelaunayTessellation(2)
    push!(tess, VoronoiDelaunay.Point2D[VoronoiDelaunay.Point(normalize(x[1]),normalize(x[2])) for x in vertices])
    tesselation = Vector{Vector{Tensors.Vec{2,Float64}}}()
    for trig in tess
        mtrig = Vector{Tensors.Vec{2,Float64}}()
        for edge in (trig._a, trig._b, trig._c)
            x = (edge._x - VoronoiDelaunay.min_coord)*maxc*width + minc
            y = (edge._y - VoronoiDelaunay.min_coord)*maxc*width + minc
            push!(mtrig, Tensors.Vec{2}((x,y)))
        end
        push!(tesselation, mtrig)
    end
    return tesselation
end

max_stars_repo_name: Paulms/jFEMToools

max_issues_repo_name: Paulms/jFEMToools

max_forks_repo_name: Paulms/jFEMToools