text:
#=
    This file handles the conversion from svg produced by `tex2svg` to Luxor commands.
    It currently misses a lot of possible commands that svg supports.
    Warnings are thrown if this happens such that one can easily grasp whether it
    should have worked or whether there is something missing
=#

"""
    float_attribute(o, name)

Get the attribute `name` of the XMLElement and parse it as a Float64
"""
float_attribute(o::LightXML.XMLElement, name) = parse(Float64, attribute(o, name))

#=
    draw_obj functions have a type which is the type of the svg element.
    i.e <rect calls the `draw_obj(::Val{:rect}, o, defs)` command or `::Val{:path}` for `<path`
=#

"""
    draw_obj(::Val{:rect}, o, defs)

Draw the rectangle defined by the object `o`.
"""
function draw_obj(::Val{:rect}, o, defs)
    width = float_attribute(o, "width")
    height = float_attribute(o, "height")
    x = float_attribute(o, "x")
    y = float_attribute(o, "y")
    # TODO: draw it anti clockwise as it gets scaled to clockwise :/
    move(Point(x + width, y))
    line(Point(x, y))
    line(Point(x, y + height))
    line(Point(x + width, y + height))
    closepath()
end

"""
    draw_obj(::Val{:g}, o, defs)

Draws a group by setting the attributes (like transformations)
and then calls `draw_obj` for all child elements.
"""
function draw_obj(::Val{:g}, o, defs)
    set_attrs(o)

    childs = collect(child_elements(o))
    for child in childs
        sym_name = Symbol(name(child))
        @layer begin
            draw_obj(Val{sym_name}(), child, defs)
        end
    end
end

"""
    draw_obj(::Val{:use}, o, defs)

Calls the command specified in `defs`.
"""
function draw_obj(::Val{:use}, o, defs)
    set_attrs(o)
    id = attribute(o, "href")[2:end]

    if haskey(defs, id)
        def = defs[id]
        sym = Symbol(name(def))
        draw_obj(Val{sym}(), def, defs)
    else
        @warn "There is no definition for $id"
    end
end

"""
    draw_obj(::Val{:path}, o, defs)

Calls the commands specified in the path data.
Currently supports only a subset of possible SVG commands.
"""
function draw_obj(::Val{:path}, o, defs)
    set_attrs(o)
    data = attribute(o, "d")
    counter = 0

    # split without loosing the command
    data_parts = split(data, r"(?=[A-Za-z])")
    # needs to keep track of the current point `c_pt` and the last point `l_pt`
    l_pt = O
    c_pt = O
    circle_pts = []
    for pi in 1:length(data_parts)
        p = data_parts[pi]
        command, args = p[1], p[2:end]
        if command != 'T'
            counter = 0
        end
        # using if else statements instead of dispatching here. Maybe it's faster :D
        if command == 'M'
            c_pt = path_move(parse.(Float64, split(args))...)
            # need to set the last control point as the current point for 'T'
            l_pt = c_pt
        elseif command == 'Q'
            l_pt, c_pt = path_quadratic(c_pt, parse.(Float64, split(args))...)
        elseif command == 'T'
            # the control point is a reflection based on the current and last point
            control_pt = l_pt + 2 * (c_pt - l_pt)
            l_pt, c_pt =
                path_quadratic(c_pt, control_pt..., parse.(Float64, split(args))...)
            push!(circle_pts, Point(parse.(Float64, split(args))...))
            counter += 1
        elseif command == 'L'
            new_pt = Point(parse.(Float64, split(args))...)
            line(new_pt)
            l_pt, c_pt = c_pt, new_pt
        elseif command == 'H'
            new_pt = Point(parse(Float64, args), c_pt.y)
            line(new_pt)
            l_pt, c_pt = c_pt, new_pt
        elseif command == 'V'
            new_pt = Point(c_pt.x, parse(Float64, args))
            line(new_pt)
            l_pt, c_pt = c_pt, new_pt
        elseif command == 'Z'
            closepath()
        else
            @warn "Couldn't parse the svg command: $command"
        end
    end
end

#=
    All kinds of commands for creating a path. The commands need to return the current point
    some of them also a previous point
=#

"""
    path_move(x,y)

Moving to the specified point
"""
function path_move(x, y)
    p = Point(x, y)
    move(p)
    p
end

"""
    path_quadratic(c_pt::Point, x,y, xe, ye)

Drawing a quadratic bezier curve by computing a cubic one as that is supported by Luxor
"""
function path_quadratic(c_pt::Point, x, y, xe, ye)
    e_pt = Point(xe, ye)
    qc = Point(x, y)
    # compute the control points of a cubic bezier curve
    c1 = c_pt + 2 / 3 * (qc - c_pt)
    c2 = e_pt + 2 / 3 * (qc - e_pt)
    curve(c1, c2, e_pt)
    return qc, e_pt
end

#=
    All kinds of setting attribute functions
=#

"""
    set_attrs(o)

Setting the attributes of the object `o` by calling `set_attr` methods.
"""
function set_attrs(o)
    for attribute in attributes(o)
        sym = Symbol(name(attribute))
        set_attr(Val{sym}(), LightXML.value(attribute))
    end
end

"""
    set_attr(::Val{:transform}, transform_strs)

Call the corresponding `set_transform` method i.e `matrix`, `scale` and `translate`
"""
function set_attr(::Val{:transform}, transform_strs)
    if transform_strs !== nothing
        transform_parts = split(transform_strs, r"(?<=[)]) ")
        for transform_str in transform_parts
            m = match(r"(.+)\((.+)\)", transform_str)
            type = Symbol(m.captures[1])
            set_transform(
                Val{type}(),
                parse.(Float64, strip.(split(m.captures[2], r"[, ]")))...,
            )
        end
    end
end

set_attr(::Val{:x}, x) = translate(parse(Float64, x), 0)
set_attr(::Val{:y}, y) = translate(0, parse(Float64, y))
set_attr(::Val{Symbol("stroke-width")}, sw) = setline(parse(Float64, sw))

set_transform(::Val{:translate}, x, y) = translate(x, y)
set_transform(::Val{:scale}, x, y = x) = scale(x, y)

"""
    set_transform(::Val{:matrix}, args...)

Multiply the new matrix with the current matrix and set it.
"""
function set_transform(::Val{:matrix}, args...)
    old = cairotojuliamatrix(getmatrix())
    update = cairotojuliamatrix(collect(args))
    new = old * update
    # only the first two rows are considered
    setmatrix(vec(new[1:2, :]))
end

#=
    General fallbacks
=#
draw_obj(::Union{Val{:title},Val{:defs}}, o, defs) = nothing
# no support for colors at this point
set_attr(t::Union{Val{:stroke},Val{:fill},Val{Symbol("aria-hidden")}}, args...) = nothing

draw_obj(t, o, defs) = @warn "Can't draw $t"
set_transform(t, args...) = @warn "Can't transform $t"
set_attr(::Val{:href}, args...) = nothing
set_attr(::Val{:d}, args...) = nothing
set_attr(::Val{:id}, args...) = nothing
set_attr(t, args...) = @warn "No attr match for $t"

"""
    svgwh(svg)

Return the width and height of svg displayed scaled to half the font size.
"""
function svgwh(svg)
    fsize = get_current_setting().fontsize
    xdoc = parse_string(svg)
    xroot = root(xdoc)
    # remove ex in the end
    ex_width = parse(Float64, attribute(xroot, "width")[1:(end - 2)])
    ex_height = parse(Float64, attribute(xroot, "height")[1:(end - 2)])
    return (fsize * (425 / 1000)) .* (ex_width, ex_height)
end

"""
    pathsvg(svg)

Convert an svg to a path using Luxor. Normally called via the `latex` command.
It handles only a subset of the full power of svg.
"""
function pathsvg(svg)
    fsize = get_current_setting().fontsize
    xdoc = parse_string(svg)
    xroot = root(xdoc)
    def_element = get_elements_by_tagname(xroot, "defs")[1]
    # create a dict for all the definitions
    defs = Dict{String,Any}()
    for def in collect(child_elements(def_element))
        defs[attribute(def, "id")] = def
    end
    x, y, width, height = parse.(Float64, split(attribute(xroot, "viewBox")))
    # remove ex in the end
    ex_width = parse(Float64, attribute(xroot, "width")[1:(end - 2)])
    ex_height = parse(Float64, attribute(xroot, "height")[1:(end - 2)])
    # everything capsulated in a layer
    @layer begin
        # unit ex: half of font size
        # such that we can scale half of a font size (size of a lower letter)
        # with the corresponding height of the svg canvas
        # and the ex_height given in it's description
        scale((fsize * (425 / 1000)) / (height / ex_height))
        translate(-x, -y)

        for child in collect(child_elements(xroot))
            sym_name = Symbol(name(child))
            @layer begin
                draw_obj(Val{sym_name}(), child, defs)
            end
        end
    end
end


text:
using AbstractTrees: AbstractShadowTree

abstract type AbstractJLBoostTree <: AbstractShadowTree end

mutable struct JLBoostTreeModel
	jlt::Vector
	loss # this should be a function with deriv defined
	target::Symbol
end

"""
	trees(jlt::JLBoostTreeModel)

Return the trees from a tree-model
"""
trees(jlt::JLBoostTreeModel) = jlt.jlt

+(v1::JLBoostTreeModel, v2::JLBoostTreeModel) = begin
	v3 = deepcopy(v1)
	v3.jlt = vcat(trees(v1), trees(v2))
	v3
end

mutable struct JLBoostTree <: AbstractJLBoostTree
    weight
	parent::Union{JLBoostTree, Nothing}
    children::Vector{AbstractJLBoostTree}
    splitfeature
    split
    gain
    JLBoostTree(w::T; parent=nothing) where {T <: Number}  = new(w, parent, JLBoostTree[], missing, missing, missing)
end

mutable struct WeightedJLBoostTree <: AbstractJLBoostTree
	tree::JLBoostTree
	eta::Number
end

Base.getproperty(jlt::WeightedJLBoostTree, sym::Symbol) = begin
	if sym == :eta
		return getfield(jlt, :eta)
	elseif sym == :tree
		return getfield(jlt, :tree)
	else
		return getproperty(getfield(jlt, :tree), sym)
	end
end

*(jlt::JLBoostTree, eta::Number) = WeightedJLBoostTree(jlt, eta)

*(eta::Number, jlt::JLBoostTree) = WeightedJLBoostTree(jlt, eta)

*(jlt::WeightedJLBoostTree, eta::Number) = begin
	jlt.eta *= eta
	jlt
end

*(eta::Number, jlt::WeightedJLBoostTree) = begin
	jlt.eta *= eta
	jlt
end

"""
	show(io, jlt, ntabs; splitfeature="")

Show a JLBoostTree
"""
function show(io, jlt::WeightedJLBoostTree, ntabs::I; kwargs...) where {I <: Integer}
	println(io, "eta = $(jlt.eta) (tree weight)")
	show(io, jlt.tree, ntabs;  kwargs...)
end

function show(io, jlt::JLBoostTree, ntabs::I; splitfeature="") where {I <: Integer}
    if ntabs == 0
        tabs = ""
    else ntabs >= 1
        tabs = reduce(*, ["  " for i = 1:ntabs])
    end
    if splitfeature != ""
        println(io, "")
        print(io, "$tabs -- $splitfeature")
        if ismissing(jlt.splitfeature)
            println(io, "")
            println(io, "  $tabs ---- weight = $(jlt.weight)")
        else
            #println(io, "  $tabs ---- $(jlt.splitfeature), weight = $(jlt.weight)")
            #println(io, "  $tabs ---- $(jlt.splitfeature)")
        end
    elseif ismissing(jlt.splitfeature)
        println(io, "$tabs weight = $(jlt.weight)")
    else
        #println("$tabs $(jlt.splitfeature), weight = $(jlt.weight)")
    end

    if length(jlt.children) == 2
        show(io, jlt.children[1], ntabs + 1; splitfeature = "$(jlt.splitfeature) <= $(jlt.split)")
        show(io, jlt.children[2], ntabs + 1; splitfeature = "$(jlt.splitfeature) > $(jlt.split)")
    end
end

"""
	show(io, jlt)

Show a JLBoostTree or a Vector{JLBoostTree}
"""
function show(io::IO, jlt::AbstractJLBoostTree)
    show(io, jlt, 0)
end

function show(io::IO, ::MIME"text/plain", jlt::AbstractVector{T}) where T <: AbstractJLBoostTree
	for (i, tree) in enumerate(jlt)
		println(io, "Tree $i")
    	show(io, tree, 0)
    	println(io, " ")
    end
end

ext: jl


text:
using Test
using PolygonArea
using PolygonArea: PolarHalfPlane, Point, corner_point, Intersection, Reunion
using PolygonArea: signed_distance, distance, outward_normal, angle
using PolygonArea: complement, rotate, translate, distance, exchange_x_and_y
using StaticArrays

@testset "Half-planes" begin

    @testset "creation and equality" begin
        @test HalfPlane(1.0, 0.0, 0.0) isa HalfPlane{Float64}
        @test PolarHalfPlane(0.0, π, center=(1.0, 1.0)) isa HalfPlane{Float64}

        @test HalfPlane(1, 0, 0) isa HalfPlane{Int}
        @test HalfPlane{Float64}(1, 0, 0) isa HalfPlane{Float64}

        @test HalfPlane((0.0, 0.0), (1.0, 0.0)) ≈ HalfPlane(0.0, 1.0, 0.0)

        @test HalfPlane(1.0, 0.0, 0.0) == HalfPlane(1.0, 0.0, 0.0)
        @test HalfPlane(1.0, 0.0, 0.0) != HalfPlane(1.0, 1e-15, 0.0)
        @test HalfPlane(1.0, 0.0, 0.0) ≈ HalfPlane(1.0, 0.0, 0.0)
    end

    @testset "contained points" begin
        @test Point(0, -1)   in HalfPlane(0, 1, 0)
        @test Point(0, -1)   in HalfPlane(0.0, 1.0, 0.0)
        @test Point(0, 0)    in HalfPlane(0, 1, 0)
        @test !(Point(0, 1)  in HalfPlane(0, 1, 0))
        @test Point(1, 0)    in HalfPlane(-1, 0, 0)
        @test !(Point(-1, 0) in HalfPlane(-1, 0, 0))

        @test (-1.0, -1.0) ∈ HalfPlane((0.0, 0.0), (1.0, 0.0), :right)
        @test (1.0, 1.0) ∈ HalfPlane((0.0, 0.0), (1.0, 0.0), :left)

        @test Point(0.5, 0.5) in PolarHalfPlane(-1.0, 0π)
        @test (0.5, 0.5) in PolarHalfPlane(-1.0, 0π)

        @test Point(0.5, 0.5) in PolarHalfPlane(-1.0, π/2)
        @test Point(0.5, 0.5) in PolarHalfPlane(-1.0, π)
        @test Point(0.5, 0.5) in PolarHalfPlane(-1.0, 3π/2)
        @test Point(1.0, 1.0) in PolarHalfPlane(-1.0, 0π)
        @test Point(0.5, 2.0) in PolarHalfPlane(-1.0, 0π)
        @test Point(2.0, 0.5) in PolarHalfPlane(-1.0, π/2)
        @test Point(0.0, 3.0) in PolarHalfPlane(-1.0, π)
    end

    @testset "properties" begin
        @test signed_distance((0.0, 0.0), PolarHalfPlane(-4.0, 2π/3)) ≈ -4.0
        @test signed_distance((0.0, 0.0), PolarHalfPlane(4.0, π/3)) ≈ 4.0

        # Works with Points (SVector) as well as tuples
        @test distance(Point(9.0, 0.0), PolarHalfPlane(0.0, 0.0)) ≈ 9.0
        @test distance((9.0, 0.0), PolarHalfPlane(0.0, 0.0)) ≈ 9.0

        # Scaling of the definition of the plane
        @test distance((1.0, 1.0), HalfPlane(-1.0, 0.0, 0.0)) == distance((1.0, 1.0), HalfPlane(-2.0, 0.0, 0.0))
        @test distance((1.0, 1.0), HalfPlane(-1.0, -1.0, 0.0)) == distance((1.0, 1.0), HalfPlane(-2.0, -2.0, 0.0))

        @test distance((9.0, 0.0), PolarHalfPlane(1.0, 0.0)) ≈ 10.0
        @test distance((9.0, -10.0), PolarHalfPlane(1.0, 0.0)) ≈ 10.0
        @test distance((9.0, 1.0), PolarHalfPlane(0.0, π/2)) ≈ 1.0
        @test distance((1.0, 1.0), PolarHalfPlane(0.0, π/4)) ≈ √2
        @test distance((2.0, 2.0), PolarHalfPlane(-√2, π/4)) ≈ √2
        @test distance((6.0, 1.0), PolarHalfPlane(0.0, π/4, center=Point(5.0, 0.0))) ≈ √2

        @test angle(PolarHalfPlane(0.0, π/8)) ≈ π/8
        @test angle(PolarHalfPlane(0.0, -π/8)) ≈ mod(-π/8, 2π)
        @test angle(PolarHalfPlane(4.0, π/3)) ≈ π/3

        @test outward_normal(HalfPlane(1.0, 0.0, 0.0)) == SVector(1.0, 0.0)
        @test outward_normal(HalfPlane(0.0, 1.0, rand(1)[1])) == SVector(0.0, 1.0)
    end

    @testset "transformations" begin
        # Rotate
        @test outward_normal(rotate(PolarHalfPlane(0π), π)) == outward_normal(PolarHalfPlane(π))
        @test outward_normal(rotate(PolarHalfPlane(0π, center=(1.0, 1.0)), π)) == outward_normal(PolarHalfPlane(π))
        @test outward_normal(rotate(PolarHalfPlane(π/2, center=(1.0, 1.0)), π)) == outward_normal(PolarHalfPlane(3π/2))

        @test (PolarHalfPlane(0.234π) |> rotate(π) |> rotate(-π)) ≈ PolarHalfPlane(0.234π)
        @test (PolarHalfPlane(0.234π) |> rotate(π) |> rotate(π)) ≈ PolarHalfPlane(0.234π)

        @test Point(1.0, 0.0) in rotate(HalfPlane(1.0, 0.0, 0.0) ∪ HalfPlane(1.0, 0.1, 0.0), π)
        @test Point(1.0, 0.0) in rotate(HalfPlane(1.0, 0.0, 0.0) ∩ HalfPlane(1.0, 0.1, 0.0), π)

        # Translate
        @test outward_normal(translate(PolarHalfPlane(0π), SVector{2}(1.0, 1.0))) == outward_normal(PolarHalfPlane(0π))
        @test distance(Point(2.0, 0.0), translate(PolarHalfPlane(0π), SVector{2}(1.0, 1.0))) == 1.0
        @test Point(1.0, 0.0) in translate(HalfPlane(1.0, 0.0, 0.0), SVector(2.0, 0.0))

        @test Point(1.0, 0.0) in translate(HalfPlane(1.0, 0.0, 0.0) ∪ HalfPlane(1.0, 0.1, 0.0), SVector(2.0, 0.0))
        @test Point(1.0, 0.0) in translate(HalfPlane(1.0, 0.0, 0.0) ∩ HalfPlane(1.0, 0.1, 0.0), SVector(2.0, 0.0))

        @test translate(HalfPlane(1.0, 0.0, 0.0), SVector(2.0, 0.0)) == (HalfPlane(1.0, 0.0, 0.0) |> translate(SVector(2.0, 0.0)))

        # Exchange x and y
        @test isapprox(exchange_x_and_y(PolarHalfPlane(0.0, 0.0)), PolarHalfPlane(0.0, π/2), atol=1e-14)
        @test isapprox(exchange_x_and_y(PolarHalfPlane(1.0, 0.0)), PolarHalfPlane(1.0, π/2), atol=1e-14)
        @test isapprox(exchange_x_and_y(PolarHalfPlane(0.0, 0.0, center=Point(1.0, 0.0))), PolarHalfPlane(0.0, π/2, center=Point(0.0, 1.0)), atol=1e-14)
    end

    @testset "unions, intersections and corners" begin
        inferior_hp = HalfPlane(0, 1, 0)
        right_hp = HalfPlane(-1, 0, 0)
        @test Point(1, -1) in intersect(inferior_hp, right_hp)
        @test Point(1, -1) in inferior_hp ∩ right_hp
        @test !(Point(1, 1) in intersect(inferior_hp, right_hp))
        @test Point(1, 1) in complement(inferior_hp ∩ right_hp)
        @test Point(1, 1) in complement(inferior_hp ∪ complement(right_hp))

        @test !(Point(-1, 1) in union(inferior_hp, right_hp))
        @test (Point(1, 1) in union(inferior_hp, right_hp))

        @test corner_point(HalfPlane(0.0, 1.0, 0.0), HalfPlane(-1.0, 0.0, 0.0)) == Point(0.0, 0.0)
        @test corner_point(HalfPlane(1.0, 0.0, 1.0), HalfPlane(0.0, -1.0, 0.5)) == Point(-1.0, 0.5)
        @test corner_point(PolarHalfPlane(0.0, 0.0), HalfPlane(0.0, 1.0, -1.0)) == Point(0.0, 1.0)
        @test corner_point(PolarHalfPlane(-1.0, π/2, center=Point(0.5, 0.5)), HalfPlane(1.0, 0.0, -1.0)) == Point(1.0, 1.5)

        # Union and intersection
        h1 = HalfPlane(1.0, 0.0, 0.0)
        @test ((h1 ∩ h1) ∪ (h1 ∩ h1)) isa Reunion{Intersection{HalfPlane{Float64}}}
        @test ((h1 ∪ h1) ∩ (h1 ∪ h1)) isa Reunion{Intersection{HalfPlane{Float64}}}

        HalfPlane(1.0, 0.0, 0.0) ∪ HalfPlane(0.0, 1.0, 0.0) ∩ (HalfPlane(1.0, 0.0, 0.0) ∪ HalfPlane(0.0, 1.0, 0.0))

        @test Point(0.0, 0.0) in ((h1 ∪ h1) ∩ (h1 ∪ h1) |> complement)
    end

    @testset "conversion" begin
        a = HalfPlane{Int}(1, 1, 1)
        b = HalfPlane{Int}(-1, 1, 1)
        @test convert(Reunion{Intersection{HalfPlane{Int}}}, a ∪ b) isa Reunion{Intersection{HalfPlane{Int}}}
    end

end


text:
function list(config::Config; show_links=false)
	init(config)

    println("Julia Versions:")
    for j in readdir(config.dir.bin)
        if j != "playground"
            file_path = joinpath(config.dir.bin, j)
            if islink(file_path) && show_links
                println("\t$(j) -> $(readlink(file_path))")
            else
                println("\t$(j)")
            end
        end
    end

    println("\nPlaygrounds:")
    for p in readdir(config.dir.store)
        file_path = joinpath(config.dir.store, p)
        if islink(file_path) && show_links
            println("\t$(p) -> $(readlink(file_path))")
        else
            println("\t$(p)")
        end
    end
end


text:
function venkat_limiter(qtilde, globaldata_itm, VL_CONST, phi)
    ds = globaldata_itm.short_distance
    epsi = VL_CONST * ds
    epsi = epsi ^ 3
    del_pos = zero(Float64)
    del_neg = zero(Float64)
    for i in 1:4
        q = globaldata_itm.q[i]
        del_neg = qtilde[i] - q
        if abs(del_neg) <= 1e-5
            phi[i] = 1.0
        elseif abs(del_neg) > 1e-5
            if del_neg > 0.0
                del_pos = globaldata_itm.max_q[i] - q
            elseif del_neg < 0.0
                del_pos = globaldata_itm.min_q[i] - q
            end
            num = (del_pos*del_pos) + (epsi*epsi)
            num = (num*del_neg) + 2 * (del_neg*del_neg*del_pos)

            den = (del_pos*del_pos) + (2 *del_neg*del_neg)
            den = den + (del_neg*del_pos) + (epsi*epsi)
            den = den*del_neg

            temp = num/den
            if temp < 1.0
                phi[i] = temp
            else
                phi[i] = 1.0
            end
        end
    end
    return nothing
end

@inline function qtilde_to_primitive(result, qtilde, gamma)
    beta = -qtilde[4]*0.5
    temp = 0.5/beta
    u1 = qtilde[2]*temp
    u2 = qtilde[3]*temp

    temp1 = qtilde[1] + beta*(u1*u1 + u2*u2)
    temp2 = temp1 - (log(beta)/(gamma-1))
    rho = exp(temp2)
    pr = rho*temp
    result[1] = u1
    result[2] = u2
    result[3] = rho
    result[4] = pr
    return  nothing
end

max_stars_repo_name: TestSubjector/HybridMeshfreeSolver

max_issues_repo_name: TestSubjector/HybridMeshfreeSolver

max_forks_repo_name: TestSubjector/HybridMeshfreeSolver


text:
## Make an animation of a random walk
using SparseArrays
using DelimitedFiles
## Load the graph
ei,ej = copy.(eachcol(Int.(readdlm("undir-graph-14.edges"))))
xy = readdlm("undir-graph-14.xy")'
##
A = sparse(ei,ej,1,14,14) # make the sparse matrix
A = A+A' # make it symmetric (we only store half)
n = size(A,1)
##
d = vec(sum(A,dims=1)) # get the degrees
Di = Diagonal(1.0./d)
##
x = zeros(n)
x[1] = 1
A'*Di*x
##
P = Matrix(A'*Di)
vals,vecs = eigen(P)
vals
vecs[:,1]/sum(vecs[:,1]) # normalize to prob.
##
d/sum(d)


text:
# Rows grouping.
# Maps row contents to the indices of all the equal rows.
# Used by groupby(), join(), nonunique()
struct RowGroupDict{T<:AbstractDataFrame}
    "source data table"
    df::T
    "number of groups"
    ngroups::Int
    "row hashes"
    rhashes::Vector{UInt}
    "hashindex -> index of group-representative row"
    gslots::Vector{Int}
    "group index for each row"
    groups::Vector{Int}
    "permutation of row indices that sorts them by groups"
    rperm::Vector{Int}
    "starts of ranges in rperm for each group"
    starts::Vector{Int}
    "stops of ranges in rperm for each group"
    stops::Vector{Int}
end

# "kernel" functions for hashrows()
# adjust row hashes by the hashes of column elements
function hashrows_col!(h::Vector{UInt},
                       n::Vector{Bool},
                       v::AbstractVector{T}) where T
    @inbounds for i in eachindex(h)
        el = v[i]
        h[i] = hash(el, h[i])
        if T >: Missing && length(n) > 0
            # el isa Missing should be redundant
            # but it gives much more efficient code on Julia 0.6
            n[i] |= (el isa Missing || ismissing(el))
        end
    end
    h
end

# should give the same hash as AbstractVector{T}
function hashrows_col!(h::Vector{UInt},
                       n::Vector{Bool},
                       v::AbstractCategoricalVector{T}) where T
    # TODO is it possible to optimize by hashing the pool values once?
    @inbounds for (i, ref) in enumerate(v.refs)
        h[i] = hash(CategoricalArrays.index(v.pool)[ref], h[i])
    end
    h
end

# should give the same hash as AbstractVector{T}
# enables efficient sequential memory access pattern
function hashrows_col!(h::Vector{UInt},
                       n::Vector{Bool},
                       v::AbstractCategoricalVector{>: Missing})
    # TODO is it possible to optimize by hashing the pool values once?
    @inbounds for (i, ref) in enumerate(v.refs)
        if ref == 0
            h[i] = hash(missing, h[i])
            length(n) > 0 && (n[i] = true)
        else
            h[i] = hash(CategoricalArrays.index(v.pool)[ref], h[i])
        end
    end
    h
end

# Calculate the vector of `df` rows hash values.
function hashrows(df::AbstractDataFrame, skipmissing::Bool)
    rhashes = zeros(UInt, nrow(df))
    missings = fill(false, skipmissing ? nrow(df) : 0)
    for col in columns(df)
        hashrows_col!(rhashes, missings, col)
    end
    return (rhashes, missings)
end

# Helper function for RowGroupDict.
# Returns a tuple:
# 1) the number of row groups in a data table
# 2) vector of row hashes
# 3) slot array for a hash map, non-zero values are
#    the indices of the first row in a group
# Optional group vector is set to the group indices of each row
function row_group_slots(df::AbstractDataFrame,
                         groups::Union{Vector{Int}, Nothing} = nothing,
                         skipmissing::Bool = false)
    rhashes, missings = hashrows(df, skipmissing)
    row_group_slots(ntuple(i -> df[i], ncol(df)), rhashes, missings, groups, skipmissing)
end

function row_group_slots(cols::Tuple{Vararg{AbstractVector}},
                         rhashes::AbstractVector{UInt},
                         missings::AbstractVector{Bool},
                         groups::Union{Vector{Int}, Nothing} = nothing,
                         skipmissing::Bool = false)
    @assert groups === nothing || length(groups) == length(cols[1])
    # inspired by Dict code from base cf. https://github.com/JuliaData/DataFrames.jl/pull/17#discussion_r102481481
    sz = Base._tablesz(length(rhashes))
    @assert sz >= length(rhashes)
    szm1 = sz-1
    gslots = zeros(Int, sz)
    # If missings are to be skipped, they will all go to group 1
    ngroups = skipmissing ? 1 : 0
    @inbounds for i in eachindex(rhashes)
        # find the slot and group index for a row
        slotix = rhashes[i] & szm1 + 1
        # Use 0 for non-missing values to catch bugs if group is not found
        gix = skipmissing && missings[i] ? 1 : 0
        probe = 0
        # Skip rows contaning at least one missing (assigning them to group 0)
        if !skipmissing || !missings[i]
            while true
                g_row = gslots[slotix]
                if g_row == 0 # unoccupied slot, current row starts a new group
                    gslots[slotix] = i
                    gix = ngroups += 1
                    break
                elseif rhashes[i] == rhashes[g_row] # occupied slot, check if miss or hit
                    if isequal_row(cols, i, g_row) # hit
                        gix = groups !== nothing ? groups[g_row] : 0
                    end
                    break
                end
                slotix = slotix & szm1 + 1 # check the next slot
                probe += 1
                @assert probe < sz
            end
        end
        if groups !== nothing
            groups[i] = gix
        end
    end
    return ngroups, rhashes, gslots
end

# Builds RowGroupDict for a given DataFrame.
# Partly uses the code of Wes McKinney's groupsort_indexer in pandas (file: src/groupby.pyx).
function group_rows(df::AbstractDataFrame, skipmissing::Bool = false)
    groups = Vector{Int}(undef, nrow(df))
    ngroups, rhashes, gslots = row_group_slots(df, groups, skipmissing)

    # count elements in each group
    stops = zeros(Int, ngroups)
    @inbounds for g_ix in groups
        stops[g_ix] += 1
    end

    # group start positions in a sorted table
    starts = Vector{Int}(undef, ngroups)
    if !isempty(starts)
        starts[1] = 1
        @inbounds for i in 1:(ngroups-1)
            starts[i+1] = starts[i] + stops[i]
        end
    end

    # define row permutation that sorts them into groups
    rperm = Vector{Int}(undef, length(groups))
    copyto!(stops, starts)
    @inbounds for (i, gix) in enumerate(groups)
        rperm[stops[gix]] = i
        stops[gix] += 1
    end
    stops .-= 1

    # drop group 1 which contains rows with missings in grouping columns
    if skipmissing
        splice!(starts, 1)
        splice!(stops, 1)
        ngroups -= 1
    end

    return RowGroupDict(df, ngroups, rhashes, gslots, groups, rperm, starts, stops)
end

# Find index of a row in gd that matches given row by content, 0 if not found
function findrow(gd::RowGroupDict,
                 df::DataFrame,
                 gd_cols::Tuple{Vararg{AbstractVector}},
                 df_cols::Tuple{Vararg{AbstractVector}},
                 row::Int)
    (gd.df === df) && return row # same table, return itself
    # different tables, content matching required
    rhash = rowhash(df_cols, row)
    szm1 = length(gd.gslots)-1
    slotix = ini_slotix = rhash & szm1 + 1
    while true
        g_row = gd.gslots[slotix]
        if g_row == 0 || # not found
            (rhash == gd.rhashes[g_row] &&
            isequal_row(gd_cols, g_row, df_cols, row)) # found
            return g_row
        end
        slotix = (slotix & szm1) + 1 # miss, try the next slot
        (slotix == ini_slotix) && break
    end
    return 0 # not found
end

# Find indices of rows in 'gd' that match given row by content.
# return empty set if no row matches
function findrows(gd::RowGroupDict,
                  df::DataFrame,
                  gd_cols::Tuple{Vararg{AbstractVector}},
                  df_cols::Tuple{Vararg{AbstractVector}},
                  row::Int)
    g_row = findrow(gd, df, gd_cols, df_cols, row)
    (g_row == 0) && return view(gd.rperm, 0:-1)
    gix = gd.groups[g_row]
    return view(gd.rperm, gd.starts[gix]:gd.stops[gix])
end

function Base.getindex(gd::RowGroupDict, dfr::DataFrameRow)
    g_row = findrow(gd, dfr.df, ntuple(i -> gd.df[i], ncol(gd.df)),
                    ntuple(i -> dfr.df[i], ncol(dfr.df)), dfr.row)
    (g_row == 0) && throw(KeyError(dfr))
    gix = gd.groups[g_row]
    return view(gd.rperm, gd.starts[gix]:gd.stops[gix])
end


text:
# =====================================
# Projectile
#
# system type: affine continuous system
# state dimension: 4
# =====================================

function projectile(; X=Universe(4))
    m = 1.0
    g = 9.81

    # variables: x, y, xdot, ydot
    A = [0 0 1 0;
         0 0 0 1;
         0 0 0 0;
         0 0 0 0.]
    X0 = Singleton([0., 1000.0, 10., 40.0])
    b = [0., 0., 0., -g]
    @ivp(x'= A*x + b, x(0) ∈ X0, x ∈ X)
    tspan = (0.0, 20.0)
    return prob, tspan
end

@taylorize function _projectile!(dx, x, p, t)
    local m = 1.0
    local g = 9.81

    # variables: x, y, xdot, ydot
    dx[1] = x[3]
    dx[2] = x[4]
    dx[3] = zero(x[3])
    dx[4] = -g * one(x[4])
    return dx
end

function projectile_taylor(; X=Universe(4))
    X0 = Singleton([0., 1000.0, 10., 40.0])
    prob = @ivp(x'= _projectile!(x), dim:4, x(0) ∈ X0, x ∈ X)
    tspan = (0.0, 20.0)
    return prob, tspan
end


text:
@testset "Add pseudo + measurements, run mixed SE" begin

    # set measurement path
    msr_path = joinpath(mktempdir(),"temp.csv")
    pseudo_path = joinpath(BASE_DIR, "test/data/extra/measurements/distr_example.csv")

    # 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);

    # solve the power flow
    pf_result = _PMD.solve_mc_pf(data, _PMD.ACPUPowerModel, ipopt_solver)
    pseudo_loads = [3, 11, 7, 16]
    cluster_list = [1, 2, 1, 3]
    _PMDSE.assign_load_pseudo_measurement_info!(data, pseudo_loads, cluster_list; time_step=1, day=1)

    # write measurements based on power flow
    _PMDSE.write_measurements_and_pseudo!(_PMD.ACPUPowerModel, data, pf_result, msr_path, distribution_info=pseudo_path, σ = 0.005)

    # read-in measurement data and set initial values
    _PMDSE.add_measurements!(data, msr_path, actual_meas = true)
    _PMDSE.assign_basic_individual_criteria!(data; chosen_criterion="rwlav")
    data["meas"]["21"]["crit"] = "mle"
    data["meas"]["22"]["crit"] = "mle"
    _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 )

    for (m, meas) in data["meas"]
        if meas["cmp"] == :load && m ∉ ["21", "22", "30", "31"]
            for phase in 1:length(meas["dst"])
                if isa(meas["dst"][phase], _DST.Normal{Float64})
                    @test mean(meas["dst"][phase]) == data["load"]["$(meas["cmp_id"])"][string(meas["var"])][phase]
                end
            end
        end
    end

    data["se_settings"] = Dict()

    # solve the state estimation
    # se_result = _PMDSE.solve_mc_se(data, _PMD.ACPUPowerModel, ipopt_solver)
    # @test se_result["termination_status"] == LOCALLY_SOLVED #CHECK WHY NUMERICAL_ERROR when result seems actually correct
    # @test isapprox( se_result["objective"], 1.41998; atol = 1e-5)
end


text:
module Interface

import Distributions: sample, Sampleable
import Random: GLOBAL_RNG, AbstractRNG

export AbstractSampler,
       AbstractTransition,
       sample_init!,
       sample_end!,
       sample,
       step!

"""
    AbstractSampler

The `AbstractSampler` type is intended to be inherited from when
implementing a custom sampler. Any persistent state information should be
saved in a subtype of `AbstractSampler`.

When defining a new sampler, you should also overload the function
`transition_type`, which tells the `sample` function what type of parameter
it should expect to receive.
"""
abstract type AbstractSampler end

"""
    AbstractTransition

The `AbstractTransition` type describes the results of a single step
of a given sampler. As an example, one implementation of an
`AbstractTransition` might include be a vector of parameters sampled from
a prior distribution.

Transition types should store a single draw from any sampler, since the
interface will sample `N` times, and store the results of each step in an
array of type `Array{Transition<:AbstractTransition, 1}`. If you were
using a sampler that returned a `NamedTuple` after each step, your
implementation might look like:

```
struct MyTransition <: AbstractTransition
    draw :: NamedTuple
end
```
"""
abstract type AbstractTransition end

"""
    sample(
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...)

    sample(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...)

A generic interface for samplers.
"""
function sample(
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    return sample(GLOBAL_RNG, ℓ, s, N)
end

"""
    sample(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

`sample` returns an `MCMCChains.Chains` object containing `N` samples from a given model and
sampler.
"""
function sample(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    # Perform any necessary setup.
    sample_init!(rng, ℓ, s, N; kwargs...)

    # Preallocate the TransitionType vector.
    ts = transitions_init(rng, ℓ, s, N; kwargs...)

    # Step through the sampler.
    for i=1:N
        if i == 1
            ts[i] = step!(rng, ℓ, s, N; kwargs...)
        else
            ts[i] = step!(rng, ℓ, s, N, ts[i-1]; kwargs...)
        end
    end

    # Wrap up the sampler, if necessary.
    sample_end!(rng, ℓ, s, N, ts; kwargs...)

    return Chains(ts)
end

"""
    sample_init!(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

Performs whatever initial setup is required for your sampler.
"""
function sample_init!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    # Do nothing.
    @warn "No sample_init! function has been implemented for objects
           of types $(typeof(ℓ)) and $(typeof(s))"
end

"""
    sample_end!(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer,
        ts::Vector{TransitionType};
        kwargs...
    )

Performs whatever finalizing the sampler requires.
"""
function sample_end!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    ts::Vector{TransitionType};
    kwargs...
) where {
    ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    TransitionType<:AbstractTransition
}
    # Do nothing.
    @warn "No sample_end! function has been implemented for objects
           of types $(typeof(ℓ)) and $(typeof(s))"
end

"""
    step!(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

Returns a single `AbstractTransition` drawn using the model and sampler type.
This is a unique step function called the first time a sampler runs.
"""
function step!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    # Do nothing.
    @warn "No step! function has been implemented for objects
           of types $(typeof(ℓ)) and $(typeof(s))"
end

"""
    step!(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer,
        t::TransitionType;
        kwargs...
    )

Returns a single `AbstractTransition` drawn using the model and sampler type.
"""
function step!(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer,
    t::TransitionType;
    kwargs...
) where {ModelType<:Sampleable,
    SamplerType<:AbstractSampler,
    TransitionType<:AbstractTransition
}
    # Do nothing.
    @warn "No step! function has been implemented for objects
           of types $(typeof(ℓ)) and $(typeof(s))"
end

"""
    transitions_init(
        rng::AbstractRNG,
        ℓ::ModelType,
        s::SamplerType,
        N::Integer;
        kwargs...
    )

Generates a vector of `AbstractTransition` types of length `N`.
"""
function transitions_init(
    rng::AbstractRNG,
    ℓ::ModelType,
    s::SamplerType,
    N::Integer;
    kwargs...
) where {ModelType<:Sampleable, SamplerType<:AbstractSampler}
    @warn "No transitions_init function has been implemented
           for objects of types $(typeof(ℓ)) and $(typeof(s))"
    return Vector(undef, N)
end

end # module Interface


text:
# Symmetric triangular distribution
immutable SymTriangularDist <: ContinuousUnivariateDistribution
    location::Float64
    scale::Float64
    function SymTriangularDist(l::Real, s::Real)
        s > zero(s) || error("scale must be positive")
        new(float64(l), float64(s))
    end
end

SymTriangularDist(location::Real) = SymTriangularDist(location, 1.0)
SymTriangularDist() = SymTriangularDist(0.0, 1.0)

@distr_support SymTriangularDist d.location - d.scale d.location + d.scale

## Properties
mean(d::SymTriangularDist) = d.location
median(d::SymTriangularDist) = d.location
mode(d::SymTriangularDist) = d.location

var(d::SymTriangularDist) = d.scale^2 / 6.0
skewness(d::SymTriangularDist) = 0.0
kurtosis(d::SymTriangularDist) = -0.6

entropy(d::SymTriangularDist) = 0.5 + log(d.scale)

## Functions
function pdf(d::SymTriangularDist, x::Float64)
    if insupport(d, x)
        (1.0 - abs(x-d.location)/d.scale) / d.scale
    else
        0.0
    end
end

function logpdf(d::SymTriangularDist, x::Float64)
    if insupport(d, x)
        log1p(-abs(x-d.location)/d.scale) -log(d.scale)
    else
        -Inf
    end
end

function cdf(d::SymTriangularDist, x::Float64)
    u = (x - d.location) / d.scale
    if u <= -1.0
        0.0
    elseif u <= 0.0
        0.5*(1.0 + u)^2
    elseif u <= 1.0
        1.0 - 0.5*(1.0 - u)^2
    else
        1.0
    end
end
function ccdf(d::SymTriangularDist, x::Float64)
    u = (x - d.location) / d.scale
    if u <= -1.0
        1.0
    elseif u <= 0.0
        1.0 - 0.5*(1.0 + u)^2
    elseif u <= 1.0
        0.5*(1.0 - u)^2
    else
        0.0
    end
end

function logcdf(d::SymTriangularDist, x::Float64)
    u = (x - d.location) / d.scale
    if u <= -1.0
        -Inf
    elseif u <= 0.0
        loghalf + 2.0*log1p(u)
    elseif u <= 1.0
        log1p(-0.5*(1.0 - u)^2)
    else
        0.0
    end
end
function logccdf(d::SymTriangularDist, x::Float64)
    u = (x - d.location) / d.scale
    if u <= -1.0
        0.0
    elseif u <= 0.0
        log1p(- 0.5*(1.0 + u)^2)
    elseif u <= 1.0
        loghalf + 2.0*log1p(-u)
    else
        -Inf
    end
end

function quantile(d::SymTriangularDist, p::Float64)
    @checkquantile p begin
        if p < 0.5
            d.location - d.scale*(1.0 - sqrt(2.0*p))
        else
            d.location + d.scale*(1.0 - sqrt(2.0*(1.0-p)))
        end
    end
end

function cquantile(d::SymTriangularDist, p::Float64)
    @checkquantile p begin
        if p > 0.5
            d.location - d.scale*(1.0 - sqrt(2.0*(1.0-p)))
        else
            d.location + d.scale*(1.0 - sqrt(2.0*p))
        end
    end
end

function invlogcdf(d::SymTriangularDist, lp::Float64)
    @checkinvlogcdf lp begin
        if lp < loghalf
            d.location + d.scale*expm1(0.5*(lp-loghalf))
        else
            d.location + d.scale*(1.0 - sqrt(-2.0*expm1(lp)))
        end
    end
end

function invlogccdf(d::SymTriangularDist, lp::Float64)
    @checkinvlogcdf lp begin
        if lp > loghalf
            d.location - d.scale*(1.0 - sqrt(-2.0*expm1(lp)))
        else
            d.location - d.scale*expm1(0.5*(lp-loghalf))
        end
    end
end

function mgf(d::SymTriangularDist, t::Real)
    a = d.scale*t
    a == zero(a) && return one(a)
    4.0*exp(d.location*t)*(sinh(0.5*a)/a)^2
end
function cf(d::SymTriangularDist, t::Real)
    a = d.scale*t
    a == zero(a) && return complex(one(a))
    4.0*exp(im*d.location*t)*(sin(0.5*a)/a)^2
end

## Sampling
function rand(d::SymTriangularDist)
    ξ1, ξ2 = rand(), rand()
    return d.location + (ξ1 - ξ2) * d.scale
end


text:
function get_ema(x; ema=[8,13,21,55])
    output = []
    for e in ema
        push!(output, Indicators.ema(float.(x), n=e))
    end
    return output
end


text:
using Poly
using Hwloc
using ISL
import ISL.API

"""
get the best tiling of arrays of n dimensions, on band_node, using l1 cache sizes only
assumes float64 arrays (worst-case sizing)
trys to find a tile size tile_1 x tile_2 x .... x tile_t where tile_1 x ... x tile_n <= l1cachesize (the array being tiled will fit in the l1 cache) and tile_i <= l1linesize for all i (all tiled dimensions will fit in a cache line)
finds square tiles

returns: isl_multi_val representation of tiling sizes
"""
function tiling_sizes(n, band::Ptr{ISL.API.isl_schedule_node}, context::Ptr{ISL.API.isl_ctx})::Ptr{ISL.API.isl_multi_val}
    size = Hwloc.l1cache_sizes()[1]/sizeof(Float64)
    linesize = Hwloc.l1cache_linesizes()[1] # note: don't divide by sizeof(float) to utilize l2 cache also
    tiles = Int[]
    tile_size = Int(floor(min(linesize, size^(1/n))))
    band_space = ISL.API.isl_schedule_node_band_get_space(band)
    dim = ISL.API.isl_space_dim(band_space, ISL.API.isl_dim_type(3))
    tiles = ISL.API.isl_val_list_alloc(context, dim)
    for i=1:dim
        ISL.API.isl_val_list_add(tiles, ISL.API.isl_val_int_from_si(context, tile_size))
    end
    tiles = ISL.API.isl_multi_val_from_val_list(band_space, tiles)
    return tiles
end

"""
shift the band to start at 0 so that full tiling occurs
"""
function shift_band_zero(band::Ptr{ISL.API.isl_schedule_node}, context::Ptr{ISL.API.isl_ctx}) ::Tuple{Ptr{ISL.API.isl_multi_union_pw_aff}, Ptr{ISL.API.isl_multi_union_pw_aff}}
    # band domain
    domain = ISL.API.isl_schedule_node_get_domain(band)
    # mupa partial schedule
    partial_schedule = ISL.API.isl_schedule_node_band_get_partial_schedule(band)
    mupa = ISL.API.isl_multi_union_pw_aff_intersect_domain(ISL.API.isl_multi_union_pw_aff_copy(partial_schedule), ISL.API.isl_union_set_copy(domain))
    # get mins of range
    n = ISL.API.isl_multi_union_pw_aff_size(mupa)
    multi_val = ISL.API.isl_multi_union_pw_aff_min_multi_val(mupa)
    for i=0:n-1
        v = ISL.API.isl_multi_val_get_val(multi_val, i)
        if ISL.API.isl_val_is_neginfty(v) == ISL.API.isl_bool_true
            # if -inf, set to 1 by default
            new_v = ISL.API.isl_val_int_from_si(context, 1)
            multi_val = ISL.API.isl_multi_val_set_val(multi_val, i, new_v)
        end
        ISL.API.isl_val_free(v)
    end
    # get shift mupa
    shift = ISL.API.isl_multi_union_pw_aff_multi_val_on_domain(domain, multi_val)
    shift_neg = ISL.API.isl_multi_union_pw_aff_neg(ISL.API.isl_multi_union_pw_aff_copy(shift))
    # shift partial schedule
    partial_schedule = ISL.API.isl_multi_union_pw_aff_add(partial_schedule, shift_neg)
    # keep shift (mupa) to shift back up later
    return partial_schedule, shift
end

"""
tile a partial schedule
taken from ISL schedule.c private function isl_multi_union_pw_aff_tile
"""
function tile_partial_schedule(sched::Ptr{ISL.API.isl_multi_union_pw_aff}, sizes::Ptr{ISL.API.isl_multi_val})::Ptr{ISL.API.isl_multi_union_pw_aff}
    ctx = ISL.API.isl_multi_val_get_ctx(sizes)
    scale = ISL.API.isl_options_get_tile_scale_tile_loops(ctx)
    n = ISL.API.isl_multi_union_pw_aff_size(sched)
    for i=0:n-1
        upa = ISL.API.isl_multi_union_pw_aff_get_union_pw_aff(sched, i)
        v = ISL.API.isl_multi_val_get_val(sizes, i)
        upa = ISL.API.isl_union_pw_aff_scale_down_val(upa, ISL.API.isl_val_copy(v))
		upa = ISL.API.isl_union_pw_aff_floor(upa)
		if scale == ISL.API.isl_bool_true
			upa = ISL.API.isl_union_pw_aff_scale_val(upa, ISL.API.isl_val_copy(v))
        end
        ISL.API.isl_val_free(v)
        sched = ISL.API.isl_multi_union_pw_aff_set_union_pw_aff(sched, i, upa)
    end
    return sched
end

"""
tile a band node with tile dimension n (max dimension of tile)
"""
function tile_band(n, band::Ptr{ISL.API.isl_schedule_node}, context::Ptr{ISL.API.isl_ctx})::Ptr{ISL.API.isl_schedule}
    # shift to zero
    partial_schedule, shift = shift_band_zero(band, context)
    # tile
    multi_val = tiling_sizes(n, band, context)
    partial_schedule = tile_partial_schedule(partial_schedule, multi_val)
    # shift back to original dims
    partial_schedule = ISL.API.isl_multi_union_pw_aff_add(partial_schedule, shift)
    # insert tiled schedule below band
    band = ISL.API.isl_schedule_node_insert_partial_schedule(band, partial_schedule)
    # get new schedule
    schedule = ISL.API.isl_schedule_node_get_schedule(band)
    return schedule
end

"""
get the max tile dimension
"""
function get_tile_dim(kernel::LoopKernel)::Int
    max_dim = 0
    for instruction in kernel.instructions
        dim = get_tile_dim(instruction.body)
        max_dim = max(max_dim, dim)
    end
    return max_dim
end

function get_tile_dim(ex::Expr)::Int
    max_dim = 0
    if ex.head == :ref
        # care about dimension of arrays, so references
        return length(ex.args[2:end]) # number of dimensions in ref
    else
        for arg in ex.args
            max_dim = max(get_tile_dim(arg), max_dim)
        end
    end
    return max_dim
end

function get_tile_dim(any)::Int
    return 0
end

"""
tile a schedule

inputs:
    kernel - LoopKernel to tile
    schedule - ISL schedule of kernel
    context - ISL context
    tile_dim - dimensions to tile (0 means don't tile)

returns:
    schedule - ISL schedule that is tiled
"""
function tile_schedule(kernel::LoopKernel, schedule::Ptr{ISL.API.isl_schedule}, context::Ptr{ISL.API.isl_ctx}, tile_dim::Int)::Ptr{ISL.API.isl_schedule}
    tile_schedule = schedule
    n = get_tile_dim(kernel)
    root = ISL.API.isl_schedule_get_root(schedule)
    node = root
    next_nodes = Set()
    # search for the outermost band node to tile
    while ISL.API.isl_schedule_node_n_children(node) > 0
        # go through all children of current node
        tiled = false
        for i=0:ISL.API.isl_schedule_node_n_children(node)-1
            band = ISL.API.isl_schedule_node_get_child(node, i)
            if ISL.API.isl_schedule_node_get_type(band) == ISL.API.isl_schedule_node_type(0) # band node
                # in a band node, tile band with dimension n
                tile_schedule = tile_band(n, band, context)
                # need to add child of band, since split band
                push!(next_nodes, ISL.API.isl_schedule_node_get_child(band, 0))
                tiled = true
                break
            else
                # add to next nodes
                push!(next_nodes, band)
            end
        end

        if tiled
            break
        end

        if length(next_nodes) == 0
            break # done, no children left
        end
        node = pop!(next_nodes) # get next node
    end

    return tile_schedule
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENCE in the project root.
# ------------------------------------------------------------------

"""
    𝓒(x)

Return closure of `x`.
"""
𝓒(x) = x ./ sum(x)

"""
    Composition{D}(parts)

A D-part composition as defined by Aitchison 1986.
"""
struct Composition{D}
  parts::SVector{D,Float64}
end

Composition(parts) = Composition{length(parts)}(parts)
Composition(parts...) = Composition(parts)

+(c₁::Composition, c₂::Composition) = Composition(𝓒(c₁.parts .* c₂.parts))

-(c::Composition) = Composition(𝓒(1. ./ c.parts))

-(c₁::Composition, c₂::Composition) = c₁ + -c₂

*(λ::Real, c::Composition) = Composition(𝓒(c.parts.^λ))

==(c₁::Composition, c₂::Composition) = 𝓒(c₁.parts) ≈ 𝓒(c₂.parts)

"""
    inner(c₁, c₂)

Inner product between compositions `c₁` and `c₂`.
"""
function inner(c₁::Composition{D}, c₂::Composition{D}) where {D}
  x = c₁.parts; y = c₂.parts
  sum(log(x[i]/x[j])*log(y[i]/y[j]) for j=1:D for i=j+1:D) / D
end

"""
    norm(c)

Aitchison norm of composition `c`.
"""
norm(c::Composition) = √inner(c,c)

"""
    distance(c₁, c₂)

Aitchison distance between compositions `c₁` and `c₂`.
"""
distance(c₁::Composition, c₂::Composition) = norm(c₁ - c₂)

# ------------
# IO methods
# ------------
function Base.show(io::IO, c::Composition{D}) where {D}
  print(io, c.parts)
end

function Base.show(io::IO, ::MIME"text/plain", c::Composition{D}) where {D}
  print(barplot(["part $i" for i in 1:D], c.parts, title="$D-part composition"))
end

function Base.show(io::IO, ::MIME"text/html", c::Composition{D}) where {D}
  print(barplot(["part $i" for i in 1:D], c.parts, title="$D-part composition"))
end


text:
using OffsetArrays

@testset "Associated Legendre Functions" begin
    nmax = 1_000
    min_length = (nmax * (nmax + 3)) ÷ 2 + 1

    for T in [BigFloat, Float64, Float32]
        ϵ = eps(T)
        P̄ = zeros(T, min_length)
        recursion_coefficients = ALFRecursionCoefficients(nmax, T)

        θs = collect(zero(T):one(T)/10:T(π)/2)

        for (i, θ) in enumerate(θs)
            expiβ = exp(θ * 1im)
            ALFcompute!(P̄, expiβ, nmax, recursion_coefficients)
            offset = -1
            for n in 0:nmax
                offset -= n
                P̄ₙ = OffsetVector(P̄, offset)
                total = sum(P̄ₙ[m]^2 for m in 0:n)
                @test total≈T(2n+1) rtol=2ϵ*n
            end
        end
    end
end


text:
@testset "bearing" begin

    start = Position([-75, 45])
    stop = Position([20, 60])

    @test bearing(start, stop, false) ≈ 37.75 atol=0.01
    @test bearing(start, stop, true) ≈ 120.01 atol=0.005
    @test rhumb_bearing(start, stop, false) ≈ 75.28 atol=0.001
    @test rhumb_bearing(start, stop, true) ≈ -104.71 atol=0.01
    @test bearing_to_azimuth(40.) == 40.
    @test bearing_to_azimuth(-105.) == 255.
    @test bearing_to_azimuth(410.) == 50.
end


text:
@testset "Kagome" begin

    @test ndims(Kagome((10,5))) == 2

    @testset "Constructors" begin
        # test default BC
        @test Kagome((10, 5)) == Kagome((10, 5), Periodic())

        for B1 in [Periodic, Open, Helical]
            @test (
                Kagome((10, 5), B1())
                == Kagome((10, 5), (B1(), B1()))
                == Kagome((10, 5), B1())
                == Kagome((10, 5), (B1(), B1()))
            )
            if B1 != Helical
                for B2 in [Periodic, Open]
                    @test (
                        Kagome((10, 5), (B1(), B2()))
                        == Kagome((10, 5), (B1(), B2()))
                    )
                end
                @test_throws ArgumentError Kagome((10, 5), (B1(), Helical()))
            end
        end
    end

    @testset "Translation Vectors" begin
        l = Kagome((10, 10))

        u1 = Set([
            Coordinate(0, 0, 2),
            Coordinate(0, 0, 3),
            Coordinate(-1, 0, 2),
            Coordinate(0, -1, 3)
        ])

        u2 = Set([
            Coordinate(0, 0, 1),
            Coordinate(0, 0, 3),
            Coordinate(1, 0, 1),
            Coordinate(1, -1, 3)
        ])

        u3 = Set([
            Coordinate(0, 0, 1),
            Coordinate(0, 0, 2),
            Coordinate(0, 1, 1),
            Coordinate(-1, 1, 2)
        ])

        @test Set(translation_vectors(l, 1, 1)) == u1
        @test Set(translation_vectors(l, 1, 2)) == u2
        @test Set(translation_vectors(l, 1, 3)) == u3
    end

    @testset "Nearest Neighbors" begin
        l = Kagome((10, 10), Periodic())
        @test Set(neighbors(l, Coordinate(10, 10, 1))) == Set([Coordinate(10, 10, 2), Coordinate(10, 10, 3), Coordinate(9, 10, 2), Coordinate(10, 9, 3)])
        @test Set(neighbors(l, Coordinate(5, 10, 1))) == Set([Coordinate(5, 10, 2), Coordinate(5, 10, 3), Coordinate(4, 10, 2), Coordinate(5, 9, 3)])
        @test Set(neighbors(l, Coordinate(10, 5, 1))) == Set([Coordinate(10, 5, 2), Coordinate(10, 5, 3), Coordinate(9, 5, 2), Coordinate(10, 4, 3)])
        @test Set(neighbors(l, Coordinate(5, 5, 1))) == Set([Coordinate(5, 5, 2), Coordinate(5, 5, 3), Coordinate(4, 5, 2), Coordinate(5, 4, 3)])

        @test Set(neighbors(l, Coordinate(10, 10, 2))) == Set([Coordinate(10, 10, 1), Coordinate(10, 10, 3), Coordinate(1, 10, 1), Coordinate(1, 9, 3)])
        @test Set(neighbors(l, Coordinate(5, 10, 2))) == Set([Coordinate(5, 10, 1), Coordinate(5, 10, 3), Coordinate(6, 10, 1), Coordinate(6, 9, 3)])
        @test Set(neighbors(l, Coordinate(10, 5, 2))) == Set([Coordinate(10, 5, 1), Coordinate(10, 5, 3), Coordinate(1, 5, 1), Coordinate(1, 4, 3)])
        @test Set(neighbors(l, Coordinate(5, 5, 2))) == Set([Coordinate(5, 5, 1), Coordinate(5, 5, 3), Coordinate(6, 5, 1), Coordinate(6, 4, 3)])

        @test Set(neighbors(l, Coordinate(10, 10, 3))) == Set([Coordinate(10, 10, 1), Coordinate(10, 10, 2), Coordinate(10, 1, 1), Coordinate(9, 1, 2)])
        @test Set(neighbors(l, Coordinate(5, 10, 3))) == Set([Coordinate(5, 10, 1), Coordinate(5, 10, 2), Coordinate(5, 1, 1), Coordinate(4, 1, 2)])
        @test Set(neighbors(l, Coordinate(10, 5, 3))) == Set([Coordinate(10, 5, 1), Coordinate(10, 5, 2), Coordinate(10, 6, 1), Coordinate(9, 6, 2)])
        @test Set(neighbors(l, Coordinate(5, 5, 3))) == Set([Coordinate(5, 5, 1), Coordinate(5, 5, 2), Coordinate(5, 6, 1), Coordinate(4, 6, 2)])

        l = Kagome((10, 10), Open())
        @test Set(neighbors(l, Coordinate(10, 10, 1))) == Set([Coordinate(10, 10, 2), Coordinate(10, 10, 3), Coordinate(9, 10, 2), Coordinate(10, 9, 3)])
        @test Set(neighbors(l, Coordinate(5, 10, 1))) == Set([Coordinate(5, 10, 2), Coordinate(5, 10, 3), Coordinate(4, 10, 2), Coordinate(5, 9, 3)])
        @test Set(neighbors(l, Coordinate(10, 5, 1))) == Set([Coordinate(10, 5, 2), Coordinate(10, 5, 3), Coordinate(9, 5, 2), Coordinate(10, 4, 3)])
        @test Set(neighbors(l, Coordinate(5, 5, 1))) == Set([Coordinate(5, 5, 2), Coordinate(5, 5, 3), Coordinate(4, 5, 2), Coordinate(5, 4, 3)])

        @test Set(neighbors(l, Coordinate(10, 10, 2))) == Set([Coordinate(10, 10, 1), Coordinate(10, 10, 3)])
        @test Set(neighbors(l, Coordinate(5, 10, 2))) == Set([Coordinate(5, 10, 1), Coordinate(5, 10, 3), Coordinate(6, 10, 1), Coordinate(6, 9, 3)])
        @test Set(neighbors(l, Coordinate(10, 5, 2))) == Set([Coordinate(10, 5, 1), Coordinate(10, 5, 3)])
        @test Set(neighbors(l, Coordinate(5, 5, 2))) == Set([Coordinate(5, 5, 1), Coordinate(5, 5, 3), Coordinate(6, 5, 1), Coordinate(6, 4, 3)])

        @test Set(neighbors(l, Coordinate(10, 10, 3))) == Set([Coordinate(10, 10, 1), Coordinate(10, 10, 2)])
        @test Set(neighbors(l, Coordinate(5, 10, 3))) == Set([Coordinate(5, 10, 1), Coordinate(5, 10, 2)])
        @test Set(neighbors(l, Coordinate(10, 5, 3))) == Set([Coordinate(10, 5, 1), Coordinate(10, 5, 2), Coordinate(10, 6, 1), Coordinate(9, 6, 2)])
        @test Set(neighbors(l, Coordinate(5, 5, 3))) == Set([Coordinate(5, 5, 1), Coordinate(5, 5, 2), Coordinate(5, 6, 1), Coordinate(4, 6, 2)])

    end

end


text:
# imlements some behavior for length-2 StaticVectors and 2x2
# StaticMatrices of Fields or FieldOps. these come up when you do
# gradients, e.g. ∇*ϕ is a vector of Fields.

# where elements are Fields or FieldOps
const FieldOrOpVector{F<:FieldOrOp}    = StaticVector{2,F}
const FieldOrOpMatrix{F<:FieldOrOp}    = StaticMatrix{2,2,F}
const FieldOrOpRowVector{F<:FieldOrOp} = Adjoint{<:Any,<:FieldOrOpVector{F}}
const FieldOrOpArray{F<:FieldOrOp}     = Union{FieldOrOpVector{F}, FieldOrOpMatrix{F}, FieldOrOpRowVector{F}}
# where elements are just Fields
const FieldVector{F<:Field}    = SVector{2,F}
const FieldMatrix{F<:Field}    = SMatrix{2,2,F,4}
const FieldRowVector{F<:Field} = FieldOrOpRowVector{F}
const FieldArray{F<:Field}     = FieldOrOpArray{F}

# non-broadcasted algebra
for f in (:/, :\, :*)
    if f != :/
        @eval ($f)(A::Field, B::FieldOrOpArray) = broadcast($f, A, B)
    end
    if f != :\
        @eval ($f)(A::FieldOrOpArray, B::Field) = broadcast($f, A, B)
    end
end

# this makes Vector{Diagonal}' * Vector{Field} work right
dot(D::DiagOp{<:Field{B}}, f::Field) where {B} = conj(D.diag) .* B(f)

# needed since v .* f is not type stable
*(v::FieldOrOpVector, f::Field) = @SVector[v[1]*f, v[2]*f]
*(v::FieldOrOpVector, w::FieldOrOpRowVector) = @SMatrix[v[1]*w[1] v[1]*w[2]; v[2]*w[1] v[2]*w[2]]
*(f::Field, v::FieldOrOpVector) = @SVector[f*v[1], f*v[2]]
*(f::Field, v::FieldOrOpRowVector) = @SVector[(f*v[1])', (f*v[2])']'
*(v::FieldOrOpRowVector, w::FieldOrOpVector) = v[1]*w[1] + v[2]*w[2]

# ffs how is something this simple broken in StaticArrays...
adjoint(L::FieldOrOpMatrix) = @SMatrix[L[1,1]' L[2,1]'; L[1,2]' L[2,2]']

# eventually replace having to do this by hand with Cassette-based solution
mul!(f::Field, v::FieldOrOpRowVector{<:Diagonal}, w::FieldVector) =
    ((@. f = v[1].diag * w[1] + v[2].diag * w[2]); f)
mul!(f::Field, v::FieldOrOpRowVector{<:Diagonal}, x::Diagonal, w::FieldVector) =
    ((@. f = x.diag * (v[1].diag * w[1] + v[2].diag * w[2])); f)
mul!(v::FieldOrOpVector{<:Diagonal}, M::FieldOrOpMatrix{<:Diagonal}, w::FieldOrOpVector{<:Diagonal}) =
    ((@. v[1].diag = M[1,1].diag*w[1].diag + M[1,2].diag*w[2].diag); (@. v[2].diag = M[2,1].diag*w[1].diag + M[2,2].diag*w[2].diag); v)
mul!(v::FieldVector, M::FieldOrOpMatrix{<:Diagonal}, w::FieldVector) =
    ((@. v[1] = M[1,1].diag*w[1] + M[1,2].diag*w[2]); (@. v[2] = M[2,1].diag*w[1] + M[2,2].diag*w[2]); v)
mul!(v::FieldVector, w::FieldOrOpVector{<:Diagonal}, f::Field) =
    ((@. v[1] = w[1].diag * f); (@. v[2] = w[2].diag * f); v)
mul!(v::FieldVector, x::Diagonal, w::FieldOrOpVector{<:Diagonal}, f::Field) =
    ((@. v[1] = x.diag * w[1].diag * f); (@. v[2] = x.diag * w[2].diag * f); v)
# only thing needed for SpinAdjoints
mul!(v::FieldVector, f::SpinAdjoint, w::FieldVector) = (mul!(v[1], f, w[1]); mul!(v[2], f, w[2]); v)

function pinv!(dst::FieldOrOpMatrix{<:Diagonal}, src::FieldOrOpMatrix{<:Diagonal})
    a,b,c,d = src
    det⁻¹ = pinv(@. a*d-b*c)
    @. dst[1,1] =  det⁻¹ * d
    @. dst[1,2] = -det⁻¹ * b
    @. dst[2,1] = -det⁻¹ * c
    @. dst[2,2] =  det⁻¹ * a
    dst
end

promote_rule(::Type{F}, ::Type{<:Scalar}) where {F<:Field} = F
arithmetic_closure(::F) where {F<:Field} = F


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:light
#     text_representation:
#       extension: .jl
#       format_name: light
#       format_version: '1.4'
#       jupytext_version: 1.2.1
#   kernelspec:
#     display_name: Julia 1.1.1
#     language: julia
#     name: julia-1.1
# ---

# # DifferentialEquations.jl Workshop Exercise 1
#
# - Exercise 1 takes the user through solving a stiff ordinary differential equation
#   and using the ModelingToolkit.jl to automatically convert the function to a
#   symbolic form to derive the analytical Jacobian to speed up the solver. The
#   same biological system is then solved with stochasticity, utilizing
#   EnsembleProblems to understand 95% bounds on the solution. Finally,
#   probabilistic programming is employed to perform Bayesian parameter estimation
#   of the parameters against data.

# # Problem 1: Investigating Sources of Randomness and Uncertainty in a Stiff Biological System (B)
#
# In this problem we will walk through the basics of simulating models with
# DifferentialEquations.jl. Let's take the
# [Oregonator model of the Belousov-Zhabotinskii chemical reaction system](https://www.radford.edu/~thompson/vodef90web/problems/demosnodislin/Demos_Pitagora/DemoOrego/demoorego.pdf).
# This system describes a classical example in non-equilibrium thermodynmics
# and is a well-known natural chemical oscillator.
#
# ## Part 1: Simulating the Oregonator ODE model
#
# When modeling, usually one starts off by investigating the deterministic model.
# The deterministic ODE formulation of the Oregonator is
# given by the equations
#
# $$\begin{align}
# \frac{dx}{dt} &= s(y-xy + x - qx^2)\\
# \frac{dy}{dt} &= (-y - xy + z)/s\\
# \frac{dz}{dt} &= w(x - z)\end{align}$$

# +
using ModelingToolkit, DifferentialEquations, Plots

@parameters t s w q
@variables x(t) y(t) z(t)
@derivatives D'~t
# -

eqs = [D(x) ~ s*(y-x*y+x-q*x^2),
       D(y) ~ (-y - x*y + z)/s,
       D(z) ~ w*(x - z)]

de = ODESystem(eqs)

generate_function(de, [x,y,z], [s, w, q])

f = ODEFunction(de, [x,y,z], [s, w, q])

# with parameter values $s=77.27$, $w=0.161$, and $q=8.375 \times 10^{-6}$, and
# initial conditions $x(0)=1$, $y(0)=2$, and $z(0)=3$. Use
# [the tutorial on solving ODEs](http://docs.juliadiffeq.org/latest/tutorials/ode_example.html)
# to solve this differential equation on the
# timespan of $t\in[0,360]$ with the default ODE solver. To investigate the result,
# plot the solution of all components over time, and plot the phase space plot of
# the solution (hint: use `vars=(1,2,3)`). What shape is being drawn in phase space?

u0 = [1.0,2.0,3.0]
tspan = (0.0,360.0)
p = [77.27,0.161,8.375e-6]
prob = ODEProblem(f,u0,tspan,p)

sol = solve(prob);

plot(sol,vars=(1,2,3))

# ## Part 2: Investigating Stiffness
#
# Because the reaction rates of `q` vs `s` is very large, this model has a "fast"
# system and a "slow" system. This is typical of ODEs which exhibit a property
# known as stiffness. Stiffness changes the ODE solvers which can handle the
# equation well. [Take a look at the ODE solver page](http://docs.juliadiffeq.org/latest/solvers/ode_solve.html)
# and investigate solving the equation using methods for non-stiff equations
# (ex: `Tsit5`) and stiff equations (ex: `Rodas5`).
#
# Benchmark using $t\in[0,50]$ using `@btime` from BenchmarkTools.jl. What
# happens when you increase the timespan?

# ## (Optional) Part 3: Specifying Analytical Jacobians (I)
#
# Stiff ODE solvers internally utilize the Jacobian of the ODE system in order
# to improve the stepsizes in the solution. However, computing and factorizing
# the Jacobian is costly, and thus it can be beneficial to provide the analytical
# solution.
#
# Use the
# [ODEFunction definition page](http://docs.juliadiffeq.org/latest/features/performance_overloads.html)
# to define an `ODEFunction` which holds both the OREGO ODE and its Jacobian, and solve using `Rodas5`.
#
# ## (Optional) Part 4: Automatic Symbolicification and Analytical Jacobian Calculations
#
# Deriving Jacobians by hand is tedious. Thankfully symbolic mathematical systems
# can do the work for you. And thankfully, DifferentialEquations.jl has tools
# to automatically convert numerical problems into symbolic problems to perform
# the analysis on!
#
# follow the [ModelingToolkit.jl README](https://github.com/JuliaDiffEq/ModelingToolkit.jl)
# to automatically convert your ODE definition
# to its symbolic form using `modelingtoolkitize` and calculate the analytical
# Jacobian. Use the compilation functions to build the `ODEFunction` with the
# embedded analytical solution.
#
# ## Part 5: Adding stochasticity with stochastic differential equations
#
# How does this system react in the presense of stochasticity? We can investigate
# this question by using stochastic differential equations. A stochastic
# differential equation formulation of this model is known as the multiplicative
# noise model, is created with:
#
# $$\begin{align}
# dx &= s(y-xy + x - qx^2)dt + \sigma_1 x dW_1\\
# dy &= \frac{-y - xy + z}{s}dt + \sigma_2 y dW_2\\
# dz &= w(x - z)dt + \sigma_3 z dW_3\end{align}$$
#
# with $\sigma_i = 0.1$ where the `dW` terms describe a Brownian motion, a
# continuous random process with normally distributed increments. Use the
# [tutorial on solving SDEs](http://docs.juliadiffeq.org/latest/tutorials/sde_example.html)
# to solve simulate this model. Then,
# [use the `EnsembleProblem`](http://docs.juliadiffeq.org/latest/features/ensemble.html)
# to generate and plot 100 trajectories of the stochastic model, and use
# `EnsembleSummary` to plot the mean and 5%-95% region over time.
#
# Try solving with the `ImplicitRKMil` and `SOSRI` methods. Notice that it isn't
# stiff every single time!
#
# (For fun, see if you can make the Euler-Maruyama `EM()` method solve this equation.
# This requires a choice of `dt` small enough to be stable. This is the "standard"
# method!)
#
# ## Part 6: Gillespie jump models of discrete stochasticity
#
# When biological models have very few particles, continuous models no longer
# make sense, and instead using the full discrete formulation can be required
# to accuracy describe the dynamics. A discrete differential equation, or
# Gillespie model, is a continuous-time Markov chain with Poisson-distributed
# jumps. A discrete description of the Oregonator model is given by a chemical
# reaction systems:

A+Y -> X+P
X+Y -> 2P
A+X -> 2X + 2Z
2X  -> A + P (note: this has rate kX^2!)
B + Z -> Y

# where reactions take place at a rate which is propoertional to its components,
# i.e. the first reaction has a rate `k*A*Y` for some `k`.
# Use the [tutorial on Gillespie SSA models](http://docs.juliadiffeq.org/latest/tutorials/discrete_stochastic_example.html)
# to implement the `JumpProblem` for this model, and use the `EnsembleProblem`
# and `EnsembleSummary` to characterize the stochastic trajectories.
#
# For what rate constants does the model give the oscillatory dynamics for the
# ODE approximation? For information on the true reaction rates, consult
# [the original paper](https://pubs.acs.org/doi/abs/10.1021/ja00780a001).
#
# ## Part 7: Probabilistic Programming / Bayesian Parameter Estimation with DiffEqBayes.jl + Turing.jl (I)
#
# In many casees, one comes to understand the proper values for their model's
# parameters by utilizing data fitting techniques. In this case, we will use
# the DiffEqBayes.jl library to perform a Bayesian estimation of the parameters.
# For our data we will the following potential output:

t = 0.0:1.0:30.0
data = [1.0 2.05224 2.11422 2.1857 2.26827 2.3641 2.47618 2.60869 2.7677 2.96232 3.20711 3.52709 3.97005 4.64319 5.86202 9.29322 536.068 82388.9 57868.4 1.00399 1.00169 1.00117 1.00094 1.00082 1.00075 1.0007 1.00068 1.00066 1.00065 1.00065 1.00065
        2.0 1.9494 1.89645 1.84227 1.78727 1.73178 1.67601 1.62008 1.56402 1.50772 1.45094 1.39322 1.33366 1.2705 1.19958 1.10651 0.57194 0.180316 0.431409 251.774 591.754 857.464 1062.78 1219.05 1335.56 1419.88 1478.22 1515.63 1536.25 1543.45 1539.98
        3.0 2.82065 2.68703 2.58974 2.52405 2.48644 2.47449 2.48686 2.52337 2.58526 2.67563 2.80053 2.9713 3.21051 3.5712 4.23706 12.0266 14868.8 24987.8 23453.4 19202.2 15721.6 12872.0 10538.8 8628.66 7064.73 5784.29 4735.96 3877.66 3174.94 2599.6]

# [Follow the exmaples on the parameter estimation page](http://docs.juliadiffeq.org/latest/analysis/parameter_estimation.html#Bayesian-Methods-1)
# to perform a Bayesian parameter estimation. What are the most likely parameters
# for the model given the posterior parameter distributions?
#
# Use the `ODEProblem` to perform the fit. If you have time, use the `EnsembleProblem`
# of `SDEProblem`s to perform a fit over averages of the SDE solutions. Note that
# the SDE fit will take significantly more computational resources! See the GPU
# parallelism section for details on how to accelerate this.
#
# ## (Optional) Part 8: Using DiffEqBiological's Reaction Network DSL
#
# DiffEqBiological.jl is a helper library for the DifferentialEquations.jl
# ecosystem for defining chemical reaction systems at a high leevel for easy
# simulation in these various forms. Use the descrption
# [from the Chemical Reaction Networks documentation page](http://docs.juliadiffeq.org/latest/models/biological.html)
# to build a reaction network and generate the ODE/SDE/jump equations, and
# compare the result to your handcoded versions.
#
#


text:
using Documenter, hotrg_frac

makedocs()


text:
### A Pluto.jl notebook ###
# v0.14.2

using Markdown
using InteractiveUtils

# ╔═╡ 29443685-beb0-4102-a46d-63c48b008742
begin
	using LinearAlgebra
	using Plots
	plotly()
end

# ╔═╡ bdeabaa3-8a91-4f9d-8cb7-4527a6385438
md"""
This notebook is a Julia translation of Josh Bevan's [intro to DG YouTube video](https://www.youtube.com/watch?v=sMfPHJUKfaI).
"""

# ╔═╡ dbbb3fea-9dea-11eb-1f7b-afb4a0c5114e
md"""
Here we attempt to solve the 1D scalar conservation eqn of the form:

$ \displaystyle \frac{\partial u}{\partial t} + \frac{\partial f(u)}{\partial x} = 0 $

where $f(u)$ is some flux function describing the "flow" of a conserved
quantity. In this simplified case $f(u) = u$ giving us a linear PDE
We will use periodic boundary conditions to examine numerical dissipation
effects
"""

# ╔═╡ 35708c02-2ef3-432f-963b-62b77f072b19
md"""
Let $u(x,t)$ be the exact solution on the domain $0 \leq x \leq 1$.

Let $u(x,0) = u_0(x) = \sin(2\pi x)$.

Let $V_h$ be the finite vector space of linear polynomials.

Let $u_h$ be the approximate numerical solution consisting of a linear
combination of basis functions ($\psi_n$ for the nth basis) in $V_h$ with
scalar coefficients $a_n$.

Let $\phi_n$ be the nth test function in the same vector space as the basis
functions ($V_h$), in this simplified case N=2.
"""

# ╔═╡ 14c235cf-7b61-45b2-85a8-3d8af45b4edb
N = 2

# ╔═╡ 6be6e442-6983-4629-b734-e6487ac8001d
md"""
Discretize the domain into $K$ elements with $K+1$ nodes.
"""

# ╔═╡ ec7aeb56-230b-4d05-a561-632f71c09e02
struct GridElement{T}
	xl :: T  # Left endpoint
	xr :: T  # Right endpoint
end

# ╔═╡ 34a301b0-58b4-4617-ad47-bd7883b09d8b
struct Grid{T, E}
	K :: Int
	Δ :: T
	elements :: E
end

# ╔═╡ d7c33092-ca55-4e1f-92be-d62aef900393
begin
	K = 32
	Δx = 1/K
	elements = [GridElement(i*Δx, (i+1)*Δx) for i in 0:K-1]
	grid = Grid(K, Δx, elements)
end

# ╔═╡ 28c5d715-c14a-4573-81b5-1e835bc4926c
md"""
Compute and store the initial condition $u_0(x) = \sin(2\pi x)$.
"""

# ╔═╡ d75f11b0-e29b-41d6-87bb-108b0d3cb468
struct FieldElement{D}
	data :: D
end

# ╔═╡ 3f8100e7-dd79-44eb-864c-45ae08c6445a
begin
	struct Field{G, E}
		grid :: G
		elements :: E
	end

	function Field(grid, f::Function)
		field_elements = []
		for grid_element in grid.elements
			field_data = [f(grid_element.xl), f(grid_element.xr)]
			field_element = FieldElement(field_data)
			push!(field_elements, field_element)
		end
		return Field(grid, field_elements)
	end
end

# ╔═╡ b52a8325-a0e0-4697-81b4-26abfbe70a98
begin
	u₀(x) = sin(2π*x)
	u = Field(grid, u₀)
end

# ╔═╡ 45016410-7e25-4f5b-82c5-aa665aaa3c7f
md"""
Let's plot ``u_0(x)`` now!
"""

# ╔═╡ 2f617585-62c1-4f6c-9088-aed874dc43e0
begin
	p = plot()
	for (ge, ue) in zip(grid.elements, u.elements)
		plot!([ge.xl, ge.xr], ue.data, linewidth=2, label="")
	end
	p
end

# ╔═╡ 4de78270-1f76-4f66-9271-34d0f148bc0c
md"""
According to Cockburn & Shu (2001, eq 2.2) let $u_h(x,0)$ be computed by

`` \displaystyle \int u_h(x,0) \phi(x) dx = \int u_0(x) \phi(x) dx ``

for each element (``x_{j-1/2} < x < x_{j+1/2}``).

We can explicitly define a formula for the value of the RHS

`` \displaystyle \int u_0(x) \phi_N(x) dx ``
"""

# ╔═╡ 98e97ab0-6623-4a9e-a9b9-fb80ca182f53
md"""
For each element in the domain $a \le x \le b$, the cell coordinate is $0 \le \xi \le 1$ where

`` x = a + (b - a) \xi  = a + \Delta \xi ``

and

`` \displaystyle \xi = \frac{x - a}{\Delta} ``

so the basis ramp function can be written as

`` \displaystyle \phi_0(\xi) = 1 - \xi \implies \phi_0(x) = \frac{b - x}{\Delta} ``

`` \displaystyle \phi_1(\xi) = \xi \implies \phi_1(x) = \frac{x - a}{\Delta} ``
"""

# ╔═╡ c22c283e-615e-4964-9ff1-f9102f6862b2
md"""
For the left endpoint

``
\displaystyle
\int_a^b u_0(x) \phi_0(x) \; dx
= \int_a^b \sin(2\pi x) \frac{b-x}{\Delta} \; dx
= \frac{2\pi\Delta \cos(2\pi a) + \sin(2\pi a) - \sin(2\pi b)}{4\pi^2\Delta}
``
"""

# ╔═╡ 520fcb4f-bffc-4d53-97b1-703dfaf0b8ca
md"""
For the right endpoint

``
\displaystyle
\int_a^b u_0(x) \phi_1(x) \; dx
= \int_a^b \sin(2\pi x) \frac{x-a}{\Delta} \; dx
= \frac{-2\pi\Delta \cos(2\pi b) + \sin(2\pi b) - \sin(2\pi a)}{4\pi^2\Delta}
``
"""

# ╔═╡ 1096111c-9281-46bc-b5a4-4178a9855cad
begin
	∫u₀ϕ₀dx(a, b, Δ) = (2π*Δ*cos(2π*a) + sin(2π*a) - sin(2π*b)) / (4π^2*Δ)
	∫u₀ϕ₁dx(a, b, Δ) = (-2π*Δ*cos(2π*b) + sin(2π*b) - sin(2π*a)) / (4π^2*Δ)

	∫u₀ϕ_elements = []
	for grid_element in grid.elements
		xₗ = grid_element.xl
		xᵣ = grid_element.xr

		e0 = ∫u₀ϕ₀dx(xₗ, xᵣ, grid.Δ)
		e1 = ∫u₀ϕ₁dx(xₗ, xᵣ, grid.Δ)

		field_element = FieldElement([e0, e1])
		push!(∫u₀ϕ_elements, field_element)
	end

	∫u₀ϕ = Field(grid, ∫u₀ϕ_elements)
end

# ╔═╡ 16e1782b-3948-4d9b-8bfd-fb3134211977
md"""
To find the initial basis weights we need to solve for $a_i$ from

``
\displaystyle
\int_a^b u_0(x) \phi_j(x) \; dx
= \sum_{i=0}^1 a_i \int_a^b \phi_i(x) \phi_j(x) \; dx
``

or from the system of linear equations

``
\begin{pmatrix}
\int u_0 \phi_0
\int u_0 \phi_1
\end{pmatrix}
``
"""

# ╔═╡ 7bc02a02-947b-4526-8bed-49007f3fa8a4
mass_matrix = [1/3 1/6;
               1/6 1/3]

# ╔═╡ 92bdd228-617f-42da-b3d9-b9e9e18a7e4b
begin
	u_basis_elements = []
	for (k, grid_element) in enumerate(grid.elements)
		e = mass_matrix \ ∫u₀ϕ.elements[k].data / grid.Δ
		u_basis_element = FieldElement(e)
		push!(u_basis_elements, u_basis_element)
	end
	u_basis = Field(grid, u_basis_elements)
end

# ╔═╡ 50142bc9-3dcd-4817-b9a6-3c4f37b292b1
begin
	p2 = plot()
	for (ge, ue) in zip(grid.elements, u_basis.elements)
		plot!([ge.xl, ge.xr], ue.data, linewidth=2, label="")
	end
	p2
end

# ╔═╡ 4a103139-e2f9-4182-89b0-9b0b3b4b2503
stiffness_matrix = [0 -1/2 -1/2;
                    0  1/2  1/2]

# ╔═╡ 5c7f8b83-1830-4466-82a2-dd2a6da4115b
upwind_flux = [1 0  0;
               0 0 -1]

# ╔═╡ 2d955948-bf0e-4f89-8bfa-2ffa6b9ac328
stencil = mass_matrix \ (stiffness_matrix + upwind_flux)

# ╔═╡ aed204c6-d68d-449a-aee7-b5fcc537b81a
Δt = 1e-3

# ╔═╡ b7d972b4-82de-4c81-8a9a-80b98850763f
linear_operator = [0 1 0; 0 0 1] + Δt * stencil / grid.Δ

# ╔═╡ 741398fd-cfb6-46af-9bd0-ad047f8f4187
gr()

# ╔═╡ ca064da4-ab88-468e-a2cb-5b76debf6e95
begin
	iterations = 1000
	anim = @animate for i in 1:iterations
		@info "Iteration $i"
		for k in 1:length(grid.elements)
			left_index = k == 1 ? grid.K : k-1

			left_element = u_basis.elements[left_index]
			current_element = u_basis.elements[k]

			u_basis_local = [left_element.data[end], current_element.data...]
			current_element.data .= linear_operator * u_basis_local
		end

		pp = plot()
		for (ge, ue) in zip(grid.elements, u_basis.elements)
			plot!(pp, [ge.xl, ge.xr], ue.data, linewidth=2, label="")
		end
	end every 5
	mp4(anim)
end

# ╔═╡ Cell order:
# ╠═29443685-beb0-4102-a46d-63c48b008742
# ╟─bdeabaa3-8a91-4f9d-8cb7-4527a6385438
# ╟─dbbb3fea-9dea-11eb-1f7b-afb4a0c5114e
# ╟─35708c02-2ef3-432f-963b-62b77f072b19
# ╟─14c235cf-7b61-45b2-85a8-3d8af45b4edb
# ╟─6be6e442-6983-4629-b734-e6487ac8001d
# ╠═ec7aeb56-230b-4d05-a561-632f71c09e02
# ╠═34a301b0-58b4-4617-ad47-bd7883b09d8b
# ╠═d7c33092-ca55-4e1f-92be-d62aef900393
# ╟─28c5d715-c14a-4573-81b5-1e835bc4926c
# ╠═d75f11b0-e29b-41d6-87bb-108b0d3cb468
# ╠═3f8100e7-dd79-44eb-864c-45ae08c6445a
# ╠═b52a8325-a0e0-4697-81b4-26abfbe70a98
# ╟─45016410-7e25-4f5b-82c5-aa665aaa3c7f
# ╠═2f617585-62c1-4f6c-9088-aed874dc43e0
# ╟─4de78270-1f76-4f66-9271-34d0f148bc0c
# ╟─98e97ab0-6623-4a9e-a9b9-fb80ca182f53
# ╟─c22c283e-615e-4964-9ff1-f9102f6862b2
# ╟─520fcb4f-bffc-4d53-97b1-703dfaf0b8ca
# ╠═1096111c-9281-46bc-b5a4-4178a9855cad
# ╟─16e1782b-3948-4d9b-8bfd-fb3134211977
# ╠═7bc02a02-947b-4526-8bed-49007f3fa8a4
# ╠═92bdd228-617f-42da-b3d9-b9e9e18a7e4b
# ╠═50142bc9-3dcd-4817-b9a6-3c4f37b292b1
# ╠═4a103139-e2f9-4182-89b0-9b0b3b4b2503
# ╠═5c7f8b83-1830-4466-82a2-dd2a6da4115b
# ╠═2d955948-bf0e-4f89-8bfa-2ffa6b9ac328
# ╠═aed204c6-d68d-449a-aee7-b5fcc537b81a
# ╠═b7d972b4-82de-4c81-8a9a-80b98850763f
# ╠═741398fd-cfb6-46af-9bd0-ad047f8f4187
# ╠═ca064da4-ab88-468e-a2cb-5b76debf6e95


text:
function makedoctree(d::Dict)
    return XNode(:doctree, [XNode(:ul, [makedoctree(name, val) for (name, val) in d])])

end

function makedoctree(name::String, path)
    link = XNode(:a, Dict(:href => "/$path"), [XLeaf(name)])
    return XNode(:li, [link])
end

function makedoctree(name::String, d::Dict)
    return XNode(:li, [XLeaf(name), XNode(:ul, [makedoctree(n, val) for (n, val) in d])])
end

function loaddoctree(p::AbstractPath)
    xdoctree = Pollen.parse(p)
    xdoctree = cata(xdoctree, SelectTag(:a) & SelectHasAttr(:href)) do x
        href = attributes(x)[:href]
        if !startswith(href, '/')
            href = "/" * href
        end
        attributes(x)[:href] = href
        return x
    end
    return XNode(:doctree, [xdoctree])
end


text:
#== # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# Description
#
#   Tests related to IAU-76/FK5 precession algorithm.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # #
#
# References
#
#   [1] Vallado, D. A (2013). Fundamentals of Astrodynamics and Applications.
#       Microcosm Press, Hawthorn, CA, USA.
#
# # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # # ==#

# File: ./src/transformations/fk5/precession.jl
# =============================================

# Function precession_fk5
# -----------------------

################################################################################
#                                 Test Results
################################################################################
#
# Scenario 01
# ===========
#
# Example 3-15: Performing IAU-76/FK5 reduction.
#
# In this example, using JD_TT = 2453101.828154745, one gets the following data
# related to the precession:
#
#   ζ = 0.0273055˚
#   Θ = 0.0237306˚
#   z = 0.0273059˚
#
# SatelliteToolbox provides the following results:
#
#   julia> (ζ,Θ,z) = precession_fk5(2453101.828154745)
#   julia> ζ*180/pi
#   0.027305539219877804
#
#   julia> Θ*180/pi
#   0.023730619896279084
#
#   julia> z*180/pi
#   0.02730593931829398
#
################################################################################

@testset "Precession" begin
    (ζ,Θ,z) = precession_fk5(2453101.828154745)

    @test ζ*180/pi ≈ 0.0273055 atol=1e-7
    @test Θ*180/pi ≈ 0.0237306 atol=1e-7
    @test z*180/pi ≈ 0.0273059 atol=1e-7
end


text:
module InstrumentOperator

using OffsetArrays
using Polynomials
using NCDatasets
using Interpolations
using Parameters
using DocStringExtensions
using Distributions
using ImageFiltering   # for convolutions
using Unitful
using OrderedCollections
using UnitfulEquivalences
using YAML

@derived_dimension PerTime inv(Unitful.𝐓)

include("types.jl")                   # All types used in this module
include("ils_tables_io.jl")           # IO functions
include("prepare_ils.jl")             # function for ILS preparation
include("instrument_convolutions.jl") # conv using ImageFiltering, can be changed later
include("noise_functions.jl")
include("io_L1_files.jl")

export FixedKernelInstrument, VariableKernelInstrument, FTSInstrument, create_instrument_kernel, conv_spectra

end


text:
module Initial

export initial_sod,initial_static,initial_uniform,initial_one_wave
"""
    this function's inputs are nx, xlim, and gamma

        nx is the number of grid points

        xlim is the limits of the domain

        gamma is the heat capacity ratio

    this function's outputs are ux and uu

        ux has 1 row
        the 1st row is x (x cooridnates)

        uu has 3 rows
        the 1st row is u (velocity)
        the 2nd row is m (mass flow rate)
        the 3rd row is ̂E (total energy per volume)

    this function initialize the Sod shock tube problem (Steger and Warming, 1981)
"""

function initial_sod(nx::Int64, xlim::Any, gamma::Float64)

    Lx = xlim[2]-xlim[1];
    Δx = Lx/(nx-1)

    # Initialization of ux and urhop
    ux = xlim[1]:Δx:xlim[2]
    urhop = Array{Float64,2}(UndefInitializer(), 3,length(ux))

    for i in 1:1:nx
    if (ux[i] < (xlim[2]-xlim[1])/2) urhop[:,i] = [0,1,1]
    else urhop[:,i] = [0,0.125,0.1] end
    end

    #From urhop make uu (u rho*u e)
    uu = Array{Float64,2}(UndefInitializer(), 3,length(ux))
    uu[1,:] = urhop[2,:];
    uu[2,:] = urhop[2,:].*urhop[1,:];
    uu[3,:] = urhop[3,:]/(gamma-1) + 0.5urhop[2,:].*urhop[1,:].^2;

    return ux,uu

end

"""
    this function's inputs are nx, xlim, and gamma

        nx is the number of grid points

        xlim is the limits of the domain

        gamma is the heat capacity ratio

    this function's outputs are ux and uu

        ux has 1 row
        the 1st row is x (x cooridnates)

        uu has 3 rows
        the 1st row is u (velocity)
        the 2nd row is m (mass flow rate)
        the 3rd row is ̂E (total energy per volume)

    this function initialize the Sod shock tube problem (Steger and Warming, 1981)
"""
function initial_static(nx::Int64, xlim::Any, gamma::Float64)

    Lx = xlim[2]-xlim[1];
    Δx = Lx/(nx-1)

    # Initialization of ux and urhop
    ux = xlim[1]:Δx:xlim[2]
    urhop = Array{Float64,2}(UndefInitializer(), 3,length(ux))

    for i in 1:1:nx
    if (ux[i] < (xlim[2]-xlim[1])/2) urhop[:,i] = [0,0.125,1]
    else urhop[:,i] = [0,0.125,1] end
    end

    #From urhop make uu (u rho*u e)
    uu = Array{Float64,2}(UndefInitializer(), 3,length(ux))
    uu[1,:] = urhop[2,:];
    uu[2,:] = urhop[2,:].*urhop[1,:];
    uu[3,:] = urhop[3,:]/(gamma-1) + 0.5urhop[2,:].*urhop[1,:].^2;

    return ux,uu

end

"""
    this function's inputs are nx, xlim, and gamma

        nx is the number of grid points

        xlim is the limits of the domain

        gamma is the heat capacity ratio

    this function's outputs are ux and uu

        ux has 1 row
        the 1st row is x (x cooridnates)

        uu has 3 rows
        the 1st row is u (velocity)
        the 2nd row is m (mass flow rate)
        the 3rd row is ̂E (total energy per volume)

    this function initialize the Sod shock tube problem (Steger and Warming, 1981)
"""
function initial_uniform(nx::Int64, xlim::Any, gamma::Float64)

    Lx = xlim[2]-xlim[1];
    Δx = Lx/(nx-1)

    # Initialization of ux and urhop
    ux = xlim[1]:Δx:xlim[2]
    urhop = Array{Float64,2}(UndefInitializer(), 3,length(ux))

    for i in 1:1:nx
    if (ux[i] < (xlim[2]-xlim[1])/2) urhop[:,i] = [1,0.125,1]
    else urhop[:,i] = [1,0.125,1] end
    end

    #From urhop make uu (u rho*u e)
    uu = Array{Float64,2}(UndefInitializer(), 3,length(ux))
    uu[1,:] = urhop[2,:];
    uu[2,:] = urhop[2,:].*urhop[1,:];
    uu[3,:] = urhop[3,:]/(gamma-1) + 0.5urhop[2,:].*urhop[1,:].^2;

    return ux,uu

end

"""
    this function's inputs are nx, xlim, and gamma

        nx is the number of grid points

        xlim is the limits of the domain

        gamma is the heat capacity ratio

    this function's outputs are ux and uu

        ux has 1 row
        the 1st row is x (x cooridnates)

        uu has 3 rows
        the 1st row is u (velocity)
        the 2nd row is m (mass flow rate)
        the 3rd row is ̂E (total energy per volume)

    this function initialize a wave
"""

    """
    still working on that
    """
function initial_one_wave(nx::Int64, xlim::Any, gamma::Float64)

    Lx = xlim[2]-xlim[1];
    Δx = Lx/(nx-1)

    # Initialization of ux and urhop
    ux = xlim[1]:Δx:xlim[2]
    urhop = Array{Float64,2}(UndefInitializer(), 3,length(ux))

    urhop[1,:].=1
    urhop[2,:].=0.001*exp.(-(collect(ux) .- 0.5.*(ux[end]-ux[1])).^2/0.1^2).+0.125
    urhop[3,:].=0.01*exp.(-(collect(ux) .- 0.5.*(ux[end]-ux[1])).^2/0.1^2).+1

    #From urhop make uu (u rho*u e)
    uu = Array{Float64,2}(UndefInitializer(), 3,length(ux))
    uu[1,:] = urhop[2,:];
    uu[2,:] = urhop[2,:].*urhop[1,:];
    uu[3,:] = urhop[3,:]/(gamma-1) + 0.5urhop[2,:].*urhop[1,:].^2;

    return ux,uu

end

"""
still working on that initial one wave at a relative location
"""
function initial_one_wave(nx::Int64, xlim::Any, gamma::Float64, xposi::Float64)

    Lx = xlim[2]-xlim[1];
    Δx = Lx/(nx-1)

    # Initialization of ux and urhop
    ux = xlim[1]:Δx:xlim[2]
    urhop = Array{Float64,2}(UndefInitializer(), 3,length(ux))

    urhop[1,:].=1
    urhop[2,:].=0.001*exp.(-(collect(ux) .- xposi.*(ux[end]-ux[1])).^2/0.1^2).+0.125
    urhop[3,:].=0.01*exp.(-(collect(ux) .- xposi.*(ux[end]-ux[1])).^2/0.1^2).+1

    #From urhop make uu (u rho*u e)
    uu = Array{Float64,2}(UndefInitializer(), 3,length(ux))
    uu[1,:] = urhop[2,:];
    uu[2,:] = urhop[2,:].*urhop[1,:];
    uu[3,:] = urhop[3,:]/(gamma-1) + 0.5urhop[2,:].*urhop[1,:].^2;

    return ux,uu

end

end


text:
# TODO Maybe take a look at:
#   https://github.com/JuliaString/InternedStrings.jl/blob/master/src/InternedStrings.jl
#
#   The main observation would be that we could potentially do away with WeakNode.
#   Cache becomes `WeakKeyDict{Node,Nothing}` (a.k.a. a weak set)
#   Base.isequal / == should check for _equality_ of ops, and _identitical equality_ of
#   parents.
#
#   Overall a clear benefit would be shrinking the size of the identity map, since we no
#   longer duplicate information.

"""
    WeakNode(parents, op)

Represent a node-like object that doesn't hold strong references to its parents.

This exists purely such that `hash` and `==` *do* allow multiple instances of
`WeakNode` to compare equal if they have the same `parents` and `op`.
"""
struct WeakNode
    parents::NTuple{N,WeakRef} where {N}
    op::NodeOp
end

# Weak nodes need to have hash & equality defined such that instances with equal
# parents and op compare equal. This will be relied upon in `obtain_node` later.
Base.hash(a::WeakNode, h::UInt) = hash(a.op, hash(a.parents, hash(:WeakNode, h)))
function Base.isequal(a::WeakNode, b::WeakNode)
    return isequal(a.parents, b.parents) && isequal(a.op, b.op)
end
Base.:(==)(a::WeakNode, b::WeakNode) = isequal(a, b)

"""
    IdentityMap

An abstract identity map.

Any implementation of this type needs to implement [`obtain_node!`](@ref).
"""
abstract type IdentityMap end

"""
    WeakIdentityMap

Represent a collection of nodes which doesn't hold strong references to any nodes.

This is useful, as it allows the existence of this cache to be somewhat transparent to the
user, and they only have to care about holding on to references for nodes that they care
about.

This structure contains nodes, but also node weak nodes -- these allow us to determine
whether we ought to create a given node.
"""
mutable struct WeakIdentityMap <: IdentityMap
    # TODO This could be wrapped up into a WeakValueDict data structure.
    weak_to_ref::Dict{WeakNode,WeakRef}
    lock::ReentrantLock
    dirty::Bool
    finalizer::Function

    function WeakIdentityMap()
        id_map = new(Dict(), ReentrantLock(), false)
        id_map.finalizer = _ -> id_map.dirty = true
        return id_map
    end
end

Base.lock(f, id_map::WeakIdentityMap) = lock(f, id_map.lock)

function _cleanup!(id_map::WeakIdentityMap)
    id_map.dirty || return nothing

    # We need to clean up stale entries from weak_to_ref.
    id_map.dirty = false

    # This is analogous to the implementation in WeakKeyDict.
    # Note that we use hidden functionality of Dict here. This is because we can no longer
    # rely on the keys to be a good indexer, since they contain weak references that may
    # have gone stale.
    idx = Base.skip_deleted_floor!(id_map.weak_to_ref)
    while idx != 0
        if id_map.weak_to_ref.vals[idx].value === nothing
            Base._delete!(id_map.weak_to_ref, idx)
        end
        idx = Base.skip_deleted(id_map.weak_to_ref, idx + 1)
    end
end

function Base.length(id_map::WeakIdentityMap)
    return lock(id_map) do
        _cleanup!(id_map)
        length(id_map.weak_to_ref)
    end
end
Base.isempty(id_map::WeakIdentityMap) = length(id_map) == 0

function all_nodes(id_map::WeakIdentityMap)
    return lock(id_map) do
        _cleanup!(id_map)
        Node[x.value for x in values(id_map.weak_to_ref)]
    end
end

"""
    _create_node!(id_map::WeakIdentityMap, parents, op, weak_node) -> Node

Create a node and insert it into `id_map`.
"""
function _create_node!(id_map::WeakIdentityMap, parents, op, weak_node::WeakNode)
    node = Node(parents, op)

    # Insert the node & its weak counterpart to the mappings.
    id_map.weak_to_ref[weak_node] = WeakRef(node)

    # Add a finalizer to the node that will declare the id_map to be dirty when it is
    # deleted. We handle this above.
    finalizer(id_map.finalizer, node)

    return node
end

"""
    obtain_node!(id_map::IdentityMap, parents::NTuple{N,Node}, op::NodeOp) -> Node

If a node with `parents` and `op` doesn't exist inside `id_map`, create and insert it.

Return either the new or existing node.
"""
function obtain_node!(
    id_map::WeakIdentityMap, parents::NTuple{N,Node}, op::NodeOp
) where {N}
    return lock(id_map) do
        # Before attempting to query or modify the id_map, ensure it is free of dangling
        # references.
        _cleanup!(id_map)

        weak_node = WeakNode(map(WeakRef, parents), op)
        node_ref = get(id_map.weak_to_ref, weak_node, nothing)
        if !isnothing(node_ref)
            # Remember that we need to unwrap the value from the WeakRef...
            node = node_ref.value

            # This is inspired by iterate(::WeakKeyDict, ...) in base/weakkeydict.jl ,
            # and included due to otherwise seeing transient issues of `node_ref.value`
            # being `nothing`.
            GC.safepoint()  # ensure `node` is now gc-rooted
            isnothing(node) ? _create_node!(id_map, parents, op, weak_node) : node
        else
            # An equivalent node does not yet exist in the id_map; create it.
            _create_node!(id_map, parents, op, weak_node)
        end
    end
end

# This is the single instance of the id_map that we want
const _GLOBAL_IDENTITY_MAP = WeakIdentityMap()

"""
    global_identity_map() -> IdentityMap

Get the global IdentityMap instance used in TimeDag.
"""
global_identity_map() = _GLOBAL_IDENTITY_MAP


text:
function WRITE(f, x)
    for i in 1:length(x)
        write(f, x[i])
    end
end

function WRITE_CHUNK(rchunk, f, ds, n, ff, delim, quotechar)
    chunk = div(n, rchunk)
    for j in 1:chunk
        WRITE(f, InMemoryDatasets.row_join(view(ds, (j-1)*rchunk+1:j*rchunk, :), ff, :, delim = delim, quotechar = quotechar))
    end
    WRITE(f, InMemoryDatasets.row_join(view(ds, chunk*rchunk+1:n, :), ff, :, delim = delim, quotechar = quotechar))
end

# basic function for writing csv files
function filewriter(path::AbstractString, ds::AbstractDataset; delim = ',', quotechar = nothing, mapformats = false, append = false, header = true)
    ncols = InMemoryDatasets.ncol(ds)
    ff = Function[]
    if mapformats
        for j in 1:ncols
            push!(ff, InMemoryDatasets.getformat(ds, j))
        end
    else
        ff = repeat([identity], ncols)
    end
    n, p = size(ds)
    f = open(path, write = true, append = append)
    if header
        allnames = names(ds)
        for j in 1:ncols
            if quotechar !== nothing
                write(f, quotechar)
            end
            write(f, allnames[j])
            if quotechar !== nothing
                write(f, quotechar)
            end
            if j == ncols
                write(f, '\n')
            else
                write(f, delim)
            end
        end
    end
    if p < 1000
        WRITE_CHUNK(1000, f, ds, n, ff, delim, quotechar)
    elseif p<10000
        WRITE_CHUNK(500, f, ds, n, ff, delim, quotechar)
    elseif p<50000
        WRITE_CHUNK(100, f, ds, n, ff, delim, quotechar)
    else
        WRITE_CHUNK(20, f, ds, n, ff, delim, quotechar)
    end
    close(f)
end


text:
#------------------------------
# Runtime and complexity analysis
# for linear ODE system with two unkown parameters
#------------------------------

#------------------------------
# Setup
#------------------------------

using Distributed

addprocs(5)

@everywhere using Pkg;
#@everywhere cd("/Users/imbi-mac-102/Desktop/ManuscriptMasterThesis/code/") # omit or change accordingly on server
@everywhere Pkg.activate(".")
Pkg.instantiate()

@everywhere using BenchmarkTools
@everywhere using DataFrames
@everywhere using Distributed
@everywhere using Distributions
@everywhere using Random
@everywhere using Flux
@everywhere using DiffEqFlux
@everywhere using OrdinaryDiffEq
@everywhere using SharedArrays

@everywhere include("../src/simulation.jl")
@everywhere include("../src/model.jl")
@everywhere include("../src/benchmarking.jl")

#------------------------------
# Define and obtain ground-truth developments
#------------------------------

# define initial condition
@everywhere true_u0 = Float32[2, 1]
# define time span on which to solve the ODE
@everywhere tspan = (0.0f0, 10.0f0)
# define parameters for the two distinct groups
@everywhere true_odeparams_group1 = Float32[-0.2, 0.00, 0.00, -0.2]
@everywhere true_odeparams_group2 = Float32[-0.2, 0.00, 0.00, 0.2]

# define corresponding ODE problems for the two groups
@everywhere prob1 = ODEProblem(linear_2d_system,true_u0,tspan,true_odeparams_group1)
@everywhere prob2 = ODEProblem(linear_2d_system,true_u0,tspan,true_odeparams_group2)

# solve ODE systems to obtain "true" underlying trajectory in each group
@everywhere dt=0.1
@everywhere sol_group1 = solve(prob1, Tsit5(), saveat = dt);
@everywhere sol_group2 = solve(prob2, Tsit5(), saveat = dt);

#------------------------------
# Train model using benchmarking
#------------------------------

# define number of observations, variables and baseline variables to try
@everywhere n_obs = [50, 100, 250, 500, 1000, 2000, 5000]
@everywhere n_vars = [10, 20, 50, 100, 200]
@everywhere n_baselinevars = [10, 20, 50, 100, 200]
@everywhere lenobs, lenvars, lenbvars = length(n_obs), length(n_vars), length(n_baselinevars)

# construct dataframe: n, p, q, time, memory, allocations
benchmarkdf = DataFrame(n_obs = cat(n_obs, fill(100, lenvars + lenbvars), dims=1),
        n_vars = cat(fill(10, lenobs), n_vars, fill(10, lenbvars), dims=1),
        n_baselinevars = cat(fill(50, lenobs + lenvars), n_baselinevars, dims=1),
        time = fill(0.0, lenobs+lenvars+lenbvars),
        gctime = fill(0.0, lenobs+lenvars+lenbvars),
        memory = fill(0, lenobs+lenvars+lenbvars),
        allocs = fill(0, lenobs+lenvars+lenbvars)
)

# get it to a shared array for distributed computing
benchmarkarray = SharedArrays.SharedMatrix{Float64}(size(Matrix(benchmarkdf)));
benchmarkarray[:,1:3] = Matrix(benchmarkdf)[:,1:3];

@everywhere eval($benchmarkarray);

# scenario 1: fixed number of time-dep and baseline variables, varying number of observations
@sync @distributed for n_ind in 1:lenobs
    # warmup (first run takes longer because of precompilation times and shouldnt be included)
    n_warmup, p_warmup, q_warmup, q_info_warmup = 100, 10, 10, 10
    xs, x_baseline, tvals, group1, group2 = generate_all(n_warmup, p_warmup, q_warmup, q_info_warmup);
    trainingdata = zip(xs, x_baseline, tvals);
    zdim = nODEparams = 2
    m = init_vae(p_warmup, q_warmup, zdim, nODEparams, prob1)
    L = loss_wrapper(m)
    ps = getparams(m)
    opt = ADAM(0.0005)
    for epoch in 1:35
        Flux.train!(L, ps, trainingdata, opt)
    end
    println("warmup done")
    # now start for real
    n, p, q = n_obs[n_ind], 10, 50
    println("n=$n, p=$p, q=$q")
    q_info = Int(q/5)
    xs, x_baseline, tvals, group1, group2 = generate_all(n, p, q, q_info);
    trainingdata = zip(xs, x_baseline, tvals);
    zdim = nODEparams = 2
    m = init_vae(p, q, zdim, nODEparams, prob1)
    b = @benchmark run_benchmark($trainingdata, $m) samples=1 evals=1
    println("training done")
    row = n_ind
    benchmarkarray[row,4] = b.times[1] # times
    benchmarkarray[row,5] = b.gctimes[1] # gctimes
    benchmarkarray[row,6] = b.memory # memory
    benchmarkarray[row,7] = b.allocs # allocations
end

# scenario 2: fixed number of observations and baseline variables, varying number of time-dependent variables
@sync @distributed for p_ind in 1:lenvars
    n, p, q = 100, n_vars[p_ind], 50
    println("n=$n, p=$p, q=$q")
    q_info = Int(q/5)
    xs, x_baseline, tvals, group1, group2 = generate_all(n, p, q, q_info);
    trainingdata = zip(xs, x_baseline, tvals);
    zdim = nODEparams = 2
    m = init_vae(p, q, zdim, nODEparams, prob1)
    b = @benchmark run_benchmark($trainingdata, $m) samples=1 evals=1
    println("training done")
    row = lenobs + p_ind
    benchmarkarray[row,4] = b.times[1] # times
    benchmarkarray[row,5] = b.gctimes[1] # gctimes
    benchmarkarray[row,6] = b.memory # memory
    benchmarkarray[row,7] = b.allocs # allocations
end

# scenario 3: fixed number of observations and time-dependent variables, varying number of baseline variables
@sync @distributed for q_ind in 1:lenbvars
    n, p, q = 100, 10, n_baselinevars[q_ind]
    println("n=$n, p=$p, q=$q")
    q_info = Int(q/5)
    xs, x_baseline, tvals, group1, group2 = generate_all(n, p, q, q_info);
    trainingdata = zip(xs, x_baseline, tvals);
    zdim = nODEparams = 2
    m = init_vae(p, q, zdim, nODEparams, prob1)
    b = @benchmark run_benchmark($trainingdata, $m) samples=1 evals=1
    println("training done")
    row = lenobs + lenvars + q_ind
    benchmarkarray[row,4] = b.times[1] # times
    benchmarkarray[row,5] = b.gctimes[1] # gctimes
    benchmarkarray[row,6] = b.memory # memory
    benchmarkarray[row,7] = b.allocs # allocations
end

# if desired: save as JLD2 file
using JLD2
JLD2.@save "benchmarkresults.jld2" benchmarkarray
# and re-load from saved
JLD2.@load "benchmarkresults.jld2"
benchmarkarray = eval(:benchmarkarray)

# copy back to dataframe, to be saved later as CSV
benchmarkdf[:,:time] = benchmarkarray[:,4]
benchmarkdf[:,:gctime] = benchmarkarray[:,5]
benchmarkdf[:,:memory] = benchmarkarray[:,6]
benchmarkdf[:,:allocs] = benchmarkarray[:,7]
benchmarkdf[:,:time_in_seconds] = round.(benchmarkarray[:,4] .* 1e-9, digits=3)

# turn memory into human-readable format (taken from BenchmarkTools.jl source code)
benchmarkdf[:,:prettymemory] = prettymemory.(benchmarkarray[:,6])

# save entire dataframe as CSV
using CSV
CSV.write("benchmarkresults.csv", benchmarkdf)

# extract tables as in manuscript appendix (and save as CSV files)

# different number of observations for fixed p (10) and q (50)
benchmark_obs = benchmarkdf[1:7, [:n_obs, :time_in_seconds, :prettymemory]]
CSV.write("benchmarkresults_obs.csv", benchmark_obs)

benchmark_vars = benchmarkdf[8:12, [:n_vars, :time_in_seconds, :prettymemory]]
CSV.write("benchmarkresults_vars.csv", benchmark_vars)

benchmark_bvars = benchmarkdf[13:17, [:n_baselinevars, :time_in_seconds, :prettymemory]]
CSV.write("benchmarkresults_baselinevars.csv", benchmark_bvars)