text:
# module ReadableNumbers

export stringpretty, showpretty

import Base: parse

#=
export
  # generating and showing prettier numeric strings
      stringpretty, showpretty,
  # span char: UTF8 char used to separate spans of contiguous digits
      betweenNums , betweenInts , betweenFlts ,
      betweenNums!, betweenInts!, betweenFlts!,
  # span size: the number of contiguous digits used to form a span
      numsSpanned , intsSpanned , fltsSpanned ,
      numsSpanned!, intsSpanned!, fltsSpanned!
=#

# module level control of numeric string formatting (span char, span size)
#   span char and span size are each the only value within a const vector
#   span char and span size are assignable for all parts of numeric strings
#     or they may be assigned to constrain
#        (a) integers & integral part of float strings
#        (b) floats of abs() < 1.0 & the fractional part of float strings

const charBetweenNums = '_'
const lengthOfNumSpan =  3

const ints_spanned = [ lengthOfNumSpan ]; intsSpanned() = ints_spanned[1]
const between_ints = [ charBetweenNums ]; betweenInts() = between_ints[1]
const flts_spanned = [ lengthOfNumSpan ]; fltsSpanned() = flts_spanned[1]
const between_flts = [ charBetweenNums ]; betweenFlts() = between_flts[1]

#  make numeric strings easier to read

stringpretty(val::Signed, group::Int, sep::Char=betweenInts()) =
    prettyInteger(val, group, sep)
stringpretty(val::Signed, sep::Char, group::Int=intsSpanned()) =
    stringpretty(val, group, sep)
function stringpretty(val::Signed)
    group, sep = intsSpanned(), betweenInts()
    stringpretty(val, group, sep)
end

stringpretty{T<:Signed}(val::Rational{T}, group::Int, sep::Char=betweenInts()) =
    string(prettyInteger(val.num, group, sep),"//",prettyInteger(val.den, group, sep))
stringpretty{T<:Signed}(val::Rational{T}, sep::Char, group::Int=intsSpanned()) =
    stringpretty(val, group, sep)
function stringpretty{T<:Signed}(val::Rational{T})
    group, sep = intsSpanned(), betweenInts()
    stringpretty(val, group, sep)
end

stringpretty(val::AbstractFloat,
        intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char) =
    prettyFloat(val, intGroup, fracGroup, intSep, fltSep)
stringpretty(val::AbstractFloat,
        intGroup::Int, fracGroup::Int, sep::Char=betweenFlts()) =
    stringpretty(val, intGroup, fracGroup, sep, sep)
stringpretty(val::AbstractFloat,
        group::Int, intSep::Char, fltSep::Char) =
    stringpretty(val, group, group, intSep, fltSep)
stringpretty(val::AbstractFloat,
        group::Int, sep::Char=betweenFlts()) =
    stringpretty(val, group, group, sep, sep)
stringpretty(val::AbstractFloat,
        intSep::Char, fltSep::Char, intGroup::Int, fracGroup::Int) =
    stringpretty(val, intGroup, fracGroup, intSep, fltSep)
stringpretty(val::AbstractFloat,
        intSep::Char, fltSep::Char, group::Int) =
    stringpretty(val, group, group, intSep, fltSep)
stringpretty(val::AbstractFloat,
        sep::Char, intGroup::Int, fracGroup::Int) =
    stringpretty(val, intGroup, fracGroup, sep, sep)
stringpretty(val::AbstractFloat,
        sep::Char, group::Int=fltsSpanned()) =
    stringpretty(val, group, group, sep, sep)
function stringpretty(val::AbstractFloat)
    group, sep = fltsSpanned(), betweenFlts()
    stringpretty(val, group, group, sep, sep)
end

function stringpretty(val::Real,
          intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char)
    if !prettyfiable(val)
       ty = typeof(val)
       throw(ErrorException("type $ty is not supported"))
    end
    prettyFloat(string(val), intGroup, fracGroup, intSep, fltSep)
end
stringpretty(val::Real, intGroup::Int, fracGroup::Int, sep::Char=betweenFlts()) =
    stringpretty(val, intGroup, fracGroup, sep, sep)
stringpretty(val::Real, group::Int, intSep::Char, fltSep::Char) =
    stringpretty(val, group, group, intSep, fltSep)
stringpretty(val::Real, group::Int, sep::Char=betweenFlts()) =
    stringpretty(val, group, group, sep, sep)
stringpretty(val::Real, intSep::Char, fltSep::Char, intGroup::Int, fracGroup::Int) =
    stringpretty(val, intGroup, fracGroup, intSep, fltSep)
stringpretty(val::Real, intSep::Char, fltSep::Char, group::Int) =
    stringpretty(val, group, group, intSep, fltSep)
stringpretty(val::Real, sep::Char, intGroup::Int, fracGroup::Int) =
    stringpretty(val, intGroup, fracGroup, sep, sep)
stringpretty(val::Real, sep::Char, group::Int=fltsSpanned()) =
    stringpretty(val, group, group, sep, sep)
function stringpretty(val::Real)
    group, sep = fltsSpanned(), betweenFlts()
    stringpretty(val, group, group, sep, sep)
end

# show easy-to-read numbers

showpretty(io::IO, val::Signed, group::Int, sep::Char=betweenInts()) =
    print(io, stringpretty(val, group, sep))
showpretty(io::IO, val::Signed, sep::Char, group::Int=intsSpanned()) =
    print(io, stringpretty(val, group, sep))
function showpretty(io::IO, val::Signed)
    group, sep = intsSpanned(), betweenInts()
    print(io, stringpretty(val, group, sep))
end

showpretty{T<:Signed}(io::IO, val::Rational{T}, group::Int, sep::Char=betweenInts()) =
    print(io,prettyInteger(val.num, group, sep),"//",prettyInteger(val.den, group, sep))
showpretty{T<:Signed}(io::IO, val::Rational{T}, sep::Char, group::Int=intsSpanned()) =
    print(io, stringpretty(val, group, sep))
function showpretty{T<:Signed}(io::IO, val::Rational{T})
    group, sep = intsSpanned(), betweenInts()
    print(io, stringpretty(val, group, sep))
end

function showpretty(io::IO, val::AbstractFloat)
    group, sep = fltsSpanned(), betweenFlts()
    print(io, stringpretty(val, group, group, sep, sep))
end
function showpretty(io::IO, val::AbstractFloat, group::Int)
    sep = betweenFlts()
    print(io, stringpretty(val, group, group, sep, sep))
end
function showpretty(io::IO, val::AbstractFloat, sep::Char)
    group = fltsSpanned()
    print(io, stringpretty(val, group, group, sep, sep))
end
showpretty(io::IO, val::AbstractFloat, prettyFormat...) =
    print(io, stringpretty(val, prettyFormat...))

function showpretty(io::IO, val::Real,
          intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char)
    if !prettyfiable(val)
       ty = typeof(val)
       throw(ErrorException("type $ty is not supported"))
    end
    print(io, stringpretty(val, intGroup, fracGroup, intSep, fltSep))
end
function showpretty(io::IO, val::Real)
    group, sep = fltsSpanned(), betweenFlts()
    showpretty(io, val, group, group, sep, sep)
end
function showpretty(io::IO, val::Real, group::Int)
    sep = betweenFlts()
    showpretty(io, val, group, group, sep, sep)
end
function showpretty(io::IO, val::Real, sep::Char)
    group = fltsSpanned()
    showpretty(io, val, group, group, sep, sep)
end
showpretty(io::IO, val::Real, prettyFormat...) =
    print(io, stringpretty(val, prettyFormat...))

# show on STDOUT

showpretty(val::Signed, group::Int, sep::Char=betweenInts()) =
    showpretty(Base.STDOUT, val, group, sep)
showpretty(val::Signed, sep::Char, group::Int=intsSpanned()) =
    showpretty(Base.STDOUT, val, group, sep)
function showpretty(val::Signed)
    group, sep = intsSpanned(), betweenInts()
    showpretty(Base.STDOUT, val, group, sep)
end

showpretty{T<:Signed}(val::Rational{T}, group::Int, sep::Char=betweenInts()) =
    string(prettyInteger(val.num, group, sep),"//",prettyInteger(val.den, group, sep))
showpretty{T<:Signed}(val::Rational{T}, sep::Char, group::Int=intsSpanned()) =
    showpretty(Base.STDOUT, val, group, sep)
function showpretty{T<:Signed}(val::Rational{T})
    group, sep = intsSpanned(), betweenInts()
    showpretty(Base.STDOUT, val, group, sep)
end

function showpretty(val::AbstractFloat, intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char)
    showpretty(Base.STDOUT, val, intGroup, fracGroup, intSep, fltSep)
end
showpretty(val::AbstractFloat, group::Int) =
    showpretty(Base.STDOUT, val, group)
showpretty(val::AbstractFloat, sep::Char)  =
    showpretty(Base.STDOUT, val, sep)
showpretty(val::AbstractFloat, prettyFormat...) =
    showpretty(Base.STDOUT, val, prettyFormat...)

function showpretty{T<:Real}(val::T, intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char)
    if !prettyfiable(val)
       throw(ErrorException("type $T is not supported"))
    end
    showpretty(Base.STDOUT, val, intGroup, fracGroup, intSep, fltSep)
end
showpretty(val::Real, group::Int) =
    showpretty(Base.STDOUT, val, group)
showpretty(val::Real, sep::Char)  =
    showpretty(Base.STDOUT, val, sep)
showpretty(val::Real, prettyFormat...) =
    showpretty(Base.STDOUT, val, prettyFormat...)

# accept integers and floats

prettyInteger{T<:Signed}(val::T, group::Int, span::Char) =
    integerString(string(val), group, span)

prettyFloat{T<:AbstractFloat}(val::T,
  intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char) =
    prettyFloat(string(val), intGroup, fracGroup, intSep, fltSep)

prettyFloat{T<:AbstractFloat}(val::T,
  intGroup::Int, fracGroup::Int, span::Char) =
    prettyFloat(string(val), intGroup, fracGroup, span, span)

prettyFloat{T<:AbstractFloat}(val::T,
  group::Int, intSep::Char, fltSep::Char) =
    prettyFloat(string(val), group, intSep, fltSep)

prettyFloat{T<:AbstractFloat}(val::T,  group::Int, span::Char) =
    prettyFloat(string(val), group, span, span)

# handle integer and float strings

if VERSION > v"0.4.999"
   splitstr(str::String, at::String) = map(String, split(str, at))
else
   splitstr(str::String, at::String) = map(bytestring, split(str, at))
end

prettyInteger(s::String, group::Int, span::Char) =
    integerString(s, group, span)

function prettyFloat(s::String, intGroup::Int, fracGroup::Int, intSep::Char, fltSep::Char)
    sinteger, sfrac =
        if contains(s,".")
           splitstr(s,".")
        else
           s, ""
        end

    istr = integerString(sinteger, intGroup, intSep)
    if sfrac == ""
       istr
    else
       fstr = fractionalString(sfrac, fracGroup, fltSep)
       string(istr, ".", fstr)
    end
end

prettyFloat(s::String, group::Int, span::Char) =
    prettyFloat(s, group, group, span, span)

prettyFloat(s::String, group::Int, intSep::Char, fltSep::Char) =
    prettyFloat(s, group, group, intSep, fltSep)

prettyFloat(s::String, intGroup::Int, fracGroup::Int, span::Char) =
    prettyFloat(s, intGroup, fracGroup, span, span)

# do the work

function nonnegIntegerString(s::String, group::Int, span::Char)
    n = length(s)
    n==0 && return "0"

    sinteger, sexponentOf2 =
        if contains(s,"e")
           strsplit(s,"e")
        else
           s, ""
        end

    n = length(sinteger)

    fullgroups, finalgroup = divrem(n, group)

    sv = convert(Vector{Char},sinteger)
    p = repeat(" ", n+(fullgroups-1)+(finalgroup!=0))
    pretty = convert(Vector{Char},p)

    sourceidx = n
    targetidx = length(pretty)
    for k in fullgroups:-1:1
        pretty[(targetidx-group+1):targetidx] = sv[(sourceidx-group+1):sourceidx]
        sourceidx -= group
        targetidx -= group
        if k > 1
            pretty[targetidx] = span
            targetidx -= 1
        end
    end

    if finalgroup > 0
        if fullgroups > 0
            pretty[targetidx] = span
            targetidx -= 1
        end
        pretty[(targetidx-finalgroup+1):targetidx] = sv[(sourceidx-finalgroup+1):sourceidx]
    end

    prettystring = convert(String, pretty)

    if length(sexponentOf2) != 0
       string(prettystring,"e",sexponentOf2)
    else
       prettystring
    end
end

function integerString(s::String, group::Int, span::Char)
    if s[1] != '-'
       nonnegIntegerString(s, group, span)
    else
       s1 = string(s[2:end])
       pretty = nonnegIntegerString(s1, group, span)
       string("-", pretty)
    end
end

function fractionalString(s::String, group::Int, span::Char)
    sfrac, sexponentOf2 =
        if contains(s,"e")
           map(String, split(s,"e"))
        else
           s, ""
        end

    pretty = reverse(nonnegIntegerString(reverse(sfrac), group, span))

    if length(sexponentOf2) != 0
       string(pretty,"e",sexponentOf2)
    else
       pretty
    end
end

# get and set shared parameters

function intsSpanned!(n::Int)
    n = max(0,n)
    ints_spanned[1]   = n
    nothing
end
intsSpanned(n::Int) = intsSpanned!(n)

function fltsSpanned!(n::Int)
    n = max(0,n)
    fltsSpanned[1] = n
    nothing
end
fltsSpanned(n::Int) = fltsSpanned!(n)

numsSpanned() = (intsSpanned() == fltsSpanned()) ? intsSpanned() : (intsSpanned(), fltsSpanned())
function numsSpanned!(n::Int)
    n = max(0,n)
    intsSpanned!(n)
    fltsSpanned!(n)
    nothing
end
numsSpanned(n::Int) = numsSpanned!(n)

betweenNums() = (betweenInts() == betweenFlts()) ? betweenInts() : (betweenInts(), betweenFlts())
function betweenNums!(ch::Char)
    betweenFlts!(ch)
    betweenInts!(ch)
    nothing
end
betweenNums!(s::String) = betweenNums!(s[1])
betweenNums(ch::Int)    = betweenNums!(ch)
betweenNums(s::String)  = betweenNums!(s)

function betweenInts!(ch::Char)
    between_ints[1] = ch
    nothing
end
betweenInts(ch::Char)  = betweenInts!(ch)
betweenInts(s::String) = betweenInts!(s[1])

function betweenFlts!(ch::Char)
    between_flts[1] = ch
    nothing
end
betweenFlts(ch::Char)  = betweenFlts!(ch)
betweenFlts(s::String) = betweenFlts!(s[1])

# is this a type that can be handled above
function prettyfiable{T<:Real}(val::T)
    try
        convert(BigFloat,val); true
    catch
        false
    end
end

# parse pretty numeric strings
parse{T<:Union{Signed,AbstractFloat}}(::Type{T}, s::String, ch::Char) =
    parse(T, join(split(s,ch),""))
parse{T<:AbstractFloat}(::Type{T}, s::String, ch1::Char, ch2::Char) =
    parse(T, join(split(s,(ch1,ch2)),""))

# end # module


text:
# Thread based scheduling for StanBase?

function intensive_computation()
    started = time_ns()
    while time_ns() - started < 1_000_000_000
    end
end

@show Threads.nthreads()
let t0 = time_ns()
    tasks = [Threads.@spawn(intensive_computation()) for _ in 1:Threads.nthreads()]
    t1 = time_ns()
    print("all ")
    t2 = time_ns()
    println("spawned")
    t3 = time_ns()
    foreach(wait, tasks)
    t4 = time_ns()
    @show (t1 - t0) / 1_000_000
    @show (t2 - t1) / 1_000_000
    @show (t3 - t2) / 1_000_000
    @show (t4 - t3) / 1_000_000
end

ext: jl


text:
using LibCUDA

struct NN
    g
end

function NN(embeds_w::Matrix{T}, embeds_c::Matrix{T}, ntags::Int;
            nlayers::Int=2, winsize_c::Int=2, winsize_w::Int=5, droprate::Float64=0.1) where T

    embeds_w = zerograd(embeds_w)
    w = lookup(Node(embeds_w), Node(name="w"))

    wsize = size(embeds_w, 1)

    embeds_c = zerograd(embeds_c)
    c = lookup(Node(embeds_c), Node(name="c"))
    batchdims_c = Node(name="batchdims_c")
    c = window1d(c, winsize_c, batchdims_c)

    csize = (winsize_c * 2 + 1) * size(embeds_c, 1)
    c = Linear(T, csize, csize)(c)
    c = maximum(c, 2, batchdims_c)

    h = concat(1, w, c)
    batchdims_w = Node(name="batchdims_w")

    hsize = wsize + csize
    wsize = (winsize_w * 2 + 1) * hsize

    for i in 1:nlayers
        h = Node(window1d, h, winsize_w, batchdims_w)
        h = Node(dropout, h, droprate)
        h = Linear(T, wsize, hsize)(h)
        h = relu(h)
    end

    h = Linear(T, hsize, ntags)(h)
    g = BACKEND(Graph(h))

    NN(g)
end

function (nn::NN)(x::Sample, train::Bool)
    Merlin.CONFIG.train = train
    z = nn.g(x.batchdims_c, x.batchdims_w, Var(x.c), Var(x.w))
    # z = nn.g(x.batchdims_c, Var(x.c), Var(x.w))
    if train
        softmax_crossentropy(Var(x.t), z)
    else
        argmax(z.data, 1)
    end
end


text:
"""
Complex populations may consist of multiple demes, where all members of a given
deme belong to the same species, but each deme may use a different species.

**Properties:**

*`demes::Vector{DemeDefinition}`, an ordered list containing specifications
  for each of the demes within this population.
"""
complex(demes::Vector{_deme_.DemeDefinition}) =
  Population(Deme[_deme_.compose!(d) for d in demes])


text:
using JuMP, Gurobi

function printAllocations(a, edges, T, Tf)
    println("------------------ Allocations ----------------------\n")
    for i in 1:size(a,1)
        for j in 1:size(a,2)
            println("Flow ",i, ", tunnel ", j, " allocated : ", a[i,j])
            print("Edges in use: ")
            for e in T[Tf[i][j]]
                print(edges[e])
            end
            println("\n")
        end
    end
end

function MaxMin(env, edges, capacity, flows, demand, k, T, Tf, explain=false)
    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end
    a = zeros(nflows, k)
    b = zeros(nflows) .- 1
    U = collect(1:nflows)
    i = 0
    while length(U) != 0
        a, amin = MaxMinLP(env, edges, capacity, flows, U, b, demand, k, T, Tf)
        Z = []
        for u in 1:length(U)
            if amin >= demand[U[u]]
                push!(Z, U[u])
                b[U[u]] = sum(a[U[u],:])
            end
        end
        U = filter(u -> u ∉ Z, U)
        i += 1
    end
    if explain
        printAllocations(a, edges, T, Tf)
    end
    return a
end

function MaxMinLP(env,
              edges,
              capacity,
              flows,
              U,
              b_fixed,
              demand,
              k,
              T,
              Tf)

    nedges = size(edges,1)
    nflows = size(flows,1)
    ntunnels = size(T,1)

    #CREATE TUNNEL EDGE MATRIX
    L = zeros(ntunnels, nedges)
    for t in 1:ntunnels
        for e in 1:nedges
            if in(e, T[t])
                L[t,e] = 1
            end
        end
    end

    model = Model(solver=GurobiSolver(env, OutputFlag=0))
    @variable(model, a[1:nflows, 1:k] >= 0, basename="a", category=:SemiCont)
    @variable(model, amin >= 0, basename="amin", category=:SemiCont)

    for e in 1:nedges
        @constraint(model, sum(a[U[u],t] * L[Tf[U[u]][t],e] for u in 1:length(U), t in 1:size(Tf[U[u]],1)) <= capacity[e])
    end

    for u in 1:length(U)
        @constraint(model, sum(a[U[u],t] for t in 1:size(Tf[U[u]],1)) <= demand[U[u]])
    end

    for f in 1:nflows
        if b_fixed[f] != -1
            @constraint(model, sum(a[f,t] for t in 1:size(Tf[f],1)) == b_fixed[f])
        end
    end

    for u in 1:length(U)
        @constraint(model, amin <= sum(a[U[u],t] for t in 1:size(Tf[U[u]],1)))
    end

    @objective(model, Max, amin)

    solve(model)
    return (getvalue(a), getvalue(amin))
end


text:
using MLStyle

function capturing_analysis(expr, out, is_literal)

    @match expr begin
        ::QuoteNode =>
             capturing_analysis(expr.value, out, true)

        if is_literal end && Expr(:$, args...) =>

            foreach(x -> capturing_analysis(x, out, false), args)

        if is_literal end && Expr(_, args...) =>

            foreach(x -> capturing_analysis(x, out, true), args)

        if is_literal end && _ => nothing

        # not literal
        ::Symbol => (push!(out, expr); nothing)

        Expr(:quote, args...) =>
            foreach(x -> capturing_analysis(x, out, true), args)

        :(Do($(args...))) =>
            foreach(args) do arg
                @match arg begin
                    Expr(:kw, key :: Symbol, value) =>
                        begin
                            push!(out, key)
                        end
                    _ => nothing
                end
            end
        :($a || $b) =>
        let out1 = Set(Symbol[]),
            out2 = Set(Symbol[])

            capturing_analysis(a, out1, false)
            capturing_analysis(b, out2, false)

            union!(out, intersect(out1, out2))
            nothing
        end
        # type pattern
        Expr(:(::), a, _) => capturing_analysis(a, out, false)
        # dict pattern
        :(Dict($(args...))) =>
            foreach(args) do arg
                @match arg begin
                    :($_ => $v) => capturing_analysis(v, out, false)
                    _ => nothing
                end
            end
        # app pattern
        :($_($(args...))) => foreach(x -> capturing_analysis(x, out, false), args)
        # other expr
        Expr(_, args...) => foreach(x -> capturing_analysis(x, out, false), args)
        # ref pattern
        Expr(:&, _)              ||
        # predicate
        Expr(:function, _...)    ||
        Expr(:if, _...)          ||
        x => nothing
    end
end

function capture(template)
    syms = Set(Symbol[])
    capturing_analysis(template, syms, true)
    out_expr = Expr(:call, Dict, (Expr(:call, =>, QuoteNode(each), each) for each in syms)...)
    let template = Expr(:quote, template)
        quote
            @λ $template -> $out_expr
        end
    end
end

macro capture(template)
    capture(template) |> esc
end
macro capture(template, ex)
    Expr(:call, capture(template), ex) |> esc
end

@info @capture f($x) :(f(1))

destruct_fn = @capture function $(fname :: Symbol)(a, $(args...)) $(body...) end
@info destruct_fn(:(
    function f(a, x, y, z)
        x + y + z
    end
))

max_stars_repo_name: thautwarm/MLStyle-Playground

max_issues_repo_name: cherichy/MLStyle-Playground

max_forks_repo_name: cherichy/MLStyle-Playground


text:
tuple_type_reverse(::Type{T}) where T<:Tuple = Tuple{reverse(tuple(T.parameters...))...}
tuple_type_reverse(::Type{Tuple{}}) = Tuple{}
tuple_type_reverse(::Type{Tuple{A}}) where A = Tuple{A}
tuple_type_reverse(::Type{Tuple{A,B}}) where {A,B} = Tuple{B,A}

abstract type MemoryLayout end
struct UnknownLayout <: MemoryLayout end
abstract type AbstractStridedLayout <: MemoryLayout end
abstract type AbstractIncreasingStrides <: AbstractStridedLayout end
abstract type AbstractColumnMajor <: AbstractIncreasingStrides end
struct DenseColumnMajor <: AbstractColumnMajor end
struct ColumnMajor <: AbstractColumnMajor end
struct IncreasingStrides <: AbstractIncreasingStrides end
abstract type AbstractDecreasingStrides <: AbstractStridedLayout end
abstract type AbstractRowMajor <: AbstractDecreasingStrides end
struct DenseRowMajor <: AbstractRowMajor end
struct RowMajor <: AbstractRowMajor end
struct DecreasingStrides <: AbstractIncreasingStrides end
struct StridedLayout <: AbstractStridedLayout end
struct ScalarLayout <: MemoryLayout end

"""
    UnknownLayout()

is returned by `MemoryLayout(A)` if it is unknown how the entries of an array `A`
are stored in memory.
"""
UnknownLayout

"""
    AbstractStridedLayout

is an abstract type whose subtypes are returned by `MemoryLayout(A)`
if an array `A` has storage laid out at regular offsets in memory,
and which can therefore be passed to external C and Fortran functions expecting
this memory layout.

Julia's internal linear algebra machinery will automatically (and invisibly)
dispatch to BLAS and LAPACK routines if the memory layout is BLAS compatible and
the element type is a `Float32`, `Float64`, `ComplexF32`, or `ComplexF64`.
In this case, one must implement the strided array interface, which requires
overrides of `strides(A::MyMatrix)` and `unknown_convert(::Type{Ptr{T}}, A::MyMatrix)`.
"""
AbstractStridedLayout

"""
    DenseColumnMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
equivalent to an `Array`, so that `stride(A,1) == 1` and
`stride(A,i) ≡ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`. In particular,
if `A` is a matrix then `strides(A) == `(1, size(A,1))`.

Arrays with `DenseColumnMajor` memory layout must conform to the `DenseArray` interface.
"""
DenseColumnMajor

"""
    ColumnMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a column major array, so that `stride(A,1) == 1` and
`stride(A,i) ≥ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`.

Arrays with `ColumnMajor` memory layout must conform to the `DenseArray` interface.
"""
ColumnMajor

"""
    IncreasingStrides()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a strided array with  increasing strides, so that `stride(A,1) ≥ 1` and
`stride(A,i) ≥ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`.
"""
IncreasingStrides

"""
    DenseRowMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a row major array with dense entries, so that `stride(A,ndims(A)) == 1` and
`stride(A,i) ≡ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`. In particular,
if `A` is a matrix then `strides(A) == `(size(A,2), 1)`.
"""
DenseRowMajor

"""
    RowMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a row major array, so that `stride(A,ndims(A)) == 1` and
stride(A,i) ≥ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`.

If `A` is a matrix  with `RowMajor` memory layout, then
`transpose(A)` should return a matrix whose layout is `ColumnMajor`.
"""
RowMajor

"""
    DecreasingStrides()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a strided array with decreasing strides, so that `stride(A,ndims(A)) ≥ 1` and
stride(A,i) ≥ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`.
"""
DecreasingStrides

"""
    StridedLayout()

is returned by `MemoryLayout(A)` if an array `A` has storage laid out at regular
offsets in memory. `Array`s with `StridedLayout` must conform to the `DenseArray` interface.
"""
StridedLayout

"""
    ScalarLayout()

is returned by `MemoryLayout(A)` if A is a scalar, which does not live in memory
"""
ScalarLayout

"""
    MemoryLayout(A)

specifies the layout in memory for an array `A`. When
you define a new `AbstractArray` type, you can choose to override
`MemoryLayout` to indicate how an array is stored in memory.
For example, if your matrix is column major with `stride(A,2) == size(A,1)`,
then override as follows:

    MemoryLayout(::MyMatrix) = DenseColumnMajor()

The default is `UnknownLayout()` to indicate that the layout
in memory is unknown.

Julia's internal linear algebra machinery will automatically (and invisibly)
dispatch to BLAS and LAPACK routines if the memory layout is compatible.
"""
@inline MemoryLayout(::Type) = UnknownLayout()

@inline MemoryLayout(::Type{<:Number}) = ScalarLayout()
@inline MemoryLayout(::Type{<:DenseArray}) = DenseColumnMajor()

@inline MemoryLayout(::Type{<:ReinterpretArray{T,N,S,P}}) where {T,N,S,P} = reinterpretedlayout(MemoryLayout(P))
@inline reinterpretedlayout(::MemoryLayout) = UnknownLayout()
@inline reinterpretedlayout(::DenseColumnMajor) = DenseColumnMajor()

MemoryLayout(::Type{<:ReshapedArray{T,N,A,DIMS}}) where {T,N,A,DIMS} = reshapedlayout(MemoryLayout(A), DIMS)
@inline reshapedlayout(_, _) = UnknownLayout()
@inline reshapedlayout(::DenseColumnMajor, _) = DenseColumnMajor()

@inline MemoryLayout(A::Type{<:SubArray{T,N,P,I}}) where {T,N,P,I} =
    sublayout(MemoryLayout(P), I)
sublayout(_1, _2) = UnknownLayout()
sublayout(_1, _2, _3)= UnknownLayout()
sublayout(::DenseColumnMajor, ::Type{<:Tuple{<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex}}}) =
    DenseColumnMajor()  # A[:] is DenseColumnMajor if A is DenseColumnMajor
sublayout(ml::AbstractColumnMajor, inds) = _column_sublayout1(ml, inds)
sublayout(::AbstractRowMajor, ::Type{<:Tuple{<:Any}}) =
    UnknownLayout()  # A[:] does not have any structure if A is AbstractRowMajor
sublayout(ml::AbstractRowMajor, inds) = _row_sublayout1(ml, tuple_type_reverse(inds))
sublayout(ml::AbstractStridedLayout, inds) = _strided_sublayout(ml, inds)

_column_sublayout1(::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Union{Int,AbstractCartesianIndex} =
    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
_column_sublayout1(::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Slice =
    DenseColumnMajor() # view(A,:,1,2) is a DenseColumnMajor vector
_column_sublayout1(::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:AbstractUnitRange{Int} =
    DenseColumnMajor() # view(A,1:3,1,2) is a DenseColumnMajor vector
_column_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Union{Int,AbstractCartesianIndex} =
    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
_column_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:AbstractUnitRange{Int} =
    DenseColumnMajor() # view(A,1:3,1,2) is a DenseColumnMajor vector
_column_sublayout1(::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:Slice =
    _column_sublayout(DenseColumnMajor(), DenseColumnMajor(), tuple_type_tail(inds))
_column_sublayout1(par::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:AbstractUnitRange{Int} =
    _column_sublayout(par, ColumnMajor(), tuple_type_tail(inds))
_column_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:AbstractUnitRange{Int} =
    _column_sublayout(par, ColumnMajor(), tuple_type_tail(inds))
_column_sublayout1(par::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _column_sublayout(par, StridedLayout(), tuple_type_tail(inds))
_column_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _column_sublayout(par, StridedLayout(), tuple_type_tail(inds))
_column_sublayout1(par, inds) = UnknownLayout()
_column_sublayout(par, ret, ::Type{<:Tuple{}}) = ret
_column_sublayout(par, ret, ::Type{<:Tuple{I}}) where I = UnknownLayout()
_column_sublayout(::DenseColumnMajor, ::DenseColumnMajor, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    DenseColumnMajor() # A[:,1:3,1,2] is DenseColumnMajor if A is DenseColumnMajor
_column_sublayout(par::DenseColumnMajor, ::DenseColumnMajor, inds::Type{<:Tuple{I, Vararg{Int}}}) where I<:Slice =
    DenseColumnMajor()
_column_sublayout(par::DenseColumnMajor, ::DenseColumnMajor, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Slice =
    _column_sublayout(par, DenseColumnMajor(), tuple_type_tail(inds))
_column_sublayout(par, ::AbstractColumnMajor, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    _column_sublayout(par, ColumnMajor(), tuple_type_tail(inds))
_column_sublayout(par, ::AbstractStridedLayout, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _column_sublayout(par, StridedLayout(), tuple_type_tail(inds))

_row_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Union{Int,AbstractCartesianIndex} =
    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
_row_sublayout1(::DenseRowMajor, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Slice =
    DenseColumnMajor() # view(A,1,2,:) is a DenseColumnMajor vector
_row_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:AbstractUnitRange{Int} =
    DenseColumnMajor() # view(A,1,2,1:3) is a DenseColumnMajor vector
_row_sublayout1(::DenseRowMajor, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:Slice =
    _row_sublayout(DenseRowMajor(), DenseRowMajor(), tuple_type_tail(inds))
_row_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:AbstractUnitRange{Int} =
    _row_sublayout(par, RowMajor(), tuple_type_tail(inds))
_row_sublayout1(par, inds::Type{<:Tuple{I,Vararg{Any}}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _row_sublayout(par, StridedLayout(), tuple_type_tail(inds))
_row_sublayout1(par, inds) = UnknownLayout()
_row_sublayout(par, ret, ::Type{<:Tuple{}}) = ret
_row_sublayout(par, ret, ::Type{<:Tuple{I}}) where I = UnknownLayout()
_row_sublayout(::DenseRowMajor, ::DenseRowMajor, inds::Type{<:Tuple{I,Vararg{Int}}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    DenseRowMajor() # A[1,2,1:3,:] is DenseRowMajor if A is DenseRowMajor
_row_sublayout(par::DenseRowMajor, ::DenseRowMajor, inds::Type{<:Tuple{I, Vararg{Int}}}) where I<:Slice =
    DenseRowMajor()
_row_sublayout(par::DenseRowMajor, ::DenseRowMajor, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Slice =
    _row_sublayout(par, DenseRowMajor(), tuple_type_tail(inds))
_row_sublayout(par::AbstractRowMajor, ::AbstractRowMajor, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    _row_sublayout(par, RowMajor(), tuple_type_tail(inds))
_row_sublayout(par::AbstractRowMajor, ::AbstractStridedLayout, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _row_sublayout(par, StridedLayout(), tuple_type_tail(inds))

_strided_sublayout(par, inds) = UnknownLayout()
_strided_sublayout(par, ::Type{<:Tuple{}}) = StridedLayout()
_strided_sublayout(par, inds::Type{<:Tuple{I, Vararg{Any}}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _strided_sublayout(par, tuple_type_tail(inds))

# MemoryLayout of transposed and adjoint matrices
struct ConjLayout{ML<:MemoryLayout} <: MemoryLayout end

conjlayout(_1, _2) = UnknownLayout()
conjlayout(::Type{<:Complex}, ::ConjLayout{ML}) where ML = ML()
conjlayout(::Type{<:Complex}, ::ML) where ML<:AbstractStridedLayout = ConjLayout{ML}()
conjlayout(::Type{<:Real}, M::MemoryLayout) = M

sublayout(::ConjLayout{ML}, t::Type{<:Tuple}) where ML = ConjLayout{typeof(sublayout(ML(), t))}()

MemoryLayout(::Type{Transpose{T,P}}) where {T,P} = transposelayout(MemoryLayout(P))
MemoryLayout(::Type{Adjoint{T,P}}) where {T,P} = adjointlayout(T, MemoryLayout(P))
transposelayout(_) = UnknownLayout()
transposelayout(::StridedLayout) = StridedLayout()
transposelayout(::ColumnMajor) = RowMajor()
transposelayout(::RowMajor) = ColumnMajor()
transposelayout(::DenseColumnMajor) = DenseRowMajor()
transposelayout(::DenseRowMajor) = DenseColumnMajor()
transposelayout(::ConjLayout{ML}) where ML = ConjLayout{typeof(transposelayout(ML()))}()
adjointlayout(::Type{T}, M::MemoryLayout) where T = transposelayout(conjlayout(T, M))

# MemoryLayout of Symmetric/Hermitian
"""
    SymmetricLayout{layout}()

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
as a symmetrized version of `layout`, where the entries used are dictated by the
`uplo`, which can be `'U'` or `L'`.

A matrix that has memory layout `SymmetricLayout(layout, uplo)` must overrided
`symmetricdata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] == B[k,j]` for `j ≥ k` if `uplo == 'U'` (`j ≤ k` if `uplo == 'L'`) and
`A[k,j] == B[j,k]` for `j < k` if `uplo == 'U'` (`j > k` if `uplo == 'L'`).
"""
struct SymmetricLayout{ML<:MemoryLayout} <: MemoryLayout end

"""
    HermitianLayout(layout, uplo)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
as a hermitianized version of `layout`, where the entries used are dictated by the
`uplo`, which can be `'U'` or `L'`.

A matrix that has memory layout `HermitianLayout(layout, uplo)` must overrided
`hermitiandata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] == B[k,j]` for `j ≥ k` if `uplo == 'U'` (`j ≤ k` if `uplo == 'L'`) and
`A[k,j] == conj(B[j,k])` for `j < k` if `uplo == 'U'` (`j > k` if `uplo == 'L'`).
"""
struct HermitianLayout{ML<:MemoryLayout} <: MemoryLayout end

MemoryLayout(::Type{Hermitian{T,P}}) where {T,P} = hermitianlayout(T, MemoryLayout(P))
MemoryLayout(::Type{Symmetric{T,P}}) where {T,P} = symmetriclayout(MemoryLayout(P))
hermitianlayout(_1, _2) = UnknownLayout()
hermitianlayout(::Type{<:Complex}, ::ML) where ML<:AbstractColumnMajor = HermitianLayout{ML}()
hermitianlayout(::Type{<:Real}, ::ML) where ML<:AbstractColumnMajor = SymmetricLayout{ML}()
hermitianlayout(::Type{<:Complex}, ::ML) where ML<:AbstractRowMajor = HermitianLayout{ML}()
hermitianlayout(::Type{<:Real}, ::ML) where ML<:AbstractRowMajor = SymmetricLayout{ML}()
symmetriclayout(_1) = UnknownLayout()
symmetriclayout(::ML) where ML<:AbstractColumnMajor = SymmetricLayout{ML}()
symmetriclayout(::ML) where ML<:AbstractRowMajor = SymmetricLayout{ML}()
transposelayout(S::SymmetricLayout) = S
adjointlayout(::Type{T}, S::SymmetricLayout) where T<:Real = S
adjointlayout(::Type{T}, S::HermitianLayout) where T = S
sublayout(S::SymmetricLayout, ::Type{<:Tuple{<:Slice,<:Slice}}) = S
sublayout(S::HermitianLayout, ::Type{<:Tuple{<:Slice,<:Slice}}) = S

symmetricdata(A::Symmetric) = A.data
symmetricdata(A::Hermitian{<:Real}) = A.data
symmetricdata(V::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = symmetricdata(parent(V))
symmetricdata(V::Adjoint{<:Real}) = symmetricdata(parent(V))
symmetricdata(V::Transpose) = symmetricdata(parent(V))
hermitiandata(A::Hermitian) = A.data
hermitiandata(V::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = hermitiandata(parent(V))
hermitiandata(V::Adjoint) = hermitiandata(parent(V))
hermitiandata(V::Transpose{<:Real}) = hermitiandata(parent(V))

symmetricuplo(A::Symmetric) = A.uplo
symmetricuplo(A::Hermitian) = A.uplo
symmetricuplo(A::Adjoint) = symmetricuplo(parent(A))
symmetricuplo(A::Transpose) = A.uplo
symmetricuplo(A::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = symmetricuplo(parent(A))

# MemoryLayout of triangular matrices
struct TriangularLayout{UPLO,UNIT,ML} <: MemoryLayout end

"""
    LowerTriangularLayout(layout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `LowerTriangular(B)` where `B` satisfies `MemoryLayout(B) == layout`.

A matrix that has memory layout `LowerTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ zero(eltype(A))` for `j > k` and
`A[k,j] ≡ B[k,j]` for `j ≤ k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `UpperTriangularLayout`.
"""
const LowerTriangularLayout{ML} = TriangularLayout{'L','N',ML}

"""
    UnitLowerTriangularLayout(ML::MemoryLayout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `UnitLowerTriangular(B)` where `B` satisfies `MemoryLayout(B) == layout`.

A matrix that has memory layout `UnitLowerTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ zero(eltype(A))` for `j > k`,
`A[k,j] ≡ one(eltype(A))` for `j == k`,
`A[k,j] ≡ B[k,j]` for `j < k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `UnitUpperTriangularLayout`.
"""
const UnitLowerTriangularLayout{ML} = TriangularLayout{'L','U',ML}

"""
    UpperTriangularLayout(ML::MemoryLayout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `UpperTriangularLayout(B)` where `B` satisfies `MemoryLayout(B) == ML`.

A matrix that has memory layout `UpperTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ B[k,j]` for `j ≥ k` and
`A[k,j] ≡ zero(eltype(A))` for `j < k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `LowerTriangularLayout`.
"""
const UpperTriangularLayout{ML} = TriangularLayout{'U','N',ML}

"""
    UnitUpperTriangularLayout(ML::MemoryLayout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `UpperTriangularLayout(B)` where `B` satisfies `MemoryLayout(B) == ML`.

A matrix that has memory layout `UnitUpperTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ B[k,j]` for `j > k`,
`A[k,j] ≡ one(eltype(A))` for `j == k`,
`A[k,j] ≡ zero(eltype(A))` for `j < k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `UnitLowerTriangularLayout`.
"""
UnitUpperTriangularLayout{ML} = TriangularLayout{'U','U',ML}

MemoryLayout(A::Type{UpperTriangular{T,P}}) where {T,P} = triangularlayout(UpperTriangularLayout, MemoryLayout(P))
MemoryLayout(A::Type{UnitUpperTriangular{T,P}}) where {T,P} = triangularlayout(UnitUpperTriangularLayout, MemoryLayout(P))
MemoryLayout(A::Type{LowerTriangular{T,P}}) where {T,P} = triangularlayout(LowerTriangularLayout, MemoryLayout(P))
MemoryLayout(A::Type{UnitLowerTriangular{T,P}}) where {T,P} = triangularlayout(UnitLowerTriangularLayout, MemoryLayout(P))
triangularlayout(_::Type{Tri}, ::MemoryLayout) where Tri = Tri{UnknownLayout}()
triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AbstractColumnMajor} = Tri{ML}()
triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AbstractRowMajor} = Tri{ML}()
triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:ConjLayout{<:AbstractRowMajor}} = Tri{ML}()
sublayout(layout::TriangularLayout, ::Type{<:Tuple{<:Union{Slice,Base.OneTo},<:Union{Slice,Base.OneTo}}}) = layout
conjlayout(::Type{<:Complex}, ::TriangularLayout{UPLO,UNIT,ML}) where {UPLO,UNIT,ML} =
    TriangularLayout{UPLO,UNIT,ConjLayout{ML}}()

for (TriLayout, TriLayoutTrans) in ((UpperTriangularLayout,     LowerTriangularLayout),
                                    (UnitUpperTriangularLayout, UnitLowerTriangularLayout),
                                    (LowerTriangularLayout,     UpperTriangularLayout),
                                    (UnitLowerTriangularLayout, UnitUpperTriangularLayout))
    @eval transposelayout(::$TriLayout{ML}) where ML = $TriLayoutTrans{typeof(transposelayout(ML()))}()
end

triangulardata(A::AbstractTriangular) = parent(A)
triangulardata(A::Adjoint) = Adjoint(triangulardata(parent(A)))
triangulardata(A::Transpose) = Transpose(triangulardata(parent(A)))
triangulardata(A::SubArray{<:Any,2,<:Any,<:Tuple{<:Union{Slice,Base.OneTo},<:Union{Slice,Base.OneTo}}}) =
    view(triangulardata(parent(A)), parentindices(A)...)

abstract type AbstractBandedLayout <: MemoryLayout end

struct DiagonalLayout{ML} <: AbstractBandedLayout end
struct SymTridiagonalLayout{ML} <: AbstractBandedLayout end

diagonallayout(_) = DiagonalLayout{UnknownLayout}()
diagonallayout(::ML) where ML<:AbstractStridedLayout = DiagonalLayout{ML}()
MemoryLayout(D::Type{Diagonal{T,P}}) where {T,P} = diagonallayout(MemoryLayout(P))
diagonaldata(D::Diagonal) = parent(D)

MemoryLayout(::Type{SymTridiagonal{T,P}}) where {T,P} = SymTridiagonalLayout{typeof(MemoryLayout(P))}()
diagonaldata(D::SymTridiagonal) = D.dv
offdiagonaldata(D::SymTridiagonal) = D.ev

transposelayout(ml::DiagonalLayout) = ml
transposelayout(ml::SymTridiagonalLayout) = ml
transposelayout(ml::ConjLayout{DiagonalLayout}) = ml

adjointlayout(::Type{<:Real}, ml::SymTridiagonalLayout) = ml

###
# Fill
####
abstract type AbstractFillLayout <: MemoryLayout end
struct FillLayout <: AbstractFillLayout end
struct ZerosLayout <: AbstractFillLayout end
struct EyeLayout <: MemoryLayout end

MemoryLayout(::Type{<:AbstractFill}) = FillLayout()
MemoryLayout(::Type{<:Zeros}) = ZerosLayout()
diagonallayout(::ML) where ML<:AbstractFillLayout = DiagonalLayout{ML}()
# all sub arrays are same
sublayout(L::AbstractFillLayout, inds::Type) = L
reshapedlayout(L::AbstractFillLayout, _) = L
adjointlayout(::Type, L::AbstractFillLayout) = L
transposelayout(L::AbstractFillLayout) = L

rowsupport(_, A, k) = axes(A,2)
""""
    rowsupport(A, k)

gives an iterator containing the possible non-zero entries in the k-th row of A.
"""
rowsupport(A, k) = rowsupport(MemoryLayout(typeof(A)), A, k)
rowsupport(A) = rowsupport(A, axes(A,1))

colsupport(_, A, j) = axes(A,1)

""""
    colsupport(A, j)

gives an iterator containing the possible non-zero entries in the j-th column of A.
"""
colsupport(A, j) = colsupport(MemoryLayout(typeof(A)), A, j)
colsupport(A) = colsupport(A, axes(A,2))

rowsupport(::ZerosLayout, A, _) = 1:0
colsupport(::ZerosLayout, A, _) = 1:0


text:
using ..CUDASupport
using ..TopOpt: @init_cuda
@init_cuda()
using ..GPUUtils
import ..TopOpt: whichdevice
using ..TopOpt: GPU

### Compliance

whichdevice(c::Compliance) = whichdevice(c.cell_comp)
whichdevice(c::MeanCompliance) = whichdevice(c.compliance)

hutch_rand!(v::CuArray) = v .= round.(CuArrays.CURAND.curand(size(v))) .* 2 .- 1

function Compliance(::GPU, problem::StiffnessTopOptProblem{dim, T}, solver::AbstractDisplacementSolver, ::Type{TI}=Int; tracing = false, logarithm = false, maxfevals = 10^8) where {dim, T, TI}
    comp = T(0)
    cell_comp = zeros(CuVector{T}, getncells(problem.ch.dh.grid))
    grad = CuVector(fill(T(NaN), length(cell_comp) - sum(problem.black) - sum(problem.white)))
    topopt_trace = TopOptTrace{T,TI}()
    reuse = false
    fevals = TI(0)
    return Compliance(problem, solver, comp, cell_comp, grad, tracing, topopt_trace, reuse, fevals, logarithm, maxfevals)
end
@define_cu(Compliance, :solver, :cell_comp, :grad, :cheqfilter)

function compute_compliance(cell_comp::CuVector{T}, grad, cell_dofs, Kes, u,
                            black, white, varind, x, penalty, xmin) where {T}

    args = (cell_comp, grad, cell_dofs, Kes, u, black, white, varind, x, penalty, xmin)
    callkernel(dev, comp_kernel1, args)
    CUDAdrv.synchronize(ctx)
    obj = compute_obj(cell_comp, x, varind, black, white, penalty, xmin)

    return obj
end

# CUDA kernels
function comp_kernel1(cell_comp::AbstractVector{T}, grad, cell_dofs, Kes, u,
                        black, white, varind, x, penalty, xmin) where {T}

    i = @thread_global_index()
    offset = @total_threads()
    @inbounds while i <= length(cell_comp)
        cell_comp[i] = zero(T)
        Ke = rawmatrix(Kes[i])
        for w in 1:size(Ke, 2)
            for v in 1:size(Ke, 1)
                if Ke isa Symmetric
                    cell_comp[i] += u[cell_dofs[v,i]]*Ke.data[v,w]*u[cell_dofs[w,i]]
                else
                    cell_comp[i] += u[cell_dofs[v,i]]*Ke[v,w]*u[cell_dofs[w,i]]
                end
            end
        end
        if !(black[i] || white[i])
            d = ForwardDiff.Dual{T}(x[varind[i]], one(T))
            if PENALTY_BEFORE_INTERPOLATION
                p = density(penalty(d), xmin)
            else
                p = penalty(density(d, xmin))
            end
            grad[varind[i]] = -p.partials[1] * cell_comp[i]
        end

        i += offset
    end
    return
end

function compute_obj(cell_comp::AbstractVector{T}, x, varind, black, white, penalty, xmin, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    result = similar(cell_comp, T, (blocksize,))
    args = (result, cell_comp, x, varind, black, white, penalty, xmin, Val(threads))
    @cuda blocks = blocksize threads = threads comp_kernel2(args...)
    CUDAnative.synchronize()
    obj = reduce(+, Array(result))
    return obj
end

function comp_kernel2(result, cell_comp::AbstractVector{T}, x, varind, black, white, penalty, xmin, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(i, length(cell_comp), +, T, LMEM, result, begin
        w_comp(cell_comp[i], x[varind[i]], black[i], white[i], penalty, xmin)
    end)

    return
end
@inline function w_comp(comp::T, x, black, white, penalty, xmin) where {T}
    if PENALTY_BEFORE_INTERPOLATION
        return ifelse(black, comp,
        ifelse(white, xmin * comp,
                        (d = ForwardDiff.Dual{T}(x, one(T));
                        p = density(penalty(d), xmin); p.value * comp)))
    else
        return ifelse(black, comp,
        ifelse(white, xmin * comp,
                        (d = ForwardDiff.Dual{T}(x, one(T));
                        p = penalty(density(d, xmin)); p.value * comp)))
    end
end

### Volume

whichdevice(v::Volume) = whichdevice(v.cellvolumes)
@define_cu(Volume, :cellvolumes, :grad, :problem) # should be optimized to avoid replicating problem

function compute_volume(cellvolumes::CuVector{T}, x, fixed_volume, varind, black, white, ::Val{blocksize} = Val(80), ::Val{threads} = Val(256)) where {T, blocksize, threads}
    result = similar(cellvolumes, T, (blocksize,))
    args = (result, cellvolumes, x, varind, black, white, Val(threads))
    @cuda blocks = blocksize threads = threads volume_kernel(args...)
    CUDAnative.synchronize()
    vol = reduce(+, Array(result)) + fixed_volume
    return vol
end

function volume_kernel(result, cellvolumes::AbstractVector{T}, x, varind, black, white, ::Val{LMEM}) where {T, LMEM}
    @mapreduce_block(i, length(cellvolumes), +, T, LMEM, result, begin
        if !(black[i]) && !(white[i])
            x[varind[i]]*cellvolumes[i]
        else
            zero(T)
        end
    end)

    return
end


text:
module BIDSTools

using JSON
using DataStructures
using CSV
using DataFrames
using PrettyTables

include("File.jl")
include("Session.jl")
include("Subject.jl")
include("Layout.jl")
# Class
export
    Layout,
    Subject,
    Session,
    File

# Utility functions
export
    total_subjects,
    total_sessions,
    total_files,
    get_metadata_path,
    parse_path,
    parse_fname,
    print_dataset_description,
    list_subject_detail,
    list_scans_detail,
    get_files,
    get_sub,
    get_ses,
    construct_fname

end # module


text:
function _function(
    ::Any,
    ::Type{MOI.VariableIndex},
    x::Vector{MOI.VariableIndex},
)
    return x[1]
end

function _function(
    ::Any,
    ::Type{MOI.VectorOfVariables},
    x::Vector{MOI.VariableIndex},
)
    return MOI.VectorOfVariables(x)
end

function _function(
    ::Type{T},
    ::Type{MOI.ScalarAffineFunction},
    x::Vector{MOI.VariableIndex},
) where {T}
    return MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(one(T), x), zero(T))
end

function _function(
    ::Type{T},
    ::Type{MOI.ScalarQuadraticFunction},
    x::Vector{MOI.VariableIndex},
) where {T}
    return MOI.ScalarQuadraticFunction(
        MOI.ScalarQuadraticTerm.(one(T), x, x),
        MOI.ScalarAffineTerm.(one(T), x),
        zero(T),
    )
end

function _function(
    ::Type{T},
    ::Type{MOI.VectorAffineFunction},
    x::Vector{MOI.VariableIndex},
) where {T}
    return MOI.VectorAffineFunction(
        MOI.VectorAffineTerm.(1:length(x), MOI.ScalarAffineTerm.(one(T), x)),
        zeros(T, length(x)),
    )
end

function _function(
    ::Type{T},
    ::Type{MOI.VectorQuadraticFunction},
    x::Vector{MOI.VariableIndex},
) where {T}
    return MOI.VectorQuadraticFunction(
        MOI.VectorQuadraticTerm.(
            1:length(x),
            MOI.ScalarQuadraticTerm.(one(T), x, x),
        ),
        MOI.VectorAffineTerm.(1:length(x), MOI.ScalarAffineTerm.(one(T), x)),
        zeros(T, length(x)),
    )
end

# Default fallback.
_set(::Any, ::Type{S}) where {S} = _set(S)

_set(::Type{T}, ::Type{MOI.LessThan}) where {T} = MOI.LessThan(one(T))
_set(::Type{T}, ::Type{MOI.GreaterThan}) where {T} = MOI.GreaterThan(one(T))
_set(::Type{T}, ::Type{MOI.EqualTo}) where {T} = MOI.EqualTo(one(T))
_set(::Type{T}, ::Type{MOI.Interval}) where {T} = MOI.Interval(zero(T), one(T))
_set(::Type{MOI.ZeroOne}) = MOI.ZeroOne()
_set(::Type{MOI.Integer}) = MOI.Integer()
function _set(::Type{T}, ::Type{MOI.Semicontinuous}) where {T}
    return MOI.Semicontinuous(zero(T), one(T))
end
function _set(::Type{T}, ::Type{MOI.Semiinteger}) where {T}
    return MOI.Semiinteger(zero(T), one(T))
end
_set(::Type{T}, ::Type{MOI.SOS1}) where {T} = MOI.SOS1(convert.(T, 1:2))
_set(::Type{T}, ::Type{MOI.SOS2}) where {T} = MOI.SOS2(convert.(T, 1:2))
_set(::Type{MOI.Zeros}) = MOI.Zeros(2)
_set(::Type{MOI.Nonpositives}) = MOI.Nonpositives(2)
_set(::Type{MOI.Nonnegatives}) = MOI.Nonnegatives(2)
_set(::Type{MOI.NormInfinityCone}) = MOI.NormInfinityCone(3)
_set(::Type{MOI.NormOneCone}) = MOI.NormOneCone(3)
_set(::Type{MOI.SecondOrderCone}) = MOI.SecondOrderCone(3)
_set(::Type{MOI.RotatedSecondOrderCone}) = MOI.RotatedSecondOrderCone(3)
_set(::Type{MOI.GeometricMeanCone}) = MOI.GeometricMeanCone(3)
_set(::Type{MOI.ExponentialCone}) = MOI.ExponentialCone()
_set(::Type{MOI.DualExponentialCone}) = MOI.DualExponentialCone()
_set(::Type{MOI.PowerCone}) = MOI.PowerCone(0.5)
_set(::Type{MOI.DualPowerCone}) = MOI.DualPowerCone(0.5)
_set(::Type{MOI.RelativeEntropyCone}) = MOI.RelativeEntropyCone(3)
_set(::Type{MOI.NormSpectralCone}) = MOI.NormSpectralCone(2, 3)
_set(::Type{MOI.NormNuclearCone}) = MOI.NormNuclearCone(2, 3)
function _set(::Type{MOI.PositiveSemidefiniteConeTriangle})
    return MOI.PositiveSemidefiniteConeTriangle(3)
end
function _set(::Type{MOI.PositiveSemidefiniteConeSquare})
    return MOI.PositiveSemidefiniteConeSquare(3)
end
_set(::Type{MOI.LogDetConeTriangle}) = MOI.LogDetConeTriangle(3)
_set(::Type{MOI.LogDetConeSquare}) = MOI.LogDetConeSquare(3)
_set(::Type{MOI.RootDetConeTriangle}) = MOI.RootDetConeTriangle(3)
_set(::Type{MOI.RootDetConeSquare}) = MOI.RootDetConeSquare(3)
_set(::Type{MOI.Complements}) = MOI.Complements(2)

function _set(::Type{MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.LessThan{Float64}}})
    return MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.LessThan(3.0))
end

function _set(
    ::Type{MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.GreaterThan{Float64}}},
)
    return MOI.Indicator{MOI.ACTIVATE_ON_ONE}(MOI.GreaterThan(3.0))
end

function _basic_constraint_test_helper(
    model::MOI.ModelLike,
    config::Config{T},
    ::Type{UntypedF},
    ::Type{UntypedS},
) where {T,UntypedF,UntypedS}
    set = _set(T, UntypedS)
    N = MOI.dimension(set)
    x = MOI.add_variables(model, N)
    constraint_function = _function(T, UntypedF, x)
    F, S = typeof(constraint_function), typeof(set)
    ###
    ### Test MOI.supports_constraint
    ###
    @requires MOI.supports_constraint(model, F, S)
    ###
    ### Test MOI.NumberOfConstraints
    ###
    @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 0
    c = MOI.add_constraint(model, constraint_function, set)
    @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 1
    _test_attribute_value_type(model, MOI.NumberOfConstraints{F,S}())
    ###
    ### Test MOI.is_valid
    ###
    @test MOI.is_valid(model, c)
    ###
    ### Test MOI.ConstraintName
    ###
    if _supports(config, MOI.ConstraintName)
        if F == MOI.VariableIndex
            @test_throws(
                MOI.VariableIndexConstraintNameError(),
                MOI.supports(model, MOI.ConstraintName(), typeof(c)),
            )
            @test_throws(
                MOI.VariableIndexConstraintNameError(),
                MOI.set(model, MOI.ConstraintName(), c, "c"),
            )
        else
            @test MOI.get(model, MOI.ConstraintName(), c) == ""
            @test MOI.supports(model, MOI.ConstraintName(), typeof(c))
            MOI.set(model, MOI.ConstraintName(), c, "c")
            @test MOI.get(model, MOI.ConstraintName(), c) == "c"
            _test_attribute_value_type(model, MOI.ConstraintName(), c)
        end
    end
    ###
    ### Test MOI.ConstraintFunction
    ###
    if _supports(config, MOI.ConstraintFunction)
        # Don't compare directly, because `f` might not be canonicalized.
        f = MOI.get(model, MOI.ConstraintFunction(), c)
        @test isapprox(MOI.Utilities.canonical(f), constraint_function, config)
        @test isapprox(
            MOI.get(model, MOI.CanonicalConstraintFunction(), c),
            constraint_function,
            config,
        )
        _test_attribute_value_type(model, MOI.ConstraintFunction(), c)
        _test_attribute_value_type(model, MOI.CanonicalConstraintFunction(), c)
    end
    ###
    ### Test MOI.ConstraintSet
    ###
    if _supports(config, MOI.ConstraintSet)
        @test MOI.get(model, MOI.ConstraintSet(), c) == set
        _test_attribute_value_type(model, MOI.ConstraintSet(), c)
    end
    ###
    ### Test MOI.ListOfConstraintIndices
    ###
    c_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
    @test length(c_indices) == 1
    ###
    ### Test MOI.add_constraints
    ###
    if F != MOI.VariableIndex && F != MOI.VectorOfVariables
        # We can't add multiple variable constraints as these are
        # interpreted as bounds etc.
        MOI.add_constraints(
            model,
            [constraint_function, constraint_function],
            [set, set],
        )
        @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) == 3
        @test length(MOI.get(model, MOI.ListOfConstraintIndices{F,S}())) == 3
        c_indices = MOI.get(model, MOI.ListOfConstraintIndices{F,S}())
        @test all(MOI.is_valid.(model, c_indices))
    end
    ###
    ### Test MOI.delete
    ###
    if _supports(config, MOI.delete)
        MOI.delete(model, c_indices[1])
        @test MOI.get(model, MOI.NumberOfConstraints{F,S}()) ==
              length(c_indices) - 1
        @test !MOI.is_valid(model, c_indices[1])
        @test_throws(
            MOI.InvalidIndex(c_indices[1]),
            MOI.delete(model, c_indices[1]),
        )
        if _supports(config, MOI.ConstraintFunction)
            @test_throws(
                MOI.InvalidIndex(c_indices[1]),
                MOI.get(model, MOI.ConstraintFunction(), c_indices[1]),
            )
        end
        if _supports(config, MOI.ConstraintSet)
            @test_throws(
                MOI.InvalidIndex(c_indices[1]),
                MOI.get(model, MOI.ConstraintSet(), c_indices[1]),
            )
        end
    end
    return
end

for s in [
    :GreaterThan,
    :LessThan,
    :EqualTo,
    :Interval,
    :ZeroOne,
    :Semicontinuous,
    :Semiinteger,
    :SOS1,
    :SOS2,
    :Zeros,
    :Nonpositives,
    :Nonnegatives,
    :NormInfinityCone,
    :NormOneCone,
    :SecondOrderCone,
    :RotatedSecondOrderCone,
    :GeometricMeanCone,
    :ExponentialCone,
    :DualExponentialCone,
    :PowerCone,
    :DualPowerCone,
    :RelativeEntropyCone,
    :NormSpectralCone,
    :NormNuclearCone,
    :PositiveSemidefiniteConeSquare,
    :PositiveSemidefiniteConeTriangle,
    :LogDetConeTriangle,
    :LogDetConeSquare,
    :RootDetConeTriangle,
    :RootDetConeSquare,
    :Complements,
]
    S = getfield(MOI, s)
    functions = if S <: MOI.AbstractScalarSet
        (:VariableIndex, :ScalarAffineFunction, :ScalarQuadraticFunction)
    else
        (:VectorOfVariables, :VectorAffineFunction, :VectorQuadraticFunction)
    end
    for f in functions
        func = Symbol("test_basic_$(f)_$(s)")
        F = getfield(MOI, f)
        @eval begin
            function $(func)(model::MOI.ModelLike, config::Config)
                _basic_constraint_test_helper(model, config, $F, $S)
                return
            end
        end
    end
end

function test_basic_VectorAffineFunction_Indicator_LessThan(
    model::MOI.ModelLike,
    config::Config,
)
    _basic_constraint_test_helper(
        model,
        config,
        MOI.VectorAffineFunction,
        MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.LessThan{Float64}},
    )
    return
end

function test_basic_VectorAffineFunction_Indicator_GreaterThan(
    model::MOI.ModelLike,
    config::Config,
)
    _basic_constraint_test_helper(
        model,
        config,
        MOI.VectorAffineFunction,
        MOI.Indicator{MOI.ACTIVATE_ON_ONE,MOI.GreaterThan{Float64}},
    )
    return
end


text:
@testset "Title Function Tests" begin

  @test isdefined(Julz, :title) == true

  @test_throws MethodError Julz.title(1)

  @test Julz.title("first") == "First"

  @test Julz.title("first class") == "First Class"

  @test Julz.title("first_class") == "First Class"

  @test Julz.title("first-class") == "First-class"

end


text:
using FourierFlows, PyPlot, JLD2, Printf, Random, FFTW

using Random: seed!

import GeophysicalFlows.TwoDTurb
import GeophysicalFlows.TwoDTurb: energy, enstrophy
import GeophysicalFlows: peakedisotropicspectrum

dev = CPU()     # Device (CPU/GPU)

# Parameters
  n = 256
  L = 2π
 nν = 2
  ν = 0.0
 dt = 1e-2
nsteps = 5000
 nsubs = 200

# Files
filepath = "."
plotpath = "./plots_McWilliams1984"
plotname = "snapshots"
filename = joinpath(filepath, "McWilliams1984.jld2")

# File management
if isfile(filename); rm(filename); end
if !isdir(plotpath); mkdir(plotpath); end

# Initialize problem
prob = TwoDTurb.Problem(; nx=n, Lx=L, ny=n, Ly=L, ν=ν, nν=nν, dt=dt, stepper="FilteredRK4", dev=dev)

sol, cl, vs, gr, filter = prob.sol, prob.clock, prob.vars, prob.grid, prob.timestepper.filter
x, y = gridpoints(gr)

# Initial condition closely following pyqg barotropic example
# that reproduces the results of the paper by McWilliams (1984)
seed!(1234)
k0, E0 = 6, 0.5
zetai  = peakedisotropicspectrum(gr, k0, E0, mask=filter)
TwoDTurb.set_zeta!(prob, zetai)

# Create Diagnostic -- energy and enstrophy are functions imported at the top.
E = Diagnostic(energy, prob; nsteps=nsteps)
Z = Diagnostic(enstrophy, prob; nsteps=nsteps)
diags = [E, Z] # A list of Diagnostics types passed to "stepforward!" will
# be updated every timestep.

# Create Output
get_sol(prob) = Array(prob.sol) # extracts the Fourier-transformed solution
get_u(prob) = Array(irfft(im*gr.l.*gr.invKrsq.*sol, gr.nx))
out = Output(prob, filename, (:sol, get_sol), (:u, get_u))
saveproblem(out)

function plot_output(prob, fig, axs; drawcolorbar=false)
  # Plot the vorticity field and the evolution of energy and enstrophy.
  TwoDTurb.updatevars!(prob)
  sca(axs[1])
  pcolormesh(x, y, vs.zeta)
  clim(-40, 40)
  axis("off")
  axis("square")
  if drawcolorbar==true
    colorbar()
  end

  sca(axs[2])
  cla()
  plot(E.t[1:E.i], E.data[1:E.i]/E.data[1])
  plot(Z.t[1:Z.i], Z.data[1:E.i]/Z.data[1])
  xlabel(L"t")
  ylabel(L"\Delta E, \, \Delta Z")

  pause(0.01)
end

# Step forward

fig, axs = subplots(ncols=2, nrows=1, figsize=(12, 4))
plot_output(prob, fig, axs; drawcolorbar=true)

startwalltime = time()
while cl.step < nsteps
  stepforward!(prob, diags, nsubs)
  saveoutput(out)

  # Message
  log = @sprintf("step: %04d, t: %d, ΔE: %.4f, ΔZ: %.4f, τ: %.2f min",
    cl.step, cl.t, E.data[E.i]/E.data[1], Z.data[Z.i]/Z.data[1], (time()-startwalltime)/60)

  println(log)
  plot_output(prob, fig, axs; drawcolorbar=false)
end
println("finished")
plot_output(prob, fig, axs; drawcolorbar=false)

savename = @sprintf("%s_%09d.png", joinpath(plotpath, plotname), cl.step)
savefig(savename, dpi=240)


text:
module Flagline

using ..Ahorn, Maple
using Random

const placements = Ahorn.PlacementDict(
    "Clothes Line" => Ahorn.EntityPlacement(
        Maple.ClothesLine,
        "line"
    ),
    "Cliffside Flags" => Ahorn.EntityPlacement(
        Maple.CliffFlags,
        "line"
    ),
)

const flagLineUnion = Union{Maple.ClothesLine, Maple.CliffFlags}

Ahorn.nodeLimits(entity::flagLineUnion) = 1, 1

function Ahorn.selection(entity::flagLineUnion)
    nodes = get(entity.data, "nodes", ())
    x, y = Ahorn.position(entity)

    res = Ahorn.Rectangle[Ahorn.Rectangle(x - 4, y - 4, 8, 8)]

    for node in nodes
        nx, ny = node

        push!(res, Ahorn.Rectangle(nx - 4, ny - 4, 8, 8))
    end

    return res
end

const droopAmount = 0.6
const selectedColor = Ahorn.colors.selection_selected_fc

const clothesLineColor = (128, 128, 163, 1) ./ (255, 255, 255, 1)
const clothesPinColor = (128, 128, 128, 1) ./ (255, 255, 255, 1)
const clothesColors = Ahorn.colorTupleType[
    (13, 46, 107, 1) ./ (255, 255, 255, 1),
    (61, 38, 136, 1) ./ (255, 255, 255, 1),
    (79, 110, 157, 1) ./ (255, 255, 255, 1),
    (71, 25, 74, 1) ./ (255, 255, 255, 1)
]
const clothesMinFlagHeight = 8
const clothesMaxFlagHeight = 20
const clothesMinFlagLength = 8
const clothesMaxFlagLength = 16
const clothesMinSpace = 2
const clothesMaxSpace = 8

const cliffLineColor = (128, 128, 163, 1) ./ (255, 255, 255, 1)
const cliffPinColor = (128, 128, 128, 1) ./ (255, 255, 255, 1)
const cliffColors = Ahorn.colorTupleType[
    (216, 95, 47, 1) ./ (255, 255, 255, 1),
    (216, 47, 99, 1) ./ (255, 255, 255, 1),
    (47, 216, 162, 1) ./ (255, 255, 255, 1),
    (216, 214, 47, 1) ./ (255, 255, 255, 1)
]
const cliffMinFlagHeight = 10
const cliffMaxFlagHeight = 10
const cliffMinFlagLength = 10
const cliffMaxFlagLength = 10
const cliffMinSpace = 2
const cliffMaxSpace = 8

const flaglineUnion = Union{Maple.ClothesLine, Maple.CliffFlags}

function renderFlagLine(ctx::Ahorn.Cairo.CairoContext, entity::flaglineUnion, lineColor::Ahorn.colorTupleType, pinColor::Ahorn.colorTupleType, colors::Array{Ahorn.colorTupleType, 1},
        minFlagHeight::Integer, maxFlagHeight::Integer, minFlagLength::Integer, maxFlagLength::Integer, minSpace::Integer, maxSpace::Integer)
    rng = Ahorn.getSimpleEntityRng(entity)
    x, y = Ahorn.position(entity)

    start = (x, y)
    stop = get(entity.data, "nodes", [start])[1]
    control = (start .+ stop) ./ 2 .+ (0, 10)

    startPoint = start[1] < stop[1] ? start : stop
    stopPoint = start[1] < stop[1] ? stop : start
    lastPoint = startPoint

    progress = 0
    curveLength = sqrt(sum((start .- stop) .^ 2))

    drawFlag = false

    baseCurve = Ahorn.SimpleCurve(startPoint, stopPoint, control)

    # Fix extra small line being wonky
    # Still looks bad for flagless longer lines
    if curveLength <= 16
        Ahorn.drawSimpleCurve(ctx, baseCurve, lineColor, thickness=1)

        return true
    end

    Ahorn.Cairo.save(ctx)

    Ahorn.set_antialias(ctx, 1)
    Ahorn.set_line_width(ctx, 1)
    Ahorn.fill_preserve(ctx)

    while progress < 1
        color = rand(rng, colors)
        hightlightColor = color .+ (0.1, 0.1, 0.1, 0)
        height = rand(rng, minFlagHeight:maxFlagHeight)
        length = rand(rng, minFlagLength:maxFlagLength)
        step = rand(rng, minSpace:maxSpace)

        progress += (drawFlag ? length : step) / curveLength

        point = Ahorn.getPoint(baseCurve, progress)

        Ahorn.move_to(ctx, lastPoint...)
        Ahorn.line_to(ctx, point...)
        Ahorn.setSourceColor(ctx, lineColor)
        Ahorn.stroke(ctx)

        if progress < 1 && drawFlag
            droop = length * droopAmount

            flagCurve = Ahorn.SimpleCurve(lastPoint, point, (lastPoint .+ point) ./ 2 .+ (0, droopAmount * length))
            prevFlagPoint = Ahorn.getPoint(flagCurve, 1 / length)

            Ahorn.setSourceColor(ctx, color)

            for i in 0:length
                flagPoint = Ahorn.getPoint(flagCurve, i / length)

                if flagPoint[1] > prevFlagPoint[1]
                    segmentWidth = min(floor(Int, flagPoint[1] - prevFlagPoint[1]), point[1] - flagPoint[1])

                    Ahorn.rectangle(ctx, flagPoint[1], flagPoint[2], segmentWidth, height - 1)
                end

                prevFlagPoint = flagPoint
            end

            Ahorn.stroke(ctx)

            Ahorn.setSourceColor(ctx, hightlightColor)
            Ahorn.rectangle(ctx, lastPoint[1], lastPoint[2], 0, height - 1)
            Ahorn.rectangle(ctx, point[1], point[2], 0, height - 1)
            Ahorn.stroke(ctx)

            Ahorn.setSourceColor(ctx, pinColor)
            Ahorn.rectangle(ctx, lastPoint[1], lastPoint[2] - 1, 0, 3)
            Ahorn.rectangle(ctx, point[1], point[2] - 1, 0, 3)
            Ahorn.stroke(ctx)
        end

        lastPoint = point
        drawFlag = !drawFlag
    end

    Ahorn.move_to(ctx, lastPoint...)
    Ahorn.line_to(ctx, stopPoint...)
    Ahorn.setSourceColor(ctx, lineColor)
    Ahorn.stroke(ctx)

    Ahorn.Cairo.restore(ctx)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ClothesLine)
    renderFlagLine(ctx, entity, selectedColor, clothesPinColor, clothesColors, clothesMinFlagHeight, clothesMaxFlagHeight,
            clothesMinFlagLength, clothesMaxFlagLength, clothesMinSpace, clothesMaxSpace)
end

function Ahorn.renderSelectedAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.CliffFlags)
    renderFlagLine(ctx, entity, selectedColor, cliffPinColor, cliffColors, cliffMinFlagHeight, cliffMaxFlagHeight,
            cliffMinFlagLength, cliffMaxFlagLength, cliffMinSpace, cliffMaxSpace)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.ClothesLine, room::Maple.Room)
    renderFlagLine(ctx, entity, clothesLineColor, clothesPinColor, clothesColors, clothesMinFlagHeight, clothesMaxFlagHeight,
            clothesMinFlagLength, clothesMaxFlagLength, clothesMinSpace, clothesMaxSpace)
end

function Ahorn.renderAbs(ctx::Ahorn.Cairo.CairoContext, entity::Maple.CliffFlags, room::Maple.Room)
    renderFlagLine(ctx, entity, cliffLineColor, cliffPinColor, cliffColors, cliffMinFlagHeight, cliffMaxFlagHeight,
            cliffMinFlagLength, cliffMaxFlagLength, cliffMinSpace, cliffMaxSpace)
end

end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
struct LCurve{T}  <: AutoRegMethod
    reg::T
end

function auto_reg(t,y,s,kernel::Function,method::LCurve;yoffset=true,tdomain=:realplus)
    creg = typeof(method.reg).name.wrapper

    if method.reg isa UnionAll
        λ_ini = 1
        creg = method.reg
    else
        λ_ini = method.reg.α
        creg = typeof(method.reg).name.wrapper
    end

    function obj(λ)
        reg=creg(λ)

        AR,yr = build_ar(t,y,s,kernel,reg,yoffset)
        #TODO: could reduce allocs , problem with forward diff
        # p=load!(o.problem, AR,yr)
        sy = get_yt(AR,yr;tdomain=tdomain)
        ρ = (AR*sy)[1:length(y)].-y |>norm
        η =  sy[1:end-1] |> norm
        return (ρ ,η)
    end
    #calculate curvatur
    function k(λ)
        dd=ForwardDiff.Dual{:a}(λ,one(1.))
        o=obj(dd)
        η =  o[2].value
        ρ =  o[1].value
        dη = o[2].partials[1]
        c= -2*η*ρ/dη*(λ^2*dη*ρ + 2*λ*η*ρ + λ^4*η*dη)/((λ^2*η^2+ ρ^2)^3/2 )
        @debug "L-curvature" κ=c λ=λ
       # return c
        #supress negative curvatures
        if c<0
            c=exp(c)
        else
            c=c+1
        end

        return c
    end

   # return k(λ_ini)
    #optimize on log scales
    λ_opt = 10. ^gradient_decent(x ->log(1/k(10.0 ^x)),ini=log10(λ_ini),d_ini=log10(λ_ini))
    @debug "found λ" λ=λ_opt
    return creg( λ_opt)
end

function gradient_decent(f;ini=1e-5,d_ini=ini,maxiters=40)
    k=f
    α_old =ini
    α_cur =ini+1
    α_new =ini+1
    dual=ForwardDiff.Dual{:c}(α_old,one(1.))
    kevel= k(dual)
    Δ_old = kevel.partials[1]
    dual=ForwardDiff.Dual{:c}(α_new,one(1.))
    kevel= k(dual)
    Δ_cur = kevel.partials[1]
    k_best= kevel.value
    k_cur=k_best
    α_opt=α_old

   for i  in 1:maxiters

        α_cur=α_new
        dual=ForwardDiff.Dual{:c}(α_cur,one(1.))
        kn= k(dual)

        Δ_cur = kn.partials[1]
        k_cur= kn.value

        # save best values
         if  k_cur< k_best
            k_best=k_cur
            α_opt=α_cur
        end
        γ=(abs(α_old-α_cur)*abs(Δ_cur-Δ_old))/((Δ_cur-Δ_old)^2+eps(Float64))#+1e-4
        α0=1
        α_new= clamp(α_cur-γ*Δ_cur,α_cur-α0,α_cur+α0)
        #α_new=clamp(α_cur-Δ_cur,α_cur-α0,α_cur+α0)
        α_old=α_cur
        Δ_old=Δ_cur
        i>1 &&γ<0.01 && abs(Δ_cur)<0.01 &&break
    end

    return α_opt
end

struct XuPei{T,S} <: AutoRegMethod
    reg::T
    λ::S  #shrinking parameter
end
XuPei(r::AbstractRegularization{T}) where {T}  = XuPei(r,0.9)

#from http://dx.doi.org/10.1016/j.flowmeasinst.2016.05.004
function auto_reg(t,y,s,kernel,method::XuPei;yoffset=true,tdomain=:realplus)

    creg = typeof(method.reg).name.wrapper

    λ =method.λ
    #define regularization vector
    if method.reg isa UnionAll
        Λ=ones(length(s)-1)*1e-3
        creg = method.reg
    else
        Λ=ones(length(s)-1)*method.reg.α
        creg = typeof(method.reg).name.wrapper
    end

    res = Inf
    for i in 1:100
        reg=creg(Λ)
        AR,yr = build_ar(t,y,s,kernel,reg,yoffset)
        sy= get_yt(AR,yr;tdomain=:real)
        sy[sy[1:end].<0][1:end-1].=0
        res_new= norm((AR*sy)[1:length(y)]-y)
        reg_new = norm((AR*sy)[length(y)+1:end-1])

        if  res < res_new
            break
        end

        res=res_new
        sy=sy[1:end-1]

        Λ[sy.>0] .=Λ[sy.>0]*λ
    end
    return creg(Λ)
end


text:
import Distributions

using Distributions, UncertainData

# Create an uncertain dataset containing both theoretical values with known parameters,
# theoretical values with fitted parameters and kernel density estimated distributions.
u1 = UncertainValue(Gamma, rand(Gamma(), 500))
u2 = UncertainValue(rand(MixtureModel([Normal(1, 0.3), Normal(-3, 3)]), 500))
uvals3 = [UncertainValue(Normal, rand(), rand()) for i = 1:11]
measurements = [u1; u2; uvals3]
d = UncertainDataset(measurements)

# Test all available constraints
constraints = [
    NoConstraint(),
    TruncateLowerQuantile(0.2),
    TruncateUpperQuantile(0.3),
    TruncateQuantiles(0.2, 0.6),
    TruncateMaximum(0.2),
    TruncateMinimum(0.2),
    TruncateRange(0.2, 0.6),
    TruncateStd(1)
]

for i = 1:length(constraints)

    # A single constraint applied to all values in the dataset
    @test constrain(d, constraints[i]) isa ConstrainedUncertainDataset

    # Element-wise application of different constraints to the values in the dataset
    @test constrain(d, [constraints[i] for k = 1:length(d)]) isa ConstrainedUncertainDataset

    # Constraining constrained datasets (might be nested several times)
    constrained_dataset = constrain(d, constraints[i])
    @test constrain(constrained_dataset, constraints[i]) isa ConstrainedUncertainDataset
    @test constrain(constrained_dataset, [constraints[i] for k = 1:length(d)]) isa ConstrainedUncertainDataset

end


text:
#static_unsigned = (SUInt128,SUInt16,SUInt32,SUInt64,SUInt8)

#static_signed = (SInt128,SInt16,SInt32,SInt64,SInt8)

#static_integer = (static_unsigned..., static_signed...)

const SNumber{V} = Union{SComplex{V},SReal{V}}

SNumber(x::BaseReal) = SReal(x)
SNumber(x::Complex) = SComplex(x)

function SNumber(val::Val{V}) where V
    if V isa Real
        if V isa Integer
            if V isa Unsigned
                if V isa UInt8
                    SUInt8{V}()
                elseif V isa UInt16
                    SUInt16{V}()
                elseif V isa UInt32
                    SUInt32{V}()
                elseif V isa UInt128
                    SUInt128{V}()
                else
                    SUInt64{V}()
                end
            elseif V isa Bool
                SBool{V}()
            elseif V isa Signed
                if V isa Int8
                    SInt8{V}()
                elseif V isa Int16
                    SInt16{V}()
                elseif V isa Int32
                    SInt32{V}()
                elseif V isa Int128
                    SInt128{V}()
                else
                    SInt64{V}()
                end
            end
        elseif V isa AbstractFloat
            if V isa Float16
                SFloat16{V}()
            elseif V isa Float32
                SFloat32{V}()
            else
                SFloat64{V}()
            end
        elseif V isa Rational
            SRational(val)
        end
    elseif V isa Complex
        SComplex(val)
    end
end

"""
    TPVal{H,L}
"""
struct TPVal{H,L}
    function TPVal{H,L}() where {H,L}
        if eltype(H) == eltype(L)
            return new{H,L}()
        else
            error("high and low precision values must be of the same type but got,
                  $(eltype(H)) and $(eltype(L)).")
        end
    end
end

const SVal{V} = Union{SNumber{V},SChar{V},TPVal{V}}

SVal(x::BaseNumber) = SNumber(x)
SVal(x::AbstractChar) = SChar(x)
SVal(x::TwicePrecision) = TPVal(x)

function SVal(val::Val{V}) where V
    if V isa Number
        SNumber(val)
    elseif V isa AbstractChar
        SChar(val)
    elseif V isa TwicePrecision
        TPVal(val)
    end
end

# static to base dict

for (SIT,BIT) in SI2BI
    BFT = float(BIT)
    SFT = BF2SF[BFT]
    @eval begin
        Base.AbstractFloat(::$SIT{X}) where X = $SFT{$BFT(X::$BIT)}()
    end
end

Base.AbstractFloat(::SBool{X}) where X = SFloat64{Float64(X::Bool)}()

for (SI,BI) in SSI2BSI
    BUI = unsigned(BI)
    SUI = BUI2SUI[BUI]
    @eval begin
        Base.unsigned(::$SI{X}) where X = $SUI{Base.bitcast($BUI, X::$BI)::$BUI}()
        Base.unsigned(::Type{<:$SI{X}}) where X = $SUI

        Base.signed(::$SUI{X}) where X = $SI{$BI(X::$BUI)::$BI}()
        Base.signed(::Type{<:$SUI{X}}) where X = $SI
    end
end

#Base.signed(x::Unsigned) = reinterpret(typeof(convert(Signed, zero(x))), x

#Base.signed(x::Unsigned) = reinterpret(typeof(convert(Signed, zero(x))), x

for (ST1,BT1) in S2B
    for (ST2,BT2) in S2B
        if BT1 != BT2
            BASERULE = promote_rule(BT1, BT2)
            if BASERULE != Union{}
                NEWRULE = B2S[BASERULE]
                @eval begin
                    promote_rule(::Type{<:$ST1{<:Any}}, ::Type{$BT2}) = $BASERULE
                    promote_rule(::Type{<:$ST2{<:Any}}, ::Type{$BT1}) = $BASERULE
                    promote_rule(::Type{<:$ST1{<:Any}}, ::Type{<:$ST2}) = $NEWRULE
                end
            end
        end
    end
end

promote_rule(::Type{<:AbstractIrrational}, ::Type{SFloat16}) = SFloat16
promote_rule(::Type{<:AbstractIrrational}, ::Type{SFloat32}) = SFloat32
promote_rule(::Type{<:AbstractIrrational}, ::Type{T}) where {T<:SReal} = promote_type(SFloat64, T)

for (ST1,BT1) in S2B
    for (ST2,BT2) in S2B
        if BT1 != BT2
            @eval begin
                @generated function (::Type{<:$ST2{<:Any}})(::$ST1{X}) where X
                    ret = $ST2{$BT2(X::$BT1)}()
                    :($ret)
                end
                (::Type{$BT2})(::$ST1{X}) where X = $BT2(X::$BT1)
                ofeltype(::Type{$BT1}, val::$ST2{V}) where V = $ST1{$BT1(V::$BT2)}()

                ofeltype(::Type{$BT1}, val::$BT2) = $BT1(val)
                ofeltype(::$BT1, val::$ST2{V}) where V = $ST1{$BT1(V::$BT2)}()
                ofeltype(::$BT1, val::$BT2) = $BT1(val)

                ofeltype(::Type{<:$ST1}, val::$ST2{V}) where V = $ST1{$BT1(V::$BT2)}()
                ofeltype(::Type{<:$ST1}, val::$BT2) = $BT1(val)
                ofeltype(::$ST1, val::$ST2{V}) where V = $ST1{$BT1(V::$BT2)}()
                ofeltype(::$ST1, val::$BT2) = $BT1(val)

                Base.flipsign(::$ST1{V1}, ::$ST2{V2}) where {V1,V2} = $ST1{flipsign(V1::$BT1,V2::$BT2)}()
                Base.copysign(::$ST1{V1}, ::$ST2{V2}) where {V1,V2} = $ST1{copysign(V1::$BT1,V2::$BT2)}()
            end
        end
    end
end


text:
include("perlin.jl")

include("gtk-compose.jl")

max_stars_repo_name: lobingera/ColorMaps


text:
using Symbolics, SymbolicUtils
using SymbolicUtils: istree, operation, arguments
using Symbolics: value, get_variables, solve_for, derivative
using SymbolicUtils.Rewriters
using SymbolicUtils.Code

include("poly.jl")
include("arith.jl")
include("roots.jl")
include("factors.jl")
include("tests.jl")


text:
using MathProgComplex, OPFInstances
using Base.Test

brokentestset = Set([])

testcases = OrderedSet([("WB2", 1),
                        ("LMBM3", 1),
                        ("WB5", 1),
                        ("case6ww", 1),
                        ("case9mod", 1),
                        ("case9", 1),
                        ("case14", 1)])
                        # ("case24_ieee_rts", 1),
                        # ("case_ieee30", 1),
                        # ("case30pwl", 1),
                        # ("case30", 1)])

@testset "SDP hierarchy ACOPFs relaxation - rank 1" begin

    @testset "$pbsolved tests" for pbsolved in [:SOSRelaxation, :MomentRelaxation]
        @testset "$instance problem - rank $d" for (instance, d) in testcases
            problem_c, point = import_from_dat(getinstancepath("Matpower", "QCQP", instance))
            problem = pb_cplx2real(problem_c)

            relax_ctx = set_relaxation(problem; hierarchykind=:Real,
                                                d = d,
                                                params = Dict(:opt_outlev=>0,
                                                                :opt_solver=>testsolver,
                                                                :opt_pbsolved=>pbsolved))

            primobj, dualobj = run_hierarchy(problem, relax_ctx);

            ε_rel = OPFsols[(instance, 1)].rel_opt_gap
            ε_abs = max(primobj, dualobj) * ε_rel

            @show testsolver, pbsolved, instance, d
            @show primobj, dualobj,  OPFsols[(instance, 1)].primal_solvalue, ε_abs, instance

            @test primobj ≈ OPFsols[(instance, 1)].primal_solvalue atol=ε_abs
            @test dualobj ≈ primobj atol=ε_abs
        end
    end
end


text:
# closed loop moving obstacle avoidance when goal is within lidar range
function goal_in_lidar()
  c=defineCase(;(:mode=>:autoBench));
  setMisc!(c;mpc_max_iter=300,tex=0.5,max_cpu_time=0.466,Nck=[10,8,6]);
  n=initializeAutonomousControl(c);
  n.s.evalConstraints=false;
  driveStraight!(n)
  for ii=2:n.mpc.max_iter
     if ((n.r.dfs_plant[end][:x][end]-c.g.x_ref)^2 + (n.r.dfs_plant[end][:y][end]-c.g.y_ref)^2)^0.5 < 1.8*n.XF_tol[1]
        n.mpc.goal_reached=true; break;
      end
      updateAutoParams!(n,c);                      # update model parameters
      status=autonomousControl!(n);                # rerun optimization
      n.mpc.t0_actual=(n.r.eval_num-1)*n.mpc.tex;  # external so that it can be updated easily in PathFollowing
      simPlant!(n)
  end
  return n.mpc.goal_reached
end


text:
files =
  [
    joinpath("BivariateNormal", "AM", "function"),
    joinpath("BivariateNormal", "AM", "pdf"),
    joinpath("BivariateNormal", "MALA", "function", "analytical"),
    joinpath("BivariateNormal", "MALA", "function", "forwarddiff"),
    joinpath("BivariateNormal", "MALA", "function", "reversediff"),
    joinpath("BivariateNormal", "MALA", "pdf", "analytical"),
    joinpath("BivariateNormal", "MALA", "pdf", "forwarddiff"),
    joinpath("BivariateNormal", "MALA", "pdf", "reversediff"),
    joinpath("BivariateNormal", "SMMALA", "analytical"),
    joinpath("BivariateNormal", "SMMALA", "forwarddiff"),
    joinpath("BivariateNormal", "SMMALA", "reversediff"),
    joinpath("BivariateNormal", "Gibbs"),
    joinpath("Gamma", "MH", "transformation"),
    joinpath("Gamma", "MH", "truncation"),
    joinpath("Gamma", "MALA"),
    joinpath("Normal", "AM", "function"),
    joinpath("Normal", "AM", "pdf"),
    joinpath("Normal", "HMC", "function", "analytical"),
    joinpath("Normal", "HMC", "function", "forwarddiff"),
    joinpath("Normal", "HMC", "pdf", "analytical"),
    joinpath("Normal", "HMC", "pdf", "forwarddiff"),
    joinpath("Normal", "MALA", "function", "analytical"),
    joinpath("Normal", "MALA", "function", "forwarddiff"),
    joinpath("Normal", "MALA", "pdf", "analytical"),
    joinpath("Normal", "MALA", "pdf", "forwarddiff"),
    joinpath("Normal", "NUTS", "function", "dualaveraging", "analytical"),
    joinpath("Normal", "NUTS", "function", "dualaveraging", "forwarddiff"),
    joinpath("Normal", "NUTS", "function", "noadaptation", "analytical"),
    joinpath("Normal", "NUTS", "function", "noadaptation", "forwarddiff"),
    joinpath("Normal", "SliceSampler"),
    joinpath("Poisson", "MH"),
    joinpath("swiss", "HMC", "dualaveraging", "analytical"),
    joinpath("swiss", "HMC", "dualaveraging", "forwarddiff"),
    joinpath("swiss", "HMC", "dualaveraging", "reversediff"),
    joinpath("swiss", "HMC", "noadaptation", "analytical"),
    joinpath("swiss", "HMC", "noadaptation", "forwarddiff"),
    joinpath("swiss", "HMC", "noadaptation", "reversediff"),
    joinpath("swiss", "NUTS", "dualaveraging", "analytical"),
    joinpath("swiss", "NUTS", "dualaveraging", "forwarddiff"),
    joinpath("swiss", "NUTS", "dualaveraging", "reversediff"),
    joinpath("swiss", "NUTS", "noadaptation", "analytical"),
    joinpath("swiss", "NUTS", "noadaptation", "forwarddiff"),
    joinpath("swiss", "NUTS", "noadaptation", "reversediff"),
    joinpath("swiss", "SMMALA", "analytical"),
    joinpath("swiss", "SMMALA", "forwarddiff"),
    joinpath("swiss", "SMMALA", "reversediff"),
    joinpath("swiss", "RAM"),
    joinpath("swiss", "SliceSampler")
  ]

println("Running examples:")

for file in files
  filejl = file*".jl"
  println("  * $(filejl) *")
  include(filejl)
end


text:
function isprime(n::Int64)
    if n <= 3
        return n > 1
    elseif n % 2 == 0 || n % 3 == 0
        return false
    end
    i = 5
    while i^2 <= n
        if n % i == 0 || n % (i+2) == 0
            return false
        end
        i += 6
    end
    return true
end

function get_smallest_divisor(n::Int64)
    k = 2
    while n % k ≠ 0
        k += 1
    end
    return k
end

function get_prime_factors(N::Int64)
    S = Vector{Int64}()
    while true
        !isprime(N) || break
        k = get_smallest_divisor(N)
        push!(S, k)
        N::Int64 = N/k
    end
    push!(S, N)
    return S
end

function reverse_int(n::Int64)
    if n <= 9
        return n
    else
        return parse(Int, reduce(string, digits(n)))
    end
end

function totient(n::Int64)
    result::Float64 = n
    p = 2
    while p*p <= n
        if n % p == 0
            while n % p == 0
                n = n//p
            end
            result *= (1.0 - (1.0 / p))
        end
        p += 1
    end

    if n > 1
        result *= (1.0 - (1.0 / n))
    end

    return result
end

function prime_sieve(N::Int64)
    iscoprime(P, i) = !any(x -> i % x == 0, P)

    primes = Int64[]
    for num = 2:N
        if iscoprime(primes, num)
            push!(primes, num)
        end
    end
    primes
end


text:
using Test
using EllZeroTrendFiltering: unsafe_minimum

# Subtraction of quadratic polynomials
@test QuadraticPolynomial(4.0, -2.0, 3.0) - QuadraticPolynomial(1.0, 1.0, 2.0) == QuadraticPolynomial(3.0,-3.0,1.0)

### Minimization of quadratic polynomials ###
@test find_minimum(QuadraticPolynomial(1 ,0 ,0)) == (0, 0)
@test find_minimum(QuadraticPolynomial(1 ,0 ,1)) == (0, 1)

# 2*x*(x+1) + 1
@test find_minimum(QuadraticPolynomial(2, 2, 1)) == (-0.5, 0.5)

# (x+1)^2 + 1
@test find_minimum(QuadraticPolynomial(1, 2, 2)) == (-1, 1)

# 2*((x+1)^2 + 1)
@test find_minimum(QuadraticPolynomial(2, 4, 4)) == (-1, 2)

# (3x+1)^2 + 1 = 9x^2 + 6x + 2
@test find_minimum(QuadraticPolynomial(9, 6, 2)) == (-1/3, 1)

# (3x+1)^2 + 1 = 9x^2 + 6x + 2
@test unsafe_minimum(QuadraticPolynomial(9.0, 6.0, 2.0)) == 1.0

@test zero(QuadraticPolynomial{Float64}) == QuadraticPolynomial(0.0, 0.0, 0.0)

@test zero(QuadraticPolynomial(4.0, -2.0, 3.0)) == QuadraticPolynomial(0.0, 0.0, 0.0)

@test zero(QuadraticPolynomial(4, -2, 3)) == QuadraticPolynomial(0, 0, 0)


text:
# These functions compute the sizes of combinations of parameters.

using Combinatorics: combinations

"""
    total_combinations(arity, n_way)

Given an array of the number of values for each parameter and the
level of coverage, return the number of n-tuples required for complete
coverage.
"""
function total_combinations(arity, n_way)
  param_cnt = length(arity)
  sum(prod(arity[key_set]) for key_set in combinations(1:param_cnt, n_way))
end

"""
    next_multiplicative!(values, arity)

Given a set of values for parameters, this returns the next possible
value. When it reaches the end, it cycles back to the start. The first
value is all ones. The last value equals the arity.
"""
function next_multiplicative!(values, arity)
    carry = 1
    for slot_idx in length(values):-1:1
        values[slot_idx] += carry
        if arity[slot_idx] < values[slot_idx]
            values[slot_idx] = 1
            carry = 1
        else
            carry = 0
        end
    end
end

"""
    all_combinations(arity, n_way)

This represents possible coverage as a matrix, one row per parameter,
zero if not used. `arity` is a list of the number of values for each parameter,
and `n_way` is the order of the combinations, most commonly 2-way.
"""
function all_combinations(arity, n_way)
    v_cnt = length(arity)
    # This returns a list of lists, so it has length, not size.
    indices = collect(combinations(1:v_cnt, n_way))
    combinations_cnt = total_combinations(arity, n_way)

    coverage = zeros(eltype(arity), v_cnt, combinations_cnt)
    idx = 1
    for indices_idx in 1:length(indices)
        offset = indices[indices_idx]
        sub_arity = arity[offset]
        sub_cnt = prod(sub_arity)
        values = copy(sub_arity)
        for sub_idx in 1:sub_cnt
            next_multiplicative!(values, sub_arity)
            coverage[offset, idx] = values
            idx += 1
        end
    end
    coverage
end

"""
    one_parameter_combinations(arity, n_way)

Generates all combinations that are nonzero for the last parameter.
This is for in-parameter-order generation, where we need
only those tuples that end with this column being nonzero.

The construction method is to leave out the given parameter
and construct all `n_way` - 1 tuples. Then copy and paste that
once for each possible value of the given parameter.
"""
function one_parameter_combinations(arity, n_way)
    param_cnt = length(arity)
    if n_way > 1
        partial = all_combinations(arity[1:(param_cnt-1)], n_way - 1)
        one_set = size(partial, 2)
        comb = zeros(eltype(arity), param_cnt, one_set * arity[param_cnt])
        for vidx in 1:(arity[param_cnt])
            col_begin = (vidx - 1) * one_set + 1
            col_end = vidx * one_set
            comb[1:size(partial, 1), col_begin:col_end] .= partial
            comb[param_cnt, col_begin:col_end] .= vidx
        end
    else  # n_way == 1
        comb = zeros(eltype(arity), param_cnt, arity[param_cnt])
        comb[param_cnt, :] .= 1:(arity[param_cnt])
    end
    comb
end

combination_number(n, m) = prod(n:-1:(n-m+1)) ÷ factorial(m)

function pairs_in_entry(entry, n_way)
    n = length(entry)
    ans = zeros(Int, combination_number(n, n_way), n)
    col_set_idx = 1
    for param_idx in combinations(1:length(entry), n_way)
        for row_idx in 1:n_way
            ans[param_idx[row_idx], col_set_idx] = entry[param_idx[row_idx]]
        end
        col_set_idx += 1
    end
    ans
end


text:
using AlgebraPDF
using Test
using Random
using Measurements

@testset "iminuit should be imported and have :Minuit" begin
    @test isiminuitimported()
    @test hasproperty(AlgebraPDF.iminuit, :Minuit)
end

@testset "Fit gauss with Minuit" begin
    d0 = Normalized(FGauss((μ=1.1, σ=0.3)), (-3, 5))
    Random.seed!(1234)
    data = randn(1000)
    fs0 = fit(d0, data, MigradAndHesse())

    @test keys(fs0.parameters) == (:μ,:σ)
    @test typeof(Tuple(fs0.measurements)) <: Tuple{Vararg{Measurement}}

    d1 = fixpar(Normalized(FGauss(Ext(μ=1.1, σ=0.9)), (-3, 5)), :σ)
    fs1 = fit(d1, data, MigradAndHesse())

    @test keys(fs1.parameters) == (:μ,)
    @test pars(fs1.best_model).σ == 0.9
end

@testset "Extended Likelihood fit" begin

    d = Normalized(FGauss((μ=1.1,σ=0.1)), (-2, 3))
    Random.seed!(4321)
    data = filter(x->inrange(x,lims(d)), randn(1000) .+ 0.5)

    fr = fit(d, data)
    @show fr.parameters.μ
    @show abs(fr.parameters.σ)
    @test 0.3 < fr.parameters.μ < 0.7
    @test 0.7 < abs(fr.parameters.σ) < 1.3

    s = FSum([d],(α=2.2,))

    @test integral(s) == 2.2
    @test s(1.1) == 2.2*d(1.1)
    @test s([1,2,3]) == 2.2*d([1,2,3])

    enll = Extended(NegativeLogLikelihood(s, data))
    fr2 = fit(enll)
    @test 0.3 < fr2.parameters.μ < 0.7
    @test 0.7 < abs(fr2.parameters.σ) < 1.3
    @test 900 < fr2.parameters.α < 1100
end
#
# Optim often fails with `isfinite(phi_c) && isfinite(dphi_c)`
# removed from the tests
#
@testset "Fit gauss with Optim" begin
    d0 = Normalized(FGauss((μ=0.01, σ=0.9)), (-3, 5))
    Random.seed!(1234)
    data = randn(1000)
    fs0 = fit(d0, data, BFGSApproxHesse())

    @test keys(fs0.parameters) == (:μ,:σ)
    @test typeof(Tuple(fs0.measurements)) <: Tuple{Vararg{Measurement}}

    d1 = fixpar(Normalized(FGauss(Ext(μ=0.01, σ=0.9)), (-3, 5)), :σ)
    fs1 = fit(d1, data, BFGSApproxHesse())

    @test keys(fs1.parameters) == (:μ,)
    @test pars(fs1.best_model).σ == 0.9
end