text:
# # Heat Equation
#
# ![](heat_square.png)
#
# *Figure 1*: Temperature field on the unit square with an internal uniform heat source
# solved with homogeneous Dirichlet boundary conditions on the boundary.
#
#-
#md # !!! tip
#md #     This example is also available as a Jupyter notebook:
#md #     [`heat_equation.ipynb`](@__NBVIEWER_ROOT_URL__/examples/heat_equation.ipynb).
#-
#
# ## Introduction
#
# The heat equation is the "Hello, world!" equation of finite elements.
# Here we solve the equation on a unit square, with a uniform internal source.
# The strong form of the (linear) heat equation is given by
#
# ```math
#  -\nabla \cdot (k \nabla u) = f  \quad x \in \Omega,
# ```
#
# where $u$ is the unknown temperature field, $k$ the heat conductivity,
# $f$ the heat source and $\Omega$ the domain. For simplicity we set $f = 1$
# and $k = 1$. We will consider homogeneous Dirichlet boundary conditions such that
# ```math
# u(x) = 0 \quad x \in \partial \Omega,
# ```
# where $\partial \Omega$ denotes the boundary of $\Omega$.
#
# The resulting weak form is given by
# ```math
# \int_{\Omega} \nabla v \cdot \nabla u \ d\Omega = \int_{\Omega} v \ d\Omega,
# ```
# where $v$ is a suitable test function.
#-
# ## Commented Program
#
# Now we solve the problem in Ferrite. What follows is a program spliced with comments.
#md # The full program, without comments, can be found in the next [section](@ref heat_equation-plain-program).
#
# First we load Ferrite, and some other packages we need
using Ferrite, SparseArrays
# We start  generating a simple grid with 20x20 quadrilateral elements
# using `generate_grid`. The generator defaults to the unit square,
# so we don't need to specify the corners of the domain.
grid = generate_grid(Quadrilateral, (20, 20));

# ### Trial and test functions
# A `CellValues` facilitates the process of evaluating values and gradients of
# test and trial functions (among other things). Since the problem
# is a scalar problem we will use a `CellScalarValues` object. To define
# this we need to specify an interpolation space for the shape functions.
# We use Lagrange functions (both for interpolating the function and the geometry)
# based on the reference "cube". We also define a quadrature rule based on the
# same reference cube. We combine the interpolation and the quadrature rule
# to a `CellScalarValues` object.
dim = 2
ip = Lagrange{dim, RefCube, 1}()
qr = QuadratureRule{dim, RefCube}(2)
cellvalues = CellScalarValues(qr, ip);

# ### Degrees of freedom
# Next we need to define a `DofHandler`, which will take care of numbering
# and distribution of degrees of freedom for our approximated fields.
# We create the `DofHandler` and then add a single field called `u`.
# Lastly we `close!` the `DofHandler`, it is now that the dofs are distributed
# for all the elements.
dh = DofHandler(grid)
push!(dh, :u, 1)
close!(dh);

# Now that we have distributed all our dofs we can create our tangent matrix,
# using `create_sparsity_pattern`. This function returns a sparse matrix
# with the correct elements stored.
K = create_sparsity_pattern(dh)

# ### Boundary conditions
# In Ferrite constraints like Dirichlet boundary conditions
# are handled by a `ConstraintHandler`.
ch = ConstraintHandler(dh);

# Next we need to add constraints to `ch`. For this problem we define
# homogeneous Dirichlet boundary conditions on the whole boundary, i.e.
# the `union` of all the face sets on the boundary.
∂Ω = union(getfaceset.((grid, ), ["left", "right", "top", "bottom"])...);

# Now we are set up to define our constraint. We specify which field
# the condition is for, and our combined face set `∂Ω`. The last
# argument is a function which takes the spatial coordinate $x$ and
# the current time $t$ and returns the prescribed value. In this case
# it is trivial -- no matter what $x$ and $t$ we return $0$. When we have
# specified our constraint we `add!` it to `ch`.
dbc = Dirichlet(:u, ∂Ω, (x, t) -> 0)
add!(ch, dbc);

# We also need to `close!` and `update!` our boundary conditions. When we call `close!`
# the dofs which will be constrained by the boundary conditions are calculated and stored
# in our `ch` object. Since the boundary conditions are, in this case,
# independent of time we can `update!` them directly with e.g. $t = 0$.
close!(ch)
update!(ch, 0.0);

# ### Assembling the linear system
# Now we have all the pieces needed to assemble the linear system, $K u = f$.
# We define a function, `doassemble` to do the assembly, which takes our `cellvalues`,
# the sparse matrix and our DofHandler as input arguments. The function returns the
# assembled stiffness matrix, and the force vector.
function doassemble(cellvalues::CellScalarValues{dim}, K::SparseMatrixCSC, dh::DofHandler) where {dim}
    # We allocate the element stiffness matrix and element force vector
    # just once before looping over all the cells instead of allocating
    # them every time in the loop.
    #+
    n_basefuncs = getnbasefunctions(cellvalues)
    Ke = zeros(n_basefuncs, n_basefuncs)
    fe = zeros(n_basefuncs)

    # Next we define the global force vector `f` and use that and
    # the stiffness matrix `K` and create an assembler. The assembler
    # is just a thin wrapper around `f` and `K` and some extra storage
    # to make the assembling faster.
    #+
    f = zeros(ndofs(dh))
    assembler = start_assemble(K, f)

    # It is now time to loop over all the cells in our grid. We do this by iterating
    # over a `CellIterator`. The iterator caches some useful things for us, for example
    # the nodal coordinates for the cell, and the local degrees of freedom.
    #+
    for cell in CellIterator(dh)
        # Always remember to reset the element stiffness matrix and
        # force vector since we reuse them for all elements.
        #+
        fill!(Ke, 0)
        fill!(fe, 0)

        # For each cell we also need to reinitialize the cached values in `cellvalues`.
        #+
        reinit!(cellvalues, cell)

        # It is now time to loop over all the quadrature points in the cell and
        # assemble the contribution to `Ke` and `fe`. The integration weight
        # can be queried from `cellvalues` by `getdetJdV`.
        #+
        for q_point in 1:getnquadpoints(cellvalues)
            dΩ = getdetJdV(cellvalues, q_point)
            # For each quadrature point we loop over all the (local) shape functions.
            # We need the value and gradient of the testfunction `v` and also the gradient
            # of the trial function `u`. We get all of these from `cellvalues`.
            #+
            for i in 1:n_basefuncs
                v  = shape_value(cellvalues, q_point, i)
                ∇v = shape_gradient(cellvalues, q_point, i)
                fe[i] += v * dΩ
                for j in 1:n_basefuncs
                    ∇u = shape_gradient(cellvalues, q_point, j)
                    Ke[i, j] += (∇v ⋅ ∇u) * dΩ
                end
            end
        end

        # The last step in the element loop is to assemble `Ke` and `fe`
        # into the global `K` and `f` with `assemble!`.
        #+
        assemble!(assembler, celldofs(cell), fe, Ke)
    end
    return K, f
end
#md nothing # hide

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

# To account for the boundary conditions we use the `apply!` function.
# This modifies elements in `K` and `f` respectively, such that
# we can get the correct solution vector `u` by using `\`.
apply!(K, f, ch)
u = K \ f;

# ### Exporting to VTK
# To visualize the result we export the grid and our field `u`
# to a VTK-file, which can be viewed in e.g. [ParaView](https://www.paraview.org/).
vtk_grid("heat_equation", dh) do vtk
    vtk_point_data(vtk, dh, u)
end

## test the result                #src
using Test                        #src
@test norm(u) ≈ 3.307743912641305 #src

#md # ## [Plain program](@id heat_equation-plain-program)
#md #
#md # Here follows a version of the program without any comments.
#md # The file is also available here: [`heat_equation.jl`](heat_equation.jl).
#md #
#md # ```julia
#md # @__CODE__
#md # ```


text:
using Plots
using Colors, StaticArrays
using StaticArrays.FixedSizeArrays

const FuncOrFuncs{F} = F

# ==========================================================
"""
    plot_funcs

shorthands for plot functions from Plots.jl
"""
plot_funcs = [
    bar,
    barh,
    barhist,
    boxplot,
    contour,
    contour3d,
    contourf,
    curves,
    density,
    heatmap,
    hexbin,
    histogram,
    histogram2d,
    hline,
    hspan,
    ohlc,
    path3d,
    plots_heatmap,
    quiver,
    scatter,
    scatter3d,
    scatterhist,
    stephist,
    sticks,
    surface,
    violin,
    vline,
    vspan,
    wireframe
]

# ==========================================================
"""
    signatures

function signatures from Plots.jl
"""
signatures = [
    # 1 argument

    # An AbstractMatrix where the type is the union of Integer and AbstractFloat.
    Tuple{AbstractMatrix{T} where T<:Union{Integer,AbstractFloat},},
    # --> supported by function convert_arguments(P, data::AbstractMatrix)

    # Tuple{Formatted{T} where T<:AbstractMatrix,},

    # An AbstractArray where the type is Number, and the dimension is 3
    Tuple{AbstractArray{T,3} where T<:Number,},
    # --> supported by function convert_arguments(MT::Type{Mesh, xyz::AbstractVector{<: VecTypes{3, T}})

    # An AbstractMatrix where the type is of a ColorTypes.Gray
    Tuple{AbstractMatrix{T} where T<:Gray,},
    # --> handled by convert_attribute

    # An AbstractMatrix where the type is of a ColorTypes.Colorant
    Tuple{AbstractMatrix{T} where T<:Colorant,},
    # --> handled by convert_attribute

    # plotting arbitrary shapes/polygons
    #Tuple{Shape,},
    #Tuple{AbstractVector{Shape},},
    #Tuple{AbstractMatrix{Shape},},

    # function without range... use the current range of the x-axis

    # A function
    Tuple{F} where F<:Function,

    # 2 arguments

    # A function and a number
    Tuple{F, Number} where F<:Function,

    # 3 arguments
    # 3d line or scatter

    # Three AbstractVector's
    Tuple{AbstractVector,AbstractVector,AbstractVector,},
    # --> supported by function convert_arguments(P, x::AbstractVector, y::AbstractVector, z::AbstractVector, i::AbstractArray{T, 3}) where T

    # Three AbstractMatrix's
    Tuple{AbstractMatrix,AbstractMatrix,AbstractMatrix,},
    # --> supported by function convert_arguments(P, x::AbstractMatrix, y::AbstractMatrix, z::AbstractMatrix)

    # Two AbstractVector's and a function
    Tuple{AbstractVector,AbstractVector,Function,},

    # Two AbstractVector's and an AbstractMatrix
    Tuple{AbstractVector,AbstractVector,AbstractMatrix,},
    # --> supported by function convert_arguments(P, x::AbstractVector, y::AbstractVector, z::AbstractMatrix)

    # Two AbstractVector's and an AbstractMatrix with the type ColorTypes.Gray
    Tuple{AbstractVector,AbstractVector,AbstractMatrix{T} where T<:Gray},
    # --> supported like above with AVeC, AVec, AMat

    # Two AbstractVector's and an AbstractMatrix with the type ColorTypes.Colorant
    Tuple{AbstractVector,AbstractVector,AbstractMatrix{T} where T<:Colorant},
    # --> supported like above with AVeC, AVec, AMat

    #parametric functions

    # A function and two numbers
    Tuple{Function, Number, Number,},

    # An AbstractArray of functions, and two numbers
    Tuple{AbstractArray{F}, Number, Number,} where F<:Function,

    # Tuple{FuncOrFuncs{F}, FuncOrFuncs{G}, AbstractVector}  where {F<:Function,G<:Function},
    # Tuple{FuncOrFuncs{F}, FuncOrFuncs{G}, Number, Number, Number} where {F<:Function,G<:Function},
    # Tuple{FuncOrFuncs{F}, FuncOrFuncs{G}, FuncOrFuncs{H}, AbstractVector} where {F<:Function,G<:Function,H<:Function},
    # Tuple{FuncOrFuncs{F}, FuncOrFuncs{G}, FuncOrFuncs{H}, Number, Number, Number} where {F<:Function,G<:Function,H<:Function},

    # Lists of tuples and Tuple{FixedSizeArrays,},
    # Tuple{Tuple},

    # (x,y) tuples

    # An AbstractVector of the Tuple of R1 and R2, where both R1 and R2 are numbers
    Tuple{AbstractVector{Tuple{R1,R2}}} where {R1<:Number,R2<:Number},

    # A Tuple of R1 and R2, where both R1 and R2 are numbers
    Tuple{Tuple{R1,R2}} where {R1<:Number,R2<:Number},

    # (x,y,z) tuples

    # An AbstractVector of the Tuple of R1, R2 and R3, where R1, R2 and R3 are numbers
    Tuple{AbstractVector{Tuple{R1,R2,R3}}} where {R1<:Number,R2<:Number,R3<:Number},

    # A Tuple of R1, R2 and R3, where R1, R2 and R3 are numbers
    Tuple{Tuple{R1,R2,R3}} where {R1<:Number,R2<:Number,R3<:Number},

    # these might be points+velocity, or OHLC or something else

    # An AbstractVector of the Tuple of R1, R2, R3 and R4, where R1, R2, R3 and R4 are numbers
    Tuple{AbstractVector{Tuple{R1,R2,R3,R4}}} where {R1<:Number,R2<:Number,R3<:Number,R4<:Number},

    # A Tuple of R1, R2, R3 and R4, where R1, R2, R3 and R4 are numbers
    Tuple{Tuple{R1,R2,R3,R4}} where {R1<:Number,R2<:Number,R3<:Number,R4<:Number},

    # 2D Tuple{FixedSizeArrays,},

    # An AbstractVector of a Vector of a fixed size, with the type Number and dimension 2
    Tuple{AbstractVector{FixedSizeArrays.Vec{2,T}}} where {T<:Number},

    # A Vector of a fixed size, with the type Number and dimension 2
    Tuple{FixedSizeArrays.Vec{2,T}} where {T<:Number},

    # 3D Tuple{FixedSizeArrays,},

    # An AbstractVector of a Vector of a fixed size, with the type Number and dimension 3
    Tuple{AbstractVector{FixedSizeArrays.Vec{3,T}}} where {T<:Number},

    # A Vector of a fixed size, with the type Number and dimension 3
    Tuple{FixedSizeArrays.Vec{3,T}} where {T<:Number}
]

# ==========================================================
# basic generator function

basic_generator(x) = sin(x)
basic_generator(x, y) = sin(x) + cos(y)
basic_generator(x, y, z) = sin(x) + cos(y) + (cos(z) * sin(z))

# convert Types and assign values
function to_values(::Type{Union{T1, T2, T3}}) where {T1, T2, T3}
    T1
end
to_values(::Type{T}) where T <: Function = basic_generator
function to_values(::Type{<: AbstractArray{T, N} where T}) where N
    rand(ntuple(i->10, N)...)
end
function to_values(::Type{<: AbstractArray{<: Function, N}}) where N
    fill(basic_generator, ntuple(i->10, N)...)
end
function to_values(::Type{<: AbstractArray{<: Function}})
    fill(basic_generator, 10)
end
# function to_values(::Type{A}) where A <: AbstractArray{<: Function, N} where N
#     rand(basic_generator, ntuple(i->10, N)...)
# end
to_values(::Type{Number}) = rand()
to_values(::Type{T}) where T <: Number = rand(T)
to_values(::Type{Tuple{}}) = ()
function to_values(::Type{T}) where T <: Tuple
    T1 = Base.tuple_type_head(T)
    Ttail = Base.tuple_type_tail(T)
    (to_values(T1), to_values(Ttail)...)
end

# for signature in signatures
#     # println(to_values(signature))
#     Plots.scatter(to_values(signature))
# end
#
# Plots.scatter(to_values(signatures[1]))
# generated_values = to_values.(signatures)

# ==========================================================
# batch testing

plots = []
for func in plot_funcs
    for signature in signatures
        generated_values = to_values(signature)
        # println("Testing $func with $signature")
        try
            plt = func(generated_values...)
            push!(plots, (func, generated_values, plt))
            println("No errors thrown for $func with $signature")
        catch
            println("FAILED: $func with $signature")
        end
    end
end

plots

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

# displaying outputs
plt = plot(plots[60:64][3]...) #the "..." means squashing

# to show plots
plot(plots[N][3])

# to show the generated values
plots[N][2]

# to show the function
plots[N][1]

ext: jl


text:
using PathLoss
using Base.Test

# write your own tests here
@test 1 == 1

# Distance test
lat1 = 38.893
lon1 = -77.037852
lat2 = 38.899147
lon2 = -77.023934
d_calc = distanceInKm(lat1,lon1,lat2,lon2)
d_real = 1.385 # km
@test_approx_eq_eps d_calc d_real 0.01 # precision 10 m

# Azimuth test
az_real = 60.42
az_calc  = azimuthLineAtoB(lat1,lon1,lat2,lon2)
@test_approx_eq_eps az_calc az_real 0.1 # precision o.1

# FreeSpaceModel Test
fs = FreeSpaceModel() # For ERP
fs.freq =800 #Mhz
p_fs1 = pathloss(fs,1.385)
p_fs2 = 93.33
@test_approx_eq_eps p_fs1 p_fs2 0.01 # precision 0.01 Db


text:
module ReadCFGPARA
    include("readCsvLib.jl")
    include("mappingLib.jl")

    CsvInfoObj = CsvInfos()

    CsvInfoObj.防区表fieldname ="定位界标1,定位界标2,xxxx主机IP,通道号,定位界标1光程,定位界标2光程,子防区编号,防区起点GPS经度坐标,防区起点GPS纬度坐标,防区终点GPS经度坐标,防区终点GPS纬度坐标,左右岸,地名信息,界标之间围栏长度,管理处,分局"
    CsvInfoObj.防区表signalStringName="定位界标1光程"
    CsvInfoObj.防区表keyFieldNames ="定位界标1,定位界标1光程,定位界标2,定位界标2光程"
    CsvInfoObj.防区表filename = "../Jl-Aux/glbCfg//防区信息.csv"

    CsvInfoObj.主机信息表fieldname = "主机编号,友好名,主机位置,主机类型"
    CsvInfoObj.主机信息表signalStringName = "主机类型"
    CsvInfoObj.主机信息表keyFieldNames = "主机编号,友好名,主机位置,主机类型"
    CsvInfoObj.主机信息表filename     = "../Jl-Aux/glbCfg/A100-xxxx主机IP信息表.csv"

    CsvInfoObj.老版参数表fieldname ="编号,通道,标牌坐标1,标牌坐标2,动态阈值时长,动态阈值系数,动态阈值下限,动态阈值上限,事件合并时长,事件滤除时间,事件合并距离,事件滤除距离,告警空间下限,告警空间上限,告警续时下限,告警续时上限,地名标识,说明"
    CsvInfoObj.老版参数表signalStringName="标牌坐标1"
    CsvInfoObj.老版参数表keyFieldNames ="通道,标牌坐标1,标牌坐标2,动态阈值时长,动态阈值系数,动态阈值下限,动态阈值上限,事件合并时长,事件滤除时间,事件合并距离,事件滤除距离,告警空间下限,告警空间上限,告警续时下限,告警续时上限"
    CsvInfoObj.老版参数表filename ="../Jl-Aux/glbCfg/xxxx-AlarmParas.csv"
    CsvInfoObj.ParasPosNames  = "通道,标牌坐标1,标牌坐标2"

    function Readin_cfgParas(;machine友好名0="GYYS001" )
        CsvInfoObj.老版参数表filename = replace(CsvInfoObj.老版参数表filename, Pair("xxxx",machine友好名0) )
        # println( " Here -----------------------00--------------------")
        basicInit( CsvInfoObj)

        convertPara_Fiberlength( CsvInfoObj,machine友好名1=machine友好名0)
        # println( " Here ------------------------11-------------------")
        return
    end

    function Paras_tofloatArray(;metersPerSeg::Float64,timeUnit::Float64=0.426)::Array
        fieldNum =  length( split(CsvInfoObj.老版参数表keyFieldNames,",") )
        rowdata = CsvInfoObj.csvInfoDefencePara.Area2FiberLen
        rowNum = length( rowdata )
     #   println( ":::rowdata len is $rowNum")
        rows = Array{Float64,2}(undef,rowNum,fieldNum)

        spacefactor1s=[2,3, 10,11,12,13 ]
        # timefactors =[4,8,9,14,15]
        timefactors =[4,14,15] # IN 0.426秒

        for row =1: rowNum
            fieldidx = 1
            rows[row,fieldidx]= parse(Float64, rowdata[row][ fieldidx ] )
            fieldidx = 2
            rows[row,fieldidx] = rowdata[row][ fieldidx ]
            # rows[row,fieldidx]
            fieldidx =3
            rows[row,fieldidx] =  rowdata[row][ fieldidx ]
            for fieldidx =4:fieldNum
                rows[row,fieldidx]= parse(Float64,rowdata[row][ fieldidx ] )

            end

            for fieldidx in spacefactor1s
                rows[row,fieldidx] = rows[row,fieldidx]/metersPerSeg +0.01
                # println( rows[row,fieldidx], "->",rowdata[row][ fieldidx ]," in unit :" ,metersPerSeg)
            end
            # println()
            # for fieldidx in timefactors
            #     rows[row,fieldidx] = rows[row,fieldidx]/ timeUnit + 0.01
            # end
            # rows[row,5]  =rows[row,5]
           # "通道,标牌坐标1,标牌坐标2,动态阈值时长,动态阈值系数,动态阈值下限,动态阈值上限,事件合并时长,事件滤除时间,事件合并距离,事件滤除距离,告警空间下限,告警空间上限,告警续时下限,告警续时上限"

        end
        # return rows
        function transfer(data::Float64)
            return floor( UInt, data)
        end

        rows1 = Array{UInt,2}(undef,rowNum,fieldNum)
        rows1 = map(transfer, rows)
        return rows1
    end

    function initParaDict( ;paraDict::Dict,rowsObj::Array,chn1SegIDRng::Array , chn2SegIDRng::Array )
        # println(chn2SegIDRng ," **************** ",chn1SegIDRng  )
        # exit()

        # println("----------------------------------")
        # rowNum,fieldNum   = size( rowsObj )
        # println("--------------#$rowNum,$fieldNum--------------------")
        # for row =1: rowNum
        #         # l=fieldNum
        #         for  idx = 1:fieldNum
        #             print( rowsObj[ row,idx] ," --")
        #         end
        #         println()
        # end
        # exit()

        rowNum,fieldNum   = size( rowsObj )
        # for row =1:rowNum
        for rowidx =1: rowNum # in rowsObj
            rowobj = rowsObj[rowidx,:]
            # println( size( rowobj))
            #
            # l1 =length(rowobj )
            # for i =1:l1
            #         println("$i,$(rowobj[i]) ")
            # end
            # println("=====================================================")
            if rowobj[ 2 ] == UInt(0) &&  rowobj[ 3 ]  == UInt(0) ## the default para

                if rowobj[ 1 ] == 1 #Channel 1
                    # println("default para Fill................1............",chn1SegIDRng[1],";",chn1SegIDRng[2])
                    for idx = chn1SegIDRng[1] : chn1SegIDRng[2]
                         paraDict[ idx ] = rowobj[4:end]
                    end
                elseif rowobj[ 1 ] == 2 #Channel 2

                    for idx = chn2SegIDRng[1] : chn2SegIDRng[2]

                        paraDict[ idx ] = rowobj[4:end]
                    end
                end

            end
        end

        for rowidx = 1: rowNum # in rowsObj
            rowobj = rowsObj[rowidx,:]
            startSegID = rowobj[ 2 ]
            endSegID   = rowobj[ 3 ]
            # println( "...............$( rowobj[ 1 ])_$startSegID,  $endSegID")
            if startSegID != 0 ||  endSegID != 0 ## NOT the default para
                if rowobj[ 1 ] == 1
                    offsetSegID = 0
                elseif rowobj[ 1 ] == 2
                    offsetSegID = chn2SegIDRng[1]-1
                else
                   ;# println( "$startSegID,  $endSegID")
                end
                begid = offsetSegID + startSegID
                endid= offsetSegID + endSegID
                rngid= begid: endid
                #println( "----$begid: $endid----" )
                for idx = rngid
                    paraDict[ idx ] = rowobj[4:end]
                end
            end
        end
    # println("======================================")
    # v=[]
    # for i in keys(paraDict )
    #     append!(v,i)
    #     # println( i )
    # end
    # a= sort(v )
    # for a1 =a
    #     println(a1)
    # end

    # exit()
    end
end #module ReadCFGPARA
function MainCfgPara(;machine友好名0="GYYS-001",metersPerSeg1::Float64,timeUnit1::Float64=0.426,Chn1Seg1::Int =1,Chn1Seg2::Int =4196, Chn2Seg1::Int =4197, Chn2Seg2::Int =8192 )
    ReadCFGPARA.Readin_cfgParas(machine友好名0=machine友好名0 )
    # println( " Here -------------------------------------------")
    rows = ReadCFGPARA.Paras_tofloatArray( metersPerSeg = metersPerSeg1,timeUnit=timeUnit1 )

    paraDictObj = Dict()
    chn1SegIDRng=[Chn1Seg1,Chn1Seg2]
    chn2SegIDRng=[Chn2Seg1,Chn2Seg2]
    ReadCFGPARA.initParaDict(paraDict=paraDictObj ,rowsObj=rows  ,chn1SegIDRng=chn1SegIDRng , chn2SegIDRng=chn2SegIDRng)
    # v=[]
    # for i in keys(paraDictObj )
    #     append!(v,i)
    #     # println( i )
    # end
    # a= sort(v )
    # for a1 =a
    #     println(a1)
    # end
    return paraDictObj
end

function TestMain(;machine友好名0="GYYS001" )
    ReadCFGPARA.Readin_cfgParas(machine友好名0=machine友好名0 )

    rows = ReadCFGPARA.Paras_tofloatArray( metersPerSeg = Float64(6.001) )
    # show(ReadCFGPARA.CsvInfoObj.csvInfoDefencePara )
    rowNum,fieldNum   = size( rows )

    for row =1:rowNum
        for fieldidx = 1:fieldNum
            print( rows[row,fieldidx],"," )
        end
        println()
    end
    paraDictObj = Dict()
    chn1SegIDRng=[1,8000]
    chn2SegIDRng=[9000,9000+8000]
    ReadCFGPARA.initParaDict(paraDict=paraDictObj ,rowsObj=rows  ,chn1SegIDRng=chn1SegIDRng , chn2SegIDRng=chn2SegIDRng)
    for key = chn1SegIDRng[1] :chn1SegIDRng[2] # 9000:17000)
        try
            valobj = paraDictObj[key]
            print(key,": ")
            for idx =1: length(valobj )
                print( valobj[idx] ,",  ")
            end
            println()
        catch( e )
            # println(e)
        end
    end

    for key in   chn2SegIDRng[1] :chn2SegIDRng[2]# 9000:17000)
        try
            valobj = paraDictObj[key]
            print(key,": ")
            for idx =1: length(valobj )
                print( valobj[idx] ,",  ")
            end
            println()
        catch( e )
            # println(e)
        end
    end

end
# TestMain(machine友好名0="GYYS001" )

max_stars_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_issues_repo_name: wangzhongtian/fiberAlarmAnalyseTool

max_forks_repo_name: wangzhongtian/fiberAlarmAnalyseTool


text:
using Test

include("")


text:
## AbstractHit ##

abstract type AbstractHit end

function Base.show(io::IO, h::H) where {H<:AbstractHit}
    x, y, z = location(h)
    print(io, H, "(", x, ", ", y, ", ", z, ", ")
    print(io, energy(h), ", ", time(h), ", ")
    print(io, particleid(h), ", ", trackid(h), ", ", trackparentid(h), ")")
end

## Hit ##

struct Hit <: AbstractHit
    x::Float64
    y::Float64
    z::Float64
    E::Float64
    t::Float64
    particleid::Int
    trackid::Int
    trackparentid::Int
end

location(h::Hit) = (h.x, h.y, h.z)
energy(h::Hit) = h.E
Base.time(h::Hit) = h.t
particleid(h::Hit) = h.particleid
trackid(h::Hit) = h.trackid
trackparentid(h::Hit) = h.trackparentid

## AbstractHitIter ##

abstract type AbstractHitIter end

Base.IteratorSize(::Type{<:AbstractHitIter}) = Base.HasLength()
Base.IteratorEltype(::Type{<:AbstractHitIter}) = Base.HasEltype()

## Abstract Event ##

abstract type AbstractEvent end

Base.IteratorSize(::Type{<:AbstractEvent}) = Base.HasLength()
Base.length(e::AbstractEvent) = hitcount(e)

Base.IteratorEltype(::Type{<:AbstractEvent}) = Base.HasEltype()

Base.iterate(e::AbstractEvent) = iterate(hits(e))
Base.iterate(e::AbstractEvent, state) = iterate(hits(e), state)

function Base.show(io::IO, e::AbstractEvent)
    print(io, "Event(")
    print(io, eventnum(e), ", ", hitcount(e), ", ", primarycount(e), ", ")
    print(IOContext(io, :limit => true, :compact => true), hits(e))
    print(")")
end

function Base.show(io::IO, mime::MIME"text/plain", e::AbstractEvent)
    print(io, typeof(e), " with ")
    print(io, "eventnum: ", eventnum(e), ", hitcount: ", hitcount(e))
    print(io, ", primarycount: ", primarycount(e), " and hits:\n")
    show(IOContext(io, :limit => true), mime, hits(e))
end

energy(event::E) where {E<:AbstractEvent} = sum(energy, hits(event))

## Event ##

struct Event{T} <: AbstractEvent
    eventnum::Int
    hitcount::Int
    primarycount::Int
    hits::T
    function Event{T}(eventnum, hitcount, primarycount, hits) where {T}
        if length(hits) == hitcount
            return new(eventnum, hitcount, primarycount, hits)
        end
        throw(ArgumentError("hitcount $hitcount must equal number of hits $(length(hits))"))
    end
end

function Event(eventnum, hitcount, primarycount, hits::T) where {T}
    Event{T}(eventnum, hitcount, primarycount, hits)
end

function Event{Vector{H}}(eventnum, hitcount, primarycount, hits::AbstractHitIter) where {H}
    Event{Vector{H}}(eventnum, hitcount, primarycount, collect(H, hits))
end

function Base.convert(::Type{Event{Vector{H}}}, e::Event{<:AbstractHitIter}) where {H}
    return Event{Vector{H}}(eventnum(e), hitcount(e), primarycount(e), hits(e))
end

Base.eltype(::Type{Event{T}}) where {T} = eltype(T)

hits(e::Event) = e.hits
hitcount(e::Event) = e.hitcount
eventnum(e::Event) = e.eventnum
primarycount(e::Event) = e.primarycount


text:
# This code is part of QuantumCircuits.
#
# (C) Copyright Rafał Pracht 2022.
#
# This code is licensed under the Apache License, Version 2.0. You may
# obtain a copy of this license in the LICENSE.txt file in the root directory
# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
#
# Any modifications or derivative works of this code must retain this
# copyright notice, and modified files need to carry a notice indicating
# that they have been altered from the originals.

using QuantumCircuits
using QuantumCircuits.QML
using QuantumCircuits.QML.Optimization
using QuantumCircuits.Execute
using QuantumCircuits.Execute: state2probability
using QuantumCircuits.QCircuits.Circuit

using Distributions: LogNormal, quantile, cdf, pdf

using PyCall

# This code was copied from classic coding competition, it is used to check solution.
py"""
from scipy.stats import lognorm
from scipy.interpolate import interp1d
import numpy as np

TARGET_ERR = 1e-2

def l2_error(pmf: np.array, x_grid: np.array, sigma=0.1):
    pmf = np.array(pmf)
    x_grid = np.array(x_grid)
    assert all(pmf >= 0)
    assert np.isclose(sum(pmf), 1)
    assert all(x_grid >= 0)
    assert all(np.diff(x_grid) > 0)
    assert len(pmf) + 1 == len(x_grid)

    n_point = 2 ** 22
    tail_value = (TARGET_ERR / 100) ** 2
    min_x = lognorm.ppf(tail_value, sigma)
    max_x = lognorm.ppf(1 - tail_value, sigma)
    x_middle = np.linspace(min_x,max_x, n_point)
    x_lower_tail = np.linspace(0, min_x, n_point//1000)
    x_upper_tail = np.linspace(max_x, x_grid[-1], n_point//1000) if x_grid[-1] > max_x else np.array([])

    x_approx = np.diff(x_grid) / 2 + x_grid[:-1]
    x_approx = np.concatenate(([x_approx[0]], x_approx, [x_approx[-1]]) )
    pdf_approx = pmf / np.diff(x_grid)
    pdf_approx = np.concatenate(([pdf_approx[0]], pdf_approx, [pdf_approx[-1]]))

    fy = interp1d(x_approx, pdf_approx, kind='nearest', assume_sorted=True, fill_value=(0, 0), bounds_error=False)
    x_full = np.concatenate((x_lower_tail[:-1], x_middle, x_upper_tail[1:]))
    approx_pdf = fy(x_full)

    full_pdf = lognorm.pdf(x_full, sigma)
    dx = np.diff(x_full)
    dx = np.append(dx, 0)

    upper_tail_err_2_approx = lognorm.sf(x_full[-1], sigma)
    main_err_2 = sum((full_pdf - approx_pdf) ** 2 * dx)
    err = (upper_tail_err_2_approx + main_err_2) ** 0.5
    return err
"""

# Target error
TARGET_ERR = 1e-2

# Definition of target distribution
μ = 0
σ = 0.1

"Cumulative normal probability distribution"
cldf(val) = cdf(LogNormal(μ, σ), val)
mydiff(v::AbstractVector) = @views v[(begin+1):end] .- v[begin:(end-1)]
pldf(val) = pdf(LogNormal(μ, σ), val)

# Qubits
N = 9

####
npoints = 2^N
tail_value = 0.00018 # 9
start_val = quantile(LogNormal(μ, σ), tail_value)
end_val = quantile(LogNormal(μ, σ), 1 - tail_value)

# Generate expected distribution manualy
diss_x = LinRange(start_val, end_val, npoints+1)
diss_p = mydiff(cldf.(diss_x))
diss_p = diss_p ./ sum(diss_p)

# Check the value
py"l2_error"(diss_p, diss_x) # 0.005128950543079309

# Generate ansact
qr = QuantumRegister(N)
cr = ClassicalRegister(N)
qc = QCircuit(qr, cr)
#add_ent_gate(i, j) = qc.cx(i, j)
add_ent_gate(i, j) = qc.u4(i, j)

#qc.x(4)
for i in 0:3
    add_ent_gate(4-i, 4-i-1)
    add_ent_gate(4+i, 4+i+1)
end

# decompose
qc = decompose(qc)
qc.measure(0:N-1, 0:N-1)

# Random parameter
params = getRandParameters(qc)
setparameters!(qc, params)

# Generate julia wrapper to call python code to check the solution
function python_l2(param, diss_x)
    start = ket"000000000"
    mat = tomatrix(qc, param)
    opt_p = qc.measures_matrix * state2probability(mat * start)

    return py"l2_error"(opt_p, diss_x)
end

# Generate expected distribution in first stage of the optimization
dx = diss_x[2] - diss_x[1]
ADDN = 2^10
const diss_y_pldf = zeros(2^N, ADDN)
for j in 1:2^N
    for i in 0:ADDN-1
        x = diss_x[j] + i/ADDN * dx
        diss_y_pldf[j, i+1] = pldf(x)
    end
end
function loss_l2_Δ(state)
    err = 0
    n = length(state)
    for j in 1:n #(x_val, p_val) in zip(diss_x[1:end-1], state)
        for i in 0:ADDN-1
            err += (diss_y_pldf[j, i+1] - state[j]/dx)^2 * dx/ADDN
        end
    end
    sqrt(err)
end

# Generate expected distribution in second stage of the optimization when the changing the discretization is alowed
function loss_l2(state, start, Δx)
    err = 0
    n = length(state)
    xx = start
    for j in 1:n
        dx = Δx[j]
        for i in 0:ADDN-1
            x = xx + i/ADDN * dx
            err += (pldf(x) - state[j]/dx)^2 * dx/ADDN
        end
        xx += dx
    end
    err += 1 - cldf(xx)
    err += cldf(start)
    sqrt(err) + max(xx - 1.75, 0) + max(0.57 - start, 0)
end

function loss_stage1(params)
    start = ket"000000000"
    mat = tomatrix(qc, params)
    state_end = state2probability(mat * start)

    return loss_l2_Δ(state_end)
end

# Loss function
function loss(params)
    Δx = abs.(params[2:513])
    model_param = params[514:end]
    start = ket"000000000"
    mat = tomatrix(qc, model_param)
    #state_end = qc.measures_matrix * state2probability(mat * start)
    state_end = state2probability(mat * start)

    #return 100*loss_mse(state_end, diss_p)
    return loss_l2(state_end, abs(params[1]), Δx)
end

# Stege 1 - only ansact
# In the beginning, we start with random parameters. During the optimization, we save in the log the best parameters
# found till now. Bat sometimes the optimization method cannot find a better point so, in that case,
# it is good to restart it from the best parameters found lastly.
# For this we use best_param, variable.
best_param = [5.679918941917893, 4.251047958788928, 3.5437338099599174, 5.441122661263436, 2.3547800750989194, 1.2375106360535937, -5.926733578722595, 5.105113492307948, 6.220803793450637, -4.731043991920058, 1.8800031306308296, 5.5737349030831105, 3.2005996610955543, 1.7527345229272375, 2.276393940822889, 6.433743641909078, 1.8804087180413855, 1.3759654438241162, 2.111986048933823, 2.646277694072601, 5.301248156653536, -0.6540576713830272, 1.9133064975706302, 6.1010852093340455, 0.6528345760801835, 3.972915237735164, 6.886297952490527, 3.8040168359747124, 3.091485664926333, 4.931157193197006, 0.553865644097088, 6.351589535169496, 5.73496498969936, 4.631313009319303, 3.515834694977789, 1.5823482363368873, 4.984427837909524, 5.867871515378414, 3.617232165066993, 3.2714246937337843, 4.64111910105604, 5.060652198181708, 3.8212182162433823, -1.4551420517097016, -0.009585091570906114, 0.995754170516813, 5.594528591259475, 3.746195066671812, 3.608551912859001, -0.8894019646716065, 1.6715257777301376, 2.2870113766942772, 5.512166184617886, 0.5714521142705539, 3.8011358974571405, 0.4677817117962828, 4.972762238325715, -0.2699783614659865, -0.25001889490376006, 4.837342114428392, 2.5394235906884934, 3.550645631398327, 0.5985420153831956, 5.611424965417167, -0.24557751621378168, 1.3413055441945323, 3.973210486309615, 4.716470372135131, -1.4402489326761916, 0.5628093958419573, -2.07108824150549, 1.0143699305757663, 5.6810691483952125, 0.31689723006715803, 8.074224906481794, 4.0008420314440345, 1.8065031699397647, 4.351954781482624, 0.27393162181003833, 1.3946594956376743, 3.068970152437304, 0.5583909059538796, -0.08708818785242389, 3.5150056229387507, 0.2860842587147386, 2.3650804786584056, 1.5834413918199446, 6.609933255262227, 5.628757478395043, 2.2812658090447404, 4.988185892342614, 2.0926004903541244, 2.7907978362966746, 1.3885592047923723, -1.7975187312537604, 0.5383417309200358, 1.9079673632105147, 5.568083164220874, 0.3042963636675245]
params = best_param

# check if loss function works
# 0.08392833431402887
loss_stage1(params)
# check if derivative of loss function works
loss_stage1'(params)

# run the optimization to find the best parameters
val, x, itr = gradientDescent(loss_stage1, loss_stage1', params, α=0.0001, maxItr=500, debug=true,
                              useBigValInc=true, argsArePeriodic=true)
python_l2(params, diss_x)

# Stage 2 - ansact and the discretization
# In the beginning, we start with parameters found in the stage 1.
# During the optimization, we save in the log the best parameters
# found till now. Bat sometimes the optimization method cannot find a better point so, in that case,
# it is good to restart it from the best parameters found lastly.
# For this we use best_param, variable.
#params = vcat([start_val], [dx for i in 1:512], best_param)
best_param = [0.6985367664022247, 0.0015250508769132493, 0.0014626137853927677, 0.0014590429186790747, 0.0014557643006992006, 0.001454049229170965, 0.0014152434090339097, 0.0014154645424014189, 0.0014274690355849293, 0.0014274598161225635, 0.0014514789740261806, 0.001452862041818621, 0.001452679344403947, 0.0014532633167190195, 0.0014534062224029842, 0.0014551139270538953, 0.0014566279666116142, 0.0014294643302957537, 0.0014292728552362352, 0.0014292486679427577, 0.001430919408368201, 0.0014541163023379018, 0.0014264683582781207, 0.0014242021400162903, 0.001422070773373309, 0.0014204928677025452, 0.0014179696641390485, 0.0014156865915134157, 0.0014124008315158862, 0.001408529509483185, 0.0014041771678095737, 0.0013973378725397319, 0.001388952298992998, 0.001570418557700881, 0.0014962322106539664, 0.001462427285624219, 0.0014459947314643237, 0.0014374257774495227, 0.001428716648156318, 0.0014245562365398049, 0.0014193533092526453, 0.0014024038001961888, 0.0013576144501029127, 0.0012866200181534816, 0.0013275133528779161, 0.0012857835770346154, 0.001227452723645235, 0.001111186406109634, 0.000983910846987461, 0.0027216896074014592, 0.0026272799861254514, 0.002329703282226615, 0.002179491224346476, 0.0020651851601318915, 0.001921146666865109, 0.001839085273350531, 0.0017395127606655075, 0.001667490450043708, 0.001589894352531833, 0.0015273881601382078, 0.0014604953385356485, 0.0014123153055424584, 0.001384503281780681, 0.0013010079121028586, 0.001263345763270071, 0.0019392437018189792, 0.0018531845905195134, 0.001765656258767773, 0.0016771608319892481, 0.0016207620536458761, 0.0015348904139311991, 0.001493174201442382, 0.0014363428176601125, 0.0013926572937170993, 0.0013445637548616521, 0.001303336750369339, 0.0012633468710517872, 0.0012145241282838202, 0.0012020632550880847, 0.0011534308481918403, 0.001131235941861219, 0.0017193376537336074, 0.0016597724888533748, 0.0015901715171418467, 0.0015368898036240673, 0.0015123767318647112, 0.0014356311396553619, 0.0014058761655063466, 0.0013590412156249114, 0.0013240567601352233, 0.0012954142560914297, 0.001261017421320202, 0.0012300604590274806, 0.0011888331122815257, 0.0011818315641344086, 0.0011398084522533415, 0.0011214758206358768, 0.0015952229340373536, 0.0015461925410152327, 0.0014995244580107022, 0.0014592495867959273, 0.0014333969000741467, 0.0013773370345070363, 0.001358835077933233, 0.001321376946287528, 0.0012995981533614211, 0.0012690075562490314, 0.001243308038767213, 0.0012173104928943585, 0.0011814401070210274, 0.0011799743246367661, 0.0011418620791040885, 0.001127538019388915, 0.0014900792572983812, 0.0014535935769283306, 0.001416737071400055, 0.0013869152062812814, 0.0013699774839846459, 0.0013234320726864385, 0.001311715064279069, 0.0012817918904826082, 0.0012662405975883608, 0.0012421022350219309, 0.001221652572814067, 0.0012010675891506494, 0.0011703450077255977, 0.0011732892508686432, 0.0011399300760484056, 0.0011295400629943355, 0.0016594652531698017, 0.0016230398622856263, 0.0015872650000226403, 0.0015574890611454932, 0.0015420416085757233, 0.0014928638191521655, 0.001483048503302747, 0.0014522684330481185, 0.0014376312813719342, 0.0014130008241249332, 0.0013931714162676737, 0.0013723542214696756, 0.001339787931259404, 0.0013457691869644605, 0.0013097756400954263, 0.001300341987439281, 0.0015507491008393008, 0.0015261637203704986, 0.0015041752896155628, 0.001484135126106822, 0.0014770029342604327, 0.0014374906455304132, 0.0014345485851729042, 0.001411519470600002, 0.0014030938952634524, 0.0013853572635658213, 0.001372876594962893, 0.001358128424334991, 0.0013313246231621044, 0.0013425723249070706, 0.0013116713644199353, 0.0013068903139337466, 0.0014971642717757075, 0.0014806648421181528, 0.0014637987052587196, 0.0014511946305855916, 0.0014507379280071727, 0.0014176945940493546, 0.0014216898840083323, 0.0014045104945161378, 0.0014028029289151278, 0.0013900960653398112, 0.0013815855798807312, 0.001371748170733806, 0.0013495241540608085, 0.0013660767411063943, 0.0013391527001339867, 0.0013388561546748396, 0.0014548163476431684, 0.0014442007458181228, 0.001431436572855298, 0.0014242463622673295, 0.001429384891909674, 0.0014021305361988597, 0.0014106402110206475, 0.0013987389556391816, 0.001401710245399672, 0.0013942765445974204, 0.0013899512841421967, 0.0013849156039217575, 0.0013671914382223681, 0.0013883239084301665, 0.001366084747485027, 0.0013705043752493574, 0.0013840065268518242, 0.0013794834232450104, 0.001374076911130199, 0.0013717731814049312, 0.0013821109336263964, 0.0013602654517263244, 0.0013734619280333395, 0.0013663365130791572, 0.0013737369401904806, 0.0013708579158846687, 0.001372215836762136, 0.00137146348633419, 0.0013568791977554935, 0.0013813002955968415, 0.0013651425292796244, 0.0013738576291753948, 0.0013966543411162767, 0.0013962886949657309, 0.0013977392050306366, 0.0014000265365980416, 0.00141440593165411, 0.001396724868848769, 0.0014139010343112094, 0.0014110635372021141, 0.001422288144436652, 0.0014238694581242843, 0.0014305376863160484, 0.0014344734044599609, 0.001425162292164838, 0.0014565362316065714, 0.0014417910054116323, 0.0014556711490415631, 0.0013787858094527255, 0.0013828199143807994, 0.0013859964322869632, 0.0013925654186953103, 0.0014113706622629332, 0.0013979391751956372, 0.001420571156880079, 0.0014221252718759232, 0.001438979791828751, 0.0014449032400958306, 0.0014549279104735925, 0.0014634873395270214, 0.0014583744222639685, 0.0014952136795143062, 0.001484861862948541, 0.0015039656551359283, 0.001354786320248275, 0.0013619791805626442, 0.0013662705364860031, 0.0013765310096760343, 0.0013988488689639967, 0.0013891921581173571, 0.0014148844074786226, 0.0014203947896928753, 0.0014415318137095939, 0.0014517868782907546, 0.0014652378648731226, 0.0014781177545579486, 0.0014775263477548076, 0.0015195358662817666, 0.0015144612819753395, 0.001538705231276322, 0.0012840137340744007, 0.0012938874230638791, 0.0013029288336586435, 0.0013156135454958054, 0.0013399901325483105, 0.0013337090722877268, 0.0013618681801404043, 0.001370192122697584, 0.0013933266340701044, 0.0014063998672897058, 0.001423895286844402, 0.0014398289402775231, 0.0014426524387608159, 0.0014869307509599626, 0.0014851529875569563, 0.0015127373214016589, 0.0012393131455630112, 0.001251561614344178, 0.0012641905073923548, 0.0012788580915241606, 0.0013054740210821098, 0.001303146783233153, 0.0013319773456067546, 0.0013437187834735053, 0.001367915289861407, 0.0013829186501197284, 0.0014047859142958473, 0.0014233761569990915, 0.0014287846546001967, 0.001476309139256843, 0.00147784427291122, 0.001508969458616347, 0.0012420577683725047, 0.0012569511149877653, 0.0012715060243388807, 0.0012884207691442088, 0.0013183224635017886, 0.001317831288422731, 0.0013513282156961545, 0.0013659313034867318, 0.0013956292353101637, 0.0014148868239542906, 0.001439609603830724, 0.0014626520104149816, 0.0014720890893140204, 0.00152547562292487, 0.0015320381466175936, 0.0015689304728814196, 0.0012389797986165198, 0.0012555562230620748, 0.0012704212005528214, 0.0012902424569442904, 0.0013223065125146135, 0.001324356830059311, 0.0013606476238089066, 0.0013781036744586422, 0.001410684671661639, 0.0014339894561811609, 0.0014617247081477665, 0.001488909420539718, 0.0015019785284637044, 0.001561516430633, 0.001572792516786633, 0.0016150328324987346, 0.0012622225879640489, 0.0012822064094558925, 0.0013016566469174224, 0.0013253110828427808, 0.0013615178314468801, 0.0013672753309123304, 0.0014087879513501422, 0.0014307191508915913, 0.0014690283419203894, 0.0014975723732589554, 0.0015310555032172716, 0.001564964969537854, 0.001584814514229899, 0.0016522516558236666, 0.0016695146756874769, 0.0017209848115561583, 0.0010936539804409181, 0.0011104691647516255, 0.0011299520645421525, 0.001149610347052481, 0.0011798660740155295, 0.0011834877789316535, -0.0012180712533561228, 0.001235185966170205, 0.0012650205838860226, 0.0012897330585531074, 0.001314951260909776, 0.0013441076570292988, 0.0013567260646748984, 0.001414987413240526, 0.0014282553961591391, 0.0014700810874366777, 0.00100665466146024, 0.001021092316286709, 0.0010392438265376797, 0.0010564071246929843, 0.001084735390328361, 0.0010880831965889336, 0.001119829764209365, 0.001136087660152228, 0.0011646361811603009, 0.0011851066184531333, 0.0012114419188506262, 0.0012351043376467381, 0.0012485553618770849, 0.0012988280639822246, 0.001309405279508703, 0.0013477824956327862, 0.001325569134400243, 0.001354820717448663, 0.0013845536604623742, 0.0014184147066623118, -0.0014669605978256798, 0.0014820134565243213, 0.0015383428117321405, 0.0015709900388268887, 0.0016249550312820381, 0.0016732266340927205, 0.0017270755345527521, 0.0017819416832398885, 0.0018221652613379811, 0.0019149042351733745, 0.0019589830951940107, 0.002041166820981876, 0.001049914441511213, 0.0010695306833592379, 0.0010890970281212508, 0.001111922325393364, 0.0011424092348095946, 0.001148122323619421, 0.001186850254733248, 0.0012085399722806661, 0.0012456339118897064, 0.0012740577277863866, 0.0013051377461435994, 0.0013437995213847341, 0.0013608435012981645, 0.0014261567018563652, 0.001442187100708035, 0.0014950194539597914, 0.00114555470342776, 0.0011690785285024658, 0.0011953819369324209, 0.0012224745957354897, 0.0012660644295784073, 0.0012790470563672606, 0.0013252699223236396, 0.0013590678356097331, 0.0013998231391182873, 0.0014374979624094744, 0.001479309383381721, 0.0015258755232421593, 0.0015638996161407464, 0.0016627081767713022, 0.001696824626273422, 0.0017678261017965403, 0.0012507901376389002, 0.0012819999359381866, 0.001297363796473682, 0.0013318723180596533, 0.001385852791421099, 0.0014046865444862135, 0.0014743521322658013, 0.0015037511646375202, 0.0015632301944286569, 0.0016086222995817566, 0.0016696504645476386, 0.0017420944477222488, -0.0017904681659136575, 0.0018925273282342004, 0.00192874644861548, 0.002018049717598483, 0.001253563775560509, 0.001282420843141195, 0.0013077340349358642, 0.001347245964585091, 0.0014012603142619766, 0.0014300265226613433, -0.0014918195742329803, -0.0015427727090240747, 0.001606053931367674, 0.0016548222695746712, 0.0017040828744978845, 0.0017688632895878444, 0.0018298503029674695, 0.001949157418540235, 0.0019938468062454706, 0.0020996307485186883, 0.0012967738588555982, 0.0013429560480351073, 0.0013816676067445646, -0.0014237371566038857, -0.001477898120439894, 0.0015058505177974486, 0.0015785162494333627, 0.0016363137735252462, -0.0017051362876069971, -0.0017896302948190803, -0.0018847679392540767, 0.0019409331455027597, -0.0020006945038912693, 0.0021727976557029778, -0.002273716859230789, 0.0023396047635150775, 0.0008859252345423805, 0.0009025228352248297, 0.000911962866244815, -0.0009435656796191546, 0.0009737825350755122, -0.0009854183131287981, 0.0010197613779826201, 0.0010449658441497784, 0.0010668650783340354, -0.0011031229007982025, -0.0011449295039070411, -0.0011707300436082097, 0.0012022595830936337, -0.0012449545990435287, -0.0012534255930387156, 0.0013137634698896727, 0.0009119468110754321, 0.0009339717288566699, 0.0009536434452721565, 0.0009827274065103246, 0.0010150434871248233, 0.001018913956270902, -0.0010662324719556522, 0.001068086371752219, 0.0011203838673951186, 0.001151533705969377, 0.0011954824226171228, -0.0012330488985407872, 0.0012865356517785733, -0.001299623727728437, -0.0013382247989299093, -0.0013900490044093006, -0.004867287466420086, 0.00556596746561182, 0.005815442338898269, -0.00690731301959757, 0.008085156694745178, 0.009150634244362169, 0.015804269929404896, 0.021626647423421942, 0.028145487136155216, 0.030391451106210242, 0.032856427720289384, 0.03301799389709809, 0.034801623581296376, 0.03886783120824117, 0.03887771050976683, 0.038342573698221444, 5.679770217925697, 4.251016852930093, 3.5437338099599174, 5.4414444305521865, 2.354307135256903, 1.2375106360535937, -5.926509974424452, 5.105579985273053, 6.220803793450637, -4.7314616203642075, 1.8798682031836234, 5.5737349030831105, 3.201179776528541, 1.753779628606808, 2.276393940822889, 6.4337424520833295, 1.8803472065474063, 1.37622694809517, 2.111479170265415, 2.647643433365171, 5.302181247182872, -0.6554348594430898, 1.9142676280384965, 6.10046184925552, 0.653441307827227, 3.972677408951232, 6.8859818579008865, 3.8038146362032093, 3.0912079807503754, 4.930723984116728, 0.5539812892044788, 6.3521240785694975, 5.734803410944227, 4.631396211030235, 3.5160160070405357, 1.5822384837667445, 4.984363180126585, 5.867871515378414, 3.6173393724017853, 3.2716835155406696, 4.64111910105604, 5.061027693997932, 3.821114354775824, -1.4551420517097016, -0.009680068910482396, 0.996787111078856, 5.594528591259475, 3.7461283499644282, 3.617111888578461, -0.8880106139940267, 1.67043976315479, 2.287238539333411, 5.511579494499221, 0.5714521142705539, 3.802725492932484, 0.47130547011305796, 4.97332952918446, -0.29953259079769096, -0.25001889490376006, 4.838712232639852, 2.5470329143628603, 3.549985645444483, 0.6010744636307404, 5.610669656210048, -0.24548376869623095, 1.3413055441945323, 4.032151127065874, 4.7104087155607095, -1.4537100582360183, 0.563020725313747, -2.07108824150549, 0.997181039950089, 5.681055401024746, 0.33144537477752495, 8.074229094966102, 4.037738446851733, 1.806980350762186, 4.351954781482624, 0.2745523724520492, 1.4089227530864126, 3.057200171651369, 0.5583868988591522, -0.08708818785242389, 3.5151991456138365, 0.2789465162342843, 2.3576112389327166, 1.581734903155837, 6.553041368749906, 5.636391301613432, 2.2812658090447404, 4.988148484777791, 2.1057427438462026, 2.8124362586259126, 1.383133441772288, -1.7975187312537604, 0.5343188007447871, 1.894759150030361, 5.568083164220874, 0.2691081343450627]
params = best_param

# check if loss function works
# 0.014844713200192563
loss(params)
# check if derivative of loss function works
loss'(params)

# run the optimization to find the best parameters
val, x, itr = gradientDescent(loss, loss', params, α=0.0000001, maxItr=500, debug=true,
                              useBigValInc=true, argsArePeriodic=true)
# Set the parameters into the circuits
setparameters!(qc, params[514:end])

################################################################################
# Generate the discretization and check the final results
Δx_param = abs.(params[2:513])
opt_x = zeros(513)
opt_x[1] = abs(params[1])
for i in 1:length(Δx_param)
    opt_x[i+1] = opt_x[i] + Δx_param[i]
end
println("<$(params[1]), $(opt_x[end])>")
python_l2(params[514:end], opt_x)
# <0.6985367664022247, 1.749949851055071>
# 0.007330783173579125

# Print the discretization
# !!!! this is copied to Python file state_preparation
println(opt_x)

# Generate the Python code from the circuits
# !!!! this is copied to Python file state_preparation
println(toPythonQiskitCode(qc))


text:
module HermiteNormalForm

using LinearAlgebra
export hnf!, hnf, ishnf

"""
    ishnf(A)

Check if matrix `A` is in Hermite norm form.
"""
function ishnf(A)
    m, n = size(A)
    for j in 1:n
        for i in 1:min(m, j-1)
            A[i, j] == 0 || return false
        end
        if m >= j
            A[j, j] > 0 || return false
        end
    end
    for i in 1:min(m, n)
        p = A[i, i]
        for j in 1:min(n, i-1)
            0 <= A[i, j] < p || return false
        end
    end
    return true
end

"""
    H, U = hnf(A)

`U` is unimodular and `H` is the Hermite normal form of `A`. If `A` is singular,
then it returns `nothing`.
"""
hnf(A) = hnf!(copy(A))

function hnf!(A)
    m, n = size(A)
    T = eltype(A)
    U = Matrix{T}(I, n, n)
    zz = zero(T)
    for k in 1:m
        # Pivot: A[k, k] should not be 0
        if n >= k && A[k, k] == zz
            pivot = k
            for j in k+1:n
                if A[k, j] != zz
                    pivot = j
                end
            end
            pivot == k && return nothing # rank deficient
            for i in 1:m
                A[i, k], A[i, pivot] = A[i, pivot], A[i, k]
            end
            for i in 1:n
                U[i, k], U[i, pivot] = U[i, pivot], U[i, k]
            end
        end
        # [A11   0 ] k-1 rows
        # [A21  A22]
        # Zero out A[k, k+1:n] === A22[1, 2:end] by multiplying
        # [p  -A[k, j]/d]
        # [q   A[k, k]/d]
        for j in k+1:n
            Akk, Akj = A[k, k], A[k, j]
            d, p, q = gcdx(Akk, Akj)
            Akkd, Akjd = div(Akk, d), div(Akj, d)
            for i in 1:m
                Aik, Aij = A[i, k], A[i, j]
                A[i, k] = Aik * p + Aij * q
                A[i, j] = -Aik * Akjd + Aij * Akkd
            end
            for i in 1:n
                Uik, Uij = U[i, k], U[i, j]
                U[i, k] = Uik * p + Uij * q
                U[i, j] = -Uik * Akjd + Uij * Akkd
            end
        end
        n >= k || continue
        # Ensure the positivity of A[k, k]
        if A[k, k] < zz
            @. A[:, k] = -A[:, k]
            @. U[:, k] = -U[:, k]
        end
        # Minimize A[k, 1:k-1] === A21[1, :]
        for j in 1:k-1
            mul = fld(A[k, j], A[k, k])
            @. A[:, j] -= mul * A[:, k]
            @. U[:, j] -= mul * U[:, k]
        end
    end
    A, U
end

end


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

mutable struct SourceVault <: SwaggerModel
    id::Any # spec type: Union{ Nothing, String } # spec name: id

    function SourceVault(;id=nothing)
        o = new()
        validate_property(SourceVault, Symbol("id"), id)
        setfield!(o, Symbol("id"), id)
        o
    end
end # type SourceVault

const _property_map_SourceVault = Dict{Symbol,Symbol}(Symbol("id")=>Symbol("id"))
const _property_types_SourceVault = Dict{Symbol,String}(Symbol("id")=>"String")
Base.propertynames(::Type{ SourceVault }) = collect(keys(_property_map_SourceVault))
Swagger.property_type(::Type{ SourceVault }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_SourceVault[name]))}
Swagger.field_name(::Type{ SourceVault }, property_name::Symbol) =  _property_map_SourceVault[property_name]

function check_required(o::SourceVault)
    true
end

function validate_property(::Type{ SourceVault }, name::Symbol, val)
end


text:
function aggregate_reports(outfile; rootdir=".", pattern=raw"report.md$")
    report_files = filter(x->x!=outfile && !isnothing(match(Regex(pattern), x)), readdir(rootdir))

    for report in report_files
        report_name = replace(strip(split(lowercase(report), "report")[1], '_'), '_'=>' ')
        @info "Processing \"$report\"..."
        report_name = titlecase(report_name) * " Report"
        open(outfile, "a") do io
            println(io, "\n<details>\n<summary> $report_name </summary>\n")
            open(report, "r") do src_io
                write(io, read(src_io))
            end
            println(io, "\n</details>\n")
        end
    end
    return report_files
end

outfile = "cml_report.md"
report_files = aggregate_reports(outfile)

if !isempty(report_files)
    run(`cml-send-comment $outfile`)
end


text:
function npqreg(y, x, tau, method=IP(); m=50, h=2, xrange=nothing)
	if xrange === nothing
		xrange = collect(range(extrema(x)..., length=m))
	end
	z = copy(x)
	Z = hcat(fill(1, length(x)), z)

	ghat = copy(xrange)
	dghat = copy(xrange)
	for i in 1:length(xrange)
	  z .= x .- xrange[i]
	  Z[:, 2] .= z

	  w = pdf.(Normal(), z./h)

	  wy = w.*y
	  wZ = w.*Z

	  r = qreg_coef(wy, wZ, tau, method)

	  ghat[i] = r[1]
	  dghat[i] = r[2]
	end

	xrange, ghat, dghat
end


text:
using Random

include("stable_bl_model.jl")

function add_perturbations!(state, localgeo)
    FT = eltype(state)
    z = localgeo.coord[3]
    if z <= FT(50) # Add random perturbations to bottom 50m of model
        state.energy.ρe += (rand() - 0.5) * state.energy.ρe / 100
    end
end

function set_clima_parameters(filename)
    eval(:(include($filename)))
end

function main(cl_args)

    surface_flux = cl_args["surface_flux"]

    FT = Float64
    config_type = AtmosLESConfigType
    # DG polynomial order
    N = 4
    # Domain resolution and size
    Δh = FT(20)
    Δv = FT(20)

    resolution = (Δh, Δh, Δv)

    # Prescribe domain parameters
    xmax = FT(100)
    ymax = FT(100)
    zmax = FT(400)

    t0 = FT(0)

    # Required simulation time == 9hours
    timeend = FT(3600 * 0.1)
    CFLmax = FT(0.4)

    C_smag_ = C_smag(param_set) #FT(0.23)

    model = stable_bl_model(
        FT,
        config_type,
        zmax,
        surface_flux;
        turbulence = SmagorinskyLilly{FT}(C_smag_),
    )

    ics = model.problem.init_state_prognostic

    # Assemble configuration
    driver_config = ClimateMachine.AtmosLESConfiguration(
        "StableBoundaryLayer",
        N,
        resolution,
        xmax,
        ymax,
        zmax,
        param_set,
        init_problem!,
        model = model,
    )

    # Choose default IMEX solver
    ode_solver_type = ClimateMachine.ExplicitSolverType()

    solver_config = ClimateMachine.SolverConfiguration(
        t0,
        timeend,
        driver_config,
        ode_solver_type = ode_solver_type,
        init_on_cpu = true,
        Courant_number = CFLmax,
    )
    dgn_config = config_diagnostics(driver_config)

    check_cons = (ClimateMachine.ConservationCheck("ρ", "1mins", FT(0.0001)),)

    result = ClimateMachine.invoke!(
        solver_config;
        diagnostics_config = dgn_config,
        check_cons = check_cons,
        check_euclidean_distance = true,
    )
end

# ArgParse in global scope to modify Clima Parameters
sbl_args = ArgParseSettings(autofix_names = true)
add_arg_group!(sbl_args, "StableBoundaryLayer")
@add_arg_table! sbl_args begin
    "--cparam-file"
    help = "specify CLIMAParameters file"
    arg_type = Union{String, Nothing}
    default = nothing

    "--surface-flux"
    help = "specify surface flux for energy and moisture"
    metavar = "prescribed|bulk|custom_sbl"
    arg_type = String
    default = "custom_sbl"
end

cl_args = ClimateMachine.init(parse_clargs = true, custom_clargs = sbl_args)
if !isnothing(cl_args["cparam_file"])
    filename = cl_args["cparam_file"]
    set_clima_parameters(filename)
end

main(cl_args)


text:
using Plots

# From the inspiring blog post by John D. Cook
# https://www.johndcook.com/blog/2022/05/09/logarithmic-spiral/

ϕ = (1 + 5^0.5)/2
y₀ = 1/(2 + ϕ)
x₀ = (2ϕ + 1) * y₀
θ₀ = atan(-y₀, -x₀)
k = 2 * log(ϕ) / π
a = -x₀ / (exp(k * θ₀) * cos(θ₀))

t = range(-20, θ₀, length = 1000)

x(t) = x₀ + a * exp(k * t) * cos(t)
y(t) = y₀ + a * exp(k * t) * sin(t)

plot(x.(t), y.(t))


text:
###############################################################################
#
# Initialize SoilOpticals
#
###############################################################################
"""
    create_soil_opticals(wls::WaveLengths{FT}) where {FT<:AbstractFloat}

Create [`SoilOpticals`](@ref) struct, given
- `wls` [`WaveLengths`](@ref) type struct
"""
function create_soil_opticals(wls::WaveLengths{FT}) where {FT<:AbstractFloat}
    @unpack iWLF, nWL = wls;

    albedo_SW     = FT(0.2) * ones(FT, nWL);
    albedo_SW_SIF = albedo_SW[iWLF];
    emsvty_SW     = 1 .- albedo_SW;

return SoilOpticals{FT}(albedo_SW,
                        albedo_SW_SIF,
                        emsvty_SW,
                        FT[0.1],
                        FT(290.0))
end

max_stars_repo_name: bbuman/Land

max_issues_repo_name: bbuman/Land

max_forks_repo_name: bbuman/Land


text:
module ShaderAbstractions

using StaticArrays, ColorTypes, FixedPointNumbers, StructArrays, Observables
import GeometryBasics, Tables

include("types.jl")
include("uniforms.jl")
include("context.jl")
include("program.jl")

end # module


text:
"""
   ograv_pk(pk, z, Ωm; kmin = 5e-4, kmax = 3d1)

Comoving density parameter of gravitational binding energy of large-scale structure of the Universe.

*Reference*: Equation (60) of Fukugita & Peebles, ApJ, 616, 643 (2004)
- With the redshift evolution factor of Equation (14) of Chiang, Makiya, Komatsu & Ménard, arXiv:2007.01679

# Arguments
- `pk`(k): a function which returns a matter power spectrum with the argument k being the comoving wavenumber.
    - **pk times k^3 must be dimensionless**. For example, if k is in units of h/Mpc, `pk` must be in units of Mpc^3/h^3.
- `z::Real`: redshift.
- `Ωm::Real`: present-day matter density parameter.

# Optional keyword arguments
- `kmin::Real=5e-4`: minimum wavenumber for integration, ``∫_{kmin}^{kmax} dk P(k)``.
- `kmax::Real=3e1`: maximum wavenumber for integration, ``∫_{kmin}^{kmax} dk P(k)``.
"""
function ograv_pk(pk, z::Real, Ωm::Real; kmin = 5e-4, kmax = 3e1)
   res, err = hquadrature(pk, kmin, kmax)
   halfW = -3 * Ωm * (1 + z) / 2998^2 / 16 / π^2 * res
   Ωgrav = halfW * Ωm
end

"""
   ograv_halo(pk, z, Ωm[, Ωcb=Ωm]; Mmin = 1e11, Mmax = 5e15, virial = false, t10MF = false)

Comoving density parameter of gravitational binding energy of gravitationally collapsed structures (halos).

*Reference*: Equation (17) of Chiang, Makiya, Komatsu & Ménard, arXiv:2007.01679

# Arguments
- `pk`(k): a function which returns a matter power spectrum with the argument k being the comoving wavenumber.
   - `pk` is in units of (Mpc/h)^3 and the argument k is in units of h/Mpc.
- `z::Real`: redshift.
- `Ωm::Real`: present-day total matter density parameter.

# Optional Arguments
- `Ωcb::Real=Ωm`: present-day baryon + cold dark matter density parameter.

# Optional keyword arguments
- `Mmin::Real=1e11`: minimum mass for integration, ``∫_{Mmin}^{Mmax} dM dn/dM Ag GM^2/R``.
- `Mmax::Real=5e15`: maximum mass for integration, ``∫_{Mmin}^{Mmax} dM dn/dM Ag GM^2/R``.
- `virial::Bool=false`: if `true`, use the virial overdensity `Δvir`. If `false` (the default), use `Δm=200`.
- `t10MF::Bool=false`: if `true`, use `tinker10MF` for the halo multiplicity function. If `false` (the default), use `tinker08MF`.
"""
function ograv_halo(
   pk,
   z::Real,
   Ωm::Real,
   Ωcb = Ωm;
   Mmin = 1e11,
   Mmax = 5e15,
   virial = false,
   t10MF = false,
)
   ρc = 2.775e11 # in units of h^2 M⊙/Mpc^3
   GN = 3 / 8π / 2998^2 / ρc # = 4.786e-20 Mpc/M⊙ is Newton's constant
   if virial # Determine Δm used by the mass function
      E2 = Ωm * (1 + z)^3 + 1 - Ωm # flat Universe
      Ωmz = Ωm * (1 + z)^3 / E2
      Δvir = 18 * π^2 + 82 * (Ωmz - 1) - 39 * (Ωmz - 1)^2
      Δm = Δvir / Ωmz
   else
      Δm = 200
   end
   nmass = 100
   func = zeros(nmass)
   lnMh = range(log(Mmin), length = nmass, log(Mmax))
   for i = 1:nmass
      func[i] = dndlnMh(pk, z, Δm, Ωcb, lnMh[i], t10MF = t10MF)
   end
   spl = Spline1D(lnMh, func)
   # %% Compute ρgrav = (1/2) * \int dlnM dn/dlnM Ag*GM^2/R, in units of h^2 M⊙/Mpc^3
   # c(M) relation is due to Equation (2) and "NFW, F, 0-2" of Table 1 in
   # Duffy et al., MNRAS, 390, L64 (2008)
   function dρdlnMh(lnMh)
      if virial # For Δ = Δvir(z)
         RΔh = cbrt(exp(lnMh) * 3 / 4π / (ρc * E2 * Δvir)) # in h^-1 Mpc
         A0, B0, C0 = 7.85, -0.081, -0.71 # Avir, Bvir, Cvir in Table 1 of Duffy et al.
      else # For Δ = Δm = 200
         RΔh = cbrt(exp(lnMh) * 3 / 4π / (ρc * Ωcb * Δm)) / (1 + z) # in h^-1 Mpc
         A0, B0, C0 = 10.14, -0.081, -1.01 # Amean, Bmean, Cmean in Table 1 of Duffy et al.
      end
      c = A0 * (exp(lnMh) / 2e12)^B0 * (1 + z)^C0
      dρdlnMh = 0.5 * spl(lnMh) * GN * exp(2 * lnMh) / RΔh * Ag(c)
   end
   ρgrav, err = hquadrature(dρdlnMh, log(Mmin), log(Mmax))
   redshift = z
   if t10MF
         fnu(x) = tinker10MF(x, redshift, Δm)
         factor, err = hquadrature(fnu, -50, 5)
         ρgrav /= factor
   end
   halfW = -ρgrav / ρc / Ωcb
   Ωgrav = halfW * Ωcb
end

"""
   Ag(c:Real; method = "poly4")

Integral of the squared Fourier transform of an NFW density profile, ``(1/π)∫_0^∞ dx |u(x,c)|^2``.

*Reference*: Equation (18) of Chiang, Makiya, Komatsu & Ménard, arXiv:2007.01679

# Arguments
- `c::Real`: concentration parameter.

# Optional keyword arguments
- `method::String="poly4"`: approximation method.
   - `poly4` (the default) or other values for 4th order polynomial fit.
   - `poly3` for 3rd order polynomial fit.
   - `exact` for numerical integration.
"""
function Ag(c::Real; method = "poly4")
   if method == "poly3" # 3rd order polynomial fit
      y = c - 5
      Ag = 1.0184 + 5.389e-2 * y - 1.514e-3 * y^2 + 3.845e-5 * y^3
   elseif method == "exact" # exact integration
      f = log(1 + c) - c / (1 + c)
      u2(x) =
         (
            sin(x / c) * (sinint((1 / c + 1) * x) - sinint(x / c)) +
            cos(x / c) * (cosint((1 / c + 1) * x) - cosint(x / c)) -
            sin(x) / (1 / c + 1) / x
         )^2
      res, err = quadde(u2, 0, 100)
      Ag = res / π / f^2
   else # 4th order polynomial fit
      y = c - 5
      Ag =
         1.0202 + 5.376e-2 * y - 1.842e-3 * y^2 + 9.489e-5 * y^3 -
         2.352e-6 * y^4
   end
end


text:
using Metatheory
using Test
using Metatheory.Library
using Metatheory.EGraphs
using Metatheory.EGraphs.Schedulers

t = [
    RewriteRule(PatTerm(:call, PatVar(:f), [PatVar(:a), PatVar(:b)], @__MODULE__), PatLiteral(:matched))
]

g = EGraph(:(foo(bar)))
saturate!(g, t)

@test !areequal(g, RewriteRule[], :(foo(bar)), :matched)

addexpr!(g, :(foo(bar, baz)))
saturate!(g, t)

display(g.classes); println()

@test areequal(g, RewriteRule[], :(foo(bar,baz)), :matched)


text:
# test_basic_queries.jl

@testset "Test DBInterface.execute" begin
    con = DBInterface.connect(DuckDB.DB)

    results = DBInterface.execute(con, "SELECT 42 a")

    # iterator
    for row in results
        @test row.a == 42
        @test row[1] == 42
    end

    # convert to DataFrame
    df = DataFrame(results)
    @test names(df) == ["a"]
    @test size(df, 1) == 1
    @test df.a == [42]

    DBInterface.close!(con)
end

@testset "Test numeric data types" begin
    con = DBInterface.connect(DuckDB.DB)

    results = DBInterface.execute(
        con,
        """
SELECT 42::TINYINT a, 42::INT16 b, 42::INT32 c, 42::INT64 d, 42::UINT8 e, 42::UINT16 f, 42::UINT32 g, 42::UINT64 h
UNION ALL
SELECT NULL, NULL, NULL, NULL, NULL, NULL, NULL, NULL
"""
    )

    df = DataFrame(results)

    @test size(df, 1) == 2
    @test isequal(df.a, [42, missing])
    @test isequal(df.b, [42, missing])
    @test isequal(df.c, [42, missing])
    @test isequal(df.d, [42, missing])
    @test isequal(df.e, [42, missing])
    @test isequal(df.f, [42, missing])
    @test isequal(df.g, [42, missing])
    @test isequal(df.h, [42, missing])

    DBInterface.close!(con)
end

@testset "Test strings" begin
    con = DBInterface.connect(DuckDB.DB)

    results = DBInterface.execute(
        con,
        """
SELECT 'hello world' s
UNION ALL
SELECT NULL
UNION ALL
SELECT 'this is a long string'
UNION ALL
SELECT 'obligatory mühleisen'
UNION ALL
SELECT '🦆🍞🦆'
"""
    )

    df = DataFrame(results)
    @test size(df, 1) == 5
    @test isequal(df.s, ["hello world", missing, "this is a long string", "obligatory mühleisen", "🦆🍞🦆"])

    DBInterface.close!(con)
end

@testset "DBInterface.execute errors" begin
    con = DBInterface.connect(DuckDB.DB)

    # parser error
    @test_throws DuckDB.QueryException DBInterface.execute(con, "SELEC")

    # binder error
    @test_throws DuckDB.QueryException DBInterface.execute(con, "SELECT * FROM this_table_does_not_exist")

    # run-time error
    @test_throws DuckDB.QueryException DBInterface.execute(
        con,
        "SELECT i::int FROM (SELECT '42' UNION ALL SELECT 'hello') tbl(i)"
    )

    DBInterface.close!(con)
end

max_stars_repo_name: lokax/duckdb

max_issues_repo_name: lokax/duckdb

max_forks_repo_name: lokax/duckdb


text:
"""
    independent_set(g, reps, MaximalIndependentSet(); parallel=:threads)

Perform [`LightGraphs.independent_set(g, MaximalIndependentSet())`](@ref) `reps` times in parallel
and return the solution with the most vertices.

### Optional Arguements
- `parallel=:threads`: If `parallel=:distributed` then the multiprocessor implementation is
used. This implementation is more efficient if `reps` is large.
"""
independent_set(g::AbstractGraph{T}, reps::Integer, alg::MaximalIndependentSet; parallel=:threads) where T <: Integer =
LightGraphs.Parallel.generate_reduce(g, (g::AbstractGraph{T})->LightGraphs.independent_set(g, alg),
(x::Vector{T}, y::Vector{T})->length(x)>length(y), reps; parallel=parallel)


text:
# Apply boundary conditions to the system

function apply_boundary_conditions_rhs_only(var, b, t)
    return apply_boundary_conditions_lhs_rhs(var, nothing, b, t)
end

function apply_boundary_conditions_lhs_rhs(var, A, b, t)
    if A === nothing
        rhs_only = true
    else
        rhs_only = false
    end
    multivar = typeof(var) <: Array;
    maxvarindex = 0;
    if multivar
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        var_to_dofs = [];
        for vi=1:length(var)
            tmp = dofs_per_node;
            dofs_per_node += var[vi].total_components;
            push!(var_to_dofs, (tmp+1):dofs_per_node);
            maxvarindex = max(maxvarindex,var[vi].index);
        end
    else
        # one variable
        dofs_per_node = var.total_components;
        maxvarindex = var.index;
    end

    # If there is a reason to zero the rhs bdry entries, uncomment this.
    # But this destroys the case NO_BC.
    #b = zero_rhs_bdry_vals(b, dofs_per_node);

    bidcount = length(grid_data.bids); # the number of BIDs
    if !rhs_only
        dirichlet_rows = zeros(0);
        neumann_Is = zeros(Int,0);
        neumann_Js = zeros(Int,0);
        neumann_Vs = zeros(0);
    end
    if dofs_per_node > 1
        if multivar
            dofind = 0;
            for vi=1:length(var)
                for compo=1:var[vi].total_components
                    dofind = dofind + 1;
                    for bid=1:bidcount
                        if prob.bc_type[var[vi].index, bid] == NO_BC
                            # do nothing
                        elseif rhs_only
                            # When doing RHS only, values are simply placed in the vector.
                            b = dirichlet_bc(b, prob.bc_func[var[vi].index, bid][compo], grid_data.bdryface[bid], t, dofind, dofs_per_node, rhs_only=true);
                        elseif prob.bc_type[var[vi].index, bid] == DIRICHLET
                            (tmprows, b) = dirichlet_bc(b, prob.bc_func[var[vi].index, bid][compo], grid_data.bdryface[bid], t, dofind, dofs_per_node, rhs_only=false);
                            append!(dirichlet_rows, tmprows);
                        elseif prob.bc_type[var[vi].index, bid] == NEUMANN
                            (tmpIs, tmpJs, tmpVs, b) = neumann_bc(b, prob.bc_func[var[vi].index, bid][compo], grid_data.bdryface[bid], bid, t, dofind, dofs_per_node);
                            append!(neumann_Is, tmpIs);
                            append!(neumann_Js, tmpJs);
                            append!(neumann_Vs, tmpVs);
                        elseif prob.bc_type[var[vi].index, bid] == ROBIN
                            printerr("Robin BCs not ready.");
                        else
                            printerr("Unsupported boundary condition type: "*prob.bc_type[var[vi].index, bid]);
                        end
                    end
                end
            end
        else
            for d=1:dofs_per_node
                dofind = d;
                for bid=1:bidcount
                    if prob.bc_type[var.index, bid] == NO_BC
                        # do nothing
                    elseif rhs_only
                        # When doing RHS only, values are simply placed in the vector.
                        b = dirichlet_bc(b, prob.bc_func[var.index, bid][d], grid_data.bdryface[bid], t, d, dofs_per_node, rhs_only=true);
                    elseif prob.bc_type[var.index, bid] == DIRICHLET
                        (tmprows, b) = dirichlet_bc(b, prob.bc_func[var.index, bid][d], grid_data.bdryface[bid], t, dofind, dofs_per_node, rhs_only=false);
                        append!(dirichlet_rows, tmprows);
                    elseif prob.bc_type[var.index, bid] == NEUMANN
                        (tmpIs, tmpJs, tmpVs, b) = neumann_bc(b, prob.bc_func[var.index, bid][d], grid_data.bdryface[bid], bid, t, dofind, dofs_per_node);
                        append!(neumann_Is, tmpIs);
                        append!(neumann_Js, tmpJs);
                        append!(neumann_Vs, tmpVs);
                    elseif prob.bc_type[var.index, bid] == ROBIN
                        printerr("Robin BCs not ready.");
                    else
                        printerr("Unsupported boundary condition type: "*prob.bc_type[var.index, bid]);
                    end
                end
            end
        end
    else
        for bid=1:bidcount
            if prob.bc_type[var.index, bid] == NO_BC
                # do nothing
            elseif rhs_only
                # When doing RHS only, values are simply placed in the vector.
                b = dirichlet_bc(b, prob.bc_func[var.index, bid][1], grid_data.bdryface[bid], t, rhs_only=true);
            elseif prob.bc_type[var.index, bid] == DIRICHLET
                #(A, b) = dirichlet_bc(A, b, prob.bc_func[var.index, bid][1], grid_data.bdry[bid], t);
                (tmprows, b) = dirichlet_bc(b, prob.bc_func[var.index, bid][1], grid_data.bdryface[bid], t, rhs_only=false);
                append!(dirichlet_rows, tmprows);
            elseif prob.bc_type[var.index, bid] == NEUMANN
                #(A, b) = neumann_bc(A, b, prob.bc_func[var.index, bid][1], grid_data.bdry[bid], bid, t);
                (tmpIs, tmpJs, tmpVs, b) = neumann_bc(b, prob.bc_func[var.index, bid][1], grid_data.bdryface[bid], bid, t);
                append!(neumann_Is, tmpIs);
                append!(neumann_Js, tmpJs);
                append!(neumann_Vs, tmpVs);
            elseif prob.bc_type[var.index, bid] == ROBIN
                printerr("Robin BCs not ready.");
            else
                printerr("Unsupported boundary condition type: "*prob.bc_type[var.index, bid]);
            end
        end
    end

    if !rhs_only
        # If there are overlaps(like corner nodes), favor dirichlet.
        if length(neumann_Is)>0
            A = insert_sparse_rows(A, neumann_Is, neumann_Js, neumann_Vs);
        end
        # display(A[1:2,1:5]);
        if length(dirichlet_rows)>0
            A = identity_rows(A, dirichlet_rows, length(b));
        end
        # display(A[1:2,1:5]);
    end

    # Reference points
    if size(prob.ref_point,1) >= maxvarindex
        if multivar
            posind = zeros(Int,0);
            vals = zeros(0);
            for vi=1:length(var)
                if prob.ref_point[var[vi].index,1]
                    eii = prob.ref_point[var[vi].index, 2];
                    tmp = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + var_to_dofs[vi][1];
                    if length(prob.ref_point[var[vi].index, 3]) > 1
                        tmp = tmp:(tmp+length(prob.ref_point[var[vi].index, 3])-1);
                    end
                    posind = [posind; tmp];
                    vals = [vals; prob.ref_point[var[vi].index, 3]];
                end
            end
            if length(vals) > 0
                if !rhs_only A = identity_rows(A, posind, length(b)); end
                b[posind] = vals;
            end

        else
            if prob.ref_point[var.index,1]
                eii = prob.ref_point[var.index, 2];
                posind = (grid_data.glbvertex[eii[1], eii[2]] - 1)*dofs_per_node + 1;
                if length(prob.ref_point[var.index, 3]) > 1
                    posind = posind:(posind+length(prob.ref_point[var[vi].index, 3])-1);
                else
                    posind = [posind];
                end
                if !rhs_only A = identity_rows(A, posind, length(b)); end
                b[posind] = prob.ref_point[var.index, 3];
            end
        end
    end

    if rhs_only
        return b;
    else
        return (A, b);
    end
end

# This is just to speed things up (huge improvement)
function is_in_rows(k, rows)
    # rows are sorted
    b = 1;
    e = length(rows);
    c = Int(round((e+1)/2));
    while e-b > 1
        if k == rows[c]
            return true;
        elseif k > rows[c]
            b=c;
        else
            e=c;
        end
        c = Int(round((b+e)/2));
    end
    return k==rows[b] || k==rows[e];
end

# zeros the rows and puts a 1 on the diagonal
function identity_rows(A, rows, N)
    if length(rows) < 1
        return A;
    end
    if issparse(A)
        (I, J, V) = findnz(A);
        eN = length(I);
        newI = zeros(Int,eN);
        newJ = zeros(Int,eN);
        newV = zeros(eN);
        newind = 0;
        (M,N) = size(A);

        sort!(rows);

        # remove duplicates
        newrows = zeros(Int, length(rows));
        newcount = 1;
        newrows[1] = rows[1];
        for i=2:length(rows)
            if rows[i] != rows[i-1] # because sorted
                newcount += 1;
                newrows[newcount] = rows[i];
            end
        end
        rows = newrows[1:newcount];

        # Remove bdry row elements
        for k in 1:length(I)
            isbdry = is_in_rows(I[k], rows);
            if !isbdry
                newind = newind+1;
                newI[newind] = I[k];
                newJ[newind] = J[k];
                newV[newind] = V[k];
            end
        end

        # Add diogonal 1s
        rn = length(rows);
        if newind + rn <= eN
            for i=1:rn
                newind = newind+1;
                newI[newind] = rows[i];
                newJ[newind] = rows[i];
                newV[newind] = 1.0;
            end
            newI = newI[1:newind];
            newJ = newJ[1:newind];
            newV = newV[1:newind];
        else
            toadd = rn - (eN-newind);
            for i=1:(eN-newind)
                newind = newind+1;
                newI[newind] = rows[i];
                newJ[newind] = rows[i];
                newV[newind] = 1.0;
            end
            append!(newI, rows[toadd:rn]);
            append!(newJ, rows[toadd:rn]);
            append!(newV, ones(toadd));
        end

        return sparse(newI,newJ,newV, M,N);

    else
        for i=1:length(rows)
            A[rows[i],:] = zeros(N);
            A[rows[i],rows[i]] = 1;
        end
        return A;
    end
end

# Inserts the rows of S in A
function insert_rows(A, S, rows, N)
    if issparse(A)
        (I, J, V) = findnz(A);
        (M,N) = size(A);
        # determine which elements to remove
        toskip = zeros(Bool, length(I));
        includen = 0;
        for k = 1:length(I)
            for i=1:length(rows)
                if I[k] == rows[i]
                    toskip[k] = 1;
                end
            end
            if toskip[k] == 0
                includen = includen + 1;
            end
        end

        # put the remaining elements in a new matrix
        newI = zeros(Int, includen);
        newJ = zeros(Int, includen);
        newV = zeros(includen);
        ind = 1;
        for k = 1:length(I)
            if toskip[k] == 0
                newI[ind] = I[k];
                newJ[ind] = J[k];
                newV[ind] = V[k];
                ind = ind + 1;
            end
        end

        # Do something similar to extract elements of S
        (I, J, V) = findnz(S);

        toskip = ones(Int, length(I));
        includen = 0;
        for k = 1:length(I)
            for i=1:length(rows)
                if I[k] == rows[i]
                    toskip[k] = 0;
                end
            end
            if toskip[k] == 0
                includen = includen + 1;
            end
        end
        newI2 = zeros(Int, includen);
        newJ2 = zeros(Int, includen);
        newV2 = zeros(includen);
        ind = 1;
        for k = 1:length(I)
            if toskip[k] == 0
                newI2[ind] = I[k];
                newJ2[ind] = J[k];
                newV2[ind] = V[k];
                ind = ind + 1;
            end
        end

        append!(newI, newI2);
        append!(newJ, newJ2);
        append!(newV, newV2);
        return sparse(newI,newJ,newV, M, N);
    else
        for i=1:length(rows)
            A[rows[i],:] = S[rows[i],:];
        end
        return A;
    end
end

function insert_sparse_rows(A, SI, SJ, SV)
    (I, J, V) = findnz(A);
    (M,N) = size(A);

    # since we need to sort the three arrays in sync, just get the sorted permutation
    sorted_order = sortperm(SI);

    SI = SI[sorted_order];
    SJ = SJ[sorted_order];
    SV = SV[sorted_order];

    # Zero existing elements in SI rows
    for k in 1:length(I)
        if is_in_rows(I[k], SI)
            V[k] = 0;
        end
    end

    # append S values
    append!(I, SI);
    append!(J, SJ);
    append!(V, SV);

    # Form sparse matrix
    return sparse(I, J, V, M, N);
end

function zero_rhs_bdry_vals(rhs, dofs)
    for bid=1:length(grid_data.bids)
        for fi=1:length(grid_data.bdryface[bid])
            face_nodes = grid_data.face2glb[:, 1, grid_data.bdryface[bid][fi]];
            for ni=1:length(face_nodes)
                rhs_ind = ((face_nodes[ni]-1)*dofs + 1):(face_nodes[ni]*dofs);
                rhs[rhs_ind] .= 0;
            end
        end
    end
    return rhs;
end

function evaluate_at_nodes(val, nodes, t=0)
    N = length(nodes);
    if typeof(val) <: Number
        result = fill(val, N);

    elseif typeof(val) == Coefficient && typeof(val.value[1]) == GenFunction
        result = zeros(N);
        if config.dimension == 1
            for i=1:N
                result[i]=val.value[1].func(grid_data.allnodes[1,nodes[i]],0,0,t,nodes[i]);
            end
        elseif config.dimension == 2
            for i=1:N
                result[i]=val.value[1].func(grid_data.allnodes[1,nodes[i]],grid_data.allnodes[2,nodes[i]],0,t,nodes[i]);
            end
        else
            for i=1:N
                result[i]=val.value[1].func(grid_data.allnodes[1,nodes[i]],grid_data.allnodes[2,nodes[i]],grid_data.allnodes[3,nodes[i]],t,nodes[i]);
            end
        end

    elseif typeof(val) == GenFunction
        result = zeros(N);
        if config.dimension == 1
            for i=1:N
                result[i]=val.func(grid_data.allnodes[1,nodes[i]],0,0,t,nodes[i]);
            end
        elseif config.dimension == 2
            for i=1:N
                result[i]=val.func(grid_data.allnodes[1,nodes[i]],grid_data.allnodes[2,nodes[i]],0,t,nodes[i]);
            end
        else
            for i=1:N
                result[i]=val.func(grid_data.allnodes[1,nodes[i]],grid_data.allnodes[2,nodes[i]],grid_data.allnodes[3,nodes[i]],t,nodes[i]);
            end
        end
    end
    return result;
end

function dirichlet_bc(b, val, bdryface, t=0, dofind = 1, totaldofs = 1; rhs_only=false)
    # Assuming every face has the same number of nodes
    nodes_per_face = size(grid_data.face2glb,1);
    nnodes = length(bdryface)*nodes_per_face;
    if !rhs_only
        bdry_rows = zeros(Int, nnodes);
    end

    for fi = 1:length(bdryface)
        # global indices for the nodes on the face
        face_nodes = grid_data.face2glb[:, 1, bdryface[fi]];

        # offset for multi dof
        if totaldofs > 1
            offsetglb = (face_nodes.-1) .* totaldofs .+ dofind;
        else
            offsetglb = face_nodes;
        end

        if !rhs_only
            # add to the list of bdry rows
            bdry_rows[((fi-1)*nodes_per_face+1):(fi*nodes_per_face)] = offsetglb;
        end

        b[offsetglb] = evaluate_at_nodes(val, face_nodes, t);

    end

    if !rhs_only
        return (bdry_rows, b);
    else
        return b;
    end
end

function neumann_bc(b, val, bdryface, bid, t=0, dofind = 1, totaldofs = 1)
    # Assuming every element has the same number of nodes
    nodes_per_element = size(grid_data.loc2glb,1);
    nodes_per_face = size(grid_data.face2glb,1);

    # To deal with the possibility of overlapping nodes between elements or other bids
    # limit it to nodes assigned to this BID and assign each node to one of the bdryfaces.
    # Need to form a list of valid nodes for each face in bdryface.
    bidnodes = grid_data.bdry[bid];
    node_available = fill(true, length(bidnodes));
    valid_face_mask = zeros(Bool, nodes_per_face, length(bdryface));
    for fi=1:length(bdryface)
        face_nodes = grid_data.face2glb[:, 1, bdryface[fi]];
        for ni=1:length(face_nodes)
            bnode_ind = indexin([face_nodes[ni]], bidnodes);
            if !(bnode_ind[1] === nothing)
                if node_available[bnode_ind[1]]
                    # This is a valid node for this face. Reserve it
                    node_available[bnode_ind[1]] = false;
                    valid_face_mask[ni,fi] = true;
                end
            end
        end
    end
    total_entries = length(bidnodes) * nodes_per_element;

    # These will be sparse rows to insert into A
    S_I = zeros(Int, total_entries);
    S_J = zeros(Int, total_entries);
    S_V = zeros(total_entries);
    next_S_ind = 1;

    # This can be precomputed for uniform grid meshes
    if config.mesh_type == UNIFORM_GRID && config.geometry == SQUARE
        precomputed = true;
        glb = grid_data.loc2glb[:,1];
        detJ = geo_factors.detJ[1];
        J = geo_factors.J[1]
        if config.dimension == 1
            RD1 = zeros(size(refel.Ddr));
            # Multiply rows of Qr and Ddr by J.rx
            for i=1:size(RD1,2) # loop over columns
                RD1[:,i] = J.rx .* refel.Ddr[:,i];
            end

        elseif config.dimension == 2
            RD1 = zeros(size(refel.Ddr));
            RD2 = zeros(size(refel.Ddr));
            for i=1:size(RD1,2)
                RD1[:,i] = J.rx .* refel.Ddr[:,i] + J.sx .* refel.Dds[:,i];
                RD2[:,i] = J.ry .* refel.Ddr[:,i] + J.sy .* refel.Dds[:,i];
            end

        elseif config.dimension == 3
            RD1 = zeros(size(refel.Ddr));
            RD2 = zeros(size(refel.Ddr));
            RD3 = zeros(size(refel.Ddr));
            for i=1:size(RD1,2)
                RD1[:,i] = J.rx .* refel.Ddr[:,i] + J.sx .* refel.Dds[:,i] + J.tx .* refel.Ddt[:,i];
                RD2[:,i] = J.ry .* refel.Ddr[:,i] + J.sy .* refel.Dds[:,i] + J.ty .* refel.Ddt[:,i];
                RD3[:,i] = J.rz .* refel.Ddr[:,i] + J.sz .* refel.Dds[:,i] + J.tz .* refel.Ddt[:,i];
            end
        end
    else
        precomputed = false;
    end

    for fi = 1:length(bdryface)
        # find the relevant element
        e = grid_data.face2element[1,bdryface[fi]] # >0 ? grid_data.face2element[1,bdryface[fi]] : grid_data.face2element[2,bdryface[fi]];
        glb = grid_data.loc2glb[:,e];       # global indices of this element's nodes
        # offset for multi dof
        if totaldofs > 1
            elglb = (glb.-1) .* totaldofs .+ dofind;
        else
            elglb = glb;
        end

        # Local indices for the nodes on the face
        flocal = refel.face2local[grid_data.faceRefelInd[1,bdryface[fi]]];

        # global indices for the nodes on the face
        face_nodes = grid_data.face2glb[:, 1, bdryface[fi]];

        # remove nodes that do not belong to this face and bid
        newflocal = zeros(Int, length(flocal));
        newface_nodes = zeros(Int, length(face_nodes));
        node_count = 0;
        for ni=1:length(face_nodes)
            if valid_face_mask[ni, fi]
                node_count += 1;
                newflocal[node_count] = flocal[ni];
                newface_nodes[node_count] = face_nodes[ni];
            end
        end
        face_nodes = newface_nodes[1:node_count];
        flocal = newflocal[1:node_count];

        # offset for multi dof
        if totaldofs > 1
            offsetface = (face_nodes.-1) .* totaldofs .+ dofind;
        else
            offsetface = face_nodes;
        end

        # Build diff matrices if needed
        if !precomputed
            detJ = geo_factors.detJ[e];
            J = geo_factors.J[e]
            if config.dimension == 1
                RD1 = zeros(size(refel.Ddr));
                # Multiply rows of Qr and Ddr by J.rx
                for i=1:size(RD1,2) # loop over columns
                    RD1[:,i] = J.rx .* refel.Ddr[:,i];
                end

            elseif config.dimension == 2
                RD1 = zeros(size(refel.Ddr));
                RD2 = zeros(size(refel.Ddr));
                for i=1:size(RD1,2)
                    RD1[:,i] = J.rx .* refel.Ddr[:,i] + J.sx .* refel.Dds[:,i];
                    RD2[:,i] = J.ry .* refel.Ddr[:,i] + J.sy .* refel.Dds[:,i];
                end

            elseif config.dimension == 3
                RD1 = zeros(size(refel.Ddr));
                RD2 = zeros(size(refel.Ddr));
                RD3 = zeros(size(refel.Ddr));
                for i=1:size(RD1,2)
                    RD1[:,i] = J.rx .* refel.Ddr[:,i] + J.sx .* refel.Dds[:,i] + J.tx .* refel.Ddt[:,i];
                    RD2[:,i] = J.ry .* refel.Ddr[:,i] + J.sy .* refel.Dds[:,i] + J.ty .* refel.Ddt[:,i];
                    RD3[:,i] = J.rz .* refel.Ddr[:,i] + J.sz .* refel.Dds[:,i] + J.tz .* refel.Ddt[:,i];
                end
            end
        end

        # Add to S matrix
        for i=1:length(elglb)
            irange = (next_S_ind + (i-1)*node_count):(next_S_ind + i*node_count - 1);
            S_I[irange] = elglb[flocal];
            S_J[irange] .= elglb[i];
        end
        if config.dimension == 1
            S1_V= RD1[flocal,:][:];

        elseif config.dimension == 2
            S1_V= RD1[flocal,:][:];
            S2_V= RD2[flocal,:][:];

        elseif config.dimension == 3
            S1_V= RD1[flocal,:][:];
            S2_V= RD2[flocal,:][:];
            S3_V= RD3[flocal,:][:];
        end

        # Add the right components of S1,S2,S3 according to normal vector
        ind_range = next_S_ind:(next_S_ind + node_count * nodes_per_element - 1);
        norm = grid_data.facenormals[:, bdryface[fi]];
        # println("bid "*string(bid)*" normal "*string(norm));
        if config.dimension == 1
            S_V[ind_range] = norm[1] .* S1_V;
        elseif config.dimension == 2
            S_V[ind_range] = norm[1] .* S1_V + norm[2] .* S2_V;
        elseif config.dimension == 3
            S_V[ind_range] = norm[1] .* S1_V + norm[2] .* S2_V + norm[3] .* S3_V;
        end

        # Do the rhs vector
        b[offsetface] = evaluate_at_nodes(val, face_nodes, t);

        # Increase S index
        next_S_ind += node_count * nodes_per_element;
    end

    if next_S_ind <= length(S_I)
        # This should not happen
        log_entry("Unexpected entry count when doing Neumann BC. expected "*string(length(S_I))*" found "*string(next_S_ind-1));
        S_I = S_I[1:next_S_ind-1];
        S_J = S_J[1:next_S_ind-1];
        S_V = S_V[1:next_S_ind-1];
    end

    # println("S for "*string(bid)*": ")
    # println(S_I);
    # println(S_J);
    # println(S_V);

    return (S_I, S_J, S_V, b);
end

# function neumann_bc_rhs_only(b, val, bdryface, bid, t=0, dofind = 1, totaldofs = 1)
#     return dirichlet_bc(b, val, bdryface, t, dofind, totaldofs, rhs_only=true); # how convenient
# end


text:
module SSA

export ssa_decompose, ssa_reconstruct

using LinearAlgebra
using Distributed
using SharedArrays

mutable struct SSA_Info
   eig_vals::Array{<:AbstractFloat, 1}
   eig_vecs::Array{<:AbstractFloat, 2}
   X::Array{<:AbstractFloat, 2}
   C::Array{<:AbstractFloat, 2}
   N::Integer
   D::Integer
   J::Integer
   M::Integer
   SSA_Info(;eig_vals, eig_vecs, X,
            C=nothing, N, D, J, M) = new(eig_vals, eig_vecs, X, C, N, D, J, M)
end

"""
Singular spectrum analysis eigendecomposition with the Broomhead–King approach
"""
function ssa_decompose(x::Array{float_type, dim},
                       M::Integer) where {float_type<:AbstractFloat} where dim
   if (dim == 1)
      # Single-channel SSA
      N = length(x)
      D = 1
      J = 1
   elseif (dim == 2)
      # Multi-channel SSA
      N, D = size(x)
      J = 1
   elseif (dim == 3)
      # Multi-channel SSA with multiple non-contiguous samples of a series
      N, D, J = size(x)
   else
      throw(ArgumentError("x must be of 1, 2, or 3 dimensions"))
   end

   if M > N
      throw(ArgumentError("M cannot be greater than N"))
   end

   N′ = N-M+1
   X = zeros(float_type, N′*J, D*M)

   for j = 1:J
      for d = 1:D
         for i = 1:N′
             X[(j-1)*N′+i, 1+M*(d-1):M*d] = x[i:i+M-1, d, j]
         end
      end
   end

   if N′*J >= D*M
      C = X'*X/(N′*J)
      eig_vals, eig_vecs = eigen(C)
   else
      # Use PCA transpose trick; see section A2 of Ghil et al. (2002)
      C = X*X'/(N′*J)
      eig_vals, eig_vecs = eigen(C)

      eig_vecs = X'*eig_vecs

      # Normalize eigenvectors
      eig_vecs = eig_vecs./mapslices(norm, eig_vecs, dims=1)
   end

   eig_vals = reverse(eig_vals)
   eig_vecs = reverse(eig_vecs, dims=2)

   return SSA_Info(eig_vals=eig_vals, eig_vecs=eig_vecs, X=X, C=C, N=N, D=D,
                   J=J, M=M)
end

"""
Reconstruct the specified modes
"""
function ssa_reconstruct(ssa_info, modes; sum_modes=false)
   eig_vecs = ssa_info.eig_vecs
   X = ssa_info.X
   M = ssa_info.M
   N = ssa_info.N
   D = ssa_info.D
   J = ssa_info.J
   r = SharedArray{eltype(X), 4}((length(modes), N, D, J))

   for (i_k, k) in enumerate(modes)
      ek = reshape(eig_vecs[:, k], M, D)
      for j = 1:J
         A = X[1+(j-1)*(N-M+1):j*(N-M+1), :]*eig_vecs
         @sync @distributed for n = 1:N
            if 1 <= n <= M - 1
               M_n = n
               L_n = 1
               U_n = n
            elseif M <= n <= N - M + 1
               M_n = M
               L_n = 1
               U_n = M
            elseif N - M + 2 <= n <= N
               M_n = N - n + 1
               L_n = n - N + M
               U_n = M
            end
            for d=1:D
               r[i_k, n, d, j] = 1/M_n*sum([A[n - m + 1, k]*ek[m, d] for m=L_n:U_n])
            end
         end
      end
   end

   if J == 1
      r = r[:, :, :, 1]

      if D == 1
         r = r[:, :, 1]
      end
   end

   if sum_modes
      r = copy(selectdim(sum(r, dims=1), 1, 1))
   end

   return r
end
end


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    HierarchicalPartitioner(first, second)

A partitioning method in which a `first` partition is applied
and then a `second` partition is applied to each subset of the
`first`.
"""
struct HierarchicalPartitioner <: AbstractPartitioner
  first::AbstractPartitioner
  second::AbstractPartitioner
end

function partition(object::AbstractSpatialObject,
                   partitioner::HierarchicalPartitioner)
  result = Vector{Vector{Int}}()

  # use first partition method
  p = partition(object, partitioner.first)

  # use second method to partition the first
  s = subsets(p)
  for (i, d) in Iterators.enumerate(p)
    q = partition(d, partitioner.second)

    for js in subsets(q)
      push!(result, s[i][js])
    end
  end

  SpatialPartition(object, result)
end

→(first, second) = HierarchicalPartitioner(first, second)


text:
julia> # Load DataFrame package

julia> require("DataFrames")

julia> using DataFrames

julia>

julia> # Load a very simple lm() function

julia> require(Pkg.dir("DataFrames", "demo", "lm.jl"))

julia>

julia> # Load a CSV file into a DataFrame

julia> df = read_table(Pkg.dir("DataFrames", "demo", "toy_example.csv"))
6x3 DataFrame:
          A       B   C
[1,]    2.5   "One" 3.0
[2,]    3.6   "One" 5.0
[3,]    3.5   "Two" 3.0
[4,]    4.5   "Two" 5.0
[5,]    4.5 "Three" 3.0
[6,]    5.5 "Three" 5.0

julia>

julia> # Run a simple linear model that predicts A using B and C

julia> lm_fit = lm(:(A ~ B + C), df)
ArgumentError("float64(String): invalid number format")
 in float64 at string.jl:1209
 in float at string.jl:1222
 in model_matrix at /Users/johnmyleswhite/.julia/DataFrames/src/formula.jl:141
 in lm at /Users/johnmyleswhite/.julia/DataFrames/demo/lm.jl:66

julia>

julia> # Print out a summary of the results

julia> print(lm_fit)
lm_fit not defined

julia>

julia> #

julia> # Behind The Scenes

julia> #

julia>

julia> # Generate a Formula object

julia> f = Formula(:(A ~ B + C))
Formula([:A],[:(+(B,C))])

julia>

julia> # Generate a ModelFrame object

julia> mf = model_frame(f, df)
ModelFrame(6x3 DataFrame:
          A       B   C
[1,]    2.5   "One" 3.0
[2,]    3.6   "One" 5.0
[3,]    3.5   "Two" 3.0
[4,]    4.5   "Two" 5.0
[5,]    4.5 "Three" 3.0
[6,]    5.5 "Three" 5.0
,[1],Formula([:A],[:(+(B,C))]))

julia>

julia> # Generate a ModelMatrix object

julia> mm = model_matrix(mf)
ArgumentError("float64(String): invalid number format")
 in float64 at string.jl:1209
 in float at string.jl:1222
 in model_matrix at /Users/johnmyleswhite/.julia/DataFrames/src/formula.jl:141

julia>


text:
# Boundary condition functions

for f in [:prescribed_surface_jump!,:prescribed_surface_average!]
  @eval function $f(dfb,t::Real,sys::ILMSystem)
          @unpack base_cache, bc, phys_params, motions = sys
          $f(dfb,t,base_cache,bc,phys_params,motions)
        end
  @eval function $f(dfb,sys::ILMSystem)
          @unpack base_cache, bc, phys_params = sys
          $f(dfb,base_cache,bc,phys_params)
        end
end

function prescribed_surface_jump!(dfb,t::Real,base_cache::BasicILMCache{N},bc,phys_params,motions) where N
    dfb .= bc["exterior"](t,base_cache,phys_params,motions)
    dfb .-= bc["interior"](t,base_cache,phys_params,motions)
    return dfb
end

function prescribed_surface_jump!(dfb,base_cache::BasicILMCache{N},bc,phys_params) where N
    dfb .= bc["exterior"](base_cache,phys_params)
    dfb .-= bc["interior"](base_cache,phys_params)
    return dfb
end

function prescribed_surface_jump!(dfb,t::Real,base_cache::BasicILMCache{0},bc,phys_params,motions)
    fill!(dfb,0.0)
    return dfb
end

function prescribed_surface_jump!(dfb,base_cache::BasicILMCache{0},bc,phys_params)
    fill!(dfb,0.0)
    return dfb
end

function prescribed_surface_average!(fb,t::Real,base_cache::BasicILMCache{N},bc,phys_params,motions) where N
    fb .= 0.5*bc["exterior"](t,base_cache,phys_params,motions)
    fb .+= 0.5*bc["interior"](t,base_cache,phys_params,motions)
    return fb
end

function prescribed_surface_average!(fb,t::Real,base_cache::BasicILMCache{0},bc,phys_params,motions)
    fill!(fb,0.0)
    return fb
end

function prescribed_surface_average!(fb,base_cache::BasicILMCache{N},bc,phys_params) where N
    fb .= 0.5*bc["exterior"](base_cache,phys_params)
    fb .+= 0.5*bc["interior"](base_cache,phys_params)
    return fb
end

function prescribed_surface_average!(fb,base_cache::BasicILMCache{0},bc,phys_params)
    fill!(fb,0.0)
    return fb
end