text:
# --------------------------------------------------------------------------
# ACE1.jl: Julia implementation of the Atomic Cluster Expansion
# Copyright (c) 2019 Christoph Ortner <christophortner0@gmail.com>
# Licensed under ASL - see ASL.md for terms and conditions.
# --------------------------------------------------------------------------

module DAG

import Base: ==, length

using Combinatorics: combinations, partitions

const BinDagNode{TI} = Tuple{TI,TI}

struct CorrEvalGraph{T, TI}
   nodes::Vector{BinDagNode{TI}}
   vals::Vector{T}
   num1::TI
   numstore::TI
end

length(dag::CorrEvalGraph) = length(dag.nodes)

==(dag1::CorrEvalGraph, dag2::CorrEvalGraph) = ACE1._allfieldsequal(dag1, dag2)

CorrEvalGraph{T, TI}() where {T, TI} =
      CorrEvalGraph(Vector{BinDagNode{TI}}(undef, 0),
                    Vector{T}(undef, 0),
                    0, 0)

# -------------- FIO

# TODO: maybe there is a cleverer way to do this, for now this is just
#       a quick hack to make sure it can be read without any ambiguity

write_dict(gr::CorrEvalGraph{T, TI}) where {T <: Real, TI <: Integer} =
   Dict( "__id__" => "ACE1_CorrEvalGraph",
         "T" => write_dict(T),
         "TI" => write_dict(TI),
         "nodes1" => [ n[1] for n in gr.nodes ],
         "nodes2" => [ n[2] for n in gr.nodes ],
         "vals" => gr.vals,
         "num1" => gr.num1,
         "numstore" => gr.numstore
      )

function read_dict(::Val{:ACE1_CorrEvalGraph}, D::Dict)
   T = read_dict(D["T"])
   TI = read_dict(D["TI"])
   @assert T <: Real
   @assert TI <: Integer
   return CorrEvalGraph{T, TI}(
      collect(zip(D["nodes1"], D["nodes2"])),
      D["vals"],
      D["num1"],
      D["numstore"]
   )
end

# ---------------------------------------------------------------------
#   partition generator

_score_partition(p) = isempty(p) ? Inf : (1e9 * length(p) + maximum(p))

function _get_ns(p, specnew, specnew_dict)
   out = Vector{Int}(undef, length(p))
   for (i, kk_) in enumerate(p)
      if haskey(specnew_dict, kk_)
         out[i] = specnew_dict[kk_]
      else
         return Int[]
      end
   end
   return out
end

function _find_partition(kk, specnew, specnew_dict)
   worstp = _get_ns([ [k] for k in kk ], specnew, specnew_dict)
   @assert worstp == kk
   bestp = worstp
   bestscore = _score_partition(bestp)

   for ip in partitions(1:length(kk))
      p = _get_ns([ kk[i] for i in ip ], specnew, specnew_dict)
      # p_old = _get_ns_old([ kk[i] for i in ip ], specnew, specnew_dict)
      # if p != p_old
      #    @infiltrate
      # end
      score = _score_partition(p)
      if !isempty(p) && score < bestscore
         bestp = p
         bestscore = score
      end
   end

   return bestp
end

# function _insert_node!(nodes, coeffsnew, specnew, specnew_dict,
#                        newnode)

# return value is the number of fake nodes added to the dag
function _insert_partition!(nodes, coeffsnew, specnew, specnew_dict,
                            kk, p,
                            ikk, coeffsN, specN,
                            TI = Int)
   if length(p) == 2
      newnode = BinDagNode{TI}((p[1], p[2]))
      # _insert_node!(nodes, coeffsnew, specnew, specnew_dict, newnode)
      push!(nodes, newnode)
      push!(coeffsnew, coeffsN[ikk])
      push!(specnew, kk)
      specnew_dict[kk] = length(specnew)
      return 0
   else
      # @show kk, p
      # @infiltrate
      # reduce the partition by pushing a new node
      push!(nodes, BinDagNode{TI}((p[1], p[2])))
      push!(coeffsnew, 0)
      kk1 = sort(vcat(specnew[p[1]], specnew[p[2]]))
      push!(specnew, kk1)
      specnew_dict[kk1] = length(specnew)
      # and now recurse with the reduced partition
      return 1 + _insert_partition!(nodes, coeffsnew, specnew, specnew_dict,
                         kk, vcat( [length(nodes)], p[3:end] ),
                         ikk, coeffsN, specN, TI)
   end
end

function get_eval_graph(inner, coeffs;   # inner::InnerPIBasis
                        filter = _->true,
                        TI = Int,
                        verbose = false)
   # make a list of all basis functions as vectors so we can search it
   # TODO: should also check the tuples are sorted lexicographically
   spec = [ inner.iAA2iA[iAA, 1:inner.orders[iAA]]
            for iAA = 1:length(inner) ]
   @assert issorted(length.(spec))
   @assert all(issorted, spec)
   # we need to separate them into 1-p and many-p
   spec1 = spec[ length.(spec) .== 1 ]
   coeffs1 = coeffs[1:length(spec1)]
   IN = (length(spec1)+1):length(spec)
   specN = spec[IN]
   coeffsN = coeffs[IN]

   # start assembling the dag
   nodes = BinDagNode{TI}[]
   sizehint!(nodes, length(inner))
   coeffsnew = Vector{eltype(coeffs)}()
   sizehint!(coeffsnew, length(inner))
   specnew = Vector{Int}[]
   specnew_dict = Dict{Vector{Int}, Int}()
   sizehint!(specnew, length(inner))

   # add the full 1-particle basis (N=1) into the dag
   num1 = maximum(inner.iAA2iA)
   for i = 1:num1
      push!(nodes, BinDagNode{TI}((i, 0)))
      push!(specnew, [i])
      specnew_dict[ [i] ] = length(specnew)
      # find that index in `spec`
      ispec = findfirst(isequal([i]), spec1)
      if isnothing(ispec)
         push!(coeffsnew, 0)
      else
         push!(coeffsnew, coeffs1[ispec])
      end
   end

   # now we can construct the rest
   extranodes = 0
   for (ikk, kk) in enumerate(specN)
      # find a good partition of kk
      p = _find_partition(kk, specnew, specnew_dict)
      extranodes += _insert_partition!(nodes, coeffsnew, specnew, specnew_dict,
                                       kk, p, ikk, coeffsN, specN, TI)
   end

   verbose && @info("Extra nodes inserted into the dag: $extranodes")
   numstore = length(nodes)
   num1old = num1

   # re-organise the dag layout to minimise numstore
   nodesfinal, coeffsfinal, num1, numstore = _reorder_dag!(nodes, coeffsnew)

   return CorrEvalGraph(nodesfinal, coeffsfinal, TI(num1), TI(numstore))
end

function _reorder_dag!(nodes::Vector{BinDagNode{TI}}, coeffs::Vector{T}) where {TI, T}
   # collect all AA indices that are used anywhere in the dag
   newinds = zeros(Int, length(nodes))
   newnodes = BinDagNode{TI}[]
   newcoeffs = T[]

   # inds2 = stage-2 indices, i.e. temporary storage
   # inds3 = stage-3 indices, i.e. no intermediate storage
   inds2 = sort(unique([[ n[1] for n in nodes ]; [n[2] for n in nodes]]))

   # first add all 1p nodes
   for i = 1:length(nodes)
      n, c = nodes[i], coeffs[i]
      if (n[2] == 0) ## && ((c != 0) || (i in inds2))
         @assert n[1] == i
         newinds[i] = i
         push!(newnodes, n)
         push!(newcoeffs, c)
      end
   end
   num1 = length(newnodes)

   # next add the remaining dependent nodes
   for i = 1:length(nodes)
      n, c = nodes[i], coeffs[i]
      # not 1p basis && dependent node
      if (n[2] != 0) && (i in inds2)
         push!(newnodes, BinDagNode{TI}((newinds[n[1]], newinds[n[2]])))
         push!(newcoeffs, c)
         newinds[i] = length(newnodes)
      end
   end
   numstore = length(newnodes)

   # now go through one more time and add the independent nodes
   for i = 1:length(nodes)
      n, c = nodes[i], coeffs[i]
      if (n[2] != 0) && (newinds[i] == 0) && (c != 0)
         push!(newnodes, BinDagNode{TI}((newinds[n[1]], newinds[n[2]])))
         push!(newcoeffs, c)
         newinds[i] = length(newnodes)
      end
   end

   return newnodes, newcoeffs, num1, numstore
end

# ---------------------------------------------------------------------
#   evaluation codes

function traverse_fwd!(AA, dag::CorrEvalGraph, A, fun)
   nodes = dag.nodes
   vals = dag.vals
   @assert length(AA) >= dag.numstore
   @assert length(A) >= dag.num1

   # Stage-1: copy the 1-particle basis into AA
   @inbounds for i = 1:dag.num1
      AA[i] = a = A[i]
      fun(vals[i], a)
   end

   # Stage-2: go through the dag and store the intermediate results we need
   @inbounds for i = (dag.num1+1):dag.numstore
      n1, n2 = nodes[i]
      AA[i] = a = AA[n1] * AA[n2]
      fun(vals[i], a)
   end

   # Stage 3:
   # continue going through the dag, but now we don't need to store the new
   # correlations since the later expressions don't depend on them
   @inbounds for i = (dag.numstore+1):length(dag)
      n1, n2 = nodes[i]
      a = AA[n1] * AA[n2]
      fun(vals[i], a)
   end

   return nothing
end

# function traverse_fwdbwd!(AA, dag::CorrEvalGraph, A, fwdfun, bwdfun)
#    # function fwdbwd_traverse!(dEs, tmpd, V::GraphPIPot{T}, ) where {T}
#    nodes = dag.nodes
#    vals = dag.vals
#
#    # FORWARD PASS
#    # --------------
#
#    for i = 1:dag.num1
#       AA[i] = a = A[i]
#       fwdfun(vals[i], a)
#    end
#
#    # go through the dag and store the intermediate results we need
#    for i = (dag.num1+1):dag.numstore
#       n1, n2 = nodes[i]
#       AA[i] = a = muladd(AA[n1], AA[n2], AA[i])
#       fwdfun(vals[i], a)
#    end
#
#    # go through the dag and store the intermediate results we need
#    for i = (dag.num1+1):dag.numstore
#       n1, n2 = nodes[i]
#       AA[i] = a = muladd(AA[n1], AA[n2], AA[i])
#       fwdfun(vals[i], a)
#    end
#
#
#    # BACKWARD PASS
#    # --------------
#    # fill the B array -> coefficients of the derivatives
#    #  AA_i = AA_{n1} * AA_{n2}
#    #  ∂AA_i = AA_{n1} * ∂AA_{n2} + AA_{n1} * AA_{n2}
#    #  c_{n1} * ∂AA_{n1} <- (c_{n1} + c_i AA_{n2}) ∂AA_{n1}
#    @inbounds @fastmath for i = length(dag):-1:(dag.numstore+1)
#       c = coeffs[i]
#       n1, n2 = nodes[i]
#       B[n1] = muladd(c, AA[n2], B[n1])
#       B[n2] = muladd(c, AA[n1], B[n2])
#    end
#    # in stage 2 c = C[i] is replaced with b = B[i]
#    @inbounds @fastmath for i = dag.numstore:-1:(dag.num1+1)
#       n1, n2 = nodes[i]
#       b = B[i]
#       B[n1] = muladd(b, AA[n2], B[n1])
#       B[n2] = muladd(b, AA[n1], B[n2])
#    end
#
#    # stage 3: get the gradients
#    fill!(dEs, zero(JVec{T}))
#    A = tmpd.A
#    dA = tmpd.dA
#    for (iR, (R, Z)) in enumerate(zip(Rs, Zs))
#       evaluate_d!(A, dA, tmpd_basis1p, basis1p, R, Z, z0)
#       iz = z2i(basis1p, Z)
#       inds = basis1p.Aindices[iz, iz0]
#       for iA = 1:length(basis1p, iz, iz0)
#          dEs[iR] += real(B[inds[iA]] * dA[inds[iA]])
#       end
#    end
#
#    return dEs
# end

end


text:
using BinaryProvider # requires BinaryProvider 0.3.0 or later

# Parse some basic command-line arguments
const verbose = "--verbose" in ARGS
const prefix = Prefix(get([a for a in ARGS if a != "--verbose"], 1, joinpath(@__DIR__, "usr")))
products = [
    LibraryProduct(prefix, String["libbkmaxflow"], :libbkmaxflow),
]

# Download binaries from hosted location
bin_prefix = "https://github.com/Gnimuc/BKMaxflowBuilder/releases/download/v3.0.1.1"

# Listing of files generated by BinaryBuilder:
download_info = Dict(
    Linux(:aarch64, :glibc) => ("$bin_prefix/BKMaxflow.v3.0.1.aarch64-linux-gnu.tar.gz", "4e34f1164d2ab6c7bfb691e6b2fb5e0cdf5f92fc61062cf3c6ac60e6e43fb8a9"),
    Linux(:aarch64, :musl) => ("$bin_prefix/BKMaxflow.v3.0.1.aarch64-linux-musl.tar.gz", "699c182be216b665310404b9b76c0f41572b613eb8110c7c4ad8848cdda62e7d"),
    Linux(:armv7l, :glibc, :eabihf) => ("$bin_prefix/BKMaxflow.v3.0.1.arm-linux-gnueabihf.tar.gz", "2405eacd8ae8869e2b274fdc25e9812d997c03ae7dd720b92a8f89d290e50c75"),
    Linux(:armv7l, :musl, :eabihf) => ("$bin_prefix/BKMaxflow.v3.0.1.arm-linux-musleabihf.tar.gz", "df5103480a48474e6a12ba71ca37b577f19952c5175ba8564c9562f9d15dd39b"),
    Linux(:i686, :glibc) => ("$bin_prefix/BKMaxflow.v3.0.1.i686-linux-gnu.tar.gz", "a44e6fcc583fe1997bba6c8229472449548844575a3e90bdd6c6e8b9ddb33637"),
    Linux(:i686, :musl) => ("$bin_prefix/BKMaxflow.v3.0.1.i686-linux-musl.tar.gz", "4c80c9a4044d3a496f8a333696be21c5ec89a87f6c4df1f84d36c9e2fd86bce8"),
    Windows(:i686) => ("$bin_prefix/BKMaxflow.v3.0.1.i686-w64-mingw32.tar.gz", "dc2ac1539661c7bbbce62ac9e3bbfb5fc1622c77bb439f7acd29aa4510a99ac0"),
    Linux(:powerpc64le, :glibc) => ("$bin_prefix/BKMaxflow.v3.0.1.powerpc64le-linux-gnu.tar.gz", "80011beb757ce1b2c37b111418b01c7f74f84d809da0acdb8909bd43b2855241"),
    MacOS(:x86_64) => ("$bin_prefix/BKMaxflow.v3.0.1.x86_64-apple-darwin14.tar.gz", "5a150601288704665cdac03cec1e59d8f38bc50616099743afdb405ae3e93fb3"),
    Linux(:x86_64, :glibc) => ("$bin_prefix/BKMaxflow.v3.0.1.x86_64-linux-gnu.tar.gz", "3f4ef5f7e99be9d6b736d8e0f2d9f7cdefcf3837e78d6e91efcc430087da6300"),
    Linux(:x86_64, :musl) => ("$bin_prefix/BKMaxflow.v3.0.1.x86_64-linux-musl.tar.gz", "c6371cc21a2926e62157eb51b683ce534035c0b58fefec7c185b290a6db1e439"),
    FreeBSD(:x86_64) => ("$bin_prefix/BKMaxflow.v3.0.1.x86_64-unknown-freebsd11.1.tar.gz", "f67d5f82a5aebc4802fea33738e15b1b60d21da4f60dceb9f35ab3ff74c975b1"),
    Windows(:x86_64) => ("$bin_prefix/BKMaxflow.v3.0.1.x86_64-w64-mingw32.tar.gz", "f36f9e187cafe396b8180fafa24c262ad4647dc4b6dffc212678c819d60ecdce"),
)

# Install unsatisfied or updated dependencies:
unsatisfied = any(!satisfied(p; verbose=verbose) for p in products)
if haskey(download_info, platform_key())
    url, tarball_hash = download_info[platform_key()]
    if unsatisfied || !isinstalled(url, tarball_hash; prefix=prefix)
        # Download and install binaries
        install(url, tarball_hash; prefix=prefix, force=true, verbose=verbose)
    end
elseif unsatisfied
    # If we don't have a BinaryProvider-compatible .tar.gz to download, complain.
    # Alternatively, you could attempt to install from a separate provider,
    # build from source or something even more ambitious here.
    error("Your platform $(triplet(platform_key())) is not supported by this package!")
end

# Write out a deps.jl file that will contain mappings for our products
write_deps_file(joinpath(@__DIR__, "deps.jl"), products)

ext: jl

max_stars_repo_name: Gnimuc/BKmaxflow

max_issues_repo_name: Gnimuc/BKmaxflow

max_forks_repo_name: Gnimuc/BKmaxflow


text:
#=

Chase Coleman (ccoleman@stern.nyu.edu)
Spencer Lyon (spencer.lyon@stern.nyu.edu)

References
----------

Judd, Kenneth L, Lilia Maliar, Serguei Maliar, and Rafael Valero. 2013.
    "Smolyak Method for Solving Dynamic Economic Models: Lagrange
    Interpolation, Anisotropic Grid and Adaptive Domain".

Krueger, Dirk, and Felix Kubler. 2004. "Computing Equilibrium in OLG
    Models with Stochastic Production." Journal of Economic Dynamics and
    Control 28 (7) (April): 1411-1436.

=#

# union type to handle isotropic/anisotropic mu parameters
typealias IntSorV Union{Int, Vector{Int}}

num_grid_points(d::Int, mu::Int) =
    mu == 1 ? 2d - 1 :
    mu == 2 ? Int(1 + 4d + 4d*(d-1)/2 ):
    mu == 3 ? Int(1 + 8d + 12d*(d-1)/2 + 8d*(d-1)*(d-2)/6) :
    error("We only know the number of grid points for mu ∈ [1, 2, 3]")

m_i(i::Int) = i < 0 ?  error("DomainError: i must be positive") :
              i == 0 ? 0 :
              i == 1 ? 1 : 2^(i - 1) + 1

function cheby2n{T<:Number}(x::Array{T}, n::Int, kind::Int=1)
    out = Array(T, size(x)..., n+1)
    cheby2n!(out, x, n, kind)
end

@noinline function cheby2n!{T<:Number,N}(out::Array{T}, x::Array{T,N}, n::Int,
                                         kind::Int)
    if size(out) != tuple(size(x)..., n+1)
        error("out must have dimensions $(tuple(size(x)..., n+1))")
    end

    R = CartesianRange(size(x))
    # fill first element with ones
    @inbounds @simd for I in R
        out[I, 1] = one(T)
        out[I, 2] = kind*x[I]
    end

    @inbounds for i in 3:n+1
        @simd for I in R
            out[I, i] = 2x[I] * out[I, i - 1] - out[I, i - 2]
        end
    end
    out
end

function s_n(n::Int)
    n < 1 && error("DomainError: n must be positive")

    if n == 1
        return [0.0]
    end

    m = m_i(n)
    j = 1:m
    pts = cos(pi * (j - 1.0) / (m - 1.0))
    @inbounds @simd for i in eachindex(pts)
        pts[i] = abs(pts[i]) < 1e-12 ? 0.0 : pts[i]
    end

    pts
end

function a_n(n::Int)
    n < 1 && error("DomainError: n must be positive")

    if n == 1
        return [0.0]
    elseif n == 2
        return [-1.0, 1.0]
    end

    sn = s_n(n)

    sn[2:2:end]
end

function a_chain(n::Int)
    sn = s_n(n)
    a = Dict{Int,Vector{Float64}}()
    sizehint!(a, n)

    # These are constant and don't follow the pattern.
    a[1] = [0.0]
    a[2] = [-1.0, 1.0]

    for i=n:-1:3
        # push!(a, sn[[2:2:size(sn, 1)]])
        a[i] = sn[2:2:end]
        sn = sn[1:2:end]
    end

    a
end

function phi_chain(n::Int)
    max_ind = m_i(n)
    phi = Dict{Int, UnitRange{Int64}}()
    phi[1] = 1:1
    phi[2] = 2:3
    low_ind = 4  # current lower index

    for i = 3:n
        high_ind = m_i(i)
        phi[i] = low_ind:high_ind
        low_ind = high_ind + 1
    end

    phi
end

## ---------------------- ##
#- Construction Utilities -#
## ---------------------- ##

# Isotropic inds
function smol_inds(d::Int, mu::Int)

    p_vals = 1:(mu+1)

    # PERF: size_hint here if it is slow
    poss_inds = Vector{Int}[]

    for el in with_replacement_combinations(p_vals, d)
        if d < sum(el) <= d + mu
            push!(poss_inds, el)
        end
    end

    # PERF: size_hint here if it is slow
    true_inds = Vector{Int}[ones(Int64, d)]  # we will always have (1, 1, ...,  1)
    for val in poss_inds
        for el in @task pmute(val)
            push!(true_inds, el)
        end
    end

    return true_inds
end

# Ansotropic inds
function smol_inds(d::Int, mu::Vector{Int})
    # Compute indices needed for anisotropic smolyak grid given number of
    # dimensions d and a vector of mu parameters mu

    length(mu) != d &&  error("ValueError: mu must have d elements.")

    mu_max = maximum(mu)
    mup1 = mu + 1

    p_vals = 1:(mu_max+1)

    poss_inds = Vector{Int}[]

    for el in with_replacement_combinations(p_vals, d)
        if d < sum(el) <= d + mu_max
            push!(poss_inds, el)
        end
    end

    true_inds = Vector{Int}[ones(Int64, d)]
    for val in poss_inds
        for el in @task pmute(val)
            if all(el .<= mup1)
                push!(true_inds, el)
            end
        end
    end

    return true_inds
end

function poly_inds(d::Int, mu::IntSorV, inds::Vector{Vector{Int}}=smol_inds(d, mu))
    phi_n = phi_chain(maximum(mu) + 1)
    vcat([cartprod([phi_n[i] for i in el]) for el in inds]...)::Matrix{Int}
end

# Build grid
function build_grid(d::Int, mu::IntSorV, inds::Vector{Vector{Int}}=smol_inds(d, mu))
    An = a_chain(maximum(mu) + 1)  # use maximum in case mu is Vector
    vcat([cartprod([An[i] for i in el]) for el in inds]...)::Matrix{Float64}
end

# Build B-matrix
function build_B!{T}(out::AbstractMatrix{T}, d::Int, mu::IntSorV,
                     pts::Matrix{Float64}, b_inds::Matrix{Int64})
    # check dimensions
    npolys = size(b_inds, 1)
    npts = size(pts, 1)
    size(out) == (npts, npolys) || error("Out should be size $((npts, npolys))")

    # fill out with ones so tensor product below works
    fill!(out, one(T))

    # compute all the chebyshev polynomials we'll need
    Ts = cheby2n(pts, m_i(maximum(mu) + 1))

    @inbounds for ind in 1:npolys, k in 1:d
        b = b_inds[ind, k]
        for i in 1:npts
            out[i, ind] *= Ts[i, k, b]
        end
    end

    return out
end

build_B(d::Int, mu::IntSorV, pts::Matrix{Float64}, b_inds::Matrix{Int64}) =
    build_B!(Array(Float64, size(pts, 1), size(b_inds, 1)), d, mu, pts, b_inds)

function dom2cube!(out::AbstractMatrix, pts::AbstractMatrix,
                   lb::AbstractVector, ub::AbstractVector)
    d = length(lb)
    n = size(pts, 1)

    size(out) == (n, d) || error("out should be $((n, d))")

    @inbounds for i_d in 1:d
        center = lb[i_d] + (ub[i_d] - lb[i_d])/2
        radius = (ub[i_d] - lb[i_d])/2
        @simd for i_n in 1:n
            out[i_n, i_d] = (pts[i_n, i_d] - center)/radius
        end
    end

    out
end

function cube2dom!(out::AbstractMatrix, pts::AbstractMatrix,
                   lb::AbstractVector, ub::AbstractVector)
    d = length(lb)
    n = size(pts, 1)

    size(out) == (n, d) || error("out should be $((n, d))")

    @inbounds for i_d in 1:d
        center = lb[i_d] + (ub[i_d] - lb[i_d])/2
        radius = (ub[i_d] - lb[i_d])/2
        @simd for i_n in 1:n
            out[i_n, i_d] = center + pts[i_n, i_d]*radius
        end
    end

    out
end

for f in [:dom2cube!, :cube2dom!]
    no_bang = symbol(string(f)[1:end-1])
    @eval $(no_bang){T}(pts::AbstractMatrix{T}, lb::AbstractVector, ub::AbstractVector) =
        $(f)(Array(T, size(pts, 1), length(lb)), pts, lb, ub)
end

## ----------------- ##
#- Type: SmolyakGrid -#
## ----------------- ##

### type SmolyakGrid
# TODO: decide if we should keep B or sparse(B)
type SmolyakGrid{Tmu<:IntSorV,Tlb<:Number,Tub<:Number}
    d::Int  # number of dimensions
    mu::Tmu  # density. Int or d element vector of Int
    lb::Vector{Tlb} # This is the lower bound for the grid (the domain)
    ub::Vector{Tub} # This is the upper bound for the grid (the domain)
    cube_grid::Matrix{Float64}  # Smolyak grid on [-1, 1]^d
    grid::Matrix{Float64} # Smolyak grid on original domain
    inds::Vector{Vector{Int}}  # Smolyak indices
    pinds::Matrix{Int64}  # Polynomial indices
    B::Matrix{Float64}  # matrix representing interpoland
    B_fact::Base.LinAlg.LU{Float64,Matrix{Float64}}  # LU factorization of B

    function SmolyakGrid(d::Int, mu::IntSorV, lb::Vector=-ones(d),
                         ub::Vector=ones(d))
        d < 2 && error("You passed d = $d. d must be greater than 1")
        mu < 1 && error("You passed mu = $mu. mu must be greater than 1")
        if length(mu) > 1
            # working on building an anisotropic grid
            if length(mu) != d
                error("mu must have d elements. It has $(length(mu))")
            end
        end

        inds = smol_inds(d, mu)
        pinds = poly_inds(d, mu, inds)
        cube_grid = build_grid(d, mu, inds)
        grid = cube2dom(cube_grid, lb, ub)
        B = build_B(d, mu, cube_grid, pinds)
        B_fact = lufact(B)

        new(d, mu, lb, ub, cube_grid, grid, inds, pinds, B, B_fact)
    end

    # Default constructor, just in case someone gets brave
    function SmolyakGrid(d, mu, lb, ub, cube_grid, grid, inds, pinds, B, B_fact)
        new(d, mu, lb, ub, cube_grid, grid, inds, pinds, B, B_fact)
    end
end

# outer constructor to infer types
function SmolyakGrid{Tmu<:IntSorV,Tlb<:Number,Tub<:Number}(d::Int,
                                                           mu::Tmu,
                                                           lb::Vector{Tlb}=-ones(d),
                                                           ub::Vector{Tub}=ones(d))
    SmolyakGrid{Tmu,Tlb,Tub}(d, mu, lb, ub)
end

function SmolyakGrid(d::Int, mu::IntSorV, lb::Number, ub::Number)
    lb = fill(lb, d)
    ub = fill(ub, d)
    SmolyakGrid(d, mu, lb, ub)
end

function SmolyakGrid(d::Int, mu::IntSorV, lub::Number)
    lub = abs(lub)
    lb = fill(-lub, d)
    ub = fill(lub, d)
    SmolyakGrid(d, mu, lb, ub)
end

function Base.show(io::IO, sg::SmolyakGrid)
    npoints = size(sg.cube_grid, 1)
    non_zero_pts = countnz(sg.B)
    pct_non_zero = (non_zero_pts / (npoints ^ 2)) * 100
    if isa(sg.mu, Array)
        mu_str = replace(strip(string(sg.mu)), '\n', " x ")
        msg = "Anisotropic Smolyak Grid:\n"
        msg *= "\td: $(sg.d)\n\tmu: $(mu_str)\n\tnpoints: $npoints"
        msg *= @sprintf "\n\tB: %.2f%% non-zero" pct_non_zero
    else
        msg = "Smolyak Grid:\n\td: $(sg.d)\n\tmu: $(sg.mu)\n\tnpoints: $npoints"
        msg *= @sprintf "\n\tB: %.2f%% non-zero" pct_non_zero
    end
    print(io, msg)
end

max_issues_repo_name: EconForge/Smolyak

max_forks_repo_name: EconForge/Smolyak


text:
"""
The `Traitor` package aims to bring a prototype of traits-based dispatch to
Julia. Users will primarily interact with `Traitor` via the `@traitor` macro.

Traits in `Traitor` follow a simple pattern. Firstly, all traits and groupings
of traits are types, and the traits corresponding to a Julia object can be
determined from its type alone. A "class" of traits is defined by an abstract
type, for example `abstract Size`. Examples of each trait are described by its
subtypes, for example `immutable Big <: Size; end` and
`immutable Small <: Size; end`.

Traits can be grouped in a couple of ways. Traits from the *same* class can be
combined in a `Union`, for example `Union{Big, Small}`. This union represents a
*larger* set of Julia objects than either `Big` or `Small` - they represent all
objects possessing either trait.

Conversely, traits from *different* classes can be combined in a `Tuple`, for
example `Tuple{Big, Smelly}` (where `Smelly <: Odor` is another trait from
the `Odor` class). These combinations represent a *smaller* set of Julia objects
than either `Big` or `Smelly` - they represent only objects that possess both
traits simultaneously.

To tell Julia that a type has a certain trait, you overload the constructor of
the trait class (the abstract supertype, e.g. `Size`). For example, we might
have `Size(::Type{Int}) = Small` and `Size{::Type{BigInt}} == Big`. Note that
types are used as both inputs and outputs. A default trait for all objects
can be given, e.g. `Size(::Any) = UnknownSize` (where `UnknownSize <: Size`).

Trait-based methods are created with the `@traitor` macro. Traits are defined
with an extra (optional) set of `::` in the function signature. For example:

    @traitor howbig(::Integer::Big) = "Huge!"
    @traitor howbig(::Integer::Union{Medium, Small}) = "Smallish"
    @traitor howbig(::Integer::Tuple{Big,Smelly}) = "That integer is really big and quite smelly"
    @traitor howbig(::Integer) = "Dunno"

The function is then called in the standard way, and `Traitor` will automatically
manage dispatch. If your code is statically typed, then this will be done at
compile time with no run-time overhead. For example:

    julia> howbig(1)
    "Smallish"

    julia> howbig(BigInt(1))
    "Huge!"

Finally, it is important to note that trait-dispatch occurs after standard
Julia multiple-dispatch. A function may contain many methods, of which
none, some or all will be split into trait-based submethods. The `@traitor`
macro will automatically take ownership of a single method, overwriting any
previous definition.

(see also @traitor, @betray and supertrait)
"""
module Traitor

import Base: @pure
export @traitor, @betray, supertrait

"""
    @betray f

The goal of the `@betray` macro is to take over the methods of `f` and ready
them to be compatible with traits-based dispatch, so that new submethods
defined by `@traitor` do not overwrite existing default methods (and used where
no more specific trait match is found).

WARNING: Please note this macro has NOT been implemented yet! (If anyone knows
how to turn a `Method` object into a new method definition, please let us know).

(See also `Traitor` and `@traitor`.)
"""
macro betray(ex)
    if !isa(ex, Symbol)
        error("Use @betray functionname")
    end

    return quote
        if !isdefined(ex)
            eval(:(function $ex end))
        else
            # TODO
            warn("""The @betray macro isn't currently implemented. The idea is to
                    take the methods of a generic function and put them into a
                    trait-dispatch table, so the user may then specialize these
                    methods further with traits while not losing all the
                    pre-existing definitions.""")
        end
    end
end

"""
    supertrait(trait)

For a simple trait type, returns `supertype(t)`, while for a `Union` of traits
is returns their common supertype (or else throws an error).
"""
supertrait(t::DataType) = supertype(t)
supertrait(t::Union) = supertrait_union(t)
Base.@pure function supertrait_union(t)
    traitclass = supertype(t.types[1])
    for j = 2:length(t.types)
        if traitclass != supertype(t.types[j])
            error("Unions of traits must be in the same class. Got $t")
        end
    end
    return traitclass
end
supertrait(t) = error("Unknown trait $t")

"""
    extract_arg_trait(ex) -> (normal_argument, trait)

Parse an argument signature to extract both the "normal"
typed-variable declaration and the trait expression.

If no trait is present in the expression, return `Tuple{}` for
`trait`.

# Examples

```julia
extract_arg_trait(:(x::Int::Big)) -> (:(x::Int), :(Tuple{Big}))
```
"""
function extract_arg_trait(ex)
    if isa(ex, Symbol)
        return Expr(:(::), ex, :Any), Tuple{} # argument with no type annotation
    end
    if ex.head == :parameters
        # TODO
        error("Keyword arguments are not supported by @traitor")
    end
    if ex.head == :kw # default values, eg f(x::T::A=1)
        # TODO currently these are ignored!!
        warn("Default values are currently ignored by @traitor...")
        arg, traittype = extract_arg_trait(ex.args[1])
        return Expr(:kw, arg, ex.args[2]), traittype
    end
    @assert ex.head == :(::)
    if length(ex.args) == 1 # of form ::Type, etc. We need an actual symbol name as we refer to it in the generated function
        ex.args = Any[gensym(), ex.args[1]]
    end
    isfourcolon = isa(ex.args[2],Expr) && ex.args[2].head == :(::) # x::::A syntax
    if isa(ex.args[1], Symbol) && !isfourcolon
        return ex, Tuple{} # argument with normal type annotation
    end
    if isfourcolon
        traittype = ex.args[2].args[1]
        if !isa(traittype, Expr) || traittype.head != :curly || traittype.args[1] != :Tuple
            traittype = Expr(:curly, :Tuple, traittype)
        end
        arg = (isa(ex.args[1], Symbol) ? Expr(:(::), ex.args[1], :Any) : ex.args[1])
        return arg, traittype
    end
    normalarg = ex.args[1]
    @assert normalarg.head == :(::)
    traittype = ex.args[2]
    if !isa(traittype, Expr) || traittype.head != :curly || traittype.args[1] != Tuple
        traittype = Expr(:curly, :Tuple, traittype)
    end
    return normalarg, traittype
end

"""
Each standard Julia method of a betrayed function may contain an arbitrary
number of trait-based submethods. This generated function returns (and creates,
if necessary, the first time it is called) the table joining the traits to
the submethods
"""
@generated function get_trait_table{Signature<:Tuple}(f::Function, ::Type{Signature})
    d = Dict{Any, Function}()
    return quote
        #$(Expr(:meta,:inline))
        ($d)
    end
end

"""
    @traitor f(x1::T1::Trait1, x2::T2::Trait, ...) =
    @traitor function f(x1::T1::Trait1, x2::T2::Trait, ...); ...; end

The @traitor takes over a method of function `f` to use trait-based dispatch.
If a prior non-@traitor method is defined, then that method will be deleted.
However, methods of `f` with signatures other than `f(::T1, ::T2, ...)` remain
unaffected.

Dispatch on traits occurs after normal Julian multiple-dispatch for further
specialization to sub-methods. The most specific matching trait-signature is
chosen, or an error is thrown in the case of ambiguities.

NOTE: Currently, default values and keyword arguments are not supported.

(See help for `Traitor` to learn how traits are to be used and defined.)
"""
macro traitor(ex)
    isa(ex, Expr) || error("trait expression must be a function method definition")

    # If it's a compact function definition `f(x) = y`, we convert it over to long form function `f(x); y; end`
    if ex.head == :(=) && isa(ex.args[1], Expr) && ex.args[1].head == :call
        body = ex.args[2]
        if !isa(body, Expr) || body.head != :block # Many one-liners don't have a block...
            body = Expr(:block, body)
        end
        ex = Expr(:function, ex.args[1], ex.args[2])
    end

    ex.head == :function || error("trait expression must be a function method definition")

    # Get the function name, signature and traits
    def = ex.args[1]
    body = ex.args[2]

    (isa(def, Expr) && def.head == :call) || error("trait expression must be a function method definition")
    funcname = def.args[1]

    # Go through all the arguments
    n_args = length(def.args) - 1
    args = Vector{Any}()
    argnames = Vector{Any}()
    quotednames = Vector{Any}()
    argtypes = Vector{Any}()
    traits = Vector{Any}()
    for i = 1:n_args
        arg, trait = extract_arg_trait(def.args[1+i])
        push!(args, arg)
        push!(quotednames, Expr(:quote, arg.args[1]))
        push!(argnames, arg.args[1])
        push!(argtypes, arg.args[2])
        push!(traits, trait)
    end

    argnames = (argnames...) # I'm conused why $((argnames...)) doesn't work in the quote below
    quotednames = (quotednames...)

    # Make a new name for this specialized function
    internalname = "_"*string(funcname)*"{"
    for i = 1:length(argtypes)
        internalname = internalname * "::" * string(argtypes[i]) * "::" * string(traits[i])
        if i < length(argtypes)
            internalname = internalname * ","
        end
    end
    internalname = internalname * "}"
    internalname = Symbol(internalname)

    # It's hard to get all of this right with nest quote blocks, AND it's hard
    # to get this right with Expr() objects... grrr...
    esc(Expr(:block,
        Expr(:stagedfunction, Expr(:call, funcname, args...), Expr(:block,
            :( dict = Traitor.get_trait_table($funcname, $(Expr(:curly, :Tuple, argtypes...))) ),
            :( f = Traitor.trait_dispatch(dict, $(Expr(:curly, :Tuple, argnames...))) ),
            Expr(:quote, Expr(:block,
                Expr(:meta, :inline),
                Expr(:call, Expr(:$, :f), argnames...)
            ))
        )),
        Expr(:function, Expr(:call, internalname, args...), body),
        :( d = Traitor.get_trait_table($funcname, $(Expr(:curly, :Tuple, argtypes...))) ),
        :( d[$(Expr(:curly, :Tuple, traits...))] = $internalname ),
    ))
end

"""
This function deterimines the dispatch on the traits. Each instance of
this generated function is specialized on the (standard) signature of the
inputs, so the only task remaining is to find the most specific matching
"trait signature".
"""
function trait_dispatch{Sig <: Tuple}(trait_dictionary::Dict{Any, Function}, ::Type{Sig})
    n_args = length(Sig.parameters)

    # First check which (if any) of our trait conditions is satisfied by Sig
    matching_traits = Vector{Any}()
    for traits ∈ keys(trait_dictionary)
        match = true
        # Iterate and search for lack of matches
        for n in 1:n_args
            traits_n = traits.parameters[n]

            # The traits for argument i should be a Tuple type
            for i in 1:length(traits_n.parameters)
                traitclass = supertrait(traits_n.parameters[i])

                if !(traitclass(Sig.parameters[n]) <: traits_n.parameters[i])
                    match = false
                    break
                end
            end
        end
        if match
            push!(matching_traits, traits)
        end
    end

    # Now see if there is exactly one most specific matching trait signature
    if length(matching_traits) == 0
        error("Can't find any matching trait method. Got signature $Sig and traits $(collect(keys(trait_dictionary)))")
    elseif length(matching_traits) == 1
        return trait_dictionary[matching_traits[1]]
    else # length(matching_traits) > 1
        # See if there is any matching trait that is more specific than *all*
        # other matching traits, otherwise throw an ambiguity error
        for i = 1:length(matching_traits)
            most_specific = true
            for j = 1:length(matching_traits)
                if i == j
                    continue
                end

                if is_more_specific_traitsig(matching_traits[j], matching_traits[i])
                    most_specific = false
                    break
                end
            end
            if most_specific
                return trait_dictionary[matching_traits[i]]
            end
        end
        error("Ambiguity. Could not determine most specific matching trait. Options are: $matching_traits")
    end
end

function is_more_specific_traitsig(traits, traits2)
    n = length(traits.parameters)

    more_specific = true
    for i = 1:n
        more_specific = more_specific & is_more_specific_traitarg(traits.parameters[i], traits2.parameters[i])
    end

    return more_specific
end

function is_more_specific_traitarg(traits, traits2)
    # Now arrange these into trait classes
    d = Dict{Any, Any}()
    for trait in traits.parameters
        traitclass = supertrait(trait)
        if haskey(d, traitclass)
            error("Encountered trait class $d more than once in $traits")
        else
            d[traitclass] = trait
        end
    end

    d2 = Dict{Any, Any}()
    for trait in traits2.parameters
        traitclass = supertrait(trait)
        if haskey(d2, traitclass)
            error("Encountered trait class $d more than once in $traits2")
        else
            d2[traitclass] = trait
        end
    end

    # For each trait, calculate what is going on
    more_specific = true

    for traitclass in union(keys(d), keys(d2))
        if haskey(d, traitclass)
            if haskey(d2, traitclass)
                if !(d[traitclass] <: d2[traitclass])
                    more_specific = false
                    break
                end
            end
        elseif haskey(d2, traitclass)
            # only in old trait
            more_specific = false
            break
        end
    end

    return more_specific
end

end # module


text:
using ZXCalculus
using YaoPlots
using Plots

function generate_example()
    zxd = ZXDiagram(4)
    push_gate!(zxd, Val{:Z}(), 1, 3//2)
    push_gate!(zxd, Val{:H}(), 1)
    push_gate!(zxd, Val{:Z}(), 1, 1//2)
    push_gate!(zxd, Val{:H}(), 4)
    push_gate!(zxd, Val{:CZ}(), 4, 1)
    push_gate!(zxd, Val{:CNOT}(), 1, 4)
    push_gate!(zxd, Val{:H}(), 1)
    push_gate!(zxd, Val{:H}(), 4)
    push_gate!(zxd, Val{:Z}(), 1, 1//4)
    push_gate!(zxd, Val{:Z}(), 4, 3//2)
    push_gate!(zxd, Val{:X}(), 4, 1//1)
    push_gate!(zxd, Val{:H}(), 1)
    push_gate!(zxd, Val{:Z}(), 4, 1//2)
    push_gate!(zxd, Val{:X}(), 4, 1//1)

    push_gate!(zxd, Val{:Z}(), 2, 1//2)
    push_gate!(zxd, Val{:CNOT}(), 3, 2)
    push_gate!(zxd, Val{:H}(), 2)
    push_gate!(zxd, Val{:CNOT}(), 3, 2)
    push_gate!(zxd, Val{:Z}(), 2, 1//4)
    push_gate!(zxd, Val{:Z}(), 3, 1//2)
    push_gate!(zxd, Val{:H}(), 2)
    push_gate!(zxd, Val{:H}(), 3)
    push_gate!(zxd, Val{:Z}(), 3, 1//2)
    push_gate!(zxd, Val{:CNOT}(), 3, 2)

    return zxd
end

zxd = generate_example()
ex_zxd = clifford_simplification(zxd)

plot(zxd)
plot(ex_zxd)

macro Name(arg)
   string(arg)
end

name = @Name ex1

#function save_plot(name, plot, label)
    #cd(dirname(@__FILE__));
    #dir = pwd();
    #println("Saving plots to the directory in $dir","\\Graphs")
    #strlabel = string(label);
    #mkpath(string(dir,"\\Graphs\\","\\",name,))
#    pngfilename = string(dir,"\\Graphs\\","\\",name,"\\","\\",name,strlabel,".png")
    #pdffilename = string(dir,"\\Graphs\\","\\",name,"_",q, " qubits\\","\\",step," steps\\",strlabel,".pdf")
    #epsfilename = string(dir,"\\Graphs\\","\\",name,"_",q, " qubits\\","\\",step," steps\\",strlabel,".eps")
    #savefig(plot, pngfilename)
    #savefig(plot, pdffilename)
    #savefig(plot, epsfilename)
#end

#title = "" *string(plot_name)* ""
#Plots.title!(title)
#lables=title

#save_plot(plot_name, zxd, label) # Saves the plots to github

#dir = pwd();
#println("Saving plots to the directory in $dir","\\Graphs")
#mkpath(string(dir,"\\Graphs\\","\\",name,))
#pngfilename = string(dir,"\\Graphs\\","\\",name,"\\","\\",name,".png")
#pdffilename = string(dir,"\\Graphs\\","\\",name,"_",q, " qubits\\","\\",step," steps\\",strlabel,".pdf")
#epsfilename = string(dir,"\\Graphs\\","\\",name,"_",q, " qubits\\","\\",step," steps\\",strlabel,".eps")
#savefig(zxd, pngfilename)


text:
"""
    immutable Selector
      dict::Dict{AbstractString, Any}
    end

Immutable represention of a Mango Selector used to query a Mango index.

Usually created using the custom string literal `q"..."` (see the `@q_str` macro),
but can be created directly from either the raw json string containing a Selector
expression or a Julia `Dict(...)` representing the same.

### Examples

    sel = q"name = bob"
    sel = Selector("{\"name\":{\"\$eq\":\"bob\"}}")
    sel = Selector(Dict("name" => Dict("\$eq" => "bob")))
    sel = and([q"name = bob", q"age > 18"])

[API reference](https://docs.cloudant.com/cloudant_query.html#selector-syntax)
"""
immutable Selector
  dict::Dict{AbstractString, Any}
end

"""
    Selector()

The empty Selector.
"""
function Selector()
  Selector(Dict{AbstractString, Any}())
end

"""
    Selector(raw_json::AbstractString)

Create a Selector from the raw json.

[API reference](https://docs.cloudant.com/cloudant_query.html#selector-syntax)
"""
function Selector(raw_json::AbstractString)
  Selector(JSON.parse(raw_json))
end

"""
    isempty(sel::Selector)

True if sel is the empty Selector.
"""
function isempty(sel::Selector)
  length(sel.dict) == 0
end

"""
    q"....."

Custom string literal for a limited Selector definition DSL.

It takes the form:

    field op data

where `field` is a field name, op is one of

    =, !=, <, <=, >, >=, in, !in, all

This allows you to write things like:

    q"name = bob"
    q"value < 5"
    q"occupation in [fishmonger, pilot, welder]"

Note that the Selector DSL only covers a fraction of the full Selector
syntax. It can be used with the boolean functions `and()`, `or()` etc
to build up more complex Selectors, e.g.

  sel = and([q"name = bob", q"age > 18"])

For more information on the actual Selector syntax, see link below.

[API reference](https://docs.cloudant.com/cloudant_query.html#selector-syntax)
"""
macro q_str(data)
  quote
    operators = Dict{String, String}(
      "=="  => "\$eq",
      "="   => "\$eq",
      "!="  => "\$ne",
      "<"   => "\$lt",
      "<="  => "\$lte",
      ">"   => "\$gt",
      ">="  => "\$gte",
      "in"  => "\$in",
      "!in" => "\$nin",
      "all" => "\$all"
    )

    m = match(r"^(.+?)\s*(==|=|!=|<|<=|>|>=|in|!in|all)\s*\[?(.+)\]?$", strip($(data)))
    if m == nothing
      error("Badly formatted selector string")
    end

    field    = m.captures[1]
    operator = m.captures[2]
    value    = m.captures[3]

    if !haskey(operators, operator)
      error("Unknown operator '$operator'")
    end

    if operator in ["in", "!in", "all"]
      value = map(strip, split(value, ","))
    end

    mydict = Dict(field => Dict(operators[operator] => value))
    Selector(mydict)
  end
end

# Boolean logic composition of Selectors:
#
# sel = and([q"name = bob", q"age > 18"])
# sel = or( [q"name = bob", q"age > 18"])
# sel = nor([q"name = bob", q"age > 18"])
#
for boolop in [:and, :or, :nor]
  boolop_str = "\$"*string(boolop)
  @eval begin
    function ($boolop)(sel::Vector{Selector})
      Selector(Dict($boolop_str => map(sel) do s
        s.dict
      end))
    end
  end
end

# Selector modifiers:
#
# sel = not(and([q"name = bob", q"age > 18"]))
#
for op in [:not]
  op_str = "\$"*string(op)
  @eval begin
    function ($op)(sel::Selector)
      Selector(Dict($op_str => sel.dict))
    end
  end
end


text:
import ArgParse, MPI, HDF5, RIOPA
import Test: @testset, @test, @test_throws

@testset "hello" begin
    RIOPA.hello(nothing)
    if HDF5.has_parallel()
        @test_exists_and_rm "hello.h5"
    else
        @test_exists_and_rm "hello_0.h5"
    end
end


text:
#=
type RobustNoCrashMDP <: RobustMDP{MLPhysicalState, MLAction}
    base::MLMDP
end
RobustMCTS.representative_mdp(rmdp::RobustNoCrashMDP) = StochasticBehaviorNoCrashMDP(rmdp.base)

function RobustMCTS.next_model(gen::RandomModelGenerator, rmdp::RobustNoCrashMDP, s::MLPhysicalState, a::MLAction)
    behaviors = Dict{Int, Union{Nothing,BehaviorModel}}(1=>nothing)
    for c in s.cars
        # TODO instead of just using sample, try to pick models that will be bad
        behaviors[c.id] = rand(gen.rng, rmdp.base.dmodel.behaviors)
    end
    return FixedBehaviorNoCrashMDP(behaviors, rmdp.base)
end
=#

abstract type EmbeddedBehaviorMDP <: MDP{MLPhysicalState, MLAction} end

POMDPs.actions(p::EmbeddedBehaviorMDP) = actions(p.base)
POMDPs.actions(p::EmbeddedBehaviorMDP, s::MLPhysicalState, as::NoCrashActionSpace) = actions(p.base, s, as)
# create_action(p::EmbeddedBehaviorMDP) = create_action(p.base)
create_state(p::EmbeddedBehaviorMDP) = MLPhysicalState(false, 0.0, 0.0, [])
POMDPs.discount(p::EmbeddedBehaviorMDP) = discount(p.base)
# reward(p::EmbeddedBehaviorMDP, s::MLPhysicalState, a::MLAction, sp::MLPhysicalState)

mutable struct FixedBehaviorNoCrashMDP <: EmbeddedBehaviorMDP
    behaviors::Dict{Int,Union{Nothing,BehaviorModel}}
    base::MLMDP
end

function POMDPs.gen(mdp::FixedBehaviorNoCrashMDP, s::MLPhysicalState, a::MLAction, rng::AbstractRNG)
    full_s = MLState(s, Vector{CarState}(undef, length(s.cars)))
    for (i,c) in enumerate(s.cars)
        full_s.cars[i] = CarState(c, mdp.behaviors[c.id])
    end
    full_sp = @gen(:sp)(mdp.base, full_s, a, rng)
    sp = create_state(mdp)
    for c in sp.cars
        if !haskey(mdp.behaviors, c.id)
            mdp.behaviors[c.id] = c.behavior
        end
    end
    return (sp=MLPhysicalState(full_sp), r=reward(mdp.base, full_s, a, full_sp), info=nothing)
end

mutable struct StochasticBehaviorNoCrashMDP <: EmbeddedBehaviorMDP
    base::MLMDP
end

function POMDPs.gen(mdp::StochasticBehaviorNoCrashMDP, s::MLPhysicalState, a::MLAction, rng::AbstractRNG)
    full_s = MLState(s, Vector{CarState}(undef, length(s.cars)))
    for i in 1:length(s.cars)
        full_s.cars[i] = CarState(s.cars[i], rand(rng, mdp.base.dmodel.behaviors))
    end
    full_sp = @gen(:sp)(mdp.base, full_s, a, rng)
    return (sp=MLPhysicalState(full_sp), r=reward(mdp.base, full_s, a, full_sp), info=nothing)
end

function POMDPs.initialstate(mdp::Union{FixedBehaviorNoCrashMDP,StochasticBehaviorNoCrashMDP})
    ImplicitDistribution() do rng
        full_state = rand(rng, initialstate(mdp.base))
        return MLPhysicalState(full_state)
    end
end

function initial_state(mdp::Union{FixedBehaviorNoCrashMDP,StochasticBehaviorNoCrashMDP}, rng::AbstractRNG)
    full_state = initial_state(mdp.base, rng)
    return MLPhysicalState(full_state)
end

#=
type RobustMLSolver <: Solver
    rsolver
end

type RobustMLPolicy <: Policy{MLState}
    rpolicy
end

function action(p::RobustMLPolicy, s::MLState, a::MLAction=MLAction(0,0))
    mdp = representative_mdp(p.rpolicy.rmdp)
    as = actions(mdp, MLPhysicalState(s), actions(mdp))
    if length(as) == 1
        return collect(as)[1]
    end
    action(p.rpolicy, MLPhysicalState(s))
end

function solve(s::RobustMLSolver, mdp::MLMDP)
    RobustMLPolicy(solve(s.rsolver, RobustNoCrashMDP(mdp::MLMDP)))
end
=#


text:
# Exercise 12-1
function sumall(args...)
    return sum(args)
end

# println(sumall(1,2,3,4,5))

# Exercise 12-2
function histogram(s)
    d = Dict()
    for c in s
        d[c] = get(d, c, 0) + 1
    end
    d
end

function keyvaluepairs(h)
    lst = Array{Tuple, 1}()
    for (k, v) in h
        push!(lst, (k, v))
    end
    return lst
end

function mostfrequent(str)
    lst = keyvaluepairs(histogram(str))
    sort!(lst, by = x -> x[2])
    return lst
end

# println(mostfrequent("brontosaurus"))

# Exercise 12-3
# Hashing by FTOA.

# Exercise 12-4
function readerhist()
    fin = open("words.txt")
    words = Dict{String, Dict}()
    for line in eachline("words.txt")
        words[line] = histogram(line)
    end
    return words
end

function invertdict(d)
    inverse = Dict()
    for key in keys(d)
        val = d[key]
        push!(get!(inverse, val, []), key)
    end
    return inverse
end

function metathesis(words)
    chartowords = invertdict(words)
    for (word, hist) in words
        if !haskey(chartowords, hist)
            continue
        end
        candidate = chartowords[hist]
        if length(candidate) > 1
            println(candidate)
            delete!(chartowords, hist)
        end
    end
end

# metathesis(readerhist())

# Exercise 12-5
# Finding longest path in DAG, toposort then DP.


text:
"""
    gewekediag(x::AbstractVector{<:Real}; first::Real=0.1, last::Real=0.5, kwargs...)

Compute the Geweke diagnostic from the `first` and `last` proportion of samples `x`.
"""
function gewekediag(x::AbstractVector{<:Real}; first::Real=0.1, last::Real=0.5, kwargs...)
    0 < first < 1 || throw(ArgumentError("`first` is not in (0, 1)"))
    0 < last < 1 || throw(ArgumentError("`last` is not in (0, 1)"))
    first + last <= 1 || throw(ArgumentError("`first` and `last` proportions overlap"))

    n = length(x)
    x1 = x[1:round(Int, first * n)]
    x2 = x[round(Int, n - last * n + 1):n]
    z =
        (Statistics.mean(x1) - Statistics.mean(x2)) /
        hypot(mcse(x1; kwargs...), mcse(x2; kwargs...))
    p = SpecialFunctions.erfc(abs(z) / sqrt(2))

    return (zscore=z, pvalue=p)
end


text:
# # Welford's algorithm
#
# * Welford, B. P. (1962). "Note on a method for calculating corrected sums of squares and products". Technometrics. 4 (3): 419–420. doi:10.2307/1266577. JSTOR 1266577
# * [Welford's online algorithm - Algorithms for calculating variance - Wikipedia](https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm)

using Transducers
using Transducers: R_, start, next, complete, inner, xform, wrap, unwrap, wrapping

struct MeanVar <: Transducer
end

function Transducers.start(rf::R_{MeanVar}, result)
    private_state = (0, 0.0, 0.0)
    return wrap(rf, private_state, start(inner(rf), result))
end

function Transducers.next(rf::R_{MeanVar}, result, input)
    wrapping(rf, result) do st, iresult
        (n, μ, M2) = st
        n += 1
        δ = input - μ
        μ += δ/n
        δ2 = input - μ
        M2 += δ*δ2
        iinput = (μ, M2 / (n-1))
        iresult = next(inner(rf), iresult, iinput)
        return (n, μ, M2), iresult
    end
end

function Transducers.complete(rf::R_{MeanVar}, result)
    _private_state, inner_result = unwrap(rf, result)
    return complete(inner(rf), inner_result)
end

collect(MeanVar(),randn(10))

@time foldl(right, MeanVar(), 1:100)


text:
import StateSpaceReconstruction:
    cembed,
    assign_bin_labels,
    groupslices, groupinds

import PerronFrobenius:
        InvariantDistribution,
        get_binvisits,
        estimate_transferoperator_from_binvisits,
        invariantmeasure

import StaticArrays:
    SVector, MVector

"""
	marginal_indices(A)

Returns a column vector with the same number of elements as there are unique
rows in A. The value of the ith element is the row indices of rows in A
matching the ith unique row.
"""
marginal_indices(A) = groupinds(groupslices(A, 1))

"""
    marginal(along_which_axes::Union{Vector{Int}, UnitRange{Int}},
            visited_bin_labels::Array{Int, 2},
            iv::PerronFrobenius.InvariantDistribution)

Compute the marginal distributions along the specified axes, taking
into account the invariant distribution over the bins.
"""
function marginal(along_which_axes::Union{Vector{Int}, UnitRange{Int}},
                visited_bin_labels::Array{Int, 2},
                iv::InvariantDistribution)

    marginal_inds::Vector{Vector{Int}} =
        marginal_indices(visited_bin_labels[:, along_which_axes])

    # Only loop over the nonzero entries of the invariant distribution
    nonzero_elements_of_dist = iv.dist[iv.nonzero_inds]
    marginal = zeros(Float64, size(marginal_inds, 1))
    @inbounds for i = 1:size(marginal_inds, 1)
        marginal[i] = sum(nonzero_elements_of_dist[marginal_inds[i]])
    end
    return marginal
end

"""
    transferentropy_transferoperator_grid(
            bins_visited_by_orbit::Array{Int, 2},
            iv::PerronFrobenius.InvariantDistribution,
            v::TEVars)

Using the invariant probability distribution obtained from the
transfer operator to the visited bins of the partitioned state
space. `bins_visited_by_orbit` are the bin labels assigned to
each element of the partition that gets visited by the orbit.

We calculate transfer entropy from the embedding
`E`, given a discretization scheme controlled by `\epsilon` and
information `v::TEVars` about which columns of the embedding to
consider for each of the marginal distributions. From these
marginal distributions, we calculate marginal entropies and
insert these into the transfer entropy expression.
"""
function transferentropy_transferoperator_grid(
            bins_visited_by_orbit::Array{Int, 2},
            iv::InvariantDistribution,
            v::TEVars; b = 2)

    # Verify that the number of dynamical variables in
    # the embedding agrees with the number of dynamical
    # variables in `v`(which indicates which variables of
    # the embedding should be considered part of the
    # transfer entropy computation). If these do not
    # agree, the marginals are computed wrongly, and
    # the transfer entropy estimate is not correct
    # (it could even be a negative value, which is
    # nonsensical).
    nvars = length(v.target_future) +
            length(v.target_presentpast) +
            length(v.source_presentpast) +
            length(v.conditioned_presentpast)

    dim = size(bins_visited_by_orbit, 1) # one variable per row

    if dim != nvars
        warn(""" The embedding dimension does not agree with the
            provided $v. The transfer entropy estimate is not
            correct.
            """)
    end

    unique_visited_bins = transpose(unique(bins_visited_by_orbit, dims = 2))
    positive_measure_bins = unique_visited_bins[iv.nonzero_inds, :]

    C = v.conditioned_presentpast
    XY = [v.target_future;      v.target_presentpast; C]
    YZ = [v.target_presentpast; v.source_presentpast; C]
    Y =  [v.target_presentpast;                       C]

    p_Y  = marginal(Y, positive_measure_bins, iv)
    p_XY = marginal(XY, positive_measure_bins, iv)
    p_YZ = marginal(YZ, positive_measure_bins, iv)

    te = StatsBase.entropy(p_YZ, b) +
        StatsBase.entropy(p_XY, b) -
        StatsBase.entropy(p_Y, b) -
        StatsBase.entropy(iv.dist[iv.nonzero_inds], b)

    return te
end

function transferentropy_transferoperator_grid(
                    pts::Vector{T},
                    ϵ::Union{Int, Float64, Vector{Float64}, Vector{Int}},
                    v::TEVars;
                    allocate_frac = 1.0, b = 2) where {T <: Union{Vector, SVector, MVector}}
    transferentropy_transferoperator_grid(hcat(pts...,), ϵ, v,
        allocate_frac = allocate_frac, b = b)
end

"""
	transferentropy_transferoperator_grid(
        E::Embeddings.AbstractEmbedding,
        ϵ::Union{Int, Float64, Vector{Float64}, Vector{Int}},
        v::TransferEntropy.TEVars;
        allocate_frac = 1, b = 2) -> Float64

Using the transfer operator to calculate probability
distributions,  calculate transfer entropy from the embedding
`E`, given a discretization scheme controlled by `ϵ` and
information `v::TEVars` about which columns of the embedding to
consider for each of the marginal distributions. From these
marginal distributions, we calculate marginal entropies and
insert these into the transfer entropy expression.
"""
function transferentropy_transferoperator_grid(
                    E::Embeddings.AbstractEmbedding,
                    ϵ::Union{Int, Float64, Vector{Float64}, Vector{Int}},
                    v::TEVars;
                    allocate_frac = 1.0, b = 2)

    # Verify that the number of dynamical variables in
    # the embedding agrees with the number of dynamical
    # variables in `v`(which indicates which variables of
    # the embedding should be considered part of the
    # transfer entropy computation). If these do not
    # agree, the marginals are computed wrongly, and
    # the transfer entropy estimate is not correct
    # (it could even be a negative value, which is
    # nonsensical).
    nvars = length(v.target_future) +
    		length(v.target_presentpast) +
    		length(v.source_presentpast) +
    		length(v.conditioned_presentpast)

    dim = size(E.points, 1) # one variable per row

    if dim != nvars
        warn(""" The embedding dimension $dim does not agree with the
            provided $v. The transfer entropy estimate is not
            correct.
            """)
    end

    # Identify which bins of the partition resulting
    # from using ϵ each point of the embedding visits.
    bins_visited_by_orbit = assign_bin_labels(E, ϵ)

    # Which are the visited bins, which points
    # visits which bin, repetitions, etc...
    binvisits = get_binvisits(bins_visited_by_orbit)

    # Use that information to estimate transfer operator
    TO = estimate_transferoperator_from_binvisits(binvisits,
                        allocate_frac = allocate_frac)

    # Calculate the invariant distribution over the bins.
    invdist = invariantmeasure(TO)

    transferentropy_transferoperator_grid(
        bins_visited_by_orbit, invdist, v, b = b)
end

"""
    transferentropy_transferoperator_grid(E::Embeddings.AbstractEmbedding,
        ϵ::Vector{Union{Int, Float64, Vector{Float64}, Vector{Int}}},
        v::TEVars; allocate_frac = 1.0, b = 2)

Compute transfer entropy over a range of bin sizes.

Using the transfer operator to calculate probability
distributions,  calculate transfer entropy from the embedding
`E`, given a discretization scheme controlled by the `ϵ`s and
information `v::TEVars` about which columns of the embedding to
consider for each of the marginal distributions. From these
marginal distributions, we calculate marginal entropies and
insert these into the transfer entropy expression.
"""
function transferentropy_transferoperator_grid(E::Embeddings.AbstractEmbedding,
        ϵ::Vector{Union{Int, Float64, Vector{Float64}, Vector{Int}}},
        v::TEVars; allocate_frac = 1.0, b = 2)
    map(ϵᵢ -> transferentropy_transferoperator_grid(
        E, ϵᵢ, v;
        allocate_frac = allocate_frac), ϵ)
end

# Shorter alias
tetogrid = transferentropy_transferoperator_grid

#############################################################################
# One should be able to just provide some points that has been pre-embedded too
#############################################################################
"""
    transferentropy_transferoperator_grid(pts::AbstractArray{T, 2},
        ϵ::Union{Int, Float64, Vector{Float64}, Vector{Int}},
        v::TransferEntropy.TEVars)

Using the transfer operator to calculate probability distributions,
calculate transfer entropy from the points `pts`,
given a discretization scheme controlled by `ϵ` and
information `v::TEVars` about which columns of the embedding to
consider for each of the marginal distributions. From these
marginal distributions, we calculate marginal entropies and
insert these into the transfer entropy expression.
The points will be embedded behind the scenes.
"""
function transferentropy_transferoperator_grid(pts::AbstractArray{T, 2},
    ϵ::Union{Int, Float64, Vector{Float64}, Vector{Int}},
    v::TEVars;
    allocate_frac = allocate_frac, b = 2) where T

    tetogrid(cembed(pts), ϵ, v;
        allocate_frac = allocate_frac, b = b)
end

tetogrid = transferentropy_transferoperator_grid

export
transferentropy_transferoperator_grid,
tetogrid


text:
set_global_client()

@testset "AUTH SELECT PING" begin
    @test_throws RedisError auth("")
    @test select(0) == "OK"
    @test ping() == "PONG"
end

@testset "GET SET" begin
    @test set("key", "value") == "OK"
    @test setnx("key", "value") == 0
    @test get("key") == "value"
    @test Jedis.keys("k*") == ["key"]
    @test del("key") == 1
    @test setnx("key", "value") == 1
    @test flushdb() == "OK" && isnothing(get("key"))
end

@testset "EXPIRE" begin
    @test setex("key", 10, "value") == "OK"
    @test ttl("key") > 0
    @test del("key") == 1 && isnothing(get("key"))
    @test set("key", "value") == "OK"
    @test ttl("key") == -1
    @test expire("key", 10) == 1
    @test ttl("key") > 0
    @test flushdb() == "OK" && isnothing(get("key"))
end

@testset "MULTI EXEC" begin
    @test multi() == "OK"
    @test set("key", "value") == "QUEUED"
    @test get("key") == "QUEUED"
    @test get("key") == "QUEUED"
    @test exec() == ["OK", "value", "value"]
    @test ["OK", "value", "value"] == multi_exec() do
        set("key", "value")
        get("key")
        get("key")
    end
    @test flushall() == "OK" && isnothing(get("key"))
end

@testset "HASH" begin
    @test hset("key", "f1", 1, "f2", 2, "f3", 3) == 3
    @test hget("key", "f2") == "2"
    @test hgetall("key") == ["f1", "1", "f2", "2", "f3", "3"]
    @test hmget("key", "f1", "f2", "doesnotexist") == ["1", "2", nothing]
    @test hdel("key", "f1") == 1
    @test isnothing(hget("key", "f1"))
    @test flushall() == "OK" && isnothing(get("key"))
end

@testset "LIST" begin
    @test lpush("mylist", 3, 2, 1) == 3
    @test rpush("mylist", 4, 5, 6) == 6
    @test llen("mylist") == 6
    @test lrange("mylist", 0, -1) == ["1", "2", "3", "4", "5", "6"]
    @test lpop("mylist") == "1"
    @test rpop("mylist") == "6"
    @test lrange("mylist", 0, -1) == ["2", "3", "4", "5"]
    flushall()
    @test begin
        task = @async blpop("mylist", "otherlist") == ["otherlist", "1"]
        other = Client()
        lpush("otherlist", 1; client=other)
        disconnect!(other)
        fetch(task)
    end
    @test begin
        task = @async brpop("mylist", "otherlist") == ["mylist", "6"]
        other = Client()
        lpush("mylist", 6; client=other)
        disconnect!(other)
        fetch(task)
    end
    flushall()
end

# @testset "QUIT" begin
#     @test quit() == "OK"
#     @test_throws Base.IOError ping()
# end


text:
export Curve

struct IntervalCurve{S<:Space,T,VT} <: SegmentDomain{T}
    curve::Fun{S,T,VT}
end

struct PeriodicCurve{S<:Space,T,VT} <: PeriodicDomain{T}
    curve::Fun{S,T,VT}
end

"""
`Curve` Represents a domain defined by the image of a Fun.  Example
usage would be

```julia
x=Fun(1..2)
Curve(exp(im*x))  # represents an arc
```
"""
const Curve{S,T} = Union{IntervalCurve{S,T},PeriodicCurve{S,T}}

==(a::Curve, b::Curve) = a.curve == b.curve
isempty(::Curve) = false

for TYP in (:IntervalCurve,:PeriodicCurve)
    @eval points(c::$TYP,n::Integer) = c.curve.(points(domain(c.curve),n))
end

checkpoints(d::Curve) = fromcanonical.(Ref(d),checkpoints(domain(d.curve)))

for op in (:(leftendpoint),:(rightendpoint),:(rand))
    @eval $op(c::Curve)=c.curve($op(domain(c.curve)))
end

canonicaldomain(c::Curve) = domain(c.curve)

fromcanonical(c::Curve{S,T},x) where {S<:Space,T<:Number} = c.curve(x)
function tocanonical(c::Curve,x)
    rts=roots(c.curve-x)
    @assert length(rts)==1
    first(rts)
end

fromcanonicalD(c::Curve,x)=differentiate(c.curve)(x)

function indomain(x,c::Curve)
    rts=roots(c.curve-x)
    if length(rts) ≠ 1
        false
    else
        in(first(rts),canonicaldomain(c))
    end
end

reverseorientation(d::IntervalCurve) = IntervalCurve(reverseorientation(d.curve))
reverseorientation(d::PeriodicCurve) = PeriodicCurve(reverseorientation(d.curve))

isambiguous(d::Curve) = ncoefficients(d.curve)==0 && isambiguous(domain(d.curve))
convert(::Type{IntervalCurve{S,T}},::AnyDomain) where {S,T}=Fun(S(AnyDomain()),[NaN])
convert(::Type{PeriodicCurve{S,T}},::AnyDomain) where {S,T}=Fun(S(AnyDomain()),[NaN])

arclength(d::Curve) = linesum(ones(d))


text:
@testset "Each k-mer" begin
    function string_eachkmer(seq::AbstractString, k, step)
        kmers = String[]
        i = 1
        for i in 1:step:length(seq) - k + 1
            subseq = seq[i:i + k - 1]
            if !in('N', subseq)
                push!(kmers, subseq)
            end
        end
        return kmers
    end

    function test_eachkmer(S, seq::AbstractString, k, step)
        xs = [convert(String, x)
              for (i, x) in collect(each(Kmer{eltype(S),k}, S(seq), step))]
        ys = [convert(String, x)
              for (i, x) in collect(eachkmer(S(seq), k, step))]
        zs = string_eachkmer(seq, k, step)
        @test xs == ys == zs
    end

    for k in [0, 1, 3, 16, 32], step in 1:3, len in [1, 2, 3, 5, 10, 100, 1000]
        test_eachkmer(BioSequence{DNAAlphabet{4}}, random_dna(len), k, step)
        test_eachkmer(BioSequence{RNAAlphabet{4}}, random_rna(len), k, step)
        test_eachkmer(ReferenceSequence, random_dna(len), k, step)

        probs = [0.25, 0.25, 0.25, 0.25, 0.00]
        test_eachkmer(BioSequence{DNAAlphabet{2}}, random_dna(len, probs), k, step)
        test_eachkmer(BioSequence{RNAAlphabet{2}}, random_rna(len, probs), k, step)
    end

    @test isempty(collect(each(DNAKmer{1}, dna"")))
    @test isempty(collect(each(DNAKmer{1}, dna"NNNNNNNNNN")))
    @test_throws Exception each(DNAKmer{-1}, dna"ACGT")
    @test_throws Exception each(DNAKmer{33}, dna"ACGT")
    @test collect(each(DNAKmer{3}, dna"AC-TGAG--TGC")) == [(4, DNACodon(DNA_T, DNA_G, DNA_A)),
                                                           (5, DNACodon(DNA_G, DNA_A, DNA_G)),
                                                           (10, DNACodon(DNA_T, DNA_G, DNA_C))]
end


text:
using Turing
using Dates

using ReverseDiffDSGE
using Random
using Zygote
using Distributions
using ComponentArrays
using SharedArrays
using DataFrames
using CSV

using DimensionalData
using PGFPlots
using DataFrames
using KernelDensity
using Statistics
using StatsBase
using LinearAlgebra

# Estimation parameters
nperiods = 100
B = [1,5,6]  # observables
σ_obs = 0.01 # sd of measurement error shocks

## Likelihood related functions
@everywhere logll(θ,ys) = kalman_ll(blanchardkahn(HerbstSchorfheide(θ)),ys,B,σ_obs)

## Generate fake data
θ = rand.(𝒫ₕₛ);
println("True parameter values:")
println("θ = $(round.(θ,digits=2))")

ts(θ,B,σ_obs,nperiods) = ReverseDiffDSGE.timeseries(blanchardkahn(HerbstSchorfheide(θ)),nperiods)[:,B] +
        σ_obs*rand(Normal(),nperiods,length(B)); # ad hoc measurement error shocks
ys = ts(θ,B,σ_obs,nperiods)

θ_axes = getaxes(θ);

## Turing Model
@model function m(y)

        τ⁻¹ ~ InverseGamma(8,8)
        κ   ~ Uniform(0.0,1.0)
        ψ_1 ~ Gamma(4,1/8)
        ψ_2 ~ Gamma(4,1/8)
        ρ_R ~ Uniform(0.5,0.9)
        ρ_g ~ Uniform(0.90,0.99)
        ρ_z ~ Uniform(0.90,0.99)
        σ_R ~ InverseGamma(4,0.32)
        σ_g ~ InverseGamma(4,2.0)
        σ_z ~ InverseGamma(4,0.5)

        ρ = ComponentArray([θ.rA,τ⁻¹,κ,ψ_1,ψ_2,ρ_R,ρ_g,ρ_z,σ_R,σ_g,σ_z],θ_axes)

        @Turing.addlogprob! logll(ρ,y)

end

ch_nuts = sample(m(ys), NUTS{Turing.ZygoteAD}(0.65;max_depth=5,Δ_max=500.0,init_ϵ=0.01), 500)
ch_mh = sample(m(ys), MH(),200_000)

params = (:τ⁻¹,:κ,:ψ_1,:ψ_2,:ρ_R,:ρ_g,:ρ_z,:σ_R,:σ_g,:σ_z)

# irf_true  = impulse_response(blanchardkahn(HerbstSchorfheide(θ)),80)

@dim Θ
@dim ε
@dim σ

function true_irf(θ,periods)
        vars = (:R,:g,:z,:y,:Δy,:π)
        shocks = (:ε_z,:ε_g,:ε_R)
        bkhs = blanchardkahn(HerbstSchorfheide(θ))
        irf_table = zeros(
                        Θ(Val(vars)),
                        Ti(Val(Tuple(1:periods))),
                        ε(Val(shocks)),
                        )
        for k in 1:length(shocks)
                irf_table[:,:,k] = impulse_response(bkhs,periods,k)'
        end
        return irf_table
end

function prior_irf(𝒫,samples,periods)
        vars = (:R,:g,:z,:y,:Δy,:π)
        shocks = (:ε_z,:ε_g,:ε_R)
        irf_table_all = zeros(
                        Θ(Val(vars)),
                        Ti(Val(Tuple(1:periods))),
                        X(Val(Tuple(1:samples))),
                        ε(Val(shocks)),
                        )
        for i in 1:samples
                θ_ = rand.(𝒫)
                bkhs = blanchardkahn(HerbstSchorfheide(θ_))
                for j in 1:length(shocks)
                        irf_table_all[:,:,i,j] = impulse_response(bkhs,periods,j)'
                end
        end
        stats = (:mean,:mode,:median,:p5,:p10,:p90,:p95)
        irf_table = zeros(
                        Θ(Val(vars)),
                        Ti(Val(Tuple(1:periods))),
                        σ(Val(stats)),
                        ε(Val(shocks)),
                        )
        for i in vars, j in 1:periods, k in shocks
                irf_table[i,j,:mean,k] = mean(irf_table_all[i,j,:,k])
                irf_table[i,j,:mode,k] = StatsBase.mode(irf_table_all[i,j,:,k])
                irf_table[i,j,:median,k] = StatsBase.median(irf_table_all[i,j,:,k])
                irf_table[i,j,:p5,k] = percentile(irf_table_all[i,j,:,k],5)
                irf_table[i,j,:p10,k] = percentile(irf_table_all[i,j,:,k],10)
                irf_table[i,j,:p90,k] = percentile(irf_table_all[i,j,:,k],90)
                irf_table[i,j,:p95,k] = percentile(irf_table_all[i,j,:,k],95)
        end
        return irf_table
end

function posterior_irf(chain,samples,periods)
        ndraws,nparams,nchains = size(chain)
        vars = (:R,:g,:z,:y,:Δy,:π)
        shocks = (:ε_z,:ε_g,:ε_R)
        irf_table_all = zeros(
                        Θ(Val(vars)),
                        Ti(Val(Tuple(1:periods))),
                        X(Val(Tuple(1:samples))),
                        ε(Val(shocks)),
                        )
        for i in 1:samples
                θi = (rand(1:ndraws),rand(1:nchains))
                θ_ = ComponentArray([[θ.rA];[chain.value[θi[1],p,θi[2]] for p in params]],θ_axes)
                bkhs = blanchardkahn(HerbstSchorfheide(θ_))
                for j in 1:length(shocks)
                        irf_table_all[:,:,i,j] = impulse_response(bkhs,periods,j)'
                end
        end
        stats = (:mean,:mode,:median,:p5,:p10,:p90,:p95)
        irf_table = zeros(
                        Θ(Val(vars)),
                        Ti(Val(Tuple(1:periods))),
                        σ(Val(stats)),
                        ε(Val(shocks)),
                        )
        for i in vars, j in 1:periods, k in shocks
                irf_table[i,j,:mean,k] = mean(irf_table_all[i,j,:,k])
                irf_table[i,j,:mode,k] = StatsBase.mode(irf_table_all[i,j,:,k])
                irf_table[i,j,:median,k] = StatsBase.median(irf_table_all[i,j,:,k])
                irf_table[i,j,:p5,k] = percentile(irf_table_all[i,j,:,k],5)
                irf_table[i,j,:p10,k] = percentile(irf_table_all[i,j,:,k],10)
                irf_table[i,j,:p90,k] = percentile(irf_table_all[i,j,:,k],90)
                irf_table[i,j,:p95,k] = percentile(irf_table_all[i,j,:,k],95)
        end
        return irf_table
end

nperiods = 40
irf_true      = true_irf(θ,nperiods)
irf_prior     = prior_irf(𝒫ₕₛ,500,nperiods)
irf_posterior_nuts = posterior_irf(ch_nuts,500,nperiods)
irf_posterior_mh = posterior_irf(ch_mh,500,nperiods)

x_labels = (
        y = L"y",
        R = L"R",
        π = L"\pi",
)

ε_labels = (
        ε_z = L"\varepsilon_z",
        ε_g = L"\varepsilon_g",
        ε_R = L"\varepsilon_R",
)

g = GroupPlot(3,3, groupStyle = "horizontal sep = 1.75cm, vertical sep = 1.5cm")

for ε in (:ε_z,:ε_g,:ε_R)
        for x in (:y,:R,:π)
            push!(g,
                PGFPlots.Axis([
                    Plots.Linear([collect(1:nperiods); collect(nperiods:-1:1)],[irf_prior[x,:,:p5,ε];reverse(irf_prior[x,:,:p95,ε])], style="olive!0, fill=black!100, fill opacity=0.1, no marks", closedCycle=true),
                    Plots.Linear(collect(1:nperiods),irf_posterior_nuts[x,:,:p5,ε], style="blue!100, very thick, no marks"),
                    Plots.Linear(collect(1:nperiods),irf_posterior_nuts[x,:,:p95,ε], style="blue!100, very thick, no marks"),
                    Plots.Linear(collect(1:nperiods),irf_posterior_mh[x,:,:p5,ε], style="red!100, dotted, very thick, no marks"),
                    Plots.Linear(collect(1:nperiods),irf_posterior_mh[x,:,:p95,ε], style="red!100, dotted, very thick, no marks"),
                    Plots.Linear(irf_true[x,:,ε], style="black, dashed, very thick, no marks"),
                ],
                title=L"%$(x_labels[x]) response to %$(ε_labels[ε])",
                ymin=minimum([irf_posterior_nuts[x,:,:p5,ε];irf_posterior_mh[x,:,:p5,ε];irf_true[x,:,ε]]),
                ymax=maximum([irf_posterior_nuts[x,:,:p95,ε];irf_posterior_mh[x,:,:p95,ε];irf_true[x,:,ε]]),
                xmin=0,
                xmax=(ε == :ε_R ? 8 : 40),
                width="5cm",
                height="5cm"))
        end
end

save("src/images/irfs_standalone.tex",g)
save("src/images/irfs.tex",g, include_preamble=false)

post_nuts = (;zip(params,(kde(ch_nuts[:,p,1].data) for p in params))...)
post_mh = (;zip(params,(kde(ch_mh[:,p,1].data) for p in params))...)

labels = (
        τ⁻¹ = L"\tau^{-1}",
        κ   = L"\kappa",
        ψ_1 = L"\psi_1",
        ψ_2 = L"\psi_2",
        ρ_R = L"\rho_R",
        ρ_g = L"\rho_G",
        ρ_z = L"\rho_z",
        σ_R = L"\sigma_R",
        σ_g = L"\sigma_g",
        σ_z = L"\sigma_z",
)

function support(var,percentile)
    upper = Distributions.percentile(𝒫ₕₛ[var],50 + percentile/2)
    lower = Distributions.percentile(𝒫ₕₛ[var],50 - percentile/2)
    return (lower,upper)
end

g = GroupPlot(3,4, groupStyle = "horizontal sep = 1.75cm, vertical sep = 1.5cm")

for v in params
    push!(g,
        PGFPlots.Axis([
            Plots.Linear(x -> pdf(𝒫ₕₛ[v],x), support(v,99.5), style="black!10, fill=black!10", closedCycle=true),
            Plots.Linear(x -> pdf(post_nuts[v],x),support(v,99.9), style="blue, thick"),
            Plots.Linear(x -> pdf(post_mh[v],x),support(v,99.9), style="red, dotted, thick"),
            Plots.Linear([θ[v],θ[v]],[0.0,1.5*maximum(post_nuts[v].density)], style="black, no marks, very thick"),
        ],ymin=0,title=labels[v],width="5cm", height="5cm"))
end

save("src/images/posteriors_standalone.tex",g)
save("src/images/posteriors.tex",g, include_preamble=false)


text:
function MCTS.estimate_value(est::BasicPOMCP.SolvedFORollout,
                             bmdp::GenerativeBeliefMDP,
                             belief,
                             d::Int64)
    sim = RolloutSimulator(est.rng, Nullable{Any}(), Nullable{Float64}(), Nullable(d))
    return simulate(sim, bmdp.pomdp, est.policy, rand(est.rng, belief))
end

function MCTS.estimate_value(est::BasicPOMCP.SolvedFOValue{P},
                             bmdp::MeanRewardBeliefMDP,
                             b::AbstractParticleBelief,
                             d::Int64) where P <: ValueIterationPolicy
    r = [estimate_value(est, bmdp.pomdp, s, d) for s in particles(b)]
    mean(r, weights(weights(b)))
end

function MCTS.estimate_value(estimator::BasicPOMCP.SolvedFOValue{P}, pomdp::POMDP, start_state, steps::Int) where P <: ValueIterationPolicy
    POMDPs.value(estimator.policy, start_state)
end

function MCTS.estimate_value(est::BasicPOMCP.SolvedFOValue{P},
                             bmdp::MeanRewardBeliefMDP,
                             b::AbstractParticleBelief,
                             d::Int64) where P <: AlphaVectorPolicy
    estimate_value(est, bmdp.pomdp, b, d)
end

function MCTS.estimate_value(estimator::BasicPOMCP.SolvedFOValue{P}, pomdp::POMDP, b, steps::Int) where P <: AlphaVectorPolicy
    POMDPs.value(estimator.policy, b)
end

Random.seed!(planner::Policy, x) = planner


text:
# # QCDNUM example program
#
# Based on `example.f` detailed [here](https://www.nikhef.nl/~h24/qcdnum-files/testjobs1701/example.f).

using QCDNUM
using Printf

# Let's define some useful inputs

xmin = Float64.([1.0e-4])
iwt = Int32.([1])
ng = 1
nxin = 100
iosp = 3
nx = 10

qq = Float64.([2e0, 1e4])
wt = Float64.([1e0, 1e0])
nq = 1
nqin = 60
ngq = 2
itype = 1

as0 = 0.364
r20 = 2.0

q2c = 3.0
q2b = 25.0
q0 = 2.0
iqt = 999

def = Float64.([0., 0., 0., 0., 0.,-1., 0., 1., 0., 0., 0., 0., 0.,
      0., 0., 0., 0.,-1., 0., 0., 0., 1., 0., 0., 0., 0.,
      0., 0., 0.,-1., 0., 0., 0., 0., 0., 1., 0., 0., 0.,
      0., 0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0.,-1., 0., 0., 0., 0., 0., 0., 0., 1., 0., 0.,
      0., 0., 1., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.,
      0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.]);

nfin = 0
x = 1.0e-3
q = 1.0e3
pdf = Array{Float64}(undef, 13)
qmz2 = 8315.25;

# Now, we can define the input PDF function for the PDF we would like to evolve. We also then have to use the
# Julia `@cfunction` macro so that we can pass it to QCDNUM through the wrapper.

function func(ipdf, x)::Float64
    i = ipdf[]
    xb = x[]
    adbar = 0.1939875
    f = 0
    if (i == 0)
        ag = 1.7
        f = ag * xb^-0.1 * (1.0-xb)^5.0
    end
    if (i == 1)
        ad = 3.064320
        f = ad * xb^0.8 * (1.0-xb)^4.0
    end
    if (i == 2)
        au = 5.107200
        f = au * xb^0.8 * (1.0-xb)^3.0
    end
    if (i == 3)
        f = 0.0
    end
    if (i == 4)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0
    end
    if (i == 5)
        f = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
    end
    if (i == 6)
        xdbar = adbar * xb^-0.1 * (1.0-xb)^6.0
        xubar = adbar * xb^-0.1 * (1.0-xb)^6.0 * (1.0-xb)
        f = 0.2 * (xdbar + xubar)
    end
    if (i == 7)
        f = 0.0
    end
    if (i == 8)
        f = 0.0
    end
    if (i == 9)
        f = 0.0
    end
    if (i == 10)
        f = 0.0
    end
    if (i == 11)
        f = 0.0
    end
    if (i == 12)
        f = 0.0
    end
    return f
end

func_c = @cfunction(func, Float64, (Ref{Int32}, Ref{Float64}))

# Next, let's initialise QCDNUM and run QCDNUM to evolve the PDF.

QCDNUM.qcinit(-6, " ")
nx = QCDNUM.gxmake(xmin, iwt, ng, nxin, iosp)
nq = QCDNUM.gqmake(qq, wt, ngq, nqin)
nw = QCDNUM.fillwt(itype)
QCDNUM.setord(3)
QCDNUM.setalf(as0, r20)
iqc = QCDNUM.iqfrmq(q2c)
iqb = QCDNUM.iqfrmq(q2b)
iqt = QCDNUM.iqfrmq(1e11)
QCDNUM.setcbt(0, iqc, iqb, 0)
iq0 = QCDNUM.iqfrmq(q0)
eps = QCDNUM.evolfg(itype, func_c, def, iq0)
pdf = QCDNUM.allfxq(itype, x, q, 0, 1)
asmz, a, b = QCDNUM.asfunc(qmz2)
csea = 2*pdf[3];

# Now, we can briefly check the results

@printf("x, q, CharmSea = %0.4e, %0.4e, %0.4e\n", x, q, csea)
@printf("as(mz2) = %0.4e", asmz)


text:
# using Base.Test
# using MIToS.Pfam
# using MIToS.PDB
# using MIToS.Information
# using PairwiseListMatrices
using ROCAnalysis

print("""

Test download from Pfam
=======================
""")

let pfam_code = "PF11591"
    @test_throws ErrorException downloadpfam("2vqc")
    filename = downloadpfam(pfam_code)
    try
        aln = read(filename, Stockholm)
        if size(aln) == (6,34)
            @test getannotfile(aln, "ID") == "2Fe-2S_Ferredox"
        end
    finally
        rm(filename)
    end
end

print("""

Test PDB code from Pfam
=======================
""")

@test getseq2pdb(read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm))["F112_SSV1/3-112"] == [("2VQC","A")]

print("""

Test Mapping PDB/Pfam
=====================
""")

let msa = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, generatemapping=true, useidcoordinates=true),
    map = msacolumn2pdbresidue(msa, "F112_SSV1/3-112", "2VQC", "A", "PF09645", ascii(joinpath(pwd(), "data", "2vqc.xml.gz"))),
    res = residuesdict(read(joinpath(pwd(), "data", "2VQC.xml"), PDBML), "1", "A", "ATOM", "*")

    #     -45              20 pdb
    #.....QTLNSYKMAEIMYKILEK  msa seq
    #     123456789012345678  msa col
    #     345678901234567890  uniprot 3-20
    #    ****              *
    #12345678901234567890123  ColMap

    @test_throws KeyError map[5]  # insert
    @test map[6]  == ""           # missing
    @test map[7]  == "4"
    @test map[8]  == "5"
    @test map[23] == "20"

    #.....QTLNSYKMAEIMYKILEKKGELTLEDILAQFEISVPSAYNIQRALKAICERHPDECEVQYKNRKTTFKWIKQEQKEEQKQEQTQDNIAKIFDAQPANFEQTDQGFIKAKQ..... msa seq
    #.....X---HHHHHHHHHHHHHHHSEE-HHHHHHHH---HHHHHHHHHHHHHHHHH-TTTEEEEE-SS-EEEEE--XXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX..... pdb ss
    #                                                                                111111111111111111111 ColMap hundreds
    #         111111111122222222223333333333444444444455555555556666666666777777777788888888889999999999000000000011111111112 ColMap tens
    #123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890123456789012345678901234567890 ColMap ones
    #                                                                           **                                     **

    @test_throws KeyError map[116]   # insert
    @test map[115] == ""             # missing
    @test map[77]  == ""             # missing
    @test map[76]  == "73"

    print("""

  Test Residues
  -------------
  """)

    msares = msaresidues(msa, res, map)

    @test_throws KeyError msares[5]   # insert
    @test_throws KeyError msares[6]   # missing
    @test msares[7].id.name  == "THR" # T
    @test msares[8].id.name  == "LEU" # L
    @test msares[23].id.name == "LYS" # K

    @test_throws KeyError msares[116]   # insert
    @test_throws KeyError msares[117]   # missing
    @test_throws KeyError msares[77]    # missing
    @test msares[76].id.name == "LYS"   # K

end

print("""

Test Contacts
=============
""")

let msa = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, generatemapping=true, useidcoordinates=true),
    map = msacolumn2pdbresidue(msa, "F112_SSV1/3-112", "2VQC", "A", "PF09645", ascii(joinpath(pwd(), "data", "2vqc.xml.gz"))),
    res = residuesdict(read(joinpath(pwd(), "data", "2VQC.xml"), PDBML), "1", "A", "ATOM", "*")

    @test msacolumn2pdbresidue(msa, "F112_SSV1/3-112", "2VQC", "A", "PF09645", ascii(joinpath(pwd(), "data", "2vqc.xml.gz")), strict=true, checkpdbname=true) == map

    @test length(res) == 70
    @test length(unique(values(map))) == 71

    #     -45              20 pdb
    #.....QTLNSYKMAEIMYKILEK  msa seq
    #     123456789012345678  msa col
    #12345678901234567890123  ColMap
    #     345678901234567890  uniprot 3-20
    #     ***              *

    @test_throws KeyError res["3"]
    @test res[map[7]].id.number == "4"

    contacts = msacontacts(msa, res, map, 6.05)
    missings = sum(isnan(contacts), 1)

    @test size(contacts) == (110,110)

    @test missings[1]  == 110              # missing
    @test missings[2]  == (110 - (70 - 1)) # 70 residues, 1 diagonal NaN
    @test missings[3]  == (110 - (70 - 1))
    @test missings[18] == (110 - (70 - 1))

    @test missings[110] == 110                 # missing
    @test missings[72]  == 110                 # missing
    @test missings[71]  == (110 - (70 - 1))

    ncontacts = sum(contacts .== 1.0, 1)

    @test ncontacts[1]  == 0
    @test ncontacts[2]  == 2
    @test ncontacts[3]  == 6

    print("""

  Test MSA contact map using PDBResidues from the MSA
  ---------------------------------------------------
  """)

    println("")

    msares = msaresidues(msa, res, map)

    ncol = ncolumns(msa)
    colmap = getcolumnmapping(msa)

    for i in 1:(ncol-1), j in (i+1):ncol
        if (colmap[i] in keys(msares)) && (colmap[j] in keys(msares))
            print(".")
            @test (contacts[i,j] .== 1.0) == contact(msares[ colmap[i] ], msares[ colmap[j] ], 6.05)
        end
    end

    println("")

    print("""

  Test hasresidues
  """)

    mask = hasresidues(msa, map)

    @test mask[1] == false
    @test mask[2] == true
    @test sum(mask) == length(msares)

end

print("""

Test AUC and Contact Masks
==========================
""")

let ntru = 90,
    nfal = 100,
    score_tru = Float16[ 2 + 2x for x in randn(ntru)],
    score_fal = Float16[-2 + 2x for x in randn(nfal)],
    msacontacts = PairwiseListMatrix(vcat(ones(Float16, ntru), zeros(Float16, nfal)), Int16[x for x in 1:20], false),
    score = PairwiseListMatrix(vcat(score_tru, score_fal), Int16[x for x in 1:20], false),
    correct = 1 - auc(roc(score_tru, score_fal))

    @test AUC(score, msacontacts) == correct

end

print("""

AUC Example using MI
--------------------
""")

let msa = read(joinpath(pwd(), "data", "PF09645_full.stockholm"), Stockholm, generatemapping=true, useidcoordinates=true),
    map = msacolumn2pdbresidue(msa, "F112_SSV1/3-112", "2VQC", "A", "PF09645", ascii(joinpath(pwd(), "data", "2vqc.xml.gz"))),
    res = residuesdict(read(joinpath(pwd(), "data", "2VQC.xml"), PDBML), "1", "A", "ATOM", "*"),
    contacts = msacontacts(msa, res, map)

    @test round( AUC(buslje09(msa, lambda=0.05, threshold=62.0, samples=0)[2], contacts) , 4) == 0.5291
end


text:
"""
    algorithm971(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="", colsumlist::AbstractString="", dim::Number=3, noversamples::Number=5, niter::Number=3, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)

Algorithm 971, which is one of randomized SVD algorithm.

Input Arguments
---------
- `input` : Julia Binary file generated by `OnlinePCA.csv2bin` function.
- `outdir` : The directory specified the directory you want to save the result.
- `scale` : {log,ftt,raw}-scaling of the value.
- `pseudocount` : The number specified to avoid NaN by log10(0) and used when `Feature_LogMeans.csv` <log10(mean+pseudocount) value of each feature> is generated.
- `rowmeanlist` : The mean of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `rowvarlist` : The variance of each row of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `colsumlist` : The sum of counts of each columns of matrix. The CSV file is generated by `OnlinePCA.sumr` functions.
- `dim` : The number of dimension of PCA.
- `noversamples` : The number of over-sampling.
- `niter` : The number of power interation.
- `initW` : The CSV file saving the initial values of eigenvectors.
- `initV` : The CSV file saving the initial values of loadings.
- `logdir` : The directory where intermediate files are saved, in every evalfreq (e.g. 5000) iteration.
- `perm` : Whether the data matrix is shuffled at random.
- `cper` : Count per X (e.g. CPM: Count per million <1e+6>)

Output Arguments
---------
- `V` : Eigen vectors of covariance matrix (No. columns of the data matrix × dim)
- `λ` : Eigen values (dim × dim)
- `U` : Loading vectors of covariance matrix (No. rows of the data matrix × dim)
- `Scores` : Principal component scores
- `ExpVar` : Explained variance by the eigenvectors
- `TotalVar` : Total variance of the data matrix
"""
function algorithm971(;input::AbstractString="", outdir::Union{Nothing,AbstractString}=nothing, scale::AbstractString="ftt", pseudocount::Number=1f0, rowmeanlist::AbstractString="", rowvarlist::AbstractString="", colsumlist::AbstractString="", dim::Number=3, noversamples::Number=5, niter::Number=3, initW::Union{Nothing,AbstractString}=nothing, initV::Union{Nothing,AbstractString}=nothing, logdir::Union{Nothing,AbstractString}=nothing, perm::Bool=false, cper::Number=1f0)
    # Initial Setting
    pca = ALGORITHM971()
    pseudocount, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar = init(input, pseudocount, dim, rowmeanlist, rowvarlist, colsumlist, initW, initV, logdir, pca, cper, scale)
    # Perform PCA
    out = algorithm971(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, noversamples, niter, logdir, pca, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, perm, cper)
    # Output
    if outdir isa String
        writecsv(joinpath(outdir, "Eigen_vectors.csv"), out[1])
        writecsv(joinpath(outdir, "Eigen_values.csv"), out[2])
        writecsv(joinpath(outdir, "Loadings.csv"), out[3])
        writecsv(joinpath(outdir, "Scores.csv"), out[4])
        writecsv(joinpath(outdir, "ExpVar.csv"), out[5])
        writecsv(joinpath(outdir, "TotalVar.csv"), out[6])
    end
    return out
end

function algorithm971(input, outdir, scale, pseudocount, rowmeanlist, rowvarlist, colsumlist, dim, noversamples, niter, logdir, pca, W, D, rowmeanvec, rowvarvec, colsumvec, N, M, TotalVar, perm, cper)
    N, M = nm(input)
    tmpN = zeros(UInt32, 1)
    tmpM = zeros(UInt32, 1)
    x = zeros(UInt32, M)
    normx = zeros(Float32, M)
    l = dim + noversamples
    @assert 0 < dim ≤ l ≤ min(N, M)
    Ω = rand(Float32, M, l)
    Y = zeros(Float32, N, l)
    L = zeros(Float32, N, l)
    Q = zeros(Float32, N, l)
    B = zeros(Float32, l, M)
    G = zeros(Float32, M, l)
    # If not 0 the calculation is converged
    n = 1
    # Each epoch s
    println("Random Projection : Y = A Ω")
    progress = Progress(N)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        # Each step n
        while(n <= N)
            next!(progress)
            # Row vector of data matrix
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            # Random Projection
            Y[n,:] .= (normx'*Ω)[1,:]
            n += 1
        end
        close(stream)
    end

    if niter > 0
        # LU factorization
        println("LU factorization : L = lu(Y)")
        L .= lu!(Y).L
        for i in 1:niter
            println("Normalized power iterations (1/3) : A' L")
            n = 1
            AtL = zeros(Float32, M, l)
            progress = Progress(N)
            open(input) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                while(n <= N)
                    next!(progress)
                    # Row vector of data matrix
                    read!(stream, x)
                    normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
                    if perm
                        normx .= normx[randperm(length(normx))]
                    end
                    AtL .+= normx*L[n,:]'
                    n += 1
                end
                close(stream)
            end
            println("L = lu(A' L)")
            G .= lu!(AtL).L

            println("Normalized power iterations (2/3) : Y = A lu(A' L)")
            n = 1
            progress = Progress(N)
            open(input) do file
                stream = ZstdDecompressorStream(file)
                read!(stream, tmpN)
                read!(stream, tmpM)
                # Each step n
                while(n <= N)
                    next!(progress)
                    # Row vector of data matrix
                    read!(stream, x)
                    normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
                    if perm
                        normx .= normx[randperm(length(normx))]
                    end
                    Y[n,:] .= (normx'*G)[1,:]
                    n += 1
                end
                close(stream)
            end
            if i < niter
                println("Normalized power iterations (3/3) : L = lu(A lu(A' L))")
                # Renormalize with LU factorization
                L .= lu!(Y).L
            else
                println("QR factorization  (3/3) : Q = qr(A lu(A' L))")
                # Renormalize with QR factorization
                Q .= Array(qr!(Y).Q)
            end
        end
    else
        println("QR factorization : Q = qr(Y)")
        # Renormalize with QR factorization
        Q .= Array(qr!(Y).Q)
    end

    println("Calculation of small matrix : B = Q' A")
    n = 1
    progress = Progress(N)
    open(input) do file
        stream = ZstdDecompressorStream(file)
        read!(stream, tmpN)
        read!(stream, tmpM)
        # Each step n
        while(n <= N)
            next!(progress)
            # Row vector of data matrix
            read!(stream, x)
            normx = normalizex(x, n, stream, scale, pseudocount, rowmeanlist, rowmeanvec, rowvarlist, rowvarvec, colsumlist, colsumvec, cper)
            if perm
                normx .= normx[randperm(length(normx))]
            end
            B .+= Q[n,:]*normx'
            n += 1
        end
        close(stream)
    end
    # SVD with small matrix
    println("SVD with small matrix : svd(B)")
    W, σ, V = svd(B)
    U = Q*W
    λ = σ .* σ ./ M
    # PC scores, Explained Variance
    Scores = zeros(Float32, M, dim)
    for n = 1:dim
        Scores[:, n] .= λ[n] .* V[:, n]
    end
    ExpVar = sum(λ) / TotalVar
    # Return
    return (V[:,1:dim], λ[1:dim], U[:,1:dim], Scores[:,1:dim], ExpVar, TotalVar)
end


text:
# File for generating data for Figure 5

using Pkg
using Random, Distributions
using Statistics
using ProgressMeter
using JLD2

include("/home/christian/UniMA/EMMA/src/RFR.jl")
include("/home/christian/UniMA/EMMA/src/cross_val.jl")
include("/home/christian/UniMA/EMMA/src/aux_functions.jl")

# Function for getting the best cv-parameter
function get_best_model(n::Int, d::Int, type::String, σ::Number, m_features::Int, parameter::Symbol)

    xtrain1, xtest1, ytrain1, ytest1 = make_data(n, d, type, σ)

    a_list = floor.(collect(LinRange(0, 50, 21)))

    d1 = Dict{Symbol, Vector{Float64}}(
        :max_features => [m_features],
        :n_trees => [30])

    if parameter == :CART
        d1[:α] = [0.0]
    else
        d1[parameter] = a_list
    end

    cv = cross_val(d1, random_state = 0)
    fit!(cv, xtrain1, ytrain1, nfolds=3)

    return(best_model(cv, "mse"))
end

function get_data(n_runs, max_d, σ, parameter::Symbol)

    # Allocate space for results
    res_mat = zeros(max_d, n_runs)

    # Setting parameters
    n = 2000
    d = 10
    m_features = d
    # m_features = Int(floor(d/3))

    # Getting the best model given error term variance
    rf = get_best_model(n, d, "friedman", σ, m_features, parameter)
    println(rf.α, rf.min_samples_leaf)

    @showprogress for r in 1:n_runs

        x_train, x_test, y_train, y_test = make_data(n, d, "friedman", σ)

        fit!(rf, x_train, y_train)

        result_arr = Array{Float64}(undef, 0, 2)

        for i in 1:length(rf.trees)

            # Save depth and splitting dimension
            res = zeros(length(rf.trees[i].depth_list), 2)
            res[:,1] = rf.trees[i].depth_list
            res[:,2] = rf.trees[i].split_dimensions

            result_arr = vcat(result_arr, res)
        end

        # Only comparing non-terminal nodes
        plot_data = result_arr[result_arr[:,2].!=-2, :]
        # Noisy dimensions are all regressors above 5
        plot_data[:,2] = plot_data[:,2].>5

        # Save figures as dataframe
        plot_data = DataFrame(plot_data, :auto)
        rename!(plot_data, ["depth", "noise"])

        # Grouping them by depth
        gdf = groupby(plot_data, :depth)
        # Getting mean of fraction split on noise
        cdf2 = combine(gdf, :noise => mean)

        res_mat[1:size(cdf2, 1), r] = cdf2[!, :noise_mean]
    end
    return res_mat, rf
end

# Set seed for reproducability
Random.seed!(68151)

n_runs=100
max_d=50 #Max depth allowed - need to set this parameter; max acieved is around 30

# Change σ here accordingly {1, 3, 8} to get figures

# Getting data
res_mat1, rf_α = get_data(n_runs, max_d, 8, :α)
res_mat1_const, rf_const = get_data(n_runs, max_d, 8, :CART)

out_list = [res_mat1, res_mat1_const]

jldsave("data/wide_sigma8.jld2"; out_list)

max_stars_repo_name: christianhilscher/emma

max_issues_repo_name: christianhilscher/emma

max_forks_repo_name: christianhilscher/emma


text:
# This file was generated, do not modify it. # hide
using MLJ, PrettyPrinting, DataFrames, Statistics, CSV, Dates
using PyPlot, HTTP
MLJ.color_off() # hide

req = HTTP.get("https://raw.githubusercontent.com/bbrandom91/KC_Housing/master/kc_house_data.csv")

df = CSV.read(req.body, missingstring="NA")
describe(df)

max_stars_repo_name: ven-k/MLJTutorials

max_issues_repo_name: ven-k/MLJTutorials

max_forks_repo_name: ven-k/MLJTutorials