text:
#===============================================================================

     standard_curve.jl

===============================================================================#

import JSON.json
import DataFrames: DataFrame, by
import Memento: debug, warn, error

## if isnull(sample) well not considered
## what if isnull(cq)

#===============================================================================
    functions >>
===============================================================================#

## called by dispatch()
## `/slope` for log (DNA copy#) is on the x-axis and Cq on the y-axis, otherwise `*slope`
function act(
    ::Type{Val{standard_curve}},
    req         ::Vector{Any};
    out_format  ::OutputFormat = pre_json_output,
    json_digits ::Int_T = JSON_DIGITS,
    qty_base    ::Real = 10,
    empty_tre   ::TargetResultEle = EMPTY_TRE,
    empty_gre   ::GroupResultEle  = EMPTY_GRE
)
    debug(logger, "at act(::Type{Val{standard_curve}})")
    #
    ## parse data
    req_df =
        try
            parse_raw_data(Val{standard_curve}, req)
        catch
            throw(ArgumentError(
                "cannot parse raw data for standard curve analysis"))
        end ## try
    #
    ## empty dataset
    if  size(req_df, 2) == 0 ||
        [:target, :cq, :qty] |> mold(symbl -> all(isnan.(req_df[symbl]))) |> any
            output = OrderedDict(:targets => [], :groups => [], :valid => false)
            return output |> out(out_format)
    end
    #
    ## target result set
    use = .!isnan.(req_df[:, :qty]) .& .!isnan.(req_df[:, :cq])
    target_result_df = by(req_df[use, :], :target) do chunk_target
        target_id = chunk_target[1, :target]
        if isnan(target_id)
            empty_tre
        else
            ## better to use GLM.jl here
            x =
                if qty_base == 10
                    log10.(chunk_target[:, :qty])
                else
                    log.(qty_base, chunk_target[:, :qty])
                end
            y = chunk_target[:, :cq]
            b0, b1 = linreg(x, y)
            eff =
                if qty_base == 10
                    exp10(-1.0 / b1) - 1.0
                else
                    exp(-log(e, qty_base) / b1) - 1.0
                end
            yhat = b0 + b1 .* x
            ybar = mean(y)
            ss_res = y .- yhat |> sumsq
            ss_tot = y .- ybar |> sumsq
            r2_ = 1.0 - ss_res / ss_tot
            ## adjusted R2 ?
            TargetResultEle(
                target_id,
                round.([b1, b0, eff, r2_], json_digits)...)
        end ## if isnan
    end ## by chunk_target

    ## group results calculation commented out
    ## as this calculation is duplicated in the front end
    ## TP 2019/03/12
    #
    ## group result set
    # group_result_df = by(req_df, :sample) do chunk_sample
    #     if isnan(chunk_sample[1, :sample])
    #         return hcat(
    #             DataFrame(target=0),
    #             DataFrame(x1=empty_gre)) ## assuming the default column name is `:x1`
    #     else
    #         return by(chunk_sample, :target) do chunk_target
    #             target_id = chunk_target[1, :target]
    #             if isnan(target_id)
    #                 return empty_gre
    #             else
    #                 cq_vec  = chunk_target[.!isnan.(chunk_target[:cq]), :cq]
    #                 qty_vec = chunk_target[.!isnan.(chunk_target[:qty]), :qty]
    #                 return GroupResultEle(
    #                     chunk_target[:well] |> unique |> sort,
    #                     chunk_target[1, :target],
    #                     round.([
    #                         mean(cq_vec),
    #                         std(cq_vec),
    #                         mean(qty_vec),
    #                         std(qty_vec)
    #                     ], json_digits)...
    #                 )
    #             end ## if
    #         end ## do chunk_target
    #     end ## if
    # end ## do chunk

    ## report results
    (out_format == full_output) && return target_result_df

    ## out_format != full_output
    ## target result set
    tre_vec = target_result_df[.!isnan.(target_result_df[:target]), 2]
    target_vec = Vector{Any}()
    for tre in tre_vec
        if isnan(tre.slope) && isnan(tre.offset)
            err_msg = "less 2 valid data points of cq and/or qty available for fitting standard curve"
            warn(logger, err_msg)
            target_result = OrderedDict(
                :target_id => tre.target_id,
                :error     => err_msg)
        else
            target_result = tre
        end
        push!(target_vec, target_result)
    end ## next tre

    ## group results reporting commented out
    ## as the calculation is duplicated in the front end
    ## TP 2019/03/12
    #
    ## group result set
    # gre_vec = group_result_df[:, 3]
    # uniq_well_combins = map(x -> x[:well], gre_vec) |> unique
    # grp_vec = Vector{OrderedDict}()
    # for well_combin in uniq_well_combins
    #     if length(well_combin) > 1
    #         grp_target_vec = Vector{OrderedDict}()
    #         gre_idc = find(gre_vec) do gre
    #             gre_wells = getfield(gre, :well)
    #             return length(gre_wells) == length(well_combin) && all(gre_wells .== well_combin)
    #         end
    #         for gre_idx in gre_idc
    #             gre = gre_vec[gre_idx]
    #             target_id = getfield(gre, :target_id)
    #             if target_id != 0
    #                 qty_mean_m, qty_mean_b = scinot(getfield(gre, :qty_mean), json_digits)
    #                 qty_sd_m, qty_sd_b = scinot(getfield(gre, :qty_sd), json_digits)
    #                 push!(grp_target_vec, OrderedDict(
    #                     :target_id              =>  target_id,
    #                     :cq                     =>  OrderedDict(
    #                         :mean               =>      getfield(gre, :cq_mean),
    #                         :standard_deviation =>      getfield(gre, :cq_sd)),
    #                     :quantity               =>  OrderedDict(
    #                         :mean               =>      OrderedDict(
    #                             :m              =>          qty_mean_m,
    #                             :b              =>          qty_mean_b),
    #                         :standard_deviation =>      OrderedDict(
    #                             :m              =>          qty_sd_m,
    #                             :b              =>          qty_sd_b))))
    #             end ## if target_id
    #         end ## do gre_i
    #     push!(grp_vec, OrderedDict(
    #         :wells   => well_combin,
    #         :targets => grp_target_vec))
    #     end ## if
    # end ## do well_combin
    #
    output = OrderedDict(
        :targets => target_vec,
        :groups  => Vector(),
        :valid   => true)
    return output |> out(out_format)
end ## act(::Type{Val{standard_curve}})

## called by act
sumsq(y) = sum(y .^ 2)

#==============================================================================#

## dependencies of `standard_curve`

## parse req into a dataframe
function parse_raw_data(::Type{Val{standard_curve}}, req ::Vector{Any})
    (length(req) == 0) && return DataFrame()
    #
    well_vec    = Vector{Int_T}()
    channel_vec = Vector{Int_T}()
    target_vec  = Vector{Real}()
    cq_vec      = Vector{Float_T}()
    qty_vec     = Vector{Float_T}()
    sample_vec  = Vector{Real}()
    #
    num_channels =
        map(req) do req_ele
            try
                length(req_ele[WELL_KEY])
            catch
                0
            end ## try
        end #= do req_ele =# |> maximum
    #
    @inline nothing2NaN(x) = thing(x) ? x : NaN_T
    for well_i in eachindex(req)
        req_ele = req[well_i]
        for channel_i in 1:num_channels
            measrmt_dict =
                try
                    req_ele[WELL_KEY][channel_i]
                catch
                    Dict{String,Any}()
                end ## try
            if length(measrmt_dict) == 0
                target = cq = qty = NaN_T
            else
                target   = nothing2NaN(measrmt_dict[TARGET_KEY])
                cq       = nothing2NaN(measrmt_dict[CQ_KEY])
                qty_dict = measrmt_dict[QUANTITY_KEY]
                qty      = nothing2NaN(qty_dict["m"]) * 10.0 ^ nothing2NaN(qty_dict["b"])
            end ## if
            push!(well_vec,     well_i)
            push!(channel_vec,  channel_i)
            push!(target_vec,   target)
            push!(cq_vec,       cq)
            push!(qty_vec,      qty)
            push!(sample_vec,   try
                    nothing2NaN(req_ele[SAMPLE_KEY])
                catch
                    NaN_T
                end) ## try
        end ## for channel_i
    end ## for well_i
    #
    return DataFrame(
        well    = well_vec,
        channel = channel_vec,
        target  = target_vec,
        cq      = cq_vec,
        qty     = qty_vec,
        sample  = sample_vec)
end

## end: dependencies of `standard_curve`

## generate standard_curve request for testing, and dependency functions;
## not to be used in production

## generate unique integers
function generate_uniq_ints(num_ints ::Integer, S, rng ::AbstractRNG = Base.GLOBAL_RNG)
    if num_ints > length(S)
        warn(logger, "num_ints($num_ints) > length(S), i.e. no enough values to choose from")
    end
    uniq_ints = unique(rand(rng, S, num_ints))
    while length(uniq_ints) < num_ints
        uniq_ints = unique(vcat(uniq_ints, rand(rng, S, num_ints - length(uniq_ints))))
    end
    return uniq_ints
end

## insert slice of the same element into array
## by searching along one dimension for positions to insert.
## user specifies element, number of insertion actions, destination array, search dimension.
function insert2ary(
    el2ins,
    num_ins             ::Integer,
    ary                 ::AbstractArray,
    seek2ins_along_dim  ::Integer = 1,
    rng                 ::AbstractRNG = Base.GLOBAL_RNG
)
    dim_len = size(ary)[seek2ins_along_dim]
    idx_range = 1:(dim_len + num_ins)
    ins_idc = generate_uniq_ints(num_ins, idx_range, rng)
    ins_vec = sort(ins_idc) .- (1:num_ins) .+ 1
    ins_mtx =
        hcat(
            vcat(1, ins_vec),
            vcat(ins_vec .- 1, dim_len))

    ary_ndims = ndims(ary)
    ary_ut = Array{Union{eltype(ary),typeof(el2ins)},ary_ndims}(ary)
    select_all_idx_vec = Vector{Any}(fill(Colon(), ary_ndims)) ## make sure eltype is a supertype of Vector{Int}
    ary_size = size(ary)
    ins_size = map(1:ary_ndims) do i
        i == seek2ins_along_dim ? 1 : ary_size[i]
    end
    ins_slice = fill(el2ins, ins_size...)

    ary_wins =
        getindex(
            cat(seek2ins_along_dim,
                map(range(1, num_ins + 1)) do i
                    select_idx_vec = copy(select_all_idx_vec)
                    setindex!(select_idx_vec, ins_mtx[i,1] : ins_mtx[i,2], seek2ins_along_dim)
                    cat(seek2ins_along_dim, getindex(ary_ut, select_idx_vec...), ins_slice)
                end...), ## do i
            setindex!(select_all_idx_vec, idx_range, seek2ins_along_dim)...)
    return ary_wins
end ## insert2ary

## generate standard_curve request
function generate_req_sc(;
    ## random unless individually specified
    ## NA not counted as a unique value for `num_uniq`

    num_wells           ::Integer = 16,
    num_channels        ::Integer = 2,

    target_vec          ::AbstractVector = [],
    num_uniq_targets    ::Integer = 10,
    lb_num_na_targets   ::Integer = 2,

    cq_vec              ::AbstractVector = [],
    cq_bounds           ::Tuple = (0.01, 40.0),
    num_na_cqs          ::Integer = 3,

    qm_vec              ::AbstractVector = [],
    qm_bounds           ::Tuple = (1, 10),
    qb_vec              ::AbstractVector = [],
    qb_bounds           ::Tuple{Int_T,Int_T} = (-20, 20),
    num_na_qtys         ::Integer = 0,

    sample_vec          ::AbstractVector = [],
    num_uniq_samples    ::Integer = 4,
    num_na_samples      ::Integer = 2,

    rng_type            ::DataType = MersenneTwister,
    seed                ::Integer = 1
)
    rng = rng_type(seed)

    debug(logger, "at generate_req_sc: num_wells=$num_wells, num_channels=$num_channels")

    if num_uniq_targets < num_channels
        num_na_channels = num_channels - num_uniq_targets
        channelwide_num_na_targets = num_na_channels * num_wells
        na_channel_idc = generate_uniq_ints(num_na_channels, 1:num_channels, rng)
        addi_num_na_targets = max(0, lb_num_na_targets - channelwide_num_na_targets)
        println("num_uniq_targets < num_channels. targets will all be na for channel(s) ",
            join(na_channel_idc, ","))
    else
        channelwide_num_na_targets = 0
        na_channel_idc = Vector{Int_T}()
        addi_num_na_targets = lb_num_na_targets
    end

    num_measurements = num_wells * num_channels
    println("num_wells and num_channels mandate num_measurements to be $num_measurements (including na)")

    num_targets = length(target_vec)
    if num_targets == 0
        println("randomly generating targets...") ## target values should not be the same across different channels for the same well
        num_nna_targets = num_measurements - channelwide_num_na_targets - addi_num_na_targets
        target_vec = insert2ary(NaN_T, addi_num_na_targets, fill(0, num_nna_targets), 1, rng)
        nna_channel_idc =
            find(1:num_channels) do channel_i
                !(channel_i in na_channel_idc)
            end
        available_targets = 1:num_uniq_targets
        num_uniq_targets_perchannel = Int_T(floor(num_uniq_targets / num_channels))
        for channel_i in nna_channel_idc
            target_idc_thischannel = ((1:num_wells) .- 1) .* num_channels .+ channel_i
            nna_target_idc = target_idc_thischannel[
                map(target_idx -> !isnan(target_vec[target_idx]), target_idc_thischannel)]
            uniq_targets_thischannel =
                generate_uniq_ints(num_uniq_targets_perchannel, available_targets, rng)
            available_targets = setdiff(available_targets, uniq_targets_thischannel)
            nna_target_vec = rand(rng, uniq_targets_thischannel, length(nna_target_idc))
            for nna_idx_i in eachindex(nna_target_idc)
                target_vec[nna_target_idc[nna_idx_i]] = nna_target_vec[nna_idx_i]
            end
        end ## for channel_i

    elseif num_targets != num_measurements
        errmsg = "target_vec not empty but length($num_targets) not same as num_measurements($num_measurements)"
        warn(logger, errmsg)
        throw(ErrorException(errmsg))
    end ## if num_targets

    num_cqs = length(cq_vec)
    if num_cqs == 0
        println("randomly generating cq values...")
        num_nna_cqs = num_measurements - num_na_cqs
        nna_cq_vec = rand(rng, num_nna_cqs) .* -(cq_bounds...) .+ cq_bounds[2] # upperbound - (0,1)seq * scaling_factor
        cq_vec = insert2ary(NaN_T, num_na_cqs, nna_cq_vec, 1, rng)
    elseif num_cqs != num_measurements
        errmsg = "cq_vec not empty but length($num_cqs) not same as num_measurements($num_measurements)"
        warn(logger, errmsg)
        throw(ErrorException(errmsg))
    end

    num_qm = length(qm_vec)
    num_qb = length(qb_vec)
    if num_qm == num_qb == 0
        println("randomly generating quantity values...")
        num_nna_qtys = num_measurements - num_na_qtys
        nna_qm_vec = rand(rng, num_nna_qtys) .* -(qm_bounds...) .+ qm_bounds[2]
        nna_qb_vec = rand(rng, qb_bounds[1]:qb_bounds[2], num_nna_qtys)
        nna_qty_mtx = hcat(nna_qm_vec, nna_qb_vec)
        qty_mtx = insert2ary(NaN_T, num_na_qtys, nna_qty_mtx, 1, rng)
        qm_vec = qty_mtx[:, 1]
        qb_vec = qty_mtx[:, 2] # not convert to integer due to NaN
    elseif num_qm != num_qb
        errmsg = "lengths of qm_vec($num_qm) and qb_vec($num_qb) not equal"
        warn(logger, errmsg)
        throw(ErrorException(errmsg))
    elseif num_qm != num_measurements
        errmsg = "qm_vec and qb_vec with equal non-0 length($num_qm) but not same as num_measurements($num_measurements)"
        warn(logger, errmsg)
        throw(ErrorException(errmsg))
    end

    num_samples = length(sample_vec)
    if num_samples == 0
        println("randomly generating samples...")
        num_nna_samples = num_wells - num_na_samples
        nna_sample_vec = rand(rng, 1:num_uniq_samples, num_nna_samples)
        sample_vec = insert2ary(NaN_T, num_na_samples, nna_sample_vec, 1, rng)
    elseif num_samples != num_wells
        errmsg = "sample_vec not empty but length($num_samples) not same as num_wells($num_wells)"
        warn(logger, errmsg)
        throw(ErrorException(errmsg))
    end

    req = map(1:num_wells) do well_i
        OrderedDict(
            :well =>
                map(1:num_channels) do channel_i
                    measrmt_i = (well_i - 1) * num_channels + channel_i
                    OrderedDict(
                        :target     => target_vec[measrmt_i],
                        :cq         => cq_vec[measrmt_i],
                        :quantity   => OrderedDict(
                            :m          => qm_vec[measrmt_i],
                            :b          => qb_vec[measrmt_i]))
                end, ## next channel_i
            :sample => sample_vec[well_i])
    end ## next well_i

    return (json(req), req)
end ## generate_req_sc

## experimental code >>

## recursive functions to simplify output

function simplify(x ::Number, out_format ::OutputFormat)
println("Number $x format $out_format")
    x
end

function simplify(tup ::Tuple, out_format ::OutputFormat)
println("Tuple $tup format $out_format")
    if out_format == full_output
        tup
    else
        Tuple([
            simplify(tup[i], out_format) for i in eachindex(tup)
        ])
    end
end

function simplify(vec ::Vector, out_format ::OutputFormat)
println("Vector $vec format $out_format")
    if out_format == full_output
        vec
    else
        Vector([
            simplify(vec[i], out_format) for i in eachindex(vec)
        ])
    end
end

function simplify(dict ::Associative, out_format ::OutputFormat)
println("Dict $dict format $out_format")
    if out_format == full_output
        dict
    else
        x = OrderedDict()
        for i in keys(dict)
            x[i] = simplify(dict[i], out_format)
        end
        x
    end
end

function simplify(df ::DataFrame, out_format ::OutputFormat)
println("DataFrame $df format $out_format")
    if out_format == full_output
        df
    else
        x = OrderedDict()
        for i in names(df)
            x[i] = simplify(df[i], out_format)
        end
        x
    end
end

function simplify(r ::StandardCurveResult, out_format ::OutputFormat)
println("Result $r format $out_format")
    if out_format == full_output
        r
    else
        x = OrderedDict()
        if (haskey(x,:well) && length(x[:well])==0)         ||  # empty gre
            (haskey(x,:target_id) && isnan(x[:target_id]))      # empty tre
            return x
        end
        for i in fieldnames(r)
            x[i] = simplify(getfield(r,i), out_format)
        end
        x
    end
end

## << experimental code

## notes
## 2018-04-21. not reproducible:
## Tuple{Float16,Float16}((0.01, 40)) gives (Float16(0.04), Float16(40.0))
## Tuple{Float32,Float16}((0.01, 40)) gives (0.04f0, Float16(40.0))
## Tuple{Float64,Float16}((0.01, 40) gives (0.04, Float16(40.0))
## Tuple{Float64,Float16}((0.001, 40)) gives (0.001, Float16(40.0))
## Tuple{Float16,Float16}((0.001, 40)) gives (Float16(0.004), Float16(40.0))
## then as expected

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr


text:
#=
The Computer Language Benchmarks Game
 https://salsa.debian.org/benchmarksgame-team/benchmarksgame/
 direct transliteration of the swift#3 program by Ralph Ganszky and Daniel Muellenborn:
 https://benchmarksgame-team.pages.debian.net/benchmarksgame/program/mandelbrot-swift-3.html
 modified for Julia 1.0 by Simon Danisch
=#
const zerov8 = ntuple(x-> 0f0, 8)

@inline function step_mandel(Zr,Zi,Tr,Ti,cr,ci)
    Zi = 2f0 .* Zr .* Zi .+ ci
    Zr = Tr .- Ti .+ cr
    Tr = Zr .* Zr
    Ti = Zi .* Zi
    return Zr,Zi,Tr,Ti
end

# Calculate mandelbrot set for one Vec8 into one byte
Base.@propagate_inbounds function mand8(cr, ci)
    Zr = zerov8
    Zi = zerov8
    Tr = zerov8
    Ti = zerov8
    t = zerov8
    i = 0

    while i<50
        for _ in 1:5
            Zr,Zi,Tr,Ti = step_mandel(Zr,Zi,Tr,Ti,cr,ci)
            i += 1
        end
        t = Tr .+ Ti
        all(x-> x > 4f0, t) && (return 0x00)
    end
    byte = 0xff
    t[1] <= 4.0 || (byte &= 0b01111111)
    t[2] <= 4.0 || (byte &= 0b10111111)
    t[3] <= 4.0 || (byte &= 0b11011111)
    t[4] <= 4.0 || (byte &= 0b11101111)
    t[5] <= 4.0 || (byte &= 0b11110111)
    t[6] <= 4.0 || (byte &= 0b11111011)
    t[7] <= 4.0 || (byte &= 0b11111101)
    t[8] <= 4.0 || (byte &= 0b11111110)
    return byte
end

function mandel_inner(rows, ci, y, N, xvals)
    @simd for x in 1:8:N
        @inbounds begin
            cr = ntuple(i-> xvals[x + (i - 1)], 8)
            rows[((y-1)*N÷8+(x-1)÷8) + 1] = mand8(cr, ci)
        end
    end
end

function mandelbrot(n = 200)
    inv_ = 2.0 / n
    N = n
    xvals = zeros(Float32, n)
    yvals = zeros(Float32, n)
    Threads.@threads for i in 0:(N-1)
        @inbounds xvals[i + 1] = i * inv_ - 1.5
        @inbounds yvals[i + 1] = i * inv_ - 1.0
    end
    rows = zeros(UInt8, n*N÷8)
    Threads.@threads for y in 1:N
        @inbounds ci = yvals[y]
        mandel_inner(rows, ci, y, N, xvals)
    end
    write(stdout, "P4\n$n $n\n")
    write(stdout, rows)
end

mandelbrot(parse(Int, ARGS[1]))

ext: jl


text:
function Base.show(io::IO, ::MIME"text/html", model::Model)

    source = """

             Model

            name
            $(model.name)

            filename
            $(model.filename)

    """

    header = "
Type	Equation
\n"

    table_lines = []
    for eq_type in keys(model.equations)
        equations = model.equations[eq_type]
        if eq_type in [:controls_ub, :controls_lb]
            continue
        end
        for i=1:length(equations)
            if i==1
                eqt = eq_type
            else
                eqt = ""
            end
            eq = equations[i]
            eq = sanitize(eq, model)
            eqtex = Dolang.latex(eq)
            eqtex = replace(eqtex, "*"=>" ")
            fmt_line = "
$eqt	\\[$eqtex\\]
"
            push!(table_lines, fmt_line)
        end
    end

    source = string(source, "
", header, table_lines..., "
")
    print(io, source)
end


text:
module Lattices

greet() = print("Hello World!")

export SquareLattice, print_lat, flip!, get_spins#, set_spin!

include("utils.jl")
include("structures.jl")

end # module


text:
using Base.Test
include("../src/Etcd.jl")
include("etcd_mock.jl")

# XXX ideally this macro would create that temp func and use that
# instead of this interface here
@etcd_mock test_machines(et) = Etcd.machines(et)
@etcd_mock test_set(et,k,v,t) = Etcd.set(et,k,v,ttl=t)
@etcd_mock test_update(et,k,v,t) = Etcd.update(et,k,v,ttl=t)
@etcd_mock test_create(et,k,v,t) = Etcd.create(et,k,v,ttl=t)
@etcd_mock test_get(et,k,s,r) = Etcd.get(et,k,sort=s,recursive=r)
@etcd_mock test_create_dir(et,k,t) = Etcd.create_dir(et,k,ttl=t)
@etcd_mock test_update_dir(et,k,t) = Etcd.update_dir(et,k,ttl=t)
@etcd_mock test_add_child(et,k,v,t) = Etcd.add_child(et,k,v,ttl=t)
@etcd_mock test_add_child_dir(et,k,t) = Etcd.add_child_dir(et,k,ttl=t)
@etcd_mock test_delete(et,k) = Etcd.delete(et,k)
@etcd_mock test_set_dir(et,k,t) = Etcd.set_dir(et,k,ttl=t)
@etcd_mock test_delete_dir(et,k,r) = Etcd.delete_dir(et,k,recursive=r)
@etcd_mock test_compare_and_delete(et,k,prvv,prvi) = Etcd.compare_and_delete(et,k,
                                                                             prev_value=prvv,
                                                                             prev_index=prvi)
@etcd_mock test_compare_and_swap(et,k,v,pv,previ,t) = Etcd.compare_and_swap(et,k,v,
                                                                            prev_value=pv,
                                                                            prev_index=previ,
                                                                            ttl=t)

function setup_etcd()
    et = Etcd.EtcdServer()
    println("Etcd server created at: ",et)
    et
end

function test_etcd_machines(et)
    mach = test_machines(et)
    @test mach == "http://127.0.0.1:4001"
end

function test_etcd_set(et)
    key = "/test"
    val = "testvalue"
    set_node = test_set(et,key,val,5)
    @test haskey(set_node,"node")
    @test haskey(set_node["node"],"key")
    @test set_node["node"]["key"] == key
    @test haskey(set_node["node"],"value")
    @test set_node["node"]["value"] == val
    @test set_node["node"]["ttl"] == 5
end

function test_etcd_update(et)
    key = "/test_update"
    val = "update-val"
    new_val = "updatED-val"

    test_set(et,key,val,5)

    update = test_update(et,key,"updatED-val",5)
    @test update["action"] == "update"
    @test update["node"]["key"] == key
    @test update["node"]["ttl"] == 5
    @test update["node"]["value"] == new_val
    @test update["prevNode"]["key"] == key
    @test update["prevNode"]["value"] == val

    # This should fail because the key does not exist.
    update = test_update(et,"/nonexistent-key","anything",5)
    @test haskey(update,"errorCode")
end

function test_etcd_create(et)
    key = "/test_create"
    val = "create-val"

    # this should succeed
    created = test_create(et,key,val,5)
    @test haskey(created,"errorCode") == false
    @test haskey(created,"prevNode") == false
    @test created["action"] == "create"
    @test created["node"]["key"] == key
    @test created["node"]["value"] == val
    @test created["node"]["ttl"] == 5

    # this should fail since the key has been already created
    created = test_create(et,key,val,6)
    @test haskey(created,"errorCode")
end

function test_etcd_set_dir(et)
    key = "/set_dir"
    value = "bar"
    ttl = 5
    # set it
    set_node = test_set(et,key,value,ttl)
    @test haskey(set_node,"errorCode") == false

    # This should succeed it should replace the key
    set_dir = test_set_dir(et,key,5)
    @test set_dir["node"]["key"] == key
    @test haskey(set_dir["node"],"value") == false
    @test set_dir["node"]["ttl"] == 5
    @test set_dir["prevNode"]["key"] == key
    @test set_dir["prevNode"]["value"] == value
    @test set_dir["prevNode"]["ttl"] == 5
end

function test_etcd_update_dir(et)
    key = "/update_dir"

    # make sure creating dir doesn't fail
    cr = test_create_dir(et,key,5)
    @test haskey(cr,"errorCode") == false

    # updating the directory should succeed
    up_dir = test_update_dir(et,key,5)
    @test up_dir["action"] == "update"
    @test up_dir["node"]["key"] == key
    @test haskey(up_dir["node"],"value") == false
    @test up_dir["node"]["ttl"] == 5
    @test up_dir["prevNode"]["key"] == key
    @test haskey(up_dir["prevNode"],"value") == false
    @test up_dir["prevNode"]["ttl"] == 5
    @test up_dir["prevNode"]["dir"] == true

    # updating a non-exitent key should fail
    up_dir = test_update_dir(et,"/nonexistent_key",5)
    @test haskey(cr,"errorCode") == false
end

function test_etcd_create_dir(et)
    key = "/create_dir"
    # make sure creating dir doesn't fail
    cr = test_create_dir(et,key,5)
    @test haskey(cr,"errorCode") == false
    @test cr["action"] == "create"
    @test cr["node"]["key"] == key
    @test haskey(cr["node"],"value") == false
    @test cr["node"]["ttl"] == 5
    @test haskey(cr,"prevNode") == false

    # This should fail, because the key is already there
    cr = test_create_dir(et,key,6)
    @test haskey(cr,"errorCode")
end

function test_etcd_get(et)
    key = "/foo"
    value = "bar"
    ttl = 5
    # set it
    set_node = test_set(et,key,value,ttl)
    # then get it
    get_node = test_get(et,key,false,false)
    @test haskey(get_node,"node")
    @test haskey(get_node["node"],"key")
    @test get_node["node"]["key"] == key
    @test haskey(get_node["node"],"value")
    @test get_node["node"]["value"] == value
end

function test_etcd_get_all(et)
    # create dir
    d_name = "/fooDir"
    ttl = 5
    dir = test_create_dir(et,d_name,ttl)
    k0 = test_set(et,"/fooDir/k0", "v0", ttl)
    k1 = test_set(et,"/fooDir/k1", "v1", ttl)

    # Return kv-pairs in sorted order
    nodes = test_get(et,d_name,true,false)
    @test length(nodes["node"]["nodes"]) == 2
    @test nodes["node"]["nodes"][1]["key"] == "/fooDir/k0"
    @test nodes["node"]["nodes"][1]["value"] == "v0"
    @test nodes["node"]["nodes"][1]["ttl"] == ttl

    @test nodes["node"]["nodes"][2]["key"] == "/fooDir/k1"
    @test nodes["node"]["nodes"][2]["value"] == "v1"
    @test nodes["node"]["nodes"][2]["ttl"] == ttl

    ch_dir = test_create_dir(et,d_name*"/childDir",ttl)
    k2 = test_set(et,d_name*"/childDir/k2", "v2", ttl)

    # recursively get kv-pairs in sorted order
    nodes = test_get(et,d_name,true,true)
    @test length(nodes["node"]["nodes"]) == 3
    @test nodes["node"]["nodes"][1]["nodes"][1]["key"] == "/fooDir/childDir/k2"
    @test nodes["node"]["nodes"][1]["nodes"][1]["value"] == "v2"
    @test nodes["node"]["nodes"][1]["nodes"][1]["ttl"] == ttl

    @test nodes["node"]["nodes"][2]["key"] == "/fooDir/k0"
    @test nodes["node"]["nodes"][2]["value"] == "v0"
    @test nodes["node"]["nodes"][2]["ttl"] == ttl

    @test nodes["node"]["nodes"][3]["key"] == "/fooDir/k1"
    @test nodes["node"]["nodes"][3]["value"] == "v1"
    @test nodes["node"]["nodes"][3]["ttl"] == ttl
end

function test_etcd_add_child(et)
    d_name = "/booDir"
    ch_dir = test_create_dir(et,d_name,5)
    c1 = test_add_child(et,d_name,"v0",5)
    c2 = test_add_child(et,d_name,"v1",5)

    nodes = test_get(et,d_name,true,false)
    @test length(nodes["node"]["nodes"]) == 2
    @test nodes["node"]["nodes"][1]["value"] == "v0"
    @test nodes["node"]["nodes"][2]["value"] == "v1"

    # Creating a child under a nonexistent directory should succeed.
    # The directory should be created.
    c3 = test_add_child(et,"/nonexistentDir","foo",5)
    @test c3["node"]["value"] == "foo"
end

function test_etcd_add_child_dir(et)
    d_name = "/looDir"
    ch_dir = test_create_dir(et,d_name,5)
    c1 = test_add_child_dir(et,d_name,5)
    c2 = test_add_child_dir(et,d_name,6)

    nodes = test_get(et,d_name,true,false)
    @test length(nodes["node"]["nodes"]) == 2
    @test nodes["node"]["nodes"][1]["ttl"] == 5
    @test nodes["node"]["nodes"][2]["ttl"] == 6

    # Creating a child under a nonexistent directory should succeed.
    # The directory should be created.
    c3 = test_add_child_dir(et,"/nonexistentDir",5)
    @test haskey(c3["node"],"key")
end

function test_etcd_delete(et)
    key = "/foo"
    value = "baz"
    ttl = 5
    set_node = test_set(et,key,value,ttl)
    # delete node
    del_node = test_delete(et,key)
    @test del_node["prevNode"]["value"] == value
end

function test_etcd_delete_dir(et)
    test_set_dir(et,"/foo",5)
    # test delete an empty dir
    del_dir = test_delete_dir(et,"/foo",false)

    @test haskey(del_dir["node"],"value") == false
    @test del_dir["prevNode"]["dir"] == true
    @test haskey(del_dir["prevNode"],"value") == false

    # test ability to not delete a non-empty directory
    d_name = "/gooDir"
    ttl = 5
    dir = test_create_dir(et,d_name,ttl)
    foo = test_set(et,d_name*"/goo", "gar", ttl)
    del_dir = test_delete_dir(et,d_name,false)
    @test haskey(del_dir,"errorCode")

    del_dir = test_delete_dir(et,d_name,true)
    @test del_dir["prevNode"]["dir"] == true
    @test haskey(del_dir["prevNode"],"value") == false
    @test haskey(del_dir["node"],"value") == false
end

function test_etcd_compare_and_delete(et)
    key = "/test_c_and_d"
    val = "testvalue"
    set_node = test_set(et,key,val,5)

    cd = test_compare_and_delete(et,key,val,nothing)
    @test cd["prevNode"]["value"] == val
    @test cd["prevNode"]["key"] == key
    @test cd["prevNode"]["ttl"] == 5

    # verify that wrong value fails
    set_node = test_set(et,key,val,6)
    index = set_node["node"]["modifiedIndex"]

    cd = test_compare_and_delete(et,key,"wrongval",nothing)
    @test haskey(cd,"errorCode")

    # should succeed due to correct prevIndex
    cd = test_compare_and_delete(et,key,nothing,index)

    @test cd["prevNode"]["value"] == val
    @test cd["prevNode"]["key"] == key
    @test cd["prevNode"]["ttl"] == 6

    # test giving an incorrect index
    set_node = test_set(et,key,val,7)

    cd = test_compare_and_delete(et,key,nothing,123456)
    @test haskey(cd,"errorCode")
end

function test_etcd_compare_and_swap(et)
    key = "/test_c_and_s"
    val = "testvalue"
    set_node = test_set(et,key,val,5)

    cs = test_compare_and_swap(et,key,"newval",val,nothing,5)
    @test cs["node"]["value"] == "newval"
    @test cs["node"]["key"] == key
    @test cs["node"]["ttl"] == 5
    @test cs["prevNode"]["value"] == val
    @test cs["prevNode"]["key"] == key
    @test cs["prevNode"]["ttl"] == 5

    # verify that wrong prev value fails
    set_node = test_set(et,key,val,6)
    index = set_node["node"]["modifiedIndex"]

    cs = test_compare_and_swap(et,key,"wrongval","wrongprevval",nothing,5)
    @test haskey(cs,"errorCode")

    # swap by index
    cs = test_compare_and_swap(et,key,"newval",nothing,index,5)
    @test cs["node"]["value"] == "newval"
    @test cs["node"]["key"] == key
    @test cs["node"]["ttl"] == 5
    @test cs["prevNode"]["value"] == val
    @test cs["prevNode"]["key"] == key
    @test cs["prevNode"]["ttl"] == 6

    # test giving an incorrect previous index
    set_node = test_set(et,key,val,7)

    cs = test_compare_and_swap(et,key,val,nothing,123456,5)
    @test haskey(cs,"errorCode")
end

function test_etcd()
    et = setup_etcd()
    test_funcs = [test_etcd_machines,
                  test_etcd_set,
                  test_etcd_get,
                  test_etcd_get_all,
                  test_etcd_add_child,
                  test_etcd_add_child_dir,
                  test_etcd_delete,
                  test_etcd_delete_dir,
                  test_etcd_compare_and_delete,
                  test_etcd_compare_and_swap,
                  test_etcd_update,
                  test_etcd_create,
                  test_etcd_set_dir,
                  test_etcd_update_dir,
                  test_etcd_create_dir]
    [f(et) for f in test_funcs]
end

test_etcd()


text:
###
### de-bruijn graph and initial workspace construction
###

# An improvement on the old algorithm:
# Builds a graph of unitigs straight away.
# Uses a sorted list (vector) of canonical kmers instead of a set. Checking kmer
# membership in the list requires binary lookup, which is slower than the O(1) of
# a set's hash, but typical sets of kmers for bigger genomes use waaaaay to much
# memory, so a list it is. But, this algorithm tries to query the list as few times
# as possible, and keeps track of all kmers used up in unitigs.
# Also, unlike the previous algorithm, when it came to making connections, an all
# v all double for loop was used, which was a bit pointless, so we avoid that
# in this version, since the two connection vectors are sorted, we can do a single
# pass of the two vectors instead.
# With decent Kmer counting, this thing should be able to at least put together
# medium size genome like arabidopsis without much trouble. Larger genomes are
# probably fine too, but might need a big machine.
# If we use a large K value that will be allowed by the Skipmer types comming
# in BioSequences v2, the unitig graphs this will produce for medium size genomes
# should be pretty darned good in the first place.

struct Kidx{M<:AbstractMer}
    kmer::M
    idx::UInt64
end

# TODO: Update BioSequences with new neighbour iterators instead. For now, these
# functions will do.
function kmer_fw_neighbours(mer::M) where {M<:AbstractMer}
    d = BioSequences.encoded_data(mer)
    base = d << 2
    return (M(base), M(base + 0x01), M(base + 0x02), M(base + 0x03))
end

function kmer_bw_neighbours(mer::M) where {M<:AbstractMer}
    k = BioSequences.ksize(M) - 1
    d = BioSequences.encoded_data(mer)
    base = d >> 2
    BT = typeof(base)
    return (M(base), M(base + (BT(1) << 2(k))), M(base + (BT(2) << 2(k))), M(base + (BT(3) << 2(k))))
end

function is_end_bw(mer::M, merlist::Vector{M}) where {M<:AbstractMer}
    @debug "Checking if kmer is end BW" mer
    next = Vector{Kidx{M}}()
    get_bw_idxs!(next, mer, merlist)
    @debug "BW neighbours:" next
    length(next) != 1 && return true
    @inbounds p = next[1].kmer
    get_fw_idxs!(next, p, merlist)
    @debug "FW neighbours of only BW neighbour:" p next
    length(next) != 1 && return true
    return false
end

function is_end_fw(mer::M, merlist::Vector{M}) where {M<:AbstractMer}
    @debug "Checking if kmer is end FW" mer
    next = Vector{Kidx{M}}()
    get_fw_idxs!(next, mer, merlist)
    @debug "FW neighbours:" next
    length(next) != 1 && return true
    @inbounds p = next[1].kmer
    get_bw_idxs!(next, p, merlist)
    @debug "BW neighbours of only FW neighbour:" p next
    length(next) != 1 && return true
    return false
end

function get_bw_idxs!(out::Vector{Kidx{M}}, kmer::M, kmerlist::Vector{M}) where {M<:AbstractMer}
    empty!(out)
    for n in kmer_bw_neighbours(kmer)
        cnext = canonical(n)
        cidx = min(searchsortedfirst(kmerlist, cnext), length(kmerlist))
        if @inbounds kmerlist[cidx] == cnext
            push!(out, Kidx{M}(n, cidx))
        end
    end
end

function get_fw_idxs!(out::Vector{Kidx{M}}, kmer::M, kmerlist::Vector{M}) where {M<:AbstractMer}
    empty!(out)
    for n in kmer_fw_neighbours(kmer)
        cnext = canonical(n)
        cidx = min(searchsortedfirst(kmerlist, cnext), length(kmerlist))
        if @inbounds kmerlist[cidx] == cnext
            push!(out, Kidx{M}(n, cidx))
        end
    end
end

const GRAPH_TYPE = Graphs.SequenceDistanceGraph{LongSequence{DNAAlphabet{4}}}

function build_unitigs!(
    sg::Graphs.SequenceDistanceGraph{LongSequence{A}},
    kmerlist::Vector{M}) where {A<:DNAAlphabet,M<:AbstractMer{DNAAlphabet{2}}}

    if !issorted(kmerlist)
        sort!(kmerlist)
    end

    @info string("Constructing unitigs from ", length(kmerlist), " ", BioSequences.ksize(M), "-mers")
    used_kmers = falses(length(kmerlist))

    for start_kmer_idx in eachindex(kmerlist)
        @debug "Considering new kmer" start_kmer_idx

        # Any kmer can only occur in one unitig.
        if used_kmers[start_kmer_idx]
            @debug "Kmer has been used" start_kmer_idx
            continue
        end

        # Check if the kmer is an end/junction of a unitig.
        start_kmer = kmerlist[start_kmer_idx]
        end_bw = is_end_bw(start_kmer, kmerlist)
        end_fw = is_end_fw(start_kmer, kmerlist)

        if !end_bw && !end_fw
            @debug "Kmer is middle of a unitig" start_kmer_idx start_kmer
            continue
        end

        if end_bw && end_fw
            @debug "Kmer is single unitig" start_kmer_idx start_kmer
            # Kmer as unitig
            s = LongSequence{A}(start_kmer)
            used_kmers[start_kmer_idx] = true
        else
            # A unitig starts on this kmer.
            # Make sure the unitig starts on FW.
            current_kmer = start_kmer
            used_kmers[start_kmer_idx] = true
            if end_fw
                current_kmer = reverse_complement(start_kmer)
                end_fw = end_bw
            end
            @debug "Start of unitig" start_kmer current_kmer end_bw end_fw
            # Start unitig
            s = LongSequence{A}(current_kmer)
            fwn = Vector{Kidx{M}}()
            while !end_fw
                # Add end nucleotide, update current kmer.
                get_fw_idxs!(fwn, current_kmer, kmerlist)
                @debug "Extending unitig" fwn
                current_kmer = first(fwn).kmer
                if used_kmers[first(fwn).idx]
                    @debug "New kmer is already used" current_kmer
                    break # Break circular contigs into lines.
                end
                used_kmers[first(fwn).idx] = true
                push!(s, last(current_kmer))
                end_fw = is_end_fw(current_kmer, kmerlist)
            end
        end
        Graphs.add_node!(sg, canonical!(s))
    end

    # Check for perfect circles now.
    next_unused = findnext(!, used_kmers, firstindex(used_kmers))
    if !isnothing(next_unused)
        @info "Perfect circles are present in the input mers"
        @info "Breaking the circles and incorporating them into the graph"
    end
    while !isnothing(next_unused)
        start_kmer = kmerlist[next_unused]
        end_fw = is_end_fw(start_kmer, kmerlist)
        fwn = Vector{Kidx{M}}()
        current_kmer = start_kmer
        used_kmers[next_unused] = true
        s = LongSequence{A}(current_kmer)
        while true
            get_fw_idxs!(fwn, current_kmer, kmerlist)
            current_kmer = first(fwn).kmer
            current_kmer == start_kmer && break
            used_kmers[first(fwn).idx] = true
            push!(s, last(current_kmer))
        end
        Graphs.add_node!(sg, canonical!(s))
        next_unused = findnext(!, used_kmers, next_unused)
    end

    @info string("Constructed ", length(Graphs.nodes(sg)), " unitigs")
    return sg
end

minus_one_k(::Type{Mer{A,K}}) where {A,K} = Mer{A,K-1}
minus_one_k(::Type{BigMer{A,K}}) where {A,K} = BigMer{A,K-1}

function find_unitig_overlaps(sg::GRAPH_TYPE, ::Type{M}) where {M<:AbstractMer}
    @info string("Identifying the ", BioSequences.ksize(M) - 1, "bp (K - 1) overlaps between ", length(Graphs.nodes(sg)), " unitigs")
    # Save the (k-1)mer in (rev on first k-1 / fw on last k-1) or out ( fw on first k-1 / bw on last k-1)
    @debug "Sorting K - 1 overlaps as `in` or `out`"

    in = Vector{Tuple{minus_one_k(M),Graphs.NodeID}}()
    out = Vector{Tuple{minus_one_k(M),Graphs.NodeID}}()
    sizehint!(in, length(Graphs.nodes(sg)))
    sizehint!(out, length(Graphs.nodes(sg)))
    for nid in eachindex(Graphs.nodes(sg))
        nodeseq = Graphs.node(sg, nid).seq
        firstmer = minus_one_k(M)(nodeseq[1:BioSequences.ksize(M) - 1])
        @debug string("Considering node ", nid) nodeseq
        if iscanonical(firstmer)
            @debug "Source overlap is canonical"
            push!(in, (firstmer, nid))
        else
            @debug "Source overlap is not canonical"
            push!(out, (reverse_complement(firstmer), nid))
        end
        lastmer = minus_one_k(M)(nodeseq[end - (BioSequences.ksize(M) - 2):end])
        if iscanonical(lastmer)
            @debug "Sink overlap is canonical"
            push!(out, (lastmer, -nid))
        else
            @debug "Sink overlap is not canonical"
            push!(in, (reverse_complement(lastmer), -nid))
        end
    end
    sort!(in)
    sort!(out)
    return in, out
end

function connect_unitigs_by_overlaps!(sg::GRAPH_TYPE, ::Type{M}) where {M<:AbstractMer}
    in, out = find_unitig_overlaps(sg, M)
    ol = length(out)
    @info string("Linking ", length(Graphs.nodes(sg)), " unitigs by their ", BioSequences.ksize(M) - 1, "bp (K - 1) overlaps")
    # Connect all out -> in for all combinations on each kmer.
    next_out_idx = 1
    for i in in
        while next_out_idx <= ol && first(out[next_out_idx]) < first(i)
            next_out_idx += 1
        end
        oidx = next_out_idx
        while oidx <= ol && first(out[oidx]) == first(i)
            Graphs.add_link!(sg, last(i), last(out[oidx]), -BioSequences.ksize(M) + 1) # No support, although we could add the DBG operation as such.
            oidx += 1
        end
    end
end

#=
function dbg!(ws::WorkSpace, counted_kmers::Vector{<:MerCount}, min_freq::Integer)
    return dbg!(graph(ws), counted_kmers, min_freq)
end
=#

#=
function dbg!(graph::GRAPH_TYPE, kmerlist::Vector{M}) where {M<:AbstractMer}
    str = string("onstructing compressed de-bruijn graph from ", length(kmerlist), ' ', BioSequences.ksize(M), "-mers")
    @info string('C', str)
    build_unitigs_from_sorted_kmers!(graph, kmerlist)
    if Graphs.n_nodes(graph) > 1
        connect_unitigs_by_overlaps!(graph, M)
    end
    @info string("Done c", str)
    return sg
end
=#

function dbg!(sdg::GRAPH_TYPE, kmers)
    if !(eltype(kmers) <: AbstractMer)
        throw(ArgumentError("Didn't provide a collection of kmers for the `kmers` parameter did ya?"))
    end
    @info dbg_message(kmers)
    build_unitigs!(sdg, kmers)
    if Graphs.n_nodes(sdg) > 1
        connect_unitigs_by_overlaps!(sdg, eltype(kmers))
    end
    @info "Done"
    return sdg
end

function dbg!(sg::GRAPH_TYPE, counted_kmers::Vector{<:MerCount}, min_freq::Integer)
    filter!(x -> freq(x) ≥ min_freq, counted_kmers)
    merlist = [mer(x) for x in counted_kmers]
    dbg!(sg, merlist)
end

@inline function dbg_message(x::Vector{M}) where {M<:AbstractMer}
    return string("Constructing a compressed de-bruijn graph from ", length(x), ' ', BioSequences.ksize(M), "-mers")
end

max_issues_repo_name: BioJulia/BioSequenceGraphs

max_forks_repo_name: BioJulia/BioSequenceGraphs


text:
using TopOpt, LinearAlgebra, StatsFuns
using TopOptMakie
using VTKDataIO
using PyCall
using VTKDataTypes
using DashVtk
using Dash
using DashHtmlComponents

#dash = pyimport("dash")
#dash_html = pyimport("dash_html_components")
dash_vtk = pyimport("dash_vtk")

E = 1.0 # Young’s modulus
v = 0.3 # Poisson’s ratio
f = 1.0 # downward force
rmin = 4.0 # filter radius
xmin = 0.0001 # minimum density
problem_size = (160, 40)
x0 = fill(1.0, prod(problem_size)) # initial design
p = 4.0 # penalty
stress_threshold = 1.0 # maximum von Mises stress
compliance_threshold = 800 # maximum compliance

problem = PointLoadCantilever(Val{:Linear}, problem_size, (1.0, 1.0), E, v, f)
#problem = HalfMBB(Val{:Linear}, problem_size, (1.0, 1.0), E, v, f)

solver = FEASolver(Displacement, Direct, problem, xmin = xmin)

cheqfilter = DensityFilter(solver, rmin = rmin)
stress = TopOpt.MicroVonMisesStress(solver)
comp = TopOpt.Compliance(problem, solver)

function obj(x)
    # minimize volume
    return sum(cheqfilter(x)) / length(x)
end
function constr1(x)
    # L-20 norm aggregation to approximate maximum function
    return norm(stress(cheqfilter(x)), 20) - stress_threshold
end
function constr2(x)
    # compliance upper-bound
    return comp(cheqfilter(x)) - compliance_threshold
end

m = Model(obj)
addvar!(m, zeros(length(x0)), ones(length(x0)))
Nonconvex.add_ineq_constraint!(m, constr1)
Nonconvex.add_ineq_constraint!(m, constr2)

options = MMAOptions(
    maxiter=1000, tol = Tolerance(kkt = 1e-4, f = 1e-4),
)
TopOpt.setpenalty!(solver, p)
@time r = Nonconvex.optimize(
    m, MMA87(dualoptimizer = ConjugateGradient()),
    x0, options = options,
);

@show obj(r.minimizer)
@show constr1(r.minimizer)
@show constr2(r.minimizer)
@show maximum(stress(cheqfilter(r.minimizer)))
topology = cheqfilter(r.minimizer);

mesh = VTKUnstructuredData(problem)
topology = round.(topology)
inds = findall(isequal(0), topology)
deleteat!(mesh.cell_connectivity, inds)
deleteat!(mesh.cell_types, inds)

pymesh = PyVTK(mesh)
dash_vtk_utils = pyimport("dash_vtk.utils")
dash_vtk_utils.to_mesh_state
state = dash_vtk_utils.to_mesh_state(pymesh)

content = vtk_view([
        vtk_geometryrepresentation([
        vtk_mesh(state=state)
    ]),
])

app = dash()

app.layout = html_div() do
    html_div(
    style = Dict(
        "width" => "100%",
        "height" => "400px"
    ),
    children=[content],
    )
end

run_server(app, "0.0.0.0", debug = true)


text:
# Model
include_model("car")

# Horizon
T = 51

# Time step
tf = 1.0
h = tf / (T - 1)

# Bounds

# ul <= u <= uu
ul, uu = control_bounds(model, T, -3.0, 3.0)

# Initial and final states
x1 = [0.0; 0.0; 0.0]
xT = [1.0; 1.0; 0.0]
xl, xu = state_bounds(model, T, x1 = x1, xT = xT)

# Objective
obj = quadratic_tracking_objective(
        [(t < T ? Diagonal(ones(model.n))
            : Diagonal(10.0 * ones(model.n))) for t = 1:T],
        [Diagonal(1.0e-1 * ones(model.m)) for t = 1:T-1],
        [xT for t = 1:T], [zeros(model.m) for t = 1:T])

# Constraints
include_constraints("stage")

circles = [(0.85, 0.3, 0.1),
           (0.375, 0.75, 0.1),
           (0.25, 0.2, 0.1),
           (0.75, 0.8, 0.1)]

function circle_obs(x, y, xc, yc, r)
    (x - xc)^2.0 + (y - yc)^2.0 - r^2.0
end

function obstacles!(c, x, u, t)
    c[1] = circle_obs(x[1], x[2], circles[1][1], circles[1][2], circles[1][3])
    c[2] = circle_obs(x[1], x[2], circles[2][1], circles[2][2], circles[2][3])
    c[3] = circle_obs(x[1], x[2], circles[3][1], circles[3][2], circles[3][3])
    c[4] = circle_obs(x[1], x[2], circles[4][1], circles[4][2], circles[4][3])
    nothing
end

n_stage = 4
n_con = n_stage * (T - 1)
t_idx = [t for t = 1:T-1]
con_obstacles = stage_constraints(obstacles!, n_stage, (1:n_stage), t_idx)

# Problem
prob = trajectory_optimization_problem(
           model,
           obj,
           T,
           h = h,
           xl = xl,
           xu = xu,
           ul = ul,
           uu = uu,
           con = con_obstacles)

# Trajectory initialization
x0 = linear_interpolation(x1, xT, T) # linear interpolation on state
u0 = [0.001 * ones(model.m) for t = 1:T-1]

# Pack trajectories into vector
z0 = pack(x0, u0, prob)

# Solve nominal problem
if false
    @time z̄, info = solve(prob, copy(z0))
    @save joinpath(@__DIR__, "sol_to.jld2") z̄
else
    println("Loading solution...")
    @load joinpath(@__DIR__, "sol_to.jld2") z̄
end


text:
#!/usr/bin/env julia --project
using ClimaCoupler
using Unitful

include("utilities/boilerplate.jl")
include("numerics/timestepper_abstractions.jl")
include("balance_law_interface_slabocean.jl")

########
# Set up parameters and initial conditions
########
include("parameters_initialconditions.jl")

FT = Float64

# Collects tracer valus/fluxes for conservation checks
mutable struct LogThatFlux{F}
    A::F
    B::F
end

########
# Set up domain
########

# A: low atmos level (ocean)
domainOcean = SphericalShell(radius = parameters.a - parameters.H, height = parameters.H)

# B: high atmos level (atmosphere)
domainAtmos = SphericalShell(radius = parameters.a, height = parameters.H)

gridOcean = DiscretizedDomain(
    domainOcean;
    elements = (vertical = 4, horizontal = 8),
    polynomial_order = (vertical = 2, horizontal = 2),
    overintegration_order = (vertical = 0, horizontal = 0),
)

gridAtmos = DiscretizedDomain(
    domainAtmos;
    elements = (vertical = 4, horizontal = 8),
    polynomial_order = (vertical = 2, horizontal = 2),
    overintegration_order = (vertical = 0, horizontal = 0),
)

# ########
# # Set up boundary conditions
# ########

oceanBC = (nothing, nothing)

# # Fixed SST:
# T_sfc(𝒫, ϕ) = 𝒫.ΔT * exp(-ϕ^2 / 2 / 𝒫.Δϕ^2) + 𝒫.Tₘᵢₙ
# AtmosSfcBC = BulkFormulaTemperature(
#     drag_coef_temperature = (params, ϕ) -> params.Cₑ,
#     drag_coef_moisture = (params, ϕ) -> params.Cₗ,
#     surface_temperature = T_sfc,
# )
AtmosSfcBC = CoupledPrimarySlabOceanBC(parameters.Cₑ, parameters.Cₗ)
atmosBC = (AtmosSfcBC, DefaultBC())#DefaultBC(), ) inner, outer

########
# Set up model physics
########

FT = Float64

ref_state = DryReferenceState(DecayingTemperatureProfile{FT}(parameters, FT(290), FT(220), FT(8e3)))

physicsOcean = Physics(orientation = SphericalOrientation(), parameters = parameters)

physicsAtmos = Physics(
    orientation = SphericalOrientation(),
    ref_state = ref_state,
    eos = MoistIdealGas(),
    lhs = (NonlinearAdvection{(:ρ, :ρu, :ρe)}(), PressureDivergence()),
    sources = (
        DeepShellCoriolis(),
        FluctuationGravity(),
        ZeroMomentMicrophysics(),
        HeldSuarezForcing(held_suarez_parameters),
        FluxAccumulator(AtmosSfcBC),
    ),
    parameters = parameters,
)

linear_physicsAtmos = Physics(
    orientation = physicsAtmos.orientation,
    ref_state = physicsAtmos.ref_state,
    eos = physicsAtmos.eos,
    lhs = (LinearAdvection{(:ρ, :ρu, :ρe)}(), LinearPressureDivergence()),
    sources = (FluctuationGravity(),),
    parameters = parameters,
)

########
# Set up model
########

modelOcean = SlabOceanModelSetup(
    physics = physicsOcean,
    boundary_conditions = oceanBC,
    initial_conditions = (T_sfc = T_sfc₀,),
    numerics = (flux = nothing, flux_second_order = nothing, direction = EveryDirection()),
)

modelAtmos = DryAtmosModel(
    physics = physicsAtmos,
    boundary_conditions = atmosBC, #atmosBC,
    initial_conditions = (ρ = ρ₀ᶜᵃʳᵗ, ρu = ρu⃗₀ᶜᵃʳᵗ, ρe = ρeᶜᵃʳᵗ, ρq = ρqᶜᵃʳᵗ),
    numerics = (flux = LMARSNumericalFlux(),), #LMARSNumericalFlux;RoeNumericalFlux
)

linear_modelAtmos = DryAtmosModel(
    physics = linear_physicsAtmos,
    boundary_conditions = (DefaultBC(), DefaultBC()), # (Insulating, FreeSlip), (Insulating, FreeSlip)
    initial_conditions = modelAtmos.initial_conditions,
    numerics = (flux = RefanovFlux(),),
)

########
# Set up time steppers (could be done automatically in simulation)
########

nstepsOcean = 1
nstepsAtmos = 1

dx = min_node_distance(gridAtmos.numerical)
cfl = 5 # 13 for 10 days, 7.5 for 200+ days
Δt = cfl * dx / 330.0 * nstepsAtmos

total_steps = 10000
start_time = 0
end_time = Δt * total_steps#30 * 24 * 3600

methodOcean = SSPRK22Heuns
methodAtmos = IMEX()

callbacks = (Info(), CFL(), VTKState(
    iteration = 100, #Int(floor(24*3600/Δt)),
    filepath = "./output/SphereSlabOcean",
), TMARCallback())

########
# Set up simulation
########

simOcean = CplSimulation(
    modelOcean;
    grid = gridOcean,
    timestepper = (method = methodOcean, timestep = Δt / nstepsOcean),
    time = (start = start_time, finish = end_time),
    nsteps = nstepsOcean,
    boundary_z = boundary_mask,
    callbacks = callbacks,
)
simAtmos = CplSimulation(
    (Explicit(modelAtmos), Implicit(linear_modelAtmos));
    grid = gridAtmos,
    timestepper = (method = methodAtmos, timestep = Δt / nstepsAtmos),
    time = (start = start_time, finish = end_time),
    nsteps = nstepsAtmos,
    boundary_z = boundary_mask,
    callbacks = callbacks,
)

## Create a Coupler State object for holding import/export fields.
coupler = CplState()
coupler_register!(
    coupler,
    :SeaSurfaceTemerature,
    deepcopy(simOcean.state.T_sfc[simOcean.boundary]),
    simOcean.grid,
    DateTime(0),
    u"K",
) # value on top of domainA for calculating upward flux into domainB
coupler_register!(
    coupler,
    :BoundaryEnergyFlux,
    deepcopy(simAtmos.state.F_ρe_accum[simAtmos.boundary]),
    simAtmos.grid,
    DateTime(0),
    u"J",
) # downward flux

compOcean = (pre_step = preOcean, component_model = simOcean, post_step = postOcean)
compAtmos = (pre_step = preAtmos, component_model = simAtmos, post_step = postAtmos)
component_list = (domainOcean = compOcean, domainAtmos = compAtmos)

cpl_solver = CplSolver(
    component_list = component_list,
    coupler = coupler,
    coupling_dt = Δt,
    t0 = FT(start_time),
    fluxlog = LogThatFlux(zeros(total_steps), zeros(total_steps)),
)

########
# Run the simulation
########

evolve!(cpl_solver, total_steps)

########
# Check conservation
########

#using Plots
fluxA = cpl_solver.fluxlog.A
fluxB = cpl_solver.fluxlog.B

fluxT = abs.(fluxA) .+ abs.(fluxB)
tme = collect(1:1:total_steps)
rel_error = [((fluxT .- fluxT[1]) / fluxT[1])]
# plot(tme .* cpl_solver.dt, rel_error, ylabel = "rel. error = (fluxT - fluxT[1]) / fluxT[1]", xlabel = "time (s)")
# plot(tme .* cpl_solver.dt, [(fluxA .- fluxA[1])  (fluxB .- fluxB[1]) ],  label = ["Ocean Energy" "Atmos Energy"], xlabel = "time (s)", ylabel = "J / m2")

max_stars_repo_name: CliMA/CouplerMachine

max_issues_repo_name: CliMA/CouplerMachine


text:
function _PDiscRefFE(::Type{T},p::Polytope,orders) where T
    order = first(orders)
    @notimplementedif any( orders .!= order ) "Anisotropic serentopity FEs not allowed"
    _PDiscRefFE(T,p,order)
end

function _PDiscRefFE(::Type{T},p::Polytope,order::Integer) where T
  D = num_cell_dims(p)
  extrusion = tfill(TET_AXIS,Val{D}())
  simplex = ExtrusionPolytope(extrusion)
  reffe = LagrangianRefFE(T,simplex,order)
  metadata = nothing
  face_nodes = [Int[] for face in 1:num_faces(p)]
  face_nodes[end] = collect(1:num_nodes(reffe))
  dofs = get_dof_basis(reffe)
  face_dofs = _generate_face_own_dofs(face_nodes,dofs.node_and_comp_to_dof)

  reffe = GenericRefFE(
    num_dofs(reffe),
    p,
    get_prebasis(reffe),
    get_dof_basis(reffe),
    L2Conformity(),
    metadata,
    face_dofs)
  GenericLagrangianRefFE(reffe,face_nodes)
end


text:
export @esc, isexpr, isline, iscall, rmlines, unblock, block, inexpr, namify, isdef,
  longdef, shortdef, @expand, makeif, prettify, combinedef, splitdef, splitarg

"""
    assoc!(d, k, v)

is the same as `d[k] = v` but returns `d` rather than `v`.
"""
assoc!(d, k, v) = (d[k] = v; d)

"""
    @esc x y

is the same as:
```julia
x = esc(x)
y = esc(y)
```
"""
macro esc(xs...)
  :($([:($x = esc($x)) for x in map(esc, xs)]...);)
end

"""
    @q [expression]

Like the `quote` keyword but doesn't insert line numbers from the construction
site. e.g. compare `@q begin end` with `quote end`. Line numbers of interpolated
expressions are preserverd.

See also: [`rmlines`](@ref)
"""
macro q(ex)
  # use esc here, so we don't rely on https://github.com/JuliaLang/julia/issues/37540 for
  # interpolation to work
  esc(Expr(:quote, striplines(ex)))
end

"""
    isexpr(x, ts...)

Convenient way to test the type of a Julia expression.
Expression heads and types are supported, so for example
you can call

    isexpr(expr, String, :string)

to pick up on all string-like expressions.
"""
isexpr(x::Expr) = true
isexpr(x) = false
isexpr(x::Expr, ts...) = x.head in ts
isexpr(x, ts...) = any(T->isa(T, Type) && isa(x, T), ts)

isline(ex) = isexpr(ex, :line) || isa(ex, LineNumberNode)

iscall(ex, f) = isexpr(ex, :call) && ex.args[1] == f

"""
    rmlines(x)

Remove the line nodes from a block or array of expressions.

Compare `quote end` vs `rmlines(quote end)`

### Examples

To work with nested blocks:

```julia
prewalk(rmlines, ex)
```

See also: [`@q`](@ref)
"""
rmlines(x) = x
function rmlines(x::Expr)
  # Do not strip the first argument to a macrocall, which is
  # required.
  if x.head === :macrocall && length(x.args) >= 2
    Expr(x.head, x.args[1], nothing, filter(x->!isline(x), x.args[3:end])...)
  else
    Expr(x.head, filter(x->!isline(x), x.args)...)
  end
end

striplines(ex) = prewalk(rmlines, ex)

"""
    unblock(expr)

Remove outer `begin` blocks from an expression, if the block is
redundant (i.e. contains only a single expression).
"""
function unblock(ex)
  isexpr(ex, :block) || return ex
  exs = rmlines(ex).args
  length(exs) == 1 || return ex
  return unblock(exs[1])
end

block(ex) = isexpr(ex, :block) ? ex : :($ex;)

"""
An easy way to get the (function/type) name out of
expressions like `foo{T}` or `Bar{T} <: Vector{T}`.
"""
namify(s::Symbol) = s
namify(ex::Expr) = namify(ex.args[1])

walk(x, inner, outer) = outer(x)
walk(x::Expr, inner, outer) = outer(Expr(x.head, map(inner, x.args)...))

"""
    postwalk(f, expr)

Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *after* they have been transformed by the walk.

See also: [`prewalk`](@ref).
"""
postwalk(f, x) = walk(x, x -> postwalk(f, x), f)

"""
    prewalk(f, expr)

Applies `f` to each node in the given expression tree, returning the result.
`f` sees expressions *before* they have been transformed by the walk, and the
walk will be applied to whatever `f` returns.

This makes `prewalk` somewhat prone to infinite loops; you probably want to try
[`postwalk`](@ref) first.
"""
prewalk(f, x)  = walk(f(x), x -> prewalk(f, x), identity)

replace(ex, s, s′) = prewalk(x -> x == s ? s′ : x, ex)

"""
    inexpr(expr, x)

Simple expression match; will return `true` if the expression `x` can be found
inside `expr`.

```julia
inexpr(:(2+2), 2) == true
```
"""
function inexpr(ex, x)
  result = false
  MacroTools.postwalk(ex) do y
    if y == x
      result = true
    end
    return y
  end
  return result
end

isgensym(s::Symbol) = occursin("#", string(s))
isgensym(s) = false

function gensymname(x::Symbol)
  m = Base.match(r"##(.+)#\d+", String(x))
  m === nothing || return m.captures[1]
  m = Base.match(r"#\d+#(.+)", String(x))
  m === nothing || return m.captures[1]
  return "x"
end

"""
    gensym_ids(expr)

Replaces gensyms with unique ids (deterministically).

```julia
julia> x, y = gensym("x"), gensym("y")
(Symbol("##x#363"), Symbol("##y#364"))

julia> MacroTools.gensym_ids(:(\$x+\$y))
:(x_1 + y_2)
```
"""
function gensym_ids(ex)
  counter = 0
  syms = Dict{Symbol, Symbol}()
  prewalk(ex) do x
    isgensym(x) ?
      get!(()->Symbol(gensymname(x), "_", counter+=1), syms, x) :
      x
  end
end

"""
    alias_gensyms(expr)

Replaces gensyms with animal names.
This makes gensym'd code far easier to follow.

```julia
julia> x, y = gensym("x"), gensym("y")
(Symbol("##x#363"), Symbol("##y#364"))

julia> MacroTools.alias_gensyms(:(\$x+\$y))
:(porcupine + gull)
```
"""
function alias_gensyms(ex)
  left = copy(animals)
  syms = Dict{Symbol, Symbol}()
  prewalk(ex) do x
    isgensym(x) ? get!(()->pop!(left), syms, x) : x
  end
end

"""
More convenient macro expansion, e.g.

```julia
@expand @time foo()
```
"""
@static if VERSION <= v"0.7.0-DEV.484"
  macro expand(ex)
    :(alias_gensyms(macroexpand($(current_module()), $(ex,)[1])))
  end
else
  macro expand(ex)
    :(alias_gensyms(macroexpand($(__module__), $(ex,)[1])))
  end
end

"Test for function definition expressions."
isdef(ex) = isshortdef(ex) || longdef1(ex) !== nothing

isshortdef(ex) = (@capture(ex, (fcall_ = body_)) &&
                  (@capture(gatherwheres(fcall)[1],
                            (f_(args__) |
                             f_(args__)::rtype_))))

function longdef1(ex)
  if @capture(ex, (arg_ -> body_))
    Expr(:function, arg isa Symbol ? :($arg,) : arg, body)
  elseif isshortdef(ex)
    @assert @capture(ex, (fcall_ = body_))
    Expr(:function, fcall, body)
  else
    ex
  end
end
longdef(ex) = prewalk(longdef1, ex)

function shortdef1(ex)
  @match ex begin
    function f_(args__) body_ end => @q $f($(args...)) = $(body.args...)
    function f_(args__) where T__ body_ end => @q $f($(args...)) where $(T...) = $(body.args...)
    function f_(args__)::rtype_ body_ end => @q $f($(args...))::$rtype = $(body.args...)
    function f_(args__)::rtype_ where T__ body_ end => @q ($f($(args...))::$rtype) where $(T...) = $(body.args...)
    function (args__,) body_ end => @q ($(args...),) -> $(body.args...)
    ((args__,) -> body_) => ex
    (arg_ -> body_) => @q ($arg,) -> $(body.args...)
    _ => ex
  end
end
shortdef(ex) = prewalk(shortdef1, ex)

"""
```julia
gatherwheres(:(f(x::T, y::U) where T where U)) == (:(f(x::T, y::U)), (:U, :T))
```
"""
function gatherwheres(ex)
  if @capture(ex, (f_ where {params1__}))
    f2, params2 = gatherwheres(f)
    (f2, (params1..., params2...))
  else
    (ex, ())
  end
end

"""    splitdef(fdef)

Match any function definition

```julia
function name{params}(args; kwargs)::rtype where {whereparams}
   body
end
```

and return `Dict(:name=>..., :args=>..., etc.)`. The definition can be rebuilt by
calling `MacroTools.combinedef(dict)`.

See also: [`combinedef`](@ref)
"""
function splitdef(fdef)
  error_msg = "Not a function definition: $(repr(fdef))"
  @assert(@capture(longdef1(fdef),
                   function (fcall_ | fcall_) body_ end),
          "Not a function definition: $fdef")
  fcall_nowhere, whereparams = gatherwheres(fcall)
  func = args = kwargs = rtype = nothing
  if @capture(fcall_nowhere, ((func_(args__; kwargs__)) |
                              (func_(args__; kwargs__)::rtype_) |
                              (func_(args__)) |
                              (func_(args__)::rtype_)))
  elseif isexpr(fcall_nowhere, :tuple)
    if length(fcall_nowhere.args) > 1 && isexpr(fcall_nowhere.args[1], :parameters)
      args = fcall_nowhere.args[2:end]
      kwargs = fcall_nowhere.args[1].args
    else
      args = fcall_nowhere.args
    end
  elseif isexpr(fcall_nowhere, :(::))
    args = Any[fcall_nowhere]
  else
    throw(ArgumentError(error_msg))
  end
  if func !== nothing
    @assert(@capture(func, (fname_{params__} | fname_)), error_msg)
    di = Dict(:name=>fname, :args=>args,
              :kwargs=>(kwargs===nothing ? [] : kwargs), :body=>body)
  else
    params = nothing
    di = Dict(:args=>args, :kwargs=>(kwargs===nothing ? [] : kwargs), :body=>body)
  end
  if rtype !== nothing; di[:rtype] = rtype end
  if whereparams !== nothing; di[:whereparams] = whereparams end
  if params !== nothing; di[:params] = params end
  di
end

"""
    combinedef(dict::Dict)

`combinedef` is the inverse of [`splitdef`](@ref). It takes a `splitdef`-like Dict
and returns a function definition.

This function approximately does the following (but more sophisticated to avoid
emitting parts that did not actually appear in the original function definition.)

```julia
rtype = get(dict, :rtype, :Any)
all_params = [get(dict, :params, [])..., get(dict, :whereparams, [])...]
:(function \$(dict[:name]){\$(all_params...)}(\$(dict[:args]...);
                                            \$(dict[:kwargs]...))::\$rtype
      \$(dict[:body])
  end)
```

"""
function combinedef(dict::Dict)
  rtype = get(dict, :rtype, nothing)
  params = get(dict, :params, [])
  wparams = get(dict, :whereparams, [])
  body = block(dict[:body])
  if haskey(dict, :name)
    name = dict[:name]
    name_param = isempty(params) ? name : :($name{$(params...)})
    # We need the `if` to handle parametric inner/outer constructors like
    # SomeType{X}(x::X) where X = SomeType{X}(x, x+2)
    if isempty(wparams)
      if rtype==nothing
        @q(function $name_param($(dict[:args]...);
                                $(dict[:kwargs]...))
          $(body.args...)
          end)
      else
        @q(function $name_param($(dict[:args]...);
                                $(dict[:kwargs]...))::$rtype
          $(body.args...)
          end)
      end
    else
      if rtype==nothing
        @q(function $name_param($(dict[:args]...);
                                $(dict[:kwargs]...)) where {$(wparams...)}
          $(body.args...)
          end)
      else
        @q(function $name_param($(dict[:args]...);
                                $(dict[:kwargs]...))::$rtype where {$(wparams...)}
          $(body.args...)
          end)
      end
    end
  else
    if isempty(dict[:kwargs])
      arg = :($(dict[:args]...),)
    else
      arg = Expr(:tuple, Expr(:parameters, dict[:kwargs]...), dict[:args]...)
    end
    if isempty(wparams)
      if rtype==nothing
        @q($arg -> $body)
      else
        @q(($arg::$rtype) -> $body)
      end
    else
      if rtype==nothing
        @q(($arg where {$(wparams...)}) -> $body)
      else
        @q(($arg::$rtype where {$(wparams...)}) -> $body)
      end
    end
  end
end

"""
    combinearg(arg_name, arg_type, is_splat, default)

`combinearg` is the inverse of [`splitarg`](@ref).
"""
function combinearg(arg_name, arg_type, is_splat, default)
    a = arg_name===nothing ? :(::$arg_type) : :($arg_name::$arg_type)
    a2 = is_splat ? Expr(:..., a) : a
    return default === nothing ? a2 : Expr(:kw, a2, default)
end

"""
    splitarg(arg)

Match function arguments (whether from a definition or a function call) such as
`x::Int=2` and return `(arg_name, arg_type, is_splat, default)`. `arg_name` and
`default` are `nothing` when they are absent. For example:

```julia
julia> map(splitarg, (:(f(a=2, x::Int=nothing, y, args...))).args[2:end])
4-element Array{Tuple{Symbol,Symbol,Bool,Any},1}:
 (:a, :Any, false, 2)
 (:x, :Int, false, :nothing)
 (:y, :Any, false, nothing)
 (:args, :Any, true, nothing)
```

See also: [`combinearg`](@ref)
"""
function splitarg(arg_expr)
    splitvar(arg) =
        (@match arg begin
            ::T_ => (nothing, T)
            name_::T_ => (name, T)
            x_ => (x, :Any)
        end)::NTuple{2,Any} # the pattern `x_` matches any expression
    (is_splat = @capture(arg_expr, arg_expr2_...)) || (arg_expr2 = arg_expr)
    if @capture(arg_expr2, arg_ = default_)
        @assert default !== nothing "splitarg cannot handle `nothing` as a default. Use a quoted `nothing` if possible. (MacroTools#35)"
        return (splitvar(arg)..., is_splat, default)
    else
        return (splitvar(arg_expr2)..., is_splat, nothing)
    end
end

function flatten1(ex)
  isexpr(ex, :block) || return ex
  #ex′ = :(;)
  ex′ = Expr(:block)
  for x in ex.args
    isexpr(x, :block) ? append!(ex′.args, x.args) : push!(ex′.args, x)
  end
  # Don't use `unblock` to preserve line nos
  return length(ex′.args) == 1 ? ex′.args[1] : ex′
end

"""
    flatten(ex)

Flatten any redundant blocks into a single block, over the whole expression.
"""
flatten(ex) = postwalk(flatten1, ex)

function makeif(clauses, els = nothing)
  @static if VERSION < v"0.7.0-"
    foldr((c, ex)->:($(c[1]) ? $(c[2]) : $ex), els, clauses)
  else
    foldr((c, ex)->:($(c[1]) ? $(c[2]) : $ex), clauses; init=els)
  end
end

unresolve1(x) = x
unresolve1(f::Function) = methods(f).mt.name

unresolve(ex) = prewalk(unresolve1, ex)

function resyntax(ex)
  prewalk(ex) do x
    @match x begin
      setfield!(x_, :f_, x_.f_ + v_) => :($x.$f += $v)
      setfield!(x_, :f_, v_) => :($x.$f = $v)
      getindex(x_, i__) => :($x[$(i...)])
      tuple(xs__) => :($(xs...),)
      adjoint(x_) => :($x')
      _ => x
    end
  end
end

"""
    prettify(ex)

Makes generated code generaly nicer to look at.
"""
prettify(ex; lines = false, alias = true) =
  ex |> (lines ? identity : striplines) |> flatten |> unresolve |> resyntax |> (alias ? alias_gensyms : identity)


text:
using Graphs, Printf

# Optional dependencies from Requires.jl:
using .GraphMakie, .GLMakie

"""
GraphMakie.graphplot(pn::TPN, x=pn.x₀, args...; nlabels_align=(:center,:bottom), nlabels_distance=25, elabels_distance=5, kwargs...)

Creates an interactive plot of the TPN `pn` in the initial state `x`. Additional `args...` and `kwargs...` are directly passed on to `graphplot(::SimpleDiGraph)`.
For a description of keyword arguments, see the `GraphMakie.jl` documentation.
"""
function GraphMakie.graphplot(pn::TPN, x=pn.x₀, args...; nlabels_align=(:center,:bottom), nlabels_distance=25, elabels_distance=5, kwargs...)
    g = SimpleDiGraph(length(pn.P)+length(pn.T))

    nlabels = ([
        string.(pn.P)
        ; string.(pn.T)])
    node_marker = [fill(:circle,length(pn.P));fill(:rect,length(pn.T))]

    elabels=Dict{Graphs.SimpleGraphs.SimpleEdge,String}()
    estyles=Dict{Graphs.SimpleGraphs.SimpleEdge,Symbol}()
    rows = rowvals(pn.ΔF)
    vals = nonzeros(pn.ΔF)
    for t_id in eachindex(pn.T)
        t_id′ = t_id + length(pn.P)
        for i in nzrange(pn.ΔF, t_id)
            p_id = rows[i]
            w = vals[i]
            e = if w > 0
                Graphs.SimpleGraphs.SimpleEdge(t_id′,p_id)
            elseif w < 0
                Graphs.SimpleGraphs.SimpleEdge(p_id,t_id′)
            else
                error("Found edge with 0 token change for id $(p_id)")
            end

            add_edge!(g,e)
            elabels[e] = abs(w)==1 ? "" : repr(Int(w))
            estyles[e] = :black
        end
    end

    # Add pure read arcs
    R = pn.C + pn.ΔF .* (pn.ΔF .< 0)
    rows = rowvals(R)
    vals = nonzeros(R)
    for t_id in eachindex(pn.T)
        t_id′ = t_id + length(pn.P)
        for i in nzrange(R, t_id)
            p_id = rows[i]
            w = vals[i]
            e = Graphs.SimpleGraphs.SimpleEdge(p_id,t_id′)

            add_edge!(g,e)
            elabels[e] = abs(w)==1 ? "" : repr(Int(w))
            estyles[e] = :gray
        end
    end

    # draw plot
    #edge_attr=(linestyle=getindex.(Ref(estyles),edges(g)),)
    f,ax,p= graphplot(g, args...; nlabels, elabels=getindex.(Ref(elabels),edges(g)), nlabels_align, nlabels_distance, elabels_distance, node_marker, node_size=0.5, edge_color=getindex.(Ref(estyles),edges(g)), node_attr=(markerspace = SceneSpace, strokecolor=:black, strokewidth=1, glowwidth=5, glowcolor=fill(:transparent, length(pn.T)+length(pn.P))), node_color=fill(:gray,length(pn.T)+length(pn.P)), kwargs...)
    ax.autolimitaspect = 1

    # add important attributes to plot
    tokens = Observable(Point2f[])
    p.attributes.tokens=tokens
    p.attributes.p_marking=Observable(copy(x.m))
    p.attributes.t_marking=Observable(copy(x.h))
    p.attributes.max_delay=lift((m,h)->max_delay(pn, TPNState(m,h)), p.attributes.p_marking,p.attributes.t_marking)
    p.attributes.selected=1

    # draw menu
    f[1,2] = sidebar = GridLayout(width=Fixed(200), tellheight=false)
    sidebar[1,1] =  Label(f, "Selected: ")
    l_idx = sidebar[1,2] =  Label(f, lift(
        idx-> (idx>length(pn.P) ? (string(pn.T[idx-length(pn.P)])*"\n(transition)") : (string(pn.P[idx]))*"\n(place)"),
        p.attributes.selected)
    )

    place_menu = GridLayout()
    transition_menu = GridLayout()
    delay_menu = GridLayout()
    save_menu = GridLayout()
    sidebar[2,1:2] = place_menu
    sidebar[3,1:2] = transition_menu
    sidebar[4,1:2] = delay_menu
    sidebar[5,1:2] = save_menu

    # save menu
    save_menu[1,1:2] = Label(f, "Save frame:")
    t_fn = save_menu[2,1] = Textbox(f, placeholder = "filename")
    b_save = save_menu[2,2] = Button(f, label="Save")

    # place menu
    l_nt_t= place_menu[1,1] = Label(f, "#tokens: ")
    l_nt = place_menu[1,2] =  Label(f, lift((idx,m)-> (idx>length(pn.P) ? "N.A." : string(Int(m[idx]))),p.attributes.selected,p.attributes.p_marking))
    b_at = place_menu[2,1] =  Button(f, label="+1")
    b_st = place_menu[2,2] =  Button(f, label="-1")

    # transition menu
    l_en_t= transition_menu[1,1] =  Label(f, "Status: ")
    l_en = transition_menu[1,2] =  Label(f, lift(
        (idx,m,h)-> if idx ≤ length(pn.P)
            "N.A."
        elseif isready(pn, TPNState(m,h), idx-length(pn.P))
            "ready"
        elseif isenabled(pn, TPNState(m,h), idx-length(pn.P))
            "enabled"
        else
            "disabled"
        end,
        p.attributes.selected,
        p.attributes.p_marking,
        p.attributes.t_marking
    ))
    l_h_t = transition_menu[2,1] = Label(f, "Clock: ")
    l_h = transition_menu[2,2] = Label(f, lift((idx,h)->if idx ≤ length(pn.P)
            "N.A."
        elseif isnothing(h[idx-length(pn.P)])
            "♮"
        else
            @sprintf "%.2f" h[idx-length(pn.P)]
        end,
        p.attributes.selected,
        p.attributes.t_marking
    ))
    b_f  = transition_menu[3,1] = Button(f,label="Fire!")
    b_fa  = transition_menu[3,2] = Button(f,label="Fire all necessary!")

    # delay menu
    delay = Observable(0.0)
    l_dm_t = delay_menu[1,1] = Label(f, "Max. delay:")
    l_dm   = delay_menu[1,2] = Label(f, lift(dm->(@sprintf "%.2f" dm),p.attributes.max_delay))
    s_d    = delay_menu[2,1:2] = Slider(f, snap=true, startvalue=1.0, range=LinRange(0.0,1.0,11))
    t_d    = delay_menu[3,1] = Textbox(f, displayed_string=lift(d->(@sprintf "%.2f" d),delay), validator=(s->(res=tryparse(Float64,s);!isnothing(res) && 0≤res≤p.attributes.max_delay[])), reset_on_defocus=false)
    b_d    = delay_menu[3,2] = Button(f, label="Elapse!")
    t_d.stored_string[] = @sprintf "%.2f" delay[]

    on(t_d.stored_string) do v
        f = parse(Float64,v)
        if p.attributes.max_delay[] ≈ 0.0
            set_close_to!(s_d, 0.0)
        else
            set_close_to!(s_d, f/p.attributes.max_delay[])
        end
    end

    on(s_d.value) do v
        delay[] = v*p.attributes.max_delay[]
    end

    # draw tokens
    scatter!(ax, p.attributes.tokens; color=:black, markersize=0.1, markerspace=SceneSpace, zindex=1)

    # change tokens with nodes, markings
    onany(p.attributes.node_pos, p.attributes.p_marking) do pos,mm
        _tokens = Point2f[]
        for (i,t) in enumerate(mm)
            if t == 0
            elseif t == 1
                push!(_tokens, pos[i])
            else
                r = 0.125
                ϕ = LinRange(0,2π,t+1)[1:end-1]
                append!(_tokens, Ref(pos[i]) .+ r.*Point2f.(sin.(ϕ), cos.(ϕ)))
            end
        end
        empty!(p.attributes.tokens[])
        append!(p.attributes.tokens[],_tokens)
        p.attributes.tokens[] = p.attributes.tokens[]
    end

    # change display based on marking
    onany(p.attributes.p_marking, p.attributes.t_marking) do m,t
        x = TPNState(p.attributes.p_marking[],p.attributes.t_marking[])
        # color enabled transitions
        colors = fill(:gray, length(pn.P)+length(pn.T))
        for t_id in eachindex(pn.T)
            if isready(pn, x, t_id)
                colors[length(pn.P)+t_id] = :red
            elseif isenabled(pn, x, t_id)
                colors[length(pn.P)+t_id] = :blue
            end
        end
        p.attributes.node_color[] .= colors
        p.attributes.node_color[] = p.attributes.node_color[]
    end

    # change display & menu when selected node changes
    on(p.attributes.selected) do idx
        p.attributes.node_attr.glowcolor[] .= :transparent
        p.attributes.node_attr.glowcolor[][idx]=:red
        p.attributes.node_attr.glowcolor[] = p.attributes.node_attr.glowcolor[]
    end

    # fire transition on button click
    on(b_f.clicks) do n
        idx = p.attributes.selected[]-length(pn.P)
        x = TPNState(p.attributes.p_marking[],p.attributes.t_marking[])
        if idx <= 0
            println("No transition selected!")
        elseif !isready(pn, x, idx)
            println("Transition not ready!")
        else
            println("Fired $(pn.T[idx])")
            fire!(pn, x, idx)
            p.attributes.p_marking[] .= x.m
            p.attributes.t_marking[] .= x.h
            p.attributes.p_marking[] = p.attributes.p_marking[]
            p.attributes.t_marking[] = p.attributes.t_marking[]
        end
    end

    # fire all transitions on button click
    on(b_fa.clicks) do n
        x = TPNState(p.attributes.p_marking[],p.attributes.t_marking[])
        transitions = fire_necessary!(pn, x)

        p.attributes.p_marking[] .= x.m
        p.attributes.t_marking[] .= x.h
        p.attributes.p_marking[] = p.attributes.p_marking[]
        p.attributes.t_marking[] = p.attributes.t_marking[]

        if isempty(transitions)
            println("No ready transitions to fire.")
        else
            println("Fired transitions "*join(transitions, ",", " and ")*".")
        end
    end

    # add or remove tokens on button click
    on(b_at.clicks) do n
        idx = p.attributes.selected[]
        if idx > length(pn.P)
            println("No place selected!")
        else
            p.attributes.p_marking[][idx] += 1
            x = TPNState(p.attributes.p_marking[],p.attributes.t_marking[])
            for t_id ∈ eachindex(pn.T)
                update_clock!(pn,x,t_id)
            end
            p.attributes.t_marking[] .= x.h
            p.attributes.p_marking[] = p.attributes.p_marking[]
            p.attributes.t_marking[] = p.attributes.t_marking[]
        end
    end
    on(b_st.clicks) do n
        idx = p.attributes.selected[]
        if idx > length(pn.P)
            println("No place selected!")
        elseif p.attributes.p_marking[][idx] <= 0
            println("Number of tokens must be ≥ 0")
        else
            p.attributes.p_marking[][idx] -= 1
            x = TPNState(p.attributes.p_marking[],p.attributes.t_marking[])
            for t_id ∈ eachindex(pn.T)
                update_clock!(pn,x,t_id)
            end
            p.attributes.t_marking[] .= x.h
            p.attributes.p_marking[] = p.attributes.p_marking[]
            p.attributes.t_marking[] = p.attributes.t_marking[]
        end
    end

    # delay time on button press
    on(b_d.clicks) do n
        if delay[] > p.attributes.max_delay[]
            println("Delay $(delay[]) exceeds max. delay of $(p.attributes.max_delay).")
        else
            x = TPNState(p.attributes.p_marking[],p.attributes.t_marking[])
            elapse!(pn, x, delay[])
            println("Elapsed τ=$(delay[]).")
            p.attributes.p_marking[] .= x.m
            p.attributes.t_marking[] .= x.h
            p.attributes.p_marking[] = p.attributes.p_marking[]
            p.attributes.t_marking[] = p.attributes.t_marking[]
        end
    end

    # save on button press
    on(b_save.clicks) do n
        save(t_fn.stored_string[], ax.scene)
    end

    # add click interaction to select node
    function select_node(idx, event, axis)
        p.attributes.selected=idx
    end
    register_interaction!(ax, :nodeclick, NodeClickHandler(select_node))

    # add drag interaction
    deregister_interaction!(ax, :rectanglezoom)
    register_interaction!(ax, :ndrag, NodeDrag(p))

    hidedecorations!(ax)

    # force execution of some updates
    p.attributes.selected[] = p.attributes.selected[]
    p.attributes.p_marking[] = p.attributes.p_marking[]
    p.attributes.t_marking[] = p.attributes.t_marking[]

    return f,ax,p
end


text:
using GridArrays
using Test, LinearAlgebra, DomainSets, Plots
using GridArrays: MaskedGrid, IndexSubGrid, randompoint, element, elements, cartesianproduct, iscomposite

function delimit(s::AbstractString)
    println()
    println("## ",s)
end
interval_grids = (EquispacedGrid, PeriodicEquispacedGrid, MidpointEquispacedGrid, ChebyshevNodes, ChebyshevExtremae, FourierGrid,
                    ChebyshevUNodes, LegendreNodes)
types = (Float64,BigFloat)

function test_grids(T)
    ## Equispaced grids
    len = 21
    a = -T(1.2)
    b = T(3.5)
    g1 = EquispacedGrid(len, Interval(a, b))
    g2 = PeriodicEquispacedGrid(len, Interval(a-1, b+1))
    g3 = g1 × g2
    g4 = g1 × g3

    test_generic_grid(g1)
    test_generic_grid(g2)
    test_generic_grid(g3)
    test_generic_grid(g4)

    # Test a subgrid
    g5 = g1[10:20]
    @test g5[1] == g1[10]
    @test g5[11] == g1[20]
    @test length(g5) == 20-10+1
    test_generic_grid(g5)
    g6 = g1[10:2:20]
    @test g6[2] == g1[12]
    @test length(g6) == 6

    g = EquispacedGrid(len, a, b)
    idx = 5
    @test g[idx] ≈ a + (idx-1) * (b-a)/(len-1)
    @test g[len] ≈ b
    @test_throws BoundsError g[len+1] == b

    # Test iterations
    (l,s) = grid_iterator1(g)
    @assert s ≈ len * (a+b)/2

    ## Periodic equispaced grids
    len = 120
    a = -T(1.2)
    b = T(3.5)
    g = PeriodicEquispacedGrid(len, a, b)

    idx = 5
    @test g[idx] ≈ a + (idx-1) * (b-a)/len
    @test g[len] ≈ b - step(g)
    @test_throws BoundsError g[len+1] == b

    (l,s) = grid_iterator1(g)
    @test s ≈ (len-1)*(a+b)/2 + a

    (l,s) = grid_iterator2(g)
    @test s ≈ (len-1)*(a+b)/2 + a

    ## Tensor product grids
    len = 11
    g1 = PeriodicEquispacedGrid(len, -one(T), one(T))
    g = g1^2
    @test isperiodic(g1)
    @test isperiodic(g)
    @test size(g) == (length(g1),length(g1))
    @test cartesianproduct(g1) ==g1

    g2 = EquispacedGrid(len, -one(T), one(T))
    @test !isperiodic(g2)
    g = g1 × g2
    @test !isperiodic(g)
    @test length(g) == length(g1) * length(g2)
    @test size(g) == (length(g1),length(g2))
    @test size(g,1) == length(g1)

    @test element(g, 1) == g1
    @test element(g, 2) == g2
    @test element(g,1:2) == g
    @test !iscomposite(g1)&& !iscomposite(g1)
    @test iscomposite(g)
    @test support(g) ≈ support(g1)×support(g2)

    idx1 = 5
    idx2 = 9
    x1 = g1[idx1]
    x2 = g2[idx2]
    x = g[idx1,idx2]
    @test x[1] ≈ x1
    @test x[2] ≈ x2

    (l,s) = grid_iterator1(g)
    @test s ≈ -len

    (l,s) = grid_iterator2(g)
    @test s ≈ -len

    # Test a tensor of a tensor
    g3 = g × g2
    idx1 = 5
    idx2 = 7
    idx3 = 4
    x = g3[idx1,idx2,idx3]
    x1 = g1[idx1]
    x2 = g2[idx2]
    x3 = g2[idx3]
    @test x[1] ≈ x1
    @test x[2] ≈ x2
    @test x[3] ≈ x3

    # Test a mapped grid
    m = interval_map(T(0), T(1), T(2), T(3))
    # Make a MappedGrid by hand because mapped_grid would simplify
    mg1 = MappedGrid(PeriodicEquispacedGrid(30, T(0), T(1)), m)
    test_generic_grid(mg1)
    # Does mapped_grid simplify?
    mg2 = mapped_grid(PeriodicEquispacedGrid(30, T(0), T(1)), m)
    @test typeof(mg2) <: PeriodicEquispacedGrid
    @test infimum(support(mg2)) ≈ T(2)
    @test supremum(support(mg2)) ≈ T(3)

    # Apply a second map and check whether everything simplified
    m2 = interval_map(T(2), T(3), T(4), T(5))
    mg3 = mapped_grid(mg1, m2)
    @test infimum(support(mg3)) ≈ T(4)
    @test supremum(support(mg3)) ≈ T(5)
    @test typeof(supergrid(mg3)) <: PeriodicEquispacedGrid

    # Scattered grid
    pts = rand(T, 10)
    sg = ScatteredGrid(pts)
    test_generic_grid(sg)
end

function test_laguerre(T)
    grid = LaguerreNodes(10,rand(T))
    @test infimum(support(LaguerreNodes(0.,rand(10)))) == 0
    @test supremum(support(LaguerreNodes(0.,rand(10)))) == Inf
    test_generic_grid(grid)
end

function test_hermite(T)
    grid = HermiteNodes(10)
    @test DomainSets.GeometricSpace{T}()== support(HermiteNodes(rand(T,10)))
    test_generic_grid(grid)
end

function test_jacobi(T)
    grid = JacobiNodes(10,rand(T),rand(T))
    @test support(JacobiNodes(T(0),T(0),rand(T,10))) == ChebyshevInterval{T}()
    test_generic_grid(grid)
    @test JacobiNodes(10,zero(T),zero(T)) ≈ LegendreNodes(10)
    @test JacobiNodes(10,T(1//2),T(1//2)) ≈ ChebyshevUNodes(10)
    @test JacobiNodes(10,T(-1//2),T(-1//2)) ≈ ChebyshevTNodes(10)
end

for T in types
    delimit(string(T))
    for GRID in interval_grids
        @testset "$(rpad(string(GRID),80))" begin
            g = instantiate(GRID,10,T)
            test_interval_grid(g)
        end
    end

    for grid in (JacobiNodes(10,rand(T),rand(T)),)
        @testset "$(rpad(string(typeof(grid)),80))" begin
            test_interval_grid(grid)
        end
    end

    @testset "HermiteNodes" begin
        test_hermite(T)
    end

    @testset "Laguerre" begin
        test_laguerre(T)
    end

    @testset "Jacobi" begin
        test_jacobi(T)
    end

    @testset "$(rpad("Specific grid tests",80))" begin
        test_grids(T)
    end
end

function test_generic_subgrid(grid, s)
    @test support(grid) ≈ s
    for x in grid
        @test x ∈ support(grid)
    end
    for x in subindices(grid)
        @test issubindex(x, grid)
    end

    cnt = 0
    for i in eachindex(supergrid(grid))
        if issubindex(i, grid)
            cnt += 1
        end
    end
    @test cnt == length(grid)
end

function test_subgrids()
    delimit("Grid functionality")
    @testset  "SubGrids" begin
    n = 20
    grid1 = EquispacedGrid(n, -1.0, 1.0)
    subgrid1 = MaskedGrid(grid1, -0.5..0.7)
    test_generic_subgrid(subgrid1, -0.5..0.7)

    subgrid2 = IndexSubGrid(grid1, 4:12)
    subgrid3 = subgrid(grid1, -0.5..0.7)
    test_generic_subgrid(subgrid3, -0.5..0.7)

    @test subgrid1 == subgrid3
    @test mask(subgrid1) == mask(subgrid3)
    @test subindices(subgrid1) == subindices(subgrid3)

    @test support(subgrid1) ∈ support(grid1)
    @test support(subgrid2) ∈ support(grid1)

    G1 = EquispacedGrid(n, -1.0, 1.0)
    G2 = EquispacedGrid(n, -1.0, 1.0)
    ProductG = G1 × G2

    C = UnitDisk()
    refgrid = MaskedGrid(ProductG, C)
    circle_grid = subgrid(ProductG, C)
    @test circle_grid isa MaskedGrid
    @test refgrid ≈ circle_grid
    test_generic_subgrid(circle_grid, C)
    @testset "Generic MaskedGrid" begin

        @test (length(circle_grid)/length(ProductG)-pi*0.25) < 0.01

        G1s = IndexSubGrid(G1,2:4)
        G2s = IndexSubGrid(G2,3:5)
        ProductGs = G1s × G2s
        @test G1s[1] == G1[2]
        @test G2s[1] == G2[3]
        @test ProductGs[1,1] == [G1[2],G2[3]]
    end

    # subgrid of a subgrid
    g1 = EquispacedGrid(10,-1,1)^2
    @test mask(subgrid(subgrid(g1,(0.0..1.0)^2),UnitSimplex{2}()))==[i+j<6 for i in 1:5 , j in 1:5]

    g1 = EquispacedGrid(10,0,1)^2
    for x in g1[BitArray([i+j<12 for i in 1:10 , j in 1:10])]
        @test x ∈ UnitSimplex{2}()
    end

    C = UnitInterval()^2
    productgrid = subgrid(ProductG, C)
    refgrid = MaskedGrid(ProductG, C)

    @test supergrid(productgrid) == ProductG
    @test productgrid isa TensorSubGrid
    test_generic_subgrid(productgrid, C)
    refgrid = MaskedGrid(ProductG, C)
    @test reshape(refgrid,10,10) == productgrid
    @test subindices(refgrid) == subindices(productgrid)

    # Generic tests for the subgrids
    @testset "result" for (grid,subgrid) in ( (grid1,subgrid1), (grid1,subgrid2), (ProductG, circle_grid))
        # print("Subgrid is ")
        # println(typeof(subgrid))
        # Count the number of elements in the subgrid
        cnt = 0
        for i in 1:length(grid)
            if issubindex(i, subgrid)
                cnt += 1
            end
        end
        @test cnt == length(subgrid)
    end

    g = subgrid(ScatteredGrid(rand(10)), Interval(0,.5))
    @test support(g) ≈ Interval(0,.5)
    for x in g
        @test x ∈ support(g)
    end

    for x in boundary(EquispacedGrid(100,-1,1)^2,UnitDisk())
        @test norm(x)≈1
    end

    for x in  boundary(subgrid(EquispacedGrid(100,-1,1)^2,UnitDisk()),UnitDisk())
        @test norm(x)≈1
    end
end
end

function test_randomgrids()
    println("Random grids:")
    @testset begin
        d = UnitDisk()
        g = randomgrid(d, 10)
        @test length(g) == 10
        @test length(eltype(g)) == dimension(d)
        @test reduce(&, [x ∈ d for x in g])

        g2 = randomgrid(UnitDisk{BigFloat}(), 10)
        @test eltype(g2[1]) == BigFloat

        g3 = randomgrid(0.0..1.0, 10)
        @test length(g3) == 10
        # 1D is a special case where we don't use vectors of length 1
        @test eltype(g3) == Float64

        box1 = UnitInterval()^2
        p1 = randompoint(box1)
        @test length(p1) == dimension(box1)
        @test p1 ∈ box1
        box2 = 0.0..1.0
        p2 = randompoint(box2)
        @test typeof(p2) == Float64
        @test p2 ∈ box2

        g3 = randompoint(UnionDomain(0.0..1.5,1.5..2.0))
        @test typeof(g3) == Float64
        @test g3 ∈ 0.0..2.0

        g3 = randompoint(IntersectionDomain(0.0..1.5,1.0..2.0))
        @test typeof(g3) == Float64
        @test g3 ∈ 1.0..1.5

        g3 = randompoint(DifferenceDomain(0.0..1.5,1.0..2.0))
        @test typeof(g3) == Float64
        @test g3 ∈ 0.0..1.0

    end
end
using StaticArrays

include("test_modcartesianindices.jl")

include("test_boundingbox.jl")
include("test_broadcast.jl")
include("test_gauss.jl")
test_subgrids()
test_randomgrids()

@testset "Plots" begin
    plot(FourierGrid(4))
    plot(FourierGrid(4)× FourierGrid(4) )
    plot(FourierGrid(4)× FourierGrid(4) × FourierGrid(4) )
    plot(FourierGrid(4)× FourierGrid(4) ,rand(4,4))
    plot(FourierGrid(4) ,rand(4))
end


text:
abstract type Joint{T,Nc} end

Base.show(io::IO, joint::Joint) = summary(io, joint)

getT(joint::Joint{T}) where T = T
getNc(joint::Joint{T,Nc}) where {T,Nc} = Nc

@inline setForce!(joint::Joint) = return
@inline setForce!(joint::Joint, body1::AbstractBody, body2::AbstractBody, Fτ::SVector{0,T}, No) where T = setForce!(joint)
@inline setForce!(joint::Joint, body1::AbstractBody, body2::AbstractBody, Fτ::Nothing, No) = setForce!(joint)

@inline minimalCoordinates(joint::Joint{T,Nc}) where {T,Nc} = @SVector zeros(T, 3 - Nc)
@inline g(joint::Joint{T,Nc}) where {T,Nc} = @SVector zeros(T, Nc)

@inline ∂g∂posa(joint::Joint{T,Nc}) where {T,Nc} = @SMatrix zeros(T, Nc, 6)
@inline ∂g∂posb(joint::Joint{T,Nc}) where {T,Nc} = @SMatrix zeros(T, Nc, 6)
@inline ∂g∂vela(joint::Joint{T,Nc}) where {T,Nc} = @SMatrix zeros(T, Nc, 6)
@inline ∂g∂velb(joint::Joint{T,Nc}) where {T,Nc} = @SMatrix zeros(T, Nc, 6)


text:
using SplitDataset
using PaddedBlocks
using Test

@testset "reading" begin
    a = reshape(1:81, 9, 9)
    mktempdir() do tdir
        ds = H5SplitDataset(tdir, a, (3, 3))
        @test all(a[1:4, 1:4] == ds[1:4, 1:4])
    end
end

@testset "incremental writing" begin
    T = UInt16
    a = T.(reshape(1:9*9*9, 9, 9, 9))
    for bs in [(1, 3, 3), (3,3,1), (9,9,1)]
        mktempdir() do tdir
            ds = H5SplitDataset(T, tdir, size(a), bs)
            bl = Blocks(size(a),bs)
            for sls in PaddedBlocks.slices(bl)
                ds[sls...] .= a[sls...]
            end
            @test all(a .== ds[:, :, :])
        end
    end
end


text:
# PREAMBLE

# PKG_SETUP

# ## Setup

using LinearAlgebra
using Plots
using DirectTrajectoryOptimization

# ## horizon
T = 101

# ## cartpole
num_state = 4
num_action = 1
num_parameter = 0

function cartpole(x, u, w)
    mc = 1.0
    mp = 0.2
    l = 0.5
    g = 9.81

    q = x[1:2]
    qd = x[3:4]

    s = sin(q[2])
    c = cos(q[2])

    H = [mc+mp mp*l*c; mp*l*c mp*l^2]
    Hinv = 1.0 / (H[1, 1] * H[2, 2] - H[1, 2] * H[2, 1]) * [H[2, 2] -H[1, 2]; -H[2, 1] H[1, 1]]

    C = [0 -mp*qd[2]*l*s; 0 0]
    G = [0, mp*g*l*s]
    B = [1, 0]

    qdd = -Hinv * (C*qd + G - B*u[1])

    return [qd; qdd]
end

function rk3_explicit(x, u, w)
    h = 0.05 # timestep

    k1 = h * cartpole(x, u, w)
    k2 = h * cartpole(x + 0.5 * k1, u, w)
    k3 = h * cartpole(x - k1 + 2.0 * k2, u, w)

    return x + (k1 + 4.0 * k2 + k3) / 6.0
end

function rk3_implicit(y, x, u, w)
    return y - rk3_explicit(x, u, w)
end

# ## model
dt = Dynamics(rk3_implicit, num_state, num_state, num_action,
    num_parameter=num_parameter)
dyn = [dt for t = 1:T-1]

# ## initialization
x1 = [0.0; 0.0; 0.0; 0.0]
xT = [0.0; π; 0.0; 0.0]

# ## objective
Q = 1.0e-2
R = 1.0e-1
Qf = 1.0e2

ot = (x, u, w) -> 0.5 * Q * dot(x - xT, x - xT) + 0.5 * R * dot(u, u)
oT = (x, u, w) -> 0.5 * Qf * dot(x - xT, x - xT)
ct = Cost(ot, num_state, num_action,
    num_parameter=num_parameter)
cT = Cost(oT, num_state, 0,
    num_parameter=num_parameter)
obj = [[ct for t = 1:T-1]..., cT]

# ## constraints
u_bnd = 3.0
bnd1 = Bound(num_state, num_action,
    action_lower=[-u_bnd],
    action_upper=[u_bnd])
bndt = Bound(num_state, num_action,
    action_lower=[-u_bnd],
    action_upper=[u_bnd])
bndT = Bound(num_state, 0)
bounds = [bnd1, [bndt for t = 2:T-1]..., bndT]

cons = [
            Constraint((x, u, w) -> x - x1, num_state, num_action),
            [Constraint() for t = 2:T-1]...,
            Constraint((x, u, w) -> x - xT, num_state, 0)
       ]

# ## problem
p = Solver(dyn, obj, cons, bounds,
    options=Options{Float64}())

# ## initialize
u_guess = [0.01 * ones(num_action) for t = 1:T-1]
x_rollout = [x1]
for t = 1:T-1
    push!(x_rollout, rk3_explicit(x_rollout[end], u_guess[t], zeros(num_parameter)))
end

initialize_states!(p, x_rollout)
initialize_controls!(p, u_guess)

# ## solve
@time solve!(p)

# ## solution
x_sol, u_sol = get_trajectory(p)

@show x_sol[1]
@show x_sol[T]

# ## state
plot(hcat(x_sol...)')

# ## control
plot(hcat(u_sol[1:end-1]..., u_sol[end-1])', linetype = :steppost)


text:
#    Copyright (c) 2021 J.A. Duffek
#    Copyright (c) 2000 D.M. Spink
#
#    This program is free software; you can redistribute it and/or modify
#    it under the terms of the GNU General Public License as published by
#    the Free Software Foundation; either version 2 of the License, or
#    (at your option) any later version.
#
#    This program is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
#    GNU General Public License for more details.
#
#    You should have received a copy of the GNU General Public License
#    along with this program; if not, see <http://www.gnu.org/licenses/>.

using Plots;

"""
    demoruled()

Demonstration of ruled surface construction.

# Examples:
```julia
julia> demoruled()
```
"""
function demoruled()
crv1 = nrbtestcrv();
crv2 = nrbtform(nrbcirc(4.0,[4.5;0.0],1.0*pi,0.0),vectrans([0.0;4.0;-4.0]));
srf = nrbruled(crv1,crv2);
nrbplot(srf,[40;20],
        title = "Ruled surface construction from two NURBS curves.",
        legend = false);
end # demoruled


text:
export funCompR1

"""
GenSPGL (this function is used for traditional spgl1)

Use:    f,g1,g2 = funCompositeR(A, r, funForward, funPenalty, params)
"""
function funCompR1(A::TA,
                       x::AbstractArray,
                       r::AbstractArray,
                       funForward::Function, funPenalty::Function,
                       nProdAt::Int64,
                       params::Dict{String,Any}) where
                         {TA<:Union{joAbstractLinearOperator,AbstractArray,Function}}

    nProdAt += one(Int64)
    f,v = funPenalty(r, params)

    if ~(params["proxy"])
        g1 = funForward(A, x, -v, params)
        g2 = [zero(eltype(g1))]
    else
        g1,g2 = funForward(A, x, -v, params)
    end

    return f,g1,g2
end


text:
function function_iteration!(g, x::AbstractArray{T}, iter,
                            tol=sqrt(eps(T)), gcur=similar(x), fcur=similar(x)) where T
    for k = 1:iter
        g(gcur, x)
        fcur .= gcur .- x
        x .= gcur
        if norm(fcur, Inf) < tol
            return x, true
        end
    end
    x, false
end
function function_iteration(g, x::AbstractArray{T}, iter,
                            tol=sqrt(eps(T))) where T
    for k = 1:iter
        gcur = g(x)
        fcur = gcur - x
        x = gcur
        if norm(fcur, Inf) < tol
            return x, true
        end
    end
    x, false
end