text:
##########   I/O   ###########

import .LibMagick
import .LibOSXNative

abstract ImageFileType

# Database of image file extension, magic bytes, and file type stubs
# Add new items to this database using add_image_file_format (below)
#
# Filename extensions do not uniquely specify the image type, so we
# have extensions pointing to a list of candidate formats. On reading,
# this list is really just a hint to improve efficiency: in the end,
# it's the set of magic bytes in the file that determine the file
# format.
fileext = Dict{ByteString, Vector{Int}}()
filemagic = Array(Vector{Uint8}, 0)
filetype = Array(Any, 0)
filesrcloaded = Array(Bool, 0)
filesrc = Array(String, 0)

function _loadformat(index::Int)
    filename = joinpath("ioformats", filesrc[index])
    if !isfile(filename)
        filename = joinpath(Pkg.dir(), "Images", "src", "ioformats", filesrc[index])
    end
    include(filename)
    filesrcloaded[index] = true
end
function loadformat{FileType<:ImageFileType}(::Type{FileType})
    indx = find(filetype .== FileType)
    if length(indx) == 1
        _loadformat(indx[1])
    elseif isempty(indx)
        error("File format ", FileType, " not found")
    else
        error("File format ", FileType, " is entered multiple times")
    end
    nothing
end

function add_image_file_format{ImageType<:ImageFileType}(ext::ByteString, magic::Vector{Uint8}, ::Type{ImageType}, filecode::ASCIIString)
    # Check to see whether these magic bytes are already in the database
    for i in 1:length(filemagic)
        if magic == filemagic[i]
            fileext[ext] = [i]  # update/add extension lookup
            filetype[i] = ImageType
            return
        end
    end
    # Add a new entry on the list
    push!(filemagic, magic)
    len = length(filemagic)
    push!(filetype, ImageType)
    if haskey(fileext, ext)
        push!(fileext[ext], len)
    else
        fileext[ext] = [len]
    end
    push!(filesrcloaded, isempty(filecode))
    push!(filesrc, filecode)
end
add_image_file_format{ImageType<:ImageFileType}(ext::ByteString, magic::Vector{Uint8}, ::Type{ImageType}) = add_image_file_format(ext, magic, ImageType, "")
add_image_file_format{ImageType<:ImageFileType}(ext::ByteString, ::Type{ImageType}, filecode::ASCIIString) = add_image_file_format(ext, b"", ImageType, filecode)

# Define our fallback file formats now, because we need them in generic imread.
# This has no extension (and is not added to the database), because it is always used as a stream.
type ImageMagick <: ImageFileType end
type OSXNative <: ImageFileType end

function imread(filename::String)
    _, ext = splitext(filename)
    ext = lowercase(ext)

    img = open(filename, "r") do stream
        # Tries to read the image using a set of potential type candidates.
        # Returns the image if successful, `nothing` else.
        function tryread(candidates)
            if (index = image_decode_magic(stream, candidates)) > 0
                # Position to end of this type's magic bytes
                seek(stream, length(filemagic[index]))
                if !filesrcloaded[index]
                    _loadformat(index)
                end
                imread(stream, filetype[index])
            end
        end

        # Use the extension as a hint to determine file type
        if haskey(fileext, ext) && (img = tryread(fileext[ext])) != nothing
            return img
        end

        # Extension wasn't helpful, look at all known magic bytes
        if (img = tryread(1:length(filemagic))) != nothing
            return img
        end
    end

    if img != nothing
        return img
    end

    @osx_only begin
        img = imread(filename, OSXNative)
        if img != nothing
            return img
        end
    end

    # There are no registered readers for this type. Try using ImageMagick if available.
    if have_imagemagick
        return imread(filename, ImageMagick)
    else
        error("Do not know how to read file ", filename, ". Is ImageMagick installed properly? See README.")
    end
end

imread{C<:ColorValue}(filename::String, ::Type{C}) = imread(filename, ImageMagick, C)

# Identify via magic bytes
function image_decode_magic{S<:IO}(stream::S, candidates::AbstractVector{Int})
    maxlen = 0
    for i in candidates
        len = length(filemagic[i])
        maxlen = (len > maxlen) ? len : maxlen
    end
    # If there are no magic bytes, simply use the file extension
    if maxlen == 0 && length(candidates) == 1
        return candidates[1]
    end

    magicbuf = zeros(Uint8, maxlen)
    for i=1:maxlen
        if eof(stream) break end
        magicbuf[i] = read(stream, Uint8)
    end
    for i in candidates
        if length(filemagic[i]) == 0
            continue
        end
        ret = ccall(:memcmp, Int32, (Ptr{Uint8}, Ptr{Uint8}, Int), magicbuf, filemagic[i], min(length(filemagic[i]), length(magicbuf)))
        if ret == 0
            return i
        end
    end
    return -1
end

function imwrite(img, filename::String; kwargs...)
    _, ext = splitext(filename)
    ext = lowercase(ext)
    if haskey(fileext, ext)
        # Write using specific format
        candidates = fileext[ext]
        index = candidates[1]  # TODO?: use options, don't default to first
        if !filesrcloaded[index]
            _loadformat(index)
        end
        imwrite(img, filename, filetype[index]; kwargs...)
    elseif have_imagemagick
        # Fall back on ImageMagick
        imwrite(img, filename, ImageMagick; kwargs...)
    else
        error("Do not know how to write file ", filename)
    end
end

function imwrite{T<:ImageFileType}(img, filename::String, ::Type{T}; kwargs...)
    open(filename, "w") do s
        imwrite(img, s, T; kwargs...)
    end
end

# only mime writeable to PNG if 2D (used by IJulia for example)
import Base.mimewritable
Base.mimewritable(::MIME"image/png", img::AbstractImage) = sdims(img) == 2 && timedim(img) == 0
# We have to disable Color's display via SVG, because both will get sent with unfortunate results.
# See IJulia issue #229
mimewritable{T<:ColorValue}(::MIME"image/svg+xml", ::AbstractMatrix{T}) = false

function _writemime(stream::IO, ::MIME"image/png", img::AbstractImage; mapi=mapinfo(ImageMagick, img))
    assert2d(img)
    if isa(img, AbstractImageIndexed)
        # For now, convert to direct
        img = convert(Image, img)
    end
    A = data(img)
    nc = ncolorelem(img)
    npix = length(A)/nc
    while npix > 1e6
        A = restrict(A, coords_spatial(img))
        npix = length(A)/nc
    end
    if eltype(A) != eltype(img)
        mapi = similar(mapi, eltype(mapi), eltype(A))
    end
    wand = image2wand(share(img, A), mapi, nothing)
    blob = LibMagick.getblob(wand, "png")
    write(stream, blob)
end
writemime(stream::IO, mime::MIME"image/png", img::AbstractImage; mapi=mapinfo_writemime(img)) =
    _writemime(stream, mime, img, mapi=mapi)

mapinfo_writemime{T}(img::AbstractImage{Gray{T}}) = mapinfo(Gray{Ufixed8},img)
mapinfo_writemime{C<:ColorValue}(img::AbstractImage{C}) = mapinfo(RGB{Ufixed8},img)
mapinfo_writemime{T}(img::AbstractImage{GrayAlpha{T}}) = mapinfo(GrayAlpha{Ufixed8},img)
mapinfo_writemime{AC<:AbstractAlphaColorValue}(img::AbstractImage{AC}) = mapinfo(RGBA{Ufixed8},img)
mapinfo_writemime(img::AbstractImage) = mapinfo(Ufixed8,img)

#### Implementation of specific formats ####

#### OSX Native readers from CoreGraphics

imread(filename::String, ::Type{OSXNative}) = LibOSXNative.imread(filename)

#### ImageMagick library

# fixed type for depths > 8
const ufixedtype = @Dict(10=>Ufixed10, 12=>Ufixed12, 14=>Ufixed14, 16=>Ufixed16)

function imread(filename::String, ::Type{ImageMagick})
    wand = LibMagick.MagickWand()
    LibMagick.readimage(wand, filename)
    LibMagick.resetiterator(wand)
    imtype = LibMagick.getimagetype(wand)
    # Determine what we need to know about the image format
    sz = size(wand)
    n = LibMagick.getnumberimages(wand)
    if n > 1
        sz = tuple(sz..., n)
    end
    havealpha = LibMagick.getimagealphachannel(wand)
    prop = @Dict("spatialorder" => ["x", "y"], "pixelspacing" => [1,1])
    cs = LibMagick.getimagecolorspace(wand)
    if imtype == "GrayscaleType" || imtype == "GrayscaleMatteType"
        cs = "Gray"
    end
    prop["IMcs"] = cs

    depth = LibMagick.getimagechanneldepth(wand, LibMagick.DefaultChannels)
    if depth <= 8
        T = Ufixed8     # always use 8-bit for 8-bit and less
    else
        T = ufixedtype[depth]
    end

    channelorder = cs
    if havealpha
        if channelorder == "sRGB" || channelorder == "RGB"
            if is_little_endian
                T, channelorder = BGRA{T}, "BGRA"
            else
                T, channelorder = ARGB{T}, "ARGB"
            end
        elseif channelorder == "Gray"
            T, channelorder = GrayAlpha{T}, "IA"
        else
            error("Cannot parse colorspace $channelorder")
        end
    else
        if channelorder == "sRGB" || channelorder == "RGB"
            T, channelorder = RGB{T}, "RGB"
        elseif channelorder == "Gray"
            T, channelorder = Gray{T}, "I"
        else
            error("Cannot parse colorspace $channelorder")
        end
    end
    # Allocate the buffer and get the pixel data
    buf = Array(T, sz...)
    LibMagick.exportimagepixels!(buf, wand, cs, channelorder)
    if n > 1
        prop["timedim"] = ndims(buf)
    end
    Image(buf, prop)
end

imread{C<:ColorType}(filename::String, ::Type{ImageMagick}, ::Type{C}) = convert(Image{C}, imread(filename, ImageMagick))

function imwrite(img, filename::String, ::Type{ImageMagick}; mapi = mapinfo(ImageMagick, img), quality = nothing)
    wand = image2wand(img, mapi, quality)
    LibMagick.writeimage(wand, filename)
end

function image2wand(img, mapi, quality)
    if isa(img, AbstractImageIndexed)
        # For now, convert to direct
        img = convert(Image, img)
    end
    imgw = map(mapi, img)
    imgw = permutedims_horizontal(imgw)
    have_color = colordim(imgw)!=0
    if ndims(imgw) > 3+have_color
        error("At most 3 dimensions are supported")
    end
    wand = LibMagick.MagickWand()
    if haskey(img, "IMcs")
        cs = img["IMcs"]
    else
        cs = colorspace(imgw)
        if in(cs, ("RGB", "RGBA", "ARGB", "BGRA"))
            cs = LibMagick.libversion > v"6.7.5" ? "sRGB" : "RGB"
        end
    end
    channelorder = colorspace(imgw)
    if channelorder == "Gray"
        channelorder = "I"
    elseif channelorder == "GrayAlpha"
        channelorder = "IA"
    end
    tmp = to_explicit(to_contiguous(data(imgw)))
    LibMagick.constituteimage(tmp, wand, cs, channelorder)
    if quality != nothing
        LibMagick.setimagecompressionquality(wand, quality)
    end
    LibMagick.resetiterator(wand)
    wand
end

# ImageMagick mapinfo client. Converts to RGB and uses Ufixed.
mapinfo{T<:Ufixed}(::Type{ImageMagick}, img::AbstractArray{T}) = MapNone{T}()
mapinfo{T<:FloatingPoint}(::Type{ImageMagick}, img::AbstractArray{T}) = ClampMinMax(Ufixed8, zero(T), one(T))
for ACV in (ColorValue, AbstractRGB,AbstractGray)
    for CV in subtypes(ACV)
        (length(CV.parameters) == 1 && !(CV.abstract)) || continue
        CVnew = CV<:AbstractGray ? Gray : RGB
        @eval mapinfo{T<:Ufixed}(::Type{ImageMagick}, img::AbstractArray{$CV{T}}) = MapNone{$CVnew{T}}()
        @eval mapinfo{T<:FloatingPoint}(::Type{ImageMagick}, img::AbstractArray{$CV{T}}) =
            Clamp{$CVnew{Ufixed8}}()
        CVnew = CV<:AbstractGray ? Gray : BGR
        for AC in subtypes(AbstractAlphaColorValue)
            (length(AC.parameters) == 2 && !(AC.abstract)) || continue
            @eval mapinfo{T<:Ufixed}(::Type{ImageMagick}, img::AbstractArray{$AC{$CV{T},T}}) = MapNone{$AC{$CVnew{T},T}}()
            @eval mapinfo{T<:FloatingPoint}(::Type{ImageMagick}, img::AbstractArray{$AC{$CV{T},T}}) = Clamp{$AC{$CVnew{Ufixed8}, Ufixed8}}()
        end
    end
end

# Make the data contiguous in memory, this is necessary for imagemagick since it doesn't handle stride.
to_contiguous(A::AbstractArray) = A
to_contiguous(A::SubArray) = copy(A)

to_explicit(A::AbstractArray) = A
to_explicit{T<:Ufixed}(A::AbstractArray{T}) = reinterpret(FixedPointNumbers.rawtype(T), A)
to_explicit{T<:Ufixed}(A::AbstractArray{RGB{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(3, size(A)...))
to_explicit{T<:FloatingPoint}(A::AbstractArray{RGB{T}}) = to_explicit(map(ClipMinMax(RGB{Ufixed8}, zero(RGB{T}), one(RGB{T})), A))
to_explicit{T<:Ufixed}(A::AbstractArray{Gray{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, size(A))
to_explicit{T<:FloatingPoint}(A::AbstractArray{Gray{T}}) = to_explicit(map(ClipMinMax(Gray{Ufixed8}, zero(Gray{T}), one(Gray{T})), A))
to_explicit{T<:Ufixed}(A::AbstractArray{GrayAlpha{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(2, size(A)...))
to_explicit{T<:FloatingPoint}(A::AbstractArray{GrayAlpha{T}}) = to_explicit(map(ClipMinMax(GrayAlpha{Ufixed8}, zero(GrayAlpha{T}), one(GrayAlpha{T})), A))
to_explicit{T<:Ufixed}(A::AbstractArray{BGRA{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(4, size(A)...))
to_explicit{T<:FloatingPoint}(A::AbstractArray{BGRA{T}}) = to_explicit(map(ClipMinMax(BGRA{Ufixed8}, zero(BGRA{T}), one(BGRA{T})), A))
to_explicit{T<:Ufixed}(A::AbstractArray{RGBA{T}}) = reinterpret(FixedPointNumbers.rawtype(T), A, tuple(4, size(A)...))
to_explicit{T<:FloatingPoint}(A::AbstractArray{RGBA{T}}) = to_explicit(map(ClipMinMax(RGBA{Ufixed8}, zero(RGBA{T}), one(RGBA{T})), A))

# Write values in permuted order
let method_cache = Dict()
global writepermuted
function writepermuted(stream, img, mapi::MapInfo, perm; gray2color::Bool = false)
    cd = colordim(img)
    key = (perm, cd, gray2color)
    if !haskey(method_cache, key)
        mapfunc = cd > 0 ? (:map1) : (:map)
        loopsyms = [symbol(string("i_",d)) for d = 1:ndims(img)]
        body = gray2color ? quote
                g = $mapfunc(mapi, img[$(loopsyms...)])
                write(stream, g)
                write(stream, g)
                write(stream, g)
            end : quote
                write(stream, $mapfunc(mapi, img[$(loopsyms...)]))
            end
        loopargs = [:($(loopsyms[d]) = 1:size(img, $d)) for d = 1:ndims(img)]
        loopexpr = Expr(:for, Expr(:block, loopargs[perm[end:-1:1]]...), body)
        f = @eval begin
            local _writefunc_
            function _writefunc_(stream, img, mapi)
                $loopexpr
            end
        end
    else
        f = method_cache[key]
    end
    f(stream, img, mapi)
    nothing
end
end

function write{T<:Ufixed}(io::IO, c::AbstractRGB{T})
    write(io, reinterpret(c.r))
    write(io, reinterpret(c.g))
    write(io, reinterpret(c.b))
end

function write(io::IO, c::RGB24)
    write(io, red(c))
    write(io, green(c))
    write(io, blue(c))
end

## PPM, PGM, and PBM ##
type PPMBinary <: ImageFileType end
type PGMBinary <: ImageFileType end
type PBMBinary <: ImageFileType end

add_image_file_format(".ppm", b"P6", PPMBinary)
#add_image_file_format(".ppm", b"P3", PPMASCII)
add_image_file_format(".pgm", b"P5", PGMBinary)
#add_image_file_format(".pgm", b"P2", PGMASCII)
add_image_file_format(".pbm", b"P4", PBMBinary)
#add_image_file_format(".pbm", b"P1", PBMASCII)

function parse_netpbm_size(stream::IO)
    szline = strip(readline(stream))
    while isempty(szline) || szline[1] == '#'
        szline = strip(readline(stream))
    end
    parseints(szline, 2)
end

function parse_netpbm_maxval(stream::IO)
    skipchars(stream, isspace, linecomment='#')
    maxvalline = strip(readline(stream))
    parseint(maxvalline)
end

function imread{S<:IO}(stream::S, ::Type{PPMBinary})
    w, h = parse_netpbm_size(stream)
    maxval = parse_netpbm_maxval(stream)
    local dat
    if maxval <= 255
        datraw = read(stream, Ufixed8, 3, w, h)
        dat = reinterpret(RGB{Ufixed8}, datraw, (w, h))
    elseif maxval <= typemax(Uint16)
        # read first as Uint16 so the loop is type-stable, then convert to Ufixed
        datraw = Array(Uint16, 3, w, h)
        # there is no endian standard, but netpbm is big-endian
        if !is_little_endian
            for indx = 1:3*w*h
                datraw[indx] = read(stream, Uint16)
            end
        else
            for indx = 1:3*w*h
                datraw[indx] = bswap(read(stream, Uint16))
            end
        end
        # Determine the appropriate Ufixed type
        T = ufixedtype[iceil(log2(maxval)/2)<<1]
        dat = reinterpret(RGB{T}, datraw, (w, h))
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
    T = eltype(dat)
    Image(dat, @Dict("spatialorder" => ["x", "y"], "pixelspacing" => [1,1]))
end

function imread{S<:IO}(stream::S, ::Type{PGMBinary})
    w, h = parse_netpbm_size(stream)
    maxval = parse_netpbm_maxval(stream)
    local dat
    if maxval <= 255
        dat = read(stream, Ufixed8, w, h)
    elseif maxval <= typemax(Uint16)
        datraw = Array(Uint16, w, h)
        if !is_little_endian
            for indx = 1:w*h
                datraw[indx] = read(stream, Uint16)
            end
        else
            for indx = 1:w*h
                datraw[indx] = bswap(read(stream, Uint16))
            end
        end
        # Determine the appropriate Ufixed type
        T = ufixedtype[iceil(log2(maxval)/2)<<1]
        dat = reinterpret(RGB{T}, datraw, (w, h))
    else
        error("Image file may be corrupt. Are there really more than 16 bits in this image?")
    end
    T = eltype(dat)
    Image(dat, @Dict("colorspace" => "Gray", "spatialorder" => ["x", "y"], "pixelspacing" => [1,1]))
end

function imread{S<:IO}(stream::S, ::Type{PBMBinary})
    w, h = parse_netpbm_size(stream)
    dat = BitArray(w, h)
    nbytes_per_row = iceil(w/8)
    for irow = 1:h, j = 1:nbytes_per_row
        tmp = read(stream, Uint8)
        offset = (j-1)*8
        for k = 1:min(8, w-offset)
            dat[offset+k, irow] = (tmp>>>(8-k))&0x01
        end
    end
    Image(dat, @Dict("spatialorder" => ["x", "y"], "pixelspacing" => [1,1]))
end

function imwrite(img, filename::String, ::Type{PPMBinary})
    open(filename, "w") do stream
        write(stream, "P6\n")
        write(stream, "# ppm file written by Julia\n")
        imwrite(img, stream, PPMBinary)
    end
end

pnmmax{T<:FloatingPoint}(::Type{T}) = 255
pnmmax{T<:Ufixed}(::Type{T}) = reinterpret(FixedPointNumbers.rawtype(T), one(T))
pnmmax{T<:Unsigned}(::Type{T}) = typemax(T)

function imwrite{T<:ColorValue}(img::AbstractArray{T}, s::IO, ::Type{PPMBinary}, mapi = mapinfo(ImageMagick, img))
    w, h = widthheight(img)
    TE = eltype(T)
    mx = pnmmax(TE)
    write(s, "$w $h\n$mx\n")
    p = permutation_horizontal(img)
    writepermuted(s, img, mapi, p; gray2color = T <: AbstractGray)
end

function imwrite{T}(img::AbstractArray{T}, s::IO, ::Type{PPMBinary}, mapi = mapinfo(ImageMagick, img))
    w, h = widthheight(img)
    cs = colorspace(img)
    in(cs, ("RGB", "Gray")) || error("colorspace $cs not supported")
    mx = pnmmax(T)
    write(s, "$w $h\n$mx\n")
    p = permutation_horizontal(img)
    writepermuted(s, img, mapi, p; gray2color = cs == "Gray")
end

function parseints(line, n)
    ret = Array(Int, n)
    pos = 1
    for i = 1:n
        pos2 = search(line, ' ', pos)
        if pos2 == 0
            pos2 = length(line)+1
        end
        ret[i] = parseint(line[pos:pos2-1])
        pos = pos2+1
        if pos > length(line) && i < n
            error("Line terminated without finding all ", n, " integers")
        end
    end
    tuple(ret...)

end

# Permute to a color, horizontal, vertical, ... storage order (with time always last)
function permutation_horizontal(img)
    cd = colordim(img)
    td = timedim(img)
    p = spatialpermutation(["x", "y"], img)
    if cd != 0
        p[p .>= cd] += 1
        insert!(p, 1, cd)
    end
    if td != 0
        push!(p, td)
    end
    p
end

permutedims_horizontal(img) = permutedims(img, permutation_horizontal(img))

### Register formats for later loading here
type Dummy <: ImageFileType; end
add_image_file_format(".dummy", b"Dummy Image", Dummy, "dummy.jl")

# NRRD image format
type NRRDFile <: ImageFileType end
add_image_file_format(".nrrd", b"NRRD", NRRDFile, "nrrd.jl")
add_image_file_format(".nhdr", b"NRRD", NRRDFile, "nrrd.jl")

# Andor Technologies SIF file format
type AndorSIF <: Images.ImageFileType end
add_image_file_format(".sif", b"Andor Technology Multi-Channel File", AndorSIF, "SIF.jl")

# Imagine file format (http://holylab.wustl.edu, "Software" tab)
type ImagineFile <: ImageFileType end
add_image_file_format(".imagine", b"IMAGINE", ImagineFile, "Imagine.jl")

# PCO b16 image format
type B16File <: ImageFileType end
add_image_file_format(".b16", b"PCO-", B16File, "b16.jl")


text:
# an example on mixed integer programming
#
#   maximize x + 2 y + 5 z
#
#   s.t.  x + y + z <= 10
#         x + 2 y + z <= 15
#
#         x is continuous: 0 <= x <= 5
#         y is integer: 0 <= y <= 10
#         z is binary
#

using CPLEX

env = CPLEX.Env()

model = CPLEX.Model(env, "mip_01")
CPLEX.set_sense!(model, :Max)
CPLEX.set_param!(env, "CPX_PARAM_MIPDISPLAY",1)
CPLEX.set_param!(env, "CPX_PARAM_MIPINTERVAL",1)

CPLEX.add_var!(model, 1., 0., 5.)  # x
CPLEX.add_var!(model, 2., 0, 10)   # y
CPLEX.add_var!(model, 5., 0, 1)    # z
CPLEX.set_vartype!(model, ['C', 'I', 'B'])

CPLEX.add_constr!(model, ones(3), '<', 10.)
CPLEX.add_constr!(model, [1., 2., 1.], '<', 15.)

CPLEX.optimize!(model)

println("sol = $(CPLEX.get_solution(model))")
println("objv = $(CPLEX.get_objval(model))")

ext: jl


text:
function sylmat(g::AbstractVector{T}, u::AbstractVector, v::AbstractVector{T}) where {T<:Number}
#
#  sylmat generates the Sylvester resultant matrix
#    of polynomial g, u, v
#
    mg = length(g) - 1
    mu = length(u) - 1
    mv = length(v) - 1
    np = mg + mu
    nq = mg + mv

    m = mg + mu + mv
    n = np + nq

    A = zeros(T, n, m)

    for j = 1:mg
       A[j:j+mu,j] = u
       A[np+j:np+j+mv,j] = v
    end

    for j = 1:mu
       A[j:j+mg,j+mg] = g
    end

    for j = 1:mv
       A[np+j:np+j+mg, j+mg+mu] = g
    end
    A
end


text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

import Base: put!, wait, isready, take!, fetch

type DictChannel <: AbstractChannel
    d::Dict
    cond_take::Condition    # waiting for data to become available
    DictChannel() = new(Dict(), Condition())
end

function put!(D::DictChannel, k, v)
    D.d[k] = v
    notify(D.cond_take)
    D
end

function take!(D::DictChannel, k)
    v=fetch(D,k)
    delete!(D.d, k)
    v
end

isready(D::DictChannel) = length(D.d) > 1
isready(D::DictChannel, k) = haskey(D.d,k)
function fetch(D::DictChannel, k)
    wait(D,k)
    D.d[k]
end

function wait(D::DictChannel, k)
    while !isready(D, k)
        wait(D.cond_take)
    end
end

# Usage:

# RemoteRef to a DictChannel on worker pid
# dc_ref=RemoteRef(()->DictChannel(), pid)

# Test if there is any data
# isready(dc_ref)

# add
# put!(dc_ref, 1, 2)

# Test if key 1 exists
# isready(dc_ref, 1)

# fetch key 1
# fetch(dc_ref, 1)

# fetch and remove key 1
# take!(dc_ref, 1)

# wait for key 3 to be added
# wait(dc_ref, 3)

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
# ---
# title: 392. Is Subsequence
# id: problem392
# author: Tian Jun
# date: 2020-11-14
# difficulty: Easy
# categories: Binary Search, Dynamic Programming, Greedy
# link: <https://leetcode.com/problems/is-subsequence/description/>
# hidden: true
# ---
#
# Given a string **s** and a string **t** , check if **s** is subsequence of
# **t**.
#
# A subsequence of a string is a new string which is formed from the original
# string by deleting some (can be none) of the characters without disturbing the
# relative positions of the remaining characters. (ie, `"ace"` is a subsequence
# of `"abcde"` while `"aec"` is not).
#
# **Follow up:**
# If there are lots of incoming S, say S1, S2, ... , Sk where k >= 1B, and you
# want to check one by one to see if T has its subsequence. In this scenario,
# how would you change your code?
#
# **Credits:**
# Special thanks to [@pbrother](https://leetcode.com/pbrother/) for adding this
# problem and creating all test cases.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "abc", t = "ahbgdc"
#     Output: true
#
#
# **Example 2:**
#
#
#
#     Input: s = "axc", t = "ahbgdc"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `0 <= s.length <= 100`
#   * `0 <= t.length <= 10^4`
#   * Both strings consists only of lowercase characters.
#
#
## @lc code=start
using LeetCode

function is_subsequence(s::String, t::String)::Bool
    if length(s) == 0
        true
    else
        i, j = 1, 1
        while j <= length(t)
            if s[i] == t[j]
                i += 1
                if i > length(s)
                    return true
                end
            end
            j += 1
        end
        false
    end
end
## @lc code=end


text:
using Base: StackFrame
# Virtual Pilot - made with GameZero.jl
# by Felipe & Ricardo (rmsrosa@gmail.com)
# sound effects from https://mixkit.co/free-sound-effects/game/?page=3

WD = 1800 # actual sreen width
HT = 1200 # actual screen heigth
WIDTH = div(WD,2) # app width needs half of that - bug?
HEIGHT = div(HT,2) # app height needs half of that - bug?
BACKGROUND = colorant"#202020"

# background stuff

## mountains
MOUNTAIN_COLOR = colorant"#757575"
num_peaks = 20
mntpos = [1; accumulate(+, rand(200:600, num_peaks))]
mnt_max = mntpos[end]
mntpeak = rand(600:1000, num_peaks)
mountains = [
    Triangle(mntpos[n], 1200,
    div(mntpos[n] + mntpos[n+1], 2), mntpeak[n],
    mntpos[n+1], 1200)
    for n in 1:num_peaks-1
]

## blocks
BRICK_COLOR = colorant"#b08070"
num_bricks = 50
brkpos = [1; accumulate(+, rand(20:100, num_bricks))]
brk_max = brkpos[end]
brkht = rand(10:200, num_bricks)
bricks = [Rect(brkpos[n], 1200, brkpos[n+1]-brkpos[n], -brkht[n]) for n in 1:num_bricks]

## grass
GRASS_COLOR = colorant"#208020"
grass = rand(1:50, WIDTH)
grass_pos = 1

# Actors

## Space pod
space_pod = Actor("modfighter.png")

## Lasers
mutable struct Laser
    actor::Actor
    frame_fired::Int
end

## Enemies
mutable struct Enemy
    actor::Actor
    xvel::Int
    yvel::Int
    mobility::Int
    health::Int
end

## Info
game_name = TextActor(
    "Piloto Virtual", # text
    "moonhouse", # font
    font_size=96)
    #color=colorant"#b08070")
game_name.pos = (div(WD, 2) - div(WD, 5), div(HT, 10))

win_text = TextActor(
    "voce venceu!!!", # text
    "moonhouse", # font
    font_size=160)
    #color=colorant"#b08070")
win_text.pos = (div(WD, 2) - div(WD, 3), div(HT, 5))

loose_text = TextActor(
    "Voce perdeu...", # text
    "moonhouse", # font
    font_size=160)
    #color=colorant"#b08070")
loose_text.pos = (div(WD, 2) - div(WD, 3), div(HT, 5))

# weighted random function
weighted_rand(items::AbstractVector, weights::AbstractVector{<:Real}) =
    items[findfirst(cumsum(weights)/sum(weights) .> rand())]

# Spawn game
function spawn_game()
    global space_pod, pod_lasers
    global num_enemies, enemies

    global game_state = :playing
    global frame_num = 0
    global num_lasers = 10

    space_pod.pos = div(HT,20), div(HT,4) # starting x, y position
    pod_lasers = [
        Laser(
            Actor("pod_laser.png"), # image
            0 # frame_fired
        )
        for _ in 1:num_lasers
    ]

    enemies = []
    num_enemies = 10
    for n in 1:num_enemies
        enemytype = weighted_rand(1:2, [2,1])
        if enemytype == 1
            imagefilename = "enemyfighter4.png"
            xvel = rand(-4:-1)
            yvel = rand(-1:1)
            mobility = rand(1:20)
            health = 5
        elseif enemytype == 2
            imagefilename = "enemyfighter2b.png"
            xvel = rand(-6:-1)
            yvel = rand(-2:2)
            mobility = rand(1:5)
            health = 10
        end
        enemy = Enemy(
            Actor(imagefilename), # image
            xvel, yvel, mobility, health,
        )
        enemy.actor.x = div(WD,2) + 200*n
        enemy.actor.y = rand(20:20:900)
        push!(enemies, enemy)
    end
end

spawn_game()

function on_key_down(g::Game, key)
    global pod_lasers, frame_num, game_state
    laser_frames = getfield.(pod_lasers, :frame_fired)
    last_fire = maximum(laser_frames)
    charged = findfirst(iszero, laser_frames)
    if key == Keys.LSHIFT && charged !== nothing && frame_num > last_fire + 3
        pod_lasers[charged].frame_fired = frame_num
        pod_lasers[charged].actor.x = space_pod.x + div(space_pod.w,2)
        pod_lasers[charged].actor.y = space_pod.y + div(3space_pod.h,4)
        play_sound("mixkit-electronic-retro-block-hit-2185.wav")
        #play_sound("mixkit-video-game-power-up-3164.wav")
    end
    if game_state != :playing && key == Keys.RETURN
        game_state = :playing
        spawn_game()
    end
    # println(key)
end

function update(g::Game)
    global frame_num, pod_lasers
    global grass_pos, enemies, game_state, brk_max
    if game_state == :playing
        frame_num += 1
        if rem(frame_num, 4) == 0
            for n in eachindex(mountains)
                mountains[n].p1[1] -= 1
                mountains[n].p2[1] -= 1
                mountains[n].p3[1] -= 1
                if mountains[n].p3[1] < 1
                    mountains[n].p1[1] += mnt_max
                    mountains[n].p2[1] += mnt_max
                    mountains[n].p3[1] += mnt_max
                end
            end
        end
        for n in 1:length(bricks)
            bricks[n].x -= 1
            if bricks[n].x < -100
                bricks[n].x += brk_max
            end
        end
        grass_pos += ifelse(grass_pos == 900, -899, 1)
        for enemy in enemies
            if rand(1:100) ≤ enemy.mobility
                enemy.yvel = -enemy.yvel
            end
            enemy.actor.x += enemy.xvel
            enemy.actor.y -= enemy.xvel * enemy.yvel
            if enemy.actor.x < -192
                enemy.actor.x = WD
            end
            enemy.actor.y = max(20, min(div(7*HT,10), enemy.actor.y))

        end
        #space_pod move
        if g.keyboard.DOWN && space_pod.y < 940
            space_pod.y += 4
        elseif g.keyboard.UP && space_pod.y > 20
            space_pod.y -= 4
        end
        if g.keyboard.RIGHT && space_pod.x < 200
            space_pod.x += 2
        elseif g.keyboard.LEFT && space_pod.x > 20
            space_pod.x -= 2
        end
        for pod_laser in pod_lasers
            if pod_laser.frame_fired > 0
                pod_laser.actor.x += 6
                if pod_laser.actor.x > WD
                    pod_laser.frame_fired = 0
                end
            end
        end
        for enemy in filter(x -> x.health > 0, enemies)
            if enemy.health == 10
                enemy.actor.image = mod(div(frame_num, 10), 2) == 0 ?
                    "enemyfighter2a.png" : "enemyfighter2b.png"
            end
            if collide(enemy.actor, space_pod)
                game_state = :loose
                play_sound("mixkit-explosion-hit-1704.wav")
            end
            for pod_laser in filter(l -> l.frame_fired > 0, pod_lasers)
                if collide(enemy.actor, pod_laser.actor) == true
                    enemy.health -= 5
                    if enemy.health > 0
                        enemy.actor.image = "enemyfighter2.png"
                    end
                    pod_laser.frame_fired = 0
                    play_sound("mixkit-shot-light-explosion-1682.wav")
                end
            end
        end
        if maximum(getfield.(enemies, :health)) == 0
            game_state = :win
        end
    end
end

function draw(g::Game)
    if frame_num < 50
        draw(game_name)
    end
    for mnt in mountains
        draw(mnt, MOUNTAIN_COLOR, fill=true)
    end
    for brick in bricks
        draw(brick, BRICK_COLOR, fill=true)
    end
    for n in 0:WIDTH
        np = mod(frame_num + n, WIDTH) + 1
        draw(Rect(2n-1, 1200, 2, -grass[np]), GRASS_COLOR, fill=true)
    end
    draw(space_pod)
    for enemy in enemies
        if enemy.health > 0
            draw(enemy.actor)
        end
    end
    for pod_laser in pod_lasers
        if pod_laser.frame_fired > 0
            draw(pod_laser.actor)
        end
    end
    if game_state == :win
        draw(win_text)
    elseif game_state == :loose
        draw(loose_text)
    end
end


text:
#= Resampling function -- takes a vector of weights that are < alpha
# to resample and does SRC resampling, then returns a vector with resampled
 weights
=#
function resample(to.resample,alpha){

  log.alpha <- log(alpha)
  log.u <- log.alpha + log( runif(1) )
  k <- 1

  while (k <= length(to.resample) ){

    # want to find u <- u - w then if u <= 0
    # u <- u + alpha
    # IN OUR CASE IF u <= w
    # find u <- u + alpha - w
    # as working with logs
    # SO first check whether log(u) < log(w)

    if ( log.u < to.resample[k] ){
      # u <- u + alpha - w
      # first find log(alpha-w) label as log.alpha.weight
      temp = c( log.alpha ,  to.resample[k] )
      c = max(temp)
      log.alpha.weight <- c + log( exp( temp[1] - c ) - exp( temp[2] - c )  )
      temp = c( log.u , log.alpha.weight )
      c = max(temp)
      log.u <- c + log( sum( exp( temp - c ) ) )
      to.resample[k] <- log.alpha
    }

    else{
      # u <- u - w
      # and w is not resampled, given weight 0 (-Inf for log(w))
      temp <- c( log.u , to.resample[k] )
      c = max(temp)
      log.u <- c + log( exp( temp[1] - c ) - exp( temp[2] - c )  )
      to.resample[k] <- - Inf

    }

    k <- k+1

  }

  return(to.resample)

}


text:
using JuliaDB, DataFrames, IterableTables, Gadfly
using JuliaDB: groupby

data = loadtable("results/bench.csv")
data = reindex(data, :Benchmark)

baseline = filter(x->x.Package=="Baseline", data)
data = filter(x->x.Package != "Baseline", data)

normalised = join(data, baseline) do a, b
  (Package = a.Package, Time = a.Time / b.Time)
end

geomean(xs) = prod(xs)^(1/length(xs))

normalised = join(normalised, groupby(geomean, normalised, :Package, select = :Time),
                  lkey = :Package)

df = sort(DataFrame(normalised), cols = [:geomean, :Time])

p = plot(df, x = :Package, y = :Time, color = :Benchmark,
         Guide.ylabel(nothing), Guide.xlabel(nothing), Scale.y_log10)

draw(SVG(joinpath("results/benchmarks.svg"), 8inch, 8inch/1.6), p)


text:
module ChemicalReactionNetworks

export
    Reaction,
    ReactionNetwork,
    SpeciesComplex,
    ==,
    hash,
    isredundant,
    stoichiometric_matrix,
    mass_action,
    mass_action!,
    reaction_current,
    jacobian,
    from_stoichiometric_matrix,
    normal_form,
    conservation_laws,
    cycles,
    n_conservation_laws,
    n_cycles,
    rand_concentrations,
    chemostatted!,
    relabel!,
    complexes,
    complex_decomposition,
    deficiency,
    complex_concentration,
    equilibrium_state,
    equilibrium_state_space,
    cycle_affinities,
    specificity,
    complex_basis_reactions,
    free_energies,
    reaction_rates,
    petri_net

include("crn.jl")
include("algebra.jl")
include("sampling.jl")
include("chemostat.jl")
include("complexes.jl")
include("thermodynamics.jl")
include("petri.jl")

"""
A Julia package for studying and simulating Chemical
Reaction Networks, of the kind that would be encountered in
the study of complex metabolic networks, self-replication,
and so on.

API Overview:
-`Reaction(...)` creates a reaction.
-`ReactionNetwork(...)` creates a network of reactions.
-`mass_action(reactions, concentration)` returns change in concentration over time.
-`chemostatted!(reactions, chemostatted_species, chemostatted_concentrations)` chemostats certain species and modifies the network.
-`equilibrium_state(reaction_net)` returns detailed-balance equilibrium state.
-`conservation_laws(reaction_net)` returns conservation laws of net
-`cycles(reaction_net)` returns cycles of net.
-`deficiency(reaction_net)` returns network deficiency.
-`complex_decomposition(reaction_net)` decomposes a CRN in terms of a set of complexes and incidence matrix.
"""

end


text:
export StateSpace, propagate, observe

"""
$(TYPEDEF)

An immutable structure representing the dynamical `f` and observation `h` operators.
The dynamical model is provided by the right hand side of the ODE to solve.
For a system of ODEs, we will prefer an in-place syntax `f(du, u, p, t)`, where `p` are parameters of the model.
We rely on `OrdinaryDiffEq` to integrate the dynamical system with the Tsitouras 5/4 Runge-Kutta method adaptive time marching.
`h` must be a function of the form `h(u, t)`, where `u` is the state vector and `t` is the time.

## Fields
$(TYPEDFIELDS)

"""
struct StateSpace
    "Propagatation f"
    f::Function

    "Observation h"
    h::Function
end

# To remove
function propagate(f::Function, X, t::Float64, Ny::Int64, Nx::Int64; P::Parallel=serial)
    Nypx, Ne = size(X)
    @assert Nypx == Ny + Nx "Wrong dimension of Ny or Nx"
    if typeof(P)==Serial
        @inbounds for i=1:Ne
        x = view(X, Ny+1:Nypx, i)
        x .= f(x, t)[2]
        end
    else
        @inbounds Threads.@threads for i=1:Ne
        x = view(X, Ny+1:Nypx, i)
        x .= f(x, t)[2]
        end
    end
end

"""
$(TYPEDSIGNATURES)

Evaluate the function `h` for the different state vectors of the `X` at time `t`, and store the results in the first `Ny` columns of `X`.
`X` is an ensemble matrix that contains the observation vectors in the first `Ny` lines, and the state vectors in the lines `Ny+1` to `Ny+Nx`.
The code can run in serial or with multithreading.
"""
function observe(h::Function, X, t::Float64, Ny::Int64, Nx::Int64; P::Parallel=serial)
    Nypx, Ne = size(X)
    @assert Nypx == Ny + Nx "Wrong dimension of Ny or Nx"
    if typeof(P)==Serial
        @inbounds for i=1:Ne
        x = view(X, Ny+1:Nypx, i)
        y = view(X, 1:Ny, i)

        y .= h(x, t)
        end
    elseif typeof(P)==Thread
        @inbounds Threads.@threads for i=1:Ne
        x = view(X, Ny+1:Nypx, i)
        y = view(X, 1:Ny, i)

        y .= h(x, t)
        end
    end
end

"""
$(TYPEDSIGNATURES)

Apply the observation operator of the `StateSpace` `F` to the ensemble matrix `X` at time `t`.
"""
observe(F::StateSpace, x::Array{Float64,1}, t::Float64) = F.h(x, t)


text:
"""
TimeTrees
---------

Basic rooted phylogenetic time tree manipulation module which includes TimeTree
    and Node classes.  The TimeTree class can be initialized from a Newick
    string and provides methods for producing visualizations using ASCII art.

* `TimeTree`: a rooted phylogenetic time tree.
* `Node`: the basic building block of a phylogenetic tree.

Read the documentation for these types for further information.
"""
module TimeTrees

export Node, isRoot, isLeaf, addChild!, edgeLength,
    getDecendentCount, getSorted, getCopy,
    TimeTree, getLeaves, getNodes, getInternalNodes,
    getLeafCount, getNodeCount,
    getNewick, plot

import Base.show

"""
The `Node` type is the basic element of any `TimeTree`. It
defines the following attributes:

* `parent`: the parent node. If the node is a root, then this
attribute refers to the current node object.

* `children`: the list of child nodes. If empty, this node is
a leaf.

* `age`: the age or age of the node (usually measured
from the most recent sample/leaf node.

* `label`: a (possibly empty) string labeling this node. Ofen
only used for leaves.

* `number`: an integer used to uniquely identify the node in a tree.
    This is set to -1 by default.

An empty node can be constructed using the `Node()` method.

Several methods on nodes are defined: `isRoot`, `isLeaf`, `addChild!`,
`edgeLength`, `getLeaves`, `getInternalNodes`, `getSorted`, `getNewick`,
`getCopy`.  Read the documentation for these methods for further information.
"""
type Node
    parent::Node
    children::Array{Node, 1}
    age::Float64
    label::AbstractString
    annotation::Dict{AbstractString,AbstractString}
    number::Int

    Node() = begin
        n = new()
        n.parent = n
        n.children = []
        n.age = 0.0
        n.label = ""
        n.annotation = Dict{AbstractString,AbstractString}()
        n.number = -1
        return n
    end
end

function show(io::IO, n::Node)
    if isRoot(n)
        print(io, "Root")
    else
        print(io, "Non-root")
    end

    print(io, " node (age: $(n.age), children: $(length(n.children))")

    if length(n.label)>0
        print(io,", label: $(n.label))")
    else
        print(io,", no label)")
    end

end

"""
`getNewick(n::Node)` retrieves the Newick representation of the subtree below `n`.
"""
function getNewick(n::Node)
    res = ""
    if length(n.children)>0
        res = string(res, "(")
        for i in 1:length(n.children)
            if i > 1
                res = string(res, ",")
            end
            res = string(res, getNewick(n.children[i]))
        end
        res = string(res, ")")
    end
    res = string(res, n.label, ":", edgeLength(n))
end

"""
`isRoot(n::Node)` returns `true` if `n` is a root node.
"""
isRoot(n::Node) = n.parent == n

"""
`isLeaf(n::Node)` returns `true` if `n` is a leaf node.
"""
isLeaf(n::Node) = length(n.children) == 0

"""
`addChild!(n::Node, c::Node)` adds node `c` as a child of node `n`.
"""
function addChild!(n::Node, c::Node)
    c.parent = n
    push!(n.children, c)
end

"""
`edgeLength(n::Node)` returns the length of the edge above node `n`.  Always
returns 0 if `n` is a root node.
"""
function edgeLength(n::Node)
    if isRoot(n)
        return 0.0
    else
        return n.parent.age - n.age
    end
end

"""
`getLeaves(n::Node)` returns array of leaf nodes below `n`.
"""
function getLeaves(n::Node)
    if isLeaf(n)
        return [n]
    else
        res = []
        for c in n.children
            res = [res; getLeaves(c)]
        end

        return res
    end
end

"""
`getInternalNodes(n::Node)` returns array of internal nodes below (and including) node `n`.
"""
function getInternalNodes(n::Node)
    if isLeaf(n)
        return []
    else
        res = [n]
        for c in n.children
            res = [res; getInternalNodes(c)]
        end

        return res
    end
end

"""
`getDecendentCount(n::Node)` returns number of decendents below (and including) `n`.
"""
function getDecendentCount(n::Node)
    res = 1
    for c in n.children
        res += getDecendentCount(c)
    end

    return res
end

"""
`getSorted(n::Node, [rev=false])` produces a copy of the clade below `n`,
sorted so that children appear in order of the number of decendents they have.
Setting `rev=true` reverses the sort.
"""
function getSorted(n::Node; rev = false)

    copy = Node()
    copy.label = n.label
    copy.age = n.age

    for c in n.children
        addChild!(copy, getSorted(c))
    end

    perm = sortperm([getDecendentCount(c) for c in copy.children], rev=rev)
    copy.children = copy.children[perm]

    return copy
end

"""
`getCopy(n::Node)` produces a deep copy of the clade below `n`.
"""
function getCopy(n::Node)
    copy = Node()
    copy.label = n.label
    copy.age = n.age

    for c in n.children
        addChild!(copy, getCopy(c))
    end

    return copy
end

"""
The `TimeTree` type represents a rooted phylogenetic tree.  It defines a
single attribute `root::Node` which specifies the root node of the tree.

`TimeTree` objects can be constructed in two ways:

1. Using the constructor `TimeTree(root::Node)` to create a tree
with the given root node, and

2. Using the constructor `TimeTree(newick::String)` which creates a
tree by parsing a Newick string given as an argument.  For instance,
`TimeTree("((A:1,B:1):1,C:2):0;")` constructs a tree with 3 leaves and
a total age of two.

The root node of a `TimeTree` is accessible via the `TimeTree`'s `root`
attribute.

There are several methods which act on trees: `getLeaves`,
`getInternalNodes`, `getNodes`, `getLeafCount`, `getNodeCount`,
`getSorted`, `getCopy`, `getNewick` and `plot`.
"""
type TimeTree
    root::Node
    nodes::Array{Node,1}
    nLeaves::Int
    nNodes::Int
end

function TimeTree(root::Node)
    leaves = getLeaves(root)
    internals = getInternalNodes(root)

    nodes = [leaves; internals]
    for (i, node) in enumerate(nodes)
        node.number = i
    end

    TimeTree(root, nodes, length(leaves), length(nodes))
end

"""
`getLeaves(t::TimeTree)` returns array of leaf nodes that `t` contains.
"""
getLeaves(t::TimeTree) = slice(t.nodes, 1:t.nLeaves)

"""
`getLeafCount(t::TimeTree)` returns the number of leaves in the tree.
"""
getLeafCount(t::TimeTree) = t.nLeaves

"""
`getNodeCount(t::TimeTree)` returns the number of nodes in the tree.
"""
getNodeCount(t::TimeTree) = t.nNodes

"""
`getInternalNodes(t::TimeTree)` returns an array of internal nodes that `t` contains.
"""
getInternalNodes(t::TimeTree) = slice(t.nodes, (t.nLeaves+1):t.nNodes)

"""
`getNodes(t::TimeTree)` returns array of nodes that `t` contains.
"""
getNodes(t::TimeTree) = t.nodes

"""
`getSorted(t::TimeTree, [rev=false])` returns a copy of `t` in which children are
sorted in order of the number of of their respective decendents.
"""
getSorted(t::TimeTree; rev = false) = TimeTree(getSorted(t.root))

"""
`getCopy(t::TimeTree)` returns a deep copy of `t`.
"""
getCopy(t::TimeTree) = TimeTree(getCopy(t.root))

function show(io::IO, t::TimeTree)
    print(io, string("A phylogenetic tree with ",
        length(getLeaves(t)), " leaves (",
        length(getNodes(t)), " nodes in total)"))
end

"""
`getNewick(t::TimeTree)` retrieves the Newick representation of `t`.
"""
getNewick(t::TimeTree) = string(getNewick(t.root), ";")

"""
Construct a `TimeTree` from a Newick string.
"""
function TimeTree(newick::AbstractString)

    ### Parser ###

    patterns = Dict{AbstractString,Regex}(
        "open_paren" => r"^\(",
        "close_paren" => r"^\)",
        "open_an" => r"^\[&",
        "close_an" => r"^\]",
        "eq" => r"^=",
        "colon" => r"^:",
        "comma" => r"^,",
        "number" => r"^\d+(\.\d*)?([eE]-?\d+)?",
        "string" => r"^((\w+)|(\"[^\"]*\")|(^'[^']*'))"
    )

    function matchToken(token::AbstractString; mustMatch = false)

        # Skip whitespace
        while i<=length(newick) && (newick[i] == ' ' || newick[i] == '\t')
            i += 1
        end

        m = i<=length(newick) ? match(patterns[token], newick[i:end]) : nothing
        if m == nothing
            if mustMatch
                c = newick[i]
                error("Expected token $token at index $i but got $c instead.")
            else
                return nothing
            end
        end

        i += length(m.match)
        return m.match
    end

    function ruleT()
       return TimeTree(ruleN())
    end

    function ruleN()
        node = Node()
        if matchToken("open_paren") != nothing
            while true
                addChild!(node, ruleN())
                if matchToken("comma") == nothing
                    break
                end
            end
            matchToken("close_paren", mustMatch = true)
        end

        node.label = ruleL()
        node.annotation = ruleA()
        node.age = ruleH()

        return node
    end

    function ruleL()
        res = matchToken("string")
        if res == nothing
            return ""
        end

        if startswith(res, "\"") || startswith(res, "'")
            return res[2:end-1]
        else
            return res
        end
    end

    function ruleA()
        res = Dict{AbstractString, AbstractString}()

        if matchToken("open_an") == nothing
            return res
        end

        while true
            key = matchToken("number")
            if key == nothing
                key = matchToken("string", mustMatch = true)
            end

            matchToken("eq", mustMatch = true)

            value = matchToken("number")
            if value == nothing
                value = matchToken("string", mustMatch = true)
            end

            res[key] = value

            if matchToken("comma") == nothing
                break
            end
        end

        matchToken("close_an", mustMatch = true)

        return res
    end

    function ruleH()
        res = matchToken("colon")
        if res == nothing
            return 0.0
        end

        return float(matchToken("number", mustMatch = true))
    end

    i = 1
    tree =  ruleT()

    ### Post-processing ###

    function getTime(node::Node, currentTime::AbstractFloat)
        currentTime += node.age

        maxTime = currentTime
        for child in node.children
            maxTime = max(maxTime, getTime(child, currentTime))
        end

        return maxTime
    end

    function branchLengthToHeight(node::Node, treeHeight::AbstractFloat, currentTime::AbstractFloat)
        currentTime += node.age
        node.age = treeHeight - currentTime

        for child in node.children
            branchLengthToHeight(child, treeHeight, currentTime)
        end
    end

    age = getTime(tree.root, 0.0)
    branchLengthToHeight(tree.root, age, 0.0)

    return tree
end

"""
`plot(t::TimeTree)` produces an ASCII representation of `t`.
Optional arguments are

* `width`: the number of columns used for the plot (default 70)
* `labelLeaves`: whether to label the leaves (default true)
* `dots`: whether to connect leaves to labels using dots (default true)
"""
function plot(t::TimeTree; width = 70, labelLeaves = true, dots = true)

    leaves = getLeaves(t)
    nodes = getNodes(t)
    pos = zeros(length(nodes))

    grid = fill(' ', length(leaves), width)

    function computePos(n::Node)
        idx = findfirst(nodes, n)
        if isLeaf(n)
            pos[idx] = findfirst(leaves, n)
        else
            for c in n.children
                pos[idx] += computePos(c)
            end
            pos[idx] /= length(n.children)
        end
        return pos[idx]
    end
    computePos(t.root)

    function ageToIdx(h::Float64)
        return round(Int, h/t.root.age*(width-1) + 1)
    end

    # Edges
    for i in 1:length(nodes)
        n = nodes[i]

        if isRoot(n)
            continue
        else
            pi = findfirst(nodes, n.parent)
            x1 = ageToIdx(nodes[i].age)
            y1 = round(Int, pos[i])
            x2 = ageToIdx(nodes[pi].age)
            y2 = round(Int, pos[pi])
            ymin = min(y1, y2)
            ymax = max(y1, y2)

            grid[y1,x1:x2] = '-'
            grid[ymin:ymax, x2] = '|'
            grid[ymin,x2] = '/'
            grid[ymax,x2] = '\\'
        end
    end

    # Nodes
    for i in 1:length(nodes)
        x = ageToIdx(nodes[i].age)
        y = round(Int, pos[i])
        if isLeaf(nodes[i])
            grid[y, x] = '*'
            if x>1 && labelLeaves && dots
                grid[y,x-1:-1:1] = '⋅'
            end
        else
            grid[y, x] = '+'
        end
    end

    # Display
    for i in 1:length(leaves)
        if labelLeaves
            println(string(join(reverse(vec(grid[i,:]))), " ", leaves[i].label))
        else
            println(string(join(reverse(vec(grid[i,:])))))
        end
    end

end

end


text:
const SIMPLE_TYPES = [
    Cint,
    Cstring,
    Cvoid,
    Int8,
    UInt8,
    Int16,
    UInt16,
    Int32,
    UInt32,
    bool,
]

function is_simple_type(fld_tp)
    return fld_tp in SIMPLE_TYPES
end

function is_simple_type(fld_tp::Type{Ptr{T}}) where T
    return is_simple_type(T)
end

function simple_type_value(fld)
    return fld
end

function simple_type_value(fld::Ptr{Cvoid})
    return fld
end

function simple_type_value(fld::Cstring)
    return fld == C_NULL ? "" : convert_cstring_to_str(fld)
end

function simple_type_value(fld::Ptr{Node})
    return fld == C_NULL ? fld : _convert_to_proper_node_type(fld)
end

function simple_type_value(fld::Ptr{T}) where T
    return fld == C_NULL ? fld : simple_type_value(unsafe_load(fld))
end

function convert_to_actual_value(x::Value)
    if x.type == T_Integer
        return convert(Int64, x.val)
    elseif x.type == T_Float || x.type == T_String || x.type == T_BitString
        val = simple_type_value(convert(Cstring, x.val))
        if x.type == T_Float
            return Base.parse(Float64, val)
        elseif x.type == T_String || x.type == T_BitString
            return val
        end
    elseif x.type == T_Null
        return nothing
    else
        throw("not implemented")
    end
end

function _convert_to_proper_node_type(node)
    node_type_enum = unsafe_load(convert(Ptr{Node}, node)).type
    if node_type_enum in [T_Integer, T_Float, T_String, T_BitString, T_Null]
        unsafe_load(convert(Ptr{Value}, node))
    else
        node_type = convert_to_node_type(Val(node_type_enum))
        unsafe_load(convert(Ptr{node_type}, node))
    end
end

############################################
# `AbstractVector{Any}` interface for `List`
############################################

function Base.iterate(lst::List, state=(1, C_NULL))
    (state_idx, state_ptr) = state
    state_idx > lst.length && return nothing
    currentNode = state_idx == 1 ? lst.head : state_ptr
    state_idx == lst.length && @assert currentNode == lst.tail "The last item should be equal to tail"
    currentListCell = unsafe_load(currentNode)

    return (_convert_to_proper_node_type(currentListCell.data.ptr_value), (state_idx+1, currentListCell.next))
end

Base.IteratorSize(::Type{List}) = Base.HasLength()
Base.IteratorEltype(::Type{List}) = Base.EltypeUnknown()
Base.length(lst::List) = lst.length
Base.size(lst::List) = (Base.length(lst),)

function Base.getindex(lst::List, i::Int)
    1 <= i <= lst.length || throw(BoundsError(lst, i))
    currentNode = lst.head
    for idx in 2:i
        currentListCell = unsafe_load(currentNode)
        currentNode = currentListCell.next
    end
    currentListCell = unsafe_load(currentNode)
    return _convert_to_proper_node_type(currentListCell.data.ptr_value)
end


text:
function constraint_active_load_gen_aggregation(pm::GenericPowerModel, n::Int, i::Int, bus_gens, pd)
    pnode = PowerModels.var(pm, n, :pnode)[i]
    pg = PowerModels.var(pm, n, :pg)

    @constraint(pm.model, pnode == sum(pg[g] for g in bus_gens) - pd)
end

function constraint_reactive_load_gen_aggregation(pm::GenericPowerModel, n::Int, i::Int, bus_gens, qd)
    qnode = PowerModels.var(pm, n, :qnode)[i]
    qg = PowerModels.var(pm, n, :qg)

    @constraint(pm.model, qnode == sum(qg[g] for g in bus_gens) - qd)
end

""
function constraint_active_load_gen_aggregation_sheddable(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, bus_gens, bus_loads)
    pnode = PowerModels.var(pm, n, cnd, :pnode, i)
    pg = PowerModels.var(pm, n, cnd, :pg)
    pl = PowerModels.var(pm, n, cnd, :pl)

    @constraint(pm.model, pnode == sum(pg[g] for g in bus_gens) - sum(pl[l] for l in bus_loads))
end

function constraint_reactive_load_gen_aggregation_sheddable(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, bus_gens, bus_loads)
    qnode = PowerModels.var(pm, n, cnd, :qnode, i)
    qg = PowerModels.var(pm, n, cnd, :qg)
    ql = PowerModels.var(pm, n, cnd, :ql)

    @constraint(pm.model, qnode == sum(qg[g] for g in bus_gens) - sum(ql[l] for l in bus_loads))
end

""
function constraint_flexible_active_load(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, prated, pref)
    pl_delta = PowerModels.var(pm, n, cnd, :pl_delta, i)
    load_ind = PowerModels.var(pm, n, cnd, :load_ind, i)

    @constraint(pm.model,   pl_delta <= 2 * prated * load_ind)
end

function constraint_flexible_reactive_load(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, qrated, qref)
    ql_delta = PowerModels.var(pm, n, cnd, :ql_delta, i)
    load_ind = PowerModels.var(pm, n, cnd, :load_ind, i)

    @constraint(pm.model,   ql_delta <= 2 * qrated * load_ind)
end
# function constraint_flexible_active_load(pm::GenericPowerModel, n::Int, i::Int, prated, pref)
#     pl_delta = PowerModels.var(pm, n, :pl_delta)[i]
#     load_ind = PowerModels.var(pm, n, :load_ind)[i]
#
#     @constraint(pm.model,   pl_delta <= 2 * prated * load_ind)
# end
#
# function constraint_flexible_reactive_load(pm::GenericPowerModel, n::Int, i::Int, qrated, qref)
#     ql_delta = PowerModels.var(pm, n, :ql_delta)[i]
#     load_ind = PowerModels.var(pm, n, :load_ind)[i]
#
#     @constraint(pm.model,   ql_delta <= 2 * qrated * load_ind)
# end

""
function constraint_fixed_active_load(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, pref::AbstractFloat)
    pl_delta = PowerModels.var(pm, n, cnd, :pl_delta, i)
    pl = PowerModels.var(pm, n, cnd, :pl, i)

    @constraint(pm.model, pl_delta == 0)
    @constraint(pm.model, pl == pref)
end

function constraint_fixed_reactive_load(pm::GenericPowerModel, n::Int,  cnd::Int, i::Int, qref::AbstractFloat)
    ql_delta = PowerModels.var(pm, n, cnd, :ql_delta, i)
    ql = PowerModels.var(pm, n, cnd, :ql, i)

    @constraint(pm.model, ql_delta == 0)
    @constraint(pm.model, ql == qref)
end

# function constraint_fixed_active_load(pm::GenericPowerModel, n::Int, i::Int, pref::AbstractFloat)
#     pl_delta = PowerModels.var(pm, n, :pl_delta)[i]
#     pl = PowerModels.var(pm, n, :pl)[i]
#
#     @constraint(pm.model, pl_delta == 0)
#     @constraint(pm.model, pl == pref)
# end
#
# function constraint_fixed_reactive_load(pm::GenericPowerModel, n::Int, i::Int, qref::AbstractFloat)
#     ql_delta = PowerModels.var(pm, n, :ql_delta)[i]
#     ql = PowerModels.var(pm, n, :ql)[i]
#
#     @constraint(pm.model, ql_delta == 0)
#     @constraint(pm.model, ql == qref)
# end
""
function constraint_flexible_active_gen(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, prated, pref)
    pg_delta = PowerModels.var(pm, n, cnd, :pg_delta, i)
    gen_ind = PowerModels.var(pm, n, cnd, :gen_ind, i)

    @constraint(pm.model,   pg_delta <= 2 * prated * gen_ind)
end

function constraint_flexible_reactive_gen(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, qrated, qref)
    qg_delta = PowerModels.var(pm, n, cnd, :qg_delta, i)
    gen_ind = PowerModels.var(pm, n, cnd, :gen_ind, i)

    @constraint(pm.model,   qg_delta <= 2 * qrated * gen_ind)
end
# function constraint_flexible_active_gen(pm::GenericPowerModel, n::Int, i::Int, prated, pref)
#     pg_delta = PowerModels.var(pm, n, :pg_delta)[i]
#     gen_ind = PowerModels.var(pm, n, :gen_ind)[i]
#
#     @constraint(pm.model,   pg_delta <= 2 * prated * gen_ind)
# end
#
# function constraint_flexible_reactive_gen(pm::GenericPowerModel, n::Int, i::Int, qrated, qref)
#     qg_delta = PowerModels.var(pm, n, :qg_delta)[i]
#     gen_ind = PowerModels.var(pm, n, :gen_ind)[i]
#
#     @constraint(pm.model,   qg_delta <= 2 * qrated * gen_ind)
# end
# function constraint_redispatch_active_power_gen(pm::GenericPowerModel, n::Int, i::Int, pref::AbstractFloat)
#     pg_delta = PowerModels.var(pm, n, :pg_delta)[i]
#     pg = PowerModels.var(pm, n, :pg)[i]
#
#     @constraint(pm.model, pg_delta >= pg - pref)
#     @constraint(pm.model, pg_delta >= -(pg - pref))
# end
#
# function constraint_redispatch_reactive_power_gen(pm::GenericPowerModel, n::Int, i::Int, qref::AbstractFloat)
#     qg_delta = PowerModels.var(pm, n, :qg_delta)[i]
#     qg = PowerModels.var(pm, n, :qg)[i]
#
#     @constraint(pm.model, qg_delta >= qg - qref)
#     @constraint(pm.model, qg_delta >= -(qg - qref))
# end
# function constraint_second_stage_redispatch_active_power_gen(pm::GenericPowerModel, n::Int, i::Int, first_stage_network_id)
#     pg_delta = PowerModels.var(pm, n, :pg_delta)[i]
#     pg = PowerModels.var(pm, n, :pg)[i]
#     pg_first_stage = PowerModels.var(pm, first_stage_network_id, :pg)[i]
#
#     @constraint(pm.model, pg_delta >= pg - pg_first_stage)
#     @constraint(pm.model, pg_delta >= -(pg - pg_first_stage))
# end
#
# function constraint_second_stage_redispatch_reactive_power_gen(pm::GenericPowerModel, n::Int, i::Int, first_stage_network_id)
#     qg_delta = PowerModels.var(pm, n, :qg_delta)[i]
#     qg = PowerModels.var(pm, n, :qg)[i]
#     qg_first_stage = PowerModels.var(pm, first_stage_network_id, :qg)[i]
#
#     @constraint(pm.model, qg_delta >= qg - qg_first_stage)
#     @constraint(pm.model, qg_delta >= -(qg - qg_first_stage))
# end
# function constraint_redispatch_active_power_load(pm::GenericPowerModel, n::Int, i::Int, pref::AbstractFloat)
#     pl_delta = PowerModels.var(pm, n, :pl_delta)[i]
#     pl = PowerModels.var(pm, n, :pl)[i]
#
#     @constraint(pm.model, pl_delta >= pl - pref)
#     @constraint(pm.model, pl_delta >= -(pl - pref))
# end
#
# function constraint_redispatch_reactive_power_load(pm::GenericPowerModel, n::Int, i::Int, qref::AbstractFloat)
#     ql_delta = PowerModels.var(pm, n, :ql_delta)[i]
#     ql = PowerModels.var(pm, n, :ql)[i]
#
#     @constraint(pm.model, ql_delta >= ql - qref)
#     @constraint(pm.model, ql_delta >= -(ql - qref))
# end
# function constraint_second_stage_redispatch_active_power_load(pm::GenericPowerModel, n::Int, i::Int, first_stage_network_id)
#     pl_delta = PowerModels.var(pm, n, :pl_delta)[i]
#     pl = PowerModels.var(pm, n, :pl)[i]
#     pl_first_stage = PowerModels.var(pm, first_stage_network_id, :pl)[i]
#
#     @constraint(pm.model, pl_delta >= pl - pl_first_stage)
#     @constraint(pm.model, pl_delta >= -(pl - pl_first_stage))
# end
#
# function constraint_second_stage_redispatch_reactive_power_load(pm::GenericPowerModel, n::Int, i::Int, first_stage_network_id)
#     ql_delta = PowerModels.var(pm, n, :ql_delta)[i]
#     ql = PowerModels.var(pm, n, :ql)[i]
#     ql_first_stage = PowerModels.var(pm, first_stage_network_id, :ql)[i]
#
#     @constraint(pm.model, ql_delta >= ql - ql_first_stage)
#     @constraint(pm.model, ql_delta >= -(ql - ql_first_stage))
# end

"""
Q = P * tan(load_angle): To avoid division by 0: Q*cos(load_angle) = P*sin(load_angle)
"""
function constraint_tan_phi_load(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, load_angle)
    ql = PowerModels.var(pm, n, cnd, :ql, i)
    pl = PowerModels.var(pm, n, cnd, :pl, i)

    @constraint(pm.model, ql * cos(load_angle) == pl * sin(load_angle))
end
# function constraint_tan_phi_load(pm::GenericPowerModel, n::Int, i::Int, load_angle)
#     ql = PowerModels.var(pm, n, :ql)[i]
#     pl = PowerModels.var(pm, n, :pl)[i]
#
#     @constraint(pm.model, ql * cos(load_angle) == pl * sin(load_angle))
# end

""
function constraint_active_power_gen_contingency(pm::GenericPowerModel, n::Int, cnd::Int, i::Int)
    pg_delta = PowerModels.var(pm, n, cnd, :pg_delta, i)
    pg = PowerModels.var(pm, n, cnd, :pg, i)

    @constraint(pm.model, pg_delta == 0)
    @constraint(pm.model, pg == 0)
end

function constraint_reactive_power_gen_contingency(pm::GenericPowerModel, n::Int, cnd::Int, i::Int)
    qg_delta = PowerModels.var(pm, n, cnd, :qg_delta, i)
    qg = PowerModels.var(pm, n, cnd, :qg, i)

    @constraint(pm.model, qg_delta == 0)
    @constraint(pm.model, qg == 0)
end
# function constraint_active_power_gen_contingency(pm::GenericPowerModel, n::Int, i::Int)
#     pg_delta = PowerModels.var(pm, n, :pg_delta)[i]
#     pg = PowerModels.var(pm, n, :pg)[i]
#
#     @constraint(pm.model, pg_delta == 0)
#     @constraint(pm.model, pg == 0)
# end
#
# function constraint_reactive_power_gen_contingency(pm::GenericPowerModel, n::Int, i::Int)
#     qg_delta = PowerModels.var(pm, n, :qg_delta)[i]
#     qg = PowerModels.var(pm, n, :qg)[i]
#
#     @constraint(pm.model, qg_delta == 0)
#     @constraint(pm.model, qg == 0)
# end

""
function constraint_active_power_branch_contingency(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, f_idx, t_idx)
    p_fr = PowerModels.var(pm, n, cnd, :p, f_idx)
    p_to = PowerModels.var(pm, n, cnd, :p, t_idx)
    va_shift_fr = PowerModels.var(pm, n, cnd, :va_shift, f_idx)
    va_shift_to = PowerModels.var(pm, n, cnd, :va_shift, t_idx)

    @constraint(pm.model, p_fr == 0)
    @constraint(pm.model, p_to == 0)
    @constraint(pm.model, va_shift_fr  == 0)
    @constraint(pm.model, va_shift_to  == 0)
end

function constraint_reactive_power_branch_contingency(pm::GenericPowerModel, n::Int, cnd::Int, i::Int, f_idx, t_idx)
    q_fr = PowerModels.var(pm, n, cnd, :q, f_idx)
    q_to = PowerModels.var(pm, n, cnd, :q, t_idx)
    vm_tap_fr = PowerModels.var(pm, n, cnd, :vm_tap, f_idx)
    vm_tap_to = PowerModels.var(pm, n, cnd, :vm_tap, t_idx)

    @constraint(pm.model, q_fr == 0)
    @constraint(pm.model, q_to == 0)
    @constraint(pm.model, vm_tap_fr  == 1)
    @constraint(pm.model, vm_tap_to  == 1)
end
# function constraint_active_power_branch_contingency(pm::GenericPowerModel, n::Int, i::Int, f_idx, t_idx)
#     p_fr = PowerModels.var(pm, n, :p)[f_idx]
#     p_to = PowerModels.var(pm, n, :p)[t_idx]
#     va_shift_fr = PowerModels.var(pm, n, :va_shift)[f_idx]
#     va_shift_to = PowerModels.var(pm, n, :va_shift)[t_idx]
#
#     @constraint(pm.model, p_fr == 0)
#     @constraint(pm.model, p_to == 0)
#     @constraint(pm.model, va_shift_fr  == 0)
#     @constraint(pm.model, va_shift_to  == 0)
# end
#
# function constraint_reactive_power_branch_contingency(pm::GenericPowerModel, n::Int, i::Int, f_idx, t_idx)
#     q_fr = PowerModels.var(pm, n, :q)[f_idx]
#     q_to = PowerModels.var(pm, n, :q)[t_idx]
#     vm_tap_fr = PowerModels.var(pm, n, :vm_tap)[f_idx]
#     vm_tap_to = PowerModels.var(pm, n, :vm_tap)[t_idx]
#
#     @constraint(pm.model, q_fr == 0)
#     @constraint(pm.model, q_to == 0)
#     @constraint(pm.model, vm_tap_fr  == 1)
#     @constraint(pm.model, vm_tap_to  == 1)
# end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
xz  = isapple ?
    "/usr/local/bin/xz" :
    "/usr/bin/xz"

gunzip = isapple ?
    "/usr/bin/gunzip" :
    "/usr/bin/gunzip"


text:
module NURBStoolbox_mixedtest01

using Test;
using MagMechFEM_Matlab2Julia.NURBStoolbox;

@testset "Mixed Tests: 01" begin

  coefs = cat([0.0 0.0; 0.0 1.0],[1.0 1.0; 0.0 1.0],dims=3);

  # TODO: this here should fail as the second knot vector is not valid!
  #       should throw argument error
  @test_skip nrbmak(coefs,[vec([0.0 0.0 1.0 1.]),vec([0.0 0.5 1.0 1.0])]);
  # TODO: this here should fail as the second knot vector is not valid!
  #       should throw argument error
  @test_skip nrbmak(coefs,[vec([0.0 0.0 1.0 1.]),vec([0.0 1.0 1.0])]);
  # TODO: this here should fail as the second knot vector is not valid!
  #       should throw argument error
  @test_skip nrbmak(coefs,[vec([0.0 0.0 1.0 1.]),vec([0.0])]);

  # TODO: this here should fail as the knot vector is not valid!
  #       should throw argument error
  @test_skip nrbmak([0.0 1.5; 0.0 3.0],vec([0.0]));

end # testset

end # module


text:
######################################################################################################
# Copy from https://github.com/xuedong/LinBAI.jl  by Xuedong Shang                                   #
######################################################################################################
using Random;
using CPUTime;
using LinearAlgebra;
include("../utilities/peps.jl");
include("../utilities/expfam.jl");
include("samplingrules.jl");

# Run the learning algorithm, paramterised by a sampling rule
# The stopping and recommendation rules are common
#
# βs must be a list of thresholds *in increasing order*

function play!(i, k, rng, pep, µ, S, N, P, Vinv)
    arm = pep.arms[k]
    Y = sample(rng, getexpfam(pep, 1), arm'µ)
    S .+= Y .* arm
    Vinv .= sherman_morrison(Vinv, arm)
    N[k] += 1
end

function runit(seed, sr, μs, pep::LinearThreshold, βs)
    # seed: random seed. UInt.
    # sr: sampling rule.
    # µs: mean vector.
    # pep: pure exploration problem.
    # βs: list of thresholds.
    convex_sr = (typeof(sr) == ConvexGame) || (typeof(sr) == LearnerK);  # test if P is needed.
    qbc_sr = typeof(sr) == FWSampling;
    βs = collect(βs) # mutable copy
    rng = MersenneTwister(seed)
    K = narms(pep, µs)
    nb_I = nanswers(pep, µs)
    dim = length(μs)
    convex_sr ? P = ones(Int64, (nb_I, K)) : P = ones(Int64, (1, 1))  # counts detailed by answer
    N = zeros(Int64, K)              # counts
    S = zeros(dim)                     # sum of samples
    Vinv = Matrix{Float64}(I, dim, dim)  # inverse of the design matrix

    if qbc_sr # initialize Vxinv matrix
        Vxinv = zeros(Float64, dim,dim);
        for k in 1:K
            Vxinv += pep.arms[k]*(pep.arms[k]')/K;
        end
    end

    baseline = CPUtime_us()
    for k = 1:K
        play!(1, k, rng, pep, µs, S, N, P, Vinv)
    end

    # start sampling rules
    if qbc_sr
        state = start(sr, N, Vxinv)
    elseif convex_sr
        state = start(sr, N, P);
    else
        state = start(sr, N);
    end

    R = Tuple{Array{Int64,1},Array{Int64,1},UInt64}[] # collect return values
    while true
        t = sum(N)
        hµ = Vinv * S  # emp. estimates
        # test stopping criterion
        Z, answer, ξ = glrt(pep, N, hμ, Vinv);

        while Z > βs[1](t)
            popfirst!(βs)
            push!(R, (copy(answer), copy(N), CPUtime_us() - baseline));
            if isempty(βs)
                return R
            end
        end
        # invoke sampling rule
        if convex_sr
            i, k = nextsample(state, pep, argmax([hμ'pep.arms[k] for k=1:K]), ξ, N, P, S, Vinv);
        else
            i, k = nextsample(state, pep, N, S, Vinv)
        end
        play!(i, k, rng, pep, µ, S, N, P, Vinv)
        convex_sr ? P[i, k] += 1 : nothing
        t += 1
    end
end


text:
a = b = 0

n = parse(Int, readline())
for i in 1:n
    sa, sb = split(readline())
    if sa < sb
        global b += 3
    elseif sa > sb
        global a += 3
    else
        a += 1
        b += 1
    end
end
println(a, " ", b)


text:
function bisection(f′, a, b, ϵ)
    if a > b; a, b = b, a; end
    ya, yb = f′(a), f′(b)
    if ya == 0; b = a; end
    if yb == 0; a = b; end
    while b-a > ϵ
        x = (a+b)/2;
        y = f′(x);
        if y == 0
            a, b = x, x;
        elseif sign(y) == sign(ya)
            a = x;
        else
            b = x;
        end
    end
    return (a, b);
end


text:
# From AnticNumberField into something

mutable struct MapDataFromAnticNumberField{T}
  prim_image::T
  isid::Bool

  function MapDataFromAnticNumberField{T}(x::Bool) where T
    @assert x
    z = new{T}()
    z.isid = true
    return z
  end

  function MapDataFromAnticNumberField{T}(x::T) where T
    z = new{T}(x, false)
    return z
  end
end

function image(f::MapDataFromAnticNumberField, y::nf_elem)
  f.isid && return y
  z = parent(defining_polynomial(parent(y)))(y)
  return evaluate(z, f.prim_image)
end

function (f::MapDataFromAnticNumberField)(y::nf_elem)
  return image(f, y)
end

map_data_type(K::AnticNumberField, L) = MapDataFromAnticNumberField{elem_type(L)}

map_data(K::AnticNumberField, L, y) = MapDataFromAnticNumberField{typeof(y)}(y)

map_data(K::AnticNumberField, ::Bool) = MapDataFromAnticNumberField{nf_elem}(true)

map_data(K::AnticNumberField, L, ::Bool) = MapDataFromAnticNumberField{elem_type(L)}(true)

# From NfRel into something
mutable struct MapDataFromNfRel{T, S}
  prim_image::T
  base_field_map_data::S
  isid::Bool

  function MapDataFromNfRel{T, S}(x::T, y::S) where {T, S}
    z = new{T, S}(x, y, false)
    return z
  end

  function MapDataFromNfRel{T, S}(x::Bool) where {T, S}
    @assert x
    z = new{T, S}()
    z.isid = true
    return z
  end
end

function image(f::MapDataFromNfRel, y)
  f.isid && return y
  z = change_base_ring(y.data, f.base_field_map_data)
  return evaluate(z, f.prim_image)
end

function (f::MapDataFromNfRel)(y)
  return image(f, y)
end

map_data_type(K::NfRel, L) = MapDataFromNfRel{elem_type(L), map_data_type(base_field(K), L)}

map_data_type(K::NfRel, L, base_field_map_data::S) where {S} = MapDataFromNfRel{elem_type(L), S}

map_data(K::NfRel, L, y, g) = MapDataFromNfRel{typeof(y), typeof(g.image_data)}(y, g.image_data)

map_data(K::NfRel, L, y) = MapDataFromNfRel{typeof(y), map_data_type(base_field(K), L)}(y, map_data_type(base_field(K), L)(true))

map_data(K::NfRel, L, ::Bool) = MapDataFromNfRel{elem_type(L), map_data_type(base_field(K), L)}(true)

function map_data(K::NfRel, y)
  z = MapDataFromNfRel{typeof(y), map_data_type(base_field(K), parent(y))}
  return z(y, map_data(base_field(K), K, true))
end

# From NfAbsNS into something
mutable struct MapDataFromNfAbsNS{T}
  images::T
  isid::Bool

  function MapDataFromNfAbsNS{T}(x::T) where {T}
    z = new{T}(x, false)
    return z
  end

  function MapDataFromNfAbsNS{T}(x::Bool) where {T}
    @assert x
    z = new{T}()
    z.isid = true
    return z
  end
end

function image(f::MapDataFromNfAbsNS, y)
  f.isid && return y
  return evaluate(y.data, f.images)
end

function (f::MapDataFromNfAbsNS)(y)
  return image(f, y)
end

map_data_type(K::NfAbsNS, L) = MapDataFromNfAbsNS{Vector{elem_type(L)}}

map_data(K::NfAbsNS, L, y) = MapDataFromNfAbsNS{typeof(y)}(y)

map_data(K::NfAbsNS, L, ::Bool) = MapDataFromNfAbsNS{Vector{elem_type(L)}}(true)

# From NfRelNS into something
mutable struct MapDataFromNfRelNS{T, S}
  images::T
  base_field_map_data::S
  isid::Bool

  function MapDataFromNfRelNS{T, S}(x::T, y::S) where {T, S}
    z = new{T, S}(x, y, false)
    return z
  end

  function MapDataFromNfRelNS{T, S}(x::Bool) where {T, S}
    @assert x
    z = new{T, S}()
    z.isid = true
    return z
  end
end

function image(f::MapDataFromNfRelNS, y)
  z = change_base_ring(y.data, f.base_field_map_data)
  return evaluate(z, f.images)
end

function (f::MapDataFromNfRelNS)(y)
  return image(f, y)
end

map_data_type(K::NfRelNS, L) = MapDataFromNfRelNS{Vector{elem_type(L)}, map_data_type(base_field(K), L)}

map_data_type(K::NfRelNS, L, base_field_map_data::S) where {S} = MapDataFromNfRelNS{Vector{elem_type(L)}, S}

# Here g is a NumFieldMor
map_data(K::NfRelNS, L, y, g) = MapDataFromNfRelNS{typeof(y), typeof(g.image_data)}(y, g.image_data)

map_data(K::NfRelNS, L, y) = MapDataFromNfRelNS{typeof(y), map_data_type(base_field(K), L)}(y, map_data_type(base_field(K), L)(true))

#map_data(K::NfRelNS, ::Bool) = MapDataFromNfRelNS{Vector{elem_type(K)}, map_data_type(base_ring(K), K)}(true)

map_data(K::NfRelNS, L, ::Bool) = MapDataFromNfRelNS{Vector{elem_type(L)}, map_data_type(base_ring(K), L)}(true)

mutable struct NumFieldMor{S, T, U, V}
  header::MapHeader{S, T}
  image_data::U
  preimage_data::V

  function NumFieldMor{S, T, U, V}() where {S, T, U, V}
    z = new{S, T, U, V}()
    return z
  end

  NumFieldMor{S, T, U, V}(h::MapHeader{S, T}, i::U, p::V) where {S, T, U, V} = new{S, T, U, V}(f, i, p)
end

function _hom(K, L, x...; preimage = nothing)
  header = MapHeader(K, L)
  image_data = map_data(K, L, x...)
  if preimage !== nothing
    preimage_data = map_data(L, preimage...)
    return NumFieldMor(header, image_data, preimage_data)
  end
  z = NumFieldMor{typeof(K), typeof(L), typeof(image_data), map_data_type(L, K)}()
  z.header = header
  z.image_data = image_data
  return z
end

function image(f::NumFieldMor, x)
  return image(f.image_data, x)
end

function preimage(f::NumFieldMor, x)
  return image(f.preimage_data, x)
end

function (f::NumFieldMor)(x)
  return image(f, x)
end

function image_primitive_element(f::NumFieldMor{AnticNumberField})
  return f.image_data.prim_img
end

function image_primitive_element(f::NumFieldMor{<:NfRel})
  return f.image_data.prim_img
end


text:
@testset "5.1 Inverse sine" begin
    include("5.1.2 (d x)^m (a+b arcsin(c x))^n.jl")
    include("5.1.4 (f x)^m (d+e x^2)^p (a+b arcsin(c x))^n.jl")
    include("5.1.5 Inverse sine functions.jl")
end

max_stars_repo_path: test/5 Inverse trig functions/5.1 Inverse sine/.jl


text:
using FSM

cn = Constants{Float64}()

ebm = EBM{Float64}(
    am=1,
    cm=1,
    dm=1,
    em=1,
    hm=1,
)

Ta = 292.3
Ua = 5.0
Ps = 87360.0
SW = 100.0
LW = 100.0
RH = 100.0

Qs = qsat(true, Ps, Ta, cn)
Qa = (RH / 100.0) * Qs

surf_props(ebm, cn, Sf)
surf_exch(ebm, cn, Ta, Ua)

surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW)

@code_warntype surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW)

time_surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW) = @time surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW)

time_surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW)
time_surf_ebal(ebm, cn, Ta, Qa, Ua, Ps, SW, LW)


text:
abstract type Item end
id(i::Item)=String(typeof(i).name.name)
text(i::Item)=langs[:name][typeof(i).name.name]
stack(::Item)=0x10 # 16
c_use(::Item)=false
one_use(::Item)=true
t_blk(::Item)=missing
function i_show(i::Item,con=dcon)
	clear_rect(con,0,0;color=RGB_(0x70,0x80,0x90))
	i_show_l(i,con::DrawContext)
end
function i_show_l(i::Item,con::DrawContext)
	iid=id(i)
	so=get_loadedimg(iid)
	fill_image(con,so,0,0)
end
function destroy_with(::Item,b::Block)
	har=hard(b)
	stre=strength(ply)
	if har>stre+15 return -1 end
	temp=har-stre+20
	return temp<0 ? 0 : temp<<3
end

struct EI<:Item end
id(i::EI)=""
text(i::EI)=""
i_show_l(::EI,::Any)=nothing

macro iwrap(ty,sy::Symbol)
	return :($sy(i::$ty)=$sy(i.wr))
end
macro iwraps(ty,sys::Symbol...)
	for sy in sys
		eval(:($sy(i::$ty)=$sy(i.wr)))
	end
end

mutable struct IFB<:Item
	wr::Block
end
@iwraps IFB id text stack
t_blk(i::IFB)=deepcopy(i.wr)
i_show(i::IFB,con=dcon)=i_show(i.wr,con::DrawContext)

mutable struct NI<:Item
	wr::Item
	name::String
end
id(::NI)="NI"
text(i::NI)=i.name
stack(::NI)=0x1
@iwraps NI c_use one_use use e_use updg t_blk
i_show(i::NI,con::DrawContext)=i_show(i.wr,con::DrawContext)

mutable struct DI<:Item
	wr::Item
	detail::Dict{Symbol,Any}
end
id(::DI)="DI"
text(i::DI)=haskey(i.detail,:name) ? i.detail[:name] : text(i.wr)
stack(::DI)=0x1
@iwraps DI c_use one_use use e_use updg t_blk
i_show(i::DI,con::DrawContext)=i_show(i.wr,con::DrawContext)

include("模具.jl")
include("书.jl")
include("桶.jl")
abstract type 零件<:Item end
id(i::零件)="$(i.mat)$(typeof(i).name.name)"
text(i::零件)=langs[:name][typeof(i.mat).name.name]*langs[:name][typeof(i).name.name]
include("薄板.jl")
mutable struct 锭<:零件 mat::Material end
mutable struct 棍<:零件 mat::Material end


text:
#  #nv-nsight-cu-cli --mode=launch julia
# #  using  Test, Revise ,CUDA

# #  includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\aPrfofiling\\profilingProcessMaskData.jl")

# #  CUDA.@profile wrapForProfile()

# using  Test, Revise,CUDA
# includet("C:\\GitHub\\GitHub\\NuclearMedEval\\test\\includeAllUseFullForTest.jl")

# using ..CUDAGpuUtils ,..IterationUtils,..ReductionUtils , ..MemoryUtils,..CUDAAtomicUtils, ..BasicStructs
# using Shuffle,..ResultListUtils, ..MetadataAnalyzePass,..MetaDataUtils,..WorkQueueUtils,..ProcessMainDataVerB,..HFUtils, ..ScanForDuplicates
# using ..MainOverlap, ..RandIndex , ..ProbabilisticMetrics , ..VolumeMetric ,..InformationTheorhetic
# using ..CUDAAtomicUtils, ..TpfpfnKernel, ..InterClassCorrKernel,..MeansMahalinobis
# using HDF5

# pathToHd5= "C:\\Users\\1\\PycharmProjects\\pythonProject3\\mytestfile.hdf5"

# const g = h5open(pathToHd5, "r+")
# not_translated=read(g["not_translated"])
# translated=read(g["translated"])

# arrGold = CuArray(not_translated)
# arrAlgo = CuArray(translated)

# sizz= size(not_translated)
# conf = ConfigurtationStruct(dice=true)
# numberToLooFor = UInt8(1)
# args,threads,blocks,metricsTuplGlobal= TpfpfnKernel.prepareForconfusionTableMetrics(arrGold    , arrAlgo    ,numberToLooFor  ,conf)

# argsB = TpfpfnKernel.getTpfpfnData!(arrGold ,arrAlgo   ,args,threads,blocks,metricsTuplGlobal)
# dice = metricsTuplGlobal[4][1]

# function getBlockTpFpFnProf(goldGPU#goldBoolGPU
#     , segmGPU#segmBoolGPU
#     #,sliceMetricsTupl
#     ,tp,tn,fp,fn#tp,tn,fp,fn
#     ,arrDims,totalNumbOfVoxels,iterLoop,pixPerSlice
#     ,numberToLooFor#numberToLooFor
#     #,metricsTuplGlobal
#     ,conf)

#     shmemSum = @cuStaticSharedMem(UInt32, (33,3))
#     shmemblockData = @cuStaticSharedMem(UInt32,(32, 32 ,3))

#     # locArr= zeros(MVector{3,UInt16})
#     @iterateLinearlyMultipleBlocks(iterLoop,pixPerSlice,totalNumbOfVoxels,
#     #inner expression
#     begin
#         # CUDA.@cuprint "i $(i)  val $(goldGPU[i])"
#         #updating variables needed to calculate means
#         boolGold = goldGPU[i]==numberToLooFor
#         boolSegm = segmGPU[i]==numberToLooFor
#         @inbounds  shmemblockData[threadIdxX(),threadIdxY(), (boolGold & boolSegm + boolSegm +1)]+=(boolGold | boolSegm)
#         #    locArr[1]=boolGold
#         #   @inbounds locArr[ (boolGold & boolSegm + boolSegm +1) ]+=(boolGold | boolSegm)
#     end)

#     # tell what variables are to be reduced and by what operation
#     @redWitAct(offsetIter,shmemSum, shmemblockData[threadIdxX(),threadIdxY(),1],+,    shmemblockData[threadIdxX(),threadIdxY(),2],+,     shmemblockData[threadIdxX(),threadIdxY(),3],+   )
#     sync_threads()
#     @addAtomic(shmemSum, fn, fp,tp)

# return
# end

# function getTpfpfnDataProf(goldGPU
#     , segmGPU
#     ,args,threads,blocks,metricsTuplGlobal)

#     for i in  1:4
#         CUDA.fill!(args[i],0)
#     end

#     for i in 1:length(metricsTuplGlobal)
#         metricsTuplGlobal[i]=0
#         #CUDA.fill!(args[11][i],0)
#     end

#     #get tp,fp,fna and slicewise results if required
#     @cuda threads=threads blocks=blocks getBlockTpFpFnProf(goldGPU, segmGPU,args...) #args[8][3]  is number of slices ...
#     #TpfpfnKernel.getMetricsCPU(args[1][1],args[3][1], args[4][1],(args[6]-(args[1][1] +args[3][1]+ args[4][1] )) ,metricsTuplGlobal,args[10],1 )
#     return args
# end#getTpfpfnData

# function wrapForProfile()

#     getTpfpfnDataProf(arrGold ,arrAlgo   ,args,threads,blocks,metricsTuplGlobal)
#     return metricsTuplGlobal[4][1]
# end

# # testArrIn = CUDA.ones(Bool,32,32,32);
# # testArrOut = CUDA.zeros(Bool,34,34,34);

# # function testKernprocessMaskDataB(testArrInn,testArrOut)
# #     resShmem =  @cuStaticSharedMem(Bool,(34,34,34))#+2 in order to get the one padding
# #     clearMainShmem(resShmem)
# #     blockBeginingX=UInt8(0)
# #     blockBeginingY=UInt8(0)
# #     blockBeginingZ=UInt8(0)
# #     isMaskFull= false
# #     isMaskEmpty= false
# #     #here we will store in registers data uploaded from mask for later verification wheather we should send it or not
# #     locArr= UInt32(0)

# #     @unroll for zIter in UInt8(1):UInt8(32)# most outer loop is responsible for z dimension
# #         processMaskData( testArrInn[(blockBeginingX+threadIdxX()),(blockBeginingY +threadIdxY()),(blockBeginingZ+zIter)], zIter, resShmem,locArr)
# #     end
# #     sync_threads()
# #         #fillGlobalFromShmem(testArrOut,resShmem)

# #     return
# # end

# # @cuda threads=(32,32) blocks=1 testKernprocessMaskDataB(testArrIn,testArrOut)

# # numb = UInt32(0)
# # #settingcorrectly
# # numb |= true << UInt8(2)
# # numb

# # numb = UInt32(0)
# # #settingcorrectly
# # numb |= 1 << 1
# # numb |= 1 << 2
# # numb |= 1 << 0
# # numb |= 1 << 0
# # numb |= 1 << 5
# # numb

# #reading...
# # numb>>1 & UInt32(1)
# # numb>>2 & UInt32(1)
# # numb>>3 & UInt32(1)
# # numb>>4 & UInt32(1)
# # numb>>5 & UInt32(1)

# # #processMaskDataB( testArrInn[threadIdxY(),1,1], zIter, resShmem ) # coalesced


text:
module PanelCutting

using Revise
using Printf
using Markdown
using Unitful
using UnitfulUS
# using UnitfulCurrency     #Problems
using Match
using DataStructures
using InteractiveUtils
using NativeSVG
using DisplayAs
using UUIDs
using Plots
using Logging
using NahaGraphs

# Work around issues with UnitfulCurrency
function money(amount)
    try
        amount * uparse("USD")
    catch e
        amount
    end
end

MoneyType = try
    Quantity{N, CURRENCY} where N
catch e
    Real
end

if MoneyType == Real
    unmoney(amount::Real) = amount
else
    eval(:(unmoney(amount::Quantity) = ustrip(Real, u"USD", amount)))
end

export money, MoneyType, unmoney

include("AllOf.jl")
include("axes.jl")
include("panels.jl")
include("cutting.jl")
include("supplier.jl")
include("search.jl")
include("report.jl")
include("graph.jl")
include("dot_styles.jl")
include("svg.jl")

#  include("./panels_notebook.jl")

end


text:
typealias Chars Union(Char,AbstractVector{Char},Set{Char})

# search in a string from the reverse direction
function _rsearch(s::String, c::Chars, i::Integer)
    l = endof(s)
    ret = search(reverse(s), c, l-i+1)
    (0 == ret) ? 0 : (l-ret+1)
end
_rsearch(s::String, c::Chars) = _rsearch(s,c,endof(s))

function _rsearch(s::String, t::String, i::Integer)
    l = endof(s)
    r = search(reverse(s), reverse(t), l-i+1)
    (0 == r.start) ? r : Range1(l-r.start-r.len+2, r.len)
end
_rsearch(s::String, t::String) = _rsearch(s,t,endof(s))

function _rsplit(str::String, splitter, limit::Integer, keep_empty::Bool)
    strs = String[]
    i = start(str)
    n = endof(str)
    r = _rsearch(str,splitter,n)
    j = first(r)-1
    k = last(r)
    while((0 <= j < n) && (length(strs) != limit-1))
        if(i <= k)
            (keep_empty || (k < n)) && unshift!(strs, str[k+1:n])
            n = j
        end
        (k <= j) && (j = prevind(str,j))
        r = _rsearch(str,splitter,j)
        j = first(r)-1
        k = last(r)
    end
    (keep_empty || (n > 0)) && unshift!(strs, str[1:n])
    return strs
end
_rsplit(s::String, spl, n::Integer) = _rsplit(s, spl, n, true)
_rsplit(s::String, spl, keep::Bool) = _rsplit(s, spl, 0, keep)
_rsplit(s::String, spl)             = _rsplit(s, spl, 0, true)

const _default_delims = [' ','\t','\n','\v','\f','\r']
_rsplit(str::String) = _rsplit(str, _default_delims, 0, false)


text:
using PkgCleanup
using Test

@testset "PkgCleanup.jl" begin
    @test isdir(PkgCleanup.logs_directory())
end


text:
using Test, ReferenceTests, BSON
using ParallelStencil
import ParallelStencil: @reset_parallel_stencil

ENV["USE_GPU"] = false
ENV["DO_VIZ"]  = false
ENV["DO_SAVE"] = false
ENV["DO_SAVE_VIZ"] = false
ENV["NX"] = 64
ENV["NY"] = 64
ENV["USE_RETURN"] = true

# Reference test using ReferenceTests.jl
comp(d1, d2) = keys(d1) == keys(d2) && all([ isapprox(v1, v2) for (v1,v2) in zip(values(d1), values(d2))])

## 2D tests
include("../scripts/diff_2D_lin.jl")
@reset_parallel_stencil()
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(yc), 12)))
d2d1  = Dict(:X=> xc[indsx], :H=>H[indsx,indsy])

include("../scripts/diff_2D_linstep.jl")
@reset_parallel_stencil()
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(yc), 12)))
d2d2  = Dict(:X=> xc[indsx], :H=>H[indsx,indsy])

include("../scripts/diff_2D_nonlin.jl")
@reset_parallel_stencil()
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(yc), 12)))
d2d3  = Dict(:X=> xc[indsx], :H=>H[indsx,indsy])

include("../scripts/diff_2D_nonlin_perf.jl")
@reset_parallel_stencil()
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(yc), 12)))
d2d4  = Dict(:X=> xc[indsx], :H=>H[indsx,indsy])

include("../scripts/diff_2D_react.jl")
@reset_parallel_stencil()
indsx = Int.(ceil.(LinRange(1, length(xc), 12)))
indsy = Int.(ceil.(LinRange(1, length(yc), 12)))
d2d5  = Dict(:X=> xc[indsx], :H=>H[indsx,indsy])

@testset "Reference-tests diffusion 2D" begin
    @test_reference "reftest-files/test_diff_2D_lin.bson" d2d1 by=comp
    @test_reference "reftest-files/test_diff_2D_linstep.bson" d2d2 by=comp
    @test_reference "reftest-files/test_diff_2D_nonlin.bson" d2d3 by=comp
    @test_reference "reftest-files/test_diff_2D_nonlin_perf.bson" d2d4 by=comp
    @test_reference "reftest-files/test_diff_2D_react.bson" d2d5 by=comp
end

@reset_parallel_stencil()


text:
# # using Revise
using Test
using Logging
using LarSurf
using SparseArrays
# using LinearAlgebraicRepresentation
# Lar = LinearAlgebraicRepresentation
# Logging.configure(level==Logging.Debug)

# include("../src/LarSurf.jl")
# include("../src/block.jl")

@testset "array compare roll independent" begin
    @test LarSurf.array_equal_roll_invariant([1,2,3,4,5], [2,3,4,5,1])

end


text:
# https://mambajl.readthedocs.io/en/latest/tutorial.html#mcmc-simulation

using Mamba, LinearAlgebra, Distributions

model = Model(y = Stochastic(1,(μ , ν) ->  MvNormal(μ , sqrt(ν)), false),
              μ  = Logical(1,(xm, β) -> xm * β, false),
              β = Stochastic(1, () -> MvNormal(2, sqrt(1000))),
              ν = Stochastic(() -> InverseGamma(0.001, 0.001)) );

## Hybrid No-U-Turn and Slice Sampling Scheme
sc1 = [NUTS(:β), Slice(:ν, 3.0)];

## No-U-Turn Sampling Scheme
sc2 = [NUTS([:β, :ν])];

## User-Defined Samplers

Gibbs_β = Sampler([:β],
   (β, ν, xm, y) ->  begin
      β_mean = mean(β.distr)
      β_invcov = invcov(β.distr)
      σ = inv(Symmetric(xm' * xm / ν + β_invcov))
      μ = σ * (xm' * y / ν  +  β_invcov * β_mean)
      rand(MvNormal(μ , σ))
   end
);

Gibbs_ν = Sampler([:ν], (μ, ν, y) -> begin
      a = length(y) / 2.0 + shape(ν.distr)
      b = sum(abs2, y - μ) / 2.0 + scale(ν.distr)
      rand(InverseGamma(a, b))
    end
);

## User-Defined Sampling Scheme
sc3 = [Gibbs_β, Gibbs_ν];

setsamplers!(model, sc3);

draw(model, filename="lineDAG.dot");

run(`ls -l lineDAG.dot`)     # Recall wildcards do NOT work, they are a function of the shell

## Data
line = Dict{Symbol, Any}(:x => [1, 2, 3, 4, 5], :y => [1, 3, 3, 3, 5] )
line[:xm] = [ones(5) line[:x]];

## Initial Values
inits = [ Dict{Symbol, Any}(
    :y => line[:y], :β => rand(Normal(0, 1), 2),:ν => rand(Gamma(1, 1))) for i in 1:3 ];

setsamplers!(model, sc1)
sim1 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

gelmandiag(sim1, mpsrf=true, transform=true) |> showall

setsamplers!(model, sc2)
sim2 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

gelmandiag(sim2, mpsrf=true, transform=true) |> showall

setsamplers!(model, sc3)
sim3 = mcmc(model, line, inits, 10000, burnin=250, thin=2, chains=3)

gelmandiag(sim3, mpsrf=true, transform=true) |> showall

sim = sim1

gewekediag(sim) |> showall

heideldiag(sim) |> showall

rafterydiag(sim) |> showall

# https://mambajl.readthedocs.io/en/latest/mcmc/chains.html#section-convergence-diagnostics

## Summary Statistics
describe(sim)

## Highest Posterior Density Intervals
hpd(sim)

## Lag-Autocorrelations
autocor(sim)

## State Space Change Rate (per Iteration)
changerate(sim)

## Deviance Information Criterion
dic(sim3)

## Subset Sampler Output
sim = sim[1000:5000, ["β[1]", "β[2]"], :]
describe(sim)

using Gadfly

p = plot(sim)

draw(p, filename="summaryplot.svg")
draw(p)

## Autocorrelation and running mean plots
p = plot(sim, [:autocor, :mean], legend=true)

draw(p, nrow=3, ncol=2, filename="autocormeanplot.svg")
draw(p, nrow=3, ncol=2)

## Pairwise contour plots
p = plot(sim, :contour)

draw(p, nrow=2, ncol=2, filename="contourplot.svg")
draw(p, nrow=2, ncol=2)


text:
module Solution

using Random
using Distributions
using ..Space, ..RNG, ..Tagging, ..Var, ..Traits
using ..Condition: Conditional
import ..Var
export solution

# # Equality Condition
# Many (but not all) inference problems are of the form `X = x` where `X` is a
# a random variable abd `x` is a concrete constant.  Problems in this form often
# permit more tractable inference.  To exploit this we use a new type EqualityCondition
# so that we can identify theses cases just from their types

"""
`solution(x)`

Solution (aka model, interpretation) to constraints on random variable.

Returns any `ω` such that `x(ω) != ⊥`
"""
function solution end

const ConstTypes = Union{Real, Array{<:Real}}
const EqualityCondition{A, B} = PwVar{Tuple{A, B}, typeof(==)} where {A, B <: ConstTypes}
tagcondition(ω, condition) = tag(ω, (condition = condition,))

"""

`solution(rng, f, Ω)`

```
μ = 1 ~ Normal(0, 1)
y = 2 ~ Normal(μ, 1)
μc = μ |ᶜ (y ==ₚ 5.0)
solution(μc)
```
"""
function solution(rng::AbstractRNG,
                  f::Conditional{X, Y},
                  Ω = defΩ()) where {X, Y <: EqualityCondition}
  ω = tagrng(Ω(), rng)
  ωc = tagcondition(ω, f.y)
  f(ωc)
  ω
end

idof(m::Member) = m.id
  idof(v::Variable) = v.f.id

function Var.prehook(::trait(Cond), d::Distribution, id, ω)
  # FIXME: Is this correct?
  matches = idof(ω.tags.condition.a) == id
  if matches
    inv = invert(d, ω.tags.condition.b)
    ω[id] = (primdist(d), inv)
  end
end

solution(f::Conditional, Ω = defΩ()) =
  solution(Random.GLOBAL_RNG, f, Ω)

end

## What's wrong
# We need to indicate that a value does not need to be modified for conditioning
# I'm not 100 percent confident in the way im checking (id matching)
# this g(nothing) seems smelly
# What if we have more than one condition??
# Can we distinguish from X == X to X == Const at type level?
# how are we actually distapplyign to do it !
# what more do we need for MH

# Changes
# change condition such that we also wat condition and add default behaviour
#   to do this erroring
# setup function barrier in this prehook such that if logpdf is there we update
# get rid of seen shit and rely on fact that rand in lazyomega is statefuk
# account for fact that we may have already done the update


text:
#Define functions for the calculation and formation of the covariance matrix S

#Compute the element of the S-matrix and store them into a vector
function S_gen(dim::Int64)

       S_vec = Complex{Float64}[]

       for i=1:dim
           for j=i:dim
      @views   fr = cov(ENWfDer[:,j+1],ENWfDer[:,i+1],corrected=false)
               push!(S_vec,fr)
             end
       end

     return S_vec

end

#Form and return the S-matrix
function Gen_S(S_vec::Array{Complex{Float64},1})

    S_mat = zeros(Complex{Float64},dim,dim)

    k = 1
    for i=1:dim
        for j=i:dim
            @views S_mat[i,j] = S_vec[k]
            @views S_mat[j,i] = conj(S_vec[k])
            k +=1
        end
    end

    return S_mat

end

max_stars_repo_name: GiamFbn/MyULTRAFAST

max_issues_repo_name: GiamFbn/MyULTRAFAST

max_forks_repo_name: GiamFbn/MyULTRAFAST


text:
# normal
solver = MCTSSolver(n_iterations=n_iter, depth=depth, exploration_constant=ec, enable_tree_vis=true)
mdp = LegacyGridWorld()

policy = solve(solver, mdp)

state = GridWorldState(1,1)

a = action(policy, state)

tree = D3Tree(policy, state)

io = IOBuffer()
show(io, MIME("text/plain"), tree)
take!(io)
show(io, MIME("text/html"), tree)

# dpw
solver = DPWSolver(n_iterations=n_iter, depth=depth, exploration_constant=ec, rng=MersenneTwister(13), tree_in_info=true)
mdp = LegacyGridWorld()

policy = solve(solver, mdp)

state = GridWorldState(1,1)

a, info = action_info(policy, state)

tree = D3Tree(policy)
tree = D3Tree(info[:tree])

io = IOBuffer()
show(io, MIME("text/plain"), tree)
take!(io)
show(io, MIME("text/html"), tree)


text:
# Reporting.jl

module Reporting

using Dates

export aging

import ..AppliAR: retrieve_unpaid_invoices, retrieve_paid_invoices

import ..AppliAR: id, name, date, meta, header, body, students, price_per_student, vat_perc

struct Aging
    id_inv::String
    csm::String
    inv_date::Date
    amount::Float64
    days::Day
end

aging(path_file_unpaid, path_file_paid) = begin
    unpaid_invoices = retrieve_unpaid_invoices(;path=path_file_unpaid)
    paid_invoices = retrieve_paid_invoices(;path=path_file_paid)

    paid = [id(x) for x in paid_invoices]
    unpaid = filter(x -> id(x) ∉ paid, unpaid_invoices)

    list = []
    for invoice in unpaid
        id_inv = id(invoice)
        csm = name(header(invoice))
        inv_date = Date(date(meta(invoice)))
        b = body(invoice)
        amount = (price_per_student(b) * length(students(b))) * (1 + vat_perc(b))
        days = Dates.today() - inv_date
        aging_item = Aging(id_inv, csm, inv_date, amount, days)
        push!(list, aging_item)
    end
    return list
end # aging

end # module

max_stars_repo_name: rbontekoe/AppliAR

max_issues_repo_name: rbontekoe/AppliAR

max_forks_repo_name: rbontekoe/AppliAR


text:
push!(LOAD_PATH, pwd())

module TestModule

using A
using B
using Test
using JLD2

x =  AType(1)

S = BType(x)

@testset "save and load" begin
	T = deepcopy(S)
	@test T == S
	fn = joinpath(mktempdir(), "test_out.jld")
	@save fn T
	T = nothing
	@load fn T
	@test T == S
end

end