text:
#' ---
#' title: Integer Programming
#' ---

#' **Originally Contributed by**: Arpit Bhatia

#' While we already know how to set a variable as integer or binary in the `@variable` macro,
#' this tutorial covers other JuMP features for integer programming along with some modelling techniques.

using JuMP

#' ## Modelling Logical Conditions
#' Generally, in a mathematical programming problem, all constraints must hold.
#' However, we might want to have conditions where we have some logical conditions between constraints.
#' In such cases, we can use binary variables for modelling logical conditions between constraints.

#' ### Disjunctive Constraints (OR)
#' Suppose that we are given two constraints $a'x \geq b$ and $c' x \geq d$,
#' in which all components of $a$ and $c$ are non-negative.
#' We would like to model a requirement that at least one of the two constraints is satisfied.
#' For this, we defined a binary variable $y$ and impose the constraints:

#' $$
#' \begin{align*}
#' a' x \geq y b \\
#' c' x \geq (1 - y) d \\
#' y \in \{0,1\}
#' \end{align*}
#' $$

a = rand(1:100, 5, 5)
c = rand(1:100, 5, 5)
b = rand(1:100, 5)
d = rand(1:100, 5)

model = Model()
@variable(model, x[1:5])
@variable(model, y, Bin)
@constraint(model, a * x .>= y .* b)
@constraint(model, c * x .>= (1 - y) .* d);

#' ### Conditional Constraints ($\implies$)
#' Suppose we want to model that a certain linear inequality must be satisfied when some other event occurs.
#' i.e. for a binary variable $z$, we want to model the implication

#' $$
#' \begin{align*}
#' z = 1 \implies a^Tx\leq b
#' \end{align*}
#' $$

#' If we know in advance an upper bound $a^Tx\leq b$. Then we can write the above as a linear inequality

#' $$
#' \begin{align*}
#' a^Tx\leq b + M(1-z)
#' \end{align*}
#' $$

a = rand(1:100, 5, 5)
b = rand(1:100, 5)
m = rand(10000:11000, 5)

model = Model()
@variable(model, x[1:5])
@variable(model, z, Bin)
@constraint(model, a * x .<=  b .+ (m .* (1 - z)));
# If z was a regular Julia variable, we would not have had to use the vectorized dot operator

#' ### Boolean Operators on Binary Variables
#' The following table is useful when we want to model boolean operators in the form of
#' linear inequalities that can be given to a solver.

#' | Boolean Expression | Constraint                           |
#' |:----------         |                           ----------:|
#' | $z=x \lor y$       | $x \leq z,  y \leq z,  z \leq x+y$   |
#' | $z=x \land y$      | $x \geq z,  y \geq z,  z+1 \geq x+y$ |
#' | $z= \neg x$        | $z = 1 − x$                          |
#' | $x \implies y$     | $x \leq y$                           |
#' | $x \iff y$         | $x = y$                              |

#' ## Modelling Integer Variables

#' ### Integer Variables using Constraints
#' We can add binary and integer restrictions to the domain of each variable using the `@constraint` macro as well.

model = Model()

@variable(model, x)
@variable(model, y)
@constraint(model, x in MOI.ZeroOne())
@constraint(model, y in MOI.Integer());

#' ### Semi-Continuous Variables
#' A semi-continuous variable is a continuous variable
#' between bounds $[l,u]$ that also can assume the value zero. ie.
#' $$
#' x \in \{0\} \cup \{l,u\}
#' $$

l = 7.45
u = 22.22
@variable(model, a)
@constraint(model, a in MOI.Semicontinuous(l, u))

#' ### Semi-Integer Variables
#' A semi-integer variable is a variable which asummes integer values
#' between bounds $[l,u]$ and can also assume the value zero. ie.

#' $$
#' x \in \{0\} \cup (\{l,u\} \cap \mathbb{Z})
#' $$

l = 5
u = 34
@variable(model, b)
@constraint(model, b in MOI.Semiinteger(l, u))

#' Note that the bounds specified in `MOI.Semiinteger` must be integral otherwise it would throw an error.

#' ## Special Ordered Sets

#' ### Special Ordered Sets of type 1 (SOS1)
#' A Special Ordered Set of type 1 is a set of variables,
#' at most one of which can take a non-zero value, all others being at 0.
#' They most frequently apply where a set of variables are actually 0-1 variables:
#' in other words, we have to choose at most one from a set of possibilities.

@variable(model, u[1:3])
@constraint(model, u in MOI.SOS1([1.0, 2.0, 3.0]))

#' Note that we have to pass MOI.SOS1 a weight vector which is essentially an ordering on the variables.
#' If the decision variables are related and have a physical ordering, then the weight vector,
#' although not used directly in the constraint, can help the solver make a better decision in the solution process.

#' ### Special Ordered Sets of type 2 (SOS2)

#' A Special Ordered Set of type 2 is a set of non-negative variables,
#' of which at most two can be non-zero, and if two are non-zero these must be consecutive in their ordering.

@variable(model, v[1:3])
@constraint(model, v in MOI.SOS2([3.0, 1.0, 2.0]))

#' The ordering provided by the weight vector is more important in this case as
#' the variables need to be consecutive according to the ordering.
#' For example, in the above constraint, the possible pairs are:
#' * Consecutive
#'   * (`x[1]` and `x[3]`) as they correspond to 3 and 2 resp. and thus can be non-zero
#'   * (`x[2]` and `x[3]`) as they correspond to 1 and 2 resp. and thus can be non-zero
#' * Non-consecutive
#'   * (`x[1]` and `x[2]`) as they correspond to 3 and 1 resp. and thus cannot be non-zero


text:
################################################################################
#
#  Parent
#
################################################################################

parent(a::AlgAssElem) = a.parent

################################################################################
#
#  Special elements
#
################################################################################

zero(A::AlgAss) = A()

one(A::AlgAss) = A(A.one)

################################################################################
#
#  Unary operations
#
################################################################################

function -(a::AlgAssElem{T}) where {T}
  coeffs = [ -a.coeffs[i] for i = 1:dim(parent(a)) ]
  return AlgAssElem{T}(parent(a), coeffs)
end

################################################################################
#
#  Binary operations
#
################################################################################

function +(a::AlgAssElem{T}, b::AlgAssElem{T}) where {T}
  parent(a) != parent(b) && error("Parents don't match.")
  coeffs = Array{T, 1}(dim(parent(a)))
  for i = 1:dim(parent(a))
    coeffs[i] = a.coeffs[i] + b.coeffs[i]
  end
  return AlgAssElem{T}(parent(a), coeffs)
end

function -(a::AlgAssElem{T}, b::AlgAssElem{T}) where {T}
  parent(a) != parent(b) && error("Parents don't match.")
  coeffs = Array{T, 1}(dim(parent(a)))
  for i = 1:dim(parent(a))
    coeffs[i] = a.coeffs[i] - b.coeffs[i]
  end
  return AlgAssElem{T}(parent(a), coeffs)
end

function *(a::AlgAssElem{T}, b::AlgAssElem{T}) where {T}
  parent(a) != parent(b) && error("Parents don't match.")

  A = parent(a)
  n = dim(A)
  c = A()
  t = base_ring(A)()
  for i = 1:n
    if iszero(a.coeffs[i])
      continue
    end
    for j = 1:n
      t = a.coeffs[i]*b.coeffs[j]
      if iszero(t)
        continue
      end
      for k = 1:n
        c.coeffs[k] += A.mult_table[i, j, k]*t
      end
    end
  end
  return c
end

function mul!(c::AlgAssElem{T}, a::AlgAssElem{T}, b::AlgAssElem{T}) where {T}
  A = parent(a)
  n = dim(A)
  t = base_ring(A)()
  s = base_ring(A)()

  if c === a || c === b
    z = parent(a)()
    mul!(z, a, b)
    return z
  end

  for k in 1:n
    c.coeffs[k] = zero!(c.coeffs[k])
  end

  for i = 1:n
    if iszero(a.coeffs[i])
      continue
    end
    for j = 1:n
      t = a.coeffs[i]*b.coeffs[j]
      if iszero(t)
        continue
      end
      for k = 1:n
        s = mul!(s, A.mult_table[i, j, k], t)
        c.coeffs[k] = add!(c.coeffs[k], c.coeffs[k], s)
        #c.coeffs[k] += A.mult_table[i, j, k]*t
      end
    end
  end
  #@assert c == a * b
  return c
end

################################################################################
#
#  Ad hoc operations
#
################################################################################

function *(a::AlgAssElem{T}, b::T) where { T <: RingElem }
  return AlgAssElem{T}(parent(a), a.coeffs.*b)
end

*(b::T, a::AlgAssElem{T}) where { T <: RingElem } = a*b

*(a::AlgAssElem{T}, b::Union{Integer, fmpz}) where {T} = a*base_ring(parent(a))(b)

*(b::Union{Integer, fmpz}, a::AlgAssElem{T}) where {T} = a*b

dot(a::AlgAssElem{T}, b::T) where {T} = a*b

dot(b::T, a::AlgAssElem{T}) where {T} = b*a

dot(a::AlgAssElem{T}, b::Union{Integer, fmpz}) where {T} = a*b

dot(b::Union{Integer, fmpz}, a::AlgAssElem{T}) where {T} = b*a

################################################################################
#
#  Exponentiation
#
################################################################################

function ^(a::AlgAssElem, b::Int)
  if b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  else
    if b < 0
      error("Not implemented yet")
      #a = inv(a)
      #b = -b
    end
    bit = ~((~UInt(0)) >> 1)
    while (UInt(bit) & b) == 0
      bit >>= 1
    end
    z = deepcopy(a)
    bit >>= 1
    while bit != 0
      z = mul!(z, z, z)
      if (UInt(bit) & b) != 0
        z = mul!(z, z, a)
      end
      bit >>= 1
    end
    return z
  end
end

function ^(a::AlgAssElem, b::fmpz)
  if nbits(b) < 64
    return a^Int(b)
  end
  if b < 0
    error("Not implemented yet")
  elseif b == 0
    return one(parent(a))
  elseif b == 1
    return deepcopy(a)
  elseif mod(b, 2) == 0
    c = a^(div(b, 2))
    return c*c
  elseif mod(b, 2) == 1
    return a^(b - 1)*a
  end
end

################################################################################
#
#  Parent object overloading
#
################################################################################

(A::AlgAss{T})() where {T} = AlgAssElem{T}(A)

function (A::AlgAss{T})(c::Array{T, 1}) where {T}
  length(c) != dim(A) && error("Dimensions don't match.")
  return AlgAssElem{T}(A, c)
end

function Base.getindex(A::AlgAss{T}, i::Int) where {T}
  (i < 1 || i > dim(A)) && error("Index must be in range $(1:dim(A))")
  n = dim(A)
  v = Vector{T}(n)
  for j in 1:n
    v[j] = zero(base_ring(A))
  end
  v[i] = one(base_ring(A))
  return A(v)
end

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

function show(io::IO, a::AlgAssElem)
  print(io, "Element of ")
  print(io, parent(a))
  print(io, " with coefficients ")
  print(io, a.coeffs)
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AlgAssElem{T}, dict::ObjectIdDict) where {T}
  b = parent(a)()
  for x in fieldnames(a)
    if x != :parent && isdefined(a, x)
      setfield!(b, x, Base.deepcopy_internal(getfield(a, x), dict))
    end
  end
  return b
end

################################################################################
#
#  Equality
#
################################################################################

function ==(a::AlgAssElem{T}, b::AlgAssElem{T}) where {T}
  parent(a) != parent(b) && return false
  return a.coeffs == b.coeffs
end

################################################################################
#
#  Minpoly
#
################################################################################

function Generic.minpoly(a::AlgAssElem)
  M = representation_matrix(a)
  R = PolynomialRing(base_ring(parent(a)), "x", cached=false)[1]
  return minpoly(R, M)
end

################################################################################
#
#  Representation matrix
#
################################################################################

function elem_to_mat_row!(M::MatElem{T}, i::Int, a::AlgAssElem{T}) where T
  for c = 1:cols(M)
    M[i, c] = deepcopy(a.coeffs[c])
  end
  return nothing
end

function elem_from_mat_row(A::AlgAss{T}, M::MatElem{T}, i::Int) where T
  a = A()
  for c = 1:cols(M)
    a.coeffs[c] = deepcopy(M[i, c])
  end
  return a
end

function representation_matrix(a::AlgAssElem, action::Symbol=:left)
  A = parent(a)
  M = zero_matrix(base_ring(A), dim(A), dim(A))
  if action==:left
    for i = 1:dim(A)
      if iszero(a.coeffs[i])
        continue
      end
      for j = 1:dim(A)
        for k = 1:dim(A)
          M[j, k] += a.coeffs[i]*A.mult_table[i, j, k]
        end
      end
    end
  elseif action==:right
    for i = 1:dim(A)
      if iszero(a.coeffs[i])
        continue
      end
      for j = 1:dim(A)
        for k = 1:dim(A)
          M[j, k] += a.coeffs[i]*A.mult_table[j, i, k]
        end
      end
    end
  else
    error("Not yet implemented")
  end
  return M
end

################################################################################
#
#  isone/iszero
#
################################################################################

isone(a::AlgAssElem) = a == one(parent(a))

iszero(a::AlgAssElem) = all(i -> i == 0, a.coeffs)

ext: jl


text:
using EGraphs
using Documenter

makedocs(;
    modules=[EGraphs],
    authors="Philip Zucer",
    repo="https://github.com/philzook58/EGraphs.jl/blob/{commit}{path}#L{line}",
    sitename="EGraphs.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://philzook58.github.io/EGraphs.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/philzook58/EGraphs.jl",
)


text:
using PkgBenchmark
using BenchmarkHelper

import Whirl: Systems, TimeMarching, Fields
using Fields
using Systems
using TimeMarching
#using SaddlePointSystems

#include(joinpath(Pkg.dir("Whirl"), "src/Whirl.jl"))
#import Whirl
#import Whirl:@get, Soln, IntFactSystems, TimeMarching, Fields

#using IntFactSystems
#include(joinpath(Pkg.dir("Whirl"), "src/systems/navier_stokes.jl"))
#import TimeMarching: ifrk!, RK31

@benchgroup "Navier Stokes" begin
  woseen(x::Tuple{Real,Real},t;Re=1.0,x0::Tuple{Real,Real}=(0,0),t0=1) =
                      exp(-((x[1]-x0[1])^2+(x[2]-x0[2])^2)/(4(t+t0)/Re))/(1+t/t0)

  nx = 260; ny = 130
  Re = 200
  U = 1.0
  U∞ = (U,0.0)
  Lx = 2.0
  Δx = Lx/(ny-1)
  Δt = min(0.5*Δx,0.5*Δx^2*Re)
  w₀ = Nodes(Dual,(nx,ny))

  xg,yg = coordinates(w₀,dx=Δx)
  x0 = (1,1)
  t0 = 1
  wexact(t) = [woseen((x,y),t;Re=Re,x0=x0.+U∞.*t,t0=t0) for x in xg, y in yg]

  sys = Systems.NavierStokes((nx,ny),Re,Δx,Δt,U∞ = U∞)

  ifrk = IFRK(w₀,sys.Δt,
            (t,w) -> plan_intfact(t,w,sys),
            (w,t) -> r₁(w,t,sys) ,rk=TimeMarching.RK31)

  t = 0.0
  w₀ .= wexact(t)
  w = deepcopy(w₀)
  @bench "Single-step NS solver" $ifrk($t,$w)

#=
    wexact(x, y, t) = exp(-((x-t)^2+y^2)/(0.2^2+4t/200))/(1+4t/(200*0.2^2))

    Re = 200
    Δx = 0.02
    x = linspace(-1-0.5Δx, 3+0.5Δx, 202)
    y = linspace(-1-0.5Δx, 1+0.5Δx, 102)
    Δt = min(0.5*Δx,0.5*Δx^2*Re)
    dims = length.((x,y))

    ns = NavierStokes(dims, Re, Δx, Δt, (1.0, 0.0))
    s = Soln(0.0, Fields.DualNodes(dims));
    s.u .= [wexact(xᵢ, yᵢ, 0)*Δx for xᵢ in x, yᵢ in y]
    s₊ = Soln(0.0, Fields.DualNodes(dims));

    @bench "Single-Step NS Solver" ifrk!($s₊, $s, $Δt, $RK31, $ns)
    =#
end

#=
@benchgroup "Fields"  begin
    celldims = (500, 500)

    nodes = DualNodes(celldims)
    nodes .= rand(Float64,size(nodes))

    edges = Edges(Primal, nodes)
    edges.u .= reshape(1:length(edges.u), size(edges.u))
    edges.v .= reshape(1:length(edges.v), size(edges.v))

    dual = Edges(Dual, nodes)

    L = Laplacian(celldims, with_inverse = true, fftw_flags = FFTW.PATIENT)

    @bench "Shift Edges to Edges" Fields.shift!($dual, $edges)
    @bench "Shift Nodes to Edges" Fields.shift!($dual, $nodes)

    @bench "Node to Edge Curl" curl!($edges, $nodes)
    @bench "Edge to Node Curl" curl!($nodes, $edges)

    out = DualNodes(celldims)
    @bench "Laplacian" A_mul_B!($out, $L, $nodes)
    @bench "Inverse Laplacian" A_ldiv_B!($out, $L, $nodes)
end
=#


text:
# No regularization
# ------------------------------------------------------------
"""
    NoRegularization

Empty functor object for running an L-shaped algorithm without regularization.

"""
struct NoRegularization <: AbstractRegularization end

function initialize_regularization!(::AbstractLShapedSolver, ::NoRegularization)
    return nothing
end

function log_regularization!(::AbstractLShapedSolver, ::NoRegularization)
    return nothing
end

function log_regularization!(::AbstractLShapedSolver, ::Integer, ::NoRegularization)
    return nothing
end

function take_step!(::AbstractLShapedSolver, ::NoRegularization)
    return nothing
end

function decision(lshaped::AbstractLShapedSolver, ::NoRegularization)
    return lshaped.x
end

function objective(lshaped::AbstractLShapedSolver, ::NoRegularization)
    return lshaped.data.Q
end

function gap(lshaped::AbstractLShapedSolver, ::NoRegularization)
    @unpack Q,θ = lshaped.data
    return abs(θ-Q)/(abs(Q)+1e-10)
end

function add_regularization_params!(::NoRegularization; kwargs...)
    return nothing
end

# API
# ------------------------------------------------------------
"""
    DontRegularize

Factory object for [`NoRegularization`](@ref). Passed by default to `regularize ` in the `LShapedSolver` factory function.

"""
struct DontRegularize <: AbstractRegularizer end

function (::DontRegularize)(::AbstractVector)
    return NoRegularization()
end

function str(::DontRegularize)
    return "L-shaped"
end


text:
#!/usr/bin/env julia

using Downloads, JSON3, Base.BinaryPlatforms, Scratch, SHA, Pkg, TOML
using Base: UUID

# Download versions.json, start iterating over Julia versions
versions_json_url = "https://julialang-s3.julialang.org/bin/versions.json"
num_concurrent_downloads = 8

@info("Downloading versions.json...")
json_buff = IOBuffer()
Downloads.download(versions_json_url, json_buff)
versions_json = JSON3.read(String(take!(json_buff)))

# Collect all versions that are >= 1.0.0, and are a stable release
versions = filter(versions_json) do (v, d)
    if VersionNumber(string(v)) < v"1.0.0"
        return false
    end
    if !d["stable"]
        return false
    end
    return true
end

# Build download URLs for each one, and then tack on the next release as well
function select_url_hash(data, host = HostPlatform())
    d = first(filter(f -> platforms_match(f.triplet, host), data.files))
    return (d.url, d.sha256)
end
version_urls = sort(select_url_hash.(values(versions)), by = pair -> pair[1])

function generate_nightly_url(jlver, host = HostPlatform())
    # Map arch
    arch_str = Dict("x86_64" => "x64", "i686" => "x86", "aarch64" => "aarch64", "armv7l" => "armv7l", "ppc64le" => "ppc64le")[arch(host)]
    # Map OS name
    os_str = Dict("linux" => "linux", "windows" => "winnt", "macos" => "mac", "freebsd" => "freebsd")[os(host)]
    # Map wordsize tag
    wordsize_str = Dict("x86_64" => "64", "i686" => "32", "aarch64" => "aarch64", "armv7l" => "armv7l", "ppc64le" => "ppc64")[arch(host)]

    # If `jlver` is nothing, we don't namespace by version and just get the absolute latest version
    ver_str = ""
    if jlver !== nothing
        ver_str = string(jlver.major, ".", jlver.minor, "/")
    end

    return string(
        "https://julialangnightlies-s3.julialang.org/bin/",
        # linux/
        os_str, "/",
        # x64/
        arch_str, "/",
        # 1.6/ (or nothing, if `jlver === nothing`)
        ver_str,
        "julia-latest-",
        # linux64
        os_str, wordsize_str,
        ".tar.gz",
    )
end
highest_release = maximum(VersionNumber.(string.(keys(versions))))
next_release = VersionNumber(highest_release.major, highest_release.minor + 1, 0)
push!(version_urls, (generate_nightly_url(next_release), ""))
push!(version_urls, (generate_nightly_url(nothing),      ""))
@info("Identified $(length(version_urls)) versions to try...")

# Next, we're going to download each of these to a scratch space
scratch_dir = @get_scratch!("julia_installers")

# Ensure we always download the nightly
rm(joinpath(scratch_dir, basename(last(version_urls)[1])); force=true)

# Helper function to print out stdlibs from a Julia installation
function get_stdlibs(scratch_dir, julia_installer_name)
    installer_path = joinpath(scratch_dir, julia_installer_name)
    mktempdir() do dir
        @info("Extracting $(julia_installer_name)")
        mount_dir = joinpath(dir, "mount_dir")
        try
            if endswith(installer_path, ".dmg")
                mkdir(mount_dir)
                # Try to mount many times, as this seems to fail randomly
                mount_cmd = `hdiutil mount $(installer_path) -mountpoint $(mount_dir)`
                tries = 0
                while !success(mount_cmd)
                    if tries > 10
                        error("Unable to mount via hdiutil!")
                    end
                    sleep(0.1)
                    tries += 1
                end
                #@show readdir(mount_dir; join=true)
                app_dir = first(filter(d -> startswith(basename(d), "Julia-"), readdir(mount_dir; join=true)))
                symlink(joinpath(app_dir, "Contents", "Resources", "julia", "bin"), joinpath(dir, "bin"))
            elseif endswith(installer_path, ".exe")
                error("This script doesn't work with `.exe` downloads")
            else
                run(`tar -C $(dir) --strip-components=1 -zxf $(installer_path)`)
            end

            jlexe = joinpath(dir, "bin", @static Sys.iswindows() ? "julia.exe" : "julia")
            jlflags = ["--startup-file=no", "-O0"]
            jlvers = VersionNumber(readchomp(`$(jlexe) $(jlflags) -e 'print(VERSION)'`))
            jlvers = VersionNumber(jlvers.major, jlvers.minor, jlvers.patch)
            @info("Auto-detected Julia version $(jlvers)")

            if jlvers < v"1.1"
                stdlibs_str = readchomp(`$(jlexe) $(jlflags) -e 'import Pkg; print(repr(Pkg.Types.gather_stdlib_uuids()))'`)
            else
                stdlibs_str = readchomp(`$(jlexe) $(jlflags) -e 'import Pkg; print(repr(Pkg.Types.load_stdlib()))'`)
            end

            # This will give us a dictionary of UUID => (name, version) mappings for all standard libraries
            if jlvers < v"1.8"
                stdlibs = Dict{Base.UUID, Tuple}(uuid => (name, nothing) for (uuid, name) in eval(Meta.parse(stdlibs_str)))

                # We're going to try and get versions for each stdlib:
                stdlib_path = readchomp(`$(jlexe) $(jlflags) -e 'import Pkg; print(Pkg.Types.stdlib_path(""))'`)
                for uuid in keys(stdlibs)
                    # If this stdlib has a `Project.toml`, try to parse it for its version field
                    name = first(stdlibs[uuid])
                    project_path = joinpath(stdlib_path, name, "Project.toml")
                    if isfile(project_path)
                        d = TOML.parsefile(project_path)
                        if haskey(d, "version")
                            stdlibs[uuid] = (name, VersionNumber(d["version"]))
                        end
                    end
                end
            else
                stdlibs = Dict{Base.UUID, Tuple}(uuid => (name, version) for (uuid, (name, version)) in eval(Meta.parse(stdlibs_str)))
            end

            return (jlvers, stdlibs)
        finally
            # Clean up mounted directories
            if isdir(mount_dir)
                unmount_cmd = `hdiutil detach $(mount_dir)`
                tries = 0
                while !success(unmount_cmd)
                    if tries > 10
                        error("Unable to unmount $(mount_dir)")
                    end
                    sleep(0.1)
                    tries += 1
                end
            end
        end
    end
end

jobs = Channel()
versions_dict = Dict()

@sync begin
    # Feeder task
    @async begin
        for (url, hash) in version_urls
            put!(jobs, (url, hash))
        end
        close(jobs)
    end

    # Consumer tasks
    for _ in 1:num_concurrent_downloads
        @async begin
            for (url, hash) in jobs
                try
                    # We might try to download two files that have the same basename
                    url_tag = bytes2hex(sha256(url))
                    fname = joinpath(scratch_dir, string(url_tag, "-", basename(url)))
                    if !isfile(fname)
                        @info("Downloading $(url)")
                        Downloads.download(url, fname)
                    end

                    if !isempty(hash)
                        calc_hash = bytes2hex(open(io -> sha256(io), fname, "r"))
                        if calc_hash != hash
                            @error("Hash mismatch on $(fname); deleting and re-downloading")
                            rm(fname; force=true)
                            Downloads.download(url, fname)
                            calc_hash = bytes2hex(open(io -> sha256(io), fname, "r"))
                            if calc_hash != hash
                                @error("Hash mismatch on $(fname); re-download failed!")
                                continue
                            end
                        end
                    end

                    version, stdlibs = get_stdlibs(scratch_dir, basename(fname))
                    versions_dict[version] = stdlibs
                catch e
                    if isa(e, InterruptException)
                        rethrow()
                    end
                    @error(e, exception=(e, catch_backtrace()))
                end
            end
        end
    end
end

# Next, drop versions that are the same as the one "before" them:
sorted_versions = sort(collect(keys(versions_dict)))
versions_to_drop = VersionNumber[]
for idx in 2:length(sorted_versions)
    if versions_dict[sorted_versions[idx-1]] == versions_dict[sorted_versions[idx]]
        push!(versions_to_drop, sorted_versions[idx])
    end
end
for v in versions_to_drop
    delete!(versions_dict, v)
end

# Next, figure out which stdlibs are actually unresolvable, because they've never been registered
all_stdlibs = Dict{UUID,Tuple}()
for (julia_ver, stdlibs) in versions_dict
    merge!(all_stdlibs, stdlibs)
end

registries = Pkg.Registry.reachable_registries()
_unregistered_stdlibs_with_vers = filter(all_stdlibs) do (uuid, name)
    return !any(haskey(reg.pkgs, uuid) for reg in registries)
end

# For each entry in `UNREGISTERED_STDLIBS`, set the version number to `nothing`
unregistered_stdlibs = Dict((k, (v[1], nothing)) for (k, v) in pairs(_unregistered_stdlibs_with_vers))

# Helper function for getting these printed out in a nicely-sorted order
function print_sorted(io::IO, d::Dict; indent::Int=0)
    println(io, "Dict{UUID,Tuple{String,Union{VersionNumber,Nothing}}}(")
    for pair in sort(collect(d), by = kv-> kv[2][1])
        println(io, " "^indent, repr(pair[1]), " => ", repr(pair[2]), ",")
    end
    print(io, " "^(max(indent - 4, 0)), ")")
end

output_fname = joinpath(dirname(dirname(@__DIR__)), "src", "HistoricalStdlibs.jl")
@info("Outputting to $(output_fname)")
sorted_versions = sort(collect(keys(versions_dict)))
open(output_fname, "w") do io
    print(io, """
    ## This file autogenerated by ext/HistoricalStdlibGenerator/generate_historical_stdlibs.jl

    using Base: UUID

    # Julia standard libraries with duplicate entries removed so as to store only the
    # first release in a set of releases that all contain the same set of stdlibs.
    const STDLIBS_BY_VERSION = [
    """)
    for v in sorted_versions
        print(io, "    $(repr(v)) => ")
        print_sorted(io, versions_dict[v]; indent=8)
        println(io, ",")
    end
    println(io, "]")

    print(io, """
    # Next, we also embed a list of stdlibs that must _always_ be treated as stdlibs,
    # because they cannot be resolved in the registry; they have only ever existed within
    # the Julia stdlib source tree, and because of that, trying to resolve them will fail.
    const UNREGISTERED_STDLIBS = """)
    print_sorted(io, unregistered_stdlibs; indent=4)
    println(io)
end


text:
"""
    RPBVISolver <: Solver

POMDP solver type using the point-based value iteration algorithm.
"""
mutable struct RPBVISolver <: Solver
    beliefpoints::Vector{Vector{Float64}}
    max_iterations::Int64
    tolerance::Float64
end

"""
    RPBVISolver(; max_iterations, tolerance)

Initialize a point-based value iteration solver with the `max_iterations` limit and desired `tolerance`.
"""
function RPBVISolver(;beliefpoints::Vector{Vector{Float64}}=[[0, 1.0], [0.2, 0.8], [0.4, 0.6], [0.6,0.4], [0.8,0.2], [1.0,0.0]], max_iterations::Int64=10, tolerance::Float64=1e-3)
    return RPBVISolver(beliefpoints, max_iterations, tolerance)
end

"""
    AlphaVec

Alpha vector type of paired vector and action.
"""
struct AlphaVec
    alpha::Vector{Float64} # alpha vector
    action::Any # action associated wtih alpha vector
end

"""
    AlphaVec(vector, action_index)

Create alpha vector from `vector` and `action_index`.
"""
AlphaVec() = AlphaVec([0.0], 0)

# define alpha vector equality
==(a::AlphaVec, b::AlphaVec) = (a.alpha,a.action) == (b.alpha, b.action)
Base.hash(a::AlphaVec, h::UInt) = hash(a.alpha, hash(a.action, h))

"""
    create_policy(solver, pomdp)

Create AlphaVectorPolicy for `solver` using immediate rewards from `pomdp`.
"""
function create_policy(solver::RPBVISolver, pomdp::Union{POMDP,RPOMDP})
    ns = n_states(pomdp)
    na = n_actions(pomdp)
    S = ordered_states(pomdp)
    A = ordered_actions(pomdp)
    alphas = [[reward(pomdp,S[i],A[j]) for i in 1:ns] for j in 1:na]
    AlphaVectorPolicy(pomdp, alphas)
end

function create_policy(solver::RPBVISolver, pomdp::Union{IPOMDP,RIPOMDP})
    ns = n_states(pomdp)
    na = n_actions(pomdp)
    S = ordered_states(pomdp)
    A = ordered_actions(pomdp)
    b = fill(1 / ns, ns)
    alphas = [[reward(pomdp,b,A[j]) for i in 1:ns] for j in 1:na]
    AlphaVectorPolicy(pomdp, alphas)
end

# Eq. (5) Osogami 2015
# from p[z,a,t] to p[t,z,s]
function minutil(prob::Union{RPOMDP,RIPOMDP}, b::Vector{Float64}, a, alphavecs::Vector{Vector{Float64}})
    nz = n_observations(prob)
    ns = n_states(prob)
    nα = length(alphavecs)
    aind = RPOMDPModels.action_index(prob, a)
    plower, pupper = dynamics(prob)
    m = Model(solver = ClpSolver())
    @variable(m, u[1:nz])
    @variable(m, p[1:ns, 1:nz, 1:ns])
    @objective(m, Min, sum(u))
    for zind = 1:nz, αind = 1:nα
        @constraint(m, u[zind] >= sum(b[sind] * p[:,zind,sind]' * alphavecs[αind] for sind = 1:ns))
    end
    @constraint(m, p .<= pupper[:,:,:,aind])
    @constraint(m, p .>= plower[:,:,:,aind])
    for sind = 1:ns
        @constraint(m, sum(p[:,:,sind]) == 1)
    end
    JuMP.solve(m)
    getvalue(u), getvalue(p)
end

function findαz(zind::Int, u::Vector{Float64}, b::Vector{Float64}, p::Array{Float64}, alphavecs::Vector{Vector{Float64}})
    TOL = 1e-6
    αz = nothing
    ns = size(p, 1)
    for α in alphavecs
        s = sum(b[sind] * p[:,zind,sind]' * α for sind = 1:ns)
        # @show u[zind]
        # @show sum(b[sind] * p[:,zind,sind]' * α for sind = 1:ns)
        if isapprox(u[zind], sum(b[sind] * p[:,zind,sind]' * α for sind = 1:ns), atol = TOL)
            αz = α
            break
        end
    end
    # @show αz
    αz
end

function findαstar(rpomdp::Union{POMDP,RPOMDP}, b, s, a, pstar::Array{Float64}, alphaveczstar::Array{Array{Float64}})
    γ = rpomdp.discount
    ns = n_states(rpomdp)
    nz = n_observations(rpomdp)
    sind = state_index(rpomdp, s)
    αstars = reward(rpomdp, s, a)
    for tind in 1:ns, zind in 1:nz
        αstars += γ * pstar[tind, zind, sind] * alphaveczstar[zind][tind]
    end
    αstars
end

function findαstar(rpomdp::Union{IPOMDP,RIPOMDP}, b, s, a, pstar::Array{Float64}, alphaveczstar::Array{Array{Float64}})
    γ = rpomdp.discount
    ns = n_states(rpomdp)
    nz = n_observations(rpomdp)
    sind = state_index(rpomdp, s)
    αstars = rewardalpha(rpomdp, b, a)[sind]
    for tind in 1:ns, zind in 1:nz
        αstars += γ * pstar[tind, zind, sind] * alphaveczstar[zind][tind]
    end
    αstars
end

"""
    robustdpupdate(αset, beliefset, rpomdp)

Robust point-based dynamic programming backup value of `αset` for `beliefset` in `rpomdp`.
"""
function robustdpupdate(Vold::Vector{AlphaVec}, beliefset::Vector{Vector{Float64}}, rp::Union{RPOMDP,RIPOMDP})
    alphaset = [avec.alpha for avec in Vold]
    Vnew = Vector{AlphaVec}(size(beliefset, 1))
    @showprogress 1 "Updating..." for (bi, b) in enumerate(beliefset)
        Vbset = Set{AlphaVec}()
        for a in ordered_actions(rp)
            u, pstar = minutil(rp, b, a, alphaset)
            αz = Array{Array{Float64}}(n_observations(rp))
            αstar = Vector{Float64}(n_states(rp))
            for (zind, z) in enumerate(ordered_observations(rp))
                αz[zind] = findαz(zind, u, b, pstar, alphaset)
            end
            for (sind, s) in enumerate(ordered_states(rp))
                αstar[sind] = findαstar(rp, b, s, a, pstar, αz)
            end
            push!(Vbset, AlphaVec(αstar,a))
        end
        αmax = nothing
        vmax = -Inf
        for α in Vbset
            if b' * α.alpha > vmax
                vmax = b' * α.alpha
                αmax = α
            end
        end
        Vnew[bi] = αmax
    end
    Vnew
end

"""
    dpupdate(αset, beliefset, pomdp)

Point-based dynamic programming backup value of `αset` for `beliefset` in `pomdp`.
"""
function dpupdate(Vold::Vector{AlphaVec}, beliefset::Vector{Vector{Float64}}, prob::Union{POMDP,IPOMDP})
    alphaset = Set([avec.alpha for avec in Vold])
    Vnew = Vector{AlphaVec}(size(beliefset, 1))
    p = dynamics(prob)
    ns = n_states(prob)
    @showprogress 1 "Updating..." for (bi, b) in enumerate(beliefset)
        Vbset = Set{AlphaVec}()
        for (aind,a) in enumerate(ordered_actions(prob))
            αz = Array{Array{Float64}}(n_observations(prob))
            αstar = Vector{Float64}(n_states(prob))
            for (zind, z) in enumerate(ordered_observations(prob))
                smax = -Inf
                for αnew in collect(alphaset)
                    s = sum(b[sind] * p[:,zind,sind,aind]' * αnew for sind = 1:ns)
                    if s > smax
                        smax = s
                        αz[zind] = αnew
                    end
                end
            end
            for (sind, s) in enumerate(ordered_states(prob))
                αstar[sind] = findαstar(prob, b, s, a, p[:,:,:,aind], αz)
            end
            push!(Vbset, AlphaVec(αstar,a))
        end
        αmax = nothing
        vmax = -Inf
        for α in Vbset
            if b' * α.alpha > vmax
                vmax = b' * α.alpha
                αmax = α
            end
        end
        Vnew[bi] = αmax
    end
    Vnew
end

"""
    diffvalue(Vnew, Vold, pomdp)

Maximum difference between new alpha vectors `Vnew` and old alpha vectors `Vold` in `pomdp`.
"""
function diffvalue(Vnew::Vector{AlphaVec}, Vold::Vector{AlphaVec}, pomdp::Union{POMDP,IPOMDP,RPOMDP,RIPOMDP})
    nv = length(Vnew)
    ns = n_states(pomdp)
    dmax = -Inf
    for bi = 1:nv
        for si = 1:ns
            diff = abs(Vnew[bi].alpha[si] - Vold[bi].alpha[si])
            dmax = max(dmax, diff)
        end
    end
    dmax
end

"""
    solve(solver::RPBVISolver, rpomdp)

AlphaVectorPolicy for `pomdp` caluclated by the incremental pruning algorithm.
"""
function solve(solver::RPBVISolver, prob::Union{RPOMDP,RIPOMDP}; verbose::Bool=false)
    # println("Solver started...")
    ϵ = solver.tolerance
    iterlimit = solver.max_iterations
    policy = create_policy(solver, prob)
    avecs = [AlphaVec(policy.alphas[i], policy.action_map[i]) for i in 1:length(policy.action_map)]
    Vold = fill(AlphaVec(zeros(n_states(prob)), ordered_actions(prob)[1]), length(solver.beliefpoints))
    Vnew = Vector{AlphaVec}()
    del = Inf
    iter = 0
    while del > ϵ && iter < iterlimit
        iter += 1
        Vnew = robustdpupdate(Vold, solver.beliefpoints, prob)
        del = diffvalue(Vnew, Vold, prob)
        Vold = copy(Vnew)
        verbose && println("\riter: $iter, del: $del")
    end
    alphas_new = [v.alpha for v in Vnew]
    actions_new = [v.action for v in Vnew]
    policy = AlphaVectorPolicy(prob, alphas_new, actions_new)
    return policy
end

function solve(solver::RPBVISolver, prob::Union{POMDP,IPOMDP}; verbose::Bool=false)
    # println("Solver started...")
    ϵ = solver.tolerance
    iterlimit = solver.max_iterations
    policy = create_policy(solver, prob)
    avecs = [AlphaVec(policy.alphas[i], policy.action_map[i]) for i in 1:length(policy.action_map)]
    Vold = fill(AlphaVec(zeros(n_states(prob)), ordered_actions(prob)[1]), length(solver.beliefpoints))
    Vnew = Vector{AlphaVec}()
    del = Inf
    iter = 0
    while del > ϵ && iter < iterlimit
        iter += 1
        Vnew = dpupdate(Vold, solver.beliefpoints, prob)
        del = diffvalue(Vnew, Vold, prob)
        Vold = copy(Vnew)
        verbose && println("\riter: $iter, del: $del")
    end
    alphas_new = [v.alpha for v in Vnew]
    actions_new = [v.action for v in Vnew]
    policy = AlphaVectorPolicy(prob, alphas_new, actions_new)
    return policy
end

policyvalue(policy::AlphaVectorPolicy, b::Vector{Float64}) = maximum(dot(policy.alphas[i], b) for i in 1:length(policy.alphas))

max_stars_repo_name: ajkeith/RobustValueIteration

max_issues_repo_name: ajkeith/RobustValueIteration

max_forks_repo_name: ajkeith/RobustValueIteration


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

using Markdown
using InteractiveUtils

# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error).
macro bind(def, element)
    quote
        local el = $(esc(element))
        global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : missing
        el
    end
end

# ╔═╡ c823a222-6e1c-11eb-18a8-01a1b673d7b0
begin
	#import Pkg;
	#Pkg.add("LinearAlgebra")
	#Pkg.add("Random")
	#Pkg.add("PlutoUI")
	#Pkg.add("Plots")
	using LinearAlgebra
	using Random
	using Plots
    using PlutoUI
	md""" **0) import packages** """
end

# ╔═╡ 53b7c092-6e9b-11eb-0fa7-57fc1385457c
md"
# _The Gaussian Process_
Given a vector of pivot points $\boldsymbol x$,
the corresponding random vector  $\boldsymbol y(\boldsymbol x)$
is drawn from a multivariate Gaussian

$\boldsymbol z \sim {\cal G}(\mu,C)$

with zero mean $\boldsymbol \mu  = \boldsymbol0$
and a covariance with matrix elements: $C_{ij} = C(x_i,x_j)$.

The different Gaussian processes differ in the kernel function $C(x,y)$

1) Straight line process: $C(x,y) = \alpha \cdot x\cdot y$

2) Wiener process: $C(x,y) = \alpha * min(x, y)$

3) Squared exponential kernel: $C(x,y) = \alpha * \exp(- \frac{(x-y)^2}{\sigma^2})$

4) Ornstein-Uhlenbeck process: $C(x,y) = \alpha * \exp(- \frac{|x-y|}{\sigma})$

GPs are magic ✨
"

# ╔═╡ 0ac101ba-6ecc-11eb-090e-6b5998c6bbfd
md"
## Setting up the mathematical details
"

# ╔═╡ dca69d60-6e8c-11eb-3305-3fa2cedd5dac
begin
	Nslide = 100
	L_𝛼  = [0 100 ; 0.01 20 ; 0.01 40 ; 0.1 10]
	L_𝜎 = [-1 -1; -1 -1; 0.1 20; 1 20]
	L_d_𝛼 = 	(L_𝛼[:,2]-L_𝛼[:,1])/Nslide
	L_d_𝜎 = 	(L_𝜎[:,2]-L_𝜎[:,1])/Nslide
	md"""**2) range of parameters**"""
end

# ╔═╡ 1edc6c32-6e1d-11eb-2f8a-b9dd2aff7b99
begin
	L_x0 = [-25:0.5:25;]
	md""" **3) pivot points** """
end

# ╔═╡ 21cff596-6ecc-11eb-093c-f97fac9ab7cb
md"
## Try it out!
"

# ╔═╡ 2837bf6e-6eb2-11eb-2b2a-190fa3669085
md"""
**Choose a kernel:** $(@bind kernel Select(["1"=>"Straight lines","2"=>"Wiener process","3"=>"Squared exponential","4"=>"Ornstein-Uhlenbeck process"]))
"""

# ╔═╡ 0bef3462-6e1d-11eb-08f6-e3b2c36a0721

begin

if kernel == "1"
		txt = "Straight line";
		L_x    = L_x0;
		K(x, y, 𝛼, sigma) = 𝛼 * x * y;
	elseif kernel == "2"
		txt = "Wiener process";
		L_x    = [x for x in L_x0 if x >= 0];
		K(x, y, 𝛼, sigma) = 𝛼 * min(x, y);
	elseif kernel == "3"
		txt = "Squared exponential";
		L_x    = L_x0;
		K(x, y, 𝛼, 𝜎) = 𝛼 * exp(- (x - y)^2/𝜎^2);
	elseif kernel == "4"
		txt = "Ornstein-Uhlenbeck process";
		L_x    = L_x0;
		K(x, y, 𝜶, 𝜎) = 𝜶 * exp(abs(x - y)/𝜎);
	else
		error("$(kernel) not supported")
	end
	n = length(L_x);
	md"""**1) kernel definitions**"""
end

# ╔═╡ bc5fc60a-6e92-11eb-0003-b58331ceddcf
ik = parse(Int32,kernel);

# ╔═╡ 6be2bb7a-6eb6-11eb-2cfb-87383e68bf54
md"
**Specify the seed value of random number generator**

1 $(@bind seed Slider(1:1000)) 1000
"

# ╔═╡ 05d7032c-6e25-11eb-0b48-7f9a594b5b88
begin
	rng = MersenneTwister(seed)
	u = randn(rng,n, 1)
	md"""**4) normal random vector zero mean, unit variance**"""
end

# ╔═╡ 4fdde52e-6e8b-11eb-3199-8112d143a31a
md"""
The seed you chose: $(seed)
"""

# ╔═╡ 431df5da-6eb4-11eb-11af-898fdf3601d5
	if L_𝜎[ik,2] > 0
		md"""**Choose parameters 𝛼:**
	$(L_𝛼[ik,1]) $(@bind 𝛼 Slider(L_𝛼[ik,1]:L_d_𝛼[ik]:L_𝛼[ik,2])) $(L_𝛼[ik,2]) **and 𝜎:**
		$(L_𝜎[ik,1]) $(@bind 𝜎 Slider(L_𝜎[ik,1]:L_d_𝜎[ik]:L_𝜎[ik,2])) $(L_𝜎[ik,2])
		"""
	else
		𝜎 = L_𝜎[ik,2]
				md"""**Choose parameter 𝛼:**
	$(L_𝛼[ik,1]) $(@bind 𝛼 Slider(L_𝛼[ik,1]:L_d_𝛼[ik]:L_𝛼[ik,2])) $(L_𝛼[ik,2])
	"""
	end

# ╔═╡ 2506c99c-6e1d-11eb-31a2-dd1b99a80d8a
begin
	C = [K(x, y, 𝛼, 𝜎) for x in L_x, y in L_x];
	F = svd(C);
	A = F.U * Diagonal(sqrt.(F.S));
	z = A * u;
	md""" **5) sample from mvG** """
end

# ╔═╡ ac913758-6e27-11eb-1649-37b1be855f5b
if L_𝜎[ik,2] > 0
	md" 	𝛼 = $(𝛼),    𝜎 = $(𝜎)		"
else
	md" 	𝛼 = $(𝛼)"
end

# ╔═╡ b3e26166-6e1e-11eb-1708-25be9e2a339f
plot(
    L_x,
    z,
	line =(1,1,:line),
    grid = true,
    title = txt,
    palette = :tab10,
    legend = :none,
    label = false,
    marker = :cross,
    xlabel = "x",
    ylabel = "y",
    ylim = (-20, 20),
)

# ╔═╡ Cell order:
# ╟─53b7c092-6e9b-11eb-0fa7-57fc1385457c
# ╟─0ac101ba-6ecc-11eb-090e-6b5998c6bbfd
# ╟─c823a222-6e1c-11eb-18a8-01a1b673d7b0
# ╟─0bef3462-6e1d-11eb-08f6-e3b2c36a0721
# ╟─dca69d60-6e8c-11eb-3305-3fa2cedd5dac
# ╟─1edc6c32-6e1d-11eb-2f8a-b9dd2aff7b99
# ╟─05d7032c-6e25-11eb-0b48-7f9a594b5b88
# ╟─2506c99c-6e1d-11eb-31a2-dd1b99a80d8a
# ╟─21cff596-6ecc-11eb-093c-f97fac9ab7cb
# ╟─2837bf6e-6eb2-11eb-2b2a-190fa3669085
# ╟─bc5fc60a-6e92-11eb-0003-b58331ceddcf
# ╟─6be2bb7a-6eb6-11eb-2cfb-87383e68bf54
# ╟─4fdde52e-6e8b-11eb-3199-8112d143a31a
# ╟─431df5da-6eb4-11eb-11af-898fdf3601d5
# ╟─ac913758-6e27-11eb-1649-37b1be855f5b
# ╟─b3e26166-6e1e-11eb-1708-25be9e2a339f

max_stars_repo_name: fonsp/bayes-imoox

max_issues_repo_name: fonsp/bayes-imoox

max_forks_repo_name: fonsp/bayes-imoox


text:
function get_irnss_orbital_elements(obj::IRNSSPropagator)
    return jcall(obj, "getIRNSSOrbitalElements", IRNSSOrbitalElements, ())
end


text:
using FastClosures

function preallocated_LinearOperator(A)
  (n, m) = size(A)
  Ap = zeros(n)
  prod = @closure p -> mul!(Ap, A, p)
  Atq = zeros(m)
  tprod = @closure q -> mul!(Atq, transpose(A), q)
  T = eltype(A)
  F1 = typeof(prod)
  F2 = typeof(tprod)
  LinearOperator{T,F1,F2,F2}(n, m, false, false, prod, tprod, tprod)
end

function nonallocating_opEye(n)
  prod = @closure v -> v
  F = typeof(prod)
  LinearOperator{Float64,F,F,F}(n, n, true, true, prod, prod, prod)
end


text:
using Test

using StructJuMP

@testset "printhook" begin

    numScen = 2
    parent = StructuredModel(num_scenarios=numScen)

    @variable(parent, 0 <= x <= 1)
    @variable(parent, 0 <= y <= 1)

    @constraint(parent, x + y == 1)
    @objective(parent, Min, x*x + y)

    for i in 1:numScen
        bl = StructuredModel(parent=parent, id=i)
        @variable(bl, w >= 0)
        @constraint(bl, w - x - y <= 1)
        @objective(bl, Min, w*w + w)
    end

    str = string(parent)
    @test sprint(print, parent) == """
Min x² + y
Subject to
 x + y = 1.0
"""
    @test sprint(show, parent) == """
A JuMP Model
Minimization problem with:
Variables: 2
Objective function type: GenericQuadExpr{Float64,StructJuMP.StructuredVariableRef}
Constraint: 1
Names registered in the model: x, y"""
    @test sprint(show, "text/latex", parent) == """
\$\$ \\begin{alignat*}{1}\\min\\quad & x^2 + y\\\\
\\text{Subject to} \\quad & x + y = 1.0\\\\
\\end{alignat*}
 \$\$"""
    #@test occursin("Child", str)
end


text:
using JSON
using UUIDs
using BinaryBuilder: jll_uuid, build_project_dict

module TestJLL end

@testset "JLLs - utils" begin
    @test jll_uuid("Zlib_jll") == UUID("83775a58-1f1d-513f-b197-d71354ab007a")
    @test jll_uuid("FFMPEG_jll") == UUID("b22a6f82-2f65-5046-a5b2-351ab43fb4e5")

    project = build_project_dict("LibFoo", v"1.3.5", [Dependency("Zlib_jll"), Dependency(PackageSpec(name = "XZ_jll", version = v"2.4.6"))])
    @test project["deps"] == Dict("JLLWrappers" => "692b3bcd-3c85-4b1f-b108-f13ce0eb3210",
                                  "Artifacts"   => "56f22d72-fd6d-98f1-02f0-08ddc0907c33",
                                  "Pkg"         => "44cfe95a-1eb2-52ea-b672-e2afdf69b78f",
                                  "Zlib_jll"    => "83775a58-1f1d-513f-b197-d71354ab007a",
                                  "Libdl"       => "8f399da3-3557-5675-b5ff-fb832c97cbdb",
                                  "XZ_jll"      => "ffd25f8a-64ca-5728-b0f7-c24cf3aae800")
    @test project["name"] == "LibFoo_jll"
    @test project["uuid"] == "b250f842-3251-58d3-8ee4-9a24ab2bab3f"
    @test project["compat"] == Dict("julia" => "1.0", "XZ_jll" => "=2.4.6", "JLLWrappers" => "1.1.0")
    @test project["version"] == "1.3.5"
    # Make sure BuildDependency's don't find their way to the project
    @test_throws MethodError build_project_dict("LibFoo", v"1.3.5", [Dependency("Zlib_jll"), BuildDependency("Xorg_util_macros_jll")])
end

@testset "JLLs - building" begin
    mktempdir() do build_path
        name = "libfoo"
        version = v"1.0.0"
        sources = [DirectorySource(build_tests_dir)]
        # Build for this platform and a platform that isn't this one for sure:
        # FreeBSD.
        platforms = [platform, Platform("x86_64", "freebsd")]
        dependencies = [Dependency("Zlib_jll")]
        # The buffer where we'll write the JSON meta data
        buff = IOBuffer()

        # First: call `get_meta_json` twice, once for each platform, and write the
        # JSON meta data.  In this way we can test that merging multiple JSON
        # objects work correctly.
        for p in platforms
            dict = get_meta_json(
                name,
                version,
                sources,
                # Use a build script depending on the target platform.
                Sys.isfreebsd(p) ? libfoo_make_script : libfoo_meson_script,
                [p],
                # The products we expect to be build
                libfoo_products,
                dependencies;
            )
            # Generate the JSON file
            println(buff, JSON.json(dict))
        end

        # Now build for real
        autobuild(
            build_path,
            name,
            version,
            sources,
            # Use the particular build script we're interested in
            libfoo_make_script,
            platforms,
            # The products we expect to be build
            libfoo_products,
            dependencies;
            # Run audit passes to make sure the library has the correct soname
            skip_audit=false,
        )

        withenv(
            "JULIA_PKG_DEVDIR" => joinpath(build_path, "devdir"),
            # Let's pretend to be in Yggdrasil, set the relevant environment
            # variables.
            "YGGDRASIL" => "true",
            "BUILD_SOURCEVERSION" => "0123456789abcdef0123456789abcdef01234567",
            "PROJECT" => "L/$(name)",
        ) do
            # What follows loosely mimics what we do to build JLL packages in
            # Yggdrasil.
            buff = IOBuffer(strip(String(take!(buff))))
            objs = []
            while !eof(buff)
                push!(objs, BinaryBuilder.JSON.parse(buff))
            end

            # Merging modifies `obj`, so let's keep an unmerged version around
            objs_unmerged = deepcopy(objs)

            # Merge the multiple outputs into one
            merged = BinaryBuilder.merge_json_objects(objs)
            BinaryBuilder.cleanup_merged_object!(merged)
            BinaryBuilder.cleanup_merged_object!.(objs_unmerged)

            # Determine build version
            name = merged["name"]
            version = merged["version"]
            # Filter out build-time dependencies that will not go into the dependencies of
            # the JLL packages.
            dependencies = Dependency[dep for dep in merged["dependencies"] if !isa(dep, BuildDependency)]
            lazy_artifacts = merged["lazy_artifacts"]
            build_version = BinaryBuilder.get_next_wrapper_version(name, version)
            repo = "JuliaBinaryWrappers/$(name)_jll.jl"
            code_dir = joinpath(Pkg.devdir(), "$(name)_jll")
            download_dir = joinpath(build_path, "products")

            # Skip init of the remote repository
            # Filter out build-time dependencies also here
            for json_obj in [merged, objs_unmerged...]
                json_obj["dependencies"] = Dependency[dep for dep in json_obj["dependencies"] if !isa(dep, BuildDependency)]
            end

            tag = "$(name)-v$(build_version)"
            upload_prefix = "https://github.com/$(repo)/releases/download/$(tag)"

            # This loop over the unmerged objects necessary in the event that we have multiple packages being built by a single build_tarballs.jl
            for (i,json_obj) in enumerate(objs_unmerged)
                from_scratch = (i == 1)

                # A test to make sure merging objects and reading them back work
                # as expected.
                if json_obj["platforms"] == [Platform("x86_64", "freebsd")]
                    @test occursin("make install", json_obj["script"])
                else
                    @test occursin("MESON_TARGET_TOOLCHAIN", json_obj["script"])
                end

                BinaryBuilder.rebuild_jll_package(json_obj; download_dir=download_dir, upload_prefix=upload_prefix, verbose=false, lazy_artifacts=json_obj["lazy_artifacts"], from_scratch=from_scratch)
            end

            env_dir = joinpath(build_path, "foo")
            mkpath(env_dir)
            Pkg.activate(env_dir)
            Pkg.develop(PackageSpec(path=code_dir))
            @eval TestJLL using libfoo_jll
            @test 6.08 ≈ @eval TestJLL ccall((:foo, libfoo), Cdouble, (Cdouble, Cdouble), 2.3, 4.5)
            @test @eval TestJLL libfoo_jll.is_available()
        end
    end
end


text:
using Coverage

coverage = process_folder()

coverage = merge_coverage_counts(coverage, filter!(
    let prefixes = (joinpath(pwd(), "src", ""),
                    joinpath(pwd(), "deps", ""))
        c -> any(p -> startswith(c.filename, p), prefixes)
    end,
    LCOV.readfolder("test")))

Codecov.submit_local(process_folder())
Codecov.clean_folder(".")


text:
"""
Module that implements the Alpha Zero learning design.

The module exports (i) implementations of various two-player board games, like
tic-tac-toe and generalizations thereof, (ii) functions for creating neural
networks that act on the states of games, and (iii) functionality to generate
datasets via the selfplay of a Monte-Carlo Tree Search player assisted by
a neural model. Aspects (ii) and (iii) of Jtac are kept orthogonal to (i), such
that it is easy to extend the library with new games.
"""
module Jtac

# -------- Packages ---------------------------------------------------------- #

using Random, Statistics, LinearAlgebra, Printf, Distributed

import ProgressMeter, Crayons, BSON, NLopt

import AutoGrad, Knet
import Knet: identity,
             relu,
             elu,
             softmax,
             tanh,
             sigm

export AutoGrad, Knet

export identity,
       relu,
       elu,
       softmax,
       tanh,
       sigm

# -------- Utilities --------------------------------------------------------- #

include("util.jl")

export prepare,
       branch

# -------- Games ------------------------------------------------------------- #

include("game.jl")

export Game,
       Status,
       ActionIndex

export status,
       current_player,
       legal_actions,
       apply_action!,
       is_action_legal,
       representation,
       policy_length,
       random_playout,
       augment,
       draw,
       is_over,
       random_turn!

# -------- CPU/GPU Elements and NN Models ------------------------------------ #

include("feature.jl")
include("element.jl")
include("layer.jl")
include("model.jl")

export Feature,
       ConstantFeature

export Model

export apply,
       apply_features,
       to_gpu,
       to_cpu,
       swap,
       on_gpu,
       base_model,
       playing_model,
       training_model,
       gametype

export Pointwise,
       Dense,
       Conv,
       Deconv,
       Pool,
       Dropout,
       Batchnorm,
       Chain,
       Stack

export valid_insize,
       outsize,
       layers,
       @chain,
       @stack,
       @residual

# -------- Specific models implementations ----------------------------------- #

include("models/toy.jl")
include("models/neural.jl")
include("models/async.jl")

export DummyModel,
       RandomModel,
       RolloutModel,
       NeuralModel,
       Shallow,
       MLP,
       ShallowConv,
       Async

# -------- MCTS, Players, and Bayesian ELO rankings -------------------------- #

include("mc.jl")
include("player.jl")
include("belo.jl")

export RandomPlayer,
       MCTSPlayer,
       IntuitionPlayer,
       HumanPlayer

export pvp,
       name,
       think,
       decide,
       turn!,
       compete

export Ranking
export summary

# -------- Training ---------------------------------------------------------- #

include("dataset.jl")
include("loss.jl")
include("learning.jl")

export DataSet

export augment,
       minibatch,
       record_self,
       record_against,
       Loss,
       loss,
       caption,
       set_optimizer!,
       train_step!,
       train!,
       train_self!,
       train_against!,
       train_from_model!,
       with_contest

# -------- Distributed creation of datasets # -------------------------------- #

include("distributed.jl")

export record_self_distributed,
       record_against_distributed,
       compete_distributed

# -------- Saving and loading games and datasets ----------------------------- #

include("io.jl")

export save_model,
       load_model,
       save_dataset,
       load_dataset

# -------- Game implementations ---------------------------------------------- #

include("games/mnkgame.jl")
include("games/metatac.jl")
include("games/nim.jl")
include("games/morris.jl")

export TicTacToe,
       MNKGame,
       MetaTac,
       Nim,
       Morris

end # module Jtac

max_stars_repo_name: roSievers/jtac


text:
using StructuredProximalOperators
using LinearAlgebra
using Test

include("check_manifold.jl")

@testset "StructuredProximalOperators.jl" begin

    ## Testing manifolds
    # @testset "Manifold $M" for M in [
    #     l1Manifold{5}([0, 0, 1, 1, 1]),
    #     FixedRankMatrices{5,6,3,ℝ}(),
    # ]
    #     check_manifold(M)

    #     check_retraction(M)
    # end

    ## Testing regularizers
    include("regularizer_l1.jl")
    include("regularizer_lnuclear.jl")

end


text:
# Conjugate

export Conjugate

"""
**Convex conjugate**

    Conjugate(f)

Returns the convex conjugate (also known as Fenchel conjugate, or Fenchel-Legendre transform) of function `f`, that is
```math
f^*(x) = \\sup_y \\{ \\langle y, x \\rangle - f(y) \\}.
```
"""
struct Conjugate{T <: ProximableFunction} <: ProximableFunction
    f::T
    function Conjugate{T}(f::T) where {T<: ProximableFunction}
        if is_convex(f) == false
            error("`f` must be convex")
        end
        new(f)
    end
end

is_prox_accurate(f::Conjugate) = is_prox_accurate(f.f)
is_convex(f::Conjugate) = true
is_cone(f::Conjugate) = is_cone(f.f) && is_convex(f.f)
is_smooth(f::Conjugate) = is_strongly_convex(f.f)
is_strongly_convex(f::Conjugate) = is_smooth(f.f)
is_quadratic(f::Conjugate) = is_strongly_convex(f.f) && is_generalized_quadratic(f.f)
is_generalized_quadratic(f::Conjugate) = is_quadratic(f.f)
is_set(f::Conjugate) = is_convex(f.f) && is_support(f.f)
is_positively_homogeneous(f::Conjugate) = is_convex(f.f) && is_set(f.f)

fun_dom(f::Conjugate) = fun_dom(f.f)

Conjugate(f::T) where {T <: ProximableFunction} = Conjugate{T}(f)

# only prox! is provided here, call method would require being able to compute
# an element of the subdifferential of the conjugate

function prox!(y::AbstractArray{R}, g::Conjugate, x::AbstractArray{R}, gamma::R=R(1)) where R <: Real
    # Moreau identity
    v = prox!(y, g.f, x/gamma, 1/gamma)
    if is_set(g)
        v = R(0)
    else
        v = dot(x, y) - gamma * dot(y, y) - v
    end
    y .= x .- gamma .* y
    return v
end

# complex case, need to cast inner products to real

function prox!(y::AbstractArray{Complex{T}}, g::Conjugate, x::AbstractArray{Complex{T}}, gamma::R=R(1)) where {R <: Real, T <: RealOrComplex{R}}
    v = prox!(y, g.f, x/gamma, 1/gamma)
    if is_set(g)
        v = R(0)
    else
        v = real(dot(x, y)) - gamma * real(dot(y, y)) - v
    end
    y .= x .- gamma .* y
    return v
end

# naive implementation

function prox_naive(g::Conjugate, x::AbstractArray{T}, gamma=R(1)) where {R, T <: RealOrComplex{R}}
    y, v = prox_naive(g.f, x/gamma, 1/gamma)
    return x - gamma * y, if is_set(g) R(0) else real(dot(x, y)) - gamma * real(dot(y, y)) - v end
end

# TODO: hard-code conjugation rules? E.g. precompose/epicompose


text:
export TicTacToe

const BLANK = Int8(3)

struct TicTacToe
    player::Int
    board::Vector{Int8}
    state::Int
end

TicTacToe() = TicTacToe(1, fill(BLANK,9), 100)

Base.isequal(x::TicTacToe, y::TicTacToe) = isequal(x.board, y.board)
Base.hash(x::TicTacToe) = hash(x.board)

function transition(ttt::TicTacToe, action::Int)
    @assert ttt.board[action] == BLANK
    player = ttt.player == 1 ? 2 : 1
    board = copy(ttt.board)
    board[action] = ttt.player
    state = getstate(player, board)
    TicTacToe(player, board, state)
end

function getstate(player::Int, board::Vector{Int8})
    pats = [(1,2,3), (4,5,6), (7,8,9), (1,4,7), (2,5,8), (3,6,9), (1,5,9), (3,5,7)]
    for p in pats
        board[p[1]] == BLANK && continue
        if board[p[1]] == board[p[2]] == board[p[3]]
            @assert player != board[p[1]]
            return -1
        end
    end
    c = count(x -> x == BLANK, board)
    c > 0 ? 100 : 0
end

function legalactions(ttt::TicTacToe)
    ttt.state == 100 || return Int[]
    actions = Int[]
    for i = 1:length(ttt.board)
        ttt.board[i] == BLANK && push!(actions,i)
    end
    actions
end

function reward(ttt::TicTacToe)
    # ttt is final
    @assert ttt.state != 100
    Float32(ttt.state)
end

function Base.print(io::IO, ttt::TicTacToe)
    s = map(ttt.board) do i
        i == 1 && return "○"
        i == 2 && return "×"
        i == BLANK && return "_"
        throw("")
    end
    print(io, "player: $(ttt.player)\n$(s[1]) $(s[2]) $(s[3])\n$(s[4]) $(s[5]) $(s[6])\n$(s[7]) $(s[8]) $(s[9])")
end

mutable struct NN
    g
end

function NN()
    T = Float32
    embeds = zerograd(Uniform(-0.001,0.001)(T,10,3))
    h = lookup(Node(embeds), Node(name="x"))
    h = Linear(T,90,50)(h)
    h = relu(h)
    p = Linear(T,50,9)(h)
    v = Linear(T,50,1)(h)
    v = tanh(v)
    NN(Graph(p,v))
end

function (nn::NN)(ttt::TicTacToe)
    if ttt.player == 1
        x = ttt.board
    else
        board = copy(ttt.board)
        for i = 1:length(board)
            board[i] == BLANK && continue
            board[i] = board[i] == 1 ? 2 : 1
        end
        x = board
    end
    nn.g(Var(x))
end


text:
#simple MNIST example with training on gpu
# once we build model , we simply push
# the loss function, and datasets to gpu to achieve this.
using Base64
using Flux,Statistics,BenchmarkTools, MLBase, Flux.Data.MNIST
using Flux: onehotbatch, onecold
using Base.Iterators: repeated
using ColorTypes
using MLJBase
using CUDA

#load training dataset
imgs = Flux.Data.MNIST.images()
print(size(imgs))
# flatten all thsimple_mniste 28x28 into 764 input [reshape(image,:) for image in imgs ]
# Reorder the layout of the data for the ANN
imagestrip(image::Matrix{<:ColorTypes.Gray}) = Float32.(reshape(image, :))
X = hcat(imagestrip.(imgs)...)
#Target to one-hot
labels = Flux.Data.MNIST.labels()
Y = onehotbatch(labels,0:9)
# defining the model ( a neural network )
m = Chain(
                Dense(784,32,relu),
                Dense(32,10),
                softmax
         ) |> gpu

loss(x, y) = Flux.Losses.crossentropy(m(X),Y)
data = Flux.Data.DataLoader((X, Y), batchsize=128) |> gpu
#data = Base.Iterators.repeated(data,200)
opt = Flux.Optimise.ADAM(0.1) |> gpu

# test dataset
tX = hcat(float.(reshape.(Flux.Data.MNIST.images(:test), :))...) |> gpu
tY = onehotbatch(Flux.Data.MNIST.labels(:test), 0:9) |> gpu
# train for 10 epochs
using IterTools: ncycle
@btime Flux.Optimise.train!(loss, Flux.params(m), ncycle(data, 2), opt, cb =  () -> println(loss(tX,tY)))
acc(x, y) = mean(onecold((m(x))) .== onecold(y))
@show acc(tX, tY)

#@show MLJBase.confusion_matrix(onecold((m(tX))),onecold(tY))


text:
using KMarkovGenSeq
using Test

@testset "counting test" begin
    function test_func(s::String)
        bases = "AGCT"
        d = Dict{String,Int}()
        seq2idx = Dict{String,Int}()
        idx = 0
        db_freq2 = zeros(Float32, 64, 1)
        for b1 in bases, b2 in bases, b3 in bases
            qs = b1 * b2 * b3
            d[qs] = 0
            idx += 1
            seq2idx[qs] = idx
        end
        for b1 in bases, b2 in bases, b3 in bases
            qs = b1 * b2 * b3
            d[qs] += length(findall(qs, s; overlap=true))
        end
        for (k, v) in d
            db_freq2[seq2idx[k]] = v
        end
        for i in 1:4:64
            db_freq2[i:(i + 3), 1] ./= sum(@view(db_freq2[i:(i + 3), 1]))
        end
        db_freq2 .= .-log.(db_freq2)
        db_freq2[isnan.(db_freq2)] .= 10
        db_freq2[isinf.(db_freq2)] .= 10

        _, db_freq1 = KMarkovGenSeq.db_freq(["test_genome.fna"], 2)

        @test db_freq2 == db_freq1[:, :, 1]

        return rm("test_genome.fna")
    end
    @testset "normal sequence" begin
        open("test_genome.fna", "w") do io
            return write(io, """
                   >seq1| This is a dummy ref seq
                   TTTATTTATTCAGAATAATAAGTTATTCATTATTTTTATACTGAAGTTTTTATACTTAGGTATTAATATA
                   AGAGTAATTAGTTTTAGTCAGTTATATTGAGTCGCAGCGCAAGCTGCAAATTAAACTTTAAATTGAAGAG
                   TTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCTTAACACATGCAAGTCGAACGGTAACAGAGAATA
                   GCTTGCTATTTTGCTGACGAGTGGCGGACGGGTGAGTAATGCTTGGGAATTTGCCTTGTTGCGGGGGACA
                   """)
        end
        s = "TTTATTTATTCAGAATAATAAGTTATTCATTATTTTTATACTGAAGTTTTTATACTTAGGTATTAATATAAGAGTAATTAGTTTTAGTCAGTTATATTGAGTCGCAGCGCAAGCTGCAAATTAAACTTTAAATTGAAGAGTTTGATCATGGCTCAGATTGAACGCTGGCGGCAGGCTTAACACATGCAAGTCGAACGGTAACAGAGAATAGCTTGCTATTTTGCTGACGAGTGGCGGACGGGTGAGTAATGCTTGGGAATTTGCCTTGTTGCGGGGGACA"
        test_func(s)
    end

    @testset "sequence with unknown (with M)" begin
        open("test_genome.fna", "w") do io
            return write(io, """
                   >seq1| This is a dummy ref seq
                   CGCAAGCTGCAAATTAAACTTTAAATTGAAGAM
                   """)
        end
        s = "CGCAAGCTGCAAATTAAACTTTAAATTGAAGAM"
        test_func(s)
    end
end


text:
#=
Associated YT tutorial: https://youtu.be/TEWy9vZcxW4
=#

using Random

Random.seed!(0)

X = [ 1    2    3    2.5;
      2.0  5.0 -1.0  2.0;
     -1.5  2.7  3.3 -0.8]

# Julia doesn't have classes, so instead we make a struct.
mutable struct Layer_Dense
    weights::Array{Float64, 2}
    biases::Vector{Float64}
end

# Make a constructor for this struct.
function Layer_Dense(n_inputs::Int, n_neurons::Int)
    weights = 0.10 * randn(n_inputs, n_neurons)
    biases = zeros(n_neurons)

    return Layer_Dense(weights, biases)
end

# And define the forward function for this struct.
function forward(layer::Layer_Dense, inputs)
    return inputs * layer.weights .+ layer.biases'
end

# Now do the example.
layer1 = Layer_Dense(4,5)
layer2 = Layer_Dense(5,2)

output1 = forward(layer1, X)
output2 = forward(layer2, output1)

println(output2)


text:
"""
    function from_bson_file(::Type{T},s::AbstractString) where T
"""
function from_bson_file(::Type{T},s::AbstractString) where T
  dict = BSON.load(s)
  from_dict(T,dict)
end

"""
    to_bson_file(object,filename)
"""
function to_bson_file(object,filename)
  dict = to_dict(object)
  bson(filename, dict)
end


text:
function updateT0!(allvar::AllVar)
    @extract allvar : alg pbf
    @extract pbf : L
    @extract alg : verbose nprint epsconv maxiter mindec μext μint
    en = Inf
    assign = zeros(Int,L,2)
    ndec = 0
    initialize_msg0(allvar)

    for it in 1:maxiter
        onesweepT0!(allvar)
        ndec, assign = assignment!(allvar,ndec,assign)
        en = compute_en(allvar)
        if mod(it, nprint) == 0 && verbose
            @printf("it = %d -- en %.3f μext %.2f μint %.2f ndec %d\n", it, en,μext, μint,ndec)
        end
        if ndec >= mindec
            if verbose
                @printf("it = %d -- en %.3f μext %.2f μint %.2f ndec %d\n", it, en,μext, μint,ndec)
            end
            return it, :converged, en
        end
    end
    if verbose
        @printf("it = %d -- en %.3f μext %.2f μint %.2f ndec %d\n", maxiter, en,μext, μint,ndec)
    end
    return maxiter, :unconverged, en
end

function initialize_msg0(allvar::AllVar)
    @extract allvar : pbf
    @extract pbf : L N P B F newP newB newF

    for marg in (newP, newB, newF)
        marg[1,0] = -Inf
        marg[1,N+1] = -Inf
        value, idx = findmax(marg)
        marg .= marg .- marg[idx]
    end
    for i in 1:L
        for marg in (P[i], F[i], B[i])
            marg[1,0] = - Inf
            marg[1,N+1] = -Inf
            if i == 1
                for n = 1:N+1
                    marg[0,n] = -Inf
                end
            end
            if i == L
                for n = 0:N
                    marg[0,n] = -Inf
                end
            end
            value, idx = findmax(marg)
            marg .= marg .- marg[idx]
        end
    end

end

function onesweepT0!(allvar::AllVar)
    @extract allvar : pbf alg
    @extract pbf : L P B F newP newB newF
    @extract alg : damp

    v = randperm(L)
    @inbounds for idx in 1:L # i-th marginal 1 ≤ i ≤ L
        i = v[idx]
        newPBFiT0!(allvar, i::Int)
        p,f,b = P[i], F[i], B[i]
        p .= newP; f .= newF; b .= newB;
        #p .= damp .* p .+ (1.0 - damp) .* newP # inplace sum
        #f .= damp .* f .+ (1.0 - damp) .* newF
        #b .= damp .* b .+ (1.0 - damp) .* newB
    end
end

function newPBFiT0!(allvar::AllVar, i::Int)
    @extract allvar : pbf
    @extract pbf : L N forward central bckward newP newB newF

    for fcb in (forward,central,bckward)
        fill!(fcb, 0.0)
    end

    centralT0!(central, allvar, i)
    i > 1 &&  forwardT0!(forward, allvar, i)
    i < L && backwardT0!(bckward, allvar, i)

    if i==1 # case site = 1 first site
        newF .= central
        newP .= central .+ bckward
        newB .= newP
    elseif 1 < i < L #internal sites
        newF .= forward .+ central
        newP .= newF .+ bckward
        newB .= central .+ bckward
    elseif i == L# case site = L
        newF .= forward .+ central
        newB .= central
        newP .= newF
    else
        error("the end of the world has come ...")
    end
    newP[1,0] = -Inf
    newB[1,0] = -Inf
    newF[1,0] = -Inf
    newP[1,N+1] = -Inf
    newB[1,N+1] = -Inf
    newF[1,N+1] = -Inf
    if i == 1
        for n in 1:N+1
            newP[0,n] = -Inf
            newB[0,n] = -Inf
            newF[0,n] = -Inf
        end
    end
    if i == L
        for n in 0:N
            newP[0,n] = -Inf
            newB[0,n] = -Inf
            newF[0,n] = -Inf
        end
    end
    for marg in (newP,newB,newF)
        v, idx = findmax(marg)
        marg .= marg .- marg[idx]
    end
end

# function myargmax(Pj,nbeg,nend)
#     #argm = -Inf
#     val = (-1,-1,-Inf)
#     @inbounds for n = nbeg:nend
#         for x = 0:1
#             val = ifelse(Pj[x, n]>val[3],(x,n,Pj[x,n]),val)
#         end
#     end
#     @assert val[1] > -1
#     @assert val[2] > -1
#     return val[1],val[2]
# end

function myargmax(Pj,nbeg,nend)
    argm = -Inf
    xj = -1
    nj = -1
    @inbounds for n = nbeg:nend
        for x = 0:1
            if Pj[x, n] > argm
                argm = Pj[x, n]
                xj = x
                nj = n
            end
        end
    end
    @assert(xj > -1)
    @assert(nj > -1)
    return xj,nj
end

function centralT0!(central::Marginal, allvar::AllVar, i::Int)
    @extract allvar:pbf jh seq alg
    @extract pbf:N P L q
    @extract jh:J h
    @extract seq:intseq
    @extract alg:μext μint

    A0 = q
    # case xᵢ = 1 (match)
    @inbounds for ni = 1:N
        Ani = intseq[ni]
        scra = h[Ani, i]
        for j = 1:i-2
            xj, nj = myargmax(P[j], 0, ni - 1)
            Anj = (nj == 0 || xj == 0) ? A0 : intseq[nj]
            scra += J[Ani, Anj, i, j]
        end
        for j = i+2:L
            xj, nj = myargmax(P[j], ni + 1, N + 1)
            Anj = (xj == 0 || nj == N + 1) ? A0 : intseq[nj]
            scra += J[Ani, Anj, i, j]
        end
        central[1, ni] = scra
    end
    # case xᵢ = 0 (gap)
    @inbounds for ni = 0:N+1
        scra = 0
        Ani = A0
        for j = 1:i-2
            xj, nj = myargmax(P[j], 0, ni)
            Anj = (nj == 0 || nj == N + 1 || xj == 0) ? A0 : intseq[nj]
            scra += J[Ani, Anj, i, j]
        end
        for j = i+2:L
            xj, nj = myargmax(P[j], ni, N + 1)
            Anj = (nj == 0 || nj == N + 1 || xj == 0) ? A0 : intseq[nj]
            scra += J[Ani, Anj, i, j]
        end
        central[0, ni] = scra
    end
    @inbounds for ni = 0:N+1
        central[0, ni] =
            (ni == 0 || ni == N + 1) ? central[0, ni] - μext + h[A0, i] :
            central[0, ni] - μint + h[A0, i]
    end
    #cannot match n = 0 pointer
    #central[1,0] = -Inf
    #mcentral = maximum(central)
    #central = central .- mcentral
end

function forwardT0!(forward::Marginal,allvar::AllVar,i::Int)
    @extract allvar : pbf jh seq alg
    @extract pbf : N P F q
    @extract jh : J
    @extract seq : intseq
    @extract alg : λo λe
    fwdm1 = F[i-1]
    A0 = q
    # case xᵢ = 1 (match)
    @inbounds for ni in 1:N
        Ani = intseq[ni]
        scra = -Inf
        @simd for nim in 0:ni-1
           Anim = nim == 0 ? A0 : intseq[nim]
           Δn = ni - nim - 1
           pen = (Δn > 0)*(nim > 0)*(- λo[i] - λe[i]*(Δn - 1))
           scra = max(scra, fwdm1[0,nim] + J[Ani, A0,   i, i-1] + pen)
           scra = max(scra, fwdm1[1,nim] + J[Ani, Anim, i, i-1] + pen)
        end
        forward[1, ni] = scra
    end
    # case xᵢ = 0 gap
    @inbounds @simd for ni in 1:N
        Ani = intseq[ni]
        forward[0,ni] = max(fwdm1[1,ni]+J[A0,Ani,i,i-1],fwdm1[0,ni]+J[A0,A0,i,i-1])
    end
    forward[0,0] = fwdm1[0,0] + J[A0,A0,i,i-1]
    scra = -Inf
    @inbounds @simd for nim = 1:N
        Anim = intseq[nim]
        scra = max(scra, fwdm1[1,nim] + J[A0,Anim,i,i-1])
    end
    forward[0,N+1] = max(scra, fwdm1[0,N+1] + J[A0,A0,i,i-1])
end

function backwardT0!(backward::Marginal,allvar::AllVar,i::Int)
    @extract allvar : pbf jh seq alg
    @extract pbf : N P B q
    @extract jh : J
    @extract seq : intseq
    @extract alg : λo λe

    bckp1 = B[i+1]
    A0 = q

    # case xᵢ = 1 (match)
    @inbounds for ni in 1:N
        # case xᵢp = 0 (gap)
        Ani = intseq[ni]
        scra = bckp1[0,ni] + J[A0,Ani,i+1,i]
        scra = max(scra, bckp1[0,N+1] + J[A0,Ani,i+1,i])
        for nip in ni+1:N
            Anip = intseq[nip]
            Δn = nip - ni- 1
	        pen = (Δn > 0)*(- λo[i+1] - λe[i+1]*(Δn - 1))
            scra = max(scra, bckp1[1,nip] + J[Anip, Ani, i+1,i] + pen)
        end
        backward[1,ni] = scra
    end
    # case xᵢ = 0 (gap)
    @inbounds for ni in 0:N+1
        scra = bckp1[0,ni] + J[A0,A0,i+1,i]
        for nip in ni+1:N
            Anip = intseq[nip]
            Δn = nip - ni - 1
	        pen = (ni > 0)*(Δn > 0)*(-λo[i+1] - λe[i+1]*(Δn - 1))
            scra = max(scra, bckp1[1,nip], J[Anip,A0,i+1,i] + pen)
        end
        backward[0,ni] = scra
    end
    #mbackward = maximum(backward)
    #backward = backward .- mbackward
end

function assignment!(allvar::AllVar, ndec::Int, assign::Matrix{Int})

    @extract allvar : pbf
    @extract pbf : L N P
    count = 0
    @inbounds for i in 1:L
        Pi = P[i]
        #v, idx = findmax(Pi)
        new_a = argmax(Pi)
        #display(new_a)
        #count how many nodes change state
        if (new_a[1] != assign[i,1]) || (new_a[2] != assign[i,2])
            count += 1
        end
        assign[i,1] = new_a[1]
        assign[i,2] = new_a[2]
    end

    if count == 0
        ndec += 1
    else
        ndec = 0
    end

    return ndec, assign
end


text:
@testset "Evalcode" begin
    # see `converter/md_blocks:convert_code_block`
    # see `converter/lx/resolve_lx_input_*`
    # --------------------------------------------
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{scripts/test1}
        done.
        """ * J.EOS |> seval

    spath = joinpath(J.PATHS[:assets], "scripts", "test1.jl")
    @test isfile(spath)
    @test isapproxstr(read(spath, String), """
        # This file was generated by JuDoc, do not modify it. # hide
        a = 5\nprint(a^2)""")

    opath = joinpath(J.PATHS[:assets], "scripts", "output", "test1.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, raw"""

Simple code:

a = 5
            print(a^2)

            then:

25

            done.
""")
end

@testset "Eval (errs)" begin
    # see `converter/md_blocks:convert_code_block`
    # --------------------------------------------
    h = raw"""
        Simple code:
        ```python:scripts/testpy
        a = 5
        print(a**2)
        ```
        done.
        """ * J.EOS |> seval

    @test isapproxstr(h, raw"""

Simple code:

a = 5
            print(a**2)

            done.
""")
end

@testset "Eval (rinput)" begin
    h = raw"""
        Simple code:
        ```julia:/scripts/test2
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{/scripts/test2}
        done.
        """ * J.EOS |> seval

    spath = joinpath(J.PATHS[:folder], "scripts", "test2.jl")
    @test isfile(spath)
    @test occursin("a = 5\nprint(a^2)", read(spath, String))

    opath = joinpath(J.PATHS[:folder], "scripts", "output", "test2.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, """

Simple code:

a = 5
            print(a^2)

            then:

25

            done.
""")

    # ------------

    J.CUR_PATH[] = joinpath(J.PATHS[:src], "pages", "pg1.md")[lastindex(J.PATHS[:src])+2:end]

    h = raw"""
        Simple code:
        ```julia:./code/test2
        a = 5
        print(a^2)
        ```
        then:
        \input{output}{./code/test2}
        done.
        """ * J.EOS |> seval

    spath = joinpath(J.PATHS[:assets], "pages", "code", "test2.jl")
    @test isfile(spath)
    @test occursin("a = 5\nprint(a^2)", read(spath, String))

    opath = joinpath(J.PATHS[:assets], "pages", "code", "output" ,"test2.out")
    @test isfile(opath)
    @test read(opath, String) == "25"

    @test isapproxstr(h, """

Simple code:

a = 5
            print(a^2)

            then:

25
  done.
""")
end

@testset "Eval (module)" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        using LinearAlgebra
        a = [5, 2, 3, 4]
        print(dot(a, a))
        ```
        then:
        \input{output}{scripts/test1}
        done.
        """ * J.EOS |> seval
    # dot(a, a) == 54
    @test occursin("then:
54
 done.", h)
end

@testset "Eval (img)" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        write(joinpath(@__DIR__, "output", "test1.png"), "blah")
        ```
        then:
        \input{plot}{scripts/test1}
        done.
        """ * J.EOS |> seval
    @test occursin("then:  done.", h)
end

@testset "Eval (throw)" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test1
        sqrt(-1)
        ```
        then:
        \input{output}{scripts/test1}
        done.
        """ * J.EOS |> seval
    # errors silently
    @test occursin("then:
There was an error running the code: DomainError", h)
end

@testset "Eval (nojl)" begin
    h = raw"""
        Simple code:
        ```python:scripts/test1
        sqrt(-1)
        ```
        done.
        """ * J.EOS

    @test (@test_logs (:warn, "Eval of non-julia code blocks is not yet supported.") h |> seval) == "
Simple code:
sqrt(-1)\n
 done.
\n"
end

# temporary fix for 186: make error appear and also use `abspath` in internal include
@testset "eval #186" begin
    h = raw"""
        Simple code:
        ```julia:scripts/test186
        fn = "tempf.jl"
        write(fn, "a = 1+1")
        println("Is this a file? $(isfile(fn))")
        include(abspath(fn))
        println("Now: $a")
        rm(fn)
        ```
        done.
        \output{scripts/test186}
        """ * J.EOS |> seval
    @test isapproxstr(h, raw"""

Simple code:
fn = "tempf.jl"
            write(fn, "a = 1+1")
            println("Is this a file? $(isfile(fn))")
            include(abspath(fn))
            println("Now: $a")
            rm(fn)
 done.
Is this a file? true
            Now: 2

            """)
end


text:
# power divergence statistic for testing goodness of fit
# Cressie and Read 1984; Read and Cressie 1988
# lambda  =  1: Pearson's Chi-square statistic
# lambda ->  0: Converges to Likelihood Ratio test stat
# lambda -> -1: Converges to minimum discrimination information statistic (Gokhale, Kullback 1978)
# lambda =  -2: Neyman Modified chi-squared statistic (Neyman 1949)
# lambda = -.5: Freeman-Tukey statistic (Freeman, Tukey 1950)
# Mirror functionality of https://github.com/JuliaStats/HypothesisTests.jl/blob/master/src/power_divergence.jl

# Under regularity conditions, their asymptotic distributions are all the same (Drost 1989)
# Chi-squared null approximation works best for lambda near 2/3
function power_divergence_statistic(x::Matrix{T};
                                                lambda=1, min_freq=5) where T<:Integer

  nrows, ncols = size(x)
  n = sum(x)

  #validate
  any( x .< 0) || any( .!isfinite.(x)) ? error("all entries must be nonnegative and finite") : nothing
  n == 0 ? error("at least one entry must be positive") : nothing
  isfinite(nrows) && isfinite(ncols) && isfinite(nrows*ncols) ? nothing : error("Invalid number of rows or columns")

  if nrows > 1 && ncols > 1
    rowsums = mapslices(sum, x, dims=2)
    colsums = mapslices(sum, x, dims=1)
    df = (nrows-1)*(ncols-1)
    xhat = rowsums * colsums / n
    if( minimum(x) < min_freq )
        @warn "Min frequecy requirement violated - returning NaN:" minimum(x)
        return NaN
    end
   else
    error("Number of rows and columns must be greater than 1")
  end

  stat = 0
  if lambda == 0
    for i in 1:length(x)
      stat += x[i]*(log(x[i]) - log(xhat[i]))
    end
    stat *= 2
  elseif lambda == -1
    for i in 1:length(x)
      stat += xhat[i]*(log(xhat[i]) - log(x[i]))
    end
    stat *= 2
  else
    for i in 1:length(x)
      stat += x[i]*( (x[i]/(xhat[i]))^lambda - 1)
    end
    stat *= 2/(lambda*(lambda+1))
  end

  return stat
end

#Chi2 Statictic
function chi2_statistic(x::AbstractMatrix{T}; min_freq=5) where T<:Integer
  power_divergence_statistic(x,lambda=1.0, min_freq=min_freq)
end

function chi2_statistic(x::AbstractVector{T}, y::AbstractVector{T},
                        levels::UnitRange{T}; min_freq=5) where T<:Integer
  d = counts(x,y,levels)
  power_divergence_statistic(d,lambda=1.0, min_freq=min_freq)
end

#Likelihood ratio (G2) Statictic
function likelihood_ratio(x::AbstractMatrix{T}; min_freq=5) where T<:Integer
  power_divergence_statistic(x,lambda=0, min_freq=min_freq)
end

function likelihood_ratio(x::AbstractVector{T}, y::AbstractVector{T},
                                      levels::UnitRange{T}; min_freq=5) where T<:Integer
  d = counts(x,y,levels)
  power_divergence_statistic(d,lambda=0, min_freq=min_freq)
end

#Minimum discrimination information
function mod_likelihood_ratio(x::AbstractMatrix{T}; min_freq=5) where T<:Integer
  power_divergence_statistic(x,lambda=-1, min_freq=min_freq)
end

function mod_likelihood_ratio(x::AbstractVector{T}, y::AbstractVector{T},
                                          levels::UnitRange{T}; min_freq=5) where T<:Integer
  d = counts(x,y,levels)
  power_divergence_statistic(d,lambda=-1, min_freq=min_freq)
end

#Neyman's Statistic
function neyman_statistic(x::AbstractMatrix{T}; min_freq=5) where T<:Integer
  power_divergence_statistic(x,lambda=-2, min_freq=min_freq)
end

function neyman_statistic(x::AbstractVector{T}, y::AbstractVector{T},
                                          levels::UnitRange{T}; min_freq=5) where T<:Integer
  d = counts(x,y,levels)
  power_divergence_statistic(d,lambda=-2, min_freq=min_freq)
end

max_stars_repo_name: bcbi/BCBIStats

max_issues_repo_name: bcbi/BCBIStats

max_forks_repo_name: bcbi/BCBIStats


text:
using Documenter, AlternatingProjections

makedocs(sitename="AlternatingProjections.jl",
    modules = [AlternatingProjections],
    format = Documenter.HTML(
        # Use clean URLs, unless built as a "local" build
        prettyurls = true,
        canonical = "https://olejorik.github.io/AlternatingProjections.jl/stable/",
        # assets = ["assets/favicon.ico"],
        # highlights = ["yaml"],
    ),
    clean = false,
    authors = "O.S.",)

# using DocumenterLaTeX
# makedocs(
# #     format = LaTeX(),
#     format = LaTeX(platform = "docker"),
#     sitename="AlternatingProjections"
# )

deploydocs(repo = "github.com/olejorik/AlternatingProjections.jl.git",
    target = "build",)


text:
#!/bin/env julia

# The MIT License (MIT)
#
# Copyright (c) 2015 Alexander Samoilov
#
# Permission is hereby granted, free of charge, to any person obtaining a copy of
# this software and associated documentation files (the "Software"), to deal in
# the Software without restriction, including without limitation the rights to
# use, copy, modify, merge, publish, distribute, sublicense, and/or sell copies of
# the Software, and to permit persons to whom the Software is furnished to do so,
# subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in all
# copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
# FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE AUTHORS OR
# COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY, WHETHER
# IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
# CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.

using DocOpt
using Printf

macro p_str(s) s end
const rexp_dir = r"\s*([\w]+)/"
#const rexp_abc = r"^ABC:\s*([\w]+)\s+(\d+.\d+)\s+(\d+.\d+)\s+(\d+)"
const transcript_rexp_suf = p"\s*([\w]+)\s+(\d+.\d+)\s+(\d+.\d+)\s+(\d+)"
const rexp_abc_prefix = p"^ABC:"
const rexp_abcdef_prefix = p"^ABCDEF:"
const rexp_abc    = Regex("$rexp_abc_prefix$transcript_rexp_suf")
const rexp_abcdef = Regex("$rexp_abcdef_prefix$transcript_rexp_suf")

mutable struct TestSummary
    name::AbstractString
    exec_time::Float64
    grad_time::Float64
    ratio::Float64

    # a constructor
    TestSummary(n::AbstractString = "") = new(n, 0, 0, 0)
end

@enum State none test_name exec_time gradient_time

function main()

    script_name = basename(@__FILE__)
    doc = """$script_name

Postprocess transcript files - get timers statistics.

Usage:
  $script_name -h | --help
  $script_name [-v | --verbose] [-s | --def] <input> [<output>]

Options:
  -h --help                  Show this screen.
  -s --DEF                   DEF mode.
"""

    if length(ARGS) < 1
        #println("$rexp_abc")
        #println("$rexp_abcdef")
        println(doc)
        exit(0)
    end

    arguments = docopt(doc)

    verbose = arguments["--verbose"]
    def     = arguments["--def"]
    fname   = arguments["<input>"]

    rexp = def ? rexp_abcdef : rexp_abc
    state = none
    summary = TestSummary[]
    open(fname,"r") do f
        d = Dict{AbstractString, UInt}()
        test_summary = TestSummary()
        for line in eachline(f)
            # print(line) ##
            test_name_m = match(rexp_dir, line)
            if test_name_m != nothing
                if verbose println("Bingo test name: $(test_name_m[1])") end
                state = test_name
                test_summary = TestSummary(test_name_m[1])
            end
            abc_time_m = match(rexp, line)
            if abc_time_m != nothing
                if verbose
                    println("Bingo abc time counter: $(abc_time_m[1]) "
                          * "elapsed time: $(abc_time_m[2]) "
                          * "real_time: $(abc_time_m[3]) "
                          * "# of repetitions:  $(abc_time_m[4])")
                end
                if abc_time_m[1] == "EXEC"
                    state = exec_time
                    if test_summary.exec_time == 0.0
                        test_summary.exec_time = parse(Float64,abc_time_m[3])
                    end
                end
                if abc_time_m[1] == "GRADIENT"
                    state = gradient_time
                    test_summary.grad_time = parse(Float64,abc_time_m[3])
                    # get all we needed - computue ratio, push to summary and switch to `none`
                    test_summary.ratio = test_summary.grad_time / test_summary.exec_time
                    tname = test_summary.name
                    if !haskey(d,tname)
                        push!(summary,test_summary)
                        #println("$(test_summary)")
                        d[tname] = 1
                    else
                        d[tname] = d[tname] + 1
                    end
                    state = none
                end
            end
        end

        sort!(summary, lt = (a,b) -> a.ratio > b.ratio)
        for s in summary
            println("$s")
        end

        # now emit csv
        f = open("$(fname).csv", "w")
        @printf(f, "%-48.48s, %7.7s, %7.7s, %5.5s%%\n", "Test name", "Exec. s", "Grad. s", "ratio")
        for s in summary
            (name, etime, gtime, ratio) = (s.name, s.exec_time, s.grad_time, s.ratio)
            if ratio < 1.0 # cut erroneous
                @printf(f, "%-48.48s, %7.1f, %7.1f, %5.2f%%\n", name, etime, gtime, ratio * 100)
            end
        end
    end
end

main()

max_stars_repo_name: alsam/scriptics-samples

max_issues_repo_name: alsam/scriptics-samples

max_forks_repo_name: alsam/scriptics-samples


text:
using Htop
using Test

@testset "Htop.jl" begin
    help = readchomp(`$(htop) --help`)
    @test contains(help, r"^htop")
    version = readchomp(`$(htop) --version`)
    @test contains(version, r"^htop [.0-9]+$")
end


text:
import Flux: gpu
import Flux.Tracker: data, collect

SequenceLengths(x::AbstractArray{T, 1} where T) = [Int32(size(x[batch], 1)) for batch in 1:size(x, 1)]

function SequenceMask(x::AbstractArray{T, 3} where T, seqLens::AbstractArray{T, 1} where T)
    (S, D, N) = size(x)
    mask = zeros(Float32, S, D, N)
    for batch in 1:N
        len = seqLens[batch]
        @inbounds mask[:, :, batch] = vcat(ones(Float32, len, D), zeros(Float32, S - len, D))
    end
    # return data(gpu(mask))
    return mask |> Flux.gpu |> Flux.data
end

# _MakePad(rows::Integer, columns::Integer) = data(gpu(zeros(Float32, rows, columns)))
_MakePad(rows::Integer, columns::Integer) = zeros(Float32, rows, columns) |> Flux.gpu |> Flux.data
_PadMerge(x::AbstractArray, pad::AbstractArray) = collect(vcat(x, pad))

function SequencePad(x::AbstractArray{T, 1} where T, maxLen::Integer)
    seqLens = SequenceLengths(x)
    return cat([_PadMerge(x[batch], _MakePad(maxLen - seqLens[batch], size(x[batch], 2))) for batch in 1:size(x, 1)]..., dims=3) |> Flux.gpu
end

function ArrayPad(x::AbstractArray{T, 1} where T, maxLen::Integer)
    seqLens = SequenceLengths(x)
    return cat([vcat(x[batch], zeros(Float32, maxLen - seqLens[batch], size(x[batch], 2))) for batch in 1:size(x, 1)]..., dims=3) |> Flux.gpu
end


text:
using Plots

include("rw_1d.jl")

nsteps = length(ARGS)>0 ? parse(Int, ARGS[1]) : 1000
prob = length(ARGS)>1 ? parse(Float64, ARGS[1]) : 0.5

result = randomwalk(nsteps,prob)

pl = plot(result, xaxis=("t"), yaxis=("x"), color=[:black], legend=false)
savefig(pl, "rw_1d.pdf")

max_stars_repo_name: JuliaTagBot/JuliaBook-Samples

max_issues_repo_name: JuliaTagBot/JuliaBook-Samples

max_forks_repo_name: JuliaTagBot/JuliaBook-Samples


text:
using Documenter, DiffEqBayes

makedocs(
    sitename="DiffEqBayes.jl",
    authors="Chris Rackauckas et al.",
    clean=true,
    doctest=false,
    modules=[DiffEqBayes],

    format=Documenter.HTML(assets=["assets/favicon.ico"],
                           canonical="https://diffeqbayes.sciml.ai/stable/"),

    pages=[
        "DiffEqBayes.jl: Bayesian Parameter Estimation for Differential Equations" => "index.md",
        "Methods" => "methods.md",
        "Examples" => "examples.md"
    ]
)

deploydocs(
    repo="github.com/SciML/DiffEqBayes.jl.git";
    push_preview=true
)


text:
function plotPaths( )

actualDir=pwd()
cd(joinpath(actualDir,"output","images"))

bsample = DataFrame(t=vec(1:T+1),Sim1=vec(apath[:,1]),Sim2=vec(apath[:,2]))
bsample[:mtindex] = 1:size(bsample,1)     # Add an identifier
datamia=stack(bsample,[2:3])              # Reshape the data so we can use colors!

plot1=plot(datamia,x="t",y="value",color="variable",Geom.line,Guide.title("Time path of assets"))
draw(PNG("pathAssets.png", 24cm, 12cm), plot1)

# -------------------------------------------------------------------------

bsample = DataFrame(t=vec(1:T),C1=vec(cpath[:,1]),Y1=vec(ypath[:,1]))
bsample[:mtindex] = 1:size(bsample,1)     # Add an identifier
datamia=stack(bsample,[2:3])              # Reshape the data so we can use colors!

plot2=plot(datamia,x="t",y="value",color="variable",Geom.line,Guide.title("Time path of income and consumption Individual 1"))
draw(PNG("pathIndiv1.png", 24cm, 12cm), plot2)

# -------------------------------------------------------------------------

bsample = DataFrame(t=vec(1:T),C2=vec(cpath[:,2]),Y2=vec(ypath[:,2]))
bsample[:mtindex] = 1:size(bsample,1)     # Add an identifier
datamia=stack(bsample,[2:3])              # Reshape the data so we can use colors!

plot3=plot(datamia,x="t",y="value",color="variable",Geom.line,Guide.title("Time path of income and consumption Individual 2"))
draw(PNG("pathIndiv2.png", 24cm, 12cm), plot3)

cd(actualDir)

end

max_stars_repo_name: floswald/ucl-econ-julia

max_issues_repo_name: floswald/ucl-econ-julia

max_forks_repo_name: floswald/ucl-econ-julia


text:
using MLJBase
using OutlierDetection
using Statistics
using Test

import OutlierDetectionInterface
const OD = OutlierDetectionInterface

struct MinimalDetectorModel <: OD.DetectorModel end
struct MinimalUnsupervisedDetector <: OD.UnsupervisedDetector end
struct MinimalSupervisedDetector <: OD.SupervisedDetector end

score(X) = dropdims(mean(X, dims=1), dims=1)
OD.fit(::MinimalUnsupervisedDetector, X::OD.Data; verbosity)::OD.Fit = MinimalDetectorModel(), score(X)
OD.fit(::MinimalSupervisedDetector, X::OD.Data, y::OD.Labels; verbosity)::OD.Fit = MinimalDetectorModel(), score(X)
OD.transform(::Union{MinimalSupervisedDetector, MinimalUnsupervisedDetector},
             model::MinimalDetectorModel, X::OD.Data)::OD.Scores = score(X)

OD.@default_frontend(MinimalUnsupervisedDetector)
OD.@default_frontend(MinimalSupervisedDetector)

include("tests.jl")


text:
module HHChannels

export Stim, run_hh

using Sundials
using SIUnits.ShortUnits

immutable HHParam
    c_m       # membrane spacific capacitance
    gnabar    # Na channel cross-membrane conductivity
    gkbar     # K channel cross-membrane conductivity
    gl        # Leak conductivity
    ena       # Na channel reversal potential
    ek        # K channel reversal potential
    el        # Leak reversal potential
    q10       # temperature dependent rate coefficient
              # (= 3^((T-T₀)/10K) with T₀ = 6.3 °C)

    # constructor with default values, corresponding
    # to a resting potential of -65 mV and temperature 6.3 °C
    HHParam(;
        # default values from HH paper

        # For reversal potentials we use those computed using
        # the Nernst equation with the following values:
        #       R   8.3144598
        #       F   96485.33289
        #       nao 140   mM
        #       nai  10   mM
        #       ko    2.5 mM
        #       ki   64.4 nM
        # We don't use the default values for ena and ek taken
        # from the HH paper:
        #   ena    = 115.0mV + -65.0mV,
        #   ek     = -12.0mV + -65.0mV,
        ena    =  63.55148117386mV,
        ek     = -74.17164678272mV,

        c_m    = 0.01F*m^-2,
        gnabar = .12S*cm^-2,
        gkbar  = .036S*cm^-2,
        gl     = .0003S*cm^-2,
        el     = -54.3mV,
        q10    = 1
    ) = new(c_m, gnabar, gkbar, gl, ena, ek, el, q10)

end

immutable Stim
    t0        # start time of stimulus
    t1        # stop time of stimulus
    i_e       # stimulus current density

    Stim() = new(0s, 0s, 0A/m^2)
    Stim(t0, t1, i_e) = new(t0, t1, i_e)
end

vtrap(x,y) = x/(exp(x/y) - 1.0)

# "m" sodium activation system
function m_lims(v, q10)
    alpha = .1mV^-1 * vtrap(-(v+40mV),10mV)
    beta =  4 * exp(-(v+65mV)/18mV)
    sum = alpha + beta
    mtau = 1ms / (q10*sum)
    minf = alpha/sum
    return mtau, minf
end

# "h" sodium inactivation system
function h_lims(v, q10)
    alpha = 0.07*exp(-(v+65mV)/20mV)
    beta = 1 / (exp(-(v+35mV)/10mV) + 1)
    sum = alpha + beta
    htau = 1ms / (q10*sum)
    hinf = alpha/sum
    return htau, hinf
end

# "n" potassium activation system
function n_lims(v, q10)
    alpha = .01mV^-1 * vtrap(-(v+55mV),10mV)
    beta = .125*exp(-(v+65mV)/80mV)
    sum = alpha + beta
    ntau = 1ms / (q10*sum)
    ninf = alpha/sum
    return ntau, ninf
end

# Choose initial conditions for the system such that the gating variables
# are at steady state for the user-specified voltage v
function initial_conditions(v, q10)
    mtau, minf = m_lims(v, q10)
    htau, hinf = h_lims(v, q10)
    ntau, ninf = n_lims(v, q10)

    return (v, minf, hinf, ninf)
end

# Given time t and state (v, m, h, n),
# return (vdot, mdot, hdot, ndot)
function f(t, state; p=HHParam(), stim=Stim())
    v, m, h, n = state

    # calculate current density due to ion channels
    gna = p.gnabar*m*m*m*h
    gk = p.gkbar*n*n*n*n

    ina = gna*(v - p.ena)
    ik = gk*(v - p.ek)
    il = p.gl*(v - p.el)

    itot = ik + ina + il

    # calculate current density due to stimulus
    if t>=stim.t0 && t<stim.t1
        itot -= stim.i_e
    end

    # calculate the voltage dependent rates for the gating variables
    mtau, minf = m_lims(v, p.q10)
    htau, hinf = h_lims(v, p.q10)
    ntau, ninf = n_lims(v, p.q10)

    return (-itot/p.c_m, (minf-m)/mtau, (hinf-h)/htau, (ninf-n)/ntau)
end

function run_hh(t_end; v0=-65mV, stim=Stim(), param=HHParam(), sample_dt=0.01ms)
    v_scale = 1V
    t_scale = 1s

    v0, m0, h0, n0 = initial_conditions(v0, param.q10)
    y0 = [ v0/v_scale, m0, h0, n0 ]

    samples = collect(0s: sample_dt: t_end)

    fbis(t, y, ydot) = begin
        vdot, mdot, hdot, ndot =
            f(t*t_scale, (y[1]*v_scale, y[2], y[3], y[4]), stim=stim, p=param)

        ydot[1], ydot[2], ydot[3], ydot[4] =
            vdot*t_scale/v_scale, mdot*t_scale, hdot*t_scale, ndot*t_scale

        return Sundials.CV_SUCCESS
    end

    # Ideally would run with vector absolute tolerance to account for v_scale,
    # but this would prevent us using the nice cvode wrapper.

    res = Sundials.cvode(fbis, y0, map(t->t/t_scale, samples), abstol=1e-6, reltol=5e-10)

    # Use map here because of issues with type deduction with arrays and SIUnits.
    return samples, map(v->v*v_scale, res[:, 1])
end

end # module HHChannels

max_stars_repo_name: kabicm/arbor

max_issues_repo_name: kabicm/arbor

max_forks_repo_name: kabicm/arbor


text:
export Cylinder

const Cylinder = Prism{Sphere{2,4,Nothing}}

# Below, if we called Cylinder(c, ...) in the function body, it would call the inner
# constructor Prism{Sphere{2,4,Nothing}}(c, ...) because Cylinder = Prism{Sphere{2,4,Nothing}},
# which is not what we want.
# To call the outer constructor of Prism, we should call Prism(c, ...) instead of Cylinder(c, ...).
Cylinder(c::SVector{3,<:Real},
         r::Real,
         h::Real=Inf,
         a::SVector{3,<:Real}=SVector(0.0,0.0,1.0),
         data=nothing) where {D} =
    (â = normalize(a); Prism(c, Sphere(SVector(0.0,0.0),r), h, [orthoaxes(â)... â], data))

Cylinder(c::AbstractVector{<:Real},  # center of cylinder
         r::Real,  # radius of base
         h::Real=Inf,  # height of cylinder
         a::AbstractVector{<:Real}=[0.0,0.0,1.0],  # axis direction of cylinder
         data=nothing) =
    Cylinder(SVector{3}(c), r, h, SVector{3}(a), data)

# Return the bounds of the center cut with respect to the prism center.
function bounds_ctrcut(s::Cylinder)
    ax = inv(s.p)  # prism axes: columns are not only unit vectors, but also orthogonal
    r = s.b.r
    el = Ellipsoid(SVector(0.0,0.0,0.0), SVector(r,r,0.0), ax)  # center is set at origin to return bounds with respect to prism center

    return bounds(el)
end


text:
using Test
using NDTensors

@testset "NDTensors" begin
  @testset "$filename" for filename in
                           ["linearalgebra.jl", "dense.jl", "blocksparse.jl", "diag.jl"]
    println("Running $filename")
    include(filename)
  end
end

nothing


text:
module Types

export Process_info, Instance_info, Call_info, Requested_task, Completed_task

type Process_info
    userhost::AbstractString
    process_id::Integer
end

type Instance_info
    userhost::AbstractString
    password::AbstractString
    process_id::Integer
    gpu_id::Integer
    error_counter::Integer
    status::AbstractString # "OK" "BUSY" "ERROR"
end

type Requested_task
    id::Integer
    command::AbstractString
end

type Completed_task
    task::Requested_task
    instance::Instance_info
end

type Call_info
    active::Bool
    process_id::Integer
    task::Requested_task
    process::Base.Process
    log_file::IOStream
    log_pos::Integer
    log_start::DateTime
    log_last_update::DateTime
    instance::Instance_info
end

Instance_info(userhost, password, gpu_id) =
        Instance_info(userhost, password, 1, gpu_id, 0, "OK")
Instance_info(userhost, password) =
        Instance_info(userhost, password, 1, 0, 0, "OK")

end # module

max_stars_repo_name: Gollor/JuliaAutomation

max_issues_repo_name: Gollor/JuliaAutomation

max_forks_repo_name: Gollor/JuliaAutomation