text:
struct StreamerInfo
    streamer
    dependencies
end

struct Streamers
    tkey::TKey
    refs::Dict{Int32, Any}
    elements::Vector{StreamerInfo}
end

Base.length(s::Streamers) = length(s.elements)

function Base.show(io::IO, s::Streamers)
    for streamer_info in s.elements
        println(io, "$(streamer_info.streamer.fName)")
        # streamer = streamer_info.streamer
        # print(io, "$(streamer.fName): fType = $(streamer.fType), ")
        # print(io, "fTypeName: $(streamer.fTypeName)")
    end
end

# Structures required to read streamers
struct TStreamerInfo
    fName
    fTitle
    fCheckSum
    fClassVersion
    fElements
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TStreamerInfo})
    preamble = Preamble(io, T)
    fName, fTitle = nametitle(io)
    fCheckSum = readtype(io, UInt32)
    fClassVersion = readtype(io, Int32)
    fElements = readobjany!(io, tkey, refs)
    endcheck(io, preamble)
    T(fName, fTitle, fCheckSum, fClassVersion, fElements)
end

safename(s::AbstractString) = replace(s, "::" => "_3a3a_")

function initialise_streamer(s::StreamerInfo)
    # FIXME Abstract is not needed when switched to autogenerated streamers
    base = Symbol(safename(s.streamer.fName))
    supername = Symbol(:Abstract, base)
    if !isdefined(@__MODULE__, supername)
        @debug "Defining abstract type $supername"
        @eval abstract type $(supername) <: ROOTStreamedObject end
    end

    name = Symbol(base, Symbol("_$(s.streamer.fClassVersion)"))
    if !isdefined(@__MODULE__, name)
        @debug "  creating versioned struct '$name <: $supername'"
        @eval struct $(name) <: $(supername) end
        # FIXME create the stream!() functions somewhere here...
        # println(name)
        # @eval struct $name <: ROOTStreamedObject
        #     data::Dict{Symbol, Any}
        # end
    else
        @debug "Not defining $name since it already has a bootstrapped version"
    end
end

"""
    function Streamers(io)

Reads all the streamers from the ROOT source.
"""
function Streamers(io)
    refs = Dict{Int32, Any}()

    start = position(io)
    tkey = unpack(io, TKey)

    if iscompressed(tkey)
        @debug "Compressed stream at $(start)"
        seekstart(io, tkey)
        compression_header = unpack(io, CompressionHeader)
        #FIXME for some reason we need to re-pack such that it ends at exact bytes.
        skipped = position(io) - start
        # notice our `TKey` size is not the same as official TKey, can't use sizeof()
        io_buf = IOBuffer(read(io, tkey.fNbytes - skipped))

        if String(compression_header.algo) == "ZL"
            stream = IOBuffer(read(ZlibDecompressorStream(io_buf), tkey.fObjlen))
        elseif String(compression_header.algo) == "XZ"
            stream = IOBuffer(read(XzDecompressorStream(io_buf), tkey.fObjlen))
        elseif String(compression_header.algo) == "L4"
            skip(io_buf, 8) #skip checksum
            stream = IOBuffer(lz4_decompress(read(io_buf), tkey.fObjlen))
        else
            error("Unsupported compression type '$(String(compression_header.algo))'")
        end

    else
        @debug "Unompressed stream at $(start)"
        stream = io
    end
    preamble = Preamble(stream, Streamers)
    skiptobj(stream)

    name = readtype(stream, String)
    size = readtype(stream, Int32)
    streamer_infos = Vector{StreamerInfo}()
    @debug "Found $size streamers, continue with parsing."
    for i ∈ 1:size
        obj = readobjany!(stream, tkey, refs)
        if typeof(obj) == TStreamerInfo
            @debug "  processing streamer info for '$(obj.fName)' (v$(obj.fClassVersion))"
            @debug "    number of dependencies: $(length(obj.fElements.elements))"
            dependencies = Set()
            for element in obj.fElements.elements
                if typeof(element) == TStreamerBase
                    @debug "      + adding dependency '$(element.fName)'"
                    push!(dependencies, element.fName)
                else
                    @debug "      - skipping dependency '$(element.fName)' with type '$(typeof(element))'"
                end
            end
            @debug "      => finishing dependency readout for: $(obj.fName)"
            push!(streamer_infos, StreamerInfo(obj, dependencies))
        else
            @debug "  not a TStreamerInfo but '$(typeof(obj))', skipping."
        end
        # FIXME why not just skip a byte?
        skip(stream, readtype(stream, UInt8))
    end

    endcheck(stream, preamble)

    streamer_infos = topological_sort(streamer_infos)

    for streamer_info in streamer_infos
        initialise_streamer(streamer_info)
    end

    Streamers(tkey, refs, streamer_infos)
end

"""
    function topological_sort(streamer_infos)

Sort the streamers with respect to their dependencies and keep only those
which are not defined already.

The implementation is based on https://stackoverflow.com/a/11564769/1623645
"""
function topological_sort(streamer_infos)
    @debug "Starting topological sort of streamers"
    provided = Set{String}()
    sorted_streamer_infos = []
    while length(streamer_infos) > 0
        remaining_items = []
        emitted = false
        @debug "  number of remaining streamers to sort: $(length(streamer_infos))"

        for streamer_info in streamer_infos
            # if all(d -> isdefined(@__MODULE__, Symbol(d)) || d ∈ provided, streamer_info.dependencies)
            #     if !isdefined(@__MODULE__, Symbol(streamer_info.streamer.fName)) && aliasfor(streamer_info.streamer.fName) === nothing
            @debug "    processing '$(streamer_info.streamer.fName)' with $(length(streamer_info.dependencies))' dependencies"
            if length(streamer_infos) ==  1 || all(d -> d ∈ provided, streamer_info.dependencies)
                if aliasfor(streamer_info.streamer.fName) === nothing
                    push!(sorted_streamer_infos, streamer_info)
                end
                push!(provided, streamer_info.streamer.fName)
                emitted = true
            else
                push!(remaining_items, streamer_info)
            end
        end

        if !emitted
            for streamer_info in streamer_infos
                filter!(isequal(streamer_info), remaining_items)
            end
        end

        streamer_infos = remaining_items
    end
    @debug "Finished the topological sort of streamers"
    sorted_streamer_infos
end

"""
    function readobjany!(io, tkey::TKey, refs)

The main entrypoint where streamers are parsed and cached for later use.
The `refs` dictionary holds the streamers or parsed data which are reused
when already available.
"""
function readobjany!(io, tkey::TKey, refs)
    beg = position(io) - origin(tkey)
    bcnt = readtype(io, UInt32)
    if Int64(bcnt) & Const.kByteCountMask == 0 || Int64(bcnt) == Const.kNewClassTag
        # New class or 0 bytes
        version = 0
        start = 0
        tag = bcnt
        bcnt = 0
    else
        version = 1
        start = position(io) - origin(tkey)
        tag = readtype(io, UInt32)
    end

    if Int64(tag) & Const.kClassMask == 0
        # reference object
        if tag == 0
            return missing
        elseif tag == 1
            error("Returning parent is not implemented yet")
        elseif !haskey(refs, tag)
            # skipping
            seek(io, origin(tkey) + beg + bcnt + 4)
            return missing
        else
            return refs[tag]
        end

    elseif tag == Const.kNewClassTag
        cname = readtype(io, CString)
        streamer = getfield(@__MODULE__, Symbol(cname))

        if version > 0
            refs[start + Const.kMapOffset] = streamer
        else
            refs[length(refs) + 1] = streamer
        end

        obj = unpack(io, tkey, refs, streamer)

        if version > 0
            refs[beg + Const.kMapOffset] = obj
        else
            refs[length(refs) + 1] = obj
        end

        return obj
    else
        # reference class, new object
        ref = Int64(tag) & ~Const.kClassMask
        haskey(refs, ref) || error("Invalid class reference.")

        streamer = refs[ref]
        obj = unpack(io, tkey, refs, streamer)

        if version > 0
            refs[beg + Const.kMapOffset] = obj
        else
            refs[length(refs) + 1] = obj
        end

        return obj
    end
end

struct TList
    preamble
    name
    size
    objects
end

Base.length(l::TList) = length(l.objects)

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TList})
    preamble = Preamble(io, T)
    skiptobj(io)

    name = readtype(io, String)
    size = readtype(io, Int32)
    objects = []
    for i ∈ 1:size
        push!(objects, readobjany!(io, tkey, refs))
        skip(io, readtype(io, UInt8))
    end

    endcheck(io, preamble)
    TList(preamble, name, size, objects)
end

struct TObjArray
    name
    low
    elements
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TObjArray})
    preamble = Preamble(io, T)
    skiptobj(io)
    name = readtype(io, String)
    size = readtype(io, Int32)
    low = readtype(io, Int32)
    elements = Vector{Any}(undef, size)
    for i in 1:size
        ele = readobjany!(io, tkey, refs)
        # !ismissing(ele) && @show ele.fName
        elements[i] = ele
     end
    endcheck(io, preamble)
    return TObjArray(name, low, elements)
end

abstract type AbstractTStreamerElement end

@premix @with_kw mutable struct TStreamerElementTemplate
    version
    fOffset
    fName
    fTitle
    fType
    fSize
    fArrayLength
    fArrayDim
    fMaxIndex
    fTypeName
    fXmin
    fXmax
    fFactor
end

@TStreamerElementTemplate mutable struct TStreamerElement end

@pour initparse begin
    fields = Dict{Symbol, Any}()
end

function parsefields!(io, fields, T::Type{TStreamerElement})
    preamble = Preamble(io, T)
    fields[:version] = preamble.version
    fields[:fOffset] = 0
    fields[:fName], fields[:fTitle] = nametitle(io)
    fields[:fType] = readtype(io, Int32)
    fields[:fSize] = readtype(io, Int32)
    fields[:fArrayLength] = readtype(io, Int32)
    fields[:fArrayDim] = readtype(io, Int32)

    n = preamble.version == 1 ? readtype(io, Int32) : 5
    fields[:fMaxIndex] = [readtype(io, Int32) for _ in 1:n]

    fields[:fTypeName] = readtype(io, String)

    if fields[:fType] == 11 && (fields[:fTypeName] == "Bool_t" || fields[:fTypeName] == "bool")
        fields[:fType] = 18
    end

    fields[:fXmin] = 0.0
    fields[:fXmax] = 0.0
    fields[:fFactor] = 0.0

    if preamble.version == 3
        fields[:fXmin] = readtype(io, Float64)
        fields[:fXmax] = readtype(io, Float64)
        fields[:fFactor] = readtype(io, Float64)
    end

    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TStreamerElement})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

@TStreamerElementTemplate mutable struct TStreamerBase
    fBaseVersion
end

function parsefields!(io, fields, T::Type{TStreamerBase})
    preamble = Preamble(io, T)
    parsefields!(io, fields, TStreamerElement)
    fields[:fBaseVersion] = fields[:version] >= 2 ? readtype(io, Int32) : 0
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TStreamerBase})
    @initparse
    parsefields!(io, fields, T)
    T(;fields...)
end

@TStreamerElementTemplate mutable struct TStreamerBasicType end

function parsefields!(io, fields, ::Type{TStreamerBasicType})
    parsefields!(io, fields, TStreamerElement)

    if Const.kOffsetL < fields[:fType] < Const.kOffsetP
        fields[:fType] -= Const.kOffsetP
    end
    basic = true
    if fields[:fType] ∈ (Const.kBool, Const.kUChar, Const.kChar)
        fields[:fSize] = 1
    elseif fields[:fType] in (Const.kUShort, Const.kShort)
        fields[:fSize] = 2
    elseif fields[:fType] in (Const.kBits, Const.kUInt, Const.kInt, Const.kCounter)
        fields[:fSize] = 4
    elseif fields[:fType] in (Const.kULong, Const.kULong64, Const.kLong, Const.kLong64)
        fields[:fSize] = 8
    elseif fields[:fType] in (Const.kFloat, Const.kFloat16)
        fields[:fSize] = 4
    elseif fields[:fType] in (Const.kDouble, Const.kDouble32)
        fields[:fSize] = 8
    elseif fields[:fType] == Const.kCharStar
        fields[:fSize] = sizeof(Int)
    else
        basic = false
    end

    if basic && fields[:fArrayLength] > 0
        fields[:fSize] *= fields[:fArrayLength]
    end

end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TStreamerBasicType})
    @initparse
    preamble = Preamble(io, T)
    parsefields!(io, fields, T)
    endcheck(io, preamble)
    T(;fields...)
end

@TStreamerElementTemplate mutable struct TStreamerBasicPointer
    fCountVersion
    fCountName
    fCountClass
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TStreamerBasicPointer})
    @initparse
    preamble = Preamble(io, T)
    parsefields!(io, fields, TStreamerElement)
    fields[:fCountVersion] = readtype(io, Int32)
    fields[:fCountName] = readtype(io, String)
    fields[:fCountClass] = readtype(io, String)
    endcheck(io, preamble)
    T(;fields...)
end

@TStreamerElementTemplate mutable struct TStreamerLoop
    fCountVersion
    fCountName
    fCountClass
end

components(::Type{TStreamerLoop}) = [TStreamerElement]

function parsefields!(io, fields, ::Type{TStreamerLoop})
    fields[:fCountVersion] = readtype(io, Int32)
    fields[:fCountName] = readtype(io, String)
    fields[:fCountClass] = readtype(io, String)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TStreamerLoop})
    @initparse
    preamble = Preamble(io, T)
    for component in components(T)
        parsefields!(io, fields, component)
    end
    parsefields!(io, fields, T)
    endcheck(io, preamble)
    T(;fields...)
end

abstract type AbstractTStreamSTL end

@TStreamerElementTemplate mutable struct TStreamerSTL <: AbstractTStreamSTL
    fSTLtype
    fCtype
end

@TStreamerElementTemplate mutable struct TStreamerSTLstring <: AbstractTStreamSTL
    fSTLtype
    fCtype
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, ::Type{T}) where T <: AbstractTStreamSTL
    @initparse
    if T == TStreamerSTLstring
        wrapper_preamble = Preamble(io, T)
    end
    preamble = Preamble(io, T)
    parsefields!(io, fields, TStreamerElement)

    fields[:fSTLtype] = readtype(io, Int32)
    fields[:fCtype] = readtype(io, Int32)

    if fields[:fSTLtype] == Const.kSTLmultimap || fields[:fSTLtype] == Const.kSTLset
        if startswith(fields[:fTypeName], "std::set") || startswith(fields[:fTypeName], "set")
            fields[:fSTLtype] = Const.kSTLset
        elseif startswith(fields[:fTypeName], "std::multimap") || startswith(fields[:fTypeName], "multimap")
            fields[:fSTLtype] = Const.kSTLmultimap
        end
    end

    endcheck(io, preamble)
    if T == TStreamerSTLstring
        endcheck(io, wrapper_preamble)
    end
    T(;fields...)
end

const TObjString = String

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{TObjString})
    preamble = Preamble(io, T)
    skiptobj(io)
    value = readtype(io, String)
    endcheck(io, preamble)
    T(value)
end

abstract type AbstractTStreamerObject end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, ::Type{T}) where T<:AbstractTStreamerObject
    @initparse
    preamble = Preamble(io, T)
    parsefields!(io, fields, TStreamerElement)
    endcheck(io, preamble)
    T(;fields...)
end

@TStreamerElementTemplate mutable struct TStreamerObject <: AbstractTStreamerObject end
@TStreamerElementTemplate mutable struct TStreamerObjectAny <: AbstractTStreamerObject end
@TStreamerElementTemplate mutable struct TStreamerObjectAnyPointer <: AbstractTStreamerObject end
@TStreamerElementTemplate mutable struct TStreamerObjectPointer <: AbstractTStreamerObject end
@TStreamerElementTemplate mutable struct TStreamerString <: AbstractTStreamerObject end

abstract type ROOTStreamedObject end

# function stream(io, ::Type{T}) where {T<:ROOTStreamedObject}
#     fields = Dict{Symbol, Any}()
#     preamble = Preamble(io, T)
#     stream!(io, fields, T{preamble.version})
#     endcheck(io, preamble)
#     T(fields)
# end

function stream!(io, fields, ::Type{T}) where {T<:ROOTStreamedObject}
    preamble = Preamble(io, T)
    streamer_name = Symbol(T, "_$(preamble.version)")
    # @show streamer_name
    mod, typename = split(String(streamer_name), ".")
    # @show mod typename
    streamer = getfield(@__MODULE__, Symbol(typename))
    # @show streamer
    readfields!(io, fields, streamer)
    endcheck(io, preamble)
end

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, ::Type{T}) where {T<:ROOTStreamedObject}
    cursor = Cursor(position(io), io, tkey, refs)
    @initparse
    preamble = Preamble(io, T)
    streamer_name = Symbol(T, "_$(preamble.version)")
    mod, typename = split(String(streamer_name), ".")
    streamer = getfield(@__MODULE__, Symbol(typename))
    readfields!(cursor, fields, streamer)
    streamer(;cursor=cursor, fields...)
end

# function stream!(io, fields, ::Type{T{V}}) where {V, T<:ROOTStreamedObject}
#     println("Don't know how to stream $T")
# end

struct TObject <: ROOTStreamedObject end
parsefields!(io, fields, ::Type{TObject}) = skiptobj(io)

struct TString <: ROOTStreamedObject end
unpack(io, tkey::TKey, refs::Dict{Int32, Any}, ::Type{TString}) = readtype(io, String)

struct Undefined <: ROOTStreamedObject
    skipped_bytes
end

Base.show(io::IO, u::Undefined) = print(io, "$(typeof(u)) ($(u.skipped_bytes) bytes)")

function unpack(io, tkey::TKey, refs::Dict{Int32, Any}, T::Type{Undefined})
    preamble = Preamble(io, T)
    bytes_to_skip = preamble.cnt - 6
    skip(io, bytes_to_skip)
    endcheck(io, preamble)
    Undefined(bytes_to_skip)
end

const TArrayD = Vector{Float64}
const TArrayI = Vector{Int32}

function readtype(io, T::Type{Vector{U}}) where U <: Union{Integer, AbstractFloat}
    size = readtype(io, eltype(T))
    [readtype(io, eltype(T)) for _ in 1:size]
end


text:
using Lindenmayer, Test

function run_all_tests()
simple                = LSystem(Dict("F" => "F[t+FoF-F]"), "F")
koch                  = LSystem(Dict("F" => "F+F-"), "F") # use turn of 60 degrees
koch                  = LSystem(Dict("F" => "F+F--F+F"), "F") # 4 # 60 degrees
koch_snowflake        = LSystem(Dict("F" => "F+F--F+F"), "F-F-F")
peano                 = LSystem(Dict("F" => "TF+F-F-toF-F+F+F+F-F"), "3F")
peano_gosper          = LSystem(Dict("X" => "X+YF++YF-tFX--FXFX-YF+", "Y" => "-FX+YFYF++YF+FX--FX-Y" ), "FX" )  # turn 60 degree
quadratic_Koch        = LSystem(Dict("F" => "tF-F+F+FFF-F-F+F"), "4F+F+F+F")
thirty_two_segment    = LSystem(Dict("F" => "t-F+F-F-F+F+FF-F+F+FF+F-F-FF+FF-FF+F+F-FF-F-F+FF-F-F+F+F-F+"), "F+F+F+F")
sierpinski_triangle   = LSystem(Dict("F" => "G+F+Gt", "G"=>"F-G-F"), "G") #  6  # 60 Degree
square_curve          = LSystem(Dict("X"=>"XF-F+F-XF+F+XtF-F+F-X"), "F+XF+F+XF")
dragon_curve          = LSystem(Dict("F"=>"F+G+t", "G"=>"-F-G"), "F") #  90 degrees
hilbert               = LSystem(Dict("L" => "+RF-LFL-cFR+", "R" => "-LF+RFR+FL-"), "1L")
hilbert_curve         = LSystem(Dict("L" => "+RF-LFL-tFR+", "R" => "-LF+RFR+FL-"), "3L") # 90deg
hilbert_curve2        = LSystem(Dict("X" => "XFYFX+F+YFXFcY-F-XFYFX", "Y" => "YFXFY-F-XFYFX+F+YFXFY"), "2X")
plant                 = LSystem(Dict("F" => "F[-F]cF[+F][F]"), "F") # use turn eg 17, 23
plant1                = LSystem(Dict("F" => "FF", "X" => "F−[[cX]+X]+F[+FX]−X"), "1X")
branch                = LSystem(Dict("F" => "FF-[F+F+Fc]+[+F-F-F][+++F+F-F---][---F+F-F---]"), "1FFFF")
penrose               = LSystem(Dict("X"  =>  "PM++QM----YM[-PM----XM]++t",
                                      "Y"  => "+PM--QM[---XM--YM]+t",
                                      "P"  => "-XM++YM[+++PM++QM]-t",
                                      "Q"  => "--PM++++XM[+QM++++YM]--YMt",
                                      "M"  => "F",
                                      "F"  => ""),
                                      "1[Y]++[Y]++[Y]++[Y]++[Y]") # 36 degrees

drawLSystem(simple, forward=50, turn=90, iterations=6, startingx=0, startingy=-150, filename="simple.pdf", showpreview=false)

drawLSystem(koch, forward=5, turn=60, iterations=6, startingx=-1800, startingy= 0, width=4000, height=4000, filename="koch.pdf", showpreview=false)

drawLSystem(quadratic_Koch, forward=10, iterations=2, turn=90, filename="quadratic_koch.pdf", showpreview=false)

drawLSystem(hilbert, forward=12, turn=90, iterations=6, startingx=-450, startingy= -450, filename="hilbert.pdf", showpreview=false)

drawLSystem(hilbert_curve, forward=25, turn=90, iterations=4, startingx=-200, startingy=-200, filename="hilbert_curve.pdf", showpreview=false)

drawLSystem(hilbert_curve2, forward=10, turn=90, iterations=4, startingx=-450, startingy= -450, filename="hilbert_curve2.pdf", showpreview=false)

drawLSystem(dragon_curve, forward=12, turn=90, iterations=10, filename="dragon_curve.pdf", showpreview=false)

drawLSystem(peano, turn=90, forward=20, iterations=3, startingx=-250, filename="peano.pdf", showpreview=false)

drawLSystem(peano_gosper, forward=10, turn=60, iterations=4, startingpen=(0, 0.8, 0.2), startingorientation=-pi/2, startingx=-200, startingy= -50, filename="peano-gosper.pdf", showpreview=false)

drawLSystem(thirty_two_segment, forward=5, iterations=2, turn=90, filename="32segments.pdf", showpreview=false)

drawLSystem(sierpinski_triangle, forward=3, startingx=-400, startingy= -350, turn=60, iterations=8, filename="sierpinski-triangle.pdf", showpreview=false)

drawLSystem(plant, forward=7, startingpen=(0, 0.8, 0.3), startingx=0, startingy= 460, startingorientation=-pi/2, turn=23, iterations=6, filename="plant.pdf", showpreview=false)

drawLSystem(plant1, forward=3, turn=13, iterations=7, startingpen=(0, 0.8, 0.2), startingorientation=-pi/2, startingx=-50, startingy= 450, filename="plant1.pdf", showpreview=false)

drawLSystem(branch, forward=12, turn=20, iterations=3, startingpen=(0, 0.9, 0.2), startingorientation=-pi/2, startingx=0, startingy= 300, filename="branch.pdf", showpreview=false)

drawLSystem(penrose, forward=25, turn=36, iterations=7, startingpen=(.5, 0.8, 0.2), startingorientation=-pi/2, startingx=0, startingy= 0, filename="penrose.pdf", showpreview=false)

end

if get(ENV, "LINDENMAYER_KEEP_TEST_RESULTS", false) == "true"
        # they changed mktempdir in v1.3
        if VERSION <= v"1.2"
            cd(mktempdir())
        else
            cd(mktempdir(cleanup=false))
        end
        @info("...Keeping the Lindenmayer test output in: $(pwd())")
        run_all_tests()
        @info("Lindenmayer test images were saved in: $(pwd())")
else
    mktempdir() do tmpdir
        cd(tmpdir) do
            @info("running Lindenmayer tests in: $(pwd())")
            @info("but not keeping the results")
            @info("because you didn't do: ENV[\"LINDENMAYER_KEEP_TEST_RESULTS\"] = \"true\"")
            run_all_tests()
            @info("Test images weren't saved. To see the test images, next time do this before running:")
            @info(" ENV[\"LINDENMAYER_KEEP_TEST_RESULTS\"] = \"true\"")
        end
    end
end

ext: jl

max_stars_repo_name: cormullion/Lindenmayer

max_issues_repo_name: cormullion/Lindenmayer

max_forks_repo_name: cormullion/Lindenmayer


text:
using CSDP

using Test
using LinearAlgebra # For Diagonal

@testset "Interact with BLAS" begin
    vec = Cdouble[1.0, 2.0, 0.0, -1.0]
    l = length(vec)
    inc = 1
    n1 = ccall((BLAS.@blasfunc(dasum_), LinearAlgebra.BLAS.libblas),
               Cdouble,
               (Ref{LinearAlgebra.BlasInt}, Ptr{Cdouble},
                Ref{LinearAlgebra.BlasInt}),
               l, vec, inc)
    @test abs(n1 - 4) < 1e-15
end

@testset "Call libcsdp.norm1" begin
    vec = Cdouble[1.0, 2.0, 0.0, -1.0]
    try
        n1 = ccall((:norm1, CSDP.libcsdp),
                   Cdouble,
                   (Cint, Ptr{Cdouble}),
                   length(vec), vec)
        @test abs(n1 - 4) < 1e-15
    catch
        println("n1 = $n1, vec=$vec, length(vec)=$(length(vec))")
        rethrow()
    end
end

@static if !Sys.iswindows() # FIXME Segfault on Windows
    @testset "Example" begin
        cd(joinpath(dirname(dirname(pathof(CSDP))), "examples")) do
            include(joinpath(pwd(), "example.jl"))
            @test size(X) == (7, 7)
            @test length(y) == 2
            @test size(Z) == (7, 7)
            X✓ = [3 3  0 0 0 0 0;
                  3 3  0 0 0 0 0;
                  0 0 16 0 0 0 0;
                  0 0  0 0 0 0 0;
                  0 0  0 0 0 0 0;
                  0 0  0 0 0 0 0;
                  0 0  0 0 0 0 0] / 24
            @test norm(convert(AbstractArray, X) - X✓) < 1e-6
            y✓ = [3, 4] / 4
            @test norm(convert(AbstractArray, y) - y✓) < 1e-6
            Z✓ = [1 -1 0 0 0 0 0;
                 -1  1 0 0 0 0 0;
                  0  0 0 0 0 0 0;
                  0  0 0 8 0 0 0;
                  0  0 0 0 8 0 0;
                  0  0 0 0 0 3 0;
                  0  0 0 0 0 0 4] / 4
            @test norm(convert(AbstractArray, Z) - Z✓) < 1e-6
        end
    end
end

@testset "Options" begin
    @test CSDP.paramstruc(Dict(:axtol => 1e-7)).axtol == 1e-7
end

@testset "MathOptInterface" begin
    include("MOI_wrapper.jl")
end


text:
module Day19v2

export day19

using CombinedParsers
using CombinedParsers.Regexp: whitespace_maybe

struct Input
    rules::Dict{Int,Union{AbstractString,Vector{Vector{Int}}}}
    messages::Vector{String}
end

function parserule(str)
    numbers = join(Repeat(Numeric(Int)), ' ')
    indirect = join(numbers, " | ")
    direct = Either("a", "b")
    q = Either(indirect, direct)
    p = Sequence(:no => Numeric(Int), ": ", :rules => q )
    str |> p
end

function parseinput(filename)
    rules = Dict{Int,Union{AbstractString,Vector{Vector{Int}}}}()
    messages = Vector{String}()
    mode = :rule
    for line in readlines(filename)
        if line == ""
            mode = :message
            continue
        end
        if mode == :rule
            result = parserule(line)
            if result[:rules] isa AbstractString
                rules[result[:no]] = result[:rules]
            else
                rules[result[:no]] = result[:rules]
            end
        else
            push!(messages, strip(line))
        end
    end
    Input(rules, messages)
end

get_re(r::AbstractString, rules) = r
get_re(r::Int, rules) = get_re(rules[r], rules)

function get_re(list, rules)
    rr = Either{Any}(Any[])
    for rule in list
        inner = [get_re(r, rules) for r in rule]
        push!(rr, Sequence(inner...))
    end
    rr
end

function day19(::Val{:a}, filename)
    input = parseinput(filename)
    re = AtStart() * get_re(input.rules[0], input.rules) * AtEnd()
    count = 0
    for message in input.messages
        if !isnothing(match(re, message))
            count += 1
        end
    end
    count
end

function day19(::Val{:b}, filename)
    input = parseinput(filename)
    re42 = get_re(input.rules[42], input.rules)
    re31 = get_re(input.rules[31], input.rules)
    count = 0
    for message in input.messages
        for i in 2:9, j in 1:i-1
                rr = AtStart() * Sequence([re42 for k in 1:i]...) * Sequence([re31 for l in 1:j]...) * AtEnd()
                if !isnothing(match(rr, message))
                    count += 1
                    break
                end
        end
    end
    count
end

end


text:
n = rand(1:10)

A = rand(n,n)

r = polarfact(A, alg =:hybrid);

# Test unitary matrix U

U = r.U
H = r.H

@test_approx_eq_eps U'*U eye(n) 1e-7

# Test Hermitian positive semifefinite matrix H

@test issym(H)

for i in eigvals(H)
    @test i >= 0.
end

@test_approx_eq_eps A U*H 1e-7

println("Hybrid method passed test...")


text:
"""
    DispersalKernel <: AbstractKernelNeighborhood

    DispersalKernel(; kw...)

Dispersal kernel for taking the dot product of the neighborhood and a matching
kernel of weights. May hold any `Neighborhood` object: the kernel
will be built to match the shape, using the `folumation`, `cellsize`
and `distancemethod`.

# Keyword Arguments

- `neighborhood`: Any DynamicGrids.jl `Neighborhood`, or an
    already constructed [`DispersalKernel`](@ref). Using this keyword means `radius` is
    ignored, and for a `DispersalKernel`, all other keywords are ignored.
- `neighborhood`: `Neighborhood` object specifying the range from the origin of the
    discretised dispersal kernal. Defaults to `Window(radius)`.
- `formulation`: kernel formulation object holding the exact form of the kernal.
    Default [`ExponentialKernel`](@ref).
- `cellsize`: the cell size of the discretised kernal (i.e. simulation grid size).
    Default is 1.0.
- `distancemethod`: [`DistanceMethod`](@ref) object for calculating distance between cells.
    The default is [`CentroidToCentroid`](@ref).
"""
struct DispersalKernel{R,L,N<:Neighborhood{R,L},K,F,C,D<:DistanceMethod} <: AbstractKernelNeighborhood{R,L}
    neighborhood::N
    kernel::K
    formulation::F
    cellsize::C
    distancemethod::D
    function DispersalKernel(
        hood::N, kernel, formulation::F, cellsize::C, distancemethod::D
    ) where {N<:Neighborhood{R,L},F,C,D<:DistanceMethod} where {R,L}
        if hood isa DispersalKernel
            hood
        else
            # Build the kernel matrix
            newkernel = scale(buildkernel(hood, formulation, distancemethod, cellsize))
            new{R,L,N,typeof(newkernel),F,C,D}(
                hood, newkernel, formulation, cellsize, distancemethod
            )
        end
    end
    function DispersalKernel{R,L,N,K,F,C,D}(
        hood::N, kernel::K, formulation::F, cellsize::C, distancemethod::D
    ) where {R,L,N,K,F,C,D}
        new{R,L,N,K,F,C,D}(hood, kernel, formulation, cellsize, distancemethod)
    end
end
function DispersalKernel(;
    radius=1,
    neighborhood=Window{radius}(),
    formulation=ExponentialKernel(),
    cellsize=1.0,
    distancemethod=CentroidToCentroid(),
)
    DispersalKernel(neighborhood, nothing, formulation, cellsize, distancemethod)
end
DispersalKernel{R}(; radius=R, kw...) where R = DispersalKernel(; radius=radius, kw...)

ConstructionBase.constructorof(::Type{<:DispersalKernel}) = DispersalKernel

function DG._setbuffer(n::DispersalKernel{R,L,<:Any,K,F,C,D}, buffer) where {R,L,K,F,C,D}
    newhood = DG._setbuffer(neighborhood(n), buffer)
    DispersalKernel{R,L,typeof(newhood),K,F,C,D}(
        newhood, kernel(n), formulation(n), cellsize(n), distancemethod(n)
    )
end

cellsize(hood::DispersalKernel) = hood.cellsize
distancemethod(hood::DispersalKernel) = hood.distancemethod
formulation(hood::DispersalKernel) = hood.formulation

function buildkernel(window::Window{R}, formulation, distancemethod, cellsize) where R
    # The radius doesn't include the center cell, so add it
    S = 2R + 1
    kernel = zeros(typeof(cellsize), S, S)
    r1 = R + one(R)
    # Paper: l. 97
    for x = 0:R, y = 0:R
        # Calculate the distance effect from the center cell to this cell
        prob = dispersalprob(formulation, distancemethod, x, y, cellsize)
        # Update the kernel value based on the formulation and distance
        kernel[ x + r1,  y + r1] = prob
        kernel[ x + r1, -y + r1] = prob
        kernel[-x + r1,  y + r1] = prob
        kernel[-x + r1, -y + r1] = prob
    end
    SMatrix{S,S}(kernel)
end
function buildkernel(window::Neighborhood{<:Any,L}, f, dm, cellsize) where L
    SVector{L}(Tuple(dispersalprob(f, dm, x, y, cellsize) for (x, y) in offsets(window)))
end

scale(x) = x ./ sum(x)


text:
mutable struct ImportExportData{GID <: Integer, PID <: Integer, LID <: Integer}
    source::BlockMap{GID, PID, LID}
    target::BlockMap{GID, PID, LID}

    permuteToLIDs::Array{LID, 1}
    permuteFromLIDs::Array{LID, 1}
    remoteLIDs::Array{LID, 1}

    exportLIDs::Array{LID, 1}
    exportPIDs::Array{PID, 1}

    numSameIDs::GID
    distributor::Distributor{GID, PID, LID}

    isLocallyComplete::Bool
end

## Constructors ##
function ImportExportData(source::BlockMap{GID, PID, LID}, target::BlockMap{GID, PID, LID})::ImportExportData{GID, PID, LID} where GID <: Integer where PID <:Integer where LID <: Integer
    ImportExportData{GID, PID, LID}(source, target, [], [], [], [], [], 0, createDistributor(getComm(source)), true)
end

## Getters ##
function sourceMap(data::ImportExportData{GID, PID, LID})::BlockMap{GID, PID, LID} where GID <: Integer where PID <:Integer where LID <: Integer
    data.source
end

function targetMap(data::ImportExportData{GID, PID, LID})::BlockMap{GID, PID, LID} where GID <: Integer where PID <:Integer where LID <: Integer
    data.target
end

function permuteToLIDs(data::ImportExportData{GID, PID, LID})::Array{LID} where GID <: Integer where PID <:Integer where LID <: Integer
    data.permuteToLIDs
end

function permuteFromLIDs(data::ImportExportData{GID, PID, LID})::Array{LID} where GID <: Integer where PID <:Integer where LID <: Integer
    data.permuteFromLIDs
end

function remoteLIDs(data::ImportExportData{GID, PID, LID})::Array{LID} where GID <: Integer where PID <: Integer where LID <: Integer
    data.remoteLIDs
end

function remoteLIDs(data::ImportExportData{GID, PID, LID}, remoteLIDs::AbstractArray{<: Integer}) where GID <: Integer where PID <: Integer where LID <: Integer
    data.remoteLIDs = remoteLIDs
end

function exportLIDs(data::ImportExportData{GID, PID, LID})::Array{LID} where GID <: Integer where PID <: Integer where LID <: Integer
    data.exportLIDs
end

function exportLIDs(data::ImportExportData{GID, PID, LID}, exportLIDs::AbstractArray{<: Integer}) where GID <: Integer where PID <: Integer where LID <: Integer
    data.exportLIDs = exportLIDs
end

function exportPIDs(data::ImportExportData{GID, PID, LID})::Array{PID} where GID <: Integer where PID <: Integer where LID <: Integer
    data.exportPIDs
end

function exportPIDs(data::ImportExportData{GID, PID, LID}, exportPIDs::AbstractArray{<: Integer}) where GID <: Integer where PID <: Integer where LID <: Integer
    data.exportPIDs = exportPIDs
end

function numSameIDs(data::ImportExportData{GID, PID, LID})::LID where GID <: Integer where PID <: Integer where LID <: Integer
    data.numSameIDs
end

function numSameIDs(data::ImportExportData{GID, PID, LID}, numSame::Integer)::LID where GID <: Integer where PID <: Integer where LID <: Integer
    data.numSameIDs = numSame
end

function distributor(data::ImportExportData{GID, PID, LID})::Distributor{GID, PID, LID} where GID <: Integer where PID <: Integer where LID <: Integer
    data.distributor
end

function isLocallyComplete(data::ImportExportData{GID, PID, LID})::Bool where GID <: Integer where PID <: Integer where LID <: Integer
    data.isLocallyComplete
end

function isLocallyComplete(data::ImportExportData{GID, PID, LID}, isLocallyComplete::Bool) where GID <: Integer where PID <: Integer where LID <: Integer
    data.isLocallyComplete = isLocallyComplete
end


text:
using Test
using DelimitedFiles, LinearAlgebra
import MLJBase
const Mlj = MLJBase
using StableRNGs
rng = StableRNG(123)
using BetaML.Perceptron

println("*** Additional testing for the Perceptron algorithms...")

println("Testing MLJ interface for Perceptron models....")

X, y                      = Mlj.@load_iris

model                     = PerceptronClassifier()
regressor                 = Mlj.machine(model, X, y)
Mlj.evaluate!(regressor, resampling=Mlj.CV(), measure=Mlj.LogLoss())

model                     = KernelPerceptronClassifier()
regressor                 = Mlj.machine(model, X, y)
Mlj.evaluate!(regressor, resampling=Mlj.CV(), measure=Mlj.LogLoss())

model                     = PegasosClassifier()
regressor                 = Mlj.machine(model, X, y)
Mlj.evaluate!(regressor, resampling=Mlj.CV(), measure=Mlj.LogLoss())


text:
module PlotReferenceImages

local_path(args...) = normpath(@__DIR__, "..", args...)

"""
    reference_file(backend::Symbol, i::Int, version::String)

Find the latest version of the reference image file for the reference image `i` and the backend `be`.
This returns a path to the file in the folder of the latest version.
If no file is found, a path pointing to the file of the folder specified by `version` is returned.
"""
function reference_file(backend, i, version)
    refdir = local_path("Plots", string(backend))
    fn = "ref$i.png"
    versions = sort(VersionNumber.(readdir(refdir)), rev = true)

    reffn = joinpath(refdir, string(version), fn)
    for v in versions
        tmpfn = joinpath(refdir, string(v), fn)
        if isfile(tmpfn)
            reffn = tmpfn
            break
        end
    end

    return reffn
end

reference_path(backend, version) = local_path("Plots", string(backend), string(version))

export reference_file, reference_path

end # module


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct Group <: Edit
        file::String
        from::Symbol
        to::Array{Symbol,1}
        input::Symbol
        output::Array{Symbol,1}
    end

Use to edit files containing data in successive dataframes with an identifying header cell or row.

# Arguments
- `file::String`: mapping .csv file name in the coremaps directory. The mapping file should correlate with the header information identifying each data group. It will be used to separate the header rows from data.
- `from::Symbol`: name of the mapping column containing input values
- `to::Array{Symbol,1}`: name of the mapping column containing output values
- `input::Symbol`: name of the input column containing
- `output::Array{Symbol,1}`: name of the output column created
"""
mutable struct Group <: Edit
    "mapping .csv file name in the coremaps directory. The mapping file should correlate with the header information identifying each data group. It will be used to separate the header rows from data."
    file::String
    "name of the mapping column containing input values"
    from::Symbol
    "name of the mapping column containing output values"
    to::Array{Symbol,1}
    "name of the input column containing"
    input::Symbol
    "name of the output column created"
    output::Array{Symbol,1}
end

function Group(; file, from, to, input, output, )
    Group(file, from, to, input, output, )
end

"""Get [`Group`](@ref) `file`."""
get_file(value::Group) = value.file
"""Get [`Group`](@ref) `from`."""
get_from(value::Group) = value.from
"""Get [`Group`](@ref) `to`."""
get_to(value::Group) = value.to
"""Get [`Group`](@ref) `input`."""
get_input(value::Group) = value.input
"""Get [`Group`](@ref) `output`."""
get_output(value::Group) = value.output

"""Set [`Group`](@ref) `file`."""
set_file!(value::Group, val) = value.file = val
"""Set [`Group`](@ref) `from`."""
set_from!(value::Group, val) = value.from = val
"""Set [`Group`](@ref) `to`."""
set_to!(value::Group, val) = value.to = val
"""Set [`Group`](@ref) `input`."""
set_input!(value::Group, val) = value.input = val
"""Set [`Group`](@ref) `output`."""
set_output!(value::Group, val) = value.output = val

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
function TAME_LRT_increasing_cliques(testing=true)

    if testing

        n_sizes = [25]
        orders_to_test = [4]
        samples = 10^3
        trials = 2
        betas=[0.0]
        alphas=[.5]
        results_path = all_results_path*"RG_DupNoise/test/"

    else

        n_sizes = [100]
        orders_to_test = [3,4,5,6,7]
        samples = 10^4
        trials = 25
        betas=[0.0,1.0]
        alphas=[.5,1.0]
        results_path = all_results_path*"RG_DupNoise/"

    end
    step_percentage = [.2]
    edge_inclusion_p = .5

    if !isdir(results_path)
        mkpath(results_path)
    end

    LambdaTAME_results = []
    for order in orders_to_test
        X = distributed_random_trials(trials,DuplicationNoise(),true;postProcessing=noPostProcessing(),
                                        n_sizes=n_sizes,edge_inclusion_p,step_percentage,
                                        method=ΛTAME_MultiMotif_M(),graph=RandomGeometric(),degreedist=LogNormal(log(5),1),
                                        betas,profile=true,orders=[order],samples=samples)
        push!(LambdaTAME_results,(order,X))
    end
    filename=filter(x-> !isspace(x),"LambdaTAME_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_orders:$(orders_to_test)_p:$(edge_inclusion_p)_sample:$(samples)_sp:$(step_percentage)_trials:$(trials)_results.json")
    save_results(results_path*filename,LambdaTAME_results)

    LowRankTAME_results = []
    for order in orders_to_test
        X = distributed_random_trials(trials,DuplicationNoise(),true;postProcessing=noPostProcessing(),
                                      n_sizes=n_sizes,edge_inclusion_p,step_percentage,method=TAME_MultiMotif_M(),
                                      graph=RandomGeometric(),degreedist=LogNormal(log(5),1),betas,profile=true,
                                      orders=[order],samples=samples)
        push!(LowRankTAME_results,(order,X))
    end
    filename=filter(x-> !isspace(x),"LowRankTAME_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_orders:$(orders_to_test)_p:$(edge_inclusion_p)_sample:$(samples)_sp:$(step_percentage)_trials:$(trials)_results.json")
    save_results(results_path*filename,LowRankTAME_results)

    TAME_results = []
    for order in orders_to_test
        X = distributed_random_trials(trials,DuplicationNoise(),true;postProcessing=noPostProcessing(),
                                      n_sizes=n_sizes,edge_inclusion_p,step_percentage,
                                      method=LowRankTAME_MultiMotif_M(),graph=RandomGeometric(),degreedist=LogNormal(log(5),1),
                                      betas,profile=true,orders=[order],samples=samples)
        push!(TAME_results,(order,X))
    end
    filename=filter(x-> !isspace(x),"TAME_RandomGeometric_degreedist:LogNormal_alpha:$(alphas)_beta:$(betas)_n:$(n_sizes)_orders:$(orders_to_test)_p:$(edge_inclusion_p)_sample:$(samples)_sp:$(step_percentage)_trials:$(trials)_results.json")
    save_results(results_path*filename,TAME_results)

end

function increased_clique_size_duplication_noise_experiments(testing=true)

    if testing
        #  --  test params  --  #
        n_sizes = [500]
        step_percentage = [.1]
        edge_inclusion_p = [.5]
        trials = 1
        samples = 1000
        maxiter = 100
        betas = [1.0]
        alphas = [.5]
        orders_to_test = [3]
        k_vals = [5]
        save = true
        result_path = all_results_path*"RG_DupNoise/test/"

    else
        edge_inclusion_p = [.5]
        n_sizes = [500]
        step_percentage = [.25]
        trials = 25
        samples = 1000000
        maxiter = 1000

        orders_to_test = [2,3,4,5,6,7,8,9]
        #orders_to_test = [3,4,5]
        k_vals = [15,30,45,60,75,90]
        betas = [0.0,1.0,10.0,100.0]
        save = true
        result_path = all_results_path*"RG_DupNoise/"

    end

    if !isdir(result_path)
        mkpath(result_path)
    end

    LambdaTAME_Klau_results = []
    for order in orders_to_test
        for k in k_vals

            X = distributed_random_trials(trials,DuplicationNoise(),true;postProcessing=KlauAlgo(k=k,maxiter=maxiter,rtype=1),n_sizes=n_sizes,edge_inclusion_p=edge_inclusion_p,step_percentage=step_percentage,
                                                            method=ΛTAME_MultiMotif_M(),graph=RandomGeometric(),degreedist=LogNormal(log(5),1),motif=Clique(),matchingMethod=ΛTAME_GramMatching(),
                                                            orders=[order],samples=samples,betas,profile=true)
            push!(LambdaTAME_Klau_results,(order,k,X))
        end
    end
    filename = filter(x-> !isspace(x),"RandomGeometric_degreedist:LogNormal_KlauAlgokvals:$(k_vals)_noiseModel:Duplication_n:$(n_sizes)_KAmiter:$(maxiter)_orders:$(orders_to_test)_p:$(edge_inclusion_p)_samples:$(samples)_sp:$(step_percentage)_trials:$(trials)_data.json")
    save_results(result_path*filename,LambdaTAME_Klau_results)

    LambdaTAME_LocalSearch_results = []
    for order in orders_to_test
        for k in k_vals

            X = distributed_random_trials(trials,DuplicationNoise(),true;postProcessing=LocalSearch(k=k),n_sizes=n_sizes,edge_inclusion_p=edge_inclusion_p, step_percentage=step_percentage,
                                                            method=ΛTAME_MultiMotif_M(),graph=RandomGeometric(),degreedist=LogNormal(log(5),1),motif=Clique(),matchingMethod=ΛTAME_GramMatching(),
                                                            orders=[order],samples=samples,betas,profile=true)
            push!(LambdaTAME_LocalSearch_results,(order,k,X))
        end
    end
    filename = filter(x-> !isspace(x),"RandomGeometric_degreedist:LogNormal_LocalSearchkvals:$(k_vals)_noiseModel:Duplication_n:$(n_sizes)_orders:$(orders_to_test)_p:$(edge_inclusion_p)_samples:$(samples)_sp:$(step_percentage)_trials:$(trials)_data.json")
    save_results(result_path*filename,LambdaTAME_LocalSearch_results)

end

max_stars_repo_name: chuckcol/LambdaTAME

max_issues_repo_name: chuckcol/LambdaTAME

max_forks_repo_name: chuckcol/LambdaTAME


text:
mutable struct AdaptiveQRData{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}}
    data::CachedMatrix{T,DM,M}
    τ::Vector{T}
    ncols::Int
end

function AdaptiveQRData(::AbstractBandedLayout, A::AbstractMatrix{T}) where T
    l,u = bandwidths(A)
    data = BandedMatrix{T}(undef,(2l+u+1,0),(l,l+u)) # pad super
    AdaptiveQRData(CachedArray(data,A), Vector{T}(), 0)
end
AdaptiveQRData(A::AbstractMatrix{T}) where T = AdaptiveQRData(MemoryLayout(typeof(A)), A)

function partialqr!(F::AdaptiveQRData{<:Any,<:BandedMatrix}, n::Int)
    if n > F.ncols
        l,u = bandwidths(F.data.array)
        resizedata!(F.data,n+l,n+u+l);
        resize!(F.τ,n);
        ñ = F.ncols
        factors = view(F.data.data,ñ+1:n+l,ñ+1:n);
        τ = view(F.τ,ñ+1:n);
        _banded_qr!(factors, τ);
        # multiply remaining columns
        n̄ = max(ñ+1,n-l-u+1) # first column interacting with extra data
        Q = QRPackedQ(view(F.data.data,n̄:n+l,n̄:n), view(F.τ,n̄:n))
        lmul!(Q',view(F.data.data,n̄:n+l,n+1:n+u+l))
        F.ncols = n
    end
    F
end

struct AdaptiveQRFactors{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: AbstractMatrix{T}
    data::AdaptiveQRData{T,DM,M}
end

struct AdaptiveLayout{M} <: MemoryLayout end
MemoryLayout(::Type{AdaptiveQRFactors{T,DM,M}}) where {T,DM,M} = AdaptiveLayout{typeof(MemoryLayout(DM))}()
triangularlayout(::Type{Tri}, ::ML) where {Tri, ML<:AdaptiveLayout} = Tri{ML}()

size(F::AdaptiveQRFactors) = size(F.data.data)
bandwidths(F::AdaptiveQRFactors) = bandwidths(F.data.data)
function colsupport(F::AdaptiveQRFactors, j)
    partialqr!(F.data, j)
    colsupport(F.data.data.data, j)
end

function rowsupport(F::AdaptiveQRFactors, j)
    partialqr!(F.data, j+bandwidth(F,2))
    rowsupport(F.data.data.data, j)
end

function getindex(F::AdaptiveQRFactors, k::Int, j::Int)
    partialqr!(F.data, j)
    F.data.data[k,j]
end

colsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = colsupport(F.factors, j)
rowsupport(F::QRPackedQ{<:Any,<:AdaptiveQRFactors}, j) = rowsupport(F.factors, j)

Base.replace_in_print_matrix(A::AdaptiveQRFactors, i::Integer, j::Integer, s::AbstractString) =
    i in colsupport(A,j) ? s : Base.replace_with_centered_mark(s)
Base.replace_in_print_matrix(A::UpperTriangular{<:Any,<:AdaptiveQRFactors}, i::Integer, j::Integer, s::AbstractString) =
    i in colsupport(A,j) ? s : Base.replace_with_centered_mark(s)

struct AdaptiveQRTau{T,DM<:AbstractMatrix{T},M<:AbstractMatrix{T}} <: AbstractVector{T}
    data::AdaptiveQRData{T,DM,M}
end

size(F::AdaptiveQRTau) = (size(F.data.data,1),)

function getindex(F::AdaptiveQRTau, j::Int)
    partialqr!(F.data, j)
    F.data.τ[j]
end

getR(Q::QR, ::Tuple{OneToInf{Int},OneToInf{Int}}) where T = UpperTriangular(Q.factors)

function _banded_qr(::NTuple{2,OneToInf{Int}}, A)
    data = AdaptiveQRData(A)
    QR(AdaptiveQRFactors(data), AdaptiveQRTau(data))
end

#########
# lmul!
#########

function lmul!(A::QRPackedQ{<:Any,<:AdaptiveQRFactors}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
    require_one_based_indexing(B)
    mA, nA = size(A.factors)
    sB = B.datasize[1]
    mB = length(B)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    Afactors = A.factors

    resizedata!(B, maximum(colsupport(A,sB)))
    b = B.data

    sB_2 = length(b)
    partialqr!(A.factors.data, sB)

    @inbounds begin
        for k = sB:-1:1
            cs = colsupport(Afactors,k)
            vBj = b[k]
            for i = (k+1:sB_2) ∩ cs
                vBj += conj(Afactors[i,k])*b[i]
            end
            vBj = A.τ[k]*vBj
            b[k] -= vBj
            for i = (k+1:sB_2) ∩ cs
                b[i] -= Afactors[i,k]*vBj
            end
        end
    end
    B
end

function lmul!(adjA::Adjoint{<:Any,<:QRPackedQ{<:Any,<:AdaptiveQRFactors}}, B::CachedVector{T,Vector{T},<:Zeros{T,1}}) where T
    COLGROWTH = 1000 # rate to grow columns
    tol = floatmin(real(T))

    require_one_based_indexing(B)
    A = adjA.parent
    mA, nA = size(A.factors)
    mB = length(B)
    if mA != mB
        throw(DimensionMismatch("matrix A has dimensions ($mA,$nA) but B has dimensions ($mB, $nB)"))
    end
    sB = length(B.data)
    jr = 1:min(COLGROWTH,nA)

    @inbounds begin
        while first(jr) < nA
            cs = colsupport(A.factors, last(jr))
            cs_max = maximum(cs)
            kr = first(jr):cs_max
            resizedata!(B, min(cs_max,mB))
            if (first(jr) > sB && maximum(abs,view(B.data,colsupport(A.factors,first(jr)))) ≤ tol)
                break
            end
            partialqr!(A.factors.data, min(cs_max,nA))
            Q_N = QRPackedQ(view(A.factors.data.data.data,kr,jr), view(A.τ.data.τ,jr))
            lmul!(Q_N', view(B.data, kr))
            jr = last(jr)+1:min(last(jr)+COLGROWTH,nA)
        end
    end
    B
end

function (*)(A::QRPackedQ{T,<:AdaptiveQRFactors}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
end

function (*)(A::Adjoint{T,<:QRPackedQ{T,<:AdaptiveQRFactors}}, x::AbstractVector{S}) where {T,S}
    TS = promote_op(matprod, T, S)
    lmul!(A, Base.copymutable(convert(AbstractVector{TS},x)))
end

function ldiv!(R::UpperTriangular{<:Any,<:AdaptiveQRFactors}, B::CachedVector{<:Any,<:Any,<:Zeros{<:Any,1}})
    n = B.datasize[1]
    partialqr!(parent(R).data, n)
    materialize!(Ldiv(UpperTriangular(view(parent(R).data.data.data,OneTo(n),OneTo(n))), view(B.data,OneTo(n))))
    B
end

ldiv!(dest::AbstractVector, F::QR{<:Any,<:AdaptiveQRFactors}, b::AbstractVector) =
    ldiv!(F, copyto!(dest, b))
ldiv!(F::QR{<:Any,<:AdaptiveQRFactors}, b::AbstractVector) = ldiv!(F.R, lmul!(F.Q',b))
\(F::QR{<:Any,<:AdaptiveQRFactors}, B::AbstractVector) = ldiv!(F.R, F.Q'B)

factorize(A::BandedMatrix{<:Any,<:Any,<:OneToInf}) = qr(A)


text:
#!/usr/bin/env julia
include("commandUtils.jl")
activatePkg()
using ArgParse
using JLD
using Distributed
s = ArgParseSettings(description="Resonant Mode Direct sweep, implements algorithm 3 in Heider paper.")
@add_arg_table s begin
    "--simFile", "-s"
        help = "Simulation file to include, must set up boundary named simBoundary, can also edit or set k or other params to avoid the need to
                set in command arguments."
        required = true
    "-n"
        help = "Discretizaton N"
        arg_type = Int
        default = 50
    "--k0"
        help = "Search for resonances close to this k (can also set k0 in simFile)."
        arg_type = Complex{Float64}
    "--procs", "-p"
        help = "Number of processes to use in parallel"
        arg_type = Int
        default = convert(Int, Sys.CPU_THREADS/2)
end
arguments = parse_args(s)

@debug arguments

addprocs(arguments["procs"])
using BoundaryIntegralMethod
using BoundaryIntegralMethod.DirectSweep
n = arguments["n"]
k0 = arguments["k0"]

simFile = arguments["simFile"]
simName = split(basename(simFile), ".")[1]
if (!startswith(simFile, "/"))
  simFile = "$(pwd())/$simFile"
end
include(simFile)
b = simBoundary

(ks, xs) = directSweep(b, n, k0)
mkpath("output")
jldopen("output/sweep-$simName-$n-$(real(k0))_$(imag(k0))im.h5", "w") do file
  write(file, "ks/full", ks)
  write(file, "ks/Re", real(ks))
  write(file, "ks/Im", imag(ks))
  write(file, "xs/full", xs)
  write(file, "xs/Re", real(xs))
  write(file, "xs/Im", imag(xs))
end
exit(0)


text:
function fit!(model::SVDDClassifier, solver)
    debug(LOGGER, "[FIT] Fitting $(typeof(model)).")
    model.state == model_created && throw(ModelStateException(model.state, model_initialized))
    status = solve!(model, solver)
    model.R, model.const_term = get_R_and_const_term(model)
    model.state = model_fitted
    debug(LOGGER, "[FIT] $(typeof(model)) is now in state $(model.state).")
    return status
end

function predict(model::SVDDClassifier, target::Array{T,2}) where T <: Real
    model.state == model_fitted || throw(ModelStateException(model.state, model_fitted))
    @assert size(model.data, 1) == size(target, 1) "Dimension mismatch between traning data and target.
        target must be a column array, e.g., transpose([1.0 2.0])"

    α = get_alpha_prime(model)
    function predict_observation(z)
        kernel(model.kernel_fct, z, z) -
            2 * sum(α[i] * kernel(model.kernel_fct, model.data[:,i], z) for i in eachindex(α)) +
            model.const_term
    end
    return vec(sqrt.(mapslices(predict_observation, target, dims=1)) .- model.R)
end

function get_R_and_const_term(model::SVDDClassifier)
    @assert length(model.alpha_values) == size(model.K, 1)

    sv_idx = get_support_vectors(model)
    # take the maximum R over all support vectors for numerical reasons
    α = get_alpha_prime(model)

    const_term = sum(α[i] * α[j]*model.K[i,j] for i in eachindex(α) for j in eachindex(α))
    R_squared = 0.0
    tmp = 0.0
    for k in sv_idx
        tmp = model.K[k,k] - 2 * sum(α[i] * model.K[i,k] for i in eachindex(α)) + const_term
        R_squared = tmp > R_squared ? tmp : R_squared
    end
    return (sqrt(R_squared), const_term)
end

function invalidate_solution!(model::SVDDClassifier)
    model.alpha_values = []
    model.const_term = -Inf
    model.R = -Inf
    model.state = model_initialized
    return nothing
end


text:
struct MixRule{M<:ThermoModel,X} <: ThermoModel
    model::M
    xn::X
end

function mol_fraction(st::MixRule{M,X}) where {M,X}
    return st.xn
end

function thermomodel(st::MixRule{M,X}) where {M,X}
    return st.model
end

abstract type AbstractMixtureRule end

"""
    LinarMixing(p)

corresponds to the mix rule Σxᵢpᵢ.

"""
struct LinearMixing{M} <: AbstractMixtureRule
    prop::M
end

function mixrule(model::T,x) where T<: LinearMixing
    return dot(model.prop,x)
end
"""
    CuadraticMixing(p,f=(x,y)/2,Aij= nothing)

corresponds to the mix rule ΣAᵢⱼxᵢxⱼf(pᵢ,pⱼ)
where f is an averaging function. That means:
- f(pᵢ,pⱼ) = f(pⱼ,pᵢ)
- f(pₖ,pₖ) = pₖ

`Aij` is the weighed interaction matrix. if not provided, A[i,j] == Aᵢⱼ == 1
"""
struct QuadraticMixing{P,F,M} <: AbstractMixtureRule
    prop::P
    func::F
    Aij::M
end

function QuadraticMixing(prop)
    return QuadraticMixing(prop,(x,y)->0.5*(x+y),nothing)
end

function QuadraticMixing(prop,f)
    return QuadraticMixing(prop,f,nothing)
end

function mixrule(model::QuadraticMixing,x)
    return mixing_rule(model.func, x, model.prop,model.Aij)
end


text:
@testset "gewekediag.jl" begin
    samples = randn(100)

    @testset "results" begin
        @test @inferred(gewekediag(samples)) isa
              NamedTuple{(:zscore, :pvalue),Tuple{Float64,Float64}}
    end

    @testset "exceptions" begin
        for x in (-0.3, 0, 1, 1.2)
            @test_throws ArgumentError gewekediag(samples; first=x)
            @test_throws ArgumentError gewekediag(samples; last=x)
        end
        @test_throws ArgumentError gewekediag(samples; first=0.6, last=0.5)
    end
end


text:
module Runtests

using Test

@testset "Test ModiaPlot" begin
    include("test_01_OneScalarSignal.jl")
    include("test_02_OneScalarSignalWithUnit.jl")
    include("test_03_OneVectorSignalWithUnit.jl")
    include("test_04_ConstantSignalsWithUnit.jl")
    include("test_05_ArraySignalsWithUnit.jl")

    include("test_06_OneScalarMeasurementSignal.jl")
    include("test_07_OneScalarMeasurementSignalWithUnit.jl")

    include("test_20_SeveralSignalsInOneDiagram.jl")
    include("test_21_VectorOfPlots.jl")
    include("test_22_MatrixOfPlots.jl")
    include("test_23_MatrixOfPlotsWithTimeLabelsInLastRow.jl")
    include("test_24_Reuse.jl")
    include("test_25_SeveralFigures.jl")
    include("test_26_TooManyLegends.jl")

    include("test_40_Warnings.jl")

    include("test_51_OneScalarMonteCarloMeasurementsSignal.jl")
    include("test_52_MonteCarloMeasurementsWithDistributions.jl")

end

end


text:
using ITensorBenchmarks

# Run all of the benchmarks over the range of BLAS threads blas_num_threads and write the results
# By default, a verbose output is written to a log file
runbenchmarks(write_results = true, blas_num_threads = 1:6:Sys.CPU_THREADS)


text:
# Create example point data
x = 8 .* rand(100) .- 4
y = 8 .* rand(100) .- 4
z = sin.(x) + cos.(y)
# Triangular surface plot
trisurf(x, y, z)


text:
import LinearAlgebra: BlasInt

export DMDParams, QRStore

# This file is available under the terms of the MIT License

@doc """
Types for describing DMD fit and storing necessary
variables
"""

"""
DMDParams - stores the description of the problem
including data and sizes. We
only provide constructors for BLAS friendly types
for now
"""

mutable struct DMDParams{T<:AbstractFloat}
    # dimensions
    m::Integer
    n::Integer
    k::Integer
    #
    # data
    X::Matrix{Complex{T}}
    x::Array{Vector{Complex{T}},1}
    t::Vector{Complex{T}}
    #
    # loss function and gradient
    lossFunc::Function
    lossGrad::Function
    #
    # variables
    a::Vector{Complex{T}}
    P::Matrix{Complex{T}}
    B::Matrix{Complex{T}}
    b::Array{Vector{Complex{T}},1}
    R::Matrix{Complex{T}}
    r::Array{Vector{Complex{T}},1}
    #
    ar::Vector{T}
    br::Array{Vector{T},1}
    # l2 direct storage
    qra::Union{Nothing,Vector{Complex{T}}}
    qrfact::Union{Nothing,QRPivoted{Complex{T},Array{Complex{T},2}}}
    qrTemp::Union{Nothing,Matrix{Complex{T}}}
    qrtemp::Union{Nothing,Array{Vector{Complex{T}},1}}

    # temp variables
    tv::Vector{Complex{T}}
    tM::Matrix{Complex{T}}
    tcnorm::Vector{T}
    idsort::Vector{S} where S <: Integer
    #
    # inner solve options
    inner_solver::Any
    inner_opts::Any
    inner_directl2::Bool

    # trimming
    nkeep::Integer
    ikeep::Any
end

# constructors
function DMDParams(k, X, t, lossFunc, lossGrad;
                   inner_directl2::Bool=false,
                   inner_opts=Optim.Options(),
                   inner_solver=BFGS(),
                   nkeep=size(X,2))
    m,n = size(X);
    x   = col_view(X);
    #
    Tx  = eltype(X);
    Tt  = eltype(t);
    if Tx != Tt
        t = convert(Vector{Tx}, t);
    end
    #
    a = zeros(Tx, k);
    P = zeros(Tx, m, k);
    B = zeros(Tx, k, n);
    R = zeros(Tx, m, n);
    #
    b = col_view(B);
    r = col_view(R);
    #
    ar = col_view_real(a);
    br = col_view_real(B);
    #
    tv = zeros(Tx, k);
    tM = zeros(Tx, k, n);
    tMkeep = zeros(Tx, k, nkeep);
    tcnorm = zeros(real(Tx),n)
    idsort = Array{Int64}(undef,n)
    #

    qra = nothing
    qrfact = nothing
    qrTemp = nothing
    qrtemp = nothing

    if inner_directl2
        qra = copy(a); qra .= NaN;
        qrfact = QRPivoted(copy(P),Array{Tx}(undef,min(m,k)),
                           zeros(BlasInt,k))
        qrTemp = copy(R)
        qrtemp = col_view(qrTemp)
    end

    ikeep = collect(1:nkeep)

    return DMDParams(m, n, k, X, x, t, lossFunc, lossGrad,
                     a, P, B, b,
                     R, r, ar, br,
                     qra, qrfact, qrTemp, qrtemp,
                     tv, tM, tcnorm, idsort,
                     inner_solver,
                     inner_opts, inner_directl2,
                     nkeep, ikeep)
end


text:
using LinearAlgebra, Statistics

export embed

export transform, scale, rotate, translate, recenter
export edge_length, ensure_embed

const DEFAULT_MARKER_SIZE = 6

# A graph embedding includes these items in the graphs' cache:
# :xy           a dictionary from vertices of [x,y] coordinates
# :vcolor       a dictionay from vertices to color specifications (default all :white)
# :vsize        marker size for vertices  (default: DEFAULT_MARKER_SIZE)
# :line_color   color to draw edges and boundaries of vertices (default :black)

"""
`_new_embed(G)` sets the graph up with a default embedding.

* :vcolor is set to :white for all vertices
* :line_color is set to :black
* :vsize is set to DEFAULT_MARKER_SIZE
* :xy is set to a circular embedding
"""
function _new_embed(G::SimpleGraph{T}) where {T}
    cache_save(G, :vsize, DEFAULT_MARKER_SIZE)

    G.cache[:vcolor] = Dict{T,Any}()
    for v in G.V
        G.cache[:vcolor][v] = :white
    end

    cache_save(G, :line_color, :black)
    n = NV(G)

    G.cache[:xy] = _circular_xy(G)

    nothing
end

"""
`ensure_embed(G)` gives `G` a default embedding if it
doesn't already have an embedding.
"""
function ensure_embed(G::SimpleGraph)
    if !cache_check(G, :xy)
        _new_embed(G)
    end
end

function private_adj(G::SimpleGraph)
    n = NV(G)
    A = zeros(Int, n, n)
    vv = vlist(G)
    for i = 1:n-1
        u = vv[i]
        for j = (i+1):n
            w = vv[j]
            if has(G, u, w)
                A[i, j] = 1
                A[j, i] = 1
            end
        end
    end
    return A, vv
end

"""
`embed(G)` creates a new embedding for `G`. The full call is
```
embed(G,algorithm;args...)
```
The `symbol` algorithm indicates the embedding algorithm.
The `args` are a collection of possible arguments to be sent
to the algorithm.

The `algorithm` defaults to `:circular` and may be one of the following:

* `:circular`: arrange the vertices evenly in a circle.
* `:random`: arrange the vertices randomly.
* `:spring`: use the `spring` layout from `GraphLayout`. Optional argument:
  * `iterations`.
* `:stress`: use the `stress` layout from `GraphLayout`.
* `:spectral`: use the `spectral` embedding. Optional arguments:
  * `xcol` -- which eigenvalue to use for the `x` coordinate.
  * `ycol` -- which eigenvalue to use for the `y` coordinate
* `:tutte` -- create a Tutte embedding using a longest face (assuming the graph has a rotation system)
  * `outside` [optional] -- a list of vertices to be the outer face of the embedding

Note that if the graph already has (say) an embedding, that embedding may
be used as the starting point for one of the algorithms.
"""
function embed(G::SimpleGraph, algorithm::Symbol = :circular; args...)
    arg_dict = list2dict(collect(args))
    n = NV(G)
    m = NE(G)

    ## FULL LIST OF POSSIBLE ARGUMENTS PRE-PARSED HERE

    verbose::Bool = true
    if haskey(arg_dict, :verbose)
        verbose = arg_dict[:verbose]
    end

    iterations::Int = 0
    if haskey(arg_dict, :iterations) && arg_dict[:iterations] > 0
        iterations = arg_dict[:iterations]
    end

    tolerance = 0.001
    if haskey(arg_dict, :tolerance) && arg_dict[:tolerance] > 0
        tolerance = arg_dict[:tolerance]
    end

    xcol::Int = 2
    if haskey(arg_dict, :xcol) && arg_dict[:xcol] > 0
        xcol = arg_dict[:xcol]
    end

    ycol::Int = 3
    if haskey(arg_dict, :ycol) && arg_dict[:ycol] > 0
        xcol = arg_dict[:ycol]
    end

    # See if there already is an xy-embedding
    if !cache_check(G, :xy)
        _new_embed(G)
        if algorithm == :circular # no need to do it again
            return nothing
        end
    end

    if algorithm == :circular
        G.cache[:xy] = _circular_xy(G)
        return nothing
    end

    if algorithm == :random
        G.cache[:xy] = _random_xy(G)
        return nothing
    end

    if algorithm == :spring
        if n < 3 || m == 0
            embed(G)
            return
        end
        if iterations <= 0
            G.cache[:xy] = _spring(G)
        else
            G.cache[:xy] = _spring(G, iterations)
        end
        return nothing
    end

    if algorithm == :stress
        if n < 3 || m == 0
            embed(G)
            return
        end
        G.cache[:xy] = _stress(G)
        return nothing
    end

    if algorithm == :combined
        if n < 3 || m == 0
            embed(G)
            return
        end
        embed(G, :spring, iterations = iterations)
        scale(G)
        embed(G, :stress)
        return nothing
    end

    if algorithm == :spectral
        _spectral(G, xcol, ycol)
        return nothing
    end

    if algorithm == :tutte
        if haskey(args, :outside)
            _tutte(G, args[:outside])
        else
            _tutte(G)
        end
        return nothing
    end

    @warn "Unknown embedding algorithm; no action taken"
    return nothing
end

"""
`embed(G,d)` specifies an embedding of the graph `G` with
a dictionary `d` that maps vertices to coordinates (as two
dimensional vectors `[x,y]`).
"""
function embed(G::SimpleGraph, xy::Dict)
    ensure_embed(G)
    G.cache[:xy] = deepcopy(xy)
    nothing
end

"""
`list2dict(list)` takes a list of `(Symbol,Any)` pairs and
converts them to a dictionary mapping the symbols to their
associated values.
"""
function list2dict(list::Vector)
    d = Dict{Symbol,Any}()
    for it in list
        x, y = it
        d[x] = y
    end
    return d
end

"""
`_circular_xy(G)` creates a standard circular embedding.
"""
function _circular_xy(G::SimpleGraph{T})::Dict where {T}
    n = NV(G)
    xy = Dict{T,Vector{Float64}}()
    if n == 0
        return xy
    end

    r = sqrt(n)
    VV = vlist(G)
    theta = 2 * pi / n
    for k = 1:n
        t = (k - 1) * theta
        x = r * sin(t)
        y = r * cos(t)
        v = VV[k]
        xy[v] = [x, y]
    end
    return xy
end

"""
`_random_xy(G)` creates a random xy-embedding for `G`
"""
function _random_xy(G::SimpleGraph{T})::Dict where {T}
    rootn = sqrt(NV(G))
    xy = Dict{T,Vector{Float64}}()
    for v in G.V
        xy[v] = [rand(), rand()] * rootn
    end
    _recenter(xy)
    return xy
end

"""
`edge_length(G)` returns an array containing the
lengths of the edges in the current embedding of `G`.
"""
edge_length(G::SimpleGraph) = [edge_length(G, ee) for ee in G.E]

"""
`edge_length(G,e)` gives the distance between the
embedded endpoints of the edge `e` in the drawing `G`.
"""
function edge_length(G::SimpleGraph{T}, v, w) where {T}
    p1 = getxy(G, v)
    p2 = getxy(G, w)
    return norm(p1 - p2)
end
edge_length(G::SimpleGraph, ee) = edge_length(G, ee[1], ee[2])

"""
`scale(G,m)` multiplies all coordinates in the graph's drawing by
`m`. If `m` is omitted, the drawing is rescaled so that the average
length of an edge equals 1.
"""
function scale(G::SimpleGraph, m::T) where {T<:Real}
    ensure_embed(G)
    _scale(G.cache[:xy], m)
end

function scale(G::SimpleGraph)
    if NE(G) == 0
        return
    end
    L = mean(edge_length(G))
    if L != 0
        _scale(G.cache[:xy], 1 / L)
    else
        @warn "Cannot scale $G"
    end
    nothing
end

include("transforms.jl")
include("getset.jl")
include("tutte.jl")
include("my_spring.jl")
include("my_stress.jl")
include("colorize.jl")
include("graffle.jl")
include("spectral.jl")
# include("geogebra.jl")
include("tikz.jl")


text:
# Tests for the linear discriminant classifier
function t_lindisc()

Ac = j4pr.DataGenerator.fish(20)   # classification dataset

tol  = 1e-6; # tolerance when comparing results
Wclass = j4pr.lindisc(rand(),rand())

# Test classification
Test.@test try
	wt1 = Ac |> Wclass
	wt2 = j4pr.strip(Ac) |> Wclass
	result = Ac |> wt1
	result2 = +Ac |> wt2
	sum(abs.(result.x - result2) .>= tol) > 0 ? false : true
	catch
		false
	end

# Test classification on a 1-D dataset
A1d = j4pr.datacell([rand(100);2*rand(100)],[zeros(100);ones(100)])
wt1 = A1d |> Wclass
Test.@test try
	wt2 = j4pr.strip(A1d) |> Wclass
	result = A1d |> wt1
	result2 = +A1d |> wt2
	sum(abs.(result.x - result2) .>= tol) > 0 ? false : true
	catch
		false
	end

# Test Printer
Test.@test try
	buf = IOBuffer()
	Base.show(buf,wt1.x.data)
	true #works
	catch
		false
	end

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder
import Pkg: PackageSpec
import Pkg.Types: VersionSpec

name = "libpolymake_julia"
upstream_version = v"0.4.1"
version = VersionNumber(upstream_version.major, upstream_version.minor, upstream_version.patch * 100 + 10 + julia_version.minor)

# Collection of sources required to build libpolymake_julia
sources = [
    ArchiveSource("https://github.com/oscar-system/libpolymake-julia/archive/v$(upstream_version).tar.gz",
                  "f7733a8eb2b2e75d63db6cb6b035b9fae24ce7f6c2992253a3df0529a914ad37"),
]

# Bash recipe for building across all platforms
script = raw"""
# remove $libdir from LD_LIBRARY_PATH as this causes issues with perl
if [[ -n "$LD_LIBRARY_PATH" ]]; then
LD_LIBRARY_PATH=$(echo -n $LD_LIBRARY_PATH | sed -e "s|[:^]$libdir\w*|:|g")
fi

cmake libpolymake-j*/ -B build \
   -DJulia_PREFIX="$prefix" \
   -DCMAKE_INSTALL_PREFIX="$prefix" \
   -DCMAKE_FIND_ROOT_PATH="$prefix" \
   -DCMAKE_TOOLCHAIN_FILE="${CMAKE_TARGET_TOOLCHAIN}" \
   -DCMAKE_BUILD_TYPE=Release

VERBOSE=ON cmake --build build --config Release --target install -- -j${nproc}

install_license $WORKSPACE/srcdir/libpolymake-j*/LICENSE.md
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = [
    Platform("x86_64", "linux"; libc="glibc"),
    Platform("x86_64", "macos"),
]
platforms = expand_cxxstring_abis(platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("polymake_run_script", :polymake_run_script),
    LibraryProduct("libpolymake_julia", :libpolymake_julia),
    FileProduct("share/libpolymake_julia/type_translator.jl",:type_translator),
]

# Dependencies that must be installed before this package can be built
dependencies = [
    BuildDependency(PackageSpec(name="libjulia_jll", version=julia_version)),
    BuildDependency(PackageSpec(name="GMP_jll", version=v"6.1.2")),
    BuildDependency(PackageSpec(name="MPFR_jll", version=v"4.0.2")),
    Dependency("CompilerSupportLibraries_jll"),
    Dependency("libcxxwrap_julia_jll", VersionNumber(0, 8, julia_version.minor)),
    Dependency("polymake_jll", v"400.400.0"; compat = "~400.400.0"),
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies;
    preferred_gcc_version=v"8",
    julia_compat = "$(julia_version.major).$(julia_version.minor)")


text:
# Non-bipartite: https://theory.stanford.edu/~jvondrak/CS369P/lec5.pdf

mutable struct PerfectBipartiteMatchingLPSolver <: PerfectBipartiteMatchingSolver
  solver

  # Optimisation model (reused by solve_linear).
  model # ::Model
  x # ::Matrix{Variable}

  function PerfectBipartiteMatchingLPSolver(solver)
    return new(solver, nothing, nothing)
  end
end

function build!(solver::PerfectBipartiteMatchingLPSolver, n_arms::Int)
  # Build the optimisation model behind solve_linear.
  solver.model = Model(solver.solver)
  indices = collect((i, j) for i in 1:n_arms, j in 1:n_arms)
  solver.x = @variable(solver.model, [indices], lower_bound=0, upper_bound=1)

  for i in 1:n_arms # Left nodes.
    @constraint(solver.model, sum(solver.x[(i, j)] for j in 1:n_arms) <= 1)
  end
  for j in 1:n_arms # Right nodes.
    @constraint(solver.model, sum(solver.x[(i, j)] for i in 1:n_arms) <= 1)
  end
end

has_lp_formulation(::PerfectBipartiteMatchingLPSolver) = true

function get_lp_formulation(solver::PerfectBipartiteMatchingLPSolver, rewards::Dict{Tuple{Int, Int}, Float64})
  return solver.model,
    sum(rewards[(i, j)] * solver.x[(i, j)] for (i, j) in keys(rewards)),
    Dict{Tuple{Int, Int}, JuMP.VariableRef}((i, j) => solver.x[(i, j)] for (i, j) in keys(rewards))
end

function solve_linear(solver::PerfectBipartiteMatchingLPSolver, rewards::Dict{Tuple{Int, Int}, Float64})
  m, obj, vars = get_lp_formulation(solver, rewards)
  @objective(m, Max, obj)

  set_silent(m)
  optimize!(m)

  if termination_status(m) != MOI.OPTIMAL
    return Tuple{Int, Int}[]
  end
  return Tuple{Int, Int}[(i, j) for (i, j) in keys(rewards) if value(vars[(i, j)]) > 0.5]
end


text:
# Ahorn's installation script. Connects to the internet using Pkg3 to install it.

if VERSION < v"1.3"
    println("""Ahorn, and thus this installer, require Julia version 1.3 or later to run.
    The version you installed is $VERSION. Please update to a more recent version.
    Press Enter to quit the installer.""")
    readline()
    exit()
end

using Pkg

install_or_update(url::String, pkg::String) = try
    if Pkg.installed()[pkg] !== nothing
        println("Updating $pkg...")
        Pkg.update(pkg)
    end
catch err
    println("Installing $pkg...")
    Pkg.add(PackageSpec(url=url))
end

if Sys.iswindows()
    installpath = joinpath(ENV["LOCALAPPDATA"], "Ahorn", "env")
else
    installpath = joinpath(get(ENV, "XDG_CONFIG_HOME", joinpath(get(ENV, "HOME", ""), ".config")) , "Ahorn", "env")
end

println("""
======================
Hey, thanks for giving Ahorn a try! Ahorn is currently still in a very early state. If you find bugs, please report them!
======================
""")

println("This installer is going to install or update Maple and Ahorn via the package manager as well as create a file in the same directory as this installer file, are you okay with that? [y/N]")

if !occursin(r"^[Yy]", readline())
    println("Exiting installer.")
    exit()
end

println("""

======================
The installer will now download and install required dependencies as well as the program itself. This will likely take a few minutes, so grab yourself a glass of juice while you wait.
We will tell you when the installation process is done, any warnings (not errors) in the console are fine and can be ignored.
======================
""")

mkpath(installpath)

println("Environment path: " * installpath)

Pkg.activate(installpath)

install_or_update("https://github.com/CelestialCartographers/Maple.git", "Maple")
install_or_update("https://github.com/kingW3/Ahorn.git", "Ahorn")

println("""

======================
Precompiling a few dependencies. This may take a while, so get yourself some cheese and a cup of fennel tea.
======================
""")

sleep(1)
Pkg.API.precompile()
import Ahorn

pkgdir = normpath(joinpath(dirname(pathof(Ahorn)), ".."))
println("Package path: " * pkgdir)

launchpath = nothing
if Sys.iswindows()
    launchpath = joinpath(@__DIR__, "ahorn.bat")
    cp(joinpath(pkgdir, "ahorn.bat"), launchpath, force=true)
else
    #symlink(joinpath(pkgdir, "src", "run_ahorn.jl"), joinpath(@__DIR__, "ahorn.jl"))
    launchpath = joinpath(@__DIR__, "ahorn.sh")
    if isfile(launchpath)
        println("Removing old ahorn.sh...")
        rm(launchpath)
    end
    #symlink(joinpath(pkgdir, "ahorn"), launchpath)
    cp(joinpath(pkgdir, "ahorn"), launchpath, force=true)
end

if Sys.iswindows()
    println("""

Done! Ahorn should be installed now. Run $launchpath to launch it.
Note that it will take quite a while to launch the first time as its dependencies compile, so please be patient.

Thanks for giving Ahorn a try!

Press Enter to quit the installer.
""")
else
    println("""

Done! Ahorn should be installed now. If Julia is in your path, just run $launchpath to launch it.
Note that it will take quite a while to launch the first time as its dependencies compile, so please be patient.

Thanks for giving Ahorn a try!

Press Enter to quit the installer.
""")
end
readline()


text:
const RATIO_TABLE = [
    # (numerator, denominator)-pair
    ("branch-misses", "branches"),
    ("L1-dcache-load-misses", "L1-dcache-loads"),
    ("L1-icache-load-misses", "L1-icache-loads"),
    ("LLC-load-misses", "LLC-loads"),
    ("dTLB-load-misses", "dTLB-loads"),
    ("iTLB-load-misses", "iTLB-loads"),
    ("stalled-cycles-frontend", "cycles"),
    ("stalled-cycles-backend", "cycles"),
]

const OTHER_KEYS = [
    # "task-clock"  TODO: float
    "context-switches",
    "cpu-migrations",
    "page-faults",
    "instructions",
]

as_property_symbol(str::AbstractString) = Symbol(lowercase(replace(str, "-" => "_")))

const ALL_KEYS = unique!(vcat(first.(RATIO_TABLE), last.(RATIO_TABLE), OTHER_KEYS))

for key in ALL_KEYS
    val = Val{as_property_symbol(key)}
    @eval getvalue(t::TrialPerf, ::$val)::Union{Int,Nothing} =
        @something(get(t, $key, nothing), return)
end

for (num, den) in RATIO_TABLE
    numval = Val{as_property_symbol(num)}
    denval = Val{as_property_symbol(den)}
    @eval function getratio(t::TrialPerf, ::$numval)::Union{Float64,Nothing}
        num = @something(getvalue(t, $numval()), return)
        den = @something(getvalue(t, $denval()), return)
        return num / den
    end
end

# perf calls them (and the ratios) "shaddow stats"
const DERIVED_NAMES = [
    :instructions_per_cycle,  # "insn per cycle"
    :stalled_cycles_per_instructions,  # "stalled cycles per insn"
]

const ALL_NAMES = (as_property_symbol.(ALL_KEYS)..., DERIVED_NAMES...)
const RATIO_NAMES = (as_property_symbol.(first.(RATIO_TABLE))...,)

struct RatioAccessor
    trial::TrialPerf
end

struct PercentAccessor
    trial::TrialPerf
end

Base.getproperty(accessor::RatioAccessor, name::Symbol) =
    getratio(getfield(accessor, :trial), Val{name}())
Base.getproperty(accessor::PercentAccessor, name::Symbol) =
    @something(getratio(getfield(accessor, :trial), Val{name}()), return) * 100

function Base.propertynames(
    accessor::Union{RatioAccessor,PercentAccessor},
    private::Bool = false,
)
    pubnames = RATIO_NAMES
    if private
        return (pubnames..., fieldnames(typeof(accessor)))
    else
        return pubnames
    end
end

getratio(::TrialPerf, @nospecialize(_val::Val{name})) where {name} =
    error(name, " is not a name of suportd ratio ", RATIO_NAMES)

function getvalue(trial::TrialPerf, ::Val{:instructions_per_cycle})
    num = @something(trial.instructions, return)
    den = @something(trial.cycles, return)
    return num / den
end

# https://github.com/torvalds/linux/blob/v5.15/tools/perf/util/stat-shadow.c#L979-L991
function getvalue(trial::TrialPerf, ::Val{:stalled_cycles_per_instructions})
    f = @something(trial.stalled_cycles_frontend, return)
    b = @something(trial.stalled_cycles_backend, return)
    num = max(f, b)
    den = @something(trial.instructions, return)
    return num / den
end


text:
function degrees(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    if issymmetric(adj)
        d = vec(sum(adj, dims=1))
    else
        if dir == :out
            d = vec(sum(adj, dims=1))
        elseif dir == :in
            d = vec(sum(adj, dims=2))
        elseif dir == :both
            d = vec(sum(adj, dims=1)) + vec(sum(adj, dims=2))
        else
            throw(DomainError(dir, "invalid argument, only accept :in, :out and :both"))
        end
    end
    d
end

function degree_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    d = degrees(adj, T, dir=dir)
    return SparseMatrixCSC(T.(diagm(0=>d)))
end

function inv_sqrt_degree_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    d = degrees(adj, T, dir=dir).^(-0.5)
    return SparseMatrixCSC(T.(diagm(0=>d)))
end

function laplacian_matrix(adj::AbstractMatrix, T::DataType=eltype(adj); dir::Symbol=:out)
    degree_matrix(adj, T, dir=dir) - SparseMatrixCSC(T.(adj))
end

function normalized_laplacian(adj::AbstractMatrix, T::DataType=eltype(adj))
    inv_sqrtD = inv_sqrt_degree_matrix(adj, T, dir=:both)
    I - inv_sqrtD * SparseMatrixCSC(T.(adj)) * inv_sqrtD
end

function neighbors(adj::AbstractMatrix, T::DataType=eltype(adj))
    n = size(adj,1)
    @assert n == size(adj,2) "adjacency matrix is not a square matrix."
    A = (adj .!= zero(T))
    if !issymmetric(adj)
        A = A .| A'
    end
    indecies = collect(1:n)
    ne = Vector{Int}[indecies[view(A, :, i)] for i = 1:n]
    return ne
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 53584796-6447-4a4b-acc6-ca4428768a93
if isdefined(Main, :PlutoRunner)
	# If this file runs under Pluto, we activate a temporary
	# environment and add all necessary packages, including PyPlot
    using Pkg
    Pkg.activate(mktempdir())
	Pkg.add(["PyPlot","Triangulate","GridVisualize","SimplexGridFactory"])
end

# ╔═╡ d432ad64-f91f-11ea-2e48-4bc7472ac64c
begin
    # Trick pyplot into not using , as floating point decimal delimiter
    # in certain language enviromnents
	ENV["LC_NUMERIC"]="C"
    ENV["MPLBACKEND"]="agg"
	using SimplexGridFactory,GridVisualize,Triangulate,PyPlot
	PyPlot.svg(true)
end;

# ╔═╡ 940b1996-fe9d-11ea-2fa4-8b72bee62b76
md"""
## Pluto notebook showing the use of SimplexGridBuilder
"""

# ╔═╡ bc6135b8-451f-11eb-106a-2f9c9a8229d7
md"""
__Prepare the environment__:
"""

# ╔═╡ f32d9f04-f923-11ea-3a4a-53cc3df5642c
html"""<hr>"""

# ╔═╡ fd27b44a-f923-11ea-2afb-d79f7e62e214
md"""
### Using the SimplexGridBuilder

__Local refinement center:__
"""

# ╔═╡ b12838f0-fe9c-11ea-2939-155ed907322d
center=[0.8,0.2]

# ╔═╡ d5d8a1d6-fe9d-11ea-0fd8-df6e81492cb5
md"""
__Local refinement callback:__
"""

# ╔═╡ aae2e82a-fe9c-11ea-0427-593f8d2c7746
function unsuitable(x1,y1,x2,y2,x3,y3,area)
        bary_x=(x1+x2+x3)/3.0
        bary_y=(y1+y2+y3)/3.0
        dx=bary_x-center[1]
        dy=bary_y-center[2]
        qdist=dx^2+dy^2
        qdist>1.0e-5 && area>0.1*qdist
end;

# ╔═╡ 1ae86964-fe9e-11ea-303b-65bb128384a5
md"""
__Set up a SimplexGridBuilder:__
"""

# ╔═╡ 511b26c6-f920-11ea-1228-51c3750f495c
function make_builder()
	factory=SimplexGridBuilder(Generator=Triangulate)

	#  Specfy points
	p1=point!(factory,0,0)
	p2=point!(factory,1,0)
	p3=point!(factory,1,1)
	p4=point!(factory,0,1.5)

	# Specify outer boundary
	facetregion!(factory,1)
	facet!(factory,p1,p2)
	facet!(factory,p2,p3)
	facet!(factory,p3,p4)
	facet!(factory,p1,p4)

	# specify interionr boundary
	facetregion!(factory,2)
	facet!(factory,p1,p3)

	# Coarse elements in upper left region
	cellregion!(factory,1)
	regionpoint!(factory,0.1,0.5)

	# Fine elements in lower right region
	cellregion!(factory,2)
	maxvolume!(factory,0.01)
	regionpoint!(factory,0.9,0.5)

	# Activate unsuitable callback
	options!(factory,unsuitable=unsuitable)
	factory
end

# ╔═╡ dd85d88e-44b3-11eb-3e31-c7dcbe07d0de
builder=make_builder()

# ╔═╡ 9735eba3-a51b-4870-8557-c3a2985b2297
md"""
__Testplot with input and output__
"""

# ╔═╡ 8f0bd5c0-f920-11ea-3b1c-db90fc95f990
builderplot(builder,Plotter=PyPlot,resolution=(600,600))

# ╔═╡ 54677c86-fa91-11ea-3518-3b8aba2c8488
md"""
These are the Triangle control flags created from the default options:
"""

# ╔═╡ 1067166c-451f-11eb-2916-c94382a220a3
SimplexGridFactory.makeflags(builder.options,:triangle)

# ╔═╡ d2129483-285b-49a2-a11d-886956146b85
md"""
__Create a simplex grid from the builder__
"""

# ╔═╡ ac93589b-6315-4677-9542-c0a2333f1755
simplexgrid(builder)

# ╔═╡ Cell order:
# ╟─940b1996-fe9d-11ea-2fa4-8b72bee62b76
# ╟─bc6135b8-451f-11eb-106a-2f9c9a8229d7
# ╠═53584796-6447-4a4b-acc6-ca4428768a93
# ╠═d432ad64-f91f-11ea-2e48-4bc7472ac64c
# ╟─f32d9f04-f923-11ea-3a4a-53cc3df5642c
# ╟─fd27b44a-f923-11ea-2afb-d79f7e62e214
# ╠═b12838f0-fe9c-11ea-2939-155ed907322d
# ╟─d5d8a1d6-fe9d-11ea-0fd8-df6e81492cb5
# ╠═aae2e82a-fe9c-11ea-0427-593f8d2c7746
# ╟─1ae86964-fe9e-11ea-303b-65bb128384a5
# ╠═511b26c6-f920-11ea-1228-51c3750f495c
# ╠═dd85d88e-44b3-11eb-3e31-c7dcbe07d0de
# ╟─9735eba3-a51b-4870-8557-c3a2985b2297
# ╠═8f0bd5c0-f920-11ea-3b1c-db90fc95f990
# ╟─54677c86-fa91-11ea-3518-3b8aba2c8488
# ╠═1067166c-451f-11eb-2916-c94382a220a3
# ╟─d2129483-285b-49a2-a11d-886956146b85
# ╠═ac93589b-6315-4677-9542-c0a2333f1755


text:
module BLASBenchmarksGPU

import BenchmarkTools
import CUDA
import DataFrames
import GemmKernels
import KernelAbstractions
import Plots
import PyPlot
import ProgressMeter
import Random
import Statistics
import StatsPlots
import Tullio

using CUDA               # this `using` statement is required to use Tullio on the GPU
using KernelAbstractions # this `using` statement is required to use Tullio on the GPU

export benchmark_result_df
export benchmark_result_types
export plotbench
export runbench

include("types.jl")

include("benchconfig.jl")
include("benchmarkresult.jl")
include("plotting.jl")
include("runbenchmark.jl")

end # end module BLASBenchmarksGPU


text:
function compute_lagrangian(data::Dict)
  x        = data["x"];
  d        = data["d"];
  s        = data["s"];
  lambda_s = data["lambda_s"];
  mu_s     = data["mu_s"];
  z        = data["z"];
  t        = data["t"];
  lambda_t = data["lambda_t"];
  mu_t     = data["mu_t"];
  beta     = data["beta"];

  numPatches = size(x, 2);
  dictSize = size(d, 1);

  loss = 0.5*vecnorm( x-d'*z ) ^2;
  spar = beta*vecnorm(t, 1);
  cons = mu_t*0.5*vecnorm( z-t )^2 + [lambda_t[:]'*( z[:]-t[:] )][1];
  cons += mu_s*0.5*vecnorm( d-s )^2 + [lambda_s[:]'*( d[:]-s[:] )][1];
  obj = loss + spar + cons;

  return (obj, loss, spar);
end

max_stars_repo_name: bkong/FSC

max_issues_repo_name: bkong/FSC

max_forks_repo_name: bkong/FSC


text:
import .StatsFuns
using .StatsFuns: xlogx, logistic, logit, log1psq, log1pexp, logsumexp

@adjoint function xlogx(x::Real)
    y = xlogx(x)
    return y, function(Δ::Real)
        return (x > zero(x) ? Δ * (log(x) + one(y)) : zero(y),)
    end
end

@adjoint function logistic(x::Real)
    y = logistic(x)
    return y, Δ->(Δ * y * (1 - y),)
end

@adjoint logit(x::Real) = logit(x), Δ->(Δ / (x * (1 - x)),)

@adjoint log1psq(x::Real) = log1psq(x), Δ->(Δ * 2x / (1 + abs2(x)),)

@adjoint function log1pexp(x::Float64)
    return log1pexp(x), Δ->(Δ * (x < 18.0 ? logistic(x) : x < 33.3 ? 1 - exp(-x) : 1),)
end

@adjoint function log1pexp(x::Float32)
    return log1pexp(x), Δ->(Δ * (x < 9f0 ? logistic(x) : x < 16f0 ? 1 - exp(-x) : 1),)
end

@adjoint function logsumexp(X::AbstractArray{<:Real}; dims=:)
    lse = logsumexp(X; dims=dims)
    return lse, Δ -> (Δ .* exp.(X .- lse),)
end


text:
using Gadfly
using RDatasets

plot(data("datasets", "iris"), x="Sepal.Length", y="Sepal.Width", Geom.point)


text:
using LinearAlgebra

struct ValueOne; end
ValueOne()

# Compute X <- a X + b I.
function matfun_axpby!(X,a,b,Y::UniformScaling)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:n
        X[i,i]+=(b isa ValueOne) ? 1 : b
    end
end

# Compute X <- a X + b Y.
function matfun_axpby!(X,a,b,Y)
    m,n=size(X)
    if ~(a isa ValueOne)
        rmul!(X,a)
    end
    @inbounds for i=1:m
        @inbounds for j=1:n
            if (b isa ValueOne)
                X[i,j]+=Y[i,j]
            else
                X[i,j]+=b*Y[i,j]
            end
        end
    end
end

@inline function exp_native_jl_rho0_95(A)
    T=promote_type(eltype(A),Float64)
    A_copy=similar(A,T); A_copy .= A;
    return exp_native_jl_rho0_95!(A_copy)
end

@inline function exp_native_jl_rho0_95!(A)
    T=promote_type(eltype(A),Float64) # Make it work for many 'bigger' types (matrices and scalars)
    # max_memslots=6
    n=size(A,1)
    # The first slots are precomputed nodes [:A]
    memslots2 = similar(A,T)
    memslots3 = similar(A,T)
    memslots4 = similar(A,T)
    memslots5 = similar(A,T)
    memslots6 = similar(A,T)
    # Assign precomputed nodes memslots
    memslots1=A # overwrite A
    # Uniform scaling is exploited.
    # No matrix I explicitly allocated.
    value_one=ValueOne()
    # Computation order: A2 A4 A6 Ua U V Z X P
    # Computing A2 with operation: mult
    mul!(memslots2,memslots1,memslots1)
    # Computing A4 with operation: mult
    mul!(memslots3,memslots2,memslots2)
    # Computing A6 with operation: mult
    mul!(memslots4,memslots2,memslots3)
    # Computing Ua = x*I+x*A2+x*A4+x*A6
    coeff1=8.64864e6
    coeff2=277200.0
    coeff3=1512.0
    coeff4=1.0
    memslots5 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4
    mul!(memslots5,true,I*coeff1,true,true)
    # Computing U with operation: mult
    mul!(memslots6,memslots5,memslots1)
    # Deallocating Ua in slot 5
    # Deallocating A in slot 1
    # Computing V = x*I+x*A2+x*A4+x*A6
    coeff1=1.729728e7
    coeff2=1.99584e6
    coeff3=25200.0
    coeff4=56.0
    # Smart lincomb recycle A2
    memslots2 .= coeff2.*memslots2 .+ coeff3.*memslots3 .+ coeff4.*memslots4
    mul!(memslots2,true,I*coeff1,true,true)
    # Deallocating A4 in slot 3
    # Deallocating A6 in slot 4
    # Computing Z = x*V+x*U
    coeff1=1.0
    coeff2=-1.0
    memslots1 .= coeff1.*memslots2 .+ coeff2.*memslots6
    # Computing X = x*V+x*U
    coeff1=1.0
    coeff2=1.0
    # Smart lincomb recycle V
    memslots2 .= coeff1.*memslots2 .+ coeff2.*memslots6
    # Deallocating U in slot 6
    # Computing P with operation: ldiv
    memslots3 .=memslots1\memslots2
    # Deallocating Z in slot 1
    # Deallocating X in slot 2
    return memslots3 # Returning P
end

max_stars_repo_name: matrixfunctions/GraphMatFunData

max_issues_repo_name: matrixfunctions/GraphMatFunData

max_forks_repo_name: matrixfunctions/GraphMatFunData


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule <: SwaggerModel
    apiGroups::Any # spec type: Union{ Nothing, Vector{String} } # spec name: apiGroups
    clusterScope::Any # spec type: Union{ Nothing, Bool } # spec name: clusterScope
    namespaces::Any # spec type: Union{ Nothing, Vector{String} } # spec name: namespaces
    resources::Any # spec type: Union{ Nothing, Vector{String} } # spec name: resources
    verbs::Any # spec type: Union{ Nothing, Vector{String} } # spec name: verbs

    function IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule(;apiGroups=nothing, clusterScope=nothing, namespaces=nothing, resources=nothing, verbs=nothing)
        o = new()
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("apiGroups"), apiGroups)
        setfield!(o, Symbol("apiGroups"), apiGroups)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("clusterScope"), clusterScope)
        setfield!(o, Symbol("clusterScope"), clusterScope)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("namespaces"), namespaces)
        setfield!(o, Symbol("namespaces"), namespaces)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("resources"), resources)
        setfield!(o, Symbol("resources"), resources)
        validate_property(IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule, Symbol("verbs"), verbs)
        setfield!(o, Symbol("verbs"), verbs)
        o
    end
end # type IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule

const _property_map_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule = Dict{Symbol,Symbol}(Symbol("apiGroups")=>Symbol("apiGroups"), Symbol("clusterScope")=>Symbol("clusterScope"), Symbol("namespaces")=>Symbol("namespaces"), Symbol("resources")=>Symbol("resources"), Symbol("verbs")=>Symbol("verbs"))
const _property_types_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule = Dict{Symbol,String}(Symbol("apiGroups")=>"Vector{String}", Symbol("clusterScope")=>"Bool", Symbol("namespaces")=>"Vector{String}", Symbol("resources")=>"Vector{String}", Symbol("verbs")=>"Vector{String}")
Base.propertynames(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }) = collect(keys(_property_map_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule))
Swagger.property_type(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule[name]))}
Swagger.field_name(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }, property_name::Symbol) =  _property_map_IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule[property_name]

function check_required(o::IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule)
    (getproperty(o, Symbol("apiGroups")) === nothing) && (return false)
    (getproperty(o, Symbol("resources")) === nothing) && (return false)
    (getproperty(o, Symbol("verbs")) === nothing) && (return false)
    true
end

function validate_property(::Type{ IoK8sApiFlowcontrolV1alpha1ResourcePolicyRule }, name::Symbol, val)
end


text:
# Fisher metric
struct FisherMetric{N} <: AbstractMetric
    I::Function
    dom::Function
end

# Constructors
FisherMetricManifold(N::Int,I::Function,dom::Function=p->true) = MetricManifold(
    Euclidean(N),
    FisherMetric{N}(I,dom)
)
FisherMetricManifold(d,dom::Function=p->true) = MetricManifold(
    Euclidean(length(fieldnames(d))),
    FisherMetric{length(fieldnames(d))}(p -> fim(eval(Expr(:call,d,p...))),dom)
)

# Methods
function local_metric(M::MetricManifold{ℝ,<:AbstractManifold,<:FisherMetric},p,B::InducedBasis)
    metric(M).I(p)
end
function check_point(M::MetricManifold{ℝ,<:AbstractManifold,<:FisherMetric{N}}, p; kwargs...) where {N}
    (size(p)) == (N,) || return DomainError(size(p),"The size of $p is not $((N,)).")
    metric(M).dom(p) || return DomainError(p,"p is not a valid point in the metric domain")
    return nothing
end
function log!(M::MetricManifold{ℝ,<:AbstractManifold,<:FisherMetric},X,p,q)
    # Use shooting
    X .= optimize(X -> norm(q - geodesic(M,p,X)(1.0)), q - p * norm(exp(M,p,X) - p) / norm(q - p)).minimizer
end
function inner(M::MetricManifold{ℝ,<:AbstractManifold,<:FisherMetric},p,X,Y)
    X' * metric(M).I(p) * Y
end

# Default basis (to make things easy)
function default_basis(M,p)
    induced_basis(M,Manifolds.RetractionAtlas(),p,TangentSpace)
end

max_stars_repo_name: ap-browning/FisherMetrics

max_issues_repo_name: ap-browning/FisherMetrics

max_forks_repo_name: ap-browning/FisherMetrics


text:
ENV["ZYGOTE_TYPED"] = true
include("../utils.jl")

using Zygote, NNlib

const loop_iters = 500

function loop(x = 2, n = loop_iters)
  r = one(x)/one(x)
  while n > 0
    r *= sin(x)
    n -= 1
  end
  return r
end

# Ref prevents LLVM from constant folding everything away
@result loop Baseline loop(2, n = $(Ref(loop_iters))[])

@result loop Zygote gradient(loop, 2)

logreg(W, b, x) = sum(σ.(W * x .+ b))

W = randn(10,28^2)
b = randn(10)
x = rand(28^2)

@result logreg Baseline logreg($W, $b, $x)

@result logreg Zygote gradient(logreg, $W, $b, $x)