text:
module Javis

using FFMPEG
using LaTeXStrings
using LightXML
import Luxor
import Luxor: Point, @layer
using Random

const FRAMES_SYMBOL = [:same]

"""
    Video

Defines the video canvas for an animation.

# Fields
- `width::Int` the width in pixel
- `height::Int` the height in pixel
- `defs::Dict{Symbol, Any}` Some definitions which should be accessible throughout the video.
"""
mutable struct Video
    width::Int
    height::Int
    defs::Dict{Symbol,Any}
end

"""
    CURRENT_VIDEO

holds the current video in an array to be declared as a constant
The current video can be accessed using CURRENT_VIDEO[1]
"""
const CURRENT_VIDEO = Array{Video,1}()

"""
    Video(width, height)

Create a video with a certain `width` and `height` in pixel.
This also sets `CURRENT_VIDEO`.
"""
function Video(width, height)
    video = Video(width, height, Dict{Symbol,Any}())
    if isempty(CURRENT_VIDEO)
        push!(CURRENT_VIDEO, video)
    else
        CURRENT_VIDEO[1] = video
    end
    return video
end

"""
    Rel

Ability to define frames in a relative fashion.

# Example
```
 Action(1:100, ground; in_global_layer=true),
 Action(1:90, :red_ball, (args...)->circ(p1, "red"), Rotation(from_rot, to_rot)),
 Action(Rel(10), :blue_ball, (args...)->circ(p2, "blue"), Rotation(2π, from_rot, :red_ball)),
 Action((video, args...)->path!(path_of_red, pos(:red_ball), "red"))
```
is the same as
```
Action(1:100, ground; in_global_layer=true),
Action(1:90, :red_ball, (args...)->circ(p1, "red"), Rotation(from_rot, to_rot)),
Action(91:100, :blue_ball, (args...)->circ(p2, "blue"), Rotation(2π, from_rot, :red_ball)),
Action(91:100, (video, args...)->path!(path_of_red, pos(:red_ball), "red"))
```

# Fields
- rel::UnitRange defines the frames in a relative fashion.
"""
struct Rel
    rel::UnitRange
end

"""
    Rel(i::Int)

Shorthand for Rel(1:i)
"""
Rel(i::Int) = Rel(1:i)

"""
    Frames

Stores the actual computed frames and the user input
which can be i.e. `:same` or `Rel(10)`.
The `frames` are computed in `javis`.
"""
mutable struct Frames{T}
    frames::Union{Nothing,UnitRange}
    user::T
end

Base.convert(::Type{Frames}, x::Union{Symbol,Rel}) = Frames(nothing, x)
Base.convert(::Type{Frames}, x::UnitRange) = Frames(x, x)

"""
    Transformation

Defines a transformation which can be returned by an action to be accessible later.
See the `circ` function inside the [`javis`](@ref) as an example.

It can be accessed by another [`Action`])(@ref) using the symbol notation
like `:red_ball` in the example.

# Fields
- `p::Point`: the translation part of the transformation
- `angle::Float64`: the angle component of the transformation (in radians)
"""
mutable struct Transformation
    p::Point
    angle::Float64
end

abstract type Transition end
abstract type InternalTransition end

abstract type AbstractAction end

"""
    SubAction <: AbstractAction

A SubAction can be used in the keyword arguments of an [`Action`](@ref) to define small
sub actions on the action function, such as [`appear`](@ref).

A SubAction should not be created by hand but instead by using one of the constructors.

# Fields
- `frames::Frames`: the frames relative to the parent [`Action`](@ref)
- `func::Function`: the function that gets called in each of those frames.
    Takes the following arguments: `video, action, subaction, rel_frame`
- `transitions::Vector{Transition}`: A list of transitions like [`Translation`](@ref)
- `internal_transitions::Vector{InternalTransition}`:
    A list of internal transitions which store the current transition for a specific frame.
"""
mutable struct SubAction <: AbstractAction
    frames::Frames
    func::Function
    transitions::Vector{Transition}
    internal_transitions::Vector{InternalTransition}
end

SubAction(transitions::Transition...) = SubAction(:same, transitions...)
SubAction(func::Function) = SubAction(:same, func)

"""
    SubAction(frames, func::Function)

A `SubAction` can be defined with frames and a function
inside the `subactions` kwarg of an [`Action`](@ref).
In the following example a filled circle with radius 50 appears in the first 20 frames,
which means the opacity is increased from 0 to 1.0.
Then it stays at full opacity and disappears the same way in the last 20 frames.

# Example
javis(demo, [
    BackgroundAction(1:100, ground),
    Action((args...)->circle(O, 50, :fill); subactions = [
        SubAction(1:20, appear(:fade)),
        SubAction(81:100, disappear(:fade))
    ])
])

# Arguments
- `frames`: A list of frames for which the function should be called.
    - The frame numbers are relative to the parent [`Action`](@ref).
- `func::Function`: The function that gets called for the frames.
    - Needs to have four arguments: `video, action, subaction, rel_frame`
    - For [`appear`](@ref) and [`disappear`](@ref) a closure exists,
      such that `appear(:fade)` works.
"""
SubAction(frames, func::Function) = SubAction(frames, func, [], [])

"""
    SubAction(frames, trans::Transition...)

A `SubAction` can also be defined this way with having a list of transitions.
This is similar to defining transitions inside [`Action`](@ref)

In the following example a circle is faded in during the first 25 frames then moves to
- `Point(100, 20)` then to `Point(120, -20)` (the translations are added)
- and then back to the origin
In the last 25 frames it disappears from the world.

# Example
```
javis(demo, [
        BackgroundAction(1:200, ground_opacity),
        Action((args...)->circle(O, 50, :fill); subactions = [
            SubAction(1:25, appear(:fade)),
            SubAction(26:75, Translation(Point(100, 20))),
            SubAction(76:100, Translation(Point(20, -40))),
            SubAction(101:175, Translation(Point(-120, 20))),
            SubAction(176:200, disappear(:fade))
        ]),
    ], tempdirectory="current/images", pathname="current/circle_square.gif")
```

# Arguments
- `frames`: A list of frames for which the function should be called.
    - The frame numbers are relative to the parent [`Action`](@ref).
- `trans::Transition...`: A list of transitions that shall be performed.
"""
SubAction(frames, trans::Transition...) =
    SubAction(frames, (args...) -> 1, collect(trans), [])

"""
    ActionSetting

The current settings of an [`Action`](@ref) which are saved in `action.current_setting`.

# Fields
- `line_width::Float64`: the current line width
- `mul_line_width::Float64`: the current multiplier for line width.
    The actual line width is then: `mul_line_width * line_width`
- `opacity::Float64`: the current opacity
- `mul_opacity::Float64`: the current multiplier for opacity.
    The actual opacity is then: `mul_opacity * opacity`
- `fontsize::Float64` the current font size
- `current_scale::Tuple{Float64, Float64}`: the current scale
- `desired_scale::Tuple{Float64, Float64}`: the new desired scale
- `mul_scale::Float64`: the multiplier for the new desired scale.
    The actual new scale is then: `mul_scale * desired_scale`
"""
mutable struct ActionSetting
    line_width::Float64
    mul_line_width::Float64 # the multiplier of line width is between 0 and 1
    opacity::Float64
    mul_opacity::Float64 # the multiplier of opacity is between 0 and 1
    fontsize::Float64
    # scale has three fields instead of just the normal two
    # current scale
    # desired scale and scale multiplier => `desired_scale*mul_scale` is the new desired scale
    # the scale change needs to be computed using `current_scale` and the desired scale
    current_scale::Tuple{Float64,Float64}
    desired_scale::Tuple{Float64,Float64}
    mul_scale::Float64 # the multiplier of scale is between 0 and 1
end

ActionSetting() = ActionSetting(1.0, 1.0, 1.0, 1.0, 10.0, (1.0, 1.0), (1.0, 1.0), 1.0)

"""
    update_ActionSetting!(as::ActionSetting, by::ActionSetting)

Set the fields of `as` to the same as `by`. Basically copying them over.
"""
function update_ActionSetting!(as::ActionSetting, by::ActionSetting)
    as.line_width = by.line_width
    as.mul_line_width = by.mul_line_width
    as.opacity = by.opacity
    as.mul_opacity = by.mul_opacity
    as.fontsize = by.fontsize
    as.current_scale = by.current_scale
    as.desired_scale = by.desired_scale
    as.mul_scale = by.mul_scale
end

"""
    Action

Defines what is drawn in a defined frame range.

# Fields
- `frames::Frames`: A range of frames for which the `Action` is called
- `id::Union{Nothing, Symbol}`: An id which can be used to save the result of `func`
- `func::Function`: The drawing function which draws something on the canvas.
    It gets called with the arguments `video, action, frame`
- `transitions::Vector{Transition}` a list of transitions
    which can be performed before the function gets called.
- `internal_transitions::Vector{InternalTransition}`:
    Similar to `transitions` but holds the concrete information whereas `Transition` can
    hold links to other actions which need to be computed first.
    See [`compute_transformation!`](@ref)
- `opts::Any` can hold any options defined by the user
"""
mutable struct Action <: AbstractAction
    frames::Frames
    id::Union{Nothing,Symbol}
    func::Function
    transitions::Vector{Transition}
    internal_transitions::Vector{InternalTransition}
    subactions::Vector{SubAction}
    current_setting::ActionSetting
    opts::Dict{Symbol,Any}
end

"""
    set_frames!(a::AbstractAction, last_frames::UnitRange)

Compute the frames based on a.frames and `last_frames`.
Save the result in `a.frames.frames` which can be accessed via [`get_frames`](@ref).
"""
function set_frames!(a::AbstractAction, last_frames::UnitRange)
    frames = a.frames.user
    a.frames.frames = get_frames(frames, last_frames)
end

"""
    get_frames(a::AbstractAction)

Return `a.frames.frames` which holds the computed frames for the AbstractAction `a`.
"""
get_frames(a::AbstractAction) = a.frames.frames

"""
    get_frames(frames::Symbol, last_frames::UnitRange)

Get the frames based on a symbol (currently only `same`) and the `last_frames`.
Throw `ArgumentError` if symbol is unknown
"""
function get_frames(frames::Symbol, last_frames::UnitRange)
    if frames === :same
        return last_frames
    else
        throw(ArgumentError("Currently the only symbol supported for defining frames is `:same`"))
    end
end

"""
    get_frames(frames::Rel, last_frames::UnitRange)

Return the frames based on a relative frames [`Rel`](@ref) object and the `last_frames`.
"""
function get_frames(frames::Rel, last_frames::UnitRange)
    start_frame = last(last_frames) + first(frames.rel)
    last_frame = last(last_frames) + last(frames.rel)
    return start_frame:last_frame
end

"""
    CURRENT_ACTION

holds the current action in an array to be declared as a constant
The current action can be accessed using CURRENT_ACTION[1]
"""
const CURRENT_ACTION = Array{Action,1}()

"""
    Action(frames, func::Function, args...)

The most simple form of an action (if there are no `args`/`kwargs`) just calls
`func(video, action, frame)` for each of the frames it is defined for.
`args` are defined it the next function definition and can be seen in action
    in this example [`javis`](@ref)
"""
Action(frames, func::Function, args...; kwargs...) =
    Action(frames, nothing, func, args...; kwargs...)

"""
    Action(frames_or_id::Symbol, func::Function, args...)

This function decides whether you wrote `Action(frames_symbol, ...)`,
    or `Action(id_symbol, ...)`
If the symbol `frames_or_id` is not a `FRAMES_SYMBOL` then it is used as an id_symbol.
"""
function Action(frames_or_id::Symbol, func::Function, args...; kwargs...)
    if frames_or_id in FRAMES_SYMBOL
        Action(frames_or_id, nothing, func, args...; kwargs...)
    else
        Action(:same, frames_or_id, func, args...; kwargs...)
    end
end

"""
    Action(func::Function, args...)

Similar to the above but uses the same as frames as the action above.
"""
Action(func::Function, args...; kwargs...) =
    Action(:same, nothing, func, args...; kwargs...)

"""
    Action(frames, id::Union{Nothing,Symbol}, func::Function,
           transitions::Transition...; kwargs...)

# Arguments
- `frames`: defines for which frames this action is called
- `id::Symbol`: Is used if the `func` returns something which
    shall be accessible by other actions later
- `func::Function` the function that is called after the `transitions` are performed
- `transitions::Transition...` a list of transitions that are performed before
    the function `func` itself is called

The keywords arguments will be saved inside `.opts` as a `Dict{Symbol, Any}`
"""
function Action(
    frames,
    id::Union{Nothing,Symbol},
    func::Function,
    transitions::Transition...;
    kwargs...,
)
    CURRENT_VIDEO[1].defs[:last_frames] = frames
    opts = Dict(kwargs...)
    subactions = SubAction[]
    if haskey(opts, :subactions)
        subactions = opts[:subactions]
        delete!(opts, :subactions)
    end
    Action(frames, id, func, collect(transitions), [], subactions, ActionSetting(), opts)
end

"""
    BackgroundAction(frames, func::Function, args...; kwargs...)

Create an Action where `in_global_layer` is set to true such that
i.e the specified color in the background is applied globally (basically a new default)
"""
function BackgroundAction(frames, func::Function, args...; kwargs...)
    Action(frames, nothing, func, args...; in_global_layer = true, kwargs...)
end

"""
    BackgroundAction(frames, id::Symbol, func::Function, args...; kwargs...)

Create an Action where `in_global_layer` is set to true and saves the return into `id`.
"""
function BackgroundAction(frames, id::Symbol, func::Function, args...; kwargs...)
    Action(frames, id, func, args...; in_global_layer = true, kwargs...)
end

mutable struct InternalTranslation <: InternalTransition
    by::Point
end

mutable struct InternalRotation <: InternalTransition
    angle::Float64
    center::Point
end

mutable struct InternalScaling <: InternalTransition
    scale::Tuple{Float64,Float64}
end

"""
    Translation <: Transition

Stores the `Point` or a link for the start and end position of the translation

# Fields
`from::Union{Point, Symbol}`: The start position or a link to the start position.
    See `:red_ball` in [`javis`](@ref)
`to::Union{Point, Symbol}`: The end position or a link to the end position
"""
struct Translation <: Transition
    from::Union{Point,Symbol}
    to::Union{Point,Symbol}
end

"""
    Translation(p::Union{Point, Symbol})

Create a `Translation(O, p)` such that a translation is done from the origin.
"""
Translation(p::Union{Point,Symbol}) = Translation(O, p)

"""
    Translation(x::Real, y::Real)

Create a `Translation(O, Point(x,y))` such that a translation is done from the origin.
Shorthand for writing `Translation(Point(x,y))`.
"""
Translation(x::Real, y::Real) = Translation(Point(x, y))

"""
    Rotation <: Transition

Stores the rotation similar to [`Translation`](@ref) with `from` and `to`
but also the rotation point.

# Fields
- `from::Union{Float64, Symbol}`: The start rotation or a link to it
- `to::Union{Float64, Symbol}`: The end rotation or a link to it
- `center::Union{Point, Symbol}`: The center of the rotation or a link to it.
"""
struct Rotation <: Transition
    from::Union{Float64,Symbol}
    to::Union{Float64,Symbol}
    center::Union{Point,Symbol}
end

"""
    Rotation(r::Union{Float64, Symbol})

Rotation as a transition from 0.0 to `r` .
Can be used as a short-hand.
"""
Rotation(r::Union{Float64,Symbol}) = Rotation(0.0, r)

"""
    Rotation(r::Union{Float64, Symbol}, center::Union{Point, Symbol})

Rotation as a transition from `0.0` to `r` around `center`.
Can be used as a short-hand for rotating around a `center` point.
"""
Rotation(r::Union{Float64,Symbol}, center::Union{Point,Symbol}) = Rotation(0.0, r, center)

"""
    Rotation(from, to)

Rotation as a transition from `from` to `to` (in radians) around the origin.
"""
Rotation(from, to) = Rotation(from, to, O)

"""
    Scaling <: Transition

Stores the scaling similar to [`Translation`](@ref) with `from` and `to`.

# Example
- Can be called with different constructors like:
```
Scaling(10) -> Scaling(CURRENT_SCALING, (10.0, 10.0))
Scaling(10, :my-scale) -> Scaling((10.0, 10.0), :my_scale)
Scaling(10, 2) -> Scaling((10.0, 10.0), (2.0, 2.0))
Scaling(10, (1,2)) -> Scaling((10.0, 10.0), (1.0, 2.0))
```

# Fields
- `from::Union{Tuple{Float64, Float64}, Symbol}`: The start scaling or a link to it
- `to::Union{Tuple{Float64, Float64}, Symbol}`: The end scaling or a link to it
- `compute_from_once::Bool`: Saves whether the from is computed for the first frame or
    every frame. Is true if from is `:_current_scale`.
"""
mutable struct Scaling <: Transition
    from::Union{Tuple{Float64,Float64},Symbol}
    to::Union{Tuple{Float64,Float64},Symbol}
    compute_from_once::Bool
end

Scaling(to::Tuple) = Scaling(:_current_scale, to, true)
Scaling(to::Real) = Scaling(:_current_scale, convert(Float64, to), true)
Scaling(to::Symbol) = Scaling(:_current_scale, to, true)

function Scaling(from::Real, to::Real, compute_from_once = false)
    from_flt = convert(Float64, from)
    to_flt = convert(Float64, to)
    Scaling((from_flt, from_flt), (to_flt, to_flt), compute_from_once)
end

function Scaling(from::Real, to, compute_from_once = false)
    flt = convert(Float64, from)
    Scaling((flt, flt), to, compute_from_once)
end

function Scaling(from, to::Real, compute_from_once = false)
    flt = convert(Float64, to)
    Scaling(from, (flt, flt), compute_from_once)
end

"""
    Line

A type to define a line by two points. Can be used i.e. in [`projection`](@ref)
We mean the mathematic definition of a continuous line and not a segment of a line.

# Fields
- `p1::Point`: start point
- `p2::Point`: second point to define the line
"""
struct Line
    p1::Point
    p2::Point
end

"""
    Base.:*(m::Array{Float64,2}, transformation::Transformation)

Convert the transformation to a matrix and multiplies m*trans_matrix.
Return a new Transformation
"""
function Base.:*(m::Array{Float64,2}, transformation::Transformation)
    θ = transformation.angle
    p = transformation.p
    trans_matrix = [
        cos(θ) -sin(θ) p.x
        sin(θ) cos(θ) p.y
        0 0 1
    ]
    res = m * trans_matrix
    return Transformation(Point(gettranslation(res)...), getrotation(res))
end

# cache such that creating svgs from LaTeX don't need to be created every time
# this is also used for test cases such that `tex2svg` doesn't need to be installed on Github Actions
const LaTeXSVG = Dict{LaTeXString,String}(
    L"\mathcal{O}(\log{n})" =>
        "<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"8.413ex\" height=\"2.843ex\" style=\"vertical-align: -0.838ex;\" viewBox=\"0 -863.1 3622.2 1223.9\" role=\"img\" focusable=\"false\" xmlns=\"http://www.w3.org/2000/svg\" aria-labelledby=\"MathJax-SVG-1-Title\">\n<title id=\"MathJax-SVG-1-Title\">script upper O left-parenthesis log n right-parenthesis</title>\n<defs aria-hidden=\"true\">\n<path stroke-width=\"1\" id=\"E1-MJCAL-4F\" d=\"M308 428Q289 428 289 438Q289 457 318 508T378 593Q417 638 475 671T599 705Q688 705 732 643T777 483Q777 380 733 285T620 123T464 18T293 -22Q188 -22 123 51T58 245Q58 327 87 403T159 533T249 626T333 685T388 705Q404 705 404 693Q404 674 363 649Q333 632 304 606T239 537T181 429T158 290Q158 179 214 114T364 48Q489 48 583 165T677 438Q677 473 670 505T648 568T601 617T528 636Q518 636 513 635Q486 629 460 600T419 544T392 490Q383 470 372 459Q341 430 308 428Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-28\" d=\"M94 250Q94 319 104 381T127 488T164 576T202 643T244 695T277 729T302 750H315H319Q333 750 333 741Q333 738 316 720T275 667T226 581T184 443T167 250T184 58T225 -81T274 -167T316 -220T333 -241Q333 -250 318 -250H315H302L274 -226Q180 -141 137 -14T94 250Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-6C\" d=\"M42 46H56Q95 46 103 60V68Q103 77 103 91T103 124T104 167T104 217T104 272T104 329Q104 366 104 407T104 482T104 542T103 586T103 603Q100 622 89 628T44 637H26V660Q26 683 28 683L38 684Q48 685 67 686T104 688Q121 689 141 690T171 693T182 694H185V379Q185 62 186 60Q190 52 198 49Q219 46 247 46H263V0H255L232 1Q209 2 183 2T145 3T107 3T57 1L34 0H26V46H42Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-6F\" d=\"M28 214Q28 309 93 378T250 448Q340 448 405 380T471 215Q471 120 407 55T250 -10Q153 -10 91 57T28 214ZM250 30Q372 30 372 193V225V250Q372 272 371 288T364 326T348 362T317 390T268 410Q263 411 252 411Q222 411 195 399Q152 377 139 338T126 246V226Q126 130 145 91Q177 30 250 30Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-67\" d=\"M329 409Q373 453 429 453Q459 453 472 434T485 396Q485 382 476 371T449 360Q416 360 412 390Q410 404 415 411Q415 412 416 414V415Q388 412 363 393Q355 388 355 386Q355 385 359 381T368 369T379 351T388 325T392 292Q392 230 343 187T222 143Q172 143 123 171Q112 153 112 133Q112 98 138 81Q147 75 155 75T227 73Q311 72 335 67Q396 58 431 26Q470 -13 470 -72Q470 -139 392 -175Q332 -206 250 -206Q167 -206 107 -175Q29 -140 29 -75Q29 -39 50 -15T92 18L103 24Q67 55 67 108Q67 155 96 193Q52 237 52 292Q52 355 102 398T223 442Q274 442 318 416L329 409ZM299 343Q294 371 273 387T221 404Q192 404 171 388T145 343Q142 326 142 292Q142 248 149 227T179 192Q196 182 222 182Q244 182 260 189T283 207T294 227T299 242Q302 258 302 292T299 343ZM403 -75Q403 -50 389 -34T348 -11T299 -2T245 0H218Q151 0 138 -6Q118 -15 107 -34T95 -74Q95 -84 101 -97T122 -127T170 -155T250 -167Q319 -167 361 -139T403 -75Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMATHI-6E\" d=\"M21 287Q22 293 24 303T36 341T56 388T89 425T135 442Q171 442 195 424T225 390T231 369Q231 367 232 367L243 378Q304 442 382 442Q436 442 469 415T503 336T465 179T427 52Q427 26 444 26Q450 26 453 27Q482 32 505 65T540 145Q542 153 560 153Q580 153 580 145Q580 144 576 130Q568 101 554 73T508 17T439 -10Q392 -10 371 17T350 73Q350 92 386 193T423 345Q423 404 379 404H374Q288 404 229 303L222 291L189 157Q156 26 151 16Q138 -11 108 -11Q95 -11 87 -5T76 7T74 17Q74 30 112 180T152 343Q153 348 153 366Q153 405 129 405Q91 405 66 305Q60 285 60 284Q58 278 41 278H27Q21 284 21 287Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-29\" d=\"M60 749L64 750Q69 750 74 750H86L114 726Q208 641 251 514T294 250Q294 182 284 119T261 12T224 -76T186 -143T145 -194T113 -227T90 -246Q87 -249 86 -250H74Q66 -250 63 -250T58 -247T55 -238Q56 -237 66 -225Q221 -64 221 250T66 725Q56 737 55 738Q55 746 60 749Z\"></path>\n</defs>\n<g stroke=\"currentColor\" fill=\"currentColor\" stroke-width=\"0\" transform=\"matrix(1 0 0 -1 0 0)\" aria-hidden=\"true\">\n <use xlink:href=\"#E1-MJCAL-4F\" x=\"0\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMAIN-28\" x=\"796\" y=\"0\"></use>\n<g transform=\"translate(1186,0)\">\n <use xlink:href=\"#E1-MJMAIN-6C\"></use>\n <use xlink:href=\"#E1-MJMAIN-6F\" x=\"278\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMAIN-67\" x=\"779\" y=\"0\"></use>\n</g>\n <use xlink:href=\"#E1-MJMATHI-6E\" x=\"2632\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMAIN-29\" x=\"3232\" y=\"0\"></use>\n</g>\n</svg>",
    L"\mathcal{O}\left(\frac{\log{x}}{2}\right)" =>
        "<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"11.183ex\" height=\"6.176ex\" style=\"vertical-align: -2.505ex;\" viewBox=\"0 -1580.7 4814.8 2659.1\" role=\"img\" focusable=\"false\" xmlns=\"http://www.w3.org/2000/svg\" aria-labelledby=\"MathJax-SVG-1-Title\">\n<title id=\"MathJax-SVG-1-Title\">script upper O left-parenthesis StartFraction log x Over 2 EndFraction right-parenthesis</title>\n<defs aria-hidden=\"true\">\n<path stroke-width=\"1\" id=\"E1-MJCAL-4F\" d=\"M308 428Q289 428 289 438Q289 457 318 508T378 593Q417 638 475 671T599 705Q688 705 732 643T777 483Q777 380 733 285T620 123T464 18T293 -22Q188 -22 123 51T58 245Q58 327 87 403T159 533T249 626T333 685T388 705Q404 705 404 693Q404 674 363 649Q333 632 304 606T239 537T181 429T158 290Q158 179 214 114T364 48Q489 48 583 165T677 438Q677 473 670 505T648 568T601 617T528 636Q518 636 513 635Q486 629 460 600T419 544T392 490Q383 470 372 459Q341 430 308 428Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-28\" d=\"M94 250Q94 319 104 381T127 488T164 576T202 643T244 695T277 729T302 750H315H319Q333 750 333 741Q333 738 316 720T275 667T226 581T184 443T167 250T184 58T225 -81T274 -167T316 -220T333 -241Q333 -250 318 -250H315H302L274 -226Q180 -141 137 -14T94 250Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-6C\" d=\"M42 46H56Q95 46 103 60V68Q103 77 103 91T103 124T104 167T104 217T104 272T104 329Q104 366 104 407T104 482T104 542T103 586T103 603Q100 622 89 628T44 637H26V660Q26 683 28 683L38 684Q48 685 67 686T104 688Q121 689 141 690T171 693T182 694H185V379Q185 62 186 60Q190 52 198 49Q219 46 247 46H263V0H255L232 1Q209 2 183 2T145 3T107 3T57 1L34 0H26V46H42Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-6F\" d=\"M28 214Q28 309 93 378T250 448Q340 448 405 380T471 215Q471 120 407 55T250 -10Q153 -10 91 57T28 214ZM250 30Q372 30 372 193V225V250Q372 272 371 288T364 326T348 362T317 390T268 410Q263 411 252 411Q222 411 195 399Q152 377 139 338T126 246V226Q126 130 145 91Q177 30 250 30Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-67\" d=\"M329 409Q373 453 429 453Q459 453 472 434T485 396Q485 382 476 371T449 360Q416 360 412 390Q410 404 415 411Q415 412 416 414V415Q388 412 363 393Q355 388 355 386Q355 385 359 381T368 369T379 351T388 325T392 292Q392 230 343 187T222 143Q172 143 123 171Q112 153 112 133Q112 98 138 81Q147 75 155 75T227 73Q311 72 335 67Q396 58 431 26Q470 -13 470 -72Q470 -139 392 -175Q332 -206 250 -206Q167 -206 107 -175Q29 -140 29 -75Q29 -39 50 -15T92 18L103 24Q67 55 67 108Q67 155 96 193Q52 237 52 292Q52 355 102 398T223 442Q274 442 318 416L329 409ZM299 343Q294 371 273 387T221 404Q192 404 171 388T145 343Q142 326 142 292Q142 248 149 227T179 192Q196 182 222 182Q244 182 260 189T283 207T294 227T299 242Q302 258 302 292T299 343ZM403 -75Q403 -50 389 -34T348 -11T299 -2T245 0H218Q151 0 138 -6Q118 -15 107 -34T95 -74Q95 -84 101 -97T122 -127T170 -155T250 -167Q319 -167 361 -139T403 -75Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMATHI-78\" d=\"M52 289Q59 331 106 386T222 442Q257 442 286 424T329 379Q371 442 430 442Q467 442 494 420T522 361Q522 332 508 314T481 292T458 288Q439 288 427 299T415 328Q415 374 465 391Q454 404 425 404Q412 404 406 402Q368 386 350 336Q290 115 290 78Q290 50 306 38T341 26Q378 26 414 59T463 140Q466 150 469 151T485 153H489Q504 153 504 145Q504 144 502 134Q486 77 440 33T333 -11Q263 -11 227 52Q186 -10 133 -10H127Q78 -10 57 16T35 71Q35 103 54 123T99 143Q142 143 142 101Q142 81 130 66T107 46T94 41L91 40Q91 39 97 36T113 29T132 26Q168 26 194 71Q203 87 217 139T245 247T261 313Q266 340 266 352Q266 380 251 392T217 404Q177 404 142 372T93 290Q91 281 88 280T72 278H58Q52 284 52 289Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-32\" d=\"M109 429Q82 429 66 447T50 491Q50 562 103 614T235 666Q326 666 387 610T449 465Q449 422 429 383T381 315T301 241Q265 210 201 149L142 93L218 92Q375 92 385 97Q392 99 409 186V189H449V186Q448 183 436 95T421 3V0H50V19V31Q50 38 56 46T86 81Q115 113 136 137Q145 147 170 174T204 211T233 244T261 278T284 308T305 340T320 369T333 401T340 431T343 464Q343 527 309 573T212 619Q179 619 154 602T119 569T109 550Q109 549 114 549Q132 549 151 535T170 489Q170 464 154 447T109 429Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-29\" d=\"M60 749L64 750Q69 750 74 750H86L114 726Q208 641 251 514T294 250Q294 182 284 119T261 12T224 -76T186 -143T145 -194T113 -227T90 -246Q87 -249 86 -250H74Q66 -250 63 -250T58 -247T55 -238Q56 -237 66 -225Q221 -64 221 250T66 725Q56 737 55 738Q55 746 60 749Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJSZ3-28\" d=\"M701 -940Q701 -943 695 -949H664Q662 -947 636 -922T591 -879T537 -818T475 -737T412 -636T350 -511T295 -362T250 -186T221 17T209 251Q209 962 573 1361Q596 1386 616 1405T649 1437T664 1450H695Q701 1444 701 1441Q701 1436 681 1415T629 1356T557 1261T476 1118T400 927T340 675T308 359Q306 321 306 250Q306 -139 400 -430T690 -924Q701 -936 701 -940Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJSZ3-29\" d=\"M34 1438Q34 1446 37 1448T50 1450H56H71Q73 1448 99 1423T144 1380T198 1319T260 1238T323 1137T385 1013T440 864T485 688T514 485T526 251Q526 134 519 53Q472 -519 162 -860Q139 -885 119 -904T86 -936T71 -949H56Q43 -949 39 -947T34 -937Q88 -883 140 -813Q428 -430 428 251Q428 453 402 628T338 922T245 1146T145 1309T46 1425Q44 1427 42 1429T39 1433T36 1436L34 1438Z\"></path>\n</defs>\n<g stroke=\"currentColor\" fill=\"currentColor\" stroke-width=\"0\" transform=\"matrix(1 0 0 -1 0 0)\" aria-hidden=\"true\">\n <use xlink:href=\"#E1-MJCAL-4F\" x=\"0\" y=\"0\"></use>\n<g transform=\"translate(963,0)\">\n <use xlink:href=\"#E1-MJSZ3-28\"></use>\n<g transform=\"translate(736,0)\">\n<g transform=\"translate(120,0)\">\n<rect stroke=\"none\" width=\"2138\" height=\"60\" x=\"0\" y=\"220\"></rect>\n<g transform=\"translate(60,726)\">\n <use xlink:href=\"#E1-MJMAIN-6C\"></use>\n <use xlink:href=\"#E1-MJMAIN-6F\" x=\"278\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMAIN-67\" x=\"779\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMATHI-78\" x=\"1446\" y=\"0\"></use>\n</g>\n <use xlink:href=\"#E1-MJMAIN-32\" x=\"819\" y=\"-687\"></use>\n</g>\n</g>\n <use xlink:href=\"#E1-MJSZ3-29\" x=\"3115\" y=\"-1\"></use>\n</g>\n</g>\n</svg>",
    L"E=mc^2" =>
        "<svg xmlns:xlink=\"http://www.w3.org/1999/xlink\" width=\"8.976ex\" height=\"2.676ex\" style=\"vertical-align: -0.338ex;\" viewBox=\"0 -1006.6 3864.5 1152.1\" role=\"img\" focusable=\"false\" xmlns=\"http://www.w3.org/2000/svg\" aria-labelledby=\"MathJax-SVG-1-Title\">\n<title id=\"MathJax-SVG-1-Title\">upper E equals m c squared</title>\n<defs aria-hidden=\"true\">\n<path stroke-width=\"1\" id=\"E1-MJMATHI-45\" d=\"M492 213Q472 213 472 226Q472 230 477 250T482 285Q482 316 461 323T364 330H312Q311 328 277 192T243 52Q243 48 254 48T334 46Q428 46 458 48T518 61Q567 77 599 117T670 248Q680 270 683 272Q690 274 698 274Q718 274 718 261Q613 7 608 2Q605 0 322 0H133Q31 0 31 11Q31 13 34 25Q38 41 42 43T65 46Q92 46 125 49Q139 52 144 61Q146 66 215 342T285 622Q285 629 281 629Q273 632 228 634H197Q191 640 191 642T193 659Q197 676 203 680H757Q764 676 764 669Q764 664 751 557T737 447Q735 440 717 440H705Q698 445 698 453L701 476Q704 500 704 528Q704 558 697 578T678 609T643 625T596 632T532 634H485Q397 633 392 631Q388 629 386 622Q385 619 355 499T324 377Q347 376 372 376H398Q464 376 489 391T534 472Q538 488 540 490T557 493Q562 493 565 493T570 492T572 491T574 487T577 483L544 351Q511 218 508 216Q505 213 492 213Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-3D\" d=\"M56 347Q56 360 70 367H707Q722 359 722 347Q722 336 708 328L390 327H72Q56 332 56 347ZM56 153Q56 168 72 173H708Q722 163 722 153Q722 140 707 133H70Q56 140 56 153Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMATHI-6D\" d=\"M21 287Q22 293 24 303T36 341T56 388T88 425T132 442T175 435T205 417T221 395T229 376L231 369Q231 367 232 367L243 378Q303 442 384 442Q401 442 415 440T441 433T460 423T475 411T485 398T493 385T497 373T500 364T502 357L510 367Q573 442 659 442Q713 442 746 415T780 336Q780 285 742 178T704 50Q705 36 709 31T724 26Q752 26 776 56T815 138Q818 149 821 151T837 153Q857 153 857 145Q857 144 853 130Q845 101 831 73T785 17T716 -10Q669 -10 648 17T627 73Q627 92 663 193T700 345Q700 404 656 404H651Q565 404 506 303L499 291L466 157Q433 26 428 16Q415 -11 385 -11Q372 -11 364 -4T353 8T350 18Q350 29 384 161L420 307Q423 322 423 345Q423 404 379 404H374Q288 404 229 303L222 291L189 157Q156 26 151 16Q138 -11 108 -11Q95 -11 87 -5T76 7T74 17Q74 30 112 181Q151 335 151 342Q154 357 154 369Q154 405 129 405Q107 405 92 377T69 316T57 280Q55 278 41 278H27Q21 284 21 287Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMATHI-63\" d=\"M34 159Q34 268 120 355T306 442Q362 442 394 418T427 355Q427 326 408 306T360 285Q341 285 330 295T319 325T330 359T352 380T366 386H367Q367 388 361 392T340 400T306 404Q276 404 249 390Q228 381 206 359Q162 315 142 235T121 119Q121 73 147 50Q169 26 205 26H209Q321 26 394 111Q403 121 406 121Q410 121 419 112T429 98T420 83T391 55T346 25T282 0T202 -11Q127 -11 81 37T34 159Z\"></path>\n<path stroke-width=\"1\" id=\"E1-MJMAIN-32\" d=\"M109 429Q82 429 66 447T50 491Q50 562 103 614T235 666Q326 666 387 610T449 465Q449 422 429 383T381 315T301 241Q265 210 201 149L142 93L218 92Q375 92 385 97Q392 99 409 186V189H449V186Q448 183 436 95T421 3V0H50V19V31Q50 38 56 46T86 81Q115 113 136 137Q145 147 170 174T204 211T233 244T261 278T284 308T305 340T320 369T333 401T340 431T343 464Q343 527 309 573T212 619Q179 619 154 602T119 569T109 550Q109 549 114 549Q132 549 151 535T170 489Q170 464 154 447T109 429Z\"></path>\n</defs>\n<g stroke=\"currentColor\" fill=\"currentColor\" stroke-width=\"0\" transform=\"matrix(1 0 0 -1 0 0)\" aria-hidden=\"true\">\n <use xlink:href=\"#E1-MJMATHI-45\" x=\"0\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMAIN-3D\" x=\"1042\" y=\"0\"></use>\n <use xlink:href=\"#E1-MJMATHI-6D\" x=\"2098\" y=\"0\"></use>\n<g transform=\"translate(2977,0)\">\n <use xlink:href=\"#E1-MJMATHI-63\" x=\"0\" y=\"0\"></use>\n <use transform=\"scale(0.707)\" xlink:href=\"#E1-MJMAIN-32\" x=\"613\" y=\"583\"></use>\n</g>\n</g>\n</svg>\n",
)

include("util.jl")
include("luxor_overrides.jl")
include("backgrounds.jl")
include("svg2luxor.jl")
include("morphs.jl")
include("subaction_animations.jl")

latex(text::LaTeXString) = latex(text, O)
latex(text::LaTeXString, pos::Point) = latex(text, pos, :stroke)
latex(text::LaTeXString, x, y) = latex(text, Point(x, y), :stroke)
@deprecate latex(text::LaTeXString, fsize::Real) begin
    fontsize(fsize)
    latex(text)
end

"""
    latex(text::LaTeXString, pos::Point, action::Symbol)

Add the latex string `text` to the top left corner of the LaTeX path.
Can be added to `Luxor.jl` graphics via [`Video`](@ref).

**NOTES:**
- **This only works if `tex2svg` is installed.**
    It can be installed using the following command (you may have to prefix this command with `sudo` depending on your installation):

        npm install -g mathjax-node-cli

- **The `latex` method must be called from within an [`Action`](@ref).**

# Arguments
- `text::LaTeXString`: a LaTeX string to render.
- `pos::Point`: position of the upper left corner of the latex text. Default: `O`
    - can be written as `x, y` instead of `Point(x, y)`
- `action::Symbol`: graphics actions defined by `Luxor.jl`. Default `:stroke`.
Available actions:
  - `:stroke` - Draws the latex string on the canvas. For more info check `Luxor.strokepath`
  - `:path` - Creates the path of the latex string but does not render it to the canvas.

# Throws
- `IOError`: mathjax-node-cli is not installed

# Example

```
using Javis
using LaTeXStrings

function ground(args...)
    background("white")
    sethue("black")
end

function draw_latex(video, action, frame)
    x = 100
    y = 120
    latex(L"\\sqrt{5}", x, y)
end

demo = Video(500, 500)
javis(demo, [BackgroundAction(1:2, ground), Action(draw_latex)],
      pathname = "latex.gif")
```

"""
function latex(text::LaTeXString, pos::Point, action::Symbol)
    translate(pos)
    # check if it's cached
    if haskey(LaTeXSVG, text)
        svg = LaTeXSVG[text]
    else
        # remove the $
        ts = text.s[2:(end - 1)]
        command = `tex2svg $ts`
        try
            svg = read(command, String)
        catch e
            @warn "Using LaTeX needs the program `tex2svg` which might not be installed"
            @info "It can be installed using `npm install -g mathjax-node-cli`"
            throw(e)
        end
        LaTeXSVG[text] = svg
    end
    pathsvg(svg)
    if action != :path
        # stroke is also fill for letters
        do_action(:fill)
    end
    translate(-pos)
end

"""
    compute_transformation!(action::AbstractAction, video::Video, frame::Int)

Update action.internal_transitions for the current frame number
"""
function compute_transformation!(action::AbstractAction, video::Video, frame::Int)
    for (trans, internal_trans) in zip(action.transitions, action.internal_transitions)
        compute_transition!(internal_trans, trans, video, action, frame)
    end
end

"""
    compute_transition!(internal_rotation::InternalRotation, rotation::Rotation, video,
                        action::AbstractAction, frame)

Computes the rotation transformation for the `action`.
If the `Rotation` is given directly it uses the frame number for interpolation.
If `rotation` includes symbols the current definition of that look up is used for computation.
"""
function compute_transition!(
    internal_rotation::InternalRotation,
    rotation::Rotation,
    video,
    action::AbstractAction,
    frame,
)
    t = (frame - first(get_frames(action))) / (length(get_frames(action)) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
    from, to, center = rotation.from, rotation.to, rotation.center

    center isa Symbol && (center = pos(center))
    from isa Symbol && (from = angle(from))
    to isa Symbol && (to = angle(to))

    internal_rotation.angle = from + t * (to - from)
    internal_rotation.center = center
end

"""
    compute_transition!(internal_translation::InternalTranslation, translation::Translation,
                        video, action::AbstractAction, frame)

Computes the translation transformation for the `action`.
If the `translation` is given directly it uses the frame number for interpolation.
If `translation` includes symbols the current definition of that symbol is looked up
and used for computation.
"""
function compute_transition!(
    internal_translation::InternalTranslation,
    translation::Translation,
    video,
    action::AbstractAction,
    frame,
)
    t = (frame - first(get_frames(action))) / (length(get_frames(action)) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
    from, to = translation.from, translation.to

    from isa Symbol && (from = pos(from))
    to isa Symbol && (to = pos(to))

    internal_translation.by = from + t * (to - from)
end

"""
    compute_transition!(internal_translation::InternalScaling, translation::Scaling,
                        video, action::AbstractAction, frame)

Computes the scaling transformation for the `action`.
If the `scaling` is given directly it uses the frame number for interpolation.
If `scaling` includes symbols, the current definition of that symbol is looked up
and used for computation.
"""
function compute_transition!(
    internal_scale::InternalScaling,
    scale::Scaling,
    video,
    action::AbstractAction,
    frame,
)
    t = (frame - first(get_frames(action))) / (length(get_frames(action)) - 1)
    # makes sense to only allow 0 ≤ t ≤ 1
    t = min(1.0, t)
    from, to = scale.from, scale.to

    if !scale.compute_from_once || frame == first(get_frames(action))
        from isa Symbol && (from = get_scale(from))
        if scale.compute_from_once
            scale.from = from
        end
    end
    to isa Symbol && (to = get_scale(to))
    internal_scale.scale = from .+ t .* (to .- from)
end

"""
    perform_transformation(action::AbstractAction)

Perform the transformations as described in action.internal_transitions
"""
function perform_transformation(action::AbstractAction)
    for trans in action.internal_transitions
        perform_transformation(trans)
    end
end

"""
    perform_transformation(trans::InternalTranslation)

Translate as described in `trans`.
"""
function perform_transformation(trans::InternalTranslation)
    translate(trans.by)
end

"""
    perform_transformation(trans::InternalRotation)

Translate and rotate as described in `trans`.
"""
function perform_transformation(trans::InternalRotation)
    translate(trans.center)
    rotate(trans.angle)
end

"""
    perform_transformation(trans::InternalScaling)

Scale as described in `trans`.
"""
function perform_transformation(trans::InternalScaling)
    scaleto(trans.scale...)
end

"""
    get_value(s::Symbol)

Get access to the value that got saved in `s` by a previous action.
If you want to access a position or angle check out [`get_position`](@ref)
and [`get_angle`](@ref).

# Returns
- `Any`: the value stored by a previous action.
"""
function get_value(s::Symbol)
    is_internal = first(string(s)) == '_'
    if is_internal
        internal_sym = Symbol(string(s)[2:end])
        if hasfield(ActionSetting, internal_sym)
            return getfield(get_current_setting(), internal_sym)
        end
    end

    defs = CURRENT_VIDEO[1].defs
    if haskey(defs, s)
        return defs[s]
    else
        error("The symbol $s is not defined.")
    end
end

"""
    val(x)

`val` is just a short-hand for [`get_value`](@ref)
"""
val(x) = get_value(x)

get_position(p::Point) = p
get_position(t::Transformation) = t.p

"""
    get_position(s::Symbol)

Get access to the position that got saved in `s` by a previous action.

# Returns
- `Point`: the point stored by a previous action.
"""
get_position(s::Symbol) = get_position(val(s))

"""
    pos(x)

`pos` is just a short-hand for [`get_position`](@ref)
"""
pos(x) = get_position(x)

# As it is just the number tuple -> return it
get_scale(x::Tuple{<:Number,<:Number}) = x

# If just the number -> return it as a tuple
get_scale(x::Number) = (x, x)

"""
    get_scale(s::Symbol)

Get access to the scaling that got saved in `s` by a previous action.

# Returns
- `Scaling`: the scale stored by a previous action.
"""
get_scale(s::Symbol) = get_scale(val(s))

"""
    scl(x)

`scl` is just a short-hand for [`get_scale`](@ref)
"""
scl(x) = get_scale(x)

get_angle(t::Transformation) = t.angle

"""
    get_angle(s::Symbol)

Get access to the angle that got saved in `s` by a previous action.

# Returns
- `Float64`: the angle stored by a previous action i.e via `return Transformation(p, angle)`
"""
get_angle(s::Symbol) = get_angle(val(s))

"""
    ang(x)

`ang` is just a short-hand for [`get_angle`](@ref)
"""
ang(x) = get_angle(x)

"""
    projection(p::Point, l::Line)

Return the projection of a point to a line.
"""
function projection(p::Point, l::Line)
    # move line to origin and describe it as a vector
    o = l.p1
    v = l.p2 - o
    # point also moved to origin
    x = p - o

    # scalar product <x,v>/<v,v>
    c = (x.x * v.x + x.y * v.y) / (v.x^2 + v.y^2)
    return c * v + o
end

"""
    create_internal_transitions!(action::AbstractAction)

For every translation an internal translation is added to `action.internal_transitions`.
Same is true for all other transitions.
"""
function create_internal_transitions!(action::AbstractAction)
    for trans in action.transitions
        if trans isa Translation
            push!(action.internal_transitions, InternalTranslation(O))
        elseif trans isa Rotation
            push!(action.internal_transitions, InternalRotation(0.0, O))
        elseif trans isa Scaling
            push!(action.internal_transitions, InternalScaling((1.0, 1.0)))
        end
    end
end

"""
    javis(
        video::Video,
        actions::Vector{AbstractAction};
        framerate=30,
        pathname="",
        tempdirectory=""
    )

Similar to `animate` in Luxor with a slightly different structure.
Instead of using actions and a video instead of scenes in a movie.

# Arguments
- `video::Video`: The video which defines the dimensions of the output
- `actions::Vector{Action}`: All actions that are performed

# Keywords
- `framerate::Int`: The frame rate of the video
- `pathname::String`: The path for the rendered gif or mp4 (i.e `output.gif` or `output.mp4`)
- `tempdirectory::String`: The folder where each frame is stored
    Defaults to a temporary directory when not set

# Example
```
function ground(args...)
    background("white")
    sethue("black")
end

function circ(p=O, color="black")
    sethue(color)
    circle(p, 25, :fill)
    return Transformation(p, 0.0)
end

from = Point(-200, -200)
to = Point(-20, -130)
p1 = Point(0,-100)
p2 = Point(0,-50)
from_rot = 0.0
to_rot = 2π

demo = Video(500, 500)
javis(demo, [
    Action(1:100, ground),
    Action(1:100, :red_ball, (args...)->circ(p1, "red"), Rotation(from_rot, to_rot)),
    Action(1:100, (args...)->circ(p2, "blue"), Rotation(to_rot, from_rot, :red_ball))
], tempdirectory="images", pathname="rotating.gif")
```

This structure makes it possible to refer to positions of previous actions
i.e :red_ball is an id for the position or the red ball which can be used in the
rotation of the next ball.
"""
function javis(
    video::Video,
    actions::Vector{AA};
    framerate = 30,
    pathname = "javis_$(randstring(7)).gif",
    tempdirectory = mktempdir(),
) where {AA<:AbstractAction}
    compute_frames!(actions)

    for action in actions
        compute_frames!(action.subactions; last_frames = get_frames(action))
    end

    # get all frames
    frames = Int[]
    for action in actions
        append!(frames, collect(get_frames(action)))
    end
    frames = unique(frames)

    # create internal transition objects
    for action in actions
        create_internal_transitions!(action)
        for subaction in action.subactions
            create_internal_transitions!(subaction)
        end
    end

    # create defs object
    for action in actions
        if action.id !== nothing
            video.defs[action.id] = Transformation(O, 0.0)
        end
    end

    if isempty(CURRENT_ACTION)
        push!(CURRENT_ACTION, actions[1])
    else
        CURRENT_ACTION[1] = actions[1]
    end

    filecounter = 1
    for frame in frames
        background_settings = ActionSetting()
        Drawing(
            video.width,
            video.height,
            "$(tempdirectory)/$(lpad(filecounter, 10, "0")).png",
        )
        origin()
        origin_matrix = cairotojuliamatrix(getmatrix())
        # this frame needs doing, see if each of the scenes defines it
        for action in actions
            # if action is not in global layer this sets the background_settings
            # from the parent background action
            update_action_settings!(action, background_settings)
            CURRENT_ACTION[1] = action
            if frame in get_frames(action)
                # check if the action should be part of the global layer (i.e BackgroundAction)
                # or in its own layer (default)
                in_global_layer = get(action.opts, :in_global_layer, false)
                if !in_global_layer
                    @layer begin
                        perform_action(action, video, frame, origin_matrix)
                    end
                else
                    perform_action(action, video, frame, origin_matrix)
                    # update origin_matrix as it's inside the global layer
                    origin_matrix = cairotojuliamatrix(getmatrix())
                end
            end
            # if action is in global layer this changes the background settings
            update_background_settings!(background_settings, action)
        end
        finish()
        filecounter += 1
    end

    isempty(pathname) && return
    path, ext = splitext(pathname)
    if ext == ".gif"
        # generate a colorpalette first so ffmpeg does not have to guess it
        ffmpeg_exe(`-loglevel panic -i $(tempdirectory)/%10d.png -vf
                    "palettegen=stats_mode=diff" -y "$(tempdirectory)/palette.bmp"`)
        # then apply the palette to get better results
        ffmpeg_exe(`-loglevel panic -framerate $framerate -i $(tempdirectory)/%10d.png -i
                    "$(tempdirectory)/palette.bmp" -lavfi
                    "paletteuse=dither=sierra2_4a" -y $pathname`)
    elseif ext == ".mp4"
        ffmpeg_exe(`-loglevel panic -framerate $framerate -i $(tempdirectory)/%10d.png  -c:v libx264 -pix_fmt yuv420p $pathname`)

    else
        @error "Currently, only gif and mp4 creation is supported. Not a $ext."
    end
    return pathname
end

function update_background_settings!(setting::ActionSetting, action::Action)
    in_global_layer = get(action.opts, :in_global_layer, false)
    if in_global_layer
        update_ActionSetting!(setting, action.current_setting)
    end
end

function update_action_settings!(action::Action, setting::ActionSetting)
    update_ActionSetting!(action.current_setting, setting)
end

"""
    perform_action(action, video, frame, origin_matrix)

Is called inside the `javis` and does everything handled for an `AbstractAction`.
It is a 4-step process:
- compute the transformation for this action (translation, rotation, scale)
- do the transformation
- call the action function
- save the result of the action if wanted inside `video.defs`
"""
function perform_action(action, video, frame, origin_matrix)
    # first compute and perform the global transformations of this action
    compute_transformation!(action, video, frame)
    perform_transformation(action)

    # relative frame number for subactions
    rel_frame = frame - first(get_frames(action)) + 1
    # call currently active subactions and their transformations
    for subaction in action.subactions
        if rel_frame in get_frames(subaction)
            subaction.func(video, action, subaction, rel_frame)
            compute_transformation!(subaction, video, rel_frame)
            perform_transformation(subaction)
        elseif rel_frame > last(get_frames(subaction))
            # call the subaction on the last frame i.e. disappeared things stay disappeared
            subaction.func(video, action, subaction, last(get_frames(subaction)))
            # have the transformation from the last active frame
            compute_transformation!(subaction, video, last(get_frames(subaction)))
            perform_transformation(subaction)
        end
    end

    # set the defaults for the frame like setline() and setopacity()
    # which can depend on the subactions
    set_action_defaults!(action)
    res = action.func(video, action, frame)
    if action.id !== nothing
        current_global_matrix = cairotojuliamatrix(getmatrix())
        # obtain current matrix without the initial matrix part
        current_matrix = inv(origin_matrix) * current_global_matrix

        # if a transformation let's save the global coordinates
        if res isa Point
            vec = current_matrix * [res.x, res.y, 1.0]
            video.defs[action.id] = Point(vec[1], vec[2])
        elseif res isa Transformation
            trans = current_matrix * res
            video.defs[action.id] = trans
        else # just save the result such that it can be used as one wishes
            video.defs[action.id] = res
        end
    end
end

"""
    set_action_defaults!(action)

Set the default action values
- line_width and calls `Luxor.setline`.
- opacity and calls `Luxor.opacity`.
- scale and calls `Luxor.scale`.
"""
function set_action_defaults!(action)
    cs = action.current_setting
    current_line_width = cs.line_width * cs.mul_line_width
    Luxor.setline(current_line_width)
    current_opacity = cs.opacity * cs.mul_opacity
    Luxor.setopacity(current_opacity)

    desired_scale = cs.desired_scale .* cs.mul_scale
    scaleto(desired_scale...)
end

const LUXOR_DONT_EXPORT = [
    :boundingbox,
    :Boxmaptile,
    :Sequence,
    :setline,
    :setopacity,
    :fontsize,
    :get_fontsize,
    :scale,
]

# Export each function from Luxor
for func in names(Luxor; imported = true)
    if !(func in LUXOR_DONT_EXPORT)
        eval(Meta.parse("import Luxor." * string(func)))
        eval(Expr(:export, func))
    end
end

export javis, latex
export Video, Action, BackgroundAction, SubAction, Rel
export Line, Translation, Rotation, Transformation, Scaling
export val, pos, ang, get_value, get_position, get_angle
export projection, morph
export appear, disappear

# custom override of luxor extensions
export setline, setopacity, fontsize, get_fontsize, scale

end


text:
#
# Slater exchange with alpha=2/3, spin-polarized case
#
function XC_x_slater_spin( rho, zeta )

    f = -1.10783814957303361 # f = -9/8*(3/pi)^(1/3)
    alpha = 2.0/3.0
    third = 1.0/3.0
    p43 = 4.0/3.0

    rho13 = ( (1.0 + zeta) * rho)^third
    exup = f * alpha * rho13
    vxup = p43 * f * alpha * rho13

    rho13 = ( (1.0 - zeta) * rho)^third
    exdw = f * alpha * rho13
    vxdw = p43 * f * alpha * rho13
    ex = 0.5 * ( (1.0 + zeta) * exup + (1.0 - zeta) * exdw)

    return ex, vxup, vxdw
end

ext: jl


text:
# double a number, the macro way!
#=
macro doubled(ex)
    return :( 2 * $(esc(ex)))
end
=#

# This code does not work.  Don't try it.
#=
macro doubled(ex)
    if typeof(ex) isa AbstractFloat
        return :( double_super_duper($(esc(ex))) )
    else
        return :( 2 * $(esc(ex)))
    end
end
=#

# generated function
@generated function doubled(x)
    return :( 2 * x )
end
doubled(2)
#=
julia> @generated function doubled(x)
           return :( 2 * x )
       end

julia> doubled(2)
4
=#

# -----------------------------------------------------------------------------
# generated function arguments

# debug
@generated function doubled(x)
    @show x
    return :( 2 * x )
end
doubled(2)
doubled(2)
#=
julia> @generated function doubled(x)
           @show x
           return :( 2 * x )
       end

julia> doubled(2)
x = Int64
4
=#

doubled(3.0)
#=
julia> doubled(3.0)
x = Float64
6.0
=#

# using type info
#=
julia> @generated function doubled(x)
           if x <: AbstractFloat
               return :( double_super_duper(x) )
           else
               return :( 2 * x )
           end
       end
=#


text:
const test_bounds = [
    Bounds.Ellipsoid,
    Bounds.MultiEllipsoid
]
const test_props = [
    Proposals.Rejection(maxiter=Int(1e6)),
    Proposals.RWalk(ratio=0.5, walks=50),
    Proposals.RStagger(ratio=0.5, walks=50),
    Proposals.Slice(slices=10),
    Proposals.RSlice(slices=10)
]

const MAXZSCORES = Dict(zip(
    Iterators.product(test_bounds, test_props),
    [4, 3, 9, 8, 6, 3, 5, 7, 4, 4] # rwalk is bad...
))

function test_logz(measured, actual, error, bound, proposal)
    diff = measured - actual
    zscore = abs(diff) / error
    @test measured ≈ actual atol=MAXZSCORES[(bound, proposal)] * error
end

@testset "$(nameof(bound)), $(nameof(typeof(proposal)))" for bound in test_bounds, proposal in test_props
    @testset "Correlated Gaussian Conjugate Prior - ndims=$D" for D in [2, 4]
        model, logz = Models.CorrelatedGaussian(D)
        # match JAXNS paper setup, generally
        sampler = Nested(D, 50D; bounds=bound, proposal=proposal)

        chain, state = sample(rng, model, sampler; dlogz=0.01)
        chain_res = sample(rng, chain, Weights(vec(chain[:weights])), length(chain))
        # test posteriors
        vals = Array(chain_res)
        means = mean(vals, dims=1)
        tols = 2std(vals, mean=means, dims=1) # 2-sigma
        μ = fill(2.0, D)
        Σ = fill(0.95, D, D)
        Σ[diagind(Σ)] .= 1
        expected = Σ * ((Σ + I) \ μ)
        @test all(@.(abs(means - expected) < tols))

        # logz
        test_logz(state.logz, logz, state.logzerr, bound, proposal)
    end

    @testset "Gaussian Shells" begin
        model, logz = Models.GaussianShells()

        sampler = Nested(2, 1000; bounds=bound, proposal=proposal)

        chain, state = sample(rng, model, sampler; dlogz=0.01)

        # logz
        test_logz(state.logz, logz, state.logzerr, bound, proposal)
    end

    @testset "Gaussian Mixture Model" begin
        σ = 0.1
        μ1 = ones(2)
        μ2 = -ones(2)
        inv_σ = diagm(0 => fill(1 / σ^2, 2))

        function logl(x)
            dx1 = x .- μ1
            dx2 = x .- μ2
            f1 = -dx1' * (inv_σ * dx1) / 2
            f2 = -dx2' * (inv_σ * dx2) / 2
            return logaddexp(f1, f2)
        end

        prior(X) = muladd.(10, X, -5)
        model = NestedModel(logl, prior)

        analytic_logz = log(4π * σ^2 / 100)

        spl = Nested(2, 1000, bounds=bound, proposal=proposal)
        chain, state = sample(rng, model, spl; dlogz=0.01)
        chain_res = sample(rng, chain, Weights(vec(chain[:weights])), length(chain))

        test_logz(state.logz, analytic_logz, state.logzerr, bound, proposal)

        xmodes = sort!(findpeaks(chain_res[:, 1, 1])[1:2])
        @test xmodes[1] ≈ -1 atol = σ
        @test xmodes[2] ≈ 1 atol = σ
        ymodes = sort!(findpeaks(chain_res[:, 2, 1])[1:2])
        @test ymodes[1] ≈ -1 atol = σ
        @test ymodes[2] ≈ 1 atol = σ
    end

    @testset "Eggbox" begin
        model, logz = Models.Eggbox()

        sampler = Nested(2, 1000; bounds=bound, proposal=proposal)

        chain, state = sample(rng, model, sampler; dlogz=0.1)

        test_logz(state.logz, logz, state.logzerr, bound, proposal)

        chain_res = sample(rng, chain, Weights(vec(chain[:weights])), length(chain))
        xmodes = sort!(findpeaks(chain_res[:, 1, 1])[1:5])
        @test all(isapprox.(xmodes, 0.1:0.2:0.9, atol=0.2))
        ymodes = sort!(findpeaks(chain_res[:, 2, 1])[1:5])
        @test all(isapprox.(ymodes, 0.1:0.2:0.9, atol=0.2))
    end
end


text: # TODO: write device test

max_stars_repo_name: unbornchikken/julia-ml-proto

max_issues_repo_name: unbornchikken/julia-ml-proto

max_forks_repo_name: unbornchikken/julia-ml-proto


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module QueryTests

using Test
using Dates

Jan = Dates.DateTime(2013, 1, 1) # Tuesday
Feb = Dates.DateTime(2013, 2, 2) # Saturday
Mar = Dates.DateTime(2013, 3, 3) # Sunday
Apr = Dates.DateTime(2013, 4, 4) # Thursday
May = Dates.DateTime(2013, 5, 5) # Sunday
Jun = Dates.DateTime(2013, 6, 7) # Friday
Jul = Dates.DateTime(2013, 7, 7) # Sunday
Aug = Dates.DateTime(2013, 8, 8) # Thursday
Sep = Dates.DateTime(2013, 9, 9) # Monday
Oct = Dates.DateTime(2013, 10, 10) # Thursday
Nov = Dates.DateTime(2013, 11, 11) # Monday
Dec = Dates.DateTime(2013, 12, 11) # Wednesday
monthnames = ["January", "February", "March", "April",
                "May", "June", "July", "August", "September",
                "October", "November", "December"]
daysofweek = [Dates.Tue, Dates.Sat, Dates.Sun, Dates.Thu, Dates.Sun, Dates.Fri,
              Dates.Sun, Dates.Thu, Dates.Mon, Dates.Thu, Dates.Mon, Dates.Wed]
dows = ["Tuesday", "Saturday", "Sunday", "Thursday", "Sunday", "Friday",
        "Sunday", "Thursday", "Monday", "Thursday", "Monday", "Wednesday"]
daysinmonth = [31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31]
@testset "Name functions" begin
    for (i, dt) in enumerate([Jan, Feb, Mar, Apr, May, Jun, Jul, Aug, Sep, Oct, Nov, Dec])
        @test Dates.month(dt) == i
        @test Dates.monthname(dt) == monthnames[i]
        @test Dates.monthname(i) == monthnames[i]
        @test Dates.monthabbr(dt) == monthnames[i][1:3]
        @test Dates.monthabbr(i) == monthnames[i][1:3]
        @test Dates.dayofweek(dt) == daysofweek[i]
        @test Dates.dayname(dt) == dows[i]
        @test Dates.dayname(Dates.dayofweek(dt)) == dows[i]
        @test Dates.dayabbr(dt) == dows[i][1:3]
        @test Dates.dayabbr(Dates.dayofweek(dt)) == dows[i][1:3]
        @test Dates.daysinmonth(dt) == daysinmonth[i]
    end
end
@testset "Customizing locale" begin
    Dates.LOCALES["french"] = Dates.DateLocale(
        ["janvier", "février", "mars", "avril", "mai", "juin",
         "juillet", "août", "septembre", "octobre", "novembre", "décembre"],
        ["janv", "févr", "mars", "avril", "mai", "juin",
         "juil", "août", "sept", "oct", "nov", "déc"],
        ["lundi", "mardi", "mercredi", "jeudi", "vendredi", "samedi", "dimanche"],
        [""],
    )
    @test Dates.dayname(Nov; locale="french") == "lundi"
    @test Dates.dayname(Jan; locale="french") == "mardi"
    @test Dates.dayname(Dec; locale="french") == "mercredi"
    @test Dates.dayname(Apr; locale="french") == "jeudi"
    @test Dates.dayname(Jun; locale="french") == "vendredi"
    @test Dates.dayname(Feb; locale="french") == "samedi"
    @test Dates.dayname(May; locale="french") == "dimanche"

    @test Dates.monthname(Jan; locale="french") == "janvier"
    @test Dates.monthname(Feb; locale="french") == "février"
    @test Dates.monthname(Mar; locale="french") == "mars"
    @test Dates.monthname(Apr; locale="french") == "avril"
    @test Dates.monthname(May; locale="french") == "mai"
    @test Dates.monthname(Jun; locale="french") == "juin"
    @test Dates.monthname(Jul; locale="french") == "juillet"
    @test Dates.monthname(Aug; locale="french") == "août"
    @test Dates.monthname(Sep; locale="french") == "septembre"
    @test Dates.monthname(Oct; locale="french") == "octobre"
    @test Dates.monthname(Nov; locale="french") == "novembre"
    @test Dates.monthname(Dec; locale="french") == "décembre"
end
@testset "leap years" begin
    @test Dates.isleapyear(Dates.DateTime(1900)) == false
    @test Dates.isleapyear(Dates.DateTime(2000)) == true
    @test Dates.isleapyear(Dates.DateTime(2004)) == true
    @test Dates.isleapyear(Dates.DateTime(2008)) == true
    @test Dates.isleapyear(Dates.DateTime(0)) == true
    @test Dates.isleapyear(Dates.DateTime(1)) == false
    @test Dates.isleapyear(Dates.DateTime(-1)) == false
    @test Dates.isleapyear(Dates.DateTime(4)) == true
    @test Dates.isleapyear(Dates.DateTime(-4)) == true

    @test Dates.daysinyear(2000) == 366
    @test Dates.daysinyear(2001) == 365
    @test Dates.daysinyear(2000) == 366
    @test Dates.daysinyear(2001) == 365

    @test Dates.daysinyear(Dates.Date(2000)) == 366
    @test Dates.daysinyear(Dates.Date(2001)) == 365
    @test Dates.daysinyear(Dates.DateTime(2000)) == 366
    @test Dates.daysinyear(Dates.DateTime(2001)) == 365
end
@testset "dayofweek" begin
    # Days of week from Monday = 1 to Sunday = 7
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 22)) == 7
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 23)) == 1
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 24)) == 2
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 25)) == 3
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 26)) == 4
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 27)) == 5
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 28)) == 6
    @test Dates.dayofweek(Dates.DateTime(2013, 12, 29)) == 7
end
@testset "daysofweekinmonth" begin
    # There are 5 Sundays in December, 2013
    @test Dates.daysofweekinmonth(Dates.DateTime(2013, 12, 1)) == 5
    # There are 4 Sundays in November, 2013
    @test Dates.daysofweekinmonth(Dates.DateTime(2013, 11, 24)) == 4

    @test Dates.dayofweekofmonth(Dates.DateTime(2013, 12, 1)) == 1
    @test Dates.dayofweekofmonth(Dates.DateTime(2013, 12, 8)) == 2
    @test Dates.dayofweekofmonth(Dates.DateTime(2013, 12, 15)) == 3
    @test Dates.dayofweekofmonth(Dates.DateTime(2013, 12, 22)) == 4
    @test Dates.dayofweekofmonth(Dates.DateTime(2013, 12, 29)) == 5
end
@testset "dayofyear" begin
    @test Dates.dayofyear(2000, 1, 1) == 1
    @test Dates.dayofyear(2004, 1, 1) == 1
    @test Dates.dayofyear(20013, 1, 1) == 1
    # Leap year
    @test Dates.dayofyear(2000, 12, 31) == 366
    # Non-leap year
    @test Dates.dayofyear(2001, 12, 31) == 365

    @test Dates.dayofyear(Dates.DateTime(2000, 1, 1)) == 1
    @test Dates.dayofyear(Dates.DateTime(2004, 1, 1)) == 1
    @test Dates.dayofyear(Dates.DateTime(20013, 1, 1)) == 1
    # Leap year
    @test Dates.dayofyear(Dates.DateTime(2000, 12, 31)) == 366
    # Non-leap year
    @test Dates.dayofyear(Dates.DateTime(2001, 12, 31)) == 365
    # Test every day of a year
    dt = Dates.DateTime(2000, 1, 1)
    for i = 1:366
        @test Dates.dayofyear(dt) == i
        dt += Dates.Day(1)
    end
    dt = Dates.DateTime(2001, 1, 1)
    for i = 1:365
        @test Dates.dayofyear(dt) == i
        dt += Dates.Day(1)
    end
end
@testset "quarterofyear" begin
    @test Dates.quarterofyear(Dates.Date(2000, 1, 1))  == 1
    @test Dates.quarterofyear(Dates.Date(2000, 1, 31))  == 1
    @test Dates.quarterofyear(Dates.Date(2000, 2, 1))  == 1
    @test Dates.quarterofyear(Dates.Date(2000, 2, 29))  == 1
    @test Dates.quarterofyear(Dates.Date(2000, 3, 1))  == 1
    @test Dates.quarterofyear(Dates.Date(2000, 3, 31))  == 1
    @test Dates.quarterofyear(Dates.Date(2000, 4, 1)) == 2
    @test Dates.quarterofyear(Dates.Date(2000, 4, 30)) == 2
    @test Dates.quarterofyear(Dates.Date(2000, 5, 1)) == 2
    @test Dates.quarterofyear(Dates.Date(2000, 5, 31)) == 2
    @test Dates.quarterofyear(Dates.Date(2000, 6, 1)) == 2
    @test Dates.quarterofyear(Dates.Date(2000, 6, 30)) == 2
    @test Dates.quarterofyear(Dates.Date(2000, 7, 1)) == 3
    @test Dates.quarterofyear(Dates.Date(2000, 7, 31)) == 3
    @test Dates.quarterofyear(Dates.Date(2000, 8, 1)) == 3
    @test Dates.quarterofyear(Dates.Date(2000, 8, 31)) == 3
    @test Dates.quarterofyear(Dates.Date(2000, 9, 1)) == 3
    @test Dates.quarterofyear(Dates.Date(2000, 9, 30)) == 3
    @test Dates.quarterofyear(Dates.Date(2000, 10, 1)) == 4
    @test Dates.quarterofyear(Dates.Date(2000, 10, 31)) == 4
    @test Dates.quarterofyear(Dates.Date(2000, 11, 1)) == 4
    @test Dates.quarterofyear(Dates.Date(2000, 11, 30)) == 4
    @test Dates.quarterofyear(Dates.Date(2000, 12, 1)) == 4
    @test Dates.quarterofyear(Dates.Date(2000, 12, 31)) == 4

    @test Dates.quarterofyear(Dates.DateTime(2000, 1, 1))  == 1
    @test Dates.quarterofyear(Dates.DateTime(2000, 1, 31))  == 1
    @test Dates.quarterofyear(Dates.DateTime(2000, 2, 1))  == 1
    @test Dates.quarterofyear(Dates.DateTime(2000, 2, 29))  == 1
    @test Dates.quarterofyear(Dates.DateTime(2000, 3, 1))  == 1
    @test Dates.quarterofyear(Dates.DateTime(2000, 3, 31))  == 1
    @test Dates.quarterofyear(Dates.DateTime(2000, 4, 1)) == 2
    @test Dates.quarterofyear(Dates.DateTime(2000, 4, 30)) == 2
    @test Dates.quarterofyear(Dates.DateTime(2000, 5, 1)) == 2
    @test Dates.quarterofyear(Dates.DateTime(2000, 5, 31)) == 2
    @test Dates.quarterofyear(Dates.DateTime(2000, 6, 1)) == 2
    @test Dates.quarterofyear(Dates.DateTime(2000, 6, 30)) == 2
    @test Dates.quarterofyear(Dates.DateTime(2000, 7, 1)) == 3
    @test Dates.quarterofyear(Dates.DateTime(2000, 7, 31)) == 3
    @test Dates.quarterofyear(Dates.DateTime(2000, 8, 1)) == 3
    @test Dates.quarterofyear(Dates.DateTime(2000, 8, 31)) == 3
    @test Dates.quarterofyear(Dates.DateTime(2000, 9, 1)) == 3
    @test Dates.quarterofyear(Dates.DateTime(2000, 9, 30)) == 3
    @test Dates.quarterofyear(Dates.DateTime(2000, 10, 1)) == 4
    @test Dates.quarterofyear(Dates.DateTime(2000, 10, 31)) == 4
    @test Dates.quarterofyear(Dates.DateTime(2000, 11, 1)) == 4
    @test Dates.quarterofyear(Dates.DateTime(2000, 11, 30)) == 4
    @test Dates.quarterofyear(Dates.DateTime(2000, 12, 1)) == 4
    @test Dates.quarterofyear(Dates.DateTime(2000, 12, 31)) == 4
end
@testset "dayofquarter" begin
    # Non-leap and leap year
    for y in (2014, 2020)
        @test Dates.dayofquarter(Dates.Date(y, 1, 1)) == 1
        @test Dates.dayofquarter(Dates.Date(y, 4, 1)) == 1
        @test Dates.dayofquarter(Dates.Date(y, 7, 1)) == 1
        @test Dates.dayofquarter(Dates.Date(y, 10, 1)) == 1
        @test Dates.dayofquarter(Dates.Date(y, 3, 31)) == 90 + isleapyear(y)
        @test Dates.dayofquarter(Dates.Date(y, 6, 30)) == 91
        @test Dates.dayofquarter(Dates.Date(y, 9, 30)) == 92
        @test Dates.dayofquarter(Dates.Date(y, 12, 31)) == 92
        @test Dates.dayofquarter(Dates.DateTime(y, 1, 1)) == 1
        @test Dates.dayofquarter(Dates.DateTime(y, 4, 1)) == 1
        @test Dates.dayofquarter(Dates.DateTime(y, 7, 1)) == 1
        @test Dates.dayofquarter(Dates.DateTime(y, 10, 1)) == 1
        @test Dates.dayofquarter(Dates.DateTime(y, 3, 31)) == 90 + isleapyear(y)
        @test Dates.dayofquarter(Dates.DateTime(y, 6, 30)) == 91
        @test Dates.dayofquarter(Dates.DateTime(y, 9, 30)) == 92
        @test Dates.dayofquarter(Dates.DateTime(y, 12, 31)) == 92
    end
end

end

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
using Restarts
using StatsBase
using Test

if !@isdefined(Returns)  # Julia < 1.7
    Returns(x) = (_args...; _kwargs...) -> x
end

@define_restart use_value
@define_restart reparse_entry

struct MalformedLogEntryError <: Exception
    text::String
end

function parse_log_entry(text)
    parts = split(text, ":", limit = 2)
    if length(parts) == 2
        return (level = Symbol(strip(parts[1])), text = strip(parts[2]))
    end

    restart(
        MalformedLogEntryError(text),
        use_value => value -> value,  # or `use_value => identity`
        reparse_entry => fixed_text -> parse_log_entry(fixed_text),
    )
end

function test_parse_log_entry()
    @test parse_log_entry("INFO: message") == (level = :INFO, text = "message")
    @test_throws MalformedLogEntryError parse_log_entry("malformed log entry")
    @test Restarts.with(reparse_entry => ex -> Some("_MALFORMED_:" * ex.text)) do
        parse_log_entry("malformed log entry")
    end == (level = :_MALFORMED_, text = "malformed log entry")
end

function countlevels(io::IO)
    Restarts.with(use_value => Returns(Some((level = :_MALFORMED_, text = "")))) do
        countmap(map(t -> parse_log_entry(t).level, eachline(io)))
    end
end

function test_count_levels()
    io = IOBuffer("""
    INFO: message 1
    INFO: message 2
    malformed
    DEBUG: message 2
    """)
    levels = countlevels(io)
    @test levels == Dict(:INFO => 2, :DEBUG => 1, :_MALFORMED_ => 1)
end


text:
using Test, JuLIP

at = bulk(:Si)

at2 = copy(at)
@test at == at2
@test at !== at2
@test at.X === at2.X

at3 = deepcopy(at)
@test at == at3
@test at !== at3
@test at.X !== at3.X
@test at.X == at3.X

rot_at = rotate!(deepcopy(at), x=[1,1,1], y=[1,-1,0])
apply_defm!(rot_at, inv(rotation_matrix(x=[1,1,1], y=[1,-1,0])))
@test rot_at ≈ at


text:
_rng() = MersenneTwister(123456)

function generate_noise_matrix(rng::AbstractRNG, N::Int)
    A = randn(rng, N, N)
    return Symmetric(A * A' + I)
end

@testset "finite_gp" begin
    @testset "convert" begin
        f = GP(sin, SqExponentialKernel())
        x = randn(10)
        Σy = generate_noise_matrix(Random.GLOBAL_RNG, 10)
        fx = f(x, Σy)
        @test fx isa AbstractGPs.FiniteGP

        dist = @inferred(convert(MvNormal, fx))
        @test dist isa MvNormal{Float64}
        @test mean(dist) ≈ mean(fx)
        @test cov(dist) ≈ cov(fx)

        dist = @inferred(convert(MvNormal{Float32}, fx))
        @test dist isa MvNormal{Float32}
        @test mean(dist) ≈ mean(fx)
        @test cov(dist) ≈ cov(fx)
    end
    @testset "statistics" begin
        rng, N, N′ = MersenneTwister(123456), 1, 9
        x, x′, Σy, Σy′ = randn(rng, N), randn(rng, N′), zeros(N, N), zeros(N′, N′)
        σ² = 1e-3
        Xmat = randn(rng, N, N′)
        f = GP(sin, SqExponentialKernel())
        fx, fx′ = f(x, Σy), f(x′, Σy′)
        @test fx isa AbstractGPs.FiniteGP
        @test fx′ isa AbstractGPs.FiniteGP

        for (x, obsdim) in ((RowVecs(Xmat), 1), (ColVecs(Xmat), 2))
            @test f(x) isa AbstractGPs.FiniteGP
            @test f(x, σ²) isa AbstractGPs.FiniteGP
            @test f(Xmat; obsdim=obsdim) == f(x)
            @test f(Xmat, σ²; obsdim=obsdim) == f(x, σ²)
        end
        @test @test_deprecated(f(Xmat)) == f(ColVecs(Xmat))

        @test mean(fx) == mean(f, x)
        @test cov(fx) == cov(f, x)
        @test var(fx) == diag(cov(fx))
        @test cov(fx, fx′) == cov(f, x, x′)
        @test mean.(marginals(fx)) == mean(f(x))
        @test var.(marginals(fx)) == var(f, x)
        @test std.(marginals(fx)) == sqrt.(var(f, x))
        let
            m, C = mean_and_cov(fx)
            @test m == mean(fx)
            @test C == cov(fx)
        end
        let
            m, c = mean_and_var(fx)
            @test m == mean(fx)
            @test c == var(fx)
        end
    end
    @testset "rand (deterministic)" begin
        rng = MersenneTwister(123456)
        N = 10
        x = randn(rng, N)
        Σy = generate_noise_matrix(rng, N)
        fx = GP(1, SqExponentialKernel())(x, Σy)
        @test fx isa AbstractGPs.FiniteGP

        # Check that samples are the correct size.
        @test length(rand(rng, fx)) == length(x)
        @test size(rand(rng, fx, 10)) == (length(x), 10)
        @test length(rand(fx)) == length(x)
        @test size(rand(fx, 10)) == (length(x), 10)

        # Check that `rand!` calls do not error
        y = similar(x)
        rand!(rng, fx, y)
        rand!(fx, y)
        ys = similar(x, length(x), 10)
        rand!(rng, fx, ys)
        rand!(fx, ys)
    end
    @testset "rand (statistical)" begin
        rng = MersenneTwister(123456)
        N = 10
        m0 = 1
        S = 100_000
        x = range(-3.0, 3.0; length=N)
        f = GP(m0, SqExponentialKernel())(x, 1e-12)
        @test f isa AbstractGPs.FiniteGP

        # Check mean + covariance estimates approximately converge for single-GP sampling.
        f̂1 = rand(rng, f, S)
        f̂2 = similar(f̂1)
        rand!(rng, f, f̂2)

        for f̂ in (f̂1, f̂2)
            @test maximum(abs.(mean(f̂; dims=2) - mean(f))) < 1e-2

            Σ′ = (f̂ .- mean(f)) * (f̂ .- mean(f))' ./ S
            @test mean(abs.(Σ′ - cov(f))) < 1e-2
        end
    end
    # @testset "rand (gradients)" begin
    #     rng, N, S = MersenneTwister(123456), 10, 3
    #     x = collect(range(-3.0, stop=3.0, length=N))
    #     Σy = 1e-12

    #     # Check that the gradient w.r.t. the samples is correct (single-sample).
    #     adjoint_test(
    #         x->rand(MersenneTwister(123456), GP(sin, EQ(), GPC())(x, Σy)),
    #         randn(rng, N),
    #         x;
    #         atol=1e-9, rtol=1e-9,
    #     )

    #     # Check that the gradient w.r.t. the samples is correct (multisample).
    #     adjoint_test(
    #         x->rand(MersenneTwister(123456), GP(sin, EQ(), GPC())(x, Σy), S),
    #         randn(rng, N, S),
    #         x;
    #         atol=1e-9, rtol=1e-9,
    #     )
    # end
    # @testset "tr_Cf_invΣy" begin
    #     N = 11
    #     x = collect(range(-3.0, 3.0; length=N))
    #     @testset "dense" begin
    #         rng = MersenneTwister(123456)
    #         A = randn(rng, N, N - 2)
    #         adjoint_test(
    #             (x, A)->begin
    #                 f = GP(sin, EQ(), GPC())
    #                 Σy = _to_psd(A)
    #                 C = cholesky(Σy)
    #                 return tr_Cf_invΣy(f(x, Σy), Σy, C)
    #             end,
    #             randn(rng), x, A,
    #         )
    #     end
    #     @testset "Diagonal" begin
    #         rng = MersenneTwister(123456)
    #         a = 0.01 .* randn(rng, N)
    #         adjoint_test(
    #             (x, a)->begin
    #                 f = GP(sin, EQ(), GPC())
    #                 Σy = Diagonal(exp.(a .+ 1))
    #                 C = cholesky(Σy)
    #                 return tr_Cf_invΣy(f(x, Σy), Σy, C)
    #             end,
    #             randn(rng), x, a,
    #         )
    #     end
    # end
    @testset "logpdf / loglikelihood" begin
        rng = MersenneTwister(123456)
        N = 10
        S = 11
        σ = 1e-1
        x = collect(range(-3.0; stop=3.0, length=N))
        f = GP(1, SqExponentialKernel())
        fx = f(x, 0)
        @test fx isa AbstractGPs.FiniteGP
        y = f(x, σ^2)
        @test y isa AbstractGPs.FiniteGP
        ŷ = rand(rng, y)

        # Check that logpdf returns the correct type and roughly agrees with Distributions.
        @test logpdf(y, ŷ) isa Real
        @test logpdf(y, ŷ) ≈ logpdf(MvNormal(Vector(mean(y)), cov(y)), ŷ)
        @test loglikelihood(y, ŷ) == logpdf(y, ŷ)
        # Check that multi-sample logpdf returns the correct type and is consistent with
        # single-sample logpdf
        Ŷ = rand(rng, y, S)
        @test logpdf(y, Ŷ) isa Vector{Float64}
        @test logpdf(y, Ŷ) ≈ [logpdf(y, Ŷ[:, n]) for n in 1:S]
        @test loglikelihood(y, Ŷ) == sum(logpdf(y, Ŷ))

        # # Check gradient of logpdf at mean is zero for `f`.
        # adjoint_test(ŷ->logpdf(fx, ŷ), 1, ones(size(ŷ)))
        # lp, back = Zygote.pullback(ŷ->logpdf(fx, ŷ), ones(size(ŷ)))
        # @test back(randn(rng))[1] == zeros(size(ŷ))

        # # Check that gradient of logpdf at mean is zero for `y`.
        # adjoint_test(ŷ->logpdf(y, ŷ), 1, ones(size(ŷ)))
        # lp, back = Zygote.pullback(ŷ->logpdf(y, ŷ), ones(size(ŷ)))
        # @test back(randn(rng))[1] == zeros(size(ŷ))

        # # Check that gradient w.r.t. inputs is approximately correct for `f`.
        # x, l̄ = randn(rng, N), randn(rng)
        # adjoint_test(
        #     x->logpdf(f(x, 1e-3), ones(size(x))),
        #     l̄, collect(x);
        #     atol=1e-8, rtol=1e-8,
        # )
        # adjoint_test(
        #     x->sum(logpdf(f(x, 1e-3), ones(size(Ŷ)))),
        #     l̄, collect(x);
        #     atol=1e-8, rtol=1e-8,
        # )

        # # Check that the gradient w.r.t. the noise is approximately correct for `f`.
        # σ_ = randn(rng)
        # adjoint_test((σ_, ŷ)->logpdf(f(x, exp(σ_)), ŷ), l̄, σ_, ŷ)
        # adjoint_test((σ_, Ŷ)->sum(logpdf(f(x, exp(σ_)), Ŷ)), l̄, σ_, Ŷ)

        # # Check that the gradient w.r.t. a scaling of the GP works.
        # adjoint_test(
        #     α->logpdf((α * f)(x, 1e-1), ŷ), l̄, randn(rng);
        #     atol=1e-8, rtol=1e-8,
        # )
        # adjoint_test(
        #     α->sum(logpdf((α * f)(x, 1e-1), Ŷ)), l̄, randn(rng);
        #     atol=1e-8, rtol=1e-8,
        # )
    end
    @testset "Type Stability - $T" for T in [Float64, Float32]
        rng = MersenneTwister(123456)
        x = randn(rng, T, 123)
        z = randn(rng, T, 13)
        f = GP(T(0), SqExponentialKernel())

        fx = f(x, T(0.1))

        y = rand(rng, fx)
        @test y isa Vector{T}
        @test logpdf(fx, y) isa T
    end
    @testset "AbstractMvNormal API" begin
        rng = MersenneTwister(424242)
        N = 5
        x = randn(rng, N)
        f = GP(SqExponentialKernel())
        fx = f(x, 0.1)
        y = rand(rng, N)
        Y = rand(rng, N, 10)
        r = rand(rng, 10)

        Distributions.TestUtils.test_mvnormal(fx, 10^6, rng)
        @test Distributions.invcov(fx) ≈ inv(cov(fx))
        @test Distributions.gradlogpdf(fx, y) ≈
            first(FiniteDifferences.grad(central_fdm(3, 1), Base.Fix1(logpdf, fx), y))
        @test Distributions.sqmahal!(r, fx, Y) ≈ Distributions.sqmahal(fx, Y)
    end
end

@testset "Docs" begin
    docstring = string(Docs.doc(logpdf, Tuple{AbstractGPs.FiniteGP,Vector{Float64}}))
    @test occursin("logpdf(f::FiniteGP, y::AbstractVecOrMat{<:Real})", docstring)
end

# """
#     simple_gp_tests(rng::AbstractRNG, f::GP, xs::AV{<:AV}, σs::AV{<:Real})

# Integration tests for simple GPs.
# """
# function simple_gp_tests(
#     rng::AbstractRNG,
#     f::GP,
#     xs::AV{<:AV},
#     isp_σs::AV{<:Real};
#     atol=1e-8,
#     rtol=1e-8,
# )
#     for x in xs, isp_σ in isp_σs

#         # Test gradient w.r.t. random sampling.
#         N = length(x)
#         adjoint_test(
#             (x, isp_σ)->rand(_rng(), f(x, exp(isp_σ)^2)),
#             randn(rng, N),
#             x,
#             isp_σ,;
#             atol=atol, rtol=rtol,
#         )
#         adjoint_test(
#             (x, isp_σ)->rand(_rng(), f(x, exp(isp_σ)^2), 11),
#             randn(rng, N, 11),
#             x,
#             isp_σ,;
#             atol=atol, rtol=rtol,
#         )

#         # Check that gradient w.r.t. logpdf is correct.
#         y, l̄ = rand(rng, f(x, exp(isp_σ))), randn(rng)
#         adjoint_test(
#             (x, isp_σ, y)->logpdf(f(x, exp(isp_σ)), y),
#             l̄, x, isp_σ, y;
#             atol=atol, rtol=rtol,
#         )

#         # Check that elbo is tight-ish when it's meant to be.
#         fx, yx = f(x, 1e-9), f(x, exp(isp_σ))
#         @test isapprox(elbo(yx, y, fx), logpdf(yx, y); atol=1e-6, rtol=1e-6)

#         # Check that gradient w.r.t. elbo is correct.
#         adjoint_test(
#             (x, ŷ, isp_σ)->elbo(f(x, exp(isp_σ)), ŷ, f(x, 1e-9)),
#             randn(rng), x, y, isp_σ;
#             atol=1e-6, rtol=1e-6,
#         )
#     end
# end

# __foo(x) = isnothing(x) ? "nothing" : x

# @testset "FiniteGP (integration)" begin
#     rng = MersenneTwister(123456)
#     xs = [collect(range(-3.0, stop=3.0, length=N)) for N in [2, 5, 10]]
#     σs = log.([1e-1, 1e0, 1e1])
#     for (k, name, atol, rtol) in vcat(
#         [
#             (EQ(), "EQ", 1e-6, 1e-6),
#             (Linear(), "Linear", 1e-6, 1e-6),
#             (PerEQ(), "PerEQ", 5e-5, 1e-8),
#             (Exp(), "Exp", 1e-6, 1e-6),
#         ],
#         [(
#             k(α=α, β=β, l=l),
#             "$k_name(α=$(__foo(α)), β=$(__foo(β)), l=$(__foo(l)))",
#             1e-6,
#             1e-6,
#         )
#             for (k, k_name) in ((EQ, "EQ"), (Linear, "linear"), (Matern12, "exp"))
#             for α in (nothing, randn(rng))
#             for β in (nothing, exp(randn(rng)))
#             for l in (nothing, randn(rng))
#         ],
#     )
#         @testset "$name" begin
#             simple_gp_tests(_rng(), GP(k, GPC()), xs, σs; atol=atol, rtol=rtol)
#         end
#     end
# end

# @testset "FiniteGP (BlockDiagonal obs noise)" begin
#     rng, Ns = MersenneTwister(123456), [4, 5]
#     x = collect(range(-5.0, 5.0; length=sum(Ns)))
#     As = [randn(rng, N, N) for N in Ns]
#     Ss = [A' * A + I for A in As]

#     S = block_diagonal(Ss)
#     Smat = Matrix(S)

#     f = GP(cos, EQ(), GPC())
#     y = rand(f(x, S))

#     @test logpdf(f(x, S), y) ≈ logpdf(f(x, Smat), y)
#     adjoint_test(
#         (x, S, y)->logpdf(f(x, S), y), randn(rng), x, Smat, y;
#         atol=1e-6, rtol=1e-6,
#     )
#     adjoint_test(
#         (x, A1, A2, y)->logpdf(f(x, block_diagonal([A1 * A1' + I, A2 * A2' + I])), y),
#         randn(rng), x, As[1], As[2], y;
#         atol=1e-6, rtol=1e-6
#     )

#     @test elbo(f(x, Smat), y, f(x)) ≈ logpdf(f(x, Smat), y)
#     @test elbo(f(x, S), y, f(x)) ≈
#         elbo(f(x, Smat), y, f(x))
#     adjoint_test(
#         (x, A, y)->elbo(f(x, _to_psd(A)), y, f(x)),
#         randn(rng), x, randn(rng, sum(Ns), sum(Ns)), y;
#         atol=1e-6, rtol=1e-6,
#     )
#     adjoint_test(
#         (x, A1, A2, y) -> begin
#             S = block_diagonal([A1 * A1' + I, A2 * A2' + I])
#             return elbo(f(x, S), y, f(x))
#         end,
#         randn(rng), x, As[1], As[2], y;
#         atol=1e-6, rtol=1e-6,
#     )
# end


text:
parent(i::Int) = i ÷ 2
left(i::Int) = 2i
right(i::Int) = 2i + 1
parent(L::AbstractVector, i::Int) = L[parent(i)]
left(L::AbstractVector, i::Int) = L[left(i)]
right(L::AbstractVector, i::Int) = L[right(i)]

function exchange!(L::AbstractVector, i::Int, j::Int)
    #     println((i,j),": ", (L[i],L[j]))
    temp = L[i]
    L[i] = L[j]
    L[j] = temp
end

function minheapify!(L::AbstractVector, i::Int)
    min_ = i
    if left(i) <= length(L) && left(L, i) |> first < L[min_] |> first
        min_ = left(i)
    end
    if right(i) <= length(L) && right(L, i) |> first < L[min_] |> first
        min_ = right(i)
    end
    if min_ != i
        exchange!(L, i, min_)
        minheapify!(L, min_)
    end
end

function buildminheap!(L::AbstractVector)
    len = length(L)
    for i in len ÷ 2:-1:1
        minheapify!(L, i)
    end
end

heapmin(L::AbstractVector) = L[1]

function heapextractmin!(L::AbstractVector)
    head = L[1]
    tail = pop!(L)
    if length(L) > 0
        L[1] = tail
        minheapify!(L, 1)
    end
    return head
end

function heapdecreasekey!(L::AbstractVector, i::Int)
    while i > 1 && first(parent(L, i)) > first(L[i])
        exchange!(L, i, parent(i))
        i = parent(i)
    end
end

function minheapinsert!(L::AbstractVector, v::T) where T
    push!(L, v)
    heapdecreasekey!(L, length(L))
end

function heapchangekey!(L::AbstractVector{T}, i::Int, v::T, oldkey) where T
    #     @assert i<=length(L)
    if first(v) == oldkey
        nothing
    elseif first(v) < oldkey
        L[i] = v
        heapdecreasekey!(L, i)
    else
        L[i] = v
        minheapify!(L, i)
    end
end

heapchangekey!(L::AbstractVector{T}, i::Int, v::T) where T = heapchangekey!(L, i, v, first(L[i]))

function buildminheap2!(L::AbstractVector)
    len = length(L)
    for i in 1:len
        heapdecreasekey!(L, i)
    end
end

# 可回滚heap
HistoryType = Vector{Tuple{Int64,Int64}}
abstract type AbstractRecHeap{T} <: AbstractVector{T} end
mutable struct RecHeap{T} <: AbstractRecHeap{T}
    vector::Vector{T}
    history::HistoryType
    unused::Int
end

RecHeap() = RecHeap(Vector(), HistoryType(), 0)
RecHeap(v::AbstractVector) = RecHeap(v, HistoryType(), 0)
Base.getindex(L::RecHeap, i) =  getindex(L.vector, i)
Base.setindex!(L::RecHeap, v, i) = setindex!(L.vector, v, i)
Base.length(L::RecHeap) = length(L.vector) - L.unused
Base.size(L::RecHeap) = (length(L),)

function exchange!(L::AbstractRecHeap, i::Int, j::Int) end
function heapextractmin!(L::AbstractRecHeap) end
function undo!(L::AbstractRecHeap, nsteps::Int) end
function historylen(L::AbstractRecHeap) end
"仅恢复堆的结构，节点的内容不恢复，即无法撤销heapchangekey!带来的节点key或者value的变化"
undo!(L::AbstractRecHeap; tostep::Int) = undo!(L, historylen(L) - tostep)

function minheapinsert!(L::AbstractRecHeap, v::T, his::HistoryType) where T
    error("RecHeap doesn't allow insert nodes.")
end

function exchange!(L::RecHeap, i::Int, j::Int)
    push!(L.history, (i, j))
    temp = L[i]
    L[i] = L[j]
    L[j] = temp
end

function heapextractmin!(L::RecHeap)
    head = L[1]
    exchange!(L, 1, length(L))
    L.unused += 1
    push!(L.history, (-1, -1))
    if length(L) > 0
        minheapify!(L, 1)
    end
    return head
end

historylen(L::RecHeap) = length(L.history)
function undo!(L::RecHeap, nsteps::Int)
    for i in 1:nsteps
        e = pop!(L.history)
        if e == (-1, -1)
            L.unused -= 1
        else
            exchange!(L.vector, e...)
        end
    end
end
# 带位置追踪的heap
mutable struct tracHeap{T} <: AbstractRecHeap{T}
    heap::RecHeap{T}
    locL::AbstractVector
end

tracHeap() = tracHeap(RecHeap(), Vector())
tracHeap(v::AbstractVector) = tracHeap(RecHeap(v), [1:length(v)...])
tracHeap(v::AbstractVector, lv::AbstractVector) = tracHeap(RecHeap(v), lv)
Base.getindex(L::tracHeap, i) =  getindex(L.heap, i)
Base.setindex!(L::tracHeap, v, i) = setindex!(L.heap, v, i)
Base.length(L::tracHeap) = length(L.heap)
Base.size(L::tracHeap) = (length(L),)

function exchange!(L::tracHeap, i::Int, j::Int)
#     @assert i!=j
    exchange!(L.locL, last(L.heap[i]), last(L.heap[j]))
    exchange!(L.heap, i, j)
end

function heapextractmin!(L::tracHeap)
    head = L[1]
    exchange!(L, 1, length(L))
    L.heap.unused += 1
    push!(L.heap.history, (-1, -1))
    if length(L) > 0
        minheapify!(L, 1)
    end
    return head
end

historylen(L::tracHeap) = length(L.heap.history)

"仅恢复堆的结构，节点的内容不恢复，即无法撤销heapchangekey!带来的节点key或者value的变化"
function undo!(L::tracHeap, nsteps::Int)
    for i in 1:nsteps
        e = pop!(L.heap.history)
        if e == (-1, -1)
            L.heap.unused -= 1
        else
            exchange!(L.locL, last(L.heap[e[1]]), last(L.heap[e[2]]))
            exchange!(L.heap.vector, e...)
        end
    end
end


text:
println("hello world")

function sphere_vol(r)
    return 4/3*pi*r^3
end

quadratic(a,sqr_term,b) = (-b+sqr_term)/2a #2a == 2*a

function quadratic2(a::Float64, b::Float64, c::Float64)

    sqr_term = sqrt(Complex(b^2-4a*c))
    r1 = quadratic(a,sqr_term,b)
    r2 = quadratic(a,-sqr_term,b)

    return r1,r2 #if no return then last term is returned

end


text:
using Rackauckaspackage
using Documenter

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

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

deploydocs(;
    repo="github.com/atantos/Rackauckaspackage.jl",
)


text:
using Unitful, Gadfly, DataFrames

a = -9.81u"m/s^2"
t = (1:0.5:10)u"s"
v = a .* t
v2 = v .+ 2u"m/s"
h = 0.5 * a .* t.^2
df = DataFrame(time=t, velocity=v, position=h,
               unitlesst=ustrip.(t), unitlessv=ustrip.(v), unitlessh=ustrip.(h),
               position2=reverse(h))

# Test that it's possible to categorize by Unitful quantities
vstack(plot(df, x=:time, y=:position, color=:velocity, Geom.point),
       plot(df, x=:time, y=:velocity, color=:position, Geom.line),
       plot(df, x=:time, y=:position, color=:unitlessv))


text:
using Random
using LightGraphs
using ProgressMeter
using Distributions

function add_and_check!(hs, hc::Number, site::Int)::Bool
    hs[site] += 1
    if hs[site] >= hc
        return true
    else
        return false
    end
end

"""
add perturbation at random lattice site and let sytem relax until no further topplings
"""
function perturb_and_relax(rng::AbstractRNG, hs::Array{Int64}, lattice::SimpleGraph{Int64}, hc::Int64, p_dis::Float64;
        start_site = missing
    )
    topple_sites_over_time = Vector{Int}[]

    # select a random lattice site to perturb by adding extra grain
    if ismissing(start_site)
        site = rand(rng, vertices(lattice))
    else
        site = start_site
    end
    if add_and_check!(hs,hc,site)
        push!(topple_sites_over_time, [site])
    end

    # iterate over all topple sites that were created in last time step
    time = 1
    while length(topple_sites_over_time) >= time
        for site in topple_sites_over_time[time]
            # here happens the topple (add one grain to every neighboring sites)
            hs[site] -= hc
            # topple occasionally dissipates (required for lattice with periodic boundary conditions)
            if rand(rng) < 1-p_dis
                for neighbor_site in outneighbors(lattice, site)
                    if add_and_check!(hs,hc,neighbor_site)
                        #if true then this mean hs will topple (possible multiple counts here...)
                        if length(topple_sites_over_time) < time+1
                            push!(topple_sites_over_time, [neighbor_site])
                        elseif ! (neighbor_site in topple_sites_over_time[time+1])
                            push!(topple_sites_over_time[time+1], neighbor_site)
                        end
                    end
                end
            end
        end
        time += 1
    end

    return topple_sites_over_time
end

function add_external_input!(rng::AbstractRNG, topple_sites::Vector{Int}, hs, hc, P_input)
    num_input = rand(rng, P_input)
    for i in 1:num_input
        site = rand(rng, 1:length(hs))
        if add_and_check!(hs,hc,site)
            if ! (site in topple_sites)
                push!(topple_sites, site)
            end
        end
    end
end

function step!(rng::AbstractRNG, topple_sites::Vector{Int}, hs, lattice::SimpleGraph{Int64}, hc::Int64, p_dis::Float64)
    topple_sites_next_step = Int[]
    for site in topple_sites
        # here happens the topple (add one grain to every neighboring sites)
        hs[site] -= hc
        # topple occasionally dissipates (required for lattice with periodic
        # boundary conditions)
        if rand(rng) < 1-p_dis
            for neighbor_site in outneighbors(lattice, site)
                if add_and_check!(hs,hc,neighbor_site)
                    #avoid multiple occurances
                    if ! (neighbor_site in topple_sites_next_step)
                        push!(topple_sites_next_step, neighbor_site)
                    end
                end
            end
        end
    end
    return topple_sites_next_step
end


text:
"""
    PostStackWEM(m,d,adj;<keyword arguments>)

Post Stack Wave Equation Migration and Demigration of 3D isotropic data.

# Arguments
- `m` : filename of image
- `d` : filename of data
- `adj` : flag for adjoint (migration), or forward (demigration) (default=true)

**Keyword arguments**

- `vel = "NULL"` : seis file containing the velocity (should have same x and z dimensions as the desired image)
- `damping=1000`: damping for deconvolution imaging condition
- `nt=1`: trace length
- `ot=0`:trace origin
- `dt=0.001` : trace sample rate
- `wav="NULL"` :seis file containing the source wavelet (in time domain)
- `gz=0` : receiver depth
- `fmin=0` : min frequency to process (Hz)
- `fmax=99999` : max frequency to process (Hz)
- `padt=2` : pad factor for the time axis
- `padx=2` : pad factor for the spatial axes
- `verbose=true` : flag for error / debugging messages

"""
function PostStackWEM(m,d,adj;vel="NULL",damping=1000,nt=1,ot=0,dt=0.001,wav="NULL",gz=0,fmin=0,fmax=99999,padt=2,padx=2,verbose=true)

    v,h = SeisRead(vel)
    min_imx = h[1].imx
    max_imx = h[end].imx
    nmx = max_imx - min_imx + 1
    min_imy = h[1].imy
    max_imy = h[end].imy
    nmy = max_imy - min_imy + 1
    dmx = h[nmy+1].mx - h[nmy].mx > 0 ? h[nmy+1].mx - h[nmy].mx : 1.
    dmy = h[2].my - h[1].my > 0 ? h[2].my - h[1].my : 1.
    nz = int(h[1].n1)
    dz = h[1].d1
    oz = h[1].o1
    if (adj)
	d1,h = SeisRead(d)
	d1 = reshape(d1,nt,nmx,nmy)
	m1 = zowem(d1,param_zowem)
	for ix = 1 : nmx*nmy
	    h[ix].n1 = nz
	    h[ix].o1 = oz
	    h[ix].d1 = dz
	end
	SeisWrite(m,m1,h)
    else
	m1,h = SeisRead(m)
	m1 = reshape(m1,nz,nmx,nmy)
	d1 = zowem(m1,adj,omx=min_imx,dmx=dmx,nmx=nmx,omy=min_imy,dmy=dmy,nmy=nmy,oz=oz,dz=dz,nz=nz,ot=ot,dt=dt,nt=nt,padt=padt,padx=padx,v=v,fmin=fmin,fmax=fmax)
	for ix = 1 : nmx*nmy
	    h[ix].n1 = nt
	    h[ix].o1 = ot
	    h[ix].d1 = dt
	end
	SeisWrite(d,d1,h)
    end

    return
end

function zowem(in,adj;omx=0,dmx=1,nmx=1,omy=0,dmy=1,nmy=1,oz=0,dz=1,nz=1,ot=0,dt=1,nt=1,padt=2,padx=2,v=[],fmin=0.,fmax=999999.)

    nf = padt*nt
    nkx = nmx > 1 ? padx*nmx : 1
    nky = nmy > 1 ? pady*nmy : 1
    dkx = 2*pi/nkx/dmx;
    dky = 2*pi/nky/dmy;
    dw = 2*pi/nf/dt
    nw = int(nf/2) + 1
    if (adj)
	d = pad3d(in,nf,nkx,nky)
	m = zeros(Float32,nz,nkx,nky)
    else
	m = pad3d(in,nz,nkx,nky)
	d = zeros(Float32,nf,nkx,nky)
    end

    # decompose slowness into layer average, and layer perturbation
    po = zeros(Float32,nz)
    pd = zeros(Float32,nz,nmx*nmy)
    for iz = 1 : nz
	po[iz] = 1/mean(v[iz,:])
	pd[iz,:] = 1 ./ v[iz,:] - po[iz]
    end
    pd = reshape(pd,nz,nmx,nmy)
    pd = pad3d(pd,nz,nkx,nky)
    param["po"] = po
    param["pd"] = pd

    if(1 <= int(fmin*dt*nf) < nw)
	iw_min = int(floor(fmin*dt*nf))
    else
	iw_min = 1
    end
    if(fmax*dt*nf < nw)
	iw_max = int(floor(fmax*dt*nf))
	else
	    iw_max = int(floor(0.5/dt))
    end
    # set up kx and ky vectors
    kx = zeros(Float32,nkx,nky)
    ky = zeros(Float32,nkx,nky)
    for ikx = 1 : nkx
	kx[ikx,:] = ikx <= int(nkx/2) ? dkx*(ikx-1) : -(dkx*(nkx-1) - dkx*(ikx-1))
    end
    for iky = 1 : nky
	ky[:,iky] = iky <= int(nky/2) ? dky*(iky-1) : -(dky*(nky-1) - dky*(iky-1))
    end
    param["kx"] = kx
    param["ky"] = ky
    if (adj)
	D = fft(d,1)
	for iw=iw_min:iw_max
	    param["w"] = iw*dw
	    x = D[iw,:,:]
	    m,x = extrap1f(m,x,param)
	    D[iw,:,:] = x
	end
	return real(m[1:nz,1:nmx,1:nmy])
	else
	    D = zeros(Complex{Float32},nf,nkx,nky)
	for iw=iw_min:iw_max
	    param["w"] = iw*dw
	    x = D[iw,:,:]
	    m,x = extrap1f(m,x,adj,w=w,dz=dz,po=po,pd=pd,nz=nz,kx=kx,ky=ky)
	    D[iw,:,:] = x
	end
		# symmetries
	for iw=nw+1:nf
	    D[iw,:,:] = conj(D[nf-iw+2,:,:])
		end
	d = ifft(D,1)
	return real(d[1:nt,1:nmx,1:nmy])
    end

end

function extrap1f(m,d,adj,w=1,dz=1,po=[],pd=[],nz=1,kx=[],ky=[])

    if (adj)
	for iz = 1 : 1 : nz
	    d = FFTOp(d,true)
	    s = (w.^2)*(po[iz]*po[iz]) - kx.^2 - ky.^2
	    s[(s.<.0)] = 0.
	    d = d.*exp(1im*sqrt(s)*dz)
	    d = fft_op(d,false)
	    #d = d.*exp(1im*w*pd[iz]*dz)
	    m[iz,:,:] = m[iz,:,:] + 2*real(d)
	end
    else
	for iz = nz : -1 : 1
	    d += m[iz,:,:]
	    #d = d.*exp(-1im*w*pd[iz]*dz)
	    d = fft_op(d,true)
	    s = (w.^2)*(po[iz]^2) - kx.^2 - ky.^2
	    s[(s.<.0)] = 0.
	    d = d.*exp(-1im*sqrt(s)*dz)
	    d = fft_op(d,false)
	end
    end
    return m,d

end

function pad3d(a,N1,N2,N3)
    n1,n2,n3 = size(a)
    b = zeros(Float32,N1,N2,N3)
    b[1:n1,1:n2,1:n3] = a
    return b
end