text:
# This file is a companion to `src/api.jl` --- it defines the raw ccall wrappers, while
# here small normalizations are made to make the calls more Julian.
# For instance, many property getters return values through pointer output arguments,
# so the methods here handle making the appropriate `Ref`s and return them (as tuples).

###
### HDF5 General library functions
###

function h5_get_libversion()
    majnum, minnum, relnum = Ref{Cuint}(), Ref{Cuint}(), Ref{Cuint}()
    h5_get_libversion(majnum, minnum, relnum)
    VersionNumber(majnum[], minnum[], relnum[])
end

function h5_is_library_threadsafe()
    is_ts = Ref{Cuint}()
    h5_is_library_threadsafe(is_ts)
    return is_ts[] > 0
end

###
### Attribute Interface
###

function h5a_get_name(attr_id)
    len = h5a_get_name(attr_id, 0, C_NULL)
    buf = StringVector(len)
    h5a_get_name(attr_id, len+1, buf)
    return String(buf)
end

function h5a_get_name_by_idx(loc_id, obj_name, idx_type, order, idx, lapl_id)
    len = h5a_get_name_by_idx(loc_id, obj_name, idx_type, order, idx, C_NULL, 0, lapl_id)
    buf = StringVector(len)
    h5a_get_name_by_idx(loc_id, obj_name, idx_type, order, idx, buf, len + 1, lapl_id)
    return String(buf)
end

###
### Dataset Interface
###

function h5d_vlen_get_buf_size(dataset_id, type_id, space_id)
    sz = Ref{hsize_t}()
    h5d_vlen_get_buf_size(dataset_id, type_id, space_id, sz)
    return sz[]
end

###
### Error Interface
###

function h5e_get_auto(error_stack)
    func = Ref{Ptr{Cvoid}}()
    client_data = Ref{Ptr{Cvoid}}()
    h5e_get_auto(error_stack, func, client_data)
    return func[], client_data[]
end

###
### File Interface
###

function h5f_get_intent(file_id)
    intent = Ref{Cuint}()
    h5f_get_intent(file_id, intent)
    return intent[]
end

function h5f_get_name(loc_id)
    len = h5f_get_name(loc_id, C_NULL, 0)
    buf = StringVector(len)
    h5f_get_name(loc_id, buf, len+1)
    return String(buf)
end

function h5f_get_obj_ids(file_id, types)
    sz = h5f_get_obj_count(file_id, types)
    hids = Vector{hid_t}(undef, sz)
    sz2 = h5f_get_obj_ids(file_id, types, sz, hids)
    sz2 != sz && resize!(hids, sz2)
    return hids
end

function h5f_get_vfd_handle(file_id, fapl)
    file_handle = Ref{Ptr{Cvoid}}()
    h5f_get_vfd_handle(file_id, fapl, file_handle)
    return file_handle[]
end

###
### Group Interface
###

function h5g_get_info(loc_id)
    ginfo = Ref{H5G_info_t}()
    h5g_get_info(loc_id, ginfo)
    return ginfo[]
end

function h5g_get_num_objs(loc_id)
    num_objs = Ref{hsize_t}()
    h5g_get_num_objs(loc_id, num_objs)
    return num_objs[]
end

###
### Identifier Interface
###

function h5i_get_name(loc_id)
    len = h5i_get_name(loc_id, C_NULL, 0)
    buf = StringVector(len)
    h5i_get_name(loc_id, buf, len+1)
    return String(buf)
end

###
### Link Interface
###

function h5l_get_info(link_loc_id, link_name, lapl_id)
    info = Ref{H5L_info_t}()
    h5l_get_info(link_loc_id, link_name, info, lapl_id)
    return info[]
end

function h5l_get_name_by_idx(loc_id, group_name, idx_type, order, idx, lapl_id)
    len = h5l_get_name_by_idx(loc_id, group_name, idx_type, order, idx, C_NULL, 0, lapl_id)
    buf = StringVector(len)
    h5l_get_name_by_idx(loc_id, group_name, idx_type, order, idx, buf, len + 1, lapl_id)
    return String(buf)
end

###
### Object Interface
###

function h5o_get_info(loc_id)
    oinfo = Ref{H5O_info_t}()
    h5o_get_info(loc_id, oinfo)
    return oinfo[]
end

###
### Property Interface
###

function h5p_get_alignment(fapl_id)
    threshold = Ref{hsize_t}()
    alignment = Ref{hsize_t}()
    h5p_get_alignment(fapl_id, threshold, alignment)
    return threshold[], alignment[]
end

function h5p_get_alloc_time(plist_id)
    alloc_time = Ref{Cint}()
    h5p_get_alloc_time(plist_id, alloc_time)
    return alloc_time[]
end

function h5p_get_char_encoding(plist_id)
    encoding = Ref{Cint}()
    h5p_get_char_encoding(plist_id, encoding)
    return encoding[]
end

function h5p_get_chunk(plist_id)
    n = h5p_get_chunk(plist_id, 0, C_NULL)
    cdims = Vector{hsize_t}(undef, n)
    h5p_get_chunk(plist_id, n, cdims)
    return cdims
end

function h5p_get_create_intermediate_group(plist_id)
    cig = Ref{Cuint}()
    h5p_get_create_intermediate_group(plist_id, cig)
    return cig[]
end

function h5p_get_dxpl_mpio(dxpl_id)
    xfer_mode = Ref{Cint}()
    h5p_get_dxpl_mpio(dxpl_id, xfer_mode)
    return xfer_mode[]
end

function h5p_get_fclose_degree(fapl_id)
    out = Ref{Cint}()
    h5p_get_fclose_degree(fapl_id, out)
    return out[]
end

function h5p_get_libver_bounds(plist_id)
    low = Ref{Cint}()
    high = Ref{Cint}()
    h5p_get_libver_bounds(plist_id, low, high)
    return low[], high[]
end

function h5p_get_local_heap_size_hint(plist_id)
    size_hint = Ref{Csize_t}()
    h5p_get_local_heap_size_hint(plist_id, size_hint)
    return size_hint[]
end

function h5p_get_obj_track_times(plist_id)
    track_times = Ref{UInt8}()
    h5p_get_obj_track_times(plist_id, track_times)
    return track_times[] != 0x0
end

function h5p_get_userblock(plist_id)
    len = Ref{hsize_t}()
    h5p_get_userblock(plist_id, len)
    return len[]
end

# Note: The following function(s) implement direct ccalls because the binding generator
# cannot (yet) do the string wrapping and memory freeing.

function h5p_get_class_name(pcid)
    pc = ccall((:H5Pget_class_name, libhdf5), Ptr{UInt8}, (hid_t,), pcid)
    if pc == C_NULL
        error("Error getting class name")
    end
    s = unsafe_string(pc)
    h5_free_memory(pc)
    return s
end

###
### Reference Interface
###

###
### Dataspace Interface
###

function h5s_get_simple_extent_dims(space_id)
    n = h5s_get_simple_extent_ndims(space_id)
    dims = Vector{hsize_t}(undef, n)
    maxdims = Vector{hsize_t}(undef, n)
    h5s_get_simple_extent_dims(space_id, dims, maxdims)
    return dims, maxdims
end

###
### Datatype Interface
###

function h5t_get_array_dims(type_id)
    nd = h5t_get_array_ndims(type_id)
    dims = Vector{hsize_t}(undef, nd)
    h5t_get_array_dims(type_id, dims)
    return dims
end

function h5t_get_fields(type_id)
    spos = Ref{Csize_t}()
    epos = Ref{Csize_t}()
    esize = Ref{Csize_t}()
    mpos = Ref{Csize_t}()
    msize = Ref{Csize_t}()
    h5t_get_fields(type_id, spos, epos, esize, mpos, msize)
    return (spos[], epos[], esize[], mpos[], msize[])
end

# Note: The following two functions implement direct ccalls because the binding generator
# cannot (yet) do the string wrapping and memory freeing.
function h5t_get_member_name(type_id, index)
    pn = ccall((:H5Tget_member_name, libhdf5), Ptr{UInt8}, (hid_t, Cuint), type_id, index)
    if pn == C_NULL
        error("Error getting name of compound datatype member #", index)
    end
    s = unsafe_string(pn)
    h5_free_memory(pn)
    return s
end

function h5t_get_tag(type_id)
    pc = ccall((:H5Tget_tag, libhdf5), Ptr{UInt8}, (hid_t,), type_id)
    if pc == C_NULL
        error("Error getting opaque tag")
    end
    s = unsafe_string(pc)
    h5_free_memory(pc)
    return s
end

###
### Optimized Functions Interface
###

###
### HDF5 Lite Interface
###

function h5lt_dtype_to_text(dtype_id)
    len = Ref{Csize_t}()
    h5lt_dtype_to_text(dtype_id, C_NULL, 0, len)
    buf = StringVector(len[] - 1)
    h5lt_dtype_to_text(dtype_id, buf, 0, len)
    return String(buf)
end

###
### Table Interface
###

###
### Filter Interface
###


text:
abstract type AbstractMode end
struct NoTransformations <: AbstractMode end
struct Squash <: AbstractMode end
struct FullSquash <: AbstractMode end
struct MDown <: AbstractMode end

struct MiniLogger{AM <: AbstractMode, IOT1 <: IO, IOT2 <: IO, DFT <: DateFormat} <: AbstractLogger
    io::IOT1
    ioerr::IOT2
    errlevel::LogLevel
    minlevel::LogLevel
    message_limits::Dict{Any,Int}
    flush::Bool
    format::Vector{Token}
    dtformat::DFT
    mode::AM
    squash_delimiter::String
    flush_threshold::Int
    lastflush::Base.RefValue{Int64}
end

getio(io, append) = io
getio(io::AbstractString, append) = open(io, append ? "a" : "w")

getflushthreshold(x::Integer) = x
getflushthreshold(x::TimePeriod) = Dates.value(Millisecond(x))

getmode(mode) = mode
getmode(mode::AbstractString) = getmode(Symbol(mode))
function getmode(mode::Symbol)
    if mode == :notransformations
        return NoTransformations()
    elseif mode == :squash
        return Squash()
    elseif mode == :fullsquash
        return FullSquash()
    elseif mode == :markdown
        return MDown()
    end
end

function Base.close(logger::MiniLogger)
    if logger.io != stdout && logger.io != stderr && isopen(logger.io)
        close(logger.io)
    end

    if logger.ioerr != stdout && logger.ioerr != stderr && isopen(logger.ioerr)
        close(logger.ioerr)
    end
end

"""
    MiniLogger(; <keyword arguments>)

MiniLogger constructor creates custom logger which can be used with usual `@info`, `@debug` commands.

Supported keyword arguments include:

* `io` (default `stdout`): IO stream which is used to output log messages below `errlevel` level. Can be either `IO` or `String`, in the latter case it is treated as a name of the output file.
* `ioerr` (default `stderr`): IO stream which is used to output log messages above `errlevel` level. Can be either `IO` or `String`, in the latter case it is treated as a name of the output file.
* `errlevel` (default `Error`): determines which output IO to use for log messages. If you want for all messages to go to `io`, set this parameter to `MiniLoggers.AboveMaxLevel`. If you want for all messages to go to `ioerr`, set this parameter to `MiniLoggers.BelowMinLevel`.
* `minlevel` (default: `Info`): messages below this level are ignored. For example with default setting `@debug "foo"` is ignored.
* `append` (default: `false`): defines whether to append to output stream or to truncate file initially. Used only if `io` or `ioerr` is a file path.
* `message_mode` (default: `:squash`): choose how message is transformed before being printed out. Following modes are supported:
    * `:notransformations`: message printed out as is, without any extra transformations
    * `:squash`: message is squashed to a single line, i.e. all `\\n` are changed to `squash_delimiter` and `\\r` are removed.
    * `:fullsquash`: all messages including error stacktraces are squashed to a single line, i.e. all `\\n` are changed to `squash_delimiter` and `\\r` are removed
    * `:markdown`: message is treated as if it is written in markdown
* `squash_delimiter`: (default: "\\t"): defines which delimiter to use in squash mode.
* `flush` (default: `true`): whether to `flush` IO stream for each log message. Flush behaviour also affected by `flush_threshold` argument.
* `flush_threshold::Union{Integer, TimePeriod}` (default: 0): if this argument is nonzero and `flush` is `true`, then `io` is flushed only once per `flush_threshold` milliseconds. I.e. if time between two consecutive log messages is less then `flush_threshold`, then second message is not flushed and will have to wait for the next log event.
* `dtformat` (default: "yyyy-mm-dd HH:MM:SS"): if `datetime` parameter is used in `format` argument, this dateformat is applied for output timestamps.
* `format` (default: "[{timestamp:func}] {level:func}: {message}"): format for output log message. It accepts following keywords, which should be provided in curly brackets:
    * `timestamp`: timestamp of the log message
    * `level`: name of log level (Debug, Info, etc)
    * `filepath`: filepath of the file, which produced log message
    * `basename`: basename of the filepath of the file, which produced log message
    * `line`: line number of the log command in the file, which produced log message
    * `group`: log group
    * `module`: name of the module, which contains log command
    * `id`: log message id
    * `message`: message itself

Each keyword accepts color information, which should be added after colon inside curly brackets. Colors can be either from `Base.text_colors` or special keyword `func`, in which case is used automated coloring. Additionaly, `bold` modifier is accepted by the `format` argument. For example: `{line:red}`, `{module:cyan:bold}`, `{group:func}` are all valid parts of the format command.

Colour information is applied recursively without override, so `{{line} {module:cyan} {group}:red}` is equivalent to `{line:red} {module:cyan} {group:red}`.

If part of the format is not a recognised keyword, then it is just used as is, for example `{foo:red}` means that output log message contain word "foo" printed in red.
"""
function MiniLogger(; io = stdout, ioerr = stderr, errlevel = Error, minlevel = Info, append = false, message_limits = Dict{Any, Int}(), flush = true, format = "[{timestamp:func}] {level:func}: {message}", dtformat = dateformat"yyyy-mm-dd HH:MM:SS", flush_threshold = 0, message_mode = Squash(), squash_delimiter = "\t")
    tio = getio(io, append)
    tioerr = io == ioerr ? tio : getio(ioerr, append)
    lastflush = Dates.value(Dates.now())
    MiniLogger(tio,
               tioerr,
               errlevel,
               minlevel,
               message_limits,
               flush,
               tokenize(format),
               dtformat,
               getmode(message_mode),
               squash_delimiter,
               getflushthreshold(flush_threshold),
               Ref(lastflush))
end

shouldlog(logger::MiniLogger, level, _module, group, id) =
    get(logger.message_limits, id, 1) > 0

min_enabled_level(logger::MiniLogger) = logger.minlevel

catch_exceptions(logger::MiniLogger) = true

# Formatting of values in key value pairs
squash(msg, logger, mode) = string(msg)
function squash(msg, logger, mode::Squash)
    smsg = string(msg)
    smsg = replace(smsg, "\r" => "")
    smsg = replace(smsg, "\n" => logger.squash_delimiter)
end

showvalue(io, msg, logger, mode) = print(io, squash(msg, logger, mode))

function showvalue(io, e::Tuple{Exception,Any}, logger, mode)
    ex, bt = e
    Base.showerror(io, ex, bt; backtrace = bt!==nothing)
end
showvalue(io, ex::Exception, logger, mode) = Base.showerror(io, ex)
showvalue(io, ex::AbstractVector{Union{Ptr{Nothing}, Base.InterpreterIP}}, logger, mode) = Base.show_backtrace(io, ex)

# Here we are fighting with multiple dispatch.
# If message is `Exception` or `Tuple{Exception, Any}` or anything else
# then we want to ignore third argument.
# But if it is any other sort of message we want to dispatch result
# on the type of message transformation
_showmessage(io, msg, logger, mode) = print(io, squash(msg, logger, mode))
_showmessage(io, msg, logger, ::MDown) = show(io, MIME"text/plain"(), Markdown.parse(msg))

showmessage(io, msg, logger, mode) = _showmessage(io, msg, logger, mode)
showmessage(io, e::Tuple{Exception,Any}, logger, mode) = showvalue(io, e, logger, mode)
showmessage(io, ex::Exception, logger, mode) = showvalue(io, ex, logger, mode)
showmessage(io, ex::AbstractVector{Union{Ptr{Nothing}, Base.InterpreterIP}}, logger, mode) = Base.show_backtrace(io, ex)

function postprocess(mode, delimiter, iobuf)
    print(iobuf, "\n")
    take!(iobuf)
end

function postprocess(mode::FullSquash, delimiter, iobuf)
    buf = take!(iobuf)
    delm = Vector{UInt8}(delimiter)
    res = similar(buf)
    L = length(res)
    j = 1
    @inbounds for (i, c) in pairs(buf)
        c == UInt8('\r') && continue
        if c == UInt8('\n')
            for c2 in delm
                if j > L
                    resize!(res, 2*L)
                    L *= 2
                end
                res[j] = c2
                j += 1
            end
            continue
        end
        res[j] = c
        j += 1
    end
    if j > L
        resize!(res, 2*L)
        L *= 2
    end
    res[j] = UInt8('\n')

    resize!(res, j)

    return res
end

tsnow(dtf) = Dates.format(Dates.now(), dtf)

function colorfunc(level::LogLevel, _module, group, id, filepath, line, element)
    level < Info  ? Color(Base.debug_color()) :
    level < Warn  ? Color(Base.info_color())  :
    level < Error ? Color(Base.warn_color())  :
                    Color(Base.error_color())
end

function extractcolor(token::Token, level, _module, group, id, filepath, line)
    if token.c.c == :func
        color = colorfunc(level, _module, group, id, filepath, line, token.val)
        return Color(color.c, color.isbold || token.c.isbold)
    else
        return token.c
    end
end

function printwcolor(iob, val, color)
    if color.c == -1
        if color.isbold
            printstyled(iob, val; bold = color.isbold)
        else
            print(iob, val)
        end
    else
        printstyled(iob, val; color = color.c, bold = color.isbold)
    end
end

function handle_message(logger::MiniLogger, level, message, _module, group, id,
                        filepath, line; maxlog=nothing, kwargs...)
    if maxlog !== nothing && maxlog isa Integer
        remaining = get!(logger.message_limits, id, maxlog)
        logger.message_limits[id] = remaining - 1
        remaining > 0 || return
    end

    io = if level < logger.errlevel
            isopen(logger.io) ? logger.io : stdout
        else
            isopen(logger.ioerr) ? logger.ioerr : stderr
        end

    buf = IOBuffer()
    iob = IOContext(buf, io)

    for token in logger.format
        c = extractcolor(token, level, _module, group, id, filepath, line)
        val = token.val
        if val == "timestamp"
            printwcolor(iob, tsnow(logger.dtformat), c)
        elseif val == "message"
            showmessage(iob, message, logger, logger.mode)

            iscomma = false
            for (key, val) in kwargs
                if string(key) == val
                    print(iob, key)
                    iscomma = false
                else
                    iscomma && print(iob, ", ")
                    print(iob, key, " = ")
                    showvalue(iob, val, logger, logger.mode)
                    iscomma = true
                end
            end
        elseif val == "level"
            levelstr = string(level)
            printwcolor(iob, levelstr, c)
        elseif val == "basename"
            printwcolor(iob, basename(filepath), c)
        elseif val == "filepath"
            printwcolor(iob, filepath, c)
        elseif val == "line"
            printwcolor(iob, line, c)
        elseif val == "group"
            printwcolor(iob, group, c)
        elseif val == "module"
            printwcolor(iob, _module, c)
        elseif val == "id"
            printwcolor(iob, id, c)
        else
            printwcolor(iob, val, c)
        end
    end
    write(io, postprocess(logger.mode, logger.squash_delimiter, buf))
    if logger.flush
        if logger.flush_threshold <= 0
            flush(io)
        else
            t = Dates.value(Dates.now())
            if t - logger.lastflush[] >= logger.flush_threshold
                logger.lastflush[] = t
                flush(io)
            end
        end
    end
    nothing
end

ext: jl


text:
"""

Helper infrastructure package to compile and sample models using Stan's `cmdstan`.
Not really intended to be called directly by a user.

# Extended help

Exports:
```Julia
* `CMDSTAN_HOME`                       : Fields common to all CmdStanModels
* `@shared_fields_stanmodels`          : Fields common to all CmdStanModels
* `CmdStanModels`                      : Supertype of CmdStanModels, see ?CmdStanModels
* `StanModelError`                     : Exception error on cmdstan compilation failure
* `HelpModel`                          : Simplests of the CmdStanModels
* `cmdline`                            : Inital portion of the cmdstan cammand line
* `stan_help`                          : Alias for calling stan_sample with a HelpModel
* `stan_sample`                        : Generic call to cmdstan, dispatched on CmdStanModel
* `read_summary`                       : Read the cmdstan summary .csv file
* `stan_summary`                       : Create the cmdstan summary .csv file
* `get_cmdstan_home``                  : Directory with cmdstan binaries
* `set_cmdstan_home!``                 : Update the cmdstan directory
* `findall`                            : Defined only for Julia versions < 1.3
```
"""
module StanBase

using DocStringExtensions: FIELDS, SIGNATURES, TYPEDEF
using Unicode, DelimitedFiles, Distributed
using Parameters

using StanDump

Int64(VERSION.minor) < 3 && include("utils/findall.jl")
include("stanmodel/cmdstan_home.jl")
include("stanmodel/shared_fields.jl")
include("stanmodel/top_level_types.jl")
include("stanmodel/help_types.jl")
include("stanmodel/HelpModel.jl")
include("stanmodel/update_model_file.jl")
include("stanmodel/number_of_chains.jl")
include("stanrun/executable.jl")
include("stanrun/cmdline.jl")
include("stanrun/stan_sample.jl")
include("stansamples/stan_summary.jl")
include("stansamples/read_summary.jl")
include("utils/par.jl")

"""
The directory which contains the cmdstan executables such as `bin/stanc` and
`bin/stansummary`.

# Extended help

Inferred from the environment variable `JULIA_CMDSTAN_HOME` or `ENV["JULIA_CMDSTAN_HOME"]`
when available.

If these are not available, use `set_cmdstan_home!` to set the value of CMDSTAN_HOME.

Example: `set_cmdstan_home!(homedir() * "/Projects/Stan/cmdstan/")`

Executing `versioninfo()` will display the value of `JULIA_CMDSTAN_HOME` if defined.
"""
CMDSTAN_HOME=""

function __init__()
  global CMDSTAN_HOME = if isdefined(Main, :JULIA_CMDSTAN_HOME)
    Main.JULIA_CMDSTAN_HOME
  elseif haskey(ENV, "JULIA_CMDSTAN_HOME")
    ENV["JULIA_CMDSTAN_HOME"]
  elseif haskey(ENV, "CMDSTAN_HOME")
    ENV["CMDSTAN_HOME"]
  else
    @warn("Environment variable CMDSTAN_HOME not set. Use set_cmdstan_home!.")
    ""
  end
end

"""

Alias for stan_sample(helpmodel,...)

# Extended help

see ?HelpModel
see ?stan_sample
"""
stan_help = stan_sample

export
  CMDSTAN_HOME,
  @shared_fields_stanmodels,
  CmdStanModels,
  StanModelError,
  HelpModel,
  cmdline,
  stan_help,
  stan_sample,
  read_summary,
  stan_summary,
  get_cmdstan_home,
  set_cmdstan_home!

end # module


text:
module TimeZones

using Base.Dates
import Base.Dates: TimeZone, AbstractTime, days, hour, minute, second, millisecond
import Compat: is_windows

export TimeZone, FixedTimeZone, VariableTimeZone, ZonedDateTime, DateTime,
    TimeError, AmbiguousTimeError, NonExistentTimeError, UnhandledTimeError,
    # discovery.jl
    timezone_names, all_timezones, timezones_from_abbr, timezone_abbrs,
    # accessors.jl
    hour, minute, second, millisecond,
    # adjusters.jl
    firstdayofweek, lastdayofweek,
    firstdayofmonth, lastdayofmonth,
    firstdayofyear, lastdayofyear,
    firstdayofquarter, lastdayofquarter,
    # Re-export from Base.Dates
    yearmonthday, yearmonth, monthday, year, month, week, day, dayofmonth,
    # conversion.jl
    now, astimezone,
    # local.jl
    localzone,
    # ranges.jl
    guess

const PKG_DIR = normpath(joinpath(dirname(@__FILE__), ".."))
const TZDATA_DIR = joinpath(PKG_DIR, "deps", "tzdata")
const COMPILED_DIR = joinpath(PKG_DIR, "deps", "compiled")
const TIME_ZONES = Dict{AbstractString,TimeZone}()

if is_windows()
    const WIN_TRANSLATION_FILE = joinpath(PKG_DIR, "deps", "windows_to_posix")
end

function __init__()
    # SLOT_RULE extension needs to happen everytime the module is loaded (issue #24)
    Base.Dates.SLOT_RULE['z'] = TimeZone
    Base.Dates.SLOT_RULE['Z'] = TimeZone

    global const ISOZonedDateTimeFormat = DateFormat("yyyy-mm-ddTHH:MM:SS.szzz")
end

include("utils.jl")
include("time.jl")
include("utcoffset.jl")
include("types.jl")
include("interpret.jl")
include("accessors.jl")
include("arithmetic.jl")
include("io.jl")
include("tzfile.jl")
include("adjusters.jl")
include("Olson.jl")
include("conversions.jl")
include("local.jl")
include("ranges.jl")
include("discovery.jl")
VERSION >= v"0.5.0-dev+5244" && include("rounding.jl")

"""
    TimeZone(name::AbstractString) -> TimeZone

Constructs a `TimeZone` instance based upon its `name`. A list of available time zones can
be determined using `timezone_names()`.

See `FixedTimeZone(::AbstractString)` for making a custom `TimeZone` instances.
"""
function TimeZone(name::AbstractString)
    return get!(TIME_ZONES, name) do
        tz_path = joinpath(COMPILED_DIR, split(name, "/")...)
        isfile(tz_path) || throw(ArgumentError("Unknown time zone named $name"))

        open(tz_path, "r") do fp
            return deserialize(fp)
        end
    end
end

end # module


text:
module Ipopt
using Libdl
using LinearAlgebra

if isfile(joinpath(dirname(@__FILE__),"..","deps","deps.jl"))
    include("../deps/deps.jl")
else
    error("Ipopt not properly installed. Please run Pkg.build(\"Ipopt\")")
end

function amplexefun(arguments::String)
    temp_env = copy(ENV)
    for var in Ipopt.amplexe_env_var
        temp_env[var] = Ipopt.amplexe_env_val
    end
    temp_dir = abspath(dirname(Ipopt.amplexe))
    proc = run(pipeline(Cmd(`$(Ipopt.amplexe) $arguments`,env=temp_env,dir=temp_dir), stdout=stdout))
    wait(proc)
    kill(proc)
    proc.exitcode
end

export createProblem, addOption
export openOutputFile, setProblemScaling, setIntermediateCallback
export solveProblem
export IpoptProblem

function __init__()
    julia_libdir = joinpath(dirname(first(filter(x -> occursin("libjulia", x), Libdl.dllist()))), "julia")
    julia_bindir = Sys.BINDIR
    ipopt_libdir = dirname(libipopt)
    ipopt_bindir = joinpath(dirname(libipopt), "..", "bin")
    pathsep = Sys.iswindows() ? ';' : ':'
    @static if Sys.isapple()
        global amplexe_env_var = ["DYLD_LIBRARY_PATH"]
        global amplexe_env_val = "$(julia_libdir)$(pathsep)$(get(ENV,"DYLD_LIBRARY_PATH",""))"
    elseif Sys.islinux()
        global amplexe_env_var = ["LD_LIBRARY_PATH"]
        global amplexe_env_val = "$(julia_libdir)$(pathsep)$(get(ENV,"LD_LIBRARY_PATH",""))"
    elseif Sys.iswindows()
        # for some reason windows sometimes needs Path instead of PATH
        global amplexe_env_var = ["PATH","Path","path"]
        global amplexe_env_val = "$(julia_bindir)$(pathsep)$(get(ENV,"PATH",""))"
    end

    # Still need this for AmplNLWriter to work until it uses amplexefun defined above
    # (amplexefun wraps the call to the binary and doesn't leave environment variables changed.)
    @static if Sys.isapple()
         ENV["DYLD_LIBRARY_PATH"] = string(get(ENV, "DYLD_LIBRARY_PATH", ""), ":", julia_libdir)
    elseif Sys.islinux()
         ENV["LD_LIBRARY_PATH"] = string(get(ENV, "LD_LIBRARY_PATH", ""), ":", julia_libdir, ":", ipopt_libdir)
    end
end

mutable struct IpoptProblem
    ref::Ptr{Cvoid}  # Reference to the internal data structure
    n::Int  # Num vars
    m::Int  # Num cons
    x::Vector{Float64}  # Starting and final solution
    g::Vector{Float64}  # Final constraint values
    mult_g::Vector{Float64} # lagrange multipliers on constraints
    mult_x_L::Vector{Float64} # lagrange multipliers on lower bounds
    mult_x_U::Vector{Float64} # lagrange multipliers on upper bounds
    obj_val::Float64  # Final objective
    status::Int  # Final status

    # Callbacks
    eval_f::Function
    eval_g::Function
    eval_grad_f::Function
    eval_jac_g::Function
    eval_h  # Can be nothing
    intermediate  # Can be nothing

    # For MathProgBase
    sense::Symbol

    function IpoptProblem(
        ref::Ptr{Cvoid}, n, m,
        eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
        prob = new(ref, n, m, zeros(Float64, n), zeros(Float64, m), zeros(Float64,m),
                   zeros(Float64,n), zeros(Float64,n), 0.0, 0,
                   eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h, nothing,
                   :Min)
        # Free the internal IpoptProblem structure when
        # the Julia IpoptProblem instance goes out of scope
        finalizer(freeProblem, prob)
        return prob
    end
end

# From Ipopt/src/Interfaces/IpReturnCodes_inc.h
ApplicationReturnStatus = Dict(
0=>:Solve_Succeeded,
1=>:Solved_To_Acceptable_Level,
2=>:Infeasible_Problem_Detected,
3=>:Search_Direction_Becomes_Too_Small,
4=>:Diverging_Iterates,
5=>:User_Requested_Stop,
6=>:Feasible_Point_Found,
-1=>:Maximum_Iterations_Exceeded,
-2=>:Restoration_Failed,
-3=>:Error_In_Step_Computation,
-4=>:Maximum_CpuTime_Exceeded,
-10=>:Not_Enough_Degrees_Of_Freedom,
-11=>:Invalid_Problem_Definition,
-12=>:Invalid_Option,
-13=>:Invalid_Number_Detected,
-100=>:Unrecoverable_Exception,
-101=>:NonIpopt_Exception_Thrown,
-102=>:Insufficient_Memory,
-199=>:Internal_Error)

###########################################################################
# Callback wrappers
###########################################################################
# Objective (eval_f)
function eval_f_wrapper(n::Cint, x_ptr::Ptr{Float64}, new_x::Cint, obj_ptr::Ptr{Float64}, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    # Calculate the new objective
    new_obj = convert(Float64, prob.eval_f(unsafe_wrap(Array,x_ptr, Int(n))))::Float64
    # Fill out the pointer
    unsafe_store!(obj_ptr, new_obj)
    # Done
    return Int32(1)
end

# Constraints (eval_g)
function eval_g_wrapper(n::Cint, x_ptr::Ptr{Float64}, new_x::Cint, m::Cint, g_ptr::Ptr{Float64}, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    # Calculate the new constraint values
    new_g = unsafe_wrap(Array,g_ptr, Int(m))
    prob.eval_g(unsafe_wrap(Array,x_ptr, Int(n)), new_g)
    # Done
    return Int32(1)
end

# Objective gradient (eval_grad_f)
function eval_grad_f_wrapper(n::Cint, x_ptr::Ptr{Float64}, new_x::Cint, grad_f_ptr::Ptr{Float64}, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    # Calculate the gradient
    new_grad_f = unsafe_wrap(Array,grad_f_ptr, Int(n))
    prob.eval_grad_f(unsafe_wrap(Array,x_ptr, Int(n)), new_grad_f)
    # Done
    return Int32(1)
end

# Jacobian (eval_jac_g)
function eval_jac_g_wrapper(n::Cint, x_ptr::Ptr{Float64}, new_x::Cint, m::Cint, nele_jac::Cint, iRow::Ptr{Cint}, jCol::Ptr{Cint}, values_ptr::Ptr{Float64}, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    # Determine mode
    mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
    x = unsafe_wrap(Array, x_ptr, Int(n))
    rows = unsafe_wrap(Array, iRow, Int(nele_jac))
    cols = unsafe_wrap(Array,jCol, Int(nele_jac))
    values = unsafe_wrap(Array,values_ptr, Int(nele_jac))
    prob.eval_jac_g(x, mode, rows, cols, values)
    # Done
    return Int32(1)
end

# Hessian
function eval_h_wrapper(n::Cint, x_ptr::Ptr{Float64}, new_x::Cint, obj_factor::Float64, m::Cint, lambda_ptr::Ptr{Float64}, new_lambda::Cint, nele_hess::Cint, iRow::Ptr{Cint}, jCol::Ptr{Cint}, values_ptr::Ptr{Float64}, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    # Did the user specify a Hessian
    if prob.eval_h === nothing
        # No Hessian provided
        return Int32(0)
    else
        # Determine mode
        mode = (values_ptr == C_NULL) ? (:Structure) : (:Values)
        x = unsafe_wrap(Array,x_ptr, Int(n))
        lambda = unsafe_wrap(Array,lambda_ptr, Int(m))
        rows = unsafe_wrap(Array,iRow, Int(nele_hess))
        cols = unsafe_wrap(Array,jCol, Int(nele_hess))
        values = unsafe_wrap(Array,values_ptr, Int(nele_hess))
        prob.eval_h(x, mode, rows, cols, obj_factor, lambda, values)
        # Done
        return Int32(1)
    end
end

# Intermediate
function intermediate_wrapper(alg_mod::Cint, iter_count::Cint, obj_value::Float64, inf_pr::Float64, inf_du::Float64, mu::Float64, d_norm::Float64, regularization_size::Float64, alpha_du::Float64, alpha_pr::Float64, ls_trials::Cint, user_data::Ptr{Cvoid})
    # Extract Julia the problem from the pointer
    prob = unsafe_pointer_to_objref(user_data)::IpoptProblem
    keepgoing = prob.intermediate(Int(alg_mod), Int(iter_count), obj_value, inf_pr, inf_du, mu, d_norm, regularization_size, alpha_du, alpha_pr, Int(ls_trials))
    # Done
    return keepgoing ? Int32(1) : Int32(0)
end

###########################################################################
# C function wrappers
###########################################################################
function createProblem(n::Int, x_L::Vector{Float64}, x_U::Vector{Float64},
    m::Int, g_L::Vector{Float64}, g_U::Vector{Float64},
    nele_jac::Int, nele_hess::Int,
    eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h = nothing)
    @assert n == length(x_L) == length(x_U)
    @assert m == length(g_L) == length(g_U)
    # Wrap callbacks
    eval_f_cb = @cfunction(eval_f_wrapper, Cint,
    (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cvoid}))
    eval_g_cb = @cfunction(eval_g_wrapper, Cint,
    (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Cvoid}))
    eval_grad_f_cb = @cfunction(eval_grad_f_wrapper, Cint,
    (Cint, Ptr{Float64}, Cint, Ptr{Float64}, Ptr{Cvoid}))
    eval_jac_g_cb = @cfunction(eval_jac_g_wrapper, Cint,
    (Cint, Ptr{Float64}, Cint, Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Float64}, Ptr{Cvoid}))
    eval_h_cb = @cfunction(eval_h_wrapper, Cint,
    (Cint, Ptr{Float64}, Cint, Float64, Cint, Ptr{Float64}, Cint, Cint, Ptr{Cint}, Ptr{Cint}, Ptr{Float64}, Ptr{Cvoid}))

    ret = ccall((:CreateIpoptProblem, libipopt), Ptr{Cvoid},
    (Cint, Ptr{Float64}, Ptr{Float64},  # Num vars, var lower and upper bounds
    Cint, Ptr{Float64}, Ptr{Float64},  # Num constraints, con lower and upper bounds
    Cint, Cint,                        # Num nnz in constraint Jacobian and in Hessian
    Cint,                              # 0 for C, 1 for Fortran
    Ptr{Cvoid}, Ptr{Cvoid},              # Callbacks for eval_f, eval_g
    Ptr{Cvoid}, Ptr{Cvoid}, Ptr{Cvoid}),  # Callbacks for eval_grad_f, eval_jac_g, eval_h
    n, x_L, x_U, m, g_L, g_U, nele_jac, nele_hess, 1,
    eval_f_cb, eval_g_cb, eval_grad_f_cb, eval_jac_g_cb, eval_h_cb)

    if ret == C_NULL
        error("IPOPT: Failed to construct problem.")
    else
        return IpoptProblem(ret, n, m, eval_f, eval_g, eval_grad_f, eval_jac_g, eval_h)
    end
end

# TODO: Not even expose this? Seems dangerous, should just destruct
# the IpoptProblem object via GC
function freeProblem(prob::IpoptProblem)
    if prob.ref != C_NULL
        ccall((:FreeIpoptProblem, libipopt), Cvoid, (Ptr{Cvoid},), prob.ref)
        prob.ref = C_NULL
    end
end

function addOption(prob::IpoptProblem, keyword::String, value::String)
    #/** Function for adding a string option.  Returns FALSE the option
    # *  could not be set (e.g., if keyword is unknown) */
    if !(isascii(keyword) && isascii(value))
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall((:AddIpoptStrOption, libipopt),
    Cint, (Ptr{Cvoid}, Ptr{UInt8}, Ptr{UInt8}),
    prob.ref, keyword, value)
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
end

function addOption(prob::IpoptProblem, keyword::String, value::Float64)
    #/** Function for adding a Number option.  Returns FALSE the option
    # *  could not be set (e.g., if keyword is unknown) */
    if !isascii(keyword)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall((:AddIpoptNumOption, libipopt),
    Cint, (Ptr{Cvoid}, Ptr{UInt8}, Float64),
    prob.ref, keyword, value)
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
end

function addOption(prob::IpoptProblem, keyword::String, value::Integer)
    #/** Function for adding an Int option.  Returns FALSE the option
    # *  could not be set (e.g., if keyword is unknown) */
    if !isascii(keyword)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall((:AddIpoptIntOption, libipopt),
    Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint),
    prob.ref, keyword, value)
    if ret == 0
        error("IPOPT: Couldn't set option '$keyword' to value '$value'.")
    end
end

function openOutputFile(prob::IpoptProblem, file_name::String, print_level::Int)
    #/** Function for opening an output file for a given name with given
    # *  printlevel.  Returns false, if there was a problem opening the
    # *  file. */
    if !isascii(file_name)
        error("IPOPT: Non ASCII parameters not supported")
    end
    ret = ccall((:OpenIpoptOutputFile, libipopt),
    Cint, (Ptr{Cvoid}, Ptr{UInt8}, Cint),
    prob.ref, file_name, print_level)
    if ret == 0
        error("IPOPT: Couldn't open output file.")
    end
end

# TODO: Verify this function even works! Trying it with 0.5 on HS071
# seems to change nothing.
function setProblemScaling(prob::IpoptProblem, obj_scaling::Float64,
    x_scaling = nothing,
    g_scaling = nothing)
    #/** Optional function for setting scaling parameter for the NLP.
    # *  This corresponds to the get_scaling_parameters method in TNLP.
    # *  If the pointers x_scaling or g_scaling are NULL, then no scaling
    # *  for x resp. g is done. */
    x_scale_arg = (x_scaling == nothing) ? C_NULL : x_scaling
    g_scale_arg = (g_scaling == nothing) ? C_NULL : g_scaling
    ret = ccall((:SetIpoptProblemScaling, libipopt),
    Cint, (Ptr{Cvoid}, Float64, Ptr{Float64}, Ptr{Float64}),
    prob.ref, obj_scaling, x_scale_arg, g_scale_arg)
    if ret == 0
        error("IPOPT: Error setting problem scaling.")
    end
end

function setIntermediateCallback(prob::IpoptProblem, intermediate::Function)
    intermediate_cb = @cfunction(intermediate_wrapper, Cint,
    (Cint, Cint, Float64, Float64, Float64, Float64,
    Float64, Float64, Float64, Float64, Cint, Ptr{Cvoid}))
    ret = ccall((:SetIntermediateCallback, libipopt), Cint,
    (Ptr{Cvoid}, Ptr{Cvoid}), prob.ref, intermediate_cb)
    prob.intermediate = intermediate
    if ret == 0
        error("IPOPT: Something went wrong setting the intermediate callback.")
    end
end

function solveProblem(prob::IpoptProblem)
    final_objval = [0.0]
    ret = ccall((:IpoptSolve, libipopt),
    Cint, (Ptr{Cvoid}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Ptr{Float64}, Any),
    prob.ref, prob.x, prob.g, final_objval, prob.mult_g, prob.mult_x_L, prob.mult_x_U, prob)
    prob.obj_val = final_objval[1]
    prob.status = Int(ret)

    return Int(ret)
end

include("MPB_wrapper.jl")
include("MOI_wrapper.jl")

end # module


text:
# __BEGIN_LICENSE__
#
# WaveOptics.jl
#
# Copyright (c) 2015, Stanford University
#
# All rights reserved.
#
# Redistribution and use in source and binary forms for academic and other
# non-commercial purposes with or without modification, are permitted provided
# that the following conditions are met:
#
# * Redistributions of source code, including modified source code, must retain
#   the above copyright notice, this list of conditions and the following
#   disclaimer.
#
# * Redistributions in binary form or a modified form of the source code must
#   reproduce the above copyright notice, this list of conditions and the
#   following disclaimer in the documentation and/or other materials provided with
#   the distribution.
#
# * Neither the name of The Leland Stanford Junior University, any of its
#   trademarks, the names of its employees, nor contributors to the source code
#   may be used to endorse or promote products derived from this software without
#   specific prior written permission.
#
# * Where a modified version of the source code is redistributed publicly in
#   source or binary forms, the modified source code must be published in a freely
#   accessible manner, or otherwise redistributed at no charge to anyone
#   requesting a copy of the modified source code, subject to the same terms as
#   this agreement.
#
# THIS SOFTWARE IS PROVIDED BY THE TRUSTEES OF THE LELAND STANFORD JUNIOR
# UNIVERSITY "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE LELAND STANFORD JUNIOR
# UNIVERSITY OR ITS TRUSTEES BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
# PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR
# BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
# CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING
# IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY
# OF SUCH DAMAGE.
#
# __END_LICENSE__

# -----------------------------------------------------------------------------
#               BACK APERTURE PLANE COORDINATE SYSTEMS
# -----------------------------------------------------------------------------

# Types for fourier_plane_meshgrid()
struct CartesianCoordinateSystem end
struct PolarCoordinateSystem end
struct BothCoordinateSystem end

abstract type OpticalCoordinateUnits end
struct SpatialFrequencyUnits <: OpticalCoordinateUnits end
struct MetricUnits <: OpticalCoordinateUnits end
struct NormalizedUnits <: OpticalCoordinateUnits end

mutable struct OpticalCoordinates
    fx::Array{Float32,2}
    fy::Array{Float32,2}
    rho::Array{Float32,2}
    theta::Array{Float32,2}
    units::OpticalCoordinateUnits
end

function _fourier_plane_meshgrid_freqs(sim_size_px::Int64, sample_period::Float64,
                                       objective_na::Float64, wavelength::Float64, focal_length::Float64,
                                       units::SpatialFrequencyUnits)
    return fftshift(fftfreq(sim_size_px, sample_period))
end

function _fourier_plane_meshgrid_freqs(sim_size_px::Int64, sample_period::Float64,
                                       objective_na::Float64, wavelength::Float64, focal_length::Float64,
                                       units::MetricUnits)
    freqs = fftshift(fftfreq(sim_size_px, sample_period))
    freqs *= focal_len * wavelength
    return freqs
end

function _fourier_plane_meshgrid_freqs(sim_size_px::Int64, sample_period::Float64,
                                       objective_na::Float64, wavelength::Float64, focal_length::Float64,
                                       units::NormalizedUnits)
    freqs = fftshift(fftfreq(sim_size_px, sample_period))
    back_aperture_radius = objective_na / wavelength  # From equation 21-4 in Gross
    freqs /= back_aperture_radius
    return freqs
end

function _fourier_plane_meshgrid_mesh(freqs, coordinate_system::CartesianCoordinateSystem)
    fx, fy = WaveOptics.meshgrid(freqs, freqs)
    return (fx, fy)
end

function _fourier_plane_meshgrid_mesh(freqs, coordinate_system::PolarCoordinateSystem)
    fx, fy = WaveOptics.meshgrid(freqs, freqs)
    r = sqrt.(fx.*fx+fy.*fy)
    theta = atan.(fy, fx)
    return (r, theta)
end

function _fourier_plane_meshgrid_mesh(freqs, coordinate_system::BothCoordinateSystem)
    fx, fy = WaveOptics.meshgrid(freqs, freqs)
    r = sqrt.(fx.*fx+fy.*fy)
    theta = atan.(fy, fx)
    return (fx, fy, r, theta)
end

"""
Returns a 2D grid of spatial frequencies (in both cartesian and polar
coordinates). Frequencies are normalized such that `R == 1` is along the edge
of the aperture defined by `objective_na / wavelength`.

Many functions will repeatedly re-use the coordinates produced by
fourier_plane_meshgrid(), so we will momoize the results so that they are
cached if the arguments to the function remain the same again and again.

**Parameters**

`sim_size_px` - the size of the simulation in the object (or image) plane

`sample_period` - the period (i.e. sample spacing) in the object (or image) plane

`objective_na, wavelength` - numerical aperture and emission wavelength

`focal_length` - focal length of the objective (or tubelens).  Used when units = "meters"

`units` - string
    one of : `[ "spatial_freq", "m", "normalized" ]`

* `"spatial_freq"` - return back aperture coordinates in terms of the spatial frequency of waves in the object plane.
* `"meters"` - returns back aperture coordinates in terms of their real
* `"normalized"` - returns back aperture coordinates in a normalized coordinate space where the aperture lies on the unit disc, |r| <= 1.

`coordinate_system` - string
    one of : `[ "cartesian", "polar", "both" ]`

**Returns**

`(fx, fy)` when `coordinate_system == "cartesian"` or `(r, theta)` when `coordinate_system == "polar"`

"""
function fourier_plane_meshgrid(sim_size_px::Int64, sample_period::Float64,
                                objective_na::Float64, wavelength::Float64,
                                focal_length::Float64, units::S) where S<:OpticalCoordinateUnits
    freqs = _fourier_plane_meshgrid_freqs(sim_size_px, sample_period, objective_na,
                                          wavelength, focal_length, units)
    (fx,fy,rho,theta) = _fourier_plane_meshgrid_mesh(freqs, BothCoordinateSystem())
    return OpticalCoordinates(fx, fy, rho, theta, units)
end


text:
# Autogenerated wrapper script for blaspp_jll for i686-w64-mingw32-cxx11
export libblaspp

using CompilerSupportLibraries_jll
using OpenBLAS32_jll
JLLWrappers.@generate_wrapper_header("blaspp")
JLLWrappers.@declare_library_product(libblaspp, "libblaspp.dll")
function __init__()
    JLLWrappers.@generate_init_header(CompilerSupportLibraries_jll, OpenBLAS32_jll)
    JLLWrappers.@init_library_product(
        libblaspp,
        "bin\\libblaspp.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
    randtree(n::Int, rooted::Bool = false, labels=string.(1:n))

Create a binary tree by the random addition of `n` tips.
"""
function _randtree_rooted(
	n::Int,
	labels::Vector{String},
	spp::Vector{Int}
	)

	nodes = Vector{Node}(undef, 2n - 1)
	origin = Node()
	nodes[1:3] .= origin, Node(labels[1]), Node(labels[2])
	link!(origin, nodes[2])
	link!(origin, nodes[3])
	nodes[2].species = spp[1]
	nodes[3].species = spp[2]

	for i ∈ 3:n
		p = Node(labels[i])
		p.species = spp[i]
		q = rand(nodes[1:(2i - 3)])
		r = rand(neighbours(q))

		s = Node()
		link!(p, s)

		graft!(s, q, r)

		nodes[(2i - 2):(2i - 1)] .= r, s
	end

	return origin
end

function _randtree_unrooted(n::Int, labels::Vector{String}, spp::Vector{Int})
	# TODO: create all the nodes in `nodes` first so that the rest of the function just creates links and sets labels.
	origin = Node()
	nodes = Vector{Node}(undef, 2n - 2)
	nodes[1:4] .= origin, Node(labels[1]), Node(labels[2]), Node(labels[3])
	nodes[2].species = spp[1]
	nodes[3].species = spp[2]
	nodes[4].species = spp[3]

	link!(origin, nodes[2])
	link!(origin, nodes[3])
	link!(origin, nodes[4])

	for i ∈ 4:n
		p = Node(labels[i])
		p.species = spp[i]
		q = rand(nodes[1:(2i - 4)])
		r = rand(neighbours(q))

		s = Node()
		link!(s, p)

		graft!(s, q, r)

		nodes[(2i - 3):(2i - 2)] .= p, s
	end

	return origin
end

"""
	randtree(n, rooted = false; labels, spp)

Create a binary tree (un`rooted` by default) by the random addition of `n` tips.

The tip `labels` can be provided as a vector of strings, or else the tips will be labelled in sequence from 1 to `n`. The tip `species` can be given as a vector of positive integers, or else they will all be 0.
"""
function randtree(
	n::Int,
	rooted::Bool=false;
	labels::Vector{String}=string.(1:n),
	species::Vector{Int}=zeros(Int, n)
	)
	@assert n == length(labels) == length(species)
	if n < 2
		@error "at least two tips are required for a rooted binary tree, or three for an unrooted binary tree."
	elseif n == 2
		rooted || @error "cannot create an unrooted binary tree with fewer than 3 tips."

		origin = Node()
		tip1 = Node(labels[1])
		tip1.species = species[1]
		tip2 = Node(labels[2])
		tip2.species = species[2]

		# Although the topologies are identical, we are going to randomise the order of the tips.
		if rand(Bool)
			link!(origin, tip1)
			link!(origin, tip2)
		else
			link!(origin, tip2)
			link!(origin, tip1)
		end

	elseif n == 3
		if rooted
			origin = _randtree_rooted(n, labels, species)
		else
			origin = Node()
			tips = map(Node, labels)
			for i in randperm(3)
				link!(origin, tips[i])
				tips[i].species = species[i]
			end
		end
	else
		origin = rooted ? _randtree_rooted(n, labels, species) :
			_randtree_unrooted(n, labels, species)
	end
	tree = Tree(origin)
	rooted && root!(tree)

	return tree
end

"""
	randtree(labels, rooted = false)

Create a binary tree (un`rooted` by default) by the random addition of tips with the given `labels`.
"""
randtree(labels::Vector{String}, rooted::Bool = false) =
	randtree(length(labels), rooted; labels = labels)

"""
	randtree(dir, rooted = false)

Create a binary tree (un`rooted` by default) by the random addition of tips with the species from a species `dir`ectory.
"""
function randtree(dir::SpeciesDirectory, rooted::Bool = false)
	tree =
		randtree(length(dir), rooted; labels = dir.list, species = collect(1:length(dir)))
	tree.dir = dir

	return tree
end


text:
AsocList{K, V} = AbstractArray{Tuple{K, V}}

function _genlex(argsym, retsym, lexer_table, reserved_words)
    gen_many_lexers = map(lexer_table) do (k, v)
        lnode = LineNumberNode(1, "lexing rule: $k")

        token_type = reserved_words === nothing ? QuoteNode(k) :
                        :(s in $reserved_words ? :reserved : $(QuoteNode(k)))
        quote
            $lnode
            s = $v($argsym, offset)
            if s !== nothing
            line_inc = count(x -> x === '\n', s)
            n = ncodeunits(s)
            push!($retsym, $Token{$token_type}(lineno, colno, offset, s, n))
            if  line_inc === 0
                colno += n
            else
                lineno += line_inc
                colno =  n - findlast(x -> x === '\n', s) + 1
            end
            offset += n
            continue
            end
        end
    end
    final = :($throw((SubString($argsym, offset), "Token NotFound")))
    Expr(:block, gen_many_lexers..., final)
end

rmlines = @λ begin
    e :: Expr           -> Expr(e.head, filter(x -> x !== nothing, map(rmlines, e.args))...)
      :: LineNumberNode -> nothing
    a                   -> a
end

struct LexerSpec{K}
    a :: K
end

function genlex(lexer_table, reserved_words)
    ex = quote
        function (tokens)
            lineno :: $Int64 = 1
            colno  :: $Int32 = 1
            offset :: $Int64  = 1
            ret = $Token[]
            N :: $Int64 = $ncodeunits(tokens)
            while offset <= N
                $(_genlex(:tokens, :ret, lexer_table, reserved_words))
            end
            ret
        end
    end
    # println(rmlines(ex))
    ex
end

macro genlex(lexer_table)
    genlex(lexer_table) |> esc
end

function mklexer(a :: LexerSpec{Char})
    quote (chars, i) -> chars[i] === $(a.a) ? String([$(a.a)]) : nothing end
end

function mklexer(a :: LexerSpec{String})
    let str = a.a, n = ncodeunits(str)
        quote (chars, i) ->
            startswith(SubString(chars, i), $str) ? $str : nothing
        end
    end
end

function mklexer(a :: LexerSpec{Regex})
    let regex = a.a
        quote
            function (chars, i)
                v = match($regex, SubString(chars, i))
                v === nothing ? nothing : v.match
            end
        end
    end
end

struct Quoted
    left :: String
    right :: String
    escape :: String
end

function mklexer(a :: LexerSpec{Quoted})
    let left = a.a.left,
        right = a.a.right,
        escape = a.a.escape
        quote
            function (chars, i)
                off = i
                n = $ncodeunits(chars)
                subs = SubString(chars, off)
                if $startswith(subs, $left)
                    off = off + $(ncodeunits(left))
                else
                    return nothing
                end
                while off <= n
                    subs = SubString(chars, off)
                    if $startswith(subs, $right)
                        off = off + $(ncodeunits(right))
                        # so tricky here for the impl of Julia string.
                        return chars[i:prevind(chars, off)]
                    end
                    if $startswith(subs, $escape)
                        off = off + $(ncodeunits(escape))
                    else
                        off = nextind(chars, off)
                    end
                end
                nothing
            end
        end
    end
end


text:
using Exercise1
using Documenter

DocMeta.setdocmeta!(Exercise1, :DocTestSetup, :(using Exercise1); recursive=true)

makedocs(;
    modules=[Exercise1],
    authors="Paul Schrimpf <paul.schrimpf@gmail.com> and contributors",
    repo="https://github.com/schrimpf/Exercise1.jl/blob/{commit}{path}#{line}",
    sitename="Exercise1.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://schrimpf.github.io/Exercise1.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Benchmarks" => "benchmarks.md"
    ],
)

deploydocs(;
    repo="github.com/schrimpf/Exercise1.jl",
    devbranch="main",
)


text:
import SpecialFunctions: erfc, erfcinv

"""
inverse of cdf.
"""
@inline function ppf(x)
    return -  sqrt(2) * erfcinv(2 * x)
end

@inline function pdf(x)
    return (1. / (sqrt(2 * pi) ) *
        exp(-(((x)) ^ 2 / 2.0)))
end

@inline function cdf(x)
    return 0.5 * erfc(-(x) / (sqrt(2)))
end


text:
using Base.Meta: isexpr
using ASTInterpreter2: JuliaStackFrame
using ASTInterpreter2: pc_expr, plain, evaluate_call!, finish_and_return!, @lookup

# Steps through the whole expression using `s`
function step_through(frame)
    state = DebuggerFramework.dummy_state([frame])
    while !isexpr(plain(pc_expr(state.stack[end])), :return)
        execute_command(state, state.stack[1], Val{:s}(), "s")
    end
    lastframe = state.stack[end]
    return @lookup(lastframe, plain(pc_expr(lastframe)).args[1])
end

# Execute a frame using Julia's regular compiled-code dispatch for any :call expressions
runframe(frame, pc=frame.pc[]) = Some{Any}(finish_and_return!(Compiled(), frame, pc))

# Execute a frame using the interpreter for all :call expressions (except builtins & intrinsics)
function runstack(frame::JuliaStackFrame, pc=frame.pc[])
    stack = JuliaStackFrame[]
    return Some{Any}(finish_and_return!(stack, frame, pc))
end


text:
abstract AbstractProposalDensity

abstract SymmetricDensity <: AbstractProposalDensity
abstract ASymmetricDensity <: AbstractProposalDensity

issymmetric{D<:SymmetricDensity}(::Union{D,Type{D}}) = true
issymmetric{D<:ASymmetricDensity}(::Union{D,Type{D}}) = false

###Factory function for known distributions (currently :normal and :lognormal)
density(s::Symbol,args...;keyargs...) = _density(Val{s},args...;keyargs...)

###For some densities, conditioning on a point can be faster than creating the density
###But in the general case, we don't know how to condition a density so throw an error
condition!(d::AbstractProposalDensity,x::AbstractVector) = throw(MethodError(condition!,(d,x)))
scale!(d::AbstractProposalDensity,s::AbstractFloat) = throw(MethodError(scale!,(d,s)))
update!(d::AbstractProposalDensity,x::AbstractVector,s::AbstractArray) = throw(MethodError(update!,(d,x,s)))

###For generic proposal densities, we don't know how to propose points, so throw an error
propose!(d::AbstractProposalDensity,v::AbstractArray) = throw(MethodError(propose!, (d,v)))

###For generic proposal densities, we don't know how to calculate the logprobability of a sample, so throw an error
logprobability!(r::AbstractVector,d::AbstractProposalDensity,v::AbstractArray) = throw(MethodError(logprobability!, (r,d,v)))
logprobability(d::AbstractProposalDensity,v::AbstractArray) = throw(MethodError(logprobability, (d,v)))

function show(io::IO,d::AbstractProposalDensity)
    println(io,"$(typeof(d)) with fields: $(fieldnames(d))")
    nothing
end


text:
mutable struct Position
    x::Int
    y::Int
end

struct Size
    width::Int
    height::Int
end

# Type of weapon
@enum Weapon Laser Missile

# A thing should have a position and size
abstract type Thing end

position(t::Thing) = t.position
size(t::Thing) = t.size
shape(t::Thing) = :unknown

struct Spaceship <: Thing
    position::Position
    size::Size
    weapon::Weapon
end
shape(s::Spaceship) = :saucer

Base.show(io::IO, s::Spaceship) =
    print(io, "Spaceship (", s.position.x, ",", s.position.y, ") ",
          s.size.width, "x", s.size.height, "/", s.weapon)

struct Asteroid <: Thing
    position::Position
    size::Size
end

Base.show(io::IO, a::Asteroid) =
    print(io, "Asteroid (", a.position.x, ",", a.position.y,
          ") ", a.size.width, "x", a.size.height)

# Test
s1 = Spaceship(Position(0,0), Size(30,5), Missile)
s2 = Spaceship(Position(10,0), Size(30,5), Laser)
a1 = Asteroid(Position(20,0), Size(20,20))
a2 = Asteroid(Position(0,20), Size(20,20))

# --------------------------------------------------------

function explode(things::AbstractVector{Any})
    for t in things
        println("Exploding $t")
    end
end

function explode(things::AbstractVector{T}) where {T <: Thing}
    for t in things
        println("Exploding thing => $t")
    end
end

# specifying abstract/concrete types in method signature
function tow(A::Spaceship, B::Thing)
    "tow 1"
end

# equivalent of parametric type
function tow(A::Spaceship, B::T) where {T <: Thing}
    "tow 2"
end

function group_anything(A::Thing, B::Thing)
    println("Grouped $A and $B")
end

function group_same_things(A::T, B::T) where {T <: Thing}
    println("Grouped $A and $B")
end

eltype(things::AbstractVector{T}) where {T <: Thing} = T


text:
if v"1.0" <= VERSION < v"1.1"
    isnothing(x) = x===nothing
end

if VERSION < v"1.4"
    pkgdir(pkg) = abspath(@__DIR__, "..")
end


text:
module HMC
using Distributions

"""
Get random momentum from a parameter.
"""
rand_momentum(x) = randn(size(x))

"""
Compute kinetic energy given momentum
"""
compute_kinetic_energy(momentum) = sum(momentum .^ 2) / 2

"""
See p. 14 of
https://arxiv.org/pdf/1206.1901.pdf
"""
# TODO: Check this
function hmc_update(curr_state,
                    logprob::Function, grad_logprob::Function,
                    num_leapfrog_steps::Int, eps::Float64;
                    momentum_sd::Real=1)
  # Position variables
  qs = copy(curr_state)

  # Momentum variables
  ps = rand_momentum(qs)

  # Current kinetic energy
  curr_K = compute_kinetic_energy(ps)

  # Potential enerdy
  U(s) = -logprob(s)

  # Current potential energy
  curr_U = U(qs)

  # Compute gradient of U
  grad_U(qs) = -grad_logprob(qs)

  # Make a half step for momentum at the beginning
  gs = grad_U(qs)
  ps .-= eps * qs / 2

  # Alternate full steps for position and momentum
  for i in 1:num_leapfrog_steps
    # Make a full step for the position
    qs .+= eps * ps

    # Make a full step for the momentum, except at end of trajectory
    if i < num_leapfrog_steps
      gs = grad_U(qs)
      ps .-= eps * qs
    end
  end

  # Make a half step for momentum at the end
  gs = grad_U(qs)
  ps .-= eps * qs / 2

  # Proposed potential energy
  cand_U = U(qs)

  # Proposed kinetic  energy
  cand_K = compute_kinetic_energy(ps)

  ### Metropolis step ###
  # log acceptance ratio.  NOTE: potential = -log_prob
  log_acceptance_ratio = curr_U + curr_K - cand_U - cand_K

  # Accept or reject
  if log_acceptance_ratio > log(rand())
    return qs, -cand_U
  else
    return curr_state, -curr_U
  end
  ### End of Metropolis step ###
end

end


text:
#!/usr/bin/env julia

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

# using RobotOS
# @rosimport barc.msg: ECU, pos_info, Vel_est
# @rosimport geometry_msgs.msg: Vector3
# @rosimport sensor_msgs.msg: Imu
# @rosimport marvelmind_nav.msg: hedge_pos
# rostypegen()
# using barc.msg
# using geometry_msgs.msg
# using sensor_msgs.msg
# using std_msgs.msg
# using marvelmind_nav.msg
# using JLD
# include("barc_lib/log_functions.jl")

# function main()

#     # Set up logging
#     buffersize      = 60000
#     gps_meas        = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
#     imu_meas        = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,9))
#     cmd_log         = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
#     cmd_pwm_log     = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
#     vel_est_log     = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,5))
#     pos_info_log    = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,20))

#     # Initialize ROS node and topics
#     init_node("mpc_traj")
#     s1  = Subscriber("ecu", ECU, ECU_callback, (cmd_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.ECU}
#     s2  = Subscriber("imu/data", Imu, IMU_callback, (imu_meas,), queue_size=1)::RobotOS.Subscriber{sensor_msgs.msg.Imu}
#     s3  = Subscriber("ecu_pwm", ECU, ECU_PWM_callback, (cmd_pwm_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.ECU}
#     s4  = Subscriber("hedge_pos", hedge_pos, GPS_callback, (gps_meas,), queue_size=1)::RobotOS.Subscriber{marvelmind_nav.msg.hedge_pos}
#     s5  = Subscriber("pos_info", pos_info, pos_info_callback, (pos_info_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.pos_info}
#     s6  = Subscriber("vel_est", Vel_est, vel_est_callback, (vel_est_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.Vel_est}

#     run_id      = get_param("run_id")

#     loop_rate   = Rate(10)
#     pub         = Publisher("ecu_pwm", ECU, queue_size=1)::RobotOS.Publisher{barc.msg.ECU}
#     cmd         = ECU()      # command type
#     cmd.motor   = 0
#     cmd.servo   = 0
#     cmd.header.stamp = get_rostime()
#     println("Starting open loop control. Waiting.")
#     rossleep(3.0)

#     t0_ros      = get_rostime()
#     t0          = to_sec(t0_ros)
#     t           = 0.0

#     gps_meas.i      = 1
#     imu_meas.i      = 1
#     cmd_log.i       = 1
#     cmd_pwm_log.i   = 1
#     vel_est_log.i   = 1
#     pos_info_log.i  = 1

#     chicane_speed = 97#1.0
#     chicane_turn  = [120,60]

#     cmd_m = 94
#     cmd_s = 80
#     t_next = 7
#     mode = 1        # 1 = acc, 2 = break, 3 = break
#     # Idea: 5 seconds acc, 3 seconds breaking.

#     # Start node
#     while ! is_shutdown()#t < 29.0              # run exactly x seconds
#         t = to_sec(get_rostime())-t0
#         if t <= 3
#             cmd.motor = 90.0
#             cmd.servo = 80.0
#         elseif t <= t_next
#             if mode == 1
#                 cmd.motor = cmd_m
#                 cmd.servo = cmd_s
#             elseif mode == 2
#                 cmd.motor = 70
#             end
#         elseif t <= 300
#             if mode == 1
#                 mode = 2
#             else
#                 cmd_s += 0
#                 cmd_m += 1
#                 mode = 1
#             end
#             t_next = t + 4
#         # **************************************
#         # CHICANE:
#         # elseif t <= 3
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[1]
#         # elseif t <= 4
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[2]
#         # elseif t <= 5
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[1]
#         # elseif t <= 6
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[2]
#         # elseif t <= 7
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[1]
#         # elseif t <= 8
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[2]
#         # elseif t <= 9
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[1]
#         # elseif t <= 10
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[2]
#         # elseif t <= 11
#         #     cmd.motor = chicane_speed
#         #     cmd.servo = chicane_turn[1]
#         # **************************************
#         # CONTINUOUS ACCELERATION:
#         # elseif t <= 123
#         #     cmd.motor = 0.0+(t-3)/20
#         #     cmd.servo = 0.0#-(t-3.0)/300-0.15
#         # elseif t <= 8                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 97.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = 100.0
#         # elseif t <= 13                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 97.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = 105.0
#         # elseif t <= 18                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 97.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = 110.0
#         # elseif t <= 23                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 97.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = 115.0
#         # elseif t <= 28                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 97.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = 120.0
#         # elseif t <= 33                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 1.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = -0.4
#         # elseif t <= 38                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 1.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = -0.45
#         # elseif t <= 43                   # CHECK TIME AND ACCELERATION !!!
#         #     cmd.motor = 1.0             # CHECK TIME AND ACCELERATION !!!
#         #     cmd.servo = -0.5
#         # elseif t <= 44
#         #     cmd.motor = -2.0
#         #     cmd.servo = 0
#         else
#             cmd.motor = 70.0
#             cmd.servo = 90.0
#         end
#         cmd.header.stamp = get_rostime()
#         publish(pub, cmd)

#         rossleep(loop_rate)
#     end

#     # Clean up buffers
#     clean_up(gps_meas)
#     clean_up(imu_meas)
#     clean_up(cmd_log)
#     clean_up(cmd_pwm_log)
#     clean_up(pos_info_log)
#     clean_up(vel_est_log)

#     println("Exiting open loop control.")
#     log_path = "$(homedir())/open_loop/output-record-$(run_id[1:4]).jld"
#     if isfile(log_path)
#         log_path = "$(homedir())/open_loop/output-record-$(run_id[1:4])-2.jld"
#         warn("Warning: File already exists.")
#     end
#     save(log_path,"gps_meas",gps_meas,"imu_meas",imu_meas,"cmd_log",cmd_log,"cmd_pwm_log",cmd_pwm_log,"pos_info",pos_info_log,"vel_est",vel_est_log)
#     println("Exiting node... Saved recorded data to $log_path.")
# end

# if ! isinteractive()
#     main()
# end

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

using RobotOS
@rosimport barc.msg: ECU, pos_info, Vel_est
@rosimport geometry_msgs.msg: Vector3
@rosimport sensor_msgs.msg: Imu
@rosimport marvelmind_nav.msg: hedge_pos
rostypegen()
using barc.msg
using geometry_msgs.msg
using sensor_msgs.msg
using std_msgs.msg
using marvelmind_nav.msg
using PyPlot
using JLD
include("barc_lib/log_functions.jl")

function main()

    # Set up logging
    buffersize      = 60000
    gps_meas        = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
    imu_meas        = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,9))
    cmd_log         = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
    cmd_pwm_log     = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,2))
    vel_est_log     = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,5))
    pos_info_log    = Measurements{Float64}(1,zeros(buffersize),zeros(buffersize),zeros(buffersize,20))

    # Initialize ROS node and topics
    init_node("mpc_traj_openloop")
    s1  = Subscriber("ecu", ECU, ECU_callback, (cmd_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.ECU}
    s2  = Subscriber("imu/data", Imu, IMU_callback, (imu_meas,), queue_size=1)::RobotOS.Subscriber{sensor_msgs.msg.Imu}
    s3  = Subscriber("ecu_pwm", ECU, ECU_PWM_callback, (cmd_pwm_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.ECU}
    s4  = Subscriber("hedge_pos", hedge_pos, GPS_callback, (gps_meas,), queue_size=1)::RobotOS.Subscriber{marvelmind_nav.msg.hedge_pos}
    s5  = Subscriber("pos_info", pos_info, pos_info_callback, (pos_info_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.pos_info}
    s6  = Subscriber("vel_est", Vel_est, vel_est_callback, (vel_est_log,), queue_size=1)::RobotOS.Subscriber{barc.msg.Vel_est}

    run_id      = get_param("run_id")

    loop_rate   = Rate(10)
    pub         = Publisher("ecu_pwm", ECU, queue_size=1)::RobotOS.Publisher{barc.msg.ECU}
    cmd         = ECU()      # command type
    cmd.motor   = 0
    cmd.servo   = 0
    cmd.header.stamp = get_rostime()
    println("Starting open loop control. Waiting.")
    rossleep(3.0)

    t0_ros      = get_rostime()
    t0          = to_sec(t0_ros)
    t           = 0.0

    t_next = 2 # SET THIS TO 4 IF YOU DO ACCELERATION TESTS

    gps_meas.i      = 1
    imu_meas.i      = 1
    cmd_log.i       = 1
    cmd_pwm_log.i   = 1
    vel_est_log.i   = 1
    pos_info_log.i  = 1

    cmd_m = 100 # PWM COMMAND FOR MOTOR

    cmd_s = 93  # PWM COMMAND FOR STEERING

    # t_cmd = zeros(10)
    # t_no_cmd = zeros(10)

    # index=1

    # time_index = 0

    # a=zeros(10)

    while ! is_shutdown()             # run exactly x seconds
        t = to_sec(get_rostime())-t0

        if t <= t_next

            cmd.motor = cmd_m
            cmd.servo = cmd_s

            #t_cmd[index] = t

        elseif t > t_next && t <= t_next + 4

            cmd.motor = 40

            #t_no_cmd[index] = t

        else
            cmd.motor = -1

            println("******FLAG*******")
        #     t_next = t + 4

        #     println("t== ",t)
        #     println("time_index== ",time_index)

        #     time_vector = time_index+0.1:.02:t-0.1

        #     v = zeros(length(time_vector))
        #     command = zeros(length(time_vector))

        #      println("v= ",v)

        #     # println("time vector= ",time_vector)

        #     for i=1:length(time_vector)

        #         v[i] = vel_est_log.z[time_vector[i].>vel_est_log.t,1][end]
        #         command[i] = cmd_pwm_log.z[time_vector[i].>cmd_pwm_log.t,1][end]
        #     end

        #     println("v= ",v)

        #     v_opt   = v[1:end]
        #     command_opt = command[1:end]
        #     v_opt   = v[command.>90]
        #     command_opt = command[command.>90]

        #     figure()
        #     plot(v_opt)

        #     # println("FLAG=",v_opt)

        #     a[index] = v_opt[end]/4

        #     println("a= ",a[index])

        #     index=index+1

        #     time_index = t
        end

        cmd.header.stamp = get_rostime()
        publish(pub, cmd)

        rossleep(loop_rate)
    end

     # Clean up buffers
    clean_up(gps_meas)
    clean_up(imu_meas)
    clean_up(cmd_log)
    clean_up(cmd_pwm_log)
    clean_up(pos_info_log)
    clean_up(vel_est_log)

    println("Exiting open loop control.")
    log_path = "$(homedir())/open_loop/output-record-$(run_id[1:4]).jld"
    if isfile(log_path)
        log_path = "$(homedir())/open_loop/output-record-$(run_id[1:4])-2.jld"
        warn("Warning: File already exists.")
    end
    save(log_path,"gps_meas",gps_meas,"imu_meas",imu_meas,"cmd_log",cmd_log,"cmd_pwm_log",cmd_pwm_log,"pos_info",pos_info_log,"vel_est",vel_est_log)#,"t_cmd",t_cmd,"t_no_cmd",t_no_cmd)
    println("Exiting node... Saved recorded data to $log_path.")
end

if ! isinteractive()
    main()
end

max_stars_repo_name: RoyaFiroozi/barc

max_issues_repo_name: RoyaFiroozi/barc

max_forks_repo_name: RoyaFiroozi/barc


text:
"""
    nullgraph!(net)

Remove all nodes and links from the network.

The first argument *net* is a FastNet structure that is be used in the simulation.

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,1,[])
Network of 0 nodes and 0 links

julia> randomgraph!(net)
Network of 1000 nodes and 2000 links

julia> nullgraph!(net)
Network of 0 nodes and 0 links
```
"""
function nullgraph!(net::FastNet)
    while countnodes_f(net)>0
        destroynode_f!(net,net.nid[1])
    end
    net
end

"""
    randomgraph!(net;<keyword arguments>)

Create an ER random graph in the network *net*.

The network isn't guaranteed to be a simple graph, but in large sparse
networks it is simple with high probability.

By default all nodes and links that the network can accommodate will be used and
all nodes will be set to state one. This behavior can be controlled by the following
keyword arguments:

- N : The number of nodes that will be used in the creation of the random graph.
  All other nodes will be removed from the network.
- K : The number of links that will be used in the creation of the random graph.
  All other links will be removed from the network.
- S : The state of the nodes. All nodes will be set to this state.

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,1,[])
Network of 0 nodes and 0 links

julia> randomgraph!(net)
Network of 1000 nodes and 2000 links

julia> nullgraph!(net)
Network of 0 nodes and 0 links

julia> randomgraph!(net,N=100,K=10)
Network of 100 nodes and 10 links
```
"""
function randomgraph!(net::FastNet; N::Int=0,K::Int=0,S::Int=1)
    nullgraph!(net)
    n=N;
    k=K;
    s=S;
    if n===0
        n=net.N
    end
    if k===0
        k=net.K
    end
    if n<1 && k>0
        throw(ArgumentError("In order to create links the net has to have at least one node"))
    end
    if n>net.N
        throw(ArgumentError("Trying to create more nodes than maximum allowed by the net"))
    end
    if k>net.K
        throw(ArgumentError("Trying to create more links than maximum allowed by the net"))
    end
    if s<1 || s>net.C-1
        msg="The net passed to randomgraph! only supports node states between 1 and $(net.C-1),"
        msg*=" but you are asking it to set nodes to state $s."
        throw(ArgumentError(msg))
    end
    makenodes!(net,n,s)
    for i=1:k
        src=randomnode_f(net,s)
        dst=randomnode_f(net,s)
        makelink_f!(net,src,dst)
    end
    net
end

"""
    configmodel!(net::FastNet,degreedist;<keyword arguments>)

Create a configuration model style network with prescribed degree distribution, *degreedist*.

The network in *net* is replaced with the new topology. The degree distribution is specified in
terms of a vector of Float64 variables such that *degreedist[k]* is specifies, pk, the probability
that a randomly drawn node has degree k. If the elements of *degreedist* add up to less than 1.0 the
remaining nodes will have degree zero.

The network generation is fast and unbiased but isn't guaranteed to result in a simple graph.
The algorithm will try to match the desired degree distribution as closely as possible,
but small discrepancies can appear if the degree distribution would result in an odd degree sum
or non integer numbers of nodes of certain degrees.

If there FastNet is not large enough to accomodate the desired number of links or nodes an argument error
will be thrown.

The keyword arguments are
- N : The number of nodes that will be used in the creation of the network
- S : The state of the nodes. All nodes will be set to this state.

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> configmodel!(net,[0.5,0.25,0.25],N=200)
Network of 200 nodes and 175 links

julia> degreedist(net)
3-element Vector{Float64}:
 0.5
 0.25
 0.25

julia> configmodel!(net,[0.5,0.25],N=200)
Network of 200 nodes and 100 links

julia> degreedist(net)
2-element Vector{Float64}:
 0.5
 0.25
```
"""
function configmodel!(net::FastNet,degreedist;N=0,S=1)
    (counts,totalstubs)=_configmodelsetup!(net,degreedist,N,S)
    nn=sum(counts);
    l=length(counts)
    nde=Array{Int,1}(undef,totalstubs)
    curnode=1
    curstub=1
    for i=1:l
        for j=1:counts[i]
            for s=1:i
                nde[curstub]=curnode
                curstub+=1
            end
            curnode+=1
        end
    end
    cutoff=totalstubs
    rng=net.rng
    while cutoff>1
        stub=rand(rng,1:cutoff)
        src=node_f(net,nde[stub])
        nde[stub]=nde[cutoff]
        cutoff-=1
        stub=rand(rng,1:cutoff)
        dst=node_f(net,nde[stub])
        nde[stub]=nde[cutoff]
        cutoff-=1
        makelink_f!(net,src,dst)
    end
    net
end

"""
    regulargraph!(net::FastNet,deg;<keyword arguments>)

Create a regular graph with node degree *deg*.

The network in *net* is replaced with the new topology in which all nodes have degree *deg*
and are randomly connected. If the number of nodes is not specified the function will
try to use all nodes allowed by net.

The network generation is fast and unbiased, but isn't guaranteed to result in a simple graph.

Note that finite regular graphs with odd node degree and odd number of nodes do not exist. Hence
either *deg* or the number of nodes must be even.

If there FastNet is not large enough to accomodate the desired number of links or nodes an argument error
will be thrown.

The keyword arguments are
- N : The number of nodes that will be used in the creation of the network
- S : The state of the nodes. All nodes will be set to this state.

# Examples
```jldoctest
julia> using Fastnet

julia>  net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> regulargraph!(net,4)
Network of 1000 nodes and 2000 links

julia> degreedist(net)
4-element Vector{Float64}:
 0.0
 0.0
 0.0
 1.0
```
"""
function regulargraph!(net::FastNet,deg;N=0,S=1)
    s=0
    n=0
    d=0
    try
        s=convert(Int,S)
    catch e
        throw(ArgumentError("regulargraph expects s to be an integer"))
    end
    try
        n=convert(Int,N)
    catch e
        throw(ArgumentError("regulargraph expects N to be an integer"))
    end
    try
        d=convert(Int,deg)
    catch e
        throw(ArgumentError("regulargraph expects deg to be an integer"))
    end
    if n==0
        n=net.N
    end
    checknodestate(net,s,"Trying to create regular graph")
    if n<1
        throw(ArgumentError("regulargraph expects number of nodes n to be positive"))
    end
    if d<0
        throw(ArgumentError("regulargraph expects nodedegree deg to be non-negative"))
    end
    if n>net.N
        throw(ArgumentError("Requested size for regulargraph exceeds max node count of underlyng FastNet"))
    end
    if (n*d)÷2>net.K
        throw(ArgumentError("Requested regulargraph exceeds max link count of underlyng FastNet"))
    end
    if isodd(n) && isodd(d)
        throw(ArgumentError("A regular graph of odd node degree must have an even number of nodes."))
    end
    nullgraph!(net)
    makenodes!(net,n,s)
    totalstubs=n*d
    nde=Array{Int,1}(undef,totalstubs)
    for i=1:totalstubs
        nde[i]=((i-1)÷deg)+1
    end
    cutoff=totalstubs
    rng=net.rng
    while cutoff>0
        srcs=rand(rng,1:cutoff)         # pick the source
        src=node_f(net,nde[srcs])
        nde[srcs]=nde[cutoff]
        cutoff-=1
        dsts=rand(rng,1:cutoff)         # pick the source
        dst=node_f(net,nde[dsts])
        nde[dsts]=nde[cutoff]
        cutoff-=1
        makelink_f!(net,src,dst)
    end
    net
end

"""
    rectlattice!(net::FastNet,dims, <keyword arguments>)

Create a rectangular lattice with given dimensions *dims*.

The network in *net* is replaced with the new topology, that is a lattice specified by dims.
*dims* can be a number, in this case it indicates the number of nodes to be arranged into a 1D lattice.
Alternatively, *dims* can be a vector of Ints. In this case the dimension of the lattice is identical to the
length of *dims* and each element of *dims* specifies the length of the lattice in one of these dimensions.

If there FastNet is not large enough to accomodate the desired number of nodes or links an argument error
will be thrown.

Keyword arguments are
- periodic : If this argument is true the lattice is generated with periodic boundary conditions in all dimensions.
  Alternatively a Vector of Bool of the same length as *dims* can be supplied. In this case the n'th argument of
  the vector specifies if the lattice is periodic in the n'th dimension.
- S : The state of the nodes. All nodes will be set to this state.

# Examples
```jldoctest
julia> using Fastnet

julia>  net=FastNet(2000,6000,2,[])
Network of 0 nodes and 0 links

julia> rectlattice!(net,[10,20,10],periodic=[true,false,true])
Network of 2000 nodes and 5900 links

julia> degreedist(net)
6-element Vector{Float64}:
 0.0
 0.0
 0.0
 0.0
 0.1
 0.9
```
"""
function rectlattice!(
        net         ::FastNet,
        dims        ::Union{Int,Tuple,AbstractVector};
        S           ::Integer=1,
        periodic    ::Union{AbstractArray,Bool,Tuple}=false
    )
    task="Trying to create a rectangular lattice"
    d=[dims...]
    nd=length(d)
    per=[false]
    if isa(periodic,Bool)
        per=fill(periodic,nd)
    else
        per=[periodic...]
    end
    if length(per)!=length(d)
        throw(ArgumentError(task*", but the number of values passed for parameter periodic does not agree with dims"))
    end
    for x in per
        if !isa(x,Bool)
            throw(ArgumentError(task*", but $x in periodic is not of type Bool"))
        end
    end
    n=1
    s=checknodestate(net,S,task)
    for i=1:nd
        n*=d[i]
    end
    nullgraph!(net)
    makenodes!(net,n,s)
    linksneeded=0
    for i=1:nd
        linksneeded+=n
        if !per[i]
            linksneeded-=n÷d[i]
        end
    end
    for i=1:n
        lowmult=1
        for j=1:nd
            if ((i-1)÷lowmult)%d[j]==d[j]-1
                if per[j]
                    src=node_f(net,i)
                    dst=node_f(net,i+lowmult-d[j]*lowmult)
                    makelink_f!(net,src,dst)
                end
            else
                src=node_f(net,i)
                dst=node_f(net,i+lowmult)
                makelink_f!(net,src,dst)
            end
            lowmult*=d[j]
        end
    end
    net
end

"""
    adjacency!(net,mat;S=1)

Create a network with given adjacency matrix.

The network in *net* is replaced with the new topology that is specified by the adjacency matrix *mat*.
If direction of links matters note that the element *mat[i,j]* corresponds to the link from j to i.

Symmeric matrices will not result in parallel links, instead the link is placed in an arbitrary direction.

Note that node *n* in the matrix will be the node in position *n* in *net* after creation, which is
not necessarily the node with ID *n*, if you need to find a particular node at a later time then it
is best to save its id using the node(net,pos) function directly after calling adjacency!(net,mat).

If *net* is not large enough to accomodate the desired number of nodes or links an argument error
will be thrown.

Keyword arguments are
- S : The state of the nodes. All nodes will be set to this state.

# Examples
```jldoctest
julia> using Fastnet

julia> net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> mat=[0 1 0; 1 0 1; 0 1 0]
3×3 Matrix{Int64}:
 0  1  0
 1  0  1
 0  1  0

julia> adjacency!(net,mat)
Network of 3 nodes and 2 links
```
"""
function adjacency!(
        net     ::FastNet,
        mat     ::AbstractMatrix;
        S       ::Integer=1
    )
    task="Trying to create topology from adjacency matrix"
    x,y=size(mat)
    if x!=y
        err=task*", but adjacency! expects a square matrix as its second argument and received a rectangular one"
        throw(ArgumentError(err))
    end
    if x>net.N
        throw(ArgumentError(task*", but the matrix is larger than the maximum number of allowed nodes by the network"))
    end
    s=checknodestate(net,S,task)
    count=0
    b1=false
    b2=false
    for i=1:x-1
        for j=i+1:x
            try
                b1=convert(Bool,mat[i,j])
                b2=convert(Bool,mat[i,j])
            catch e
                throw(ArgumentError(task*", but was unable to convert an element of mat to Bool"))
            end
            if b1||b2
                count+=1
            end
        end
    end
    if count>net.K
        throw(ArgumentError(task*", but the matrix contains more links than permitted by net"))
    end
    nullgraph!(net)
    makenodes!(net,x,s)
    for i=1:x
        dst=node_f(net,i)
        for j=i+1:x
            src=node_f(net,j)
            linked::Bool = mat[i,j]
            if linked
                makelink_f!(net,src,dst)
            else
                linked = mat[j,i]
                if linked
                    makelink_f!(net,dst,src)
                end
            end
        end
    end
    net
end

### WIP
function configmodel_DG!(net::FastNet,degreedist,N::Int=0,S::Int=1)
    (counts,totalstubs)=_configmodelsetup!(net,degreedist,N,S)

    println("On track")
end

"""
    randomgeometricgraph!(net::FastNet; N::Int,K::Int,S::Int, dim::Int, deg::Float64)

Create a random geometric graph with mean degree *deg*.

The network in *net* is replaced with the new topology in which each node has *dim* vector entries uniformly drawn
from (0,1). A maximal connection distance is calculated from the given mean degree and nodes get connected,
when their euclidean distance is smaller than the maximal connection distance.

If there FastNet is not large enough to accomodate the desired number of links or nodes an argument error
will be thrown.

Additionaly, if the calculated maximal connection distance results in more links being generated than supported by the net,
an argument error will be thrown as well.

The keyword arguments are
- N : The number of nodes that will be used in the creation of the network
- S : The state of the nodes. All nodes will be set to this state.
- dim : Dimensions of the random geometric graph.
- deg : The graph´s mean degree which is used to calculate the maximal connection distance between nodes.

# Examples
```jldoctest
julia> using Fastnet

julia>  net=FastNet(1000,2000,2,[])
Network of 0 nodes and 0 links

julia> randomgeometricgraph!(net,dim=2,meandegree=2.0)
Network of 1000 nodes and 970 links

```
"""
function randomgeometricgraph!(net::FastNet; N::Int=0,K::Int=0,S::Int=1, dim::Int=2, deg::Float64=2.0)
    nullgraph!(net)
    n=N;
    k=K;
    s=S;
    if n===0
        n=net.N
    end
    if k===0
        k=net.K
    end
    if n<1 && k>0
        throw(ArgumentError("In order to create links the net has to have at least one node"))
    end
    if n>net.N
        throw(ArgumentError("Trying to create more nodes than maximum allowed by the net"))
    end
    if k>net.K
        throw(ArgumentError("Trying to create more links than maximum allowed by the net"))
    end
    if s<1 || s>net.C-1
        msg="The net passed to randomgeometricgraph! only supports node states between 1 and $(net.C-1),"
        msg*=" but you are asking it to set nodes to state $s."
        throw(ArgumentError(msg))
    end

    r = (1/2)*((deg/n) * double_factorial(dim) / (pi/2)^(floor(dim/2)))^(1/dim)
    makenodes!(net,n,s)
    loc = rand(n,dim)
    for i = 1:n
        for j = i+1:n
            dist = sqrt(sum((loc[i,:].-loc[j,:]).^2))
            if (dist < r)
                if countlinks_f(net)>=k
                    throw(ArgumentError("For the given mean degree, you would create more links than allowed by the net"))
                end
                makelink_f!(net,i,j)
            end
        end
    end
    net
end

function double_factorial(n::Int)
    if isodd(n)
        k=Int((n+1)/2)
        return factorial(2*k-1)/(2^(k-1)*factorial(k-1))
    else
        k=Int(n/2)
        return 2^k*factorial(k)
    end
end


text:
#= Reference data from
- Fischer, C. (1977). The Hartree-Fock Method for Atoms : A Numerical
  Approach. New York: Wiley.
=#

# # Table 6.2
# exact_Yᵏs = Dict(((1,0),(1,0),0) => r -> 1 - exp(-2r)*(1 + r),
#                  ((2,0),(2,0),0) => r -> 1 - exp(-r)*(1 + 3/4*r + 1/4*r^2 + 1/8*r^3),
#                  ((2,1),(2,1),0) => r -> 1 - exp(-r)*(1 + 3/4*r + 1/4*r^2 + 1/24*r^3),
#                  ((2,1),(2,1),2) => r -> 30/r^2*(1 - exp(-r)*(1 + r + 1/2*r^2 + 1/6*r^3 +
#                                                             1/24*r^4 + 1/144*r^5)))

# Reference data generated using analytic-integrals.jl
include("exact_Yks.jl")

# Table 6.3–6.3B
exact_Fᵏs = Dict(((1,0),(1,0),0)=>5//8,
                 ((2,0),(1,0),0)=>17//81,
                 ((2,0),(2,0),0)=>77//512,
                 ((2,1),(1,0),0)=>59//243,
                 ((2,1),(2,0),0)=>83//512,
                 ((2,1),(2,1),0)=>93//512,
                 ((4,0),(4,0),0)=>19541//524288,
                 ((4,0),(4,1),0)=>19943//524288,
                 ((4,0),(4,2),0)=>20693//524288,
                 ((4,0),(4,3),0)=>21743//524288,
                 ((4,1),(4,1),0)=>20413//524288,
                 ((4,1),(4,2),0)=>21239//524288,
                 ((4,1),(4,3),0)=>22373//524288,
                 ((4,2),(4,2),0)=>22373//524288,
                 ((4,2),(4,3),0)=>23759//524288,
                 ((4,3),(4,3),0)=>26333//524288)

# Table 6.3
exact_Gᵏs = Dict(((2,0),(1,0),0)=>16//729,
                 ((2,1),(1,0),1)=>112//2187,
                 ((2,1),(2,0),1)=>45//512)

#= Table 1 of

- Golden, L. (1978). Exact Slater Integrals. Computer Physics
  Communications, 14(3-4),
  255–260. http://dx.doi.org/10.1016/0010-4655(78)90018-8
=#

exact_Fᵏs[(((1,0), (1,0), 0))] = (5//8)
exact_Gᵏs[(((1,0), (1,0), 0))] = (5//8)
exact_Fᵏs[(((1,0),(2,0), 0))] = (17//81)
exact_Gᵏs[(((1,0),(2,0), 0))] = (16//729)
exact_Fᵏs[(((1,0),(3,0), 0))] = (815//8192)
exact_Gᵏs[(((1,0), (3,0), 0))] = (189//32768)
exact_Fᵏs[(((2,0),(2,0), 0))] = (77//512)
exact_Gᵏs[(((2,0), (2,0), 0))] = (77//512)
exact_Fᵏs[(((2,0), (3,0), 0))] = (197142//2343750)
exact_Gᵏs[(((2,0), (3,0), 0))] = (73008//9765625)
exact_Fᵏs[(((3,0), (3,0), 0))] = (17//256)
exact_Gᵏs[(((3,0), (3,0), 0))] = (17//256)
exact_Fᵏs[(((1,0),(2,1), 0))] = (59//243)
exact_Gᵏs[(((1,0), (2,1), 1))] = (112//2187)
exact_Fᵏs[(((1,0),(3,1), 0))] = (1783//16384)
exact_Gᵏs[(((1,0), (3,1), 1))] = (891//65536)
exact_Fᵏs[(((2,0), (2,1), 0))] = (83//512)
exact_Gᵏs[(((2,0), (2,1), 1))] = (45//512)
exact_Fᵏs[(((2,0), (3,1), 0))] = (35161//390625)
exact_Gᵏs[(((2,0), (3,1), 1))] = (82944//9765625)
exact_Fᵏs[(((3,0), (3,1), 0))] = (317//4608)
exact_Gᵏs[(((3,0), (3,1), 1))] = (65//1536)
exact_Fᵏs[(((3,0), (2,1), 0))] = (33889//390625)
exact_Gᵏs[(((3,0), (2,1), 1))] = (92016//9765625)
exact_Fᵏs[(((1,0), (3,2), 0))] = (1819//16384)
exact_Gᵏs[(((1,0), (3,2), 2))] = (81//65536)
exact_Fᵏs[(((2,0), (3,2), 0))] = (202301//1953125)
exact_Gᵏs[(((2,0), (3,2), 2))] = (1783296//48828125)
exact_Fᵏs[(((3,0),(3,2), 0))] = (337//4608)
exact_Gᵏs[(((3,0), (3,2), 2))] = (35//1536)
exact_Fᵏs[(((2,1),(2,1), 0))] = (93//512)
exact_Gᵏs[(((2,1),(2,1), 0))] = (93//512)
exact_Fᵏs[(((2,1), (3,1), 0))] = (36697//390625)
exact_Gᵏs[(((2,1), (3,1), 0))] = (96768//9765625)
exact_Fᵏs[(((3,1), (3,1), 0))] = (1987//27648)
exact_Gᵏs[(((3,1), (3,1), 0))] = (1987//27648)
exact_Fᵏs[(((2,1), (2,1), 2))] = (45//512)
exact_Gᵏs[(((2,1), (2,1), 2))] = (45//512)
exact_Fᵏs[(((2,1), (3,1), 2))] = (2022//78125)
exact_Gᵏs[(((2,1), (3,1), 2))] = (110592//9765625)
exact_Fᵏs[(((3,1), (3,1), 2))] = (995//27648)
exact_Gᵏs[(((3,1), (3,1), 2))] = (995//27648)
exact_Fᵏs[(((2,1), (3,2), 0))] = (207677//1953125)
exact_Gᵏs[(((2,1), (3,2), 1))] = (1824768//48828125)
exact_Fᵏs[(((3,1),(3,2), 0))] = (709//9216)
exact_Gᵏs[(((3,1), (3,2), 1))] = (105//3072)
exact_Fᵏs[(((2,1), (3,2), 2))] = (14502//390625)
exact_Gᵏs[(((2,1), (3,2), 3))] = (1064448//48828125)
exact_Fᵏs[(((3,1), (3,2), 2))] = (37//1024)
exact_Gᵏs[(((3,1), (3,2), 3))] = (1995//82944)
exact_Fᵏs[(((3,2), (3,2), 0))] = (793//9216)
exact_Gᵏs[(((3,2), (3,2), 0))] = (793//9216)
exact_Fᵏs[(((3,2), (3,2), 2))] = (2093//46080)
exact_Gᵏs[(((3,2), (3,2), 2))] = (2093//46080)
exact_Fᵏs[(((3,2), (3,2), 4))] = (91//3072)
exact_Gᵏs[(((3,2), (3,2), 4))] = (91//3072)


text:
using AxisArrays
using JuMP
using MathOptInterface


text:
struct Mystruct
    data::Ptr{Cdouble}
    size::Csize_t
    Mystruct() = new(Ptr{Cdouble}(), 0)
end

function print_struct(v::Ref{Mystruct})
    ccall((:recvstruct, "./libfoo.so"), Cvoid, (Ptr{Mystruct},), v)
end

function update_struct()
    st = Ref(Mystruct())

    # Use Ptr but not Ref since the C is managing its memory and
    # we do not want Julia manages the memory (garbage collecting
    # and producing a segfault)
    ccall((:updatestruct, "./libfoo.so"), Cvoid, (Ptr{Mystruct},), st)

    # Display the C array
    print_struct(st)

    # Convert the C array into a Julia vector
    V = Vector{Cdouble}(undef, st[].size)
    for i in 1:st[].size
        V[i] = unsafe_load(st[].data, i)
    end

    return V
end

V = update_struct()
show(V)

max_stars_repo_name: Lecrapouille/BacASable

max_issues_repo_name: Lecrapouille/BacASable

max_forks_repo_name: Lecrapouille/BacASable


text:
using ReinforcementLearning, Base.Test

env = GymEnv("Breakout-v0")

@test length(env.actions) == 4

start = getInitialState(env)
@test isTerminal(start, env) == false
@test size(start.data) == (210, 160, 3)
@test (sample(env) in env.actions) == true

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
import Base: sum, findmax, findmin, maximum, argmax
import LinearAlgebra: norm
export average

function sum(x::CuArray; dims=())
    if isempty(dims)
        x = vec(x)
        dims = (1,)
    else
        isa(dims,Int) && (dims = (dims,))
    end
    for d in dims
        x = CUDNN.reduce(x, d, CUDNN.CUDNN_REDUCE_TENSOR_ADD)[1]
    end
    x
end

function argmax(x::CuArray, dims::Int)
    findmax(x, dims=dims)[2]
end

function findmax(x::CuArray; dims::Int)
    dim = dims
    y, idx = CUDNN.reduce(x, dim, CUDNN.CUDNN_REDUCE_TENSOR_MAX)
    idx = reshape(idx, Base.setindex(size(x),1,dim))
    y, idx
end

function findmin(x::CuArray; dims::Int)
    CUDNN.reduce(x, dims, CUDNN.CUDNN_REDUCE_TENSOR_MIN)
end

average(x::CuArray; dims::Int) = CUDNN.reduce(x, dims, CUDNN.CUDNN_REDUCE_TENSOR_AVG)[1]

function norm(x::CuArray, dim::Int, p::Int)
    if p == 1
        CUDNN.reduce(x, dim, CUDNN.CUDNN_REDUCE_TENSOR_NORM1)[1]
    elseif p == 2
        CUDNN.reduce(x, dim, CUDNN.CUDNN_REDUCE_TENSOR_NORM2)[1]
    else
        throw("Not supported. Valid p: 1 or 2.")
    end
end

#=
function argmax(x::Array, dim::Int)
    _, index = findmax(x, dim)
    y = ind2sub(size(x), vec(index))[dim]
    dims = ntuple(i -> i==dim ? 1 : size(x,i), ndims(x))
    reshape(y, dims)
end
=#


text:
export _3PointBendingNotched, _3PointBending, _4PointBendingNotched, _4PointBending

function _3PointBendingNotched(args...; kwargs...)
    _3PointBending(args...; notched=true, kwargs...)
end

function _4PointBendingNotched(args...; kwargs...)
    _4PointBending(args...; notched=true, kwargs...)
end

function _4PointBending(args...; kwargs...)
    _3PointBending(args...; _4point=true, kwargs...)
end

function _3PointBending(; gen_mat=nothing, spc_mat=nothing, resolution=0.1, solver_=:qs, notched=false, _4point=true, file_prefix="TensileBar")
    if solver_ in [:qs]
        Steps = 100
        sargs = (Steps, 0.01)
        skwargs = Dict(:max_iter=>500)
        fwf = 1
    else
        Steps = 1000
        sargs = (Steps, )
        skwargs = Dict()
        fwf = 10
    end
    solver = PDBenchmark.NameParam(solver_, sargs, Dict(:filewrite_freq=>fwf, :file_prefix=>file_prefix, :start_at=>0, skwargs...))

    geom = []
    reso = resolution

    obj = () -> move(Disk(0.5, 1.0), by=[1.0,  -0.5, 0.0])
    push!(geom, PDBenchmark.NameParam(obj, (), Dict(:resolution=>reso, :type=>2)))

    obj = () -> move(Disk(0.5, 1.0), by=[11.0, -0.5, 0.0])
    push!(geom, PDBenchmark.NameParam(obj, (), Dict(:resolution=>reso, :type=>2)))

    obj = () -> move(Disk(0.5, 1.0), by=[6.00,  1.5, 0.0])
    push!(geom, PDBenchmark.NameParam(obj, (), Dict(:resolution=>reso, :type=>3)))

    if notched
        obj = () -> NotchedBar(bounds=[0 12; 0 1; 0 1], notch=[5.75 6.25; 0.0 0.25; 0.0 1.0])
    else
        obj = () -> Bar(bounds=[0 12; 0 1; 0 1])
    end

    push!(geom, PDBenchmark.NameParam(obj, (), Dict(:resolution=>reso, :type=>1)))

    gen_mat = [gen_mat, gen_mat, gen_mat, gen_mat]

    mat = getproperty_(PeriDyn, spc_mat.name)(spc_mat.args...; spc_mat.kwargs...)
    spc_mat0 = PDBenchmark.NameParam(:SkipSpecific, (mat.density, ), Dict())
    spc_mat = [spc_mat0, spc_mat0, spc_mat0, spc_mat]

    if solver_ in [:qs]
        bc_f1 = env -> (vec(env.type .== 2), )
        bc_f2 = env -> (vec(env.type .== 3), [0.0, -0.01, 0])
        bc = [  PDBenchmark.NameParam(:FixBC, (bc_f1,)),
                PDBenchmark.NameParam(:MoveBC, (bc_f2,)),
                ]
    else
        bc_f1 = env -> (vec(env.type .== 2), )
        bc_f2 = env -> (vec(env.type .== 3), [0.0, -0.001, 0])
        bc = [  PDBenchmark.NameParam(:FixBC, (bc_f1,)),
                PDBenchmark.NameParam(:MoveBC, (bc_f2,))
                ]
    end

    RM = []
    push!(RM, PDBenchmark.NameParam(:LinearRepulsionModel, (100.0,), Dict(:blocks=>(4,), :distanceX=>3, :max_neighs=>200)))
    push!(RM, PDBenchmark.NameParam(:LinearRepulsionModel, (100.0,), Dict(:blocks=>(4,1), :distanceX=>3, :max_neighs=>200)))
    push!(RM, PDBenchmark.NameParam(:LinearRepulsionModel, (100.0,), Dict(:blocks=>(4,2), :distanceX=>3, :max_neighs=>200)))
    push!(RM, PDBenchmark.NameParam(:LinearRepulsionModel, (100.0,), Dict(:blocks=>(4,3), :distanceX=>3, :max_neighs=>200)))

    f = env -> begin
        pin_ = (env.type .== 3)
        pin_max = maximum(env.y[2, pin_])
        env.Params = Dict("pin" => pin_, "pin_max"=>pin_max)
        env.Out = Dict("Force" => zeros(3, Steps), "Displacement" => zeros(1, Steps))
        env.Collect! = (env, step) -> begin
                env.Out["Force"][:, step] = sum(env.f[:,env.Params["pin"]], dims=2)
                env.Out["Displacement"][1, step] = -(maximum(env.y[1,env.Params["pin"]]) - env.Params["pin_max"])
                if step%fwf==0
                    println("∂P/∂y = 48EI/L^3 = 4.0e-3 * E = $(env.Out["Force"][2, step]/env.Out["Displacement"][1, step])")
                end
            end
    end

    test = PDBenchmark.Test(solver, geom, gen_mat, spc_mat, bc, RM, f; names=["LSupport", "RSupport", "Pin", "Beam"])

    return test
end

max_stars_repo_name: ravinderbhattoo/PDBenchmark

max_issues_repo_name: ravinderbhattoo/PDBenchmark

max_forks_repo_name: ravinderbhattoo/PDBenchmark