text:
# Attributes that can be associated with parameters, function results, or the function
# itself.

export Attribute,
       EnumAttribute, StringAttribute,
       kind, value

abstract type Attribute end
reftype(::Type{T}) where {T<:Attribute} = API.LLVMAttributeRef

@checked struct EnumAttribute <: Attribute
    ref::reftype(Attribute)
end

@checked struct StringAttribute <: Attribute
    ref::reftype(Attribute)
end

# TODO: make @reftype's identify mechanism flexible enough to cover cases like this one,
#       and not only Value and Type

function Attribute(ref::API.LLVMAttributeRef)
    ref == C_NULL && throw(UndefRefError())
    if convert(Core.Bool, API.LLVMIsEnumAttribute(ref))
        return EnumAttribute(ref)
    elseif convert(Core.Bool, API.LLVMIsStringAttribute(ref))
        return StringAttribute(ref)
    else
        error("unknown attribute kind")
    end
end

function Base.show(io::IO, attr::T) where T<:Attribute
    print(io, "$T $(kind(attr))=$(value(attr))")
end

## enum attribute

# NOTE: the AttrKind enum is not exported in the C API,
#       so we don't expose a way to construct EnumAttribute from its raw enum value
#       (which also would conflict with the inner ref constructor)
function EnumAttribute(kind::String, value::Integer=0, ctx::Context=GlobalContext())
    enum_kind = API.LLVMGetEnumAttributeKindForName(kind, Csize_t(length(kind)))
    return EnumAttribute(API.LLVMCreateEnumAttribute(ref(ctx), enum_kind, UInt64(value)))
end

kind(attr::EnumAttribute) = API.LLVMGetEnumAttributeKind(ref(attr))

value(attr::EnumAttribute) = API.LLVMGetEnumAttributeValue(ref(attr))

## string attribute

StringAttribute(kind::String, value::String="", ctx::Context=GlobalContext()) =
    StringAttribute(API.LLVMCreateStringAttribute(ref(ctx), kind, length(kind),
                                                  value, length(value)))

function kind(attr::StringAttribute)
    len = Ref{Cuint}()
    data = API.LLVMGetStringAttributeKind(ref(attr), len)
    return unsafe_string(convert(Ptr{Int8}, data), len[])
end

function value(attr::StringAttribute)
    len = Ref{Cuint}()
    data = API.LLVMGetStringAttributeValue(ref(attr), len)
    return unsafe_string(convert(Ptr{Int8}, data), len[])
end


text:
module PandorasBoxEntityActivator

using ..Ahorn, Maple

@mapdef Entity "pandorasBox/entityActivator" EntityActivator(x::Integer, y::Integer, mode::String="ActivateInsideDeactivateOutside", activationMode::String="OnEnter", targets::String="", useTracked::Bool=true, flag::String="", updateInterval::Number=-1.0)

const placements = Ahorn.PlacementDict(
    "Entity Activator (Pandora's Box)" => Ahorn.EntityPlacement(
        EntityActivator,
        "rectangle"
    )
)

const modes = String[
    "ActivateInsideDeactivateOutside",
    "ActivateInside",
    "DeactivateInside",
    "ActivateOutside",
    "DeactivateOutside",
    "ActivateOnScreenDeactivateOffScreen"
]

const activationModes = String[
    "OnEnter",
    "OnStay",
    "OnLeave",
    "OnFlagActive",
    "OnFlagInactive",
    "OnFlagActivated",
    "OnFlagDeactivated",
    "OnUpdate",
    "OnAwake",
    "OnCameraMoved"
]

Ahorn.editingOptions(entity::EntityActivator) = Dict{String, Any}(
    "mode" => modes,
    "activationMode" => activationModes
)

Ahorn.minimumSize(entity::EntityActivator) = 8, 8
Ahorn.resizable(entity::EntityActivator) = true, true

Ahorn.selection(entity::EntityActivator) = Ahorn.getEntityRectangle(entity)

function Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::EntityActivator, room::Maple.Room)
    x = Int(get(entity.data, "x", 0))
    y = Int(get(entity.data, "y", 0))

    width = Int(get(entity.data, "width", 32))
    height = Int(get(entity.data, "height", 32))

    Ahorn.drawRectangle(ctx, 0, 0, width, height, (0.7, 1.0, 0.7, 0.4), (0.7, 1.0, 0.7, 1.0))
end

end

ext: jl

max_stars_repo_name: DemoJameson/PandorasBox

max_issues_repo_name: DemoJameson/PandorasBox

max_forks_repo_name: DemoJameson/PandorasBox


text:
export IsobeKakinuma

"""
    IsobeKakinuma(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
the Isobe-Kakinuma model proposed by [Isobe](https://doi.org/10.1061/9780784400890.023).

# Argument
`param` is of type `NamedTuple` (or a collection `NamedTuple`s) of and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.

## Optional keyword arguments
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `iterative`: solve the elliptic problem through GMRES if `true`, LU decomposition if `false` (default is `true`);
- `precond`: use a (left) preconditioner for GMRES if `true` (default), choose `precond` as the preconditioner if provided;
- `gtol`: relative tolerance of the GMRES algorithm (default is `1e-14`);
- `restart`: the corresponding option of the GMRES algorithm (default is `100`);
- `maxiter`: the corresponding option of GMRES (default is `nothing`);
- `ktol`: tolerance of the Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"Isobe-Kakinuma"`);

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `IsobeKakinuma.f!` to be called in explicit time-integration solvers;
2. a function `IsobeKakinuma.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `IsobeKakinuma.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct IsobeKakinuma <: AbstractModel

	label   :: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info	:: String

    function IsobeKakinuma(param::NamedTuple;
				mesh = Mesh(param),
				dealias = 0,
				ktol	= 0,
				iterate = true,
				gtol	= 1e-14,
				precond = true,
				restart	= nothing,
				maxiter	= nothing,
				label	= "Isobe-Kakinuma"
				)

		# Set up
		μ 	= param.μ
		ϵ 	= param.ϵ

		if isnothing(maxiter) maxiter = mesh.N end
		if isnothing(restart) restart = min(20,mesh.N) end

		# Print information
		info = "Isobe-Kakinuma model of order 2.\n"
		info *= "├─Shallowness parameter μ=$μ, nonlinearity parameter ϵ=$ϵ.\n"
		if dealias == 0
			info *= "├─No dealiasing. "
		else
			info *= "├─Dealiasing with Orszag's rule adapted to power $(dealias + 1) nonlinearity. "
		end
		if ktol == 0
			info *= "No Krasny filter. "
		else
			info *= "Krasny filter with tolerance $ktol."
		end
		if iterate == true
			if precond == false out="out" else out="" end
			info *= "\n└─Elliptic problem solved with GMRES method with$out preconditioning, \
			tolerance $gtol, maximal number of iterations $maxiter, restart after $restart iterations \
			(consider `iterate=false` for non-iterative method). "
		else
			info *= "\n└─Elliptic problem solved with standard LU factorization \
			(consider `iterate=true` for faster results). "
		end
		info *= "\nDiscretized with $(mesh.N) collocation points on [$(mesh.xmin), $(mesh.xmax)]."

		# Pre-allocate useful data
		k = mesh.k
		x 	= mesh.x
		x₀ = mesh.x[1]

		∂ₓ	=  1im * k
		if dealias == 0
			Π⅔ 	= ones(size(k)) # No dealiasing (Π⅔=Id)
		else
			K = (mesh.kmax-mesh.kmin)/(2+dealias)
			Π⅔ 	= abs.(k) .<= K # Dealiasing low-pass filter
		end
		FFT = exp.(-1im*k*(x.-x₀)');
        IFFT = exp.(1im*k*(x.-x₀)')/length(x);
        Id = Diagonal(ones(size(x)));

		if precond == true
			Precond = lu([  Id   μ*Id  ;
					1/2*Diagonal( ∂ₓ.^2) (Id + μ/10 * Diagonal(∂ₓ.^2 .* Π⅔) ) ])
		elseif precond == false
			Precond = Diagonal( ones( size( [k;k])) )
		else
			Precond = precond
		end

		h = zeros(Complex{Float64}, mesh.N)
		u, fftv, fftη, u, w, G = (similar(h),).*ones(6)
		C = similar([h ; h])
		guess = 0*C
		fftϕ = reshape(C,:,2)
		L = similar([FFT FFT ; FFT FFT])

		# Evolution equations are ∂t U = f(U)
		function f!(U)
			fftη .= U[:,1]
			h .= 1 .+ ϵ*ifft(fftη)
			fftv .= U[:,2]
			C .= [ zero(fftv) ; -1/2 * ∂ₓ.* fftv]
			if iterate == false
				L .= [Id    μ*(FFT * Diagonal( h.^2 ) * IFFT .* Π⅔)  ;
						1/2*Diagonal( ∂ₓ.^2 .* Π⅔) (Id + μ/10 * Diagonal(Π⅔) * FFT * Diagonal( h.^2 ) * IFFT * Diagonal( ∂ₓ.^2 .* Π⅔)) ]

				fftϕ .= reshape(L \ C,:,2)
			elseif iterate == true # does not work yet
				guess .= fftϕ[:]  # bof, le guess n'ameliore pas tellement les performances
		        function LL(fftϕ)
		            [fftϕ[1:end÷2] + μ*fft( (h.^2) .* ifft(  Π⅔ .* fftϕ[end÷2+1:end]))   ;
					1/2*∂ₓ.^2 .* Π⅔.* fftϕ[1:end÷2] .+ fftϕ[end÷2+1:end] .+ μ/10*Π⅔.*fft( (h.^2) .* ifft( ∂ₓ.^2 .* Π⅔ .* fftϕ[end÷2+1:end] ) ) ]
				end
				fftϕ .= reshape( gmres!( guess, LinearMap(LL, 2*length(x); issymmetric=false, ismutating=false) , C ;
						restart = restart, maxiter = maxiter, Pl = Precond, reltol = gtol ) , : ,2)
			end
			u .= ifft( ∂ₓ.* fftϕ[:,1] .+ fftv) .+ μ * (h.^2)  .* ifft(∂ₓ.* fftϕ[:,2])
			w .= 2*h.* ifft( fftϕ[:,2])
			G .= -∂ₓ.*Π⅔.* fft( h.* ifft( ∂ₓ.* fftϕ[:,1] .+ fftv) .+ μ*(h.^3)/3  .* ifft(∂ₓ.* fftϕ[:,2]) )

		   	U[:,1] .= G
			U[:,2] .= -∂ₓ .* (fftη .+ ϵ * Π⅔ .* fft( μ * w .* ifft(-G)
								.+ 1/2 * (u.^2 .+ μ * w.^2 ) ) )
			U[abs.(U).< ktol ].=0
		end

		# Build raw data from physical data.
		# Discrete Fourier transform with, possibly, dealiasing and Krasny filter.
		function mapto(data::InitialData)
			U = [Π⅔ .* fft(data.η(x)) Π⅔ .*fft(data.v(x))]
			U[abs.(U).< ktol ].=0
			return U
		end

		# Reconstruct physical variables from raw data
		# Return `(η,v,x)`, where
		# - `η` is the surface deformation;
		# - `v` is the derivative of the trace of the velocity potential;
		# - `x` is the vector of collocation points
		function mapfro(U)
			real(ifft(U[:,1])),real(ifft(U[:,2])),mesh.x
		end

        new(label, f!, mapto, mapfro, info)
    end
end


text:
# build test case for using Pose2Pose2 to landmarks -- to help AprilTag

using RoME, Distributions
using RoMEPlotting, Gadfly

# const TU = TransformUtils

# true landmark locations
global gtl = Dict{Symbol, Vector{Float64}}()
gtl[:l1] = [1.0;-1.0;0]
gtl[:l2] = [2.0;-1.0;0]
gtl[:l3] = [3.0;-1.0;0]
gtl[:l4] = [4.0;-1.0;0]

gtl[:l11] = [1.0;1.0;0]
gtl[:l12] = [2.0;1.0;0]
gtl[:l13] = [3.0;1.0;0]
gtl[:l14] = [4.0;1.0;0]

# true pose locations
global gtp = Dict{Symbol,Vector{Float64}}()
gtp[:x0] = [0.0;0;0]
gtp[:x1] = [0.5;0;0]
gtp[:x2] = [1.0;0;0]
gtp[:x3] = [2.0;0;0]
gtp[:x4] = [2.5;0;0]
gtp[:x5] = [3.0;0;0]

# calculate true measurements
global gtpt = Dict{Symbol, Vector{Float64}}()
for (ps,gp) in gtp, (ls,gl) in gtl
  gtpt[Symbol("$(ps)$(ls)")] = se2vee(SE2(gp) \ SE2(gl))
end

# Figure export folder
global currdirtime = now()
global imgdir = joinpath(ENV["HOME"], "Pictures", "testimgs", "$(currdirtime)")
mkdir(imgdir)

# construct the factor graph

global N=75
global fg = initfg()

addVariable!(fg, :x0, Pose2)
addFactor!(fg, [:x0], PriorPose2(MvNormal(zeros(3),Matrix(Diagonal([0.01;0.01;0.001].^2)))))

addVariable!(fg, :l1, Pose2)
addFactor!(fg, [:x0;:l1], Pose2Pose2(MvNormal(gtpt[:x0l1],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l11, Pose2)
addFactor!(fg, [:x0;:l11], Pose2Pose2(MvNormal(gtpt[:x0l11],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 0
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x1
addVariable!(fg, :x1, Pose2)
addFactor!(fg, [:x0;:x1], Pose2Pose2(MvNormal(zeros(3),Matrix(Diagonal([0.4;0.1;0.4].^2)))) )

addFactor!(fg, [:x1;:l1], Pose2Pose2(MvNormal(gtpt[:x1l1],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addFactor!(fg, [:x1;:l11], Pose2Pose2(MvNormal(gtpt[:x1l11],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l2, Pose2)
addFactor!(fg, [:x1;:l2], Pose2Pose2(MvNormal(gtpt[:x1l2],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l12, Pose2)
addFactor!(fg, [:x1;:l12], Pose2Pose2(MvNormal(gtpt[:x1l12],Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 1
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x2
addVariable!(fg, :x2, Pose2)
addFactor!(fg, [:x1;:x2], Pose2Pose2(MvNormal(zeros(3), Matrix(Diagonal([0.4;0.1;0.4].^2)))) )

addFactor!(fg, [:x2;:l2], Pose2Pose2(MvNormal(gtpt[:x2l2], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addFactor!(fg, [:x2;:l12], Pose2Pose2(MvNormal(gtpt[:x2l12], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l3, Pose2)
addFactor!(fg, [:x2;:l3], Pose2Pose2(MvNormal(gtpt[:x2l3], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 2
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# gtp[:x3] = [2.0;0;0]
addVariable!(fg, :x3, Pose2)
addFactor!(fg, [:x2;:x3], Pose2Pose2(MvNormal(zeros(3), Matrix(Diagonal([0.4;0.1;0.4].^2)))) )

addFactor!(fg, [:x3;:l3], Pose2Pose2(MvNormal(gtpt[:x3l3], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l13, Pose2)
addFactor!(fg, [:x3;:l13], Pose2Pose2(MvNormal(gtpt[:x3l13], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 3
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x4
# gtp[:x4] = [2.5;0;0]
addVariable!(fg, :x4, Pose2)
addFactor!(fg, [:x3;:x4], Pose2Pose2(MvNormal(zeros(3), Matrix(Diagonal([0.4;0.1;0.4].^2)))) )

addFactor!(fg, [:x4;:l3], Pose2Pose2(MvNormal(gtpt[:x4l3], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addFactor!(fg, [:x4;:l13], Pose2Pose2(MvNormal(gtpt[:x4l13], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l4, Pose2)
addFactor!(fg, [:x4;:l4], Pose2Pose2(MvNormal(gtpt[:x4l4], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addVariable!(fg, :l14, Pose2)
addFactor!(fg, [:x4;:l14], Pose2Pose2(MvNormal(gtpt[:x4l14], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 4
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

# pose :x5
# gtp[:x5] = [3.0;0;0]
addVariable!(fg, :x5, Pose2)
addFactor!(fg, [:x4;:x5], Pose2Pose2(MvNormal(zeros(3), Matrix(Diagonal([0.4;0.1;0.4].^2)))) )

addFactor!(fg, [:x5;:l4], Pose2Pose2(MvNormal(gtpt[:x5l4], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )
addFactor!(fg, [:x5;:l14], Pose2Pose2(MvNormal(gtpt[:x5l14], Matrix(Diagonal([0.1;0.1;0.01].^2)))) )

# writeGraphPdf(fg)

tree, smt, hist = solveTree!(fg)
# global tree = wipeBuildNewTree!(fg)
# inferOverTree!(fg, tree, N=N)

global psid = 5
global pl = drawPosesLandms(fg, spscale=0.1, drawhist=false)#,   meanmax=:mean,xmin=-3,xmax=6,ymin=-5,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"x$(psid).png"),30cm, 25cm),pl)
global pl = drawPosesLandms(fg, spscale=0.1)#,   meanmax=:mean,xmin=-3,xmax=3,ymin=-2,ymax=2);
Gadfly.draw(PNG(joinpath(imgdir,"hist_x$(psid).png"),30cm, 25cm),pl)

#


text:
export page, row, cell

"""
    `function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "", style::String = "",
                    channel::String = Genie.config.webchannels_default_route , head_content::String = "", kwargs...)`

Generates the HTML code corresponding to an SPA (a single page application), defining the root element of the Vue app.

### Example

```julia
julia> page(:elemid, [
        span("Hello", @text(:greeting))
        ])
"<!DOCTYPE html>\n<html><head><title></title><meta name=\"viewport\" content=\"width=device-width, initial-scale=1, maximum-scale=1, user-scalable=no, minimal-ui\" /></head><body class style><link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.min.css\" rel=\"stylesheet\" /><link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" /><div id=elemid><span v-text='greeting'>Hello</span></div><script src=\"/js/channels.js?v=1.17.1\"></script><script src=\"/js/stipple/underscore-min.js\"></script><script src=\"/js/stipple/vue.js\"></script><script src=\"/js/stipple/quasar.umd.min.js\"></script>\n<script src=\"/js/stipple/apexcharts.min.js\"></script><script src=\"/js/stipple/vue-apexcharts.min.js\"></script><script src=\"/js/stipple/stipplecore.js\" defer></script><script src=\"/js/stipple/vue_filters.js\" defer></script></body></html>"
```
"""
function page(elemid, args...; partial::Bool = false, title::String = "", class::String = "", style::String = "",
              channel::String = Genie.config.webchannels_default_route , head_content::String = "", core_theme::Bool = true,
              kwargs...)
  Stipple.Layout.layout(Genie.Renderer.Html.div(id = elemid, args...; class = class, kwargs...), partial = partial, title = title,
                        style = style, head_content = head_content, channel = channel, core_theme = core_theme)
end

"""
    `function row(args...; kwargs...)`

Creates a `div` HTML element with a CSS class named `row`. This works with Stipple's Twitter Bootstrap to create the
responsive CSS grid of the web page. The `row` function creates rows which should include `cell`s.

### Example

```julia
julia> row(span("Hello"))
"<div class=\"row\"><span>Hello</span></div>"
```
"""
function row(args...; kwargs...)
  kwargs = NamedTuple(Dict{Symbol,Any}(kwargs...), :class, "row")

  Genie.Renderer.Html.div(args...; kwargs...)
end

"""
    `function cell(args...; size::Int=0, kwargs...)`

Creates a `div` HTML element with CSS classes named `col col-12` and `col-sm-$size`.
This works with Stipple's Twitter Bootstrap to create the responsive CSS grid of the web page. The `cell`s should be
included within `row`s.

### Example

```julia
julia> row(cell(size=2, span("Hello")))
"<div class=\"row\"><div class=\"col col-12 col-sm-2\"><span>Hello</span></div></div>"
```
"""
function cell(args...; size::Int=0, kwargs...)
  kwargs = NamedTuple(Dict{Symbol,Any}(kwargs...), :class, "col col-12 col-sm$(size > 0 ? "-$size" : "")")

  Genie.Renderer.Html.div(args...; kwargs...)
end


text:
"""
    echo_time(x)
    echo_time!(x, val)

The echo time (TE) for the acquisition. This parameter is REQUIRED if corresponding
fieldmap data is present or the data comes from a multi echo sequence.
"""
@defprop EchoTime{:echo_time}::F64Second

"""
    inversion_time(x)
    inversion_time!(x, val)

Returns the inversion time (TI) for the acquisition, specified in seconds.
Inversion time is the time after the middle of inverting RF pulse to middle of
excitation pulse to detect the amount of longitudinal magnetization.
"""
@defprop InversionTime{:inversion_time}::F64Second

"""
    slice_timing(x)
    slice_timing!(x, val)

The time at which each slice was acquired within each volume (frame) of the
acquisition. Slice timing is not slice order -- rather, it is a list of times
(in JSON format) containing the time (in seconds) of each slice acquisition in
relation to the beginning of volume acquisition. The list goes through the
slices along the slice axis in the slice encoding dimension (see below). Note
that to ensure the proper interpretation of the `slice_timing` field, it is
important to check if the OPTIONAL `slice_encoding_direction` exists. In
particular, if `slice_encoding_direction` is negative, the entries in
`slice_timing` are defined in reverse order with respect to the slice axis
(i.e., the final entry in the `slice_timing` list is the time of acquisition of
slice 0). This parameter is REQUIRED for sparse sequences that do not have the
`delay_time` field set. In addition without this parameter slice time correction
will not be possible.
"""
@defprop SliceTiming{:slice_timing}::Vector{F64Second}

"""
    dwell_time(x)
    dwell_time!(x, val)

Actual dwell time (in seconds) of the receiver per point in the readout
direction, including any oversampling. For Siemens, this corresponds to DICOM
field (0019,1018) (in ns). This value is necessary for the optional readout
distortion correction of anatomicals in the HCP Pipelines. It also usefully
provides a handle on the readout bandwidth, which isn’t captured in the other
metadata tags. Not to be confused with `effective_echo_spacing`, and the frequent
mislabeling of echo spacing (which is spacing in the phase encoding direction)
as "dwell time" (which is spacing in the readout direction).
"""
@defprop DwellTime{:dwell_time}::F64Second

"""
    delay_time(x)
    delay_time!(x, val)

Returns the user specified time (in seconds) to delay the acquisition of data for
the following volume. If the field is not present it is assumed to be set to zero.
Corresponds to Siemens CSA header field ldelay_timeInTR. This field is REQUIRED
for sparse sequences using the `repetition_time` field that do not have the
slice_timing field set to allowed for accurate calculation of "acquisition time".
This field is mutually exclusive with `volume_timing`.
"""
@defprop DelayTime{:delay_time}::F64Second

"""
    acquisition_duration(x)
    acquisition_duration!(x, val)

Duration (in seconds) of volume acquisition. This field is REQUIRED for
sequences that are described with the `volume_timingfield` and that do not have the
`slice_timing` field set to allowed for accurate calculation of "acquisition time".
This field is mutually exclusive with `repetition_time`.
"""
@defprop AcquisitionDuration{:acquisition_duration}::F64Second

"""
    delay_after_trigger(x)
    delay_after_trigger!(x, val)

Returns duration (in seconds) from trigger delivery to scan onset. This delay is
commonly caused by adjustments and loading times. This specification is entirely
independent of `nvol_discarded_by_scanner` or `nvol_discarded_by_user`,
as the delay precedes the acquisition.
"""
@defprop DelayAfterTrigger{:delay_after_trigger}::F64Second

"""
    volume_timing(x)
    volume_timing!(x, val)

Returns the time at which each volume was acquired during the acquisition. It is
described using a list of times (in JSON format) referring to the onset of each
volume in the BOLD series. The list must have the same length as the BOLD
series, and the values must be non-negative and monotonically increasing. This
field is mutually exclusive with repetition_time and delay_time. If defined, this
requires acquisition time (TA) be defined via either slice_timing or
acquisition_duration be defined.
"""
@defprop VolumeTiming{:volume_timing}::Vector{F64Second}

"""
    repetition_time(x)
    repetition_time!(x, val)

Returns the time in seconds between the beginning of an acquisition of one volume
and the beginning of acquisition of the volume following it (TR). Please note that
this definition includes time between scans (when no data has been acquired) in
case of sparse acquisition schemes. This value needs to be consistent with the
pixdim[4] field (after accounting for units stored in xyzt_units field) in the
NIfTI header. This field is mutually exclusive with `volume_timing` and is derived
from DICOM Tag 0018, 0080 and converted to seconds.
"""
@defprop RepetitionTime{:repetition_time}::F64Second


text:
import StatsBase: sample

# Model selection

# Taken from https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2017.0009
function AIC(k::Int64, X::AbstractArray, Y::AbstractArray; likelyhood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return 2*k - 2*log(likelyhood(X, Y))
end
# Taken from https://royalsocietypublishing.org/doi/pdf/10.1098/rspa.2017.0009
function AICC(k::Int64, X::AbstractArray, Y::AbstractArray; likelyhood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return AIC(k, X, Y, likelyhood = likelyhood)+ 2*(k+1)*(k+2)/(size(X)[2]-k-2)
end

# Double check on that
# Taken from https://www.immagic.com/eLibrary/ARCHIVES/GENERAL/WIKIPEDI/W120607B.pdf
function BIC(k::Int64, X::AbstractArray, Y::AbstractArray; likelyhood = (X,Y) -> sum(abs2, X-Y))
    @assert size(X) == size(Y) "Dimensions of trajectories should be equal !"
    return - 2*log(likelyhood(X, Y)) + k*log(size(X)[2])
end

# Optimal Shrinkage for data in presence of white noise
# See D. L. Donoho and M. Gavish, "The Optimal Hard Threshold for Singular
# Values is 4/sqrt(3)", http://arxiv.org/abs/1305.5870
# Code taken from https://github.com/erichson/optht

function optimal_svht(m::Int64, n::Int64; known_noise::Bool = false)
    @assert m/n > 0
    @assert m/n <= 1

    β = m/n
    ω = (8*β) / (β+1+sqrt(β^2+14β+1))
    c = sqrt(2*(β+1)+ω)

    if known_noise
        return c
    else
        median = median_marcenko_pastur(β)
        return c / sqrt(median)
    end
end

function marcenko_pastur_density(t, lower, upper, beta)
    sqrt((upper-t).*(t-lower))./(2π*beta*t)
end

function incremental_marcenko_pastur(x, beta, gamma)
    @assert beta <= 1
    upper = (1+sqrt(beta))^2
    lower = (1-sqrt(beta))^2

    @inline marcenko_pastur(x) = begin
        if (upper-x)*(x-lower) > 0
            return marcenko_pastur_density(x, lower, upper, beta)
        else
            return zero(eltype(x))
        end
    end

    if gamma ≈ zero(eltype(gamma))
        i, ϵ = quadgk(x->(x^gamma)*marcenko_pastur(x), x, upper)
        return i
    else
        i, ϵ = quadgk(x->marcenko_pastur(x), x, upper)
        return i
    end
end

function median_marcenko_pastur(beta)
    @assert 0 < beta <= 1
    upper = (1+sqrt(beta))^2
    lower = (1-sqrt(beta))^2
    change = true
    x = ones(eltype(upper), 5)
    y = similar(x)
    while change && (upper - lower > 1e-5)
        x = range(lower, upper, length = 5)
        for (i,xi) in enumerate(x)
            y[i] = one(eltype(x)) - incremental_marcenko_pastur(xi, beta, 0)
        end
        any(y .< 0.5) ? lower = maximum(x[y .< 0.5]) : change = false
        any(y .> 0.5) ? upper = minimum(x[y .> 0.5]) : change = false
    end
    return (lower+upper)/2
end

function optimal_shrinkage(X::AbstractArray{T, 2}) where T <: Number
    m,n = minimum(size(X)), maximum(size(X))
    U, S, V = svd(X)
    τ = optimal_svht(m,n)
    inds = S .>= τ*median(S)
    return U[:, inds]*Diagonal(S[inds])*V[:, inds]'
end

function optimal_shrinkage!(X::AbstractArray{T, 2}) where T <: Number
    m,n = minimum(size(X)), maximum(size(X))
    U, S, V = svd(X)
    τ = optimal_svht(m,n)
    inds = S .>= τ*median(S)
    X .= U[:, inds]*Diagonal(S[inds])*V[:, inds]'
    return
end

function savitzky_golay(x::AbstractVector{T}, windowSize::Integer, polyOrder::Integer; deriv::Integer=0, dt::Real=1.0, crop::Bool = true) where T <: Number
	# Polynomial smoothing with the Savitzky Golay filters
	# Adapted from: https://github.com/BBN-Q/Qlab.jl/blob/master/src/SavitskyGolay.jl
	# More information: https://pdfs.semanticscholar.org/066b/7534921b308925f6616480b4d2d2557943d1.pdf
	# Requires LinearAlgebra and DSP modules loaded.

	# Some error checking
	@assert isodd(windowSize) "Window size must be an odd integer."
	@assert polyOrder < windowSize "Polynomial order must be less than window size."

	# Calculate filter coefficients
	filterCoeffs = calculate_filterCoeffs(windowSize, polyOrder, deriv, dt)

	# Pad the signal with the endpoints and convolve with filter
	halfWindow = Int(ceil((windowSize - 1)/2))
	paddedX = [x[1]*ones(halfWindow); x; x[end]*ones(halfWindow)]
	y = conv(filterCoeffs[end:-1:1], paddedX)

	if !crop
		# Return the valid midsection
		return y[2*halfWindow+1:end-2*halfWindow]
	else
		# Return cropped data. Excluding borders, where the estimation is less accurate
		return x[halfWindow+2:end-halfWindow-1], y[3*halfWindow+2:end-3*halfWindow-1]
	end
end

function savitzky_golay(x::AbstractMatrix{T}, windowSize::Integer, polyOrder::Integer; deriv::Integer=0, dt::Real=1.0, crop::Bool = true) where T <: Number
	# Polynomial smoothing with the Savitzky Golay filters
	# Adapted from: https://github.com/BBN-Q/Qlab.jl/blob/master/src/SavitskyGolay.jl
	# More information: https://pdfs.semanticscholar.org/066b/7534921b308925f6616480b4d2d2557943d1.pdf
	# Requires LinearAlgebra and DSP modules loaded.

	# Some error checking
	@assert isodd(windowSize) "Window size must be an odd integer."
	@assert polyOrder < windowSize "Polynomial order must be less than window size."

	# Calculate filter coefficients
	filterCoeffs = calculate_filterCoeffs(windowSize, polyOrder, deriv, dt)

	# Apply filter to each component
	halfWindow = Int(ceil((windowSize - 1)/2))

	if !crop
		y = similar(x)
		for (i, xi) in enumerate(eachrow(x))
			paddedX = [xi[1]*ones(halfWindow); xi; xi[end]*ones(halfWindow)]
			y₀ = conv(filterCoeffs[end:-1:1], paddedX)
			y[i,:] = y₀[2*halfWindow+1:end-2*halfWindow]
		end
		return y
	else
		cropped_x = x[:,halfWindow+2:end-halfWindow-1]
		y = similar(cropped_x)
		for (i, xi) in enumerate(eachrow(x))
			paddedX = [xi[1]*ones(halfWindow); xi; xi[end]*ones(halfWindow)]
			y₀ = conv(filterCoeffs[end:-1:1], paddedX)
			y[i,:] = y₀[3*halfWindow+2:end-3*halfWindow-1]
		end
		return cropped_x, y
	end
end

function calculate_filterCoeffs(windowSize::Integer, polyOrder::Integer, deriv::Integer, dt::Real)
	# Some error checking
	@assert isodd(windowSize) "Window size must be an odd integer."
	@assert polyOrder < windowSize "Polynomial order must be less than window size."

	# Form the design matrix A
	halfWindow = Int(ceil((windowSize - 1)/2))
	A = zeros(windowSize, polyOrder+1)
	for order = 0:polyOrder
		A[:, order+1] = (-halfWindow:halfWindow).^(order)
	end

	# Compute the required column of the inverse of A'*A
	# and calculate filter coefficients
	ei = zeros(polyOrder+1)
	ei[deriv+1] = 1.0
	inv_col = (A'*A) \ ei
	return A*inv_col * factorial(deriv) ./(dt^deriv)
end

@inline function burst_sampling(x::AbstractArray, samplesize::Int64, bursts::Int64)
    @assert size(x)[end] >= samplesize*bursts "Length of data array too small for subsampling of size $size!"
    inds = sample(1:size(x)[end]-samplesize, bursts, replace = false)
    inds = sort(unique(vcat([collect(i:i+samplesize) for i in inds]...)))
    return resample(x, inds)
end

@inline function burst_sampling(x::AbstractArray, y::AbstractArray, samplesize::Int64, bursts::Int64)
    @assert size(x)[end] >= samplesize*bursts "Length of data array too small for subsampling of size $size!"
    @assert size(x)[end] == size(y)[end]
    inds = sample(1:size(x)[end]-samplesize, bursts, replace = false)
    inds = sort(unique(vcat([collect(i:i+samplesize) for i in inds]...)))
    return resample(x, inds), resample(y, inds)
end

@inline function burst_sampling(x::AbstractArray, t::AbstractVector, period::T, bursts::Int64) where T <: AbstractFloat
    @assert period > zero(typeof(period)) "Sampling period has to be positive."
    @assert size(x)[end] == size(t)[end] "Provide consistent data."
    @assert bursts >= 1 "Number of bursts has to be positive."
    @assert t[end]-t[1]>= period*bursts "Bursting impossible. Please provide more data or reduce bursts."
    t_ids = zero(eltype(t)) .<= t .- period  .<= t[end] .- 2*period
    samplesize = Int64(floor(period/(t[end]-t[1])*length(t)))
    inds = sample(collect(1:length(t))[t_ids], bursts, replace = false)
    inds = sort(unique(vcat([collect(i:i+samplesize) for i in inds]...)))
    return resample(x, inds), resample(t, inds)
end

@inline function subsample(x::AbstractVector, frequency::Int64)
    @assert frequency > 1
    return x[1:frequency:end]
end

@inline function subsample(x::AbstractArray, frequency::Int64)
    @assert frequency > 1
    return x[:, 1:frequency:end]
end

@inline function subsample(x::AbstractArray, t::AbstractVector, period::T) where T <: AbstractFloat
    @assert period > zero(typeof(period)) "Sampling period has to be positive."
    @assert size(x)[end] == size(t)[end] "Provide consistent data."
    @assert t[end]-t[1]>= period "Subsampling impossible. Sampling period exceeds time window."
    idx = Int64[1]
    t_now = t[1]
    @inbounds for (i, t_current) in enumerate(t)
        if t_current - t_now >= period
            push!(idx, i)
            t_now = t_current
        end
    end
    return resample(x, idx), resample(t, idx)
end

@inline function resample(x::AbstractArray{T,1}, indx::AbstractArray{Int64}) where T <: Number
    @assert maximum(indx) <= length(x)
    @assert minimum(indx) >= 1
    return x[indx]
end

@inline function resample(x::AbstractArray{T,2}, indx::AbstractArray{Int64}) where T <: Number
    @assert maximum(indx) <= size(x, 2)
    @assert minimum(indx) >= 1
    return x[:, indx]
end


text:
module Freq

using ..LIKWID: LibLikwid

_ptr2str(ptr) = ptr == C_NULL ? "" : unsafe_string(ptr)

init() = LibLikwid.freq_init()

finalize() = LibLikwid.freq_finalize()

get_cpu_clock_current(cpuid::Integer) = LibLikwid.freq_getCpuClockCurrent(cpuid)
get_cpu_clock_max(cpuid::Integer) = LibLikwid.freq_getCpuClockMax(cpuid)
get_cpu_clock_min(cpuid::Integer) = LibLikwid.freq_getCpuClockMin(cpuid)
get_conf_cpu_clock_max(cpuid::Integer) = LibLikwid.freq_getConfCpuClockMax(cpuid)
get_conf_cpu_clock_min(cpuid::Integer) = LibLikwid.freq_getConfCpuClockMin(cpuid)

function set_cpu_clock_max(cpuid::Integer, freq::Integer)
    return LibLikwid.freq_setCpuClockMax(cpuid, freq)
end
function set_cpu_clock_min(cpuid::Integer, freq::Integer)
    return LibLikwid.freq_setCpuClockMin(cpuid, freq)
end

get_governor(cpuid::Integer) = _ptr2str(LibLikwid.freq_getGovernor(cpuid))
set_governor(cpuid::Integer, g::String) = LibLikwid.freq_setGovernor(cpuid, g)

get_avail_freq(cpuid::Integer) = _ptr2str(LibLikwid.freq_getAvailFreq(cpuid))
get_avail_govs(cpuid::Integer) = LibLikwid.freq_getAvailGovs(cpuid)

get_uncore_clock_min(s::Integer) = LibLikwid.freq_getUncoreFreqMin(s) * 1000000
set_uncore_clock_min(s::Integer, f::Integer) = LibLikwid.freq_setUncoreFreqMin(s, f)
get_uncore_clock_max(s::Integer) = LibLikwid.freq_getUncoreFreqMax(s) * 1000000
set_uncore_clock_max(s::Integer, f::Integer) = LibLikwid.freq_setUncoreFreqMax(s, f)

end # module


text:
using Stipple

export theme

const THEMES = Function[]

"""
    `function theme() :: String

Provides theming support for Stipple apps and pages. It includes Stipple's default CSS files and additional elements,
  in the form of HTML tags, can be injected by pushing to the `Stipple.Layout.THEMES` collection.

### Example

```julia
julia> theme()
"<link href=\"https://fonts.googleapis.com/css?family=Material+Icons\" rel=\"stylesheet\" /><link href=\"https://fonts.googleapis.com/css2?family=Lato:ital,wght@0,400;0,700;0,900;1,400&display=swap\" rel=\"stylesheet\" /><link href=\"/css/stipple/stipplecore.css\" rel=\"stylesheet\" />"

julia> StippleUI.theme()
"<link href=\"/css/stipple/quasar.min.css\" rel=\"stylesheet\" />"

julia> push!(Stipple.Layout.THEMES, StippleUI.theme)
```
"""
function theme(; core_theme::Bool = true) :: String
  output = ""

  if core_theme
    if ! Genie.Assets.external_assets(Stipple.assets_config)
      Genie.Router.route(Genie.Assets.asset_path(Stipple.assets_config, :css, file="stipplecore")) do
        Genie.Renderer.WebRenderable(
          Genie.Assets.embedded(Genie.Assets.asset_file(cwd=normpath(joinpath(@__DIR__, "..", "..")), type="css", file="stipplecore")),
          :css) |> Genie.Renderer.respond
      end
    end

    output *= string(
      stylesheet("https://fonts.googleapis.com/css?family=Material+Icons"),
      stylesheet(Genie.Assets.asset_path(Stipple.assets_config, :css, file="stipplecore"))
    )
  end

  string(output, join([f() for f in THEMES], "\n")) |> ParsedHTMLString
end


text:
using Q
import Q: kpointer, ktypecode, asarray, K_new, K_None, K_Ref
using Q._k, Q._k.GOT_Q
using Base.Test
using Base.Dates.AbstractTime
using DataFrames

include("test-utils.jl")
const side = GOT_Q ? "server" : "client"

@testset "Q tests" begin
  include("lowlevel-tests.jl")
  include("temporal-tests.jl")
  include("new-tests.jl")
  include("refcount-tests.jl")
  include("atom-tests.jl")
  include("vector-tests.jl")
  include("serialize-tests.jl")
  include("text-tests.jl")
  include("list-tests.jl")
  include("table-tests.jl")
  include("key-table-tests.jl")
  include("parser-tests.jl")
  include("$side-tests.jl")
  include("repl-tests.jl")
  include("expr-tests.jl")
  include("q-macro-tests.jl")
end


text:
# Cross-correlation tracker not present in AnnealingAPI
"""
    CCTracker(phase1, phase2)

Descriptor for cross-corelation function for the phases `phase1` and
`phase2`.

**NB:** Cross-correlation function does not commute,
i.e. `CCTracker(phase1, phase2)` and `CCTracker(phase2, phase1)` track
different functions.
"""
struct CCTracker{T} <: AbstractTracker{T}
    phase1 :: T
    phase2 :: T
end

(tracked :: CCTracker{T})(system :: AbstractArray{T}; kwargs...) where T =
    Directional.s2(system, Directional.SeparableIndicator(x -> x == tracked.phase1,
                                                          x -> x == tracked.phase2);
                   kwargs...)

const SimpleTracker{T}  = Union{L2Tracker{T}, S2Tracker{T}, CCTracker{T}}
const SurfaceTracker{T} = Union{SSTracker{T}, SVTracker{T}}
const CorrdataDict{T}   = Dict{AbstractTracker{T}, Directional.CorrelationData}

# Is gradient update needed?
update_gradient_p(:: SSTracker)       = true
update_gradient_p(:: SVTracker)       = true
update_gradient_p(:: AbstractTracker) = false

function gradient(array :: AbstractArray)
    deltas = imgradients(array, KernelFactors.sobel)
    return map((x...) -> norm(x), deltas...)
end

struct CorrelationTracker{T, N, A} <: AbstractArray{T, N}
    system     :: A
    periodic   :: Bool
    corrdata   :: CorrdataDict{T}
    grad       :: Array{Float64, N}
    fft_plans  :: Directional.S2FTPlans

    # For quick access
    corrlen    :: Int
    directions :: Vector{Symbol}
end

CDUpdateInfo{T} = Dict{Tuple{AbstractTracker{T}, Symbol}, Vector}
struct RollbackToken{T, N} <: AbstractRollbackToken
    val      :: T
    index    :: CartesianIndex{N}
    corrdata :: CDUpdateInfo{T}
    grad     :: Array{Float64, N}
end

@doc raw"""
    default_trackers(T)

Construct a vector of correlation functions which are tracked by
default (that is $S_2^1(x)$, $L_2^1(x)$ and $L_2^0(x)$). `T` is the
type of `x`.
"""
default_trackers(T :: Type) = [S2Tracker{T}(0), L2Tracker{T}(1), L2Tracker{T}(0)]

"""
    CorrelationTracker(system   :: AbstractArray{T, N};
                       tracking = default_trackers(T),
                       periodic = false[, directions][, kwargs...])

Create correlation functions tracker.

Create correlation tracker for the array `system`. `tracking` is a
vector of `AbstractTracker` types which specify correlation functions
you wish to track. `periodic` and `direction` have the same meaning as
in the most functions in `CorrelationFunctions.jl` package. Additional
arguments such as `len` may be passed in `kwargs`.

Returned tracker supports interface of `AbstractArray` (e.g. you can
perform element-wise read and write operations).

# Examples
```jldoctest
julia> let
       system = rand(MersenneTwister(35), 0:1, (30, 10))
       tracker = CorrelationTracker(system)
       end
30×10 CorrelationTracker{Int64, 2, Matrix{Int64}}:
 0  1  0  1  1  0  0  1  1  0
 1  1  1  0  0  0  0  0  1  1
 0  0  0  0  0  0  1  1  0  1
 1  1  1  0  1  1  1  0  1  0
 0  1  0  0  1  0  0  1  1  1
 0  0  0  0  0  0  1  0  1  1
 0  0  1  0  1  1  0  1  0  1
 1  0  0  1  0  0  1  0  1  0
 0  1  1  0  0  1  1  1  1  1
 0  0  1  1  1  1  0  0  0  0
 0  0  1  1  0  0  1  1  1  0
 0  1  0  0  0  1  0  0  1  0
 1  0  0  1  0  0  1  1  0  1
 0  1  0  1  0  0  1  1  1  0
 1  1  0  1  1  1  0  1  0  1
 1  1  1  0  0  0  0  1  0  1
 1  0  0  1  0  0  1  1  1  0
 0  0  0  1  0  0  0  1  1  0
 1  0  1  0  1  0  0  0  1  0
 1  0  0  1  0  0  0  0  0  1
 1  1  1  0  1  0  1  0  1  1
 0  1  0  1  1  0  0  1  0  1
 0  0  0  1  0  0  1  1  1  1
 0  0  1  1  1  1  0  1  1  0
 1  0  1  1  0  0  0  0  0  1
 1  1  0  1  0  1  1  0  1  0
 0  1  1  0  0  1  1  0  1  0
 0  1  0  0  1  0  0  1  0  0
 1  1  0  0  1  1  0  0  0  1
 0  0  1  1  0  1  1  1  1  0
```
"""
function CorrelationTracker(system     :: AbstractArray{T, N};
                            tracking   :: Vector{<:AbstractTracker{T}} = default_trackers(T),
                            periodic   :: Bool           = false,
                            directions :: Vector{Symbol} = system |> Directional.default_directions,
                            kwargs...) where {T, N}
    corrdata =
        CorrdataDict{T}(data => data(system;
                                     periodic   = periodic,
                                     directions = directions,
                                     kwargs...)
                        for data in tracking)
    len = length(first(corrdata)[2])
    # FIXME: What about multiphase systems?
    return CorrelationTracker{T, N, typeof(system)}(
        copy(system), periodic,
        corrdata, gradient(system),
        Directional.S2FTPlans(system, periodic),
        len, directions)
end

function update_cf(tracker :: CorrelationTracker{T, N},
                   data    :: SimpleTracker{T},
                   val,
                   index   :: CartesianIndex{N}) where {T, N}
    corrdata = tracker.corrdata[data]
    dict     = Dict{Symbol, Vector}()

    for direction in tracker.directions
        slice = get_slice(tracker.system,
                          tracker.periodic,
                          Tuple(index), direction)
        scorr = data(slice;
                     plans    = tracker.fft_plans,
                     periodic = tracker.periodic,
                     len      = tracker.corrlen)
        dict[direction] = scorr.success[:x]
    end

    return dict
end

function update_at_point(tracker   :: CorrelationTracker{T, N},
                         index     :: CartesianIndex{N},
                         direction :: Symbol,
                         _         :: SSTracker{T}) where {T, N}
    slice = get_slice(tracker.grad,
                      tracker.periodic,
                      Tuple(index), direction)
    s2 = Directional.s2(slice, Directional.SeparableIndicator(identity);
                        periodic = tracker.periodic,
                        plans    = tracker.fft_plans,
                        len      = tracker.corrlen)
    return s2.success[:x]
end

function update_at_point(tracker   :: CorrelationTracker{T, N},
                         index     :: CartesianIndex{N},
                         direction :: Symbol,
                         _         :: SVTracker{T}) where {T, N}
    slice_surface = get_slice(tracker.grad,
                              tracker.periodic,
                              Tuple(index), direction)
    slice_system = get_slice(tracker.system,
                             tracker.periodic,
                             Tuple(index), direction)

    χ1(x) = slice_system[x] == 0
    χ2(x) = slice_surface[x]
    s2 = Directional.s2(CartesianIndices(slice_system),
                        Directional.SeparableIndicator(χ1, χ2);
                        periodic = tracker.periodic,
                        plans    = tracker.fft_plans,
                        len      = tracker.corrlen)
    return s2.success[:x]
end

function update_cf(tracker :: CorrelationTracker{T, N},
                   data    :: SurfaceTracker{T},
                   val,
                   index   :: CartesianIndex{N}) where {T, N}
    corrdata = tracker.corrdata[data]
    dict     = Dict{Symbol, Vector}()

    indices    = CartesianIndices(tracker)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(fidx)

    for direction in tracker.directions
        segment_seeds = Vector{CartesianIndex{N}}(undef, 0)
        for idx in max(index - uidx, fidx):min(index + uidx, lidx)
            if !any(x -> same_line_p(x, idx, direction), segment_seeds)
                push!(segment_seeds, idx)
            end
        end

        success = mapreduce(.+, segment_seeds) do idx
            update_at_point(tracker, idx, direction, data)
        end

        dict[direction] = success
    end

    return dict
end

function update_gradient!(tracker  :: CorrelationTracker{T, N},
                          index    :: CartesianIndex{N}) where {T, N}
    grad   = tracker.grad
    system = tracker.system

    indices    = CartesianIndices(grad)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(index)

    # Select 5x5 square (or 5x5x5 cube) used to update gradient
    gradstart = max(index - 2uidx, fidx)
    gradstop  = min(index + 2uidx, lidx)
    subsys    = system[gradstart:gradstop]
    subgrad   = gradient(subsys)

    # Index of the updated element in subgrad
    sindex       = index - gradstart + uidx
    sindices     = CartesianIndices(subgrad)
    sfidx, slidx = first(sindices), last(sindices)
    substart     = max(sindex - uidx, sfidx)
    substop      = min(sindex + uidx, slidx)

    # Update gradient
    upstart = max(index - uidx, fidx)
    upstop  = min(index + uidx, lidx)

    oldgrad = grad[upstart:upstop]
    grad[upstart:upstop] .= subgrad[substart:substop]

    return oldgrad
end

function rollback_gradient!(tracker :: CorrelationTracker{T, N},
                            index   :: CartesianIndex{N},
                            subgrad :: Array{Float64, N}) where {T, N}
    grad       = tracker.grad
    indices    = CartesianIndices(grad)
    fidx, lidx = first(indices), last(indices)
    uidx       = oneunit(index)
    upstart    = max(index - uidx, fidx)
    upstop     = min(index + uidx, lidx)

    grad[upstart:upstop] .= subgrad

    return nothing
end

function AnnealingAPI.update_corrfns!(tracker :: CorrelationTracker{T,N},
                                      val,
                                      index   :: CartesianIndex{N}) where {T, N}
    trackers = keys(tracker.corrdata)
    update_info = CDUpdateInfo{T}()
    directions = tracker.directions

    # Do everything we can before updating the value in the underlying array
    for tr in trackers
        success = update_cf(tracker, tr, val, index)
        for direction in directions
            update_info[(tr, direction)] = -success[direction]
        end
    end

    # Change the value
    oldval = tracker.system[index]
    tracker.system[index] = val

    # Update gradient if needed
    grad = Array{Float64, N}(undef, (0 for n in 1:N)...)
    if any(update_gradient_p, trackers)
        grad = update_gradient!(tracker, index)
    end

    # Do everything else
    for tr in trackers
        corrdata = tracker.corrdata[tr]
        success = update_cf(tracker, tr, val, index)
        for direction in directions
            update_info[(tr, direction)] .+= success[direction]
            corrdata.success[direction] .+= update_info[(tr, direction)]
        end
    end

    return RollbackToken{T, N}(oldval, index, update_info, grad)
end

function AnnealingAPI.rollback!(tracker  :: CorrelationTracker{T, N},
                                rollback :: RollbackToken{T, N}) where {T, N}
    trackers = keys(tracker.corrdata)
    index    = rollback.index
    val      = rollback.val

    # Roll back correlation functions
    for tr in trackers
        corrdata = tracker.corrdata[tr]

        for direction in Directional.directions(corrdata)
            corrdata.success[direction] .-= rollback.corrdata[(tr, direction)]
        end
    end

    # Roll back the value of the changed element
    tracker.system[index] = val

    # Roll back gradient
    if any(update_gradient_p, trackers)
        rollback_gradient!(tracker, index, rollback.grad)
    end

    return tracker
end

# Interface
AnnealingAPI.tracked_data(x :: CorrelationTracker) = x.corrdata |> keys
AnnealingAPI.tracked_length(x :: CorrelationTracker) = x.corrlen
AnnealingAPI.tracked_directions(x :: CorrelationTracker) = x.directions

# Array interface
Base.size(x :: CorrelationTracker) = size(x.system)
Base.getindex(x :: CorrelationTracker, idx :: Vararg{Int}) = getindex(x.system, idx...)
Base.setindex!(x :: CorrelationTracker{T}, val, idx :: Vararg{Int}) where T =
    AnnealingAPI.update_corrfns!(x, val, CartesianIndex(idx))


text:
@testset "Mixed measurements (vs pf)" begin

    # set measurement path
    msr_path = joinpath(mktempdir(),"temp.csv")

    @testset "ACP-rwlav" begin
        # set model
        crit     = "rwlav"
        model    = _PMD.ACPUPowerModel #this is going to be used for the SE

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
         data = _PMD.transform_data_model(data);

        for data_model in [_PMD.ACRUPowerModel,_PMD.IVRUPowerModel]
            # solve the power flow
            pf_result = _PMD.solve_mc_pf(data, data_model, ipopt_solver)

            # write measurements based on power flow
            _PMDSE.write_measurements!(data_model, data, pf_result, msr_path, exclude = ["vi","vr"])

            # read-in measurement data and set initial values
            _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
            _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
            _PMDSE.assign_start_to_variables!(data)
            _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
            # set se settings
            data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                               "rescaler" => 1)

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max, 0.0; atol = 5e-5)
            @test isapprox(avg, 0.0; atol = 1e-7)
        end
    end
    @testset "ACP-WLS" begin
        # set model
        crit     = "wls"
        model    = _PMD.ACPUPowerModel
        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        for data_model in [_PMD.ACRUPowerModel,_PMD.IVRUPowerModel]
            # solve the power flow
            pf_result = _PMD.solve_mc_pf(data, data_model, ipopt_solver)

            # write measurements based on power flow
            _PMDSE.write_measurements!(data_model, data, pf_result, msr_path, exclude = ["vi","vr"])

            # read-in measurement data and set initial values
            _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
            _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
            _PMDSE.assign_start_to_variables!(data)
            _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
            # set se settings
            data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                               "rescaler" => 1)

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max, 0.0; atol = 5e-5)
            @test isapprox(avg, 0.0; atol = 2e-6)
        end
    end
    @testset "ACR-rwlav" begin
        # set model
        crit     = "rwlav"
        model    = _PMD.ACRUPowerModel
        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        for data_model in [_PMD.ACPUPowerModel,_PMD.IVRUPowerModel]
            # solve the power flow
            pf_result = _PMD.solve_mc_pf(data, data_model, ipopt_solver)

            # write measurements based on power flow
            _PMDSE.write_measurements!(data_model, data, pf_result, msr_path)

            # read-in measurement data and set initial values
            _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
            _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
            _PMDSE.assign_start_to_variables!(data)
            _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )

            # set se settings
            data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                               "rescaler" => 1)

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max, 0.0; atol = 5e-5)
            @test isapprox(avg, 0.0; atol = 1e-7)
        end
    end
    @testset "ACR-WLS" begin
        # set model
        crit     = "wls"
        model    = _PMD.ACRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        for data_model in [_PMD.ACPUPowerModel,_PMD.IVRUPowerModel]
            # solve the power flow
            pf_result = _PMD.solve_mc_pf(data, data_model, ipopt_solver)

            # write measurements based on power flow
            _PMDSE.write_measurements!(data_model, data, pf_result, msr_path)

            # read-in measurement data and set initial values
            _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
            _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
            _PMDSE.assign_start_to_variables!(data)
            _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
            # set se settings
            data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                               "rescaler" => 10)

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max, 0.0; atol = 5e-5)
            @test isapprox(avg, 0.0; atol = 1e-7)
        end
    end
    @testset "IVR-rwlav" begin
        # set model
        crit     = "rwlav"
        model    = _PMD.IVRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        for data_model in [_PMD.ACPUPowerModel,_PMD.ACRUPowerModel]
            # solve the power flow
            pf_result = _PMD.solve_mc_pf(data, data_model, ipopt_solver)

            # write measurements based on power flow
            _PMDSE.write_measurements!(data_model, data, pf_result, msr_path)

            # read-in measurement data and set initial values
            _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
            _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
            _PMDSE.assign_start_to_variables!(data)
            _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
            # set se settings
            data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                               "rescaler" => 1)

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max, 0.0; atol = 5e-5)
            @test isapprox(avg, 0.0; atol = 1e-7)
        end
    end
    @testset "IVR-WLS" begin
        # set model
        crit     = "wls"
        model    = _PMD.IVRUPowerModel

        # load data
        data = _PMD.parse_file(_PMDSE.get_enwl_dss_path(ntw, fdr))
        if rm_transfo _PMDSE.rm_enwl_transformer!(data) end
        if rd_lines   _PMDSE.reduce_enwl_lines_eng!(data) end

        # insert the load profiles
        _PMDSE.insert_profiles!(data, season, elm, pfs, t = time_step)

        # transform data model
        data = _PMD.transform_data_model(data);

        for data_model in [_PMD.ACPUPowerModel,_PMD.ACRUPowerModel]
            # solve the power flow
            pf_result = _PMD.solve_mc_pf(data, data_model, ipopt_solver)

            # write measurements based on power flow
            _PMDSE.write_measurements!(data_model, data, pf_result, msr_path)

            # read-in measurement data and set initial values
            _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
            _PMDSE.add_voltage_measurement!(data, pf_result, 0.005)
            _PMDSE.assign_start_to_variables!(data)
            _PMDSE.update_all_bounds!(data; v_min = 0.8, v_max = 1.2, pg_min=-1.0, pg_max = 1.0, qg_min=-1.0, qg_max=1.0, pd_min=-1.0, pd_max=1.0, qd_min=-1.0, qd_max=1.0 )
            # set se settings
            data["se_settings"] = Dict{String,Any}("criterion" => crit,
                                               "rescaler" => 10)

            # solve the state estimation
            se_result = _PMDSE.solve_mc_se(data, model, ipopt_solver)

            # tests
            delta, max, avg = _PMDSE.calculate_voltage_magnitude_error(se_result, pf_result)
            @test isapprox(max, 0.0; atol = 5e-5)
            @test isapprox(avg, 0.0; atol = 4e-7)
        end
    end
end


text:
##
# @author b33rNc0d3
# @date 2020-07-05
# Linear Regression in Julia with Turing for Bayesian Statistics: Techniques and Models
##

using Turing
using Distributions


text:
struct Game{GameState <: AbstractGameState, Renderer <: AbstractRenderer, InputHandler <: AbstractInputHandler, Config <: AbstractConfig}
    gamestate::GameState
    renderer::Renderer
    inputhandler::InputHandler
    config::Config
end

function init_game()
    # do basic shit first
    config = NoConfig()
    gamestate = init_game_state()
    renderer = init_renderer()
    inputhandler = generate_input_handler(scene(renderer).events, defaultkeyboardmapping())
    Game(gamestate, renderer, inputhandler, config)
end

function run_game(game::Game)
    newlantern = DistanceGoalLightHouse(Complex(0.7, 0.7), Complex(-1.4, 0.0))
    lasttime = time()

    stoprenderloopindex = game.inputhandler.playerintentindex[stoprenderloop]
    plantlighthouseindex = game.inputhandler.playerintentindex[plantlighthouse]
    alreadychangedlanatern = false
    while true
        if game.inputhandler.currentintent[stoprenderloopindex].intended
            break
        else
            if game.inputhandler.currentintent[plantlighthouseindex].intended && !alreadychangedlanatern
                newlantern = DistanceGoalLightHouse(Complex(randn(Float64, 2)...), Complex(-1.4, 0.0))
                @show newlantern.startpoint
                alreadychangedlanatern = true
            elseif !game.inputhandler.currentintent[plantlighthouseindex].intended
                alreadychangedlanatern = false
            end
            nowtime = time()
            deltatime = nowtime - lasttime
            lasttime = nowtime
            timestep_game_state(game.gamestate, game.inputhandler.currentintent, game.gamestate.worldstate.current_transform, deltatime)

            render_game(game.renderer, game.gamestate)
            game.gamestate.worldstate.current_transform = evolvetransform(newlantern, game.gamestate.worldstate.current_transform, game.gamestate.player_state)
        end
    end
    game
end

max_stars_repo_name: zoemcc/fractrixzeta

max_issues_repo_name: zoemcc/fractrixzeta

max_forks_repo_name: zoemcc/fractrixzeta


text:
using Test
using Forecast
import Forecast: fvar

@testset "forecast_ar" begin

    Φ,Σ2,n = .991,1.,10
    @test fvar(Φ,Σ2,n) ≈ [1.0, 1.982081, 2.946564090561, 3.8937646086222375, 4.823992240600336, 5.737551123641018, 6.634739945056495, 7.515852039981027, 8.381175487276607, 9.230993203720098]

    Φ,Σ,n = [1],1,10
    @test fvar(Φ,Σ,n) == [1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]

    Φ,Σ2,n = [.1,.2,.3,.4],1.,10
    @test fvar(Φ,Σ2,n) ≈ [1.0, 1.3, 1.3900000000000001, 1.4170000000000003, 1.4251, 1.42753, 1.4282590000000002, 1.4284777000000002, 1.4285433100000002, 1.4285629930000001]

    Φ,Σ2,n = [.1 .2;.3 .4],[1. 0.; 0. 1.],10
    @test fvar(Φ,Σ2,n) ≈ [1.0 1.0; 1.05 1.25; 1.0649 1.3209; 1.069185 1.341385; 1.07042201 1.34729681; 1.0707790225 1.3490030505; 1.070882061449 1.3494954953290001; 1.07091179993665 1.34963762183705; 1.07092038288254 1.3496786415463762; 1.0709228600414522 1.3496904804114789]

    Φ,Σ2,n = reshape([[.1 .2; .3 .4] [.5 .6; .7 .8]],2,2,2),[.9 .10; .10 .9],10
    @test fvar(Φ,Σ2,n) ≈ [0.9 0.9; 1.896 3.3360000000000003; 5.307360000000001 11.574240000000001; 16.897152000000002 39.571737600000006; 56.28080870400001 134.71025894400003; 190.11105110016 458.00206718976006; 644.8818211024897 1556.5853271158787; 2190.2459254584123 5289.699343417738; 7441.573347348584 17975.25606597047; 25286.196251506655 61082.256992946226]

    # ar(Φ,n)
    Φ,n = .5,n
    Φ,Φ0,x0,n = .1,.2,.3,10
    xar = ar(arsim(Φ,Φ0,x0,n),1)
    fxar = forecast(xar,10)
    @test size(fxar.mean,1) == 10

end


text:
# output functions

@doc """
### SimpleODEMod.printSolution

  This function prints the solution vector u to a file named solution.dat.
  Each line contains: element number, local node number, and the 4 conservative
  variable values.

  Inputs:
   * mesh :: AbstractMesh
   * u  : vector of values
"""->
function printSolution(mesh::AbstractMesh, u::AbstractVector)
# print solution to file
# format = for each element, node, print u rho*u rho*v E

myrank = mesh.myrank
println("entered printSolution")
f = open("solution_$myrank.dat", "w")

for i=1:mesh.numEl
  dofnums_i = getGlobalNodeNumbers(mesh, i)

  for j=1:3
    u_vals = u[dofnums_i[:, j]]
    str = @sprintf("%d %d %16.15e \n", i, j, u_vals[1])  # print element number

    print(f, str)
  end
#  print(f, "\n")
end

close(f)
return nothing

end

@doc """
### SimpleODEMod.printCoordinates

  This function prints the mesh vertex coordinates in the same format
  as printSolution prints the solution vector.

  Inputs:
    * mesh :: AbstractMesh
"""->
function printCoordinates(mesh::AbstractMesh)
# print solution to file
# format = for each element, node, print u rho*u rho*v E

myrank = mesh.myrank
println("entered printCoordinates")
writedlm("coords_output_$myrank.dat", mesh.coords)

return nothing
end

#TODO: automate adding _$myrank to file names
@doc """
### SimpleODEMod.printSolution

  This function prints the solution vector to the file with the given name.
  The solution vector is printed one value per line.

  Inputs:
    * name : AbstractString naming file, including extension
    * u  : vector to print
"""->
function printSolution(name::AbstractString, u::AbstractVector)

  f = open(name, "a+")

  for i=1:length(u)
    write(f, string(u[i], "\n"))
  end

  close(f)

  return nothing
end

@doc """
### SimpleODEMod.printMatrix

  This function prints a 2D matrix of Float64s to a file with the given name

  Inputs:
    * name  : AbstractString naming file, including extension
    * u  : matrix to be printed
"""->
function printMatrix(name::AbstractString, u::AbstractArray{Float64, 2})
# print a matrix to a file, in a readable format

(m,n) = size(u)

f = open(name, "a+")

for i=1:m
  for j=1:n
    str = @sprintf("%16.15e", u[i,j])
    print(f, str)
    if j < n  # don't put space at end of line
      print(f, " ")
    end
  end

  print(f, "\n")
end

close(f)
return nothing

end

function printMatrix(name::AbstractString, u::AbstractArray{T, 3}) where T
# print a matrix to a file, in a quasireadable format

#println("printing 3d matrix")
(tmp, m,n) = size(u)
#println("numel = ", n)
#println("nnodes = ", m)

f = open(name, "a+")

for i=1:n  # loop over elements
  for j=1:m  # loop over nodes on element
#    str = @sprintf("%16.15e", u[i,j])
    println(f, "el ", i, " node ", j, " flux = ", u[:, j, i] )
#    if j < n  # don't put space at end of line
#      print(f, " ")
#    end
  end

#  print(f, "\n")
end

close(f)
return nothing

end


text:
export Boussinesq

"""
    Boussinesq(param;kwargs)

Define an object of type `AbstractModel` in view of solving the initial-value problem for
`abcd`-Boussinesq models (with `b=d` and `c=0`).
See [Bona, Chen, and Saut](https://doi.org/10.1007/s00332-002-0466-4)

# Argument
`param` is of type `NamedTuple` and must contain
- dimensionless parameters `ϵ` (nonlinearity) and `μ` (dispersion);
- numerical parameters to construct the mesh of collocation points, if `mesh` is not provided as a keyword argument.

## Optional keyword arguments
- two parameters `a` (default is `-1/3`) and `b` (default is `+1/3`) which determine the model solved. You need `a+2*b=1/3` for validity as a long wave model (without surface tension).
- `mesh`: the mesh of collocation points. By default, `mesh = Mesh(param)`;
- `ktol`: tolerance of the low-pass Krasny filter (default is `0`, i.e. no filtering);
- `dealias`: dealiasing with Orlicz rule `1-dealias/(dealias+2)` (default is `0`, i.e. no dealiasing);
- `label`: a label for future references (default is `"Boussinesq"`);

# Return values
Generate necessary ingredients for solving an initial-value problem via `solve!`:
1. a function `Boussinesq.f!` to be called in explicit time-integration solvers;
2. a function `Boussinesq.mapto` which from `(η,v)` of type `InitialData` provides the raw data matrix on which computations are to be executed;
3. a function `Boussinesq.mapfro` which from such data matrix returns the Tuple of real vectors `(η,v,x)`, where
    - `η` is the values of surface deformation at collocation points `x`;
    - `v` is the derivative of the trace of the velocity potential at `x`.

"""
mutable struct Boussinesq <: AbstractModel

	label 	:: String
	f!		:: Function
	mapto	:: Function
	mapfro	:: Function
	info	:: String

    function Boussinesq(param::NamedTuple;
						mesh = Mesh(param),
						a=-1//3,b=1//3,
						dealias=0,ktol=0,
						label="Boussinesq"
						)

		m=WhithamBoussinesq(param;Boussinesq=true,
							mesh=mesh,
							a=a,b=b,
							dealias=dealias,ktol=ktol,
							label=label)

		new(m.label, m.f!, m.mapto, m.mapfro, m.info)
    end
end


text:
"""
    InternalTranslation <: InternalTransition

Saves a translation as described by [`Translation`](@ref) for the current frame.
Is part of the [`Action`](@ref) struct.
"""
mutable struct InternalTranslation <: InternalTransition
    by::Point
end

"""
    InternalRotation <: InternalTransition

Saves a rotation as described by [`Rotation`](@ref) for the current frame.
Is part of the [`Action`](@ref) struct.
"""
mutable struct InternalRotation <: InternalTransition
    angle::Float64
    center::Point
end

"""
    InternalScaling <: InternalTransition

Saves a scaling as described by [`Scaling`](@ref) for the current frame.
Is part of the [`Action`](@ref) struct.
"""
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


text:
###########################################################################################
#Forecast using random forest regression with lagged aggregate variable
#and lagged disaggregate variables.
#This is the inner function to compute the RMSE for each λ
###########################################################################################
@everywhere function rfr_inner_RMSE(J, h, h_lag, J_cv, Y, Y_lag, X_lag, n_est, m, oos_i)
    println("m = $m")
    ## Vector Index for oos forecast goes from J_cv - (h_lag + h) to J - (1 + h_lag + h)
    ŷ_cv = zeros(Float64, size(J_cv - (h_lag + h):J + (oos_i - 1) - (1 + h_lag + h)))
    model = RandomForestRegressor(n_est, max_depth = m, random_state = 0)::PyCall.PyObject
    for (n, vec_î) in enumerate(J_cv - (h_lag + h):J + (oos_i - 1) - (1 + h_lag + h))
        LHS = Y[1:vec_î - h]
        sample_size = length(LHS)
        RHS = [Y_lag[1:sample_size] X_lag[1:sample_size, :]]
        fit!(model, RHS, LHS)
        RHS_crnt = [Y_lag[vec_î]; X_lag[vec_î, :]]'
        ŷ_cv[n] = ScikitLearn.predict(model, RHS_crnt)[1]::Float64
    end
    RMSE = sqrt(mean((Y[J_cv - (h_lag + h):J + (oos_i - 1) - (1 + h_lag + h)] - ŷ_cv) .^ 2))
    return [RMSE, m]
end

###########################################################################################
#Forecast using random forest regression with lagged aggregate variable
#and lagged disaggregate variables.
#This is the outer function to compute the RMSE for each λ
###########################################################################################
function rfr_oos_forc(J, h, h_lag, J_cv, Y, Y_lag, X_lag, m_vec, n_est, oos_i)
    # Step1 : pick the best m
    RMSE_cv = zeros(length(m_vec), 2)
    result = pmap(m -> rfr_inner_RMSE(J, h, h_lag, J_cv, Y, Y_lag, X_lag, n_est, m, oos_i),
                  m_vec)::Array{Array{Float64,1},1}
    for i in 1:length(m_vec)
        RMSE_cv[i,:] = result[i]
    end
    index = findmin(RMSE_cv[:, 1])[2]::Int64
    m = m_vec[index]::Int64
    model = RandomForestRegressor(n_est, max_depth = m, random_state = 0)::PyCall.PyObject
    # Step 2: oos forecasting
    ## Vector Index for oos forecast goes from J - (h_lag + h) to length(Y)
    vec_î = (J - (h_lag + h):length(Y))[oos_i]
    LHS = Y[1:vec_î - h]
    sample_size = length(LHS)
    RHS = [Y_lag[1:sample_size] X_lag[1:sample_size, :]]
    fit!(model, RHS, LHS)
    RHS_crnt = [Y_lag[vec_î]; X_lag[vec_î, :]]'
    ŷ = ScikitLearn.predict(model, RHS_crnt)[1]
    return m, ŷ, RMSE_cv
end

max_stars_repo_name: jungangl/forecast-machine-learning

max_issues_repo_name: jungangl/forecast-machine-learning

max_forks_repo_name: jungangl/forecast-machine-learning