text:
################################
# BIBLIOTECA JULIA CalcNumCTEC #
################################

using LinearAlgebra
using Printf
using Base

function Derivada_MDFCentrada(f::Function, x₀::Number, i::Int; h::Number=10^-4, n::Int=9)
    x = x₀-((n-1)/2)*h:h:x₀+((n-1)/2)*h

    A = zeros(n,n)
    for i in range(1,stop=n)
        for j in range(1,stop=n)
            A[i,j] = x[i]^(n-j)
        end
    end

    coefs = A \ f.(x)

    function dif_f_interp(x_in)
        res = 0
        for k in 1:n-i
            res += coefs[k] * ((factorial(n-k))/(factorial(n-k-i))) * x_in^(n-k-i)
        end
        return res
    end

    return dif_f_interp(x₀)
end

# Cálculo de zeros de funções ------------------------------------------------------------------------------

function Zeros_Bissecao(f::Function, a::Number, b::Number;
    tol::Number=10^-12, klim::Number=10^6)
    if f(a) * f(b) < 0
        k = 0
        x = (a + b) / 2

        while abs(f(x)) > tol
            if k == klim
                @assert false "ERRO: O método não convergiu para $klim iterações..."
            else
                k += 1
            end

            x = (a + b) / 2
            if f(x) * f(a) < 0
                b = x
            else
                a = x
            end
        end

        return x
    else
        @assert false "\n\nERRO NA FUNÇÃO ZeroBissecao:\nOs valores da função nos extremos do intervalo dado\nprecisam ter valores de sinais opostos entre si\npara que a convergência do método seja garantida..."
    end
end

function Zeros_Cordas(f::Function, a::Number, b::Number;
    tol::Float64=10^-12, klim::Int64=10^6)
    if f(a) * f(b) < 0
        k = 0
        x = (a * f(b) - b * f(a)) / (f(b) - f(a))

        while abs(f(x)) > tol
            if k == klim
                @assert false "ERRO: O método não convergiu para $klim iterações..."
            else
                k += 1
            end

            x = (a + b) / 2
            if f(x) * f(a) < 0
                b = x
            else
                a = x
            end
        end

        return x
    else
        @assert false "ERRO NA FUNÇÃO ZeroCordas:\nOs valores da função nos extremos do intervalo dado",
        "\nprecisam ter valores de sinais opostos entre si",
        "\npara que a convergência do método seja garantida..."
    end
end

function Zeros_NR(f::Function, x::Number;
    tol::Float64=10^-12, klim::Int64=10^6)
    k = 0
    while abs(f(x)) > tol
        if k == klim
            @assert false "ERRO NA FUNÇÃO ZeroNR:\nO método não convergiu para $klim iterações..."
        else
            k += 1
        end

        if abs(Derivada_MDFCentrada(f, x, 1, n=5)) < tol
            @assert false "ERRO NA FUNÇÃO ZeroNR:\nO ponto x=$x possui derivada nula..."
        end
        x = x - (f(x) / Derivada_MDFCentrada(f, x, 1, n=5))
    end

    return x
end

# Resolução de sistemas de equações lineares --------------------------------------------------------------

function SEL_EliminGauss(A::Matrix, b::Vector; tol::Float64=10^-12)
    (m,n) = size(A)

    if m==n && m == length(b)
        if abs(det(A)) < tol
            @assert false "ERRO na função EliminGauss:\nA matriz é singular..."
        else
            function test(M)
                aux = 0
                for j in range(1, stop = n)
                    for i in range(1, stop = n)
                        if i > j
                            aux += abs(M[i,j])
                        end
                    end
                end
                return aux
            end

            AG = copy(A)
            bG = copy(b)

            for j in 1:n
                i = j
                while A[j,j] == 0
                    i+=1

                    aux_A = AG[j,:]
                    AG[j,:] = AG[i,:]
                    AG[i,:] = aux_A

                    aux_b = b[j]
                    bG[j] = bG[i]
                    bG[i] = aux_b
                    if (i,j) == (m,n)
                        @assert false "ERRO na função EliminGauss:\nA matriz é singular..."
                    end
                end
            end

            i = 1
            j = 1
            k = 0
            while test(AG) > tol
                k += 1
                if i == n
                    j += 1
                    i = j
                end
                i += 1

                bG[i] -= bG[j] * (AG[i,j] / AG[j,j])
                AG[i,:] -= AG[j,:] * (AG[i,j] / AG[j,j])
            end

            x = zeros(n)
            for k in n:-1:1
                x[k] = (bG[k] - dot(AG[k,:], x)) / AG[k,k]
            end

            return x
        end
    else
        @assert false "ERRO NA FUNÇÃO ElimGauss:\nInconsistência nas dimensões da matriz e do vetor..."
    end
end

function SEL_CritConverg(A::Matrix, modo::Char)
    (m,n) = size(A)

    if m == n
        if uppercase(modo) == 'L'
            beta=zeros(n)
            for i in range(1, stop = n)
                for j in range(1, stop = n)
                    beta[i] += abs(A[i,j])
                end
                beta[i] /= abs(A[i,i])
            end
            return max(beta...)
        elseif uppercase(modo) == 'C'
            beta=zeros(n)
            for j in range(1, stop = n)
                for i in range(1, stop = n)
                    beta[j] += abs(A[i,j])
                end
                beta[j] /= abs(A[j,j])
            end
            return max(beta...)
        elseif uppercase(modo) == 'S'
            beta = ones(n)
            for i in range(1, stop = n)
                for j in range(1, stop = n)
                    beta[i] += abs(A[i,j]) * beta[j]
                end
                beta[i] -= 1
                beta[i] /= abs(A[i,i])
            end
            return max(beta...)
        end
    else
        @assert false "\n\nERRO NA FUNÇÃO CritConverg:\nA matriz dada não é quadrada..."
    end
end

function SEL_GaussJacobi(A::Matrix, b::Vector; x_in::Vector=zeros(length(b)),
    tol::Float64=1e-12, klim::Int=10^6)
    # Rseolve o sistema linear A*x=b através de iterações pelo Método da Gauss-Jacobi.
    # Tolerância padrão: 10^(-12)
    # Métrica de erro: eₐ = ||x_out - x_in|| / ||x_out||
    # Retorno da (x_out, err)

    (m,n) = size(A)

    if m==n && m == length(b)
        if abs(det(A)) < tol
            @assert false "ERRO na função GaussJacobi:\nA matriz é singular..."
        else
            for j in 1:n
                i = j
                while A[j,j] == 0
                    i+=1

                    aux_A = A[j,:]
                    A[j,:] = A[i,:]
                    A[i,:] = aux_A

                    aux_b = b[j]
                    b[j] = b[i]
                    b[i] = aux_b
                    if (i,j) == (m,n)
                        @assert false "ERRO na função EliminGauss:\nA matriz é singular..."
                    end
                end
            end

            function err(x1, x2)
                return norm(x2 - x1) / norm(x2)
            end

            k = 1
            x_out = zeros(n)
            for i in range(1, stop = n)
                x_out[i] = (b[i] - dot(A[i,1:i-1], x_in[1:i-1]) - dot(A[i,i+1:n], x_in[i+1:n])) / A[i,i]
            end

            while err(x_in, x_out) > tol
                if k == klim
                    @assert false "ERRO NA FUNÇÃO GaussJacobi:\nO método não convergiu para $klim iterações\ne tolerância de $tol..."
                else
                    k += 1
                end

                x_in = copy(x_out)

                for i in range(1, stop = n)
                    x_out[i] = (b[i] - dot(A[i,1:i-1], x_in[1:i-1]) - dot(A[i,i+1:n], x_in[i+1:n])) / A[i,i]
                end
            end

            return x_out
        end
    else
        @assert false "ERRO NA FUNÇÃO GaussJacobi:\nInconsistência nas dimensões da matriz e do vetor..."
    end
end

function SEL_GaussSeidel(A::Matrix, b::Vector; x_in::Vector=zeros(length(b)),
    tol::Float64=1e-12, klim::Int=10^6)
    # Rseolve o sistema linear A*x=b através de iterações pelo Método da Gauss-Seidel.
    # Tolerância padrão: 10^(-12)
    # Métrica de erro: eₐ = ||x_out - x_in|| / ||x_out||
    # Retorno da (x_out, err)

    (m,n) = size(A)

    if m==n && m == length(b)
        if abs(det(A)) < tol
            @assert false "\n\nERRO na função GaussJacobi:\nA matriz é singular...\n\n"
        else
            for j in 1:n
                i = j
                while A[j,j] == 0
                    i+=1

                    aux_A = A[j,:]
                    A[j,:] = A[i,:]
                    A[i,:] = aux_A

                    aux_b = b[j]
                    b[j] = b[i]
                    b[i] = aux_b
                    if (i,j) == (m,n)
                        @assert false "ERRO na função GaussSeidel:\nA matriz é singular..."
                    end
                end
            end

            function err(x1, x2)
                return norm(x2 - x1) / norm(x2)
            end

            k = 1
            x_out = zeros(n)
            for i in range(1, stop = n)
                x_out[i] = (b[i] - dot(A[i,1:i-1], x_out[1:i-1]) - dot(A[i,i+1:end], x_in[i+1:end])) / A[i,i]
            end

            while err(x_in, x_out) > tol
                if k == klim
                    @assert false "ERRO NA FUNÇÃO GaussSeidel:\nO método não convergiu para $klim iterações\ne tolerância de $tol..."
                else
                    k += 1
                end

                x_in = copy(x_out)

                for i in range(1, stop = n)
                    x_out[i] = (b[i] - dot(A[i,1:i-1], x_out[1:i-1]) - dot(A[i,i+1:end], x_in[i+1:end])) / A[i,i]
                end
            end

            return x_out
        end
    else
        printstyled("ERRO NA FUNÇÃO GaussJacobi:\nInconsistência nas dimensões da matriz e do vetor...\n\n", color=:bold)
    end
end

# Resolução de sistemas de equações não lineares ---------------------------------------------------------

function SENL_NRSecante(F::Function, x_in::Vector;
    tol::Float64=1e-6, klim::Int64=10^3)
    # Rseolve o sistema não linear F(x)=0 através de iterações pelo Método de Newton-Raphson Secante.
    # Ou seja, a Jacobiana é calculada de maneira aproximada em função da tolerância.
    # Tolerância padrão: 10^(-12)

    function J_aprox(x, δ=tol)
        n = length(x)
        Jap = zeros(n,n)

        for j in range(1, stop = n)
            δx = zeros(n)
            δx[j] = δ * x[j]
            for i in range(1, stop = n)
                Jap[i,j] = (F(x + δx)[i] - F(x)[i]) / δx[j]
            end
        end

        return Jap
    end

    if length(F(x_in)) == length(x_in)
        err(x_in, x_out) = norm(x_out - x_in) / norm(x_out)

        k = 1
        x_out = x_in - J_aprox(x_in) \ F(x_in)

        while err(x_in, x_out) > tol
            if k==klim
                @assert false "\n\nERRO na função NR_Secante:\nO método não convergiu para $klim iterações\ne tolerância de $tol...\n\n"
            else
                k+=1
            end

            x_in = copy(x_out)
            x_out = x_in - J_aprox(x_in) \ F(x_in)
        end

        return x_out
    else
        printstyled("\n\nERRO na função NR_Secante:\nOs vetores possuem dimensões diferentes...\n\n", color=:bold)
        return nothing
    end
end

# Interpolação e ajuste ------------------------------------------------------------------------------

function Interpolacao_Polinomial(x::Vector, y::Vector)
    # Calcula e retorna o polinômio interpolador para os dados de entrada via polinômios de Lagrange.
    # O grau do polinômio é determinado pelas dimensões dos vetores x e y.

    n = Int(length(x))

    if n==length(y)
        function f_interpoladora(x_in)
            res = 0
            for i in range(1, stop=n)
                aux = 1
                for j in range(1, stop=n)
                    if i != j
                        aux *= (x_in - x[j])/(x[i] - x[j])
                    end
                end
                res += y[i] * aux
            end
            return res
        end

        return f_interpoladora
    else
        @assert false "ERRO NA FUNÇÃO Interpolacao_Polinomial:\nOs vetores de entrada possuem dimensões diferentes..."
    end
end

function Interpolacao_Polinomial(f::Function, x::Vector)
    # Calcula e retorna o polinômio interpolador para os dados de entrada via polinômios de Lagrange.
    # O grau do polinômio é determinado pela dimensão do vetor x.

    n = Int(length(x))

    function f_interpoladora(x_in)
        res = 0
        for i in range(1, stop=n)
            aux = 1
            for j in range(1, stop=n)
                if i != j
                    aux *= (x_in - x[j])/(x[i] - x[j])
                end
            end
            res += f(x[i]) * aux
        end
        return res
    end

    return f_interpoladora
end

function Interpolacao_Bilinear(x::Vector, y::Vector, z::Vector)
    if length(x) == length(y) == length(z) == 4
        M = [ones(n) x y x.*y]
        b = M \ z
        g(x,y) = dot(b,[1 x y x*y])
        return g
    else
        @assert false "\n\nERRO NA FUNÇÃO Interpolacao_Bilinear:\nOs vetores de entrada precisam ter dimensões de 4 elementos..."
    end
end

function Ajuste_MinimosQuadrados(x::Vector, y::Vector; grau::Int64=1)
    # Calcula e retorna o polinômio ajustado para os dados de entrada via Método dos Mínimos Quadrados.
    # O grau padrão do polinômio é 1.
    # Retorno da função: (polinomio, r²)
    if length(x) != length(y)
        @assert false "ERRO NA FUNÇÃO MinimosQuadrados:\nOs vetores de entrada possuem dimensões diferentes..."
    else
        A = zeros(grau+1, grau+1)
        for i in 0:grau
            for j in 0:grau
                A[i+1,j+1] = sum(x.^(i+j))
            end
        end

        b = zeros(grau+1)
        for i in 0:grau
            b[i+1] = sum(y.*(x.^i))
        end

        coefs = A \ b
    end

    function f_ap(x_in)
        xᵏ = []
        for i in 0:grau
            push!(xᵏ, x_in^i)
        end
        return dot(coefs, xᵏ)
    end

    Sr = sum((y - f_ap.(x)).^2)
    ym = sum(y)/length(y)
    St = sum((y .- ym).^2)
    r² = 1 - (Sr/St)

    return f_ap, r²
end

# Integração Numérica ------------------------------------------------------------------

# OBS: COMPLETAR O ERRO COM A DIFERENCIAÇÃO NUMÉRICA

function Integral_Trapezio(f::Function, x::Tuple; tol::Float64 = 1e-12)
    if length(x) == 2
        a = x[1]
        b = x[2]
        c = (a+b)/2

        T(a,b) = (b-a) * (f(a) + f(b)) / 2

        if T(a,b) - T(a,c) - T(c,b) < tol
            return T(a,b)
        else
            return Integral_Trapezio(f, (a,c), tol=tol) + Integral_Trapezio(f, (c,b),tol=tol)
        end
    else
        @assert false "ERRO NA FUNÇÃO Intergral_Trapezio:\nO intervalo precisa ter apenas 2 valores..."
    end
end

function Integral_Trapezio(x::Vector, y::Vector)
    n = length(x)

    if n == length(y)
        I = 0
        for i in 1:n-1
            I += ((x[i+1]-x[i])/2)*(y[i] + y[i+1])
        end

        #err = abs(-(b-a)^4 * (IntegSimpson(dif_f4, a, b, n=10^3))[1] / (180*n^4))

        return I#, err
    else
        @assert false "ERRO NA FUNÇÃO Integral_Trapezio:\nOs vetores possuem dimensões diferentes..."
    end
end

function IntegralDupla_Trapezio(f::Function, x::Tuple, y::Tuple; tol::Float64=1e-12)
    if length(x) == length(y) == 2
        x₀, x₁ = x
        y₀, y₁ = y
        xₘ = (x₀ + x₁) / 2
        yₘ = (y₀ + y₁) / 2

        function T(x,y)
            return Integral_Trapezio([y...], [Integral_Trapezio([x...], f.([x...], [y[1],y[1]])),Integral_Trapezio([x...], f.([x...], [y[2],y[2]]))])
        end

        if T(x,y) - T((x₀,xₘ),y) - T((xₘ,x₁),y) < tol
            if T(x,y) - T(x,(y₀,yₘ)) - T(x,(yₘ,y₁)) < tol
                println("teste1")
                return T(x,y)
            else
                println("teste2")
                return IntegralDupla_Trapezio(f, x, (y₀,yₘ), tol=tol) + IntegralDupla_Trapezio(f, x, (yₘ,y₁), tol=tol)
            end
        elseif T(x,y) - T(x,(y₀,yₘ)) - T(x,(yₘ,y₁)) < tol
            println("teste3")
            return IntegralDupla_Trapezio(f, (x₀,xₘ), y, tol=tol) + IntegralDupla_Trapezio(f, (xₘ,x₁), y, tol=tol)
        else
            println("teste4")
            return IntegralDupla_Trapezio(f, (x₀,xₘ), (y₀,yₘ), tol=tol) + IntegralDupla_Trapezio(f, (xₘ,x₁), (y₀,yₘ), tol=tol) + IntegralDupla_Trapezio(f, (x₀,xₘ), y, tol=tol) + IntegralDupla_Trapezio(f, (xₘ,x₁), y, tol=tol)
        end
    else
        @assert false "ERRO NA FUNÇÃO IntegralDupla_Trapezio:\nOs intervalos precisam ter apenas 2 valores..."
    end
end

function Integral_Simpson(f::Function, x::Tuple; tol::Float64 = 1e-12)
    if length(x) == 2
        a = x[1]
        b = x[2]
        c = (a+b)/2

        S(a,b) = (b-a) * (f(a) + 4*f((a+b)/2) + f(b)) / 6

        if S(a,b) - S(a,c) - S(c,b) < tol
            return S(a,b)
        else
            return Integral_Simpson(f, (a,c), tol=tol) + Integral_Simpson(f, (c,b), tol=tol)
        end
    else
        @assert false "ERRO NA FUNÇÃO Integral_Simpson:\nO intervalo precisa ter apenas 2 valores..."
    end
end

function Integral_Simpson(x::Vector, y::Vector)
    n = length(x)

    if n == length(y)
        if n%2==1
            I = 0
            for i in 1:2:n-2
                I += ((x[i+2]-x[i])/6)*(y[i] + 4*y[i+1] + y[i+2])
            end

            return I
        else
            @assert false "ERRO NA FUNÇÃO Integral_Simpson:\nO número de elementos dos vetores deve ser ímpar..."
        end
    else
        @assert false "ERRO NA FUNÇÃO IntegSimpson:\nOs vetores possuem dimensões diferentes..."
    end
end

# function IntegralDupla_Simpson(f::Function, x1::Number, x2::Number,
#     y1::Number, y2::Number; nx::Number=(b-a)*10^3, ny::Number=(d-c)*10^3)

#     hx = (x2-x1)/2nx
#     hy = (y2-y1)/2ny

#     M = zeros(2nx+1, 2ny+1)
#     for (j,y) in enumerate(y1:hy:y2)
#         for (i,x) in enumerate(x1:hx:x2)
#             M[i,j] = f(x,y)
#         end
#     end

#     V = zeros(2ny+1)
#     for k in 1:2ny+1
#         V[k] += Integral_Simpson(Vector(x1:hx:x2), M[:,k])
#     end

#     return Integral_Simpson(Vector(y1:hy:y2), V)
# end

# function IntegralDupla_Simpson(f::Function, x₀::Number, x₁::Number, y₀::Number, y₁::Number; tol::Float64=1e-12)
#     S(a,b) = (b-a) * (f(a) + 4*f((a+b)/2) + f(b)) / 6

#     xₘ = (x₀ + x₁) / 2
#     yₘ = (y₀ + y₁) / 2

#     if

#     else

#     end
# end

function IntegralTripla_Simpson(f::Function, x1::Number, x2::Number, y1::Number, y2::Number,
    z1::Number, z2::Number; nx::Number=(x2-x1)*100, ny::Number=(y2-y1)*100, nz::Number=(z2-z1)*100)

    hx = (x2-x1)/2nx
    hy = (y2-y1)/2ny
    hz = (z2-z1)/2nz

    I = []

    M = Matrix{Number}(undef, 2ny+1, 2nz+1)
    for (k,z) in enumerate(z1:hz:z2)
        for (j,y) in enumerate(y1:hy:y2)
            for (i,x) in enumerate(x1:hx:x2)
                M[i,j] = f(x,y,z)
            end
        end
        # println(M)
        V = zeros(2nx+1)
        for k in 1:2nx+1
            V[k] = Integral_Simpson(Vector(x1:hx:x2), M[:,k])
        end
        # println(V)
        push!(I, Integral_Simpson(Vector(y1:hy:y2), V))
    end
    res = Integral_Simpson(Vector(z1:hz:z2), I)
    return res

    # M = Matrix{Number}(undef, 2nx+1, 2ny+1)
    # for (j,y) in enumerate(y1:hy:y2)
    #     for (i,x) in enumerate(x1:hx:x2)
    #         M[i,j] = Integral_Simpson(Vector(x1:hx:x2), S[i,j,:])
    #     end
    # end

    # V = zeros(2ny+1)
    # for k in 1:2ny+1
    #     V[k] = Integral_Simpson(Vector(x1:hx:x2), M[:,k])
    # end

    # return Integral_Simpson(Vector(y1:hy:y2), V)
end

function IntegralTripla_Trapezio(f::Function, x1::Number, x2::Number, y1::Number, y2::Number,
    z1::Number, z2::Number; nx::Int=10^3, ny::Int=10^3, nz::Int=10^3)

    hx = (x2-x1)/nx
    hy = (y2-y1)/ny
    hz = (z2-z1)/nz

    I = []

    M = Matrix{Number}(undef, ny+1, nz+1)
    for (k,z) in enumerate(z1:hz:z2)
        for (j,y) in enumerate(y1:hy:y2)
            for (i,x) in enumerate(x1:hx:x2)
                M[i,j] = f(x,y,z)
            end
        end
        # println(M)
        V = zeros(nx+1)
        for k in 1:nx+1
            V[k] = Integral_Trapezio(Vector(x1:hx:x2), M[:,k])
        end
        # println(V)
        push!(I, Integral_Trapezio(Vector(y1:hy:y2), V))
    end
    return Integral_Trapezio(Vector(z1:hz:z2), I)

    # M = Matrix{Number}(undef, 2nx+1, 2ny+1)
    # for (j,y) in enumerate(y1:hy:y2)
    #     for (i,x) in enumerate(x1:hx:x2)
    #         M[i,j] = Integral_Simpson(Vector(x1:hx:x2), S[i,j,:])
    #     end
    # end

    # V = zeros(2ny+1)
    # for k in 1:2ny+1
    #     V[k] = Integral_Simpson(Vector(x1:hx:x2), M[:,k])
    # end

    # return Integral_Simpson(Vector(y1:hy:y2), V)
end

# Geometria -----------------------------------------------------------------------------------

struct Vertex2d
    x::Number
    y::Number

    function Vertex2d(x::Number, y::Number)

    end

    function Base.print(v::Vertex2d)
        print("(", v.x, ", ", v.y, ")")
    end
    function Base.println(v::Vertex2d)
        println("(", v.x, ", ", v.y, ")")
    end
end

struct Vertex3d
    x::Number
    y::Number
    z::Number

    function Vertex3d(v::Vertex2d,z::Number)
        return Vertex3d(v.x, v.y, z)
    end

    function Base.print(v::Vertex3d)
        print("(", v.x, ", ", v.y, ", ", v.z, ")")
    end
    function Base.println(v::Vertex3d)
        println("(", v.x, ", ", v.y, ", ", v.z, ")")
    end
end

"CalcNumCTEC incluída com sucesso!"

max_stars_repo_name: Lucashcr/CalcNumCTEC

max_issues_repo_name: Lucashcr/CalcNumCTEC

max_forks_repo_name: Lucashcr/CalcNumCTEC


text:
module CUDNN

import CUDAdrv: CUDAdrv, CuContext

using ..CuArrays
using ..CuArrays: libcudnn, active_context, configured

include("libcudnn_types.jl")

const _handles = Dict{CuContext,cudnnHandle_t}()
const _handle = Ref{cudnnHandle_t}(C_NULL)

function handle()
    if _handle[] == C_NULL
        @assert isassigned(active_context) # some other call should have initialized CUDA
        _handle[] = get!(_handles, active_context[]) do
            context = active_context[]
            handle = cudnnCreate()
            atexit(()->CUDAdrv.isvalid(context) && cudnnDestroy(handle))
            handle
        end
    end

    return _handle[]
end

include("error.jl")
include("helpers.jl")
include("libcudnn.jl")
include("nnlib.jl")

function __init__()
    configured || return

    global CUDNN_VERSION = convert(Int, ccall((:cudnnGetVersion,libcudnn), Csize_t, ()))
end

end

ext: jl


text:
function scalsq(A::AbstractArray{T,2}, b::AbstractVector{T}, w = nothing) where {T<:Number}
#
#  Solving the scaled least squares problem
#          W(Ax - b) = 0
#      with iterative refinement.
#
#  input  A --- the matrix
#         b --- the right-hand side vector
#         w --- the scaling vector (diagonal entries of W)
#
#
    m, n = size(A)
	E = real(one(T))

    if w == nothing
        w = ones(T, m)
        for j = 1:m
            abj = abs(b[j])
            if abj > E
                w[j] = E / abj
            end
        end
    end

    for j = 1:m
        A[j,1:end] = A[j,1:end] * w[j]
    end

    b = b .* w

    Q, S = qr(A)

    d = Q' * b
    x = backsub(S, d)

    # one step refinement
    bb = [b; zeros(T, n)]
	B = [I A; A' zeros(T, n, n)]
    r = b - A * x

    #for j = 1:3
    rr = bb - B * [r; x]; #disp(norm(rr));

    s = Q' * rr[1:m]
    c = forsub(S', rr[m+1:m+n])
	c2 = backsub(S, s[1:n] - c)
    #c1 = Q * [c; s[n+1:m]]
    #r = r + c1
    x = x + c2
    #%end;
    x
end


text:
function tg_sum(tg::TreeGradient)
    rv = zero(tg.∇)
    for (_, node) in tg
        rv .+= node.∇
    end
    rv
end

function tg_cmp(a, b)
    key_a, node_a = a
    key_b, node_b = b
    @test key_a == key_b
    @test node_a.∇ ≈ node_b.∇
    @test node_a.isinitialized == node_b.isinitialized
    return
end

∇ = fill(1.0, 2)
n = 10
tg = TreeGradient(∇, n)
ans = collect(tg)
@test length(ans) == 1
tg_cmp(ans[1], (1, TreeGradientNode(fill(1.0, 2), true)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 10

# make a series of insertions, and test that the state is as expected
i, j = 2, 8
v = 2.0
∇i = fill(v, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 3
tg_cmp(ans[1], (1, TreeGradientNode(zeros(2), false)))
tg_cmp(ans[2], (2, TreeGradientNode(fill(v, 2), true)))
tg_cmp(ans[3], (9, TreeGradientNode(zeros(2), false)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 7

i, j = 2, 8
v = 2.1
∇i = fill(v, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 3
tg_cmp(ans[1], (1, TreeGradientNode(zeros(2), false)))
tg_cmp(ans[2], (2, TreeGradientNode(fill(v, 2), true)))
tg_cmp(ans[3], (9, TreeGradientNode(zeros(2), false)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 7

i, j = 2, 10
v = 3.0
∇i = fill(v, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 2
tg_cmp(ans[1], (1, TreeGradientNode(zeros(2), false)))
tg_cmp(ans[2], (2, TreeGradientNode(fill(v, 2), true)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 9

i, j = 4, 6
v = 4.0
∇i = fill(v, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 4
tg_cmp(ans[1], (1, TreeGradientNode(zeros(2), false)))
tg_cmp(ans[2], (2, TreeGradientNode(zeros(2), false)))
tg_cmp(ans[3], (4, TreeGradientNode(fill(v, 2), true)))
tg_cmp(ans[4], (7, TreeGradientNode(zeros(2), false)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 3

i, j = 1, 10
v = 5.0
∇i = fill(v, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 1
tg_cmp(ans[1], (1, TreeGradientNode(fill(v, 2), true)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 10

i, j = 1, 3
∇i = fill(1.0, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
i, j = 4, 4
∇i = fill(3.0, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
i, j = 5, 6
∇i = fill(4.0, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
i, j = 2, 8
∇i = fill(8.0, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 3
tg_cmp(ans[1], (1, TreeGradientNode(zeros(2), false)))
tg_cmp(ans[2], (2, TreeGradientNode(fill(8.0, 2), true)))
tg_cmp(ans[3], (9, TreeGradientNode(zeros(2), false)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 7

i, j = 1, 1
∇j = fill(2.0, 2)
insert!(tg, i, j, ∇j)
∇j .= 0
ans = collect(tg)
@test length(ans) == 3
tg_cmp(ans[1], (1, TreeGradientNode(fill(2.0, 2), true)))
tg_cmp(ans[2], (2, TreeGradientNode(fill(8.0, 2), true)))
tg_cmp(ans[3], (9, TreeGradientNode(zeros(2), false)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 8

i, j = 2, 10
∇i = fill(3.0, 2)
insert!(tg, i, j, ∇i)
∇i .= 0
ans = collect(tg)
@test length(ans) == 2
tg_cmp(ans[1], (1, TreeGradientNode(fill(2.0, 2), true)))
tg_cmp(ans[2], (2, TreeGradientNode(fill(3.0, 2), true)))
@test tg.∇ ≈ tg_sum(tg)
@test tg.ninit == 10

# test that new memory was allocated
∇j .= 1.0
tg_cmp(ans[1], (1, TreeGradientNode(fill(2.0, 2), true)))

# test input validation
i, j = 0, 3
@test_throws ArgumentError insert!(tg, i, j, ∇i)
i, j = 1, n+1
@test_throws ArgumentError insert!(tg, i, j, ∇i)
i, j = 1, 2
∇i = zeros(3)
@test_throws DimensionMismatch insert!(tg, i, j, ∇i)


text:
# ---
# title: 1531. String Compression II
# id: problem1531
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Dynamic Programming
# link: <https://leetcode.com/problems/string-compression-ii/description/>
# hidden: true
# ---
#
# [Run-length encoding](http://en.wikipedia.org/wiki/Run-length_encoding) is a
# string compression method that works by replacing consecutive identical
# characters (repeated 2 or more times) with the concatenation of the character
# and the number marking the count of the characters (length of the run). For
# example, to compress the string `"aabccc"` we replace `"aa"` by `"a2"` and
# replace `"ccc"` by `"c3"`. Thus the compressed string becomes `"a2bc3"`.
#
# Notice that in this problem, we are not adding `'1'` after single characters.
#
# Given a string `s` and an integer `k`. You need to delete **at most**  `k`
# characters from `s` such that the run-length encoded version of `s` has
# minimum length.
#
# Find the _minimum length of the run-length encoded  version of _`s` _after
# deleting at most_`k` _characters_.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "aaabcccd", k = 2
#     Output: 4
#     Explanation: Compressing s without deleting anything will give us "a3bc3d" of length 6. Deleting any of the characters 'a' or 'c' would at most decrease the length of the compressed string to 5, for instance delete 2 'a' then we will have s = "abcccd" which compressed is abc3d. Therefore, the optimal way is to delete 'b' and 'd', then the compressed version of s will be "a3c3" of length 4.
#
# **Example 2:**
#
#
#
#     Input: s = "aabbaa", k = 2
#     Output: 2
#     Explanation: If we delete both 'b' characters, the resulting compressed string would be "a4" of length 2.
#
#
# **Example 3:**
#
#
#
#     Input: s = "aaaaaaaaaaa", k = 0
#     Output: 3
#     Explanation: Since k is zero, we cannot delete anything. The compressed string is "a11" of length 3.
#
#
#
#
# **Constraints:**
#
#   * `1 <= s.length <= 100`
#   * `0 <= k <= s.length`
#   * `s` contains only lowercase English letters.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
module C172

using LinearAlgebra
using StaticArrays
using UnPack
using Unitful

using Flight.Modeling
import Flight.Modeling: init_x0, init_y0, init_u0, init_d0, f_cont!, f_disc!
using Flight.Plotting
import Flight.Plotting: plots
using Flight.Misc

using Flight.Attitude

using Flight.Terrain
using Flight.Airdata
using Flight.Kinematics
using Flight.Dynamics

using Flight.Components
import Flight.Components: MassTrait, WrenchTrait, AngularMomentumTrait, get_mass_properties

using Flight.Propulsion
using Flight.Aerodynamics
using Flight.LandingGear

using Flight.Aircraft
import Flight.Aircraft: assign_joystick_inputs!

using Flight.Input

export C172Aircraft, C172Pwp, C172Ldg

struct C172ID <: AbstractAircraftID end

############################## C172Powerplant ################################

struct C172Pwp <: SystemGroupDescriptor
    left::EThruster
    right::EThruster
end

WrenchTrait(::System{<:C172Pwp}) = HasWrench()
AngularMomentumTrait(::System{<:C172Pwp}) = HasAngularMomentum()

function C172Pwp()

    prop = SimpleProp(kF = 4e-3, J = 0.25)

    left = EThruster(propeller = prop, motor = ElectricMotor(α = CW))
    right = EThruster(propeller = prop, motor = ElectricMotor(α = CCW))

    C172Pwp(left, right)

end

############################ C172LandingGear ############################

struct C172Ldg{L <: LandingGearUnit, R <: LandingGearUnit,
    C <: LandingGearUnit} <: SystemGroupDescriptor
    left::L
    right::R
    center::C
end

WrenchTrait(::System{<:C172Ldg}) = HasWrench()
AngularMomentumTrait(::System{<:C172Ldg}) = HasNoAngularMomentum()

function C172Ldg()

    mlg_damper = SimpleDamper(k_s = 25000, k_d_ext = 1000, k_d_cmp = 1000)
    nlg_damper = SimpleDamper(k_s = 25000, k_d_ext = 1000, k_d_cmp = 1000)

    left = LandingGearUnit(
        strut = Strut(
            t_bs = FrameTransform(r = [-1, -1.25, 1], q = RQuat() ),
            l_0 = 0.0,
            damper = mlg_damper),
        braking = DirectBraking())

    right = LandingGearUnit(
        strut = Strut(
            t_bs = FrameTransform(r = [-1, 1.25, 1], q = RQuat() ),
            l_0 = 0.0,
            damper = mlg_damper),
        braking = DirectBraking())

    center = LandingGearUnit(
        strut = Strut(
            t_bs = FrameTransform(r = [2, 0, 1] , q = RQuat()),
            l_0 = 0.0,
            damper = nlg_damper),
        steering = DirectSteering())

    C172Ldg(left, right, center)

end

############################### C172 Aerodynamics ###########################

C172_aero() = SimpleDrag()

############################## C172Airframe #################################

Base.@kwdef struct C172Airframe{ Aero <: SystemDescriptor, Pwp <: SystemDescriptor,
                        Ldg <: SystemDescriptor} <: SystemGroupDescriptor
    aero::Aero = C172_aero()
    pwp::Pwp = C172_pwp()
    ldg::Ldg = C172_ldg()
end

MassTrait(::System{<:C172Airframe}) = HasMass()
WrenchTrait(::System{<:C172Airframe}) = HasWrench()
AngularMomentumTrait(::System{<:C172Airframe}) = HasAngularMomentum()

function get_mass_properties(::System{<:C172Airframe})

    #for an aircraft implementing a fuel system here we could compute the actual
    #mass properties from the different fuel tank contents

    MassProperties(
        #upreferred(2650u"lb") |> ustrip,
        m = 1202.0197805,
        #upreferred.([948, 1346, 1967]u"m") |> ustrip |> diagm |> SMatrix{3,3,Float64},
        J_Ob_b = SA[948.0 0 0; 0 1346.0 0; 0 0 1967.0],
        r_ObG_b = zeros(SVector{3}))
end

#get_wr_b and get_hr_b use the fallback for SystemGroups, which will in turn
#call get_wr_b and get_hr_b on aero, pwp and ldg

#uses default SystemGroup f_disc! implementation

############################## C172Controls #################################

struct C172Controls <: SystemDescriptor end

Base.@kwdef mutable struct C172ControlsU
    throttle::Bounded{Float64, 0, 1} = 0.0
    yoke_x::Bounded{Float64, -1, 1} = 0.0 #ailerons (+ bank right)
    yoke_y::Bounded{Float64, -1, 1} = 0.0 #elevator (+ pitch up)
    pedals::Bounded{Float64, -1, 1} = 0.0 #rudder and nose wheel (+ yaw right)
    brake_left::Bounded{Float64, 0, 1} = 0.0 #[0, 1]
    brake_right::Bounded{Float64, 0, 1} = 0.0 #[0, 1]
end

#const is essential when declaring type aliases!
struct C172ControlsY
    throttle::Float64
    yoke_x::Float64
    yoke_y::Float64
    pedals::Float64
    brake_left::Float64
    brake_right::Float64
end

init_u0(::C172Controls) = C172ControlsU()

init_y0(::C172Controls) = C172ControlsY(zeros(SVector{6})...)

###################### Continuous update functions ##########################

function assign_component_inputs!(afr::System{<:C172Airframe},
    ctl::System{<:C172Controls})

    @unpack throttle, yoke_x, yoke_y, pedals, brake_left, brake_right = ctl.u
    @unpack aero, pwp, ldg = afr.subsystems

    pwp.u.left.throttle = throttle
    pwp.u.right.throttle = throttle
    ldg.u.center.steering[] = pedals
    ldg.u.left.braking[] = brake_left
    ldg.u.right.braking[] = brake_right

    return nothing
end

function f_cont!(afr::System{<:C172Airframe}, ctl::System{<:C172Controls},
                kin::KinData, air::AirData, trn::AbstractTerrain)

    @unpack aero, pwp, ldg = afr.subsystems

    #could this go in the main aircraft f_cont!?
    assign_component_inputs!(afr, ctl)
    # f_cont!(srf, air) #update surface actuator continuous state & outputs
    f_cont!(ldg, kin, trn) #update landing gear continuous state & outputs
    f_cont!(pwp, kin, air) #update powerplant continuous state & outputs
    f_cont!(aero, air, kin, trn) #requires previous srf update
    # f_cont!(aero, air, kin, srf, trn) #requires previous srf update

    afr.y = (aero = aero.y, pwp = pwp.y, ldg = ldg.y)

end

function f_cont!(ctl::System{<:C172Controls}, ::System{<:C172Airframe},
                ::KinData, ::AirData, ::AbstractTerrain)

    #here, controls do nothing but update their output state. in a more complex
    #aircraft this could even be a complete autopilot implementation
    @unpack throttle, yoke_x, yoke_y, pedals, brake_left, brake_right = ctl.u
    return C172ControlsY(throttle, yoke_x, yoke_y, pedals, brake_left,
    brake_right)
    # return (throttle = 0.0, yoke_x = 0.0, yoke_y = 0.0, pedals = 0.0,
    #         brake_left = 0.0, brake_right = 0.0)

end

######################## Discrete Update Functions ########################

function f_disc!(::System{<:C172Controls}, ::System{<:C172Airframe})
    #this does nothing, but could also be used to implement open loop or closed
    #loop control laws, predefined maneuvers, etc. by calling an additional
    #function we could call as part of
    return false
end

function f_disc!(afr::System{<:C172Airframe}, ::System{<:C172Controls})
    #fall back to the default SystemGroup implementation, the f_disc! for the
    #C172 components don't have to deal with the C172Controls
    return f_disc!(afr)
end

#these should eventually be declared as constants:

C172_kin() = KinLTF()

function C172Aircraft(; kin = C172_kin(), aero = C172_aero(), pwp = C172Pwp(),
                        ldg = C172Ldg())
    AircraftBase(
        C172ID(),
        kin,
        C172Airframe(aero, pwp, ldg),
        C172Controls())
end

######################## Joystick Input Interface ########################

yoke_curve(x) = exp_axis_curve(x, strength = 0.5, deadzone = 0.05)
pedal_curve(x) = exp_axis_curve(x, strength = 1.5, deadzone = 0.1)
brake_curve(x) = exp_axis_curve(x, strength = 0, deadzone = 0.05)

function assign_joystick_inputs!(ac::System{<:AircraftBase{C172ID}}, joystick::XBoxController)

    if get_button_change(joystick, :dpad_up) === button_released
        ac.u.throttle = Float64(ac.u.throttle) + 0.1
    elseif get_button_change(joystick, :dpad_down) === button_released
        ac.u.throttle = Float64(ac.u.throttle) - 0.1
    end

    ac.u.yoke_x = get_axis_data(joystick, :right_analog_x) |> yoke_curve
    ac.u.yoke_y = get_axis_data(joystick, :right_analog_y) |> yoke_curve
    ac.u.pedals = get_axis_data(joystick, :left_analog_x) |> pedal_curve
    ac.u.brake_left = get_axis_data(joystick, :left_trigger) |> brake_curve
    ac.u.brake_right = get_axis_data(joystick, :right_trigger) |> brake_curve
end

end #module


text:
# Simple example of JuMP extension used in the tests to check that JuMP works well with extensions
# The main difference between `JuMP.Model` and `JuMPExtension.MathPrgm` is the fact that in `add_variable` (resp. `add_constraint`),
# `JuMP.Model` applies the modification to its `moi_backend` field while
# `JuMPExtension.MathPrgm` stores the `AbstractVariable` (resp. `AbstractConstraint`) in a list.

using JuMP

function remove!(a, item)
    deleteat!(a, findall(x -> x==item, a))
end

struct ConstraintIndex
    value::Int # Index in `model.constraints`
end

" Mathematical Programm representation "
mutable struct MathPrgm <: JuMP.AbstractModel
    emp
    vars::Vector{Int}
    equs::Vector{Tuple{Int, Bool}}
    matching::Dict{Int,Tuple{Int,Bool}}
    objectivesense::Union{Symbol,MOI.OptimizationSense}
    mps::Vector{MathPrgm}
    equils::Vector{Vector{MathPrgm}}
    solverobj::Ptr{}

    objective_function::Union{JuMP.AbstractJuMPScalar,JuMP._NonlinearExprData}
    obj_dict::Dict{Symbol, Any}                     # Same that JuMP.Model's field `obj_dict`
    function MathPrgm(m)
        mp = new(m, Vector{Int}(), Vector{Tuple{Int, Bool}}(),
                      Dict{Int,Tuple{Int,Bool}}(), MOI.FEASIBILITY_SENSE,
                      Vector{MathPrgm}(), Vector{Vector{MathPrgm}}(),
                      C_NULL, zero(JuMP.GenericAffExpr{Float64, JuMP.VariableRef}),
                      Dict{Symbol, Any}())
        push!(m.mps, mp)
        return mp
    end
end

Base.broadcastable(model::MathPrgm) = Ref(model)

JuMP.object_dictionary(model::MathPrgm) = model.obj_dict

# Variables
#Base.copy(v::MyVariableRef) = v
#Base.copy(v::MyVariableRef, new_model::MathPrgm) = MyVariableRef(new_model, v.idx)

#Base.:(==)(v::MyVariableRef, w::MyVariableRef) = v.model === w.model && v.idx == w.idx
#Base.broadcastable(v::MyVariableRef) = Ref(v)
#JuMP.isequal_canonical(v::MyVariableRef, w::MyVariableRef) = v == w
JuMP.variable_type(::MathPrgm) = JuMP.VariableRef
function JuMP.add_variable(model::MathPrgm, v::JuMP.AbstractVariable, name::String="")
    vref = JuMP.add_variable(model.emp.backend, v, name)
    push!(model.vars, vref.index.value)
    vref
end
function JuMP.add_variable(model::MathPrgm, variable::JuMP.VariableConstrainedOnCreation, name::String)
    var_ref = JuMP.add_variable(model, variable.scalar_variable, name)
    JuMP.add_constraint(model, JuMP.ScalarConstraint(var_ref, variable.set))
    return var_ref
end
function JuMP.add_variable(model::MathPrgm, variable::JuMP.VariablesConstrainedOnCreation, names)
    var_refs = JuMP.add_variable.(model, variable.scalar_variables,
                                  JuMP.vectorize(names, variable.shape))
    JuMP.add_constraint(model, JuMP.VectorConstraint(var_refs, variable.set))
    return JuMP.reshape_vector(var_refs, variable.shape)
end
function JuMP.delete(model::MathPrgm, vref::JuMP.VariableRef)
    @assert JuMP.is_valid(model.emp.backend, vref)
    remove!(model.vars, vref.index.value)
    delete(model.emp.backend, vref)
end
function JuMP.delete(model::MathPrgm, vrefs::Vector{JuMP.VariableRef})
    JuMP.delete.(model, vrefs)
end
JuMP.is_valid(model::MathPrgm, vref::JuMP.VariableRef) = JuMP.is_valid(model.emp.backend, vref)
JuMP.num_variables(model::MathPrgm) = length(model.vars)

# Constraints
JuMP.constraint_type(::MathPrgm) = JuMP.ConstraintRef
function JuMP.add_constraint(model::MathPrgm, c::JuMP.AbstractConstraint,
                             name::String="")
    cref = JuMP.add_constraint(model.emp.backend, c, name)
    push!(model.equs, cref.index.value)
    cref
end
function JuMP.delete(model::MathPrgm, cref::JuMP.ConstraintRef)
    @assert JuMP.is_valid(model, constraint_ref)
    remove!(model.equs, cref.index.value)
    delete(model.emp.backend, cref)
end
function JuMP.delete(model::MathPrgm, con_refs::Vector{<:JuMP.ConstraintRef})
    JuMP.delete.(model, con_refs)
end
function JuMP.is_valid(model::MathPrgm, constraint_ref::JuMP.ConstraintRef)
    return (model === constraint_ref.model &&
            constraint_ref.index in keys(model.constraints))
end
function JuMP.num_constraints(model::MathPrgm,
    F::Type{<:JuMP.AbstractJuMPScalar},
    S::Type{<:MOI.AbstractSet})
    crefs = all_constraints(model.emp.backend, F, S)
    return count(cref -> cref.index.value in model.equs, crefs)
end
function JuMP.num_constraints(model::MathPrgm,
    T::Type{<:Vector{F}},
    S::Type{<:MOI.AbstractSet}) where F<:JuMP.AbstractJuMPScalar
    crefs = all_constraints(model.emp.backend, T, S)
    return count(cref -> cref.index.value in model.equs, crefs)
end

# Objective
function JuMP.set_objective_function(model::MathPrgm, f::JuMP.AbstractJuMPScalar)
    model.objective_function = f
end
function JuMP.set_objective_function(model::MathPrgm, f::Real)
    model.objective_function = JuMP.GenericAffExpr{Float64, JuMP.VariableRef}(f)
end
function JuMP.set_objective_function(model::MathPrgm, f::JuMP._NonlinearExprData)
    model.objective_function = f
end
JuMP.objective_sense(model::MathPrgm) = model.objectivesense
function JuMP.set_objective_sense(model::MathPrgm, sense)
    model.objectivesense = sense
end
JuMP.objective_function_type(model::MathPrgm) = typeof(model.objective_function)
JuMP.objective_function(model::MathPrgm) = model.objective_function
function JuMP.objective_function(model::MathPrgm, FT::Type)
    # InexactError should be thrown, this is needed in `objective.jl`
    if !(model.objective_function isa FT)
        throw(InexactError(:objective_function, FT,
                           typeof(model.objective_function)))
    end
    return model.objective_function::FT
end

# Names
function JuMP.variable_by_name(model::MathPrgm, name::String)
    return variable_by_name(model.emp.backend, name)
end
function JuMP.constraint_by_name(model::MathPrgm, name::String)
    return constraint_by_name(model.emp.backend, name)
end

# Show
function JuMP.show_backend_summary(io::IO, model::MathPrgm) end
function JuMP.show_objective_function_summary(io::IO, model::MathPrgm)
    println(io, "Objective function type: ",
            JuMP.objective_function_type(model))
end
function JuMP.objective_function_string(print_mode, model::MathPrgm)
    return JuMP.function_string(print_mode, JuMP.objective_function(model))
end
_plural(n) = (isone(n) ? "" : "s")
function JuMP.show_constraints_summary(io::IO, model::MathPrgm)
    n = length(model.equs)
    print(io, "Constraint", _plural(n), ": ", n)
end
function JuMP.constraints_string(print_mode, model::MathPrgm)
    strings = String[]
    # Sort by creation order, i.e. ConstraintIndex value
    crefs = all_constraints(model.emp.backend)
    constraints = sort(crefs, by = c -> c.first.index.value)
    for (index, constraint) in constraints
        push!(strings, JuMP.constraint_string(print_mode, constraint))
    end
    return strings
end

JuMP.all_variables(model::MathPrgm) = JuMP.all_variables(model.emp.backend)
#########################################################################
# EMP add
#########################################################################
JuMP._parse_NL_expr_runtime(model::MathPrgm, a::Any, b::Any, c::Any, d::Any) = JuMP._parse_NL_expr_runtime(model.emp.backend, a, b, c, d)


text:
module SDFReader

export header, Header,
    file_summary,
    cached_read,
    labels,
    Stagger, CellCentre, FaceX, FaceY, FaceZ, EdgeX, EdgeY, EdgeZ, Vertex,
    AbstractBlockHeader,
    PlainVariableBlockHeader,
    PointVariableBlockHeader,
    PlainMeshBlockHeader,
    PointMeshBlockHeader

using Unitful

include("constants.jl")
include("sdf_header.jl")
include("read_header.jl")
include("read_data.jl")
include("cache.jl")
include("units.jl")
include("utils.jl")

@generated function typemap(data_type::Val{N}) :: DataType where N
    if N == DATATYPE_NULL
        Nothing
    elseif N == DATATYPE_INTEGER4
        Int32
    elseif N == DATATYPE_INTEGER8
        Int64
    elseif N == DATATYPE_REAL4
        Float32
    elseif N == DATATYPE_REAL8
        Float64
    elseif N == DATATYPE_REAL16
        Nothing # unsupported
    elseif N == DATATYPE_CHARACTER
        Char
    elseif N == DATATYPE_LOGICAL
        Bool
    else
        Nothing
    end
end

file_summary(filename) = open(file_summary, filename)[2]

function file_summary(f::IOStream)
    h = header(f)
    blocks = Dict{Symbol,AbstractBlockHeader}()

    block_start = h.summary_location
    for i in Base.OneTo(h.nblocks)
        block = BlockHeader(f, block_start, h.string_length, h.block_header_length)
        blocks = push!(blocks, Symbol(block.id) => read(f, block))
        block_start = block.next_block_location
    end

    h, blocks
end

function Base.read(f, block::AbstractBlockHeader{T, D}) where {T, D}
    raw_data = read!(f, block)
    apply_normalization!(raw_data, block)

    return add_units(raw_data, block)
end

end # module


text:
precompile(Tuple{typeof(Base.similar), Array{Base.Grisu.Bignums.Bignum, 1}})
precompile(Tuple{typeof(Base.length), Array{Base.Grisu.Bignums.Bignum, 1}})
precompile(Tuple{typeof(Base.deepcopy_internal), Any, Base.IdDict{Any, Any}})
precompile(Tuple{typeof(Base.deepcopy_internal), Array{UInt32, 1}, Base.IdDict{Any, Any}})
precompile(Tuple{typeof(Base.create_expr_cache), String, String, Array{Base.Pair{Base.PkgId, UInt64}, 1}, Base.UUID})
precompile(Tuple{Type{NamedTuple{(:stderr,), T} where T<:Tuple}, Tuple{Base.TTY}})
precompile(Tuple{getfield(Base, Symbol("##pipeline#553")), Nothing, Nothing, Base.TTY, Bool, typeof(Base.pipeline), Base.Cmd})
precompile(Tuple{getfield(Base, Symbol("#kw##pipeline")), NamedTuple{(:stderr,), Tuple{Base.TTY}}, typeof(Base.pipeline), Base.Cmd})
precompile(Tuple{getfield(Base, Symbol("##open#563")), Bool, Bool, typeof(Base.open), Base.CmdRedirect, Base.TTY})
precompile(Tuple{typeof(Base.open), Base.CmdRedirect, String, Base.TTY})
precompile(Tuple{typeof(Base.convert), Type{IO}, Base.PipeEndpoint})
precompile(Tuple{typeof(Base.getproperty), Base.Process, Symbol})
precompile(Tuple{typeof(Base.show), Base.GenericIOBuffer{Array{UInt8, 1}}, Array{String, 1}})
precompile(Tuple{typeof(Base.write), Base.PipeEndpoint, String})
precompile(Tuple{typeof(Base.print), Base.GenericIOBuffer{Array{UInt8, 1}}, Base.UUID})
precompile(Tuple{getfield(Base, Symbol("##sprint#339")), Nothing, Int32, typeof(Base.sprint), Function, UInt64})
precompile(Tuple{typeof(Base.show), Base.GenericIOBuffer{Array{UInt8, 1}}, UInt64})
precompile(Tuple{typeof(Base.write), Base.Process, String})
precompile(Tuple{typeof(Base.show_delim_array), Base.GenericIOBuffer{Array{UInt8, 1}}, Tuple{UInt64, UInt64}, Char, Char, Char, Bool, Int32, Int32})
precompile(Tuple{typeof(Base.print), Base.GenericIOBuffer{Array{UInt8, 1}}, Tuple{UInt64, UInt64}})
precompile(Tuple{getfield(Base, Symbol("##open#271")), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(Base.open), getfield(Base, Symbol("#698#700")), String, Vararg{String, N} where N})
precompile(Tuple{typeof(Requires.__init__)})
precompile(Tuple{typeof(Requires.loadpkg), Base.PkgId})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Requires.listenpkg), Any, Base.PkgId})
precompile(Tuple{typeof(ImageCore.__init__)})
precompile(Tuple{typeof(Base.rehash!), Base.Dict{Base.PkgId, Array{Function, 1}}, Int32})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Base.copyto!), Base.IndexLinear, Array{Function, 1}, Base.IndexLinear, Array{Any, 1}})
precompile(Tuple{typeof(Base.get!), getfield(Requires, Symbol("#1#2")), Base.Dict{Base.PkgId, Array{Function, 1}}, Base.PkgId})
precompile(Tuple{typeof(Base.push!), Array{Function, 1}, Function})
precompile(Tuple{typeof(CompilerSupportLibraries_jll.__init__)})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##load_overrides#3")), Bool, typeof(Pkg.Artifacts.load_overrides)})
precompile(Tuple{typeof(Pkg.Artifacts.process_overrides), Base.Dict{String, Any}, Base.UUID})
precompile(Tuple{typeof(Pkg.Artifacts.do_artifact_str), String, Base.Dict{String, Any}, String, Module})
precompile(Tuple{typeof(Base.getindex), Base.Dict{Symbol, Base.Dict{K, V} where V where K}, Symbol})
precompile(Tuple{typeof(Base.haskey), Base.Dict{Base.UUID, Base.Dict{String, Union{Base.SHA1, String}}}, Base.UUID})
precompile(Tuple{Type{Base.Generator{I, F} where F where I}, getfield(Pkg.Artifacts, Symbol("#21#22")){String, String}, Array{Base.Dict{String, Any}, 1}})
precompile(Tuple{typeof(Pkg.Artifacts.unpack_platform), Base.Dict{String, Any}, String, String})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("#21#22")){String, String}, Base.Dict{String, Any}})
precompile(Tuple{Type{Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}}, Base.Generator{Array{Base.Dict{String, Any}, 1}, getfield(Pkg.Artifacts, Symbol("#21#22")){String, String}}})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("#nosym#16")), String})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("#nover#17")), String})
precompile(Tuple{Type{NamedTuple{(:libgfortran_version, :libstdcxx_version, :cxxstring_abi), T} where T<:Tuple}, Tuple{Base.VersionNumber, Nothing, Nothing}})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), NamedTuple{(:libgfortran_version, :libstdcxx_version, :cxxstring_abi), Tuple{Base.VersionNumber, Nothing, Nothing}}, Type{Pkg.BinaryPlatforms.CompilerABI}})
precompile(Tuple{Type{NamedTuple{(:libc, :compiler_abi), T} where T<:Tuple}, Tuple{Symbol, Pkg.BinaryPlatforms.CompilerABI}})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#Linux#13#14")), Symbol, Nothing, Pkg.BinaryPlatforms.CompilerABI, Type{Pkg.BinaryPlatforms.Linux}, Symbol})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), NamedTuple{(:libc, :compiler_abi), Tuple{Symbol, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.Linux}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.Linux, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{Type{NamedTuple{(:libc, :compiler_abi), T} where T<:Tuple}, Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#FreeBSD#19#20")), Nothing, Nothing, Pkg.BinaryPlatforms.CompilerABI, Type{Pkg.BinaryPlatforms.FreeBSD}, Symbol})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.FreeBSD}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.FreeBSD, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Pkg.BinaryPlatforms.FreeBSD})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.FreeBSD})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#Windows#17#18")), Nothing, Nothing, Pkg.BinaryPlatforms.CompilerABI, Type{Pkg.BinaryPlatforms.Windows}, Symbol})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.Windows}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.Windows, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Pkg.BinaryPlatforms.Windows})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.Windows})
precompile(Tuple{getfield(Core, Symbol("#kw#Type")), NamedTuple{(:libc, :compiler_abi), Tuple{Nothing, Pkg.BinaryPlatforms.CompilerABI}}, Type{Pkg.BinaryPlatforms.MacOS}, Symbol})
precompile(Tuple{Type{Base.Pair{A, B} where B where A}, Pkg.BinaryPlatforms.MacOS, Base.Dict{String, Any}})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Pkg.BinaryPlatforms.MacOS})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Base.Dict{String, Any}, Pkg.BinaryPlatforms.MacOS})
precompile(Tuple{typeof(Base.filter), Function, Base.KeySet{Pkg.BinaryPlatforms.Platform, Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}}})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#39#41")){Pkg.BinaryPlatforms.Linux}, Pkg.BinaryPlatforms.Windows})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#39#41")){Pkg.BinaryPlatforms.Linux}, Pkg.BinaryPlatforms.MacOS})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#39#41")){Pkg.BinaryPlatforms.Linux}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{getfield(Pkg.BinaryPlatforms, Symbol("#39#41")){Pkg.BinaryPlatforms.Linux}, Pkg.BinaryPlatforms.FreeBSD})
precompile(Tuple{typeof(Base.ht_keyindex2!), Base.Dict{Pkg.BinaryPlatforms.Platform, Nothing}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{typeof(Base.setindex!), Base.Dict{Pkg.BinaryPlatforms.Platform, Nothing}, Nothing, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{typeof(Base.push!), Base.Set{Pkg.BinaryPlatforms.Platform}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{typeof(Base.ht_keyindex), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{typeof(Base.getindex), Base.Dict{Pkg.BinaryPlatforms.Platform, Base.Dict{String, Any}}, Pkg.BinaryPlatforms.Linux})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("#kw##ensure_artifact_installed")), NamedTuple{(:platform,), Tuple{Pkg.BinaryPlatforms.Linux}}, typeof(Pkg.Artifacts.ensure_artifact_installed), String, Base.Dict{String, Any}, String})
precompile(Tuple{typeof(Pkg.Artifacts.query_override), Base.SHA1})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##artifact_paths#9")), Bool, typeof(Pkg.Artifacts.artifact_paths), Base.SHA1})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##artifact_exists#11")), Bool, typeof(Pkg.Artifacts.artifact_exists), Base.SHA1})
precompile(Tuple{typeof(Base.convert), Type{Nothing}, Base.SubString{String}})
precompile(Tuple{typeof(Base.convert), Type{Base.Pair{Nothing, String}}, Base.Pair{Base.SubString{String}, String}})
precompile(Tuple{typeof(Base.setindex!), Array{Base.Pair{Nothing, String}, 1}, Base.Pair{Base.SubString{String}, String}, Int32})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Base.Pair{Nothing, String}, 1}, Base.Pair{Base.SubString{String}, String}, Base.Generator{Array{Array{Union{Nothing, Base.SubString{String}}, 1}, 1}, getfield(Pkg.PlatformEngines, Symbol("#88#90"))}, Int32})
precompile(Tuple{typeof(Base.collect), Base.Generator{Array{Array{Union{Nothing, Base.SubString{String}}, 1}, 1}, getfield(Pkg.PlatformEngines, Symbol("#88#90"))}})
precompile(Tuple{typeof(Pkg.PlatformEngines.list_tarball_symlinks), String})
precompile(Tuple{typeof(Base.write), Base.IOStream, Nothing})
precompile(Tuple{getfield(Pkg.PlatformEngines, Symbol("##unpack#91")), Bool, typeof(Pkg.PlatformEngines.unpack), String, String})
precompile(Tuple{getfield(Pkg.PlatformEngines, Symbol("##download_verify_unpack#96")), Nothing, Bool, Bool, Bool, typeof(Pkg.PlatformEngines.download_verify_unpack), String, String, String})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("#kw##artifact_paths")), NamedTuple{(:honor_overrides,), Tuple{Bool}}, typeof(Pkg.Artifacts.artifact_paths), Base.SHA1})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##artifact_path#10")), Bool, typeof(Pkg.Artifacts.artifact_path), Base.SHA1})
precompile(Tuple{typeof(Pkg.Artifacts.create_artifact), getfield(Pkg.Artifacts, Symbol("#39#40")){Bool, String, String}})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##download_artifact#38")), Bool, typeof(Pkg.Artifacts.download_artifact), Base.SHA1, String, String})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##ensure_artifact_installed#42")), Pkg.BinaryPlatforms.Platform, Bool, typeof(Pkg.Artifacts.ensure_artifact_installed), String, Base.Dict{String, Any}, String})
precompile(Tuple{typeof(Pkg.Artifacts.artifact_path), Base.SHA1})
precompile(Tuple{getfield(Pkg.Artifacts, Symbol("##query_override#7")), Base.Dict{Symbol, Base.Dict{K, V} where V where K}, typeof(Pkg.Artifacts.query_override), Base.SHA1})
precompile(Tuple{typeof(Libdl.dlopen), String})
precompile(Tuple{typeof(OpenSpecFun_jll.__init__)})
precompile(Tuple{typeof(Base.foreach), getfield(OpenSpecFun_jll, Symbol("#7#9")), Tuple{Array{String, 1}}})
precompile(Tuple{typeof(Base.foreach), getfield(OpenSpecFun_jll, Symbol("#8#10")), Tuple{Array{String, 1}}})
precompile(Tuple{typeof(VideoIO.loglevel)})
precompile(Tuple{typeof(VideoIO._read_packet), Ptr{VideoIO.AVInput{I} where I}, Ptr{UInt8}, Int32})
precompile(Tuple{typeof(VideoIO.AVFormat.av_find_input_format), String})
precompile(Tuple{typeof(Base.setindex!), VideoIO.AVDict, Int32, String})
precompile(Tuple{typeof(VideoIO.__init__)})
precompile(Tuple{typeof(Glob._match_bracket), Base.SubString{String}, Char, Int32, Char, Char})
precompile(Tuple{typeof(Glob._match), Base.SubString{String}, Int32, Char, Bool, Bool})
precompile(Tuple{Type{Glob.GlobMatch}, String})
precompile(Tuple{typeof(Glob.glob), String, String})
precompile(Tuple{typeof(Base.getindex), Array{Union{Base.SubString{String}, Glob.FilenameMatch{Base.SubString{String}}}, 1}, Int32})
precompile(Tuple{typeof(Base.occursin), Glob.FilenameMatch{Base.SubString{String}}, String})
precompile(Tuple{typeof(Glob._glob!), Array{String, 1}, Glob.FilenameMatch{Base.SubString{String}}})
precompile(Tuple{typeof(VideoIO.AVUtil.av_dict_set), Base.RefValue{Ptr{VideoIO.AVUtil.AVDictionary}}, String, String, Int32})
precompile(Tuple{typeof(Base.stat), Base.Libc.RawFD})
precompile(Tuple{typeof(Base.isvalid_file_crc), Base.IOStream})
precompile(Tuple{typeof(Base.stale_cachefile), String, String})
precompile(Tuple{typeof(Base._include_from_serialized), String, Array{Any, 1}})
precompile(Tuple{typeof(Base._tryrequire_from_serialized), Base.PkgId, UInt64, String})
precompile(Tuple{typeof(Base._require_search_from_serialized), Base.PkgId, String})
precompile(Tuple{typeof(Base.create_expr_cache), String, String, Array{Base.Pair{Base.PkgId, UInt64}, 1}, Nothing})
precompile(Tuple{typeof(Base.compilecache), Base.PkgId, String})
precompile(Tuple{typeof(Base._tryrequire_from_serialized), Base.PkgId, UInt64, Nothing})
precompile(Tuple{typeof(Base._require_from_serialized), String})
precompile(Tuple{typeof(Base._require), Base.PkgId})
precompile(Tuple{typeof(Base.require), Base.PkgId})
precompile(Tuple{typeof(Base.require), Module, Symbol})
precompile(Tuple{typeof(Base.create_expr_cache), String, String, Array{Base.Pair{Base.PkgId, UInt64}, 1}, Base.UUID})
precompile(Tuple{typeof(Base._crc32c), Base.IOStream, UInt32})
precompile(Tuple{getfield(Base, Symbol("##open#271")), Base.Iterators.Pairs{Union{}, Union{}, Tuple{}, NamedTuple{(), Tuple{}}}, typeof(Base.open), getfield(Base, Symbol("#698#700")), String, Vararg{String, N} where N})
precompile(Tuple{typeof(Base.isassigned), Core.SimpleVector, Int32})
precompile(Tuple{typeof(ImageInTerminal.__init__)})
precompile(Tuple{typeof(ImageInTerminal.use_24bit)})
precompile(Tuple{typeof(Core.Compiler.eltype), Type{Array{Array{VideoIO.StreamInfo, 1}, 1}}})
precompile(Tuple{typeof(Base.array_subpadding), Type{T} where T, Type{T} where T})
precompile(Tuple{Type{Base.CyclePadding{P} where P}, DataType})
precompile(Tuple{typeof(Base.padding), Type{T} where T})
precompile(Tuple{typeof(VideoIO.opencamera)})
precompile(Tuple{typeof(Base.collect), Base.Generator{Core.SimpleVector, getfield(VideoIO.AVCodecs, Symbol("#1#2"))}})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Array{VideoIO.StreamInfo, 1}, 1}, Array{VideoIO.StreamInfo, 1}, Base.Generator{Base.UnitRange{Int32}, getfield(VideoIO, Symbol("#2#3"))}, Int32})
precompile(Tuple{typeof(VideoIO.AVFormat.avformat_open_input), Array{Ptr{VideoIO.AVFormat.AVFormatContext}, 1}, String, Ptr{VideoIO.AVFormat.AVInputFormat}, VideoIO.AVDict})
precompile(Tuple{getfield(VideoIO, Symbol("##AVInput#1")), Int32, Type{VideoIO.AVInput{I} where I}, String, Ptr{VideoIO.AVFormat.AVInputFormat}, VideoIO.AVDict})
precompile(Tuple{typeof(Base.divgcd), Int32, Int32})
precompile(Tuple{typeof(Base.collect), Base.Generator{Core.SimpleVector, getfield(VideoIO.AVUtil, Symbol("#1#2"))}})
precompile(Tuple{typeof(VideoIO.AVUtil.fieldposition), Type{VideoIO.AVUtil.AVFrame}, Symbol})
precompile(Tuple{Type{VideoIO.AVUtil.AVFrame}})
precompile(Tuple{getfield(Base, Symbol("#thrownonint#206")){Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, Ptr{UInt8}}, Type{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}}, Type{Ptr{UInt8}}, Int32})
precompile(Tuple{typeof(Base.unalias), Array{Ptr{UInt8}, 1}, Base.ReinterpretArray{Ptr{UInt8}, 1, Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, Array{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, 1}}})
precompile(Tuple{typeof(Base.throw_boundserror), Base.ReinterpretArray{Ptr{UInt8}, 1, Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, Array{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, 1}}, Tuple{Int32}})
precompile(Tuple{typeof(Base.throw_boundserror), Array{Ptr{UInt8}, 1}, Tuple{Base.OneTo{Int32}}})
precompile(Tuple{typeof(Base.copyto!), Array{Ptr{UInt8}, 1}, Base.ReinterpretArray{Ptr{UInt8}, 1, Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, Array{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, 1}}})
precompile(Tuple{getfield(Base, Symbol("#thrownonint#206")){Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Int32}, Type{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}}, Type{Int32}, Int32})
precompile(Tuple{typeof(Base.unalias), Array{Int32, 1}, Base.ReinterpretArray{Int32, 1, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Array{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, 1}}})
precompile(Tuple{typeof(Base.throw_boundserror), Base.ReinterpretArray{Int32, 1, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Array{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, 1}}, Tuple{Int32}})
precompile(Tuple{typeof(Base.copyto!), Array{Int32, 1}, Base.ReinterpretArray{Int32, 1, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Array{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, 1}}})
precompile(Tuple{Type{VideoIO.VideoTranscodeContext}, Ptr{VideoIO.SWScale.SwsContext}, Int32, Int32, Int32, Int32, Int32, Array{UInt8, 3}, Array{VideoIO.AVUtil.AVFrame, 1}, Base.ReinterpretArray{Ptr{UInt8}, 1, Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, Array{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, 1}}, Base.ReinterpretArray{Int32, 1, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Array{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, 1}}})
precompile(Tuple{getfield(VideoIO, Symbol("##VideoReader#4")), Bool, Int32, Int32, Type{VideoIO.VideoReader{transcode} where transcode}, VideoIO.AVInput{String}, Int32})
precompile(Tuple{typeof(VideoIO.opencamera), String, Ptr{VideoIO.AVFormat.AVInputFormat}, VideoIO.AVDict})
precompile(Tuple{getfield(VideoIO.AVCodecs, Symbol("#1#2")), Type{T} where T})
precompile(Tuple{typeof(Base._array_for), Type{Ptr{Nothing}}, Core.SimpleVector, Base.HasLength})
precompile(Tuple{typeof(Base.collect_to!), Array{Ptr{Nothing}, 1}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVCodecs, Symbol("#1#2"))}, Int32, Int32})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Ptr{Nothing}, 1}, Ptr{Nothing}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVCodecs, Symbol("#1#2"))}, Int32})
precompile(Tuple{typeof(Base.copyto!), Array{Any, 1}, Int32, Array{Ptr{Nothing}, 1}, Int32, Int32})
precompile(Tuple{typeof(Base.setindex_widen_up_to), Array{Ptr{Nothing}, 1}, Int64, Int32})
precompile(Tuple{typeof(Base.collect_to!), Array{Any, 1}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVCodecs, Symbol("#1#2"))}, Int32, Int32})
precompile(Tuple{Type{VideoIO.AVCodecs.AVPacket}, Ptr{Nothing}, Int64, Int64, Ptr{Nothing}, Int32, Int32, Int32, Ptr{Nothing}, Int32, Int64, Int64, Int64})
precompile(Tuple{Type{VideoIO.AVInput{String}}, String, Array{Ptr{VideoIO.AVFormat.AVFormatContext}, 1}, Array{Ptr{VideoIO.AVFormat.AVIOContext}, 1}, Int32, Array{VideoIO.AVCodecs.AVPacket, 1}, Array{VideoIO.StreamInfo, 1}, Array{VideoIO.StreamInfo, 1}, Array{VideoIO.StreamInfo, 1}, Array{VideoIO.StreamInfo, 1}, Array{VideoIO.StreamInfo, 1}, Array{VideoIO.StreamInfo, 1}, Base.Set{Int32}, Array{VideoIO.StreamContext, 1}, Bool})
precompile(Tuple{typeof(Base.unsafe_convert), Type{Ref{Ptr{VideoIO.AVUtil.AVDictionary}}}, Base.RefValue{Ptr{VideoIO.AVUtil.AVDictionary}}})
precompile(Tuple{getfield(VideoIO.AVUtil, Symbol("#1#2")), Type{T} where T})
precompile(Tuple{typeof(Base.zero), Type{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}}})
precompile(Tuple{typeof(Base._array_for), Type{Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}}, Core.SimpleVector, Base.HasLength})
precompile(Tuple{typeof(Base.collect_to!), Array{Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, 1}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVUtil, Symbol("#1#2"))}, Int32, Int32})
precompile(Tuple{typeof(Base.collect_to_with_first!), Array{Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, 1}, Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVUtil, Symbol("#1#2"))}, Int32})
precompile(Tuple{typeof(Base.zero), Type{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}}})
precompile(Tuple{typeof(Base.copyto!), Array{Tuple{Any, Any, Any, Any, Any, Any, Any, Any}, 1}, Int32, Array{Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, 1}, Int32, Int32})
precompile(Tuple{typeof(Base.setindex_widen_up_to), Array{Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, 1}, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Int32})
precompile(Tuple{typeof(Base.collect_to!), Array{Tuple{Any, Any, Any, Any, Any, Any, Any, Any}, 1}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVUtil, Symbol("#1#2"))}, Int32, Int32})
precompile(Tuple{typeof(Base.copyto!), Array{Any, 1}, Int32, Array{Tuple{Any, Any, Any, Any, Any, Any, Any, Any}, 1}, Int32, Int32})
precompile(Tuple{typeof(Base.setindex_widen_up_to), Array{Tuple{Any, Any, Any, Any, Any, Any, Any, Any}, 1}, Ptr{Nothing}, Int32})
precompile(Tuple{typeof(Base.collect_to!), Array{Any, 1}, Base.Generator{Core.SimpleVector, getfield(VideoIO.AVUtil, Symbol("#1#2"))}, Int32, Int32})
precompile(Tuple{typeof(Base.zero), Type{VideoIO.AVUtil.AVRational}})
precompile(Tuple{typeof(Base.zero), Type{Tuple{UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64}}})
precompile(Tuple{typeof(Base.zero), Type{Tuple{Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}, Ptr{VideoIO.AVUtil.AVBufferRef}}}})
precompile(Tuple{typeof(Base.findfirst), Base.Fix2{typeof(Base.isequal), Symbol}, Tuple{Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol, Symbol}})
precompile(Tuple{Type{VideoIO.AVUtil.AVFrame}, Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Ptr{Nothing}, Int32, Int32, Int32, Int32, Int32, UInt32, VideoIO.AVUtil.AVRational, Int64, Int64, Int64, Int32, Int32, Int32, Ptr{Nothing}, Tuple{UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64, UInt64}, Int32, Int32, Int32, Int32, Int64, Int32, UInt64, Tuple{Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}}, Ptr{Nothing}, Int32, Ptr{Nothing}, Int32, Int32, UInt32, UInt32, UInt32, UInt32, UInt32, Int64, Int64, Int64, Ptr{Nothing}, Int32, Int32, Int32, Ptr{Nothing}, Int32, Int32, Ptr{Nothing}, Ptr{Nothing}, Ptr{Nothing}, UInt32, UInt32, UInt32, UInt32, Ptr{Nothing}})
precompile(Tuple{Type{VideoIO.VideoReader{true}}, VideoIO.AVInput{String}, VideoIO.StreamInfo, Int32, Ptr{VideoIO.AVCodecs.AVCodecContext}, Ptr{VideoIO.AVCodecs.AVCodec}, Array{VideoIO.AVUtil.AVFrame, 1}, Array{Int32, 1}, Int32, Int32, Int32, Base.Rational{Int32}, Base.Rational{Int32}, Array{Array{UInt8, 1}, 1}, VideoIO.VideoTranscodeContext})
precompile(Tuple{Type{Base.PermutedDimsArrays.PermutedDimsArray{T, N, perm, iperm, AA} where AA<:(AbstractArray{T, N} where N where T) where iperm where perm where N where T}, Array{ColorTypes.Gray{FixedPointNumbers.Normed{UInt8, 8}}, 2}, Tuple{Int32, Int32}})
precompile(Tuple{Type{Base.PermutedDimsArrays.PermutedDimsArray{T, N, perm, iperm, AA} where AA<:(AbstractArray{T, N} where N where T) where iperm where perm where N where T}, Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8, 8}}, 2}, Tuple{Int32, Int32}})
precompile(Tuple{typeof(VideoIO.retrieve), VideoIO.VideoReader{true}})
precompile(Tuple{typeof(Base.read), VideoIO.VideoReader{true}})
precompile(Tuple{typeof(VideoIO.pump), VideoIO.AVInput{String}})
precompile(Tuple{typeof(Base.getindex), Array{VideoIO.StreamContext, 1}, Int32})
precompile(Tuple{typeof(VideoIO.decode_packet), VideoIO.VideoReader{true}, Array{VideoIO.AVCodecs.AVPacket, 1}})
precompile(Tuple{typeof(Base._collect), Base.UnitRange{Int32}, Base.Generator{Base.UnitRange{Int32}, getfield(Base.PermutedDimsArrays, Symbol("#1#2")){(2, 1), (2, 1)}}, Base.EltypeUnknown, Base.HasShape{1}})
precompile(Tuple{Type{Base.PermutedDimsArrays.PermutedDimsArray{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8, 8}}, 2, (2, 1), (2, 1), Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8, 8}}, 2}}}, Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8, 8}}, 2}})
precompile(Tuple{typeof(VideoIO.SWScale.sws_scale), Ptr{VideoIO.SWScale.SwsContext}, Base.ReinterpretArray{Ptr{UInt8}, 1, Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, Array{Tuple{Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}, Ptr{UInt8}}, 1}}, Base.ReinterpretArray{Int32, 1, Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, Array{Tuple{Int32, Int32, Int32, Int32, Int32, Int32, Int32, Int32}, 1}}, Int32, Int32, Array{Ptr{UInt8}, 1}, Array{Int32, 1}})
precompile(Tuple{getfield(Base, Symbol("#throwbits#204")), Type{Array{UInt8, 1}}, Type{Ptr{UInt8}}, Type{Array{UInt8, 1}}})
precompile(Tuple{typeof(Base.reinterpret), Type{Ptr{UInt8}}, Array{Array{UInt8, 1}, 1}})
precompile(Tuple{typeof(VideoIO.retrieve!), VideoIO.VideoReader{true}, Base.PermutedDimsArrays.PermutedDimsArray{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8, 8}}, 2, (2, 1), (2, 1), Array{ColorTypes.RGB{FixedPointNumbers.Normed{UInt8, 8}}, 2}}})
precompile(Tuple{typeof(Base.empty!), Base.Dict{Int32, Nothing}})
precompile(Tuple{typeof(VideoIO.fieldposition), Type{VideoIO.AVFormat.AVIOContext}, Symbol})
precompile(Tuple{typeof(Base.close), VideoIO.AVInput{String}})
precompile(Tuple{typeof(VideoIO._close), VideoIO.VideoReader{true}})


text:
using Sockets
using Printf
using Dates
#####################              ########################
function writeInt16( sock,num::UInt16)
    # println( num )
    a1=UInt8(num & 0x00FF)
    a2 =UInt8( num>>8)
    write( sock, a1,a2 )
    return 2
end
#####################################################################
function SendHeadID(sock)
    writeInt16( sock ,0xEB90)
    return writeInt16( sock ,0xEB90)
end

function readHeadID(sock) :: UInt32
    return read( sock ,UInt32)
end
# unsigned short type ; /*01表示原始数据，02表示滤波后数据，03表示过零数据，04表示能量数据    */
function sendTypeCode(sock ,type::UInt16)
    return writeInt16( sock ,type)
    # unsigned short checkCode ;
end
function readTypeCode(sock ) :: UInt16
    return read( sock ,UInt16)
    # unsigned short checkCode ;
end

function sendcheckCode( sock )
    return writeInt16( sock ,0xEB90)
end
function readcheckCode( sock )  :: UInt16
    return read( sock ,UInt16)
end
# char machineName[20];// null terminated string，固定长度
function sendmachineName( sock,machineName::String)
    len = length( machineName  )
    if len >= 20
        temstr = machineName[1:20]
        return write( sock ,temstr)
    else
        temstr = @sprintf( "%s",machineName)
        # println( "--",temstr,"--")
        write( sock ,temstr)
        for i=len+1:20
            write(sock,UInt8(0x00))
        end
        return 20
    end

end

function readmachineName( sock ) :: String
    machinename =""
    terminated = false
    # position( )
    for i=1:20
        c= read( sock ,UInt8)
        # print(string(c))
        if c != 0x00 && terminated == false
            # print(string("-",c,"-"))
            c= Char(c)
            machinename = string( machinename, c)
        else
            terminated = true
            continue
        end
    end
    machinename =strip( machinename)
    # println(" machinename is: --$machinename----")
    return machinename
end

function sendDT( sock ,DT::String)
    len = length( DT  )
    if len > 14
        temstr = DT[1:14]
    else
        temstr = @sprintf( "%-14s",DT)
    end
    # println( "--------------",temstr)
    return write( sock ,temstr)
    # return writeInt16( sock ,0xEB90)
end

function readdDT( sock)::String
    DT  =""
    for i=1:14
        c= read( sock ,Char)
        # if c != 0x00
        DT = string(DT,c)
        # else
        #     break
        # end
    end
    return DT
end

function sendparas(sock ,CfgData::Array{UInt16,1})
    for i in CfgData
        writeInt16(sock,i)
    end
end

function readCfg(fi1) :: Array{UInt16,1}
    cfgs=Array{UInt16,1}(undef,7)
    for idx =1: length( cfgs )
        cfgs[idx] = read( fi1,UInt16)
    end
    return cfgs
end

function senddataType(sock )
    # unsigned short dataType;
    # //01表示char，02表示uchar，03表示int16，04表示uint16，05表示int32，06表示uint32，07表示float，08表示double
    writeInt16(sock,0x0004)
end

function readdataType(sock )
    # unsigned short dataType;
    # //01表示char，02表示uchar，03表示int16，04表示uint16，05表示int32，06表示uint32，07表示float，08表示double
    return read(sock,UInt16)
end

function sendDataPayLoad(sock ,data)
    for i =1:length( data )
        writeInt16(sock,data[i] )
    end

end

function sendDataPayLoad(sock ,frameLen::UInt,data::Array{UInt16,1})
    Ar = reinterpret(UInt8,data)
    # println( "------",length(Ar) ," ",frameLen )
    a=0x0000
    for i =1:frameLen
       yield()
       a+=write( sock , Ar)
    end
    # println(a)
end

function readindatas(fi1,dataBuf::Array{UInt16,1} ) :: Array{UInt16,1}
    cfgs=dataBuf #Array{UInt16,1}(undef,8192)
    for idx =1: length( cfgs )
        cfgs[idx] = read( fi1,UInt16)
    end
    return cfgs
end

function readindatas(fi1,frameLen,dataBuf::Array{UInt16,1}) #:: Array{UInt16,1}
    cfgs =@view dataBuf[ 1:8192*frameLen ]
    for idx =1: length( cfgs )
        cfgs[idx] = read( fi1,UInt16)
    end
    # return cfgs
end

function readindatas!(fi1,frameLen,dataBuf::Array{UInt16,1}) #:: Array{UInt16,1}
    cfgs =@view dataBuf[ 1:8192* frameLen]
    dr = reinterpret(UInt8, cfgs)
    unsafe_read( fi1 ,pointer(dr ), length(dr ))
end

function  sendtailCheckCode(sock )
    # unsigned int tailCheckCode; //FFFF FFFF
    writeInt16(sock,0xFFFF)
    writeInt16(sock,0xFFFF)
end
function  readtailCheckCode(sock )
    return read(sock,UInt32)
end

######################

function sendTestFrames(sock,fi1)
    cfgdata =readCfg(fi1)
    dataBuf = Array{UInt16,1}(undef,8192 )
    dataBuf.=0x0000
    while ! eof(fi1)

        pos1 = position( fi1 )
        sendTestFrame(sock,fi1,UInt16(1),cfgdata,dataBuf  )

        seek( fi1,pos1 )
        sendTestFrame(sock,fi1,UInt16(2),cfgdata, dataBuf )

        seek( fi1,pos1 )
        sendTestFrame(sock,fi1,UInt16(3),cfgdata , dataBuf  )

        seek( fi1,pos1 )
        sendTestFrame(sock,fi1,UInt16(4),cfgdata , dataBuf  )
        # println(pos1)
    end
end

function sendTestFrame(sock,fi1,typeCode::UInt16,cfgdata::Array{UInt16,1} ,dataBuf::Array{UInt16,1} )

    SendHeadID(sock)
    sendTypeCode(sock,typeCode)

    sendcheckCode(sock)
    a= sendmachineName(sock,"TEST-567890")

    DT = readdDT(fi1)
    a= sendDT(sock,DT)
    sendparas(sock,cfgdata)

    readdataType(fi1)
    senddataType(sock)

    readindatas(fi1,1,dataBuf)
    # Datas =@view dataBuf[1:8192]
    if typeCode ==1 || typeCode == 2
        # for i=1:512
        @time  sendDataPayLoad(sock,UInt(512),dataBuf)
        # end
    else
        @time  sendDataPayLoad(sock,UInt(1),dataBuf)
    end
    # sendDataPayLoad(sock,Datas)

    sendtailCheckCode(sock)
end

# Base.include(Main,"ClientTsk.jl")
# Base.include(Main,"SvrTsk.jl")

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
function reorder_bytes(bvec::T, elsize) where T<:AbstractVector{I} where I<:Union{UInt8,Int8}
    l = length(bvec)
    h = Int(l//elsize)
    m = reshape(bvec, elsize, h)
    m1 = m[end:-1:begin, :]  # having @views result in a slow collect with a lot of allocations downstream
    return reshape(m1, l)
end

function reorder_bytes(vc::T) where T<:AbstractArray{N} where N<:Number
    elsize = sizeof(eltype(vc))
    l = length(vc)
    v = vec(vc)
    s = elsize * l
    bvec = reinterpret(UInt8, v)
    return reorder_bytes(bvec, elsize)
end


text:
"""
# module Substitution

# Examples

```jldoctest
julia>
```
"""
module Substitution

using LinearAlgebra: rank, UpperTriangular, LowerTriangular

using APMAE4001.Matrix: issquare

export Forward,
    Backward,
    substitution

struct Forward end
struct Backward end

function substitution(::Backward, A::AbstractMatrix, b::AbstractVector)
    issquare(A) || throw(DimensionMismatch("The matrix `A` is not a square matrix!"))
    m = size(A, 1)
    solutions = zeros(m)
    solutions[m] = b[m] / A[m, m]
    for i in (m - 1):-1:1
        solutions[i] = (b[i] - sum(A[i, j] * solutions[j] for j in (i + 1):m)) / A[i, i]
    end
    solutions
end
function substitution(::Forward, A::AbstractMatrix, b::AbstractVector)
    issquare(A) || throw(DimensionMismatch("The matrix `A` is not a square matrix!"))
    m = size(A, 1)
    solutions = zeros(m)
    solutions[1] = b[1] / A[1, 1]
    for i in 2:m
        solutions[i] = (b[i] - sum(A[i, j] * solutions[j] for j in 1:(i - 1))) / A[i, i]
    end
    solutions
end
substitution(A::UpperTriangular, b::AbstractVector) = substitution(Backward(), A.data, b)
substitution(A::LowerTriangular, b::AbstractVector) = substitution(Forward(), A.data, b)

end


text:
function parsubset(differ_group, start_val)
    differ_group = hcat(differ_group...)
    npar, ngroup = size(differ_group)

    parunique = Vector{Int64}(undef, npar)
    for i in 1:npar
        parunique[i] = length(unique(differ_group[i, :]))
    end

    npar_effective = sum(parunique)

    locations = zeros(Int64, ngroup, npar)
    for i in 1:ngroup
        locations[i, :] .=
            cumsum(parunique) .- parunique .+ differ_group[:, i]
    end

    locations = [locations[i, :] for i in 1:size(locations, 1)]

    start_val_long = vcat(fill.(start_val, parunique)...)

    # parsubsets = falses(ngroup, npar_effective)
    # parsubsets = [parsubsets[i, :] for i in 1:size(parsubsets, 1)]

    # for i in 1:ngroup
    #     parsubsets[i][locations[i, :]] .= true
    # end

    return locations, start_val_long
end

function MGSem(
    data, start_val, differ_group, rowind,
    observed, imply, loss, diff;
    obs_args = (),
    imply_agrs = (),
    loss_args = (),
    diff_args = ())

    #diff = diff(diff_args...)
    #imply = imply(imply_args...)
    #loss = loss(loss_args...)
    #observed = observed(obs_args...)

    sem_vec = semvec(observed, imply, loss, diff)

    par_subsets, start_val_long = parsubset(differ_group, start_val)

    return MGSem(sem_vec, par_subsets), start_val_long
end

function (semmg::MGSem)(par)
    F = zero(eltype(par))
    for i in 1:length(semmg.sem_vec)
        F += (semmg.sem_vec[i].observed.n_obs - 1)*
            semmg.sem_vec[i](view(par, semmg.par_subsets[i]))
    end
    return F
end


text:
module AtomicLevels

using UnicodeFun
using Formatting
using Parameters
using BlockBandedMatrices
using WignerSymbols
macro hi_str(s)
    parse(HalfInteger, s)
end
using Combinatorics

if VERSION < v"1.1-DEV"
    isnothing(::Nothing) = true
    isnothing(::Any) = false
end

include("common.jl")
include("parity.jl")
include("orbitals.jl")
include("configurations.jl")
include("excited_configurations.jl")
include("terms.jl")
include("allchoices.jl")
include("jj_terms.jl")
include("couple_terms.jl")
include("csfs.jl")
include("jj2lsj.jl")
include("levels.jl")

module Utils
include("utils/print_states.jl")
end

end # module


text:
# Tests for Atomistic.jl

using Atomistic

using AtomsBase
using InteratomicPotentials
using Molly
using NBodySimulator
using StaticArrays
using Test
using Unitful

import Atomistic: ElementMassBody, AugmentedAtomData, UnimplementedError

Base.:(==)(a1::AtomsBase.Atom, a2::AtomsBase.Atom) = propertynames(a1) == propertynames(a2) && all(getproperty(a1, p) == getproperty(a2, p) for p ∈ propertynames(a1))
Base.:(==)(b1::ElementMassBody, b2::ElementMassBody) = all(getproperty(b1, p) == getproperty(b2, p) for p ∈ propertynames(b1))
Base.:(==)(a1::AugmentedAtomData, a2::AugmentedAtomData) = propertynames(a1) == propertynames(a2) && all(getproperty(a1, p) == getproperty(a2, p) for p ∈ propertynames(a1))

@testset "Atomistic.jl" begin
    include("unit_convention.jl")
    include("initialization.jl")
    @testset "Atomistic API Unit Tests" begin
        include("api/exceptions.jl")
        include("api/md_result.jl")
        include("api/md_simulator.jl")
    end
    @testset "Implementations Unit Tests" begin
        @testset "NBodySimulator Implementation" begin
            include("implementations/nbodysimulator/nbodysimulator_base.jl")
            include("implementations/nbodysimulator/nbodysimulator_simulator.jl")
            include("implementations/nbodysimulator/nbodysimulator_result.jl")
        end
        @testset "Molly Implementation" begin
            include("implementations/molly/molly_base.jl")
            include("implementations/molly/molly_simulator.jl")
            include("implementations/molly/molly_result.jl")
        end
    end
    @testset "Analysis Unit Tests" begin
        include("analysis/plotting.jl")
    end
end


text:
# This is an example script demonstrating how PARDISO works on a
# medium-sized Hermitian positive definite matrix.
using Pardiso
using LinearAlgebra # for norm
using SparseArrays
using Random

# Script parameters.
# -----------------
verbose = false
n       = 100
lambda  = 3

# Create the Hermitian positive definite matrix A and the vector b in the
# linear system Ax = b.
e = ones(n)
e2 = ones(n-1)
A = spdiagm(-1 => im*e2, 0 => lambda*e, 1 => -im*e2)
b = rand(n) + im * zeros(n)

# Initialize the PARDISO internal data structures.
# ps = PardisoSolver()
ps = MKLPardisoSolver()

if verbose
    set_msglvl!(ps, Pardiso.MESSAGE_LEVEL_ON)
end

# If we want, we could just solve the system right now.
# Pardiso.jl will automatically detect the correct matrix type,
# solve the system and free the data
X1 = solve(ps, A, b)

# First set the matrix type to handle general complex
# hermitian positive definite matrices
set_matrixtype!(ps, Pardiso.COMPLEX_HERM_POSDEF)

# Initialize the default settings with the current matrix type
pardisoinit(ps)

# Get the correct matrix to be sent into the pardiso function.
# :N for normal matrix, :T for transpose, :C for conjugate
A_pardiso = get_matrix(ps, A, :N)

# Analyze the matrix and compute a symbolic factorization.
set_phase!(ps, Pardiso.ANALYSIS)
pardiso(ps, A_pardiso, b)
@printf("The factors have %d nonzero entries.\n", get_iparm(ps, 18))

# Compute the numeric factorization.
set_phase!(ps, Pardiso.NUM_FACT)
pardiso(ps, A_pardiso, b)

# Compute the solutions X using the symbolic factorization.
set_phase!(ps, Pardiso.SOLVE_ITERATIVE_REFINE)
x = similar(b) # Solution is stored in X
pardiso(ps, x, A_pardiso, b)
@printf("PARDISO performed %d iterative refinement steps.\n", get_iparm(ps, 7))

# Compute the residuals.
r = abs.(A*x - b)
@printf("The maximum residual for the solution is %0.3g.\n",maximum(r))
@test norm(R) < 1e-10

# Free the PARDISO data structures.
set_phase!(ps, Pardiso.RELEASE_ALL)
pardiso(ps)


text:
# TODO: tests with multiD variables

Random.seed!(1234)

p = 2
P1 = randn(p)
P1 *= P1'
mu = zeros(p)
mvg = PDSampler.MvGaussianCanon(mu,P1)
chain = chaingraph([ Factor( (x,v)->PDSampler.nextevent_bps(mvg, x, v),
                              x   ->PDSampler.gradloglik(mvg, x),
                              i) for i ∈ 1:3])

Random.seed!(1234)

# test_local_definegraph defines a chain.
nvars = chain.structure.nvars
x0 = randn(nvars)
v0 = randn(nvars)
v0 /= norm(v0)
T = abs(randn())

### TEST LOCALSIMULATION
v0a = randn(nvars+1)
x0a = randn(nvars-1)
@test_throws AssertionError LocalSimulation(chain,x0,v0a,T, 1000, 0.01)
@test_throws AssertionError LocalSimulation(chain,x0a,v0,T, 1000, 0.01)
@test_throws AssertionError LocalSimulation(chain,x0,v0,-0.1,1000,0.01)
@test_throws AssertionError LocalSimulation(chain,x0,v0,T,-10,0.01)
@test_throws AssertionError LocalSimulation(chain,x0,v0,T,1000,-0.1)

ls = LocalSimulation(chain, x0, v0, T, 1000, 0.001)
@test ls.x0 == x0 && ls.v0 == v0 && ls.T == T && ls.maxnevents == 1000 &&
        ls.lambdaref == 0.001 && ls.fg.structure.flist == chain.structure.flist

ls2 = LocalSimulation(
        factorgraph=chain,
        x0=x0, v0=v0, T=T,
        maxnevents=1000, lambdaref=0.001)
@test ls2.x0 == x0 &&
        ls2.v0 == v0 &&
        ls2.T == T &&
        ls2.maxnevents == 1000 &&
        ls2.lambdaref == 0.001 &&
        ls2.fg.structure.flist == chain.structure.flist

### TEST LS_INIT
Random.seed!(140)
(start, all_evlist, pq, tref) = PDSampler.ls_init(ls2)

@test time()-start > 0.0 # and expected to be small...
i = rand(1:chain.structure.nvars)
@test   all_evlist.evl[i].xs[1] == x0[i] &&
        all_evlist.evl[i].vs[1] == v0[i] &&
        all_evlist.evl[i].ts[1] == 0.0
Random.seed!(140)
@test tref == randexp()/0.001

# pq for i <> NOTE it's hard to test precisely
# in parallel setting because would need to know
# how many "rand()" have been called before actual call..
@test pq[i] > 0.0

### TEST LS_RESHAPE
v = Vector{PDSampler.AllowedVarType}(undef, 5)
v = [randn(), randn(7), randn(), randn(50), randn(10)]
w = Vector{PDSampler.AllowedVarType}(undef, 5)
for i in 1:length(v)
    w[i] = v[i]*0.
end
u = vcat(v...)

@test PDSampler.ls_reshape(u, w) == v

all_evlist = AllEventList(Float64, nvars)

for i in 1:nvars
    evi = Event(x0[i], v0[i], 0.0)
    pushevent!(all_evlist.evl[i], evi)
end

i   = rand(1:nvars)
j   = rand(1:nvars)
k   = rand(1:nvars)
t   = rand()
t1  = t+rand()
t2  = t1+rand()
t3  = t2+rand()
ev1 = Event(randn(), randn(), t1)
ev2 = Event(randn(), randn(), t2)
ev3 = Event(randn(), randn(), t3)
pushevent!(all_evlist.evl[i], ev1)
pushevent!(all_evlist.evl[j], ev2)
pushevent!(all_evlist.evl[k], ev3)

### Testing LS_RETRIEVE (1D) (without reflection)
fidx = 2
(xf, vf, g, vars) = PDSampler.ls_retrieve(chain, fidx, all_evlist, 0.0)
@test   xf == x0[assocvariables(chain, fidx)] &&
        vf == v0[assocvariables(chain, fidx)] &&
        g  == PDSampler.gradloglik(mvg, vcat(xf...))

fidx = 1
t    = t3+rand()
(xf, vf, g, vars) = PDSampler.ls_retrieve(chain, fidx, all_evlist, t, false)
aev1x = all_evlist.evl[1].xs[end]
aev1t = all_evlist.evl[1].ts[end]
aev1v = all_evlist.evl[1].vs[end]
aev2x = all_evlist.evl[2].xs[end]
aev2t = all_evlist.evl[2].ts[end]
aev2v = all_evlist.evl[2].vs[end]

@test   vars == [1,2] &&
        isapprox(vcat(xf...), [aev1x+(t-aev1t)*aev1v; aev2x+(t-aev2t)*aev2v]) &&
        isapprox(vcat(vf...), [aev1v; aev2v]) &&
        isapprox(g, chain.factors[fidx].gll(vcat(xf...)))

### Testing LS_RETRIEVE (1D) (with reflection)
fidx = 2
(xf, vf, g, vars) = PDSampler.ls_retrieve(chain, fidx, all_evlist, t, true)
aev1x = all_evlist.evl[2].xs[end]
aev1t = all_evlist.evl[2].ts[end]
aev1v = all_evlist.evl[2].vs[end]
aev2x = all_evlist.evl[3].xs[end]
aev2t = all_evlist.evl[3].ts[end]
aev2v = all_evlist.evl[3].vs[end]

@test   vars == [2,3] &&
        isapprox(vcat(xf...), [aev1x+(t-aev1t)*aev1v; aev2x+(t-aev2t)*aev2v]) &&
        isapprox(vcat(vf...), reflect_bps!(g, vcat(aev1v, aev2v))) &&
        isapprox(g, chain.factors[fidx].gll(vcat(xf...)))

### Testing LS_SAVEUPDATE!
PDSampler.ls_saveupdate!(all_evlist, vars, xf, vf, t)
ii = rand(1:length(vars))
vi = vars[ii]

@test getlastevent(all_evlist.evl[vi]) == Event(xf[ii], vf[ii], t)

### Testing LS_UPDATEPQ!
pq = PDSampler.PriorityQueue()
Random.seed!(123); PDSampler.ls_updatepq!(pq, chain, fidx, xf, vf, g, t)
# It's gaussian so no need for thinning
Random.seed!(123); bounce = chain.factors[fidx].nextevent(vcat(xf...), vcat(vf...))

@test pq[fidx] == t .+ bounce.tau

### Testing LS_RANDOM
evl1 = EventList(Float64)
push!(evl1.xs, randn())
push!(evl1.vs, randn())
push!(evl1.ts, abs(randn()))

evl2 = EventList(Vector{Float64})
push!(evl2.xs, randn(2))
push!(evl2.vs, randn(2))
push!(evl2.ts, abs(randn()))

Random.seed!(321); a1 = PDSampler.ls_random(evl1)
Random.seed!(321); a2 = PDSampler.ls_random(evl2)
Random.seed!(321); aa = randn()
Random.seed!(321); ab = randn(2)

@test a1==aa && a2==ab

### Testing LS_REFRESHMENT
t  = randexp()

Random.seed!(123); pq = PDSampler.ls_refreshment(chain, t, all_evlist)
Random.seed!(123)
### this is a bit of a silly test but ls_refreshment is the
# composite of functions that have all been tested, so there's not much
# else to test than to just check it "works"
v = Vector{PDSampler.AllowedVarType}(undef, chain.structure.nvars)
for i in 1:length(v)
    v[i] = PDSampler.ls_random(all_evlist.evl[i])
end
pq2 = PDSampler.PriorityQueue{Int, Real}()
for fidx in 1:chain.structure.nfactors
    (xf, vf, g, vars) = PDSampler.ls_retrieve(chain, fidx, all_evlist, t)
    PDSampler.ls_updatepq!(pq2, chain, fidx, xf, v[vars], g, t)
end
@test pq == pq2

#### TEST LS_FIRSTBRANCH!
fidx = 2
all_evlist_copy  = deepcopy(all_evlist)
pq_copy          = deepcopy(pq)
Random.seed!(9)
(all_evlist, pq) = PDSampler.ls_firstbranch!(chain, fidx, all_evlist, pq, t)
## reproduce part of the action
Random.seed!(9)
(xf, vf, g, vars) = PDSampler.ls_retrieve(chain, fidx, all_evlist_copy, t, true)
PDSampler.ls_saveupdate!(all_evlist_copy, vars, xf, vf, t)
PDSampler.ls_updatepq!(pq_copy, chain, fidx, xf, vf, g, t)
# update a single attached factor
fpidx = linkedfactors(chain, fidx)[1]
(xfp, vfp, gp) = PDSampler.ls_retrieve(chain, fpidx, all_evlist_copy, t)
PDSampler.ls_updatepq!(pq_copy, chain, fpidx, xfp, vfp, gp, t)

@test   getlastevent(all_evlist.evl[1]).x ==
        getlastevent(all_evlist_copy.evl[1]).x
@test pq_copy[fpidx] == pq[fpidx]


text:
@sumProductRule(:node_type     => Probit,
                :outbound_type => Message{Bernoulli},
                :inbound_types => (Nothing, Message{Gaussian}),
                :name          => SPProbitOutNG)

@sumProductRule(:node_type     => Probit,
                :outbound_type => Message{Function},
                :inbound_types => (Message{PointMass}, Nothing),
                :name          => SPProbitIn1PN)

@sumProductRule(:node_type     => Probit,
                :outbound_type => Message{Function},
                :inbound_types => (Message{Bernoulli}, Nothing),
                :name          => SPProbitIn1BN)

@expectationPropagationRule(:node_type     => Probit,
                            :outbound_type => Message{GaussianWeightedMeanPrecision},
                            :inbound_types => (Message{Bernoulli}, Message{Gaussian}),
                            :outbound_id   => 2,
                            :name          => EPProbitIn1BG)

@expectationPropagationRule(:node_type     => Probit,
                            :outbound_type => Message{GaussianWeightedMeanPrecision},
                            :inbound_types => (Message{Categorical}, Message{Gaussian}),
                            :outbound_id   => 2,
                            :name          => EPProbitIn1CG)

@expectationPropagationRule(:node_type     => Probit,
                            :outbound_type => Message{GaussianWeightedMeanPrecision},
                            :inbound_types => (Message{PointMass}, Message{Gaussian}),
                            :outbound_id   => 2,
                            :name          => EPProbitIn1PG)


text:
using Gadfly, Colors

plot(ones(10,3).*[1 2 3], x=Row.index, y=Col.value, color=Col.index,
        Geom.line, Scale.color_discrete_hue(x->linspace(RGB(1,0,0),RGB(0,0,1),x)))


text:
## for development

using DSP
using HDF5

if !isdefined(:WarpedArray)
    include("types.jl")
end

include("rasta.jl")
include("mfccs.jl")
include("warpedarray.jl")
include("feacalc.jl")
include("io.jl")


text:
# Angular quadrature defined on the unit sphere octant in the upper right, closest to the viewer.
# The angles and weights are transformed to the other octants by symmetry.
#     +----+----+
#    /    /    /|
#   +----+----+ |
#  /    /    /| +
# +----+----+ |/|
# |    |this| + |
# |    | one|/| +
# +----+----+ |/
# |    |    | +
# |    |    |/
# +----+----+
#
# The spherical coordinates are defined in the following manner
# Ω̂ = (Ω_i, Ω_j, Ω_k) = (cos(θ),  sin(θ)cos(γ),  sin(θ)sin(γ))
#                     = (     μ, √(1-μ²)cos(γ), √(1-μ²)sin(γ))
#
#        j
#        ^
#        |   θ is the polar angle about the i-axis (x-direction)
#        |   γ is the azimuthal angle in the j-k plane, from the j-axis
#        |
#        |
#        |
#       /|
#      (γ|
#       \|--------------------> i
#       / \θ)
#      /   \
#     /     \
#    /       \ Ω̂
#   /         v
#  𝘷
#  k

# TODO: Function to convert product into general quadrature
#struct GeneralAngularQuadrature{T <: AbstractFloat} <: AngularQuadrature
#    Ω̂::Tuple{Point{T}} # Points on the unit sphere satisfying θ ∈ (0, π/2), γ ∈ (0, π)
#    w::Tuple{T} # Weights for each point
#end

struct ProductAngularQuadrature{nγ, nθ} <: AngularQuadrature
    γ::SVector{nγ, Float64}    # Azimuthal angles, γ ∈ (0, π)
    w_γ::SVector{nγ, Float64}  # Weights for the azimuthal angles
    θ::SVector{nθ, Float64}    # Polar angles, θ ∈ (0, π/2)
    w_θ::SVector{nθ, Float64}  # Weights for the polar angles
end

function generate_chebyshev_angular_quadrature(M::Int64)
    # A Chebyshev-type quadrature for a given weight function is a quadrature formula with equal
    # weights. This function produces evenly spaced angles with equal weights.
    angles = [(π*(2m-1)/(4M)) for m = 1:M]
    weights = zeros(M) .+ 1/M
    return angles, weights
end

# nγ and nθ are azimuthal and polar angles per octant
function generate_angular_quadrature(quadrature_type::String, nγ::Int, nθ::Int)
    if quadrature_type == "Chebyshev-Chebyshev"
        (azi_angles, azi_weights) = generate_chebyshev_angular_quadrature(nγ)
        (pol_angles, pol_weights) = generate_chebyshev_angular_quadrature(nθ)
        append!(azi_angles, reverse(π .- azi_angles))
        azi_weights = azi_weights./2
        append!(azi_weights, azi_weights)
        quadrature = ProductAngularQuadrature(SVector{2nγ}(azi_angles), SVector{2nγ}(azi_weights),
                                              SVector{nθ}(pol_angles), SVector{nθ}(pol_weights))
    else
        @error "Unsupported quadrature type"
    end
    return quadrature
end

max_stars_repo_name: KyleVaughn/MOCNeutronTransport

max_issues_repo_name: KyleVaughn/MOCNeutronTransport

max_forks_repo_name: KyleVaughn/MOCNeutronTransport