text:
"""
    IRISRequest

An abstract type representing requests conforming to the IRIS Web Services specification.

Current subtypes of `IRISRequest`:
- `IRISTimeSeries`
"""
abstract type IRISRequest <: SeisRequest end

protocol_string(::IRISRequest) = "irisws"

"Available processing filters and allowable types for them"
const IRISTimeSeries_PROCESSING_FIELDS = Dict(
    :taper => Real,
    :taper_type => AbstractString,
    :envelope => Bool,
    :lpfilter => Real,
    :hpfilter => Real,
    :bpfilter => Tuple{Number,Number},
    :demean => Bool,
    :scale => Union{Real,AbstractString},
    :divscale => Real,
    :correct => Bool,
    :freqlimits => Tuple{Number,Number,Number,Number},
    :autolimits => Union{AbstractString,Tuple{Number,Number}},
    :units => AbstractString,
    :diff => Bool,
    :int => Bool,
    :decimate => Real)

const IRISTimeSeries_OUTPUT_TYPES = ("ascii1", "ascii2", "ascii", "geocsv", "geocsv.tspair",
    "geocsv.slist", "audio", "miniseed", "plot", "saca", "sacbb", "sacbl")

const IRISTimeSeries_FORMAT_TYPES = ("ascii1", "ascii2", "ascii", "geocsv", "geocsv.tspair",
    "geocsv.slist", "audio", "miniseed", "plot", "sac.zip")

"""
    IRISTimeSeries

Create an data query which can be sent to a datacentre which implements the IRIS
Web Services timeseries specification.

## Available options

### Channel Options (required)
The four SCNL parameters (Station – Channel – Network – Location) are used to determine
the channel of interest, and are all required. Wildcards are not accepted.

- `network`: Seismic network name
- `station`: Station name
- `location`: Location code. Use `loc=--` for empty location codes
- `channel`: Channel Code

This can instead be combined into the single option:

- `code`: SEED-style channel code of the form `"⟨network⟩.⟨station⟩.⟨location⟩.⟨channel⟩"`.

Note that `code` cannot be combined with the separate SCNL keywords.

### Date-Range Options (required)
Time ranges must be defined by specifying a start time with either an end time or a duration
in seconds.

For example, the following two time specifications are valid, and are equivalent.

#### Examples

    IRISTimeSeries([SCNL]..., starttime=Dates.DateTime(2005, 02, 23, 12), endtime=Dates.DateTime(2005, 02, 23, 12, 05))
    IRISTimeSeries([SCNL]..., starttime=Dates.DateTime("2005-02-23T12:00:00"), duration=300)

- `starttime`: Start time
- `endtime`: End time
- `duration`: Duration of requested data, in seconds

### Time Series Processing Options
The following parameters act as filters upon the timeseries. Parameter order matters
because each filter operation is performed in the order given.

#### Examples

    IRISTimeSeries(..., demean=true, lpfilter=2.0) # will demean and then apply a low-pass filter
    IRISTimeSeries(..., lpfilter=2.0, demean=true) # will apply a low-pass filter, and then demean

- `taper`: Apply a time domain symmetric tapering function to the timeseries data. The
  width is specified as a fraction of the trace length from 0 to 0.5.
- `taper_type`: Form of taper to apply.  Supported types: "HANNING" (default), "HAMMING", "COSINE".
  Note that all input is uppercased automatically.
- `envelope`: Calculate the envelope of the time series. This calculation uses a Hilbert
  transform approximated by a time domain filter
- `lpfilter`: Low-pass filter the time-series using an IIR 4th order filter, using this value
  (in Hertz) as the cutoff
- `hpfilter`: High-pass filter the time-series using an IIR 4th order filter, using this value
  (in Hertz) as the cutoff
- `bpfilter`: Band pass frequencies, in Hz.  This may be given as a length-2 tuple or AbstractArray,
  or as an AbstractString. (If the latter, then two values must be separated by either `-`, `/`, `,` or `;`.)
- `demean`: Remove mean value from data. `true` or `false` (default `false`).
- `scale`: Scale data samples by specified factor. When `scale="AUTO"` scales by the stage-zero gain.
  Mutually exclusive with `divscale`.
- `divscale`: Scale data samples by the inverse of the specified factor.  Mutually exclusive with
  `scale`.
- `correct`: Apply instrument correction to convert to earth units. `true` or `false`. Uses
  either deconvolution3,4 or polynomial response correction. (default `false`).
  Mutually exclusive with `scale="AUTO"`.
- `freqlimits`: Specify an envelope for a spectrum taper for deconvolution. Frequencies are
  specified in Hertz. This cosine taper scales the spectrum from 0 to 1 between f1 and f2 and
  from 1 to 0 between f3 and f4. Can only be used with the correct option. Cannot be used in
  combination with the autolimits option.
- `autolimits`: Automatically determine frequency limits for deconvolution. A pass band is
  determined for all frequencies with the lower and upper corner cutoffs defined in terms of
  dB down from the maximum amplitude. This algorithm is designed to work with flat responses,
  i.e. a response in velocity for an instrument which is flat to velocity. Other combinations
  will likely result in unsatisfactory results. Cannot be used in combination with the
  `freqlimits` option.
- `units`: Specify output units. Can be "DIS", "VEL", "ACC" or "DEF", where "DEF" results in no unit
  conversion.  `units` can only be used with `correct`.
- `diff`: Differentiate using 2 point (uncentered) method. `true` or `false` (default `false`)
- `int`: Integrate using trapezoidal (midpoint) method. `true` or `false` (default `false`)
- `decimate`: Sample-rate to decimate to. See [online help](http://service.iris.edu/irisws/timeseries/docs/1/help/#deci)
  for more details. A linear-phase, anti-alias filter is applied during decimation.

### Plot Options
- `antialias`: If true, the created image will have anti-aliasing used. Can only be specified
  with the `output=plot` option
- `width`: Width of the output plot (pixels). Can only be specified with the `output=plot` option
- `height`: Height of the output plot (pixels). Can only be specified with the `output=plot` option

### Audio Options
- `audiosamplerate`: The sample rate of the output wav file in Hz. Defaults to 16000.
  Can only be specified with the output=audio option
- `audiocompress`: Apply dynamic range compression to waveform data. This makes more of the
  signal audible. Can only be specified with the output=audio option

## Format options (required)

Provide one of the following to determine the output format with the `format` option:

|Value|Description|
|:----|:----------|
|`ascii1`|ASCII data format, 1 column (values)|
|`ascii2`|ASCII data format, 2 columns (time, value)|
|`ascii`|Same as `ascii2`|
|`geocsv`, `geocsv.tspair`|ASCII GeoCSV data format, 2 columns (time, values)|
|`geocsv.slist`|ASCII geocsv data format, 1 column (value)|
|`audio`|Audio WAV file|
|`miniseed`|FDSN miniSEED format [default]|
|`plot`|A simple plot of the timeseries in PNG format|
|`sac.zip`|SAC Zipped binary traces|

The following are deprecated and can only be used with the `output` option, which cannot
be combined with the `format` option:

|Value|Description|
|:----|:----------|
|`saca`|SAC – ASCII format [deprecated and used with `output`]|
|`sacbb`|SAC – binary big-endian format [deprecated and used with `output`]|
|`sacbl`|SAC – binary little-endian format [deprecated and used with `output`]|

As these outputs are deprecated, they may stop working in the future.

## Specifying dates
Both `Dates.DateTime` objects and `String`s can be passed as values to the
`starttime` and `endtime` keywords so long as they are valid.
Examples include `"2000-01-01"` and `"2008-02-03T00:00:30"`.
"""
struct IRISTimeSeries <: IRISRequest
    network::String
    station::String
    location::String
    channel::String
    starttime::DateTime
    endtime::MDateTime
    duration::MFloat
    antialias::MBool
    width::MInt
    height::MInt
    audiosamplerate::MInt
    audiocompress::MBool
    output::MString
    format::MString
    process::OrderedDict{Symbol,Any} # Because order matters for processing

    function IRISTimeSeries(network, station, location, channel, starttime, endtime, duration,
                            antialias, width, height, audiosamplerate, audiocompress, output,
                            format, process)
        location in ("  ", "") && (location = "--")
        !ismissing(starttime) && (starttime = DateTime(starttime))
        !ismissing(endtime) && (endtime = DateTime(endtime))
        !ismissing(starttime) && !ismissing(endtime) && starttime > endtime &&
            throw(ArgumentError("`starttime` must be before endtime"))
        coalesce(duration, 1) > 0 || throw(ArgumentError("`duration` must be positive"))
        0 <= get(process, :taper, 0.1) <= 0.5 ||
            throw(ArgumentError("`taper` length must be between 0 and 0.5"))
        # Check types for processing
        for (f, T) in IRISTimeSeries_PROCESSING_FIELDS
            if haskey(process, f)
                try
                    if T <: Tuple
                        # Convert arrays
                        process[f] = convert(T, tuple(process[f]...))
                    else
                        process[f] = convert(T, process[f])
                    end
                catch err
                    throw(ArgumentError("field `$f` must be of type $T"))
                    rethrow(err)
                end
            end
        end
        if haskey(process, :taper_type)
            haskey(process, :taper) || throw(ArgumentError("`taper_type` must be specficied with taper"))
            taper_type = uppercase(process[:taper_type])
            taper_type in ("HANNING", "HAMMING", "COSINE") ||
            throw(ArgumentError("taper_type must be one of \"HANNING\", \"HAMMING\" or \"COSINE\""))
        end
        haskey(process, :lpfilter) && process[:lpfilter] <= 0 &&
            throw(ArgumentError("Low-pass filter cutoff must be greater than 0 Hz"))
        haskey(process, :hpfilter) && process[:hpfilter] <= 0 &&
            throw(ArgumentError("High-pass filter cutoff must be greater than 0 Hz"))
        if haskey(process, :bpfilter)
            0 < process[:bpfilter][1] < process[:bpfilter][2] ||
                throw(ArgumentError("Low-pass corner of band-pass filter must be above high-pass corner"))
        end
        if haskey(process, :scale)
            haskey(process, :divscale) &&
                        throw(ArgumentError("Cannot specify both `scale` and `divscale`"))
            if process[:scale] isa AbstractString
                process[:scale] = uppercase(process[:scale])
                process[:scale] == "AUTO" || throw(ArgumentError("`scale` must be a float or \"AUTO\""))
                get(process, :correct, false) &&
                    throw(ArgumentError("Cannot specify both `correct` and `scale=\"AUTO\"`"))
            end
        end
        if haskey(process, :freqlimits)
            get(process, :correct, false) ||
                throw(ArgumentError("`freqlimits` can only be used when `correct`=`true`"))
            f1, f2, f3, f4 = process[:freqlimits]
            0 < f1 < f2 < f3 < f4 ||
                throw(ArgumentError("`freqlimits` must increase monotonically and all be positive"))
        end
        if haskey(process, :units)
            get(process, :correct, false) ||
                throw(ArgumentError("`units` can only be used when `correct` is `true`"))
            process[:units] in ("DIS", "VEL", "ACC", "DEF") ||
                throw(ArgumentError("`units` must be one of \"DIS\", \"VEL\", \"ACC\" or \"DEF\""))
        end
        haskey(process, :autolimits) && haskey(process, :freqlimits) &&
            throw(ArgumentError("Cannot specify both `autolimits` and `freqlimits`"))

        # Disallowed combinations
        coalesce(audiocompress, false) || !ismissing(audiosamplerate) && output != "audio" &&
            throw(ArgumentError("Cannot specify `audiocompress` if `output` is not \"audio\""))
        count(ismissing, (output, format)) == 1 ||
            throw(ArgumentError("One and only one of `output` and `format` must be specified"))

        # Disallowed values
        (!ismissing(output) && output ∉ IRISTimeSeries_OUTPUT_TYPES) &&
            throw(ArgumentError("`output` ($output) must be one of $(IRISTimeSeries_OUTPUT_TYPES)"))
        (!ismissing(format) && format ∉ IRISTimeSeries_FORMAT_TYPES) &&
            throw(ArgumentError("`format` ($output) must be one of $(IRISTimeSeries_FORMAT_TYPES)"))

        new(network, station, location, channel, starttime, endtime, duration,
            antialias, width, height, audiosamplerate, audiocompress, output, format, process)
    end
end

function IRISTimeSeries(;
                        code=missing,
                        network=missing, station=missing, location=missing, channel=missing,
                        starttime=missing, endtime=missing, duration=missing, antialias=missing,
                        width=missing, height=missing, audiosamplerate=missing,
                        audiocompress=missing, output=missing, format=missing, kwargs...)
    if !ismissing(code)
        any(!ismissing, (network, station, location, channel)) &&
            throw(ArgumentError("`code` cannot be provided with any of " *
                                "`network`, `station`, `location`, or `channel`"))
        network, station, location, channel = split_channel_code(code)
    end
    any(ismissing, (network, station, location, channel, starttime)) &&
        throw(ArgumentError("network, station, location, channel and starttime, plus one of " *
                            "endtime or duration, must all be specified"))
    count(ismissing, (endtime, duration)) == 1 ||
        throw(ArgumentError("One and only one of `endtime` or `duration` must be specified"))
    if ismissing(output) && ismissing(format)
        format = "miniseed"
    end
    process = OrderedDict{Symbol,Any}()
    for (k, v) in kwargs
        k in keys(IRISTimeSeries_PROCESSING_FIELDS) ||
            throw(ArgumentError("Field `$k` is not a valid IRISTimeSeries processing field"))
        process[k] = v
    end
    IRISTimeSeries(network, station, location, channel, starttime, endtime, duration,
                   antialias, width, height, audiosamplerate, audiocompress, output,
                   format, process)
end

function Base.:(==)(r1::IRISTimeSeries, r2::IRISTimeSeries)
    for f in fieldnames(IRISTimeSeries)
        v1, v2 = getfield.((r1, r2), f)
        # Special case as OrderedDicts compare as equal even with different
        # orders of the keys!
        if f === :process
            length(v1) == length(v1) || return false
            for ((key1, val1), (key2,val2)) in zip(v1, v2)
                key1 === key2 || return false
                val1 === val2 || return false
            end
        else
            if v1 === missing || v2 === missing
                v1 === v2 || return false
            else
                v1 == v2 || return false
            end
        end
    end
    true
end

Base.broadcastable(request::IRISTimeSeries) = Ref(request)

service_string(::IRISTimeSeries) = "timeseries"

function request_uri(r::IRISTimeSeries; server=DEFAULT_SERVER)
    server = server in keys(SERVERS) ? SERVERS[server] : server
    protocol = protocol_string(r)
    service = service_string(r)
    version = version_string(r)
    uri = join((server, protocol, service, version, "query?"), "/")
    firstfield = true
    for f in fieldnames(typeof(r))
        if f == :process
            for (field, v) in r.process
                field_val = if field == :taper
                    value = "$v"
                    if haskey(r.process, :taper_type)
                        value = join((value, uppercase(r.process[:taper_type])), ",")
                    end
                    "$field=$value"
                elseif field == :taper_type
                    continue
                elseif field in (:bpfilter, :freqlimits, :autolimits)
                    "$field=$(join(v, "-"))"
                else
                    "$field=$v"
                end
                @assert !firstfield # process is the last field
                uri = join((uri, field_val), "&")
            end
        else
            v = getfield(r, f)
            if !ismissing(v)
                if firstfield
                    uri *= "$f=$v"
                    firstfield = false
                else
                    uri = join((uri, "$f=$v"), "&")
                end
            end
        end
    end
    uri
end


text:
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
# ==============================================================================
#
#   Compute the beta angle of a satellite.
#
# References
# ==============================================================================
#
# [1] Mortari, D., Wilkins, M. P., and Bruccoleri, C.  On Sun-Synchronous Orbits
#     and Associated Constellations
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #

export beta_angle

"""
    beta_angle(jd₀::Number, a::Number, e::Number, i::Number, Ω::Number, Δt::Integer; kwargs...)

Compute the beta angle of an orbit with semi-major axis `a` [m], eccentricity
`e`, inclination `i` [rad], and initial right ascension of the ascending node
`Ω₀` [rad]. The orbit epoch, which is also the day in which the analysis will
begin, is `jd₀` [Julian Day]. The analysis will be performed for each day during
`Δt` days. The algorithm was obtained from **[1]**.

!!! note
    It is expected that the input elements are represented in the TOD reference
    frame. If it is not the case, they can be converted using the function
    `change_oe_frame` of **SatelliteToolbox.jl**.

# Keywords

- `pert::Symbol`: Select the perturbation terms that must be used when
    propagating the right ascencion of the ascending node. The possible values
    are:
    - `:J0`: Consider a Keplerian orbit.
    - `:J2`: Consider the perturbation terms up to J2.
    - `:J4`: Consider the perturbation terms J2, J4, and J2².
    (**Default**: `:J2`)

# Returns

- An array with the beta angle [rad] for each day.

# References

- **[1]**: Mortari, D., Wilkins, M. P., and Bruccoleri, C.  On Sun-Synchronous
           Orbits and Associated Constellations

# Extended help

The beta angle is the angle between the orbit plane and the Sun. The positive
direction is defined as that of the orbit angular momentum.

The beta angle is useful when computing the mean amount of solar radiation a
satellite receives in a particular orbit.

## Examples

```julia
julia> using SatelliteAnalysis, UnicodePlots

julia> jd₀ = SatelliteAnalysis.date_to_jd(2021, 1, 1, 0, 0, 0)

julia> β = beta_angle(
    jd₀,
    7130.982e3,
    0.001111,
    deg2rad(98.405),
    SatelliteAnalysis.compute_RAAN_lt(jd₀, 22.5),
    365
)

julia> lineplot(0:1:364, rad2deg.(β), xlabel = "Day", ylabel = "Beta angle [°]")
                     ┌────────────────────────────────────────┐
                  30 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⣠⠞⠉⠙⢦⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⡴⠁⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⣰⠃⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠹⡄⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⡰⠃⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠹⡀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⢀⣀⡀⠀⠀⠀⠀⠀⠀⢀⠞⠁⠀⠀⠀⠀⠀│
   Beta angle [°]    │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡀⠀⠀⠀⠀⠀⠀⠀⠀⢀⠴⠋⠉⠀⠉⠙⠒⠤⠤⠤⠖⠉⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠱⡄⠀⠀⠀⠀⠀⢀⠔⠁⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠘⠦⣀⣀⣠⠔⠉⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                  10 │⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀│
                     └────────────────────────────────────────┘
                     ⠀0⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀400
                     ⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀⠀Day
```
"""
function beta_angle(
    jd₀::Number,
    a::Number,
    e::Number,
    i::Number,
    Ω₀::Number,
    Δt::Integer;
    pert::Symbol = :J2
)
    # Check the input.
    Δt ≤ 0 && error("The number of days for the analysis `Δt` must be equal or higher than 0")

    # Vector of the days in which the beta angle will be computed.
    days = 1:1:Δt

    # Output vector.
    β = Vector{Float64}(undef, Δt)

    # RAAN rotation rate [rad/day].
    δΩ = 86400 * dRAAN(a, e, i, pert)

    # Loop
    @inbounds for d in days
        # Compute the RAAN at the day d.
        Ω = Ω₀ + δΩ * d

        # Compute the unit vector aligned with the orbit normal `n` represented
        # in the TOD reference frame.
        D_tod_o = angle_to_dcm(-i, -Ω, 0, :XZX)
        n̄_tod   = D_tod_o * (@SVector [0, 0, 1])

        # Compute the Sun position at noon (UT) represented in the TOD reference
        # frame.
        s_mod     = sun_position_i(jd₀ + d)
        D_tod_mod = r_eci_to_eci(TOD(), jd₀, MOD(), jd₀)
        s̄_tod     = D_tod_mod * s_mod / norm(s_mod)

        # Compute the beta angle, which is the angle between the Sun vector and
        # the orbit plane.
        β[d] = π / 2 - acos(n̄_tod ⋅ s̄_tod)
    end

    return β
end

ext: jl


text:
using brainflow

# enable all possible logs from all three libs
brainflow.enable_dev_brainflow_logger(brainflow.BOARD_CONTROLLER)
brainflow.enable_dev_brainflow_logger(brainflow.DATA_HANDLER)
brainflow.enable_dev_brainflow_logger(brainflow.ML_MODULE)

params = BrainFlowInputParams()
board_shim = brainflow.BoardShim(brainflow.SYNTHETIC_BOARD, params)
sampling_rate = brainflow.get_sampling_rate(brainflow.SYNTHETIC_BOARD)
nfft = brainflow.get_nearest_power_of_two(sampling_rate)

brainflow.prepare_session(board_shim)
brainflow.start_stream(board_shim)
sleep(5)
brainflow.stop_stream(board_shim)
data = brainflow.get_board_data(board_shim)
brainflow.release_session(board_shim)

eeg_channels = brainflow.get_eeg_channels(brainflow.SYNTHETIC_BOARD)

bands = brainflow.get_avg_band_powers(data, eeg_channels, sampling_rate, true)
feature_vector = vcat(bands[1], bands[2])

# calc concentration
model_params = BrainFlowModelParams(metric = "concentration", classifier = "KNN")
brainflow.prepare(model_params)
print(brainflow.predict(feature_vector, model_params))
brainflow.release(model_params)

# calc relaxation
model_params = BrainFlowModelParams(metric = "relaxation", classifier = "regression")
brainflow.prepare(model_params)
print(brainflow.predict(feature_vector, model_params))
brainflow.release(model_params)

max_stars_repo_name: matthijscox/brainflow

max_issues_repo_name: matthijscox/brainflow

max_forks_repo_name: matthijscox/brainflow


text:
# updated version of nested_optim_d_only.jl
# run align_prep.jl first.

include("../src/NMRCalibrate.jl")
import .NMRCalibrate

PyPlot.close("all")
fig_num = 1

Random.seed!(25)
PyPlot.matplotlib["rcParams"][:update](["font.size" => 22, "font.family" => "serif"])

### optim all regions. # 900 secs.
region_select = 3
y_cost = y[cost_inds_set[region_select]]
U_cost = U_y[cost_inds_set[region_select]]
P_cost = P_y[cost_inds_set[region_select]]

# y_cost = y_cost_all
# U_cost = U_cost_all
# P_cost = P_cost_all

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P_y, real.(y), label = "data spectrum")
PyPlot.plot(P_cost, real.(y_cost), "^", label = "positions")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("positions against data spectrum, r = $(region_select), real part")

### optim.
N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )
N_β = sum( NMRCalibrate.getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
LS_inds = 1:length(U_cost)

# inputs.
κ_lb_default = 1e-2
κ_ub_default = 1e2
p_initial = zeros(N_d)
p_β = zeros(N_β)

# prepare.
U_rad_cost = U_cost .* (2*π)
Es = collect( NMRSpectraSimulator.καFIDModelType(Bs[i]) for i = 1:length(Bs) )

a_setp, b_setp, minxs,
    rets = NMRCalibrate.setupitpab(0.1, 10, 0.7; optim_algorithm = :LN_BOBYQA)
#
β_optim_algorithm = :GN_DIRECT_L
β_max_iters = 500
β_xtol_rel = 1e-9
β_ftol_rel = 1e-9
β_maxtime = Inf
q, updatedfunc, getshiftfunc, N_vars_set,
run_optim, obj_func_β, E_BLS, κ_BLS, b_BLS, updateβfunc, updateκfunc,
q_β = NMRCalibrate.setupcostnesteddwarp(Es, Bs, As, fs, SW, LS_inds, U_rad_cost,
    y_cost, Δsys_cs, a_setp, b_setp, κs_β_DOFs, κs_β_orderings;
    w = w,
    optim_algorithm = β_optim_algorithm,
    κ_lb_default = κ_lb_default,
    κ_ub_default = κ_ub_default,
    max_iters = β_max_iters,
    xtol_rel = β_xtol_rel,
    ftol_rel = β_ftol_rel,
    maxtime = β_maxtime)
#
obj_func = pp->NMRCalibrate.costnestedd(U_rad_cost, y_cost, updatedfunc, updateκfunc, pp,
Es, Bs, κs_β_orderings, κs_β_DOFs, q, run_optim, E_BLS, κ_BLS, b_BLS, p_β)

grad_func = xx->FiniteDiff.finite_difference_gradient(obj_func, xx)

#NMRCalibrate.evaldesignmatrixκ!(E_BLS, U_rad_cost, Es, As, w)

#optim_algorithm = :LN_BOBYQA
optim_algorithm = :GN_ESCH
#optim_algorithm = :GN_ISRES
#optim_algorithm = :GN_DIRECT_L # larger objective!
opt = NLopt.Opt(optim_algorithm, length(p_initial))

println("obj_func(p_initial) = ", obj_func(p_initial))
println()

# for display.
U = LinRange(u_min, u_max, 50000)
P = hz2ppmfunc.(U)
U_rad = U .* (2*π)

q_U = q.(U_rad)

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P_y, real.(y), label = "y")
PyPlot.plot(P, real.(q_U), label = "q")
#PyPlot.plot(P, real.(q3_star_U), "--", label = "run_optim")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("r = $(region_select). initial vs data")

# ## test updates.
# N_d = sum( NMRCalibrate.getNd(Bs[n]) for n = 1:length(Bs) )
# N_β = sum( NMRCalibrate.getNβ(κs_β_DOFs[n], Bs[n]) for n = 1:length(Bs) )
#
# #p = ones(N_d)
# p = -ones(N_d) .* 0.8
#
# st_ind = 1
# NMRCalibrate.updatemixtured!(Bs, p, st_ind, fs, SW, Δsys_cs)
#
# p2 = similar(p)
# NMRCalibrate.applywarptoshifts!(p2, Bs, p, st_ind, a_setp, b_setp)
#
# p = ones(N_β)
# p[end-1] = 0.23
# p[end] = 1.23
# st_ind = 1
# NMRCalibrate.updateβ!(Bs, κs_β_orderings, κs_β_DOFs, p, st_ind)

println("Starting multi-start optim")
N_starts = 100
local_optim_algorithm = NLopt.LN_BOBYQA
xtol_rel = 1e-3
maxeval = 50 # 855 sec.        obj_func(minx) = 0.01690795832619695
# maxeval = 100 # 900 seconds. obj_func(minx) = 0.01690795832619695
maxtime = Inf

shift_lb = -ones(T, N_d)
shift_ub = ones(T, N_d)
shift_initial = zeros(T, N_d)

p_lb = shift_lb
p_ub = shift_ub
p_initial = shift_initial # not used.

prob = MultistartOptimization.MinimizationProblem(obj_func, p_lb, p_ub)

local_method = MultistartOptimization.NLoptLocalMethod(; algorithm = local_optim_algorithm,
xtol_rel = xtol_rel,
maxeval = maxeval,
maxtime = maxtime)

multistart_method = MultistartOptimization.TikTak(N_starts)
@time ret_mo = MultistartOptimization.multistart_minimization(multistart_method,
    local_method, prob)
#
minf = ret_mo.value
minx = ret_mo.location
ret = ret_mo.ret
# 1249 secs.
#
#minx = [-0.09; -0.1; 1.1; 1.1] # a good answer.

minf2 = obj_func(minx)
println("obj_func(minx) = ", minf2)
println()

q2 = uu->NMRSpectraSimulator.evalitpproxymixture(uu, As, Es; w = w)

q_U = q2.(U_rad)
# q0_U = q.(U_rad)
# norm(q_U-q0_U)

PyPlot.figure(fig_num)
fig_num += 1

PyPlot.plot(P_y, real.(y), label = "y")
PyPlot.plot(P, real.(q_U), label = "q")
#PyPlot.plot(P, real.(q3_star_U), "--", label = "run_optim")

PyPlot.legend()
PyPlot.xlabel("ppm")
PyPlot.ylabel("real")
PyPlot.title("r = $(region_select). nested optim result vs data")


text:
using ComputedFieldTypes
using Test

include("readme_examples.jl")
include("inheritance.jl")


text:
module PolynomialBases

using LinearAlgebra

using ArgCheck
using FastGaussQuadrature: gausslegendre, gausslobatto, gaussjacobi
using FastTransforms: clenshawcurtisnodes, clenshawcurtisweights, chebyshevmoments1
using Requires: @require
using UnPack: @unpack
using SpecialFunctions

# types
abstract type AbstractDomain{Dim} end
struct Line <: AbstractDomain{1} end

abstract type AbstractBasis{Domain} end
abstract type NodalBasis{Domain} <: AbstractBasis{Domain} end

# source files
include("canonical_mappings.jl")
include("interpolation.jl")
include("integration.jl")
include("derivative.jl")
include("nodal_bases.jl")
include("legendre.jl")
include("gegenbauer.jl")
include("jacobi.jl")
include("hermite.jl")
include("laguerre.jl")
include("hahn.jl")

# export
## types
export NodalBasis, LobattoLegendre, GaussLegendre, GaussJacobi, ClosedNewtonCotes,
        ClenshawCurtis

## mappings
export map_to_canonical, map_to_canonical!, map_from_canonical, map_from_canonical!

## interpolation
export interpolate, interpolate!, interpolation_matrix, interpolation_matrix!,
        change_basis, change_basis!

## evaluation of coefficients
export compute_coefficients, compute_coefficients!,
        evaluate_coefficients, evaluate_coefficients!

## derivative
export derivative_at, derivative_at!, derivative_matrix, derivative_matrix!

## integration
export integrate

## Continuous orthogonal polynomials
export legendre, legendre_vandermonde, legendre_D, legendre_M,
       gegenbauer,
       jacobi, jacobi_vandermonde,
       hermite,
       laguerre

## Discrete orthogonal polynomials
export hahn

## other utilities
export utility_matrices, includes_boundaries, satisfies_sbp

function __init__()
    # interpolation_matrix
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
        function interpolation_matrix!(mat, dest, src::AbstractVector{SymPy.Sym}, baryweights)
            @boundscheck begin
                @assert length(src) == length(baryweights)
                @assert size(mat,1) == length(dest)
                @assert size(mat,2) == length(src)
            end
            fill!(mat, zero(eltype(mat)))
            @inbounds for k in 1:size(mat,1)
                row_has_match = false
                for j in 1:size(mat,2)
                    # no method matching rtoldefault(::Type{SymPy.Sym})
                    if dest[k] == src[j]
                        row_has_match = true
                        mat[k,j] = 1
                    end
                end
                if row_has_match == false
                    s = zero(eltype(mat))
                    for j in 1:size(mat,2)
                        t = baryweights[j] / (dest[k] - src[j])
                        mat[k,j] = t
                        s = s + t
                    end
                    for j in 1:size(mat,2)
                        mat[k,j] /= s
                    end
                end
            end
            @inbounds for idx in eachindex(mat)
                mat[idx] = SymPy.simplify(mat[idx])
            end
            nothing
        end
    end

    @require SymEngine="123dc426-2d89-5057-bbad-38513e3affd8" begin
        function interpolation_matrix!(mat, dest, src::AbstractVector{SymEngine.Basic}, baryweights)
            @boundscheck begin
                @assert length(src) == length(baryweights)
                @assert size(mat,1) == length(dest)
                @assert size(mat,2) == length(src)
            end
            fill!(mat, zero(eltype(mat)))
            @inbounds for k in 1:size(mat,1)
                row_has_match = false
                for j in 1:size(mat,2)
                    # no method matching rtoldefault(::Type{SymEngine.Basic})
                    if dest[k] == src[j]
                        row_has_match = true
                        mat[k,j] = 1
                    end
                end
                if row_has_match == false
                    s = zero(eltype(mat))
                    for j in 1:size(mat,2)
                        t = baryweights[j] / (dest[k] - src[j])
                        mat[k,j] = t
                        s = s + t
                    end
                    for j in 1:size(mat,2)
                        mat[k,j] /= s
                    end
                end
            end
            @inbounds for idx in eachindex(mat)
                mat[idx] = SymEngine.expand(mat[idx])
            end
            nothing
        end
    end

    # LobattoLegendre
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
        function LobattoLegendre(p::Int, T::Type{SymPy.Sym})
            if p == 0
                nodes   = T[0]
                weights = T[2]
            elseif p == 1
                nodes   = T[-1, 1]
                weights = T[1, 1]
            elseif p == 2
                nodes   = T[-1, 0, 1]
                weights = T[1//3, 4//3, 1//3]
            elseif p == 3
                sqrt_1_5 = sqrt(one(T) / 5)
                nodes   = T[-1, -sqrt_1_5, sqrt_1_5, 1]
                weights = T[1//6, 5//6, 5//6, 1//6]
            elseif p == 4
                sqrt_3_7 = sqrt(T(3) / 7)
                nodes   = T[-1, -sqrt_3_7, 0, sqrt_3_7, 1]
                weights = T[1//10, 49//90, 32//45, 49//90, 1//10]
            elseif p == 5
                sqrt_7 = sqrt(T(7))
                sqrt_m = sqrt(one(T)/3 - 2*sqrt_7/21)
                sqrt_p = sqrt(one(T)/3 + 2*sqrt_7/21)
                nodes   = T[-1, -sqrt_p, -sqrt_m, sqrt_m, sqrt_p, 1]
                weights = T[1//15, (14-sqrt_7)/30, (14+sqrt_7)/30, (14+sqrt_7)/30, (14-sqrt_7)/30, 1//15]
            elseif p == 6
                sqrt_5_3 = sqrt(T(5) / 3)
                sqrt_15 = sqrt(T(15))
                sqrt_m = sqrt((5 - 2*sqrt_5_3)/11)
                sqrt_p = sqrt((5 + 2*sqrt_5_3)/11)
                nodes   = T[-1, -sqrt_p, -sqrt_m, 0, sqrt_m, sqrt_p, 1]
                weights = T[1//21, (124-7*sqrt_15)/350, (124+7*sqrt_15)/350, 256//525, (124+7*sqrt_15)/350, (124-7*sqrt_15)/350, 1//21]
            else
                throw(ArgumentError("Polynomial degree p = $p not implemented yet."))
            end
            baryweights = SymPy.simplify.(barycentric_weights(nodes))
            D = SymPy.simplify.(derivative_matrix(nodes, baryweights))
            LobattoLegendre(nodes, weights, baryweights, D)
        end
    end

    @require SymEngine="123dc426-2d89-5057-bbad-38513e3affd8" begin
        function LobattoLegendre(p::Int, T::Type{SymEngine.Basic})
            if p == 0
                nodes   = T[0]
                weights = T[2]
            elseif p == 1
                nodes   = T[-1, 1]
                weights = T[1, 1]
            elseif p == 2
                nodes   = T[-1, 0, 1]
                weights = T[1//3, 4//3, 1//3]
            elseif p == 3
                sqrt_1_5 = sqrt(one(T) / 5)
                nodes   = T[-1, -sqrt_1_5, sqrt_1_5, 1]
                weights = T[1//6, 5//6, 5//6, 1//6]
            elseif p == 4
                sqrt_3_7 = sqrt(T(3) / 7)
                nodes   = T[-1, -sqrt_3_7, 0, sqrt_3_7, 1]
                weights = T[1//10, 49//90, 32//45, 49//90, 1//10]
            elseif p == 5
                sqrt_7 = sqrt(T(7))
                sqrt_m = sqrt(one(T)/3 - 2*sqrt_7/21)
                sqrt_p = sqrt(one(T)/3 + 2*sqrt_7/21)
                nodes   = T[-1, -sqrt_p, -sqrt_m, sqrt_m, sqrt_p, 1]
                weights = T[1//15, (14-sqrt_7)/30, (14+sqrt_7)/30, (14+sqrt_7)/30, (14-sqrt_7)/30, 1//15]
            elseif p == 6
                sqrt_5_3 = sqrt(T(5) / 3)
                sqrt_15 = sqrt(T(15))
                sqrt_m = sqrt((5 - 2*sqrt_5_3)/11)
                sqrt_p = sqrt((5 + 2*sqrt_5_3)/11)
                nodes   = T[-1, -sqrt_p, -sqrt_m, 0, sqrt_m, sqrt_p, 1]
                weights = T[1//21, (124-7*sqrt_15)/350, (124+7*sqrt_15)/350, 256//525, (124+7*sqrt_15)/350, (124-7*sqrt_15)/350, 1//21]
            else
                throw(ArgumentError("Polynomial degree p = $p not implemented yet."))
            end
            baryweights = SymEngine.expand.(barycentric_weights(nodes))
            D = SymEngine.expand.(derivative_matrix(nodes, baryweights))
            LobattoLegendre(nodes, weights, baryweights, D)
        end
    end

    # GaussLegendre
    @require SymPy="24249f21-da20-56a4-8eb1-6a02cf4ae2e6" begin
        function GaussLegendre(p::Int, T::Type{SymPy.Sym})
            if p == 0
                nodes   = T[0]
                weights = T[2]
            elseif p == 1
                sqrt_1_3 = sqrt(one(T) / 3)
                nodes   = T[-sqrt_1_3, sqrt_1_3]
                weights = T[1, 1]
            elseif p == 2
                sqrt_3_5 = sqrt(T(3) / 5)
                nodes   = T[-sqrt_3_5, 0, sqrt_3_5]
                weights = T[5//9, 8//9, 5//9]
            elseif p == 3
                sqrt_30 = sqrt(T(30))
                sqrt_6_5 = sqrt(T(6)/5)
                sqrt_m = sqrt((3-2*sqrt_6_5)/7)
                sqrt_p = sqrt((3+2*sqrt_6_5)/7)
                nodes   = T[-sqrt_p, -sqrt_m, sqrt_m, sqrt_p]
                weights = T[(18-sqrt_30)/36, (18+sqrt_30)/36, (18+sqrt_30)/36, (18-sqrt_30)/36]
            elseif p == 4
                sqrt_70 = sqrt(T(70))
                sqrt_10_7 = sqrt(T(10)/7)
                sqrt_m = sqrt(5 - 2*sqrt_10_7) / 3
                sqrt_p = sqrt(5 + 2*sqrt_10_7) / 3
                w_m = (322-13*sqrt_70) / 900
                w_p = (322+13*sqrt_70) / 900
                nodes   = T[-sqrt_p, -sqrt_m, 0, sqrt_m, sqrt_p]
                weights = T[w_m, w_p, 128//225, w_p, w_m]
            else
                throw(ArgumentError("Polynomial degree p = $p not implemented yet."))
            end
            baryweights = SymPy.simplify.(barycentric_weights(nodes))
            D = SymPy.simplify.(derivative_matrix(nodes, baryweights))
            R = interpolation_matrix([-1, 1], nodes, baryweights)
            GaussLegendre(nodes, weights, baryweights, D, R[1,:], R[2,:])
        end
    end

    @require SymEngine="123dc426-2d89-5057-bbad-38513e3affd8" begin
        function GaussLegendre(p::Int, T::Type{SymEngine.Basic})
            if p == 0
                nodes   = T[0]
                weights = T[2]
            elseif p == 1
                sqrt_1_3 = sqrt(one(T) / 3)
                nodes   = T[-sqrt_1_3, sqrt_1_3]
                weights = T[1, 1]
            elseif p == 2
                sqrt_3_5 = sqrt(T(3) / 5)
                nodes   = T[-sqrt_3_5, 0, sqrt_3_5]
                weights = T[5//9, 8//9, 5//9]
            elseif p == 3
                sqrt_30 = sqrt(T(30))
                sqrt_6_5 = sqrt(T(6)/5)
                sqrt_m = sqrt((3-2*sqrt_6_5)/7)
                sqrt_p = sqrt((3+2*sqrt_6_5)/7)
                nodes   = T[-sqrt_p, -sqrt_m, sqrt_m, sqrt_p]
                weights = T[(18-sqrt_30)/36, (18+sqrt_30)/36, (18+sqrt_30)/36, (18-sqrt_30)/36]
            elseif p == 4
                sqrt_70 = sqrt(T(70))
                sqrt_10_7 = sqrt(T(10)/7)
                sqrt_m = sqrt(5 - 2*sqrt_10_7) / 3
                sqrt_p = sqrt(5 + 2*sqrt_10_7) / 3
                w_m = (322-13*sqrt_70) / 900
                w_p = (322+13*sqrt_70) / 900
                nodes   = T[-sqrt_p, -sqrt_m, 0, sqrt_m, sqrt_p]
                weights = T[w_m, w_p, 128//225, w_p, w_m]
            else
                throw(ArgumentError("Polynomial degree p = $p not implemented yet."))
            end
            baryweights = SymEngine.expand.(barycentric_weights(nodes))
            D = SymEngine.expand.(derivative_matrix(nodes, baryweights))
            R = interpolation_matrix([-1, 1], nodes, baryweights)
            GaussLegendre(nodes, weights, baryweights, D, R[1,:], R[2,:])
        end
    end
end

end # module


text:
module CompatEntryUtilities

import Pkg

export semver_spec_string

@static if Base.VERSION >= v"1.7-"
    const PKG_VERSIONS = Pkg.Versions
else
    const PKG_VERSIONS = Pkg.Types
end

function _check_result(spec::PKG_VERSIONS.VersionSpec, str::String)
    spec_from_str = try
        PKG_VERSIONS.semver_spec(str)
    catch
        @error("Encountered an error while checking the result", spec, str)
        rethrow()
    end
    result_is_correct = spec == spec_from_str
    if !result_is_correct
        msg = "`Pkg.Versions.semver_spec(str)` is not equal to the original `spec`"
        @error(msg, spec, str, spec_from_str)
        throw(ErrorException(msg))
    end
    return nothing
end

function _fix_zero(s::String)
    if s === "0.0.0"
        return "0"
    end
    return s
end

function _semver_spec_string(r::PKG_VERSIONS.VersionRange)::String
    m, n = r.lower.n, r.upper.n

    if (m, n) == (0, 0)
        return "≥0"
    end

    if ((m, n) == (3, 3)) && (r.lower == r.upper)
        # "=a.b.c"
        return string("=", _fix_zero(join(r.lower.t, ".")))
    end

    if ((m, n) == (3, 1)) && (r.lower.t[1] == r.upper.t[1]) && (r.lower.t[1] ≥ 1)
        # "^a.b.c", where a ≥ 1
        #  "a.b.c", where a ≥ 1
        return _fix_zero(join(r.lower.t[1:m], "."))
    end

    if (m !=0 ) && (n == 0)
        # "≥a.b.c"
        return string("≥", _fix_zero(join(r.lower.t[1:m], ".")))
    end

    if (m !=0 ) && (n != 0)
        # "a.b.c - x.y.z"
        return string(join(r.lower.t[1:m], "."), " - ", join(r.upper.t[1:n], "."))
    end

    msg = "This version range cannot be represented using SemVer notation"
    @error(msg, r, r.lower, r.upper)
    throw(ArgumentError(msg))
end

"""
    semver_spec_string(spec::Pkg.Versions.VersionSpec)::String

Returns `str::AbstractString` such that `Pkg.Versions.semver_spec(str) == spec`.
"""
function semver_spec_string(spec::PKG_VERSIONS.VersionSpec)
    if isempty(spec.ranges)
        result_string = "1 - 0"
    else
        spec_strings = _semver_spec_string.(spec.ranges)
        result_string = join(spec_strings, ", ")
    end
    _check_result(spec, result_string)
    return result_string
end

end # module


text:
using DiffEqFinancial, StochasticDiffEq
using Test

# write your own tests here
u0 = [1.;0.5]
σ =  0.25
prob = HestonProblem(1.,1.,σ,1.,1.,u0,(0.0,1.0))
sol = solve(prob,SRIW1(),adaptive=false,dt=1/10)

prob = BlackScholesProblem((t)->t^2,(t,u)->1.,σ,0.5,(0.0,1.0))
sol = solve(prob,SRIW1())

max_stars_repo_name: JuliaDiffEq/FinancialModels

max_issues_repo_name: JuliaDiffEq/FinancialModels

max_forks_repo_name: JuliaDiffEq/FinancialModels


text:
module GenieTables

using Genie, SearchLight

export AbstractColumn, AbstractTable, render, header
export Column, FieldColumn, LinkColumn

################################################################################
####    Abstract Column                                                     ####
################################################################################
abstract type AbstractColumn end

function render(column::AbstractColumn, object::AbstractModel) end
function header(column::AbstractColumn) end

################################################################################
####    Column                                                              ####
################################################################################
struct Column <: AbstractColumn
    header::String
    template::String
end

function header(column::Column)
    column.header
end

function render(column::Column, object::T) where T <: AbstractModel
    partial(joinpath(Genie.config.path_resources, "tables", "views", column.template); object = object)
end

################################################################################
####    Field column                                                        ####
################################################################################
struct FieldColumn <: AbstractColumn
    header::String
    template::String
    accessor::Symbol
end

function render(column::FieldColumn, object::T) where T <: AbstractModel
    partial(joinpath(Genie.config.path_resources, "tables", "views", column.template);
            object = object, accessor = accessor)
end

################################################################################
####    Link column                                                         ####
################################################################################
struct LinkColumn <: AbstractColumn
    header::String
    template::String
    accessor::Symbol
    view_name::Symbol
end

function render(col::LinkColumn, object::T) where T <: AbstractModel
    partial(joinpath(Genie.config.path_resources, "tables", "views", column.template);
            object = object, accessor = accessor, view_name = view_name)
end

################################################################################
####    Tables
################################################################################
abstract type AbstractTable end

struct AdminTable{T} <: AbstractTable where T
    objects::T
    columns::Vector{<:AbstractColumn}
end

function install(dest::String; force = false) :: Nothing
    src = abspath(normpath(joinpath(@__DIR__, "..", Genie.Plugins.FILES_FOLDER)))

    for f in readdir(src)
        isfile(f) && continue
        isdir(f) || mkpath(joinpath(src, f))

        Genie.Plugins.install(joinpath(src, f), dest, force = force)
    end

    nothing
end

end


text:
module StarGraph
    export H0, Hnoise, Hamiltonian
    """
        H0(N; ν=1.)

    Constructs the noiseless Hamiltonian for the star graph with ``N`` nodes
    """
    function H0(N; ν=1.)
        H = spzeros(N,N)
        H[1,2:end] = ν
        H = H + H'
        H = spdiagm(0 => sum(H,1)[:]) - H
        return H
    end

    """
        Hnoise(N, noise_id)

    Returns the noisy part of the Hamiltonian for the 1D lattice with `N` nodes
    for a specific state of the noise fluctuators specified by `noise_id`.
    """
    function Hnoise(N, noise_id)
        if noise_id < 0 || noise_id > 2^(N-1)
            throw(ArgumentError("noise_id must be between 1 and 2^(N-1)"))
        end

        if noise_id ==0
            return spzeros(N,N)
        end

        noise = 2 * digits(noise_id - 1, base=2, pad=N-1) - 1

        H = spzeros(N,N)
        H[1, 2:end] = noise
        H = H + H'
        H = spdiagm(0 => sum(H,1)[:]) - H
        return H
    end

    """
        Hamiltonian(N, ν, ν0=1.)

    Returns an array of Hamiltonians for all possible configurations of the noise

    # Arguments
        * `N::Integer`: the dimension of the Hilbert space
        * `ν::Float`: the strength of the noise
        * `ν0::Float`: The rate of the noiseless Hamiltonian (ν0 = 1.)
    """
    function Hamiltonian(N::Integer, ν::Float64, ν0=1.)
        return [H0(N;ν=ν0) + ν*Hnoise(N, id) for id = 1:2^(N-1)]
    end

    function H0_spatial_search(N; ν=1., target::Int=1)
        H = spzeros(N,N)
        H[1,2:end] = ν
        H = H + H'
        H = spdiagm(0 => sum(H,1)[:]) - H
        H[target,target] = -1.
        return H
    end

    function Hamiltonian_spatial_search(N::Integer, ν::Float64; ν0=.3, target::Int=1)
        return [H0_spatial_search(N;ν=ν0) + ν*Hnoise(N, id) for id = 1:2^(N-1)]
    end
end


text:
# test multi kde plot, issue 65

using RoME, RoMEPlotting
using Test
# Gadfly.set_default_plot_size(35cm, 25cm)

@testset "test plotKDE array (issue 65)" begin

fg = initfg()

addVariable!(fg, :x0, ContinuousScalar)
addVariable!(fg, :x1, ContinuousScalar)

addFactor!(fg, [:x0], Prior(Normal()))
addFactor!(fg, [:x0;:x1], LinearRelative(Normal()))

ensureAllInitialized!(fg)

plotKDE(fg, ls(fg))

end


text: sqr(x) = x ^ 2

max_stars_repo_path: Module 2/Chapter06/anynomous.jl

max_issues_repo_path: Module 2/Chapter06/anynomous.jl

max_forks_repo_path: Module 2/Chapter06/anynomous.jl


text:
# # clip-11-16.jl

# Load Julia packages (libraries) needed  for the snippets in chapter 0

using StatisticalRethinking, Optim
#gr(size=(600,600));

# ### snippet 3.11

p_grid = range(0, step=0.001, stop=1)
prior = ones(length(p_grid))
likelihood = [pdf(Binomial(3, p), 3) for p in p_grid]
posterior = likelihood .* prior
posterior = posterior / sum(posterior)

# Draw 10000 samples from this posterior distribution

N = 10000
samples = sample(p_grid, Weights(posterior), N);

# In StatisticalRethinkingJulia samples will always be stored
# in an MCMCChains.Chains object.

chn = MCMCChains.Chains(reshape(samples, N, 1, 1), ["toss"]);

# ### snippet 3.12

MCMCChains.describe(chn)

# ### snippet 3.13

MCMCChains.hpd(chn, alpha=0.5)

# ### snippet 3.14

mode(samples)

# ### snippet 3.15

mean(samples)

# ### snippet 3.16

median(samples)

# End of `03/clip-11-16.jl`


text:
function data_gen_unit_tests()

    generate_symmetric_distance_matrix_tests()
    generate_random_pd_points_tests()
    generate_data_random_tsp_tests()
    generate_random_unit_prices()
    generate_random_depot_tests()
    generate_random_vehicle_category_tests()
    generate_random_vehicle_sets_tests()
    generate_full_data_tests()

end

function generate_symmetric_distance_matrix_tests()
    coords = [RVRP.Coord(rand(1:20), rand(1:20)) for i in 1:5]
    matrix_1 = RVRP.generate_symmetric_distance_matrix(coords)
    matrix_2 = RVRP.generate_symmetric_distance_matrix(5)
    @test size(matrix_1) == (5,5)
    @test size(matrix_2) == (5,5)
    for i in 1:5
        @test matrix_1[i,i] == 0.0
        @test matrix_2[i,i] == 0.0
        for j in i+1:5
            @test matrix_1[i,j] == matrix_1[j,i]
            @test matrix_2[i,j] == matrix_2[j,i]
        end
    end
end

function generate_random_pd_points_tests()
    ps = RVRP.generate_random_pd_points(RVRP.PickupPoint, 3, 2)
    @test length(ps) == 3
    @test eltype(ps) == RVRP.PickupPoint
    for p in ps
        @test p.location.index >= 2
        @test p.location.index <= 5
    end
    ps = RVRP.generate_random_pd_points(RVRP.DeliveryPoint, 3, 2)
    @test length(ps) == 3
    @test eltype(ps) == RVRP.DeliveryPoint
    for p in ps
        @test p.location.index >= 2
        @test p.location.index <= 5
    end
end

function generate_data_random_tsp_tests()
    data = RVRP.generate_data_random_tsp(5)
    @test typeof(data) == RVRP.RvrpInstance
    @test data.id[1:11] == "tsp_random_"
    @test length(data.vehicle_categories) == 1
    @test length(data.vehicle_sets) == 1
    @test length(data.depot_points) == 1
    @test length(data.pickup_points) == 4
    @test data.delivery_points == RVRP.DeliveryPoint[]
    @test length(data.requests) == 4
    for i in 1:length(data.requests)
        req = data.requests[i]
        @test req.index == i
        @test req.is_optional == false
        @test req.price_reward == 0.0
        @test req.compartment_capacity_consumption == 0.0
        @test req.split_fulfillment == false
        @test req.precedence_restriction == 0
        @test req.alternative_pickup_point_ids == [string("PickupPoint_", i+1)]
        @test req.alternative_delivery_point_ids == String[]
        @test req.pickup_service_time == 0.0
        @test req.delivery_service_time == 0.0
        @test req.max_duration == typemax(Int32)
        @test req.product_id == "unique_specific_product"
    end
    @test data.travel_time_matrix == Array{Float64,2}(undef, 0, 0)
    @test size(data.travel_distance_matrix) == (5,5)
    for i in 1:5
        @test data.travel_distance_matrix[i,i] == 0.0
        for j in i+1:5
            @test data.travel_distance_matrix[i,j] == data.travel_distance_matrix[j,i]
        end
    end
    @test length(data.product_categories) == 1
    @test length(data.products) == 1
end

function generate_random_unit_prices()
    costs = RVRP.generate_random_unit_prices()
    @test typeof(costs) == RVRP.UnitPrices
end

function generate_random_depot_tests()
    depot = RVRP.generate_random_depot(3, 2)
    @test typeof(depot) == RVRP.DepotPoint
    @test depot.id == "depot_2"
    @test depot.index == 2
    @test depot.location.index == 3
    @test depot.location.id == "depot_3"
end

function generate_random_vehicle_category_tests()
    v_cats = RVRP.generate_random_vehicle_category(5)
    @test eltype(v_cats) == RVRP.VehicleCategory
    @test length(v_cats) == 5
end

function generate_random_vehicle_sets_tests()
    v_cats = RVRP.generate_random_vehicle_category(2)
    depot_points = [RVRP.generate_random_depot(i+1, i) for i in 1:2]
    vs = RVRP.generate_random_vehicle_sets(3, v_cats, depot_points)
    @test length(vs) == 3
    @test eltype(vs) == RVRP.HomogeneousVehicleSet
    for v in vs
        @test v.departure_depot_index >= 1
        @test v.departure_depot_index <= 2
        @test length(v.arrival_depot_indices) == 1
        @test v.departure_depot_index == v.arrival_depot_indices[1]
    end
end

function generate_full_data_tests()
    data = RVRP.generate_full_data_random(3)
    @test typeof(data) == RVRP.RvrpInstance
    @test length(data.pickup_points) == 3
    @test length(data.delivery_points) == 3
    @test length(data.requests) == 6
    for i in 1:length(data.requests)
        @test data.requests[i].index == i
    end
end


text:
isdefined(Base, :__precompile__) && __precompile__()
module PowerSimulations

#################################################################################
# Exports

# Base Models
export Stage
export Simulation
export OperationsProblem
export OperationsProblemTemplate
export InitialCondition
export SimulationSequence

#Network Relevant Exports
export StandardPTDFModel
export CopperPlatePowerModel

######## Device Models ########
export DeviceModel
export FixedOutput
######## Service Models ########
export ServiceModel
export RangeReserve
######## Branch Models ########
export StaticLine
export StaticLineBounds
export StaticLineUnbounded
export StaticTransformer
export StaticTransformerBounds
export StaticTransformerUnbounded
export FlowMonitoredLine
export HVDCLossless
export HVDCDispatch
export HVDCUnbounded
#export VoltageSourceDC
######## Load Models ########
export StaticPowerLoad
export InterruptiblePowerLoad
export DispatchablePowerLoad
######## Renewable Formulations ########
export RenewableFullDispatch
export RenewableConstantPowerFactor
######## Hydro Formulations ########
export HydroDispatchRunOfRiver
export HydroDispatchReservoirFlow
export HydroDispatchReservoirStorage
# export HydroCommitmentRunOfRiver
# export HydroCommitmentReservoirFlow
# export HydroCommitmentReservoirStorage
######## Renewable Formulations ########
export BookKeeping
export BookKeepingwReservation
######## Thermal Formulations ########
export ThermalStandardUnitCommitment
export ThermalBasicUnitCommitment
export ThermalDispatch
export ThermalRampLimited
export ThermalDispatchNoMin

# feedforward chrons
export RecedingHorizon
export Synchronize
export Consecutive

# feedforward models
export UpperBoundFF
export SemiContinuousFF
export RangeFF
export IntegralLimitFF

# InitialConditions chrons
export InterStageChronology
export IntraStageChronology

# Initial Conditions Quantities
export DevicePower
export DeviceStatus
export TimeDurationON
export TimeDurationOFF
export EnergyLevel

# cache_models
export TimeStatusChange
export StoredEnergy

#operation_models
export GenericOpProblem
export UnitCommitmentProblem
export EconomicDispatchProblem
#export OptimalPowerFlow

# Functions
## Construction Exports
export construct_device!
export construct_network!
## Op Model Exports
export solve!
export get_initial_conditions
export set_transmission_model!
export set_devices_template!
export set_branches_template!
export set_services_template!
export set_device_model!
export set_branch_model!
export set_device_model!
## Sim Model Exports
export build!
export execute!
export make_references
## Template Exports
export template_economic_dispatch
export template_unit_commitment
export EconomicDispatchProblem
export UnitCommitmentProblem
export run_economic_dispatch
export run_unit_commitment
## Results interfaces
export get_duals

## Utils Exports
export SimulationResultsReference
export write_results
export check_file_integrity
export load_operation_results
export load_simulation_results
export write_to_CSV
export get_all_constraint_index
export get_all_var_index
export get_con_index
export get_results_variable
export get_var_index
export get_result_variable
export configure_logging
export show_recorder_events
export list_simulation_events
export show_simulation_events

#################################################################################
# Imports
import Logging
import Serialization
#Modeling Imports
import JuMP
# so that users do not need to import JuMP to use a solver with PowerModels
import JuMP: optimizer_with_attributes
export optimizer_with_attributes
import MathOptInterface
import ParameterJuMP
import LinearAlgebra
import PowerSystems
import InfrastructureSystems
# so that users have access to IS.Results interfaces
import InfrastructureSystems:
    get_base_power,
    get_variables,
    get_total_cost,
    get_optimizer_log,
    get_time_stamp,
    write_results
export get_base_power
export get_variables
export get_dual_values
export get_total_cost
export get_optimizer_log
export get_time_stamp
export write_results
import PowerModels
import TimerOutputs

#TimeStamp Management Imports
import Dates
import TimeSeries

#I/O Imports
import DataFrames
import Feather
import JSON
import CSV
import SHA

include("core/definitions.jl")

################################################################################
# Includes

include("logging.jl")
include("utils.jl")

#Models and constructors
include("core/abstract_types.jl")
include("core/aux_structs.jl")

include("services_models/services_model.jl")
include("devices_models/device_model.jl")
include("network_models/networks.jl")

include("core/parameters.jl")
include("core/cache.jl")
include("core/initial_condition_types.jl")
include("core/initial_condition.jl")
include("core/initial_conditions.jl")
include("core/update_initial_conditions.jl")
include("core/operations_problem_template.jl")
include("core/psi_settings.jl")
include("core/psi_container.jl")
include("core/operations_problem_results.jl")
include("core/operations_problem.jl")
include("core/simulation_stages.jl")
include("core/simulation_sequence.jl")
include("core/simulation.jl")
include("core/feedforward.jl")
include("core/simulation_results.jl")
include("core/recorder_events.jl")

#Device Modeling components
include("devices_models/devices/common.jl")
include("devices_models/devices/renewable_generation.jl")
include("devices_models/devices/thermal_generation.jl")
include("devices_models/devices/electric_loads.jl")
include("devices_models/devices/AC_branches.jl")
include("devices_models/devices/DC_branches.jl")
include("devices_models/devices/storage.jl")
include("devices_models/devices/hydro_generation.jl")

#Services Models
include("services_models/reserves.jl")
include("services_models/services_constructor.jl")

#Network models
include("network_models/copperplate_model.jl")
include("network_models/powermodels_interface.jl")
include("network_models/ptdf_model.jl")
include("network_models/network_slack_variables.jl")

#Device constructors
include("devices_models/device_constructors/common/constructor_validations.jl")
include("devices_models/device_constructors/thermalgeneration_constructor.jl")
include("devices_models/device_constructors/hydrogeneration_constructor.jl")
include("devices_models/device_constructors/branch_constructor.jl")
include("devices_models/device_constructors/renewablegeneration_constructor.jl")
include("devices_models/device_constructors/load_constructor.jl")
include("devices_models/device_constructors/storage_constructor.jl")

#Network constructors
include("network_models/network_constructor.jl")

#Templates
include("operations_problems_templates.jl")
# Printing
include("printing.jl")

end


text:
using LinearAlgebra
# positions
x⃗ = [-1, 0, 1] .* 1.0

# position operator
x̂ = Diagonal(x⃗)

# construct differentiation matrix / position operator
X⃗ = [x⃗[i]^(j-1) for i in eachindex(x⃗), j in eachindex(x⃗)]
dX⃗ = [(j-1) * x⃗[i]^(j-2) for i in eachindex(x⃗), j in eachindex(x⃗)]
dX⃗[:,1] .= 0.0 # get rid of NaNs
# p̂ * X⃗ = dX⃗, so just invert to find p̂
p̂ = dX⃗ / X⃗

# fake identity
# won't be identiy for polynomial orders ≥ length(x⃗)-1
I_s =  p̂ * x̂ - x̂ * p̂

# and check that its correct to machine precision
check = 3 .+ x⃗
I_s * check - check

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
using NamedDimensions
using FactCheck, FunctionalData

facts("basic") do
    data = [1 2 3;4 5 6]
    n = named(data, :a, :b)
    r = named(n, :a => 2)
    @fact r.data --> [4,5,6]
    @fact r.names --> [:b]
    @fact array(r) --> [4,5,6]

    @fact ndims(n) --> 2
    @fact size(n) --> (2,3)
    @fact size(n,1) --> 2
    @fact size(n,:a) --> 2
    @fact size(n,2) --> 3
    @fact size(n,:b) --> 3
    @fact length(n) --> 6
    @fact len(n) --> 3

    r = @p named n :a
    @fact r.names --> [:b, :a]
    @fact size(r.data) --> (3,2)

    n = named(data, :c)
    @fact n.names --> [:dimA, :c]

    b = ones(2,3,4)
    n = named(b, :a, :b, :c)
    r = @p named n :a :b
    @fact r.names --> [:c, :a, :b]
    @fact size(r.data) --> (4,2,3)

    r = @p named n :a :b | named :c :a
    @fact r.names --> [:b, :c, :a]
    @fact size(r.data) --> (3,4,2)

    r = @p named data :a :b | named (:b=>2:3)
    @fact r.names --> [:a, :b]
    @fact r.data --> [2 3; 5 6]
    r = @p named data :a :b | named (:a=>1:2)
    @fact r.names --> [:b, :a]
    @fact r.data --> data'

    n = named(data, :a, :b)
    @fact n[:a] --> at(n,:a)
    @fact (@p n :b (2:3)) --> named(n, :b => 2:3)
    @fact (@p n :b (2:3) :a 2) --> named([5,6], :b)

    context("linear indexing") do
        @fact n[1] --> 1
        @fact n[2:3] --> [4,2]
    end

    buf = IOBuffer()
    showinfo(buf,n)
    @fact takebuf_string(buf) --> "----  2 a x 3 b\n    type: Array{Int64,2}   size: (2,3)\n    min:  1   max: 6\n    mean: 3.5   median: 3.5\n"
    showinfo(buf, n, "n in runtests")
    @fact takebuf_string(buf) --> "n in runtests  --  2 a x 3 b\n    type: Array{Int64,2}   size: (2,3)\n    min:  1   max: 6\n    mean: 3.5   median: 3.5\n"

    @fact minimum(n) --> 1
    @fact maximum(n) --> 6
    @fact mean(n) --> 3.5
    @fact minimum(n,1) --> named([1,2,3],:b)
    @fact minimum(n,:a) --> named([1,2,3],:b)
    @fact minimum(n,2) --> named([1,4],:a)
    @fact minimum(n,:b) --> named([1,4],:a)
    @fact maximum(n,2) --> named([3,6],:a)
    @fact std(n,2) --> named(vec(std(n.data,2)), :a)
    @fact std(n,:a) --> named(vec(std(n.data,1)), :b)
    @fact n .+ n --> named([2 4 6;8 10 12], :a, :b)
    @fact n .+ [1 2 3] --> named([2 4 6; 5 7 9], :a, :b)
    @fact array(named(ones(2,3,4)) .+ ones(2)) --> 2*ones(2,3,4)
    @fact array(named(ones(2,3,4)) .+ ones(2,1)) --> 2*ones(2,3,4)
    @fact array(named(ones(2,3,4)) .+ ones(2,1,1)) --> 2*ones(2,3,4)
    @fact n .+ [1 2]' --> named([2 3 4; 6 7 8], :a, :b)
    @fact [1 2 3] .+ n --> named([2 4 6; 5 7 9], :a, :b)
    @fact [1 2]' .+ n --> named([2 3 4; 6 7 8], :a, :b)
    @fact n .* data --> named(data.^2, :a, :b)
    @fact n * data' --> named([14 32; 32 77], :a, :dimB)
    @fact (@p minimum n :b | minimum :a) --> 1

    @fact NamedDimensions.dimnames(named(rand(2,3),:a,:b), 1) --> [:dimA]
    @fact NamedDimensions.dimnames(named(rand(2,3),:dimA,:b), 1) --> [:dimB]
    @fact NamedDimensions.dimnames(named(rand(2,3),:dimA,:dimB), 1) --> [:dimC]
    @fact NamedDimensions.dimnames(named(rand(2,3),:dimC,:dimB), 1) --> [:dimA]
    @fact NamedDimensions.dimnames(named(rand(2,3),:dimC,:dimB), 3) --> [:dimA, :dimD, :dimE]

end

facts("FunctionalData") do
    data = [1 2 3;4 5 6]
    n = named(data, :a, :b)

    @fact concat(n,n) --> named(concat(data,data), :a, :b)
    @fact flatten([n,n]) --> named(concat(data,data), :a, :b)
    @fact concat(n,n) --> named(concat(data,data), :a, :b)
    @fact flatten([n,n]) --> named(concat(data,data), :a, :b)
    @fact stack([n,n], :c) --> named(stack(Any[data,data]), :a, :b, :c)
    @fact unstack(stack([n,n], :c)) --> [named(data, :a, :b), named(data, :a, :b)]

    @fact fst(n) --> named([1,4], :a)
    @fact last(n) --> named([3,6], :a)
    @fact last(n,:b) --> named([3,6], :a)
    @fact last(n,:a) --> named([4,5,6], :b)
    @fact take(n,2) --> named([1 2; 4 5], :a, :b)
    @fact takelast(n,2) --> named([2 3; 5 6], :a, :b)
    @fact part(n,3:3) --> named(col([3;6]), :a, :b)
    @fact drop(n,2) --> part(n, 3:3)

    @fact named(n,:a=>:end) --> named([4,5,6],:b)
    @fact named(n,:b=>:end) --> named([3,6], :a)

    # @fact (@p map n id) --> n
    # @fact (@p map n array) --> n.data

    # @fact (@p map2 n n plus) --> named(data+data, :a, :b)
    # @fact (@p map2 n data plus) --> named(data+data, :a, :b)
    # @fact (@p map2 data n plus) --> named(data+data, :a, :b)
end

s = getstats()
if s["nNonSuccessful"] == 0
    print("   PASSED!")
else
    print("   FAILED: $(s["nFailures"]) failures and $(s["nErrors"]) errors")
end
println("   ( $(s["nNonSuccessful"]+s["nSuccesses"]) tests for runtests.jl $(join(ARGS, " ")))")
exit(s["nNonSuccessful"])


text:
# module for representing agent
# decide control order to robot
# get action based on policy
# by dynamic programming
# using augmented mdp(amdp)

mutable struct AmdpPolicyAgent
  speed
  yaw_rate
  delta_time
  estimator
  prev_spd
  prev_yr
  puddle_coef
  puddle_depth
  total_reward
  in_goal
  final_value
  goal
  pose_min
  pose_max
  widths
  index_nums
  dev_borders
  policy_data
  stop_timer

  # init
  function AmdpPolicyAgent(;delta_time::Float64=0.1,
                           estimator=nothing,
                           goal=nothing,
                           puddle_coef=100,
                           widths=[0.2, 0.2, pi/18],
                           lower_left=[-4.0, -4.0],
                           upper_right=[4.0, 4.0],
                           dev_borders=[0.1, 0.2, 0.4, 0.8])
    self = new()
    self.speed = 0.0
    self.yaw_rate = 0.0
    self.delta_time = delta_time
    self.estimator = estimator
    self.prev_spd = 0.0
    self.prev_yr = 0.0
    self.puddle_coef = puddle_coef
    self.puddle_depth = 0.0
    self.total_reward = 0.0
    self.in_goal = false
    self.final_value = 0.0
    self.goal = goal

    self.pose_min = [lower_left[1], lower_left[2], 0.0]
    self.pose_max = [upper_right[1], upper_right[2], 2*pi]
    self.widths = widths
    index_nums_3d = round.(Int64, (self.pose_max - self.pose_min)./self.widths)
    self.index_nums = [index_nums_3d[1], index_nums_3d[2], index_nums_3d[3], length(dev_borders)]
    self.dev_borders = dev_borders
    self.policy_data = init_belief_policy(self)
    self.stop_timer = 0.0

    return self
  end
end

function init_belief_policy(self::AmdpPolicyAgent)
  tmp = zeros(Tuple([self.index_nums[1], self.index_nums[2], self.index_nums[3], self.index_nums[4], 2]))

  txt_path = "src/decision_making/partially_observable_mdp/policy_amdp.txt"
  open(joinpath(split(@__FILE__, "src")[1], txt_path), "r") do fp
    for line in eachline(fp)
      d = split(line) # [i_x, i_y, i_theta, i_sigma, speed, yaw_rate]
      tmp[parse(Int64, d[1])+1, parse(Int64, d[2])+1, parse(Int64, d[3])+1, parse(Int64, d[4])+1, :] .= [parse(Float64, d[5]), parse(Float64, d[6])]
    end
  end

  return tmp
end

function cov_to_index(self::AmdpPolicyAgent, cov)
  sigma = det(cov)^(1.0/6)

  for (i, e) in enumerate(self.dev_borders)
    if sigma < e
      return i-1
    end
  end

  return length(self.dev_borders)-1
end

function reward_per_sec(self::AmdpPolicyAgent)
  return (-1.0 - self.puddle_depth*self.puddle_coef)
end

function to_index(self::AmdpPolicyAgent, pose)
  index = Int64.(floor.((pose - self.pose_min)./self.widths))

  # normalize direction index
  index[3] = (index[3] + self.index_nums[3]*1000)%self.index_nums[3]

  for i in 1:2
    if index[i] < 0
      index[i] = 0
    elseif index[i] >= self.index_nums[i]
      index[i] = self.index_nums[i] - 1
    end
  end

  return index
end

function policy(self::AmdpPolicyAgent, pose, cov)
  pose_index = to_index(self, pose)
  belief_index = cov_to_index(self, cov)

  action = self.policy_data[pose_index[1]+1, pose_index[2]+1, pose_index[3]+1, belief_index+1, :]

  return action[1], action[2]
end

function draw_decision!(self::AmdpPolicyAgent, observation)
  # reached goal
  if self.in_goal == true
    return 0.0, 0.0 # stop
  end

  if self.estimator !== nothing
    motion_update(self.estimator, self.prev_spd, self.prev_yr, self.delta_time)
    observation_update(self.estimator, observation)

    self.total_reward += self.delta_time * reward_per_sec(self)

    self.speed, self.yaw_rate = policy(self, self.estimator.estimated_pose, self.estimator.estimated_cov)
    self.prev_spd, self.prev_yr = self.speed, self.yaw_rate

    draw!(self.estimator)
  end

  return self.speed, self.yaw_rate
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
@with_kw mutable struct AdversarialOffPolicySolver <: Solver
    𝒮_pro::OffPolicySolver # Solver parameters for the protagonist
    𝒮_ant::OffPolicySolver # Solver parameters for the antagonist
    px::PolicyParams # Nominal disturbance policy
    train_pro_every::Int = 1
    train_ant_every::Int = 1
    log::Union{Nothing, LoggerParams} = nothing # The logging parameters
    i::Int = 0 # The current number of environment interactions
end

function POMDPs.solve(𝒮::AdversarialOffPolicySolver, mdp)
    # make some assertions on the parameters that should be shared
    @assert 𝒮.𝒮_pro.buffer == 𝒮.𝒮_ant.buffer
    @assert 𝒮.𝒮_pro.buffer_init == 𝒮.𝒮_ant.buffer_init
    @assert 𝒮.𝒮_pro.required_columns == 𝒮.𝒮_ant.required_columns
    @assert 𝒮.𝒮_pro.N == 𝒮.𝒮_ant.N
    @assert 𝒮.𝒮_pro.S == 𝒮.𝒮_ant.S
    @assert 𝒮.𝒮_pro.max_steps == 𝒮.𝒮_ant.max_steps

    # define shared parameters for easy reference
    buffer = 𝒮.𝒮_pro.buffer
    buffer_init = 𝒮.𝒮_pro.buffer_init
    ΔN = 𝒮.𝒮_pro.ΔN + 𝒮.𝒮_ant.ΔN
    N = 𝒮.𝒮_pro.N
    S = 𝒮.𝒮_pro.S
    max_steps = 𝒮.𝒮_pro.max_steps

    # Construct the training buffern for the protagonist and antagonist
    𝒟_pro = buffer_like(buffer, capacity=𝒮.𝒮_protagonist.c_opt.batch_size, device=device(𝒮.𝒮_protagonist.π))
    𝒟_ant = buffer_like(buffer, capacity=𝒮.𝒮_protagonist.c_opt.batch_size, device=device(𝒮.𝒮_protagonist.π))

    # Get the discount factor
    γ = Float32(discount(mdp))

    # Two samplers, one for the traditional ap
    s_pro = Sampler(mdp, 𝒮.𝒮_pro.agent, adversary=𝒮.px, S=S, max_steps=max_steps, required_columns=extra_columns(buffer))
    s_ant = Sampler(mdp, 𝒮.𝒮_pro.agent, adversary=𝒮.𝒮_ant.agent, S=S, max_steps=max_steps, required_columns=extra_columns(buffer))

    s_nom = Sampler(mdp_nom, 𝒮.π, S=𝒮.S, A=𝒮.A, max_steps=𝒮.max_steps, required_columns=setdiff(extra_columns(𝒮.buffer), [:x]))
    push!(𝒮.log.fns, (;kwargs...) -> Dict("nominal_undiscounted_return" => undiscounted_return(s_nom, Neps=10)))

    isnothing(𝒮.log.sampler) && (𝒮.log.sampler = s_pro)

    # Log the pre-train performance
    log(𝒮.log, 𝒮.i)

    # Fill the buffer with initial observations before training
    𝒮.i += fillto!(buffer, s_pro, buffer_init, i=𝒮.i, explore=true)
    𝒮.𝒮_pro.i = 𝒮.i
    𝒮.𝒮_ant.i = 𝒮.i

    # Loop over the desired number of environment interactions
    for 𝒮.i in range(𝒮.i, stop=𝒮.i + N - ΔN, step=ΔN)
        # Update the per solver iteration
        𝒮.𝒮_pro.i = 𝒮.i
        𝒮.𝒮_ant.i = 𝒮.i

        # Sample transitions into the replay buffer
        push!(buffer, steps!(s_pro, Nsteps=𝒮.𝒮_pro.ΔN, explore=true, i=𝒮.i))
        push!(buffer, steps!(s_ant, Nsteps=𝒮.𝒮_ant.ΔN, explore=true, i=𝒮.i))

        # callback for potentially updating the buffer
        𝒮.post_experience_callback(buffer)

        # Train the networks
        infos = []
        elapsed(𝒮.i + 1:𝒮.i + ΔN, 𝒮.train_pro_every) && merge!(infos, train_step(𝒮.𝒮_pro, 𝒟_pro, γ))
        elapsed(𝒮.i + 1:𝒮.i + ΔN, 𝒮.train_ant_every) && merge!(infos, train_step(𝒮.𝒮_ant, 𝒟_ant, γ))

        # Log the results
        log(𝒮.log, 𝒮.i + 1:𝒮.i + ΔN, aggregate_info(infos))
    end
    𝒮.i += 𝒮.ΔN
    𝒮.𝒮_pro.i = 𝒮.i
    𝒮.𝒮_ant.i = 𝒮.i

    𝒮.𝒮_pro.π
end


text:
export Category2, FreeCategory2, Hom2, Hom2Expr, composeV, composeH, ⋆,
  DoubleCategory, FreeDoubleCategory, HomH, HomV, HomHExpr, HomVExpr,
  left, right, top, bottom, idH, idV, id2, id2V, id2H,
  MonoidalDoubleCategory, SymmetricMonoidalDoubleCategory,
  FreeSymmetricMonoidalDoubleCategory, braidH, braidV, σH, σV

abstract type Hom2Expr{T} <: CategoryExpr{T} end
abstract type HomVExpr{T} <: CategoryExpr{T} end
abstract type HomHExpr{T} <: CategoryExpr{T} end

# 2-category
############

""" Theory of (strict) *2-categories*
"""
@signature Category2{Ob,Hom,Hom2} <: Category{Ob,Hom} begin
  """ 2-morphism in a 2-category """
  Hom2(dom::Hom(A,B), codom::Hom(A,B))::TYPE ⊣ (A::Ob, B::Ob)
  @op begin
    (⇒) := Hom2
    (⋆) := composeH
  end

  # Hom categories: vertical composition
  id(f)::(f ⇒ f) ⊣ (A::Ob, B::Ob, f::(A ⇒ B))
  compose(α::(f ⇒ g), β::(g ⇒ h))::(f ⇒ h) ⊣
    (A::Ob, B::Ob, f::(A → B), g::(A → B), h::(A → B))

  # Horizontal compostion
  composeH(α::(f ⇒ g), β::(h ⇒ k))::((f ⋅ h) ⇒ (g ⋅ k)) ⊣
    (A::Ob, B::Ob, C::Ob, f::(A → B), g::(A → B), h::(B → C), k::(B → C))
end

# Convenience constructors
composeH(αs::AbstractVector) = foldl(composeH, αs)
composeH(α, β, γ, αs...) = composeH([α, β, γ, αs...])

""" Syntax for a 2-category.

Checks domains of morphisms but not 2-morphisms.
"""
@syntax FreeCategory2{ObExpr,HomExpr,Hom2Expr} Category2 begin
  compose(f::Hom, g::Hom) = associate_unit(new(f,g; strict=true), id)
  compose(α::Hom2, β::Hom2) = associate_unit(new(α,β), id)
  composeH(α::Hom2, β::Hom2) = associate(new(α,β))
end

function show_unicode(io::IO, expr::Hom2Expr{:compose}; kw...)
  Syntax.show_unicode_infix(io, expr, "⋅"; kw...)
end
function show_unicode(io::IO, expr::Hom2Expr{:composeH}; kw...)
  Syntax.show_unicode_infix(io, expr, "*"; kw...)
end

function show_latex(io::IO, expr::Hom2Expr{:compose}; kw...)
  Syntax.show_latex_infix(io, expr, "\\cdot"; kw...)
end
function show_latex(io::IO, expr::Hom2Expr{:composeH}; kw...)
  Syntax.show_latex_infix(io, expr, "*"; kw...)
end

# Double category
#################

""" Theory of (strict) *double categories*
"""
@theory DoubleCategory{Ob,HomV,HomH,Hom2} begin
  Ob::TYPE
  """ Vertical morphism in a double category """
  HomV(dom::Ob, codom::Ob)::TYPE
  """ Horizontal morphism in a double category """
  HomH(dom::Ob, codom::Ob)::TYPE
  """ 2-cell in a double category """
  Hom2(top::HomH(A,B), bottom::HomH(C,D),
       left::HomV(A,C), right::HomV(B,D))::TYPE ⊣ (A::Ob, B::Ob, C::Ob, D::Ob)
  @op begin
    (→) := HomH
    (↓) := HomV
    (⇒) := Hom2
    (⋆) := composeH
    (⋅) := composeV
  end

  idH(A::Ob)::(A → A) ⊣ (A::Ob)
  idV(A::Ob)::(A ↓ A) ⊣ (A::Ob)
  composeH(f::(A → B), g::(B → C))::(A → C) ⊣ (A::Ob, B::Ob, C::Ob)
  composeV(f::(A ↓ B), g::(B ↓ C))::(A ↓ C) ⊣ (A::Ob, B::Ob, C::Ob)

  # Category axioms for horizontal morphisms
  ((f ⋆ g) ⋆ h == f ⋆ (g ⋆ h)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob, f::(A → B), g::(B → C), h::(C → D)))
  f ⋆ idH(B) == f ⊣ (A::Ob, B::Ob, f::(A → B))
  idH(A) ⋆ f == f ⊣ (A::Ob, B::Ob, f::(A → B))

  # Category axioms for vertical morphisms
  ((f ⋅ g) ⋅ h == f ⋅ (g ⋅ h)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob, f::(A ↓ B), g::(B ↓ C), h::(C ↓ D)))
  f ⋅ idV(B) == f ⊣ (A::Ob, B::Ob, f::(A ↓ B))
  idV(A) ⋅ f == f ⊣ (A::Ob, B::Ob, f::(A ↓ B))

  # Identity 2-cell on object
  id2(X::Ob)::Hom2(idH(X), idH(X), idV(X), idV(X)) ⊣ (X::Ob)
  # Identity 2-cell for vertical composition
  id2V(f::(X→Y))::Hom2(f, f, idV(X), idV(Y)) ⊣ (X::Ob, Y::Ob)
  # Identity 2-cell for horizontal composition
  id2H(f::(X↓Y))::Hom2(idH(X), idH(Y), f, f) ⊣ (X::Ob, Y::Ob)

  # Vertical composition of 2-cells
  composeV(α::Hom2(t,b,l,r), β::Hom2(b,b′,l′,r′))::Hom2(t, b′, l⋅l′, r⋅r′) ⊣
    (A::Ob, B::Ob, X::Ob, Y::Ob, C::Ob, D::Ob,
     t::(A→B), b::(X→Y), l::(A↓X), r::(B↓Y),
     b′::(C→D), l′::(X↓C), r′::(Y↓D))

  # Horizontal composition of 2-cells
  composeH(α::Hom2(t,b,l,r), β::Hom2(t′,b′,r,r′))::Hom2(t⋆t′, b⋆b′, l, r′) ⊣
    (A::Ob, B::Ob, X::Ob, Y::Ob, C::Ob, D::Ob,
     t::(A→X), b::(B→Y), l::(A↓B), r::(X↓Y),
     t′::(X→C), b′::(Y→D), r′::(C↓D))
end

# Convenience constructors
composeV(αs::AbstractVector) = foldl(composeV, αs)
composeV(α, β, γ, αs...) = composeV([α, β, γ, αs...])

""" Syntax for a double category.

Checks domains of morphisms but not 2-morphisms.
"""
@syntax FreeDoubleCategory{ObExpr,HomVExpr,HomHExpr,Hom2Expr} DoubleCategory begin
  composeV(f::HomV, g::HomV) = associate_unit(new(f,g; strict=true), idV)
  composeH(f::HomH, g::HomH) = associate_unit(new(f,g; strict=true), idH)
  composeV(α::Hom2, β::Hom2) = associate_unit(new(α,β), id2V)
  composeH(α::Hom2, β::Hom2) = associate_unit(new(α,β), id2H)
end

function show_unicode(io::IO, expr::Hom2Expr{:composeV}; kw...)
  Syntax.show_unicode_infix(io, expr, "⋅"; kw...)
end

function show_latex(io::IO, expr::Hom2Expr{:composeV}; kw...)
  Syntax.show_latex_infix(io, expr, "\\cdot"; kw...)
end

# Monoidal double category
##########################

""" Theory of *monoidal double categories*

To avoid associators and unitors, we assume the monoidal double category is
*strict* in both the horizontal and vertical directions. Apart from assuming
strictness, this theory follows the definition of a monoidal double category in
(Shulman, 2010, Constructing symmetric monoidal bicategories) and other recent
papers, starting from an internal category (S,T: D₁ ⇉ D₀, U: D₀ → D₁,
⋆: D₁ ×_{D₀} D₁ → D₁) in Cat where

- the objects of D₀ are objects
- the morphisms of D₀ are vertical 1-cells
- the objects of D₁ are horizontal 1-cells
- the morphisms of D₁ are 2-cells.

The top and bottom of a 2-cell are given by domain and codomain in D₁ and the
left and right are given by the functors S,T. In a monoidal double category, D₀
and D₁ are each required to be monoidal categories, subject to further axioms
such as S and T being strict monoidal functors.

Despite the apparent asymmetry in this setup, the definition of a monoidal
double category unpacks to be nearly symmetric with respect to horizontal and
vertical, except that the monoidal unit I of D₀ induces the monoidal unit of D₁
as U(I), which I think has no analogue in the vertical direction.
"""
@theory MonoidalDoubleCategory{Ob,HomV,HomH,Hom2} <: DoubleCategory{Ob,HomV,HomH,Hom2} begin
  otimes(A::Ob, B::Ob)::Ob
  otimes(f::(A → B), g::(C → D))::((A ⊗ C) → (B ⊗ D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)
  otimes(f::(A ↓ B), g::(C ↓ D))::((A ⊗ C) ↓ (B ⊗ D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)
  otimes(f::Hom2(t,b,l,r),g::Hom2(t′,b′,l′,r′))::Hom2(t⊗t′,b⊗b′,l⊗l′,r⊗r′) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob, G::Ob, H::Ob,
     t::(A → B), b::(C → D), l::(A ↓ C), r::(B ↓ D),
     t′::(E → F), b′::(G → H), l′::(E ↓ G), r′::(F ↓ H))

  @op (⊗) := otimes
  munit()::Ob

  # Monoid axioms, vertical.
  (A ⊗ B) ⊗ C == A ⊗ (B ⊗ C) ⊣ (A::Ob, B::Ob, C::Ob)
  A ⊗ munit() == A ⊣ (A::Ob)
  munit() ⊗ A == A ⊣ (A::Ob)
  (f ⊗ g) ⊗ h == f ⊗ (g ⊗ h) ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
                                f::(A ↓ X), g::(B ↓ Y), h::(C ↓ Z))

  # Monoid axioms, horizontal.
  (f ⊗ g) ⊗ h == f ⊗ (g ⊗ h) ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
                                f::(A → X), g::(B → Y), h::(C → Z))
  f ⊗ idH(munit()) == f ⊣ (A::Ob, B::Ob, f::(A → B))
  idH(munit()) ⊗ f == f ⊣ (A::Ob, B::Ob, f::(A → B))
  (α ⊗ β) ⊗ γ == α ⊗ (β ⊗ γ) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob,
     G::Ob, H::Ob, I::Ob, J::Ob, K::Ob, L::Ob,
     t1::(A → B), b1::(C → D), l1::(A ↓ C), r1::(B ↓ D),
     t2::(E → F), b2::(G → H), l2::(E ↓ G), r2::(F ↓ H),
     t3::(I → J), b3::(K → L), l3::(I ↓ K), r3::(J ↓ L),
     α::Hom2(t1, b1, l1, r1), β::Hom2(t2, b2, l2, r2), γ::Hom2(t3, b3, l3, r3))

  # Functorality axioms.
  ((f ⊗ g) ⋆ (h ⊗ k) == (f ⋆ h) ⊗ (g ⋆ k)
    ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
       f::(A → B), h::(B → C), g::(X → Y), k::(Y → Z)))
  ((f ⊗ g) ⋅ (h ⊗ k) == (f ⋅ h) ⊗ (g ⋅ k)
    ⊣ (A::Ob, B::Ob, C::Ob, X::Ob, Y::Ob, Z::Ob,
       f::(A ↓ B), h::(B ↓ C), g::(X ↓ Y), k::(Y ↓ Z)))
  ((α ⊗ β) ⋅ (γ ⊗ δ) == (α ⋅ γ) ⊗ (β ⋅ δ)
    ⊣ (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob, G::Ob, H::Ob,
       I::Ob, J::Ob, K::Ob, L::Ob, M::Ob, N::Ob, O::Ob, P::Ob,
       t1::(A → B), b1::(C → D), l1::(A ↓ C), r1::(B ↓ D),
       t2::(E → F), b2::(G → H), l2::(E ↓ G), r2::(F ↓ H),
       t3::(I → J), b3::(K → L), l3::(I ↓ K), r3::(J ↓ L),
       t4::(M → N), b4::(O → P), l4::(M ↓ O), r4::(N ↓ P),
       α::Hom2(t1, b1, l1, r1), β::Hom2(t2, b2, l2, r2),
       γ::Hom2(t3, b3, l3, r3), δ::Hom2(t4, b4, l4, r4)))
  idH(A ⊗ B) == idH(A) ⊗ idH(B) ⊣ (A::Ob, B::Ob)
  idV(A ⊗ B) == idV(A) ⊗ idV(B) ⊣ (A::Ob, B::Ob)
  id2(A ⊗ B) == id2(A) ⊗ id2(B) ⊣ (A::Ob, B::Ob)
  id2H(f ⊗ g) == id2H(f) ⊗ id2H(g) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                      f::(A↓C), g::(B↓D))
  id2V(f ⊗ g) == id2V(f) ⊗ id2V(g) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                      f::(A→C), g::(B→D))
end

""" Theory of (strict) *symmetric monoidal double categories*

Unlike classical double categories, symmetric monoidal double categories do not
treat the vertical and horizontal directions on an equal footing, even in the
strict case. See [`MonoidalDoubleCategory`](@ref) for details and references.
"""
@theory SymmetricMonoidalDoubleCategory{Ob,HomV,HomH,Hom2} <: MonoidalDoubleCategory{Ob,HomV,HomH,Hom2} begin
  braidV(A::Ob, B::Ob)::((A ⊗ B) ↓ (B ⊗ A))
  braidH(f::(A → C), g::(B → D))::Hom2((f ⊗ g), (g ⊗ f), σV(A,B), σV(C,D)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob)
  @op (σV) := braidV
  @op (σH) := braidH

  # Involutivity axioms.
  σV(A,B) ⋅ σV(B,A) == idV(A ⊗ B) ⊣ (A::Ob, B::Ob)
  σH(f,g) ⋅ σH(g,f) == id2V(f ⊗ g) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                      f::(A → C), g::(B → D))

  # Naturality axioms.
  (f⊗g) ⋅ σV(C,D) == σV(A,B) ⋅ (g⊗f) ⊣ (A::Ob, B::Ob, C::Ob, D::Ob,
                                        f::(A ↓ C), g::(B ↓ D))
  ((α⊗β) ⋅ σH(h,k) == σH(f,g) ⋅ (β⊗α) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob, G::Ob, H::Ob,
     f::(A → C), g::(B → D), h::(E → G), k::(F → H),
     ℓ1::(A ↓ E), r1::(C ↓ G), ℓ2::(B ↓ F), r2::(D ↓ H),
     α::Hom2(f,h,ℓ1,r1), β::Hom2(g,k,ℓ2,r2)))

  # Coherence axioms.
  σV(A,B⊗C) == (σV(A,B) ⊗ idV(C)) ⋅ (idV(B) ⊗ σV(A,C)) ⊣ (A::Ob, B::Ob, C::Ob)
  σV(A⊗B,C) == (idV(A) ⊗ σV(B,C)) ⋅ (σV(A,C) ⊗ idV(B)) ⊣ (A::Ob, B::Ob, C::Ob)
  (σH(f,g⊗h) == (σH(f,g) ⊗ id2V(h)) ⋅ (id2V(g) ⊗ σH(f,h)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob,
     f::(A → D), g::(B → E), h::(C → F)))
  (σH(f⊗g,h) == (id2V(f) ⊗ σH(g,h)) ⋅ (σH(f,h) ⊗ id2V(g)) ⊣
    (A::Ob, B::Ob, C::Ob, D::Ob, E::Ob, F::Ob,
     f::(A → D), g::(B → E), h::(C → F)))
end

@syntax FreeSymmetricMonoidalDoubleCategory{ObExpr,HomVExpr,HomHExpr,Hom2Expr} SymmetricMonoidalDoubleCategory begin
  otimes(A::Ob, B::Ob) = associate_unit(new(A,B), munit)
  otimes(f::HomV, g::HomV) = associate(new(f,g))
  otimes(f::HomH, g::HomH) = associate(new(f,g))
  otimes(f::Hom2, g::Hom2) = associate(new(f,g))
  composeV(f::HomV, g::HomV) = associate_unit(new(f,g; strict=true), idV)
  composeH(f::HomH, g::HomH) = associate_unit(new(f,g; strict=true), idH)
  composeV(α::Hom2, β::Hom2) = associate_unit(new(α,β), id2V)
  composeH(α::Hom2, β::Hom2) = associate_unit(new(α,β), id2H)
end

function show_unicode(io::IO, expr::HomVExpr{:braidV}; kw...)
  Syntax.show_unicode_infix(io, expr, "σV"; kw...)
end
function show_unicode(io::IO, expr::Hom2Expr{:braidH}; kw...)
  Syntax.show_unicode_infix(io, expr, "σH"; kw...)
end

function show_latex(io::IO, expr::HomVExpr{:braidV}; kw...)
  Syntax.show_latex_script(io, expr, "\\sigmaV")
end
function show_latex(io::IO, expr::Hom2Expr{:braidH}; kw...)
  Syntax.show_latex_script(io, expr, "\\sigmaH")
end


text:
julia> x = [1, 2, 3]
julia> ptr = pointer_from_objref(x)
Ptr{Void} @0x000000010282e4a0
julia> unsafe_pointer_to_objref(ptr)
3-element Array{Int64,1}:
 1
 2
 3

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
macro localscribe(f::Symbol, dict=:(IdDict()))
    sname = scribename(f)
    esc(quote
        !@isdefined($sname) && ($sname = Anamnesis.Scribe($f, $dict))
        $sname
    end)
end

macro scribe(f::Symbol, dict=:(IdDict()))
    sname = scribename(f)
    esc(quote
        !@isdefined($sname) && (global $sname = Anamnesis.Scribe($f, $dict))
        $sname
    end)
end

macro localmem(expr::Expr)
    @capture(expr, f_(args__))
    esc(:(@localscribe($f)($(args...))))
end

macro mem(expr::Expr)
    @capture(expr, f_(args__))
    esc(:(@scribe($f)($(args...))))
end

argname(x)::Symbol = splitarg(x)[1]

function _translate_kwarg(ex::Expr)
    arg = argname(ex)
    Expr(:kw, arg, arg)
end

macro anamnesis(expr::Expr)
    origdef = splitdef(expr)
    fname = origdef[:name]

    rawfuncdef = copy(origdef)
    rawname = rawfuncname(fname)
    rawfuncdef[:name] = rawname

    scrname = scribename(rawname)

    funcdef = copy(origdef)
    kwargs = _translate_kwarg.(funcdef[:kwargs])
    funcdef[:body] = quote
        $scrname($(map(argname, funcdef[:args])...); $(kwargs...))
    end

    esc(quote
        $(MacroTools.combinedef(rawfuncdef))
        $scrname = Anamnesis.Scribe($rawname,
                                    IdDict{Any, $(get(origdef, :rtype, :Any))}())
        $(MacroTools.combinedef(funcdef))
    end)
end

macro anamnesis(expr::Expr, funcs::Symbol...)
    expr = MacroTools.postwalk(expr) do ex
        if @capture(ex, f_(args__)) && (f ∈ funcs)
            :(@mem $ex)
        else
            ex
        end
    end
    esc(expr)
end


text:
using PNML, PNML.PnmlTypes
#using PrettyPrinting
using Documenter

using PNML
using PNML: AbstractHLCore,
    AbstractLabel,
    AbstractPnmlTool,
    Arc,
    Condition,
    Coordinate,
    AbstractDeclaration,
    AbstractSort,
    BuiltInSort,
    MultisetSort,
    ProductSort,
    UserSort,
    Term,
    Variable,
    AbstractOperator,
    BuiltInOperator,
    BuiltInConst,
    MultiSetOperator,
    PnmlTuple,
    UserOperator,
    Fill,
    Font,
    Graphics,
    HLCore,
    HLAnnotation,
    HLInscription,
    HLMarking,
    HLNet,
    HLPetriNet,
    IDRegistry,
    Line,
    MalformedException,
    Maybe,
    MissingIDException,
    Name,
    ObjectCommon,
    OpenNet,
    PNML,
    PTInscription,
    PTMarking,
    PTNet,
    PT_HLPNG,
    Page,
    PetriNet,
    Place,
    PnmlCfg,
    PnmlCore,
    PnmlDict,
    PnmlException,
    PnmlLabel,
    PnmlModel,
    PnmlNet,
    PnmlNode,
    PnmlObject,
    PnmlType,
    PnmlTypes,
    RefPlace,
    RefTransition,
    SimpleNet,
    StochasticNet,
    SymmetricNet,
    TimedNet,
    TokenGraphics,
    ToolInfo,
    Transition,
    XMLNode,
    _harvest_any!,
    _match,
    add_label!,
    add_toolinfo!,
    all_arcs,
    allchildren,
    append_page!,
    arc,
    arc_ids,
    arcs,
    condition,
    conditions,
    convert,
    deref!,
    deref_place,
    deref_transition,
    duplicate_id_action,
    find_net,
    find_nets,
    first_net,
    firstchild,
    firstpage,
    flatten_pages!,
    get_label,
    get_toolinfo,
    has_arc,
    has_graphics,
    has_label,
    has_labels,
    has_name,
    has_place,
    has_refP,
    has_refT,
    has_structure,
    has_text,
    has_toolinfo,
    has_tools,
    has_transition,
    has_xml,
    in_out,
    inc_indent,
    indent,
    indent_width,
    infos,
    initialMarking,
    ins,
    inscription,
    isregistered,
    marking,
    nets,
    number_value,
    outs,
    parse_and,
    parse_arbitraryoperator,
    parse_arbitrarysort,
    parse_arc,
    parse_bool,
    parse_booleanconstant,
    parse_condition,
    parse_declaration,
    parse_declarations,
    parse_equality,
    parse_file,
    parse_graphics,
    parse_graphics_coordinate,
    parse_graphics_fill,
    parse_graphics_font,
    parse_graphics_line,
    parse_hlinitialMarking,
    parse_hlinscription,
    parse_imply,
    parse_inequality,
    parse_initialMarking,
    parse_inscription,
    parse_label,
    parse_mulitsetsort,
    parse_name,
    parse_namedoperator,
    parse_net,
    parse_node,
    parse_not,
    parse_or,
    parse_page,
    parse_place,
    parse_pnml,
    parse_pnml_common!,
    parse_pnml_label_common!,
    parse_pnml_node_common!,
    parse_productsort,
    parse_refPlace,
    parse_refTransition,
    parse_sort,
    parse_str,
    parse_structure,
    parse_term,
    parse_text,
    parse_tokengraphics,
    parse_tokenposition,
    parse_toolspecific,
    parse_transition,
    parse_tuple,
    parse_type,
    parse_unparsed,
    parse_useroperator,
    parse_usersort,
    parse_variable,
    parse_variabledecl,
    pid,
    place,
    place_ids,
    places,
    pnml_common_defaults,
    pnml_label_defaults,
    pnml_namespace,
    pnml_node_defaults,
    pnml_ns,
    rate,
    rates,
    refplace,
    refplace_ids,
    refplaces,
    reftransition,
    reftransition_ids,
    reftransitions,
    register_id!,
    reset_registry!,
    show_common,
    show_page_field,
    source,
    src_arcs,
    structure,
    tag,
    tagmap,
    target,
    text,
    tgt_arcs,
    transition,
    transition_function,
    transition_ids,
    transitions,
    nettype,
    unclaimed_label,
    update_maybe!,
    xmlnode

# Makie.jl is a source of many of these good ideas. (Bad ones are mine?)

################################################################################
#                              Utility functions                               #
################################################################################

################################################################################
#                                    Setup                                     #
################################################################################

pathroot   = normpath(@__DIR__, "..")
docspath   = joinpath(pathroot, "docs")
srcpath    = joinpath(docspath, "src")
buildpath  = joinpath(docspath, "build")
genpath    = joinpath(srcpath,  "generated")
srcgenpath = joinpath(docspath, "src_generation")

# Eventually we plan on generating pictures, et al in genpath.

mkpath(genpath) #TODO where should initialization happen?

################################################################################
#                          Syntax highlighting theme                           #
################################################################################

#TODO

################################################################################
#                      Automatic Markdown page generation                      #
################################################################################

#TODO

################################################################################
#                 Building HTML documentation with Documenter                  #
################################################################################

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

for m ∈ [PNML]
    for i ∈ propertynames(m)
       xxx = getproperty(m, i)
       println(xxx)
    end
 end

@info("Running `makedocs` from make.jl.")

makedocs(;
         clean = true,
         doctest=true,
         modules=[PNML], #, PNML.PnmlTypes],
         authors="Jeff Hurst <strangehurst@users.noreply.github.com>",
         repo="https://github.com/strangehurst/PNML.jl/blob/{commit}{path}#{line}",
         #repo="/home/jeff/PNML/{path}",
         checkdocs=:all,

         format=Documenter.HTML(;
                                # CI means publish documentation on GitHub.
                                prettyurls=get(ENV, "CI", nothing) == "true",
                                canonical="https://strangehurst.github.io/PNML.jl",
                                assets=String[],
                                sidebar_sitename=true,
                                prerender=false,
                                #no highlight.js
                                ),
         sitename="PNML.jl",
         pages=[
            "Petri Net Markup Language" => "pnml.md",
            "Intermediate Representation" => "IR.md",
            "Interfaces" => "interface.md",
            "API" => [
                "PNML"      => "API/library.md",
                "PnmlTypes" => "API/pnmltypes.md"
            ],
            "Examples"   => "examples.md",
            "Index" => "index.md",
            "acknowledgments.md",
          ],
         )

################################################################################
#                           Deploying documentation                            #
################################################################################

if !isempty(get(ENV, "DOCUMENTER_KEY", ""))
    deploydocs(;
               repo="github.com/strangehurst/PNML.jl",
               devbranch = "main",
               push_preview = true,
               )
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -2.9922365365158297 <= q <= 1.3959867501971397)

                     add_NL_constraint(m, :(log(1 + exp(-0.7172866311917714 + 0.8530094465195428*log(1 + exp(0.49975691401319544 + 0.817850285997455*log(1 + exp(-0.9044530867886094 + -0.3504940249833268*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + -0.925387475816613*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + -0.15011340339301915*log(1 + exp(-0.49866226232151334 + 0.31623027699891626*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + 0.1823838704363938*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))) + -0.2302663292873568*log(1 + exp(0.6056562038214484 + 0.017530899400252764*log(1 + exp(-0.9044530867886094 + -0.3504940249833268*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + -0.925387475816613*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + 0.9141090581590978*log(1 + exp(-0.49866226232151334 + 0.31623027699891626*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + 0.1823838704363938*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))))) + log(1 + exp(0.11802092343800519 + -0.23198020992542823*log(1 + exp(0.49975691401319544 + 0.817850285997455*log(1 + exp(-0.9044530867886094 + -0.3504940249833268*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + -0.925387475816613*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + -0.15011340339301915*log(1 + exp(-0.49866226232151334 + 0.31623027699891626*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + 0.1823838704363938*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))) + 0.3694431087223533*log(1 + exp(0.6056562038214484 + 0.017530899400252764*log(1 + exp(-0.9044530867886094 + -0.3504940249833268*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + -0.925387475816613*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))) + 0.9141090581590978*log(1 + exp(-0.49866226232151334 + 0.31623027699891626*log(1 + exp(0.9239332442221326 + 0.41157080299867177*$(x[1]) + 0.6856618842893334*$(x[2]) + 0.23250449051531152*$(x[3]) + -0.8171094215290706*$(x[4]))) + 0.1823838704363938*log(1 + exp(-0.1441575648301785 + -0.7324773919132506*$(x[1]) + 0.1666248349273194*$(x[2]) + 0.5034013066229974*$(x[3]) + -0.3351866390106042*$(x[4]))))))))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions