text:
#!/usr/bin/julia
#
# Plotting functions.
#
using RecipesBase;
using Plots: @layout, grid;
#using IJulia; # Required to prevent plotting window opening and crashing.

export check_path, gen_filename, format_plot_title_variable,
       GeneralPlotOptions, LanczosPlotOptions, SplitStepPlotOptions, SimPlot;

"""
Implementations contain mainly boolean objects specifying what should be
included within plots.

All implementations should at least implement all the properties of
`GeneralPlotOptions`.

See also: [`GeneralPlotOptions`](@ref).
"""
abstract type PlotOptions end;

"Plotting options for a simulation agnostic plot"
struct GeneralPlotOptions <: PlotOptions
  """
  If the real and imaginary parts of the wave function should be plotted in
  addition to the absolute square value.
  """
  components::Bool;

  """
  If the value of the wave function is less than the value of this property at
  the start or end of the plot range, then they will not be plotted.
  """
  extremaϵ::Number;

end
GeneralPlotOptions(;components::Bool=false, extremaϵ::Number=0.1) =
  GeneralPlotOptions(components, extremaϵ);

"Plotting options for split step simulations."
struct SplitStepPlotOptions <: PlotOptions
  components::Bool;
  extremaϵ::Number;
end
SplitStepPlotOptions(;components::Bool=false, extremaϵ::Number=0.1) =
  SplitStepPlotOptions(components, extremaϵ);

"Plotting options for Lanczos simulations."
struct LanczosPlotOptions <: PlotOptions
  components::Bool;
  extremaϵ::Number;
  "Plot basis functions."
  basis::Bool;
  "Plot eigenfunctions of the reduced Hamiltonian."
  eigenfct::Bool;
end

abstract type PlotType end;
struct SimPlot<:PlotType end;

"""
    check_path(path::String)

Check the path given exists and create it if it does not.
"""
function check_path(path::String)
  if (!ispath(path))
    mkpath(path)
  end
end

# Parameters common to all simulation methods.
function gen_filename(p::SimParams, grid::UniformGrid{S,1}) where S
  return string(real(p.Δt) == 0 ? "ti" : "td",
    round(typeof(grid.Δ), grid.Δ, sigdigits=4), "d",
    round(typeof(grid.max), grid.max, sigdigits=4), "mx",
    round(typeof(p.Δt), abs(p.Δt), sigdigits=4), "dt-", p.steps, "s");
end
# Simulation variation method.
function gen_filename(var::Vararg{SimVariation})
  name = string();

  # Iterate through simulation variations and append parameters to filename.
  for i ∈ var
    if i isa AdaptiveΔt
      name *= string("-", round(i.ϵmax, sigdigits=4), "emx",
                     round(i.ϵmin, sigdigits=4), "emn",
                     round(i.Δtdown, sigdigits=4), "tdn",
                     round(i.Δtup, sigdigits=4), "tup",
                     round(i.Δtmax, sigdigits=4), "tmx",
                     round(i.Δtmin, sigdigits=4), "tmn",
                     i.maxeval, "mxev");
      if i.maxevallimit != nothing
        name *= string(i.maxevallimit, "mxevlm");
      end
    elseif i isa ConvergenceVariation
      name *= string("-", i.ϵ, "ecnv");
    end
  end

  return name;
end
# Common plotting options.
function gen_filename(plt::PlotOptions)
  if plt.components
    return string("pcmpnt");
  end
  return string();
end
# Simulation specific methods.
function gen_filename(sim::Lanczos, plt::LanczosPlotOptions)
  name = string("-sil-", sim.basisno, "b", sim.basissteps, "bs-" );

  if plt.basis name *= string("pb"); end
  if plt.eigenfct name *= string("pefct"); end

  name *= gen_filename(plt);

  return name;
end
function gen_filename(sim::SplitStep, plt::SplitStepPlotOptions)
  name = string("-ss-");

  name *= gen_filename(plt);

  return name
end
"""
    gen_filename(p::SimParams, grid::UniformGrid{S,1},
                      sim::SimMethod, plt::PlotOptions,
                      var::Union{Vararg{SimVariation},Nothing}=nothing) where S

Generate an appropriate filename reflecting the parameters used in the
simulation.
"""
function gen_filename(p::SimParams, grid::UniformGrid{S,1},
                      sim::SimMethod, plt::PlotOptions,
                      var::Vararg{Union{Nothing,SimVariation}}=nothing) where S
  name = gen_filename(p, grid);
  name *= gen_filename(sim, plt);
  if (var[1] != nothing) name *= gen_filename(var...); end

  return name;
end

"""
    format_plot_title_variable(variable::Number, decimalplaces::Integer)

Format `variable` with a fixed number of decimal places given by `decimalplaces`
such that it is in a suitable format to be printed in a GIF frame.
"""
function format_plot_title_variable(variable::Unitful.Quantity, decimalplaces::Integer)
  return string(format_plot_title_variable(ustrip(variable), decimalplaces),
                " ", unit(variable));
end
function format_plot_title_variable(variable::Number, decimalplaces::Integer)
  # Determine number of digits before the decimal place.
  predecimal = findfirst(".", string(variable));
  if predecimal === nothing predecimal = length(string(variable));
  else predecimal = predecimal[1]-1 end

  return string(rpad(round(variable, digits=decimalplaces),
                     predecimal+decimalplaces+1, '0'));
end

"""
    get_effective_minmax(v::AbstractVector, ϵ::Float64=0.1)

Return the minimum and maximum indices of `v` below and above which the value of
`v` does not exceed `ϵ`.
"""
function get_effective_minmax(v::AbstractVector, ϵ::Number=0.1)
  len = length(v);
  min = 1;
  max = len;
  # Find the first and last indices of v where v > ϵ.
  for i ∈ 1:len
    if abs(v[i]) > ϵ
      min = i;
      break;
    end;
  end
  for i ∈ len:-1:1
    if abs(v[i]) > ϵ
      max = i;
      break;
    end
  end

  return (min, max);
end
function get_effective_minmax(v::Vararg{AbstractVector}; ϵ=0.1)
  len = length(v);
  mins = Vector{Int}(undef, len);
  maxs = Vector{Int}(undef, len);

  for i ∈ 1:len
    (mins[i], maxs[i]) = get_effective_minmax(v[i], ϵ);
  end

  return (mins, maxs);
end
"""
    get_effective_minmax(m::AbstractMatrix, ϵ::Number=0.1)

Treat the columns of the matrix `m` as vectors.
"""
function get_effective_minmax(m::AbstractMatrix, ϵ::Number=0.1)
  numv = size(m)[2];
  mins = Vector{Int}(undef, numv);
  maxs = Vector{Int}(undef, numv);

  for i ∈ 1:numv
    (mins[i], maxs[i]) = get_effective_minmax(m[:,i], ϵ);
  end

  return (mins, maxs);
end

"""
    get_global_effective_minmax(v::Vararg{AbstractVector}; ϵ=0.1)

Returns the maximum and minimum values found by the respective
`get_effective_minmax()` function.

See also: [`get_effective_minmax`](@ref)
"""
function get_global_effective_minmax(v::Vararg{AbstractVector}; ϵ=0.1)
  (xmins, xmaxs) = get_effective_minmax(v..., ϵ=ϵ);
  return (minimum(xmins), maximum(xmaxs));
end
function get_global_effective_minmax(m::AbstractMatrix, ϵ::Number=0.1)
  (xmins, xmaxs) = get_effective_minmax(m, ϵ);
  return (minimum(xmins), maximum(xmaxs));
end

"""
    get_ranged_minmax(min::Int, max::Int, v::Vararg{AbstractVector})

Find the minimum and maximum of the `v` in the index range `min:max`.

Note: values of the `v` should be real to allow size comparison.
"""
function get_ranged_minmax(min::Int, max::Int, v::Vararg{AbstractVector})
  ymin = 0.0;
  ymax = 0.0;
  for i ∈ 2:numψ
    w = v[i][xmin:xmax];
    t = minimum(w);
    if (t < ymin || i == 1) ymin = t; end
    t = maxmum(w);
    if (t > ymmax || i ==1) ymax = t; end
  end

  return (ymin, ymax);
end

"""
    subscript(i::Integer)

Generate a string with a subscript number `i`. From Stack Overflow user
improbable https://stackoverflow.com/a/46709534.
"""
subscript(i::Integer) = i<0 ? error("$i is negative") : join('₀'+d for d in reverse(digits(i)));

"""
    generate_wavefunction_labels(num::Int)

Generate labels to identify `num` wave functions.
"""
function generate_wavefunction_labels(num::Int)
  labels = Vector{String}(undef, num);
  if (num == 1) labels[1] = "|ψ|²";
  else
    for i ∈ 1:num labels[i] = string("|ψ", subscript(i), "|²"); end
  end
  return reshape(labels, 1, num);
end

"Set appropriate axis labels"
@recipe function f(::Type{O}, x::Vector, y::Array) where
                   {T,N,S<:ConfigurationSpace,O<:WaveObject{T,N,S,1}}
  xguide --> "r";
  yguide --> "Probability density";
  (x, y)
end
@recipe function f(::Type{O}, x::Vector, y::Array) where
                   {T,N,S<:MomentumSpace,O<:WaveObject{T,N,S,1}}
  xguide --> "k";
  yguide --> "Probability density";
  (x, y)
end

"""
    vectors_to_matrix(v::Tuple{Vararg{Vector}}, min::Int=1, max::Int=0)

Construct a matrix from the vectors `v` between the indices `min` and `max`.
"""
function vectors_to_matrix(v::Tuple{Vararg{AbstractVector}}, min::Int=1, max::Int=0)
  num = length(v);
  if (max == 0) max = num; end

  m = Matrix{Number}(undef, max - min + 1, num);
  for i ∈ 1:num m[:,i] = abs2.(v[i][min:max]); end

  return m;
end

"Plot absolute square of wave function across relevant range."
@recipe function f(grid::Grid1D, ψ::Vararg{WaveFunction{T,1}};
                   opts=GeneralPlotOptions()) where T
  (opts, grid, ψ...);
end
# Keyword arguments are not dispatched properly in the recipe chain,
# so define this alternative order of the arguments to implement the function.
# 2020/12/01
@recipe function f(opts::PlotOptions, grid::Grid1D,
                   ψ::Vararg{WaveFunction{T,1}}) where T
  label --> generate_wavefunction_labels(length(ψ));

  (min, max) = get_global_effective_minmax(ψ..., ϵ=opts.extremaϵ);

  (typeof(ψ[1]), grid[min:max], vectors_to_matrix(ψ, min, max));
end

"""
Plot potential function and absolute square of wave function across interesting
wave function range.
"""
@recipe function f(grid::Grid1D, H::SeparableHamiltonian,
                   ψ::Vararg{WaveFunction{T,1}};
                   opts=GeneralPlotOptions()) where T
  (opts, grid, H, ψ...);
end
# Workaround for same issue as above.
@recipe function f(opts::PlotOptions, grid::Grid1D, H::SeparableHamiltonian,
                   ψ::Vararg{WaveFunction{T,1}}) where T
  label --> hcat(["Potential"], generate_wavefunction_labels(length(ψ)));

  (min, max) = get_global_effective_minmax(ψ..., ϵ=opts.extremaϵ);

  (typeof(ψ[1]), grid[min:max], hcat(real.(H.V[min:max]),
                                     vectors_to_matrix(ψ, min, max)));
end

"""
Plot absolute square of wave function components of wave matrix across relevant
range.
"""
@recipe function f(grid::G, ϕ::WaveMatrix{T,S,G},
                   opts::PlotOptions=GeneralPlotOptions()) where {T,S,G}
  label --> generate_wavefunction_labels(size(ϕ)[2]);

  (min, max) = get_global_effective_minmax(ϕ, opts.extremaϵ);

  (typeof(ϕ), grid[min:max], abs2.(ϕ[min:max,:]))
end

"""
Plot potential function and absolute square of wave function components of wave
matrix across relevant range.
"""
@recipe function f(grid::G, H::SeparableHamiltonian, ϕ::WaveMatrix{T,S,G},
                   opts::PlotOptions=GeneralPlotOptions()) where {T,S,G}
  label --> hcat(["Potential"], generate_wavefunction_labels(size(ϕ)[2]));

  (min, max) = get_global_effective_minmax(ϕ, opts.extremaϵ);

  (typeof(ϕ), grid[min:max], hcat(real.(H.V[min:max]),
                                  abs2.(ϕ[min:max,:])))
end

"Plot potential functions included in the passed Hamiltonians."
@recipe function f(grid::Grid1D, H::Vararg{SeparableHamiltonian})
  potentials = Vector{Vector{real(eltype(H[1]))}}(undef, length(H));
  for i ∈ 1:length(H) potentials[i] = real.(H[i].V); end
  (grid, potentials)
end

"Plot real and imaginary parts of wave functions alongside main plot."
@recipe function f(p::SimPlot, opts::PlotOptions, g::G, ψ::WaveFunction{T,1,S,G},
                   H::Union{SeparableHamiltonian,Nothing}=nothing) where {T,S,G}
  size := (800, 450);
  dpi := 300;

  if opts.components
    #layout := @layout [ a{0.6w} grid(2,1) ]
    layout := @layout [ a{0.5w} b ]

    @series begin
      subplot := 2;
      title := "Wave function decomposition";
      ylims := (-1,1);
      realψ = real.(ψ);
      (remin, remax) = get_effective_minmax(realψ, opts.extremaϵ);
      imagψ = imag.(ψ);
      (immin, immax) = get_effective_minmax(imagψ, opts.extremaϵ);
      min = Base.min(immin, remin);
      max = Base.max(immax, remax);
      label := hcat(["Re(ψ)"], ["Im(ψ)"]);
      (g[min:max], hcat(realψ[min:max], imagψ[min:max]))
    end
    # Allow negative y values on subplots.
    #=
    try
      subylims = ylims();
      if subylims[1] >= 0
        subylims = (-subylims[2], subylims[2]);
      end
    catch
    end=#

    #=@series begin
      subplot := 2;
      title := "Re(ψ)";
      ylims := (-1,1);
      realψ = real.(ψ);
      (min, max) = get_effective_minmax(realψ, opts.extremaϵ);
      (g[min:max], realψ[min:max])
    end

    @series begin
      subplot := 3;
      title := "Im(ψ)";
      ylims := (-1,1);
      imagψ = imag.(ψ);
      (min, max) = get_effective_minmax(imagψ, opts.extremaϵ);
      (g[min:max], imagψ[min:max])
    end=#
  end

  @series begin
    subplot := 1;
    if (H === nothing) return (opts, g, ψ);
    else return (opts, g, H, ψ); end
  end
end

"Plot Krylov basis and reduced Hamiltonian eigenfunctions from Lanczos method."
@recipe function f(p::SimPlot, opts::LanczosPlotOptions, g::G, ψ::WaveFunction{T,1,S,G},
                   Φ::KrylovBasis{T,S,G}, H::Union{SeparableHamiltonian,Nothing}=nothing) where {T,S,G}
  if opts.components && opts.basis && opts.eigenfct
    layout := @layout [ a{0.7w} grid(2,1); grid(1,2){0.4h}]
    size := (1400, 750);
  elseif !opts.components && opts.basis && opts.eigenfct
    layout := @layout [ a; grid(1,2)]
    size := (1400, 750);
  elseif !opts.components && !opts.basis && !opts.eigenfct
    size := (1400, 750);
  elseif opts.components && !opts.basis && !opts.eigenfct
    layout := @layout [ a{0.6w} grid(2,1) ]
    size := (1400, 750);
  else
    error("the requested plotting format has not yet been implemented.");
  end

  if opts.components
    @series begin
      subplot := 2;
      title := "Re(ψ)";
      ylims := (-1,1);
      realψ = real.(ψ);
      (min, max) = get_effective_minmax(realψ, opts.extremaϵ);
      (g[min:max], realψ[min:max])
    end
    @series begin
      subplot := 3;
      title := "Im(ψ)";
      ylims := (-1,1);
      imagψ = imag.(ψ);
      (min, max) = get_effective_minmax(imagψ, opts.extremaϵ);
      (g[min:max], imagψ[min:max])
    end
  end
  if opts.basis && opts.eigenfct
    # Plot basis functions.
    @series begin
      if !opts.components subplot := 2;
      else subplot := 4; end
      title := "Krylov Basis"
      (g, Φ, opts)
    end

    # Plot eigenfunctions.
    @series begin
      if !opts.components subplot := 3;
      else subplot := 5; end

      # Label the eigenstates by their energies.
      len = length(Φ.λ);
      labels = Vector{String}(undef, len)
      for i ∈ 1:len labels[i] = string("|ψ|² where E = ", round(Φ.λ[i], sigdigits=4)); end
      label := reshape(labels, 1, len);

      title := "Reduced Hamiltonian Eigenstates"
      (g, KrylovBasis{T,S,G}(Φ*Φ.v, Φ.λ, Φ.v, Φ.Δ), opts)
    end
  end

  @series begin
    subplot := 1;
    if (H === nothing) return (opts, g, ψ);
    else return (opts, g, H, ψ); end
  end
end


text:
"""
    BoundingBox(pMin::Point3, pMax::Point3)

A bounding box is defined by its extreme corners `pMin` and `pMax`.

"""
struct BoundingBox <: BoundingVolume
    pMin::Point3
    pMax::Point3
    function BoundingBox(p1::Point3, p2::Point3)
        pMin = min(p1, p2)
        pMax = max(p1, p2)
        return new(pMin, pMax)
    end
end

"""
    BoundingBox(boxes::Array{BoundingBox,1})

Construct the bounding box of a list of bounding boxes.

"""
function BoundingBox(boxes::Array{BoundingBox,1})
    pMin = Point3(Inf, Inf, Inf)
    pMax = Point3(-Inf, -Inf, -Inf)
    for bb in boxes
        pMin = min(pMin, bb.pMin)
        pMax = max(pMax, bb.pMax)
    end
    return BoundingBox(pMin, pMax)
end

"""
    BoundingBox(p1::Point3, p2::Point3, p3::Point3)

Construct the bounding box of three points.

"""
function BoundingBox(p1::Point3, p2::Point3, p3::Point3)
    pMin = min(p1, p2, p3)
    pMax = max(p1, p2, p3)
    return BoundingBox(pMin, pMax)
end

"""
    BoundingBox(contents::Array{Point3,1})

Construct the bounding box of a list of points.

"""
function BoundingBox(contents::Array{Point3,1})
    pMin = Point3(Inf, Inf, Inf)
    pMax = Point3(-Inf, -Inf, -Inf)
    for p in contents
        pMin = min(pMin, p)
        pMax = max(pMax, p)
    end
    return BoundingBox(pMin, pMax)
end

"""
    BoundingBox(BB::BoundingBox, p::Point3)

Return a bounding box that covers both the box `BB` and the point `p`.

"""
BoundingBox(BB::BoundingBox, p::Point3) = BoundingBox(min(BB.pMin, p), max(BB.pMax, p))

"""
    BoundingBox(p::Point3, BB::BoundingBox)

Return a bounding box that covers both the box `BB` and the point `p`.

"""
BoundingBox(p::Point3, BB::BoundingBox) = BoundingBox(BB, p)

"""
    BoundingBox(BB1::BoundingBox, BB2::BoundingBox)

Return a bounding box that covers the two bounding boxes given.

"""
BoundingBox(BB1::BoundingBox, BB2::BoundingBox) = BoundingBox(min(BB1.pMin, BB2.pMin), max(BB1.pMax, BB2.pMax))

"""
    (transform::Transformation)(BB::BoundingBox)

Apply a transformation on a bounding box (probably not in the most efficient way).

"""
function (transform::Transformation)(BB::BoundingBox)
    corners = Point3[
        Point3(BB.pMin.x, BB.pMin.y, BB.pMin.z),
        Point3(BB.pMin.x, BB.pMin.y, BB.pMax.z),
        Point3(BB.pMin.x, BB.pMax.y, BB.pMin.z),
        Point3(BB.pMin.x, BB.pMax.y, BB.pMax.z),
        Point3(BB.pMax.x, BB.pMin.y, BB.pMin.z),
        Point3(BB.pMax.x, BB.pMin.y, BB.pMax.z),
        Point3(BB.pMax.x, BB.pMax.y, BB.pMin.z),
        Point3(BB.pMax.x, BB.pMax.y, BB.pMax.z)
    ]
    return BoundingBox([transform(p) for p in corners])
end

"""
    intersect(R::Ray, BB::BoundingBox)

Ray intersection with a bounding box.

Returns `(hit::Bool, tmin::Float64, tmax::Float64)`, where tmin and tmax are the near and far
hit locations on the box. They will be NaN if hit is false.

"""
function intersect(R::Ray, BB::BoundingBox)
    t0 = R.tmin
    t1 = R.tmax
    for i = 1:3
        tNear = (BB.pMin[i] - R.origin[i]) / R.direction[i]
        tFar = (BB.pMax[i] - R.origin[i]) / R.direction[i]
        tNear, tFar = min(tNear, tFar), max(tNear, tFar)
        t0 = tNear > t0 ? tNear : t0
        t1 = tFar < t1 ? tFar : t1
        if t0 > t1
            return false, NaN, NaN
        end
    end
    return true, t0, t1
end

"""
    intersectP(R::Ray, BB::BoundingBox)

True/false intersection test with bounding box.

"""
function intersectP(R::Ray, BB::BoundingBox)
    t0 = R.tmin
    t1 = R.tmax
    for i = 1:3
        invd = 1.0 / R.direction[i]
        tNear = (BB.pMin[i] - R.origin[i]) * invd
        tFar = (BB.pMax[i] - R.origin[i]) * invd
        tNear, tFar = min(tNear, tFar), max(tNear, tFar)
        t0 = tNear > t0 ? tNear : t0
        t1 = tFar < t1 ? tFar : t1
        if t0 > t1
            return false
        end
    end
    return true
end

"""
    area(BB::BoundingBox)

Return the surface area of the bounding box.

"""
area(BB::BoundingBox) = 2 * ((BB.pMax[1] - BB.pMin[1]) * (BB.pMax[2] - BB.pMin[2])
                           + (BB.pMax[2] - BB.pMin[2]) * (BB.pMax[3] - BB.pMin[3])
                           + (BB.pMax[3] - BB.pMin[3]) * (BB.pMax[1] - BB.pMin[1]))

"""
    max_extent(BB::BoundingBox)

Get the index of the axis where the bounding box is widest. X is 1, Y is 2 and Z is 3.

"""
function max_extent(BB::BoundingBox)
    dx = BB.pMax.x - BB.pMin.x
    dy = BB.pMax.y - BB.pMin.y
    dz = BB.pMax.z - BB.pMin.z
    if dx > dy && dx > dz
        return 1
    elseif dy > dx && dy > dz
        return 2
    else
        return 3
    end
end

##### BOUNDING SPHERE

"""
    BoundingSphere(center::Point3, radius::Float64)

A spherical bounding volume.

"""
struct BoundingSphere <: BoundingVolume
    center::Point3
    radius::Float64
end

"""
    BoundingSphere(BB:BoundingBox ; r_eps=1e-15)

Construct a `BoundingSphere` which circumscribes the given `BoundingBox`.

The small value `r_eps` is added to the computed radius to avoid clipping with the corners
of the box.

"""
function BoundingSphere(BB::BoundingBox ; r_eps=1e-15)
    c = mean(BB.pMin, BB.pMax)
    r = norm(c - BB.pMax) + r_eps
    return BoundingSphere(c, r)
end

BoundingSphere(stuff ; r_eps=1e-15) = BoundingSphere(BoundingBox(stuff) ; r_eps = r_eps)

"""
    intersectP(R::Ray, BS::BoundingSphere)

True/false intersection test between ray and given `BoundingSphere`.

"""
function intersectP(R::Ray, BS::BoundingSphere)
    A = dot(R.direction, R.direction)
    B = 2.0 * dot(R.direction, R.origin)
    C = dot(R.origin, R.origin) - BS.radius^2
    hit, tnear, tfar = quadratic(A, B, C)

    if !hit || tnear > R.tmax || tfar < R.tmin || (tnear < R.tmin && tfar > R.tmax)
        return false
    else
        return true
    end
end

@inline cross_section(B::BoundingSphere) = π * B.radius^2

ext: jl

max_stars_repo_name: dronir/GRAYTR

max_issues_repo_name: dronir/GRAYTR

max_forks_repo_name: dronir/GRAYTR


text:
using FileIO, DelimitedFiles
using Images
# using NPZ
using LinearAlgebra
import Rotations

"""
time stamps are in μs
"""
mutable struct TumData <:RgbdReader
    k::Int64
    frame_num::Int64
    folder::String
    ids::Vector{Int64}
    timestamps::Vector{Int64}
    id2rgb::Dict{Int64, String}
    id2depth::Dict{Int64, String}
    cam::CameraModel.RgbdCamParams
end

"""
at k = 0 means it has not started yet
"""
function TumData(folder::String; k::Int64=0, use_ros_cam::Bool=false)
    idx_file = joinpath(folder, "index.txt")
    idx_data = readdlm(idx_file)
    n = size(idx_data, 1)

    id2rgb = Dict{Int64, String}()
    id2depth = Dict{Int64, String}()

    ids = Vector{Int64}(undef, n)
    timestamps::Vector{Int64} = Vector{Int64}(undef, n)
    time_factor = 1000000
    for i=1:n
        id = idx_data[i, 1]
        ids[i] = id
        timestamps[i] = Int64(round(idx_data[i, 2] * time_factor))
        id2rgb[id] = joinpath(folder, idx_data[i, 3])
        id2depth[id] = joinpath(folder, idx_data[i, 4])
    end
    # ids start from 0. timestamp starts from 1
    cam = if use_ros_cam
        CameraModel.TumCamROS()
    elseif occursin("freiburg1", folder)

        fx = 517.306408

        # IR projector baseline times fx (aprox.)
        bf = 40.0

        CameraModel.RgbdCamParams(
            width=640,
            height=480,
            fx = fx,
            fy = 516.469215,
            cx = 318.643040,
            cy = 255.313989,
            depth_factor=5000.0,
            baseline=bf/fx,
            k1 = 0.262383,
            k2 = -0.953104,
            p1 = -0.005358,
            p2 = 0.002628,
            k3 = 1.163314
        )

    elseif occursin("freiburg2", folder)

        fx = 520.908620

        # IR projector baseline times fx (aprox.)
        bf = 40.0

        CameraModel.RgbdCamParams(
            width=640,
            height=480,
            fx = fx,
            fy = 521.007327,
            cx = 325.141442,
            cy = 249.701764,
            depth_factor=5208.0,
            baseline=bf/fx,
            k1 = 0.231222,
            k2 = -0.784899,
            p1 = -0.003257,
            p2 = -0.000105,
            k3 = 0.917205,
        )

    elseif occursin("freiburg3", folder)

        fx=535.4

        # IR projector baseline times fx (aprox.)
        bf = 40.0

        CameraModel.RgbdCamParams(
            width=640,
            height=480,
            fx = fx,
            fy=539.2,
            cx=320.1,
            cy=247.6,
            depth_factor=5000.0,
            baseline=bf/fx,
            k1=0.0,
            k2=0.0,
            p1=0.0,
            p2=0.0,
        )
    else
        @warn "no camera type detected"
    end

    return TumData(k, n, folder, ids, timestamps, id2rgb, id2depth, cam)
end

function read(tum_data::TumData)
    tum_data.k += 1
    return readk(tum_data, tum_data.k)
    # id = tum_data.ids[tum_data.k]
    # img = load(tum_data.id2rgb[id])
    # depth = load(tum_data.id2depth[id])
    # depth_view = channelview(depth)
    # depth_uint16 = copy(reinterpret(UInt16, depth_view[1,:,:]))
    #
    # return img, depth_uint16
end

function readk(tum_data::TumData, k::Int64)

    id = tum_data.ids[k]
    return readid(tum_data, id)
    # img = load(tum_data.id2rgb[id])
    # depth = load(tum_data.id2depth[id])
    # depth_view = channelview(depth)
    # depth_uint16 = copy(reinterpret(UInt16, depth_view[1,:,:]))
    # return img, depth_uint16
end

function readlabelid(tum_data::TumData, id::Int64)::Matrix{Int32}
    rgb_fname = tum_data.id2rgb[id]
    label_fname = "$(splitext(rgb_fname)[1]).npy"
    return npzread(label_fname)
end

function loaddepth(fname::String)
    depth = load(fname)
    depth_view = channelview(depth)
    depth_uint16 = copy(reinterpret(UInt16, depth_view[1,:,:]))
    return depth_uint16
end

function readid(tum_data::TumData, id::Int64)
    # id = tum_data.ids[k]
    img = FileIO.load(tum_data.id2rgb[id])
    depth = FileIO.load(tum_data.id2depth[id])
    depth_view = Images.channelview(depth)
    # depth_uint16 = copy(reinterpret(UInt16, depth_view[1,:,:]))
    depth_uint16 = if ndims(depth_view) ==  3
        copy(reinterpret(UInt16, depth_view[1,:,:]))
    else
        copy(reinterpret(UInt16, depth_view))
    end
    return img, depth_uint16
end

"""
readRGBD(reader::Sun3DReader; id::Int64)
read rgbd images with id as input
"""
readRGBD(tum_data::TumData; id::Int64) = readid(tum_data, id)

function fastforward!(tum_data::TumData, Δk::Int64)
    tum_data.k += Δk
end

function fastbackward!(tum_data::TumData, Δk::Int64)
    tum_data.k -= Δk
end

function id2time(tum_data::TumData, id::Int64)
    idx = get_id_index(tum_data, id)
    return tum_data.timestamps[idx]
end

"""
mutable struct Trajectory
    n::Int64
    fixed_id::Int64
    ids::Vector{Int64}
    poses::Dict{Int64, Matrix{Float64}}
end
tum data use float number as time, we apply a time factor to make it an integer
x_idx: the index for the x colum with 1 start index
"""
function load_tum_traj(f_path::String; time_factor::Int64=1000000, x_idx=2)
    # ids = Vector{Int64}(undef, 0)
    id2poses = Dict{Int64, Matrix{Float64}}()
    ids = Vector{Int64}(undef, 0)
    open(f_path) do file
        # do stuff with the open file
        while !eof(file)
            aline = readline(file)
            split_strs = split(aline, (' ', '\t'))
            if split_strs[1] == "#"
                # skip comment
                continue
            end
            # @show split_strs
            raw_time = parse(Float64, split_strs[1])
            id = Int64(round(raw_time * time_factor))#parse(Int64, split_strs[1])
            xyzqxqyqzqw = [parse(Float64, x) for x = split_strs[x_idx:x_idx+6]]
            qx, qy, qz, qw = xyzqxqyqzqw[4:7]
            # convert to transformation matrix
            R = Rotations.Quat(qw, qx, qy, qz)
            T = Matrix{Float64}(LinearAlgebra.I, 4, 4)
            T[1:3, 1:3] = R
            T[1:3, 4] = xyzqxqyqzqw[1:3]
            # display(T)

            # @show id, num2, num3
            # @show mat4
            # @show id, xyzqxqyqzqw
            # save res
            push!(ids, id)
            id2poses[id] = T
        end
    end
    n = length(ids)
    fixed_id = minimum(ids)
    is_pose_Twc = true
    traj = Trajectory(n, fixed_id, ids, id2poses, is_pose_Twc)
    return traj
end

"""
The format of each line is 'timestamp tx ty tz qx qy qz qw'
https://vision.in.tum.de/data/datasets/rgbd-dataset/file_formats
timestamp is loaded as us in Int64 type
"""
function save_tum_traj(
    f_path::String,
    traj_dict::Dict{Int64, Matrix{FT}};
    id2time::Dict{Int64, Int64}=Dict{Int64, Int64}(),
    time_factor::Int64=1000000,
    do_inv::Bool=false
) where FT<:Number

    traj_keys = sort!(collect(keys(traj_dict)))
    open(f_path, write=true) do file
        for traj_key = traj_keys
            T = traj_dict[traj_key]

            time_stamp = if length(id2time) > 0
                id2time[traj_key]
            else
                traj_key
            end
            T = if do_inv
                inv(T)
            else
                T
            end

            tx, ty, tz, = T[1:3, 4]
            q = Rotations.Quat(T[1:3, 1:3])
            qx = q.x
            qy = q.y
            qz = q.z
            qw = q.w
            out_string = if time_factor == 1
                out_string = Printf.@sprintf(
                    "%d %f %f %f %f %f %f %f",
                    time_stamp,
                    tx,
                    ty,
                    tz,
                    qx,
                    qy,
                    qz,
                    qw
                )
            else
                out_string = Printf.@sprintf(
                    "%.6f %f %f %f %f %f %f %f",
                    time_stamp/time_factor,
                    tx,
                    ty,
                    tz,
                    qx,
                    qy,
                    qz,
                    qw
                )
                # out_string = "$(time_stamp/time_factor) $tx $ty $tz $qx $qy $qz $qw"
            end
            println(file, out_string)
        end

    end

end

"""
bf index is in a format similar to TUM index
"""
function load_boundlefusion_index(
    idx_fpath::String;
    seq_folder::String="",
    only_keep_filename::Bool=true
)
    # idx_file = joinpath(folder, "index.txt")
    idx_data = readdlm(idx_fpath)
    n = size(idx_data, 1)

    id2rgb = Dict{Int64, String}()
    id2depth = Dict{Int64, String}()

    ids = Vector{Int64}(undef, n)
    # timestamps::Vector{Int64} = Vector{Int64}(undef, n)
    # time_factor = 1000000
    for i=1:n
        id = idx_data[i, 1]
        ids[i] = id
        # timestamps[i] = Int64(round(idx_data[i, 2] * time_factor))
        rgb_fname, depth_fname = if only_keep_filename
            splitpath(String(idx_data[i, 2]))[end], splitpath(String(idx_data[i, 3]))[end]
        else
            idx_data[i, 2], idx_data[i, 3]
        end
        id2rgb[id] = joinpath(seq_folder, rgb_fname)
        id2depth[id] = joinpath(seq_folder, depth_fname)
    end
    return ids, id2rgb, id2depth
end


text:
using Copulas
using Documenter

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

makedocs(;
    modules=[Copulas],
    authors="Oskar Laverny",
    repo="https://github.com/lrnv/Copulas.jl/blob/{commit}{path}#{line}",
    sitename="Copulas.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://lrnv.github.io/Copulas.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
        "Elliptical Copulas" => "ellipticals.md",
        "Archimedean Copulas" => "archimedeans.md"
    ],
)

deploydocs(;
    repo="github.com/lrnv/Copulas.jl",
    devbranch="main",
)


text:
# test interoperation between LaguerrePolynomial and LaguerreFunction.
# in particular, a product of LaguerrePolynomial and LaguerreFunction
# can be expanded into a linear combination of LaguerreFunctions.
@testset "Interoperation" begin
    @testset "Test throw" begin
        l1 = LaguerrePolynomial(2, 2.0)
        l2 = LaguerreFunction(3, 1.0)
        @test_throws ErrorException l1 + l2
        @test_throws ErrorException l1 - l2
    end
end


text:
mutable struct Person
	id::Int
	age::Int
	occupation::Int

	movement_state::Int
	work_school_place::Int
	transport::Int
	time_spent_out::Int

	function Person(id::Int, age::Int, occupation::Int, work_school_place::Int, transport::Int)
		new(id,age,occupation,1,work_school_place,transport,0)
	end
end

mutable struct Household
	id::Int
	workers::Array{Person,1}
	students::Array{Person,1}
	unemployeds::Array{Person,1}
	num_components::Int

	function Household(id::Int,num_components::Int)
		new(id,[],[],[],num_components)
	end

end

mutable struct Company
	id::Int
	max_size::Int
	current_size::Int

	function Company(id::Int,max_size::Int)
		new(id,max_size,0)
	end
end

mutable struct Contact_simulation_options
	work_place_num::Int
	school_place_num::Int
	transport_place_num::Int
	leisure_place_num::Int
	days::Int
end

max_stars_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_issues_repo_name: Luigi-Crisci/Hypergraphs-Epidemics

max_forks_repo_name: Luigi-Crisci/Hypergraphs-Epidemics


text:
#
# This file is a part of MolecularGraph.jl
# Licensed under the MIT License http://opensource.org/licenses/MIT
#

export
    descendants, ancestors,
    topologicalsort

descendants(graph::DirectedGraph, idx::Int) = reachablenodes(graph, idx)
ancestors(graph::DirectedGraph, idx::Int) = reversereachablenodes(graph, idx)

function topologicalsort(graph::DirectedGraph)
    result = Int[]
    stack = [i for i in nodeset(graph) if indegree(graph, i) == 0]
    used = Set{Int}()
    while !isempty(stack)
        n = pop!(stack)
        push!(result, n)
        for (ninc, nadj) in outneighbors(graph, n)
            ninc in used && continue
            push!(used, ninc)
            if isempty(setdiff(in_incidences(graph, nadj), used))
                push!(stack, nadj)
            end
        end
    end
    edgecount(graph) == length(used) || throw(ErrorException("cycle found"))
    return result
end


text:
module Counter

using Redux

# actions
struct IncrementAction <: AbstractSyncAction end
struct DecrementAction <: AbstractSyncAction end

const INCREMENT = IncrementAction()
const DECREMENT = DecrementAction()

# state
struct State <: AbstractImmutableState
    counter::Int
end

# reducers
counter(state::AbstractState, action::AbstractAction) = state
counter(state::Vector{<:AbstractState}, action::AbstractAction) = state
counter(state::State, action::IncrementAction) = State(state.counter + 1)
counter(state::State, action::DecrementAction) = State(state.counter - 1)

end # module


text:
using ITensorsBenchmarks

# Run all of the block sparse benchmarks using all of the available Julia threads and write the results
# By default, a verbose output is written to a log file
runbenchmarks(write_results = true, blocksparse_num_threads = Threads.nthreads())


text:
function sum_of_sqrt1(xs)
    total = zero(eltype(xs))
    for i in eachindex(xs)
        total += sqrt(xs[i])
    end
    return total
end

function sum_of_sqrt2(xs)
    total = zero(eltype(xs))
    for i in eachindex(xs)
        try
            total += sqrt(xs[i])
        catch
            # ignore error intentionally
        end
    end
    return total
end

using BenchmarkTools
x = rand(100_000);
@btime sum_of_sqrt1($x);
@btime sum_of_sqrt2($x);

#=
julia> x = rand(100_000);

julia> @btime sum_of_sqrt1($x);
  482.818 μs (0 allocations: 0 bytes)

julia> @btime sum_of_sqrt2($x);
  2.455 ms (0 allocations: 0 bytes)
=#

function sum_of_sqrt3(xs)
    total = zero(eltype(xs))
    for i in eachindex(xs)
        if xs[i] >= 0.0
            total += sqrt(xs[i])
        end
    end
    return total
end

@btime sum_of_sqrt3($x);
#=
julia> @btime sum_of_sqrt3($x);
  545.789 μs (2 allocations: 781.33 KiB)
=#


text:
import Base: size, iterate

export Point, Vec2D

struct Point{T <: Number}
  x::T
  y::T
end

struct Vec2D{T <: Number}
  dx::T
  dy::T
end

include("point.jl")
include("vec2d.jl")
include("rect.jl")

# Allows you to do things such as:
#  x, y, w, h = Rect(2, 4, 20, 30)
function iterate(t::Union{Point, Vec2D, Size, Rect}, i::Int=1)
    if 1 <= i <= fieldcount(typeof(t))
        (getfield(t, i), i+1)
    else
        nothing
    end
end


text:
using JetPack, Jets, Test, FFTW

n1,n2 = 101,201

@testset "JopRemoveDC, correctness test, T=$T" for T in (Float32,Float64,Complex{Float32},Complex{Float64})
    A = JopRemoveDC(JetSpace(T,n1,n2))
    x = -1 .+ 2 .* rand(domain(A)) .+ rand(T)
    y2 = A * x
    X = fft(x,(1,))
    X[1,:] .= 0
    y1 = ifft(X,(1,))
    @test y1 ≈ y2
end

@testset "JopRemoveDC, linearity test, T=$T" for T in (Float32,Float64,Complex{Float32},Complex{Float64})
    A = JopRemoveDC(JetSpace(T,n1,n2))
    lhs,rhs = linearity_test(A)
    @test lhs ≈ rhs
end

@testset "JopRemoveDC, dot product test, T=$T" for T in (Float32,Float64,Complex{Float32},Complex{Float64})
    #=
    I observe the following bug on windows:
    x = rand(Complex{Float64},101,201)
    using LinearAlgebra
    dot(x,x) # crashes Julia, presumably due to an issue with BLAS on windows.
    =#
    if Sys.iswindows() && T == Complex{Float64}
        @warn "dot product test fails on Windows due to the LinearAlgebra.dot method not working for Complex{Float64}"
    else
        A = JopRemoveDC(JetSpace(T,n1,n2))
        lhs,rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
        @test lhs ≈ rhs
    end
end


text:
import .ExpressionExplorer
import .ExpressionExplorer: join_funcname_parts, SymbolsState, FunctionNameSignaturePair

"Return a copy of `old_topology`, but with recomputed results from `cells` taken into account."
function updated_topology(old_topology::NotebookTopology, notebook::Notebook, cells)

	updated_codes = Dict{Cell,ExprAnalysisCache}()
	updated_nodes = Dict{Cell,ReactiveNode}()
	unresolved_cells = copy(old_topology.unresolved_cells)
	for cell in cells
		old_code = old_topology.codes[cell]
		if old_code.code !== cell.code
			new_code = updated_codes[cell] = ExprAnalysisCache(notebook, cell)
			new_symstate = new_code.parsedcode |>
				ExpressionExplorer.try_compute_symbolreferences
			new_reactive_node = ReactiveNode(new_symstate)

			updated_nodes[cell] = new_reactive_node
		elseif old_code.forced_expr_id !== nothing
			# reset computer code
			updated_codes[cell] = ExprAnalysisCache(old_code; forced_expr_id=nothing, function_wrapped=false)
		end

		new_reactive_node = get(updated_nodes, cell, old_topology.nodes[cell])
		if !isempty(new_reactive_node.macrocalls)
			# The unresolved cells are the cells for wich we cannot create
			# a ReactiveNode yet, because they contains macrocalls.
			push!(unresolved_cells, cell)
		else
			pop!(unresolved_cells, cell, nothing)
		end
	end
	new_codes = merge(old_topology.codes, updated_codes)
	new_nodes = merge(old_topology.nodes, updated_nodes)

	for removed_cell in setdiff(keys(old_topology.nodes), notebook.cells)
		delete!(new_nodes, removed_cell)
		delete!(new_codes, removed_cell)
		delete!(unresolved_cells, removed_cell)
	end

	NotebookTopology(nodes=new_nodes, codes=new_codes, unresolved_cells=unresolved_cells)
end


text:
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\Materials.jl")
include("C:\\Users\\warner\\Desktop\\My Stuff\\Personal Learning\\Monte Carlo\\NautilusMC\\NautilusMC\\get_XS.jl")

mutable struct Particle
    type::String
    x::Vector{Float64}
    v::Vector{Float64}
    E::Float64
    cell::Cell
    weight::Float64
    alive::Bool
end

function transport(particle)
    material = particle.cell.material
    max_distance = cell.max_distance(x, v)
    # Get total cross sections
    nuclide_idxs = findall(ZA_T -> (ZA_T[1] in material.nuclides) & (ZA_T[2] == material.temperature), xs_nuclides)
    map_lo = grid_index(E)
    idxs_lo = grid_map[nuclide_idxs, map_lo]
    idxs_hi = grid_map[nuclide_idxs, map_lo + 1] .+ 1
    idxs = find_idx.(xs_data[nuclide_idxs], E, idxs_lo, idxs_hi)
    interp_fracs = interp_frac_log.(xs_data[nuclide_idxs], E, idxs)
    σ_T0 = σ_T_idx.(xs_data[nuclide_idxs], idxs)
    σ_T1 = σ_T_idx.(xs_data[nuclide_idxs], idxs.+1)
    Σ_Ts = material.atomic_dens.*exp.(log.(σ_T0) + log.(σ_T1./σ_T0).*interp_fracs) * 1.0E-24
    Σ_T = sum(Σ_T)
    # Move particle
    distance = -log(rand()) / Σ_T
    particle.x += particle.v * distance
    # Determine collision nuclide
    nuclide_probs = cumsum(Σ_Ts) / Σ_T
    collision_nuclide = material.nuclides[findlast(p -> p < nuclide_probs)]
    # Determine collision reaction
    println(collision_nuclide)
    # Do the collision
end

max_stars_repo_name: wamcghee/NautilusMC

max_issues_repo_name: wamcghee/NautilusMC

max_forks_repo_name: wamcghee/NautilusMC


text:
export Pairwise, pairwise, collision_waypoints_itr
export SteeringEdge

# Pairwise Iterator
struct Pairwise{I}
    itr::I
end
pairwise(itr) = Pairwise(itr)
Base.length(p::Pairwise) = length(p.itr) - 1
Base.eltype(::Type{Pairwise{I}}) where {I} = Tuple{eltype(I),eltype(I)}
@inline function Base.iterate(p::Pairwise, state=iterate(p.itr))
    state === nothing && return nothing
    first_item, itr_state = state
    state = iterate(p.itr, itr_state)
    state === nothing && return nothing
    (first_item, state[1]), state
end

# SteeringBVP Collision Checking
collision_waypoints_itr(f::DifferentialDynamics, x::State, controls) = waypoints_itr(f, x, controls, 8)
function is_free_edge(CC::CollisionChecker, bvp::SteeringBVP, x0::State, xf::State, controls::Nothing=nothing)
    is_free_edge(CC, bvp.dynamics, x0, bvp(x0, xf).controls)
end
function is_free_edge(CC::CollisionChecker, bvp::SteeringBVP, x0::State, xf::State, controls)
    is_free_edge(CC, bvp.dynamics, x0, controls)
end
function is_free_edge(CC::CollisionChecker, f::DifferentialDynamics, x0::State, controls)
    CC.edge_count[] += 1
    all(is_free_motion(CC, x1, x2) for (x1, x2) in pairwise(collision_waypoints_itr(f, x0, controls)))
end

# Plotting
struct SteeringEdge{BVP,S,S2C,U}
    state2config::S2C
    bvp::BVP
    x0::S
    xf::S
    controls::U
end
SteeringEdge(state2config, bvp, x0, xf) = SteeringEdge(state2config, bvp, x0, xf, bvp(x0, xf).controls)
@recipe function f(E::SteeringEdge; dims=(1, 2), edge_waypoints=10, plot_endpoints=true) # , edge_color=:grey, edge_alpha=1, edge_markershape=:none
    # color :=  edge_color
    # alpha :=  edge_alpha
    # label --> ""
    x, y = dims
    plot_endpoints && @series begin
        seriestype  := :scatter
        # markershape := edge_markershape
        q0, qf = E.state2config(E.x0), E.state2config(E.xf)
        SVector(q0[x], qf[x]), SVector(q0[y], qf[y])
    end
    @series begin
        pts = waypoints(E.bvp.dynamics, E.x0, E.controls, edge_waypoints)
        [p[x] for p in pts], [p[y] for p in pts]
    end
end

max_stars_repo_name: pbpf/clean

max_issues_repo_name: pbpf/clean

max_forks_repo_name: pbpf/clean


text:
tol = 1e-13

pmax = 4

### Lagrangian basis functions
println( "Lagrangian basis functions")
for ii = 1:pmax

  ploc,   = luteos.genlocal(ii)
  ploc1D = fill( 0.0, ii+1, 1 )
  for jj = 1:ii+1
    ploc1D[jj] = (jj-1)/ii
  end
  temp = ploc1D[2:ii]
  ploc1D[1:2] = ploc1D[ [1,ii+1] ]
  ploc1D[3:ii+1] = temp

  master = Master2D( ii; typeb = "lag" )

  phi1D, dphi1D = luteos.basisFuncLineLag( Val{ii}, ploc1D )
  sz    = size( phi1D )
  for jj = 1:sz[1], kk = 1:sz[2]
    if jj == kk
      @test abs(phi1D[jj,kk] - 1.0) < tol
    else
      @test abs(phi1D[jj,kk]) < tol
    end
  end

  phi, dphi = luteos.basisFuncTriangleLag( Val{ii}, ploc[:,2], ploc[:,3] )

  sz  = size( phi )
  for jj = 1:sz[1], kk = 1:sz[2]
    if jj == kk
      @test abs(phi[jj,kk] - 1.0) < tol
    else
      @test abs(phi[jj,kk]) < tol
    end
  end

end

# test interpolation

for ii = 1:pmax

  (ploc,) = luteos.genlocal( ii )

  (ϕ,∇ϕ) = luteos.basisFuncTriangleLag( Val{ii}, ploc[:,2], ploc[:,3] )

  sz  = size( ϕ )
  for jj = 1:sz[1], kk = 1:sz[2]
    if jj == kk
      @test abs(ϕ[jj,kk] - 1.0) < tol
    else
      @test abs(ϕ[jj,kk]) < tol
    end
  end

  # check evaluation for random point
  xloc   = [ 0.11698  0.2678;
             0.06905  0.1376 ]
  (ϕ,∇ϕ) = luteos.basisFuncTriangleLag( Val{ii}, xloc[:,1], xloc[:,2] )
  for jj in 1:2, dd in 1:2
    xloctest = ϕ[:,jj]' * ploc[:,dd+1]
    @test abs( xloctest - xloc[jj,dd] ) < tol

  end

end

# ### Legendre basis functions
# println( "Legendre basis functions")
# for ii = 1:pmax
#
#   @printf( "  Test 2D master element at %i\n", ii )
#
#   master = Master2D( ii; typeb = "leg" )
#
#   println( "      1D")
#
#   res1D = master.ϕ1d' * diagm(master.gwts1d) * master.ϕ1d
#   sz    = size( res1D )
#   for jj = 1:sz[1], kk = 1:sz[2]
#     if jj == kk
#       @test abs(res1D[jj,kk] - 1.0) < tol
#     else
#       @test abs(res1D[jj,kk]) < tol
#     end
#   end
#
#   println( "      2D")
#
#   res = master.ϕ' * diagm(master.gwts) * master.ϕ
#   sz  = size( res )
#   for jj = 1:sz[1], kk = 1:sz[2]
#     if jj == kk
#       @test abs(res[jj,kk] - 1.0) < tol
#     else
#       @test abs(res[jj,kk]) < tol
#     end
#   end
#
# end


text:
BigInt(10)^1000   # 大整数
BigFloat(10)^1000 # 大浮点数，如需改变默认精度，可参考文档
big(1.5)          # 对应大数类型的值，在此例是 BigFloat
1 + 1im           # 复数
123//456          # 使用 // 算符写分数

max_stars_repo_name: zhaiyusci/Julia-Notes

max_issues_repo_name: zhaiyusci/Julia-Notes

max_forks_repo_name: zhaiyusci/Julia-Notes


text:
## This file contains mostly meval, infseval.

# Choose infinite or single evaluation.
# The test suite assumes infseval is used.
#doeval(x) = infseval(x)  # infinite evaluation

## export this for users
"""
    symeval(expr::Any)

send `expr` through the Symata evaluation sequence. `expr` is
an `Mxpr` or number, symbol, etc.

In particular, an `Expr` (i.e. not translated to Symata) Symata-evaluated (whic in this case means returned unchanged) by the Symata
evaluation sequence.
"""
symeval(args...) = doeval(args...)

"""
    doeval(expr::Any)

is the main entry point to the Symata evaluation sequence.
"""
function doeval(x)
    if is_throw()
        return x
    end
    return infseval(x)
end

#@inline doeval(x) = meval(x)   # single evaluation

# Enable or disable hashing expressions here.
# We hash expressions so that only unique copies are stored.
# This slows code. It could be faster in some circumstances.
# Eg. if we continue to store the dependent variables of an
# expression as metadata, then we don't need to regenerate them,
# if we find the cached copy of an expression.
#
lcheckhash(x) = x      # do nothing

# We should generate our own  (why did I write this ?)

"""
    get_localized_symbol(s::Symbol)

return a new symbol, like gensym. But, also set `Attribute` `Temporary` for
the new symbol.
"""
function get_localized_symbol(s::Symbol)
    increment_localization_counter()
    localized_symbol = Symbol(s, "!", get_localization_counter())
    set_attribute(localized_symbol, :Temporary)
    return localized_symbol
end

@mkapprule ClearTemporary  nargs => 0

@sjdoc ClearTemporary """
    ClearTemporary()

remove temporary symbols, ie all beginnig with "##", from the symbol table.
"""
# The Temporary attribute is not working. the symbols that escape are just gensysms
@doap function ClearTemporary()
    syms = usersymbols()
    for sym in syms
        ss = string(sym)
#        if length(ss) > 2 && ss[1:2] == "##"
        if  occursin(r"^##", ss)
            delete_sym(Symbol(sym))
        end
    end
    Null
end

@sjdoc Evaluate """
    Evaluate(expr)

specify that `expr` should be evaluated even if it appears as an argument protected by `HoldAll`, etc.

Example:

```
expr = 1 + x^2
f = Compile(Evaluate(expr))
```

!!! note
    `Evaluate` is only implemented for `HoldAll`; not `HoldRest`, etc.
"""

## Mma allows 0 or more args
@mkapprule Evaluate  nodefault => true

# This only handles trivial cases, where Evaluate is not protected by HoldXXX. The cases in which it
# *is* protected are handled below in meval_arguments()
@doap Evaluate(x) = x
@doap Evaluate(args...) = mxpr(:Sequence,args...)
@doap Evaluate() = mxpr(:Sequence)

# @mkapprule Unevaluated
# @doap function Unevaluated(arg)
#     arg
# end

## Macro for translation and evaluation, at repl or from file

# Read a line of user input, translate Expr to Mxpr, but don't evaluate result
macro exnoeval(ex)
    mx = extomx(ex)
    :(($(esc(mx))))
end

# This is called from the REPL
macro ex(ex)   # use this macro from the julia prompt
    mx = symataevaluate(ex)
    :(($(esc(mx))))
end

# The same macro, but we export it for users

"""
    @sym expr

Embed symata expression `expr` in Julia code.

Read and evaluate `expr` embedded in Julia code.

```
julia> a = 1       # Julia symbol `a`
julia> @sym a = 2  # Symata symbol
julia> a
1
julia> @sym a
2
```
"""
macro sym(ex)   # use this macro from the julia prompt
    mx = _sym_inner(ex, EvaluateJuliaSyntaxSimple())
    :(($(esc(mx))))
end

macro symfull(ex) # use this macro for Symata prompt.
    mx = _sym_inner(ex)
    :(($(esc(mx))))
end

function _sym_inner(ex, options=EvaluateJuliaSyntax())
    mx = symataevaluate(ex, options)
    if isa(mx, Symbol)
        mx = QuoteNode(mx)
    end
    return mx
end

# Simply evaluate. Do not print "Out", or count things, etc.
"""
    symtranseval(expr::Any)

Translate `expr` from `Expr` to `Mxpr` and send to the top level of the Symata evaluation sequence. `expr`
may also be a number, symbol, string, etc.
"""
function symtranseval(expr)
    symataevaluate(expr, EvaluateJuliaSyntaxSimple())
end

"""
    symparseeval(s::String)

Parse `s` into Julia expressions, translates them to Symata expressions, and Symata-evaluates each one,
returning the value returned by the final evaluation.

See `symparsestring`, `symtranseval`, and `@sym`.
"""
function symparseeval(s::String)
    mxprs = symparsestring(s)
    local res
    for mx in mxprs
        res = doeval(mx)
    end
    res
end

"""
    macro exsimple(ex::Expr)

translates `ex` to Symata expression and evaluates the result. Most code for interactive sessions disabled.
`@exsimple` is more conventient than `@sym` at the Julia prompt.
"""
macro exsimple(ex)   # use this macro from the julia prompt
    mx = symataevaluate(ex, EvaluateJuliaSyntaxSimple())
    :(($(esc(mx))))
end

"""
    debugmxpr(s::String)

parse and evaluate `s` as a Symata expression. From Julia, this prints the stack trace when an error is raised.
If the same expression entered at the Symata command line, the stack trace will not be printed.
"""
function debugmxpr(s::String)
    symataevaluate(Meta.parse(s))
end

"""
    number_of_Os

is the number of previous outputs to bind to `O`, `OO`, etc.
"""
const number_of_Os = 10

const Os = Array{SJSym}(undef, 0)
for i in 1:number_of_Os
    push!(Os, Symbol("O"^i))
end

# Put this in the System symbol table
# so that they are imported into Global
for i in 1:number_of_Os
    set_system_symval(Os[i], Null)
end

# Bind O to O(n), where n is the most recent line number.
# Bind OO to O(n-1), etc.  When O evaluates to Out(n), the Out rule evaluate
# the expression bound to Out(n).
macro bind_Os()
    expr = Expr(:block)
    for i in 1:number_of_Os
        sym =  string(Os[i])
        newex = :(
                  if (length(Output) - $i + 1) >= 1
                       oexp = mxpr(:Out, get_line_number() - $i + 1)
                       set_system_symval(Meta.parse($sym), oexp)
                       set_sysattributes($sym)
                  end
                     )
        push!(expr.args, newex)
    end
    expr
end

abstract type AbstractEvaluateOptions end

mutable struct EvaluateJuliaSyntax <: AbstractEvaluateOptions
end

function maybe_print_Out_label(opt::EvaluateJuliaSyntax)
    # do_we_print_Out_label is for IJulia
    if (! simple(opt)) && isinteractive() && do_we_print_Out_label
        print("Out(" * string(get_line_number()) * ") = ")
    end
    nothing
end

simple(opt::EvaluateJuliaSyntax) = false

mutable struct EvaluateJuliaSyntaxSimple <: AbstractEvaluateOptions
end

simple(opt::EvaluateJuliaSyntaxSimple) = true
maybe_print_Out_label(opt::EvaluateJuliaSyntaxSimple) = nothing

# """
#     NullExFuncOptions

# options for `symataevaluate` that cut out most of the code for interactive sessions.
# """
# const SimpleExFuncOptions = ExFuncOptions(true)

# """
#     NullExFuncOptions

# options for `symataevaluate` that do not modify its behavior
# """
# const NullExFuncOptions = ExFuncOptions(false)

"""
    symataevaluate(ex::Any, options)

Translate `ex` to Symata and evaluate the result. `ex` is typically an `Expr`,
or else an "elementary" type, such as a `Number` or `String`. If `ex` is already of type `Mxpr`, then
the translation is the identity. `ex` may be an object of any type.
"""
function symataevaluate(ex, options=EvaluateJuliaSyntax())
    if ! simple(options)
        check_doc_query(ex) && return nothing  # Asking for doc? Currently this is:  ?, SomeHead
    end
    res = extomx(ex)  # Translate to Mxpr
    local mx
    if ! simple(options)
        reset_meval_count()
        reset_try_downvalue_count()
        reset_try_upvalue_count()
        if is_timing()  #  && is_sjinteractive()  #  Disallow when reading  code. For now.
            @time mx = trysymataevaluate(res)
            println("tryrule count: downvalue ", get_try_downvalue_count(),", upvalue ", get_try_upvalue_count())
        else
            mx = trysymataevaluate(res)
        end
    else
        mx = trysymataevaluate(res)
    end
    symval(mx) == Null && return nothing
    maybe_print_Out_label(options)
    if (! simple(options) ) && isinteractive()  # && is_sjinteractive()
        increment_line_number()
        set_system_symval(:ans,mx)  # Like Julia and matlab, not Mma
        set_symata_prompt(get_line_number())
        if getkerneloptions(:history_length) > 0 push_output(mx) end
        @bind_Os
    end
    if is_throw()
        if is_Mxpr(mx,:Throw)
            @warn("Uncaught Throw")
            clear_throw()
        else
            @symwarn("Throw flag set, but expression is not throw.",  mx) # This should probably be an assert
            clear_throw()
            return mx
        end
    end
    return mx
end

function trysymataevaluate(mxin)
    try
        doeval(mxin)
    catch e
        if isa(e,ArgCheckErr)
            @warn(e.msg)
            return mxin
        elseif isa(e,RecursionLimitError)
            @warn(e.msg)
            e.mx
        else
            println(e)
            rethrow(e)  # It would be nice to get better error information.
        end
    end
end

#################################################################################
#                                                                               #
#  infseval                                                                     #
#  repeats meval till we reach a fixed point                                    #
#                                                                               #
#################################################################################

# We use infinite or fixed point evaluation: the Mxpr is evaled repeatedly until it does
# not change. Actually Mma, and Symata try to detect and avoid more evaluations.
# Also try to detect if the expression is simplified, (fixed or canonical).
# This is also complicated by infinite evaluation because whether an expression is
# simplified depends on the current environment. We try to solve this with lists of 'free' symbols.
# Note: lcheckhash is the identity (ie disabled)
# doeval is infseval: ie, we use 'infinite' evaluation. Evaluate till expression does not change.

# We could make this user settable somehow
recursion_limit() =  1024

# Diagnostic. Count number of exits from points in infseval
global const exitcounts = Int[0,0,0,0]

"""
    infseval(mxin::Mxpr)

Apply the evaluation `meval` to `mxin` repeatedly to either a fixed point,
or the recursion limit.
"""
function infseval(mxin::Mxpr)
    if is_throw()
        println("Caught throw at toplevel infseval() for Mxpr")
        return mxin  # TODO check where this test is actually used. We have it in three places
    end
    @mdebug(2, "infseval ", mxin)
    neval = 0
    if checkdirtysyms(mxin)                   # is timestamp on any free symbol in mxin more recent than timestamp on mxin ?
        unsetfixed(mxin)                      # Flag mxin as not being at its fixed point in this environment.
    end                                       # This might be good for iterating over list of args in Mxpr.
    is_fixed(mxin) && return lcheckhash(mxin) # If mxin was already fixed and none of its free vars changed, just return.
    mx = meval(mxin)                          # Do the first evaluation
    if is_Mxpr(mx)
        is_fixed(mx) && return lcheckhash(mx) # The first meval may have set the fixed point flag. Eg, an Mxpr with only numbers, does not need another eval.
        if mx == mxin       # meval did not set fixed flag, but we see that it is at fixed point.
            setfixed(mx)    # They may be equal but we need to set fixed bit in mx !
            setfixed(mxin)  # Do we need to do this to both ?
            return lcheckhash(mx)
        end
    end
    local mx1
    while true              # After 1 eval fixed bit is not set and input not equal to result of first eval
        mx1 = meval(mx)     # So, we do another eval.
        if (is_Mxpr(mx1) && is_fixed(mx1))  || mx1 == mx  # The most recent eval was enough, we are done
            mx = mx1
            break
        end
        neval += 1
        if neval > 10  # Defining Subdivide via Subdivide(a_,b_) := Subdivide(0,a,b), etc. causes neval == 3, which was too many
            setfixed(mx)
            throw(RecursionLimitError("infseval: Too many, $neval, evaluations. Expression still changing", mxprcf(:Hold,mx)))
        end
        mx = mx1
    end
    if is_Mxpr(mx) && mx == mxin
        setfixed(mxin)
        setfixed(mx)
    end
    return lcheckhash(mx)  # checking hash code is disbled.
end

function infseval(qs::Qsym)
    mx = meval(qs)
    return mx == qs ? qs : infseval(mx)
end

function infseval(s::SJSym)
    mx = meval(s)
    return mx == s ? s : infseval(mx)
end
# Any type that other than SJSym (ie Symbol) or Mxpr is not meval'd.
@inline infseval(x) = x
@inline infseval(x::Complex{T}) where {T<:Real} = x.im == zero(x.im) ? x.re : x

infseval(x::Complex{Rational{T}}) where {T<:Integer} = meval(x)
infseval(x::Rational{T}) where {T<:Integer} = meval(x)

meval(x::Float64) = x == Inf ? Infinity : x == -Inf ? MinusInfinity : x

#################################################################################
#                                                                               #
#  meval  Evaluation of Mxpr                                                    #
#  main evaluation routine. Call doeval on head                                 #
#  and some of the arguments. Then apply rules and other things on the result.  #
#                                                                               #
#################################################################################

# These are normally not called, but rather are caught by infseval.
@inline meval(x::Complex{T}) where {T<:Real} = x.im == 0 ? x.re : x
meval(x::Rational{T}) where {T<:Integer} = x.den == 1 ? x.num : x.den == 0 ? ComplexInfinity : x
meval(x::Nothing) = Null
meval(x) = x
meval(s::SJSym) = symval(s) # this is where var subst happens

function meval(qs::Qsym)
    res = symval(qs)
    res
end

function meval(mx::Mxpr)
    increment_meval_count()
    if get_meval_count() > recursion_limit()
        res = mxprcf(:Hold,mx)
        deepsetfixed(res)
        reset_meval_count()
        reset_try_downvalue_count()
        reset_try_upvalue_count()
        throw(RecursionLimitError("Recursion depth of " * string(recursion_limit()) *  " exceeded.", res))
    end
    local ind::String = ""  # some places get complaint that its not defined. other places no !?
    if is_meval_trace()
        ind = " " ^ (get_meval_count() - 1)
        println(ind,">>", get_meval_count(), " " , wrapout(mx))
    end
    nmx::Mxpr = meval_arguments(mx)
    @mdebug(2, "meval: done meval_args ", nmx)
    setfreesyms(nmx,revisesyms(mx)) # set free symbol list in nmx
    @mdebug(2, "meval: done setfreesyms ")
    res = meval_apply_all_rules(nmx)
    @mdebug(2, "meval: done meval_apply_all_rules ", res)
    is_meval_trace() && println(ind,get_meval_count(), "<< ", wrapout(res))
    decrement_meval_count()
    @mdebug(2, "meval: returning ", res)
    return res
end

# ?? why the first test  ! is_canon(nmx). This must apparently always be satisfied.
function meval_apply_all_rules(nmx::Mxpr)
    if  ! is_canon(nmx)
        if isFlat(nmx) nmx = flatten(nmx) end
        if isListable(nmx) nmx = threadlistable(nmx) end
        res = canonexpr!(nmx)
    end
    @mdebug(2, "meval_apply_all_rules: entering apprules: ", res)
    res = apprules(res)           # apply "builtin" rules
    @mdebug(2, "meval_apply_all_rules: exited apprules ", res)
    is_Mxpr(res) || return res
    @mdebug(2, "meval_apply_all_rules: entering ev_upvalues ", res)
    res = ev_upvalues(res)
    @mdebug(2, "meval_apply_all_rules: exited ev_upvalues ", res)
    res = ev_downvalues(res)
    @mdebug(2, "meval_apply_all_rules: entering merge_args_if_empty_syms")
    merge_args_if_emtpy_syms(res) # merge free symbol lists from arguments
    @mdebug(2, "meval_apply_all_rules: exiting")
    return res
end

@inline function argeval(arg)
    if isa(arg,Mxpr{:Unevaluated}) && length(arg) > 0
        ## We may want to strip Unevaluated here, to make it exactly like Mma
        ## But, then it will not prevent Sequence substitution below
        ## And it is not held is some other forms as well.
        ## So, these would need to be fixed.
        return sjcopy(arg)
    else
        return doeval(arg)
    end
end

# Evaluate arguments of mx, construct and return new Mxpr
function meval_arguments(mx::Mxpr)
    nhead = doeval(mhead(mx))
    local nargs::MxprArgs
    mxargs::MxprArgs = margs(mx)
    len::Int = length(mxargs)
    if len == 0
        return mxpr(nhead)
    elseif isHoldFirst(nhead)
        nargs = newargs(len)
        nargs[1] = mxargs[1]
        @inbounds for i in 2:length(mxargs)
            nargs[i] = argeval(mxargs[i])
        end
    elseif isHoldAll(nhead)
        nargs = copy(mxargs)
        for i=1:length(nargs)
            if isa(nargs[i],Mxpr{:Evaluate}) && length(nargs[i]) > 0
                nargs[i] = doeval(nargs[i][1])
            end
        end
    elseif isHoldAllComplete(nhead)
        nargs = copy(mxargs)
    elseif isHoldRest(nhead)
        nargs = copy(mxargs)
        nargs[1] = argeval(nargs[1])
    else      # Evaluate all arguments
        nargs = newargs(len)
        @inbounds for i in 1:len
            nargs[i] = argeval(mxargs[i])
        end
    end
    if  (! isSequenceHold(nhead))  &&  (! isHoldAllComplete(nhead))
        splice_sequences!(nargs)
        len = length(nargs)
    end
    for i=1:len
        if (isa(nargs[i],Mxpr{:Unevaluated}) && length(nargs[i]) > 0) &&
            ! (isHoldAll(nhead) || isHoldAllComplete(nhead))  && !( (isHoldFirst(nhead) && i == 1) || (isHoldRest(nhead) && i > 1))
            nargs[i] = nargs[i][1]
        end
    end
    nmx::Mxpr = mxpr(nhead,nargs)   # new expression with evaled args
end

## We do meval_arguments specially for List, in order to remove occurrances of Nothing
## TODO We probably need to worry about Unevaluated, etc. here.
function meval_arguments(mx::Mxpr{:List})
    nhead = :List
    local nargs::MxprArgs
    mxargs::MxprArgs = margs(mx)
    len::Int = length(mxargs)
    # maybe using this flag is efficient. 'Nothing' is relatively rare.
    got_nothing::Bool = false
    if len == 0
        return mxpr(:List)
    else
        nargs = newargs(len)
        @inbounds for i in 1:len
            nargs[i] = doeval(mxargs[i])
            if nargs[i] == :Nothing
                got_nothing = true
            end
        end
    end
    if got_nothing
        ninds = Array{Int}(undef, 0)
        for i in 1:len
            if nargs[i] == :Nothing
                push!(ninds,i)
            end
        end
        deleteat!(nargs,ninds)
    end
    splice_sequences!(nargs)
    nmx::Mxpr = mxpr(nhead,nargs)   # new expression with evaled args
end

# Similar to checkdirtysyms. The original input Mxpr had a list of free symbols.
# That input has been mevaled at least once and the result is mx, the argument
# to revisesyms. Here, we make a free-symbol list for mx. We look at its current
# free symbol list, which is inherited, and identify those that are no longer
# free. Eg. the environment changed. E.g The user set a = 1. Or 'a' may
# evaluate to an expression with other symbols.
#
# move this to mxpr_type ??
function revisesyms(mx::Mxpr)
    s::FreeSyms = getfreesyms(mx)
    mxage::UInt64 = getage(mx)
    nochange::Bool = true      # Don't create a new symbol list if nothing changed
    for sym in keys(s)         # Check if changes. Does this save or waste time ?
        if symage(sym) > mxage
            nochange = false
            break
        end
    end
    # Need to return a copy here, or Table(x^i + x*i + 1,[i,10]) shows a bug.
    nochange == true && return copy(s)
    nsyms::FreeSyms = newsymsdict()
    for sym in keys(s)
        if symage(sym) > mxage
            mergesyms(nsyms,symval(sym))
        else
            mergesyms(nsyms,sym)  # just copying from the old symbol list
        end
    end
    return nsyms
end

## Try applying downvalues

@inline function ev_downvalues(mx::Mxpr)
    if has_downvalues(mx)
        return applydownvalues(mx)
    else
        return mx
    end
end
@inline ev_downvalues(x) = x

## Applying upvalues. This has to be efficient, we must not iterate over args.
#  Instead, we check free-symbol list.

@inline function ev_upvalues(mx::Mxpr)
    merge_args_if_emtpy_syms(mx) # do upvalues are for free symbols in mx.
    for s in listsyms(mx)
        if has_upvalues(s)
            mx = applyupvalues(mx,s)
            break      # I think we are supposed to only apply one rule
        end
    end
    return mx
end
@inline ev_upvalues(x) = x

## Build list of free syms in Mxpr if list is empty.
@inline function merge_args_if_emtpy_syms(mx::Mxpr)
    if isempty(getfreesyms(mx))    # get free symbol lists from arguments
        mergeargs(mx)              # This is costly if it is not already done.
        add_nothing_if_no_syms(mx) # If there are no symbols, add :nothing, so this is not called again.
    end
end
@inline merge_args_if_emtpy_syms(x) = nothing

#### Thread Listable over Lists

# If any arguments to mx are lists, thread over them. If there are
# more than one list, they must be of the same length.  Eg.
# f([a,b,c],d) -> [f(a,d),f(b,d),f(c,d)].
#
# This is general, but is not the most efficient way. Some Specific
# cases, or classes should be handled separately. Eg Adding two large
# lists of numbers, is slow this way. The problem is that
# threadlistable is called early in the evaluation sequence, as per
# Mma. So we can't add lists of numbers at that point.

function threadlistable(mx::Mxpr)
    pos = Array{Int}(undef, 0)      # FIXME should avoid this
    lenmx = length(mx)
    lenlist::Int = -1
    h = mhead(mx)
    @inbounds for i in 1:lenmx
        if is_Mxpr(mx[i],:List)
            nlen = length(mx[i])
            if lenlist >= 0 && nlen != lenlist
                error("Can't thread over lists of different lengths.")
            end
            lenlist = nlen
            push!(pos,i)
        end
    end
    lenp = length(pos)
    lenp == 0 && return mx      # Nothing to do. return input array
    largs = newargs(lenlist)
    @inbounds for i in 1:lenlist
        nargs = newargs(lenmx)
        p = 1
        @inbounds for j in 1:lenmx
            if p <= lenp && pos[p] == j
                nargs[j] = mx[j][i]
                p += 1
            else
                nargs[j] = mx[j]
            end
        end
        largs[i] = mxpr(h,nargs)
    end
    nmx = mxpr(:List,largs)
    return nmx
end

#### Splice expressions with head Sequence into argument list

## f(a,b,Sequence(c,d),e,f) -> f(a,b,c,d,e,f)
## The following is broken. returns a non-symata object
## FIXME!: f(Sequence([i1,3],[i2,4],[i3,2],[i4,10]))
## No, the sequence is ok. something is broken with Out(n), or ...
## args are args of an Mxpr
function splice_sequences!(args)
    length(args) == 0 && return
    i::Int = 1
    while true
        if is_Mxpr(args[i], :Sequence)
            sargs = margs(args[i])
            splice!(args,i,sargs)
            i += length(sargs)   # skip over new arguments.
        end                      # splicing in lower levels should be done already.
        i += 1
        i > length(args) && break
    end
end


text:
module Positions

import JLD2
using Printf: @sprintf
using ..SimLib
using ..SimLib: Maybe, FArray
using XXZNumerics

export PositionDataDescriptor, PositionData, load_positions

## Data structure

"""
    struct PositionDataDescriptor

The important bits of information needed to specify a set positions are:
 - geometry
 - dimension
 - system_size
 - shots
 - densities ρs
For `load`ing data the last 2 may be omitted.
"""
struct PositionDataDescriptor <: SimLib.AbstractDataDescriptor
    geometry::Symbol
    dimension::Int
    system_size::Int
    shots::Maybe{Int}
    ρs::Maybe{FArray{1}}
    pathdata::SaveLocation
    function PositionDataDescriptor(geom, dimension, system_size, shots, ρs, pathdata::SaveLocation)
        geom ∈ SimLib.GEOMETRIES || error("Unknown geometry: $geom")
        new(geom, dimension, system_size, shots, ismissing(ρs) ? missing : unique!(sort(collect(ρs))), pathdata)
    end
end

function PositionDataDescriptor(geom, dimension, system_size, shots=missing, ρs=missing; prefix=path_prefix(), suffix="")
    PositionDataDescriptor(geom, dimension, system_size, shots, ρs, SaveLocation(prefix, suffix))
end

function Base.:(==)(d1::PositionDataDescriptor, d2::PositionDataDescriptor)
    all(getfield(d1, f) == getfield(d2, f) for f in [:geometry, :dimension, :system_size, :shots, :ρs])
end

SimLib._filename(desc::PositionDataDescriptor) = filename(desc.geometry, desc.dimension, desc.system_size)
filename(geometry, dimension, system_size) = @sprintf("positions/%s_%id_N_%02i", geometry, dimension, system_size)

load_positions(geometry, dimension, system_size, location=SaveLocation(); prefix=location.prefix, suffix=location.suffix) = load(PositionDataDescriptor(geometry, dimension, system_size; prefix, suffix))

"""
    struct PositionData

Stores the actual data for the positions specified by the descriptor [`PositionDataDescriptor`](@ref).
The indices mean:
 - coordinate index
 - particle
 - shot number
 - density ρ
Thus the data as dimensions: [D, N, #SHOT, #ρ]

The default save directory is "positions".
"""
struct PositionData <: SimLib.AbstractSimpleData
    descriptor::PositionDataDescriptor
    # [xyz, N, shot, ρ]
    data::Array{Float64, 4}
end

PositionData(desc::PositionDataDescriptor) = PositionData(desc, zeros(Float64, desc.dimension, desc.system_size, desc.shots, length(desc.ρs)))
PositionData(args...; kwargs...) = PositionData(PositionDataDescriptor(args..., kwargs...))

## main function
function create_positions!(empty_posdata; fail_rate=0.3)
    logmsg("Generating positions for $(empty_posdata.descriptor)")
    dim = empty_posdata.dimension
    N = empty_posdata.system_size
    shots = empty_posdata.shots

    max_misses = shots*fail_rate
    for (i, ρ) in enumerate(empty_posdata.ρs)
        sampler = geometry_from_density(empty_posdata.geometry, ρ, N, dim)
        misses = 0
        for j in 1:shots
            positions = Vector{Vector{Float64}}()
            while length(positions) == 0 && misses < max_misses
                try
                    positions = sample_blockaded(sampler, N)
                catch e;
                    misses += 1
                    if !(e isa ErrorException)
                        logmsg("Got unexpected Error: $e\nCurrent ρ=$ρ, dim=$dim, shot=$j, N=$N")
                    end
                end
            end
            if misses >= max_misses
                error("sampler $sampler did fail to converge too often (fail rate $(round(fail_rate*100; digits=1))%! Current ρ=$ρ, dim=$dim, shot=$j, N=$N")
            end
            if dim == 1
                # for spin chains, sorting the positions is very helpful!
                # it enables f.e. easier computation of half-chain entropy, tracing out subsystems...
                positions = sort!(positions)
            end
            empty_posdata[:,:,j,i] = hcat(positions...)
        end
    end
    empty_posdata # now filled
end

SimLib.create(desc::PositionDataDescriptor) = create_positions!(PositionData(desc))

end #module


text:
module TestExamplesKPP

using Test
using Trixi

include("test_trixi.jl")

EXAMPLES_DIR = joinpath(examples_dir(), "tree_2d_dgsem")

@testset "KPP" begin
  @trixi_testset "elixir_kpp.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_kpp.jl"),
      l2   = [0.36563290910786106],
      linf = [9.116732052340398],
      max_refinement_level = 6,
      tspan = (0.0, 0.01),
      atol = 1e-6,
      rtol = 1e-6,
      skip_coverage = true)
  end
end

end # module


text:
using LinearAlgebra
using LinearMaps
using ArnoldiMethod
using SphericalOperators
import SphericalOperators: ord

struct InverseMap{M} <: LinearMap{M}
    factorization
end
InverseMap(A::AbstractMatrix) = InverseMap{eltype(A)}(factorize(A))

Base.size(IM::InverseMap) = size(IM.factorization)

function LinearMaps.A_mul_B!(B, IM::InverseMap, A)
    B[:] = IM.factorization\A
end

function eigenstates(H, σ=-10.0; kwargs...)
    Hs = H - sparse(I,size(H)...)*σ
    Hf = InverseMap(Hs)

    ee, = partialschur(Hf; kwargs...)

    λ = 1 ./ diag(ee.R) .+ σ

    λ,ee.Q
end

function ground_state(H, args...; kwargs...)
    λ,ϕ = eigenstates(H, args...; kwargs...)
    real(λ[1]),real.(ϕ[:,1])
end

function ground_state(basis::FEDVR.Basis, L::AbstractSphericalBasis,
                      ℓ₀::Integer, m₀::Integer=0,
                      ::Type{O}=SphericalOperators.LexicalOrdering,
                      args...; kwargs...) where {O<:SphericalOperators.Ordering}
    ℓ₀ ∉ eachℓ(L) && error("Requested initial partial wave $(ℓ₀) not in range $(eachℓ(L))")
    m₀ ∉ eachm(L,ℓ₀) && error("Requested initial projection quantum number $(m₀) not in range $(eachm(L,ℓ₀))")
    Hℓ₀ = hamiltonian(basis, ℓ₀)

    nᵣ = basecount(basis.grid)
    M = prod(size(L))
    ψ₀ = zeros(M)
    ψ₀[ord(L,O,ℓ₀,m₀,1:nᵣ)] = ground_state(Hℓ₀, args...; kwargs...)[2]
    ψ₀
end

export eigenstates, ground_state


text:
"""
    struct ZeroMeanFESpace <: SingleFieldFESpace
      # private fields
    end
"""
struct ZeroMeanFESpace{CA,S} <: SingleFieldFESpace
  space::FESpaceWithConstantFixed{CA,S}
  vol_i::Vector{Float64}
  vol::Float64
end

"""
"""
function ZeroMeanFESpace(space::SingleFieldFESpace,dΩ::Measure)
  _space = FESpaceWithConstantFixed(space,true,num_free_dofs(space))
  vol_i = assemble_vector(v->∫(v)*dΩ,space)
  vol = sum(vol_i)
  ZeroMeanFESpace(_space,vol_i,vol)
end

# Genuine functions

function TrialFESpace(f::ZeroMeanFESpace)
  U = TrialFESpace(f.space)
  ZeroMeanFESpace(U,f.vol_i,f.vol)
end

function FEFunction(
  f::ZeroMeanFESpace,
  free_values::AbstractVector,
  dirichlet_values::AbstractVector)
  c = _compute_new_fixedval(
    free_values,
    dirichlet_values,
    f.vol_i,
    f.vol,
    f.space.dof_to_fix)
  fv = lazy_map(+,free_values,Fill(c,length(free_values)))
  dv = dirichlet_values .+ c
  FEFunction(f.space,fv,dv)
end

function EvaluationFunction(f::ZeroMeanFESpace,free_values)
  FEFunction(f.space,free_values)
end

function _compute_new_fixedval(fv,dv,vol_i,vol,fixed_dof)
  @assert length(fv) + 1 == length(vol_i)
  @assert length(dv) == 1
  c = zero(eltype(vol_i))
  for i=1:fixed_dof-1
    c += fv[i]*vol_i[i]
  end
  c += first(dv)*vol_i[fixed_dof]
  for i=fixed_dof+1:length(vol_i)
    c += fv[i-1]*vol_i[i]
  end
  c = -c/vol
  c
end

# Delegated functions

get_triangulation(f::ZeroMeanFESpace) = get_triangulation(f.space)

ConstraintStyle(::Type{ZeroMeanFESpace{CA,S}}) where {CA,S} = ConstraintStyle(S)

CellField(t::ZeroMeanFESpace,cell_vals) = CellField(t.space,cell_vals)

get_cell_isconstrained(f::ZeroMeanFESpace) = get_cell_isconstrained(f.space)

get_cell_constraints(f::ZeroMeanFESpace) = get_cell_constraints(f.space)

get_dirichlet_values(f::ZeroMeanFESpace) = get_dirichlet_values(f.space)

get_cell_shapefuns(f::ZeroMeanFESpace) = get_cell_shapefuns(f.space)

get_cell_shapefuns_trial(f::ZeroMeanFESpace) = get_cell_shapefuns_trial(f.space)

get_cell_dof_basis(f::ZeroMeanFESpace) = get_cell_dof_basis(f.space)

get_free_dof_ids(f::ZeroMeanFESpace) = get_free_dof_ids(f.space)

zero_free_values(f::ZeroMeanFESpace) = zero_free_values(f.space)

get_vector_type(f::ZeroMeanFESpace) = get_vector_type(f.space)

get_dof_value_type(f::ZeroMeanFESpace) = get_dof_value_type(f.space)

get_cell_dof_ids(f::ZeroMeanFESpace) = get_cell_dof_ids(f.space)

get_dirichlet_dof_ids(f::ZeroMeanFESpace) = get_dirichlet_dof_ids(f.space)

zero_dirichlet_values(f::ZeroMeanFESpace) = zero_dirichlet_values(f.space)

num_dirichlet_tags(f::ZeroMeanFESpace) = num_dirichlet_tags(f.space)

get_dirichlet_dof_tag(f::ZeroMeanFESpace) = get_dirichlet_dof_tag(f.space)

scatter_free_and_dirichlet_values(f::ZeroMeanFESpace,fv,dv) = scatter_free_and_dirichlet_values(f.space,fv,dv)

gather_free_and_dirichlet_values(f::ZeroMeanFESpace,cv) = gather_free_and_dirichlet_values(f.space,cv)

gather_free_and_dirichlet_values!(fv,dv,f::ZeroMeanFESpace,cv) = gather_free_and_dirichlet_values!(fv,dv,f.space,cv)

gather_dirichlet_values(f::ZeroMeanFESpace,cv) = gather_dirichlet_values(f.space,cv)

gather_dirichlet_values!(dv,f::ZeroMeanFESpace,cv) = gather_dirichlet_values!(dv,f.space,cv)

gather_free_values(f::ZeroMeanFESpace,cv) = gather_free_values(f.space,cv)

gather_free_values!(fv,f::ZeroMeanFESpace,cv) = gather_free_values!(fv,f.space,cv)


text:
# Julia program to implement Circle Sort

# Performs recursive circular swaps and returns true if atleast one swap occurs
function rec_sort(arr, high, low)
    isSwap = false

    # If concerned array is empty, Return False
    if (high == low)
        return false
    end

    # Storing the values of beg, end to later use in the recursive call
    highA = high
    lowA = low
    while( highA < lowA)
        if (arr[highA] > arr[lowA])
            arr[highA], arr[lowA] = arr[lowA], arr[highA]
            isSwap = true;
        end
        highA = highA + 1
        lowA = lowA - 1
    end

    # If the array has odd number of elements
    if (highA == lowA)
        if (arr[highA] > arr[lowA + 1])
            arr[high], arr[lowA+1] = arr[lowA+1], arr[high]
            isSwap = true
        end
    end

    mid = (low - high) ÷ 2
    isSwapA = rec_sort(arr, high, high+mid)
    isSwapB = rec_sort(arr, high+mid+1, low)

    return (isSwap || isSwapA || isSwapB)
end

function circle_sort(arr, n)
    while (rec_sort(arr, 1, n))
    end
end

print("How many numbers do you want to sort? ")
n = readline()
n = parse(Int64, n)
if (n <= 0)
    println("No numbers to sort!!!")
    exit()
end
arr = Int[]
print("Enter the numbers: ")
arr = [parse(Int, num) for num in split(readline())]
circle_sort(arr, n)
print("The numbers in sorted order is ")
print(arr)

"""
Time Complexity - O(n * log(n)), where 'n' is the size of the array
Space Complexity - O(n)

SAMPLE INPUT AND OUTPUT

SAMPLE I

How many numbers do you want to sort? 6
Enter the numbers: 76 243 89 25 8337 72
The numbers in sorted order is [25, 72, 76, 89, 243, 8337]

SAMPLE II

How many numbers do you want to sort? -3
No numbers to sort!!!
"""


text:
using PowerSystems
using NLsolve
const PSY = PowerSystems

############### Data Network ########################
include(joinpath(dirname(@__FILE__), "dynamic_test_data.jl"))
include(joinpath(dirname(@__FILE__), "data_utils.jl"))
############### Data Network ########################
threebus_file_dir = joinpath(dirname(@__FILE__), "ThreeBusNetwork.raw")
threebus_sys = System(
    PowerModelsData(threebus_file_dir),
    runchecks = false,
    unit_system = "device_base",
)
add_source_to_ref(threebus_sys)
#Reduce generator output
for g in get_components(Generator, threebus_sys)
    set_active_power!(g, 0.75)
end
res = solve_powerflow!(threebus_sys)

function dyn_gen_five_mass_shaft_order(generator)
    return PSY.DynamicGenerator(
        generator,
        1.0, # ω_ref,
        machine_oneDoneQ(), #machine
        shaft_fivemass(), #shaft
        avr_type1(), #avr
        tg_none(),
        pss_none(),
    ) #pss
end

function dyn_gen_first_order(generator)
    return PSY.DynamicGenerator(
        generator,
        1.0, # ω_ref,
        machine_oneDoneQ(), #machine
        shaft_damping(), #shaft
        avr_type1(), #avr
        tg_none(), #tg
        pss_none(),
    ) #pss
end

for g in get_components(Generator, threebus_sys)
    if get_number(get_bus(g)) == 103
        case_gen = dyn_gen_five_mass_shaft_order(g)
        add_component!(threebus_sys, case_gen)
    elseif get_number(get_bus(g)) == 102
        case_inv = dyn_gen_first_order(g)
        add_component!(threebus_sys, case_inv)
    end
end
#Compute Y_bus after fault
fault_branches = deepcopy(collect(get_components(Branch, threebus_sys))[2:end])
Ybus_fault = PSY.Ybus(fault_branches, get_components(Bus, threebus_sys))[:, :]


text:
using Ranges
using Test

@testset "positional range" begin
    for length in [2, 10, 100]
        @test range(length) == Base.OneTo(length)
        for start in [-1, 0, 1, 10]
            @test range(start, length) == start:start+length-1
            for stop in [-2, 0, 2, 100]
                @test range(start, stop, length) == Base.range(start, stop; length)
                for step in [-1, 1, 1]
                    @test_throws ArgumentError range(start, step, stop, length)
                    @test range(start, step, stop, nothing) == Base.range(start, stop; step)
                    @test range(start, step, nothing, length) == Base.range(start; step, length)
                    @test range(start, nothing, stop, length) == Base.range(start, stop; length)
                    @test (range(nothing, step, stop, length),) .|> (Base.step, last, Base.length) == (step, stop, length)
                end
            end
        end
    end
end

@testset "Pair range and PartialRange" begin
    for start in [-1, 0, 1, 10]
        for stop in [-2, 0, 2, 100]
            r = range(start => stop)
            @test isa(r,Ranges.PartialRange)
            for length in [2, 10, 100]
                @test range(start => stop, length) == Base.range(start, stop; length)
                @test range(length, start => stop) == Base.range(start, stop; length)
                @test r(length) == Base.range(start, stop; length)
                @test r(;length) == Base.range(start, stop; length)
            end
            for step in [-1, 1, 1]
                @test r(;step) == Base.range(start, stop; step)
            end
        end
    end
end

@testset "Keyword range" begin
    for length in [2, 10, 100]
        @test range(;length) == Base.OneTo(length)
        @test range(:bels;l=length, b=1) == Base.OneTo(length)
        for start in [-1, 0, 1, 10]
            @test range(;start, length) == start:start+length-1
            @test range(:bels;b=start, l=length) == start:start+length-1
            for stop in [-2, 0, 2, 100]
                @test range(;start, stop, length) == Base.range(start, stop; length)
                @test range(:bels;b=start, e=stop, l=length) == Base.range(start, stop; length)
            end
            for step in [-1, 1, 1]
                @test range(;start, length, step) == Base.range(start; step, length)
                @test range(:bels;b=start, l=length, s=step) == Base.range(start; step, length)
            end
        end
        for stop in [-2, 0, 2, 100]
            @test range(;length, stop) == Base.range(1; length, stop)
            @test range(:bels;l=length, e=stop, b=1) == Base.range(1; length, stop)
        end
        for step in [-1, 1, 1]
            @test range(;length, step) == Base.range(1; length, step)
            @test range(:bels;l=length, s=step, b=1) == Base.range(1; length, step)
        end
    end
end

@testset "Abbreviated keyword range" begin
    b = 1
    e = 5
    l = 5
    s = 1
    bels = [:b, :e, :l, :s]
    dict = Dict(:b => b, :e => e, :l => l, :s => s)
    for sym in bels
        range(sym, dict[sym])
        s1 = setdiff(bels, (sym,))
        for sym2 in s1
            range(Symbol(sym,sym2), dict[sym], dict[sym2])
            s2 = setdiff(s1, (sym2,))
            for sym3 in s2
                @test range(Symbol(sym,sym2,sym3), dict[sym], dict[sym2], dict[sym3]) == range(1, 5; length=5)
                @test range(Symbol(sym,sym2,sym3), dict[sym], dict[sym2], dict[sym3]) == range(:bels; sym => dict[sym], sym2 => dict[sym2], sym3 => dict[sym3])
                if sym === :l
                    @test range(Symbol(sym,sym2,sym3), dict[sym], dict[sym2]) == range(Symbol(sym,sym2), dict[sym], dict[sym2])
                    @test range(Symbol(sym,sym2,sym3), dict[sym]) == range(sym, dict[sym])
                end
            end
        end
    end
end

@testset "Length changing" begin
    @test length(1:5, 11) == range(1,5; length=11)
    @test length(1 => 5.5, 11) == range(1,5.5; length=11)
end

@testset "Step changing" begin
    @test step(1:5, 2) == range(1,5; step=2)
    @test step(1 => 5.5, 2) == range(1,5.5; step=2)
end

@testset "First and Last changing" begin
    @test first(1:5, 3) == 3:5
    @test last(1:5, 3) == 1:3
end


text:
VERSION
cd("C:/Users/Markus/.julia/dev/Hausch/")
@__DIR__()
pwd()
using Pkg
Pkg.activate(".")
# now add dependencies in the REPL via "add Module"
dirname(@__FILE__)
using Hausch
using Documenter
using DocumenterLaTeX

cd("C:/Users/Markus/.julia/dev/Hausch/docs/")

makedocs(modules = [Hausch],
        format = Documenter.HTML(assets = ["assets/style.css"], prettyurls = true),
        doctest = false,
        clean = false, # ?
        sitename = "Hausch.jl",
        authors = "Markus Hauschel",
        pages = [
            "Home" => "index.md",
            "Examples" => "examples.md",
            "Getting Started" => "getting_started.md",
            "Types" => "types.md",
            "Functions" => "functions.md",
            "References" => "references.md"
        ])

deploydocs(repo = "C:/Users/Markus/.julia/dev/Hausch/git",
                julia = "1.2.0", osname = "windows")

makedocs(modules = [Hausch],
        format = LaTeX(),
        doctest = false,
        clean = true, # ?
        sitename = "Hausch.jl",
        authors = "Markus Hauschel",
        pages = [
            "Home" => "index.md",
            "Examples" => "examples.md",
            "Getting Started" => "getting_started.md",
            "Types" => "types.md",
            "Functions" => "functions.md",
            "References" => "references.md"
        ])


text:
module GenML

abstract MLAlgorithm{F} #general type for any machine learning algorithm.
abstract Layer{F} <: MLAlgorithm{F} #a generalized "layer" type

# interface definition for MLAlgorithm.

# basics.  input and output vector size count.
inputs{T <: MLAlgorithm}(::Type{T})     = throw(ArgumentError("inputs for $T not properly implemented"))
outputs{T <: MLAlgorithm}(::Type{T})    = throw(ArgumentError("outputs for $T not properly implemented"))
parameters{T <: MLAlgorithm}(::Type{T}) = throw(ArgumentError("parameters for $T not properly implemented"))

inputs{T <: MLAlgorithm}(::T)     = inputs(T)
outputs{T <: MLAlgorithm}(::T)    = outputs(T)
parameters{T <: MLAlgorithm}(::T) = parameters(T)

@generated function flatten!{F}(::Vector{F}, mla::MLAlgorithm{F}, offset::Integer = 0)
  T = mla
  :(throw(ArgumentError("flatten! for $T not propertly implemented")))
end
@generated function unflatten!{F}(mla::MLAlgorithm{F}, ::Vector{F}, offset::Integer = 0)
  T = mla
  :(throw(ArgumentError("unflatten! for $T not propertly implemented")))
end

################################################################################
## nonparametric data storage.

doc"""
  `GenML.Storage(::Type{MLAlgorithm}, N = :v)`
  `GenML.Storage(::MLAlgorithm, N = :v)`

  takes an ML algorithm and allocates all non-parametric memory required to
  evaluate the algorithm.
"""
abstract Storage{T <: MLAlgorithm, N}

doc"""
  `GenML.BackpropStorage(::Type{MLAlgorithm})`
  `GenML.BackpropStorage(::MLAlgorithm)`

  takes an ML algorithm and allocates all non-parametric memory required to
  evaluate AND backpropagate the algorithm.  Implementations of backprop storage
  should be written in such a way that it be passed as a parameter to the
  forward evaluation which will persistently populate the layer data for the
  backwards propagation.
"""
abstract BackpropStorage{T, N} <: Storage{T, N}

#storage and backpropstorage for layer types is generally not allowed.
(::Type{Storage}){T <: Layer}(::Type{T})         = throw(ArgumentError("storage of layer types not allowed"))
(::Type{Storage}){T <: Layer}(::T)               = throw(ArgumentError("storage of layer types not allowed"))
(::Type{BackpropStorage}){T <: Layer}(::Type{T}) = throw(ArgumentError("storage of layer types not allowed"))
(::Type{BackpropStorage}){T <: Layer}(::T)       = throw(ArgumentError("storage of layer types not allowed"))

################################################################################
## evaluation

# MLAlgorithms should implement evaluate() methods.
@generated function evaluate!{F}(::AbstractArray{F}, mla::MLAlgorithm{F}, ::AbstractArray, ::Storage)
  if mla <: Layer
    :(throw(ArgumentError("Layer types don't implement storage evaluation.")))
  else
    :(throw(MethodError(evaluate!, (AbstractArray{F}, MLAlgorithm{F}, Array, Storage))))
  end
end

@generated function evaluate!{F}(output::AbstractArray{F}, mla::MLAlgorithm{F}, input::AbstractArray)
  #for the situation where the memory has not been preallocated, do the following.
  if input <: AbstractVector
    #catch mismatched input/output type parameters at the compile stage.
    (output <: AbstractVector) || throw(ArgumentError("vector/matrix status of input and output must match."))
    quote
      #generate the temporary storage
      sbuf = Storage(mla, :v)
      #go!
      evaluate!(output, mla, input, sbuf)
    end
  elseif input <: AbstractMatrix
    #catch mismatched input/output type parameters at the compile stage.
    (output <: AbstractMatrix) || throw(ArgumentError("vector/matrix status of input and output must match."))
    quote
      #do a second check to make sure the size of the batches are the same.
      size(output, 2) == size(input, 2) || throw(ArgumentError("input and output matrices must have matching batch sizes"))
      #generate the temporary storage with the correct batch size.
      sbuf = Storage(mla, size(input, 2))
      #go!
      evaluate!(output, mla, input, sbuf)
    end
  end
end

@generated function ml_call{F}(mla::MLAlgorithm{F}, input::Array)
  if input <: Vector
    quote
      output = Array{F,1}(outputs(mla))
      evaluate!(output, mla, input)
      output
    end
  elseif input <: Matrix
    quote
      output = Array{F,2}(outputs(mla), size(input, 2))
      evaluate!(output, mla, input)
      output
    end
  end
end

#optionally, an MLAlgorithm can implement backpropagation.
hasbackpropagation{T <: MLAlgorithm}(::Type{T}) = false  #default does not implement backpropagation

#to give backpropagation algorithms the option of not returning the next deltas for the last layer.
typealias VoidableDeltas{F} Union{AbstractArray{F}, Void}

#compound ml algorithms with backpropagation should implement this function.
@generated function backpropagate!{F}(mla::MLAlgorithm{F},                   #algorithm to be backpropagated.
                                      inp::AbstractArray,                    #input values to be trained
                                      oup::AbstractArray{F},                 #expected output values for training
                                      ouδ::AbstractArray{F},                 #deltas for these output values
                                      sto::BackpropStorage,                  #backprop storage data.
                                                                             #input deltas for passing to the next layer.
                                      inδ::VoidableDeltas{F} = nothing)
  T = mla
  :(throw(ArgumentError("backpropagate! for $T not properly implemented.")))
end

#layers should implement backpropagation without storage function.
function backpropagate!{F}(mla::Layer{F},
                           inp::AbstractArray,
                           oup::AbstractArray{F},
                           ouδ::AbstractArray{F},
                           inδ::VoidableDeltas{F} = nothing)
  T = mla
  :(throw(ArgumentError("backpropagate! for $T not properly implemented.")))
end

macro import_interface()
  quote
    #key functions in the interface
    import ..inputs; import ..outputs; import ..parameters; import ..flatten!; import ..unflatten!
    import ..evaluate!; import ..hasbackpropagation; import ..backpropagate!
    #key type definitions in the interface
    import ..MLAlgorithm; import ..Layer; import ..VoidableDeltas
    import ..Storage; import ..BackpropStorage
    #calling convenience function
    import ..ml_call
    #possibly useful modules.
    import ..TF; import ..CF
    #differentiation interface
    import ..d; import ..dxasy
  end
end

#mathematical utilities
#tools for differentiation
include("./differentiation.jl")
#transfer functions.
include("./transferfunctions.jl")
#include general cost functions
include("./costfunctions.jl")

#subtypes
include("./FCL/FCL.jl") #fully connected layer
include("./MLP/MLP.jl") #multilayer perceptrons

#include general optimizers
include("optimizers/optimizers.jl")

end # module

max_stars_repo_name: interplanetary-robot/GenML

max_issues_repo_name: interplanetary-robot/GenML

max_forks_repo_name: interplanetary-robot/GenML


text:
"""
The MovingWall boundary condition applies a bounce back rule and adds a velocity
to each distribution that is bounced back
"""
struct MovingWall{D <: Direction, Ints, T <: Real, VT <: AbstractVector{<:Real}} <:
       BoundaryCondition
    direction::D
    xs::Ints
    ys::Ints
    u::VT
    ρ::T
    T::T
end
MovingWall(direction, xs, ys, u::VT) where {T <: Real, VT <: AbstractVector{T}} =
    MovingWall(direction, xs, ys, u, one(T), one(T))

function apply!(bc::MovingWall{<:North}, q::Quadrature, f_new, f_old)
    cs = q.speed_of_sound_squared

    a_1_eq = equilibrium_coefficient(Val{1}, q, bc.ρ, bc.u, bc.T)
    nx, ny, nf = size(f_new)
    for f_idx in Base.OneTo(size(f_new, 3))
        H_1 = hermite(Val{1}, q.abscissae[:, f_idx], q)
        a_1 = q.weights[f_idx] * cs * dot(a_1_eq, H_1)

        opposite_f_idx = opposite(q, f_idx)

        for y_idx in Base.OneTo(ny)
            if y_idx + q.abscissae[2, opposite_f_idx] <= ny
                continue
            end
            for x_idx in Base.OneTo(nx)
                f_new[x_idx, y_idx, f_idx] = f_old[x_idx, y_idx, opposite_f_idx] + 2 * a_1
            end
        end

    end
end