text:
include("utilities.jl")

#selectivity curve using spearman correlation
function spearman_selectivity_filters_curves( data, predictions::Vector{Vector{Float64}};
    plot::Bool=true,
    selectivity=data.θ,
    filter_feature=-data.Δ,
    thresh=compute_thresh(filter_feature),
    labels=[string(i) for i=1:length(predictions)])

    Sf=length(filter_feature)
    n_predictions=size(predictions)[1]

    v_thr_cor=zeros(length(thresh),n_predictions+1)
    for t=1:length(thresh)

            filter_idx = [s for s = 1 : Sf if filter_feature[s] ≥ thresh[t]];

            if !isempty(filter_idx)
                v_thr_cor[t,1] = length(filter_idx)/Sf;
                for k=1:n_predictions
                    v_thr_cor[t,2] = corspearman(predictions[k][filter_idx],selectivity[filter_idx])
                end
            end
    end

    if plot

        for k=1:n_predictions
            semilogx(v_thr_cor[:,1],v_thr_cor[:,k+1],label=labels[k])
        end
        ylabel("Spearman correlation")
        xlabel("fraction of sequences")
        legend()
        title("correlation with mean selectivity")
    else
        return v_thr_cor
    end

end

#plots logselectivity vs energy curves for all train-set splits of the crossvalidation
function plot_fold_curves(data::Vector{Any}, energies::Vector{Vector{Float64}}, range=log_range(1/100,1,20); metric::Symbol = :pearson,
        xlab::String, ylab::String, tit::String)
    nfolds = length(data)
    sel_curves = (metric == :pearson) ? selectivity_filters_curves : spearman_selectivity_filters_curves
    corr_train = [sel_curves(data[k][1], [-energies[k]],
        thresh = compute_thresh(-data[k][1].Δ;
        fraction=range), plot = false) for k in 1:nfolds];
    for k in 1:nfolds
        plot(corr_train[k][:,1], corr_train[k][:,2])
    end
    xlabel(xlab)
    ylabel(ylab)
    title(tit)
    xscale(:log)
    return corr_train
end

export spearman_selectivity_filters_curves, plot_fold_curves


text:
#=

  Euler #36 in Julia.

  Problem 36:
  """
  The decimal number, 585 = 1001001001_(2) (binary), is palindromic
  in both bases.

  Find the sum of all numbers, less than one million, which are palindromic
  in base 10 and base 2.

  (Please note that the palindromic number, in either base, may not
   include leading zeros.)
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

# 0.06355465s
function euler36a()
    s = 0
    for n in 1:999_999
        if palindromic_number(n) &&
            palindromic_list(digits(n,base=2))
            s += n;
        end
    end
    return s
end

# 0.19043701s
function euler36b()
    return ((1:999999)
        .|>n->n*palindromic_number(n)
        .|>n->n*palindromic_list(digits(n,base=2))
        )|>sum
end

function euler36c()
    s = 0
    for n in 1:999_999
        if palindromic_list(digits(n,base=10)) &&
            palindromic_list(digits(n,base=2))
            s += n;
        end
    end
    return s
end

run_euler(euler36a)
# run_euler(euler36b)
# run_euler(euler36a)

ext: jl

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# observation features for the OccludedCrosswalk
# can be either Lidar (defined in ADM)
# or an image representation of the scene

# uncomment for Lidar
function POMDPs.convert_o(::Type{Vector{Float64}}, o::OCObs, pomdp::OCPOMDP)
    return [o.ranges/pomdp.sensor.max_range, o.range_rates/pomdp.env.params.speed_limit]
end

# not used now
mutable struct MapFeatures
    width::Int64
    height::Int64
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
################################################################################
#  Copyright 2020, Tom Van Acker                                               #
################################################################################
# MultiStateSystems.jl                                                         #
# A Julia package to solve multi-state system models.                          #
# See http://github.com/timmyfaraday/MultiStateSystems.jl                      #
################################################################################

"""
Example considering a tree-like wind farm with four wind turbines.

The output generated by the wind turbines depend on a single stochastic source:
the wind speed. To enforce this, the argument `dep` must be passed as true.

  WT₁         WT₂           n₃
   |----///----|----///----|
   n₁    c₁     n₂    c₂     |
                           |--⋅ PCC
  WT₄         WT₃           |
   |----///----|----///----|
   n₅    c₄     n₄    c₃
Fig. 1: Tree-like wind farm with four wind turbines.
"""

# load pkgs
using Unitful
using MultiStateSystems

# include the state-transition diagrams for the wind turbines and cables
stdʷᵗᵒ = solvedSTD(prob = [0.3,0.7], power = [0.0u"MW",2.0u"MW"])
stdᶜᵇˡ = solvedSTD(prob = [0.1,0.9], power = [0.0u"MW",4.0u"MW"])

# initialize the network corresponding to the wind farm
ntwʷᶠ = Network()

# add the user, sources and components to the network
add_user!(ntwʷᶠ, node = 3)
add_sources!(ntwʷᶠ, node = [1,2,4,5], std = stdʷᵗᵒ, dep = true)
add_components!(ntwʷᶠ, edge = [(1,2),(2,3),(4,3),(5,4)],
                       std  = stdᶜᵇˡ)

# return the network
return ntwʷᶠ


text:
module ImageDistortion

include("distortion.jl")

export
    distortpoints,
    undistortpoints,
    getrectangles,
    get_optimalnewcameramatrix,
    init_undistortrectifymap,
    remap

"""
Algorithms:
    - Distortion: `distortpoints`, `undistortpoints`, `get_optimalnewcameramatrix`, `init_undistortrectifymap`, `remap`
"""
end


text:
# ---
# title: 1544. Make The String Great
# id: problem1544
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String, Stack
# link: <https://leetcode.com/problems/make-the-string-great/description/>
# hidden: true
# ---
#
# Given a string `s` of lower and upper case English letters.
#
# A good string is a string which doesn't have  **two adjacent characters**
# `s[i]` and `s[i + 1]` where:
#
#   * `0 <= i <= s.length - 2`
#   * `s[i]` is a lower-case letter and `s[i + 1]` is the same letter but in upper-case or **vice-versa**.
#
# To make the string good, you can choose **two adjacent** characters that make
# the string bad and remove them. You can keep doing this until the string
# becomes good.
#
# Return _the string_ after making it good. The answer is guaranteed to be
# unique under the given constraints.
#
# **Notice** that an empty string is also good.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "leEeetcode"
#     Output: "leetcode"
#     Explanation: In the first step, either you choose i = 1 or i = 2, both will result "leEeetcode" to be reduced to "leetcode".
#
#
# **Example 2:**
#
#
#
#     Input: s = "abBAcC"
#     Output: ""
#     Explanation: We have many possible scenarios, and all lead to the same answer. For example:
#     "abBAcC" --> "aAcC" --> "cC" --> ""
#     "abBAcC" --> "abBA" --> "aA" --> ""
#
#
# **Example 3:**
#
#
#
#     Input: s = "s"
#     Output: "s"
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 100`
#   * `s` contains only lower and upper case English letters.
#
#
## @lc code=start
using LeetCode

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


text:
#
# Custom vector types for representing data in columnar form.
#

using Tables

import Base:
    IndexStyle,
    OneTo,
    convert,
    eltype,
    getindex,
    iterate,
    show,
    size,
    summary,
    &, |, ~

#
# Vector of tuples in columnar form.
#

# Constructors.

"""
    TupleVector([lbls::Vector{Symbol},] len::Int, cols::Vector{AbstractVector})
    TupleVector([lbls::Vector{Symbol},] idxs::AbstractVector{Int}, cols::Vector{AbstractVector})
    TupleVector(lcols::Pair{<:Union{Symbol,AbstractString},<:AbstractVector}...)

Vector of tuples stored as a collection of column vectors.

- `cols` is a vector of columns; optional `lbls` is a vector of column labels.
  Alternatively, labels and columns could be provided as a list of pairs
  `lcols`.
- `len` is the vector length, which must coincide with the length of all the
  columns.  Alternatively, the vector could be constructed from a subset of the
  column data using a vector of indexes `idxs`.
"""
struct TupleVector{I<:AbstractVector{Int}} <: AbstractVector{Any}
    lbls::Vector{Symbol}            # isempty(lbls) means plain tuples
    idxs::I
    cols::Vector{AbstractVector}
    icols::Vector{AbstractVector}   # [col[idxs] for col in cols]

    @inline function TupleVector{I}(lbls::Vector{Symbol}, idxs::I, cols::Vector{AbstractVector}) where {I<:AbstractVector{Int}}
        @boundscheck _checktuple(lbls, idxs, cols)
        icols = Vector{AbstractVector}(undef, length(cols))
        new{I}(lbls, idxs, cols, icols)
    end

    @inline function TupleVector{I}(lbls::Vector{Symbol}, len::Int, cols::Vector{AbstractVector}) where {I<:OneTo{Int}}
        @boundscheck _checktuple(lbls, len, cols)
        idxs = OneTo(len)
        new{I}(lbls, idxs, cols, cols)
    end
end

@inline TupleVector(lbls::Vector{Symbol}, idxs::I, cols::Vector{AbstractVector}) where {I<:AbstractVector{Int}} =
    TupleVector{I}(lbls, idxs, cols)

@inline TupleVector(lbls::Vector{Symbol}, len::Int, cols::Vector{AbstractVector}) =
    TupleVector{OneTo{Int}}(lbls, len, cols)

@inline TupleVector(idxs::AbstractVector{Int}, cols::Vector{AbstractVector}) =
    TupleVector(Symbol[], idxs, cols)

@inline TupleVector(len::Int, cols::Vector{AbstractVector}) =
    TupleVector(Symbol[], len, cols)

@inline TupleVector(len::Int) =
    TupleVector(Symbol[], len, AbstractVector[])

function TupleVector(lcol1::Pair{<:Union{Symbol,AbstractString},<:AbstractVector},
                     more::Pair{<:Union{Symbol,AbstractString},<:AbstractVector}...)
    len = length(last(lcol1))
    lcols = Any[lcol1, more...]
    lbls = Symbol[Symbol(first(lcol)) for lcol in lcols]
    cols = AbstractVector[last(lcol) for lcol in lcols]
    TupleVector(lbls, len, cols)
end

function _checklabels(lbls::Vector{Symbol}, width::Int)
    if !isempty(lbls)
        length(lbls) == width || error("number of labels ≠ number of columns")
        seen = Set{Symbol}()
        for lbl in lbls
            !(lbl in seen) || error("duplicate column label $(repr(lbl))")
            push!(seen, lbl)
        end
    end
end

function _checktuple(lbls::Vector{Symbol}, idxs::AbstractVector{Int}, cols::Vector{AbstractVector})
    _checklabels(lbls, length(cols))
    if !isempty(idxs)
        m = maximum(idxs)
        for col in cols
            length(col) >= m || error("insufficient column height")
        end
    end
end

function _checktuple(lbls::Vector{Symbol}, len::Int, cols::Vector{AbstractVector})
    _checklabels(lbls, length(cols))
    for col in cols
        length(col) == len || error("unexpected column height")
    end
end

# Printing.

show(io::IO, tv::TupleVector) =
    show_vectortree(io, tv)

show(io::IO, ::MIME"text/plain", tv::TupleVector) =
    display_vectortree(io, tv)

# Properties.

@inline labels(tv::TupleVector) = tv.lbls

@inline width(tv::TupleVector) = length(tv.cols)

function columns(tv::TupleVector)
    for j = 1:length(tv.cols)
        if !isassigned(tv.icols, j)
            @inbounds tv.icols[j] = tv.cols[j][tv.idxs]
        end
    end
    tv.icols
end

function column(tv::TupleVector, j::Int)
    if !isassigned(tv.icols, j)
        @inbounds tv.icols[j] = tv.cols[j][tv.idxs]
    end
    tv.icols[j]
end

@inline column(tv::TupleVector, lbl::Symbol) =
    column(tv, findfirst(isequal(lbl), tv.lbls))

@inline locate(tv::TupleVector, j::Int) =
    1 <= j <= length(tv.cols) ? j : nothing

@inline locate(tv::TupleVector, lbl::Symbol) =
    findfirst(isequal(lbl), tv.lbls)

# Vector interface.

@inline size(tv::TupleVector) = size(tv.idxs)

IndexStyle(::Type{<:TupleVector}) = IndexLinear()

eltype(tv::TupleVector) =
    if !isempty(tv.lbls)
        NamedTuple{(tv.lbls...,),Tuple{Any[eltype(col) for col in tv.cols]...}}
    else
        Tuple{Any[eltype(col) for col in tv.cols]...}
    end

@inline function getindex(tv::TupleVector, k::Int)
    @boundscheck checkbounds(tv, k)
    @inbounds k′ = tv.idxs[k]
    @inbounds t = tuple(Any[col[k′] for col in tv.cols]...)
    if !isempty(tv.lbls)
        NamedTuple{(tv.lbls...,)}(t)
    else
        t
    end
end

"""
    (::TupleVector)[ks::AbstractVector{Int}] :: TupleVector

Returns a new `TupleVector` with a subset of rows specified by indexes `ks`.
"""
@inline function getindex(tv::TupleVector, ks::AbstractVector)
    @boundscheck checkbounds(tv, ks)
    @inbounds idxs′ = tv.idxs[ks]
    @inbounds tv′ = TupleVector(tv.lbls, idxs′, tv.cols)
    tv′
end

# Tables.jl export interface.

Tables.istable(tv::TupleVector) = !isempty(tv.lbls)

Tables.columnaccess(::TupleVector) = true

Tables.columns(tv::TupleVector) =
    NamedTuple{(labels(tv)...,)}(Any[collect(col) for col in columns(tv)])

#
# Cardinality of a data block.
#

"""
    x1to1::Cardinality
    x0to1::Cardinality
    x1toN::Cardinality
    x0toN::Cardinality

Cardinality constraints on a block of data.
"""
Cardinality

@enum Cardinality::UInt8 x1to1 x0to1 x1toN x0toN

convert(::Type{Cardinality}, c::Symbol) =
    c == :x1to1 ? x1to1 : c == :x0to1 ? x0to1 : c == :x1toN ? x1toN : c ==:x0toN ? x0toN :
    error("invalid cardinality literal: $c")

quoteof(c::Cardinality) =
    c == x1to1 ? :x1to1 : c == x0to1 ? :x0to1 : c == x1toN ? :x1toN : :x0toN

syntaxof(c::Cardinality) =
    c == x1to1 ? :(1:1) : c == x0to1 ? :(0:1) : c == x1toN ? :(1:N) : :(0:N)

# Bitwise operations.

(~)(c::Cardinality) =
    Base.bitcast(Cardinality, (~UInt8(c))&UInt8(x0toN))

(|)(c1::Cardinality, c2::Cardinality) =
    Base.bitcast(Cardinality, UInt8(c1)|UInt8(c2))

(&)(c1::Cardinality, c2::Cardinality) =
    Base.bitcast(Cardinality, UInt8(c1)&UInt8(c2))

# Predicates.

ismandatory(c::Cardinality) =
    c & x0to1 != x0to1

issingular(c::Cardinality) =
    c & x1toN != x1toN

#
# Vector of data blocks in a columnar form.
#

# Constructors.

"""
    BlockVector(offs::AbstractVector{Int}, elts::AbstractVector, card::Cardinality=x0toN)
    BlockVector(:, elts::AbstractVector, card::Cardinality=x1to1)

Vector of data blocks stored as a vector of elements partitioned by a
vector of offsets.

- `elts` is a continuous vector of block elements.
- `offs` is a vector of indexes that subdivide `elts` into separate blocks.
  Should be monotonous with `offs[1] == 1` and `offs[end] == length(elts)+1`.
  Use `:` if the offset vector is a unit range.
- `card` is the cardinality constraint on the blocks.
"""
struct BlockVector{CARD,O<:AbstractVector{Int},E<:AbstractVector} <: AbstractVector{Any}
    offs::O
    elts::E

    @inline function BlockVector{CARD,O,E}(offs::O, elts::E) where {CARD,O<:AbstractVector{Int},E<:AbstractVector}
        @boundscheck _checkblock(length(elts), offs, CARD)
        new{CARD,O,E}(offs, elts)
    end
end

@inline BlockVector{CARD}(offs::O, elts::E) where {CARD,O<:AbstractVector{Int},E<:AbstractVector} =
    BlockVector{CARD,O,E}(offs, elts)

@inline BlockVector{CARD}(::Colon, elts::AbstractVector) where {CARD} =
    BlockVector{CARD}(OneTo{Int}(length(elts)+1), elts)

@inline BlockVector(offs::AbstractVector{Int}, elts::AbstractVector, card::Cardinality=x1toN|x0to1) =
    BlockVector{card}(offs, elts)

@inline BlockVector(::Colon, elts::AbstractVector, card::Cardinality=x1to1) =
    BlockVector{card}(:, elts)

function _checkblock(len::Int, offs::OneTo{Int}, ::Cardinality)
    !isempty(offs) || error("offsets must be non-empty")
    offs[end] == len+1 || error("offsets must enclose the elements")
end

function _checkblock(len::Int, offs::AbstractVector{Int}, card::Cardinality)
    !isempty(offs) || error("offsets must be non-empty")
    @inbounds off = offs[1]
    off == 1 || error("offsets must start with 1")
    for k = 2:lastindex(offs)
        @inbounds off′ = offs[k]
        off′ >= off || error("offsets must be monotone")
        !issingular(card) || off′ <= off+1 || error("singular blocks must have at most one element")
        !ismandatory(card) || off′ >= off+1 || error("mandatory blocks must have at least one element")
        off = off′
    end
    off == len+1 || error("offsets must enclose the elements")
end

# Printing.

show(io::IO, bv::BlockVector) =
    show_vectortree(io, bv)

show(io::IO, ::MIME"text/plain", bv::BlockVector) =
    display_vectortree(io, bv)

# Properties.

@inline offsets(bv::BlockVector) = bv.offs

@inline elements(bv::BlockVector) = bv.elts

@inline cardinality(bv::BlockVector{CARD}) where {CARD} = CARD

@inline cardinality(::Type{<:BlockVector{CARD}}) where {CARD} = CARD

# Vector interface.

@inline size(bv::BlockVector) = (length(bv.offs)-1,)

IndexStyle(::Type{<:BlockVector}) = IndexLinear()

eltype(bv::BlockVector) =
    Vector{eltype(bv.elts)}

eltype(bv::BlockVector{x0to1}) =
    Union{eltype(bv.elts),Missing}

eltype(bv::BlockVector{x1to1}) =
    eltype(bv.elts)

@inline function getindex(bv::BlockVector, k::Int)
    @boundscheck checkbounds(bv, k)
    @inbounds rng = bv.offs[k]:bv.offs[k+1]-1
    @inbounds elt = bv.elts[rng]
    elt
end

@inline function getindex(bv::BlockVector{x0to1}, k::Int)
    @boundscheck checkbounds(bv, k)
    @inbounds rng = bv.offs[k]:bv.offs[k+1]-1
    @inbounds elt = !isempty(rng) ? bv.elts[rng.start] : missing
    elt
end

@inline function getindex(bv::BlockVector{x1to1}, k::Int)
    @boundscheck checkbounds(bv, k)
    @inbounds elt = bv.elts[k]
    elt
end

"""
    (::BlockVector)[ks::AbstractVector{Int}] :: BlockVector

Returns a new `BlockVector` with a selection of blocks specified by indexes
`ks`.
"""
@inline function getindex(bv::BlockVector, ks::AbstractVector)
    @boundscheck checkbounds(bv, ks)
    _getindex(bv, ks)
end

function _getindex(bv::BlockVector{CARD}, ks::AbstractVector) where {CARD}
    offs′ = Vector{Int}(undef, length(ks)+1)
    @inbounds offs′[1] = top = 1
    i = 1
    @inbounds for k in ks
        l = bv.offs[k]
        r = bv.offs[k+1]
        offs′[i+1] = top = top + r - l
        i += 1
    end
    perm = Vector{Int}(undef, top-1)
    j = 1
    @inbounds for k in ks
        l = bv.offs[k]
        r = bv.offs[k+1]
        copyto!(perm, j, l:r-1)
        j += r - l
    end
    @inbounds elts′ = bv.elts[perm]
    @inbounds bv′ = BlockVector{CARD}(offs′, elts′)
    bv′
end

function _getindex(bv::BlockVector{CARD,OneTo{Int}}, ks::AbstractVector) where {CARD}
    offs′ = OneTo(length(ks)+1)
    @inbounds elts′ = bv.elts[ks]
    @inbounds bv′ = BlockVector{CARD}(offs′, elts′)
    bv′
end

function _getindex(bv::BlockVector{CARD}, ks::OneTo) where {CARD}
    len = length(ks)
    if len == length(bv.offs)-1
        return bv
    end
    @inbounds offs′ = bv.offs[OneTo(len+1)]
    @inbounds top = bv.offs[len+1]
    @inbounds elts′ = bv.elts[OneTo(top-1)]
    @inbounds bv′ = BlockVector{CARD}(offs′, elts′)
    bv′
end

function _getindex(bv::BlockVector{CARD,OneTo{Int}}, ks::OneTo) where {CARD}
    len = length(ks)
    if len == length(bv.offs)-1
        return bv
    end
    offs′ = OneTo(len+1)
    @inbounds elts′ = bv.elts[ks]
    @inbounds bv′ = BlockVector{CARD}(offs′, elts′)
    bv′
end

# Allocation-free view.

mutable struct BlockCursor{T,O<:AbstractVector{Int},E<:AbstractVector{T}} <: AbstractVector{T}
    pos::Int
    l::Int
    r::Int
    offs::O
    elts::E

    @inline BlockCursor{T,O,V}(bv::BlockVector) where {T,O<:AbstractVector{Int},V<:AbstractVector{T}} =
        new{T,O,V}(0, 1, 1, bv.offs, bv.elts)

    @inline function BlockCursor{T,O,V}(pos, bv::BlockVector) where {T,O<:AbstractVector{Int},V<:AbstractVector{T}}
        @boundscheck checkbounds(bv.offs, pos:pos+1)
        @inbounds cr = new{T,O,V}(pos, bv.offs[pos], bv.offs[pos+1], bv.offs, bv.elts)
        cr
    end
end

BlockCursor(bv::BlockVector{CARD,O,E}) where {CARD,T,O<:AbstractVector{Int},E<:AbstractVector{T}} =
    BlockCursor{T,O,E}(bv)

BlockCursor(pos, l, r, bv::BlockVector{CARD,O,E}) where {CARD,T,O<:AbstractVector{Int},E<:AbstractVector{T}} =
    BlockCursor{T,V}(pos, l, r, bv)

# Cursor interface for block vector.

@inline cursor(bv::BlockVector) =
    BlockCursor(bv)

@inline cursor(bv::BlockVector, pos::Int) =
    BlockCursor(pos, bv)

@inline function next!(cr::BlockCursor)
    cr.pos += 1
    cr.l = cr.r
    cr.pos < length(cr.offs) || return false
    cr.r = cr.offs[cr.pos+1]
    true
end

# Vector interface for cursor.

@inline size(cr::BlockCursor) = (cr.r - cr.l,)

IndexStyle(::Type{<:BlockCursor}) = IndexLinear()

@inline function getindex(cr::BlockCursor, k::Int)
    @boundscheck checkbounds(cr, k)
    @inbounds elt = cr.elts[cr.l + k - 1]
    elt
end

@inline function setindex!(cr::BlockCursor, elt, k::Int)
    @boundscheck checkbounds(cr, k)
    @inbounds cr.elts[cr.l + k - 1] = elt
    cr
end

#
# Printing columnar vectors.
#

summary(io::IO, v::Union{TupleVector,BlockVector}) =
    pprint(io, summary_layout(v))

summary_layout(v::AbstractVector) =
    pair_layout(literal("@VectorTree"),
                tile_expr(Expr(:call, :×, length(v), syntaxof(shapeof(v)))),
                sep=" of ")

if VERSION < v"1.5.0-DEV"
Base.typeinfo_prefix(io::IO, cv::Union{TupleVector,BlockVector}) =
    if get(io, :typeinfo, nothing) === nothing
        "@VectorTree $(syntaxof(shapeof(cv))) "
    else
        ""
    end
else
Base.typeinfo_prefix(io::IO, cv::Union{TupleVector,BlockVector}) =
    if get(io, :typeinfo, nothing) === nothing
        "@VectorTree $(syntaxof(shapeof(cv))) ", false
    else
        "", false
    end
end

show_vectortree(io::IO, v::AbstractVector) =
    Base.show_vector(io, v)

function display_vectortree(io::IO, v::AbstractVector)
    if get(io, :compact, false)
        return show_vectortree(io, v)
    end
    summary(io, v)
    !isempty(v) || return
    println(io, ":")
    io = IOContext(io, :typeinfo => eltype(v), :compact => true)
    Base.print_array(io, v)
end

#
# @VectorTree constructor.
#

"""
    @VectorTree sig vec

Constructs a tree of columnar vectors from a plain vector literal.

The first parameter, `sig`, describes the tree structure.  It is defined
recursively:

- Julia type `T` indicates a regular vector of type `T`.
- Tuple `(col₁, col₂, ...)` indicates a `TupleVector` instance.
- Named tuple `(lbl₁ = col₁, lbl₂ = col₂, ...)` indicates a `TupleVector`
  instance with the given labels.
- Prefixes `(0:N)`, `(1:N)`, `(0:1)`, `(1:1)` indicate a `BlockVector` instance
  with the respective cardinality constraints (no constraints, mandatory,
  singular, mandatory+singular).

The second parameter, `vec`, is a vector literal in row-oriented format:

- `TupleVector` data is specified either by a matrix or by a vector of (regular
  or named) tuples.
- `BlockVector` data is specified by a vector of vectors.  A one-element block
  could be represented by its element; an empty block by `missing` literal.
"""
macro VectorTree(sig, vec)
    mk = sig2mk(sig)
    ex = vectorize(mk, vec)
    return esc(ex)
end

const CARD_MAP = Dict(:(0:N) => x0toN,
                      :(0:n) => x0toN,
                      :(1:N) => x1toN,
                      :(1:n) => x1toN,
                      :(0:1) => x0to1,
                      :(1:1) => x1to1)

function sig2mk(sig)
    if Meta.isexpr(sig, :tuple)
        lbls = Symbol[]
        col_mks = MakeAbstractVector[]
        for arg in sig.args
            if Meta.isexpr(arg, (:(=), :kw), 2) && arg.args[1] isa Union{Symbol,String}
                push!(lbls, Symbol(arg.args[1]))
                push!(col_mks, sig2mk(arg.args[2]))
            else
                push!(col_mks, sig2mk(arg))
            end
        end
        return MakeTupleVector(lbls, col_mks)
    elseif Meta.isexpr(sig, :vect, 1)
        elts_mk = sig2mk(sig.args[1])
        return MakeBlockVector(elts_mk, x0toN)
    elseif Meta.isexpr(sig, :call, 3) &&
           sig.args[1] in (:×, :*) && sig.args[2] in keys(CARD_MAP)
        elts_mk = sig2mk(sig.args[3])
        card = CARD_MAP[sig.args[2]]
        return MakeBlockVector(elts_mk, card)
    elseif Meta.isexpr(sig, :call) && length(sig.args) >= 1 && sig.args[1] in keys(CARD_MAP)
        elts_mk = sig2mk(Expr(:tuple, sig.args[2:end]...))
        card = CARD_MAP[sig.args[1]]
        return MakeBlockVector(elts_mk, card)
    else
        ty = sig == :Bottom ? :(Union{}) : sig
        return MakeVector(ty)
    end
end

abstract type MakeAbstractVector end

mutable struct MakeTupleVector <: MakeAbstractVector
    lbls::Vector{Symbol}
    col_mks::Vector{MakeAbstractVector}
    len::Int

    MakeTupleVector(lbls, col_mks) = new(lbls, col_mks, 0)
end

mutable struct MakeBlockVector <: MakeAbstractVector
    elts_mk::MakeAbstractVector
    offs::Vector{Int}
    card::Cardinality
    top::Int

    MakeBlockVector(elts_mk, card) = new(elts_mk, [1], card, 1)
end

mutable struct MakeVector <: MakeAbstractVector
    ty::Any
    vals::Vector{Any}

    MakeVector(ty) = new(ty, [])
end

function vectorize(mk::MakeAbstractVector, ex)
    if Meta.isexpr(ex, (:vect, :vcat))
        for arg in ex.args
            _rearrange!(mk, arg)
        end
        return _reconstruct(mk)
    else
        error("expected a vector literal; got $(repr(ex))")
    end
end

function _rearrange!(mk::MakeTupleVector, ex)
    if Meta.isexpr(ex, (:tuple, :row))
        if length(ex.args) == length(mk.col_mks)
            for (j, (arg, col_mk)) in enumerate(zip(ex.args, mk.col_mks))
                if Meta.isexpr(arg, :(=), 2)
                    if j <= length(mk.lbls) && arg.args[1] == mk.lbls[j]
                        arg = arg.args[2]
                    elseif j < length(mk.lbls)
                        error("expected label $(repr(mk.lbls[j])); got $(repr(arg))")
                    else
                        error("expected no label; got $(repr(arg))")
                    end
                end
                _rearrange!(col_mk, arg)
            end
        else
            error("expected $(length(mk.col_mks)) column(s); got $(repr(ex))")
        end
        mk.len += 1
    elseif length(mk.col_mks) == 1
        _rearrange!(mk.col_mks[1], ex)
        mk.len += 1
    else
        error("expected a tuple or a row literal; got $(repr(ex))")
    end
end

function _rearrange!(mk::MakeBlockVector, ex)
    if Meta.isexpr(ex, (:vect, :vcat))
        for arg in ex.args
            _rearrange!(mk.elts_mk, arg)
            mk.top += 1
        end
    elseif ex !== :missing
        _rearrange!(mk.elts_mk, ex)
        mk.top += 1
    end
    push!(mk.offs, mk.top)
end

function _rearrange!(mk::MakeVector, ex)
    push!(mk.vals, ex)
end

_reconstruct(mk::MakeTupleVector) =
    Expr(:call, TupleVector,
                mk.lbls,
                mk.len,
                Expr(:ref, AbstractVector, Any[_reconstruct(col_mk) for col_mk in mk.col_mks]...))

_reconstruct(mk::MakeBlockVector) =
    Expr(:call, Expr(:curly, BlockVector, mk.card),
                mk.offs == (1:length(mk.offs)) ? :(:) : mk.offs,
                _reconstruct(mk.elts_mk))

_reconstruct(mk::MakeVector) =
    Expr(:ref, :(try $(mk.ty)::Type catch; error("expected a type; got $($(repr(mk.ty)))") end), mk.vals...)


text:
# License for this file: MIT (expat)
# Copyright 2017-2018, DLR Institute of System Dynamics and Control

"""
    module ModiaMath.NonlinearEquations

Solve nonlinear algebraic equations

# Main developer
Martin Otter, [DLR - Institute of System Dynamics and Control](https://www.dlr.de/sr/en)
"""
module NonlinearEquations

export KINSOL, NonlinearEquationsInfo, solveNonlinearEquations!

# using/imports
import ModiaMath

# include code
include("sundials_kinsol.jl")

end


text:
# ------------------------------------------------------------------------------------------
# Follow the video on https://www.youtube.com/watch?v=vAp6nUMrKYg.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# # Autodiff:
 Calculus  from another angle
# (and the special role played by Julia's multiple dispatch and compiler technology)
#
#
#    At the heart of modern machine learning, so popular in (2018),  is an optimization
# problem.  Optimization means gradients, so suddenly differentiation, especially automatic
# differentiation, is exciting.
#
#
#   The first time one  hears about automatic differentiation, it is easy to imagine what it
# is.  Surely it is  straightforward symbolic differentiation applied to code.  One imagines
# automatically doing what is  learned  in a calculus class.
#
#   .... and anyway if it is not that, then it must be finite differences, like one learns
# in a numerical computing class.
#
#
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Babylonian sqrt
#
# We start with a simple example, the computation of sqrt(x), where  how autodiff works
# comes as both a mathematical surprise, and a computing wonder.  The example is  the
# Babylonian algorithm, known to mankind for millenia, to compute sqrt(x):
#
#
#  > Repeat $ t \leftarrow  (t+x/t) / 2 $ until $t$ converges to $\sqrt{x}$.
#
#  Each iteration has one add and two divides. For illustration purposes, 10 iterations
# suffice.
# ------------------------------------------------------------------------------------------

function Babylonian(x; N = 10)
    t = (1+x)/2
    for i = 2:N; t=(t + x/t)/2  end
    t
end

# ------------------------------------------------------------------------------------------
# Check that it works:
# ------------------------------------------------------------------------------------------

α = π
Babylonian(α), √α

x=2; Babylonian(x),√x  # Type \sqrt+<tab> to get the symbol

# Pkg.add(plots)
# Pkg.add(plotly)
using Plots
plotly()
#gr()
#pyplot()

## Warning first plots load packages, takes time
i = 0:.01:49

plot([x->Babylonian(x,N=i) for i=1:5],i,label=["Iteration $j" for i=1:1,j=1:5])

plot!(sqrt,i,c="black",label="sqrt",
      title = "Those Babylonians really knew how to √")

# ------------------------------------------------------------------------------------------
# ## ...and now the derivative, almost by magic
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Eight lines of Julia!  No mention of 1/2 over sqrt(x).
# D for "dual number", invented by the famous algebraist Clifford in 1873.
# ------------------------------------------------------------------------------------------

struct D <: Number  # D is a function-derivative pair
    f::Tuple{Float64,Float64}
end

# ------------------------------------------------------------------------------------------
# Sum Rule: (x+y)' = x' + y'

# Quotient Rule: (x/y)' = (yx'-xy') / y^2
# ------------------------------------------------------------------------------------------

import Base: +, /, convert, promote_rule
+(x::D, y::D) = D(x.f .+ y.f)
/(x::D, y::D) = D((x.f[1]/y.f[1], (y.f[1]*x.f[2] - x.f[1]*y.f[2])/y.f[1]^2))
convert(::Type{D}, x::Real) = D((x,zero(x)))
promote_rule(::Type{D}, ::Type{<:Number}) = D

# ------------------------------------------------------------------------------------------
# The same algorithm with no rewrite at all computes properly
# the derivative as the check shows.
# ------------------------------------------------------------------------------------------

x=49; Babylonian(D((x,1))), (√x,.5/√x)

x=π; Babylonian(D((x,1))), (√x,.5/√x)

# ------------------------------------------------------------------------------------------
# ## It just works!
#
# How does it work?  We will explain in a moment.  Right now marvel that it does.  Note we
# did not
# import any autodiff package.  Everything is just basic vanilla Julia.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## The assembler
#
# Most folks don't read assembler, but one can see that it is short.
# The shortness is a clue that suggests speed!
# ------------------------------------------------------------------------------------------

@inline function Babylonian(x; N = 10)
    t = (1+x)/2
    for i = 2:N; t=(t + x/t)/2  end
    t
end
@code_native(Babylonian(D((2,1))))

# ------------------------------------------------------------------------------------------
# ## Symbolically
#
# We haven't yet explained how it works, but it may be of some value to understand that the
# below is mathematically
# equivalent, though not what the computation is doing.
#
# Notice in the below that Babylonian works on SymPy symbols.
#
# Note: Python and Julia are good friends.  It's not a competition!  Watch how nicely we can
# use the same code now with SymPy.
# ------------------------------------------------------------------------------------------

#Pkg.add("SymPy")
using SymPy

x = symbols("x")
display("Iterations as a function of x")
for k = 1:5
 display( simplify(Babylonian(x,N=k)))
end

display("Derivatives as a function of x")
for k = 1:5
 display(simplify(diff(simplify(Babylonian(x,N=k)),x)))
end

# ------------------------------------------------------------------------------------------
# The code is computing answers mathematically equivalent to the functions above, but not
# symbolically, numerically.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## How autodiff is getting the answer
# Let us by hand take the "derivative" of the Babylonian iteration with respect to x.
# Specifically t′=dt/dx.  This is the old fashioned way of a human rewriting code.
# ------------------------------------------------------------------------------------------

function dBabylonian(x; N = 10)
    t = (1+x)/2
    t′ = 1/2
    for i = 1:N;
        t = (t+x/t)/2;
        t′= (t′+(t-x*t′)/t^2)/2;
    end
    t′

end

# ------------------------------------------------------------------------------------------
# See this rewritten code gets the right answer.  So the trick is for the computer system to
# do it for you, and without any loss of speed or convenience.
# ------------------------------------------------------------------------------------------

x = π; dBabylonian(x), .5/√x

# ------------------------------------------------------------------------------------------
# What just happened?  Answer: We created an iteration by hand for t′ given our iteration
# for t. Then we ran the iteration alongside the iteration for t.
# ------------------------------------------------------------------------------------------

Babylonian(D((x,1)))

# ------------------------------------------------------------------------------------------
# How did this work?  It created the same derivative iteration that we did by hand, using
# very general rules that are set once and need not be written by hand.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Important:: The derivative is substituted before the JIT compiler, and thus efficient
# compiled code is executed.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# ## Dual Number Notation
#
# Instead of D(a,b) we can write a + b ϵ, where ϵ satisfies ϵ^2=0.  (Some people like to
# recall imaginary numbers where an i is introduced with i^2=-1.)
#
# Others like to think of how engineers just drop the O(ϵ^2) terms.
#
# The four rules are
#
# $ (a+b\epsilon) \pm (c+d\epsilon) = (a+c) \pm (b+d)\epsilon$
#
# $ (a+b\epsilon) * (c+d\epsilon) = (ac) + (bc+ad)\epsilon$
#
# $ (a+b\epsilon) / (c+d\epsilon) = (a/c) + (bc-ad)/d^2 \epsilon $
# ------------------------------------------------------------------------------------------

Base.show(io::IO,x::D) = print(io,x.f[1]," + ",x.f[2]," ϵ")

# Add the last two rules
import Base: -,*
-(x::D, y::D) = D(x.f .- y.f)
*(x::D, y::D) = D((x.f[1]*y.f[1], (x.f[2]*y.f[1] + x.f[1]*y.f[2])))

D((1,0))

D((0,1))^2

D((2,1)) ^2

ϵ = D((0,1))
@code_native(ϵ^2)

ϵ * ϵ

ϵ^2

1/(1+ϵ)  # Exact power series:  1-ϵ+ϵ²-ϵ³-...

(1+ϵ)^5 ## Note this just works (we didn't train powers)!!

# ------------------------------------------------------------------------------------------
# ## Generalization to arbitrary roots
# ------------------------------------------------------------------------------------------

function nthroot(x, n=2; t=1, N = 10)
    for i = 1:N;   t += (x/t^(n-1)-t)/n; end
    t
end

nthroot(2,3), ∛2 # take a cube root

nthroot(2+ϵ,3)

nthroot(7,12), 7^(1/12)

x = 2.0
nthroot( x+ϵ,3), ∛x, 1/x^(2/3)/3

# ------------------------------------------------------------------------------------------
# ## Forward Diff
# Now that you understand it, you can use the official package
# ------------------------------------------------------------------------------------------

Pkg.add("ForwardDiff")

using ForwardDiff

ForwardDiff.derivative(sqrt, 2)

ForwardDiff.derivative(Babylonian, 2)

@which ForwardDiff.derivative(sqrt, 2)

# ------------------------------------------------------------------------------------------
# ## Close Look at Convergence with big floats
# the -log10 gives the number of correct digits.  Watch the quadratic convergence right
# before your eyes.
# ------------------------------------------------------------------------------------------

setprecision(3000)
round.(Float64.(log10.([Babylonian(BigFloat(2),N=k) for k=1:10] - √BigFloat(2))),3)

struct D1{T} <: Number  # D is a function-derivative pair
    f::Tuple{T,T}
end

z = D((2.0,1.0))
z1 = D1((BigFloat(2.0),BigFloat(1.0)))

import Base: +, /, convert, promote_rule
+(x::D1, y::D1) = D1(x.f .+ y.f)
/(x::D1, y::D1) = D1((x.f[1]/y.f[1], (y.f[1]*x.f[2] - x.f[1]*y.f[2])/y.f[1]^2))
convert(::Type{D1{T}}, x::Real) where {T} = D1((convert(T, x), zero(T)))
promote_rule(::Type{D1{T}}, ::Type{S}) where {T,S<:Number} = D1{promote_type(T,S)}

A = randn(3,3)

x = randn(3)

ForwardDiff.gradient(x->x'A*x,x)

(A+A')*x

n = 4
Strang = SymTridiagonal(2*ones(n),-ones(n-1))

# ------------------------------------------------------------------------------------------
# ##  But wait there's more!
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Many packages need to be taught how to compute autodiffs of matrix factorications such as
# the svd or lu.  Julia will "just do it," no
# teaching necessary for reasons such as the above.  This is illustrated in another
# notebook, not included here.
# ------------------------------------------------------------------------------------------

max_stars_repo_path: .nbexports/introductory-tutorials/broader-topics-and-ecosystem/intro-to-ml/20. Automatic Differentiation in 10 Minutes.jl

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
# Autogenerated wrapper script for FriBidi_jll for x86_64-w64-mingw32
export fribidi, libfribidi

JLLWrappers.@generate_wrapper_header("FriBidi")
JLLWrappers.@declare_executable_product(fribidi)
JLLWrappers.@declare_library_product(libfribidi, "libfribidi-0.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_executable_product(
        fribidi,
        "bin\\fribidi.exe",
    )

    JLLWrappers.@init_library_product(
        libfribidi,
        "bin\\libfribidi-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
import Compat.String
using DataFrames

type Html
    html::String
end

import Base.writemime
Base.writemime(io::IO, ::MIME"text/html", x::Html) = write(io, x.html)

script(expr) =
    Html(string("<script>", expr, "</script>"))

alert(level, text) =
    Dict("alert" => level, "msg" => text)

teeprint(x) =
    begin println(x); x end # A useful debugging function

get_response_data(x) =
    x.body |> takebuf_string |> JSON.parse

set_metadata(question, obj) =
    script("window.Homework && Homework.set_meta(" * JSON.json(question) * ", " * JSON.json(obj) * ")")

function make_score_dataframe(data, field="score")
    @assert field in ["score", "attempts"]
    d = Dict()
    meta = Dict()

    for q in data["questions"]
        meta[q["id"]] = string(q["max_" * field])
        for s in q["students"]
            student_scores = get(d, s["id"], Dict())
            student_scores[q["id"]] = string(s[field])
            d[s["id"]] = student_scores
        end
    end

    df = DataFrame()

    students = sort(collect(keys(d)))
    names = DataFrame(Names=["MAX", students])

    num(x) = (try int(matchall(r"^[0-9]*", x)[1]) catch 0 end, x)
    questions = sort(collect(keys(meta)), lt=(x, y) -> num(x) < num(y))

    for q in questions
        df[symbol(q)] = [meta[q], [get(get(d, s, Dict()), q, 0) for s in students]]
    end
    if field == "score"
        total = DataFrame(total=[sum([float(v) for (k, v) in row]) for row in DataFrames.eachrow(df)])
    end
    hcat(names, total, df)
end


text:
# Array subindexing:
function Fgetcol(x::Array,incol)
    col=Int(incol[1])
    if col>0
        tmp=zeros(size(x,1),1)
        tmp[:]=x[:,col]
        return (tmp, zeros(size(x)))
    else
        return (zeros(size(x,1)), zeros(size(x)))
    end
end

function Fgetcol_inplace(value,auxvalue,x,incol)
    col=Int(incol[1])
    if col>0
        tmp=zeros(size(x,1),1)
        tmp[:]=x[:,col]
        copy!(value,tmp)
        copy!(auxvalue,zeros(size(x)))
    else
        return
        copy!(value,zeros(size(x,1)))
        copy!(auxvalue,zeros(size(x)))
    end
end

function Dgetcol(derivativeIDX,f_c,faux_c,grad_c,grad_n,x,incol)
    col=Int(incol[1]) # this is ugly since currently we force all values stored on the graph to be real arrays.
    if derivativeIDX==1
        if col>0
            faux_c[:,col]=grad_c
            axpy!(1.0,faux_c,grad_n)
        end
    end
end

#TODO: gpu version

export Fgetcol

Derivative[Fgetcol]=Dgetcol
Inplace[Fgetcol]=Fgetcol_inplace

import Base.*
getindex(A::ADnode,:,i::Int)=ADnode(Fgetcol,[A ADint(i)])
getcol(A::ADnode,i)=ADnode(Fgetcol,[A ADint(i)])

export getindex
export getcol

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
module collision_newtons_cradle

using Modia3D
using Modia3D.StaticArrays
import Modia3D.ModiaMath

vmatGraphics = Modia3D.Material(color="LightBlue" , transparency=0.5)    # material of Graphics
vmatSolids = Modia3D.Material(color="Red" , transparency=0.5)         # material of solids

#c = 1e9, d = 100.0
cmat = "Steel"
Lx = 0.3
Ly = 4.2
Lz = 0.1

@assembly NewtonsCradle begin
  world = Modia3D.Object3D(visualizeFrame=true)
  box    = Modia3D.Object3D(world, Modia3D.Box(Lx, Ly, Lz; material=vmatGraphics) )

  frame1 = Modia3D.Object3D(world; r=[0.0, -Ly/2, 0.0])

  frame2 = Modia3D.Object3D(box; r=[0.0, -Ly/4, 0.0])
  frame3 = Modia3D.Object3D(box; r=[0.0, 0.0, 0.0])
  frame4 = Modia3D.Object3D(box; r=[0.0, Ly/4, 0.0])
  frame5 = Modia3D.Object3D(box; r=[0.0, Ly/2, 0.0])

  pendulum1 = Pendulum()
  pendulum2 = Pendulum()
  pendulum3 = Pendulum()
  pendulum4 = Pendulum()
  pendulum5 = Pendulum()

  rev1   = Modia3D.Revolute(frame1, pendulum1.frame1; axis = 1 , phi_start =  -pi/3)
  rev2   = Modia3D.Revolute(frame2, pendulum2.frame1; axis = 1 , phi_start =  -pi/3)
  rev3   = Modia3D.Revolute(frame3, pendulum3.frame1; axis = 1 ) #, phi_start =  -pi/3)
  rev4   = Modia3D.Revolute(frame4, pendulum4.frame1; axis = 1 , phi_start =  pi/3 ) #
  rev5   = Modia3D.Revolute(frame5, pendulum5.frame1; axis = 1 , phi_start =  pi/3 ) #
end

diameter = 1.0
lengthWire = 4.0

@assembly Pendulum() begin
  wire = Modia3D.Object3D( Modia3D.Cylinder(diameter/5, lengthWire; material=vmatGraphics) )

  frame1 = Modia3D.Object3D(wire; r=[0.0, 0.0, lengthWire/2])
  frame2 = Modia3D.Object3D(wire; r=[0.0, 0.0, -lengthWire/2])

  cyl    = Modia3D.Cylinder(0.15,Lx; material=vmatGraphics)
  cylinder   = Modia3D.Object3D(frame1, cyl; R=ModiaMath.rot2(-pi/2), visualizeFrame=false)

  sphere = Modia3D.Object3D(frame2, Modia3D.Solid(Modia3D.SolidSphere(diameter) , "Steel", vmatSolids; contactMaterial = cmat);   fixed=true ) # , R=ModiaMath.rot2(-pi/2) )

end

gravField = Modia3D.UniformGravityField(g=9.81, n=[0,0,-1])
newton = NewtonsCradle(sceneOptions=Modia3D.SceneOptions(gravityField=gravField,visualizeFrames=true, defaultFrameLength=0.7,nz_max = 100, enableContactDetection=true, visualizeContactPoints=true, visualizeSupportPoints=true))

# Modia3D.visualizeAssembly!( newton )

model = Modia3D.SimulationModel( newton )
result = ModiaMath.simulate!(model; stopTime=10.0, tolerance=1e-8,interval=0.001, log=true)

println("... success of collision_newtons_cradle.jl!")
end


text:
using Test

include("robot-name.jl")

# Random names means a risk of collisions.
history = String[]

isname(x) = occursin(r"^[A-Z]{2}[0-9]{3}$", x)

@testset "one robot" begin
    global r1 = Robot()
    push!(history, r1.name)

    @testset "name of robot is valid" begin
        @test isname(r1.name)
    end

    @testset "names of robot instance are valid and unique in history" for i=1:100
        reset!(r1)
        @test isname(r1.name)
        @test !in(r1.name, history)
        push!(history, r1.name)
    end
end

@testset "two robots" begin
    global r2 = Robot()
    global r3 = Robot()

    @testset "names of robots are valid" begin
        @test isname(r2.name)
        @test isname(r3.name)
    end

    @testset "names of robots are not equal" begin
        @test r2.name != r3.name
    end

    @testset "names of both robots are unique in history" begin
        @test !in(r2.name, history)
        @test !in(r3.name, history)
    end

    push!(history, r2.name)
    push!(history, r3.name)
end

@testset "many robots" begin
    global robots = Robot[]

    for i=1:10
        push!(robots, Robot())

        @testset "name of robot is valid and unique in history" begin
            @test isname(robots[i].name)
            @test !in(robots[i].name, history)
        end

        push!(history, robots[i].name)
    end

    @testset "fresh names of reset robots are valid and unique in history" for r in robots
        reset!(r)
        @test isname(r.name)
        @test !in(r.name, history)
        push!(history, r.name)
    end
end


text:
using Statistics
using Random
using Distributions
import Future
import EvaluationOfRLAlgs

using Plots
gr()

# sample uniformly on natural log scale
function logRand(low, high, rng)
    X = exp(rand(rng, Uniform(log(low), log(high))))
    logp = -log(log(high) - log(high) - log(X))
    return X, logp
end

function quickeval()
    num_sweeps = 25
    results = zeros((num_sweeps,2))
    rng = Random.MersenneTwister(0)  # change 0 to change random seed

	# use this if you have set Blas to single threaded
    # rngs = [rng; accumulate(Future.randjump, fill(big(10)^20, num_sweeps-1), init=rng)]  # use different random number generator for each trial
	# Threads.@threads for s in 1:num_sweeps
    #     results[s, :] .= run_alg(rngs[s])
    # end

	# use this to get results sequentially
	for s in 1:num_sweeps
		results[s, :] .= run_alg(rng)
	end

	# get results for empirical quantile plots
    lrets = sort(results[:, 1])
	frets = sort(results[:, 2])
	xpts = collect(Float64, 1:num_sweeps) / num_sweeps
	# plot empirical quantile plot for both distributions
	p = plot(xpts, lrets, label="Average Return", color=:dodgerblue)
	p = plot!(p, xpts, frets, label="Final Return", color=:crimson)
	ylabel!(p, "Return")
	xlabel!(p, "Probability")
	display(p)
end

function run_alg(rng)
	env = MountainCar(Float64, Int, false) # takes about 20 seconds
	# env = PinBall(Float64, "pinball_medium.cfg") # takes about 2 minutes
    num_trials = 1
    num_episodes = 100

    state_ranges = get_state_desc(env)
    num_actions = get_action_desc(env)
    ϕ = FourierBasis(state_ranges, 5,6,false)  		# create a fourier basis function that use both sine and cosine functions with dependent order of 6 and independent order of 7
    p = LinearFunction(Float64, ϕ, num_actions)  	# create a linear function using the fourier basis to use as the policy
    policy = LinearSoftmaxPolicy(p)  				# make the policy a softmax of the function p
    vf = LinearFunction(Float64, ϕ, 1)  			# create a linear function for the value function
	qf = EvaluationOfRLAlgs.LinearFunction(Float64, ϕ, num_actions)  # creates a q function for use in Sarsa and Q learning algorithms
    λ = rand(rng, Uniform(0., 1.0))  				# sample a random elgibility trace decay
    γ = get_gamma(env)  							# get discount factor from environment (most are 1.0 so setting this lower could be useful)
	ϵ = rand(rng, Uniform(0.0, 0.1))  				# sample the #ϵ-greedy exploration parameter
	αp = logRand(1e-3, 1.0, rng)[1] / get_num_params(policy) # scale the learning rate by the number of parameters. Helps for setting linear rates of linear functions.
	#αv = logRand(1e-3, 1.0, rng)[1] / get_num_params(vf)  # learning rate for value functions when not using Parl2 (i.e., REINFORCE)

	# create the optimization functions
    popt = AccumulatingTraceOptimizer(policy, αp, γ, λ)
    vopt = Parl2AccumulatingTraceOptimizer(vf, γ, λ)
	# qopt = Parl2AccumulatingTraceOptimizer(qf, γ, λ)

	# create the algorithm
    agent = ActorCritic(policy, vf, popt, vopt, γ)
    # agent = NACTD(policy, vf, wf, popt, vopt, wopt, γ, true)  # exercise: create linear functon wf using IdentityBasis
    # agent = REINFORCE(policy, vf, αp, αv, γ, false, true)
    # agent = Sarsa(qf, qopt, γ, ϵ)
    # agent = QLearning(qf, qopt, γ, ϵ)

	# finite difference based search methods (black box optimization)
    # agent = BasicRandomSearch(policy, 0.01, 1.0, 1.0, 3, rng, num_episodes-10)
    # agent = AugmentedRandomSearch(policy, 0.015, 1.0, 0.5, 10, 5, rng, num_episodes-10)

    # returns = run_agent_parallel(env, agent, num_trials, num_episodes, rng)  # run multiple trials using several threads
    # returns = run_agent(env, agent, num_trials, num_episodes, rng) # run multiple trials using one thread
    learning_returns = run_agent!(env, agent, num_episodes, rng) # run one trial. the ! denotes that agent will be update. The two functions above do not update the agent.
	# returns = run_agent!(env, agent, num_episodes, rng, true)  # the true flag displays the agent during leraning. PyPlot backend is slow so use gr(). Do not run this in parallel or if you want speed.
	final_returns = eval_agent(env, agent, 30, rng) # evaluate updated agent after learning

	return mean(learning_returns), mean(final_returns)

end

max_stars_repo_name: ScottJordan/EvaluationOfRLAlgs

max_issues_repo_name: ScottJordan/EvaluationOfRLAlgs

max_forks_repo_name: ScottJordan/EvaluationOfRLAlgs


text:
###############################################################################
#
#   Rational.jl : Additional Nemo functionality for Julia Rationals
#
###############################################################################

###############################################################################
#
#   Data type and parent object methods
#
###############################################################################

JuliaQQ = Rationals{BigInt}()

qq = Rationals{Int}()

parent(a::Rational{T}) where T <: Integer = Rationals{T}()

elem_type(::Type{Rationals{T}}) where T <: Integer = Rational{T}

parent_type(::Type{Rational{T}}) where T <: Integer = Rationals{T}

base_ring(a::Rational{Int}) = zz

base_ring(a::Rational{BigInt}) = JuliaZZ

base_ring(a::Rationals{Int}) = zz

base_ring(a::Rationals{BigInt}) = JuliaZZ

base_ring(a::Rationals{T}) where T <: Integer = Integers{T}()

base_ring(a::Rational{T}) where T <: Integer = Integers{T}()

isexact_type(::Type{Rational{T}}) where T <: Integer = true

isdomain_type(::Type{Rational{T}}) where T <: Integer = true

###############################################################################
#
#   Basic manipulation
#
###############################################################################

zero(::Rationals{T}) where T <: Integer = Rational{T}(0)

one(::Rationals{T}) where T <: Integer = Rational{T}(1)

if VERSION < v"0.7.0-DEV.1144"
isone(a::Rational) = a == 1
end

isunit(a::Rational) = a != 0

canonical_unit(a::Rational)  = a

###############################################################################
#
#   String I/O
#
###############################################################################

function show(io::IO, R::Rationals)
   print(io, "Rationals over ")
   show(io, base_ring(R))
end

needs_parentheses(::Rational) = false

isnegative(a::Rational) = a < 0

show_minus_one(::Type{Rational{T}}) where T <: Integer = false

###############################################################################
#
#   Divides
#
###############################################################################

function divides(a::Rational{T}, b::Rational{T}) where T <: Integer
   return true, a//b
end

###############################################################################
#
#   Exact division
#
###############################################################################

divexact(a::Rational, b::Integer) = a//b

divexact(a::Integer, b::Rational) = a//b

divexact(a::Rational, b::Rational) = a//b

###############################################################################
#
#   GCD
#
###############################################################################

function gcd(p::Rational{T}, q::Rational{T}) where T <: Integer
   a = p.num*q.den
   b = p.den*q.num
   n = gcd(a, b)
   d = p.den*q.den
   if d != 1 && n != 0
      g = gcd(n, d)
      n = divexact(n, g)
      d = divexact(d, g)
   end
   if n == 0
      return Rational{T}(n, T(1))
   else
      return Rational{T}(n, d)
   end
end

###############################################################################
#
#   Exponential
#
###############################################################################

function exp(a::Rational{T}) where T <: Integer
   a != 0 && throw(DomainError())
   return Rational{T}(1)
end

###############################################################################
#
#   Unsafe functions
#
###############################################################################

function zero!(a::Rational{T}) where T <: Integer
   n = a.num
   n = zero!(n)
   if a.den == 1
      return Rational{T}(n, a.den)
   else
      return Rational{T}(n, T(1))
   end
end

function mul!(a::Rational{T}, b::Rational{T}, c::Rational{T}) where T <: Integer
   n = a.num
   d = a.den
   n = mul!(n, b.num, c.num)
   d = mul!(d, b.den, c.den)
   if d != 1 && n != 0
      g = gcd(n, d)
      n = divexact(n, g)
      d = divexact(d, g)
   end
   if n == 0
      return Rational{T}(n, T(1))
   else
      return Rational{T}(n, d)
   end
end

function add!(a::Rational{T}, b::Rational{T}, c::Rational{T}) where T <: Integer
   if a === b
      return addeq!(a, c)
   elseif a == c
      return addeq!(a, b)
   else # no aliasing
      n = a.num
      d = a.den
      d = mul!(d, b.den, c.den)
      n = mul!(n, b.num, c.den)
      n = addmul!(n, b.den, c.num)
      if d != 1 && n != 0
         g = gcd(n, d)
         n = divexact(n, g)
         d = divexact(d, g)
      end
      if n == 0
         return Rational{T}(n, T(1))
      else
         return Rational{T}(n, d)
      end
   end
end

function addeq!(a::Rational{T}, b::Rational{T}) where T <: Integer
   if a === b
      if iseven(a.den)
         return Rational{T}(a.num, div(b.den, 2))
      else
         return Rational{T}(2*a.num, b.den)
      end
   else
      n = a.num
      n = mul!(n, n, b.den)
      n = addmul!(n, b.num, a.den)
      d = a.den
      d = mul!(d, d, b.den)
      if d != 1 && n != 0
         g = gcd(n, d)
         n = divexact(n, g)
         d = divexact(d, g)
      end
      if n == 0
         return Rational{T}(n, T(1))
      else
         return Rational{T}(n, d)
      end
   end
end

function addmul!(a::Rational{T}, b::Rational{T}, c::Rational{T}, d::Rational{T}) where T <: Integer
   d = mul!(d, b, c)
   a = addeq!(a, d)
   return a
end

###############################################################################
#
#   Random generation
#
###############################################################################

function rand(R::Rationals{T}, n::UnitRange{Int}) where T <: Integer
   d = T(0)
   while d == 0
      d = T(rand(n))
   end
   n = T(rand(n))
   return R(n, d)
end

###############################################################################
#
#   Parent object call overload
#
###############################################################################

function (R::Rationals{T})() where T <: Integer
   return Rational{T}(0)
end

function (R::Rationals{T})(b) where T <: Integer
   return Rational{T}(b)
end

function (R::Rationals{T})(b::Integer, c::Integer) where T <: Integer
   return Rational{T}(b, c)
end

###############################################################################
#
#   FractionField constructor
#
###############################################################################

FractionField(R::Integers{T}) where T <: Integer = Rationals{T}()


text:
# LZ78 Prediction Implementation

"""
    LZ78{SymbolType}()

Creates a LZ78 predictor for `SymbolType`. `SymbolType` can be any valid type including
`Char`, `Int` etc.,

## Examples
```julia-repl
julia> p = LZ78{Int}()
DiscretePredictors.LZ78{Int}([*] (0)
, Int[], Int[])

julia> p = LZ78{Char}()
DiscretePredictors.LZ78{Char}([*] (0)
, Char[], Char[])

julia> p = LZ78{String}()
DiscretePredictors.LZ78{String}([*] (0)
, String[], String[])
```

Reference:
Ziv, Jacob, and Abraham Lempel. "Compression of individual sequences via variable-rate coding." IEEE transactions on Information Theory 24.5 (1978): 530-536.
"""
mutable struct LZ78{T} <: BasePredictor{T}
    model::Trie{T,Int}
    phrase::Vector{T}
    context::Vector{T}

    # Constructor
    LZ78{T}() where {T} = new( Trie{T,Int}(), Vector{T}(), Vector{T}() )
    # (::Type{LZ78{T}}){T}() = new{T}( Trie{T,Int}(), Vector{T}(), Vector{T}() );
end

function add!( p::LZ78{T}, sym::T ) where {T}
    p.model.value += 1

    # Update model
    push!( p.phrase, sym )
    if( p.phrase ∈ keys(p.model) ) # Already in dictionary
        # Meww!!
        p.context   = p.phrase[1:end]
        p.model[p.phrase] += 1
    else
        # keyBuffer   = copy( p.phrase );
        p.model[p.phrase]  = 1 # Add this phrase to trie
        # keyBuffer = p.phrase[1:end-1];  # Suffix
        # while !isempty(keyBuffer)
        #     p.model[ keyBuffer ] += 1;
        #     pop!( keyBuffer );
        # end
        # println( "Spliting phrase: ", p.phrase )
        p.context = p.phrase[2:end]
        p.phrase  = Vector{T}()
    end

    nothing
end

function predict( p::LZ78{T} ) where {T}
    # Create a dictionary with symbols
    # println( "Predicting on context: ", p.context )
    children_of_root    = children(p.model,Vector{T}())
    if isempty(children_of_root) # No symbol seen
        return Dict{T,Float64}()   # No prediction
    end
    root_sum            = mapreduce( c->p.model[[c]], +, keys(children_of_root) )
    symbols = Dict( k => (p.model[[k]]/root_sum) for k ∈ keys(children_of_root) )

    window  = p.context[1:end]
    buffer  = Vector{T}()
    for i = length(window):-1:1
        pushfirst!( buffer, window[i] )
        list_of_children    = children( p.model, buffer )
        # Get sum of all children values
        s   = isempty(list_of_children) ?
                0 :
                mapreduce( k->list_of_children[k].value, +, keys(list_of_children) );
        # Apply escape probability
        if haskey(p.model,buffer)
            esc_prob = (p.model[buffer]-s)/p.model[buffer]
        else
            esc_prob = 1.0;
        end
        # println( "    escape from cxt: ", buffer, " with ", esc_prob )
        for symbol in keys(symbols)
            symbols[symbol] *= esc_prob
        end
        # Get each child probability
        for k ∈ keys(list_of_children)
            # println( "        ", k, " -> ", list_of_children[k].value )
            if haskey(symbols,k)    # Not an orphan symbol
                symbols[k] += list_of_children[k].value/p.model[buffer]
            else
                symbols[k] = 0.00
            end
        end
    end

    return symbols;
end

function info_string( p::LZ78{T} ) where {T}
    return @sprintf( "LZ78" )
end

function unique_string( p::LZ78{T} ) where {T}
    return @sprintf( "LZ78" )
end


text:
#Getting overlapping heatmaps
using GLMakie, Colors

xs = LinRange(0, 10, 100)
ys = LinRange(0, 15, 100)
zs = [cos(x) * sin(y) for x in xs, y in ys]

heatmap(xs,ys,zs)

as = LinRange(0, 10, 100)
bs = LinRange(0, 15, 100)
cs = [exp(-((x-5)^2/2 + (y-7.5)^2/2)) for x in xs, y in ys]

c1 = RGBA(0,0,0,0.0);
c2 = RGBA(1,0,0,1.0);

cm = range(c1, stop=c2, length=100);

heatmap!(as,bs,cs,colormap = cm)


text:
@testset "Models" begin
    @testset "Scheffer" begin
        setup = lake_initial_state(1.7, 50.0, 30.0, Scheffer)
        @test setup[1] == [50.0, 30.0]
        @test setup[2].nutrients == 1.7

        lp = LakeParameters(Scheffer, 0.9)
        @test typeof(lp) <: LimnoSES.SchefferParameters
        @test lp.nutrients == 0.9

        @test LimnoSES.preset_conditions(Clear, Scheffer) == (0.7, [20.0, 1.8])
        @test LimnoSES.preset_conditions(Turbid, Scheffer) == (2.5, [84.0, 0.04])

        @test sprint(show, "text/plain", lp) ==
              "Parameters for lake dynamics (Scheffer) model):\nNutrient level: 0.9\nImmigration rate (g⋅m⁻²⋅day⁻¹) for bream: 2.0e-5, pike: 2.0e-5\nGrowth rate (day⁻¹) for bream: 0.0075, Predation rate (day⁻¹) of pike: 0.05\nHalf sturation constants: H₁ 0.5, H₂ 10.0%, H₃ 20.0 g⋅m⁻², H₄ 15.0 g⋅m⁻²\nIntraspecific competition constant (g⋅m⁻²⋅day⁻¹) for bream: 7.5e-5, pike: 0.000275\nPike food conversion efficiency to growth: 0.1, Mortalitiy rate (day⁻¹): 0.00225\nMaximum vegetation coverage: 100.0%\n"

        du = [0.0, 0.0]
        LimnoSES.lake_dynamics!(du, [20.0, 1.8], lp, 0.0)
        @test du == [0.008848571428571449, -0.0001209999999999985]
    end
    @testset "Martin" begin
        setup = lake_initial_state(1.7, 50.0, 30.0, 80.0, Martin)
        @test setup[1] == [50.0, 30.0, 80.0]
        @test setup[2].nutrients == 1.7

        lp = LakeParameters(Martin, 0.9)
        @test typeof(lp) <: LimnoSES.MartinParameters
        @test lp.nutrients == 0.9

        @test LimnoSES.preset_conditions(Clear, Martin) == (0.7, [20.5172, 1.7865, 56.8443])
        @test LimnoSES.preset_conditions(Turbid, Martin) ==
              (2.5, [83.0128, 0.0414705, 6.40048])
        @test LimnoSES.preset_conditions(X1, Martin) == (2.2, [36.730, 2.87725, 26.6798])
        @test LimnoSES.preset_conditions(X2, Martin) == (1.05, [59.0606, 0.819124, 12.0023])
        @test LimnoSES.preset_conditions(X3, Martin) == (1.05, [64.0559, 0.374008, 10.3631])
        @test LimnoSES.preset_conditions(S1, Martin) == (2.0, [79.597, 0.050, 6.928])
        @test LimnoSES.preset_conditions(S2, Martin) == (3.5, [85.436, 0.0373, 6.061])
        @test LimnoSES.preset_conditions(S3, Martin) == (3.5, [20.0, 3.5, 50.0])
        @test LimnoSES.preset_conditions(T1, Martin) == (0.9, [21.647, 2.050, 53.726])
        @test LimnoSES.preset_conditions(T2, Martin) == (2.0, [26.901, 2.809, 41.530])
        @test LimnoSES.preset_conditions(T3, Martin) == (3.5, [20.0, 3.5, 50.0])

        @test sprint(show, "text/plain", lp) ==
              "Parameters for lake dynamics (Martin) model):\nNutrient level: 0.9\nImmigration rate (g⋅m⁻²⋅day⁻¹) for bream: 2.0e-5, pike: 2.0e-5\nGrowth rate (day⁻¹) for bream: 0.0075, Predation rate (day⁻¹) of pike: 0.05\nHalf sturation constants: H₁ 0.5, H₂ 11.0%, H₃ 20.0 g⋅m⁻², H₄ 15.0 g⋅m⁻²\nIntraspecific competition constant (g⋅m⁻²⋅day⁻¹) for bream: 7.5e-5, pike: 0.000275\nPike food conversion efficiency to growth: 0.1, Mortalitiy rate (day⁻¹): 0.00225\nVegetation rates. Growth: 0.007 day⁻¹, mortality: 0.007 day⁻¹, competition: 6.0e-5 m²\n"

        du = [0.0, 0.0, 0.0]
        LimnoSES.lake_dynamics!(du, [20.0, 1.8, 56.8], lp, 0.0)
        @test du == [0.008848571428571449, -9.551327433628197e-5, 0.005225600000000025]
    end
end


text:
@testset "Model Access" begin
@testset "Accessing model" begin
    model = ABM(Agent0; properties = Dict(:a => 2, :b => "test"))
    for i in 1:5
        add_agent!(model)
    end
    @test model.scheduler == schedule_fastest
    @test typeof(model.agents) <: Dict
    @test model.space === nothing
    @test model.properties == Dict(:a => 2, :b => "test")
    a = model[1]
    @test a isa Agent0
    @test a.id == 1
    @test model.a == 2
    @test model.b == "test"

    model.a = 7
    model.b = "Changed"
    @test model.a == 7
    @test model.b == "Changed"
    @test_throws ErrorException model.c = 5

    newa = Agent0(6)
    model[6] = newa
    @test model[6] == newa
    @test_throws ArgumentError model[7] = newa

    prop2 = Agent2(1, 0.5)
    model2 = ABM(Agent0; properties = prop2)
    @test model2.weight == 0.5

    mutable struct Properties
        par1::Int
        par2::Float64
        par3::String
    end
    properties = Properties(1,1.0,"Test")
    model = ABM(Agent0; properties = properties)
    @test model.par1 == 1
    @test model.par2 == 1.0
    @test model.par3 == "Test"

    model.par1 = 7
    model.par2 = 7
    model.par3 = "Changed"
    @test model.par1 == 7
    @test model.par2 == 7.0
    @test model.par3 == "Changed"
    @test_throws ErrorException model.par4 = 5

    model = ABM(Agent0)
    @test_throws ErrorException model.a = 5
end

@testset "model access typestability" begin
    prop1 = Dict(:a => 0.5)
    prop2 = Agent2(1, 0.5)
    model1 = ABM(Agent0; properties = prop1)
    model2 = ABM(Agent0; properties = prop2)

    test1(model) = model.a
    test2(model) = model.weight
    @test_nowarn @inferred test1(model1)
    @test_nowarn @inferred test2(model2)
end

@testset "Core methods" begin
    model = ABM(Agent3, GridSpace((5,5)))
    @test Agents.agenttype(model) == Agent3
    @test Agents.spacetype(model) <: GridSpace
    @test size(model.space) == (5,5)
    @test all(isempty(p, model) for p in positions(model))
end

@testset "Display" begin
    model = ABM(Agent0)
    @test sprint(show, model)[47:54] == "no space"
    model = ABM(Agent3, GridSpace((5,5)))
    @test sprint(show, model)[1:29] == "AgentBasedModel with 0 agents"
    @test sprint(show, model.space) == "GridSpace with size (5, 5), metric=chebyshev, periodic=true"
    model = ABM(Agent6, ContinuousSpace((1,1), 0.1))
    @test sprint(show, model.space) == "periodic continuous space with 10×10 divisions"
    model = ABM(Agent5, GraphSpace(path_graph(5)))
    @test sprint(show, model.space) == "GraphSpace with 5 positions and 4 edges"
end

end


text:
@testset " IncidenceTable" begin
    ex1 = :(foo = bing + bong)
    ex2 = :(cowboy = yee - haw! + foo)
    ex3 = :(x = foo(-3) + bing(1))
    exs = [ex1, ex2, ex3]

    # visiting symbols
    it = Dolang.IncidenceTable()
    Dolang.visit!(it, :x, 1, 0)
    want_eq = OrderedDict(1=>Dict(:x=>Set(0)))
    @test want_eq == it.by_eq
    @test Dict(:x=>Set(0)) == it.by_var

    Dolang.visit!(it, :x, 1, 1)
    want_eq[1][:x] = Set([0, 1])
    @test want_eq == it.by_eq
    @test Dict(:x=>Set([0, 1])) == it.by_var

    Dolang.visit!(it, :x, 2, 3)
    want_eq[1][:x] = Set([0, 1])
    want_eq[2] = Dict(:x=>Set(3))
    @test want_eq == it.by_eq
    @test Dict(:x=>Set([0, 1, 3])) == it.by_var

    Dolang.visit!(it, :y, 1, -1)
    want_eq[1][:y] = Set([-1])
    @test want_eq == it.by_eq
    @test Dict(:x=>Set([0, 1, 3]), :y=>Set([-1])) == it.by_var

    # visiting numbers shouldn't change anything
    for i in 1:100
        Dolang.visit!(it, i, 1, 1)
    end
    @test want_eq == it.by_eq
    @test Dict(:x=>Set([0, 1, 3]), :y=>Set([-1])) == it.by_var

    # reset `it` before working on expressions
    it = Dolang.IncidenceTable()
    want_eq = deepcopy(it.by_eq)
    want_bv = deepcopy(it.by_var)

    s0 = Set(0)
    Dolang.visit!(it, ex1, 1, 0)
    want_eq[1] = Dict(:foo=>s0, :bing=>s0, :bong=>s0)
    want_bv = deepcopy(want_eq[1])
    @test want_eq == it.by_eq
    @test want_bv == it.by_var

    Dolang.visit!(it, ex2, 2, 0)
    want_eq[2] = Dict(:cowboy=>s0, :yee=>s0, :haw! =>s0, :foo=>s0)
    # want_bv[:cowboy] =s0; want_bv[:yee] = s0; want_bv[:haw!] =s0
    merge!(want_bv, want_eq[2])
    @test want_eq == it.by_eq
    @test want_bv == it.by_var

    Dolang.visit!(it, ex3, 3, 0)
    want_eq[3] = Dict(:x=>s0, :foo=>Set([-3]), :bing=>Set(1))
    want_bv[:x] = s0
    want_bv[:foo] = Set([-3, 0])
    want_bv[:bing] = Set([0, 1])
    @test want_eq == it.by_eq
    @test want_bv == it.by_var

    Dolang.visit!(it, ex3, 4, 1)
    want_eq[4] = Dict(:x=>Set(1), :foo=>Set([-2]), :bing=>Set(2))
    want_bv[:x] = Set([0, 1])
    want_bv[:foo] = Set([-3, -2, 0])
    want_bv[:bing] = Set([0, 1, 2])
    @test want_eq == it.by_eq
    @test want_bv == it.by_var

    # using expressions
    it1 = Dolang.IncidenceTable(ex1)
    want = OrderedDict(1=>Dict(:foo=>Set(0), :bing=>Set(0),
                               :bong=>Set(0)))
    @test it1.by_eq == want

    want = Dict(:foo=>Set(0), :bing=>Set(0), :bong=>Set(0))
    @test it1.by_var == want

    # getindex
    @test it[1] == it.by_date[1]
    @test it[:foo] == it.by_var[:foo]
end

max_issues_repo_name: EconForge/Dolang


text:
"""
    This might not be the right way to go about it to solve the problem for small firms.
"""
function residuals_unc(x; par::Pars_v3, type::Symbol)
    l_P = x[1]^2
    l_G = x[2]^2
    D_E = x[3]^2
    D_I = x[4]^2

    D = bundle(D_I, D_E; par, type)
    l = labor_demand(l_P, D; par, type)

    #! FOCs
    res = zeros(4)
    res[1] = FOC_data_buying(l_P, D, l, D_E; par, type)
    res[2] = FOC_data_gen(l_G, l_P, D, l, D_I; par, type)
    res[3] = FOC_data_proc(l_P, D, l; par, type)
    res[4] = FOC_data_sharing(l_P, D, D_I, l; par, type)

    return sum(abs.(res)) # Try to force the algorithm to come up with something more precise
end

function residuals_unc_nl(x, res; par::Pars_v3, type::Symbol, sec_period::Symbol = :no)
    l_P = x[1]^2
    l_G = x[2]^2
    D_E = x[3]^2
    D_I = x[4]^2

    D = bundle(D_I, D_E; par, type)
    l = labor_demand(l_P, D; par, type)

    #! FOCs
    res[1] = FOC_data_buying(l_P, D, l, D_E; par, type, sec_period)
    res[2] = FOC_data_gen(l_G, l_P, D, l, D_I; par, type, sec_period)
    res[3] = FOC_data_proc(l_P, D, l; par, type, sec_period)
    res[4] = FOC_data_sharing(l_P, D, D_I, l; par, type, sec_period)

    return res # Try to force the algorithm to come up with something more precise
end

function residuals_unc_nl2(x, res; par::Pars_v3, type::Symbol, sec_period::Symbol = :no)
    #! FOCs
    if sec_period == :yes
        l_P = x[1]^2
        l_G = x[2]^2
        D_E = x[3]^2
        D_I = x[4]^2

        D = bundle(D_I, D_E; par, type)
        l = labor_demand(l_P, D; par, type, sec_period)
        res[1] = FOC_data_buying(l_P, D, l, D_E; par, type, sec_period)
        res[2] = FOC_data_gen(l_G, l_P, D, l, D_I; par, type, sec_period)
        res[3] = FOC_data_proc(l_P, D, l; par, type, sec_period)
        res[4] = FOC_data_sharing(l_P, D, D_I, l; par, type, sec_period)
    elseif sec_period == :no
        l_P = x[1]^2
        l_G = x[2]^2
        D_E = x[3]^2
        D_I = x[4]^2

        l_P2 = x[5]^2
        l_G2 = x[6]^2
        D_E2 = x[7]^2
        D_I2 = x[8]^2
        D = bundle(D_I, D_E; par, type)
        D2 = bundle(D_I2, D_E2; par, type)
        l = labor_demand(l_P, D; par, type)
        l2 = labor_demand(l_P2, D2; par, type)
        #! This is the solution to the decision in the second period
        res[1] = FOC_data_buying(l_P2, D2, l2, D_E2; par, type, sec_period=:yes)
        res[2] = FOC_data_gen(l_G2, l_P2, D2, l2, D_I2; par, type, sec_period=:yes)
        res[3] = FOC_data_proc(l_P2, D2, l2; par, type, sec_period=:yes)
        res[4] = FOC_data_sharing(l_P2, D2, D_I2, l2; par, type, sec_period=:yes)

        #! This is the solution to the decision in the first period
        res[5] = FOC_data_buying2(l_P, D, l, D_E, l_P2, D2, l2; par, type)
        res[6] = FOC_data_gen2(l_G, l_P, D, l, D_I, l_P2, D2, l2; par, type)
        res[7] = FOC_data_proc2(l_P, D, l, l_P2, D2, l2; par, type)
        res[8] = FOC_data_sharing2(l_P, D, D_I, l, l_P2, D2, l2; par, type)
    end
    return res # Try to force the algorithm to come up with something more precise
end

"""
    Constrained residuals, D_S = 0, no data selling.
"""
function residuals_con(x; par::Pars_v3, type::Symbol)
    l_P = x[1]^2
    l_G = x[2]^2
    D_E = x[3]^2

    D_I = data_gen(l_G; par, type)
    D = bundle(D_I, D_E; par, type)
    l = labor_demand(l_P, D; par, type)

    #! FOCs
    res = zeros(3)
    res[1] = FOC_data_buying(l_P, D, l, D_E; par, type)
    res[2] = FOC_data_gen(l_G, l_P, D, l, D_I; par, type)
    res[3] = FOC_data_proc(l_P, D, l; par, type)

    return sum(abs.(res)) # Try to force the algorithm to come up with something more precise
end

"""
    Constrained residuals, D_S = 0, no data selling.
"""
function residuals_con_nl(x, res; par::Pars_v3, type::Symbol)
    l_P = x[1]^2
    l_G = x[2]^2
    D_E = x[3]^2

    D_I = data_gen(l_G; par, type)
    D = bundle(D_I, D_E; par, type)
    l = labor_demand(l_P, D; par, type)

    #! FOCs
    res[1] = FOC_data_buying(l_P, D, l, D_E; par, type)
    res[2] = FOC_data_gen(l_G, l_P, D, l, D_I; par, type)
    res[3] = FOC_data_proc(l_P, D, l; par, type)

    return res
end

"""
    Main residuals for finding the general equilibrium steady state
"""
function residuals_SS(x; par::Pars_v3)
    w_guess = x[1]^2
    w_P_guess = x[2]^2
    w_G_guess = x[3]^2

    @set! par.w = w_guess
    @set! par.w_P = w_P_guess
    @set! par.w_G = w_G_guess
    @set! par.p_D = find_p_D(par)

    l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = sol_f_problem(par; type = :BS)
    l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = sol_f_problem(par; type = :LS)
    l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = sol_f_problem(par; type = :LU)
    l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = sol_f_problem(par; type = :BU)

    Π_BS = firm_revenue(l_P_BS, D_BS, l_BS; par, type = :BS)
    Π_LS = firm_revenue(l_P_LS, D_LS, l_LS; par, type = :LS)
    Π_LU = firm_revenue(l_P_LU, D_LU, l_LU; par, type = :LU)
    Π_BU = firm_revenue(l_P_BU, D_BU, l_BU; par, type = :BU)

    D = [D_BS, D_LS, D_LU, D_BU]
    D_I = [D_I_BS, D_I_LS, D_I_LU, D_I_BU]
    Π = [Π_BS, Π_LS, Π_LU, Π_BU]
    l_G = [l_G_BS, l_G_LS, l_G_LU, l_G_BU]

    w = wage_int_goods(Π; par)
    w_G = wage_gen_data(Π, D, l_G, D_I; par)
    w_P = wage_gen_proc(Π, D; par)

    res = zeros(3)

    res[1] = w_guess - w
    res[2] = w_P_guess - w_P
    res[3] = w_G_guess - w_G

    return sum(res .^ 2)
end

"""
    Including a guess on Y.
"""
function residuals_SS2(x; par::Pars_v3)
    Y_guess = x[1]^2
    w_guess = x[2]^2
    w_P_guess = x[3]^2
    w_G_guess = x[4]^2
    #p_D_guess = x[5]^2

    @set! par.Y = Y_guess
    @set! par.w = w_guess
    @set! par.w_P = w_P_guess
    @set! par.w_G = w_G_guess
    @set! par.p_D = find_p_D(par)

    sol_BS = solve_firm_problem(par; type = :BS)
    sol_LS = solve_firm_problem(par; type = :LS)
    sol_LU = solve_firm_problem(par; type = :LU)
    sol_BU = solve_firm_problem(par; type = :BU)

    l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = get_solution(sol_BS; par, type = :BS)
    l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = get_solution(sol_LS; par, type = :LS)
    l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = get_solution(sol_LU; par, type = :LU)
    l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = get_solution(sol_BU; par, type = :BU)

    #D_BS = bundle(D_I_BS, D_E_BS; par)
    #D_LS = bundle(D_I_LS, D_E_LS; par)
    #D_LU = bundle(D_I_LU, D_E_LU; par)
    #D_BU = bundle(D_I_BU, D_E_BU; par)

    Y_BS = firm_output(l_P_BS, D_BS, l_BS; par, type = :BS)
    Y_LS = firm_output(l_P_LS, D_LS, l_LS; par, type = :LS)
    Y_LU = firm_output(l_P_LU, D_LU, l_LU; par, type = :LU)
    Y_BU = firm_output(l_P_BU, D_BU, l_BU; par, type = :BU)

    Y = aggregate_output(Y_BS, Y_LS, Y_LU, Y_BU; par)

    res = zeros(4)

    res[1] = (Y - Y_guess)
    res[2] = labor_market_clearing(l_BS, l_LS, l_LU, l_BU; par)
    res[3] = labor_market_clearing(l_P_BS, l_P_LS, l_P_LU, l_P_BU; par)
    res[4] = labor_market_clearing(l_G_BS, l_G_LS, l_G_LU, l_G_BU; par)

    # l_G = [l_G_BS, l_G_LS, l_G_LU, l_G_BU]
    # D_E = [D_E_BS, D_E_LS, D_E_LU, D_E_BU]
    # D_I = [D_I_BS, D_I_LS, D_I_LU, D_I_BU]

    # res[5] = data_market_clearing(l_G, D_E, D_I; par)

    return sum(res .^ 2)
end

"""
    Data market clearing condition
"""
function data_market_clearing(l_G, D_E, D_I; par::Pars_v3)
    l_G_BS = l_G[1]
    l_G_LS = l_G[2]
    l_G_LU = l_G[3]
    l_G_BU = l_G[4]

    D_E_BS = D_E[1]
    D_E_LS = D_E[2]
    D_E_LU = D_E[3]
    D_E_BU = D_E[4]

    D_I_BS = D_I[1]
    D_I_LS = D_I[2]
    D_I_LU = D_I[3]
    D_I_BU = D_I[4]

    D_S_BS = sold_data(D_I_BS, l_G_BS; par, type = :BS)
    D_S_LS = sold_data(D_I_LS, l_G_LS; par, type = :LS)
    D_S_LU = sold_data(D_I_LU, l_G_LU; par, type = :LU)
    D_S_BU = sold_data(D_I_BU, l_G_BU; par, type = :BU)

    w_BS = weight_of_type(:BS; par)
    w_LS = weight_of_type(:LS; par)
    w_LU = weight_of_type(:LU; par)
    w_BU = weight_of_type(:BU; par)

    D_S = w_BS * D_S_BS + w_LS * D_S_LS + w_LU * D_S_LU + w_BU * D_S_BU
    D_B = w_BS * D_E_BS + w_LS * D_E_LS + w_LU * D_E_LU + w_BU * D_E_BU
    #! 1-τ times all data sold is equal to what is bought
    return D_S * (1 - par.τ) - D_B
end

function residuals_p_D(p_D; par::Pars_v3)

    @set! par.p_D = p_D

    l_P_BS, l_G_BS, l_BS, D_E_BS, D_I_BS, D_BS = sol_f_problem_nl(par, type = :BS)
    l_P_LS, l_G_LS, l_LS, D_E_LS, D_I_LS, D_LS = sol_f_problem_nl(par, type = :LS)
    l_P_LU, l_G_LU, l_LU, D_E_LU, D_I_LU, D_LU = sol_f_problem_nl(par, type = :LU)
    l_P_BU, l_G_BU, l_BU, D_E_BU, D_I_BU, D_BU = sol_f_problem_nl(par, type = :BU)

    l_G = [l_G_BS, l_G_LS, l_G_LU, l_G_BU]
    D_E = [D_E_BS, D_E_LS, D_E_LU, D_E_BU]
    D_I = [D_I_BS, D_I_LS, D_I_LU, D_I_BU]

    data_market_clearing(l_G, D_E, D_I; par)
end


text:
#=
is_unique_bit_vector:
- Julia version: 1.5.0
- Author: shisa
- Date: 2020-09-23
=#

module UniqueCharBitChecker
    # methods
    function check_unique(str)
        println("")
        println("Check string $(str) is unique or not")

        checker = 0

        for i in 1:length(str)
            val = str[i] - 'a'
            if (checker & (1 << val)) > 0
                println("Not unique string")
                return false
            end
            checker |= (1 << val)
        end

        println("Unique string")
        return true
    end

    function main()
        @time println(check_unique("abcdefg"))
        @time println(check_unique("aaaaaaaaaaaa"))
        @time println(check_unique("axgilnrewqczpoh"))
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    using .UniqueCharBitChecker
    UniqueCharChecker.main()
end

max_stars_repo_name: ShisatoYano/JuliaPractice

max_issues_repo_name: ShisatoYano/JuliaPractice

max_forks_repo_name: ShisatoYano/JuliaPractice


text:
export Product, LibraryProduct, FileProduct, ExecutableProduct, satisfied,
       locate, write_deps_file, variable_name
import Base: repr

"""
A `Product` is an expected result after building or installation of a package.

Examples of `Product`s include `LibraryProduct`, `ExecutableProduct` and
`FileProduct`.  All `Product` types must define the following minimum set of
functionality:

* `locate(::Product)`: given a `Product`, locate it within the wrapped `Prefix`
  returning its location as a string

* `satisfied(::Product)`: given a `Product`, determine whether it has been
  successfully satisfied (e.g. it is locateable and it passes all callbacks)

* `variable_name(::Product)`: return the variable name assigned to a `Product`

* `repr(::Product)`: Return a representation of this `Product`, useful for
  auto-generating source code that constructs `Products`, if that's your thing.
"""
abstract type Product end

"""
    satisfied(p::Product; platform::Platform = platform_key_abi(),
              verbose::Bool = false, isolate::Bool = false)

Given a `Product`, return `true` if that `Product` is satisfied, e.g. whether
a file exists that matches all criteria setup for that `Product`.  If `isolate`
is set to `true`, will isolate all checks from the main Julia process in the
event that `dlopen()`'ing a library might cause issues.
"""
function satisfied(p::Product; platform::Platform = platform_key_abi(),
                               verbose::Bool = false, isolate::Bool = false)
    return locate(p; platform=platform, verbose=verbose, isolate=isolate) != nothing
end

"""
    variable_name(p::Product)

Return the variable name associated with this `Product` as a string
"""
function variable_name(p::Product)
    return string(p.variable_name)
end

"""
A `LibraryProduct` is a special kind of `Product` that not only needs to exist,
but needs to be `dlopen()`'able.  You must know which directory the library
will be installed to, and its name, e.g. to build a `LibraryProduct` that
refers to `"/lib/libnettle.so"`, the "directory" would be "/lib", and the
"libname" would be "libnettle".  Note that a `LibraryProduct` can support
multiple libnames, as some software projects change the libname based on the
build configuration.
"""
struct LibraryProduct <: Product
    dir_path::Union{String, Nothing}
    libnames::Vector{String}
    variable_name::Symbol
    prefix::Union{Prefix, Nothing}

    """
        LibraryProduct(prefix::Prefix, libname::AbstractString,
                       varname::Symbol)

    Declares a `LibraryProduct` that points to a library located within the
    `libdir` of the given `Prefix`, with a name containing `libname`.  As an
    example, given that `libdir(prefix)` is equal to `usr/lib`, and `libname`
    is equal to `libnettle`, this would be satisfied by the following paths:

        usr/lib/libnettle.so
        usr/lib/libnettle.so.6
        usr/lib/libnettle.6.dylib
        usr/lib/libnettle-6.dll

    Libraries matching the search pattern are rejected if they are not
    `dlopen()`'able.
    """
    function LibraryProduct(prefix::Prefix, libname::AbstractString,
                            varname::Symbol)
        return LibraryProduct(prefix, [libname], varname)
    end

    function LibraryProduct(prefix::Prefix, libnames::Vector{S},
                            varname::Symbol) where {S <: AbstractString}
        @show prefix, libnames, varname
        return new(nothing, libnames, varname, prefix)
    end

    """
        LibraryProduct(dir_path::AbstractString, libname::AbstractString,
                       varname::Symbol)

    For finer-grained control over `LibraryProduct` locations, you may directly
    pass in the `dir_path` instead of auto-inferring it from `libdir(prefix)`.
    """
    function LibraryProduct(dir_path::AbstractString, libname::AbstractString,
                            varname::Symbol)
        return LibraryProduct(dir_path, [libname], varname)
    end

    function LibraryProduct(dir_path::AbstractString, libnames::Vector{S},
                            varname::Symbol) where {S <: AbstractString}
        @show dir_path, libnames, varname
       return new(dir_path, libnames, varname, nothing)
    end
end

function repr(p::LibraryProduct)
    libnames = repr(p.libnames)
    varname = repr(p.variable_name)
    if p.prefix === nothing
        return "LibraryProduct($(repr(p.dir_path)), $(libnames), $(varname))"
    else
        return "LibraryProduct(prefix, $(libnames), $(varname))"
    end
end

"""
locate(lp::LibraryProduct; verbose::Bool = false,
        platform::Platform = platform_key_abi())

If the given library exists (under any reasonable name) and is `dlopen()`able,
(assuming it was built for the current platform) return its location.  Note
that the `dlopen()` test is only run if the current platform matches the given
`platform` keyword argument, as cross-compiled libraries cannot be `dlopen()`ed
on foreign platforms.
"""
function locate(lp::LibraryProduct; verbose::Bool = false,
                platform::Platform = platform_key_abi(), isolate::Bool = false)
    dir_path = lp.dir_path
    if dir_path === nothing
        dir_path = libdir(lp.prefix, platform)
    end

    if !isdir(dir_path)
        if verbose
            @info("Directory $(dir_path) does not exist!")
        end
        return nothing
    end
    for f in readdir(dir_path)
        # Skip any names that aren't a valid dynamic library for the given
        # platform (note this will cause problems if something compiles a `.so`
        # on OSX, for instance)
        if !valid_dl_path(f, platform)
            continue
        end

        if verbose
            @info("Found a valid dl path $(f) while looking for $(join(lp.libnames, ", "))")
        end

        # If we found something that is a dynamic library, let's check to see
        # if it matches our libname:
        for libname in lp.libnames
            if startswith(basename(f), libname)
                dl_path = abspath(joinpath(dir_path), f)
                if verbose
                    @info("$(dl_path) matches our search criteria of $(libname)")
                end

                # If it does, try to `dlopen()` it if the current platform is good
                if platforms_match(platform, platform_key_abi())
                    if isolate
                        # Isolated dlopen is a lot slower, but safer
                        dl_esc_path = replace(dl_path, "\\"=>"\\\\")
                        if success(`$(Base.julia_cmd()) -e "import Libdl; Libdl.dlopen(\"$(dl_esc_path)\")"`)
                            return dl_path
                        end
                    else
                        hdl = Libdl.dlopen_e(dl_path)
                        if !(hdl in (C_NULL, nothing))
                            Libdl.dlclose(hdl)
                            return dl_path
                        end
                    end

                    if verbose
                        try
                            dlopen(dl_path)
                        catch dlopen_result
                            @info("$(dl_path) cannot be dlopen'ed",dlopen_result)
                        end
                    end
                else
                    # If the current platform doesn't match, then just trust in our
                    # cross-compilers and go with the flow
                    return dl_path
                end
            end
        end
    end

    if verbose
        @info("Could not locate $(join(lp.libnames, ", ")) inside $(dir_path)")
    end
    return nothing
end

"""
An `ExecutableProduct` is a `Product` that represents an executable file.

On all platforms, an ExecutableProduct checks for existence of the file.  On
non-Windows platforms, it will check for the executable bit being set.  On
Windows platforms, it will check that the file ends with ".exe", (adding it on
automatically, if it is not already present).
"""
struct ExecutableProduct <: Product
    path::AbstractString
    variable_name::Symbol
    prefix::Union{Prefix, Nothing}

    """
    `ExecutableProduct(prefix::Prefix, binname::AbstractString,
                       varname::Symbol)`

    Declares an `ExecutableProduct` that points to an executable located within
    the `bindir` of the given `Prefix`, named `binname`.
    """
    function ExecutableProduct(prefix::Prefix, binname::AbstractString,
                               varname::Symbol)
        return new(joinpath(bindir(prefix), binname), varname, prefix)
    end

    """
    `ExecutableProduct(binpath::AbstractString, varname::Symbol)`

    For finer-grained control over `ExecutableProduct` locations, you may directly
    pass in the full `binpath` instead of auto-inferring it from `bindir(prefix)`.
    """
    function ExecutableProduct(binpath::AbstractString, varname::Symbol)
        return new(binpath, varname, nothing)
    end
end

function repr(p::ExecutableProduct)
    varname = repr(p.variable_name)
    if p.prefix === nothing
        return "ExecutableProduct($(repr(p.path)), $(varname))"
    else
        rp = relpath(p.path, bindir(p.prefix))
        return "ExecutableProduct(prefix, $(repr(rp)), $(varname))"
    end
end

"""
`locate(fp::ExecutableProduct; platform::Platform = platform_key_abi(),
                               verbose::Bool = false, isolate::Bool = false)`

If the given executable file exists and is executable, return its path.

On all platforms, an ExecutableProduct checks for existence of the file.  On
non-Windows platforms, it will check for the executable bit being set.  On
Windows platforms, it will check that the file ends with ".exe", (adding it on
automatically, if it is not already present).
"""
function locate(ep::ExecutableProduct; platform::Platform = platform_key_abi(),
                verbose::Bool = false, isolate::Bool = false)
    # On windows, we always slap an .exe onto the end if it doesn't already
    # exist, as Windows won't execute files that don't have a .exe at the end.
    path = if platform isa Windows && !endswith(ep.path, ".exe")
        "$(ep.path).exe"
    else
        ep.path
    end

    if !isfile(path)
        if verbose
            @info("$(ep.path) does not exist, reporting unsatisfied")
        end
        return nothing
    end

    # If the file is not executable, fail out (unless we're on windows since
    # windows doesn't honor these permissions on its filesystems)
    @static if !Sys.iswindows()
        if uperm(path) & 0x1 == 0
            if verbose
                @info("$(path) is not executable, reporting unsatisfied")
            end
            return nothing
        end
    end

    return path
end

"""
A `FileProduct` represents a file that simply must exist to be satisfied.
"""
struct FileProduct <: Product
    path::AbstractString
    variable_name::Symbol
    prefix::Union{Prefix, Nothing}

    """
        FileProduct(prefix::Prefix, relative_path::AbstractString,
                                    varname::Symbol)`

    Declares a `FileProduct` that points to a file located relative to a the
    root of a `Prefix`.
    """
    function FileProduct(prefix::Prefix, relative_path::AbstractString,
                                         varname::Symbol)
        file_path = joinpath(prefix.path, relative_path)
        return new(file_path, varname, prefix)
    end

    """
        FileProduct(file_path::AbstractString, varname::Symbol)

    For finer-grained control over `FileProduct` locations, you may directly
    pass in the full `file_pathpath` instead of defining it in reference to
    a root `Prefix`.
    """
    function FileProduct(file_path::AbstractString, varname::Symbol)
        return new(file_path, varname, nothing)
    end
end

function repr(p::FileProduct)
    varname = repr(p.variable_name)
    if p.prefix === nothing
        return "FileProduct($(repr(p.path)), $(varname))"
    else
        rp = relpath(p.path, p.prefix.path)
        return "FileProduct(prefix, $(repr(rp)), $(varname))"
    end
end

"""
locate(fp::FileProduct; platform::Platform = platform_key_abi(),
                        verbose::Bool = false, isolate::Bool = false)

If the given file exists, return its path.  The platform argument is ignored
here, but included for uniformity.
"""
function locate(fp::FileProduct; platform::Platform = platform_key_abi(),
                                 verbose::Bool = false, isolate::Bool = false)
    # Limited variable expansion capabilities
    mappings = Dict()

    for (var, val) in [("target", triplet(platform)), ("nbits", wordsize(platform))]
        mappings["\$$(var)"] = string(val)
        mappings["\${$(var)}"] = string(val)
    end

    expanded = fp.path
    for (old, new) in mappings
        expanded = replace(expanded, old => new)
    end

    if isfile(expanded)
        if verbose
            @info("FileProduct $(fp.path) found at $(realpath(expanded))")
        end
        return expanded
    end
    if verbose
        @info("FileProduct $(fp.path) not found")
    end
    return nothing
end

"""
    write_deps_file(depsjl_path::AbstractString, products::Vector{Product};
                    verbose::Bool = false)

Generate a `deps.jl` file that contains the variables referred to by the
products within `products`.  As an example, running the following code:

    fooifier = ExecutableProduct(..., :foo_exe)
    libbar = LibraryProduct(..., :libbar)
    write_deps_file(joinpath(@__DIR__, "deps.jl"), [fooifier, libbar])

Will generate a `deps.jl` file that contains definitions for the two variables
`foo_exe` and `libbar`.  If any `Product` object cannot be satisfied (e.g.
`LibraryProduct` objects must be `dlopen()`-able, `FileProduct` objects must
exist on the filesystem, etc...) this method will error out.  Ensure that you
have used `install()` to install the binaries you wish to write a `deps.jl`
file for.

The result of this method is a `deps.jl` file containing variables named as
defined within the `Product` objects passed in to it, holding the full path to the
installed binaries.  Given the example above, it would contain code similar to:

    global const foo_exe = "<pkg path>/deps/usr/bin/fooifier"
    global const libbar = "<pkg path>/deps/usr/lib/libbar.so"

This `deps.jl` file is intended to be `include()`'ed from within the top-level
source of your package.  Note that all files are checked for consistency on
package load time, and if an error is discovered, package loading will fail,
asking the user to re-run `Pkg.build("package_name")`.
"""
function write_deps_file(depsjl_path::AbstractString, products::Vector{P};
                         verbose::Bool=false, isolate::Bool=true) where {P <: Product}
    # helper function to escape paths
    escape_path = path -> replace(path, "\\" => "\\\\")

    # Grab the package name as the name of the top-level directory of a package
    # Auto-detect being placed within the `~/.julia/packages/XXXXX` directory
    if basename(dirname(dirname(dirname(dirname(depsjl_path))))) == "packages"
        package_name = basename(dirname(dirname(dirname(depsjl_path))))
    else
        # otherwise assume same as folder name
        # TODO: check Project.toml
        package_name = basename(dirname(dirname(depsjl_path)))
    end
    # We say this a couple of times
    rebuild = strip("""
    Please re-run Pkg.build(\\\"$(package_name)\\\"), and restart Julia.
    """)

    # Begin by ensuring that we can satisfy every product RIGHT NOW
    for p in products
        if !satisfied(p; verbose=verbose, isolate=isolate)
            error("$p is not satisfied, cannot generate deps.jl!")
        end
    end

    # If things look good, let's generate the `deps.jl` file
    open(depsjl_path, "w") do depsjl_file
        # First, dump the preamble
        println(depsjl_file, strip("""
        ## This file autogenerated by BinaryProvider.write_deps_file().
        ## Do not edit.
        ##
        ## Include this file within your main top-level source, and call
        ## `check_deps()` from within your module's `__init__()` method

        if isdefined((@static VERSION < v"0.7.0-DEV.484" ? current_module() : @__MODULE__), :Compat)
            import Compat.Libdl
        elseif VERSION >= v"0.7.0-DEV.3382"
            import Libdl
        end
        """))

        # Next, spit out the paths of all our products
        for product in products
            # Escape the location so that e.g. Windows platforms are happy with
            # the backslashes in a string literal
            product_path = locate(product, platform=platform_key_abi(),
                                           verbose=verbose)
            product_path = relpath(product_path, dirname(depsjl_path))
            product_path = escape_path(product_path)
            vp = variable_name(product)
            println(depsjl_file, strip("""
            const $(vp) = joinpath(dirname(@__FILE__), \"$(product_path)\")
            """))
        end

        # Next, generate a function to check they're all on the up-and-up
        println(depsjl_file, "function check_deps()")

        for product in products
            varname = variable_name(product)

            # Add a `global $(name)`
            println(depsjl_file, "    global $(varname)");

            # Check that any file exists
            println(depsjl_file, """
                if !isfile($(varname))
                    error("\$($(varname)) does not exist, $(rebuild)")
                end
            """)

            # For Library products, check that we can dlopen it:
            if typeof(product) <: LibraryProduct
                println(depsjl_file, """
                    if Libdl.dlopen_e($(varname)) in (C_NULL, nothing)
                        error("\$($(varname)) cannot be opened, $(rebuild)")
                    end
                """)
            end
        end

        # If any of the products are `ExecutableProduct`s, we need to add Julia's
        # library directory onto the end of {DYLD,LD}_LIBRARY_PATH
        @static if !Sys.iswindows()
            if any(p isa ExecutableProduct for p in products)
                dllist = Libdl.dllist()
                libjulia = filter(x -> occursin("libjulia", x), dllist)[1]
                julia_libdir = repr(joinpath(dirname(libjulia), "julia"))
                envvar_name = @static if Sys.isapple()
                    "DYLD_LIBRARY_PATH"
                else Sys.islinux()
                    "LD_LIBRARY_PATH"
                end
                envvar_name = repr(envvar_name)

                println(depsjl_file, """
                    libpaths = split(get(ENV, $(envvar_name), ""), ":")
                    if !($(julia_libdir) in libpaths)
                        push!(libpaths, $(julia_libdir))
                    end
                    ENV[$(envvar_name)] = join(filter(!isempty, libpaths), ":")
                """)
            end
        end

        # Close the `check_deps()` function
        println(depsjl_file, "end")
    end
end


text:
import Cassette
Cassette.@context OmegaCtx

# if there's already a context append to it
# @inline (x::Intervened)(ω) = Cassette.overdub(OmegaCtx(metadata = (intervention = x.i,)), x.x, ω)

# function Cassette.overdub(ctx::OmegaCtx{NamedTuple{K, V}}, x, ω::MaybeTagged{<:AbstractΩ}) where {K, V}
#   # @show K
#   # if :intervention in K
#   #   :(
#       handleintervention(ctx, ctx.metadata.intervention, x, ω)
#       # )
#   # else
#   #   :(Cassette.recurse(ctx, x, ω))
#   # end
# end

# @inline function handleintervention(ctx, i::Intervention{X, V}, x::X, ω) where {X, V}
#   if i.x == x
#     i.v
#   else
#     Casstte.recurse(ctx, x, ω)
#   end
# end

# @inline handleintervention(ctx, ::Intervention, x, ω) = Cassette.recurse(ctx, x, ω)

# # # Syntactic Sugar