text:
module Periods

import Base: -, *, /, +, get, isapprox, show, zero, eltype, (:), isless

export TimeUnit, Second, Minute, Hour, Day, Year, Century,
    seconds, minutes, hours, days, years, centuries,
    Period, -, *, /, +, value, unit,
    SECONDS_PER_MINUTE,
    SECONDS_PER_HOUR,
    SECONDS_PER_DAY,
    SECONDS_PER_YEAR,
    SECONDS_PER_CENTURY,
    MINUTES_PER_HOUR,
    MINUTES_PER_DAY,
    MINUTES_PER_YEAR,
    MINUTES_PER_CENTURY,
    HOURS_PER_DAY,
    HOURS_PER_YEAR,
    HOURS_PER_CENTURY,
    DAYS_PER_YEAR,
    DAYS_PER_CENTURY,
    YEARS_PER_CENTURY

const SECONDS_PER_MINUTE   = 60.0
const SECONDS_PER_HOUR     = 60.0 * 60.0
const SECONDS_PER_DAY      = 60.0 * 60.0 * 24.0
const SECONDS_PER_YEAR     = 60.0 * 60.0 * 24.0 * 365.25
const SECONDS_PER_CENTURY  = 60.0 * 60.0 * 24.0 * 365.25 * 100.0

const MINUTES_PER_HOUR     = 60.0
const MINUTES_PER_DAY      = 60.0 * 24.0
const MINUTES_PER_YEAR     = 60.0 * 24.0 * 365.25
const MINUTES_PER_CENTURY  = 60.0 * 24.0 * 365.25 * 100.0

const HOURS_PER_DAY        = 24.0
const HOURS_PER_YEAR       = 24.0 * 365.25
const HOURS_PER_CENTURY    = 24.0 * 365.25 * 100.0

const DAYS_PER_YEAR        = 365.25
const DAYS_PER_CENTURY     = 365.25 * 100.0

const YEARS_PER_CENTURY    = 100.0

"""
All time units are subtypes of the abstract type `TimeUnit`.
The following time units are defined:

* `Second`
* `Minute`
* `Hour`
* `Day`
* `Year`
* `Century`
"""
abstract type TimeUnit end

struct Second <: TimeUnit end
struct Minute <: TimeUnit end
struct Hour <: TimeUnit end
struct Day <: TimeUnit end
struct Year <: TimeUnit end
struct Century <: TimeUnit end

const seconds = Second()
const minutes = Minute()
const hours = Hour()
const days = Day()
const years = Year()
const centuries = Century()

Base.broadcastable(u::TimeUnit) = Ref(u)

"""
    Period{U, T}(unit, Δt) where {U<:TimeUnit, T}

A `Period` object represents a time interval of `Δt` with a [`TimeUnit`](@ref) of `unit`.
Periods should be constructed via the shorthand syntax shown in the examples below.

### Examples ###

```jldoctest
julia> 3.0seconds
3.0 seconds

julia> 1.0minutes
1.0 minutes

julia> 12hours
12 hours

julia> days_per_year = 365
365
julia> days_per_year * days
365 days

julia> 10.0years
10.0 years

julia> 1centuries
1 century
```
"""
struct Period{U<:TimeUnit, T}
    unit::U
    Δt::T
end

value(p::Period) = p.Δt
unit(p::Period) = p.unit
zero(p::Period) = Period(unit(p), zero(value(p)))
zero(p::Type{<:Period{U}}) where {U} = Period(U(), 0.0)
zero(p::Type{<:Period{U,T}}) where {U, T} = Period(U(), zero(T))
eltype(p::Period) = typeof(value(p))
eltype(p::Type{<:Period{U,T}}) where {U, T} = T

name(::Second, val::Integer) = ifelse(val == one(val), "second", "seconds")
name(::Second, ::Any) = "seconds"
name(::Minute, val::Integer) = ifelse(val == one(val), "minute", "minutes")
name(::Minute, ::Any) = "minutes"
name(::Hour, val::Integer) = ifelse(val == one(val), "hour", "hours")
name(::Hour, ::Any) = "hours"
name(::Day, val::Integer) = ifelse(val == one(val), "day", "days")
name(::Day, ::Any) = "days"
name(::Year, val::Integer) = ifelse(val == one(val), "year", "years")
name(::Year, ::Any) = "years"
name(::Century, val::Integer) = ifelse(val == one(val), "century", "centuries")
name(::Century, ::Any) = "centuries"

function show(io::IO, p::Period)
    u = unit(p)
    v = value(p)
    print(io, v, " ", name(u, v))
end

*(Δt::T, unit::TimeUnit) where {T<:Number} = Period(unit, Δt)
*(unit::TimeUnit, Δt::T) where {T<:Number} = Period(unit, Δt)
*(A::TimeUnit, B::AbstractArray) = broadcast(*, A, B)
*(A::AbstractArray, B::TimeUnit) = broadcast(*, A, B)

(+)(p1::Period{U}, p2::Period{U}) where {U}= Period(unit(p1), p1.Δt + p2.Δt)
(-)(p1::Period{U}, p2::Period{U}) where {U}= Period(unit(p1), p1.Δt - p2.Δt)
(-)(p::Period) = Period(unit(p), -p.Δt)
(*)(x, p::Period) = Period(unit(p), p.Δt * x)
(*)(p::Period, x) = Period(unit(p), p.Δt * x)
(/)(p::Period, x) = Period(unit(p), p.Δt / x)

isless(p1::Period{U}, p2::Period{U}) where {U} = isless(value(p1), value(p2))
isapprox(p1::Period{U}, p2::Period{U}) where {U} = value(p1) ≈ value(p2)

(:)(start::Period{U,T}, stop::Period{U,T}) where {U,T} = (:)(start, one(T) * U(), stop)

function (:)(start::Period{U}, step::Period{U}, stop::Period{U}) where {U}
    step = start < stop ? step : -step
    StepRangeLen(start, step, floor(Int, value(stop-start)/value(step))+1)
end

Period{U,T}(p::Period{U,T}) where {U,T} = p

Base.step(r::StepRangeLen{<:Period}) = r.step

(::Second)(p::Period{Second})   = p
(::Second)(p::Period{Minute})   = Period(seconds, p.Δt * SECONDS_PER_MINUTE)
(::Second)(p::Period{Hour})     = Period(seconds, p.Δt * SECONDS_PER_HOUR)
(::Second)(p::Period{Day})      = Period(seconds, p.Δt * SECONDS_PER_DAY)
(::Second)(p::Period{Year})     = Period(seconds, p.Δt * SECONDS_PER_YEAR)
(::Second)(p::Period{Century})  = Period(seconds, p.Δt * SECONDS_PER_CENTURY)

(::Minute)(p::Period{Second})   = Period(minutes, p.Δt / SECONDS_PER_MINUTE)
(::Minute)(p::Period{Minute})   = p
(::Minute)(p::Period{Hour})     = Period(minutes, p.Δt * MINUTES_PER_HOUR)
(::Minute)(p::Period{Day})      = Period(minutes, p.Δt * MINUTES_PER_DAY)
(::Minute)(p::Period{Year})     = Period(minutes, p.Δt * MINUTES_PER_YEAR)
(::Minute)(p::Period{Century})  = Period(minutes, p.Δt * MINUTES_PER_CENTURY)

(::Hour)(p::Period{Second})     = Period(hours, p.Δt / SECONDS_PER_HOUR)
(::Hour)(p::Period{Minute})     = Period(hours, p.Δt / MINUTES_PER_HOUR)
(::Hour)(p::Period{Hour})       = p
(::Hour)(p::Period{Day})        = Period(hours, p.Δt * HOURS_PER_DAY)
(::Hour)(p::Period{Year})       = Period(hours, p.Δt * HOURS_PER_YEAR)
(::Hour)(p::Period{Century})    = Period(hours, p.Δt * HOURS_PER_CENTURY)

(::Day)(p::Period{Second})      = Period(days, p.Δt / SECONDS_PER_DAY)
(::Day)(p::Period{Minute})      = Period(days, p.Δt / MINUTES_PER_DAY)
(::Day)(p::Period{Hour})        = Period(days, p.Δt / HOURS_PER_DAY)
(::Day)(p::Period{Day})         = p
(::Day)(p::Period{Year})        = Period(days, p.Δt * DAYS_PER_YEAR)
(::Day)(p::Period{Century})     = Period(days, p.Δt * DAYS_PER_CENTURY)

(::Year)(p::Period{Second})     = Period(years, p.Δt / SECONDS_PER_YEAR)
(::Year)(p::Period{Minute})     = Period(years, p.Δt / MINUTES_PER_YEAR)
(::Year)(p::Period{Hour})       = Period(years, p.Δt / HOURS_PER_YEAR)
(::Year)(p::Period{Day})        = Period(years, p.Δt / DAYS_PER_YEAR)
(::Year)(p::Period{Year})       = p
(::Year)(p::Period{Century})    = Period(years, p.Δt * YEARS_PER_CENTURY)

(::Century)(p::Period{Second})  = Period(centuries, p.Δt / SECONDS_PER_CENTURY)
(::Century)(p::Period{Minute})  = Period(centuries, p.Δt / MINUTES_PER_CENTURY)
(::Century)(p::Period{Hour})    = Period(centuries, p.Δt / HOURS_PER_CENTURY)
(::Century)(p::Period{Day})     = Period(centuries, p.Δt / DAYS_PER_CENTURY)
(::Century)(p::Period{Year})    = Period(centuries, p.Δt / YEARS_PER_CENTURY)
(::Century)(p::Period{Century}) = p

end


text:
# This implements Adaptive Generalized Approximate Message Parsing (adGAMP)
# for compressed sensing / regularized regression as described in:
#   Kamilov et al, "Approximate Message Passing with Consistent Parameter
#   Estimation and Applications to Sparse Learning", arxiv.org, 2012,
#   http://arxiv.org/pdf/1207.3859v3.pdf
#

# Unpack the adGamp parameters and options.
#
# Output:
#  Tuple of parameters
#
# Input:
#  A - measurement matrix (M*N) where M is the number of cases and N the number of covariates
#  y - noisy measurements (M*1) assuming y = A * x + err
#
# Options:
#  x0 - starting point for search (if used in sequential scenarios etc)
#  maxIterations -
#  maxParamIterations -
function unpack_adGamp_params(A, y;
  adapt = true,
  maxIterations      = 100,
  maxParamIterations = 300,
  x0 = nothing,       # A starting point, will be set to zero if not given

  f::Function = scad, # Defining function used in tresholding, by default SCAD. Should be soft?

  tolerance = 1e-6,   # Convergence level, i.e. terminate algorithm if norm change is less than this between iterations

  )

  M, N = size(A)
  xhat = x0 || 1

end

  a = 1

immutable AdGampOptions
  adapt::Bool
  maxIterations::Int64
  maxParamIterations::Int64
end

function adaptiveGamp(A, y; kws...)
  A, y, denoiseFunc, estimatingFuncs, xhat, vx, theta, options = unpack_adGamp_params(A, y; kws...)

  m, n = size(A)
  trA = A'
  hadamardA = A .* A
  trhadamardA = hadamardA'
  xhatprev = 1e6 * ones(n, 1)
  shat = zeros(m, 1)
  numBelowTolerance = 0

  for i in 1:maxIterations
    vp = hadamardA * vx
    phat = A * xhat .- (vp .* shat)
    vs = 1 ./ (vp + noiseVariance)
    shat = vs .* (y .- phat)

    # Update variables
    vr = 1 ./ (trhadamardA * vs)
    rhat = xhat .+ (vr .* (trA * shat))

    if options.adapt
      theta = update_distributional_params(rhat, vr, theta, estimatingFuncs, options)
    end

    xhat, vx = denoiseFunc(rhat, vr, theta)

    # If the change in estimated vector is less than the tolerance level
    # we count it as a violation. We only allow a limited number of such
    # violations.
    if(norm(xhat - xhatprev, 1) < options.tolerance)
      numBelowTolerance += 1
      if numBelowTolerance > options.maxToleranceViolations
        break
      end
    end

  end

end

ext: jl


text:
using PiADS1x15
using Test

@testset "PiADS1x15.jl" begin
    # Write your own tests here.
end


text:
# generates the table that defines interval boundaries / heights and
# the various constants required for Chopin's algorithm to sample from
# a truncated normal.

using Roots, Distributions

# returns residual area above N intervals of size 'area' starting at x = 0,
# works for d symmetric around 0 with descreasing pdf above x > 0
function arearesidual(d::ContinuousUnivariateDistribution, area::Float64, N::Int)
    x = 0.0
    for i = 1:N
        x += area / pdf(d, x)  # s.t. area = dx * pdf(min(xi, xi+1))
    end
    area - cdf(d, -x)   # cdf(d, x) = 1 - cdf(d, -x) for symmetric d
end

# returns list of interval boundaries for given number N of intervals above 0
# and the lowest interval of interest, xmin
function genxtable(N::Int, xmin::Real=-2.0)
    stdn = Normal(0, 1)
    # find area such that all intervals cover the same area
    area = fzero(x -> arearesidual(stdn, x, N), 0.5 / N, 0.7 / N)
    # find interval boundaries for x >= 0
    x = Array(Float64, N+1)
    x[1] = 0.0
    for i = 1:N
        xi = x[i]
        x[i+1] = xi + area / pdf(stdn, xi)
    end
    # do not store boundaries below amin
    j = N + 1
    while x[j] > -xmin
        j -= 1
    end  # j points at interval that contains amin
    [-x[j+1:-1:2], x]
end

# returns the other interval tables given the list of interval boundaries. Here,
# yl is the smaller of the two pdf(x)'s' at the interval boundaries
# yu is the larger of the two pdf(x)'s' at the interval boundaries
# d is the interval width
# delta is d yu / yl for each interval
function genydtables(x::Vector{Float64})
    stdn = Normal(0, 1)
    nd = length(x) - 1  # number of intervals [x[i], x[i+1]]
    d = Array(Float64, nd)
    delta = Array(Float64, nd)
    yl = Array(Float64, nd)
    yu = Array(Float64, nd)
    for i = 1:nd
        xi, xi1 = x[i], x[i+1]
        if xi < 0
            yl[i] = pdf(stdn, xi)
            yu[i] = pdf(stdn, xi1)
        else
            yl[i] = pdf(stdn, xi1)
            yu[i] = pdf(stdn, xi)
        end
        d[i] = xi1 - xi
        delta[i] = d[i] * yu[i] / yl[i]
    end
    yl, yu, d, delta
end

# returns the required constants, xmin, xmax, i0, invh
function genconsts(x, d)
    h = minimum(d)    # minimum interval range
    xmin, xmax = x[1], x[end]
    i0 = 1 - ifloor(xmin / h)
    xmin, xmax, i0, 1/h
end

# returns list of indices such that, for some xmin <= a < xmax,
#   k = ncell[i0 + ifloor(a * invh)]
# guarantees
#   x[k] <= a < x[k+2]
function genncell(x, xmin, xmax, i0, invh)
    nb = length(x) - 1
    p = 1 + ifloor(xmax * invh) - ifloor(xmin * invh)
    ncell = Array(Uint16, p)
    k = 1
    for i = 1:p  # ith interval [h(i-i0), h(i-i0+1))
        while k < nb && x[k] < (i - i0) / invh  # ensure x[k] just above h(i-i0)
            k += 1
        end
        ncell[i] = max(k - 1, 1)
    end
    ncell
end

# write vector of numbers in Julia-valid format to io
function writevector{T <: Number}(io::IO, x::Vector{T}, xname, fmtstr::String, xperline::Int)
    n = length(x)
    @eval dofmt(x) = @sprintf($fmtstr, x)
    println(io, "const ", xname, " = ", string(T), "[")
    for i = 1:xperline:n
        print(io, "    ")
        for j = i:min(n, i+xperline-1)
            print(io, dofmt(x[j]), j == n ? "]" : ", ")
        end
        println(io)
    end
end
writevector(io::IO, x::Vector{Float64}, xname, xperline::Int=3) = writevector(io, x, xname, "%.16e", xperline)
writevector(io::IO, x::Vector{Uint16}, xname, xperline::Int=10) = writevector(io, x, xname, "%5d", xperline)

# writes full table to io
function writetable(io::IO, N::Int, xmin::Real=-2.0)
    xmin_ini = xmin
    x = genxtable(N, xmin)
    yl, yu, d, delta = genydtables(x)
    xmin, xmax, i0, invh = genconsts(x, d)
    ncell = genncell(x, xmin, xmax, i0, invh)

    println(io, "# value tables for sampling from a truncated normal distribution")
    println(io, "# generated by https://github.com/jdrugo/fasttn/src/gentables.jl")
    println(io, "# called with N=", N, ", xmin=", xmin_ini)
    println(io, "# number of intervals: ", length(x) - 1)
    println(io)

    # constants
    println(io, "# constants")
    println(io, "const tnzig_xmin = ", xmin)
    println(io, "const tnzig_xmax = ", xmax)
    println(io, "const tnzig_i0 = ", i0)
    println(io, "const tnzig_invh = ", invh)
    println(io)

    # tables
    println(io, "# inverval tables")
    writevector(io, x, "tnzig_x")
    println(io)
    writevector(io, yl, "tnzig_yl")
    println(io)
    writevector(io, yu, "tnzig_yu")
    println(io)
    writevector(io, d, "tnzig_d")
    println(io)
    writevector(io, delta, "tnzig_delta")
    println(io)
    writevector(io, ncell, "tnzig_ncell")
    println(io)
end

# create tables
f = open("randtn_zigtables.jl", "w")
writetable(f, 2048)
close(f)

max_stars_repo_name: jdrugo/fasttn

max_issues_repo_name: jdrugo/fasttn

max_forks_repo_name: jdrugo/fasttn


text:
@testset "1041.robot-bounded-in-circle.jl" begin
    @test is_robot_bounded("GGLLGG") == true
    @test is_robot_bounded("GG") == false
    @test is_robot_bounded("GL") == true
end


text:
## ----------------  Individual files

# Convenience abbreviation
rf(ds :: DataSettings, baseFn, groups, selfTrans, area, momentType) =
    RawDataFile(selfTrans, area, momentType,
        file_name(ds, baseFn, groups), ds);

# AFQT percentile by quality
raw_afqt_pct_qual(ds :: DataSettings; momentType :: Symbol = :mean) =
    RawDataFile(:selfReport, :freshmen, momentType,
        file_name(ds, "afqt_pctile", :qual), ds);

raw_frac_local_qual(ds :: DataSettings; momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :freshmen, momentType,
        file_name(ds, "frac_loc_50", [:qual]), ds);

# ----------  By quality / gpa

rf_qg(ds :: DataSettings, baseFn, selfTrans, area, momentType) =
    rf(ds, baseFn, [:qual, :afqt], selfTrans, area, momentType);

# Conditional on entry, fraction of students in each [quality, gpa] cell
raw_entry_qual_gpa(ds :: DataSettings; momentType :: Symbol = :mean) =
    rf_qg(ds, "jointdist", :transcript, :freshmen, momentType);
        # file_name(ds, "jointdist", [:qual, :afqt]), ds);

# By year
raw_frac_drop_qual_gpa(ds :: DataSettings, year :: Integer;
    momentType :: Symbol = :mean) =
    rf_qg(ds, "drop_rate_y$year", :transcript, :progress, momentType);
        # file_name(ds, "drop_rate_y$year", [:qual, :afqt]), ds);

raw_grad_rate_qual_gpa(ds :: DataSettings; momentType :: Symbol = :mean) =
    rf_qg(ds, "grad_rate", :transcript, :progress, momentType)
        # file_name(ds, "grad_rate", [:quality, :afqt]), ds);

raw_study_time_qual_gpa(ds :: DataSettings; momentType :: Symbol = :mean) =
    rf_qg(ds, "studytime", :selfReport, :progress, momentType);
        # file_name(ds, "studytime", [:quality, :afqt]), ds);

raw_class_time_qual_gpa(ds :: DataSettings; momentType :: Symbol = :mean) =
    rf_qg(ds, "classtime", :selfReport, :progress, momentType);
    #    file_name(ds, "classtime", [:quality, :afqt]), ds);

# Transcript time to drop contains a 0 in one cell
raw_time_to_drop_4y_qual_gpa(ds :: DataSettings;
    momentType :: Symbol = :mean) =
    rf_qg(ds, "time_to_drop_4Y", :transcript, :progress, momentType);
        # file_name(ds, "time_to_drop_4Y", [:qual, :afqt]), ds);

raw_time_to_grad_4y_qual_gpa(ds :: DataSettings;
    momentType :: Symbol = :mean) =
    rf_qg(ds, "time_to_grad_4Y", :transcript, :progress, momentType);
        # file_name(ds, "time_to_grad_4Y", [:qual, :afqt]), ds);

raw_work_hours_qual_gpa(ds :: DataSettings, year :: Integer;
    momentType :: Symbol = :mean) =
    RawDataFile(:selfReport, :finance, momentType,
        file_name(ds, "hours_y$year", [:qual, :afqt]), ds);

raw_net_price_qual_gpa(ds :: DataSettings, year :: Integer;
    momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :finance, momentType,
        file_name(ds, "net_price_y$year", [:qual, :afqt]), ds);

raw_credits_taken_qual_gpa(ds :: DataSettings, year :: Integer;
    momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :progress, momentType,
        file_name(ds, "creds_att_y$year", [:qual, :afqt]), ds);

# Can load for selected percentiles giving `percentile` input (e.g. 90)
function raw_cum_loans_qual_year(
    ds :: DataSettings, year :: Integer;
    momentType :: Symbol = :mean, percentile = nothing
    )
    fn = file_name(ds, "cumloans_y$year", [:qual, :afqt];
        percentile = percentile);
    return RawDataFile(:transcript, :finance, momentType, fn, ds);
end

# Conditional on entry, fraction of students in each [quality, gpa] cell
raw_mass_entry_qual_gpa(ds :: DataSettings; momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :freshmen, momentType,
        file_name(ds, "jointdist", [:qual, :afqt]), ds);

# -------  By quality / grad outcome

raw_credits_taken_qual_grad_year(ds :: DataSettings, year :: Integer;
    momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :progress, momentType,
        file_name(ds, "creds_att_y$year", [:qual, :gradDrop]), ds);

# -----  By quality / parental

rf_qp(ds :: DataSettings, baseFn, selfTrans, area, momentType) =
    rf(ds, baseFn, [:qual, :yp], selfTrans, area, momentType);

raw_work_hours_qual_parental(ds :: DataSettings;
    momentType :: Symbol = :mean) =
    rf_qp(ds, "hours_y1", :selfReport, :finance, momentType);

# Conditional on entry, fraction of students in each [quality, parental] cell
raw_entry_qual_parental(ds :: DataSettings; momentType :: Symbol = :mean) =
    rf_qp(ds, "jointdist", :transcript, :freshmen, momentType);

raw_frac_grad_qual_parental(ds :: DataSettings;
    momentType :: Symbol = :mean) =
    rf_qp(ds, "grad_rate", :transcript, :progress, momentType);

raw_time_to_drop_4y_qual_parental(ds :: DataSettings;
    momentType :: Symbol = :mean) =
    rf_qp(ds, "time_to_drop_4Y", :transcript, :progress, momentType);

raw_time_to_grad_4y_qual_parental(ds :: DataSettings;
    momentType :: Symbol = :mean) =
    rf_qp(ds, "time_to_grad_4Y", :transcript, :progress, momentType);

# ------  by [gpa, parental]
# Data files have this transposed as [parental, gpa]

# Entry rates [gpa, parental], but TRANSPOSED in data files!
raw_entry_gpa_parental(ds :: DataSettings; momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :hsGrads, momentType,
        file_name(ds, "entrants", [:yp, :afqt]), ds);

# Fraction by quality, by [gpa, parental], TRANSPOSED in data files.
# Conditional on entry.
raw_qual_entry_gpa_parental(ds :: DataSettings, iCollege;
    momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :freshmen, momentType,
        file_name(ds, "prob_ent_q$(iCollege)", [:yp, :afqt]), ds);

# Mass of HSG by [parental, hs gpa]
raw_mass_gpa_parental(ds :: DataSettings; momentType :: Symbol = :mean) =
    RawDataFile(:transcript, :hsGrads, momentType,
        file_name(ds, "jointdist_hsgrads", [:yp, :afqt]), ds);

# -----  Regressions

# Wage regression for graduates with college quality coefficients
# raw_wage_regr_grads(ds :: DataSettings) =
#     RawDataFile(:selfReport, :none, :regression, "loginc_reg2.dat", ds);

raw_transfer_regr(ds :: DataSettings; momentType = nothing) =
    RawDataFile(:transcript, :none, :regression, "parental_transfers1_reg1.dat", ds);

raw_tuition_regr(ds :: DataSettings; momentType = nothing) =
    RawDataFile(:transcript, :none, :regression, "net_price1_reg1.dat", ds);

# ---------


text:
@testset "AbsoluteValue2MILPBridge: $(fct_type), $(T)" for fct_type in ["vector of variables", "vector affine function"], dim in [2, 3], T in [Int, Float64]
    mock = MOIU.MockOptimizer(MILPModel{T}())
    model = COIB.AbsoluteValue2MILP{T}(mock)

    if T == Int
        @test MOI.supports_constraint(model, MOI.VariableIndex, MOI.Integer)
    end
    @test MOI.supports_constraint(
        model,
        MOI.ScalarAffineFunction{T},
        MOI.LessThan{T},
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorAffineFunction{T},
        CP.AbsoluteValue,
    )
    @test MOIB.supports_bridging_constraint(
        model,
        MOI.VectorOfVariables,
        CP.AbsoluteValue,
    )

    if T == Int
        x, _ = MOI.add_constrained_variable(model, MOI.Integer())
        x_abs, _ = MOI.add_constrained_variable(model, MOI.Integer())
    elseif T == Float64
        x = MOI.add_variable(model)
        x_abs = MOI.add_variable(model)
    end

    fct = if fct_type == "vector of variables"
        MOI.VectorOfVariables([x_abs, x])
    elseif fct_type == "vector affine function"
        MOIU.vectorize([x_abs, x])
    else
        @assert false
    end

    @test_throws AssertionError MOI.add_constraint(model, fct, CP.AbsoluteValue())

    MOI.add_constraint(model, x, MOI.LessThan(5 * one(T)))
    MOI.add_constraint(model, x, MOI.GreaterThan(2 * one(T)))

    c = MOI.add_constraint(model, fct, CP.AbsoluteValue())

    @test MOI.is_valid(model, x)
    @test MOI.is_valid(model, x_abs)
    @test MOI.is_valid(model, c)

    bridge = MOIBC.bridges(model)[MOI.ConstraintIndex{MOI.VectorOfVariables, CP.AbsoluteValue}(-1)]

    @testset "Bridge properties" begin
        @test MOIBC.concrete_bridge_type(typeof(bridge), MOI.VectorOfVariables, CP.AbsoluteValue) == typeof(bridge)
        @test MOIB.added_constrained_variable_types(typeof(bridge)) == [(MOI.ZeroOne), (MOI.GreaterThan{T})]
        @test MOIB.added_constraint_types(typeof(bridge)) == [
            (MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}),
            (MOI.ScalarAffineFunction{T}, MOI.LessThan{T}),
        ]

        @test MOI.get(bridge, MOI.NumberOfVariables()) == 3
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.ZeroOne}()) == 1
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.VariableIndex, MOI.GreaterThan{T}}()) == 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == 2
        @test MOI.get(bridge, MOI.NumberOfConstraints{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == 2

        @test MOI.get(bridge, MOI.ListOfVariableIndices()) == [bridge.var_bin, bridge.var_pos, bridge.var_neg]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.ZeroOne}()) == [bridge.var_bin_con]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.VariableIndex, MOI.GreaterThan{T}}()) == [bridge.var_pos_con, bridge.var_neg_con]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.EqualTo{T}}()) == [bridge.con_original_var, bridge.con_abs_var]
        @test MOI.get(bridge, MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{T}, MOI.LessThan{T}}()) == [bridge.con_pos_var_big_m, bridge.con_neg_var_big_m]
    end

    @testset "Set of variables" begin
        @test MOI.is_valid(model, bridge.var_bin)
        @test MOI.is_valid(model, bridge.var_pos)
        @test MOI.is_valid(model, bridge.var_neg)
    end

    @testset "Set of constraints" begin
        @test MOI.is_valid(model, bridge.var_bin_con)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.var_bin_con)
        @test f == bridge.var_bin
        @test MOI.get(model, MOI.ConstraintSet(), bridge.var_bin_con) == MOI.ZeroOne()

        @test MOI.is_valid(model, bridge.var_pos_con)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.var_pos_con)
        @test f == bridge.var_pos
        @test MOI.get(model, MOI.ConstraintSet(), bridge.var_pos_con) == MOI.GreaterThan(zero(T))

        @test MOI.is_valid(model, bridge.var_neg_con)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.var_neg_con)
        @test f == bridge.var_neg
        @test MOI.get(model, MOI.ConstraintSet(), bridge.var_neg_con) == MOI.GreaterThan(zero(T))

        @test MOI.is_valid(model, bridge.con_original_var)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_original_var)
        @test f.constant === zero(T)
        @test length(f.terms) == 3
        t1 = f.terms[1]
        @test t1.coefficient === one(T)
        @test t1.variable == x
        t2 = f.terms[2]
        @test t2.coefficient === -one(T)
        @test t2.variable == bridge.var_pos
        t3 = f.terms[3]
        @test t3.coefficient === one(T)
        @test t3.variable == bridge.var_neg
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_original_var) == MOI.EqualTo(zero(T))

        @test MOI.is_valid(model, bridge.con_abs_var)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_abs_var)
        @test f.constant === zero(T)
        @test length(f.terms) == 3
        t1 = f.terms[1]
        @test t1.coefficient === one(T)
        @test t1.variable == x_abs
        t2 = f.terms[2]
        @test t2.coefficient === -one(T)
        @test t2.variable == bridge.var_pos
        t3 = f.terms[3]
        @test t3.coefficient === -one(T)
        @test t3.variable == bridge.var_neg
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_abs_var) == MOI.EqualTo(zero(T))

        @test MOI.is_valid(model, bridge.con_pos_var_big_m)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_pos_var_big_m)
        @test abs(f.constant) === zero(T)
        @test length(f.terms) == 2
        t1 = f.terms[1]
        @test t1.coefficient === -5 * one(T)
        @test t1.variable == bridge.var_bin
        t2 = f.terms[2]
        @test t2.coefficient === one(T)
        @test t2.variable == bridge.var_pos
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_pos_var_big_m) == MOI.LessThan(zero(T))

        @test MOI.is_valid(model, bridge.con_neg_var_big_m)
        f = MOI.get(model, MOI.ConstraintFunction(), bridge.con_neg_var_big_m)
        @test abs(f.constant) === 5 * one(T)
        @test length(f.terms) == 2
        t1 = f.terms[1]
        @test t1.coefficient === 5 * one(T)
        @test t1.variable == bridge.var_bin
        t2 = f.terms[2]
        @test t2.coefficient === one(T)
        @test t2.variable == bridge.var_neg
        @test MOI.get(model, MOI.ConstraintSet(), bridge.con_neg_var_big_m) == MOI.LessThan(zero(T))
    end
end


text:
#=
Afft.jl
2019-07-07, Jeff Fessler, University of Michigan
=#

export Afft

# using MIRT: embed
using LinearMaps: LinearMap
using FFTW: fft, ifft
using Test: @test

"""
`A = Afft(samp::AbstractArray{Bool}; T::DataType = ComplexF32)`
Make a LinearMap object for (under-sampled) FFT, of type `T`,
using given sampling pattern `samp`.
Especially for compressed sensing MRI with Cartesian sampling.
"""
function Afft(samp::AbstractArray{Bool}; T::DataType = ComplexF32)
	dim = size(samp)
	return LinearMap{T}(
		x -> fft(reshape(x,dim))[samp],
		y -> prod(dim) * ifft(embed(y,samp)),
		sum(samp), prod(dim))
end

"""
`A = Afft(:test)`
self test
"""
function Afft(test::Symbol)
	test != :test && throw(ArgumentError("test $test"))
	samp = trues(3,2); samp[2] = false
	A = Afft(samp)
	@test isapprox(Matrix(A)', Matrix(A'))
	@test A * ones(3,2)[:] == [6, 0, 0, 0, 0]
	@test A' * [1, 0, 0, 0, 0] == ones(3,2)[:]
	true
end

# Afft(:test)


text:
using CUDA.CUDNN: cudnnPoolingMode_t, CUDNN_POOLING_MAX,
                  CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING,
                  cudnnPoolingForward!, pooldims, cudnnPoolingBackward

import CUDA.CUDNN: cudnnPoolingDescriptor

function cudnnPoolingDescriptor(pdims::PoolDims, x::DenseCuArray{T}, mode::cudnnPoolingMode_t) where T
    pdims, x = fix1d(pdims), fix1d(x)
    window, padding, stride = NNlib.kernel_size(pdims), nnlibPadding(pdims), NNlib.stride(pdims)
    nanOpt = CUDNN_NOT_PROPAGATE_NAN
    cudnnPoolingDescriptor(mode, nanOpt, Cint(max(2,ndims(x)-2)), pooldims(window,size(x)), pooldims(padding,size(x)), pooldims(stride,size(x)))
end

function maxpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_MAX)
    cudnnPoolingForward!(y, x, d)
end

function ∇maxpool!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    xDesc, yDesc = cudnnTensorDescriptor.((x, y))
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_MAX)
    alpha, beta = scalingParameter(T,1), scalingParameter(T,0)
    cudnnPoolingBackward(handle(), d, alpha, yDesc, y, yDesc, dy, xDesc, x, beta, xDesc, dx)
    return dx
end

function meanpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING)
    cudnnPoolingForward!(y, x, d)
end

function ∇meanpool!(dx::DenseCuArray{T}, dy::DenseCuArray{T}, y::DenseCuArray{T}, x::DenseCuArray{T}, pdims::PoolDims) where T<:CUDNNFloat
    xDesc, yDesc = cudnnTensorDescriptor.((x, y))
    d = cudnnPoolingDescriptor(pdims, x, CUDNN_POOLING_AVERAGE_COUNT_INCLUDE_PADDING)
    alpha, beta = scalingParameter(T,1), scalingParameter(T,0)
    cudnnPoolingBackward(handle(), d, alpha, yDesc, y, yDesc, dy, xDesc, x, beta, xDesc, dx)
    return dx
end

function maxpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, k; pad=map(_->0,k), stride=k) where {T<:CUDNNFloat}
    pdims = PoolDims(x, k; padding=pad, stride=stride)
    return maxpool!(y, x, pdims)
end

function meanpool!(y::DenseCuArray{T}, x::DenseCuArray{T}, k; pad=map(_->0,k), stride=k) where {T<:CUDNNFloat}
    pdims = PoolDims(x, k; padding=pad, stride=stride)
    return meanpool!(y, x, pdims)
end

# Deprecated methods
function cudnnPoolingForward(y::DenseCuArray{T,N}, x::DenseCuArray{T,N}, pdims::NNlib.PoolDims;
                             alpha=1, beta=0, mode=CUDNN_POOLING_MAX) where {T,N}
    # @warn "`cudnnPoolingForward(y,x,d::PoolDims)` is deprecated, please use one of the methods in `@doc cudnnPoolingForward`." maxlog=1
    cudnnPoolingForward!(y, x; window=NNlib.kernel_size(pdims), padding=nnlibPadding(pdims), stride=NNlib.stride(pdims), mode, alpha, beta)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule ConnectFourSolver_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("ConnectFourSolver")
JLLWrappers.@generate_main_file("ConnectFourSolver", UUID("5ab91c8a-5c16-538d-97cb-84915cefd4b0"))
end  # module ConnectFourSolver_jll


text:
"""
    CDFT(Yc, Xc, Xp)

"""
struct CDFT <: ContinuousUnivariateDistribution
    Yc::ContinuousUnivariateDistribution  # Distribution des données observées pour la période de calibration
    Xc::ContinuousUnivariateDistribution  # Distribution des données du modèle pour la période de calibration
    Xp::ContinuousUnivariateDistribution  # Distribution des données du modèle pour la période projetée
end

"""
    T_transform(m::CDFT, x::Real)

"""
function T_transform(m::CDFT, x::Real)

    F_Xp = cdf(m.Xp, x)  # Évaluation en x de la CDF des données du modèle pour la période projeté

    Fm1_Xc = quantile(m.Xc, F_Xp)  # Quantile correspondant de la distribution des données du modèle pour la période de calibration

    F_Yp = cdf(m.Yc,  Fm1_Xc)  # Évaluation de la CDF des données observées pour la période de calibration au point correspondant

    return F_Yp  # Valeur de la CDF évaluée au point x
end

#### Evaluation

cdf(m::CDFT, x::Real) = T_transform(m, x)

function pdf(m::CDFT, x::Real)
    # très approximatif... à revoir
    # Inf en zéro (?)

    δ = 0.00001
    x₁ = x - δ
    x₂ = x + δ

    return (cdf(m, x₂) - cdf(m, x₁)) / (x₂ - x₁)
end

function quantile(m::CDFT,
                  p::Real;
                  interval::Array{Float64,1}=[0.0, 500.0])

    @assert zero(p)<p<one(p) "the quantile level should be between 0 and 1."

    res = optimize(θ -> (p - cdf(m, θ))^2, interval[1], interval[2])

    return Optim.minimizer(res)[1]
end

#### Sampling

function rand(rng::AbstractRNG, m::CDFT)
    # Generate a Float64 random number uniformly in (0,1].
    u = 1 - rand(rng)

    return quantile(m, u)
end

sampler(m::CDFT) = m    # à vérifier


text:
mutable struct Marks
    n_vert::Int64
    m::Int64
    get_vertex::Dict{Int64,Int64}
end

function Marks(v::Int64,m::Int64)::Marks

    return Marks(v, m , Dict{Int64, Int64}())
end

function Base.iterate(mark::Marks, c::Int64 = 0)::Union{Nothing, Tuple{Marks, Int64}}

    if mark.m == 0 && c == 1
        return nothing
    end

    if isempty(mark.get_vertex)
        mark.get_vertex = Dict([i for i in 1:mark.m].=> 1)
        return mark, 1
    end

    index = findfirst(x -> mark.get_vertex[x] < mark.n_vert, 1:mark.m)
    if index === nothing
        return nothing
    end

    next_mark!(mark, index)
    return mark, 1
end

function next_mark!(mark::Marks, index::Int64)::Nothing

    mark.get_vertex[index] += 1
    for i in 1:(index-1)
        mark.get_vertex[i] = 1
    end

    return nothing
end
#=
"""
Take the marks of a graph and return a dictionary Dict{Int64,Vector{Int64}} which associate to each vertex the set of marks on the vertex.
"""=#
function invert_marks(mark::Marks)::Dict{Int64,Vector{Int64}}
    new_D = Dict{Int64,Vector{Int64}}()
    D = mark.get_vertex
    for k in keys(D)
        D[k] in keys(new_D) ? push!(new_D[D[k]],k) : new_D[D[k]] = [k]
    end
    for v in 1:mark.n_vert
        if !(v in keys(new_D))
            new_D[v]=Int64[]
        end
    end
    return new_D
end
#=
"""
Return the number of marks assigned to some vertex v.
# Arguments
- `mark::Marks`: the marks assigned to a graph.
- `v::Int64`: the vertex.
"""=#
function num_marks(mark::Marks, v::Int64)::Int64

    mark.m == 0 && return 0

    return count(t -> mark.get_vertex[t] == v, 1:mark.m)
end


text:
function get_double_value(obj::KeyValue)
    return jcall(obj, "getDoubleValue", jdouble, ())
end

function get_integer_value(obj::KeyValue)
    return jcall(obj, "getIntegerValue", jint, ())
end

function get_key(obj::KeyValue)
    return jcall(obj, "getKey", JString, ())
end

function get_keyword(obj::KeyValue)
    return jcall(obj, "getKeyword", Keyword, ())
end

function get_value(obj::KeyValue)
    return jcall(obj, "getValue", JString, ())
end


text:
export FastDownwardPlanner

"Wrapper to the FastDownward planning system."
@kwdef mutable struct FastDownwardPlanner <: Planner
    search::String = "astar"
    heuristic::String = "add"
    h_params::Dict{String, String} = Dict()
    timeout::Float64 = 10
    verbose::Bool = false
end

"Calls the FastDownward planning system to produce a plan."
function solve(planner::FastDownwardPlanner,
               domain::Domain, state::State, spec::Specification)
    if !haskey(ENV, "FD_PATH")
        error("FD_PATH not set to location of fast_downward.py")
    end
    @unpack search, heuristic, h_params, timeout, verbose = planner
    # Write temporary domain and problem files
    goals = get_goal_terms(spec)
    metric = hasfield(typeof(spec), :metric) ? spec.metric : nothing
    if metric != nothing metric = (-1, metric) end
    problem = Problem(state, Compound(:and, goals), metric; domain=domain.name)
    domain_path = save_domain(tempname(), domain)
    problem_path = save_problem(tempname(), problem)
    # Set up shell command to fast_downward.py
    fd_path = ENV["FD_PATH"]
    py_cmd = get(ENV, "PYTHON", "python")
    h_params = join(["$key=$val" for (key, val) in h_params], ", ")
    search_params = "$search($heuristic($h_params))"
    cmd = `$py_cmd $fd_path $domain_path $problem_path --search $search_params`
    # Run command up to timeout
    out = Pipe()
    proc = run(pipeline(cmd, stdout=out); wait=false)
    cb() = process_exited(proc)
    timedwait(cb, float(timeout))
    if process_running(proc)
        @debug "Planner timed out."
        Base.Filesystem.rm(domain_path); Base.Filesystem.rm(problem_path)
        kill(proc); close(out.in)
        return NullSolution()
    end
    # Read output and check if solution was found
    close(out.in)
    output = read(out, String)
    Base.Filesystem.rm(domain_path); Base.Filesystem.rm(problem_path)
    if verbose println(output) end
    if !occursin("Solution found", output)
        if occursin("aborting after translate", output)
            error("Could not translate domain and problem files.")
        end
        return NullSolution()
    end
    # Read plan from file
    plan = readlines("./sas_plan")[1:end-1]
    Base.Filesystem.rm("./sas_plan")
    plan = parse_pddl.(plan)
    return OrderedPlan(plan)
end


text:
# mcmc.jl

export timemod!, etr, mcmc!, adapt!, reset_adapt!;

## from Arthur Lui
function deepcopyFields(state::T, fields::Vector{Symbol}) where T
  substate = Dict{Symbol, Any}()

  for field in fields
    substate[field] = deepcopy(getfield(state, field))
  end

  return substate
end

function postSimsInit(n_keep::Int, init_state::Union{State_GPMTD},
    monitor::Vector{Symbol}=[:intercept, :lλ, :μ, :σ2, :κ,
        :κ_hypers, :corParams, :corHypers, :fx])

    L = length(init_state.mixcomps)

    monitor_outer = intersect(monitor, fieldnames(typeof(init_state)))
    monitor_mixcomps = intersect(monitor, fieldnames(typeof(init_state.mixcomps[1])))

    state = deepcopyFields(init_state, monitor_outer)
    state[:mixcomps] = [ deepcopyFields(init_state.mixcomps[j], monitor_mixcomps) for j = 1:L ]

    state[:llik] = 0.0

    sims = [ deepcopy(state) for i = 1:n_keep ]

    return sims
end

## MCMC timing for benchmarks
function timemod!(n::Int64, model::Union{Model_GPMTD}, niter::Int, outfilename::String; n_procs=1, save=false, tol_posdef=1.0e-9)
    outfile = open(outfilename, "a+")
    write(outfile, "\ntiming for $(niter) iterations each with $(n_procs) parallel processes:\n")
    for i in 1:n
        tinfo = @timed mcmc!(model, niter, n_procs=n_procs, save=save, tol_posdef=tol_posdef)
        write(outfile, "trial $(i), elapsed: $(tinfo[2]) seconds, allocation: $(tinfo[3]/1.0e6) Megabytes\n")
    end
    close(outfile)
end

## estimate time remaining
function etr(timestart::DateTime, n_iter_timed::Int, n_keep::Int, thin::Int, outfilename::String)
    timeendburn = now()
    durperiter = (timeendburn - timestart).value / float(n_iter_timed) # in milliseconds
    milsecremaining = durperiter * (n_keep * thin)
    estimatedfinish = now() + Dates.Millisecond(Int64(round(milsecremaining)))
    report_file = open(outfilename, "a+")
    write(report_file, "Completed burn-in at $(durperiter/1.0e3*1000.0) seconds per 1000 iterations \n
      $(durperiter/1.0e3/60.0*1000.0) minutes per 1000 iterations \n
      $(durperiter/1.0e3/60.0/60.0*1000.0) hours per 1000 iterations \n
      estimated completion time $(estimatedfinish)")
    close(report_file)
end

function mcmc!(model::Model_GPMTD, n_keep::Int;
    save::Bool=true,
    thin::Int=1,
    n_procs::Int=1,
    report_filename::String="out_progress.txt",
    report_freq::Int=10000,
    update::Vector{Symbol}=[:intercept, :lλ, :ζ, :μ, :σ2,
        :κ, :κ_hypers, :corParams, :corHypers, :fx],
    monitor::Vector{Symbol}=[:intercept, :lλ, :μ, :σ2,
        :κ, :κ_hypers, :corParams, :corHypers, :fx],
    tol_posdef=1.0e-9)

    ## output files
    report_file = open(report_filename, "a+")
    write(report_file, "Commencing MCMC at $(Dates.now()) for $(n_keep * thin) iterations.\n")

    ## split update parameters
    update_outer = intersect(update, fieldnames(typeof(model.state)))
    update_mixcomps = intersect(update, setdiff( fieldnames(typeof(model.state.mixcomps[1])), [:κ_hypers, :corHypers]))
    up_mixcomps = length(update_mixcomps) > 0
    update_hypers = intersect(update, [:κ_hypers, :corHypers])

    ## collect posterior samples
    if save
        sims = postSimsInit(n_keep, model.state, monitor)
        monitor_outer = intersect(monitor, fieldnames(typeof(model.state)))
        monitor_mixcomps = intersect(monitor, fieldnames(typeof(model.state.mixcomps[1])))
    end

    ## for calculating Metropolis acceptance rates
    start_accpt = deepcopy(model.state.accpt)
    prev_accpt = deepcopy(model.state.accpt)
    start_iter = deepcopy(model.state.iter)

    ## sampling
    for i in 1:n_keep
        for j in 1:thin

            if (:ζ in update_outer)
                update_ζ!(model.y, model.state, model.state.mixcomps[1])
            end

            if (:intercept in update_outer)
                update_interc_params!(model.state.intercept,
                    model.prior.intercept, model.y, model.state.ζ)
            end

            if up_mixcomps
                update_mixcomps!(model.state, model.prior, model.y, update_mixcomps, n_procs=n_procs, tol_posdef=tol_posdef)
            end

            if (:lλ in update_outer)
                update_lλ!(model.state, model.prior.λ)
            end

            if (:κ_hypers in update_hypers) || (:corHypers in update_hypers)
                update_cov_hypers!(model.state, model.prior.covhyper)
            end

            model.state.iter += 1
            if model.state.iter % report_freq == 0
                write(report_file, "Iter $(model.state.iter) at $(Dates.now())\n")
                model.state.llik = llik(model)
                write(report_file, "Log-likelihood $(model.state.llik)\n")
                write(report_file, "Current Metropolis acceptance rates: $(float((model.state.accpt - prev_accpt) / report_freq))\n\n")
                prev_accpt = deepcopy(model.state.accpt)
            end

        end

        if save
            for field in monitor_outer
                sims[i][field] = deepcopy(getfield(model.state, field))
            end
            for j = 1:model.L
                for field in monitor_mixcomps
                    sims[i][:mixcomps][j][field] = deepcopy(getfield(model.state.mixcomps[j], field))
                end
            end
            sims[i][:llik] = llik(model)
        end

    end

    model.state.llik = llik(model)

    close(report_file)

    accptr = float(model.state.accpt - start_accpt) ./ float(model.state.iter - start_iter)

    if save
        return (sims, accptr)
    else
        return (model.state.iter, accptr)
    end

end

function adjust_from_accptr(accptr::T, target::T, adjust_bnds::Array{T,1}) where T <: Real
    if accptr < target
        out = (1.0 - adjust_bnds[1]) * accptr / target + adjust_bnds[1]
    else
        out = (adjust_bnds[2] - 1.0) * (accptr - target) / (1.0 - target) + 1.0
    end

    return out
end

function adapt!(model::Model_GPMTD;
    n_iter_collectSS::Int=2000, n_iter_scale::Int=500,
    accptr_bnds::Vector{T}=[0.23, 0.44],
    adjust_bnds::Vector{T}=[0.01, 10.0],
    maxtries::Int=50,
    report_filename::String="out_progress.txt",
    update::Vector{Symbol}=[:intercept, :lλ, :ζ, :μ, :σ2,
            :κ, :κ_hypers, :corParams, :corHypers, :fx],
    tol_posdef=1.0e-9) where T <: Real

    target = (accptr_bnds[2] - accptr_bnds[1]) / 2.0
    d = Int(2)
    collect_scale = 2.38^2 / float(d)

    ## initial runs
    report_file = open(report_filename, "a+")
    write(report_file, "\n\nBeginning Adaptation Phase 1 of 4 (initial scaling) at $(Dates.now())\n\n")
    close(report_file)

    model.state.adapt = false
    reset_adapt!(model)
    tries = 0

    fails = trues(model.L)

    while any(fails)
        tries += 1
        if tries > maxtries
            println("Exceeded maximum adaptation attempts, Phase 1\n")
            report_file = open(report_filename, "a+")
            write(report_file, "\n\nExceeded maximum adaptation attempts, Phase 1\n\n")
            close(report_file)
            break
        end

        iter, accptr = mcmc!(model, n_iter_scale, save=false,
            report_filename=report_filename,
            report_freq=n_iter_scale, update=update, tol_posdef=tol_posdef)

        for j = 1:model.L
            fails[j] = (accptr[j] < accptr_bnds[1])
            if fails[j]
                model.state.mixcomps[j].cSig = PDMat_adj(adjust_from_accptr(accptr[j], target, adjust_bnds) * model.state.mixcomps[j].cSig.mat)
            end
        end

    end

    ## local scaling
    report_file = open(report_filename, "a+")
    write(report_file, "\n\nBeginning Adaptation Phase 2 of 4 (local scaling) at $(Dates.now())\n\n")
    close(report_file)

    model.state.adapt = false
    reset_adapt!(model)

    localtarget = target * 1.0

    for ii = 1:3

        for k = 1:d

            tries = 0
            fails = trues(model.L)
            while any(fails)
                tries += 1
                if tries > maxtries
                    println("Exceeded maximum adaptation attempts, Phase 2\n")
                    report_file = open(report_filename, "a+")
                    write(report_file, "\n\nExceeded maximum adaptation attempts, Phase 2\n\n")
                    close(report_file)
                    break
                end

                iter, accptr = mcmc!(model, n_iter_scale, save=false,
                    report_filename=report_filename,
                    report_freq=n_iter_scale, update=update, tol_posdef=tol_posdef)

                for j = 1:model.L
                    too_low = accptr[j] < (accptr_bnds[1] * 0.5)
                    too_high = accptr[j] > (accptr_bnds[2])

                    if too_low || too_high

                        fails[j] = true

                        tmp = Matrix(model.state.mixcomps[j].cSig)
                        σ = sqrt.(LinearAlgebra.diag(tmp))
                        ρ = StatsBase.cov2cor(tmp, σ)

                        σ[k] *= adjust_from_accptr(accptr[j], localtarget, adjust_bnds)
                        tmp = StatsBase.cor2cov(ρ, σ)
                        tmp += Diagonal(fill(0.1*minimum(σ), size(tmp,1)))
                        model.state.mixcomps[j].cSig = PDMat_adj(tmp)

                    else
                        fails[j] = false
                    end

                end

            end

        end
    end

    ## cΣ collection
    report_file = open(report_filename, "a+")
    write(report_file, "\n\nBeginning Adaptation Phase 3 of 4 (covariance collection) at $(Dates.now())\n\n")
    close(report_file)

    reset_adapt!(model)
    model.state.adapt = true

    iter, accptr = mcmc!(model, n_iter_collectSS, save=false,
        report_filename=report_filename,
        report_freq=1000, update=update, tol_posdef=tol_posdef)

    for j = 1:model.L
        Sighat = model.state.mixcomps[j].runningSS_Met / float(model.state.adapt_iter)
        # minSighat = minimum(abs.(Sighat))
        # SighatPD = Sighat + Matrix(Diagonal(fill(0.1*minSighat, d)))
        model.state.mixcomps[j].cSig = PDMat_adj(collect_scale * Sighat)
    end

    ## final scaling
    report_file = open(report_filename, "a+")
    write(report_file, "\n\nBeginning Adaptation Phase 4 of 4 (final scaling) at $(Dates.now())\n\n")
    close(report_file)

    model.state.adapt = false
    reset_adapt!(model)
    tries = 0

    fails = trues(model.L)

    while any(fails)
        tries += 1
        if tries > maxtries
            println("Exceeded maximum adaptation attempts, Phase 4\n")
            report_file = open(report_filename, "a+")
            write(report_file, "\n\nExceeded maximum adaptation attempts, Phase 4\n\n")
            close(report_file)
            break
        end

        iter, accptr = mcmc!(model, n_iter_scale, save=false,
            report_filename=report_filename,
            report_freq=n_iter_scale, update=update, tol_posdef=tol_posdef)

        for j = 1:model.L
            too_low = accptr[j] < accptr_bnds[1]
            too_high = accptr[j] > accptr_bnds[2]

            if too_low || too_high
                fails[j] = true
                model.state.mixcomps[j].cSig = PDMat_adj(adjust_from_accptr(accptr[j], target, adjust_bnds) * model.state.mixcomps[j].cSig.mat)
            else
                fails[j] = false
            end

        end

    end

    report_file = open(report_filename, "a+")
    write(report_file, "\n\nAdaptation concluded at $(Dates.now())\n\n")
    close(report_file)

    ## note that mcmc! also closes the report file
    reset_adapt!(model)
    model.state.adapt = false

    return nothing
end

function reset_adapt!(model::Model_GPMTD, nparams::Int=2)
    model.state.adapt_iter = 0
    L = length(model.state.mixcomps)
    for j = 1:L
        model.state.mixcomps[j].adapt_iter = 0
        model.state.mixcomps[j].runningsum_Met = zeros( Float64, nparams )
        model.state.mixcomps[j].runningSS_Met = zeros( Float64, nparams, nparams )
    end

    return nothing
end


text:
@testset "mapmany" begin
    @test mapmany(i->1:i, 1:3) == [1, 1,2, 1,2,3]
end

@testset "flatten" begin
    @test flatten([1:1, 1:2, 1:3]) == [1, 1,2, 1,2,3]
end


text:
function normalize(X)
    # Normalizes each column by the standard normal
    n,p = size(X)
    X_n = zeros(n,p)
    means = [sum(X[:,i])/n for i=1:p];
    vars = [sum((X[:,i] - means[i]*ones(n,1)).^2)/n for i=1:p];
    for i=1:p
        X_n[:,i] = (X[:,i] - means[i]*ones(n,1))./vars[i]^(0.5);
    end
    return X_n, means, vars
end

function threshold(Y, val)
    # Thresholds Y and turns into binary categories
    return float(Y .>= val)
end

function augment(X, fns, fn_names=nothing)
    # Augmenting the data with a set of fn's, and returns associated functions
    # This method prefers that fns represent operations on columns, not matrices!
    n,p = size(X);
    augX = deepcopy(X);
    fn_array = [x -> x[:,i] for i=1:p];
    for i=1:size(fns,1)
        next = fns[i](X);
        if typeof(X) == DataFrame && !isnothing(fn_names)
            names!(next, fn_names[i])
        end
        augX = hcat(augX, next);
        fn_array= vcat(fn_array, [fns[i]]);
    end
    return augX
end

function quad_fns(x)
    n,p = size(x);
    fns = [x->x[:,i].^2 for i = 1:p];
    fn_names = [];
    if typeof(x) == DataFrame
        fn_names = [string("x_", i, "^2") for i =1:p];
    end
    return fns, fn_names
end

function bilinear_fns(x)
    n,p = size(x);
    fn_arr = [x->x[:,i].^2 for i = 1:p];
    fn_names = [];
    if typeof(x) == DataFrame
        fn_names = [string("x_", i, "^2") for i =1:p];
    end
    for i = 1:p
        fn_arr = vcat(fn_arr, [x->x[:,i].*x[:,j] for j=i+1:p])
        if typeof(x) == DataFrame
            fn_names = vcat(fn_names, [string("x_", i,"*", "x_", j ) for j=i+1:p]);
        end
    end
    return fn_arr, fn_names
end

function sixth_order_even_fns(x)
    n,p = size(x);
    fn_arr = [x->x[:,i].^2 for i = 1:p];
    fn_arr = vcat(fn_arr,[x->x[:,i].^4 for i = 1:p]);
    fn_arr = vcat(fn_arr, [x->x[:,i].^6 for i = 1:p];)
    fn_names = [];
    if typeof(x) == DataFrame
        fn_names = [string("x_", i, "^2") for i =1:p];
        fn_names = vcat(fn_names,[string("x_", i, "^4") for i =1:p]);
        fn_names = vcat(fn_names, [string("x_", i, "^6") for i =1:p]);
    end
    for i = 1:p
        fn_arr = vcat(fn_arr, [x->x[:,i].^2 .* x[:,j].^4 for j=1:p])
        if typeof(x) == DataFrame
            fn_names = vcat(fn_names, [string("x_", i,"^2*", "x_", j, "^4*") for j=1:p]);
        end
    end
    return fn_arr, fn_names
end

function exps(x)
    n,p = size(x)
    return [x->exp.(x[:,i]) for i=1:p]
end


text:
using ZipFile

v = "0.4.2"

url = "https://github.com/tailwindlabs/heroicons/archive/v$v.zip"

rm(joinpath(@__DIR__, "heroicons-$v/"), force=true, recursive=true)

download(url, joinpath(@__DIR__, "heroicons-$v.zip"))

zipfile = ZipFile.Reader(joinpath(@__DIR__, "heroicons-$v.zip"))

for file in zipfile.files
    if endswith(file.name, '/')
        mkdir(joinpath(@__DIR__, file.name))
    elseif endswith(file.name, ".svg")
        open(joinpath(@__DIR__, file.name), "w") do io
            print(io, read(file, String))
        end
    end
end


text:
module OTSIN

using LinearAlgebra

include("optimaltransport.jl")
include("metrics.jl")
include("fitting.jl")
include("utils.jl")
include("perturbation.jl")
include("posterior.jl")
include("plotting.jl")

export isprobability, marginals
export optimaltransport
export entropy, relative_entropy, KL, utility
export uniform, report_utility, normalize
export fitM, fitMsgd, r, global_loss
export Regularization, L1, L2, ElasticNet, Exp, Entropic
export perturbation, pert_utility, pert_entropy
export parallel_cross_entropy, parallel_cross_entropy_A, parallel_cross_entropy_B,
        parallel_cross_entropy_free
export laplace_posterior
export marginalsheatmap

end # module

max_stars_repo_name: PoisotLab/OTSIN

max_issues_repo_name: PoisotLab/OTSIN

max_forks_repo_name: PoisotLab/OTSIN


text:
using PQEcolaPoint
using Dates
using DataFrames

# simulation dictionary -
simulation_archive = Dict{NamedTuple,DataFrame}()

# build array of models -
models = Array{CRRLatticeModel,1}()
Sₒ = 81.68
σ = Sₒ*(0.50)*sqrt(37/365)
L = Sₒ - 2.96*σ
U = Sₒ + 2.96*σ
DTEₒ = 37
underlying_price_range = range(L,stop=U,length = 10) |> collect
IV_array = range(0.1, stop = 1.0, step = 0.1) |> collect
strike_array = range(50, stop = 125, step = 5) |> collect
dte_array = range(DTEₒ, stop = 1.0, step=-1) |> collect

# build contract array -
contracts = Array{PutContractModel,1}()
for K ∈ strike_array

    # what is the contract?
    short_put_contract = PutContractModel()
    short_put_contract.ticker = "XYZ"
    short_put_contract.expiration_date = Date(2022, 05, 06)
    short_put_contract.strike_price = K
    short_put_contract.premium = 0.0
    short_put_contract.number_of_contracts = 1
    short_put_contract.direction = 1
    push!(contracts, short_put_contract)
end

for (i,underlying_price) ∈ enumerate(underlying_price_range)
    for (j,dte_value) ∈ enumerate(dte_array)
        for IV ∈ IV_array

            # build a binary tree with these market conditions -
            model = build(CRRLatticeModel; Sₒ = underlying_price, number_of_levels = 80, σ = IV, T = (dte_value / 365), μ = 0.0045)
            push!(models, model)
        end

        # compute grid -
        df = premium(contracts, models)

        # build the key -
        key_tuple = (
            i = i,
            j = j
        )

        # cache -
        simulation_archive[key_tuple] = df

        # clean out the model array -
        empty!(models)

        @info (underlying_price, dte_value)
    end
end


text:
using DataSets.DataREPL: complete, parse_data_repl_cmd

@testset "repl completions" begin
    @test complete("") == (["list ", "show ", "stack ", "help "], "", true)
    @test complete("s") == (["show ", "stack "], "s", true)
    @test complete("stack ") == (["push ", "pop ", "list "], "", true)

    cd(@__DIR__) do
        if Sys.iswindows()
            @test complete("stack push da") == (["data\\"], "da", true)
        else
            @test complete("stack push da") == (["data/"], "da", true)
        end
    end
end

@testset "repl commands" begin
    @test eval(parse_data_repl_cmd("help")) === DataSets.DataREPL._data_repl_help
    @test eval(parse_data_repl_cmd("?")) === DataSets.DataREPL._data_repl_help
    empty!(DataSets.PROJECT)
    @test eval(parse_data_repl_cmd("stack push $(@__DIR__)")) === DataSets.PROJECT
    @test length(DataSets.PROJECT.projects) == 1
    @test eval(parse_data_repl_cmd("stack pop")) === DataSets.PROJECT
    @test isempty(DataSets.PROJECT.projects)
end

@testset "data show utils" begin
    @test sprint(DataSets.DataREPL.hexdump, UInt8.(0:70)) == raw"""
    0000: 0001 0203 0405 0607 0809 0a0b 0c0d 0e0f  ................
    0010: 1011 1213 1415 1617 1819 1a1b 1c1d 1e1f  ................
    0020: 2021 2223 2425 2627 2829 2a2b 2c2d 2e2f   !"#$%&'()*+,-./
    0030: 3031 3233 3435 3637 3839 3a3b 3c3d 3e3f  0123456789:;<=>?
    0040: 4041 4243 4445 46                        @ABCDEF
    """
end


text:
using Documenter, Plots

makedocs(modules=[Plots], doctest=false)

deploydocs(
    deps   = Deps.pip("mkdocs",
    "mkdocs-material" ,"python-markdown-math", "pygments", "pymdown-extensions"),
    repo   = "github.com/JuliaPlots/PlotDocs.jl.git",
    julia  = "0.6",
    osname = "linux"
)


text:
X = collect(Float64, reshape(1:40, 10, 4))
x = rand(10) * 10
w = collect(Float64, 1:10)

D = DataFrame(A=rand(10), B=1:10, C=string.(11:20; base=16))
D_NA = allowmissing!(copy(D))
D_NA[1, :A] = missing
@testset "Array" begin
    scaler = FixedRangeScaler(X)
    XX = transform(X, scaler)
    @test mean(XX[:,end]) ≈ 1
    @test mean(XX[:,1]) ≈ 0
    @test maximum(XX) == 1
    @test minimum(XX) == 0

    XX, scaler = fit_transform(FixedRangeScaler, X)
    @test mean(XX[:,end]) ≈ 1
    @test mean(XX[:,1]) ≈ 0
    @test maximum(XX) == 1
    @test minimum(XX) == 0

    XX = deepcopy(X)
    scaler = fit_transform!(FixedRangeScaler, XX)
    @test mean(XX[:,end]) ≈ 1
    @test mean(XX[:,1]) ≈ 0
    @test maximum(XX) == 1
    @test minimum(XX) == 0

    scaler = FixedRangeScaler(X, -2, 2)
    XX = transform(X, scaler)
    @test mean(XX[:,end]) ≈ 2
    @test mean(XX[:,1]) ≈ -2
    @test maximum(XX) == 2
    @test minimum(XX) == -2

    scaler = fit(FixedRangeScaler, X)
    XX = transform(X, scaler)
    @test mean(XX[:,end]) ≈ 1
    @test mean(XX[:,1]) ≈ 0
    @test maximum(XX) == 1
    @test minimum(XX) == 0

    scaler = fit(FixedRangeScaler, X, obsdim=1)
    XX = transform(X, scaler)
    @test mean(XX[1,:]) ≈ 0
    @test mean(XX[end,:]) ≈ 1
    @test maximum(XX) == 1
    @test minimum(XX) == 0

    scaler = fit(FixedRangeScaler, X, -2, 2)
    XX = transform(X, scaler)
    @test mean(XX[:,end]) ≈ 2
    @test mean(XX[:,1]) ≈ -2
    @test maximum(XX) == 2
    @test minimum(XX) == -2

    XX, scaler = fit_transform(FixedRangeScaler, X, -2, 2)
    @test mean(XX[:,end]) ≈ 2
    @test mean(XX[:,1]) ≈ -2
    @test maximum(XX) == 2
    @test minimum(XX) == -2

    XX = deepcopy(X)
    scaler = fit_transform!(FixedRangeScaler, XX, -2, 2)
    @test mean(XX[:,end]) ≈ 2
    @test mean(XX[:,1]) ≈ -2
    @test maximum(XX) == 2
    @test minimum(XX) == -2

    scaler = fit(FixedRangeScaler, X, -2, 2, obsdim=1)
    XX = transform(X, scaler)
    @test mean(XX[1,:]) ≈ -2
    @test mean(XX[end,:]) ≈ 2
    @test maximum(XX) == 2
    @test minimum(XX) == -2

    scaler = fit(FixedRangeScaler, X, -2, 2, obsdim=2)
    XX = transform(X, scaler)
    @test mean(minimum(XX; dims=2)) ≈ -2
    @test mean(maximum(XX; dims=2)) ≈ 2
    @test maximum(XX) == 2
    @test minimum(XX) == -2

    scaler = fit(FixedRangeScaler, X, -2, 2, obsdim=1, operate_on=[1,2])
    XX = transform(X, scaler)
    @test mean(minimum(XX[:,[1,2]], dims=1)) ≈ -2
    @test mean(maximum(XX[:,[1,2]], dims=1)) ≈ 2
    ww = deepcopy(w)
    transform!(ww, scaler)
    @test ww[1] == -2
    @test all(w[3:end] .== ww[3:end])

    scaler = fit(FixedRangeScaler, X, -2, 2, obsdim=2, operate_on=[1,2])
    XX = transform(X, scaler)
    @test mean(minimum(XX[[1,2],:], dims=2)) ≈ -2
    @test mean(maximum(XX[[1,2],:], dims=2)) ≈ 2

    Xi = round.(Int, X)
    transform(Xi, scaler)
    @test mean(minimum(XX[[1,2],:], dims=2)) ≈ -2
    @test mean(maximum(XX[[1,2],:], dims=2)) ≈ 2

    scaler = fit(FixedRangeScaler, X, -2, 2, obsdim=2, operate_on=[1,2])
    XX = deepcopy(X)
    transform!(XX, scaler)
    @test mean(minimum(XX[[1,2],:], dims=2)) ≈ -2
    @test mean(maximum(XX[[1,2],:], dims=2)) ≈ 2

    XX = deepcopy(X)
    fixedrange!(XX)
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)))
    @test all(minimum(XX; dims=2) .== zeros(size(XX, 1)))

    XX = deepcopy(X)
    fixedrange!(XX, ObsDim.Last(), collect(1:size(X,1)))
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)))
    @test all(minimum(XX; dims=2) .== zeros(size(XX, 1)))

    XX = deepcopy(X)
    fixedrange!(XX, -2, 2)
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)) .* 2)
    @test all(minimum(XX; dims=2) .== ones(size(XX, 1)) .* -2)

    XX = deepcopy(X)
    fixedrange!(XX, 0, 1, ObsDim.First(), collect(1:size(X, 2)))
    @test all(maximum(XX; dims=1) .== ones(size(XX, 2)))
    @test all(minimum(XX; dims=1) .== zeros(size(XX, 2)))

    XX = deepcopy(X)
    fixedrange!(XX, 0, 1, ObsDim.Last(), collect(1:size(X,1)))
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)))
    @test all(minimum(XX; dims=2) .== zeros(size(XX, 1)))

    XX = deepcopy(X)
    fixedrange!(XX, 0, 1, vec(minimum(XX; dims=2)), vec(maximum(XX; dims=2)))
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)))
    @test all(minimum(XX; dims=2) .== zeros(size(XX, 1)))

    XX = deepcopy(X)
    fixedrange!(XX, 0, 1, vec(minimum(XX; dims=2)), vec(maximum(XX; dims=2)), ObsDim.Last(), collect(1:size(XX, 1)))
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)))
    @test all(minimum(XX; dims=2) .== zeros(size(XX, 1)))

    XX = deepcopy(X)
    fixedrange!(XX, 0, 1, vec(minimum(XX; dims=2)), vec(maximum(XX; dims=2)), ObsDim.Constant{2}(), collect(1:size(XX, 1)))
    @test all(maximum(XX; dims=2) .== ones(size(XX, 1)))
    @test all(minimum(XX; dims=2) .== zeros(size(XX, 1)))

    xx = deepcopy(x)
    fixedrange!(xx)
    @test minimum(xx) == 0
    @test maximum(xx) == 1

    xx = deepcopy(x)
    fixedrange!(xx, -1, 1)
    @test minimum(xx) == -1
    @test maximum(xx) == 1

    xx = deepcopy(x)
    xmin = minimum(x; dims=2) .- 1
    xmax = maximum(x; dims=2)

    xx = deepcopy(x)
    fixedrange!(xx, -1, 1, xmin, xmax, ObsDim.First(), collect(1:length(x)))
    @test all(minimum(xx; dims=2) .== ones(length(x)))
    @test all(maximum(xx; dims=2) .== ones(length(x)))

    xx = deepcopy(x)
    fixedrange!(xx, -1, 1, xmin, xmax, ObsDim.Last(), collect(1:length(x)))
    @test all(minimum(xx; dims=2) .== ones(length(x)))
    @test all(maximum(xx; dims=2) .== ones(length(x)))
end

@testset "DataFrame" begin
    scaler = FixedRangeScaler(D)
    DD = transform(D, scaler)
    @test minimum(DD[:A]) == 0
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == 0
    @test maximum(DD[:B]) == 1

    scaler = fit(FixedRangeScaler, D)
    DD = transform(D, scaler)
    @test minimum(DD[:A]) == 0
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == 0
    @test maximum(DD[:B]) == 1

    DD, scaler = fit_transform(FixedRangeScaler, D)
    @test minimum(DD[:A]) == 0
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == 0
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    scaler = fit_transform!(FixedRangeScaler, DD)
    @test minimum(DD[:A]) == 0
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == 0
    @test maximum(DD[:B]) == 1

    scaler = fit(FixedRangeScaler, D, -1, 1)
    DD = transform(D, scaler)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD, scaler = fit_transform(FixedRangeScaler, D, -1, 1)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    scaler = fit_transform!(FixedRangeScaler, DD, -1, 1)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    scaler = FixedRangeScaler(D, -1, 1)
    DD = transform(D, scaler)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    scaler = fit(FixedRangeScaler, D, -1, 1, operate_on=[:A])
    DD = transform(D, scaler)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == minimum(D[:B])
    @test maximum(DD[:B]) == maximum(D[:B])

    scaler = fit(FixedRangeScaler, D, -1, 1)
    DD = transform(D_NA, scaler)
    @test ismissing(DD[1,:A])
    @test DD[end,:A] == D_NA[end,:A]
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    scaler = fit(FixedRangeScaler, D, -1, 1, operate_on=[:A, :B])
    DD = transform(D_NA, scaler)
    @test ismissing(DD[1,:A])
    @test DD[end,:A] == D_NA[end,:A]
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    scaler = fit(FixedRangeScaler, D, -1, 1, operate_on=[:A, :B, :C])
    DD = transform(D_NA, scaler)
    @test ismissing(DD[1,:A])
    @test DD[end,:A] == D_NA[end,:A]
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    scaler = fit(FixedRangeScaler, D, -1, 1, operate_on=[:A, :B])
    transform!(DD, scaler)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    fixedrange!(DD)
    @test minimum(DD[:A]) == 0
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == 0
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    fixedrange!(DD, -1, 1)
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    fixedrange!(DD, -1, 1, [:A, :B, :C])
    @test minimum(DD[:A]) == -1
    @test maximum(DD[:A]) == 1
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D_NA)
    fixedrange!(DD, -1, 1, [:A, :B, :C])
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1

    DD = deepcopy(D)
    fixedrange!(DD, -1, 1, [0,1], [1,10])
    @test minimum(DD[:B]) == -1
    @test maximum(DD[:B]) == 1
end


text:
"""
    retrieving!(actr, chunk, args...; kwargs...)

Sets the declarative memory module as busy and Submits a request for a chunk and registers
a new event for the retrieval

# Arguments

- `actr`: an ACT-R model object
- `request...`: a variable list of slot-value pairs
"""
function retrieving!(actr; request...)
    actr.declarative.state.busy = true
    description = "Retrieve"
    type = "model"
    id = get_name(actr)
    cur_time = get_time(actr)
    chunk = retrieve(actr, cur_time; request...)
    tΔ = compute_RT(actr, chunk)
    register!(actr, retrieve!, after, tΔ , actr, chunk; id, description, type)
    return tΔ
end

"""
    retrieve!(actr, chunk, args...; kwargs...)

Completes a memory retrieval by adding chunk to declarative memory buffer and setting
busy = false and empty = false. Error is set to true if retrieval failure occurs.

# Arguments

- `actr`: an ACT-R model object
- `request...`: a variable list of slot-value pairs
"""
function retrieve!(actr, chunk)
    actr.declarative.state.busy = false
    if isempty(chunk)
        actr.declarative.state.error = true
    else
        actr.declarative.state.empty = false
        add_to_buffer!(actr.declarative, chunk[1])
    end
    return nothing
end


text:
const _cl_error_codes = Dict{Int, Symbol}(
     +0 => :CL_SUCCESS,
     -1 => :CL_DEVICE_NOT_FOUND,
     -2 => :CL_DEVICE_NOT_AVAILABLE,
     -3 => :CL_COMPILER_NOT_AVAILABLE,
     -4 => :CL_MEM_OBJECT_ALLOCATION_FAILURE,
     -5 => :CL_OUT_OF_RESOURCES,
     -6 => :CL_OUT_OF_HOST_MEMORY,
     -7 => :CL_PROFILING_INFO_NOT_AVAILABLE,
     -8 => :CL_MEM_COPY_OVERLAP,
     -9 => :CL_IMAGE_FORMAT_MISMATCH,
    -10 => :CL_IMAGE_FORMAT_NOT_SUPPORTED,
    -11 => :CL_BUILD_PROGRAM_FAILURE,
    -12 => :CL_MAP_FAILURE,
    -13 => :CL_MISALIGNED_SUB_BUFFER_OFFSET,
    -14 => :CL_EXEC_STATUS_ERROR_FOR_EVENTS_IN_WAIT_LIST,
    -15 => :CL_COMPILE_PROGRAM_FAILURE,
    -15 => :CL_LINKER_NOT_AVAILABLE,
    -17 => :CL_LINK_PROGRAM_FAILURE,
    -18 => :CL_DEVICE_PARTITION_FAILED,
    -19 => :CL_KERNEL_ARG_INFO_NOT_AVAILABLE,

    -30 => :CL_INVALID_VALUE,
    -31 => :CL_INVALID_DEVICE_TYPE,
    -32 => :CL_INVALID_PLATFORM,
    -33 => :CL_INVALID_DEVICE,
    -34 => :CL_INVALID_CONTEXT,
    -35 => :CL_INVALID_QUEUE_PROPERTIES,
    -36 => :CL_INVALID_COMMAND_QUEUE,
    -37 => :CL_INVALID_HOST_PTR,
    -38 => :CL_INVALID_MEM_OBJECT,
    -39 => :CL_INVALID_IMAGE_FORMAT_DESCRIPTOR,
    -40 => :CL_INVALID_IMAGE_SIZE,
    -41 => :CL_INVALID_SAMPLER,
    -42 => :CL_INVALID_BINARY,
    -43 => :CL_INVALID_BUILD_OPTIONS,
    -44 => :CL_INVALID_PROGRAM,
    -45 => :CL_INVALID_PROGRAM_EXECUTABLE,
    -46 => :CL_INVALID_KERNEL_NAME,
    -47 => :CL_INVALID_KERNEL_DEFINITION,
    -48 => :CL_INVALID_KERNEL,
    -49 => :CL_INVALID_ARG_INDEX,
    -50 => :CL_INVALID_ARG_VALUE,
    -51 => :CL_INVALID_ARG_SIZE,
    -52 => :CL_INVALID_KERNEL_ARGS,
    -53 => :CL_INVALID_WORK_DIMENSION,
    -54 => :CL_INVALID_WORK_GROUP_SIZE,
    -55 => :CL_INVALID_WORK_ITEM_SIZE,
    -56 => :CL_INVALID_GLOBAL_OFFSET,
    -57 => :CL_INVALID_EVENT_WAIT_LIST,
    -58 => :CL_INVALID_EVENT,
    -59 => :CL_INVALID_OPERATION,
    -60 => :CL_INVALID_GL_OBJECT,
    -61 => :CL_INVALID_BUFFER_SIZE,
    -62 => :CL_INVALID_MIP_LEVEL,
    -63 => :CL_INVALID_GLOBAL_WORK_SIZE,
    -64 => :CL_INVALID_PROPERTY,
    -65 => :CL_INVALID_IMAGE_DESCRIPTOR,
    -66 => :CL_INVALID_COMPILER_OPTIONS,
    -67 => :CL_INVALID_LINKER_OPTIONS,
    -68 => :CL_INVALID_DEVICE_PARTITION_COUNT,
    -69 => :CL_INVALID_PIPE_SIZE,
    -70 => :CL_INVALID_DEVICE_QUEUE,

    -1000 => :CL_INVALID_GL_SHAREGROUP_REFERENCE_KHR,
    -1001 => :CL_PLATFORM_NOT_FOUND_KHR,
    -1002 => :CL_INVALID_D3D10_DEVICE_KHR,
    -1003 => :CL_INVALID_D3D10_RESOURCE_KHR,
    -1004 => :CL_D3D10_RESOURCE_ALREADY_ACQUIRED_KHR,
    -1005 => :CL_D3D10_RESOURCE_NOT_ACQUIRED_KHR,
    -1006 => :CL_INVALID_D3D11_DEVICE_KHR,
    -1007 => :CL_INVALID_D3D11_RESOURCE_KHR,
    -1008 => :CL_D3D11_RESOURCE_ALREADY_ACQUIRED_KHR,
    -1009 => :CL_D3D11_RESOURCE_NOT_ACQUIRED_KHR,
    -1010 => :CL_INVALID_DX9_MEDIA_ADAPTER_KHR,
    -1011 => :CL_INVALID_DX9_MEDIA_SURFACE_KHR,
    -1012 => :CL_DX9_MEDIA_SURFACE_ALREADY_ACQUIRED_KHR,
    -1013 => :CL_DX9_MEDIA_SURFACE_NOT_ACQUIRED_KHR,

    -1057 => :CL_DEVICE_PARTITION_FAILED_EXT,
    -1058 => :CL_INVALID_PARTITION_COUNT_EXT,
    -1059 => :CL_INVALID_PARTITION_NAME_EXT,

    -1092 => :CL_EGL_RESOURCE_NOT_ACQUIRED_KHR,
    -1093 => :CL_INVALID_EGL_OBJECT_KHR,
)

const _cl_err_desc = Dict{Integer, String}(
    CL_INVALID_CONTEXT =>
    "Context is not a valid context.",

    CL_INVALID_BUFFER_SIZE =>
    "Buffer size is 0",

    CL_INVALID_EVENT =>
    "Event objects specified in event_list are not valid event objects",

    CL_INVALID_HOST_PTR =>
    string("If host_ptr is NULL and CL_MEM_USE_HOST_PTR or ",
           "CL_MEM_COPY_HOST_PTR are set in flags or if host_ptr is not NULL but ",
           "CL_MEM_COPY_HOST_PTR or CL_MEM_USE_HOST_PTR are not set in flags."),

    CL_MEM_OBJECT_ALLOCATION_FAILURE =>
    "Failure to allocate memory for buffer object.",

    CL_OUT_OF_RESOURCES =>
    "Failure to allocate resources required by the OpenCL implementation on the device.",

    CL_OUT_OF_HOST_MEMORY =>
    "Failure to allocate resources required by the OpenCL implementation on the host",

    CL_INVALID_PROGRAM =>
    "Program is not a valid program object.",

    CL_INVALID_VALUE =>
    "CL_INVALID_VALUE: this one should have been caught by julia!",

    CL_INVALID_DEVICE =>
    "OpenCL devices listed in device_list are not in the list of devices associated with program.",

    CL_INVALID_BINARY =>
    string("program is created with clCreateWithProgramBinary and devices listed in ",
           "device_list do not have a valid program binary loaded."),

    CL_INVALID_BUILD_OPTIONS =>
    "The build options specified by options are invalid.",

    CL_INVALID_OPERATION =>
    string("The build of a program executable for any of the devices listed in device_list by a ",
           "previous call to clBuildProgram for program has not  completed."),

    CL_COMPILER_NOT_AVAILABLE =>
    "Program is created with clCreateProgramWithSource and a compiler is not available",

    CL_BUILD_PROGRAM_FAILURE =>
    string("Failure to build the program executable. ",
           "This error will be returned if clBuildProgram ",
           "does not return until the build has completed"),

    CL_INVALID_OPERATION =>
    "There are kernel objects attached to program.",

    CL_OUT_OF_HOST_MEMORY =>
    "if there is a failure to allocate resources required by the OpenCL implementation on the host.",

    CL_EXEC_STATUS_ERROR_FOR_EVENTS_IN_WAIT_LIST =>
    "The execution status of any of the events in event_list is a negative integer value",

    CL_INVALID_PROGRAM_EXECUTABLE =>
    "there is no successfully built executable for program",

    CL_INVALID_KERNEL_NAME =>
    "kernel_name is not found in program.",

    CL_INVALID_KERNEL_DEFINITION =>
    string("The function definition for __kernel  function ",
           "given by kernel_name such as the number of arguments, the argument types are not the ",
           "same for all devices for which the program executable has been built"),

    CL_PROFILING_INFO_NOT_AVAILABLE =>
    string("The CL_QUEUE_PROFILING_ENABLE flag ",
           "is not set for the command-queue, if the execution status of the command identified by ",
           "event is not CL_COMPLETE or if event is a user event objec"),
)

struct CLMemoryError <: Exception
    msg::String
end

Base.show(io::IO, err::CLMemoryError) = Base.print(io, "OpenCL.MemObject: $(err.msg)")

struct OpenCLException <: Exception
    msg::String
end

Base.show(io::IO, err::OpenCLException) = Base.print(io, "OpenCL Error: $(err.msg)")

struct CLError <: Exception
    code::CL_int
    desc::Symbol

    function CLError(c::Integer)
        new(c, get(_cl_error_codes, Int(c), :CL_UNKNOWN_ERROR_CODE))
    end
end

Base.show(io::IO, err::CLError) =
        Base.print(io, "CLError(code=$(err.code), $(err.desc))")

function error_description(err::CLError)
    get(_cl_err_desc, err.code, "no description for error $(err.code)")
end


text:
# # Filtered Oscillator

# ## Model

module filtered_oscillator      #jl

using ReachabilityAnalysis

one_loop_iteration = false
n0 = 4
n1 = (one_loop_iteration ? n0 + 1 : n0)
n = n1 + 2
z = zeros(n1)

## common flow
A = zeros(n, n)
A[1,1], A[2,2] = -2., -1.
A[3,1], A[3,3] = 5., -5.
for i = 4 : n-1
    A[i,i-1], A[i,i] = 5., -5.
end

function mode1(z)
    b = [1.4; -0.7; z]
    X = HPolyhedron([HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                     HalfSpace([1.0; 0.0; z], 0.0)])  # x <= 0
    @system(x' = Ax + b, x ∈ X)
end

function mode2(z)
    b = [-1.4; 0.7; z]
    X = HPolyhedron([HalfSpace([1.0; 0.0; z], 0.0),  # x <= 0
                     HalfSpace([0.714286; 1.0; z], 0.0)])  # 0.714286*x + y <= 0
    @system(x' = Ax + b, x ∈ X)
end

function mode3(z)
    b = [1.4; -0.7; z]
    X = HPolyhedron([HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                     HalfSpace([-0.714286; -1.0; z], 0.0)])  # 0.714286*x + y >= 0
    @system(x' = Ax + b, x ∈ X)
end

function mode4(z, one_loop_iteration)
    b = [-1.4; 0.7; z]
    X = HPolyhedron([HalfSpace([0.714286; 1.0; z], 0.0),  # 0.714286*x + y <= 0
                     HalfSpace([-1.0; 0.0; z], 0.0)])  # x >= 0

    if one_loop_iteration
        ## k <= 2 (2.1 to avoid numerical issues)
        addconstraint!(X, HalfSpace([zeros(n-1); 1.], 2.1))
    end
    @system(x' = Ax + b, x ∈ X)
end

function filtered_oscillator_hybrid(n0, one_loop_iteration)

    n1 = (one_loop_iteration ? n0 + 1 : n0)
    n = n1 + 2
    z = zeros(n1)

    ## transition graph (automaton)
    a = LightAutomaton(4)
    add_transition!(a, 3, 4, 1)
    add_transition!(a, 4, 2, 2)
    add_transition!(a, 2, 1, 3)
    add_transition!(a, 1, 3, 4)

    mode1 = mode1(z)
    mode2 = mode2(z)
    mode3 = mode3(z)
    mode4 = mode4(z, one_loop_iteration)
    m = [mode1, mode2, mode3, mode4]

    ## transitions

    ## transition l3 -> l4
    X_l3l4 = HPolyhedron([HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                          HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                          HalfSpace([0.714286; 1.0; z], 0.0)])  # 0.714286*x + y <= 0
    if one_loop_iteration
        A_trans_34 = Matrix(1.0I, n, n)
        A_trans_34[n, n] = 2.  # k' = k * 2
        r1 = ConstrainedLinearMap(A_trans_34, X_l3l4)
    else
        r1 = ConstrainedIdentityMap(n, X_l3l4)
    end

    ## transition l4 -> l2
    X_l4l2 = HPolyhedron([HalfSpace([0.714286; 1.0; z], 0.0),  # 0.714286*x + y <= 0
                          HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                          HalfSpace([1.0; 0.0; z], 0.0)])  # x <= 0
    r2 = ConstrainedIdentityMap(n, X_l4l2)

    ## transition l2 -> l1
    X_l2l1 = HPolyhedron([HalfSpace([1.0; 0.0; z], 0.0),  # x <= 0
                          HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                          HalfSpace([0.714286; 1.0; z], 0.0)])  # 0.714286*x + y <= 0
    r3 = ConstrainedIdentityMap(n, X_l2l1)

    ## transition l1 -> l3
    X_l1l3 = HPolyhedron([HalfSpace([-0.714286; -1.0; z], 0.0),  # 0.714286*x + y >= 0
                          HalfSpace([-1.0; 0.0; z], 0.0),  # x >= 0
                          HalfSpace([1.0; 0.0; z], 0.0)])  # x <= 0
    r4 = ConstrainedIdentityMap(n, X_l1l3)

    r = [r1, r2, r3, r4]

    ## switchings
    s = [HybridSystems.AutonomousSwitching()]

    return HybridSystem(a, m, r, s)
end

function model(X0; n0::Int=4, one_loop_iteration::Bool=false)
    H = filtered_oscillator_hybrid(n0, one_loop_iteration)
    return IVP(H, [(1, X0)])
end

end #jl

# ## References

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
# TODO: generalize to K-integrator setting
struct NPlayer2DDoubleIntegratorCost{TIC,TICR,TSC, TSCR,TPC,TGC} <: NPlayerNavigationCost
    "the index of the player this cost applies to"
    player_id::Int
    inputcost::TIC
    inputconstr::TICR
    statecost::TSC
    stateconstr::TSCR
    proximitycost::TPC
    goalcost::TGC
end

function NPlayer2DDoubleIntegratorCost(player_id, xg, t_final, np;
         inputcost::TIC=QuadCost(SMatrix{2,2}([1. 0.; 0. 1.]) * 10),
         inputconstr::TICR=tuple(),
         statecost::TSC=QuadCost(SMatrix{4,4}(diagm([0., 1., 0., 1.])) * 40),
         # TODO: think of a way to enforce v > 0
         stateconstr::TSCR=(SoftConstr(2, -2., 2., 50),
                            SoftConstr(4, -2., 2., 50)),
         proximitycost::TPC=ProximityCost(2.0, 50.0, np),
         goalcost::TGC=GoalCost(t_final, xg, SMatrix{4,4}(diagm([1.,0.,1.,0.])) *
                                500)) where {TXI,TUI,TIC,TICR,TSC,TSCR,TPC,TGC}

    return NPlayer2DDoubleIntegratorCost(player_id, inputcost, inputconstr,
                                         statecost, stateconstr, proximitycost,
    goalcost)
end

"-----------------Implementing the NPlayerNavigationCost interface-----------------"

player_id(c::NPlayer2DDoubleIntegratorCost) = c.player_id
inputcost(c::NPlayer2DDoubleIntegratorCost) = c.inputcost
inputconstr(c::NPlayer2DDoubleIntegratorCost) = c.inputconstr
statecost(c::NPlayer2DDoubleIntegratorCost) = c.statecost
stateconstr(c::NPlayer2DDoubleIntegratorCost) = c.stateconstr
proximitycost(c::NPlayer2DDoubleIntegratorCost) = c.proximitycost
goalcost(c::NPlayer2DDoubleIntegratorCost) = c.goalcost


text:
"""
    randBall(nDims, nPts)

Draw nPts random points uniformely distributed in an nDims-dimensional ball.

# Implementation
Samples uniformely in a unit cube, then throws away unnecessary points.
Efficient only for low number of dimensions (nDims=1,2,3).

# Examples
```julia-repl
julia> using FPRS; using Plots;
julia> x = randBall(2,50);
julia> figure(); scatter(x[1,:], x[2,:]);
```
"""
function randBall(nDims::Integer, nPts::Integer)

    res = zeros(nDims, nPts)
    numInBall = 0

    while numInBall < nPts
        newPoint = 2.0 * rand(nDims) .- 1.0
        if norm(newPoint) <= 1
            numInBall += 1
            res[1:end, numInBall] = newPoint
        end
    end

    return res
end

"""
    randRotationMatrix(deflection=1.0; randSeed)

Return a random rotation matrix in 3d.

Sources:
-http://blog.lostinmyterminal.com/python/2015/05/12/random-rotation-matrix.html
-http://www.realtimerendering.com/resources/GraphicsGems/gemsiii/rand_rotation.c

# Arguments
- `deflection`: number between 0 and 1 controls the size of the perturbation.
                For 0, no rotation; for 1, competely random rotation.
                Small deflection => small perturbation.

# Examples
```julia-repl
julia> mx = randRotationMatrix();
```
"""
function randRotationMatrix(deflection=1.0)
    # Random.seed!(seedVal) # do not reseed rng

    theta = 2.0*deflection* pi * rand() # Rotation about the pole (Z)
    phi = 2.0* pi * rand() # For direction of pole deflection
    z = 2.0*deflection * rand() # For magnitude of pole deflection

    # Compute a vector V used for distributing points over the sphere
    # via the reflection I - V Transpose(V).  This formulation of V
    # will guarantee that if x[1] and x[2] are uniformly distributed,
    # the reflected points will be uniform on the sphere.  Note that V
    # has length sqrt(2) to eliminate the 2 in the Householder matrix.

    r =  sqrt(z)
    Vx = sin(phi) * r
    Vy = cos(phi) * r
    Vz = sqrt(2.0 - z)

    # Compute the row vector S = Transpose(V) * R, where R is a simple
    # rotation by theta about the z-axis.  No need to compute Sz since
    # it's just Vz.

    st = sin(theta)
    ct = cos(theta)
    Sx = Vx * ct - Vy * st
    Sy = Vx * st + Vy * ct

    # Construct the rotation matrix  ( V Transpose(V) - I ) R, which
    # is equivalent to V S - R.

    M = zeros(3,3)
    M[1,1] = Vx * Sx - ct;
    M[1,2] = Vx * Sy - st;
    M[1,3] = Vx * Vz;

    M[2,1] = Vy * Sx + st;
    M[2,2] = Vy * Sy - ct;
    M[2,3] = Vy * Vz;

    M[3,1] = Vz * Sx;
    M[3,2] = Vz * Sy;
    M[3,3] = 1.0 - z;   # This equals Vz * Vz - 1.0

    return M
end

max_stars_repo_name: tsipenyuk/FPRS

max_issues_repo_name: tsipenyuk/FPRS

max_forks_repo_name: tsipenyuk/FPRS


text:
# See https://en.wikipedia.org/wiki/Derivative
# Examples from https://github.com/JuliaSymbolics/Symbolics.jl

using Symbolics

@variables x y z

expressions = [
    :(Symbolics.derivative(x + x^2, x))
    :(Symbolics.derivative(x + 2x^2 + 6x^3 - 5x^4, x))
    :(Symbolics.derivative(cos(x) * sin(x), x))
    :(Symbolics.derivative(exp(x), x))
    :(Symbolics.derivative(log(x), x))
    :(Symbolics.jacobian([x + x*y, x^4 + y], [x, y]))
    :(Symbolics.hessian(x^4 + z, [x, y, z]))
]

for expr in expressions
    println("[] :julia-answer \"\"\"", escape_string(string(expr)), " = ", eval(expr), "\"\"\".")
end

max_stars_repo_name: pmoura/eye

max_issues_repo_name: pmoura/eye

max_forks_repo_name: pmoura/eye


text:
const OCTAVIAN_NUM_TASKS = Ref(1)
_nthreads() = OCTAVIAN_NUM_TASKS[]

@generated function calc_factors(::Union{Val{nc},StaticInt{nc}} = num_cores()) where {nc}
    t = Expr(:tuple)
    for i ∈ nc:-1:1
        d, r = divrem(nc, i)
        iszero(r) && push!(t.args, (i, d))
    end
    t
end
# const CORE_FACTORS = calc_factors()

MᵣW_mul_factor(::True) = StaticInt{4}()
MᵣW_mul_factor(::False) = StaticInt{9}()
MᵣW_mul_factor() = MᵣW_mul_factor(has_feature(Val(:x86_64_avx512f)))

W₁Default(::True) = StaticFloat64{0.006089395198610773}()
W₂Default(::True) = StaticFloat64{0.7979822724696168}()
R₁Default(::True) = StaticFloat64{0.5900561503730485}()
R₂Default(::True) = StaticFloat64{0.762152930709678}()

W₁Default_arch(::Val{:znver1}) = StaticFloat64{0.053918949422353986}()
W₂Default_arch(::Val{:znver1}) = StaticFloat64{0.3013238122374886}()
R₁Default_arch(::Val{:znver1}) = StaticFloat64{0.6077103834481342}()
R₂Default_arch(::Val{:znver1}) = StaticFloat64{0.8775382433240162}()

W₁Default_arch(::Union{Val{:znver2},Val{:znver3}}) = StaticFloat64{0.1}()
W₂Default_arch(::Union{Val{:znver2},Val{:znver3}}) = StaticFloat64{0.993489411720157}()
R₁Default_arch(::Union{Val{:znver2},Val{:znver3}}) = StaticFloat64{0.6052218809954467}()
R₂Default_arch(::Union{Val{:znver2},Val{:znver3}}) = StaticFloat64{0.7594052633561165}()

W₁Default_arch(_) = StaticFloat64{0.1}()
W₂Default_arch(_) = StaticFloat64{0.15989396641218157}()
R₁Default_arch(_) = StaticFloat64{0.4203583148344484}()
R₂Default_arch(_) = StaticFloat64{0.8775382433240162}()

W₁Default(::False) = W₁Default_arch(VectorizationBase.cpu_name())
W₂Default(::False) = W₂Default_arch(VectorizationBase.cpu_name())
R₁Default(::False) = R₁Default_arch(VectorizationBase.cpu_name())
R₂Default(::False) = R₂Default_arch(VectorizationBase.cpu_name())

W₁Default() = W₁Default(has_feature(Val(:x86_64_avx512f)))
W₂Default() = W₂Default(has_feature(Val(:x86_64_avx512f)))
R₁Default() = R₁Default(has_feature(Val(:x86_64_avx512f)))
R₂Default() = R₂Default(has_feature(Val(:x86_64_avx512f)))

_first_cache(::StaticInt{1}) = StaticInt{1}()
_first_cache(::StaticInt) = StaticInt{2}()
first_cache() = _first_cache(VectorizationBase.num_l2cache())
second_cache() = first_cache() + One()

_first_cache_size(fcs::StaticInt) = ifelse(eq(first_cache(), StaticInt(2)) & cache_inclusive(StaticInt(2)), fcs - cache_size(One()), fcs)
_first_cache_size(::Nothing) = StaticInt(262144)
first_cache_size() = _first_cache_size(cache_size(first_cache()))

_second_cache_size(scs::StaticInt, ::True) = scs - cache_size(first_cache())
_second_cache_size(scs::StaticInt, ::False) = scs
_second_cache_size(::StaticInt{0}, ::Nothing) = StaticInt(3145728)
function second_cache_size()
    sc = second_cache()
    _second_cache_size(cache_size(sc), cache_inclusive(sc))
end

first_cache_size(::Type{T}) where {T} = first_cache_size() ÷ static_sizeof(T)
second_cache_size(::Type{T}) where {T} = second_cache_size() ÷ static_sizeof(T)

bcache_count() = VectorizationBase.num_cache(second_cache())

const BCACHEPTR = Ref{Ptr{Cvoid}}(C_NULL)
const BCACHE_LOCK = Threads.Atomic{UInt}(zero(UInt))

if Sys.WORD_SIZE ≤ 32
    const ACACHEPTR = Ref{Ptr{Cvoid}}(C_NULL)
end


text:
m6 = let
    dat = MixedModels.dataset(:mrk17_exp1)
    form = @formula(1000 / rt ~ 1 + F*P*Q*lQ*lT + (1+F+P+Q+lQ+lT|subj) + (1+P+Q+lQ+lT|item))
    contr = Dict(
        :F => EffectsCoding(),
        :P => EffectsCoding(),
        :Q => EffectsCoding(),
        :lQ => EffectsCoding(),
        :lT => EffectsCoding(),
        :subj => Grouping(),
        :item => Grouping(),
    )
    fit(MixedModel, form, dat; contrasts = contr)
end


text:
"""
Read CSV-formatted posterior output.

# API

See [`StanCSV.read_chain`](@ref) and [`StanCSV.matching_files`], neither is exported.

# File format

Posterior samples are supposed to be available in a CSV file which follows the following
conventions:

1. Whitespace and content after `#` is ignored. No fields are quoted or escaped. If a line
has a `#`, it cannot have content.

2. The first non-ignored line contains a comma-separated list of variable names, of the
format `var[.i1.i2.…]` where optional indexes follow dots. Indexes for the same variable are
contiguous, and have a column-major layout.

3. Subsequent lines contain the same number of comma-separated floating-point values.

`cmdstan` outputs posterior samples in this format (diagnostics are ignored for now), hence
the name.
"""
module StanCSV

using ArgCheck: @argcheck
using DocStringExtensions: FUNCTIONNAME, SIGNATURES
import ..Chains

####
#### Header parsing building blocks
####

"""
$(SIGNATURES)

Parse a Stan variable name as `name::Symbol => indexes::Tuple`.
"""
function parse_variable_name(variable_name::AbstractString)
    parts = split(variable_name, '.')
    @argcheck !isempty(first(parts))
    name = Symbol(strip(first(parts)))
    index = tuple(parse.(Ref(Int), parts[2:end])...)
    name => index
end

"""
$(SIGNATURES)

Given a vector of `name => index` pairs and a starting index `i`, which should be `1`s, find
the last index for this variable (ie the dimensions) and check that that intermediate
indices are contiguous (in column major ordering).

Return `name`, `dimensions`, and the position for the continuation.
"""
function collapse_contiguous_dimensions(names_indexes, i)
    name, index = names_indexes[i]
    @argcheck all(index .== 1) "Indexes for $(name) don't start with ones."
    n = length(names_indexes)
    j = i + 1
    while j ≤ n
        if names_indexes[j][1] ≠ name
            break
        else
            j += 1
        end
    end
    dimensions = names_indexes[j - 1][2]
    indexes = CartesianIndices(dimensions)
    for (k, ix) in zip(i:(j-1), CartesianIndices(dimensions))
        index_k = names_indexes[k][2]
        @argcheck CartesianIndex(index_k) == ix "Non-contiguous index $(k) for $(name)."
    end
    name, dimensions, j
end

"""
$(SIGNATURES)

Parse a vector of `names => indexes` to a column schema.
"""
function parse_schema(names_indexes)
    n = length(names_indexes)
    i = 1
    schema = Pair{Symbol,Tuple{Vararg{Int}}}[]
    while i ≤ n
        name, dimensions, next_i = collapse_contiguous_dimensions(names_indexes, i)
        @argcheck all(s -> first(schema) ≠ name, schema) "Duplicate name $(name) in column $(i)."
        push!(schema, name => dimensions)
        i = next_i
    end
    Chains.IndexSchema(NamedTuple{Tuple(first.(schema))}(last.(schema)))
end

####
#### CSV file reading.
####

"""
$(SIGNATURES)

Test if `line` is a comment line.
"""
is_comment_line(line) = occursin(r"^ *#", line)

"""
$(SIGNATURES)

Find the first non-comment line and read it as a `Chains.IndexSchema`. When there is no
such line, throw an `EOFError`.
"""
function read_schema(io::IO)
    while true                  # will throw an EOFError if not fount
        line = readline(io; keep = false)
        if !is_comment_line(line)
            return parse_schema(parse_variable_name.(split(line, ',')))
        end
    end
end

"""
$(SIGNATURES)

Read `nfields` fields from `io`, parse as `Float64`, and push into `buffer`, returning
`true`. Fewer fields than `nfields` results in an error.

If a line starts with a `'#'`, do no parsing, return `false`.
"""
function read_csv_line!(buffer, io::IO, nfields::Integer)
    line = readline(io; keep = false)
    is_comment_line(line) && return false
    pos = 1
    last_pos = lastindex(line)
    for _ in 1:nfields
        @argcheck pos ≤ last_pos "Fewer than $(nfields) fields in line."
        delim_pos = something(findnext(isequal(','), line, pos), last_pos + 1)
        push!(buffer, parse(Float64, SubString(line, pos, delim_pos - 1)))
        pos = delim_pos + 1
    end
    true
end

function read_csv_flat_data!(buffer, io, nfields)
    row_count = 0
    while !eof(io)
        row_count += read_csv_line!(buffer, io, nfields)
    end
    row_count
end

function csv_buffer_to_array(buffer, row_count, dims)
    n = length(dims)
    permutedims(reshape(buffer, dims..., row_count), vcat([1 + n], 1:n))
end

function read_chain(io::IO)
    sch = read_schema(io)
    nfields = length(sch)
    buffer = Vector{Float64}()
    row_count = read_csv_flat_data!(buffer, io, nfields)
    sample_matrix = collect(permutedims(reshape(buffer, nfields, row_count)))
    Chains.Chain(sch, sample_matrix)
end

"""
$(SIGNATURES)

Read data from a CSV file, that uses the output format from Stan.

Return a vector of `name => array` pairs, where `name::Symbol` is the variable name, and the
`array` is the data for that variable.

Consecutive columns for the same variable name with `n`-dimensional indexes are assembled
into an `n+1` dimensional array, with index `var[row, i1, i2, …]` containing the results
for `var.i1.i2.…` in the given `row`.
"""
read_chain(filename::AbstractString) = open(read_chain, filename, "r")

"""
$(SIGNATURES)

Return a vector `id::Int => filename` where `filename` is `prefix_<digits>.csv` and the path
is a file. Prefix should be a path, eg

```julia
$(FUNCTIONNAME)("/tmp/samples_")
```

will return `[1 => "/tmp/samples_1.csv", 2 => "/tmp/samples_2.csv", …]` if the files exist.

Recommended use: obtaining all the filenames in a directory.
"""
function matching_files(prefix::AbstractString)
    dir = dirname(prefix)
    base = basename(prefix)
    pattern = Regex("^" * base * raw"(\d+)\.csv$")
    ids_files = Pair{Int,String}[]
    for filename in readdir(dir)
        m = match(pattern, filename)
        p = joinpath(dir, filename)
        if m ≢ nothing && isfile(p)
            push!(ids_files, parse(Int, m.captures[1]) => p)
        end
    end
    @argcheck allunique(first.(ids_files)) "Non-unique file ids, perhaps because of 0-padding?"
    sort!(ids_files, by = first)
    ids_files
end

end


text:
function read_boarding_passes( filename )

    open(filename) do f
        return readlines(f)
    end

end

tickets = read_boarding_passes( "input05.txt")
tickets .= replace.( tickets,  r"[RB]" => "1")
tickets .= replace.( tickets,  r"[FL]" => "0")
seat_id = parse.(Int, tickets, base=2)

println(maximum(seat_id))
seat = setdiff(minimum(seat_id):maximum(seat_id), seat_id)
println(seat)

max_stars_repo_name: pnavaro/AdventOfCode

max_issues_repo_name: pnavaro/AdventOfCode

max_forks_repo_name: pnavaro/AdventOfCode


text:
module Models

using ImageBase
using ImageBase.ImageCore.MappedArrays: of_eltype
using ImageBase.FiniteDiff

# Introduced in ColorVectorSpace v0.9.3
# https://github.com/JuliaGraphics/ColorVectorSpace.jl/pull/172
using ImageBase.ImageCore.ColorVectorSpace.Future: abs2

"""
This submodule provides predefined image-related models and its solvers that can be reused
by many image processing tasks.

- solve the Rudin Osher Fatemi (ROF) model using the primal-dual method: [`solve_ROF_PD`](@ref) and [`solve_ROF_PD!`](@ref)
"""
Models

export solve_ROF_PD, solve_ROF_PD!

##### implementation details

"""
    solve_ROF_PD([T], img::AbstractArray, λ; kwargs...)

Return a smoothed version of `img`, using Rudin-Osher-Fatemi (ROF) filtering, more commonly
known as Total Variation (TV) denoising or TV regularization. This algorithm is based on the
primal-dual method.

This function applies to generic N-dimensional colorant array and is also CUDA-compatible.
See also [`solve_ROF_PD!`](@ref) for the in-place version.

# Arguments

- `T`: the output element type. By default it is `float32(eltype(img))`.
- `img`: the input image, usually a noisy image.
- `λ`: the regularization coefficient. Larger `λ` results in more smoothing.

# Parameters

- `num_iters::Int`: The number of iterations before stopping.

# Examples

```julia
using ImageFiltering
using ImageFiltering.Models: solve_ROF_PD
using ImageQualityIndexes
using TestImages

img_ori = float.(testimage("cameraman"))
img_noisy = img_ori .+ 0.1 .* randn(size(img_ori))
assess_psnr(img_noisy, img_ori) # ~20 dB

img_smoothed = solve_ROF_PD(img_noisy, 0.015, 50)
assess_psnr(img_smoothed, img_ori) # ~27 dB

# larger λ produces over-smoothed result
img_smoothed = solve_ROF_PD(img_noisy, 5, 50)
assess_psnr(img_smoothed, img_ori) # ~21 dB
```

# Extended help

Mathematically, this function solves the following ROF model using the primal-dual method:

```math
\\min_u \\lVert u - g \\rVert^2 + \\lambda\\lvert\\nabla u\\rvert
```

# References

- [1] Chambolle, A. (2004). "An algorithm for total variation minimization and applications". _Journal of Mathematical Imaging and Vision_. 20: 89–97
- [2] [Wikipedia: Total Variation Denoising](https://en.wikipedia.org/wiki/Total_variation_denoising)
"""
solve_ROF_PD(img::AbstractArray{T}, args...) where T = solve_ROF_PD(float32(T), img, args...)
function solve_ROF_PD(::Type{T}, img::AbstractArray, args...) where T
    u = similar(img, T)
    buffer = preallocate_solve_ROF_PD(T, img)
    solve_ROF_PD!(u, buffer, img, args...)
end

# non-exported helper
preallocate_solve_ROF_PD(img::AbstractArray{T}) where T = preallocate_solve_ROF_PD(float32(T), img)
function preallocate_solve_ROF_PD(::Type{T}, img) where T
    div_p = similar(img, T)
    p = ntuple(i->similar(img, T), ndims(img))
    ∇u = ntuple(i->similar(img, T), ndims(img))
    ∇u_mag = similar(img, eltype(T))
    return div_p, p, ∇u, ∇u_mag
end

"""
    solve_ROF_PD!(out, buffer, img, λ, num_iters)

The in-place version of [`solve_ROF_PD`](@ref).

It is not uncommon to use ROF solver in a higher-level loop, in which case it makes sense to
preallocate the output and intermediate arrays to make it faster.

!!! note "Buffer"
    The content and meaning of `buffer` might change without any notice if the internal
    implementation is changed. Use `preallocate_solve_ROF_PD` helper function to avoid
    potential changes.

# Examples

```julia
using ImageFiltering.Models: preallocate_solve_ROF_PD

out = similar(img)
buffer = preallocate_solve_ROF_PD(img)
solve_ROF_PD!(out, buffer, img, 0.2, 30)
```

"""
function solve_ROF_PD!(
        out::AbstractArray{T},
        buffer::Tuple,
        img::AbstractArray,
        λ::Real,
        num_iters::Integer) where T
    # seperate a stub method to reduce latency
    FT = float32(T)
    if FT == T
        solve_ROF_PD!(out, buffer, img, Float32(λ), Int(num_iters))
    else
        solve_ROF_PD!(out, buffer, FT.(img), Float32(λ), Int(num_iters))
    end
end
function solve_ROF_PD!(
        out::AbstractArray,
        (div_p, p, ∇u, ∇u_mag)::Tuple,
        img::AbstractArray,
        λ::Float32,
        num_iters::Int)
    # Total Variation regularized image denoising using the primal dual algorithm
    # Implement according to reference [1]
    τ = 1//4 # see 2nd remark after proof of Theorem 3.1.

    # use the same symbol in the paper
    u, g = out, img

    fgradient!(p, g)
    # This iterates Eq. (9) of [1]
    # TODO(johnnychen94): set better stop criterion
    for _ in 1:num_iters
        fdiv!(div_p, p)
        # multiply term inside ∇ by -λ. Thm. 3.1 relates this to `u` via Eq. 7.
        @. u = g - λ*div_p
        fgradient!(∇u, u)
        _l2norm_vec!(∇u_mag, ∇u) # |∇(g - λdiv p)|
        # Eq. (9): update p
        for i in 1:length(p)
            @. p[i] = (p[i] - (τ/λ)*∇u[i])/(1 + (τ/λ) * ∇u_mag)
        end
    end
    return u
end

function _l2norm_vec!(out, Vs::Tuple)
    all(v->axes(out) == axes(v), Vs) || throw(ArgumentError("All axes of input data should be the same."))
    @. out = abs2(Vs[1])
    for v in Vs[2:end]
        @. out += abs2(v)
    end
    @. out = sqrt(out)
    return out
end

end # module


text: # LDLᵀ factorization for block-symmetric systems