text:
__precompile__()

module BotCoreLCMTypes

using LCMCore
using StaticArrays: SVector

export atlas_command_t,
       force_torque_t,
       gps_data_t,
       gps_satellite_info_list_t,
       gps_satellite_info_t,
       image_metadata_t,
       image_sync_t,
       image_t,
       images_t,
       ins_t,
       joint_angles_t,
       joint_state_t,
       kvh_raw_imu_batch_t,
       kvh_raw_imu_t,
       planar_lidar_t,
       pointcloud2_t,
       pointcloud_t,
       pointfield_t,
       pose_t,
       position_3d_t,
       quaternion_t,
       raw_t,
       rigid_transform_t,
       robot_state_t,
       robot_urdf_t,
       sensor_status_t,
       six_axis_force_torque_array_t,
       six_axis_force_torque_t,
       system_status_t,
       twist_t,
       utime_t,
       vector_3d_t,
       viewer_command_t,
       viewer_draw_t,
       viewer_geometry_data_t,
       viewer_link_data_t,
       viewer_load_robot_t

include("atlas_command_t.jl")
include("force_torque_t.jl")
include("gps_data_t.jl")
include("gps_satellite_info_t.jl")
include("gps_satellite_info_list_t.jl")
include("image_metadata_t.jl")
include("image_sync_t.jl")
include("image_t.jl")
include("images_t.jl")
include("ins_t.jl")
include("joint_angles_t.jl")
include("joint_state_t.jl")
include("kvh_raw_imu_t.jl")
include("kvh_raw_imu_batch_t.jl")
include("planar_lidar_t.jl")
include("pointfield_t.jl")
include("pointcloud2_t.jl")
include("pointcloud_t.jl")
include("pose_t.jl")
include("vector_3d_t.jl")
include("quaternion_t.jl")
include("position_3d_t.jl")
include("twist_t.jl")
include("raw_t.jl")
include("rigid_transform_t.jl")
include("robot_state_t.jl")
include("robot_urdf_t.jl")
include("sensor_status_t.jl")
include("six_axis_force_torque_t.jl")
include("six_axis_force_torque_array_t.jl")
include("system_status_t.jl")
include("utime_t.jl")
include("viewer_command_t.jl")
include("viewer_draw_t.jl")
include("viewer_geometry_data_t.jl")
include("viewer_link_data_t.jl")
include("viewer_load_robot_t.jl")
end


text:
using C3D, Test

cd(@__DIR__)
datadir = joinpath("..", "data")

include("identical.jl")
include("bigdata.jl")
include("singledata.jl")
include("invalid.jl")
include("blanklabels.jl")
include("badformats.jl")

ext: jl


text:
using Pkg.Pkg.TOML
using Printf

# directory where Julia packages are stored
pkg_dir = "/Users/mengliu/.julia/packages/"

# this variable needs to be changed based on your own installation directory
path = "/Users/mengliu/.julia/registries/General"
packages_dict = TOML.parsefile(joinpath(path,"Registry.toml"))["packages"]
# this variable needs to be changed based on your own installation directory
const STDLIB_DIR = "/Applications/Julia-1.3.app/Contents/Resources/julia/share/julia/stdlib/v1.3/"
const STDLIBS = readdir(STDLIB_DIR)

# find all Julia standard packages
std_pkg_names = []
for (i, stdlib) in enumerate(STDLIBS)
    if isfile(joinpath(STDLIB_DIR, stdlib, "Project.toml"))
        proj = TOML.parsefile(joinpath(STDLIB_DIR, stdlib, "Project.toml"))
        push!(std_pkg_names,proj["name"])
    end
end

# next, we use the first package as an example
std_pkg_names[1]

joinpath(STDLIB_DIR,std_pkg_names[1])

# grab the first source file
curr_path = joinpath(joinpath(STDLIB_DIR,std_pkg_names[1]),"src")
src_file = readdir(curr_path)[1]
curr_file = joinpath(curr_path,src_file)

# read source file as string
src = read(curr_file,String)

# parse source file
src_expr = Meta.parse(src)

# this function will go thourgh all the expressions in the source file
found_functions = Set()
function find_functions(expr::Union{GlobalRef,
    LineNumberNode,String,Symbol,Bool,QuoteNode},
    found_functions) end # do nothing
function find_functions(expr::Expr,found_functions)
    for arg in expr.args
        # expr.head will store which function is called in this expression
        if typeof(arg) == Expr
            push!(found_functions, expr.head)
        end
        find_functions(arg, found_functions)
    end
end

find_functions(src_expr,found_functions)

## TODO
# run find_functions on all source files of all packages
# build a mapping between all possible function names and node indices
# build co-calling graph

max_stars_repo_name: empet/more-graphs

max_issues_repo_name: empet/more-graphs

max_forks_repo_name: empet/more-graphs


text:
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie

using ForwardDiff
Base.adjoint(f::Function) = x -> ForwardDiff.derivative(f, float(x))

function tangent_line(f=nothing, a=0, b=pi)

    if f == nothing
        f = x -> sin(x)
    end

descr = """
The tangent line has a slope approximated by the slope of secant lines.
This demo allows the points c and c+h to be adjusted to see the two lines"""

## From http://juliaplots.org/MakieReferenceImages/gallery//edit_polygon/index.html:
function add_move!(scene, points, pplot)
    idx = Ref(0); dragstart = Ref(false); startpos = Base.RefValue(Point2f0(0))
    on(events(scene).mousedrag) do drag
        if ispressed(scene, Mouse.left)
            if drag == Mouse.down
                plot, _idx = mouse_selection(scene)
                if plot == pplot
                    idx[] = _idx; dragstart[] = true
                    startpos[] = to_world(scene, Point2f0(scene.events.mouseposition[]))
                end
            elseif drag == Mouse.pressed && dragstart[] && checkbounds(Bool, points[], idx[])
                pos = to_world(scene, Point2f0(scene.events.mouseposition[]))

                # we work with components not vector
                x,y = pos

                x = clamp(x, a, b)
                y = f(x)

                points[][idx[]] = [x,y]
                points[] = points[]
            end
        else
            dragstart[] = false
        end
        return
    end
end

c, h = pi/4, .5
points = Node(Point2f0[(c, f(c)), (c+h, f(c+h))])

# where we lay our scene:
scene, layout = layoutscene()
layout.halign = :left
layout.valign = :top

p = layout[1,1:2] = LScene(scene)

rowsize!(layout, 1, Auto(1))
colsize!(layout, 1, Auto(1))

layout[2,1:2] = LText(scene, descr)

secline = lift(points) do pts
    c, ch = pts
    x0, y0 = c
    x1, y1 = ch
    m = (y1 - y0)/(x1 - x0)
    sl = x -> y0 + m * (x - x0)
    Point2f0[(a, sl(a)), (b, sl(b))]
end

tangentline = lift(points) do pts
    c, ch = pts
    x0, y0 = c
    m = f'(x0)
    tl = x -> y0 + m * (x - x0)
    Point2f0[(a, tl(a)), (b, tl(b))]
end

lines!(p.scene, a..b, f, strokecolor=:red, strokewidth=15)
lines!(p.scene, secline, color = :blue, strokewidth = 10, raw=true)
lines!(p.scene, tangentline, color = :red, strokewidth = 10, raw=true)
scatter!(p.scene, points, color = :white, strokewidth = 10, markersize = 0.05, strokecolor = :black, raw = true)
xlims!(p.scene, (a, b))
#ylims!(p.scene, (0, 1.5))

add_move!(p.scene, points, p.scene[end])

scene

    end

tangent_line()


text:
using Documenter
using TensorCast

makedocs(
    sitename = "TensorCast",
    modules = [TensorCast],
    pages = [
        "Home" => "index.md",
        "Basics" => "basics.md",
        "Reduction" => "reduce.md",
        "Options" => "options.md",
        "Docstrings" => "docstrings.md",
    ],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
)

deploydocs(
    repo = "github.com/mcabbott/TensorCast.jl"
)

# julia --color=yes make.jl


text:
"""
    SeisLinearEvents(; <keyword arguments>)

Generate five dimensional data `d` consisting of linear events.

# Arguments

**Keyword arguments**

* `ot=0.0`: first sample for the time axis in secs.
* `dt=0.004`: sampling interval in secs.
* `nt=500`: number of time samples.
* `ox1=0.0`: first sample for the first spatial dimension in meters.
* `dx1=10.0`: sample interval for the first spatial dimension in meters.
* `nx1=100`: number of samples for the first spatial dimension.
* `ox2=0.0`: first sample for the second spatial dimension in meters.
* `dx2=10.0`: sample interval for the second spatial dimension in meters.
* `nx2=1`: number of samples for the second spatial dimension.
* `ox3=0.0`: second sample for the third spatial dimension in meters.
* `dx3=10.0`: sample interval for the third spatial dimension in meters.
* `nx3=1`: number of samples for the third spatial dimension.
* `ox4=0.0`: third sample for the fourth spatial dimension in meters.
* `dx4=10.0`: sample interval for the fourth spatial dimension in meters.
* `nx4=1`:number of samples for the fourth spatial dimension.
* `tau=[1.0, 1.6]`: intercept traveltimes for each event.
* `p1=[0.0000,-0.0001]`
* `p2=[0.0003, 0.0002]`
* `p3=[-0.0001,-0.0001]`
* `p4=[0.0001,-0.0000]`
* `amp=[1.0,-1.0]`: amplitudes for each linear event.
* `wavelet="ricker"`: wavelet used to model the linear events.
* `f0=[20.0]`: central frequency of wavelet for each linear event.

# Example
```julia
julia> d,extent = SeisLinearEvents(); SeisPlot(d);
```

Credits: Aaron Stanton, 2015
"""
function SeisLinearEvents(; ot=0.0, dt=0.004, nt=500, ox1=0.0, dx1=10.0,
                          nx1=100, ox2=0.0, dx2=10.0, nx2=1, ox3=0.0, dx3=10.0,
                          nx3=1, ox4=0.0, dx4=10.0, nx4=1, tau=[1.0,1.5],
                          p1=[0.0001,-0.0003],p2=[0.,0.],p3=[0.,0],p4=[0.,0.],
                          amp=[1.0,-1.0], f0=20.0)

    w = Ricker(dt=dt,f0=f0)
    nf = nextpow2(nt)
    nw = length(w)
    t_delay = (nw-1)*dt/2
    w = vcat(w, zeros(nf-nw))
    W = fft(w)
    x1 =ox1+collect(0:1:nx1-1)*dx1
    x2 =ox2+collect(0:1:nx2-1)*dx2
    x3 =ox3+collect(0:1:nx3-1)*dx3
    x4 =ox4+collect(0:1:nx4-1)*dx4
    n_events = length(p1)
    d = zeros(Float64, nt, nx1, nx2, nx3, nx4)
    D = zeros(Complex{Float64}, nf, nx1, nx2, nx3, nx4)
    nfh = round(Int, floor(nf/2)) + 1
    wrs = collect(0:1:nfh-1)*2*pi/(nf*dt)     # Frequency in rad/sec
    for ie = 1:n_events
     for ix4 = 1:nx4
      for ix3 = 1:nx3
       for ix2 = 1:nx2
        for ix1 = 1:nx1
          for iw = 2:nfh-1
              phase = wrs[iw]*(tau[ie] + p1[ie]*x1[ix1] + p2[ie]*x2[ix2]
                               + p3[ie]*x3[ix3] + p4[ie]*x3[ix4] - t_delay)
              D[iw:iw, ix1,ix2,ix3,ix4] += W[iw]*amp[ie]*exp(-im*phase)
              D[nf-iw+2,ix1,ix2,ix3,ix4] = conj(D[iw,ix1,ix2,ix3,ix4])
          end
        end
       end
      end
     end
    end
    d = ifft(D,1)
    d = real(d[1:nt,:,:,:,:]);
    # Add extent header
    extent = Seismic.Extent(convert(Int32, nt), convert(Int32, nx1),
                            convert(Int32, nx2), convert(Int32, nx3),
                            convert(Int32, nx4), convert(Float32, ot),
                            convert(Float32, ox1), convert(Float32, ox2),
                            convert(Float32, ox3), convert(Float32, ox4),
                            convert(Float32, dt), convert(Float32, dx1),
                            convert(Float32, dx2), convert(Float32, dx3),
                            convert(Float32,dx4), "Time", "ix1", "ix2", "ix3",
                            "ix4", "s", "index", "index", "index", "index", "")
    if extent.n5 == 1 && extent.n4 == 1 && extent.n3 == 1 && extent.n2 == 1
        d = reshape(d, round(Int, extent.n1))
    elseif extent.n5 == 1 && extent.n4 == 1 && extent.n3 == 1
        d = reshape(d, round(Int, extent.n1), round(Int, extent.n2))
    elseif extent.n5 == 1 && extent.n4 == 1
        d = reshape(d, round(Int, extent.n1), round(Int, extent.n2),
                    round(Int, extent.n3))
    elseif extent.n5 == 1
        d = reshape(d, round(Int, extent.n1), round(Int, extent.n2),
                    round(Int, extent.n3), round(Int, extent.n4))
    else
        d = reshape(d, round(Int, extent.n1), round(Int, extent.n2),
                    round(Int, extent.n3), round(Int, extent.n4),
                    round(Int, extent.n5))
    end
    return d, extent

end


text:
###
  Name     : c8_24_01.jl
  Book     : Hands-on Data Science with Anaconda )
  Publisher: Packt Publishing Ltd.
  Author   : Yuxing Yan and James Yan
  Date     : 3/25/2018
  email    : yany@canisius.edu
             paulyxy@hotmail.com
###

using Base.Test
using ClusterAnalysis

lbls = vec([1 1 1 2 2 2])
expected = vec([2 2 2 3 3 3])

@test 1.0 == ClusterAnalysis.adjusted_rand_index(lbls, expected)

lbls = vec([1 1 1 2 2 2])
expected = vec([1 2 1 1 2 1])


text:
# DigilentWaveForms.jl: wrapper for Digilent WaveForms SDK
# Copyright (C) 2019 Ian Butterworth

# wrap.jl: automatically wrap C API using Clang.jl
# Builded by wrap.jl within Spinnaker.jl by Samuel Powell

@info("NOTE MANUAL CHANGES in deps/wrap.jl")

using Clang
const LIBDWF_INCLUDE = joinpath("/Library", "Frameworks", "dwf.framework", "Headers") |> normpath
const LIBDWF_HEADERS = [joinpath(LIBDWF_INCLUDE, header) for header in ["dwf.h"]]

wc = init(; headers = LIBDWF_HEADERS,
            output_file = joinpath(@__DIR__, "libdwf_api.jl"),
            common_file = joinpath(@__DIR__, "libdwf_common.jl"),
            clang_includes = vcat(LIBDWF_INCLUDE, CLANG_INCLUDE),
            clang_args = ["-I", joinpath(LIBDWF_INCLUDE, "..")],
            header_wrapped = (root, current)->root == current,
            header_library = x->"libdwf",
            clang_diagnostics = true,
            )

run(wc)

open(joinpath(@__DIR__, "libdwf_fixes.jl"), "w") do io
    write(io, "# manual fixes\n\n")
    write(io, "const TRUE = true #fix for first line in libdwf_common.jl\n\n")
end

# Manually wrapping the consts, as dwf.h uses consts instead of enums
include("wrapconsts.jl")
wrapconsts(LIBDWF_HEADERS[1], joinpath(@__DIR__, "libdwf_consts.jl"))


text:
using MD5: md5

part1(input) = md5search(input, 0x00f0ffff)
part2(input) = md5search(input, 0x00ffffff)

function md5search(input, mask)
    salt = readchomp(input)
    i = 1
    while first(reinterpret(UInt32, md5(string(salt, i)))) & mask != 0
        i += 1
    end
    return i
end


text:
module BinaryAnalysis

using MachO

import Base.show
export get_dltree

abstract node
abstract leaf <: node
type branch <: node
    path::String
    children::Array{node,1}
end

type good_leaf <: leaf
    path::String
end

type bad_leaf <: leaf
    path::String
    error::String
end

nodemap = Dict{String,node}()

function get_dltree(path, rpaths=String[])
    # Return cached results if we've got them
    if haskey(nodemap, path)
        return nodemap[path]
    end

    if !isfile(path)
        return bad_leaf(path, "File not found")
    end

    lcs = nothing
    h = readmeta(path)
    if isa(h, MachO.FatMachOHandle)
        lcs = [collect(LoadCmds(h[1])), collect(LoadCmds(h[2]))]
    else
        lcs = collect(LoadCmds(h))
    end

    # Get rpaths from the binary, if they exist
    new_rpath_cmds = filter(x -> isa(x, MachO.LoadCmd{MachO.rpath_command}), lcs)
    new_rpaths = [x.cmd.path for x in new_rpath_cmds]

    # Canonicalize the new_rpaths and add them to rpaths
    new_rpaths = [abspath(replace(x, "@executable_path", dirname(path))) for x in new_rpaths]
    prepend!(new_rpaths, rpaths)

    # Get only LOAD_DYLIB lc's
    dylib_cmds = filter(x -> isa(x, MachO.LoadCmd{MachO.dylib_command}) && x.cmd_id == MachO.LC_LOAD_DYLIB, lcs)
    dylibs = [x.cmd.name for x in dylib_cmds]

    # For any dylib paths that might have @rpath in them, search through rpaths for them
    dylibs = map(dylibs) do search_path
        if contains(search_path, "@rpath")
            for rpath in new_rpaths
                poss_dylib_path = abspath(replace(search_path, "@rpath", rpath))

                # See if we can descend down the tree this way
                poss_dylib = get_dltree(poss_dylib_path, new_rpaths)
                if isa(poss_dylib, good_leaf) || isa(poss_dylib, branch)
                    return poss_dylib
                end
            end

            # If none of our rpaths worked out, complain
            return bad_leaf(search_path, "File not found")
        else
            return get_dltree(search_path)
        end
    end

    # dylibs is an array of nodes
    ret = good_leaf(path)
    if length(dylibs) > 0
        ret = branch(path, dylibs)
    end

    nodemap[path] = ret
    return ret
end

# Convenience function for showing a
function show(io::IO, n::node, tablevel=0)
    if tablevel >= 2
        print(io, "| "^(tablevel - 1))
    end
    if tablevel >= 1
        print(io, "+-")
    end

    if isa(n,branch)
        println(io, n.path, ":")
        for c in n.children
            show(io, c, tablevel + 1 )
        end
    else
        if isa(n,good_leaf)
            println(io, n.path)
        else
            println(io, n.path, ": ", n.error)
        end
    end
end

end #module


text:
include("../src/CholeskyExperiment.jl")
using .CholeskyExperiment: ExperimentManager,
                           DecompositionDB
const db = DecompositionDB

function main()
    manager = ExperimentManager()
    db.delete_similar_decomposition(manager.decompositions; verbose=true)
end
main()

max_stars_repo_name: charlyalizadeh/OPFdecgen

max_issues_repo_name: charlyalizadeh/OPFdecgen

max_forks_repo_name: charlyalizadeh/OPFdecgen


text:
using Pkg
Pkg.activate("")
using DataFrames, Dates, JLD2

list_problems = "list_problems_1000000.dat"
# list_problems = "list_problems_test.dat"
problems = readlines(list_problems)

# problems = ["ARGLINA", "ARGLINB", "ARGLINC"]
npb = length(problems)

for solver in ["ARC", "TRU"]
    name_table = String[] #
    time_table = Float64[] # Vector{Float64}(undef, npb)
    neval_obj_table = Int[] # Vector{Float64}(undef, npb)
    neval_grad_table = Int[] # Vector{Float64}(undef, npb)
    neval_hprod_table = Int[] # Vector{Float64}(undef, npb)
    objective_table = Float64[] # Vector{Float64}(undef, npb)
    dual_feas_table = Float64[] # Vector{Float64}(undef, npb)
    iter_table = Int[]

    for pb in problems
        push!(name_table, pb)
        result = readlines("$(pb)/$(solver)/" * pb * ".dat")
        neval_obj = NaN
        neval_grad = NaN
        neval_hprod = NaN
        objective = NaN
        dual_feas = NaN
        iter = NaN
        for line in result
            if occursin("Total time", line)
                push!(
                    time_table,
                    parse(Float64, replace(split(line, "=")[2], "seconds" => "")),
                )
            end
            if occursin("function evaluations", line)
                neval_obj = parse(Int, split(line, "=")[2])
            end
            if occursin("gradient evaluations", line)
                neval_grad = parse(Int, split(line, "=")[2])
            end
            if occursin("Hessian evaluations", line)
                neval_hprod = parse(Int, split(line, "=")[2])
            end
            if occursin("objective value", line)
                objective = parse(Float64, split(line, "=")[2])
            end
            if occursin("gradient norm", line)
                dual_feas = parse(Float64, split(line, "=")[2])
            end
            if occursin("major  iterations", line)
                iter = parse(Int, split(line, "=")[2])
            end
        end
        push!(neval_obj_table, neval_obj)
        push!(neval_grad_table, neval_grad)
        push!(neval_hprod_table, neval_hprod)
        push!(objective_table, objective)
        push!(dual_feas_table, dual_feas)
        push!(iter_table, iter)
    end
    df = DataFrame(
        [
            name_table,
            time_table,
            objective_table,
            iter_table,
            dual_feas_table,
            neval_obj_table,
            neval_grad_table,
            neval_hprod_table,
        ],
        [
            :name,
            :elapsed_time,
            :objective,
            :iter,
            :dual_feas,
            :neval_obj,
            :neval_grad,
            :neval_hprod,
        ],
    )

    open("$(today())_$(solver)_cutest_$(npb).dat", "w") do io
        println(io, df)
    end

    using JLD2
    @save "$(today())_$(solver)_cutest_$(npb).jld2" df
end


text:
export htmldoc, htmlroot, htmlhead, htmlbody, attributes, attribute_names, attribute_values,
    pushfirst_sibling!, push_sibling!, pushfirst_child!, push_child!,
    HTMLDocument, HTMLElement, HTMLNode, HTMLText, NullNode, attrs, getattr, setattr!, tag, text

import Base: string

const EmptyHtmlDoc = Gumbo.parsehtml("")
htmldoc() = deepcopy(EmptyHtmlDoc)

push_sibling!(elem::Gumbo.HTMLElement,val) = push!(elem, val)
pushfirst_sibling!(elem::Gumbo.HTMLElement,val) = pushfirst!(elem, val)
push_child!(elem::Gumbo.HTMLElement,val) = push!(elem.children, val)
pushfirst_child!(elem::Gumbo.HTMLElement,val) = pushfirst!(elem.children, val)

const HTMLTextTag = "(_)"
Gumbo.tag(x::Gumbo.HTMLText) = HTMLTextTag

htmldoctype(x::Gumbo.HTMLDocument) = x.doctype
htmlroot(x::Gumbo.HTMLDocument) = x.root
htmlhead(x::Gumbo.HTMLDocument) = Gumbo.children(x.root)[1]
htmlbody(x::Gumbo.HTMLDocument) = Gumbo.children(x.root)[2]

EzXML.parentelement(x::Gumbo.HTMLElement) = x.parent

# subtypes(HTMLNode) == [ HTMLElement, HTMLText, NullNode ]

typeofnode(x::HTMLElement{T}) where {T} = T
typeofnode(x::HTMLText) = String
typeofnode(x::NullNode) = Nothing

attributes(x::HTMLElement{T}) where {T} = x.attributes
# attribute_pairs(x::HTMLElement{T}) where {T} =
attribute_names(x::HTMLElement{T}) where {T}  = keys(attributes(x))
attribute_values(x::HTMLElement{T}) where {T} = values(attributes(x))

function string(::Type{HTMLElement{:HTML}}, x, lang="en")
  doctypestr = "<!DOCTYPE html>"
  htmlstr = string("")
  headstr = string(HTMLElement{:head}, x[1])
  bodystr = string(HTMLElement{:body}, x[2])
  return join((doctypestr, htmlstr, headstr, "", bodystr, "", "\n"), "\n")
end

function string(::Type{HTMLElement{:head}}, x)

end

function string(::Type{HTMLElement{:body}}, x)
    atrs = attrs(x)
    if isempty(atrs)
        result = "\n"
    else
        result = string("\n")
    end
    for child in children(x)
        result = string(result, "  ", string(typeof(child), child))
    end
    result = string(result, "\n")
    return result
end

function string(::Type{HTMLElement{T}}, x) where {T}
    attrs = attributes(x)
    isempty(attrs) && return stringcontent(HTMLElement{T})
    childs = children(x)
    while !isempty(childs)
       return string(typeof(childs), childs)
    end
    return ""
end

function stringcontent(::Type{HTMLElement{T}}, x) where {T}
  content =
  return string("<",T,"></",T,">")
end

function string_of_attributes(attributes::Dict) where {T}
    k, v = collect(keys(attributes)), collect(values(attributes))
    vstrs = map(s->string("\"",s,"\""), v)
    z = zip(k,vstrs)
    return join(join.(collect(z), "="), " ")
end

#=
prettyhtml(x::AbstractString) = prettyhtml(Gumbo.parsehtml(x))

function prettyhtml(x::Gumbo.HTMLDocument)
    io = IOBuffer()
    Gumbo.prettyprint(io, x)
    s = String(take!(io))
    return s
end
=#

#=
htmlstr = """
<!DOCTYPE html>

    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1" />
    <title>Title</title>
    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/bulma/0.7.5/css/bulma.min.css"></link>

      content

"""

htmlstr = """
<!DOCTYPE html>

    <meta charset="utf-8" />
    <meta name="viewport" content="width=device-width, initial-scale=1" />
    <title>Title</title>
    <link rel="stylesheet" href="https://cdnjs.cloudflare.com/ajax/libs/bulma/0.7.5/css/bulma.min.css" />

          these words

"""

htmldoc = Gumbo.parsehtml(htmlstr)
htmlroot = htmldoc.root
htmlhead = htmlroot[1]
htmlbody = htmlroot[2]

julia> text(r[2][2][1][1])
"these words"

julia> source = "
first div

second div
"
"
first div

second div
"

julia> substitute(match) = match[1] * "\n" * match[2]
substitute (generic function with 1 method)

julia> replace(source, r">." => substitute)
"
\nfirst div
\n
\nsecond div
"

julia> print(ans)

first div

second div

=#


text:
struct ShaderSource
    code::Vector{UInt8}
    language::ShaderLanguage
    stage::Vk.ShaderStageFlag
    entry_points::Vector{Symbol}
end

function Base.show(io::IO, source::ShaderSource)
    print(io, "ShaderSource(", source.language, ", ", source.stage, ", ", length(source.code), " bytes)")
end

function pad_shader_code!(code::Vector{UInt8})
    size = cld(length(code), 4)
    rem = size * 4 - length(code)
    if rem ≠ 0
        resize!(code, size * 4)
        code[end - rem + 1:end] .= 0
    end
    @assert length(code) % 4 == 0
    code
end

struct Shader
    source::ShaderSource
    shader_module::Vk.ShaderModule
    entry_point::Symbol
    push_constant_ranges::Vector{Vk.PushConstantRange}
    specialization_constants::Vector{Vk.SpecializationInfo}
end

device(shader::Shader) = shader.shader_module.device

Shader(source, shader_module, entry_point) = Shader(source, shader_module, entry_point, [], [])

struct ShaderCache
    device::Vk.Device
    compiled::Dictionary{String,ShaderSource}
    shaders::Dictionary{ShaderSource,Shader}
end

Base.hash(source::ShaderSource, h::UInt64) = objectid(source.code) + h

ShaderCache(device) = ShaderCache(device, Dictionary(), Dictionary())

function find_source!(cache::ShaderCache, spec::ShaderSpecification)
    file = string(spec.source_file)
    if haskey(cache.compiled, file)
        source = cache.compiled[file]
        if spec.entry_point in source.entry_points
            return source
        end
    else
        source = ShaderSource(spec)
        if source.language ≠ SPIR_V
            source = compile(source)
        end
        insert!(cache.compiled, file, source)
        source
    end
end

function ShaderSource(spec::ShaderSpecification)
    ShaderSource(read(spec.source_file), spec.language, spec.stage, [spec.entry_point])
end

function find_shader!(cache::ShaderCache, source::ShaderSource, entry_point::Symbol)
    if source.language ≠ SPIR_V
        source = compile(source)
    end
    find_shader!(cache, source, ShaderSpecification("", false, entry_point, SPIR_V))
end

find_shader!(cache::ShaderCache, spec::ShaderSpecification) = find_shader!(cache, find_source!(cache, spec), spec)

function find_shader!(cache::ShaderCache, source::ShaderSource, spec::ShaderSpecification)
    if haskey(cache.shaders, source)
        cache.shaders[source]
    else
        shader_module = Vk.ShaderModule(cache.device, source)
        shader = Shader(source, shader_module, spec.entry_point)
        insert!(cache.shaders, source, shader)
        shader
    end
end

"""
Retrieve a shader from the provided cache.

Note that the cache will be modified if it does not contain the requested shader.
"""
Shader(cache::ShaderCache, spec::ShaderSpecification) = find_shader!(cache, spec)


text:
function findRow(A::Array,b::Array)
 # finds the first row in A such that A[k,:] = b

 N = round(Int128,maximum(A[:]))
 Av = arrayToUniqueInt(A,N)
 bv = arrayToUniqueInt(b,N)
 return indexin(bv,Av)

end

function arrayToUniqueInt(A::Array{Int64,2},N::Int128)

	k = zeros(Int128,size(A,1))
	for i=1:size(A,1)
		for j=1:size(A,2)
		 	k[i] += A[i,j]*N^(j-1)
		end
	end
	return k
end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
function measure_components(properties::Tuple{MeasurementProperties, Vararg{MeasurementProperties}}, labels::AbstractArray; N::Int = 0)
    N = N == 0 ? maximum(labels) : N
    # Set primary key as the component label.
    t = table(Base.OneTo(N); names = [:l], pkey=:l)
    # Some measurements require that we have computed so-called `Bit Quads`.
    usecodes = any(map(x->typeof(x) <: Union{Measurement, BasicTopology}, properties))
    t = usecodes ? append_bitcodes(t, labels, N) : t
    for i = 1:length(properties)
        t = measure_feature(properties[i], t, labels, N)
    end
    t
end

function append_bitcodes(t::IndexedTable, labels::AbstractArray,N::Int)
    # Stores counts of Bit Quad patterns for each component.
    𝓠₀ = zeros(Int, N)
    𝓠₁ = zeros(Int, N)
    𝓠₂ = zeros(Int, N)
    𝓠₃ = zeros(Int, N)
    𝓠₄ = zeros(Int, N)
    𝓠ₓ = zeros(Int, N)
    rows, cols = axes(labels)
    labels_padded = padarray(labels, Fill(0, (2,2)))
    rows, cols = axes(labels_padded)
    for r = first(rows):last(rows)-1
        for c = first(rows):last(rows)-1
            for n = 1:N
                p₁ = labels_padded[r, c] == n ? 1 : 0
                p₂ = labels_padded[r + 1, c] == n ? 1 : 0
                p₃ = labels_padded[r, c + 1] == n ? 1 : 0
                p₄ = labels_padded[r + 1, c + 1] == n ? 1 : 0
                𝓠₀[n] = Q₀(p₁, p₂, p₃, p₄) ? 𝓠₀[n] + 1 : 𝓠₀[n]
                𝓠₁[n] = Q₁(p₁, p₂, p₃, p₄) ? 𝓠₁[n] + 1 : 𝓠₁[n]
                𝓠₂[n] = Q₂(p₁, p₂, p₃, p₄) ? 𝓠₂[n] + 1 : 𝓠₂[n]
                𝓠₃[n] = Q₃(p₁, p₂, p₃, p₄) ? 𝓠₃[n] + 1 : 𝓠₃[n]
                𝓠₄[n] = Q₄(p₁, p₂, p₃, p₄) ? 𝓠₄[n] + 1 : 𝓠₄[n]
                𝓠ₓ[n] = Qₓ(p₁, p₂, p₃, p₄) ? 𝓠ₓ[n] + 1 : 𝓠ₓ[n]
            end
        end
    end
    pushcol(t, :Q₀ => 𝓠₀, :Q₁ => 𝓠₁, :Q₂ => 𝓠₂, :Q₃ => 𝓠₃, :Q₄ => 𝓠₄, :Qₓ => 𝓠ₓ)
end

function Q₀(p₁, p₂, p₃, p₄)
    c₁ = (p₁ == 0 && p₂ == 0 && p₃ == 0 &&  p₄ == 0) ? true : false
end

function Q₁(p₁, p₂, p₃, p₄)
    c₁ = (p₁ == 1 && p₂ == 0 && p₃ == 0 &&  p₄ == 0) ? true : false
    c₂ = (p₁ == 0 && p₂ == 0 && p₃ == 1 &&  p₄ == 0) ? true : false
    c₃ = (p₁ == 0 && p₂ == 0 && p₃ == 0 &&  p₄ == 1) ? true : false
    c₄ = (p₁ == 0 && p₂ == 1 && p₃ == 0 &&  p₄ == 0) ? true : false
    c₁ || c₂ || c₃ || c₄
end

function Q₂(p₁, p₂, p₃, p₄)
    c₁ = (p₁ == 1 && p₂ == 0 && p₃ == 1 &&  p₄ == 0) ? true : false
    c₂ = (p₁ == 0 && p₂ == 0 && p₃ == 1 &&  p₄ == 1) ? true : false
    c₃ = (p₁ == 0 && p₂ == 1 && p₃ == 0 &&  p₄ == 1) ? true : false
    c₄ = (p₁ == 1 && p₂ == 1 && p₃ == 0 &&  p₄ == 0) ? true : false
    c₁ || c₂ || c₃ || c₄
end

function Q₃(p₁, p₂, p₃, p₄)
    c₁ = (p₁ == 1 && p₂ == 0 && p₃ == 1 &&  p₄ == 1) ? true : false
    c₂ = (p₁ == 0 && p₂ == 1 && p₃ == 1 &&  p₄ == 1) ? true : false
    c₃ = (p₁ == 1 && p₂ == 1 && p₃ == 0 &&  p₄ == 1) ? true : false
    c₄ = (p₁ == 1 && p₂ == 1 && p₃ == 1 &&  p₄ == 0) ? true : false
    c₁ || c₂ || c₃ || c₄
end

function Q₄(p₁, p₂, p₃, p₄)
    c₁ = (p₁ == 1 && p₂ == 1 && p₃ == 1 &&  p₄ == 1) ? true : false
end

function Qₓ(p₁, p₂, p₃, p₄)
    c₁ = (p₁ == 1 && p₂ == 0 && p₃ == 0 &&  p₄ == 1) ? true : false
    c₂ = (p₁ == 0 && p₂ == 1 && p₃ == 1 &&  p₄ == 0) ? true : false
    c₁ || c₂
end


text:
N = [[1 2 3 4], [5 6], [7 8 9 ]]
plateList = [[2,3,4,8],[5,6,9]]

nPlate = length(plateList)

nNodes = 0
for i in 1:length(N)
    nNodes += length(N[i])
end

nodeList = collect(1:nNodes)

plateLevelList = []

for i in 1: nPlate
    h = []
    for k in 1:length(plateList[i])
        for j in 1: length(N)
            if plateList[i][k] ∈ N[j]
                push!(h,j)
            end
        end
    end
    push!(plateLevelList, h)
end

rangeList=[]
for i in 1: nPlate
    h = []
    push!(h,minimum(plateLevelList[i]), maximum(plateLevelList[i]))
    push!(rangeList, h)
end

#model = Model(with_optimizer(Mosek.Optimizer, QUIET=false, INTPNT_CO_TOL_DFEAS=1e-7))
model = Model(GLPK.Optimizer)

nX = nNodes^2 - nNodes

nA = 0
for j in 1:nPlate
    for i in rangeList[j][1] : rangeList[j][2]
        nA += length(setdiff(N[i], plateList[j]))
    end
end

@variable(model, A[1:nA] >= 0)
@variable(model, B[1:nA], Bin)
@variable(model, X[1:nNodes,1:nNodes], Int)
#@variable(model, Y[1:14,1:2], Bin)

@objective(model, Min, sum(A[i] for i in 1:nA))

w = 1
for j in 1:nPlate
    for i in rangeList[j][1] : rangeList[j][2]
        for k in setdiff(N[i], plateList[j])
            @constraint(model, sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) <=
                A[w] + length(plateList[j]))
            @constraint(model, -sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) <=
                A[w] + length(plateList[j]))
            @constraint(model, sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) >=
                length(plateList[j]) - A[w] - 1000*B[w])
            @constraint(model, sum(X[k,plateList[j][h]] for h in 1:length(plateList[j])) <=
                1000 - length(plateList[j]) + A[w] - 1000*B[w])
            w = w + 1
        end
    end
end

for i in 1:(nNodes-1)
    for j in (i+1):nNodes
        @constraint(model, X[i,j] + X[j,i] == 0)
    end
end

for i in 1:nNodes
        @constraint(model, X[i,i] == 0)
end

for i in 1:nNodes
    for j ∈ nodeList[setdiff(1:nNodes, i)]
        for k ∈ nodeList[setdiff(1:nNodes, [i,j])]
            @constraint(model, X[i,j] + X[j,k] - X[i,k]  <= 1)
        end
    end
end

@constraint(model, -1 .<= X .<= 1)

optimize!(model)
@show value.(X);

            adj_list = Vector{Int}[
                [2,4], #node1
                [5], #node2
                [5, 7, 8], #node3
                [], #node4
                [8], #node 5
                [7], #node 6
                [8], #node 7
                []
            ]

            labels = ["1", "2", "3", "4", "5", "6", "7", "8"]

            layoutTree(adj_list,labels,filename="tree_1.svg",
            cycles = false, ordering = :barycentric, coord = :optimal,
            xsep = 50, x_start = 100, ysep = 100, scale = 0.2, labelpad = 1.2,
            background = "white")

            N = [[1 2 3 4], [5 6], [7 8  ]]
            plateList = [[2,3,4,8],[5,6]]
            Pos = getpos(N, plateList)
            hierviz(N, plateList, Pos)


text:
# 6 Bus Model based on http://motor.ece.iit.edu/data/6bus_Data_ES.pdf

# OPTIONS for caseID:
# wind_reserves: wind generator can provide reserves. this is the default
# no_wind_reserves: wind generator cannot provide reserves
# nondispatchable_wind: wind generators are not dispatchable

const caseID = "no_wind_reserves"

# include basic model parameters
include("params.jl")

#*========== BASIC TOPOLOGY ===========

# number of buses, generators (regular and wind), loads
const numBuses = 6
const numGen = 3
const numRegGen = 2
const numWindGen = 1
const numLoads = 3

# map regular generators to buses
const regGenToBus = [1,6]

# map wind generators to buses
const windGenToBus = [2]

# map loads to buses
const loadToBus = [3,4,5]

# lines between buses. assume that the graph is connected
const lines = ((1,2), (1,4), (2,3), (2,4), (3,6), (4,5), (5,6))

# if the wind generator is nondispatchable, then treat it like a negative load
if(caseID == "nondispatchable_wind")
	numLoads = 4
	numWindGen = 0
	numGen = 2
	loadToBus = [2,3,4,5]
end

#*========== GENERATOR BOUNDS ===========

# lower and upper power bounds for regular generators (pmin and pmax)
const regGenMinPow = zeros(Float64,numRegGen)
const regGenMaxPow = zeros(Float64,numRegGen)
# set high enough upper bound for the most expensive (second) generator to ensure relatively complete recourse
regGenMaxPow[1] = 220.0
regGenMaxPow[2] = 1000.0

# lower power bounds for wind generators. upper bounds are random
const windGenMinPow = zeros(Float64,numWindGen)

#*========== COST COEFFICIENTS ===========

# regular generator costs including cost of reserves
# only use linear cost coefficients
const regGenCost = [13.5,17.7]

# scale reserves cost by suitable factor
const regGenPosResCost = reservesCostFactor*regGenCost
const regGenNegResCost = 1.0*regGenPosResCost

# wind generator costs including cost of reserves (assume zero marginal costs)
const windGenCost = zeros(Float64,numWindGen)

# set cost of reserves for each wind generator to be a scaled version of minimum regular generator reserves cost
const minRegGenPosResCost = minimum(regGenPosResCost)
const minRegGenNegResCost = minimum(regGenNegResCost)

const windGenPosResCost = minRegGenPosResCost*windReservesCostFactor*ones(Float64,numWindGen)
const windGenNegResCost = minRegGenNegResCost*windReservesCostFactor*ones(Float64,numWindGen)

# allow wind spill (with participation factor 0.1*min_alpha) in the no_wind_reserves case
if(caseID == "no_wind_reserves")
	windGenNegResCost = zeros(Float64,numWindGen)
end

#*========== RESERVES INFORMATION ===========

# which regular generators are ready to provide reserves?
const regGenRes = trues(numRegGen)

# which regular generators are ALWAYS READY to provide reserves?
# make sure that the generator with large enough max generation limit always provides reserves to ensure relatively complete recourse holds
# we usually assume that all generators provide guaranteed reserves with a minimum participation factor min_alpha (that isn't too small) so that numerical issues don't crop up while running the stochastic approximation method
const regGenGuarRes = trues(numRegGen)

# which wind generators are ready to provide reserves?
const windGenRes = trues(numWindGen)
if(caseID == "no_wind_reserves" || caseID == "nondispatchable_wind")
	windGenRes = falses(numWindGen)
end

# which wind generators are ALWAYS READY to provide reserves?
const windGenGuarRes = trues(numWindGen)
if(caseID == "no_wind_reserves" || caseID == "nondispatchable_wind")
	windGenGuarRes = falses(numWindGen)
end

# minimum participation factors for generators ALWAYS READY to provide reserves
# reset this min participation factor to not restrict the space of solutions too much, if necessary
const min_alpha = 0.001
const numGenProvRes = sum(windGenGuarRes) + sum(regGenGuarRes)
if(numGenProvRes > 10)
	min_alpha = 1.0/(100.0*numGenProvRes)
end

#*========== LINE PARAMETERS ===========

# these line parameters are assumed to be in the same order as "lines"
# line susceptances (beta)
const beta = 1./[0.170,0.258,0.037,0.197,0.018,0.037,0.140]

# line flow limits
const lineFlowFactor = 1.0
const flow_max = [200.0,100.0,100.0,100.0,100.0,100.0,100.0]*lineFlowFactor

#*========== LOAD PARAMETERS ===========

# power demands/loads
const demandFactor = 1.0

# average demands
const totalMeanDemand = 280.0*demandFactor
const meanDemand = zeros(Float64,numLoads)
meanDemand[1] = 0.2*totalMeanDemand
meanDemand[2] = 0.4*totalMeanDemand
meanDemand[3] = 0.4*totalMeanDemand

#*========== WIND PARAMETERS ===========

# Wind prediction data
const windFactor = 1.4
const windPredMean = [100.0*windFactor]

# in the nondispatchable case, treat wind generator like a negative load
if(caseID == "nondispatchable_wind")
	meanDemand = [-windPredMean;meanDemand[1:3]]
end

#*========== UNCERTAIN PARAMETERS ===========

# standard deviation of loads/demands
const demandStdevFactor = 0.1
const demandStdev = demandStdevFactor*meanDemand

# randomly generated correlation matrix for the demands
const demandCorrMat = [[1.0 0.5827 0.8654];[0.5827 1.0 0.1761];[0.8654 0.1761 1.0]]

# standard deviation of wind power output
const windPredStdevFactor = 0.1
const windPredStdev = windPredStdevFactor*windPredMean

if(caseID == "nondispatchable_wind")
	# use smaller standard deviation to ensure that relatively complete recourse holds
	demandStdevFactor = 0.05
	demandStdev = demandStdevFactor*meanDemand
	demandCorrMat = [[1.0 0.0 0.0 0.0];[0.0 1.0 0.5827 0.8654];[0.0 0.5827 1.0 0.1761];[0.0 0.8654 0.1761 1.0]]
end

#*========== SCENARIO GENERATION ===========

# store the Cholesky factorization of the demand covariance matrix for numerical efficiency
const demandCovMat = Diagonal(demandStdev)*demandCorrMat*Diagonal(demandStdev)
const demandCovMat_chol = (cholfact(demandCovMat))[:L]

# tailored implementation of multivariate normal distribution for demand scenarios
function generateDemandScenarios(numSamples::Int64)

	omega_d = zeros(Float64,numLoads,numSamples)

	omega_d_tmp = rand(Normal(0.0,1.0),numLoads,numSamples)
	for samp = 1:numSamples
		omega_d[:,samp] = demandCovMat_chol*omega_d_tmp[:,samp] + meanDemand
	end

	return omega_d
end

# Generate demand and wind scenarios
function generateScenarios(numSamples::Int64)

	omega_d = generateDemandScenarios(numSamples)

	omega_w = Array{Float64}(numWindGen,numSamples)
	for i = 1:numWindGen
		if(windPredStdev[i] > 0)
			omega_w[i,:] = windPredMean[i] + rand(Normal(0,windPredStdev[i]), numSamples)
		else
			omega_w[i,:] = windPredMean[i]*ones(Float64,numSamples)
		end
	end

	return omega_d, omega_w
end

max_stars_repo_name: rohitkannan/DCOPF-reserve-saturation

max_issues_repo_name: rohitkannan/DCOPF-reserve-saturation

max_forks_repo_name: rohitkannan/DCOPF-reserve-saturation


text:
module NCP_EM

include("../Misc.jl")

using .Misc
using Distributions, SpecialFunctions
import Base.Iterators: product

export standard_EM

function standard_EM(X::Array{ℜ,N}, R::Ƶ; a=1.0, b=a/nansum(X), EPOCHS::Ƶ=1, ϵ=1e-16) where {ℜ<:Real, Ƶ<:Integer, N}
    I = size(X)
    α_I = (a/R) ./ I
    α_R = a/R

    S_R, S₊ = zeros(ℜ,R), sum(X)
    S_I = [Array{ℜ}(undef,R,Iₙ) for Iₙ in I]

    log_λ = log(max(S₊ + a - 1, ϵ)) - log(b + 1)
    log_θ_R = log.(rand(Dirichlet(R,α_R + 0.1)))
    log_θ_I = [log.(rand(Dirichlet(Iₙ,α_Iₙ + 0.1),R)) for (Iₙ,α_Iₙ) in zip(I,α_I)]

    s, log_p = Array{ℜ}(undef,R), Array{ℜ}(undef,R)

    for eph=1:EPOCHS
        # initialization
        for n ∈ 1:N
            S_I[n] .= 0
        end
        S_R .= 0

        # expectation step
        for i=product(map(Iₙ -> 1:Iₙ, I)...)
            log_p .= log_θ_R + sum(n -> log_θ_I[n][i[n],:],1:N)
            log_p .-= logsumexp(log_p)
            s .= X[i...] .* exp.(log_p)

            for (n,iₙ) ∈ enumerate(i)
                S_I[n][:,iₙ] .+= s
            end
            S_R .+= s
        end

        # maximization step
        for (n,Iₙ) ∈ enumerate(I)
            log_θ_I[n] .= log.(max.(S_I[n] .+ α_I[n] .- 1, ϵ))'
            log_θ_I[n] .-= logsumexp(log_θ_I[n],dims=1)
        end
        log_θ_R .= log.(max.(S_R .+ α_R .- 1, ϵ))
        log_θ_R .-= logsumexp(S_R)
    end
    return log_λ, (log_θ_R, log_θ_I...)
end

end

max_stars_repo_name: mehmetburakkurutmaz/bam

max_issues_repo_name: mehmetburakkurutmaz/bam

max_forks_repo_name: atcemgil/bam


text:
function equals(obj::BigRealField, arg0::Object)
    return jcall(obj, "equals", jboolean, (Object,), arg0)
end

function get_instance(::Type{BigRealField})
    return jcall(BigRealField, "getInstance", BigRealField, ())
end

function get_one(obj::BigRealField)
    return jcall(obj, "getOne", BigReal, ())
end

function get_runtime_class(obj::BigRealField)
    return jcall(obj, "getRuntimeClass", Class, ())
end

function get_zero(obj::BigRealField)
    return jcall(obj, "getZero", BigReal, ())
end

function hash_code(obj::BigRealField)
    return jcall(obj, "hashCode", jint, ())
end


text:
############################################################
## get methods
## - reorganise result_type here?

Base.map(f::Function, p::TextParse.AbstractToken) =
    map(f, parser(p))

function Base.get(parser::WrappedParser, sequence, till, after, i, state)
    get(parser.parser, sequence, till, after, i, state)
end

"""
    Base.get(parser::Bytes{N,T}, sequence::Vector{UInt8})

!!! note

    Endianness can be achieved by just mapping `bswap`
    ```jldoctest
    julia> map(bswap, Bytes(2,UInt16))([0x16,0x11])
    0x1611

    julia> Bytes(2,UInt16)([0x16,0x11])
    0x1116
    ```
"""
function Base.get(parser::Bytes{N,T}, sequence::Vector{UInt8}, till,
                  after, i, state) where {N,T}
    if isbitstype(T)
        reinterpret(T,sequence[i:after-1])[1]
    else
        T(sequence[i:after-1])
    end
end

function Base.get(parser::AbstractTokenParser,
                  sequence, till,
                  after, i, state)
    state.state
end

function Base.get(parser::Union{NIndexParser{1}, ConstantParser{Char}}, sequence, till, after, i, state)
    @inbounds sequence[i]
end

function Base.get(parser::ConstantParser, sequence, till, after, i, state)
    parser.parser
end

# Base.get(parser::ConstantParser, sequence, till, after, i, state) = parser.parser

function Base.get(parser::FlatMap, sequence, till, after, i, state)
    li = rightof(sequence,i,parser.left,tuple_pos(state))
    get(right_parser(state),
        sequence, till,
        after, li,
        right_state(state))
end

function Base.get(parser::Either,
                  sequence, till,
                  after, i, state)
    get(parser.options[either_state_option(state)],
        sequence, till,
        after, i, either_state_state(state))
end

function Base.get(parser::Either{<:AbstractTrie},
                  sequence, till,
                  after, i, state)
    get(state.state)
end

if VERSION<v"1.5"
    @warn "julia < 1.5: defining ismutable function for creating copies of default Optional values (only on Union{AbstractDict, AbstractVector, AbstractSet})."
    ismutable(x::Union{AbstractDict, AbstractVector, AbstractSet}) = true
    ismutable(x) = false
end

"""
    _copy(x)

`copy(x)` iif `ismutable(x)`; used when [`defaultvalue`](@ref) of [`Optional`](@ref) results in [`get`](@ref).
"""
_copy(x) =
    if ismutable(x)
        copy(x)
    else
        x
    end

# a String is mutable in julia
# ismutable("a")
_copy(x::AbstractString) = x

Base.get(parser::Optional, sequence, till, after, i, state::NoMatch) =
    _copy(parser.default)

Base.get(parser::Optional, sequence, till, after, i, state) =
    get(parser.parser,sequence, till, after, i, state)

function Base.get(parser::Repeat, sequence, till, after, i, state::Vector)
    r = Vector{result_type(parser.parser)}(undef,length(state))
    i_ = i
    for (p,s) in enumerate(state)
        after_ = rightof(sequence,i_,parser.parser,s)
        @inbounds r[p] = get(parser.parser, sequence, till, after_, i_, s)
        i_=after_
    end
    r
end

function Base.get(parser::Repeat, sequence, till, after, i, state::Int)
    r = Vector{result_type(parser.parser)}(undef,state)
    i_=i
    s=MatchState()
    for p in 1:state
        after_ = rightof(sequence,i_,parser.parser,s)
        @inbounds r[p] = get(parser.parser,sequence, till, after_, i_, s)
        i_=after_
    end
    r
end

# Base.get(parser::Sequence, sequence, till, after, i, state::MatchState) =
#     get(parser, sequence, till, after, i, ( MatchState() for i in 1:length(parser.parts)) )
# function Base.get(parser::Sequence, sequence, till::Int, after::Int, i::Int, state)
#     r = Vector{Any}(undef,length(parser.parts))
#     i_::Int = i
#     for (p,s) in enumerate(state)
#         after_ = rightof(sequence,i_,parser.parts[p],s)
#         r[p] = get(parser.parts[p],sequence, till, after_, i_, s)
#         i_=after_
#     end
#     1
#     tuple(r...)
# end
@generated function get(parser::Sequence{pts,sts}, sequence, till::Int, after::Int, posi::Int, states) where {pts,sts}
    fpts = fieldtypes(pts)
    spts = Type[ Union{Nothing,state_type(t)} for t in fpts ]
    n = length(fpts)
    subresult = Symbol[ gensym(:r) for i in fpts ]
    part = Symbol[ gensym(:part) for i in fpts ]
    pposi = Symbol[ gensym(:pos) for i in 1:(n+1) ]
    substate = Symbol[ gensym(:s) for i in fpts ]
    init = if states<:MatchState
        [
            quote
            $(substate[i])::MatchState = MatchState()
            @inbounds $(part[i])::$p = parser.parts[$i]
            $(pposi[i])::Int = 0
            end
            for (i,(t,p)) in enumerate(zip(spts,fpts))
        ]
    else
        [
            quote
            $(substate[i])::$t = states[$i]
            @inbounds $(part[i])::$p = parser.parts[$i]
            $(pposi[i])::Int = 0
            end
            for (i,(t,p)) in enumerate(zip(spts,fpts))
        ]
    end
    parseparts = [
        quote
        $(pposi[i+1]) = rightof(sequence,$(pposi[i]),$(part[i]),$(substate[i]))
        $(subresult[i]) = get($(part[i]),sequence, till, $(pposi[i+1]), $(pposi[i]), $(substate[i]))
        end
        for (i,t) in enumerate(fpts)
    ]
    R = quote
        $(pposi[end])::Int = after
        $(init...)
        $(pposi[1]) = posi
        $(parseparts...)
        tuple($(subresult...))
    end
    R
end

## Base.get(x::Palimdrome, str, state) = SubString(str,state...)
## TODO: dispatch Base.get(x::CombinedParser, str, till::Int, after::Int, posi::Int, state) = Base.get(x, str, after, posi, state)
## TODO: dispatch Base.get(x::CombinedParser, str, after::Int, posi::Int, state) = Base.get(x, str, posi, state)
## TODO: dispatch Base.get(x::CombinedParser, str, posi::Int, state) = Base.get(x, str, state)


text:
"""
	grdedit(cmd0::String="", arg1=[], kwargs...)

Reads the header information in a binary 2-D grid file and replaces the information with
values provided on the command line.

Full option list at [`grdedit`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html)

Parameters
----------

- **A** : **adjust** : -- Bool --

    If necessary, adjust the file’s x_inc, y_inc to be compatible with its domain.
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#a)
- **C** : **adjust** : -- Bool --

    Clear the command history from the grid header.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#c)
- **D** : **header** : -- Str --    Flags = [+xxname][+yyname][+zzname][+sscale][+ooffset][+ninvalid][+ttitle][+rremark

    Change these header parameters.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#d)
- **E** : **header** : -- Str --    Flags = [a|h|l|r|t|v]

    Transform the grid in one of six ways and (for l|r|t) interchange the x and y information
    [`-E`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#e)
- **G** : **outgrid** : -- Str --

    Output grid file name. Note that this is optional and to be used only when saving
    the result directly on disk. Otherwise, just use the G = grdedit(....) form.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#g)
- $(GMT.opt_J)
- **N** : **replace** : -- Str or Mx3 array --

    Read the ASCII (or binary) file table and replace the corresponding nodal values in the
    grid with these x,y,z values. Alternatively, provide a Mx3 matrix with values to be changed.
    [`-N`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#n)
- $(GMT.opt_R)
- **S** : **wrap** : -- Bool --

    For global, geographical grids only. Grid values will be shifted longitudinally according to
    the new borders given in ``limits`` (R option).
    [`-S`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#s)
- **T** : **toggle** : -- Bool --

    Make necessary changes in the header to convert a gridline-registered grid to a pixel-registered
    grid, or vice-versa.
    [`-T`](http://gmt.soest.hawaii.edu/doc/latest/grdedit.html#t)
- $(GMT.opt_V)
- $(GMT.opt_bi)
- $(GMT.opt_di)
- $(GMT.opt_e)
- $(GMT.opt_f)
- $(GMT.opt_swap_xy)
"""
function grdedit(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("grdedit", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)

	cmd, = parse_R("", d)
	cmd, = parse_J(cmd, d)
	cmd = parse_V_params(cmd, d)
	cmd, = parse_bi(cmd, d)
	cmd, = parse_di(cmd, d)
	cmd, = parse_e(cmd, d)
	cmd, = parse_f(cmd, d)
	cmd, = parse_swap_xy(cmd, d)

	cmd = add_opt(cmd, 'A', d, [:A :adjust])
	cmd = add_opt(cmd, 'C', d, [:C :clear_history])
	cmd = add_opt(cmd, 'D', d, [:D :header])
	cmd = add_opt(cmd, 'E', d, [:E :flip])
	cmd = add_opt(cmd, 'G', d, [:G :outgrid])
	cmd = add_opt(cmd, 'N', d, [:N :replace])
	cmd = add_opt(cmd, 'S', d, [:S :wrap])
	cmd = add_opt(cmd, 'T', d, [:T :toggle])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "grdedit", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
grdedit(arg1=[], cmd0::String=""; kw...) = grdedit(cmd0, arg1; kw...)


text:
using Symbolics, LinearAlgebra

n = 2
m1 = 1.0
m2 = 1.0
gr = 9.81
l1 = 1.0
l2 = 1.0
k1(q) = [0.5 * l1 * cos(q[1]); 0.5 * l1 * sin(q[1])]
k2(q) = k1(q) + [0.5 * l2 * cos(q[1] + q[2]); 0.5 * l1 * sin(q[1] + q[2])]

@variables q[1:n] q̇[1:n]

# kinematics
k1s = k1(q)
k2s = k2(q)
ks = [k1(q); k2(q)]

# kinematics Jacobians
j1s = Symbolics.jacobian(k1s, q)
j2s = Symbolics.jacobian(k2s, q)
js = Symbolics.jacobian(ks, q)
j1 = eval(Symbolics.build_function(j1s, q)[1])
j2 = eval(Symbolics.build_function(j2s, q)[1])
j = eval(Symbolics.build_function(js, q)[1])

# direct approach
@variables λ1[1:n] λ2[1:n] λ3[1:n] x[1:n]

a = js * λ1
∂a∂q = Symbolics.jacobian(a, q)

# b = transpose(\lamba    ∂a∂q) * λ2
# ∂b∂q = Symbolics.jacobian(b, q)

# mass matrix
mass_matrix_max = Diagonal([m1; m1; m2; m2])
Ms = transpose(js) * mass_matrix_max * js
mass_matrix = eval(Symbolics.build_function(Ms, q)[1])

# mass matrix derivative
Mx = Ms * x
∂Mx∂q = Symbolics.jacobian(Mx, q)
mass_matrix_deriv = eval(Symbolics.build_function(∂Mx∂q, q, x)[1])

# dynamics bias
Cs = ∂Mx∂q * x
dynamics_bias = eval(Symbolics.build_function(Cs, q, x)[1])

# dynamics bias derivative
∂C∂q̇ = 2.0 * transpose(js) * mass_matrix_max * ∂a∂q
∂C∂q = transpose(∂a∂q) * mass_matrix_max * ∂a∂q

# differentiation approach
function lagrangian(q, q̇)
    L = 0.0

    v1 = j1(q) * q̇
    v2 = j2(q) * q̇

    L += 0.5 * m1 * transpose(v1) * v1
    # L -= m1 * gr * k1(q)[2]

    L += 0.5 * m2 * transpose(v2) * v2
    # L -= m2 * gr * k2(q)[2]

    return L
end

L = lagrangian(q, q̇)

ddL = Symbolics.hessian(L, [q; q̇])
dLq = Symbolics.gradient(L, q)
ddLq̇q = ddL[n .+ (1:n), 1:n]

M = ddL[n .+ (1:n), n .+ (1:n)]
C = ddLq̇q * q̇ - dLq

mass_matrix_lag = eval(Symbolics.build_function(M, q)[1])
dynamics_bias_lag = eval(Symbolics.build_function(C, q, q̇)[1])

# test
q0 = rand(n)
q̇0 = rand(n)

norm(mass_matrix(q0) - mass_matrix_lag(q0), Inf)
norm(dynamics_bias(q0, q̇0) - dynamics_bias_lag(q0, q̇0), Inf)

dynamics_bias(q0, q̇0)
dynamics_bias_lag(q0, q̇0)


text:
# Convenience function to allow for things like Array(tp) or CuArray(tp)
# Not sure if this counts as type piracy...
(::Type{T})(p::TensorProto) where T = array(p) |> T
(::Type{Ref{T}})(p::TensorProto) where T = array(p) |> T |> Ref

"""
    array(p::TensorProto)

Return `p` as an reshaped and reinterpreted array.
"""
function array(p::TensorProto)
    # Copy pasted from jl
    # Can probably be cleaned up a bit
    # TODO: Add missing datatypes...
    if p.data_type === TensorProto_DataType.INT64
        if isdefined(p, :int64_data) && !isempty(p.int64_data)
            return reshape(reinterpret(Int64, p.int64_data), reverse(p.dims)...)
        end
        return reshape(reinterpret(Int64, p.raw_data), reverse(p.dims)...)
    end

    if p.data_type === TensorProto_DataType.INT32
        if isdefined(p, :int32_data) && !isempty(p.int32_data)
            return reshape(p.int32_data , reverse(p.dims)...)
        end
        return reshape(reinterpret(Int32, p.raw_data), reverse(p.dims)...)
    end

    if p.data_type === TensorProto_DataType.INT8
        return reshape(reinterpret(Int8, p.raw_data), reverse(p.dims)...)
    end

    if p.data_type === TensorProto_DataType.DOUBLE
        if isdefined(p, :double_data) && !isempty(p.double_data)
            return reshape(p.double_data , reverse(p.dims)...)
        end
        return reshape(reinterpret(Float64, p.raw_data), reverse(p.dims)...)
    end

    if p.data_type === TensorProto_DataType.FLOAT
        if isdefined(p,:float_data) && !isempty(p.float_data)
            return reshape(reinterpret(Float32, p.float_data), reverse(p.dims)...)
        end
        return reshape(reinterpret(Float32, p.raw_data), reverse(p.dims)...)
    end

    if p.data_type === TensorProto_DataType.FLOAT16
        return reshape(reinterpret(Float16, p.raw_data), reverse(p.dims)...)
    end
end

Base.size(vip::ValueInfoProto) = size(vip._type)
Base.size(tp::TypeProto) = size(tp.tensor_type)
Base.size(tp_t::TypeProto_Tensor) = size(tp_t.shape)
Base.size(tsp::TensorShapeProto) = size.(Tuple(reverse(tsp.dim)))
Base.size(tsp_d::TensorShapeProto_Dimension) = isdefined(tsp_d, :dim_value) ? tsp_d.dim_value : missing

"""
    attribute(p::AttributeProto)

Return attribute in `p` as a name => value pair.
"""
function attribute(p::AttributeProto)
    # Copy paste from ONNX.jl
    if (p._type != 0)
        field = [:f, :i, :s, :t, :g, :floats, :ints, :strings, :tensors, :graphs][p._type]
        if field === :s
            return Symbol(p.name) => String(getfield(p, field))
        elseif  field === :strings
            return Symbol(p.name) => String.(getfield(p, field))
        end
        return Symbol(p.name) => getfield(p, field)
    end
end

Base.Dict(pa::AbstractVector{AttributeProto}) = Dict(attribute(p) for p in pa)


text:
"""
    rect(xmin, ymin, w, h; action=:none)
    rect(xmin, ymin, w, h, action)

Create a rectangle with one corner at (`xmin`/`ymin`) with
width `w` and height `h`, and add it to the current path.
Then apply `action`.

Returns a tuple of two points, the corners of a bounding box
that encloses the rectangle.

See `box()` for more ways to do similar things, such as
supplying two opposite corners, placing by centerpoint and
dimensions.
"""
function rect(xmin::Real, ymin::Real, w::Real, h::Real;
        action = :none)
    if action != :path
        newpath()
    end
    Cairo.rectangle(get_current_cr(), xmin, ymin, w, h)
    do_action(action)
    return Point(xmin, ymin), Point(xmin + w, ymin + h)
end

rect(xmin::Real, ymin::Real, w::Real, h::Real, action::Symbol) = rect(xmin, ymin, w, h, action = action)

"""
    rect(cornerpoint, w, h; action = none, reversepath=false,
        vertices=false)
    rect(cornerpoint, w, h, action; reversepath=false,
        vertices=false)

Create a rectangle with one corner at `cornerpoint` with
width `w` and height `h`, and add it to the current path.
Then apply `action`.

Use `vertices=true` to return an array of the four corner
points: bottom left, top left, top right, bottom right.

`reversepath` reverses the direction of the path (and
returns points in the order: bottom left, bottom right, top
right, top left).

Returns the four corner vertices.
"""
function rect(cornerpoint::Point, w::Real, h::Real;
        action = :none,
        reversepath = false,
        vertices = false)

    pts =  [
        Point(cornerpoint.x,     cornerpoint.y + h),
        Point(cornerpoint.x,     cornerpoint.y),
        Point(cornerpoint.x + w, cornerpoint.y),
        Point(cornerpoint.x + w, cornerpoint.y + h)
        ]
    if reversepath == true
        pts = pts[[1, 4, 3, 2]]
    end
    if vertices == false && action != :none
        move(pts[1])
        line(pts[2])
        line(pts[3])
        line(pts[4])
        closepath()
        do_action(action)
    end
    return pts
end

rect(cornerpoint::Point, w::Real, h::Real, action::Symbol;
        reversepath = false,
        vertices = false) = rect(cornerpoint, w, h;
                action = action,
                reversepath = reversepath,
                vertices = vertices)

"""
    box(cornerpoint1, cornerpoint2; action=:none, vertices=false, reversepath=false)
    box(cornerpoint1, cornerpoint2, action; vertices=false, reversepath=false)

Create a box (rectangle) between two points and add it to
the current path. Then apply `action`.

Use `vertices=true` to return an array of the four corner
points: bottom left, top left, top right, bottom right
rather than execute action.

`reversepath` reverses the direction of the path (and
returns points in the order: bottom left, bottom right, top
right, top left).
"""
function box(corner1::Point, corner2::Point;
        action::Symbol=:none,
        reversepath=false,
        vertices=false)
    # rearrange to get topleft -> bottom right
    # I'm not sure whether this is worth doing just
    # to get the order of vertices correct
    xmin, xmax = extrema((corner1.x, corner2.x))
    ymin, ymax = extrema((corner1.y, corner2.y))

    c1 = Point(xmin, ymin)
    c2 = Point(xmax, ymax)
    pts = [
        Point(c1.x, c2.y),
        Point(c1.x, c1.y),
        Point(c2.x, c1.y),
        Point(c2.x, c2.y)
        ]
    if reversepath == true
        pts = pts[[1, 4, 3, 2]]
    end
    if vertices == false && action != :none
        move(pts[1])
        line(pts[2])
        line(pts[3])
        line(pts[4])
        closepath()
        do_action(action)
    end
    return pts
end

box(corner1::Point, corner2::Point, action::Symbol;
        reversepath=false,
        vertices=false) = box(corner1, corner2;
                action=action,
                reversepath=reversepath,
                vertices=vertices)

"""
    box(points::Array; action=:none,
        reversepath=reversepath,
        vertices=vertices)
    box(points::Array; action=:none,
        reversepath=reversepath,
        vertices=vertices)

Create a box/rectangle using the first two points of an
array of Points to defined opposite corners, and add it to
the current path. Then apply `action`.

Use `vertices=true` to return an array of the four corner
points: bottom left, top left, top right, bottom right
rather than execute action.
"""
box(bbox::Array, action::Symbol; kwargs...) =
    box(bbox[1], bbox[2], action; kwargs...)

box(bbox::Array; kwargs...) =
    box(bbox[1], bbox[2]; kwargs...)

"""
    box(pt::Point, width, height; action=:none, vertices=false)
    box(pt::Point, width, height, action=:none; vertices=false)

Create a box/rectangle centered at point `pt` with width and
height. Use `vertices=true` to return an array of the four
corner points rather than apply the action.

`reversepath` reverses the direction of the path.
"""
function box(pt::Point, width, height;
        action = :none,
        reversepath = false,
        vertices = false)
    pts = [
        Point(pt.x - width/2, pt.y + height/2),
        Point(pt.x - width/2, pt.y - height/2),
        Point(pt.x + width/2, pt.y - height/2),
        Point(pt.x + width/2, pt.y + height/2)
        ]
    if reversepath == true
        pts = pts[[1, 4, 3, 2]]
    end
    if vertices == false && action != :none
        move(pts[1])
        line(pts[2])
        line(pts[3])
        line(pts[4])
        closepath()
        do_action(action)
    end
    return pts
end

box(pt::Point, width, height, action::Symbol;
        reversepath = false,
        vertices = false) = box(pt, width, height;
                action = action,
                reversepath = reversepath,
                vertices = vertices)

box(x::Real, y::Real, width::Real, height::Real, action::Symbol) =
    rect(x - width/2.0, y - height/2.0, width, height, action)

box(x::Real, y::Real, width::Real, height::Real; action::Symbol=:none) =
    rect(x - width/2.0, y - height/2.0, width, height, action)

"""
    box(pt, width, height, cornerradius, action=:none)
    box(pt, width, height, cornerradius; action=:none)

Draw a box/rectangle centered at point `pt` with `width` and `height` and
round each corner by `cornerradius`.

The constructed path consists of arcs and straight lines.
"""
box(centerpoint::Point, width, height, cornerradius::Real, action::Symbol) =
    box(centerpoint, width, height, fill(cornerradius, 4), action=action)

box(centerpoint::Point, width, height, cornerradius::Real; action=:none) =
    box(centerpoint, width, height, fill(cornerradius, 4), action=action)

"""
    box(pt, width, height, cornerradii::Array; action=:none)
    box(pt, width, height, cornerradii::Array, action=:none)

Draw a box/rectangle centered at point `pt` with `width` and `height` and
round each corner by the corresponding value in the array `cornerradii`.

The constructed path consists of arcs and straight lines.

The first corner is the one at the bottom left, the second at
the top left, and so on.

## Example

```
@draw begin
    box(O, 120, 120, [0, 20, 40, 60], :fill)
end
```
"""
function box(centerpoint::Point, width, height, cornerradii::Array;
        action=:none)
    gsave()
    translate(centerpoint)

    length(cornerradii) != 4 && throw(error("box() must have four values to specify rounded corners"))

    # bottom left
    p1start  = Point(O.x - width/2 + cornerradii[1], O.y + height/2)
    p1center = Point(O.x - width/2 + cornerradii[1], O.y + height/2 - cornerradii[1])
    p1end    = Point(O.x - width/2, O.y + height/2 - cornerradii[1])

    # top left
    p2start  = Point(O.x - width/2, O.y - height/2 + cornerradii[2])
    p2center = Point(O.x - width/2 + cornerradii[2], O.y - height/2 + cornerradii[2])
    p2end    = Point(O.x - width/2 + cornerradii[2], O.y - height/2)

    # top right
    p3start  = Point(O.x + width/2 - cornerradii[3], O.y - height/2)
    p3center = Point(O.x + width/2 - cornerradii[3], O.y - height/2 + cornerradii[3])
    p3end    = Point(O.x + width/2, O.y - height/2 + cornerradii[3])

    # bottom right
    p4start  = Point(O.x + width/2, O.y + height/2 - cornerradii[4])
    p4center = Point(O.x + width/2 - cornerradii[4], O.y + height/2 - cornerradii[4])
    p4end    = Point(O.x + width/2 - cornerradii[4], O.y + height/2)

    #  start at bottom center then bottomleft→topleft→topright→bottomright
    newpath()
    move(Point(O.x, O.y + height/2))

    line(p1start)
    arc(p1center, cornerradii[1], π/2, π, :none)
    line(p1end)

    line(p2start)
    arc(p2center, cornerradii[2], π, (3π)/2, :none)
    line(p2end)

    line(p3start)
    arc(p3center, cornerradii[3], 3π/2, 0, :none)
    line(p3end)

    line(p4start)
    arc(p4center, cornerradii[4], 0, π/2, :none)
    line(p4end)

    closepath()
    grestore()
    do_action(action)
    return Point(centerpoint.x - width/2, centerpoint.y - height/2),
           Point(centerpoint.x + width/2, centerpoint.y + height/2)
end

box(centerpoint::Point, width, height, cornerradii::Array, action::Symbol) =
     box(centerpoint::Point, width, height, cornerradii::Array, action=action)

"""
    ngon(x, y, radius, sides=5, orientation=0;
        action = :none,
        vertices = false,
        reversepath = false)

Draw a regular polygon centered at point `Point(x,y)`.
"""
function ngon(x::Real, y::Real, radius::Real, sides::Int=5, orientation=0.0;
            action=:none,
            vertices=false,
            reversepath=false)
    ptlist = [Point(x+cos(orientation + n * 2pi/sides) * radius,
                    y+sin(orientation + n * 2pi/sides) * radius) for n in 1:sides]
    if !vertices
        poly(ptlist, action=action, close=true, reversepath=reversepath)
    end
    if reversepath
        return reverse(ptlist)
    else
        return ptlist
    end
end

ngon(x::Real, y::Real, radius::Real, sides::Int, action::Symbol) =
    ngon(x, y, radius, sides, action=action, vertices=false, reversepath=false)

ngon(x::Real, y::Real, radius::Real, sides::Int, orientation::Real, action::Symbol) =
    ngon(x, y, radius, sides, orientation, action=action, vertices=false, reversepath=false)

# Point version
"""
    ngon(centerpos, radius, sides=5, orientation=0;
        action=:none,
        vertices=false,
        reversepath=false)

Draw a regular polygon centered at point `centerpos`.

Find the vertices of a regular n-sided polygon centered at `x`, `y` with
circumradius `radius`.

The polygon is constructed counterclockwise, starting with the first vertex
drawn below the positive x-axis.

If you just want the raw points, use keyword argument `vertices=true`, which
returns the array of points. Compare:

```julia
ngon(0, 0, 4, 4, 0, vertices=true) # returns the polygon's points:

4-element Array{Luxor.Point, 1}:
Luxor.Point(2.4492935982947064e-16, 4.0)
Luxor.Point(-4.0, 4.898587196589413e-16)
Luxor.Point(-7.347880794884119e-16, -4.0)
Luxor.Point(4.0, -9.797174393178826e-16)
```

whereas

```
ngon(0, 0, 4, 4, 0, :close) # draws a polygon
```
"""
ngon(centerpoint::Point, radius, sides::Int=5, orientation=0.0;
        action=:none, vertices=false, reversepath=false) = ngon(centerpoint.x, centerpoint.y, radius, sides, orientation,
        action=action, vertices=vertices, reversepath=reversepath)

# action as argument

# this is a bit untidy?
ngon(centerpoint::Point, radius::Real, sides::Int, orientation::Real, a::Symbol;
        action=:none, vertices=false, reversepath=false) = ngon(centerpoint, radius, sides, orientation,
        action=a, vertices=vertices, reversepath=reversepath)

"""
    ngonside(centerpoint::Point, sidelength::Real, sides::Int=5, orientation=0;
            action=:none,
            vertices=false,
            reversepath=false)
    ngonside(centerpoint::Point, sidelength::Real, sides::Int, orientation, action;
            vertices=false,
            reversepath=false)

Draw a regular polygon centered at `centerpoint` with `sides` sides of length `sidelength`.
"""
function ngonside(centerpoint::Point, sidelength::Real, sides::Int=5, orientation=0;
        action=:none,
        vertices=false,
        reversepath=false)
    radius = 0.5 * sidelength * csc(pi/sides)
    ngon(centerpoint, radius, sides, orientation, action=action,
        vertices=vertices,
        reversepath=reversepath)
end

ngonside(centerpoint::Point, sidelength::Real, sides::Int, orientation::Real, action::Symbol;
        vertices=false,
        reversepath=false) = ngonside(centerpoint, sidelength, sides, orientation, action=action,
        vertices=vertices, reversepath=reversepath)

ngonside(centerpoint::Point, sidelength::Real, sides::Int, action::Symbol;
        vertices=false,
        reversepath=false) = ngonside(centerpoint, sidelength, sides, action=action,
        vertices=vertices, reversepath=reversepath)

function star(x::Real, y::Real, radius::Real, npoints::Int=5, ratio::Real=0.5, orientation=0;
        action=:none,
        vertices = false,
        reversepath=false)
    outerpoints = [Point(x+cos(orientation + n * 2pi/npoints) * radius,
                         y+sin(orientation + n * 2pi/npoints) * radius) for n in 1:npoints]
    innerpoints = [Point(x+cos(orientation + (n + 1/2) * 2pi/npoints) * (radius * ratio),
                         y+sin(orientation + (n + 1/2) * 2pi/npoints) * (radius * ratio)) for n in 1:npoints]
    result = Point[]
    for i in eachindex(outerpoints)
        push!(result, outerpoints[i])
        push!(result, innerpoints[i])
    end
    if reversepath
        result = reverse(result)
    end
    if !vertices
        poly(result, action=action, close=true)
    end
    return result
end

"""
    star(center, radius, npoints, ratio=0.5, orientation, action=:none;
        vertices = false, reversepath=false)
    star(center, radius, npoints, ratio=0.5, orientation=0.0;
        action=:none, vertices = false, reversepath=false)

Make a star centered at `center` with `npoints` sections oriented by `orientation`.
`ratio` specifies the height of the smaller radius of the
star relative to the larger.

Returns the vertices of the star.

Use `vertices=true` to only return the vertices of a star
instead of making it.

## Examples

```julia
star(O, 120, 5, 0.5, 0.0, :fill,
    vertices = false,
    reversepath=false)

star(O, 220, 5, 0.5;
    action=:stroke,
    vertices = false,
    reversepath=false)
```
"""
function star(centerpoint::Point, radius::Real, npoints::Int=5, ratio::Real=0.5, orientation=0;
              action=:none,
              vertices=false,
              reversepath=false)
    star(centerpoint.x, centerpoint.y, radius, npoints, ratio, orientation, action=action,
         vertices = vertices, reversepath=reversepath)
end

star(x::Real, y::Real, radius::Real, npoints::Int, ratio::Real, orientation, a::Symbol;
    action=:none,
    vertices = false,
    reversepath=false) = star(Point(x, y), radius, npoints, ratio, orientation,
    action=a, vertices = vertices, reversepath=reversepath)

star(pt::Point, radius::Real, npoints::Int, ratio::Real, orientation, a::Symbol;
    action=:none,
    vertices = false,
    reversepath=false) = star(pt, radius, npoints, ratio, orientation,
    action=a, vertices = vertices, reversepath=reversepath)

"""
    cropmarks(center, width, height)

Draw cropmarks (also known as trim marks). Use current color.
"""
function cropmarks(center, width, height)
    gap = 5
    crop = 15
    gsave()
    setline(0.5)
    setdash("solid")
    # horizontal top left
    line(Point(-width/2 - gap - crop, -height/2),
         Point(-width/2 - gap, -height/2),
         :stroke)

    # horizontal bottom left
    line(Point(-width/2 - gap - crop, height/2),
         Point(-width/2 - gap, height/2),
         :stroke)

    # horizontal top right
    line(Point(width/2 + gap, -height/2),
         Point(width/2 + gap + crop, -height/2),
         :stroke)

    # horizontal bottom right
    line(Point(width/2 + gap, height/2),
         Point(width/2 + gap + crop, height/2),
         :stroke)

    # vertical top left
    line(Point(-width/2, -height/2 - gap - crop),
         Point(-width/2, -height/2 - gap),
         :stroke)

    # vertical bottom left
    line(Point(-width/2, height/2 + gap),
         Point(-width/2, height/2 + gap + crop),
         :stroke)

    # vertical top right
    line(Point(width/2, -height/2 - gap - crop),
         Point(width/2, -height/2 - gap),
         :stroke)

    # vertical bottom right
    line(Point(width/2, height/2 + gap),
         Point(width/2, height/2 + gap + crop),
         :stroke)

    grestore()
end

"""
    polycross(pt::Point, radius, npoints::Int, ratio=0.5, orientation=0.0;
        action      = :none,
        splay       = 0.5,
        vertices    = false,
        reversepath = false)
    polycross(pt::Point, radius, npoints::Int, ratio=0.5, orientation=0.0, action;
        splay       = 0.5,
        vertices    = false,
        reversepath = false)

Make a cross-shaped polygon with `npoints` arms to fit
inside a circle of radius `radius` centered at `pt`.

`ratio` specifies the ratio of the two sides of each arm.
`splay` makes the arms ... splayed.

Use `vertices=true` to return the vertices of the shape
instead of executing the action.

(Adapted from Compose.jl.xgon()))

## Examples

```julia
polycross(O, 100, 5,
    action = :fill,
    splay  = 0.5)

polycross(O, 120, 5, 0.5, 0.0, :stroke,
    splay  = 0.5)
```
"""
function polycross(pt::Point, radius::Real, npoints::Int, ratio=0.5, orientation=0.0;
        action      = :none,
        splay       = .5,
        vertices    = false,
        reversepath = false)
    # adapted from:    Compose.jl, https://github.com/GiovineItalia/Compose.jl/src/form.jl
    # original author: mattriks

    ratio = clamp(ratio, 0.0, 1.0)

    θ₁ = range(π/2 + orientation + 0, stop = π/2 + orientation + 2π, length = npoints + 1)[1:end-1]

    width = 2radius * ratio * sin(π/npoints)

    dₒ = asin(clamp(splay * width/radius, -1.0, 1.0))
    dᵢ = asin(0.5 * width/(radius * ratio))

    r₂ = repeat([radius * ratio, radius, radius], outer = npoints)
    θ₂ = vec([mod2pi(θ + x) for x in [-dᵢ, -dₒ, dₒ], θ in θ₁])

    pts = @. Point.(pt.x .+ r₂ .* cos.(θ₂), pt.y .+ r₂ .* sin.(θ₂))
    if reversepath
        reverse!(pts)
    end
    if !vertices
        poly(pts, action=action, close=true)
    end
    return pts
end

polycross(pt::Point, radius::Real, npoints::Int, ratio, orientation::Real, action::Symbol;
    splay       = .5,
    vertices    = false,
    reversepath = false) = polycross(pt, radius, npoints, ratio, orientation,
        action = action,
        splay       = splay,
        vertices    = vertices,
        reversepath = reversepath)


text:
@testset "send_retx.jl" begin
    @testset "Retx SYN twice at the right times, then ack" begin
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        tick!(sender, rto - 1)
        expect_no_seg(sender)
        tick!(sender, 1)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_state(sender, JLSponge.SYN_SENT)
        @test sender.bytes_in_flight == 1

        tick!(sender, 2rto - 1)
        expect_no_seg(sender)
        tick!(sender, 1)
        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_state(sender, JLSponge.SYN_SENT)
        @test sender.bytes_in_flight == 1
        ack_received_test(sender, WrappingInt32(isn+1))
        @test sender.bytes_in_flight == 0
        tick!(sender, 2rto - 1)
        @test sender.ticks == 0
    end

    @testset "Retx SYN until too many retransmissions" begin
        MAX_RETX_ATTEMPTS = 10
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        expect_state(sender, JLSponge.SYN_SENT)

        for i in 0:MAX_RETX_ATTEMPTS-1
            tick!(sender, (Int(rto) << i) - 1)
            @test consecutive_retransmissions(sender) <= MAX_RETX_ATTEMPTS
            tick!(sender, 1)
            expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
            expect_no_seg(sender)
            expect_state(sender, JLSponge.SYN_SENT)
            @test sender.bytes_in_flight == 1
        end
        tick!(sender, (Int(rto) << MAX_RETX_ATTEMPTS) - 1)
        @test sender.consecutive_retransmissions <= MAX_RETX_ATTEMPTS
        tick!(sender, 1)
        @test sender.consecutive_retransmissions > MAX_RETX_ATTEMPTS
    end

    @testset "Send some data, the retx and succeed, then retx till limit" begin
        MAX_RETX_ATTEMPTS = 10
        isn = rand(UInt32(0):typemax(UInt32))
        rto = rand(UInt16(30):UInt16(10000))
        window_size = rand(UInt16(50000):UInt16(63000))

        sender = TCPSender(fixed_isn=WrappingInt32(isn), retx_timeout = rto)
        fill_window!(sender)

        expect_seg(sender, no_flag=true, syn=true, payload_size=0, seqno=WrappingInt32(isn))
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn+1))
        expect_state(sender, JLSponge.SYN_ACKED)
        write_bytes!(sender, "abcd")
        expect_seg(sender, payload_size=4)
        expect_no_seg(sender)
        ack_received_test(sender, WrappingInt32(isn+5))
        @test sender.bytes_in_flight == 0
        write_bytes!(sender, "efgh")
        expect_seg(sender, payload_size=4)
        expect_no_seg(sender)
        @test sender.consecutive_retransmissions <= MAX_RETX_ATTEMPTS
        ack_received_test(sender, WrappingInt32(isn+9))
        @test sender.bytes_in_flight == 0

        write_bytes!(sender, "efgh")
        expect_seg(sender, payload_size=4, seqno=WrappingInt32(isn+9))
        expect_no_seg(sender)

        for i in 0:MAX_RETX_ATTEMPTS-1
            tick!(sender, (Int(rto) << i) - 1)
            @test sender.consecutive_retransmissions <= MAX_RETX_ATTEMPTS
            expect_no_seg(sender)
            tick!(sender, 1)
            expect_seg(sender, payload_size=4, seqno=WrappingInt32(isn + 9))
            expect_no_seg(sender)
            expect_state(sender, JLSponge.SYN_ACKED)
            @test sender.bytes_in_flight == 4
        end

        tick!(sender, (Int(rto) << MAX_RETX_ATTEMPTS) - 1)
        @test sender.consecutive_retransmissions <= MAX_RETX_ATTEMPTS
        tick!(sender, 1)
        @test sender.consecutive_retransmissions > MAX_RETX_ATTEMPTS
    end
end

max_stars_repo_name: AquaIndigo/JLSponge

max_issues_repo_name: AquaIndigo/JLSponge

max_forks_repo_name: AquaIndigo/JLSponge


text:
using DSGE, Test

# Initialize model object
m = AnSchorfheide(testing = true)
m <= Setting(:forecast_horizons, 12)

# Define dummy scenario
scen = Scenario(:testscen, "Test Scenario", [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF")

# Check compute_scenario_system throws error before zeroing out measurement error
@testset "Check compute_scenario_system error throwing" begin
    @test_throws ErrorException compute_scenario_system(m, scen)
end

# Zero out measurement error (or else filtering shocks won't work exactly)
for para in [:e_y, :e_π, :e_R]
    m[para].value = 0
end

# Test compute_scenario_system
@testset "Test compute_scenario_system, filter_shocks!, and shock_scaling" begin
    global sys = compute_scenario_system(m, scen)
    @test all(x -> x == 0, sys[:CCC])
    @test all(x -> x == 0, sys[:DD])
    @test all(x -> x == 0, sys[:DD_pseudo])
    for shock in keys(m.exogenous_shocks)
        ind = m.exogenous_shocks[shock]
        if shock in scen.instrument_names
            @test sys[:QQ][ind, ind] != 0
        else
            @test sys[:QQ][ind, ind] == 0
        end
    end

    # Test filter_shocks!
    for var in scen.target_names
        scen.targets[var] = rand(forecast_horizons(m))
    end
    global forecastshocks = filter_shocks!(m, scen, sys)
    for var in scen.instrument_names
        ind = m.exogenous_shocks[var]
        @test scen.instruments[var] == forecastshocks[ind, :]
    end

    s_T = zeros(n_states_augmented(m))
    global _, forecastobs,_ ,_  = forecast(m, sys, s_T, shocks = forecastshocks)
    for var in scen.target_names
        ind = m.observables[var]
        @test scen.targets[var] ≈ forecastobs[ind, :]
    end

    # Test scenario with shock scaling
    global scale = Scenario(:scaledscen, "Test Shock Scaling Scenario",
                            [:obs_gdp, :obs_cpi], [:g_sh, :rm_sh], "REF",
                            shock_scaling = 2.0)
    scale.targets = scen.targets

    global forecastshocks = filter_shocks!(m, scale, sys)
    for var in scale.target_names
        ind = m.observables[var]
        @test scale.targets[var] == scen.targets[var]
    end
    for var in scale.instrument_names
        ind = m.exogenous_shocks[var]
        @test scale.instruments[var] == forecastshocks[ind, :]
    end
end


text:
module TestChannels

using Reagents
using Reagents.Internal: @trace
using Test
using ..Utils: concurrently

function test_simple_swap()
    @testset for spawn in [false, true]
        test_simple_swap(spawn)
    end
end

function test_simple_swap(spawn::Bool)
    c1, c2 = Reagents.channel(Int, Char)
    local y1, y2
    concurrently(function f1()
        y1 = c1(1)
    end, function f2()
        y2 = c2('2')
    end; spawn = spawn)
    @test (y1, y2) == ('2', 1)
end

function test_multi_swap()
    @testset for spawn in [false, true]
        test_multi_swap(spawn)
    end
end

function test_multi_swap(spawn::Bool, nrepeat = spawn ? 1000 : 1)
    c1, c2 = Reagents.channel(Int)
    c3, c4 = Reagents.channel(Int)
    local y1, y2, y3
    ok1 = ok2 = ok3 = true
    failed = Threads.Atomic{Bool}(false)
    concurrently(
        function f1()
            for i in 1:nrepeat
                failed[] && break
                @trace label = :f1 i taskid = objectid(current_task())
                y1 = (c1 ⨟ c3)(1 + 100i)
                ok1 &= y1 == 3 + 100i
                if !ok1
                    failed[] = true
                    error("f1 failed: got $(y1)")
                    break
                end
            end
        end,
        function f2()
            for i in 1:nrepeat
                failed[] && break
                @trace label = :f2 i taskid = objectid(current_task())
                y2 = c2(2 + 100i)
                ok2 &= y2 == 1 + 100i
                if !ok1
                    failed[] = true
                    error("f2 failed: got $(y2)")
                    break
                end
            end
        end,
        function f3()
            for i in 1:nrepeat
                failed[] && break
                @trace label = :f3 i taskid = objectid(current_task())
                y3 = c4(3 + 100i)
                ok3 &= y3 == 2 + 100i
                if !ok1
                    failed[] = true
                    error("f3 failed: got $(y3)")
                    break
                end
            end
        end;
        spawn = spawn,
    )
    @test ok1 == ok2 == ok3 == true
    @test (y1, y2, y3) == (3, 1, 2) .+ 100nrepeat
end

function test_simple_choice()
    @testset for spawn in [false, true]
        test_simple_choice(spawn)
    end
end

function test_simple_choice(spawn::Bool, nrepeat = spawn ? 1000 : 1)
    a2b, b2a = Reagents.channel(Int, Nothing)
    c2d, d2c = Reagents.channel(Int, Nothing)
    ok = true
    concurrently(
        function receiver()
            for i in 1:nrepeat
                y = (b2a | d2c)()
                ok &= y == i
                ok || error("receiver failed: got $y")
            end
        end,
        function sender()
            for i in 1:nrepeat
                if rand() > 0.5
                    a2b(i)
                else
                    c2d(i)
                end
            end
        end,
        spawn = spawn,
    )
    @test ok
end

end  # module


text:
# Interactive plots with Observables

"""
    vlslice(meta, var; normal=:y, axisunit=SI, comp=0)

Interactive 2D slice of 3D `var` in `normal` direction.
"""
function vlslice(meta::MetaVLSV, var; normal=:y, axisunit=SI, comp=0)
   dir, str1, str2 =
      if normal == :x
         1, "y", "z"
      elseif normal == :y
         2, "x", "z"
      else
         3, "x", "y"
      end

   unitx = axisunit == EARTH ? " [Re]" : " [m]"

   dx = meta.dcoord[dir] / 2^meta.maxamr

   pArgs = Vlasiator.set_args(meta, var, axisunit; normal, origin=0.0)
   x, y = Vlasiator.get_axis(axisunit, pArgs.plotrange, pArgs.sizes)

   nsize = meta.ncells[dir]
   depth = nsize*2^meta.maxamr

   fig = Figure()
   ax = Axis(fig[1, 1], aspect=DataAspect())
   ax.xlabel = str1*unitx
   ax.ylabel = str2*unitx

   lsgrid = SliderGrid(fig[2, 1],
      (label = "location in normal direction $(String(normal))",
       range=1:depth,
       format=x -> "$(x) cells"),
   )

   sliderobservables = [s.value for s in lsgrid.sliders]

   slice = lift(sliderobservables...) do slvalues...
      begin
         origin = (slvalues[1]-1)*dx + meta.coordmin[dir]
         pArgs = Vlasiator.set_args(meta, var, axisunit; normal, origin)
         Vlasiator.prep2dslice(meta, var, normal, comp, pArgs)
      end
   end

   heatmap!(ax, x, y, slice, colormap=:turbo)

   set_close_to!(lsgrid.sliders[1], .5depth)

   fig
end

"""
    vdfslices(meta, location; fmin=1f-16, species="proton", unit=SI, verbose=false)

Three orthogonal slices of VDFs from `meta` at `location`.
# Optional Arguments
- `fmin`: minimum VDF threshold for plotting.
- `species`: name of particle.
- `unit`: unit of input `location`, `SI` or `EARTH`.
"""
function vdfslices(meta, location; fmin=1f-16, species="proton", unit=SI, verbose=false)
   if haskey(meta.meshes, species)
      vmesh = meta.meshes[species]
   else
      throw(ArgumentError("Unable to detect population $species"))
   end

   unit == EARTH && (location .*= Vlasiator.RE)

   # Calculate cell ID from given coordinates
   cidReq = getcell(meta, location)
   cidNearest = getnearestcellwithvdf(meta, cidReq)

   cellused = getcellcoordinates(meta, cidNearest)

   if verbose
      @info "Original coordinates : $location"
      @info "Original cell        : $(getcellcoordinates(meta, cidReq))"
      @info "Nearest cell with VDF: $cellused"
      let
         x, y, z = getcellcoordinates(meta, cidNearest)
         @info "cellid $cidNearest, x = $x, y = $y, z = $z"
      end
   end

   vcellids, vcellf = readvcells(meta, cidNearest; species)

   f = Vlasiator.reconstruct(vmesh, vcellids, vcellf)

   fig = Figure()
   ax = Axis3(fig[1, 1], aspect=(1,1,1), title = "VDF at $cellused in log scale")
   ax.xlabel = "vx [m/s]"
   ax.ylabel = "vy [m/s]"
   ax.zlabel = "vz [m/s]"

   x = LinRange(vmesh.vmin[1], vmesh.vmax[1], vmesh.vblock_size[1]*vmesh.vblocks[1])
   y = LinRange(vmesh.vmin[2], vmesh.vmax[2], vmesh.vblock_size[2]*vmesh.vblocks[2])
   z = LinRange(vmesh.vmin[3], vmesh.vmax[3], vmesh.vblock_size[3]*vmesh.vblocks[3])

   lsgrid = SliderGrid(fig[2, 1],
      (label="vx", range=1:length(x), format=i -> "$(round(x[i], digits=2))"),
      (label="vy", range=1:length(y), format=i -> "$(round(y[i], digits=2))"),
      (label="vz", range=1:length(z), format=i -> "$(round(z[i], digits=2))"),
   )

   for i in eachindex(f)
      if f[i] < fmin f[i] = fmin end
   end
   data = log10.(f)

   plt = volumeslices!(ax, x, y, z, data, colormap=:viridis)
   #TODO: improve on colormap setup!
   cbar = Colorbar(fig, plt,
      label="f(v)",
      minorticksvisible=true)

   fig[1, 2] = cbar

   # connect sliders to volumeslices update methods
   sl_yz, sl_xz, sl_xy = lsgrid.sliders

   on(sl_yz.value) do v; plt[:update_yz][](v) end
   on(sl_xz.value) do v; plt[:update_xz][](v) end
   on(sl_xy.value) do v; plt[:update_xy][](v) end

   set_close_to!(sl_yz, .5length(x))
   set_close_to!(sl_xz, .5length(y))
   set_close_to!(sl_xy, .5length(z))

   fig
end


text:
"""
    setMutationOperator(method::Function, varArgs...; genj::GeneticJuliaStructure = GenJ,
                        probability::AbstractFloat = 0.3)

Sets the mutation operator of parents for `genj`, receiving a function that
mutates an individual and returns the genotype of the mutated child. This function
constructs a structure of type `MutationOperator` and adds it to the main structure.

# Arguments
- `method::Function`: function that receives, first, a parent as argument which
    is an individual's genotype and `"varArgs"`, which is a list of arguments,
    as own arguments, and returns the genotype of the mutated child.
    (For information about representation look [`Genotype`](@ref)). See `"Examples"`.
- `varArgs`: list of arguments for the mutation method. See `"Examples"`.

# Keyword Arguments
- `genj::GeneticJuliaStructure = GenJ`: the main structure. For `code user`: don't modify
    the default value unless multiple experiments are going to be run.
- `probability::AbstractFloat = 0.3`: a floating point number, between `0` and `1`,
    that represents the probability of doing each single mutation operation.

# Examples
```jldoctest
julia> #Function without extra arguments
julia> function ExampleMutation(Individual::Genotype)
    #···Code···
    return children::Genotype # Type expressed here for user awareness
end

julia> setMutationOperator(ExampleMutation, probability = 0.1)
```

```jldoctest
julia> #Function with extra arguments
julia> function AlsoExampleMutation(Individual::Genotype, arg1, arg2)
    #···Code···
    return children::Genotype # Type expressed here for user awareness
end

julia> setMutationOperator(AlsoExampleMutation, arg1, arg2, probability = 0.1)
```

See also: [`MutationOperator`](@ref)
"""
function setMutationOperator(method::Function, varArgs...; genj::GeneticJuliaStructure = GenJ,
                             probability::AbstractFloat = 0.3)

    if 0 <= probability <= 1
        varArgs = [i for i in varArgs]
        genj._mutationOp = MutationOperator(method, probability, varArgs)
    else
        error("Mutation probability must be within 0 and 1")
    end
end # function

"""
    setNullMutation()
Makes the mutation non-existent for the evolutionary proccess.
This function constructs a structure of type `CrossoverOperator` and adds it to the main structure.

At least, one between crossover and mutation (or both) must exist in the experiment.
"""
function setNullMutation(; genj::GeneticJuliaStructure = GenJ)
    noFunction() = nothing
    genj._mutationOp = MutationOperator(noFunction, 0.0, [])
end # function

max_stars_repo_name: SergioGmezM/GeneticJulia

max_issues_repo_name: SergioGmezM/GeneticJulia

max_forks_repo_name: SergioGmezM/GeneticJulia


text:
using ClassicalOrthogonalPolynomials, BlockArrays, LazyBandedMatrices, FillArrays, ContinuumArrays, StaticArrays, Test
import ClassicalOrthogonalPolynomials: PiecewiseInterlace, SetindexInterlace, plotgrid

@testset "Interlace" begin
    @testset "Piecewise" begin
        @testset "expansion" begin
            T1,T2 = chebyshevt(-1..0), chebyshevt(0..1)
            T = PiecewiseInterlace(T1, T2)
            @test T[-0.1,1:2:10] ≈ T1[-0.1,1:5]
            @test T[0.1,2:2:10] ≈ T2[0.1,1:5]
            @test T[0.0,1:2:10] ≈ T1[0.0,1:5]
            @test T[0.0,2:2:10] ≈ T2[0.0,1:5]

            x = axes(T,1)
            u = T / T \ exp.(x)
            @test u[-0.1] ≈ exp(-0.1)
            @test u[0.1] ≈ exp(0.1)
            @test u[0.] ≈ 2
        end

        @testset "two-interval ODE" begin
            T1,T2 = chebyshevt(-1..0), chebyshevt(0..1)
            U1,U2 = chebyshevu(-1..0), chebyshevu(0..1)
            T = PiecewiseInterlace(T1, T2)
            U = PiecewiseInterlace(U1, U2)

            @test copy(T) == T

            D = U \ (Derivative(axes(T,1))*T)
            C = U \ T

            A = BlockVcat(T[-1,:]',
                            BlockBroadcastArray(vcat,unitblocks(T1[end,:]),-unitblocks(T2[begin,:]))',
                            D-C)
            N = 20
            M = BlockArray(A[Block.(1:N+1), Block.(1:N)])

            u = M \ [exp(-1); zeros(size(M,1)-1)]
            x = axes(T,1)

            F = factorize(T[:,Block.(Base.OneTo(N))])
            @test F \ exp.(x) ≈ (T \ exp.(x))[Block.(1:N)] ≈ u
        end

        @testset "two-interval p-FEM" begin
            P1,P2 = jacobi(1,1,-1..0), jacobi(1,1,0..1)
            W1,W2 = Weighted(P1), Weighted(P2)

            W = PiecewiseInterlace(W1, W2)

            x = axes(W,1)
            D = Derivative(x)
            Δ = -((D*W)'*(D*W))

            @test (W'exp.(x))[1:2:10] ≈ (W1'*exp.(axes(W1,1)))[1:5]
            @test (W'exp.(x))[2:2:10] ≈ (W2'*exp.(axes(W2,1)))[1:5]
            M = W'W

            # Δ \ (W'exp.(x))
        end

        @testset "three-interval ODE" begin
            d = (-1..0),(0..1),(1..2)
            T = PiecewiseInterlace(chebyshevt.(d)...)
            U = PiecewiseInterlace(chebyshevu.(d)...)

            D = U \ (Derivative(axes(T,1))*T)
            C = U \ T

            A = BlockVcat(T[-1,:]',
                            BlockBroadcastArray(vcat,unitblocks(T.args[1][end,:]),-unitblocks(T.args[2][begin,:]),Zeros((unitblocks(axes(T.args[3],2)),)))',
                            BlockBroadcastArray(vcat,Zeros((unitblocks(axes(T.args[1],2)),)),unitblocks(T.args[2][end,:]),-unitblocks(T.args[3][begin,:]))',
                            D-C)
            N = 20
            M = BlockArray(A[Block.(1:N+2), Block.(1:N)])

            u = M \ [exp(-1); zeros(size(M,1)-1)]
            x = axes(T,1)

            F = factorize(T[:,Block.(Base.OneTo(N))])
            @test F \ exp.(x) ≈ (T \ exp.(x))[Block.(1:N)] ≈ u
        end

        @testset "plot" begin
            T1,T2 = chebyshevt(-1..0), chebyshevt(0..1)
            T = PiecewiseInterlace(T1, T2)
            @test plotgrid(T[:,1:5]) == sort([plotgrid(T1[:,1:3]); plotgrid(T2[:,1:3])])
        end
    end

    @testset "SetindexInterlace" begin
        @testset "Chebyshev" begin
            T = ChebyshevT(); U = ChebyshevU();
            V = SetindexInterlace(SVector(0.0,0.0), T, U)
            C² = Ultraspherical(2)
            C³ = Ultraspherical(3)
            M = SetindexInterlace(zero(SMatrix{2,2,Float64}), T, U, C², C³)
            x = axes(T,1)

            @test axes(V,1) == x

            @testset "evaluation" begin
                @test V[0.1,1:2:6] == vcat.(T[0.1,1:3],0)
                @test V[0.1,2:2:6] == vcat.(0,U[0.1,1:3])
                @test M[0.1,1:4:12] == hvcat.(2, T[0.1,1:3], 0, 0, 0)
                @test M[0.1,2:4:12] == hvcat.(2, 0, 0, U[0.1,1:3], 0)
                @test M[0.1,3:4:12] == hvcat.(2, 0, C²[0.1,1:3], 0, 0)
                @test M[0.1,4:4:12] == hvcat.(2, 0, 0, 0, C³[0.1,1:3])
            end

            @testset "expansion" begin
                f = broadcast(x -> SVector(1,x),x)
                c = V[:,Block.(Base.OneTo(5))] \ f
                @test c[1] ≈ 1
                @test V[:,1:5] \ f ≈ [1; 0; 0; 0.5; 0]
                c = V \ f
                u = V * c
                @test u[0.1] ≈ f[0.1]

                F = broadcast(x -> SMatrix{2,2}(1,x,exp(x),cos(x)),x)
                u = M / M \ F
                @test u[0.1] ≈ F[0.1]

                u = V / V \ SVector.(exp.(x), cos.(x))
                @test u[0.1] ≈ [exp(0.1),cos(0.1)]
            end

            @testset "Operators" begin
                W = SetindexInterlace(SVector(0.0,0.0), U, C²)
                R = W \ V
                @test (W * R * (V \ broadcast(x -> SVector(exp(x),cos(x-1)),x)))[0.1] ≈ [exp(0.1),cos(0.1-1)]
                D = W\Derivative(x)*V
                @test (W * D * (V \ broadcast(x -> SVector(exp(x),cos(x-1)),x)))[0.1] ≈ [exp(0.1),-sin(0.1-1)]
            end
        end
        @testset "Fourier" begin
            F = Fourier()
            V = SetindexInterlace{SVector{2,Float64}}(F, F)
            θ = axes(V,1)
            f = broadcast(θ -> SVector(sin.(θ),cos.(θ)),θ)
            u = V / V \ f
            @test u[0.1] ≈ [sin(0.1),cos(0.1)]
        end
        @testset "Fill" begin
            d = 10
            T = ChebyshevT()
            V = SetindexInterlace(zeros(d), Fill(T,d))
            x = axes(V,1)
            V_N = V[:,Block.(Base.OneTo(50))]
            U_N = V_N / V_N \ broadcast(x -> cos.((1:10) .* x), x)
            @test U_N[0.1] ≈ cos.((1:10) .* 0.1)
            U = V / V \ broadcast(x -> cos.((1:10) .* x), x)
            @test U[0.1] ≈ cos.((1:10) .* 0.1)
        end
    end
end