text:
### ARS holds the fields of the acceptance-rejection sampler

immutable ARS <: MCSampler
  logproposal::Function # Possibly unnormalized log-proposal
  proposalscale::Float64 # Scale factor to ensure the scaled-up logproposal covers target
  jumpscale::Float64 # Scale factor for adapting the jump size

  function ARS(l::Function, ps::Float64, js::Float64)
    @assert typeof(l) == Function "logproposal should be a function."
    @assert js > 0 "Scale factor for adapting the jump size is not positive."
    new(l, ps, js)
  end
end

ARS(l::Function) = ARS(l, 1., 1.)
ARS(l::Function, ps::Float64) = ARS(l, ps, 1.)

ARS(l::Function; proposalscale::Float64 = 1.0, jumpscale::Float64=1.0) = ARS(l, proposalscale, jumpscale)

### ARSHeap type holds the internal state ("local variables") of the ARS sampler

type ARSHeap <: MCHeap{MCBaseSample}
  instate::MCState{MCBaseSample} # Monte Carlo state used internally by the sampler
  outstate::MCState{MCBaseSample} # Monte Carlo state outputted by the sampler
  tune::MCTune
  count::Int # Current number of iterations
  jumpscale::Vector{Float64}
  logproposal::Float64
  weight::Float64
end

ARSHeap() =
  ARSHeap(MCState(MCBaseSample(), MCBaseSample()), MCState(MCBaseSample(), MCBaseSample()), VanillaMCTune(), 0,
  Float64[], NaN, NaN)

ARSHeap(l::Int) =
  ARSHeap(MCState(MCBaseSample(l), MCBaseSample(l)), MCState(MCBaseSample(l), MCBaseSample(l)), VanillaMCTune(), 0,
  Float64[], NaN, NaN)

### Initialize Metropolis-Hastings sampler

function initialize_heap(m::MCModel, s::ARS, r::MCRunner, t::MCTuner)
  heap::ARSHeap = ARSHeap(m.size)

  heap.instate.current = MCBaseSample(copy(m.init))
  logtarget!(heap.instate.current, m.eval)
  @assert isfinite(heap.instate.current.logtarget) "Initial values out of model support."

  heap.jumpscale = m.scale.*s.jumpscale  # rescale model scale by sampler scale

  heap.tune = VanillaMCTune()

  heap.count = 1

  heap
end

function reset!(heap::ARSHeap, x::Vector{Float64}, m::MCModel)
  heap.instate.current = MCBaseSample(copy(x))
  logtarget!(heap.instate.current, m.eval)
end

function initialize_task!(heap::ARSHeap, m::MCModel, s::ARS, r::MCRunner, t::MCTuner)
  # Hook inside Task to allow remote resetting
  task_local_storage(:reset, (x::Vector{Float64})->reset!(heap, x, m))

  while true
    iterate!(heap, m, s, r, t, produce)
  end
end

### Perform iteration for ARS sampler

function iterate!(heap::ARSHeap, m::MCModel, s::ARS, r::MCRunner, t::MCTuner, send::Function)
  if t.verbose
    heap.tune.proposed += 1
  end

  heap.instate.successive = MCBaseSample(heap.instate.current.sample+randn(m.size).*heap.jumpscale)
  logtarget!(heap.instate.successive, m.eval)
  heap.logproposal = s.logproposal(heap.instate.successive.sample)

  heap.weight = heap.instate.successive.logtarget.-s.proposalscale.-heap.logproposal
  if heap.weight > log(rand())
    heap.outstate = MCState(heap.instate.successive, heap.instate.current, Dict{Any, Any}("accept" => true))
    heap.instate.current = deepcopy(heap.instate.successive)

    if t.verbose
      heap.tune.accepted += 1
    end
  else
    heap.outstate = MCState(heap.instate.current, heap.instate.current, Dict{Any, Any}("accept" => false))
  end

  if t.verbose && heap.count <= r.burnin && mod(heap.count, t.period) == 0
    rate!(heap.tune)
    println("Burnin iteration $(heap.count) of $(r.burnin): ", round(100*heap.tune.rate, 2), " % acceptance rate")
  end

  heap.count += 1

  send(heap.outstate)
end


text:
import RowEchelon
import LinearAlgebra

#=
Homework 7
Section 3.2
Problems 29
=#

# Create B
B = [1 0 1 ; 1 1 2 ; 1 2 1]
# Create B^5
B5 = B^5
# Get detB5
detB5 = LinearAlgebra.det(B5)

println("Original Matrix, B:")
println(B)
println("B^5:")
println(B5)
println("Determinate of B^5")
println(detB5)

ext: jl


text:
# ---
# title: 1108. Defanging an IP Address
# id: problem1108
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: String
# link: <https://leetcode.com/problems/defanging-an-ip-address/description/>
# hidden: true
# ---
#
# Given a valid (IPv4) IP `address`, return a defanged version of that IP
# address.
#
# A _defanged  IP address_ replaces every period `"."` with `"[.]"`.
#
#
#
# **Example 1:**
#
#
#
#     Input: address = "1.1.1.1"
#     Output: "1[.]1[.]1[.]1"
#
#
# **Example 2:**
#
#
#
#     Input: address = "255.100.50.0"
#     Output: "255[.]100[.]50[.]0"
#
#
#
#
# **Constraints:**
#
#   * The given `address` is a valid IPv4 address.
#
#
## @lc code=start
using LeetCode

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


text:
module Utils

#include("SparseArrays.jl")
include("BoundingBoxes.jl")
include("FakeSegmentations.jl")
include("SynapseTables.jl")
include("Mathes.jl")
include("VectorClouds.jl")
include("RangeIndexingArrays.jl")

end # module


text:
# Note (1). See
# https://discourse.julialang.org/t/weight-regularisation-which-iterates-params-m-in-flux-mutating-arrays-is-not-supported/64314

""" Penalizer(λ, α)

Returns a callable object `penalizer` for evaluating regularization
penalties associated with some numerical array. Specifically,
`penalizer(A)` returns

   λ*(α*L1 + (1 - α)*L2),

where `L1` is the sum of absolute values of the elments of `A` and
`L2` is the sum of squares of those elements.
"""
struct Penalizer{T}
    lambda::T
    alpha::T
    function Penalizer(lambda, alpha)
        lambda == 0 && return new{Nothing}(nothing, nothing)
        T = promote_type(typeof.((lambda, alpha))...)
        return new{T}(lambda, alpha)
    end
end

(::Penalizer{Nothing})(::Any) = 0
function (p::Penalizer)(A)
    λ = p.lambda
    α = p.alpha
    # avoiding broadcasting; see Note (1) above
    L2 = sum(abs2, A)
    L1 = sum(abs,  A)
    return  λ*(α*L1 + (1 - α)*L2)
end

"""
    Penalty(model)

Returns a callable object `p`, for returning the regularization
penalty `p(w)` associated with some collection of parameters `w`. For
example, `w = params(chain)` where `chain` is some Flux
model. Here `model` is an MLJFlux model ("model" in the MLJ
sense, not the Flux sense). Specifically, `p(w)` returns

   sum(Penalizer(λ, α).w)

where `α = model.alpha`, `λ = model.lambda`.

See also [`Penalizer`](@ref)

"""
struct Penalty{P}
    penalizer::P
    function Penalty(model)
        penalizer = Penalizer(model.lambda, model.alpha)
        return new{typeof(penalizer)}(penalizer)
    end
end
(p::Penalty{Penalizer{Nothing}})(w) = zero(Float32)
(p::Penalty)(w) = sum(p.penalizer(wt) for wt in w)


text:
# ---
# title: 1456. Maximum Number of Vowels in a Substring of Given Length
# id: problem1456
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: String, Sliding Window
# link: <https://leetcode.com/problems/maximum-number-of-vowels-in-a-substring-of-given-length/description/>
# hidden: true
# ---
#
# Given a string `s` and an integer `k`.
#
# Return _the maximum number of vowel letters_ in any substring of `s` with
# length `k`.
#
# **Vowel letters** in English are (a, e, i, o, u).
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abciiidef", k = 3
#     Output: 3
#     Explanation: The substring "iii" contains 3 vowel letters.
#
#
# **Example 2:**
#
#
#
#     Input: s = "aeiou", k = 2
#     Output: 2
#     Explanation: Any substring of length 2 contains 2 vowels.
#
#
# **Example 3:**
#
#
#
#     Input: s = "leetcode", k = 3
#     Output: 2
#     Explanation: "lee", "eet" and "ode" contain 2 vowels.
#
#
# **Example 4:**
#
#
#
#     Input: s = "rhythms", k = 4
#     Output: 0
#     Explanation: We can see that s doesn't have any vowel letters.
#
#
# **Example 5:**
#
#
#
#     Input: s = "tryhard", k = 4
#     Output: 1
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 10^5`
#   * `s` consists of lowercase English letters.
#   * `1 <= k <= s.length`
#
#
## @lc code=start
using LeetCode

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


text:
using Images, WNNMDenoise, TestImages
using Random

img = imresize(Float64.(load("house.png")) * 255, ratio=0.5);
img = Float32.(load("house.png")) * 255;

clean_img = img;
σₙ = 40;
Random.seed!(0)
noisy_img = clean_img .+ σₙ .* randn(Float32, size(clean_img))

assess_psnr(noisy_img, img, 255)

out = copy(noisy_img)
f = WNNM(σₙ)
@time f(out, noisy_img; clean_img=clean_img)


text:
"""
    SubgradientMaster

Implementation of projected subgradient method
"""
mutable struct SubgradientMaster <: AbstractLagrangeMaster
    num_vars::Int
    num_functions::Int
    eval_f::Union{Nothing,Function}

    iter::Int # current iteration count
    maxiter::Int
    obj_limit::Float64

    f::Float64
    best_f::Float64
    x::Vector{Float64}
    best_x::Vector{Float64}
    ∇f::Vector{Float64}
    α::Float64 # step size

    step_size::Function

    constraint_matrix::Union{Nothing,SparseMatrixCSC{Float64,Int}}

    iteration_time::Vector{Float64}

    function SubgradientMaster()
        sg = new()
        sg.num_vars = 0
        sg.num_functions = 0
        sg.eval_f = nothing
        sg.iter = 0
        sg.maxiter = 1000
        sg.obj_limit = +Inf
        sg.f = -Inf
        sg.best_f = -Inf
        sg.x = []
        sg.best_x = []
        sg.∇f = []
        sg.α = 0.1
        # sg.step_size = (method) -> method.α / method.iter # square summable but not summable
        sg.step_size = (method) -> method.α / sqrt(method.iter) # Nonsummable diminishing
        sg.constraint_matrix = nothing
        sg.iteration_time = []
        return sg
    end
end

function load!(method::SubgradientMaster, num_coupling_variables::Int, num_blocks::Int, eval_function::Function, init_sol::Vector{Float64})
    method.num_vars = num_coupling_variables
    method.num_functions = num_blocks
    method.eval_f = eval_function
    method.x = init_sol
    method.best_x = init_sol
end

function add_constraints!(LD::AbstractLagrangeDual, method::SubgradientMaster)
    I = Int[]; J = Int[]; V = Float64[]
    row = 1
    for (id, vars) in LD.block_model.variables_by_couple
        for v in vars
            push!(I, row)
            push!(J, index_of_λ(LD, v))
            push!(V, 1.0)
        end
        row += 1
    end
    method.constraint_matrix = sparse(I, J, V)
end

function run!(method::SubgradientMaster)

    # projection matrix
    m, n = size(method.constraint_matrix)
    I = sparse(SparseArrays.I, n, n)
    A = method.constraint_matrix
    proj = transpose(A) * ((A * transpose(A)) \ A)

    total_time = 0.0
    total_master_time = 0.0
    total_stime = time()

    method.∇f = zeros(method.num_vars)

    while method.iter < method.maxiter
        method.iter += 1
        if method.iter % 100 == 1
            @printf("%6s", "Iter")
            @printf("\t%13s", "fbest")
            @printf("\t%13s", "f")
            @printf("\t%13s", "|∇f|")
            @printf("\t%13s", "α")
            # @printf("\t%13s", "res")
            @printf("\t%7s", "sg time")
            @printf("\t%8s", "tot time")
            @printf("\n")
        end
        f, g = method.eval_f(method.x)
        method.f = -sum(f)
        if method.best_f < method.f
            method.best_f = method.f
            copy!(method.best_x, method.x)
        end

        master_stime = time()

        fill!(method.∇f, 0.0)
        for (k, gk) in g
            method.∇f .+= gk
        end

        # update
        α = method.step_size(method)
        method.x .-= α * method.∇f

        # projection
        method.x .-= proj * method.x
        # res = method.constraint_matrix * method.x

        # clock iteration time
        push!(method.iteration_time, time() - master_stime)
        total_master_time += time() - master_stime
        total_time = time() - total_stime

        @printf("%6d", method.iter)
        @printf("\t%+6e", method.best_f)
        @printf("\t%+6e", method.f)
        @printf("\t%+6e", norm(method.∇f))
        @printf("\t%+6e", α)
        # @printf("\t%+6e", norm(res))
        @printf("\t%7.2f", total_master_time)
        @printf("\t%8.2f", total_time)
        @printf("\n")

        if method.obj_limit < method.best_f
            break
        end
    end
end

get_objective(method::SubgradientMaster) = method.best_f
get_solution(method::SubgradientMaster) = method.best_x
get_times(method::SubgradientMaster)::Vector{Float64} = method.iteration_time
function set_obj_limit!(method::SubgradientMaster, val::Float64)
    method.obj_limit = val
end


text:
using Base: sin, cos

using ..Geometry

export
  AffineTransform,
  copy,
  to_identity!, set!, transform!,
  translate!, make_translate!,
  scale!, make_scale!,
  rotate!, set_rotate!, make_rotate!,
  multiply!, invert!,
  transpose!

# A minified affine transform.
# Column major
#     x'   |a c tx|   |x|
#     y' = |b d ty| x |y|                  <=== Post multiply
#     1    |0 0  1|   |1|
# or
# Row major
#                           |a  b   0|
#     |x' y' 1| = |x y 1| x |c  d   0|     <=== Pre multiply
#                           |tx ty  1|
#
mutable struct AffineTransform{T <: AbstractFloat}
    a::T
    b::T
    c::T
    d::T
    tx::T
    ty::T

    function AffineTransform{T}() where {T <: AbstractFloat}
        new(1.0, 0.0, 0.0, 1.0, 0.0, 0.0)
    end

    function AffineTransform{T}(a::T, b::T, c::T, d::T, tx::T, ty::T) where {T <: AbstractFloat}
        new(a, b, c, d, tx, ty)
    end
end

# copy specializations
copy(t::AffineTransform) = AffineTransform{Float64}(t.a, t.b, t.c, t.d, t.tx, t.ty)

# initializations
function to_identity!(aft::AffineTransform)
    aft.a = aft.d = 1.0;
    aft.b = aft.c = aft.tx = aft.ty = 0.0;
end

# transforms
function transform!(aft::AffineTransform{T}, v::Geometry.Point{T}) where {T <: AbstractFloat}
    x = (aft.a * v.x) + (aft.c * v.y) + aft.tx
    y = (aft.b * v.x) + (aft.d * v.y) + aft.ty
    v.x = x
    v.y = y;
end

function transform!(aft::AffineTransform{T}, in::Geometry.Point{T}, out::Geometry.Point{T}) where {T <: AbstractFloat}
    out.x = (aft.a * in.x) + (aft.c * in.y) + aft.tx
    out.y = (aft.b * in.x) + (aft.d * in.y) + aft.ty;
end

function transform!(aft::AffineTransform{T}, x::T, y::T, out::Geometry.Point{T}) where {T <: AbstractFloat}
    out.x = (aft.a * x) + (aft.c * y) + aft.tx
    out.y = (aft.b * x) + (aft.d * y) + aft.ty;
end

function transform!(aft::AffineTransform{T}, vx::T, vy::T) where {T <: AbstractFloat}
    ((aft.a * vx) + (aft.c * vy) + aft.tx,
    (aft.b * vx) + (aft.d * vy) + aft.ty)
end

function transform!(aft::AffineTransform{T}, polygon::Geometry.Polygon{T}) where {T <: AbstractFloat}
    for v in polygon.mesh.vertices
        Geometry.set!(v, (aft.a * vx) + (aft.c * vy) + aft.tx, (aft.b * vx) + (aft.d * vy) + aft.ty)
    end
end

# setters/getters
function set!(aft::AffineTransform{T}, a::T, b::T, c::T, d::T, tx::T, ty::T) where {T <: AbstractFloat}
    aft.a = a
    aft.b = b
    aft.c = c
    aft.d = d
    aft.tx = tx
    aft.ty = ty;
end

function set!(to::AffineTransform{T}, from::AffineTransform{T}) where {T <: AbstractFloat}
    to.a = from.a
    to.b = from.b
    to.c = from.c
    to.d = from.d
    to.tx = from.tx
    to.ty = from.ty;
end

# Concatenate translation
function translate!(aft::AffineTransform{T}, tx::T, ty::T) where {T <: AbstractFloat}
    aft.tx += (aft.a * tx) + (aft.c * ty)
    aft.ty += (aft.b * tx) + (aft.d * ty);
end

function make_translate!(aft::AffineTransform{T}, tx::T, ty::T) where {T <: AbstractFloat}
    aft.a = aft.d = 1.0
    aft.b = aft.c = 0.0
    aft.tx = tx
    aft.ty = ty;
end

function scale!(aft::AffineTransform{T}, sx::T, sy::T) where {T <: AbstractFloat}
    aft.a *= sx
    aft.b *= sx
    aft.c *= sy
    aft.d *= sy
end

function make_scale!(aft::AffineTransform{T}, sx::T, sy::T) where {T <: AbstractFloat}
    aft.a = sx
    aft.d = sy
    aft.b = aft.c = aft.tx = aft.ty = 0.0;
end

# Concatinate a rotation (radians) onto this transform.
#
# Rotation is just a matter of perspective. A CW rotation can be seen as
# CCW depending on what you are talking about rotating. For example,
# if the coordinate system is thought as rotating CCW then objects are
# seen as rotating CW, and that is what the 2x2 matrix below represents.
#
# It is also the frame of reference we use. In this library +Y axis is downward
#     |cos  -sin|   object appears to rotate CW.
#     |sin   cos|
#
# In the matrix below the object appears to rotate CCW.
#     |cos  sin|
#     |-sin cos|
#
#     |a  c|    |cos  -sin|
#     |b  d|  x |sin   cos|
#
function rotate!(aft::AffineTransform{T}, radians::T) where {T <: AbstractFloat}
    rsin = sin(radians)
    rcos = cos(radians)
    a = aft.a
    b = aft.b
    c = aft.c
    d = aft.d

    aft.a = a * rcos + c * rsin
    aft.b = b * rcos + d * rsin
    aft.c = c * rcos - a * rsin
    aft.d = d * rcos - b * rsin;
end

function set_rotate!(aft::AffineTransform{T}, radians::T) where {T <: AbstractFloat}
    rsin = sin(radians)
    rcos = cos(radians)
    aft.a = rcos
    aft.b = rsin
    aft.c = -rsin
    aft.d = rcos;
end

function make_rotate!(aft::AffineTransform{T}, radians::T) where {T <: AbstractFloat}
    set_rotate!(aft, radians)
    aft.tx = 0.0
    aft.ty = 0.0;
end

# n = n * m
function multiply_pre!(m::AffineTransform{T}, n::AffineTransform{T}) where {T <: AbstractFloat}
    a = n.a
    b = n.b
    c = n.c
    d = n.d
    tx = n.tx
    ty = n.ty

    n.a = a * m.a + b * m.c
    n.b = a * m.b + b * m.d
    n.c = c * m.a + d * m.c
    n.d = c * m.b + d * m.d
    n.tx = (a, m.tx) + (c * m.ty) + m.tx
    n.ty = (b, m.tx) + (d * m.ty) + m.ty;
end

# n = m * n
function multiply_post!(m::AffineTransform{T}, n::AffineTransform{T}) where {T <: AbstractFloat}
    a = m.a
    b = m.b
    c = m.c
    d = m.d
    tx = m.tx
    ty = m.ty

    n.a = a * n.a + b * n.c
    n.b = a * n.b + b * n.d
    n.c = c * n.a + d * n.c
    n.d = c * n.b + d * n.d
    n.tx = (a, n.tx) + (c * n.ty) + n.tx
    n.ty = (b, n.tx) + (d * n.ty) + n.ty;
end

# out = m * n
function multiply!(m::AffineTransform{T}, n::AffineTransform{T}, out::AffineTransform{T}) where {T <: AbstractFloat}
    out.a = m.a * n.a + m.b * n.c
    out.b = m.a * n.b + m.b * n.d
    out.c = m.c * n.a + m.d * n.c
    out.d = m.c * n.b + m.d * n.d
    out.tx = (m.tx * n.a) + (m.ty * n.c) + n.tx
    out.ty = (m.tx * n.b) + (m.ty * n.d) + n.ty;
end

function invert!(aft::AffineTransform{T}) where {T <: AbstractFloat}
    determinant = 1.0 / (aft.a * aft.d - aft.b * aft.c)
    a = aft.a
    b = aft.b
    c = aft.c
    d = aft.d
    tx = aft.tx
    ty = aft.ty
    aft.a = determinant * d
    aft.b = -determinant * b
    aft.c = -determinant * c
    aft.d = determinant * a
    aft.tx = determinant * (c * ty - d * tx)
    aft.ty = determinant * (b * tx - a * ty);
end

function invert!(aft::AffineTransform{T}, out::AffineTransform{T}) where {T <: AbstractFloat}
    determinant = 1.0 / (aft.a * aft.d - aft.b * aft.c)
    out.a = determinant * aft.d
    out.b = -determinant * aft.b
    out.c = -determinant * aft.c
    out.d = determinant * aft.a
    out.tx = determinant * (aft.c * aft.ty - aft.d * aft.tx)
    out.ty = determinant * (aft.b * aft.tx - aft.a * aft.ty);
end

# Converts either from or to pre or post multiplication.
#     a c
#     b d
# to
#     a b
#     c d
function transpose!(aft::AffineTransform{T}) where {T <: AbstractFloat}
    c = aft.c
    aft.c = aft.b
    aft.b = c;
end

max_stars_repo_name: wdevore/Ranger-Julia-SDL

max_issues_repo_name: wdevore/Ranger-Julia-SDL

max_forks_repo_name: wdevore/Ranger-Julia-SDL


text:
# This example is used to generate plots in the paper, "Reflection from a multi-species material and its transmitted effective wavenumber." Proc. R. Soc. (2018): 20170864.

# Everything related to ploting has been commented so that this example can run with requiring the Plots package. That way, this example can be checked everytime the package is modified.

using EffectiveWaves

# using LaTeXStrings
# using Plots
# height=500
 # pyplot(linewidth=3, size=(2.6*height,height))
 # Plots.scalefontsizes(1.8)

filename="fluid"
mediumname = "water"
## choose material

  # fluid
  medium     = WaterDistilled
  inclusion1 = Glycerol # surfactant ?
  inclusion2 = Hexadecane
  # inclusion2 = SodiumNitrate
  # ωs = 2.0*pi*LinRange(1.0e1,1.0e7,100)

  ωfactor = 4000;
  ω = ωfactor*real(medium.c)/2.0

  vols = 0.01:0.003:0.4
  vol_prop = 0.5
  r1 = 0.1/ωfactor; r2 = 1.0/ωfactor;

  # True wavenumber
  kTs = map(vols) do v
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = v*vol_prop)
    sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = v*(1.0 - vol_prop))
    wavenumber_low_volfrac(ω, medium, [sp1,sp2])
  end
  # low volfrac wavenumber
  kT2s = map(vols) do v
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = v*vol_prop)
    sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = v*(1.0 - vol_prop))
    wavenumber_very_low_volfrac(ω, medium, [sp1,sp2])
  end

  # Approximate wavenumber
  kTLSs = map(vols) do v
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = v*vol_prop)
    sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = v*(1.0 - vol_prop))
    two_species_approx_wavenumber(ω, medium, [sp1,sp2])
  end
  # Approximate challis
  kTCs = map(vols) do v
    sp1 = Specie(ρ=inclusion1.ρ ,r=r1, c=inclusion1.c, volfrac = v*vol_prop)
    sp2 = Specie(ρ=inclusion2.ρ ,r=r2, c=inclusion2.c, volfrac = v*(1.0 - vol_prop))
    wavenumber_challis(ω, medium, [sp1,sp2])
  end

  kTs_arr = [kTs,kT2s, kTLSs,kTCs];
  speed_arr = [ ω./real(ks) for ks in kTs_arr];
  push!(speed_arr, 0 .*real(kTs) .+ real(medium.c))
  atten_arr = imag.(kTs_arr)

  styles = [:solid :dash :dashdot :dashdot :dot]
  ys_arr = speed_arr;
  xs = vols;
  # labs = [L"k_{*}" L"k_{*0}" L"k_{*LS}" L"k_{*C}" mediumname]
  # xlabel = "Total volume fraction"
  # m =5;
  # p1 = plot(xs, ys_arr, xlabel=xlabel, ylabel="sound speed (m/s)", labels=labs
                #  , line = styles, ylims = (minimum(minimum.(ys_arr)), maximum(ys_arr[1])*1.001));

  ys_arr = atten_arr;
  styles = [:solid :dash :dashdot :dashdot]
  # labs = [L"k_{*}" L"k_{*0}" L"k_{*LS}" L"k_{*C}"]
  # p2 = plot(xs, ys_arr, labels=labs, xlabel=xlabel, ylabel="attenuation (1/m)"
                    #   , line = styles )
                #  , ylims = (minimum(ys_arr[1]), maximum(ys_arr[1])*1.04));
  # plot(p1,p2)
  # gui()
  # try mkdir("media") end
  # savefig("media/compare_$(filename)_volfrac.png")
  # savefig("media/compare_$(filename)_volfrac.pdf")

# Plots.scalefontsizes(1/1.8)


text:
using NNlib: conv, maxpool, meanpool, relu, softmax
using ChainRulesCore

# Initialization utils

xavier_uniform(T::Type, dims...; fan_in, fan_out) = sqrt(T(6) / (fan_in + fan_out)) * (2 * rand(T, dims...) .- 1)
xavier_uniform(dims...; fan_in, fan_out) = xavier_uniform(Float32, dims..., fan_in=fan_in, fan_out=fan_out)

# Linear layers

struct Linear{TW, Tb} <: Model
    W::TW
    b::Tb
end

function Linear(T::Type, dims::Pair)
    W = xavier_uniform(T, dims[2], dims[1], fan_in=dims[1], fan_out=dims[2])
    b = zeros(T, dims[2])
    return Linear(W, b)
end

Linear(dims::Pair) = Linear(Float32, dims)

function (l::Linear)(x)
    size_x = size(x)
    x_reshaped = reshape(x, size_x[1], :)
    out_reshaped = l.W * x_reshaped .+ l.b
    out = reshape(out_reshaped, :, size_x[2:end]...)
    return out
end

# Convolutional layers

struct Conv{Tw, Tb} <: Model
    w::Tw
    b::Tb
end

function Conv(T::Type, num_channels::Pair, filter_size::Tuple)
    fan_in = prod(filter_size) * num_channels[1]
    fan_out = prod(filter_size) * num_channels[2]
    w = xavier_uniform(T, filter_size..., num_channels..., fan_in=fan_in, fan_out=fan_out)
    b = zeros(T, 1, 1, num_channels[2])
    return Conv(w, b)
end

Conv(num_channels::Pair, filter_size::Tuple) = Conv(Float32, num_channels, filter_size)

(c::Conv)(x) = conv(x, c.w) .+ c.b

# Composition of layers

struct Compose{T <: Tuple} <: Model
    layers::T
end

Compose(first, rest...) = Compose((first, rest...))

function (m::Compose)(x)
    y = x
    for layer in m.layers
        y = layer(y)
    end
    return y
end

# Parameter-free layers (activations, pooling, reshaping...)

struct ReLU <: Model end

(::ReLU)(x) = relu.(x)

struct SoftMax{D} <: Model end

SoftMax(d=1) = SoftMax{d}()

(::SoftMax{D})(x) where D = softmax(x, dims=D)

struct MaxPool{S} <: Model end

MaxPool(s) = MaxPool{s}()

(p::MaxPool{S})(x) where S = maxpool(x, S)

struct MeanPool{S} <: Model end

MeanPool(s) = MeanPool{s}()

(p::MeanPool{S})(x) where S = meanpool(x, S)

struct Reshape{S} <: Model end

Reshape(s) = Reshape{s}()

(r::Reshape{S})(x) where S = reshape(x, S)

# Dropout

_dropout_shape(s, ::Colon) = size(s)
_dropout_shape(s, dims) = tuple((i in dims ? si : 1 for (i, si) in enumerate(size(s)))...)

function _dropout_mask(x::AbstractArray{T}, p; dims=:) where T
    y = rand(eltype(x), _dropout_shape(x, dims)...)
    y .= ifelse.(y .> p, T(1/(1 - p)), T(0))
    return y
end

function dropout(x, p; dims=:)
    if !is_training_mode()
        return x
    end
    y = _dropout_mask(x, p, dims=dims)
    return x .* y
end

function ChainRulesCore.rrule(::typeof(dropout), x, p; dims=:)
    if !is_training_mode()
        return x, c -> (ChainRulesCore.NoTangent(), c, ChainRulesCore.ZeroTangent())
    end
    y = _dropout_mask(x, p, dims=dims)
    return x .* y, c -> (ChainRulesCore.NoTangent(), c .* y, ChainRulesCore.ZeroTangent())
end

struct Dropout{P, D} <: Model end

Dropout(p, dims=:) = Dropout{p, dims}()

(d::Dropout{P, D})(x) where {P, D} = dropout(x, P, dims=D)


text:
"""
```
mutable struct MeansBands
```

Stores the means and bands of results for a particular set of outputs from the forecast step.

Specifically, forecasts can be made for any element in the Cartesian product of 4 sets:

1. `input_type`: some subset of the parameter draws from the estimation
   step. See `forecast_one` for all possible options.

2. `cond_type`: conditional type. See `forecast_one` for all possible options.

3. *product*: a particular result computed in the forecast. This could be one of
   the following:

```
  - `hist`: smoothed histories
  - `forecast`: forecasted values
  - `shockdec`: shock decompositions
  - `irf`: impulse responses
```

4. variable *class*: the category in which a particular variable, like `:y_t`,
   falls. Options are:

```
  - `state`: state (from `m.endogenous_states` or `m.endogenous_states_augmented`)
  - `obs`: observable (from `m.observables`)
  - `pseudo`: pseudoobservable (from `pseudo_measurement` equation)
  - `shock`: shock (from `m.exogenous_shocks`)
```

Note that the Cartesian product (product x class) is the set of options for
`output_vars` in the `forecast_one` function signature.

### Fields

- `metadata::Dict{Symbol,Any}`: Contains metadata keeping track of the
  `input_type`, `cond_type`, product (history, forecast, shockdec,
  etc), and variable class (observable, pseudoobservable, state, etc)
  stored in this `MeansBands` structure.
- `means::DataFrame`: a `DataFrame` of the mean of the time series
- `bands::Dict{Symbol,DataFrame}`: a `Dict` mapping variable names to
  `DataFrame`s containing confidence bands for each variable. See
  `find_density_bands` for more information.
"""
mutable struct MeansBands
    metadata::Dict{Symbol,Any}
    means::DataFrame
    bands::Dict{Symbol,DataFrame}

    function MeansBands(key, means, bands)

        if !isempty(bands)

            # assert that means and bands fields have the same keys (provide info for same products)
            @assert sort(setdiff(propertynames(means),[:date])) == sort(collect(keys(bands)))

            # check to make sure that # of periods in all dataframes are the same
            n_periods_means = size(means,1)
            for df in values(bands)
                n_periods_bands = size(df,1)
                @assert(n_periods_means == n_periods_bands,
                        "means and bands must have same number of periods")
            end
        end

        new(key, means, bands)
    end
end

# A dummy MeansBands object
function MeansBands()
    metadata   = Dict(:class => :none, :product => :none,
                      :cond_type => :none, :para => :none,
                      :date_inds => OrderedDict{Date, Int}(Date(1) => 1),
                      :indices => Dict{Symbol, Int}(:none => 1))

    means = DataFrame(date = [Dates.Date(0)], none = [0.0])
    bands = Dict{Symbol,DataFrame}(:none => DataFrame(date = [Dates.Date(0)]))

    MeansBands(metadata, means, bands)
end

function Base.show(io::IO, mb::MeansBands)
    @printf io "MeansBands\n"
    @printf io "  class: %s\n"   get_class(mb)
    @printf io "  product: %s\n" get_product(mb)
    @printf io "  cond: %s\n"    get_cond_type(mb)
    @printf io "  para: %s\n"    get_para(mb)
    if mb.metadata[:product] != :trend && mb.metadata[:product] != :irf
        if isempty(mb.metadata[:date_inds])
            @printf io "  dates: []\n"
        else
            @printf io "  dates: %s - %s\n" startdate_means(mb) enddate_means(mb)
        end
    end
    @printf io "  # of variables: %s\n" n_vars_means(mb)
    @printf io "  bands: %s\n" which_density_bands(mb, uniquify=true)
    if haskey(mb.metadata, :scenario_key)
        @printf io "  scenario: %s\n" get_scenario_key(mb)
    end
end

"""
```
Base.isempty(mb::MeansBands)
```

Returns whether the `mb` object in question is a dummy.
"""
function Base.isempty(mb::MeansBands)

    return get_class(mb) == :none && get_product(mb) == :none &&
        startdate_means(mb) == Dates.Date(0) &&
        collect(propertynames(mb.means)) ==  [:date, :none] &&
        collect(keys(mb.bands)) ==  [:none]
end

"""
```
Base.cat(mb1::MeansBands, mb2::MeansBands; out_product = Symbol(),
    forecast_string = "")
```

Concatenate 2 compatible `MeansBands` objects together by date. 2
`MeansBands` objects are defined to be compatible if the class of
variables is the same, the conditional type is the same, and the
input_type for the forecast used to create the two `MeansBands` object
is the same. Furthermore, we require that the dates covered by each
`MeansBands` object form a continguous interval.

### Inputs

- `mb1::MeansBands`
- `mb2::MeansBands`

Note that the dates in `mb1` should come chronologically first, with the dates
in `mb2` beginning 1 period after the final period in `mb1`.

### Keyword Arguments

- `out_product::Symbol`: desired product of the resulting concatenated
`MeansBands` object. This argument is required if `mb1` and `mb2` do not
represent a history and a forecast respectively
- `forecast_string::String`: desired `forecast_string` of the resulting
  concatenated `MeansBands`. This argument is recommended (but not required) if
  the `forecast_string`s of `mb1` and `mb2` do not match. If this is the case
  but `forecast_string` is not provided, `mb1`'s `forecast_string` will be used.
"""
function Base.cat(mb1::MeansBands, mb2::MeansBands;
                  out_product::Symbol = Symbol(),
                  forecast_string::String = "")

    # If either mb1 or mb2 is empty, return just the other one
    if isempty(mb1)
        return mb2
    elseif isempty(mb2)
        return mb1
    end

    # Assert class, cond type, para, and scenario info are the same
    @assert get_class(mb1) == get_class(mb2)
    @assert get_cond_type(mb1) == get_cond_type(mb2)
    @assert get_para(mb1) == get_para(mb2)
    @assert haskey(mb1.metadata, :scenario_key) == haskey(mb2.metadata, :scenario_key) ==
        haskey(mb1.metadata, :scenario_vint) == haskey(mb2.metadata, :scenario_vint)
    if haskey(mb1.metadata, :scenario_key)
        @assert mb1.metadata[:scenario_key] == mb2.metadata[:scenario_key]
        @assert mb1.metadata[:scenario_vint] == mb2.metadata[:scenario_vint]
    end

    # Assert dates are contiguous
    last_mb1_date  = enddate_means(mb1)
    first_mb2_date = startdate_means(mb2)
#    @assert iterate_quarters(last_mb1_date, 1) == first_mb2_date

    # compute means field
    means = vcat(mb1.means, mb2.means)
    # na2nan!(means) # Removed b/c DataFrames uses missing instead of NA, and missings are already handled

    # compute bands field
    bands = Dict{Symbol, DataFrame}()

    mb1vars = collect(keys(mb1.bands))
    mb2vars = collect(keys(mb2.bands))
    nperiods_mb1 = length(mb1.metadata[:date_inds])
    nperiods_mb2 = length(mb2.metadata[:date_inds])

    bothvars = intersect(mb1vars, mb2vars)
    for var in setdiff(union(keys(mb1.bands), keys(mb2.bands)), [:date])
        bands[var] = if var in bothvars
            vcat(mb1.bands[var], mb2.bands[var])
        elseif var in setdiff(mb1vars, mb2vars)
            vcat(mb1vars[var], fill(NaN, nperiods_mb2))
        else
            vcat(fill(NaN, nperiods_mb1), mb2vars[var])
        end
        # na2nan!(bands[var]) # Removed b/c DataFrames uses missing instead of NA, and missings are already handled
    end

    # compute metadata
    # product
    mb1_product = get_product(mb1)
    mb2_product = get_product(mb2)
    product = if mb1_product == :hist
        if mb2_product == :forecast
            :histforecast
        elseif mb2_product == :bddforecast
            :bddforecast
        end
    elseif mb1_product == :histut
        if mb2_product == :forecastut
            :histforecastut
        elseif mb2_product == :bddforecast
            :bddforecastut
        end
    elseif mb1_product == :hist4q
        if mb2_product == :forecast4q
            :histforecast4q
        elseif mb2_product == :bddforecast4q
            :bddhistforecast4q
        end
    elseif mb1_product == mb2_product
        mb1_product
    else
        @assert !isempty(out_product) "Please supply a product name for the output MeansBands"
        out_product
    end

    # date indices
    date_indices = Dict(d::Dates.Date => i::Int for (i, d) in enumerate(means[!, :date]))

    # variable indices
    indices = Dict(var::Symbol => i::Int for (i, var) in enumerate(propertynames(means)))

    # forecast string
    if isempty(forecast_string) && (mb1.metadata[:forecast_string] != mb2.metadata[:forecast_string])
        @warn "No forecast_string provided: using $(mb1.metadata[:forecast_string])"
    end
    forecast_string = mb1.metadata[:forecast_string]

    metadata = Dict{Symbol, Any}(
                   :para            => get_para(mb1),
                   :cond_type       => get_cond_type(mb1),
                   :class           => get_class(mb1),
                   :product         => product,
                   :indices         => indices,
                   :forecast_string => forecast_string,
                   :date_inds       => sort(date_indices, by = x -> date_indices[x]))
    if haskey(mb1.metadata, :scenario_key)
        metadata[:scenario_key] = mb1.metadata[:scenario_key]
        metadata[:scenario_vint] = mb1.metadata[:scenario_vint]
    end

    # construct the new MeansBands object and return
    MeansBands(metadata, means, bands)
end

###################################
## METADATA
###################################

"""
```
get_class(mb::MeansBands)
```
Returns the class of variables contained this `MeansBands` object (observables, pseudoobservables)
"""
get_class(mb::MeansBands) = mb.metadata[:class]

"""
```
get_product(mb::MeansBands)
```
Returns the product stored in this `MeansBands` object (history, forecast, shock decompostion, etc).
"""
get_product(mb::MeansBands) = mb.metadata[:product]

"""
```
get_cond_type(mb::MeansBands)
```
Returns the conditional type of the forecast used to create this `MeansBands` object.
"""
get_cond_type(mb::MeansBands) = mb.metadata[:cond_type]

"""
```
get_para(mb::MeansBands)
```
Returns the `input_type` from the forecast used to create this `MeansBands` object.
"""
get_para(mb::MeansBands) = mb.metadata[:para]

"""
```
get_shocks(mb::MeansBands)
```

Returns a list of shock names that are used for the shock
decomposition stored in a shock decomposition or irf MeansBands object `mb`.
"""
function get_shocks(mb::MeansBands)
    @assert get_product(mb) in [:shockdec, :irf] "Function only for shockdec or irf MeansBands objects"
    varshocks = setdiff(propertynames(mb.means), [:date])
    unique(map(x -> Symbol(split(string(x), DSGE_SHOCKDEC_DELIM)[2]), varshocks))
end

"""
```
parse_mb_colname(s::Symbol)
```

`MeansBands` column names are saved in the format
`\$var\$DSGE_SHOCKDEC_DELIM\$shock`. `parse_mb_colname` returns (`var`,
`shock`).
"""
function parse_mb_colname(s::Symbol)
    map(Symbol, split(string(s), DSGE_SHOCKDEC_DELIM))
end

"""
```
get_variables(mb::MeansBands)
```

Returns a list of variable names that are used for the shock
decomposition stored in a shock decomposition or irf MeansBands object `mb`.
"""
function get_variables(mb::MeansBands)
    @assert get_product(mb) in [:shockdec, :irf] "Function only for shockdec or irf MeansBands objects"
    varshocks = setdiff(propertynames(mb.means), [:date])
    unique(map(x -> Symbol(split(string(x), DSGE_SHOCKDEC_DELIM)[1]), varshocks))
end

"""
```
get_scenario_key(mb::MeansBands)
```

If `mb` is an alternative scenario `MeansBands`, return the scenario key.
"""
function get_scenario_key(mb::MeansBands)
    @assert haskey(mb.metadata, :scenario_key) "Function only for scenario MeansBands objects"
    mb.metadata[:scenario_key]
end

###################################
## MEANS
###################################

"""
```
n_vars_means(mb::MeansBands)
````

Get number of variables (`:y_t`, `:OutputGap`, etc) in `mb.means`.
"""
function n_vars_means(mb::MeansBands)
    length(get_vars_means(mb))
end

"""
```
get_vars_means(mb::MeansBands)
````

Get variables (`:y_t`, `:OutputGap`, etc) in `mb.means`. Note that
`mb.metadata[:indices]` is an `OrderedDict`, so the keys will be in the correct
order.
"""
function get_vars_means(mb::MeansBands)
    collect(keys(mb.metadata[:indices]))
end

"""
```
n_periods_means(mb::MeansBands)
```

Get number of periods in `mb.means`.
"""
n_periods_means(mb::MeansBands) = size(mb.means,1)

"""
```
startdate_means(mb::MeansBands)
```

Get first period in`mb.means`. Assumes `mb.means[product]` is already sorted by
date.
"""
startdate_means(mb::MeansBands) = mb.means[!, :date][1]

"""
```
enddate_means(mb::MeansBands)
```

Get last period for which `mb` stores means. Assumes `mb.means[product]` is
already sorted by date.
"""
enddate_means(mb::MeansBands) = mb.means[!, :date][end]

"""
```
get_shockdec_means(mb::MeansBands, var::Symbol;
    shocks::Vector{Symbol} = Vector{Symbol}())
```

Return the mean value of each shock requested in the shock decomposition of a
particular variable. If `shocks` is empty, returns all shocks.
"""
function get_shockdec_means(mb::MeansBands, var::Symbol; shocks::Vector{Symbol} = Vector{Symbol}())

    # Extract the subset of columns relating to the variable `var` and the shocks listed in `shocks.`
    # If `shocks` not provided, give all the shocks
    var_cols = collect(propertynames(mb.means))[findall([contains(string(col), string(var)) for col in propertynames(mb.means)])]
    if !isempty(shocks)
        var_cols = [col -> contains(string(col), string(shock)) ? col : nothing for shock in shocks]
    end

    # Make a new DataFrame with column the column names
    out = DataFrame()
    for col in var_cols
        shockname = split(string(col), DSGE_SHOCKDEC_DELIM)[2]
        out[Symbol(shockname)] = mb.means[col]
    end

    return out
end

###################################
## BANDS
###################################

"""
```
n_vars_bands(mb::MeansBands)
```

Get number of variables (`:y_t`, `:OutputGap`, etc) for which `mb`
stores bands for the specified `product` (`hist`, `forecast`, `shockdec`, etc).
"""
n_vars_bands(mb::MeansBands) = length(mb.bands)

"""
```
n_periods_bands(mb::MeansBands)
```

Get number of periods for which `mb` stores bands for the specified
product` (`hist`, `forecast`, `shockdec`, etc).
"""
function n_periods_bands(mb::MeansBands)
    size(mb.bands[collect(keys(mb.bands))[1]],1)
end

"""
```
startdate_bands(mb::MeansBands)
```

Get first period for which `mb` stores bands. Assumes `mb.bands` is already sorted by date.
Also assumes the startdate is same for all observables in MeansBands
"""
startdate_bands(mb::MeansBands) = if isempty(mb)
    Date(0000, 1, 1)
else
    first(mb.bands)[2][:, :date][1]
end

"""
```
enddate_bands(mb::MeansBands)
```

Get last period in `mb.bands`. Assumes `mb.bands` is already sorted by date.
Also assumes the startdate is same for all observables in MeansBands
"""
enddate_bands(mb::MeansBands) = if isempty(mb)
    Date(0000, 1, 1)
else
    first(mb.bands)[2][:, :date][end]
end

"""
```
which_density_bands(mb, uniquify=false)
```

Return a list of the bands stored in mb.bands. If `uniquify = true`,
strips \"upper\" and \"lower\" band tags and returns unique list of percentage values.
"""
function which_density_bands(mb::MeansBands; uniquify=false, ordered=true)

    if isempty(mb.bands)
        return String[]
    else
        # extract one of the keys in mb.bands
        var  = collect(keys(mb.bands))[1]

        # get all the columns in the corresponding dataframe that aren't dates
        strs = map(string,propertynames(mb.bands[var]))
        strs = setdiff(strs, ["date"])

        lowers = strs[map(occursin, repeat([r"LB"], outer=length(strs)), strs)]
        uppers = strs[map(occursin, repeat([r"UB"], outer=length(strs)), strs)]

        # sort
        if ordered
            sort!(lowers, rev=true)
            sort!(uppers)
        end

        # return both upper and lower bands, or just percents, as desired
        strs = if uniquify
            sort([convert(String, split(x, " ")[1]) for x in lowers])
        else
            [lowers; uppers]
        end
        return strs
    end
end

"""
```
get_shockdec_bands(mb, var; shocks = Vector{Symbol}(), bands = Vector{Symbol}())
```

Return a `Dict{Symbol,DataFrame}` mapping shock names to bands for a particular
variable.

### Inputs

- `mb::MeansBands`
- `var::Symbol`: the variable of interest (eg the state `:y_t`, or observable
  `:obs_hours`)

### Keyword Arguments

- `shocks::Vector{Symbol}`: subset of shock names for which to return bands. If
  empty, `get_shockdec_bands` returns all bands
- `bands::Vector{Symbol}`: subset of bands stored in the DataFrames of
  `mb.bands` to return
"""
function get_shockdec_bands(mb::MeansBands, var::Symbol;
                            shocks::Vector{Symbol} = Vector{Symbol}(),
                            bands::Vector{Symbol} = Vector{Symbol}())

    @assert get_product(mb) == :shockdec

    # Extract the subset of columns relating to the variable `var` and the shocks listed in `shocks.`
    # If `shocks` not provided, give all the shocks
    var_cols = collect(keys(mb.bands))[findall([contains(string(col), string(var)) for col in keys(mb.bands)])]
    if !isempty(shocks)
        var_cols = [col -> contains(string(col), string(shock)) ? col : nothing for shock in shocks]
    end

    # Extract the subset of bands we want to return. Return all bands if `bands` not provided.
    bands_keys = if isempty(bands)
        propertynames(mb.bands[var_cols[1]])
    else
        [[Symbol("$(100x)% LB") for x in bands]; [Symbol("$(100x)% UB") for x in bands]]
    end

    # Make a new dictionary mapping shock names to bands
    out = Dict{Symbol, DataFrame}()
    for col in var_cols
        shockname = parse_mb_colname(col)[2]
        out[shockname] = mb.bands[col][bands_keys]
    end

    return out
end

################################################
## EXTRACTING VARIABLES
################################################

"""
```
prepare_meansbands_table_timeseries(mb, var)
```

Returns a `DataFrame` of means and bands for a particular time series variable
(either `hist` or `forecast` of some type). Columns are sorted such that the
bands are ordered from smallest to largest, and the means are on the far
right. For example, a `MeansBands` containing 50\\% and 68\\% bands would be
ordered as follows: [68\\% lower, 50\\% lower, 50\\% upper, 68\\% upper, mean].

### Inputs

- `mb::MeansBands`: time-series MeansBands object
- `var::Symbol`: an economic variable stored in `mb`. If `mb` stores
  observables, `var` would be an element of `propertynames(m.observables)`. If
  it stores pseudo-observables, `var` would be the name of a
  pseudo-observable defined in the pseudo-measurement equation.

### Keyword Arguments

- `bands_pcts::Vector{String}`: vector of (uniquified) band percentiles to
  include in the table
"""
function prepare_meansbands_table_timeseries(mb::MeansBands, var::Symbol;
                                             bands_pcts::Vector{String} = which_density_bands(mb, uniquify = true))

    @assert get_product(mb) in [:hist, :histut, :hist4q, :forecast, :forecastut, :forecast4q,
                                :bddforecast, :bddforecastut, :bddforecast4q,
                                :histforecast, :histforecastut, :histforecast4q,
                                :bddhistforecast, :bddhistforecastut, :bddhistforecast4q,
                                :trend, :dettrend] "prepare_meansbands_table_timeseries can only be used for time-series products"
    @assert var in get_vars_means(mb) "$var is not stored in this MeansBands object"

    # Get bands
    uppers = sort!([pct * " UB" for pct in bands_pcts], rev = true)
    lowers = sort!([pct * " LB" for pct in bands_pcts])
    my_bands = map(Symbol, vcat(lowers, uppers))

    # Extract this variable from Means and bands
    means = mb.means[:, [:date, var]]
    bands = mb.bands[var][:, [:date; my_bands]]

    # Join so mean is on far right and date is on far left
    df = isdefined(DataFrames, :innerjoin) ? innerjoin(bands, means, on = :date) : join(bands, means, on = :date)
    rename!(df, var => Symbol("mean"))

    return df
end

"""
```
prepare_means_table_shockdec(mb_shockdec, mb_trend, mb_dettrend, var;
    shocks = get_shocks(mb_shockdec), mb_forecast = MeansBands(),
    mb_hist = MeansBands(), detexify = true, groups = [])
```

Returns a `DataFrame` representing a detrended shock decompostion for
the variable `var`. The columns of this dataframe represent the
contributions of each shock in `shocks` (or all shocks, if the keyword
argument is omitted) and the deterministic trend.

### Inputs

- `mb_shockdec::MeansBands`: a `MeansBands` object for a shock decomposition
- `mb_trend::MeansBands`: a `MeansBands` object for a trend  product.
- `mb_dettrend::MeansBands`: a `MeansBands` object for a deterministic trend
  product.
- `var::Symbol`: name of economic variable for which to return the means and bands table

### Keyword Arguments

- `shocks::Vector{Symbol}`: If `mb` is a shock decomposition, this is
  an optional list of shocks to print to the table. If omitted, all
  shocks will be printed.
- `mb_forecast::MeansBands`: a `MeansBands` object for a forecast.
- `mb_hist::MeansBands`: a `MeansBands` object for smoothed states.
- `detexify::Bool`: whether to remove Unicode characters from shock names
- `groups::Vector{ShockGroup}`: if provided, shocks will be grouped accordingly
"""
function prepare_means_table_shockdec(mb_shockdec::MeansBands, mb_trend::MeansBands,
                                      mb_dettrend::MeansBands, var::Symbol;
                                      shocks::Vector{Symbol} = get_shocks(mb_shockdec),
                                      mb_forecast::MeansBands = MeansBands(),
                                      mb_hist::MeansBands = MeansBands(),
                                      detexify_shocks::Bool = true,
                                      groups::Vector{ShockGroup} = ShockGroup[])

    @assert get_product(mb_shockdec) == :shockdec "The first argument must be a MeansBands object for a shockdec"
    @assert get_product(mb_trend)    == :trend    "The second argument must be a MeansBands object for a trend"
    @assert get_product(mb_dettrend) == :dettrend "The third argument must be a MeansBands object for a deterministic trend"

    # Get the variable-shock combinations we want to print
    varshocks = [Symbol("$var" * DSGE_SHOCKDEC_DELIM * "$shock") for shock in shocks]

    # Fetch the columns corresponding to varshocks
    df_shockdec = mb_shockdec.means[!,union([:date], varshocks)]
    df_trend    = mb_trend.means[!,[:date, var]]
    df_dettrend = mb_dettrend.means[!,[:date, var]]

    # Line up dates between trend, dettrend and shockdec
    has_ij = isdefined(DataFrames, :innerjoin)
    if has_ij
        df_shockdec = innerjoin(df_shockdec, df_trend, on = :date)
        rename!(df_shockdec, var => :trend)
        df_shockdec = innerjoin(df_shockdec, df_dettrend, on = :date)
        rename!(df_shockdec, var => :dettrend)
    else
        df_shockdec = join(df_shockdec, df_trend, on = :date, kind = :inner)
        rename!(df_shockdec, var => :trend)
        df_shockdec = join(df_shockdec, df_dettrend, on = :date, kind = :inner)
        rename!(df_shockdec, var => :dettrend)
    end

    # Add each shock's contribution and deterministic trend to output DataFrame
    df = DataFrame(date = df_shockdec[!,:date])
    for col in setdiff(propertynames(df_shockdec), [:date, :trend])
        df[!,col] = df_shockdec[!,col]
    end
    df_shockdec[!,:dettrend] = df_shockdec[!,:dettrend]

    # Rename columns to just the shock names
    map(x -> rename!(df, x => parse_mb_colname(x)[2]), setdiff(propertynames(df), [:date, :trend, :dettrend]))

    # If mb_forecast and mb_hist are passed in, add the detrended time series
    # mean of var to the table
    if !(isempty(mb_forecast) && isempty(mb_hist))
        mb_timeseries = cat(mb_hist, mb_forecast)

        # Truncate to just the dates we want
        startdate = df[1, :date]
        enddate   = df[end, :date]
        df_mean   = mb_timeseries.means[startdate .<= mb_timeseries.means[!,:date] .<= enddate, [:date, var]]

        df_shockdec = has_ij ? innerjoin(df_shockdec, df_mean, on = :date) :
            join(df_shockdec, df_mean, on = :date, kind = :inner)
        df[!,:detrendedMean] = df_shockdec[!,var] - df_shockdec[!,:trend]
    end

    # Group shocks if desired
    nperiods = size(df, 1)
    v0 = zeros(nperiods)
    for group in groups
        # Sum shock values for each group
        shock_vectors = [df[!,shock] for shock in group.shocks]
        shock_sum = reduce(+, shock_vectors; init = v0)
        df[!,Symbol(group.name)] = shock_sum

        # Delete original (ungrouped) shocks from df
        select!(df, Not(setdiff(group.shocks, [Symbol(group.name)])))
    end

    # Remove Unicode characters from shock names
    if detexify_shocks
        for x in setdiff(propertynames(df), [:date, :trend, :dettrend])
            x_detexed = detexify(x)
            if x != x_detexed
                rename!(df, x => x_detexed)
            end
        end
    end

    return df
end

#="""
```
prepare_meansbands_table_irf(mb, var, shock)
```

Returns a `DataFrame` of means and bands for a particular impulse
response function of variable (observable, pseudoobservable, or state)
`v` to shock `s`. Columns are sorted such that the bands are ordered from
smallest to largest, and the means are on the far right. For example,
a MeansBands object containing 50\\% and 68\\% bands would be ordered as
follows: [68\\% lower, 50\\% lower, 50\\% upper, 68\\% upper, mean].

### Inputs
- `mb::MeansBands`: time-series MeansBands object
- `var::Symbol`: an economic variable stored in `mb`. If `mb` stores
  observables, `var` would be an element of `propertynames(m.observables)`. If
  it stores pseudoobservables, `var` would be the name of a
  pseudoobservable defined in the pseudomeasurement equation.
"""
function prepare_meansbands_table_irf(mb::MeansBands, shock::Symbol, var::Symbol)

    @assert get_product(mb) in [:irf] "prepare_meansbands_table_irf can only be used for irfs"
    @assert var in get_vars_means(mb) "$var is not stored in this MeansBands object"

    # get the variable-shock combination we want to print
    # varshock = Symbol["$var" * DSGE_SHOCKDEC_DELIM * "$shock" for var in vars]
    varshock = Symbol("$var" * DSGE_SHOCKDEC_DELIM * "$shock")

    # extract the means and bands for this irf
    df = mb.bands[varshock][map(Symbol, which_density_bands(mb))]
    df[:mean] = mb.means[varshock]

    return df
end

function prepare_meansbands_table_irf(mb::MeansBands, shock::Symbol, vars::Vector{Symbol})

    # Print all vars by default
    if isempty(vars)
        vars = get_variables(mb)
    end

    # Make dictionary to return
    irfs = Dict{Symbol, DataFrame}()
    df = DataFrame()

    # Make tables for each irf
    for var in vars
        @show df
        if isempty(df)
            df = prepare_meansbands_table_irf(mb, shock, var)
        else
            join(df, prepare_meansbands_table_irf(mb, shock, var), on = :date)
        end
    end

    return df
end=#


text:
using Plots

blue = "#8DC0FF"
red = "#FFAEA6"

t = 0:10:3040
ϵ₁ = randn(length(t))*15     # spread on Dolphin IQ
ϵ₂ = randn(length(t))*20     # spreak on Human IQ

Y₁ = dolphinsIQ = t/12 + ϵ₁
Y₂ = humanIQ = t/20 + ϵ₂

scatter(t,Y₁; label="Dolphins", color=blue,
  ylabel="IQ (-)", xlabel ="Time (year BC)", legend=:topleft)
scatter!(t,Y₂; label="Humans", color=red)

# Solution
function ridgeRegression(X,Y;λ=0)
  I = one(Y*Y')
  β = (X*X' + λ*I)^-1*X'Y
end

# Solution
function ridgeRegression2(X,Y;λ=0)
  I = one(Y*Y')
  β = (X*X' + λ*I)^-1*X'Y
end

β₁ =  ridgeRegression(t,Y₁;λ=1e-3)
β₂ =  ridgeRegression(t,Y₂;λ=1e-3)

plot!(time,β₁*time;label="Prediction Dolphins",color=blue)
plot!(time,β₂*time;label="Prediction Human", color=red)


text:
function range_to_ind(range)
    length(range) == 1 || throw(ArgumentError("length(range) != 1"))
    first(range)
end


text:
# simple.jl
# A collection of simple domains.

###########################
# The unit ball and sphere
###########################

"""
The unit ball (of radius 1) in `N` dimensions.
"""
struct UnitHyperBall{N,T} <: EuclideanDomain{N,T} end

UnitHyperBall{N}() where N = UnitHyperBall{N,Float64}()

const UnitDisk{T} = UnitHyperBall{2,T}
const UnitBall{T} = UnitHyperBall{3,T}

indomain(x, ::UnitHyperBall) = norm(x) <= 1

approx_indomain(x, ::UnitHyperBall, tolerance) = norm(x) <= 1+tolerance

isempty(::UnitHyperBall) = false

show(io::IO, d::UnitHyperBall{N}) where {N} = print(io, "the $(N)-dimensional unit ball")

# We choose the origin here
point_in_domain(d::UnitHyperBall) = zero(eltype(d))

###########################
# An n-dimensional simplex
###########################

struct UnitSimplex{N,T} <: EuclideanDomain{N,T} end
UnitSimplex{N}() where N = UnitSimplex{N,Float64}()

indomain(x, ::UnitSimplex) = mapreduce( t-> t >= 0, &, x) && norm(x,1) <= 1

approx_indomain(x, ::UnitSimplex, tolerance) = mapreduce( t-> t >= -tolerance, &, x) && norm(x,1) <= 1+tolerance

isempty(::UnitSimplex) = false

# We pick the origin, because it belongs to the domain regardless of what T is
point_in_domain(d::UnitSimplex) = zero(eltype(d))

##############
# A cylinder
##############

cylinder(::Type{T} = Float64) where {T} = UnitDisk{T}() × UnitInterval{T}()

cylinder(radius::T, length::T) where {T} = radius * UnitDisk{T}() × (0 .. length)

"""
    Point(x)

represents a single point at `x`.
"""
struct Point{T} <: Domain{T}
    x::T
end

convert(::Type{Number}, d::Point{<:Number}) = d.x
convert(::Type{N}, d::Point{<:Number}) where N<:Number = convert(N, convert(Number, d.x))
Number(d::Point) = convert(Number, d)

convert(::Type{Domain{T}}, d::Point{T}) where T = d
convert(::Type{Domain{T}}, d::Point) where T = Point(T(d.x))
convert(::Type{Domain}, c::Number) = Point(c)
convert(::Type{Domain{T}}, c::Number) where T = Point(convert(T,c))
convert(::Type{Domain}, s::Set) = UnionDomain(map(Domain,collect(s)))
Domain(d) = convert(Domain, d)

==(a::Point,b::Point) = a.x == b.x
indomain(x, d::Point) = x == d.x
isempty(::Point) = false

approx_indomain(x, d::Point, tolerance) = norm(x-d.x) <= tolerance

isopen(d::Point) = false
isclosed(d::Point) = true

point_in_domain(d::Point) = d.x

for op in (:*,:+,:-)
    @eval begin
        $op(c::Number, d::Point)  = Point($op(c,d.x))
        $op(d::Point,  c::Number) = Point($op(d.x,c))
    end
end

/(d::Point,  c::Number) = Point(d.x/c)
\(c::Number, d::Point)  = Point(c\d.x)

for op in (:+,:-)
    @eval $op(a::Point, b::Point) = Point($op(a.x,b.x))
end

for op in (:*,:+)
    @eval begin
        $op(a::Point, v::AbstractVector) = map(y->$op(a,y),v)
        $op(v::AbstractVector, a::Point) = map(y->$op(y,a),v)
    end
end

function setdiff(d::Interval{L,R,T}, p::Point{T}) where {L,R,T}
    a = leftendpoint(d)
    b = rightendpoint(d)

    a == p.x && return Interval{:open,R}(a,b)
    a < p.x < b && return UnionDomain(Interval{L,:open}(a,p.x)) ∪ UnionDomain(Interval{:open,R}(p.x,b))
    b == p.x && return Interval{L,:open}(a,b)

    return d
end

issubset(p::Point, d::Domain) = p.x ∈ d


text:
module DatagenCopulaBased
  using Distributions
  using Combinatorics
  using HypothesisTests
  using HCubature
  using StatsBase
  using QuadGK
  using Roots
  using SpecialFunctions
  using Random
  using LinearAlgebra
  using Distributed
  using SharedArrays

  include("sampleunivdists.jl")
  include("archcopcorrelations.jl")
  include("archcopulagendat.jl")
  include("nestedarchcopulagendat.jl")
  include("chaincopulagendat.jl")
  include("subcopulasgendat.jl")
  include("copulagendat.jl")
  include("marshallolkincopcor.jl")

  export archcopulagen, chaincopulagen, nestedarchcopulagen
  export cormatgen, cormatgen_constant, cormatgen_toeplitz, convertmarg!, gcop2arch
  export cormatgen_constant_noised, cormatgen_toeplitz_noised, cormatgen_rand
  export cormatgen_two_constant, cormatgen_two_constant_noised
  export tstudentcopulagen, gausscopulagen
  export frechetcopulagen, marshallolkincopulagen, chainfrechetcopulagen
  export gcop2tstudent, gcop2frechet, gcop2marshallolkin
end


text:
# Author: Ritchie Lee, ritchie.lee@sv.cmu.edu
# Date: 11/21/2014

#ACASX implementation: based on CCAS interface to libCAS

module ACASXCommonImpl

export
    addObserver,

    ACASXInput,
    ACASXOutput,
    ACASXCoordRecord

    #use explicit form for these
    #initialize
    #step
    #update_from_coord!
    #update_to_coord!

using AbstractCollisionAvoidanceSystemImpl
using AbstractCollisionAvoidanceSystemInterfaces
using CommonInterfaces

using AbstractCASCoordImpl

using CASCoordination
using CASInterface
using RLESUtils, Observers

import CASInterface.reset! #extend

typealias ACASXInput CASInterface.Input
typealias ACASXOutput CASInterface.Output

type ACASXCoordRecordIntruder
  id::UInt32
  cvc::UInt8
  vrc::UInt8
  vsb::UInt8

  ACASXCoordRecordIntruder(id::Int) = new(id, 0x0, 0x0, 0x0)
end

type ACASXCoordRecord
  my_id::UInt32
  modes::UInt32
  equipage::EQUIPAGE
  quant::UInt8
  sensitivity_index::UInt8
  protection_mode::UInt8
  intruders::Vector{ACASXCoordRecordIntruder}
end
function ACASXCoordRecord(my_id::Int64, equipage::EQUIPAGE, quant::Int64, max_intruders::Int64)
  ACASXCoordRecord(UInt32(my_id), equipage, UInt8(quant), max_intruders)
end
function ACASXCoordRecord(my_id::UInt32, equipage::EQUIPAGE, quant::UInt8, max_intruders::Int64)
  ACASXCoordRecord(my_id, my_id, equipage, quant, 0x0, 0x0,
                   [ACASXCoordRecordIntruder(i) for i = 1:max_intruders])
end

getListId(list_owner_id::Integer,id::Integer) = id < list_owner_id ? id : id - 1

function update_from_coord!(input::ACASXInput, coord::AbstractCASCoord, my_id::Int64)
  #Grab from coordination object
  coord_records = CASCoordination.getAll(coord)

  #Augment the input with info from coord
  #note: looping over intruders skips self
  for intruder in input.intruders

    my_list_id = getListId(intruder.id, my_id)
    record = coord_records[intruder.id]

    #intruder-shared
    intruder.equipage             = record.equipage
    intruder.quant                = record.quant
    intruder.sensitivity_index    = record.sensitivity_index
    intruder.protection_mode      = record.protection_mode

    #quantize sensor reading depending on equipage
    #TODO: should these be moved to Sensor?
    intruder.z = quantize(intruder.z, Float64(intruder.quant))

    #intruder-specific
    intruder_self = record.intruders[my_list_id] #self in intruders' record

    @assert my_id == intruder_self.id #sanity check the record

    intruder.cvc = intruder_self.cvc
    intruder.vrc = intruder_self.vrc
    intruder.vsb = intruder_self.vsb
  end
end

#update coordination records with output
function update_to_coord!(coord::AbstractCASCoord, my_id::Int64, output::ACASXOutput)
  my_record = getRecord(coord, my_id)

  #update my record in coordination object
  #intruder-shared
  my_record.sensitivity_index = output.sensitivity_index
  #others (id,modes,equipage,quant) don't need to be updated...

  #intruder-specific
  for i = 1:endof(my_record.intruders)
    my_record.intruders[i].id           = output.intruders[i].id
    my_record.intruders[i].cvc          = output.intruders[i].cvc
    my_record.intruders[i].vrc          = output.intruders[i].vrc
    my_record.intruders[i].vsb          = output.intruders[i].vsb
  end
  setRecord(coord, my_id, my_record) #push back to coord

  return output
end

function initialize(cas)
  reset!(cas.input)
  reset!(cas.output)
  reset!(cas, getRecord(cas.coord, cas.my_id))
end

function reset!(cas, rec::ACASXCoordRecord)
  rec.my_id = cas.my_id
  rec.modes = cas.my_id
  rec.equipage = cas.equipage
  rec.quant = cas.quant
  rec.sensitivity_index = 0x0
  rec.protection_mode = 0x0

  for i = 1:endof(rec.intruders)
    rec.intruders[i].id = UInt32(i)
    rec.intruders[i].cvc = 0x0
    rec.intruders[i].vrc = 0x0
    rec.intruders[i].vsb = 0x0
  end
end

function quantize(x::AbstractFloat, b::AbstractFloat)
  # quantize x to the nearest multiple of b
  d, r = divrem(x, b)
  return b * (d + round(r / b))
end

end #module

max_stars_repo_name: sisl/RLES-SISLES

max_issues_repo_name: sisl/RLES-SISLES

max_forks_repo_name: sisl/RLES-SISLES


text:
module AdobeGlyphList

export  aglfn,
        agl,
        zapfdingbats

using DelimitedFiles
using Pkg

load_file(filename) = readdlm(filename, ';', String, '\n', comments=true, comment_char='#')

to_uint16(x) = parse(UInt16, x, base=16)

"""
```
    aglfn() -> Matrix{Any}
```
    Returns the mapping of glyph code to unicode character code for new fonts.
"""
function aglfn()
    path = joinpath(@__DIR__, "..", "aglfn.txt")
    a = load_file(path)
    b = similar(a, Any)
    b[:, 2:3] = a[:, 2:3]
    b[:, 1]   = map(x->Char(to_uint16(x)), a[:, 1])
    return b
end

"""
```
    agl() -> Matrix{Any}
```
    Returns the mapping of glyph code to unicode character code.
"""
function agl()
    path = joinpath(@__DIR__, "..", "glyphlist.txt")
    a = load_file(path)
    b = similar(a, Any)
    b[:, 1] = a[:, 1]
    b[:, 2] = map(a[:, 2]) do x
        y = split(x, ' ')[1]
        return Char(to_uint16(y))
    end
    return b
end

"""
```
    zapfdingbats() -> Matrix{Any}
```
    Returns the mapping of glyph code to unicode character code.
"""
function zapfdingbats()
    path = joinpath(@__DIR__, "..", "zapfdingbats.txt")
    a = load_file(path)
    b = similar(a, Any)
    b[:, 1] = a[:, 1]
    b[:, 2] = map(x->Char(to_uint16(x)), a[:, 2])
    return b
end

end


text:
using ..MultiInterface

export
    importance_sample,
    support_quality,
    support_quality_rank,
    support_quality_from_rank,
    __OptVec,
    __Opt

# These are the new operator definitions where the type signature is specified

"""__OptVec{T} = Union{Vector{Union{}}, Vector{T}}"""
__OptVec{T} = Union{Vector{Union{}}, Vector{T}}
"""__Opt{T} = Union{Nothing, T}"""
__Opt{T} = Union{Nothing, T}

# to support
MultiInterface.get_imp(::Nothing, args...) = nothing

@interface is_deterministic(sf::SFunc)::Bool
@interface sample(sf::SFunc{I,O}, i::I)::O where {I,O}
@interface sample_logcpdf(sf::SFunc{I,O}, i::I)::Tuple{O, AbstractFloat} where {I,O}
@interface invert(sf::SFunc{I,O}, o::O)::I where {I,O}
@interface lambda_msg(sf::SFunc{I,O}, i::SFunc{<:__Opt{Tuple{}}, O})::SFunc{<:__Opt{Tuple{}}, I} where {I,O}
@interface marginalize(sf::SFunc{I,O}, i::SFunc{<:__Opt{Tuple{}}, I})::SFunc{<:__Opt{Tuple{}}, O} where {I,O}
@interface logcpdf(sf::SFunc{I,O}, i::I, o::O)::AbstractFloat where {I,O}
@interface cpdf(sf::SFunc{I,O}, i::I, o::O)::AbstractFloat where {I,O}
@interface log_cond_prob_plus_c(sf::SFunc{I,O}, i::I, o::O)::AbstractFloat where {I,O}
@interface f_expectation(sf::SFunc{I,O}, i::I, fn::Function) where {I,O}
@interface expectation(sf::SFunc{I,O}, i::I)::O where {I,O}
@interface variance(sf::SFunc{I,O}, i::I)::O where {I,O}
@interface get_score(sf::SFunc{Tuple{I},O}, i::I)::AbstractFloat where {I,O}
@interface get_log_score(sf::SFunc{Tuple{I},O}, i::I)::AbstractFloat where {I,O}

@impl begin
    struct SFuncExpectation end

    function expectation(sf::SFunc{I,O}, i::I) where {I,O}
        return f_expectation(sf, i, x -> x)
    end
end

@interface support(sf::SFunc{I,O},
                    parranges::NTuple{N,Vector},
                    size::Integer,
                    curr::Vector{<:O}) where {I,O,N}

function importance_sample end

"""
    support_quality_rank(sq::Symbol)

Convert the support quality symbol into an integer for comparison.
"""
function support_quality_rank(sq::Symbol)
    if sq == :CompleteSupport return 3
    elseif sq == :IncrementalSupport return 2
    else return 1 end
end

"""
    support_quality_from_rank(rank::Int)

Convert the rank back into the support quality.
"""
function support_quality_from_rank(rank::Int)
    if rank == 3 return :CompleteSupport
    elseif rank == 2 return :IncrementalSupport
    else return :BestEffortSupport() end
end

@interface support_quality(sf::SFunc, parranges)

@impl begin
    struct SFuncSupportQuality end

    function support_quality(s::SFunc, parranges)
            :BestEffortSupport
    end
end

@interface bounded_probs(sf::SFunc{I,O},
                         range::__OptVec{<:O},
                         parranges::NTuple{N,Vector})::Tuple{Vector{<:AbstractFloat}, Vector{<:AbstractFloat}} where {I,O,N}

@interface make_factors(sf::SFunc{I,O},
                        range::__OptVec{<:O},
                        parranges::NTuple{N,Vector},
                        id,
                        parids::Tuple)::Tuple{Vector{<:Scruff.Utils.Factor}, Vector{<:Scruff.Utils.Factor}} where {I,O,N}

#= Statistics computation not included in the release
@interface initial_stats(sf::SFunc)

# TODO create an abstract type Stats{I,O}
# (range, parranges, pi's, lambda's)
@interface expected_stats(sf::SFunc{I,O},
                          range::__OptVec{<:O},
                          parranges::NTuple{N,Vector},
                          pis::NTuple{M,Dist},
                          child_lambda::Score{<:O}) where {I,O,N,M}

@interface accumulate_stats(sf::SFunc, existing_stats, new_stats)
@interface maximize_stats(sf::SFunc, stats)
=#
@interface compute_bel(sf::SFunc{I,O},
                      range::__OptVec{<:O},
                      pi::Dist{<:O},
                      lambda::Score)::Dist{<:O} where {I,O}

@interface compute_lambda(sf::SFunc,
                          range::__OptVec,
                          lambda_msgs::Vector{<:Score})::Score

@interface send_pi(sf::SFunc{I,O},
                   range::__OptVec{O},
                   bel::Dist{O},
                   lambda_msg::Score)::Dist{<:O} where {I,O}

@interface outgoing_pis(sf::SFunc,
                        range::__OptVec,
                        bel::Dist,
                        incoming_lambdas::__OptVec{<:Score})::Vector{<:Dist}

@interface outgoing_lambdas(sf::SFunc{I,O},
                  lambda::Score,
                  range::__OptVec,
                  parranges::NTuple{N,Vector},
                  incoming_pis::Tuple)::Vector{<:Score} where {N,I,O}

@interface compute_pi(sf::SFunc{I,O},
                     range::__OptVec{O},
                     parranges::NTuple{N,Vector},
                     incoming_pis::Tuple)::Dist{<:O} where {N,I,O}

@interface send_lambda(sf::SFunc{I,O},
                       lambda::Score,
                       range::__OptVec,
                       parranges::NTuple{N,Vector},
                       incoming_pis::Tuple,
                       parent_idx::Integer)::Score where {N,I,O}


text:
using ModelingToolkit, Test

@testset "Parsing Test" begin include("variable_parsing.jl") end
@testset "Differentiation Test" begin include("derivatives.jl") end
@testset "Simplify Test" begin include("simplify.jl") end
@testset "Direct Usage Test" begin include("direct.jl") end
@testset "System Construction Test" begin include("system_construction.jl") end
@testset "Distributed Test" begin include("distributed.jl") end


text:
#=
    pt="C:\\Users\\bernhard.konig\\Documents\\ASync\\home\\Code\\Julia\\DecisionTrees.jl"
    using Pkg
    Pkg.activate(pt) =#
import CSV
import DataFrames
import DataFrames: DataFrame

using DecisionTrees
using Profile
using Pkg

pkgdir = Pkg.dir("DecisionTrees") # how long will this remain supported?
testdir = joinpath(pkgdir, "test")
datadir = joinpath(pkgdir, "data")

# Read the data
datafile = joinpath(datadir, "freMTPL2\\freMTPL2.csv")
@assert isfile(datafile);
@time fullData = CSV.read(datafile, rows_for_type_detect=100000, allowmissing=:none, categorical=false);

areasSorted = sort(unique(fullData[:Area]))
AreaInteger = map(x->findall((in)([x]), areasSorted)[1], fullData[:Area])
fullData[:AreaInteger] = AreaInteger

# correct for unreasonable observations
for i = 1:size(fullData, 1)
    fullData[:ClaimNb][i] = min(4, fullData[:ClaimNb][i])
    fullData[:Exposure][i] = min(1.0, fullData[:Exposure][i])
end

# set independent variables
selected_explanatory_vars = ["Area","AreaInteger","VehPower","VehAge","DrivAge","BonusMalus","VehBrand","VehGas","Density","Region"]

# Prepare the data
##############################

dtmtable, sett, dfprepped = prepare_dataframe_for_dtm!(fullData, keycol="IDpol", trnvalcol="trnTest", numcol="ClaimNb", denomcol="Exposure", weightcol="Exposure", independent_vars=selected_explanatory_vars);

##############################
# Define model SETTINGS

updateSettingsMod!(sett,minWeight=-0.03,model_type="boosted_tree",iterations=2,learningRate=0.1,subsampling_features_prop=.7,calculatePoissonError=true)

##############################
# Run single Boosting Model
resultingFiles, resM = dtm(dtmtable, sett)

# measure allocation
updateSettingsMod!(sett,
iterations=40,
model_type="boosted_tree",
nScores="1000",
writeStatistics="false",
writeSasCode="false",
writeIterationMatrix="false",
writeResult="false",
writeCsharpCode="false",
writeVbaCode="false",
saveJLDFile="false",
saveResultAsJLDFile="false",

showProgressBar_time="false",
prroduceEstAndLeafMatrices="false",
write_dot_graph="false",

calculatePoissonError="true",
performanceMeasure="Average Poisson Error Val"
)

Profile.clear_malloc_data()
strs, resm2 = dtm(dtmtable, sett)
@info("allocation run done.")
quit();


text:
fname = joinpath(dirname(@__FILE__), "..", "data", "test.avr")
sname = joinpath(dirname(@__FILE__), "..", "data", "tmp", "same.avr")

a, b = read_avr(fname)

write_avr(sname, a, b, 8192)

a2, b2 = read_avr(sname)

@fact a --> a2
@fact b --> b2


text:
using Documenter, FIGlet, DocumenterMarkdown

cp(joinpath(@__DIR__, "../README.md"), joinpath(@__DIR__, "./src/index.md"), force=true, follow_symlinks=true)

makedocs(
         sitename="FIGlet.jl documentation",
         format = Markdown()
        )

deploydocs(;
    repo="github.com/kdheepak/FIGlet.jl",
    deps = Deps.pip(
                   "mkdocs==0.17.5",
                   "mkdocs-material==2.9.4",
                   "python-markdown-math",
                   "pygments",
                   "pymdown-extensions",
                   ),
    make = () -> run(`mkdocs build`),
    target = "site",
)


text:
using Dash, DashHtmlComponents, DashCoreComponents

app = dash()

dropdown_options = [
    Dict("label" => "New York City", "value" => "NYC"),
    Dict("label" => "Montreal", "value" => "MTL"),
    Dict("label" => "San Francisco", "value" => "SF"),
]
app.layout = html_div(style = Dict("columnCount" => 2)) do
    html_label("Dropdown"),
    dcc_dropdown(options = dropdown_options, value = "MTL"),
    html_label("Multi-Select Dropdown"),
    dcc_dropdown(
        options = dropdown_options,
        value = ["MTL", "SF"],
        multi = true,
    ),
    html_label("Radio Items"),
    dcc_radioitems(options = dropdown_options, value = "MTL"),
    html_label("Checkboxes"),
    dcc_checklist(options = dropdown_options, value = ["MTL", "SF"]),
    html_label("Text Input"),
    dcc_input(value = "MTL", type = "text"),
    html_label("Slider"),
    dcc_slider(
        min = 0,
        max = 9,
        marks = Dict([i => (i == 1 ? "Label $(i)" : "$(i)") for i = 1:6]),
        value = 5,
    )
end

run_server(app, "0.0.0.0", 8000, debug = true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
# we write a function that computes x factorial recursively
function recursiveFactorial(x)

	# we check that x is positive
	if x < 0
		error("factorials of non-negative numbers only")
	end

	# we check if x is either 0 or 1
	if x <= 1

		# we return the asnwer 1
		return 1
	else

		# we recurse
		return x * recursiveFactorial(x - 1)
	end
end


text:
## linear.jl : implementation of linear covariance function

## Linear ##
struct Linear{T} <: IsotropicCovarianceStructure{T}
    λ::T
    σ::T
    p::T

    function Linear{T}(λ::T, σ::T, p::T) where T
        λ > 0 || throw(DomainError(λ, "correlation length λ of linear covariance cannot be negative or zero"))
        σ > 0 || throw(DomainError(σ, "marginal standard deviation σ of linear covariance cannot be negative or zero"))
        p >= 1 || throw(DomainError(p, "in p-norm, p must be greater than or equal to 1"))
        isinf(p) && throw(DomainError(p, "in p-norm, p cannot be infinity"))

        new{T}(λ, σ, p)
    end
end

"""
    Linear(λ, [σ = 1], [p = 2])

Linear covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm, defined as

``C(x, y) = \\begin{cases} σ \\left(1 - \\displaystyle\\frac{ρ}{λ}\\right) & \\text{if }ρ ≤ λ\\\\ 0 & \\text{if }ρ>λ\\end{cases}``

with ``ρ = ||x - y||_p``.

# Examples
```jldoctest
julia> Linear(0.1)
linear (λ=0.1, σ=1.0, p=2.0)

julia> Linear(1.0, σ=2)
linear (λ=1.0, σ=2.0, p=2.0)

```
See also: [`Exponential`](@ref), [`Spherical`](@ref), [`Whittle`](@ref), [`Gaussian`](@ref), [`SquaredExponential`](@ref), [`Matern`](@ref)
"""
Linear(λ::Real; σ::Real=1.0, p::Real=2) = Linear{promote_type(typeof(λ),typeof(σ),typeof(p))}(promote(λ, σ, p)...)

# evaluate exponential covariance
apply(l::Linear, x::Real) = max(zero(x), 1 - x / l.λ)

# short name
shortname(::Linear) = "linear"


text:
using Krotov
using Documenter

# Generate examples
include("generate.jl")

DocMeta.setdocmeta!(Krotov, :DocTestSetup, :(using Krotov); recursive=true)

println("Starting makedocs")

makedocs(;
    modules=[Krotov],
    authors="Michael Goerz <mail@michaelgoerz.net> and contributors",
    repo="https://github.com/JuliaQuantumControl/Krotov.jl/blob/{commit}{path}#{line}",
    sitename="Krotov.jl",
    format=Documenter.HTML(;
        prettyurls=true,
        canonical="https://juliaquantumcontrol.github.io/Krotov.jl",
        assets=String[],
        mathengine=KaTeX(),
    ),
    pages=[
        "Home" => "index.md",
        "Overview" => "overview.md",
        "Examples" => [
            "List of Examples" => "examples/index.md",
            "Example 1 (TLS)" => "examples/simple_state_to_state.md",
            "Example 2 (RWA)" => "examples/state_to_state_rwa.md",
            "Example 3 (Diss. Gate)" => "examples/rho_3states.md",
            "Example 4 (Parametrization)" => "examples/state_to_state_parametrizations.md",
        ],
        "API" => "api.md",
        "History" => "history.md",
    ],
)

println("Finished makedocs")

rm(joinpath(@__DIR__, "build", "examples", ".gitignore"))

deploydocs(;
    repo="github.com/JuliaQuantumControl/Krotov.jl",
)


text:
"""Counts the number of each type of example in a dataset."""
function class_counts(rows::Data)::Dict{AbstractString, Int64}
    counts = Dict()  # a dictionary of label -> count.
    for row in rows
        # in our dataset format, the label is always the last column
        label = row[end]
        if !(label in keys(counts)) # ! binds to the label not the result of in expr
            counts[label] = 0
        end
        counts[label] += 1
    end
    counts
end


text:
using DataFrames

"""
List all columns that only consists of `NA`

    all_nacols(df::AbstractDataFrames)

"""
function all_nacols(df::AbstractDataFrame)
    nacols = Array{Symbol}([])
    for n in names(df)
        allna(df[n]) && push!(nacols, n)
    end
    nacols
end

"""
Drop all columns that only consists of `NA`

    drop_all_nacols!(df::AbstracDataFrame)

"""
function drop_all_nacols!(df::AbstractDataFrame)
    nacols = all_nacols(df)
    df = delete!(df, nacols)
    df
end


text:
# tests based on AMLCR 2nd ed,
@testset "AMLCR Math" begin

    # §5.7 payment frequency and certain (Tables 5.1 and 5.2)
    function f(age)
        life = SingleLife(mort = AMLCR_std, issue_age = age)
        int = Yields.Constant(0.05)
        lc = [
            AnnuityImmediate(life,int)
            AnnuityImmediate(life,int,frequency=4)
            AnnuityDue(life,int)
            AnnuityDue(life,int,frequency=4)
            AnnuityImmediate(life,int,n=10)
            AnnuityImmediate(life,int,n=10,frequency=4)
            AnnuityDue(life,int,n=10,frequency=4)
        ]

        return present_value.(lc)
    end

    # atol only to two decimals b/c AMLCR uses rounded rates?
    @test all(isapprox.(f(20), [18.966, 19.338, 19.966, 19.588, 7.711, 7.855, 7.952],atol=.01))
    @test all(isapprox.(f(80), [ 7.548,  7.917,  8.548,  8.167, 6.128, 6.373, 6.539],atol=.01))

    life = SingleLife(mort = AMLCR_std, issue_age = 20)
    int = Yields.Constant(0.05)

    # property tests

    @test present_value(AnnuityDue(life,int)) > present_value(AnnuityDue(life,int,start_time=5))

    @test present_value(AnnuityDue(life,int,n=5,certain=5,start_time=5)) == sum([1 for t in 5:9] .* [1/1.05 ^t for t in 5:9])

    # relation on pg 124
    @test issorted(present_value.([
        AnnuityImmediate(life,int)
        AnnuityImmediate(life,int,frequency=4)
        AnnuityDue(life,int,frequency=4)
        AnnuityDue(life,int)
    ]))

end#


text:
include("liftedfrombindeps.jl")

# define current version:
mskvmajor = "8"
mskvminor = "0"

mskplatform,distroext =
  if Sys.ARCH == :i386 || Sys.ARCH == :i686
    if     is_linux()   "linux32x86",  ".tar.bz2"
    elseif is_windows() "win32x86",  ".zip"
    else   error("Platform not supported")
    end
  elseif Sys.ARCH == :x86_64
    if     is_linux()   "linux64x86",".tar.bz2"
    elseif is_apple()   "osx64x86",  ".tar.bz2"
    elseif is_windows() "win64x86",  ".zip"
    else   error("Platform not supported")
    end
  else
    error("Platform not supported")
  end

libalternatives =
  if     mskplatform == "linux32x86"
                                       [ ("libmosek.so.8.0",     "libmosekscopt8_0.so"),    ]
  elseif mskplatform == "linux64x86"
                                       [ ("libmosek64.so.8.0",   "libmosekscopt8_0.so"),    ]
  elseif mskplatform == "osx64x86"
                                       [ ("libmosek64.8.0.dylib","libmosekscopt8_0.dylib"), ]
  elseif mskplatform == "win32x86"
                                       [ ("mosek8_0.dll",        "mosekscopt8_0.dll"),      ]
  elseif mskplatform == "win64x86"
                                       [ ("mosek64_8_0.dll",     "mosekscopt8_0.dll"),      ]
  else   error("Platform not supported")
  end

bindepsdir = dirname(@__FILE__)

function versionFromBindir(bindir ::String)
    try
        mosekbin = if is_windows() "mosek.exe" else "mosek" end
        txt = readstring(`$bindir/$mosekbin`)
        m = match(r"\s*MOSEK Version ([0-9]+\.[0-9]+\.[0-9]+\.[0-9]+)",txt)
        if m == nothing
            return nothing
        else
            return m.captures[1]
        end
    catch
        return nothing
    end
end

function bindirIsCurrentVersion(bindir)
    ver = versionFromBindir(bindir)
    return ver != nothing && ver[1] == mskvmajor && ver[2] == mskvminor
end

# Detect previous installation method:
#   "internal" -> the MOSEK distro was downloaded and installed by the installer
#   "external" -> the MOSEK distro was detected in some other location and not downloaded by the installer
# If the previous installation was internal, we will not attempt to locate a MOSEK installation
#
# This ensures that you can rebuild an internal installation without
# the installer checking for other installations.
instmethod =
    try
        open(joinpath(bindepsdir,"inst_method"),"r") do f
            strip(readstring(f))
        end
    catch
        nothing
    end

usepreinstalled = ! haskey(ENV,"MOSEKJL_FORCE_DOWNLOAD") && instmethod != "internal"

mskbindir =
# 1. Is MOSEKBINDIR set? If so this must point to the binaries dir in the MOSEK DISTRO
    if  usepreinstalled &&
        haskey(ENV,"MOSEKBINDIR")  &&
        bindirIsCurrentVersion(ENV["MOSEKBINDIR"])

        instmethod = "external"

        ENV["MOSEKBINDIR"]
    elseif ! usepreinstalled && haskey(ENV,"MOSEK_8_0_BINDIR") && bindirIsCurrentVersion(ENV["MOSEKBINDIR"])
        instmethod = "external"
        ENV["MOSEK_8_0_BINDIR"]
# 2a. Otherwise, use the default installation path (Linux)
    elseif usepreinstalled &&
        haskey(ENV,"HOME") &&
        bindirIsCurrentVersion(joinpath(ENV["HOME"],"mosek","8","tools","platform",mskplatform,"bin"))

        instmethod = "external"

        joinpath(ENV["HOME"],"mosek","8","tools","platform",mskplatform,"bin")
# 2b. Windows default install path
    elseif usepreinstalled &&
        haskey(ENV,"HOMEDRIVE") &&
        haskey(ENV,"HOMEPATH") &&
        bindirIsCurrentVersion(joinpath(string(ENV["HOMEDRIVE"],ENV["HOMEPATH"]),"mosek","8","tools","platform",mskplatform,"bin"))

        home = string(ENV["HOMEDRIVE"],ENV["HOMEPATH"])

        instmethod = "external"

        joinpath(home,"mosek","8","tools","platform",mskplatform,"bin")
# 3. Otherwise, fetch the MOSEK distro and unpack it
    else
        srcdir   = joinpath(bindepsdir,"src")
        dldir    = joinpath(bindepsdir,"downloads")
        archname = "mosektools$(mskplatform)$(distroext)"
        verurl   = "http://download.mosek.com/stable/$mskvmajor/version"

        mkpath(dldir)

        cur_version =
            if isfile(joinpath(bindepsdir,"version"))
                open(joinpath(bindepsdir,"version"),"r") do f
                    cur_version = strip(readstring(f))
                end
            else
                nothing
            end

        info("Get latest MOSEK version (http://download.mosek.com/stable/version)")
        success(download_cmd(verurl, joinpath(dldir,"new_version"))) || error("Failed to get MOSEK version")

        new_version =
            open(joinpath(dldir,"new_version"),"r") do f
                strip(readstring(f))
            end
        info("Latest MOSEK version = $new_version, currently installed = $cur_version")

        instmethod = "internal"

        if cur_version == nothing || cur_version != new_version
            archurl = "http://download.mosek.com/stable/$(new_version)/$archname"
            info("Download MOSEK distro ($archurl)")

            basename,ext,sndext = splittarpath(archname)

            success(download_cmd(archurl, joinpath(dldir,archname))) || error("Failed to download MOSEK distro")

            mkpath(srcdir)
            info("Unpack MOSEK distro ($dldir/$archname -> $srcdir)")
            success(unpack_cmd(joinpath(dldir,archname),srcdir, ext, sndext)) || error("Failed to unpack MOSEK distro")

            joinpath(srcdir,"mosek",mskvmajor,"tools","platform",mskplatform,"bin")
        else
            info("Update not necessary")
            joinpath(srcdir,"mosek",mskvmajor,"tools","platform",mskplatform,"bin")
        end
    end

mskbindir = replace(mskbindir,"\\","/")

idxs = reverse(find(libs -> all(lib -> isfile(joinpath(mskbindir,lib)), libs), libalternatives))
if length(idxs) == 0
    error("Failed to find any usable MOSEK libraries")
else
    version = versionFromBindir(mskbindir)

    if version == nothing
        error("MOSEK package is broken")
    else
        open(joinpath(bindepsdir,"version"),"w") do f
            write(f,version)
        end
    end

    if instmethod == "external"
        info("""Found MOSEK $version at $mskbindir""")
    end

    open(joinpath(bindepsdir,"inst_method"),"w") do f
        write(f,instmethod)
    end

    libmosek,libmosekscopt = libalternatives[idxs[1]]

    libmosekpath = escape_string(normpath("$mskbindir/$libmosek"))
    libscoptpath = escape_string(normpath("$mskbindir/$libmosekscopt"))

    open(joinpath(bindepsdir,"deps.jl"),"w") do f
        write(f,"""# This is an auto-generated file; do not edit\n""")
        write(f,"""# Macro to load a library\n""")
        write(f,"""macro checked_lib(libname, path)\n""")
        write(f,"""    (Libdl.dlopen_e(path) == C_NULL) && error("Unable to load \\n\\n\$libname (\$path)\\n\\nPlease re-run Pkg.build(package), and restart Julia.")\n""")
        write(f,"""    quote const \$(esc(libname)) = \$path end\n""")
        write(f,"""end\n""")

        write(f,"""\n""")
        write(f,"""# Load dependencies\n""")
        write(f,string("""@checked_lib libmosek      \"$libmosekpath\"\n"""))
        write(f,string("""@checked_lib libmosekscopt \"$libscoptpath\"\n"""))
    end
end


text:
mutable struct Env
    ptr::Ptr{Cvoid}
    num_models::Int
    finalize_called::Bool

    function Env()
      stat = Vector{Cint}(undef, 1)
      tmp = @cpx_ccall(openCPLEX, Ptr{Cvoid}, (Ptr{Cint},), stat)
      if tmp == C_NULL
          error("CPLEX: Error creating environment")
      end
      env = new(tmp, 0, false)
      function env_finalizer(env)
          if env.num_models == 0
              close_CPLEX(env)
          else
              env.finalize_called = true
          end
      end
      @compat finalizer(env_finalizer, env)
      return env
    end
end

convert(ty::Type{Ptr{Cvoid}}, env::Env) = env.ptr::Ptr{Cvoid}
unsafe_convert(ty::Type{Ptr{Cvoid}}, env::Env) = convert(ty, env)

function is_valid(env::Env)
    env.ptr != C_NULL
end

function notify_new_model(env::Env)
    env.num_models += 1
end

function notify_freed_model(env::Env)
    @assert env.num_models > 0
    env.num_models -= 1
    if env.num_models <= 0 && env.finalize_called
        close_CPLEX(env)
    end
end

function close_CPLEX(env::Env)
    tmp = Ptr{Cvoid}[env.ptr]
    stat = @cpx_ccall(closeCPLEX, Cint, (Ptr{Cvoid},), tmp)
    env.ptr = C_NULL
    if stat != 0
        throw(CplexError(env, stat))
    end
end

function set_logfile(env::Env, filename::String)
  @assert isascii(filename)
  fp = @cpx_ccall(fopen, Ptr{Cvoid}, (Ptr{Cchar}, Ptr{Cchar}), filename, "w")
  if fp == C_NULL
    error("CPLEX: Error setting logfile")
  end
  stat = @cpx_ccall(setlogfile, Cint, (Ptr{Cvoid}, Ptr{Cvoid}), env, fp)
  if stat != 0
    throw(CplexError(env, stat))
  end
end

function get_error_msg(env::Env, code::Number)
    @assert env.ptr != C_NULL
    buf = Vector{Cchar}(undef, 4096) # minimum size for Cplex to accept
    errstr = @cpx_ccall(geterrorstring, Ptr{Cchar}, (Ptr{Cvoid}, Cint, Ptr{Cchar}), env.ptr, convert(Cint, code), buf)
    if errstr != C_NULL
      return unsafe_string(pointer(buf))
    else
      error("CPLEX: error getting error message(!)")
    end
end

function version(env::Env = Env())
    charptr = @cpx_ccall(version, Ptr{Cchar}, (Ptr{Cvoid},), env.ptr)
    if charptr != C_NULL
        return unsafe_string(charptr)
    else
        return error("CPLEX: error getting version")
    end
end

struct CplexError <: Exception
  code::Int
  msg::String

  function CplexError(env::Env, code::Integer)
    new(convert(Cint, code), get_error_msg(env, code))
  end
end


text:
## exponential.jl : implementation of exponential covariance function

## Exponential ##
struct Exponential{T} <: IsotropicCovarianceStructure{T}
    λ::T
    σ::T
    p::T

    function Exponential{T}(λ::T, σ::T, p::T) where T
        λ > 0 || throw(DomainError(λ, "correlation length λ of exponential covariance cannot be negative or zero"))
        σ > 0 || throw(DomainError(σ, "marginal standard deviation σ of exponential covariance cannot be negative or zero"))
        p >= 1 || throw(DomainError(p, "in p-norm, p must be greater than or equal to 1"))
        isinf(p) && throw(DomainError(p, "in p-norm, p cannot be infinity"))

        new{T}(λ, σ, p)
    end
end

"""
    Exponential(λ, [σ = 1], [p = 2])

Exponential covariance structure with correlation length `λ`, (optional) marginal standard deviation `σ` and (optional) `p`-norm defined as

``C(x, y) = σ \\exp\\left(-\\displaystyle\\frac{ρ}{λ}\\right)``

with ``ρ = ||x - y||_p``.

# Examples
```jldoctest
julia> Exponential(0.1)
exponential (λ=0.1, σ=1.0, p=2.0)

julia> Exponential(1.0, σ=2)
exponential (λ=1.0, σ=2.0, p=2.0)

```
See also: [`Linear`](@ref), [`Spherical`](@ref), [`Whittle`](@ref), [`Gaussian`](@ref), [`SquaredExponential`](@ref), [`Matern`](@ref)
"""
Exponential(λ::Real; σ::Real=1.0, p::Real=2) = Exponential{promote_type(typeof(λ),typeof(σ),typeof(p))}(promote(λ, σ, p)...)

# evaluate exponential covariance
apply(e::Exponential, x::Real) = exp(-x / e.λ)

# short name
shortname(::Exponential) = "exponential"


text:
using OPCSPs
using OPCSPs.MVNTools

r = [0, 10.1, 5, 5, 5, 0]
cov = Float64[ 0 0 0 0 0 0; 0 0 0 0 0 0; 0 0 2 2 -2 0; 0 0 2 4 0 0; 0 0 -2 0 4 0; 0 0 0 0 0 0]

dist = MVN(r,cov)

# BAD NOTATION delta is not what it usually is
M = 1000000
rewards = Array(Float64, M)
rng = MersenneTwister(1)
for i = 1:M
    delta = rand(dist)
    if delta[3] >= 5.0
        rewards[i] = delta[3] + delta[4]
    else
        rewards[i] = delta[3] + delta[5]
    end
end

#=
using Gadfly
p = plot(rewards, Geom.histogram)
draw(SVG("not_gaussian.svg", 6inch, 3inch), p)
=#

import PyPlot
PyPlot.plt[:clf]()
PyPlot.plt[:hist](rewards, 100)

# not gaussian


text:
@testset "Non-negative CANDECOMP" begin
Random.seed!(12345)

r = 2
T = rand_kruskal3(r, (10, 20, 30), true)

@time factors = nncp(T, r, compute_error=true)
@test length(factors.factors) == ndims(T)
@test size.(factors.factors) == ntuple(i -> (size(T, i), r), 3)
@test rank(factors) == r
@test rel_residue(factors) < 0.05

end


text:
using Test

@testset "all" begin
    @testset "FED" begin include("FED.jl") end
    @testset "FEDW" begin include("FEDW.jl") end
    @testset "Perturbation" begin include("Perturbation.jl") end
    #@testset "RPA" begin include("RPA.jl") end
end


text:
using MittagLefflerFunctions
import MittagLefflerFunctions.Quadrature.MLQuad
using RationalApproximations
using PyPlot
using Printf

G(t) = (1-t)/(1+t)
α_vals = collect(1:5) ./ 5
Eα = [ MLQuad(α, 1.0, 14, :hyperbola) for α in α_vals ]
fα = Vector{Function}(undef, 5)
for col = 1:5
    fα[col] = ( t -> Eα[col](-G(t)) )
end
t = range(-1, 1, length=800)
max_m = 10
tol = 1e-15

rows = 10
cols = length(α_vals)
errors = fill(NaN, rows, cols)
for col = 1:cols
    α = α_vals[col]
    println("α = ", α)
    F = fα[col].(t)
    r, err = aaa_v2(F, t, max_m, tol)
    display(err)
    M = length(r)
    for row = 1:M
        errors[row,col] = err[row]
    end
end

@printf("%5s ", "m")
for col = 1:cols
    @printf("&    α = %3.1f", α_vals[col])
end
@printf("\n\n")

for row = 1:rows
    @printf("%5d ", row)
    for col = 1:cols
        @printf("& %10.2e", errors[row,col])
    end
    @printf("\\\\\n")
end


text:
#!/usr/bin/env julia

# Best run in the REPL until I work out how to get unicodeplots to print to stdout when in an `include`
# `env JULIA_NUM_THREADS=(nproc) julia --project`

using JustJoshing

import Plots
Plots.unicodeplots()

# This should look like Figure 5.3, page 121 in Joshi (it does)
# NB: looks like Joshi's graph is mislabelled - should be spot price
p = Plots.plot(); for t in 0:0.2499:1; Plots.plot!(p,x->C(x,t;K=100,σ=0.3,r=0),60:140); end; p

Plots.plot(x->C(1,x),0:0.001:0.9999) # Value of at-the-money approaches zero as time approaches maturity (Fig 2.1, page 35)

Plots.plot(x->C(0.5,0;σ=x),0:0.01:1) # Call options with volatile underlyings are more expensive (Fig 3.8, page 66)

# Simulated stock price
Plots.plot(x->B(x,σ=0.02),0:0.01:1)

# Monte-Carlo validation of contracts from payoffs

# max(S-K,0): call option
let S_t=100, t=0.01, K=100, r=0.02, σ=0.05, T=1, trials=100_000_000
    @show bs = C(S_t,t;T=T,K=K,r=r,σ=σ)

    @show mc = mc_pricer((S,K,t,T,r)->max(S-K,0),S_t;t=t,T=T,K=S_t,r=r,σ=σ) |> first

    @assert isapprox(mc, bs, rtol=1e-3)
end

# S-K: forward contract
let S_t=100, t=0.01, K=100, r=0.02, σ=0.05, T=1, trials=100_000_000
    @show exact = S_t - K*exp(-r*(T-t))

    @show mc = mc_pricer((S,K,t,T,r)->S-K,S_t;t=t,T=T,K=K,r=r,σ=σ) |> first

    @assert isapprox(mc, exact, rtol=1e-3)
end

# Int(S>E): binary option
let S_t=100, t=0.01, E=100, r=0.02, σ=0.05, T=1, trials=100_000_000
    @show exact = binary(S_t,t;E=E,r=r,σ=σ,T=T)

    @show mc = mc_pricer((S,E,t,T,r)->Int(S>E),S_t;t=t,T=T,K=E,r=r,σ=σ) |> first

    @assert isapprox(mc, exact, rtol=1e-3)
end

# Asian option, fixed strike
# TODO: compare prices with European options, look at impact volatility has
s = 80:1:120; t = [0; round.(-1*10 .^(-2:0.2:-1);sigdigits = 2)]
a = (x->mc_pricer_pathdep((S,K,t,T,r)->max(mean(S)-100,0),x[1],trials=100_000,r=0,t=x[2]:0.1:0,T=0,σ=0.5)).(Iterators.product(s,t))
plot(s,first.(a),labels=string.(t'))


text:
######################################################################
# plot_trajectory.jl
#
# Contains files for plotting the trajectory
######################################################################

# what other stuff do we need here?
# only marks, colors, etc
function plot_trajectory(xd::VVF;
                         c::String="b",
                         ls::String="-",
                         lw::Real=1.5,
                         mew::Real=1,
                         mfc::String="w",
                         ms::Real=10,
                         dims::Int=2
                        )

    N = length(xd)
    xvals = zeros(N)
    yvals = zeros(N)
    zvals = zeros(N)
    for i = 1:N
        xvals[i] = xd[i][1]
        yvals[i] = xd[i][2]
    end

	# if in 3 dimensions
    if dims == 3
        for i = 1:N
            zvals[i] = xd[i][3]
        end
        plot3D(xvals, yvals, zvals, ls=ls, marker=m, ms=ms, lw=lw)
        xlabel(L"x")
        ylabel(L"y")
        zlabel(L"z")
        return
    end

    m = "."
    plot(xvals, yvals, c, linestyle=ls, marker=m, lw=lw, ms=ms, mfc=mfc)

end

# plotting for multiple agents
# vtm is a vector of trajectory managers
function plot_trajectory(xd::VVF, vtm::VTM;
                         lw::Real=1,
                         ms::Real=9,
                         dims::Int=2
                        )

    num_agents = length(vtm)
    xds = vvf2vvvf(xd, vtm)
    for j = 1:num_agents
        plot_trajectory(xds[j], lw=lw, ms=ms,onlyMarks=onlyMarks,dims=dims)
    end
end