text:
using YAML

include("extract_sprites_meta.jl")

const atlases = Dict{String, Dict{String, SpriteHolder}}()

const defaultBlankColor = (0.0, 0.0, 0.0, 0.0)
const defaultBlackColor = (0.0, 0.0, 0.0, 1.0)
const defaultWhiteColor = (1.0, 1.0, 1.0, 1.0)

# For missing resources
# Provides a simple image indicating the resource is missing
const fileNotFoundSurface = Assets.missingImage
const fileNotFoundWidth = width(fileNotFoundSurface)
const fileNotFoundHeight = height(fileNotFoundSurface)
const fileNotFoundSpriteHolder = SpriteHolder(
    Sprite(
        0,
        0,

        fileNotFoundWidth,
        fileNotFoundHeight,

        0,
        0,
        fileNotFoundWidth,
        fileNotFoundHeight,

        fileNotFoundSurface,
        ""
    ),

    time(),
    ""
)

function getResourceFromZip(filename::String, resource::String, atlas::String, allowRetry=true)
    fileCache = get(spriteZipCache, filename, nothing)

    if fileCache !== nothing
        return get(fileCache, "$atlas/$resource", nothing)

    else
        if allowRetry
            cacheZipContent(filename)

            return getResourceFromZip(filename, resource, atlas, false)

        else
            return nothing
        end
    end
end

const metaYamlResourceNotFound = Dict{String, Any}()

function getSpriteSurface(resource::String, filename::String, atlas="Gameplay")
    # If we have a filename use that, otherwise look for the resource
    filename = isempty(filename) ? something(findExternalSprite(resource), "") : filename

    if hasExt(filename, ".png")
        return open(Cairo.read_from_png, filename), filename

    elseif hasExt(filename, ".zip")
        surface = getResourceFromZip(filename, resource * ".png", atlas)

        return something(surface, fileNotFoundSurface), filename
    end

    return fileNotFoundSurface, filename
end

function getSpriteMeta(resource::String, filename::String, atlas="Gameplay")
    if !get(config, "load_image_meta_yaml", false)
        return false, false
    end

    # If we have a filename use that, otherwise look for the resource
    filename = isempty(filename) ? something(findExternalSprite(resource), "") : filename
    yamlFilename = splitext(filename)[1] * ".meta.yaml"

    if hasExt(filename, ".png")
        if isfile(yamlFilename)
            try

                return true, open(YAML.load, yamlFilename)

            catch
                return false, "Invalid YAML data"
            end

        else
            return false, false
        end

    elseif hasExt(filename, ".zip")
        data = getResourceFromZip(resource, yamlFilename, atlas)

        return data !== nothing, something(data, metaYamlResourceNotFound)
    end

    return fileNotFoundSurface, filename
end

function addSprite!(resource::String, filename::String=""; atlas="Gameplay", force=false)::Sprite
    spriteAtlas = getAtlas(atlas)
    spriteHolder = get(spriteAtlas, resource, nothing)

    if force || spriteHolder === nothing || spriteHolder.sprite.surface === fileNotFoundSurface || spriteHolder.sprite.width == 0 || spriteHolder.sprite.height == 0 || mtime(spriteHolder.path) > spriteHolder.readtime || mtime(spriteHolder.path) == 0
        surface, filename = getSpriteSurface(resource, filename)
        success, meta = getSpriteMeta(resource, filename, atlas)
        hasMeta = isa(meta, Dict)

        if !success && isa(meta, String)
            println(Base.stderr, "Problem with YAML file for resource '$resource', at file '$filename' - '$meta'")
        end

        imWidth = Int(width(surface))
        imHeight = Int(height(surface))

        realWidth = hasMeta ? get(meta, "Width", imWidth) : imWidth
        realHeight = hasMeta ? get(meta, "Height", imHeight) : imHeight

        offsetX = hasMeta ? -get(meta, "X", 0) : 0
        offsetY = hasMeta ? -get(meta, "Y", 0) : 0

        spriteAtlas[resource] = SpriteHolder(
            Sprite(
                0,
                0,

                imWidth,
                imHeight,

                offsetX,
                offsetY,
                realWidth,
                realHeight,

                surface,
                filename
            ),

            time(),
            filename
        )
    end

    return spriteAtlas[resource].sprite
end

function getAtlas(atlas="Gameplay")::Dict{String, SpriteHolder}
    return get!(Dict{String, SpriteHolder}, atlases, atlas)
end

function getAtlasResourceNames(atlas="Gameplay")
    atlas = getAtlas(atlas)

    return collect(keys(atlas))
end

# Get a resource if loaded, otherwise returns the "missing image" sprite
function getSprite(resource::String, atlas="Gameplay")::Sprite
    spriteAtlas = getAtlas(atlas)
    spriteHolder::SpriteHolder = get(spriteAtlas, resource, fileNotFoundSpriteHolder)

    return spriteHolder.sprite
end

function loadAllExternalSprites!(force=false)
    externalSprites = findAllExternalSprites()

    for (texture, raw, atlas) in externalSprites
        addSprite!(fixTexturePath(texture), raw, atlas=atlas, force=force)
    end
end

# Always forced loads, we know they have changed
function loadChangedExternalSprites!()
    externalSprites = findChangedExternalSprites()
    unloadedZips = Set{String}()

    for (texture, raw, atlas) in externalSprites
        if !(raw in unloadedZips) && hasExt(raw, ".zip")
            uncacheZipContent(raw)
        end

        addSprite!(fixTexturePath(texture), raw, atlas=atlas, force=true)
    end
end

const fixedTexturePathsCache = Dict{String, String}()

# Remove .png and convert to unix paths
function fixTexturePath(texture::String)::String
    return get!(fixedTexturePathsCache, texture) do
        replace(last(texture, 4) == ".png" ? texture[1:end - 4] : texture, '\\' => '/')
    end
end

function getTextureSprite(texture::String, atlas="Gameplay")::Sprite
    return getSprite(fixTexturePath(texture), atlas)
end

mutable struct DrawingAlpha
    alpha::Float64
end

const drawingAlpha = DrawingAlpha(1.0)

function setGlobalAlpha!(alpha=1.0)
    drawingAlpha.alpha = alpha
end

getGlobalAlpha() = drawingAlpha.alpha

# This is probably one of the stupid ways to do it.
# Look into this if it turns out being slow

# Image based #

function drawImage(ctx::Cairo.CairoContext, surface::Cairo.CairoSurface, x, y, quadX, quadY, width, height; alpha=nothing, tint=nothing)
    if ctx.ptr != C_NULL
        if tint !== nothing
            tinted = CairoImageSurface(zeros(UInt32, height, width), Cairo.FORMAT_ARGB32)
            tintedCtx = getSurfaceContext(tinted)

            # Don't pass tint or alpha
            drawImage(tintedCtx, surface, 0, 0, quadX, quadY, width, height)
            tinted = tintSurfaceMultiplicative(tinted, tint)

            drawImage(ctx, tinted, x, y, 0, 0, width, height, alpha=alpha)

            deleteSurface(tinted)

        else
            drawingAlpha = something(alpha, getGlobalAlpha())

            Cairo.save(ctx)

            rectangle(ctx, x, y, width, height)
            clip(ctx)

            set_source_surface(ctx, surface, x - quadX, y - quadY)
            pattern_set_filter(get_source(ctx), Cairo.FILTER_NEAREST)

            paint_with_alpha(ctx, drawingAlpha)

            restore(ctx)
        end
    end
end

function drawImage(ctx::Cairo.CairoContext, surface::Cairo.CairoSurface, x, y; alpha=nothing, tint=nothing)
    drawImage(ctx, surface, x, y, 0, 0, Int(surface.width), Int(surface.height); alpha=alpha, tint=tint)
end

function drawImage(canvas::Gtk.GtkCanvas, name::String, x, y, quadX, quadY, width, height; alpha=nothing, tint=nothing, atlas="Gameplay")
    drawImage(canvas, getSprite(name, atlas), x, y, quadX, quadY, width, height, alpha=alpha, tint=tint)
end

function drawImage(ctx::Cairo.CairoContext, name::String, x, y, quadX, quadY, width, height; alpha=nothing, tint=nothing, atlas="Gameplay")
    drawImage(ctx, getSprite(name, atlas), x, y, quadX, quadY, width, height, alpha=alpha, tint=tint)
end

function drawImage(canvas::Gtk.GtkCanvas, sprite::Sprite, x, y, quadX, quadY, width, height; alpha=nothing, tint=nothing)
    drawImage(Gtk.getgc(canvas), sprite.surface, x, y, quadX, quadY, width, height, alpha=alpha, tint=tint)
end

function drawImage(ctx::Cairo.CairoContext, sprite::Sprite, x, y, quadX=0, quadY=0, width=sprite.width, height=sprite.height; alpha=nothing, tint=nothing)
    if ctx.ptr != C_NULL
        if tint !== nothing
            tinted = CairoImageSurface(zeros(UInt32, height, width), Cairo.FORMAT_ARGB32)
            tintedCtx = getSurfaceContext(tinted)

            # Don't pass tint or alpha
            drawImage(tintedCtx, sprite, 0, 0, quadX, quadY, width, height)
            tinted = tintSurfaceMultiplicative(tinted, tint)

            drawImage(ctx, tinted, x, y, 0, 0, width, height, alpha=alpha)

            deleteSurface(tinted)

        else
            drawingAlpha = something(alpha, getGlobalAlpha())

            Cairo.save(ctx)

            rectangle(ctx, x, y, width, height)
            clip(ctx)

            set_source_surface(ctx, sprite.surface, x - quadX - sprite.x, y - quadY - sprite.y)
            pattern_set_filter(get_source(ctx), Cairo.FILTER_NEAREST)

            paint_with_alpha(ctx, drawingAlpha)

            restore(ctx)
        end
    end
end

function drawImage(canvas::Gtk.GtkCanvas, name::String, x, y; alpha=nothing, tint=nothing, atlas="Gameplay")
    drawImage(canvas, getSprite(name, atlas), x, y, alpha=alpha, tint=tint)
end

function drawImage(ctx::Cairo.CairoContext, name::String, x, y; alpha=nothing, tint=nothing, atlas="Gameplay")
    drawImage(ctx, getSprite(name, atlas), x, y, alpha=alpha, tint=tint)
end

function drawImage(canvas::Gtk.GtkCanvas, sprite::Sprite, x, y; alpha=nothing, tint=nothing)
    drawImage(Gtk.getgc(canvas), sprite, x, y, alpha=alpha, tint=tint)
end

# TODO - Only scale 1/-1 works as intended
function drawSprite(ctx::Cairo.CairoContext, texture::String, x, y; jx=0.5, jy=0.5, sx=1, sy=1, rot=0, alpha=nothing, tint=nothing, atlas="Gameplay")
    if ctx.ptr != C_NULL
        sprite = getTextureSprite(texture, atlas)

        width, height = sprite.realWidth, sprite.realHeight
        drawX, drawY = floor(Int, x - width * jx - sprite.offsetX) , floor(Int, y - height * jy - sprite.offsetY)

        Cairo.save(ctx)

        scale(ctx, sx, sy)
        translate(ctx, sx > 0 ? 0 : -2 * x, sy > 0 ? 0 : -2 * y)
        translate(ctx, drawX, drawY)
        rotate(ctx, rot)
        drawImage(ctx, sprite, 0, 0, alpha=alpha, tint=tint)

        restore(ctx)
    end
end

# Create image from matrix
# Bitmap colors are 0 index based (false = 0, true = 1)
function matrixToSurface(bitmap::Array{T, 2}, colors::Array{rgbaTupleType, 1}=rgbaTupleType[defaultBlackColor, defaultWhiteColor]) where T <: Integer
    height, width = size(bitmap)
    pixels = width * height

    data = Array{Cairo.ARGB32, 1}(undef, pixels)

    ptr = 1
    for x in 1:width, y in 1:height
        color = colors[bitmap[y, x] + 1]
        data[ptr] = Cairo.ARGB32(color[1], color[2], color[3], color[4])

        ptr += 1
    end

    data = permutedims(reshape(data, (height, width)))

    return CairoImageSurface(data)
end

# Image tinting, super inefficient

function tintSurfaceAdditive(surface::Cairo.CairoSurface, color::rgbaTupleType, width=width(surface), height=height(surface))
    for x in 1:floor(Int, width), y in 1:floor(Int, height)
        pixelColor = argb32ToRGBATuple(surface.data[x, y]) ./ 255
        surface.data[x, y] = normRGBA32Tuple2ARGB32(min.(1.0, pixelColor .+ color))
    end

    return surface
end

function tintSurfaceMultiplicative(surface::Cairo.CairoSurface, color::rgbaTupleType, width=width(surface), height=height(surface))
    for x in 1:floor(Int, width), y in 1:floor(Int, height)
        pixelColor = argb32ToRGBATuple(surface.data[x, y]) ./ 255
        surface.data[x, y] = normRGBA32Tuple2ARGB32(pixelColor .* color)
    end

    return surface
end

# Shapes #

function setSourceColor(ctx, c::rgbTupleType)
    set_source_rgb(ctx, c[1], c[2], c[3])
end

function setSourceColor(ctx, c::rgbaTupleType)
    set_source_rgba(ctx, c[1], c[2], c[3], c[4])
end

function paintSurface(ctx::Cairo.CairoContext, c=defaultWhiteColor)
    if ctx.ptr != C_NULL
        setSourceColor(ctx, c)
        paint(ctx)
    end
end

paintSurface(surface::Cairo.CairoSurface, c=defaultWhiteColor) = paintSurface(getSurfaceContext(surface), c)

function clearArea(ctx::Cairo.CairoContext, x, y, width, height)
    if ctx.ptr != C_NULL
        Cairo.save(ctx)

        set_operator(ctx, Cairo.OPERATOR_CLEAR)
        rectangle(ctx, x, y, width, height)
        fill(ctx)

        restore(ctx)
    end
end

clearArea(surface::Cairo.CairoSurface, x, y, width, height) = clearArea(getSurfaceContext(surface), x, y, width, height)

function clearSurface(ctx::Cairo.CairoContext)
    return clearArea(ctx, 0, 0, width(ctx), height(ctx))
end

clearSurface(surface::Cairo.CairoSurface) = clearSurface(getSurfaceContext(surface))

# With border
function drawRectangle(ctx::Cairo.CairoContext, x, y, w, h, fc=defaultBlackColor, rc=defaultBlankColor)
    if ctx.ptr != C_NULL
        Cairo.save(ctx)

        setSourceColor(ctx, fc)
        rectangle(ctx, x, y, w, h)
        fill_preserve(ctx)
        setSourceColor(ctx, rc)
        stroke(ctx)

        restore(ctx)
    end
end

function drawRectangle(ctx::Cairo.CairoContext, rect::Rectangle, fc=defaultBlackColor, rc=defaultBlankColor)
    return drawRectangle(ctx, rect.x, rect.y, rect.w, rect.h, fc, rc)
end

function drawLines(ctx::Cairo.CairoContext, nodes, sc=defaultBlackColor; thickness=2, filled=false, fc=sc)
    if ctx.ptr != C_NULL
        if length(nodes) < 2
            return
        end

        Cairo.save(ctx)

        set_antialias(ctx, thickness)
        set_line_width(ctx, thickness)

        move_to(ctx, nodes[1][1], nodes[1][2])

        for node in nodes[2:end]
            line_to(ctx, node[1], node[2])
        end

        if filled
            setSourceColor(ctx, fc)
            fill_preserve(ctx)
        end

        setSourceColor(ctx, sc)
        stroke(ctx)
        restore(ctx)
    end
end

function drawSimpleCurve(ctx::Cairo.CairoContext, curve::SimpleCurve, c=defaultBlackColor; resolution=25, thickness=2)
    points = [getPoint(curve, i / resolution) for i in 0:resolution]

    drawLines(ctx, points, c, thickness=thickness)
end

function drawArc(ctx::Cairo.CairoContext, x, y, r, alpha, beta, c=defaultBlackColor)
    if ctx.ptr != C_NULL
        Cairo.save(ctx)

        setSourceColor(ctx, c)
        arc(ctx, x, y, r, alpha, beta)
        stroke(ctx)

        restore(ctx)
    end
end

function drawCircle(ctx::Cairo.CairoContext, x, y, r, c=defaultBlackColor)
    drawArc(ctx, x, y, r, 0, 2 * pi, c)
end

function drawArrow(ctx::Cairo.CairoContext, x1, y1, x2, y2, c=defaultBlackColor; filled=true, headLength=15, headTheta=pi / 6)
    theta = atan(y2 - y1, x2 - x1)
    len = sqrt((x1 - x2)^2 + (y1 - y2)^2)

    shaftLen = len - headLength
    sideLen = tan(headTheta) * headLength

    shaftX = shaftLen * cos(theta) + x1
    shaftY = shaftLen * sin(theta) + y1

    dx = sideLen * cos(theta + pi / 2)
    dy = sideLen * sin(theta + pi / 2)

    nodes = Tuple{Float64, Float64}[
        (x1, y1),
        (shaftX, shaftY),
        (shaftX + dx, shaftY + dy),
        (x2, y2),
        (shaftX - dx, shaftY - dy),
        (shaftX, shaftY),
    ]

    return drawLines(ctx, nodes, c, filled=filled)
end

max_stars_repo_name: asdfsdfadfasfasd/Ahorn

max_issues_repo_name: asdfsdfadfasfasd/Ahorn

max_forks_repo_name: asdfsdfadfasfasd/Ahorn


text:
#!/usr/bin/env julia

"""
    function detect_polyA_in_a_string(fq_seq::String,minimum_polyA_length::Int64,maximum_non_A_symbols::Int64;maximum_search_fragment_length::Int64=50,debug=false)

    Checks if a string of a read contains polyA strechA starting seach from the 3'.
    Once a g search once a strech is found. A first and last symbol of a strech must be A.

    # Arguments
    - `fq_seq::String`: string for testing.
    - `minimum_polyA_length::Int64`: minimum length of polyA strech.
    - `maximum_non_A_symbols`::Int64: maximum numer of nonA symbols in polyA strech.
    - `re:Regex`: regex to search for polyA.

    # Keyword Arguments
    - `maximum_search_fragment_length::Int64=50`: fragment length from 3' end.
"""
function detect_polyA_in_a_string(
    fq_seq::String,
    minimum_polyA_length::Int64,
    maximum_non_A_symbols::Int64,
    re::Regex;
    maximum_search_fragment_length::Int64=50,
    debug=false
    )::Bool

    if maximum_non_A_symbols == 1
        return detect_polyA_in_a_string(fq_seq,minimum_polyA_length,maximum_search_fragment_length,re,debug)
    end
    has_poly_a = false
    fq_length = length(fq_seq)
    window_position_from_3_end = minimum_polyA_length
    startini = minimum_polyA_length

    while (window_position_from_3_end < maximum_search_fragment_length) && (startini  != fq_length)
        window_position_from_3_end += 1
        start = fq_length - startini + 1
        # skip searching for polyA far from the  right end
        substring = String(fq_seq[start:start+minimum_polyA_length-1])

        if substring[1] == 'A' && substring[minimum_polyA_length] == 'A' #first and last symbol of a polyA strech must be A
            ct = 0
            for symb in substring
                if symb != 'A'
                    ct += 1
                end
            end

            if ct <= maximum_non_A_symbols
                has_poly_a = true

                if debug
                    println("substring having at most than $maximum_non_A_symbols non A symbols")
                    println(substring)
                end
                break
            end
        end
        startini += 1
    end
    return has_poly_a
end

function detect_polyA_in_a_string(
    fq_seq::String,
    minimum_polyA_length::Int64,
    maximum_search_fragment_length::Int64,
    re::Regex,
    debug::Bool
    )::Bool

    fq_seq = uppercase(fq_seq)
    fq_length = length(fq_seq)

    if fq_length <= maximum_search_fragment_length
        if fq_length < minimum_polyA_length
            if debug
                println(STDERR,"ERROR! Sequence too short to detect polyA")
            end
            return false
        else
            maximum_search_fragment_length = fq_length - 1
        end
    end
    fq_seq = reverse(fq_seq[end-maximum_search_fragment_length:end])

    for m in eachmatch(re,fq_seq, true)
        if m[1] == nothing
            len1 = 0
        else
            len1 = length(m[1])
        end

        if m[3] == nothing
            len3 = 0
        else
            len3 = length(m[3])
        end

        len2 = length(m[2])

        if len2 >= minimum_polyA_length
            return true
        elseif len1 + len2 >= minimum_polyA_length || len3 + len2 >= minimum_polyA_length
            return true
        end
    end
    return false
end

"""
    extend_poly_A(fq_seq1::FASTQ.Record,fq_seq2::FASTQ.Record)

    Tries to merge and extend forward read.

    # Arguments
    - `fq_seq1::FASTQ.Record`: string of forward read.
    - `fq_seq2::FASTQ.Record`: string of reverse read reverse complement.
"""
function extend_poly_A(
    fq_seq1::FASTQ.Record,
    fq_seq2::FASTQ.Record)::FASTQ.Record

    return fq_seq1
end

"""
    check_polyA_prefixes(fqo_trimmed::FASTQ.Record,prefixes::FMIndexes.FMIndex{7,},maximum_distance_with_prefix_database::Int64,minimum_not_polyA::Int64)

    Check if the read is in pefixes list of naturall polyA sreches.

    # Arguments
    - `fqo_trimmed::FASTQ.Record,`: fastq record.
    - `prefixes::FMIndexes.FMIndex{7,}`: array of natural prefixes.
    - `maximum_distance_with_prefix_database::Int64`: maximum Levenstain distance.
    - `minimum_not_polyA::Int64`: Minimum of polyA lenght.
"""
function check_polyA_prefixes(
    fqo_trimmed::FASTQ.Record,
    prefixes::FMIndexes.FMIndex{7,},
    maximum_distance_with_prefix_database::Int64,
    minimum_not_polyA::Int64
    )::Bool

    read_substring_for_check=String(FASTQ.sequence(fqo_trimmed))[(end-minimum_not_polyA+1):end]

    has_no_match=(count(read_substring_for_check,prefixes)<1)
    return  has_no_match
end

"""
    trim_polyA_3end(seq::String, minimum_poly_A_between::Int64)

    Function that trimmes 3' end of polyA having sequence from non A symbols that might originate fue to reamins of adapters or sequencing artefacts

    # Arguments
    - `seq::String`: read sequence in String
    - `minimum_poly_A_between::Int64`: minimum length of polyA strech that might occure between non A symbols forming a fragment to be trimmed off (starting from the very 3' end)
"""
function trim_polyA_3end(seq::String,
    minimum_poly_A_between::Int64
    )::Tuple{String,Int64}

    length_sequence = length(seq)
    i = length(seq)
    trimming_position = 0
    polyA_strech_length = 0

    while (polyA_strech_length <= minimum_poly_A_between) && (i>0)
        symbol=seq[i]

        if symbol != 'A'
            trimming_position = i
            polyA_strech_length = 0
        else
            polyA_strech_length += 1
        end
        i -= 1
    end

    if trimming_position == 0
        out = (seq,0)
    else
        out = (seq[1:trimming_position-1],length_sequence-trimming_position+1)
    end
    return out
end

"""
    count_nona(seq::String)

    Counts not A in a string.
"""
function count_nona(seq::String)::Int64

    ct = 0

    for s in seq
        if s != 'A'
            ct += 1
        end
    end
    return ct
end

"""
    first_nona(seq::String)

    Finds first not A from the end (3')
"""
function first_nona(seq::String)::Int64

    ct = 0
    l = length(seq)

    for i in 1:l
        pos = l - ct

        if seq[pos] != 'A'
            break
        end
        ct += 1
    end
    return ct
end

"""
    Function that trimmes 3 polyA tail of a given fastq entry and outputs trimmed read with number of polyA indicated in read's name

    # Arguments
    - `minimum_not_polyA::Int64`: minimum length of not polyA strech in a read
    - `minimum_polyA_length::Int64`: minimum length of polyA strech
    - `minimum_poly_A_between::Int64`: minimum length of polyA strech that might occure between non A symbols forming a fragment to be trimmed off (starting from the very 3' end)

    # Keywords arguments
    - `debug::Bool`: turns out additional output
    - `window_length::Int64`: windows size for polyA trimming
    - `max_nonA_in_window::Int64`: tollerated nonA number per window

    # The algorithms works like
    1. Trimmes offs remains of adapters at the 3'
    2. Search for a window wich exceeds nonA limit and has no other high A content window towards 5' end
    3. Ajusts the A position in the window
"""
function trim_polyA_from_fastq_record(fq::FASTQ.Record,
    minimum_not_polyA::Int64,
    minimum_polyA_length::Int64,
    minimum_poly_A_between::Int64;
    debug=false,
    window_length=10,
    max_nonA_in_window=1,
    limit_for_3end_trimming=20
    )::Tuple{FASTQ.Record,Bool,Bool}

    has_proper_polyA = false
    polyA_detected = true
    seq_ini = String(FASTQ.sequence(fq),)
    seq,trimmed_to_polyA = trim_polyA_3end(seq_ini,minimum_poly_A_between)
    seq_len = length(seq)
    name = FASTQ.identifier(fq)

    if FASTQ.hasdescription(fq)
        description = FASTQ.description(fq)
    else
        description = nothing
    end

    quality = FASTQ.quality(fq)[1:seq_len]
    i = seq_len
    polyA_start = 0
    first_nona_position = 0

    if debug
        println("-----------------------------------------------------------------------------------------------------------------------------\nSequence for trimming:\n$seq_ini")
        println("Sequence with chopped 3'end adapter remains:\n$seq")

        if !(polyA_detected)
            println("WARNING! this sequence will be concidered as not having polyA tail due to trimming of 3end adapter remains too much (> $limit_for_3end_trimming bp)")
        end
    end

    # search for a polyA start with a sliding window from the 3' end
    #Analysis limit to the cases where we meet a new non A symbol (going from 3' end)
    initial_position = 0 # position of polyA start based on windows.
    first_nona_position = 0

    while i > 0 & initial_position == 0
        s = seq[i]

        if s != 'A'
            fragment_start = i
            fragment_end = i + window_length - 1
            # println("fragment_start:fragment_end-ini seq_len $fragment_start $fragment_end $seq_len")
            # println(seq)

            if fragment_end > seq_len #handle the cases when there are non A at the 3' end
                fragment_end_ajusted = seq_len
                fragment_start_ajusted = i - (fragment_end-seq_len)
                fragment_start = fragment_start_ajusted
                fragment_end = fragment_end_ajusted
            end

            if fragment_start < 1
                fragment_start = 1
            end
            # println("fragment_start:fragment_end $fragment_start $fragment_end")
            fragment = seq[fragment_start:fragment_end]
            number_of_nonA = count_nona(fragment)

            if debug
                println("Fragment with non A, starting position $i, with $number_of_nonA nonA symbols")
                println(fragment)
            end

            if max_nonA_in_window < number_of_nonA
                #test if a window towards 5' also contains non A sequences
                to_5end_frag_end=fragment_start - 1
                to_5end_frag_start=to_5end_frag_end - window_length + 1

                if to_5end_frag_start < 1
                     to_5end_frag_start = 1
                end

                if to_5end_frag_start == 1 || (count_nona(seq[to_5end_frag_start:to_5end_frag_end]) > max_nonA_in_window)
                    initial_position = fragment_end
                    first_nona_position = initial_position - first_nona(seq[1:initial_position])

                    if debug
                        println("Initial non A sequence:")
                        println(seq[1:initial_position])
                        println("Adjusted non A sequence:")
                        println(seq[1:first_nona_position])
                    end
                    break
                end
            end
        end
        i -= 1
    end

    polyA_start = first_nona_position + 1
    polyA_length = seq_len-polyA_start + 1

    if (polyA_start > 0) & (polyA_start > minimum_not_polyA)
        not_poly_a_seq = seq[1:polyA_start-1]
        nname = string(polyA_length) * ":A:" * name
        nquality = quality[1:polyA_start-1]

        if description != nothing
            fqo = FASTQ.Record(nname, description, not_poly_a_seq, nquality)
        else
            fqo = FASTQ.Record(nname, not_poly_a_seq, nquality)
        end

        has_proper_polyA = true
    end

    if (trimmed_to_polyA > limit_for_3end_trimming) || (polyA_length < minimum_polyA_length)
        polyA_detected = false
    end

    if has_proper_polyA
        return fqo,has_proper_polyA,polyA_detected
    else
        return fq,has_proper_polyA,polyA_detected
    end
end

"""
    get_polyA_prefixes(fasta_record::FASTA.Record,minimum_not_polyA::Int64,minimum_polyA_length::Int64)

    Returns minimum kmer streches of not polyA ins a supplied transcript sequences

    # Arguments
    - `sequence::FASTA.Record`: string with polyA streches (type - BioSequences reference seq)
    - `minimum_not_polyA::Int64`: minimum length of not polyA strech
    - `minimum_polyA_length::Int64`: minimum length of polyA strech
"""
function get_polyA_prefixes(fasta_record::FASTA.Record,
    minimum_not_polyA::Int64,
    minimum_polyA_length::Int64,
    re::Regex)::Array{String,1}

    output = Array{String,1}()

    for m in eachmatch(re, String(sequence(fasta_record)))
        push!(output,m[1])
    end
    return output
end

"""
    get_polyA_prefixes(fasta_record::FASTA.Record,minimum_not_polyA::Int64,minimum_polyA_length::Int64,counter::SharedArray)

    Returns minimum kmer streches of not polyA ins a supplied transcript sequences

    # Arguments
    - `sequence::FASTA.Record`: string with polyA streches (type - BioSequences reference seq)
    - `minimum_not_polyA::Int64`: minimum length of not polyA strech
    - `minimum_polyA_length::Int64`: minimum length of polyA strech
    - `counter`: shared array for tracking processed reads - length matches julia workers
"""
function get_polyA_prefixes(fasta_record::FASTA.Record,
    minimum_not_polyA::Int64,
    minimum_polyA_length::Int64,
    re::Regex,
    counter::SharedArray)::Array{String,1}

    output = get_polyA_prefixes(fasta_record,minimum_not_polyA,minimum_polyA_length,re)
    counter[(myid()-1)] += 1
    total_jobs = sum(counter)

    if mod(total_jobs,10000) == 0
      print(STDERR,"Number of parsed transcripts = $total_jobs\r")
    end
    return output
end

"""
    trim_polyA_from_fastq_pair(fastq1::FASTQ.Record,fastq2::FASTQ.Record,prefixes::FMIndexes.FMIndex{7,},minimum_not_polyA::Int64,minimum_polyA_length::Int64,maximum_non_A_symbols::Int64,maximum_distance_with_prefix_database::Int64,minimum_poly_A_between::Int64)

    Finds and trims polyA having reads from a pair of fastq records

    # Arguments
    - `fastq1::FASTQ.Record`: FASTQ record
    - `fastq2::FASTQ.Record`: FASTQ record
    - `prefixes::FMIndexes.FMIndex{7,}`: array of prefixes of natural polyA
    - `minimum_not_polyA::Int64`: minimum length of not polyA strech in a read
    - `minimum_polyA_length::Int64`: minimum length of polyA strech
    - `maximum_non_A_symbols::Int64`: maximum numer of nonA symbols in polyA strech
    - `maximum_distance_with_prefix_database::Int64`: allowed Levenshtein distance between prefix of natural polyA in transcripts and the read
    - `minimum_poly_A_between::Int64`: minimum length of polyA strech that might occure between non A symbols forming a fragment to be trimmed off (starting from the very 3' end)
"""
function trim_polyA_from_fastq_pair(
    fastq1::FASTQ.Record,
    fastq2::FASTQ.Record,
    prefixes::FMIndexes.FMIndex{7,},
    minimum_not_polyA::Int64,
    minimum_polyA_length::Int64,
    maximum_non_A_symbols::Int64,
    maximum_distance_with_prefix_database::Int64,
    minimum_poly_A_between::Int64,
    re::Regex;
    debug_id="None",
    debug=false
    )

    polyA_detected = false
    has_proper_polyA = false
    seq_for_read = String(FASTQ.sequence(fastq1))
    revseq_rev_read = String(reverse_complement!(FASTQ.sequence(fastq2)))

    if debug
        println("for_read:", "\n", fastq1, "\n")
        println("rev_read:", "\n", fastq2, "\n")
        println("seq_for_read", seq_for_read,"\n")
        println("revseq_rev_read", revseq_rev_read,"\n")
        println("Is polyA in seq_for_read: ", detect_polyA_in_a_string(seq_for_read,minimum_polyA_length,maximum_non_A_symbols,re))
        println("Is polyA in revseq_rev_read: ", detect_polyA_in_a_string(revseq_rev_read,minimum_polyA_length,maximum_non_A_symbols,re))
        println("used parameter: minimum_polyA_length $minimum_polyA_length, maximum_non_A_symbols $maximum_non_A_symbols")
    end

    if  detect_polyA_in_a_string(seq_for_read,minimum_polyA_length,maximum_non_A_symbols,re) &&
        detect_polyA_in_a_string(revseq_rev_read,minimum_polyA_length,maximum_non_A_symbols,re)
            #tries sto esxtend polyA
            cosensus_extended_fwd_read=extend_poly_A(fastq1,fastq2)
            fqo_trimmed, has_proper_polyA, polyA_detected = trim_polyA_from_fastq_record(cosensus_extended_fwd_read,
                minimum_not_polyA,
                minimum_polyA_length,
                minimum_poly_A_between,debug=debug)

            if debug
                println("Output of $trim_polyA_from_fastq_record")
                println("Input seq",cosensus_extended_fwd_read)
                println(fqo_trimmed, has_proper_polyA, polyA_detected)
                println(minimum_not_polyA)
                println(minimum_polyA_length)
                println(minimum_poly_A_between)
            end

            #check if the read is in pefixes list of natural polyA streches
            if has_proper_polyA  && (length(prefixes) != 1 )
                has_proper_polyA = check_polyA_prefixes(fqo_trimmed,prefixes,maximum_distance_with_prefix_database,minimum_not_polyA)

                if !has_proper_polyA
                    fqo_trimmed = cosensus_extended_fwd_read
                end
            end
    else
        fqo_trimmed = fastq1
    end
    return fqo_trimmed, has_proper_polyA, polyA_detected
end

function trim_polyA_from_fastq_pair_pararell(
        chunk_for_channels::Array{Tuple{NTuple{4,String},NTuple{4,String}},1},
        fastqo_1_2_s_d::RemoteChannel{Channel{NTuple{4,String}}},
        ct_all::SharedArray{Int64,1},
        ct_pair_without_polyA::SharedArray{Int64,1},
        ct_pair_with_proper_polyA::SharedArray{Int64,1},
        ct_pair_with_discarded_polyA::SharedArray{Int64,1},
        ct_finished::SharedArray{Int64,1},
        ct_output_chunks::SharedArray{Int64,1},
        prefixes::FMIndexes.FMIndex{7,},
        minimum_not_polyA::Int64,
        minimum_polyA_length::Int64,
        maximum_non_A_symbols::Int64,
        maximum_distance_with_prefix_database::Int64,
        minimum_poly_A_between::Int64,
        include_polyA::Bool,
        re::Regex;
        debug_id="None",
        debug=false,
        )

        # mark that task is not finished on this process
        ct_finished[(myid()-1)] = 0
        # initialise buffers for collect of processed read
        fastq1_b = IOBuffer()
        fastq2_b = IOBuffer()
        fastq_s_b = IOBuffer()
        fastq_d_b = IOBuffer()
        #counrer for processed fastq entries per a process
        ct_local = 0

        for chunk_for_channel in chunk_for_channels
            fastq1_4l=chunk_for_channel[1][1]*chunk_for_channel[1][2]*chunk_for_channel[1][3]*chunk_for_channel[1][4]
            fastq2_4l=chunk_for_channel[2][1]*chunk_for_channel[2][2]*chunk_for_channel[2][3]*chunk_for_channel[2][4]
            fastq1 = FASTQ.Record(fastq1_4l)
            fastq2 = FASTQ.Record(fastq2_4l)
            ct_local += 1
            ct_all[(myid()-1)] += 1
            has_proper_polyA = false
            polyA_detected = false

            if debug_id == FASTQ.identifier(fastq1)
                debug = true
            else
                debug = false
            end
            fqo_trimmed, has_proper_polyA, polyA_detected = trim_polyA_from_fastq_pair(
                                                            fastq1,
                                                            fastq2,
                                                            prefixes,
                                                            minimum_not_polyA,
                                                            minimum_polyA_length,
                                                            maximum_non_A_symbols,
                                                            maximum_distance_with_prefix_database,
                                                            minimum_poly_A_between,
                                                            re,
                                                            debug_id=debug_id,
                                                            debug=debug
                                                            )

            if debug
                println("trim_polyA_from_fastq_pair output for $debug_id")
                println("fqo_trimmed:\n $fqo_trimmed \n, has_proper_polyA: $has_proper_polyA, polyA_detected: $polyA_detected")
            end

            if !polyA_detected
                ct_pair_without_polyA[(myid()-1)] += 1
                println(fastq1_b,fastq1)
                println(fastq2_b,fastq2)
            else
                a = 1

                if has_proper_polyA
                    println(fastq_s_b,fqo_trimmed)
                    ct_pair_with_proper_polyA[(myid()-1)] += 1
                    #get reverse complement of the polyA read

                    if include_polyA
                        name = FASTQ.identifier(fqo_trimmed)

                        if FASTQ.hasdescription(fqo_trimmed)
                            description = FASTQ.description(fqo_trimmed)
                        else
                            description = nothing
                        end

                        rev_quality = reverse(FASTQ.quality(fqo_trimmed))
                        rev_seq = reverse_complement!(FASTQ.sequence(fqo_trimmed))

                        if description != nothing
                            fqo_trimmed_rev = FASTQ.Record(name, description, rev_seq, rev_quality)
                        else
                            fqo_trimmed_rev = FASTQ.Record(name, rev_seq, rev_quality)
                        end

                        println(fastq1_b,fqo_trimmed)
                        println(fastq2_b,fqo_trimmed_rev)
                    end
                else
                    println(fastq_d_b,fastq1)
                    ct_pair_with_discarded_polyA[(myid()-1)] += 1
                end
            end
            #write out
        end
    #put to output pipe the final fraction
    fastq1_b_s = String(fastq1_b)
    fastq2_b_s = String(fastq2_b)
    fastq_s_b_s = String(fastq_s_b)
    fastq_d_b_s = String(fastq_d_b)

    if fastq1_b_s != ""   fastq1_b_z=String(transcode(CodecZlib.GzipCompressor,fastq1_b_s)) else fastq1_b_z="" end
    if fastq2_b_s != ""   fastq2_b_z=String(transcode(CodecZlib.GzipCompressor,fastq2_b_s)) else fastq2_b_z="" end
    if fastq_s_b_s != ""   fastq_s_b_z=String(transcode(CodecZlib.GzipCompressor,fastq_s_b_s)) else fastq_s_b_z="" end
    if fastq_d_b_s != ""   fastq_d_b_z=String(transcode(CodecZlib.GzipCompressor,fastq_d_b_s)) else fastq_d_b_z="" end
    ct_output_chunks[(myid()-1)] += 1
    put!(fastqo_1_2_s_d,(fastq1_b_z,fastq2_b_z,fastq_s_b_z,fastq_d_b_z))
    close(fastq1_b)
    close(fastq2_b)
    close(fastq_s_b)
    close(fastq_d_b)
    #close input buffers
    ct_finished[(myid()-1)] = 1

end

ext: jl

max_stars_repo_name: laurabiggins/polyAterminus

max_issues_repo_name: laurabiggins/polyAterminus

max_forks_repo_name: laurabiggins/polyAterminus


text:
# Test statistical measures of traces
using Test
using SAC

@testset "Statistics" begin
    @testset "Window average" begin
        # Trace of constant absolute value of 1 with random polarity
        let delta = 0.1, npts = 1001
            local s = SACtr(rand([-1, 1], npts), delta)
            @test window_average(s, s.b, s.delta*(s.npts-1)) ≈ 1.0
            @test_throws ErrorException window_average(s, s.b-s.delta, s.delta*(s.npts))
            @test_throws ErrorException window_average(s, s.b, -1)
        end
        # Trace of amplitude 1 at 1 sample after zero
        let delta = 1.0, npts = 101, b = -50.0
            local s = SACtr(zeros(npts), delta, b)
            s.t[52] = 1
            @test window_average(s, s.b, 40delta) ≈
                window_average(s, 2delta, 40delta) atol=eps(SAC.SACFloat)
            @test window_average(s, 0, 2delta) ≈ 1/3
            @test window_average(s, delta, 2delta, true) ≈ 1/3
        end
    end

    @testset "STA/LTA" begin
        # Step function from 1 to 2 one sample after 0 time
        let delta = 1.0, npts = 101, b = -50.0, t1 = 2, t2 = 4
            local s = SACtr(ones(npts), delta, b)
            s.t[52:end] .= 2
            @test stalta(s, delta, 10, 10) ≈ 2.0
            @test stalta(s, -delta, 2delta, 0.5, 0.5) ≈ [1, 1, 2, 2]
            @test typeof(stalta(s, 0.1, 1)) == SACtr
            local s′ = stalta(s, t1, t2)
            local n = npts - t1 - t2 - 3
            @test s′.npts == n
            @test time(s′)[argmax(s′.t)] ≈ 1.0
        end
    end
end


text:
"""
    jgraph_morph(g::JGraph, layout_to, lims_to; frames=nothing)

This functions allows to animate the change from one layout to another for a graph.

# Arguments
- `g::JGraph` the `JGraph` to animate.
- `layout_to` a `NetworkLayout.AbstractLayout` or a `Vector{NetworkLayout.AbstractLayout}` that will be applied to nodes and edges
to determine the configuration to reach.
- `lims_to` a tuple or a vector of tuples each holding width and height to scale the positions of nodes and edges in the animation.

# Keywords
- `frames` a range of frames during which the morph will take place.
If not specified the animation will last as long aas the JGraph itself.
- `starting_positions` specifies the starting positions of each node, defaults
to [`get_starting_positions(g)`](@ref)
"""
function jgraph_morph(
    g::JGraph,
    layout_to::Vector,
    lims_to::Vector;
    frames = nothing,
    starting_positions = get_starting_positions(g),
    closed=false,
)

    frames = isnothing(frames) ? g.data.frames : frames
    if closed
        push!(layout_to, g.data.layout)
        push!(lims_to[1], g.data.width)
        push!(lims_to[2], g.data.height)
    end

    frames_stops =
        floor.(Int, range(first(frames) - 1, last(frames), length = length(layout_to) + 1))

    for i in 1:length(layout_to)
        starting_positions = jgraph_morph(
            g,
            layout_to[i],
            lims_to[i],
            frames = (frames_stops[i] + 1):frames_stops[i + 1],
            starting_positions = starting_positions,
        )
    end
    return starting_positions
end

function jgraph_morph(
    g::JGraph,
    layout_to,
    lims_to;
    frames = nothing,
    starting_positions = get_starting_positions(g),
)
    Jnodes = jnodes(g)
    frames = isnothing(frames) ? g.data.frames : frames
    positions = _adapt_to_lims(layout_to(g.data.graph), lims_to...)
    for (node, pos_from, pos_to) in zip(Jnodes, starting_positions, positions)
        act!(node, Action(frames, anim_translate(pos_from, pos_to)))
    end
    return positions
end


text:
using ActiveAppearanceModels
using FaceDatasets
# using PiecewiseAffineTransforms

function smoke_cootes()
    imgs = load_images(FaceDatasets.CootesDataset)
    shapes = load_shapes(FaceDatasets.CootesDataset)
    test_img_idx = 8
    one_out = [1:test_img_idx-1; test_img_idx+1:length(imgs)]
    m = AAModel()
    @time train(m, imgs[one_out], shapes[one_out])
    img = imgs[test_img_idx]
    init_shape = shapes[18]
    @time fitted_shape, fitted_app = fit(m, img, init_shape, 50);
    # triplot(img, fitted_shape, m.wparams.trigs)
end

smoke_cootes()


text:
module CaratheodorySolver

"""
This is a sample function
"""
greet() = print("Hello CaratheodorySolver!!")

end # module


text:
module MeinBeispiel

import ForwardDiff

greet() = println("Hallo du Pfeife")
# Write your package code here.

include("meineFunktion.jl")

export meineFunktion, derivativat_von_meineFunktion

end


text:
using Documenter, ScalarEquations, DocumenterLaTeX, DocumenterTools
push!(LOAD_PATH,"../src/")
makedocs(sitename="ScalarEquations.jl",
authors="C. T. Kelley",
format = Documenter.HTML(
               prettyurls = get(ENV, "CI", nothing) == "true"
           ),
pages = Any[
     "Home" => "index.md",
     "Scalar Equations" => Any[
        "Scalar.md"
      ],
     "Scalar Equations Functions" => Any[
       "functions/functions.md",
       "functions/sptc.md",
      ]
]
)
deploydocs(
     repo="github.com/ctkelley/ScalarEquations.jl.git"
)


text:
using Soundex

# Check basic examples
@assert soundex("Robert") == "R163"
@assert soundex("Rupert") == "R163"
@assert soundex("Rubin") == "R150"

# Verify double-coding gives single output when 's' and 'c' are separated by a 'h' or 'w'
@assert soundex("Ashcraft") == "A261"
@assert soundex("Ashcroft") == "A261"

# Verify double-coding gives double output when 'z' and 'k' are separated by a vowel
@assert soundex("Tymczak") == "T522"

# Verify double-coding for initial 'P' and 'f'
@assert soundex("Pfister") == "P236"

# Verify 0-padding is done correctly
@assert soundex("a") == "A000"
@assert soundex("ab") == "A100"
@assert soundex("abc") == "A120"
@assert soundex("abcd") == "A123"


text:
# Julia wrapper for header: core.h
# Automatically generated using Clang.jl

function sodium_init()
    ccall((:sodium_init, libsodium), Cint, ())
end

function sodium_set_misuse_handler(handler)
    ccall((:sodium_set_misuse_handler, libsodium), Cint, (Ptr{Cvoid},), handler)
end

function sodium_misuse()
    ccall((:sodium_misuse, libsodium), Cvoid, ())
end
# Julia wrapper for header: export.h
# Automatically generated using Clang.jl

# Julia wrapper for header: crypto_hash_sha512.h
# Automatically generated using Clang.jl

function crypto_hash_sha512_statebytes()
    ccall((:crypto_hash_sha512_statebytes, libsodium), Csize_t, ())
end

function crypto_hash_sha512_bytes()
    ccall((:crypto_hash_sha512_bytes, libsodium), Csize_t, ())
end

function crypto_hash_sha512(out, in, inlen)
    ccall((:crypto_hash_sha512, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}, Culonglong), out, in, inlen)
end

function crypto_hash_sha512_init(state)
    ccall((:crypto_hash_sha512_init, libsodium), Cint, (Ptr{crypto_hash_sha512_state},), state)
end

function crypto_hash_sha512_update(state, in, inlen)
    ccall((:crypto_hash_sha512_update, libsodium), Cint, (Ptr{crypto_hash_sha512_state}, Ptr{Cuchar}, Culonglong), state, in, inlen)
end

function crypto_hash_sha512_final(state, out)
    ccall((:crypto_hash_sha512_final, libsodium), Cint, (Ptr{crypto_hash_sha512_state}, Ptr{Cuchar}), state, out)
end
# Julia wrapper for header: crypto_sign_ed25519.h
# Automatically generated using Clang.jl

function crypto_sign_ed25519ph_statebytes()
    ccall((:crypto_sign_ed25519ph_statebytes, libsodium), Csize_t, ())
end

function crypto_sign_ed25519_bytes()
    ccall((:crypto_sign_ed25519_bytes, libsodium), Csize_t, ())
end

function crypto_sign_ed25519_seedbytes()
    ccall((:crypto_sign_ed25519_seedbytes, libsodium), Csize_t, ())
end

function crypto_sign_ed25519_publickeybytes()
    ccall((:crypto_sign_ed25519_publickeybytes, libsodium), Csize_t, ())
end

function crypto_sign_ed25519_secretkeybytes()
    ccall((:crypto_sign_ed25519_secretkeybytes, libsodium), Csize_t, ())
end

function crypto_sign_ed25519_messagebytes_max()
    ccall((:crypto_sign_ed25519_messagebytes_max, libsodium), Csize_t, ())
end

function crypto_sign_ed25519(sm, smlen_p, m, mlen, sk)
    ccall((:crypto_sign_ed25519, libsodium), Cint, (Ptr{Cuchar}, Ptr{Culonglong}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}), sm, smlen_p, m, mlen, sk)
end

function crypto_sign_ed25519_detached(sig, siglen_p, m, mlen, sk)
    ccall((:crypto_sign_ed25519_detached, libsodium), Cint, (Ptr{Cuchar}, Ptr{Culonglong}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}), sig, siglen_p, m, mlen, sk)
end

function crypto_sign_ed25519_verify_detached(sig, m, mlen, pk)
    ccall((:crypto_sign_ed25519_verify_detached, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}, Culonglong, Ptr{Cuchar}), sig, m, mlen, pk)
end

function crypto_sign_ed25519_keypair(pk, sk)
    ccall((:crypto_sign_ed25519_keypair, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}), pk, sk)
end

function crypto_sign_ed25519_seed_keypair(pk, sk, seed)
    ccall((:crypto_sign_ed25519_seed_keypair, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}, Ptr{Cuchar}), pk, sk, seed)
end

function crypto_sign_ed25519_pk_to_curve25519(curve25519_pk, ed25519_pk)
    ccall((:crypto_sign_ed25519_pk_to_curve25519, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}), curve25519_pk, ed25519_pk)
end

function crypto_sign_ed25519_sk_to_curve25519(curve25519_sk, ed25519_sk)
    ccall((:crypto_sign_ed25519_sk_to_curve25519, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}), curve25519_sk, ed25519_sk)
end

function crypto_sign_ed25519_sk_to_seed(seed, sk)
    ccall((:crypto_sign_ed25519_sk_to_seed, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}), seed, sk)
end

function crypto_sign_ed25519_sk_to_pk(pk, sk)
    ccall((:crypto_sign_ed25519_sk_to_pk, libsodium), Cint, (Ptr{Cuchar}, Ptr{Cuchar}), pk, sk)
end

function crypto_sign_ed25519ph_init(state)
    ccall((:crypto_sign_ed25519ph_init, libsodium), Cint, (Ptr{crypto_sign_ed25519ph_state},), state)
end

function crypto_sign_ed25519ph_update(state, m, mlen)
    ccall((:crypto_sign_ed25519ph_update, libsodium), Cint, (Ptr{crypto_sign_ed25519ph_state}, Ptr{Cuchar}, Culonglong), state, m, mlen)
end

function crypto_sign_ed25519ph_final_create(state, sig, siglen_p, sk)
    ccall((:crypto_sign_ed25519ph_final_create, libsodium), Cint, (Ptr{crypto_sign_ed25519ph_state}, Ptr{Cuchar}, Ptr{Culonglong}, Ptr{Cuchar}), state, sig, siglen_p, sk)
end

function crypto_sign_ed25519ph_final_verify(state, sig, pk)
    ccall((:crypto_sign_ed25519ph_final_verify, libsodium), Cint, (Ptr{crypto_sign_ed25519ph_state}, Ptr{Cuchar}, Ptr{Cuchar}), state, sig, pk)
end


text:
"""
$(TYPEDEF)

A floating joint type that uses a unit quaternion representation for orientation.

Floating joints are 6-degree-of-freedom joints that are in a sense degenerate,
as they impose no constraints on the relative motion between two bodies.

The full, 7-dimensional configuration vector of a `QuaternionFloating` joint
type consists of a unit quaternion representing the orientation that rotates
vectors from the frame 'directly after' the joint to the frame 'directly before'
it, and a 3D position vector representing the origin of the frame after the
joint in the frame before the joint.

The 6-dimensional velocity vector of a `QuaternionFloating` joint is the twist
of the frame after the joint with respect to the frame before it, expressed in
the frame after the joint.
"""
struct QuaternionFloating{T} <: JointType{T} end

Base.show(io::IO, jt::QuaternionFloating) = print(io, "Quaternion floating joint")
Random.rand(::Type{QuaternionFloating{T}}) where {T} = QuaternionFloating{T}()

num_positions(::Type{<:QuaternionFloating}) = 7
num_velocities(::Type{<:QuaternionFloating}) = 6
has_fixed_subspaces(jt::QuaternionFloating) = true
isfloating(::Type{<:QuaternionFloating}) = true

@propagate_inbounds function rotation(jt::QuaternionFloating, q::AbstractVector, normalize::Bool = true)
    quat = Quat(q[1], q[2], q[3], q[4], normalize)
    quat
end

@propagate_inbounds function set_rotation!(q::AbstractVector, jt::QuaternionFloating, rot::Rotation{3})
    T = eltype(rot)
    quat = convert(Quat{T}, rot)
    q[1] = quat.w
    q[2] = quat.x
    q[3] = quat.y
    q[4] = quat.z
    nothing
end

@propagate_inbounds function set_rotation!(q::AbstractVector, jt::QuaternionFloating, rot::AbstractVector)
    q[1] = rot[1]
    q[2] = rot[2]
    q[3] = rot[3]
    q[4] = rot[4]
    nothing
end

@propagate_inbounds translation(jt::QuaternionFloating, q::AbstractVector) = SVector(q[5], q[6], q[7])
@propagate_inbounds set_translation!(q::AbstractVector, jt::QuaternionFloating, trans::AbstractVector) = copyto!(q, 5, trans, 1, 3)

@propagate_inbounds angular_velocity(jt::QuaternionFloating, v::AbstractVector) = SVector(v[1], v[2], v[3])
@propagate_inbounds set_angular_velocity!(v::AbstractVector, jt::QuaternionFloating, ω::AbstractVector) = copyto!(v, 1, ω, 1, 3)

@propagate_inbounds linear_velocity(jt::QuaternionFloating, v::AbstractVector) = SVector(v[4], v[5], v[6])
@propagate_inbounds set_linear_velocity!(v::AbstractVector, jt::QuaternionFloating, ν::AbstractVector) = copyto!(v, 4, ν, 1, 3)

@propagate_inbounds function set_configuration!(q::AbstractVector, joint::Joint{<:Any, <:QuaternionFloating}, config::Transform3D)
    @boundscheck check_num_positions(joint, q)
    @framecheck config.from frame_after(joint)
    @framecheck config.to frame_before(joint)
    set_rotation!(q, joint_type(joint), rotation(config))
    set_translation!(q, joint_type(joint), translation(config))
    q
end

@propagate_inbounds function set_velocity!(v::AbstractVector, joint::Joint{<:Any, <:QuaternionFloating}, twist::Twist)
    @boundscheck check_num_velocities(joint, v)
    @framecheck twist.base frame_before(joint)
    @framecheck twist.body frame_after(joint)
    @framecheck twist.frame frame_after(joint)
    set_angular_velocity!(v, joint_type(joint), angular(twist))
    set_linear_velocity!(v, joint_type(joint), linear(twist))
    v
end

@propagate_inbounds function joint_transform(jt::QuaternionFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D, q::AbstractVector)
    Transform3D(frame_after, frame_before, rotation(jt, q, false), translation(jt, q))
end

@inline function motion_subspace(jt::QuaternionFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector)
    S = promote_eltype(jt, q)
    angular = hcat(one(SMatrix{3, 3, S}), zero(SMatrix{3, 3, S}))
    linear = hcat(zero(SMatrix{3, 3, S}), one(SMatrix{3, 3, S}))
    GeometricJacobian(frame_after, frame_before, frame_after, angular, linear)
end

@inline function constraint_wrench_subspace(jt::QuaternionFloating, joint_transform::Transform3D)
    S = promote_eltype(jt, joint_transform)
    WrenchMatrix(joint_transform.from, zero(SMatrix{3, 0, S}), zero(SMatrix{3, 0, S}))
end

@inline function bias_acceleration(jt::QuaternionFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    S = promote_eltype(q, v)
    zero(SpatialAcceleration{S}, frame_after, frame_before, frame_after)
end

@propagate_inbounds function configuration_derivative_to_velocity!(v::AbstractVector, jt::QuaternionFloating,
        q::AbstractVector, q̇::AbstractVector)
    quat = rotation(jt, q, false)
    quatdot = SVector(q̇[1], q̇[2], q̇[3], q̇[4])
    ω = angular_velocity_in_body(quat, quatdot)
    posdot = translation(jt, q̇)
    linear = inv(quat) * posdot
    set_angular_velocity!(v, jt, ω)
    set_linear_velocity!(v, jt, linear)
    nothing
end

@propagate_inbounds function configuration_derivative_to_velocity_adjoint!(fq, jt::QuaternionFloating, q::AbstractVector, fv)
    quatnorm = sqrt(q[1]^2 + q[2]^2 + q[3]^2 + q[4]^2) # TODO: make this nicer
    quat = Quat(q[1] / quatnorm, q[2] / quatnorm, q[3] / quatnorm, q[4] / quatnorm, false)
    rot = (velocity_jacobian(angular_velocity_in_body, quat)' * angular_velocity(jt, fv)) ./ quatnorm
    trans = quat * linear_velocity(jt, fv)
    set_rotation!(fq, jt, rot)
    set_translation!(fq, jt, trans)
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative!(q̇::AbstractVector, jt::QuaternionFloating,
        q::AbstractVector, v::AbstractVector)
    quat = rotation(jt, q, false)
    ω = angular_velocity(jt, v)
    linear = linear_velocity(jt, v)
    quatdot = quaternion_derivative(quat, ω)
    transdot = quat * linear
    set_rotation!(q̇, jt, quatdot)
    set_translation!(q̇, jt, transdot)
    nothing
end

@propagate_inbounds function velocity_to_configuration_derivative_jacobian(jt::QuaternionFloating, q::AbstractVector)
    quat = rotation(jt, q, false)
    vj = velocity_jacobian(quaternion_derivative, quat)
    R = RotMatrix(quat)
    # TODO: use hvcat once it's as fast
    @SMatrix(
        [vj[1] vj[5] vj[9]  0    0    0;
         vj[2] vj[6] vj[10] 0    0    0;
         vj[3] vj[7] vj[11] 0    0    0;
         vj[4] vj[8] vj[12] 0    0    0;
         0     0     0      R[1] R[4] R[7];
         0     0     0      R[2] R[5] R[8];
         0     0     0      R[3] R[6] R[9]])
end

@propagate_inbounds function configuration_derivative_to_velocity_jacobian(jt::QuaternionFloating, q::AbstractVector)
    quat = rotation(jt, q, false)
    vj = velocity_jacobian(angular_velocity_in_body, quat)
    R_inv = RotMatrix(inv(quat))
    # TODO: use hvcat once it's as fast
    @SMatrix(
        [vj[1] vj[4] vj[7] vj[10] 0        0        0;
         vj[2] vj[5] vj[8] vj[11] 0        0        0;
         vj[3] vj[6] vj[9] vj[12] 0        0        0;
         0     0     0     0      R_inv[1] R_inv[4] R_inv[7];
         0     0     0     0      R_inv[2] R_inv[5] R_inv[8];
         0     0     0     0      R_inv[3] R_inv[6] R_inv[9]])
end

@propagate_inbounds function zero_configuration!(q::AbstractVector, jt::QuaternionFloating)
    T = eltype(q)
    set_rotation!(q, jt, one(Quat{T}))
    set_translation!(q, jt, zero(SVector{3, T}))
    nothing
end

@propagate_inbounds function rand_configuration!(q::AbstractVector, jt::QuaternionFloating)
    T = eltype(q)
    set_rotation!(q, jt, rand(Quat{T}))
    set_translation!(q, jt, rand(SVector{3, T}) - 0.5)
    nothing
end

@propagate_inbounds function joint_twist(jt::QuaternionFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector)
    S = promote_eltype(jt, q, v)
    angular = convert(SVector{3, S}, angular_velocity(jt, v))
    linear = convert(SVector{3, S}, linear_velocity(jt, v))
    Twist(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function joint_spatial_acceleration(jt::QuaternionFloating, frame_after::CartesianFrame3D, frame_before::CartesianFrame3D,
        q::AbstractVector, v::AbstractVector, vd::AbstractVector)
    S = promote_eltype(jt, q, v, vd)
    angular = convert(SVector{3, S}, angular_velocity(jt, vd))
    linear = convert(SVector{3, S}, linear_velocity(jt, vd))
    SpatialAcceleration(frame_after, frame_before, frame_after, angular, linear)
end

@propagate_inbounds function joint_torque!(τ::AbstractVector, jt::QuaternionFloating, q::AbstractVector, joint_wrench::Wrench)
    set_angular_velocity!(τ, jt, angular(joint_wrench))
    set_linear_velocity!(τ, jt, linear(joint_wrench))
    nothing
end

# uses exponential coordinates centered around q0
@propagate_inbounds function local_coordinates!(ϕ::AbstractVector, ϕ̇::AbstractVector,
        jt::QuaternionFloating, q0::AbstractVector, q::AbstractVector, v::AbstractVector)
    # anonymous helper frames # FIXME
    frame_before = CartesianFrame3D()
    frame0 = CartesianFrame3D()
    frame_after = CartesianFrame3D()

    quat0 = rotation(jt, q0, false)
    quat = rotation(jt, q, false)
    p0 = translation(jt, q0)
    p = translation(jt, q)

    quat0inv = inv(quat0)
    δquat = quat0inv * quat
    δp = quat0inv * (p - p0)
    relative_transform = Transform3D(frame_after, frame0, δquat, δp)

    twist = joint_twist(jt, frame_after, frame0, q, v) # (q_0 is assumed not to change)
    ξ, ξ̇ = log_with_time_derivative(relative_transform, twist)

    copyto!(ϕ, 1, angular(ξ), 1, 3)
    copyto!(ϕ, 4, linear(ξ), 1, 3)

    copyto!(ϕ̇, 1, angular(ξ̇), 1, 3)
    copyto!(ϕ̇, 4, linear(ξ̇), 1, 3)

    nothing
end

@propagate_inbounds function global_coordinates!(q::AbstractVector, jt::QuaternionFloating, q0::AbstractVector, ϕ::AbstractVector)
    # anonymous helper frames #FIXME
    frame_before = CartesianFrame3D()
    frame0 = CartesianFrame3D()
    frame_after = CartesianFrame3D()

    t0 = joint_transform(jt, frame0, frame_before, q0)
    ξrot = SVector(ϕ[1], ϕ[2], ϕ[3])
    ξtrans = SVector(ϕ[4], ϕ[5], ϕ[6])
    ξ = Twist(frame_after, frame0, frame0, ξrot, ξtrans)
    relative_transform = exp(ξ)
    t = t0 * relative_transform
    set_rotation!(q, jt, rotation(t))
    set_translation!(q, jt, translation(t))
    nothing
end

@propagate_inbounds normalize_configuration!(q::AbstractVector, jt::QuaternionFloating) = set_rotation!(q, jt, rotation(jt, q, true))

@propagate_inbounds function is_configuration_normalized(jt::QuaternionFloating, q::AbstractVector, rtol, atol)
    isapprox(quatnorm(rotation(jt, q, false)), one(eltype(q)); rtol = rtol, atol = atol)
end

@propagate_inbounds function principal_value!(q::AbstractVector, jt::QuaternionFloating)
    quat = rotation(jt, q, false)
    set_rotation!(q, jt, principal_value(quat))
    nothing
end


text:
header("NODES")

@testset "place" begin
    node = xml"""
      <place id="place1">
        <name> <text>Some place</text> </name>
        <initialMarking> <text>100</text> </initialMarking>
      </place>
    """
    @test_call parse_node(node; reg = PNML.IDRegistry())
    n = @inferred Place parse_node(node; reg = PNML.IDRegistry())
    printnode(n)
    @test typeof(n) <: Place
    @test_call PNML.has_xml(n)
    @test !PNML.has_xml(n)
    @test @inferred(pid(n)) === :place1
    @test @inferred PNML.has_name(n)
    @test @inferred(PNML.name(n)) == "Some place"
    @test_call marking(n)
    @test marking(n) == 100
end

@testset "no text" begin
    node = xml"""
      <place id="place1">
        <name> <text>Some place</text> </name>
        <initialMarking>100</initialMarking>
      </place>
    """
    @test_call parse_node(node; reg = PNML.IDRegistry())
    n = @inferred Place parse_node(node; reg = PNML.IDRegistry())
    printnode(n)
    @test typeof(n) <: Place
    @test_call PNML.has_xml(n)
    @test !PNML.has_xml(n)
    @test @inferred(pid(n)) === :place1
    @test @inferred PNML.has_name(n)
    @test @inferred(PNML.name(n)) == "Some place"
    @test_call marking(n)
    @test marking(n) == 100
end

@testset "transition" begin
    node = xml"""
      <transition id="transition1">
        <name> <text>Some transition</text> </name>
        <condition> <structure>100</structure> </condition>
      </transition>
    """
    n = @inferred Transition parse_node(node; reg = PNML.IDRegistry())
    printnode(n)
    @test typeof(n) <: PNML.Transition
    @test !PNML.has_xml(n)
    @test pid(n) === :transition1
    @test PNML.has_name(n)
    @test PNML.name(n) == "Some transition"
    @test condition(n) == 100

    node = xml"""<transition id ="t1"> <condition><text>test</text></condition></transition>"""
    @test_throws ErrorException parse_node(node; reg = PNML.IDRegistry())
    node = xml"""<transition id ="t2"> <condition/> </transition>"""
    @test_throws ErrorException parse_node(node; reg = PNML.IDRegistry())
    node = xml"""<transition id ="t3"> <condition><structure/></condition> </transition>"""
    t = parse_node(node; reg = PNML.IDRegistry())
    @test t isa PNML.Transition
    @test_call condition(t)
    @test condition(t) === true
end

@testset "arc" begin
    node = xml"""
      <arc source="transition1" target="place1" id="arc1">
        <name> <text>Some arc</text> </name>
        <inscription> <text>6</text> </inscription>
      </arc>
    """
    n = @inferred Arc parse_node(node; reg = PNML.IDRegistry())
    printnode(n)
    @test typeof(n) <: PNML.Arc
    @test !PNML.has_xml(n)
    @test pid(n) === :arc1
    @test PNML.has_name(n)
    @test PNML.name(n) == "Some arc"
    @test_call inscription(n)
    @test PNML.inscription(n) == 6
end

@testset "ref Trans" begin
    node = xml"""
        <referenceTransition id="rt1" ref="t1"/>
    """
    n = parse_node(node; reg = PNML.IDRegistry())
    printnode(n)
    @test typeof(n) <: PNML.RefTransition
    @test !PNML.has_xml(n)
    @test pid(n) === :rt1
    @test n.ref === :t1
end

@testset "ref Place" begin
    n1 = (node = xml"""
    <referencePlace id="rp2" ref="rp1"/>
    """,
    id="rp2", ref="rp1" )
    n2 = (node = xml"""
    <referencePlace id="rp1" ref="Sync1">
        <graphics>
          <position x="734.5" y="41.5"/>
          <dimension x="40.0" y="40.0"/>
        </graphics>
    </referencePlace>
    """,
    id="rp1", ref="Sync1")
    @testset for s in [n1, n2]
        n = parse_node(s.node; reg = PNML.IDRegistry())
        printnode(n)
        @test typeof(n) <: PNML.RefPlace
        @test !PNML.has_xml(n)
        @test typeof(n.id) == Symbol
        @test typeof(n.ref) == Symbol
        @test n.id === Symbol(s.id)
        @test n.ref === Symbol(s.ref)
    end
end


text:
"""
$(SIGNATURES)

Helper function to process an individual block when it's a `HFun` such as
`{{ fill author }}`. See also [`convert_html`](@ref).
"""
function convert_html_fblock(β::HFun)::String
    fun = Symbol("hfun_" * lowercase(β.fname))
    ex  = isempty(β.params) ? :($fun()) : :($fun($β.params))
    # see if a hfun was defined in utils
    if isdefined(Main, :Utils) && isdefined(Main.Utils, fun)
        res = Core.eval(Main.Utils, ex)
        return string(res)
    end
    # see if a hfun was defined internally
    isdefined(Franklin, fun) && return eval(ex)
    # if zero parameters, see if can fill (case: {{vname}})
    if isempty(β.params) &&
        (!isnothing(locvar(β.fname)) || β.fname in UTILS_NAMES)
        return hfun_fill([β.fname])
    end
    # if we get here, then the function name is unknown, warn and ignore
    @warn "I found a function block '{{$(β.fname) ...}}' but I don't " *
          "recognise the function name. Ignoring."
    # returning empty
    return ""
end

"""
$(SIGNATURES)

H-Function of the form `{{ fill vname }}` or `{{ fill vname rpath}}` to plug in
the content of a franklin-var `vname` (assuming it can be represented as a
string).
"""
function hfun_fill(params::Vector{String})::String
    # check params
    if length(params) > 2 || isempty(params)
        throw(HTMLFunctionError("{{fill ...}} should have one or two " *
                                "($(length(params)) given). Verify."))
    end
    # form the fill
    repl  = ""
    vname = params[1]
    if length(params) == 1
        if vname in UTILS_NAMES
            repl = string(getfield(Main.Utils, Symbol(vname)))
        else
            tmp_repl = locvar(vname)
            if isnothing(tmp_repl)
                @warn "I found a '{{fill $vname}}' but I do not know the " *
                      "variable '$vname'. Ignoring."
            else
                repl = string(tmp_repl)
            end
        end
    else # two parameters, look in a path
        rpath = params[2]
        tmp_repl = pagevar(rpath, vname)
        if isnothing(tmp_repl)
            @warn "I found a '{{fill $vname $rpath}}' but I do not know the " *
                  "variable '$vname' or the path '$rpath'. Ignoring."
        else
            repl = string(tmp_repl)
        end
    end
    return repl
end

"""
$(SIGNATURES)

H-Function of the form `{{ insert fpath }}` to plug in the content of a file at
`fpath`. Note that the base path is assumed to be `PATHS[:src_html]`
(`< v"0.2"`) and `PATHS[:layout]` otherwise and so paths have to be expressed
relative to that.
"""
function hfun_insert(params::Vector{String})::String
    # check params
    if length(params) != 1
        throw(HTMLFunctionError("I found a {{insert ...}} with more than one parameter. Verify."))
    end
    # apply
    repl   = ""
    layout = path(ifelse(FD_ENV[:STRUCTURE] < v"0.2", :src_html, :layout))
    fpath  = joinpath(layout, split(params[1], "/")...)
    if isfile(fpath)
        repl = convert_html(read(fpath, String))
    else
        @warn "I found an {{insert ...}} block and tried to insert '$fpath' but I couldn't find the file. Ignoring."
    end
    return repl
end

"""
$(SIGNATURES)

H-Function of the form `{{href ... }}`.
"""
function hfun_href(params::Vector{String})::String
    # check params
    if length(params) != 2
        throw(HTMLFunctionError("I found an {{href ...}} block and expected 2 parameters" *
                                "but got $(length(params)). Verify."))
    end
    # apply
    repl = "??"
    dname, hkey = params[1], params[2]
    if params[1] == "EQR"
        haskey(PAGE_EQREFS, hkey) || return repl
        repl = html_ahref_key(hkey, PAGE_EQREFS[hkey])
    elseif params[1] == "BIBR"
        haskey(PAGE_BIBREFS, hkey) || return repl
        repl = html_ahref_key(hkey, PAGE_BIBREFS[hkey])
    else
        @warn "Unknown dictionary name $dname in {{href ...}}. Ignoring"
    end
    return repl
end

"""
$(SIGNATURES)

H-Function of the form `{{toc min max}}` (table of contents). Where `min` and
`max` control the minimum level and maximum level of  the table of content.
The split is as follows:

* key is the refstring
* f[1] is the title (header text)
* f[2] is irrelevant (occurence, used for numbering)
* f[3] is the level
"""
function hfun_toc(params::Vector{String})::String
    if length(params) != 2
        throw(HTMLFunctionError("I found a {{toc ...}} block and expected 2 " *
                              "parameters but got $(length(params)). Verify."))
    end
    isempty(PAGE_HEADERS) && return ""

    # try to parse min-max level
    min = 0
    max = 100
    try
        min = parse(Int, params[1])
        max = parse(Int, params[2])
    catch
        throw(HTMLFunctionError("I found a {{toc min max}} but couldn't " *
                                "parse min/max. Verify."))
    end

    inner   = ""
    headers = filter(p -> min ≤ p.second[3] ≤ max, PAGE_HEADERS)
    baselvl = minimum(h[3] for h in values(headers)) - 1
    curlvl  = baselvl
    for (rs, h) ∈ headers
        lvl = h[3]
        if lvl ≤ curlvl
            # Close previous list item
            inner *= "
"
            # Close additional sublists for each level eliminated
            for i = curlvl-1:-1:lvl
                inner *= "

"
            end
            # Reopen for this list item
            inner *= "
"
        elseif lvl > curlvl
            # Open additional sublists for each level added
            for i = curlvl+1:lvl
                inner *= "

"
            end
        end
        inner *= html_ahref_key(rs, h[1])
        curlvl = lvl
        # At this point, number of sublists (

) open equals curlvl
    end
    # Close remaining lists, as if going down to the base level
    for i = curlvl-1:-1:baselvl
        inner *= "

"
    end
    toc = "
" * inner * "
"
end

"""
$(SIGNATURES)

H-Function of the form `{{taglist}}`.
"""
function hfun_taglist()::String
    tag = locvar(:fd_tag)

    c = IOBuffer()
    write(c, "
")

    rpaths = globvar("fd_tag_pages")[tag]
    sorter(p) = begin
        pvd = pagevar(p, "date")
        if isnothing(pvd)
            return Date(Dates.unix2datetime(stat(p * ".md").ctime))
        end
        return pvd
    end
    sort!(rpaths, by=sorter, rev=true)

    (:hfun_list, "tag: $tag, loop over $rpaths") |> logger

    for rpath in rpaths
        title = pagevar(rpath, "title")
        if isnothing(title)
            title = "/$rpath/"
        end
        url = get_url(rpath)
        write(c, "
$title
")
    end
    write(c, "
")
    return String(take!(c))
end


text:
function build_jlData()
    fields = Tuple{Symbol,Type}[(:d, mjData), (:cptr, Ptr{mjData})]
    ctorargs = []

    getpropbody = :(getfield(value, name))
    setpropbody = :(setfield!(value, name, x))
    propnames = Symbol[]
    skippedfields = [:buffer, :stack]

    sinfo = STRUCT_INFO[:mjData]
    for (name, T) in zip(fieldnames(mjData), fieldtypes(mjData))
        if name in skippedfields
            continue
        end
        push!(propnames, name)
        if T <: Ptr
            if haskey(sinfo, name)
                finfo = sinfo[name]
                N = length(finfo.mjsize)
                T = Array{eltype(T),N}
                push!(fields, (name, T))
                push!(ctorargs, quote
                    ptr = getfield(d, $(QuoteNode(name)))
                    size = lookupsize(:mjData, m, $(QuoteNode(name)))
                    Base.unsafe_wrap(Array, ptr, size)
                end)
                if !isnothing(finfo.dsize)
                    getpropbody = Expr(
                        :elseif,
                        :(name === $(QuoteNode(name))),
                        quote
                            A = getfield(value, $(QuoteNode(name)))
                            d = getfield(value, :d)
                            sz = lookup_dsize(
                                Val{:mjData}(),
                                d,
                                Val{$(QuoteNode(name))}(),
                            )
                            UnsafeArray(pointer(A), sz)
                        end,
                        getpropbody,
                    )
                end
            else
                @warn "Found Ptr $name but no size, skipping"
            end
        elseif T <: StaticArray
            N = ndims(T)
            elT = eltype(T)
            sz = size(T)
            Tnew = Array{eltype(T),N}
            foffset = fieldoffset(mjData, Base.fieldindex(mjData, name))

            push!(fields, (name, Tnew))
            push!(ctorargs, quote
                ptr = Ptr{$elT}(pd + $foffset)
                Base.unsafe_wrap(Array, ptr, $sz)
            end)
        else
            # forward
            getpropbody = Expr(
                :elseif,
                :(name === $(QuoteNode(name))),
                :(getfield(getfield(value, :d), $(QuoteNode(name)))),
                getpropbody,
            )
            setpropbody = Expr(
                :elseif,
                :(name === $(QuoteNode(name))),
                :(setfield!(getfield(value, :d), $(QuoteNode(name)), convert($T, x))),
                setpropbody,
            )
        end
    end
    structdef = build_structdef(:jlData, fields, mutable = true, supertype = AbstractMJData)
    structexpr = combinestructdef(structdef)

    getpropbody = Expr(:if, getpropbody.args...)
    setpropbody = Expr(:if, setpropbody.args...)

    quote
        $structexpr
        function Base.getproperty(value::jlData, name::Symbol)
            $getpropbody
        end
        function Base.setproperty!(value::jlData, name::Symbol, x)
            $setpropbody
        end
        function Base.propertynames(::Union{jlData,Type{jlData}})
            $propnames
        end
        function jlData(
            pm::Ptr{mjModel},
            pd::Ptr{mjData};
            ownmodel::Bool = false,
            owndata::Bool = false,
        )
            m = unsafe_load(pm)
            d = unsafe_load(pd)
            jd = jlData(d, pd, $(ctorargs...))
            mj_resetData(pm, jd)
            mj_forward(pm, jd)
            # must free C-allocated Ptr/pd!
            owndata && Base.finalizer(jd -> MJCore.mj_deleteData(jd.cptr), jd)
            ownmodel && mj_deleteModel(pm)
            jd
        end
        function jlData(pm::Ptr{mjModel}; own::Bool = false)
            jlData(pm, mj_makeData(pm), ownmodel = own, owndata = true)
        end
        jlData(jm::jlModel) = jlData(jm.cptr, own = false)
        Base.cconvert(::Type{Ptr{mjData}}, jd::jlData) = Base.cconvert(Ptr{mjData}, jd.d)
        Base.deepcopy_internal(jd::jlData, ::IdDict) = Base.copy(jd)
        function Base.copy(jd::jlData)
            error("jlData is not deep copiable. See `mj_copyData` and `mj_makeData` instead.")
        end

    end
end

@eval $(build_jlData())


text:
"""
    transfer_from_left(X, [O,] A, [B=A])

```
   Λ
  ╱ ╲  ┌───────┐
┌1 X 2─1   A   2─(2)
│ ╲ ╱  └───3───┘
│  V       │
│         .2.
│        ( O )
│         `1'
│          │
│      ┌───3───┐
└──────1conj(B)2─┐
       └───────┘ │
(1)──────────────┘
```
"""
function transfer_from_left(X::AbstractMatrix, A::AbstractTensor3, B::AbstractTensor3 = A)
    @tensoropt (p = 1, (br, kr, bl, kl) = χ) transfer_from_left[br, kr] :=
        X[bl, kl] * A[kl, kr, p] * conj(B)[bl, br, p]
end
function transfer_from_left(
    X::UniformScaling{T},
    A::AbstractTensor3,
    B::AbstractTensor3 = A,
) where {T}
    @tensoropt (p = 1, (br, kr, l) = χ) transfer_from_left[br, kr] :=
        one(T) * A[l, kr, p] * conj(B)[l, br, p]
end
function transfer_from_left(
    X::AbstractMatrix,
    O::AbstractMatrix,
    A::AbstractTensor3,
    B::AbstractTensor3 = A,
)
    @tensoropt ((bp, kp) = 1, (br, kr, bl, kl) = χ) transfer_from_left[br, kr] :=
        X[bl, kl] * A[kl, kr, kp] * conj(B)[bl, br, bp] * O[bp, kp]
end
function transfer_from_left(
    X::UniformScaling{T},
    O::AbstractMatrix,
    A::AbstractTensor3,
    B::AbstractTensor3 = A,
) where {T}
    @tensoropt ((bp, kp) = 1, (br, kr, l) = χ) transfer_from_left[br, kr] :=
        one(T) * A[l, kr, kp] * conj(B)[l, br, bp] * O[bp, kp]
end

"""
    transfer_from_left(X, [O,] (AL, AR), [(BL, BR)=(AL, AR)])

```
   Λ
  ╱ ╲  ┌────────┐  ┌────────┐
┌1 X 2─1   AL   2──1   AR   2─(2)
│ ╲ ╱  └────3───┘  └────3───┘
│  V        │           │
│          .3───────────4.
│         (       O       )
│          `1───────────2'
│           │           │
│      ┌────3───┐  ┌────3───┐
└──────1conj(BL)2──1conj(BR)2─┐
       └────────┘  └────────┘ │
(1)───────────────────────────┘
```
"""
function transfer_from_left(
    X::AbstractMatrix,
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
)
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((lp, rp) = 1, (br, kr, bl, kl, b, k) = χ) transfer_from_left[br, kr] :=
        X[bl, kl] *
        AL[kl, k, lp] *
        AR[k, kr, rp] *
        conj(BL)[bl, b, lp] *
        conj(BR)[b, br, rp]
end
function transfer_from_left(
    X::UniformScaling{T},
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
) where {T}
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((lp, rp) = 1, (br, kr, l, b, k) = χ) transfer_from_left[br, kr] :=
        one(T) * AL[l, k, lp] * AR[k, kr, rp] * conj(BL)[l, b, lp] * conj(BR)[b, br, rp]
end
function transfer_from_left(
    X::AbstractMatrix,
    O::AbstractTensor4,
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
)
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((blp, brp, klp, krp) = 1, (br, kr, bl, kl, b, k) = χ) transfer_from_left[
        br,
        kr,
    ] :=
        X[bl, kl] *
        AL[kl, k, klp] *
        AR[k, kr, krp] *
        conj(BL)[bl, b, blp] *
        conj(BR)[b, br, brp] *
        O[blp, brp, klp, krp]
end
function transfer_from_left(
    X::UniformScaling{T},
    O::AbstractTensor4,
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
) where {T}
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((blp, brp, klp, krp) = 1, (br, kr, l, b, k) = χ) transfer_from_left[
        br,
        kr,
    ] :=
        one(T) *
        AL[l, k, klp] *
        AR[k, kr, krp] *
        conj(BL)[l, b, blp] *
        conj(BR)[b, br, brp] *
        O[blp, brp, klp, krp]
end

"""
    transfer_from_right(X, [O,] A, [B=A])

```
                Λ
    ┌───────┐  ╱ ╲
(1)─1   A   2─1 X 2─┐
    └───3───┘  ╲ ╱  │
        │       V   │
       .2.          │
      ( O )         │
       `1'          │
        │           │
    ┌───3───┐       │
┌───1conj(B)2───────┘
│   └───────┘
└───────────────────(2)
```
"""
function transfer_from_right(X::AbstractMatrix, A::AbstractTensor3, B::AbstractTensor3 = A)
    @tensoropt (p = 1, (br, kr, bl, kl) = χ) transfer_from_right[kl, bl] :=
        A[kl, kr, p] * conj(B)[bl, br, p] * X[kr, br]
end
function transfer_from_right(
    X::UniformScaling{T},
    A::AbstractTensor3,
    B::AbstractTensor3 = A,
) where {T}
    @tensoropt (p = 1, (r, bl, kl) = χ) transfer_from_right[kl, bl] :=
        one(T) * A[kl, r, p] * conj(B)[bl, r, p]
end
function transfer_from_right(
    X::AbstractMatrix,
    O::AbstractMatrix,
    A::AbstractTensor3,
    B::AbstractTensor3 = A,
)
    @tensoropt ((bp, kp) = 1, (br, kr, bl, kl) = χ) transfer_from_right[kl, bl] :=
        A[kl, kr, kp] * conj(B)[bl, br, bp] * X[kr, br] * O[bp, kp]
end
function transfer_from_right(
    X::UniformScaling{T},
    O::AbstractMatrix,
    A::AbstractTensor3,
    B::AbstractTensor3 = A,
) where {T}
    @tensoropt ((bp, kp) = 1, (r, bl, kl) = χ) transfer_from_right[kl, bl] :=
        one(T) * A[kl, r, kp] * conj(B)[bl, r, bp] * O[bp, kp]
end

"""
    transfer_from_right(X, [O,] (AL, AR), [(BL, BR)=(AL, AR)])

```
                             Λ
    ┌────────┐  ┌────────┐  ╱ ╲
(1)─1   AL   2──1   AR   2─1 X 2┐
    └────3───┘  └────3───┘  ╲ ╱ │
         │           │       V  │
        .3───────────4.         │
       (       O       )        │
        `1───────────2'         │
         │           │          │
    ┌────3───┐  ┌────3───┐      │
┌───1conj(BL)2──1conj(BR)2──────┘
│   └────────┘  └────────┘
└───────────────────────────────(2)
```
"""
function transfer_from_right(
    X::AbstractMatrix,
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
)
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((lp, rp) = 1, (br, kr, bl, kl, b, k) = χ) transfer_from_right[kl, bl] :=
        AL[kl, k, lp] *
        AR[k, kr, rp] *
        conj(BL)[bl, b, lp] *
        conj(BR)[b, br, rp] *
        X[kr, br]
end
function transfer_from_right(
    X::UniformScaling{T},
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
) where {T}
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((lp, rp) = 1, (r, bl, kl, b, k) = χ) transfer_from_right[kl, bl] :=
        one(T) * AL[kl, k, lp] * AR[k, r, rp] * conj(BL)[bl, b, lp] * conj(BR)[b, r, rp]
end
function transfer_from_right(
    X::AbstractMatrix,
    O::AbstractTensor4,
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
)
    (AL, AR), (BL, BR) = A, B
    @tensoropt ((blp, brp, klp, krp) = 1, (br, kr, bl, kl, b, k) = χ) transfer_from_right[
        kl,
        bl,
    ] :=
        AL[kl, k, klp] *
        AR[k, kr, krp] *
        conj(BL)[bl, b, blp] *
        conj(BR)[b, br, brp] *
        X[kr, br] *
        O[blp, brp, klp, krp]
end
function transfer_from_right(
    X::UniformScaling{T},
    O::AbstractTensor4,
    A::Tuple{<:AbstractTensor3,<:AbstractTensor3},
    B::Tuple{<:AbstractTensor3,<:AbstractTensor3} = A,
) where {T}

    (AL, AR), (BL, BR) = A, B
    @tensoropt ((blp, brp, klp, krp) = 1, (r, bl, kl, b, k) = χ) transfer_from_right[
        kl,
        bl,
    ] :=
        one(T) *
        AL[kl, k, klp] *
        AR[k, r, krp] *
        conj(BL)[bl, b, blp] *
        conj(BR)[b, r, brp] *
        O[blp, brp, klp, krp]
end


text:
using LAP_julia: average_dicts, print_dict

# prepare new dataframe
calculate_flow_df = DataFrame(
    index = Int[],
    mode = Symbol[],
    img_size = Int[],
    whs = Int[],
    timer = TimerOutput[],
    results = Dict{String,Float64}[]
    )

window_half_sizes = cat(collect(1:51), collect(61:10:101), collect(141:40:381), dims=1);
img_sizes = [50, 100, 200, 400, 800, 1600];
modes = [:known_value; :calculation];

##
# fill dataframe
df = calculate_flow_df
let index = 0
    for img_size in img_sizes
        img, imgw, flow = gen_init(:chess, chess_args=[25, img_size/25])

        # limit the max window half size to 1/4 the image
        whs_limit = img_size/4
        whs_modified = filter(x -> x <= whs_limit, window_half_sizes)

        for whs in whs_modified
            for mode in modes

                window = [whs * 2 + 1, whs * 2 + 1]
                timer = TimerOutput("reg alg: lap")
                # run once:
                flow_est, source_reg, timer, results = test_registration_alg(lap, img, imgw, flow, [whs, window],
                    Dict(:timer => timer, :new_feature => mode == :calculation ? false : true), timer=timer, display=true)

                index = index + 1
                println("at index: ", index, " img_size: ", img_size, " whs: ", whs)
                push!(df, Dict(:index => index,
                               :mode => mode,
                               :img_size => img_size,
                               :whs => whs,
                               :timer => timer,
                               :results => results))
            end
        end
    end
end

new_avg_dicts = []
old_avg_dicts = []

for img_size in img_sizes

    n = count(k -> k < (img_size/4), window_half_sizes)

    new_dicts = df[(df.mode .== :known_value) .& (df.img_size .== img_size), :][:, :results][1:n]
    old_dicts = df[(df.mode .== :calculation) .& (df.img_size .== img_size), :][:, :results][1:n]

    new_avg = average_dicts(new_dicts)
    old_avg = average_dicts(old_dicts)

    append!(new_avg_dicts, new_avg)
    append!(old_avg_dicts, old_avg)

    print_dict(new_avg, "new_avg, with img_size: " * string(img_size))
    print_dict(old_avg, "old_avg, with img_size: " * string(img_size))
end

img_size = 50
n = count(k -> k < (img_size/4), window_half_sizes)
new_dicts = df[(df.mode .== :known_value) .& (df.img_size .== img_size), :][:, :results][1:n]
old_dicts = df[(df.mode .== :calculation) .& (df.img_size .== img_size), :][:, :results][1:n]


text:
# https://gitlab.com/ase/ase/blob/master/doc/ase/calculators/lammpsrun.rst

using ASEPotential

a = atoms([atom("Na", [0,0,0]), atom("Cl", [0,0,2.3608])], [6.5, 6.5, 7.7])
parameters = LAMMPSCalculatorParameters()

@time println(get_potential_energy(a, parameters))
@time println(get_forces(a, parameters))


text:
using Overseer: Indices, @indices_in, @reverse_indices_in

@testset "Construction, collect" begin
    data_in = (1,5,100)
    s = Indices(data_in)
    data_out = collect(s)
    @test all(map(d->in(d,data_out), data_in))
    @test length(data_out) == length(data_in)
end

@testset "eltype, empty" begin
    @test eltype(Indices()) == Int
    @test eltype(typeof(Indices())) == Int
    @test isequal(empty(Indices([1,2,3])), Indices())
end

@testset "Core Functionality" begin
    s = Indices([1,2,10,20,200,300,1000,10000,10002])
    @test last(s) == 10002
    @test first(s) == 1
    @test length(s) == 9
    @test pop!(s) == 10002
    @test length(s) == 8
    @test popfirst!(s) == 1
    @test length(s) == 7
    @test !in(1,s)
    @test !in(10002,s)
    @test in(10000,s)
    @test_throws ArgumentError first(Indices())
    @test_throws ArgumentError last(Indices())
    t = copy(s)
    s = Indices()
    push!(s, 1, 2, 100)
    @test 1 in s
    @test !(3 in s)
    @test 2 in s
    @test 100 in s
    @test !(101 in s)
    @test !(1000 in s)
    @test first(s) == 1
    @test last(s) == 100
    @test s == Indices([1, 2, 100])
    push!(s, 1000)
    @test [i for i in s] == [1, 2, 100, 1000]
    @test pop!(s) == 1000
    @test s == Indices([1, 2, 100])
    push!(s, 5000)
    push!(s, 2000)
    pop!(s, 5000)
    @test s.reverse[end] === Overseer.NULL_INT_PAGE
    b = 1:1000
    s = Indices(b)
    @test collect(s) == collect(b)
    @test length(s) == length(b)
    @test pop!(s, 100) == 100
    @test_throws BoundsError pop!(s, 100)
    @test pop!(s, 100, 1) == 1
    @test pop!(s, 99, 1) == 99
    @test !in(500000, s)
    @test !in(99, s)
end

@testset "setdiff / symdiff" begin
    @test setdiff(Indices([1, 2, 3, 4]), Indices([2, 4, 5, 6])) == Indices([1, 3])
end

@testset "setdiff!" begin
    s2 = Indices([1, 2, 3, 4])
    setdiff!(s2, Indices([2, 4, 5, 6]))

    @test s2 == Indices([1, 3])
end

@testset "issue #7851" begin
    @test_throws DomainError Indices(-1)
    @test !(-1 in Indices(1:10))
end
@testset "Copy, copy!, empty" begin
    s1 = Indices([1,2,3])
    s2 = empty(s1)
    push!(s2, 10000)
    @test !in(10000, s1)
    copy!(s2, s1)
    @test !in(10000, s2)
    push!(s2, 10000)
    @test !in(10000, s1)
    s3 = copy(s2)
    push!(s3, 1000)
    @test !in(1000, s2)
    pop!(s3, 1000)
    pop!(s2, 10000)
    @test in(10000, s3)
    pop!(s3, 10000)
    @test s3 == s2 == s1
    @test collect(s3) == collect(s2) == [1,2,3]

end

@testset "Push, union" begin
    # Push, union
    s1 = Indices()
    @test_throws DomainError push!(s1, -1)
    push!(s1, 1, 10, 100, 1000)
    @test collect(s1) == [1, 10, 100, 1000]
    push!(s1, 606)
    @test collect(s1) == [1, 10, 100, 1000, 606]
    s2 = Indices()
    @test s2 === union!(s2, s1)
    s3 = Indices([1, 10, 100])
    union!(s3, [1, 606, 1000])
    s4 = union(Indices([1, 100, 1000]), Indices([10, 100, 606]))
    @test s1 == s2 == s3 == s4
end

@testset "pop!, delete!" begin
    s = Indices(1:2:10)
    @test pop!(s, 1) == 1
    @test !(1 in s)
    @test_throws BoundsError pop!(s, 1)
    @test_throws ArgumentError pop!(s, -1)
    @test_throws ArgumentError pop!(s, -1, 1)
    @test pop!(s, 1, 0) == 0
    is = copy(s.packed)
    for i in is; pop!(s, i); end
    @test isempty(s)
    push!(s, 1:2:10...)
    @test pop!(s) == 9
    @test pop!(s) == 7
    @test popfirst!(s) == 1
    @test collect(s) == [5, 3]
    empty!(s)
    @test isempty(s)
end

@testset "Intersect" begin
    @test isempty(intersect(Indices()))
    @test isempty(intersect(Indices(1:10), Indices()))
    @test isempty(intersect(Indices(), Indices(1:10)))

    @test intersect(Indices([1,2,3])) == Indices([1,2,3])

    @test intersect(Indices(1:7), Indices(3:10)) ==
          intersect(Indices(3:10), Indices(1:7)) == Indices(3:7)

    @test intersect!(Indices(1:10), Indices(1:4), 1:5, [1,2,10]) == Indices(1:2)
end

@testset "Setdiff" begin
    s1 = Indices(1:100)
    setdiff!(s1, Indices(1:2:100))
    s2 = setdiff(Indices(1:100), Indices(1:2:100))
    @test s1 == s2 == Indices(2:2:100)

    s1 = Indices(1:10)
    s2 = Indices([1:2; 6:100])
    @test setdiff(s1, s2) == setdiff(s1, [1:2; 6:100]) == Indices(3:5)
end

@testset "Subsets, equality" begin
    @test Indices(2:2:10) < Indices(1:10)
    @test !(Indices(2:2:10) < Indices(2:2:10))
    @test Indices(2:2:10) <= Indices(2:10)
    @test Indices(2:2:10) <= Indices(2:2:10)
end

@testset "indices_in" begin
    a = Indices(2:2:10)
    b = Indices(10:3:20)
    c = [Indices(3:10)]

    d = Indices(1)
    t = 0
    it = @indices_in(((a && c[1]) || d) && !b)
    for e in it
        t+=e
    end
    @test t == 4+6+8+1
    @test length(it) == 4
    t = 0
    it = @reverse_indices_in(((a && c[1]) || d) && !b)
    for e in it
        t+=e
    end
    @test t == 4+6+8+1
    @test length(it) == 4
    @test collect(it) == [1, 8, 6, 4]
end


text:
function sampling(;input::Input, interf::Interferometer, distinguishability::Real, reflectivity::Real, exact=nothing)

    if exact == nothing || exact == true

        if distinguishability == 1 && reflectivity == 1
            return cliffords_sampler(input=input, interf=interf)
        elseif distinguishability == 0 && reflectivity == 1
            return classical_sampler(input=input, interf=interf)
        else
            throw(ArgumentError("need approximation"))
        end

    else

        if distinguishability == 1 && reflectivity == 1
                  n = input.r.n
                  m = input.r.m
                  U = copy(interf.U)

                  # generate a collisionless state as a starting point
                  starting_state = iterate_until_collisionless(() -> random_occupancy(n,m))

                  known_pdf(state) = process_probability_distinguishable(U, input_state, state)
                  target_pdf(state) = process_probability(U, input_state, state)
                  known_sampler = () -> iterate_until_collisionless(() -> classical_sampler(U = U, m = m, n = n)) # gives a classical sampler

                  return metropolis_sampler(;target_pdf = target_pdf, known_pdf = known_pdf , known_sampler = known_sampler , starting_state = starting_state, n_iter = 100)
        else
            return noisy_sampling(input=input, distinguishability=distinguishability, reflectivity=reflectivity, interf=interf)
        end

    end

end


text:
struct UnitCommitment<:AbstractOperationsProblem end

function UnitCommitment(sys::PSY.System, transmission::Type{S}; optimizer::Union{Nothing, JuMP.OptimizerFactory}=nothing, kwargs...) where {S<:PM.AbstractPowerModel}

    devices = Dict{Symbol, DeviceModel}(:ThermalGenerators => DeviceModel(PSY.ThermalGen, ThermalStandardUnitCommitment),
                                            :RenewableGenerators => DeviceModel(PSY.RenewableGen, RenewableFullDispatch),
                                            :Loads => DeviceModel(PSY.PowerLoad, StaticPowerLoad))

branches = Dict{Symbol, DeviceModel}(:Lines => DeviceModel(PSY.Branch, SeriesLine))
services = Dict{Symbol, ServiceModel}(:Reserves => ServiceModel(PSY.Reserve, AbstractReservesFormulation))
    return OperationsProblem(UnitCommitment,
                                   transmission,
                                    devices,
                                    branches,
                                    services,
                                    system;
                                    optimizer = optimizer, kwargs...)
end


text:
#=
This file is auto-generated. Do not edit.
=#
"""
    mutable struct StaticReserve{T <: ReserveDirection} <: Reserve{T}
        name::String
        available::Bool
        timeframe::Float64
        requirement::Float64
        ext::Dict{String, Any}
        internal::InfrastructureSystemsInternal
    end

Data Structure for a proportional reserve product for system simulations.

# Arguments
- `name::String`
- `available::Bool`
- `timeframe::Float64`: the relative saturation timeframe, validation range: (0, nothing), action if invalid: error
- `requirement::Float64`: the static value of required reserves, validation range: (0, nothing), action if invalid: error
- `ext::Dict{String, Any}`
- `internal::InfrastructureSystemsInternal`: power system internal reference, do not modify
"""
mutable struct StaticReserve{T <: ReserveDirection} <: Reserve{T}
    name::String
    available::Bool
    "the relative saturation timeframe"
    timeframe::Float64
    "the static value of required reserves"
    requirement::Float64
    ext::Dict{String, Any}
    "power system internal reference, do not modify"
    internal::InfrastructureSystemsInternal
end

function StaticReserve{T}(name, available, timeframe, requirement, ext=Dict{String, Any}(), ) where T <: ReserveDirection
    StaticReserve{T}(name, available, timeframe, requirement, ext, InfrastructureSystemsInternal(), )
end

function StaticReserve{T}(; name, available, timeframe, requirement, ext=Dict{String, Any}(), ) where T <: ReserveDirection
    StaticReserve{T}(name, available, timeframe, requirement, ext, )
end

# Constructor for demo purposes; non-functional.
function StaticReserve{T}(::Nothing) where T <: ReserveDirection
    StaticReserve{T}(;
        name="init",
        available=false,
        timeframe=0.0,
        requirement=0.0,
        ext=Dict{String, Any}(),
    )
end

InfrastructureSystems.get_name(value::StaticReserve) = value.name
"""Get StaticReserve available."""
get_available(value::StaticReserve) = value.available
"""Get StaticReserve timeframe."""
get_timeframe(value::StaticReserve) = value.timeframe
"""Get StaticReserve requirement."""
get_requirement(value::StaticReserve) = value.requirement
"""Get StaticReserve ext."""
get_ext(value::StaticReserve) = value.ext
"""Get StaticReserve internal."""
get_internal(value::StaticReserve) = value.internal

InfrastructureSystems.set_name!(value::StaticReserve, val::String) = value.name = val
"""Set StaticReserve available."""
set_available!(value::StaticReserve, val::Bool) = value.available = val
"""Set StaticReserve timeframe."""
set_timeframe!(value::StaticReserve, val::Float64) = value.timeframe = val
"""Set StaticReserve requirement."""
set_requirement!(value::StaticReserve, val::Float64) = value.requirement = val
"""Set StaticReserve ext."""
set_ext!(value::StaticReserve, val::Dict{String, Any}) = value.ext = val
"""Set StaticReserve internal."""
set_internal!(value::StaticReserve, val::InfrastructureSystemsInternal) = value.internal = val


text:
"""
	$(SIGNATURES)

Study times permitted in a given college. In model units.

Study times can be per course or total.

Grid points INclude fixed time per course.

The grid (and the outside) never knows about the min time per course. It just chooses how much to study per course. Everything else is up to the h production function.
"""
mutable struct StudyTimeGrid
    studyTimeV :: Vector{Double}
    # fixedTimePerCourse :: Double
    # minTimePerCourse :: Double
    sTimeGridPerCourse :: Bool
end

function make_test_study_time_grid(perCourse :: Bool)
    return StudyTimeGrid([0.2, 0.4, 0.6], perCourse);
end

function validate_stg(stg :: StudyTimeGrid)
    isValid = all(stg.studyTimeV .> 0.0)  &&  all(diff(stg.studyTimeV) .> 0.0);
    return isValid
end

"""
	$(SIGNATURES)

Number of allowed study time grid points.
"""
n_study_times(stg :: StudyTimeGrid) = length(stg.studyTimeV);

"""
	$(SIGNATURES)

Is the study time grid per course (including any fixed time costs per course) or total?
"""
stime_grid_per_course(stg :: StudyTimeGrid) = stg.sTimeGridPerCourse;

"""
	$(SIGNATURES)

Study times (total) available, given no of courses tried.
"""
study_time_grid(stg :: StudyTimeGrid, nTried :: Integer) =
    study_time_from_grid(stg, stg.studyTimeV, nTried);

study_time_grid(stg :: StudyTimeGrid, nTried :: Integer, idxV) =
    study_time_from_grid(stg, stg.studyTimeV[idxV], nTried);

"""
	$(SIGNATURES)

Max total study time for this no of courses tried.
"""
max_study_time(stg :: StudyTimeGrid, nTried :: Integer) =
    last(study_time_grid(stg, nTried));

# Bound from above. ++++++
# Study time total, from study time per course (INCL fixed time).
function study_time_total(stg :: StudyTimeGrid, sTimePerCourseV, nTriedV)
    sTimeV = sTimePerCourseV .* nTriedV;
    return sTimeV
end

# Inputs are grid points (which can be total or per course)
function study_time_from_grid(stg :: StudyTimeGrid, sTimeGridV, nTriedV)
    if stime_grid_per_course(stg)
        return study_time_total(stg, sTimeGridV, nTriedV);
    else
        return sTimeGridV
    end
end

# fixed_time_per_course(stg :: StudyTimeGrid) = stg.fixedTimePerCourse;

## Study time per course (INCL fixed time), from total study time.
# function study_time_per_course(stg :: StudyTimeGrid, sTimeTotalV, nTriedV)
#     return max.(min_time_per_course(stg),  study_time_per_course_no_min(stg, sTimeTotalV, nTriedV))
# end

# study_time_per_course_no_min(stg :: StudyTimeGrid, sTimeTotalV, nTriedV) =
#     sTimeTotalV ./ nTriedV .- fixed_time_per_course(stg);

# stime_per_course_grid(stg :: StudyTimeGrid, nTried :: Integer) =
#     stime_per_course_from_grid(stg, stg.studyTimeV, nTried);

# # Inputs are study time grid points.
# function stime_per_course_from_grid(stg :: StudyTimeGrid,
#     sTimeGridV, nTried :: Integer)
#     if stime_grid_per_course(stg)
#         sTimeV = max.(min_time_per_course(stg),
#             sTimeGridV .- fixed_time_per_course(stg));
#     else
#         sTimeV = study_time_per_course(stg, sTimeGridV, nTried);
#     end
#     return sTimeV
# end

# # Study time total, from study time per course (INCL fixed time).
# function study_time_total(stg :: StudyTimeGrid, sTimePerCourseV, nTriedV)
#     sTimeV = sTimePerCourseV .* nTriedV;
#     return sTimeV
# end

# # Inputs are grid points
# function stime_total_from_grid(stg :: StudyTimeGrid, sTimeGridV, nTriedV)
#     if stime_grid_per_course(stg)
#         return study_time_total(stg, sTimeGridV, nTriedV);
#     else
#         return sTimeGridV
#     end
# end

# -------------


text:
matrixdata =
"""
"random symmetric matrix"
function randsym{T}(::Type{T}, n)
 A = zeros(T, n, n)
  for j = 1:n
      for i = j:n
          A[i,j] = randn()
      end
  end
  A = A + tril(A, -1)'
  return A
end
randsym(n) = randsym(Float64, n)
include_generator(FunctionName, "randsym", randsym)
include_generator(Group, "random", randsym)
include_generator(Group, "symmetric", randsym)
"""
user_dir = joinpath(dirname(@__FILE__), "..", "myMatrixDepot")

open(string(user_dir, "/generator.jl"), "w") do f
    write(f, matrixdata)
end
workspace()
using MatrixDepot
using Base.Test

n = rand(1:8)
A = matrixdepot("randsym", n)
matrixdepot("randsym")
@test "randsym" in matrixdepot("random")
@test "randsym" in matrixdepot("symmetric")

user_dir = joinpath(dirname(@__FILE__), "..", "myMatrixDepot")
rm(user_dir, recursive = true)
workspace()
using MatrixDepot
using Base.Test