text:
using ModelingToolkit, Unitful, OrdinaryDiffEq, DiffEqJump, IfElse
using Test
MT = ModelingToolkit
@parameters τ [unit = u"ms"] γ
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
D = Differential(t)

#This is how equivalent works:
@test MT.equivalent(u"MW", u"kJ/ms")
@test !MT.equivalent(u"m", u"cm")
@test MT.equivalent(MT.get_unit(P^γ), MT.get_unit((E / τ)^γ))

# Basic access
@test MT.get_unit(t) == u"ms"
@test MT.get_unit(E) == u"kJ"
@test MT.get_unit(τ) == u"ms"
@test MT.get_unit(γ) == MT.unitless
@test MT.get_unit(0.5) == MT.unitless
@test MT.get_unit(MT.SciMLBase.NullParameters()) == MT.unitless

# Prohibited unit types
@parameters β [unit = u"°"] α [unit = u"°C"] γ [unit = 1u"s"]
@test_throws MT.ValidationError MT.get_unit(β)
@test_throws MT.ValidationError MT.get_unit(α)
@test_throws MT.ValidationError MT.get_unit(γ)

# Non-trivial equivalence & operators
@test MT.get_unit(τ^-1) == u"ms^-1"
@test MT.equivalent(MT.get_unit(D(E)), u"MW")
@test MT.equivalent(MT.get_unit(E / τ), u"MW")
@test MT.get_unit(2 * P) == u"MW"
@test MT.get_unit(t / τ) == MT.unitless
@test MT.equivalent(MT.get_unit(P - E / τ), u"MW")
@test MT.equivalent(MT.get_unit(D(D(E))), u"MW/ms")
@test MT.get_unit(IfElse.ifelse(t > t, P, E / τ)) == u"MW"
@test MT.get_unit(1.0^(t / τ)) == MT.unitless
@test MT.get_unit(exp(t / τ)) == MT.unitless
@test MT.get_unit(sin(t / τ)) == MT.unitless
@test MT.get_unit(sin(1u"rad")) == MT.unitless
@test MT.get_unit(t^2) == u"ms^2"

eqs = [D(E) ~ P - E / τ
       0 ~ P]
@test MT.validate(eqs)
@named sys = ODESystem(eqs)

@test !MT.validate(D(D(E)) ~ P)
@test !MT.validate(0 ~ P + E * τ)

# Array variables
@variables t [unit = u"s"] x[1:3](t) [unit = u"m"]
@parameters v[1:3]=[1, 2, 3] [unit = u"m/s"]
D = Differential(t)
eqs = D.(x) .~ v
ODESystem(eqs, name = :sys)

# Difference equation
@parameters t [unit = u"s"] a [unit = u"s"^-1]
@variables x(t) [unit = u"kg"]
δ = Differential(t)
D = Difference(t; dt = 0.1u"s")
eqs = [
    δ(x) ~ a * x,
]
de = ODESystem(eqs, t, [x], [a], name = :sys)

# Nonlinear system
@parameters a [unit = u"kg"^-1]
@variables x [unit = u"kg"]
eqs = [
    0 ~ a * x,
]
@named nls = NonlinearSystem(eqs, [x], [a])

# SDE test w/ noise vector
@parameters τ [unit = u"ms"] Q [unit = u"MW"]
@variables t [unit = u"ms"] E(t) [unit = u"kJ"] P(t) [unit = u"MW"]
D = Differential(t)
eqs = [D(E) ~ P - E / τ
       P ~ Q]

noiseeqs = [0.1u"MW",
    0.1u"MW"]
@named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])

# With noise matrix
noiseeqs = [0.1u"MW" 0.1u"MW"
            0.1u"MW" 0.1u"MW"]
@named sys = SDESystem(eqs, noiseeqs, t, [P, E], [τ, Q])

# Invalid noise matrix
noiseeqs = [0.1u"MW" 0.1u"MW"
            0.1u"MW" 0.1u"s"]
@test !MT.validate(eqs, noiseeqs)

# Non-trivial simplifications
@variables t [unit = u"s"] V(t) [unit = u"m"^3] L(t) [unit = u"m"]
@parameters v [unit = u"m/s"] r [unit = u"m"^3 / u"s"]
D = Differential(t)
eqs = [D(L) ~ v,
    V ~ L^3]
@named sys = ODESystem(eqs)
sys_simple = structural_simplify(sys)

eqs = [D(V) ~ r,
    V ~ L^3]
@named sys = ODESystem(eqs)
sys_simple = structural_simplify(sys)

@variables V [unit = u"m"^3] L [unit = u"m"]
@parameters v [unit = u"m/s"] r [unit = u"m"^3 / u"s"] t [unit = u"s"]
eqs = [V ~ r * t,
    V ~ L^3]
@named sys = NonlinearSystem(eqs, [V, L], [t, r])
sys_simple = structural_simplify(sys)

eqs = [L ~ v * t,
    V ~ L^3]
@named sys = NonlinearSystem(eqs, [V, L], [t, r])
sys_simple = structural_simplify(sys)

#Jump System
@parameters β [unit = u"(mol^2*s)^-1"] γ [unit = u"(mol*s)^-1"] t [unit = u"s"] jumpmol [
    unit = u"mol",
]
@variables S(t) [unit = u"mol"] I(t) [unit = u"mol"] R(t) [unit = u"mol"]
rate₁ = β * S * I
affect₁ = [S ~ S - 1 * jumpmol, I ~ I + 1 * jumpmol]
rate₂ = γ * I
affect₂ = [I ~ I - 1 * jumpmol, R ~ R + 1 * jumpmol]
j₁ = ConstantRateJump(rate₁, affect₁)
j₂ = VariableRateJump(rate₂, affect₂)
js = JumpSystem([j₁, j₂], t, [S, I, R], [β, γ], name = :sys)

affect_wrong = [S ~ S - jumpmol, I ~ I + 1]
j_wrong = ConstantRateJump(rate₁, affect_wrong)
@test_throws MT.ValidationError JumpSystem([j_wrong, j₂], t, [S, I, R], [β, γ], name = :sys)

rate_wrong = γ^2 * I
j_wrong = ConstantRateJump(rate_wrong, affect₂)
@test_throws MT.ValidationError JumpSystem([j₁, j_wrong], t, [S, I, R], [β, γ], name = :sys)

# mass action jump tests for SIR model
maj1 = MassActionJump(2 * β / 2, [S => 1, I => 1], [S => -1, I => 1])
maj2 = MassActionJump(γ, [I => 1], [I => -1, R => 1])
@named js3 = JumpSystem([maj1, maj2], t, [S, I, R], [β, γ])

#Test unusual jump system
@parameters β γ t
@variables S(t) I(t) R(t)

maj1 = MassActionJump(2.0, [0 => 1], [S => 1])
maj2 = MassActionJump(γ, [S => 1], [S => -1])
@named js4 = JumpSystem([maj1, maj2], t, [S], [β, γ])


text:
using FinElt

include("params.jl")
successive_refine("rect", 2, refinements, 0.4)

ext: jl


text:
##########
##########
# Get environment ready
##########
##########
# AUTHOR: Cole B. Brookson
# DATE OF CREATION: 2022-01-06
##########
##########

# simple set up ================================================================

import Pkg
Pkg.instantiate()

max_stars_repo_name: colebrookson/trait-based-rewiring

max_issues_repo_name: colebrookson/trait-based-rewiring

max_forks_repo_name: colebrookson/trait-based-rewiring


text:
function cross_attack(occ::UInt64, ui::UInt64)
    free_squares = Array{UInt64,1}()
    edge_squares = Array{UInt64,1}()

    free_squares, edge_squares = slide(free_squares, edge_squares,
        occ, ui, DIAGONALS, 9)
    free_squares, edge_squares = slide(free_squares, edge_squares,
        occ, ui, ANTIDIAGONALS, 7)

    return free_squares, edge_squares
end

function cross_attack(free_squares::Array{UInt64,1},
    edge_squares::Array{UInt64,1}, occ::UInt64, ui::UInt64)

    free_squares, edge_squares = slide(free_squares, edge_squares,
        occ, ui, DIAGONALS, 9)
    free_squares, edge_squares = slide(free_squares, edge_squares,
        occ, ui, ANTIDIAGONALS, 7)

    return free_squares, edge_squares
end

function gen_all_diagonal_masks()
    diag_masks = Dict{UInt64, UInt64}()
    for i = 1:64
        ui = INT2UINT[i]
        diag_mask = EMPTY
        for j = 1:15
            if ui & DIAGONALS[j] != EMPTY
                diag_mask |= DIAGONALS[j]
                break
            end
        end
        for k = 1:15
            if ui & ANTIDIAGONALS[k] != EMPTY
                diag_mask |= ANTIDIAGONALS[k]
                break
            end
        end
        diag_masks[ui] = diag_mask
    end
    return diag_masks
end
const DIAG_MASKS = gen_all_diagonal_masks()


text:
using Documenter
using EzXMLPatched

makedocs(
    sitename="EzXMLPatched.jl",
    modules=[EzXMLPatched],
    pages=["index.md", "manual.md", "reference.md", "devnotes.md"])

deploydocs(
    repo="github.com/bicycle1885/EzXMLPatched.jl.git",
    target="build")


text:
module SmolkaFormat

using CairoMakie
include("themes.jl")
export monohalfpage

using Colors, ColorSchemes
include("colorschemes.jl")
include("colors.jl")
export cmap_tum, cmap_tum2
export colors_tum

end


text:
function _parse_julia_file(filename::String)
    return Base.parse_input_line(read(filename, String); filename=filename)
end

function lowercase_drive(a)
    if length(a) >= 2 && a[2]==':'
        return lowercase(a[1]) * a[2:end]
    else
        return a
    end
end

max_stars_repo_name: caleb-allen/julia-vscode

max_issues_repo_name: caleb-allen/julia-vscode

max_forks_repo_name: caleb-allen/julia-vscode


text:
# By default, Julia/LLVM does not use fused multiply-add operations (FMAs).
# Since these FMAs can increase the performance of many numerical algorithms,
# we need to opt-in explicitly.
# See https://ranocha.de/blog/Optimizing_EC_Trixi for further details.
@muladd begin

abstract type AbstractVolumeIntegral end

get_element_variables!(element_variables, u, mesh, equations,
                       volume_integral::AbstractVolumeIntegral, dg, cache) = nothing

"""
    VolumeIntegralStrongForm

The classical strong form volume integral type for FD/DG methods.
"""
struct VolumeIntegralStrongForm <: AbstractVolumeIntegral end

"""
    VolumeIntegralWeakForm

The classical weak form volume integral type for DG methods as explained in standard
textbooks such as
- Kopriva (2009)
  Implementing Spectral Methods for Partial Differential Equations:
  Algorithms for Scientists and Engineers
  [doi: 10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)
"""
struct VolumeIntegralWeakForm <: AbstractVolumeIntegral end

create_cache(mesh, equations, ::VolumeIntegralWeakForm, dg, uEltype) = NamedTuple()

"""
    VolumeIntegralFluxDifferencing

Volume integral type for DG methods based on SBP operators and flux differencing using
symmetric two-point volume fluxes. Based upon the theory developed by
- LeFloch, Mercier, Rohde (2002)
  Fully Discrete, Entropy Conservative Schemes of Arbitrary Order
  [doi: 10.1137/S003614290240069X](https://doi.org/10.1137/S003614290240069X)
- Fisher, Carpenter (2013)
  High-order entropy stable finite difference schemes for nonlinear
  conservation laws: Finite domains
  [doi: 10.1016/j.jcp.2013.06.014](https://doi.org/10.1016/j.jcp.2013.06.014)
- Hendrik Ranocha (2017)
  Comparison of Some Entropy Conservative Numerical Fluxes for the Euler Equations
  [arXiv: 1701.02264](https://arxiv.org/abs/1701.02264)
  [doi: 10.1007/s10915-017-0618-1](https://doi.org/10.1007/s10915-017-0618-1)
- Chen, Shu (2017)
  Entropy stable high order discontinuous Galerkin methods with suitable
  quadrature rules for hyperbolic conservation laws
  [doi: 10.1016/j.jcp.2017.05.025](https://doi.org/10.1016/j.jcp.2017.05.025)
"""
struct VolumeIntegralFluxDifferencing{VolumeFlux} <: AbstractVolumeIntegral
  volume_flux::VolumeFlux
end

function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralFluxDifferencing)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    setup = [
            "volume flux" => integral.volume_flux
            ]
    summary_box(io, "VolumeIntegralFluxDifferencing", setup)
  end
end

"""
    VolumeIntegralShockCapturingHG

Shock-capturing volume integral type for DG methods proposed by
- Hennemann, Gassner (2020)
  "A provably entropy stable subcell shock capturing approach for high order split form DG"
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)
"""
struct VolumeIntegralShockCapturingHG{VolumeFluxDG, VolumeFluxFV, Indicator} <: AbstractVolumeIntegral
  volume_flux_dg::VolumeFluxDG # symmetric, e.g. split-form or entropy-conservative
  volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
  indicator::Indicator
end

function VolumeIntegralShockCapturingHG(indicator; volume_flux_dg=flux_central,
                                                   volume_flux_fv=flux_lax_friedrichs)
  VolumeIntegralShockCapturingHG{typeof(volume_flux_dg), typeof(volume_flux_fv), typeof(indicator)}(
    volume_flux_dg, volume_flux_fv, indicator)
end

function Base.show(io::IO, mime::MIME"text/plain", integral::VolumeIntegralShockCapturingHG)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    summary_header(io, "VolumeIntegralShockCapturingHG")
    summary_line(io, "volume flux DG", integral.volume_flux_dg)
    summary_line(io, "volume flux FV", integral.volume_flux_fv)
    summary_line(io, "indicator", integral.indicator |> typeof |> nameof)
    show(increment_indent(io), mime, integral.indicator)
    summary_footer(io)
  end
end

"""
    VolumeIntegralPureLGLFiniteVolume

A volume integral that only uses the subcell finite volume scheme from the paper
- Hennemann, Gassner (2020)
  "A provably entropy stable subcell shock capturing approach for high order split form DG"
  [arXiv: 2008.12044](https://arxiv.org/abs/2008.12044)

This gives a formally O(1)-accurate finite volume scheme on an LGL-type subcell
mesh (LGL = Legendre-Gauss-Lobatto).

!!! warning "Experimental implementation"
    This is an experimental feature and may change in future releases.
"""
struct VolumeIntegralPureLGLFiniteVolume{VolumeFluxFV} <: AbstractVolumeIntegral
  volume_flux_fv::VolumeFluxFV # non-symmetric in general, e.g. entropy-dissipative
end
# TODO: Figure out if this can also be used for Gauss nodes, not just LGL, and adjust the name accordingly

function Base.show(io::IO, ::MIME"text/plain", integral::VolumeIntegralPureLGLFiniteVolume)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    setup = [
            "FV flux" => integral.volume_flux_fv
            ]
    summary_box(io, "VolumeIntegralPureLGLFiniteVolume", setup)
  end
end

function get_element_variables!(element_variables, u, mesh, equations,
                                volume_integral::VolumeIntegralShockCapturingHG, dg, cache)
  # call the indicator to get up-to-date values for IO
  volume_integral.indicator(u, mesh, equations, dg, cache)
  get_element_variables!(element_variables, volume_integral.indicator, volume_integral)
end

abstract type AbstractSurfaceIntegral end

"""
    SurfaceIntegralWeakForm(surface_flux=flux_central)

The classical weak form surface integral type for DG methods as explained in standard
textbooks such as
- Kopriva (2009)
  Implementing Spectral Methods for Partial Differential Equations:
  Algorithms for Scientists and Engineers
  [doi: 10.1007/978-90-481-2261-5](https://doi.org/10.1007/978-90-481-2261-5)

See also [`VolumeIntegralWeakForm`](@ref).
"""
struct SurfaceIntegralWeakForm{SurfaceFlux} <: AbstractSurfaceIntegral
  surface_flux::SurfaceFlux
end

SurfaceIntegralWeakForm() = SurfaceIntegralWeakForm(flux_central)

function Base.show(io::IO, ::MIME"text/plain", integral::SurfaceIntegralWeakForm)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    setup = [
            "surface flux" => integral.surface_flux
            ]
    summary_box(io, "SurfaceIntegralWeakForm", setup)
  end
end

"""
    SurfaceIntegralStrongForm(surface_flux=flux_central)

The classical strong form surface integral type for FD/DG methods.

See also [`VolumeIntegralStrongForm`](@ref).
"""
struct SurfaceIntegralStrongForm{SurfaceFlux} <: AbstractSurfaceIntegral
  surface_flux::SurfaceFlux
end

SurfaceIntegralStrongForm() = SurfaceIntegralStrongForm(flux_central)

function Base.show(io::IO, ::MIME"text/plain", integral::SurfaceIntegralStrongForm)
  @nospecialize integral # reduce precompilation time

  if get(io, :compact, false)
    show(io, integral)
  else
    setup = [
            "surface flux" => integral.surface_flux
            ]
    summary_box(io, "SurfaceIntegralStrongForm", setup)
  end
end

"""
    DG(; basis, mortar, surface_integral, volume_integral)

Create a discontinuous Galerkin method.
If [`basis isa LobattoLegendreBasis`](@ref LobattoLegendreBasis),
this creates a [`DGSEM`](@ref).
"""
struct DG{Basis, Mortar, SurfaceIntegral, VolumeIntegral}
  basis::Basis
  mortar::Mortar
  surface_integral::SurfaceIntegral
  volume_integral::VolumeIntegral
end

function Base.show(io::IO, dg::DG)
  @nospecialize dg # reduce precompilation time

  print(io, "DG{", real(dg), "}(")
  print(io,       dg.basis)
  print(io, ", ", dg.mortar)
  print(io, ", ", dg.surface_integral)
  print(io, ", ", dg.volume_integral)
  print(io, ")")
end

function Base.show(io::IO, mime::MIME"text/plain", dg::DG)
  @nospecialize dg # reduce precompilation time

  if get(io, :compact, false)
    show(io, dg)
  else
    summary_header(io, "DG{" * string(real(dg)) * "}")
    summary_line(io, "basis", dg.basis)
    summary_line(io, "mortar", dg.mortar)
    summary_line(io, "surface integral", dg.surface_integral |> typeof |> nameof)
    show(increment_indent(io), mime, dg.surface_integral)
    summary_line(io, "volume integral", dg.volume_integral |> typeof |> nameof)
    if !(dg.volume_integral isa VolumeIntegralWeakForm)
      show(increment_indent(io), mime, dg.volume_integral)
    end
    summary_footer(io)
  end
end

Base.summary(io::IO, dg::DG) = print(io, "DG(" * summary(dg.basis) * ")")

@inline Base.real(dg::DG) = real(dg.basis)

function get_element_variables!(element_variables, u, mesh, equations, dg::DG, cache)
  get_element_variables!(element_variables, u, mesh, equations, dg.volume_integral, dg, cache)
end

const MeshesDGSEM = Union{TreeMesh, StructuredMesh, UnstructuredMesh2D, P4estMesh}

@inline ndofs(mesh::MeshesDGSEM, dg::DG, cache) = nelements(cache.elements) * nnodes(dg)^ndims(mesh)

# TODO: Taal performance, 1:nnodes(dg) vs. Base.OneTo(nnodes(dg)) vs. SOneTo(nnodes(dg)) for DGSEM
@inline eachnode(dg::DG) = Base.OneTo(nnodes(dg))
@inline nnodes(dg::DG)   = nnodes(dg.basis)

# This is used in some more general analysis code and needs to dispatch on the
# `mesh` for some combinations of mesh/solver.
@inline nelements(mesh, dg::DG, cache) = nelements(dg, cache)

@inline eachelement(dg::DG, cache)   = Base.OneTo(nelements(dg, cache))
@inline eachinterface(dg::DG, cache) = Base.OneTo(ninterfaces(dg, cache))
@inline eachboundary(dg::DG, cache)  = Base.OneTo(nboundaries(dg, cache))
@inline eachmortar(dg::DG, cache)    = Base.OneTo(nmortars(dg, cache))
@inline eachmpiinterface(dg::DG, cache) = Base.OneTo(nmpiinterfaces(dg, cache))

@inline nelements(dg::DG, cache)   = nelements(cache.elements)
@inline nelementsglobal(dg::DG, cache) = mpi_isparallel() ? cache.mpi_cache.n_elements_global : nelements(dg, cache)
@inline ninterfaces(dg::DG, cache) = ninterfaces(cache.interfaces)
@inline nboundaries(dg::DG, cache) = nboundaries(cache.boundaries)
@inline nmortars(dg::DG, cache)    = nmortars(cache.mortars)
@inline nmpiinterfaces(dg::DG, cache) = nmpiinterfaces(cache.mpi_interfaces)

# The following functions assume an array-of-structs memory layout
# We would like to experiment with different mamory layout choices
# in the future, see
# - https://github.com/trixi-framework/Trixi.jl/issues/88
# - https://github.com/trixi-framework/Trixi.jl/issues/87
# - https://github.com/trixi-framework/Trixi.jl/issues/86
@inline function get_node_coords(x, equations, solver::DG, indices...)
  SVector(ntuple(@inline(idx -> x[idx, indices...]), Val(ndims(equations))))
end

@inline function get_node_vars(u, equations, solver::DG, indices...)
  # There is a cut-off at `n == 10` inside of the method
  # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17
  # in Julia `v1.5`, leading to type instabilities if
  # more than ten variables are used. That's why we use
  # `Val(...)` below.
  # We use `@inline` to make sure that the `getindex` calls are
  # really inlined, which might be the default choice of the Julia
  # compiler for standard `Array`s but not necessarily for more
  # advanced array types such as `PtrArray`s, cf.
  # https://github.com/JuliaSIMD/VectorizationBase.jl/issues/55
  SVector(ntuple(@inline(v -> u[v, indices...]), Val(nvariables(equations))))
end

@inline function get_surface_node_vars(u, equations, solver::DG, indices...)
  # There is a cut-off at `n == 10` inside of the method
  # `ntuple(f::F, n::Integer) where F` in Base at ntuple.jl:17
  # in Julia `v1.5`, leading to type instabilities if
  # more than ten variables are used. That's why we use
  # `Val(...)` below.
  u_ll = SVector(ntuple(@inline(v -> u[1, v, indices...]), Val(nvariables(equations))))
  u_rr = SVector(ntuple(@inline(v -> u[2, v, indices...]), Val(nvariables(equations))))
  return u_ll, u_rr
end

@inline function set_node_vars!(u, u_node, equations, solver::DG, indices...)
  for v in eachvariable(equations)
    u[v, indices...] = u_node[v]
  end
  return nothing
end

@inline function add_to_node_vars!(u, u_node, equations, solver::DG, indices...)
  for v in eachvariable(equations)
    u[v, indices...] += u_node[v]
  end
  return nothing
end

# Use this function instead of `add_to_node_vars` to speed up
# multiply-and-add-to-node-vars operations
# See https://github.com/trixi-framework/Trixi.jl/pull/643
@inline function multiply_add_to_node_vars!(u, factor, u_node, equations, solver::DG, indices...)
  for v in eachvariable(equations)
    u[v, indices...] = u[v, indices...] + factor * u_node[v]
  end
  return nothing
end

# Used for analyze_solution
SolutionAnalyzer(dg::DG; kwargs...) = SolutionAnalyzer(dg.basis; kwargs...)

AdaptorAMR(mesh, dg::DG) = AdaptorL2(dg.basis)

# General structs for discretizations based on the basic principle of
# DGSEM (discontinuous Galerkin spectral element method)
include("dgsem/dgsem.jl")

function allocate_coefficients(mesh::AbstractMesh, equations, dg::DG, cache)
  # We must allocate a `Vector` in order to be able to `resize!` it (AMR).
  # cf. wrap_array
  zeros(eltype(cache.elements), nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache))
end

@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DGSEM, cache)
  @boundscheck begin
    @assert length(u_ode) == nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)
  end
  # We would like to use
  #     reshape(u_ode, (nvariables(equations), ntuple(_ -> nnodes(dg), ndims(mesh))..., nelements(dg, cache)))
  # but that results in
  #     ERROR: LoadError: cannot resize array with shared data
  # when we resize! `u_ode` during AMR.
  #
  # !!! danger "Segfaults"
  #     Remember to `GC.@preserve` temporaries such as copies of `u_ode`
  #     and other stuff that is only used indirectly via `wrap_array` afterwards!

  # Currently, there are problems when AD is used with `PtrArray`s in broadcasts
  # since LoopVectorization does not support `ForwardDiff.Dual`s. Hence, we use
  # optimized `PtrArray`s whenever possible and fall back to plain `Array`s
  # otherwise.
  if LoopVectorization.check_args(u_ode)
    # This version using `PtrArray`s from StrideArrays.jl is very fast and
    # does not result in allocations.
    #
    # !!! danger "Heisenbug"
    #     Do not use this code when `@threaded` uses `Threads.@threads`. There is
    #     a very strange Heisenbug that makes some parts very slow *sometimes*.
    #     In fact, everything can be fast and fine for many cases but some parts
    #     of the RHS evaluation can take *exactly* (!) five seconds randomly...
    #     Hence, this version should only be used when `@threaded` is based on
    #     `@batch` from Polyester.jl or something similar. Using Polyester.jl
    #     is probably the best option since everything will be handed over to
    #     Chris Elrod, one of the best performance software engineers for Julia.
    PtrArray(pointer(u_ode),
             (StaticInt(nvariables(equations)), ntuple(_ -> StaticInt(nnodes(dg)), ndims(mesh))..., nelements(dg, cache)))
            #  (nvariables(equations), ntuple(_ -> nnodes(dg), ndims(mesh))..., nelements(dg, cache)))
  else
    # The following version is reasonably fast and allows us to `resize!(u_ode, ...)`.
    unsafe_wrap(Array{eltype(u_ode), ndims(mesh)+2}, pointer(u_ode),
                (nvariables(equations), ntuple(_ -> nnodes(dg), ndims(mesh))..., nelements(dg, cache)))
  end
end

# General fallback
@inline function wrap_array(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DG, cache)
  wrap_array_native(u_ode, mesh, equations, dg, cache)
end

# Like `wrap_array`, but guarantees to return a plain `Array`, which can be better
# for interfacing with external C libraries (MPI, HDF5, visualization),
# writing solution files etc.
@inline function wrap_array_native(u_ode::AbstractVector, mesh::AbstractMesh, equations, dg::DG, cache)
  @boundscheck begin
    @assert length(u_ode) == nvariables(equations) * nnodes(dg)^ndims(mesh) * nelements(dg, cache)
  end
  unsafe_wrap(Array{eltype(u_ode), ndims(mesh)+2}, pointer(u_ode),
              (nvariables(equations), ntuple(_ -> nnodes(dg), ndims(mesh))..., nelements(dg, cache)))
end

function compute_coefficients!(u, func, t, mesh::AbstractMesh{1}, equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for i in eachnode(dg)
      x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i, element)
      u_node = func(x_node, t, equations)
      set_node_vars!(u, u_node, equations, dg, i, element)
    end
  end
end

function compute_coefficients!(u, func, t, mesh::AbstractMesh{2}, equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for j in eachnode(dg), i in eachnode(dg)
      x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, element)
      u_node = func(x_node, t, equations)
      set_node_vars!(u, u_node, equations, dg, i, j, element)
    end
  end
end

function compute_coefficients!(u, func, t, mesh::AbstractMesh{3}, equations, dg::DG, cache)

  @threaded for element in eachelement(dg, cache)
    for k in eachnode(dg), j in eachnode(dg), i in eachnode(dg)
      x_node = get_node_coords(cache.elements.node_coordinates, equations, dg, i, j, k, element)
      u_node = func(x_node, t, equations)
      set_node_vars!(u, u_node, equations, dg, i, j, k, element)
    end
  end
end

# Discretizations specific to each mesh type of Trixi
# If some functionality is shared by multiple combinations of meshes/solvers,
# it is defined in the directory of the most basic mesh and solver type.
# The most basic solver type in Trixi is DGSEM (historic reasons and background
# of the main contributors).
# We consider the `TreeMesh` to be the most basic mesh type since it is Cartesian
# and was the first mesh in Trixi. The order of the other mesh types is the same
# as the include order below.
include("dgsem_tree/dg.jl")
include("dgsem_structured/dg.jl")
include("dgsem_unstructured/dg.jl")
include("dgsem_p4est/dg.jl")

# Finite difference methods using summation by parts (SBP) operators
# These methods are very similar to DG methods since they also impose interface
# and boundary conditions weakly. Thus, these methods can re-use a lot of
# functionality implemented for DGSEM.
include("fdsbp_tree/fdsbp_2d.jl")

end # @muladd


text:
# * Diagonal matrices

function diagmuladd!(α, A::Diagonal, B::AbstractVector, β, C::AbstractVector)
    mA, nA = size(A)
    mB = length(B)
    nA == mB || throw(DimensionMismatch("Dimensions must match"))
    length(C) == mA || throw(DimensionMismatch("Dimensions must match"))

    isone(β) || lmul!(β, C)
    (nA == 0 || mB == 0) && return C

    @inbounds for k = 1:mB
        C[k] += α * A.diag[k] * B[k]
    end

    C
end

function LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T,
                                            <:RadialOperator{<:Any, <:Any, <:Diagonal},
                                            Source, Dest}) where {T,Source<:RadialOrbital,Dest<:RadialOrbital}
    diagmuladd!(ma.α, ma.A.args[2], ma.B.args[2], ma.β, ma.C.args[2])
    ma.C
end

LazyArrays.default_blasmul!(α, A::Diagonal, B::AbstractVector, β, C::AbstractVector) =
    diagmuladd!(α, A, B, β, C)

# * SymTridiagonal matrices

function symtridiagmuladd!(α, A::SymTridiagonal, B::AbstractVector, β, C::AbstractVector)
    mA, nA = size(A)
    m = length(B)
    mA == nA == m || throw(DimensionMismatch("Dimensions must match"))
    length(C) == m || throw(DimensionMismatch("Dimensions must match"))

    isone(β) || lmul!(β, C)
    (nA == 0 || m == 0) && return C

    d = A.dv
    e = A.ev

    @inbounds begin
        x₊ = B[1]
        x₀ = zero(x₊)
        # If m == 1 then e[1] is out of bounds
        e₀ = m > 1 ? zero(e[1]) : zero(eltype(e))
        for i = 1:m - 1
            x₋, x₀, x₊ = x₀, x₊, B[i + 1]
            e₋, e₀ = e₀, e[i]
            C[i] += α * (e₋*x₋ + d[i]*x₀ + e₀*x₊)
        end
        C[m] += α * (e₀*x₀ + d[m]*x₊)
    end

    C
end

function LazyArrays.materialize!(ma::MulAdd{<:Any, <:Any, <:Any, T,
                                            <:RadialOperator{<:Any, <:Any, <:SymTridiagonal},
                                            Source, Dest}) where {T,Source<:RadialOrbital,Dest<:RadialOrbital}
    symtridiagmuladd!(ma.α, ma.A.args[2], ma.B.args[2], ma.β, ma.C.args[2])
    ma.C
end

LazyArrays.default_blasmul!(α, A::SymTridiagonal, B::AbstractVector, β, C::AbstractVector) =
    symtridiagmuladd!(α, A, B, β, C)


text:
using Iterators

# type definitions
immutable ProblemDescription{T<:FloatingPoint}
    timestep::T

    # organized in the order of x-y-z
    boundary::Array{T,1}
    discretesteps::Array{Int16,1}

    tolerence::T

    nstates::Int8
end

type SystemDescription{T<:FloatingPoint}
    cres::Array{T, 1}
    coordindex::Iterators.Product
    pcoordinategrid::Array{Array{T,1}, 1}
    mcoordinategrid::Array{Array{T,1}, 1}

    fftplan!::Function
    ifftplan!::Function

    normenergies::Array{Complex{T},1}
    energies::Array{Complex{T}, 1}

    fftnorm::Complex{T}
end

type StateZero
end

type State{T<:Complex, N}
    components::Array{T, N}
    inpositionbasis::Bool
    frozen::Bool
end

typealias State1D{T<:Complex} State{T, 1}
typealias State2D{T<:Complex} State{T, 2}
typealias State3D{T<:Complex} State{T, 3}

type StateSet{T<:Complex, N}
    states::Array{State{T,N},1}
end

typealias StateSet1D{T<:Complex} StateSet{T, 1}
typealias StateSet2D{T<:Complex} StateSet{T, 2}
typealias StateSet3D{T<:Complex} StateSet{T, 3}

immutable QuantumOperator{T<:Complex, N}
    doperator::Array{T, N}
    positionbasis::Bool
end

immutable CompositeQuantumOperator{T<:Complex, N}
    queue::Array{QuantumOperator{T,N},1}
end

QOperator = Union(QuantumOperator, CompositeQuantumOperator)
### Outer constructors
State{T<:Complex}(::Type{T},vars...) = State(rand(T, vars...), true, false)

StateSet{T<:Complex}(::Type{T}, number, vars...) =
                                    StateSet([State(T,vars...) for i=1:number])
StateSet{T<:Complex,N}(x::Array{Array{T,N},1}) = StateSet([State(i, true, false) for i in x])

### Method definitions
Base.call{T<:Complex}(::Type{State1D}, ::Type{T}, a) = State{T,1}(a)
Base.call{T<:Complex}(::Type{State2D}, ::Type{T}, a, b) = State{T,2}(a, b)
Base.call{T<:Complex}(::Type{State3D}, ::Type{T}, a, b, c) = State{T,3}(a, b, c)
Base.call{T<:Complex}(::Type{StateSet1D}, ::Type{T}, n, a) =
                                                StateSet{T,1}(n, a)
Base.call{T<:Complex}(::Type{StateSet2D}, ::Type{T}, n, a, b) =
                                                StateSet{T,2}(n, a, b)
Base.call{T<:Complex}(::Type{StateSet3D}, ::Type{T}, n, a, b, c) =
                                                StateSet{T,3}(n, a, b, c)

Base.call{T,N}(x::State{T,N}, y::Array{T,N}) = begin
    if(size(x)==size(y))
        x.components=y
    else
        error("Cannot change the dimensions of a state.")
    end
end
Base.call{T,N}(x::StateSet{T,N}, a::Array{Array{T,N},1}) = begin
    for i=1:length(x)
        x.states[i](a[i])
    end
end
Base.call{T,N}(x::StateSet{T,N}, a::Array{State{T,N},1}) = (x.states = a)

Base.getindex(x::State, vars...) = x.components[vars...]
Base.getindex(x::StateSet, n) = x.states[n]
Base.getindex(x::StateSet, n, vars...) = x.states[n].components[vars...]
Base.getindex(x::QuantumOperator, vars...) = x.doperator[vars...]

Base.size{T<:State}(x::T) = size(x.components)
Base.size{T<:StateSet}(x::T) = [length(x.states), size(x.states[1])...]

Base.length{T<:State}(x::T) = length(x.components)
Base.length{T<:StateSet}(x::T) = length(x.states)

Base.convert{T,N}(::Type{Array{T,N}}, x::State{T,N}) = x.components
Base.convert{T,N}(::Type{Array{State{T,N},1}}, x::StateSet{T,N}) = x.states
Base.convert{T,N}(::Type{Array{Array{T,N},1}}, x::StateSet{T,N}) =
                                            [i.components for i in x.states]
Base.convert{T,N}(::Type{StateSet{T,N}}, x::Array{Any, 1}) = begin
    a::Array{State{T,N}, 1} = [i for i in x]
    StateSet(a)
end

Base.start{T<:StateSet}(x::T) = 1
Base.done{T<:StateSet}(x::T, s) = ((length(x)+1)==s)
Base.next{T<:StateSet}(x::T, s) = x[s], s+1

Base.zero(::Type(Any)) = StateZero()
### Function definitions

function generatesystemdescriptor{T<:FloatingPoint}(pdesc::ProblemDescription{T})
    cres = zeros(T, length(pdesc.boundary))
    cgrid = product([[i for i=1:pdesc.discretesteps[j]] for j=1:length(pdesc.discretesteps)])
    for i=1:length(pdesc.boundary)
        cres[i] = 2*pdesc.boundary[i]/pdesc.discretesteps[i]
    end
    csys = [[-pdesc.boundary[j]+(i-1)*cres[j] for i=1:(pdesc.discretesteps[j])] for j=1:length(pdesc.discretesteps)]
    kcoord = [[(i<=(pdesc.discretesteps[j]/2)?(i-1):(i-pdesc.discretesteps[j]-1))*π/pdesc.boundary[j] for i=1:pdesc.discretesteps[j]] for j=1:length(pdesc.discretesteps)]

    fplan = plan_fft!(rand(Complex{T},pdesc.discretesteps...))
    ifplan = plan_bfft!(rand(Complex{T},pdesc.discretesteps...))
    nenergies = zeros(Complex{T},pdesc.nstates)
    energies = zeros(Complex{T},pdesc.nstates)

    csys1 = convert(Array{Array{T,1},1},[convert(Array{T,1},[p...]) for p in product(csys...)])
    kcoord1 = convert(Array{Array{T,1},1},[convert(Array{T,1},[p...]) for p in product(kcoord...)])

    SystemDescription(cres, cgrid, csys1, kcoord1,fplan, ifplan, nenergies, energies, convert(Complex{T},prod(pdesc.discretesteps)^0.5))
end

### too expensive to assert the same dimensions every single time?
function assertsamedimensions{T<:State}(x::T, y::T)
    size(x) == size(y) ? size(x): error("Incompatible States")
end

function discretizeOperators{T<:FloatingPoint}(operator::Function,
                    syd::SystemDescription{T}, pdesc::ProblemDescription{T},
                                positionbasis::Bool=true)
    if positionbasis==true
        grid = syd.pcoordinategrid
    else
        grid = syd.mcoordinategrid
    end

    t = [convert(Complex{T}, operator(p...)) for p in grid]
    QuantumOperator(convert(Array{Complex{T}},reshape(t, convert(Array{Int64,1},pdesc.discretesteps)...)), positionbasis)
end

(+){T<:State}(a::StateZero, b::T) = b
(+){T<:State}(a::T, b::T) = begin
    if a.inpositionbasis==b.inpositionbasis
        c = deepcopy(a)
        c.components = a.components + b.components
        c
    else
        error("States represented in two different bases cannot be added.")
    end
end
(+){T<:StateSet}(a::T, b::T) = begin
    convert(T, map(+, a, b))
end

(+){T<:QuantumOperator}(a::T, b::T) = begin
    if a.positionbasis == b.positionbasis
        QuantumOperator(a.doperator+b.doperator, a.positionbasis)
    else
        error("Two operators are expressed in different basis.")
    end
end

(-){T<:State}(a::T, b::T) = a + (-1*b)
(-){T<:StateSet}(a::T, b::T) = a + (-1*b)
(-){T<:QuantumOperator}(a::T, b::T) = a + (-1*b)

(*){T<:Complex,Z<:Number,N}(a::Z, b::State{T,N}) = begin
    c = deepcopy(b)
    c.components = convert(T,a) * b.components
    c
end
(*){T<:Complex,Z<:Number,N}(a::Z, b::StateSet{T,N}) = begin
    c::Array{State{T,N},1} = [a*i for i in b]
    StateSet(c)
end
(*){T<:Complex,N,Z<:Number}(a::Array{Z,2}, b::StateSet{T,N}) = begin
    convert(StateSet{T,N}, a*b.states)
end

(*){T<:Complex,N}(O::QuantumOperator{T, N}, s::State{T,N}) = begin
    s1 = deepcopy(s)
    multiply!(O, s1)
    s1
end

(*){T<:Complex,N}(O::QuantumOperator{T,N}, s::StateSet{T,N}) = begin
    s1 = deepcopy(s)
    multiply!(O, s1)
    s1
end

(*){T<:State}(s1::T, s2::T) = innerproduct(s1, s2)

(*){T<:QuantumOperator,Z<:Number}(a::Z, O::T) = begin
    c = a*O.doperator
    QuantumOperator(c, O.positionbasis)
end

(*){T<:QuantumOperator}(a::T, b::T) = begin
    CompositeQuantumOperator([b,a])
end

(*){T<:Complex,N}(a::QuantumOperator{T,N},b::CompositeQuantumOperator{T,N}) = begin
    c = deepcopy(b)
    push!(c.queue, a)
    c
end

(*){T<:Complex,N}(a::CompositeQuantumOperator{T,N}, b::QuantumOperator{T,N}) = begin
    c = deepcopy(a)
    insert!(c.queue, 1, b)
    c
end

(*){T<:Complex,N}(a::CompositeQuantumOperator{T,N}, b::State{T,N}) = begin
    c = deepcopy(b)
    o = deepcopy(a)
    while !isempty(o.queue)
        c  = pop!(o.queue)*c
    end
    c
end

(/){T<:Complex,N,Z<:Number}(x::State{T,N}, a::Z) = begin
    c = deepcopy(x)
    c.components = x.components / convert(T,a)
    c
end

(/){T<:QuantumOperator,Z<:Number}(a::Z, O::T) = begin
    c = O.doperator/a
    QuantumOperator(c, O.positionbasis)
end

Base.exp{T<:QuantumOperator}(O::T) = begin
    c = exp(O.doperator)
    QuantumOperator(c, O.positionbasis)
end

function multiply!{T<:Complex,N}(O::QuantumOperator{T,N}, b::State{T,N})
    if O.positionbasis==b.inpositionbasis
        for ind in eachindex(b.components)
            b.components[ind...] = O[ind...]*b[ind...]
        end
    else
        error("Operator represented in a  different basis.")
    end
end

function multiply!{T<:Complex,N}(O::CompositeQuantumOperator{T,N}, b::State{T,N})
    Op = deepcopy(O)
    while !isempty(Op.queue)
        multiply!(pop!(Op.queue), b)
    end
    b
end

function multiply!{T<:Complex,N}(O::QuantumOperator{T,N}, b::StateSet{T,N})
    for c in b
        multiply!(O, c)
    end
end

function multiply!{T<:Complex,N}(O::CompositeQuantumOperator{T,N}, b::StateSet{T,N})
    for c in b
        multiply!(O, c)
    end
end

function multiply!{T<:FloatingPoint,N}(O::QuantumOperator{Complex{T},N}, s::State{Complex{T},N}, syd::SystemDescription{T})
    if O.positionbasis==s.inpositionbasis
        multiply!(O, s)
    else
        if O.positionbasis
            positionspace!(s, syd)
            for ind in eachindex(s.components)
                s.components[ind...] = O[ind...]*s[ind...]
            end
            momentumspace!(s, syd)
        else
            momentumspace!(s, syd)
            for ind in eachindex(s.components)
                s.components[ind...] = O[ind...]*s[ind...]
            end
            positionspace!(s, syd)
        end
    end
end

function multiply!{T<:FloatingPoint,N}(O::CompositeQuantumOperator{Complex{T},N}, s::State{Complex{T},N}, syd::SystemDescription{T})
    Op = deepcopy(O)
    while !isempty(Op.queue)
        multiply!(pop!(Op.queue), s, syd)
    end
end

function multiply!{T<:FloatingPoint,N}(O::QuantumOperator{Complex{T},N}, s::StateSet{Complex{T},N}, syd::SystemDescription{T})
    for c in s
        multiply!(O,c,syd)
    end
end

function multiply!{T<:FloatingPoint,N}(O::CompositeQuantumOperator{Complex{T},N}, s::StateSet{Complex{T},N}, syd::SystemDescription{T})
    for c in s
        multiply!(O,c,syd)
    end
end

function multiply!{T<:Complex,Z<:Number,N}(a::Z, b::State{T,N})
    b.components = convert(T,a)*b.components
end

function multiply!{T<:Complex,N,Z<:Number}(a::Array{Z,2}, b::StateSet{T,N})
    b.states = a*b.states
end

function add!{T<:State}(a::T, b::T)
    a.components = a.components + b.components
end

function divide!{T<:Complex,Z<:Number,N}(a::State{T,N}, b::Z)
    a.components = a.components / convert(T, b)
end

# Does not check if they are of the same basis
function innerproduct{T<:Complex,N}(x::State{T,N}, y::State{T,N})
    s = assertsamedimensions(x,y)
    n = length(x.components)
    k = zero(T)
    for i=1:n
        indices = ind2sub(s, i)
        k += conj(x[indices...])*y[indices...]
    end
    k
end

function innerproduct{T<:Complex,N}(x::StateSet{T,N}, y::StateSet{T,N})
    convert(Array{T,1}, map(innerproduct, x, y))
end

LinAlg.norm{T<:State}(x::T) = sqrt(innerproduct(x,x))
LinAlg.norm{T<:StateSet}(x::T) = [norm(i) for i in x]

function overlapmatrix{T<:Complex,N}(x::StateSet{T,N})
    temp = zeros(T, length(x), length(x))
    for i=1:length(x), j=1:length(x)
        temp[i,j] = innerproduct(x[i], x[j])
    end
    Hermitian(temp)
end

function normalize!{T,N}(x::State{T,N})
    divide!(x,norm(x))
end

function normalize!{T,N}(x::StateSet{T,N})
    for i in x
        normalize!(i)
    end
end

function orthonormalize!{T,N}(s::StateSet{T,N})
    overlap = overlapmatrix(s)
    eigfactors = eigfact(overlap)
    multiply!(eigfactors[:vectors].', s)
    for i=1:length(s)
        divide!(s[i], sqrt(complex(eigfactors[:values][i])))
    end
    eigfactors[:values]
end

function orthonormalize{T,N}(s::StateSet{T,N})
    overlap = overlapmatrix(s)
    eigfactors = eigfact(overlap)
    l = eigfactors[:vectors].' * s
    for i=1:length(l)
        divide!(l[i], sqrt(complex(eigfactors[:values][i])))
    end
    l, eigfactors[:values]
end

function positionspace!{T<:FloatingPoint,N}(s::State{Complex{T},N}, syd::SystemDescription{T})
    if s.inpositionbasis==true
        return
    end
    syd.ifftplan!(s.components)
    divide!(s,syd.fftnorm)
    s.inpositionbasis = true
end

function positionspace!{T<:FloatingPoint,N}(s::StateSet{Complex{T}, N}, syd::SystemDescription{T})
    for st in s
        positionspace!(st, syd)
    end
end

function positionspace!{T<:FloatingPoint,N}(s::State{Complex{T},N})
    if s.inpositionbasis==true
        return
    end
    bfft!(s.components)
    divide!(s,sqrt(prod(size(s))))
    s.inpositionbasis = true
end

function momentumspace!{T<:FloatingPoint,N}(s::State{Complex{T},N} , syd::SystemDescription{T})
    if s.inpositionbasis==false
        return
    end

    syd.fftplan!(s.components)
    divide!(s,syd.fftnorm)
    s.inpositionbasis = false
end

function momentumspace!{T<:FloatingPoint,N}(s::State{Complex{T},N})
    if s.inpositionbasis==false
        return
    end

    fft!(s.components)
    divide!(s, sqrt(prod(size(s))))
    s.inpositionbasis = false
end

function momentumspace!{T<:FloatingPoint,N}(s::StateSet{Complex{T},N}, syd::SystemDescription{T})
    for st in s
        momentumspace!(st, syd)
    end
end

function expectationvalue{T<:Complex,N}(s::State{T,N}, o::QuantumOperator{T})
    innerproduct(s, o*s)
end

function expectationvalue{T<:Complex,N}(s::StateSet{T,N}, o::QuantumOperator{T})
    innerproduct(s, o*s)
end


text:
function datasets(package_name::AbstractString)
	d = datasets()
	d[findfirst(isequal(package_name), d[:Package]), :]
end

function datasets()
    if RDatasets.__datasets === nothing
        path = joinpath(dirname(@__FILE__), "..", "doc", "datasets.csv")
        global __datasets = CSV.read(path)
    end
    return RDatasets.__datasets::DataFrame
end

max_stars_repo_name: SCXFPROTOCOL/GameRewards-System

max_issues_repo_name: SCXFPROTOCOL/GameRewards-System

max_forks_repo_name: SCXFPROTOCOL/GameRewards-System


text:
#####
# Example evaluation script for TintiNet.jl
# Model IGBT (InceptGOR-Bert): 3-state SS Classifier and PHI, PSI and ACC regressor.
#####

#####
# 0. Loading Packages and defining Hyperparameters
#####

using CUDA, Flux, Transformers, JSON, BSON, Random, Dates, Printf
using IterTools: ncycle;
using Transformers.Basic, Transformers.Pretrain, Transformers.Datasets, Transformers.Stacks, Transformers.BidirectionalEncoder
using BSON: @load, @save
using Flux: @functor, onehot, onecold, mse, shuffle, sigmoid, gradient, unsqueeze
using Flux.Optimise: update!
using Statistics: mean, std

using TintiNet
Random.seed!(102528)

hyperpar = Dict(
    "max_seq_len"       => 128,     # Maximum sequence length on the model (default 192)
    "batch_size"        => 128,     # Batch size for training (default 128)
)

inputfile_path = ARGS[1] ## First argument is the location of the sequence input file
output_path = ARGS[2] ## Second argument is the directory to write the outputs
#print_header_prediction()

#####
# 1. Obtaining the data
#####

headers, sequences = read_sequences_from_file(inputfile_path, hyperpar["max_seq_len"]; filetype="fasta")

#####
# 2. Loading the saved models
#####

@load "./example_classifier_model.bson" checkpoint_model
classifier_model = checkpoint_model |> gpu
Flux.testmode!(classifier_model)

@load "./example_regressor_model.bson" checkpoint_model
regressor_model = checkpoint_model |> gpu
Flux.testmode!(regressor_model)

#####
# 3. Model evaluation
#####

ss_predictions = compute_classifier_predictions(
    classifier_model,
    sequences;
    batched=true,
    eval_batch_size=64,
    use_gpu=true,
    sleep_time_seconds=0.01)

(phi_predictions, psi_predictions, acc_predictions) = compute_regressor_predictions(
    regressor_model,
    sequences;
    batched=false,
    eval_batch_size=64,
    use_gpu=true,
    sleep_time_seconds=0.01)

#####
# 4. Saving Prediction Results
#####

write_csv_predictions(output_path, headers, sequences, ss_predictions, phi_predictions, psi_predictions, acc_predictions)


text:
__precompile__
import FastTransforms
using FastGaussQuadrature

struct Quadrature
    Nq::Int;
    w::Array{Float64,1};
    wTens::Array{Float64,1};
    xi::Array{Float64,1};

   function Quadrature(NqVal,quadType)
        if quadType == "Gauss"
            #xi,w = QuadNodesWeights(NqVal,-1.0,1.0)
            xi,w = gausslegendre(NqVal)
        elseif quadType == "ClenshawCurtis"
            xi,w = QuadNodesWeightsClCu(NqVal,-1.0,1.0)
        end
        xi = xi[end:-1:1];
        w = w[end:-1:1];

        wTens = zeros(NqVal^2);

        for k = 1:NqVal
            for q = 1:NqVal
                wTens[(q-1)*NqVal+k] = w[k]*w[q];
            end
        end

        new(NqVal,w,wTens,xi)
   end

end

function QuadNodesWeightsClCu(Nq::Int64,a::Float64,b::Float64) # transformation to a,b not tested!
    xiQuad,wQuad = FastTransforms.clenshawcurtis(Nq,0.,0.)
    xiQuad=(a*(1-xiQuad)+b*(1+xiQuad))/2.0;
    wQuad = wQuad*(2.0/(b-a));
    return xiQuad,wQuad;
end

function Integral(q::Quadrature, f::Function)
    dot( q.w, f.(q.xi) );
end

function Integral(q::Quadrature, f::Function, a::Float64, b::Float64)
    w = q.w*(b-a)*0.5;
    xi = (a*(ones(Float64,length(w)).-q.xi)+b*(ones(Float64,length(w))+q.xi))/2.0;
    #return q.w'f.(xi);
    return sum(w'f.(xi))
end

function IntegralVec(q::Quadrature, fVec::Array{Float64,1}, a::Float64, b::Float64)
    w = q.w*(b-a)*0.5;
    return sum(w'fVec)
end

function Integral(q::Quadrature, fVec::Array{Float64,1})
    return sum(q.wTens'fVec)
end

function IntegralMat(q::Quadrature, fVec::Array{Float64,2}, a::Float64, b::Float64)
    w = q.w*(b-a)*0.5;
    return fVec*w;
end

function IntegralMat(q::Quadrature, fVec::Array{Float64,2})
    return fVec*q.w;
end


text:
"Perform post-processing to azimuthally-weight RT matrices <<Suniti>>"
function postprocessing_vza!(iμ₀, pol_type, composite_layer, vza, qp_μ, m, vaz, μ₀, weight, nSpec, SFI, R, R_SFI, T, T_SFI)
    # idx of μ0 = cos(sza)
    st_iμ0, istart0, iend0 = get_indices(iμ₀, pol_type);

    # Convert these to Arrays (if CuArrays), so they can be accessed by index
    R⁻⁺ = Array(composite_layer.R⁻⁺);
    T⁺⁺ = Array(composite_layer.T⁺⁺);
    J₀⁺ = Array(composite_layer.J₀⁺);
    J₀⁻ = Array(composite_layer.J₀⁻);

    # Loop over all viewing zenith angles
    for i = 1:length(vza)

        # Find the nearest quadrature point idx
        iμ = nearest_point(qp_μ, cosd(vza[i]));
        st_iμ, istart, iend = get_indices(iμ, pol_type);

        # Compute bigCS
        cos_m_phi, sin_m_phi = (cosd(m * vaz[i]), sind(m * vaz[i]));
        bigCS = weight * Diagonal([cos_m_phi, cos_m_phi, sin_m_phi, sin_m_phi][1:pol_type.n]);

        # Accumulate Fourier moments after azimuthal weighting
        for s = 1:nSpec

            if SFI
                R_SFI[i,:,s] += bigCS * J₀⁻[istart:iend,1, s];
                T_SFI[i,:,s] += bigCS * J₀⁺[istart:iend,1, s];
            else
                R[i,:,s] += bigCS * (R⁻⁺[istart:iend, istart0:iend0, s] / μ₀) * pol_type.I₀;
                T[i,:,s] += bigCS * (T⁺⁺[istart:iend, istart0:iend0, s] / μ₀) * pol_type.I₀;
            end

        end
    end
end


text:
using Distributions;
using Statistics

### Utils
include("utils.jl")

### Ito Processes
include("brownian_motion.jl")
include("brownian_motion_aug.jl")
include("brownian_motion_prescribed.jl")
include("brownian_motion_sobol.jl")
include("geometric_brownian_motion.jl")
include("geometric_brownian_motion_aug.jl")
include("black_scholes.jl")
include("heston.jl")
include("heston_aug.jl")
include("heston_sobol.jl")
include("log_normal_mixture.jl")
include("shifted_log_normal_mixture.jl")

### Finite Activity Levy Processes
include("fa_levy.jl")
include("kou.jl")
include("merton.jl")

### Infinite Activity Levy Processes
include("subordinated_brownian_motion.jl")
include("subordinated_brownian_motion_aug.jl")
include("variance_gamma.jl")
include("normal_inverse_gaussian.jl")

### MultiDimensional
include("nvariate.jl")
include("nvariate_log.jl")


text:
#-----------------------------------------------------------------------------# StatWrapper
abstract type StatWrapper{T} <: OnlineStat{T} end
nobs(o::StatWrapper) = nobs(o.stat)
value(o::StatWrapper) = value(o.stat)
_merge!(a::StatWrapper{T}, b::StatWrapper{T}) where {T} = _merge!(a.stat, b.stat)
name(o::T, args...) where {T<:StatWrapper} = name(typeof(o), args...) * "($(name(o.stat, args...)))"

#-----------------------------------------------------------------------------# CountMissing
"""
    CountMissing(stat)

Calculate a `stat` along with the count of `missing` values.

# Example

    o = CountMissing(Mean())
    fit!(o, [1, missing, 3])
"""
mutable struct CountMissing{T, O<:OnlineStat{T}} <: StatWrapper{Union{Missing,T}}
    stat::O
    nmissing::Int
end
CountMissing(stat::OnlineStat) = CountMissing(stat, 0)
value(o::CountMissing) = (nmissing=o.nmissing, stat=o.stat)
nobs(o::CountMissing) = nobs(o.stat) + o.nmissing

_fit!(o::CountMissing, x) = _fit!(o.stat, x)
_fit!(o::CountMissing, ::Missing) = (o.nmissing += 1)

_merge!(a::CountMissing, b::CountMissing) = (merge!(a.stat, b.stat); a.nmissing += b.nmissing)

#-----------------------------------------------------------------------------# FilterTransform
"""
    FilterTransform(stat::OnlineStat{S}, T = S; filter = x->true, transform = identity)
    FilterTransform(T => filter => transform => stat)

Wrapper around an OnlineStat that the filters and transforms its input.  Note that, depending on
your transformation, you may need to specify the type of a single observation (`T`).

# Examples

    o = FilterTransform(Mean(), Union{Missing,Number}, filter=!ismissing)
    fit!(o, [1, missing, 3])

    o = FilterTransform(String => (x->true) => (x->parse(Int,x)) => Mean())
    fit!(o, "1")
"""
struct FilterTransform{S, T, O<:OnlineStat{T},F,F2} <: StatWrapper{S}
    stat::O
    filter::F
    transform::F2
end
FilterTransform(intype::DataType, stat::OnlineStat; kw...) = FilterTransform(stat, intype; kw...)
function FilterTransform(stat::OnlineStat{T}, intype=T; filter=always_true, transform=identity) where {T}
    FilterTransform{intype, T, typeof(stat), typeof(filter), typeof(transform)}(stat, filter, transform)
end
function FilterTransform(p::Pair{DataType, <:Pair{<:Function, <:Pair{<:Function, <:OnlineStat}}})
    FilterTransform(p[1], p[2][2][2]; filter=p[2][1], transform=p[2][2][1])
end

_fit!(o::FilterTransform, y) = o.filter(y) && _fit!(o.stat, o.transform(y))

additional_info(o::FilterTransform) = (; filter=o.filter, transform=o.transform)

always_true(x) = true

#-----------------------------------------------------------------------------# SkipMissing
"""
    SkipMissing(stat)

Wrapper around an OnlineStat that will skip over `missing` values.

# Example

    o = SkipMissing(Mean())

    fit!(o, [1, missing, 3])
"""
struct SkipMissing{T, O<:OnlineStat{T}} <: StatWrapper{Union{Missing,T}}
    stat::O
    SkipMissing(stat::OnlineStat{T}) where {T} = new{T, typeof(stat)}(stat)
end
_fit!(o::SkipMissing, x::Missing) = nothing
_fit!(o::SkipMissing, x) = _fit!(o.stat, x)
Base.skipmissing(o::OnlineStat) = SkipMissing(o)

#-----------------------------------------------------------------------------# TryCatch
"""
    TryCatch(stat; error_limit=1000, error_message_limit=90)

Wrap each call to `fit!` in a `try`-`catch` block and track the errors encountered (via [`CountMap`](@ref)).  Errors will stop
being tracked after `error_limit` unique errors are encountered.  Only the first `error_message_limit`
characters of each error message will be recorded.

# Example

    o = TryCatch(Mean())

    fit!(o, [1, missing, 3])

    OnlineStatsBase.errors(o)
"""
struct TryCatch{T, O<:OnlineStat{T}} <: StatWrapper{T}
    stat::O
    errors::CountMap{String}
    error_limit::Int
    error_message_limit::Int
end
function TryCatch(stat::OnlineStat; error_limit=1000, error_message_limit=90)
    TryCatch(stat, CountMap(String), error_limit, error_message_limit)
end

errors(o::TryCatch) = value(o.errors)

function additional_info(o::TryCatch)
    ex = errors(o)
    nex = length(ex)
    msg = length(ex) ≥ o.error_limit ? "$nex (limit reached)" : nex
    nex == 0 ? () : (; errors=msg)
end

function handle_error!(o::TryCatch, ex)
    io = IOBuffer()
    Base.showerror(io, ex)
    s = String(take!(io))
    lim = o.error_message_limit
    s = length(s) > lim ? s[1:lim] * "..." : s
    length(value(o.errors)) < o.error_limit && _fit!(o.errors, s)
end

function fit!(o::TryCatch{T}, y::T) where {T}
    try
        _fit!(o.stat, y)
    catch ex
        handle_error!(o, ex)
    end
    o
end

function fit!(o::TryCatch{I}, y::T) where {I, T}
    try
        T == eltype(y) && error("The input for $(name(o,false,false)) is $I.  Found $T.")
        for yi in y
            fit!(o, yi)
        end
    catch ex
        handle_error!(o, ex)
    end
    o
end


text:
"""
    Assignable

`Assignable` represents any data type that can be associated with traits.
It essentially covers all data types including parametric types e.g. `AbstractArray`
"""
const Assignable = Union{UnionAll, DataType}

"""
    Contract{T <: DataType, F <: Function, N}

A contract refers to a function defintion `func` that is required to satisfy
the Can-type of a trait. The function `func` must accepts `args` and returns `ret`.

# Fields
- `can_type`: can-type of a trait e.g. `CanFly`
- `func`: function that must be implemented to satisfy this trait
- `args`: argument types of the function `func`
- `kwargs`: keyword argument names of the function `func`
- `ret`: return type of the function `func`
"""
struct Contract{T <: DataType, F <: Function}
    can_type::T
    func::F
    args::Tuple
    kwargs::Tuple
    ret::Union{DataType,Nothing}
end

function Base.show(io::IO, c::Contract)
    # prefix the can-type with a different symbol because it's not subtype relation
    args = [x == c.can_type ? "🚀$x" : "::$x" for x in c.args]
    args = string("(", join(args, ", "))
    if length(c.kwargs) > 0
        args = string(args, "; ", join(c.kwargs, ", "))
    end
    args = string(args, ")")
    trait = supertype(c.can_type)
    print(io, "$(trait): $(c.can_type) ⇢ $(c.func)$(args)")
    c.ret !== nothing && print(io, "::$(c.ret)")
end

"""
    InterfaceReview

An InterfaceReview object contains the validation results of an interface.

# Fields
- `data_type`: the type being checked
- `result`: true if the type fully implements all required contracts
- `implemented`: an array of implemented contracts
- `misses`: an array of unimplemented contracts
"""
@Base.kwdef struct InterfaceReview
    data_type::Assignable
    result::Bool
    implemented::Vector{Contract}
    misses::Vector{Contract}
end

function Base.show(io::IO, ir::InterfaceReview)
    irtype = ir.data_type
    if length(ir.implemented) == length(ir.misses) == 0
        print(io, "✅ $(irtype) has no interface contract requirements.")
        return
    end
    if length(ir.implemented) > 0
        println(io, "✅ $(irtype) has implemented:")
        for (i, c) in enumerate(ir.implemented)
            println(io, "$(i). $c")
        end
    end
    if length(ir.misses) > 0
        println(io, "❌ $(irtype) is missing these implementations:")
        for (i, c) in enumerate(ir.misses)
            println(io, "$(i). $c")
        end
    end
end

"""
    SyntaxError

Syntax error for macros.
"""
struct SyntaxError <: Exception
    msg
end

const MyDict = VERSION < v"1.2.0" ? Dict : IdDict

"""
    TraitsMap

Map a data type to the Can-type of its assigned traits.
For example, `Dog => Set([CanSwim, CanRun])`.
"""
const TraitsMap = MyDict{Assignable,Set{DataType}}

"""
    InterfaceMap

Map a Can-type to a set of interface contracts.  See [`Contract`](@ref).
"""
const InterfaceMap = MyDict{DataType,Set{Contract}}

"""
    CompositeTraitMap

Maps a composite can-type to a set of its underlying can-types.
e.g. `CanFlySwim => Set([CanFly, CanSwim])`.
"""
const CompositeTraitMap = MyDict{DataType,Set{DataType}}

"""
    TraitsStorage

Keeps all traits-related dynamic data.

# Fields
- `traits_map`: see [`TraitsMap`](@ref)
- `interface_map`: see [`InterfaceMap`](@ref)
- `composite_map`: see [`CompositeTraitMap`](@ref)
"""
struct TraitsStorage
    traits_map::TraitsMap
    interface_map::InterfaceMap
    composite_map::CompositeTraitMap
    TraitsStorage() = new(TraitsMap(), InterfaceMap(), CompositeTraitMap())
end


text:
"""
    resample(rng, x, log_weights, ndraws) -> (draws, psis_result)
    resample(rng, x, ndraws) -> draws

Draw `ndraws` samples from `x`, with replacement.

If `log_weights` is provided, perform Pareto smoothed importance resampling.
"""
function resample(rng, x, log_ratios, ndraws)
    psis_result = PSIS.psis(log_ratios)
    weights = psis_result.weights
    pweights = StatsBase.ProbabilityWeights(weights, one(eltype(weights)))
    return StatsBase.sample(rng, x, pweights, ndraws; replace=true), psis_result
end
resample(rng, x, ndraws) = StatsBase.sample(rng, x, ndraws; replace=true)


text:
# For stochastic Runge-Kutta methods (SIRK and WIRK)
# x - the solution vector to be updated
# V - the matrix containing the drift vector evaluated at the internal stages v(Q_i) (SIRK) or v(Q0_i) (WIRK)
# B - the array containing the diffusion matrix evaluated at the internal stages B(Q_i) (SIRK) or B(Q1^(l)_i) (WIRK)
# bdrift - the Runge-Kutta coefficients for the drift part
# bdiff  - the Runge-Kutta coefficients for the diffusion part
# Δt - the time step
# ΔW - the increments of the Brownian motion (SFIRK) or the increments represented by the random variables \hat I^(k) (WFIRK)
function update_solution!(x::SolutionVector{T}, V::Vector{Vector{T}}, B::Vector{Matrix{T}},
                            bdrift::Vector{T}, bdiff::Vector{T}, Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}

    @assert length(bdrift) == length(bdiff) == length(V) == length(B)

    for i in eachindex(V, B)
        @assert length(x) == length(V[i]) == size(B[i], 1)
        @assert length(ΔW)== size(B[i], 2)
    end

    local Δx::T

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0
        for i in eachindex(V)
            Δx += bdrift[i] * V[i][k]
        end
        x[k] += Δt * Δx
    end

    # Contribution from the diffusion part
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(B)
            for l in eachindex(Δy)
                Δy[l] += bdiff[i] * B[i][k,l]
            end
        end
        x[k] += dot(Δy,ΔW)
    end
end

# For stochastic Runge-Kutta methods
# x - the solution vector to be updated
# V - the matrix containing the drift vector evaluated at the internal stages v(Q_i)
# B - the array containing the diffusion matrix evaluated at the internal stages B(Q_i)
# bdrift - the Runge-Kutta coefficients for the drift part
# bdiff - the Runge-Kutta coefficients for the ΔW terms of the diffusion part
# bdiff2- the Runge-Kutta coefficients for the ΔZ terms of the diffusion part
# Δt - the time step
# ΔW - the increments of the Brownian motion
# ΔZ - the integrals of the increments of the Brownian motion
function update_solution!(x::SolutionVector{T}, V::Vector{Vector{T}}, B::Vector{Matrix{T}},
                            bdrift::Vector{T}, bdiff::Vector{T}, bdiff2::Vector{T},
                            Δt::T, ΔW::Vector{T}, ΔZ::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}

    @assert length(bdrift) == length(bdiff) == length(bdiff2) == length(V) == length(B)

    for i in eachindex(V, B)
        @assert length(x) == length(V[i]) == size(B[i], 1)
        @assert length(ΔW) == length(ΔZ) == size(B[i], 2)
    end

    local Δx::T

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0
        for i in eachindex(V)
            Δx += bdrift[i] * V[i][k]
        end
        x[k] += Δt * Δx
    end

    # Contribution from the diffusion part (ΔW terms)
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(B)
            for l in eachindex(Δy)
                Δy[l] += bdiff[i] * B[i][k,l]
            end
        end
        x[k] += dot(Δy,ΔW)
    end

    # Contribution from the diffusion part (ΔZ terms)
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(B)
            for l in eachindex(Δy)
                Δy[l] += bdiff2[i] * B[i][k,l]
            end
        end
        x[k] += dot(Δy,ΔZ)/Δt
    end
end

# For stochastic partitioned Runge-Kutta methods
# q, p - the solution vector to be updated
# V, F - the matrix containing the drift vectors evaluated at the internal stages v(Q_i), f(Q_i)
# B, G - the array containing the diffusion matrices evaluated at the internal stages B(Q_i), G(Q_i)
# bqdrift, bpdrift - the Runge-Kutta coefficients for the drift parts of the q and p equations
# bqdiff, bpdiff - the Runge-Kutta coefficients for the diffusion parts of the q and p equations
# Δt - the time step
# ΔW - the increments of the Brownian motion
function update_solution!(q::SolutionVector{T}, p::SolutionVector{T},
                          V::Vector{Vector{T}}, F::Vector{Vector{T}},
                          B::Vector{Matrix{T}}, G::Vector{Matrix{T}},
                          bqdrift::Vector{T}, bqdiff::Vector{T},
                          bpdrift::Vector{T}, bpdiff::Vector{T},
                          Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW), Δz::Vector{T}=zero(ΔW)) where {T}

    @assert length(bqdrift) == length(bqdiff) == length(bpdrift) == length(bpdiff) == length(V) == length(F) == length(B) == length(G)

    for i in eachindex(V, F, B, G)
        @assert length(q) == length(p) == length(V[i]) == length(F[i]) == size(B[i], 1) == size(G[i], 1)
        @assert length(ΔW) == size(B[i], 2) == size(G[i], 2)
    end

    local Δq::T
    local Δp::T

    # Contribution from the drift part
    for k in eachindex(q, p)
        Δq = 0
        Δp = 0
        for i in eachindex(V, F)
            Δq += bqdrift[i] * V[i][k]
            Δp += bpdrift[i] * F[i][k]
        end
        q[k] += Δt * Δq
        p[k] += Δt * Δp
    end

    # Contribution from the diffusion part
    for k in eachindex(q, p)
        Δy .= 0
        Δz .= 0

        for i in eachindex(B, G)
            for l in eachindex(Δy,Δz)
                Δy[l] += bqdiff[i] * B[i][k,l]
                Δz[l] += bpdiff[i] * G[i][k,l]
            end
        end

        q[k] += dot(Δy,ΔW)
        p[k] += dot(Δz,ΔW)
    end
end

# For stochastic split partitioned Runge-Kutta methods
# q, p - the solution vector to be updated
# V, F1, F2 - the matrix containing the drift vectors evaluated at the internal stages v(Q_i), fi(Q_i)
# B, G1, G2 - the array containing the diffusion matrices evaluated at the internal stages B(Q_i), Gi(Q_i)
# bqdrift, bpdrift1, bpdrift2 - the Runge-Kutta coefficients for the drift parts of the q and p equations
# bqdiff, bpdiff1, bpdiff2 - the Runge-Kutta coefficients for the diffusion parts of the q and p equations
# Δt - the time step
# ΔW - the increments of the Brownian motion
function update_solution!(q::SolutionVector{T}, p::SolutionVector{T},
                            V::Vector{Vector{T}}, F1::Vector{Vector{T}}, F2::Vector{Vector{T}},
                            B::Vector{Matrix{T}}, G1::Vector{Matrix{T}}, G2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, bqdiff::Vector{T},
                            bpdrift1::Vector{T}, bpdrift2::Vector{T},
                            bpdiff1::Vector{T}, bpdiff2::Vector{T},
                            Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW), Δz::Vector{T}=zero(ΔW)) where {T}

   @assert length(bqdrift) == length(bqdiff) == length(bpdrift1) == length(bpdrift2) == length(bpdiff1) == length(bpdiff2) == length(V) == length(F1) == length(F2) == length(B) == length(G1) == length(G2)

   for i in eachindex(V, F1, F2, B, G1, G2)
       @assert length(q) == length(p) == length(V[i]) == length(F1[i]) == length(F2[i]) == size(B[i], 1) == size(G1[i], 1) == size(G2[i], 1)
       @assert length(ΔW) == size(B[i], 2) == size(G1[i], 2) == size(G2[i], 2)
   end

   local Δq::T
   local Δp::T

   # Contribution from the drift part
   for k in eachindex(q, p)
       Δq = 0
       Δp = 0
       for i in eachindex(V, F1, F2)
           Δq += bqdrift[i] * V[i][k]
           Δp += bpdrift1[i] * F1[i][k] + bpdrift2[i] * F2[i][k]
       end
       q[k] += Δt * Δq
       p[k] += Δt * Δp
   end

   # Contribution from the diffusion part
   for k in eachindex(q, p)
       Δy .= 0
       Δz .= 0

       for i in eachindex(B, G1, G2)
           for l in eachindex(Δy,Δz)
               Δy[l] += bqdiff[i]  * B[i][k,l]
               Δz[l] += bpdiff1[i] * G1[i][k,l] + bpdiff2[i] * G2[i][k,l]
           end
       end

       q[k] += dot(Δy,ΔW)
       p[k] += dot(Δz,ΔW)
   end
end

# For weak Runge-Kutta methods WERK
# x - the solution vector to be updated
# V - the matrix containing the drift vector evaluated at the internal stages v(Q_i)
# B1 - the array containing the diffusion matrix evaluated at the internal stages H^(l)_i, such that B1[:,l,i] is evaluated at H^(l)_i
# B2 - the array containing the diffusion matrix evaluated at the internal stages \hat H^(l)_i, such that B2[:,l,i] is evaluated at \hat H^(l)_i
# α - the Runge-Kutta coefficients for the drift part
# β1 - the Runge-Kutta coefficients for the diffusion term with the random increments
# β2 - the Runge-Kutta coefficients for the second diffusion term
# Δt - the time step
# ΔW - the increments of the Brownian motion represented by the random variables \hat I^(k)
function update_solution!(x::SolutionVector{T}, V::Vector{Vector{T}},
                            B1::Vector{Matrix{T}}, B2::Vector{Matrix{T}}, α::Vector{T},
                            β1::Vector{T}, β2::Vector{T}, Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}
    @assert length(α) == length(β1) == length(β2) == length(V) == length(B1) == length(B2)

    for i in eachindex(V, B1, B2)
        @assert length(x) == length(V[i]) == size(B1[i], 1) == size(B2[i], 1)
        @assert length(ΔW)== size(B1[i], 2) == size(B2[i], 2)
    end

    local Δx::T

    # Contribution from the drift part
    for k in eachindex(x)
        Δx = 0
        for i in eachindex(V)
            Δx += α[i] * V[i][k]
        end
        x[k] += Δt * Δx
    end

    # Contribution from the diffusion term with the random variables I^(k)_i
    for k in eachindex(x)
        Δy .= 0
        for i in eachindex(B1)
            for l in eachindex(Δy)
                Δy[l] += β1[i] * B1[i][k,l]
            end
        end
        x[k] += dot(Δy,ΔW)
    end

    # Contribution from the second diffusion term
    for k in eachindex(x)
        Δx = 0
        for i in eachindex(B2)
            for l in axes(B2[i], 2)
                Δx += β2[i] * B2[i][k,l]
            end
        end
        x[k] += sqrt(Δt)*Δx
    end
end

# For stochastic Runge-Kutta methods
function update_solution!(x::SolutionVector{T}, V::Vector{Vector{T}}, B::Vector{Matrix{T}},
                            bdrift::Vector{T}, b̂drift::Vector, bdiff::Vector{T}, b̂diff::Vector,
                            Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}
    update_solution!(x, V, B, bdrift, bdiff, Δt, ΔW, Δy)
    update_solution!(x, V, B, b̂drift, b̂diff, Δt, ΔW, Δy)
end

# For stochastic partitioned Runge-Kutta methods
function update_solution!(q::SolutionVector{T}, p::SolutionVector{T},
                            V::Vector{Vector{T}}, F::Vector{Vector{T}}, B::Vector{Matrix{T}}, G::Vector{Matrix{T}},
                            bqdrift::Vector{T}, b̂qdrift::Vector, bqdiff::Vector{T}, b̂qdiff::Vector,
                            bpdrift::Vector{T}, b̂pdrift::Vector, bpdiff::Vector{T}, b̂pdiff::Vector,
                            Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW), Δz::Vector{T}=zero(ΔW)) where {T}
    update_solution!(q, p, V, F, B, G, bqdrift, bqdiff, bpdrift, bpdiff, Δt, ΔW, Δy, Δz)
    update_solution!(q, p, V, F, B, G, b̂qdrift, b̂qdiff, b̂pdrift, b̂pdiff, Δt, ΔW, Δy, Δz)
end

# For stochastic split partitioned Runge-Kutta methods
function update_solution!(q::SolutionVector{T}, p::SolutionVector{T},
                            V::Vector{Vector{T}}, F1::Vector{Vector{T}}, F2::Vector{Vector{T}},
                            B::Vector{Matrix{T}}, G1::Vector{Matrix{T}}, G2::Vector{Matrix{T}},
                            bqdrift::Vector{T}, b̂qdrift::Vector, bqdiff::Vector{T}, b̂qdiff::Vector,
                            bpdrift1::Vector{T}, b̂pdrift1::Vector, bpdrift2::Vector{T}, b̂pdrift2::Vector,
                            bpdiff1::Vector{T}, b̂pdiff1::Vector, bpdiff2::Vector{T}, b̂pdiff2::Vector,
                            Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW), Δz::Vector{T}=zero(ΔW)) where {T}
    update_solution!(q, p, V, F1, F2, B, G1, G2, bqdrift, bqdiff, bpdrift1, bpdrift2, bpdiff1, bpdiff2, Δt, ΔW, Δy, Δz)
    update_solution!(q, p, V, F1, F2, B, G1, G2, b̂qdrift, b̂qdiff, b̂pdrift1, b̂pdrift2, b̂pdiff1, b̂pdiff2, Δt, ΔW, Δy, Δz)
end

 # For weak Runge-Kutta methods WERK
 function update_solution!(x::SolutionVector{T}, V::Vector{Vector{T}}, B1::Vector{Matrix{T}}, B2::Vector{Matrix{T}},
                             α::Vector{T}, α̂::Vector, β1::Vector{T}, β̂1::Vector, β2::Vector{T}, β̂2::Vector,
                             Δt::T, ΔW::Vector{T}, Δy::Vector{T}=zero(ΔW)) where {T}
     update_solution!(x, V, B1, B2, α, β1, β2, Δt, ΔW, Δy)
     update_solution!(x, V, B1, B2, α̂, β̂1, β̂2, Δt, ΔW, Δy)
 end


text:
function tenter(d, tg, av)
    @unpack T = d
    av["step"] = "enter"
    pathkey = join(av["path"])
    if av["txt"] == "Pasirinkti"
        for k in keys(T)
            if contains(k, pathkey)
                if !any(av["id"] .== T[k]["dav_id"])
                    push!(T[k]["dav_id"], av["id"])
                end
            end
        end
        if !(av["path"][end]=="0")
            av["path"] = av["path"][1:end-1]
        end
    elseif av["txt"] == "Nuimti žymę"
        for k in keys(T)
            if contains(k,pathkey)
                i = findlast(av["id"] .== T[k]["dav_id"])
                if i !== nothing
                    popat!(T[k]["dav_id"], i)
                end
            end
        end
        if !(av["path"][end]=="0")
            av["path"] = av["path"][1:end-1]
        end
    elseif any(T[pathkey]["children"] .== av["txt"])
        av["path"] = vcat(av["path"], av["txt"])
    elseif !(av["txt"] == "Atgal | Повернутися")
        #@warn "Netinkama įvestis"
    end
    tree(tg, av, T)

    return
end


text:
import Revise
import GLRenderer as GL
import ThreeDP3 as T
import Images as I
import MiniGSG as S
import Rotations as R
import PoseComposition: Pose
using Gen
import Serialization
import MeshCatViz

import MeshCatViz
MeshCatViz.setup_visualizer()

YCB_DIR = "/home/nishadg/mcs/ThreeDVision.jl/data/ycbv2"
world_scaling_factor = 100.0
id_to_cloud, id_to_shift, id_to_box  = T.load_ycbv_models_adjusted(YCB_DIR, world_scaling_factor);
all_ids = sort(collect(keys(id_to_cloud)));

# +
IDX =6
# IDX =301
IDX =2100
IDX =100

@show T.get_ycb_scene_frame_id_from_idx(YCB_DIR,IDX)

# +
gt_poses, ids, gt_rgb_image, gt_depth_image, cam_pose, camera = T.load_ycbv_scene_adjusted(
    YCB_DIR, IDX, world_scaling_factor, id_to_shift
);
dense_poses, dense_ids = T.load_ycbv_dense_fusion_predictions_adjusted(YCB_DIR, IDX, world_scaling_factor, id_to_shift);
dense_poses_reordered = [dense_poses[findfirst(dense_ids .== id)]   for id in ids];
inf_poses = Serialization.deserialize("../assets/results/$(lpad(IDX,4,"0")).poses").poses_after
inf_poses = T.adjust_poses(inf_poses, ids, world_scaling_factor, id_to_shift)

img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))
I.colorview(I.RGB, permutedims(Float64.(gt_rgb_image)./255.0, (3,1,2)))
# -

import Plots as P
import Measures

maximum(gt_depth_image)

P.heatmap(clamp.(gt_depth_image, 0.0,120.0), yflip=true, aspect_ratio=:equal,colorbar=false, xaxis=false, yaxis=false,xticks=false, yticks=false)
# P.savefig("depth.png")

size(gt_rgb_image)

img = I.colorview(I.Gray, gt_depth_image ./ maximum(gt_depth_image))

# +
resolution = 0.5

renderer_color = GL.setup_renderer(camera, GL.RGBMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer_color, v, n, f)
end
# -

colors = [
    I.colorant"yellow", I.colorant"cyan", I.colorant"lightgreen",
    I.colorant"red", I.colorant"purple", I.colorant"orange"
]

renderer_color.gl_instance.lightpos = [0.0, 0.0, -30.0]
images = [
    let
    rot = R.RotY(ang) * R.RotX(pi/2)
    gap = 15.0
    poses = [Pose([gap*i - 3.5*gap, 60.0, 0.0], rot) for (i,_) in enumerate(ids)]
    rgb_image, depth_image = GL.gl_render(
        renderer_color, ids, poses, colors[1:length(ids)], Pose([0.0,0.0,-160.0],R.RotX(-pi/10)))
    I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))
    end
    for ang in 0:0.1:(8*pi)
]

ang = 0.0
gap = 14.0
rot = R.RotY(ang) * R.RotX(pi/2)
p = popat!(poses, 4)
insert!(poses,2,p)
p = popat!(ids, 4)
insert!(ids,2,p)
p = popat!(colors, 4)
insert!(colors,2,p)
poses = [Pose([gap*i - 3.0*gap, 60.0, 0.0], rot) for (i,_) in enumerate(ids)]

# +

poses[5] = Pose(poses[5].pos, R.RotZ(pi/2))
poses[2] = Pose(poses[2].pos, R.RotZ(0.0))

rgb_image, depth_image = GL.gl_render(
    renderer_color, ids, poses, colors[1:length(ids)], Pose([0.0,0.0,-160.0],R.RotX(-pi/10)))
I.colorview(I.RGBA, permutedims(rgb_image,(3,1,2)))

# -

images[1]

I.colorview(I.RGB, permutedims(Float64.(gt_rgb_image)./255.0, (3,1,2)))

# +
resolution = 0.5

renderer = GL.setup_renderer(camera, GL.RGBMode())
for id in all_ids
    cloud = id_to_cloud[id]
    v,n,f = GL.mesh_from_voxelized_cloud(GL.voxelize(cloud, resolution), resolution)
    GL.load_object!(renderer, v, n, f)
end

# -

ids

new_poses = copy(gt_poses)
new_poses = [T.get_c_relative_to_a(cam_pose,p) for p in new_poses]
old_pose = new_poses[1]
new_poses[1] = Pose(old_pose.pos .- [0,0,2.0],old_pose.orientation)

renderer.gl_instance.lightpos = cam_pose.pos

rgb_image_dense, depth_image = GL.gl_render(
    renderer, ids, new_poses, colors, cam_pose)
rgb_image_dense= I.colorview(I.RGBA, permutedims(rgb_image_dense,(3,1,2)))

import Plots

import Pkg;Pkg.add("Measures")

using Measures

low, high = 56.0, 130.0
d = clamp.(depth_image, low, high)
@show minimum(d), maximum(d)
p = Plots.heatmap(d; c=:thermal,
    clim=(53.0, 130.0),
    ylim=(0, 480),
    xlim=(0, 640),
    aspect_ratio=:equal, yflip=true,  legend = :none, yticks=false, xticks=false, xaxis=false, yaxis=false, margin = 0mm)
Plots.savefig(p, "regen.pdf")
p

d = clamp.(gt_depth_image,low, high)
@show minimum(d), maximum(d)
p = Plots.heatmap(d; c=:thermal,
    clim=(53.0, 130.0),
    ylim=(0, 480),
    xlim=(0, 640),
    aspect_ratio=:equal, yflip=true,  legend = :none, yticks=false, xticks=false, xaxis=false, yaxis=false, margin = 0mm)
Plots.savefig(p, "real.pdf")
p

# +
focus_idx = 3
alpha = 0.4

rgb_image_dense, depth_image = GL.gl_render(
    renderer_color, [ids[focus_idx]], [dense_poses_reordered[focus_idx]], [I.colorant"cyan"], T.IDENTITY_POSE)
rgb_image_dense= I.colorview(I.RGBA, permutedims(rgb_image_dense,(3,1,2)))

new_gt_rgb_image = cat(Float64.(gt_rgb_image) ./ 255.0, ones(size(gt_rgb_image)[1:2]),dims=3)
new_gt_rgb_image = I.colorview(I.RGBA, permutedims(new_gt_rgb_image,(3,1,2)))
@show size(new_gt_rgb_image)
mask = (rgb_image_dense .!= I.RGBA(1.0, 1.0, 1.0, 1.0))
new_gt_rgb_image[mask] .= alpha.*new_gt_rgb_image[mask] +(1-alpha).*rgb_image_dense[mask]
new_gt_rgb_image
# +

rgb_image_3dp3, depth_image = GL.gl_render(
    renderer_color, [ids[focus_idx]], [inf_poses[focus_idx]], [I.colorant"limegreen"], T.IDENTITY_POSE)
rgb_image_3dp3= I.colorview(I.RGBA, permutedims(rgb_image_3dp3,(3,1,2)))

new_gt_rgb_image = cat(Float64.(gt_rgb_image) ./ 255.0, ones(size(gt_rgb_image)[1:2]),dims=3)
new_gt_rgb_image = I.colorview(I.RGBA, permutedims(new_gt_rgb_image,(3,1,2)))
@show size(new_gt_rgb_image)
@show size(rgb_image_3dp3)
mask = (rgb_image_3dp3 .!= I.RGBA(1.0, 1.0, 1.0, 1.0))
new_gt_rgb_image[mask] .= alpha.*new_gt_rgb_image[mask] +(1-alpha).*rgb_image_3dp3[mask]
new_gt_rgb_image
# -


text:
function draw!(camera::CameraModel, scene::AbstractScene)
    optical_center = Point3f0(camera.𝐜)
    image_width = camera.image_width
    image_height = camera.image_height
    f = camera.focal_length
    𝐞₁ = camera.𝐞₁
    𝐞₂ = camera.𝐞₂
    𝐞₃ = camera.𝐞₃
    bottom_right = optical_center + Point3f0((image_width/2)     * 𝐞₁ + (image_height/2)  * 𝐞₂ + f*𝐞₃)
    top_right =  optical_center   + Point3f0((image_width/2)     * 𝐞₁ + (-image_height/2) * 𝐞₂ + f*𝐞₃)
    top_left = optical_center     + Point3f0((-image_width/2)    * 𝐞₁ + (-image_height/2) * 𝐞₂ + f*𝐞₃)
    bottom_left = optical_center  + Point3f0((-image_width/2)    * 𝐞₁ + (image_height/2)  * 𝐞₂ + f*𝐞₃)

    centroid2film = [
        optical_center  => bottom_right;
        optical_center  => top_right;
        optical_center  => top_left;
        optical_center   => bottom_left;
    ]

    film = [
            bottom_right => top_right;
            top_right => top_left;
            top_left => bottom_left;
            bottom_left =>  bottom_right;
            ]

    scale = 20.0f0
    coordinate_system = [
        optical_center => optical_center + Point3f0(scale*𝐞₁);
        optical_center => optical_center + Point3f0(scale*𝐞₂);
        optical_center => optical_center + Point3f0(scale*𝐞₃);
    ]
    linesegments!(scene, centroid2film, color = :black, linewidth = 2)
    linesegments!(scene, film, color = :black, linewidth = 2)
    linesegments!(scene, coordinate_system, color = [:red, :green, :blue ], linewidth = 2)

end


text:
@inline Base.@propagate_inbounds Base.getindex(a::AbstractTen,I::Integer) = getfield(a,I)
@inline Base.@propagate_inbounds Base.getindex(a::AbstractTen,I::Integer,J::Integer) = getfield(a, I + 2*(J-1))

struct Ten2D{T<:Number} <: AbstractTen{T}
    xx::T
    yx::T
    xy::T
    yy::T
end

Base.size(a::Ten2D) =
    (2,2)

Base.length(a::Ten2D) =
    4

Base.linearindices(a::Ten2D) =
    Base.OneTo(4)

Base.indices(a::Ten2D) =
    (Base.OneTo(2),Base.OneTo(2))

Base.zero(a::Type{Ten2D{T}}) where {T} =
    Ten2D{T}(zero(T),zero(T),zero(T),zero(T))

Ten2D(xx,yx,xy,yy) =
    Ten2D(promote(xx,yx,xy,yy)...)

Base.:+(a::Ten2D{T},b::Ten2D{T2}) where {T,T2} =
    Ten2D{promote_type(T,T2)}(a.xx+b.xx, a.yx+b.yx, a.xy+b.xy, a.yy+b.yy)

Base.:-(a::Ten2D{T},b::Ten2D{T2}) where {T,T2} =
    Ten2D{promote_type(T,T2)}(a.xx-b.xx, a.yx-b.yx, a.xy-b.xy, a.yy-b.yy)

Base.:*(a::T,b::Ten2D{T2}) where {T<:Number,T2<:Number} =
    Ten2D{promote_type(T,T2)}(a*b.xx, a*b.yx, a*b.xy, a*b.yy)

Base.:*(b::Ten2D{T},a::T2) where {T<:Number,T2<:Number} =
    Ten2D{promote_type(T,T2)}(a*b.xx, a*b.yx, a*b.xy, a*b.yy)

Base.:/(b::Ten2D{T},a::T2) where {T<:Number,T2<:Number} =
    Ten2D{promote_type(T,T2)}(b.xx/a, b.yx/a, b.xy/a, b.yy/a)

Base.dot(a::Ten2D,b::Ten2D) =
    muladd(a.xx, b.xx, muladd(a.yx, b.yx, muladd(a.xy, b.xy, a.yy*b.yy)))

Base.norm(a::Ten2D) =
    sqrt(2dot(a,a))

Base.dot(a::Ten2D{T},b::Vec2D{T2}) where {T,T2} =
    Vec2D{promote_type(T,T2)}(a.xx*b.x + a.xy*b.y, a.yx*b.x + a.yy * b.y)

Base.dot(b::Vec2D{T},a::Ten2D{T2}) where {T,T2} =
    Vec2D{promote_type(T,T2)}(a.xx*b.x + a.yx*b.y, a.xy*b.x + a.yy * b.y)

Base.:*(a::Vec2D{T},b::Vec2D{T2}) where {T,T2} =
    Ten2D{promote_type(T,T2)}(a.x*b.x, a.y*b.x, a.x*b.y, a.y*b.y)


text:
@testset "cuda" begin
D = 10
Dcut = 5
d = 2
sites = 5

T = Float32
@testset "creation of random MPS on cuda" begin
    ψ = CUDA.randn(CuMPS{T}, sites, D, d)

    @test length(ψ) == sites
    @test ψ == ψ
    @test ψ ≈ ψ
    @test eltype(ψ) == T
    @test typeof(ψ.tensors[1]) <: CuArray
    @test copy(ψ) == ψ
end

@testset "random MPO on cuda" begin
    H = CUDA.randn(CuMPO{T}, sites, D, d)
    @test typeof(H.tensors[1]) <: CuArray
end

@testset "contraction on cuda" begin
    ψ = CUDA.randn(CuMPS{T}, sites, D, d)
    H = CUDA.randn(CuMPO{T}, sites, D, d)

    @test typeof((H*ψ).tensors[1]) <: CuArray
end
end


text:
module Dim1
import mFEM

# calc gauss_quad on interval [-1,1]
# u: (x::Float64)::Float64
function _gauss(u::Function)
  # c = sqrt(3/5)
  # xw = [-c 5/9; 0 8/9; c 5/9]

  a = 0.9061798
  b = 0.5384693
  c = 0
  wa = 0.2369269
  wb = 0.4786287
  wc = 0.5688889

  xw = [a wa; -a wa; b wb; -b wb; c wc]
  xwn,_ = size(xw)

  D = map(1:xwn) do i
    x,w = xw[i,:]
    return w * u(x)
  end
  return sum(D)
end

function uniform_mesh(a,b,N)::mFEM.Mesh
  @assert N > 1
  @assert a < b

  Nm = N + 1
  h = (b - a)/N

  P = [a + (i - 1)*h for i=1:Nm]
  T = zeros(Int,N,2)
  for i=1:N
    T[i,:] = [i,i+1]
  end

  BNode = [1,Nm]
  BEage = []
  BSurface = []

  # calc gauss_quad on the nth element of mesh
  # u: (x::Tuple)::Float64
  # n index of basis element
  function gauss_quad(u::Function,n::Int)
    a = P[T[n,1],1]
    b = P[T[n,2],1]

    @assert a < b
    return (b-a) * _gauss(x-> u(( (b-a)*x/2 + (a+b)/2, )) ) / 2
  end

  return mFEM.Mesh(P,T,BNode,BEage,BSurface,gauss_quad)
end

function Dϕ_reference_linear(k::Int, xt::Tuple, α::Int)
  @assert α ∈ [1,2]
  @assert k ∈ [0,1]

  local x = xt[1]

  if α == 1
    if k == 0
      return 1 - x
    elseif k == 1
      return  -1
    end
  elseif  α == 2
    if k == 0
      return x
    elseif k == 1
      return 1
    end
  end
end

function Dϕ_reference_quadratic(k::Int, xt::Tuple, α::Int)
  @assert α ∈ [1,2,3]
  @assert k ∈ [0,1,2]

  local x = xt[1]

  if α == 1
    if k == 0
      return 2x^2 - 3x + 1
    elseif k == 1
      return 4x - 3
    elseif k ==2
      return 4
    end
  elseif  α == 2
    if k == 0
      return 2x^2 - x
    elseif k == 1
      return 4x - 1
    elseif k ==2
      return 4
    end
  elseif  α == 3
    if k == 0
      return -4x^2 + 4x
    elseif k == 1
      return -8x + 4
    elseif k ==2
      return -8
    end
  end
end

function Dϕ_reference_cubic(k::Int, xt::Tuple, α::Int)
  @assert α ∈ 1:4
  @assert k ∈ 0:2

  local x = xt[1]

  if α == 1
    if k == 0
      return -4.5x^3 + 9x^2 -5.5x + 1
    elseif k == 1
      return -3*4.5x^2 + 2*9x -5.5
    elseif k == 2
      return -2*3*4.5x + 2*9
    end
  elseif α == 2
    if k == 0
      return 4.5x^3 - 4.5x^2 + x
    elseif k == 1
      return 3*4.5x^2 - 2*4.5x + 1
    elseif k == 2
      return 2*3*4.5x - 2*4.5
    end
  elseif α == 3
    if k == 0
      return 13.5x^3 - 22.5x^2 + 9x
    elseif k == 1
      return 3*13.5x^2 - 2*22.5x + 9
    elseif k == 2
      return 2*3*13.5x - 2*22.5
    end
  elseif α == 4
    if k == 0
      return -13.5x^3 + 18x^2 - 4.5x
    elseif k == 1
      return -3*13.5x^2 + 2*18x - 4.5
    elseif k == 2
      return -2*3*13.5x + 2*18
    end
  end
end

function reference2local(Dϕ_reference::Function;Pb::Array,Tb::Array,Nlb::Int)::Function
  # kth deriavtive of basis function at x
  function Dϕ(k::Int, xt::Tuple, n::Int, α::Int)
    @assert α ∈ 1:Nlb
    @assert k ∈ [0,1,2]
    @assert length(xt) == 1

    x = xt[1]

    xn1 = Pb[Tb[n,1],1]
    xn2 = Pb[Tb[n,2],1]

    x_hat = (x - xn1) / (xn2 - xn1)
    xt_hat = (x_hat,)

    if k == 0
      return Dϕ_reference(k,xt_hat,α)
    elseif k == 1
      return Dϕ_reference(k,xt_hat,α) / (xn2 - xn1)
    elseif k == 2
      @error "TODO"
    end
  end

  return Dϕ
end

function linear_elem(mesh::mFEM.Mesh)::mFEM.Basis
  Pb = copy(mesh.P)
  Tb = copy(mesh.T)
  Nlb = 2
  BNode = copy(mesh.BNode)
  BEage= copy(mesh.BEage)
  BSurface= copy(mesh.BSurface)
  BNode = copy(mesh.BNode)

  return mFEM.Basis(Nlb,Pb,Tb,BNode,BEage,BSurface,
    reference2local(Dϕ_reference_linear;Pb=Pb,Tb=Tb,Nlb=Nlb))
end

function quadatic_elem(mesh::mFEM.Mesh)::mFEM.Basis
  Pb = zeros(Float64,size(mesh.P,1) * 2 - 1,1)
  Tb = zeros(Int,size(mesh.T,1),3)
  Nlb = 3

  for i in 1:size(mesh.P,1)
    Pb[2i - 1] = mesh.P[i]
  end

  for i in 1:(size(mesh.P,1) - 1)
    Pb[2i] = (Pb[2i-1] + Pb[2i+1]) / 2
  end

  for i in 1:size(mesh.T,1)
    Tb[i,1] = 2 * mesh.T[i,1] - 1
    Tb[i,2] = 2 * mesh.T[i,2] - 1
    Tb[i,3] = (Tb[i,1] + Tb[i,2]) / 2
  end

  BNode = zeros(Int,size(mesh.BNode,1),1)
  for i in 1:size(BNode,1)
    BNode[i] = 2 * mesh.BNode[i] - 1
  end

  BEage = copy(mesh.BEage)
  BSurface = copy(mesh.BSurface)

  return mFEM.Basis(Nlb,Pb,Tb,BNode,BEage,BSurface,
    reference2local(Dϕ_reference_quadratic;Pb=Pb,Tb=Tb,Nlb=Nlb)
  )
end

function cubic_elem(mesh::mFEM.Mesh)::mFEM.Basis
  Pb = zeros(Float64,size(mesh.P,1) * 3 - 2,1)
  Tb = zeros(Int,size(mesh.T,1),4)
  Nlb = 4

  for i in 1:size(mesh.P,1)
    Pb[3i - 2] = mesh.P[i]
  end
  for i in 1:(size(mesh.P,1) - 1)
    Pb[3i - 1] = (2/3)* Pb[3i-2] + (1/3) * Pb[3i + 1]
    Pb[3i] = (1/3)* Pb[3i-2] + (2/3) * Pb[3i + 1]
  end

  for i in 1:size(mesh.T,1)
    Tb[i,1] = 3 * mesh.T[i,1] - 2
    Tb[i,2] = 3 * mesh.T[i,2] - 2
    Tb[i,3] = Tb[i,1] + 1
    Tb[i,4] = Tb[i,1] + 2
  end

 BNode = zeros(Int,size(mesh.BNode,1),1)
  for i in 1:size(BNode,1)
    BNode[i] = 3 * mesh.BNode[i] - 2
  end

  BEage = copy(mesh.BEage)
  BSurface = copy(mesh.BSurface)

  return mFEM.Basis(Nlb,Pb,Tb,BNode,BEage,BSurface,
    reference2local(Dϕ_reference_cubic;Pb=Pb,Tb=Tb,Nlb=Nlb)
  )
end

function Possion(mesh::mFEM.Mesh,c::Function,f::Function,g::Function;
    trial::mFEM.Basis = linear_elem(mesh),test::mFEM.Basis = linear_elem(mesh))
  # generate boundary nodes
  boundarynodes = zeros(Int,size(trial.BNode,1),2)
  for i=1:size(boundarynodes,1)
    boundarynodes[i,:] = [-1,trial.BNode[i]]
  end

  A = mFEM.stiffness_mat(mesh,trial,test,c)
  b = mFEM.load_vec(mesh,test,f)

  mFEM.set_boundary_cond!(trial,boundarynodes,A,b,g)
  x = A \ b
  return x
end

function test(n,;finite_elem::Function=linear_elem)
  u(x) = x * cos(x)
  c(x) = exp(x)
  f(x) = - exp(x) *(cos(x) - 2sin(x) - x*cos(x) - x*sin(x))
  g(x,i) = u(x[1])

  function tuplfy1d(f::Function)::Function
    return (x::Tuple) -> f(x[1])
  end

  m = uniform_mesh(0,1,n)
  trial = finite_elem(m)
  test= finite_elem(m)

  U = Possion(m,tuplfy1d(c),tuplfy1d(f),g;trial=trial,test=test)
  return maximum(@.abs(U-map(u,trial.Pb)))
end

test_quadatic(n) = test(n,finite_elem=quadatic_elem)
test_linear(n) = test(n,finite_elem=linear_elem)
test_cubic(n) = test(n,finite_elem=cubic_elem)

end # Dim1


text:
using PowerModels
using JuMP

# post a JuMP model with loads as variables and additional constraints given in the dictionaries A and B
function post_ac_opf_withload(data, model, A, B)
    @assert !(data["multinetwork"])
    ref = PowerModels.build_ref(data)[:nw][0]

    @variable(model, pd[i in keys(ref[:bus])], start = ref[:bus][i]["pd"])
    @variable(model, qd[i in keys(ref[:bus])], start = ref[:bus][i]["qd"])

    @variable(model, va[i in keys(ref[:bus])])
    @variable(model, ref[:bus][i]["vmin"] <= vm[i in keys(ref[:bus])] <= ref[:bus][i]["vmax"], start=1.0)

    @variable(model, ref[:gen][i]["pmin"] <= pg[i in keys(ref[:gen])] <= ref[:gen][i]["pmax"])
    @variable(model, ref[:gen][i]["qmin"] <= qg[i in keys(ref[:gen])] <= ref[:gen][i]["qmax"])

    @variable(model, -ref[:branch][l]["rate_a"] <= p[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])
    @variable(model, -ref[:branch][l]["rate_a"] <= q[(l,i,j) in ref[:arcs]] <= ref[:branch][l]["rate_a"])

    @variable(model, ref[:arcs_dc_param][a]["pmin"] <= p_dc[a in ref[:arcs_dc]] <= ref[:arcs_dc_param][a]["pmax"])
    @variable(model, ref[:arcs_dc_param][a]["qmin"] <= q_dc[a in ref[:arcs_dc]] <= ref[:arcs_dc_param][a]["qmax"])

    from_idx = Dict(arc[1] => arc for arc in ref[:arcs_from_dc])
    @objective(model, Min,
        sum(gen["cost"][1]*pg[i]^2 + gen["cost"][2]*pg[i] + gen["cost"][3] for (i,gen) in ref[:gen]) +
        sum(dcline["cost"][1]*p_dc[from_idx[i]]^2 + dcline["cost"][2]*p_dc[from_idx[i]] + dcline["cost"][3] for (i,dcline) in ref[:dcline])
    )

    for (i,bus) in ref[:ref_buses]
        # Refrence Bus
        @constraint(model, va[i] == 0)
    end

    for (i,bus) in ref[:bus]
        # Bus KCL
        @constraint(model,
            sum(p[a] for a in ref[:bus_arcs][i]) +
            sum(p_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) ==
            sum(pg[g] for g in ref[:bus_gens][i]) - pd[i] - bus["gs"]*vm[i]^2
        )
        @constraint(model,
            sum(q[a] for a in ref[:bus_arcs][i]) +
            sum(q_dc[a_dc] for a_dc in ref[:bus_arcs_dc][i]) ==
            sum(qg[g] for g in ref[:bus_gens][i]) - qd[i] + bus["bs"]*vm[i]^2
        )
    end

    for (i,branch) in ref[:branch]
        f_idx = (i, branch["f_bus"], branch["t_bus"])
        t_idx = (i, branch["t_bus"], branch["f_bus"])

        p_fr = p[f_idx]
        q_fr = q[f_idx]
        p_to = p[t_idx]
        q_to = q[t_idx]

        vm_fr = vm[branch["f_bus"]]
        vm_to = vm[branch["t_bus"]]
        va_fr = va[branch["f_bus"]]
        va_to = va[branch["t_bus"]]

        # Line Flow
        g, b = PowerModels.calc_branch_y(branch)
        tr, ti = PowerModels.calc_branch_t(branch)
        c = branch["br_b"]
        tm = branch["tap"]^2

        # AC Line Flow Constraints
        @NLconstraint(model, p_fr == g/tm*vm_fr^2 + (-g*tr+b*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-b*tr-g*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )
        @NLconstraint(model, q_fr == -(b+c/2)/tm*vm_fr^2 - (-b*tr-g*ti)/tm*(vm_fr*vm_to*cos(va_fr-va_to)) + (-g*tr+b*ti)/tm*(vm_fr*vm_to*sin(va_fr-va_to)) )

        @NLconstraint(model, p_to == g*vm_to^2 + (-g*tr-b*ti)/tm*(vm_to*vm_fr*cos(va_to-va_fr)) + (-b*tr+g*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )
        @NLconstraint(model, q_to == -(b+c/2)*vm_to^2 - (-b*tr+g*ti)/tm*(vm_to*vm_fr*cos(va_fr-va_to)) + (-g*tr-b*ti)/tm*(vm_to*vm_fr*sin(va_to-va_fr)) )

        # Phase Angle Difference Limit
        @constraint(model, va_fr - va_to <= branch["angmax"])
        @constraint(model, va_fr - va_to >= branch["angmin"])

        # Apparent Power Limit, From and To
        @constraint(model, p[f_idx]^2 + q[f_idx]^2 <= branch["rate_a"]^2)
        @constraint(model, p[t_idx]^2 + q[t_idx]^2 <= branch["rate_a"]^2)
    end

    for (i,dcline) in ref[:dcline]
        # DC Line Flow Constraint
        f_idx = (i, dcline["f_bus"], dcline["t_bus"])
        t_idx = (i, dcline["t_bus"], dcline["f_bus"])

        @constraint(model, (1-dcline["loss1"])*p_dc[f_idx] + (p_dc[t_idx] - dcline["loss0"]) == 0)
    end

    # adding additional constraints as contained in the dictionaries A and B
    for i in keys(ref[:bus])
        a = A[i]
        b = B[i]
        @constraint(model, a*[pd[i],qd[i]] .<= b)
    end

    var_refs = Dict{String,Any}()
    var_refs["pg"] = pg
    var_refs["qg"] = qg

    return model, var_refs
end

max_stars_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_issues_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed

max_forks_repo_name: bishnubhattarai/PFO-ADC-DER-Testbed


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:4] <= 1)

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :(swish(0.7464435003747418 + -0.7676428720221384*swish(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.5937285886485864*swish(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.70936536775501*swish(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.9647291374973435*swish(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.527810271934432*swish(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + swish(-0.35636357579153266 + 0.8665784119884141*swish(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.03175667478129274*swish(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.7799384758999897*swish(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.4104070349373172*swish(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + -0.7416668360909755*swish(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + swish(0.5038932713674513 + 0.7223377417739489*swish(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.17110213406275054*swish(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.7726394930743403*swish(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.6118790711423934*swish(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.16250007841990932*swish(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + swish(-0.7190092552996492 + 0.6952008118484776*swish(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + 0.975536624348428*swish(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + 0.03367235210941821*swish(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.8059867473566515*swish(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.5874603423029647*swish(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) + swish(-0.5715134646750415 + 0.6325466436773857*swish(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4])) + -0.2633742432221431*swish(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4])) + -0.3786600783719858*swish(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4])) + 0.5923927477083897*swish(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4])) + 0.9114166052727355*swish(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions