text:
#=
    High-level Overview

    Hyperscript is a package for representing and rendering HTML and SVG DOM trees.

    The primary abstraction is a `Node` — with a tag, attributes, and children.

    Each `Node` also has a `Context` associated with it which defines the way in which
    the basic Hyperscript functions operate on that node.

    The basic pipeline is as follows:

    A `Node` is created with some tag, attributes, and children.

    The tag, attributes, and children are normalized and then validated, using the
    functions [normalize|validate][tag|attr|child] which can dispatch on context to
    handle different node types according to their needs.

    Normalization is the conversion to a canonical form (e.g. kebab- or camel-casing) and
    validation is checking invariants and throwing an error if they are not met (e.g. no
    spaces in attribute names).

    Each of these functions takes a `ctx` argument, allowing the context of a node to fully
    dictate the terms of normalization, validation, escapining, and rendering.

    Nodes handle escaping of their contents upon render using functions named
    escape[tag|attrname|attrvalue|child].

    The full pipeline over the lifecycle of a `Node` involves all of these functions
    in order: normalize => validate => escape upon render.

    ===

    Node extensibility points

    Note that these defaults are not defined — it makes more sense to fully define the relevant
    behaviors independently for each context. The values below are just 'no-op' examples.

    # Return the normalized property value
    normalizetag(ctx, tag) = tag
    normalizeattr(ctx, tag, attr) = attr
    normalizechild(ctx, tag, child) = child

    # Return the property value or throw a validation error
    validatetag(ctx, tag) = tag
    validateattr(ctx, tag, attr) = attr
    validatechild(ctx, tag, child) = child

    # Return a Dict{Char, String} specifying escape replacements
    escapetag(ctx) = NO_ESCAPES
    escapeattrname(ctx) = NO_ESCAPES
    escapeattrvalue(ctx) = NO_ESCAPES
    escapechild(ctx) = NO_ESCAPES
=#

module Hyperscript

export @tags, @tags_noescape, m, css, Style, styles, render, Pretty, savehtml, savesvg

# Units
include(joinpath(@__DIR__, "cssunits.jl"))
export px, pt, em, #= (this one conflicts with Base) rem, =# vh, vw, vmin, vmax, pc

## Basic definitions

abstract type Context end

struct CSS <: Context
    allow_nan_attr_values::Bool
end

struct HTMLSVG <: Context
    allow_nan_attr_values::Bool
    noescape::Bool
end

abstract type AbstractNode{T} end

struct Node{T<:Context} <: AbstractNode{T}
    context::T
    tag::String
    children::Vector{Any}
    attrs::Dict{String, Any}
end

function Node(ctx::T, tag::AbstractString, children, attrs) where T <: Context
    tag = validatetag(ctx, normalizetag(ctx, tag))
    Node{T}(
        ctx,
        tag,
        processchildren(ctx, tag, children),
        processattrs(ctx, tag, attrs)
    )
end

function (node::Node{T})(cs...; as...) where T
    ctx = context(node)
    Node{T}(
        ctx,
        tag(node),
        isempty(cs) ? children(node) : prepend!(processchildren(ctx, tag(node), cs), children(node)),
        isempty(as) ? attrs(node)    : merge(attrs(node), processattrs(ctx, tag(node), as))
    )
end

tag(x::Node)      = Base.getfield(x, :tag)
attrs(x::Node)    = Base.getfield(x, :attrs)
children(x::Node) = Base.getfield(x, :children)
context(x::Node)  = Base.getfield(x, :context)

function Base.:(==)(x::Node, y::Node)
    context(x)  == context(y)  && tag(x)   == tag(y) &&
    children(x) == children(y) && attrs(x) == attrs(y)
end

## Node utils

function processchildren(ctx, tag, children)
    # Any[] for type-stability Node construction (children::Vector{Any})
    Any[validatechild(ctx, tag, normalizechild(ctx, tag, child)) for child in flat(children)]
end

# A single attribute is allowed to normalize to multiple attributes,
# for example when normalizing CSS attribute names into vendor-prefixed versions.
function processattrs(ctx, tag, attrs)
    Dict{String, Any}(
        validateattr(ctx, tag, attr′)
        for attr in attrs
        for attr′ in flat(normalizeattr(ctx, tag, attr))
    )
end

function flat(xs::Union{Base.Generator, Tuple, AbstractArray})
    out = Any[] # Vector{Any} for node children and attribute values
    for x in xs
        append!(out, flat(x))
    end
    out
end
flat(x) = (x,)
flat(x::AbstractRange) = (x,)

## Rendering

struct RenderContext
    pretty::Bool
    indent::String
    level::Int
end
const DEFAULT_RCTX = RenderContext(false, "  ", 0)

"""
Wrapper struct for pretty-printing `Node`s: `Pretty(node)`.
Line feeds are added along with indentation controlled by `indent`.
"""
struct Pretty{T <: AbstractNode}
    node::T
    rctx::RenderContext
    Pretty(node::T; indent="  ") where {T} = new{T}(node, RenderContext(true, indent, 0))
    Pretty(node::T, rctx::RenderContext) where {T} = new{T}(node, rctx)
end

# Top-level nodes render in their own node context.
render(io::IO, rctx::RenderContext, x::Node) = render(io, rctx, context(x), x)
render(io::IO, x::Pretty) = render(io, x.rctx, x.node)
render(io::IO, x::Node) = render(io, DEFAULT_RCTX, x)
render(x::AbstractNode) = sprint(render, x)

Base.show(io::IO, node::Union{Node, Pretty}) = render(io, node)
Base.show(io::IO, m::MIME"text/html", node::Union{Node, Pretty}) = render(io, node)

printescaped(io::IO, x::AbstractString, escapes) = for c in x
    print(io, get(escapes, c, c))
end

printescaped(io::IO, x::AbstractChar, escapes) = printescaped(io, string(x), escapes)

# todo: turn the above into something like an escaping IO pipe to avoid string
# allocation via sprint. future use: sprint(printescaped, x, escapes))
printescaped(io::IO, x, escapes) = printescaped(io, sprint(print, x), escapes)

# pass numbers through untrammelled
kebab(camel::String) = join(islowercase(c) || isnumeric(c) || c == '-' ? c : '-' * lowercase(c) for c in camel)

## HTMLSVG

function render(io::IO, rctx::RenderContext, ctx::HTMLSVG, node::Node{HTMLSVG})
    etag = escapetag(ctx)
    eattrname = escapeattrname(ctx)
    eattrvalue = escapeattrvalue(ctx)
    if rctx.pretty && rctx.level > 0
        print(io, "\n", rctx.indent ^ rctx.level, "<")
    else
        print(io, "<")
    end
    printescaped(io, tag(node), etag)
    for (name, value) in pairs(attrs(node))
        print(io, " ")
        printescaped(io, name, eattrname)
        if value != nothing
            print(io, "=\"")
            printescaped(io, value, eattrvalue)
            print(io, "\"")
        end
    end

    if isvoid(tag(node))
        @assert isempty(children(node))
        print(io, " />")
    else
        print(io, ">")
        for child in children(node)
            renderdomchild(io, RenderContext(rctx.pretty, rctx.indent, rctx.level + 1), ctx, child)
        end
        if rctx.pretty && any(x -> isa(x, AbstractNode), children(node))
            print(io, "\n", rctx.indent ^ rctx.level, "</")
        else
            print(io, "</")
        end
        printescaped(io, tag(node), etag)
        print(io, ">")
    end
end

const VOID_TAGS = Set([
    "track", "hr", "col", "embed", "br", "circle", "input", "base",
    "use", "source", "polyline", "param", "ellipse", "link", "img",
    "path", "wbr", "line", "stop", "rect", "area", "meta", "polygon"
])
isvoid(tag) = tag ∈ VOID_TAGS

# Rendering HTMLSVG child nodes in their own context
renderdomchild(io, rctx::RenderContext, ctx::HTMLSVG, node::AbstractNode{HTMLSVG}) = render(io, rctx, node)

# Do nothing for `nothing`; this is similar to using `nothing` as an attribute value for valueless attributes.
renderdomchild(io, rctx::RenderContext, ctx, x::Nothing) = nothing

# Render and escape other HTMLSVG children, including CSS nodes, in the parent context.
# If a child is `showable` with text/html, render with that using `repr`.
renderdomchild(io, rctx::RenderContext, ctx, x) =
    showable(MIME("text/html"), x) ? show(io, MIME("text/html"), x) : printescaped(io, x, escapechild(ctx))

# All camelCase attribute names from HTML 4, HTML 5, SVG 1.1, SVG Tiny 1.2, and SVG 2
const HTML_SVG_CAMELS = Dict(lowercase(x) => x for x in [
    "preserveAspectRatio", "requiredExtensions", "systemLanguage",
    "externalResourcesRequired", "attributeName", "attributeType", "calcMode",
    "keySplines", "keyTimes", "repeatCount", "repeatDur", "requiredFeatures",
    "requiredFonts", "requiredFormats", "baseFrequency", "numOctaves", "stitchTiles",
    "focusHighlight", "lengthAdjust", "textLength", "glyphRef", "gradientTransform",
    "gradientUnits", "spreadMethod", "tableValues", "pathLength", "clipPathUnits",
    "stdDeviation", "viewBox", "viewTarget", "zoomAndPan", "initialVisibility",
    "syncBehavior", "syncMaster", "syncTolerance", "transformBehavior", "keyPoints",
    "defaultAction", "startOffset", "mediaCharacterEncoding", "mediaContentEncodings",
    "mediaSize", "mediaTime", "maskContentUnits", "maskUnits", "baseProfile",
    "contentScriptType", "contentStyleType", "playbackOrder", "snapshotTime",
    "syncBehaviorDefault", "syncToleranceDefault", "timelineBegin", "edgeMode",
    "kernelMatrix", "kernelUnitLength", "preserveAlpha", "targetX", "targetY",
    "patternContentUnits", "patternTransform", "patternUnits", "xChannelSelector",
    "yChannelSelector", "diffuseConstant", "surfaceScale", "refX", "refY",
    "markerHeight", "markerUnits", "markerWidth", "filterRes", "filterUnits",
    "primitiveUnits", "specularConstant", "specularExponent", "limitingConeAngle",
    "pointsAtX", "pointsAtY", "pointsAtZ", "hatchContentUnits", "hatchUnits"])

normalizetag(ctx::HTMLSVG, tag) = strip(tag)

# The simplest normalization — kebab-case and don't pay attention to the tag.
# Allows both squishcase and camelCase for the attributes above.
# If the attribute name is a string and not a Symbol (using the Node constructor),
# then no normalization is performed — this way you can pass any attribute you'd like.
function normalizeattr(ctx::HTMLSVG, tag, (name, value)::Pair{Symbol, <:Any})
    name = string(name)
    get(() -> kebab(name), HTML_SVG_CAMELS, lowercase(name)) => value
end

function normalizeattr(ctx::HTMLSVG, tag, attr::Pair{<:AbstractString, <:Any})
    # Note: This must implementation must change if we begin to normalize attr values above.
    # Right now we only normalize attr names.
    attr
end

normalizechild(ctx::HTMLSVG, tag, child) = child

# Nice printing in errors
stringify(ctx::HTMLSVG, tag, attr::String=" ") = "<$tag>$attr $(isvoid(tag) ? " />" : ">")"
stringify(ctx::HTMLSVG, tag, (name, value)::Pair) = stringify(ctx, tag, " $name=$value")

function validatetag(ctx::CSS, tag)
    isempty(tag) && error("Tag cannot be empty.")
    tag
end

function validateattr(ctx::HTMLSVG, tag, attr)
    (name, value) = attr
    if !ctx.allow_nan_attr_values && typeof(value) <: AbstractFloat && isnan(value)
        error("NaN values are not allowed for HTML or SVG nodes: $(stringify(ctx, tag, attr))")
    end
    if any(isspace, name)
        error("Spaces are not allowed in HTML or SVG attribute names: $(stringify(ctx, tag, attr))")
    end
    attr
end

function validatechild(ctx::HTMLSVG, tag, child)
    if isvoid(tag)
        error("Void tags are not allowed to have children: $(stringify(ctx, tag))")
    end
    child
end

# Creates an HTML or SVG escaping dictionary
chardict(chars) = Dict(c => "&#$(Int(c));" for c in chars)

# See: https://stackoverflow.com/questions/7753448/how-do-i-escape-quotes-in-html-attribute-values
const ATTR_VALUE_ESCAPES = chardict("&<>\"\n\r\t")

# See: https://stackoverflow.com/a/9189067/1175713
const HTML_ESCAPES = chardict("&<>\"'`!@\$%()=+{}[]")

# Used for CSS nodes, as well as children of tag nodes defined with @tags_noescape
const NO_ESCAPES = Dict{Char, String}()

escapetag(ctx::HTMLSVG) = HTML_ESCAPES
escapeattrname(ctx::HTMLSVG) = HTML_ESCAPES
escapeattrvalue(ctx::HTMLSVG) = ATTR_VALUE_ESCAPES
escapechild(ctx::HTMLSVG) = ctx.noescape ? NO_ESCAPES : HTML_ESCAPES

# Concise CSS class shorthand
addclass(attrs, class) = haskey(attrs, "class") ? string(attrs["class"], " ", class) : class
Base.getproperty(x::Node{HTMLSVG}, class::Symbol) = x(class=addclass(attrs(x), kebab(String(class))))
Base.getproperty(x::Node{HTMLSVG}, class::String) = x(class=addclass(attrs(x), class))

const DEFAULT_HTMLSVG_CONTEXT = HTMLSVG(false, false)
const NOESCAPE_HTMLSVG_CONTEXT = HTMLSVG(false, true)
m(tag::AbstractString, cs...; as...) = Node(DEFAULT_HTMLSVG_CONTEXT, tag, cs, as)
m(ctx::Context, tag::AbstractString, cs...; as...) = Node(ctx, tag, cs, as)

# HTML/SVG tags macros
macro tags(args::Symbol...)
    blk = Expr(:block)
    for tag in args
        push!(blk.args, quote
            const $(esc(tag)) = m($(string(tag)))
        end)
    end
    push!(blk.args, nothing)
    blk
end

macro tags_noescape(args::Symbol...)
    blk = Expr(:block)
    for tag in args
        push!(blk.args, quote
            const $(esc(tag)) = m(NOESCAPE_HTMLSVG_CONTEXT, $(string(tag)))
        end)
    end
    push!(blk.args, nothing)
    blk
end

## CSS

ismedia(node::Node{CSS}) = startswith(tag(node), "@media")
nestchildren(node::Node{CSS}) = startswith(tag(node), "@")

function render(io::IO, rctx::RenderContext, ctx::CSS, node::Node)
    @assert ctx == context(node)

    etag = escapetag(ctx)
    eattrname = escapeattrname(ctx)
    eattrvalue = escapeattrvalue(ctx)

    if rctx.pretty
        print(io, "\n", rctx.indent ^ rctx.level)
    end
    printescaped(io, tag(node), etag)
    print(io, " {") # \n

    for (name, value) in pairs(attrs(node))
        printescaped(io, name, eattrname)
        print(io, ": ")
        printescaped(io, value, eattrvalue)
        print(io, ";") # \n
    end

    for child in children(node)
        @assert typeof(child) <: Node{CSS}  "CSS child elements must be `Node`s."
    end

    nest = nestchildren(node)
    nest && for child in children(node)
        render(io, rctx, child)
    end

    print(io, "}") # \n

    !nest && for child in children(node)
        childctx = context(child)
        render(io, rctx, Node{typeof(childctx)}(childctx, tag(node) * " " * tag(child), children(child), attrs(child)))
    end
end

renderdomchild(io, rctx::RenderContext, ctx::Context, node::AbstractNode{CSS}) = render(io, rctx, node)

normalizetag(ctx::CSS, tag) = strip(tag)

stringify(ctx::CSS, tag, (name, value)::Pair) = "$tag { $name: $value; }"

function validatetag(ctx::HTMLSVG, tag)
    isempty(tag) && error("Tag cannot be empty.")
    tag
end

function validateattr(ctx::CSS, tag, attr)
    name, value = attr
    last(attr) == nothing && error("CSS attribute value may not be `nothing`: $(stringify(ctx, tag, attr))")
    last(attr) == "" && error("CSS attribute value may not be the empty string: $(stringify(ctx, tag, attr))")
    if !ctx.allow_nan_attr_values && typeof(value) <: AbstractFloat && isnan(value)
        error("NaN values are not allowed for CSS nodes: $(stringify(ctx, tag, attr))")
    end
    attr
end

function validatechild(ctx::CSS, tag, child)
    typeof(child) <: Node{CSS} || error("CSS nodes may only have Node{CSS} children. Found $(typeof(child)): $child")
    child
end

normalizeattr(ctx::CSS, tag, attr::Pair) = kebab(string(first(attr))) => last(attr)
normalizechild(ctx::CSS, tag, child) = child

escapetag(ctx::CSS) = NO_ESCAPES
escapeattrname(ctx::CSS) = NO_ESCAPES
escapeattrvalue(ctx::CSS) = NO_ESCAPES

const DEFAULT_CSS_CONTEXT = CSS(false)
css(tag, children...; attrs...) = Node(DEFAULT_CSS_CONTEXT, tag, children, attrs)

## Scoped CSS

# A `Styled` node results from the application of a `Style` to a `Node`.
# It serves as a cascade barrier — parent styles do not bleed into nested styled nodes.
struct Styled <: AbstractNode{HTMLSVG}
    node::Node{HTMLSVG}
    style
end

# delegate
tag(x::Styled) = tag(x.node)
attrs(x::Styled) = attrs(x.node)
children(x::Styled) = children(x.node)
context(x::Styled) = context(x.node)
(x::Styled)(cs...; as...) = Styled(x.node((augmentdom(x.style.id, c) for c in  cs)...; as...), x.style)
render(io::IO, x::Styled) = render(io, x.node)
render(io::IO, rctx::RenderContext, x::Styled) = render(io, rctx, x.node)

Base.show(io::IO, x::Styled) = show(io, x.node)
Base.show(io::IO, m::MIME"text/html", x::Styled) = show(io, x.node)

struct Style
    id::Int
    styles::Vector{Node{CSS}}
    augmentcss(id, node) = Node{CSS}(
        context(node),
        isempty(attrs(node)) || ismedia(node) ? tag(node) : tag(node) * "[v-style$id]",
        augmentcss.(id, children(node)),
        attrs(node)
    )
    Style(id::Int, styles) = new(id, [augmentcss(id, node) for node in styles])
end

style_id = 0
function Style(styles...)
    global style_id
    Style(style_id += 1, styles)
end

styles(x::Style) = x.styles

render(io::IO, rctx::RenderContext, x::Style) = for node in x.styles
    render(io, rctx, node)
end

augmentdom(id, x) = x # Literals and other non-HTML/SVG objects
augmentdom(id, x::Styled) = x # `Styled` nodes act as cascade barriers
augmentdom(id, node::Node{T}) where {T} = Node{T}(
    context(node),
    tag(node),
    augmentdom.(id, children(node)),
    push!(copy(attrs(node)), "v-style$id" => nothing) # note: makes a defensive copy
)
(s::Style)(x::Node) = Styled(augmentdom(s.id, x), s)

## Useful utilities for generating HTML files

# todo: can these two functions be merged?
# todo: verify that this doesn't double-wrap...
# I think I've seen it happen in the wild.
function wraphtml(dom)
    node = if dom isa Node && tag(dom) == "html"
        dom
    else
        # todo: head, body?
        m("html", dom)
    end
    preamble = "<!doctype html>\n" # HTML 5
    string(preamble, sprint(show, node))
end
savehtml(filename, children...) = write(filename, wraphtml(children))

function wrapsvg(dom)
    preamble = """
    <?xml version="1.0" standalone="yes"?>
    <!DOCTYPE svg PUBLIC "-//W3C//DTD SVG 1.1//EN"
      "http://www.w3.org/Graphics/SVG/1.1/DTD/svg11.dtd">
    """ # SVG 1.1
    node = tag(dom) == "svg" ? dom : m("svg", dom)
    haskey(attrs(node), "xmlns")   || (node = node(xmlns="http://www.w3.org/2000/svg"))
    haskey(attrs(node), "version") || (node = node(version="1.1"))
    string(preamble, sprint(show, node))
end
savesvg(filename, children) = write(filename, wrapsvg(children))

#=
future enhancements
    - some way to import many common HTML/SVG tags at once
    - when applying a Style to a node only add the `v-style` marker to those nodes that may be affected by a style selector.
    - add linting validations for e.g. <circle x=... />
    - autoprefix css attributes based on some criterion, perhaps from caniuse.com
    process m(k => v) as m(; k=v)? this would break the "all normal arguments are children" invariant.
    - rename save[x] to write[x]
=#

end # module


text:
# Autogenerated wrapper script for PTSCOTCH_jll for armv6l-linux-gnueabihf
export libptesmumps, libptscotch, libptscotcherr, libptscotcherrexit, libptscotchparmetis

using Zlib_jll
using MPICH_jll
using SCOTCH_jll
JLLWrappers.@generate_wrapper_header("PTSCOTCH")
JLLWrappers.@declare_library_product(libptesmumps, "libptesmumps.so")
JLLWrappers.@declare_library_product(libptscotch, "libptscotch.so")
JLLWrappers.@declare_library_product(libptscotcherr, "libptscotcherr.so")
JLLWrappers.@declare_library_product(libptscotcherrexit, "libptscotcherrexit.so")
JLLWrappers.@declare_library_product(libptscotchparmetis, "libptscotchparmetis.so")
function __init__()
    JLLWrappers.@generate_init_header(Zlib_jll, MPICH_jll, SCOTCH_jll)
    JLLWrappers.@init_library_product(
        libptesmumps,
        "lib/libptesmumps.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libptscotch,
        "lib/libptscotch.so",
        nothing,
    )

    JLLWrappers.@init_library_product(
        libptscotcherr,
        "lib/libptscotcherr.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libptscotcherrexit,
        "lib/libptscotcherrexit.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@init_library_product(
        libptscotchparmetis,
        "lib/libptscotchparmetis.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()

ext: jl


text:
# To add support for loading/storing...
@inline function load(p::Ptr{UInt}, ::Type{T}) where {T<:NativeTypes}
    __vload(Base.unsafe_convert(Ptr{T}, p), False(), register_size())
end
@inline function load(p::Ptr{UInt}, ::Type{T}) where {T<:Union{Ptr,Core.LLVMPtr}}
    reinterpret(T, __vload(p, False(), register_size()))
end
# @inline function load(p::Ptr{UInt}, ::Type{T}) where {T<:NativeTypes}
#     __vload(reinterpret(Core.LLVMPtr{T,0}, p), False(), register_size())
# end
# @inline function load(p::Ptr{UInt}, ::Type{T}) where {T<:Union{Ptr,Core.LLVMPtr}}
#     reinterpret(T, __vload(reinterpret(Core.LLVMPtr{UInt,0}, p), False(), register_size()))
# end
@inline load(p::Ptr{UInt}, ::Type{T}) where {T} = unsafe_load(Base.unsafe_convert(Ptr{T}, p))
@inline function store!(p::Ptr{UInt}, x::T) where {T <: Union{Ptr,Core.LLVMPtr}}
    __vstore!(p, reinterpret(UInt, x), False(), False(), False(), register_size())
end
@inline function store!(p::Ptr{UInt}, x::T) where {T <: NativeTypes}
    __vstore!(Base.unsafe_convert(Ptr{T}, p), x, False(), False(), False(), register_size())
end
# @inline function store!(p::Ptr{UInt}, x::T) where {T <: Union{Ptr,Core.LLVMPtr}}
#     __vstore!(reinterpret(Core.LLVMPtr{UInt,0}, p), reinterpret(UInt, x), False(), False(), False(), register_size())
# end
# @inline function store!(p::Ptr{UInt}, x::T) where {T <: NativeTypes}
#     __vstore!(reinterpret(Core.LLVMPtr{T,0}, p), x, False(), False(), False(), register_size())
# end
@inline store!(p::Ptr{UInt}, x::T) where {T} = (unsafe_store!(Base.unsafe_convert(Ptr{T}, p), x); nothing)

@inline load(p::Ptr{UInt}, ::Type{StaticInt{N}}, i) where {N} = i, StaticInt{N}()
@inline store!(p::Ptr{UInt}, ::StaticInt, i) = i

@generated function load(p::Ptr{UInt}, ::Type{StridedPointer{T,N,C,B,R,X,O}}, i) where {T,N,C,B,R,X,O}
    q = quote
        $(Expr(:meta,:inline))
        i, ptr = load(p, Ptr{$T}, i)
    end
    xt = Expr(:tuple)
    Xp = X.parameters
    for n ∈ 1:N
        x = Symbol(:x_,n)
        push!(xt.args, x)
        push!(q.args, :((i, $x) = load(p, $(Xp[n]), i)))
    end
    ot = Expr(:tuple)
    Op = O.parameters
    for n ∈ 1:N
        o = Symbol(:o_,n)
        push!(ot.args, o)
        push!(q.args, :((i, $o) = load(p, $(Op[n]), i)))
    end
    push!(q.args, :((i, StridedPointer{$T,$N,$C,$B,$R}(ptr, $xt, $ot))))
    q
end
@generated function store!(p::Ptr{UInt}, ptr::StridedPointer{T,N,C,B,R,X,O}, i) where {T,N,C,B,R,X,O}
    q = quote
        $(Expr(:meta,:inline))
        i = store!(p, pointer(ptr), i)
        strd = strides(ptr)
        offs = offsets(ptr)
    end
    for n ∈ 1:N
        push!(q.args, :(i = store!(p, strd[$n], i)))
    end
    for n ∈ 1:N
        push!(q.args, :(i = store!(p, offs[$n], i)))
    end
    push!(q.args, :i)
    q
end

@inline function load(p::Ptr{UInt}, ::Type{T}, i) where {T}
    i + sizeof(T), load(p + i, T)
end
@inline function store!(p::Ptr{UInt}, x, i)
    store!(p + i, x)
    i + sizeof(x)
end

@generated function load(p::Ptr{UInt}, ::Type{T}, i) where {T<:Tuple}
    isbitstype(T) || return :(i + $(sizeof(T)), load(p + i, $T))
    q = Expr(:block, Expr(:meta,:inline))
    tup = Expr(:tuple)
    for (i,t) ∈ enumerate(T.parameters)
        ln = Symbol(:l_,i)
        push!(tup.args, ln)
        push!(q.args, :((i,$ln) = load(p, $t, i)))
    end
    push!(q.args, :(i, $tup))
    q
end
@generated store!(p::Ptr{UInt}, tup::T, i) where {T<:Tuple} = isbitstype(T) ? :(_store!(p, tup, i)) : :(store!(p+i, tup); i + sizeof(tup))
@inline function _store!(p::Ptr{UInt}, tup::Tuple{A,B,Vararg{Any,N}}, i) where {A,B,N}
    i = store!(p, first(tup), i)
    _store!(p, Base.tail(tup), i)
end
@inline function _store!(p::Ptr{UInt}, tup::Tuple{A}, i) where {A}
    store!(p, first(tup), i)
end
@inline store!(p::Ptr{UInt}, tup::Tuple{}, i) = i
@inline store!(p::Ptr{UInt}, tup::Nothing, i) = i


text:
# # [Advanced Recording Example](@id RecordingTutorial)
#
# The recording and debug possiblities make it possible to record nearly any data during the iterations.
# This tutorial illustrates how to
# * record one value during the iterations
# * record multiple values during the iterations and access them afterwards
# * define an own [`RecordAction`](@ref) to perform individual recordings.
#
# Several predefined recordings exist, for example [`RecordCost`](@ref) or [`RecordGradient`](@ref), depending on the solver used.
# For fields of the [`Options`](@ref) this can be directly done using the [`RecordEntry`](@ref).
# For others, an own [`RecordAction`](@ref) can be defined.
#
# We illustrate these using the gradient descent used in the introductionary [Get Started: Optimize!](@ref Optimize) tutorial example of computing the Riemannian Center of mass and refer to that tutorial for the mathematical details.
using Manopt, Manifolds, Random
Random.seed!(42)
m = 30
M = Sphere(m)
n = 800
σ = π / 8
x = zeros(Float64, m + 1)
x[2] = 1.0
data = [exp(M, x, random_tangent(M, x, Val(:Gaussian), σ)) for i in 1:n];
F(M, y) = sum(1 / (2 * n) * distance.(Ref(M), Ref(y), data) .^ 2)
gradF(M, y) = sum(1 / n * grad_distance.(Ref(M), data, Ref(y)))
nothing #hide
# ## Plain examples
#
# For the high level interfaces of the solvers, like [`gradient_descent`](@ref) we have to set `return_options` to `true` to obtain the whole options structure and not only the resulting resulting minimizer.
#
R = gradient_descent(M, F, gradF, data[1]; record=:Cost, return_options=true)
#
# You can attach different recorders to some operations (`:Start`. `:Stop`, `:Iteration` at time of writing), where `:Iteration` is the default, so the following is the same as `get_record(R, :Iteration)`. We get
#
get_record(R)
#
# To record more than one value, you can pass a array of a mix of symbols and [`RecordAction`](@ref) which gets mapped to a [`RecordGroup`](@ref)
R = gradient_descent(M, F, gradF, data[1]; record=[:Iteration, :Cost], return_options=true)
#
# Here, the Symbol `:Cost` is mapped to using the [`RecordCost`](@ref) action. The same holds for `:Iteration` and `:Iterate` and any member field of the current [`Options`](@ref).
# To access these you can first extract the group of records (of the `:Iteration` action) and then access the `:Cost`
ra = get_record_action(R)[:Cost]
#
# Or similarly
#
get_record(R, :Iteration, :Cost)
#
# Note that the first symbol again refers to the point where we record (not to the thing we record).
# We can also pass a Tuple as second argument to have our own order (not that now the second `:Iteration` refers to the recorded iteratons)
get_record(R, :Iteration, (:Cost, :Iteration))
#
# ## A more complex example
#
# To illustrate a complicated example let's record
# * the iteration number, cost and gradient field, but only every sixth iteration
# * the iteration at which we stop
#
# We first generate the problem and the options
#
p = GradientProblem(M, F, gradF)
o = GradientDescentOptions(
    M,
    copy(data[1]);
    stopping_criterion=StopAfterIteration(200) | StopWhenGradientNormLess(10.0^-9),
)
#
# and now decorate these with [`RecordOptions`](@ref)
#
rI = RecordEvery(
    RecordGroup([
        :Iteration => RecordIteration(),
        :Cost => RecordCost(),
        :Gradient => RecordEntry(similar(data[1]), :gradient),
    ]),
    6,
)
sI = RecordIteration()
r = RecordOptions(o, Dict(:Iteration => rI, :Stop => sI))
r2 = solve(p, r)
#
# and we see
#
get_record(r2, :Stop)
#
# as well as
#
get_record(r2, :Iteration, (:Iteration, :Cost))
#
# Here it is interesting to see, that a meta-record like [`RecordEvery`](@ref) just passes the tuple further on, so we can again also do
#
get_record_action(r2, :Iteration)[:Gradient]
#
# ## Writing an own [`RecordAction`](@ref)s
#
# Let's investigate where we want to count the number of function evaluations, again just to illustrate, since for the gradient this is just one evaluation per iteration.
# We first define a cost, that counts it's own calls.
#
mutable struct MyCost{T}
    data::T
    count::Int
end
MyCost(data::T) where {T} = MyCost{T}(data, 0)
function (c::MyCost)(M, x)
    c.count += 1
    return sum(1 / (2 * length(c.data)) * distance.(Ref(M), Ref(x), c.data) .^ 2)
end
nothing #hide
#
# and we define the following RecordAction,
#
mutable struct RecordCount <: RecordAction
    recorded_values::Vector{Int}
    RecordCount() = new(Vector{Int}())
end
function (r::RecordCount)(p::Problem, ::Options, i)
    if i > 0
        push!(r.recorded_values, p.cost.count)
    elseif i < 0 # reset if negative
        r.recorded_values = Vector{Int}()
    end
end
nothing #hide
#
# And now we can initialize the new cost and call the gradient descent.
# Note that this illustrates also the last use case – you can pass symbol-Action pairs into the `record=`array.
#
F2 = MyCost(data)
R = gradient_descent(
    M,
    F2,
    gradF,
    data[1];
    record=[:Iteration, :Count => RecordCount(), :Cost],
    return_options=true,
)
#
# We can again access the whole sets of records
#
get_record(R)
#
# this is equivalent to calling `R[:Iteration`.
# Note that since we introduced `:Count` we can also access a single recorded value using
#
R[:Iteration, :Count]


text:
using XEB
using Test

@testset "XEB.jl" begin
    # Write your tests here.
end


text:
@testset "Inflation with MixedCountryStructure" begin

    GB = 200
    PER = 120
    baseindex = rand(100:0.5:110, GB)

    v = rand(PER, GB) .- 0.25
    w = rand(GB)
    w = w / sum(w)

    basedate = Date(2001,1)
    dates = getdates(basedate, PER)

    # Base with different base indexes per product
    vcpi_mixed = VarCPIBase(v, w, dates, baseindex)
    mcs = MixedCountryStructure(vcpi_mixed)

    # Basic inflation function test
    totalfn = InflationTotalCPI()
    @test totalfn(mcs) isa Vector
    @test infl_periods(mcs) == (PER-11)

    # Mix with base with same base index
    GB = 315
    PER = 120
    baseindex = 100.0

    v2 = rand(PER, GB) .- 0.25
    w2 = rand(GB)
    w2 = w2 / sum(w)

    basedate = Date(2011,1)
    dates2 = getdates(basedate, PER)

    vcpi_uniform = VarCPIBase(v2, w2, dates2, baseindex)
    mcs = MixedCountryStructure(vcpi_mixed, vcpi_uniform)

    @test totalfn(mcs) isa Vector
    @test infl_periods(mcs) == (240-11)

    # Test indexing
    @test mcs[1] == vcpi_mixed
    @test mcs[2] == vcpi_uniform

end


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

using Markdown
using InteractiveUtils

# ╔═╡ aed9a998-1ee8-11eb-1d55-31076d75bfb6
using Pkg, DrWatson

# ╔═╡ a818b7b6-1ee8-11eb-167a-f3d0a4cb736a
md"## Listing 2.6"

# ╔═╡ 8fce84c2-1ee5-11eb-38fc-9f2c2720cb71
begin
	@quickactivate "StatisticsWithJuliaPlutoNotebooks"
end

# ╔═╡ b7052e7e-1ee8-11eb-3f69-af0a0c82375f
begin
	A = Set([2,7,2,3])
	B = Set(1:6)
	omega = Set(1:10)
end;

# ╔═╡ 49ce4878-1ee9-11eb-03bd-59a4e141e22c
begin
	AunionB = union(A, B)
	AintersectionB = intersect(A, B)
	BdifferenceA = setdiff(B,A)
	Bcomplement = setdiff(omega,B)
	AsymDifferenceB = union(setdiff(A,B),setdiff(B,A))
end;

# ╔═╡ 5fa63066-1ee9-11eb-033f-79d6fc16652a
("A = $A, B = $B")

# ╔═╡ 49ce9e86-1ee9-11eb-1371-37dac8fbaeab
("A union B = $AunionB")

# ╔═╡ 49cf15dc-1ee9-11eb-31ba-eb7bd106c7a8
("A intersection B = $AintersectionB")

# ╔═╡ 6cc4003e-1ee9-11eb-259d-39c42a956b07
("B diff A = $BdifferenceA")

# ╔═╡ 49db80e0-1ee9-11eb-34af-515e8414ea68
("B complement = $Bcomplement")

# ╔═╡ 74f756e8-1ee9-11eb-1edb-53c58d2386d1
("A symDifference B = $AsymDifferenceB")

# ╔═╡ 49dc3514-1ee9-11eb-33af-01a3ab419dee
("The element '6' is an element of A: $(in(6,A))")

# ╔═╡ 7c6cc9da-1ee9-11eb-337f-c1e274f93d10
("Symmetric difference and intersection are subsets of the union: ",
        issubset(AsymDifferenceB,AunionB),", ", issubset(AintersectionB,AunionB))

# ╔═╡ 8c95b6f0-1ee9-11eb-0fd8-cd7a30b3ca3c
md"## End of listing 2.6"

# ╔═╡ Cell order:
# ╟─a818b7b6-1ee8-11eb-167a-f3d0a4cb736a
# ╠═aed9a998-1ee8-11eb-1d55-31076d75bfb6
# ╠═8fce84c2-1ee5-11eb-38fc-9f2c2720cb71
# ╠═b7052e7e-1ee8-11eb-3f69-af0a0c82375f
# ╠═49ce4878-1ee9-11eb-03bd-59a4e141e22c
# ╠═5fa63066-1ee9-11eb-033f-79d6fc16652a
# ╠═49ce9e86-1ee9-11eb-1371-37dac8fbaeab
# ╠═49cf15dc-1ee9-11eb-31ba-eb7bd106c7a8
# ╠═6cc4003e-1ee9-11eb-259d-39c42a956b07
# ╠═49db80e0-1ee9-11eb-34af-515e8414ea68
# ╠═74f756e8-1ee9-11eb-1edb-53c58d2386d1
# ╠═49dc3514-1ee9-11eb-33af-01a3ab419dee
# ╠═7c6cc9da-1ee9-11eb-337f-c1e274f93d10
# ╟─8c95b6f0-1ee9-11eb-0fd8-cd7a30b3ca3c


text:
using StatsFuns

function log_likelihood_CIR(δ, γ, σ, λ, data)
    times = data |> keys |> collect |> sort

    α = δ/2#Alternative parametrisation

    res = Dict()

    res[times[1]] = logμπh_param_δγ(0, γ/σ^2, δ, data[times[1]])

    #Prior
    Λ_prime = [0]
    θ_prime = γ/σ^2

    #1st update
    Λ = next_Λ_from_Λ_prime(Λ_prime, data[times[1]], t_CIR)
    θ = θ_from_θ_prime(data[times[1]], θ_prime, T_CIR)
    log_wms = [0]

    #1st prediction
    Λ_prime = Λ_prime_1D(Λ)
    log_wms_prime = next_log_wms_prime_from_log_wms1D(log_wms, Λ, times[2]-times[1], θ, γ, σ)
    θ_prime = θ_prime_from_θ_CIR(θ, times[2]-times[1], γ, σ)

    res[times[2]] = res[times[1]] + logμν_i_minus_1(Λ_prime, log_wms_prime, θ_prime, data[times[2]], δ)

    for i in 2:(length(times)-1)
        t = times[i]
        next_t = times[i+1]
        Δt = next_t-t

        #Update
        Λ = next_Λ_from_Λ_prime(Λ_prime, data[t], t_CIR)
        θ = θ_from_θ_prime(data[t], θ_prime, T_CIR)
        log_wms = next_logwms_from_log_wms_prime1D(log_wms_prime, Λ_prime, data[t], θ_prime, α)

         #Prediction
        Λ_prime = Λ_prime_1D(Λ)
        θ_prime = θ_prime_from_θ_CIR(θ, Δt, γ, σ)
        log_wms_prime = next_log_wms_prime_from_log_wms1D(log_wms, Λ, Δt, θ, γ, σ)

        res[next_t] = res[t] + logμν_i_minus_1(Λ_prime, log_wms_prime, θ_prime, data[next_t], δ)

    end

    return res
end

logμπh_param_δγ(m, θ, δ, y; λ = 1) = logμπh(m, θ, δ/2, y; λ = λ)

function log_likelihood_CIR_keep_fixed_number(δ, γ, σ, λ, data, fixed_number::Int64; silence = false)
    # println("filter_WF_mem2")

    prune_keeping_fixed_number_log_wms_(Λ_of_t, log_wms_of_t) =
        prune_keeping_fixed_number_log_wms(Λ_of_t, log_wms_of_t, fixed_number)

    log_likelihood_pruning(δ, γ, σ, λ, data, prune_keeping_fixed_number_log_wms_; silence = silence)

end

function log_likelihood_CIR_keep_above_threshold(δ, γ, σ, λ, data, ε::Float64; silence = false)
    # println("filter_WF_mem2")

    prune_keeping_above_threshold_log_wms_(Λ_of_t, log_wms_of_t) = prune_keeping_above_threshold_log_wms(Λ_of_t, log_wms_of_t, ε)

    log_likelihood_pruning(δ, γ, σ, λ, data, prune_keeping_above_threshold_log_wms_; silence = silence)

end

function log_likelihood_CIR_fixed_fraction(δ, γ, σ, λ, data, fraction::Float64; silence = false)
    # println("filter_WF_mem2")
    prune_keeping_fixed_fraction_log_wms_(Λ_of_t, log_wms_of_t) = prune_keeping_fixed_fraction_log_wms(Λ_of_t, log_wms_of_t, fraction)

    log_likelihood_pruning(δ, γ, σ, λ, data, prune_keeping_fixed_fraction_log_wms_; silence = silence)

end

function log_likelihood_pruning(δ, γ, σ, λ, data, do_the_pruning_log_wms::Function; silence = false)
    times = data |> keys |> collect |> sort

    α = δ/2#Alternative parametrisation

    res = Dict()

    res[times[1]] = logμπh_param_δγ(0, γ/σ^2, δ, data[times[1]])

    #Prior
    Λ_prime = [0]
    θ_prime = γ/σ^2

    #1st update
    Λ = next_Λ_from_Λ_prime(Λ_prime, data[times[1]], t_CIR)
    θ = θ_from_θ_prime(data[times[1]], θ_prime, T_CIR)
    log_wms = [0]

    pruned_Λ, pruned_log_wms = do_the_pruning_log_wms(Λ, log_wms)

    #1st prediction
    Λ_prime = Λ_prime_1D(pruned_Λ)
    log_wms_prime = next_log_wms_prime_from_log_wms1D(pruned_log_wms, pruned_Λ, times[2]-times[1], θ, γ, σ)
    θ_prime = θ_prime_from_θ_CIR(θ, times[2]-times[1], γ, σ)

    res[times[2]] = res[times[1]] + logμν_i_minus_1(Λ_prime, log_wms_prime, θ_prime, data[times[2]], δ)

    for i in 2:(length(times)-1)
        t = times[i]
        next_t = times[i+1]
        Δt = next_t-t

        #Update
        Λ = next_Λ_from_Λ_prime(Λ_prime, data[t], t_CIR)
        θ = θ_from_θ_prime(data[t], θ_prime, T_CIR)
        log_wms = next_logwms_from_log_wms_prime1D(log_wms_prime, Λ_prime, data[t], θ_prime, α)

        #Pruning
        pruned_Λ, pruned_log_wms = do_the_pruning_log_wms(Λ, log_wms)

         #Prediction
        Λ_prime = Λ_prime_1D(pruned_Λ)
        θ_prime = θ_prime_from_θ_CIR(θ, Δt, γ, σ)
        log_wms_prime = next_log_wms_prime_from_log_wms1D(pruned_log_wms, pruned_Λ, Δt, θ, γ, σ)

        res[next_t] = res[t] + logμν_i_minus_1(Λ_prime, log_wms_prime, θ_prime, data[next_t], δ)

    end

    return res
end

function t_CIR(y::Number, m::Integer)
    return m + y
end
function t_CIR(y::Number, m::Array{U, 1}) where U<:Integer
    return m .+ y
end
function t_CIR(y::Array{T, 1}, m::Integer) where T<:Number
    return m + sum(y)
end
function t_CIR_inside(y::Array{T, 1}, m) where T<:Number
    return m .+ sum(y)
end
function t_CIR(y::Array{T, 1}, m::Array{U, 1}) where {T<:Number, U<:Integer}
    return t_CIR_inside(y::Array{T, 1}, m)
end
function t_CIR(y::Array{T, 1}, m::UnitRange{Int64}) where T<:Number
    return t_CIR_inside(y::Array{T, 1}, m)
end

function Λ_prime_1D(Λ::Union{Array{T,1}, UnitRange{Int64}}) where T<:Integer
    return 0:maximum(Λ)
end

function next_Λ_from_Λ_prime(Λ_prime, y, t)
    return t(y, Λ_prime)
end

function θ_prime_from_θ_CIR(θ, Δt, γ, σ)
    return γ/σ^2 * θ*exp(2*γ*Δt)/(θ*exp(2*γ*Δt)+γ/σ^2-θ)
end

function θ_from_θ_prime(y, θ_prime, T)
    return T(y, θ_prime)
end

T_CIR(y, θ; λ = 1) = θ + λ*length(y)

function next_logwms_from_log_wms_prime1D(log_wms_prime, Λ_prime, y, θ_prime, α)#update
    #Make sure we deal correctly with weights equal to 0
    #Probably harmonise what can be improved in the filtering algorithm
    unnormalised_log_wms = log_wms_prime .+ logμπh(Λ_prime, θ_prime, α, y)
    return unnormalised_log_wms .- StatsFuns.logsumexp(unnormalised_log_wms)
end

function next_log_wms_prime_from_log_wms1D(log_wms, Λ, Δt, θ, γ, σ)
    #The speed gain seems like a modest 20%, but much less allocations though.
    nΛ = length(Λ)
    maxΛ = maximum(Λ)
    wms_prime = zeros(maxΛ+1)
    all_u = zeros(maxΛ+1)
    all_u .= -Inf
    p = γ/σ^2*1/(θ*exp(2*γ*Δt) + γ/σ^2 - θ)
    @inbounds for k in eachindex(Λ)
        m = Λ[k]
        for n in 0:m
            idx = n+1

            x = log_wms[k] + logpdf(Binomial(m, p), n)
            if x <= all_u[idx]
                wms_prime[idx] += exp(x - all_u[idx])
            else
                wms_prime[idx] *= exp(all_u[idx] - x)
                wms_prime[idx] += 1.0
                all_u[idx] = x
            end
        end
    end
    # return @. log(wms_prime) + all_u
    return log.(wms_prime) .+ all_u
end

function logμν_i_minus_1(Λ_prime_i_minus_1, log_wms_prime_i_minus_1, θ_prime_i_minus_1, yi, δ)
    StatsFuns.logsumexp(log_wms_prime_i_minus_1 .+ logμπh_param_δγ(Λ_prime_i_minus_1, θ_prime_i_minus_1, δ, yi))
end

max_stars_repo_name: konkam/DualOptimalFiltering

max_issues_repo_name: konkam/DualOptimalFiltering

max_forks_repo_name: konkam/DualOptimalFiltering


text:
"""
`CircStats` contains routines for investigating the statistics of circular
data.

By default, all input and output from routines is in radians, but in general passing
`true` as the last argument to a routine will change this to degrees.
"""
module CircStats

using Printf

import Distributions

using StatsBase: mean

export
    # Summary stats
    cdist,
    cmean,
    cmedian,
    cresultant,
    cstd,
    cvariance,
    # Distributions
    von_mises_cdf,
    von_mises_pdf,
    # Fitting/estimation
    fit_vonMises,
    # Hypothesis testing
    V_test,
    watson_U2,
    watson_U2n

include("Datasets.jl")
using .Datasets

include("hypothesistests.jl")

"""
    cdist(a, b, degrees::Bool=false) -> angle

Return the angular distance from `a` to `b` (b - a) in the forward
direction; hence if `b` is 'behind' `a`, `distance` will be negative.

Angles are confined to be the smaller possible angle, so are in the range
[-π:π], or [-180°:180°].

Angles are in radians, unless `degrees` == true.
"""
cdist(a, b) = mod(b - a + pi, 2π) - pi
cdist(a, b, degrees::Bool) =
    degrees ? rad2deg(cdist(deg2rad(a), deg2rad(b))) : cdist(a, b)

"""
    cmean(a::Array, degrees::Bool=false) -> mean

Return the mean angle from the set of angles in `a`.

Angles are in radians, unless `degrees` == true.
"""
cmean(a) = atan(sum(sin.(a)), sum(cos.(a)))
cmean(a, degrees::Bool) = degrees ? rad2deg(cmean(deg2rad.(a))) : cmean(a)

"""
    cmedian(a::Array, degrees::Bool=false, axial=false) -> median

Return the median angle from the set of angles `a`.

Angles are in radians, unless `degrees` == true.
"""
function cmedian(a, degrees::Bool=false, axial::Bool=false)
    medians = Vector{Float64}()
    A = Float64.(sort(a))
    n = length(A)
    degrees && (A .= deg2rad.(A))
    axial && (A .= 2A)
    # Compute trial bisectors, which are either the points themselves or adjacent means
    p = Array{Float64}(undef, n)
    if iseven(n)
        for i = 1:n
            j = (i + 1 <= n) ? i + 1 : 1
            p[i] = cmean([A[i], A[j]])
        end
    else
        p[:] .= A[:]
    end
    # Try all possible diameters
    for i = 1:n
        # Count points on either side of diameter
        n_plus = sum(cdist.(p[i], A) .> 0)
        n_minus = sum(cdist.(p[i], A) .< 0)
        if n_plus == n_minus
            # Determine which side of circle is correct direction by counting the number
            # within pi/2 of each of the two opposing possible medians
            if sum(abs.(cdist.(p[i], A)) .<= pi/2) > sum(abs.(cdist.(p[i], A)) .> pi/2)
                push!(medians, A[i])
            else
                push!(medians, (A[i] + 2pi)%2pi - pi)
            end
        end
    end
    # If there is more than one median, take the mean thereof (Otieno & Anderson-Cook
    # (2003), Journal of Modern Applied Statistical Methods, 2(1), 168-176)
    median = if length(medians) > 1
        cmean(medians)
    elseif length(medians) == 1
        medians[1]
    else
        error("Zero medians found.  Are data axial but axial!=true?")
    end
    degrees && (median = rad2deg(median))
    axial && (median /= 2)
    median
end

"""
    cresultant(a, degrees=false) -> R
    cresultant(a, w, degrees=false) -> Rc

Return the resultant vector length, `R`, from a set of angles, `a`.

If data are binned by binwidth `w`, an unbiased estimate of `R`, `Rc` is returned
when `w` is supplied.

Angles are in radians, unless `degrees` == true.
"""
cresultant(a, degrees::Bool=false) = degrees ?
    sqrt(sum(sin.(deg2rad.(a)))^2 + sum(cos.(deg2rad.(a)))^2)/length(a) :
    sqrt(sum(sin.(a))^2 + sum(cos.(a))^2)/length(a)
cresultant(a, w::Real, degrees::Bool=false) = degrees ?
    cresultant(a, true)*deg2rad(w)/(2sin(deg2rad(w)/2)) : cresultant(a)*w/(2sin(w/2))

"""
    cstd(a, degrees=false) -> σ

Return the standard deviation, `σ`, for a set of angles `a`.

Angles are in radians, unless `degrees == true`.
"""
cstd(a, degrees::Bool=false) = sqrt(-2*log(cresultant(a, degrees)))

"""
    cvariance(a, degrees=false) -> σ²

Return the circular variance, `σ²`, of a set of angles `a`.

Angles are in radians, unless `degrees == true`.
"""
function cvariance(a, degrees::Bool=false)
    a_mean = cmean(a, degrees)
    if degrees
        1 - sum(cos.(deg2rad.(a .- a_mean)))/length(a)
    else
        1 - sum(cos.(a .- a_mean))/length(a)
    end
end

"""
    von_mises_pdf(a, µ, κ, degrees=false) -> p

Return the Von Mises probability density function at `a`, for a von Mises
distribution with circular mean `µ` and concentration `κ`.

Angles are in radians, unless `degrees` == true.
"""
von_mises_pdf(a, mu::Real, k::Real, degrees::Bool=false) = degrees ?
                                    von_mises_pdf(deg2rad(a), deg2rad(mu), k) :
                                    Distributions.pdf(Distributions.VonMises(mu, k), a)

"""
    von_mises_cdf(a, μ, κ, degrees=false) -> cdf

Return the von Mises cumulative distribution function at
`a` for a distribution with mean `μ` and concentration `κ`.

Angles are in radians, unless `degrees` == true.
"""
function von_mises_cdf(a, μ, κ, degrees=false)
    degrees && (μ = deg2rad(μ))
    d = Distributions.VonMises(μ, κ)
    Distributions.cdf(d, degrees ? deg2rad(a) : a)
end

"""
    fit_vonMises(θ, degrees=false; axial=false) -> μ, κ

Return the maximum-likelihood values of the mean `μ` and concentration
`κ` for the von Mises distribution which fits the set of angles `θ`.

If `axial` is `true`, then the data are assumed to have a π (or 180°)
ambiguity.

Angles are in radians, unless `degrees` == true.
"""
function fit_vonMises(θ, degrees=false; axial=false)
    degrees && (θ = deg2rad.(θ))
    axial && (θ = 2θ)
    μ = cmean(θ)
    κ = estimate_kappa_vonMises(θ, μ)
    degrees ? (rad2deg(μ), κ) : (μ, κ)
end

"""
    estimate_kappa_vonMises(θ, μ) -> κ

Return the maximum likelihood estimation of the concentration
parameter `κ` for the best-fitting von Mises distribution, given
a set of angles `θ` in radians.

Uses the polynomial approximation given by Best and Fisher (1981).
**N.B.** This may not be reliable when R̄ is small (e.g., < 0.7).
"""
function estimate_kappa_vonMises(θ, μ)
    R̄ = mean(cos.(θ .- μ))
    if 0 <= R̄ < 0.53
        2R̄ + R̄^3 + 5R̄^5/6
    elseif 0.53 <= R̄ < 0.85
        -0.4 + 1.39R̄ + 0.43/(1 - R̄)
    elseif R̄ >= 0.85
        1/(R̄^3 - 4R̄^2 + 3R̄)
    end
end

end # module


text:
################################################################################
# Stamps
################################################################################

################################################################################
# Type Definition
################################################################################

abstract type AbstractStamp
end

mutable struct Stamp <: AbstractStamp
	prob::Symbol # name of the problem
	i0::Int      # index of the problem
	n1::Symbol   # name of the variable
	i1::Int      # index of the variable
	v1::Int      # value of the variable
	n2::Symbol   # name of the variable
	i2::Int      # index of the variable
	v2::Int      # value of the variable
end

mutable struct VStamp <: AbstractStamp
	prob::Symbol # name of the problem
	i0::Int      # index of the problem
	n1::Symbol   # name of the variable
	i1::Int      # index of the variable
	v1::Int      # value of the variable
end

mutable struct HStamp <: AbstractStamp
	n2::Symbol   # name of the variable
	i2::Int      # index of the variable
	v2::Int      # value of the variable
end

################################################################################
# Base functions
################################################################################

import Base.hash
import Base.==
import Base.isequal

Base.hash(stamp::Stamp, h::UInt)  = hash(stamp.prob, hash(stamp.i0, hash(stamp.n1, hash(stamp.i1, hash(stamp.v1, hash(stamp.n2, hash(stamp.i2, hash(stamp.v2, h))))))))
Base.hash(stamp::VStamp, h::UInt) = hash(stamp.prob, hash(stamp.i0, hash(stamp.n1, hash(stamp.i1, hash(stamp.v1, h)))))
Base.hash(stamp::HStamp, h::UInt) = hash(stamp.n2, hash(stamp.i2, hash(stamp.v2, h)))

function (==)(stamp1::AbstractStamp, stamp2::AbstractStamp)
    out = true
    T1 = typeof(stamp1)
    T2 = typeof(stamp2)
    T1 == T2 ? nothing : return false
    for name in fieldnames(T1)
        out &= getfield(stamp1, name) == getfield(stamp2, name)
    end
    return out
end

function (==)(stamp1::Stamp, stamp2::Stamp)
    out = true
	for name in Stamp.name.names
        out &= getfield(stamp1, name) == getfield(stamp2, name)
    end
    return out
end

function (==)(stamp1::VStamp, stamp2::VStamp)
    out = true
	for name in VStamp.name.names
        out &= getfield(stamp1, name) == getfield(stamp2, name)
    end
    return out
end

function (==)(stamp1::HStamp, stamp2::HStamp)
    out = true
	for name in HStamp.name.names
        out &= getfield(stamp1, name) == getfield(stamp2, name)
    end
    return out
end

Base.isequal(stamp1::Stamp, stamp2::Stamp) = stamp1 == stamp2
Base.isequal(stamp1::VStamp, stamp2::VStamp) = stamp1 == stamp2
Base.isequal(stamp1::HStamp, stamp2::HStamp) = stamp1 == stamp2

################################################################################
# Methods
################################################################################

function Stamp()
	return Stamp(:x, 0, :x, 0, 0, :x, 0, 0)
end

function VStamp()
	return VStamp(:x, 0, :x, 0, 0)
end

function HStamp()
	return HStamp(:x, 0, 0)
end

function stampify(name::Symbol, step::Int)
	if name == :x_1
		name = :x
		step -= 1
	elseif name == :x1
		name = :x
		step += 1
	end
	return name, step
end

function stampify(prob::Symbol, ind_p::Int, name_i::Symbol, ind_i::Int, step_i::Int, name_j::Symbol, ind_j::Int, step_j::Int)
	name_i, step_i = stampify(name_i, step_i)
	name_j, step_j = stampify(name_j, step_j)
	return Stamp(prob, ind_p, name_i, ind_i, step_i, name_j, ind_j, step_j)
end

function stampify!(stamp::Stamp, prob::Symbol, ind_p::Int, name_i::Symbol, ind_i::Int, step_i::Int, name_j::Symbol, ind_j::Int, step_j::Int)
	name_i, step_i = stampify(name_i, step_i)
	name_j, step_j = stampify(name_j, step_j)
	stamp.prob = prob
	stamp.i0 = ind_p
	stamp.n1 = name_i
	stamp.i1 = ind_i
	stamp.v1 = step_i
	stamp.n2 = name_j
	stamp.i2 = ind_j
	stamp.v2 = step_j
	return nothing
end

function stampify(prob::Symbol, ind_p::Int, name_i::Symbol, ind_i::Int, step_i::Int)
	name_i, step_i = stampify(name_i, step_i)
	return VStamp(prob, ind_p, name_i, ind_i, step_i)
end

function stampify!(stamp::VStamp, prob::Symbol, ind_p::Int, name_i::Symbol, ind_i::Int, step_i::Int)
	name_i, step_i = stampify(name_i, step_i)
	stamp.prob = prob
	stamp.i0 = ind_p
	stamp.n1 = name_i
	stamp.i1 = ind_i
	stamp.v1 = step_i
	return nothing
end

function stampify(name_j::Symbol, ind_j::Int, step_j::Int)
	name_j, step_j = stampify(name_j, step_j)
	return HStamp(name_j, ind_j, step_j)
end

function stampify!(stamp::HStamp, name_j::Symbol, ind_j::Int, step_j::Int)
	name_j, step_j = stampify(name_j, step_j)
	stamp.n2 = name_j
	stamp.i2 = ind_j
	stamp.v2 = step_j
	return nothing
end

################################################################################
# Stamp Validity
################################################################################

function valid(s::Stamp, N::Int, p::Int)
	return valid(s.prob, s.i0, s.n1, s.i1, s.v1, s.n2, s.i2, s.v2, N, p)
end

function valid(prob::Symbol, i0::Int, n1::Symbol, i1::Int, v1::Int, n2::Symbol, i2::Int, v2::Int, N::Int, p::Int)
	b1 = valid(prob, i0, n1, i1, v1, N, p) # stamp 1 is valid
	b2 = false # stamp 2 is valid

	if prob == :opt && i0 ∈ (1:p)
		if n2 == :u && i2 ∈ (1:p) && v2 ∈ (1:N-1) # uj1, ...N-1
			b2 = true
		elseif n2 == :λ && i2 == i0 && v2 ∈ (1:N-1) # λi1 ...N-1
			b2 = true
		elseif n2 == :x && i2 == 1 && v2 ∈ (2:N) # x2...xN
			b2 = true
		end
	end
	if prob == :dyn && i0 == 1
		if n2 == :u && i2 ∈ (1:p) && v2 ∈ (1:N-1) # uj1, ...N-1
			b2 = true
		elseif n2 == :x && i2 == 1 && v2 ∈ (2:N) # x2...xN
			b2 = true
		end
	end
	return b1 && b2
end

function valid(s::VStamp, N::Int, p::Int)
	return valid(s.prob, s.i0, s.n1, s.i1, s.v1, N, p)
end

function valid(prob::Symbol, i0::Int, n1::Symbol, i1::Int, v1::Int, N::Int, p::Int)
	b1 = false # stamp 1 is valid
	if prob == :opt && i0 ∈ (1:p)
		if n1 == :u && i1 == i0 && v1 ∈ (1:N-1) # u1...uN-1
			b1 = true
		elseif n1 == :x && i1 == 1 && v1 ∈ (2:N) # x2...xN
			b1 = true
		end
	end
	if prob == :dyn && i0 == 1
		if n1 == :x && i1 == 1 && v1 ∈ (1:N-1) # x1...xN-1
			b1 = true
		end
	end
	return b1
end

function valid(s::HStamp, N::Int, p::Int)
	return valid(s.n2, s.i2, s.v2, N, p)
end

function valid(n2::Symbol, i2::Int, v2::Int, N::Int, p::Int)
	b2 = false # stamp 2 is valid
	if n2 == :u && i2 ∈ (1:p) && v2 ∈ (1:N-1) # uj1, ...N-1
		b2 = true
	elseif n2 == :λ && i2 ∈ (1:p) && v2 ∈ (1:N-1) # λi1 ...N-1
		b2 = true
	elseif n2 == :x && i2 == 1 && v2 ∈ (2:N) # x2...xN
		b2 = true
	end
	return b2
end


text:
using VLModelParametersDB
using JSON
using CSV

# load the default JSON file -
path_to_json_file = "/Users/jeffreyvarner/Desktop/julia_work/VLModelParametersDB.jl/src/database/data/default_human.json"
biophysical_json_dict = JSON.parsefile(path_to_json_file)["biophysical_constants"]

# convert -
df = transform(biophysical_json_dict, "HUMAN")

# write -
path_to_csv_file = "/Users/jeffreyvarner/Desktop/julia_work/VLModelParametersDB.jl/src/database/data/human.csv"
CSV.write(path_to_csv_file, df)


text:
push!(LOAD_PATH, "/mnt/g/Python Programme/BALL.jl/src/" )
using BALL
using BioStructures
using DataFrames
using BenchmarkTools

struc = System("1ake.pdb")
println(struc)

#profilehtml
#= Buggy code: Dataframes hashing doesnt work
internal = System("1,2-Benzodiazepine.pdb")
dfs = DataFrameSystem(internal)

graph = MolecularGraph(dfs)
=#

internal = System("/mnt/g/Python Programme/BALL.jl/1ake.pdb")
internal2 = deepcopy(internal)
dfs = DataFrameSystem(internal)
#ramachandranangles(dfs)
dfs2 = deepcopy(dfs)
struc = read("/mnt/g/Python Programme/BALL.jl/1ake.pdb", PDB)
struc2 = deepcopy(struc)

for sys in struc2 for ch in sys for res in ch for at in res
	isdisorderedatom(at) && continue
    coords!(at, at.coords .+= 1/sqrt(3))
end end end end

for row in eachrow(dfs2.atoms)
    row.position_ = getPosition(row) .+ 1/sqrt(3)
end

for at in AtomIterator(internal2)
    at.position_ = getPosition(at) .+ 1/sqrt(3)
end

println(BALL.rmsd(dfs,dfs2))
println(BioStructures.rmsd(struc,struc2,superimpose=false))
println(BALL.rmsd(internal, internal2))

@show @benchmark BALL.rmsd(dfs,dfs2)
@show @benchmark BioStructures.rmsd(struc,struc2,superimpose=false)
@show @benchmark BALL.rmsd(internal, internal2)

@show @benchmark BALL.rmsd(dfs,dfs2, calphaonly=false)
@show @benchmark BioStructures.rmsd(struc,struc2,allselector,superimpose=false)
@show @benchmark BALL.rmsd(internal, internal2,false)


text:
# Table 2:
# 2014 Dec 18 57,009.13 Gemini/GPI H 14.161 ± 0.021   0.17 ± 0.14    454.24 ± 1.88   171.22 ± 0.23
# 2015 Jan 30 57,052.06 Gemini/GPI J 14.161 ± 0.021   0.17 ± 0.14    451.81 ± 2.06   170.01 ± 0.26
# 2015 Jan 31 57,053.06 Gemini/GPI H 14.161 ± 0.021   0.17 ± 0.14    456.80 ± 2.57   170.19 ± 0.30
# 2015 Feb 1 57,054.25 Keck/NIRC2  L′ 9.952 ± 0.002 −0.252 ± 0.009    461.5 ± 23.9   170.4  ± 3.0
# 2015 Sep 1 57,266.41 Gemini/GPI H  14.161 ± 0.021   0.17 ± 0.14    455.10 ± 2.23   167.30 ± 0.26
# 2015 Nov 6 57,332.23 Gemini/GPI K1 14.161 ± 0.021   0.21 ± 0.23    452.88 ± 5.41   166.12 ± 0.57
# 2015 Dec 18 57,374.19 Gemini/GPI 2 14.161 ± 0.021   0.21 ± 0.23    455.91 ± 6.23   165.66 ± 0.57
# 2015 Dec 20 57,376.17 Gemini/GPI H 14.161 ± 0.021   0.21 ± 0.23    455.01 ± 3.03   165.69 ± 0.43
# 2016 Jan 28 57,415.05 Gemini/GPIK1 14.161 ± 0.021   0.21 ± 0.23    454.46 ± 6.03   165.94 ± 0.51
# 2016 Sep 18 57,649.39 Gemini/GPI H 14.161 ± 0.021   0.32 ± 0.15    454.81 ± 2.02   161.80 ± 0.26
# 2016 Sep 21 57,652.38 Gemini/GPI J 14.161 ± 0.021   0.32 ± 0.15    451.43 ± 2.67   161.73 ± 0.31
# 2016 Dec 17 57,739.13 Gemini/GPI J 14.161 ± 0.021   0.32 ± 0.15    449.39 ± 2.15   160.06 ± 0.27
# 2017 Nov 11 58,068.26 Gemini/GPI H 14.161 ± 0.021   0.28 ± 0.19    447.54 ± 3.02   155.23 ± 0.39
# 2018 Nov 20 58,442.21 Gemini/GPI H 14.161 ± 0.021   0.45 ± 0.11    434.22 ± 2.01   149.64 ± 0.23
using Measurements
θ = [
    171.22 ± 0.23,
    170.01 ± 0.26,
    170.19 ± 0.30,
    170.4  ± 3.0,
    167.30 ± 0.26,
    166.12 ± 0.57,
    165.66 ± 0.57,
    165.69 ± 0.43,
    165.94 ± 0.51,
    161.80 ± 0.26,
    161.73 ± 0.31,
    160.06 ± 0.27,
    155.23 ± 0.39,
    149.64 ± 0.23,
]
ρ = [
    454.24 ± 1.88,
    451.81 ± 2.06,
    456.80 ± 2.57,
    461.5 ± 23.9 ,
    455.10 ± 2.23,
    452.88 ± 5.41,
    455.91 ± 6.23,
    455.01 ± 3.03,
    454.46 ± 6.03,
    454.81 ± 2.02,
    451.43 ± 2.67,
    449.39 ± 2.15,
    447.54 ± 3.02,
    434.22 ± 2.01,
]

dates = [
    "2014-12-18",
    "2015-01-30",
    "2015-01-31",
    "2015-02-01",
    "2015-09-01",
    "2015-11-06",
    "2015-12-18",
    "2015-12-20",
    "2016-01-28",
    "2016-09-18",
    "2016-09-21",
    "2016-12-17",
    "2017-11-11",
    "2018-11-20",
]
mjds = mjd.(dates)

ra1  = ρ .* sind.(θ)
dec1 = ρ .* cosd.(θ)

##
ra = Measurements.value.(ra1)
dec = Measurements.value.(dec1)
σ_ra = Measurements.uncertainty.(ra1)
σ_dec = Measurements.uncertainty.(dec1)

hcat(mjds, ra, dec, σ_ra, σ_dec)
##
# 57009.0  -448.917   69.3356   1.8787    1.82477
# 57052.0  -444.96    78.3783   2.05971   2.05054
# 57053.0  -450.121   77.8303   2.565     2.39716
# 57054.0  -455.037   76.9638  23.9074   24.1568
# 57266.0  -443.966  100.052    2.22231   2.07345
# 57332.0  -439.656  108.641    5.36208   4.56234
# 57374.0  -441.705  112.918    6.13954   4.65729
# 57376.0  -440.892  112.464    3.0549    3.39255
# 57415.0  -440.845  110.406    5.93133   4.18857
# 57649.0  -432.057  142.053    2.02432   2.05962
# 57652.0  -428.673  141.521    2.6485    2.46576
# 57739.0  -422.449  153.258    2.14627   2.12148
# 58068.0  -406.365  187.509    3.02463   3.04171
# 58442.0  -374.674  219.468    1.9453    1.815

max_stars_repo_name: sefffal/teps-problems

max_issues_repo_name: sefffal/teps-problems

max_forks_repo_name: sefffal/teps-problems


text:
function jldopen(f::Function, args...; kws...)
    jld = jldopen(args...; kws...)
    try
        return f(jld)
    finally
        close(jld)
    end
end

"""
    @save filename var1 [var2 ...]

Write one or more variables `var1,...` from the current scope to a JLD2 file
`filename`.

For interactive use you can save all variables in the current module's global
scope using `@save filename`. More permanent code should prefer the explicit
form to avoid saving unwanted variables.

# Example

To save the string `hello` and array `xs` to the JLD2 file example.jld2:

    hello = "world"
    xs = [1,2,3]
    @save "example.jld2" hello foo
"""
macro save(filename, vars...)
    if isempty(vars)
        # Save all variables in the current module
        quote
            let
                m = $(__module__)
                f = jldopen($(esc(filename)), "w")
                wsession = JLDWriteSession()
                try
                    for vname in names(m; all=true)
                        s = string(vname)
                        if !occursin(r"^_+[0-9]*$", s) # skip IJulia history vars
                            v = getfield(m, vname)
                            if !isa(v, Module)
                                try
                                    write(f, s, v, wsession)
                                catch e
                                    if isa(e, PointerException)
                                        @warn("skipping $vname because it contains a pointer")
                                    else
                                        rethrow(e)
                                    end
                                end
                            end
                        end
                    end
                finally
                    close(f)
                end
            end
        end
    else
        writeexprs = Vector{Expr}(undef, length(vars))
        for i = 1:length(vars)
            writeexprs[i] = :(write(f, $(string(vars[i])), $(esc(vars[i])), wsession))
        end

        quote
            jldopen($(esc(filename)), "w") do f
                wsession = JLDWriteSession()
                $(Expr(:block, writeexprs...))
            end
        end
    end
end

"""
    @load filename var1 [var2 ...]

Load one or more variables `var1,...` from JLD2 file `filename` into the
current scope and return a vector of the loaded variable names.

For interactive use, the form `@load "somefile.jld2"` will load all variables
from `"somefile.jld2"` into the current scope. This form only supports literal
file names and should be avoided in more permanent code so that it's clear
where the variables come from.

# Example

To load the variables `hello` and `foo` from the file example.jld2, use

    @load "example.jld2" hello foo
"""
macro load(filename, vars...)
    if isempty(vars)
        if isa(filename, Expr)
            throw(ArgumentError("filename argument must be a string literal unless variable names are specified"))
        end
        # Load all variables in the top level of the file
        readexprs = Expr[]
        vars = Symbol[]
        f = jldopen(filename)
        try
            for n in keys(f)
                if !isgroup(f, lookup_offset(f.root_group, n))
                    push!(vars, Symbol(n))
                end
            end
        finally
            close(f)
        end
    end
    return quote
        ($([esc(x) for x in vars]...),) = jldopen($(esc(filename))) do f
            ($([:(read(f, $(string(x)))) for x in vars]...),)
        end
        $(Symbol[v for v in vars]) # convert to Array
    end
end

# Save all the key-value pairs in the dict as top-level variables of the JLD
function save(f::File{format"JLD2"}, dict::AbstractDict; kwargs...)
    jldopen(FileIO.filename(f), "w"; kwargs...) do file
        wsession = JLDWriteSession()
        for (k,v) in dict
            if !isa(k, AbstractString)
                throw(ArgumentError("keys must be strings (the names of variables), got $k"))
            end
            write(file, String(k), v, wsession)
        end
    end
end

# Or the names and values may be specified as alternating pairs
function save(f::File{format"JLD2"}, name::AbstractString, value, pairs...; kwargs...)
    if isodd(length(pairs)) || !isa(pairs[1:2:end], Tuple{Vararg{AbstractString}})
        throw(ArgumentError("arguments must be in name-value pairs"))
    end
    jldopen(FileIO.filename(f), "w"; kwargs...) do file
        wsession = JLDWriteSession()
        write(file, String(name), value, wsession)
        for i = 1:2:length(pairs)
            write(file, String(pairs[i]), pairs[i+1], wsession)
        end
    end
end

save(f::File{format"JLD2"}, value...; kwargs...) = error("must supply a name for each variable")

function loadtodict!(d::Dict, g::Union{JLDFile,Group}, prefix::String="")
    for k in keys(g)
        v = g[k]
        if v isa Group
            loadtodict!(d, v, prefix*k*"/")
        else
            d[prefix*k] = v
        end
    end
    return d
end

# load with just a filename returns a dictionary containing all the variables
function load(f::File{format"JLD2"}; kwargs...)
    jldopen(FileIO.filename(f), "r"; kwargs...) do file
        loadtodict!(Dict{String,Any}(), file)
    end
end

# When called with explicitly requested variable names, return each one
function load(f::File{format"JLD2"}, varname::AbstractString; kwargs...)
    jldopen(FileIO.filename(f), "r"; kwargs...) do file
        read(file, varname)
    end
end

load(f::File{format"JLD2"}, varnames::AbstractString...; kwargs...) =
    load(f, varnames; kwargs...)

function load(f::File{format"JLD2"}, varnames::Tuple{Vararg{AbstractString}}; kwargs...)
    jldopen(FileIO.filename(f), "r"; kwargs...) do file
        map((var)->read(file, var), varnames)
    end
end


text:
using Eyeball
using Eyeball.FoldingTrees
using REPL
using REPL.TerminalMenus
using InteractiveUtils
using Test

# From REPL/test/TerminalMenus/dynamic_menu.jl
matchback(str) = match(r"\e\[(\d+)A", str)

# from https://github.com/JuliaCollections/FoldingTrees.jl/blob/e4d3c37b18a6ee2c3ad9a5049f6fcc4b79486886/test/runtests.jl#L47-L60
function linesplitter(str)
    strs = split(str, '\n')
    s1 = strs[1]
    m = matchback(s1)
    if m === nothing
        nback, startidx = 0, 1
    else
        nback = parse(Int, m.captures[1])
        startidx = m.offset+length(m.match)
    end
    strs[1] = s1[startidx:end]  # discard the portion that moves the terminal
    @test all(s->startswith(s, "\e[2K"), strs)
    return nback, replace.(map(s->s[5:end], strs), ('\r'=>"",))
end

@testset "Basics" begin
    a = (h=zeros(5), e=:(5sin(pi*t)), f=sin, c=33im, set=Set((:a, 9, 1:5, 8)), b=(c=1,d=9,e=(i=9,f=0)), x=9 => 99:109, d=Dict(1=>2, 3=>4), ds=Dict(:s=>4,:t=>7), dm=Dict(1=>9, "x"=>8))
    root = eye(a, 1, interactive = false)
    @test count_open_leaves(root) == 11
    io = IOBuffer()
    foreach(unfold!, nodes(root))
    FoldingTrees.print_tree(io, root)
    s = String(take!(io))
    @test s[1:100] == "  : \e[32mNamedTuple{(:h, :e, :f, :c, :set, :b, :x, :d, :ds, :dm), Tuple{Vector{Float64}, Expr, typeo"
    root = eye(a, 2, interactive = false)
    @test count_open_leaves(root) == 26
    root = eye(a, interactive = false)
    @test count_open_leaves(root) == 40
end

@testset "TreeMenu" begin
    a = (h=zeros(5), e=:(5sin(pi*t)), f=sin, c=33im, set=Set((:a, 9, 1:5, 8)), b=(c=1,d=9,e=(i=9,f=0)), x=9 => 99:109, d=Dict(1=>2, 3=>4), ds=Dict(:s=>4,:t=>7), dm=Dict(1=>9, "x"=>8))
    root = eye(a, 1, interactive = false)
    menu = TreeMenu(root, pagesize = 20, dynamic = true, maxsize = 30)
    # menu = TreeMenu(root, pagesize = 20, dynamic = true, maxsize = 30, keypress = Eyeball.keypress)
    @test TerminalMenus.numoptions(menu) == 11
    io = IOBuffer()
    state = TerminalMenus.printmenu(io, menu, 2; init=true)
    str = String(take!(io))
    nback, lines = linesplitter(str)
    @test nback == 0
    @test lines[2] == " > +  \e[36mh\e[39m: \e[32mVector{Float64}\e[39m \e[35m(5,)\e[39m \e[33m40\e[39m [0.0, 0.0, 0.0, 0.0, 0.0]"
    @test lines[8] == "   +  \e[36mx\e[39m: \e[32mPair{Int64, UnitRange{Int64}}\e[39m  9=>99:109"
    @test length(lines) == 11
    menu.cursoridx = 11
    # TerminalMenus.keypress(menu, UInt32('f'))
    # @test TerminalMenus.numoptions(menu) == 19
    # node = FoldingTrees.setcurrent!(menu, menu.cursoridx)
    # o = node.data.obj
    # @test o isa Dict
end


text:
#=
Tools for working with Markov Chains

@author : Spencer Lyon, Zac Cranko

@date: 07/10/2014

References
----------

https://lectures.quantecon.org/jl/finite_markov.html

=#
import LightGraphs: DiGraph, period, attracting_components,
                    strongly_connected_components, is_strongly_connected

@inline check_stochastic_matrix(P) = maximum(abs, sum(P, dims = 2) .- 1) < 5e-15 ? true : false

"""
Finite-state discrete-time Markov chain.

Methods are available that provide useful information such as the stationary
distributions, and communication and recurrent classes, and allow simulation of
state transitions.

##### Fields

- `p::AbstractMatrix` : The transition matrix. Must be square, all elements must be nonnegative, and all rows must sum to unity.
- `state_values::AbstractVector` : Vector containing the values associated with the states.
"""
mutable struct MarkovChain{T, TM<:AbstractMatrix{T}, TV<:AbstractVector}
    p::TM # valid stochastic matrix
    state_values::TV

    function MarkovChain{T,TM,TV}(p::AbstractMatrix, state_values) where {T,TM,TV}
        n, m = size(p)

        n != m &&
            throw(DimensionMismatch("stochastic matrix must be square"))

        minimum(p) <0 &&
            throw(ArgumentError("stochastic matrix must have nonnegative elements"))

        !check_stochastic_matrix(p) &&
            throw(ArgumentError("stochastic matrix rows must sum to 1"))

        length(state_values) != n &&
            throw(DimensionMismatch("state_values should have $n elements"))

        return new{T,TM,TV}(p, state_values)
    end
end

# Provide constructor that infers T from eltype of matrix
MarkovChain(p::AbstractMatrix, state_values=1:size(p, 1)) =
    MarkovChain{eltype(p), typeof(p), typeof(state_values)}(p, state_values)

Base.eltype(mc::MarkovChain{T,TM,TV}) where {T,TM,TV} = eltype(TV)

"Number of states in the Markov chain `mc`"
n_states(mc::MarkovChain) = size(mc.p, 1)

function Base.show(io::IO, mc::MarkovChain{T,TM}) where {T,TM}
    println(io, "Discrete Markov Chain")
    println(io, "stochastic matrix of type $TM:")
    print(io, mc.p)
end

@doc doc"""
This routine computes the stationary distribution of an irreducible Markov
transition matrix (stochastic matrix) or transition rate matrix (generator
matrix) ``A``.

More generally, given a Metzler matrix (square matrix whose off-diagonal
entries are all nonnegative) ``A``, this routine solves for a nonzero solution
``x`` to ``x (A - D) = 0``, where ``D`` is the diagonal matrix for which the rows of
``A - D`` sum to zero (i.e., ``D_{ii} = \sum_j A_{ij}`` for all ``i``). One (and only
one, up to normalization) nonzero solution exists corresponding to each
reccurent class of ``A``, and in particular, if ``A`` is irreducible, there is a
unique solution; when there are more than one solution, the routine returns the
solution that contains in its support the first index ``i`` such that no path
connects ``i`` to any index larger than ``i``. The solution is normalized so that
its 1-norm equals one. This routine implements the Grassmann-Taksar-Heyman
(GTH) algorithm (Grassmann, Taksar, and Heyman 1985), a numerically stable
variant of Gaussian elimination, where only the off-diagonal entries of ``A`` are
used as the input data. For a nice exposition of the algorithm, see Stewart
(2009), Chapter 10.

##### Arguments

- `A::Matrix{T}` : Stochastic matrix or generator matrix. Must be of shape n x
  n.

##### Returns

- `x::Vector{T}` : Stationary distribution of ``A``.

##### References

- W. K. Grassmann, M. I. Taksar and D. P. Heyman, "Regenerative Analysis and
  Steady State Distributions for Markov Chains, " Operations Research (1985),
  1107-1116.
- W. J. Stewart, Probability, Markov Chains, Queues, and Simulation, Princeton
  University Press, 2009.

"""
gth_solve(A::Matrix{T}) where {T<:Real} = gth_solve!(copy(A))

# convert already makes a copy, hence gth_solve!
gth_solve(A::Matrix{T}) where {T<:Integer} =
    gth_solve!(convert(Matrix{Rational{T}}, A))

"""
Same as `gth_solve`, but overwrite the input `A`, instead of creating a copy.

"""
function gth_solve!(A::Matrix{T}) where T<:Real
    n = size(A, 1)
    x = zeros(T, n)

    @inbounds for k in 1:n-1
        scale = sum(A[k, k+1:n])
        if scale <= zero(T)
            # There is one (and only one) recurrent class contained in
            # {1, ..., k};
            # compute the solution associated with that recurrent class.
            n = k
            break
        end
        A[k+1:n, k] /= scale

        for j in k+1:n, i in k+1:n
            A[i, j] += A[i, k] * A[k, j]
        end
    end

    # backsubstitution
    x[n] = 1
    @inbounds for k in n-1:-1:1, i in k+1:n
        x[k] += x[i] * A[i, k]
    end

    # normalisation
    x /= sum(x)

    return x
end

"""
Find the recurrent classes of the Markov chain `mc`.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.

##### Returns

- `::Vector{Vector{Int}}` : Vector of vectors that describe the recurrent
  classes of `mc`.

"""
recurrent_classes(mc::MarkovChain) = attracting_components(DiGraph(mc.p))

"""
Find the communication classes of the Markov chain `mc`.

#### Arguments

- `mc::MarkovChain` : MarkovChain instance.

### Returns

- `::Vector{Vector{Int}}` : Vector of vectors that describe the communication
  classes of `mc`.

"""
communication_classes(mc::MarkovChain) = strongly_connected_components(DiGraph(mc.p))

"""
Indicate whether the Markov chain `mc` is irreducible.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.

##### Returns

- `::Bool`

"""
is_irreducible(mc::MarkovChain) =  is_strongly_connected(DiGraph(mc.p))

"""
Indicate whether the Markov chain `mc` is aperiodic.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.

##### Returns

- `::Bool`

"""
is_aperiodic(mc::MarkovChain) = period(mc) == 1

"""
Return the period of the Markov chain `mc`.

##### Arguments

- `mc::MarkovChain` : MarkovChain instance.

##### Returns

- `::Int` : Period of `mc`.

"""
function period(mc::MarkovChain)
    g = DiGraph(mc.p)
    recurrent = attracting_components(g)

    d = 1
    for r in recurrent
        pd = period(g[r])
        d *= div(pd, gcd(pd, d))
    end

    return d
end

for (S, ex_T, ex_gth) in ((Real, :(T), :(gth_solve!)),
                          (Integer, :(Rational{T}), :(gth_solve)))
    @eval function stationary_distributions(mc::MarkovChain{T}) where T<:$S
        n = n_states(mc)
        rec_classes = recurrent_classes(mc)
        T1 = $ex_T
        stationary_dists = Vector{Vector{T1}}(undef, length(rec_classes))

        for (i, rec_class) in enumerate(rec_classes)
            dist = zeros(T1, n)
            # * mc.p[rec_class, rec_class] is a copy, hence gth_solve!
            #   while gth_solve for Int matrix
            # * todense is to convert a sparse matrix to a dense matrix
            #   with eltype T1
            A = todense(T1, mc.p[rec_class, rec_class])
            dist[rec_class] = $ex_gth(A)
            stationary_dists[i] = dist
        end

        return stationary_dists
    end
end

@doc doc"""
Compute stationary distributions of the Markov chain `mc`, one for each
recurrent class.

##### Arguments

- `mc::MarkovChain{T}` : MarkovChain instance.

##### Returns

- `stationary_dists::Vector{Vector{T1}}` : Vector of vectors that represent
  stationary distributions, where the element type `T1` is `Rational` if `T` is
  `Int` (and equal to `T` otherwise).

"""
stationary_distributions

# From base/sparse/sparsematrix.jl
"""Custom version of `full`, which allows convertion to type `T`"""
function todense(T::Type, S::SparseMatrixCSC)
    A = zeros(T, S.m, S.n)
    for Sj in 1:S.n
        for Sk in nzrange(S, Sj)
            Si = S.rowval[Sk]
            Sv = S.nzval[Sk]
            A[Si, Sj] = Sv
        end
    end
    return A
end

"""If `A` is already dense, return `A` as is"""
todense(::Type, A::Array) = A

mutable struct MCIndSimulator{T<:MarkovChain,S}
    mc::T
    len::Int
    init::Int
    drvs::S
end

function MCIndSimulator(mc::MarkovChain, len::Int, init::Int)
    # NOTE: ensure dense array and transpose before slicing the array. Then
    #       when passing to DiscreteRV use `sub` to avoid allocating again
    p = Matrix(mc.p)'
    drvs = [DiscreteRV(view(p, :, i)) for i in 1:size(mc.p, 1)]
    MCIndSimulator(mc, len, init, drvs)
end

# Base.start(mcis::MCIndSimulator) = (mcis.init, 0)
# function Base.next(mcis::MCIndSimulator, state::Tuple{Int,Int})
#     ix, t = state
#     (ix, (rand(mcis.drvs[ix]), t+1))
# end
# Base.done(mcis::MCIndSimulator, s::Tuple{Int,Int}) = s[2] >= mcis.len

function Base.iterate(mcis::MCIndSimulator, state::Tuple{Int,Int}=(mcis.init, 0))
    ix, t = state
    if t >= mcis.len
        return nothing
    end
    (ix, (rand(mcis.drvs[ix]), t+1))
end

Base.length(mcis::MCIndSimulator) = mcis.len
Base.eltype(mcis::MCIndSimulator) = Int
Base.IteratorSize(mcis::MCIndSimulator) = Base.HasLength()

mutable struct MCSimulator{T<:MCIndSimulator}
    mcis::T
end

function MCSimulator(mc::MarkovChain, len::Int, init::Int)
    MCSimulator(MCIndSimulator(mc, len, init))
end

# only need to implement next and eltype differently...
function Base.iterate(mcs::MCSimulator, state::Tuple{Int,Int}=(mcs.mcis.init, 0))
    output = iterate(mcs.mcis, state)
    if output === nothing
        return output
    end
    ix, new_state = output
    (mcs.mcis.mc.state_values[ix], new_state)
end
Base.eltype(mcs::MCSimulator) = eltype(mcs.mcis.mc)

# ...the rest of the interface can derive from mcis
Base.length(mcs::MCSimulator) = length(mcs.mcis)
Base.IteratorSize(mcs::MCSimulator) = Base.IteratorSize(mcs.mcis)

"""
Simulate one sample path of the Markov chain `mc`.
The resulting vector has the state values of `mc` as elements.

### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of simulation
- `;init::Int=rand(1:n_states(mc))` : Initial state

### Returns

- `X::Vector` : Vector containing the sample path, with length
  ts_length
"""
function simulate(mc::MarkovChain, ts_length::Int;
                  init::Int=rand(1:n_states(mc)))
    X = Vector{eltype(mc)}(undef, ts_length)
    simulate!(X, mc; init=init)
end

"""
Fill `X` with sample paths of the Markov chain `mc` as columns.
The resulting matrix has the state values of `mc` as elements.

### Arguments

- `X::Matrix` : Preallocated matrix to be filled with sample paths
of the Markov chain `mc`. The element types in `X` should be the
same as the type of the state values of `mc`
- `mc::MarkovChain` : MarkovChain instance.
- `;init=rand(1:n_states(mc))` : Can be one of the following
    - blank: random initial condition for each chain
    - scalar: same initial condition for each chain
    - vector: cycle through the elements, applying each as an
      initial condition until all columns have an initial condition
      (allows for more columns than initial conditions)
"""
function simulate!(X::Union{AbstractVector,AbstractMatrix},
                   mc::MarkovChain; init=rand(1:n_states(mc), size(X, 2)))
    mcs = MCSimulator(mc, size(X, 1), init[1])

    for (i, init) in enumerate(take(cycle(init), size(X, 2)))
        mcs.mcis.init = init
        copyto!(view(X, :, i), mcs)
    end
    X
end

# ------------------------ #
# simulate_indices methods #
# ------------------------ #

"""
Simulate one sample path of the Markov chain `mc`.
The resulting vector has the indices of the state values of `mc` as elements.

### Arguments

- `mc::MarkovChain` : MarkovChain instance.
- `ts_length::Int` : Length of simulation
- `;init::Int=rand(1:n_states(mc))` : Initial state

### Returns

- `X::Vector{Int}` : Vector containing the sample path, with length
  ts_length
"""
function simulate_indices(mc::MarkovChain, ts_length::Int;
                          init::Int=rand(1:n_states(mc)))
    X = Vector{Int}(undef, ts_length)
    simulate_indices!(X, mc; init=init)
end

"""
Fill `X` with sample paths of the Markov chain `mc` as columns.
The resulting matrix has the indices of the state values of `mc` as elements.

### Arguments

- `X::Matrix{Int}` : Preallocated matrix to be filled with indices
of the sample paths of the Markov chain `mc`.
- `mc::MarkovChain` : MarkovChain instance.
- `;init=rand(1:n_states(mc))` : Can be one of the following
    - blank: random initial condition for each chain
    - scalar: same initial condition for each chain
    - vector: cycle through the elements, applying each as an
      initial condition until all columns have an initial condition
      (allows for more columns than initial conditions)
"""
function simulate_indices!(X::Union{AbstractVector{T},AbstractMatrix{T}},
               mc::MarkovChain; init=rand(1:n_states(mc), size(X, 2))) where T<:Integer
    mcis = MCIndSimulator(mc, size(X, 1), init[1])

    for (i, init) in enumerate(take(cycle(init), size(X, 2)))
        mcis.init = init
        copyto!(view(X, :, i), mcis)
    end
    X
end

@doc doc"""
Simulate one sample path of the non-stationary vector of Markov chains `mcs`.
The resulting vector has the state values of `mcs` as elements.

### Arguments

- `mcs::Vector{MarkovChain}` : Vector of MarkovChains.
- `;init::Int` : Index of initial state

### Returns

- `X::Vector` : Vector containing the sample path, with length
  length(mcs)+1 (includes initial state)
"""
function simulate_ns(mcs::Vector, init::Int)
    ind_0 = init
    ind = zeros(Int, length(mcs))
    X = zeros(Real, length(mcs))

    for t in eachindex(mcs)
        X[t] = simulate(mcs[t], 2, init = ind_0)[2]
        ind[t] = findmin( abs.(mcs[t].state_values .- X[t]) )[2]
        if t < length(mcs)
            ind_0 = ind[t] #findmin( abs.(mcs[t].state_values .- X[t]) )[2]
        end
    end
    ind = [init; ind]
    X = [mcs[1].state_values[init]; X]
end


text:
using Test
using CapacityExpansion
using Clp

@testset "CEP" begin
    @testset "TX_1" begin #Compare to Merrick Testcase
        # load data
        ts_input_data = load_timeseries_data_provided("TX_1"; T=24, years=[2008])
        cep_data = load_cep_data_provided("TX_1")
        # run optimization
        model = run_opt(ts_input_data,cep_data,Clp.Optimizer;optimizer_config=Dict{Symbol,Any}(:LogLevel => 0))
        # compare to exact result
        exact_res=[70540.26439790576;0.0;0.0;8498.278397905757;80132.88454450261;0.0]
        @test exact_res ≈ model.variables["CAP"].data[:,1,1] atol=1
    end
    scenarios=Dict{String,OptResult}()
    #Test the workflow for single node scenarios for Germany and California
    @testset "workflow $state" for (state, years) in [["GER_1", [2016]],["CA_1", [2016]]] begin
            # laod data
            ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
            cep_data = load_cep_data_provided(state)
            ## CLUSTERING ##
            ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
            ts_full_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=30)
            ## OPTIMIZATION ##
            optimizer=Clp.Optimizer
            optimizer_config=Dict{Symbol,Any}(:LogLevel => 0)
            scenarios["$state-$years-co2"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="co2",limit_emission=Dict{String,Number}("CO2/electricity"=>1000),optimizer_config=optimizer_config)
            scenarios["$state-$years-slack"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="slack",lost_load_cost=Dict{String,Number}("electricity"=>1e6), lost_emission_cost=Dict{String,Number}("CO2"=>700),optimizer_config=optimizer_config)
            scenarios["$state-$years-ex"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="ex",infrastructure=Dict{String,Array}("existing"=>["all"]),optimizer_config=optimizer_config)
            scenarios["$state-$years-simple"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="simple storage",storage_type="simple",conversion=true,optimizer_config=optimizer_config)
            scenarios["$state-$years-seasonal"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="seasonal storage",storage_type="seasonal",conversion=true,optimizer_config=optimizer_config)
            design_result=run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="des&op",optimizer_config=optimizer_config)
            scenarios["$state-$years-des&op"] = run_opt(ts_full_data.clust_data,cep_data,design_result.config,get_cep_design_variables(design_result),optimizer;lost_load_cost=Dict{String,Number}("electricity"=>1e6), lost_emission_cost=Dict{String,Number}("CO2"=>700))
        end
    end
    #Test transmission for a multi-node scenario
    @testset "workflow $state" for (state, years) in [["GER_18", [2016]],["CA_14", [2016]]] begin
           # laod data
           ts_input_data = load_timeseries_data_provided(state; T=24, years=years) #CEP
           cep_data = load_cep_data_provided(state)
           ## CLUSTERING ##
           ts_clust_data = run_clust(ts_input_data;method="hierarchical",representation="centroid",n_init=1,n_clust=3)
           ## OPTIMIZATION ##
           optimizer=Clp.Optimizer
           scenarios["$state-$years-trans"] = run_opt(ts_clust_data.clust_data,cep_data,optimizer;descriptor="trans",transmission=true,optimizer_config=Dict{Symbol,Any}(:LogLevel => 0))
       end
   end
   #Test exact values for each of the previously calculated scenarios by comparison with exact scenarios
   @load normpath(joinpath(dirname(@__FILE__),"cep_exact_data.jld2")) exact_scenarios
   @testset "exact $k" for (k,v) in scenarios begin
            #Test for each variable within this scenario
            @testset "variable $kv" for (kv,vv) in v.variables begin
                    @test sum((round.(round.(exact_scenarios[k].variables[kv].data; digits=4); sigdigits=4).== round.(round.(vv.data;digits=4); sigdigits=4)).== false) == 0
                end
            end
        end
    end
end


text:
"""
    Tullio.@einsum  A[i,j] += B[i] * C[j]

Since this package is almost superset of `Einsum.jl`, you can probable drop that and
write `using Tullio: @einsum` to use the new macro under the old name.

Differences:
* Updating `A` with weird things like `*=` uses an awful hack which may be less efficient,
  but does make tests pass!
* Options `threads=false, avx=false, grad=false` are selected for you.
"""
macro einsum(ex::Expr)
    if ex.head in [:(:=), :(=), :(+=)]
        _tullio(ex, :(avx=false), :(threads=false), :(grad=false); mod=__module__)

    elseif ex.head in [:(-=), :(*=), :(/=)]
        @gensym tmp

        if @capture_(ex.args[1], Z_[ijk__]) # array *= ...
            act = Expr(:(:=), :($tmp[$(ijk...)]), ex.args[2:end]...)
            res = _tullio(act, :(avx=false), :(threads=false), :(grad=false); mod=__module__).args[1]
            Expr(Symbol(string(".", ex.head)), Z, res) |> esc

        elseif ex.args[1] isa Symbol # scalar case
            Z = ex.args[1]
            act = Expr(:(:=), :($tmp), ex.args[2:end]...)
            res = _tullio(act, :(avx=false), :(threads=false), :(grad=false); mod=__module__).args[1]
            Expr(ex.head, Z, res) |> esc

        end
    end
end


text:
function _mat_str(str; repl=false)
    lhs, rhs = _split_expression(str)
    lhs_parsed = _try_parse(lhs)

    if rhs == ""
        lhs = interpolate(lhs)
        nargout = Int(!repl)
        return quote $(_convert_py_output(eng.eval(lhs, nargout=nargout))) end

    elseif lhs_parsed isa Expr && lhs_parsed.head == :($)
        rhs = _convert_py_output(eng.eval(rhs, nargout=1))
        return esc(:($(lhs_parsed.args[1]) = $rhs))

    elseif lhs_parsed isa Expr && (lhs_parsed.head == :vect || lhs_parsed.head == :hcat)
        # rhs = _convert_py_output.(py"eng.eval($rhs, nargout=$(length(lhs_parsed.args)))")
        rhs = eng.eval(rhs, nargout=length(lhs_parsed.args))

        block = map(zip(lhs_parsed.args, rhs)) do (arg, expr)
            if arg isa Expr && arg.head ==:($)
                esc(:($(arg.args[1]) = $(_convert_py_output(expr))))

            else
                str = string(arg) * "=" * _string_string(expr)
                quote $(eng.eval(str, nargout=0)) end
            end
        end

        return Expr(:block, block...)

    else
        str = lhs * "=" * rhs
        return quote $(eng.eval(str, nargout=0)) end
    end
end

"""
    mat"...MATLAB code..."
    out = mat"...MATLAB code..."

Runs MATLAB code in open MATLAB engine session

## Examples
```julia
julia> using SharedMATLABEngine

julia> connect_matlab();
REPL mode MATLAB initialized. Press > to enter and backspace to exit.

>> a = magic(3)

a =

     8     1     6
     3     5     7
     4     9     2

>> \$a21 = a(2,1)
3.0

julia> a21 .+ mat"a"
3×3 Array{Float64,2}:
 11.0   4.0   9.0
  6.0   8.0  10.0
  7.0  12.0   5.0
```
"""
macro mat_str(expr)
    _mat_str(expr)
end

macro matrepl_str(expr)
    _mat_str(expr; repl=true)
end


text:
# __precompile__(false)

module Sunny

using Requires
using LinearAlgebra
using StaticArrays
using OffsetArrays
using SpecialFunctions
using FFTW
using Tullio
using ProgressMeter
using Printf
using Random: rand!, randn!

# Specific to Symmetry/
using FilePaths: Path
using CrystalInfoFramework
import Spglib

# TODO: Remove in Julia 1.7
using Parameters: @unpack

const Vec3 = SVector{3, Float64}
const Mat3 = SMatrix{3, 3, Float64, 9}

# Boltzmannn factor k_B in units of meV/K
const meV_per_K = 0.086173332621451774

# Bohr magneton in units of meV / T
const BOHR_MAGNETON = 0.057883818060738013331

# Vacuum permability in units of T^2 Å^3 / meV
const VACUUM_PERM = 201.33545383470705041

include("Symmetry/Symmetry.jl")
export Crystal, subcrystal, nbasis, cell_volume, cell_type
export lattice_vectors, lattice_params
export Bond, displacement, distance, coordination_number
export print_bond, print_bond_table
export reference_bonds, basis_for_symmetry_allowed_couplings
export all_symmetry_related_bonds, all_symmetry_related_bonds_for_atom
export all_symmetry_related_couplings, all_symmetry_related_couplings_for_atom

include("Util.jl")

include("Lattice.jl")

include("Interactions.jl")
export heisenberg, exchange, dm_interaction
export easy_axis, easy_plane, single_ion_anisotropy
export external_field, dipole_dipole
export SiteInfo

include("PairInteractions.jl")

include("Ewald.jl")

include("FourierAccel.jl")

include("Hamiltonian.jl")

include("Systems.jl")
export ChargeSystem, SpinSystem, rand!, randflips!, energy, field, field!

include("Metropolis.jl")
export MetropolisSampler, IsingSampler, set_temp!, get_temp, get_system
export sample!, thermalize!, anneal!
export running_energy, running_mag, reset_running_energy!, reset_running_mag!

include("Integrators.jl")
export HeunP, LangevinHeunP, SphericalMidpoint, evolve!
export LangevinSampler

include("StructureFactors.jl")
export StructureFactor, update!, apply_dipole_factor, zero!
export dynamic_structure_factor, static_structure_factor

include("WangLandau/BinnedArray.jl")
export BinnedArray, filter_visited, reset!

include("WangLandau/WangLandau.jl")
export WangLandau, spherical_cap_update, init_bounded!, run!

function __init__()
    @require GLMakie="e9467ef8-e4e7-5192-8a1a-b1aee30e663a" begin
        include("Plotting.jl")
        export plot_lattice, plot_spins, plot_bonds, plot_all_bonds
        export anim_integration, live_integration, live_langevin_integration
    end

    @require MPI="da04e1cc-30fd-572f-bb4f-1f8673147195" begin
        include("ReplicaExchangeMC.jl")
        export init_MPI, xyz_to_file, Replica, run_REMC!, run_FBO!
    end
end

end


text:
module IfsFractal

export FractalCanvas,Orbit
export randcomplex,mandelbrot,draw!,render

using Images
using Colors

@inline function randcomplex(rng::MersenneTwister,maxnorm::Float64=2.0)
	c=Complex(maxnorm*(rand(rng)*2-1),maxnorm*(rand(rng)*2-1))
	while distsqrd(c,zero(Complex)) > maxnorm*maxnorm
		c=Complex(maxnorm*(rand(rng)*2-1),maxnorm*(rand(rng)*2-1))
	end
	return c
end
type Pixel
	x::Int
	y::Int
end
type FractalCanvas
	width::Int
	height::Int
	center::Complex
	zoom::Float64
	data::Matrix{Float64}
end
FractalCanvas(width::Int,height::Int,center::Complex,zoom::Float64) =
	FractalCanvas(width,height,center,zoom,zeros(Float64,height,width))

@inline project(canvas::FractalCanvas,pnt::Complex) =
	Pixel(round(Int,((real(pnt)-real(canvas.center))*canvas.zoom+0.5)*canvas.width),
		round(Int,((imag(pnt)-imag(canvas.center))*canvas.zoom+0.5)*canvas.height))
@inline inbounds(canvas::FractalCanvas,pixel::Pixel) =
	pixel.x>=1 && pixel.x <=canvas.width && pixel.y >= 1 && pixel.y <= canvas.height

function writepixel!(canvas::FractalCanvas,pixel::Pixel)
	if inbounds(canvas,pixel)
		canvas.data[pixel.y,pixel.x]+=1.0
	end
end

type SamplerData
	cs::Vector{Complex}
	contribs::Vector{Float64}
end
SamplerData()=SamplerData(Vector{Complex}(),Vector{Float64}())
function add!(sd::SamplerData,c::Complex,contrib::Float64)
	push!(sd.cs,c)
	push!(sd.contribs,contrib)
end

type Orbit
	c::Complex
	zs::Vector{Complex}
	escaped::Bool
end
Orbit(c::Complex)=Orbit(c,Vector{Complex}(),false)
const MAX_ORBIT_LENGTH=50000
const DEFAULT_NUM_MANDELBROT_ITERS=200
@inline function mandelbrot(c::Complex,num_iters::Int=DEFAULT_NUM_MANDELBROT_ITERS)
	z=zero(Complex)
	orbit=Orbit(c)
	# const threshold=2^16
	const threshold=4
	for i=1:num_iters
		z=z*z+orbit.c
		push!(orbit.zs,z)
		if real(z*z') > threshold
			orbit.escaped=true
			return orbit
		end
		if length(orbit.zs) >= MAX_ORBIT_LENGTH
			break
		end
	end
	return orbit
end

@inline function contrib(canvas::FractalCanvas,orbit::Orbit)
	ctrib=0
	for i=1:length(orbit.zs)
		pixel=project(canvas,orbit.zs[i])
		if inbounds(canvas,pixel)
			ctrib+=1
		end
	end
	return ctrib/length(orbit.zs)
end

function generate_sample_points(rng::MersenneTwister,canvas::FractalCanvas)
	const num_samplers=30
	const num_sampling_iters=50000
	samplerdata=SamplerData()
	for i=1:num_samplers
		c=find_initial_point(rng,canvas,Complex(0.0,0.0),2.0)
		if isinf(c)
			warn("Unable to find initial poitn for sampler $i")
			continue
		end
		mandelbrot(c,num_sampling_iters)
		cont=contrib(canvas,orbit)
		add!(samplerdata,c,cont)
	end
	return samplerdata
end

distsqrd(c1::Complex,c2::Complex)=(real(c2)-real(c1))*(real(c2)-real(c1))+
	(imag(c2)-imag(c1))*(imag(c2)-imag(c1))

function find_initial_point(rng::MersenneTwister,canvas::FractalCanvas,pnt::Complex,radius::Float64,
		depth::Int=0)
	const num_search_attempts=200
	const num_search_iters=50000
	const max_depth=500
	if depth > max_depth
		return Complex(Inf,Inf)
	end
	closest=Inf
	nextseed=pnt
	for i=1:num_search_attempts
		c=randcomplex(rng,radius)+pnt
		orbit=mandelbrot(c,num_search_iters)
		if !orbit.escaped
			continue
		end
		contribution=contrib(canvas,orbit)
		if contribution > 0.0
			return c
		end
		for j=1:length(orbit.zs)
			ds=distsqrd(orbit.zs[j],canvas.center)
			if ds < closest
				closest=ds
				nextseed=c
			end
		end
	end
	return find_initial_point(rng,canvas,nextseed,radius/2.0,depth+1)
end

function burnin!(canvas::FractalCanvas,samplerdata::SamplerData)
	for i=1:length(samplerdata.cs)
		const warmup_length=10000
		for j=1:warmup_length
		end
	end
end

function loop!(canvas::FractalCanvas)
end

function main()
	rng=MersenneTwister(1)
	canvas=FractalCanvas(512,512,zero(Complex),0.25)
	samplerdata=generate_sample_points(rng,canvas)
	burnin!(canvas,samplerdata)
	loop!(canvas)
	f=open("canvas.metro.dat","w")
	serialize(f,canvas)
	close(f)
	render(canvas,"out.metro.dat")
end

function draw!(canvas::FractalCanvas,orbit::Orbit)
	for z in orbit.zs
		px=project(canvas,z)
		# @show orbit.zs[i],px,inbounds(canvas,px)
		writepixel!(canvas,px)
	end
end

const DEFAULT_GAIN_COEF=0.2
@inline compute_bias_coef(gain_coef::Float64)=log(1.0-gain_coef)/log(0.5)
@inline bias(val::Float64,bias_coef::Float64) = val>0.0 ? val^bias_coef : 0.0
@inline gain(val::Float64,bias_coef::Float64) =
	0.5*(val<0.5 ? bias(2.0*val,bias_coef) : 2.0-bias(2.0-2.0*val,bias_coef))
@inline clamp01(val::Float64) = clamp(val,0.0,1.0)

function compute_pixels!(imagedata::Matrix{RGB{Float64}},canvas::FractalCanvas,
		gain_coef::Float64,bgcolor::RGB{Float64},fgcolor::RGB{Float64})
	maxval=maximum(canvas.data)
	minval=minimum(canvas.data)
	@show maxval,minval,mean(canvas.data)
	bias_coef=compute_bias_coef(gain_coef)
	for j=1:canvas.width
		for i=1:canvas.height
			# pixel=clamp01(gain(canvas.data[i,j]/maxval,bias_coef)*2.0)
			pixel=clamp01(2.0*(canvas.data[i,j]-minval)/(maxval-minval))^0.5
			imagedata[i,j]=RGB{Float64}(
				(1.0-pixel)*bgcolor.r+pixel*fgcolor.r,
				(1.0-pixel)*bgcolor.g+pixel*fgcolor.g,
				(1.0-pixel)*bgcolor.b+pixel*fgcolor.b)
		end
	end
end
function render(canvas::FractalCanvas,filename::AbstractString,
		gain_coef::Float64=DEFAULT_GAIN_COEF,
		bgcolor::RGB{Float64}=RGB{Float64}(1.0,1.0,1.0),
		fgcolor::RGB{Float64}=RGB{Float64}(0.0,0.0,0.0))
	# @show canvas.data
	imagedata=fill(bgcolor,(canvas.height,canvas.width))
	compute_pixels!(imagedata,canvas,gain_coef,bgcolor,fgcolor)
	img=Image(imagedata)
	save(filename,img)
end

end

max_stars_repo_name: jblondin/Fractal

max_issues_repo_name: jblondin/Fractal

max_forks_repo_name: jblondin/Fractal