text:
using Colors
export plaintext,
       fontfamily,
       fontsize,
       fontweight,
       fontcolor,
       fontstyle,
       fonttype,
       fontcase,
       serif,
       sansserif,
       slabserif,
       monospace,
       normal,
       slanted,
       italic,
       ucase,
       lcase,
       xxsmall,
       xsmall,
       small,
       medium,
       large,
       xlarge,
       xxlarge,
       bold,
       bolder,
       lighter,
       textalign,
       raggedright,
       raggedleft,
       justify,
       centertext,
       lineheight,
       letterspacing,
       heading,
       h1, h2, h3, h4,
       paragraph,
       emph,
       codeblock,
       code,
       blockquote,
       caption,
       title

plaintext(x) = Leaf(Elem(:div, string(x)))

abstract FontSize
# Note: I do not like these terms much. Could be irrational fear.
@terms FontSize begin
    xxsmall => XXSmall
    xsmall => XSmall
    small => Small
    medium => Medium
    large => Large
    xlarge => XLarge
    xxlarge => XXLarge
end
classes(::XXSmall) = "font-xx-small"
classes(::XSmall) = "font-x-small"
classes(::Small) = "font-small"
classes(::Medium) = "font-medium"
classes(::Large) = "font-large"
classes(::XLarge) = "font-x-large"
classes(::XXLarge) = "font-xx-large"

@api fontsize => (WithFontSize{T <: (@compat Union{FontSize, Length})} <: Tile) begin
    doc("Set the font size of text in one or more tiles")
    arg(size::T, doc="The font size")
    curry(tiles::TileList, doc="A tile or a list of tiles")
end
WithFontSize{T}(size::T, tiles) = WithFontSize{T}(size, tiles)

render{T <: Length}(t::WithFontSize{T}, state) =
    wrapmany(t.tiles, :span, state) & style(@d(:fontSize => t.size))
render{T <: FontSize}(t::WithFontSize{T}, state) =
    addclasses(wrapmany(t.tiles, :span, state), classes(t.size))

abstract FontWeight
# These terms add to the explosion as well.
@terms FontWeight begin
    bold => Bold
    bolder => Bolder
    lighter => Lighter
end
classes(::Bold) = "font-bold"
classes(::Bolder) = "font-bolder"
classes(::Lighter) = "font-lighter"

const allowed_font_weights = 100:100:900
@api fontweight => (WithFontWeight{T <: (@compat Union{Int, FontWeight})} <: Tile) begin
    doc("Set the font weight of text in one or more tiles.")
    arg(
        weight::T,
        doc="Font weight. Valid font weights are multiplies of 100 between 100 and 900."
    )
    curry(
        tiles::TileList,
        doc="A tile or a vector of tiles."
    )
end
WithFontWeight{T}(weight::T, tile) = WithFontWeight{T}(weight, tile)

render(t::WithFontWeight, state) = begin
    if !(t.weight in allowed_font_weights)
        error(string(t.weight, " is not an allowed font weight"))
    end
    wrapmany(t.tiles, :span, state) & style(@d(:fontWeight => t.weight))
end

render{T <: FontWeight}(t::WithFontWeight{T}, state) =
    addclasses(wrapmany(t.tiles, :span, state), classes(t.weight))

@api fontcolor => (FontColor <: Tile) begin
    doc("Set the font color.")
    arg(color::Color, doc="The color.")
    curry(tiles::TileList, doc="A tile or a vector of tiles.")
end

fontcolor(c::AbstractString) = fontcolor(parse(Colorant, c))
fontcolor(c::AbstractString, tiles) = fontcolor(parse(Colorant, c), tiles)

render(t::FontColor, state) =
    wrapmany(t.tiles, :span, state) & style(@d(:color => render_color(t.color)))

@api fontfamily => (FontFamily <: Tile) begin
    doc("Set the font family.")
    arg(family::AbstractString, doc="The font family")
    curry(tile::TileList, doc="A tile or a vector of tiles")
end
render(t::FontFamily, state) =
    wrapmany(t.tile, :span, state) & style(@d(:fontFamily => t.family))

abstract FontType
# TODO: Add serif and slab fonts
@terms FontType begin
    serif => Serif
    sansserif => SansSerif
    slabserif => SlabSerif
    monospace => Monospace
end
@api fonttype => (WithFontType <: Tile) begin
    doc("Set the font type.")
    arg(
        typ::FontType,
        doc=md"""The font type. Valid font types are `serif`, `sanserif`,
                 `slabserif` and `monospace`."""
    )
    curry(tiles::TileList, doc="A tile or a vector of tiles.")
end
classes(::Serif) = "font-serif"
classes(::SansSerif) = "font-sansserif"
classes(::SlabSerif) = "font-serif font-slab"
classes(::Monospace) = "font-monospace"

render(t::WithFontType, state) =
    addclasses(wrapmany(t.tiles, :span, state), classes(t.typ))

abstract FontStyle

@terms FontStyle begin
    normal => Normal
    slanted => Slanted
    italic => Italic
end
classes(::Normal) = "font-normal"
classes(::Slanted) = "font-slanted"
classes(::Italic) = "font-italic"

@api fontstyle => (WithFontStyle <: Tile) begin
    doc("Set the font style.")
    arg(
        style::FontStyle,
        doc=md"The font style. Valid font styles are `normal`, `slanted` and `italic`."
    )
    curry(tiles::TileList, doc="A tile or a vector of tiles.")
end
render(t::WithFontStyle, state) =
    addclasses(wrapmany(t.tiles, :span, state), classes(t.style))

abstract FontCase

# This exists purely for themability. (as opposed to using string case functions from Julia)
# e.g. .font-uppercase could sometimes use a different letter-spacing in a theme.
@terms FontCase begin
    ucase => Uppercase
    lcase => Lowercase
end
classes(::Uppercase) = "font-uppercase"
classes(::Lowercase) = "font-lowercase"

@api fontcase => (WithFontCase <: Tile) begin
    doc("Set the font case.")
    arg(case::FontCase, doc=md"The case. Valid cases are `ucase` and `lcase`")
    curry(tiles::TileList, doc="A tile or a vector of tiles.")
end

render(t::WithFontCase, state) =
    addclasses(wrapmany(t.tiles, :span, state), classes(t.case))

abstract TextAlignment

@terms TextAlignment begin
    raggedright => RaggedRight
    raggedleft => RaggedLeft
    justifytext => JustifyText
    centertext => CenterText
end

@api textalign => (AlignText{T <: TextAlignment} <: Tile) begin
    doc("Set the text alignment.")
    typedarg(
        alignment::T,
        doc=md"""The alignment. Valid alignments are `raggedright`, `raggedleft`,
                 `justifytext` and `centertext`."""
    )
    curry(tile::Tile, doc="A tile.")
end

render(t::AlignText{RaggedRight}, state) =
    render(t.tile, state) & style(@d(:textAlign => :left))
render(t::AlignText{RaggedLeft}, state) =
    render(t.tile, state) & style(@d(:textAlign => :right))
render(t::AlignText{JustifyText}, state) =
    render(t.tile, state) & style(@d(:textAlign => :justify))
render(t::AlignText{CenterText}, state) =
    render(t.tile, state) & style(@d(:textAlign => :center))

@api lineheight => (LineHeight <: Tile) begin
    doc("Set the height of lines of text.")
    arg(height::Length, doc="The height.")
    curry(tiles::TileList, doc="A tile or a vector of tiles.")
end
render(t::LineHeight, state) =
    wrapmany(t.tiles, :span, state) & style(@d(:lineHeight => t.height))

@api letterspacing => (LetterSpacing <: Tile) begin
    doc("Set the spacing between letters in text.")
    arg(space::Length, doc="The spacing.")
    curry(tiles::TileList, doc="A tile or a vector of tiles.")
end
render(t::LetterSpacing, state) =
    wrapmany(t.tiles, :span, state) & style(@d(:letterSpacing => t.space))

# Themable fonts

heading(n::Int, txt) = class("heading-$n", txt, forcewrap=true, wrap=:h1)
heading(n::Int) = t -> heading(n,t)

title(n::Int, txt) = class("title-$n", txt)

@apidoc title => (Class <: Tile) begin
    doc(md"Create a title.")
    arg(
        level::Int,
        doc="Title level. More is bigger. Valid values are integers 1 to 4."
    )
    curry(tile::TileList, doc="A tile or a vector of tiles.")
end

h1(txt) = heading(1, txt)
h2(txt) = heading(2, txt)
h3(txt) = heading(3, txt)
h4(txt) = heading(4, txt)

@apidoc heading => (Class <: Tile) begin
    doc(md"Create a heading. You can use `h1`, `h2`, `h3`, `h4` for brevity.")
    arg(level::Int, doc="Heading level.")
    curry(tile::TileList, doc="A tile or a vector of tiles.")
end

blockquote(txt) = class("blockquote", txt, forcewrap=true, wrap=:blockquote)
@apidoc blockquote => (Class <: Tile) begin
    doc("Create a quote block.")
    arg(tile::Tile, doc="The quote.")
end

caption(txt) = class("caption", txt, wrap=:span)
@apidoc caption => (Class <: Tile) begin
    doc("Create a caption.")
    arg(tile::Tile, doc="The caption.")
end

emph(txt) = class("emph", txt, forcewrap=true, wrap=:em)
@apidoc emph => (Class <: Tile) begin
    doc("Create a emphasized text.")
    arg(tile::Tile, doc="The text to be emphasized.")
end

@api code => (Code <: Tile) begin
    doc("Stylize text as code.")
    arg(code::Any, doc="The code.")
    kwarg(language::AbstractString="julia", doc="The language for syntax highlighting.")
end

## AbstractString will go inside a span, is this OK?
render(x::Code, state) = Elem(:code, render(x.code, state))

# colsize


text:
# Autogenerated wrapper script for psurface_jll for x86_64-w64-mingw32-cxx03
export libpsurface, psurface_convert, psurface_simplify, psurface_smooth

JLLWrappers.@generate_wrapper_header("psurface")
JLLWrappers.@declare_library_product(libpsurface, "libpsurface-0.dll")
JLLWrappers.@declare_executable_product(psurface_convert)
JLLWrappers.@declare_executable_product(psurface_simplify)
JLLWrappers.@declare_executable_product(psurface_smooth)
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libpsurface,
        "bin\\libpsurface-0.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_executable_product(
        psurface_convert,
        "bin\\psurface-convert.exe",
    )

    JLLWrappers.@init_executable_product(
        psurface_simplify,
        "bin\\psurface-simplify.exe",
    )

    JLLWrappers.@init_executable_product(
        psurface_smooth,
        "bin\\psurface-smooth.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# test simple solving
for i in 1:10
    A = rand(20, 20)
    sA = sparse(A)
    b = rand(20)

    @test fnnls(sA, b) ≈ fnnls(A, b)
    @test nnls(sA, b) ≈ nnls(A, b)
    @test pivot(sA, b) ≈ pivot(A, b)
    @test nonneg_lsq(sA,b;alg=:pivot, variant=:cache) ≈ nonneg_lsq(A,b;alg=:pivot, variant=:cache)
end


text:
"""
    iLQRSolver

A fast solver for unconstrained trajectory optimization that uses a Riccati recursion
to solve for a local feedback controller around the current trajectory, and then
simulates the system forward using the derived feedback control law.

# Constructor
    Altro.iLQRSolver(prob, opts; kwarg_opts...)
"""
struct iLQRSolver{T,I<:QuadratureRule,L,O,n,n̄,m,L1,C} <: UnconstrainedSolver{T}
    # Model + Objective
    model::L
    obj::O

    # Problem info
    x0::MVector{n,T}
    xf::MVector{n,T}
    tf::T
    N::Int

    opts::SolverOptions{T}
    stats::SolverStats{T}

    # Primal Duals
    Z::Traj{n,m,T,KnotPoint{T,n,m,L1}}
    Z̄::Traj{n,m,T,KnotPoint{T,n,m,L1}}

    # Data variables
    # K::Vector{SMatrix{m,n̄,T,L2}}  # State feedback gains (m,n,N-1)
    K::Vector{SizedMatrix{m,n̄,T,2,Matrix{T}}}  # State feedback gains (m,n,N-1)
    d::Vector{SizedVector{m,T,Vector{T}}}  # Feedforward gains (m,N-1)

    D::Vector{DynamicsExpansion{T,n,n̄,m}}  # discrete dynamics jacobian (block) (n,n+m+1,N)
    G::Vector{SizedMatrix{n,n̄,T,2,Matrix{T}}}        # state difference jacobian (n̄, n)

	quad_obj::TO.CostExpansion{n,m,T}  # quadratic expansion of obj
	S::TO.CostExpansion{n̄,m,T}         # Cost-to-go expansion
    E::TO.CostExpansion{n̄,m,T}         # cost expansion
    Q::TO.CostExpansion{n̄,m,T}         # Action-value expansion
    Qprev::TO.CostExpansion{n̄,m,T}     # Action-value expansion from previous iteration

    # Q_tmp::TO.QuadraticCost{n̄,m,T,SizedMatrix{n̄,n̄,T,2,Matrix{T}},SizedMatrix{m,m,T,2,Matrix{T}}}
    Q_tmp::TO.Expansion{n̄,m,T}
	Quu_reg::SizedMatrix{m,m,T,2,Matrix{T}}
	Qux_reg::SizedMatrix{m,n̄,T,2,Matrix{T}}
    ρ::Vector{T}   # Regularization
    dρ::Vector{T}  # Regularization rate of change

    cache::FiniteDiff.JacobianCache{Vector{T}, Vector{T}, Vector{T}, UnitRange{Int}, Nothing, Val{:forward}(), T}
    exp_cache::C
    grad::Vector{T}  # Gradient

    logger::SolverLogger
end

function iLQRSolver(
        prob::Problem{QUAD,T},
        opts::SolverOptions=SolverOptions(),
        stats::SolverStats=SolverStats(parent=solvername(iLQRSolver));
        kwarg_opts...
    ) where {QUAD,T}
    set_options!(opts; kwarg_opts...)

    # Init solver results
    n,m,N = size(prob)
    n̄ = RobotDynamics.state_diff_size(prob.model)

    x0 = prob.x0
    xf = prob.xf

    Z = prob.Z
    # Z̄ = Traj(n,m,Z[1].dt,N)
    Z̄ = copy(prob.Z)

	K = [zeros(T,m,n̄) for k = 1:N-1]
    d = [zeros(T,m)   for k = 1:N-1]

	D = [DynamicsExpansion{T}(n,n̄,m) for k = 1:N-1]
	G = [SizedMatrix{n,n̄}(zeros(n,n̄)) for k = 1:N+1]  # add one to the end to use as an intermediate result

    E = TO.CostExpansion{T}(n̄,m,N)
    quad_exp = TO.CostExpansion(E, prob.model)
    Q = TO.CostExpansion{T}(n̄,m,N)
    Qprev = TO.CostExpansion{T}(n̄,m,N)
    S = TO.CostExpansion{T}(n̄,m,N)

    # Q_tmp = TO.QuadraticCost{T}(n̄,m)
    Q_tmp = TO.Expansion{T}(n̄,m)
	Quu_reg = SizedMatrix{m,m}(zeros(m,m))
	Qux_reg = SizedMatrix{m,n̄}(zeros(m,n̄))
    ρ = zeros(T,1)
    dρ = zeros(T,1)

    cache = FiniteDiff.JacobianCache(prob.model)
    exp_cache = TO.ExpansionCache(prob.obj)
    grad = zeros(T,N-1)

    logger = SolverLogging.default_logger(opts.verbose >= 2)
	L = typeof(prob.model)
	O = typeof(prob.obj)
    solver = iLQRSolver{T,QUAD,L,O,n,n̄,m,n+m,typeof(exp_cache)}(
        prob.model, prob.obj, x0, xf,
		prob.tf, N, opts, stats,
        Z, Z̄, K, d, D, G, quad_exp, S, E, Q, Qprev, Q_tmp, Quu_reg, Qux_reg, ρ, dρ,
        cache, exp_cache, grad, logger)

    reset!(solver)
    return solver
end

# Getters
Base.size(solver::iLQRSolver{<:Any,<:Any,<:Any,<:Any,n,<:Any,m}) where {n,m} = n,m,solver.N
@inline TO.get_trajectory(solver::iLQRSolver) = solver.Z
@inline TO.get_objective(solver::iLQRSolver) = solver.obj
@inline TO.get_model(solver::iLQRSolver) = solver.model
@inline get_initial_state(solver::iLQRSolver) = solver.x0
@inline TO.integration(solver::iLQRSolver{<:Any,Q}) where Q = Q
solvername(::Type{<:iLQRSolver}) = :iLQR

log_level(::iLQRSolver) = InnerLoop

function reset!(solver::iLQRSolver{T}) where T
    reset_solver!(solver)
    solver.ρ[1] = 0.0
    solver.dρ[1] = 0.0
    return nothing
end


text:
function GaussIntegratorFactory()
    return GaussIntegratorFactory(())
end

function hermite(obj::GaussIntegratorFactory, arg0::jint)
    return jcall(obj, "hermite", SymmetricGaussIntegrator, (jint,), arg0)
end

function laguerre(obj::GaussIntegratorFactory, arg0::jint)
    return jcall(obj, "laguerre", GaussIntegrator, (jint,), arg0)
end

function legendre(obj::GaussIntegratorFactory, arg0::jint)
    return jcall(obj, "legendre", GaussIntegrator, (jint,), arg0)
end

function legendre(obj::GaussIntegratorFactory, arg0::jint, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "legendre", GaussIntegrator, (jint, jdouble, jdouble), arg0, arg1, arg2)
end

function legendre_high_precision(obj::GaussIntegratorFactory, arg0::jint)
    return jcall(obj, "legendreHighPrecision", GaussIntegrator, (jint,), arg0)
end

function legendre_high_precision(obj::GaussIntegratorFactory, arg0::jint, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "legendreHighPrecision", GaussIntegrator, (jint, jdouble, jdouble), arg0, arg1, arg2)
end


text:
"""

  function procrustes(x,y) -> xnew

  aligns two structures [sets of points in 3D space]. Solves
  the "Procrustes" problem. Structures are expected to be of the same size, and the
  correspondence is assumed from the vector indices.

  Returns x aligned, by performing the rigid body transformation [rotation
  and translation that minimizes the RMSD between x and y].

  x, y, and xnew (return) are matrices of dimensions (n,3)
  (n is the number of points, 3 is the dimension of the space).

  L. Martinez, Institute of Chemistry - University of Campinas
  Jan 04, 2019

"""

using LinearAlgebra

function procrustes( x :: Matrix{}, y :: Matrix{} )

  n = size(x,1)

  # Computing centroid

  cmx = zeros(3)
  cmy = zeros(3)
  for i in 1:n
    for j in 1:3
      cmx[j] = cmx[j] + x[i,j]
      cmy[j] = cmy[j] + y[i,j]
    end
  end
  cmx = cmx / n
  cmy = cmy / n

  # Translating both sets to the origin

  for i in 1:n
    for j in 1:3
      x[i,j] = x[i,j] - cmx[j]
      y[i,j] = y[i,j] - cmy[j]
    end
  end

  # Computing the quaternion matrix

  xm = Vector{Float64}(undef,n)
  ym = Vector{Float64}(undef,n)
  zm = Vector{Float64}(undef,n)
  xp = Vector{Float64}(undef,n)
  yp = Vector{Float64}(undef,n)
  zp = Vector{Float64}(undef,n)
  for i in 1:n
    xm[i] = y[i,1] - x[i,1]
    ym[i] = y[i,2] - x[i,2]
    zm[i] = y[i,3] - x[i,3]
    xp[i] = y[i,1] + x[i,1]
    yp[i] = y[i,2] + x[i,2]
    zp[i] = y[i,3] + x[i,3]
  end

  q = zeros(4,4)
  for i in 1:n
    q[1,1] = q[1,1] + xm[i]^2 + ym[i]^2 + zm[i]^2
    q[1,2] = q[1,2] + yp[i]*zm[i] - ym[i]*zp[i]
    q[1,3] = q[1,3] + xm[i]*zp[i] - xp[i]*zm[i]
    q[1,4] = q[1,4] + xp[i]*ym[i] - xm[i]*yp[i]
    q[2,2] = q[2,2] + yp[i]^2 + zp[i]^2 + xm[i]^2
    q[2,3] = q[2,3] + xm[i]*ym[i] - xp[i]*yp[i]
    q[2,4] = q[2,4] + xm[i]*zm[i] - xp[i]*zp[i]
    q[3,3] = q[3,3] + xp[i]^2 + zp[i]^2 + ym[i]^2
    q[3,4] = q[3,4] + ym[i]*zm[i] - yp[i]*zp[i]
    q[4,4] = q[4,4] + xp[i]^2 + yp[i]^2 + zm[i]^2
  end
  q[2,1] = q[1,2]
  q[3,1] = q[1,3]
  q[3,2] = q[2,3]
  q[4,1] = q[1,4]
  q[4,2] = q[2,4]
  q[4,3] = q[3,4]

  # Computing the eigenvectors 'v' of the q matrix

  v = LinearAlgebra.eigvecs(q)

  # Compute rotation matrix

  u = Matrix{Float64}(undef,3,3)
  u[1,1] = v[1,1]^2 + v[2,1]^2 - v[3,1]^2 - v[4,1]^2
  u[1,2] = 2. * ( v[2,1]*v[3,1] + v[1,1]*v[4,1] )
  u[1,3] = 2. * ( v[2,1]*v[4,1] - v[1,1]*v[3,1] )
  u[2,1] = 2. * ( v[2,1]*v[3,1] - v[1,1]*v[4,1] )
  u[2,2] = v[1,1]^2 + v[3,1]^2 - v[2,1]^2 - v[4,1]^2
  u[2,3] = 2. * ( v[3,1]*v[4,1] + v[1,1]*v[2,1] )
  u[3,1] = 2. * ( v[2,1]*v[4,1] + v[1,1]*v[3,1] )
  u[3,2] = 2. * ( v[3,1]*v[4,1] - v[1,1]*v[2,1] )
  u[3,3] = v[1,1]^2 + v[4,1]^2 - v[2,1]^2 - v[3,1]^2

  # Rotate vector x [will be stored in xnew], and restore y

  xnew = zeros(n,3)
  for i in 1:n
    for j in 1:3
      for k in 1:3
        xnew[i,j] = xnew[i,j] + u[j,k] * x[i,k]
      end
    end
  end

  # Translate vector to the centroid of y [and restore x and y]

  for i in 1:n
    for j in 1:3
      xnew[i,j] = xnew[i,j] + cmy[j]
      y[i,j] = y[i,j] + cmy[j]
      x[i,j] = x[i,j] + cmx[j]
    end
  end

  return xnew

end


text:
#####
##### Sources
#####

eq_tends(pv::AbstractPrognosticVariable, m::LandModel, tt::Source) =
    (m.source_dt[pv]...,)

#####
##### First order fluxes
#####

eq_tends(::PV, ::LandModel, ::Flux{FirstOrder}) where {PV} = ()

#####
##### Second order fluxes
#####

# Empty by default
eq_tends(pv::PV, ::AbstractSoilComponentModel, ::Flux{SecondOrder}) where {PV} =
    ()

eq_tends(pv::PV, land::LandModel, tt::Flux{SecondOrder}) where {PV} =
    (eq_tends(pv, land.soil.heat, tt)..., eq_tends(pv, land.soil.water, tt)...)

eq_tends(
    pv::PV,
    land::LandModel,
    tt::Flux{SecondOrder},
) where {PV <: VolumetricInternalEnergy} =
    (eq_tends(pv, land.soil.heat, land.soil.water, tt)...,)

# TODO: move to soi_heat.jl?
eq_tends(
    ::VolumetricInternalEnergy,
    ::SoilHeatModel,
    ::SoilWaterModel,
    ::Flux{SecondOrder},
) = (DiffHeatFlux(), DarcyDrivenHeatFlux())

# TODO: move to soi_water.jl?
eq_tends(
    ::VolumetricInternalEnergy,
    ::SoilHeatModel,
    ::PrescribedWaterModel,
    ::Flux{SecondOrder},
) = (DiffHeatFlux(),)

eq_tends(::VolumetricLiquidFraction, ::SoilWaterModel, ::Flux{SecondOrder}) =
    (DarcyFlux(),)


text:
# implementation of efficient direct methods for symmetric Toeplitz matrices
using LinearAlgebra
using LinearAlgebra: checksquare
using ToeplitzMatrices

# using ScalableGaussianProcesses
# using FastGaussianProcesses

# computes y = K \ (-r) where K = SymmetricToeplitz([1, r[1:end-1]])
# K is assumed to be positive definite
# r is of length n
durbin(r::AbstractVector) = durbin!(zero(r), r)
# z is a temporary that's necessary for the computation
function durbin!(y::AbstractVector, r::AbstractVector)
    n = length(r)
    length(y) == n || throw(DimensionMismatch("length(y) = $(length(y)) ≠ $n = length(r)"))
    y[1] = -r[1]
    α, β = -r[1], 1
    for k in 1:n-1
        β *= (1-α^2)
        r_k, y_k = @views r[1:k], y[1:k]
        α = -(r[k+1] + reverse_dot(r_k, y_k)) / β
        reverse_increment!(y_k, y_k, α)
        y[k+1] = α
    end
    return y
end

# Trench's algorithm (see page 213 of Golub & van Loan)
# computes inverse of symmetric positive definite Toeplitz matrix
function trench(T::SymmetricToeplitz)
    n = checksquare(T)
    B = zeros(eltype(T), n, n)
    trench!(B, T)
end

function trench!(B::AbstractMatrix, T::SymmetricToeplitz)
    r_0 = T.vc[1]
    r = @view T.vc[2:end]
    if r_0 == 1 # levinson implementation assumes normalization of diagonal
        r /= r_0
    end
    trench!(B, r)
    B *= r_0
end

# computes inverse of K = SymmetricToeplitz(vcat(1, r))
function trench(r::AbstractVector)
    n = length(r) + 1
    B = zeros(eltype(r), n, n)
    trench!(B, r)
end

# computes inverse of K = SymmetricToeplitz(vcat(1, r))
# uses B to store the inverse
# NOTE: only populates entries on upper triangle since the inverse is symmetric
function trench!(B::AbstractMatrix, r::AbstractVector)
    n = checksquare(B)
    n == length(r) + 1 || throw(DimensionMismatch())
    y = durbin(r)
    γ = inv(1 + dot(r, y))
    ν = γ * y[end:-1:1]
    B[1, 1] = γ
    @. B[1, 2:n] = γ * y
    for j in 2:n
        for i in 2:j
            @inbounds B[i, j] = B[i-1, j-1] + (ν[n+1-j] * ν[n+1-i] - ν[i-1] * ν[j-1]) / γ
        end
    end
    return Symmetric(B)
end

# solves K \ b where K_{ij} = r[abs(i-j)] and by assumption r[0] = 1
# i.e. K = SymmetricToeplitz(vcat(1, r))
# also, K is assumed to be positive definite
levinson(r::AbstractVector, b::AbstractVector) = levinson!(zero(b), r, b)
function levinson!(x::AbstractVector, r::AbstractVector, b::AbstractVector,
                   y::AbstractVector = zero(x))
    n = length(r) + 1
    length(x) == n || throw(DimensionMismatch("length(x) = $(length(x)) ≠ $n = length(r) + 1"))
    length(b) == n || throw(DimensionMismatch("length(b) = $(length(b)) ≠ $n = length(r) + 1"))
    y[1] = -r[1]
    x[1] = b[1]
    α, β = -r[1], 1
    @inbounds for k in 1:n-1
        β *= (1-α^2)
        r_k, x_k, y_k  = @views r[1:k], x[1:k], y[1:k]
        μ = (b[k+1] - reverse_dot(r_k, x_k)) / β
        reverse_increment!(x_k, y_k, μ)
        x[k+1] = μ
        if k < n-1
            α = -(r[k+1] + reverse_dot(r_k, y_k)) / β
            reverse_increment!(y_k, y_k, α) # computes y_k += α * reverse(y_k)
            y[k+1] = α
        end
    end
    return x
end

function levinson(T::SymmetricToeplitz, b::AbstractVector)
    r_0 = T.vc[1]
    r = @view T.vc[2:end]
    if r_0 == 1 # levinson implementation assumes normalization of diagonal
        r /= r_0
    end
    x = levinson(r, b)
    if r_0 != 1
        x ./= r_0
    end
    return x
end

# computes dot(x, reverse(y)) efficiently
function reverse_dot(x::AbstractArray, y::AbstractArray)
    n = length(x)
    n == length(y) || throw(DimensionMismatch())
    d = zero(promote_type(eltype(x), eltype(y)))
    @inbounds @simd for i in 1:n
        d += x[i] * y[n-i+1]
    end
    return d
end
# computes x += α * reverse(y) efficiently without temporary
# NOTE: works without allocations even when x === y
function reverse_increment!(x::AbstractArray, y::AbstractArray = x, α::Number = 1)
    n = length(x)
    n == length(y) || throw(DimensionMismatch())
    if x === y
        @inbounds @simd for i in 1:n÷2
            y_i = y[i] # important to store these for the iteration so that they aren't mutated
            z_i = y[n-i+1]
            x[i] += α * z_i
            x[n-i+1] += α * y_i
        end
        if isodd(n)
            i = (n÷2) + 1 # midpoint
            x[i] += α * y[i]
        end
    else
        @inbounds @simd for i in 1:n
            x[i] += α * y[n-i+1]
        end
    end
    return x
end

# via Bareiss algorithm
# function LinearAlgebra.det(T::SymmetricToeplitz)
#     cholesky()
# end


text:
module Ensembles

import ..ED
using ..Levels
using ..OPDiag
using ..EON
using ..SimLib
using XXZNumerics

using Printf: @sprintf
using Statistics: mean
using LinearAlgebra
import JLD2

export ENSEMBLE_INDICES, EnsembleDataDescriptor, EnsembleData, ensemble_predictions, load_ensemble

## Data structure

"""
    struct EnsembleDataDescriptor

Carries the information to construct a [`EDDataDescriptor`](!ref) object.
 - geometry
 - dimension
 - system_size
 - α
 - shots
 - densities ρs
 - field strengths
 - scaling of field strengths
 - symmetries to respect
For `load`ing data only the first 4 fields are required.
Can also be constructed from a `PositionDataDescriptor` by supplying a the missing bits (α, fields).
"""
struct EnsembleDataDescriptor <: ED.EDDerivedDataDescriptor
    derivedfrom::ED.EDDataDescriptor
    EnsembleDataDescriptor(args...; kwargs...) = new(EDDataDescriptor(args...; kwargs...))
end

const ENSEMBLE_INDICES = (; microcanonical = 1, canonical = 2, diagonal = 3)

"""
    struct EnsembleData

Stores the actual data for the ensemble's predictions specified by the descriptor [`EnsembleDataDescriptor`](@ref).

# Index order
The indices mean:
 - shot number
 - field index
 - density ρ
 - ensemble
 For access use the keys :microcanonical, :canonical and :diagonal either on this struct or on ENSEMBLE_INDICES to map to the index)

The default save directory is "ensemble".
"""
struct EnsembleData{N} <: SimLib.AbstractSimpleData
    descriptor::EnsembleDataDescriptor
    # [shot, h, rho, ensemble]
    # ensemble: 1=microcanonical, 2=canonical, 3=diag
    data::Array{Float64,N}
end

_data_indices(arr) = CartesianIndices(axes(arr)[1:end-1])

function Base.getproperty(ensdata::EnsembleData, s::Symbol)
    if hasfield(typeof(ensdata), s)
        getfield(ensdata, s)
    elseif haskey(ENSEMBLE_INDICES, s)
        @view ensdata.data[_data_indices(ensdata.data), ENSEMBLE_INDICES[s]]
    else
        getproperty(ensdata.descriptor, s)
    end
end

## Saving/Loading
ED._default_folder(::EnsembleDataDescriptor) = "ensemble"

load_ensemble(args...; kwargs...) = load(EnsembleDataDescriptor(args...; kwargs...))
## main

# eon [eigen_state, shot, h, rho]
# eev [eigen_state, shot, h, rho]
# out              [shot, h, rho]

function _canonical!(out, eon, eev, evals)
    # out[i,j,k] = canonical_weights*eev[a,i,j,k]
    @views for I in CartesianIndices(out)
        E_0 = dot(eon[:, I], evals[:, I])
        occ = canonical_ensemble(evals[:, I], E_0) # ΔE default: 0.5% of total spectral width
        out[I] = @views dot(occ, eev[:, I])
    end
    out
end

function _microcan!(out, eon, eev, evals)
    # out[i,j,k] = microcanonical_weights*eev[a,i,j,k]
    @views for I in CartesianIndices(out)
        E_0 = dot(eon[:, I], evals[:, I])
        if E_0 < evals[1,I] || E_0 > evals[end,I] || !issorted(evals[1,I])
            logmsg("Something off for index $I")
            @show issorted(evals[:,I])
            @show norm(sort(evals[:, I]) - evals[:, I])
            @show E_0, evals[1,I], evals[end,I], minimum(evals[:,I]), maximum(evals[:,I])
            @show eltype(eon)
            @show sum(eon[:, I])
            println()
        end
        occ = microcanonical_ensemble(evals[:, I], E_0)
        out[I] = dot(occ, eev[:, I])
    end
    out
end

function _diagonal!(out, eon, eev, ignored=nothing)
    # out[i,j,k] = eon[a,i,j,k]*eev[a,i,j,k]
    for I in CartesianIndices(out)
        out[I] = @views dot(eon[:, I], eev[:, I])
    end
    out
end

function ensemble_predictions(evals, eon, eev)
    ensemble_data = zeros(Float64, size(eon)[2:end]..., 3)
    ensemble_predictions!(ensemble_data, evals, eon, eev)
    return ensemble_data
end

function ensemble_predictions!(ensemble_data, evals, eon, eev)
    if !(eltype(eon) <: Real)
        eon = abs2.(eon)
        logmsg("Ensemble prediction: Got amplitudes -> squaring.")
    end
    I = _data_indices(evals)
    _microcan!( @view(ensemble_data[I, 1]), eon, eev, evals)
    _canonical!(@view(ensemble_data[I, 2]), eon, eev, evals)
    _diagonal!( @view(ensemble_data[I, 3]), eon, eev, evals) # evals not needed here but looks nicer ;)
end

function ensemble_predictions(levels::LevelData, eon::EONData, eev::OPDiagData)
    data = ensemble_predictions(levels.data, eon.data, eev.data)
    return EnsembleData(EnsembleDataDescriptor(descriptor(levels).derivedfrom), data)
end

function SimLib.create(desc::EnsembleDataDescriptor)
    eddesc = desc.derivedfrom
    logmsg("Computing ensemble predictions for $(eddesc)")
    levels = load(LevelDataDescriptor(eddesc))
    eon = load(EONDataDescriptor("xpol", eddesc))
    eev = load(OPDiagDataDescriptor("xmag", eddesc))
    return ensemble_predictions(levels, eon, eev)
end

end #module


text:
"""
    DiscreteSystemSetup(resampling = NoResampling, dt::Int = 1, Ttr::Int = 0,
        n_pts = 100, source, target)

Setup for causality analysis of discrete dynamical systems.

## Mandatory keyword arguments

- **`source`**:  The variable(s) to use as `source` when calling `causality`. Usually integer(s).

- **`target`**:  The variable(s) to use as `target` when calling `causality`. Usually integer(s).

## Optional keywords

- **`resampling`**: An instance of a resampling scheme. Defaults to `NoResampling()`.

- **`dt::Int`**: The time step when iterating the system.

- **`Ttr::Int`**: The number of transient iterations before starting sampling.

- **`n_pts::Int`**: The number of points in the orbit.
"""
Base.@kwdef struct DiscreteSystemSetup{R} <: DynamicalSystemSetup where R
    """ An instance of a resampling scheme (defaults to `NoResampling()`)"""
    resampling::R = NoResampling()

    """ The time step """
    dt::Int = 1

    """ The number of transient iterations before starting sampling """
    Ttr::Int = 0

    """ The number of points in the orbit """
    n_pts = 100

    """ The variable(s) to use as `source` """
    source

    """ The variable(s) to use as `target` """
    target
end

export DiscreteSystemSetup


text:
###DEPRECATED###VALID FUNCTIONS MOVE TO PREPROCESSING###WILL BE REMOVED SOON###

"""
    function adjust_model(database, variation_tres =.5, dGr_dict =Dict())

DEV NOTE: Necessary for MinFlux/MinRxn procedure, currently work in progress. Function might change severely.
Reduces model to usable reactions. It excludes exchange reactions, those where no Gibbs free energy could be calculated or with too large standard deviation, depending on the given coefficient of variation threshold (`variation_tres`).
"""
function adjust_model(database, variation_tres = 0.9)
    dGr_dict, dGr_bounds = collect_dGr_bounds(database)

    for rxn in reactions(database)
        if !reaction_mass_balanced(database, rxn) #excluding exchanges
            database = remove_reaction(database, rxn)
        end
    end
    mod_bounds = OrderedDict()
    for (k, v) in dGr_bounds
        if k in reactions(database)
            mod_bounds[k] = v
        end
    end
    for rxn in keys(mod_bounds)
        if typeof(dGr_dict[rxn]) == String #excluding reactions without calculated Gibbs energies
            mod_bounds[rxn] = (0.0, 0.0)
        else
            try
                if dGr_dict[rxn].val.err / dGr_dict[rxn].val.val > variation_tres || abs(dGr_dict[rxn].val.val) == Inf
                    mod_bounds[rxn] = (0.0, 0.0)
                end
            catch
                mod_bounds[rxn] = (0.0, 0.0)
            end
        end
    end

    return database, mod_bounds
end

max_stars_repo_name: VictorVeraFrazao/OptStoic

max_issues_repo_name: VictorVeraFrazao/OptStoic

max_forks_repo_name: VictorVeraFrazao/OptStoic


text:
#TODO deleting of solutions does not delete from h.

"""
An "implementation" of a `ContinuationMethod` using a predictor-corrector-method with
step-size adaption. Look up the theoretical documentation for more details.
"""
type PC <: ContinuationMethod
	parent::Session
	dataGUI::Dict{AbstractString, Any}
	h::Dict{Any,Float64}
	PC(ses::Session) = new(ses, Dict{AbstractString, Any}(), Dict{Any,Float64}())
end

function Base.show(C::PC)
	win = @Window("Continuation Controls", 256, 100, false)
	# setproperty!(win, :resizable, false)

	g = @Grid() #TODO obsolete
	push!(win, g)

	buttonStep = @Button("Step")
	buttonRun = @ToggleButton("Run")

	gridPC = mkControlGrid(C.dataGUI, Array[
		#[ ("ϵ", Float64, 1e-4, 1e-6:1e-8:1e-2), Void ],
		[ ("κ", Float64, .4, 1e-3:1e-3:1.0) ],
		[ ("δ", Float64, .5, 1e-2:1e-2:20.0) ],
		[ ("α", Float64, 10.0, .1:.1:40.0) ],
		[ ("Perturb", Bool, false); ("Strength", Float64, .0, -.1:1e-3:.1) ],
		[ buttonStep, buttonRun ]
	])
	g[1:2,2] = gridPC

	signal_connect(buttonStep, "clicked") do w
		@schedule begin
			C.parent.P.activeSolution == Void && return Void
			setproperty!(w, :sensitive, false)
			try
				step(C)
			finally
				setproperty!(w, :sensitive, true)
			end
		end
	end

	signal_connect(buttonRun, "toggled") do w
		if getproperty(w, :active, Bool)
			@schedule while getproperty(w, :active, Bool)
				try step(C)
				catch e
					setproperty!(w, :active, false)
					return
				end
				sleep(.2)
			end
		end
	end

	showall(win)
end

# can throw error, eg when h=0
function Base.step(C::PC)
	V = C.parent.P.activeSolution
	V == Void && return Void #TODO elaborate..

	local B, flagNewBranch = false
	if typeof(V) == Solution && (V == V.parent[1] || V == V.parent[end])
		B = V.parent
	else # internal, single or orphaned solution
		flagNewBranch = true
		B = Branch(C.parent.P)
		V = Solution(V, B)
		B.solutions = Solution[V]
	end

	htmp = get!(C.h, V, 1.0)
	V==B[1] && length(B)>1 && (htmp=-htmp)

	Htmp = C.dataGUI["Perturb"] ? V -> H(C.parent.core)(V) + C.dataGUI["Strength"] : H(C.parent.core)
	W,htmp = continuationStep(Htmp, J(C.parent.core), V.data, htmp, 1e-8, C.dataGUI["κ"], C.dataGUI["δ"], C.dataGUI["α"]) #fetch?

	W = Solution(W, B)
	C.h[W] = abs(htmp)

	flagNewBranch && push!(C.parent.P, B) #NOTE no branch pushed in case of error in continuationStep
	(htmp>0?push!:unshift!)(B, W)
	setActiveSolution(C.parent.P, W)

	return Void
end

# numerical core
function continuationStep(H, J, V, h, ϵ, κ, δ, α)
	# tangent-vector for matrix A
	function tang(A)
		t = [A; ones(1,size(A,2))] \ [zeros(size(A,1)); 1] # nullspace(A)[:,1] #takes up >90% processing time
		flipsign(t / norm(t), det([A; t']))
	end

	# newton corrector step
	function Δ(v)
		Hv, Jv = H(v), J(v)
		Hv, Jv, Jv \ Hv
	end

	local v₁, dir = nothing
	while true
		tJV = tang(J(V))
		v₀ = V + h * tJV

		# step size control
		Hv₀, Jv₀, Δv₀ = Δ(v₀)
		v₁ = v₀ - Δv₀

		Hv₁, Jv₁, Δv₁ = Δ(v₁)
		v₁ = v₁ - Δv₁
		κ′ = norm(Δv₁) / norm(Δv₀)
		δ′ = norm(Δv₀)
		α′ = acos( clamp(dot(tJV, tang(Jv₀)), -1, 1) ) #prevent domain exception through numerical error

		f = max(sqrt(κ′/κ), sqrt(δ′/δ), (α′/pi*180)/α)

		# prevent infinite loop: only allow changes of h in one direction
		dir == nothing && (dir = (f < 1.0))
		dir ≠ (f < 1.0) && break

		h /= clamp(f, .5, 2)
		.5 ≤ f ≤ 2.0 && break
		abs(h) ≤ 1e-8 && error("h≈0 !") #TODO remove?
	end

	its = 0
	v₁ = newton(H, J, v₁, predEps(ϵ) ∧ predCount(50); callback=(H,J,v)-> (its+=1)) #TODO integrate max count
	its ≥ 50 && warn("no convergence; error: ", norm(H(v₁)))
	return v₁, h
end


text:
module Heading

import Genie, Stipple

export heading

function heading(title::String = "";
                  class = "",
                  img = Genie.Renderer.Html.img(class="st-logo", src="/img/st-logo.svg"),
                  h = Genie.Renderer.Html.h1(title, class="st-header__title text-h3"),
                  kwargs...)

  Genie.Renderer.Html.header(; class="$(class) st-header q-pa-sm", kwargs...) do; [
    img
    h
  ]end
end

end


text:
module Circuits

import QXZoo

export create_test_circuit, create_qft_circuit, create_ghz_circuit
export create_rqc_circuit, gate_matrix, gate_qubits

"""
Functions to create circuits
"""

"""
    create_test_circuit()

Function to create a 3 qubit QXZoo ghz state preparation circuit useful for testing
"""
function create_test_circuit()
    circ = QXZoo.Circuit.Circ(3)
    QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.h(1))
    QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.c_x(1, 2))
    QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.c_x(2, 3))
    circ
end

"""
    create_tescreate_ghz_circuitt_circuit(n::Int64)

Function to create an n qubit QXZoo ghz state preparation circuit
"""
function create_ghz_circuit(n::Int64)
    circ = QXZoo.Circuit.Circ(n)
    QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.h(1))
    for i in 1:n-1
        QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.c_x(i, i+1))
    end
    circ
end

"""
    create_qft_circuit(n::Integer)

Function to create a qft circuit with n qubits
"""
function create_qft_circuit(n::Integer)
    circ = QXZoo.Circuit.Circ(n)
    QXZoo.QFT.apply_qft!(circ, collect(0:n-1))
    QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.c_x(0, 1))
    QXZoo.Circuit.add_gatecall!(circ, QXZoo.DefaultGates.c_x(1, 2))
    circ
end

"""
    create_rqc_circuit(rows::Int,
                       cols::Int,
                       depth::Int,
                       seed::Union{Int, Nothing}=nothing;
                       final_h::Bool=false)

Function to create a qft circuit with n qubits
"""
function create_rqc_circuit(rows::Int,
                            cols::Int,
                            depth::Int,
                            seed::Union{Int, Nothing}=nothing;
                            final_h::Bool=false)
    return QXZoo.RQC.create_RQC(rows, cols, depth, seed; final_Hadamard_layer=final_h)
end

"""
    gate_matrix(gate::QXZoo.GateOps.AGateCall)

Function to get the matrix for a gate
"""
function gate_matrix(gate::QXZoo.GateOps.AGateCall)
    convert(Array{ComplexF64}, QXZoo.GateMap.gates[gate.gate_symbol]())
end

"""
    gate_qubits(gate::QXZoo.GateOps.AGateCall)

Function to get the qubit indices that the gate acts on
"""
function gate_qubits(gate::QXZoo.GateOps.AGateCall)
    if typeof(gate) <: QXZoo.GateOps.GateCall1
        qubits = [gate.target]
    elseif typeof(gate) <: QXZoo.GateOps.GateCall2
        qubits = [gate.target, gate.ctrl]
    end
    qubits
end

end  # end module


text:
using Revise
using ADCME
using NNFEM
using JLD2
using PyPlot
using PyCall
using JLD2
using MAT
using Statistics

rcParams = PyPlot.PyDict(PyPlot.matplotlib."rcParams")
font0 = Dict(
        "font.size" => 16,
        "axes.labelsize" => 16,
        "xtick.labelsize" => 16,
        "ytick.labelsize" => 16,
        "legend.fontsize" => 16,
)
merge!(rcParams, font0)
np = pyimport("numpy")
mpl = pyimport("tikzplotlib")

reset_default_graph()
include("nnutil.jl")
t_scale, s_scale = 1000.0, 1000.0 #ms->s, MPa->GPa
markevery = 5
Nite = 10
max_idx = 4
max_nn_init_id = 10
cols = ["r", "y", "b", "g"]
nnname=nntype
if nntype == "piecewise"
    nnname = "CholNN"
elseif nntype == "ae_scaled"
    nnname = "NN1"
elseif nntype == "ae_scaled2"
    nnname = "NN2"
end

function plot_loss(nntype::String)
    close("all")
    Method = "Min"  #"Median" #"Min"
    min_nn_init_ids = zeros(Int64, max_idx)

    for idx = 1:max_idx
        loss_array = zeros(Float64, max_nn_init_id)
        # minLoss = 1e100
        for nn_init_id = 1:max_nn_init_id
            file="Data/$(nntype)/learned_nn$(idx)_loss_$(nn_init_id).txt"
            vars = matread(file)
            train_loss = vars["loss"]

            # if train_loss[end] < minLoss
            #     minLoss = train_loss[end]
            #     min_nn_init_ids[idx] = nn_init_id
            # end
            loss_array[nn_init_id] = train_loss[end]

            N = length(train_loss)
            NN = 1:div(N,100):N
            if nn_init_id == 1
                semilogy(NN, train_loss[NN],"--$(cols[idx])", label="$(idx)-layer $(nnname)")
            else
                semilogy(NN, train_loss[NN],"--$(cols[idx])")
            end
        end

        if Method == "Min"
            min_nn_init_ids[idx] = argmin(loss_array)
        elseif Method == "Median"
            min_nn_init_ids[idx] = argmin(abs.(loss_array .- median(loss_array)))
        else
            error("Method is ", Method)
        end

    end

    legend()
    xlabel("Iteration")
    ylabel("Loss")

    PyPlot.tight_layout()
    savefig("nnlearn_$(nntype)_loss.pdf")
    mpl.save("nnlearn_$(nntype)_loss.tex")
    @show min_nn_init_ids
    return min_nn_init_ids
end

function plot_nnlearn_p2p_strain_stress(tid::Int64, nntype::String, min_nn_init_ids::Array{Int64,1})
    close("all")
    sid = 8
    strain, stress = read_strain_stress("Data/$(tid).dat")
    X, Y = prepare_strain_stress_data1D(strain, stress )
    NT = Int64(length(Y)/sid)

    plot(X[NT*(sid-1)+1:markevery:NT*sid, 1], Y[NT*(sid-1)+1:markevery:NT*sid]/s_scale, "o",fillstyle="none", label="Reference")

    for idx = 1:max_idx

        reset_default_graph()
        sess = Session(); init(sess);

        # reset config in nnutil.jl

        if idx == 1
            config = [20,1]
        elseif idx == 2
            config = [20,20,1]
        elseif idx == 3
            config = [20,20,20,1]
        elseif idx == 4
            config = [20,20,20,20,1]
        elseif idx == 5
            config = [20,20,20,20,20,1]
        else
            error(idx, "idx <= 5")
        end

        function nn(ε, ε0, σ0)
            # make sure all inputs are 2d matrix
            @show ε, ε0, σ0

            local y
            if nntype=="linear"
                y = E0 * ε
            elseif nntype=="ae_scaled"
                x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
                y = ae(x, config, "ae_scaled")*stress_scale
            elseif nntype=="ae_scaled2"
                x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
                y = ae(x, config, "ae_scaled2")*stress_scale + σ0
            elseif nntype=="piecewise"
                x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
                H = ae(x, config, "piecewise")^2
                s = σ0^2
                i = sigmoid((s - 0.01e6)/1.0e3)

                y = ( H.* i + E0 * (1-i) ) .* (ε-ε0)*stress_scale/strain_scale + σ0

            elseif nntype=="piecewise2"
                x = [ε/strain_scale ε0/strain_scale σ0/stress_scale]
                H = ae(x, config, "piecewise2")
                s = σ0^2
                i = sigmoid((s - 0.01e6))

                y = ( H.* i + E0 * (1-i) ) .* (ε-ε0)*stress_scale/strain_scale + σ0
            else
                error("nntype must be specified.")
            end

            return y

        end

        y = squeeze(nn(constant(X[:,1:1]), constant(X[:,2:2]), constant(X[:,3:3])))

        ADCME.load(sess, "Data/$(nntype)/learned_nn$(idx)_ite$(Nite)_$(min_nn_init_ids[idx]).mat")

        out = run(sess, y)
        plot(X[NT*(sid-1)+1:NT*sid,1], out[NT*(sid-1)+1:NT*sid]/s_scale, "--$(cols[idx])", label="$(idx)-layer $(nnname)")

    end
    ylim((-0.1,0.6))
    xlim((-0.00125, 0.015))

    xlabel("Strain")
    ylabel("Stress (GPa)")
    legend()
    PyPlot.tight_layout()

    savefig("nnlearn_$(nntype)_p2p_strain_stress_tid$(tid).pdf")
    mpl.save("nnlearn_$(nntype)_p2p_strain_stress_tid$(tid).tex")
end

function run_nnlearn_test(tid::Int64, nntype::String, min_nn_init_ids::Array{Int64,1})

    #Data/$(nntype)/learned_nn$(idx)_ite$(i)_$(nn_init_id).mat

    for idx = 1:max_idx
        testcommand = `julia Test_NNTrussPull_Plasticity.jl $(nntype) $idx $(min_nn_init_ids[idx])`
        @show testcommand
        run(testcommand)
    end
end

function plot_nnlearn_fem_strain_stress(tid::Int64, nntype::String)
    close("all")

    @load "Data/domain$(tid).jld2" domain
    sid = 8
    strain = hcat(domain.history["strain"]...)
    stress = hcat(domain.history["stress"]...)
    plot(strain[sid,1:markevery:end], stress[sid,1:markevery:end]/s_scale, "o",fillstyle="none", label="Reference")

    for idx = 1:max_idx
        @load "Data/$(nntype)/nn$(idx)_domain_te$(tid).jld2" domain_te

        strain = hcat(domain_te.history["strain"]...)
        stress = hcat(domain_te.history["stress"]...)
        plot(strain[sid,1:markevery:end], stress[sid,1:markevery:end]/s_scale, "--$(cols[idx])", label="$(idx)-layer $(nnname)")
    end

    xlabel("Strain")
    ylabel("Stress (GPa)")
    legend()
    ylim((-0.1,0.6))
    xlim((-0.00125, 0.015))
    PyPlot.tight_layout()

    savefig("nnlearn_$(nntype)_fem_strain_stress$tid.pdf")
    mpl.save("nnlearn_$(nntype)_fem_strain_stress$tid.tex")
end

function plot_nnlearn_fem_disp(tid::Int64, nntype::String)
    close("all")
    @load "Data/domain$(tid).jld2" domain

    u1 = hcat(domain.history["state"]...)
    t1 = vcat([0.0], domain.history["time"]...)

    plot(t1[1:markevery:end]/t_scale, u1[5,1:markevery:end], "o",fillstyle="none", label="Reference")

    for idx = 1:max_idx
        @load "Data/$(nntype)/nn$(idx)_domain_te$(tid).jld2" domain_te
        u2 = hcat(domain_te.history["state"]...)
        t2 = vcat([0.0], domain_te.history["time"]...)

        @show idx, u2[5,:]
        plot(t2[1:markevery:end]/t_scale, u2[5,1:markevery:end], "--$(cols[idx])", label="$(idx)-layer $(nnname)")
    end
    xlabel("Time (s)")
    ylabel("Displacement (m)")
    legend()
    ylim((-0.005,0.06))
    PyPlot.tight_layout()
    savefig("nnlearn_$(nntype)_fem_disp$tid.pdf")
    mpl.save("nnlearn_$(nntype)_fem_disp$tid.tex")

end

tid = 3

min_nn_init_ids = plot_loss(nntype)

plot_nnlearn_p2p_strain_stress(tid, nntype, min_nn_init_ids)
run_nnlearn_test(tid, nntype, min_nn_init_ids)
plot_nnlearn_fem_strain_stress(tid, nntype)
plot_nnlearn_fem_disp(tid, nntype)


text:
export @skipif

immutable Skipped end

macro skipif(pred::Union{Expr, Bool, Symbol}, test::Expr)
  const name = test.args[1].args[1]
  quote
    if $(esc(pred))
      const $(esc(name)) = Skipped()
    else
      $(esc(test))
    end
  end
end

push_test!(tests::Vector{Tuple{AbstractString, Function}}, name::AbstractString, test::Skipped) = push!(tests, (name, () -> (true, :yellow, "$name SKIPPED")))


text:
# ziptextfiles.jl
#
#   Zip a text file or list of text files
#
# Jeremy Rüffer
# Thünen Institut
# Institut für Agrarklimaschutz
# Junior Research Group NITROSPHERE
# Julia 1.4.0
# 08.12.2016
# Last Edit: 01.04.2020

# Reference: https://zipfilejl.readthedocs.io/en/latest/

# include("/home/jeremy/AD-Home/Code/Julia/Jeremy/ziptextfiles.jl")

"# ziptextfiles(dest::String,files)

Zip a text file or a list of text files

`ziptextfiles(destination,files)`\n
* **destination**::String = Destination directory and file name of the resulting ZIP file
* **files**::Array{String,1} or String = File or list of files\n\n

---\n

#### Keywords:\n
* verbose::Bool = Display information as the function runs, TRUE is default
\n\n"
function ziptextfiles(dest::String,files::String;verbose::Bool=false)
	ziptextfiles(dest,[files],verbose=verbose)

	nothing
end

function ziptextfiles(dest::String,files::Array{String,1};verbose::Bool=false)
	verbose ? println(dest) : nothing

	# Zip each file
	tempFile = ""
	for i=1:1:length(files)
		tempFile = files[i]
		verbose ? println("\t" * tempFile) : nothing

		read(`$exe7z a -tzip $dest $tempFile`,String);
	end

	verbose ? println("Complete") : nothing
	nothing
end

function zipExtractAll(source::String,destinationFolder::String)
	##############
	##  Checks  ##
	##############
	isfile(source) ? nothing : error("ZIP file does not exist")
	isdir(destinationFolder) ? nothing : error("Destination folder does not exist")

	read(`$exe7z x $source -y -o$destinationFolder`,String);
end

function zipExtractFile(source::String,destinationFolder::String,fileToExtract::String)
	##############
	##  Checks  ##
	##############
	isfile(source) ? nothing : error("ZIP file does not exist")
	isdir(destinationFolder) ? nothing : error("Destination folder does not exist")

	read(`$exe7z x $source -y -o$destinationFolder $fileToExtract`,String);
end

function zipList(file::String)
	isfile(file) ? nothing : error("zipList input must be a file string.") # Check

	# Load the raw list and parse it into lines
	listRaw = read(`$exe7z l $file`,String);
	list = readlines(IOBuffer(listRaw))

	# Find where the file list starts
	fileListStart = 0
	for i=1:1:length(list)
		if occursin("Date",list[i]) & occursin("Time",list[i]) & occursin("Attr",list[i]) & occursin("Size",list[i]) & occursin("Compressed",list[i]) & occursin("Name",list[i])
			fileListStart = i + 2
		end
	end

	# List just the file names themselves
	fileList = Array{String}(undef,length(list)- 1 - fileListStart)
	for i=1:length(fileList)
		# Isolate the file name from the whole line
		fileList[i] = list[fileListStart + i - 1][findlast(' ',list[fileListStart + i - 1])+1:end]
	end

	return fileList
end

max_stars_repo_name: JeremyRueffer/ClimateDataIO

max_issues_repo_name: JeremyRueffer/ClimateDataIO


text:
rand()               # Defaults to n = 1 random number

rand()               # Different output each time

rand(3)              # Array of 3 random numbers

using Random
Random.seed!(1234)      # Seed number 1234
println(rand(3))
println(rand(5))
Random.seed!(1234)      # Re-seed with same number - will give the same sequence of random numbers
println(rand(2))
println(rand(6))

y = 2rand(5) .- 1     # 5 uniform random numbers on [-1,1]

z = 7rand(3) .+ 4    # 3 uniform random numbers on [4,11]

rand(1:10, 3)        # 3 random integers between 1 and 10

rand([π, exp(1), -1.0], 5)

randn()

μ = 2
σ = 5
σ*randn() + μ


text:
# These caves seem to be lava tubes. Parts are even still volcanically active;
# small hydrothermal vents release smoke into the caves that slowly settles
# like rain.
#
# If you can model how the smoke flows through the caves, you might be able to
# avoid it and be that much safer. The submarine generates a heightmap of the
# floor of the nearby caves for you (your puzzle input).
#
# Smoke flows to the lowest point of the area it's in. For example, consider
# the following heightmap:
#
# 2199943210
# 3987894921
# 9856789892
# 8767896789
# 9899965678
#
# Each number corresponds to the height of a particular location, where 9 is
# the highest and 0 is the lowest a location can be.
#
# Your first goal is to find the low points - the locations that are lower than
# any of its adjacent locations. Most locations have four adjacent locations
# (up, down, left, and right); locations on the edge or corner of the map have
# three or two adjacent locations, respectively. (Diagonal locations do not
# count as adjacent.)
#
# In the above example, there are four low points, all highlighted: two are in
# the first row (a 1 and a 0), one is in the third row (a 5), and one is in the
# bottom row (also a 5). All other locations on the heightmap have some lower
# adjacent location, and so are not low points.
#
# The risk level of a low point is 1 plus its height. In the above example, the
# risk levels of the low points are 2, 1, 6, and 6. The sum of the risk levels
# of all low points in the heightmap is therefore 15.
#
# Find all of the low points on your heightmap. What is the sum of the risk
# levels of all low points on your heightmap?

# read input
const width = 100
const height = 100
field = Array{UInt8}(undef, height, width)
foreach(eachline("input.txt"), Iterators.countfrom()) do l,y
  field[y,:] = parse.(UInt8, split(l, ""))
end

# find local minima
result = 0
for y in 1:height, x in 1:width
  v = field[y, x]
  ((x > 1 && v >= field[y, x-1]) ||
    (x < width && v >= field[y, x+1]) ||
    (y > 1 && v >= field[y-1, x]) ||
    (y < height && v >= field[y+1, x])) &&
      continue
  global result += 1 + v
end

println("Result: ", result)


text:
# test

s = "Hello World"
println(s)

abstract type Person   # type abstrait
end

abstract type Musician <: Person    # "sub-type"
end

mutable struct Rockstar <: Musician   # type concret non const (mutable)
	name::String
	instrument::String
	bandName::String
	headbandColor::String
	instrumentsPlayed::Int
end

struct ClassicMusician <: Musician   # type concret "const"
	name::String
	instrument::String
end

mutable struct Physicist <: Person
	name::String
	sleepHours::Float64
	favouriteLanguage::String
end

aure = Physicist("Aurelio", 6, "Julia")
aure.sleepHours = 8

function introduceMe(person::Person)
    println("Hello, my name is $(person.name).")
end

introduceMe(aure)

function f(x)
    2*x*x
end

print(f(2))

max_stars_repo_name: rboman/progs

max_issues_repo_name: rboman/progs

max_forks_repo_name: rboman/progs


text:
# Load Julia packages (libraries) needed.

using StatisticalRethinking
using CSV, DataFrames
using StanSample, MonteCarloMeasurements

ProjDir = @__DIR__

df = DataFrame!(CSV.File(joinpath(ProjDir, "sports.csv"), delim=','));
scale!(df, [:TM, :WUE, :FL])
println()

m2 = "
data {
 int < lower = 1 > N; // Sample size
 vector[N] WUE; // Outcome
 vector[N] TM; // Predictor
 vector[N] FL; // Predictor
}

parameters {
 real a; // Intercept
 real bFL; // Slope (regression coefficients)
 real bTM; // Slope (regression coefficients)
 real < lower = 0 > sigma;    // Error SD
}

model {
  vector[N] mu;               // mu is a vector
  a ~ normal(0, 0.2);         //Priors
  bFL ~ normal(0, 0.5);
  bTM ~ normal(0, 0.5);
  sigma ~ exponential(1);
  mu = a + bFL * FL + bTM * TM;
  WUE ~ normal(mu , sigma);     // Likelihood
}
";

# Define the SampleModel and set the output format to :mcmcchains.

m2s = SampleModel("m2", m2);

# Input data for cmdstan

m2_data = Dict("N" => size(df, 1), "FL" => df[:, :FL_s],
    "TM" => df[:, :TM_s], "WUE" => df[:, :WUE_s]);

# Sample using StanSample

rc = stan_sample(m2s, data=m2_data);

if success(rc)

  # Describe the draws

  dfa2 = read_samples(m2s; output_format=:dataframe)

  p2 = Particles(dfa2)

end


text:
# This is part of Quiz.jl. License is MIT.

# Moodle Engine

"""
`abstract type MOO <: qEngine end`\n
Moodle quiz engine.
"""
abstract type MOO <: qEngine end

export MOO


text:
module HDF5Arrays

using HDF5: DataFile, HDF5Dataset, h5open, d_create, datatype, dataspace
using TiledIteration: TileIterator

import Base: size, getindex, setindex!, iterate, to_indices, firstindex, lastindex
import HDF5: get_chunk

export HDF5Array

struct HDF5Array{T, N, D, C<:Union{Nothing, NTuple{N, Int}}} <: AbstractArray{T, N}
    data::D
    chunk::C
end

const HDF5Vector{T, D, C} = HDF5Array{T, 1, D, C}

size(arr::HDF5Array) = size(arr.data)

getindex(arr::HDF5Array{T}, I::Integer...) where {T} = getindex(arr.data, I...)::T

getindex(arr::HDF5Array{T}, I::Union{Colon, Integer, AbstractRange}...) where {T} =
    getindex(arr.data, I...)::Array{T}

setindex!(arr::HDF5Array, v, I::Union{Colon, Integer, AbstractRange}...) =
    setindex!(arr.data, v, I...)

# Read in one operation, rather than iterating index by index
to_indices(arr::HDF5Array, I::Tuple{<:CartesianIndices}) = I[1].indices

function iterate(arr::HDF5Array{T, N, HDF5Dataset, NTuple{N, Int}}) where {T, N}
    # init chunk axes
    cax_itr = TileIterator(axes(arr), arr.chunk)

    @debug "Chunk axes itr: $(cax_itr)"

    nextcax = iterate(cax_itr)
    isnothing(nextcax) && return nothing
    cax, cax_s = nextcax

    @debug "First chunk axes: $(cax)"

    buf::Array{T, N} = arr[cax...]

    @debug "First buffer: $(buf)"

    nextb = iterate(buf)
    isnothing(nextb) && return nothing

    b, b_s = nextb

    return b, ((buf, b_s), (cax_itr, cax_s))
end

function iterate(arr::HDF5Array{T, N, HDF5Dataset, NTuple{N, Int}}, state) where {T, N}
    (buf, b_s), (cax_itr, cax_s) = state

    nextb = iterate(buf, b_s)

    if isnothing(nextb)
        # get next buffer
        nextcax = iterate(cax_itr, cax_s)
        isnothing(nextcax) && return nothing
        cax, cax_s = nextcax

        @debug "Next chunk axes: $(cax)"

        buf::Array{T, N} = arr[cax...]

        @debug "Next buffer: $(buf)"

        nextb = iterate(buf)
    end

    b, b_s = nextb

    return b, ((buf, b_s), (cax_itr, cax_s))
end

firstindex(arr::HDF5Array) = first(eachindex(arr))
lastindex(arr::HDF5Array) = last(eachindex(arr))

include("construct.jl")

end


text:
function run_wide_tree_agnostic_sysid(n_leaves::Int64; value_aware=true)
    rng = MersenneTwister(1)
    widetree = WideTreeEnv(n_leaves, "real", rng)
    agent = WideTreeAgnosticSysIDAgent(widetree, rng)
    run(agent, rng; value_aware=value_aware)
end


text:
using Glob
import Oceananigans.Fields: set!

using Oceananigans.Fields: offset_data

"""
    Checkpointer{I, T, P} <: AbstractOutputWriter

An output writer for checkpointing models to a JLD2 file from which models can be restored.
"""
mutable struct Checkpointer{T, P} <: AbstractOutputWriter
      schedule :: T
           dir :: String
        prefix :: String
    properties :: P
         force :: Bool
       verbose :: Bool
       cleanup :: Bool
end

"""
    Checkpointer(model; schedule,
                        dir = ".",
                     prefix = "checkpoint",
                      force = false,
                    verbose = false,
                    cleanup = false,
                 properties = [:architecture, :boundary_conditions, :grid, :clock, :coriolis,
                               :buoyancy, :closure, :velocities, :tracers, :timestepper]
                 )

Construct a `Checkpointer` that checkpoints the model to a JLD2 file every so often as
specified by `schedule`. The `model.clock.iteration` is included
in the filename to distinguish between multiple checkpoint files.

Note that extra model `properties` can be safely specified, but removing crucial properties
such as `:velocities` will make restoring from the checkpoint impossible.

The checkpoint file is generated by serializing model properties to JLD2. However,
functions cannot be serialized to disk with JLD2. So if a model property
contains a reference somewhere in its hierarchy it will not be included in the checkpoint
file (and you will have to manually restore them).

Keyword arguments
=================
- `schedule` (required): Schedule that determines when to checkpoint.

- `dir`: Directory to save output to. Default: "." (current working directory).

- `prefix`: Descriptive filename prefixed to all output files. Default: "checkpoint".

- `force`: Remove existing files if their filenames conflict. Default: `false`.

- `verbose`: Log what the output writer is doing with statistics on compute/write times
             and file sizes. Default: `false`.

- `cleanup`: Previous checkpoint files will be deleted once a new checkpoint file is written.
             Default: `false`.

- `properties`: List of model properties to checkpoint. Some are required.
"""
function Checkpointer(model; schedule,
                             dir = ".",
                          prefix = "checkpoint",
                           force = false,
                         verbose = false,
                         cleanup = false,
                      properties = [:architecture, :grid, :clock, :coriolis,
                                    :buoyancy, :closure, :velocities, :tracers,
                                    :timestepper, :particles]
                     )

    # Certain properties are required for `restore_from_checkpoint` to work.
    required_properties = (:grid, :architecture, :velocities, :tracers, :timestepper, :particles)

    for rp in required_properties
        if rp ∉ properties
            @warn "$rp is required for checkpointing. It will be added to checkpointed properties"
            push!(properties, rp)
        end
    end

    for p in properties
        p isa Symbol || error("Property $p to be checkpointed must be a Symbol.")
        p ∉ propertynames(model) && error("Cannot checkpoint $p, it is not a model property!")

        if (p ∉ required_properties) && has_reference(Function, getproperty(model, p))
            @warn "model.$p contains a function somewhere in its hierarchy and will not be checkpointed."
            filter!(e -> e != p, properties)
        end
    end

    mkpath(dir)

    return Checkpointer(schedule, dir, prefix, properties, force, verbose, cleanup)
end

#####
##### Checkpointer utils
#####

""" Returns the full prefix (the `superprefix`) associated with `checkpointer`. """
checkpoint_superprefix(prefix) = prefix * "_iteration"

"""
    checkpoint_path(iteration::Int, c::Checkpointer)

Returns the path to the `c`heckpointer file associated with model `iteration`.
"""
checkpoint_path(iteration::Int, c::Checkpointer) =
    joinpath(c.dir, string(checkpoint_superprefix(c.prefix), iteration, ".jld2"))

# This is the default name used in the simulation.output_writers ordered dict.
defaultname(::Checkpointer, nelems) = :checkpointer

""" Returns `filepath`. Shortcut for `run!(simulation, pickup=filepath)`. """
checkpoint_path(filepath::AbstractString, checkpointers) = filepath

function checkpoint_path(pickup, checkpointers)
    length(checkpointers) == 0 && error("No checkpointers found: cannot pickup simulation!")
    length(checkpointers) > 1 && error("Multiple checkpointers found: not sure which one to pickup simulation from!")
    return checkpoint_path(pickup, first(checkpointers))
end

"""
    checkpoint_path(pickup::Bool, checkpointer)

For `pickup=true`, parse the filenames in `checkpointer.dir` associated with
`checkpointer.prefix` and return the path to the file whose name contains
the largest iteration.
"""
function checkpoint_path(pickup::Bool, checkpointer::Checkpointer)
    filepaths = glob(checkpoint_superprefix(checkpointer.prefix) * "*.jld2", checkpointer.dir)

    if length(filepaths) == 0 # no checkpoint files found
        return nothing
    else
        return latest_checkpoint(checkpointer, filepaths)
    end
end

function latest_checkpoint(checkpointer, filepaths)
    filenames = basename.(filepaths)
    leading = length(checkpoint_superprefix(checkpointer.prefix))
    trailing = 5 # length(".jld2")
    iterations = map(name -> parse(Int, name[leading+1:end-trailing]), filenames)
    latest_iteration, idx = findmax(iterations)
    return filepaths[idx]
end

#####
##### Writing checkpoints
#####

function write_output!(c::Checkpointer, model)
    filepath = checkpoint_path(model.clock.iteration, c)
    c.verbose && @info "Checkpointing to file $filepath..."

    t1 = time_ns()
    jldopen(filepath, "w") do file
        file["checkpointed_properties"] = c.properties
        serializeproperties!(file, model, c.properties)
    end

    t2, sz = time_ns(), filesize(filepath)
    c.verbose && @info "Checkpointing done: time=$(prettytime((t2-t1)/1e9)), size=$(pretty_filesize(sz))"

    @show c.cleanup
    c.cleanup && cleanup_checkpoints(c)

    return nothing
end

function cleanup_checkpoints(checkpointer)
    filepaths = glob(checkpoint_superprefix(checkpointer.prefix) * "*.jld2", checkpointer.dir)
    latest_checkpoint_filepath = latest_checkpoint(checkpointer, filepaths)
    @show filepaths
    @show latest_checkpoint_filepath
    [rm(filepath) for filepath in filepaths if filepath != latest_checkpoint_filepath]
    return nothing
end

#####
##### set! for checkpointer filepaths
#####

"""
    set!(model, filepath::AbstractString)

Set data in `model.velocities`, `model.tracers`, `model.timestepper.Gⁿ`, and
`model.timestepper.G⁻` to checkpointed data stored at `filepath`.
"""
function set!(model, filepath::AbstractString)

    jldopen(filepath, "r") do file

        # Validate the grid
        checkpointed_grid = file["grid"]
        model.grid == checkpointed_grid ||
            error("The grid associated with $filepath and model.grid are not the same!")

        # Set model fields and tendency fields
        model_fields = merge(model.velocities, model.tracers)

        for name in propertynames(model_fields)
            # Load data for each model field
            address = name ∈ (:u, :v, :w) ? "velocities/$name" : "tracers/$name"
            parent_data = file[address * "/data"]

            model_field = model_fields[name]
            copyto!(model_field.data.parent, parent_data)

            # Load tendency data
            #
            # Note: this step is unecessary for models that use RungeKutta3TimeStepper and
            # tendency restoration could be depcrecated in the future.

            # Tendency "n"
            parent_data = file["timestepper/Gⁿ/$name/data"]

            tendencyⁿ_field = model.timestepper.Gⁿ[name]
            copyto!(tendencyⁿ_field.data.parent, parent_data)

            # Tendency "n-1"
            parent_data = file["timestepper/G⁻/$name/data"]

            tendency⁻_field = model.timestepper.G⁻[name]
            copyto!(tendency⁻_field.data.parent, parent_data)
        end

        if !isnothing(model.particles)
            copyto!(model.particles.properties, file["particles"])
        end

        checkpointed_clock = file["clock"]

        # Update model clock
        model.clock.iteration = checkpointed_clock.iteration
        model.clock.time = checkpointed_clock.time

    end

    return nothing
end

#####
##### Nuke!
#####

function restore_if_not_missing(file, address)
    if !ismissing(file[address])
        return file[address]
    else
        @warn "Checkpoint file does not contain $address. Returning missing. " *
              "You might need to restore $address manually."
        return missing
    end
end

function restore_field(file, address, arch, grid, loc, kwargs)
    field_address = file[address * "/location"]
    data = offset_data(convert_to_arch(arch, file[address * "/data"]), grid, loc)

    # Extract field name from address. We use 2:end so "tracers/T "gets extracted
    # as :T while "timestepper/Gⁿ/T" gets extracted as :Gⁿ/T (we don't want to
    # apply the same BCs on T and T tendencies).
    field_name = split(address, "/")[2:end] |> join |> Symbol

    # If the user specified a non-default boundary condition through the kwargs
    # in restore_from_checkpoint, then use them when restoring the field. Otherwise
    # restore the BCs from the checkpoint file as long as they're not missing.
    if :boundary_conditions in keys(kwargs) && field_name in keys(kwargs[:boundary_conditions])
        bcs = kwargs[:boundary_conditions][field_name]
    else
        bcs = restore_if_not_missing(file, address * "/boundary_conditions")
    end

    return Field(field_address, arch, grid, bcs, data)
end

const u_location = (Face, Center, Center)
const v_location = (Center, Face, Center)
const w_location = (Center, Center, Face)
const c_location = (Center, Center, Center)

"""
    restore_from_checkpoint(filepath; kwargs=Dict())

Restore a model from the checkpoint file stored at `filepath`. `kwargs` can be passed
to the model constructor, which can be especially useful if you need to manually
restore forcing functions or boundary conditions that rely on functions.
"""
function restore_from_checkpoint(filepath; kwargs...)
    kwargs = length(kwargs) == 0 ? Dict{Symbol,Any}() : Dict{Symbol,Any}(kwargs)

    file = jldopen(filepath, "r")
    cps = file["checkpointed_properties"]

    if haskey(kwargs, :architecture)
        arch = kwargs[:architecture]
        filter!(p -> p ≠ :architecture, cps)
    else
        arch = file["architecture"]
    end

    grid = file["grid"]

    # Restore velocity fields
    kwargs[:velocities] = (u = restore_field(file, "velocities/u", arch, grid, u_location, kwargs),
                           v = restore_field(file, "velocities/v", arch, grid, v_location, kwargs),
                           w = restore_field(file, "velocities/w", arch, grid, w_location, kwargs))

    filter!(p -> p ≠ :velocities, cps) # pop :velocities from checkpointed properties

    # Restore tracer fields
    tracer_names         = Tuple(Symbol.(keys(file["tracers"])))
    tracer_fields        = Tuple(restore_field(file, "tracers/$c", arch, grid, c_location, kwargs) for c in tracer_names)
    tracer_fields_kwargs = NamedTuple{tracer_names}(tracer_fields)

    kwargs[:tracers] = TracerFields(tracer_names, arch, grid; tracer_fields_kwargs...)

    filter!(p -> p ≠ :tracers, cps) # pop :tracers from checkpointed properties

    # Restore time stepper tendency fields
    field_names = (:u, :v, :w, tracer_names...) # field names
    locs = (u_location, v_location, w_location, Tuple(c_location for c in tracer_names)...) # name locations

    G⁻_fields = Tuple(restore_field(file, "timestepper/G⁻/$(field_names[i])", arch, grid, locs[i], kwargs) for i = 1:length(field_names))
    Gⁿ_fields = Tuple(restore_field(file, "timestepper/Gⁿ/$(field_names[i])", arch, grid, locs[i], kwargs) for i = 1:length(field_names))

    G⁻_tendency_field_kwargs = NamedTuple{field_names}(G⁻_fields)
    Gⁿ_tendency_field_kwargs = NamedTuple{field_names}(Gⁿ_fields)

    # Restore time stepper
    kwargs[:timestepper] =
        QuasiAdamsBashforth2TimeStepper(arch, grid, tracer_names;
                                        G⁻ = TendencyFields(arch, grid, tracer_names; G⁻_tendency_field_kwargs...),
                                        Gⁿ = TendencyFields(arch, grid, tracer_names; Gⁿ_tendency_field_kwargs...))

    filter!(p -> p ≠ :timestepper, cps) # pop :timestepper from checkpointed properties

    # Restore the remaining checkpointed properties
    for p in cps
        if !haskey(kwargs, p)
            kwargs[p] = file["$p"]
        end
    end

    close(file)

    model = IncompressibleModel(; kwargs...)

    return model
end


text:
include("annulus.jl")
using GLMakie

lw = 3
Nθ = 100
θ = collect(range(0, 2π, length = Nθ))
xp = cos.(θ)
yp = sin.(θ)
fig = Figure(resolution = (1450, 1350))
ax = Axis(fig[1, 1], xlabel = "x", ylabel = "y", ylabelsize = 32,
    xlabelsize = 32, xgridstyle = :dash, ygridstyle = :dash, xtickalign = 1,
    xticksize = 30, ytickalign = 1, yticksize = 30, xlabelpadding = -10,
    xticklabelsize = 30, yticklabelsize = 30)
lines!(ax, xp .* r2, yp .* r2, color = :black, linewidth = lw)
lines!(ax, xp .* r1, yp .* r1, color = :black, linewidth = lw)

θs = collect(0:nn) * 2π / nn
for θ in θs
    ray_x = [0.8, 1] * cos(π / 2 + θ + φ)
    ray_y = [0.8, 1] * sin(π / 2 + θ + φ)
    lines!(ax, ray_x, ray_y, color = :black, linewidth = lw)
    ray_x = [0.8, 1] * cos(π / 2 + θ - φ)
    ray_y = [0.8, 1] * sin(π / 2 + θ - φ)
    lines!(ax, ray_x, ray_y, color = :black, linewidth = lw)
end

scatter!(ax, x_positions[2:end], y_positions[2:end], color = :red, markersize = 30)
scatter!(ax, x_positions[1:1], y_positions[1:1], color = :red, markersize = 60, marker = :star5)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
using Distributions, StatsBase, HypothesisTests, Plots, Random; pyplot()
Random.seed!(0)

n = 25
N = 10^4
xGrid = -10:0.001:10
kGrid = 0:0.01:5
dist1, dist2 = Exponential(1), Normal()

function ksStat(dist)
    data = rand(dist,n)
    Fhat = ecdf(data)
    sqrt(n)*maximum(abs.(Fhat.(xGrid) - cdf.(dist,xGrid)))
end

kStats1 = [ksStat(dist1) for _ in 1:N]
kStats2 = [ksStat(dist2) for _ in 1:N]

p1 = stephist(kStats1, bins=50,
	c=:blue, label="KS stat (Exponential)", normed=true)
p1 = plot!(kGrid, pdf.(Kolmogorov(),kGrid),
	c=:red, label="Kolmogorov PDF", xlabel="K", ylabel="Density")

p2 = stephist(kStats2, bins=50,
	c=:blue, label="KS stat (Normal)", normed=true)
p2 = plot!(kGrid, pdf.(Kolmogorov(),kGrid),
	c=:red, label="Kolmogorov PDF", xlabel="K", ylabel="Density")

plot(p1, p2, xlims=(0,2.5), ylims=(0,1.8), size=(800, 400))

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
@testset "tagpages" begin
    fs()
    write(joinpath("_layout", "head.html"), "")
    write(joinpath("_layout", "foot.html"), "")
    write(joinpath("_layout", "page_foot.html"), "")
    write("config.md", "")
    write("index.md", """
        Hello
        """)
    F.def_GLOBAL_VARS!()
    F.def_LOCAL_VARS!()
    write("pg1.md", "@def title = \"Page1\"")
    write("pg2.md", "@def title = \"Page2\"")
    F.set_var!(F.GLOBAL_VARS, "fd_page_tags", F.DTAG(("pg1" => Set(["aa", "bb"]),)))
    F.globvar("fd_page_tags")["pg2"] = Set(["bb", "cc"])

    @test Set(keys(F.globvar("fd_page_tags"))) == Set(["pg1", "pg2"])
    @test Set(union(values(F.globvar("fd_page_tags"))...)) == Set(["aa", "bb", "cc"])

    F.generate_tag_pages()

    @test F.globvar("fd_tag_pages")["aa"] == ["pg1"]
    @test F.globvar("fd_tag_pages")["bb"] == ["pg1","pg2"]
    @test F.globvar("fd_tag_pages")["cc"] == ["pg2"]

    @test isdir(F.path(:tag))
    @test isfile(joinpath(F.path(:tag), "aa", "index.html"))
    @test isfile(joinpath(F.path(:tag), "bb", "index.html"))
    @test isfile(joinpath(F.path(:tag), "cc", "index.html"))

    tagbb = read(joinpath(F.path(:tag), "bb", "index.html"), String)
    tagcc = read(joinpath(F.path(:tag), "cc", "index.html"), String)

    @test occursin("""

Page1

Page2

""", tagbb)
    @test occursin("""

Page2

""", tagcc)

    F.clear_dicts()
end

# ======= INTEGRATION ============

@testset "tags" begin
    fs()
    write(joinpath("_layout", "head.html"), "")
    write(joinpath("_layout", "foot.html"), "")
    write(joinpath("_layout", "page_foot.html"), "")
    write("config.md", "")
    write("index.md", """
        Hello
        """)
    F.def_GLOBAL_VARS!()
    F.def_LOCAL_VARS!()
    isdir("blog") && rm("blog", recursive=true)
    mkdir("blog")
    write(joinpath("blog", "pg1.md"), """
        @def tags = ["aa", "bb", "cc"]
        @def date = Date(2000, 01, 01)
        @def title = "Page 1"
        """)
    write(joinpath("blog", "pg2.md"), """
        @def tags = ["bb", "cc"]
        @def date = Date(2001, 01, 01)
        @def title = "Page 2"
        """)
    write(joinpath("blog", "pg3.md"), """
        @def tags = ["bb", "cc", "ee", "dd"]
        @def date = Date(2002, 01, 01)
        @def title = "Page 3"
        """)
    write(joinpath("blog", "pg4.md"), """
        @def tags = ["aa", "dd", "ee"]
        @def date = Date(2003, 01, 01)
        @def title = "Page 4"
        """)
    serve(clear=true, single=true, cleanup=false, nomess=true)
    @test isdir(joinpath("__site", "tag"))
    for tag in ("aa", "bb", "cc", "dd", "ee")
        local p
        p = joinpath("__site", "tag", tag, "index.html")
        @test isfile(p)
    end

    p = joinpath("__site", "tag", "aa", "index.html")
    c = read(p, String)
    @test occursin("""Tag: aa""", c)
    @test occursin("""
""", c)
    @test occursin("""

Page 4
""", c)
    @test occursin("""
Page 1

""", c)

    # Now remove some pages; we use bash commands so that they're blocking
    success(`rm $(joinpath("blog", "pg4.md"))`)
    success(`rm $(joinpath("blog", "pg3.md"))`)

    serve(clear=true, single=true, nomess=true)
    c = read(p, String)
    @test occursin("""
""", c)
    @test occursin("""

Page 1

""", c)

    F.clear_dicts()
end


text:
module Bucket

using ..Ahorn, Maple

displayName = "Bucket"
group = "Brushes"

drawingLayers = Ahorn.Layer[]

toolsLayer = nothing
targetLayer = nothing
material = '0'

bucketPosition = nothing
bucketBrush = nothing

function drawFill(x::Number, y::Number, tiles::Maple.Tiles, layer::Ahorn.Layer)
    h, w = size(tiles.data)

    if 1 <=x <= w && 1 <= y <= h
        pixels, ox, oy = Ahorn.shrinkMatrix(Ahorn.findFill(tiles.data, x, y))
        global bucketBrush = Ahorn.Brush("Bucket", pixels, (ox, oy))

        Ahorn.drawBrush(bucketBrush, layer, 1, 1)
    end
end

function drawBucket(layer::Ahorn.Layer, room::Ahorn.Room)
    if bucketPosition != nothing
        x, y = bucketPosition
        tiles = Ahorn.roomTiles(targetLayer, room)

        drawFill(x, y, tiles, layer)
    end
end

function applyFill!(x::Number, y::Number, layer::Ahorn.Layer, material::Char)
    if bucketBrush !== nothing
        Ahorn.History.addSnapshot!(Ahorn.History.RoomSnapshot("Bucket ($material)", Ahorn.loadedState.room))

        Maple.updateTileSize!(Ahorn.loadedState.room, '0', '0')
        tiles = Ahorn.roomTiles(layer, Ahorn.loadedState.room)
        Ahorn.applyBrush!(bucketBrush, tiles, material, 1, 1)
    end
end

function cleanup()
    global bucketPosition = nothing
    global bucketBrush = nothing

    Ahorn.redrawLayer!(toolsLayer)
end

function setMaterials!(layer::Ahorn.Layer)
    validTiles = Ahorn.validTiles(layer)
    tileNames = Ahorn.tileNames(layer)

    Ahorn.setMaterialList!([tileNames[mat] for mat in validTiles], row -> row[1] == get(tileNames, material, nothing))
end

function toolSelected(subTools::Ahorn.ListContainer, layers::Ahorn.ListContainer, materials::Ahorn.ListContainer)
    layerName = Ahorn.layerName(targetLayer)
    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_$(layerName)_material", tileNames["Air"])[1]

    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")
    Ahorn.updateLayerList!(["fgTiles", "bgTiles"], row -> row[1] == wantedLayer)

    Ahorn.redrawingFuncs["tools"] = drawBucket
    Ahorn.redrawLayer!(toolsLayer)
end

function layerSelected(list::Ahorn.ListContainer, materials::Ahorn.ListContainer, selected::String)
    global targetLayer = Ahorn.getLayerByName(drawingLayers, selected)
    Ahorn.persistence["brushes_layer"] = selected

    tileNames = Ahorn.tileNames(targetLayer)
    global material = get(Ahorn.persistence, "brushes_material_$(selected)", tileNames["Air"])[1]
    setMaterials!(targetLayer)
end

function materialSelected(list::Ahorn.ListContainer, selected::String)
    tileNames = Ahorn.tileNames(targetLayer)
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = tileNames[selected]
    global material = tileNames[selected]
end

function materialFiltered(list::Ahorn.ListContainer)
    tileNames = Ahorn.tileNames(targetLayer)
    Ahorn.selectRow!(list, row -> row[1] == get(tileNames, material, '0'))
end

function layersChanged(layers::Array{Ahorn.Layer, 1})
    wantedLayer = get(Ahorn.persistence, "brushes_layer", "fgTiles")

    global drawingLayers = layers
    global toolsLayer = Ahorn.getLayerByName(layers, "tools")
    global targetLayer = Ahorn.selectLayer!(layers, wantedLayer, "fgTiles")
end

function mouseMotion(x::Number, y::Number)
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    h, w = size(tiles.data)

    if bucketBrush === nothing || !get(bucketBrush.pixels, (y - bucketBrush.offset[2] + 1, x - bucketBrush.offset[1] + 1), false)
        global bucketPosition = (x, y)

        Ahorn.redrawLayer!(toolsLayer)
    end

    if x < 1 || x > w || y < 1 || y > h
        global bucketBrush = nothing
    end
end

function middleClick(x::Number, y::Number)
    tiles = Ahorn.roomTiles(targetLayer, Ahorn.loadedState.room)
    tileNames = Ahorn.tileNames(targetLayer)
    target = get(tiles.data, (y, x), '0')

    global material = target
    layerName = Ahorn.layerName(targetLayer)
    Ahorn.persistence["brushes_material_$(layerName)"] = material
    Ahorn.selectMaterialList!(tileNames[target])
end

function leftClick(x::Number, y::Number)
    applyFill!(x, y, targetLayer, material)

    Ahorn.redrawLayer!(targetLayer)
end

end


text:
# ---
# title: 643. Maximum Average Subarray I
# id: problem643
# author: Indigo
# date: 2021-02-04
# difficulty: Easy
# categories: Array
# link: <https://leetcode.com/problems/maximum-average-subarray-i/description/>
# hidden: true
# ---
#
# Given an array consisting of `n` integers, find the contiguous subarray of
# given length `k` that has the maximum average value. And you need to output
# the maximum average value.
#
# **Example 1:**
#
#
#
#     Input: [1,12,-5,-6,50,3], k = 4
#     Output: 12.75
#     Explanation: Maximum average is (12-5-6+50)/4 = 51/4 = 12.75
#
#
#
#
# **Note:**
#
#   1. 1 <= `k` <= `n` <= 30,000.
#   2. Elements of the given array will be in the range [-10,000, 10,000].
#
#
#
#
## @lc code=start
using LeetCode

function find_max_average(nums::Vector{Int}, K::Int)
    res = sum(@view(nums[1:K]))
    tmp = res
    i, j = 1, K + 1
    while j ≤ length(nums)
        tmp += nums[j] - nums[i]
        res = max(tmp, res)
        i += 1
        j += 1
    end
    res / K
end
## @lc code=end


text:
include("default_consts.jl")

##########################################################################################
# Node defaults
##########################################################################################

#=
function set_node_defaults(osmplot)
    osm = osmplot.osm[]
    gpk = osmplot.graphplotkwargs

    node_color = @lift(get($gpk, :node_color, :default))
    node_size = @lift(get($gpk, :node_size, :default))
    nlabels = @lift(show_nlabels($gpk, $(osmplot.hide_nlabels), osm))
    nlabels_textsize = @lift(get($gpk, :nlabels_textsize, :default))

    return (; node_color, node_size, nlabels, nlabels_textsize)
end
=#

#=
It would be nicer to just use the function above but NamedTuples are automatically converted
to Attributes, which are not a subtype of AbstractObservable. Hence we cannot lift from
them. Finding a workaround to this would allow for better reactivity of the plot attributes.
So far everything I've tried was unsuccessful. Maybe at some point I will be smarter or the
underlying system will have changed.

Anyways, the following workaround has to suffice for now although it only results in
static plots. The user has to recreate the plots to reflect any changes of graphplotkwargs.
=#

function set_node_defaults(osmplot, edge_width, edge_color)
    osm = osmplot.osm[]
    gpk = osmplot.graphplotkwargs
    sorted_edges = osmplot.sorted_edges[]

    node_size, maxs = size_nodes(gpk, osm, sorted_edges, edge_width)
    node_color = color_nodes(gpk, osm, edge_color, maxs, sorted_edges, edge_width)
    nlabels = @lift(show_nlabels(gpk, $(osmplot.hide_nlabels), osm))
    nlabels_textsize = haskey(gpk, :nlabels_textsize) ? gpk.nlabels_textsize : 9

    return (; node_color, node_size, nlabels, nlabels_textsize)
end

function size_nodes(gpk, osm, sorted_edges, edge_width)
    if haskey(gpk, :node_size)
        return fill(gpk.node_size[], Graphs.nv(osm.graph)), nothing
    else
        if edge_width isa Vector
            sizes = fill(0, Graphs.nv(osm.graph))
            maxs = fill(0, Graphs.nv(osm.graph))
            for i in eachindex(sizes)
                # get indices for each connected edge of current vertex
                edge_indices = findall(x -> first(x) == i, sorted_edges)
                if !isempty(edge_indices)
                    # set vertex size to max width of connected edges and save its edge index
                    sizes[i], m = findmax(edge_width[ei] for ei in edge_indices)
                    maxs[i] = edge_indices[m]
                end
            end
            return sizes, maxs
        else
            return edge_width, nothing
        end
    end
end

function color_nodes(gpk, osm, edge_color, maxs, sorted_edges, edge_width)
    if haskey(gpk, :node_color)
        return fill(gpk.node_color[], Graphs.nv(osm.graph))
    else
        if edge_color isa Vector
            colors = fill(colorant"#444", Graphs.nv(osm.graph))
            # set vertex color to color of connected edge with max width
            if !isnothing(maxs) # use maxs if previously filled in size_nodes
                for i in eachindex(maxs)
                    if !iszero(maxs[i])
                        colors[i] = edge_color[maxs[i]]
                    end
                end
            else # similar procedure to size_nodes
                for i in eachindex(colors)
                    edge_indices = findall(x -> first(x) == i, sorted_edges)
                    if !isempty(edge_indices)
                        _, m = findmax(edge_width[ei] for ei in edge_indices)
                        colors[i] = edge_color[m]
                    end
                end
            end
            return colors
        else
            return edge_color
        end
    end
end

function show_nlabels(gpk, hide_nlabels, osm)
    labels = if hide_nlabels
        nothing
    else
        haskey(gpk, :nlabels) ? gpk.nlabels[] : repr.(keys(osm.nodes))
    end

    return labels
end

##########################################################################################
# Edge defaults
##########################################################################################

function set_edge_defaults(osmplot)
    osm = osmplot.osm[]
    gpk = osmplot.graphplotkwargs
    i2w = osmplot.index_to_way[]
    sorted_edges = osmplot.sorted_edges[]
    n2i = osm.node_to_index
    ways = osm.highways

    edge_color = color_streets(gpk, i2w)
    edge_width = width_streets(gpk, i2w)
    elabels = @lift(show_elabels(gpk, $(osmplot.hide_elabels), sorted_edges, n2i, ways))
    elabels_textsize = haskey(gpk, :elabels_textsize) ? gpk.elabels_textsize : 11
    arrow_attr = (;
        markersize = arrows_streets(gpk, i2w, edge_width),
        color = :white
    )

    return (; edge_color, edge_width, elabels, elabels_textsize, arrow_attr)
end

function color_streets(gpk, i2w)
    if haskey(gpk, :edge_color)
        return fill(gpk.edge_color[], length(i2w))
    else
        colors = fill(colorant"#444", length(i2w))

        for (index, way) in pairs(i2w)
            waytype = get(way.tags, "highway", nothing)
            if !isnothing(waytype) && haskey(WAYTYPECOLORS, waytype)
                colors[index] = WAYTYPECOLORS[waytype]
            end
        end

        return colors
    end
end

function width_streets(gpk, i2w)
    if haskey(gpk, :edge_width)
        return fill(gpk.edge_width[],  length(i2w))
    else
        widths = fill(BASEWIDTH, length(i2w))

        for (index, way) in pairs(i2w)
            waytype = get(way.tags, "highway", nothing)
            if !isnothing(waytype) && haskey(WAYTYPEWIDTHS, waytype)
                widths[index] *= WAYTYPEWIDTHS[waytype]
            end
        end

        return widths
    end
end

function label_streets(sorted_edges, n2i, ways)
    labels = fill("", length(sorted_edges))

    # replace edge labels but try to avoid repetitions
    for way in values(ways)
        if haskey(way.tags, "name")
            node_indices = begin
                l = length(way.nodes)
                if l == 2 # pick first two nodes of a way
                    n1 = 1
                    n2 = 2
                elseif iseven(l) # pick middle node and the next one
                    n1 = l ÷ 2
                    n2 = n1 + 1
                else # pick the two nodes next to the middle
                    n1 = floor(Int, l / 2)
                    n2 = ceil(Int, l / 2)
                end
                [n2i[way.nodes[n1]], n2i[way.nodes[n2]]]
            end
            i = only(indexin([node_indices], sorted_edges))
            labels[i] = way.tags["name"]
        end
    end

    return labels
end

function show_elabels(gpk, hide_elabels, sorted_edges, n2i, ways)
    labels = if hide_elabels
        nothing
    else
        haskey(gpk, :elabels) ? gpk.elabels[] : label_streets(sorted_edges, n2i, ways)
    end
    return labels
end

function arrows_streets(gpk, i2w, edge_width)
    if haskey(gpk, :arrow_attr) && haskey(gpk, :arrow_size)
        return gpk.arrow_attr.arrow_size[]
    else
        markersizes = fill(0, length(i2w))

        for (index, way) in pairs(i2w)
            if get(way.tags, "oneway", false)
                markersizes[index] = edge_width[index]
            end
        end

        return markersizes
    end
end

##########################################################################################
# Building defaults
##########################################################################################

function get_building_polys(buildings)
    building_polys = Polygon[]

    for b in values(buildings)
        pntsvec = Vector{Point2f}[]
        for bp in b.polygons
            pnts = Point2f[(node.location.lon, node.location.lat) for node in bp.nodes]
            push!(pntsvec, pnts)
        end
        push!(building_polys, Polygon(pntsvec[1], pntsvec[2:end]))
    end

    return building_polys
end


text:
module TestPkg

# __precompile__(false)

using TaylorIntegration
# using BenchmarkTools
using InteractiveUtils

order = 25
N = 200
x0 = 10randn(2N)
μ = 1e-7rand(N)
x0[N] = - sum(x0[1:N-1].*μ[1:N-1])/μ[N]
x0[2N] = - sum(x0[N+1:2N-1].*μ[1:N-1])/μ[N]

@taylorize xdot1(x, p, t) = b1-x^2

@taylorize function f1!(dq, q, params, t)
    local N = Int(length(q)/2)
    local _eltype_q_ = eltype(q)
    local μ = params
    X = Array{_eltype_q_}(undef, N, N)
    accX = Array{_eltype_q_}(undef, N) #acceleration
    for j in 1:N
        accX[j] = zero(q[1])
        dq[j] = q[N+j]
    end
    #compute accelerations
    for j in 1:N
        for i in 1:N
            if i == j
            else
                X[i,j] = q[i]-q[j]
                temp_001 = accX[j] + (μ[i]*X[i,j])
                accX[j] = temp_001
            end #if i != j
        end #for, i
    end #for, j
    for i in 1:N
        dq[N+i] = accX[i]
    end
    nothing
end

function f2!(dq, q, params, t)
    local N = Int(length(q)/2)
    local _eltype_q_ = eltype(q)
    local μ = params
    X = Array{_eltype_q_}(undef, N, N)
    accX = Array{_eltype_q_}(undef, N) #acceleration
    for j in 1:N
        accX[j] = zero(q[1])
        dq[j] = q[N+j]
    end
    #compute accelerations
    for j in 1:N
        for i in 1:N
            if i == j
            else
                X[i,j] = q[i]-q[j]
                temp_001 = accX[j] + (μ[i]*X[i,j])
                accX[j] = temp_001
            end #if i != j
        end #for, i
    end #for, j
    for i in 1:N
        dq[N+i] = accX[i]
    end
    nothing
end

function TaylorIntegration.jetcoeffs!(::Val{f2!}, t::Taylor1{_T}, q::AbstractVector{Taylor1{_S}}, dq::AbstractVector{Taylor1{_S}}, params) where {_T <: Real, _S <: Number}
    order = t.order
    local N = Int(length(q) / 2)
    local _eltype_q_ = eltype(q)
    local μ = params
    X = Array{_eltype_q_}(undef, N, N)
    accX = Array{_eltype_q_}(undef, N)
    for j = 1:N
        accX[j] = Taylor1(zero(constant_term(q[1])), order)
        dq[j] = Taylor1(identity(constant_term(q[N + j])), order)
    end
    tmp337 = Array{Taylor1{_S}}(undef, size(X))
    tmp337 .= Taylor1(zero(_S), order)
    temp_001 = Array{Taylor1{_S}}(undef, size(tmp337))
    temp_001 .= Taylor1(zero(_S), order)
    for j = 1:N
        for i = 1:N
            if i == j
            else
                X[i, j] = Taylor1(constant_term(q[i]) - constant_term(q[j]), order)
                tmp337[i, j] = Taylor1(constant_term(μ[i]) * constant_term(X[i, j]), order)
                temp_001[i, j] = Taylor1(constant_term(accX[j]) + constant_term(tmp337[i, j]), order)
                accX[j] = Taylor1(identity(constant_term(temp_001[i, j])), order)
            end
        end
    end
    for i = 1:N
        dq[N + i] = Taylor1(identity(constant_term(accX[i])), order)
    end
    for __idx = eachindex(q)
        (q[__idx]).coeffs[2] = (dq[__idx]).coeffs[1]
    end
    for ord = 1:order - 1
        ordnext = ord + 1
        for j = 1:N
            TaylorSeries.zero!(accX[j], q[1], ord)
            TaylorSeries.identity!(dq[j], q[N + j], ord)
        end
        for j = 1:N
            for i = 1:N
                if i == j
                else
                    TaylorSeries.subst!(X[i, j], q[i], q[j], ord)
                    TaylorSeries.mul!(tmp337[i, j], μ[i], X[i, j], ord)
                    TaylorSeries.add!(temp_001[i, j], accX[j], tmp337[i, j], ord)
                    TaylorSeries.identity!(accX[j], temp_001[i, j], ord)
                end
            end
        end
        for i = 1:N
            TaylorSeries.identity!(dq[N + i], accX[i], ord)
        end
        for __idx = eachindex(q)
            (q[__idx]).coeffs[ordnext + 1] = (dq[__idx]).coeffs[ordnext] / ordnext
        end
    end
    return nothing
end

ex1 = :(function f3!(dq, q, params, t)
    local N = Int(length(q)/2)
    local _eltype_q_ = eltype(q)
    local μ = params
    X = Array{_eltype_q_}(undef, N, N)
    accX = Array{_eltype_q_}(undef, N) #acceleration
    for j in 1:N
        accX[j] = zero(q[1])
        dq[j] = q[N+j]
    end
    #compute accelerations
    for j in 1:N
        for i in 1:N
            if i == j
            else
                X[i,j] = q[i]-q[j]
                temp_001 = accX[j] + (μ[i]*X[i,j])
                accX[j] = temp_001
            end #if i != j
        end #for, i
    end #for, j
    for i in 1:N
        dq[N+i] = accX[i]
    end
    nothing
end)

ex2 = :(xdot2(x, p, t) = b1-x^2)

ex3 = :(function harm_osc!(dx, x, p, t)
    local ω = p[1]
    local ω2 = ω^2
    dx[1] = x[2]
    dx[2] = - (ω2 * x[1])
    return nothing
end)

nex1, narr1 = TaylorIntegration._make_parsed_jetcoeffs(ex1)
nex2, narr2 = TaylorIntegration._make_parsed_jetcoeffs(ex2)
nex3, narr3 = TaylorIntegration._make_parsed_jetcoeffs(ex3)

greet(f, parse_eqs) = begin
    t = Taylor1(order)
    x = Taylor1.(x0, t.order)
    dx = similar(x)

    @show @which TaylorIntegration.__jetcoeffs!(Val(parse_eqs), f, t, x, dx, similar(x), μ)
    @show @which TaylorIntegration.jetcoeffs!(Val(f), t, x, dx, μ)
    @show methods(TaylorIntegration.jetcoeffs!)
    # @btime TaylorIntegration.__jetcoeffs!(Val($parse_eqs), $f, $t, $x, $dx, similar($x), μ)
end

end # module


text:
"""
    AlleeExtinction <: CellRule

    AlleeExtinction(minfounders)
    AlleeExtinction{R}(; minfounders=5.0)
    AlleeExtinction{R,W}(minfounders)

Causes extinction in a cell when a population is below a minimum number of individuals.

- `minfounders`: minimum founding individuals required for a viable population.
    Must be a type that can be compared to the grid values using `isless`.

Pass grid `Symbol`s to `R` or both `R` and `W` type parameters to use to specific grids.
"""
struct AlleeExtinction{R,W,MF} <: CellRule{R,W}
    minfounders::MF
end
AlleeExtinction{R,W}(; minfounders=Param(5.0, bounds=(1.0, 200.0))) where {R,W} =
    AlleeExtinction{R,W}(minfounders)

@inline function applyrule(data, rule::AlleeExtinction, N, I)
    f = get(data, rule.minfounders, I...)

    return N >= f ? N : zero(N)
end


text:
using OrdinaryDiffEq, DiffEqDevTools, ParameterizedFunctions

f = @ode_def LotkaVolterra begin
  dx = a*x - b*x*y
  dy = -c*y + d*x*y
end a b c d
p = [1.5,1.0,3.0,1.0]
prob = ODEProblem(f,[1.0;1.0],(0.0,10.0),p)

abstols = 1 ./ 10 .^ (6:13)
reltols = 1 ./ 10 .^ (3:10);
sol = solve(prob,Vern7(),abstol=1/10^14,reltol=1/10^14)
test_sol = TestSolution(sol)

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=10000)
using Plots; gr()
DIAGRAMS["Lotka Volterra Low Order"] = plot(wp);

setups = [Dict(:alg=>DP5())
          Dict(:alg=>Tsit5())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,numruns=1000,maxiters=10000,error_estimate=:L2,dense_errors=true)
DIAGRAMS["Lotka Volterra Low Order Interpolation"] = plot(wp);

setups = [Dict(:alg=>DP8())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,save_everystep=false,numruns=1000,maxiters=1000)
DIAGRAMS["Lotka Volterra High Order"] = plot(wp)

setups = [Dict(:alg=>DP8())
          Dict(:alg=>Vern7())
          Dict(:alg=>Vern8())
          Dict(:alg=>Vern6())
]
wp = WorkPrecisionSet(prob,abstols,reltols,setups;appxsol=test_sol,dense=true,numruns=1000,maxiters=1000,error_estimate=:L2)
DIAGRAMS["Lotka Volterra High Order Interpolation"] = plot(wp)


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using ViewerGL
GL = ViewerGL
using LinearAlgebra

# input of primitive shapes
V,(VV,EV,FV) = Lar.simplex(2, true)
W,(WW,EW,FW) = Lar.cuboid([1,1], true)
triangle = (V,EV)
square = (W,EW)
# hierarchical assembly
model = Lar.Struct([
            Lar.Struct([
                square,
                Lar.t(.175,.175), Lar.s(.5,.5),
                triangle
            ]),
            Lar.t(.825,.825), Lar.s(-.5,-.5),
            triangle
        ])
# visualization of generated wire-frame model
V,EV = Lar.struct2lar(model)
GL.VIEW(GL.numbering(.4)((V,[VV, EV])));

triangles = Lar.triangulate2d(V, EV)
GL.VIEW([ GL.GLGrid(V,triangles) ]);

# generate edges
ev = map(sort,cat([[[u,v], [v,w], [w,u]] for (u,v,w) in triangles]))
# remove duplicated edges from triangulation
nodupev = collect(Set(ev))

GL.VIEW(GL.numbering(.3)((V,[VV, nodupev])));
GL.VIEW([ GL.GLGrid(V,triangles) ]);
GL.VIEW(GL.numbering(.25)((V,[VV, nodupev, triangles]),GL.COLORS[1],0.2));


text:
"""
`Scale <: ImageTransformation`

Description
============

Multiplies the image height and image width by individually specified
constant factors. This means that the size of the output image
depends on the size of the input image. If one wants to resize
each dimension by the same factor, use :class:`Zoom` instead.

Usage
======

    Scale()

    Scale(width = 1, height = 1)

Arguments
==========

- **`width`** : The constant factor that should be used to scale the
width of the output image

- **`height`** : The constant factor that should be used to scale the
height of the output image

Methods
========

- **`transform`** : Applies the transformation to the given Image
or set of images

Author(s)
==========

- Christof Stocker (Github: https://github.com/Evizero)

Examples
========

    using Augmentor
    using TestImages

    # load an example image
    img = testimage("lena")

    # Scale the image to half its width and twice its height
    img_new = transform(Scale(0.5, 2), img)

see also
=========

`ImageTransformation`, `ProbableOperation`, `transform`
"""
immutable Scale <: ImageTransformation
    width::Float64
    height::Float64

    function Scale(width::Real, height::Real)
        width > 0 || throw(ArgumentError("width has to be greater than 0"))
        height > 0 || throw(ArgumentError("height has to be greater than 0"))
        new(Float64(width), Float64(height))
    end
end

Scale(; width = 1, height = 1) = Scale(width, height)

Base.showcompact(io::IO, op::Scale) = print(io, "Scale by $(op.width) x $(op.height).")
multiplier(::Scale) = 1

function transform{T<:AbstractImage}(op::Scale, img::T)
    s1, s2 = size(img)
    if isxfirst(img)
        x = ceil(Int, s1 * op.width)
        y = ceil(Int, s2 * op.height)
        result = copyproperties(img, Images.imresize(img, (x, y)))
        _log_operation!(result, op)
    else
        x = ceil(Int, s2 * op.width)
        y = ceil(Int, s1 * op.height)
        result = copyproperties(img, Images.imresize(img, (y, x)))
        _log_operation!(result, op)
    end
end