text:
"""
    struct ErpenbeckThoss{T<:AbstractFloat} <: DiabaticModel

1D two-state diabatic system capable of modelling a molecule adsorbed on a metal surface
or a single-molecule junction.

In the two references, all of the parameters are identical except for the particle mass `m`
and the vertical shift `c` applied to the ϵ₀ state.
Both references modify the shift to ensure the quantum ground-state has an energy of 0 eV.
Note that the mass `m` is specified in atomic mass units (amu) **not** atomic units.
We calculate the offset automatically in the constructor from the Morse potential
zero-point energy.

# References

- PHYSICAL REVIEW B 97, 235452 (2018)
- J. Chem. Phys. 151, 191101 (2019)
"""
struct ErpenbeckThoss{T<:AbstractFloat} <: DiabaticModel
    Γ::T
    morse::AdiabaticModels.Morse{T}
    D₁::T
    D₂::T
    x₀′::T
    a′::T
    c::T
    V∞::T
    q::T
    ã::T
    x̃::T
    V̄ₖ::T
end

function ErpenbeckThoss(;
    Γ,

    m   = austrip(1.0u"u"),
    Dₑ  = austrip(3.52u"eV"),
    x₀  = austrip(1.78u"Å"),
    a   = austrip(1.7361u"Å^-1"),
    D₁  = austrip(4.52u"eV"),
    D₂  = austrip(0.79u"eV"),
    x₀′ = x₀,
    a′  = austrip(1.379u"Å^-1"),
    V∞  = austrip(-1.5u"eV"),

    q   = 0.05,
    ã   = austrip(0.5u"Å"),
    x̃   = austrip(3.5u"Å"),
    V̄ₖ  = sqrt(austrip(Γ)/2π)
)
    morse = AdiabaticModels.Morse(;Dₑ, x₀, a, m)
    c = -NQCModels.AdiabaticModels.eigenenergy(morse, 0) # Set c to offset zero-point energy

    return ErpenbeckThoss(austrip(Γ), morse, D₁, D₂, x₀′, a′, c, V∞, q, ã, x̃, V̄ₖ)
end

function NQCModels.potential(model::ErpenbeckThoss, r::Real)
    (;morse, c, D₁, D₂, x₀′, a′, V∞) = model
    ϵ₀(x) = NQCModels.potential(morse, x) + c
    ϵ₁(x) = D₁*exp(-2a′*(x-x₀′)) - D₂*exp(-a′*(x-x₀′)) + V∞

    (;q, ã, x̃, V̄ₖ) = model
    Vₖ(x) = V̄ₖ * ((1-q)/2*(1 - tanh((x-x̃)/ã)) + q)

    V11 = ϵ₀(r)
    V22 = ϵ₁(r)
    V12 = Vₖ(r)
    return Hermitian(SMatrix{2,2}(V11, V12, V12, V22))
end

function NQCModels.derivative(model::ErpenbeckThoss, r::Real)
    (;morse, D₁, D₂, x₀′, a′) = model
    ∂ϵ₀(x) = NQCModels.derivative(morse, x)
    ∂ϵ₁(x) = -2a′*D₁*exp(-2a′*(x-x₀′)) + a′*D₂*exp(-a′*(x-x₀′))

    (;q, ã, x̃, V̄ₖ) = model
    ∂Vₖ(x) = -V̄ₖ * (1-q)/2 * sech((x-x̃)/ã)^2 / ã

    D11 = ∂ϵ₀(r)
    D22 = ∂ϵ₁(r)
    D12 = ∂Vₖ(r)
    return Hermitian(SMatrix{2,2}(D11, D12, D12, D22))
end

NQCModels.nstates(::ErpenbeckThoss) = 2
NQCModels.ndofs(::ErpenbeckThoss) = 1


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "ORC"
version = v"4.31.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://gstreamer.freedesktop.org/src/orc/orc-0.4.31.tar.xz", "a0ab5f10a6a9ae7c3a6b4218246564c3bf00d657cbdf587e6d34ec3ef0616075")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/orc-0.4.31/
mkdir build && cd build
meson .. --cross-file=${MESON_TARGET_TOOLCHAIN}
ninja -j${nproc}
ninja install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()

# The products that we will ensure are always built
products = [
    LibraryProduct(["liborc-test-0.4", "liborc-test-0"], :liborc_test),
    ExecutableProduct("orcc", :orcc),
    LibraryProduct(["liborc-0.4", "liborc-0"], :liborc),
    ExecutableProduct("orc-bugreport", :orc_bugreport)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

ext: jl

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
################################################################################
#
# NfOrd/FacElem.jl : Factored elements over number fields
#
# This file is part of hecke.
#
# Copyright (c) 2015: Claus Fieker, Tommy Hofmann
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without
# modification, are permitted provided that the following conditions are met:
# * Redistributions of source code must retain the above copyright notice, this
#   list of conditions and the following disclaimer.
#
# * Redistributions in binary form must reproduce the above copyright notice,
#   this list of conditions and the following disclaimer in the documentation
#   and/or other materials provided with the distribution.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
# AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
# IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
# DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE
# FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
# DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
# SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
# CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
# OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
#
#
#  Copyright (C) 2015, 2016 Tommy Hofmann
#
################################################################################

# Get FacElem from ClassGrpCtx
function FacElem(x::ClassGrpCtx, y::fmpz_mat, j::Int)
  return FacElem(x.R, [ deepcopy(y[j, i]) for i in 1:ncols(y) ])
end

function FacElem(x::ClassGrpCtx, y::Array{fmpz, 1})
  return FacElem(x.R, [ deepcopy(y[i]) for i in 1:length(y) ])
end

# Get the trivial factored element from an ordinary element
function FacElem(x::nf_elem)
  z = FacElem{nf_elem, AnticNumberField}()
  z.fac[x] = fmpz(1)
  z.parent = FacElemMon(parent(x)::AnticNumberField)::FacElemMon{AnticNumberField}
  return z
end

function istorsion_unit(x::FacElem{T}, checkisunit::Bool = false, p::Int = 16) where T
  @vprint :UnitGroup 2 "Checking if factored element is torsion\n"

  if checkisunit
    _isunit(x) ? nothing : return false, p
  end

  K = base_ring(x)
  d = degree(K)
  r, s = signature(K)

  @vprint :UnitGroup 2 "Precision is now $(p) \n"
  l = 0
  @vprint :UnitGroup 2 "Computing conjugates ... \n"
  @v_do :UnitGroup 2 pushindent()
  A = ArbField(p, false)
  B = log(A(1) + A(1)//A(6) * log(A(d))//A(d^2))
  p = Int(upper_bound(-log(B)/log(A(2)), fmpz)) + 2

  cx = conjugates_arb_log(x, p) #this guarantees the result with an abs. error
                                # of 2^-p
  @v_do :UnitGroup 2 popindent()
  @vprint :UnitGroup 2 "Conjugates log are $cx\n"
  for i in 1:r
    k = abs(cx[i])
    if ispositive(k)
      return false, p
    elseif isnonnegative(B - k)
      l = l + 1
    else
      println("fail 1")
    end
  end
  for i in 1:s
    k = cx[r + i]//2
    if ispositive(k)
      return false, p
    elseif isnonnegative(B - k)
      l = l + 1
    else
      println("fail 2")
    end
  end

  if l == r + s
    return true, p
  end
  error("precision was not sufficient")
end

function factored_norm(x::FacElem{nf_elem, AnticNumberField})
  b = fmpq[]
  c = fmpz[]
  for (a, e) in x.fac
    n = norm(a)
    d = numerator(n)
    if !isone(d)
      push!(b, d)
      push!(c, e)
    end
    d = denominator(n)
    if !isone(d)
      push!(b, d)
      push!(c, -e)
    end
  end
  if length(b) == 0
    push!(b, fmpq(1))
    push!(c, 0)
  end
  f = FacElem(b, c)
  simplify!(f)
  return f
end

function norm(x::FacElem{nf_elem})
  return evaluate(factored_norm(x))
end

_base_ring(x::nf_elem) = parent(x)::AnticNumberField

_base_ring(x::FacElem{nf_elem}) = base_ring(x)::AnticNumberField

*(x::FacElem{nf_elem}, y::NfOrdElem) = x*elem_in_nf(y)

*(x::NfOrdElem, y::FacElem{nf_elem}) = y*x

function _conjugates_arb_log(A::FacElemMon{AnticNumberField}, a::nf_elem, abs_tol::Int)
  abs_tol = 1<<nbits(abs_tol)
  if haskey(A.conj_log_cache, abs_tol)
    if haskey(A.conj_log_cache[abs_tol], a)
      return A.conj_log_cache[abs_tol][a]::Array{arb, 1}
    else
      z = conjugates_arb_log(a, abs_tol)
      A.conj_log_cache[abs_tol][a] = z
      return z::Array{arb, 1}
    end
  else
    A.conj_log_cache[abs_tol] = Dict{nf_elem, arb}()
    z = conjugates_arb_log(a, abs_tol)
    A.conj_log_cache[abs_tol][a] = z
    return z::Array{arb, 1}
  end
end

function conjugates_arb(x::FacElem{nf_elem, AnticNumberField}, abs_tol::Int)
  d = degree(_base_ring(x))
  res = Array{acb}(undef, d)

  first = true

  for (a, e) in x.fac
    z = conjugates_arb(a, abs_tol)
    if first
      for j in 1:d
        res[j] = z[j]^e
      end
      first = false
    else
      for j in 1:d
        res[j] = res[j] * z[j]^e
      end
    end
    i = i + 1
  end

  return res
end

function conjugates_arb_log(x::FacElem{nf_elem, AnticNumberField}, abs_tol::Int)
  K = _base_ring(x)
  r1, r2 = signature(K)
  d = r1 + r2
  res = Array{arb}(undef, d)

  target_tol = abs_tol
  pr = abs_tol + nbits(maximum(abs, values(x.fac))) + nbits(length(x.fac))

  while true
    prec_too_low = false
    first = true
    for (a, e) in x.fac
      z = _conjugates_arb_log(parent(x), a, pr)
      if first
        for j in 1:d
          res[j] = parent(z[j])()::arb
          muleq!(res[j], z[j], e)
          if !radiuslttwopower(res[j], -target_tol) || !isfinite(res[j])
            prec_too_low = true
            break
          end
          #res[j] = x.fac[a] * z[j]
        end
        first = false
      else
        for j in 1:d
          addmul!(res[j], z[j], e)
          #res[j] = res[j] + x.fac[a]*z[j]
          if !radiuslttwopower(res[j], -target_tol) || !isfinite(res[j])
            prec_too_low = true
            break
          end
        end
      end
      if prec_too_low
        break
      end
    end
    if prec_too_low
      pr *= 2
      continue
    end

    for j in 1:d
      expand!(res[j], -target_tol)
    end
    return res
  end
end

function conjugates_arb_log(x::FacElem{nf_elem, AnticNumberField}, R::ArbField)
  z = conjugates_arb_log(x, -R.prec)
  return map(R, z)
end

@doc Markdown.doc"""
    valuation(a::FacElem{nf_elem, AnticNumberField}, P::NfOrdIdl) -> fmpz
The valuation of $a$ at $P$.
"""
function valuation(a::FacElem{nf_elem, AnticNumberField}, P::NfOrdIdl)
  val = fmpz(0)
  for (a, e) = a.fac
    val += valuation(a, P)*e
  end
  return val
end

@doc Markdown.doc"""
    valuation(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, p::NfOrdIdl)
    valuation(A::FacElem{NfOrdIdl, NfOrdIdlSet}, p::NfOrdIdl)
The valuation of $A$ at $P$.
"""
function valuation(A::FacElem{NfOrdIdl, NfOrdIdlSet}, p::NfOrdIdl)
  return sum(valuation(I, p)*v for (I, v) = A.fac)
end

function valuation(A::FacElem{NfOrdFracIdl, NfOrdFracIdlSet}, p::NfOrdIdl)
  return sum(valuation(I, p)*v for (I, v) = A.fac)
end

@doc Markdown.doc"""
     ideal(O::NfOrd, a::FacElem{nf_elem, AnticNumberField)
The factored fractional ideal $a*O$.
"""
function ideal(O::NfOrd, a::FacElem{nf_elem, AnticNumberField})
  de = Dict{NfOrdFracIdl, fmpz}()
  for (e, k) = a.fac
    I = ideal(O, e)
    add_to_key!(de, I, k)
    #if haskey(de, I)
    #  de[I] += k
    #else
    #  de[I] = k
    #end
  end
  return FacElem(de)
end

#the normalise bit ensures that the "log" vector lies in the same vector space
#well, the same hyper-plane, as the units
@doc Markdown.doc"""
    conjugates_arb_log_normalise(x::nf_elem, p::Int = 10)
    conjugates_arb_log_normalise(x::FacElem{nf_elem, AnticNumberField}, p::Int = 10)
The "normalised" logarithms, ie. the array $c_i\log |x^{(i)}| - 1/n\log|N(x)|$,
so the (weighted) sum adds up to zero.
"""
function conjugates_arb_log_normalise(x::nf_elem, p::Int = 10)
  K = parent(x)
  r,s = signature(K)
  c = conjugates_arb_log(x, p)
  n = sum(c)//degree(K)
  for i=1:r
    c[i] -= n
  end
  for i=r+1:r+s
    c[i] -= n
    c[i] -= n
  end
  return c
end

#the normalise bit ensures that the "log" vector lies in the same vector space
#well, the same hyper-plane, as the units
function conjugates_arb_log_normalise(x::FacElem{nf_elem, AnticNumberField}, p::Int = 10)
  K = base_ring(x)
  r,s = signature(K)
  c = conjugates_arb_log(x, p)
  n = sum(c)//degree(K)
  for i=1:r
    c[i] -= n
  end
  for i=r+1:r+s
    c[i] -= n
    c[i] -= n
  end
  return c
end

function _conj_arb_log_matrix_normalise_cutoff(u::Array{T, 1}, prec::Int = 32)::arb_mat where T
  z = conjugates_arb_log_normalise(u[1], prec)
  A = ArbMatSpace(parent(z[1]), length(u), length(z)-1, false)()
  for i=1:length(z)-1
    A[1,i] = z[i]
  end

  for j=2:length(u)
    z = conjugates_arb_log_normalise(u[j], prec)
    for i=1:length(z)-1
      A[j,i] = z[i]
    end
  end
  return A
end


text:
# Solves subsonic airfoil
#   ARGS[1] = path to input file

# create temporary directory to do the solve in
const dirname = "./work"
if isdir(dirname)
  error("$dirname already exists, please delete it before proceeding")
end

mkdir(dirname)
cd(dirname)

# load PDESolver and solve the PDE
using PDESolver
using EulerEquationMod

solvePDE(joinpath("../", ARGS[1]))


text: #TODO - add test


text:
using WackyOptions
using Documenter

makedocs(;
    modules=[WackyOptions],
    authors="tester",
    repo="https://x.com/tester/WackyOptions.jl/blob/{commit}{path}#L{line}",
    sitename="WackyOptions.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="http://example.com",
        assets=[
            "assets/static.txt",
        ],
    ),
    pages=[
        "Home" => "index.md",
    ],
    bar="baz",
    foo="bar",
)


text:
mutable struct SVD{P} <: Transformer
    predict::P
    iter::Int64
    function SVD(maxiter=1000)
        predict(A, r) = -svd(A, r, maxiter = iter)
        iter = maxiter
        new{typeof(predict)}(predict, iter)
    end
    function _svd(A, r; iter = 100)
        V = randn(size(A,2),r) # random initialization
        for _ = 1:maxiter
            W = A * V
            Z = A' * W
            V, R = mmids_gramschmidt(Z)
        end
        W = A * V
        S = [norm(W[:, i]) for i=1:size(W,2)] # singular values
        U = reduce(hcat,[W[:,i]/S[i] for i=1:size(W,2)]) # left singular vectors
        return U * S * V'
end
end


text:
# Transform escaped characters in JPaths back to their original value.
function unescape_jpath(raw::String)
    ret = replace(replace(raw, "~0" => "~"), "~1" => "/")
    m = match(r"%([0-9A-F]{2})", ret)
    if m !== nothing
        for c in m.captures
            ret = replace(ret, "%$(c)" => Char(parse(UInt8, "0x$(c)")))
        end
    end
    return ret
end

function type_to_dict(x)
    return Dict(name => getfield(x, name) for name in fieldnames(typeof(x)))
end

function update_id(uri::HTTP.URI, s::String)
    id2 = HTTP.URI(s)
    if !isempty(id2.scheme)
        return id2
    end
    els = type_to_dict(uri)
    delete!(els, :uri)
    els[:fragment] = id2.fragment
    if !isempty(id2.path)
        oldpath = match(r"^(.*/).*$", uri.path)
        els[:path] = oldpath == nothing ? id2.path : oldpath.captures[1] * id2.path
    end
    return HTTP.URI(; els...)
end

function get_element(schema, path::AbstractString)
    for element in split(path, "/"; keepempty = false)
        schema = _recurse_get_element(schema, unescape_jpath(String(element)))
    end
    return schema
end

function _recurse_get_element(schema::Any, ::String)
    error("unmanaged type in ref resolution $(typeof(schema)): $(schema).")
end

function _recurse_get_element(schema::AbstractDict, element::String)
    if !haskey(schema, element)
        error("missing property '$(element)' in $(schema).")
    end
    return schema[element]
end

function _recurse_get_element(schema::Vector, element::String)
    index = tryparse(Int, element)  # Remember that `index` is 0-indexed!
    if index === nothing
        error("expected integer array index instead of '$(element)'.")
    elseif index >= length(schema)
        error("item index $(index) is larger than array $(schema).")
    end
    return schema[index + 1]
end

function get_remote_schema(uri::HTTP.URI)
    r = HTTP.get(uri)
    if r.status != 200
        error("Unable to get remote schema at $uri. HTTP status = $(r.status)")
    end
    return Schema(JSON.parse(String(r.body)))
end

function find_ref(
    uri::HTTP.URI, id_map::AbstractDict, path::String, parent_dir::String
)
    if path == "" || path == "#"  # This path refers to the root schema.
        return id_map[string(uri)]
    elseif startswith(path, "#/")  # This path is a JPointer.
        return get_element(id_map[string(uri)], path[3:end])
    end
    uri = update_id(uri, path)
    els = type_to_dict(uri)
    delete!.(Ref(els), [:uri, :fragment])
    uri2 = HTTP.URI(; els...)
    is_file_uri = startswith(uri2.scheme, "file") || isempty(uri2.scheme)
    if is_file_uri && !isabspath(uri2.path)
        # Normalize a file path to an absolute path so creating a key is consistent.
        uri2 = HTTP.URIs.merge(uri2; path = abspath(joinpath(parent_dir, uri2.path)))
    end
    if !haskey(id_map, string(uri2))
        # id_map doesn't have this key so, fetch the ref and add it to id_map.
        id_map[string(uri2)] = if startswith(uri2.scheme, "http")
            @info("fetching remote ref $(uri2)")
            get_remote_schema(uri2).data
        else
            @assert is_file_uri
            @info("loading local ref $(uri2)")
            Schema(JSON.parsefile(uri2.path); parent_dir = dirname(uri2.path)).data
        end
    end
    return get_element(id_map[string(uri2)], uri.fragment)
end

# Recursively find all "$ref" fields and resolve their path.

resolve_refs!(::Any, ::HTTP.URI, ::AbstractDict, ::String) = nothing

function resolve_refs!(
    schema::Vector,
    uri::HTTP.URI,
    id_map::AbstractDict,
    parent_dir::String
)
    for s in schema
        resolve_refs!(s, uri, id_map, parent_dir)
    end
    return
end

function resolve_refs!(
    schema::AbstractDict,
    uri::HTTP.URI,
    id_map::AbstractDict,
    parent_dir::String
)
    if haskey(schema, "id") && schema["id"] isa String
        # This block is for draft 4.
        uri = update_id(uri, schema["id"])
    end
    if haskey(schema, "\$id") && schema["\$id"] isa String
        # This block is for draft 6+.
        uri = update_id(uri, schema["\$id"])
    end
    for (k, v) in schema
        if k == "\$ref" && v isa String
            # This ref has not been resolved yet (otherwise it would not be a String).
            # We will replace the path string with the schema element pointed at, thus
            # marking it as resolved. This should prevent infinite recursions caused by
            # self referencing.
            schema["\$ref"] = find_ref(uri, id_map, v, parent_dir)
        else
            resolve_refs!(v, uri, id_map, parent_dir)
        end
    end
    return
end

function build_id_map(schema::AbstractDict)
    id_map = Dict{String, Any}("" => schema)
    build_id_map!(id_map, schema, HTTP.URI())
    return id_map
end

build_id_map!(::AbstractDict, ::Any, ::HTTP.URI) = nothing

function build_id_map!(id_map::AbstractDict, schema::Vector, uri::HTTP.URI)
    build_id_map!.(Ref(id_map), schema, Ref(uri))
    return
end

function build_id_map!(id_map::AbstractDict, schema::AbstractDict, uri::HTTP.URI)
    if haskey(schema, "id") && schema["id"] isa String
        # This block is for draft 4.
        uri = update_id(uri, schema["id"])
        id_map[string(uri)] = schema
    end
    if haskey(schema, "\$id") && schema["\$id"] isa String
        # This block is for draft 6+.
        uri = update_id(uri, schema["\$id"])
        id_map[string(uri)] = schema
    end
    for value in values(schema)
        build_id_map!(id_map, value, uri)
    end
    return
end

"""

    Schema(schema::AbstractDict; parent_dir::String = abspath("."))

Create a schema but with `schema` being a parsed JSON created with `JSON.parse()` or
`JSON.parsefile()`.

`parent_dir` is the path with respect to which references to local schemas are resolved.

## Examples

    my_schema = Schema(JSON.parsefile(filename))
    my_schema = Schema(JSON.parsefile(filename); parent_dir = "~/schemas")
"""
struct Schema
    data::Union{AbstractDict, Bool}

    Schema(schema::Bool; kwargs...) = new(schema)

    function Schema(
        schema::AbstractDict;
        parent_dir::String = abspath("."),
        parentFileDirectory = nothing,
    )
        if parentFileDirectory !== nothing
            @warn("kwarg `parentFileDirectory` is deprecated. Use `parent_dir` instead.")
            parent_dir = parentFileDirectory
        end
        schema = deepcopy(schema)  # Ensure we don't modify the user's data!
        id_map = build_id_map(schema)
        resolve_refs!(schema, HTTP.URI(), id_map, parent_dir)
        return new(schema)
    end
end

"""
    Schema(schema::String; parent_dir::String = abspath("."))

Create a schema for document validation by parsing the string `schema`.

`parent_dir` is the path with respect to which references to local schemas are resolved.

## Examples

    my_schema = Schema(\"\"\"{
        \"properties\": {
            \"foo\": {},
            \"bar\": {}
        },
        \"required\": [\"foo\"]
    }\"\"\")

    # Assume there exists `~/schemas/local_file.json`:
    my_schema = Schema(
        \"\"\"{
            "\$ref": "local_file.json"
        }\"\"\",
        parent_dir = "~/schemas"
    )
"""
Schema(schema::String; kwargs...) = Schema(JSON.parse(schema); kwargs...)

Base.show(io::IO, ::Schema) = print(io, "A JSONSchema")


text:
"""
Default PowerSimulations Emulation Problem Type
"""
struct GenericEmulationProblem <: EmulationProblem end

"""
    EmulationModel(::Type{M},
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...) where {M<:EmulationProblem,
                      T<:PM.AbstractPowerFormulation}

This builds the optimization problem of type M with the specific system and template.

# Arguments

  - `::Type{M} where M<:EmulationProblem`: The abstract Emulation model type
  - `template::ProblemTemplate`: The model reference made up of transmission, devices,
    branches, and services.
  - `sys::PSY.System`: the system created using Power Systems
  - `jump_model::Union{Nothing, JuMP.Model}`: Enables passing a custom JuMP model. Use with care

# Output

  - `model::EmulationModel`: The Emulation model containing the model type, built JuMP model, Power
    Systems system.

# Example

template = ProblemTemplate(CopperPlatePowerModel, devices, branches, services)
OpModel = EmulationModel(MockEmulationProblem, template, system)

# Accepted Key Words

  - `optimizer`: The optimizer that will be used in the optimization model.
  - `warm_start::Bool`: True will use the current Emulation point in the system to initialize variable values. False initializes all variables to zero. Default is true
  - `system_to_file::Bool:`: True to create a copy of the system used in the model. Default true.
  - `export_pwl_vars::Bool`: True to export all the pwl intermediate variables. It can slow down significantly the solve time. Default is false.
  - `allow_fails::Bool`: True to allow the simulation to continue even if the optimization step fails. Use with care, default to false.
  - `optimizer_solve_log_print::Bool`: True to print the optimizer solve log. Default is false.
  - `direct_mode_optimizer::Bool` True to use the solver in direct mode. Creates a [JuMP.direct_model](https://jump.dev/JuMP.jl/dev/reference/models/#JuMP.direct_model). Default is false.
  - `initial_time::Dates.DateTime`: Initial Time for the model solve
  - `time_series_cache_size::Int`: Size in bytes to cache for each time array. Default is 1 MiB. Set to 0 to disable.
"""
mutable struct EmulationModel{M <: EmulationProblem} <: OperationModel
    name::Symbol
    template::ProblemTemplate
    sys::PSY.System
    internal::ModelInternal
    store::EmulationModelStore # might be extended to other stores for simulation
    ext::Dict{String, Any}

    function EmulationModel{M}(
        template::ProblemTemplate,
        sys::PSY.System,
        settings::Settings,
        jump_model::Union{Nothing, JuMP.Model}=nothing;
        name=nothing,
    ) where {M <: EmulationProblem}
        if name === nothing
            name = Symbol(typeof(template))
        elseif name isa String
            name = Symbol(name)
        end
        _, ts_count, _ = PSY.get_time_series_counts(sys)
        if ts_count < 1
            error(
                "The system does not contain Static TimeSeries data. An Emulation model can't be formulated.",
            )
        end
        internal = ModelInternal(
            OptimizationContainer(sys, settings, jump_model, PSY.SingleTimeSeries),
        )
        new{M}(name, template, sys, internal, EmulationModelStore(), Dict{String, Any}())
    end
end

function EmulationModel{M}(
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    name=nothing,
    optimizer=nothing,
    warm_start=true,
    system_to_file=true,
    initialize_model=true,
    initialization_file="",
    deserialize_initial_conditions=false,
    export_pwl_vars=false,
    allow_fails=false,
    calculate_conflict=false,
    optimizer_solve_log_print=false,
    detailed_optimizer_stats=false,
    direct_mode_optimizer=false,
    check_numerical_bounds=true,
    initial_time=UNSET_INI_TIME,
    time_series_cache_size::Int=IS.TIME_SERIES_CACHE_SIZE_BYTES,
) where {M <: EmulationProblem}
    settings = Settings(
        sys;
        initial_time=initial_time,
        optimizer=optimizer,
        time_series_cache_size=time_series_cache_size,
        warm_start=warm_start,
        system_to_file=system_to_file,
        initialize_model=initialize_model,
        initialization_file=initialization_file,
        deserialize_initial_conditions=deserialize_initial_conditions,
        export_pwl_vars=export_pwl_vars,
        allow_fails=allow_fails,
        calculate_conflict=calculate_conflict,
        optimizer_solve_log_print=optimizer_solve_log_print,
        detailed_optimizer_stats=detailed_optimizer_stats,
        direct_mode_optimizer=direct_mode_optimizer,
        check_numerical_bounds=true,
        horizon=1,
    )
    return EmulationModel{M}(template, sys, settings, jump_model; name=name)
end

"""
    EmulationModel(::Type{M},
    template::ProblemTemplate,
    sys::PSY.System,
    optimizer::MOI.OptimizerWithAttributes,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...) where {M <: EmulationProblem}

This builds the optimization problem of type M with the specific system and template

# Arguments

  - `::Type{M} where M<:EmulationProblem`: The abstract Emulation model type
  - `template::ProblemTemplate`: The model reference made up of transmission, devices,
    branches, and services.
  - `sys::PSY.System`: the system created using Power Systems
  - `jump_model::Union{Nothing, JuMP.Model}`: Enables passing a custom JuMP model. Use with care

# Output

  - `Stage::EmulationProblem`: The Emulation model containing the model type, unbuilt JuMP model, Power
    Systems system.

# Example

template = ProblemTemplate(CopperPlatePowerModel, devices, branches, services)
problem = EmulationModel(MyOpProblemType template, system, optimizer)

# Accepted Key Words

  - `initial_time::Dates.DateTime`: Initial Time for the model solve
  - `warm_start::Bool` True will use the current Emulation point in the system to initialize variable values. False initializes all variables to zero. Default is true
  - `export_pwl_vars::Bool` True will write the results of the piece-wise-linear intermediate variables. Slows down the simulation process significantly
  - `allow_fails::Bool` True will allow the simulation to continue if the optimizer can't find a solution. Use with care, can lead to unwanted behaviour or results
  - `optimizer_solve_log_print::Bool` Uses JuMP.unset_silent() to print the optimizer's log. By default all solvers are set to `MOI.Silent()`
  - `name`: name of model, string or symbol; defaults to the type of template converted to a symbol
"""
function EmulationModel(
    ::Type{M},
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...,
) where {M <: EmulationProblem}
    return EmulationModel{M}(template, sys, jump_model; kwargs...)
end

function EmulationModel(
    template::ProblemTemplate,
    sys::PSY.System,
    jump_model::Union{Nothing, JuMP.Model}=nothing;
    kwargs...,
)
    return EmulationModel{GenericEmulationProblem}(template, sys, jump_model; kwargs...)
end

"""
EmulationModel(directory::AbstractString)

Construct an EmulationProblem from a serialized file.

# Arguments

  - `directory::AbstractString`: Directory containing a serialized model.
  - `optimizer::MOI.OptimizerWithAttributes`: The optimizer does not get serialized.
    Callers should pass whatever they passed to the original problem.
  - `jump_model::Union{Nothing, JuMP.Model}` = nothing: The JuMP model does not get
    serialized. Callers should pass whatever they passed to the original problem.
  - `system::Union{Nothing, PSY.System}`: Optionally, the system used for the model.
    If nothing and sys_to_file was set to true when the model was created, the system will
    be deserialized from a file.
"""
function EmulationModel(
    directory::AbstractString,
    optimizer::MOI.OptimizerWithAttributes;
    jump_model::Union{Nothing, JuMP.Model}=nothing,
    system::Union{Nothing, PSY.System}=nothing,
    kwargs...,
)
    return deserialize_problem(
        EmulationModel,
        directory;
        jump_model=jump_model,
        optimizer=optimizer,
        system=system,
    )
end

function get_current_time(model::EmulationModel)
    execution_count = get_internal(model).execution_count
    initial_time = get_initial_time(model)
    resolution = get_resolution(model.internal.store_parameters)
    return initial_time + resolution * execution_count
end

function init_model_store_params!(model::EmulationModel)
    num_executions = get_executions(model)
    system = get_system(model)
    interval = resolution = PSY.get_time_series_resolution(system)
    base_power = PSY.get_base_power(system)
    sys_uuid = IS.get_uuid(system)
    model.internal.store_parameters = ModelStoreParams(
        num_executions,
        1,
        interval,
        resolution,
        base_power,
        sys_uuid,
        get_metadata(get_optimization_container(model)),
    )
    return
end

function build_pre_step!(model::EmulationModel)
    TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Build pre-step" begin
        if !is_empty(model)
            @info "EmulationProblem status not BuildStatus.EMPTY. Resetting"
            reset!(model)
        end
        # TODO-PJ: Temporary while are able to switch from PJ to POI
        container = get_optimization_container(model)
        container.built_for_recurrent_solves = true

        @info "Initializing Optimization Container For an EmulationModel"
        init_optimization_container!(
            get_optimization_container(model),
            get_network_formulation(get_template(model)),
            get_system(model),
        )
        @info "Initializing ModelStoreParams"
        init_model_store_params!(model)

        @info "Mapping Service Models"
        populate_aggregated_service_model!(get_template(model), get_system(model))
        populate_contributing_devices!(get_template(model), get_system(model))
        add_services_to_device_model!(get_template(model))

        handle_initial_conditions!(model)
        set_status!(model, BuildStatus.IN_PROGRESS)
    end
    return
end

"""
Implementation of build for any EmulationProblem
"""
function build!(
    model::EmulationModel{<:EmulationProblem};
    executions=1,
    output_dir::String,
    recorders=[],
    console_level=Logging.Error,
    file_level=Logging.Info,
    disable_timer_outputs=false,
)
    mkpath(output_dir)
    set_output_dir!(model, output_dir)
    set_console_level!(model, console_level)
    set_file_level!(model, file_level)
    TimerOutputs.reset_timer!(BUILD_PROBLEMS_TIMER)
    disable_timer_outputs && TimerOutputs.disable_timer!(BUILD_PROBLEMS_TIMER)
    file_mode = "w"
    add_recorders!(model, recorders)
    register_recorders!(model, file_mode)
    logger = configure_logging(model.internal, file_mode)
    try
        Logging.with_logger(logger) do
            try
                set_executions!(model, executions)
                TimerOutputs.@timeit BUILD_PROBLEMS_TIMER "Problem $(get_name(model))" begin
                    build_impl!(model)
                end
                set_status!(model, BuildStatus.BUILT)
                @info "\n$(BUILD_PROBLEMS_TIMER)\n"
            catch e
                set_status!(model, BuildStatus.FAILED)
                bt = catch_backtrace()
                @error "EmulationModel Build Failed" exception = e, bt
            end
        end
    finally
        unregister_recorders!(model)
        close(logger)
    end
    return get_status(model)
end

"""
Default implementation of build method for Emulation Problems for models conforming with  DecisionProblem specification. Overload this function to implement a custom build method
"""
function build_model!(model::EmulationModel{<:EmulationProblem})
    container = get_optimization_container(model)
    system = get_system(model)
    build_impl!(container, get_template(model), system)
    return
end

function reset!(model::EmulationModel{<:EmulationProblem})
    if built_for_recurrent_solves(model)
        set_execution_count!(model, 0)
    end
    model.internal.container = OptimizationContainer(
        get_system(model),
        get_settings(model),
        nothing,
        PSY.SingleTimeSeries,
    )
    model.internal.ic_model_container = nothing
    empty_time_series_cache!(model)
    empty!(get_store(model))
    set_status!(model, BuildStatus.EMPTY)
    return
end

function update_parameters!(model::EmulationModel, store::EmulationModelStore)
    update_parameters!(model, store.data_container)
    return
end

function update_parameters!(model::EmulationModel, data::DatasetContainer{DataFrameDataset})
    for key in keys(get_parameters(model))
        update_parameter_values!(model, key, data)
    end
    if !is_synchronized(model)
        update_objective_function!(get_optimization_container(model))
    end
    return
end

function update_initial_conditions!(
    model::EmulationModel,
    source::EmulationModelStore,
    ::InterProblemChronology,
)
    for key in keys(get_initial_conditions(model))
        update_initial_conditions!(model, key, source)
    end
    return
end

function update_model!(
    model::EmulationModel,
    source::EmulationModelStore,
    ini_cond_chronology,
)
    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Parameter Updates" begin
        update_parameters!(model, source)
    end
    TimerOutputs.@timeit RUN_SIMULATION_TIMER "Ini Cond Updates" begin
        update_initial_conditions!(model, source, ini_cond_chronology)
    end
    return
end

function update_model!(model::EmulationModel)
    update_model!(model, get_store(model), InterProblemChronology())
    return
end

function run_impl!(
    model::EmulationModel;
    optimizer=nothing,
    enable_progress_bar=progress_meter_enabled(),
    kwargs...,
)
    _pre_solve_model_checks(model, optimizer)
    internal = get_internal(model)
    # Temporary check. Needs better way to manage re-runs of the same model
    if internal.execution_count > 0
        error("Call build! again")
    end
    prog_bar = ProgressMeter.Progress(internal.executions; enabled=enable_progress_bar)
    initial_time = get_initial_time(model)
    for execution in 1:(internal.executions)
        TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Run execution" begin
            solve_impl!(model)
            current_time = initial_time + (execution - 1) * PSI.get_resolution(model)
            write_results!(get_store(model), model, execution, current_time)
            write_optimizer_stats!(get_store(model), get_optimizer_stats(model), execution)
            advance_execution_count!(model)
            update_model!(model)
            ProgressMeter.update!(
                prog_bar,
                get_execution_count(model);
                showvalues=[(:Execution, execution)],
            )
        end
    end
    return
end

"""
Default run method for problems that conform to the requirements of
EmulationModel{<: EmulationProblem}

This will call [`build!`](@ref) on the model if it is not already built. It will forward all
keyword arguments to that function.

# Arguments

  - `model::EmulationModel = model`: Emulation model
  - `optimizer::MOI.OptimizerWithAttributes`: The optimizer that is used to solve the model
  - `executions::Int`: Number of executions for the emulator run
  - `export_problem_results::Bool`: If true, export ProblemResults DataFrames to CSV files.
  - `output_dir::String`: Required if the model is not already built, otherwise ignored
  - `enable_progress_bar::Bool`: Enables/Disable progress bar printing
  - `serialize::Bool`: If true, serialize the model to a file to allow re-execution later.

# Examples

status = run!(model; optimizer = GLPK.Optimizer, executions = 10)
status = run!(model; output_dir = ./model_output, optimizer = GLPK.Optimizer, executions = 10)
"""
function run!(
    model::EmulationModel{<:EmulationProblem};
    export_problem_results=false,
    console_level=Logging.Error,
    file_level=Logging.Info,
    disable_timer_outputs=false,
    serialize=true,
    kwargs...,
)
    build_if_not_already_built!(
        model;
        console_level=console_level,
        file_level=file_level,
        disable_timer_outputs=disable_timer_outputs,
        kwargs...,
    )
    set_console_level!(model, console_level)
    set_file_level!(model, file_level)
    TimerOutputs.reset_timer!(RUN_OPERATION_MODEL_TIMER)
    disable_timer_outputs && TimerOutputs.disable_timer!(RUN_OPERATION_MODEL_TIMER)
    file_mode = "a"
    register_recorders!(model, file_mode)
    logger = configure_logging(model.internal, file_mode)
    try
        Logging.with_logger(logger) do
            try
                initialize_storage!(
                    get_store(model),
                    get_optimization_container(model),
                    model.internal.store_parameters,
                )
                TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Run" begin
                    run_impl!(model; kwargs...)
                    set_run_status!(model, RunStatus.SUCCESSFUL)
                end
                if serialize
                    TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Serialize" begin
                        optimizer = get(kwargs, :optimizer, nothing)
                        serialize_problem(model, optimizer=optimizer)
                        serialize_optimization_model(model)
                    end
                end
                TimerOutputs.@timeit RUN_OPERATION_MODEL_TIMER "Results processing" begin
                    results = ProblemResults(model)
                    serialize_results(results, get_output_dir(model))
                    export_problem_results && export_results(results)
                end
                @info "\n$(RUN_OPERATION_MODEL_TIMER)\n"
            catch e
                @error "Emulation Problem Run failed" exception = (e, catch_backtrace())
                set_run_status!(model, RunStatus.FAILED)
            end
        end
    finally
        unregister_recorders!(model)
        close(logger)
    end
    return get_run_status(model)
end

"""
Default solve method for an EmulationModel used inside of a Simulation. Solves problems that conform to the requirements of DecisionModel{<: DecisionProblem}

# Arguments

  - `step::Int`: Simulation Step
  - `model::OperationModel`: operation model
  - `start_time::Dates.DateTime`: Initial Time of the simulation step in Simulation time.
  - `store::SimulationStore`: Simulation output store

# Accepted Key Words

  - `exports`: realtime export of output. Use wisely, it can have negative impacts in the simulation times
"""
function solve!(
    step::Int,
    model::EmulationModel{<:EmulationProblem},
    start_time::Dates.DateTime,
    store::SimulationStore;
    exports=nothing,
)
    # Note, we don't call solve!(decision_model) here because the solve call includes a lot of
    # other logic used when solving the models separate from a simulation
    solve_impl!(model)
    @assert get_current_time(model) == start_time
    if get_run_status(model) == RunStatus.SUCCESSFUL
        advance_execution_count!(model)
        write_results!(
            store,
            model,
            get_execution_count(model),
            start_time;
            exports=exports,
        )
        write_optimizer_stats!(store, model, get_execution_count(model))
    end
    return get_run_status(model)
end


text:
# Tests for HINT layer (Kruse et al, 2020)
# Author: Philipp Witte, pwitte3@gatech.edu
# Date: January 2020

using InvertibleNetworks, LinearAlgebra, Test, Random

Random.seed!(11)

#######################################################################################################################
# Test invertibility
# Input
nx = 16
ny = 16
n_channel = 4
n_hidden = 8
batchsize = 2

function test_inv(nx, ny, n_channel, n_hidden, batchsize, permute, logdet, rev)
    # Input image
    X = randn(Float32, nx, ny, n_channel, batchsize)

    # HINT layer w/o logdet
    HL = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet)
    rev && (HL = reverse(HL))
    # Test
    Y = logdet ? HL.forward(X)[1] : HL.forward(X)
    X_ = HL.inverse(Y)
    @test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-5)

    X_ = HL.backward(0f0.*Y, Y)[2]
    @test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-5)
end

function loss(HL, X)
    Y = HL.forward(X)
    f = .5*norm(Y)^2
    ΔY = copy(Y)
    ΔX, X_ = HL.backward(ΔY, Y)
    return f, ΔX, HL.CL[1].RB.W1.grad, X_
end

function loss_logdet(HL, X)
    Y, logdet = HL.forward(X)
    f = -log_likelihood(Y) - logdet
    ΔY = -∇log_likelihood(Y)
    ΔX, X_ = HL.backward(ΔY, Y)
    return f, ΔX, HL.CL[1].RB.W1.grad, X_
end

function grad_test_X(nx, ny, n_channel, n_hidden, batchsize, permute, logdet, rev)
    logdet ? lossf = loss_logdet : lossf = loss
    # Input image
    X0 = rand(Float32, nx, ny, n_channel, batchsize)
    dX = randn(Float32, nx, ny, n_channel, batchsize)
    # Test for input X
    HL = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet)
    rev && (HL = reverse(HL))

    f0, gX, X_ = lossf(HL, X0)[[1,2,4]]
    @test isapprox(norm(X_ - X0)/norm(X0), 0f0; atol=1f-5)

    maxiter = 5
    h = 0.1f0
    err1 = zeros(Float32, maxiter)
    err2 = zeros(Float32, maxiter)

    print("\nGradient test ΔX for permute=$(permute), reverse=$(rev), logdet=$(logdet)\n")
    for j=1:maxiter
        f = lossf(HL, X0 + h*dX)[1]
        err1[j] = abs(f - f0)
        err2[j] = abs(f - f0 - h*dot(dX, gX))
        print(err1[j], "; ", err2[j], "\n")
        h = h/2f0
    end

    @test isapprox(err1[end] / (err1[1]/2^(maxiter-1)), 1f0; atol=1f0)
    @test isapprox(err2[end] / (err2[1]/4^(maxiter-1)), 1f0; atol=1f0)
end

function grad_test_layer(nx, ny, n_channel, n_hidden, batchsize, permute, logdet, rev)
    logdet ? lossf = loss_logdet : lossf = loss
    # Test for weights
    HL = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet)
    HL0 = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet)
    rev && (HL = reverse(HL))
    rev && (HL0 = reverse(HL0))

    X = glorot_uniform(nx, ny, n_channel, batchsize)
    HL.CL[1].RB.W1.data *= 4f0
    HL0.CL[1].RB.W1.data *= 4f0 # make weights larger
    HLini = deepcopy(HL0)
    dW = HL.CL[1].RB.W1.data - HL0.CL[1].RB.W1.data

    f0, gX, gW, X_ = lossf(HL0, X)
    @test isapprox(norm(X_ - X)/norm(X), 0f0; atol=1f-5)

    maxiter = 5
    h = 0.5f0
    err3 = zeros(Float32, maxiter)
    err4 = zeros(Float32, maxiter)

    print("\nGradient test weights for permute=$(permute), reverse=$(rev), logdet=$(logdet)\n")
    for j=1:maxiter
        HL0.CL[1].RB.W1.data = HLini.CL[1].RB.W1.data + h*dW
        f = lossf(HL0, X)[1]
        err3[j] = abs(f - f0)
        err4[j] = abs(f - f0 - h*dot(gW, dW))
        print(err3[j], "; ", err4[j], "\n")
        h = h/2f0
    end

    @test isapprox(err3[end] / (err3[1]/2^(maxiter-1)), 1f0; atol=1f1)
    @test isapprox(err4[end] / (err4[1]/4^(maxiter-1)), 1f0; atol=1f1)
end

# Loop over all possible permutation options
options = ["none", "lower", "both", "full"]

for permute in options
    for logdet in [true, false]
        for rev in [true, false]
            test_inv(nx, ny, n_channel, n_hidden, batchsize, permute, logdet, rev)
            grad_test_layer(nx, ny, n_channel, n_hidden, batchsize, permute, logdet, rev)
            grad_test_X(nx, ny, n_channel, n_hidden, batchsize, permute, logdet, rev)
        end
    end
end

###################################################################################################
# Jacobian-related tests

# Gradient test

# Initialization
logdet=true
# logdet=false
permute="full"
# permute="both"
# permute="none"
# permute="lower"
HL = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet)
θ = deepcopy(get_params(HL))
HL0 = CouplingLayerHINT(n_channel, n_hidden; permute=permute, logdet=logdet)
θ0 = deepcopy(get_params(HL0))
X = randn(Float32, nx, ny, n_channel, batchsize)

# Perturbation (normalized)
dθ = θ-θ0; dθ .*= norm.(θ0)./(norm.(dθ).+1f-10)
dX = randn(Float32, nx, ny, n_channel, batchsize); dX *= norm(X)/norm(dX)

# Jacobian eval
logdet ? ((dY, Y, lgdet, GNdθ) = HL.jacobian(dX, dθ, X)) : ((dY, Y) = HL.jacobian(dX, dθ, X))

# Test
print("\nJacobian test\n")
h = 0.1f0
maxiter = 5
err5 = zeros(Float32, maxiter)
err6 = zeros(Float32, maxiter)
for j=1:maxiter
    set_params!(HL, θ+h*dθ)
    logdet ? ((Y_loc, _) = HL.forward(X+h*dX)) : (Y_loc = HL.forward(X+h*dX))
    err5[j] = norm(Y_loc - Y)
    err6[j] = norm(Y_loc - Y - h*dY)
    print(err5[j], "; ", err6[j], "\n")
    global h = h/2f0
end

@test isapprox(err5[end] / (err5[1]/2^(maxiter-1)), 1f0; atol=1f1)
@test isapprox(err6[end] / (err6[1]/4^(maxiter-1)), 1f0; atol=1f1)

# Adjoint test

set_params!(HL, θ)
logdet ? ((dY, Y, _, _) = HL.jacobian(dX, dθ, X)) : ((dY, Y) = HL.jacobian(dX, dθ, X))
dY_ = randn(Float32, size(dY))
logdet ? ((dX_, dθ_, _, _) = HL.adjointJacobian(dY_, Y)) : ((dX_, dθ_, _) = HL.adjointJacobian(dY_, Y))
a = dot(dY, dY_)
b = dot(dX, dX_)+dot(dθ, dθ_)
@test isapprox(a, b; rtol=1f-3)


text:
module MapTransform

    mutable struct Relation

        lat      :: AbstractArray{Float64, 2}
        area     :: AbstractArray{Float64, 2}
        mask     :: AbstractArray{Float64, 2}

        lat_bnd  :: AbstractArray{Float64, 1}
        ∂a       :: AbstractArray{Float64, 1}
        A        :: Float64
        indices  :: Any

        function Relation(;
            lat     :: AbstractArray{Float64, 2},
            area    :: AbstractArray{Float64, 2},
            mask    :: AbstractArray{Float64, 2},
            lat_bnd :: AbstractArray{Float64, 1},
        )
            Δlat = lat_bnd[2:end] - lat_bnd[1:end-1]

            if any(Δlat .<= 0)
                throw(ErrorException("lat_bnd not monotonically increasing."))
            end

            indices = []
            mask_is_one = mask .== 1.0
            push!(indices, (lat_bnd[1] .<= lat .<= lat_bnd[2]) .& mask_is_one)
            for i = 2:length(lat_bnd)-1
                push!(indices, (lat_bnd[i] .< lat .<= lat_bnd[i+1]) .& mask_is_one )
            end

            ∂a = zeros(Float64, length(lat_bnd)-1)
            for i = 1:length(∂a)
                ∂a[i] = sum(area[indices[i]])
            end

            area_chk = sum(area[mask.==1.0])
            if abs((sum(∂a) - area_chk) / area_chk) > 1e-10
                print(abs((sum(∂a))), "; area_chk=", area_chk)
                throw(ErrorException("Sum area not equal."))
            end

            return new(
                lat,
                area,
                mask,
                lat_bnd,
                ∂a,
                sum(∂a),
                indices,
            )

        end

    end

    function f∂a(
        r :: Relation,
        f :: AbstractArray{Float64, 2},
    )
        f∂a = zeros(Float64, length(r.lat_bnd)-1)
        for i = 1:length(r.∂a)
            index = r.indices[i]
            f∂a[i] = sum(r.area[index] .* f[index])
        end

        return f∂a
    end

    function ∫∂a(
        r :: Relation,
        f :: AbstractArray{Float64, 2},
    )
        _∫f∂a = zeros(Float64, length(r.lat_bnd))
        _f∂a = f∂a(r, f)

        for i = 1:length(_f∂a)
            _∫f∂a[i+1] = _∫f∂a[i] + _f∂a[i]
        end

        return _∫f∂a

    end

    function transform(
        r :: Relation,
        f :: AbstractArray{Float64, 2},
    )
        _f∂a = f∂a(r, f)

        return _f∂a ./ r.∂a

    end

    function mean(
        r :: Relation,
        f :: AbstractArray{Float64, 2},
    )
        return ∫∂a(r, f)[end] / r.A
    end
end

max_stars_repo_name: meteorologytoday/CESM-diagnostic

max_issues_repo_name: meteorologytoday/CESM-diagnostic

max_forks_repo_name: meteorologytoday/CESM-diagnostic


text:
module Metropolis

export metropolis

using Distributions

function metropolis(; P::Function, x₀::Real = 0.0, Δ::Real = 1.0, Steps::Integer = 1000000)
    Aₙ = 0
    XList = Vector{Float64}(undef, Steps)
    x = x₀
    for n ∈ 1:Steps
        XList[n] = x
        y = x + rand(Uniform(-Δ, Δ))
        if rand() < P(y) / P(x)
            x = y
            Aₙ += 1
        end
    end
    return XList, Aₙ / Steps
end

end


text:
using Plots
using BenchmarkTools

include("random_sky.jl")

dir=false

mappa = zeros(NN, NN)
hit   = zeros(NN, NN)
N_baselines = 15

dataset_baselines = Array{Main.strip2.TOD_fake_pointing, 1}(undef, N_baselines)

for i = 1:N_baselines
    dataset = strip2.observe_sky(NN, conv_map; noise=true, direction_l_r=dir)
    m, h = strip2.binned_map(NN, dataset.time_order_data, dataset.pointing)
    mappa .+= m
    hit   .+= h
    dataset_baselines[i] = dataset
    global dir =  true ⊻ dir
end

tod = Array{Float64, 1}(undef, 0)
puntamenti = Array{Int64}(undef, 0)
for i in dataset_baselines
    append!(tod, i.time_order_data)
    append!(puntamenti, i.pointing[:,1])
end

for i in dataset_baselines
    append!(puntamenti, i.pointing[:,2])
end

puntamenti = reshape(puntamenti, NN*NN*10, 2)

#obs_map  = strip2.Plot_Sky(NN, pix_size, mappa ./ hit, cmin=-600, cmax=600)
#true_map = strip2.Plot_Sky(NN, pix_size, conv_map)

#plot(obs_map, true_map, layout = grid(2,1),size=(590,940))

map_des, a, b = strip2.destriper(NN, dataset_baselines, tolerance=1e-6)
map_d = reshape(map_des, NN, NN)
residui =  (map_d .- (sum(map_d)/(NN*NN))) .- (conv_map .- (sum(conv_map)/(NN*NN)))

#p2 = strip2.Plot_Sky(NN, pix_size, conv_map .- (sum(conv_map)/(NN*NN)), cmin = -400, cmax = 400, lab_title="T Map")
p3 = strip2.Plot_Sky(NN, pix_size, map_d .- (sum(map_d)/(NN*NN)) .-  (conv_map .- (sum(conv_map)/(NN*NN))), cmin=-20, cmax=20, lab_title="Destriped Map")
#p1 = strip2.Plot_Sky(NN, pix_size, map_d .- (sum(map_d)/(NN*NN)) .- (conv_map .- (sum(conv_map)/(NN*NN))) , cmin=-8, cmax=8, lab_title="Differences between destriped and true maps")

#plot(p1, p2, layout = grid(2,1), size=(590,940))

#strip2.Plot_Sky(NN, pix_size, map_d, cmin = -400, cmax = 400)


text:
"""
    addedge!(tree, parent, child, length=missing)

Create and add to `tree` a new edge from `parent` to `child`,
with an appropriate number.
"""
function addedge!(tree::Tree{T}, parent, child, length=missing) where T
    push!(tree.edge, Edge(parent, child, length))
    return nothing
end

const TR = Tree{T} where T

"""
    path2root(tree, node_index)
    path2root!(list_of_node_indices, tree, node_index)

Find the list of edges from the root to the `node_index`.
The second version modifies the list of edges in place and calls
itself until we find the root (a node with no parent).

Note: this function could run forever if the tree has a directed cycle,
but the function checks for no more edges in the path than edges in the tree
to avoid an infinite path cycling around itself.
"""
function path2root(tree::TR, i::Int)
    res = Edge[] # vector of Edges, of length 0
    path2root!(res, tree::TR, i::Int)
end
function path2root!(res, tree::TR, i::Int)
    pei = findfirst([e.child == i for e in tree.edge]) # parent edge index
    if pei === nothing # i is not the child of anyone: we are at the root
        return res     # ends the recursion
    end
    # pei was something, otherwise we would have quit the function
    parentedge = tree.edge[pei]
    pushfirst!(res, parentedge)
    # pushfirst! instead of push!, above, to have a path that starts at the root
    if length(res) > length(tree.edge)
        error("more nodes in the path to the root than nodes in the tree: $res")
    end
    path2root!(res, tree::TR, parentedge.parent::Int)
end


text:
using FFTW, Luxor, Images, TestImages, Colors

function fourierfill(w, h, F)
    @layer begin
        a = CartesianIndices(F)
        for i in a
            re, im = F[i].re, F[i].im
            sethue(LCHab(70, 80, rescale(re, -8, 8, 0, 255)))
            r, c = Tuple(i)
            pt = Point(-w/2 + c, -h/2 + r)
            t = mask(pt, O, w, h, easingfunction=easeoutquad)
            box(K * pt, K * t, K * t, :fill)
        end
        box(O, K , K , :fill) # it was a bug, Dave
    end
end

const w, h = 128, 128
const K = 16

function main(fname)
    img = testimage("cameraman")
    img1 = imresize(img, w, h)
    timg = fft(channelview(img1))
    # F = timg
    F = fftshift(timg)
    Drawing((K * w), (K * h), fname)
    origin()
    background("grey5")
    setline(0.25)

    @layer begin
        setopacity(0.5)
        fourierfill(w, h, F)
    end

    fontface("Archivo-Black")
    fontsize(450)
    sethue("orange")
    setline(5)

    # VIZCON
    textoutlines("VIZCON", O, halign=:center,:path)
    clip()
    fourierfill(w, h, F)
    clipreset()
    textoutlines("VIZCON", O, halign=:center,:stroke)

    # BERLIN
    textoutlines("BERLIN", O + (0, 340), halign=:center,:path)
    clip()
    fourierfill(w, h, F)
    clipreset()
    textoutlines("BERLIN", O + (0, 340), halign=:center,:stroke)

    # Julia
    @layer begin
        setopacity(1.0)
        translate(0, -650)
        scale(2.5)
        translate(-165, -125)
        julialogo()
    end

    # date
    sethue("white")
    fontsize(100)
    setopacity(0.75)
    texttrack("MARCH 12-16 2020", O + (-900, 600), 500, 100)

    # url
    fontsize(60)
    fontface("ArchivoNarrow-Bold")
    setopacity(1)
    te = textextents("Visit discourse.julialang.org/tag/plotting")
    texttrack("Visit discourse.julialang.org/tag/plotting", O + (-te[3], 900), 365, 60)

    finish()
    preview()
end

main("/tmp/vizcon.svg")

max_stars_repo_path: vizcon 2020/vizcon.jl

max_stars_repo_name: radovankavicky/graphics

max_issues_repo_path: vizcon 2020/vizcon.jl

max_issues_repo_name: radovankavicky/graphics

max_forks_repo_path: vizcon 2020/vizcon.jl

max_forks_repo_name: cormullion/graphics


text:
module DifferentialDynamicProgramming
using LinearTimeVaryingModelsBase, Requires, ValueHistories, LinearAlgebra, Statistics, Printf
const DEBUG = false # Set this flag to true in order to print debug messages

export QPTrace, boxQP, demoQP, iLQG,iLQGkl, demo_linear, demo_linear_kl, demo_pendcart, GaussianPolicy

eye(n) = Matrix{Float64}(I,n,n)

function __init__()
    @require Plots="91a5bcdd-55d7-5caf-9e0b-520d859cae80" begin
    @eval LinearAlgebra.adjoint(x::String) = x
    @eval function plotstuff_linear(x,u,cost,totalcost)
        p = Plots.plot(layout=(2,2))
        Plots.plot!(p,x', title="State Trajectories", xlabel="Time step",legend=false, subplot=1, show=false)
        Plots.plot!(p,cost,c=:black,linewidth=3, title="Cost", xlabel="Time step", subplot=2, show=false)
        Plots.plot!(p,u',title="Control signals", xlabel="Time step", subplot=3, show=false)
        Plots.plot!(p,totalcost,title="Total cost", xlabel="Iteration", subplot=4, show=false)
        display(Plots.current())
    end
    @eval function plotstuff_pendcart(x00, u00, x,u,cost00,cost,otrace)
        cp = Plots.plot(layout=(1,3))
        sp = Plots.plot(x00',title=["\$x_$(i)\$" for i=1:size(x00,1)]', lab="Simulation", layout=(2,2))
        Plots.plot!(sp,x', title=["\$x_$(i)\$" for i=1:size(x00,1)]', lab="Optimized", xlabel="Time step", legend=true)

        Plots.plot!(cp,cost00, title="Cost", lab="Simulation", subplot=2)
        Plots.plot!(cp,u', legend=true, title="Control signal",lab="Optimized", subplot=1)
        Plots.plot!(cp,cost[2:end], legend=true, title="Cost",lab="Optimized", xlabel="Time step", subplot=2, yscale=:log10)
        iters = sum(cost .> 0)
        filter!(x->x>0,cost)
        Plots.plot!(cp, get(otrace, :cost)[2], yscale=:log10,xscale=:log10, title="Total cost", xlabel="Iteration", legend=false, subplot=3)
        Plots.plot(sp,cp)
        display(Plots.current())
    end
end
end

dir(paths...) = joinpath(@__DIR__, "..", paths...)
include("boxQP.jl")
include("iLQG.jl")
include("iLQGkl.jl")
include("forward_pass.jl")
include("backward_pass.jl")
include("demo_linear.jl")
include("system_pendcart.jl")

function debug(x)
    DEBUG && printstyled(string(x),"\n", color=:blue)
end

end # module


text:
module Examples

end # module


text:
push!(LOAD_PATH, pwd())
using Day3_module

wire_paths = readlines("input.txt")
#wire_paths = ['R75,D30,R83,U83,L12,D49,R71,U7,L72',
#              'U62,R66,U55,R34,D71,R55,D58,R83']
#wire_paths = ['R98,U47,R26,D63,R33,U87,L62,D20,R33,U53,R51',
#              'U98,R91,D20,R16,D67,R40,U7,R15,U6,R7']

function get_steps(wire_path)
    return split(wire_path, ',')
end

t0 = time()
pos_a, pos_b = [get_positions(get_steps(path)) for path in wire_paths]
crossings = intersect(Set(pos_a), Set(pos_b))

function manhattan_distance(pos)
    sum(map(abs, pos))
end

manhattan_distances = [manhattan_distance(pos) for pos in crossings]
println(minimum(manhattan_distances))

t1 = time()
println(t1 - t0)

println(minimum(indexin(crossings, pos_b) + indexin(crossings, pos_a)))

t2 = time()
println(t2 - t1)

max_stars_repo_name: st--/advent-of-code

max_issues_repo_name: st--/advent-of-code

max_forks_repo_name: st--/advent-of-code


text:
using ImagineProcedures, ImagineInterface
using Test

rt = rigtemplate("ocpi-2")
isigs = rt[1:2]

#returns 2 ImagineSignals
f = x -> isigs
sg = SignalGenerator(f, (0,))

@test isequal(description(sg), "")
@test func(sg) == f
@test func_args(sg) == (0,)
@test call!(sg) == isigs

af = x->x[1]
ip = ImagineProcedure("p", sg, af)
@test isa(repr(ip), AbstractString)
@test isequal(description(ip), "p")
@test ImagineProcedures.signal_generator(ip) == sg
@test ImagineProcedures.analyze_function(ip) == af
@test outputs(ip) == isigs
@test process(ip) == isigs[1]


text:
# redundant
struct QMDPState{B,S}
    isstate::Bool
    b::Nullable{B}
    s::Nullable{S}

    function QMDPState{B,S}(isstate::Bool, x) where {B,S}
        if isstate
            return new(isstate, Nullable{B}(), Nullable(x))
        else
            return new(isstate, Nullable(x), Nullable{S}())
        end
    end
end

struct QMDPWrapper{M,B,S,A} <: MDP{QMDPState{B,S}, A}
    mdp::M
end

function QMDPWrapper(mdp::MDP, B::Type)
    S = state_type(mdp)
    A = POMDPs.action_type(mdp)
    return QMDPWrapper{typeof(mdp), B, S, A}(mdp)
end

state(m::QMDPWrapper{<:Any,B,S,<:Any}, s::S) where {B,S} = QMDPState{B,S}(true, s)
state(m::QMDPWrapper{<:Any,B,S,<:Any}, b::B) where {B,S} = QMDPState{B,S}(false, b)

function isterminal(m::QMDPWrapper, s::QMDPState)
    if s.isstate
        return isterminal(m, get(s.s))
    else
        return isterminal(m, get(s.b))
    end
end

function generate_sr(m::QMDPWrapper, s::QMDPState, a, rng::AbstractRNG)
    if s.isstate
        sp, r = generate_sr(m.mdp, get(s.s), a, rng)
        return state_type(m)(true, sp), r
    else
        sp, r = generate_sr(m.mdp, rand(rng, get(s.b)), a, rng)
        return state_type(m)(true, sp), r
    end
end

discount(m::QMDPWrapper) = discount(m.mdp)
function actions(m::QMDPWrapper, s::QMDPState)
    if s.isstate
        return actions(m.mdp, get(s.s))
    else
        # Gallium.@enter actions(m.mdp, get(s.b))
        return actions(m.mdp, get(s.b))
    end
end

struct GenQMDPSolver <: Solver
    solver
end

struct GenQMDPPolicy{P<:Policy, Q<:QMDPWrapper} <: Policy
    policy::P
    qmdp::Q
end

solve(sol::GenQMDPSolver, qmdp::QMDPWrapper) = GenQMDPPolicy(solve(sol.solver, qmdp), qmdp)

function action_info(p::GenQMDPPolicy, b)
    # XXX if this ever makes it into the toolbox, need to get the mdp some other way
    s = QMDPState{typeof(b), state_type(p.qmdp.mdp)}(false, b)
    return action_info(p.policy, s)
end

action(p::GenQMDPPolicy, b) = first(action_info(p, b))


text:
using LinearAlgebra
using PolynomialBasis
using ImplicitDomainQuadrature
using Revise
using CutCell
include("../../../../test/useful_routines.jl")

function bulk_modulus(l, m)
    return l + 2m / 3
end

function lame_lambda(k, m)
    return k - 2m / 3
end

function analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
    a = zeros(3, 3)
    a[1, 1] = inradius
    a[1, 2] = -inradius
    a[1, 3] = -1.0 / inradius
    a[2, 1] = 2 * (lc + mc)
    a[2, 2] = -2 * (ls + ms)
    a[2, 3] = 2ms / inradius^2
    a[3, 2] = 2(ls + ms)
    a[3, 3] = -2ms / outradius^2
    return a
end

function analytical_coefficient_rhs(ls, ms, theta0)
    r = zeros(3)
    Ks = bulk_modulus(ls, ms)
    r[2] = -Ks * theta0
    r[3] = Ks * theta0
    return r
end

struct AnalyticalSolution
    inradius::Any
    outradius::Any
    center::Any
    A1c::Any
    A1s::Any
    A2s::Any
    ls::Any
    ms::Any
    lc::Any
    mc::Any
    theta0::Any
    function AnalyticalSolution(inradius, outradius, center, ls, ms, lc, mc, theta0)
        a = analytical_coefficient_matrix(inradius, outradius, ls, ms, lc, mc)
        r = analytical_coefficient_rhs(ls, ms, theta0)
        coeffs = a \ r
        new(
            inradius,
            outradius,
            center,
            coeffs[1],
            coeffs[2],
            coeffs[3],
            ls,
            ms,
            lc,
            mc,
            theta0,
        )
    end
end

function radial_displacement(A::AnalyticalSolution, r)
    if r <= A.inradius
        return A.A1c * r
    else
        return A.A1s * r + A.A2s / r
    end
end

function (A::AnalyticalSolution)(x)
    relpos = x - A.center
    r = sqrt(relpos' * relpos)
    ur = radial_displacement(A, r)
    if ur ≈ 0.0
        [0.0, 0.0]
    else
        costheta = (x[1] - A.center[1]) / r
        sintheta = (x[2] - A.center[2]) / r
        u1 = ur * costheta
        u2 = ur * sintheta
        return [u1, u2]
    end
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

function pressure(stress)
    return -(stress[1, :] + stress[2, :] + stress[4, :]) / 3
end

function onleftboundary(x, L, W)
    return x[1] ≈ 0.0
end

function onbottomboundary(x, L, W)
    return x[2] ≈ 0.0
end

function onrightboundary(x, L, W)
    return x[1] ≈ L
end

function ontopboundary(x, L, W)
    return x[2] ≈ W
end

function onboundary(x, L, W)
    return onleftboundary(x, L, W) ||
           onbottomboundary(x, L, W) ||
           onrightboundary(x, L, W) ||
           ontopboundary(x, L, W)
end

function solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    theta0,
    penalty,
    analyticalsolution,
)

    bilinearforms = CutCell.BilinearForms(basis, cellquads, stiffness, cutmesh)
    interfacecondition = CutCell.incoherent_interface_condition(
        basis,
        interfacequads,
        stiffness,
        cutmesh,
        penalty,
    )

    displacementbc = CutCell.DisplacementCondition(
        analyticalsolution,
        basis,
        facequads,
        stiffness,
        cutmesh,
        x -> onboundary(x, width, width),
        penalty,
    )

    sysmatrix = CutCell.SystemMatrix()
    sysrhs = CutCell.SystemRHS()

    CutCell.assemble_bilinear_form!(sysmatrix, bilinearforms, cutmesh)
    CutCell.assemble_interface_condition!(sysmatrix, interfacecondition, cutmesh)
    CutCell.assemble_bulk_transformation_linear_form!(
        sysrhs,
        transfstress,
        basis,
        cellquads,
        cutmesh,
    )
    CutCell.assemble_incoherent_interface_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        interfacequads,
        cutmesh,
    )
    CutCell.assemble_penalty_displacement_bc!(sysmatrix, sysrhs, displacementbc, cutmesh)
    CutCell.assemble_penalty_displacement_transformation_rhs!(
        sysrhs,
        transfstress,
        basis,
        facequads,
        cutmesh,
        x -> onboundary(x, width, width),
    )

    matrix = CutCell.make_sparse(sysmatrix, cutmesh)
    rhs = CutCell.rhs(sysrhs, cutmesh)

    nodaldisplacement = matrix \ rhs

    return nodaldisplacement
end

function spatial_closest_points(refclosestpoints, refclosestcellids, mesh)
    dim, npts = size(refclosestpoints)
    spclosestpoints = zeros(dim, npts)
    for i = 1:npts
        cellmap = CutCell.cell_map(mesh, refclosestcellids[i])
        spclosestpoints[:, i] .= cellmap(refclosestpoints[:, i])
    end
    return spclosestpoints
end

function angular_position(points)
    cpoints = points[1, :] + im * points[2, :]
    return rad2deg.(angle.(cpoints))
end

K1, K2 = 247.0, 192.0    # Pa
mu1, mu2 = 126.0, 86.0   # Pa
lambda1 = lame_lambda(K1, mu1)
lambda2 = lame_lambda(K2, mu2)
stiffness = CutCell.HookeStiffness(lambda1, mu1, lambda2, mu2)

rho1 = 1.0   # Kg/m^3
rho2 = 1.0   # Kg/m^3

V1 = inv(rho1)
V2 = inv(rho2)
theta0 = -0.067
ΔG = -1022.0    # J/mol
molarmass = 0.147   # Kg/mol
diffG0 = ΔG / molarmass # J/Kg

width = 1.0
penaltyfactor = 1e2

polyorder = 3
numqp = required_quadrature_order(polyorder) + 2
nelmts = 21
center = [width / 2, width / 2]
inradius = width / 4
outradius = width

analyticalsolution =
    AnalyticalSolution(inradius, outradius, center, lambda1, mu1, lambda2, mu2, theta0)

transfstrain = CutCell.plane_transformation_strain(theta0)
transfstress = CutCell.plane_strain_transformation_stress(lambda1, mu1, theta0)

dx = width / nelmts
meanmoduli = 0.5 * (lambda1 + lambda2 + mu1 + mu2)
penalty = penaltyfactor / dx * meanmoduli

basis = TensorProductBasis(2, polyorder)
mesh = CutCell.Mesh([0.0, 0.0], [width, width], [nelmts, nelmts], basis)
levelset = InterpolatingPolynomial(1, basis)
levelsetcoeffs =
    CutCell.levelset_coefficients(x -> -circle_distance_function(x, center, inradius), mesh)

cutmesh = CutCell.CutMesh(levelset, levelsetcoeffs, mesh)
cellquads = CutCell.CellQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
interfacequads = CutCell.InterfaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)
facequads = CutCell.FaceQuadratures(levelset, levelsetcoeffs, cutmesh, numqp)

nodaldisplacement = solve_for_displacement(
    basis,
    cellquads,
    interfacequads,
    facequads,
    cutmesh,
    levelset,
    levelsetcoeffs,
    stiffness,
    theta0,
    penalty,
    analyticalsolution,
)

refseedpoints, spatialseedpoints, seedcellids =
    CutCell.seed_zero_levelset_with_interfacequads(interfacequads, cutmesh)
nodalcoordinates = CutCell.nodal_coordinates(cutmesh)

tol = 1e-8
boundingradius = 3.0
refclosestpoints, refclosestcellids, refgradients =
    CutCell.closest_reference_points_on_zero_levelset(
        nodalcoordinates,
        refseedpoints,
        spatialseedpoints,
        seedcellids,
        levelset,
        levelsetcoeffs,
        mesh,
        tol,
        boundingradius,
    )
spclosestpoints = spatial_closest_points(refclosestpoints,refclosestcellids,cutmesh)

relspatialpoints = spclosestpoints .- center
angularposition = angular_position(relspatialpoints)
sortidx = sortperm(angularposition)
angularposition = angularposition[sortidx]

refclosestpoints = refclosestpoints[:,sortidx]
refclosestcellids = refclosestcellids[sortidx]
refgradients = refgradients[:,sortidx]

invjac = CutCell.inverse_jacobian(cutmesh)
normals = diagm(invjac) * refgradients
CutCell.normalize_normals!(normals)

product_stress_at_cp = CutCell.product_stress_at_reference_points(
    refclosestpoints,
    refclosestcellids,
    basis,
    stiffness,
    transfstress,
    theta0,
    nodaldisplacement,
    cutmesh,
)
parent_stress_at_cp = CutCell.parent_stress_at_reference_points(
    refclosestpoints,
    refclosestcellids,
    basis,
    stiffness,
    nodaldisplacement,
    cutmesh,
)

productpressure = CutCell.pressure_at_points(product_stress_at_cp)
parentpressure = CutCell.pressure_at_points(parent_stress_at_cp)

productdevstress = CutCell.deviatoric_stress_at_points(product_stress_at_cp,productpressure)
parentdevstress = CutCell.deviatoric_stress_at_points(parent_stress_at_cp,parentpressure)

productdevnorm = stress_inner_product_over_points(productdevstress)
parentdevnorm = stress_inner_product_over_points(parentdevstress)

productnormaldevstress = normal_stress_component_over_points(productdevstress,normals)
parentnormaldevstress = normal_stress_component_over_points(parentdevstress,normals)

productspecificvolume = V1*(1.0 .- productpressure/K1)
parentspecificvolume = V2*(1.0 .- parentpressure/K2)

productp1 = productspecificvolume .* productpressure
productp2 = V1/(2K1)*(productpressure.^2)
productp3 = -(productspecificvolume .* productnormaldevstress)
productp4 = V1/(4mu1)*productdevnorm
productpotential = productp1+productp2+productp3+productp4

parentp1 = parentspecificvolume .* parentpressure
parentp2 = V2/(2K2)*(parentpressure.^2)
parentp3 = -(parentspecificvolume .* parentnormaldevstress)
parentp4 = V2/(4mu2)*parentdevnorm
parentpotential = parentp1+parentp2+parentp3+parentp4

potentialdifference = productpotential - parentpotential

folderpath = "examples/transformation-strain/cylindrical-particle/cylindrical-bc/"

# using PyPlot
# fig,ax = PyPlot.subplots()
# ax.plot(angularposition,productpressure,label="product")
# # ax.plot(angularposition,parentpressure,label="parent")
# ax.grid()
# ax.set_title("Pressure along interface")
# ax.set_ylabel("Pressure (Pa)")
# ax.set_xlabel("Angular position (deg)")
# ax.legend()
# fig.tight_layout()
# fig
# fig.savefig(folderpath*"pressure.png")

# using PyPlot
# fig,ax = PyPlot.subplots()
# ax.plot(angularposition,parentp1,label="term1")
# ax.plot(angularposition,parentp2,label="term2")
# ax.plot(angularposition,parentp3,label="term3")
# ax.plot(angularposition,parentp4,label="term4")
# ax.plot(angularposition,parentpotential,label="potential")
# ax.set_xlabel("Angular position (deg)")
# ax.set_ylabel("Energy density (J/Kg)")
# ax.set_title("Potential on parent side of interface")
# # ax.set_ylim(-1e6,2e6)
# ax.grid()
# ax.legend()
# fig.tight_layout()
# fig
# fig.savefig(folderpath*"parent-potential.png")
#
#
#
# fig,ax = PyPlot.subplots()
# ax.plot(angularposition,productp1,label="term1")
# ax.plot(angularposition,productp2,label="term2")
# ax.plot(angularposition,productp3,label="term3")
# ax.plot(angularposition,productp4,label="term4")
# ax.plot(angularposition,productpotential,label="potential")
# ax.set_xlabel("Angular position (deg)")
# ax.set_ylabel("Energy density (J/Kg)")
# # ax.set_ylim(-1e6,2e6)
# ax.set_title("Potential on product side of interface")
# ax.grid()
# ax.legend()
# fig.tight_layout()
# fig
# fig.savefig(folderpath*"product-potential.png")
#
#
#
fig,ax = PyPlot.subplots()
ax.plot(angularposition,potentialdifference)
ax.grid()
ax.set_xlabel("Angular position (deg)")
ax.set_ylabel("Potential difference (J/Kg)")
ax.set_title("Potential difference along interface circumference")
fig.tight_layout()
# ax.set_ylim(0.2,0.3)
fig
# fig.savefig(folderpath*"potential-difference.png")


text:
"""
    AbstractRuleset <: ModelParameters.AbstractModel

Abstract supertype for [`Ruleset`](@ref) objects and variants.
"""
abstract type AbstractRuleset <: AbstractModel end

# Getters
ruleset(rs::AbstractRuleset) = rs
rules(rs::AbstractRuleset) = rs.rules
settings(rs::AbstractRuleset) = rs.settings
boundary(rs::AbstractRuleset) = boundary(settings(rs))
proc(rs::AbstractRuleset) = proc(settings(rs))
opt(rs::AbstractRuleset) = opt(settings(rs))
cellsize(rs::AbstractRuleset) = cellsize(settings(rs))
timestep(rs::AbstractRuleset) = timestep(settings(rs))
radius(set::AbstractRuleset) = radius(rules(set))

Base.step(rs::AbstractRuleset) = timestep(rs)

# ModelParameters interface
Base.parent(rs::AbstractRuleset) = rules(rs)
ModelParameters.setparent!(rs::AbstractRuleset, rules) = rs.rules = rules
ModelParameters.setparent(rs::AbstractRuleset, rules) = @set rs.rules = rules

"""
    Rulseset <: AbstractRuleset

    Ruleset(rules...; kw...)
    Ruleset(rules, settings)

A container for holding a sequence of `Rule`s and simulation
details like boundary handing and optimisation.
Rules will be run in the order they are passed, ie. `Ruleset(rule1, rule2, rule3)`.

# Keywords

- `proc`: a [`Processor`](@ref) to specificy the hardware to run simulations on,
    like [`SingleCPU`](@ref), [`ThreadedCPU`](@ref) or [`CuGPU`](@ref) when
    KernelAbstractions.jl and a CUDA gpu is available.
- `opt`: a [`PerformanceOpt`](@ref) to specificy optimisations like
    [`SparseOpt`](@ref). Defaults to [`NoOpt`](@ref).
- `boundary`: what to do with boundary of grid edges.
    Options are `Remove()` or `Wrap()`, defaulting to [`Remove`](@ref).
- `cellsize`: size of cells.
- `timestep`: fixed timestep where this is required for some rules.
    eg. `Month(1)` or `1u"s"`.
"""
mutable struct Ruleset{S} <: AbstractRuleset
    # Rules in Ruleset are intentionally not type-stable.
    # But they are when rebuilt in a StaticRuleset later
    rules::Tuple{Vararg{<:Rule}}
    settings::S
end
Ruleset(rule1, rules::Rule...; kw...) = Ruleset((rule1, rules...); kw...)
Ruleset(rules::Tuple; kw...) = Ruleset(rules, SimSettings(; kw...))
Ruleset(rs::AbstractRuleset) = Ruleset(rules(rs), settings(rs))
function Ruleset(; rules=(), settings=nothing, kw...)
    settings1 = settings isa Nothing ? SimSettings(; kw...) : settings
    Ruleset(rules, settings1)
end

struct StaticRuleset{R<:Tuple,S} <: AbstractRuleset
    rules::R
    settings::S
end
StaticRuleset(rule1, rules::Rule...; kw...) = StaticRuleset((rule1, rules...); kw...)
StaticRuleset(rules::Tuple; kw...) = StaticRuleset(rules, SimSettings(; kw...))
StaticRuleset(rs::AbstractRuleset) = StaticRuleset(rules(rs), settings(rs))
function StaticRuleset(; rules=(), settings=nothing, kw...)
    settings1 = settings isa Nothing ? SimSettings(; kw...) : settings
    StaticRuleset(rules, settings1)
end

const SRuleset = StaticRuleset


text:
include("../intcode.jl")

orig = set_memory!(Tape(parse_data(readline(open("input.txt")))), 100000)

function getval(x, y)
    tape = copy(orig)
    @async run_program!(tape)
    take!(tape.outch)
    put!(tape.inch, x)
    take!(tape.outch)
    put!(tape.inch, y)
    out = take!(tape.outch)
    take!(tape.outch)
    out
end

function getview(xsize, ysize, xoffset, yoffset)
    res = zeros(Int64, xsize, ysize)
    for i=1:xsize
        for j=1:ysize
            res[i,j] = getval(i-1+xoffset,j-1+yoffset)
        end
    end
    res
end

function showview(res)
    for i=1:size(res,1)
        for j=1:size(res,2)
            print(res[i,j])
        end
        println()
    end
end

res1 = getview(100, 100, 0, 0)
println(sum(res1))


text:
# PlotlyJS Examples
using PlotlyJS

# NOTE: By default, here is the PlotyJS color palette. This
# may prove useful when hard-coding colors is desired.
# #1F77B4
# #FF7F0E
# #2CA02C
# #D62728
# #9575D2
# #8C564B
# #E377C0
# #7F7F7F
# #BCBD22
# #17BECF

function hist2{T<:Real}(x1::AbstractArray{T,1}; col = "#87CEEB", bargap = 0.05, alpha = 0.99, title = "", plot_bgcolor = "white")

    trace1 = histogram(x = x1,
                       opacity = alpha,
                       marker_color = col,
                       nbinsx = 50)

    layout = Layout(barmode = "overlay",
                    bargap = bargap,
                    title = title,
                    plot_bgcolor = plot_bgcolor,
                    xaxis = attr(gridcolor = "#E2E2E2"),
                    yaxis = attr(gridcolor = "#E2E2E2"))
    plot(trace1, layout)
end

## Example use
x = randn(50_000)
hist2(x)

function hist2{T<:Real, S<:Real}(x1::AbstractArray{T,1}, x2::AbstractArray{S,1}; bargap = 0.05, alpha = 0.8, plot_bgcolor = "white")

    trace1 = histogram(x = x1,
                       opacity = alpha,
                       marker_color = "#1F77B4",
                       nbinsx = 50)
    trace2 = histogram(x = x2,
                       opacity = alpha,
                       marker_color = "#FF7F0E",
                       nbinsx = 50)
    data = [trace1, trace2]

    layout = Layout(barmode = "overlay",
                    bargap = bargap,
                    plot_bgcolor = plot_bgcolor,
                    xaxis = attr(gridcolor = "#E2E2E2"),
                    yaxis = attr(gridcolor = "#E2E2E2"))

    plot(data, layout)
end

## Example use
hist2(randn(3000), randn(3000) .+ 2)

function scatterplot{T<:Real}(y1::AbstractArray{T,1}, col = "#87CEEB")
    trace1 = scatter(y = y1,
                     mode = "markers",
                     marker_color = col)

    layout1 = Layout(plot_bgcolor = "white",
                     xaxis = attr(gridcolor = "#E2E2E2"),
                     yaxis = attr(gridcolor = "#E2E2E2"))

    plot(trace1, layout1)
end

## Example use
# scatterplot(randn(1000));

function scatterplot{T<:Real, S<:Real}(x1::AbstractArray{T,1}, y1::AbstractArray{S,1}, col = "#87CEEB")
    trace1 = scatter(x = x1,
                     y = y1,
                     mode = "markers",
                     marker_color = col)

    layout1 = Layout(plot_bgcolor = "white",
                     xaxis = attr(gridcolor = "#E2E2E2"),
                     yaxis = attr(gridcolor = "#E2E2E2"))

    plot(trace1, layout1)
end

## Example use
x = randn(1000);
y = 2.0 + 1.5x .+ randn(1000);
#
scatterplot(x, y);

function scatterplot{T<:Real, S<:Real, U<:Real}(x1::AbstractArray{T,1}, y1::AbstractArray{S,1}, z1::AbstractArray{U,1}, col = "#87CEEB")
    trace1 = scatter3d(x = x1,
                     y = y1,
                     z = z1,
                     mode = "markers",
                     marker_color = col,
                     marker_size = 2)

    layout1 = Layout(scene = attr(
                     xaxis = attr(gridcolor = "#E2E2E2", backgroundcolor = "white",
                                  showbackground = true),
                     yaxis = attr(gridcolor = "#E2E2E2", backgroundcolor = "white",
                                  showbackground = true),
                     zaxis = attr(gridcolor = "#E2E2E2", backgroundcolor = "white",
                                  showbackground = true)))
    plot(trace1, layout1)
end

## Example use
n = 5_000
x2 = randn(n);
y2 = 1.5x2 .+ randn(n);
z2 = 0.5 + 1.1x2 .+ 0.5y2 .+ randn(n)

scatterplot(x2, y2, z2)


text:
export AbstractLearner, Approximator

import Flux
import Functors

abstract type AbstractLearner end

(L::AbstractLearner)(env) = env |> state |> send_to_device(L) |> L |> send_to_device(env)

Base.@kwdef mutable struct Approximator{M,O}
    model::M
    optimiser::O
end

Functors.functor(x::Approximator) = (model = x.model,), y -> Approximator(y.model, x.state)

(A::Approximator)(x) = A.model(x)

RLBase.optimise!(A::Approximator, gs) =
    Flux.Optimise.update!(A.optimiser, Flux.params(A), gs)

max_stars_repo_name: jbrea/ReinforcementLearning

max_issues_repo_name: jbrea/ReinforcementLearning

max_forks_repo_name: jbrea/ReinforcementLearning


text:
# The Computer Language Benchmarks Game
# https://salsa.debian.org/benchmarksgame-team/benchmarksgame/

# contributed by Jarret Revels and Alex Arslan
# based on the Java version

module NBody

using Printf
using LinearAlgebra

# Constants
const solar_mass = 4π^2
const days_per_year = 365.24

# Use a 4-tuple here to get better SIMD instructions.
# This is a space-time tradeoff, but this benchmark is well within L1 cache limits.
struct Vec3
    x::NTuple{4, Float64}
end
Vec3(x, y, z) = Vec3((x,y,z,0.0))
Base.:/(v::Vec3, n::Number) = Vec3(1/n .* v.x)
Base.:*(v::Vec3, n::Number) = Vec3(n .* v.x)
Base.:-(v1::Vec3, v2::Vec3) = Vec3(v1.x .- v2.x)
Base.:+(v1::Vec3, v2::Vec3) = Vec3(v1.x .+ v2.x)
# Todo, prettify
squarednorm(v1::Vec3) = v1.x[1]^2 + v1.x[2]^2 + v1.x[3]^2
Base.muladd(x::Vec3, y::Number, z::Vec3) = Vec3(muladd.(x.x, y, z.x))

# A heavenly body in the system
mutable struct Body
    pos::Vec3
    v::Vec3
    mass::Float64
end

function offset_momentum!(b::Body, p::Vec3)
    b.v -= p / solar_mass
end

function init_sun!(bodies::Vector{Body})
    p = Vec3(0.0, 0.0, 0.0)
    for b in bodies
        p += b.v * b.mass
    end
    offset_momentum!(bodies[1], p)
end

function advance(bodies::Vector{Body}, dt::Number)
    @inbounds for i = 1:length(bodies)
        bi = bodies[i]
        for j = i+1:length(bodies)
            bj = bodies[j]
            delta = bi.pos - bj.pos
            dsq = squarednorm(delta)
            distance = sqrt(dsq)
            mag = dt / (dsq * distance)
            bi.v = muladd(delta, -(bj.mass * mag), bi.v)
            bj.v = muladd(delta, (bi.mass * mag), bj.v)
        end
    end

    for b in bodies
        b.pos = muladd(b.v, dt, b.pos)
    end
end

function energy(bodies::Vector{Body})
    e = 0.0
    @inbounds for i = 1:length(bodies)
        bi = bodies[i]
        e += 0.5 * bi.mass * squarednorm(bi.v)
        for j = i+1:length(bodies)
            bj = bodies[j]
            delta = bi.pos - bj.pos
            distance = sqrt(squarednorm(delta))
            dinv = 1.0 / distance
            e = muladd((bi.mass * bj.mass), -dinv, e)
        end
    end
    return e
end

function perf_nbody(N::Int=1000)
    jupiter = Body( Vec3(4.84143144246472090e+00,                  # pos[1] = x
                         -1.16032004402742839e+00,                 # pos[2] = y
                         -1.03622044471123109e-01),                # pos[3] = z
                   Vec3(1.66007664274403694e-03 * days_per_year,   # v[1] = vx
                        7.69901118419740425e-03 * days_per_year,   # v[2] = vy
                        -6.90460016972063023e-05 * days_per_year), # v[3] = vz
                   9.54791938424326609e-04 * solar_mass)       # mass

    saturn = Body(Vec3(8.34336671824457987e+00,
                       4.12479856412430479e+00,
                      -4.03523417114321381e-01),
                  Vec3(-2.76742510726862411e-03 * days_per_year,
                       4.99852801234917238e-03 * days_per_year,
                       2.30417297573763929e-05 * days_per_year),
                  2.85885980666130812e-04 * solar_mass)

    uranus = Body(Vec3(1.28943695621391310e+01,
                    -1.51111514016986312e+01,
                    -2.23307578892655734e-01),
                  Vec3(2.96460137564761618e-03 * days_per_year,
                       2.37847173959480950e-03 * days_per_year,
                       -2.96589568540237556e-05 * days_per_year),
                  4.36624404335156298e-05 * solar_mass)

    neptune = Body(Vec3(1.53796971148509165e+01,
                        -2.59193146099879641e+01,
                        1.79258772950371181e-01),
                   Vec3(2.68067772490389322e-03 * days_per_year,
                        1.62824170038242295e-03 * days_per_year,
                       -9.51592254519715870e-05 * days_per_year),
                   5.15138902046611451e-05 * solar_mass)

    sun = Body(Vec3(0.0, 0.0, 0.0), Vec3(0.0, 0.0, 0.0), solar_mass)

    bodies = [sun, jupiter, saturn, uranus, neptune]

    init_sun!(bodies)
    @printf("%.9f\n", energy(bodies))
    for i = 1:N
        advance(bodies, 0.01)
    end
    @printf("%.9f\n", energy(bodies))
end

end # module

n = parse(Int,ARGS[1])
NBody.perf_nbody(n)


text:
@testset "862.shortest-subarray-with-sum-at-least-k.jl" begin
    @test shortest_subarray([2], 2) == 1
    @test shortest_subarray([2, -1, 2], 3) == 3
    @test shortest_subarray([1, 2], 4) == -1
    @test shortest_subarray([25, -23, 21, 24, 27, 9, -14, -11, 27, -6, 12, 9, -15, 1, -30,
                             -16, -6, 11, 19, -14, 21, -12, 0, -28, -24, 16, -1, -18, 22,
                             12, 4, 8, -11, 28, -18, -9, -14, 5, 26, -10, -29, 21, 21, -3,
                             4, -3, -19, 28, -10, 13, 30, 14, 23, 25, 10, 19, 25, -13, -28,
                             -2, -22, -13, 26, -25, 3, 1, 5, -18, 18, -28, 8, 26, -21, -26,
                             27, 8, 15, -4, 21, -1, 27, -1, -24, -7, 15, 25, 28, -21, -9,
                             21], 70) == 3
    @test shortest_subarray([25, -23, 21, 24, 27, 9, -14, -11, 27, -6, 12, 9, -15, 1, -30,
                             -16, -6, 11, 19, -14, 21, -12, 0, -28, -24, 16, -1, -18, 22,
                             12, 4, 8, -11, 28, -18, -9, -14, 5, 26, -10, -29, 21, 21, -3,
                             4, -3, -19, 28, -10, 13, 30, 14, 23, 25, 10, 19, 25, -13, -28,
                             -2, -22, -13, 26, -25, 3, 1, 5, -18, 18, -28, 8, 26, -21, -26,
                             27, 8, 15, -4, 21, -1, 27, -1, -24, -7, 15, 25, 28, -21, -9,
                             21], 86) == 4
    @test shortest_subarray([25, -23, 21, 24, 27, 9, -14, -11, 27, -6, 12, 9, -15, 1, -30,
                             -16, -6, 11, 19, -14, 21, -12, 0, -28, -24, 16, -1, -18, 22,
                             12, 4, 8, -11, 28, -18, -9, -14, 5, 26, -10, -29, 21, 21, -3,
                             4, -3, -19, 28, -10, 13, 30, 14, 23, 25, 10, 19, 25, -13, -28,
                             -2, -22, -13, 26, -25, 3, 1, 5, -18, 18, -28, 8, 26, -21, -26,
                             27, 8, 15, -4, 21, -1, 27, -1, -24, -7, 15, 25, 28, -21, -9,
                             21], 30) == 1
end


text:
import Pkg
Pkg.instantiate()

using ArgParse
import JSON

### analysis functions

function groups(data; init=Tuple{String,Any}[])
  if "events" ∈ keys(data)
    events = data["events"]
    n = length(init)
    for e ∈ events
      groups(e; init)
    end
    length(init) == n && push!(init, (data["group"], events))
  end
  init
end

function findtrace(cache, ev)
  if "threadID" ∈ keys(ev)
    id = ev["threadID"]
    id ∈ keys(cache) && return cache[id]
  end
  if "stimulus" ∈ keys(ev)
    id = ev["stimulus"]["messageID"]
    id ∈ keys(cache) && return cache[id]
  end
  0
end

function addtocache!(cache, ev, trace)
  if "threadID" ∈ keys(ev)
    id = ev["threadID"]
    cache[id] = trace
  end
  if "stimulus" ∈ keys(ev)
    id = ev["stimulus"]["messageID"]
    cache[id] = trace
  end
  if "response" ∈ keys(ev)
    id = ev["response"]["messageID"]
    cache[id] = trace
  end
end

function splitcomponent(s)
  m = match(r"^(?<name>.*)::(?<clazz>[^/]*)/?(?<node>.*)$", s)
  m === nothing && return "", "", ""
  m[:name], m[:clazz], m[:node]
end

function prettymessage(msg)
  clazz = msg["clazz"]
  p = findlast('.', clazz)
  p === nothing || (clazz = clazz[p+1:end])
  sender = "sender" ∈ keys(msg) ? msg["sender"] : ""
  recipient = "recipient" ∈ keys(msg) ? msg["recipient"] : ""
  "$clazz ⟦ $sender → $recipient ⟧"
end

function analyze(events)
  traces = Tuple{Int64,String,Vector{Any}}[]
  cache = Dict{String,Int64}()
  for ev ∈ events
    trace = findtrace(cache, ev)
    if trace == 0
      cname, cclazz, cnode = splitcomponent(ev["component"])
      cnode != "" && (cnode = "[$cnode] ")
      desc = cnode * cname
      if "stimulus" ∈ keys(ev)
        desc = cnode * prettymessage(ev["stimulus"])
      elseif "response" ∈ keys(ev)
        desc = cnode * prettymessage(ev["response"])
      end
      push!(traces, (ev["time"], desc, [ev]))
      trace = length(traces)
    else
      push!(traces[trace][3], ev)
    end
    addtocache!(cache, ev, trace)
  end
  traces
end

function capture!(actors, msgs, origin, node, agent, t, msg, stimulus)
  "recipient" ∉ keys(msg) && return
  recipient = msg["recipient"]
  if startswith(recipient, '#')
    stimulus || return
    recipient = agent
  end
  onode = node
  okey = msg["messageID"] * "→" * recipient
  if okey ∈ keys(origin)
    n = origin[okey]
    onode == n && return
    onode = n
  end
  stimulus && "sender" ∉ keys(msg) && return
  sender = "sender" ∈ keys(msg) ? msg["sender"] : agent
  if node != ""
    sender = sender * "/" * onode
    recipient = recipient * "/" * node
  end
  clazz = msg["clazz"]
  if clazz == "org.arl.fjage.Message"
    clazz = msg["performative"]
  else
    p = findlast('.', clazz)
    p === nothing || (clazz = clazz[p+1:end])
  end
  origin[okey] = node
  sender == recipient && return
  push!(actors, (onode, sender))
  push!(actors, (node, recipient))
  push!(msgs, (t, clazz, sender, recipient))
end

function sequence(events)
  actors = Tuple{String,String}[]
  msgs = Tuple{Int64,String,String,String}[]
  origin = Dict{String,String}()
  for ev ∈ events
    cname, cclazz, cnode = splitcomponent(ev["component"])
    "stimulus" ∈ keys(ev) && capture!(actors, msgs, origin, cnode, cname, ev["time"], ev["stimulus"], true)
    "response" ∈ keys(ev) && capture!(actors, msgs, origin, cnode, cname, ev["time"], ev["response"], false)
  end
  unique!(sort!(actors; by = x -> x[1]))
  actors, msgs
end

function mermaid(actors, msgs)
  println("sequenceDiagram")
  for a ∈ actors
    id = replace(a[2], '/' => '_')
    println("  participant $id as $(a[2])")
  end
  for m ∈ msgs
    id1 = replace(m[3], '/' => '_')
    id2 = replace(m[4], '/' => '_')
    print("  $id1")
    m[2] == "AGREE" && print('-')
    println("->>$id2: $(m[2])")
  end
end

### main

aps = ArgParseSettings()
@add_arg_table! aps begin
  "--group", "-g"
    help = "event group number"
    arg_type = Int
    default = 0
  "--trace", "-t"
    help = "trace number"
    arg_type = Int
    default = 0
  "filename"
    help = "JSON trace file"
    default = "trace.json"
end
args = parse_args(ARGS, aps)

if !isfile(args["filename"])
  println(args["filename"], " not found")
  exit(1)
end

grps = groups(JSON.parsefile(args["filename"]))

if length(grps) < 1
  println("no traces found")
  exit(1)
end

g = args["group"]
length(grps) == 1 && g == 0 && (g = 1)

if g == 0
  println("Specify an event group:")
  for (i, g1) ∈ enumerate(grps)
    println(i, ": ", g1[1], " (", length(g1[2]), " events)")
  end
  exit(1)
end

if g < 1 || length(grps) < g
  println("invalid event group number")
  exit(1)
end

traces = analyze(grps[g][2])

t = args["trace"]
length(traces) == 1 && t == 0 && (t = 1)

if t == 0
  println("Specify a trace:")
  for (i, t1) ∈ enumerate(traces)
    println(i, ": ", t1[1], " ", t1[2], " (", length(t1[3]), " events)")
  end
  exit(1)
end

if t < 1 || length(traces) < t
  println("invalid trace number")
  exit(1)
end

mermaid(sequence(traces[t][3])...)


text:
include("../src/CspExamples.jl")
using Test

fill_channel(ch, input) = for v in input
    put!(ch, v)
end

@testset "S31_COPY" begin
    @testset "Copy pre-filled buffer" begin
        # The input channel is a 1-character buffer as described in the paper.
        west = Channel(ch->(fill_channel(ch, "hello")), csize=1, ctype=Char, sticky=false)
        east = Channel(ch->CspExamples.S31_COPY(west, ch), ctype=Char, sticky=false)
        @test String(collect(east)) == "hello"
    end
    @testset "test incremental concurrent copying" begin
        # TODO: Should these be unbuffered channels?
        west = Channel{Char}(1)
        east = Channel(ch->CspExamples.S31_COPY(west, ch), ctype=Char, csize=Inf, sticky=false)
        put!(west, 'w'); put!(west, 'o'); put!(west, 'r'); put!(west, 'l'); put!(west, 'd');
        close(west)
        @test String(collect(east)) == "world"
    end
end

@testset "S32_SQUASH" begin
    # The input channel is a 1-character buffer as described in the paper.
    west = Channel(csize=1, ctype=Char, sticky=false) do ch
        fill_channel(ch, "hello*to**the*world")
    end
    east = Channel(ch->CspExamples.S32_SQUASH(west, ch), ctype=Char, sticky=false)
    @test String(collect(east)) == "hello*to↑the*world"
end

# Convenience function to create and fill a channel.
make_fill_close_chnl(input; csize=0, ctype=eltype(input)) = Channel(csize=csize, ctype=ctype, sticky=false) do ch
    fill_channel(ch, input)
end

@testset "S32_SQUASH_EXT" begin
    # Test 1 * at end
    west = make_fill_close_chnl("hello**world*")
    east = Channel(ch->CspExamples.S32_SQUASH_EXT(west, ch), ctype=Char, sticky=false)
    @test String(collect(east)) == "hello↑world*"

    # Test 3 *s at end
    west = make_fill_close_chnl("hello**world***", csize=1)
    east = Channel(ch->CspExamples.S32_SQUASH_EXT(west, ch), ctype=Char, sticky=false)
    @test String(collect(east)) == "hello↑world↑*"
end

@testset "S33_DISASSEMBLE" begin
    cardfile = Channel(ctype=String, sticky=false) do ch
        for s in ["hey", "buddy"]; put!(ch, s); end
    end
    @test "hey buddy " ==
            String(collect(Channel(X->CspExamples.S33_DISASSEMBLE(cardfile, X), ctype=Char, sticky=false)))
end

@testset "S34_ASSEMBLE" begin
    linelength = 3
    inputstr = rand('a':'z', linelength*2)
    X = Channel(ctype=Char, sticky=false) do ch
        # Pass two characters too few.
        for s in inputstr[1:end-2]; put!(ch, s); end
    end
    lineprinter = Channel(Inf)
    CspExamples.S34_ASSEMBLE(X, lineprinter, linelength)
    close(lineprinter)

    # Expect to get back the inputstr split into 125-character lines
    expected = inputstr
    expected[end-1:end] .= ' '
    expected = [String(expected[1+(i-1)*linelength:i*linelength]) for i in 1:2]
    @test collect(lineprinter) == expected
end

@testset "S35_Reformat" begin
    reformatted = Channel(sticky=false) do ch
        CspExamples.S35_Reformat(make_fill_close_chnl(["hello", "world"]), ch, 4)
    end
    @test collect(reformatted) == ["hell", "o wo", "rld "]
end
@testset "S35_Reformat2" begin
    reformatted = Channel(sticky=false) do ch
        CspExamples.S35_Reformat2(make_fill_close_chnl(["hello", "world"]), ch, 4)
    end
    @test collect(reformatted) == ["hell", "o wo", "rld "]
end
@testset "S35_Reformat_non_concurrent" begin
    # While the implementation of S35_Reformat_non_concurrent is more convoluted than
    # the concurrent version, the callsite is of course simpler because we don't have to
    # consider concurrency at all.
    #
    # A best-of-both-worlds approach might be to expose a "normal" interface like this one,
    # as a wrapper around a concurrent implementation.
    reformatted = CspExamples.S35_Reformat_non_concurrent(["hello", "world"], 4)
    @test reformatted == ["hell", "o wo", "rld "]
end

@testset "S36_Conway" begin
    reformatted = Channel(sticky=false) do ch
        CspExamples.S36_Conway(make_fill_close_chnl(["**hello", "world**"]), ch, 4)
    end
    @test collect(reformatted) == ["↑hel", "lo w", "orld", "↑   "]
end


text:
include("setting.jl")

# jacobian of lorenz96
function jac_lorenz(x::Vector{Float64})::Matrix{Float64}
    # jac[i, j] = -1 when j = i
    # jac[i, j] = x[i-1] when j = i+1
    # jac[i, j] = x[i+1] - x[i-2] when j = i-1
    # jac[i, j] = -x[i-1] when j = i-2
    len = length(x)
    lenm = len - 1;
    jac = zeros(len, len)

    for i = 1:len
        ip1 = (i+1+lenm)%len+1
        im1 = (i-1+lenm)%len+1
        im2 = (i-2+lenm)%len+1
        @inbounds jac[i, i] = -1
        @inbounds jac[i, ip1] = x[im1]
        @inbounds jac[i, im1] = x[ip1] - x[im2]
        @inbounds jac[i, im2] = -x[im1]
    end
    jac
end

a = repeat([1., 3, 6, 4, 5], 3)
da = repeat([.2, .2, .2, .2, .2], 3)
b = a + da
fa = lorenz96(a, Missing, Missing)
fb = lorenz96(b, Missing, Missing)
jfa = jac_lorenz(a) * da
display(jac_lorenz(a))
println()
@show fa + jfa - fb
@assert all(abs.(fa+jfa-fb) .< 1e-6)

max_stars_repo_name: Typas/Data-Assimilation-Project

max_issues_repo_name: Typas/Data-Assimilation-Project

max_forks_repo_name: Typas/Data-Assimilation-Project


text:
function ack(m,n)
    if m == 0
        return n + 1
    elseif n == 0
        return ack(m-1,1)
    else
        return ack(m-1,ack(m,n-1))
    end
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode