text:
#=
ARPACK wrapper for use in Pseudospectra.jl

This file is part of Pseudospectra.jl.
Portions derived from the Base library of Julia (later moved to Arpack.jl).

Copyright (c) 2017 Ralph A. Smith
Copyright (c) 2009-2016: Jeff Bezanson, Stefan Karpinski, Viral B. Shah,
 and other contributors

SPDX-License-Identifier: MIT-Expat
License-Filename: LICENSES/MIT-Expat
=#

import LinearAlgebra: BlasInt, checksquare
import Arpack: ARPACKException, naupd, saupd, eupd_wrapper

"""
    xeigs(A, B, channel=nothing; nev=6, ncv=max(20,2*nev+1), which=:LM,
          tol=0.0, maxiter=300, sigma=nothing, ritzvec=true, v0=zeros((0,)),
          wantH=true, options=nothing)

Compute (generalized) eigenpairs of `A x=λ B x` and projections.

Modified version of `eigs()` (from the Arpack package, q.v.) which optionally
1. provides the projection matrices,
2. provides intermediate states (typically for plotting).

For intermediates, caller must provide a `Channel` argument;
in this case `xeigs` is implemented as a producer which fills `channel`.
When finished, it `put`s `(:finale, d,[v,],nconv,niter,nmult,resid[,H,V])` to
`channel`.

The type of `A` must provide methods for `mul!`, `issymmetric`,
`eltype` and `size`. Some variants are not
implemented for the case where `A` is not a factorable matrix.
"""
function xeigs(A, B, chnl = nothing;
               nev::Integer=6, ncv::Integer=max(20,2*nev+1), which=:LM,
               tol=0.0, maxiter::Integer=300, sigma=nothing,
               v0::Vector=zeros((0,)), # zeros(eltype(A),(0,)),
               ritzvec::Bool=true, wantH::Bool=true,
               options = Dict{Symbol,Any}())

    n = checksquare(A)

    T = eltype(A)
    iscmplx = T <: Complex
    sym = issymmetric(A) && issymmetric(B) && !iscmplx

    isgeneral = B !== I
    nevmax=sym ? n-1 : n-2
    if nevmax <= 0
        throw(ArgumentError("Input matrix A is too small. Use eigfact instead."))
    end
    if nev > nevmax
        @warn("Adjusting nev from $nev to $nevmax")
        nev = nevmax
    end
    if nev <= 0
        throw(ArgumentError("requested number of eigenvalues (nev) must be ≥ 1, got $nev"))
    end
    ncvmin = nev + (sym ? 1 : 2)
    if ncv < ncvmin
        @warn("Adjusting ncv from $ncv to $ncvmin")
        ncv = ncvmin
    end
    ncv = BlasInt(min(ncv, n))
    bmat = isgeneral ? "G" : "I"
    isshift = sigma !== nothing

    if isa(which,AbstractString)
        @warn("Use symbols instead of strings for specifying which eigenvalues to compute")
        which=Symbol(which)
    end
    if (which != :LM && which != :SM && which != :LR && which != :SR &&
        which != :LI && which != :SI && which != :BE)
        throw(ArgumentError("which must be :LM, :SM, :LR, :SR, :LI, :SI, or :BE, got $(repr(which))"))
    end
    if which == :BE && !sym
        throw(ArgumentError("which=:BE only possible for real symmetric problem"))
    end
    isshift && which == :SM && @warn("use of :SM in shift-and-invert mode is not recommended, use :LM to find eigenvalues closest to sigma")

    if which==:SM && !isshift # transform into shift-and-invert method with sigma = 0
        isshift=true
        sigma=zero(T)
        which=:LM
    end

    if sigma !== nothing && !iscmplx && isa(sigma,Complex)
        throw(ArgumentError("complex shifts for real problems are not yet supported"))
    end
    sigma = isshift ? convert(T,sigma) : zero(T)

    if !isempty(v0)
        if length(v0) != n
            throw(DimensionMismatch())
        end
        if eltype(v0) != T
            throw(ArgumentError("starting vector must have element type $T, got $(eltype(v0))"))
        end
    end

    whichstr = "LM"
    if which == :BE
        whichstr = "BE"
    end
    if which == :LR
        whichstr = (!sym ? "LR" : "LA")
    end
    if which == :SR
        whichstr = (!sym ? "SR" : "SA")
    end
    if which == :LI
        if !sym
            whichstr = "LI"
        else
            throw(ArgumentError("largest imaginary is meaningless for symmetric eigenvalue problems"))
        end
    end
    if which == :SI
        if !sym
            whichstr = "SI"
        else
            throw(ArgumentError("smallest imaginary is meaningless for symmetric eigenvalue problems"))
        end
    end

    # Refer to ex-*.doc files in ARPACK/DOCUMENTS for calling sequence

    if isshift && !methodexists(factorize,(typeof(A),))
        throw(ArgumentError("shift mode only implemented for factorable matrices"))
    end

    matvecA!(y, x) = mul!(y, A, x)

    if !isgeneral           # Standard problem
        matvecB = x -> x
        if !isshift         #    Regular mode
            mode       = 1
            solveSI = x->x
        else                #    Shift-invert mode
            mode       = 3
            F = factorize(A - UniformScaling(sigma))
            solveSI = x -> F \ x
        end
    else                    # Generalized eigenproblem
        matvecB = x -> B * x
        if !isshift         #    Regular inverse mode
            mode       = 2
            F = factorize(B)
            solveSI = x -> F \ x
        else                #    Shift-invert mode
            mode       = 3
            F = factorize(A - sigma*B)
            solveSI = x -> F \ x
        end
    end

    # Compute the Ritz values and Ritz vectors
    (resid, v, ldv, iparam, ipntr, workd, workl, lworkl, rwork, TOL) =
       xaupd_wrapper(T, matvecA!, matvecB, solveSI, n, sym, iscmplx, bmat, nev, ncv, whichstr, tol, maxiter, mode, v0, options, chnl)

    V = copy(v)

    # Postprocessing to get eigenvalues and eigenvectors
    output = eupd_wrapper(T, n, sym, iscmplx, bmat, nev, whichstr, ritzvec, TOL,
                                 resid, ncv, v, ldv, sigma, iparam, ipntr, workd, workl, lworkl, rwork)

    # Issue 10495, 10701: Check that all eigenvalues are converged
    nev = length(output[1])
    nconv = output[ritzvec ? 3 : 2]
    nev ≤ nconv || @warn("not all wanted Ritz pairs converged. Requested: $nev, converged: $nconv")

    if wantH
        H = extract_mtx(matvecA!, V,ipntr,workl,ncv,iscmplx,sym)
        output = (output...,H,V)
    end
    output = (:finale, output...)

    if chnl == nothing
        return (output[2:end]...,)
    else
        put!(chnl,output)
    end
end # xeigs

function extract_mtx(matvecA!, V,ipntr,workl,ncv,iscmplx,sym)
    ih = ipntr[5]

    if !iscmplx
        H = zeros(ncv,ncv)
        if sym
            H = (diagm(workl[ih+1:ih-1+ncv],1)
                 + diagm(workl[ih+ncv:ih-1+2*ncv],0)
                 + diagm(workl[ih+1:ih-1+ncv],-1))
        else
            H[:] = workl[ih:ih-1+ncv^2]
        end
    else
        H = zeros(ncv,ncv) .+ 0im
        # CHECKME: workl is complex, but does ARPACK set ih for a real array?
        H[:] = workl[ih:ih-1+ncv^2]
    end

    if !iscmplx && sym
        β = workl[ih]
    else
        w = zeros(size(V,1))
        matvecA!(w,V[:,end])
        h = (w' * V)'
        mtxprod = w - V*h
        β = norm(mtxprod)
    end

    # form a rectangular Hessenberg matrix
    H = vcat(H,hcat(zeros(1,ncv-1),β))
end

"""
modified version of aupd_wrapper with optional plotting of intermediate results.
"""
function xaupd_wrapper(T, matvecA!::Function, matvecB::Function, solveSI::Function, n::Integer,
                      sym::Bool, cmplx::Bool, bmat::String,
                      nev::Integer, ncv::Integer, which::String,
                      tol::Real, maxiter::Integer, mode::Integer, v0::Vector,
                       options::Dict,chnl)

    lworkl = cmplx ? ncv * (3*ncv + 5) : (sym ? ncv * (ncv + 8) :  ncv * (3*ncv + 6) )
    TR = cmplx ? T.types[1] : T
    TOL = Ref{TR}(tol)

    v     = Matrix{T}(undef, n, ncv)
    workd = Vector{T}(undef, 3*n)
    workl = Vector{T}(undef, lworkl)
    rwork = cmplx ? Vector{TR}(undef, ncv) : Vector{TR}()

    if isempty(v0)
        resid = Vector{T}(undef, n)
        info  = Ref{BlasInt}(0)
    else
        resid = deepcopy(v0)
        info  = Ref{BlasInt}(1)
    end
    iparam = zeros(BlasInt, 11)
    ipntr  = zeros(BlasInt, (sym && !cmplx) ? 11 : 14)
    ido    = Ref{BlasInt}(0)

    iparam[1] = BlasInt(1)       # ishifts
    iparam[3] = BlasInt(maxiter) # maxiter
    iparam[7] = BlasInt(mode)    # mode

    zernm1 = 0:(n-1)

    all_shifts = []
    done_one_already = false
    H = []
    V = []
    iterx = 0

    while true
        if cmplx
            naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, rwork, info)
        elseif sym
            saupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, info)
        else
            naupd(ido, bmat, n, which, nev, TOL, resid, ncv, v, n,
                  iparam, ipntr, workd, workl, lworkl, info)
        end
        if info[] == 1
            @warn("incomplete convergence. Try a different"
                 * " starting vector or increase maxiter or ncv.")
            break
#            throw(ARPACKException("incomplete convergence. Try a different"
#                        * " starting vector or increase maxiter or ncv."))
        elseif info[] != 0
            throw(ARPACKException(info[]))
        end

        x = view(workd, ipntr[1] .+ zernm1)
        y = view(workd, ipntr[2] .+ zernm1)
        if mode == 1  # corresponds to dsdrv1, dndrv1 or zndrv1
            if ido[] == 1
                matvecA!(y, x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        elseif mode == 3 && bmat == "I" # corresponds to dsdrv2, dndrv2 or zndrv2
            if ido[] == -1 || ido[] == 1
                y[:] = solveSI(x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        elseif mode == 2 # corresponds to dsdrv3, dndrv3 or zndrv3
            if ido[] == -1 || ido[] == 1
                matvecA!(y, x)
                if sym
                    x[:] = y    # overwrite as per Remark 5 in dsaupd.f
                end
                y[:] = solveSI(y)
            elseif ido[] == 2
                y[:] = matvecB(x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        elseif mode == 3 && bmat == "G" # corresponds to dsdrv4, dndrv4 or zndrv4
            if ido[] == -1
                y[:] = solveSI(matvecB(x))
            elseif  ido[] == 1
                y[:] = solveSI(view(workd,ipntr[3] .+ zernm1))
            elseif ido[] == 2
                y[:] = matvecB(x)
            elseif ido[] == 99
                break
            else
                throw(ARPACKException("unexpected behavior"))
            end
        else
            throw(ArgumentError("ARPACK mode ($mode) not yet supported"))
        end

        # This is just a call counter since
        # standard ARPACK doesn't expose the real iteration nr
        iterx += 1
        # TODO: displayRitzValues

        if done_one_already
            # arpackgui_update
            if !cmplx
                if sym
                    dispvec = workl[ipntr[6] .+ (0:ncv-1)]
                    if which == :BE
                        # roughly nev large and small ews
                        the_ews = dispvec[end-2*nev+1:end]
                        the_shifts = dispvec[1:end-2*nev]
                    else
                        # nev eigenvalues
                        the_ews = dispvec[end-nev+1:end]
                        the_shifts = dispvec[1:end-nev]
                    end
                else # non-sym real
                # println("ew/shifts:",workl[ipntr[6]+(0:2*ncv-1)])
                    dispvec = (workl[ipntr[6] .+ (0:ncv-1)]
                               + workl[ipntr[7] .+ (0:ncv-1)]*im)
#                               + workl[ipntr[6]+(ncv:2*ncv-1)]*im)
                end
                # nev+1 ews (keep conjugate pairs together)
                the_ews = dispvec[end-nev:end]
                the_shifts = dispvec[1:end-nev-1]
            else # complex
                dispvec = workl[ipntr[6] .+ (0:ncv-1)]
                the_ews = dispvec[end-nev+1:end]
                the_shifts = dispvec[1:end-nev]
            end
            if isempty(all_shifts)
                all_shifts = copy(the_shifts)
            else
                append!(all_shifts,the_shifts)
            end
            # we want to use the previous H for PSA/RVec computation
            # since the current one is computed after the implicit restart
            # while the R vals and shifts apply before it.
            prevH = H
            prevV = V
            isempty(V) && (V = copy(v))
            H = extract_mtx(matvecA!, V,ipntr,workl,ncv,cmplx,sym)
            V = copy(v)

            # DEVNOTE: upstream has logic for analyzing a Ritz value here
            # if user invoked "pause"

            if ((chnl != nothing) && get(options,:ARPACK_plots,true)
                && (get(options,:ProgRVals,true) || get(options,:ProgPSA,false))
                && (mod(iterx, get(options,:ARPACK_plotinterval,ncv-nev)) == 0))
                # we only want to update about once per restart, but ARPACK
                # doesn't report restarts, so default to approximate
                # interval (ncv-nev)

                # DEVNOTE: upstream plots Pseudospectra
                # here if options[:ProgPSA] is set

                the_str = @sprintf("dim = %d, iter %d",n,iterx/(ncv-nev))
                stuff = (:media,dispvec,the_str,the_ews,
                         get(options,:ProgAllShifts,false) ? all_shifts
                                                           : the_shifts)
                put!(chnl,stuff)
            end
        end #arpackgui_update
        done_one_already = true

    end # iteration loop

    return (resid, v, n, iparam, ipntr, workd, workl, lworkl, rwork, TOL)
end


text:
using HierarchicalUtils: encode, stringify

"""
    AbstractMillModel

Supertype for any model defined in `Mill.jl`.
"""
abstract type AbstractMillModel end

const MillFunction = Union{Dense, Chain, Function}

_make_array_model(m::Union{MillFunction, Function}) = ArrayModel(m)
_make_array_model(m) = m

include("arraymodel.jl")
include("bagmodel.jl")
include("productmodel.jl")
include("lazymodel.jl")
include("reflectinmodel.jl")

ext: jl


text:
function pmap(pool::WorkerPool, f, c; output::IO=STDERR, dt::Real=0.5)
    meters = Array{Progress}(nworkers(pool))
    channels = [RemoteChannel(() -> Channel{Any}(16), w) for w in 1:nworkers()]
    remaining = length(c)

    g(x) = f(x, channels[myid() - 1])
    g = remote(pool, g)

    h(x) = (remaining -= 1; g(x))

    done = false
    lastPrint = time()

    #clear space for the progress meters
    print(output, "\n"^(nworkers(pool) + 2))

    @async while true
        for (i, channel) in enumerate(channels)
            if isready(channel)
                update = take!(channel)

                if typeof(update) == Tuple{String, Int64}
                    meters[i] = Progress(update[2], desc=update[1], do_not_print=true)
                elseif update < 0
                    finish!(meters[i])
                else
                    update!(meters[i], update)
                end
            end
        end

        t = time()
        if t > lastPrint + dt
            move_cursor_up_while_clearing_lines(output, nworkers() + 2)

            println(output, remaining, " tasks remaining to start")
            println(output, "-------------------------------")
            for i in 1:length(meters)
                if isassigned(meters, i)
                    printover(output, progress_meter_string(meters[i])*"\n", meters[i].color)
                else
                    println(output)
                end
            end

            lastPrint = t
        end

        yield()

        done && break
    end

    results = asyncmap(h, c, ntasks = () -> nworkers(pool))

    done = true
    sleep(0.1)

    return results
end
pmap(f, c) = pmap(default_worker_pool(), f, c)

macro info(str)
    if myid() == 1
        esc(:(Base.print_with_color(:cyan, "INFO: ", $str, "\n")))
    end
end

macro parallelprogress(args...)
    if length(args) < 2
        throw(ArgumentError("@paralleprogress requires at least two arguments"))
    end
    metersym = args[1]
    progressargs = args[2:end-1]
    loop = args[end]

    if loop.head != :for
        throw(ArgumentError("@paralleprogress only works on for loops"))
    end

    progressString = "Progress: "
    if length(progressargs) > 0
        progressString = progressargs[1]
    end

    loopIttr = loop.args[1].args[2]
    loopVar = loop.args[1].args[1]
    loopBody = loop.args[2].args
    loopLength = :( length($loopIttr) )

    if myid() == 1
        setup = esc(:($metersym = ProgressMeter.Progress($loopLength, desc=$progressString) ))
        update = :(ProgressMeter.update!($metersym, $loopVar))
        finish = esc(:(ProgressMeter.finish!($metersym)))
    else
        setup = esc(:(put!($metersym, ($progressString, $loopLength)); lastUpdateTime = time() ))
        update = quote
                    t = time()
                    if t > lastUpdateTime + 0.5
                        @async put!($metersym, $loopVar)
                        lastUpdateTime = t
                    end

                    yield()
                end
        finish = esc(:(put!($metersym, -1)))
    end

    push!(loop.args[2].args, update)
    loop = esc(loop)

    return quote
        $setup
        $loop
        $finish
    end
end


text:
abstract type Relocation{T<:ObjectHandle} end
abstract type RelocationRef{T<:ObjectHandle} end

struct LOIByName
    addrs::Dict{Symbol, UInt64}
end
function getSectionLoadAddress(LOI::LOIByName, x::Union{Symbol, AbstractString})
    return LOI.addrs[symbol(x)]
end
function getSectionLoadAddress(LOI::LOIByName, sec)
    return getSectionLoadAddress(LOI, bytestring(sectionname(sec)))
end


text:
"""
Base type for an MDP/POMDP solver
"""
abstract type Solver end

"""
    solve(solver::Solver, problem::POMDP)

Solves the POMDP using method associated with solver, and returns a policy.
"""
function solve end


text:
@testset "accessors" begin

    amb = LaxZonedDateTime(DateTime(2015,11,1,1),winnipeg)
    amb_first = LaxZonedDateTime(ZonedDateTime(2015,11,1,1,winnipeg,1))
    amb_last = LaxZonedDateTime(ZonedDateTime(2015,11,1,1,winnipeg,2))

    @test amb_first < amb_last
    @test !(amb_first < amb)
    @test !(amb < amb_first)
    @test !(amb_last < amb)
    @test !(amb < amb_last)

    non_existent = LaxZonedDateTime(DateTime(2015,3,8,2),winnipeg)
    @test ZonedDateTime(2015,3,8,1,winnipeg) < non_existent < ZonedDateTime(2015,3,8,3,winnipeg)

    @test hour(amb) == hour(amb_first) == hour(amb_last) == 1
    @test hour(non_existent) == 2

    @test month(amb) == 11

    null = LaxZonedDateTime()
    @test_throws Exception DateTime(null)
    @test_throws Exception DateTime(null, UTC)
    @test_throws Exception hour(null)

    @test !isvalid(null)
    @test !isvalid(amb)
    @test isvalid(amb_first)
    @test isvalid(amb_last)
    @test !isvalid(non_existent)

    @test !isinvalid(null)
    @test isinvalid(amb)
    @test !isinvalid(amb_first)
    @test !isinvalid(amb_last)
    @test isinvalid(non_existent)

    @test !isrepresentable(null)
    @test isrepresentable(amb)
    @test isrepresentable(amb_first)
    @test isrepresentable(amb_last)
    @test isrepresentable(non_existent)

    @test !isambiguous(null)
    @test isambiguous(amb)
    @test !isambiguous(amb_first)
    @test !isambiguous(amb_last)
    @test !isambiguous(non_existent)

    @test !isnonexistent(null)
    @test !isnonexistent(amb)
    @test !isnonexistent(amb_first)
    @test !isnonexistent(amb_last)
    @test isnonexistent(non_existent)

    @test amb_last - amb_first == Hour(1)
    @test amb - amb_first === nothing
    @test non_existent - amb_first === nothing
    @test null - amb_first === nothing

end


text:
using Plots
default(legend = false)
x = y = range(-5, 5, length = 40)
zs = zeros(0, 40)
n = 100

@gif for i in range(0, stop = 2π, length = n)
    f(x, y) = sin(x + 10sin(i)) + cos(y)

    # create a plot with 3 subplots and a custom layout
    l = @layout [a{0.7w} b; c{0.2h}]
    p = plot(x, y, f, st = [:surface, :contourf], layout = l)

    # induce a slight oscillating camera angle sweep, in degrees (azimuth, altitude)
    plot!(p[1], camera = (10 * (1 + cos(i)), 40))

    # add a tracking line
    fixed_x = zeros(40)
    z = map(f, fixed_x, y)
    plot!(p[1], fixed_x, y, z, line = (:black, 5, 0.2))
    vline!(p[2], [0], line = (:black, 5))

    # add to and show the tracked values over time
    global zs = vcat(zs, z')
    plot!(p[3], zs, alpha = 0.2, palette = cgrad(:blues).colors)
end

max_stars_repo_path: julia/sine and cosine wave.jl

max_stars_repo_name: abhishek-parashar/Right-From-Scratch

max_issues_repo_path: julia/sine and cosine wave.jl

max_issues_repo_name: abhishek-parashar/Right-From-Scratch

max_forks_repo_path: julia/sine and cosine wave.jl

max_forks_repo_name: abhishek-parashar/Right-From-Scratch


text:
# Copyright (c) 2020-2022, Bank of Canada
# All rights reserved.

using TimeSeriesEcon

"""
    readfame(db, args...;
        namecase=lowercase,
        prefix=nothing, glue="_",
        collect=[],
        wc_options...)

Read data from FAME database into Julia. The data is returned in a
[`Workspace`](@ref TimeSeriesEcon.Workspace).

`db` is a [`FameDatabase`](@ref) or a `String`. If `db` is a `String`, the
database will be opened in `:readonly` mode and closed after loading the data.

If `db` is the only argument, then all objects in the database will be loaded.
Arguments and options can be used to restrict which objects will be loaded.

Each positional argument after the first one should be a string or a Symbol.
* If it is a string that contains a wildcard character (`'?'` or `'^'`, see Fame
  help about wildcards) then we call [`listdb`](@ref) to obtain a list of
  objects matching the given wildcard. In this case, we pass the given
  `wc_options...` to [`listdb`](@ref). You can use them to limit the wildcard
  search to specific class, type, frequency, etc. See [`listdb`](@ref) for
  details.
* Otherwise (Symbol or a string not containing wildcard characters), an object
  with the given name will be loaded. `wc_options...` are ignored, meaning that
  the object will be loaded no matter its class, type, frequency, etc.

The following options can be used to modify how the Julia identifiers are
constructed from the FAME names.

* `namecase` - FAME identifiers are case-insensitive and FAME always returns
  them in upper case. By default we convert the names to lower case. You can
  pass any function that takes a string and returns a string as the value of the
  `namecase` option. The default is `lowercase`. For example, this may be a good
  idea if your database contains names with symbols that are not allowed in
  Julia identifiers. Currently, we just call `Symbol(namecase(fo.name))`, but you
  may want to substitute such symbols for something else, e.g.
  `namecase=(x->lowercase(replace(x, "@"=>"_")))`.
* `prefix`, if given, will be stripped (together with `glue`) from the beginning
  of the FAME name. If the name does not begin with the given `prefix` then it
  will remain unchanged. If you want to load only names starting with the given
  prefix you must use the appropriate wildcard. Default is `prefix=nothing`,
  which disables this functionality. Note that `prefix=nothing` and `prefix=""`
  are not the same.
* `collect` can be a string/symbol, a vector whose elements are strings/symbols
  or vectors of strings/symbols, etc. If the name begins with one of the given
  `collect` values (together with the `glue`), then the object will be loaded in
  a nested [`Workspace`](@ref). The idea is that nested Workspaces written to
  the database would be loaded in the same structure by providing the list of
  names of the nested Workspaces in the `collect` option value. See the examples
  below.
* `glue` is used to join the `prefix` or the `collect` strings to the rest of
  the name. Use the same value as in [`writefame`](@ref) in order for this to
  work.

## Examples
```
julia> w = Workspace(; a = 1, b=TSeries(2020Q1, randn(10)),
       s = MVTSeries(2020M1, (:q, :p), randn(24,2)),
       c = Workspace(; alpha = 0.1, beta = 0.8,
       n = Workspace(; s = "Hello World")
       ))
Workspace with 4-variables
  a ⇒ 1
  b ⇒ 10-element TSeries{Quarterly} with range 2020Q1:2022Q2
  s ⇒ 24×2 MVTSeries{Monthly} with range 2020M1:2021M12 and variables (q,p)
  c ⇒ Workspace with 3-variables

julia> writefame("data.db", w); listdb("data.db")
7-element Vector{FameObject}:
 A: scalar,numeric,undefined,0:0,0:0
 B: series,precision,quarterly_december,2020:1,2022:2
 C_ALPHA: scalar,precision,undefined,0:0,0:0
 C_BETA: scalar,precision,undefined,0:0,0:0
 C_N_S: scalar,string,undefined,0:0,0:0
 S_P: series,precision,monthly,2020:1,2021:12
 S_Q: series,precision,monthly,2020:1,2021:12

julia> # read only variables in the list
julia> readfame("data.db", "a", "b")
Workspace with 2-variables
  a ⇒ 1.0
  b ⇒ 10-element TSeries{Quarterly} with range 2020Q1:2022Q2

julia> # read everything as it is in the database
julia> readfame("data.db")
Workspace with 7-variables
        a ⇒ 1.0
        b ⇒ 10-element TSeries{Quarterly} with range 2020Q1:2022Q2
  c_alpha ⇒ 0.1
   c_beta ⇒ 0.8
    c_n_s ⇒ 11-codeunit String
      s_p ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12
      s_q ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12

julia> # prefix is stripped where it appears (still loading everything)
julia> readfame("data.db", prefix="c")
Workspace with 7-variables
      a ⇒ 1.0
      b ⇒ 10-element TSeries{Quarterly} with range 2020Q1:2022Q2
  alpha ⇒ 0.1
   beta ⇒ 0.8
    n_s ⇒ 11-codeunit String
    s_p ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12
    s_q ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12

julia> # wildcard search, no prefix (name remains unchanged)
julia> readfame("data.db", "s?")
Workspace with 2-variables
  s_p ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12
  s_q ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12

julia> # prefix (stripped) with matching wildcard search
julia> readfame("data.db", "s?", prefix="s")
Workspace with 2-variables
  p ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12
  q ⇒ 24-element TSeries{Monthly} with range 2020M1:2021M12

julia> # collect with matching wildcard
julia> readfame("data.db", "c?", collect="c")
Workspace with 1-variables
  c ⇒ Workspace with 3-variables

julia> # nested collect - matches the original structure of nested Workspaces
julia> readfame("data.db", collect=["c"=>["n"], "s"])
Workspace with 4-variables
  a ⇒ 1.0
  b ⇒ 10-element TSeries{Quarterly} with range 2020Q1:2022Q2
  c ⇒ Workspace with 3-variables
  s ⇒ Workspace with 2-variables

```

"""
function readfame end
export readfame

# The workhorse - reads one FameObject and writes it into the
function _readfame!(ret, ident, fo, db)
    try
        ret[ident] = unfame(do_read!(fo, db))
    catch e
        @info "Didn't load $(fo.name): $(sprint(showerror, e))"
    end
    return nothing
end

# If db is given as a string
@inline readfame(dbname::AbstractString, args...; kwargs...) =
    opendb(dbname) do db
        readfame(db, args...; kwargs...)
    end
# if list of names is not given
readfame(db::FameDatabase; kwargs...) = readfame(db, "?"; kwargs...)
# The other cases
function readfame(db::FameDatabase, args...; namecase = lowercase, prefix = nothing, collect = [], glue = "_", kwargs...)
    if collect isa Union{AbstractString,Symbol,Pair{<:Union{AbstractString,Symbol},<:Any}}
        collect = [collect]
    end
    ret = Workspace()
    for wc in args
        fos = _iswildcard(wc) ? listdb(db, wc; kwargs...) : [quick_info(db, wc)]
        for fo in fos
            dest, name = _destination(ret, fo.name, glue, prefix, collect...)
            _readfame!(dest, Symbol(namecase(name)), fo, db)
        end
    end
    return ret
end

_iswildcard(wc) = occursin('?', wc) || occursin('^', wc)

_remove_prefix(name, pref) = startswith(name, pref) ? replace(name, pref => ""; count = 1) : name
# Handle the prefix argument
_destination(ret, name, glue, prefix::AbstractString, args...) = _destination(ret, _remove_prefix(name, uppercase(string(prefix) * glue)), glue, nothing, args...)
# Recursion on the collect arguments
_destination(ret, name, glue, ::Nothing) = (ret, name)
_destination(ret, name, glue, ::Nothing, W::Any, args...) = error("Invalid type of W: $(typeof(W))")
_destination(ret, name, glue, ::Nothing, W::Union{Symbol,AbstractString}, args...) = _destination(ret, name, glue, nothing, W => [], args...)
@inline function _destination(ret, name, glue, ::Nothing, (Wpref, Wargs)::Pair{<:Union{Symbol,AbstractString},<:Any}, args...)
    pref = uppercase(string(Wpref) * glue)
    if startswith(name, pref)
        return _destination(get!(ret, Symbol(Wpref), Workspace()), replace(name, pref => ""; count = 1), glue, nothing, Wargs...)
    else
        return _destination(ret, name, glue, nothing, args...)
    end
end

"""
    unfame(fo::FameObject)

Convert a `FameObject` to a Julia type.

* PRECISION SCALAR => `Float64`
* NUMERIC SCALAR => `Float32`
* STRING SCALAR => `String`
* BOOLEAN SCALAR => `Bool`
* NAMELIST => `String` (Formatted as "{NAME1,NAME2,ETC}", see Fame help for
  details.)
* DATE SCALAR => [`MIT`](@ref TimeSeriesEcon.MIT) (CASE becomes `MIT{Unit}`,
  Frequencies not supported by TimeSeriesEcon throw an `ErrorException`)
* PRECISION SERIES => [`TSeries`](@ref TimeSeriesEcon.TSeries)
* NUMERIC SERIES => [`TSeries`](@ref TimeSeriesEcon.TSeries) with element type
  `Float32`
* BOOLEAN SERIES => [`TSeries`](@ref TimeSeriesEcon.TSeries) with element type
  `Bool`
* DATE SERIES => [`TSeries`](@ref TimeSeriesEcon.TSeries) with element type
  [`MIT`](@ref TimeSeriesEcon.MIT)
* STRING SERIES => `Vector{String}` (the time series metadata is lost)
"""
function unfame end
export unfame

unfame(fo::FameObject{:scalar,:precision}) = _unmissing!(fo.data, Val(:precision))[]
unfame(fo::FameObject{:scalar,:numeric}) = _unmissing!(fo.data, Val(:numeric))[]
unfame(fo::FameObject{:scalar,:string}) = fo.data[]
unfame(fo::FameObject{:scalar,:boolean}) = fo.data[] != 0
unfame(fo::FameObject{:scalar,:namelist}) = fo.data[]
unfame(fo::FameObject{:scalar,TY}) where {TY} = _date_to_mit(TY, fo.data[])

unfame(fo::FameObject{:series,:precision,FR}) where {FR} =
    TSeries(_date_to_mit(FR, fo.first_index[]), _unmissing!(fo.data, Val(:precision)))
unfame(fo::FameObject{:series,:numeric,FR}) where {FR} =
    TSeries(_date_to_mit(FR, fo.first_index[]), _unmissing!(fo.data, Val(:numeric)))
unfame(fo::FameObject{:series,:boolean,FR}) where {FR} =
    TSeries(_date_to_mit(FR, fo.first_index[]), [d != 0 for d in fo.data])
unfame(fo::FameObject{:series,TY,FR}) where {TY,FR} =
    TSeries(_date_to_mit(FR, fo.first_index[]), [_date_to_mit(TY, d) for d in fo.data])
unfame(fo::FameObject{:series,:string,FR}) where {FR} =
    copy(fo.data) # error("Can't handle string series")

@inline _freq_from_fame(fr) =
    (
        fr = string(val_to_symbol(fr, fame_freq));
        startswith(fr, "quarterly") ? Quarterly :
        startswith(fr, "monthly") ? Monthly :
        startswith(fr, "annual") ? Yearly :
        startswith(fr, "case") ? Unit :
        error("Cannot convert FAME frequency $fr to TimeSeriesEcon.")
    )

function _date_to_mit(fr, date)
    F = _freq_from_fame(fr)
    F == Unit && return MIT{F}(date)
    y = Ref{Cint}()
    p = Ref{Cint}()
    if _ismissing(date, Val(fr))
        return typenan(MIT{F})
    end
    @fame_call_check(fame_index_to_year_period,
        (Cint, FameIndex, Ref{Cint}, Ref{Cint}),
        val_to_int(fr, fame_freq), date, y, p)
    return MIT{F}(y[], p[])
end

@inline _freq_to_fame(F) =
    (
        F == Unit ? :case :
        F == Quarterly ? :quarterly_december :
        F == Monthly ? :monthly :
        F == Yearly ? :annual_december :
        error("Cannot convert TimeSeriesEcon frequency $F to a FAME frequency.")
    )

function _mit_to_date(x::MIT{F}) where {F<:Frequency}
    F == Unit && return (:case, FameIndex(Int(x)))
    fr = _freq_to_fame(F)
    if istypenan(x)
        return (fr, FAME_INDEX_NC)
    end
    yr, p = mit2yp(x)
    index = Ref{FameIndex}(-1)
    @fame_call_check(fame_year_period_to_index,
        (Cint, Ref{FameIndex}, Cint, Cint),
        val_to_int(fr, fame_freq), index, yr, p)
    return (fr, index[])
end

"""
    writefame(db, data; options)

Write Julia data to a FAME database.

### Arguments:
* `db` can be a string containing the path to a FAME .db file or an instance of
  [`FameDatabase`](@ref).
* `data` is a collection of data which will be written to the database.
  - If `data` is an `MVTSeries`, each series will be written to the database.
  - If `data` is a `Workspace`, each element will be written as a separate FAME
    object. In this case any nested `Workspace` of `MVTSeries` objects will be
    written recursively by prepending the names.

### Options:
* `mode::Symbol` - if `db` is a string, the database is opened with the
  given `mode`. The default is :overwrite.
* `prefix::Union{Nothing,String}` - the given `prefix` will be prepended to the
  name of each object written in the database. The default is `nothing`, i.e.
  nothing will be prepended. NOTE: `prefix=nothing` and `prefix=""` are not the
  same.
* `glue::String` - the `glue` is used to join the prefix to the name. The
  default is `"_"`.

### Examples:
```
julia> w = Workspace(; a=1, b=TSeries(2020Q1,ones(10)), s=MVTSeries(2020M1, collect("ab"), randn(24,2)))
Workspace with 2-variables
  a ⇒ 1.0
  b ⇒ 10-element TSeries{Quarterly} with range 2020Q1:2022Q2
  s ⇒ 24×2 MVTSeries{Monthly} with range 2020M1:2021M12 and variables …

julia> writefame("w.db", w)

julia> listdb("w.db")
2-element Vector{FAME.FameObject}:
 A: scalar,precision,undefined,0:0,0:0
 B: series,precision,quarterly_december,2020:1,2022:2
 S_A: series,precision,monthly,2020:1,2021:12
 S_B: series,precision,monthly,2020:1,2021:12

```
"""
function writefame end
export writefame

# The workhorse - write an iterable of (name, value) pairs
function _writefame(db, iterable, prefix, glue)
    for (name, value) in iterable
        if prefix !== nothing
            name = Symbol(prefix, glue, name)
        end
        if value isa Union{Workspace,MVTSeries}
            _writefame(db, pairs(value), name, glue) #= glue= =#
        else
            try
                fo = refame(name, value)
                do_write(fo, db)
            catch e
                @info "Failed to write $(name): $(sprint(showerror, e))" e
                continue
            end
        end
    end
end

# If database is given as a string
@inline writefame(dbname::AbstractString, data...; mode = :overwrite, kwargs...) =
    opendb(dbname, mode) do db
        writefame(db, data...; kwargs...)
    end

const _FameWritable = Union{MVTSeries,Workspace}

# # write a single Workspace or MVTSeries
# @inline writefame(db::FameDatabase, data::_FameWritable; prefix=nothing, glue="_") =
#     _writefame(db, pairs(data), prefix, glue)

# write a list of MVTSeries and Workspace
writefame(db::FameDatabase, data::_FameWritable...; kwargs...) = writefame(db, data; kwargs...)
@inline function writefame(db::FameDatabase, data::Tuple{_FameWritable,Vararg{_FameWritable}}; prefix = nothing, glue = "_")
    for d in data
        _writefame(db, pairs(d), prefix, glue)
    end
end

"""
    refame(name, value)

Convert the given value to a [`FameObject`](@ref)

* `Real` => NUMERIC SCALAR (includes integers)
* `Float64` => PRECISION SCALAR
* `Bool` => BOOLEAN SCALAR
* [`MIT`](@ref TimeSeriesEcon.MIT) => DATE SCALAR
* `String` => NAMELIST, if it is in the form "{name1,name2,...}", otherwise
  STRING SCALAR
* `Vector{String}` => CASE SERIES of STRING
* [`TSeries`](@ref TimeSeriesEcon.TSeries) => a SERIES of the same frequency and
  type. The type conversions are the same as for scalars.
"""
function refame end
export refame

@inline refame(name::Symbol, value::T) where {T<:Real} = refame(name, convert(promote_type(T, Float32), value))

@inline function refame(name::Symbol, value::Float64)
    if value == typenan(value)
        value = FAME.FPRCNC
    end
    return FameObject{:scalar,:precision,:undefined,Ref{Float64}}(
        string(name), Ref(0), Ref(0), Ref(Float64(value)))
end

@inline function refame(name::Symbol, value::Float32)
    if value == typenan(value)
        value = FAME.FNUMNC
    end
    return FameObject{:scalar,:numeric,:undefined,Ref{Float32}}(
        string(name), Ref(0), Ref(0), Ref(value))
end

@inline function refame(name::Symbol, value::Bool)
    return FameObject{:scalar,:boolean,:undefined,Ref{Int32}}(
        string(name), Ref(0), Ref(0), Ref(Int32(value)))
end

@inline function refame(name::Symbol, value::MIT)
    (fr, ind) = _mit_to_date(value)
    return FameObject{:scalar,fr,:undefined,Ref{FameIndex}}(
        string(name), Ref(0), Ref(0), Ref(ind))
end

function refame(name::Symbol, value::String)
    ty = :string
    if length(value) > 1 && value[1] == '{' && value[end] == '}'
        ty = :namelist
    end
    return FameObject{:scalar,ty,:undefined,Ref{String}}(
        string(name), Ref(0), Ref(0), Ref(value))
end

function refame(name::Symbol, value::AbstractVector{<:AbstractString})
    return FameObject{:series,:string,:case,Vector{String}}(
        string(name), Ref(1), Ref(length(value)), String[value...]
    )
end

@inline refame(name::Symbol, value) =
    error("Cannot write type $(typeof(value)) to FAME database.")

function refame(name::Symbol, value::TSeries)
    if isempty(value)
        fr = _freq_to_fame(frequencyof(value))
        find = lind = FAME_INDEX_NC
    else
        (fr, find) = _mit_to_date(firstdate(value))
        lind = find + length(value) - 1
    end
    ElType = eltype(value)
    if ElType == Float32
        ty = :numeric
        val = [istypenan(v) ? FNUMNC : v for v in value.values]
    elseif ElType == Float64
        ty = :precision
        val = [istypenan(v) ? FPRCNC : v for v in value.values]
    elseif ElType <: MIT
        ty = _freq_to_fame(frequencyof(ElType))
        val = [_mit_to_date(x)[2] for x in value.values]
    elseif ElType == Bool
        ty = :boolean
        val = [Int32(x) for x in value.values]
    else
        ty = :precision
        val = map(value.values) do v
            v = convert(Float64, v)
            istypenan(v) ? FPRCNC : v
        end
    end
    return FameObject{:series,ty,fr,typeof(val)}(
        string(name), Ref(find), Ref(lind), val)
end

const HNMVAL = Int32(0) # Normal value; not missing or magic
const HNCVAL = Int32(1) # Missing NC - Not Computable
const HNAVAL = Int32(2) # Missing NA - Not Available
const HNDVAL = Int32(3) # Missing ND - Not Defined
const HMGVAL = Int32(4) # Magic value; for internal use only

function _ismissing(x, ::Val{:precision})
    mt = Ref{Int32}(-1)
    @cfm_call_check(cfmispm, (Cdouble, Ref{Cint}), Float64(x), mt)
    return mt[] != HNMVAL
end

function _ismissing(x, ::Val{:numeric})
    mt = Ref{Int32}(-1)
    @cfm_call_check(cfmisnm, (Cfloat, Ref{Cint}), Float32(x), mt)
    return mt[] != HNMVAL
end

function _ismissing(x, ::Val{:boolean})
    mt = Ref{Int32}(-1)
    @cfm_call_check(cfmisbm, (Cint, Ref{Cint}), Int32(x), mt)
    return mt[] != HNMVAL
end

function _ismissing(x, ::Val{:string})
    mt = Ref{Int32}(-1)
    @cfm_call_check(cfmissm, (Cstring, Ref{Cint}), String(strip(x, '\0')), mt)
    return mt[] != HNMVAL
end

function _ismissing(x, ::Val{FR}) where {FR}
    mt = Ref{Int32}(-1)
    @fame_call_check(fame_date_missing_type, (Clonglong, Ref{Cint}), Int64(x), mt)
    return mt[] != HNMVAL
end

function _unmissing!(x::Ref, v::Val{T}) where {T}
    nan = typenan(eltype(x))
    if _ismissing(x[], v)
        x[] = nan
    end
    return x
end

function _unmissing!(x::Vector, v::Val{T}) where {T}
    nan = typenan(eltype(x))
    for i = eachindex(x)
        if _ismissing(x[i], v)
            x[i] = nan
        end
    end
    return x
end


text:
# this file contains functions for glue between this package and POMDPs.jl

function reweight!(wm, m::POMDP, b, a, pm, o)
    for i in 1:n_particles(b)
        s = particle(b, i)
        if isterminal(m, s)
            wm[i] = 0.0
        else
            sp = pm[i]
            wm[i] = obs_weight(m, s, a, sp, o)
        end
    end
end

function predict!(pm, m::POMDP, b, a, rng)
    all_terminal = true
    for i in 1:n_particles(b)
        s = particle(b, i)
        if !isterminal(m, s)
            all_terminal = false
            sp = @gen(:sp)(m, s, a, rng)
            pm[i] = sp
        end
    end
    if all_terminal
        error("Particle filter update error: all states in the particle collection were terminal.")
    end
end

particle_memory(pmodel::POMDP) = statetype(pmodel)[]

# all of the lines with resample should be changed to rand(rng, d, 3) when that becomes part of the POMDPs.jl standard

function initialize_belief(up::BasicParticleFilter, b::AbstractParticleBelief)
    return resample(ImportanceResampler(up.n_init), b, up.rng) # resample is needed to control number of particles
end

function initialize_belief(up::BasicParticleFilter, b::ParticleCollection)
    if n_particles(b) == up.n_init
        return b
    else
        return resample(ImportanceResampler(up.n_init), b, up.rng) # resample is needed to control number of particles
    end
end

function initialize_belief(up::BasicParticleFilter, b::AbstractVector)
    pc = ParticleCollection(b)
    return resample(ImportanceResampler(up.n_init), pc, up.rng)
end

function initialize_belief(up::BasicParticleFilter, d::D) where D
    # using weighted iterator here is more likely to be order n than just calling rand() repeatedly
    # but, this implementation is problematic and may change in the future
    try
        if @implemented(support(::D)) &&
            @implemented(iterate(::typeof(support(d)))) &&
            @implemented(pdf(::D, ::typeof(first(support(d)))))
                S = typeof(first(support(d)))
                particles = S[]
                weights = Float64[]
                for (s, w) in weighted_iterator(d)
                    push!(particles, s)
                    push!(weights, w)
                end
                return resample(ImportanceResampler(up.n_init), WeightedParticleBelief(particles, weights), up.rng)
        end
    catch ex
        if ex isa MethodError
            @warn("""
                Suppressing MethodError in initialize_belief in ParticleFilters.jl. Please file an issue here:

                https://github.com/JuliaPOMDP/ParticleFilters.jl/issues/new

                The error was

                $(sprint(showerror, ex))
                  """, maxlog=1)
        else
            rethrow(ex)
        end
    end

    return ParticleCollection(collect(rand(up.rng, d) for i in 1:up.n_init))
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

# NOTE: this function is a hack for Zygote.jl compatbility; see lib/zygote.jl
function warn_kmm_julialib(B, ϵ)
  # warn user that closed-form solution does
  # not consider probability simplex constraints
  isinf(B) || @warn "B parameter ignored when optlib=JuliaLib"
  iszero(ϵ) || @warn "ϵ parameter ignored when optlib=JuliaLib"
end

function _kmm_ratios(K, κ, dre::KMM, optlib::Type{JuliaLib})
  # retrieve parameters
  @unpack B, ϵ = dre

  warn_kmm_julialib(B, ϵ) # warn ignored parameters

  # density ratio via solve
  K \ vec(κ)
end


text:
struct DroneEnvParams{T}
    m::T
    L::T
    Iₘ::Matrix{T}
    Iₘ_inv::Matrix{T}
    k_m::T
    k_F::T
    k_M::T
    ω_min::T
    ω_max::T
end

Base.show(io::IO, params::DroneEnvParams) = print(
    io,
    join(["$p=$(getfield(params, p))" for p in fieldnames(DroneEnvParams)], ","),
)

mutable struct DroneEnv{A,T,R<:AbstractRNG} <: AbstractEnv
    params::DroneEnvParams{T}
    action_space::A
    observation_space::Space{Vector{ClosedInterval{T}}}
    state::Vector{T}
    done::Bool
    t::Int
    dt::T
    δt::T
    track::Track # Need to adjust this to a 3D track with the ability to add obstacles
    rng::R
end

"""
    DroneEnv(;kwargs...)

# Keyword arguments
- `T = Float64`
- `m = 0.25, 			                        # Mass (kg)`
- `L = 0.1,                                     # Length of the moment arm from the body to a motor`
- `I_m = diagm([2.5e-4, 2.5e-4, 1.2e-3]),     	# Moment of inertia in body frame (kg m^2)`
- `k_m = 20,	 			                    # Motor constant (1/s)`
- `k_F = 6.11e-8,   	                        # Motor force constant (N/rpm²)`
- `k_M = 1.5e-9,                                # Motor moment constant (N/rpm²)`
- `ω_min = 1200.0,                              # Motor min rotation speed (rpm)`
- `ω_max = 7800.0,                              # Motor max rotation speed (rpm)`
- `dt = 0.1`,                                   # Time step between actions
- `δt = 0.01,                                   # Time step used for integration`
- `track = string(@__DIR__, "/car_racing_tracks/curve.csv"),   # Track to load`
- `lane_width = 5.0,                            # Track lane_width`
- `rng = Random.GLOBAL_RNG`
"""
function DroneEnv(;
    T = Float64,
    m = 0.25,
    L = 0.1,
    I_m = diagm([2.5e-4, 2.5e-4, 1.2e-3]),
    k_m = 20,
    k_F = 6.11e-8,
    k_M = 1.5e-9,
    ω_min = 1200.0,
    ω_max = 7800.0,
    dt = 0.1,
    δt = 0.01,
    track = string(@__DIR__, "/car_racing_tracks/curve.csv"),
    lane_width = 5.0,
    rng = Random.GLOBAL_RNG,
)

    params = DroneEnvParams{T}(
        m,
        L,
        I_m,
        I_m^(-1),
        k_m,
        k_F,
        k_M,
        ω_min,
        ω_max,
    )

    return DroneEnv(params, T=T, dt=dt, δt=δt, track=track, rng=rng)
end

"""
    DroneEnv(params::DroneEnvParams;kwargs...)

# Keyword arguments
- `T = Float64`
- `dt = 0.1`,                                   # Time step between actions
- `δt = 0.01,                                   # Time step used for integration`
- `track = string(@__DIR__, "/car_racing_tracks/curve.csv"),   # Track to load`
- `lane_width = 5.0,                            # Track lane_width`
- `rng = Random.GLOBAL_RNG`
"""
function DroneEnv(params::DroneEnvParams;
    T = Float64,
    dt = 0.1,
    δt = 0.01,
    track = string(@__DIR__, "/car_racing_tracks/curve.csv"),
    lane_width = 5.0,
    rng = Random.GLOBAL_RNG,
)

    action_space = ClosedInterval{Vector{T}}(
        [-1.0, -1.0, -1.0, -1.0],
        [ 1.0,  1.0,  1.0,  1.0],
        )

    observation_space = Space([
        -Inf..Inf,                          # X position
        -Inf..Inf,                          # Y position
        -Inf..Inf,                          # Z position
        -Inf..Inf,                          # Vx (translational velocity)
        -Inf..Inf,                          # Vy
        -Inf..Inf,                          # Vz
        -π..π,                              # Euler angle 1 (roll, about x, ϕ)
        -π/2..π/2,                          # Euler angle 2 (pitch, about y, θ)
        -π..π,                              # Euler angle 3 (yaw, about z, ψ)
        -Inf..Inf,                          # Body rate 1 (roll rate, p)
        -Inf..Inf,                          # Body rate 2 (pitch rate, q)
        -Inf..Inf,                          # Body rate 4 (yaw rate, r)
        0.0..Inf,                           # Motor 1 speed (ω₁)
        0.0..Inf,                           # Motor 2 speed (ω₂)
        0.0..Inf,                           # Motor 3 speed (ω₃)
        0.0..Inf,                           # Motor 4 speed (ω₄)
        ])
    state = zeros(T,length(observation_space))
    state[end-3:end] = ones(T, 4) * 1200.0

    env = DroneEnv(
        params,
        action_space,
        observation_space,
        state,
        false,
        0,
        dt,
        δt,
        Track(track, width=lane_width),
        rng,
        )

    reset!(env)
    env
end

DroneEnv{T}(; kwargs...) where {T} = DroneEnv(; T = T, kwargs...)

Random.seed!(env::DroneEnv, seed) = Random.seed!(env.rng, seed)

RLBase.action_space(env::DroneEnv) = env.action_space
RLBase.state_space(env::DroneEnv{T}) where {T} = env.observation_space
RLBase.is_terminated(env::DroneEnv) = env.done
RLBase.state(env::DroneEnv) = env.state

# function within_track(env::DroneEnv)
#     return within_track(env.track, env.state[1:2])
# end

"""
    reward(env::DroneEnv)

"""
function RLBase.reward(env::DroneEnv{T}) where {T}
    if any(isnan.(env.state))
        return -100000.0
    end
    rew = 0.0

    # Go towards a point
    desired_pos = [10.0, 10.0, 10.0]
    dist_from_desired = norm(desired_pos - env.state[1:3])
    rew += -dist_from_desired

    # Penalize for more time spent away from the point
    if dist_from_desired >= 5
        rew += -1
    end

    return rew
end

function RLBase.reset!(env::DroneEnv{A,T}) where {A,T}
    state = zeros(T,length(env.state))
    state[end-3:end] = ones(T, 4) * 1500.0
    env.state = state
    env.t = 0
    env.done = false
    nothing
end

function RLBase.reset!(env::DroneEnv{A,T}, state::Vector{T}) where {A,T}
    env.state = state
    env.t = 0
    env.done = false
    nothing
end

"""
    a::Vector{Float64}
    size = (4,)
    a[1] = Desired motor speed for motor 1 (ω₁, rpm)
    a[2] = Desired motor speed for motor 2 (rpm)
    a[3] = Desired motor speed for motor 3 (rpm)
    a[4] = Desired motor speed for motor 4 (rpm)
"""
function (env::DroneEnv{<:ClosedInterval})(a::Vector{Float64})
    a in env.action_space || error("Action is not in action space")
    length(a) == 4 || error("Incorrect action size")
    _step!(env, a)
end

function (env::DroneEnv{<:ClosedInterval})(a::Vector{Int})
    env(Float64.(a))
end

function (env::DroneEnv{<:ClosedInterval})(a::Matrix{Float64})
    size(a)[2] == 1 || error("Only implented for one step")
    env(vec(a))
end

"""
Quadrotor dynamics NOT considering complications due to varying thrust based
on angle of attack, blade flapping, or airflow disruption.

N. Michael, D. Mellinger, Q. Lindsey, V. Kumar, The GRASP Multiple Micro UAV Testbed
"""
function _step!(env::DroneEnv, a::Vector{Float64})
    env.t += 1
    r = env.state[1:3]     # Position
    r_dot = env.state[4:6]  # Velocity
    Φ = env.state[7:9]      # Euler angles
    ω = env.state[10:12]    # body rates
    Ω = env.state[13:16]    # motor speeds

    dt = env.dt
    δt = env.δt

    # Scaled input controls from -1 to 1. Tis adjusts
    #  them to the appropriate rate (-1 = min, +1 = max)
    Ωᵈ = env.params.ω_min .+ (a .+1 ) .* ((env.params.ω_max-env.params.ω_min)/2)

    integration_steps = round(Int, dt/δt)
    for _ in 1:integration_steps

        Ω_dot = env.params.k_m .* (Ωᵈ - Ω)
        Ω += Ω_dot .* δt # Updated motor speeds

        F = env.params.k_F .* Ω.^2
        M = env.params.k_M .* Ω.^2

        temp_vec = [
            env.params.L * (F[2] - F[4]),
            env.params.L * (F[3] - F[1]),
            M[1] - M[2] + M[3] - M[4],
        ]
        ω_dot = env.params.Iₘ_inv * (temp_vec - cross(ω, env.params.Iₘ*ω))

        ω += ω_dot .* δt # Updated body rates

        temp_mat = [
            cos(Φ[2])  0 -cos(Φ[1])*sin(Φ[2]);
                0      1      sin(Φ[1])      ;
            sin(Φ[2])  0  cos(Φ[1])*cos(Φ[2])
        ]
        Φ_dot = temp_mat^(-1) * ω
        Φ += Φ_dot .* δt # Updated Euler angles
        Φ = mod.(Φ, 2π)

        R_11 = cos(Φ[3])*cos(Φ[2]) - sin(Φ[1])*sin(Φ[3])*sin(Φ[2])
        R_12 = -cos(Φ[1])*sin(Φ[3])
        R_13 = cos(Φ[3])*sin(Φ[2]) + cos(Φ[2])*sin(Φ[1])*sin(Φ[3])

        R_21 = cos(Φ[2])*sin(Φ[3]) + cos(Φ[3])*sin(Φ[1])*sin(Φ[2])
        R_22 = cos(Φ[1])*cos(Φ[3])
        R_23 = sin(Φ[3])*sin(Φ[2]) - cos(Φ[3])*cos(Φ[2])*sin(Φ[1])

        R_31 = -cos(Φ[1])*sin(Φ[2])
        R_32 = sin(Φ[1])
        R_33 = cos(Φ[1])*cos(Φ[2])

        R = [
             R_11 R_12 R_13;
             R_21 R_22 R_23;
             R_31 R_32 R_33
        ]

        r_ddot = [0,0,-9.81] + R*[0,0,sum(F)/env.params.m]
        r_dot += r_ddot .* δt
        r += r_dot .* δt
    end

    env.state[1:3] = r
    env.state[4:6] = r_dot
    env.state[7:9] = Φ
    env.state[10:12] = ω
    env.state[13:16] = Ω
    return env
end

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
const BOX_DRAWINGS_LIGHT_HORIZONTAL = '─'
const BOX_DRAWINGS_HEAVY_HORIZONTAL = '━'
const BOX_DRAWINGS_LIGHT_VERTICAL = '│'
const BOX_DRAWINGS_HEAVY_VERTICAL = '┃'
const BOX_DRAWINGS_LIGHT_TRIPLE_DASH_HORIZONTAL = '┄'
const BOX_DRAWINGS_HEAVY_TRIPLE_DASH_HORIZONTAL = '┅'
const BOX_DRAWINGS_LIGHT_TRIPLE_DASH_VERTICAL = '┆'
const BOX_DRAWINGS_HEAVY_TRIPLE_DASH_VERTICAL = '┇'
const BOX_DRAWINGS_LIGHT_QUADRUPLE_DASH_HORIZONTAL = '┈'
const BOX_DRAWINGS_HEAVY_QUADRUPLE_DASH_HORIZONTAL = '┉'
const BOX_DRAWINGS_LIGHT_QUADRUPLE_DASH_VERTICAL = '┊'
const BOX_DRAWINGS_HEAVY_QUADRUPLE_DASH_VERTICAL = '┋'
const BOX_DRAWINGS_LIGHT_DOWN_AND_RIGHT = '┌'
const BOX_DRAWINGS_DOWN_LIGHT_AND_RIGHT_HEAVY = '┍'
const BOX_DRAWINGS_DOWN_HEAVY_AND_RIGHT_LIGHT = '┎'
const BOX_DRAWINGS_HEAVY_DOWN_AND_RIGHT = '┏'
const BOX_DRAWINGS_LIGHT_DOWN_AND_LEFT = '┐'
const BOX_DRAWINGS_DOWN_LIGHT_AND_LEFT_HEAVY = '┑'
const BOX_DRAWINGS_DOWN_HEAVY_AND_LEFT_LIGHT = '┒'
const BOX_DRAWINGS_HEAVY_DOWN_AND_LEFT = '┓'
const BOX_DRAWINGS_LIGHT_UP_AND_RIGHT = '└'
const BOX_DRAWINGS_UP_LIGHT_AND_RIGHT_HEAVY = '┕'
const BOX_DRAWINGS_UP_HEAVY_AND_RIGHT_LIGHT = '┖'
const BOX_DRAWINGS_HEAVY_UP_AND_RIGHT = '┗'
const BOX_DRAWINGS_LIGHT_UP_AND_LEFT = '┘'
const BOX_DRAWINGS_UP_LIGHT_AND_LEFT_HEAVY = '┙'
const BOX_DRAWINGS_UP_HEAVY_AND_LEFT_LIGHT = '┚'
const BOX_DRAWINGS_HEAVY_UP_AND_LEFT = '┛'
const BOX_DRAWINGS_LIGHT_VERTICAL_AND_RIGHT = '├'
const BOX_DRAWINGS_VERTICAL_LIGHT_AND_RIGHT_HEAVY = '┝'
const BOX_DRAWINGS_UP_HEAVY_AND_RIGHT_DOWN_LIGHT = '┞'
const BOX_DRAWINGS_DOWN_HEAVY_AND_RIGHT_UP_LIGHT = '┟'
const BOX_DRAWINGS_VERTICAL_HEAVY_AND_RIGHT_LIGHT = '┠'
const BOX_DRAWINGS_DOWN_LIGHT_AND_RIGHT_UP_HEAVY = '┡'
const BOX_DRAWINGS_UP_LIGHT_AND_RIGHT_DOWN_HEAVY = '┢'
const BOX_DRAWINGS_HEAVY_VERTICAL_AND_RIGHT = '┣'
const BOX_DRAWINGS_LIGHT_VERTICAL_AND_LEFT = '┤'
const BOX_DRAWINGS_VERTICAL_LIGHT_AND_LEFT_HEAVY = '┥'
const BOX_DRAWINGS_UP_HEAVY_AND_LEFT_DOWN_LIGHT = '┦'
const BOX_DRAWINGS_DOWN_HEAVY_AND_LEFT_UP_LIGHT = '┧'
const BOX_DRAWINGS_VERTICAL_HEAVY_AND_LEFT_LIGHT = '┨'
const BOX_DRAWINGS_DOWN_LIGHT_AND_LEFT_UP_HEAVY = '┩'
const BOX_DRAWINGS_UP_LIGHT_AND_LEFT_DOWN_HEAVY = '┪'
const BOX_DRAWINGS_HEAVY_VERTICAL_AND_LEFT = '┫'
const BOX_DRAWINGS_LIGHT_DOWN_AND_HORIZONTAL = '┬'
const BOX_DRAWINGS_LEFT_HEAVY_AND_RIGHT_DOWN_LIGHT = '┭'
const BOX_DRAWINGS_RIGHT_HEAVY_AND_LEFT_DOWN_LIGHT = '┮'
const BOX_DRAWINGS_DOWN_LIGHT_AND_HORIZONTAL_HEAVY = '┯'
const BOX_DRAWINGS_DOWN_HEAVY_AND_HORIZONTAL_LIGHT = '┰'
const BOX_DRAWINGS_RIGHT_LIGHT_AND_LEFT_DOWN_HEAVY = '┱'
const BOX_DRAWINGS_LEFT_LIGHT_AND_RIGHT_DOWN_HEAVY = '┲'
const BOX_DRAWINGS_HEAVY_DOWN_AND_HORIZONTAL = '┳'
const BOX_DRAWINGS_LIGHT_UP_AND_HORIZONTAL = '┴'
const BOX_DRAWINGS_LEFT_HEAVY_AND_RIGHT_UP_LIGHT = '┵'
const BOX_DRAWINGS_RIGHT_HEAVY_AND_LEFT_UP_LIGHT = '┶'
const BOX_DRAWINGS_UP_LIGHT_AND_HORIZONTAL_HEAVY = '┷'
const BOX_DRAWINGS_UP_HEAVY_AND_HORIZONTAL_LIGHT = '┸'
const BOX_DRAWINGS_RIGHT_LIGHT_AND_LEFT_UP_HEAVY = '┹'
const BOX_DRAWINGS_LEFT_LIGHT_AND_RIGHT_UP_HEAVY = '┺'
const BOX_DRAWINGS_HEAVY_UP_AND_HORIZONTAL = '┻'
const BOX_DRAWINGS_LIGHT_VERTICAL_AND_HORIZONTAL = '┼'
const BOX_DRAWINGS_LEFT_HEAVY_AND_RIGHT_VERTICAL_LIGHT = '┽'
const BOX_DRAWINGS_RIGHT_HEAVY_AND_LEFT_VERTICAL_LIGHT = '┾'
const BOX_DRAWINGS_VERTICAL_LIGHT_AND_HORIZONTAL_HEAVY = '┿'
const BOX_DRAWINGS_UP_HEAVY_AND_DOWN_HORIZONTAL_LIGHT = '╀'
const BOX_DRAWINGS_DOWN_HEAVY_AND_UP_HORIZONTAL_LIGHT = '╁'
const BOX_DRAWINGS_VERTICAL_HEAVY_AND_HORIZONTAL_LIGHT = '╂'
const BOX_DRAWINGS_LEFT_UP_HEAVY_AND_RIGHT_DOWN_LIGHT = '╃'
const BOX_DRAWINGS_RIGHT_UP_HEAVY_AND_LEFT_DOWN_LIGHT = '╄'
const BOX_DRAWINGS_LEFT_DOWN_HEAVY_AND_RIGHT_UP_LIGHT = '╅'
const BOX_DRAWINGS_RIGHT_DOWN_HEAVY_AND_LEFT_UP_LIGHT = '╆'
const BOX_DRAWINGS_DOWN_LIGHT_AND_UP_HORIZONTAL_HEAVY = '╇'
const BOX_DRAWINGS_UP_LIGHT_AND_DOWN_HORIZONTAL_HEAVY = '╈'
const BOX_DRAWINGS_RIGHT_LIGHT_AND_LEFT_VERTICAL_HEAVY = '╉'
const BOX_DRAWINGS_LEFT_LIGHT_AND_RIGHT_VERTICAL_HEAVY = '╊'
const BOX_DRAWINGS_HEAVY_VERTICAL_AND_HORIZONTAL = '╋'
const BOX_DRAWINGS_LIGHT_DOUBLE_DASH_HORIZONTAL = '╌'
const BOX_DRAWINGS_HEAVY_DOUBLE_DASH_HORIZONTAL = '╍'
const BOX_DRAWINGS_LIGHT_DOUBLE_DASH_VERTICAL = '╎'
const BOX_DRAWINGS_HEAVY_DOUBLE_DASH_VERTICAL = '╏'
const BOX_DRAWINGS_DOUBLE_HORIZONTAL = '═'
const BOX_DRAWINGS_DOUBLE_VERTICAL = '║'
const BOX_DRAWINGS_DOWN_SINGLE_AND_RIGHT_DOUBLE = '╒'
const BOX_DRAWINGS_DOWN_DOUBLE_AND_RIGHT_SINGLE = '╓'
const BOX_DRAWINGS_DOUBLE_DOWN_AND_RIGHT = '╔'
const BOX_DRAWINGS_DOWN_SINGLE_AND_LEFT_DOUBLE = '╕'
const BOX_DRAWINGS_DOWN_DOUBLE_AND_LEFT_SINGLE = '╖'
const BOX_DRAWINGS_DOUBLE_DOWN_AND_LEFT = '╗'
const BOX_DRAWINGS_UP_SINGLE_AND_RIGHT_DOUBLE = '╘'
const BOX_DRAWINGS_UP_DOUBLE_AND_RIGHT_SINGLE = '╙'
const BOX_DRAWINGS_DOUBLE_UP_AND_RIGHT = '╚'
const BOX_DRAWINGS_UP_SINGLE_AND_LEFT_DOUBLE = '╛'
const BOX_DRAWINGS_UP_DOUBLE_AND_LEFT_SINGLE = '╜'
const BOX_DRAWINGS_DOUBLE_UP_AND_LEFT = '╝'
const BOX_DRAWINGS_VERTICAL_SINGLE_AND_RIGHT_DOUBLE = '╞'
const BOX_DRAWINGS_VERTICAL_DOUBLE_AND_RIGHT_SINGLE = '╟'
const BOX_DRAWINGS_DOUBLE_VERTICAL_AND_RIGHT = '╠'
const BOX_DRAWINGS_VERTICAL_SINGLE_AND_LEFT_DOUBLE = '╡'
const BOX_DRAWINGS_VERTICAL_DOUBLE_AND_LEFT_SINGLE = '╢'
const BOX_DRAWINGS_DOUBLE_VERTICAL_AND_LEFT = '╣'
const BOX_DRAWINGS_DOWN_SINGLE_AND_HORIZONTAL_DOUBLE = '╤'
const BOX_DRAWINGS_DOWN_DOUBLE_AND_HORIZONTAL_SINGLE = '╥'
const BOX_DRAWINGS_DOUBLE_DOWN_AND_HORIZONTAL = '╦'
const BOX_DRAWINGS_UP_SINGLE_AND_HORIZONTAL_DOUBLE = '╧'
const BOX_DRAWINGS_UP_DOUBLE_AND_HORIZONTAL_SINGLE = '╨'
const BOX_DRAWINGS_DOUBLE_UP_AND_HORIZONTAL = '╩'
const BOX_DRAWINGS_VERTICAL_SINGLE_AND_HORIZONTAL_DOUBLE = '╪'
const BOX_DRAWINGS_VERTICAL_DOUBLE_AND_HORIZONTAL_SINGLE = '╫'
const BOX_DRAWINGS_DOUBLE_VERTICAL_AND_HORIZONTAL = '╬'
const BOX_DRAWINGS_LIGHT_ARC_DOWN_AND_RIGHT = '╭'
const BOX_DRAWINGS_LIGHT_ARC_DOWN_AND_LEFT = '╮'
const BOX_DRAWINGS_LIGHT_ARC_UP_AND_LEFT = '╯'
const BOX_DRAWINGS_LIGHT_ARC_UP_AND_RIGHT = '╰'
const BOX_DRAWINGS_LIGHT_DIAGONAL_UPPER_RIGHT_TO_LOWER_LEFT = '╱'
const BOX_DRAWINGS_LIGHT_DIAGONAL_UPPER_LEFT_TO_LOWER_RIGHT = '╲'
const BOX_DRAWINGS_LIGHT_DIAGONAL_CROSS = '╳'
const BOX_DRAWINGS_LIGHT_LEFT = '╴'
const BOX_DRAWINGS_LIGHT_UP = '╵'
const BOX_DRAWINGS_LIGHT_RIGHT = '╶'
const BOX_DRAWINGS_LIGHT_DOWN = '╷'
const BOX_DRAWINGS_HEAVY_LEFT = '╸'
const BOX_DRAWINGS_HEAVY_UP = '╹'
const BOX_DRAWINGS_HEAVY_RIGHT = '╺'
const BOX_DRAWINGS_HEAVY_DOWN = '╻'
const BOX_DRAWINGS_LIGHT_LEFT_AND_HEAVY_RIGHT = '╼'
const BOX_DRAWINGS_LIGHT_UP_AND_HEAVY_DOWN = '╽'
const BOX_DRAWINGS_HEAVY_LEFT_AND_LIGHT_RIGHT = '╾'
const BOX_DRAWINGS_HEAVY_UP_AND_LIGHT_DOWN = '╿'
const UPPER_HALF_BLOCK = '▀'
const LOWER_ONE_EIGHTH_BLOCK = '▁'
const LOWER_ONE_QUARTER_BLOCK = '▂'
const LOWER_THREE_EIGHTHS_BLOCK = '▃'
const LOWER_HALF_BLOCK = '▄'
const LOWER_FIVE_EIGHTHS_BLOCK = '▅'
const LOWER_THREE_QUARTERS_BLOCK = '▆'
const LOWER_SEVEN_EIGHTHS_BLOCK = '▇'
const FULL_BLOCK = '█'
const LEFT_SEVEN_EIGHTHS_BLOCK = '▉'
const LEFT_THREE_QUARTERS_BLOCK = '▊'
const LEFT_FIVE_EIGHTHS_BLOCK = '▋'
const LEFT_HALF_BLOCK = '▌'
const LEFT_THREE_EIGHTHS_BLOCK = '▍'
const LEFT_ONE_QUARTER_BLOCK = '▎'
const LEFT_ONE_EIGHTH_BLOCK = '▏'
const RIGHT_HALF_BLOCK = '▐'
const LIGHT_SHADE = '░'
const MEDIUM_SHADE = '▒'
const DARK_SHADE = '▓'
const UPPER_ONE_EIGHTH_BLOCK = '▔'
const RIGHT_ONE_EIGHTH_BLOCK = '▕'
const QUADRANT_LOWER_LEFT = '▖'
const QUADRANT_LOWER_RIGHT = '▗'
const QUADRANT_UPPER_LEFT = '▘'
const QUADRANT_UPPER_LEFT_AND_LOWER_LEFT_AND_LOWER_RIGHT = '▙'
const QUADRANT_UPPER_LEFT_AND_LOWER_RIGHT = '▚'
const QUADRANT_UPPER_LEFT_AND_UPPER_RIGHT_AND_LOWER_LEFT = '▛'
const QUADRANT_UPPER_LEFT_AND_UPPER_RIGHT_AND_LOWER_RIGHT = '▜'
const QUADRANT_UPPER_RIGHT = '▝'
const QUADRANT_UPPER_RIGHT_AND_LOWER_LEFT = '▞'
const QUADRANT_UPPER_RIGHT_AND_LOWER_LEFT_AND_LOWER_RIGHT = '▟'


text:
struct TupleVector{T, V <: Tuple} <: AbstractVector{T}
    vectors::V
    function TupleVector{T, V}(v) where {T, V <: Tuple}
        new{T, V}(v)
    end
end

Base.size(tv::TupleVector{T, V}) where {T, V} =
    min([length(vector) for vector in tv.vectors]...)

@generated function TupleVector{T}(n::I) where {T, I}
    types = fieldtypes(T)
    V = Tuple{(Vector{type} for type in types)...}
    args = [:(Vector{$type}(undef, n)) for type in types]
    :(TupleVector{T, $V}(($(args...),)))
end

@generated function Base.getindex(tv::TupleVector{T}, i::I) where {T, I}
    types = fieldtypes(T)
    members = [:(tv.vectors[$index][i])
               for (index, type) in enumerate(types)]
    Expr(:block, Expr(:meta, :inline), Expr(:new, T, members...))
end

@generated function Base.setindex!(tv::TupleVector{T}, t::T, i::I) where {T, I}
    types = fieldtypes(T)
    expressions = [:(tv.vectors[$index][i] = getfield(t, $index))
                   for (index, type) in enumerate(types)]
    :($(Expr(:block, Expr(:meta, :inline), expressions...)), t)
end

@noinline function resizeend!(tv::TupleVector{T}, n::I) where {T, I}
    if n > length(tv)
        for (i, _) in enumerate(tv.vectors)
            Base._growend!(tv.vectors[i], n - length(tv))
        end
    elseif n < length(tv)
        for (i, _) in enumerate(tv.vectors)
            Base._deleteend!(tv.vectors[i], length(tv) - n)
        end
    end
end
@noinline function resizebegin!(tv::TupleVector{T}, n::I) where {T, I}
    if n > length(tv)
        for (i, _) in enumerate(tv.vectors)
            Base._growbeg!(tv.vectors[i], n - length(tv))
        end
    elseif n < length(tv)
        for (i, _) in enumerate(tv.vectors)
            Base._deletebeg!(tv.vectors[i], length(tv) - n)
        end
    end
end


text:
@testset "Convergence failure" begin

    using DataFrames
    using DiscreteChoiceModels
    using Test
    using StatsBase

    data = DataFrame(
        choice = [2, 2, 2, 2, 2, 2, 2],
        pred =   [0, 0, 1, 1, 0, 0, 0],  # perfectly predicts outcome
        otherx = [0, 1, 0, 1, 1, 1, 0]
    )

    @test_skip 1
    # Despite my best efforts I can't get this to fail to converge

    #=
    @test_throws ConvergenceException multinomial_logit(
        # asc2 ML estimate is ∞
        @utility(begin
            1 ~ 0
            2 ~ :asc2
        end),
        data.choice,
        data
    )
    =#
end


text:
using Base.Test
using Attributes

g = 5;
b = BinaryAttributes(g)
nodes = 3;
attr = get_attributes(b, nodes)

@test size(attr) == (nodes, g) #test rozmiaru
@test all(abs.(attr) .== 1) #test wartosci powinny byc +-1

attr = ones(nodes, g)
attr[1, :] = -1
attr[2, 1:2] = -1
#attr is: [-1, -1, -1, -1, -1;
#          -1, -1, 1, 1, 1;
#           1, 1, 1, 1, 1]
#This should give following weights:
# [0.0  -0.2  -1.0
#  0.0   0.0   0.2
#  0.0   0.0   0.0]
exp_weights = zeros(3, 3)
exp_weights[1, 2] = -0.2;
exp_weights[1, 3] = -1.0;
exp_weights[2, 3] = 0.2;

weights = get_attribute_layer_weights(b, attr)

@test size(weights) == (nodes, nodes)
@test weights == exp_weights
# @test false

max_stars_repo_name: pjgorski/PolarizationFramework

max_issues_repo_name: pjgorski/PolarizationFramework

max_forks_repo_name: pjgorski/PolarizationFramework


text:
using BarycentricLagrange
using BenchmarkTools

# Chebyshev nodes
n = 100
x = cos.((0:n)*pi/n)
w = bclag_interp_weights(x)

F = x -> sin(x)
f = F.(x)

xx = collect(-1:0.01:1)
ref = F.(xx)

# Benchmarks
println("Direct eval")
@btime bclag_interp_eval(x, f, xx)
println("Compute weights")
@btime bclag_interp_weights(x)
println("Eval with weights")
@btime bclag_interp_eval(x, f, xx, w)

println("Direct matrix")
@btime bclag_interp_matrix(x, xx)
println("Matrix with weights")
@btime bclag_interp_matrix(x, xx, w)
M = bclag_interp_matrix(x, xx)
println("Matrix apply")
@btime M*f


text:
using Dates
using GRIB

"""
    $(typeof(GRIB_STEP_UNITS_TO_SECONDS))

Array used to convert the grib step units to seconds. As Julia is 1-indexed,
not 0 like Python, you should take care to correctly access the array,
typically  just +1 to the step units before using it as an index.

Taken from eccodes `stepUnits.table`.
"""
const GRIB_STEP_UNITS_TO_SECONDS = [
    60,
    3600,
    86400,
    missing,
    missing,
    missing,
    missing,
    missing,
    missing,
    missing,
    10800,
    21600,
    43200,
    1,
    900,
]

"""
    $(typeof(DEFAULT_EPOCH))

Default epoch used for `from_` and `to_` methods, set to `$DEFAULT_EPOCH.`
"""
const DEFAULT_EPOCH = DateTime(1970, 1, 1, 0, 0)

"""
Returns the integer seconds from epoch to the given date and time.

$(METHODLIST)
"""
function from_grib_date_time end

""
function from_grib_date_time(
    date::Int,
    time::Int;
    epoch::DateTime=DEFAULT_EPOCH
)::Int
    hour = time ÷ 100
    minute = time % 100
    year = date ÷ 10000
    month = date ÷ 100 % 100
    day = date % 100

    data_datetime = DateTime(year, month, day, hour, minute)

    return Dates.value(Dates.Second(data_datetime - epoch))
end

"""
Pulls out the date and time from given keys and passes them to
[`from_grib_date_time(::Int, ::Int)`](@ref from_grib_date_time(::Int, ::Int))
"""
function from_grib_date_time(
    message::GRIB.Message;
    date_key="dataDate",
    time_key="dataTime",
    epoch::DateTime=DEFAULT_EPOCH
)::Union{Int,Missing}
    if !haskey(message, date_key) || !haskey(message, time_key)
        return missing
    end

    date = message[date_key]
    time = message[time_key]

    return from_grib_date_time(date, time, epoch=epoch)
end

function to_grib_date_time(args...; kwargs...)
    throw(ErrorException("Unimplemented"))
end

"""
Returns the `step_key` value in hours.

Uses [`GRIB_STEP_UNITS_TO_SECONDS`](@ref GRIB_STEP_UNITS_TO_SECONDS) to convert
the step values to seconds, then divides by `3600.0` to get hours.
"""
function from_grib_step(
    message::GRIB.Message,
    step_key::String="endStep",
    step_unit_key::String="stepUnits"
)::Float64
    #  +1 as Julia is 1-indexed not 0
    to_seconds = GRIB_STEP_UNITS_TO_SECONDS[message[step_unit_key] + 1]
    return message[step_key] * to_seconds / 3600.0
end

function to_grib_step(args...; kwargs...)
    throw(ErrorException("Unimplemented"))
end

"""
Returns the integer seconds from the epoch to the verifying month value in the
GRIB message.
"""
function from_grib_month(
    message::GRIB.Message,
    verifying_month_key::String="verifyingMonth",
    epoch::DateTime=DEFAULT_EPOCH
)::Union{Int,Missing}
    if !haskey(message, verifying_month_key)
        return missing
    end

    date = message[verifying_month_key]
    year = date ÷ 100
    month = date % 100
    data_datetime = DateTime(year, month)

    return Dates.value(Dates.Second(data_datetime - epoch))
end

"""
Returns a pair of `(dims, data)` based on the type of input
"""
function build_valid_time end

"""
```jldoctest
julia> CfGRIB.build_valid_time(10, 10)
((), 36010)
```
"""
function build_valid_time(
    time::Int, step::Int
)::Tuple{Tuple{},Int64}
    step_s = step * 3600

    data = time + step_s
    dims = ()

    return dims, data
end

"""
```jldoctest
julia> CfGRIB.build_valid_time([10], 10)
(("time",), [36010])
```
"""
function build_valid_time(
    time::Array{Int,1}, step::Int
)::Tuple{Tuple{String},Array{Int64,1}}
    step_s = step * 3600

    data = time .+ step_s
    dims = ("time",)

    return dims, data
end

"""
```jldoctest
julia> CfGRIB.build_valid_time(1, [10])
(("step",), [36001])
```
"""
function build_valid_time(
    time::Int, step::Array{Int,1}
)::Tuple{Tuple{String},Array{Int64,1}}
    step_s = step * 3600

    data = time .+ step_s
    dims = ("step",)

    return dims, data
end

"""
```jldoctest
julia> CfGRIB.build_valid_time([10, 10], [10, 10])
(("time", "step"), [36010 36010; 36010 36010])
```

```jldoctest
julia> CfGRIB.build_valid_time([10], [10])
((), 36010)
```
"""
function build_valid_time(
    time::Array{Int,1}, step::Array{Int,1}
)::Union{
    Tuple{Tuple{},Int64},
    Tuple{Tuple{String,String},Array{Int64,2}}
}
    step_s = step * 3600

    if length(time) == 1 && length(step) == 1
        return build_valid_time(time[1], step[1])
    end

    data = time' .+ step_s
    dims = ("time", "step")
    return dims, data
end

"""
Dictionary which maps a key to a conversion method. The first function is the
'to' conversion, the second is 'from'.

TODO: Actually applying the `from_grib_step` function results in different values
to cfgrib.py, so `step -> (from_grib_step, to_grib_step)` is currently disabled.

Currently converts:

```
    "time" => (from_grib_date_time, to_grib_date_time)

    "valid_time" => (
        message -> from_grib_date_time(message, date_key="validityDate", time_key="validityTime"),
        message -> to_grib_date_time(message, date_key="validityDate", time_key="validityTime"),
    )

    "verifying_time" => (from_grib_month, m -> throw(ErrorException("Unimplemented")))

    "indexing_time" => (
        message -> from_grib_date_time(message, date_key="indexingDate", time_key="indexingTime"),
        message -> to_grib_date_time(message, date_key="indexingDate", time_key="indexingTime"),
    )
```

# Example

A GRIB message containing `20160501` as the date key and `0` as the time key
would end up calling:

```jldoctest
julia> CfGRIB.COMPUTED_KEYS["time"](20160501, 0)
1462060800
```
"""
COMPUTED_KEYS = Dict(
    "time" => (from_grib_date_time, to_grib_date_time),
    #  TODO: Actually applying the from_grib_step function results in different
    #  values to cfgrib.py...?
    # "step" => (from_grib_step, to_grib_step),
    "valid_time" => (message -> from_grib_date_time(message, date_key="validityDate", time_key="validityTime"),
        message -> to_grib_date_time(message, date_key="validityDate", time_key="validityTime"),),
    "verifying_time" => (from_grib_month, m -> throw(ErrorException("Unimplemented"))),
    "indexing_time" => (message -> from_grib_date_time(message, date_key="indexingDate", time_key="indexingTime"),
        message -> to_grib_date_time(message, date_key="indexingDate", time_key="indexingTime"),),
)

"""
Reads a specific key from a GRIB.jl message. Attempts to convert the raw value
associated with that key using the [`COMPUTED_KEYS`](@ref COMPUTED_KEYS) mapping
to `from_grib_*` functions.
"""
function read_message(message::GRIB.Message, key::String)
    value = missing

    if key in keys(COMPUTED_KEYS)
        value = COMPUTED_KEYS[key][1](message)
    end

    if ismissing(value)
        value = haskey(message, key) ? message[key] : missing
    end

    value = value isa Array ? Tuple(value) : value

    return value
end

#  TODO: implement other conversion methods, but some seem unused, should these
#  be implemented as well:
#   - to_grib_date_time
#   - from_grib_step
#   - to_grib_step


text:
if basename(pwd()) == "aoc"
    cd("2016/14")
end

max_stars_repo_name: stellartux/advent-of-code

max_issues_repo_name: stellartux/advent-of-code

max_forks_repo_name: stellartux/advent-of-code


text:
deltaH = -93.2
deltaS = -.1983
T = 300

gibbsFreeEnergy = deltaH - T*deltaS
gibbsFreeEnergy2 = -93.2 - (300*-.1983)

print("Gibbs Free Energy: ")
print("$gibbsFreeEnergy2 kJ/mol")

max_stars_repo_name: j-hayes/intro-to-computing-sci-eng

max_issues_repo_name: j-hayes/intro-to-computing-sci-eng

max_forks_repo_name: j-hayes/intro-to-computing-sci-eng


text:
struct BinaryTreeNode
    left
    right
end

Base.show(io::IO, blk::BinaryTreeNode) = show(io, "plain/text", blk)
function Base.show(io::IO, ::MIME"plain/text", blk::BinaryTreeNode)
    print(io, "Contraction Tree: ")
    print_tree(io, blk)
end

function print_tree(io::IO, blk::BinaryTreeNode, print_level=0)
    print(io, "(")
    print_tree(io, blk.left, print_level+1)
    print(io, " ↔ ")
    print_tree(io, blk.right, print_level+1)
    print(io, ")")
end

function print_tree(io::IO, blk, print_level=0)
    print(io, blk)
end

Base.getindex(t::BinaryTreeNode, i::Int) = i==1 ? t.left : (i==2 ? t.right : throw(BoundsError(t, i)))
Base.iterate(t::BinaryTreeNode, args...) = iterate((t.left, t.right), args...)


text:
function buildgti(t)
    n=length(t)
    tstart =Float64[]
    tstop = Float64[]
    push!(tstart,t[1] - 1.1886703073978424)
    for i=1:n-1
        if (t[i+1] - t[i]) > 2.4
            push!(tstop,t[i] + 1.1886703073978424)
            println(t[i+1] - t[i])
            push!(tstart,t[i+1] - 1.1886703073978424)
        end
    end
    return tstart[1:end-1],  tstop
end

max_stars_repo_name: gulabd/SXTMerger

max_issues_repo_name: gulabd/SXTMerger


text:
import JSON
import MbedTLS
using Dates, Base64
using HTTP, HTTP.Messages

const API_ROOT = "https://www.googleapis.com"
const SCOPE_ROOT = "$API_ROOT/auth"
const AUD_ROOT = "$API_ROOT/oauth2/v4/token"

"""
    isurl(path)

Return true if `path` is a URL and false a path fragment.
"""
isurl(path::AbstractString) = occursin(r"^https?://", path)

"""
    SHA256withRSA(message, key)

Sign message using private key with RSASSA-PKCS1-V1_5-SIGN algorithm.
"""
SHA256withRSA(message, key::MbedTLS.PKContext) = MbedTLS.sign(key, MbedTLS.MD_SHA256,
    MbedTLS.digest(MbedTLS.MD_SHA256, message), MbedTLS.MersenneTwister(0)
)

"""
Session(...)

OAuth 2.0 session for Google using provided credential.

Caches authorisation tokens up to expiry.

```julia
sess = Session(Credential(expanduser("~/auth.json")), ["devstorage.full_control"])
```
"""
mutable struct Session{T <: Credential}
    credential::T
    scopes::Vector{String}
    authorization::Dict{Symbol, Any}
    expiry::DateTime
    """
        Session(credential, scopes)

    Set up session with credential and OAuth scopes.
    """
    function Session(credential::T, scopes::AbstractVector{<: AbstractString}) where {T <: Credential}
        scopes = [isurl(scope) ? scope : "$SCOPE_ROOT/$scope" for scope in scopes]
        new{T}(credential, scopes, Dict{String, String}(), DateTime(1))
    end
end

function Session(credential::Union{AbstractString, Nothing}, scopes::AbstractVector{<: AbstractString}=String[])
    if credential === nothing
        credential = ""
    end
    if isfile(credential)
        credential = Credential(credential)
    end
    Session(credential, scopes)
end

function Session(scopes::AbstractVector{<: AbstractString})
    Session(get(ENV, "GOOGLE_APPLICATION_CREDENTIALS", ""), scopes)
end

function Base.print(io::IO, x::Session)
    println(io, "scopes: $(x.scopes)")
    println(io, "authorization: $(x.authorization)")
    println(io, "expiry: $(x.expiry)")
end
Base.show(io::IO, x::Session) = print(io, x)

"""
    unixseconds(x)

Convert date-time into unix seconds.
"""
unixseconds(x::DateTime) = trunc(Int, datetime2unix(x))

"""
    JWTHeader

JSON Web Token header.
"""
struct JWTHeader
    algorithm::String
end
function Base.string(x::JWTHeader)
    base64encode(JSON.json(Dict(:alg => x.algorithm, :typ => "JWT")))
end
Base.print(io::IO, x::JWTHeader) = print(io, string(x))

"""
    JWTClaimSet

JSON Web Token claim-set.
"""
struct JWTClaimSet
    issuer::String
    scopes::Vector{String}
    assertion::DateTime
    expiry::DateTime
    function JWTClaimSet(issuer::AbstractString, scopes::AbstractVector{<: AbstractString}, assertion::DateTime=now(UTC), expiry::DateTime=now(UTC) + Hour(1))
        new(issuer, scopes, assertion, expiry)
    end
end
function Base.string(x::JWTClaimSet)
    base64encode(JSON.json(Dict(
        :iss => x.issuer,
        :scope => join(x.scopes, " "),
        :aud => AUD_ROOT,
        :iat => unixseconds(x.assertion),
        :exp => unixseconds(x.expiry)
    )))
end
Base.print(io::IO, x::JWTClaimSet) = print(io, string(x))

"""
    JWS(credential, claimset)

Construct the Base64-encoded JSON Web Signature based on the JWT header, claimset
and signed using the private key provided in the Google JSON service-account key.
"""
function JWS(credential::Credential, claimset::JWTClaimSet, header::JWTHeader=JWTHeader("RS256"))
    payload = "$header.$claimset"
    signature = base64encode(SHA256withRSA(payload, credential.private_key))
    "$payload.$signature"
end

function token(credential::Credential, scopes::AbstractVector{<: AbstractString}; kws...)
    # construct claim-set from service account email and requested scopes
    claimset = JWTClaimSet(credential.client_email, scopes)
    data = HTTP.URIs.escapeuri(Dict{Symbol, String}(
        :grant_type => "urn:ietf:params:oauth:grant-type:jwt-bearer",
        :assertion => JWS(credential, claimset)
    ))
    headers = Dict{String, String}("Content-Type" => "application/x-www-form-urlencoded")
    res = HTTP.post("$AUD_ROOT", headers, data; kws...)
    if res.status != 200
        throw("Unable to obtain authorization: $(read(res, String))")
    end
    authorization = JSON.parse(payload(res, String); dicttype=Dict{Symbol, Any})
    authorization, claimset.assertion
end

"""
    authorize(session; cache=true)

Get OAuth 2.0 authorisation token from Google.

If `cache` set to `true`, get a new token only if the existing token has expired.
"""
function authorize(session::Session; cache::Bool=true, kws...)
    # don't get a new token if a non-expired one exists
    if cache && (session.expiry >= now(UTC)) && !isempty(session.authorization)
        return session.authorization
    end

    authorization, assertion = try token(session.credential, session.scopes; kws...) catch e
        session.expiry = DateTime(1)
        empty!(session.authorization)
        rethrow(e)
    end

    # cache authorization if required
    if cache
        session.expiry = assertion + Second(authorization[:expires_in] - 30)
        session.authorization = authorization
    else
        session.expiry = DateTime(1)
        empty!(session.authorization)
    end
    authorization
end

max_stars_repo_name: ProfLeao/mais

max_issues_repo_name: ProfLeao/mais

max_forks_repo_name: ProfLeao/mais


text:
function linsolve_lu(A::Matrix{T}, b::Vector{T}) where {T}

    (n_row, n_col) = size(A)
    if n_row != n_col
        error("This function only works with the same dimensions")
    end
    x = zeros(eltype(b), n_row)

    (L, U, p) = lu_decomp_partialpivot(A)
    b_ = b[p]
    y_ = triangular_solve_lower_colwise(L, b_)
    x = triangular_solve_upper_colwise(U, y_)
    # y_ = triangular_solve_lower_rowwise(L, b_)
    # x = triangular_solve_upper_rowwise(U, y_)

    return x
end

function lu_decomp_partialpivot(A::Matrix{T}) where {T}

    (n_row, n_col) = size(A)

    M = copy(A)
    p = collect(1:n_row)

    for jj = 1:(n_col-1)

        @inbounds (max_elem, maxidx) = util.maximum( abs, @view( M[jj:n_row, jj] ) )
        if maxidx > 1
            maxidx = maxidx + jj - 1
            # swap the row
            @inbounds p[jj], p[maxidx] = p[maxidx], p[jj]
            @simd for cc = 1:n_col
                @inbounds M[jj, cc], M[maxidx, cc] = M[maxidx, cc], M[jj, cc]
            end

        end

        # calculate the L
        @fastmath @simd for rr = (jj+1):n_row
            M[rr, jj] /= M[jj, jj]
        end

        for cc = (jj+1):n_col
            @inbounds m_ = M[jj, cc]
            @simd for rr = (jj+1):n_row
                @fastmath @inbounds M[rr, cc] -= M[rr, jj] .* m_
            end
        end

    end

    L = UnitLowerTriangular(M);
    U = UpperTriangular(M);
    return (L, U, p)
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
using OrdinaryDiffEq, StochasticDiffEq, DelayDiffEq, Flux,
      Zygote, DiffEqFlux, Test

x = Float32[2.; 0.]
xs = Float32.(hcat([0.; 0.], [1.; 0.], [2.; 0.]))
tspan = (0.0f0, 1.0f0)
fastdudt = FastChain(FastDense(4, 50, tanh), FastDense(50, 4))
fastdudt2 = FastChain(FastDense(4, 50, tanh), FastDense(50, 4))
fastdudt22 = FastChain(FastDense(4, 50, tanh), FastDense(50, 16), (x, p) -> reshape(x, 4, 4))
fastddudt = FastChain(FastDense(12, 50, tanh), FastDense(50, 4))

# Augmented Neural ODE
anode = AugmentedNDELayer(
    NeuralODE(fastdudt, tspan, Tsit5(), save_everystep=false, save_start=false), 2
)
anode(x)

grads = Zygote.gradient(() -> sum(anode(x)), Flux.params(x, anode))
@test ! iszero(grads[x])
@test ! iszero(grads[anode.p])

# Augmented Neural DSDE
andsde = AugmentedNDELayer(
    NeuralDSDE(fastdudt, fastdudt2, (0.0f0, 0.1f0), SOSRI(), saveat=0.0:0.01:0.1), 2
)
andsde(x)

grads = Zygote.gradient(() -> sum(andsde(x)), Flux.params(x, andsde))
@test ! iszero(grads[x])
@test ! iszero(grads[andsde.p])

# Augmented Neural SDE
asode = AugmentedNDELayer(
    NeuralSDE(fastdudt, fastdudt22,(0.0f0, 0.1f0), 4, LambaEM(), saveat=0.0:0.01:0.1), 2
)
asode(x)

grads = Zygote.gradient(() -> sum(asode(x)), Flux.params(x, asode))
@test ! iszero(grads[x])
@test ! iszero(grads[asode.p])

# Augmented Neural CDDE
adode = AugmentedNDELayer(
    NeuralCDDE(fastddudt, (0.0f0, 2.0f0), (p, t) -> zeros(Float32, 4), (1f-1, 2f-1),
               MethodOfSteps(Tsit5()), saveat=0.0:0.1:2.0), 2
)
adode(x)

grads = Zygote.gradient(() -> sum(adode(x)), Flux.params(x, adode))
@test ! iszero(grads[x])
@test ! iszero(grads[adode.p])


text:
__precompile__()

module StaticPolynomials

    import MultivariatePolynomials
    const MP = MultivariatePolynomials
    import StaticArrays: SVector, MVector, SMatrix
    import LinearAlgebra
    if VERSION < v"1.0-"
        import LinearAlgebra: gradient
    end

    include("helpers.jl")
    include("sexponents.jl")
    include("polynomial.jl")
    include("show.jl")

    include("codegen_helpers.jl")
    include("evalpoly.jl")
    include("evaluate_codegen.jl")
    include("gradient_codegen.jl")
    include("evaluation.jl")

    include("polynomial_system.jl")
end # module


text:
module CatElements
export ThElements, AbstractElements, Elements, elements

using DataStructures: OrderedDict

using ..CSets, ..FinSets
using ...Present, ...Theories
using ...Theories: Category, ob, hom, dom_nums, codom_nums

@present ThElements(FreeSchema) begin
  (El, Arr, Ob, Hom)::Ob
  Name::AttrType
  src::Hom(Arr, El)
  tgt::Hom(Arr, El)
  dom::Hom(Hom, Ob)
  cod::Hom(Hom, Ob)
  πₑ::Hom(El, Ob)
  πₐ::Hom(Arr, Hom)
  nameo::Attr(Ob, Name)
  nameh::Attr(Hom, Name)

  src⋅πₑ == πₐ⋅dom
  tgt⋅πₑ == πₐ⋅cod
end

@abstract_acset_type AbstractElements
@acset_type Elements(ThElements, index=[:src, :tgt, :πₑ, :πₐ]) <: AbstractElements

"""    elements(X::AbstractACSet)

construct the category of elements from an ACSet. This only correctly handles the CSet part.
This transformation converts an instance of C into a Graph homomorphism. The codomain of the
homomorphism is a graph shaped like the schema. This is one half of the isomorphism between
databases and knowledge graphs.
"""
function elements(X::StructACSet{S}) where S
  Y = Elements{Symbol}()

  obs = ob(S)
  add_parts!(Y, :Ob, length(obs), nameo=obs)
  els = map(enumerate(obs)) do (i,c)
    add_parts!(Y, :El, nparts(X, c), πₑ = i)
  end

  add_parts!(Y, :Hom, length(hom(S)), dom=dom_nums(S), cod=codom_nums(S), nameh=hom(S))
  map(enumerate(zip(hom(S), dom_nums(S), codom_nums(S)))) do (i, (f, ci, di))
    c, d = obs[ci], obs[di]
    nc = nparts(X, c)
    add_parts!(Y, :Arr, nparts(X, c), src=els[ci], tgt=view(els[di], X[f]), πₐ=i)
  end
  return Y
end

"""    presentation(X::AbstractElements)

convert a category of elements into a new schema. This is useful for generating large schemas
that are defined as the category of elements of a specified C-Set. For example, the schema for
Bipartite Graphs is the category of elements of the graph with 2 vertices and 1 edge. The 2 vertices
will get mapped to elements `V_1, V_2` and the one edge will be `E_1` and the source and target maps
will be mapped to `src_E_1, tgt_E_1`.
"""
function presentation(X::AbstractElements)
  P = Presentation(FreeSchema)
  obs = OrderedDict{Tuple{Symbol, Int}, Any}()
  for a in 1:nparts(X, :Ob)
    class = X[a, :nameo]
    for (i,j) in enumerate(incident(X, a, :πₑ))
      ob = Ob(FreeSchema, Symbol("$(class)_$i"))
      obs[(class, j)] = ob
    end
  end
  add_generators!(P, values(obs))
  homs = OrderedDict{Tuple{Symbol, Int}, Any}()
  for f in 1:nparts(X, :Hom)
    class = X[f, :nameh]
    for (i,j) in enumerate(incident(X, f, :πₐ))
      domkey = (X[j, [:src, :πₑ, :nameo]], X[j, :src])
      dom = obs[domkey]
      codkey = (X[j, [:tgt, :πₑ, :nameo]], X[j, :tgt])
      codom = obs[codkey]
      homs[(class, j)] = Hom(Symbol("$(class)_$dom"), dom, codom)
    end
  end
  add_generators!(P, values(homs))
  return P, obs, homs
end
end


text:
module MountainCarLearnPolicy

using Curiosity, Statistics, Plots, FileIO, BSON
using Random, StatsBase, RollingFunctions
using GVFHordes, ProgressMeter

const MCU = Curiosity.MountainCarUtils

# Cumulant
struct Reward end

(c::Reward)(;kwargs...) = kwargs[:r]

struct ConstantDiscount{F}
    val::F
end

(d::ConstantDiscount)(;kwargs...) = d.val

# Environment
function construct_env()
    normalized = true
    env = MountainCar(0.0, 0.0, normalized)
end

# Agent
function construct_agent(numtilings, numtiles, lu_str, α, λ, ϵ, γ, policy_name)
    obs_size = 2
    fc = Curiosity.SparseTileCoder(numtilings, numtiles, obs_size)
    feature_size = size(fc)

    lu = if lu_str == "ESARSA"
        ESARSA(lambda=λ, opt=Curiosity.Descent(α), trace=ReplacingTraces())
    elseif lu_str == "SARSA"
        SARSA(lambda=λ, opt=Curiosity.Descent(α))
    elseif lu_str == "TB"
        TB(lambda=λ, opt=Curiosity.Descent(α))
    else
        throw(ArgumentError("$(lu_str) Not a valid behaviour learning update"))
    end
#     (update, num_features, num_actions, num_demons, w_init)
    w_init = 0 / numtilings
    learner = LinearQLearner(lu, feature_size, 3, 1,w_init)
    exploration = EpsilonGreedy(ϵ)
    # cumulant = Reward()
    # discount = ConstantDiscount(γ)
    # discount = GVFParamFuncs.ConstantDiscount(γ)

    cumulant,discount = if policy_name == "Goal"
        cumulant = MCU.step_to_goal_cumulant()
        discount = MCU.goal_pseudoterm(γ)
        cumulant,discount
    elseif policy_name == "Wall"
        cumulant = MCU.step_to_wall_cumulant()
        discount = MCU.wall_pseudoterm(γ)
        cumulant,discount
    elseif policy_name == "Task"
        cumulant = MCU.leave_cumulant()
        discount = MCU.goal_pseudoterm(γ)
        cumulant,discount
    else
        throw(ArgumentError("What Policy do you want to learn?? "))
    end

    b_gvf = make_behaviour_gvf(learner, discount, cumulant, fc, exploration)
    b_demons = Horde([b_gvf])

    Curiosity.PolicyLearner(learner,
                            fc,
                            exploration,
                            discount,
                            cumulant,
                            zeros(2),
                            0,
                            b_demons)
end

function make_behaviour_gvf(behaviour_learner, discount, cumulant, fc, exploration_strategy)
    function b_π(state_constructor, learner, exploration_strategy; kwargs...)
        s = state_constructor(kwargs[:state_t])
        preds = learner(s)
        return exploration_strategy(preds)[kwargs[:action_t]]
    end
    GVF_policy = GVFParamFuncs.FunctionalPolicy((;kwargs...) -> b_π(fc, behaviour_learner, exploration_strategy; kwargs...))
    BehaviourGVF = GVF(cumulant,discount, GVF_policy)
    return BehaviourGVF
end

function main_experiment(policy_name="tmp";progress=true, file_name = nothing)
    # Learn the policy

    if file_name isa Nothing
        file_name = policy_name
    end

    seed = 1029
    Random.seed!(seed)
    numtilings, numtiles = 8, 8
    lu_str = "ESARSA"
    α = 0.1/numtilings
    λ = 0.9
    ϵ = 0.1
    γ = 0.99

    info = Dict(
        "seed"=>seed,
        "numtilings"=>numtilings,
        "numtiles"=>numtiles,
        "lu"=>"lu_str",
        "α"=>α,
        "λ"=>λ,
        "ϵ"=>ϵ,
        "γ"=>γ,
        "rew"=>"Env"
    )

    env = construct_env()

    agent = construct_agent(numtilings, numtiles, lu_str, α, λ, ϵ, γ, policy_name)

    steps = Int[]
    ret = Float64[]
    max_num_steps = 500000
    eps = 0
    prg_bar = ProgressMeter.Progress(max_num_steps, "Step: ")

    while sum(steps) < max_num_steps
        is_terminal = false

        max_episode_steps = min(max_num_steps - sum(steps), 2000)
        s = start!(env)
        a = start!(agent, s)
        stp = 0
        a = 0
        tr = 0.0

        while !is_terminal && stp <= max_episode_steps
            s, r, is_terminal = MinimalRLCore.step!(env, a)
            a = MinimalRLCore.step!(agent, s, r, is_terminal)
            tr += r
            stp += 1

            if progress
                next!(prg_bar)
            end
        end
        push!(steps, stp)
        push!(ret, tr)

        eps += 1
    end
    Curiosity.save(agent, "./src/data/MC_learned_policies/$(policy_name).bson", info)

    steps,ret
end

end

max_stars_repo_name: MatthewMcLeod/curiosity

max_issues_repo_name: MatthewMcLeod/curiosity

max_forks_repo_name: MatthewMcLeod/curiosity


text:
function interpolate!(ex,env)
    if isa(ex,Symbol)
        if haskey(env, ex)
            return env[ex]
        end
    elseif isa(ex,Expr)
        for (i,x) in enumerate(ex.args)
            if isa(x,Symbol)
                if haskey(env, x)
                    ex.args[i] = env[x]
                end
            elseif isa(x,Expr)
                ex.args[i] = interpolate!(ex.args[i], env)
            end
        end
    end
    return ex
end

function index2coord(index::Integer, L::AbstractArray)
    D = length(L)
    coord = zeros(Int,D)
    for d in 1:D
        coord[d] = index%L[d]
        index ÷= L[d]
    end
    return coord
end

function coord2index(coord::AbstractArray, L::AbstractArray)
    D = length(L)
    index = 0
    offset = 0
    c = mod.(coord, L)
    for d in D:-1:1
        index *= L[d]
        index += c[d]
    end
    return index
end

function generatelattice(param)
    latticedict = get(param, "LatticeDict", stdlattices)
    lat = latticedict[param["Lattice"]]
    return generator(lat)(param)
end

function generatelattice_std(param)
    bravaisdict = get(param, "BravaisDict", stdbravais)
    unitcelldict = get(param, "UnitcellDict", stdunitcells)
    latticedict = get(param, "LatticeDict", stdlattices)
    lat = latticedict[param["Lattice"]]
    bra = bravaisdict[lat.bravais]
    ucell = unitcelldict[lat.unitcell]
    @assert lat.dimension == bra.dimension == ucell.dimension
    D = lat.dimension
    bparams = Dict(name => get(param,String(name),default) for (name, default) in bra.parameters)
    for (name, default) in lat.parameters
        bparams[name] = get(param,String(name),default)
    end

    bc = get(param, "Periodic Boudary Condition", lat.periodic)

    if isa(param["L"], Vector)
        L = param["L"]
    else
        L = [param["L"]]
    end
    while length(L) < D
        push!(L,L[end])
    end

    numcell = prod(L)

    latvec = interpolate!(bra.basis, bparams) |> eval

    usites = ucell.sites
    numsites_in_cell = length(usites)
    numsites = numcell * numsites_in_cell

    let x = Int[]
        for site in usites
            push!(x, site.id)
        end
        sort!(x)
        @assert x == collect(1:numsites_in_cell)
    end

    ubonds = ucell.bonds

    numbonds = 0
    numsitetypes = 0
    numbondtypes = 0
    nss = Int[]
    nbs = Int[]

    sites = Site[]
    bonds = Bond[]
    ib = 0

    use_index_as_sitetype = get(param, "Use Indecies as Site Types", false)
    use_index_as_bondtype = get(param, "Use Indecies as Bond Types", false)

    for icell in 0:(numcell-1)
        cellcoord = index2coord(icell, L)
        for site in usites
            id = numsites_in_cell * icell + site.id
            coord = latvec * (cellcoord .+ site.coord)
            s = Site(id, ifelse(use_index_as_sitetype, id, site.sitetype), Int[], Int[], coord, site.id, cellcoord)
            push!(sites,s)
        end
        for bond in ubonds
            for d in 1:D
                if !(bc[d] || bond.source.offset[d] == bond.source.offset[d])
                    if !(0 <= cellcoord[d] + (bond.target.offset[d] - bond.source.offset[d]) < L[d] )
                        continue
                    end
                end
            end
            source = numsites_in_cell * coord2index(cellcoord .+ bond.source.offset, L) + bond.source.id
            target = numsites_in_cell * coord2index(cellcoord .+ bond.target.offset, L) + bond.target.id
            dir = latvec * ((bond.target.offset .+ usites[bond.target.id].coord )
                            .- (bond.source.offset .+ usites[bond.source.id].coord ))
            ib += 1
            b = Bond(ib, ifelse(use_index_as_bondtype, ib, bond.bondtype), source, target, dir)
            push!(bonds, b)
        end
    end

    return Lattice(latvec, L, sites, bonds)
end

function generate_fully_connected_graph(param)
    N = param["N"]
    nb = div(N*(N-1),2)
    sites = Site[]
    bonds = Bond[]
    ib = 0
    for s in 1:N
        push!(sites, Site(s, 1, Int[], Int[], [0.0], 1, [0]))
        for t in (s+1):N
            ib += 1
            push!(bonds, Bond(ib, 1, s, t, [0.0]))
        end
    end
    latvec = zeros(1,1)
    return Lattice(latvec, [N], sites, bonds)
end


text:
macro _const(ex)
    ex = esc(ex)
    if VERSION < v"1.8.0-DEV.1148"
        return ex
    else
        return Expr(:const, ex)
    end
end

const var"@const" = var"@_const"

@def_singleton TRUE
@def_singleton FALSE
const BOOL = Union{typeof(TRUE),typeof(FALSE)}
const AnyBool = Union{BOOL,Bool}

boolean(b::Bool) = b
boolean(::typeof(TRUE)) = true
boolean(::typeof(FALSE)) = false

macro return_if_something(ex, before_return = nothing)
    quote
        ans = $(esc(ex))
        if ans === nothing
        else
            $(esc(before_return))
            return something(ans)
        end
    end
end

function duplicated_items_in_sorted_iterable(xs)
    dups = nothing
    y = iterate(xs)
    y === nothing && return dups
    prev, state = y
    while true
        y = iterate(xs, state)
        y === nothing && return dups
        curr, state = y
        if isless(prev, curr)
            prev = curr
        elseif isless(curr, prev)
            error("Input is not sorted")
        else
            dups = @something(dups, eltype(xs)[])
            push!(dups, curr)
        end
    end
end

Historic.@define Debug
using .Debug: @record

"""
    @yield_unsafe expression

Document that `expression` must not yield to the Julia scheduler.

This macro is used purely for documentation.  The `expression` is evaluated
as-is.
"""
macro yield_unsafe(ex)
    esc(ex)
end

function contextualized_buffer(io::IO)
    buffer = IOBuffer()
    ctx = IOContext(buffer, io)
    ctx = IOContext(ctx, :displaysize => displaysize(io))
    return (ctx, buffer)
end

""" `printing(f, io::IO)` is same as `f(io)` but tries to print it in "one shot" """
function printing(f, io::IO)
    (ctx, buffer) = contextualized_buffer(io)
    y = f(ctx)
    seekstart(buffer)
    write(io, take!(buffer))
    return y
end

struct LazyException <: Exception
    f::Any
    data::NamedTuple
end

function Base.showerror(io::IO, e::LazyException)
    printing(io) do io
        e.f(io, e.data)
    end
end

errorwith(f; data...) = throw(LazyException(f, values(data)))

donothing(_args...; _kwargs...) = nothing

function waitall(waitables; onerror = donothing, logerror = nothing)
    exception = nothing
    for w in waitables
        try
            Base.wait(w)
        catch err
            if logerror !== nothing
                @error "$logerror" exception = (err, catch_backtrace())
            end
            exception = @something(exception, CompositeException())
            push!(exception, err)
            onerror(w, err)
        end
    end
    exception === nothing || throw(exception)
end

function pause()
    ccall(:jl_cpu_pause, Cvoid, ())
    GC.safepoint()
end

const EMPTY_VECTOR = Union{}[]

struct DontCallMe end
const wait = DontCallMe()
const schedule = DontCallMe()
# const yield = DontCallMe()


text:
#Official site: https://julialang.org/

println("Hello world!")

#From Julia REPL
# julia> include("/path/to/Julia.jl")
#From terminal
# $ /path/to/julia Julia.jl


text:
module PlasmoSolverInterface

using Plasmo
using Libdl
using MPI
using Distributed

import Plasmo: solve, AbstractPlasmoSolver

export PipsSolver,pipsnlp_solve,dsp_solve

include("wrapped_solvers.jl")

include("plasmoPipsNlpInterface.jl")

include("DspCInterface.jl")

include("plasmoDspInterface.jl")

using .PlasmoPipsNlpInterface

using .PlasmoDspInterface

end # module


text:
# ---
# title: 1260. Shift 2D Grid
# id: problem1260
# author: Tian Jun
# date: 2020-11-14
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/shift-2d-grid/description/>
# hidden: true
# ---
#
# Given a 2D `grid` of size `m x n` and an integer `k`. You need to shift the
# `grid` `k` times.
#
# In one shift operation:
#
#   * Element at `grid[i][j]` moves to `grid[i][j + 1]`.
#   * Element at `grid[i][n - 1]` moves to `grid[i + 1][0]`.
#   * Element at `grid[m - 1][n - 1]` moves to `grid[0][0]`.
#
# Return the _2D grid_ after applying shift operation `k` times.
#
#
#
# **Example 1:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/05/e1.png)
#
#
#
#     Input: grid = [[1,2,3],[4,5,6],[7,8,9]], k = 1
#     Output: [[9,1,2],[3,4,5],[6,7,8]]
#
#
# **Example 2:**
#
# ![](https://assets.leetcode.com/uploads/2019/11/05/e2.png)
#
#
#
#     Input: grid = [[3,8,1,9],[19,7,2,5],[4,6,11,10],[12,0,21,13]], k = 4
#     Output: [[12,0,21,13],[3,8,1,9],[19,7,2,5],[4,6,11,10]]
#
#
# **Example 3:**
#
#
#
#     Input: grid = [[1,2,3],[4,5,6],[7,8,9]], k = 9
#     Output: [[1,2,3],[4,5,6],[7,8,9]]
#
#
#
#
# **Constraints:**
#
#   * `m == grid.length`
#   * `n == grid[i].length`
#   * `1 <= m <= 50`
#   * `1 <= n <= 50`
#   * `-1000 <= grid[i][j] <= 1000`
#   * `0 <= k <= 100`
#
#
## @lc code=start
using LeetCode

function shift_grid(grid::Matrix{Int}, k::Int)::Matrix{Int}
    A = grid'
    res = Matrix(undef, size(A)...)
    LI = LinearIndices(A)

    for i in LI
        i′ = mod1(i + k, length(LI))
        res[i′] = A[i]
    end
    return res'
end
## @lc code=end


text:
using DiffEqOptimalControl
using Base.Test