text:
using SX, Reachability, MathematicalSystems

SUITE["Build"] = BenchmarkGroup()

file = @relpath "Building_more_decimals.xml"
H = readsxmodel(file, ST=ConstrainedLinearControlContinuousSystem)

# ================
# BLDF01 - BDS01
# ================

n = size(H.modes[1].A, 1)-1
A = H.modes[1].A[1:n, 1:n]
B = eye(n)
U = Hyperrectangle(low=[0.8], high=[1.0])
X, Uin = nothing, ConstantInput(H.modes[1].B[1:n, 1] * U)
S = ConstrainedLinearControlContinuousSystem(A, B, X, Uin)
initially = " x1 >= 0.0002000 & x1 <= 0.0002500 & x2 >= 0.0002000 & x2 <= 0.0002500 & x3 >= 0.0002000 & x3 <= 0.0002500 & x4 >= 0.0002000 & x4 <= 0.0002500 & x5 >= 0.0002000 & x5 <= 0.0002500 & x6 >= 0.0002000 & x6 <= 0.0002500 & x7 >= 0.0002000 & x7 <= 0.0002500 & x8 >= 0.0002000 & x8 <= 0.0002500 & x9 >= 0.0002000 & x9 <= 0.0002500 & x10 >= 0.0002000 & x10 <= 0.0002500 & x11 >= 0.0000000 & x11 <= 0.0000000 & x12 >= 0.0000000 & x12 <= 0.0000000 & x13 >= 0.0000000 & x13 <= 0.0000000 & x14 >= 0.0000000 & x14 <= 0.0000000 & x15 >= 0.0000000 & x15 <= 0.0000000 & x16 >= 0.0000000 & x16 <= 0.0000000 & x17 >= 0.0000000 & x17 <= 0.0000000 & x18 >= 0.0000000 & x18 <= 0.0000000 & x19 >= 0.0000000 & x19 <= 0.0000000 & x20 >= 0.0000000 & x20 <= 0.0000000 & x21 >= 0.0000000 & x21 <= 0.0000000 & x22 >= 0.0000000 & x22 <= 0.0000000 & x23 >= 0.0000000 & x23 <= 0.0000000 & x24 >= 0.0000000 & x24 <= 0.0000000 & x25 >= -0.0001000 & x25 <= 0.0001000 & x26 >= 0.0000000 & x26 <= 0.0000000 & x27 >= 0.0000000 & x27 <= 0.0000000 & x28 >= 0.0000000 & x28 <= 0.0000000 & x29 >= 0.0000000 & x29 <= 0.0000000 & x30 >= 0.0000000 & x30 <= 0.0000000 & x31 >= 0.0000000 & x31 <= 0.0000000 & x32 >= 0.0000000 & x32 <= 0.0000000 & x33 >= 0.0000000 & x33 <= 0.0000000 & x34 >= 0.0000000 & x34 <= 0.0000000 & x35 >= 0.0000000 & x35 <= 0.0000000 & x36 >= 0.0000000 & x36 <= 0.0000000 & x37 >= 0.0000000 & x37 <= 0.0000000 & x38 >= 0.0000000 & x38 <= 0.0000000 & x39 >= 0.0000000 & x39 <= 0.0000000 & x40 >= 0.0000000 & x40 <= 0.0000000 & x41 >= 0.0000000 & x41 <= 0.0000000 & x42 >= 0.0000000 & x42 <= 0.0000000 & x43 >= 0.0000000 & x43 <= 0.0000000 & x44 >= 0.0000000 & x44 <= 0.0000000 & x45 >= 0.0000000 & x45 <= 0.0000000 & x46 >= 0.0000000 & x46 <= 0.0000000 & x47 >= 0.0000000 & x47 <= 0.0000000 & x48 >= 0.0000000 & x48 <= 0.0000000 & t==0"
f = i -> SX.parse_sxmath(initially)[i].args[3]
X0 = Hyperrectangle(low=[f(i) for i in 1:2:96], high=[f(i) for i in 2:2:96])
problem = InitialValueProblem(S, X0)

δ_max = 0.0009
time_horizon = 20.0
pBDS01 = LinearConstraintProperty(sparsevec([25], [1.0], 48), 0.0051) # x25 <= 0.0051

sol = solve(problem, :property => pBDS01, :T=>time_horizon, :δ=>δ_max, :vars=>[25], :mode=>"check")
@assert sol.satisfied
SUITE["Build"]["BLDF01-BDS01", "dense"] = @benchmarkable solve($problem, :property => $pBDS01, :T=>$time_horizon, :δ=>$δ_max, :vars=>[25], :mode=>"check")

sol = solve(problem, :approx_model=>"nobloating", :property=>pBDS01, :T=>time_horizon, :δ=>δ_max, :vars=>[25], :mode=>"check")
@assert sol.satisfied
SUITE["Build"]["BLDF01-BDS01", "discrete"] = @benchmarkable solve($problem, :approx_model=>"nobloating", :property => $pBDS01, :T=>$time_horizon, :δ=>$δ_max, :vars=>[25], :mode=>"check")

# ================
# BLDC01 - BDS01
# ================

n = size(H.modes[1].A, 1)-1
A = H.modes[1].A[1:n, 1:n]
A = Reachability.add_dimension(A) # add an extra zero row and column
S = LinearContinuousSystem(A)
X0 = X0 * U
problemConst = InitialValueProblem(S, X0);

pBLDC01 = LinearConstraintProperty(sparsevec([25], [1.0], n+1), 0.0051) # x25 <= 0.0051

sol = solve(problemConst, :T=>time_horizon, :δ=>δ_max, :vars=>[25], :mode=>"check", :property => pBLDC01, :assume_homogeneous=>true)
@assert sol.satisfied
SUITE["Build"]["BLDC01-BDS01", "dense"] = @benchmarkable solve($problemConst, :T=>$time_horizon, :δ=>$δ_max, :vars=>[25], :mode=>"check", :property => $pBLDC01, :assume_homogeneous=>true)

sol = solve(problemConst, :T=>time_horizon, :δ=>δ_max, :vars=>[25], :mode=>"check", :property => pBLDC01, :assume_homogeneous=>true)
@assert sol.satisfied
SUITE["Build"]["BLDC01-BDS01", "discrete"] = @benchmarkable solve($problemConst, :approx_model=>"nobloating", :T=>$time_horizon, :δ=>$δ_max, :vars=>[25], :mode=>"check", :property => $pBLDC01, :assume_homogeneous=>true)


text:
function plotlabel(string::String, fontSize, plotPadding)
    plot(x=[0],y=[0],label=[string], Geom.label(position=:centered),
         Guide.xlabel(""), Guide.ylabel(""), Guide.xticks(label=false), Guide.yticks(label=false),
         Theme(point_label_font_size=fontSize, plot_padding=[plotPadding], grid_color=RGBA(1,1,1,0)))
end

function pairwisescatterplot(X::AbstractMatrix, groupBelow, groupColorsBelow, groupAbove, groupColorsAbove;
                             padding = 10.0mm/size(X,2), point_size=0.9mm, highlight_width=point_size/3)
    N = size(X,2)
    P = Matrix{Plot}(undef,N,N)
    for j=1:N
        for i=1:N
            if i==j
                P[i,j] = plotlabel(string(i), 120mm/N, padding)
            else
                group,groupColors = i>j ? (groupBelow,groupColorsBelow) : (groupAbove,groupColorsAbove)
                x,y = X[:,i],-X[:,j]
                df = DataFrame(x=x, y=y, group=group)
                P[i,j] = plot(df, x=:x, y=:y, color=:group, groupColors,
                              Guide.xlabel(""), Guide.ylabel(""),
                              Guide.xticks(ticks=[extrema(x)...],label=false), Guide.yticks(ticks=[extrema(y)...],label=false),
                              Theme(key_position=:none,  # disable legend
                                    plot_padding=[padding],
                                    point_size=point_size,
                                    highlight_width=highlight_width))
            end
        end
    end
    P
end
pairwisescatterplot(X, group, groupColors) = pairwisescatterplot(X, group, groupColors, group, groupColors)

ext: jl


text:
using GaussianProcesses
using GaussianProcesses: update_mll_and_dmll!, update_cK!
using BenchmarkTools
using DataFrames
using JLD

const d = 10        # Input observation dimension
const nt = 3000     # Number of training points

kerns = Dict(
    "se" => SEIso(0.0,0.0),
    "mat12" => Mat12Iso(0.0,0.0),
    "rq" => RQIso(0.0,0.0,0.0),
    "se+rq" => SEIso(0.0,0.0) + RQIso(0.0,0.0,0.0),
    "se*rq" => SEIso(0.0,0.0) * RQIso(0.0,0.0,0.0),
    "se+se2+rq" => SEIso(0.0,0.0) + SEIso(0.5,0.5) + RQIso(0.0,0.0,0.0),
    "(se+se2)*rq" => (SEIso(0.0,0.0) + SEIso(0.5,0.5)) * RQIso(0.0,0.0,0.0),
    "mask(se, [1])" => Masked(SEIso(0.0,0.0), [1]),
    "mask(se, [1])+mask(rq, [2:10])" =>  Masked(SEIso(0.0,0.0), [1]) +  Masked(RQIso(0.0,0.0,0.0), collect(2:10)),
    "fix(se, σ)" => fix(SEIso(0.0,0.0), :lσ)
    )

function benchmark_kernel(group, kern)
    srand(1)
    X = randn(d, nt) # Training data
    Y = randn(nt)
    buf1=Array{Float64}(nt,nt)
    buf2=Array{Float64}(nt,nt)
    gp = GP(X, Y, MeanConst(0.0), kern, log(1.0))
    group["cK"] = @benchmarkable update_cK!($gp)
    group["mll_and_dmll"] = @benchmarkable update_mll_and_dmll!($gp, $buf1, $buf2)
end

SUITE = BenchmarkGroup()

for (label, k) in kerns
    SUITE[label] = BenchmarkGroup([label])
    benchmark_kernel(SUITE[label], k)
end
;

results = run(SUITE, verbose=false, seconds=1000, samples=10, evals=1)
;

knames = sort(collect(keys(kerns)))
times = [time(results[k]["mll_and_dmll"])/10^6 for k in knames]

df = DataFrame(kernel = knames, times=times)
print(df)

writetable("bench_results/GaussianProcesses_jl.csv", df, header=false)

srand(1)
X = randn(d, nt) # Training data
Y = randn(nt)
XY_df = DataFrame(X')
XY_df[:Y] = Y
XY_df
writetable("simdata.csv", XY_df, header=true)

buf1=Array{Float64}(nt,nt)
buf2=Array{Float64}(nt,nt)
gp = GPE(X, Y, MeanConst(0.0), kerns["se"], log(1.0))
update_mll_and_dmll!(gp, buf1, buf2)
gp.mll, gp.dmll # SE kernel

buf1=Array{Float64}(nt,nt)
buf2=Array{Float64}(nt,nt)
gp = GPE(X, Y, MeanConst(0.0), kerns["rq"], log(1.0))
Profile.clear()
update_mll_and_dmll!(gp, buf1, buf2)
@profile update_mll_and_dmll!(gp, buf1, buf2)
@profile update_mll_and_dmll!(gp, buf1, buf2)

Profile.print(mincount=10)


text:
using SparseArrays
using FastGaussQuadrature

include("SplineModel.jl")
include("MeshUtils.jl")

struct Boundary1D
    type::String #can be "Dirichlet, Neumann or Robin"
    x_val::Float64 #boundary point in the physical space
    u_val::Float64 #boundary point in the parameter space
    α::Number #parameter value for Robin boundary condition (αu+u')
    op_val::Number #boundary operator value
end

function Boundary1D(type::String, x_val::Float64, u_val::Float64, op_val::Number)
    return Boundary1D(type, x_val, u_val, 0., op_val)
end

struct Problem1D
    rhs_f::Function
    a0::Function
    a1::Function
    boundary_cond::Array{Boundary1D}
    domain::NURBS
end

struct GaussQuad
    nodes::Array{Float64,1}
    weights::Array{Float64,1}
end

"""
Generates the Gauss-Legendre points of order n
"""
function genGaussLegendre(n)
    nodes, weights = gausslegendre(n)
    gauss_rule = GaussQuad(nodes, weights)
    return gauss_rule
end

"""
Aseembles the stiffness matrix K_ij = ∫ a0(x)ϕ_i'(x)ϕ_j'(x) dΩ
"""
function assemble_stiff(mesh::Mesh, a0::Function, gauss_rule)
    B, dB = bernsteinBasis(gauss_rule.nodes, degP)
    domainLength = 0
    stiff = spzeros(mesh.numBasis, mesh.numBasis)
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 2]
        Jac_ref_par = (uMax-uMin)/2

        #compute the (B-)spline basis functions and derivatives with Bezier extraction
        N_mat = B * mesh.C[iElem]'
        dN_mat = dB * mesh.C[iElem]'/Jac_ref_par

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1, curNodes]
        wgts = mesh.weights[curNodes]
        localStiff = zeros(numNodes, numNodes)
        for iGauss = 1:length(gauss_rule.nodes)
            #compute the rational basis
            RR = N_mat[iGauss,:].* wgts
            dR = dN_mat[iGauss,:].* wgts
            w_sum = sum(RR)
            dw_xi = sum(dR)
            dR = dR/w_sum - RR*dw_xi/w_sum^2

            #compute the derivatives w.r.t to the physical space
            dxdxi = dR' * cpts
            dR = dxdxi\dR
            Jac_par_phys = det(dxdxi)

            RR /= w_sum
            phys_pt = RR'*cpts
            localStiff += Jac_par_phys * Jac_ref_par * a0(phys_pt) * (dR*dR') * gauss_rule.weights[iGauss]
            domainLength += Jac_par_phys * Jac_ref_par * gauss_rule.weights[iGauss]

        end
        stiff[curNodes, curNodes] += localStiff

    end
    @show domainLength
    return stiff
end

"""
Aseembles the mass matrix M_ij = ∫ a1(x)ϕ_i(x)ϕ_j(x) dΩ
"""
function assemble_mass(mesh::Mesh, a1::Function, gauss_rule)
    B, dB = bernsteinBasis(gauss_rule.nodes, degP)
    domainLength = 0
    mass = spzeros(Complex{Float64}, mesh.numBasis, mesh.numBasis)
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 2]
        Jac_ref_par = (uMax-uMin)/2

        #compute the (B-)spline basis functions and derivatives with Bezier extraction
        N_mat = B * mesh.C[iElem]'
        dN_mat = dB * mesh.C[iElem]'/Jac_ref_par

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1, curNodes]
        wgts = mesh.weights[curNodes]
        localMass = zeros(numNodes, numNodes)
        for iGauss = 1:length(gauss_rule.nodes)
            #compute the rational basis
            RR = N_mat[iGauss,:].* wgts
            dR = dN_mat[iGauss,:].* wgts
            w_sum = sum(RR)
            dw_xi = sum(dR)
            dR = dR/w_sum - RR*dw_xi/w_sum^2

            #compute the Jacobian of the transformation from parameter to physical space
            dxdxi = dR' * cpts
            Jac_par_phys = det(dxdxi)

            RR /= w_sum
            phys_pt = RR'*cpts
            localMass += a1(phys_pt) * (RR*RR') * Jac_par_phys * Jac_ref_par * gauss_rule.weights[iGauss]
            domainLength += Jac_par_phys * Jac_ref_par * gauss_rule.weights[iGauss]

        end
        #@show localMass
        #readline(stdin)
        mass[curNodes, curNodes] += localMass
    end
    @show domainLength
    return mass
end

"""
Assembles the RHS vector corresponding to the body force
RHS[i]=∫_Ω ϕ_i(x)*f(x) dΩ
"""
function assemble_rhs(mesh::Mesh, f::Function, gauss_rule)
    B, dB = bernsteinBasis(gauss_rule.nodes, degP)
    domainLength = 0
    rhs = zeros(Complex{Float64}, mesh.numBasis)
    for iElem = 1:mesh.numElem
        uMin = mesh.elemVertex[iElem, 1]
        uMax = mesh.elemVertex[iElem, 2]
        Jac_ref_par = (uMax-uMin)/2
        #@show Jac_ref_par
        #compute the (B-)spline basis functions and derivatives with Bezier extraction
        N_mat = B * mesh.C[iElem]'
        dN_mat = dB * mesh.C[iElem]'/Jac_ref_par

        #compute the rational spline basis
        curNodes = mesh.elemNode[iElem]
        numNodes = length(curNodes)
        cpts = mesh.controlPoints[1, curNodes]
        wgts = mesh.weights[curNodes]
        localRhs = zeros(numNodes)
        for iGauss = 1:length(gauss_rule.nodes)
            #compute the rational basis
            RR = N_mat[iGauss,:].* wgts
            dR = dN_mat[iGauss,:].* wgts
            w_sum = sum(RR)
            dw_xi = sum(dR)
            dR = dR/w_sum - RR*dw_xi/w_sum^2

            #compute the derivatives w.r.t to the physical space
            dxdxi = dR' * cpts
            Jac_par_phys = det(dxdxi)
            RR /= w_sum

            phys_pt = RR'*cpts
            localRhs += Jac_par_phys * Jac_ref_par * f(phys_pt) * RR * gauss_rule.weights[iGauss]
            domainLength += Jac_par_phys * Jac_ref_par * gauss_rule.weights[iGauss]
        end
        rhs[curNodes] += localRhs
    end
    @show domainLength
    return rhs
end

function uniquetol(arr, tol=1e-10)
    sort!(arr)
    delete_index = Array{Int64,1}(undef, 0)
    for i=1:length(arr)-1
        if abs(arr[i+1]-arr[i])<tol
            push!(delete_index, i+1)
        end
    end
    deleteat!(arr, delete_index)
end

"""
Helper function to compute the setdifference within a tolerance
"""
function setdifftol(arrA, arrB, tol=1e-10)
    delete_index = Array{Int64,1}(undef, 0)
    for i=1:length(arrA)
        for j=1:length(arrB)
            if abs(arrA[i]-arrB[j])<tol
                push!(delete_index, i)
                break
            end
        end
    end
    deleteat!(arrA, delete_index)
end

"""
Apply the boundary conditions to a linear system for a 1D
nurbs geometry
"""
function applyBCnurbs(bound_cond, stiff, mass, rhs, nrb)
    #collect the dofs and values corresponding to the boundary
    lhs = stiff+mass
    bcdof = Array{Int64,1}(undef, 0)
    bcval = Array{Float64,1}(undef, 0)
    for i=1:length(bound_cond)

        if bound_cond[i].type=="Dirichlet"
            bcdof = vcat(bcdof, findall(nrb.coefs[1,:].==bound_cond[i].x_val))
            bcval = vcat(bcval, bound_cond[i].op_val)
        elseif bound_cond[i].type=="Neumann"
            bcdof_neu = findall(nrb.coefs[1,:].==bound_cond[i].x_val)
            rhs[bcdof_neu] .+= bound_cond[i].op_val
        elseif bound_cond[i].type=="Robin"
            #Assume that the basis functions have value 1 on the Robin boundary
            bcdof_rob = findall(nrb.coefs[1,:].==bound_cond[i].x_val)
            lhs[bcdof_rob, bcdof_rob] .+= bound_cond[i].α
            rhs[bcdof_rob] .+= bound_cond[i].op_val
        else
            error("Boundary condition type not implemented")
        end
    end
    rhs = rhs - lhs[:,bcdof]*bcval
    rhs[bcdof] = bcval
    lhs[bcdof, :] .= 0.
    lhs[:, bcdof] .= 0.
    lhs[bcdof, bcdof] = sparse(I, length(bcdof), length(bcdof))
    return lhs, rhs
end


text:
using LinearAlgebra

@testset "System" begin
    Impl.@polyvar x y z

    f = Polynomial{Float64}(x^2+3y+z*x+y^2+12*x^7*z^5)
    g = Polynomial{Float64}(4*y*x^2-3z+z*x^4+y+y^2-y*x)

    ∇f = [differentiate(f, i) for i=1:3]
    ∇g = [differentiate(g, i) for i=1:3]

    w = rand(3)
    u = zeros(3)

    @test_throws MethodError JacobianConfig(f, w)

    u = zeros(2)
    F = System([f, g])
    cfg = config(F, w)
    @test [f(w), g(w)] ≈ evaluate(F, w, cfg)
    @test [f(w), g(w)] ≈ evaluate!(u, F, w, cfg)
    @test [f(w), g(w)] ≈ u

    @test_throws BoundsError evaluate(F, rand(2), cfg)

    u = zeros(2)
    cfg = JacobianConfig(F, w)
    jacobian(F, w, cfg)
    @test [f(w), g(w)] ≈ evaluate(F, w, cfg, true)
    @test [f(w), g(w)] ≈ evaluate!(u, F, w, cfg, true)
    @test [f(w), g(w)] ≈ u

    cfg2 = deepcopy(cfg)
    @test cfg2 !== cfg

    U = zeros(2, 3)

    DF = [∇f[1](w) ∇f[2](w) ∇f[3](w); ∇g[1](w) ∇g[2](w) ∇g[3](w)]
    evaluate(F, w, cfg)
    @test DF ≈ jacobian(F, w, cfg, true)
    @test DF ≈ jacobian!(U, F, w, cfg, true)
    @test DF ≈ U
    @test_throws BoundsError jacobian(F, rand(2), cfg)
    @test_throws BoundsError jacobian!(U[1:end-1,:],F, w, cfg)
    @test_throws BoundsError jacobian!(U[:,1:end-1],F, w, cfg)

    r = JacobianDiffResult(cfg)
    jacobian!(r, F, w, cfg)
    @test DF ≈ jacobian(r)
    @test [f(w), g(w)] ≈ value(r)

    u, U = evaluate_and_jacobian(F, w, cfg)
    @test U ≈ jacobian(r)
    @test u ≈ value(r)

    v = zeros(3)
    V = zeros(3, 3)

    r = JacobianDiffResult((@view v[1:2]), (@view V[1:2, :]))
    jacobian!(r, F, w, cfg)
    @test DF ≈ jacobian(r)
    @test [f(w), g(w)] ≈ value(r)
end


text:
export @dropdims

using MacroTools

"""
    @dropdims sum(A; dims=1)

Macro which wraps such reductions in `dropdims(...; dims=1)`.
Allows `sum(A; dims=1) do x stuff end`,
and works on whole blocks of code like `@views`.
Does not handle other keywords, like `reduce(...; dims=..., init=...)`.
"""
macro dropdims(ex)
    _dropdims(ex)
end

function _dropdims(ex)
    out = MacroTools.postwalk(ex) do x
        if @capture(x, red_(args__, dims=d_)) || @capture(x, red_(args__; dims=d_))
            :( dropdims($x; dims=$d) )
        elseif @capture(x, dropdims(red_(args__, dims=d1_); dims=d2_) do z_ body_ end) ||
               @capture(x, dropdims(red_(args__; dims=d1_); dims=d2_) do z_ body_ end)
            :( dropdims($red($z -> $body, $(args...); dims=$d1); dims=$d2) )
        else
            x
        end
    end
    esc(out)
end


text:
# This file was generated, do not modify it. # hide
#hideall
println("~~~") # start HTML
println("
")
for f in readdir(joinpath("src", "pages"))
    fn = splitext(splitdir(f)[2])[1]
    println("
")
    println("", fn, "")
    println("
")
end
println("
")
println("~~~") # end HTML

max_stars_repo_name: abhishalya/julialang-judoc

max_issues_repo_name: abhishalya/julialang-judoc

max_forks_repo_name: abhishalya/julialang-judoc


text:
using DataStructures
using XMLDict

function FunctionDescriptionXML()
    unsafe_string(ccall(("TA_FunctionDescriptionXML", libta_lib), Cstring, ()))
end

function _list2ordereddict(lst_ta_func)
    d = OrderedDict{Symbol,Any}()
    for func_info in lst_ta_func
        funcname = func_info["Abbreviation"]
        delete!(func_info, "Abbreviation")
        l = []  # length of RequiredInputArguments, OptionalInputArguments, OutputArguments
        for key in ["RequiredInputArgument", "OptionalInputArgument", "OutputArgument"]
            if haskey(func_info, key * "s")
                if typeof(func_info[key * "s"][key]) <: AbstractDict  # if it's a dict like
                    func_info[key * "s"] = [func_info[key * "s"][key]]  # list of only ONE element
                else
                    func_info[key * "s"] = func_info[key * "s"][key]
                end
            else
                func_info[key * "s"] = []
            end

            for arg in func_info[key * "s"]
                if !haskey(d_typ_to_c, arg["Type"])
                    error("$(arg["Type"]) is not a supported type")
                end
            end

            push!(l, length(func_info[key * "s"]))

        end
        #func_info["Length"] = l
        d[Symbol(funcname)] = func_info
    end
    d
end

#=
function get_dict_of_ta_func(filename::String)
    if !isfile(filename)
        error("$filename doesn't exist")
    end
    lst_ta_func = JSON.parsefile(filename)["FinancialFunctions"]["FinancialFunction"]
    _list2ordereddict(lst_ta_func)
end
=#

function get_dict_of_ta_func()
    s_xml = FunctionDescriptionXML()
    d_xml = xml_dict(s_xml)["FinancialFunctions"]["FinancialFunction"]
    _list2ordereddict(d_xml)
end

function get_ta_func_constants()
    #indicators = get_dict_of_ta_func(joinpath(basepath(), "generated", "ta_func_api.json"))
    d_indicators = get_dict_of_ta_func()
    indicators = Symbol[func for func in keys(d_indicators)]

    d_groups = Dict{String, Array{Symbol,1}}()
    for (func, val) in d_indicators
        key = val["GroupId"]
        if !haskey(d_groups, key)
            d_groups[key] = Symbol[func]
        else
            push!(d_groups[key], func)
        end
    end
    d_indicators, indicators, d_groups
end


text:
using SmartShockFinder

struct Traindata_Settings
    polydegrees::Vector{Int64}
    meshes::Vector{Any}
    length_data::Integer
end

function Traindata_Settings(polydegrees, meshes)
    # Calculate an upper bound for the number of entries
    # length for 1 mesh = 2 * #elements in mesh * #polydegrees * #functions
    # the # of functions is per default 26
    length_data = 0
    for i in eachindex(meshes)
        length_data += length(leaf_cells(meshes[i].tree)) * length(polydegrees) * 26
    end
    length_data = max(length_data,25000)
    Traindata_Settings(polydegrees, meshes, 2 * length_data)
end

# 1D implementation to generate dataset
include("traindata1d.jl")
include("validdata1d.jl")

# 2D implementation
include("traindata2d.jl")
include("validdata2d.jl")


text:
function GenericTimeStampedCache(arg0::jint, arg1::jint, arg2::jdouble, arg3::jdouble, arg4::TimeStampedGenerator)
    return GenericTimeStampedCache((jint, jint, jdouble, jdouble, TimeStampedGenerator), arg0, arg1, arg2, arg3, arg4)
end

function get_earliest(obj::GenericTimeStampedCache)
    return jcall(obj, "getEarliest", TimeStamped, ())
end

function get_entries(obj::GenericTimeStampedCache)
    return jcall(obj, "getEntries", jint, ())
end

function get_generate_calls(obj::GenericTimeStampedCache)
    return jcall(obj, "getGenerateCalls", jint, ())
end

function get_generator(obj::GenericTimeStampedCache)
    return jcall(obj, "getGenerator", TimeStampedGenerator, ())
end

function get_get_neighbors_calls(obj::GenericTimeStampedCache)
    return jcall(obj, "getGetNeighborsCalls", jint, ())
end

function get_latest(obj::GenericTimeStampedCache)
    return jcall(obj, "getLatest", TimeStamped, ())
end

function get_max_slots(obj::GenericTimeStampedCache)
    return jcall(obj, "getMaxSlots", jint, ())
end

function get_max_span(obj::GenericTimeStampedCache)
    return jcall(obj, "getMaxSpan", jdouble, ())
end

function get_neighbors(obj::GenericTimeStampedCache, arg0::AbsoluteDate)
    return jcall(obj, "getNeighbors", Stream, (AbsoluteDate,), arg0)
end

function get_neighbors_size(obj::GenericTimeStampedCache)
    return jcall(obj, "getNeighborsSize", jint, ())
end

function get_new_slot_quantum_gap(obj::GenericTimeStampedCache)
    return jcall(obj, "getNewSlotQuantumGap", jdouble, ())
end

function get_slots(obj::GenericTimeStampedCache)
    return jcall(obj, "getSlots", jint, ())
end

function get_slots_evictions(obj::GenericTimeStampedCache)
    return jcall(obj, "getSlotsEvictions", jint, ())
end


text:
dv1 = DataVector[9,1,8,NA,3,3,7,NA]
dv2 = 1.0 * dv1
dv3 = DataArray([1:8])
pdv1 = PooledDataArray(dv1)
idv1 = IndexedVector(dv1)

@assert sortperm(dv1) == sortperm(dv2)
@assert sortperm(dv1) == sortperm(pdv1)
@assert sortperm(dv1) == sortperm(idv1)
@assert isequal(sort(dv1), DataArray(sort(dv1)))
@assert isequal(sort(dv1), DataArray(sort(pdv1)))
@assert isequal(sort(dv1), DataArray(sort(idv1)))

d = @DataFrame(dv1 => dv1, dv2 => dv2, dv3 => dv3, pdv1 => pdv1, idv1 => idv1)

@assert sortperm(d) == sortperm(dv1)
@assert sortperm(d[["dv3","dv1"]]) == sortperm(dv3)
@assert sortby(d, "dv1")["dv3"] == sortperm(dv1)
@assert sortby(d, "dv2")["dv3"] == sortperm(dv1)
@assert sortby(d, "pdv1")["dv3"] == sortperm(dv1)
@assert sortby(d, "idv1")["dv3"] == sortperm(dv1)
@assert sortby(d, ["dv1","pdv1"])["dv3"] == sortperm(dv1)
@assert sortby(d, ["dv1","dv3"])["dv3"] == sortperm(dv1)


text:
using MethodOfLines, ModelingToolkit, DomainSets
using Combinatorics: permutations

@testset "Test 00: recognize relevant variable for equations, time defined" begin
    @parameters x, t
    @variables u(..), v(..), w(..)

    Dx(d) = Differential(x)^d

    Dt = Differential(t)

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 20.0

    dx = 1.0

    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]

    pde = [Dt(u(t,x)) ~ +Dx(1)(u(t,x))+w(t,x),
           Dt(v(t,x)) ~ -Dx(2)(v(t,x))+w(t,x),
           Dt(w(t,x)) ~ -Dx(3)(w(t,x))+u(t,x)+v(t,x)]
    bcs = [u(0,x) ~ 0, u(t,0) ~ 0, u(t,Float64(π)) ~ 0]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x), v(t,x), w(t,x)])

    # Test centered order
    disc = MOLFiniteDifference([x=>dx], t)

    depvar_ops = map(x->operation(x.val),pdesys.depvars)
    depvars = MethodOfLines.get_all_depvars(pdesys, depvar_ops)
    indvars = MethodOfLines.remove(collect(reduce(union,filter(xs->!isequal(xs, [t]), map(arguments, depvars)))),t)

    s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, disc)

    m = MethodOfLines.buildmatrix(pde, s)
    if VERSION >= v"1.7"
        @test m == [2 0 1; 0 3 1; 1 1 4] # Test the matrix is the identity matrix
    else
        @test m == [1 2 0; 1 0 3; 4 1 1]
    end

end

@testset "Test 00a: recognize relevant variable for equations, time undefined, multiple choices" begin
    @parameters x, t
    @variables u(..), v(..), w(..)

    Dx(d) = Differential(x)^d

    Dt = Differential(t)

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 20.0

    dx = 1.0

    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]

    pde = [Dt(u(t,x)) ~ +Dx(1)(u(t,x))+Dx(1)(w(t,x)),
           Dt(v(t,x)) ~ -Dx(2)(v(t,x))+Dx(2)(w(t,x)),
           Dt(w(t,x)) ~ -Dx(3)(w(t,x))+Dx(3)(u(t,x))+Dx(3)(v(t,x))]
    bcs = [u(0,x) ~ 0, u(t,0) ~ 0, u(t,Float64(π)) ~ 0]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x), v(t,x), w(t,x)])

    # Test centered order
    disc = MOLFiniteDifference([x=>dx, t=>dx])

    depvar_ops = map(x->operation(x.val),pdesys.depvars)
    depvars = MethodOfLines.get_all_depvars(pdesys, depvar_ops)
    indvars = collect(reduce(union,filter(xs->!isequal(xs, [t]), map(arguments, depvars))))

     s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, disc)

    m = MethodOfLines.buildmatrix(pde, s)
    if VERSION >= v"1.7"
        @test m == [2 0 2; 0 3 3; 4 4 4] # Test the matrix is the identity matrix
    else
        @test m == [2 2 0; 3 0 3; 4 4 4]
    end
end
#
@testset "Test 00b: recognize relevant variable for equations, time undefined, mixed derivatives, multiple choices" begin
    @parameters x, t
    @variables u(..), v(..), w(..)

    Dx(d) = Differential(x)^d

    Dt = Differential(t)

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 20.0

    dx = 1.0

    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]

    pde = [Dt(u(t,x)) ~ +Dx(1)(u(t,x))+w(t,x),
           Dt(v(t,x)) ~ -Dx(2)(v(t,x))+Dx(2)(w(t,x)),
           Dt(w(t,x)) ~ -Dx(3)(w(t,x))+Dx(3)(u(t,x)+v(t,x))]
    bcs = [u(0,x) ~ 0, u(t,0) ~ 0, u(t,Float64(π)) ~ 0]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x), v(t,x), w(t,x)])

    # Test centered order
    disc = MOLFiniteDifference([x=>dx, t=>1.0])

    depvar_ops = map(x->operation(x.val),pdesys.depvars)
    depvars = MethodOfLines.get_all_depvars(pdesys, depvar_ops)
    indvars = collect(reduce(union,filter(xs->!isequal(xs, [t]), map(arguments, depvars))))

    s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, disc)

    m = MethodOfLines.buildmatrix(pde, s)
    if VERSION >= v"1.7"
        @test m == [2 0 1; 0 3 3; 5 5 4] # Test the matrix is the identity matrix
    else
        @test m == [1 2 0; 3 0 3; 4 5 5]
    end

end

@testset "Test 00c: recognize relevant variable for equations, time undefined, ranking important" begin
    @parameters x, t
    @variables u(..), v(..), w(..)

    Dx = Differential(x)

    Dt = Differential(t)

    t_min= 0.
    t_max = 2.0
    x_min = 0.
    x_max = 20.0

    dx = 1.0

    domains = [t ∈ Interval(t_min, t_max), x ∈ Interval(x_min, x_max)]

    pde  = [Dt(u(t,x)) ~ Dx(v(t,x)),
            u(t,x) ~ w(t,x)*w(t,x),
            v(t,x) ~ w(t,x)*Dx(w(t,x))]

    bcs = [u(0,x) ~ 0, u(t,0) ~ 0, u(t,Float64(π)) ~ 0]

    @named pdesys = PDESystem(pde,bcs,domains,[t,x],[u(t,x), v(t,x), w(t,x)])

    # Test centered order
    disc = MOLFiniteDifference([x=>dx, t=>1.0])

    depvar_ops = map(x->operation(x.val),pdesys.depvars)
    depvars = MethodOfLines.get_all_depvars(pdesys, depvar_ops)
    indvars = collect(reduce(union,filter(xs->!isequal(xs, [t]), map(arguments, depvars))))

    s = MethodOfLines.DiscreteSpace(domains, depvars, indvars, disc)

    m = MethodOfLines.buildmatrix(pde, s)
    if VERSION >= v"1.7"
        @test m == [2 2 0; 1 0 1; 0 1 2]
    else
        @test m == [0 2 2; 1 1 0; 2 0 1]
    end

end

@testset "Test 01a: Build variable mapping - one right choice simple" begin
    m = hcat([0, 1, 0],
             [0, 0, 1],
             [1, 0, 0])
    pdes = ["a", "b", "c"]
    vars = ["x", "y", "z"]

    @test Dict(MethodOfLines.build_variable_mapping(m, vars, pdes)) == Dict([
        "a" => "z",
        "b" => "x",
        "c" => "y"
    ])
end

@testset "Test 01b: Build variable mapping - one right choice complex" begin
    m = hcat([0, 1, 1],
             [1, 0, 1],
             [0, 1, 0])
    pdes = ["a", "b", "c"]
    vars = ["x", "y", "z"]

    @test Dict(MethodOfLines.build_variable_mapping(m, vars, pdes)) == Dict([
        "a" => "y",
        "b" => "z",
        "c" => "x"
    ])
end

@testset "Test 01c: Build variable mapping - two right choices" begin
    m = hcat([0, 1, 1],
             [1, 0, 1],
             [1, 1, 0])
    pdes = ["a", "b", "c"]
    vars = ["x", "y", "z"]
    out = Dict(MethodOfLines.build_variable_mapping(m, vars, pdes))
    @test out == Dict([
        "a" => "y",
        "b" => "z",
        "c" => "x"
    ]) || out == Dict([
        "a" => "z",
        "b" => "x",
        "c" => "y"
    ])
end

@testset "Test 01d: Build variable mapping - any choice correct" begin
    m = hcat([1, 1, 1],
             [1, 1, 1],
             [1, 1, 1])
    pdes = ["a", "b", "c"]
    vars = ["x", "y", "z"]

    perms = permutations(vars)
    dict = Dict(MethodOfLines.build_variable_mapping(m, vars, pdes))
    @test any(map(perm -> dict == Dict(["a" => perm[1],
                                        "b" => perm[2],
                                        "c" => perm[3]]), perms))
end

@testset "Test 01e: Build variable mapping - no correct mapping" begin
    m = hcat([1, 1, 0],
             [1, 1, 0],
             [1, 1, 0])

    pdes = ["a", "b", "c"]
    vars = ["x", "y", "z"]
    try
        MethodOfLines.build_variable_mapping(m, vars, pdes)
        @test false
    catch e
        @test e isa AssertionError
    end
end


text:
using Base.Test

reload("ForwardBackwardOptim")
m = ForwardBackwardOptim

tests = [
    "optims"
]

for t in tests
    tfile = string(t, ".jl")
    println("  * $tfile ...")
    include(tfile)
end

println("Finished testing.")


text:
#=
# The final symbolic form will be an Expr expression tree.
# The leaf nodes are all SymEntitys as described below.
# The parents will be math operators.
# All used entities will also be kept in a list.
#
# A SymEntity has
# - name, the name of this variable, coefficient, etc. (or number for constants)
# - index, the component of this, such as a vector component (-1 for constants)
# - derivs
# - flags
=#
export SymExpression, SymEntity, build_symexpressions

struct SymExpression
    tree                # The expression tree
    entities::Array     # The SymEntitys corresponding to the leaf nodes of the tree
end

struct SymEntity
    name::Union{Float64, String}    # The string or number
    index::Union{Int64, Array}     # vector component, 1 for scalars, -1 for numbers, "INDEXEDBY..." for indexed vars
    derivs::Array                   # any derivatives applied
    flags::Array                    # any non-derivative modifiers or flags
end

Base.copy(x::SymEntity) = SymEntity(x.name, x.index, copy(x.derivs), copy(x.flags));

# For printing, write the entity in a more readable format
Base.show(io::IO, x::SymEntity) = print(io, symentity_string(x));
function symentity_string(a::SymEntity)
    result = "";
    for i=1:length(a.derivs)
        result *= "D"*string(a.derivs[i]) * "(";
    end

    for i=1:length(a.flags)
        result *= a.flags[i];
    end

    result *= a.name;
    if typeof(a.index) == Int && a.index > 0
        result *= "_"*string(a.index);
    elseif typeof(a.index) <: Array && length(a.index) > 0
        if length(a.index) > 1
            result *= "_{"*a.index[1];
            for i=2:length(a.index)
                result *= ","*a.index[i];
            end
            result *= "}";
        elseif length(a.index) == 1
            result *= "_"*a.index[1];
        end
    end

    for i=1:length(a.derivs)
        result *= ")";
    end

    return result;
end

# Make a Latex string for the expression
# Recursively traverses the tree turning each piece into latex.
function symexpression_to_latex(ex)
    result = "";

    if typeof(ex) <: Array
        if length(ex) < 1
            return "";
        end

        if typeof(ex[1]) <: Array
            result = "\\left[";
            for i=1:length(ex)
                if i > 1
                    result *= ", ";
                end
                result *= symexpression_to_latex(ex[i]);
            end
            result *= "\\right]";

        else # The innermost array is an array of added terms
            result = symexpression_to_latex(ex[1]);
            for i=2:length(ex)
                result *= " + " * symexpression_to_latex(ex[i]);
            end
        end

    elseif typeof(ex) == SymExpression
        result = symexpression_to_latex(ex.tree);

    elseif typeof(ex) == Expr
        if ex.head === :call
            if ex.args[1] in [:+, :.+, :-, :.-, :*, :.*] && length(ex.args) > 2
                if typeof(ex.args[2]) == Expr
                    result *= "\\left(";
                end
                result *= symexpression_to_latex(ex.args[2]);
                if typeof(ex.args[2]) == Expr
                    result *= "\\right)";
                end
                for i=3:length(ex.args)
                    result *= " "*string(ex.args[1])[end] * " ";
                    if typeof(ex.args[i]) == Expr
                        result *= "\\left(";
                    end
                    result *= symexpression_to_latex(ex.args[i]);
                    if typeof(ex.args[i]) == Expr
                        result *= "\\right)";
                    end
                end

            elseif ex.args[1] in [:-, :.-] && length(ex.args) == 2 # negative
                result *= "-";
                if typeof(ex.args[2]) == Expr
                    result *= "\\left(";
                end
                result *= symexpression_to_latex(ex.args[2]);
                if typeof(ex.args[2]) == Expr
                    result *= "\\right)";
                end

            elseif ex.args[1] in [:/, :./]
                top = symexpression_to_latex(ex.args[2]);
                if length(ex.args) > 3
                    # multiply them first
                    newex = :(a*b);
                    newex.args = [:*];
                    append!(newex.args, ex.args[3:end]);
                    bottom = symexpression_to_latex(newex);
                else
                    bottom = symexpression_to_latex(ex.args[3]);
                end

                result = "\\frac{" * top * "}{" * bottom * "}";

            elseif ex.args[1] in [:^, :.^]
                bottom = symexpression_to_latex(ex.args[2]);
                power = symexpression_to_latex(ex.args[3]);
                result = bottom * "^{" * power * "}";

            elseif ex.args[1] in [:sqrt]
                result = "\\sqrt{" * symexpression_to_latex(ex.args[2]) * "}";

            elseif ex.args[1] in [:sin, :cos, :tan]
                result = "\\"*string(ex.args[1]) * "\\left(" * symexpression_to_latex(ex.args[2]) * "\\right)";

            elseif ex.args[1] === :abs
                result = "\\left|" * symexpression_to_latex(ex.args[2]) * "\\right|";

            elseif ex.args[1] === :conditional_function
                result =  "conditional\\left(" * symexpression_to_latex(ex.args[2]);
                for i=3:length(ex.args)
                    result *= ", " * symexpression_to_latex(ex.args[i]);
                end
                result *= "\\right)";

            else
                result = string(ex.args[1]) * "\\left(" * symexpression_to_latex(ex.args[2]);
                for i=3:length(ex.args)
                    result *= ", " * symexpression_to_latex(ex.args[i]);
                end
                result *= "\\right)";
            end

        elseif ex.head === :.
            # This is a broadcast symbol. args[1] is the operator, args[2] is a tuple for the operator
            # eg. change :(sin.(a)) to :(sin(a))
            newex = Expr(:call, ex.args[1]);
            append!(newex.args, ex.args[2].args);
            result = symexpression_to_latex(newex);

        else
            # There are some other possible symbols, but they could be ignored
            result = symexpression_to_latex(ex.args[1]);
        end

    elseif typeof(ex) == SymEntity
        # Apply derivatives
        dimchars = ["x","y","z"];
        name = ex.name;
        if name == "FACENORMAL1"
            name = "normal";
        elseif name == "FACENORMAL2"
            name = "-normal";
        end
        if typeof(ex.index) == Int && ex.index > 0
            result = name * "_{"*string(ex.index);
            if "DGSIDE1" in ex.flags || "CELL1" in ex.flags
                result *= "+";
            elseif "DGSIDE2" in ex.flags || "CELL2" in ex.flags
                result *= "-";
            end
            if length(ex.derivs) > 0
                result *= ","
                for i=1:length(ex.derivs)
                    result *= dimchars[ex.derivs[i]];
                end
            end
            result *= "}";
        else
            result = name;
        end

    elseif typeof(ex) <: Number
        result = string(ex);

    else
        result = string(ex);
    end

    return result;
end

# Builds all of them one at a time.
function build_symexpressions(var, lhsvol, rhsvol, lhssurf=nothing, rhssurf=nothing; remove_zeros=false)
    lhsvolsym = build_symexpression(var, lhsvol, "L", "V", remove_zero_in_array=remove_zeros);
    rhsvolsym = build_symexpression(var, rhsvol, "R", "V", remove_zero_in_array=remove_zeros);
    lhssurfsym = build_symexpression(var, lhssurf, "L", "S", remove_zero_in_array=remove_zeros);
    rhssurfsym = build_symexpression(var, rhssurf, "R", "S", remove_zero_in_array=remove_zeros);

    if lhssurf === nothing && rhssurf === nothing
        return (lhsvolsym, rhsvolsym);
    else
        return (lhsvolsym, rhsvolsym, lhssurfsym, rhssurfsym);
    end
end

# Builds one symexpression from a SymEngine object or an Expr
function build_symexpression(var, ex, lorr, vors; remove_zero_in_array=false)
    if ex === nothing
        return nothing
    end

    if typeof(ex) <: Array
        exarray = [];
        for i=1:length(ex)
            if !(remove_zero_in_array && ex[i] == 0)
                push!(exarray, build_symexpression(var, ex[i], lorr, vors, remove_zero_in_array=remove_zero_in_array));
            end
        end
        return exarray;

    elseif typeof(ex) <: SymEngine.Basic
        # ex is a SymEngine.Basic object
        # First parse it into an Expr
        jex = Meta.parse(string(ex));

        # Then traverse it to find SymTerm leaves and build a SymExpression
        (newex, sents) = extract_symentities(var, jex);

        # Replace some special operators with their real expressions
        newex = replace_special_ops(newex);

        # If using FEM, reorder factors for easier generation
        if config.solver_type == CG || config.solver_type == DG
            newex = order_expression_for_fem(var, test_functions, newex);
        end

        symex = SymExpression(newex, sents);

        return symex;

    elseif typeof(ex) == Expr || typeof(ex) == Symbol || typeof(ex) <: Number
        # Traverse it to find SymTerm leaves and build a SymExpression
        (newex, sents) = extract_symentities(var, ex);

        # Replace some special operators with their real expressions
        newex = replace_special_ops(newex);

        # If using FEM, reorder factors for easier generation
        if config.solver_type == CG || config.solver_type == DG
            newex = order_expression_for_fem(var, test_functions, newex);
        end

        symex = SymExpression(newex, sents);

        return symex;

    end
end

# Traverse the tree and replace appropriate symbols with SymEntitys.
# Also, keep a list of symentities used.
function extract_symentities(var, ex)
    sents = [];
    if typeof(ex) == Expr
        # Recursively traverse the tree
        newex = copy(ex);
        if ex.head === :call # an operation like +(a,b)
            for i=2:length(ex.args)
                (subex, subsents) = extract_symentities(var, ex.args[i]);
                newex.args[i] = subex;
                append!(sents, subsents);
            end
        elseif ex.head === :if # a conditional like if a<b c else d end
            # This has three expr: one for a<b and a block for c and a block for d
            for i=1:length(ex.args)
                (subex, subsents) = extract_symentities(var, ex.args[i]);
                newex.args[i] = subex;
                append!(sents, subsents);
            end

        else # There are a few other things it could be. Just process the args.
            for i=1:length(ex.args)
                (subex, subsents) = extract_symentities(var, ex.args[i]);
                newex.args[i] = subex;
                append!(sents, subsents);
            end
        end

    elseif typeof(ex) == Symbol
        sen = build_symentity(var, ex);
        push!(sents, sen);
        newex = sen;

    else # probably a number
        newex = ex;
    end

    return (newex, sents);
end

# Builds a SymEntity struct from a symbolic name such as:
# DGSIDE1_D2__u_1 -> name=u, index=1, derivs=D2, flags=DGSIDE1
function build_symentity(var, ex)
    if typeof(ex) <: Number
        return SymEntity(ex, 0, [], []);
    elseif typeof(ex) == Symbol
        # Extract the index, derivs, flags
        (ind, symb, mods) = extract_entity_parts(ex);
        # Separate the derivative mods from the other flags
        (derivs, others) = get_derivative_mods(mods);

        return SymEntity(symb, ind, derivs, others);

    else
        printerr("unexpected type in build_symentity: "*string(typeof(ex))*" : "*string(ex));
    end
end

# Extract the indexing symbols from the index string
function extract_index_symbols_from_index_string(ind_str)
    indices = [];
    b = 10;
    e = 10;
    for i=10:length(ind_str)
        if ind_str[i-1:i] == "BY"
            e = i-2;
            push!(indices, ind_str[b:e]);
            b = i+1;
        elseif i == length(ind_str)
            push!(indices, ind_str[b:end]);
        end
    end

    return indices;
end

# Extract meaning from the symbolic object name
# The format of the input symbol should look like this
#   MOD1_MOD2_..._var_n
# There could be any number of mods, _var_n is the symvar symbol (n is the vector index, or 1 for scalar, or INDEXEDBY... for array type)
# Returns (n, var, [MOD1,MOD2,...]) all as strings
function extract_entity_parts(ex)
    str = string(ex);
    index = 0;
    var = "";
    mods = [];
    l = lastindex(str);
    e = l; # end of variable name
    b = l; # beginning of variable name

    # dt is a special symbol that will be passed as a number value in the generated function.
    if str == "dt"
        return(-1, str, []);
    end

    # Extract the index
    for i=l:-1:1
        if str[i] == '_' && i < l
            index_str = str[i+1:l];
            # This could be an integer or INDEXEDBY...
            if occursin("INDEXED", index_str)
                index = extract_index_symbols_from_index_string(index_str);
            else
                try
                    index = parse(Int, index_str);
                catch
                    return (-1,str,[]); # an unexpected index, perhaps it is a special symbol
                end
            end

            e = i-1; # end of variable name
            break;
        end
    end

    # Extract the name
    for i=e:-1:1
        if str[i] == '_'
            var = str[i+1:e];
            b = i-1;
            break;
        end
    end

    # extract the modifiers like D1_ etc. separated by underscores
    if b>1
        e = b-1;
        for i=e:-1:1
            if str[i] == '_'
                b = i+1;
                push!(mods, str[b:e]);
                e = b-2;

            elseif i == 1
                push!(mods, str[1:e]);
            end
        end
    end

    return (index, var, mods);
end

# Separate the derivative mods(Dn_) from the other mods and return an array
# of derivative indices and the remaining mods.
function get_derivative_mods(mods)
    derivs = [];
    others = [];

    for i=1:length(mods)
        if length(mods[i]) == 2 && mods[i][1] == 'D'
            try
                index = parse(Int, mods[i][2])
                push!(derivs, index);
            catch
                printerr("Unexpected modifier: "*mods[i]*" see get_derivative_mods() in symexpression.jl")
                push!(others, mods[i]);
            end
        else
            push!(others, mods[i]);
        end
    end

    return (derivs, others);
end

# changes special operators into the intended expressions
function replace_special_ops(ex)
    if typeof(ex) == Expr
        newex = copy(ex);
        for i=1:length(ex.args)
            newex.args[i] = replace_special_ops(newex.args[i]);
        end

    elseif typeof(ex) == Symbol
        # check against these special ones that need to be replaced
        # This will eventually contain more complex things.
        specials = [:symbolmin, :symbolmax, :isgreaterthan, :islessthan, :conditional];
        replacements = Dict([(:symbolmin, :min), (:symbolmax, :max), (:isgreaterthan, :>), (:islessthan, :<), (:conditional, :conditional_function)]);
        if ex in specials
            newex = replacements[ex];
        elseif occursin("CALLBACK_", string(ex))
            # It's a callback function. Figure out the index in callback_functions
            fun_name = string(ex)[10:end];
            for i=1:length(callback_functions)
                if callback_functions[i].name == fun_name
                    newex = :(Finch.callback_functions[$i].func);
                end
            end
        else
            newex = ex;
        end
    else
        newex = ex;
    end

    return newex;
end

function conditional_function(a, b, c)
    if a
        return b;
    else
        return c;
    end
end

# Order the factors in each term so that the test function is first and the unknown(when present) is second
function order_expression_for_fem(var, test, ex)
    # Assuming the upper levels of the tree look like t1+t2-t3+...
    # also possible: t1+(t2-t3)+... or just t1
    # Handle recursively in case there is something like t1+(t2+(t3+(t4...)))
    if typeof(ex) == Expr
        newex = copy(ex);
        if newex.head === :call
            if newex.args[1] === :+ || newex.args[1] === :.+ || newex.args[1] === :- || newex.args[1] === :.-
                # It is a + or - operation. Recursively order each term
                for i=2:length(newex.args)
                    newex.args[i] = order_expression_for_fem(var, test, newex.args[i])
                end
            elseif newex.args[1] === :* || newex.args[1] === :.*
                # search for the test function
                test_part = 0;
                var_part = 0;
                for i=2:length(newex.args)
                    if has_specific_part(newex.args[i], test)
                        test_part = i;
                    elseif has_specific_part(newex.args[i], var)
                        var_part = i;
                    end
                end

                # Place test first, var second
                if test_part > 2
                    tmp = newex.args[2];
                    newex.args[2] = newex.args[test_part];
                    newex.args[test_part] = tmp;
                    if var_part == 2
                        var_part = test_part;
                    end
                    if var_part > 3
                        tmp = newex.args[3];
                        newex.args[3] = newex.args[var_part];
                        newex.args[var_part] = tmp;
                    end

                elseif var_part > 2
                    tmp = newex.args[2];
                    newex.args[2] = newex.args[var_part];
                    newex.args[var_part] = tmp;
                end
            end

        end

        return newex;

    else # ex is not an Expr
        # It is a single symentity. Do nothing.
        return ex;
    end

end

# Returns true if a.name matches any of the parts
function has_specific_part(a, parts)
    if typeof(a) == SymEntity
        if typeof(parts) <: Array
            for i=1:length(parts)
                if has_specific_part(a, parts[i])
                    return true;
                end
            end
        elseif typeof(parts) == Variable
            return a.name == string(parts.symbol);
        elseif typeof(parts) == Coefficient
            return a.name == string(parts.symbol);
        end
    end
    return false;
end

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using Logging

if haskey(ENV, "JULIAFEM_LOGLEVEL")
    ENV["JULIAFEM_LOGLEVEL"] == "DEBUG" && Logging.configure(level=DEBUG)
end

"""
This is JuliaFEM -- Finite Element Package
"""
module JuliaFEM

importall Base

module Testing

    if VERSION >= v"0.5-"
        using Base.Test
    else
        using BaseTestNext
    end

    export @test, @testset, @test_throws

end

include("io.jl")
export Xdmf, h5file, xmffile, has_child, get_child, new_dataitem

include("fields.jl")
export Field, DCTI, DVTI, DCTV, DVTV, CCTI, CVTI, CCTV, CVTV, Increment
include("types.jl")  # data types: Point, IntegrationPoint, ...
export AbstractPoint, Point, IntegrationPoint, IP, Node

### ELEMENTS ###
include("elements.jl") # common element routines
export Node, AbstractElement, Element, update!, get_connectivity, get_basis,
       get_dbasis, inside, get_local_coordinates
include("elements_lagrange.jl") # Continuous Galerkin (Lagrange) elements
export get_reference_coordinates, get_interpolation_polynomial
export Poi1,
       Seg2, Seg3,
       Tri3, Tri6, Tri7,
       Quad4, Quad8, Quad9,
       Tet4, Tet10,
       Wedge6,
       Hex8, Hex20, Hex27

include("elements_nurbs.jl")
export NSeg, NSurf, NSolid, is_nurbs

#include("hierarchical.jl") # P-elements

include("integrate.jl")  # default integration points for elements
export get_integration_points

include("sparse.jl")
export add!, SparseMatrixCOO, get_nonzero_rows, get_nonzero_columns

include("problems.jl") # common problem routines
export Problem, AbstractProblem, FieldProblem, BoundaryProblem,
       get_unknown_field_dimension, get_gdofs, Assembly,
       get_parent_field_name, get_elements

include("problems_elasticity.jl")
export Elasticity

include("problems_dirichlet.jl")
export Dirichlet

include("problems_heat.jl")
export Heat

export assemble!, postprocess!

function assemble!(problem::Problem, element::Element, time=0.0)
    assemble!(problem.assembly, problem, element, time)
end

### ASSEMBLY + SOLVE ###
include("assembly.jl")
include("solvers_utils.jl")
include("solvers.jl")
export AbstractSolver, Solver, Nonlinear, NonlinearSolver, Linear, LinearSolver,
       get_unknown_field_name, get_formulation_type, get_problems,
       get_field_problems, get_boundary_problems,
       get_field_assembly, get_boundary_assembly,
       initialize!, create_projection, eliminate_interior_dofs,
       is_field_problem, is_boundary_problem
include("solvers_modal.jl")
export Modal

include("optics.jl")
export find_intersection, calc_reflection, calc_normal

### Mortar methods ###
include("problems_mortar.jl")
include("problems_mortar_2d.jl")
include("problems_mortar_3d.jl")
include("problems_mortar_2d_autodiff.jl")
export calculate_normals, calculate_normals!, project_from_slave_to_master,
       project_from_master_to_slave, Mortar, get_slave_elements,
       get_polygon_clip

### Mortar methods, contact mechanics extension ###
include("problems_contact.jl")
include("problems_contact_2d.jl")
include("problems_contact_3d.jl")
include("problems_contact_2d_autodiff.jl")
export Contact

module Preprocess
include("preprocess.jl")
export create_elements, Mesh, add_node!, add_nodes!, add_element!,
       add_elements!, add_element_to_element_set!, add_node_to_node_set!,
       find_nearest_nodes, reorder_element_connectivity!
include("preprocess_abaqus_reader.jl")
export parse_abaqus, parse_section, parse_element_section,
       abaqus_read_mesh, abaqus_read_model
include("preprocess_aster_reader.jl")
export aster_create_elements, parse_aster_med_file, is_aster_mail_keyword,
       parse_aster_header, aster_parse_nodes, aster_renumber_nodes!,
       aster_renumber_elements!, aster_combine_meshes, aster_read_mesh,
       filter_by_element_set, filter_by_element_id, MEDFile, aster_read_data,
       aster_read_mesh_names, aster_read_node_sets, aster_read_nodes, RMEDFile
end

function get_mesh(mesh_name::AbstractString, args...; kwargs...)
    return get_mesh(Val{Symbol(mesh_name)}, args...; kwargs...)
end

function get_model(model_name::AbstractString, args...; kwargs...)
    return get_model(Val{Symbol(model_name)}, args...; kwargs...)
end

export get_mesh, get_model

module Postprocess

include("postprocess_utils.jl")
export calc_nodal_values!, get_nodal_vector, get_nodal_dict, copy_field!,
       calculate_area, calculate_center_of_mass,
       calculate_second_moment_of_mass, extract

include("postprocess_xdmf.jl")
export XDMF, xdmf_new_result!, xdmf_save_field!, xdmf_save!, DataFrame

end
export Postprocessor

module Abaqus
include("abaqus.jl")
export abaqus_read_model, abaqus_run_model, abaqus_open_results
end

end #


text:
################################################################################
# setup
################################################################################
# using Revise
using RCall, Plots, StatPlots
# using DataBenchs
using uCSV, Feather, BenchmarkTools, CSV
using TextParse, FileIO, JLD, IterableTables
using JuliaDB, IndexedTables, uCSV, JLD2, DataFrames

using BenchmarkTools
using GR
using StatPlots
using FileIO
using Plots
using GR
gr()

outpath = "d:/tmp/"

srand(1);
N = 1_000_000;  K = 100;
# @time df = DataBench.createSynDataFrame(N, K); #31 #40
pool = "id".*dec.(1:K,3);
pool1 = "id".*dec.(1:N÷K,10);
nums = round.(rand(100).*100, 4);

############################ python package test
using PyCall
using Conda
# Conda.add("pandas") # need to run if runs into error
# Conda.add_channel("conda-forge")
# Conda.add("feather-format")
@pyimport pandas as pd
# @pyimport feather

include("benchmark/bench_read_write_df_1_code.jl")

################################################################################
# benchmark write
# ZJ: from my own testing, writing of 1m rows is sufficient to assess the
# relativities
################################################################################
println("starting testing")
if true
    julres = bench_df_write_read(1_000_000, 100, outpath, true )
    rres=rreadwrite(outpath)
    res1m = (julres, rres)
    FileIO.save(outpath*"df_write_bench_1m.jld", "res1m", res1m)

    # show read and write
    data = FileIO.load(outpath*"/df_write_bench_1m.jld")
    res1m = data["res1m"]
    (julres, rres) = res1m
    plot_bench_df_read_write(julres, rres, 1_000_000, true)
    Plots.savefig("benchmark/results/1m.png")

    wtimes = vcat([mean(a.times)/1e9 for a in julres[1]], rres[1:2:5])
    rtimes = vcat([mean(a.times)/1e9 for a in julres[2]], rres[2:2:6])

    df = DataFrame(pkg = ["Feather.jl","CSV.jl","TextParse.jl","JLD.jl","JLD2.jl", "Pandas","fst","data.table","R feather"],
    wtimes = wtimes,
    rtimes = rtimes)

    sort!(df, cols = [:rtimes])

    julres = bench_df_write_read(10_000_000, 100, outpath )
    rres=rreadwrite(outpath)
    res1m = (julres, rres)
    FileIO.save(outpath*"df_write_bench_10m.jld", "res1m", res1m)

    # show read and write
    data = FileIO.load(outpath*"/df_write_bench_10m.jld")
    res1m = data["res1m"]
    (julres, rres) = res1m
    plot_bench_df_read_write(julres, rres, 10_000_000)
    Plots.savefig("benchmark/results/10m.png")

    julres = bench_df_write_read(100_000_000, 100, outpath )
        rres=rreadwrite(outpath)
        res1m = (julres, rres)
        FileIO.save(outpath*"df_write_bench_100m.jld", "res1m", res1m)

    # show read and write
    data = FileIO.load(outpath*"/df_write_bench_100m.jld")
    res1m = data["res1m"]
    (julres, rres) = res1m
    plot_bench_df_read_write(julres, rres, 100_000_000)

    # #  show writing test
    # data = FileIO.load("e:/df_write_bench_1m_more.jld")
    # res1m = data["res1m"]
    # (julres, rres) = res1m

    # plot_bench_df_read_write(julres, rres, 1_000_000, false)

end

# if true
#     julres = bench_df_write_read(1_000_000, 100, outpath, false )
#     rres=rreadwrite(outpath)
#     res1m = (julres, rres)
#     FileIO.save("d:/tmp/df_write_bench_1m_more.jld", "res1m", res1m)
# end


text:
export LogisticBinomial, HalfCauchy, HalfNormal, EqualMix

struct HalfCauchy{T<:Real} <: Distribution{Univariate,Continuous}
    σ::T
end

HalfCauchy(σ::Integer = 1) = HalfCauchy(float(σ))

Dists.params(d::HalfCauchy) = (d.σ,)

Dists.logpdf(d::HalfCauchy, x::Real) = log(2) + logdensity(Cauchy(0, d.σ), x)

Dists.pdf(d::HalfCauchy, x) = 2 * pdf(Cauchy(0, d.σ), x)

Base.rand(rng::AbstractRNG, d::HalfCauchy) = abs(rand(rng, Cauchy(0, d.σ)))

Dists.quantile(d::HalfCauchy, p) = quantile(Cauchy(0, d.σ), (p + 1) / 2)

Dists.support(::HalfCauchy) = RealInterval(0.0, Inf)

struct HalfNormal{T<:Real} <: Distribution{Univariate,Continuous}
    σ::T
end

HalfNormal(σ::Integer = 1) = HalfNormal(float(σ))

Dists.params(d::HalfNormal) = (d.σ,)

Dists.logpdf(d::HalfNormal, x::Real) = log(2) + logdensity(Normal(0, d.σ), x)

Dists.pdf(d::HalfNormal, x) = 2 * pdf(Normal(0, d.σ), x)

Base.rand(rng::AbstractRNG, d::HalfNormal) = abs(rand(rng, Normal(0, d.σ)))

Dists.quantile(d::HalfNormal, p) = quantile(Normal(0, d.σ), (p + 1) / 2)

Dists.support(::HalfNormal) = RealInterval(0.0, Inf)

# Binomial distribution, parameterized by logit(p)
LogisticBinomial(n, x) = Binomial(n, logistic(x))

struct EqualMix{T}
    components::Vector{T}
end

Dists.logpdf(m::EqualMix, x) = logsumexp(map(d -> logdensity(d, x), m.components))

rand(m::EqualMix) = rand(GLOBAL_RNG, m)

Base.rand(rng::AbstractRNG, m::EqualMix) = rand(rng, rand(rng, m.components))

xform(d::EqualMix, _data) = xform(d.components[1], _data)

# StudentT(ν, μ = 0.0, σ = 1.0) = LocationScale(μ, σ, TDist(ν))

xform(d::Dists.Dirichlet, _data) = UnitSimplex(length(d.alpha))


text:
function combination_model(loss_vec)

    loss = sqrt(sum(loss_vec.^2))
    return loss

end


text:
export dimension
function dimension(expr :: Expr)
    # @show expr
    MLStyle.@match expr begin
        :($f |> $g)           => dimension(:($g($f)))
        # :(For($js) do $j $dist) => getTransform(:(For($j -> $dist, $js)))
        :(MixtureModel($d,$(args...))) => dimension(d)
        :(iid($n)($dist))     => n * dimension(dist)
        :(iid($n, $dist))     => n * dimension(dist)
        :(Dirichlet($k,$a))   => k
        :(Dirichlet($a))      => length(a)
        :($dist($(args...)))  => dimension(dist)
        d                     => throw(MethodError(dimension, d))
    end
end


text:
# FIXME: Move this somewhere more appropriate
convert{T <: Real, R <: Real}(::Type{SymbolicRandVar{T}}, x::R) =
  ConstantRandVar{T}(convert(T, x))

# Exponential
"exponentially distributed random variable"
exponential{T<:Real}(λ::Lift{T}, i::Id = genint()) = (-log(1-omega_component(i)))/λ

# Logistic
"Logistic exponentially distributed random variable"
logistic{T<:Real}(μ::Lift{T}, s::Lift{T}, i::Id = genint()) =
  μ + s*log(omega_component(i)/(1-omega_component(i)))

## Uniform
## =======
"univariate uniformly distributed random variable"
function uniform{T1<:Real, T2<:Real}(a::Lift{T1}, b::Lift{T2}, id::Id = genint())
  UniformRandVar{Float64, Float64, Float64}(id,
    convert(SymbolicRandVar{Float64},a),
    convert(SymbolicRandVar{Float64},b))
end

"Standard uniform a = 0.0, b = 1.0"
uniform(i::Id = genint()) = uniform(0.0, 1.0, i)

## Normal
## ======
# "Constructs Normally distributed random variable constructor"
# function normal{T1<:Real, T2<:Real}(μ::Lift{T1}, σ::Lift{T2}, id::Id = genint())
#   NormalRandVar{Float64, Float64, Float64}(id,
#     convert(SymbolicRandVar{Float64},μ),
#     convert(SymbolicRandVar{Float64},σ))
# end
#
# "Standard normal μ = 0.0 σ = 1.0"
# normal(i::Id = genint()) = normal(0.0, 1.0, i)

## Discrete Distributions
## ======================

# FIXME: Why is this causing errors
# "Bernoulli distributed random variable"
# function flip{T<:Real}(p::Lift{T}, id::Int = genint())
#   BernoulliRandVar{Bool, Float64}(id, convert(SymbolicRandVar{Float64}, p))
# end

"Bernoulli distributed random variable"
flip{T<:Real}(p::Lift{T}, id::Int = genint()) = p >= omega_component(id)

"Standard Bernoulli p = 0.5"
flip(id::Id = genint()) = flip(0.5, id)

# "Poisson distributed random variable constructor"
# function poisson{T<:Real}(λ::SymbolicRandVar{T}, id::Id = genint())
#   PoissonRandVar{Int64, Float64}(id, convert(SymbolicRandVar{Float64}, λ))
# end
#
# "Standard Poisson λ = 1.0"
# poisson(id::Id = genint()) = poisson(1.0, id)

"Poisson distributed random variable constructor"
function discreteuniform(a::Lift{Int64}, b::Lift{Int64}, id::Id = genint())
  DiscreteUniformRandVar{Int64, Int64, Int64}(id,
    convert(SymbolicRandVar{Int64}, a),
    convert(SymbolicRandVar{Int64}, b))
end


text:
module ContextAware

using Agents.BaseAgent
using Params
using Traffic.Simple
using Channel.Gilbert
using Util
using Params.QParams

export ContextQ

# Q-Learning state, may contains some fields in
# AgentState
type StateT
    chan :: Int
    buflevel :: Int
    energysaving :: EnergySavingMode
end

type ContextQ <: Agent
    s :: AgentState
    Q :: Array{Float64, 4}
    visit :: Array{Float64, 4}
    a :: Int
    P_tx :: Float64
    bitrate :: Float64
    status :: Status
    expertness :: Float64
    buf_interval :: Int
    beta_idle :: Float64
    state :: StateT
    δ :: Float64
    trustQ :: Float64
    buflevels :: Int # needed for bufer-level-baed computation in feedback
end

const n_p_levels = length(Params.P_levels)
const idle_action = Params.n_channel * n_p_levels + 1

function ContextQ(i, P, pos)
    Q = rand(EnergySaving.n, int(Params.n_channel), P.buf_levels + 1, idle_action)
    Q *= Params.P_tx * Params.t_slot # a good initial randomization
    visit = zeros(Int, EnergySaving.n, Params.n_channel, P.buf_levels + 1, idle_action)
    agent = ContextQ(AgentState(i, P, pos), Q, visit, 0, 0, 0, Initialized, 0, div(Params.B + 1, P.buf_levels), P.beta_idle, StateT(0, 0, MaxThroughput), P.δ, P.trustQ, P.buf_levels + 1)
    agent.state.energysaving = agent.s.energysaving
    agent
end

function policy!(a :: ContextQ)
    if rand() < epsilon
        a.a = rand(1:idle_action)
    else
        a.a = indmax(a.Q[a.state.energysaving.n, a.s.chan, div(a.s.B_empty, a.buf_interval) + 1, :])
    end
    nothing
end

function BaseAgent.act(a :: ContextQ, env, t)
    # check whether we have packets first
    if a.status == Initialized
        if a.s.B_max - a.s.B_empty == 0
            a.a = -1
            return idle(a)
        end

        policy!(a)
        # infer state
        a.state.chan = a.s.chan
        a.state.buflevel = div(a.s.B_empty, a.buf_interval) + 1
        a.state.energysaving = a.s.energysaving
        # update state visit count
        a.visit[a.state.energysaving.n, a.state.chan, a.state.buflevel, a.a] += 1

        if a.a == idle_action
            return idle(a)
        end

        chan = int8(fld(a.a - 1, n_p_levels) + 1)
        return switch(a, chan)
    elseif a.status == Switched
        return sense(a, env, detect_traffic)
    elseif a.status == Sensed
        if detect_traffic(a.s, env.traffics[a.s.chan], a.s.t_remaining)
            return idle(a)
        else
            a.P_tx = Params.P_levels[(a.a - 1) % n_p_levels + 1]
            a.bitrate = capacity(env.channels[a.s.chan], a.P_tx, a.s.pos)
            pkgs_to_send = min(a.s.B_max - a.s.B_empty, floor(a.s.t_remaining * a.bitrate / Params.pkt_size))

            if pkgs_to_send == 0
                return idle(a)
            end
            return transmit!(a, a.P_tx, env, pkgs_to_send)
        end
    elseif a.status == Transmitted
        return idle(a)
    end
end

function alpha(a :: ContextQ)
    return 0.4 + 0.6 / (1 + a.visit[a.state.energysaving.n, a.state.chan, a.state.buflevel, a.a])
end

function BaseAgent.feedback(a :: ContextQ, res :: Result, idle :: Bool = false, n_pkt :: Int = 0)
    feedback(a.s, res, idle, n_pkt)
    # if we didn't take any action that's worth learning, return immediately
    if a.a == -1
        return nothing
    end
    r :: Float64 = 0
    if idle
        r = - a.beta_idle * a.bitrate * Params.t_slot / a.s.E_slot
        if res == BufOverflow # If we were idle and overflow occurred, get some extra punishment
            r = r * Params.beta_overflow / a.beta_idle
        else
            if a.state.energysaving == MaxThroughput
                r *= 2 # double beta_idle in max throughput mode
            end
            new_buf_level = div(a.s.B_empty, a.buf_interval) + 1
            r *= 1 + 0.5 * (new_buf_level - a.buflevels / 2) / a.buflevels # increase idleness penalty with buffer level
        end
    elseif res == Success
        K = 1 # a.P_tx ^ 2 * Params.t_slot / a.bitrate
        r = K * Params.pkt_size * n_pkt / a.s.E_slot
    elseif res == Collision
        r = - beta_md * a.bitrate * Params.t_slot / a.s.E_slot
    elseif res == LostInChannel
        r = - beta_loss * a.bitrate * Params.t_slot / a.s.E_slot
    else
        return nothing
    end
    # UPDATE Q
    max_Q_next = maximum(a.Q[a.state.energysaving.n, a.s.chan, floor(a.s.B_empty / a.buf_interval) + 1])
    Q_now = a.Q[a.state.energysaving.n, a.state.chan, a.state.buflevel, a.a]
    a.Q[a.state.energysaving.n, a.state.chan, a.state.buflevel, a.a] += alpha(a) * (r + discount * max_Q_next - Q_now)
    # Update expertness
    if r > 0
        a.expertness += r # sum of positive reinforcements for now
    end
    nothing
end

# Coordinator object for Cooperative Q
type Coordinator
    US :: Array{Float64, 3} # U*S where U*S*V' = svd(Q) for each agent
    Vt :: Array{Float64, 3} # V' where U*S*V' = svd(Q) for each agent
    function Coordinator(P :: ParamT)
        c = new()
        c.US = zeros(P.n_agent, int(Params.n_channel) * (P.buf_levels + 1), P.d_svd)
        c.Vt = zeros(P.n_agent, P.d_svd, Params.n_channel * length(Params.P_levels) + 1)
    end
end

function BaseAgent.initcoordinator(:: Type{ContextQ}, P :: ParamT)
    return Coordinator(P)
end

function expertweights(a, expertness, agents, i)
    #expertness -= min(expertness)
    # convert weights to values in [0, 1]
    distances = [norm([a.s.x, a.s.y] + rand(a.s.location_error)) for a in agents]
    distances = abs(distances - distances[i]) + 5 # +5 part is for accuracy
    expertness ./= distances .^ a.δ
    expertness .*= Params.QParams.cap_weights[[a.cap_level for a in agents]]
    weights = logistic(expertness, mean(expertness), 1, expertness[i])
    # introduce distances into weights
    #temp = weights[i]
    #weights[i] = 0
    # weights[expertness .<= expertness[i]] = 0
    # make weights one-sum
    weights = weights ./ sum(weigths)
    # normalize weights to use trustQ
    weights .*= a.trustQ / sum(weights[1:endof(weights) .!= i])
    weights[i] = 1 - a.trustQ
end

function BaseAgent.cooperate(agents :: Vector{ContextQ}, P :: ParamT, coordinator :: Coordinator, t)
    const n_agent = int(P.n_agent)
    const sharingperiod = P.sharingperiod
    # size of US + size of Vt, used for data sharing time and energy computation
    const sizeQ = P.d_svd * 64 * (EnergySaving.n * n_channel * (P.buf_levels + 1) + (n_channel * length(P_levels) + 1))
    # time slots required for an agent to send/receive Q matrix by sending/receiving US & Vt
    const timeQ = int(ceil(sizeQ ./ Params.controlcapacity ./ t_slot))
    # time required for an agent to send/receive Q matrix by sending/receiving US & Vt
    const rawtimeQ0 = sizeQ ./ Params.controlcapacity ./ t_slot
    # Shape of the Q matrix (Shape is S x A where S and A are discrete enumerated sets)
    const shapeQ = (EnergySaving.n * int(Params.n_channel) * (P.buf_levels + 1), Params.n_channel * length(Params.P_levels) + 1)
    if n_agent > 1
        US = coordinator.US
        Vt = coordinator.Vt
        # join Cooperative Q learners
        if t > Params.t_saturation && t % sharingperiod < 2 * n_agent * timeQ
            tt = t % sharingperiod
            i = fld(tt, timeQ) + 1
            # sending phase (from CR)
            if tt < n_agent * timeQ
                # receiving complete, update combined Q matrix
                # using only positive expertness ones
                if tt - (i-1) * timeQ == timeQ - 1
                    u, s, v = svd(reshape(agents[i].expertness .* agents[i].Q, shapeQ))
                    US[i, :, :] = u[:,1:P.d_svd]*diagm(s[1:P.d_svd])
                    Vt[i, :, :] = v'[1:P.d_svd, :]
                    expertness[i] = agents[i].expertness
                end
                if tt == (i - 1) * timeQ
                    # initialize control channel usage overhead
                    rawtimeQ = rawtimeQ0
                end
                # Calculate control channel energy overhead
                if rawtimeQ > 1
                    agents[i].s.E_slot += Params.P_tx * t_slot # add control channel energy overhead
                else
                    agents[i].s.E_slot += Params.P_tx * t_slot * rawtimeQ # add control channel energy overhead
                end
                rawtimeQ -= 1 # doesn't matter when rawtimeQ < 1, it'll be reset anyways
            else # receiving phase (to CR)
                i -= n_agent
                # sending complete, update agent's Q matrix
                if tt - (i+n_agent-1) * timeQ == timeQ - 1
                    #fill!(agents[i].Q, 0)
                    # use up-to-date Q function on CR's side
                    # learning from experts (LE)
                    weights = expertweights(agents[i], expertness, agents, i)
                    agents[i].Q *= weights[i]
                    for j=1:n_agent
                        if i == j || weights[j] == 0 # || !experts[j]
                            continue
                        else
                            us = slice(US, (j, 1:size(US)[2], 1:size(US)[3]))
                            vt = slice(Vt, (j, 1:size(Vt)[2], 1:size(Vt)[3]))
                            agents[i].Q += weights[j] * reshape(us * vt, size(agents[i].Q))
                            fill!(agents[i].visit, 0)
                        end
                    end
                    agents[i].expertness = 0 # *= 1 - a.trustQ
                end
            end
        end
    end
end

end


text:
#    Main Module of the FourierAnalysis Package for julia language
#    v 1.2.0 - last update September 2021
#
#   MIT License
#   Copyright (c) 2019-2022,
#   Marco Congedo, CNRS, Grenobe, France:
#   https://sites.google.com/site/marcocongedo/home

# __precompile__()

module  FourierAnalysis

using   Base.Threads,
        LinearAlgebra,
        Statistics,
        FFTW,
        AbstractFFTs,
        DSP,
        RecipesBase

# # # # # Special instructions and variables # # # # #
# N.B. multi-threaded planners are not yet exported by FFTW
BLAS.set_num_threads(Sys.CPU_THREADS)

# # # # # constants # # # # #
#const qroot2        = 2^0.25
const titleFont     = "\x1b[36m"
const separatorFont = "\x1b[33;1m"
const defaultFont   = "\x1b[0m"
const greyFont      = "\x1b[90m"
const 📌            = titleFont*"FourierAnalysis"*defaultFont
# # # # # types # # # # #

IntOrReal    = Union{Int, Real}
RealOrComplex= Union{Real, Complex}

# Frequency ranges (in Hz) for the mean function applying to (Cross)Spcetra objects
fInterval    = Union{IntOrReal, Tuple{IntOrReal, IntOrReal}, Colon}

# time ranges (in samples) for the mean function applying to TFx objects
tInterval    = Union{Int, Tuple{Int, Int}, Colon}

# smoother kind; see tools.jl
@enum Smoother begin
    noSmoother       = 1
    hannSmoother     = 2
    hammingSmoother  = 3
    blackmanSmoother = 4
end

# see spectra.jl
SpectraData      = AbstractArray{T} where T<:Real
struct Spectra
    y           :: SpectraData  # spectra
    sr          :: Int          # sampling rate
    wl          :: Int          # FFT window length
    DC          :: Bool         # DC component in first bin if true, not present otherwise
    taper       :: String       # time-domain tapering type as a string (with parameters for Slepians)
    flabels     :: Vector{T} where T<:Union{Real, Int} # Fourier Discrete Frequencies
    func        :: Function
    smoothing   :: Smoother # enum type, see above
end
SpectraVector   = Vector{Spectra}

# see crossspectra.jl
CrossSpectraData = Union{Vector{LowerTriangular}, Vector{Hermitian}}
struct CrossSpectra
    y           :: CrossSpectraData
    sr          :: Int          # sampling rate
    wl          :: Int          # FFT window length
    DC          :: Bool         # DC component in first bin if true, not present otherwise
    taper       :: String       # time-domain tapering type as a string (with parameters for Slepians)
    flabels     :: Vector{T} where T<:Union{Real, Int} # Fourier Discrete Frequencies
    nonlinear   :: Bool
    smoothing   :: Smoother
    tril        :: Bool
end
CrossSpectraVector = Vector{CrossSpectra}

# see coherence.jl
CoherenceData = Union{Vector{LowerTriangular}, Vector{Hermitian}}
struct Coherence
    y           :: CoherenceData
    sr          :: Int          # sampling rate
    wl          :: Int          # FFT window length
    DC          :: Bool         # DC component in first bin if true, not present otherwise
    taper       :: String       # time-domain tapering type as a string (with parameters for Slepians)
    flabels     :: Vector{T} where T<:Union{Real, Int} # Fourier Discrete Frequencies
    nonlinear   :: Bool
    smoothing   :: Smoother
    tril        :: Bool
end
CoherenceVector = Vector{Coherence}
CoherenceVector₂ = Vector{CoherenceVector}

# NB data field `y` in CrossSpectra and Coherence structures is one of these types:
# Array{LowerTriangular,1}, Array{Hermitian,1}
# CrossSpectraVector and CoherenceVector are of one of these types:
# Array{Array{LowerTriangular,1},1}, Array{Array{Hermitian,1},1}
# All these types are straightforwardly compatible with package PosDefMaifold.

# see timefrequency.jl
TFAnalyticSignalData = Matrix{T} where T<:Complex
struct TFAnalyticSignal
    y          :: TFAnalyticSignalData # Analytic Signal, dim1=freq, dim2=time
    bandwidth  :: IntOrReal
    flabels    :: Vector{S} where S<:Real # center frequency axis labels (in Hz)
    nonlinear  :: Bool # true if nonlinear, false otherwise
    fsmoothing :: Smoother
    tsmoothing :: Smoother
end
TFAnalyticSignalVector = Vector{TFAnalyticSignal}

TFAmplitudeData = Matrix{T} where T<:Real
struct TFAmplitude
    y          :: TFAmplitudeData # Analytic (Instantaneous) Amplitude, dim1=freq, dim2=time
    bandwidth  :: IntOrReal
    flabels    :: Vector{S} where S<:Real # center frequency axis labels (in Hz)
    fsmoothing :: Smoother
    tsmoothing :: Smoother
    func       :: Function
end
TFAmplitudeVector = Vector{TFAmplitude}

TFPhaseData = Matrix{T} where T<:Real
struct TFPhase
    y          :: TFPhaseData # Analytic (Instantaneous) Phase, dim1=freq, dim2=time
    bandwidth  :: IntOrReal
    flabels    :: Vector{S} where S<:Real # center frequency axis labels (in Hz)
    nonlinear  :: Bool # true if nonlinear, false otherwise
    fsmoothing :: Smoother
    tsmoothing :: Smoother
    unwrapped  :: Bool
    func       :: Function
end
TFPhaseVector = Vector{TFPhase}

# all frequency domain objects
FDobjects=Union{Spectra, CrossSpectra, Coherence}
FDobjectsVector=Union{SpectraVector, CrossSpectraVector, CoherenceVector}
# all Time-frequency domain objects
TFobjects=Union{TFAnalyticSignal, TFAmplitude, TFPhase}
TFobjectsVector=Union{TFAnalyticSignalVector, TFAmplitudeVector, TFPhaseVector}

#import DSP: dpss
import Statistics.mean
import Base: show, *, conj, real, imag #,
import DSP: coherence

export

# From this module
fInterval,
tInterval,
Smoother,
    noSmoother,
    hannSmoother,
    hammingSmoother,
    blackmanSmoother,

SpectraType,            Spectra,            SpectraVector,
CrossSpectraType,       CrossSpectra,       CrossSpectraVector,
CoherenceType,          Coherence,          CoherenceVector,
TFAnalyticSignalType,   TFAnalyticSignal,   TFAnalyticSignalVector,
TFAmplitudeType,        TFAmplitude,        TFAmplitudeVector,
TFPhaseType,            TFPhase,            TFPhaseVector,

FDobjects, FDobjectsVector,
TFobjects, TFobjectsVector,

# from fftw.jl
plan_estimate,
plan_measure,
plan_patient,
plan_exhaustive,
plan_conserve_memory,
plan_wisdom_only,
plan_unaligned,
plan_preserve_input,
Planner,

# From tapers.jl
TaperKind,
    rectangular,
    triangular,
    hann,
    hamming,
    blackman,
    harris4,
    riesz,
    parzen,
    slepian,
taper,
slepians,
taperinfo,

# from tools.jl
sinusoidal,
fres,
f2b,
b2f,
fdf,
brange,
bbands,
fbands,
decibel,
amplitude,
phase,
polar,
unwrapPhase,
smooth,
extract, extr,
mean,
bands,
sameParams,
isLinear,
isNonLinear,
isUnwrapped,

# from goertzel.jl
goertzel,
goertzel_fast,
goertzel2,

# from spectra.jl
spectra,

# from crossspectra.jl
crossSpectra,

# from coherence.jl
coherence,

# from hilbert.jl
analyticsignal,

# from filters.jl
FilterDesign,
filterbank,

# from timefrequency.jl
TFanalyticsignal,
TFamplitude,
TFphase,

# from timefrequencyuni.jl
meanAmplitude, mamp,
concentration, con,
meanDirection, mdir,

# from timefrequencybi.jl
comodulation, com,
coherence, coh,

# from recipes.jl
tfAxes

include("fftw.jl")
include("tapers.jl")
include("tools.jl")
include("goertzel.jl")
include("spectra.jl")
include("crossspectra.jl")
include("coherence.jl")
include("hilbert.jl")
include("filters.jl")
include("timefrequency.jl")
include("timefrequencyuni.jl")
include("timefrequencybi.jl")
include("recipes.jl")

# welcome message
println("\n⭐ "," Welcome to the ", 📌," package", " ⭐\n")
@info(" ")
println(" Your Machine `",gethostname(),"` (",Sys.MACHINE, ")")
println(" runs on kernel ",Sys.KERNEL," with word size ",Sys.WORD_SIZE,".")
println(" CPU  Threads: ",Sys.CPU_THREADS)
# Sys.BINDIR # julia bin directory
println(" Base.Threads: ", "$(Threads.nthreads())")
println(" BLAS Threads: ", "$(Sys.CPU_THREADS)", "\n")
end # module


text:
using ITensorGaussianMPS
using LinearAlgebra
using ITensors

# Half filling
N = 50
Nf = N÷2

@show N, Nf

# Maximum MPS link dimension
_maxlinkdim = 100

@show _maxlinkdim

# DMRG cutoff
_cutoff = 1e-12

# Hopping
t = 1.0

# Electron-electron on-site interaction
U = 1.0

@show t, U

# Free fermion Hamiltonian
ampo = AutoMPO()
for n in 1:N-1
  ampo .+= -t, "Cdag", n, "C", n+1
  ampo .+= -t, "Cdag", n+1, "C", n
end

# Hopping Hamiltonian with N spinless fermions
h = hopping_hamiltonian(ampo)

# Get the Slater determinant
Φ = slater_determinant_matrix(h, Nf)

# Create an mps for the free fermion ground state
s = siteinds("Fermion", N; conserve_qns = true)
println("Making free fermion starting MPS")
@time ψ0 = slater_determinant_to_mps(s, Φ; eigval_cutoff = 1e-4,
                                           cutoff = _cutoff,
                                           maxdim = _maxlinkdim)
@show maxlinkdim(ψ0)

# Make an interacting Hamiltonian
for n in 1:N-1
  ampo .+= U, "N", n, "N", n+1
end
H = MPO(ampo, s)

# Random starting state
ψr = randomMPS(s, n -> n ≤ Nf ? "1" : "0")

println("\nRandom state starting energy")
@show flux(ψr)
@show inner(ψr, H, ψr)

println("\nFree fermion starting energy")
@show flux(ψ0)
@show inner(ψ0, H, ψ0)

println("\nRun dmrg with random starting state")
sweeps = Sweeps(20)
maxdim!(sweeps,10,20,40,_maxlinkdim)
cutoff!(sweeps,_cutoff)
@time dmrg(H, ψr, sweeps)

println("\nRun dmrg with free fermion starting state")
sweeps = Sweeps(4)
maxdim!(sweeps,_maxlinkdim)
cutoff!(sweeps,_cutoff)
@time dmrg(H, ψ0, sweeps)

nothing


text:
module QMPS

# using Yao, Yao.ConstGate
using Yao
# using StatsBase

include("Diff.jl")
include("MPSC.jl")
include("CircuitBuilder.jl")

end


text:
#!/usr/bin/env julia --project="../"
using RobotOS
@rosimport std_msgs.msg.String
rostypegen()
using .std_msgs.msg: StringMsg
function talker()
    println("Inside talker")
    init_node("talker") # node name
    pub = Publisher{StringMsg}("chatter", queue_size=10) # topic name
    rate = Rate(10) # 10 Hz
    while !is_shutdown()
        hello_str = "hello world $(to_sec(get_rostime()))"
        loginfo(hello_str)
        publish(pub, StringMsg(hello_str))
        rossleep(rate)
    end
end
if !isinteractive()
    println("Running in non-interactive mode")
    talker()
end

max_stars_repo_name: mattborghi/ros-julia-examples

max_issues_repo_name: mattborghi/ros-julia-examples

max_forks_repo_name: mattborghi/ros-julia-examples


text:
module FGES

#Packages + Reason required
using LinearAlgebra #for subarrays (?)
using Statistics #mean, std
using Random #set seeds in generation of random graph
using Combinatorics #loop through combinations and permutations
using Memoization, LRUCache #Caching scores
using Dates #record the current time
using DataFrames

#Base functions to overload
import Base: show, view, iterate, eltype, length, == #functions from base to extend

#Files to include
include("graphStructure.jl")
include("graphAlgorithms.jl")
include("mainAlgorithm.jl")

export

#graphStructure.jl
    PDAG,
    Edge,

    nv,
    ne,

    isneighbor,
    isadjacent,
    isparent,
    ischild,
    isdescendent,
    isclique,
    isblocked,
    isdirected,

    addedge!,
    remedge!,

    neighbors,
    neighbors_in,
    neighbors_out,
    neighbors_undirect,
    parents,
    children,
    descendent,

    countNeighbors,
    countNeighbors_in,
    countNeighbors_out,
    countNeighbors_undirect,
    countParents,
    countChildren,

    allundirected,
    alldirected,
    alledges,
    allpairs,

    edges,
    vertices,

    orientedge!,
    calcNAyx,
    calcT,
    degreeAverage,

    saveGraph,
    loadGraph,
    edgetable,

#graphAlgorithms.jl
    graphVStructure!,
    categorizeNeighbors,
    setCategory!,
    meekRules!,

#mainAlgorithm.jl
    fges,
    ParseData,
    Step,
    Search!,
    Insert!,
    Delete!,
    findNextEquivClass!,
    findBestInsert,
    findBestDelete,
    statusUpdate

end


text:
using SciMLBase
using SteadyStateDiffEq

# Test find_steady_states
res = find_steady_states(trajectories=10)

@test res isa EnsembleSolution
@test res[1] isa SteadyStateSolution
@test length(res) > 0


text:
push!(LOAD_PATH,"../src")
using Documenter: Documenter, makedocs, deploydocs, doctest, DocMeta
using CombinedParsers
using CombinedParsers.Regexp
using Test
using Literate

docdir = joinpath(dirname(pathof(CombinedParsers)),"../docs/src/")
mandir = joinpath(docdir,"man")
DocMeta.setdocmeta!(CombinedParsers, :DocTestSetup, quote
    using TextParse
    using CombinedParsers
    using CombinedParsers.Regexp
    using CombinedParsers.BNF
end; recursive=true)

for f in [ "pcre.jl", "example-number-ranges.jl", "example-either-trie.jl" ]
    Literate.markdown(joinpath(mandir,f), mandir,
                      repo_root_url="https://github.com/gkappler/CombinedParsers.jl/docs",
                      codefence = "```@repl session" => "```")
end

for f in [ "example-palindromes.jl", "bson.jl"]
    Literate.markdown(joinpath(mandir,f), mandir,
                      repo_root_url="https://github.com/gkappler/CombinedParsers.jl/docs")
end

makedocs(;
         source=docdir,
         modules=[CombinedParsers],
         authors="Gregor Kappler",
         repo="https://github.com/gkappler/CombinedParsers.jl/blob/{commit}{path}#L{line}",
         sitename="CombinedParsers.jl",
         format=Documenter.HTML(;
                                prettyurls=get(ENV, "CI", "false") == "true",
                                canonical="https://gkappler.github.io/CombinedParsers.jl",
                                assets=String[],
                                ),
         pages=[
             "Home" => "index.md",
             "Manual" => [
                 "Overview" => "man/guide.md",
                 "User Guide" => "man/user.md",
                 "Prefix-Tree Matching" => "man/example-either-trie.md",
                 "PCRE Compliance" => "man/pcre-compliance.md",
             ],
             "Examples" => [
                 "Names (What is Parsing?)" => "man/example-person.md",
                 "Number sequences (What are Representations?)" => "man/example-number-ranges.md",
                 "Arithmetics (Evaluation)" => "man/example-arithmetics.md",
                 "Regular Expressions" => "man/pcre.md",
                 "JSON (recursion)" => "man/json.md",
                 "BSON (`Vector{Unit8}` parsing)" => "man/bson.md",
                 "struct Palindrome<:CombinedParser" => "man/example-palindromes.md",
             ],
             "Library" => Any[
                 "Matching/Parsing" => "lib/public.md",
                 "Templates" => "lib/parsers.md",
                 "Constructors" => "lib/constructors.md",
                 "Regexp" => "lib/regexp.md",
                 "EBNF" => "man/bnf.md",
                 "Transformations" => "lib/transformation.md",
                 "Internals" => "lib/internals.md"
             ],
             # "Developer Guide" => "developer.md"
         ],
         )

deploydocs(;
    repo="github.com/gkappler/CombinedParsers.jl",
)


text:
immutable TrivialModel{T}
    α::T
    β::T
end

@testset "Collocation" begin
    fam = Chebyshev(10)
    θ = fit(fam, x->x^2)
    model = TrivialModel(2.0, 3.0)
    residual_function(model::TrivialModel, f, x) = f(x-model.α)-model.β
    cres = CollocationResidual(model, fam, residual_function)
    r = cres(θ)
    @test r ≈ (points(fam) - model.α).^2-model.β
end

immutable TrivialModel2{T}
    α::T
end

@testset "Collocation solver" begin
    (m::TrivialModel2)(x) = m.α*x
    model0 = TrivialModel2(2.0)
    residualf(model::TrivialModel2, f, x) = f(x)-model(x)
    res = CollocationResidual(model0, Chebyshev(10), residualf)
    f_sol, o = solve_collocation(res, zero)
    @test norm([f_sol(x)-model0(x) for x in linspace(domain(f_sol), 100)]) ≤ 1e-14
end