text:
#NOTE: We do the error handling of nnp_initialize while loading NNPACK
function nnp_initialize()
    ccall((:nnp_initialize, libnnpack), nnp_status, (),)
end

function nnp_deinitialize()
    @check ccall((:nnp_deinitialize, libnnpack), nnp_status, (),)
end

function pthreadpool_create(n = 0)
    ccall((:pthreadpool_create, libnnpack), Ptr{Cvoid}, (Csize_t,), n)
end

function nnp_relu_output(batch_size, channels, input, output, negative_slope, threadpool)
    @check ccall((:nnp_relu_output, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Cfloat, pthreadpool_t), batch_size, channels, input, output, negative_slope, threadpool)
end

function nnp_relu_output(x::AbstractArray{Float32,N}, y::AbstractArray{Float32,N}; negative_slope::AbstractFloat = 0.0, threadpool = shared_threadpool[]) where {N}
    # Investigate why the channel and batch dims need to specified like this
    nnp_relu_output(prod(size(x)[N-1:N]), prod(size(x)[1:N-2]), x, y, negative_slope, threadpool)
    y
end

function nnp_relu_input_gradient(batch_size, channels, grad_output, input, grad_input, negative_slope, threadpool)
    @check ccall((:nnp_relu_input_gradient, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Cfloat, pthreadpool_t), batch_size, channels, grad_output, input, grad_input, negative_slope, threadpool)
end

function nnp_relu_input_gradient(x::AbstractArray{Float32,N}, dy::AbstractArray{Float32,N}, dx::AbstractArray{Float32,N}; negative_slope::AbstractFloat = 0.0, threadpool = shared_threadpool[]) where {N}
    # Investigate why the channel and batch dims need to specified like this
    nnp_relu_input_gradient(Csize_t(prod(size(x)[N-1:N])), prod(size(x)[1:N-2]), dy, x, dx, negative_slope, threadpool)
    dx
end

function nnp_softmax_output(batch_size, channels, input, output, threadpool)
    @check ccall((:nnp_softmax_output, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t), batch_size, channels, input, output, threadpool)
end

function nnp_softmax_output(x::AbstractVecOrMat{Float32}, y::AbstractVecOrMat{Float32}; threadpool = shared_threadpool[])
    nnp_softmax_output(ndims(x) == 2 ? size(x, 2) : 1, size(x, 1), x, y, threadpool)
    y
end

#FIXME: Output of fully connected not consistent with `kernel * input`
#NOTE: This most likely due to nnpack being row major. Investigate this.

function nnp_fully_connected_output(batch_size, input_channels, output_channels, input, kernel, output, threadpool, profile)
    @check ccall((:nnp_fully_connected_output, libnnpack), nnp_status, (Csize_t, Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t, Ptr{Cvoid}), batch_size, input_channels, output_channels, input, kernel, output, threadpool, C_NULL)
end

function nnp_fully_connected_output(x::AbstractArray{Float32,2}, w::AbstractArray{Float32,2}, y::AbstractArray{Float32,2}; profile = nothing, threadpool = shared_threadpool[])
    profile = profile == nothing ? nnp_profile() : profile
    nnp_fully_connected_output(size(x, 2), size(x, 1), size(w, 1), x, w, y, threadpool, profile)
    y
end

function nnp_fully_connected_inference_f16f32(input_channels, output_channels, input, kernel, output, threadpool)
    @check ccall((:nnp_fully_connected_inference_f16f32, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cvoid}, Ptr{Cfloat}, pthreadpool_t), input_channels, output_channels, input, kernel, output, threadpool)
end

nnp_fully_connected_inference_f16f32(x::AbstractVector{Float32}, w::AbstractArray{Float16,2}, y::AbstractVector{Float32}; threadpool = shared_threadpool[]) =
    nnp_fully_connected_inference(reshape(x, size(x), 1), w, reshape(y, size(y), 1), threadpool = threadpool)

function nnp_fully_connected_inference_f16f32(x::AbstractMatrix{Float32}, w::AbstractArray{Float16,2}, y::AbstractMatrix{Float32}; threadpool = shared_threadpool[])
    nnp_fully_connected_inference(size(x, 1), size(y, 1), x, w, y, threadpool)
    y
end

function nnp_fully_connected_inference(input_channels, output_channels, input, kernel, output, threadpool)
    @check ccall((:nnp_fully_connected_inference, libnnpack), nnp_status, (Csize_t, Csize_t, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t), input_channels, output_channels, input, kernel, output, threadpool)
end

nnp_fully_connected_inference(x::AbstractVector{Float32}, w::AbstractArray{Float32,2}; threadpool = shared_threadpool[]) =
    nnp_fully_connected_inference(reshape(x, size(x), 1), w, threadpool = threadpool)

function nnp_fully_connected_inference(x::AbstractMatrix{Float32}, w::AbstractMatrix{Float32}, y::AbstractMatrix{Float32}; threadpool = shared_threadpool[])
    nnp_fully_connected_inference(size(x, 1), size(y, 1), x, w, y, threadpool)
    y
end

function nnp_max_pooling_output(batch_size, channels, input_size, input_padding, pooling_size, pooling_stride, input, output, threadpool)
    @check ccall((:nnp_max_pooling_output, libnnpack), nnp_status, (Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, pthreadpool_t), batch_size, channels, input_size, input_padding, pooling_size, pooling_stride, input, output, threadpool)
end

function nnp_max_pooling_output(x::AbstractArray{Float32,4}, y::AbstractArray{Float32,4}, kernel::Tuple; padding = 0, stride = 1, threadpool = shared_threadpool[])
    input_size = nnp_size(Csize_t.((size(x, 1), size(x, 2)))...)
    pooling_size = nnp_size(Csize_t.(kernel)...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    pooling_stride = nnp_size(Csize_t.(padtuple(x, stride))...)
    nnp_max_pooling_output(size(x, 4), size(x, 3), input_size, input_padding, pooling_size, pooling_stride, x, y, threadpool)
    y
end

#TODO: Add wrapper for convolution inference

function nnp_convolution_input_gradient(algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, grad_output, kernel, grad_input, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, profile)
    @check ccall((:nnp_convolution_input_gradient, libnnpack), nnp_status, (nnp_convolution_algorithm, Csize_t, Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cvoid}, Csize_t, nnp_activation, Ptr{Cvoid}, pthreadpool_t, Ptr{Cvoid}), algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, grad_output, kernel, grad_input, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, C_NULL)
end

function nnp_convolution_input_gradient(dx::AbstractArray{Float32,4}, x::AbstractArray{Float32,4}, dy::AbstractArray{Float32,4}, w::AbstractArray{Float32,4}; algo::nnp_convolution_algorithm = UInt32(0), workspace_buffer = nothing, workspace_size = 0, padding = 0, stride = 1, threadpool = shared_threadpool[], profile = nothing)
    input_size = nnp_size(Csize_t.((size(x,1), size(x,2)))...)
    kernel_size = nnp_size(Csize_t.((size(w,1),size(w,2)))...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    profile = profile == nothing ? nnp_profile() : profile
    workspace_buffer = workspace_buffer === nothing ? C_NULL : workspace_buffer
    nnp_convolution_input_gradient(UInt32(algo), size(x,4), size(x,3), size(w,4), input_size, input_padding, kernel_size, dy, w, dx, workspace_buffer, workspace_size, UInt32(0), C_NULL, threadpool, profile)
    dx
end

function nnp_convolution_kernel_gradient(algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, grad_output, grad_kernel, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, profile)
    @check ccall((:nnp_convolution_kernel_gradient, libnnpack), nnp_status, (nnp_convolution_algorithm, Csize_t, Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cvoid}, Csize_t, nnp_activation, Ptr{Cvoid}, pthreadpool_t, Ptr{Cvoid}), algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, grad_output, grad_kernel, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, C_NULL)
end

function nnp_convolution_kernel_gradient(dw::AbstractArray{Float32,4}, x::AbstractArray{Float32,4}, dy::AbstractArray{Float32,4}, w::AbstractArray{Float32,4}; algo::nnp_convolution_algorithm = UInt32(0), workspace_buffer = nothing, workspace_size = 0, padding = 0, stride = 1, threadpool = shared_threadpool[], profile = nothing)
    input_size = nnp_size(Csize_t.((size(x,1), size(x,2)))...)
    kernel_size = nnp_size(Csize_t.((size(w,1),size(w,2)))...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    profile = profile == nothing ? nnp_profile() : profile
    workspace_buffer = workspace_buffer === nothing ? C_NULL : workspace_buffer
    nnp_convolution_kernel_gradient(UInt32(algo), size(x,4), size(x,3), size(w,4), input_size, input_padding, kernel_size, x, dy, dw, workspace_buffer, workspace_size, UInt32(0), C_NULL, threadpool, profile)
    dw
end

function nnp_convolution_output(algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, kernel, bias, output, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, profile)
    @check ccall((:nnp_convolution_output, libnnpack), nnp_status, (nnp_convolution_algorithm, Csize_t, Csize_t, Csize_t, nnp_size, nnp_padding, nnp_size, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cfloat}, Ptr{Cvoid}, Csize_t, nnp_activation, Ptr{Cvoid}, pthreadpool_t, Ptr{Cvoid}), algorithm, batch_size, input_channels, output_channels, input_size, input_padding, kernel_size, input, kernel, bias, output, workspace_buffer, workspace_size, activation, activation_parameters, threadpool, C_NULL)
end

function nnp_convolution_output(y::AbstractArray{Float32,4}, x::AbstractArray{Float32,4}, w::AbstractArray{Float32,4}, b::AbstractArray{Float32,1}; algo::nnp_convolution_algorithm = UInt32(0), workspace_buffer = nothing, workspace_size = 0, padding = 0, stride = 1, threadpool = shared_threadpool[], profile = nothing)
    input_size = nnp_size(Csize_t.((size(x,1), size(x,2)))...)
    kernel_size = nnp_size(Csize_t.((size(w,1),size(w,2)))...)
    input_padding = nnp_padding(Csize_t(padding[2]), Csize_t(padding[1]), Csize_t(padding[2]), Csize_t(padding[1]))
    profile = profile == nothing ? nnp_profile() : profile
    workspace_buffer = workspace_buffer === nothing ? C_NULL : workspace_buffer
    nnp_convolution_output(UInt32(algo), size(x,4), size(x,3), size(w,4), input_size, input_padding, kernel_size, x, w, b, y, workspace_buffer, workspace_size, UInt32(0), C_NULL, threadpool, profile)
    y
end


text:
# There are exactly ten ways of selecting three from five, 12345:
#
# 123, 124, 125, 134, 135, 145, 234, 235, 245, and 345
#
# In combinatorics, we use the notation, (5 3)=10.
#
# In general, (n r)=n!/(r!(n−r)!), where r≤n, n!=n×(n−1)×...×3×2×1, and 0!=10!=1.
#
# It is not until n=23 that a value exceeds one-million: (23 10)=1144066.
#
# How many, not necessarily distinct, values of (n r) for 1≤n≤100, are greater
# than one-million?

using ProjectEulerSolutions

# We note that (n r) = (n (n-r)), and that when r = n/2, the value of (n r) is
# maximal.  Thus we calculate from (n r/2-k) for k until we get a value of more
# than 1 million.  We save more time by further constraining the range of k
# values.  This is mostly a book-keeping exercise.
function p053solution(n_max::Integer=3, thresh=1_000)::Integer
    n_min = 23
    r = 3

    total = 0
    for n in n_max:-1:n_min
        # The value of r where we surpass the threshold will only get larger
        while binomial(n, r) < thresh
            r += 1
        end
        total += (n - r * 2  + 1)
    end
    return total
end

p053 = Problems.Problem(p053solution)

Problems.benchmark(p053, 100, 1_000_000)

ext: jl

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
module GMSHreader
#=
Read GMSH files:

Supported MeshFormat 4, ASCII

Supported File Groups:
    - MeshFomrat
    - PhysicalNames
    - Entities
    - Nodes
    - Elements

GMSH Element Types:
1:  2-node line.
2:  3-node triangle.
3:  4-node quadrangle.
4:  4-node tetrahedron.
5:  8-node hexahedron.
6:  6-node prism.
7:  5-node pyramid.
8:  3-node second order line (2 nodes associated with the vertices and 1
    with the edge).
9:  6-node second order triangle (3 nodes associated with the vertices and 3
    with the edges).
10: 9-node second order quadrangle (4 nodes associated with the vertices,
    4 with the edges and 1 with the face).
11: 10-node second order tetrahedron (4 nodes associated with the vertices and
    6 with the edges).
12: 27-node second order hexahedron (8 nodes associated with the vertices,
    12 with the edges, 6 with the faces and 1 with the volume).
13: 18-node second order prism (6 nodes associated with the vertices,
    9 with the edges and 3 with the quadrangular faces).
14: 14-node second order pyramid (5 nodes associated with the vertices,
    8 with the edges and 1 with the quadrangular face).
15: 1-node point.
16: 8-node second order quadrangle (4 nodes associated with the vertices and
    4 with the edges).
17: 20-node second order hexahedron (8 nodes associated with the vertices and
    12 with the edges).
18: 15-node second order prism (6 nodes associated with the vertices and 9
    with the edges).
19: 13-node second order pyramid (5 nodes associated with the vertices
    and 8 with the edges).
20: 9-node third order incomplete triangle (3 nodes associated with the
    vertices, 6 with the edges)
21: 10-node third order triangle (3 nodes associated with the vertices,
    6 with the edges, 1 with the face)
22: 12-node fourth order incomplete triangle (3 nodes associated with the
    vertices, 9 with the edges)
23: 15-node fourth order triangle (3 nodes associated with the vertices,
    9 with the edges, 3 with the face)
24: 15-node fifth order incomplete triangle (3 nodes associated with the vertices,
    12 with the edges)
25: 21-node fifth order complete triangle (3 nodes associated with the vertices,
    12 with the edges, 6 with the face)
26: 4-node third order edge (2 nodes associated with the vertices, 2 internal
    to the edge)
27: 5-node fourth order edge (2 nodes associated with the vertices, 3 internal
    to the edge)
28: 6-node fifth order edge (2 nodes associated with the vertices, 4 internal
    to the edge)
29: 20-node third order tetrahedron (4 nodes associated with the vertices,
    12 with the edges, 4 with the faces)
30: 35-node fourth order tetrahedron (4 nodes associated with the vertices,
    18 with the edges, 12 with the faces, 1 in the volume)
31: 56-node fifth order tetrahedron (4 nodes associated with the vertices,
    24 with the edges, 24 with the faces, 4 in the volume)
92: 64-node third order hexahedron (8 nodes associated with the vertices,
    24 with the edges, 24 with the faces, 8 in the volume)
93: 125-node fourth order hexahedron (8 nodes associated with the vertices,
    36 with the edges, 54 with the faces, 27 in the volume)
=#

export readGMSH

include("GMSHtypes.jl")

function readGMSH(filename::String)
    io = open(filename, "r")
    blockTypes, blocks = loopGMSHFileBlockSections(io)
    gmsh = evaluatGMSHBlockSections(blockTypes, blocks)

    return gmsh
end

function loopGMSHFileBlockSections(io)

    blocks = String[]
    blockTypes = String[]
    iobuff = IOBuffer(UInt8[], read=true, write=true, append=true)

    for line in eachline(io)
        if line[1] =='$'

            if line[1:4] == "\$End"
                push!(blocks, read(iobuff,String))
            else
                push!(blockTypes, line)
            end

            flush(iobuff)
        else
            write(iobuff, line * "\n")
        end
    end

    return blockTypes, blocks
end

function evaluatGMSHBlockSections(blockTypes::Array{String,1}, blocks)
    fileInfo=[]
    gmshPhysicalEntities=[]
    entities=[]
    nodes=[]
    elementGroups=[]

    for i=1:size(blockTypes,1)
        println(blockTypes[i])

        blockLines = split(blocks[i], "\n")
        blockLines = blockLines[blockLines .!= ""]

        if blockTypes[i] == "\$MeshFormat"
            fileInfo = evalMeshFormat(blockLines)
            println("File Version:",fileInfo.version_number)
            println("File Type:",fileInfo.file_type)
            println("Data Size:",fileInfo.data_size)

        elseif blockTypes[i] == "\$PhysicalNames"
            gmshPhysicalEntities = evalPhysicalNames(blockLines, fileInfo)

        elseif blockTypes[i] == "\$Entities"
            entities = evalEntities(blockLines, fileInfo)

        elseif blockTypes[i] == "\$Nodes"
            nodes = evalNodesBlock(blockLines, fileInfo)

        elseif blockTypes[i] == "\$Elements"
            elementGroups = evalElementsBlock(blockLines, fileInfo)

        end
    end

    # may not be defined in every gmsh file
    if typeof(gmshPhysicalEntities) == Array{Any,1}
        gmshPhysicalEntities = PhysicalEntities(0, nothing)
    end

    gmsh = GMSH(fileInfo, gmshPhysicalEntities, entities, nodes, elementGroups)

    return gmsh
end

function evalMeshFormat(blockContent::Array{SubString{String},1})

    if size(blockContent,1) <= 1
        s = split(blockContent[1], " ")

        fileInfo = GMSHFileInfo(parse(Float64, s[1]), parse(Int64, s[2]), parse(Int64, s[3]))

    else
        @error "Uknown size of MeshFormat Block in GMSHreader!"
    end

    return fileInfo
end

function evalPhysicalNames(blockContent::Array{SubString{String},1},
        fileinfo::GMSHFileInfo)

    count = parse(Int64,blockContent[1])
    entitiesBlock = blockContent[2:end]
    entities = Array{PhysicalEntity,1}(undef, size(entitiesBlock,1));

    for i = 1 : size(entitiesBlock,1)
        s = split(entitiesBlock[i], r" ")

        if size(s,1) > 3
            name = join(s[3:size(s,1)])
        else
            name = s[3];
        end

        entities[i] = PhysicalEntity(parse(UInt8, s[1]), parse(Int64, s[2]), name)
    end

    gmshPhysicalEntities = PhysicalEntities(parse(Int64, blockContent[1]) ,entities)
    return gmshPhysicalEntities
end

function evalEntityBlock(BlockContent::Array{SubString{String},1}, type::Int64, fileInfo::GMSHFileInfo)
    #=
    Entity Block String eval:
    type:
    1: Point
    2: Curve
    3: Surface
    4: volume
    =#

    if type == gmsh_tag_type["Point"]
        Entity = Array{Point,1}(undef, size(BlockContent,1))
    elseif type == gmsh_tag_type["Curve"]
        Entity = Array{Curve,1}(undef, size(BlockContent,1))
    elseif type == gmsh_tag_type["Surface"]
        Entity = Array{Surface,1}(undef, size(BlockContent,1))
    elseif type == gmsh_tag_type["Volume"]
        Entity = Array{Volume,1}(undef, size(BlockContent,1))
    else
        @error "unknow type"
    end

    for i=1:1:size(BlockContent,1)
        s = split(BlockContent[i], " ")
        tag = parse(Int64, s[1])

        if type == gmsh_tag_type["Point"]

            x = parse(Float64, s[2])
            y = parse(Float64, s[3])
            z = parse(Float64, s[4])

            if isapprox(fileInfo.version_number,4.0)
                # 8 = Bug in GMSH Doku oder Implementierung für Points...
                # Doku:
                # http://gmsh.info/doc/texinfo/gmsh.html#MSH-file-format
                #=
                pointTag(int) X(double) Y(double) Z(double)
                numPhysicalTags(size_t) physicalTag(int) ...
                ...
                Allerdings sind x y z doppelt besetzt ...
                =#
                numPhysicalTags = parse(Int64, s[8])
                physicalTag = Array{Int,1}(undef, numPhysicalTags)

            elseif isapprox(fileInfo.version_number,4.1)
                numPhysicalTags = parse(Int64, s[5])
                physicalTag = Array{Int,1}(undef, numPhysicalTags)
            else
                error("Unknown gmsh file version")
            end
        else
            xmin = parse(Float64, s[2])
            ymin = parse(Float64, s[3])
            zmin = parse(Float64, s[4])
            xmax = parse(Float64, s[5])
            ymax = parse(Float64, s[6])
            zmax = parse(Float64, s[7])
        end

        if type != gmsh_tag_type["Point"]
            numPhysicalTags = parse(Int64, s[8])
            physicalTag = Array{Int,1}(undef, numPhysicalTags)

            # TODO: what todo when points have physical tags?
            if numPhysicalTags > 0
                m = 1
                for k=9:1:(8+numPhysicalTags)
                    physicalTag[m] = parse(Int64, s[k])
                    m += 1
                end
            end

            numBoundingElements = parse(Int64, s[9+numPhysicalTags])
            boundingElementTag = Array{Int64,1}(undef, numBoundingElements)

            if numBoundingElements > 0
                m = 1
                for k=(10+numPhysicalTags):1:(9+numPhysicalTags+numBoundingElements)
                    boundingElementTag[m] = parse(Int64, s[k])
                    m += 1
                end
            end
        end

        if type == gmsh_tag_type["Point"]
            Entity[i] = Point(tag, x, y, z, numPhysicalTags, physicalTag)
        elseif type == gmsh_tag_type["Curve"]
            Entity[i] = Curve(tag, xmin, ymin, zmin, xmax, ymax, zmax,
                        numPhysicalTags, physicalTag, numBoundingElements,
                        boundingElementTag)
        elseif type == gmsh_tag_type["Surface"]
            Entity[i] = Surface(tag, xmin, ymin, zmin, xmax, ymax, zmax,
                        numPhysicalTags, physicalTag, numBoundingElements,
                        boundingElementTag)

        elseif type == gmsh_tag_type["Volume"]
            Entity[i] = Volume(tag, xmin, ymin, zmin, xmax, ymax, zmax,
                        numPhysicalTags, physicalTag, numBoundingElements,
                        boundingElementTag)
        end
    end

    return Entity
end

function evalEntities(blockContent::Array{SubString{String},1},fileinfo::GMSHFileInfo)

    s = split(blockContent[1], " ")

    if size(s,1) == 4
        numPoints = parse(Int64, s[1])
        numCurves = parse(Int64, s[2])
        numSurfaces = parse(Int64, s[3])
        numVolumes = parse(Int64, s[4])

        numAllEntities = numPoints+numCurves+numSurfaces+numVolumes

        if size(blockContent,1) == (1+numAllEntities)
            Points = evalEntityBlock(blockContent[2:(numPoints+1)],
                                        gmsh_tag_type["Point"], fileinfo)
            Curves = evalEntityBlock(
                            blockContent[(numPoints+2):(numPoints+numCurves+1)],
                            gmsh_tag_type["Curve"], fileinfo )
            Surfaces = evalEntityBlock(
                            blockContent[(numPoints+numCurves+2):(numPoints
                                                 +numCurves+numSurfaces+1)],
                            gmsh_tag_type["Surface"], fileinfo)
            Volumes = evalEntityBlock(blockContent[(numPoints+numCurves+
                                            numSurfaces+2):(numPoints+
                                            numCurves+numSurfaces+
                                            numVolumes+1)],
                                        gmsh_tag_type["Volume"], fileinfo)

            entities = MeshEntities(Points, Curves, Surfaces, Volumes)
            return entities
        else
            @error "Wrong number of lines in entities block"
        end

        println(numPoints, numCurves, numSurfaces, numVolumes)

    else
        @error "Unknown size of first line in entities block"
    end

end

function evalNodesBlock(blockContent::Array{SubString{String},1},fileinfo::GMSHFileInfo)
    println("evalNodes");
    #=
        size(blockContent[2:end]) - numEntities == numNodes
        that means that each Node is exclusive to its entity no double entities for nodes!
    =#
    s = split(blockContent[1], " ")
    numEntities = parse(Int64, s[1])
    numNodes = parse(Int64, s[2])

    println("numEntities: ", numEntities)
    println("numNodes: ", numNodes)

    x = zeros(Real, numNodes)
    y = zeros(Real, numNodes)
    z = zeros(Real, numNodes)
    nodenum = zeros(Int64, numNodes)
    nenteties = zeros(Int64, numNodes)
    ckecked = falses(numNodes)
    entityTags = zeros(Int64, numNodes)
    entityTypeIDs = zeros(Int64, numNodes)

    offset = 2;
    ni = 1
    for i = 1:numEntities
        si = split(blockContent[offset], " ")
        entityTag = parse(Int64, si[1])
        entityTypeID = parse(Int64, si[2])
        numEntityNodes = parse(Int64, si[end])

        if isapprox(fileinfo.version_number, 4.0)
            blockNodes=blockContent[(offset+1):(offset+numEntityNodes)]
            println("Node entity ", i, " has ", size(blockNodes,1), " nodes")

            for si in blockNodes
                sj = split(si, " ")
                node = parse(Int64, sj[1])

                if ckecked[ni]
                    @warn "dual cells checked, should not happen"
                else
                    ckecked[ni] = true
                end

                nodenum[ni] = node
                x[ni] = parse(Float64, sj[2])
                y[ni] = parse(Float64, sj[3])
                z[ni]= parse(Float64, sj[4])
                nenteties[ni] += 1
                entityTypeIDs[ni] = entityTypeID
                entityTags[ni] = entityTag

                ni += 1
            end

            offset = offset + numEntityNodes + 1

        elseif isapprox(fileinfo.version_number, 4.1)

            blockNodesNum = blockContent[(offset+1):(offset+numEntityNodes)]
            blockNodesCoords = blockContent[(offset+1+numEntityNodes):(offset+numEntityNodes+numEntityNodes)]
            println("Node entity ", i, " has ", size(blockNodesNum,1), " nodes")

            ni_temp = ni
            for si in blockNodesNum
                sj = split(si, " ")

                node = parse(Int64, sj[1])

                if ckecked[ni_temp]
                    @warn "dual cells checked, should not happen"
                else
                    ckecked[ni_temp] = true
                end

                nodenum[ni_temp] = node

                nenteties[ni_temp] += 1
                entityTypeIDs[ni_temp] = entityTypeID
                entityTags[ni_temp] = entityTag

                ni_temp += 1
            end

            for si in blockNodesCoords
                sj = split(si, " ")

                x[ni] = parse(Float64, sj[1])
                y[ni] = parse(Float64, sj[2])
                z[ni]= parse(Float64, sj[3])

                ni += 1
            end

            offset = offset + numEntityNodes + numEntityNodes + 1
        end
    end

    p = sortperm(nodenum)

    nodes = Nodes(nodenum[p], x[p], y[p], z[p], nenteties[p], entityTypeIDs[p],
                    entityTags[p])

    return nodes
end

function evalElementLine(blockContent::SubString{String}, elementTypeID::Int)
    s = split(blockContent, " ")
    s = s[s .!= ""]

    no_nodes =  elementTypeID_NodeCount[elementTypeID]
    number = parse(Int64, s[1])

    if no_nodes == size(s[2:end],1)
        nodes = Array{Int,1}(undef, no_nodes)

        for i = 1:no_nodes
            nodes[i] = parse(Int64,s[i+1])
        end

        element = Element(number, nodes)
        return element
    else
        @error "Element Type Node Count not equal Found Nodes"
    end

end

function evalElementBlock(blockContent::Array{SubString{String},1},entityTag::Int,
                     entityTypeID::Int, numEntityElements::Int, elementTypeID::Int64)

    elements = Array{Element, 1}(undef, numEntityElements)

    for i=1:1:numEntityElements
        elements[i] = evalElementLine(blockContent[i], elementTypeID)
    end

    elementGroup = ElementGroup(entityTag, entityTypeID, elementTypeID,
                                                numEntityElements, elements)

    return elementGroup
end

function evalElementsBlock(blockContent::Array{SubString{String},1},fileinfo::GMSHFileInfo)
    println("evalElements");

    s = split(blockContent[1], " ")
    numEntities = parse(Int64, s[1])
    numElements = parse(Int64, s[2])

    println("numEntities: ", numEntities)
    println("numElements: ", numElements)

    elementGroups = Array{ElementGroup,1}(undef, numEntities)

    offset = 2;
    for i = 1:numEntities

        si = split(blockContent[offset], " ")
        entityTag = parse(Int64, si[1])
        entityTypeID = parse(Int64, si[2])
        elementTypeID = parse(Int64, si[3])
        numEntityElements = parse(Int64, si[end])

        println("Elements entity ", i, " has ", numEntityElements, " elements")

        elementGroups[i] = evalElementBlock(blockContent[offset+1:offset+numEntityElements],
                            entityTag, entityTypeID, numEntityElements, elementTypeID)

        offset = offset + numEntityElements + 1
    end

    return elementGroups

end

end

max_stars_repo_name: c-schubert/linearMeshIO

max_issues_repo_name: c-schubert/linearMeshIO

max_forks_repo_name: c-schubert/linearMeshIO


text:
# Avoid ambiguity
function update_batch_vertex(g::GATConv, M::AbstractMatrix, X::CuMatrix, u)
    M = convert(typeof(X), M)
    update_batch_vertex(g, M, X, u)
end

function update_batch_vertex(g::GATConv, M::CuMatrix, X::CuMatrix, u)
    M = M .+ g.bias
    if !g.concat
        N = size(M, 2)
        M = reshape(mean(reshape(M, :, g.heads, N), dims=2), :, N)
    end
    return M
end


text:
export GoldState;
export push_cost, pop_cost, update_state;
export HEAD_ASSIGNED, HEAD_IN_BUFFER, HEAD_IN_STACK, HEAD_UNKNOWN;

const FORBIDDEN_COST = 9000

const HEAD_UNKNOWN = 0
const HEAD_ASSIGNED = 1
const HEAD_IN_STACK = 2
const HEAD_IN_BUFFER = 3

mutable struct GoldState
  gold_tree::DependencyTree
  system::ArcEagerSystem
  root_dependents_in_stack_count::Integer
  root_dependents_in_buffer_count::Integer
  dependents_in_stack_count::Vector{Integer}
  dependents_in_buffer_count::Vector{Integer}
  heads_states::Vector{Integer}
  config::Configuration

  function GoldState(gold_tree::DependencyTree, config::Configuration, system::ArcEagerSystem)
    state = new(gold_tree, system)
    update_state(state, config)

    state
  end
end

function update_state(state::GoldState, config::Configuration)
  state.config = config
  state.dependents_in_stack_count = zeros(state.gold_tree.length)
  state.dependents_in_buffer_count = zeros(state.gold_tree.length)
  state.root_dependents_in_buffer_count = 0
  state.root_dependents_in_stack_count = 0
  state.heads_states = zeros(state.gold_tree.length)

  # calculate dependents for root separatly
  foreach(config.buffer) do dependent
    if arc_present(state.gold_tree, 0, dependent)
      state.root_dependents_in_buffer_count += 1
    end
  end

  foreach(config.stack) do dependent
    if arc_present(state.gold_tree, 0, dependent)
      state.root_dependents_in_stack_count += 1
    end
  end

  # calculate dependents and head for each node
  for head = 1:state.gold_tree.length
    foreach(config.buffer) do dependent
      if arc_present(state.gold_tree, head, dependent)
        state.dependents_in_buffer_count[head] += 1
      end
    end

    foreach(config.stack) do dependent
      if arc_present(state.gold_tree, head, dependent)
        state.dependents_in_stack_count[head] += 1
      end
    end

    state.heads_states[head] = if head_in_stack(state, config, head)
      HEAD_IN_STACK
    elseif head_in_buffer(state, config, head)
      HEAD_IN_BUFFER
    elseif has_head(state.gold_tree, head)
      HEAD_ASSIGNED
    else
      HEAD_UNKNOWN
    end
  end
end

head_in_stack(state::GoldState, config::Configuration, word::Integer) = state.gold_tree.nodes[word].head_id in config.stack
head_in_buffer(state::GoldState, config::Configuration, word::Integer) = state.gold_tree.nodes[word].head_id in config.buffer

function dependents_in_buffer_count(state::GoldState, word_id::Integer)
  word_id == 0 ? state.root_dependents_in_buffer_count : state.dependents_in_buffer_count[word_id]
end

function dependents_in_stack_count(state::GoldState, word_id::Integer)
  word_id == 0 ? state.root_dependents_in_stack_count : state.dependents_in_stack_count[word_id]
end

function head_state(state::GoldState, word_id::Integer)
  word_id == 0 && return HEAD_UNKNOWN

  state.heads_states[word_id]
end

function label_correct(gold_state::GoldState, word_id::Integer, label::String)
  if 1 <= word_id <= gold_state.gold_tree.length
    node = gold_state.gold_tree.nodes[word_id]
    return node.label == label
  end

  false
end

max_stars_repo_name: LarsWl/DependencyParser

max_issues_repo_name: LarsWl/DependencyParser

max_forks_repo_name: LarsWl/DependencyParser


text:
module baobzi

mutable struct baobzi_t
end

mutable struct baobzi_input_t
    func::Ptr{Cvoid}
    data::Ptr{Cvoid}
    dim::Cint
    order::Cint
    tol::Cdouble
    minimum_leaf_fraction::Cdouble
    split_multi_eval::Cint
end

function init(fin, dim, order, center, half_length, tol, minimum_leaf_fraction=0.0, split_multi_eval=1)
    fanon = (x, p) -> fin(x)
    fbind = @cfunction($fanon, Cdouble, (Ptr{Cdouble}, Ptr{Cvoid},))
    input = baobzi_input_t(fbind.ptr, C_NULL, dim, order, tol, minimum_leaf_fraction, split_multi_eval)

    output_ptr = ccall(
        (:baobzi_init, :libbaobzi),
        Ptr{baobzi_t},
        (Ref{baobzi_input_t}, Ptr{Cdouble}, Ptr{Cdouble},),
        input, center, half_length
    )

    return output_ptr
end

function eval(func::Ptr{baobzi_t}, x)
    return ccall(
        (:baobzi_eval, :libbaobzi),
        Cdouble,
        (Ptr{Cvoid}, Ptr{Cdouble},),
        func, x
    )
end

function eval_multi(func::Ptr{baobzi_t}, x)
    n_points = Int32(length(x) // 2)
    println(n_points)
    res = Array{Float64}(undef, n_points)
    ccall(
        (:baobzi_eval_multi, :libbaobzi),
        Cvoid,
        (Ptr{Cvoid}, Ptr{Cdouble}, Ptr{Cdouble}, Cint,),
        func, x, res, n_points
    )
    return res
end

function stats(func::Ptr{baobzi_t})
    ccall(
        (:baobzi_stats, :libbaobzi),
        Ptr{Cvoid},
        (Ptr{Cvoid},),
        func
    )
end

function free(func::Ptr{baobzi_t})
    ccall(
        (:baobzi_free, :libbaobzi),
        Ptr{Cvoid},
        (Ptr{Cvoid},),
        func
    )
end

function save(func::Ptr{baobzi_t}, filename::String)
    ccall(
        (:baobzi_save, :libbaobzi),
        Ptr{Cvoid},
        (Ptr{Cvoid}, Cstring),
        func, filename
    )
end

function restore(filename::String)
    output_ptr = ccall(
        (:baobzi_restore, :libbaobzi),
        Ptr{baobzi_t},
        (Cstring,),
        filename
    )

    return output_ptr
end

end

max_stars_repo_name: blackwer/baobzi

max_issues_repo_name: blackwer/baobzi

max_forks_repo_name: blackwer/baobzi


text:
@testset "Reweighing" begin
	X, y, _ = @load_toydata;
	model = ConstantClassifier()
	wrappedModel = ReweighingWrapper(model; grp=:Sex)
	mach = machine(wrappedModel, X, y)
	fit!(mach)
	ŷ = predict(mach, X[1:5, :])
	@test all(x-> typeof(x) <: UnivariateFinite, ŷ)
	@test length(ŷ)==5
end

@testset "ReweighingSampling" begin
	X, y, _ = @load_toydata;
	model = ConstantClassifier()
	wrappedModel = ReweighingSamplingWrapper(model; grp=:Sex)
	mach = machine(wrappedModel, X, y)
	fit!(mach)
	ŷ = predict(mach, X[1:5, :])
	@test all(x-> typeof(x) <: UnivariateFinite, ŷ)
	@test length(ŷ)==5
end


text:
"""
    S = fditspec(sysr::FDFilterIF; FDfreq = missing, block = false, poleshift = false,
                 FDtol, FDStol, atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true)

Compute the weak or strong binary structure matrix `S` of the transfer function matrix `Rf(λ)` of
the transfer channel from the fault inputs to residuals of
a fault detection filter internal form object `sysr::FDFilterIF`.
For a filter `sysr` with `q` residual outputs and `mf` fault inputs,
`Rf(λ)` is the `q x mf` transfer function matrix of the fault inputs channel with the descriptor system representation
`sysr.sys[:,sysr.faults] := (Af-lambda*Ef,Bf,Cf,Df)`.

If `FDfreq = missing` (default), then `S` contains the weak structure matrix of `Rf(λ)`.
For `block = false`, `S` is determined as a `q x mf`
binary matrix (BitMatrix), whose `(i,j)`-th element is `S[i,j] = 1`, if the `(i,j)`-th element of `Rf(λ)` is
nonzero, and otherwise, `S[i,j] = 0`.
For `block = true`, `S` is determined as a `1 x mf` binary matrix, whose `(1,j)`-th element is `S[1,j] = 1`,
if the `j`-th column of `Rf(λ)` is nonzero, and otherwise, `S[1,j] = 0`.

If `FDfreq = freq` specifies a vector `freq` of `nf` real frequencies
which characterize the classes of persistent fault signals, then
for a suitable proper and invertible `M(λ)` (see below),
`S` contains the strong structure matrix of `M(λ)*Rf(λ)` with respect to a set of `nf` complex frequencies `Ω`, defined as follows:
if `f` is a real frequency in `freq`, then the corresponding complex frequency in `Ω`
is `λ := im*f`, for a continuous-time system,
or `λ := exp(im*f*abs(Ts))`, for a discrete-time system with sampling-time `Ts`.

`FDtol = tol1` specifies an absolute threshold `tol1` for the magnitudes of nonzero elements in the system matrices
`Bf` and `Df` and is used to determine the weak structure matrix.
Its default value is `tol1 = 0.0001*max(1, norm(Bf,1), norm(Df,1))`.

`FDStol = tol2` specifies an absolute  threshold `tol2` for the magnitudes of nonzero elements in the system matrices
`Af`, `Ef`, `Bf`, `Cf` and `Df` and is used to determine the strong structure matrix.
Its default value is
`tol2 = epsm*max(1, norm(Ef,1), norm(Af,1), norm(Bf,1), norm(Cf,Inf), norm(Df,1)))`,
where `epsm` is the working machine precision.

For `block = false`, then, if `poleshift = true`, `M(λ)` is chosen diagonal such that `M(λ)*Rf(λ)`
has no poles in `Ω` and if `poleshift = false` (default), `M(λ) = I` is used and
an error is issued if `Rf(λ)` has poles in `Ω`.
`S` is determined as a `q x mf` binary matrix, whose `(i,j)`-th element is `S[i,j] = 1`,
if the `(i,j)`-th element of `M(λ)*Rf(λ)`
evaluated for all frequencies in `Ω` is nonzero, and otherwise, `S[i,j] = 0`.

For `block = true`, then, if `poleshift = true`, `M(λ)` is chosen such that `M(λ)*Rf(λ)`
as no poles in `Ω` and if `poleshift = false` (default), `M(λ) = I` is used and
an error is issued if `Rf(λ)` has poles in `Ω`.
`S` is determined as an `1 x mf` binary matrix, whose `(1,j)`-th element is `S[1,j] = 1`,
if the `j`-th column of `M(λ)*Rf(λ)` evaluated for all frequencies in `Ω`
is nonzero and otherwise `S[1,j] = 0`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Af`, `Bf`, `Cf`, `Df`, the absolute tolerance for the nonzero elements of `Ef`,
and the relative tolerance for the nonzero elements of `Af`, `Bf`, `Cf`, `Df` and `Ef`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

_Method:_ For the definition of the structure matrix, see [1]. For the
determination of the weak structure matrix, minimal realizations
are determined for each column of `Rf(λ)` if `block = true` or for
each element of `Rf(λ)` if `block = false` and the nonzero columns or
elements in each column are identified  (see Corollary 7.1 of [1]).
For the determination of the strong structure matrix, minimal realizations
are determined for each column of `M(λ)*Rf(λ)` if `block = true` or for
each element of `M(λ)*Rf(λ)` if `block = false` and  the full rank of the
corresponding system matrix is checked for all frequencies in `FDfreq`
(see Corollary 7.2 in [1]) (i.e., the lack of zeros in all frequencies).

_References:_

[1] Varga A. Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017; sec.3.4.
"""
function fditspec(sysr::Union{FDFilterIF{T},FDIModel{T}};
    block::Bool = false, FDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing,
    poleshift::Bool = false, FDtol::Real = 0., FDStol::Real = 0.,
    atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol,
    rtol::Real = 0, fast::Bool = true) where T

    S = fditspec_(sysr.sys[:,sysr.faults]; block, FDfreq, poleshift, FDtol, FDStol,
                  atol1, atol2, rtol, fast)
    ismissing(FDfreq) && (return S)
    for j = 1:size(S,2)
        for i = 1:size(S,1)
            S[i,j,1] = all(view(S,i,j,:))
        end
    end
    return S[:,:,1]
end
"""
    S = fditspec(sysr::FDIFilterIF; FDfreq = missing, poleshift = false,
                 FDtol, FDStol, atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true)

Compute the weak or strong binary structure matrix `S` of the global transfer function matrix `Rf(λ)` of
the transfer channel from the fault inputs to residuals of
a fault detection and isolation filter internal form object `sysr::FDIFilterIF`.
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
the fault to residual channel of the `i`-th filter `sysr.sys[i][:,sysr.faults]`
has `qi` residual outputs and `mf` fault inputs, has the descriptor system representation
`sysr.sys[i][:,sysr.faults] := (Afi-lambda*Efi,Bfi,Cfi,Dfi)` and
`Rfi(λ)` is the corresponding `qi x mf` transfer function matrix.
The global transfer function matrix `Rf(λ)` is formed by row concatenation of the
transfer function matrices of the `N` individual filters, i.e.,  `Rf(λ) := [ Rf1(λ); Rf2(λ); ...; RfN(λ)]`.
For the evaluation of the strong structure matrix, the structure matrix of the stable
transfer function matrix `M(λ)*Rf(λ)` is determined, with a `M(λ)` block-diagonal
`M(λ) = block-diag(M1(λ), M2(λ), ..., MN(λ))`, where `Mi(λ)` is a suitable square and invertible
transfer function matrix (see below).

`FDtol = tol1` specifies an absolute threshold `tol1` for the magnitudes of nonzero elements in the system matrices
`Bf` and `Df` and is used to determine the weak structure matrix.
Its default value is `tol1 = 0.0001*max(1, norm(Bf,1), norm(Df,1))`.

`FDStol = tol2` specifies an absolute  threshold `tol2` for the magnitudes of nonzero elements in the system matrices
`Af`, `Ef`, `Bf`, `Cf` and `Df` and is used to determine the strong structure matrix.
Its default value is
`tol2 = epsm*max(1, norm(Ef,1), norm(Af,1), norm(Bf,1), norm(Cf,Inf), norm(Df,1)))`,
where `epsm` is the working machine precision.

If `FDfreq = missing` (default), then `S` contains the weak structure matrix of `Rf(λ)`.
`S` is determined as a `N x mf` binary matrix, whose `(i,j)`-th element is `S[i,j] = 1`,
if the `j`-th column of `Rfi(λ)` is nonzero, and otherwise, `S[i,j] = 0`.

If `FDfreq = freq` specifies a vector `freq` of `nf` real frequencies
which characterize the classes of persistent fault signals, then
for a suitable proper and invertible `M(λ)` (see below),
`S` contains the strong structure matrix of `M(λ)*Rf(λ)` with respect to a set of `nf` complex frequencies `Ω`, defined as follows:
if `f` is a real frequency in `freq`, then the corresponding complex frequency in `Ω`
is `λ := im*f`, for a continuous-time system,
or `λ := exp(im*f*abs(Ts))`, for a discrete-time system with sampling-time `Ts`.

`S` is determined as a `N x mf` binary matrix, whose `(i,j)`-th element is `S[i,j] = 1`,
if the `j`-th column of `Mi(λ)*Rfi(λ)` is nonzero for all frequencies in `Ω`, and otherwise, `S[i,j] = 0`.
If `poleshift = true`, `Mi(λ)` is chosen such that `Mi(λ)*Rfi(λ)` has no poles in `Ω` and
if `poleshift = false` (default), `Mi(λ) = I` is used and an error is issued if any `Rfi(λ)` has poles in `Ω`.

The keyword arguments `atol1`, `atol2`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Af`, `Bf`, `Cf`, `Df`, the absolute tolerance for the nonzero elements of `Ef`,
and the relative tolerance for the nonzero elements of `Af`, `Bf`, `Cf`, `Df` and `Ef`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol` and `atol2 = atol`.

_Method:_ For the definition of the structure matrix, see [1]. For the
determination of the weak structure matrix, minimal realizations
are determined for each column of `Rfi(λ)` and the nonzero columns are identified  (see Corollary 7.1 of [1]).
For the determination of the strong structure matrix, minimal realizations
are determined for each column of `Mi(λ)*Rfi(λ)` and  the full rank of the
corresponding system matrix is checked for all frequencies in `Ω`
(see Corollary 7.2 in [1]) (i.e., the lack of zeros in all frequencies in `Ω`).

_References:_

[1] Varga A. Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017; sec.3.4.
"""
function fditspec(sysr::FDIFilterIF{T}; FDfreq::Union{AbstractVector{<:Real},Real,Missing} = missing,
                  poleshift::Bool = false, FDtol::Real = 0., FDStol::Real = 0.,
                  atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol,
                  rtol::Real = 0, fast::Bool = true) where T
    N = length(sysr.sys)
    mf = length(sysr.faults)
    S = trues(N,mf)
    if ismissing(FDfreq)
       for i = 1:N
           S[i,:] = fditspec_(sysr.sys[i][:,sysr.faults]; block = true, FDtol, atol1, atol2, fast,
                                         rtol = (size(sysr.sys[i].A,1)+1)*eps(float(one(real(T))))*iszero(max(atol1,atol2,rtol)))
       end
    else
       isa(FDfreq,Vector) || (FDfreq = [FDfreq])
       lfreq = length(FDfreq);
       for i = 1:N
           Sc  = fditspec_(sysr.sys[i][:,sysr.faults]; block = true, FDfreq, poleshift, FDtol, FDStol, atol1, atol2, fast,
                                         rtol = (size(sysr.sys[i].A,1)+1)*eps(float(one(real(T))))*iszero(max(atol1,atol2,rtol)))
           lfreq == 1 ? (S[i,:] = Sc[1,:,1]) : ([S[i,j] = all(view(Sc,1,j,:)) for j = 1:mf])
       end
    end
    return S
end
"""
     S = fdisspec(sysr::FDFilterIF, freq; block = false, stabilize = false, FDGainTol = 0.01,
                     atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true)

Compute, for a given set of real frequencies `freq`,
the strong binary structure matrix `S` of the stable transfer function matrix `M(λ)*Rf(λ)`,
where `Rf(λ)` is the transfer function matrix of the
transfer channel from the fault inputs to residuals of the
fault detection filter internal form object `sysr::FDFilterIF` and
`M(λ)` is a suitable proper and invertible stabilizing transfer function matrix (see below).
For a filter `sysr` with `q` residual outputs and `mf` fault inputs,
`Rf(λ)` is the `q x mf` transfer function matrix of the fault inputs channel with the descriptor system representation
`sysr.sys[:,sysr.faults] := (Af-lambda*Ef,Bf,Cf,Df)`.

`freq` must contain a real frequency value or a vector of `nf` real frequencies
which characterize the classes of persistent fault signals
(default: `freq = 0`, i.e., characterizing constant faults).
`S` contains the strong
structure matrix of `M(λ)*Rf(λ)` with respect to a set of `nf` complex frequencies `Ω`, defined as follows:
if `f` is a real frequency in `freq`, then the corresponding complex frequency in `Ω`
is `λ := im*f`, for a continuous-time system,
or `λ := exp(im*f*abs(Ts))`, for a discrete-time system with sampling-time `Ts`.

`FDGainTol = tol` specifies an absolute  threshold `tol` for the nonzero magnitudes of
the frequency response gains (default: `tol = 0.01`).

For `block = false`, then, if `stabilize = true`, `M(λ)` is chosen diagonal such that `M(λ)*Rf(λ)`
has only stable poles and if `stabilize = false` (default), `M(λ) = I` is used and
an error is issued if `Rf(λ)` has poles in `Ω`.
`S` is determined as a `q x mf` binary matrix, whose `(i,j)`-th element is `S[i,j] = 1`,
if the `(i,j)`-th element of `M(λ)*Rf(λ)`
evaluated for all frequencies in `freq` is larger than or equal to `tol`, and otherwise, `S[i,j] = 0`.

For `block = true`, then, if `stabilize = true`, `M(λ)` is chosen such that `M(λ)*Rf(λ)`
has only stable poles and if `stabilize = false` (default), `M(λ) = I` is used and
an error is issued if `Rf(λ)` has poles in `Ω`.
`S` is determined as an `1 x mf` binary matrix, whose `(1,j)`-th element is `S[1,j] = 1`,
if the `j`-th column of `M(λ)*Rf(λ)` evaluated for all frequencies in `Ω`
is larger than or equal to `tol` and otherwise, `S[1,j] = 0`.

The keyword arguments `atol1`, `atol2`, `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Af`, `Bf`, `Cf`, `Df`, the absolute tolerance for the nonzero elements of `Ef`,
the absolute tolerance for the nonzero elements of `Cf`,
and the relative tolerance for the nonzero elements of `Af`, `Bf`, `Cf`, `Df` and `Ef`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`.

The computation of minimal realizations of individual input-output channels relies on pencil manipulation algorithms,
which employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

_Method:_ `S` is evaluated using the definition of the strong structure matrix in [1].

_References:_

[1] Varga A. Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017; sec. 3.4.

"""
function fdisspec(sysr::Union{FDFilterIF{T},FDIModel{T}}, freq::Union{AbstractVector{<:Real},Real} = 0; kwargs...) where T
    S = fdisspec_(sysr.sys[:,sysr.faults], freq; kwargs...)[1]
    length(freq) == 1 && (return S[:,:,1])
    for j = 1:size(S,2)
        for i = 1:size(S,1)
            S[i,j,1] = all(view(S,i,j,:))
        end
    end
    return S[:,:,1]
end
"""
     S = fdisspec(sysr::FDIFilterIF, freq; stabilize = false, FDGainTol = 0.01,
                     atol = 0, atol1 = atol, atol2 = atol, rtol, fast = true)

Compute, for a given set of real frequencies `freq`,
the strong binary structure matrix `S`
of the stable transfer function matrix `M(λ)*Rf(λ)`, where `Rf(λ)` is the global
transfer function matrix of the transfer channel from the fault inputs to residuals of the
fault detection and isolation filter internal form object `sysr::FDIFilterIF` and
`M(λ)` is a suitable block-diagonal proper and invertible stabilizing transfer function matrix (see below).
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
the fault to residual channel of the `i`-th filter `sysr.sys[i][:,sysr.faults]`
has `qi` residual outputs and `mf` fault inputs, has the descriptor system representation
`sysr.sys[i][:,sysr.faults] := (Afi-lambda*Efi,Bfi,Cfi,Dfi)` and
`Rfi(λ)` is the corresponding `qi x mf` transfer function matrix.
The global transfer function matrix `Rf(λ)` is formed by row concatenation of the
transfer function matrices of the `N` individual filters, i.e.,  `Rf(λ) := [ Rf1(λ); Rf2(λ); ...; RfN(λ)]`.
`M(λ) = block-diag(M1(λ), M2(λ), ..., MN(λ))`, where `Mi(λ)` is square and invertible
and chosen such that `Mi(λ)Rfi(λ)` is stable (see below).

`freq` must contain a real frequency value or a vector of `nf` real frequencies
which characterize the classes of persistent fault signals
(default: `freq = 0`, i.e., characterizing constant faults).
`S` contains the strong
structure matrix of `M(λ)*Rf(λ)` with respect to a set of `nf` complex frequencies `Ω`, defined as follows:
if `f` is a real frequency in `freq`, then the corresponding complex frequency in `Ω`
is `λ := im*f`, for a continuous-time system,
or `λ := exp(im*f*abs(Ts))`, for a discrete-time system with sampling-time `Ts`.

`FDGainTol = tol` specifies an absolute  threshold `tol` for the nonzero magnitudes of
the frequency response gains (default: `tol = 0.01`).

If `stabilize = true`, `Mi(λ)` is chosen such that `Mi(λ)*Rfi(λ)` has only stable poles and
if `stabilize = false` (default), `Mi(λ) = I` is used and an error is issued
if any `Rfi(λ)` has poles in `Ω`.

`S` is determined as a `N x mf`
binary matrix, whose `(i,j)`-th element is `S[i,j] = 1`, if the norm of the
`j`-th column of `Mi(λ)*Rfi(λ)` evaluated for all frequencies in `Ω`
is larger than or equal to `tol`, and otherwise, `S[i,j] = 0`.

The keyword arguments `atol1`, `atol2`, `atol3`, and `rtol`, specify, respectively, the absolute tolerance for the
nonzero elements of matrices `Afi`, `Bfi`, `Cfi`, `Dfi`, the absolute tolerance for the nonzero elements of `Efi`,
the absolute tolerance for the nonzero elements of `Cfi`,
and the relative tolerance for the nonzero elements of `Afi`, `Bfi`, `Cfi`, `Dfi` and `Efi`.
The keyword argument `atol` can be used to simultaneously set `atol1 = atol`, `atol2 = atol` and `atol3 = atol`.

The computation of minimal realizations of individual input-output channels relies on pencil manipulation algorithms,
which employ rank determinations based on either the use of
rank revealing QR-decomposition with column pivoting, if `fast = true`, or the SVD-decomposition.
The rank decision based on the SVD-decomposition is generally more reliable, but the involved computational effort is higher.

_Method:_ `S` is evaluated using the definition of the strong structure matrix in [1].

_References:_

[1] Varga A. Solving Fault Diagnosis Problems - Linear Synthesis Techniques. Springer Verlag, 2017; sec. 3.4.

"""
function fdisspec(sysr::FDIFilterIF{T}, freq::Union{AbstractVector{<:Real},Real} = 0; stabilize::Bool = false,
                  FDGainTol::Real = 0.01, atol::Real = zero(float(real(T))), atol1::Real = atol, atol2::Real = atol,
                  rtol::Real = 0, fast::Bool = true) where T
    isa(freq,Vector) || (freq = [freq])
    lfreq = length(freq);
    N = length(sysr.sys)
    mf = length(sysr.faults)
    S = trues(N,mf)
    for i = 1:N
        Sc = fdisspec_(sysr.sys[i][:,sysr.faults], freq; block = true, stabilize, FDGainTol,
                       atol1, atol2, fast, rtol = (size(sysr.sys[i].A,1)+1)*eps(float(one(real(T))))*iszero(max(atol1,atol2,rtol)))[1]
        lfreq == 1 ? (S[i,:] = Sc[1,:,1]) : ([S[i,j] = all(view(Sc,1,j,:)) for j = 1:mf])
    end
    return S
end
"""
     fdiscond(sysr::FDFilterIF,freq) -> (scond, β, γ)

Compute for the stable transfer function matrix `Rf(λ)` of the
transfer channel from the fault inputs to residuals of the
fault detection filter internal form object `sysr::FDFilterIF` the quantities:
`β` - the H∞- index of `Rf(λ)`, `γ` - the maximum of the columns norms of `Rf(λ)` and
the fault detection sensitivity condition `scond` evaluated as `scond := β/γ`.
If `freq` is a vector of real frequency values, then `β` and `γ`
are evaluated over the frequencies contained in `freq`.
"""
function fdiscond(sysr::FDFilterIF{T}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T
    fdiscond_(sysr.sys[:,sysr.faults], freq)
end
"""
     fdiscond(sysr::FDFilterIF,SFDI,freq) -> (scond, β, γ)

Compute the detection and isolation sensitivity condition
`scond` (and related quatities `β` and `γ`) for the binary structure vector `SFDI` associated to
the stable transfer function matrix `Rf(λ)` of the
transfer channel from the fault inputs to residuals of the
fault detection filter internal form object `sysr::FDFilterIF`.
If  `Rff(λ)` is the transfer function matrix formed of those `j`-th columns of `Rf(λ)`
for which `SFDI[j] = 1`, then:
`β` - the H∞- index of `Rff(λ)`, `γ` - the maximum of the columns norms of `Rff(λ)` and
the fault detection sensitivity condition `scond` evaluated as `scond := β/γ`.
If `freq` is a vector of real frequency values, then `β` and `γ`
are evaluated over the frequencies contained in `freq`.
"""
function fdiscond(sysr::FDFilterIF{T}, SFDI::Union{BitVector,Vector{Bool}}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T
    fdiscond_(sysr.sys[:,sysr.faults[SFDI]], freq)
end
"""
     fdiscond(sysr::FDFilterIF, SFDI, freq) -> (scond, β, γ)

Compute the detection and isolation sensitivity condition vector
`scond` (and related quatities `β` and `γ`) for the `q × mf` binary structure matrix `SFDI` associated to
the stable transfer function matrix `Rf(λ)` of the
transfer channel from the fault inputs to residuals of the
fault detection filter internal form object `sysr::FDFilterIF`.
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
`sysr.sys[i][:,sysr.faults]` is the fault to residual channel of the `i`-th filter and
`Rfi(λ)` is the corresponding transfer function matrix.
The global transfer function matrix `Rf(λ)` is formed by row concatenation of the
transfer function matrices of the `N` individual filters, i.e.,  `Rf(λ) := [ Rf1(λ); Rf2(λ); ...; RfN(λ)]`.
The i-th element of the vectors `scond`, `β` and `γ` contain the quantities:
`β[i]` - the H∞- index of the nonzero columns of `Rffi(λ)`, `γ` - the maximum of the nonzero columns norms of `Rfi(λ)` and
the correspomding fault detection sensitivity condition `scond[i]` evaluated as `scond[i] := β[i]/γ[i]`,
where `Rffi(λ)` is formed of those `j`-th  columns of the `i`-th row of `Rf(λ)` for which `S[i,j] = 1` and
`Rfi(λ)` is the `i`-th row of `Rf(λ)`.
It is assumed that for each `j` such that `SFDI[i,j] = 1`, the `j`-th column of `Rfi(λ)` is nonzero  and
for each `j` such that `SFDI[i,j] = 0`, the `j`-th column of `Rfi(λ)` is zero.
If `freq` is a vector of real frequency values, then `β[i]` and `γ[i]`
are evaluated over the frequencies contained in `freq`.
"""
function fdiscond(sysr::FDFilterIF{T}, SFDI::Union{BitMatrix,Array{Bool,2}}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T
    q = size(sysr.sys,1)
    inpf = sysr.faults
    mf = length(inpf)
    (q, mf ) == size(SFDI) || error("missmatch between the number of fault inputs and column dimension of SFDI")

    mf == 0 && (return zeros(T,q,0), zeros(T,q,0), zeros(T,q,0))
    scond = similar(Array{T,1},q)
    β = similar(Array{T,1},q)
    γ = similar(Array{T,1},q)
    for i = 1:q
        scond[i], β[i], γ[i] = fdiscond_(sysr.sys[i,inpf[view(SFDI,i,:)]], freq)
    end
    return scond, β, γ
end
"""
     fdiscond(sysr::FDIFilterIF, SFDI, freq) -> (scond, β, γ)

Compute the detection and isolation sensitivity condition
`scond` (and related quatities `β` and `γ`) for the `N × mf` binary structure matrix `SFDI` associated to
the stable global transfer function matrix `Rf(λ)` of the
transfer channel from the fault inputs to residuals of the
fault detection and isolation filter internal form object `sysr::FDIFilterIF`.
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
`sysr.sys[i][:,sysr.faults]` is the fault to residual channel of the `i`-th filter and
`Rfi(λ)` is the corresponding transfer function matrix.
The global transfer function matrix `Rf(λ)` is formed by row concatenation of the
transfer function matrices of the `N` individual filters, i.e.,  `Rf(λ) := [ Rf1(λ); Rf2(λ); ...; RfN(λ)]`.
It is assumed that for each `j` such that `SFDI[i,j] = 1`, the `j`-th column of `Rfi(λ)` is nonzero  and
for each `j` such that `SFDI[i,j] = 0`, the `j`-th column of `Rfi(λ)` is zero.
The i-th element of the vectors `scond`, `β` and `γ` contain the quantities:
`β[i]` - the H∞- index of the nonzero columns of `Rfi(λ)`, `γ` - the maximum of the nonzero columns norms of `Rf(λ)` and
the fault detection sensitivity condition `scond` evaluated as `scond[i] := β[i]/γ[i]`.
If `freq` is a vector of real frequency values, then `β[i]` and `γ[i]`
are evaluated over the frequencies contained in `freq`.
"""
function fdiscond(sysr::FDIFilterIF{T}, SFDI::Union{BitMatrix,BitVector,Array{Bool,2},Array{Bool,1}}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T
    N = length(sysr.sys)
    isa(SFDI,Union{BitVector,Array{Bool,1}}) ? (nb = 1; mf = length(SFDI); SFDI = reshape(SFDI,nb,mf) ) :
                                               (nb = size(SFDI,1); mf = size(SFDI,2))
    nb == N || error("missmatch between number of filters $N and number of specifications $nb")
    inpf = sysr.faults
    mf == length(inpf) || error("missmatch between number of faults and column dimension/length of SFDI")
    scond = similar(Array{T,1},N)
    β = similar(Array{T,1},N)
    γ = similar(Array{T,1},N)
    for i = 1:N
        scond[i], β[i], γ[i] = fdiscond_(sysr.sys[i][:,inpf[view(SFDI,i,:)]], freq)
    end
    return scond, β, γ
end
fdiscond(sysr::FDIFilterIF{T}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T = fdiscond(sysr,trues(1,length(sysr.faults)),freq)

"""
     fdif2ngap(sysr::FDFilterIF, freq) -> (gap, β, γ)

Compute the fault-to-noise gap `gap` (and the related quantities `β` and `γ`)
for the stable fault detection filter internal form object `sysr::FDFilterIF`.
For the fault to residual channel of the filter `sysr.sys[:,sysr.faults]`
with the corresponding transfer function matrix `Rf(λ)` and
the noise to residual channel of the filter `sysr.sys[:,sysr.noise]`
with the corresponding transfer function matrix `Rw(λ)`,
`β` is the H∞- index of `Rf(λ)`, `γ` is the H∞-norm of `Rw(λ)` and
`gap` is the fault-to-noise gap evaluated as `gap := β/γ`.
If `freq` is a vector of real frequency values, then `β` and `γ`
are evaluated over the frequencies contained in `freq`.
`gap = ∞` if there are no noise inputs and `gap = 0` if there are no fault inputs.
"""
function fdif2ngap(sysr::FDFilterIF{T}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T
   mf = length(sysr.faults)
   mw = length(sysr.noise)

   mf+mw == 0 && (return T[], T[], T[])

   isstable(sysr.sys) || error("the system is unstable")

   β = mf == 0 ? 0 : fdhinfminus(sysr.sys[:,sysr.faults],freq)[1]
   γ = mw == 0 ? 0 : ghinfnorm(sysr.sys[:,sysr.noise])[1]
   return β/γ, β, γ
end
"""
     fdif2ngap(sysr::FDFilterIF, SFDI, freq) -> (gap, β, γ)

Compute the fault-to-noise gap `gap` (and the related quantities `β` and `γ`)
for the stable fault detection filter internal form object `sysr::FDFilterIF`
and the associated binary structure vector `SFDI`.
`sysr.sys[:,sysr.faults]` is the fault to residual channel of the filter
with the corresponding transfer function matrix `Rf(λ)` and
`sysr.sys[:,sysr.noise]` is the noise to residual channel of the filter
with the corresponding transfer function matrix `Rw(λ)`.
If  `Rff(λ)` is the transfer function matrix formed of those `j`-th columns of `Rf(λ)`
for which `SFDI[j] = 1` and `Rdf(λ)` is the transfer function matrix formed of
those `j`-th columns of `Rf(λ)` for which `SFDI[j] = false`, then:
`β` is the H∞- index of `Rff(λ)`, `γ` is the H∞-norm of `[Rdf(λ) Rw(λ)]` and
`gap` is the fault-to-noise gap evaluated as `gap := β/γ`.
If `freq` is a vector of real frequency values, then `β` and `γ`
are evaluated over the frequencies contained in `freq`.
`gap = ∞` if `[Rdf(λ) Rw(λ)] = 0` and `gap = 0` if `Rff(λ) = 0`.
"""
function fdif2ngap(sysr::FDFilterIF{T}, SFDI::Union{BitVector,Vector{Bool}}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing) where T
   inpf = sysr.faults
   mf = length(inpf)
   mw = length(sysr.noise)

   mf+mw == 0 && (return T[], T[], T[])

   isstable(sysr.sys) || error("the system is unstable")

   β = mf == 0 ? 0 : fdhinfminus(sysr.sys[:,inpf[SFDI]],freq)[1]
   γ = ghinfnorm(sysr.sys[:,[inpf[.!SFDI]; sysr.noise]])[1]
   return β/γ, β, γ
end
"""
      fdif2ngap(sysr::FDFilterIF, SFDI, freq; atol = √ϵ) -> (gap, β, γ)

Compute the fault-to-noise gap
`gap` (and the related quatities `β` and `γ`) for the `q × mf` binary structure matrix `SFDI` associated to
the stable transfer function matrices `Rf(λ)` and `Rw(λ)` of the
transfer channels from the fault inputs to residuals and noise inputs to residuals, respectively,  of the
fault detection filter internal form object `sysr::FDFilterIF`.
The `i`-th element of the vectors `gap`, `β` and `γ` contain the quantities:
`β[i]` - the H∞- index of `Rffi(λ)`, `γ[i]` - the H∞-norm of `[Rdfi(λ) Rwi(λ)]` and
the fault-to-noise gap `gap` evaluated as `gap[i] := β[i]/γ[i]`, where
`Rffi(λ)` is formed of those `j`-th  columns of the `i`-th row of `Rf(λ)` for which `S[i,j] = 1` and
`Rdfi(λ)` is formed of those `j`-th  columns of the `i`-th row of `Rf(λ)` for which `S[i,j] = 0`.
`gap[i] = ∞` if `[Rdfi(λ) Rwi(λ)] = 0` and `gap[i] = 0` if `Rffi(λ) = 0`.
If `freq` is a vector of real frequency values, then `β[i]` and `γ[i]`
are evaluated over the frequencies contained in `freq`.
`atol` is an absolute tolerance for the norms `Rwi(λ)`, such that norm values less than or equal to `atol` are
considered zero (default:  `√ϵ`, where `ϵ` is the working machine precision.)
"""
function fdif2ngap(sysr::FDFilterIF{T}, SFDI::Union{BitMatrix,Matrix{Bool}}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing;
                    atol::Real = sqrt(eps(T))) where T
    q = size(sysr.sys,1)
    inpf = sysr.faults
    mf = length(inpf)
    (q, mf ) == size(SFDI) || error("missmatch between the number of fault inputs and column dimension of SFDI")
    mw = length(sysr.noise)

    mf+mw == 0 && (return zeros(T,q,0), zeros(T,q,0), zeros(T,q,0))

    isstable(sysr.sys) || error("the system is unstable")
    β = similar(Array{T,1},q)
    γ = similar(Array{T,1},q)
    for i = 1:q
        β[i] = mf == 0 ? 0 : fdhinfminus(sysr.sys[i,inpf[SFDI[i,:]]],freq)[1]
        t = ghinfnorm(sysr.sys[i,[inpf[.!SFDI[i,:]]; sysr.noise]])[1]
        γ[i] = t <= atol ? 0 : t
    end
    return β ./ γ, β, γ
 end
 """
     fdif2ngap(sysr::FDIFilterIF, SFDI, freq; atol = √ϵ) -> (gap, β, γ)

Compute the fault-to-noise gap
`gap` (and the related quatities `β` and `γ`) for the `N × mf` binary structure matrix `SFDI` associated to
the stable global transfer function matrices `Rf(λ)` and `Rw(λ)` of the
transfer channels from the fault inputs to residuals and noise inputs to residuals, respectively,  of the
fault detection and isolation filter internal form object `sysr::FDIFilterIF`.
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
`sysr.sys[i][:,sysr.faults]` is the fault to residual channel of the `i`-th filter
with the corresponding transfer function matrix `Rfi(λ)` and `sysr.sys[i][:,sysr.noise]` is the noise to residual channel of the `i`-th filter
with the corresponding transfer function matrix `Rwi(λ)`.
The global transfer function matrices `Rf(λ)` and `Rw(λ)` are formed by row concatenation of the
transfer function matrices of the `N` individual filters, i.e.,  `Rf(λ) := [ Rf1(λ); Rf2(λ); ...; RfN(λ)]` and
`Rw(λ) := [ Rw1(λ); Rw2(λ); ...; RwN(λ)]`
Let  `Rffi(λ)` be the transfer function matrix formed of those `j`-th columns of `Rfi(λ)`
for which `SFDI[i,j] = 1` and let `Rdfi(λ)` be the transfer function matrix formed of those `j`-th columns of `Rfi(λ)`
for which `SFDI[i,j] = 0`.
The `i`-th element of the vectors `gap`, `β` and `γ` contain the quantities:
`β[i]` - the H∞- index of `Rffi(λ)`, `γ[i]` - the H∞-norm of `[Rdfi(λ) Rwi(λ)]` and
the fault-to-noise gap `gap` evaluated as `gap[i] := β[i]/γ[i]`.
`gap[i] = ∞` if `[Rdfi(λ) Rwi(λ)] = 0` and `gap[i] = 0` if `Rffi(λ) = 0`.
If `freq` is a vector of real frequency values, then `β[i]` and `γ[i]`
are evaluated over the frequencies contained in `freq`.
`atol` is an absolute tolerance for the norms `Rwi(λ)`, such that norm values less than or equal to `atol` are
considered zero (default:  `√ϵ`, where `ϵ` is the working machine precision.)
"""
function fdif2ngap(sysr::FDIFilterIF{T}, SFDI::Union{BitMatrix,BitVector,Array{Bool,2},Array{Bool,1}}, freq::Union{AbstractVector{<:Real},Real,Missing} = missing;
                   atol::Real = sqrt(eps(T))) where T
    N = length(sysr.sys)
    isa(SFDI,Union{BitVector,Array{Bool,1}}) ? (nb = 1; mf = length(SFDI); SFDI = reshape(SFDI,nb,mf) ) :
                                               (nb = size(SFDI,1); mf = size(SFDI,2))
    nb == N || error("missmatch between the number of filters $N and the number of specifications $nb")
    inpf = sysr.faults
    mf == length(inpf) || error("missmatch between the number of faults and the column dimension/length of SFDI")
    mw = length(sysr.noise)
    β = mf > 0 ? similar(Array{T,1},N) : zeros(T,N)
    γ = mf+mw > 0 ? similar(Array{T,1},N) : zeros(T,N)
    if mf > 0
        for i = 1:N
            β[i] = fdhinfminus(sysr.sys[i][:,inpf[view(SFDI,i,:)]],freq)[1]
        end
    end
    if mf+mw > 0
        for i = 1:N
            t = ghinfnorm(sysr.sys[i][:,[inpf[.!SFDI[i,:]]; sysr.noise]])[1]
            γ[i] = t <= atol ? 0 : t
        end
    end

    return β./γ, β, γ
end
"""
     γ = fdimmperf(sysr::FDFilterIF[, nrmflag])

Compute the model-matching performance `γ` of the fault detection filter internal form object `sysr::FDFilterIF`.
If `Rw(λ)` is the transfer function matrix of the transfer channel from the noise inputs to residuals
`sysr.sys[:,sysr.noise]`, then
`γ` is the  H∞-norm of `Rw(λ)`, if `nrmflag = Inf` (default) and the  H2-norm of `Rw(λ)`, if `nrmflag = 2`.
The value of `γ` is infinite for an unstable filter or if `nrmflag = 2` and the transfer function matrix
`Rw(λ)` of a continuous-time system is not strictly proper.
"""
function fdimmperf(sysr::FDFilterIF{T}, nrmflag::Real = Inf) where T

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")

    length(sysr.noise) == 0 && (return T(0))

    return isinf(nrmflag) ? ghinfnorm(sysr.sys[:,sysr.noise])[1] :
                            gh2norm(sysr.sys[:,sysr.noise])
end
"""
     γ = fdimmperf(sysr::FDFilterIF, SFDI[, nrmflag])

Compute the model-matching performance `γ` of the fault detection filter internal form object `sysr::FDFilterIF`
for a given binary structure vector `SFDI`. If `Rf(λ)` is the transfer function matrix of the
transfer channel from the fault inputs to residuals `sysr.sys[:,sysr.faults]` and
`Rw(λ)` is the transfer function matrix of the transfer channel from the noise inputs to residuals
`sysr.sys[:,sysr.noise]`, then
`γ` is the  H∞-norm of `[Rdf(λ) Rw(λ)]`, if `nrmflag = Inf` (default) and the  H2-norm of `[Rdf(λ) Rw(λ)]`,
if `nrmflag = 2`, where `Rdf(λ)` is the transfer function matrix formed by those `j`-th columns of `Rf(λ)` for which
`SFDI[j] = 0`.
The value of `γ` is infinite for an unstable filter or if `nrmflag = 2` and the transfer function matrix
`[Rdf(λ) Rw(λ)]` of a continuous-time system is not strictly proper.
"""
function fdimmperf(sysr::FDFilterIF{T}, SFDI::Union{BitVector,Vector{Bool}}, nrmflag::Real = Inf) where T

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")
    inpf = sysr.faults
    length(SFDI) == length(inpf) || error("missmatch between the number of faults and the column dimension/length of SFDI")

    return isinf(nrmflag) ? ghinfnorm(sysr.sys[:,[inpf[.!SFDI]; sysr.noise]])[1] :
                            gh2norm(sysr.sys[:,[inpf[.!SFDI]; sysr.noise]])
end
"""
     γ = fdimmperf(sysr::FDFilterIF, SFDI[, nrmflag])

Compute the model-matching performance `γ` of the fault detection filter internal form object `sysr::FDFilterIF`
for a given binary structure matrix `SFDI`. If `Rf(λ)` is the transfer function matrix of the
transfer channel from the fault inputs to residuals `sysr.sys[:,sysr.faults]` and
`Rw(λ)` is the transfer function matrix of the transfer channel from the noise inputs to residuals
`sysr.sys[:,sysr.noise]`, then
`γ` is the  H∞-norm of `[Rdf(λ) Rw(λ)]`, if `nrmflag = Inf` (default) and the  H2-norm of `[Rdf(λ) Rw(λ)]`,
if `nrmflag = 2`, where `Rdf(λ) = .!SFDI .* Rf(λ)` (i.e., the element-wise product of `.!SFDI` and `Rf(λ)`.
The value of `γ` is infinite for an unstable filter or if `nrmflag = 2` and the transfer function matrix
`[Rdf(λ) Rw(λ)]` of a continuous-time system is not strictly proper.
"""
function fdimmperf(sysr::FDFilterIF{T}, SFDI::Union{BitMatrix,Matrix{Bool}}, nrmflag::Real = Inf) where T

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")

    return isinf(nrmflag) ? ghinfnorm([dssubsel(sysr.sys[:,sysr.faults], .!SFDI) sysr.sys[:,sysr.noise]])[1] :
                            gh2norm([dssubsel(sysr.sys[:,sysr.faults], .!SFDI) sysr.sys[:,sysr.noise]])
end

"""
     γ = fdimmperf(sysr::FDIFilterIF[, nrmflag])

Compute the model-matching performance `γ` of the stable fault detection and isolation filter internal form object `sysr::FDIFilterIF`.
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
`sysr.sys[i][:,sysr.noise]` is the noise to residual channel of the `i`-th filter
with the corresponding transfer function matrix `Rwi(λ)`. Then,
`γ` is an `N`-dimensional vector whose `i`-th component is the  H∞-norm of `Rwi(λ)`, if `nrmflag = Inf` (default)
and the  H2-norm of `Rwi(λ)`, if `nrmflag = 2`.
The `i`-th component of `γ` is infinite for an unstable filter or if `nrmflag = 2` and the transfer function matrix
`Rwi(λ)` of a continuous-time system is not strictly proper.
"""
function fdimmperf(sysr::FDIFilterIF{T}, nrmflag::Real = Inf) where T
    N = length(sysr.sys)

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")

    length(sysr.noise) == 0 && (return zeros(T,N))

    γ = similar(Array{T,1},N)
    for i = 1:N
        γ[i] = isinf(nrmflag) ? ghinfnorm(sysr.sys[i][:,sysr.noise])[1] :
                                gh2norm(sysr.sys[i][:,sysr.noise])
    end
    return γ
end
"""
     γ = fdimmperf(sysr::FDIFilterIF, SFDI[, nrmflag])

Compute the model-matching performance `γ` of the stable fault detection and isolation
filter internal form object `sysr::FDIFilterIF` and the associated binary structure matrix `SFDI`.
The filter `sysr` consists of `N` individual FDI filters `sysr.sys[i]`, for `i = 1, ..., N`, where
`sysr.sys[i][:,sysr.faults]` is the fault to residual channel of the `i`-th filter
with the corresponding transfer function matrix `Rfi(λ)` and
`sysr.sys[i][:,sysr.noise]` is the noise to residual channel of the `i`-th filter
with the corresponding transfer function matrix `Rwi(λ)`. Then,
`γ` is an `N`-dimensional vector whose `i`-th component is the  H∞-norm of `[Rfdi(λ) Rwi(λ)]`, if `nrmflag = Inf` (default)
or the  H2-norm of `[Rfdi(λ) Rwi(λ)]`, if `nrmflag = 2`, where `Rfdi(λ)` is the transfer function matrix whose `j`-th column is
zero if `SFDI[i,j] = 1` and is equal to the `j`-th column of `Rfi(λ)` if `SFDI[i,j] = 0`.
The `i`-th component of `γ` is infinite for an unstable filter or if `nrmflag = 2` and the transfer function matrix
`[Rfdi(λ) Rwi(λ)]` of a continuous-time system is not strictly proper.
"""
function fdimmperf(sysr::FDIFilterIF{T}, SFDI::Union{BitMatrix,BitVector,Array{Bool,2},Array{Bool,1}}, nrmflag::Real = Inf) where T
    N = length(sysr.sys)

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")
    isa(SFDI,Union{BitVector,Array{Bool,1}}) ? (nb = 1; mf = length(SFDI); SFDI = reshape(SFDI,nb,mf) ) :
                                               (nb = size(SFDI,1); mf = size(SFDI,2))
    nb == N || error("missmatch between the number of filters $N and the number of specifications $nb")
    inpf = sysr.faults
    mf == length(inpf) || error("missmatch between the number of faults and the column dimension/length of SFDI")

    γ = similar(Array{T,1},N)
    for i = 1:N
        inddfw = [inpf[.!SFDI[i,:]]; sysr.noise]
        γ[i] = isinf(nrmflag) ? ghinfnorm(sysr.sys[i][:,inddfw])[1] :
                                gh2norm(sysr.sys[i][:,inddfw])
    end
    return γ
end
"""
     γ = fdimmperf(sysr::FDFilterIF, sysref::Union{FDFilterIF,FDIModel}[, nrmflag])

Compute the model-matching performance `γ` of the fault detection filter internal form object `sysr::FDFilterIF`
with respect to the fault detection reference filter internal form  `sysref::FDFilterIF`.
If `R(λ)` is the transfer function matrix of the fault detection filter internal form
`sysr.sys` and `Mr(λ)` is the transfer function matrix of the fault detection reference filter internal form
`sysref.sys`, then
`γ` is the  H∞-norm of `R(λ)-Mr(λ)`, if `nrmflag = Inf` (default) or the  H2-norm of `R(λ)-Mr(λ)`, if `nrmflag = 2`.
The value of `γ` is infinite for an unstable difference `R(λ)-Mr(λ)` or if `nrmflag = 2` and the transfer function matrix
`R(λ)-Mr(λ)` of a continuous-time system is not strictly proper. In general, `R(λ)` and `Mr(λ)` are partitioned as
`R(λ) = [ Ru(λ) Rd(λ) Rf(λ) Rw(λ) Ra(λ) ]` and `Mr(λ) = [ Mru(λ) Mrd(λ) Mrf(λ) Mrw(λ) Mra(λ) ]` in accordance with
the partitioning of the inputs in control inputs, disturbance inputs, fault inputs, noise inputs and auxiliary inputs.
Void components of `Mr(λ)` corresponding to non-void components in `R(λ)` are assumed to be zero.
"""
function fdimmperf(sysr::FDFilterIF{T}, sysref::Union{FDFilterIF,FDIModel}, nrmflag::Real = Inf; atolinf::Real = 0) where T

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")
    DescriptorSystems.promote_Ts(sysr.sys.Ts, sysref.sys.Ts)  # check sampling time

    inpu = sysr.controls; mu = length(inpu)
    inpd = sysr.disturbances; md = length(inpd)
    inpf = sysr.faults; mf = length(inpf)
    inpw = sysr.noise;  mw = length(inpw)
    inpaux = sysr.aux;  maux = length(inpaux)
    mru = length(sysref.controls);
    mru == 0 || mru == mu ||  error("Incompatible control input groups in sysr and sysref")
    mrd = length(sysref.disturbances);
    mrd == 0 || mrd == md ||  error("Incompatible disturbance input groups in sysr and sysref")
    mrf = length(sysref.faults);
    mrf == 0 || mrf == mf ||  error("Incompatible fault input groups in sysr and sysref")
    mrw = length(sysref.noise);
    mrw == 0 || mrw == mw ||  error("Incompatible noise input groups in sysr and sysref")
    mra = length(sysref.aux);
    mra > 0 || mra == maux ||  error("Incompatible auxiliary input groups in sysr and sysref")

    m = mu+md+mf+mw+maux;       # total number of inputs

    rinp = zeros(0,m);
    mru == 0 || (rinp = [rinp; eye(mu,m)])
    mrd == 0 || (rinp = [rinp; zeros(md,mu) eye(md,m-mu)])
    mrf == 0 || (rinp = [rinp; zeros(mf,mu+md) eye(mf,m-mu-md)])
    mrw == 0 || (rinp = [rinp; zeros(mw,mu+md+mf) eye(mw,m-mu-md-mf)])
    mra == 0 || (rinp = [rinp; zeros(maux,m-maux) eye(maux)])
    return isinf(nrmflag) ? ghinfnorm(sysr.sys-sysref.sys*rinp)[1] :
                            gh2norm(sysr.sys-sysref.sys*rinp; atolinf)
end
"""
     γ = fdimmperf(sysr::FDIFilterIF, sysref::FDIFilterIF[, nrmflag])

Compute the model-matching performance `γ` of the fault detection and isolation filter internal form object `sysr::FDIFilterIF`
with respect to the fault detection and isolation reference filter internal form  `sysref::FDIFilterIF`.
If `Ri(λ)` is the transfer function matrix of the `i`-th fault detection and isolation filter internal form
`sysr.sys[i]` and `Mri(λ)` is the transfer function matrix of the `i`-th fault detection and isolation reference filter internal form
`sysref.sys[i]`, then `γ` is a vector whose `i`-th component
`γ[i]` is the  H∞-norm of `Ri(λ)-Mri(λ)`, if `nrmflag = Inf` (default) or the  H2-norm of `Ri(λ)-Mri(λ)`, if `nrmflag = 2`.
The value of `γ[i]` is infinite for an unstable difference `Ri(λ)-Mri(λ)` or if `nrmflag = 2` and the transfer function matrix
`Ri(λ)-Mri(λ)` of a continuous-time system is not strictly proper. In general, `Ri(λ)` and `Mri(λ)` are partitioned as
`Ri(λ) = [ Rui(λ) Rdi(λ) Rfi(λ) Rwi(λ) Rai(λ) ]` and `Mri(λ) = [ Mrui(λ) Mrdi(λ) Mrfi(λ) Mrwi(λ) Mrai(λ) ]` in accordance with
the partitioning of the inputs in control inputs, disturbance inputs, fault inputs, noise inputs and auxiliary inputs.
Void components of `Mri(λ)` corresponding to non-void components in `Ri(λ)` are assumed to be zero.
"""
function fdimmperf(sysr::FDIFilterIF{T1}, sysref::FDIFilterIF{T2}, nrmflag::Real = Inf) where {T1,T2}
    N = length(sysr.sys)
    N == length(sysref.sys) || error("sysr and sysref must have the same number of component filters")
    DescriptorSystems.promote_Ts(sysr.sys[1].Ts, sysref.sys[1].Ts)  # check sampling time

    (nrmflag == Inf || nrmflag == 2) || error("only H∞- and H2-norms supported")
    inpu = sysr.controls; mu = length(inpu)
    inpd = sysr.disturbances; md = length(inpd)
    inpf = sysr.faults; mf = length(inpf)
    inpw = sysr.noise;  mw = length(inpw)
    inpaux = sysr.aux;  maux = length(inpaux)
    mru = length(sysref.controls);
    mru == 0 || mru == mu ||  error("Incompatible control input groups in sysr and sysref")
    mrd = length(sysref.disturbances);
    mrd == 0 || mrd == md ||  error("Incompatible disturbance input groups in sysr and sysref")
    mrf = length(sysref.faults);
    mrf == 0 || mrf == mf ||  error("Incompatible fault input groups in sysr and sysref")
    mrw = length(sysref.noise);
    mrw == 0 || mrw == mw ||  error("Incompatible noise input groups in sysr and sysref")
    mra = length(sysref.aux);
    mra > 0 || mra == maux ||  error("Incompatible auxiliary input groups in sysr and sysref")

    m = mu+md+mf+mw+maux;       # total number of inputs
    rinp = zeros(0,m);
    mru == 0 || (rinp = [rinp; eye(mu,m)])
    mrd == 0 || (rinp = [rinp; zeros(md,mu) eye(md,m-mu)])
    mrf == 0 || (rinp = [rinp; zeros(mf,mu+md) eye(mf,m-mu-md)])
    mrw == 0 || (rinp = [rinp; zeros(mw,mu+md+mf) eye(mw,m-mu-md-mf)])
    mra == 0 || (rinp = [rinp; zeros(maux,m-maux) eye(maux)])

    γ = similar(Array{T1,1},N)
    for i = 1:N
        γ[i] = isinf(nrmflag) ? ghinfnorm(sysr.sys[i]-sysref.sys[i]*rinp)[1] :
                                gh2norm(sysr.sys[i]-sysref.sys[i]*rinp)
    end
    return γ
end

max_stars_repo_name: andreasvarga/FaultDetectionTools

max_issues_repo_name: andreasvarga/FaultDetectionTools

max_forks_repo_name: andreasvarga/FaultDetectionTools


text:
@testset "Test CG poisson example" begin

using HDiscontinuousGalerkin
using Tensors
using SparseArrays

mesh = rectangle_mesh(TriangleCell, (10,10), Vec{2}((0.0,0.0)), Vec{2}((1.0,1.0)))

# ### Initiate function Spaces
dim = 2
Wh = ScalarFunctionSpace(mesh, ContinuousLagrange{dim,RefTetrahedron,1}(); update_face_data = false)

# Declare variables
u_h = TrialFunction(Wh)

# ### Degrees of freedom
dh = DofHandler([u_h],mesh)
K = create_sparsity_pattern(dh);

# ### Boundary conditions
dbc = Dirichlet(u_h, dh, "boundary", x -> 0)

# ### RHS function
f(x::Vec{dim}) = 2*π^2*sin(π*x[1])*sin(π*x[2])

# ### Assembling the linear system
# Now we have all the pieces needed to assemble the linear system, $K u = f$.
function doassemble(Wh, K::SparseMatrixCSC, dh::DofHandler)
    # We allocate the element stiffness matrix and element force vector
    # just once before looping over all the cells instead of allocating
    # them every time in the loop.
    n_basefuncs = getnbasefunctions(Wh)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)
    b = zeros(ndofs(dh))
    cell_dofs = Vector{Int}(undef, ndofs_per_cell(dh))

    # Next we define the global force vector `f` and
    # create an assembler. The assembler
    # is just a thin wrapper around `f` and `K` and some extra storage
    # to make the assembling faster.
    assembler = start_assemble(K, b)

    # It is now time to loop over all the cells in our grid.
    @inbounds for (cellcount, cell) in enumerate(CellIterator(mesh))
        # Always remember to reset the element stiffness matrix and
        # force vector since we reuse them for all elements.
        fill!(Ke, 0)
        fill!(fe, 0)
        reinit!(Wh, cell)

        # It is now time to loop over all the quadrature points in the cell and
        # assemble the contribution to `Ke` and `fe`. The integration weight
        # can be queried from `cellvalues` by `getdetJdV`.
        for q_point in 1:getnquadpoints(Wh)
            dΩ = getdetJdV(Wh, q_point)
            fh = function_value(f, Wh, cell, q_point)
            # For each quadrature point we loop over all the (local) shape functions.
            # We need the value and gradient of the testfunction `v` and also the gradient
            # of the trial function `u`. We get all of these from `cellvalues`.
            for i in 1:n_basefuncs
                v  = shape_value(Wh, q_point, i)
                ∇v = shape_gradient(Wh, q_point, i)
                fe[i] += fh*v * dΩ
                for j in 1:n_basefuncs
                    ∇u = shape_gradient(Wh, q_point, j)
                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                end
            end
        end
        # The last step in the element loop is to assemble `Ke` and `fe`
        # into the global `K` and `f` with `assemble!`.
        celldofs!(cell_dofs, dh, cell)
        assemble!(assembler, cell_dofs, fe, Ke)
    end
    return K, b
end
#md nothing # hide

# ### Solution of the system
# The last step is to solve the system. First we call `doassemble`
# to obtain the global stiffness matrix `K` and force vector `f`.
K, b = doassemble(Wh, K, dh);

apply!(K,b,dbc)
u = K \ b;

reconstruct!(u_h, u, dh)

# ### Compute errors
u_ex(x::Vec{dim}) = sin(π*x[1])*sin(π*x[2])
Etu_h = errornorm(u_h, u_ex)
@test Etu_h <= 0.0002

end


text:
@doc """  download genecode file to data_dir
""" ->
function download_gencode()
    gz_genecode = string(genecode_fl, ".gz")
    if !isdir(data_dir)
        mkdir(data_dir)
    end
    info("Downloading gencode.v19.annotation.gtf now")
    cur_dir = pwd()
    cd(data_dir)
    cmd = `wget -c
    ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz`
    run(cmd)
    gz_genecode = string(genecode_fl, ".gz")
    run(`gunzip $(gz_genecode)`)
    cd(cur_dir)
end

@doc """ read gtf file with out header info
""" ->
function read_gtf_gene()
    if !isfile(genecode_fl)
        download_gencode()
    end
    num_line = countlines(genecode_fl)
    info("gtf total lines $num_line")
    records =  Array{ASCIIString,2}(num_line,9)
    idx = 0
    open(genecode_fl) do file
        while !eof(file)
            line = readline(file)
            if line[1] == '#'
                continue
            end
            fields = split(strip(line,'\n'),"\t")
            if fields[3] != "gene"
                continue
            end
            if contains(fields[9],"gene_type \"pseudogene\"")
                continue
            end
            idx += 1
            idx % 1000_000 == 0 && info("processing $idx line now")
            num_field = length(fields)
            @assert num_field == 9
            records[idx,:] = fields
        end
    end
    records[1:idx,:]
end

@doc """ get gene_name, chr, start, end  given a record
         mgp means mutation genome position
""" ->
function get_gene_mgp_1(fields::Array{ASCIIString,1})
    @assert length(fields) == 9
    chr = convert(ASCIIString, fields[1])
    st  = convert(ASCIIString, fields[4])
    ed  = convert(ASCIIString, fields[5])
    gene_names = filter(subfield->contains(subfield, "gene_name"), split(fields[end],";"))
    @assert length(gene_names) == 1
    gene_name = split(gene_names[1]," ")[end]
    gene_name = split(gene_name,"\"")[2]
    gene_name = convert(ASCIIString, gene_name)
    gene_name,chr,st,ed
end

@doc """ get chr, start, end, gene_name given records
""" ->
function get_gene_mgp_2(data::Array{ASCIIString,2})
    info("begin get_gene_mgp")
    @assert size(data,2) == 9
    n = size(data,1)

    gene_location = Dict{ASCIIString,Tuple{ASCIIString,ASCIIString,ASCIIString}}()
    data = data'
    for i = 1:n
        gene_name,chr,st,ed = get_gene_mgp_1(data[:,i])
        gene_location[gene_name] = (chr,st,ed)
    end
    # gene_chrsted
    save(gene_chrsted_fl, "gene_chrsted_dict", gene_location)

    #(chr,st,ed)=>gene
    location_gene = sort(map(reverse, collect(gene_location)), by=x->x[1])

    #chr=>(st,ed)=>gene_name
    chr_sted_gene = Dict{ASCIIString,Dict{Tuple{UInt64,UInt64},ASCIIString}}()
    last_chr = location_gene[1][1][1]
    tmp_dict = Dict{Tuple{UInt64,UInt64},ASCIIString}()
    for chrsted_gene in location_gene
        chrsted,gene = chrsted_gene
        chr,st,ed    = chrsted
        if chr != last_chr
            chr_sted_gene[last_chr] = copy(tmp_dict)
            empty!(tmp_dict)
            last_chr = chr
        end
        stn = parse(UInt64, st)
        edn = parse(UInt64, ed)
        tmp_dict[(stn,edn)] = gene
    end
    chr_sted_gene[last_chr] = copy(tmp_dict)

    save(chr_sted_gene_fl, "chr_sted_gene_dict", chr_sted_gene)
    info("chr_sted_gene_pair is done")

    nothing
end

@doc """ build index pos_gene dict
""" ->
function build_index_geneloc()
    info("reading gtf file to memory")
    gtf = read_gtf_gene()
    info("Getting gene position from gtf data")
    get_gene_mgp_2(gtf)
    info("built gene_location successfully")
    true
end


text:
using Test
using GRASP
import AtomicLevels: @ro_str, @c_str

@testset "rcsfs.jl" begin

@testset "kappa" begin
    import GRASP: kappa
    @test kappa(ro"1s")  == -1
    @test kappa(ro"2s")  == -1
    @test kappa(ro"2p-") ==  1
    @test kappa(ro"2p")  == -2
    @test kappa(ro"3d-") ==  2
    @test kappa(ro"3d")  == -3
end

import GRASP: CSF, Symmetries
@testset "CSF" begin
    orb1 = ro"1s"  # 1s(1/2)
    orb1b = ro"1s" # 1s(1/2)
    orb2 = ro"3d-" # 3d(3/2)
    orb3 = ro"4d"  # 4d(5/2)

    angmom = Symmetries.AngularMomentum[0, 0]

    # CSF(total2J, parity, orbs, coupled2Js)
    csf1 = CSF([orb1, orb2], [2, 4], angmom, angmom, Symmetries.even)
    @test string(csf1) == "1s(2|0|0) 3d-(4|0|0) | 0+"
    @test nelectrons(csf1) == 6

    csf1b = CSF([orb1b, orb2], [2, 4], angmom, angmom, Symmetries.even)
    @test csf1 == csf1b

    csf2 = CSF([orb1, orb3], [2, 4], angmom, angmom, Symmetries.even)
    @test string(csf2) == "1s(2|0|0) 4d(4|0|0) | 0+"
    @test nelectrons(csf2) == 6

    @test nexcitations(csf1, csf1) == 0
    @test nexcitations(csf2, csf2) == 0
    @test nexcitations(csf1, csf2) == 4

    # Test the ::CSF iterator
    @test length(csf1) == 2
    let cs = collect(csf1), zero_angmom = Symmetries.AngularMomentum(0)
        @test cs[1][1] == orb1
        @test cs[1][2] == 2
        @test cs[1][3] == zero_angmom
        @test cs[1][4] == zero_angmom

        @test cs[2][1] == orb2
        @test cs[2][2] == 4
        @test cs[2][3] == zero_angmom
        @test cs[2][4] == zero_angmom
    end
end

import GRASP: kappa2rso
@testset "kappa2rso" begin
    @test kappa2rso(-1) == "s"
    @test kappa2rso(-2) == "p"
    @test kappa2rso(-3) == "d"
    @test kappa2rso(1)  == "p-"
    @test kappa2rso(2)  == "d-"
end

import GRASP: Symmetries, AngularMomentum, angularmomentum, parity
@testset "parse_rcsf" begin
    let csfbs = GRASP.parse_rcsf(joinpath(@__DIR__, "grasp/csls/example1.c"))
        @test isa(csfbs, Vector{GRASP.CSFBlock})
        @test length(csfbs) == 4

        @test parity(csfbs[1]) == Symmetries.even
        @test angularmomentum(csfbs[1]) == AngularMomentum(0)

        @test angularmomentum(csfbs[2]) == AngularMomentum(1)
        @test angularmomentum(csfbs[3]) == AngularMomentum(2)
        @test angularmomentum(csfbs[4]) == AngularMomentum(3)
    end

    let csfbs = GRASP.parse_rcsf(joinpath(@__DIR__, "grasp/csls/example2.c"))
        @test isa(csfbs, Vector{GRASP.CSFBlock})
        @test length(csfbs) == 4

        @test parity(csfbs[1]) == Symmetries.odd
        @test angularmomentum(csfbs[1]) == AngularMomentum(1//2)

        @test angularmomentum(csfbs[2]) == AngularMomentum(3//2)
        @test angularmomentum(csfbs[3]) == AngularMomentum(5//2)
        @test angularmomentum(csfbs[4]) == AngularMomentum(7//2)

        # Test that we interpret the in-shell and inter-shell couplings correctly
        let csf = csfbs[1].csfs[1]
            #@test string(csf) == "1s(1|1/2|1/2) 2s(1|1/2|1) 2p(3|3/2|1/2) | 1/2-"
            @test csf.terms[1] == 1//2
            @test csf.terms[2] == 1
            @test csf.terms[3] == 1//2
        end
        let csf = csfbs[1].csfs[5]
            #@test string(csf) == "1s(1|1/2|1/2) 2s(1|1/2|1) 2p⁻(2|0|1) 2p(1|3/2|1/2) | 1/2-"
            @test csf.terms[1] == 1//2
            @test csf.terms[2] == 1
            @test csf.terms[3] == 1
            @test csf.terms[4] == 1//2
        end
    end

    let csfbs = GRASP.parse_rcsf(joinpath(@__DIR__, "grasp/csls/example3.c"))
        @test isa(csfbs, Vector{GRASP.CSFBlock})
        @test length(csfbs) == 6

        @test parity(csfbs[1]) == Symmetries.even
        @test angularmomentum(csfbs[1]) == AngularMomentum(1//2)
        @test angularmomentum(csfbs[2]) == AngularMomentum(3//2)
        @test angularmomentum(csfbs[3]) == AngularMomentum(5//2)
        @test angularmomentum(csfbs[4]) == AngularMomentum(7//2)
        @test angularmomentum(csfbs[5]) == AngularMomentum(9//2)
        @test angularmomentum(csfbs[6]) == AngularMomentum(11//2)
    end

    # example4.c generated with:
    #
    #   * ! Orbital order
    #             4  ! Selected core
    #   4f(4,*)
    #   *
    #   4s,4p,4d,4f
    #             0          50  ! Lower and higher 2*J
    #             0  ! Number of excitations
    #   n
    let csfbs = GRASP.parse_rcsf(joinpath(@__DIR__, "grasp/csls/example4.c"))
        @test isa(csfbs, Vector{GRASP.CSFBlock})
        @test length(csfbs) == 11

        @test parity(csfbs[1]) == Symmetries.even
        @test angularmomentum(csfbs[1])  == AngularMomentum(0)
        @test angularmomentum(csfbs[2])  == AngularMomentum(1)
        @test angularmomentum(csfbs[3])  == AngularMomentum(2)
        @test angularmomentum(csfbs[4])  == AngularMomentum(3)
        @test angularmomentum(csfbs[5])  == AngularMomentum(4)
        @test angularmomentum(csfbs[6])  == AngularMomentum(5)
        @test angularmomentum(csfbs[7])  == AngularMomentum(6)
        @test angularmomentum(csfbs[8])  == AngularMomentum(7)
        @test angularmomentum(csfbs[9])  == AngularMomentum(8)
        @test angularmomentum(csfbs[10]) == AngularMomentum(9)
        @test angularmomentum(csfbs[11]) == AngularMomentum(10)

        # The 3rd block (J=2) contains seniority numbers. But we don't actually
        # parse store them at the moment.
        @test length(csfbs[3].csfs) == 17
    end
end

end # @testset "rcsfs.jl"


text:
@testset "FourierNeuralOperator" begin
    m = FourierNeuralOperator()

    𝐱, 𝐲 = get_burgers_data()
    𝐱, 𝐲 = Float32.(𝐱), Float32.(𝐲)
    @test size(m(𝐱)) == size(𝐲)

    loss(𝐱, 𝐲) = sum(abs2, 𝐲 .- m(𝐱)) / size(𝐱)[end]
    data = [(𝐱[:, :, 1:5], 𝐲[:, 1:5])]
    Flux.train!(loss, params(m), data, Flux.ADAM())
end


text:
"""
Main module for `GeophysicalFlows.jl` -- a collection of solvers for geophysical fluid dynamics
problems in periodic domains on CPUs and GPUs. All modules use Fourier-based pseudospectral
methods and leverage the functionality of `FourierFlows.jl` ecosystem.
"""
module GeophysicalFlows

using
  CUDA,
  Statistics,
  SpecialFunctions,
  Reexport,
  DocStringExtensions

@reexport using FourierFlows

include("utils.jl")
include("twodnavierstokes.jl")
include("singlelayerqg.jl")
include("multilayerqg.jl")
include("surfaceqg.jl")
include("barotropicqgql.jl")

@reexport using GeophysicalFlows.TwoDNavierStokes
@reexport using GeophysicalFlows.SingleLayerQG
@reexport using GeophysicalFlows.MultiLayerQG
@reexport using GeophysicalFlows.SurfaceQG
@reexport using GeophysicalFlows.BarotropicQGQL

end # module


text:
#===============================================================================

    PeakIndices.jl

    iterator for peaks in melting curve
    and temperature consistency experiments

    Author: Tom Price
    Date: June 2019

===============================================================================#

import Base: start, next, done, eltype, collect, iteratorsize, SizeUnknown

#===============================================================================
    structs >>
===============================================================================#

struct PeakIndices
    summit_heights  ::Vector{Float_T}
    summit_idc      ::Vector{Int_T}
    nadir_idc       ::Vector{Int_T}
    len_summit_idc  ::Int_T
    len_nadir_idc   ::Int_T
    PeakIndices(h,s,n) = new(vcat(h,0.0), vcat(s,0), vcat(n,0), length(s), length(n))
end

struct PeakIndicesState
    left_nadir_ii   ::Int_T
    summit_ii       ::Int_T
    right_nadir_ii  ::Int_T
end

struct PeakIndicesElement
    left_nadir_idx  ::Int_T
    summit_idx      ::Int_T
    right_nadir_idx ::Int_T
end

#===============================================================================
    methods >>
===============================================================================#

## PeakIndices methods
## iterator functions to find peaks and flanking nadirs

Base.start(iter ::PeakIndices) = PeakIndicesState(0, 0, 0)

Base.done(iter ::PeakIndices, state ::Void) = true
Base.done(iter ::PeakIndices, state ::PeakIndicesState) =
    state.left_nadir_ii > iter.len_nadir_idc

Base.iteratorsize(::PeakIndices) = SizeUnknown()

Base.eltype(iter ::PeakIndices) = Type{PeakIndicesState}

Base.collect(iter ::PeakIndices) =
    [peak for peak in iter if thing(peak)]

## fail if state === nothing
Base.next(iter ::PeakIndices, state ::Void) = (nothing, nothing)

function Base.next(iter ::PeakIndices, state ::PeakIndicesState)
    left_nadir_ii, summit_ii, right_nadir_ii =
        map(f -> getfield(state, f), fieldnames(state))
    ## next summit
    while (summit_ii < iter.len_summit_idc)
        ## summit_ii < iter.len_summit_idc
        ## increment the summit index
        summit_ii += 1
        ## extend nadir range to the right
        while
            (right_nadir_ii < iter.len_nadir_idc)
                right_nadir_ii += 1
                (iter.summit_idc[summit_ii] < iter.nadir_idc[right_nadir_ii]) && break
        end
        ## decrease nadir range to the left, if possible
        while
            (left_nadir_ii < iter.len_nadir_idc) &&
            (iter.nadir_idc[left_nadir_ii + 1] < iter.summit_idc[summit_ii])
                left_nadir_ii += 1
        end
        ## if there is a nadir to the left, break out of loop
        (left_nadir_ii > 0) && break
        ## otherwise try the next summit
    end
    ## fail if no more summits or no flanking nadirs
    if  (summit_ii > iter.len_summit_idc) ||
       !(iter.nadir_idc[left_nadir_ii] < iter.summit_idc[summit_ii] < iter.nadir_idc[right_nadir_ii])
            return (nothing, nothing)
    end
    ## find duplicate summits
    right_summit_ii = summit_ii
    while
        (right_summit_ii < iter.len_summit_idc) &&
        (iter.summit_idc[right_summit_ii + 1] < iter.nadir_idc[right_nadir_ii])
            right_summit_ii += 1
    end
    ## remove duplicate summits by choosing highest summit
    if right_summit_ii > summit_ii
        summit_ii = (iis -> iis[indmax(iter.summit_heights[iis])])(summit_ii:right_summit_ii)
    end
    ## return value
    newstate = PeakIndicesState(left_nadir_ii, summit_ii, right_nadir_ii)
    element  = PeakIndicesElement(iter, newstate)
    if newstate==state
        return (nothing, nothing)
    end
    return (element, PeakIndicesState(left_nadir_ii, right_summit_ii, right_nadir_ii))
end ## next()

## constructor
PeakIndicesElement(iter ::PeakIndices, state ::PeakIndicesState) =
    PeakIndicesElement( iter.nadir_idc[  state.left_nadir_ii  ],
                        iter.summit_idc[ state.summit_ii      ],
                        iter.nadir_idc[  state.right_nadir_ii ] )

max_stars_repo_name: MakerButt/chaipcr

max_issues_repo_name: MakerButt/chaipcr

max_forks_repo_name: MakerButt/chaipcr