text:
"restriction enzyme used for digesting genome [name, cut-sequence, cut-site, overhang]"
mutable struct restrictionEnzyme
    name::String # name of the restriction enzyme
    cutSeq::Array{String} # recognition sequence
    cutPlace::Array{Int64} # where it gets cut
    overhang::Array{String} # overhang of the restriction enzyme if exists
    restrictionEnzyme(name, cutSeq, cutPlace, overhang) =
        new(name, cutSeq, cutPlace, overhang)
end

## list of restraiction enzymes
## single-digest
ApeKI = restrictionEnzyme("ApeKI", ["GCAGC", "GCTGC"], [2, 2], ["CTG", "CAG"]); # Elshire et al. 2011
SbfI = restrictionEnzyme("SbfI", ["CCTGCAGG"], [3], ["TGCA"]);# Elshire et al. 2011
PstI = restrictionEnzyme("PstI", ["CTGCAG"], [2], ["TGCA"]); # Poland et al. 2012
MspI = restrictionEnzyme("MspI", ["CCGG"], [2], ["CGG"]); # Poland et al. 2012

## double-digest
PstI_MspI = restrictionEnzyme("PstI-MspI", ["CTGCAG", "CCGG"], [2, 2], ["TGCA", "CGG"]); # Poland et al. 2012
SbfI_MspI = restrictionEnzyme("SbfI-MspI", ["CCTGCAGG", "CCGG"], [2, 2], ["TGCA", "CGG"]); # Poland et al. 2012

## nucleotides and its complements
"nucleotides"
nucl = Dict(
    'A' => 'A',
    'C' => 'C',
    'G' => 'G',
    'T' => 'T',
    'N' => 'N',
    'R' => rand(['A', 'G']),
    'Y' => rand(['C', 'T']),
    'S' => rand(['G', 'C']),
    'W' => rand(['A', 'T']),
    'K' => rand(['G', 'T']),
    'M' => rand(['A', 'C']),
    'B' => rand(['C', 'G', 'T']),
    'D' => rand(['A', 'G', 'T']),
    'H' => rand(['A', 'C', 'T']),
    'V' => rand(['A', 'C', 'G']),
); # randomly assign IUPAC codes for uncertain ones in the assembly

"complements"
comp = Dict('A' => 'T', 'C' => 'G', 'G' => 'C', 'T' => 'A', 'N' => 'N'); # complements of DNA nucleotides

"split sequence by any specific sequence"
function splitSequence(seq, cutSeq, startSeq, endSeq)
    cuts = split(uppercase(seq), cutSeq) # split DNA sequence
    if length(cuts) > 1
        cuts[2:length(cuts)] = startSeq .* cuts[2:length(cuts)] # re-attach the strat seq to the fragment after cut-place
        cuts[1:(length(cuts)-1)] = cuts[1:(length(cuts)-1)] .* endSeq # re-attach the end seq to the fragment before cut-place
    else
        cuts
    end
    cuts
end

"digest sequence using specified restriction enzyme"
function digest(sequence, re::restrictionEnzyme)
    cutSeq = re.cutSeq
    cutPlace = re.cutPlace
    # overhang = re.overhang
    startSeq = [SubString(cutSeq[i], cutPlace[i]) for i = 1:length(cutSeq)]
    endSeq = [SubString(cutSeq[i], 1, cutPlace[i] - 1) for i = 1:length(cutSeq)] # [SubString(cutSeq[i], 1, cutPlace[i] - 1) * overhang[i] for i = 1:length(cutSeq)]
    for k = 1:length(cutSeq)
        cuts = sequence
        sequence = []
        for s in cuts
            sequence = [sequence; splitSequence(s, cutSeq[k], startSeq[k], endSeq[k])]
        end
        sequence
    end
    # if length(sequence) > 2
    #    sequence[2:length(sequence)]
    # else
    # end
    sequence
end

"double-digestion (or enzyme has more than one recongnition site)"
function digestSecond(sequence, re::restrictionEnzyme)
    cutSeq = re.cutSeq
    cutPlace = re.cutPlace
    # overhang = re.overhang
    startSeq = [SubString(cutSeq[i], cutPlace[i]) for i = 1:length(cutSeq)]
    endSeq = [SubString(cutSeq[i], 1, cutPlace[i] - 1) for i = 1:length(cutSeq)] #[SubString(cutSeq[i], 1, cutPlace[i] - 1) .* overhang[i] for i = 1:length(cutSeq)]
    goodSeq = []
    for k = 1:length(sequence)
        s = sequence[k]
        sp = splitSequence(uppercase(s), cutSeq[1], startSeq[1], endSeq[1])
        if length(sp) > 1
            if k == 1
                s2 = startSeq .* sp[length(sp)]
                goodSeq = [goodSeq; s2]
            elseif k == length(sequence)
                s1 = sp[1] .* endSeq
                goodSeq = [goodSeq; s1]
            else
                s1 = sp[1] .* endSeq
                s2 = startSeq .* sp[length(sp)]
                goodSeq = [goodSeq; s1; s2]
            end
        else
            goodSeq = [goodSeq; sp]
        end
    end
    goodSeq
end

"""
	digestGenome(genofile, re, useChr, useChrLen, lower ,upper, plotOutput, writeOutput)

Perform virtual digestion using restriction enzyme and generate GBS fragments.

This function uses specified restriction enzyme(s) to digest genome and therefore generate GBS fragments. Fragment size-selection step is also included.

# Arguments
- `genofile`: file containing the reference genome
- `re`: restriction enzyme(s) to be used
- `useChr`: either the number of chromosome or a set of chromosome(s) to be simulated
- `useChrLen`: length of chromosome in cM to be used in simulation, otherwise using entire chromosome
- `lower`: lower threshold of fragment size-selection
- `upper`: upper threshold of fragment size-selection
- `winSize`: size of window used for calculating average genomic coverage
- `plotOutput`: set to true if graphical outputs are required
- `writeOutput`: set to true if text outputs are required
...

# Examples
```julia
julia> digestGenome("ref.fa.gz", [SimGBS.ApeKI], [1], Array{Float64}(undef,0), 65 ,195, 1000000, false, true)
```
"""
function digestGenome(
    genofile::String,
    re::Array{restrictionEnzyme,1},
    useChr::Array{Int64},
    useChrLen::Array{Float64},
    lowerThresh::Int64,
    upperThresh::Int64,
    winSize::Int64,
    plotOutput::Bool,
    writeOutput::Bool,
)
    ## module 1: read genome and cut it into fragments
    genome = readdlm(GZip.open(genofile), '\t') # load genome
    index = LinearIndices(genome) # index genome
    header = index[findall(x -> x[1] == '>', genome)] # identify headers within each FASTA file
    global chrLen = Array{Float64}(undef, 0) # empty array to store estimated chromosome length
    seqFrag = Array{String}(undef, 0) # empty array to store the DNA sequence
    chrFrag = Array{Int64}(undef, 0) # empty array to store chromosome
    idFrag = Array{Int64}(undef, 0) # empty array to store the ID
    posFrag = Array{Int64}(undef, 0) # empty array to store the position
    lenFrag = Array{Int64}(undef, 0) # empty array to store the position
    if length(useChr) <= length(header)
        global numChr = length(useChr) # useChr is considered to be a set of chromosomes if a vector is supplied
        starts = (header.+1)[useChr] # select chromosomes (specified by [useChr])
        ends = [header[2:end] .- 1; size(genome)[1]][useChr] # subset genome
        println("INFO: Using chromsome $(useChr) in the simulation!")
    else
        error("Wrong chromosome number/set specified, please check your parameters!")
    end
    for c = 1:numChr
        seq = [join(genome[starts[c]:ends[c]])] # read genomic sequence
        ## virtual digestion
        if size(re, 1) == 1
            frag = digest(seq, re[1]) # single-digestion
            reName = re[1].name
        else
            raw = digest(seq, re[1])
            frag = digestSecond(raw, re[2]) # double-digestion if more than one enzymes is supplied
            reName = [re[1].name re[2].name]
        end
        numFrag = length(frag) # number of fragments
        println(
            "CHROMOSOME $(useChr[c]): $numFrag GBS fragments are generated through virtual digestion using $(reName).",
        )
        len = length.(frag) # length of fragments
        pos = [0; 1 .+ cumsum(len)[1:(end-1)]] # starting position of each fragment
        chrLenEst = pos[end] / 1e6 # estimate chromosome length
        println(
            "CHROMOSOME $(useChr[c]): Estimated chromosome length equals to $chrLenEst Mb.",
        )
        ## (optional) subsetting chromosome
        if useChrLen == [] || size(useChrLen, 1) < numChr
            keep = [1:length(frag)...] # keep all fragments if subsetting is not required
            chrLen = [chrLen; chrLenEst] # store chromosome lengths
            println(
                "CHROMOSOME $(useChr[c]): Subsetting chromosome is not required. Keeping all $numFrag available GBS fragments.",
            )
        elseif size(useChrLen, 1) == numChr || useChrLen[c] <= estLen
            keep = findall(x -> x < useChrLen[c] * 1e6, pos) # keep only fragments found within specified length (varied at each chromosome)
            chrLen = useChrLen # store user-defined chromosome lengths
            println(
                "CHROMOSOME $(useChr[c]) (up to $(useChrLen[c]) cM): $(length(keep)) out of $numFrag available GBS fragments are kept.",
            )
        elseif size(useChrLen, 1) == 1 # keep only fragments found within specified length (uniform across chromosomes)
            keep = findall(x -> x < useChrLen[1] * 1e6, pos)
            chrLen = [chrLen; useChrLen[1]]
            println(
                "CHROMOSOME $(useChr[c]) (up to $(useChrLen[c]) cM): $(length(keep)) out of $numFrag available GBS fragments are kept.",
            )
        else
            error(
                "CHROMOSOME $(useChr[c]): Wrong chromosome length(s) specified, please check your parameters!",
            )
        end
        idFrag = [idFrag; [1:length(keep)...]] # identifier of each fragment
        chrFrag = [chrFrag; fill(useChr[c], length(keep))] # origin of chromosome per fragment
        posFrag = [posFrag; pos[keep]] # starting position per fragment
        lenFrag = [lenFrag; len[keep]] # length of each fragment
        seqFrag = [seqFrag; frag[keep]] # genomic sequence of each fragment
    end
    ## module 2: fragment size-selection
    fragRaw = [idFrag chrFrag posFrag lenFrag seqFrag] # dataset stores all fragments
    numFragTotal = length(seqFrag) # total number of GBS fragments
    selected = findall(x -> lowerThresh <= x <= upperThresh, lenFrag) # fragment size-selection
    numFragSelected = length(selected) # number of fragments selected in size-selection
    println(
        "INFO: $numFragSelected out of $numFragTotal GBS fragments are selected after size-selection with lower and upper thresholds equals to [$lowerThresh, $upperThresh].",
    )
    fragGBS = fragRaw[selected, :] # dataset conatins GBS fragments (after size-selection)
    global GBSFrag = (
        frag = fragGBS[:, 5],
        len = fragGBS[:, 4],
        pos = fragGBS[:, 3],
        chr = fragGBS[:, 2],
        id = fragGBS[:, 1],
    )
    ## module 3: estimate coverage of selected GBS fragments (per Mb)
    GBSCoverChr = Array{Int64}(undef, 0)
    GBSFragWinStatrs = Array{Int64}(undef, 0)
    GBSFragWinEnds = Array{Int64}(undef, 0)
    GBSFragCover = Array{Float64}(undef, 0)
    for c = 1:numChr
        numWin = Int(round(chrLen[c]) .* 1e6 / winSize)
        winStarts = Int.([0:(numWin-1)...] .* winSize)
        winEnds = Int.([1:numWin...] .* winSize)
        GBSFragChr = findall(x -> x == useChr[c], GBSFrag.chr)
        GBSFragPos = GBSFrag.pos[GBSFragChr]
        GBSFragLen = GBSFrag.len[GBSFragChr]
        GBSFragCoverWin = Array{Float64}(undef, numWin)
        for w = 1:numWin
            GBSFragWin = findall(x -> x <= winEnds[w] && x >= winStarts[w], GBSFragPos)
            if size(GBSFragWin, 1) != 0
                GBSFragCoverWin[w] = sum(GBSFragLen[GBSFragWin]) / winSize
            else
                GBSFragCoverWin[w] = 0
            end
        end
        GBSFragCover = [GBSFragCover; GBSFragCoverWin]
        GBSFragWinEnds = [GBSFragWinEnds; winEnds]
        GBSFragWinStatrs = [GBSFragWinStatrs; winStarts]
        GBSCoverChr = [GBSCoverChr; fill(c, numWin)]
    end
    GBSCover = [GBSCoverChr GBSFragWinStatrs GBSFragWinEnds GBSFragCover]
    println(
        "INFO: Expected sequencing coverage based on $numFragSelected selected GBS fragments is approximately $(round(mean(GBSFragCover)*100;digits=2))%.",
    )
    ## module 4: writing out and plotting fragment-szie distribution and GBS coverage if required
    if plotOutput == true
        histogram(
            fragRaw[:, 4],
            normalize = :probability,
            title = "Size Distribution of Raw GBS Fragments",
            label = "",
            fmt = :png,
        )
        savefig("RawFragSize")
        histogram(
            fragGBS[:, 4],
            normalize = :probability,
            title = "Size Distribution of GBS Fragments (after Size-Selection)",
            label = "",
            fmt = :png,
        )
        savefig("GBSFragSize")
        histogram(
            GBSFragCover,
            normalize = :probability,
            title = "Coverage of GBS Fragments per Mb",
            label = "",
            fmt = :png,
        )
        savefig("GBSCoverage")
    end
    if writeOutput == true
        writedlm("RawFrag.txt", fragRaw)
        writedlm("GBSFrag.txt", fragGBS)
        writedlm("GBSCoverage.txt", GBSCover)
    end
end


text:
module SymTensors

using LinearAlgebra
using Random

import Base: convert, size, show, isless, isequal, ==, isapprox, *, conj, zero, inv
import Base: eltype, similar, copyto!
import Base: fill, fill!, rand, sum
import LinearAlgebra: mul!, rmul!, axpy!, axpby!, dot, norm, normalize!, svd, diag
import Base.intersect

#using DataStructures
include("sorteddict.jl")
export SortedDict, onlysemitokens

include("group.jl")
include("vspace.jl")
include("sector.jl")
#include("productspace.jl")
include("symtensor.jl")
include("trivialtensor.jl")
include("symreleg.jl")
include("contract.jl")
include("symtools.jl")

export +, *, zero, inv
export isequal, ==
export isapprox
export issimilar
export AbstractSymTensor, AbstractSymMatrix
export VectorSpace#, ProductSpace
export AbstractCharge
export Trivial
export U1Space, U1, U1Charge
export Z2Charge
export Sector
export SymTensor, SymMatrix, SymDiagonal, SymVector
export eye

export vtype

export SectorArray, SectorDiagonal
export _allsectorsandsizes
export _sectorisallowed

export array
export dim, dims, charges
export space
#export alldims, accdims, getdim, fulldims, signs
export rank, sectors, blocks
export fuse

#export _sectors_sortperm
export fuselegs
export splitleg
export permutelegs
export dropdummyleg
export trimleg!
export contract

export memoryrepr
export dual, isdual
export mapcharges
#export invlegs

export index_sector
export get_sector!
export set_sector!

export _svd_

export eltype, similar, copyto!
export mul!, rmul!, axpy!, axpby!, dot, norm
export normalize!
export fill, fill!
export fill_linearindex

export fermionswapgate

export isleftisometry, isrightisometry, isunitary

end

ext: jl

max_stars_repo_name: amiragha/mps

max_issues_repo_name: amiragha/mps

max_forks_repo_name: amiragha/mps


text:
float(::Type{ArbFloat})   = ArbFloat
float(::Type{ArbReal})    = ArbFloat
float(::Type{ArbComplex}) = ArbFloat

float(::Type{ArbFloat{P}})   where {P} = ArbFloat{P}
float(::Type{ArbReal{P}})    where {P} = ArbFloat{P}
float(::Type{ArbComplex{P}}) where {P} = ArbFloat{P}

float(x::ArbFloat{P})   where {P} = x
float(x::ArbReal{P})    where {P} = ArbFloat{P}(x)
float(z::ArbComplex{P}) where {P} = ArbFloat{P}(real(z))

real(::Type{ArbFloat})   = ArbFloat
real(::Type{ArbReal})    = ArbReal
real(::Type{ArbComplex}) = ArbReal

real(::Type{ArbFloat{P}})   where {P} = ArbFloat{P}
real(::Type{ArbReal{P}})    where {P} = ArbReal{P}
real(::Type{ArbComplex{P}}) where {P} = ArbReal{P}

real(x::ArbFloat{P})   where {P} = x
real(x::ArbReal{P})    where {P} = x
function real(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_get_real), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x)
    return z
end

imag(::Type{ArbFloat})   = ArbFloat
imag(::Type{ArbReal})    = ArbReal
imag(::Type{ArbComplex}) = ArbReal

imag(::Type{ArbFloat{P}})   where {P} = ArbFloat{P}
imag(::Type{ArbReal{P}})    where {P} = ArbReal{P}
imag(::Type{ArbComplex{P}}) where {P} = ArbReal{P}

imag(x::ArbFloat{P})   where {P} = x
imag(x::ArbReal{P})    where {P} = x
function imag(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_get_imag), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}), z, x)
    return z
end

complex(::Type{ArbFloat})   = ArbComplex
complex(::Type{ArbReal})    = ArbComplex
complex(::Type{ArbComplex}) = ArbComplex

complex(::Type{ArbFloat{P}})   where {P} = ArbComplex{P}
complex(::Type{ArbReal{P}})    where {P} = ArbComplex{P}
complex(::Type{ArbComplex{P}}) where {P} = ArbComplex{P}

complex(x::ArbFloat{P})   where {P} = ArbComplex{P}(x)
complex(x::ArbReal{P})    where {P} = ArbComplex{P}(x)
complex(z::ArbComplex{P}) where {P} = z

# other parts

angle(x::ArbFloat{P}) where {P} = zero(ArbFloat{P})
angle(x::ArbReal{P})  where {P} = zero(ArbReal{P})
function angle(x::ArbComplex{P}) where {P}
    z = ArbReal{P}()
    ccall(@libarb(acb_arg), Cvoid, (Ref{ArbReal}, Ref{ArbComplex}, Clong), z, x, P)
    return z
end

magnitude(x::T) where {T<:Complex} where {P} = hypot(reim(x)...)

magnitude(x::ArbFloat{P})   where {P} = x
magnitude(x::ArbReal{P})    where {P} = x
magnitude(x::ArbComplex{P}) where {P} = hypot(real(x), imag(x))


text:
struct WithArray
    a::Array{Float32,2}
end

struct WithCodeInstance
    a::Core.CodeInstance
end

struct WithDataType
    a::DataType
end

struct WithExpr
    a::Expr
end

struct WithString
    a::String
end

struct WithMethod
    a::Method
end

struct WithMethodInstance
    a::Core.MethodInstance
end

struct WithMethodTable
    a::Core.MethodTable
end

struct WithModule
    a::Module
end

struct WithSimpleVector
    a::Core.SimpleVector
end

struct WithSymbol
    a::Symbol
end

struct WithTask
    a::Task
end

struct WithTypeName
    a::Core.TypeName
end

struct WithTypeVar
    a::TypeVar
end

struct WithTypeMapEntry
    a::Core.TypeMapEntry
end

struct WithTypeMapLevel
    a::Core.TypeMapLevel
end

struct WithUnion
    a::Union
end

struct WithUnionAll
    a::UnionAll
end

@testset "Structs with builtin fields" begin
    @test begin
        b = JlrsReflect.reflect([WithArray])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithArray] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithArray")]
        pub struct WithArray<'frame, 'data> {
            pub a: ::jlrs::wrappers::ptr::ArrayRef<'frame, 'data>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithCodeInstance])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithCodeInstance] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithCodeInstance")]
        pub struct WithCodeInstance<'frame> {
            pub a: ::jlrs::wrappers::ptr::CodeInstanceRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithDataType])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithDataType] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithDataType")]
        pub struct WithDataType<'frame> {
            pub a: ::jlrs::wrappers::ptr::DataTypeRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithExpr])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithExpr] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithExpr")]
        pub struct WithExpr<'frame> {
            pub a: ::jlrs::wrappers::ptr::ExprRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithString])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithString] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithString")]
        pub struct WithString<'frame> {
            pub a: ::jlrs::wrappers::ptr::StringRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithMethod])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithMethod] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithMethod")]
        pub struct WithMethod<'frame> {
            pub a: ::jlrs::wrappers::ptr::MethodRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithMethodInstance])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithMethodInstance] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithMethodInstance")]
        pub struct WithMethodInstance<'frame> {
            pub a: ::jlrs::wrappers::ptr::MethodInstanceRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithMethodTable])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithMethodTable] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithMethodTable")]
        pub struct WithMethodTable<'frame> {
            pub a: ::jlrs::wrappers::ptr::MethodTableRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithModule])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithModule] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithModule")]
        pub struct WithModule<'frame> {
            pub a: ::jlrs::wrappers::ptr::ModuleRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithSimpleVector])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithSimpleVector] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithSimpleVector")]
        pub struct WithSimpleVector<'frame> {
            pub a: ::jlrs::wrappers::ptr::SimpleVectorRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithSymbol])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithSymbol] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithSymbol")]
        pub struct WithSymbol<'frame> {
            pub a: ::jlrs::wrappers::ptr::SymbolRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithTask])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithTask] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithTask")]
        pub struct WithTask<'frame> {
            pub a: ::jlrs::wrappers::ptr::TaskRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithTypeMapEntry])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithTypeMapEntry] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithTypeMapEntry")]
        pub struct WithTypeMapEntry<'frame> {
            pub a: ::jlrs::wrappers::ptr::TypeMapEntryRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithTypeMapLevel])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithTypeMapLevel] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithTypeMapLevel")]
        pub struct WithTypeMapLevel<'frame> {
            pub a: ::jlrs::wrappers::ptr::TypeMapLevelRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithTypeName])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithTypeName] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithTypeName")]
        pub struct WithTypeName<'frame> {
            pub a: ::jlrs::wrappers::ptr::TypeNameRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithTypeVar])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithTypeVar] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithTypeVar")]
        pub struct WithTypeVar<'frame> {
            pub a: ::jlrs::wrappers::ptr::TypeVarRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithUnion])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithUnion] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithUnion")]
        pub struct WithUnion<'frame> {
            pub a: ::jlrs::wrappers::ptr::UnionRef<'frame>,
        }"""
    end

    @test begin
        b = JlrsReflect.reflect([WithUnionAll])
        sb = JlrsReflect.StringWrappers(b)

        sb[WithUnionAll] === """#[repr(C)]
        #[derive(Clone, Debug, Unbox, ValidLayout, Typecheck)]
        #[jlrs(julia_type = "Main.WithUnionAll")]
        pub struct WithUnionAll<'frame> {
            pub a: ::jlrs::wrappers::ptr::UnionAllRef<'frame>,
        }"""
    end
end


text:
include("gs.jl")

using .GS

Lx = 2 # Comprimento em x
Ly = 1 # Comprimento em y
nx = 10 # Número de nós na direção x
ny = 5 # Número de nós na direção y
nivel = 4 # Nível de conectividade dos nós da ground structure
material = Material(1, 7e7)

estrutura = gerar_trelica_hiperconectada(Lx, Ly, nx, ny, material, nivel)
# Inserção de dados sobre as cargas e restrições dos nós.
# O nó mais próximo das coordenadas x, y receberá as propriedades
carregar_apoiar_no!(estrutura.nos, 0, 0; rx=true, ry=true)
carregar_apoiar_no!(estrutura.nos, 0, 1; rx=true, ry=true)
carregar_apoiar_no!(estrutura.nos, 2, 0; fy=-10)

vol = 0.0003

otm = Otimizador(vol, estrutura)
otimizar!(otm)
plotar_estrutura_otimizada(otm)
#plotar_estrutura_inicial(estrutura)

max_stars_repo_name: lorranfoliveira/ge

max_issues_repo_name: lorranfoliveira/ge

max_forks_repo_name: lorranfoliveira/ge


text:
using Distributed, LinearAlgebra
addprocs(2)
@everywhere using DistributedArrays, LinearAlgebra
##
function filldarray(a::AbstractArray)
    r = Array{Future, 1}(undef, nworkers())
    for (ipid, pid) in enumerate(workers())
        r[ipid] = @spawnat pid a
    end
    DArray(r)
end

function fillarray()
    r = [@spawnat pid myid()*ones(2) for pid in workers()]
    DArray(r)
end

function testclosure(d::DArray, a)
    function myclosure(a)
        localpart(d)[:] *= a
    end

    for(idx, pid) in enumerate(workers())
       remote_do(myclosure, pid, a)
    end
end

function fillarraywithnoarray()
    r = [@spawnat pid myid() for pid in workers()]
    DArray(r)
end

##
r = [@spawnat pid myid()*rand(100,100) for pid in workers()]
d = DArray(r, (200,100))

@everywhere fetchvalsvd(r::Future) = svd(fetch(r))

function dosvd(r::Array{Future, 1})
    s = Array{Any}(undef, nworkers())
    @sync for (i, rr) in enumerate(r)
        @async s[i] = remotecall_fetch(fetchvalsvd, rr.where, rr)
    end
    s
end

function dosvd(d::DArray)
    rr = Array{Future, 1}(undef, nworkers())
    @sync for (ipid, pid) in enumerate(procs(d))
        rr[ipid] = @spawnat pid svd(localpart(d))
    end
    fetch.(rr)
end

@everywhere fetchlocalsvd(d::DArray) = svd(localpart(d))

function dosvd2(d::DArray)
    s = Array{Any}(undef, nworkers())
    @sync for (ipid, pid) in enumerate(procs(d))
        @async s[ipid] = remotecall_fetch(fetchlocalsvd, pid, d)
    end
    s
end

##
function testnp(np::Int, sleeptime::Float64)
    rmprocs(workers())
    addprocs(np)
    @time @sync for (ipid, pid) in enumerate(workers())
        @async remotecall_fetch(sleep, pid, sleeptime)
    end
end


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

function _add_production_piecewise_linear_eqs!(
    model::JuMP.Model,
    g::Unit,
    formulation_prod_vars::Gar1962.ProdVars,
    formulation_pwl_costs::KnuOstWat2018.PwlCosts,
    formulation_status_vars::Gar1962.StatusVars,
)::Nothing
    eq_prod_above_def = _init(model, :eq_prod_above_def)
    eq_segprod_limit_a = _init(model, :eq_segprod_limit_a)
    eq_segprod_limit_b = _init(model, :eq_segprod_limit_b)
    eq_segprod_limit_c = _init(model, :eq_segprod_limit_c)
    segprod = model[:segprod]
    gn = g.name
    K = length(g.cost_segments)
    T = model[:instance].time

    # Gar1962.ProdVars
    prod_above = model[:prod_above]

    # Gar1962.StatusVars
    is_on = model[:is_on]
    switch_on = model[:switch_on]
    switch_off = model[:switch_off]

    for t in 1:T
        for k in 1:K
            # Pbar^{k-1)
            Pbar0 =
                g.min_power[t] +
                (k > 1 ? sum(g.cost_segments[ell].mw[t] for ell in 1:k-1) : 0.0)
            # Pbar^k
            Pbar1 = g.cost_segments[k].mw[t] + Pbar0

            Cv = 0.0
            SU = g.startup_limit   # startup rate
            if Pbar1 <= SU
                Cv = 0.0
            elseif Pbar0 < SU # && Pbar1 > SU
                Cv = Pbar1 - SU
            else # Pbar0 >= SU
                # this will imply that we cannot produce along this segment if
                # switch_on = 1
                Cv = g.cost_segments[k].mw[t]
            end
            Cw = 0.0
            SD = g.shutdown_limit  # shutdown rate
            if Pbar1 <= SD
                Cw = 0.0
            elseif Pbar0 < SD # && Pbar1 > SD
                Cw = Pbar1 - SD
            else # Pbar0 >= SD
                Cw = g.cost_segments[k].mw[t]
            end

            if g.min_uptime > 1
                # Equation (46) in Kneuven et al. (2020)
                eq_segprod_limit_a[gn, t, k] = @constraint(
                    model,
                    segprod[gn, t, k] <=
                    g.cost_segments[k].mw[t] * is_on[gn, t] -
                    Cv * switch_on[gn, t] -
                    (t < T ? Cw * switch_off[gn, t+1] : 0.0)
                )
            else
                # Equation (47a)/(48a) in Kneuven et al. (2020)
                eq_segprod_limit_b[gn, t, k] = @constraint(
                    model,
                    segprod[gn, t, k] <=
                    g.cost_segments[k].mw[t] * is_on[gn, t] -
                    Cv * switch_on[gn, t] -
                    (t < T ? max(0, Cv - Cw) * switch_off[gn, t+1] : 0.0)
                )

                # Equation (47b)/(48b) in Kneuven et al. (2020)
                eq_segprod_limit_c[gn, t, k] = @constraint(
                    model,
                    segprod[gn, t, k] <=
                    g.cost_segments[k].mw[t] * is_on[gn, t] -
                    max(0, Cw - Cv) * switch_on[gn, t] -
                    (t < T ? Cw * switch_off[gn, t+1] : 0.0)
                )
            end

            # Definition of production
            # Equation (43) in Kneuven et al. (2020)
            eq_prod_above_def[gn, t] = @constraint(
                model,
                prod_above[gn, t] == sum(segprod[gn, t, k] for k in 1:K)
            )

            # Objective function
            # Equation (44) in Kneuven et al. (2020)
            add_to_expression!(
                model[:obj],
                segprod[gn, t, k],
                g.cost_segments[k].cost[t],
            )

            # Also add an explicit upper bound on segprod to make the solver's
            # work a bit easier
            set_upper_bound(segprod[gn, t, k], g.cost_segments[k].mw[t])
        end
    end
end


text:
module Mocking

using Compat: mergewith
using ExprTools: splitdef, combinedef

export @patch, @mock, Patch, apply

include("expr.jl")
include("dispatch.jl")
include("patch.jl")
include("mock.jl")
include("deprecated.jl")

# Create the initial definition of `activated` which defaults Mocking to be disabled
activated() = false

"""
    Mocking.activate()

Enable `@mock` call sites to allow for calling patches instead of the original function.
"""
function activate()
    # Avoid redefining `activated` when it's already set appropriately
    !activated() && @eval activated() = true
    return nothing
end

"""
    Mocking.deactivate()

Disable `@mock` call sites to only call the original function.
"""
function deactivate()
    # Avoid redefining `activated` when it's already set appropriately
    activated() && @eval activated() = false
    return nothing
end

const NULLIFIED = Ref{Bool}(false)

"""
    Mocking.nullify()

Force any packages loaded after this point to treat the `@mock` macro as a no-op. Doing so
will maximize performance by eliminating any runtime checks taking place at the `@mock` call
sites but will break any tests that require patches to be applied.

Note to ensure that all `@mock` macros are inoperative be sure to call this function before
loading any packages which depend on Mocking.jl.
"""
function nullify()
    global NULLIFIED[] = true
    return nothing
end

end # module


text:
#!/usr/bin/julia

# Trizen
# 27 April 2017
# https://github.com/trizen

# Complex transform of an image, by mapping each pixel position to a complex function.

# usage:
#   julia complex_transform.jl [image]

using Images
#using SpecialFunctions

function map_val(value, in_min, in_max, out_min, out_max)
    (value - in_min) * (out_max - out_min) / (in_max - in_min) + out_min
end

function complex_transform(file)

    img = load(file)
    height, width = size(img)

    function transform(x, y)

        z = complex(
            (2 * x - width ) / width,
            (2 * y - height) / height
        )

        # Complex function
        z = (cos(z) + sin(z)*im) / im

        (real(z), imag(z))
    end

    matrix = zeros(height, width, 2)

    min_x, min_y = (Inf, Inf)
    max_x, max_y = (-Inf, -Inf)

    for y in 1:height, x in 1:width
        new_x, new_y = transform(x, y)

        matrix[y,x,1] = new_x
        matrix[y,x,2] = new_y

        if (new_x < min_x)
            min_x = new_x
        end
        if (new_y < min_y)
            min_y = new_y
        end
        if (new_x > max_x)
            max_x = new_x
        end
        if (new_y > max_y)
            max_y = new_y
        end
    end

    println("X: [$min_x, $max_x]")
    println("Y: [$min_y, $max_y]")

    out_img = zeros(RGB{N0f8}, height, width)

    for y in 1:height, x in 1:width
        out_img[y,x] = RGB{N0f8}(0,0,0)
    end

    for y in 1:height, x in 1:width

        new_x = map_val(matrix[y,x,1], min_x, max_x, 1, width)
        new_y = map_val(matrix[y,x,2], min_y, max_y, 1, height)

        if (abs(new_x) == Inf || isnan(new_x) || abs(new_y) == Inf || isnan(new_y))
            println("Skipping one pixel...")
            continue
        end

        new_x = round(Int64, new_x)
        new_y = round(Int64, new_y)

        out_img[new_y,new_x] = img[y,x]
    end

    return out_img
end

inputfile = length(ARGS) > 0 ? ARGS[1] : "input.png"
save("complex_transform.png", complex_transform(inputfile))

max_stars_repo_name: trizen/julia-scripts

max_issues_repo_name: trizen/julia-scripts

max_forks_repo_name: trizen/julia-scripts


text:
using cascadeddecay, QuantumOptics, meanfield
using PyCall
@pyimport matplotlib.pyplot as plt

const γ = 1.
const e_dipole = [0,0,1.]
const T = [0:0.1:2.5]
const D = exp(range(log(0.1), stop=log(100), length=100))

plt.figure("trace_distance", figsize=(6,4))
plt.figure("sz_distance", figsize=(6,4))

for atomnumber=2:6
    meanfield_error = Float64[]
    correlation_error = Float64[]
    meanfield_sz_error = Float64[]
    correlation_sz_error = Float64[]

    for d=D
        println(d)
        #atoms = Atom[Atom(position, e_dipole) for position=cube(0.627)]
        #atoms = Atom[Atom(position, e_dipole) for position=cube(d)]
        atoms = Atom[Atom(position, e_dipole) for position=chain(d,atomnumber)]
        system = CascadedDecaySystem(atoms, γ)
        basis_atom = atoms[1].basis

        # System operators
        H = Hint(system)
        J = JumpOperators(system)
        Jdagger =map(dagger,J)

        # Non-hermitian Hamiltonian
        Hnh = H - 0.5im*sum([Jdagger[i]*J[i] for i=1:length(atoms)])
        Hnh_dagger = dagger(Hnh)

        # Sparse system operators
        J_sparse = map(operators_sparse.SparseOperator, J)
        Jdagger_sparse = map(operators_sparse.SparseOperator, Jdagger)

        Hnh_sparse = operators_sparse.SparseOperator(Hnh)
        Hnh_dagger_sparse = operators_sparse.SparseOperator(Hnh_dagger)

        # Initial conditions
        I = identityoperator(basis_atom)
        R = 1/sqrt(2)*(I + 1im*sigmay)

        Ψ₀_atom = R*basis_ket(basis_atom,1)
        Ψ₀ = reduce(⊗, [Ψ₀_atom for i=1:length(atoms)])
        ρ₀ = Ψ₀⊗dagger(Ψ₀)

        # Meanfield timeevolution
        N = length(system.atoms)
        sz0 = zeros(ComplexF64, N)
        sp0 = zeros(ComplexF64, N) .- 0.5

        T_mf, sz_mf, sp_mf = meanfield_timeevolution(T, system, sz0, sp0)

        sm_mf = conj(sp_mf)
        sx_mf = real(sp_mf + sm_mf)
        sy_mf = real(1im*(sp_mf - sm_mf))

        # Correlation included timeevolution
        N = length(system.atoms)
        sz0 = zeros(ComplexF64, N)
        sp0 = zeros(ComplexF64, N) .- 0.5
        Cpm0 = zeros(ComplexF64, N, N)
        Cpz0 = zeros(ComplexF64, N, N)
        Cpp0 = zeros(ComplexF64, N, N)
        Czz0 = zeros(ComplexF64, N, N)

        T_cor, sz_cor, sp_cor, C_cor = correlation_timeevolution(T, system, sz0, sp0, Cpm0, Cpz0, Cpp0, Czz0)

        sm_cor = conj(sp_cor)
        sx_cor = real(sp_cor + sm_cor)
        sy_cor = real(1im*(sp_cor - sm_cor))

        # Master time evolution
        tout, ρ_t = timeevolution.master_nh(T, ρ₀, Hnh_sparse, J_sparse)

        error_mf = 0
        error_cor = 0
        error_sz_mf = 0
        error_sz_cor = 0

        for i=1:length(tout)
            ρmf = meanfield.densityoperator(sx_mf[i], sy_mf[i], real(sz_mf[i]))
            ρcor = meanfield.densityoperator(sx_cor[i], sy_cor[i], real(sz_cor[i]), C_cor[i])
            exp_sz = real(expect(embed(system.basis, 1, sigmaz), ρ_t[i]))

            error_sz_mf = max(error_sz_mf, abs(exp_sz - real(sz_mf[i][1])))
            error_sz_cor = max(error_sz_cor, abs(exp_sz - real(sz_cor[i][1])))
            error_mf = max(error_mf, tracedistance(ρ_t[i] ,ρmf))
            error_cor = max(error_cor, tracedistance(ρ_t[i] ,ρcor))
        end
        push!(meanfield_error, error_mf)
        push!(correlation_error, error_cor)
        push!(meanfield_sz_error, error_sz_mf)
        push!(correlation_sz_error, error_sz_cor)
    end

    plt.figure("trace_distance")
    plt.loglog(D, meanfield_error, "b", alpha=0.1*atomnumber, label="Meanfield ($atomnumber)")
    plt.loglog(D, correlation_error, "g", alpha=0.1*atomnumber, label="Correlations ($atomnumber)")

    plt.figure("sz_distance")
    plt.loglog(D, meanfield_sz_error, "b", alpha=0.1*atomnumber, label="Meanfield ($atomnumber)")
    plt.loglog(D, correlation_sz_error, "g", alpha=0.1*atomnumber, label="Correlations ($atomnumber)")
end

plt.figure("trace_distance")
plt.ylim(0, 1)
plt.xlabel("\$d/\\lambda\$")
plt.ylabel("Error (Trace distance)")
#plt.legend()
plt.savefig("error_analysis_atomnumber_chain6.pdf")

plt.figure("sz_distance")
plt.ylim(0, 1)
plt.xlabel("\$d/\\lambda\$")
plt.ylabel("Error (sz)")
#plt.legend()
plt.savefig("error_analysis_atomnumber_sz_chain6.pdf")
plt.show()

# # Select only first atom
# first(xvec) = Float64[real(x[1]) for x=xvec]

# name = "cube_atoms=8_d=015"
# # name = "chain_atoms=8_d=07"

# plt.figure(figsize=(6,4))
# plt.plot(T, -sz_full, lw=3.0, color="darkorange", label="Collective")
# plt.plot(T, exp(-γ*T)-1, lw=3.0, color="gray", label="Independent")
# #plt.plot(T_smf, -sz_smf, label="symmetric meanfield")
# plt.plot(T, -first(sz_mf), "-", lw=1.5, color="navy", label="Mean-field")
# plt.plot(T, -first(sz_cor), "-", lw=1.5, color="green", label="Correlations")
# plt.ylabel("\$-\\langle\\sigma_z\\rangle\$")
# plt.ylim(-1,1)
# plt.legend()
# plt.savefig("$(name)_sigmaz_timeevolution.pdf")

# plt.figure(figsize=(6,4))
# plt.plot(T, sx_full, lw=3.0, color="darkorange", label="Collective")
# plt.plot(T, -exp(-0.5*γ*T), lw=3.0, color="gray", label="Independent")
# #plt.plot(T_smf, sx_smf, label="symmetric meanfield")
# plt.plot(T, first(sx_mf), "-", lw=1.5, color="navy", label="Mean-field")
# plt.plot(T, first(sx_cor), "-", lw=1.5, color="green", label="Correlations")
# plt.ylabel("\$\\langle\\sigma_x\\rangle\$")
# plt.ylim(-1,1)
# plt.legend()
# plt.savefig("$(name)_sigmax_timeevolution.pdf")

# plt.figure(figsize=(6,4))
# plt.plot(T, sy_full, lw=3.0, color="darkorange", label="Collective")
# plt.plot(T, 0.*T, lw=3.0, color="gray", label="Independent")
# #plt.plot(T_smf, sy_smf, label="symmetric meanfield")
# plt.plot(T, first(sy_mf), "-", lw=1.5, color="navy", label="Mean-field")
# plt.plot(T, first(sy_cor), "-", lw=1.5, color="green", label="Correlations")
# plt.ylabel("\$\\langle\\sigma_y\\rangle\$")
# plt.ylim(-1,1)
# plt.legend()
# plt.savefig("$(name)_sigmay_timeevolution.pdf")

# plt.figure(figsize=(6,4))
# plt.plot(T, error_mf, lw=1.2, color="navy", label="Mean-field")
# plt.plot(T, error_cor, lw=1.2, color="green", label="Correlation")
# plt.xlabel("Time [\$\\gamma^{-1}\$]")
# plt.ylabel("Error (Trace distance)")
# plt.ylim(0,1)
# plt.legend()
# plt.savefig("$(name)_error_timeevolution.pdf")

# plt.show()
#println("⦑H⟩ = ", trunc(real(exp), 3))
#println("⦑H⟩ = ", trunc(real(exp_mean), 3))


text:
try
    using RCall
catch
end

# same as May 21, but use different target, namely posterior expectation

using Distributions
using StatsBase
using Iterators
using Gurobi
using JLD

using EmpiricalBayes

nreps =25
k_max = 16 #in total 25*16 replications

true_dists = [MixtureModel([ Normal(-0.3,.2), Normal(0,0.9)],[0.8, 0.2])]
xs = collect(-4.5:0.5:4.5)
minimax_σs = [0.2]
ms = [20_000]

combs = []

for (true_dist, x, σ, m) in product(true_dists, xs, minimax_σs, ms)
    push!(combs, Dict(:dist=>true_dist, :x=>x, :σ=>σ, :m=>m) )
end

i = parse(Int64, ARGS[1])

comb = combs[i]

m = comb[:m]
σ = comb[:σ]
true_dist = comb[:dist]
x = comb[:x]

marginal_grid = collect(range(-6, stop=6, length=1001));
prior_grid = collect(range(-3, stop=3, length=121));
marginal_h = marginal_grid[2]-marginal_grid[1];

d_true = NormalConvolutionProblem(true_dist, marginal_grid)
ds = MixingNormalConvolutionProblem(Normal, σ, prior_grid, marginal_grid);

target = PosteriorTarget(PosteriorMeanNumerator(x))

dry_run = CEB_ci(rand(d_true,Int(m/2)), rand(d_true, Int(m/2)), ds, target)

for k=1:k_max
    sim_array = []
    for l = 1:nreps
        @show l
        Xs = rand(d_true,m)
        brad_non_null =  fit(BradDeconvolveR, Xs; deltaAt = :nothing)
        brad_null =  fit(BradDeconvolveR, Xs)

        n_total = length(Xs)
        n_half = ceil(Int, n_total/2)
        idx_test = sample(1:n_total, n_half, replace=false)
        idx_train = setdiff(1:n_total, idx_test)
        Xs_train = Xs[idx_train]
        Xs_test = Xs[idx_test]

        ceb_res = CEB_ci(Xs_train, Xs_test, ds, target, dry_run[2], dry_run[3], dry_run[4].C_bias)
        push!(sim_array, (ceb_res, brad_null, brad_non_null))
    end
    res = (i, comb, sim_array)
    save("/scratch/users/ignat/sims/May21b/mysim_$(i)_$(k).jld", "res", res)
end


text:
@testset "VCF" begin
    metainfo = VCF.MetaInfo()
    @test !isfilled(metainfo)
    @test occursin(r"^GeneticVariation.VCF.MetaInfo: <not filled>", repr(metainfo))
    @test_throws ArgumentError metainfotag(metainfo)

    metainfo = VCF.MetaInfo(Vector{UInt8}("##source=foobar1234"))
    @test isfilled(metainfo)
    @test metainfotag(metainfo) == "source"
    @test metainfoval(metainfo) == "foobar1234"

    metainfo = VCF.MetaInfo("##source=foobar1234")
    @test isfilled(metainfo)
    @test metainfotag(metainfo) == "source"
    @test metainfoval(metainfo) == "foobar1234"

    metainfo = VCF.MetaInfo(metainfo)
    @test isa(metainfo, VCF.MetaInfo)
    metainfo = VCF.MetaInfo(metainfo, tag="date")
    @test metainfotag(metainfo) == "date"
    metainfo = VCF.MetaInfo(metainfo, value="2017-01-30")
    @test metainfoval(metainfo) == "2017-01-30"
    metainfo = VCF.MetaInfo(metainfo, tag="INFO", value=["ID"=>"DP", "Number"=>"1", "Type"=>"Integer", "Description"=>"Total Depth"])
    @test metainfo["ID"] == "DP"
    @test metainfo["Number"] == "1"
    @test metainfo["Type"] == "Integer"
    @test metainfo["Description"] == "Total Depth"
    @test metainfotag(metainfo) == "INFO"
    @test metainfoval(metainfo) == """<ID=DP,Number=1,Type=Integer,Description="Total Depth">"""

    record = VCF.Record()
    @test !isfilled(record)
    @test occursin(r"^GeneticVariation.VCF.Record: <not filled>", repr(record))
    @test_throws ArgumentError VCF.chrom(record)

    record = VCF.Record("20\t302\t.\tT\tTA\t999\t.\t.\tGT")
    @test isfilled(record)
    @test VCF.haschrom(record)
    @test VCF.chrom(record) == "20"
    @test VCF.haspos(record)
    @test VCF.pos(record) == 302
    @test !VCF.hasid(record)
    @test_throws MissingFieldException VCF.id(record)
    @test VCF.hasref(record)
    @test VCF.ref(record) == "T"
    @test VCF.hasalt(record)
    @test VCF.alt(record) == ["TA"]
    @test VCF.hasqual(record)
    @test VCF.qual(record) == 999
    @test !VCF.hasfilter(record)
    @test_throws MissingFieldException VCF.filter(record)
    @test VCF.infokeys(record) == String[]
    @test !VCF.hasinfo(record)
    @test_throws MissingFieldException VCF.info(record)
    @test VCF.hasformat(record)
    @test VCF.format(record) == ["GT"]

    # empty data is not a valid VCF record
    @test_throws ArgumentError VCF.Record("")
    @test_throws ArgumentError VCF.Record(Vector{UInt8}(""))

    record = VCF.Record(Vector{UInt8}(".\t.\t.\t.\t.\t.\t.\t.\t"))
    @test isfilled(record)
    @test !VCF.haschrom(record)
    @test !VCF.haspos(record)
    @test !VCF.hasid(record)
    @test !VCF.hasref(record)
    @test !VCF.hasalt(record)
    @test !VCF.hasqual(record)
    @test !VCF.hasfilter(record)
    @test !VCF.hasinfo(record)

    record = VCF.Record(record)
    @test isa(record, VCF.Record)
    record = VCF.Record(record, chrom="chr1")
    @test VCF.chrom(record) == "chr1"
    record = VCF.Record(record, pos=1234)
    @test VCF.pos(record) == 1234
    record = VCF.Record(record, id="rs1111")
    @test VCF.id(record) == ["rs1111"]
    record = VCF.Record(record, ref="A")
    @test VCF.ref(record) == "A"
    record = VCF.Record(record, alt=["AT"])
    @test VCF.alt(record) == ["AT"]
    record = VCF.Record(record, qual=11.2)
    @test VCF.qual(record) == 11.2
    record = VCF.Record(record, filter="PASS")
    @test VCF.filter(record) == ["PASS"]
    record = VCF.Record(record, info=Dict("DP" => 20, "AA" => "AT", "DB"=>nothing))
    @test VCF.info(record, "DP") == "20"
    @test VCF.info(record, "AA") == "AT"
    @test VCF.info(record, "DB") == ""
    @test VCF.infokeys(record) == ["DP", "AA", "DB"]
    record = VCF.Record(record, genotype=[Dict("GT" => "0/0", "DP" => [10,20])])
    @test VCF.format(record) == ["DP", "GT"]
    @test VCF.genotype(record) == [["10,20", "0/0"]]

    let header = VCF.Header()
        @test isempty(header)
        push!(header, "##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta")
        @test !isempty(header)
        @test length(header) == 1
        pushfirst!(header, "##fileformat=VCFv4.3")
        @test length(header) == 2
        @test collect(header) == [
            VCF.MetaInfo("##fileformat=VCFv4.3"),
            VCF.MetaInfo("##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta")]
        @test startswith(repr(header), "GeneticVariation.VCF.Header:")
    end

    let header = VCF.Header(["##fileformat=VCFv4.3"], ["Sample1"])
        @test !isempty(header)
        @test length(header) == 1
        @test header.sampleID == ["Sample1"]
        @test first(header) == VCF.MetaInfo("##fileformat=VCFv4.3")
    end

    # minimum header
    data = Vector{UInt8}("""
    ##fileformat=VCFv4.3
    #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO
    """)
    reader = VCF.Reader(BufferedInputStream(data))
    @test isa(header(reader), VCF.Header)
    let header = header(reader)
        @test length(header.metainfo) == 1
        @test metainfotag(header.metainfo[1]) == "fileformat"
        @test metainfoval(header.metainfo[1]) == "VCFv4.3"
        @test isempty(header.sampleID)
    end

    # realistic header
    data = Vector{UInt8}("""
    ##fileformat=VCFv4.2
    ##fileDate=20090805
    ##source=myImputationProgramV3.1
    ##reference=file:///seq/references/1000GenomesPilot-NCBI36.fasta
    ##contig=<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>
    ##phasing=partial
    ##INFO=<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
    ##INFO=<ID=AA,Number=1,Type=String,Description="Ancestral Allele">
    #CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	NA00001	NA00002	NA00003
    """)
    reader = VCF.Reader(BufferedInputStream(data))
    @test isa(header(reader), VCF.Header)

    let header = header(reader)
        @test length(header.metainfo) == 10

        let metainfo = header.metainfo[1]
            @test metainfotag(metainfo) == "fileformat"
            @test metainfoval(metainfo) == "VCFv4.2"
            @test_throws ArgumentError keys(metainfo)
            @test_throws ArgumentError values(metainfo)
        end
        @test length(findall(header, "fileformat")) == 1
        @test first(findall(header, "fileformat")) == header.metainfo[1]

        let metainfo = header.metainfo[2]
            @test metainfotag(metainfo) == "fileDate"
            @test metainfoval(metainfo) == "20090805"
            @test_throws ArgumentError keys(metainfo)
            @test_throws ArgumentError values(metainfo)
        end

        let metainfo = header.metainfo[5]
            @test metainfotag(metainfo) == "contig"
            @test metainfoval(metainfo) == """<ID=20,length=62435964,assembly=B36,md5=f126cdf8a6e0c7f379d618ff66beb2da,species="Homo sapiens",taxonomy=x>"""
            @test keys(metainfo) == ["ID", "length", "assembly", "md5", "species", "taxonomy"]
            @test values(metainfo) == ["20", "62435964", "B36", "f126cdf8a6e0c7f379d618ff66beb2da", "Homo sapiens", "x"]
            @test metainfo["ID"] == "20"
            @test metainfo["md5"] == "f126cdf8a6e0c7f379d618ff66beb2da"
            @test metainfo["taxonomy"] == "x"
        end

        let metainfo = header.metainfo[7]
            @test metainfotag(metainfo) == "INFO"
            @test metainfoval(metainfo) == """<ID=NS,Number=1,Type=Integer,Description="Number of Samples With Data">"""
            @test keys(metainfo) == ["ID", "Number", "Type", "Description"]
            @test values(metainfo) == ["NS", "1", "Integer", "Number of Samples With Data"]
            @test metainfo["ID"] == "NS"
            @test metainfo["Type"] == "Integer"
        end
        @test length(findall(header, "INFO")) == 4

        @test header.sampleID == ["NA00001", "NA00002", "NA00003"]
    end

    data = Vector{UInt8}("""
    ##fileformat=VCFv4.3
    ##contig=<ID=chr1>
    ##contig=<ID=chr2>
    ##INFO=<ID=DP,Number=1,Type=Integer,Description="Total Depth">
    ##INFO=<ID=AF,Number=A,Type=Float,Description="Allele Frequency">
    ##FORMAT=<ID=GT,Number=1,Description="Genotype">
    #CHROM\tPOS\tID\tREF\tALT\tQUAL\tFILTER\tINFO\tFORMAT\tNA00001\tNA00002
    chr1\t1234\trs001234\tA\tC\t30\tPASS\tDP=10;AF=0.3\tGT\t0|0\t0/1
    chr2\t4\t.\tA\tAA,AAT\t.\t.\tDP=5\tGT:DP\t0|1:42\t0/1
    """)
    reader = VCF.Reader(BufferedInputStream(data))
    record = VCF.Record()

    @test read!(reader, record) === record
    @test VCF.chrom(record) == "chr1"
    @test VCF.pos(record) === 1234
    @test VCF.id(record) == ["rs001234"]
    @test VCF.ref(record) == "A"
    @test VCF.alt(record) == ["C"]
    @test VCF.qual(record) === 30.0
    @test VCF.filter(record) == ["PASS"]
    @test VCF.info(record) == ["DP" => "10", "AF" => "0.3"]
    @test VCF.info(record, "DP") == "10"
    @test VCF.info(record, "AF") == "0.3"
    @test VCF.format(record) == ["GT"]
    @test VCF.genotype(record) == [["0|0"], ["0/1"]]
    @test VCF.genotype(record, 1) == VCF.genotype(record)[1]
    @test VCF.genotype(record, 2) == VCF.genotype(record)[2]
    @test VCF.genotype(record, 1, "GT") == "0|0"
    @test VCF.genotype(record, 2, "GT") == "0/1"
    @test VCF.genotype(record, 1:2, "GT") == ["0|0", "0/1"]
    @test VCF.genotype(record, :, "GT") == VCF.genotype(record, 1:2, "GT")
    @test occursin(r"^GeneticVariation.VCF.Record:\n.*", repr(record))

    @test read!(reader, record) === record
    @test VCF.chrom(record) == "chr2"
    @test VCF.pos(record) == 4
    @test !VCF.hasid(record)
    @test VCF.ref(record) == "A"
    @test VCF.alt(record) == ["AA", "AAT"]
    @test !VCF.hasqual(record)
    @test !VCF.hasfilter(record)
    @test VCF.info(record) == ["DP" => "5"]
    @test VCF.info(record, "DP") == "5"
    @test_throws KeyError VCF.info(record, "AF")
    @test VCF.format(record) == ["GT", "DP"]
    @test VCF.genotype(record) == [["0|1", "42"], ["0/1", "."]]
    @test VCF.genotype(record, 1) == VCF.genotype(record)[1]
    @test VCF.genotype(record, 2) == VCF.genotype(record)[2]
    @test VCF.genotype(record, 1, "GT") == "0|1"
    @test VCF.genotype(record, 1, "DP") == "42"
    @test VCF.genotype(record, 2, "GT") == "0/1"
    @test VCF.genotype(record, 2, "DP") == "."
    @test VCF.genotype(record, 1:2, "GT") == ["0|1", "0/1"]
    @test VCF.genotype(record, 1:2, "DP") == ["42", "."]
    @test VCF.genotype(record, :, "DP") == VCF.genotype(record, 1:2, "DP")
    @test_throws KeyError VCF.genotype(record, :, "BAD")

    @test_throws EOFError read!(reader, record)

    # round-trip test
    vcfdir = joinpath(fmtdir, "VCF")
    for specimen in YAML.load_file(joinpath(vcfdir, "index.yml"))
        filepath = joinpath(vcfdir, specimen["filename"])
        records = VCF.Record[]
        reader = open(VCF.Reader, filepath)
        output = IOBuffer()
        writer = VCF.Writer(output, header(reader))
        for record in reader
            write(writer, record)
            push!(records, record)
        end
        close(reader)
        flush(writer)

        records2 = VCF.Record[]
        for record in VCF.Reader(IOBuffer(take!(output)))
            push!(records2, record)
        end
        @test records == records2
    end
end

function parsehex(str)
    return map(x -> parse(UInt8, x, base = 16), split(str, ' '))
end


text:
###########################################################################
# Scaled Least-squares
###########################################################################
function bearing_sls(angles, gains, ref_angles, ref_gains)
	@assert length(angles) == length(gains)
	num_gains = length(gains)
	best_err = Inf
	best_shift = 0
	best_scale = 0
	temp_gains = copy(gains)
	for scale_factor = 0:100
		temp_gains += 1
		for shift = 0:359
			err = 0.0
			for j = 1:num_gains
				if validgain(gains[j])
					ref_gain = interp_linear(angles[j], shift, ref_angles, ref_gains)
					err += (temp_gains[j] - ref_gain)^2
				end
			end
			if err < best_err
				best_err = err
				best_shift = shift
				best_scale = scale_factor
			end
		end
	end
	return 360 - best_shift
end

function bearing_sls(dir_file::AbstractString, ref_file::AbstractString; samples::Int64=0, sense::Int64=1)
	test_gp = GainPattern(dir_file)
	ref_gp = GainPattern(ref_file)

	num_samples = length(test_gp.angles)
	num_exp = length(test_gp.samples[1])

	# Allowing user to specify count
	if (samples > 0) && (samples < num_samples)
		num_samples = samples
	end

	aoa_array = zeros(num_exp)
	angles = zeros(num_samples)
	gains = zeros(num_samples)
	for i = 1:num_exp
		if sense > 0
			for j = 1:num_samples
				# Loop through all samples in this experiment...
				angles[j] = test_gp.angles[j]
				gains[j] = test_gp.samples[j][i]
			end
		elseif sense == -1
			# backwards loop, starting from last measurement
			for j = 1:num_samples
				angles[j] = test_gp.angles[end - j + 1]
				gains[j] = test_gp.samples[end - j + 1][i]
			end
		elseif sense == -2
			# backwards loop, starting from origin
			for j = 2:num_samples
				angles[j] = test_gp.angles[end - j + 2]
				gains[j] = test_gp.samples[end - j + 2][i]
			end
			angles[1] = test_gp.angles[1]
			gains[1] = test_gp.samples[1][i]
		end
		aoa_array[i] = bearing_sls(angles, gains, ref_gp.angles, ref_gp.meangains)
	end
	return aoa_array
end


text:
using Qutilities
using Test

using LinearAlgebra: diagm, I, tr

@testset "Qutilities" begin

@testset "sigma_x, sigma_y, sigma_z" begin
    @test sigma_x^2 == I
    @test sigma_y^2 == I
    @test sigma_z^2 == I
    @test -im*sigma_x*sigma_y*sigma_z == I
end

@testset "ptrace, ptranspose" begin
    let A = I(2),
        B = reshape(1:16, 4, 4),
        C = [[0, 0, 1] [0, 1, 0] [1, 0, 0]],
        ABC = kron(A, B, C),
        dims = (2, 4, 3)

        @test ptrace(ABC, dims, 1) == tr(A) * kron(B, C)
        @test ptranspose(ABC, dims, 1) == kron(A', B, C)
        @test ptrace(ABC, dims, 2) == tr(B) * kron(A, C)
        @test ptranspose(ABC, dims, 2) == kron(A, B', C)
        @test ptrace(ABC, dims, 3) == tr(C) * kron(A, B)
        @test ptranspose(ABC, dims, 3) == kron(A, B, C')
    end

    let M = reshape(1.0:16.0, 4, 4),
        MT1 = [[1.0, 2.0, 9.0, 10.0] [5.0, 6.0, 13.0, 14.0] [3.0, 4.0, 11.0, 12.0] [7.0, 8.0, 15.0, 16.0]],
        MT2 = [[1.0, 5.0, 3.0, 7.0] [2.0, 6.0, 4.0, 8.0] [9.0, 13.0, 11.0, 15.0] [10.0, 14.0, 12.0, 16.0]]

        @test ptrace(M, (1, 4), 1) == M
        @test ptranspose(M, (1, 4), 1) == M
        @test ptrace(M, (1, 4), 2) == fill(tr(M), 1, 1)
        @test ptranspose(M, (1, 4), 2) == transpose(M)
        @test ptrace(M, (2, 2), 1) == [[12.0, 14.0] [20.0, 22.0]]
        @test ptranspose(M, (2, 2), 1) == MT1
        @test ptrace(M, (2, 2), 2) == [[7.0, 11.0] [23.0, 27.0]]
        @test ptranspose(M, (2, 2), 2) == MT2
        @test ptrace(M, (4, 1), 1) == fill(tr(M), 1, 1)
        @test ptranspose(M, (4, 1), 1) == transpose(M)
        @test ptrace(M, (4, 1), 2) == M
        @test ptranspose(M, (4, 1), 2) == M

        @test ptrace(M, 1) == ptrace(M, (2, 2), 1)
        @test ptrace(M, 2) == ptrace(M, (2, 2), 2)
        @test ptranspose(M, 1) == ptranspose(M, (2, 2), 1)
        @test ptranspose(M, 2) == ptranspose(M, (2, 2), 2)

        @test ptrace(M) == ptrace(M, (2, 2), 2)
        @test ptranspose(M) == ptranspose(M, (2, 2), 2)
    end

    let M = [[1.0, im] [im, 1.0]]

        @test ptrace(M, (1, 2), 1) == M
        @test ptranspose(M, (1, 2), 1) == M
        @test ptrace(M, (1, 2), 2) == fill(tr(M), 1, 1)
        @test ptranspose(M, (1, 2), 2) == transpose(M)
    end
end

@testset "binent" begin
    @test binent(0.0) == 0.0
    @test binent(0.5) == 1.0
    @test binent(1.0) == 0.0
end

@testset "purity, S_vn, S_renyi" begin
    let rho1 = I(4) / 4.0,
        rho2 = [[2.0, im] [-im, 2.0]] / 2.0,
        eigs2 = [1.0, 3.0] / 2.0

        @test purity(rho1) == 0.25
        @test purity(rho2) == sum(eigs2.^2)

        @test S_renyi(rho1, 0) == 2.0
        @test S_renyi(rho2, 0) == 1.0

        @test S_vn(rho1) == 2.0
        @test S_vn(rho2) == -sum(eigs2 .* log2.(eigs2))

        @test S_renyi(rho1) == 2.0
        @test S_renyi(rho2) == -log2(sum(eigs2.^2))

        @test S_renyi(rho1, Inf) == 2.0
        @test S_renyi(rho2, Inf) == -log2(maximum(eigs2))
    end
end

@testset "mutinf" begin
    let rho = diagm([3, 2, 1, 2]) / 8.0

        @test isapprox(mutinf(rho), (1.5 - 5.0 * log2(5.0) / 8.0))
        @test isapprox(mutinf(rho, S_renyi), (1.0 + 2.0 * log2(3.0) - log2(17.0)))
    end
end

@testset "spinflip, concurrence, concurrence_lb, formation, negativity" begin
    let rho = reshape(1.0:16.0, 4, 4),
        rho_f = [[16.0, -15.0, -14.0, 13.0] [-12.0, 11.0, 10.0, -9.0] [-8.0, 7.0, 6.0, -5.0] [4.0, -3.0, -2.0, 1.0]]

        @test spinflip(rho) == rho_f
    end

    let rho = zeros(4, 4)

        rho[1, 1] = 0.5
        rho[4, 4] = 0.5

        let C = concurrence(rho)

            @test spinflip(rho) == rho
            @test C == 0.0
            @test concurrence_lb(rho) == 0.0
            @test formation(C) == 0.0
            @test negativity(rho) == 0.0
        end
    end

    let rho = zeros(4, 4)

        rho[1, 1] = 0.125
        for i in 2:3, j in 2:3
            rho[i, j] = 0.375
        end
        rho[4, 4] = 0.125

        let C = concurrence(rho),
            a = (2 + sqrt(3.0)) / 4.0,
            b = (2 - sqrt(3.0)) / 4.0

            @test spinflip(rho) == rho
            @test C == 0.5
            @test concurrence_lb(rho) == sqrt(3.0)/4.0
            @test formation(C) == -a * log2(a) - b * log2(b)
            @test negativity(rho) == log2(1.5)
        end
    end

    let rho = zeros(4, 4)

        for i in 2:3, j in 2:3
            rho[i, j] = 0.5
        end

        let C = concurrence(rho)

            @test spinflip(rho) == rho
            @test C == 1.0
            @test concurrence_lb(rho) == 1.0
            @test formation(C) == 1.0
            @test negativity(rho) == 1.0
        end
    end

    let rho = Matrix{ComplexF64}(undef, 4, 4)

        for i in 1:4
            for j in 1:4
                if i < j
                    rho[i, j] = 0.0625im
                elseif i == j
                    rho[i, j] = 0.25
                else
                    rho[i, j] = -0.0625im
                end
            end
        end

        let C = concurrence(rho),
            rho_f = copy(transpose(rho))

            rho_f[1, 4] *= -1
            rho_f[2, 3] *= -1
            rho_f[3, 2] *= -1
            rho_f[4, 1] *= -1

            @test spinflip(rho) == rho_f
            @test C == 0.0
            @test concurrence_lb(rho) == 0.0
            @test formation(C) == 0.0
            @test isapprox(negativity(rho), 0.0, atol=1e-15)
        end
    end
end

end


text:
############################################################
# FactCheck.jl
# A testing framework for Julia
# http://github.com/JuliaLang/FactCheck.jl
# MIT Licensed
############################################################
# Assertion helpers
# - not
# - anything
# - truthy, falsey/falsy
# - exactly
# - roughly
############################################################

# not: logical not for values and functions
not(x) = isa(x, Function) ? (y) -> !x(y) :
                            (y) -> x != y

# anything: anything but nothing
anything(x) = (x != nothing)

# truthy: not `nothing`, false (== 0)
# falsy/falsey: not truthy
truthy(x) = (x != nothing) && (x != false)
falsey(x) = not(truthy(x))
falsy = falsey

# exactly: tests object/function equality (i.e. ===)
exactly(x) = (y) -> is(x, y)

# approx/roughly: Comparing numbers approximately
roughly(x::Number; kvtols...) = (y::Number) -> isapprox(y, x; kvtols...)

roughly(A::AbstractArray; kvtols...) = (B::AbstractArray) -> begin
    size(A) != size(B) && return false
    for i in 1:length(A)
        if !isapprox(A[i], B[i]; kvtols...)
            return false
        end
    end
    return true
end


text:
using AbstractPlotting
using AbstractPlotting.MakieLayout
using GLMakie

using Roots

## Assumes f(a), f(b) are zero
## only 1 or 2 solutions for f(x) = f(c) for any c in [a,b]
f(x) = 1 - x^4
a = -1
b = 1

    descr = """
Adjust the point (c, f(c)) with c > 0 to find the inscribed rectangle with maximal area
"""

function _inscribed_area(c)
    zs = fzeros(u -> f(u) - f(c), a, b)
    length(zs) <= 1 ? 0 : abs(zs[1] - zs[2]) * f(c)
end
D(f) = c -> (f(c + 1e-4) - f(c))/1e-4
function answer()
    h = 1e-4
    zs = fzeros(D(_inscribed_area), 0, b-h)
    a,i = findmax(_inscribed_area.(zs))
    a
end

Answer = answer()

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

                if idx[] == 3
                    pos = to_world(scene, Point2f0(scene.events.mouseposition[]))

                    # we work with components not vector
                    x,y = pos
                    c = clamp(x, a, b)
                    zs = fzeros(u -> f(u) - f(c), a , b)

                    if length(zs) == 1
                        c′ = c = first(zs)
                    else
                        c′, c = zs
                    end

                    points[][1] = [c′, 0]
                    points[][2] = [c′, f(c′)]
                    points[][3] = [c, f(c)]
                    points[][4] = [c, 0]
                    points[] = points[]
                end
            end
        else
            dragstart[] = false
        end
        return
    end
end

c = b/2
c′ = first(fzeros(u -> f(u) - f(c), a , b))
area = round(abs(c - c′) * f(c), digits=4)

points = Node(Point2f0[(c′,0), (c′, f(c′)), (c, f(c)), (c, 0)])

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

p = layout[0,0] = LScene(scene)

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

label = layout[end, 1] = LText(scene, "Area = $area")
layout[end+1,1] = LText(scene, chomp(descr))

polycolor = lift(points) do pts
    c′ = pts[1][1]
    c = pts[3][1]
    area = round(abs(c - c′)*f(c), digits=4)

    lbl = "Area = $area"
    label.text[] = lbl

    if abs(area - Answer) <= 1e-4
        :green
    else
        :gray75
    end

end

lines!(p.scene, a..b, f, strokecolor=:red, strokewidth=15)
lines!(p.scene, a..b, zero, strokecolor=:black, strokewidth=10)
poly!(p.scene, points, color = polycolor)
scatter!(p.scene, points, color = :white, strokewidth = 10, markersize = 0.05, strokecolor = :black, raw = true)
xlims!(p.scene, (a, b))

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

scene


text:
"""
    function pix_index_min_max(x::T, hsml::T, n_pixels::Int64) where T

Calculates the minimum and maximum pixel to which a particle contributes.
"""
function pix_index_min_max(x::T, hsml::T, n_pixels::Int64) where T

    iMin = max(floor(Integer, x - hsml), 0)
    iMax = min(floor(Integer, x + hsml), n_pixels-1)

    return iMin, iMax
end

"""
    function get_dxyz(x::Real, hsml::Real, i::Integer)

Calculates the extent of the current particle size in units of pixels.
"""
function get_dxyz(x::Real, hsml::Real, i::Integer)
    min(x + hsml, i + 1) - max(x - hsml, i)
end

"""
    function get_x_dx(x, hsml, i)

"""
@inline function get_x_dx(x, hsml, i)

    dx = get_dxyz(x, hsml, i)
    x_dist = x - i - 0.5

    return x_dist, dx
end

"""
    function get_xyz( pos, hsml, k::Integer,
                      par::mappingParameters)

Calculates `x, y, z` position in units of pixels and performs periodic mapping, if required.
"""
function get_xyz( pos, hsml, k::Integer,
                  par::mappingParameters)

    if par.periodic

        # perform periodic mapping of particle
        x, y, z = find_position_periodic(pos, k, par.boxsize)

        # check if the particle is still contained in the image frame
        # after the periodic mapping
        if !particle_in_image(x, y, z, hsml, par.halfsize)

            # if not return skip_particle as true
            return x, y, z, true
        end

    else
        # unpack array
        x, y, z = pos
    end

    # convert from code units to pixels
    x *= par.len2pix
    y *= par.len2pix

    x += 0.5 * par.Npixels[1]
    y += 0.5 * par.Npixels[2]

    return x, y, z, false

end

"""
    function update_image( image::Real, w_image::Real,
                           wk::Real, bin_q::Real,
                           geometry_norm::Real )

Applies the different contributions to the image and the weight image.
"""
function update_image( image::Float64, w_image::Float64,
                        wk::Float64, bin_q::Float64,
                        geometry_norm::Float64 )

    if wk > 0.0
        pix_weight  = geometry_norm * wk
        image      += bin_q * pix_weight
        w_image    += pix_weight
    end

    return image, w_image
end


text:
#############################
# Internal Structures
#############################

struct LoggedStaticPool <: AbstractThreadPool
    tids :: Vector{Int}
    lck  :: ReentrantLock
    recs :: Vector{ThreadRecord}
    t0   :: UInt64
    log  :: ThreadLog
end

#############################
# Constructors / Finalizer
#############################

"""
    LoggedStaticPool(init_thrd=1, nthrds=Threads.nthreads())

The main LoggedStaticPool object.
"""
function LoggedStaticPool(init_thrd::Integer=1, nthrds::Integer=Threads.nthreads())
    thrd0 = min(init_thrd, Threads.nthreads())
    thrd1 = min(thrd0+nthrds-1, Threads.nthreads())
    return LoggedStaticPool(thrd0:thrd1, ReentrantLock(), ThreadRecord[], time_ns(), ThreadLog())
end

function Base.close(pool::LoggedStaticPool)
    _recordstolog!(pool.log, pool.recs, pool. t0)
end

#############################
# ThreadPool Interface
#############################

function tmap(pool::LoggedStaticPool, fn::Function, itr)::Vector{_detect_type(fn, itr)}
    N = length(itr)
    sizehint!(pool.recs, N)
    result = Vector{_detect_type(fn, itr)}(undef, N)
    nts = length(pool.tids)
    n = div(N,nts)
    r = N % nts

    _fn = (tind) -> begin
        n0 = (tind-1)*n + 1 + (nts-tind+1 > r ? 0 : tind-nts+1)
        n1 = n0-1 + n + (nts-tind+1 <= r ? 1 : 0)
        tid = Threads.threadid()
        for i in n0:n1
            lock(pool.lck)
            push!(pool.recs, (i, tid, true, time_ns()))
            unlock(pool.lck)
            @inbounds result[i] = fn(Base.unsafe_getindex(itr, i))
            lock(pool.lck)
            push!(pool.recs, (i, tid, false, time_ns()))
            unlock(pool.lck)
        end
    end

    Threads.@threads for tid in 1:Threads.nthreads()
        ind = findfirst(t->t==tid, pool.tids)
        isnothing(ind) || _fn(ind)
    end

    return result
end


text:
# TODO: include matype options in macd, rsi, and adx functions

@doc doc"""
momentum{Float64}(x::Vector{Float64}, n::Int64=1)

Momentum indicator (price now vs price `n` periods back)
""" ->
function momentum{Float64}(x::Vector{Float64}, n::Int64=1)
    @assert n>0 "Argument n must be positive."
    return diffn(x, n)
end

@doc """
roc{Float64}(x::Vector{Float64}, n::Int64=1)

Rate of change indicator (percent change between i'th observation and (i-n)'th observation)
""" ->
function roc{Float64}(x::Vector{Float64}, n::Int64=1)
    @assert n<size(x,1) && n>0 "Argument n out of bounds."
    out = zeros(x)
    @inbounds for i = n:size(x,1)
        out[i] = x[i]/x[i-n] - 1.0
    end
    out[1:n] = NaN
    return out * 100.0
end

@doc doc"""
macd{Float64}(x::Vector{Float64}, nfast::Int64=12, nslow::Int64=26, nsig::Int64=9)

Moving average convergence-divergence

`Output:`

- Column 1: MACD
- Column 2: MACD Signal Line
- Column 3: MACD Histogram
""" ->
function macd{Float64}(x::Vector{Float64}, nfast::Int64=12, nslow::Int64=26, nsig::Int64=9;
                       fastMA::Function=ema, slowMA::Function=ema, signalMA::Function=sma)
    out = zeros(Float64, (length(x),3))
    out[:,1] = fastMA(x, nfast) - slowMA(x, nslow)
    out[:,2] = signalMA(out[:,1], nsig)
    out[:,3] = out[:,1] - out[:,2]
    return out
end

@doc doc"""
rsi{Float64}(x::Vector{Float64}, n::Int64=14; ma::Function=ema, args...)

Relative strength index
""" ->
function rsi{Float64}(x::Vector{Float64}, n::Int64=14; ma::Function=ema, args...)
    @assert n<size(x,1) && n>0 "Argument n is out of bounds."
    N = size(x,1)
    ups = zeros(N)
    dns = zeros(N)
    zro = 0.0
    dx = [NaN; diff(x)]
    @inbounds for i=2:N
        if dx[i] > zro
            ups[i] = dx[i]
        elseif dx[i] < zro
            dns[i] = -dx[i]
        end
    end
    rs = [NaN; ma(ups[2:end], n; args...) ./ ma(dns[2:end], n; args...)]
    return 100.0 - 100.0 ./ (1.0 + rs)
end

@doc doc"""
adx{Float64}(hlc::Array{Float64}, n::Int64=14; wilder=true)

Average directional index

`Output:`

- Column 1: DI+
- Column 2: DI-
- Column 3: ADX
""" ->
function adx{Float64}(hlc::Array{Float64}, n::Int64=14; ma::Function=ema, args...)
    @assert n<size(hlc,1) && n>0 "Argument n is out of bounds."
    if size(hlc,2) != 3
        error("HLC array must have three columns")
    end
    N = size(hlc,1)
    updm = zeros(N)
    dndm = zeros(N)
    updm[1] = dndm[1] = NaN
    @inbounds for i = 2:N
        upmove = hlc[i,1] - hlc[i-1,1]
        dnmove = hlc[i-1,2] - hlc[i,2]
        if upmove > dnmove && upmove > 0.0
            updm[i] = upmove
        elseif dnmove > upmove && dnmove > 0.0
            dndm[i] = dnmove
        end
    end
    dip = [NaN; ma(updm[2:N], n; args...)] ./ atr(hlc, n) * 100.0
    dim = [NaN; ma(dndm[2:N], n; args...)] ./ atr(hlc, n) * 100.0
    dmx = abs(dip-dim) ./ (dip+dim)
    adx = [fill(NaN,n); ma(dmx[n+1:N], n; args...)] * 100.0
    return [dip dim adx]
end

@doc doc"""
psar{Float64}(hl::Array{Float64,2}, af::Float64=0.02, af_max::Float64=0.2, af_min::Float64=0.02)

Parabolic stop and reverse (SAR)

- hl: 2D array of high and low prices in first and second columns respectively
- af_mi: starting/initial value for acceleration factor
- af_max: maximum acceleration factor (accel factor capped at this value)
- af_inc: increment to the acceleration factor (speed of increase in accel factor)
""" ->
function psar{Float64}(hl::Array{Float64}, af_min::Float64=0.02, af_max::Float64=0.2, af_inc::Float64=af_min)
    @assert af_min<1.0 && af_min>0.0 "Argument af_min must be in [0,1]."
    @assert af_max<1.0 && af_max>0.0 "Argument af_max must be in [0,1]."
    @assert af_inc<1.0 && af_inc>0.0 "Argument af_inc must be in [0,1]."
    @assert size(hl,2) == 2 "Argument hl must have 2 columns."
    ls0 = 1
    ls = 0
    af0 = af_min
    af = 0.0
    ep0 = hl[1,1]
    ep = 0.0
    maxi = 0.0
    mini = 0.0
    sar = zeros(Float64,size(hl,1))
    sar[1] = hl[1,2] - std(hl[:,1]-hl[:,2])
    @inbounds for i = 2:size(hl,1)
        ls = ls0
        ep = ep0
        af = af0
        mini = min(hl[i-1,2], hl[i,2])
        maxi = max(hl[i-1,1], hl[i,1])
        # Long/short signals and local extrema
        if (ls == 1)
            ls0 = hl[i,2] > sar[i-1] ? 1 : -1
            ep0 = max(maxi, ep)
        else
            ls0 = hl[i,1] < sar[i-1] ? -1 : 1
            ep0 = min(mini, ep)
        end
        # Acceleration vector
        if ls0 == ls  # no signal change
            sar[i] = sar[i-1] + af*(ep-sar[i-1])
            af0 = (af == af_max) ? af_max : (af + af_inc)
            if ls0 == 1  # current long signal
                af0 = (ep0 > ep) ? af0 : af
                sar[i] = min(sar[i], mini)
            else  # current short signal
                af0 = (ep0 < ep) ? af0 : af
                sar[i] = max(sar[i], maxi)
            end
        else  # new signal
            af0 = af_min
            sar[i] = ep0
        end
    end
    return sar
end

@doc doc"""
kst{Float64}(x::Vector{Float64},
                      nroc::Vector{Int64}=[10,15,20,30], navg::Vector{Int64}=[10,10,10,15],
                      wgts::Vector{Int64}=collect(1:length(nroc)); ma::Function=sma)

KST (Know Sure Thing) -- smoothed and summed rates of change
""" ->
function kst{Float64}(x::Vector{Float64}, nroc::Vector{Int64}=[10,15,20,30], navg::Vector{Int64}=[10,10,10,15];
                      wgts::Vector{Int64}=collect(1:length(nroc)), ma::Function=sma)
    @assert length(nroc) == length(navg)
    @assert all(nroc.>0) && all(nroc.<size(x,1))
    @assert all(navg.>0) && all(navg.<size(x,1))
    N = length(x)
    k = length(nroc)
    out = zeros(x)
    @inbounds for j = 1:k
        out += ma(roc(x, nroc[j]), navg[j]) * wgts[j]
    end
    return out
end

@doc doc"""
wpr{Float64}(hlc::Array{Float64,2}, n::Int64=14)

Williams %R
""" ->
function wpr{Float64}(hlc::Array{Float64,2}, n::Int64=14)
    hihi = runmax(hlc[:,1], n, false)
    lolo = runmin(hlc[:,2], n, false)
    return -100 * (hihi - hlc[:,3]) ./ (hihi - lolo)
end

@doc doc"""
cci{Float64}(hlc::Array{Float64,2}, n::Int64=20, c::Float64=0.015; ma::Function=sma)

Commodity channel index
""" ->
function cci{Float64}(hlc::Array{Float64,2}, n::Int64=20, c::Float64=0.015; ma::Function=sma, args...)
    tp = (hlc[:,1] + hlc[:,2] + hlc[:,3]) / 3.0
    dev = runmad(tp, n, false, fun=mean)
    avg = ma(tp, n; args...)
    return (tp - avg) ./ (c * dev)
end

@doc doc"""
stoch{Float64}(hlc::Array{Float64,2}, nK::Int64=14, nD::Int64=3;
               kind::ByteString="fast", ma::Function=sma, args...)

Stochastic oscillator (fast or slow)
""" ->
function stoch{Float64}(hlc::Array{Float64,2}, nK::Int64=14, nD::Int64=3;
                        kind::Symbol=:fast, ma::Function=sma, args...)
    @assert kind == :fast || kind == :slow "Argument `kind` must be either :fast or :slow"
    @assert nK<size(hlc,1) && nK>0 "Argument `nK` out of bounds."
    @assert nD<size(hlc,1) && nD>0 "Argument `nD` out of bounds."
    hihi = runmax(hlc[:,1], nK, false)
    lolo = runmin(hlc[:,2], nK, false)
    out = zeros(Float64, (size(hlc,1),2))
    out[:,1] = (hlc[:,3]-lolo) ./ (hihi-lolo) * 100.0
    out[:,2] = ma(out[:,1], nD; args...)
    if kind == "slow"
        out[:,1] = out[:,2]
        out[:,2] = ma(out[:,1], nD; args...)
    end
    return out
end

@doc doc"""
smi{Float64}(hlc::Array{Float64,2}, n::Int64=13, nFast::Int64=2, nSlow::Int64=25, nSig::Int64=9;
             maFast::Function=ema, maSlow::Function=ema, maSig::Function=sma)

SMI (stochastic momentum oscillator)
""" ->
function smi{Float64}(hlc::Array{Float64,2}, n::Int64=13, nFast::Int64=2, nSlow::Int64=25, nSig::Int64=9;
                      maFast::Function=ema, maSlow::Function=ema, maSig::Function=sma)
    hihi = runmax(hlc[:,1], n, false)
    lolo = runmin(hlc[:,2], n, false)
    hldif = hihi-lolo
    delta = hlc[:,3] - (hihi+lolo) / 2.0
    numer = maSlow(maFast(delta, nFast), nSlow)
    denom = maSlow(maFast(hldif, nFast), nSlow) / 2.0
    out = zeros(Float64, (size(hlc,1),2))
    out[:,1] = 100.0*(numer./denom)
    out[:,2] = maSig(out[:,1])
    return out
end


text:
"""
    bloch_redfield_tensor(H, a_ops; J=[], use_secular=true, secular_cutoff=0.1)

Create the super-operator for the Bloch-Redfield master equation such that ``\\dot ρ = R ρ`` based on the QuTiP implementation.

See QuTiP's documentation (http://qutip.org/docs/latest/guide/dynamics/dynamics-bloch-redfield.html) for more information and a brief derivation.

# Arguments
* `H`: Hamiltonian.
* `a_ops`: Nested list of [interaction operator, callback function] pairs for the Bloch-Redfield type processes where the callback function describes the environment spectrum for the corresponding interaction operator.
           The spectral functions must take the angular frequency as their only argument.
* `J=[]`: Vector containing the jump operators for the Linblad type processes (optional).
* `use_secular=true`: Specify whether or not to use the secular approximation.
* `secular_cutoff=0.1`: Cutoff to allow a degree of partial secularization. Terms are discarded if they are greater than (dw\\_min * secular cutoff) where dw\\_min is the smallest (non-zero) difference between any two eigenenergies of H.
                        This argument is only taken into account if use_secular=true.
"""
function bloch_redfield_tensor(H::AbstractOperator, a_ops::Array; J=[], use_secular=true, secular_cutoff=0.1)

    # use the eigenbasis
    H_evals, transf_mat = eigen(DenseOperator(H).data)
    H_ekets = [Ket(H.basis_l, transf_mat[:, i]) for i in 1:length(H_evals)]::Array{Ket{typeof(H.basis_l), Array{Complex{Float64}, 1}}, 1}

    #Define function for transforming to Hamiltonian eigenbasis
    function to_Heb(op)
        #Copy oper
        oper = copy(op)
        #Transform underlying array
        oper.data = inv(transf_mat) * oper.data * transf_mat
        return oper
    end

    N = length(H_evals)
    K = length(a_ops)

    #If only Lindblad collapse terms
    if K==0
        Heb = to_Heb(H)
        L = liouvillian(Heb, to_Heb.(J))
        L = sparse(L)
        return L, H_ekets
    end

    #Transform interaction operators to Hamiltonian eigenbasis
    A = Array{Complex{Float64}}(undef, N, N, K)
    for k in 1:K
        A[:, :, k] = to_Heb(a_ops[k][1]).data
    end

    # Trasition frequencies between eigenstates
    W = transpose(H_evals) .- H_evals

    #Array for spectral functions evaluated at transition frequencies
    Jw = Array{Complex}(undef, N, N, K)
    for k in 1:K
       # do explicit loops here
       for n in 1:N
           for m in 1:N
               Jw[m, n, k] = a_ops[k][2](W[n, m])
           end
       end
    end

    #Calculate secular cutoff scale
    W_flat = reshape(W, N*N)
    dw_min = minimum(abs.(W_flat[W_flat .!= 0.0]))

    #Pre-calculate mapping between global index I and system indices a,b
    Iabs = Array{Int}(undef, N*N, 3)
    indices = CartesianIndices((N,N))
    for I in 1:N*N
        Iabs[I, 1] = I
        Iabs[I, 2:3] = [indices[I].I...]
    end

    # Calculate Liouvillian for Lindblad temrs (unitary part + dissipation from J (if given)):
    Heb = to_Heb(H)
    L = liouvillian(Heb, to_Heb.(J))
    L = sparse(L)

    # Main Bloch-Redfield operators part
    rows = Int[]
    cols = Int[]
    data = Complex[]
    Is = view(Iabs, :, 1)
    As = view(Iabs, :, 2)
    Bs = view(Iabs, :, 3)

    for (I, a, b) in zip(Is, As, Bs)

        if use_secular
            Jcds = Int.(zeros(size(Iabs)))
            for (row, (I2, a2, b2)) in enumerate(zip(Is, As, Bs))
                if abs.(W[a, b] - W[a2, b2]) < dw_min * secular_cutoff
                    Jcds[row, :] = [I2 a2 b2]
                end
            end
            Jcds = transpose(Jcds)
            Jcds = Jcds[Jcds .!= 0]
            Jcds = reshape(Jcds, 3, Int(length(Jcds)/3))
            Jcds = transpose(Jcds)

            Js = view(Jcds, :, 1)
            Cs = view(Jcds, :, 2)
            Ds = view(Jcds, :, 3)

        else
            Js = Is
            Cs = As
            Ds = Bs
        end

        for (J, c, d) in zip(Js, Cs, Ds)

            elem = 0.5 * sum(view(A, c, a, :) .* view(A, b, d, :) .* (view(Jw, c, a, :) + view(Jw, d, b, :) ))

            if b == d
                elem -= 0.5 * sum(transpose(view(A, a, :, :)) .* transpose(view(A, c, :, :)) .* transpose(view(Jw, c, :, :) ))
            end

            if a == c
            elem -= 0.5 * sum(transpose(view(A, d, :, :)) .* transpose(view(A, b, :, :)) .* transpose(view(Jw, d, :, :) ))
            end

            #Add element
            if abs(elem) != 0.0
                push!(rows, I)
                push!(cols, J)
                push!(data, elem)
            end
        end
    end

    """
    Need to be careful here since sparse function is happy to create rectangular arrays but we dont want this behaviour so need to make square explicitly.
    """
    if !any(rows .== N*N) || !any(cols .== N*N) #Check if row or column list has (N*N)th element, if not then add one
        push!(rows, N*N)
        push!(cols, N*N)
        push!(data, 0.0)
    end

    R = sparse(rows, cols, data) #Careful with rows/cols ordering...

    #Add Bloch-Redfield part to Linblad Liouvillian calculated earlier
    L.data = L.data + R

    return L, H_ekets

end #Function

"""
    timeevolution.master_bloch_redfield(tspan, rho0, R, H; <keyword arguments>)

Time-evolution according to a Bloch-Redfield master equation.

# Arguments
* `tspan`: Vector specifying the points of time for which output should
        be displayed.
* `rho0`: Initial density operator. Can also be a state vector which is
        automatically converted into a density operator.
* `H`: Arbitrary operator specifying the Hamiltonian.
* `R`: Bloch-Redfield tensor describing the time-evolution ``\\dot ρ = R ρ`` (see timeevolution.bloch\\_redfield\\_tensor).
* `fout=nothing`: If given, this function `fout(t, rho)` is called every time
        an output should be displayed. ATTENTION: The given state rho is not
        permanent! It is still in use by the ode solver and therefore must not
        be changed.
* `kwargs...`: Further arguments are passed on to the ode solver.
"""
function master_bloch_redfield(tspan::Vector{Float64},
        rho0::T, L::SuperOperator{Tuple{B,B},Tuple{B,B}},
        H::AbstractOperator{B,B}; fout::Union{Function,Nothing}=nothing,
        kwargs...) where {B<:Basis,T<:DenseOperator{B,B}}

    #Prep basis transf
    evals, transf_mat = eigen(dense(H).data)
    transf_op = DenseOperator(rho0.basis_l, transf_mat)
    inv_transf_op = DenseOperator(rho0.basis_l, inv(transf_mat))

    # rho as Ket and L as DataOperator
    basis_comp = rho0.basis_l^2
    rho0_eb = Ket(basis_comp, (inv_transf_op * rho0 * transf_op).data[:]) #Transform to H eb and convert to Ket
    L_ = isa(L, DenseSuperOperator) ? DenseOperator(basis_comp, L.data) : SparseOperator(basis_comp, L.data)

    # Derivative function
    dmaster_br_(t::Float64, rho::T2, drho::T2) where T2<:Ket = dmaster_br(drho, rho, L_)

    return integrate_br(tspan, dmaster_br_, rho0_eb, transf_op, inv_transf_op, fout; kwargs...)
end
master_bloch_redfield(tspan::Vector{Float64}, psi::Ket, args...; kwargs...) = master_bloch_redfield(tspan::Vector{Float64}, dm(psi), args...; kwargs...)

# Derivative ∂ₜρ = Lρ
function dmaster_br(drho::T, rho::T, L::DataOperator{B,B}) where {B<:Basis,T<:Ket{B}}
    QuantumOpticsBase.gemv!(1.0, L, rho, 0.0, drho)
end

# Integrate if there is no fout specified
function integrate_br(tspan::Vector{Float64}, dmaster_br::Function, rho::T,
                transf_op::T2, inv_transf_op::T2, ::Nothing;
                kwargs...) where {T<:Ket,T2<:DenseOperator}
    # Pre-allocate for in-place back-transformation from eigenbasis
    rho_out = copy(transf_op)
    tmp = copy(transf_op)
    tmp2 = copy(transf_op)

    # Define fout
    function fout(t::Float64, rho::T)
        tmp.data[:] = rho.data
        QuantumOpticsBase.gemm!(1.0, transf_op, tmp, 0.0, tmp2)
        QuantumOpticsBase.gemm!(1.0, tmp2, inv_transf_op, 0.0, rho_out)
        return copy(rho_out)
    end

    return integrate(tspan, dmaster_br, copy(rho.data), rho, copy(rho), fout; kwargs...)
end

# Integrate with given fout
function integrate_br(tspan::Vector{Float64}, dmaster_br::Function, rho::T,
                transf_op::T2, inv_transf_op::T2, fout::Function;
                kwargs...) where {T<:Ket,T2<:DenseOperator}
    # Pre-allocate for in-place back-transformation from eigenbasis
    rho_out = copy(transf_op)
    tmp = copy(transf_op)
    tmp2 = copy(transf_op)

    # Perform back-transfomration before calling fout
    function fout_(t::Float64, rho::T)
        tmp.data[:] = rho.data
        QuantumOpticsBase.gemm!(1.0, transf_op, tmp, 0.0, tmp2)
        QuantumOpticsBase.gemm!(1.0, tmp2, inv_transf_op, 0.0, rho_out)
        return fout(t, rho_out)
    end

    return integrate(tspan, dmaster_br, copy(rho.data), rho, copy(rho), fout_; kwargs...)
end


text:
"""
    PosixPath()
    PosixPath(str)

Represents any posix path (e.g., `/home/user/docs`)
"""
struct PosixPath <: AbstractPath
    segments::Tuple{Vararg{String}}
    root::String
end

PosixPath() = PosixPath(tuple(), "")

function PosixPath(str::AbstractString)
    str = string(str)
    root = ""

    if isempty(str)
        return PosixPath(tuple("."))
    end

    tokenized = split(str, POSIX_PATH_SEPARATOR)
    if isempty(tokenized[1])
        root = POSIX_PATH_SEPARATOR
    end
    return PosixPath(tuple(map(String, filter!(!isempty, tokenized))...), root)
end

ispathtype(::Type{PosixPath}, str::AbstractString) = Sys.isunix()
isabs(fp::PosixPath) = !isempty(fp.root)

function Base.expanduser(fp::PosixPath)
    p = fp.segments

    if p[1] == "~"
        if length(p) > 1
            return PosixPath(tuple(homedir(), p[2:end]...))
        else
            return PosixPath(tuple(homedir()))
        end
    end

    return fp
end


text:
using JSON, ZMQ
using Test
# using RoMEPlotting
using Caesar, Caesar.ZmqCaesar
using Distributions, RoME, IncrementalInference
using Unmarshal

# @testset "ZMQ Integration Test (Hexagonal Solver)" begin

addOdo2DJson = "{\n  \"covariance\": [\n    [\n      0.1,\n      0.0,\n      0.1\n    ],\n    [\n      0.1,\n      0.0,\n      0.1\n    ],\n    [\n      0.1,\n      0.0,\n      0.1\n    ]\n  ],\n  \"measurement\": [\n    10.0,\n    0.0,\n    1.0471975511965976\n  ],\n  \"robot_id\": \"Hexagonal\",\n  \"session_id\": \"cjz002\",\n  \"type\": \"addOdometry2D\"\n}";
addLandmark2DJson = "{\n  \"landmark_id\": \"l1\",\n  \"robot_id\": \"Hexagonal\",\n  \"session_id\": \"cjz002\",\n  \"type\": \"addLandmark2D\"\n}";
addOdo2DCmd = JSON.parse(addOdo2DJson)
addLandmark2DCmd = JSON.parse(addLandmark2DJson)

fg = initfg()
config = Dict{String, String}()

# Send a command locally.
function sendCmd(config, fg, cmd::Dict{String, Any})::String
    callback = getfield(Caesar.ZmqCaesar, Symbol(cmd["request"]))
    resp = callback(config, fg, cmd)
    return JSON.json(resp)
end

# Add an initial variable.
addVariableCmd = Dict{String, Any}("request" => "addVariable", "payload" => JSON.parse(JSON.json(VariableRequest("x0", "Pose2", 100, ["TEST"]))))
@test sendCmd(config, fg, addVariableCmd) == "{\"status\":\"OK\",\"id\":\"x0\"}"

# ls to see that it is added to the graph.
lsCmd = Dict{String, Any}("request" => "ls", "payload" => JSON.parse(JSON.json(lsRequest("true", "true"))))
@test sendCmd(config, fg, lsCmd) == "{\"variables\":[{\"factors\":[],\"id\":\"x0\"}]}"

# Add a prior factor
dist = Dict{String, Any}("distType" => "MvNormal", "mean" => Array{Float64}(undef,3), "cov" => [1.0,0,0,0,1.0,0,0,0,1.0])
factor = Dict{String, Any}("measurement" => [dist])
priorFactCmd = Dict{String, Any}("request" => "addFactor", "payload" => JSON.parse(JSON.json(FactorRequest(["x0"], "Prior", factor))))
@test sendCmd(config, fg, priorFactCmd) == "{\"status\":\"OK\",\"id\":\"x0f1\"}"

@test sendCmd(config, fg, lsCmd) == "{\"variables\":[{\"factors\":[\"x0f1\"],\"id\":\"x0\"}]}"

# Now let's add all 6 hexagonal nodes
odo = Dict{String, Any}(
    "distType" => "MvNormal",
    "mean" => [10.0,0,pi/3.0],
    "cov" => [0.1,0,0,0,0.1,0,0,0,0.1])
for i in 1:6
    # Add variable
    addVariableCmd = Dict{String, Any}("request" => "addVariable", "payload" => JSON.parse(JSON.json(VariableRequest("x$i", "Pose2", 100, ["TEST"]))))
    @test sendCmd(config, fg, addVariableCmd) == "{\"status\":\"OK\",\"id\":\"x$i\"}"

    # Now adding odo factor
    factor = Dict{String, Any}("measurement" => [odo])
    @show odoFactCmd = Dict{String, Any}("request" => "addFactor", "payload" => JSON.parse(JSON.json(FactorRequest(["x$(i-1)", "x$i"], "Pose2Pose2", factor))))
    @test sendCmd(config, fg, odoFactCmd) == "{\"status\":\"OK\",\"id\":\"x$(i-1)x$(i)f1\"}"
end

jResponse = JSON.parse(sendCmd(config, fg, lsCmd))
# TODO: Test the content of this response.
#@test sendCmd(config, fg, lsCmd) == "{\"variables\":[{\"factors\":[\"x0f1\",\"x0x1f1\"],\"id\":\"x0\"},{\"factors\":[\"x0x1f1\",\"x1x2f1\"],\"id\":\"x1\"},{\"factors\":[\"x1x2f1\",\"x2x3f1\"],\"id\":\"x2\"},{\"factors\":[\"x2x3f1\",\"x3x4f1\"],\"id\":\"x3\"},{\"factors\":[\"x3x4f1\",\"x4x5f1\"],\"id\":\"x4\"},{\"factors\":[\"x4x5f1\",\"x5x6f1\"],\"id\":\"x5\"},{\"factors\":[\"x5x6f1\"],\"id\":\"x6\"}]}"

# Set all variables solvable
setSolvableCmd = Dict{String, Any}("request" => "setSolvable", "payload" => JSON.parse(JSON.json(setSolvableRequest(nothing, 1))))
@test sendCmd(config, fg, setSolvableCmd) == "{\"status\":\"OK\"}"

# Call batch solve
solveTreeCmd = Dict{String, Any}("request" => "solveTree")
result = sendCmd(config, fg, solveTreeCmd)
@test JSON.parse(result)["status"] == "OK"

# Just for fun
getVariableCmd = Dict{String, Any}("request" => "getVariable", "payload" => "x0")
x0Ret = JSON.parse(sendCmd(config, fg, getVariableCmd))
getFactorCmd = Dict{String, Any}("request" => "getFactor", "payload" => "x0x1f1")
# TODO: Fix
# x0x1f1Ret = JSON.parse(sendCmd(config, fg, getFactorCmd))

# Test a factor that doesn't exist, should produce legible error
dist = Dict{String, Any}("distType" => "MvNormal", "mean" => Array{Float64}(undef,3), "cov" => [1.0,0,0,0,1.0,0,0,0,1.0])
factor = Dict{String, Any}("measurement" => [dist])
priorFactCmd = Dict{String, Any}("request" => "addFactor", "payload" => JSON.parse(JSON.json(FactorRequest(["x0"], "IDontExistAsAFactor", factor))))
@test_throws ErrorException sendCmd(config, fg, priorFactCmd)

# Get the KDEs
getMAPCmd = Dict{String, Any}("request" => "getVarMAPKDE", "payload" => "x0")
mapResult = sendCmd(config, fg, getMAPCmd)
getMAPMeanCmd = Dict{String, Any}("request" => "getVarMAPMean", "payload" => "x0")
mapResult = sendCmd(config, fg, getMAPMeanCmd)

#######
### STATE OF THE ART. BEYOND HERE THERE BE KRAKENS (AND EXCEPTIONS, YAAARGH!)
#######

# Get some plots
# drawPosesCmd = Dict{String, Any}("type" => "drawPoses", "plotParams" => Dict{String, Any}("widthPx" => 1024, "heightPx" => 768, "encoding" => "none"))
# plotResult = sendCmd(config, fg, drawPosesCmd)

# end


text:
using Genie.Router
using TasksController
using BooksController

route("/") do
  serve_static_file("/welcome.html")
end

route("/tasks", TasksController.index, named = :get_tasks)
route("/tasks/new", TasksController.new, named = :new_task)
route("/tasks/create", TasksController.create, method = GET, named = :create_task)
route("/tasks/update", TasksController.update, method = GET, named = :update_task)

max_stars_repo_name: yshutaro/GenieTodoList

max_issues_repo_name: yshutaro/GenieTodoList

max_forks_repo_name: yshutaro/GenieTodoList