text:
using Statistics
struct FEMGrid{BC<:AbstractBC, S<:AbstractShape, T} <: AbstractGrid{BC, T}
  N::Int
  L::Float64
  bases::Vector{BasisFunction{S, T}}
  partitionunityweights::Vector{T}
  buckets::Dict{BasisFunction{TopHatShape,Float64},Set{Int}}
  function FEMGrid{BC}(N::Int, L::Float64, bases::Vector{BasisFunction{S,T}},
      puw=zeros(T, N)) where {BC<:AbstractBC, S<:AbstractShape, T}
    bc = BC(L)
    puw = 1 ./ [sum(v(centre(u), bc) for u ∈ bases) for v ∈ bases]
    @assert !any(iszero, puw)

    zero!.(bases)

    # Split the grid up into buckets;
    # Any basisfunction that does not intersect with a key
    # will also not intersect with any element in the associated value vector.
    # Therefore to find all bases on the grid intersecting with a basis, B,
    # simply find the key(s), k, that B intersects with and loop through the
    # associated value, which itselfis a vector of grid bases interesecting
    # with the key, k.
    # N.B. TopHatShape is used simply as an interval
    nbins = Int(round(sqrt(N))) # make the buckets approx √N "cells" wide
    bins = collect(0:nbins) ./ nbins * L
    buckets = Dict{BasisFunction{TopHatShape,Float64}, Set{Int}}()
    bc = PeriodicGridBC(L)
    for i ∈ 2:length(bins)
      l, u = bins[i-1], bins[i]
      bin = BasisFunction{TopHatShape,Float64}(TopHatShape(u - l), (l + u)/2, NaN) # weight is unimportant
      contents = Set{Int}()
      for (index, item) in enumerate(bases)
        in(translate(item, bin, bc)...) && push!(contents, index)
      end
      buckets[bin] = contents
    end

    return new{BC, S, T}(N, L, bases, puw, buckets)
  end
end

function FEMGrid(N::Int, L::Float64, shape::S, ::Type{BC}=PeriodicGridBC
    ) where {BC<:AbstractBC, S<:AbstractShape}
  width(shape) > L && @error ArgumentError "Shapes must not be wider than the grid"
  Δ = L / N
  bases = [BasisFunction(shape, (i-0.5) * Δ) for i ∈ 1:N]
  return FEMGrid{BC}(N, L, bases)
end

function (l::FEMGrid{BC})(x) where {BC}
  # This is equivalent to return sum(i(x, BC(0.0, l.L)) * weight(i) for i ∈ l)
  # but is faster
  bc = BC(l.L)
  processedindices = Set{Int}()
  amount = 0.0
  for (k, v) ∈ l.buckets
    in(x, k) || continue
    for index ∈ filter(x->!(x ∈ processedindices), v)
      item = l.bases[index]
      amount += item(x, bc) * weight(item)
      push!(processedindices, index)
    end
  end
  return amount
end
Base.size(f::FEMGrid) = (f.N,)
Base.iterate(f::FEMGrid) = iterate(f.bases)
Base.iterate(f::FEMGrid, state) = iterate(f.bases, state)
Base.getindex(l::FEMGrid, i) = l.bases[i]
bases(l::FEMGrid) = l.bases

function Base.setindex!(l::FEMGrid, v, i::Integer)
  zero!(l.bases[i])
  l.bases[i] += v
  return l
end

lower(l::FEMGrid) = 0.0
upper(l::FEMGrid) = l.L
zero!(f::FEMGrid) = map(zero!, f)
domainsize(l::FEMGrid) = l.L
numberofunknowns(l::FEMGrid) = length(l.bases)

#function solve(l::FEMGrid{BC}, f::F) where {BC, F}
#  bc = BC(lower(l), upper(l))
#  zero!(l)
#  A = zeros(length(l), length(l))
#  for (j, v) ∈ enumerate(l.bases), (i, u) ∈ enumerate(l.bases)
#    A[i, j] = integral(u, v, bc)
#  end
#  b = zeros(length(l))
#  for (j, v) ∈ enumerate(l.bases)
#    b[j] = integral(v, f, bc)
#  end
#  x = A \ b
#  return x
#end

# Pull a special trick to solve circulant matrix that pops out of periodic BCs
function solve(A::AbstractMatrix, b::AbstractVector)
  x = if Circulant(A) ≈ A
    c = fft(A[:, 1])
    c[1] = 1 # 1st mode (eigenvalue) is zero; will * by false later
    o = fft(b) ./ c
    o[1] *= false # PeriodicGridBC means we want zero DC offset
    for i ∈ eachindex(o) # incase we are dividing by any other zeros
      isnan(o[i]) && (o[i] *= false)
    end
    real.(ifft(o))
  elseif rank(Matrix(A)) == size(A, 1) # is non-singular
    demean!(A \ b)
  else
    u, v = eigen(Matrix(A))
    idu = 1 ./ u
    idu[end] *= false # just like ignoring 0th Fourier mode?
    A⁻¹ = real.(v * diagm(idu) * inv(v))
    demean!(A⁻¹ * b)
  end
  @assert all(isfinite, x)
  return x
end

"""
    fastdeposit!(l::FEMGrid{BC},particle)where{BC<:AbstractBC}

description

...
# Arguments
- `l::FEMGrid{BC}`:
- `particlewhere{BC<:AbstractBC}`:
...

# Equivalent to
```julia
function deposit!(l::FEMGrid{BC}, particle) where {BC<:AbstractBC}
  bc = BC(l.L)
  qw = charge(particle) * weight(particle)
  for (index, item) ∈ enumerate(l)
    item += integral(item, basis(particle), bc) * qw *
      l.partitionunityweights[index]
  end
  return l
end
```
"""
function fastdeposit!(l::FEMGrid{BC}, particle) where {BC<:AbstractBC}
  bc = BC(l.L)
  qw = charge(particle) * weight(particle)
  # Only deposit from particle into each basis function once
  # Record which basis functions have recevied deposition in this set
  processedindices = Set{Int}()
  for (k, v) ∈ (l.buckets)
    overlap(k, basis(particle), bc) || continue
    for index ∈ filter(x->!(x ∈ processedindices), v)
      item = l.bases[index]
      item += integral(item, basis(particle), bc) * qw *
        l.partitionunityweights[index]
      push!(processedindices, index)
    end
  end
  return l
end

"""
    fastantideposit(l::FEMGrid{BC},particle)where{BC<:AbstractBC}

description

...
# Arguments
- `l::FEMGrid{BC}`:
- `particlewhere{BC<:AbstractBC}`:
...

# Example
```julia
function antideposit(l::FEMGrid{BC}, particle) where {BC<:AbstractBC}
  bc = BC(l.L)
  amount = 0.0
  for (index, item) ∈ enumerate(l)
    amount += integral(item, basis(particle), bc) * weight(item)
  end
  return amount
end
```
"""
function fastantideposit(l::FEMGrid{BC}, particle) where {BC<:AbstractBC}
  bc = BC(l.L)
  # Only antideposit from each basis function once
  # Record which basis functions have been processed already
  processedindices = Set{Int}()
  amount = 0.0
  for (k, v) ∈ l.buckets
    overlap(k, basis(particle), bc) || continue
    for index ∈ filter(x->!(x ∈ processedindices), v)
      item = l.bases[index]
      amount += integral(item, basis(particle), bc) * weight(item)
      push!(processedindices, index)
    end
  end
  return amount
end

function deposit!(l::FEMGrid{BC}, particle) where {BC<:AbstractBC}
  return fastdeposit!(l, particle)
  #bc = BC(l.L)
  #qw = charge(particle) * weight(particle)
  #for (index, item) ∈ enumerate(l)
  #  item += integral(item, basis(particle), bc) * qw *
  #    l.partitionunityweights[index]
  #end
  #return l
end

function antideposit(l::FEMGrid{BC}, particle) where {BC<:AbstractBC}
  return fastantideposit(l, particle)
  #bc = BC(l.L)
  #amount = 0.0
  #for (index, item) ∈ enumerate(l)
  #  amount += integral(item, basis(particle), bc) * weight(item)
  #end
  #return amount
end

function update!(l::FEMGrid{BC}, f::F) where {BC, F}
  bc = BC(0.0, l.L)

  A = zeros(length(l), length(l))
  ThreadsX.foreach(enumerate(l)) do (j, v)
    for (i, u) in enumerate(l)
      A[i, j] = v(centre(u), bc)
    end
  end

  b = zeros(length(l))
  ThreadsX.foreach(enumerate(l)) do (j, v)
    b[j] = f(centre(v))
  end

  x = A \ b

  @avx for i ∈ eachindex(x)
    setindex!(l, x[i], i)
  end

  return l
end

abstract type AbstractFEMField{BC} <: AbstractField{BC} end

"""
    GalerkinFEMField

Using basis function, u to represent both the charge density the
the electric field, solve for the latter from the former via:
    ∫ u' v'dx = ∫ u' v dx.
"""
struct LSFEMField{BC, S<:AbstractShape, T} <: AbstractFEMField{BC}
  chargedensity::FEMGrid{BC,S,T}
  electricfield::FEMGrid{BC,S,T}
  function LSFEMField(chargedensity::L, electricfield::L
      ) where {BC,S,T,L<:FEMGrid{BC,S,T}}
    @assert chargedensity.N == electricfield.N
    @assert chargedensity.L == electricfield.L
    return new{BC,S,T}(chargedensity, electricfield)
  end
end
LSFEMField(a::FEMGrid) = LSFEMField(a, deepcopy(a))
LSFEMField(N::Int, L::Real, shape::S) where {S<:AbstractShape} = LSFEMField(FEMGrid(N, L, shape))

numberofunknowns(l::AbstractFEMField) = numberofunknowns(l.chargedensity)

zero!(f::AbstractFEMField) = map(zero!, (f.chargedensity, f.electricfield))

"""
    GalerkinFEMField

Using basis function, u, to represent the charge density and v to represent the
the electric field, solve for the latter from the former via:
    -∫ u' v'dx = ∫ u v' dx.
"""
struct GalerkinFEMField{BC, S1<:AbstractShape, S2<:AbstractShape, T
    } <: AbstractFEMField{BC}
  chargedensity::FEMGrid{BC,S1,T}
  electricfield::FEMGrid{BC,S2,T}
  function GalerkinFEMField(chargedensity::FEMGrid{BC,S1,T}, electricfield::FEMGrid{BC,S2,T}
      ) where {BC,S1,S2,T}
    @assert chargedensity.L == electricfield.L
    return new{BC,S1,S2,T}(chargedensity, electricfield)
  end
end
function GalerkinFEMField(N::Int, L::Real, chargeshape::S1, efieldshape::S2
  ) where {S1<:AbstractShape, S2<:AbstractShape}
  return GalerkinFEMField(FEMGrid(N, L, chargeshape), FEMGrid(N,L, efieldshape))
end

function solve(field::AbstractFEMField{PeriodicGridBC})
  A = massmatrix(field)
  b = forcevector(field)
  return solve(A, b)
end

function solve(f::AbstractFEMField{BC}) where {BC}
  A = massmatrix(f)
  b = forcevector(f)
  p = AMGPreconditioner{SmoothedAggregation}(A)
  return IterativeSolvers.gmres(A, b, Pl=p)
end

postsolve!(x, ::Type{AbstractBC}) = x
postsolve!(x, ::Type{PeriodicGridBC}) = demean!(x)

"""
    solve!(f::AbstractFEMField{PeriodicGridBC})

Solve for the electric field from the chargedensity.
The periodic boundary condition creates a circulant Symmetric matrix,
which can be solved via FFTs (n log n).

Other solution methods are eigen decomposition (n^3) and regularisation.

...
# Arguments
- `field::AbstractFEMField{PeriodicGridBC}`:
...

# Example
```julia
```
"""
function solve!(f::AbstractFEMField{BC}) where {BC}
  x = try
    solve(f)
  catch err
    A = massmatrix(f)
    b = forcevector(f)
    @save "solve!_$(typeof(f)).jld2" f err
    @warn "Caught $err, saved field struct, and exiting."
    rethrow(err)
  end
  f.electricfield .= postsolve!(x, BC)
  return f
end

function quadrature(integrand::F, intervals) where {F}
  return QuadGK.quadgk(integrand, intervals..., order=7, rtol=2eps())[1]
end

function lsfemmassmatrixintegral(u::BasisFunction{GaussianShape},
                                 v::BasisFunction{GaussianShape}, cache)
  if !(sigma(shape(u)) == sigma(shape(v)))
    return _lsfemmassmatrixintegral(u, v, cache)
  end
  Δ = abs(centre(v) - centre(u))
  σ = sigma(shape(u))
  return exp(-(Δ/σ)^2 / 2) * (1 - (Δ/σ)^2) / sqrt(2pi) / σ^3
end
function lsfemmassmatrixintegral(u::BasisFunction, v::BasisFunction, cache)
  return _lsfemmassmatrixintegral(u, v, cache)
end
function _lsfemmassmatrixintegral(u, v, cache)
  integrand(x) = ForwardDiff.derivative(u, x) * ForwardDiff.derivative(v, x)
  return get!(()->quadrature(integrand, knots(u, v)), cache, getkey(u, v))
end

function lsfemstiffnessmatrixintegral(u, v, cache)
  integrand(x) = ForwardDiff.derivative(u, x) * v(x)
  return get!(()->quadrature(integrand, knots(u, v)), cache, getkey(u, v))
end

function galerkinmassmatrixintegral(u::BasisFunction,
                                    v::BasisFunction, cache)
  return  _galerkinmassmatrixintegral(u, v, cache)
end
function _galerkinmassmatrixintegral(u, v, cache)
  integrand(x) = -ForwardDiff.derivative(u, x) * ForwardDiff.derivative(v, x)
  return get!(()->quadrature(integrand, knots(u, v)), cache, getkey(u, v))
end

function galerkinstiffnessmatrixintegral(u, v, cache)
  integrand(x) = u(x) * ForwardDiff.derivative(v, x)
  return get!(()->quadrature(integrand, knots(u, v)), cache, getkey(u, v))
end

massmatrixintegral(::LSFEMField) = lsfemmassmatrixintegral
massmatrixintegral(::GalerkinFEMField) = galerkinmassmatrixintegral
stiffnessmatrixintegral(::LSFEMField) = lsfemstiffnessmatrixintegral
stiffnessmatrixintegral(::GalerkinFEMField) = galerkinstiffnessmatrixintegral
massmatrixmatrixtype(::LSFEMField) = Symmetric
massmatrixmatrixtype(_) = SparseMatrixCSC
@memoize function massmatrix(f::AbstractFEMField)
  return matrix(f.chargedensity, f.electricfield, massmatrixintegral(f),
                massmatrixmatrixtype(f))
end
@memoize function normalisedstiffnessmatrix(f::AbstractFEMField)
  return matrix(f.chargedensity, f.chargedensity, stiffnessmatrixintegral(f))
end

function forcevector(f::AbstractFEMField)
  return normalisedstiffnessmatrix(f) * weight.(f.chargedensity)
end

function matrix(a::FEMGrid{BC}, b::FEMGrid{BC}, integral::F,
    ::Type{M}=SparseMatrixCSC) where {BC<:PeriodicGridBC, F, M<:AbstractMatrix}
  bc = BC(lower(a), upper(a))
  As = [spzeros(length(a), length(a)) for _ ∈ 1:Threads.nthreads()]
  caches = [IdDict{UInt64,Any}() for _ ∈ 1:Threads.nthreads()]
  foreach(enumerate(b)) do (j, v̄)
    cache = caches[Threads.threadid()]
    A = As[Threads.threadid()]
    for (i, ū) ∈ enumerate(a)
      M == Symmetric && j > i && (A[i, j] = A[j, i]; continue)
      u, v = translate(ū, v̄, bc)
      u ∈ v || continue
      A[i, j] = integral(u, v, cache)
    end
  end
  return M(sum(As))
end


text:
module BasisFunctionsModule
export PrimitiveGaussianBasisFunction, ContractedGaussianBasisFunction

using ..BaseModule

type PrimitiveGaussianBasisFunction
  center::Position
  exponent::Real
  mqn::MQuantumNumber
end

type ContractedGaussianBasisFunction
  coefficients::Array{Real,1}
  primitiveBFs::Array{PrimitiveGaussianBasisFunction,1}
end

end

ext: jl

max_stars_repo_name: tkelman/QuantumLab

max_issues_repo_name: tkelman/QuantumLab

max_forks_repo_name: tkelman/QuantumLab


text:
const poi_type_dialog = collect(Dialog() for _ in 1:10)
format_experiment(i) = string(i[:name], ": ", i[:description])
format_level(j) = string(j[:factor][:name], ": ", j[:level][:value])
format_run(i) = string(i[:runDate], ", ", join(format_level.(i[:factorLevels]), ", "))
const datafile = joinpath(tempdir(), "ab7bcb9c-a5ad-11e9-2f73-4ded39bd6e4d")

const experiment_menu = RadioMenu(["", " "])
empty!(experiment_menu.options)
experiment_menu.pagesize -= 2
if !isfile(datafile)
    write(datafile, "")
else
    data = read(datafile)
    append!(experiment_menu.options, format_experiment.(data))
    experiment_menu.pagesize = min(length(data), 11)
end

function newpoi(i, id, existing_types)
    @label poitypel
    _poitype = requestᵐ("What POI is in column $i (e.g. nest, feeder, track, etc.)?", poi_type_dialog[i])
    if isempty(_poitype)
        @warn "the POI must have some name"
        @goto poitypel
    end
    poitype = Symbol(_poitype)
    if poitype ∈ existing_types
        @warn "there cannot be two identical POIs in the same run" poitype
        @goto poitypel
    end
    temporal = newinterval(poitype == :track)
    i = if length(calibration_menu.options) == 1
        [1]
    else
        requestᵐ("Which calibration/s calibrates this POI?", calibration_menu)
    end
    poi = POI(calibrations[collect(i)], temporal)
    return poitype, poi
end

function newrun(designations)
    @label desig
    dialog = Dialog()
    designation = request("Give a designation for this new run", dialog)
    if designation ∈ designations
        @warn "run $designation is already taken, try again…"
        @goto desig
    end
    isempty(designation) && @goto desig
    designation
end

function register_pois(ids)
    pois = Dict{Symbol, POI}()
    for (i, id) in enumerate(ids)
        poitype, poi = newpoi(i, id, keys(pois))
        newpoiᵒ[] = "$id.csv" => poitype => poi
        pois[poitype] = poi
    end

    runs = JLSO.load(joinpath(sourcefolder, "runs.jlso"))
    if haskey(runs, "data") && length(runs) == 1
        delete!(runs, "data")
    end
    if isempty(runs)
        @info "registering a new run"
        k = newrun(keys(runs))
        runs[k] = pois
    else
        options = collect(keys(runs))
        pushfirst!(options, "Register a new run")
        run_menu = RadioMenu(options, min(1, length(options)))
        @label whichrun
        i = requestᵐ("Which run includes these POIs?", run_menu)
        if i == 1
            k = newrun(options)
            runs[k] = pois
        else
            k = options[i]
            for (kk, p) in pois
                if haskey(runs[k], kk)
                    @warn "run $k already has POI $kk"
                    @goto whichrun
                end
                runs[k][kk] = p
            end
        end
    end
    JLSO.save(joinpath(sourcefolder, "runs.jlso"), runs)
end

    #=token = gettoken()
    newdata = getexperimentsruns(token)
    olddata = read(datafile)
    data = if olddata ≠ newdata
        @info "updating local list of experiments…"
        empty!(experiment_menu.options)
        j = JSON3.read(newdata)
        newdata = j[:data][:experiments]
        append!(experiment_menu.options, format_experiment.(newdata))
        experiment_menu.pagesize = min(length(newdata), 11)
        newdata
    else
        olddata
    end
    experimenti = requestᵐ("Which experiment/s include these POIs?", experiment_menu)
    # experimentId = data[experimenti][:objectId]

    options = format_run.(data[experimenti][:runs])
    run_menu = RadioMenu(options)
    runi = requestᵐ("Which run/s include these POIs?", menu)
    runId = UUID(data[experimenti][:runs][runi][:objectId])

    AppRun(runID, pois)
end=#


text:
module Pythia

using Base: Integer, Float64, sign_mask
using Statistics, Distributions

export sglasso, glasso, lasso
export MeanForecast, NaiveForecast, SES, Holt, HoltWinters
export fit, predict

include("algorithms/sglasso.jl")
include("algorithms/basicMethods.jl")
include("algorithms/ETS.jl")

end


text:
function step!(hs :: HardSpheres, pc :: PeriodicCollisions)

    n = hs.n
    ϵ = hs.ϵ

    dtmin, fant, i1, i2 = dt_min_position(pc)

    if fant > 0

       qa = [hs.q[i2] + offset[k] + dtmin*hs.v[i2] for k in 1:9]
       qb = hs.q[i1] + dtmin .* hs.v[i1]
       vr = argmin([norm(qx-qb) for qx in qa])

       hs.q[i2] = qa[vr]
       hs.q[i1] = qb

       J = (dot(hs.v[i2]-hs.v[i1],hs.q[i2]-hs.q[i1]))/ 2ϵ
       hs.v[i1] = hs.v[i1] + J * (hs.q[i2]-hs.q[i1]) / 2ϵ
       hs.v[i2] = hs.v[i2] - J * (hs.q[i2]-hs.q[i1]) / 2ϵ
       hs.q[i1] = mod.(hs.q[i1],1)
       hs.q[i2] = mod.(hs.q[i2],1)

    end

    for i in [1:min(i1,i2)-1;min(i1,i2)+1:max(i1,i2)-1;max(i1,i2)+1:n]
        qnew = hs.q[i] .+ dtmin .* hs.v[i]
        hs.q[i] = mod.(qnew, 1.)
    end

    pc.dt .-= dtmin
    reset!( pc, i1)
    reset!( pc, i2)

    compute_dt!(pc, i1, hs)
    compute_dt!(pc, i2, hs)

    return dtmin

end


text:
### A Pluto.jl notebook ###
# v0.15.1

using Markdown
using InteractiveUtils

# ╔═╡ 7a206e18-4562-469d-b580-f2d0cc2f80df
using Pkg; Pkg.develop(path=".."); Pkg.activate("..")

# ╔═╡ 7c28cf4c-cd0c-11eb-0507-b10f2a883e8c
begin
    using RDatasets
    using LinearAlgebra
    using Plots
    using StatsPlots
    plotly(size=(680, 500))
end

# ╔═╡ 93477f7b-aef4-44e5-ba3b-8285a6261500
md"
# K-Medoids

JingYu Ning
"

# ╔═╡ 64f0d05d-9c4e-4935-bf33-a5ee5835f99a
md"
## IRIS dataset
"

# ╔═╡ 7f59bdb6-c6de-4a25-9ff0-8a81457a4242
iris = dataset("datasets", "iris")

# ╔═╡ 63073482-7add-4234-ab34-6c8068b2f5a6
begin
    cols = names(iris)

    function gen_pots()
        plots = Matrix{Plots.Plot}(undef, 4, 4)
        for i in 1:4, j in 1:4
            if i == j
                plots[i, j] = @df(iris, histogram(cols(i)))
            elseif i < j
                plots[i, j] = @df(iris, histogram2d(
                    cols(i), cols(j),
                    nbins=20, color=:coolwarm, colorbar=false
                ))
            else
                plots[i, j] =  @df(iris, scatter(cols(i), cols(j)))
            end

            (i==1) && (plots[i, j] = plot!(ylabel=cols[j])) # add x label
            (j==4) && (plots[i, j] = plot!(xlabel=cols[i])) # add y label
        end

        return plots
    end
end

# ╔═╡ fde3d98a-be02-47fe-a8c1-293092298659
plot(gen_pots()..., layout=(4, 4), legend=false)

# ╔═╡ 30f8c663-e43b-44e2-9c91-b3e8dd760e5f
md"
## Initial center points
"

# ╔═╡ b1b771a2-3e94-4666-a668-01dee57f9832
md"
**Data size and numbers of classes**
"

# ╔═╡ 52e9ada2-cff4-4fe4-abb0-ec0ebc2d3fa1
begin
    n = size(iris, 1)
    n_classes = 3
end;

# ╔═╡ cbb9a475-36d9-4c9f-a12a-d2cc5c339be7
md"
**Initial three center points**
"

# ╔═╡ ad7bfec2-3fe8-428c-ac92-ad56183cbcdd
# cs_i = rand(1:n, n_classes)
cs_i = [77, 90, 13]

# ╔═╡ 2c6511ca-e070-47ac-8f2e-c1308f491f98
let
    plots = gen_pots()
    for i in 2:4, j in 1:(i-1)
        plots[i, j] = @df iris[cs_i, :] scatter!(plots[i, j], cols(i), cols(j))
    end
    plot(plots..., layout=(4, 4), legend=false)
end

# ╔═╡ bb476b47-ba34-41d7-a32f-fc204f4a328c
md"
## Training
"

# ╔═╡ f558898f-d426-4a36-baed-fa4dd07cf7b2
let
	data = randn(2, 1000)
	data_o = hcat(data, [9, -9])

	function d²(data::Matrix, i::Integer, j::Integer)
		return sum(abs.(data[:, j]-data[:, i]))
	end

	xlim=(-10, 10); ylim=(-10, 10)

	pd = scatter(data[1, :], data[2, :], xlim=xlim, ylim=ylim)
	id = argmin([sum(d²(data, i, j) for j in 1:1000) for i in 1:1000])
	scatter!(pd, [data[1, id]], [data[2, id]])

	po = scatter(data_o[1, :], data_o[2, :], xlim=xlim, ylim=ylim)
	io = argmin([sum(d²(data_o, i, j) for j in 1:1000) for i in 1:1000])
	scatter!(po, [data_o[1, io]], [data_o[2, io]])
	plot(pd, po, legend=false)
end

# ╔═╡ e6f38859-fd04-4dd6-84f3-991fdc72317d
md"
**Manhattan distance between two points**
"

# ╔═╡ 9ce95067-b54f-4a9d-8158-8e25b4a2480f
function d(data::DataFrame, p₁_i::Integer, p₂_i::Integer)
	return sum(abs.(Vector(data[p₁_i, 1:4]) - Vector(data[p₂_i, 1:4])))
end

# ╔═╡ 9cf597b7-9ea2-4c5d-b810-83b5d027d758
begin
    signal = "⚡"

    iris[!, :Predict] = zeros(n)
    for _ in 1:100
        iris[!, :Predict] = map(i->argmin([d(iris, i, ci) for ci in cs_i]), 1:n)

		iris[!, :SumInClassD] = [
			sum((iris[i, :Predict]==iris[j, :Predict] && d(iris, i, j)) for j in 1:n)
			for i in 1:n
		]
		for i in 1:n
			if iris[cs_i[iris[i, :Predict]], :SumInClassD] > iris[i, :SumInClassD]
				cs_i[iris[i, :Predict]] = i
			end
		end
    end
end

# ╔═╡ 7b5ea310-f5c2-4f9d-91b0-2ebae7b8639e
md"
**Prediction**
"

# ╔═╡ e0b73f20-18ba-4dff-9b75-0fe82d072167
let
    signal
    plots = gen_pots()
    for i in 2:4, j in 1:(i-1)
        plots[i, j] = @df iris scatter!(plots[i, j], cols(i), cols(j), group=:Predict)
        plots[i, j] = @df iris[cs_i, :] scatter!(plots[i, j], cols(i), cols(j))
    end
    plot(plots..., layout=(4, 4), legend=false)
end

# ╔═╡ 413aa844-5ce1-4349-a715-b62c8db9a1f8
md"
**Ground truth**
"

# ╔═╡ 1ae9328e-72c1-46e4-ba1e-742205c4bd88
let
    plots = gen_pots()
    for i in 2:4, j in 1:(i-1)
        plots[i, j] = @df iris scatter!(plots[i, j], cols(i), cols(j), group=:Species)
    end
    plot(plots..., layout=(4, 4), legend=false)
end

# ╔═╡ Cell order:
# ╟─93477f7b-aef4-44e5-ba3b-8285a6261500
# ╟─7a206e18-4562-469d-b580-f2d0cc2f80df
# ╠═7c28cf4c-cd0c-11eb-0507-b10f2a883e8c
# ╟─64f0d05d-9c4e-4935-bf33-a5ee5835f99a
# ╠═7f59bdb6-c6de-4a25-9ff0-8a81457a4242
# ╟─63073482-7add-4234-ab34-6c8068b2f5a6
# ╟─fde3d98a-be02-47fe-a8c1-293092298659
# ╟─30f8c663-e43b-44e2-9c91-b3e8dd760e5f
# ╟─b1b771a2-3e94-4666-a668-01dee57f9832
# ╠═52e9ada2-cff4-4fe4-abb0-ec0ebc2d3fa1
# ╟─cbb9a475-36d9-4c9f-a12a-d2cc5c339be7
# ╠═ad7bfec2-3fe8-428c-ac92-ad56183cbcdd
# ╟─2c6511ca-e070-47ac-8f2e-c1308f491f98
# ╟─bb476b47-ba34-41d7-a32f-fc204f4a328c
# ╟─f558898f-d426-4a36-baed-fa4dd07cf7b2
# ╟─e6f38859-fd04-4dd6-84f3-991fdc72317d
# ╠═9ce95067-b54f-4a9d-8158-8e25b4a2480f
# ╠═9cf597b7-9ea2-4c5d-b810-83b5d027d758
# ╟─7b5ea310-f5c2-4f9d-91b0-2ebae7b8639e
# ╟─e0b73f20-18ba-4dff-9b75-0fe82d072167
# ╟─413aa844-5ce1-4349-a715-b62c8db9a1f8
# ╟─1ae9328e-72c1-46e4-ba1e-742205c4bd88


text:
struct SOSPolynomialBridge{
    T, F <: MOI.AbstractVectorFunction,
    DT <: SemialgebraicSets.AbstractSemialgebraicSet,
    UMCT <: Union{Vector{<:MOI.ConstraintIndex{MOI.VectorOfVariables}},
                  MOI.ConstraintIndex{MOI.VectorOfVariables}},
    UMST, MCT,
    GB <: Union{
        Vector{<:Vector{<:MultivariateBases.AbstractPolynomialBasis}}, # Symmetry
        Vector{<:MultivariateBases.AbstractPolynomialBasis},           # Sparsity
        MultivariateBases.AbstractPolynomialBasis,                     # No reduction
    },
    ZB <: MultivariateBases.AbstractPolynomialBasis,
    CT <: SOS.Certificate.AbstractIdealCertificate,
    MT <: MP.AbstractMonomial,
    MVT <: AbstractVector{MT}} <: MOIB.Constraint.AbstractBridge

    Q::Union{Vector{Vector{MOI.VariableIndex}},
             Vector{MOI.VariableIndex}} # Vector{Vector{MOI.VariableIndex}} for sparse SOS
    cQ::UMCT
    gram_basis::GB
    zero_constraint::MOI.ConstraintIndex{F, PolyJuMP.ZeroPolynomialSet{DT, ZB, MT, MVT}}
    domain::DT
    monomials::MVT
    certificate::CT
end

function MOI.Bridges.Constraint.bridge_constraint(
    ::Type{SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT}},
    model::MOI.ModelLike, f::MOI.AbstractVectorFunction,
    s::SOS.SOSPolynomialSet{<:SemialgebraicSets.AbstractAlgebraicSet}) where {T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT}

    @assert MOI.output_dimension(f) == length(s.monomials)
    # MOI does not modify the coefficients of the functions so we can modify `p`.
    # without altering `f`.
    # The monomials may be copied by MA however so we need to copy it.
    p = MP.polynomial(MOIU.scalarize(f), copy(s.monomials))
    # As `*(::MOI.ScalarAffineFunction{T}, ::S)` is only defined if `S == T`, we
    # need to call `changecoefficienttype`. This is critical since `T` is
    # `Float64` when used with JuMP and the coefficient type is often `Int` if
    # `set.domain.V` is `FullSpace` or `FixedPolynomialsSet`.
    # FIXME convert needed because the coefficient type of `r` is `Any` otherwise if `domain` is `AlgebraicSet`
    r = SOS.Certificate.reduced_polynomial(s.certificate, p, MP.changecoefficienttype(s.domain, T))
    gram_basis = SOS.Certificate.gram_basis(
        s.certificate,
        SOS.Certificate.with_variables(r, s.domain),
    )
    g, Q, cQ = SOS.add_gram_matrix(model, MCT, gram_basis, T)
    # MOI does not modify the coefficients of the functions so we can modify `r`.
    # without altering `f`.
    q = MA.operate!!(-, r, g)
    set = PolyJuMP.ZeroPolynomialSet(s.domain, SOS.Certificate.zero_basis(s.certificate), MP.monomials(q))
    coefs = MOIU.vectorize(MP.coefficients(q))
    zero_constraint = MOI.add_constraint(model, coefs, set)
    return SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT}(
        Q, cQ, gram_basis, zero_constraint, s.domain, s.monomials, s.certificate)
end

function MOI.supports_constraint(::Type{SOSPolynomialBridge{T}},
                                 F::Type{<:MOI.AbstractVectorFunction},
                                 ::Type{<:SOS.SOSPolynomialSet{<:SemialgebraicSets.AbstractAlgebraicSet}}) where T
    return MOIU.is_coefficient_type(F, T)
end
function MOIB.added_constrained_variable_types(::Type{<:SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT}}) where {T, F, DT, UMCT, UMST, MCT}
    return constrained_variable_types(MCT)
end
function MOIB.added_constraint_types(::Type{SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT}}) where {T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT}
    return [(F, PolyJuMP.ZeroPolynomialSet{DT, ZB, MT, MVT})]
end
function MOIB.Constraint.concrete_bridge_type(
    ::Type{<:SOSPolynomialBridge{T}},
    F::Type{<:MOI.AbstractVectorFunction},
    ::Type{<:SOS.SOSPolynomialSet{DT, MT, MVT, CT}}) where {T, DT<:SemialgebraicSets.AbstractAlgebraicSet, MT, MVT, CT}
    # promotes VectorOfVariables into VectorAffineFunction, it should be enough
    # for most use cases
    G = MOIU.promote_operation(-, T, F, MOI.VectorOfVariables)
    MCT = SOS.matrix_cone_type(CT)
    UMCT = union_constraint_types(MCT)
    UMST = union_set_types(MCT)
    GB = SOS.Certificate.gram_basis_type(CT)
    ZB = SOS.Certificate.zero_basis_type(CT)
    return SOSPolynomialBridge{T, G, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT}
end

# Attributes, Bridge acting as an model
_num_variables(Q::Vector{MOI.VariableIndex}) = length(Q)
_num_variables(Q::Vector{Vector{MOI.VariableIndex}}) = mapreduce(length, +, Q, init = 0)
function MOI.get(bridge::SOSPolynomialBridge, ::MOI.NumberOfVariables)
    _num_variables(bridge.Q)
end
_list_variables(Q::Vector{MOI.VariableIndex}) = Q
_list_variables(Q::Vector{Vector{MOI.VariableIndex}}) = Iterators.flatten(Q)
function MOI.get(bridge::SOSPolynomialBridge, ::MOI.ListOfVariableIndices)
    _list_variables(bridge.Q)
end
_num_constraints(cQ::Vector, ::Type{C}) where C = count(ci -> ci isa C, cQ)
_num_constraints(cQ::C, ::Type{C}) where C = 1
_num_constraints(cQ, ::Type) = 0
function MOI.get(bridge::SOSPolynomialBridge{T, F, DT, UMCT, UMST},
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where {T, F, DT, UMCT, UMST, S<:UMST}
    return _num_constraints(bridge.cQ, MOI.ConstraintIndex{MOI.VectorOfVariables, S})
end
_list_constraints(cQ::Vector, ::Type{C}) where C = filter(ci -> ci isa C, cQ)
_list_constraints(cQ::C, ::Type{C}) where C = [cQ]
_list_constraints(cQ, C::Type) = C[]
function MOI.get(bridge::SOSPolynomialBridge{T, F, DT, UMCT, UMST},
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where {T, F, DT, UMCT, UMST, S<:UMST}
    return _list_constraints(bridge.cQ, MOI.ConstraintIndex{MOI.VectorOfVariables, S})
end
function MOI.get(::SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT},
                 ::MOI.NumberOfConstraints{F, PolyJuMP.ZeroPolynomialSet{DT, ZB, MT, MVT}}) where {
        T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT
    }
    return 1
end
function MOI.get(b::SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT},
                 ::MOI.ListOfConstraintIndices{F, PolyJuMP.ZeroPolynomialSet{DT, ZB, MT, MVT}}) where {
        T, F, DT, UMCT, UMST, MCT, GB, ZB, CT, MT, MVT
    }
    return [b.zero_constraint]
end

# Indices
_delete_variables(model, Q::Vector{MOI.VariableIndex}) = MOI.delete(model, Q)
function _delete_variables(model, Qs::Vector{Vector{MOI.VariableIndex}})
    for Q in Qs
        MOI.delete(model, Q)
    end
end
function MOI.delete(model::MOI.ModelLike, bridge::SOSPolynomialBridge)
    # First delete the constraints in which the Gram matrix appears
    MOI.delete(model, bridge.zero_constraint)
    # Now we delete the Gram matrix
    _delete_variables(model, bridge.Q)
end

# Attributes, Bridge acting as a constraint
function MOI.get(model::MOI.ModelLike,
                 ::MOI.ConstraintSet,
                 bridge::SOSPolynomialBridge)
    return SOS.SOSPolynomialSet(bridge.domain, bridge.monomials, bridge.certificate)
end
function MOI.get(::MOI.ModelLike,
                 ::MOI.ConstraintPrimal,
                 ::SOSPolynomialBridge)
    throw(SOS.ValueNotSupported())
end

function MOI.get(model::MOI.ModelLike, attr::MOI.ConstraintDual,
                 bridge::SOSPolynomialBridge{T}) where T
    dual = MOI.get(model, attr, bridge.zero_constraint)
    set = MOI.get(model, MOI.ConstraintSet(), bridge.zero_constraint)
    μ = MultivariateMoments.measure(dual, set.monomials)
    function reduced(mono)
        p = MP.polynomial(mono, T)
        domain = MP.changecoefficienttype(bridge.domain, T)
        return SOS.Certificate.reduced_polynomial(bridge.certificate, p, domain)
    end
    return [dot(reduced(mono), μ) for mono in bridge.monomials]
end
function MOI.get(model::MOI.ModelLike,
                 attr::MOI.ConstraintDual,
                 bridge::SOSPolynomialBridge{T, <:MOI.AbstractVectorFunction, SemialgebraicSets.FullSpace}) where {T}
    return MOI.get(model, attr, bridge.zero_constraint)
end
function MOI.get(model::MOI.ModelLike, attr::PolyJuMP.MomentsAttribute,
                 bridge::SOSPolynomialBridge)
    return MOI.get(model, attr, bridge.zero_constraint)
end

function MOI.get(::MOI.ModelLike, ::SOS.CertificateBasis,
                 bridge::SOSPolynomialBridge)
    return bridge.gram_basis
end
function _gram(f::Function, Q::Vector{MOI.VariableIndex}, gram_basis, T::Type, MCT)
    return SOS.build_gram_matrix(convert(Vector{T}, f(Q)), gram_basis, MCT, T)
end
function _gram(f::Function, Qs::Vector{Vector{MOI.VariableIndex}}, gram_bases, T::Type, MCT)
    return SOS.build_gram_matrix(gram_bases, MCT, T) do i
        return convert(Vector{T}, f(Qs[i]))
    end
end
function MOI.get(model::MOI.ModelLike,
                 attr::SOS.GramMatrixAttribute,
                 bridge::SOSPolynomialBridge{T, F, DT, UMCT, UMST, MCT}) where {T, F, DT, UMCT, UMST, MCT}
    return _gram(Q -> MOI.get(model, MOI.VariablePrimal(attr.N), Q),
                 bridge.Q, bridge.gram_basis, T::Type, MCT)
end
function MOI.get(model::MOI.ModelLike,
                 attr::SOS.MomentMatrixAttribute,
                 bridge::SOSPolynomialBridge)
    if bridge.cQ isa Vector{<:MOI.ConstraintIndex}
        return SOS.build_moment_matrix(bridge.gram_basis) do i
            MOI.get(model, MOI.ConstraintDual(attr.N), bridge.cQ[i])
        end
    else
        return SOS.build_moment_matrix(
            MOI.get(model, MOI.ConstraintDual(attr.N), bridge.cQ),
            bridge.gram_basis,
        )
    end
end

max_stars_repo_name: blegat/SumOfSquares

max_issues_repo_name: blegat/SumOfSquares

max_forks_repo_name: blegat/SumOfSquares


text:
using Flux, Optimisers, XLA
using Flux.Data.MNIST
using Flux: onehotbatch, crossentropy

images = reduce(hcat, [vec(Float64.(im)) for im in MNIST.images()])
labels = Float64.(onehotbatch(MNIST.labels(), 0:9))

opt = Optimisers.Descent(0.1)

function step(m)
  m̄, = gradient(m) do m
    @show crossentropy(m(images), labels)
  end
  return opt(m, m̄)
end

m = Chain(Dense(28^2, 32, relu), Dense(32, 10), softmax) |> f64

xstep = xla(step)

for i = 1:100
  global m = xstep(m)
end


text:
# adjusters.jl
using Dates: UTD, value

### truncation
Base.trunc(qdt::QDate, ::Type{Year}) = firstdayofyear(qdt)
Base.trunc(qdt::QDate, ::Type{Month}) = firstdayofmonth(qdt)
Base.trunc(qdt::QDate, ::Type{Day}) = qdt

firstdayofyear(qdt::QDate) = QDate(UTD(value(qdt) - dayofyear(qdt) + 1))
function lastdayofyear(qdt::QDate)
    # cqdate = _qref(qdt)
    # return QDate(UTD(value(qdt) + daysinyear(qdt) - cqdate[3]))
    qdinfo = QREF.qref(qdt)
    jdn2qdate(QREF.lastjdninyear(qdinfo))
end
@inline Dates.firstdayofyear(qdt::QDate) = firstdayofyear(qdt)
@inline Dates.lastdayofyear(qdt::QDate) = lastdayofyear(qdt)

firstdayofmonth(qdt::QDate) = QDate(UTD(value(qdt) - day(qdt) + 1))
@inline Dates.firstdayofmonth(qdt::QDate) = firstdayofmonth(qdt)

function lastdayofmonth(qdt::QDate)
    y, m, l, d = yearmonthleapday(qdt)
    return QDate(UTD(value(qdt) + daysinmonth(y, m, l) - d))
end
@inline Dates.lastdayofmonth(qdt::QDate) = lastdayofmonth(qdt)

# Return the next TimeType that falls on dow
ISQDAYOFWEEK = Dict(先勝 => Dates.DateFunction(is先勝, Base.typemin(QDate)),
                    友引 => Dates.DateFunction(is友引, Base.typemin(QDate)),
                    先負 => Dates.DateFunction(is先負, Base.typemin(QDate)),
                    仏滅 => Dates.DateFunction(is仏滅, Base.typemin(QDate)),
                    大安 => Dates.DateFunction(is大安, Base.typemin(QDate)),
                    赤口 => Dates.DateFunction(is赤口, Base.typemin(QDate)))

Dates.tonext(qdt::QDate, dow::Int; same::Bool=false) = Dates.adjust(ISQDAYOFWEEK[dow], same ? qdt : qdt + Day(1), Day(1), 12)

Dates.toprev(qdt::QDate, dow::Int; same::Bool=false) = Dates.adjust(ISQDAYOFWEEK[dow], same ? qdt : qdt + Day(-1), Day(-1), 12)

function Dates.tofirst(qdt::QDate, dow::Int; of::Union{Type{Year},Type{Month}}=Month)
    qdt = of <: Month ? firstdayofmonth(qdt) : firstdayofyear(qdt)
    return Dates.adjust(ISQDAYOFWEEK[dow], qdt, Day(1), 385)
end

function Dates.tolast(qdt::QDate, dow::Int; of::Union{Type{Year},Type{Month}}=Month)
    qdt = of <: Month ? lastdayofmonth(qdt) : lastdayofyear(qdt)
    return Dates.adjust(ISQDAYOFWEEK[dow], qdt, Day(-1), 385)
end


text:
using Distributions
using DataFrames
using Test
using VariantsNMF

function test_convergence(V, β, local_params)
    params = copy(local_params)
    params.β = β
    res = nmf_FI(V, global_params, params)
    issorted(res.metrics.absolute_cost[1:end-global_params.stopping_iter], rev=true)
end

function test_reg_H(V, β, αs, local_params)
    params = copy(local_params)
    params.β = β
    s_H = Array{Float64,1}()
    for α_H in αs
        params.α_H = α_H
        res = nmf_FI(V, global_params, params)
        push!(s_H, VariantsNMF.matr_avg_row_sparseness(res.H))
    end
    issorted(s_H)
end

function test_reg_W(V, β, αs, local_params)
    params = copy(local_params)
    params.β = β
    s_W = Array{Float64,1}()
    for α_W in αs
        params.α_W = α_W
        res = nmf_FI(V, global_params, params)
        push!(s_W, VariantsNMF.matr_avg_col_sparseness(res.W))
    end
    issorted(s_W)
end

function run_FI_tests(global_params, local_params, βs, αs)
    # record test results
    df_test_results = DataFrame(
        β = Float64[],
        noise = String[],
        test_cvg_passed = Bool[],
        test_reg_H_passed = Bool[],
        test_reg_W_passed = Bool[]
    )

    # get simulations of factorisations for test
    noise_Vs, noise_names = get_many_simulated_V()

    for β in βs
        for (V, name) in zip(noise_Vs, noise_names)
            # cvg test: check monotonicity of the cost function
            # reg_H test: check monotonicity of matrix sparseness
            # reg_W test: check monotonicity of matrix sparseness

            if β ≤ 0 && minimum(V) == 0
                cvg_passed = test_convergence(V .+ 1e-4, β, local_params)
                reg_H_passed = test_reg_H(V .+ 1e-4, β, αs, local_params)
                reg_W_passed = test_reg_W(V .+ 1e-4, β, αs, local_params)
            else
                cvg_passed = test_convergence(V, β, local_params)
                reg_H_passed = test_reg_H(V, β, αs, local_params)
                reg_W_passed = test_reg_W(V, β, αs, local_params)
            end
            @test cvg_passed
            push!(df_test_results, [β name cvg_passed reg_H_passed reg_W_passed])
        end
    end
    println(df_test_results)
end

#### parameters
global_params = NMFParams(
    init          = :random,
    dist          = Uniform(0, 1),
    max_iter      = 1_000,
    stopping_crit = :rel_cost,
    stopping_tol  = 1e-6,
    verbose       = 0,
    seed          = 0
)

#### # 1. MM ALGORITHM
####################################################################################################################

@testset "nmf_MM" begin
    βs = [-1, -0.5, 0, 0.5, 1, 1.5, 2, 3]
    αs = [0, 0.1, 1, 10, 100]

    local_params = FIParams(
        scale_W_iter = false,
        scale_W_last = false,
        l₁ratio_H    = 1,
        l₁ratio_W    = 1,
        α_W          = 0,
        α_H          = 0,
        alg          = :mm
    )

    run_FI_tests(global_params, local_params, βs, αs)
end

#### # 2. H ALGORITHM
####################################################################################################################

@testset "nmf_H" begin
    βs = [-1, -0.5, 0, 0.5, 1, 1.5, 2, 3]
    αs = [0, 0.1, 1, 10, 100]

    local_params = FIParams(
        scale_W_iter = false,
        scale_W_last = false,
        l₁ratio_H    = 1,
        l₁ratio_W    = 1,
        α_W          = 0,
        α_H          = 0,
        alg          = :h
    )

    run_FI_tests(global_params, local_params, βs, αs)
end

#### # 3. ME ALGORITHM
####################################################################################################################

@testset "nmf_ME" begin
    βs = [-1, 0, 0.5, 1.5, 2]
    αs = [0, 0.1, 1, 10, 100]

    local_params = FIParams(
        scale_W_iter = false,
        scale_W_last = false,
        l₁ratio_H    = 1,
        l₁ratio_W    = 1,
        α_W          = 0,
        α_H          = 0,
        alg          = :h
    )

    run_FI_tests(global_params, local_params, βs, αs)
end


text:
struct Grid{FT<:Real}
	xmin :: FT
	xmax :: FT
	nx   :: Int
	δx   :: FT
	xc   :: Array{FT}
	xf   :: Array{FT}
	ymin :: FT
	ymax :: FT
	ny   :: Int
	δy   :: FT
	yc   :: Array{FT}
	yf   :: Array{FT}
end

function GenerateGrid(
	FT = Float64;
	size, x, y,
)

	nx,ny = size
	xmin,xmax = x
	ymin,ymax = y

	xv = collect(range(xmin,xmax,length=nx+1)); δx = (xv[end] - xv[1]) / nx
	xf = xv[2:end]
	xc = (xv[1:(end-1)] .+ xv[2:end]) / 2
    yf = collect(range(ymin,ymax,length=ny+1)); δy = (yf[end] - yf[1]) / ny
	yc = (yf[1:(end-1)] .+ yf[2:end]) / 2

	return Grid{FT}(
		xmin,xmax,nx,δx,xc,xf,
		ymin,ymax,ny,δy,yc,yf
	)

end

function show(io::IO, g::Grid{FT}) where FT <: Real
    print(
		io,
		"The Grid{$FT} is defined as follows:\n",
		"                  domain : x ∈ [$(g.xmin), $(g.xmax)], y ∈ [$(g.ymin), $(g.ymin)]\n",
		"      resolution (nx,ny) : ", (g.nx,g.ny), '\n',
		"    grid spacing (δx,δy) : ", (g.δx,g.δy), '\n'
	)
end


text:
using PyPlot
using Triangulation_diagrams

N = [ -1  -1  -1  0   0  1  0
       2   0  -2  0  -2  0  2 ]

T = [ 6  7  1  2  3  4
      7  2  2  3  5  5
      4  4  7  4  4  6 ]

E = [ 7  1  2  3  5  6
      1  2  3  5  6  7 ]

QN = 4

tri = Triangulation(N, T, E, QN)

figure(1)
draw_triangles(tri)
enumerate_triangles(tri)
enumerate_vertices(tri)
enumerate_bdry_edges(tri, 0.1)
mark_first_local_vertex(tri, 0.3)
axis("off")
savefig("example_triangulation.pdf")

max_stars_repo_name: billmclean/ComputationalMathsNotes

max_issues_repo_name: billmclean/ComputationalMathsNotes

max_forks_repo_name: billmclean/ComputationalMathsNotes


text:
using Powseeker
using POMDPs
using POMDPToolbox
using Plots

mdp = PowseekerMDP()
policy = Downhill(mdp)
hr = HistoryRecorder(max_steps=100)
is = SkierState(1, [1200, 0], 0)

hist = simulate(hr, mdp, policy, is)

gr()
plot(mdp)
plot!(hist)
gui()


text:
using JuMP
using Ipopt
using CPLEX
#using Gurobi
using Distributions

include("prob.jl")
include("cut.jl")

function master_problem(; file = "../data/nesta_case24_ieee_rts_nk.m", k = 4, solver = IpoptSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "AC", initial_cut_file = "cascades_300/IEEE_300_2.txt")

    initial_cut_branches = []
    (initial_cut_file != "none") && (initial_cut_branches = readdlm(initial_cut_file, ',', Int64))
    data = PowerModels.parse_file(file)
    buses = Dict([(parse(Int, k), v) for (k, v) in data["bus"]])
    num_buses = length(buses)
    buses = filter((i, bus) -> bus["bus_type"] != 4, buses)
    branches = Dict([(parse(Int, k), v) for (k, v) in data["branch"]])
    branches = filter((i, branch) -> branch["br_status"] == 1 && branch["f_bus"] in keys(buses) && branch["t_bus"] in keys(buses), branches)
    branch_indexes = collect(keys(branches))
    loads = Dict([(parse(Int, k), v) for (k, v) in data["load"]])
    loads = filter((i, load) -> load["status"] == 1 && load["load_bus"] in keys(buses), loads)

    log_p = Dict([(i, log(branches[i]["prob"])) for i in branch_indexes])

    totalload = 0.0
    for (i,load) in loads
        totalload += load["pd"]
    end
    # println("total load: $totalload ...")

    quit()
    #master_problem = Model(solver=GurobiSolver(LogToConsole=0))
    master_problem = Model(solver=CplexSolver())
    @variable(master_problem, 1e-6 <= eta <= 1e6)
    @variable(master_problem, x[i in branch_indexes], Bin)
    @variable(master_problem, p <= 0)
    @variable(master_problem, y <= 1e6)

    @constraint(master_problem, sum(x) == k)
    @constraint(master_problem, p == sum(x[i]*log_p[i] for i in branch_indexes))

    # adding cuts to remove some previously found solutions
    for i in 1:size(initial_cut_branches)[1]
        expr = zero(AffExpr)
        for j in initial_cut_branches[i,:]
            append!(expr, x[j])
        end
        @constraint(master_problem, expr <= k-1)
    end

    @objective(master_problem, Max, y)
    # add outer approximation of constraint y ⩽ p + log η using lazy constraint callback:  y ⩽ p + log η₀ + 1/η₀(η - η₀)
    # outer approximation of y ⩽ f(x) at x = xᵏ : y ⩽ f(xᵏ) + ∇f(xᵏ)⋅(x-xᵏ)

    function outer_approximate(cb)
        eta_val = getvalue(eta)
        p_val = getvalue(p)
        if eta_val > 0.0
            inv_eta_val = 1/eta_val
            @lazyconstraint(cb, y <= p + log(eta_val) + inv_eta_val*(eta - eta_val))
        end
    end

    addlazycallback(master_problem, outer_approximate)

    # anonymous function
    x_val = i -> getvalue(x[i])

    zlb = -1e10
    zub = 1e10
    eps = 1e-6

    time_taken = 0
    tic()
    solve(master_problem)

    iteration = 0
    final_xvals = Dict{Any,Any}
    final_lines = Any[]

    if (model_constructor == DCPPowerModel)
        cut_constructor = "DC"
    end

    while (zub - zlb > 0.01)
        iteration += 1
        current_xvals = Dict(i => x_val(i) for i in branch_indexes)
        current_lines = collect(keys(filter( (i, xval) -> xval > 1e-4, current_xvals)))
        p_val = getvalue(p)
        for index in current_lines
            data["branch"][string(index)]["br_status"] = 0
        end
        pm = model_constructor(data; solver = solver)
        post_pfls(pm)
        status, solve_time = solve(pm)
        for index in current_lines
            data["branch"][string(index)]["br_status"] = 1
        end
        result = PowerModels.build_solution(pm, status, solve_time)
        alpha = get_cut_coefficients(pm, cut_constructor = cut_constructor)
        sub_objective = result["objective"]
        sub_lb = p_val + log(sub_objective)
        if (sub_lb - zlb > eps)
            zlb = sub_lb
            final_xvals = current_xvals
            final_lines = current_lines
        end
        expr = zero(AffExpr)
        for i in branch_indexes
            if (current_xvals[i] > 0.9)
                append!(expr, x[i])
            end
        end
        @constraint(master_problem, expr <= k-1)
        @constraint(master_problem, eta <= sub_objective + sum(alpha[i]*x[i] for i in collect(keys(alpha))))
        solve(master_problem)
        zub = getobjectivevalue(master_problem)
        time_taken += toq()

        min_value = eps + abs(zlb)
        rel_gap = abs(zub-zlb)
        rel_gap = (min_value > eps) ? rel_gap/min_value : rel_gap
        print("ub: $zub, lb: $zlb, rel_gap: $((zub-zlb)/(eps+abs(zlb))), abs_gap: $(zub-zlb), ")
        p_val = 0.0
        for i in final_lines
            p_val += log_p[i]
        end
        print("solution: $(collect(final_lines)), ")
        log_eta = zlb - p_val
        actual_prob = exp(p_val)
        eta_val = exp(log_eta)
        print("prob: $actual_prob, ")
        print("load_shed: $eta_val, ")
        print("time_taken: $time_taken. \n")
        tic()

        if time_taken >= 86400
            break
        end

        if rel_gap <= 0.05 && sign(zlb)*sign(zub) > 0
            break
        end
    end

    p_val = 0.0
    for i in final_lines
        p_val += log_p[i]
    end
    println("###################################")
    println("k ... $k")
    println("solution = x* ... $(collect(final_lines))")
    println("model: $(model_constructor), cut_type: $(cut_constructor)")
    println("total number of iterations ... $iteration")
    println("solution cost = logprob(x*) + log(η(x*)) ... $zlb")
    println("probability = logprob(x*) ... $p_val")
    log_eta = zlb - p_val
    actual_prob = exp(p_val)
    eta_val = exp(log_eta)
    println("probability = prob(x*) ... $actual_prob")
    println("load shed = η* ... $eta_val")
    println("expected load shed = prob(x*)⋅η* ... $(actual_prob*eta_val)")
    time_taken += toq()
    println("time ... $time_taken")
    # solving the AC load shed model on the final set of lines
    for index in final_lines
        data["branch"][string(index)]["br_status"] = 0
    end
    pm = ACPPowerModel(data; solver = IpoptSolver(print_level=0))
    post_pfls(pm)
    status, solve_time = solve(pm)
    result = PowerModels.build_solution(pm, status, solve_time)
    println("full AC load shed ... $(result["objective"])")
    println("expected AC load shed ... $(actual_prob*result["objective"])")
    println("###################################")
    solution = Dict{AbstractString,Any}("k" => k, "lines" => collect(final_lines),
                                        "model" => model_constructor, "time" => time_taken,
                                        "cut" => cut_constructor, "prob" => actual_prob,
                                        "load_shed" => eta_val,
                                        "expected_load_shed" => actual_prob*eta_val,
                                        "ac_load_shed" => result["objective"],
                                        "expected_ac_load_shed" => actual_prob*result["objective"],
                                        "iterations" => iteration)
    push!(solution_vector, solution)

    return num_buses
end

# solution_dict = Dict{AbstractString,Any}()
solution_vector = Any[]

args = Dict{Any,Any}()
println(ARGS[1])
args["file"] = String(ARGS[1])
args["initial_nk_cut_file"] = String(ARGS[2])
args["k"] = parse(Int, ARGS[3])
args["nf"] = parse(Int, ARGS[4]) # 0 - all three, 1 - only nf, 2 - no nf, 3 - only soc (hurricane case/uniform case)
(args["initial_nk_cut_file"] == "none") && (assert(args["nf"] == 3))

if args["nf"] == 0
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = NFPowerModel, cut_constructor = "DC")
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = DCPPowerModel, cut_constructor = "DC")
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "DC")
elseif args["nf"] == 1
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = NFPowerModel, cut_constructor = "DC")
elseif args["nf"] == 3
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "DC", initial_cut_file=args["initial_nk_cut_file"])
else
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = DCPPowerModel, cut_constructor = "DC")
    num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "DC")
end
# num_buses = master_problem(file = args["file"], k = args["k"], solver = CplexSolver(), model_constructor = SOCWRPowerModel, cut_constructor = "AC")

println(solution_vector)

#= original n-k Networks paper output-writes
if args["nf"] == 0
    f = open("./output_files/$(num_buses)_$(args["k"])", "w")
elseif args["nf"] == 1
    f = open("./output_files/$(num_buses)_$(args["k"])_nf", "w")
elseif args["nf"] == 3 # only use for probability distribution study where case hi or ui is appended to file name
    a = args["file"]
    if args["file"] != "../data/nesta_case1354_pegase_nk.m"
        case = split(split(a, "_")[length(split(a, "_"))], ".")[1]
        f = open("./output_files/$(num_buses)_$(args["k"])_soc_$(case)", "w")
    else
        f = open("./output_files/$(num_buses)_$(args["k"])_soc", "w")
    end
else
    f = open("./output_files/$(num_buses)_$(args["k"])_nonf", "w")
end

write(f, "case model k time iterations probability load_shed expected_load_shed ac_load_shed expected_ac_load_shed lines\n")
for solution in solution_vector
    if solution["model"] == PowerModels.GenericPowerModel{StandardNFForm}
        write(f, "$(num_buses) nf $(solution["k"]) $(solution["time"]) $(solution["iterations"]) $(solution["prob"]) $(solution["load_shed"]) $(solution["expected_load_shed"]) $(solution["ac_load_shed"]) $(solution["expected_ac_load_shed"]) $(solution["lines"]) \n")
    end
    if solution["model"] == PowerModels.GenericPowerModel{PowerModels.StandardDCPForm}
        write(f, "$(num_buses) dc $(solution["k"]) $(solution["time"]) $(solution["iterations"]) $(solution["prob"]) $(solution["load_shed"]) $(solution["expected_load_shed"]) $(solution["ac_load_shed"]) $(solution["expected_ac_load_shed"]) $(solution["lines"]) \n")
    end
    if solution["model"] == PowerModels.GenericPowerModel{PowerModels.SOCWRForm}
        write(f, "$(num_buses) soc $(solution["k"]) $(solution["time"]) $(solution["iterations"]) $(solution["prob"]) $(solution["load_shed"]) $(solution["expected_load_shed"]) $(solution["ac_load_shed"]) $(solution["expected_ac_load_shed"]) $(solution["lines"]) \n")
    end
end

close(f)
=#

#= Commented out before n-k networks paper final runs
solution_dict["solution"] = solution_vector
json_string = JSON.json(solution_dict)
if endswith(args["file"], "api.m")
    write("./json_files/stoch_$(args["k"])_$(num_buses)_api.json", json_string)
elseif endswith(args["file"], "sad.m")
    write("./json_files/stoch_$(args["k"])_$(num_buses)_sad.json", json_string)
else
    write("./json_files/stoch_$(args["k"])_$(num_buses).json", json_string)
end
=#

max_stars_repo_name: kaarthiksundar/acopf-n-k

max_issues_repo_name: kaarthiksundar/acopf-n-k

max_forks_repo_name: kaarthiksundar/acopf-n-k


text:
using FFTW

export OneDPoisson

"""
    OneDPoisson

Derived type to solve the Poisson equation on 1d regular
cartesian mesh with periodic boundary conditions
- kx   : wave number in x
- nc_x : cells number in x
- dx   : x step size
- rht  : fft(rho)
"""
struct OneDPoisson

    grid::OneDGrid
    modes::Array{Float64,1}
    rht::Array{ComplexF64,1}

    function OneDPoisson(grid::OneDGrid)

        nx = grid.nx
        rht = zeros(ComplexF64, nx)
        k =  2π / grid.dimx
        modes = [1.0; k .* vcat(1:nx÷2-1,-nx÷2:-1)...]

        new(grid, modes, rht)

    end

end

export compute_e_from_rho!

"""
    compute_e_from_rho!( ex, poisson, rho )

computes electric field `e` from `rho` by solving Poisson equation.

"""
function compute_e_from_rho!(ex, poisson::OneDPoisson, rho)

   poisson.rht .= fft(rho) ./ poisson.modes
   ex .= real(ifft(-1im .* poisson.rht ))

end

export OneDPoissonPIC

"""
    OneDPoissonPIC( poisson, kernel )

- kernel : Kernel smoother taking care of charge deposition and field evaluation
- poisson : Poisson solver
- rho : Coefficients of expansion of rho (MPI global version)
- efield : Coefficients of expansion of electric field
- phi : Coefficients of expansion of potential
"""
struct OneDPoissonPIC

    ndofs::Int
    kernel::ParticleMeshCoupling2D
    poisson::OneDPoisson
    rho::Vector{Float64}
    efield::Vector{Vector{Float64}}
    phi::Vector{Float64}

    function OneDPoissonPIC(poisson::OneDPoisson, kernel::ParticleMeshCoupling1D)

        ndofs = poisson.grid.nx

        rho = zeros(ndofs)
        efield = [zeros(ndofs), zeros(ndofs)]

        new(ndofs, kernel, poisson, rho, efield)

    end

end

"""
    add_charge!(pic, position, marker_charge)

Add charge from one particle
- self : Pic Poisson solver object
- position : Position of the particle
- marker_charge : Particle weight times charge
"""
function add_charge!(pic::OneDPoissonPIC, position::Float64, marker_charge)

    GEMPIC.add_charge!(pic.rho, pic.kernel, position, marker_charge)

end

"""
    evaluate_rho!(pic, position)

Evaluate charge density at rho at one position
- self : Pic Poisson solver object
- position : Position of the particle
- func_value : Value of rho at given position
"""
function evaluate_rho!(pic :: OneDPoissonPIC, position)

    GEMPIC.evaluate!(kernel, position, pic.rho)

end

export OneDSplittingOperator

"""
Operator splitting type for 2d2v Vlasov-Poisson
- pic :: PIC poisson solver
- pg :: Particle group
"""
struct OneDSplittingOperator

    pic::OneDPoissonPIC
    pg::ParticleGroup{1,1}

end

"""
    operator_t(split, dt)
Push x
- split :: time splitting object
- dt   :: time step
"""
function operator_t!(split::OneDSplittingOperator, dt)

    for i_part = 1:split.pg.n_particles
        split.pg.array[1, i_part] += dt * split.pg.array[2, i_part]
    end

end

function operator_v!(split::OneDSplittingOperator, dt)

    qm = split.pg.q_over_m

    for i_part = 1:split.pg.n_particles

        # Evaluate efields at particle position

        xi = split.pg.array[1, i_part]

        ex = GEMPIC.evaluate(split.pic.kernel, xi, split.pic.efield)

        vx_new = split.pg.array[1, i_part]

        vx_new += dt * qm * ex

        split.pg.array[2, i_part] = vx_new

    end

end

export charge_deposition!

function charge_deposition!(split::OneDSplittingOperator)

    fill!(split.pic.rho_dofs, 0.0)
    for i_part = 1:split.pg.n_particles
        xp = split.pg.array[1, i_part]
        wp = split.pg.array[3, i_part]
        GEMPIC.add_charge!(split.pic.rho_dofs, split.pic.kernel, xp, wp)
    end

end

export solve_fields!

"""
    solve_fields!( pic )

Solve efields from rho
"""
function solve_fields!(pic :: OneDPoissonPIC)

    compute_e_from_rho!(pic.efield[1], pic.poisson, pic.rho)

end

"""
Solve Poisson's equation for the electric field
"""
solve_fields!(split::OneDSplittingOperator) = solve_fields!(split.pic)

export strang_splitting!

"""
    strang_splitting!(split, dt)

Strang splitting
- split :: time splitting object
- dt   :: time step
"""
function strang_splitting!(split::OneDSplittingOperator, dt)

    operator_t!(split, 0.5dt)
    charge_deposition!(split)
    solve_fields!(split)
    operator_v!(split, dt)
    operator_t!(split, 0.5dt)

end

export compute_field_energy

"""
    compute_field_energy

Compute the squared l2 norm of electric field
"""
compute_field_energy(pic::OneDPoissonPIC, component::Int) =
    sum(pic.efield[component] .^ 2) * pic.poisson.grid.dx


text:
function linear_pp(;nᵣ::Int64=3,
        kwargs...)
    linear_pp(nᵣ; kwargs...)
end

"""
    For the random demand, the idea is following:
    1) get random intensity for different products between 0 and 1
    2) calcultae resource size of each product.
    3) multiply selling_period*intensity*resources -> expected number of resources to be requested
    4) scale back the intensity so that the actual number fits the total expected resources expected_res
    For the approximation to make sense, at any timestep, the combined arrival of
    request should be low, i.e. <<0.5
"""
function rand_demand(P::AbstractArray{<:Product}, expected_res::Number, rnd::AbstractRNG)
    prod_selling_period_end = [p.selling_period_end for p in P]
    prod_resources = [sum(p) for p in P]
    prod_demand_intensity = [rand(rnd) for i in 1:length(P)]
    prod_resource_requests = prod_selling_period_end .* prod_demand_intensity .* prod_resources
    resource_requests = sum(prod_resource_requests)
    prod_probs = prod_demand_intensity .* (expected_res / resource_requests)

    # display(prod_probs)
    @assert sum(prod_probs)<0.5 "Product probs $prod_probs are too high with given T=$(maximum(prod_selling_period_end)) and expected resource requests $expected_res"

    D = StaggeredBernoulliScheme(prod_selling_period_end, prod_probs)
end

"""
For the actions, we want to cover the range of the budgets with some space around
"""
function action_space(P::AbstractArray{<:Product}, res_budget_μ::Number)
    prod_resources = [sum(p) for p in P]
    min_b = minimum(prod_resources)*res_budget_μ
    max_b = maximum(prod_resources)*res_budget_μ

    min_a = min_b - res_budget_μ
    max_a = max_b + res_budget_μ
    step_a = res_budget_μ/2
    A = collect(min_a:step_a:max_a)
end

function linear_pp(nᵣ::Int64;
        c::Int64=3,
        T::Int64=10,
        expected_res::Float64=3.,
        res_budget_μ::Float64=5.,
        objective::Symbol=:revenue)

    """
    Initial capacity is the same for every resource
    """
    c₀ = [c for v in 1:nᵣ]

    """
    Resource selling period ends should end one by one at each timestep.
    """
    resource_selling_period_ends = [T for i in 1:nᵣ]
    for i::Int64 in 1:ceil(nᵣ/2)
        resource_selling_period_ends[1:(end-i)].-=1
    end

    """
    Products
    """
    P = PMDPs.create_continuous_linear_products(resource_selling_period_ends)

    """
    Demand
    """
    rnd = Xorshift1024Plus()
    D = rand_demand(P, expected_res, rnd)

    """
    For the budgets, we assume normal distribution per resource.
    """
    B = PMDPs.normal_budgets_per_resource(P, res_budget_μ, res_budget_μ/2)

    """
    Actions
    """
    A = action_space(P, res_budget_μ)

    PMDPs.PMDPProblem(P, c₀, D, B, A, objective)
end


text:
ls = `ls -la`

# run print stdout directory, returns nothing
run(ls)

# read capture
re = read(ls)
println(length(re))

less = `less`
open(less, "w", stdout) do io
   for i = 1:3
       println(io, i)
   end
end

max_stars_repo_name: mkanenobu/trashbox

max_issues_repo_name: mkanenobu/trashbox

max_forks_repo_name: mkanenobu/trashbox


text:
mutable struct SingleStepBuffer{O,A} <: AbstractBuffer
    transition::Union{Nothing,Transition{O,A}}

    SingleStepBuffer{O,A}() where {O,A} = new(nothing)
end

function add!(buffer::SingleStepBuffer{O,A}, obs::O, action::A, reward::Float64) where {O,A}
    buffer.transition = Transition(obs, action, reward)
end

ready(buffer::SingleStepBuffer) = !isnothing(buffer.transition)

reset!(buffer::SingleStepBuffer) = buffer.transition = nothing

function Base.iterate(buffer::SingleStepBuffer, state=1)
    state != 1 && return nothing

    return buffer.transition, 0
end

Base.length(buffer::SingleStepBuffer) = isnothing(buffer.transition) ? 0 : 1

Base.eltype(buffer::SingleStepBuffer{O,A}) where {O,A} = Transition{O,A}


text:
export iauPmat06
"""
Precession matrix (including frame bias) from GCRS to a specified
date, IAU 2006 model.

This function is part of the International Astronomical Union's
SOFA (Standards Of Fundamental Astronomy) software collection.

Status:  support function.

Given:
   date1,date2  double          TT as a 2-part Julian Date (Note 1)

Returned:
   rbp          double[3][3]    bias-precession matrix (Note 2)

Notes:

   1. The TT date date1+date2 is a Julian Date, apportioned in any
   convenient way between the two arguments.  For example,
   JD(TT)=2450123.7 could be expressed in any of these ways,
   among others:

         date1          date2

      2450123.7           0.0       (JD method)
      2451545.0       -1421.3       (J2000 method)
      2400000.5       50123.2       (MJD method)
      2450123.5           0.2       (date & time method)

   The JD method is the most natural and convenient to use in
   cases where the loss of several decimal digits of resolution
   is acceptable.  The J2000 method is best matched to the way
   the argument is handled internally and will deliver the
   optimum resolution.  The MJD method and the date & time methods
   are both good compromises between resolution and convenience.

   2. The matrix operates in the sense V(date) = rbp * V(GCRS), where
   the p-vector V(GCRS) is with respect to the Geocentric Celestial
   Reference System (IAU, 2000) and the p-vector V(date) is with
   respect to the mean equatorial triad of the given date.

Called:
   iauPfw06     bias-precession F-W angles, IAU 2006
   iauFw2m      F-W angles to r-matrix

References:

   Capitaine, N. & Wallace, P.T., 2006, Astron.Astrophys. 450, 855

   Wallace, P.T. & Capitaine, N., 2006, Astron.Astrophys. 459, 981

This revision:  2013 June 18

SOFA release 2018-01-30

Copyright (C) 2018 IAU SOFA Board.  See notes at end.
"""

# void iauPmat06(double date1, double date2, double rbp[3][3])
function iauPmat06(date1::Real, date2::Real)
   rbp = zeros(Float64, 3, 3)

   ccall((:iauPmat06, libsofa_c), Cvoid,
       (Cdouble, Cdouble, Ptr{Cdouble}),
       convert(Float64, date1),
       convert(Float64, date2),
       rbp)

   return SMatrix{3,3}(rbp')
end


text:
module EmpiricalCovarianceOperators
    using LowRankApprox: AbstractLinearOperator, LinearOperator, pheigfact, PartialHermitianEigen
    import Base: convert, size, transpose, show
    using Nullables
    import LinearAlgebra
    using LinearAlgebra: UniformScaling, I, Diagonal
    using Statistics: mean

    include("covarianceoperator.jl")
    include("utilities.jl")

    export CovarianceOperator, CrossCovarianceOperator
    export covmul!, crosscovmul!
    export wbinv

end # module


text:
"""
    add_numerical_flux_internal!(fn, dydt, args...)

Add the numerical flux at the internal faces of the spectral space mesh.

The numerical flux is determined by evaluating

    fn(normal, argvals⁻, argvals⁺)

where:
 - `normal` is the unit normal vector, pointing from the "minus" side to the "plus" side
 - `argvals⁻` is the tuple of values of `args` on the "minus" side of the face
 - `argvals⁺` is the tuple of values of `args` on the "plus" side of the face
and should return the net flux from the "minus" side to the "plus" side.

For consistency, it should satisfy the property that

    fn(normal, argvals⁻, argvals⁺) == -fn(-normal, argvals⁺, argvals⁻)

See also:
- [`CentralNumericalFlux`](@ref)
- [`RusanovNumericalFlux`](@ref)
"""
function add_numerical_flux_internal!(fn, dydt, args...)
    space = axes(dydt)
    Nq = Spaces.Quadratures.degrees_of_freedom(space.quadrature_style)
    topology = space.topology

    for (iface, (elem⁻, face⁻, elem⁺, face⁺, reversed)) in
        enumerate(Topologies.interior_faces(topology))

        internal_surface_geometry_slab =
            slab(space.internal_surface_geometry, iface)

        arg_slabs⁻ = map(arg -> slab(Fields.todata(arg), elem⁻), args)
        arg_slabs⁺ = map(arg -> slab(Fields.todata(arg), elem⁺), args)

        dydt_slab⁻ = slab(Fields.field_values(dydt), elem⁻)
        dydt_slab⁺ = slab(Fields.field_values(dydt), elem⁺)

        for q in 1:Nq
            sgeom⁻ = internal_surface_geometry_slab[q]

            i⁻, j⁻ = Topologies.face_node_index(face⁻, Nq, q, false)
            i⁺, j⁺ = Topologies.face_node_index(face⁺, Nq, q, reversed)

            numflux⁻ = fn(
                sgeom⁻.normal,
                map(
                    slab -> slab isa DataSlab2D ? slab[i⁻, j⁻] : slab,
                    arg_slabs⁻,
                ),
                map(
                    slab -> slab isa DataSlab2D ? slab[i⁺, j⁺] : slab,
                    arg_slabs⁺,
                ),
            )

            dydt_slab⁻[i⁻, j⁻] = dydt_slab⁻[i⁻, j⁻] ⊟ (sgeom⁻.sWJ ⊠ numflux⁻)
            dydt_slab⁺[i⁺, j⁺] = dydt_slab⁺[i⁺, j⁺] ⊞ (sgeom⁻.sWJ ⊠ numflux⁻)
        end
    end
end

"""
    CentralNumericalFlux(fluxfn)

Evaluates the central numerical flux using `fluxfn`.
"""
struct CentralNumericalFlux{F}
    fluxfn::F
end

function (fn::CentralNumericalFlux)(normal, argvals⁻, argvals⁺)
    Favg =
        RecursiveApply.rdiv(fn.fluxfn(argvals⁻...) ⊞ fn.fluxfn(argvals⁺...), 2)
    return RecursiveApply.rmap(f -> f' * normal, Favg)
end

"""
    RusanovNumericalFlux(fluxfn, wavespeedfn)

Evaluates the Rusanov numerical flux using `fluxfn` with wavespeed `wavespeedfn`
"""
struct RusanovNumericalFlux{F, W}
    fluxfn::F
    wavespeedfn::W
end

function (fn::RusanovNumericalFlux)(normal, argvals⁻, argvals⁺)
    y⁻ = argvals⁻[1]
    y⁺ = argvals⁺[1]
    Favg =
        RecursiveApply.rdiv(fn.fluxfn(argvals⁻...) ⊞ fn.fluxfn(argvals⁺...), 2)
    λ = max(fn.wavespeedfn(argvals⁻...), fn.wavespeedfn(argvals⁺...))
    return RecursiveApply.rmap(f -> f' * normal, Favg) ⊞ (λ / 2) ⊠ (y⁻ ⊟ y⁺)
end

function add_numerical_flux_boundary!(fn, dydt, args...)
    space = axes(dydt)
    Nq = Spaces.Quadratures.degrees_of_freedom(space.quadrature_style)
    topology = space.topology

    for (iboundary, boundarytag) in
        enumerate(Topologies.boundary_tags(topology))
        for (iface, (elem⁻, face⁻)) in
            enumerate(Topologies.boundary_faces(topology, boundarytag))
            boundary_surface_geometry_slab =
                surface_geometry_slab =
                    slab(space.boundary_surface_geometries[iboundary], iface)

            arg_slabs⁻ = map(arg -> slab(Fields.todata(arg), elem⁻), args)
            dydt_slab⁻ = slab(Fields.field_values(dydt), elem⁻)
            for q in 1:Nq
                sgeom⁻ = boundary_surface_geometry_slab[q]
                i⁻, j⁻ = Topologies.face_node_index(face⁻, Nq, q, false)
                numflux⁻ = fn(
                    sgeom⁻.normal,
                    map(
                        slab -> slab isa DataSlab2D ? slab[i⁻, j⁻] : slab,
                        arg_slabs⁻,
                    ),
                )
                dydt_slab⁻[i⁻, j⁻] =
                    dydt_slab⁻[i⁻, j⁻] ⊟ (sgeom⁻.sWJ ⊠ numflux⁻)
            end
        end
    end
    return dydt
end


text:
# Autogenerated wrapper script for libwatcher_jll for x86_64-w64-mingw32-cxx11
export libwatcher

using LibUV_jll
JLLWrappers.@generate_wrapper_header("libwatcher")
JLLWrappers.@declare_library_product(libwatcher, "libwatcher.dll")
function __init__()
    JLLWrappers.@generate_init_header(LibUV_jll)
    JLLWrappers.@init_library_product(
        libwatcher,
        "bin\\libwatcher.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
"""
`FDPoisson2DProblem` - Initalize the data structure
"""

function FDPoisson2DProblem(a, b, c, d, n, f)

    zero_function = (x,y)-> 0;
    zero_bc = DirichletBC(zero_function);
    return FDPoisson2DProblem(a, b, c, d, n, n, f, zero_bc, zero_bc, zero_bc, zero_bc);
end

function FDPoisson2DProblem(a, b, c, d, nx, ny, f, bct::TBC, bcb::TBC, bcl::TBC, bcr::TBC) where {TBC<:DirichletBC}
    Δx = (b-a)/(nx+1)
    Δy = (d-c)/(ny+1)
    x = LinRange(a, b, nx+2)[2:end-1];
    y = LinRange(c, d, ny+2)[2:end-1];

    # allocate the matrix and the rhs
    A = spzeros(nx * ny, nx * ny);
    rhs = zeros(nx*ny);

    return FDPoisson2DProblem(a, b, c, d, Δx, nx, Δy, ny, x, y, f, A, rhs, bct, bcb, bcl, bcr)
end


text:
include("utils.jl")

@testset "Positive Numbers" begin
    M = PositiveNumbers()
    @testset "Positive Numbers Basics" begin
        @test repr(M) == "PositiveNumbers()"
        @test repr(PositiveVectors(2)) == "PositiveVectors(2)"
        @test repr(PositiveMatrices(2, 3)) == "PositiveMatrices(2, 3)"
        @test repr(PositiveArrays(2, 3, 4)) == "PositiveArrays(2, 3, 4)"
        @test representation_size(M) == ()
        @test manifold_dimension(M) == 1
        @test !is_point(M, -1.0)
        @test_throws DomainError is_point(M, -1.0, true)
        @test is_vector(M, 1.0, 0.0)
        @test vector_transport_to(M, 1.0, 3.0, 2.0, ParallelTransport()) == 6.0
        @test retract(M, 1.0, 1.0) == exp(M, 1.0, 1.0)
        @test isinf(injectivity_radius(M))
        @test isinf(injectivity_radius(M, -2.0))
        @test isinf(injectivity_radius(M, -2.0, ExponentialRetraction()))
        @test isinf(injectivity_radius(M, ExponentialRetraction()))
        @test project(M, 1.5, 1.0) == 1.0
        @test zero_vector(M, 1.0) == 0.0
        X = similar([1.0])
        zero_vector!(M, X, 1.0)
        @test X == [0.0]
    end
    types = [Float64]
    TEST_FLOAT32 && push!(types, Float32)
    for T in types
        @testset "Type $T" begin
            pts = convert.(Ref(T), [1.0, 4.0, 2.0])
            test_manifold(
                M,
                pts,
                test_forward_diff=false,
                test_reverse_diff=false,
                test_vector_spaces=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                test_vee_hat=false,
                is_mutating=false,
            )
        end
    end
    @testset "Power of Positive Numbers" begin
        M2 = PositiveVectors(2)
        for T in types
            pts2 = [convert.(Ref(T), v) for v in [[1.0, 1.1], [3.0, 3.3], [2.0, 2.2]]]
            test_manifold(
                M2,
                pts2,
                test_forward_diff=false,
                test_reverse_diff=false,
                test_vector_spaces=false,
                test_project_tangent=true,
                test_musical_isomorphisms=true,
                test_default_vector_transport=true,
                test_vee_hat=false,
            )
        end
    end
end


text:
using Telescope
using Test

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


text: print("Hello, World")


text:
function ExtendWaveform(x::AbstractArray, dt::Real, endTime::Real)::AbstractArray
    return vcat(x, [x[end] for i in 1:(endTime/dt-length(x))])
end

function LowPassFilter(x::AbstractArray, RC::Real, dt::Real)::AbstractArray
    alpha = RC/(RC+dt)

    y = zeros(length(x))
    y[1] = alpha*x[1]
    for i = 2:length(x)
        y[i] = y[i-1] + alpha*(x[i]-y[i-1])
    end
    return y
end

function HighPassFilter(x::AbstractArray, RC::Real, dt::Real)::AbstractArray
    alpha = RC/(RC+dt)

    y = zeros(length(x))
    y[1] = x[1]
    for i = 2:length(x)
        y[i] = alpha*(y[i-1] + x[i] - x[i-1])
    end
    return y
end

function CustomFilter(x::AbstractArray, dt::Real)::AbstractArray
    y1 = HighPassFilter(x, 5e-6, dt)
    y2 = LowPassFilter(y1, 7e-9, dt)
    y3 = LowPassFilter(y2, 7e-9, dt)
    #y3 = HighPassFilter(y2, 8e-9, dt)
    #y4 = LowPassFilter(y3, 10e-9, dt)

    return y3
end

function Integrate(x::AbstractArray, dt::Real)::AbstractArray
    y = zeros(length(x))
    for i = 2:length(x)
        y[i] = y[i-1] + x[i]*dt
    end
    return y
end

function TestFilters()
    x = zeros(4000)
    offset = 10
    x[offset:offset+100] .= 1.0
    x[offset:offset+20] .= 1.5

    y = Integrate(x, 0.1)
    wf = CustomFilter(y, 1e-10)

    return x, y, wf
end

max_stars_repo_name: rjtayl/NDP

max_issues_repo_name: rjtayl/NDP

max_forks_repo_name: rjtayl/NDP


text:
function OrientFaceSphere(n,m)
Orient=dot(n,m);
return Orient
end


text:
module TestExamples2DEulerMulticomponent

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "tree_2d_dgsem")

@testset "MHD Multicomponent" begin

  @trixi_testset "elixir_mhdmulti_ec.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [0.04301155595653797, 0.0429973578727621, 0.025745530869947734, 0.16206102676791534, 0.01745438427233918, 0.01745523378100093, 0.026879482381500168, 0.00020380087569638123, 0.012094208262809793, 0.024188416525619587],
      linf = [0.3156206778985403, 0.3094169692980939, 0.2116756351925415, 0.9688251298546025, 0.09076254289155039, 0.09160589769497807, 0.15698032974768816, 0.006131914796912256, 0.0783928755595098, 0.1567857511190196])
  end

  @trixi_testset "elixir_mhdmulti_es.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_es.jl"),
      l2   = [0.042511651022704613, 0.04250604885998708, 0.023846653305641972, 0.11554598944911021, 0.01636256247880442, 0.01636284413777288, 0.025815885276624734, 0.00023041275897137022, 0.010819647132572942, 0.021639294265145884],
      linf = [0.23483836992423804, 0.23500711333514657, 0.11900732096172402, 0.5334761163947968, 0.06168598081763488, 0.06170710018441339, 0.0957218767073531, 0.004187945734585529, 0.06209058954605229, 0.12418117909210458])
  end

  @trixi_testset "elixir_mhdmulti_convergence.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_convergence.jl"),
      l2   = [0.00026086814345575545, 0.0002608681434558476, 0.00048401141419704093, 0.0003246315763252851, 0.0003655995669606298, 0.00036559956696067316, 0.0005007680779079482, 0.0003855610470758256, 5.258783225290134e-5, 0.00010517566450580267, 0.00021035132901160535],
      linf = [0.0015430152095171022, 0.0015430152095179052, 0.0023883101781443854, 0.0010097026785059748, 0.001880514662696009, 0.00188051466269612, 0.0026096395828504393, 0.0010868614907852264, 0.0002056455895989573, 0.0004112911791979146, 0.0008225823583958292])
  end

  @trixi_testset "elixir_mhdmulti_rotor.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_rotor.jl"),
      l2   = [0.7176533615514566, 0.7229955636121396, 0.0, 0.834554078472726, 0.04749874413702431, 0.07911052149635948, 0.0, 0.001984674887414827, 0.19160481586887373, 0.09580240793443687],
      linf = [11.13806940653242, 11.018014654552669, 0.0, 12.446432961705831, 0.6353833312104487, 1.0674585003641777, 0.0, 0.06416407003456842, 3.489953221760654, 1.744976610880327],
      tspan = (0.0, 0.01))

  # TODO: nonconservative terms, remove
  @trixi_testset "elixir_mhdmulti_ec.jl with old nonconservative stuff" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhdmulti_ec.jl"),
      l2   = [0.04301155595653797, 0.0429973578727621, 0.025745530869947734, 0.16206102676791534, 0.01745438427233918, 0.01745523378100093, 0.026879482381500168, 0.00020380087569638123, 0.012094208262809793, 0.024188416525619587],
      linf = [0.3156206778985403, 0.3094169692980939, 0.2116756351925415, 0.9688251298546025, 0.09076254289155039, 0.09160589769497807, 0.15698032974768816, 0.006131914796912256, 0.0783928755595098, 0.1567857511190196],
      volume_flux  = flux_derigs_etal,
      surface_flux = flux_derigs_etal)
  end
end

end

end # module


text:
#!/usr/bin/env julia
# Slow code example
using Statistics

@time begin
  data = Float64[]
  for _ in 1:(10^6)
    group = Float64[]
    for _ in 1:(5 * 10^2)
      push!(group, rand())
    end
    push!(data, mean(group))
  end
  println(
    "98% of the means lie in the estimated range: ",
    (quantile(data, 0.01), quantile(data, 0.99)),
  )
end


text:
module Indicomb

export get_events_catnum_name

using JSON3, HTTP, Dates, SHA
import DataStructures: OrderedDict
import HTTP: escapeuri
import Base.Threads: @threads

_format_dt(dt) = Dates.format(dt, dateformat"y-m-d")

"""
    get_events_catnum_name(baseurl, cat_num, evt_name; params...)

Main API for basic users, given base url (Indico server domain), category number and event title filter,
return a list of (`detail=subcontributions`) events JSON object.

# Example
```
julia> get_events_catnum_name("https://indico.cern.ch", 1X3X, "XXXXX"; from="2021-06-10", to="2021-06-30", apikey=".....", secretkey="....");
JSON3.Object{Vector{UInt8}, SubArray{UInt64, 1, Vector{UInt64}, Tuple{UnitRange{Int64}}, true}} with 29 entries:
  :_type            => "Conference"
  :id               => "10521XX"
  :title            => "XXXXX group meeting"
  :description      => ""
  :startDate        => {…
  :timezone         => "Europe/Zurich"
  :endDate          => {…
  :room             => ""
  :location         => ""
  :address          => ""
  ...
```

"""
function get_events_catnum_name(
    baseurl,
    cat_num,
    evt_name;
    params...,
)

    # get category JSON which has events in simplr form
    cat = get_indico_page(baseurl, "/export/categ/$cat_num.json"; params...)
    simple_evts = filter_events(evt_name, cat)
    res = Vector{JSON3.Object}(undef, length(simple_evts))

    # add event specific headers to into original `params` to get full info
    detail_params = (pretty="yes", detail="subcontributions", params...)
    @threads for idx in eachindex(simple_evts)
        id = simple_evts[idx][:id]
        evt = get_indico_page(baseurl, "/export/event/$id.json"; detail_params...)
        # stupid Indico always return bloated/nested JSON, unpact here!
        res[idx] = only(evt[:results]) #`only` because we're quering a single event
    end
    return res
end
"""
    indico_request(path; params...)

Build a `path?params` url object (without the base domain). The key feature is to construct `signature=` query
based on the SHA1 HMAC algo. The construction of `signature` only cares about path, not base-site domain.
"""
function indico_request(path; params...)
    d = OrderedDict{Symbol,Any}(params) # don't care about performance
    if :secretkey ∈ keys(d)
        # key is `secret_key` and data is entire url except signature=
        _now = Dates.value(now(UTC)) - Dates.UNIXEPOCH
        d[:timestamp] = _now ÷ 1000# ms -> s
        sort!(d)
        temp_url = "$path?$(escapeuri(d))"
        sig = bytes2hex(sha1_hmac(d[:secretkey], temp_url))
        d[:signature] = sig
    end
    return "$path?$(escapeuri(d))"
end

function get_indico_page(baseurl, path; params...)
    _u = "$baseurl$(indico_request(path; params...))"
    return JSON3.read(HTTP.get(_u).body)
end

function filter_events(f::Function, category)
    return filter(f, category[:results])
end

function filter_events(title_name::Union{AbstractString,Regex}, category)
    return filter_events(x -> contains(x[:title], title_name), category)
end

# performe a one-time update and then return digest
"""
    sha1_hmac(key, data)
Return a 1-cycle daupted `digest` given key and data.

# Example
```
julia> bytes2hex(Indicomb.sha1_hmac("123", "Julia"))
"eca18118a715a32ab5f340f57d917c1e2eec96f5"
```
!!! note
    Indico uses *sorted*, query parameters as data. i.e.:

    Given: `"https://cern.ch/export/categ/1135.json?from=2020&apikey=8"`, `data` into this function will be
    `"apikey=8&from=2020"` (string from comes from `HTTP.escapeuri`).
"""
function sha1_hmac(key, data)
    h = HMAC_CTX(SHA1_CTX(), Vector{UInt8}(key))
    update!(h, Vector{UInt8}(data))
    return digest!(h)
end

end


text:
# HellwarthTheory.jl

#####
"""
	HellwarthBScheme(LO)

    Multiple phonon mode reduction to a single effective frequency.
	Hellwarth et al. 1999 PRB, 'B scheme'; the athermal method.
    Averaging procedure is constructed by considering the average effect of the action of multiple branches.

	Follows Eqn (58) in this paper, assuming typo on LHS, should actually be W_e.
"""
function HellwarthBScheme(LO)
    println("Hellwarth B Scheme... (athermal)")
    H58 = sum( (LO[:,2].^2)./ LO[:,1].^2 )
    println("Hellwarth (58) summation: ",H58)

    H59 = sum( LO[:,2].^2 ) # sum of total ir activity squarred
    println("Hellwarth (59) summation (total ir activity ^2): ", H59)
    println("Hellwarth (59) W_e (total ir activity ): ", sqrt(H59))

    omega = sqrt(H59 / H58)
    println("Hellwarth (61) Omega (freq): ",omega)

	return(omega)
end

# More complex scheme, involving thermodynamic Beta
# Hellwarth(50), RHS
"""
	HellwarthAScheme(LO,T=295)

    Multiple phonon mode reduction to a single effective frequency.
	Temperature dependent, defaults to T=295 K.

    Follows Hellwarth et al. 1999 PRB 'A' scheme, Eqn 50 RHS.

	UNTESTED AND UNCERTAIN CODE.
"""
function HellwarthAScheme(LO; T=295)
    println("Hellwarth A scheme...T=$T K")
    omega=LO[:,1] # No unit conversion
    #omega=LO[:,1].*2*pi*1E12 # THz --> Hz
    #omega=LO[:,1].*2*pi*0.02998*1E12 # cm^-1 --> Hz
    β=(omega.*ħ)/(kB*T) #assuming units SI
    H50 = sum( ((LO[:,2].^2).*coth.(β))./omega )
    println("Hellwarth (50) summation: ",H50)

    H51= sum( LO[:,2].^2 ) # sum of total ir activity squarred
    println("Hellwarth (51) summation (total ir activity ^2): ", H51)
    println("Hellwarth (51) W_e (total ir activity ): ", sqrt(H51))

    # OK; so this is deriving Omega / coth(Beta/2)
    omegacoth=H51/H50
    println("omegacoth: ",omegacoth)

    solnfreq=0.0 #required for Julia 0.5 so it realises variable still required for return.
    # Very primitive manner to decouple Omega from both sides of the eqn.
	# Should really rewrite as a bisection (at least!)
    maxfreq=maximum(omega)
    for freq in 0.1:maxfreq/2500:maxfreq
        pseudo_omega=omegacoth/coth(freq * ħ/(2*kB*T))
        if freq>pseudo_omega
            println("freq: $freq pseudo-omega: $pseudo_omega")
            println("freq: ",freq/(2*pi*0.02998*1E12), " pseudo-omega: ",pseudo_omega/(2*pi*0.02998*1E12))
            solnfreq=freq
            break
        end
    end
	return solnfreq
end


text:
@testset "SymbolicExpressionDeriv1" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    env[:x] = 5.6
    expr = symbolic(:(4.0 * x + 1.6))
    y = seval(expr, :x, env, cache)
    @test y == 4.0
end

@testset "SymbolicExpressionDeriv2" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    env[:x] = 5.6
    expr = symbolic(:(-4.0 * x - 1.6 * x))
    y = seval(expr, :x, env, cache)
    @test y == -4-1.6
end

@testset "SymbolicExpressionDeriv3" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    env[:x] = 5.6
    expr = symbolic(:(-4.0 * x * x))
    println(expr)
    y = seval(expr, :x, env, cache)
    @test y == -8 * x
end

@testset "SymbolicExpressionDeriv4" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    y = 10.1
    env[:x] = x
    env[:y] = y
    expr = symbolic(:(-4.0 * x / y))
    res = seval(expr, :x, env, cache)
    @test res == -4 / y
end

@testset "SymbolicExpressionDeriv5" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    y = 10.1
    env[:x] = x
    env[:y] = y
    expr = symbolic(:(-4.0 * x / y))
    res = seval(expr, :y, env, cache)
    @test res == 4 * x / y^2
end

@testset "SymbolicExpressionDeriv6" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    y = 10.1
    env[:x] = x
    env[:y] = y
    expr = symbolic(:(-4 * x * y))
    res = seval(expr, :x, env, cache)
    @test res == -4 * y
end

@testset "SymbolicExpressionDeriv7" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    env[:x] = x
    expr = symbolic(:(-4 * x ^ 3))
    res = seval(expr, :x, env, cache)
    @test res == -12 * x^2
end

@testset "SymbolicExpressionDeriv8" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    y = 10.1
    env[:x] = x
    env[:y] = y
    expr = -log(symbolic(:(x*y))) * symbolic(:y)
    res = seval(expr, :x, env, cache)
    @test res == - y / (x * y) * y
    res = seval(expr, :y, env, cache)
    @test res == - y / (x * y) * x - log(x * y)
    println(cache)
end

@testset "SymbolicExpressionDeriv2_1" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    x = 5.6
    y = 10.1
    env[:x] = x
    env[:y] = y
    expr = symbolic(:(x^2 + 10*x*y + 4*y^2))
    res = seval(expr, (:x,:y), env, cache)
    @test res == 10.0
end

@testset "lam1" begin
    env = SymbolicEnv()
    cache = SymbolicCache()
    @bind env begin
        t = rand()
        lam1 = rand()
        lam2 = rand()
    end
    expr = exp(-lam1*t) - (lam1/lam2) * exp(-lam1*t) * (exp(-lam2*t) - 1)
    res1 = seval(expr, lam1, env, cache)
    @test res1 ≈ seval(-t*exp(-lam1*t) - (1/lam2)*exp(-lam1*t)*(exp(-lam2*t)-1) + t*(lam1/lam2)*exp(-lam1*t)*(exp(-lam2*t)-1), env)
    res2 = seval(expr, lam2, env, cache)
    @test res2 ≈ seval((lam1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) + t*(lam1/lam2)*exp(-lam1*t)*exp(-lam2*t), env)
    res3 = seval(expr, (lam2,lam1), env, cache)
    res4 = seval(expr, (:lam1,:lam2), env, cache)
    @test isapprox(res3, res4)
    @test isapprox(seval((1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) - t*(lam1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) + t*(1/lam2)*exp(-lam1*t)*exp(-lam2*t) - t^2*(lam1/lam2)*exp(-lam1*t)*exp(-lam2*t), env),
                    seval((1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) + t*(1/lam2)*exp(-lam1*t)*exp(-lam2*t) - t*(lam1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) - t^2*(lam1/lam2)*exp(-lam1*t)*exp(-lam2*t), env))
    @test res3 ≈ seval((1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) - t*(lam1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) + t*(1/lam2)*exp(-lam1*t)*exp(-lam2*t) - t^2*(lam1/lam2)*exp(-lam1*t)*exp(-lam2*t), env)
    @test res4 ≈ seval((1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) + t*(1/lam2)*exp(-lam1*t)*exp(-lam2*t) - t*(lam1/lam2^2)*exp(-lam1*t)*(exp(-lam2*t)-1) - t^2*(lam1/lam2)*exp(-lam1*t)*exp(-lam2*t), env)
end


text:
path = Pkg.dir("Dolo")

import Dolo

########
######## RBC with a Markov shock
########

# fn = joinpath(path,"examples","models","LAMP.yaml")
model_mc = Dolo.yaml_import("rbc_mc.yaml")

drc = Dolo.ConstantDecisionRule(model_mc.calibration[:controls])
# @time dr0, drv0 = Dolo.solve_policy(model_mc, drc) #, verbose=true, maxit=10000 )
@time dr = Dolo.time_iteration(model_mc, verbose=false, maxit=10000, details=false)
@time drv = Dolo.evaluate_policy(model_mc, dr, verbose=false) #, verbose=true, maxit=10000)
@time drd = Dolo.time_iteration_direct(model_mc, dr, verbose=false, details=false) #, maxit=500)

sim = Dolo.simulate(model_mc, dr) #; N=100, T=20)

########
######## RBC with an iid shock
########

model = Dolo.Model("rbc_iid.yaml")

@time dr = Dolo.perturbate(model)

drc = Dolo.ConstantDecisionRule(model.calibration[:controls])

@time dr = Dolo.time_iteration(model, maxit=100, verbose=false, details=false)
@time res = Dolo.time_iteration(model, dr; maxit=100, details=true)
#
# @time dr0, drv0 = Dolo.solve_policy(model, drc) #;, verbose=true, maxit=1000 )
# @time res = Dolo.solve_policy(model, drc; details=true) #;, verbose=true, maxit=1000 )

@time drd = Dolo.time_iteration_direct(model, details=false) #, maxit=1000, verbose=true)

@time dr = Dolo.time_iteration_direct(model, drd, details=false) #, maxit=500, verbose=true)
@time res = Dolo.time_iteration_direct(model, drc; details=true)
@time drv = Dolo.evaluate_policy(model, dr; verbose=false)

Dolo.simulate(model, dr)

s0 = model.calibration[:states]+0.1
sim = Dolo.simulate(model, dr, s0)
Dolo.simulate(model, dr; N=10)

res = Dolo.response(model.exogenous, [0.01])

irf = Dolo.response(model, dr, :e_z)
irf[:k]

########
######## RBC with an AR1 driving process
########

# AR1 model: this one should be exactly equivalent to rbc_ar1
model = Dolo.Model("rbc_ar1.yaml")
dp = Dolo.discretize(model.exogenous)

@time drp = Dolo.perturbate(model)

cdr =Dolo.CachedDecisionRule(dr,dp)
@time dr = Dolo.time_iteration(model, cdr, details=false)

@time dr = Dolo.time_iteration(model, details=false)
@time drv = Dolo.evaluate_policy(model, dr, verbose=false, maxit=10000)
@time drd = Dolo.time_iteration_direct(model, dr, details=false) #, verbose=true) #, maxit=500)
#
Dolo.simulate(model, dr, N=10)

Dolo.response(model.exogenous, [0.01])

irf = Dolo.response(model, dr, :z)

irf[:z]


text:
# ------------ Module alias and APIs ------------ #

const Module = MLIR.API.MlirModule

create_empty_module(l::Location) = MLIR.API.mlirModuleCreateEmpty(l)
create_parse(ctx::Context, m::Module) = MLIR.API.mlirModuleCreateParse(ctx, m)
get_context(m::Module) = MLIR.API.mlirModuleGetContext(m)
get_body(m::Module) = MLIR.API.mlirModuleGetBody(m)
is_null(m::Module) = MLIR.API.mlirModuleIsNull(m)
destroy!(m::Module) = MLIR.API.mlirModuleDestroy(m)
get_operation(m::Module) = Operation(MLIR.API.mlirModuleGetOperation(m))

# Constructor.
Module(l::Location) = create_empty_module(l)

@doc(
"""
const Module = MLIR.API.MlirModule
""", Module)


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/optimize/optimize.jl
# Load subproblem type optimization routines and define general optimize!
# throw_optimize_hook!.
#############################################################################

include("optimize_lp_cone.jl")
include("optimize_convex.jl")
include("optimize_nonconvex.jl")

throw_optimize_hook!(m::Optimizer) = optimize_hook!(m.ext_type, m)

function MOI.optimize!(m::Optimizer)

    m._start_time = time()

    # Runs the branch and bound routine
    if !m.enable_optimize_hook

        # Need to eliminate fixed variables prior to checking
        # for convex quadratic constraints
        initial_parse!(m)

        # Determines if the problem is an LP, MILP, SOCP, MISCOP,
        # CONVEX, OF MINCVX PROBLEM TYPE
        parse_classify_problem!(m)

        m._parse_time = m._start_time - time()

        # Throws the problem to the appropriate solution routine
        optimize!(Val{m._working_problem._problem_type}(), m)
    else

        # throws to user-defined optimization hook
        throw_optimize_hook!(m)

    end

    return nothing
end


text:
# this line because the const objects are not being exported
include(joinpath(dirname(@__FILE__), "..", "src/.timeseriesrc.jl"))

facts("const values are set the package defaults") do

    context("SHOWINT") do
        @fact SHOWINT --> true
    end

    context("DECIMALS") do
        @fact DECIMALS --> 4
    end

    context("MISSING") do
        @fact MISSING --> NAN
    end
end

facts("const values are correct") do

    context("NAN") do
        @fact NAN --> "NaN"
    end

    context("NA") do
        @fact NA --> "NA"
    end

    context("BLACKHOLE") do
        @fact BLACKHOLE --> "\u2B24"
    end

    context("DOTCIRCLE") do
        @fact DOTCIRCLE --> "\u25CC"
    end

    context("QUESTION") do
        @fact QUESTION --> "\u003F"
    end
end


text:
module Frames

using ..Bodies:
    ALL_NAMES

using ItemGraphs: ItemGraph, SimpleGraph, add_edge!, items

export AbstractFrame, InertialFrame, RotatingFrame, ICRF, icrf, from_sym, @frame

abstract type AbstractFrame end
abstract type InertialFrame <: AbstractFrame end
abstract type RotatingFrame <: AbstractFrame end

Base.show(io::IO, frame::AbstractFrame) = print(io, string(nameof(typeof(frame))))

const FRAMES = ItemGraph{Symbol}(SimpleGraph())
from_sym(sym::Symbol) = from_sym(Val(sym))

struct ICRF <: InertialFrame end
const icrf = ICRF()
from_sym(::Val{:ICRF}) = icrf

function path_frames(frame1, frame2)
    name1 = nameof(typeof(frame1))
    name2 = nameof(typeof(frame2))
    items(FRAMES, name1, name2)
end

macro frame(name::Symbol, cname::Symbol, parent::Symbol, typ::Symbol)
    quote
        struct $(esc(name)) <: $(esc(typ)) end
        const $(esc(cname)) = $(esc(name))()
        add_edge!(FRAMES, $(Meta.quot(name)), $(Meta.quot(parent)))
        Frames.from_sym(::Val{$(Meta.quot(name))}) = $(esc(cname))
    end
end

include("rotations.jl")
include("iau.jl")
include("iers.jl")

function __init__()
    # IAU frames
    for name in ALL_NAMES
        endswith(name, "Barycenter") && continue

        frame = Symbol("IAU", name)
        cname = Symbol("iau_", lowercase(name))
        add_edge!(FRAMES, frame, :ICRF)
    end

    # IERS frames
    add_edge!(FRAMES, :ICRF, :CIRF)
    add_edge!(FRAMES, :CIRF, :TIRF)
    add_edge!(FRAMES, :TIRF, :ITRF)
end

end