text:
# A unit fraction contains 1 in the numerator. The decimal representation of
# the unit fractions with denominators 2 to 10 are given:

# 1/2 =   0.5
# 1/3 =   0.(3)
# 1/4 =   0.25
# 1/5 =   0.2
# 1/6 =   0.1(6)
# 1/7 =   0.(142857)
# 1/8 =   0.125
# 1/9 =   0.(1)
# 1/10    =   0.1
#
# Where 0.1(6) means 0.166666..., and has a 1-digit recurring cycle. It can be
# seen that 1/7 has a 6-digit recurring cycle.
#
# Find the value of d < 1000 for which 1/d contains the longest recurring cycle
# in its decimal fraction part.

using ProjectEulerSolutions

# Solve using greedy regex to pull out the longest repeated sequence of digits.
# Could be made slightly faster to consider only prime values between 2 and n.
function p026solution_regex(n::Integer=3)::Integer
    numdigits = 2500
    startvalue = BigInt(10)^numdigits

    regexstring = r"(\d+?)\1"
    lengths = zeros(Integer, n)

    @simd for i = 3:2:n
        stringnum = string(fld(startvalue, i))
        @inbounds lengths[i] = mapreduce(m->length(m.captures[1]), max, eachmatch(regexstring, stringnum))
    end

    return findmax(lengths)[2]
end

# Cycles through a sequence of digits building up and storing larger numbers
# until it sees one it has seen before.  Could be made slightly faster to
# consider only prime values between 2 and n.
function p026solution_repeats(n::Integer=3)::Integer
    totalmax = 0
    maxind = 1
    for i = 3:2:n
        seennumbers = Set{Integer}()
        r = 10
        j = 0
        while ! in(r, seennumbers)
            push!(seennumbers, r)
            r = 10 * (r % i)
            j += 1
        end

        if j > totalmax
            totalmax = j
            maxind = i
        end
    end

    return maxind
end

p026 = Problems.Problem(Dict("Regex" => p026solution_regex,
                             "Repeats" => p026solution_repeats))

Problems.benchmark(p026, 1000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text: # define a spatial dependency graph (internals, working)

ext: jl


text:
#   Unit "timefrequencyuni" of the FourierAnalysis Package for julia language
#
#   MIT License
#   Copyright (c) 2019, Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# ? CONTENTS :
#   This unit implements time-frequency univariate measure based
#   on timefrequency.jl and tools.jl

################### UNIVARIATE MEASURE ##########################
# All measures are weighted versions of measures described in:
# Congedo(2018): https://hal.archives-ouvertes.fr/hal-01868538v2/document

# Mean Amplitude (MAmp), Eq. 0.1 therein: <|z|> = <|re^i𝜑|> = <r> if w=[].
# If w is a positive weight vector a weighted version <wr>/<w> is computed.
# NB: w must be normalized so that the denominator is 1.0
"""
```
(1)
function meanAmplitude( 𝐀      :: TFAmplitudeVector,
                        frange :: fInterval,
                        trange :: tInterval;
                    mode   :: Function = extract,
                    func   :: Function = identity,
                    w      :: Vector   = [],
                    check  :: Bool     = true)

(2)
function meanAmplitude( 𝐙 :: TFAnalyticSignalVector,
                    < same arguments as method (1) >

(3)
function meanAmplitude( 𝐱      :: Vector{Vector{T}},
                        sr     :: Int,
                        wl     :: Int,
                        frange :: fInterval,
                        trange :: tInterval,
                     bandwidth :: IntOrReal    = 2;
                 mode          :: Function     = extract,
                 func          :: Function     = identity,
                 w             :: Vector       = [],
                 check         :: Bool         = true,
                 filtkind      :: FilterDesign = Butterworth(2),
                 fmin          :: IntOrReal    = bandwidth,
                 fmax          :: IntOrReal    = sr÷2,
                 fsmoothing    :: Smoother     = noSmoother,
                 tsmoothing    :: Smoother     = noSmoother,
                 planner       :: Planner      = getplanner,
                 ⏩           :: Bool         = true) where T<:Real

```
**alias**: mamp

(1)

Given a [TFAmplitudeVector](@ref) object, estimate the
[(weighted) mean amplitude](@ref) measure across those objects.
The time-frequency planes of all the objects in `𝐀` should be congruent.

**arguments**:

`frange` and `trange` define a time-frequency region on which the estimation
is sought. `frange` is a [fInterval](@ref) type
and delimits center frequencies of the filter bank. `trange` is a
[tInterval](@ref) type and delimits time samples. To obtain the estimation
in the whole time-frequency plane use the column (:) sign for both arguments.

**optional keyword arguments**

If `mode=extract` is passed (default), the measure will be computed
for all points in the chosen time-frequency region. If `mode=mean` is passed,
it will be computed on the mean of these points (grand-average).
The [`extract`](@ref) and [`mean`](@ref) functions are
[generic methods](@ref) of *FourierAnalysis*.

**Note**: with `mode=mean` the output of the function is always a real number,
whereas with `mode=extract` the output may be a real number, a real
row or column vector or a real matrix, depending on the shape of the
chosen time-frequency region.

Passing a function with the `func` argument you can
[derive your own time-frequency measures](@ref).

`w` may be a vector of non-negative real weights associated to
each input object. By default the unweighted version of the measure
is computed.

If `check` is true (default), check that if the column sign is passed
- as `frange` argument, all input objects have the same number of rows (center frequencies);
- as `trange` argument, all input objects have the same number of columns (time samples).
If either check fails, print an error message and return `Nothing`.
No other range checks are performed.

(2)

Given a [TFAnalyticSignalVector](@ref) object, compute the amplitude
of all objects in `𝐙` and estimate the [(weighted) mean amplitude](@ref)
measure across those objects as per method (1). In addition,
since using this method all [TFAnalyticSignal](@ref) in `𝐙` must be `linear`,
if `check` is true (default) and this is not the case,
print an error and return `Nothing`. The checks on `frange` and `trange`
performed by method (1) are also performed by this method.

(3)

Estimate the amplitude of all data vectors in `𝐱` calling the
[`TFamplitude`](@ref) constructor and then estimate the
[(weighted) mean amplitude](@ref) measure across the constructed
amplitude objects as per method (1).

`frange`, `trange`, `mode`, `func`, `w` and `check` have the same meaning as
in method (1). The other arguments are passed to the
[`TFamplitude`](@ref) constructor, to which the reader is referred
for their meaning.

**See also**: [`concentration`](@ref), [`meanDirection`](@ref),
[timefrequencybi.jl](@ref).

**Examples**:
```
using FourierAnalysis

# generate 100 vectors of data
sr, t, bandwidth=128, 512, 2
h=taper(harris4, t)
𝐱=[sinusoidal(2, 10, sr, t, 0).*h.y+randn(t) for i=1:100]

𝐘=TFanalyticsignal(𝐱, sr, t, bandwidth; fmax=32)
𝐀=TFamplitude(𝐘)

# mean amplitude in a TF region from a TFAnalyticSignalVector object
MAmp=meanAmplitude(𝐘, (4, 16), :)
heatmap(MAmp; c=:fire) # y axis labels are not correct

# mean amplitude in a TF region from a TFAmplitudeVector object
MAmp=meanAmplitude(𝐀, (4, 16), :)

# mean amplitude in a TF region directly from data
MAmp=meanAmplitude(𝐱, sr, t, (4, 16), :, bandwidth)

# NB: in the above, the analytic signal objects must all
# be linear, since meanAmplitude is computed from amplitude
# and the amplitude of non-linear analytic signal is uniformy equal to 1.

# All these computations can be obtained averaging in a TF region, e.g.,
MAmp=meanAmplitude(𝐘, (4, 16), :; mode=mean) # output a real number

# and can be obtained on smoothed Amplitude, e.g.,
MAmp=meanAmplitude(𝐱, sr, t, (4, 16), :;
                   fsmoothing=blackmanSmoother,
                   tsmoothing=blackmanSmoother)
# or, equivalently, and using the alias `mamp`,
MAmp=mamp(smooth(blackmanSmoother, blackmanSmoother, 𝐀), (4, 16), :)

# A similar syntax is used for the other univariate measures, e.g.,
# concentration averaging in a TF region from a TFAnalyticSignalVector object
ConM=concentration(𝐘, (4, 16), (128, 384); mode=mean)

# concentration in a TF region directly from data (using the alias `con`)
ConE=con(𝐱, sr, t, (4, 16), (128, 384), bandwidth; mode=extract)
heatmap(Con; c=:fire) # y axis labels are not correct

NB: ConM is not at all equivalent to mean(ConE) !

# mean direction averaging in a TF region directly from data
MDir=meanDirection(𝐱, sr, t, (4, 16), :, bandwidth; mode=mean)

# mean direction in a TF region from a TFAnalyticSignalVector object
MDir=meanDirection(𝐘, (4, 16), :)

# and for the non-linear counterpart:
# phase concentration in a TF region directly from data
Con=concentration(𝐱, sr, t, (8, 12), :; nonlinear=true)

# phase concentration at a single TF point
Con=concentration(𝐱, sr, t, 10, 256; nonlinear=true)

# phase mean direction averaging in a TF region directly from data
# and using the alias `mdir`
MDir=mdir(𝐱, sr, t, (8, 12), :; mode=mean, nonlinear=true)

# If you try to compute a non-linear measure from a linear
# AnalyticSignal object you will get en error (see the REPL), e.g.,
Con=con(𝐘, (8, 12), (1, 512); mode=mean, nonlinear=true)

# In order to compute non-linear measures from analytic signal objects
# first we need to compute non-linear analytic signal objects:
𝐘=TFanalyticsignal(𝐱, sr, t, bandwidth; fmax=32, nonlinear=true)

# then, we can obtain for example the phase concentration
Con=con(𝐘, (8, 12), :; mode=mean, nonlinear=true)

# and the phase mean direction
MDir=meanDirection(𝐘, (8, 12), :; nonlinear=true)
```
"""
meanAmplitude( 𝐀      :: TFAmplitudeVector,
               frange :: fInterval,
               trange :: tInterval;
           mode   :: Function = extract,
           func   :: Function = identity,
           w      :: Vector   = [],
           check  :: Bool     = true) =
    if !check || _checkUniData(𝐀, frange, trange, false, true, "meanAmplitude")
        return mean(func(mode(𝐀, frange, trange; w=w, check=check)))
    end

meanAmplitude( 𝐙      :: TFAnalyticSignalVector,
               frange :: fInterval,
               trange :: tInterval;
           mode   :: Function = extract,
           func   :: Function = identity,
           w      :: Vector   = [],
           check  :: Bool     = true) =
    if !check || _checkUniData(𝐙, frange, trange, false, true, "meanAmplitude")
        return mean(func(mode(TFamplitude(𝐙), frange, trange; w=w, check=check)))
    end

# NB  if `t` > 2^14 then `t` is set to 2^10.
meanAmplitude( 𝐱      :: Vector{Vector{T}},
               sr     :: Int,
               wl     :: Int,
               frange :: fInterval,
               trange :: tInterval,
            bandwidth :: IntOrReal    = 2;  # <>
        mode          :: Function     = extract,
        func          :: Function     = identity,
        w             :: Vector       = [],
        check         :: Bool         = true,
        filtkind      :: FilterDesign = Butterworth(2),
        fmin          :: IntOrReal    = bandwidth,
        fmax          :: IntOrReal    = sr÷2,
        fsmoothing    :: Smoother     = noSmoother,
        tsmoothing    :: Smoother     = noSmoother,
        planner       :: Planner      = getplanner,
        ⏩           :: Bool         = true) where T<:Real =

    meanAmplitude(TFamplitude(𝐱,
                              sr,
                              wl,
                              bandwidth;
                        fmin       = fmin,
                        fmax       = fmax,
                        fsmoothing = fsmoothing,
                        tsmoothing = tsmoothing,
                        filtkind   = filtkind,
                        planner    = planner,
                        ⏩        = ⏩),
                  frange,
                  trange;
            mode  = mode,
            func  = func,
            w     = w,
            check = check)

mamp = meanAmplitude

"""
```
(1)
function concentration( 𝐙       :: TFAnalyticSignalVector,
                        frange  :: fInterval,
                        trange  :: tInterval;
                    nonlinear :: Bool     = false,
                    mode      :: Function = extract,
                    func      :: Function = identity,
                    w         :: Vector   = [],
                    check     :: Bool     = true)

(2)
function concentration( 𝐱      :: Vector{Vector{T}},
                        sr     :: Int,
                        wl     :: Int,
                        frange :: fInterval,
                        trange :: tInterval,
                     bandwidth :: IntOrReal    = 2;
                 nonlinear  :: Bool         = false,
                 mode       :: Function     = extract,
                 func       :: Function     = identity,
                 w          :: Vector       = [],
                 check         :: Bool      = true,
                 filtkind   :: FilterDesign = Butterworth(2),
                 fmin       :: IntOrReal    = bandwidth,
                 fmax       :: IntOrReal    = sr÷2,
                 fsmoothing :: Smoother     = noSmoother,
                 tsmoothing :: Smoother     = noSmoother,
                 planner    :: Planner      = getplanner,
                 ⏩        :: Bool         = true) where T<:Real

```

**alias**: con

If optional keyword parameter `nonlinear` is false (default),
estimate the [(weighted) concentration](@ref) measure,
otherwise estimate the [(weighted) phase concentration](@ref) measure.

(1)
The desired measure is obtained averaging across the
[TFAnalyticSignal](@ref) objects in `𝐙`. Since this method uses
pre-computed analytic signal objects, their `.nonlinear` field
must agree with the `nonlinear` argument passed to
this function.

`frange`, `trange`, `w`, `mode` and `func` have the same meaning as
in the [`meanAmplitude`](@ref) function, however keep in mind that
the two possible `mode` functions, i.e.,
[`extract`](@ref) and [`mean`](@ref), in this function operate on
complex numbers.

The checks performed in the [`meanAmplitude`](@ref) function
are performed here too. In addition, if `check` is true, also check that
- if `nonlinear` is true, all objects in `𝐙` are nonlinear;
- if `nonlinear` is false, all objects in `𝐙` are linear.
If either check fails, print an error message and return `Nothing`.

(2)
Estimate the analytic signal of all data vectors in `𝐱` calling the
[`TFanalyticsignal`](@ref) constructor and then use method (1)
to obtained the desired measure.

`frange`, `trange`, `mode`, `func`, `w` and `check` have the same meaning as
in the [`meanAmplitude`](@ref) function.
The other arguments are passed to the [`TFanalyticsignal`](@ref) constructor,
to which the reader is referred for understanding their action.

**See also**: [`meanAmplitude`](@ref), [`meanDirection`](@ref),
[timefrequencybi.jl](@ref).

**Examples**: see the examples of [`meanAmplitude`](@ref) function.
"""
concentration( 𝐙         :: TFAnalyticSignalVector,
               frange    :: fInterval,
               trange    :: tInterval;
               nonlinear :: Bool     = false,
               mode      :: Function = extract,
               func      :: Function = identity,
               w         :: Vector   = [],
               check     :: Bool     = true) =
    if !check || _checkUniData(𝐙, frange, trange, nonlinear, !nonlinear, "concentration")
        return amplitude(mean(func(mode(𝐙, frange, trange; w=w, check=check))))
    end

# NB  if `t` > 2^14 then `t` is set to 2^10.
concentration( 𝐱      :: Vector{Vector{T}},
               sr     :: Int,
               wl     :: Int,
               frange :: fInterval,
               trange :: tInterval,
            bandwidth :: IntOrReal    = 2;
        nonlinear     :: Bool         = false,
        mode          :: Function     = extract,
        func          :: Function     = identity,
        w             :: Vector       = [],
        check         :: Bool         = true,
        filtkind      :: FilterDesign = Butterworth(2),
        fmin          :: IntOrReal    = bandwidth,
        fmax          :: IntOrReal    = sr÷2,
        fsmoothing    :: Smoother     = noSmoother,
        tsmoothing    :: Smoother     = noSmoother,
        planner       :: Planner      = getplanner,
        ⏩           :: Bool         = true) where T<:Real =

    concentration(TFanalyticsignal(𝐱,
                                   sr,
                                   wl,
                                   bandwidth;
                              fmin       = fmin,
                              fmax       = fmax,
                              fsmoothing = fsmoothing,
                              tsmoothing = tsmoothing,
                              filtkind   = filtkind,
                              nonlinear  = nonlinear,
                              planner    = planner,
                              ⏩        = ⏩),
                  frange,
                  trange;
              nonlinear = nonlinear,
              mode      = mode,
              func      = func,
              w         = w,
              check     = check)

con = concentration

"""
```
(1)
function meanDirection( 𝐙         :: TFAnalyticSignalVector,
                        frange    :: fInterval,
                        trange    :: tInterval;
                    nonlinear :: Bool       = false,
                    mode      :: Function   = extract,
                    func      :: Function   = identity,
                    w         :: Vector     = [],
                    check     :: Bool       = true)

(2)
function meanDirection( 𝐱      :: Vector{Vector{T}},
                        sr     :: Int,
                        wl     :: Int,
                        frange :: fInterval,
                        trange :: tInterval,
                     bandwidth :: IntOrReal    = 2;
                 nonlinear     :: Bool         = false,
                 mode          :: Function     = extract,
                 func          :: Function     = identity,
                 w             :: Vector       = [],
                 check         :: Bool         = true,
                 filtkind      :: FilterDesign = Butterworth(2),
                 fmin          :: IntOrReal    = bandwidth,
                 fmax          :: IntOrReal    = sr÷2,
                 fsmoothing    :: Smoother     = noSmoother,
                 tsmoothing    :: Smoother     = noSmoother,
                 planner       :: Planner      = getplanner,
                 ⏩           :: Bool         = true) where T<:Real
```
This function features two methods that use exactly the same syntax
as the two corresponding methods of the [`concentration`](@ref) function.
All arguements have exactly the same meaning as therein.
Only the output differs:

if optional keyword parameter `nonlinear` is false (default),
estimate the [(weighted) mean direction](@ref) measure,
otherwise estimate the [(weighted) phase mean direction](@ref) measure.

**alias**: mdir

**See also**: [`meanAmplitude`](@ref), [`concentration`](@ref),
[timefrequencybi.jl](@ref).

**Examples**: see the examples of [`meanAmplitude`](@ref).
"""
meanDirection( 𝐙         :: TFAnalyticSignalVector,
               frange    :: fInterval,
               trange    :: tInterval;
           nonlinear :: Bool       = false,
           mode      :: Function   = extract,
           func      :: Function   = identity,
           w         :: Vector     = [],
           check     :: Bool       = true) =
    if !check || _checkUniData(𝐙, frange, trange, nonlinear, !nonlinear, "meanDirection")
        return phase(mean(func(mode(𝐙, frange, trange; w=w, check=check))))
    end

# NB  if `t` > 2^14 then `t` is set to 2^10.
meanDirection( 𝐱      :: Vector{Vector{T}},
               sr     :: Int,
               wl     :: Int,
               frange :: fInterval,
               trange :: tInterval,
            bandwidth :: IntOrReal    = 2;
        nonlinear     :: Bool         = false,
        mode          :: Function     = extract,
        func          :: Function     = identity,
        w             :: Vector       = [],
        check         :: Bool         = true,
        filtkind      :: FilterDesign = Butterworth(2),
        fmin          :: IntOrReal    = bandwidth,
        fmax          :: IntOrReal    = sr÷2,
        fsmoothing    :: Smoother     = noSmoother,
        tsmoothing    :: Smoother     = noSmoother,
        planner       :: Planner      = getplanner,
        ⏩           :: Bool         = true) where T<:Real =

    meanDirection(TFanalyticsignal(𝐱,
                                   sr,
                                   wl,
                                   bandwidth;
                               fmin       = fmin,
                               fmax       = fmax,
                               fsmoothing = fsmoothing,
                               tsmoothing = tsmoothing,
                               filtkind   = filtkind,
                               nonlinear  = nonlinear,
                               planner    = planner,
                               ⏩        = ⏩),
                  frange,
                  trange;
              nonlinear = nonlinear,
              mode      = mode,
              func      = func,
              w         = w,
              check     = check)

mdir = meanDirection

# Check that:
# - if `frange` is a Colon (:), then the data in all objects in `𝒀` have the same number of rows,
# - if `trange` is a Colon (:), then the data in all objects in `𝒀` have the same number of columns,
# - if `mustBeAllNonLinear` is `true`, then the `nonlinear` field in all objects in `𝒀` is set to `true`,
# - if `mustBeAllLinear` is `true`, then the `nonlinear` field in all objects in `𝒀` is set to `false`.
# Return `true` if all checks are positive, `false` otherwise.
function _checkUniData(𝒀::TFobjectsVector, frange::fInterval, trange::tInterval,
                       mustBeAllNonLinear::Bool, mustBeAllLinear::Bool, funcname::String)
    OK = true
    frange isa Colon && !_allsame([size(Y.y, 1) for Y ∈ 𝒀]) ?
        (@error 📌*", "*funcname*" function: the number of filter banks is not the same in all input TF objects. Don't use a colon for the `frange` or make sure the input TF objects are of similar size"; OK = false) : nothing
    trange isa Colon && !_allsame([size(Y.y, 2) for Y ∈ 𝒀]) ?
        (@error 📌*", "*funcname*" function: the number of samples is not the same in all input TF objects. Don't use a colon for the `trange` or make sure the input TF objects are of similar size"; OK = false) : nothing
    mustBeAllNonLinear && !isNonLinear(𝒀) ?
        (@error 📌*", all input TF objects must be nonlinear in order to compute "*funcname; OK = false) : nothing
    mustBeAllLinear && !isLinear(𝒀) ?
        (@error 📌*", all input TF objects must be linear in order to compute "*funcname; OK = false) : nothing
    return OK
end


text:
# module Octo

import Base: ==, !=, <, >, <=, >=, -, +, *, /
const PlainTypes = Union{Number, String, Symbol, <:DatePeriod, <:TimePeriod}

for op in (:(==), :(!=), :(<), :(>), :(<=), :(>=), :(-), :(+), :(*), :(/))
    @eval begin
        ($op)(left::SQLElement, right::SQLElement) = Predicate(($op), left, right)
        ($op)(left::SQLElement, right::PlainTypes) = Predicate(($op), left, right)
        ($op)(left::PlainTypes, right::SQLElement) = Predicate(($op), left, right)
        ($op)(left::Union{SQLElement, PlainTypes}, right::Type{PlaceHolder}) = Predicate(($op), left, right)
        ($op)(left::Type{PlaceHolder}, right::Union{SQLElement, PlainTypes}) = Predicate(($op), left, right)
    end
end

function Base.:*(left::SQLKeyword, right::SQLKeyword)
    KeywordAllKeyword(left, right)
end

# module Octo


text:
import Base.Dates: Second

warsaw = resolve("Europe/Warsaw", tzdata["europe"]...)
fixed = FixedTimeZone("Fixed", -7200, 3600)

# ZonedDateTime accessors
zdt = ZonedDateTime(DateTime(2014,6,12,23,59,58,57), fixed)
@test TimeZones.localtime(zdt) == DateTime(2014,6,12,23,59,58,57)
@test TimeZones.utc(zdt) == DateTime(2014,6,13,0,59,58,57)

@test TimeZones.days(zdt) == 735396
@test TimeZones.hour(zdt) == 23
@test TimeZones.minute(zdt) == 59
@test TimeZones.second(zdt) == 58
@test TimeZones.millisecond(zdt) == 57

# Make sure that Base.Dates accessors work with ZonedDateTime.
@test Dates.year(zdt) == 2014
@test Dates.month(zdt) == 6
@test Dates.week(zdt) == 24
@test Dates.day(zdt) == 12
@test Dates.dayofmonth(zdt) == 12
@test Dates.yearmonth(zdt) == (2014, 6)
@test Dates.monthday(zdt) == (6, 12)
@test Dates.yearmonthday(zdt) == (2014, 6, 12)

# Vectorized accessors
# Note: repmat is used over broadcast to test for size and equality.
arr = repmat([zdt], 10)
@test TimeZones.hour(arr) == repmat([23], 10)
@test TimeZones.minute(arr) == repmat([59], 10)
@test TimeZones.second(arr) == repmat([58], 10)
@test TimeZones.millisecond(arr) == repmat([57], 10)

@test Dates.year(arr) == repmat([2014], 10)
@test Dates.month(arr) == repmat([6], 10)
@test Dates.day(arr) == repmat([12], 10)
@test Dates.dayofmonth(arr) == repmat([12], 10)
@test Dates.yearmonth(arr) == repmat([(2014, 6)], 10)
@test Dates.monthday(arr) == repmat([(6, 12)], 10)
@test Dates.yearmonthday(arr) == repmat([(2014, 6, 12)], 10)


text:
using BayesianOptimization, GaussianProcesses
const BO = BayesianOptimization
using Test, Random
include("acquisitionfunctions.jl")
include("acquisition.jl")
include("branin.jl")
include("warmstart.jl")
include("utils.jl")


text:
using Documenter
using DynamicPPL
using DynamicPPL: AbstractPPL

# Doctest setup
DocMeta.setdocmeta!(DynamicPPL, :DocTestSetup, :(using DynamicPPL); recursive=true)

makedocs(;
    sitename="DynamicPPL",
    format=Documenter.HTML(),
    modules=[DynamicPPL],
    pages=["Home" => "index.md", "API" => "api.md"],
    strict=true,
    checkdocs=:exports,
)

deploydocs(; repo="github.com/TuringLang/DynamicPPL.jl.git", push_preview=true)

max_stars_repo_name: TuringLang/MicroPPL

max_issues_repo_name: TuringLang/MicroPPL

max_forks_repo_name: TuringLang/MicroPPL


text:
"""
    bearing_to_azimuth(bearing::Real)::Real

Convert any bearing angle from the north line direction (positive clockwise)
and return an angle between 0-360 degrees (positive clockwise), 0 being the north line.

# Examples
```jldoctest
julia> bearing_to_azimuth(40)
40

julia> bearing_to_azimuth(-105)
255

julia> bearing_to_azimuth(410)
50
```
"""
function bearing_to_azimuth(bearing::Real)::Real
    angle = bearing % 360
    if angle < 0
        angle += 360
    end
    return angle
end

"""
    rhumb_bearing(start::Position, stop::Position, final::Bool=false)

Take two Positions and finds the bearing angle between them along a Rhumb line
i.e. the angle measured in degrees start the north line (0 degrees).

# Examples
```jldoctest
julia> start = Position([-75, 45])
2-element Array{Float64,1}:
 -75.0
  45.0

julia> stop = Position([20, 60])
2-element Array{Float64,1}:
 20.0
 60.0

julia> rhumb_bearing(start, stop)
75.28061364784332
```
"""
function rhumb_bearing(start::Position, stop::Position, final::Bool=false)
    bear360 = nothing

    if final === true
        bear360 = calculateRhumbBearing(stop, start)
    else
        bear360 = calculateRhumbBearing(start, stop)
    end

    bear180 = (bear360 > 180) ? - (360 - bear360) : bear360

    return bear180
end

rhumb_bearing(start::Point, stop::Point, final::Bool=false) = rhumb_bearing(start.coordinates, stop.coordinates, final)

function calculateRhumbBearing(a::Position, b::Position)
    ϕ1 = deg2rad(a[2])
    ϕ2 = deg2rad(b[2])

    Δλ = deg2rad((b[1] - a[1]))

    if Δλ > pi Δλ -= 2 * pi end
    if Δλ < -pi Δλ += 2* pi end

    Δψ = log(tan(ϕ2 / 2 + pi / 4) / tan(ϕ1 / 2 + pi / 4))

    θ = atan(Δλ, Δψ)

    return (rad2deg(θ) + 360) % 360
end

"""
    bearing(start::Position, stop::Position, final::Bool=false)

Take two Points or Positions and finds the geographic bearing between them,
i.e. the angle measured in degrees from the north line (0 degrees).

# Examples
```jldoctest
julia> start = Position([-75, 45])
2-element Array{Float64,1}:
 -75.0
  45.0

julia> stop = Position([20, 60])
2-element Array{Float64,1}:
 20.0
 60.0

julia> bearing(start, stop)
37.75495852601734
```
"""
function bearing(start::Position, stop::Position, final::Bool=false)
    if final === true
        bear = bearing(stop, start, false)
        return (bear + 180) % 360
    end

    lon1 = deg2rad(start[1])
    lon2 = deg2rad(stop[1])
    lat1 = deg2rad(start[2])
    lat2 = deg2rad(stop[2])

    a = sin(lon2 - lon1) * cos(lat2)
    b = cos(lat1) * sin(lat2) - sin(lat1)  * cos(lat2) * cos(lon2 - lon1)

    return rad2deg(atan(a, b))
end

bearing(start::Point, stop::Point, final::Bool=false) = bearing(start.coordinates, stop.coordinates, final)


text:
@testset "Test files with only one source of data (point or analog)" begin

    @testset "Point only files" begin
        samp29_1 = @test_nowarn readc3d(joinpath(datadir, "sample29", "Facial-Sing.c3d"))
        @test isempty(samp29_1.analog)
        samp29_2 = @test_nowarn readc3d(joinpath(datadir, "sample29", "OptiTrack-IITSEC2007.c3d"))
        @test isempty(samp29_1.analog)
        samp34_1 = @test_nowarn readc3d(joinpath(datadir, "sample34", "AlbertoINT.c3d"))
        @test isempty(samp34_1.analog)
        samp34_2 = @test_nowarn readc3d(joinpath(datadir, "sample34", "AlbertoREAL.c3d"))
        @test isempty(samp34_2.analog)
        samp34_3 = @test_nowarn readc3d(joinpath(datadir, "sample34", "Basketball.c3d"))
        @test isempty(samp34_3.analog)
    end

    @testset "Analog only files" begin
        samp35 = @test_nowarn readc3d(joinpath(datadir, "sample35", "Mega Electronics Isokinetic EMG Angle Torque Sample File.c3d"))
        @test isempty(samp35.point)
    end

end


text:
using NLopt
using Optim
using ADCME
using Statistics
using LinearAlgebra
using PyCall
using SparseArrays
using Test
using Random

conda = get_conda()
run_with_env(`$conda install -y -c anaconda matplotlib`)
using PyPlot

if has_gpu()
    use_gpu(0)
end
doctor()
sess = Session()

include("layers.jl")
include("sparse.jl")
include("random.jl")
include("io.jl")
include("variable.jl")
include("ops.jl")
include("core.jl")
include("extra.jl")
include("newton.jl")
include("optim.jl")
include("ot.jl")
include("ode.jl")
include("flow.jl")
include("mpi.jl")
include("toolchain.jl")
include("rbf.jl")
include("pcl.jl")
# The default matplotlib backend does not work for MacOSX
if !Sys.isapple()
    include("kit.jl")
end


text:
"""
A Julia suite for chaos and nonlinear dynamics
"""
module DynamicalSystems

using StaticArrays

export SVector, SMatrix, @SVector, @SMatrix, Size

using Reexport

@reexport using DelayEmbeddings
@reexport using DynamicalSystemsBase
@reexport using ChaosTools
@reexport using RecurrenceAnalysis

display_update = true
update_name = "update_v1.3.0"

if display_update
if !isfile(joinpath(@__DIR__, update_name))
printstyled(stdout,
"""
\nUpdate message: DynamicalSystems v1.3

A method that estimates the predictability properties of a
dynamical system has been implemented, following the work of:

Wernecke, H., Sándor, B. & Gros, C.
*How to test for partially predictable chaos*.

See the function `predictability`.\n
"""; color = :light_magenta)
touch(joinpath(@__DIR__, update_name))
end
end

end # module


text:
const ENLARGEMENT = 10

struct Image
    data::Any
    width::Int
    height::Int
    is_mirrored::Bool

    function Image(_data = nothing, _is_mirrored = false)

        data = _data
        if !isnothing(data)
            width = data.width
            height = data.height
        else
            width = 0
            height = 0
        end
        is_mirrored = _is_mirrored

        new(
            data,
            width,
            height,
            is_mirrored,
        )
    end
end

function mirror(image::Image)

    return Image(image.data, !image.is_mirrored)
end

function show(image::Image, rectangle)

    rectangle = rectangle * ENLARGEMENT
    left_top = rectangle[2]
    w = rectangle[3].x - left_top.x
    h = rectangle[1].y - left_top.y

    # Show rectangle in which image is placed
    #setcolor("black")
    #rect(left_top, w, h, :stroke)

    if image.is_mirrored
        right_top = Point(left_top.x + w, left_top.y)
        translate(right_top)
        scale(-1 * w / image.width, h / image.height)
    else
        translate(left_top)
        scale(w / image.width, h / image.height)
    end

    placeimage(image.data, Point(0, 0))

    # Reverse previous translation and scaling
    origin(Point(0, 0))
    return nothing
end


text:
import Flux
using Flux
using Flux: params, flatten
import Zygote
using Statistics
import LinearAlgebra
import BSON
using BSON: @save
using BSON: @load
using Dates
using Serialization

function ResNetBlock(depth; filter_size=(3,3), momentum=0.1f0)
    pad = filter_size .÷ 2
    layers = Chain(
        Conv(filter_size, depth=>depth, pad=pad),
        BatchNorm(depth, relu; momentum=momentum),
        Conv(filter_size, depth=>depth, pad=pad),
        BatchNorm(depth; momentum=momentum))
    Chain(
        SkipConnection(layers, +),
        x -> relu.(x))
end

mutable struct Incremental <: Network
    config::MuZeroConfig
    saved_at::Union{Nothing, String}
    step::Integer
    test_mode::Bool
    force_cpu::Bool
    representation_network::Chain
    prediction_common_network::Chain
    prediction_policy_network::Chain
    prediction_value_network::Chain
    dynamics_hidden_state_network::Chain
    dynamics_reward_network::Chain
    function Incremental(config; test_mode=false, force_cpu=false)
        game_type = config.game_type
        hidden_size = config.hidden_state_size
        observation_resolution = config.observation_shape
        initial_filter_size = (3,3)
        input_depth = config.observation_shape[3]
        pad = initial_filter_size .÷ 2
        num_filters = config.num_filters
        flat_size = observation_resolution[1] * observation_resolution[2] * num_filters

        rep_layers = config.representation_resnet_layers
        representation_net = Chain(
            Conv(initial_filter_size, input_depth=>num_filters, pad=pad),
            BatchNorm(num_filters, relu, momentum=0.1f0),
            [ResNetBlock(num_filters) for _ in 1:rep_layers]...,
        )

        dynamics_layers = config.dynamics_resnet_layers
        dynamics_hidden_state_net = Chain(
            Conv(initial_filter_size, (num_filters+1)=>num_filters, pad=pad),
            BatchNorm(num_filters, relu, momentum=0.1f0),
            [ResNetBlock(num_filters) for _ in 1:dynamics_layers]...,
        )

        reward_layers = config.reward_fc_layers
        reward_net = Chain(
            Conv((1,1), num_filters=>num_filters),
            BatchNorm(num_filters, relu, momentum=0.1f0),
            flatten,
            [Dense(flat_size, flat_size, relu) for _ in 1:reward_layers]...,
            Dense(flat_size, full_support_size(config), relu),
            Dense(full_support_size(config), 1, tanh))

        prediction_common = Chain(
            ResNetBlock(num_filters))

        policy_layers = config.policy_fc_layers
        prediction_policy = Chain(
            Conv((1,1), num_filters=>num_filters),
            BatchNorm(num_filters, relu, momentum=0.1f0),
            flatten,
            [Dense(flat_size, flat_size, relu) for _ in 1:policy_layers]...,
            Dense(flat_size, action_space_size(game_type)))

        value_layers = config.value_fc_layers
        prediction_value = Chain(
            Conv((1,1), num_filters=>num_filters),
            BatchNorm(num_filters, relu, momentum=0.1f0),
            flatten,
            [Dense(flat_size, flat_size, relu) for _ in 1:value_layers]...,
            Dense(flat_size, full_support_size(config), relu),
            Dense(full_support_size(config), 1, tanh))

        device(net) = force_cpu ? cpu(net) : gpu(net)
        x = new(
            config,
            nothing,
            0,
            test_mode,
            force_cpu,
            device(representation_net),
            device(prediction_common),
            device(prediction_policy),
            device(prediction_value),
            device(dynamics_hidden_state_net),
            device(reward_net),
        )
        enable_test_mode_if_needed!(x)
        x
    end
end

function get_nets(x::Incremental)
    [getfield(x, n) for n in get_net_names(x)]
end

function enable_test_mode_if_needed!(n::Incremental)
    testmode!.(get_nets(n), n.test_mode)
end

zeropadder(padto) = (x) -> Matrix{Float32}(I, padto, size(x, 1)) * x
device(n::Incremental, x) = n.force_cpu ? cpu(x) : gpu(x)

function even_1_odd_n1(x)
    res = x / Float32(2)
    floor(res) ≈ res ? Float32(1) : Float32(-1)
end

function prediction(n::Incremental, hidden_state::AbstractArray{Float32})
    common = n.prediction_common_network(hidden_state)
    policy_logits = n.prediction_policy_network(common)
    state_values = n.prediction_value_network(common)
    policy_logits, state_values
end

function initial_inference(n::Incremental, games::AbstractArray{MuZeroGame})
    initial_inference(
        n,
        current_image.(getgame.(games))
    )
end

initial_inference(n::Incremental, game::MuZeroGame) = initial_inference(n, game.game)
initial_inference(n::Incremental, game::AbstractEnv) = initial_inference(n, current_image(game))
initial_inference(n::Incremental, game::Vector{<:AbstractEnv}) = initial_inference(n, current_image.(game))
initial_inference(n::Incremental, image::AbstractImage) = initial_inference(n, tensorize(image))
initial_inference(n::Incremental, image::Vector{<:AbstractImage}) = initial_inference(n, tensorize(image))

function initial_inference(n::Incremental, s::AbstractEnvState)
    initial_inference(
        n,
        get_image(s),
    )
end

#t = TTT(0)
#out = initial_inference(pttt, current_image(t)) |> SingleNetworkOutput

function initial_inference(n::Incremental, images::AbstractArray{Float32, 4})
    hidden_state = n.representation_network(device(n, images))
    policy_logits, value = prediction(n, hidden_state)
    reward = n.dynamics_reward_network(hidden_state)
    value, reward, policy_logits, hidden_state
end

function dynamics(n::Incremental, hidden_state::AbstractArray{Float32, 4}, actions::AbstractArray{Float32, 4})
    #example: size(actions) = (3,3,1,9)
    hidden_state_plus_actions = cat(hidden_state, actions; dims=3)
    next_hidden_state = n.dynamics_hidden_state_network(hidden_state_plus_actions)
    rewards = n.dynamics_reward_network(next_hidden_state)
    next_hidden_state, rewards
end

# out = initial_inference(pttt, current_image(t)) |> SingleNetworkOutput
# next_states, rewards = dynamics(pttt, hidden_state(out), [2])
# dynamics(pttt, images, [2,2])

function recurrent_inference(n::Incremental, hidden_states::Vector{<:AbstractArray{Float32, 3}}, actions::Vector{<:Integer})
    recurrent_inference(
        n,
        reduce((x,y)->cat(x,y; dims=4), Flux.unsqueeze.(hidden_states,4)),
        reduce(
            (x,y)->cat(x,y; dims=4),
            Flux.unsqueeze.(
                map(a->get_action_image(n.config.game_type, a).channels, actions),
                4
            )
        ),
    )
end

function recurrent_inference(n::Incremental, hidden_state::AbstractArray{Float32, 3}, action::Integer)
    action_tensor = Flux.unsqueeze(get_action_image(n.config.game_type, action).channels, 4)
    recurrent_inference(n, Flux.unsqueeze(hidden_state, 4), action_tensor)
end

function recurrent_inference(n::Incremental, hidden_state::AbstractArray{Float32, 4}, action::AbstractArray{Float32, 4})
    next_hidden_state, reward = dynamics(n, hidden_state, action)
    policy_logits, value = prediction(n, next_hidden_state)
    value, reward, policy_logits, next_hidden_state
end

get_weights(x::Incremental) = params(get_nets(x)...)

# function save!(n::Incremental, optimizer, path::String)
#     ensure_pathto_exists(path)
#     saved_at = string(now())
#     n.saved_at = saved_at
#     opt                     = cpu(optimizer)
#     representation          = cpu(n.representation_network)
#     prediction_common       = cpu(n.prediction_common_network)
#     prediction_policy       = cpu(n.prediction_policy_network)
#     prediction_value        = cpu(n.prediction_value_network)
#     dynamics_hidden_state   = cpu(n.dynamics_hidden_state_network)
#     dynamics_reward         = cpu(n.dynamics_reward_network)
#     @save expanduser(path) opt saved_at representation prediction_common prediction_policy prediction_value dynamics_hidden_state dynamics_reward
# end

function copy!(from::Incremental, to::Incremental)
    to.saved_at                         = from.saved_at
    to.step                             = from.step
    # @todo - need to make sure device actually *copies* for both cpu and gpu
    to.representation_network           = device(to, from.representation_network)
    to.prediction_common_network        = device(to, from.prediction_common_network)
    to.prediction_policy_network        = device(to, from.prediction_policy_network)
    to.prediction_value_network         = device(to, from.prediction_value_network)
    to.dynamics_hidden_state_network    = device(to, from.dynamics_hidden_state_network)
    to.dynamics_reward_network          = device(to, from.dynamics_reward_network)
    enable_test_mode_if_needed!(to)
    return
end

function save!(n::Incremental, path::String)
    ensure_pathto_exists(path)
    #     d = Dict()
    #     for x in get_net_names(n)
    #         d[x] = params(getfield(n, x))
    #     end
    #     @show keys(d)
    #     BSON.bson(expanduser(path), d)
    saved_at = string(now())
    n.saved_at = saved_at
    step = n.step
    representation          = cpu(n.representation_network)
    prediction_common       = cpu(n.prediction_common_network)
    prediction_policy       = cpu(n.prediction_policy_network)
    prediction_value        = cpu(n.prediction_value_network)
    dynamics_hidden_state   = cpu(n.dynamics_hidden_state_network)
    dynamics_reward         = cpu(n.dynamics_reward_network)
    @save expanduser(path) saved_at step representation prediction_common prediction_policy prediction_value dynamics_hidden_state dynamics_reward
    dir = join(splitpath(expanduser(path))[1:(end-1)], "/")
    serialize("$dir/step.bin", step)
    serialize("$dir/saved_at.bin", saved_at)
    serialize("$dir/representation.bin", representation)
    serialize("$dir/prediction_common.bin", prediction_common)
    serialize("$dir/prediction_policy.bin", prediction_policy)
    serialize("$dir/prediction_value.bin", prediction_value)
    serialize("$dir/dynamics_hidden_state.bin", dynamics_hidden_state)
    serialize("$dir/dynamics_reward.bin", dynamics_reward)
end

get_net_names(n::Incremental) = get_net_names(typeof(n))

function get_net_names(::Type{Incremental})
    [
        :representation_network,
        :prediction_common_network,
        :prediction_policy_network,
        :prediction_value_network,
        :dynamics_reward_network,
        :dynamics_hidden_state_network,
    ]
end

# function load!(n::Incremental, d::Dict)
#     n.saved_at                      = d[:saved_at]
#     n.representation_network        = device(n, d[:representation])
#     n.prediction_common_network     = device(n, d[:prediction_common])
#     n.prediction_policy_network     = device(n, d[:prediction_policy])
#     n.prediction_value_network      = device(n, d[:prediction_value])
#     n.dynamics_reward_network       = device(n, d[:dynamics_reward])
#     n.dynamics_hidden_state_network = device(n, d[:dynamics_hidden_state])
#     enable_test_mode_if_needed!(n)
#     nothing
# end

function load!(n::Incremental, path::String)
    # load!(
    #     n,
    #     BSON.load(expanduser(path), @__MODULE__)
    # )
    dir = join(splitpath(expanduser(path))[1:(end-1)], "/")
    n.step                          = device(n, deserialize("$dir/step.bin"))
    n.saved_at                      = device(n, deserialize("$dir/saved_at.bin"))
    n.representation_network        = device(n, deserialize("$dir/representation.bin"))
    n.prediction_common_network     = device(n, deserialize("$dir/prediction_common.bin"))
    n.prediction_policy_network     = device(n, deserialize("$dir/prediction_policy.bin"))
    n.prediction_value_network      = device(n, deserialize("$dir/prediction_value.bin"))
    n.dynamics_hidden_state_network = device(n, deserialize("$dir/dynamics_hidden_state.bin"))
    n.dynamics_reward_network       = device(n, deserialize("$dir/dynamics_reward.bin"))
    enable_test_mode_if_needed!(n)
end

se(ŷ, y) = (ŷ .- y) .^ 2

function sum_se_scale_gradient(ŷ, y, n, scale)
    squared_error = se(ŷ, y)
    @assert size(squared_error) == size(scale)
    sum(scale_gradient.(squared_error, scale)) / Float32(n)
end

logitcrossentropy(ŷ, y; dims = 1) = .-sum(y .* logsoftmax(ŷ; dims = dims); dims = dims)
function sum_logitcrossentropy_scale_gradient(ŷ, y, n, scale)
    logit_cross_entropy = reshape(logitcrossentropy(ŷ, y), :)
    sum(scale_gradient.(logit_cross_entropy, scale)) / Float32(n)
end

function my_mse(ŷ, y, n, gradient_scale)
    sum(
        #scale_gradient.(
            (mse(ŷ, y; agg=collect) / Float32(n))
            #gradient_scale
        #)
    )
end
function my_crossentropy(ŷ, y, n, gradient_scale)
    sum(
        #scale_gradient.(
            (crossentropy(ŷ, y; agg=collect) / Float32(n))#,
            #gradient_scale
        #)
    )
end
#     scale_gradient.(
#         crossentropy(preds, targets; agg=collect) ./ Float32(size(targets,2)),
#         gradient_scale
#     ) |> sum
# end

function calc_target_loss(preds::Tuple, target, num_targets, gradient_scale=gpu(ones(Float32, 1, size(target, 2))))
    value_preds, reward_preds, policy_logits_preds, _ = preds

    target_values = view(target, 1:1, :)
    value_loss = sum_se_scale_gradient(value_preds, target_values, num_targets, gradient_scale)

    target_rewards = view(target, 2:2, :)
    reward_loss = sum_se_scale_gradient(reward_preds, target_rewards, num_targets, gradient_scale)

    target_policies = view(target, 3:size(target, 1), :)
    #policy_preds = softmax(policy_logits_preds)
    policy_loss = sum_logitcrossentropy_scale_gradient(policy_logits_preds, target_policies, num_targets, gradient_scale) #my_crossentropy

    #policy_preds = softmax(policy_logits_preds)
    #policy_loss = my_crossentropy(policy_preds, target_policies, num_targets, gradient_scale)

    TargetLoss(policy_loss, reward_loss, value_loss)
end

function calc_initial_loss(network::Incremental, batch::TrainingBatchTensor)
    num_targets = calc_num_initial_inference_targets(batch)
    net_output = initial_inference(network, images(batch))
    loss = calc_target_loss(net_output, image_targets(batch), num_targets)
    num_games = size(image_targets(batch), 2)
    (loss, net_output[4], 1:num_games)
end

function getloss(network::Incremental, batch::TrainingBatchTensor)
    #target_losses = calc_initial_loss(network, batch)[1]
    calc_recurrent_loss = recurrent_loss_fn(network, batch)
    target_loss, _, _ = reduce(calc_recurrent_loss, action_indices(batch); init=calc_initial_loss(network, batch))
    #regularization_loss = weight_decay * sum(l2_loss(p) for p in all_params)
    Loss(target_loss, Float32(0))
end

# @todo make action space size a config param
function recurrent_loss_fn(network::Incremental, batch::TrainingBatchTensor)
    action_targets = batch.action_targets
    actions = batch.actions
    num_targets = calc_num_recurrent_inference_targets(batch)

    return (prev_output, action_idx) -> begin

        prev_loss_sum, prev_hidden_states, prev_game_indices = prev_output

        action_idx_filter = action_targets[2,:] .== Float32(action_idx)

        game_indices = Zygote.dropgrad(view(action_targets, 1, action_idx_filter))

        remaining_games_filter = index_in(cpu(game_indices), cpu(prev_game_indices))
        remaining_hidden_states = view(prev_hidden_states, :, :, :, remaining_games_filter)

        @assert size(actions, 4) == size(action_targets, 2)
        next_actions = view(actions, :, :, :, action_idx_filter)

        # @todo make this more robust. Perhaps
        # put policy indices at the end of the vector.
        gradient_idx = size(action_targets, 1)
        curr_targets = view(action_targets, 4:(gradient_idx-1), action_idx_filter)

        net_output = recurrent_inference(network, remaining_hidden_states, next_actions)

        gradient_scale = view(action_targets, gradient_idx:gradient_idx, action_idx_filter)
        curr_loss = calc_target_loss(net_output, curr_targets, num_targets, gradient_scale)
        hidden_states = scale_gradient(net_output[4], Float32(0.5))

        (prev_loss_sum + curr_loss, hidden_states, game_indices)
    end
end


text:
module Mixtape

# Experimental re-write based upon:
# https://github.com/Keno/Compiler3.jl/blob/master/src/extracting_interpreter.jl

using LLVM
using LLVM.Interop
using LLVM_full_jll
using MacroTools: @capture,
                  postwalk,
                  rmlines,
                  unblock
using CodeInfoTools
using CodeInfoTools: resolve
using Core: MethodInstance,
            CodeInstance,
            CodeInfo
using Core.Compiler: WorldView,
                     NativeInterpreter,
                     InferenceResult,
                     coverage_enabled,
                     copy_exprargs,
                     convert_to_ircode,
                     slot2reg,
                     compact!,
                     ssa_inlining_pass!,
                     getfield_elim_pass!,
                     adce_pass!,
                     type_lift_pass!,
                     verify_ir,
                     verify_linetable

using GPUCompiler: cached_compilation,
                   FunctionSpec

import GPUCompiler: AbstractCompilerTarget,
                    NativeCompilerTarget,
                    AbstractCompilerParams,
                    CompilerJob,
                    julia_datalayout,
                    llvm_machine,
                    llvm_triple

import Core.Compiler: InferenceState,
                      InferenceParams,
                      AbstractInterpreter,
                      OptimizationParams,
                      InferenceState,
                      OptimizationState,
                      get_world_counter,
                      get_inference_cache,
                      lock_mi_inference,
                      unlock_mi_inference,
                      code_cache,
                      optimize,
                      may_optimize,
                      may_compress,
                      may_discard_trees,
                      add_remark!

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

export CompilationContext,
       NoContext,
       allow,
       transform,
       optimize!,
       OptimizationBundle,
       get_ir,
       get_state,
       julia_passes!,
       jit,
       emit,
       @ctx,
       @load_abi

include("context.jl")

#####
##### Utilities
#####

function get_methodinstance(@nospecialize(sig))
    ms = Base._methods_by_ftype(sig, 1, Base.get_world_counter())
    @assert length(ms) == 1
    m = ms[1]
    mi = ccall(:jl_specializations_get_linfo,
               Ref{MethodInstance}, (Any, Any, Any),
               m[3], m[1], m[2])
    return mi
end

function infer(interp, fn, t::Type{T}) where T <: Tuple
    mi = get_methodinstance(Tuple{typeof(fn), t.parameters...})
    src = Core.Compiler.typeinf_ext_toplevel(interp, mi)
    @assert(haskey(interp.code, mi))
    ci = getindex(interp.code, mi)
    if ci !== nothing && ci.inferred === nothing
        ci.inferred = src
    end
    return mi
end

#####
##### Interpreter
#####

# Based on: https://github.com/Keno/Compiler3.jl/blob/master/exploration/static.jl

# Holds its own cache.
struct StaticInterpreter{I, C} <: AbstractInterpreter
    code::Dict{MethodInstance, CodeInstance}
    inner::I
    messages::Vector{Tuple{MethodInstance, Int, String}}
    optimize::Bool
    ctx::C
end

function StaticInterpreter(; ctx = NoContext(), opt = false)
    return StaticInterpreter(Dict{MethodInstance, CodeInstance}(),
                             NativeInterpreter(),
                             Tuple{MethodInstance, Int, String}[],
                             opt,
                             ctx)
end

InferenceParams(si::StaticInterpreter) = InferenceParams(si.inner)
OptimizationParams(si::StaticInterpreter) = OptimizationParams(si.inner)
get_world_counter(si::StaticInterpreter) = get_world_counter(si.inner)
get_inference_cache(si::StaticInterpreter) = get_inference_cache(si.inner)
lock_mi_inference(si::StaticInterpreter, mi::MethodInstance) = nothing
unlock_mi_inference(si::StaticInterpreter, mi::MethodInstance) = nothing
code_cache(si::StaticInterpreter) = si.code
Core.Compiler.get(a::Dict, b, c) = Base.get(a, b, c)
Core.Compiler.get(a::WorldView{<:Dict}, b, c) = Base.get(a.cache,b,c)
Core.Compiler.haskey(a::Dict, b) = Base.haskey(a, b)
Core.Compiler.haskey(a::WorldView{<:Dict}, b) =
Core.Compiler.haskey(a.cache, b)
Core.Compiler.setindex!(a::Dict, b, c) = setindex!(a, b, c)
Core.Compiler.may_optimize(si::StaticInterpreter) = si.optimize
Core.Compiler.may_compress(si::StaticInterpreter) = false
Core.Compiler.may_discard_trees(si::StaticInterpreter) = false
function Core.Compiler.add_remark!(si::StaticInterpreter, sv::InferenceState, msg)
    push!(si.messages, (sv.linfo, sv.currpc, msg))
end

#####
##### Pre-inference
#####

function resolve_generic(a)
    if a <: Function && isdefined(a, :instance)
        return a.instance
    else
        return resolve(a)
    end
end

function custom_pass!(interp::StaticInterpreter, result::InferenceResult, mi::Core.MethodInstance, src)
    src === nothing && return src
    mi.specTypes isa UnionAll && return src
    sig = Tuple(mi.specTypes.parameters)
    as = map(resolve_generic, sig)
    if allow(interp.ctx, mi.def.module, as...)
        src = transform(interp.ctx, src, sig)
    end
    return src
end

function InferenceState(result::InferenceResult, cached::Bool, interp::StaticInterpreter)
    src = Core.Compiler.retrieve_code_info(result.linfo)
    mi = result.linfo
    src = custom_pass!(interp, result, mi, src)
    src === nothing && return nothing
    Core.Compiler.validate_code_in_debug_mode(result.linfo, src, "lowered")
    return InferenceState(result, src, cached, interp)
end

#####
##### Julia optimization pipeline
#####

@static if VERSION >= v"1.7.0-DEV.662"
    using Core.Compiler: finish as _finish
else
    function _finish(interp::AbstractInterpreter,
            opt::OptimizationState,
            params::OptimizationParams, ir, @nospecialize(result))
        return Core.Compiler.finish(opt, params, ir, result)
    end
end

struct OptimizationBundle
    ir::Core.Compiler.IRCode
    sv::OptimizationState
end
get_ir(b::OptimizationBundle) = b.ir
get_state(b::OptimizationBundle) = b.sv

@doc(
"""
    struct OptimizationBundle
        ir::Core.Compiler.IRCode
        sv::OptimizationState
    end
    get_ir(b::OptimizationBundle) = b.ir
    get_state(b::OptimizationBundle) = b.sv

Object which holds inferred `ir::Core.Compiler.IRCode` and a `Core.Compiler.OptimizationState`. Provided to the user through [`optimize!`](@ref), so that the user may plug in their own optimizations.
""", OptimizationBundle)

function julia_passes!(ir::Core.Compiler.IRCode, ci::CodeInfo,
        sv::OptimizationState)
    ir = compact!(ir)
    ir = ssa_inlining_pass!(ir, ir.linetable, sv.inlining, ci.propagate_inbounds)
    ir = compact!(ir)
    ir = getfield_elim_pass!(ir)
    ir = adce_pass!(ir)
    ir = type_lift_pass!(ir)
    ir = compact!(ir)
    return ir
end

julia_passes!(b::OptimizationBundle) = julia_passes!(b.ir, b.sv.src, b.sv)

function optimize(interp::StaticInterpreter, opt::OptimizationState,
        params::OptimizationParams, @nospecialize(result))
    nargs = Int(opt.nargs) - 1
    mi = opt.linfo
    meth = mi.def
    preserve_coverage = coverage_enabled(opt.mod)
    ir = convert_to_ircode(opt.src, copy_exprargs(opt.src.code), preserve_coverage, nargs, opt)
    ir = slot2reg(ir, opt.src, nargs, opt)
    b = OptimizationBundle(ir, opt)
    ir :: Core.Compiler.IRCode = optimize!(interp.ctx, b)
    verify_ir(ir)
    verify_linetable(ir.linetable)
    return _finish(interp, opt, params, ir, result)
end

#####
##### FunctionGraph
#####

abstract type FunctionGraph end

struct StaticSubGraph <: FunctionGraph
    code::Dict{MethodInstance, Any}
    instances::Vector{MethodInstance}
    entry::MethodInstance
end

entrypoint(ssg::StaticSubGraph) = ssg.entry

get_codeinstance(ssg::StaticSubGraph, mi::MethodInstance) = getindex(ssg.code, mi)

function get_codeinfo(code::Core.CodeInstance)
    ci = code.inferred
    if ci isa Vector{UInt8}
        return Core.Compiler._uncompressed_ir(code, ci)
    else
        return ci
    end
end

function get_codeinfo(graph::StaticSubGraph,
        cursor::MethodInstance)
    return get_codeinfo(get_codeinstance(graph, cursor))
end

function analyze(@nospecialize(f), tt::Type{T};
        ctx = NoContext(), opt = false) where T <: Tuple
    si = StaticInterpreter(; ctx = ctx, opt = opt)
    mi = infer(si, f, tt)
    si, StaticSubGraph(si.code, collect(keys(si.code)), mi)
end

#####
##### LLVM optimization pipeline
#####

function optimize!(tm, mod::LLVM.Module)
    ModulePassManager() do pm
        add_library_info!(pm, triple(mod))
        add_transform_info!(pm, tm)
        propagate_julia_addrsp!(pm)
        scoped_no_alias_aa!(pm)
        type_based_alias_analysis!(pm)
        basic_alias_analysis!(pm)
        cfgsimplification!(pm)
        scalar_repl_aggregates!(pm) # SSA variant?
        mem_cpy_opt!(pm)
        always_inliner!(pm)
        alloc_opt!(pm)
        instruction_combining!(pm)
        cfgsimplification!(pm)
        scalar_repl_aggregates!(pm) # SSA variant?
        instruction_combining!(pm)
        jump_threading!(pm)
        instruction_combining!(pm)
        reassociate!(pm)
        early_cse!(pm)
        alloc_opt!(pm)
        loop_idiom!(pm)
        loop_rotate!(pm)
        lower_simdloop!(pm)
        licm!(pm)
        loop_unswitch!(pm)
        instruction_combining!(pm)
        ind_var_simplify!(pm)
        loop_deletion!(pm)
        alloc_opt!(pm)
        scalar_repl_aggregates!(pm) # SSA variant?
        instruction_combining!(pm)
        gvn!(pm)
        mem_cpy_opt!(pm)
        sccp!(pm)
        instruction_combining!(pm)
        jump_threading!(pm)
        dead_store_elimination!(pm)
        alloc_opt!(pm)
        cfgsimplification!(pm)
        loop_idiom!(pm)
        loop_deletion!(pm)
        jump_threading!(pm)
        aggressive_dce!(pm)
        instruction_combining!(pm)
        barrier_noop!(pm)
        lower_exc_handlers!(pm)
        gc_invariant_verifier!(pm, false)
        late_lower_gc_frame!(pm)
        final_lower_gc!(pm)
        lower_ptls!(pm, #=dump_native=# false)
        cfgsimplification!(pm)
        instruction_combining!(pm) # Extra for Enzyme
        run!(pm, mod)
    end
end

#####
##### Entry codegen with cached compilation
#####

# Interpreter holds the cache.
function cache_lookup(si::StaticInterpreter, mi::MethodInstance,
        min_world, max_world)
    Base.get(si.code, mi, nothing)
end

# Mostly from GPUCompiler.
# In future, try to upstream any requires changes.
function codegen(job::CompilerJob)
    f = job.source.f
    tt = job.source.tt
    opt = job.params.opt
    si, ssg = analyze(f, tt;
                      ctx = job.params.ctx, opt = opt) # Populate local cache.
    world = get_world_counter(si)
    λ_lookup = (mi, min, max) -> cache_lookup(si, mi, min, max)
    lookup_cb = @cfunction($λ_lookup, Any, (Any, UInt, UInt))
    params = Base.CodegenParams(;
                                track_allocations = false,
                                code_coverage     = false,
                                prefer_specsig    = true,
                                gnu_pubnames      = false,
                                lookup            = Base.unsafe_convert(Ptr{Nothing}, lookup_cb))

    GC.@preserve lookup_cb begin
        native_code = ccall(:jl_create_native,
                            Ptr{Cvoid},
                            (Vector{MethodInstance},
                             Base.CodegenParams, Cint),
                            [ssg.entry],
                            params, 1) # = extern policy = #
        @assert native_code != C_NULL
        llvm_mod_ref = ccall(:jl_get_llvm_module,
                             LLVM.API.LLVMModuleRef,
                             (Ptr{Cvoid},),
                             native_code)
        @assert llvm_mod_ref != C_NULL
        llvm_mod = LLVM.Module(llvm_mod_ref)
    end
    code = cache_lookup(si, ssg.entry, world, world)
    llvm_func_idx = Ref{Int32}(-1)
    llvm_specfunc_idx = Ref{Int32}(-1)
    ccall(:jl_get_function_id,
          Nothing,
          (Ptr{Cvoid}, Any, Ptr{Int32}, Ptr{Int32}),
          native_code, code, llvm_func_idx, llvm_specfunc_idx)
    @assert llvm_specfunc_idx[] != -1
    @assert llvm_func_idx[] != -1
    llvm_func_ref = ccall(:jl_get_llvm_function,
                          LLVM.API.LLVMValueRef,
                          (Ptr{Cvoid}, UInt32),
                          native_code,
                          llvm_func_idx[] - 1)
    @assert llvm_func_ref != C_NULL
    llvm_func = LLVM.Function(llvm_func_ref)
    llvm_specfunc_ref = ccall(:jl_get_llvm_function,
                              LLVM.API.LLVMValueRef,
                              (Ptr{Cvoid}, UInt32),
                              native_code,
                              llvm_specfunc_idx[] - 1)
    @assert llvm_specfunc_ref != C_NULL
    llvm_specfunc = LLVM.Function(llvm_specfunc_ref)
    triple!(llvm_mod, llvm_triple(job.target))
    if julia_datalayout(job.target) !== nothing
        datalayout!(llvm_mod, julia_datalayout(job.target))
    end
    return (Any, llvm_specfunc, llvm_func, llvm_mod)
end

struct MixtapeCompilerTarget <: AbstractCompilerTarget end

llvm_triple(::MixtapeCompilerTarget) = Sys.MACHINE

function get_llvm_optlevel(opt_level::Int)
    if opt_level < 2
        optlevel = LLVM.API.LLVMCodeGenLevelNone
    elseif opt_level == 2
        optlevel = LLVM.API.LLVMCodeGenLevelDefault
    else
        optlevel = LLVM.API.LLVMCodeGenLevelAggressive
    end
    return optlevel
end

struct MixtapeCompilerParams <: AbstractCompilerParams
    opt::Bool
    optlevel::Int
    ctx::CompilationContext
end

function MixtapeCompilerParams(; opt = false,
        optlevel = Base.JLOptions().opt_level,
        ctx = NoContext())
    return MixtapeCompilerParams(opt, optlevel, ctx)
end

function llvm_machine(::MixtapeCompilerTarget,
        params::MixtapeCompilerParams)
    optlevel = get_llvm_optlevel(params.optlevel)
    tm = LLVM.JITTargetMachine(; optlevel=optlevel)
    LLVM.asm_verbosity!(tm, true)
    return tm
end

struct Entry{F, RT, TT}
    f::F
    specfunc::Ptr{Cvoid}
    func::Ptr{Cvoid}
end

@generated function (entry::Entry{F, RT, TT})(args...) where {F, RT, TT}
    # Slow ABI -- requires array allocation and unpacking. But stable.
    expr = quote
        args = Any[args...]
        ccall(entry.func,
              Any,
              (Any, Ptr{Any}, Int32),
              entry.f, args, length(args))
    end
    return expr
end

const jit_compiled_cache = Dict{UInt, Any}()

function jit(@nospecialize(f), tt::Type{T};
        ctx = NoContext(), opt = true,
        optlevel = Base.JLOptions().opt_level) where T <: Tuple
    fspec = FunctionSpec(f, tt, false, nothing) #=name=#
    job = CompilerJob(MixtapeCompilerTarget(),
                      fspec,
                      MixtapeCompilerParams(;
                                            opt = opt,
                                            ctx = ctx,
                                            optlevel = optlevel))
    return cached_compilation(jit_compiled_cache, job, _jit, _jitlink)
end

@doc(
"""
    jit(f::F, tt::Type{T}; ctx = NoContext(),
        opt = true,
        optlevel = Base.JLOptions().opt_level) where {F, T <: Type}

Compile and specialize a method instance for signature `Tuple{f, tt.parameters...}` with pipeline parametrized by `ctx::CompilationContext`.

Returns a callable instance of `Entry{F, RT, TT}` where `RT` is the return type of the instance after inference.

The user can configure the pipeline with optional arguments:

- `ctx::CompilationContext` -- configure [`transform`](@ref) and [`optimize!`](@ref).
- `opt::Bool` -- configure whether or not the Julia optimizer is run (including [`optimize!`](@ref)).
- `optlevel::Int > 0` -- configure the LLVM optimization level.
""", jit)

function _jitlink(job::CompilerJob, (rt, llvm_mod, func_name, specfunc_name))
    fspec = job.source
    tm = llvm_machine(job.target, job.params)
    orc = LLVM.OrcJIT(tm)
    atexit() do
        return LLVM.dispose(orc)
    end
    optimize!(tm, llvm_mod)
    jitted_mod = compile!(orc, llvm_mod)
    specfunc_addr = addressin(orc, jitted_mod, specfunc_name)
    specfunc_ptr = pointer(specfunc_addr)
    func_addr = addressin(orc, jitted_mod, func_name)
    func_ptr = pointer(func_addr)
    @assert(!(specfunc_ptr === C_NULL || func_ptr === C_NULL))
    return Entry{typeof(fspec.f), rt, fspec.tt}(fspec.f, specfunc_ptr, func_ptr)
end

function _jit(job::CompilerJob)
    rt, llvm_specfunc, llvm_func, llvm_mod = codegen(job)
    specfunc_name = LLVM.name(llvm_specfunc)
    func_name = LLVM.name(llvm_func)
    linkage!(llvm_func, LLVM.API.LLVMExternalLinkage)
    linkage!(llvm_specfunc, LLVM.API.LLVMExternalLinkage)
    return (rt, llvm_mod, func_name, specfunc_name)
end

#####
##### Emission and call interface
#####

function _emit(job::CompilerJob)
    f = job.source.f
    tt = job.source.tt
    opt = job.params.opt
    si, ssg = analyze(f, tt;
                      ctx = job.params.ctx, opt = opt) # Populate local cache.
    return get_codeinfo(ssg, entrypoint(ssg))
end

identity(job::CompilerJob, src) = src

const emit_compiled_cache = Dict{UInt, Any}()

function emit(@nospecialize(f), tt::Type{T};
        ctx = NoContext(), opt = false) where {F <: Function, T <: Tuple}
    fspec = FunctionSpec(f, tt, false, nothing) #=name=#
    optlevel = Base.JLOptions().opt_level
    job = CompilerJob(MixtapeCompilerTarget(),
                      fspec,
                      MixtapeCompilerParams(;
                                            ctx = ctx,
                                            opt = opt,
                                            optlevel = optlevel))
    return cached_compilation(emit_compiled_cache, job, _emit, identity)
end

@doc(
"""
    emit(@nospecialize(f), tt::Type{T};
        ctx = NoContext(), opt = false) where {F <: Function, T <: Tuple}

Emit typed (and optimized if `opt = true`) `CodeInfo` using the Mixtape pipeline. The user can configure the pipeline with optional arguments:

- `ctx::CompilationContext` -- configure [`transform`](@ref) and [`optimize!`](@ref).
- `opt::Bool` -- configure whether or not the Julia optimizer is run (including [`optimize!`](@ref)).
""", emit)

macro load_abi()
    expr = quote
        function cached_call(entry::Mixtape.Entry{F, RT, TT},
                args...) where {F, RT, TT}

            # TODO: Fast ABI.
            #ccall(entry.func, Any, (Any, $(nargs...), ), entry.f, $(_args...))

            # Slow ABI. Requires an array allocation.
            expr = quote
                vargs = Any[args...]
                ccall($(entry.func), Any, (Any, Ptr{Any}, Int32), $(entry.f), vargs,
                      $(length(args)))
            end
            return expr
        end

        @generated function _call(ctx::CompilationContext,
                optlevel::Val{T}, f::Function, args...) where T
            TT = Tuple{args...}
            entry = jit(f.instance, TT;
                        ctx = ctx(), opt = true, optlevel = T)
            return cached_call(entry, args...)
        end

        function call(f::T, args...;
                ctx = NoContext(),
                optlevel = Base.JLOptions().opt_level) where T <: Function
            _call(ctx, Val(optlevel), f, args...)
        end
    end
    return esc(expr)
end

@doc(
"""
    @load_abi()
    ...expands...
    call(f::T, args...; ctx = NoContext(),
        optlevel = Base.JLOptions().opt_level) where T <: Function

A macro which expands to define an ABI function `call` into the scope of the calling module. `call` wraps an `@generated` function which is called with signature argument types `Tuple{f <: Function, args...}`. The underlying `@generated` function creates a new instance of `ctx` (thus, a nullary constructor is an implicit requirement of your own subtypes of [`CompilationContext`](@ref) for usage with `call`) and calls [`jit`](@ref) -- it then caches a `ccall` which calls a function pointer to the [GPUCompiler](https://github.com/JuliaGPU/GPUCompiler.jl)-compiled LLVM module.

The `call` interface currently uses a slow ABI `ccall` -- which costs an array allocation for each toplevel `call`. This allocation is required to construct a `Vector{Any}` for the arguments and pass a pointer to it over the line, where the call unboxes each argument.
""", :(@Mixtape.load_abi))

end # module


text:
function (+){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftSum2(a,b)
    eftRound(hi, lo, rounding)
end

function (-){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftDiff2(a,b)
    eftRound(hi, lo, rounding)
end

function (*){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftProd2(a,b)
    eftRound(hi, lo, rounding)
end

function (/){T<:AbstractFloat, R<:RoundingMode}(a::T, b::T, rounding::R)
    hi,lo = eftDiv2Approx(a,b)
    eftRound(hi, lo, rounding)
end

function (square){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSquare(a)
    eftRound(hi, lo, rounding)
end

function (sqrt){T<:AbstractFloat, R<:RoundingMode}(a::T, rounding::R)
    hi,lo = eftSqrtApprox(a)
    eftRound(hi, lo, rounding)
end

for op in (:(+), :(-), :(*), :(/))
    @eval ($op){T<:Integer, R<:RoundingMode}(a::Rational{T}, b::Rational{T}, rounding::R) = ($op)(a,b)
end

if isdefined(Main,:FloatFloat)
FF=Main.FloatFloat.FF
for op in (:(+), :(-), :(*), :(/))
    @eval begin
        function ($op){T<:FF, R<:RoundingMode}(a::T, b::T, rounding::R)
            c = zero(T)
            with_rounding(Float64,rounding) do
                c = ($op)(a, b)
            end
            c
        end
    end
end
end
function (+){T<:Real, R<:RoundingMode}(a::T, b::T, rounding::R)
    c = zero(T)
    with_rounding(T,R) do
        c = a + b
    end
    c
end

function (-){T<:Real, R<:RoundingMode}(a::T, b::T, rounding::R)
    c = zero(T)
    with_rounding(T,R) do
        c = a - b
    end
    c
end

function (*){T<:Real, R<:RoundingMode}(a::T, b::T, rounding::R)
    c = zero(T)
    with_rounding(T,R) do
        c = a * b
    end
    c
end

function (/){T<:Real, R<:RoundingMode}(a::T, b::T, rounding::R)
    c = zero(T)
    with_rounding(T,R) do
        c = a / b
    end
    c
end

function (square){T<:Real, R<:RoundingMode}(a::T, rounding::R)
    c = zero(T)
    with_rounding(T,R) do
        c = a*a
    end
    c
end

function (sqrt){T<:Real, R<:RoundingMode}(a::T, rounding::R)
    c = zero(T)
    with_rounding(T,R) do
        c = sqrt(a)
    end
    c
end


text:
#####
##### Evaluating the mistake rate of an AlphaZero agent
##### on the standard benchmark from Pascal Pons
#####

const SOLVER_ERROR =
  "To use this script, you should install Pascal Pons' Connect Four Solver.\n" *
  "See games/connect-four/solver/README.md for instructions."

const SESSION_DIR = "sessions/connect-four"
const SAVE_FILE = "pons-benchmark-results.json"
const PLOT_FILE = "pons-benchmark-results.png"
const ITC_STRIDE = 5
const NUM_WORKERS = 128
const DEBUG_MODE = false # Launches a quick run on a tiny dataset to help debug

# A benchmark to evaluate connect-four agents is available at:
#   http://blog.gamesolver.org/solving-connect-four/02-test-protocol/
# It features 6000 positions along with their expected scores.

# Let `ne` the number of elapsed moves. A game stage is either
# `begin` (ne <= 14), `middle` (14 < ne <= 28) or `end` (ne > 28)
const STAGES = [:beginning, :middle, :end]

# Let `nr` the number of remaining moves. A difficulty level is either
# `easy` (nr < 14), `medium` (14 <= nr < 28) or `hard` (nr > 28).
const DIFFICULTIES = [:easy, :medium, :hard]

# Return either `nothing` or a `(stage, level)` tuple
function parse_test_filename(name)
  try
    m = match(r"^Test_L(\d)_R(\d)$", name)
    L = parse(Int, m.captures[1])
    R = parse(Int, m.captures[2])
    return (STAGES[L], DIFFICULTIES[R])
  catch e
    rethrow(e)
    return nothing
  end
end

struct Bench
  stage :: Symbol
  difficulty :: Symbol
  entries :: Vector{Tuple{String, Int}}
end

function load_benchmarks(dir)
  @assert isdir(dir)
  benchmarks = Bench[]
  files = readdir(dir)
  for bf in files
    meta = parse_test_filename(bf)
    f = joinpath(dir, bf)
    if !isnothing(meta) && isfile(f)
      stage, difficulty = meta
      entries = []
      for L in readlines(f)
        try
          L = split(L)
          push!(entries, (L[1], parse(Int, L[2])))
        catch e
          rethrow(e)
        end
      end
    push!(benchmarks, Bench(stage, difficulty, entries))
    end
  end
  rank(b) = (
    findfirst(==(b.difficulty), DIFFICULTIES),
    findfirst(==(b.stage), STAGES))
  sort!(benchmarks, by=rank)
  return benchmarks
end

const BENCHMARKS_DIR = joinpath(@__DIR__, "..", "benchmark")
const BENCHMARKS = load_benchmarks(BENCHMARKS_DIR)

####
#### Testing code
####

using AlphaZero
include("../main.jl")
using .ConnectFour: Game, Solver, Training
using ProgressMeter
using Formatting
using Statistics: mean

function state_of_string(str)
  g = Game()
  for c in str
    a = parse(Int, c)
    GI.play!(g, a)
  end
  return g
end

function optimal_on(solver, player, e)
  s = state_of_string(e[1])
  as, π = AlphaZero.think(player, s)
  api = argmax(π)
  qs = [Solver.qvalue(solver, s, a) for a in as]
  return sign(qs[api]) == maximum(sign, qs)
end

# TODO: the Solver is not thread safe.
function test_player_on(make_player, oracle, benchmark, progress)
  # Split benchmark entries into batches for workers
  entries = DEBUG_MODE ?
    Iterators.take(benchmark.entries, 2) :
    benchmark.entries
  batch_size = max(length(entries) ÷ NUM_WORKERS, 1)
  batches = Iterators.partition(entries, batch_size) |> collect
  nworkers = length(batches)
  # Spawn nworkers workers (nworkers <= NUM_WORKERS)
  spawn_oracle, done = AlphaZero.batchify_oracles(oracle, false, nworkers)
  solver = Solver.Player()
  results = AlphaZero.Util.threads_pmap(batches) do batch
    player = make_player(spawn_oracle())
    res = map(batch) do e
      err = optimal_on(solver, player, e) ? 0. : 1.
      next!(progress)
      return err
    end
    done()
    return res
  end
  return mean(Iterators.flatten(results) |> collect)
end

function test_player(make_player, oracle)
  errs = Float64[]
  for bench in BENCHMARKS
    p = Progress(length(bench.entries))
    err = test_player_on(make_player, oracle, bench, p)
    push!(errs, err)
    err_str =
    println("($(bench.stage), $(bench.difficulty)): $(fmt(".2f", 100 * err))%")
  end
  return errs
end

function load_alphazero_env(dir)
  return AlphaZero.UserInterface.load_env(
    Game, Training.Network{Game},
    AlphaZero.Log.Logger(devnull), dir, params=Training.params)
end

####
#### Running all the tests
####

const Errors = Vector{Float64} # One error rate for each benchmark

# Type gathering all produced data.
struct Results
  minmax :: Errors # Minmax baseline
  alphazero :: Errors #
  alphazero_training :: Vector{Tuple{Int, Errors}}
end

function test_alphazero(env)
  mcts_params = MctsParams(env.params.arena.mcts,
    temperature=ConstSchedule(0.),
    dirichlet_noise_ϵ=0.)
  net = Network.copy(env.bestnn, on_gpu=true, test_mode=true)
  return test_player(net) do net
    return MctsPlayer(net, mcts_params)
  end
end

function generate_data(session_dir)
  println("Testing the minmax baseline")
  minmax_errs = test_player(nothing) do _
     MinMax.Player{Game}(depth=5, amplify_rewards=true, τ=0)
  end
  println("")

  println("Testing AlphaZero")
  env = load_alphazero_env(session_dir)
  az_errs = test_alphazero(env)
  println("")

  println("Testing AlphaZero during training")
  az_training_errs = []
  itc = 0
  while true
    itdir = joinpath(session_dir, "iterations", "$itc")
    AlphaZero.UserInterface.valid_session_dir(itdir) || break
    DEBUG_MODE && itc > 10 && break
    println("Iteration $itc")
    env = load_alphazero_env(itdir)
    push!(az_training_errs, (itc, test_alphazero(env)))
    itc += ITC_STRIDE
  end
  return Results(minmax_errs, az_errs, az_training_errs)
end

#####
##### Main
#####

if !isdir(Solver.DEFAULT_SOLVER_DIR)
  println(stderr, SOLVER_ERROR)
  exit()
end

using JSON2
using JSON3

JSON3.StructType(::Type{Results}) = JSON3.Struct()

# Regenerate data only if needed as it can take a while.
if isfile(SAVE_FILE)
  data = open(SAVE_FILE, "r") do io
    JSON3.read(io, Results)
  end
else
  data = generate_data(SESSION_DIR)
  open(SAVE_FILE, "w") do io
    JSON2.pretty(io, JSON3.write(data))
  end
end

import Plots

# i is the benchmark number
function sub_plot(benchs, results, n)
  bench = benchs[n]
  sname(s) = uppercase(string(s)[1]) * string(s)[2:end]
  title = "$(sname(bench.stage)) - $(sname(bench.difficulty))"
  azt = results.alphazero_training
  xs = [itc for (itc, _) in azt]
  ys = [100 * errs[n] for (itc, errs) in azt]
  y_minmax = 100 * results.minmax[n]
  #ymax = max(5, maximum(ys))
  ymax = max(maximum(ys), y_minmax)
  plt = Plots.plot(
    xs, ys, ylims=(0, ymax),
    # xlabel="Training Iteration",
    ylabel="Error Rate (in %)",
    legend=nothing,
    title=title,
    titlefontsize=10,
    xguidefontsize=8,
    yguidefontsize=8)
  Plots.hline!(plt, [y_minmax])
  return plt
end

function plot_results(results)
  plts = [sub_plot(BENCHMARKS, results, n) for n in eachindex(BENCHMARKS)]
  return Plots.plot(plts..., layout=(3, 2))
end

plot = plot_results(data)
Plots.savefig(plot, PLOT_FILE)


text:
module PerfTest
using BenchmarkTools, InteractiveUtils
using LinearAlgebra, Profile, Random, Test, WiSER

@info "generate data"
Random.seed!(123)
# dimensions
m  = 6000 # number of individuals
ns = rand(5:11, m) # numbers of observations per individual
p  = 5    # number of fixed effects, including intercept
q  = 3    # number of random effects, including intercept
l  = 5    # number of WS variance covariates, including intercept
obsvec = Vector{WSVarLmmObs{Float64}}(undef, m)
# true parameter values
βtrue = [0.1; 6.5; -3.5; 1.0; 5]
τtrue = [-1.5; 1.5; -0.5; zeros(l - 3)]
Σγ    = Matrix(Diagonal([2.0; 1.2; rand(q - 2)])) # full rank case
δγω   = [0.2; 0.1; rand(q - 2) ./ 10]
σω    = [1.0]
Σγω   = [Σγ δγω; δγω' σω]
Lγω   = cholesky(Symmetric(Σγω), check = false).L
Lγ    = Lγω[1:q, 1:q]
lγω   = Lγω[q + 1, 1:q]
lω    = Lγω[q + 1, q + 1]
# generate data
γω = Vector{Float64}(undef, q + 1)
z  = similar(γω) # hold vector of iid std normal
for i in 1:m
    # first column intercept, remaining entries iid std normal
    X = Matrix{Float64}(undef, ns[i], p)
    X[:, 1] .= 1
    @views randn!(X[:, 2:p])
    # first column intercept, remaining entries iid std normal
    Z = Matrix{Float64}(undef, ns[i], q)
    Z[:, 1] .= 1
    @views randn!(Z[:, 2:q])
    # first column intercept, remaining entries iid std normal
    W = Matrix{Float64}(undef, ns[i], l)
    W[:, 1] .= 1
    @views randn!(W[:, 2:l])
    # generate random effects: γω = Lγω * z
    mul!(γω, Lγω, randn!(z))
    # generate y
    μy = X * βtrue + Z * γω[1:q]
    @views ysd = exp.(0.5 .* (W * τtrue .+ dot(γω[1:q], lγω) .+ γω[end]))
    y = ysd .* randn(ns[i]) .+ μy
    # form a VarLmmObs instance
    obsvec[i] = WSVarLmmObs(y, X, Z, W)
end
# form VarLmmModel
vlmm = WSVarLmmModel(obsvec);

@testset "init_ls!" begin
# least squares starting point
init_ls!(vlmm)
@show vlmm.β
@show vlmm.τ
@show vlmm.Lγ
# @info "type stability"
# @code_warntype init_ls!(vlmm)
@info "benchmark"
bm = @benchmark init_ls!($vlmm, gniters=5)
display(bm); println()
@test_skip allocs(bm) == 0
# @info "profile"
# Profile.clear()
# @profile @btime init_ls!(vlmm; gniters=5)
# Profile.print(format=:flat)
end

@testset "update_wtmat!" begin
# set parameter values to be the truth
copy!(vlmm.β, βtrue)
copy!(vlmm.τ, τtrue)
vlmm.τ[1] = τtrue[1] + 0.5(abs2(lω) + abs2(norm(lγω + Lγ'lγω)))
vlmm.Lγ  .= Lγ
@show vlmm.β
@show vlmm.τ
@show vlmm.Lγ
# update weight matrix
update_wtmat!(vlmm)
@show vlmm.β
@show vlmm.∇β
# @show vlmm.Hββ
# @info "type stability"
# @code_warntype update_wtmat!(vlmm)
@info "benchmark"
bm = @benchmark update_wtmat!($vlmm)
display(bm); println()
@test allocs(bm) == 0
# @info "profile"
# Profile.clear()
# @profile @btime update_wtmat!(vlmm)
# Profile.print(format=:flat)
end

@testset "nlsv_obj! (unweighted)" begin
# set parameter values to be the truth
copy!(vlmm.β, βtrue)
copy!(vlmm.τ, τtrue)
vlmm.τ[1] = τtrue[1] + 0.5(abs2(lω) + abs2(norm(lγω + Lγ'lγω)))
vlmm.Lγ  .= Lγ
@show vlmm.β
@show vlmm.τ
@show vlmm.Lγ
# evaluate objective (at truth)
@info "obj/grad/hessian at true parameter values"
vlmm.iswtnls[1] = false
@show nlsv_obj!(vlmm, true, true, true)
@show vlmm.∇β
@show vlmm.∇τ
@show vlmm.∇Lγ
H = [vlmm.Hττ vlmm.HτLγ; vlmm.HτLγ' vlmm.HLγLγ]
# display(H); println()
@test norm(H - transpose(H)) / norm(H) < 1e-8
@test all(eigvals(Symmetric(H)) .≥ 0)
# @info "type stability"
# @code_warntype nlsv_obj!(vlmm.data[1], vlmm.β, vlmm.τ, vlmm.Lγ, true)
# @code_warntype nlsv_obj!(vlmm, true)
@info "benchmark"
# bm = @benchmark nlsv_obj!($vlmm.data[1], $vlmm.β, $vlmm.τ, $vlmm.Lγ, true)
# display(bm)
# @test allocs(bm) == 0
bm = @benchmark nlsv_obj!($vlmm, true, true, true)
display(bm); println()
@test allocs(bm) == 0
# @info "profile"
# Profile.clear()
# @profile @btime nlsv_obj!($vlmm, true, true)
# Profile.print(format=:flat)
end

@testset "nlsv_obj! (weighted)" begin
# set parameter values to be the truth
copy!(vlmm.β, βtrue)
copy!(vlmm.τ, τtrue)
vlmm.τ[1] = τtrue[1] + 0.5(abs2(lω) + abs2(norm(lγω + Lγ'lγω)))
vlmm.Lγ  .= Lγ
@show vlmm.β
@show vlmm.τ
@show vlmm.Lγ
# update_wtmat then evaluate at the truth
update_wtmat!(vlmm)
vlmm.iswtnls[1] = true
@show nlsv_obj!(vlmm, true, true, true)
@show vlmm.∇β
@show vlmm.∇τ
@show vlmm.∇Lγ
H = [vlmm.Hττ vlmm.HτLγ; vlmm.HτLγ' vlmm.HLγLγ]
# display(H); println()
@test norm(H - transpose(H)) / norm(H) < 1e-8
@test all(eigvals(Symmetric(H)) .≥ 0)
# @info "type stability"
# @code_warntype nlsv_obj!(vlmm.data[1], vlmm.β, vlmm.τ, vlmm.Lγ, true)
# @code_warntype nlsv_obj!(vlmm, true)
@info "benchmark"
bm = @benchmark nlsv_obj!($vlmm, true, true, true)
display(bm); println()
@test allocs(bm) == 0
# @info "profile"
# Profile.clear()
# @profile @btime nlsv_obj!($vlmm, true, true, true)
# Profile.print(format=:flat)
end

@testset "nlsv_obj! (weighted) - parallel" begin
# set parameter values to be the truth
copy!(vlmm.β, βtrue)
copy!(vlmm.τ, τtrue)
vlmm.τ[1] = τtrue[1] + 0.5(abs2(lω) + abs2(norm(lγω + Lγ'lγω)))
vlmm.Lγ  .= Lγ
@show vlmm.β
@show vlmm.τ
@show vlmm.Lγ
# update_wtmat then evaluate at the truth
update_wtmat!(vlmm)
vlmm.iswtnls[1] = true
vlmm.ismthrd[1] = true
@show nlsv_obj!(vlmm, true, true, true)
@show vlmm.∇β
@show vlmm.∇τ
@show vlmm.∇Lγ
H = [vlmm.Hττ vlmm.HτLγ; vlmm.HτLγ' vlmm.HLγLγ]
# display(H); println()
@test norm(H - transpose(H)) / norm(H) < 1e-8
@test all(eigvals(Symmetric(H)) .≥ 0)
# @info "type stability"
# @code_warntype nlsv_obj!(vlmm.data[1], vlmm.β, vlmm.τ, vlmm.Lγ, true)
# @code_warntype nlsv_obj!(vlmm, true)
@info "benchmark"
vlmm.iswtnls[1] = true
vlmm.ismthrd[1] = true
bm = @benchmark nlsv_obj!($vlmm, true, true, true)
display(bm); println()
# @info "profile"
# Profile.clear()
# @profile @btime nlsv_obj!($vlmm, true, true, true)
# Profile.print(format=:flat)
end
end


text:
using DifferentialEquations
f(u,p,t) = 0.98u
u0 = 1.0
tspan = (0.0,1.0)
prob = ODEProblem(f,u0,tspan)

sol = solve(prob)

using Plots; gr()
plot(sol)

plot(sol,linewidth=5,title="Solution to the linear ODE with a thick line",
     xaxis="Time (t)",yaxis="u(t) (in μm)",label="My Thick Line!") ## legend=false

plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

sol.t

sol.u

[t+u for (u,t) in tuples(sol)]

sol

sol(0.45)

sol = solve(prob,abstol=1e-8,reltol=1e-8)

plot(sol)
plot!(sol.t, t->1.0*exp(0.98t),lw=3,ls=:dash,label="True Solution!")

sol = solve(prob,saveat=0.1)

sol = solve(prob,saveat=[0.2,0.7,0.9])

sol = solve(prob,saveat=[0.2,0.7,0.9],save_start = false, save_end = false)

sol = solve(prob,dense=false)

sol = solve(prob,save_everystep=false)

sol = solve(prob,alg_hints=[:stiff])

sol = solve(prob,Tsit5(),reltol=1e-6)

function lorenz!(du,u,p,t)
    σ,ρ,β = p
    du[1] = σ*(u[2]-u[1])
    du[2] = u[1]*(ρ-u[3]) - u[2]
    du[3] = u[1]*u[2] - β*u[3]
end

u0 = [1.0,0.0,0.0]

p = (10,28,8/3) ## we could also make this an array, or any other type!

tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan,p)

sol = solve(prob)

sol.t[10],sol[10]

sol[2,10]

A = convert(Array,sol)

plot(sol)

plot(sol,vars=(1,2,3))

plot(sol,vars=(1,2,3),denseplot=false)

plot(sol,vars=(0,2))

function lotka_volterra!(du,u,p,t)
  du[1] = p[1]*u[1] - p[2]*u[1]*u[2]
  du[2] = -p[3]*u[2] + p[4]*u[1]*u[2]
end

lv! = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d

u0 = [1.0,1.0]
p = (1.5,1.0,3.0,1.0)
tspan = (0.0,10.0)
prob = ODEProblem(lv!,u0,tspan,p)
sol = solve(prob)
plot(sol)

lv!.Jex

A  = [1. 0  0 -5
      4 -2  4 -3
     -4  0  0  1
      5 -2  2  3]
u0 = rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

sol[3]

big_u0 = big.(u0)

prob = ODEProblem(f,big_u0,tspan)
sol = solve(prob)

sol[1,3]

prob = ODEProblem(f,big_u0,big.(tspan))
sol = solve(prob)

using StaticArrays
A  = @SMatrix [ 1.0  0.0 0.0 -5.0
                4.0 -2.0 4.0 -3.0
               -4.0  0.0 0.0  1.0
                5.0 -2.0 2.0  3.0]
u0 = @SMatrix rand(4,2)
tspan = (0.0,1.0)
f(u,p,t) = A*u
prob = ODEProblem(f,u0,tspan)
sol = solve(prob)

sol[3]

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
@testset "node" begin
    @testset "constructors" begin
        n = AVLTrees.Node(10, 10)
        @test n.key == 10
        @test n.data == 10
        @test isnothing(n.parent)

        n1 = AVLTrees.Node(12, 12, n)
        @test n1.parent == n
        @test isnothing(n1.parent.parent)
    end
end


text:
using Example
if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# @test hello("Julia") == "Hello, Julia"
# @test domath(2.0) ≈ 7.0

### SDM
K = 5
α = ones(K)
β = 100.0
d = SparseDirMix(α, β)
rand(d)

### SBM
a, b = shape_Dir2GenDir(α)

nsim = 10e3
xD = hcat([rDirichlet(α) for i = 1:nsim ]...)
xGD = hcat([rGenDirichlet(a,b)[1] for i = 1:nsim ]...)

using Statistics
sum(abs.(mean(xD, dims=2) .- mean(xGD, dims=2)))

η = 1.0e3
π_small = 0.5
π_large = 0.1
γ = 1.5
δ = 1.5
p = SBMprior(K, η, π_small, π_large, γ, δ)

rand(p)

x = [0, 0, 2, 1, 0]

pp = SBM_multinom_post(p, x)

rand(pp)


text:
using ImageHistogram
@static if VERSION < v"0.7.0-DEV.2005"
    using Base.Test
else
    using Test
end

# just a simple test to pass
@test (1+1) == 2


text:
round3(x) = round(x, digits=3)
round2(x) = round(x, digits=2)

max_stars_repo_name: lihebi/DAG-EQ

max_issues_repo_name: lihebi/DAG-EQ

max_forks_repo_name: lihebi/DAG-EQ


text:
#!/usr/bin/julia
using Pkg.Artifacts
using Pkg.PlatformEngines
using Pkg.GitTools
using SHA

NAME = "jars"
URL = "https://github.com/JuliaSpace/orekit-artifacts/archive/10.2.tar.gz"

function sha256sum(tarball_path)
    return open(tarball_path, "r") do io
        return bytes2hex(sha256(io))
    end
end

function add_artifact!(
    artifacts_toml::String,
    name::String,
    tarball_url::String;
    clear=true,
    options...,
)
    probe_platform_engines!()

    tarball_path = download(tarball_url)
    sha256 = sha256sum(tarball_path)

    git_tree_sha1 = create_artifact() do artifact_dir
        unpack(tarball_path, artifact_dir)
    end

    rm(tarball_path)
    clear && remove_artifact(git_tree_sha1)

    bind_artifact!(
        artifacts_toml,
        name,
        git_tree_sha1;
        download_info = [(tarball_url, sha256)],
        options...,
    )

    return git_tree_sha1
end

add_artifact!(joinpath(@__DIR__, "..", "Artifacts.toml"), NAME, URL, force=true)


text:
MOIU.@model(
    DualizableModel,
    (),
    (MOI.EqualTo, MOI.GreaterThan, MOI.LessThan),
    (
        MOI.Reals,
        MOI.Zeros,
        MOI.Nonnegatives,
        MOI.Nonpositives,
        MOI.SecondOrderCone,
        MOI.RotatedSecondOrderCone,
        MOI.ExponentialCone,
        MOI.DualExponentialCone,
        MOI.PositiveSemidefiniteConeTriangle,
    ),
    (MOI.PowerCone, MOI.DualPowerCone),
    (),
    (MOI.ScalarAffineFunction,),
    (MOI.VectorOfVariables,),
    (MOI.VectorAffineFunction,)
)

mutable struct PrimalDualMap{T}
    constrained_var_idx::Dict{VI,Tuple{CI,Int}}
    constrained_var_dual::Dict{CI,CI}
    constrained_var_zero::Dict{
        CI,
        Union{MOI.VectorAffineFunction{T},MOI.ScalarAffineFunction{T}},
    }
    primal_var_dual_con::Dict{VI,CI}
    primal_con_dual_var::Dict{CI,Vector{VI}}
    primal_con_dual_con::Dict{CI,CI}
    primal_con_constants::Dict{CI,Vector{T}}

    primal_parameter::Dict{VI,VI} # in the future we could have Parameters, as in ParameterJuMP
    primal_var_dual_quad_slack::Dict{VI,VI}

    function PrimalDualMap{T}() where {T}
        return new(
            Dict{VI,Tuple{CI,Int}}(),
            Dict{CI,CI}(),
            Dict{
                CI,
                Union{MOI.VectorAffineFunction{T},MOI.ScalarAffineFunction{T}},
            }(),
            Dict{VI,CI}(),
            Dict{CI,Vector{VI}}(),
            Dict{CI,CI}(),
            Dict{CI,Vector{T}}(),
            Dict{VI,VI}(),
            Dict{VI,VI}(),
        )
    end
end

function is_empty(primal_dual_map::PrimalDualMap{T}) where {T}
    return isempty(primal_dual_map.constrained_var_idx) &&
           isempty(primal_dual_map.constrained_var_dual) &&
           isempty(primal_dual_map.constrained_var_zero) &&
           isempty(primal_dual_map.primal_var_dual_con) &&
           isempty(primal_dual_map.primal_con_dual_var) &&
           isempty(primal_dual_map.primal_con_dual_con) &&
           isempty(primal_dual_map.primal_con_constants) &&
           isempty(primal_dual_map.primal_parameter) &&
           isempty(primal_dual_map.primal_var_dual_quad_slack)
end

function empty!(primal_dual_map::PrimalDualMap)
    Base.empty!(primal_dual_map.constrained_var_idx)
    Base.empty!(primal_dual_map.constrained_var_dual)
    Base.empty!(primal_dual_map.constrained_var_zero)
    Base.empty!(primal_dual_map.primal_var_dual_con)
    Base.empty!(primal_dual_map.primal_con_dual_var)
    Base.empty!(primal_dual_map.primal_con_dual_con)
    Base.empty!(primal_dual_map.primal_con_constants)
    Base.empty!(primal_dual_map.primal_parameter)
    Base.empty!(primal_dual_map.primal_var_dual_quad_slack)
    return
end

struct DualProblem{T,OT<:MOI.ModelLike}
    dual_model::OT #It can be a model from an optimizer or a DualizableModel{T}
    primal_dual_map::PrimalDualMap{T}

    function DualProblem{T}(
        dual_optimizer::OT,
        pdmap::PrimalDualMap{T},
    ) where {T,OT<:MOI.ModelLike}
        return new{T,OT}(dual_optimizer, pdmap)
    end
end

function DualProblem{T}(dual_optimizer::OT) where {T,OT<:MOI.ModelLike}
    return DualProblem{T}(dual_optimizer, PrimalDualMap{T}())
end

function DualProblem(dual_optimizer::OT) where {OT<:MOI.ModelLike}
    return DualProblem{Float64}(dual_optimizer)
end

# Empty DualProblem cosntructor
function DualProblem{T}() where {T}
    return DualProblem{T}(DualizableModel{T}(), PrimalDualMap{T}())
end