text:
using URIParser

jlhome() = ccall(:jl_get_julia_home, Any, ())

function basepath(files...)
  srcdir = joinpath(jlhome(), "..", "..")
  releasedir = joinpath(jlhome(), "..", "share", "julia")
  normpath(joinpath(isdir(srcdir) ? srcdir : releasedir, files...))
end

function pkgrootpath(pkg)
  pkgpath = Base.find_package(pkg)

  # package not installed
  isfile(pkgpath) || return nothing

  return normpath(joinpath(dirname(pkgpath), ".."))
end

"""
    docsdir(pkg) -> String

Find the directory conataining the documenatation for package `pkg`. Will fall back to
returning a package's README.md. Returns an empty `String` if no indication of documentation
is found.
"""
function docsdir(pkg)
  # sepcial case base
  lowercase(pkg) == "base" && return joinpath(basepath("doc"), "src")

  pkgpath = pkgrootpath(pkg)

  pkgpath === nothing && return ""

  # Documenter.jl default:
  docpath = joinpath(pkgpath, "docs", "src")
  isdir(docpath) && return docpath

  # other possibility:
  docpath = joinpath(pkgpath, "doc", "src")
  isdir(docpath) && return docpath

  # fallback to readme
  readmepath = joinpath(pkgpath, "README.md")
  return isfile(readmepath) ? readmepath : ""
end

function readmepath(pkg)
  (lowercase(pkg) == "base") && return ""

  pkgpath = pkgrootpath(pkg)
  # package not installed
  pkgpath === nothing && return ""

  joinpath(pkgpath, "README.md")
end

"""
    docsurl(pkg) -> String

Return the most likely candidate for a package's online documentation or an empty string.
"""
docsurl(pkg) = baseURL(finddocsURL(pkg))

"""
    baseURL(links::Vector{Markdown.Link}) -> String

Find the most common host and return the first URL in `links` with that host.
"""
function baseURL(links::Vector{Markdown.Link})
  isempty(links) && return ""

  length(links) == 1 && return links[1].url

  # find most common host
  urls = map(x -> URI(x.url), links)
  hosts = String[url.host for url in urls]
  perm = sortperm([(host, count(x -> x == host, hosts)) for host in unique(hosts)], lt = (x,y) -> x[2] > y[2])

  # TODO: better heuristic for choosing the right path
  links[perm[1]].url
end

"""
    finddocsURL(pkg) -> Vector{Markdown.Link}

Search `pkg`s readme for links to documentation.
"""
function finddocsURL(pkg)
  lowercase(pkg) == "base" && return [Markdown.Link("", "https://docs.julialang.org")]
  pkgpath = pkgrootpath(pkg)

  doclinks = Markdown.Link[]
  pkgpath === nothing && return doclinks

  readmepath = joinpath(pkgpath, "README.md")
  isfile(readmepath) || return doclinks

  md = Markdown.parse(String(read(joinpath(pkgpath, "README.md"))))
  links = findlinks(md)

  isempty(links) && (links = findplainlinks(md))

  for link in links
    if isdoclink(link)
      push!(doclinks, link)
    end
  end
  doclinks
end

function findplainlinks(md)
  text = Markdown.plain(md)
  [Markdown.Link(link, link) for link in matchall(r"(https?:\/\/[^\s]+)\b", text)]
end

function isdoclink(link::Markdown.Link)
  p = lowercase(Markdown.plaininline(link))
  # TODO: could be a bit smarter about this
  contains(p, "docs") || contains(p, "documentation") ||
    contains(p, "/stable") || contains(p, "/latest")
end

function findlinks(mdobj)
  doclinks = Markdown.Link[]
  for obj in mdobj.content
    findlinks(obj, doclinks)
  end
  doclinks
end

function findlinks(mdobj::Markdown.Paragraph, links)
  for obj in mdobj.content
    findlinks(obj, links)
  end
end

findlinks(mdobj, links) = nothing
findlinks(mdobj::Markdown.Link, links) = push!(links, mdobj)


text:
@testset "macro methods" begin
    @testset "tensorexpressions" begin
        using TensorOperations: isassignment, isdefinition, getlhs, getrhs, isindex,
            istensor, isgeneraltensor, istensorexpr, isscalarexpr, hastraceindices,
            hastraceindices, getindices, getallindices,
            normalizeindex, instantiate_scalar, instantiate_eltype,
            decomposetensor, decomposegeneraltensor

        @test isassignment(:(a[-1,-2,-3] = b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isassignment(:(a[-1,-2,-3] += b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isassignment(:(a[-1,-2,-3] -= b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isassignment(:(a[-1,-2,-3] := b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isassignment(:(a[-1,-2,-3] ≔ b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isassignment(:(b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))

        @test !isdefinition(:(a[-1,-2,-3] = b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isdefinition(:(a[-1,-2,-3] += b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isdefinition(:(a[-1,-2,-3] -= b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isdefinition(:(a[-1,-2,-3] := b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test isdefinition(:(a[-1,-2,-3] ≔ b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))
        @test !isdefinition(:(b[-1,-2,1]*c[1,-3]+d[-2,-3,-1]))

        lhs = :(a[-1,-2,-3])
        rhs = :(b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])
        for (getside,side) in ((getlhs, lhs), (getrhs, rhs))
            @test getside(:(a[-1,-2,-3] = b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == side
            @test getside(:(a[-1,-2,-3] += b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == side
            @test getside(:(a[-1,-2,-3] -= b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == side
            @test getside(:(a[-1,-2,-3] := b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == side
            @test getside(:(a[-1,-2,-3] ≔ b[-1,-2,1]*c[1,-3]+d[-2,-3,-1])) == side
        end

        @test isindex(:a)
        @test normalizeindex(:a) == :a
        @test isindex(:(a'))
        @test normalizeindex(:(a')) == :(a′)
        @test isindex(:(a'''))
        @test normalizeindex(:(a''')) == :(a′′′)
        @test isindex(:(β))
        @test normalizeindex(:(β)) == :(β)
        @test isindex(:(β'))
        @test normalizeindex(:(β')) == :(β′)
        @test isindex(:(3))
        @test normalizeindex(:(3)) == 3
        @test isindex(:(-5))
        @test normalizeindex(:(-5)) == -5
        @test !isindex(:('a'))
        @test !isindex(:(a+b))
        @test !isindex(:("x"))
        @test !isindex(:(5.1))

        @test istensor(:(a[1,2,3]))
        @test decomposetensor(:(a[1,2,3])) == (:a, Any[1,2,3], Any[])
        @test istensor(:(a[5][a b c]))
        @test decomposetensor(:(a[5][a b c])) == (:(a[5]), Any[:a,:b,:c], Any[])
        @test istensor(:(cos(y)[a b c; 1 2 3]))
        @test decomposetensor(:(cos(y)[a b c; 1 2 3])) == (:(cos(y)), Any[:a,:b,:c], Any[1,2,3])
        @test istensor(:(x[(); (1,2,3)]))
        @test decomposetensor(:(x[(); (1,2,3)])) == (:x, Any[], Any[1,2,3])
        @test !istensor(:(2*a[1,2,3]))
        @test !istensor(:(a[1 2 3; 4 5 6; 7 8 9]))
        @test !istensor(:(conj(a[5][a b c])))
        @test !istensor(:(cos(y)[a b c; 1 2 3]*b[4,5]))
        @test !istensor(:(3+5))

        @test isgeneraltensor(:(conj(a[1,2,3])))
        @test decomposegeneraltensor(:(conj(a[1,2,3]))) == (:a, [1,2,3], [], instantiate_scalar(:(conj(1))), true)
        @test isgeneraltensor(:(x*a[5][a b c]))
        @test decomposegeneraltensor(:(x*a[5][a b c])) == (:(a[5]), [:a,:b,:c], [], instantiate_scalar(:(x*1)), false)
        @test isgeneraltensor(:(3*conj(a*cos(y)[a b c; 1 2 3])))
        @test decomposegeneraltensor(:(3*conj(a*cos(y)[a b c; 1 2 3]))) == (:(cos(y)),  Any[:a,:b,:c],  Any[1,2,3], instantiate_scalar(:(3*conj(a*1))), true)
        @test !isgeneraltensor(:(1/a[1,2,3]))
        @test !isgeneraltensor(:(a[1 2 3; 4 5 6]\x))
        @test !isgeneraltensor(:(cos(y)[a b c; 1 2 3]*b[4,5]))
        @test !isgeneraltensor(:(3+5))

        @test hastraceindices(:(a[x,y,z,1,x]))
        @test hastraceindices(:(a[x y z;1 x]))
        @test !hastraceindices(:(a[-1,y,z,1,x]))
        @test !hastraceindices(:(a[-1 y z;1 x]))

        @test isscalarexpr(:a)
        @test isscalarexpr(2.)
        @test isscalarexpr(:(2. + im*3))
        @test isscalarexpr(:(3*a+c))
        @test isscalarexpr(:(sin(x)+exp(-y)))
        @test isscalarexpr(:(scalar(a[x,y]*b[y,x])))
        @test !isscalarexpr(:(a[x,y]*b[y,x]))
        @test !isscalarexpr(:(3*scalar(a[x,y]*b[y,x]) + conj(c[z])))

        @test instantiate_eltype(:(a[1,2,3]*b[3,4,5]+c[1,2,4,5])) == :(promote_type(promote_type(eltype(a), eltype(b)), eltype(c)))
    end

    @testset "parsecost" begin
        using TensorOperations: parsecost, Power
        @test parsecost(:(3/5)) === 3/5
        @test parsecost(:(5+2)) === 5+2
        @test parsecost(:(Int128(2*8))) === Int128(2*8)
        @test parsecost(:(float(2*8))) === float(2*8)
        @test parsecost(:(big(2*8))) isa BigInt
        @test parsecost(:(x)) == Power{:x}(1,1)
        @test parsecost(:(3*x^3)) == Power{:x}(3,3)
        @test parsecost(:(3*x^3+2)) == Power{:x}(3,3)+2
        @test parsecost(:((3*x^3-2)/5)) == (Power{:x}(3,3)-2)/5
    end
end

@testset "methods for tensoropt" begin
    @testset "set methods" for T in (UInt32,UInt64,UInt128,BitSet,BitVector)
        maxint = 32
        A = Set(rand(1:maxint, 20))
        B = Set(rand(1:maxint, 20))
        U = union(A,B)
        I = intersect(A,B)
        S = setdiff(A,B)

        A2 = TensorOperations.storeset(T, A, maxint)
        B2 = TensorOperations.storeset(T, B, maxint)
        U2 = TensorOperations.storeset(T, U, maxint)
        I2 = TensorOperations.storeset(T, I, maxint)
        S2 = TensorOperations.storeset(T, S, maxint)

        @test TensorOperations._union(A2, B2) == U2
        @test TensorOperations._intersect(A2, B2) == I2
        @test TensorOperations._setdiff(A2, B2) == S2
        @test TensorOperations._isemptyset(TensorOperations.storeset(T, [], maxint))
        @test !TensorOperations._isemptyset(TensorOperations.storeset(T, [1], maxint))
    end

    @testset "poly" begin
        using TensorOperations: Power, degree, Poly
        x = Power{:x}(1,1)
        x⁰ = Power{:x}(1,0)
        x⁵= Power{:x}(1,5)
        @test x^5 == x⁵ == x*x*x*x*x
        @test x^0 == x⁰ == one(x)
        @test zero(x) == 0*x
        @test (5*x)/5 == x
        @test degree(x) == 1
        @test x + x^3 + x^5 isa Poly
        p = x - x^3 + x^5
        coeffs = [0,1,0,-1,0,1]
        @test p == Poly{:x}(coeffs)
        @test 2.5*p == Poly{:x}(2.5*coeffs) == p*2.5
        @test p/2 == Poly{:x}(coeffs/2) == 2\p
        @test p*p == x^10+x^6+x^2-2*x^8+2*x^6-2*x^4 == (x^10+x^6+x^2+1)-2*(x^8-x^6+x^4)-1
        p = x + x^3 + x^5
        for i = 0:6
            @test p[i] == isodd(i)
        end
        @test p + p/2 == 1.5*p
        @test p*one(p) == p == p*one(x)
        @test p*zero(p) == 0 == p*zero(x)
        @test p <= p*x
        @test x*p == p*x
    end

    @testset "cost methods" begin

    end
end

ext: jl


text:
@testset "bound" begin

    for platform in testplatforms

        #By local update
        m = @agent(
            3,

            [sx,sy,sz]::Local, #Steps

            [s,sMin,sMax]::Local, #Stop bounds
            [b,bMin,bMax]::Local, #Bounce bounds
            [r,rMin,rMax]::Local, #Hard bounds

            UpdateLocal = begin
                x += sx
                y += sy
                z += sz
            end
        )
        space = SimulationFree(m,
            box=[
                Bound(:x,-1,1,
                    stop=[:s],stopMin=[:sMin],stopMax=[:sMax],
                    bounce=[:b],bounceMin=[:bMin],bounceMax=[:bMax],
                    reflect=[:r],reflectMin=[:rMin],reflectMax=[:rMax]
                ),
                Bound(:y,-1,1,
                    stop=[:s],stopMin=[:sMin],stopMax=[:sMax],
                    bounce=[:b],bounceMin=[:bMin],bounceMax=[:bMax],
                    reflect=[:r],reflectMin=[:rMin],reflectMax=[:rMax]
                ),
                Bound(:z,-1,1,
                    stop=[:s],stopMin=[:sMin],stopMax=[:sMax],
                    bounce=[:b],bounceMin=[:bMin],bounceMax=[:bMax],
                    reflect=[:r],reflectMin=[:rMin],reflectMax=[:rMax]
                )
            ]
        )
        model = compile(m,space,platform=platform)
        #println(prettify(model.program))
        com = Community(model,N=6)
        com.x .= [-.9,.9,0,0,0,0]
        com.y .= [0,0,-.9,.9,0,0]
        com.z .= [0,0,0,0,-.9,.9]

        com.sx .= [-.2,.2,0,0,0,0]
        com.sy .= [0,0,-.2,.2,0,0]
        com.sz .= [0,0,0,0,-.2,.2]

        com.s .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.sMin .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.sMax .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]

        com.b .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.bMin .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.bMax .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]

        com.r .= 1
        com.rMin .= 1
        com.rMax .= 1

        comt = model.evolve(com,dt=1,tMax=1.1)

        @test all(abs.(comt[end].s .- [-1.,1.,-1.,1.,-1.,1.]) .< 0.0001)
        @test all(abs.(comt[end].sMin .- [-1.,1.2,-1.,1.2,-1.,1.2]) .< 0.0001)
        @test all(abs.(comt[end].sMax .- [-1.2,1.,-1.2,1.,-1.2,1.]) .< 0.0001)

        @test all(abs.(comt[end].b .- [-.8,.8,-.8,.8,-.8,.8]) .< 0.0001)
        @test all(abs.(comt[end].bMin .- [-.8,1.2,-.8,1.2,-.8,1.2]) .< 0.0001)
        @test all(abs.(comt[end].bMax .- [-1.2,.8,-1.2,.8,-1.2,.8]) .< 0.0001)

        @test all(abs.(comt[end].r .- [-1,-1,-1,-1,-1,-1]) .< 0.0001)
        @test all(abs.(comt[end].rMin .- [-1,1,-1,1,-1,1]) .< 0.0001)
        @test all(abs.(comt[end].rMax .- [1,-1,1,-1,1,-1]) .< 0.0001)

        m = @agent(
            3,

            [x1,y1,z1]::Local,
            [sx,sy,sz]::Local, #Steps

            UpdateLocal = begin
                x += sx
                y += sy
                z += sz
            end
        )
        space = SimulationFree(m,
            box=[
                Periodic(:x,-1,1,additional=[:x1]),
                Periodic(:y,-1,1,additional=[:y1]),
                Periodic(:z,-1,1,additional=[:z1])
            ]
        )
        model = compile(m,space,platform=platform)
        #println(prettify(model.program))

        com = Community(model,N=6)
        com.x .= [-.9,.9,0,0,0,0]
        com.y .= [0,0,-.9,.9,0,0]
        com.z .= [0,0,0,0,-.9,.9]

        com.sx .= [-.2,.2,0,0,0,0]
        com.sy .= [0,0,-.2,.2,0,0]
        com.sz .= [0,0,0,0,-.2,.2]

        com.x1 .= [-.9,.9,0,0,0,0]
        com.y1 .= [0,0,-.9,.9,0,0]
        com.z1 .= [0,0,0,0,-.9,.9]

        comt = model.evolve(com,dt=1,tMax=1.1)

        @test all(abs.(comt[end].x .- [.9,-.9,0,0,0,0]) .< 0.0001)
        @test all(abs.(comt[end].y .- [0,0,.9,-.9,0,0]) .< 0.0001)
        @test all(abs.(comt[end].z .- [0,0,0,0,.9,-.9]) .< 0.0001)

        @test all(abs.(comt[end].x1 .- [1.1,-1.1,0,0,0,0]) .< 0.0001)
        @test all(abs.(comt[end].y1 .- [0,0,1.1,-1.1,0,0]) .< 0.0001)
        @test all(abs.(comt[end].z1 .- [0,0,0,0,1.1,-1.1]) .< 0.0001)

        #Integration
        m = @agent(
            3,

            [sx,sy,sz]::Local, #Steps

            [s,sMin,sMax]::Local, #Stop bounds
            [b,bMin,bMax]::Local, #Bounce bounds
            [r,rMin,rMax]::Local, #Hard bounds

            Equation = begin
                d_x = sx*dt
                d_y = sy*dt
                d_z = sz*dt
            end
        )
        space = SimulationFree(m,
            box=[
                Bound(:x,-1,1,
                    stop=[:s],stopMin=[:sMin],stopMax=[:sMax],
                    bounce=[:b],bounceMin=[:bMin],bounceMax=[:bMax],
                    reflect=[:r],reflectMin=[:rMin],reflectMax=[:rMax]
                ),
                Bound(:y,-1,1,
                    stop=[:s],stopMin=[:sMin],stopMax=[:sMax],
                    bounce=[:b],bounceMin=[:bMin],bounceMax=[:bMax],
                    reflect=[:r],reflectMin=[:rMin],reflectMax=[:rMax]
                ),
                Bound(:z,-1,1,
                    stop=[:s],stopMin=[:sMin],stopMax=[:sMax],
                    bounce=[:b],bounceMin=[:bMin],bounceMax=[:bMax],
                    reflect=[:r],reflectMin=[:rMin],reflectMax=[:rMax]
                )
            ]
        )
        model = compile(m,space,platform=platform)
        #println(prettify(model.program))

        com = Community(model,N=6)
        com.x .= [-.9,.9,0,0,0,0]
        com.y .= [0,0,-.9,.9,0,0]
        com.z .= [0,0,0,0,-.9,.9]

        com.sx .= [-.2,.2,0,0,0,0]
        com.sy .= [0,0,-.2,.2,0,0]
        com.sz .= [0,0,0,0,-.2,.2]

        com.s .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.sMin .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.sMax .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]

        com.b .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.bMin .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]
        com.bMax .= [-1.2,1.2,-1.2,1.2,-1.2,1.2]

        com.r .= 1
        com.rMin .= 1
        com.rMax .= 1

        comt = model.evolve(com,dt=1,tMax=1.1)

        @test all(abs.(comt[end].s .- [-1.,1.,-1.,1.,-1.,1.]) .< 0.0001)
        @test all(abs.(comt[end].sMin .- [-1.,1.2,-1.,1.2,-1.,1.2]) .< 0.0001)
        @test all(abs.(comt[end].sMax .- [-1.2,1.,-1.2,1.,-1.2,1.]) .< 0.0001)

        @test all(abs.(comt[end].b .- [-.8,.8,-.8,.8,-.8,.8]) .< 0.0001)
        @test all(abs.(comt[end].bMin .- [-.8,1.2,-.8,1.2,-.8,1.2]) .< 0.0001)
        @test all(abs.(comt[end].bMax .- [-1.2,.8,-1.2,.8,-1.2,.8]) .< 0.0001)

        @test all(abs.(comt[end].r .- [-1,-1,-1,-1,-1,-1]) .< 0.0001)
        @test all(abs.(comt[end].rMin .- [-1,1,-1,1,-1,1]) .< 0.0001)
        @test all(abs.(comt[end].rMax .- [1,-1,1,-1,1,-1]) .< 0.0001)

        m = @agent(
            3,

            [x1,y1,z1]::Local,
            [sx,sy,sz]::Local, #Steps

            Equation = begin
                d_x = sx*dt
                d_y = sy*dt
                d_z = sz*dt
            end
        )
        space = SimulationFree(m,
            box=[
                Periodic(:x,-1,1,additional=[:x1]),
                Periodic(:y,-1,1,additional=[:y1]),
                Periodic(:z,-1,1,additional=[:z1])
            ]
        )
        model = compile(m,space,platform=platform)
        #println(prettify(model.program))

        com = Community(model,N=6)
        com.x .= [-.9,.9,0,0,0,0]
        com.y .= [0,0,-.9,.9,0,0]
        com.z .= [0,0,0,0,-.9,.9]

        com.sx .= [-.2,.2,0,0,0,0]
        com.sy .= [0,0,-.2,.2,0,0]
        com.sz .= [0,0,0,0,-.2,.2]

        com.x1 .= [-.9,.9,0,0,0,0]
        com.y1 .= [0,0,-.9,.9,0,0]
        com.z1 .= [0,0,0,0,-.9,.9]

        comt = model.evolve(com,dt=1,tMax=1.1)

        @test all(abs.(comt[end].x .- [.9,-.9,0,0,0,0]) .< 0.0001)
        @test all(abs.(comt[end].y .- [0,0,.9,-.9,0,0]) .< 0.0001)
        @test all(abs.(comt[end].z .- [0,0,0,0,.9,-.9]) .< 0.0001)

        @test all(abs.(comt[end].x1 .- [1.1,-1.1,0,0,0,0]) .< 0.0001)
        @test all(abs.(comt[end].y1 .- [0,0,1.1,-1.1,0,0]) .< 0.0001)
        @test all(abs.(comt[end].z1 .- [0,0,0,0,1.1,-1.1]) .< 0.0001)

    end

end

max_stars_repo_name: gatocor/embryogenesisJulia

max_issues_repo_name: gatocor/embryogenesisJulia

max_forks_repo_name: gatocor/embryogenesisJulia


text:
module JuliaCLI

import JSONRPC
using Base: Process
using Parameters: @unpack
using Pkg: TOML

include("rpcutils.jl")
include("utils.jl")
include("workers.jl")
include("broker.jl")
include("scripts.jl")

end # module


text:
export SimulatedMotorParams, SimulatedMotor

Base.@kwdef struct SimulatedMotorParams <: DeviceParams

end
SimulatedMotorParams(dict::Dict) = params_from_dict(SimulatedMotorParams, dict)

Base.@kwdef mutable struct SimulatedMotor <: Motor
  "Unique device ID for this device as defined in the configuration."
  deviceID::String
  "Parameter struct for this devices read from the configuration."
  params::SimulatedGaussMeterParams
  "Flag if the device is optional."
	optional::Bool = false
  "Flag if the device is present."
  present::Bool = false
  "Vector of dependencies for this device."
  dependencies::Dict{String, Union{Device, Missing}}

  direction::MotorDirection
  speed::typeof(1.0u"1/s")
end

function init(motor::SimulatedMotor)
  @debug "Initializing simulated motor unit with ID `$(motor.deviceID)`."

  motor.present = true
end

checkDependencies(motor::SimulatedMotor) = true

Base.close(motor::SimulatedMotor) = nothing


text:
export dualize

# MOI dualize
function dualize(
    primal_model::MOI.ModelLike;
    dual_names::DualNames = EMPTY_DUAL_NAMES,
    variable_parameters::Vector{VI} = VI[],
    ignore_objective::Bool = false,
    consider_constrained_variables::Bool = true,
)
    # Creates an empty dual problem
    dual_problem = DualProblem{Float64}()
    return dualize(
        primal_model,
        dual_problem,
        dual_names,
        variable_parameters,
        ignore_objective,
        consider_constrained_variables,
    )
end

function dualize(
    primal_model::MOI.ModelLike,
    dual_problem::DualProblem{T};
    dual_names::DualNames = EMPTY_DUAL_NAMES,
    variable_parameters::Vector{VI} = VI[],
    ignore_objective::Bool = false,
    consider_constrained_variables::Bool = true,
) where {T}
    # Dualize with the optimizer already attached
    return dualize(
        primal_model,
        dual_problem,
        dual_names,
        variable_parameters,
        ignore_objective,
        consider_constrained_variables,
    )
end

function dualize(
    primal_model::MOI.ModelLike,
    dual_problem::DualProblem{T},
    dual_names::DualNames,
    variable_parameters::Vector{VI},
    ignore_objective::Bool,
    consider_constrained_variables::Bool,
) where {T}
    # Throws an error if objective function cannot be dualized
    supported_objective(primal_model)

    # Query all constraint types of the model
    con_types = MOI.get(primal_model, MOI.ListOfConstraintTypesPresent())
    supported_constraints(con_types) # Errors if constraint cant be dualized

    # Set the dual model objective sense
    set_dual_model_sense(dual_problem.dual_model, primal_model)

    # Get primal objective in quadratic form
    # terms already split considering parameters
    primal_objective =
        get_primal_objective(primal_model, variable_parameters, T)

    # Cache information of which primal variables are `constrained_variables`
    # creating a map: constrained_var_idx, from original primal vars to original
    # constrains and their internal index (if vector constrains), 1 otherwise.
    # Also, initializes the map: `constrained_var_dual`, from original primal ci
    # to the dual constraint (latter is initilized as empty at this point).
    # If the Set constant of a VI-in-Set constraint is non-zero, the respective
    # primal variable will not be a constrained variable (with respect to that
    # constraint).
    if consider_constrained_variables
        add_constrained_variables(
            dual_problem,
            primal_model,
            variable_parameters,
        )
    end

    # Add variables to the dual model and their dual cone constraint.
    # Return a dictionary from dual variables to primal constraints
    # constants (obj coef of dual var)
    # Loops through all constraints that are not the constraint of a
    # constrained variable (defined by `add_constrained_variables`).
    # * creates the dual variable associated with the primal constraint
    # * fills `dual_obj_affine_terms`, since we are already looping through
    #   all constraints that might have constants.
    # * fills `primal_con_dual_var` mapping the primal constraint and the dual
    #   variable
    # * fills `primal_con_dual_con` to map the primal constraint to a
    #   constraint in the dual variable (if there is such constraint the dual
    #   dual variable is said to be constrained). If the primal constraint's set
    #   is EqualTo or Zeros, no constraint is added in the dual variable (the
    #   dual variable is said to be free).
    # * fills `primal_con_constants` mapping primal constraints to their
    #   respective constants, which might be inside the set.
    #   this map is used in `MOI.get(::DualOptimizer,::MOI.ConstraintPrimal,ci)`
    #   that requires extra information in the case that the scalar set constrains
    #   a constant (EqualtTo, GreaterThan, LessThan)
    dual_obj_affine_terms = add_dual_vars_in_dual_cones(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        con_types,
    )

    # Creates variables in the dual problem that represent parameters in the
    # primal model.
    # Fills `primal_parameter` mapping parameters in the primal to parameters
    # in the dual model.
    add_primal_parameter_vars(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        variable_parameters,
        primal_objective,
        ignore_objective,
    )

    # Add dual slack variables that are associated to the primal quadratic terms
    # All primal variables that appear in the objective products will have an
    # associated dual slack variable that is created here.
    # also, `primal_var_dual_quad_slack` is filled, mapping primal variables
    # (that appear in quadritc objective terms) to dual "slack" variables.
    add_quadratic_slack_vars(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        primal_objective,
    )

    # Add dual constraints
    # that will be equality if associated to "free variables"
    # but will be constrained in the dual set of the associated primal
    # constrained variable if such variable is not "free"
    # Also, fills the link dictionary.
    # returns `scalar_affine_terms`
    # because the terms associated to variables that are parameters will be used
    # in `get_dual_objective`
    scalar_affine_terms = add_dual_equality_constraints(
        dual_problem.dual_model,
        primal_model,
        dual_problem.primal_dual_map,
        dual_names,
        primal_objective,
        con_types,
        variable_parameters,
    )

    if ignore_objective
        # do not add objective
    else
        # Fill Dual Objective Coefficients Struct
        dual_objective = get_dual_objective(
            dual_problem,
            dual_obj_affine_terms,
            primal_objective,
            scalar_affine_terms,
            variable_parameters,
        )
        # Add dual objective to the model
        set_dual_objective(dual_problem.dual_model, dual_objective)
    end
    return dual_problem
end

# JuMP dualize
function dualize(model::JuMP.Model; dual_names::DualNames = EMPTY_DUAL_NAMES)
    # Create an empty JuMP model
    JuMP_model = JuMP.Model()

    if JuMP.mode(model) != JuMP.AUTOMATIC # Only works in AUTOMATIC mode
        error(
            "Dualization does not support solvers in $(model.moi_backend.mode) mode",
        )
    end
    # Dualize and attach to the model
    dualize(
        backend(model),
        DualProblem(backend(JuMP_model));
        dual_names = dual_names,
    )
    fill_obj_dict_with_variables!(JuMP_model)
    fill_obj_dict_with_constraints!(JuMP_model)
    return JuMP_model
end

function dualize(
    model::JuMP.Model,
    optimizer_constructor;
    dual_names::DualNames = EMPTY_DUAL_NAMES,
)
    # Dualize the JuMP model
    dual_JuMP_model = dualize(model; dual_names = dual_names)
    # Set the optimizer
    JuMP.set_optimizer(dual_JuMP_model, optimizer_constructor)
    return dual_JuMP_model
end

# dualize docs
"""
    dualize(args...; kwargs...)

The `dualize` function works in three different ways. The user can provide:

* A `MathOptInterface.ModelLike`

The function will return a `DualProblem` struct that has the dualized model
and `PrimalDualMap{Float64}` for users to identify the links between primal
and dual model. The `PrimalDualMap{Float64}` maps variables and constraints
from the original primal model into the respective objects of the dual model.

* A `MathOptInterface.ModelLike` and a `DualProblem{T}`

* A `JuMP.Model`

The function will return a JuMP model with the dual representation of the problem.

* A `JuMP.Model` and an optimizer constructor

The function will return a JuMP model with the dual representation of the problem with
the optimizer constructor attached.

On each of these methods, the user can provide the following keyword arguments:

* `dual_names`: of type `DualNames` struct. It allows users to set more intuitive names
for the dual variables and dual constraints created.

* `variable_parameters`: A vector of MOI.VariableIndex containing the variables that
should not be considered model variables during dualization. These variables will behave
like constants during dualization. This is specially useful for the case of bi-level
modelling, where the second level depends on some decisions from the upper level.

* `ignore_objective`: a boolean indicating if the objective function should be
added to the dual model. This is also useful for bi-level modelling, where the second
level model is represented as a KKT in the upper level model.

"""
function dualize end

function fill_obj_dict_with_variables!(model::JuMP.Model)
    list = MOI.get(model, MOI.ListOfVariableAttributesSet())
    if !(MOI.VariableName() in list)
        return
    end
    all_indices =
        MOI.get(model, MOI.ListOfVariableIndices())::Vector{MOI.VariableIndex}
    for vi in all_indices
        name = MOI.get(backend(model), MOI.VariableName(), vi)
        if !isempty(name)
            model[Symbol(name)] = VariableRef(model, vi)
        end
    end
    return
end

function fill_obj_dict_with_constraints!(model::JuMP.Model)
    con_types = MOI.get(model, JuMP.MOI.ListOfConstraintTypesPresent())
    for (F, S) in con_types
        fill_obj_dict_with_constraints!(model, F, S)
    end
    return
end

function fill_obj_dict_with_constraints!(model::JuMP.Model, F::Type, S::Type)
    list = MOI.get(model, MOI.ListOfConstraintAttributesSet{F,S}())
    if !(MOI.ConstraintName() in list)
        return
    end
    for ci in MOI.get(backend(model), MOI.ListOfConstraintIndices{F,S}())
        name = MOI.get(backend(model), MOI.ConstraintName(), ci)
        if !isempty(name)
            con_ref = JuMP.constraint_ref_with_index(model, ci)
            model[Symbol(name)] = con_ref
        end
    end
end


text:
################################################################################
#
# Generate synthetic data for the Gaussian lens model
#
################################################################################

### Module loading

using LinearAlgebra, JUDI.TimeModeling, JLD2, PyPlot
push!(LOAD_PATH, string(pwd(), "/src/"))
using TWRIdual

### Set true model

n = (201, 201)
d = (10f0, 10f0)
o = (0f0, 0f0)
val = 0.5f0
x = reshape((collect(1:n[1]).-1f0).*d[1], :, 1)
z = reshape((collect(1:n[2]).-1f0).*d[2], 1, :)
v = 2f0.-val*exp.(-((x./1000f0.-1f0)/0.25f0).^2f0.-((z./1000f0.-1f0)/0.5f0).^2f0)
m_true = 1f0./v.^2f0
model_true = Model(n, d, o, m_true)

### Acquisition geometry

# Sources
ix_src = 3:14:199
nsrc = length(ix_src)
iz_src = 3
x_src = convertToCell((ix_src.-1)*d[1])
y_src = convertToCell(range(0f0, stop = 0f0, length = nsrc))
z_src = convertToCell(range((iz_src-1)*d[2], stop = (iz_src-1)*d[2], length = nsrc))

# Source wavelet
dt = 1f0
T = 3000f0 # total recording time [ms]
freq_peak = 0.006f0
wavelet = ricker_wavelet(T, dt, freq_peak)

# Receivers
ix_rcv = 3:199
iz_rcv = 199
nrcv = length(ix_rcv)
x_rcv = (ix_rcv.-1)*d[1]
y_rcv = 0f0
z_rcv = range((iz_rcv-1)*d[2], stop = (iz_rcv-1)*d[2], length = nrcv)

# Geometry structures
src_geom = Geometry(x_src, y_src, z_src; dt = dt, t = T)
rcv_geom = Geometry(x_rcv, y_rcv, z_rcv; dt = dt, t = T, nsrc = nsrc)

# Source function
fsrc = judiVector(src_geom, wavelet)

### Modeling synthetic data

dt_comp = Array{Any, 1}(undef, nsrc); dt_comp .= dt
@time dat = gendata(model_true, fsrc, rcv_geom; opt = Options(), dt_comp = dt_comp)

### Saving data

@save string("./data/GaussLens/GaussLens_data.jld") model_true fsrc dat


text:
module UUID

VERSION >= v"0.4-" && warn("UUID.jl is deprecated, please use Base.Random.uuid1(), Base.Random.uuid4(), and Base.Random.UUID instead.")

using Compat

import Base.show, Base.print, Base.string, Base.int, Base.uint, Base.hex

# uuid layout and byteorder from RFC 4122
# 0                   1                   2                   3
# 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1 2 3 4 5 6 7 8 9 0 1
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# |                          time_low                             |
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# |       time_mid                |         time_hi_and_version   |
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# |clk_seq_hi_res |  clk_seq_low  |         node (0-1)            |
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+
# |                         node (2-5)                            |
# +-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+

type TimestampUuid
    ts_low::Uint32
    ts_mid::Uint16
    ts_hi_and_ver::Uint16
end

type ClockSeqUuid
    clock_seq_hi_and_rsvd::Uint8
    clock_seq_low::Uint8
end

# regex pattern to parse uuid
const regex_uuid = r"^(urn\:uuid\:)?\{?([a-z0-9]{8})-([a-z0-9]{4})-([1-5][a-z0-9]{3})-([a-z0-9]{4})-([a-z0-9]{12})\}?$"

immutable Uuid
    ts::TimestampUuid
    cseq::ClockSeqUuid
    node::Vector{Uint8}

    Uuid(ts, cseq, node) =
         length(node) != 6 ? error("Node vector size should be 6") :
                             new(TimestampUuid(ts.ts_low,ts.ts_mid,ts.ts_hi_and_ver),
                                 ClockSeqUuid(cseq.clock_seq_hi_and_rsvd,
                                              cseq.clock_seq_low),
                                 node)
    function Uuid(uuid::String)
        uuid_match = match(regex_uuid,uuid)
        if uuid_match == nothing
            error("Not a valid UUID string")
        end
        new(TimestampUuid(uint32(parseint(uuid_match.captures[2],16)),
                          uint16(parseint(uuid_match.captures[3],16)),
                          uint16(parseint(uuid_match.captures[4],16))),
            ClockSeqUuid(uint8(parseint(uuid_match.captures[5][1:2],16)),
                         uint8(parseint(uuid_match.captures[5][3:end],16))),
            [uint8(parseint(uuid_match.captures[6][1:2],16)),
             uint8(parseint(uuid_match.captures[6][3:4],16)),
             uint8(parseint(uuid_match.captures[6][5:6],16)),
             uint8(parseint(uuid_match.captures[6][7:8],16)),
             uint8(parseint(uuid_match.captures[6][9:10],16)),
             uint8(parseint(uuid_match.captures[6][11:12],16))])
    end
end

function build_timestamp()
    timestamp = get_timestamp()
    ts_low::Uint32 = timestamp & typemax(Uint32)
    ts_mid ::Uint16 = (timestamp >> 32) & typemax(Uint16)
    # timestamp is 60 bits so cutoff upper nibble
    ts_hi_and_ver::Uint16 = ((timestamp >> 48) & 0x0FFF)
    TimestampUuid(ts_low,ts_mid,ts_hi_and_ver)
end

function build_version(ts,ver)
    ts.ts_hi_and_ver |= (ver << 12)
    ts
end

function build_clock_sequence(clock_seq)
    if clock_seq == 0
        # generate a random clock sequence
        clock_seq = rand(1:((1<<14) - 1))
    end
    clock_seq_low::Uint8 = clock_seq & 0xFF
    clock_seq_hi_and_rsvd::Uint8 = ((clock_seq & 0x3F00) >> 8) | 0x80
    ClockSeqUuid(clock_seq_hi_and_rsvd,clock_seq_low)
end

# as specified in the RFC 4.5 we can generate a random node id
# and set unicast/multicast bit to 1 (the least significant first octet
# bit)
function build_node()
    node = get_rand_narray(6)
    # set unicast/multicast bit
    node[2] |= 0x01
    node
end

function get_rand_narray(n)
    rand(Uint8,n)
end

# Per RFC 4122:
# Get the current time as a 60-bit count of 100-nanosecond intervals
# since 00:00:00.00, 15 October 1582
# 0x01b21dd213814000 is the number of 100-ns intervals between the
# UUID epoch 1582-10-15 00:00:00 and the Unix epoch 1970-01-01 00:00:00
const number_100ns_intervals = 0x01b21dd213814000

function get_timestamp()
    # 1e9 /100 for 100ns intervals
    round(Integer, time() * 1e7) + number_100ns_intervals
end

function v1(; clock_seq = 0)
    ts = build_timestamp()
    ts = build_version(ts,1)
    clock_seq = build_clock_sequence(clock_seq)
    node = build_node()
    Uuid(ts, clock_seq, node)
end

# version 4 is a random uuid
function v4()
    ts = TimestampUuid(rand(Uint32),
                       rand(Uint16),
                       rand(Uint16) & 0x0FFF) # clear version
    build_version(ts,4)
    Uuid(ts,
         ClockSeqUuid(rand(Uint8), rand(Uint8)),
         get_rand_narray(6))
end

function asString(uuid::Uuid, joiner::String="")
    string(hex(uuid.ts.ts_low, 8), joiner,
           hex(uuid.ts.ts_mid, 4), joiner,
           hex(uuid.ts.ts_hi_and_ver, 4), joiner,
           hex(uuid.cseq.clock_seq_hi_and_rsvd, 2),
           hex(uuid.cseq.clock_seq_low, 2),
           joiner,
           hex(uuid.node[1],2), hex(uuid.node[2],2),
           hex(uuid.node[3],2), hex(uuid.node[4],2),
           hex(uuid.node[5],2), hex(uuid.node[6],2))
end

hex(uuid::Uuid) = asString(uuid)
string(uuid::Uuid) = asString(uuid,"-")

function uint(uuid::Uuid)
    (convert(Uint128,uuid.ts.ts_low)                  << 96) |
    (convert(Uint128,uuid.ts.ts_mid)                  << 80) |
    (convert(Uint128,uuid.ts.ts_hi_and_ver)           << 64) |
    (convert(Uint128,uuid.cseq.clock_seq_hi_and_rsvd) << 56) |
    (convert(Uint128,uuid.cseq.clock_seq_low)         << 48) |
    (convert(Uint128,uuid.node[1])                    << 40) |
    (convert(Uint128,uuid.node[2])                    << 32) |
    (convert(Uint128,uuid.node[3])                    << 24) |
    (convert(Uint128,uuid.node[4])                    << 16) |
    (convert(Uint128,uuid.node[5])                    <<  8) |
    (convert(Uint128,uuid.node[6])                         )
end

int(uuid::Uuid) = reinterpret(Int128, uint(uuid))

function get_version(uuid::Uuid)
    uuid.ts.ts_hi_and_ver & 0xF000 >> 12 |> int
end

# override default show and print for UUID type
show(io::IO, uuid::Uuid) = print(io,string("Uuid(\'",uuid,"\')"))

print(io::IO, uuid::Uuid) = print(io, string(uuid))

end


text:
# This file is a part of TypeDBClient.  License is MIT: https://github.com/Humans-of-Julia/TypeDBClient.jl/blob/main/LICENSE

#
# package typedb.client.api.concept.thing;
#
# import typedb.client.api.TypeDBTransaction;
# import typedb.client.api.concept.type.AttributeType;
# import typedb.client.api.concept.type.ThingType;
#
# import javax.annotation.CheckReturnValue;
# import java.time.LocalDateTime;
# import java.util.stream.Stream;
#
# public interface Attribute<VALUE> extends Thing {
#
#     @Override
#     @CheckReturnValue
#     AttributeType getType();
#
#     @CheckReturnValue
#     VALUE getValue();
#
#     @Override
#     @CheckReturnValue
#     default boolean isAttribute() {
#         return true;
#     }
#
#     @CheckReturnValue
#     default boolean isBoolean() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isLong() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isDouble() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isString() {
#         return false;
#     }
#
#     @CheckReturnValue
#     default boolean isDateTime() {
#         return false;
#     }
#
#     @CheckReturnValue
#     Attribute.Boolean asBoolean();
#
#     @CheckReturnValue
#     Attribute.Long asLong();
#
#     @CheckReturnValue
#     Attribute.Double asDouble();
#
#     @CheckReturnValue
#     Attribute.String asString();
#
#     @CheckReturnValue
#     Attribute.DateTime asDateTime();
#
#     @Override
#     @CheckReturnValue
#     Attribute.Remote<VALUE> asRemote(TypeDBTransaction transaction);
#
#     interface Remote<VALUE> extends Thing.Remote, Attribute<VALUE> {
#
#         @CheckReturnValue
#         Stream<? extends Thing> getOwners();
#
#         @CheckReturnValue
#         Stream<? extends Thing> getOwners(ThingType ownerType);
#
#         @Override
#         @CheckReturnValue
#         Attribute.Remote<VALUE> asAttribute();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Boolean.Remote asBoolean();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Long.Remote asLong();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Double.Remote asDouble();
#
#         @Override
#         @CheckReturnValue
#         Attribute.String.Remote asString();
#
#         @Override
#         @CheckReturnValue
#         Attribute.DateTime.Remote asDateTime();
#     }
#
#     interface Boolean extends Attribute<java.lang.Boolean> {
#
#         @Override
#         @CheckReturnValue
#         default boolean isBoolean() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Boolean getType();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Boolean.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends Attribute.Boolean, Attribute.Remote<java.lang.Boolean> {}
#     }
#
#     interface Long extends Attribute<java.lang.Long> {
#
#         @Override
#         @CheckReturnValue
#         default boolean isLong() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Long getType();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Long.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends Attribute.Long, Attribute.Remote<java.lang.Long> {}
#     }
#
#     interface Double extends Attribute<java.lang.Double> {
#
#         @Override
#         @CheckReturnValue
#         default boolean isDouble() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.Double getType();
#
#         @Override
#         @CheckReturnValue
#         Attribute.Double.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends Attribute.Double, Attribute.Remote<java.lang.Double> {}
#     }
#
#     interface String extends Attribute<java.lang.String> {
#
#         @Override
#         @CheckReturnValue
#         default boolean isString() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.String getType();
#
#         @Override
#         @CheckReturnValue
#         Attribute.String.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends Attribute.String, Attribute.Remote<java.lang.String> {}
#     }
#
#     interface DateTime extends Attribute<LocalDateTime> {
#
#         @Override
#         @CheckReturnValue
#         default boolean isDateTime() {
#             return true;
#         }
#
#         @Override
#         @CheckReturnValue
#         AttributeType.DateTime getType();
#
#         @Override
#         @CheckReturnValue
#         Attribute.DateTime.Remote asRemote(TypeDBTransaction transaction);
#
#         interface Remote extends Attribute.DateTime, Attribute.Remote<LocalDateTime> {}
#     }
# }


text:
using DataFrames, LowRankModels, Compat

# boolean example with only entries greater than threshold t observed
# ie, censored data
# example with only entries greater than threshold t observed
m,n,k,ktrue = 100,100,1,1
A = rand(m,ktrue)*rand(ktrue,n)
println("max value of A is ",maximum(maximum(A))," which is less than $ktrue")
B = round(Int, ktrue*rand(m,n) .>= A) # Bernoulli samples with probability proportional to A
losses = fill(QuadLoss(),n)
r = QuadReg(.1)
obs = @compat Array(Tuple{Int,Int}, 0)
for i=1:m
    for j=1:n
        if B[i,j] == 1
            push!(obs, (i,j))
        end
    end
end

(train_observed_features, train_observed_examples, test_observed_features,  test_observed_examples) =
    get_train_and_test(obs, m, n, .2)

train_glrm = GLRM(B,losses,r,r,k, observed_features=train_observed_features, observed_examples=train_observed_examples)

train_error, test_error, prec_at_k, train_time, reg_params, solution =
    precision_at_k(train_glrm, test_observed_features, params=Params(1,200,.00001,.01),
                                 reg_params=logspace(2,-2,9))
df = DataFrame(train_error = train_error, prec_at_k = prec_at_k,
                   train_time = train_time, reg_param = reg_params, solution_1norm = [s[2] for s in solution])
println(df)


text:
# Command: Execute Script
"""
    script!(session::Session, script::AbstractString, args...; async::Bool = false)

Executes JavaScript (sync / async)

!!! note

    This command uses the JSON Wire Protocol instead of the current W3c WebDriver API.
"""
function script!(session::Session, _script::AbstractString, args...; async::Bool = false)
    @unpack addr, id = session
    response = HTTP.post(
        "$addr/session/$id/execute$(async ? "_async" : "")",
        [("Content-Type" => "application/json; charset=utf-8")],
        JSON3.write(Dict(
            "script" => _script,
            "args" => [isa(arg, Element) ? "ELEMENT" => arg.id : arg for arg in args],
        )),
    )
    @assert response.status == 200
    JSON3.read(response.body).value
end


text:
# Utilities for filtering Quantities with window operations

"""
Filters the quantity q with a window function w

w::AbstractArray{T,1} -> T
"""
function filter_window{T<:Quantity}(q::T,w,k)
    ts,qs = unzip(q)
    N = length(qs)
    ps = zeros(qs)
    kk = isodd(k) ? k : k+1
    k2 = kk>>1
    for i in eachindex(ps)
        i1 = max(1,i-k2)
        i2 = min(N,i+k2)
        ps[i] = w(qs[i1:i2])
    end
    T(ts,ps)
end


text:
"""
    AbstractScaling <: Transform

Linearly scale the data according to some statistics.
"""
abstract type AbstractScaling <: Transform end
cardinality(::AbstractScaling) = OneToOne()

"""
    IdentityScaling <: AbstractScaling

Represents the no-op scaling which simply returns the `data` it is applied on.
"""
struct IdentityScaling <: AbstractScaling end
IdentityScaling(args...) = IdentityScaling()

@inline _apply(x, ::IdentityScaling; kwargs...) = x

"""
    MeanStdScaling(μ, σ) <: AbstractScaling

Linearly scale the data by the statistical mean `μ` and standard deviation `σ`.
This is also known as standardization, or the Z score transform.

# Keyword arguments to `apply`
* `inverse=true`: inverts the scaling (e.g. to reconstruct the unscaled data).
* `eps=1e-3`: used in place of all 0 values in `σ` before scaling (if `inverse=false`).
"""
struct MeanStdScaling <: AbstractScaling
    μ::Real
    σ::Real

    """
        MeanStdScaling(A::AbstractArray; dims=:, inds=:) -> MeanStdScaling
        MeanStdScaling(table, [cols]) -> MeanStdScaling

    Construct a [`MeanStdScaling`](@ref) transform from the statistics of the given data.
    By default _all the data_ is considered when computing the mean and standard deviation.
    This can be restricted to certain slices via the keyword arguments (see below).

    Since `MeanStdScaling` is a stateful transform, i.e. the parameters depend on the data
    it's given, you should define it independently before applying it so you can keep the
    information for later use. For instance, if you want to invert the transform or apply it
    to a test set.

    # `AbstractArray` keyword arguments
    * `dims=:`: the dimension along which to take the `inds` slices. Default uses all dims.
    * `inds=:`: the indices to use in computing the statistics. Default uses all indices.

    # `Table` keyword arguments
    * `cols`: the columns to use in computing the statistics. Default uses all columns.

    !!! note
        If you want the `MeanStdScaling` to transform your data consistently you should use
        the same `inds`, `dims`, or `cols` keywords when calling `apply`. Otherwise, `apply`
        might rescale the wrong data or throw an error.
    """
    function MeanStdScaling(A::AbstractArray; dims=:, inds=:)
        dims == Colon() && return new(compute_stats(A)...)
        return new(compute_stats(selectdim(A, dims, inds))...)
    end

    function MeanStdScaling(table; cols=_get_cols(table))
        Tables.istable(table) || throw(MethodError(MeanStdScaling, table))
        columntable = Tables.columns(table)
        data = reduce(vcat, [getproperty(columntable, c) for c in _to_vec(cols)])
        return new(compute_stats(data)...)
    end
end

compute_stats(x) = (mean(x), std(x))

function _apply(A::AbstractArray, scaling::MeanStdScaling; inverse=false, eps=1e-3, kwargs...)
    inverse && return scaling.μ .+ scaling.σ .* A
    # Avoid division by 0
    # If std is 0 then data was uniform, so the scaled value would end up ≈ 0
    # Therefore the particular `eps` value should not matter much.
    σ_safe = max(scaling.σ, eps)
    return (A .- scaling.μ) ./ σ_safe
end


text:
#
# Autogenerated by Thrift Compiler (1.0.0-dev)
#
# DO NOT EDIT UNLESS YOU ARE SURE THAT YOU KNOW WHAT YOU ARE DOING

type _enum_TProtocolVersion
  HIVE_CLI_SERVICE_PROTOCOL_V1::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V2::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V3::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V4::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V5::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V6::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V7::Int32
  HIVE_CLI_SERVICE_PROTOCOL_V8::Int32
end
const TProtocolVersion = _enum_TProtocolVersion(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

type _enum_TTypeId
  BOOLEAN_TYPE::Int32
  TINYINT_TYPE::Int32
  SMALLINT_TYPE::Int32
  INT_TYPE::Int32
  BIGINT_TYPE::Int32
  FLOAT_TYPE::Int32
  DOUBLE_TYPE::Int32
  STRING_TYPE::Int32
  TIMESTAMP_TYPE::Int32
  BINARY_TYPE::Int32
  ARRAY_TYPE::Int32
  MAP_TYPE::Int32
  STRUCT_TYPE::Int32
  UNION_TYPE::Int32
  USER_DEFINED_TYPE::Int32
  DECIMAL_TYPE::Int32
  NULL_TYPE::Int32
  DATE_TYPE::Int32
  VARCHAR_TYPE::Int32
  CHAR_TYPE::Int32
  INTERVAL_YEAR_MONTH_TYPE::Int32
  INTERVAL_DAY_TIME_TYPE::Int32
end
const TTypeId = _enum_TTypeId(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7), Int32(8), Int32(9), Int32(10), Int32(11), Int32(12), Int32(13), Int32(14), Int32(15), Int32(16), Int32(17), Int32(18), Int32(19), Int32(20), Int32(21))

type _enum_TStatusCode
  SUCCESS_STATUS::Int32
  SUCCESS_WITH_INFO_STATUS::Int32
  STILL_EXECUTING_STATUS::Int32
  ERROR_STATUS::Int32
  INVALID_HANDLE_STATUS::Int32
end
const TStatusCode = _enum_TStatusCode(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4))

type _enum_TOperationState
  INITIALIZED_STATE::Int32
  RUNNING_STATE::Int32
  FINISHED_STATE::Int32
  CANCELED_STATE::Int32
  CLOSED_STATE::Int32
  ERROR_STATE::Int32
  UKNOWN_STATE::Int32
  PENDING_STATE::Int32
end
const TOperationState = _enum_TOperationState(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7))

type _enum_TOperationType
  EXECUTE_STATEMENT::Int32
  GET_TYPE_INFO::Int32
  GET_CATALOGS::Int32
  GET_SCHEMAS::Int32
  GET_TABLES::Int32
  GET_TABLE_TYPES::Int32
  GET_COLUMNS::Int32
  GET_FUNCTIONS::Int32
  UNKNOWN::Int32
end
const TOperationType = _enum_TOperationType(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5), Int32(6), Int32(7), Int32(8))

type _enum_TGetInfoType
  CLI_MAX_DRIVER_CONNECTIONS::Int32
  CLI_MAX_CONCURRENT_ACTIVITIES::Int32
  CLI_DATA_SOURCE_NAME::Int32
  CLI_FETCH_DIRECTION::Int32
  CLI_SERVER_NAME::Int32
  CLI_SEARCH_PATTERN_ESCAPE::Int32
  CLI_DBMS_NAME::Int32
  CLI_DBMS_VER::Int32
  CLI_ACCESSIBLE_TABLES::Int32
  CLI_ACCESSIBLE_PROCEDURES::Int32
  CLI_CURSOR_COMMIT_BEHAVIOR::Int32
  CLI_DATA_SOURCE_READ_ONLY::Int32
  CLI_DEFAULT_TXN_ISOLATION::Int32
  CLI_IDENTIFIER_CASE::Int32
  CLI_IDENTIFIER_QUOTE_CHAR::Int32
  CLI_MAX_COLUMN_NAME_LEN::Int32
  CLI_MAX_CURSOR_NAME_LEN::Int32
  CLI_MAX_SCHEMA_NAME_LEN::Int32
  CLI_MAX_CATALOG_NAME_LEN::Int32
  CLI_MAX_TABLE_NAME_LEN::Int32
  CLI_SCROLL_CONCURRENCY::Int32
  CLI_TXN_CAPABLE::Int32
  CLI_USER_NAME::Int32
  CLI_TXN_ISOLATION_OPTION::Int32
  CLI_INTEGRITY::Int32
  CLI_GETDATA_EXTENSIONS::Int32
  CLI_NULL_COLLATION::Int32
  CLI_ALTER_TABLE::Int32
  CLI_ORDER_BY_COLUMNS_IN_SELECT::Int32
  CLI_SPECIAL_CHARACTERS::Int32
  CLI_MAX_COLUMNS_IN_GROUP_BY::Int32
  CLI_MAX_COLUMNS_IN_INDEX::Int32
  CLI_MAX_COLUMNS_IN_ORDER_BY::Int32
  CLI_MAX_COLUMNS_IN_SELECT::Int32
  CLI_MAX_COLUMNS_IN_TABLE::Int32
  CLI_MAX_INDEX_SIZE::Int32
  CLI_MAX_ROW_SIZE::Int32
  CLI_MAX_STATEMENT_LEN::Int32
  CLI_MAX_TABLES_IN_SELECT::Int32
  CLI_MAX_USER_NAME_LEN::Int32
  CLI_OJ_CAPABILITIES::Int32
  CLI_XOPEN_CLI_YEAR::Int32
  CLI_CURSOR_SENSITIVITY::Int32
  CLI_DESCRIBE_PARAMETER::Int32
  CLI_CATALOG_NAME::Int32
  CLI_COLLATION_SEQ::Int32
  CLI_MAX_IDENTIFIER_LEN::Int32
end
const TGetInfoType = _enum_TGetInfoType(Int32(0), Int32(1), Int32(2), Int32(8), Int32(13), Int32(14), Int32(17), Int32(18), Int32(19), Int32(20), Int32(23), Int32(25), Int32(26), Int32(28), Int32(29), Int32(30), Int32(31), Int32(32), Int32(34), Int32(35), Int32(43), Int32(46), Int32(47), Int32(72), Int32(73), Int32(81), Int32(85), Int32(86), Int32(90), Int32(94), Int32(97), Int32(98), Int32(99), Int32(100), Int32(101), Int32(102), Int32(104), Int32(105), Int32(106), Int32(107), Int32(115), Int32(10000), Int32(10001), Int32(10002), Int32(10003), Int32(10004), Int32(10005))

type _enum_TFetchOrientation
  FETCH_NEXT::Int32
  FETCH_PRIOR::Int32
  FETCH_RELATIVE::Int32
  FETCH_ABSOLUTE::Int32
  FETCH_FIRST::Int32
  FETCH_LAST::Int32
end
const TFetchOrientation = _enum_TFetchOrientation(Int32(0), Int32(1), Int32(2), Int32(3), Int32(4), Int32(5))

typealias TTypeEntryPtr Int32

typealias TIdentifier Compat.UTF8String

typealias TPattern Compat.UTF8String

typealias TPatternOrIdentifier Compat.UTF8String

type TTypeQualifierValue
  i32Value::Int32
  stringValue::Compat.UTF8String
  TTypeQualifierValue() = (o=new(); fillunset(o); o)
end # type TTypeQualifierValue
meta(t::Type{TTypeQualifierValue}) = meta(t, Symbol[:i32Value,:stringValue], Int[], Dict{Symbol,Any}())

type TTypeQualifiers
  qualifiers::Dict{Compat.UTF8String,TTypeQualifierValue}
  TTypeQualifiers() = (o=new(); fillunset(o); o)
end # type TTypeQualifiers

type TPrimitiveTypeEntry
  _type::Int32
  typeQualifiers::TTypeQualifiers
  TPrimitiveTypeEntry() = (o=new(); fillunset(o); o)
end # type TPrimitiveTypeEntry
meta(t::Type{TPrimitiveTypeEntry}) = meta(t, Symbol[:typeQualifiers], Int[], Dict{Symbol,Any}())

type TArrayTypeEntry
  objectTypePtr::TTypeEntryPtr
  TArrayTypeEntry() = (o=new(); fillunset(o); o)
end # type TArrayTypeEntry

type TMapTypeEntry
  keyTypePtr::TTypeEntryPtr
  valueTypePtr::TTypeEntryPtr
  TMapTypeEntry() = (o=new(); fillunset(o); o)
end # type TMapTypeEntry

type TStructTypeEntry
  nameToTypePtr::Dict{Compat.UTF8String,TTypeEntryPtr}
  TStructTypeEntry() = (o=new(); fillunset(o); o)
end # type TStructTypeEntry

type TUnionTypeEntry
  nameToTypePtr::Dict{Compat.UTF8String,TTypeEntryPtr}
  TUnionTypeEntry() = (o=new(); fillunset(o); o)
end # type TUnionTypeEntry

type TUserDefinedTypeEntry
  typeClassName::Compat.UTF8String
  TUserDefinedTypeEntry() = (o=new(); fillunset(o); o)
end # type TUserDefinedTypeEntry

type TTypeEntry
  primitiveEntry::TPrimitiveTypeEntry
  arrayEntry::TArrayTypeEntry
  mapEntry::TMapTypeEntry
  structEntry::TStructTypeEntry
  unionEntry::TUnionTypeEntry
  userDefinedTypeEntry::TUserDefinedTypeEntry
  TTypeEntry() = (o=new(); fillunset(o); o)
end # type TTypeEntry

type TTypeDesc
  types::Vector{TTypeEntry}
  TTypeDesc() = (o=new(); fillunset(o); o)
end # type TTypeDesc

type TColumnDesc
  columnName::Compat.UTF8String
  typeDesc::TTypeDesc
  position::Int32
  comment::Compat.UTF8String
  TColumnDesc() = (o=new(); fillunset(o); o)
end # type TColumnDesc
meta(t::Type{TColumnDesc}) = meta(t, Symbol[:comment], Int[], Dict{Symbol,Any}())

type TTableSchema
  columns::Vector{TColumnDesc}
  TTableSchema() = (o=new(); fillunset(o); o)
end # type TTableSchema

type TBoolValue
  value::Bool
  TBoolValue() = (o=new(); fillunset(o); o)
end # type TBoolValue
meta(t::Type{TBoolValue}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TByteValue
  value::UInt8
  TByteValue() = (o=new(); fillunset(o); o)
end # type TByteValue
meta(t::Type{TByteValue}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TI16Value
  value::Int16
  TI16Value() = (o=new(); fillunset(o); o)
end # type TI16Value
meta(t::Type{TI16Value}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TI32Value
  value::Int32
  TI32Value() = (o=new(); fillunset(o); o)
end # type TI32Value
meta(t::Type{TI32Value}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TI64Value
  value::Int64
  TI64Value() = (o=new(); fillunset(o); o)
end # type TI64Value
meta(t::Type{TI64Value}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TDoubleValue
  value::Float64
  TDoubleValue() = (o=new(); fillunset(o); o)
end # type TDoubleValue
meta(t::Type{TDoubleValue}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TStringValue
  value::Compat.UTF8String
  TStringValue() = (o=new(); fillunset(o); o)
end # type TStringValue
meta(t::Type{TStringValue}) = meta(t, Symbol[:value], Int[], Dict{Symbol,Any}())

type TColumnValue
  boolVal::TBoolValue
  byteVal::TByteValue
  i16Val::TI16Value
  i32Val::TI32Value
  i64Val::TI64Value
  doubleVal::TDoubleValue
  stringVal::TStringValue
  TColumnValue() = (o=new(); fillunset(o); o)
end # type TColumnValue

type TRow
  colVals::Vector{TColumnValue}
  TRow() = (o=new(); fillunset(o); o)
end # type TRow

type TBoolColumn
  values::Vector{Bool}
  nulls::Vector{UInt8}
  TBoolColumn() = (o=new(); fillunset(o); o)
end # type TBoolColumn

type TByteColumn
  values::Vector{UInt8}
  nulls::Vector{UInt8}
  TByteColumn() = (o=new(); fillunset(o); o)
end # type TByteColumn

type TI16Column
  values::Vector{Int16}
  nulls::Vector{UInt8}
  TI16Column() = (o=new(); fillunset(o); o)
end # type TI16Column

type TI32Column
  values::Vector{Int32}
  nulls::Vector{UInt8}
  TI32Column() = (o=new(); fillunset(o); o)
end # type TI32Column

type TI64Column
  values::Vector{Int64}
  nulls::Vector{UInt8}
  TI64Column() = (o=new(); fillunset(o); o)
end # type TI64Column

type TDoubleColumn
  values::Vector{Float64}
  nulls::Vector{UInt8}
  TDoubleColumn() = (o=new(); fillunset(o); o)
end # type TDoubleColumn

type TStringColumn
  values::Vector{Compat.UTF8String}
  nulls::Vector{UInt8}
  TStringColumn() = (o=new(); fillunset(o); o)
end # type TStringColumn

type TBinaryColumn
  values::Vector{Vector{UInt8}}
  nulls::Vector{UInt8}
  TBinaryColumn() = (o=new(); fillunset(o); o)
end # type TBinaryColumn

type TColumn
  boolVal::TBoolColumn
  byteVal::TByteColumn
  i16Val::TI16Column
  i32Val::TI32Column
  i64Val::TI64Column
  doubleVal::TDoubleColumn
  stringVal::TStringColumn
  binaryVal::TBinaryColumn
  TColumn() = (o=new(); fillunset(o); o)
end # type TColumn

type TRowSet
  startRowOffset::Int64
  rows::Vector{TRow}
  columns::Vector{TColumn}
  TRowSet() = (o=new(); fillunset(o); o)
end # type TRowSet
meta(t::Type{TRowSet}) = meta(t, Symbol[:columns], Int[], Dict{Symbol,Any}())

type TStatus
  statusCode::Int32
  infoMessages::Vector{Compat.UTF8String}
  sqlState::Compat.UTF8String
  errorCode::Int32
  errorMessage::Compat.UTF8String
  TStatus() = (o=new(); fillunset(o); o)
end # type TStatus
meta(t::Type{TStatus}) = meta(t, Symbol[:infoMessages,:sqlState,:errorCode,:errorMessage], Int[], Dict{Symbol,Any}())

type THandleIdentifier
  guid::Vector{UInt8}
  secret::Vector{UInt8}
  THandleIdentifier() = (o=new(); fillunset(o); o)
end # type THandleIdentifier

type TSessionHandle
  sessionId::THandleIdentifier
  TSessionHandle() = (o=new(); fillunset(o); o)
end # type TSessionHandle

type TOperationHandle
  operationId::THandleIdentifier
  operationType::Int32
  hasResultSet::Bool
  modifiedRowCount::Float64
  TOperationHandle() = (o=new(); fillunset(o); o)
end # type TOperationHandle
meta(t::Type{TOperationHandle}) = meta(t, Symbol[:modifiedRowCount], Int[], Dict{Symbol,Any}())

type TOpenSessionReq
  client_protocol::Int32
  username::Compat.UTF8String
  password::Compat.UTF8String
  configuration::Dict{Compat.UTF8String,Compat.UTF8String}
  TOpenSessionReq() = (o=new(); fillunset(o); o)
end # type TOpenSessionReq
meta(t::Type{TOpenSessionReq}) = meta(t, Symbol[:username,:password,:configuration], Int[], Dict{Symbol,Any}(:client_protocol => Int32(7)))

type TOpenSessionResp
  status::TStatus
  serverProtocolVersion::Int32
  sessionHandle::TSessionHandle
  configuration::Dict{Compat.UTF8String,Compat.UTF8String}
  TOpenSessionResp() = (o=new(); fillunset(o); o)
end # type TOpenSessionResp
meta(t::Type{TOpenSessionResp}) = meta(t, Symbol[:sessionHandle,:configuration], Int[], Dict{Symbol,Any}(:serverProtocolVersion => Int32(7)))

type TCloseSessionReq
  sessionHandle::TSessionHandle
  TCloseSessionReq() = (o=new(); fillunset(o); o)
end # type TCloseSessionReq

type TCloseSessionResp
  status::TStatus
  TCloseSessionResp() = (o=new(); fillunset(o); o)
end # type TCloseSessionResp

type TGetInfoValue
  stringValue::Compat.UTF8String
  smallIntValue::Int16
  integerBitmask::Int32
  integerFlag::Int32
  binaryValue::Int32
  lenValue::Int64
  TGetInfoValue() = (o=new(); fillunset(o); o)
end # type TGetInfoValue

type TGetInfoReq
  sessionHandle::TSessionHandle
  infoType::Int32
  TGetInfoReq() = (o=new(); fillunset(o); o)
end # type TGetInfoReq

type TGetInfoResp
  status::TStatus
  infoValue::TGetInfoValue
  TGetInfoResp() = (o=new(); fillunset(o); o)
end # type TGetInfoResp

type TExecuteStatementReq
  sessionHandle::TSessionHandle
  statement::Compat.UTF8String
  confOverlay::Dict{Compat.UTF8String,Compat.UTF8String}
  runAsync::Bool
  TExecuteStatementReq() = (o=new(); fillunset(o); o)
end # type TExecuteStatementReq
meta(t::Type{TExecuteStatementReq}) = meta(t, Symbol[:confOverlay,:runAsync], Int[], Dict{Symbol,Any}(:runAsync => false))

type TExecuteStatementResp
  status::TStatus
  operationHandle::TOperationHandle
  TExecuteStatementResp() = (o=new(); fillunset(o); o)
end # type TExecuteStatementResp
meta(t::Type{TExecuteStatementResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetTypeInfoReq
  sessionHandle::TSessionHandle
  TGetTypeInfoReq() = (o=new(); fillunset(o); o)
end # type TGetTypeInfoReq

type TGetTypeInfoResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetTypeInfoResp() = (o=new(); fillunset(o); o)
end # type TGetTypeInfoResp
meta(t::Type{TGetTypeInfoResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetCatalogsReq
  sessionHandle::TSessionHandle
  TGetCatalogsReq() = (o=new(); fillunset(o); o)
end # type TGetCatalogsReq

type TGetCatalogsResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetCatalogsResp() = (o=new(); fillunset(o); o)
end # type TGetCatalogsResp
meta(t::Type{TGetCatalogsResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetSchemasReq
  sessionHandle::TSessionHandle
  catalogName::TIdentifier
  schemaName::TPatternOrIdentifier
  TGetSchemasReq() = (o=new(); fillunset(o); o)
end # type TGetSchemasReq
meta(t::Type{TGetSchemasReq}) = meta(t, Symbol[:catalogName,:schemaName], Int[], Dict{Symbol,Any}())

type TGetSchemasResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetSchemasResp() = (o=new(); fillunset(o); o)
end # type TGetSchemasResp
meta(t::Type{TGetSchemasResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetTablesReq
  sessionHandle::TSessionHandle
  catalogName::TPatternOrIdentifier
  schemaName::TPatternOrIdentifier
  tableName::TPatternOrIdentifier
  tableTypes::Vector{Compat.UTF8String}
  TGetTablesReq() = (o=new(); fillunset(o); o)
end # type TGetTablesReq
meta(t::Type{TGetTablesReq}) = meta(t, Symbol[:catalogName,:schemaName,:tableName,:tableTypes], Int[], Dict{Symbol,Any}())

type TGetTablesResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetTablesResp() = (o=new(); fillunset(o); o)
end # type TGetTablesResp
meta(t::Type{TGetTablesResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetTableTypesReq
  sessionHandle::TSessionHandle
  TGetTableTypesReq() = (o=new(); fillunset(o); o)
end # type TGetTableTypesReq

type TGetTableTypesResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetTableTypesResp() = (o=new(); fillunset(o); o)
end # type TGetTableTypesResp
meta(t::Type{TGetTableTypesResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetColumnsReq
  sessionHandle::TSessionHandle
  catalogName::TIdentifier
  schemaName::TPatternOrIdentifier
  tableName::TPatternOrIdentifier
  columnName::TPatternOrIdentifier
  TGetColumnsReq() = (o=new(); fillunset(o); o)
end # type TGetColumnsReq
meta(t::Type{TGetColumnsReq}) = meta(t, Symbol[:catalogName,:schemaName,:tableName,:columnName], Int[], Dict{Symbol,Any}())

type TGetColumnsResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetColumnsResp() = (o=new(); fillunset(o); o)
end # type TGetColumnsResp
meta(t::Type{TGetColumnsResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetFunctionsReq
  sessionHandle::TSessionHandle
  catalogName::TIdentifier
  schemaName::TPatternOrIdentifier
  functionName::TPatternOrIdentifier
  TGetFunctionsReq() = (o=new(); fillunset(o); o)
end # type TGetFunctionsReq
meta(t::Type{TGetFunctionsReq}) = meta(t, Symbol[:catalogName,:schemaName], Int[], Dict{Symbol,Any}())

type TGetFunctionsResp
  status::TStatus
  operationHandle::TOperationHandle
  TGetFunctionsResp() = (o=new(); fillunset(o); o)
end # type TGetFunctionsResp
meta(t::Type{TGetFunctionsResp}) = meta(t, Symbol[:operationHandle], Int[], Dict{Symbol,Any}())

type TGetOperationStatusReq
  operationHandle::TOperationHandle
  TGetOperationStatusReq() = (o=new(); fillunset(o); o)
end # type TGetOperationStatusReq

type TGetOperationStatusResp
  status::TStatus
  operationState::Int32
  sqlState::Compat.UTF8String
  errorCode::Int32
  errorMessage::Compat.UTF8String
  TGetOperationStatusResp() = (o=new(); fillunset(o); o)
end # type TGetOperationStatusResp
meta(t::Type{TGetOperationStatusResp}) = meta(t, Symbol[:operationState,:sqlState,:errorCode,:errorMessage], Int[], Dict{Symbol,Any}())

type TCancelOperationReq
  operationHandle::TOperationHandle
  TCancelOperationReq() = (o=new(); fillunset(o); o)
end # type TCancelOperationReq

type TCancelOperationResp
  status::TStatus
  TCancelOperationResp() = (o=new(); fillunset(o); o)
end # type TCancelOperationResp

type TCloseOperationReq
  operationHandle::TOperationHandle
  TCloseOperationReq() = (o=new(); fillunset(o); o)
end # type TCloseOperationReq

type TCloseOperationResp
  status::TStatus
  TCloseOperationResp() = (o=new(); fillunset(o); o)
end # type TCloseOperationResp

type TGetResultSetMetadataReq
  operationHandle::TOperationHandle
  TGetResultSetMetadataReq() = (o=new(); fillunset(o); o)
end # type TGetResultSetMetadataReq

type TGetResultSetMetadataResp
  status::TStatus
  schema::TTableSchema
  TGetResultSetMetadataResp() = (o=new(); fillunset(o); o)
end # type TGetResultSetMetadataResp
meta(t::Type{TGetResultSetMetadataResp}) = meta(t, Symbol[:schema], Int[], Dict{Symbol,Any}())

type TFetchResultsReq
  operationHandle::TOperationHandle
  orientation::Int32
  maxRows::Int64
  fetchType::Int16
  TFetchResultsReq() = (o=new(); fillunset(o); o)
end # type TFetchResultsReq
meta(t::Type{TFetchResultsReq}) = meta(t, Symbol[:fetchType], Int[], Dict{Symbol,Any}(:orientation => Int32(0), :fetchType => Int16(0)))

type TFetchResultsResp
  status::TStatus
  hasMoreRows::Bool
  results::TRowSet
  TFetchResultsResp() = (o=new(); fillunset(o); o)
end # type TFetchResultsResp
meta(t::Type{TFetchResultsResp}) = meta(t, Symbol[:hasMoreRows,:results], Int[], Dict{Symbol,Any}())

type TGetDelegationTokenReq
  sessionHandle::TSessionHandle
  owner::Compat.UTF8String
  renewer::Compat.UTF8String
  TGetDelegationTokenReq() = (o=new(); fillunset(o); o)
end # type TGetDelegationTokenReq

type TGetDelegationTokenResp
  status::TStatus
  delegationToken::Compat.UTF8String
  TGetDelegationTokenResp() = (o=new(); fillunset(o); o)
end # type TGetDelegationTokenResp
meta(t::Type{TGetDelegationTokenResp}) = meta(t, Symbol[:delegationToken], Int[], Dict{Symbol,Any}())

type TCancelDelegationTokenReq
  sessionHandle::TSessionHandle
  delegationToken::Compat.UTF8String
  TCancelDelegationTokenReq() = (o=new(); fillunset(o); o)
end # type TCancelDelegationTokenReq

type TCancelDelegationTokenResp
  status::TStatus
  TCancelDelegationTokenResp() = (o=new(); fillunset(o); o)
end # type TCancelDelegationTokenResp

type TRenewDelegationTokenReq
  sessionHandle::TSessionHandle
  delegationToken::Compat.UTF8String
  TRenewDelegationTokenReq() = (o=new(); fillunset(o); o)
end # type TRenewDelegationTokenReq

type TRenewDelegationTokenResp
  status::TStatus
  TRenewDelegationTokenResp() = (o=new(); fillunset(o); o)
end # type TRenewDelegationTokenResp

abstract TCLIServiceClientBase


text:
"""
Identify the driver that can open a raster file.

This function will try to identify the driver that can open the passed filename
by invoking the Identify method of each registered `GDALDriver` in turn. The
first driver that successful identifies the file name will be returned. If all
drivers fail then `NULL` is returned.

It is possible to give an optional list of files. This is the list of all files
at the same level in the file system as the target file, including the target
file. The filenames will not include any path components, and are essentially
just the output of `CPLReadDir()` on the parent directory. If the target object
does not have filesystem semantics then the file list should be `NULL`.

### Parameters
* `pszFilename`   the name of the file to access. In the case of exotic drivers
this may not refer to a physical file, but instead contain information for the
driver on how to access a dataset.
* `papszFileList` an array of strings, whose last element is the `NULL`
pointer. These strings are filenames that are auxiliary to the main filename.
The passed value may be `NULL`.

### Returns
A `GDALDriverH` handle or `NULL` on failure. For C++ applications this handle can
be cast to a `GDALDriver *`.
"""
_identifydriver(pszFilename::Ptr{UInt8}, papszFileList::Ptr{Ptr{UInt8}}) =
    GDALIdentifyDriver(pszFilename, papszFileList)::GDALDriverH

function identifydriver(filename::String,
                        filelist::Vector{String} = Vector{String}())
    driver = _identifydriver(pointer(filename),
                             Ptr{Ptr{UInt8}}(pointer(filelist)))
    (driver == C_NULL) && error("Could not identify driver")
    driver
end

"""
Fetch a driver based on the short name (such as `GTiff`).

Returns `NULL` if no match is found.
"""
_getdriverbyname(drivername::Ptr{UInt8}) =
    GDALGetDriverByName(drivername)::GDALDriverH

function driverbyname(drivername::String)
    driver = _getdriverbyname(pointer(drivername))
    (driver == C_NULL) && error("Could not find driver $drivername")
    driver
end

"Fetch the number of registered drivers."
_getdrivercount() = GDALGetDriverCount()

drivercount() = _getdrivercount()

"""
Fetch driver by index (from `0` to `GetDriverCount()-1`).

Throws an error if the index is invalid.
"""
_getdriver(i::Integer) = GDALGetDriver(i)

function driverbyindex(i::Integer)
    driver = _getdriver(i-1)
    (driver == C_NULL) && error("driver index $i is invalid")
    driver
end

"""
Return the short name of a driver (e.g. "GTiff")

This name can be passed to the `GDALGetDriverByName()` function.
The returned string should not be freed and is owned by the driver.
"""
_getdrivershortname(ptr::GDALDriverH) = GDALGetDriverShortName(ptr)::Ptr{UInt8}

drivershortname(ptr::GDALDriverH) = unsafe_string(_getdrivershortname(ptr))
drivershortname(i::Integer) = drivershortname(driverbyindex(i))

"""
Return the long name of a driver (e.g. "GeoTIFF"), or empty string.

The returned string should not be freed and is owned by the driver.
"""
_getdriverlongname(ptr::GDALDriverH) = GDALGetDriverLongName(ptr)

driverlongname(ptr::GDALDriverH) = unsafe_string(_getdriverlongname(ptr))
driverlongname(i::Integer) = driverlongname(driverbyindex(i))

"""
Fetch the driver that the dataset was created with `GDALOpen()`/`GDALCreate()`.
"""
_getdatasetdriver(dataset::GDALDatasetH) = GDALGetDatasetDriver(dataset)

"""
Destroy a `GDALDriver`.

This is roughly equivelent to deleting the driver, but is guaranteed to take
place in the GDAL heap. It is important this that function not be called on a
driver that is registered with the `GDALDriverManager`.

### Parameters
`hDriver`     the driver to destroy.
"""
_destroydriver(hDriver::GDALDriverH) = GDALDestroyDriver(hDriver)

"Register a driver for use."
_registerdriver(driver::GDALDriverH) = GDALRegisterDriver(driver)::Cint

"Deregister the passed driver."
_deregisterdriver(driver::GDALDriverH) = GDALDeregisterDriver(driver)

# """
# Destroy the driver manager. Incidently unloads all managed drivers.

# NOTE: This function is not thread safe. It should not be called while other
# threads are actively using GDAL.
# """
# _destroydrivermanager() = GDALDestroyDriverManager()

"""
Return the URL to the help that describes the driver.

That URL is relative to the GDAL documentation directory.
For the GeoTIFF driver, this is "frmt_gtiff.html"

### Parameters
* `hDriver`     the handle of the driver

### Returns
the URL to the help that describes the driver or `NULL`. The returned string
should not be freed and is owned by the driver.
"""
_getdriverhelptopic(driver::GDALDriverH) =
    GDALGetDriverHelpTopic(driver)::Ptr{UInt8}

driverhelptopic(driver::GDALDriverH) = unsafe_string(_getdriverhelptopic(driver))

"""
Return the list of creation options of the driver used by `Create()` and
`CreateCopy()` as an XML string

### Parameters
* `hDriver`     the handle of the driver

### Returns
an XML string that describes the list of creation options or empty string.
The returned string should not be freed and is owned by the driver.
"""
_getdrivercreationoptionlist(driver::GDALDriverH) =
    GDALGetDriverCreationOptionList(driver)::Ptr{UInt8}

drivercreationoptionlist(driver::GDALDriverH) =
    unsafe_string(_getdrivercreationoptionlist(driver))


text:
using Test
using LatticeTools

@testset "Lattice1D" begin
    unitcell = makeunitcell(2.0; SiteType=String)
    addsite!(unitcell, "Spin", FractCoord([0], [0.0]))
    lattice = make_lattice(unitcell, 3)
    scale_matrix = 3 * ones(Int, (1,1))
    @test lattice.unitcell == unitcell
    #@test lattice.hypercube == HypercubicLattice(scale_matrix)
    @test lattice.hypercube == Hypercube(scale_matrix)
    @test numsite(lattice.supercell) == 3

    @test_throws DimensionMismatch make_lattice(unitcell, [3 0])
    @test_throws DimensionMismatch make_lattice(unitcell, [3 0; 0 3])
end

@testset "Lattice" begin
    unitcell = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
    addsite!(unitcell, "d", FractCoord([0,0], [0.0, 0.0]))
    addsite!(unitcell, "px", FractCoord([0,0], [0.5, 0.0]))
    addsite!(unitcell, "py", FractCoord([0,0], [0.0, 0.5]))

    @test_throws DimensionMismatch make_lattice(unitcell, 3)
    @test_throws DimensionMismatch make_lattice(unitcell, [1 0 0; 0 1 0; 0 0 1])
    @test_throws DimensionMismatch make_lattice(unitcell, [1 0 0; 0 1 0; 0 0 1], [1 0; 0 1])

    @test_throws DimensionMismatch make_lattice(unitcell, [3 0; 0 3], ones(Int, 1, 1))
    @test_throws DimensionMismatch make_lattice(unitcell, [3 0; 0 3], [1 0 0; 0 1 0; 0 0 1])
    @test_throws ArgumentError make_lattice(unitcell, [4 0; 0 3], [1 1; 0 1]) # generators not orthogonal
    @test_throws ArgumentError make_lattice(unitcell, [4 0; 0 4], [1 0; 0 2]) # generator not unimodular

    lattice = make_lattice(unitcell, [2 0; 0 2])

    @test lattice.unitcell == unitcell
    @test lattice.hypercube.shape_matrix == [2 0; 0 2]
    @test numsite(lattice.supercell) == 3 * 2 * 2
    @test lattice.supercell.latticevectors == [2.0 0.0; 0.0 2.0]

    @test lattice == make_lattice(unitcell, [2 0; 0 2])
    @test lattice != make_lattice(unitcell, [3 0; 0 2])

    let lattice2 = make_lattice(unitcell, [2 0; 0 2], [1 0; 0 1])
        @test lattice == lattice2
    end
    let lattice3 = make_lattice(unitcell, [2 0; 0 2], [0 1; 1 0])
        @test lattice != lattice3
        @test lattice.unitcell == lattice3.unitcell
        @test lattice.hypercube == lattice3.hypercube
        @test lattice.bravais_coordinates != lattice3.bravais_coordinates
    end
    let unitcell2 = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String)
        addsite!(unitcell2, "d", FractCoord([0,0], [0.0, 0.0]))
        @test lattice != make_lattice(unitcell2, [2 0; 0 2])
    end
end

@testset "LatticeWithOrbital" begin
    unitcell = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String, OrbitalType=String)
    addsite!(unitcell, "Cu", FractCoord([0,0], [0.0, 0.0]))
    addsite!(unitcell, "Ox", FractCoord([0,0], [0.5, 0.0]))
    addsite!(unitcell, "Oy", FractCoord([0,0], [0.0, 0.5]))

    addorbital!(unitcell, "d", "Cu")
    addorbital!(unitcell, "px", "Ox")
    addorbital!(unitcell, "py", 3)

    @test_throws DimensionMismatch make_lattice(unitcell, 3)
    @test_throws DimensionMismatch make_lattice(unitcell, [1 0 0; 0 1 0; 0 0 1])
    @test_throws DimensionMismatch make_lattice(unitcell, [1 0 0; 0 1 0; 0 0 1], [1 0; 0 1])

    @test_throws DimensionMismatch make_lattice(unitcell, [3 0; 0 3], ones(Int, 1, 1))
    @test_throws DimensionMismatch make_lattice(unitcell, [3 0; 0 3], [1 0 0; 0 1 0; 0 0 1])
    @test_throws ArgumentError make_lattice(unitcell, [4 0; 0 3], [1 1; 0 1]) # generators not orthogonal
    @test_throws ArgumentError make_lattice(unitcell, [4 0; 0 4], [1 0; 0 2]) # generator not unimodular

    lattice = make_lattice(unitcell, [2 0; 0 2])

    @test lattice.unitcell == unitcell
    @test lattice.hypercube.shape_matrix == [2 0; 0 2]
    @test numsite(lattice.supercell) == 3 * 2 * 2
    @test lattice.supercell.latticevectors == [2.0 0.0; 0.0 2.0]

    @test lattice == make_lattice(unitcell, [2 0; 0 2])
    @test lattice != make_lattice(unitcell, [3 0; 0 2])

    let lattice2 = make_lattice(unitcell, [2 0; 0 2], [1 0; 0 1])
        @test lattice == lattice2
    end
    let lattice3 = make_lattice(unitcell, [2 0; 0 2], [0 1; 1 0])
        @test lattice != lattice3
        @test lattice.unitcell == lattice3.unitcell
        @test lattice.hypercube == lattice3.hypercube
        @test lattice.bravais_coordinates != lattice3.bravais_coordinates
    end
    let unitcell2 = makeunitcell([1.0 0.0; 0.0 1.0]; SiteType=String, OrbitalType=String)
        addsite!(unitcell2, "d", FractCoord([0,0], [0.0, 0.0]))
        @test lattice != make_lattice(unitcell2, [2 0; 0 2])
    end

    @test numsite(lattice.supercell) == 12
    for (siteindex, (sitename, sitecoord)) in enumerate(lattice.supercell.sites)
        @test hassite(lattice.unitcell, sitename[1])
        @test isa(sitename[2], Vector{Int})
        @test getsite(lattice.supercell, siteindex) == (sitename, sitecoord)
        @test getsitename(lattice.supercell, siteindex) == sitename
        @test getsitecoord(lattice.supercell, siteindex) == sitecoord
    end

    @test numorbital(lattice.supercell) == 12
    for (orbitalindex, (orbitalname, siteindex)) in enumerate(lattice.supercell.orbitals)
        @test hasorbital(lattice.unitcell, orbitalname[1])
        @test isa(orbitalname[2], Vector{Int})
        small_orbitalindex = getorbitalindex(lattice.unitcell, orbitalname[1])
        small_siteindex = getorbitalsiteindex(lattice.unitcell, small_orbitalindex)
        small_sitename = getsitename(lattice.unitcell, small_siteindex)
        @test small_sitename == getsitename(lattice.supercell, siteindex)[1]
        @test getsitename(lattice.supercell, siteindex) == (small_sitename, orbitalname[2])
    end
end


text:
module CardUtils

using StatsBase

export ps2pdf, gen_rand_array

# Converts PostScript files to PDF files.
# Requires GhostScript to be installed and the `ps2pdf` utility to be in the
# program's $PATH.
ps2pdf(ps) = open(`ps2pdf - -`, read=true, write=true) do proc
  write(proc.in, ps)
  close(proc.in)
  read(proc.out, String)
end

# TODO: Write tests
gen_rand_col(rng, yl, replace) = (r) -> sample(rng, r, yl, replace=replace)

# TODO: Write tests
gen_rand_array(rng; ranges=[01:15, 16:30, 31:45, 46:65, 66:75], yl=5,
               replace=false) = hcat(map(gen_rand_col(rng, yl, replace), ranges)...)

end

max_stars_repo_name: tuckertwo/ECbingo

max_issues_repo_name: tuckertwo/ECbingo

max_forks_repo_name: tuckertwo/ECbingo


text:
include("Elements.jl");

# numerical fluxes
central(ρ⁻, ρ⁺) = (ρ⁻ + ρ⁺) / 2
right(ρ⁻, ρ⁺) = ρ⁺
left(ρ⁻, ρ⁺) = ρ⁻

# DG D4 operators
function D_DG(ρ, JI, LIFT, D, numflux)
  K = size(ρ, 2)
  # element derivative
  Dρ = JI * D * ρ

  # Face continuity
  ρ⁻ = @view ρ[end, :]
  ρ⁺ = @view ρ[ 1, mod1.(2:K+1, K)]
  Dρ[end, :] += LIFT * (numflux(ρ⁻, ρ⁺) - ρ⁻)

  ρ⁻ = @view ρ[1, :]
  ρ⁺ = @view ρ[end, mod1.(0:K-1, K)]
  Dρ[1, :] -= LIFT * (numflux(ρ⁺, ρ⁻) - ρ⁻)

  return Dρ
end

function dg_data(N, K, start::FT, stop::FT;
                 numflux = (central, central, central, central)) where FT
  elems = range(start, stop = stop, length=K + 1)

  # Grid spacing
  Δx = FT(elems.step)

  # LGL points and weights
  ξ, ω = lglpoints(FT, N)

  # Get the derivative matrix
  D = spectralderivative(ξ)

  # Create the degree of freedom mesh
  x = Δx * (ξ .+ 1) / 2 .+ elems[1:end-1]'

  # Inverse of the Jacobian determinant
  JI = 2 / Δx

  # lift operator
  LIFT = JI / ω[1]

  D4 = zeros(FT, (N+1) * K, (N+1) * K)

  v = zeros(FT, N+1, K)

  for n = 1:length(v)
    v[n] = 1
    Dv = D_DG(v, JI, LIFT, D, numflux[1])
    D2v = D_DG(Dv, JI, LIFT, D, numflux[2])
    D3v = D_DG(D2v, JI, LIFT, D, numflux[3])
    D4[:, n] = D_DG(D3v, JI, LIFT, D, numflux[4])[:]
    v[n] = 0
  end

  return (D4 = D4, x = x, M = ω  / JI)
end


text:
function miscCombineTimeSteps(base_expr::Basic, base_radius::Int,
                            expr::Basic, radius::Int)

    b_sym, b_off = CreateSymOffMatrices(base_radius)
    sym,off = CreateSymOffMatrices(radius)
    new_sym, new_off = CreateSymOffMatrices(radius+base_radius)

    new_eq = 0

    for i = 1:(2*radius+1), j = 1:(2*radius+1),k = 1:(2*radius+1)
        c1 = coeff(expr, sym[i,j,k])
        (c1 == 0) && continue
        for ii = 1:(2*base_radius+1), jj = 1:(2*base_radius+1),kk = 1:(2*base_radius+1)
            c2 = coeff(base_expr, b_sym[ii,jj,kk])
            (c2 == 0) && continue
            off2 = b_off[ii,jj,kk] + [i;j;k] .+ base_radius
            new_eq += c1*c2*new_sym[off2...]
        end
    end
    return new_eq
end

function countComp(expr, radius, m_step)
    count_arr = zeros(Int, m_step)
    r = radius
    new_sym = expr
    count_arr[1] = count_coeffs(expr, radius)
    for i = 2:m_step
        new_sym = CombineTimeSteps(expr, r, new_sym, r*i)
        count_arr[i] = count_coeffs(new_sym, r*i)
    end
    return count_arr
end

function count_coeffs(expr, radius)
    Dsym = [symbols("w$(i)_$(j)_$(k)") for i = 1:(2radius+1),
                                        j = 1:(2radius+1),
                                        k = 1:(2radius+1)]
    s = 0
    for d in Dsym
        c = coeff(expr, d)
        if c  != 0
            s += 1
        end
    end
    return s

end