text:
import Base: ==, map, AbstractVector, map, +, *, length, size, zero, iszero
import Statistics: mean

export CustomMean, ZeroMean, ConstantMean, mean

abstract type MeanFunction end
abstract type BaseMeanFunction <: MeanFunction end

"""
    AbstractVector(μ::MeanFunction)

Convert `μ` into an `AbstractVector` if such a representation exists.
"""
function AbstractVector(μ::MeanFunction)
    @assert isfinite(length(μ))
    return map(μ, eachindex(μ))
end

eachindex(μ::BaseMeanFunction) = throw(ErrorException("Cannot construct indices for $μ"))
length(::BaseMeanFunction) = Inf
size(μ::MeanFunction) = (length(μ),)

_map_fallback(f::MeanFunction, X::AV) = [f(x) for x in X]
_map(f::MeanFunction, X::AV) = _map_fallback(f, X)
map(f::MeanFunction, X::AV) = _map(f, X)
map(f::MeanFunction, X::BlockData) = BlockVector([map(f, x) for x in blocks(X)])
map(f::MeanFunction, ::Colon) = map(f, eachindex(f))

"""
    CustomMean <: BaseMeanFunction

A user-defined mean function. `f(x)` should return a scalar for whatever type of `x` this is
intended to work with.
"""
struct CustomMean{T} <: BaseMeanFunction
    f::T
end
@inline (f::CustomMean)(x) = f.f(x)

"""
    ZeroMean <: BaseMeanFunction

Returns zero (of the appropriate type) everywhere.
"""
struct ZeroMean{T<:Real} <: BaseMeanFunction end
@inline (::ZeroMean{T})(x) where T = zero(T)
@inline _map(z::ZeroMean{T}, D::AbstractVector) where T = Zeros{T}(length(D))
==(::ZeroMean, ::ZeroMean) = true

const ZM = ZeroMean{Float64}
zero(::Type{<:MeanFunction}, N::Int) = FiniteZeroMean(1:N)
zero(::Type{<:MeanFunction}) = ZeroMean{Float64}()
zero(μ::MeanFunction) = length(μ) < Inf ? FiniteZeroMean(eachindex(μ)) : ZM()

"""
    ConstantMean{T} <: BaseMeanFunction

Returns `c` (of the appropriate type) everywhere.
"""
struct ConstantMean{T<:Real} <: BaseMeanFunction
    c::T
end
@inline (μ::ConstantMean)(x) = μ.c
@inline _map(μ::ConstantMean, D::AbstractVector) = Fill(μ.c, length(D))
==(μ::ConstantMean, μ′::ConstantMean) = μ.c == μ′.c

+(μ::ConstantMean, μ′::ConstantMean) = ConstantMean(μ.c + μ′.c)
*(μ::ConstantMean, μ′::ConstantMean) = ConstantMean(μ.c * μ′.c)

"""
    EmpiricalMean <: BaseMeanFunction

A finite-dimensional mean function specified by a vector of values `μ`.
"""
struct EmpiricalMean{T<:Real, Tμ<:AbstractVector{T}} <: BaseMeanFunction
    μ::Tμ
    EmpiricalMean(μ::Tμ) where {T<:Real, Tμ<:AbstractVector{T}} = new{T, Tμ}(μ)
end
@inline (μ::EmpiricalMean)(n) = μ.μ[n]
==(μ1::EmpiricalMean, μ2::EmpiricalMean) = μ1.μ == μ2.μ
@inline length(μ::EmpiricalMean) = length(μ.μ)
@inline eachindex(μ::EmpiricalMean) = eachindex(μ.μ)

@inline map(μ::EmpiricalMean, ::Colon) = μ.μ
function _map(μ::EmpiricalMean, X::AV)
    if X == eachindex(μ)
        return μ.μ
    else
        return μ[X]
    end
end
AbstractVector(μ::EmpiricalMean) = μ.μ

+(x::ZeroMean, x′::ZeroMean) = zero(x)
function +(μ::MeanFunction, μ′::MeanFunction)
    @assert size(μ) == size(μ′)
    if iszero(μ)
        return μ′
    elseif iszero(μ′)
        return μ
    else
        return CompositeMean(+, μ, μ′)
    end
end
function *(μ::MeanFunction, μ′::MeanFunction)
    @assert size(μ) == size(μ′)
    return iszero(μ) || iszero(μ′) ? zero(μ) : CompositeMean(*, μ, μ′)
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    Closure()

The transform that applies the closure operation 𝒞
to all rows of the input table. The rows of the
output table sum to one.
"""
struct Closure <: Stateless end

isrevertible(::Type{Closure}) = true

assertions(::Type{Closure}) = [TT.assert_continuous]

function apply(transform::Closure, table)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # original column names
  names = Tables.columnnames(table)

  # table as matrix and get the sum acros dims 2
  X = Tables.matrix(table)
  S = sum(X, dims=2)

  # divides each row by its sum (closure operation)
  Z = X ./ S

  # table with the old columns and the new values
  𝒯 = (; zip(names, eachcol(Z))...)
  newtable = 𝒯 |> Tables.materializer(table)

  newtable, S
end

function revert(::Closure, newtable, cache)
  # transformed column names
  names = Tables.columnnames(newtable)

  # table as matrix
  Z = Tables.matrix(newtable)

  # retrieve cache
  S = cache

  # undo operation
  X = Z .* S

  # table with original columns
  𝒯 = (; zip(names, eachcol(X))...)
  𝒯 |> Tables.materializer(newtable)
end

ext: jl


text:
doc"""
    validate(g, parent_, ij, search_key)
validate `parent_` info for bfs on graph `g` with key `search_key` based on Graph500 result
validation(http://graph500.org/?page_id=12#sec-8)
Return `1` for no error
"""
function validate(
  g::Graph{T},
  parent_::Vector{T},
  ij::Array{T,2},
  search_key::T
  ) where T<:Integer

    out = 1 # Default: no error.

    if parent_[search_key] != search_key
      # root must be the parent of itself
      out = 0
      return out
    end

    N = nv(g)
    slice = findall(x->(x > zero(T)), parent_)

    #Compute levels and check for cycles.
    level_ = zeros(T, length(parent_))
    level_[slice] .= 1
    P = parent_[slice]
    mask = (P .!= search_key)
    k = 0
    while any(mask)
      level_[slice[mask]] .+= 1
      P = parent_[P]
      mask = (P .!= search_key)
      k = k + 1
      if k > N
        #There must be a cycle in the tree.
        out = -3
        return out
      end
    end

    #Check that there are no edges with only one end in the tree.
    #This also checks the component condition.
    lij = level_[ij]
    neither_in = zeros(Bool, size(lij)[2])
    not_neither_in_or_both_in = zeros(Bool, size(lij)[2])
    for i in 1:size(lij)[2]
      if ij[1,i] != ij[2,i]
        neither_in[i] = (lij[1,i] == 0 && lij[2,i] == 0)
        not_neither_in_or_both_in[i] = !(neither_in[i] || (lij[1,i] > 0 && lij[2,i] > 0))
      else
        neither_in[i] = true
        not_neither_in_or_both_in[i] = false
      end
    end
    if any(not_neither_in_or_both_in)
      out = -4
      return out
    end

    #Validate the distances/levels.
    respects_tree_level = trues(1, size(ij)[2])
    not_neither_in_or_respects_tree_level = zeros(Bool,size(ij)[2])
    level_difference = one(T);
    @simd for i in 1:size(ij)[2]
        # To prevent Unsigned integer overflows
        if(lij[1,i] > lij[2,i])
            respects_tree_level[i] = (lij[1,i] .- lij[2,i]) <= level_difference
        else
            respects_tree_level[i] = (lij[2,i] .- lij[1,i]) <= level_difference
        end

        not_neither_in_or_respects_tree_level[i] = !(neither_in[i] || respects_tree_level[i])
    end
    if any(not_neither_in_or_respects_tree_level)
      out = -5
      return out
    end

    return out
  end


text:
# ---
# title: 902. Numbers At Most N Given Digit Set
# id: problem902
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: Math, Dynamic Programming
# link: <https://leetcode.com/problems/numbers-at-most-n-given-digit-set/description/>
# hidden: true
# ---
#
# Given an array of `digits` which is sorted in **non-decreasing** order. You
# can write numbers using each `digits[i]` as many times as we want. For
# example, if `digits = ['1','3','5']`, we may write numbers such as `'13'`,
# `'551'`, and `'1351315'`.
#
# Return _the number of positive integers that can be generated_ that are less
# than or equal to a given integer `n`.
#
#
#
# **Example 1:**
#
#
#
#     Input: digits = ["1","3","5","7"], n = 100
#     Output: 20
#     Explanation:
#     The 20 numbers that can be written are:
#     1, 3, 5, 7, 11, 13, 15, 17, 31, 33, 35, 37, 51, 53, 55, 57, 71, 73, 75, 77.
#
#
# **Example 2:**
#
#
#
#     Input: digits = ["1","4","9"], n = 1000000000
#     Output: 29523
#     Explanation:
#     We can write 3 one digit numbers, 9 two digit numbers, 27 three digit numbers,
#     81 four digit numbers, 243 five digit numbers, 729 six digit numbers,
#     2187 seven digit numbers, 6561 eight digit numbers, and 19683 nine digit numbers.
#     In total, this is 29523 integers that can be written using the digits array.
#
#
# **Example 3:**
#
#
#
#     Input: digits = ["7"], n = 8
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= digits.length <= 9`
#   * `digits[i].length == 1`
#   * `digits[i]` is a digit from `'1'` to `'9'`.
#   * All the values in `digits` are **unique**.
#   * `digits` is sorted in  **non-decreasing** order.
#   * `1 <= n <= 109`
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
#=
if occursin(r"RvSpectMLEcoSystem$", pwd())   cd("SunAsAStar")   end
if occursin(r"SunAsAStar", pwd())   cd("test")   end
 using Pkg
 Pkg.activate(".")
=#

using SunAsAStar
using Test
#using CSV
using DataFrames

make_plots = false
if make_plots
   using Plots
end

@testset "Compare to Collier-Cameron table" begin

   HARPSN_df = SunAsAStar.CSV.read("data/Sun_harpsn_qualflag.csv",DataFrame)
   jd = HARPSN_df.JD .+ 2400000
   HARPSN_df.Delta = HARPSN_df.FWHMobs.^2 .- HARPSN_df.FWHMsid.^2

   # Test scalar
   # Equatorial
   @test SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs(jd[1], obs=:HARPSN) ≈ -0.5969944883008037
   # Ecliptic
   # @test SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs(jd[1], obs=:HARPSN) ≈ -0.5967647214702758

   # Test on array of times
   Delta = SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs.(jd[1:50:end],obs=:HARPSN)

#Delta = SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs(jd,obs=:HARPSN)
#jdsmall = copy(jd[1:50:end])
#Delta = SolarRotation.CalculateFWHMDifference_SolarRotation_from_obs(jdsmall,obs=:HARPSN)

   if make_plots
      using Plots
      plot(jd,HARPSN_df.Delta, label="ACC")
      scatter!(jd[1:50:end],Delta, ms=2,label="Julia")
      xlabel!("Time (d)")
      ylabel!("FWHM_obs^2 - FWHM_sid^2")

   end

   if make_plots
      plot(jd[1:50:end],Delta.-HARPSN_df.Delta[1:50:end], label="Julia-ACC")
      xlabel!("Time (d)")
      ylabel!("Residuals in Δ(FWHM^2)")
   end

   if make_plots && false
      plot(jd[1:50:end],Delta./HARPSN_df.Delta[1:50:end], label="Julia/ACC")
      xlabel!("Time (d)")
      ylabel!("Ratio of Δ(FWHM^2)'s")
   end
   @test maximum(abs.(HARPSN_df.Delta[1:50:end].-Delta)) < 0.01
end


text:
"""
    $(SIGNATURES)
Structure for a graph summary.
"""
struct AbstractDFGSummary
  variables::Dict{Symbol, DFGVariableSummary}
  factors::Dict{Symbol, DFGFactorSummary}
  userId::String
  robotId::String
  sessionId::String
end


text:
const Arr = Union{Vector, Matrix}

type SVM
    X::Matrix
    y::Vector
    C::Float64
    tol::Float64
    max_iter::Integer
    kernel::String
    degree::Integer
    gamma::Float64
    alpha::Vector
    b::Float64
    sv_indx::Vector
    K::Matrix
end

function svm(X::Matrix,
             y::Vector;
             C::Float64 = 1.0,
             kernel::String = "linear",
             max_iter::Integer = 100,
             tol::Float64 = 1e-3,
             degree::Integer = 2,
             gamma::Float64 = 0.1,
             alpha::Vector = zeros(10),
             b::Float64 = 0.0)
    n = size(X,1)
    alpha = zeros(n)
    K = zeros(n,n)
    sv_indx = collect(1:n)
    return SVM(X,y,C,tol,max_iter,kernel, degree,gamma,alpha,b,sv_indx,K)
end

function predict(model::SVM,
                 x::Arr)
    n = size(x,1)
    res = zeros(n)
    if n == 1
        res[1] = predict_row(x,model)
    else
        for i = 1:n
            res[i] = predict_row(x[i,:],model)
        end
    end
    return res
end

function train!(model::SVM)
    n_sample = size(model.X,1)
    model.K = zeros(n_sample,n_sample)
    for i in 1:n_sample
        model.K[:,i] = kernel_c(model.X,model.X[i,:],model)
    end
    # start training

    iters = 0
    while iters < model.max_iter
        iters += 1
       # println("Processing $(iters)/$(model.max_iter)")
        alpha_prev = copy(model.alpha)
        for j = 1:n_sample
            i = rand(1:n_sample)
            eta = 2.0 * model.K[i, j] - model.K[i, i] - model.K[j, j]
            if eta >= 0
                continue
            end
            L, H = count_bounds(i, j,model)

            # Error for current examples
            e_i, e_j = error_(i,model), error_(j,model)

            # Save old alphas
            alpha_io, alpha_jo = model.alpha[i], model.alpha[j]

            # Update alpha
            model.alpha[j] -= (model.y[j] * (e_i - e_j)) / eta
            model.alpha[j] = clamp(model.alpha[j], L, H)

            model.alpha[i] = model.alpha[i] + model.y[i] * model.y[j] * (alpha_jo - model.alpha[j])

            # Find intercept
            b1 = model.b - e_i - model.y[i] * (model.alpha[i] - alpha_jo) * model.K[i, i] -
                 model.y[j] * (model.alpha[j] - alpha_jo) * model.K[i, j]
            b2 = model.b - e_j - model.y[j] * (model.alpha[j] - alpha_jo) * model.K[j, j] -
                 model.y[i] * (model.alpha[i] - alpha_io) * model.K[i, j]
            if 0 < model.alpha[i] < model.C
                model.b = b1
            elseif 0 < model.alpha[j] < model.C
                model.b = b2
            else
                model.b = 0.5 * (b1 + b2)
            end

            # Check convergence
            diff = norm(model.alpha - alpha_prev)
            if diff < model.tol
                break
            end
        end
    end

    #println("Convergence has reached after $(iters). for $(model.kernel)")

    # Save support vectors index
    model.sv_indx = find(model.alpha .> 0)

end

function kernel_c(X::Matrix,
                y::Vector,
                model::SVM)
    if model.kernel == "linear"
        return X * y
    elseif model.kernel == "poly"
        return (X * y).^model.degree
    elseif model.kernel == "rbf"
        n = size(X,1)
        res = zeros(n)
        for i = 1:n
            res[i] = e^(-model.gamma*sumabs2(X[i,:]-y))
        end
        return res
    end
end

function count_bounds(i,j,model)
    if model.y[i] != model.y[j]
        L = max(0, model.alpha[j] - model.alpha[i])
        H = min(model.C, model.C - model.alpha[i] + model.alpha[j])
    else
        L = max(0, model.alpha[i] + model.alpha[j] - model.C)
        H = min(model.C, model.alpha[i] + model.alpha[j])
    end
    return L, H
end

function predict_row(x,model)
    res = kernel_c(model.X,x,model)
    return sign(res' * (model.alpha .* model.y) + model.b)[1]
end

function error_(i,model)
    return predict_row(model.X[i,:],model) - model.y[i]
end

function test_svm()
    X_train, X_test, y_train, y_test = make_cla(n_features = 14)
    predictions = 0
    for kernel in ["linear", "rbf"]
        model = svm(X_train, y_train, max_iter=500, kernel=kernel, C=0.6)
        train!(model)
        predictions = predict(model,X_test)
        println("Classification accuracy $(kernel): $(accuracy(y_test, predictions))")

    end
    #PCA
    pca_model = PCA()
    train!(pca_model, X_test)
    plot_in_2d(pca_model, X_test, predictions, "svm")
end


text:
@testset "Sticky ZigZag 1d" begin
    Random.seed!(1)
    t0 = 0.0
    x0 = [1.0]
    θ0 = [0.8]
    σ = sqrt(0.5)
    μ = 0.9
    ∇ϕ(x, i) = (x[] .- μ)/σ^2

    Γ = sparse(fill(1.0, 1, 1))
    c = 20ones(1)

    Z = ZigZag(Γ, x0*0)
    T = 1000.0

    #sspdmp(∇ϕ, t0, x0, θ0, T, c, F::ZigZag, κ, args...; strong_upperbounds = false,  factor=1.5, adapt=false)
    κ = [0.7]
    trace, _, acc = sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ)
    @show acc[1]/acc[2]
    dt = 0.2
    global ts, xs = sep(collect(discretize(trace, dt)))

    w = sqrt(2π)*σ/(sqrt(2π)*σ + exp(-0.5*μ^2/σ^2)/κ[]) # compute weight

    @test abs(mean(getindex.(xs) .!= 0) - w) < 2.5/sqrt(T) # P(X = 0) = w
    @test abs(mean(getindex.(xs)) - w*μ) < 2/sqrt(T) # E X = wμ
    @test abs(mean(getindex.(xs).^2) - w*(σ^2 + μ^2)) < 2.5/sqrt(T) # E X^2 = w(σ^2 + μ^2)
end

@testset "Sticky SZigZag" begin
    global Γ
    d = size(Γ, 1)

    ∇ϕ(x, i, Γ) = ZigZagBoomerang.idot(Γ, i, x) # sparse computation
    ∇ϕmoving(t, x, θ, i, t′, F, Γ) = ZigZagBoomerang.idot_moving!(Γ, i, t, x, θ, t′, F) # sparse computation

    t0 = 0.0
    x0 = rand(d)
    κ = 1000.0*ones(d) # dont stop, actually
    θ0 = rand([-1.0,-0.5,0.5,1.0], d)

    c = .7*[norm(Γ[:, i], 2) for i in 1:d]

    Z = ZigZag(0.9Γ, x0*0)
    T = 1000.0

    trace, _, acc = @time sspdmp(∇ϕ, t0, x0, θ0, T, c, Z, κ, Γ)
    dt = 0.5
    ts, xs = sep(collect(discretize(trace, dt)))

    @test mean(abs.(mean(xs))) < 2/sqrt(T)
    @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 2.5/sqrt(T)
end

@testset "Sticky Boomerang" begin
    global Γ
    d = size(Γ, 1)

    μ = rand(d)
    ∇ϕ!(y, x, Γ, μ) = mul!(y, Γ, x-μ)

    t0 = 0.0
    x0 = rand(d)
    κ = 1000.0*ones(d) # dont stop, actually
    θ0 = rand([-1.0,-0.5,0.5,1.0], d)

    c = 10.0

    B = Boomerang(sparse(I(d)), μ, 0.5, 0.95, sparse(I(d)))
    T = 1000.0

    trace, _, acc = @time sspdmp(∇ϕ!, t0, x0, θ0, T, c, B, κ, Γ, μ)
    dt = 0.5
    ts, xs = sep(collect(discretize(trace, dt)))

    @test mean(abs.(mean(xs) - μ)) < 2/sqrt(T)
    @test mean(abs.(cov(xs) - inv(Matrix(Γ)))) < 2.5/sqrt(T)
end


text:
"""
    Sanitize a model.

    If your model is stored in `m` and your data are stored in `x1`,
    `x2`, `x3`, etc. then you can sanitize your model with:
    ```julia
    sanitize!(Model(M), Data(x1), Data(x2), Data(x3), ...)
    ```
"""
function sanitize! end

function sanitize!(m::Model{T}; kwargs...)::Model{T} where T
    sanitize!(m, Vector{Data}(undef, 0); kwargs...)
    return m
end

function sanitize!(m::Model{T}, varargs...; kwargs...)::Model{T} where T
    sanitize!(m, convert(Vector{Data}, collect(varargs)); kwargs...)
    return m
end

function sanitize!(m::Model{T}, data::Vector{Data}; kwargs...)::Model{T} where T
    _sanitize!(m, data, _elements(data); kwargs...)
    return m
end

function sanitize!(location::ForceSanitize, varargs...; kwargs...)
    sanitize!(convert(Vector{ForceSanitize}, vcat(location, collect(varargs))); kwargs...)
    return location, varargs...
end

function sanitize!(locations::Vector{ForceSanitize}; kwargs...)
    _sanitize!(locations; kwargs...)
    return locations
end

function _sanitize!(locations::Vector{ForceSanitize}; kwargs...)
    for i = 1:length(locations)
        if isassigned(locations, i)
            _sanitize!(locations[i]; kwargs...)
        end
    end
    return locations
end

function _sanitize!(location::ForceSanitize; kwargs...)
    zero!(location.loc; kwargs...)
    return location
end

function _sanitize!(m::T, data::Vector{Data}, elements::_DataElements; kwargs...)::T where T
    _sanitize_fields!(m, data, elements)
    _sanitize_iterable!(m, data, elements)
    _sanitize_indexable!(m, data, elements)
    return m
end

function _sanitize_fields!(m::T, data::Vector{Data}, elements::_DataElements; kwargs...)::T where T
    for field in fieldnames(T)
        _sanitize!(_get_property(m, field), data, elements)
    end
    return m
end

function _sanitize_iterable!(m::T, data::Vector{Data}, elements::_DataElements; kwargs...)::T where T
    if _is_iterable(T)
        try
            for object in m
                try
                    _sanitize!(object, data, elements)
                catch ex_inner
                    # showerror(stderr, ex)
                    # Base.show_backtrace(stderr, catch_backtrace())
                    # @debug("Ignoring exception [inner]", exception=(ex, catch_backtrace()))
                end
            end
        catch ex_outer
            # showerror(stderr, ex)
            # Base.show_backtrace(stderr, catch_backtrace())
            # @debug("Ignoring exception [outer]", exception=(ex, catch_backtrace()))
        end
    end
    return m
end

function _sanitize_indexable!(m::T, data::Vector{Data}, elements::_DataElements; kwargs...)::T where T
    if _has_isassigned(T)
        _sanitize_indexable_with_check_assigned!(m, data, elements)
    else
        _sanitize_indexable_without_check_assigned!(m, data, elements)
    end
    return m
end

function _sanitize_indexable_with_check_assigned!(m::T, data::Vector{Data}, elements::_DataElements; kwargs...)::T where T
    if _is_indexable(T)
        for i = 1:length(m)
            if isassigned(m, i)
                _sanitize!(m[i], data, elements)
            end
        end
    end
    return m
end

function _sanitize_indexable_without_check_assigned!(m::T, data::Vector{Data}, elements::_DataElements; kwargs...) where T
    if _is_indexable(T)
        for i = 1:length(m)
            try
                _sanitize!(m[i], data, elements)
            catch
            end
        end
    end
    return m
end


text:
using ApproxFun, MultivariateOrthogonalPolynomials, Plots
 import MultivariateOrthogonalPolynomials: DirichletTriangle, Vec

# Neumann
d = Triangle()
S = DirichletTriangle{1,1,1}(d)
Dx = Derivative(S,[1,0])
Dy = Derivative(S,[0,1])
B₁ = I : S → Legendre(Vec(0,0) .. Vec(0,1))
B₂ = I : S → Legendre(Vec(0,0) .. Vec(1,0))
B₃ = I : S → Legendre(Vec(0,1) .. Vec(1,0))

L = [0 B₁ 0;
 0  0  B₂;
 0  B₃ B₃;
 Dx -eye(S) 0;
  Dy 0 -eye(S);
   0.1eye(S) Dx Dy]

f = Fun((x,y) -> exp(x*cos(y)), d)

u = \(L , [0;0; 0; 0;0; f]; tolerance=1E-5)
plot(u[1])

d = Triangle() , Triangle(Vec(1,1),Vec(1,0),Vec(0,1))
S₁,S₂ = S = DirichletTriangle{1,1,1}.(d)
Dx₁ = Derivative(S₁,[1,0])
Dy₁ = Derivative(S₁,[0,1])
Dx₂ = Derivative(S₂,[1,0])
Dy₂ = Derivative(S₂,[0,1])
B₁₁ = I : S[1] → Legendre(Vec(0,0) .. Vec(0,1))
B₁₂ = I : S[1] → Legendre(Vec(0,0) .. Vec(1,0))
B₁₃ = I : S[1] → Legendre(Vec(0,1) .. Vec(1,0))
B₂₁ = I : S[2] → Legendre(Vec(1,1) .. Vec(0,1))
B₂₂ = I : S[2] → Legendre(Vec(1,1) .. Vec(1,0))
B₂₃ = I : S[2] → Legendre(Vec(0,1) .. Vec(1,0))
I₁ = eye(S₁)
I₂ = eye(S₂)

L = [B₁₁  0    0    0    0    0
     B₁₂  0    0    0    0    0
     0    0    0    B₂₁  0    0
     0    0    0    B₂₂  0    0
     B₁₃  0    0   -B₂₃  0    0
     0    B₁₃  B₁₃  0   -B₂₃ -B₂₃
     Dx₁ -I₁   0    0    0    0
     Dy₁  0   -I₁   0    0    0
     0    0    0    Dx₂ -I₂   0
     0    0    0    Dy₂  0   -I₂
     0    Dx₁  Dy₁  0    0    0
     0    0    0    0    Dx₂  Dy₂]

rs = rangespace(L)
f = vcat(map(r->Fun((x,y) -> real(exp(x+im*y)), r,20), rs[1:4])...)

u = \(L, [f; zeros(8)...]; tolerance=1E-5)

plot(u[1])
 plot!(u[4])

plot(u[4])
@which ApproxFun.plotptsvals(u[4])

ncoefficients(u[4])
@which points(u[4].space,210)
points(u[4].space,ncoefficients(u[4]))
pts=points(JacobiTriangle(0,0,0),210)
fromcanonical(S,pts[1])
K = S
n = 210
map(Vec,map(vec,points(ProductTriangle(K),round(Int,sqrt(n)),round(Int,sqrt(n))))...)
S = JacobiTriangle(0,0,0,d[2])
fromcanonical(S, Vec(0.1,0.2))
u[1](0.1,0.2)-real(exp(0.1+0.2im)) # ~ 4.5E-16
u[4](0.6,0.7)-real(exp(0.6+0.7im)) # ~ 6.7E-16


text:
using DZLinearAlgebra: find_saturation_threshold
using MultiFloats
using Printf

const INVERSE_FACTORIALS_F64 = setprecision(() ->
    [Float64x{20}(inv(BigFloat(factorial(BigInt(i))))) for i = 1 : 170],
    BigFloat, 1200)

const LOG2_F64 = setprecision(() ->
    Float64x{20}(log(BigFloat(2))),
    BigFloat, 1200)

function test_multifloat_exp(x::Float64x{N}, num_terms::Int=170,
        reduction_power::Int=0, ::Val{B}=Val(false)) where {N,B}
    exponent_f = Base.rint_llvm(x._limbs[1] / 0.6931471805599453)
    exponent_i = Base.fptosi(Int, exponent_f)
    y = ldexp(Float64x{N}(Float64x{N+!B}(x) -
            exponent_f * Float64x{N+!B}(LOG2_F64)), -reduction_power)
    terms = Vector{Float64x{N}}(undef, num_terms)
    @inbounds terms[1] = y
    for i = 2 : num_terms
        @inbounds terms[i] = terms[div(i, 2)] * terms[div(i+1, 2)]
    end
    exp_y = terms[num_terms] * Float64x{N}(INVERSE_FACTORIALS_F64[num_terms])
    for i = num_terms-1 : -1 : 1
        exp_y += terms[i] * Float64x{N}(INVERSE_FACTORIALS_F64[i])
    end
    exp_y += 1.0
    for _ = 1 : reduction_power
        exp_y *= exp_y
    end
    return ldexp(exp_y, exponent_i)
end

function test_exp_precision(x::Float64x{N}, args...) where {N}
    setprecision(3000) do
        exp_x = test_multifloat_exp(x, args...)
        exp_exact = exp(BigFloat(x))
        -Float64(log2(abs((BigFloat(exp_x) - exp_exact) / exp_exact)))
    end
end

function run_test(num_trials::Int, max_reduction_power::Int, ::Val{B}) where {B}
    for N = 1 : 8
        for reduction_power = 0 : max_reduction_power
            data = [(Float64(num_terms), minimum(test_exp_precision(
                    100 * rand(Float64x{N}) - 50,
                    num_terms, reduction_power, Val{B}())
                for _ = 1 : num_trials)) for num_terms = 1 : 170]
            threshold = round(Int, find_saturation_threshold(data)[1]) + 1
            lost_bits = precision(Float64x{N}) - floor(Int,
                minimum(p[2] for p in data[threshold:end]))
            @printf("Float64x%d (power %2d): %3d terms; %3d bits\n",
                N, reduction_power, threshold, lost_bits)
        end
        println()
    end
end

num_trials = 1000
max_reduction_power = 20

use_clean_multifloat_arithmetic(9)
println("CLEAN (EXTRA BIT)")
run_test(num_trials, max_reduction_power, Val{false}())
println("CLEAN (NO EXTRA)")
run_test(num_trials, max_reduction_power, Val{true}())

use_standard_multifloat_arithmetic(9)
println("STANDARD (EXTRA BIT)")
run_test(num_trials, max_reduction_power, Val{false}())
println("STANDARD (NO EXTRA)")
run_test(num_trials, max_reduction_power, Val{true}())

use_sloppy_multifloat_arithmetic(9)
println("SLOPPY (EXTRA BIT)")
run_test(num_trials, max_reduction_power, Val{false}())
println("SLOPPY (NO EXTRA)")
run_test(num_trials, max_reduction_power, Val{true}())


text:
function compute_matrix(L::LinearMap)
    M = L.M
    N = L.N
    cols = [L*[k==j ? 1.0 : 0.0 for k=1:N] for j=1:N]
    return cat(cols...; dims=2)
end


text:
function test_integration(f1, f2, angles)
    reconstructed = @. f2(f1(angles))
    @test reconstructed ≈ angles
end

@testset "integration - $f1/$f2" for (f1, f2) in zip(
    [deg2hms, deg2dms, deg2ha], [hms2deg, dms2deg, ha2deg]
)
    angles = randdegree(rng, 1000)
    test_integration(f1, f2, angles)
end

@testset "integration - $f1/$f2" for (f1, f2) in zip(
    [rad2hms, rad2dms, rad2ha], [hms2rad, dms2rad, ha2rad]
)
    angles = randrad(rng, 1000)
    test_integration(f1, f2, angles)
end

@testset "integration - $f1/$f2" for (f1, f2) in zip(
    [ha2hms, ha2dms], [hms2ha, dms2ha]
)
    angles = randha(rng, 1000)
    test_integration(f1, f2, angles)
end


text:
# A piece of Data is sent from a lens
immutable Data
  procid::Int # The process that executed the lens
  lensname::Symbol
  data::Dict{Symbol, Any} # Variable ids, hmm
end

getindex(x::Data,i::Symbol) = x.data[i]

# A listener is a named function
# It has a name so that we can remove it if desired easily
type Listener
  name::Symbol
  f::Function
  enabled::Bool
  kwargs::Bool #Does the function take its argument as keywords
end

## Lenses
## ======
# global mutable directory connecting lenses to listeners
lens_to_listeners = Dict{Symbol, Set{Listener}}()

# Two listeners are considered equal if they have the same name
hash(a::Listener) = hash(a.name)
# ==(a::Listener,b::Listener) = a.name == b.name
isequal(a::Listener, b::Listener) = a.name == b.name

## Lens KW
function lens(lensname::Symbol; data...)
#   println("in lens",myid(),data, lens_to_listeners)
  listeners = getlisteners(lensname)
  if !isempty(listeners)
    for listener in listeners
      datum = Data(myid(), lensname, Dict{Symbol, Any}(data))
      listener.enabled && listener.kwargs && listener.f(datum)
      # When listener does not use kwargs just pass in the values
      listener.enabled && !listener.kwargs && listener.f([d[2] for d in data]...)
    end
  end
  v = [d[2] for d in data]
  length(v) == 1 ? v[1] : v
end

lens(lensname::Symbol, data...) =
  lens(lensname; [(symbol("x$i"),data[i]) for i = 1:length(data)]...)

## Enabling/Disabling Listeners
## ==========================
# Register a lens to a listener.
function register!(lensname::Symbol, f::Listener)
#   println("I'm her eand registering!!!\n")
  # DefaultDict: add listener to list or create singleton listener vec
  if haskey(Lens.lens_to_listeners, lensname)
    push!(Lens.lens_to_listeners[lensname],f)
  else
    Lens.lens_to_listeners[lensname] = Set([f])
  end
end

# Create/register a Listener with function and f and name listenername to lensname
register!(f::Function, listenername::Symbol, lensname::Symbol, kwargs::Bool) =
  register!(lensname, Listener(listenername, f, true, kwargs))

# Clearing is permanent
clear_all_listeners!() = global lens_to_listeners = Dict{Symbol, Set{Listener}}()

function delete_listener!(listener_name::Symbol)
  for l in values(lens_to_listeners), f in l
    if f.name == listener_name delete!(l,f) end
  end
end

function enable_listener!(watch_name::Symbol, listener_name::Symbol)
  for f in lens_to_listeners[watch_name]
    if f.name == listener_name f.enabled = true end
  end
end

function disable_listener!(lensname::Symbol, listener_name::Symbol)
  for f in lens_to_listeners[lensname]
    if f.name == listener_name f.enabled = false end
  end
end

enable_all_listeners!() = for fset in values(lens_to_listeners), f in fset f.enabled = true end
disable_all_listeners!() = for fset in values(lens_to_listeners), f in fset f.enabled = false end
# What listeners are registed to lens name
getlisteners(name::Symbol) = haskey(lens_to_listeners, name) ? lens_to_listeners[name] : Set{Listener}()
nlisteners() = (i=0; for fset in values(lens_to_listeners), f in fset i += 1 end; i)


text:
export clenshaw_curtis

"""
    clenshaw_curtis(N)

Compute the nodes `x` and weights `w` for integrating a
continuous functions from [-1,1] using the Clenshaw-Curtis
integration rule with order `N`.
"""
function clenshaw_curtis(N::Int64)

    @assert N>0 "Order must be greater than 0"

    x = zeros(N)
    w = ones(N)

    if N==1
        #x[1]=0.0
        w[1] = 2.0
        return x, w
    end

    # Compute x values
    @inbounds for i=1:N
        x[i] = cos((N-i)*π/(N-1))
    end

    n = N
    x[1] = -1.0

    if mod(n,2)==1
        x[ceil(Int64, (n+1)/2)] = 0.0
    end

    x[n] = 1.0

    @inbounds for i=1:n
        θ = (i-1)*π/(n-1)
        for j=1:ceil(Int64,(n-1)/2)
            if 2*j == n-1
                b = 1.0
            else
                b = 2.0
            end

            w[i] -= b*cos(2*j*θ)/(4*j*j-1)
        end
    end

    w[1] *= 1/(n-1)
    w[2:n] .*= 2.0/(n-1)
    w[n] *= 1/(n-1)

    return x, w
end


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

export
    OBJC_NEW_PROPERTIES,
    AVFILTER_FLAG_DYNAMIC_INPUTS,
    AVFILTER_FLAG_DYNAMIC_OUTPUTS,
    AVFILTER_FLAG_SLICE_THREADS,
    AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC,
    AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL,
    AVFILTER_FLAG_SUPPORT_TIMELINE,
    AVFILTER_THREAD_SLICE,
    AVFILTER_CMD_FLAG_ONE,
    AVFILTER_CMD_FLAG_FAST,
    AVFilterPad,
    AVFilter,
    AVFilterGraphInternal,
    AVFilterGraph,
    AVFilterInternal,
    AVFilterContext,
    AVFilterFormats,
    ANONYMOUS_1,
    AVFILTER_AUTO_CONVERT_ALL,
    AVFILTER_AUTO_CONVERT_NONE,
    AVFilterInOut,
    LIBAVFILTER_VERSION_MAJOR,
    LIBAVFILTER_VERSION_MINOR,
    LIBAVFILTER_VERSION_MICRO,
    LIBAVFILTER_BUILD,
    AV_BUFFERSINK_FLAG_PEEK,
    AV_BUFFERSINK_FLAG_NO_REQUEST,
    AVBufferSinkParams,
    AVABufferSinkParams,
    ANONYMOUS_2,
    AV_BUFFERSRC_FLAG_NO_CHECK_FORMAT,
    AV_BUFFERSRC_FLAG_PUSH,
    AV_BUFFERSRC_FLAG_KEEP_REF

const OBJC_NEW_PROPERTIES  =  1
const AVFILTER_FLAG_DYNAMIC_INPUTS  =  1 << 0
const AVFILTER_FLAG_DYNAMIC_OUTPUTS  =  1 << 1
const AVFILTER_FLAG_SLICE_THREADS  =  1 << 2
const AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC  =  1 << 16
const AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL  =  1 << 17
const AVFILTER_FLAG_SUPPORT_TIMELINE  =  AVFILTER_FLAG_SUPPORT_TIMELINE_GENERIC | AVFILTER_FLAG_SUPPORT_TIMELINE_INTERNAL
const AVFILTER_THREAD_SLICE  =  1 << 0
const AVFILTER_CMD_FLAG_ONE  =  1
const AVFILTER_CMD_FLAG_FAST  =  2

const AVFilterPad = Void

immutable AVFilter
    name::Cstring
    description::Cstring
    inputs::Ptr{AVFilterPad}
    outputs::Ptr{AVFilterPad}
    priv_class::Ptr{AVClass}
    flags::Cint
    init::Ptr{Void}
    init_dict::Ptr{Void}
    uninit::Ptr{Void}
    query_formats::Ptr{Void}
    priv_size::Cint
    next::Ptr{AVFilter}
    process_command::Ptr{Void}
    init_opaque::Ptr{Void}
end

const AVFilterGraphInternal = Void

immutable AVFilterGraph
    av_class::Ptr{AVClass}
    filters::Ptr{Ptr{AVFilterContext}}
    nb_filters::UInt32
    scale_sws_opts::Cstring
    resample_lavr_opts::Cstring
    thread_type::Cint
    nb_threads::Cint
    internal::Ptr{AVFilterGraphInternal}
    opaque::Ptr{Void}
    execute::Ptr{avfilter_execute_func}
    aresample_swr_opts::Cstring
    sink_links::Ptr{Ptr{AVFilterLink}}
    sink_links_count::Cint
    disable_auto_convert::UInt32
end

const AVFilterInternal = Void

immutable AVFilterContext
    av_class::Ptr{AVClass}
    filter::Ptr{AVFilter}
    name::Cstring
    input_pads::Ptr{AVFilterPad}
    inputs::Ptr{Ptr{AVFilterLink}}
    nb_inputs::UInt32
    output_pads::Ptr{AVFilterPad}
    outputs::Ptr{Ptr{AVFilterLink}}
    nb_outputs::UInt32
    priv::Ptr{Void}
    graph::Ptr{AVFilterGraph}
    thread_type::Cint
    internal::Ptr{AVFilterInternal}
    command_queue::Ptr{AVFilterCommand}
    enable_str::Cstring
    enable::Ptr{Void}
    var_values::Ptr{Cdouble}
    is_disabled::Cint
end

const AVFilterFormats = Void

# begin enum ANONYMOUS_1
const ANONYMOUS_1 = Cint
const AVFILTER_AUTO_CONVERT_ALL  =  (Int32)(0)
const AVFILTER_AUTO_CONVERT_NONE  =  (Int32)(-1)
# end enum ANONYMOUS_1

immutable AVFilterInOut
    name::Cstring
    filter_ctx::Ptr{AVFilterContext}
    pad_idx::Cint
    next::Ptr{AVFilterInOut}
end

const LIBAVFILTER_VERSION_MAJOR  =  6
const LIBAVFILTER_VERSION_MINOR  =  31
const LIBAVFILTER_VERSION_MICRO  =  100

# Skipping MacroDefinition: LIBAVFILTER_VERSION_INT AV_VERSION_INT ( LIBAVFILTER_VERSION_MAJOR , LIBAVFILTER_VERSION_MINOR , LIBAVFILTER_VERSION_MICRO )
# Skipping MacroDefinition: LIBAVFILTER_VERSION AV_VERSION ( LIBAVFILTER_VERSION_MAJOR , LIBAVFILTER_VERSION_MINOR , LIBAVFILTER_VERSION_MICRO )

const LIBAVFILTER_BUILD  =  LIBAVFILTER_VERSION_INT

# Skipping MacroDefinition: LIBAVFILTER_IDENT "Lavfi" AV_STRINGIFY ( LIBAVFILTER_VERSION )
# Skipping MacroDefinition: FF_API_OLD_FILTER_OPTS ( LIBAVFILTER_VERSION_MAJOR < 7 )
# Skipping MacroDefinition: FF_API_OLD_FILTER_OPTS_ERROR ( LIBAVFILTER_VERSION_MAJOR < 7 )
# Skipping MacroDefinition: FF_API_AVFILTER_OPEN ( LIBAVFILTER_VERSION_MAJOR < 7 )
# Skipping MacroDefinition: FF_API_AVFILTER_INIT_FILTER ( LIBAVFILTER_VERSION_MAJOR < 7 )
# Skipping MacroDefinition: FF_API_OLD_FILTER_REGISTER ( LIBAVFILTER_VERSION_MAJOR < 7 )
# Skipping MacroDefinition: FF_API_NOCONST_GET_NAME ( LIBAVFILTER_VERSION_MAJOR < 7 )

const AV_BUFFERSINK_FLAG_PEEK  =  1
const AV_BUFFERSINK_FLAG_NO_REQUEST  =  2

immutable AVBufferSinkParams
    pixel_fmts::Ptr{AVPixelFormat}
end

immutable AVABufferSinkParams
    sample_fmts::Ptr{AVSampleFormat}
    channel_layouts::Ptr{Int64}
    channel_counts::Ptr{Cint}
    all_channel_counts::Cint
    sample_rates::Ptr{Cint}
end

# begin enum ANONYMOUS_2
const ANONYMOUS_2 = UInt32
const AV_BUFFERSRC_FLAG_NO_CHECK_FORMAT  =  (UInt32)(1)
const AV_BUFFERSRC_FLAG_PUSH  =  (UInt32)(4)
const AV_BUFFERSRC_FLAG_KEEP_REF  =  (UInt32)(8)
# end enum ANONYMOUS_2

max_stars_repo_name: asafasban/AAVideoIO

max_issues_repo_name: asafasban/AAVideoIO

max_forks_repo_name: asafasban/AAVideoIO


text:
function holtslag_psis(xi)
    #println("xi")
    ##println(xi)
    #println("-------------")
if (xi>=0) & (xi<=1.)
    y=1+5*xi;
elseif xi>=1
    y=5+xi
else xi<0
    y=(1-15*xi)^(-1/2)
end
return y
end


text:
module spring_pendulum
using ReachabilityAnalysis
@taylorize function spring_pendulum!(du, u, p, t)
    local g, k, L = 9.8, 2.0, 1.0
    du[1] = u[3]
    du[2] = u[4]
    du[3] = ((u[1] * (u[4]*u[4])) + g*cos(u[2])) - k*(u[1]-L)
    du[4] = -((2*u[3]*u[4]) + g*sin(u[2]))/(u[1])

    ## change of variables: r <- r - L
    #du[3] = (((u[1]+L) * (u[4]*u[4])) + g*cos(u[2])) - k*(u[1])
    #du[4] = -((2*u[3]*u[4]) + g*sin(u[2]))/(u[1]+L)

    return du
end

function model(X0)
    S = @system(x' = spring_pendulum!(x), dim:4)
    return IVP(S, X0)
end
end # module

max_stars_repo_name: mforets/ReachabilityLibrary

max_issues_repo_name: mforets/ReachabilityLibrary

max_forks_repo_name: mforets/ReachabilityLibrary


text:
using ClusterManagers
using Distributed

OnCluster = false #set to false to run locally
addWorkers = false #set to false to run serially
println("OnCluster = $(OnCluster)")

# Current number of workers
#--------------------------
currentWorkers = nworkers()
println("Initial number of workers = $(currentWorkers)")

# Increase the number of workers available
#-----------------------------------------
maxNumberWorkers = 1
if addWorkers == true
    if OnCluster == true
        addprocs(SlurmManager(maxNumberWorkers))
    else
        addprocs(maxNumberWorkers)
    end
end

# Sanity checks
#-------------
hosts = []
pids = []
for i in workers()
    host, pid = fetch(@spawnat i (gethostname(), getpid()))
    println("Hello I am worker $(i), my host is $(host)")
    push!(hosts, host)
    push!(pids, pid)
end

currentWorkers = nworkers()
println("Number of workers = $(currentWorkers)")

using Plots
using LaTeXStrings
using ParallelDataTransfer

using MSM
using DataStructures
using OrderedCollections
using Distributions
using Random
using DataStructures
using Statistics
using LinearAlgebra

# Generate simulated data
Random.seed!(1234)  #for replicability reasons
T = 10000         #number of periods
P = 2               #number of dependent variables
beta0 = rand(P)     #choose true coefficients by drawing from a uniform distribution on [0,1]
alpha0 = rand(1)[]  #intercept
theta0 = 0.0        #coefficient to create serial correlation in the error terms
println("True intercept = $(alpha0)")
println("True coefficient beta0 = $(beta0)")
println("Serial correlation coefficient theta0 = $(theta0)")

# Generation of error terms
# row = individual dimension
# column = time dimension
U = zeros(T)
d = Normal()
U[1] = rand(d, 1)[] #first error term
# loop over time periods
for t = 2:T
    U[t] = rand(d, 1)[] + theta0 * U[t-1]
end

# Let's simulate the dependent variables x_t
x = zeros(T, P)

d = Uniform(0, 5)
for p = 1:P
    x[:, p] = rand(d, T)
end

# Let's calculate the resulting y_t
y = zeros(T)

for t = 1:T
    y[t] = alpha0 + x[t, 1] * beta0[1] + x[t, 2] * beta0[2] + U[t]
end

# Visualize data
p1 = scatter(
    x[1:100, 1],
    y[1:100],
    xlabel = "x1",
    ylabel = "y",
    legend = :none,
    smooth = true,
)
p2 = scatter(
    x[1:100, 2],
    y[1:100],
    xlabel = "x2",
    ylabel = "y",
    legend = :none,
    smooth = true,
)
plot(p1, p2)

# Define locally
optionsSMM = MSMOptions(
    maxFuncEvals = 2000,
    globalOptimizer = :dxnes,
    localOptimizer = :NelderMead,
)
myProblem = MSMProblem(options = optionsSMM);

# Priors
dictPriors = OrderedDict{String,Array{Float64,1}}()
dictPriors["alpha"] = [0.5, 0.001, 1.0]
dictPriors["beta1"] = [0.5, 0.001, 1.0]
dictPriors["beta2"] = [0.5, 0.001, 1.0]

# Empirical moments
dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
dictEmpiricalMoments["mean"] = [mean(y)] #informative on the intercept
dictEmpiricalMoments["mean_x1y"] = [mean(x[:, 1] .* y)] #informative on betas
dictEmpiricalMoments["mean_x2y"] = [mean(x[:, 2] .* y)] #informative on betas
dictEmpiricalMoments["mean_x1y^2"] = [mean((x[:, 1] .* y) .^ 2)] #informative on betas
dictEmpiricalMoments["mean_x2y^2"] = [mean((x[:, 2] .* y) .^ 2)] #informative on betas

W = Matrix(1.0 .* I(length(dictEmpiricalMoments)))#initialization
#Special case: diagonal matrix
#(you may choose something else)
for (indexMoment, k) in enumerate(keys(dictEmpiricalMoments))
    W[indexMoment, indexMoment] = 1.0 / (dictEmpiricalMoments[k][1])^2
end

set_priors!(myProblem, dictPriors)
set_empirical_moments!(myProblem, dictEmpiricalMoments)
set_weight_matrix!(myProblem, W)

# x[1] corresponds to the intercept, x[2] corresponds to beta1, x[3] corresponds to beta2
function functionLinearModel(
    x;
    uniform_draws::Array{Float64,1},
    simX::Array{Float64,2},
    nbDraws::Int64 = length(uniform_draws),
    burnInPerc::Int64 = 0,
)
    T = nbDraws
    P = 2       #number of dependent variables

    alpha = x[1]
    beta = x[2:end]
    theta = 0.0     #coefficient to create serial correlation in the error terms

    # Creation of error terms
    # row = individual dimension
    # column = time dimension
    U = zeros(T)
    d = Normal()
    # Inverse cdf (i.e. quantile)
    gaussian_draws = quantile.(d, uniform_draws)
    U[1] = gaussian_draws[1] #first error term

    # loop over time periods
    for t = 2:T
        U[t] = gaussian_draws[t] + theta * U[t-1]
    end

    # Let's calculate the resulting y_t
    y = zeros(T)

    for t = 1:T
        y[t] = alpha + simX[t, 1] * beta[1] + simX[t, 2] * beta[2] + U[t]
    end

    # Get rid of the burn-in phase:
    #------------------------------
    startT = max(1, Int(nbDraws * (burnInPerc / 100)))

    # Moments:
    #---------
    output = OrderedDict{String,Float64}()
    output["mean"] = mean(y[startT:nbDraws])
    output["mean_x1y"] = mean(simX[startT:nbDraws, 1] .* y[startT:nbDraws])
    output["mean_x2y"] = mean(simX[startT:nbDraws, 2] .* y[startT:nbDraws])
    output["mean_x1y^2"] = mean((simX[startT:nbDraws, 1] .* y[startT:nbDraws]) .^ 2)
    output["mean_x2y^2"] = mean((simX[startT:nbDraws, 2] .* y[startT:nbDraws]) .^ 2)
    return output
end

# Let's freeze the randomness during the minimization
d_Uni = Uniform(0, 1)
nbDraws = T #number of draws in the simulated data
uniform_draws = rand(d_Uni, nbDraws)
simX = zeros(length(uniform_draws), 2)
d = Uniform(0, 5)
for p = 1:2
    simX[:, p] = rand(d, length(uniform_draws))
end

# Construct the objective function everywhere
set_simulate_empirical_moments!(
    myProblem,
    x -> functionLinearModel(x, uniform_draws = uniform_draws, simX = simX),
)
construct_objective_function!(myProblem)

println("Global Algorithm")
# Choose a global optimizer that supports parallel evaluations (e.g. xnes or dxnes)
# (see the documentation: https://github.com/robertfeldt/BlackBoxOptim.jl)
msm_optimize!(myProblem, verbose = false)

# Access the results using best_candidate
minimizer = msm_minimizer(myProblem)
minimum_val = msm_minimum(myProblem)

println("Minimum objective function = $(minimum_val)")
println("Estimated value for alpha = $(minimizer[1]). True value for beta1 = $(alpha0[1]) \n")
println("Estimated value for beta1 = $(minimizer[2]). True value for beta1 = $(beta0[1]) \n")
println("Estimated value for beta2 = $(minimizer[3]). True value for beta2 = $(beta0[2]) \n")

println("Multistart Algorithm")
# Start several local optimization algorithms in parallel
# Choose algorithms from the package Optim.jl (https://github.com/JuliaNLSolvers/Optim.jl)
# The "global" mininimum is the minimum of the local minima.
msm_multistart!(myProblem, nums = nworkers(), verbose = false)
minimizer_multistart = msm_multistart_minimizer(myProblem)
minimum_multistart = msm_multistart_minimum(myProblem)

println("Minimum objective function = $(minimum_multistart)")
println("Estimated value for alpha = $(minimizer_multistart[1]). True value for beta1 = $(alpha0[1]) \n")
println("Estimated value for beta1 = $(minimizer_multistart[2]). True value for beta1 = $(beta0[1]) \n")
println("Estimated value for beta2 = $(minimizer_multistart[3]). True value for beta2 = $(beta0[2]) \n")

# Empirical Series
#-----------------
X = zeros(T, length(dictEmpiricalMoments))
X[:, 1] = y
X[:, 2] = (x[:, 1] .* y)
X[:, 3] = (x[:, 2] .* y)
X[:, 4] = (x[:, 1] .* y) .^ 2
X[:, 5] = (x[:, 2] .* y) .^ 2

# "Distance Matrix" (see Duffie and Singleton, 1993)
Sigma0 = cov(X)

# Heatmap to visualize correlation
xs = [string("x", i) for i = 1:length(dictEmpiricalMoments)]
ys = [string("x", i) for i = 1:length(dictEmpiricalMoments)]
z = cor(X)
hh = heatmap(xs, ys, z, aspect_ratio = 1)

set_Sigma0!(myProblem, Sigma0)
# nbDraws = number of draws in the simulated data
# To decrease standard errors, increase nbDraws
calculate_Avar!(myProblem, minimizer_multistart, tau = T / nbDraws)

df = summary_table(myProblem, minimizer_multistart, T, 0.05)
println(df)

# Check the rank condition
# Local identification requires D to be full column rank (in a neighborhood of the solution)
D = calculate_D(myProblem, minimizer_multistart)
println("number of parameters: $(size(D,2))")
println("rank of D is: $(rank(D))")

# Slices
vXGrid, vYGrid = msm_slices(myProblem, minimizer_multistart, nbPoints = 7);
p1 = plot(vXGrid[:, 1],vYGrid[:, 1],title = L"\alpha", label = "",linewidth = 3, xrotation = 45)
plot!(p1, [minimizer_multistart[1]], seriestype = :vline, label = "",linewidth = 1)
p2 = plot(vXGrid[:, 2],vYGrid[:, 2],title = L"\beta_1", label = "",linewidth = 3, xrotation = 45)
plot!(p2, [minimizer_multistart[2]], seriestype = :vline, label = "",linewidth = 1)
p3 = plot(vXGrid[:, 3],vYGrid[:, 3],title = L"\beta_2", label = "",linewidth = 3, xrotation = 45)
plot!(p3, [minimizer_multistart[3]], seriestype = :vline, label = "",linewidth = 1)
plot_combined = plot(p1, p2, p3)
savefig(plot_combined, "slices.png")
display(plot_combined)

# Efficient GMM
#---------------
myProblem.W = inv(Sigma0)
msm_multistart!(myProblem, nums = nworkers(), verbose = false)
minimizer_multistart = msm_multistart_minimizer(myProblem)
minimum_multistart = msm_multistart_minimum(myProblem)

# J-test
J, c = J_test(myProblem, minimizer_multistart, T, nbDraws, 0.05)
println("J-statistic: $(J)")
println("Critical value: $(c)")

# Compare results with GLM
using DataFrames, GLM
data = DataFrame(x1 = x[:, 1], x2 = x[:, 2], y = y[:]);
ols = lm(@formula(y ~ x1 + x2), data)
println(ols)


text:
"""
    module GenerateMediumDict

Module that is used to generate the MediumDict and store it on file.

This package is currently under development.
"""
module GenerateMediumDict

### Importing packages -------------------------------------------------------------------------------
import ModiaMedia
using  JSON
using  StaticArrays
import Serialization

const path    = dirname(dirname(@__FILE__))          # Absolute path of package directory
const Version = "0.1.0-dev from 2018-11-13 21:34"
const dict    = Dict{AbstractString,ModiaMedia.AbstractMedium}()

println(" \nImporting GenerateMediumDict version ", Version)

### Utility functions
"""
    filled_obj = fillobj(dict,obj)

Store the values of all (string) keys of dict in struct obj and return the filled obj.
"""
function fillobj(dict, obj)
    for (key, value) in dict
        # println("... value = ", value, ", typeof(value) = ", typeof(value))
        if typeof(value) <: Array{Any,1}
           setfield!(obj, Symbol(key), SVector{length(value),Float64}(value))
        else
           setfield!(obj, Symbol(key), value)
        end
    end
    return obj
end

### Including data --------------------------------------------------------
include("SimpleMedium.jl")
include("SimpleIdealGasMedium.jl")
include("SingleGasesNasa.jl")

### Write data to file ----------------------------------------------------
file   = "$path/src/Media/media.julia_serializer"
nmedia = length(dict)
println("... Write media dictionary ($nmedia media) as serialized binary object to file:\n",
        "    ", file)
open("$path/src/Media/media.julia_serializer", "w") do f
    Serialization.serialize(f, dict)
end

end # module


text:
"""
    arrow(p::Point, f::F; kwargs...) where F <: Force

Draw a single force arrow with filled head.

    p        starting position.
    f        scalar force quantity.
    α = 0°   is direction to the right, positive values rotate around z axis.

"""
arrow(p::Point, f::F; kwargs...) where F <: Force = arrow(p, f, 0 * f; kwargs...)

"""
    arrow(p::Point, f::ForceTuple; kwargs...)

Draw a single velocity arrow - no filled arrow head.

p     starting position (at the normal between vanes tip and arrow spine)
F     is a force tuple referred to an x-axis which can optionally be rotated around z by
      keyword argument α.
components = false: No components forces on the x and y axes are drawn.
"""
function arrow(p::Point, f::ForceTuple; components = true, kwargs...)
    if !components
        αcalc = atan(f[2], f[1]) |> rad |> °
        F = hypot(f)
        arrow(p, F, 0.0 * f[2]; kwargs..., α = αcalc)
    else
        arrow(p, f[1], f[2]; kwargs...)
    end
end

"""
    arrow(p::Point, Fx::F, Fy::F;
        α::Angle = 0°,
        linewidth = PT,
        backgroundcolor = colorant"white",
        labellength::Bool = false) where {F <: Force}

Draw a force arrow with filled head resultant, and perpendicular components in dashed style.
p                  starting position.
Fx, Fy             scalar velocity perpendicular components.
α                  defines the direction of the x component.
backgroundcolor    is used to pick colors for components and label outlines
labellength = true adds a label with arrow length (i.e. vector norm) and unit
"""
function arrow(p::Point, Fx::F, Fy::F;
                α::Angle = 0°,
                linewidth = PT,
                backgroundcolor = colorant"white",
                labellength::Bool = false) where {F <: Force}
    ΔFx = Point(Fx * cos(α), Fx *sin(α))
    ΔFy = Point(-Fy * sin(α), Fy *cos(α))
    ΔF = ΔFx + ΔFy

    # Prepare color selection
    luminback = lumin(backgroundcolor)
    luminfront = get_current_lumin()
    avglumin = (luminback + luminfront) / 2
    curcol = get_current_RGB()
    avgcol = color_with_lumin(curcol, avglumin)
    # Draw the components, but only if the resultant is not coincident with one of them.
    if hypot(ΔFx) > 0 && hypot(ΔFy) > 0
        gsave()
        setdash("longdashed")
        deltalumin = luminback - luminfront
        lesscontrast = color_with_lumin(get_current_RGB(), luminfront + deltalumin)
        sethue(lesscontrast)
        shaftlength = hypot(ΔFx)
        arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
        arrow_nofill(p, p + ΔFx, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT / 2)
        shaftlength = hypot(ΔFy)
        arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
        arrow_nofill(p, p + ΔFy, arrowheadlength = arrowheadlength, arrowheadangle = π / 12, linewidth = PT / 2)
        grestore()
    end

    # Draw the resultant arrow
    endpoint = p + ΔF
    shaftlength = hypot(ΔF)
    arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3
    if !isapprox(shaftlength, 0.0)
        arrow(p, endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
    else
        # with shaft length approximately zero, draw two opposing arrows with minimum possible shaft lengths
        tiny = oneunit(shaftlength) # almost zero
        arrow(p + ( tiny, 0), endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
        arrow(p + (-tiny, 0), endpoint, arrowheadlength = arrowheadlength, arrowheadangle = pi/12, linewidth = PT)
    end

    # add an optional label with the euclidean length of the resultant, at the tip of it.
    if labellength
        sethue(avgcol)
        fontface("Calibri bold")
        setlabel(endpoint, Fx, Fy, α)
        fontface("Calibri")
        sethue(curcol)
        setlabel(endpoint, Fx, Fy, α)
    end
end

"""
arrow(p::Point, v::V; kwargs...) where V<:Velocity

Draw a single velocity arrow - no filled arrow head.

p        starting position (at the normal between vanes tip and arrow spine).
v        scalar velocity
α = 0°   is direction to the right, positive values rotate around z axis.
"""
arrow(p::Point, v::V; kwargs...) where V<:Velocity = arrow(p, v, 0 * v; kwargs...)

"""
    arrow(p::Point, v::VelocityTuple; kwargs...)

Draw a single velocity arrow - no filled arrow head.

p     starting position (at the normal between vanes tip and arrow spine)
v     is a 2-tuple here, but you can optionally rotate the x-axis by keyword argument α.
"""
arrow(p::Point, v::VelocityTuple; kwargs...) = arrow(p, v[1], v[2]; kwargs...)

"""
    arrow(p::Point, vx::V, vy::V;
               α::Angle = 0°,
               linewidth = PT,
               backgroundcolor = colorant"white",
               labellength::Bool = false) where {V <: Velocity}

Draw a single velocity arrow - no filled arrow head.

p                   starting position (at the normal between vanes tip and arrow spine).
vx, vy              scalar velocity perpendicular components.
α                   defines the direction of the x component.
backgroundcolor     is used to pick a color for an outline.
labellength = true  label at the point with arrow length and unit
"""
function arrow(p::Point, vx::V, vy::V;
               α::Angle = 0°,
               linewidth = PT,
               backgroundcolor = colorant"white",
               labellength::Bool = false) where {V <: Velocity}
    Δvx = Point(vx * cos(α), vx *sin(α))
    Δvy = Point(-vy * sin(α), vy *cos(α))
    Δv = Δvx + Δvy
    endpoint = p + Δv
    luminback = lumin(backgroundcolor)
    luminfront = get_current_lumin()
    avglumin = (luminback + luminfront) / 2
    curcol = get_current_RGB()
    avgcol = color_with_lumin(backgroundcolor, avglumin)
    gsave()
    setline(linewidth * 2)
    sethue(avgcol)
    curved_point_arrow_with_vanes(p, Δv, endpoint, linewidth)
    sethue(curcol)
    setline(linewidth)
    curved_point_arrow_with_vanes(p, Δv, endpoint, linewidth)
    # add a label with the euclidean length of the arrow, at the tip of it.
    if labellength
        sethue(avgcol)
        fontface("Calibri bold")
        setlabel(endpoint, vx, vy, α)
        fontface("Calibri")
        sethue(curcol)
        setlabel(endpoint, vx, vy, α)
    end
    grestore()
    return nothing
end

function curved_point_arrow_with_vanes(p, Δv, endpoint, linewidth)
    arrowheadangle = pi/24

    isapprox(Δv, Point(0,0)) && throw(error("can't draw velocity arrow between two identical points"))
    shaftlength = hypot(Δv)
    arrowheadlength = shaftlength > 3 * EM ? EM :  shaftlength / 3

    shaftangle = atan(p.y - endpoint.y, p.x - endpoint.x)

    # shorten the length so that lines
    # stop before we get to the arrow
    # thus wide shafts won't stick out through the head of the arrow.
    max_undershoot = shaftlength - ((linewidth/2) / tan(arrowheadangle))
    ratio = max_undershoot / shaftlength
    tox = p.x + Δv.x * ratio
    toy = p.y + Δv.y * ratio
    fromx = p.x
    fromy = p.y

    # draw the shaft of the arrow
    newpath()
    line(Point(fromx, fromy), Point(tox, toy), :stroke)

    # draw the arrowhead
    arrowheadtopsideangle = shaftangle + arrowheadangle
    arrowheadtopsidemidangle = shaftangle + 0.67arrowheadangle
    topx = endpoint.x + cos(arrowheadtopsideangle) * arrowheadlength
    topy = endpoint.y + sin(arrowheadtopsideangle) * arrowheadlength
    topmidx = endpoint.x + cos(arrowheadtopsidemidangle) * arrowheadlength / 2
    topmidy = endpoint.y + sin(arrowheadtopsidemidangle) * arrowheadlength / 2

    arrowheadbottomsideangle = shaftangle - arrowheadangle
    arrowheadbottomsidemidangle = shaftangle - 0.67arrowheadangle
    botx = endpoint.x + cos(arrowheadbottomsideangle) * arrowheadlength
    boty = endpoint.y + sin(arrowheadbottomsideangle) * arrowheadlength
    botmidx = endpoint.x + cos(arrowheadbottomsidemidangle) * arrowheadlength / 2
    botmidy = endpoint.y + sin(arrowheadbottomsidemidangle) * arrowheadlength / 2

    poly([Point(topx, topy), Point(topmidx, topmidy), endpoint,
        Point(botmidx, botmidy), Point(botx, boty)], :stroke)

    # draw the rearmost vane
    vaneangle = shaftangle + π / 3
    vanelength = arrowheadlength / 2
    vanex = fromx + cos(vaneangle) * vanelength
    vaney = fromy + sin(vaneangle) * vanelength

    # draw the mid vane
    from1x = p.x + Δv.x * (1 - ratio)
    from1y = p.y + Δv.y * (1 - ratio)
    vane1x = from1x + cos(vaneangle) * vanelength
    vane1y = from1y + sin(vaneangle) * vanelength

    # draw the front vane
    from2x = p.x + Δv.x * 2 * (1 - ratio)
    from2y = p.y + Δv.y * 2 * (1 - ratio)
    vane2x = from2x + cos(vaneangle) * vanelength
    vane2y = from2y + sin(vaneangle) * vanelength

    #poly([Point(vanex, vaney), Point(fromx, fromy)] , :stroke)
    #poly([Point(vane1x, vane1y), Point(from1x, from1y)] , :stroke)
    #poly([Point(vane2x, vane2y), Point(from2x, from2y)] , :stroke)

    poly([Point(vanex, vaney), Point(fromx, fromy), Point(from1x, from1y),
        Point(vane1x, vane1y), Point(from1x, from1y), Point(from2x, from2y),
        Point(vane2x, vane2y)] , :stroke)

end

"""
Same as Luxor.arrow, except for leaving the arrow head unfilled
"""
function arrow_nofill(startpoint::Point, endpoint::Point;
    linewidth=1.0,
    arrowheadlength= EM / 4.4,
    arrowheadangle = π / 8,
    decoration = 0.5,
    decorate = () -> ())
    gsave()
    setlinejoin("butt")
    setline(linewidth)

    isapprox(startpoint, endpoint) && throw(error("can't draw arrow between two identical points"))
    shaftlength = distance(startpoint, endpoint)
    shaftangle = atan(startpoint.y - endpoint.y, startpoint.x - endpoint.x)

    arrowheadtopsideangle = shaftangle + arrowheadangle

    # shorten the length so that lines
    # stop before we get to the arrow
    # thus wide shafts won't stick out through the head of the arrow.
    max_undershoot = shaftlength - ((linewidth/2) / tan(arrowheadangle))
    ratio = max_undershoot/shaftlength
    tox = startpoint.x + (endpoint.x - startpoint.x) * ratio
    toy = startpoint.y + (endpoint.y - startpoint.y) * ratio
    fromx = startpoint.x
    fromy = startpoint.y

    # draw the shaft of the arrow
    newpath()
    line(Point(fromx, fromy), Point(tox, toy), :stroke)

    # prepare to add decorations at point along shaft
    for decpos in decoration
        decpoint = between(startpoint, endpoint, decpos)
        slp = slope(startpoint, endpoint)
        @layer begin
            translate(decpoint)
            rotate(slp)
            decorate()
        end
    end
    # draw the arrowhead
    draw_arrowhead(endpoint, shaftangle;
        arrowheadangle,
        arrowheadlength, filled = false)
    grestore()
end
"""
    draw_arrowhead(endpoint, shaftangle;
        arrowheadangle = π / 8,
        arrowheadlength  = EM / 4.4; filled = false)
"""
function draw_arrowhead(endpoint, shaftangle;
    arrowheadangle = π / 8,
    arrowheadlength  = EM / 4.4, filled = false)
    arrowheadtopsideangle = shaftangle + arrowheadangle
    # draw the arrowhead
    topx = endpoint.x + cos(arrowheadtopsideangle) * arrowheadlength
    topy = endpoint.y + sin(arrowheadtopsideangle) * arrowheadlength
    arrowheadbottomsideangle = shaftangle - arrowheadangle
    botx = endpoint.x + cos(arrowheadbottomsideangle) * arrowheadlength
    boty = endpoint.y + sin(arrowheadbottomsideangle) * arrowheadlength
    poly([Point(topx, topy), endpoint, Point(botx, boty)], filled ? :fill : :stroke)
end

"""
    arrow(centerpos::Point, radius, α_sta::Angle, α_end::Angle;
            linewidth=1.0,
            arrowheadlength = EM / 4.4,
            arrowheadangle = π / 8,
            decoration = 0.5,
            decorate = () -> (),
            clockwise = false)

Draw a curved arrow, an arc centered at `centerpos` starting at `α_sta` and
ending at `α_end`.

The arrowhead at the end is optional. Angles are measured counter-clockwise
from the positive x-axis.

Arrows don't use the current linewidth setting (`setline()`); you can specify
the linewidth.

The `decorate` keyword argument accepts a function that can execute code at
locations on the arrow's shaft. The inherited graphic environment is centered at
points on the curve between 0 and 1 given by scalar or vector `decoration`, and
the x-axis is aligned with the direction of the curve at that point.
"""
function arrow(centerpos::Point, radius, α_sta::Angle, α_end::Angle;
        linewidth=1.0,
        arrowheadlength = EM / 4.4,
        arrowheadangle = π / 8,
        decoration = 0.5,
        decorate = () -> (),
        clockwise = false,
        filled = false)
    @assert radius >= zero(radius) " radius < 0 not supported. Flip keyword argument clockwise?"
    gsave()
    setlinejoin("butt")
    setline(linewidth)
    # Radius in points, unitless
    r = scale_to_pt(oneunit(radius)) * radius / unit(radius)
    while α_sta > α_end
        α_end += unit(α_end)(360°)
    end
    # Full arc ccw angle span. Negative values for clockwise pointing arrows
    Δα = clockwise ? α_sta - α_end : α_end - α_sta

    # End arrow arc ccw angle span. Negative values for clockwise pointing arrows.
    Δα_head = unit(α_end)(arrowheadlength / r) * (clockwise ? -1 : 1)
    α_head = α_end - Δα_head

    # Is there room for one arrow head?
    if sign(Δα_head) != sign(Δα) && abs(Δα_head) < abs(Δα)
        @warn "Arrow head too large, try negative keyword argument arrowheadlength"
        return
    end

    # In case we are drawing a sharp arrow head point and a wide shaft, we want to end the arc somewhere inside of the head.
    # Sign convention: Positive values for cutting off ccw arcs.
    Δα_cutoff = unit(α_end)((linewidth / 2) / tan(arrowheadangle) / r) * sign(Δα)
    # Draw the arc to
    α_cutoff = unit(α_end)(α_end - Δα_cutoff)

    # Local origo at centre
    translate(centerpos)
    # Start, end, head crossing points and end of drawn arc in local pt coordinates, where +y is down
    p0 = Point(r∙cos(α_sta), -r∙sin(α_sta))
    p1 = Point(r∙cos(α_end), -r∙sin(α_end))
    p1h = Point(r∙cos(α_end - Δα_head), -r∙sin(α_end - Δα_head))
    p1c = Point(r∙cos(α_end - Δα_cutoff), -r∙sin(α_end - Δα_cutoff))

    # Unitless, clockwise radian angles from origo
    β_sta = rad(-α_sta) / 1rad
    β_head = rad(-α_head) / 1rad
    β_cutoff = rad(-α_cutoff) / 1rad
    β_end = rad(-α_end) / 1rad

    # Draw the cut-off shaft arc
    newpath()
    move(p0.x, p0.y)
    if clockwise
        arc(0, 0, r, β_sta, β_cutoff, :stroke)
    else
        carc(0, 0, r,  β_sta, β_cutoff, :stroke)
    end
    closepath()
    # prepare to add decorations at points along shaft

    for decpos in decoration
        decorationangle = rescale(decpos, 0, 1, β_sta, β_end)
        decorationpoint = Point(r * cos(decorationangle), r * sin(decorationangle))
        perp = perpendicular(decorationpoint)
        @layer begin
            translate(decorationpoint)
            rotate(slope(decorationpoint, perp))
            decorate()
        end
    end

    # Head
    draw_arrowhead(p1, β_cutoff + (clockwise ? -π / 2 : π / 2 );
        arrowheadangle, arrowheadlength, filled)
    grestore()
    return
end


text:
background_human_zero_markov()   = [0.2953, 0.2047, 0.2047, 0.2953]
background_xenopus_zero_markov() = [0.2996, 0.2004, 0.2004, 0.2996]

### Load meme into named tuple
function loadmeme(file, w = 1e-2, background=background_human_zero_markov())

    motifname = ""
    motifid   = ""
    mode = 1
    data = Vector{Vector{Float64}}()
    for line in eachline(file)
        if mode == 1
            if startswith(line, "MOTIF")
                fields = split(line)
                motifid   = fields[2]
                motifname = fields[3]
            elseif startswith(line, "letter-probability")
                mode = 2
            end
        else
            if startswith(line, "URL")
                mode = 1
                break
            end
            v = parse.(Float64, split(line))
            push!(data, v)
        end
    end
    pwm = reduce(hcat, data)
    pwm .+= w
    pwm ./= sum(pwm, dims=1)

    pbg = log2.(pwm) .- log2.(background)
    (name=motifname, id=motifid, pwm=pwm, pbg=pbg)
end

function loadmemelibrary(file, w = 1e-2, background=background_human_zero_markov())

    motifnames = String[]
    motifids = String[]
    motiflengths = Int[]
    nsites = Int[]
    mode = :notreading
    currentmotif = Vector{Vector{Float64}}()
    pwms = Vector{Matrix{Float64}}()

    for line in eachline(file)
        if mode == :notreading
            if startswith(line, "MOTIF")

                fields = split(line)
                motifid   = fields[2]
                motifname = fields[3]

                push!(motifnames, motifname)
                push!(motifids, motifid)

                mode = :reading
            end
        elseif mode == :reading

            if startswith(line, "letter-probability")
                currentmotif = Vector{Vector{Float64}}()
                ldata = replace(line, "letter-probability matrix: " => "")
                fields = split.(split(replace(ldata, "= " => "=")), "=")

                dd = Dict(k => v for (k, v) in fields)
                push!(motiflengths, parse(Int, dd["w"]))
                push!(nsites, parse(Int, dd["nsites"]))

            elseif startswith(line, "URL")
                mode = :notreading
                pwm = reduce(hcat, currentmotif)
                pwm .+= w
                pwm ./= sum(pwm, dims=1)
                push!(pwms, pwm)

            else
                v = parse.(Float64, split(line))
                push!(currentmotif, v)
            end
        end
    end

    motifmeta = DataFrame(MotifName=motifnames, MotifID=motifids, Lengths=motiflengths, NumSites=nsites)
    pbgs = [log2.(pwm) .- log2.(background) for pwm in pwms]

    motifdata = [(name=motifname, id=id, pwm=pwm, pbg=pbg) for (motifname, id, pwm, pbg) in zip(motifnames, motifids, pwms, pbgs)]

    motifmeta, motifdata

end

function loadhomer(file, w = 1e-2, background=background_human_zero_markov())

    io = open(file)
    line = readline(io)

    fields = split(line)
    consensus = fields[1][2:end]
    longname = fields[2]
    name = first(split(longname, "/"))
    sct = parse(Float64, fields[3])

    data = Vector{Vector{Float64}}()
    for line in eachline(io)
        v = parse.(Float64, split(line))
        push!(data, v)
    end

    pwm = reduce(hcat, data)
    pwm .+= w
    pwm ./= sum(pwm, dims=1)

    pbg = log2.(pwm) .- log2.(background)
    (name=name, id=longname, pwm=pwm, pbg=pbg, sct=sct, consensus=consensus)
end


text:
using Metatheory
using Metatheory.EGraphs
using Metatheory.Util
# Description here:
# https://www.philipzucker.com/metatheory-progress/

Metatheory.options.verbose = false
Metatheory.options.printiter = false

# https://github.com/AlgebraicJulia/Catlab.jl/blob/ce2fde9c63a8aab65cf2a7697f43cd24e5e00b3a/src/theories/Monoidal.jl#L127

cat_rules = @theory begin
    f ⋅ id(b) => f
    id(a) ⋅ f => f
    f == f ⋅ id(cod(type(f)))
    f == id(dom(type(f))) ⋅ f

    a ⊗ₒ munit() == a
    munit() ⊗ₒ a == a

    f ⋅ (g ⋅ h) == (f ⋅ g) ⋅ h
end

monoidal_rules = @theory begin
    id(munit()) ⊗ₘ f => f
    f ⊗ₘ id(munit()) => f
    a ⊗ₒ (b ⊗ₒ c) == (a ⊗ₒ b) ⊗ₒ c
    f ⊗ₘ (h ⊗ₘ j) == (f ⊗ₘ h) ⊗ₘ j
    id(a ⊗ₒ b) == id(a) ⊗ₘ id(b)
    (f ⋅ g) ⊗ₘ (p ⋅ q) => (f ⊗ₘ p) ⋅ (g ⊗ₘ q)
end

push!(monoidal_rules,
    MultiPatRewriteRule(@pat((f ⊗ₘ p) ⋅ (g ⊗ₘ q)), @pat((f ⋅ g) ⊗ₘ (p ⋅ q)),
        [PatEquiv(@pat(cod(type(f))), @pat(dom(type(g)))), PatEquiv(@pat(cod(type(p))), @pat(dom(type(q))))])
)

sym_rules = @theory begin
    σ(a, b) ⋅ σ(b, a) == id(a ⊗ₒ b)
    (σ(a, b) ⊗ₘ id(c)) ⋅ (id(b) ⊗ₘ σ(a, c)) == σ(a, b ⊗ₒ c)
    (id(a) ⊗ₘ σ(b, c)) ⋅ (σ(a, c) ⊗ₘ id(b)) == σ(a ⊗ₒ b, c)

    # these rules arer not catlab
    σ(a, munit()) == id(a)
    σ(munit(), a) == id(a)
    σ(munit(), munit()) => id(munit() ⊗ₒ munit())

end

push!(sym_rules,
    MultiPatRewriteRule(@pat((f ⊗ₘ h) ⋅ σ(a, b)), @pat(σ(dom(type(f)), dom(type(h))) ⋅ (h ⊗ₘ f)),
        [PatEquiv(@pat(cod(type(f))), @pat(a)), PatEquiv(@pat(cod(type(h))), @pat(b))]),

    MultiPatRewriteRule(@pat(σ(c, d) ⋅ (h ⊗ₘ f)), @pat((f ⊗ₘ h) ⋅ σ(cod(type(f)), cod(type(h)))),
    [PatEquiv(@pat(dom(type(f))), PatVar(:c)), PatEquiv(@pat(dom(type(f))), PatVar(:d))])
)

diag_rules = @theory begin
    Δ(a) ⋅ (⋄(a) ⊗ₘ id(a)) == id(a)
    Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a)) == id(a)
    Δ(a) ⋅ σ(a, a) == Δ(a)

    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ σ(a, b) ⊗ₘ id(b)) == Δ(a ⊗ₒ b)

    Δ(a) ⋅ (Δ(a) ⊗ₘ id(a)) == Δ(a) ⋅ (id(a) ⊗ₘ Δ(a))
    ⋄(a ⊗ₒ b) == ⋄(a) ⊗ₘ ⋄(b)

    Δ(munit()) == id(munit())
    ⋄(munit()) == id(munit())
end

cart_rules = @theory begin

    pair(f, k) == Δ(dom(type(f))) ⋅ (f ⊗ₘ k)
    proj1(a, b) == id(a) ⊗ₘ ⋄(b)
    proj2(a, b) == ⋄(a) ⊗ₘ id(b)
    f ⋅ ⋄(b) => ⋄(dom(type(f)))
    # Has to invent f. Hard to fix.
    # ⋄(b) => f ⋅ ⋄(b)

    f ⋅ Δ(b) => Δ(dom(type(f))) ⋅ (f ⊗ₘ f)
    Δ(a) ⋅ (f ⊗ₘ f) => f ⋅ Δ(cod(type(f)))
end

push!(cart_rules,
MultiPatRewriteRule(@pat(Δ(a) ⋅ (f ⊗ₘ k)), @pat(pair(f,k)),
[PatEquiv(@pat(dom(type(f))), @pat(dom(type(k))))])
)

typing_rules = @theory begin
    dom(hom(a, b)) => a
    cod(hom(a, b)) => b
    type(id(a)) => hom(a, a)
    type(f ⋅ g) => hom(dom(type(f)), cod(type(g)))
    type(f ⊗ₘ g) => hom(dom(type(f)) ⊗ₒ dom(type(g)), cod(type(f)) ⊗ₒ cod(type(g)))
    type(a ⊗ₒ b) => :ob
    type(munit()) => :ob
    type(σ(a, b)) => hom(a ⊗ₒ b, b ⊗ₒ a)
    type(⋄(a)) => hom(a, munit())
    type(Δ(a)) => hom(a, a ⊗ₒ a)
    type(pair(f, g)) => hom(dom(type(f)), cod(type(f)) ⊗ₒ cod(type(g)))
    type(proj1(a, b)) => hom(a ⊗ₒ b, a)
    type(proj2(a, b)) => hom(a ⊗ₒ b, b)
end

rules = typing_rules ∪ cat_rules ∪ monoidal_rules ∪ sym_rules ∪ diag_rules ∪ cart_rules ∪ typing_rules

# A goofy little helper macro
# Taking inspiration from Lean/Dafny/Agda
using Metatheory.Util
using Metatheory.EGraphs.Schedulers
macro calc(e...)
    theory = eval(e[1])
    e = rmlines(e[2])
    @assert e.head == :block

    trues = Bool[]

    for (a, b) in zip(e.args[1:end-1], e.args[2:end])
        # println(a, " =? ", b)
        params = SaturationParams(
            timeout=12,
            eclasslimit=12000,
            scheduler=BackoffScheduler
            )
        g = EGraph()
        ta = addexpr!(g, :(type(a)))
        tao = addexpr!(g, :(:ob))
        merge!(g, ta.id, tao.id)

        eq = @time areequal(g, theory, a, b; params=params)
        push!(trues, (eq !== missing) && eq)
        println(eq)
        #  WOULD WORK IF COST FUNCTION IS SIMILARITY TO OTHER FUN
        # if !eq
        #     i = 0
        #     while !eq && i < 4
        #         ga = EGraph(a); gb = EGraph(b)
        #         ga_extr = addanalysis!(ga, ExtractionAnalysis, astsize)
        #         gb_extr = addanalysis!(gb, ExtractionAnalysis, astsize)
        #         @time saturate!(ga, theory; timeout = 9)
        #         @time saturate!(gb, theory; timeout = 9)
        #
        #         new_a = extract!(ga, ga_extr)
        #         new_b = extract!(gb, gb_extr)
        #         println("i = $i \nnew a = $new_a \nnew b = $new_b")
        #         eq = @time areequal(theory, new_a, new_b; timeout = 9)
        #         i += 1
        #     end
        #
        #     if !eq && i == 4
        #         return false
        #     end
        # end
    end
    all(trues)
end

@calc rules begin
    id(a ⊗ₒ b)
    id(a) ⊗ₘ id(b)
    (Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a))) ⊗ₘ id(b)
    (Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a))) ⊗ₘ (Δ(b) ⋅ (⋄(b) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ ⋄(a)) ⊗ₘ (⋄(b) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ (⋄(a) ⊗ₘ ⋄(b)) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ ((⋄(a) ⊗ₘ ⋄(b)) ⋅ σ(munit(), munit())) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⋅ id(a)) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ σ(a, b) ⊗ₘ id(b)) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b)))
    Δ(a ⊗ₒ b) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b))
    Δ(a ⊗ₒ b) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b))
    Δ(a ⊗ₒ b) ⋅ (proj1(a, b) ⊗ₘ proj2(a, b))
    pair(proj1(a, b), proj2(a, b))
end

# shorter proof also accepted
@calc rules begin
    id(a ⊗ₒ b)
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    pair(proj1(a, b), proj2(a, b))
end

# shorter proof not quite there
@calc rules begin
    id(a ⊗ₒ b)
    pair(proj1(a, b), proj2(a, b))
end

@calc rules begin
    id(a ⊗ₒ b)
    id(a) ⊗ₘ id(b)
    (Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a))) ⊗ₘ id(b)
    (Δ(a) ⋅ (id(a) ⊗ₘ ⋄(a))) ⊗ₘ (Δ(b) ⋅ (⋄(b) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ ⋄(a)) ⊗ₘ (⋄(b) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ (⋄(a) ⊗ₘ ⋄(b)) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ (id(a) ⊗ₘ ((⋄(a) ⊗ₘ ⋄(b)) ⋅ σ(munit(), munit())) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b)))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⋅ id(a)) ⊗ₘ (σ(a, b) ⋅ (⋄(b) ⊗ₘ ⋄(a))) ⊗ₘ id(b))
    (Δ(a) ⊗ₘ Δ(b)) ⋅ ((id(a) ⊗ₘ σ(a, b) ⊗ₘ id(b)) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b)))
    Δ(a ⊗ₒ b) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b))
    Δ(a ⊗ₒ b) ⋅ (id(a) ⊗ₘ (⋄(b) ⊗ₘ ⋄(a)) ⊗ₘ id(b))
    Δ(a ⊗ₒ b) ⋅ (proj1(a, b) ⊗ₘ proj2(a, b))
    pair(proj1(a, b), proj2(a, b))
end

Metatheory.options.verbose = true
Metatheory.options.printiter = true
Metatheory.options.multithreading = false

G = EGraph( :(pair(proj1(a, b), proj2(a, b))))
params = SaturationParams(timeout=10)
saturate!(G, rules, params )
ex = extract!(G, astsize)

G.classes


text:
module DataEnvelopmentAnalysis

    """
        DataEnvelopmentAnalysis
    A Julia package for efficiency and productivity measurement using Data Envelopment Analysis (DEA).
    [DataEnvelopmentAnalysis repository](https://github.com/javierbarbero/DataEnvelopmentAnalysis.jl).
    """

    using JuMP
    using GLPK
    using Ipopt
    using SparseArrays
    using LinearAlgebra

    using Printf: @sprintf
    using Statistics: std
    using StatsBase: CoefTable

    import StatsBase: nobs, mean

    export
    # optimizer
    DEAOptimizer,
    newdeamodel,

    # Types
    AbstractDEAModel,
    DEAModel,
    AbstractDEAPeers, AbstractDEAPeersDMU,
    DEAPeers, DEAPeersDMU,
    AbstractTechnicalDEAModel, AbstractRadialDEAModel,
    TechnicalDEAModel,
    RadialDEAModel, AdditiveDEAModel, DirectionalDEAModel, GeneralizedDFDEAModel,
    RussellDEAModel, EnhancedRussellGraphDEAModel, ModifiedDDFDEAModel,
    AbstractHolderDEAModel, HolderL1DEAModel, HolderL2DEAModel, HolderLInfDEAModel,
    ReverseDDFDEAModel,
    AbstractEconomicDEAModel,
    AbstractCostDEAModel, AbstractRevenueDEAModel, AbstractProfitDEAModel, AbstractProfitabilityDEAModel,
    CostDEAModel, RevenueDEAModel, ProfitDEAModel, ProfitabilityDEAModel,
    AbstractProductivityDEAModel,
    MalmquistDEAModel,
    # All models
    nobs, ninputs, noutputs,
    # Peers
    peers, peersmatrix, ispeer,
    # Technical models
    dea, deaadd, deaaddweights, deaddf, deagdf,
    dearussell, deaerg, deamddf, deaholder, dearddf,
    efficiency, slacks,
    # Economic models
    deamincost, deamaxrevenue, deamaxprofit,
    deacost, dearevenue, deaprofit, deaprofitability,
    normfactor, ismonetary,
    # Common technical and economic models
    targets,
    # Productivity models
    malmquist,
    nperiods, prodchange

    # Include code of functions
    include("optimizer.jl")

    include("model.jl")
    include("peers.jl")

    include("technical.jl")
    include("dea.jl")
    include("deaadd.jl")
    include("deaddf.jl")
    include("deagdf.jl")
    include("dearussell.jl")
    include("deaerg.jl")
    include("deamddf.jl")
    include("deaholder.jl")
    include("dearddf.jl")

    include("economic.jl")
    include("econoptim.jl")
    include("deacost.jl")
    include("dearevenue.jl")
    include("deaprofit.jl")
    include("deaprofitability.jl")

    include("productivity.jl")
    include("malmquist.jl")

    function __init__()

        nothing
    end

end # module


text:
#####
##### Methods for dealing with PMEM to DRAM ratios
#####

footprint(datum::Dict) = convert(Int, datum[:pmem_alloc_size] + datum[:dram_alloc_size])
footprint(f::nGraph.NFunction) =
    convert(Int, nGraph.get_pmem_pool_size(f) + nGraph.get_temporary_pool_size(f))

function getratio(datum::Dict)
    r = get(datum, :ratio, nothing)
    isnothing(r) || return r

    pmem = convert(Int, datum[:pmem_alloc_size])
    dram = convert(Int, datum[:dram_alloc_size])
    return pmem // dram
end

getratio(fex::nGraph.FluxExecutable) = getratio(fex.ex.ngraph_function)
function getratio(f::nGraph.NFunction)
    pmem = convert(Int, nGraph.get_pmem_pool_size(f))
    dram = convert(Int, nGraph.get_temporary_pool_size(f))
    return pmem // dram
end

ratio_string(x::Rational) = "$(x.num):$(x.den)"
compare_ratio(a, b) = iszero(b.den) ? inv(a) : a - b

max_stars_repo_name: darchr/AutoTM

max_issues_repo_name: darchr/AutoTM

max_forks_repo_name: darchr/AutoTM


text:
using Statistics, SpecialFunctions

"""
    xerfcx_pi(x)

Computes 1/√π - x * erfcx(x), without cancellation.
"""
function xerfcx_pi(x::Real)
    if xerfcx_asym_thresh(x, 3)
        return xerfcx_asym_pi(x, 3)
    else
        return 1/√π - x * erfcx(x)
    end
end

"""
    xerfcx_asym(x, m)

Evaluates x * erfcx(x) by its asymptotic series for large x > 0,
up to m'th order.
"""
function xerfcx_asym(x::Real, m::Integer)
    #= NIST Handbook of Mathematical Functions (2010) 1st Edition, 7.12.1 =#
    @assert m ≥ 0 && x > 0
    s = one(x / 1)
    for j = m : -1 : 0
        #s = 1 - (j + 1/2) / x^2 * s
        s = muladd((j + 1/2) / x / x, -s, 1)
    end
    return s / √π
end

"""
    xerfcx_asym(x, m)

Estimated reminder of xerfcx_asym(x, m).
"""
function xerfcx_asym_rem(x::Real, m::Integer)
    #= NIST Handbook of Mathematical Functions (2010) 1st Edition, 7.12.1. =#
    @assert m ≥ 0 && x > 0
    return 1 / √π * exp(lgamma(m + 1/2) - (2m + 1) * x) / erfcx(x)
end

"""
    xerfcx_asym_pi(x, m)

Computes 1/√π - x * erfcx(x) by its asymptotic series for large x > 0, up to
m'th order.
"""
function xerfcx_asym_pi(x::Real, m::Integer)
    #= NIST Handbook of Mathematical Functions (2010) 1st Edition, 7.12.1 =#
    @assert m ≥ 0 && x > 0
    s = one(x / 1)
    for j = m : -1 : 1
        #s = 1 - s * (j + 1/2) / x^2
        s = muladd((j + 1/2) / x / x, -s, 1)
    end
    s = s / 2 / x / x
    return s / √π
end

"""
    xerfcx_asym_thresh(x, m)

Returns true if the asymptotic series of x * erfcx(x) to order m can be used
for this value of x without appreciable error.
"""
xerfcx_asym_thresh(x::Real, m::Integer) = xerfcx_asym_thresh(x, Val(m))
function xerfcx_asym_thresh(x::T, ::Val{m}) where {T <: Real, m}
    @assert m isa Integer
    @assert m > 0
    return x > (√π * eps(float(T)))^(-1/2m) * exp(lgamma(m + 1/2) / 2m)
end
xerfcx_asym_thresh(x::Union{Float64,Integer}, ::Val{3}) = x > 452
xerfcx_asym_thresh(x::Union{Float32}, ::Val{3}) = x > 16
xerfcx_asym_thresh(x::Union{Float16}, ::Val{3}) = x > 4

"""
    two

Returns 2 of appropriate type. Like `one`, but `two`.
"""
two(::Union{Type{T}, T}) where {T} = convert(T, 2)


text:
function samplepoints(streamlinelength, λ_min , λ_max, n)
    x0 = rand() * 1000.0m
    xs = range(x0, x0 + streamlinelength, length = n)
    no = noise_between_wavelengths(λ_min, λ_max, xs, normalize = false)
end

"""
    z_transform_contributions(vs::T, z::Complex)
        -> Complex
Unilateral z-transform
where
    vs are the samples
    z is a point in the complex plane

"""
function z_transform_contributions(vs::T, z::Complex) where {T<:Union{AbstractRange, Vector}}
    halfn = div(length(vs), 2)
    [s * z^-(i - halfn) for (i, s) in enumerate(vs)]
end
z_transform_contributions_points(vs, z) = [EM * Point(real(ŝ), imag(ŝ)) for ŝ in z_transform_contributions(vs, z)]


text:
append!(empty!(LOAD_PATH), Base.DEFAULT_LOAD_PATH)
using Pkg

# ################################################################################
# # Generate notebooks
# # ################################################################################
# exampledir = joinpath(@__DIR__, "..", "examples")
# Pkg.activate(exampledir)
# Pkg.instantiate()
# include(joinpath(exampledir, "generate_notebooks.jl"))

################################################################################
# Build simulation environments
################################################################################
robotdir = joinpath(@__DIR__, "..", "src/robots")
Pkg.activate(robotdir)
Pkg.instantiate()

using JLD2
using Symbolics
using LinearAlgebra
using Scratch

include("../src/utils.jl")
include("../src/simulator/model.jl")
include("../src/simulator/dimensions.jl")
include("../src/simulator/residual.jl")
include("../src/simulator/codegen.jl")

include("../src/robots/integrator.jl")
include("../src/robots/codegen.jl")

include("../src/solver/cones.jl")

include("../src/robots/hopper/model.jl")  # hopper special case
include("../src/robots/hopper/simulator.jl")

include("../src/robots/particle/model.jl")  # particle special case
include("../src/robots/particle/simulator.jl")

include("../src/robots/centroidal_quadruped/model.jl")  # centroidal quadruped special case
include("../src/robots/centroidal_quadruped/simulator.jl")

include("../src/robots/centroidal_quadruped_ic/model.jl")  # centroidal quadruped special case
include("../src/robots/centroidal_quadruped_ic/simulator.jl")

pkgdir = joinpath(@__DIR__, "..")
Pkg.activate(pkgdir)
path_robots = @get_scratch!("robots")

robots = [
        :hopper, :hopper1, :biped, :quadruped, :box, :particle, :centroidal_quadruped, :centroidal_quadruped_ic,
        :halfcheetah, :halfcheetah4,
        ]

for robot in robots
    # robot model
    include(joinpath("../src/robots", String(robot), "model.jl"))

    # expr path
    path_expr = joinpath(path_robots, String(robot) * ".jld2")

    # kinematics
    contact_kinematics = eval(Symbol(String(robot) * "_contact_kinematics"))
    contact_kinematics_jacobians = eval(Symbol(String(robot) * "_contact_kinematics_jacobians"))

    # codegen
    r_model, rz_model, rθ_model = codegen_residual(eval(robot), codegen_dynamics(eval(robot))..., contact_kinematics, contact_kinematics_jacobians)

    @save path_expr r_model rz_model rθ_model
end

# # robot model
# robot = :box
# include(joinpath("../src/robots", String(robot), "model.jl"))

# # expr path
# path_expr = joinpath(path_robots, String(robot) * ".jld2")

# # kinematics
# contact_kinematics = eval(Symbol(String(robot) * "_contact_kinematics"))
# contact_kinematics_jacobians = eval(Symbol(String(robot) * "_contact_kinematics_jacobians"))

# # codegen
# mass_matrix, dynamics_bias = codegen_dynamics(eval(robot))
# r_model, rz_model, rθ_model = codegen_residual(eval(robot), mass_matrix, dynamics_bias, contact_kinematics, contact_kinematics_jacobians)

# @save path_expr r_model rz_model rθ_model

# # robot model
# robot = :centroidal_quadruped
# include(joinpath("../src/robots", String(robot), "model.jl"))
# include(joinpath("../src/robots", String(robot), "simulator.jl"))

# # expr path
# path_expr = joinpath(path_robots, String(robot) * ".jld2")

# # kinematics
# contact_kinematics = eval(Symbol(String(robot) * "_contact_kinematics"))
# contact_kinematics_jacobians = eval(Symbol(String(robot) * "_contact_kinematics_jacobians"))

# # codegen
# r_model, rz_model, rθ_model = codegen_residual(eval(robot),
#     q -> mass_matrix(eval(robot), q),
#     (q, q̇) -> dynamics_bias(eval(robot), q, q̇),
#     contact_kinematics, contact_kinematics_jacobians)

# @save path_expr r_model rz_model rθ_model