text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

module Resolve

using ..Versions

using Printf
using Random
using UUIDs

export resolve, sanity_check, Graph, pkgID

####################
# Requires / Fixed #
####################
const Requires = Dict{UUID,VersionSpec}

struct Fixed
    version::VersionNumber
    requires::Requires
end
Fixed(v::VersionNumber) = Fixed(v, Requires())

Base.:(==)(a::Fixed, b::Fixed) = a.version == b.version && a.requires == b.requires
Base.hash(f::Fixed, h::UInt) = hash((f.version, f.requires), h + (0x68628b809fd417ca % UInt))

Base.show(io::IO, f::Fixed) = isempty(f.requires) ?
    print(io, "Fixed(", repr(f.version), ")") :
    print(io, "Fixed(", repr(f.version), ",", f.requires, ")")

struct ResolverError <: Exception
    msg::AbstractString
    ex::Union{Exception,Nothing}
end
ResolverError(msg::AbstractString) = ResolverError(msg, nothing)

function Base.showerror(io::IO, pkgerr::ResolverError)
    print(io, pkgerr.msg)
    if pkgerr.ex !== nothing
        pkgex = pkgerr.ex
        if isa(pkgex, CompositeException)
            for cex in pkgex
                print(io, "\n=> ")
                showerror(io, cex)
            end
        else
            print(io, "\n")
            showerror(io, pkgex)
        end
    end
end

const uuid_julia = UUID("1222c4b2-2114-5bfd-aeef-88e4692bbb3e")

include("versionweights.jl")
include("fieldvalues.jl")
include("graphtype.jl")
include("maxsum.jl")

"Resolve package dependencies."
function resolve(graph::Graph)
    # attempt trivial solution first
    greedy_ok, sol = greedysolver(graph)

    greedy_ok && @goto solved

    log_event_global!(graph, "greedy solver failed")

    # trivial solution failed, use maxsum solver
    maxsum_ok, sol, staged = maxsum(graph)

    maxsum_ok && @goto solved

    log_event_global!(graph, "maxsum solver failed")

    # the problem is unsat, force-trigger a failure
    # in order to produce a log - this will contain
    # information about the best that the solver could
    # achieve
    trigger_failure!(graph, sol, staged)

    @label solved

    # verify solution (debug code) and enforce its optimality
    @assert verify_solution(sol, graph)
    greedy_ok || enforce_optimality!(sol, graph)

    log_event_global!(graph, "the solver found $(greedy_ok ? "an optimal" : "a feasible") configuration")

    # return the solution as a Dict mapping UUID => VersionNumber
    return compute_output_dict(sol, graph)
end

"""
Scan the graph for (explicit or implicit) contradictions. Returns a list of problematic
(package,version) combinations.
"""
function sanity_check(graph::Graph, sources::Set{UUID} = Set{UUID}(), verbose::Bool = true)
    req_inds = graph.req_inds
    fix_inds = graph.fix_inds

    id(p) = pkgID(p, graph)

    isempty(req_inds) || @warn("sanity check called on a graph with non-empty requirements")
    if !any(is_julia(graph, fp0) for fp0 in fix_inds)
        @warn("sanity check called on a graph without julia requirement, adding it")
        add_fixed!(graph, Dict(uuid_julia=>Fixed(VERSION)))
    end
    if length(fix_inds) ≠ 1
        @warn("sanity check called on a graph with extra fixed requirements (besides julia)")
    end

    isources = isempty(sources) ?
        Set{Int}(1:graph.np) :
        Set{Int}(graph.data.pdict[p] for p in sources)

    simplify_graph!(graph, isources)

    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    data = graph.data
    pkgs = data.pkgs
    pdict = data.pdict
    pvers = data.pvers
    eq_classes = data.eq_classes

    problematic = Tuple{String,VersionNumber}[]

    np == 0 && return problematic

    vers = [(pkgs[p0],pvers[p0][v0]) for p0 = 1:np for v0 = 1:(spp[p0]-1)]
    sort!(vers, by=pv->(-length(gadj[pdict[pv[1]]])))

    nv = length(vers)

    svdict = Dict{Tuple{UUID,VersionNumber},Int}(vers[i] => i for i = 1:nv)

    checked = falses(nv)

    last_str_len = 0

    for (i,(p,vn)) in enumerate(vers)
        if verbose
            frac_compl = i / nv
            print("\r", " "^last_str_len, "\r")
            progr_msg = @sprintf("%.3i/%.3i (%i%%) — problematic so far: %i", i, nv, round(Int, 100 * frac_compl), length(problematic))
            print(progr_msg)
            last_str_len = length(progr_msg)
        end

        length(gadj[pdict[p]]) == 0 && break
        checked[i] && continue

        push_snapshot!(graph)

        # enforce package version
        # TODO: use add_reqs! instead...
        p0 = graph.data.pdict[p]
        v0 = graph.data.vdict[p0][vn]
        fill!(graph.gconstr[p0], false)
        graph.gconstr[p0][v0] = true
        push!(graph.req_inds, p0)

        ok = false
        try
            simplify_graph_soft!(graph, Set{Int}([p0]), log_events = false)
        catch err
            isa(err, ResolverError) || rethrow()
            @goto done
        end

        ok, sol = greedysolver(graph)
        ok && @goto done
        ok, sol = maxsum(graph)

        @label done

        if !ok
            for vneq in eq_classes[p][vn]
                push!(problematic, (id(p), vneq))
            end
        else
            @assert verify_solution(sol, graph)
            sol_dict = compute_output_dict(sol, graph)
            for (sp,svn) in sol_dict
                j = svdict[sp,svn]
                checked[j] = true
            end
        end

        # state reset
        empty!(graph.req_inds)
        pop_snapshot!(graph)
    end
    if verbose
        print("\r", " "^last_str_len, "\r")
        println("found $(length(problematic)) problematic versions")
    end
    return sort!(problematic)
end

"""
Translate the solver output (a Vector{Int} of package states) into a Dict which
associates a VersionNumber to each installed package UUID.
"""
function compute_output_dict(sol::Vector{Int}, graph::Graph)
    np = graph.np
    spp = graph.spp
    fix_inds = graph.fix_inds
    pkgs = graph.data.pkgs
    pvers = graph.data.pvers
    pruned = graph.data.pruned

    want = Dict{UUID,VersionNumber}()
    for p0 = 1:np
        p0 ∈ fix_inds && continue
        p = pkgs[p0]
        s0 = sol[p0]
        s0 == spp[p0] && continue
        vn = pvers[p0][s0]
        want[p] = vn
    end
    for (p,vn) in pruned
        @assert !haskey(want, p)
        want[p] = vn
    end

    return want
end

"""
Preliminary solver attempt: tries to maximize each version; bails out as soon as
some non-trivial requirement is detected.
"""
function greedysolver(graph::Graph)
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    np = graph.np

    push_snapshot!(graph)
    gconstr = graph.gconstr

    # initialize solution: all uninstalled
    sol = [spp[p0] for p0 = 1:np]

    # packages which are not allowed to be uninstalled
    # (NOTE: this is potentially a superset of graph.req_inds,
    #        since it may include implicit requirements)
    req_inds = Set{Int}(p0 for p0 = 1:np if !gconstr[p0][end])

    # set up required packages to their highest allowed versions
    for rp0 in req_inds
        # look for the highest version which satisfies the requirements
        rv0 = findlast(gconstr[rp0])
        @assert rv0 ≠ 0 && rv0 ≠ spp[rp0]
        sol[rp0] = rv0
        fill!(gconstr[rp0], false)
        gconstr[rp0][rv0] = true
    end

    # propagate the requirements
    try
        simplify_graph_soft!(graph, req_inds, log_events = false)
    catch err
        err isa ResolverError || rethrow()
        pop_snapshot!(graph)
        return (false, Int[])
    end

    # we start from required packages and explore the graph
    # following dependencies
    staged = req_inds
    seen = copy(staged)

    while !isempty(staged)
        staged_next = Set{Int}()
        for p0 in staged
            s0 = sol[p0]
            @assert s0 < spp[p0]

            # scan dependencies
            for (j1,p1) in enumerate(gadj[p0])
                msk = gmsk[p0][j1]
                # look for the highest version which satisfies the requirements
                v1 = findlast(msk[:,s0] .& gconstr[p1])
                v1 == spp[p1] && continue # p1 is not required by p0's current version
                # if we found a version, and the package was uninstalled
                # or the same version was already selected, we're ok;
                # otherwise we can't be sure what the optimal configuration is
                # and we bail out
                if v1 > 0 && (sol[p1] == spp[p1] || sol[p1] == v1)
                    sol[p1] = v1
                else
                    pop_snapshot!(graph)
                    return (false, Int[])
                end

                p1 ∈ seen || push!(staged_next, p1)
            end
        end
        union!(seen, staged_next)
        staged = staged_next
    end

    pop_snapshot!(graph)

    for p0 = 1:np
        log_event_greedysolved!(graph, p0, sol[p0])
    end

    return true, sol
end

"""
Verifies that the solver solution fulfills all hard constraints
(requirements and dependencies). This is intended as debug code.
"""
function verify_solution(sol::Vector{Int}, graph::Graph)
    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    gconstr = graph.gconstr

    @assert length(sol) == np
    @assert all(sol .> 0)

    # verify constraints and dependencies
    for p0 = 1:np
        s0 = sol[p0]
        gconstr[p0][s0] || (@warn("gconstr[$p0][$s0] fail"); return false)
        for (j1,p1) in enumerate(gadj[p0])
            msk = gmsk[p0][j1]
            s1 = sol[p1]
            msk[s1,s0] || (@warn("gmsk[$p0][$p1][$s1,$s0] fail"); return false)
        end
    end
    return true
end

"""
Push the given solution to a local optimium if needed: keeps increasing
the states of the given solution as long as no constraints are violated.
It also removes unnecessary parts of the solution which are unconnected
to the required packages.
"""
function enforce_optimality!(sol::Vector{Int}, graph::Graph)
    np = graph.np
    spp = graph.spp
    gadj = graph.gadj
    gmsk = graph.gmsk
    gconstr = graph.gconstr
    pkgs = graph.data.pkgs

    # keep a track for the log
    why = Union{Symbol,Int}[0 for p0 = 1:np]

    restart = true
    while restart
        restart = false
        for p0 = 1:np
            s0 = sol[p0]
            s0 == spp[p0] && (why[p0] = :uninst; continue) # the package is not installed

            # check if bumping to the higher version would violate a constraint
            gconstr[p0][s0+1] || (why[p0] = :constr; continue)

            # check if bumping to the higher version would violate a constraint
            viol = false
            for (j1,p1) in enumerate(gadj[p0])
                s1 = sol[p1]
                msk = gmsk[p0][j1]
                if !msk[s1, s0+1]
                    viol = true
                    why[p0] = p1
                    break
                end
            end
            viol && continue

            # So the solution is non-optimal: we bump it manually
            sol[p0] += 1
            restart = true
        end
    end

    # Finally uninstall unneeded packages:
    # start from the required ones and keep only
    # the packages reachable from them along the graph.
    # (These should have been removed in the previous step, but in principle
    # an unconnected yet self-sustaining cycle may have survived.)
    uninst = trues(np)
    staged = Set{Int}(p0 for p0 = 1:np if !gconstr[p0][end])
    seen = copy(staged)

    while !isempty(staged)
        staged_next = Set{Int}()
        for p0 in staged
            s0 = sol[p0]
            @assert s0 < spp[p0]
            uninst[p0] = false
            for (j1,p1) in enumerate(gadj[p0])
                gmsk[p0][j1][end,s0] && continue # the package is not required by p0 at version s0
                p1 ∈ seen || push!(staged_next, p1)
            end
        end
        union!(seen, staged_next)
        staged = staged_next
    end

    for p0 in findall(uninst)
        sol[p0] = spp[p0]
        why[p0] = :uninst
    end

    @assert verify_solution(sol, graph)

    for p0 = 1:np
        log_event_maxsumsolved!(graph, p0, sol[p0], why[p0])
    end
end

function apply_maxsum_trace!(graph::Graph, sol::Vector{Int})
    gconstr = graph.gconstr

    for (p0,s0) in enumerate(sol)
        s0 == 0 && continue
        gconstr0 = gconstr[p0]
        old_constr = copy(gconstr0)
        @assert old_constr[s0]
        fill!(gconstr0, false)
        gconstr0[s0] = true
        gconstr0 ≠ old_constr && log_event_maxsumtrace!(graph, p0, s0)
    end
end

function trigger_failure!(graph::Graph, sol::Vector{Int}, staged::Tuple{Int,Int})
    apply_maxsum_trace!(graph, sol)
    simplify_graph_soft!(graph, Set(findall(sol .> 0)), log_events = true) # this may throw an error...

    np = graph.np
    gconstr = graph.gconstr
    p0, v0 = staged

    @assert gconstr[p0][v0]
    fill!(gconstr[p0], false)
    gconstr[p0][v0] = true
    log_event_maxsumtrace!(graph, p0, v0)
    simplify_graph!(graph) # this may throw an error...
    outdict = resolve(graph) # ...otherwise, this MUST throw an error
    open(io->showlog(io, graph, view=:chronological), "logchrono.errresolve.txt", "w")
    error("this is not supposed to happen... $(Dict(pkgID(p, graph) => vn for (p,vn) in outdict))")
end

end # module


text:
function convert_klines(klines)
    t     = []
    start = []
    high  = []
    low   = []
    out   = []
    for i in 1:1:size(klines,1)
        push!(t, unix2datetime(klines[i, 7]/1000))
        push!(start, parse(Float64, klines[i,2]))
        push!(high, parse(Float64, klines[i,3]))
        push!(low, parse(Float64, klines[i,4]))
        push!(out, parse(Float64, klines[i,5]))
    end
    return t, out, high, low
end

function klines2table(klines; time_warp=3600)
    open_time  = []
    open       = []
    high       = []
    low        = []
    close      = []
    volume     = []
    close_time = []
    for i in 1:1:size(klines,1)
        push!(open_time, unix2datetime(klines[i, 1]/1000 + time_warp))
        push!(open, parse(Float64, klines[i,2]))
        push!(high, parse(Float64, klines[i,3]))
        push!(low, parse(Float64, klines[i,4]))
        push!(close, parse(Float64, klines[i,5]))
        push!(volume, parse(Float64, klines[i,6]))
        push!(close_time, unix2datetime(klines[i, 7]/1000 + time_warp))
    end
    return Table(
        open_time = Array{DateTime,1}(open_time),
        open = Array{Float64,1}(open),
        high = Array{Float64,1}(high),
        low = Array{Float64,1}(low),
        close = Array{Float64,1}(close),
        close_time = Array{DateTime,1}(close_time)
    )
end

function get_klines(client, symbol, kline, datapoints; time_warp=3600)
    timespans = Dict()
    timespans["1m"] = [1, "minute"]
    timespans["3m"] = [3, "minute"]
    timespans["5m"] = [5, "minute"]
    timespans["15m"] = [15, "minute"]
    timespans["30m"] = [30, "minute"]

    timespans["1h"] = [1, "hour"]
    timespans["2h"] = [2, "hour"]
    timespans["4h"] = [4, "hour"]
    timespans["6h"] = [6, "hour"]
    timespans["8h"] = [8, "hour"]
    timespans["12h"] = [12, "hour"]

    timespans["1d"] = [1, "day"]
    timespans["3d"] = [3, "day"]

    timespans["1w"] = [1, "week"]

    timespans["1M"] = [1, "month"]
    kline_type = split(kline, "_")[end]
    kline_historical = string(datapoints * timespans[kline_type][1]) * " " * timespans[kline_type][2]
    sleep(2)
    klines = client.get_historical_klines(symbol, kline_type, kline_historical * " ago UTC");
    return Pinguu.klines2table(klines; time_warp=time_warp)
end

ext: jl


text:
let
    using Random
    Random.seed!(0)

    n = 50
    p = 2
    t = (1:n)/n
    b = randn(n)
    Ut,Vt = egrss.generators(t,p);
    K = egrss.full(Ut,Vt)

    Lref = cholesky(K).L
    Wt = egrss.potrf(Ut,Vt)
    xref = Lref\b
    x = egrss.trsv(Ut,Wt,b)
    @test ( isapprox(x,xref) )
    x = egrss.trsv(Ut,Wt,b,'N')
    @test ( isapprox(x,xref) )
    xref = Lref'\b
    x = egrss.trsv(Ut,Wt,b,'T')
    @test ( isapprox(x,xref) )

    d = collect((1:n)/n) .+ 1e-3
    Lref = cholesky(egrss.full(Ut,Vt,d)).L;
    Wt,c = egrss.potrf(Ut,Vt,d)
    xref = Lref\b
    x = egrss.trsv(Ut,Wt,c,b)
    @test ( isapprox(x,xref) )
    x = egrss.trsv(Ut,Wt,c,b,'N');
    @test ( isapprox(x,xref) )
    xref = Lref'\b;
    x = egrss.trsv(Ut,Wt,c,b,'T');
    @test ( isapprox(x,xref) )
end

max_stars_repo_name: martinandersen/EGRSS

max_issues_repo_name: martinandersen/EGRSS

max_forks_repo_name: martinandersen/EGRSS


text:
module MatrixLH

export round_to_grid

"""
accumarray for linear indices

Due to Tim Holy
"https://groups.google.com/forum/#!topic/julia-users/QwdFHbbasis"

IN
    subs :: Vector
"""
function accumarrayTH(subs, val, sz=(maximum(subs),))
    A = zeros(eltype(val), sz...)
    for i = 1:length(val)
        A[subs[i]] += val[i]
    end
    A
end

"""
accumarray for arrays

"https://groups.google.com/forum/#!topic/julia-users/QwdFHbbasis"
"""
function accumarray(subs, val :: Array{T,1}, fun=sum,
    sz=maximum(subs, dims=1))  where T <: Number

   counts = Dict()
   for i = 1:size(subs,1)
       # Append value
       counts[subs[i,:]] = [get(counts, subs[i,:], []); val[i...]]
   end

   A = zeros(T, sz...);
   for j = keys(counts)
        A[j...] = fun(counts[j])
   end
   return A
end

"""
round_to_grid()

Round a matrix to the nearest points on a grid
"""
function round_to_grid(xM :: Array{T}, gridV :: Vector{T}) where T <: AbstractFloat
    idxM = similar(xM, Int);
    for (i1, x) in enumerate(xM)
        iMin = findmin(abs.(gridV .- x));
        idxM[i1] = iMin[2];
    end
    return idxM
end

end # module


text:
module AffineFEOperatorsTests

using Test
using Gridap.Arrays
using Gridap.Algebra
using Gridap.TensorValues
using Gridap.ReferenceFEs
using Gridap.Geometry
using Gridap.Fields
using Gridap.FESpaces
using LinearAlgebra
using Gridap.CellData

domain =(0,1,0,1)
partition = (2,2)
model = CartesianDiscreteModel(domain,partition)
order = 1
reffe = ReferenceFE(lagrangian,Float64,order)
V = FESpace(model,reffe;dirichlet_tags=[1,10])
U = V

trian = get_triangulation(model)
degree = 4
quad = CellQuadrature(trian,degree)

#f(x) = x[2]
#v = get_cell_shapefuns(V)
#u = get_cell_shapefuns_trial(U)
#
#cellmat = integrate(∇(v)⊙∇(u),quad)
#cellvec = integrate(v*f,quad)
#cellids = collect(1:num_cells(trian))
#
#assem = SparseMatrixAssembler(U,V)
#matdata = ([cellmat],[cellids],[cellids])
#vecdata = ([cellvec],[cellids])
#A =  assemble_matrix(assem,matdata)
#b =  assemble_vector(assem,vecdata)
#
#op = AffineFEOperator(U,V,A,b)
#@test A === get_matrix(op)
#@test b === get_vector(op)
#
#x = ones(length(b))
#r = A*x - b
#
#test_fe_operator(op,x,r,≈,jac=A)
#@test isa(get_algebraic_operator(op), AffineOperator)

#

tol = 1.0e-9

u_sol(x) = x[1]+x[2]
f_fun(x) = 0

V = FESpace(model,reffe;dirichlet_tags="boundary")
U = TrialFESpace(V,u_sol)

dΩ = Measure(quad)

a(u,v) = ∫(∇(v)⊙∇(u))*dΩ
l(v) = ∫(v*f_fun)*dΩ

assem = SparseMatrixAssembler(U,V)

op = AffineFEOperator(U,V,assem) do u,v
  ∫(∇(v)⊙∇(u))*dΩ, ∫(v*f_fun)*dΩ
end
uh = solve(op)
e = u_sol - uh

@test sum(∫(e*e)*dΩ) < tol
@test ∑(∫(e*e)*dΩ) < tol

op = AffineFEOperator(U,V) do u,v
  ∫(∇(v)⊙∇(u))*dΩ, ∫(v*f_fun)*dΩ
end
uh = solve(op)
e = u_sol - uh
@test sum(∫(e*e)*dΩ) < tol

op = AffineFEOperator(a,l,U,V)
uh = solve(op)
e = u_sol - uh
@test sum(∫(e*e)*dΩ) < tol

end # module


text:
# tests for find_zeros

using Roots
using Compat.Test

mutable struct CallableFunction
f
n
end
(F::CallableFunction)(x) = (F.n+=1; F.f(x))

@testset "find_zeros" begin

    function azero(f, x)
        fx = f(x)
        iszero(fx) && return true
        sign(fx) * sign(f(nextfloat(x))) < 0 && return true
        sign(fx) * sign(f(prevfloat(x))) < 0 && return true
        abs(fx) <= 8eps(x) && return true
        false
    end

    ## easier ones, counting steps
    F = CallableFunction(x -> exp(x) - x^4, 0)
    xrts = find_zeros(F, -5, 20)
    @test length(xrts) == 3
    @test all(azero.((F,), xrts))
    @test F.n <= 3000

    F = CallableFunction(x -> cos(x) + cos(2x), 0)
    xrts = find_zeros(F, 0, 4pi)
    @test length(xrts) == 6
    @test all(azero.((F,), xrts))
    @test F.n <= 5000

    T11(x) = 1024x^11 -2816x^9 + 2816x^7 -1232x^5 + 220x^3 - 11x
    U9(x) = 512x^9 - 1024x^7 + 672x^5 - 160x^3 +10x

    F = CallableFunction(T11, 0)
    xrts = find_zeros(F, -1, 1)
    @test length(xrts) == 11
    @test all(azero.((F,), xrts))
    @test F.n <= 10_000

    F = CallableFunction(U9, 0)
    xrts = find_zeros(F, -1, 1)
    @test length(xrts) == 9
    @test all(azero.((F,), xrts))
    @test F.n <= 10_000

    W(n) = x -> prod(x-i for i in 1:n)
    Wi(n) = x -> prod((x-i)^i for i in 1:n)

    F = CallableFunction(W(20), 0)
    xrts = find_zeros(F, -1, 21)
    @test length(xrts) == 20
    @test all(azero.((F,), xrts))
    @test F.n <= 20_000

    F = CallableFunction(Wi(6), 0)
    xrts = find_zeros(F, -1, 7)
    @test length(xrts) == 6
    @test all(azero.((F,), xrts))
    @test F.n <= 10_000

    ## Harder ones
    f1(x) = 2*exp(.5*x) *(sin(5*x) + sin(101*x))
    tiger_tail(x) = f1(x) - round(f1(x)) # (-2,1) then filter

    f2(x) = (x-0.5)*(x-0.5001)*(x-4)*(x-4.05)*(x-9.3) # (0,10)
    f3(x) = (x-3)^2*(x-4)^2 # (0,5)
    delta = .001
    f4(x) = (x - 0.5)^3 * (x - (0.5+delta)) * (x - 1)
    f5(x) = (x - 0.5)^3 * (x - (0.5+delta))^3 * (x-4)*(x-(4+delta))*(x-4.2)^2
    M(n) =  x -> prod((x - (0.5 - (1/10)^i)) for i in 1:n)
    f6 = M(4)
    f7 = M(5) # too much

    xrts = find_zeros(tiger_tail, -2.0, 1.0)
    xrts = filter(u -> -1/4 < tiger_tail(u) < 1/4, xrts)
    @test 345 <= length(xrts) <= 345
    @test all(azero.((tiger_tail,), xrts))

    xrts = find_zeros(f2, 0.0, 10.0)
    @test length(xrts) == 5
    @test all(azero.((f2,), xrts))

    xrts = find_zeros(f3, 0.0, 5.0)
    @test length(xrts) == 2
    @test all(azero.((f3,), xrts))

    xrts = find_zeros(f4, 0.0, 1.5)
    @test length(xrts) == 3
    @test all(azero.((f4,), xrts))

    xrts = find_zeros(f5, 0.0, 5.0)
    @test length(xrts) >= 3          # too hard to get 5 w/o luck, as with no_pts=21/k=4
    @test all(azero.((f5,), xrts))

    xrts = find_zeros(f6, 0.0, 10.0)
    @test length(xrts) == 4
    @test all(azero.((f6,), xrts))

    xrts = find_zeros(f7, 0.0, 10.0)  # too sensitive to interval
    @test length(xrts) >= 3           # should be 5
    @test all(azero.((f7,), xrts))
end


text:
# Define some points and data in 2D
points = permutedims([0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0], (2,1))
points_adjoint = [0.0 0.0; 0.0 1.0; 1.0 0.0; 1.0 1.0]'
data = [0.0; 0.5; 0.5; 1.0]

@testset "NearestNeighbor" begin

    @testset "Evaluation" for p in (points, points_adjoint)

        itp = interpolate(NearestNeighbor(), p, data)

        # Check that we get back the original data at the sample points
        ev = evaluate(itp, points)
        @test ev ≈ data

        # Check for a different point, should be closest to the last point
        ev = evaluate(itp, [2.0; 2.0])
        @test ev[1] ≈ data[end]
    end
end


text:
# simple vertex shader (GLSL source) for simple graphics pipeline

vert_source="""

#define LIGHTING_ENABLED        arg(LIGHTING_ENABLED)
#define COLOR_ATTRIBUTE_ENABLED arg(COLOR_ATTRIBUTE_ENABLED)

uniform mat4 u_modelview_matrix;
uniform mat4 u_projection_matrix;
uniform vec4 u_color;

attribute  vec4 a_position;

#if LIGHTING_ENABLED
attribute  vec3 a_normal;
#endif

#if COLOR_ATTRIBUTE_ENABLED
attribute vec4 a_color;
#endif

#if LIGHTING_ENABLED
uniform mat3 u_normal_matrix;
uniform vec3 u_light_position;
varying vec3 v_normal;
varying vec3 v_light_dir;
varying vec3 v_eye_vec;
#endif

#if COLOR_ATTRIBUTE_ENABLED
varying vec4 v_color;
#endif

void main()
{
	vec4 eye_pos= u_modelview_matrix * a_position;

#if LIGHTING_ENABLED
	v_normal = u_normal_matrix * a_normal;
	vec3 vVertex = vec3(u_modelview_matrix * a_position);
	v_light_dir  = normalize(u_light_position - vVertex);
	v_eye_vec    = normalize(-vVertex);
#endif

#if COLOR_ATTRIBUTE_ENABLED
	v_color=a_color;
#endif

	gl_Position = u_projection_matrix * eye_pos;
}
"""

# simple fragment shader (GLSL source) for simple graphics pipeline

frag_source="""

#define LIGHTING_ENABLED        arg(LIGHTING_ENABLED)
#define COLOR_ATTRIBUTE_ENABLED arg(COLOR_ATTRIBUTE_ENABLED)

uniform vec4 u_color;

#if LIGHTING_ENABLED
varying vec3 v_normal;
varying vec3 v_light_dir;
varying vec3 v_eye_vec;
#endif

#if COLOR_ATTRIBUTE_ENABLED
varying vec4 v_color;
#endif

void main()
{
	vec4 frag_color=u_color;

  #if LIGHTING_ENABLED
	vec3 N = normalize(v_normal   );
	vec3 L = normalize(v_light_dir);
	vec3 E = normalize(v_eye_vec  );

	vec4  u_material_ambient  = vec4(0.2,0.2,0.2,1.0);
	vec4  u_material_diffuse  = vec4(0.8,0.8,0.8,1.0);
	vec4  u_material_specular = vec4(0.1,0.1,0.1,1.0);
	float u_material_shininess=100.0;

	if(gl_FrontFacing)
	{
		frag_color = u_material_ambient;
		float NdotL = abs(dot(N,L));
		if (NdotL>0.0)
			{
			vec3 R = reflect(-L, N);
			float NdotHV = abs(dot(R, E));
			frag_color += u_material_diffuse * NdotL;
			frag_color += u_material_specular * pow(NdotHV,u_material_shininess);
		}
	}
	else
	{
		frag_color = u_material_ambient;
		float NdotL = abs(dot(-N,L));
		if (NdotL>0.0);
		{
			vec3 R = reflect(-L, -N);
			float NdotHV=abs(dot(R, E));
			frag_color += u_material_diffuse * NdotL;
			frag_color += u_material_specular * pow(NdotHV,u_material_shininess);
		}
	}
#endif

#if COLOR_ATTRIBUTE_ENABLED
	frag_color =v_color;
#endif

	gl_FragColor = frag_color;
}
"""

"""
	GLPhongShader(lighting_enabled,color_attribute_enabled)

Enable Phong's color shading and lighting shading.
Vertex and fragment shaders are given in GLSL.
"""
function GLPhongShader(lighting_enabled,color_attribute_enabled)

	v=vert_source
	f=frag_source

	v=replace(v, "arg(LIGHTING_ENABLED)"=>lighting_enabled ? "1" : "0")
	f=replace(f, "arg(LIGHTING_ENABLED)"=>lighting_enabled ? "1" : "0")

	v=replace(v, "arg(COLOR_ATTRIBUTE_ENABLED)"=>color_attribute_enabled ? "1" : "0")
	f=replace(f, "arg(COLOR_ATTRIBUTE_ENABLED)"=>color_attribute_enabled ? "1" : "0")

	# this is needed for #version 330
	# v=replace(v, "attribute"=>"in")
	# f=replace(f, "attribute"=>"in")

	# v=replace(v, "varying"=>"out")
	# f=replace(f, "varying"=>"out")

	# v=string("#version 120\n",v)
	# f=string("#version 120\n",f)

	return GLShader(v,f)
end


text:
"""
    to_symbol(input) -> String

Cleans up `input` string and converts it into a symbol, needed so that channel
names work with AxisArrays.
"""
function to_symbol(input::String)
    fixed = replace(input, r"[^\w\ \-\_]"=>"")
    fixed = replace(fixed, r"[\ \-\_]+"=>"_")
    Symbol(replace(fixed, r"^[\d]"=>s"_\g<0>"))
end

"""
    myendian() -> UInt16

Returns the endianness of the host machine
"""
function myendian()
    if ENDIAN_BOM == 0x04030201
        return 0x4949
    elseif ENDIAN_BOM == 0x01020304
        return 0x4d4d
    end
end

"""
    check_bswap(io::Union{Stream, IOStream})

Check endianness of TIFF file to see if we need to swap bytes
"""
function check_bswap(io::Union{Stream, IOStream})
    seekstart(io)
    endianness = read(io, UInt16)
    # check if we need to swap byte order
    need_bswap = endianness != myendian()

    tiff_version = Array{UInt8}(undef, 2)
    read!(io, tiff_version)
    (tiff_version != [0x2a, 0x00] && tiff_version != [0x00, 0x2a]) && error("Big-TIFF files aren't supported yet")

    return need_bswap
end

"""
    extract_filename(io) -> String

Extract the name of the file backing a stream
"""
extract_filename(io::IOStream) = split(io.name, " ")[2][1:end-1]
extract_filename(io::Stream) = io.filename

"""Corresponding Julian types for OME-XML types"""
type_mapping = Dict(
    "uint8" => (UInt8, N0f8),
    "uint16" => (UInt16, N0f16),
    "uint32" => (UInt32, N0f32),
    "float" => (Float32, Float32),
    "double" => (Float64, Float64),
    "int8" => (Int8, N0f8)
)

function getstream(fmt, io, name)
    # adapted from https://github.com/JuliaStats/RDatasets.jl/pull/119/
    if isdefined(FileIO, :action)
        # FileIO >= 1.6
        return Stream{fmt}(io, name)
    else
        # FileIO < 1.6
        return Stream(fmt, io, name)
    end
end

getstream(fmt, io::IOBuffer) = getstream(fmt, io, "")
getstream(fmt, io::IOStream) = getstream(fmt, io, extract_filename(io))
# assume OMETIFF if no format given
getstream(io) = getstream(format"OMETIFF", io)


text:
platforms_and_sources = Dict(
    Platform("x86_64", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.11.4/agnostic/x64/nccl_2.11.4-1+cuda10.2_x86_64.txz",
                      "d2f59b9385cb2026d9a8d4dee607f9aa93356761257fd57487f81de8563c8da1")],
    Platform("powerpc64le", "linux") => [
        ArchiveSource("https://developer.nvidia.com/compute/machine-learning/nccl/secure/2.11.4/agnostic/ppc64le/nccl_2.11.4-1+cuda10.2_ppc64le.txz",
                      "78b9b08d0d943239706c82db12c23e098ee2094108462055e03d74bc1a364601")],
)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
module MeshTests

using Test
using Jalerkin

p = Point{2}(0,0)
@test p == [0,0]

c = Cell{2}(1,2)
@test c == [1,2]

# 1D mesh
m = Mesh1D(0, 1, 10)

@test get_cell_vector(m) == [[i, i+1] for i=1:10]
@test get_coord_vector(m) ≈ [[0.1*i] for i=0:10]

# Quadrangular 2D mesh
C = Cell{4}
P = Point{2}
m = Mesh{C,P}( [C(1,2,3,4), C(3,4,5,6)],
        [P(0,0),P(0,1), P(0.5,0),P(0.5,1), P(1,0),P(1,1)])

end #module

max_stars_repo_name: rrgalvan/Jalerkin

max_issues_repo_name: rrgalvan/Jalerkin

max_forks_repo_name: rrgalvan/Jalerkin


text:
struct Button{T} <: AbstractWidget{Bool}
    content::Observable{T}
    value::Observable{Bool}
    attributes::Dict{Symbol, Any}
end

function Button(content; kw...)
    return Button(
        Observable(content), Observable(false), Dict{Symbol, Any}(kw)
    )
end

struct TextField <: AbstractWidget{String}
    value::Observable{String}
    attributes::Dict{Symbol, Any}
end

function TextField(value::String; kw...)
    TextField(Observable(value), Dict{Symbol, Any}(kw))
end

struct NumberInput <: AbstractWidget{Float64}
    value::Observable{Float64}
    attributes::Dict{Symbol, Any}
end

function NumberInput(value::Float64; kw...)
    NumberInput(Observable(value), Dict{Symbol, Any}(kw))
end

struct Slider{T <: AbstractRange, ET} <: AbstractWidget{T}
    range::Observable{T}
    value::Observable{ET}
    attributes::Dict{Symbol, Any}
end

function Slider(range::T, value = first(range); kw...) where T <: AbstractRange
    Slider{T, eltype(range)}(
        to_node(range),
        to_node(value),
        Dict{Symbol, Any}(kw)
    )
end

@enum Orientation vertical horizontal

struct RangeSlider{T <: AbstractRange, ET <: AbstractArray} <: AbstractWidget{T}
    attributes::Dict{Symbol, Any}
    range::Observable{T}
    value::Observable{ET}
    connect::Observable{Bool}
    orientation::Observable{Orientation}
    tooltips::Observable{Bool}
    ticks::Observable{Dict{String, Any}}
end

function RangeSlider(range::T; value = [first(range)], kw...) where T <: AbstractRange
    RangeSlider{T, typeof(value)}(
        Dict{Symbol, Any}(kw),
        Observable(range),
        Observable(value),
        Observable(true),
        Observable(horizontal),
        Observable(false),
        Observable(Dict{String, Any}()),
    )
end

struct Checkbox <: AbstractWidget{Bool}
    value::Observable{Bool}
    attributes::Dict{Symbol, Any}
end

function Checkbox(value::Bool; kw...)
    return Checkbox(Observable(value), Dict{Symbol, Any}(kw))
end


text:
mutable struct EM
    probA::Float32
    probB::Float32
    probC::Float32
end

EM() = EM(0.5, 0.5, 0.5)
"""
    E_step(self::EM, iter::Int32) -> Float32

返回一个浮点数μ^{i+1}
"""
function E_step(self::EM, data::Bool)
    probFromB = self.probA * self.probB^data * (1-self.probB)^(1-data)
    probSum = probFromB +
                (1-self.probA) * self.probC^data * (1-self.probC)^(1-data)
    return probFromB / probSum
end

function M_step(self::EM, datas::Array{Bool, 1}; maxIter=4)
    for i in 1:maxIter
    muProbs = [E_step(self, data) for data in datas]
    self.probA = 1/length(datas) * sum(muProbs)
    self.probB = sum([mu*data for (mu, data) in zip(muProbs, datas)]) / sum(muProbs)
    self.probC = sum([(1-mu)*data for (mu, data) in zip(muProbs, datas)]) / sum(map(mu -> 1-mu, muProbs))
    println("第$(i)轮结果:")
    println("A的概率为$(self.probA)")
    println("B的概率为$(self.probB)")
    println("C的概率为$(self.probC)\n")
    end
end

function EMbase()
    boolArray(x) = map(Bool, x)
    datas = boolArray([1, 1, 0, 1, 0, 0, 1, 0, 1, 1])
    #M_step(EM(), datas)
    M_step(EM(0.4, 0.6, 0.7), datas)
end

mutable struct GaussModel
    alpha::Float32
    mu::Float32
    sigma::Float32
end

gaussProb(gm::GaussModel, x::AbstractFloat) = (1/sqrt(gm.sigma)) * exp(-(x-gm.mu)^2 / (2gm.sigma^2))

gausses = [GaussModel(0.3, 1.2, 2), GaussModel(0.5, 1., 0.4), GaussModel(0.2, 0.8, 1.3)]

using Distributions
"""
    GMMgenerator() -> Float32
产生基于gausses的数据
"""
function GMMgenerator()
    randNum, accumulate = rand(), 0.
    for (i, gauss) in enumerate(gausses)
        if randNum > accumulate
            return rand(Normal(gauss.mu, gauss.sigma))
        else
            accumulate += gauss.alpha
        end
    end
end

function GMM(_K=3, maxIter=10)
    _K = length(gausses)
    println("生成数据的Gauss分布为：")
    println(gausses, "\n")
    # Init
    datas = [GMMgenerator() for _ in 1:30000]
    alphas = [rand() for _ in 1:_K]; alphas ./= sum(alphas)
    fakeGausses = [GaussModel(alphas[i], rand(), rand()) for i in 1:_K]

    # 判断循环终止条件
    for _ in 1:maxIter
    println(fakeGausses)

    # E_step
    response = zeros(_K, length(datas))
    for j in 1:length(datas)
        dominators = zeros(1, _K)
        for k in 1:_K
            response[k, j] = gaussProb(fakeGausses[k], datas[j])
        response[:, j] ./= sum(response[:, j])
        end
    end

    # M_step
    for k in 1:_K
        old_mu = fakeGausses[k].mu
        fakeGausses[k].mu   = sum( map(x -> x[1]*x[2] ,zip(response[k, :], datas)) ) / sum(response[k, :])
        fakeGausses[k].sigma= sum( map(x -> x[1]*(x[2]-old_mu)^2, zip(response[k, :], datas)) ) / sum(response[k, :])
        fakeGausses[k].alpha= sum( response[k, :] ) / length(datas)
    end

    end

    println("GaussMixtureModel运行完毕")
    println("共$(_K)个高斯分布")
    for k in 1:_K
        println("第$(k)个高斯分布参数：")
        println("权重：$(fakeGausses[k].alpha), 均值：$(fakeGausses[k].mu), 方差：$(fakeGausses[k].sigma)")
    end

end

function main()
    GMM()
end

main()

max_stars_repo_name: gajanlee/toys

max_issues_repo_name: gajanlee/toys

max_forks_repo_name: gajanlee/toys


text:
immutable Grouped{T,I,L}
    src::T
    group_indices::I
    groupbys::Vector{Symbol}
    group_levels::L
    metadata::Dict{Symbol, Any}
end


text:
# k-NN Searcher
# =============

"""
A set of k-nearest neighboring cells.

Fields
------

- `indexes`: cell indexes (k x #queries).
- `hammingdistances`: Hamming distances (k x #queries).
"""
struct NearestCells
    indexes::Matrix{Int}
    hammingdistances::Matrix{Int16}
    # some scores?

    NearestCells(k, n) = new(zeros(Int, k, n), zeros(Int16, k, n))
end

"""
    similarities(ncs::NearestCells, index::CellIndex) -> Matrix{Float32}

Estimate cosine similarities between queries and neighbors from the Hamming distance.
"""
function similarities(ncs::NearestCells, index::CellIndex)
    T = nbits(index) * length(index.lshashes)
    return cos.(ncs.hammingdistances .* (Float32(pi) ./ T))
end

"""
    findneighbors(
        k::Integer,
        counts::AbstractMatrix,
        featurenames::AbstractVector{String},
        index::CellIndex;
        inferstats::Symbol=:both
    ) -> NearestCells

Find `k`-nearest neighboring cells from `index`.

The `inferstats` parameter specifies the way to infer feature (gene)
statistics. If it is `:both` (default), they are inferred both from the query
and the database cells. If it is `:query`, they are inferred only from the
query cells.  If it is `:database`, they are inferred only from the database
cells. No other values are allowed.
"""
function findneighbors(
        k::Integer,
        counts::AbstractMatrix,
        featurenames::AbstractVector{String},
        index::CellIndex;
        inferstats::Symbol=:both)
    return findneighbors(k, ExpressionMatrix(counts, featurenames), index; inferstats=inferstats)
end

function findneighbors(k::Integer, Y::ExpressionMatrix, index::CellIndex; inferstats::Symbol=:both)
    if k < 0
        throw(ArgumentError("negative k"))
    end
    if inferstats ∉ (:both, :query, :database)
        throw(ArgumentError("invalid value for the inferstats parameter"))
    end
    n = size(Y, 2)
    L = length(index.lshashes)
    T = bitvectype(first(index.lshashes))
    @assert L ≥ 1
    rtstats = index.rtstats
    tic!(rtstats)
    # preprocess
    X = preprocess(index.preproc, Y, inferstats)
    # allocate temporary memories
    neighbors = Matrix{Int}(undef, k * L, n)
    neighbors_tmp = Matrix{Int}(undef, k, n)
    Z = Matrix{T}(undef, L, n)
    Z_tmp = Vector{T}(undef, n)
    for l in 1:L
        rtstats.preproc_time += tic!(rtstats)
        findneighbors!(neighbors_tmp, Z_tmp, X, index.lshashes[l])
        neighbors[k*(l-1)+1:k*l,:] = neighbors_tmp
        Z[l,:] = Z_tmp
        rtstats.search_time += tic!(rtstats)
    end
    @assert all(neighbors .> 0)
    # Rank neighbor cells.
    ncs = rankcells!(neighbors, Z, index.lshashes)
    rtstats.rank_time += tic!(rtstats)
    return ncs
end

function findneighbors(k::Integer, I::AbstractVector{Int}, index::CellIndex)
    if k < 0
        throw(ArgumentError("negative k"))
    end
    n = length(I)
    L = length(index.lshashes)
    @assert L ≥ 1
    T = CellFishing.bitvectype(first(index.lshashes))
    neighbors = zeros(Int, k * L, n)
    Z = Matrix{T}(undef, L, n)
    nns = HammingIndexes.NearestNeighbors(k)
    for l in 1:L
        lshash = index.lshashes[l]
        for j in 1:n
            h = lshash.hammingindex[I[j]]
            Z[l,j] = h
            HammingIndexes.findknn!(nns, h, lshash.hammingindex)
            for i in 1:k
                neighbors[k*(l-1)+i,j] = nns.indexes[i]
            end
        end
    end
    @assert all(neighbors .> 0)
    return CellFishing.rankcells!(neighbors, Z, index.lshashes)
end

function rankcells!(neighbors::Matrix{Int}, Z::Matrix{T}, lshashes::Vector{LSHash{T}}) where T
    L, N = size(Z)
    k = div(size(neighbors, 1), L)
    @assert size(neighbors) == (k * L, N)
    ncs = NearestCells(k, N)
    fill!(ncs.hammingdistances, L * bitsof(T) + 1)
    for j in 1:N
        last = 0
        sort!(@view neighbors[:,j])
        for i in 1:k*L
            x = neighbors[i,j]
            if x == last
                continue
            end
            d = 0
            for l in 1:L
                d += count_ones(Z[l,j] ⊻ lshashes[l].hammingindex[x])
            end
            if ncs.indexes[k,j] == 0 || d < ncs.hammingdistances[k,j]
                ncs.indexes[k,j] = x
                ncs.hammingdistances[k,j] = d
                l = k - 1
                while l ≥ 1 && ncs.hammingdistances[l,j] > ncs.hammingdistances[l+1,j]
                    ncs.indexes[l,j], ncs.indexes[l+1,j] = ncs.indexes[l+1,j], ncs.indexes[l,j]
                    ncs.hammingdistances[l,j], ncs.hammingdistances[l+1,j] = ncs.hammingdistances[l+1,j], ncs.hammingdistances[l,j]
                    l -= 1
                end
            end
            last = x
        end
        #@assert issorted(ncs.hammingdistances[:,j])
    end
    return ncs
end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

                     @variable(m, -1 <= x[i=1:3] <= 1)

                     @variable(m, -12.75118444054646 <= q <= 9.039664285972906)

                     add_NL_constraint(m, :(swish(0.715796138346545 + -0.1516695982843861*swish(-0.317423658088571 + -0.8843064794993589*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + -0.0820632819191216*swish(-0.8960501996212993 + -0.7682728024607406*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + -0.0571486344530765*swish(0.054021626349418206 + -0.05819989275612114*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) + swish(-0.8464220616292248 + -0.9331432460947604*swish(-0.317423658088571 + -0.8843064794993589*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.425513343659472*swish(-0.8960501996212993 + -0.7682728024607406*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.1421259709107252*swish(0.054021626349418206 + -0.05819989275612114*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) + swish(-0.5205460270098605 + 0.13812073353255405*swish(-0.317423658088571 + -0.8843064794993589*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.6219431368805965*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + -0.24757066957523488*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.7749706249828683*swish(-0.8960501996212993 + -0.7682728024607406*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + -0.3378399703859083*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.4153867033969405*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3]))) + 0.5007112790542019*swish(0.054021626349418206 + -0.05819989275612114*swish(0.5491785730840042 + -0.947600132880416*$(x[1]) + -0.5975167017748926*$(x[2]) + -0.008573066442256305*$(x[3])) + 0.8875139231924525*swish(0.2436152968029921 + 0.15213749662876674*$(x[1]) + -0.7774236630629989*$(x[2]) + 0.8200325525995154*$(x[3])) + 0.11280184637312685*swish(-0.7915430864078541 + -0.04671642019396449*$(x[1]) + 0.6799142288942424*$(x[2]) + 0.9251799998319439*$(x[3])))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
# Hamiltonian operator
h(φ::Float64, Δφ::Float64, v::Float64) = -Δφ/2 + v*φ

function h!(dst::AbstractArray, φ::AbstractArray, Δφ::AbstractArray, v::AbstractArray)
    @assert length(dst) == length(φ) == length(Δφ) == length(v)
    for i in 1:length(dst)
        dst[i] = h(φ[i],Δφ[i],v[i])
    end
end

function h!(dst::AbstractMatrix, φ::AbstractMatrix, Δφ::AbstractMatrix, v::AbstractArray)
    @assert size(dst) == size(φ) == size(Δφ)
    @assert size(dst,1) == length(v)
    for i in 1:size(dst,2)
        for j in 1:size(dst,1)
            dst[j,i] = h(φ[j,i],Δφ[j,i],v[j])
        end
    end
end

function h(φ::AbstractArray, Δφ::AbstractArray, v::AbstractArray)
    dst = similar(φ)
    h!(dst,φ,Δφ,v)
    dst
end

function h!(dst::SparseMatrixCSC, φ::SparseMatrixCSC, Δφ::SparseMatrixCSC, v::AbstractArray)
    @assert size(dst) == size(φ) == size(Δφ)
    @assert size(dst,1) == length(v)
    for i in 1:length(φ.nzval)
        j = φ.rowval[i]
        dst.nzval[i] = h(φ.nzval[i],Δφ.nzval[i],v[j])
    end
end

function h(φ::SparseMatrixCSC, Δφ::SparseMatrixCSC, v::AbstractArray)
    dst = SparseMatrixCSC(φ.m,φ.n,φ.colptr,φ.rowval,similar(φ.nzval))
    h!(dst,φ,Δφ,v)
    dst
end

max_stars_repo_name: johncsnyder/Element


text:
Random.seed!(0)
@testset "SparseVector: Constructors B=$B T=$T" for B in [1, 2, 4, 8],
                                                    T in [Int32, Float32,
                                                     Complex{Float32},
                                                     SMatrix{2,2,Float32,2 * 2}]

    len = 10 * B
    n = bspzeros(Val(B), T, len)
    @test length(n) == len
    @test eltype(n) == T

    x = bsprand(Val(B), T, len, 0.3)
    @test length(x) == len
    @test eltype(x) == T
    if !(T <: Integer)
        y = bsprandn(Val(B), T, len, 0.3)
        @test length(y) == len
        @test eltype(y) == T
    end

    xf = Vector(x)
    @test eltype(xf) == T
    @test length(xf) == length(x)
    x′ = blocksparse(Val(B), xf)
    @test length(x′) == length(xf)
    @test eltype(x′) == T
    @test x′ == x

    xs = SparseVector(x)
    @test eltype(xs) == T
    @test length(xs) == len
    x′ = blocksparse(Val(B), xs)
    @test eltype(x′) == T
    @test x′ == x

    zf = rand(T, len)::Vector{T}
    z = blocksparse(Val(B), zf)
    @test eltype(z) == T
    @test size(z) == size(zf)
    zf′ = Array(z)
    @test eltype(zf′) == T
    @test zf′ == zf

    zs = sprand(T, len, 0.3)::SparseVector{T}
    z = blocksparse(Val(B), zs)
    @test eltype(z) == T
    @test size(z) == size(zs)
    zs′ = SparseVector(z)
    @test eltype(zs′) == T
    @test zs′ == zs
end

Random.seed!(0)
@testset "SparseVector: Collection B=$B T=$T" for B in [1, 2, 4, 8],
                                                  T in [Int32, Float32,
                                                   Complex{Float32},
                                                   SMatrix{2,2,Float32,2 * 2}]

    len = 10 * B
    x = bsprand(Val(B), T, len, 0.3)
    xs = SparseVector(x)
    @test x == xs
end

Random.seed!(0)
@testset "SparseVector: Vector space B=$B T=$T" for B in [1, 2, 4, 8],
                                                    T in [Int32, Float32,
                                                     Complex{Float32},
                                                     SMatrix{2,2,Float32,2 * 2}]

    len = 10 * B

    n = bspzeros(Val(B), T, len)
    x = bsprand(Val(B), T, len, 0.3)
    y = bsprand(Val(B), T, len, 0.3)
    z = bsprand(Val(B), T, len, 0.3)
    a = rand(eltype(T))
    b = rand(eltype(T))

    @test x == x
    @test y ≠ x

    @test +x == x
    @test -(-x) == x
    @test n + x == x
    @test x + n == x
    @test n - x == -x
    @test x - n == x
    @test x + y == y + x
    @test (x + y) + z == x + (y + z)

    @test zero(a) * x == n
    @test one(a) * x == x
    @test (-one(a)) * x == -x
    @test a * x == x * a
    @test a * (b * x) ≈ (a * b) * x

    if T <: Union{AbstractFloat,Complex{<:AbstractFloat}}
        @test inv(a) * (a * x) ≈ x
        @test x / a ≈ x * inv(a)
        @test a \ x ≈ inv(a) * x
    end

    @test (a + b) * x ≈ a * x + b * x
    @test a * (x + y) ≈ a * x + a * y
end


text:
module HPD

using Base: @kwdef
using Julog
using AutoHashEquals
using Reexport

include(joinpath("parser", "parser.jl"))

@reexport using .Parser

end


text:
module CMBLensing

using Adapt
using Base.Broadcast: AbstractArrayStyle, ArrayStyle, Broadcasted, broadcasted,
    DefaultArrayStyle, preprocess_args, Style, result_style
using Base.Iterators: flatten, product, repeated, cycle, countfrom
using Base.Threads
using Base: @kwdef, @propagate_inbounds, Bottom, OneTo, showarg, show_datatype,
    show_default, show_vector, typed_vcat
using Combinatorics
using DataStructures
using DelimitedFiles
using Distributed: pmap, nworkers, myid, workers, addprocs, @everywhere, remotecall_wait, @spawnat
using FileIO
using FFTW
using InteractiveUtils
using IterTools: flagfirst
using JLD2: jldopen, JLDWriteSession
import KahanSummation
using Loess
using LinearAlgebra
using LinearAlgebra: diagzero, matprod, promote_op
using MacroTools: @capture, combinedef, isdef, isexpr, postwalk, splitdef
using Match
using Markdown
using Measurements
using Memoization
using OptimKit
using Pkg
using Printf
using ProgressMeter
using Random
using Random: seed!, AbstractRNG
using Roots
using Requires
using Setfield
using SparseArrays
using StaticArrays: @SMatrix, @SVector, SMatrix, StaticArray, StaticArrayStyle,
    StaticMatrix, StaticVector, SVector, SArray, SizedArray
using Statistics
using StatsBase
using UnPack
using Zygote
using Zygote: unbroadcast, Numeric, @adjoint, @nograd

import Adapt: adapt_structure
import Base: +, -, *, \, /, ^, ~, ≈, <, <=, |, &, ==,
    abs, adjoint, all, any, axes, broadcast, broadcastable, BroadcastStyle, conj, copy, convert,
    copy, copyto!, eltype, eps, fill!, getindex, getproperty, hash, hcat, hvcat, inv, isfinite,
    iterate, keys, lastindex, length, literal_pow, mapreduce, materialize!,
    materialize, one, permutedims, print_array, promote, promote_rule,
    promote_rule, promote_type, propertynames, real, setindex!, setproperty!, show,
    show_datatype, show_vector, similar, size, sqrt, string, sum, summary,
    transpose, zero
import Base.Broadcast: instantiate, preprocess
import LinearAlgebra: checksquare, diag, dot, isnan, ldiv!, logdet, mul!, norm,
    pinv, StructuredMatrixStyle, structured_broadcast_alloc, tr
import Measurements: ±
import Statistics: std

export
    @BandpowerParamOp, @ismain, @namedtuple, @repeated, @unpack, animate,
    argmaxf_lnP, BandPassOp, BaseDataSet, batch, batchindex, batchsize, beamCℓs, cache,
    CachedLenseFlow, camb, cov_to_Cℓ, cpu, Cℓ_2D, Cℓ_to_Cov, DataSet, DerivBasis,
    diag, Diagonal, DiagOp, dot, EBFourier, EBMap, expnorm, Field, FieldArray, fieldinfo,
    FieldMatrix, FieldOrOpArray, FieldOrOpMatrix, FieldOrOpRowVector,
    FieldOrOpVector, FieldRowVector, FieldTuple, FieldVector, FieldVector,
    firsthalf, fixed_white_noise, Flat, FlatEB, FlatEBFourier, FlatEBMap, FlatField,
    FlatFieldFourier, FlatFieldMap, FlatFourier, FlatIEBCov, FlatIEBFourier, FlatIEBMap,
    FlatIQUFourier, FlatIQUMap, FlatMap, FlatQU, FlatQUFourier, FlatQUMap, FlatS0,
    FlatS02, FlatS2, FlatS2Fourier, FlatS2Map, Fourier, fourier∂, FuncOp, get_Cℓ,
    get_Cℓ, get_Dℓ, get_αℓⁿCℓ, get_ρℓ, get_ℓ⁴Cℓ, gradhess, gradient, HighPass,
    IEBFourier, IEBMap, InterpolatedCℓs, IQUFourier, IQUMap,
    lasthalf, LazyBinaryOp, LenseBasis, LenseFlow, LinOp, lnP, load_camb_Cℓs,
    load_chains, load_sim, LowPass, make_mask, Map, MAP_joint, MAP_marg,
    map∂, mean_std_and_errors, MidPass, mix, nan2zero, new_dataset, noiseCℓs,
    OuterProdOp, ParamDependentOp, pixwin, PowerLens, QUFourier, QUMap, resimulate!,
    resimulate, RK4Solver, S0, S02, S2, sample_joint, seed_for_storage!, shiftℓ,
    simulate, SymmetricFuncOp, symplectic_integrate, Taylens, toCℓ, toDℓ,
    tuple_adjoint, ud_grade, unbatch, unmix, Ð, Ł, δfϕ_δf̃ϕ, ℓ², ℓ⁴, ∇, ∇², ∇¹, ∇ᵢ,
    ∇⁰, ∇ⁱ, ∇₀, ∇₁, ⋅, ⨳

# generic stuff
include("util.jl")
include("util_fft.jl")
include("numerical_algorithms.jl")
include("generic.jl")
include("cls.jl")
include("field_tuples.jl")
include("field_vectors.jl")
include("specialops.jl")
include("flowops.jl")

# lensing
include("lenseflow.jl")
include("powerlens.jl")

# flat-sky maps
include("flat_fftgrid.jl")
include("flat_s0.jl")
include("flat_s2.jl")
include("flat_s0s2.jl")
include("flat_generic.jl")
include("flat_batch.jl")
include("masking.jl")
include("taylens.jl")
include("bilinearlens.jl")

# plotting
function animate end
@init @require PyPlot="d330b81b-6aea-500a-939a-2ce795aea3ee" include("plotting.jl")

# sampling and maximizing the posteriors
include("dataset.jl")
include("posterior.jl")
include("maximization.jl")
include("sampling.jl")
include("chains.jl")

# other estimates
include("quadratic_estimate.jl")

# curved-sky
include("healpix.jl")

include("autodiff.jl")

# gpu
is_gpu_backed(x) = false
@init @require CUDA="052768ef-5323-5732-b1bb-66c8b64840ba" include("gpu.jl")

# misc init
# see https://github.com/timholy/ProgressMeter.jl/issues/71 and links therein
@init if ProgressMeter.@isdefined ijulia_behavior
    ProgressMeter.ijulia_behavior(:clear)
end

end


text:
using LinearAlgebra

function atom_under_mouse(E::Ensamble)
    obj, _ = mouse_selection(E.scene)
    if obj !== nothing
        for i = eachindex(E.atoms)
            if obj === E.atoms[i].plot
                return E.atoms[i]
            end
        end
    else
        return nothing
    end
end

function closest_atom(E::Ensamble, xyz::Point3f0)
    dist = zeros(Float32, length(E.atoms))
    for i = eachindex(dist)
        dist[i] = √sum((E.atoms[i].center .- xyz).^2)
    end
    return E.atoms[findmin(dist)[2]]
end

function create_atom_under_mouse(E::Ensamble)
    ori, dir = Makie.view_ray(E.scene)
    if length(E.atoms) > 0
        A = closest_atom(E, Point3f0(ori))
        xyz = ori .+ dir * norm(A.center .- ori)
        CreateAtom!(E, "N", xyz...)
    else
        xyz = ori .- dir * 10
        CreateAtom!(E, "N", xyz...)
    end
end

function respond_mouse_event(E::Ensamble, M::SelectionMode, event)
    # Event involving left mouse
    if event.button == Mouse.left
        # On click
        if event.action == Mouse.press
            A = atom_under_mouse(E)
            if A === nothing
                clear_selection!(E)
                return false
            end
            select!(E, A)
            return false
        end
    end
    return false
end

function respond_mouse_event(E::Ensamble, M::InsertMode, event)
    # Event involving left mouse
    if event.button == Mouse.left
        # On click
        if event.action == Mouse.press
            A = atom_under_mouse(E)
            # If no atom under mouse, create one
            if A === nothing
                create_atom_under_mouse(E)
            end
            return true
        # On release
        else
            return false
        end
    end
    return false
end


text:
push!(LOAD_PATH,"../src/")

using Documenter, PDFIO

makedocs(
    format = :html,
    sitename = "PDFIO",
    pages = [
        "intro.md",
        "API Structure and Design" => "index.md"
    ]
)

deploydocs(
    repo   = "github.com/sambitdash/PDFIO.jl.git",
    target = "build",
    branch = "gh-pages",
    julia  = "0.6",
    deps   = nothing,
    make   = nothing
)


text:
# The most basic terminals have a set of 8 different colors

@inline function black(str)
    return "\u001b[30m$str\u001b[0m"
end

@inline function red(str)
    return "\u001b[31m$str\u001b[0m"
end

@inline function green(str)
    return "\u001b[32m$str\u001b[0m"
end

@inline function yellow(str)
    return "\u001b[33m$str\u001b[0m"
end

@inline function blue(str)
    return "\u001b[34m$str\u001b[0m"
end

@inline function magenta(str)
    return "\u001b[35m$str\u001b[0m"
end

@inline function cyan(str)
    return "\u001b[36m$str\u001b[0m"
end

@inline function white(str)
    return "\u001b[37m$str\u001b[0m"
end


text:
function overprint(str)
  print("\u1b[1F")
  print(str)
  print("\u1b[0K")
  println()
end

function print_info_dataset(data)
    println("Dataset name : ", data[1, :name])
    println("Number of vertices : ", data[1, :nodes])
    println("Number of edges (required/not-required) : ", data[1, :req_edges], " / ", data[1, :noreq_edges])
    println("Number of arcs (required/not-required) : ", data[1, :req_arcs], " / ", data[1, :noreq_arcs])
    println("Max capacity of vehicles : ", data[1, :capa])
    println("Depot(s) : ", string(data[1, :depots])[2:end-1])
    println("Number of vehicle(s) : ", data[1, :nb_vehicles])
end

function remove_multiple_service(tours, graph)

  for arc in 1:nrow(graph["req_arcs"])
    servicing_tours_indices = []
    for (i,tour) in enumerate(tours)
      if tour["X"]["req_arcs"][arc] == 1
        push!(servicing_tours_indices, i)
      end
    end
    if length(servicing_tours_indices) != 0
      (val,ind) = findmin([tours[i]["cost"] for i in servicing_tours_indices])
      for index in servicing_tours_indices
        if index != servicing_tours_indices[ind]
          tours[index]["X"]["req_arcs"][arc] = 0
          tours[index]["Y"]["req_arcs"][arc] += 1
          tours[index]["cost"] += (graph["req_arcs"][arc, :trav_cost] - graph["req_arcs"][arc, :serv_cost])
          tours[index]["capacity"] -= graph["req_arcs"][arc, :demand]
        end
      end
    end
  end

  for arc in 1:nrow(graph["req_edges"])
    servicing_tours_indices = []
    for (i,tour) in enumerate(tours)
      if (tour["X"]["req_edges"][arc]==1) | (tour["X"]["req_edges_swap"][arc]==1)
        push!(servicing_tours_indices, i)
      end
    end
    if length(servicing_tours_indices) != 0
      (val,ind) = findmin([tours[i]["cost"] for i in servicing_tours_indices])
      for index in servicing_tours_indices
        if index != servicing_tours_indices[ind]
          if tours[index]["X"]["req_edges"][arc] == 1
            tours[index]["X"]["req_edges"][arc] = 0
            tours[index]["Y"]["req_edges"][arc] += 1
          else
            tours[index]["X"]["req_edges_swap"][arc] = 0
            tours[index]["Y"]["req_edges_swap"][arc] += 1
          end
          tours[index]["cost"] += (graph["req_edges"][arc, :trav_cost] - graph["req_edges"][arc, :serv_cost])
          tours[index]["capacity"] -= graph["req_edges"][arc, :demand]
        end
      end
    end
  end

  return(tours)
end

max_stars_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_issues_repo_name: TheoGuyard/Groomer-Optimization-Problem

max_forks_repo_name: TheoGuyard/Groomer-Optimization-Problem


text:
"""
    module MaxwellOperators
"""
module MaxwellOperators

export Maxwell

files = (
    "1D/MaxwellOperators1D.jl",
    # "2D/MaxwellOperators2D.jl",
    # "3D/MaxwellOperators3D.jl"
    )

import ..PRINTED_COLOR_DARK
import ..PRINTED_COLOR_VARIABLE
import ..PRINTED_COLOR_INSTRUCTION

using ..VectorFields
using ..Simulations
using ..Dispersions
using LinearAlgebra
using SparseArrays

const I3 = sparse(I,3,3)

"""
    maxwell(simulation{1}; ky=0, kz=0) -> m
"""
struct Maxwell{N,TΣ,TSIM}
    A::SparseMatrixCSC{ComplexF64,Int}
    D²::SparseMatrixCSC{ComplexF64,Int}
    αεpFχ::SparseMatrixCSC{Complex{Float64},Int}
    curlcurl::SparseMatrixCSC{ComplexF64,Int}
    Σs::TΣ
    αε::SparseMatrixCSC{ComplexF64,Int}
    Fχ::SparseMatrixCSC{ComplexF64,Int}
    Fχs::Vector{SparseMatrixCSC{ComplexF64,Int}}
    dFχdψr::Array{ComplexF64,3}
    dFχdψi::Array{ComplexF64,3}
    dFχdω::Array{ComplexF64,3}
    dFχdϕ::Array{ComplexF64,3}
    simulation::TSIM
    kx::Float64
    ky::Float64
    kz::Float64
    ka::Float64
    kb::Float64
    kc::Float64
end

# load 1D, 2D, 3D
foreach(include,files)

# Pretty Printing
function Base.show(io::IO,::Maxwell{N}) where N
    print(io,N,"D ")
    printstyled(io,"Maxwell ",color=PRINTED_COLOR_DARK)
    print(io,"Operator ")
    printstyled(io,"(call w/ args ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ω",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,", [",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ωs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,",",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"ψs",color=PRINTED_COLOR_VARIABLE)
    printstyled(io,"]) -> ",color=PRINTED_COLOR_INSTRUCTION)
    printstyled(io,"A",color=PRINTED_COLOR_VARIABLE)
end

function Base.getproperty(m::Maxwell,sym::Symbol)
    if sym==:sim
        return getfield(m,:simulation)
    else
        return getfield(m,sym)
    end
end

function Base.propertynames(::Maxwell{1},private=false)
    if private
        return fieldnames(Maxwell)
    else
        (:simulation,:ky,:kz)
    end
end

end #module


text:
module TestWithSequential

using FoldsTapir
using Folds.Testing: test_with_sequential

executors = [TapirEx(), TapirEx(basesize = 1), TapirEx(basesize = 3)]

test_with_sequential(executors)

end


text:
@test first(evaluate5(A♣, A♡, A♢, K♠, K♣)) == 167    # Aces Full over Kings
@test first(evaluate5(A♣, A♡, A♢, Q♠, Q♣)) == 168    # Aces Full over Queens
@test first(evaluate5(A♣, A♡, A♢, J♠, J♣)) == 169    # Aces Full over Jacks
@test first(evaluate5(A♣, A♡, A♢, T♠, T♣)) == 170    # Aces Full over Tens
@test first(evaluate5(A♣, A♡, A♢, 9♠, 9♣)) == 171    # Aces Full over Nines
@test first(evaluate5(A♣, A♡, A♢, 8♠, 8♣)) == 172    # Aces Full over Eights
@test first(evaluate5(A♣, A♡, A♢, 7♠, 7♣)) == 173    # Aces Full over Sevens
@test first(evaluate5(A♣, A♡, A♢, 6♠, 6♣)) == 174    # Aces Full over Sixes
@test first(evaluate5(A♣, A♡, A♢, 5♠, 5♣)) == 175    # Aces Full over Fives
@test first(evaluate5(A♣, A♡, A♢, 4♠, 4♣)) == 176    # Aces Full over Fours
@test first(evaluate5(A♣, A♡, A♢, 3♠, 3♣)) == 177    # Aces Full over Treys
@test first(evaluate5(A♣, A♡, A♢, 2♠, 2♣)) == 178    # Aces Full over Deuces
@test first(evaluate5(K♣, K♡, K♢, A♠, A♣)) == 179    # Kings Full over Aces
@test first(evaluate5(K♣, K♡, K♢, Q♠, Q♣)) == 180    # Kings Full over Queens
@test first(evaluate5(K♣, K♡, K♢, J♠, J♣)) == 181    # Kings Full over Jacks
@test first(evaluate5(K♣, K♡, K♢, T♠, T♣)) == 182    # Kings Full over Tens
@test first(evaluate5(K♣, K♡, K♢, 9♠, 9♣)) == 183    # Kings Full over Nines
@test first(evaluate5(K♣, K♡, K♢, 8♠, 8♣)) == 184    # Kings Full over Eights
@test first(evaluate5(K♣, K♡, K♢, 7♠, 7♣)) == 185    # Kings Full over Sevens
@test first(evaluate5(K♣, K♡, K♢, 6♠, 6♣)) == 186    # Kings Full over Sixes
@test first(evaluate5(K♣, K♡, K♢, 5♠, 5♣)) == 187    # Kings Full over Fives
@test first(evaluate5(K♣, K♡, K♢, 4♠, 4♣)) == 188    # Kings Full over Fours
@test first(evaluate5(K♣, K♡, K♢, 3♠, 3♣)) == 189    # Kings Full over Treys
@test first(evaluate5(K♣, K♡, K♢, 2♠, 2♣)) == 190    # Kings Full over Deuces
@test first(evaluate5(Q♣, Q♡, Q♢, A♠, A♣)) == 191    # Queens Full over Aces
@test first(evaluate5(Q♣, Q♡, Q♢, K♠, K♣)) == 192    # Queens Full over Kings
@test first(evaluate5(Q♣, Q♡, Q♢, J♠, J♣)) == 193    # Queens Full over Jacks
@test first(evaluate5(Q♣, Q♡, Q♢, T♠, T♣)) == 194    # Queens Full over Tens
@test first(evaluate5(Q♣, Q♡, Q♢, 9♠, 9♣)) == 195    # Queens Full over Nines
@test first(evaluate5(Q♣, Q♡, Q♢, 8♠, 8♣)) == 196    # Queens Full over Eights
@test first(evaluate5(Q♣, Q♡, Q♢, 7♠, 7♣)) == 197    # Queens Full over Sevens
@test first(evaluate5(Q♣, Q♡, Q♢, 6♠, 6♣)) == 198    # Queens Full over Sixes
@test first(evaluate5(Q♣, Q♡, Q♢, 5♠, 5♣)) == 199    # Queens Full over Fives
@test first(evaluate5(Q♣, Q♡, Q♢, 4♠, 4♣)) == 200    # Queens Full over Fours
@test first(evaluate5(Q♣, Q♡, Q♢, 3♠, 3♣)) == 201    # Queens Full over Treys
@test first(evaluate5(Q♣, Q♡, Q♢, 2♠, 2♣)) == 202    # Queens Full over Deuces
@test first(evaluate5(J♣, J♡, J♢, A♠, A♣)) == 203    # Jacks Full over Aces
@test first(evaluate5(J♣, J♡, J♢, K♠, K♣)) == 204    # Jacks Full over Kings
@test first(evaluate5(J♣, J♡, J♢, Q♠, Q♣)) == 205    # Jacks Full over Queens
@test first(evaluate5(J♣, J♡, J♢, T♠, T♣)) == 206    # Jacks Full over Tens
@test first(evaluate5(J♣, J♡, J♢, 9♠, 9♣)) == 207    # Jacks Full over Nines
@test first(evaluate5(J♣, J♡, J♢, 8♠, 8♣)) == 208    # Jacks Full over Eights
@test first(evaluate5(J♣, J♡, J♢, 7♠, 7♣)) == 209    # Jacks Full over Sevens
@test first(evaluate5(J♣, J♡, J♢, 6♠, 6♣)) == 210    # Jacks Full over Sixes
@test first(evaluate5(J♣, J♡, J♢, 5♠, 5♣)) == 211    # Jacks Full over Fives
@test first(evaluate5(J♣, J♡, J♢, 4♠, 4♣)) == 212    # Jacks Full over Fours
@test first(evaluate5(J♣, J♡, J♢, 3♠, 3♣)) == 213    # Jacks Full over Treys
@test first(evaluate5(J♣, J♡, J♢, 2♠, 2♣)) == 214    # Jacks Full over Deuces
@test first(evaluate5(T♣, T♡, T♢, A♠, A♣)) == 215    # Tens Full over Aces
@test first(evaluate5(T♣, T♡, T♢, K♠, K♣)) == 216    # Tens Full over Kings
@test first(evaluate5(T♣, T♡, T♢, Q♠, Q♣)) == 217    # Tens Full over Queens
@test first(evaluate5(T♣, T♡, T♢, J♠, J♣)) == 218    # Tens Full over Jacks
@test first(evaluate5(T♣, T♡, T♢, 9♠, 9♣)) == 219    # Tens Full over Nines
@test first(evaluate5(T♣, T♡, T♢, 8♠, 8♣)) == 220    # Tens Full over Eights
@test first(evaluate5(T♣, T♡, T♢, 7♠, 7♣)) == 221    # Tens Full over Sevens
@test first(evaluate5(T♣, T♡, T♢, 6♠, 6♣)) == 222    # Tens Full over Sixes
@test first(evaluate5(T♣, T♡, T♢, 5♠, 5♣)) == 223    # Tens Full over Fives
@test first(evaluate5(T♣, T♡, T♢, 4♠, 4♣)) == 224    # Tens Full over Fours
@test first(evaluate5(T♣, T♡, T♢, 3♠, 3♣)) == 225    # Tens Full over Treys
@test first(evaluate5(T♣, T♡, T♢, 2♠, 2♣)) == 226    # Tens Full over Deuces
@test first(evaluate5(9♣, 9♡, 9♢, A♠, A♣)) == 227    # Nines Full over Aces
@test first(evaluate5(9♣, 9♡, 9♢, K♠, K♣)) == 228    # Nines Full over Kings
@test first(evaluate5(9♣, 9♡, 9♢, Q♠, Q♣)) == 229    # Nines Full over Queens
@test first(evaluate5(9♣, 9♡, 9♢, J♠, J♣)) == 230    # Nines Full over Jacks
@test first(evaluate5(9♣, 9♡, 9♢, T♠, T♣)) == 231    # Nines Full over Tens
@test first(evaluate5(9♣, 9♡, 9♢, 8♠, 8♣)) == 232    # Nines Full over Eights
@test first(evaluate5(9♣, 9♡, 9♢, 7♠, 7♣)) == 233    # Nines Full over Sevens
@test first(evaluate5(9♣, 9♡, 9♢, 6♠, 6♣)) == 234    # Nines Full over Sixes
@test first(evaluate5(9♣, 9♡, 9♢, 5♠, 5♣)) == 235    # Nines Full over Fives
@test first(evaluate5(9♣, 9♡, 9♢, 4♠, 4♣)) == 236    # Nines Full over Fours
@test first(evaluate5(9♣, 9♡, 9♢, 3♠, 3♣)) == 237    # Nines Full over Treys
@test first(evaluate5(9♣, 9♡, 9♢, 2♠, 2♣)) == 238    # Nines Full over Deuces
@test first(evaluate5(8♣, 8♡, 8♢, A♠, A♣)) == 239    # Eights Full over Aces
@test first(evaluate5(8♣, 8♡, 8♢, K♠, K♣)) == 240    # Eights Full over Kings
@test first(evaluate5(8♣, 8♡, 8♢, Q♠, Q♣)) == 241    # Eights Full over Queens
@test first(evaluate5(8♣, 8♡, 8♢, J♠, J♣)) == 242    # Eights Full over Jacks
@test first(evaluate5(8♣, 8♡, 8♢, T♠, T♣)) == 243    # Eights Full over Tens
@test first(evaluate5(8♣, 8♡, 8♢, 9♠, 9♣)) == 244    # Eights Full over Nines
@test first(evaluate5(8♣, 8♡, 8♢, 7♠, 7♣)) == 245    # Eights Full over Sevens
@test first(evaluate5(8♣, 8♡, 8♢, 6♠, 6♣)) == 246    # Eights Full over Sixes
@test first(evaluate5(8♣, 8♡, 8♢, 5♠, 5♣)) == 247    # Eights Full over Fives
@test first(evaluate5(8♣, 8♡, 8♢, 4♠, 4♣)) == 248    # Eights Full over Fours
@test first(evaluate5(8♣, 8♡, 8♢, 3♠, 3♣)) == 249    # Eights Full over Treys
@test first(evaluate5(8♣, 8♡, 8♢, 2♠, 2♣)) == 250    # Eights Full over Deuces
@test first(evaluate5(7♣, 7♡, 7♢, A♠, A♣)) == 251    # Sevens Full over Aces
@test first(evaluate5(7♣, 7♡, 7♢, K♠, K♣)) == 252    # Sevens Full over Kings
@test first(evaluate5(7♣, 7♡, 7♢, Q♠, Q♣)) == 253    # Sevens Full over Queens
@test first(evaluate5(7♣, 7♡, 7♢, J♠, J♣)) == 254    # Sevens Full over Jacks
@test first(evaluate5(7♣, 7♡, 7♢, T♠, T♣)) == 255    # Sevens Full over Tens
@test first(evaluate5(7♣, 7♡, 7♢, 9♠, 9♣)) == 256    # Sevens Full over Nines
@test first(evaluate5(7♣, 7♡, 7♢, 8♠, 8♣)) == 257    # Sevens Full over Eights
@test first(evaluate5(7♣, 7♡, 7♢, 6♠, 6♣)) == 258    # Sevens Full over Sixes
@test first(evaluate5(7♣, 7♡, 7♢, 5♠, 5♣)) == 259    # Sevens Full over Fives
@test first(evaluate5(7♣, 7♡, 7♢, 4♠, 4♣)) == 260    # Sevens Full over Fours
@test first(evaluate5(7♣, 7♡, 7♢, 3♠, 3♣)) == 261    # Sevens Full over Treys
@test first(evaluate5(7♣, 7♡, 7♢, 2♠, 2♣)) == 262    # Sevens Full over Deuces
@test first(evaluate5(6♣, 6♡, 6♢, A♠, A♣)) == 263    # Sixes Full over Aces
@test first(evaluate5(6♣, 6♡, 6♢, K♠, K♣)) == 264    # Sixes Full over Kings
@test first(evaluate5(6♣, 6♡, 6♢, Q♠, Q♣)) == 265    # Sixes Full over Queens
@test first(evaluate5(6♣, 6♡, 6♢, J♠, J♣)) == 266    # Sixes Full over Jacks
@test first(evaluate5(6♣, 6♡, 6♢, T♠, T♣)) == 267    # Sixes Full over Tens
@test first(evaluate5(6♣, 6♡, 6♢, 9♠, 9♣)) == 268    # Sixes Full over Nines
@test first(evaluate5(6♣, 6♡, 6♢, 8♠, 8♣)) == 269    # Sixes Full over Eights
@test first(evaluate5(6♣, 6♡, 6♢, 7♠, 7♣)) == 270    # Sixes Full over Sevens
@test first(evaluate5(6♣, 6♡, 6♢, 5♠, 5♣)) == 271    # Sixes Full over Fives
@test first(evaluate5(6♣, 6♡, 6♢, 4♠, 4♣)) == 272    # Sixes Full over Fours
@test first(evaluate5(6♣, 6♡, 6♢, 3♠, 3♣)) == 273    # Sixes Full over Treys
@test first(evaluate5(6♣, 6♡, 6♢, 2♠, 2♣)) == 274    # Sixes Full over Deuces
@test first(evaluate5(5♣, 5♡, 5♢, A♠, A♣)) == 275    # Fives Full over Aces
@test first(evaluate5(5♣, 5♡, 5♢, K♠, K♣)) == 276    # Fives Full over Kings
@test first(evaluate5(5♣, 5♡, 5♢, Q♠, Q♣)) == 277    # Fives Full over Queens
@test first(evaluate5(5♣, 5♡, 5♢, J♠, J♣)) == 278    # Fives Full over Jacks
@test first(evaluate5(5♣, 5♡, 5♢, T♠, T♣)) == 279    # Fives Full over Tens
@test first(evaluate5(5♣, 5♡, 5♢, 9♠, 9♣)) == 280    # Fives Full over Nines
@test first(evaluate5(5♣, 5♡, 5♢, 8♠, 8♣)) == 281    # Fives Full over Eights
@test first(evaluate5(5♣, 5♡, 5♢, 7♠, 7♣)) == 282    # Fives Full over Sevens
@test first(evaluate5(5♣, 5♡, 5♢, 6♠, 6♣)) == 283    # Fives Full over Sixes
@test first(evaluate5(5♣, 5♡, 5♢, 4♠, 4♣)) == 284    # Fives Full over Fours
@test first(evaluate5(5♣, 5♡, 5♢, 3♠, 3♣)) == 285    # Fives Full over Treys
@test first(evaluate5(5♣, 5♡, 5♢, 2♠, 2♣)) == 286    # Fives Full over Deuces
@test first(evaluate5(4♣, 4♡, 4♢, A♠, A♣)) == 287    # Fours Full over Aces
@test first(evaluate5(4♣, 4♡, 4♢, K♠, K♣)) == 288    # Fours Full over Kings
@test first(evaluate5(4♣, 4♡, 4♢, Q♠, Q♣)) == 289    # Fours Full over Queens
@test first(evaluate5(4♣, 4♡, 4♢, J♠, J♣)) == 290    # Fours Full over Jacks
@test first(evaluate5(4♣, 4♡, 4♢, T♠, T♣)) == 291    # Fours Full over Tens
@test first(evaluate5(4♣, 4♡, 4♢, 9♠, 9♣)) == 292    # Fours Full over Nines
@test first(evaluate5(4♣, 4♡, 4♢, 8♠, 8♣)) == 293    # Fours Full over Eights
@test first(evaluate5(4♣, 4♡, 4♢, 7♠, 7♣)) == 294    # Fours Full over Sevens
@test first(evaluate5(4♣, 4♡, 4♢, 6♠, 6♣)) == 295    # Fours Full over Sixes
@test first(evaluate5(4♣, 4♡, 4♢, 5♠, 5♣)) == 296    # Fours Full over Fives
@test first(evaluate5(4♣, 4♡, 4♢, 3♠, 3♣)) == 297    # Fours Full over Treys
@test first(evaluate5(4♣, 4♡, 4♢, 2♠, 2♣)) == 298    # Fours Full over Deuces
@test first(evaluate5(3♣, 3♡, 3♢, A♠, A♣)) == 299    # Treys Full over Aces
@test first(evaluate5(3♣, 3♡, 3♢, K♠, K♣)) == 300    # Treys Full over Kings
@test first(evaluate5(3♣, 3♡, 3♢, Q♠, Q♣)) == 301    # Treys Full over Queens
@test first(evaluate5(3♣, 3♡, 3♢, J♠, J♣)) == 302    # Treys Full over Jacks
@test first(evaluate5(3♣, 3♡, 3♢, T♠, T♣)) == 303    # Treys Full over Tens
@test first(evaluate5(3♣, 3♡, 3♢, 9♠, 9♣)) == 304    # Treys Full over Nines
@test first(evaluate5(3♣, 3♡, 3♢, 8♠, 8♣)) == 305    # Treys Full over Eights
@test first(evaluate5(3♣, 3♡, 3♢, 7♠, 7♣)) == 306    # Treys Full over Sevens
@test first(evaluate5(3♣, 3♡, 3♢, 6♠, 6♣)) == 307    # Treys Full over Sixes
@test first(evaluate5(3♣, 3♡, 3♢, 5♠, 5♣)) == 308    # Treys Full over Fives
@test first(evaluate5(3♣, 3♡, 3♢, 4♠, 4♣)) == 309    # Treys Full over Fours
@test first(evaluate5(3♣, 3♡, 3♢, 2♠, 2♣)) == 310    # Treys Full over Deuces
@test first(evaluate5(2♣, 2♡, 2♢, A♠, A♣)) == 311    # Deuces Full over Aces
@test first(evaluate5(2♣, 2♡, 2♢, K♠, K♣)) == 312    # Deuces Full over Kings
@test first(evaluate5(2♣, 2♡, 2♢, Q♠, Q♣)) == 313    # Deuces Full over Queens
@test first(evaluate5(2♣, 2♡, 2♢, J♠, J♣)) == 314    # Deuces Full over Jacks
@test first(evaluate5(2♣, 2♡, 2♢, T♠, T♣)) == 315    # Deuces Full over Tens
@test first(evaluate5(2♣, 2♡, 2♢, 9♠, 9♣)) == 316    # Deuces Full over Nines
@test first(evaluate5(2♣, 2♡, 2♢, 8♠, 8♣)) == 317    # Deuces Full over Eights
@test first(evaluate5(2♣, 2♡, 2♢, 7♠, 7♣)) == 318    # Deuces Full over Sevens
@test first(evaluate5(2♣, 2♡, 2♢, 6♠, 6♣)) == 319    # Deuces Full over Sixes
@test first(evaluate5(2♣, 2♡, 2♢, 5♠, 5♣)) == 320    # Deuces Full over Fives
@test first(evaluate5(2♣, 2♡, 2♢, 4♠, 4♣)) == 321    # Deuces Full over Fours
@test first(evaluate5(2♣, 2♡, 2♢, 3♠, 3♣)) == 322    # Deuces Full over Treys


text:
module LabelNumerals

import Base: ==, isless, <, >, <=, +, -, max, min,
             convert, promote_rule, print, show, string, hash, parse

include("LabelNumeral.jl")
include("AlphaNumeral.jl")
include("ExtNumerals.jl")
include("LookupNumeral.jl")
include("AlphaNumNumeral.jl")

export LabelNumeral, findLabels,
       @an_str, AlphaNumeral,
       @ln_str, LookupNumeral, registerLookupNumerals,
       @ann_str, AlphaNumNumeral

allNumerals = [AlphaNumeral, RomanNumeral, Int, LookupNumeral, AlphaNumNumeral]

"""
```
    findLabels(label::String, ::Vector{DataType} = allNumerals; pfxList=Vector{String}=[""])
        -> Vector{Tuple{LabelNumeral,Type}}}
```
Given `allNumerals =[AlphaNumeral, RomanNumeral, Int, LookupNumeral, AlphaNumNumeral]`

Finds the `LabelNumeral` that is most suitably matching to the input `String`.
`pfxList` provides one or more label prefix values.

The function returns an array of all the matching `LabelNumeral` and the `Type` of numeral
that best matches its internal composition.
"""
function findLabels(str::String, types::Vector{DataType} = allNumerals; pfxList=Vector{String}=[""])
    arr::Vector{Tuple{LabelNumeral,Type}} = Vector{Tuple{LabelNumeral,Type}}()
    for pfx in pfxList
        if(startswith(str, pfx))
            tstr = str[length(pfx)+1:end]
            for T in types
                try
                    ln = LabelNumeral(T, tstr; prefix=pfx)
                    push!(arr, (ln, T))
                catch
                end
            end
        end
    end
    return arr
end

end # module


text:
using Materialis
using Test

@testset "Test: FEM2Grid interpolation matrix" begin

    # generate the FEM mesh
    Lx = 0.5
    Ly = 1.0
    Lz = 1.2

    testNodes = [ 0    0    0 ;
                  0    0   Lz ;
                  0   Ly   Lz ;
                  0   Ly    0 ;
                  Lx   0    0 ;
                  Lx   0   Lz ;
                  Lx  Ly   Lz ;
                  Lx  Ly    0 ]

    testConec = [ [ 1, 4, 2, 6 ],
                  [ 6, 2, 3, 4 ],
                  [ 4, 3, 6, 7 ],
                  [ 4, 1, 5, 6 ],
                  [ 4, 6, 5, 8 ],
                  [ 4, 7, 6, 8 ] ]

    # generate the image grid
    iniG = -0.2
    endG = 1.0
    numVox = [ 10, 5, 8]
    testGrid = createGrid( iniG*ones(3), endG*ones(3), numVox )

    # generate synthetic intensity in FEM mesh
    Ix, Iy, Iz = [.2, .3, -.4]
    magnitudeFEM = [ testNodes[i,1]*Ix+testNodes[i,2]*Iy+testNodes[i,3]*Iz  for i in (1:8) ]

    # interpolate intensity
    interpMatrix = computeFEM2GridInterpMatrix( testNodes, testConec, testGrid, 1, 1 )
    magnitudeGrid = interpMatrix * magnitudeFEM

    # paraview output of unstructured FEM mesh data
    vtkUnsGridPlot( testConec, testNodes, magnitudeFEM, "testunstrgrid" )

    # paraview output of structured grid data
    magnitudes = reshape( magnitudeGrid, ( numVox[1], numVox[2], numVox[3] ) )
    vtkStrGridPlot( testGrid, [magnitudes], "teststrgrid" )

    numericalInte = magnitudes[3,2,2]
    analyticaInte = sum( ( testGrid.voxelWidths .* [2,1,1] .+ testGrid.startVox ) .* [Ix,Iy,Iz] )

    @test abs( numericalInte - analyticaInte ) < ( 1.0e-8*abs(analyticaInte) )
end


text:
export MatrixOp

"""
`MatrixOp(domainType=Float64::Type, dim_in::Tuple, A::AbstractMatrix)`

`MatrixOp(A::AbstractMatrix)`

`MatrixOp(A::AbstractMatrix, n_colons)`

Creates a `LinearOperator` which, when multiplied with a vector `x::AbstractVector`, returns the product `A*x`.

The input `x` can be also a matrix: the number of columns must be given either in the second entry of `dim_in::Tuple` or using the constructor `MatrixOp(A::AbstractMatrix, n_colons)`.

```julia
julia> MatrixOp(Float64,(10,),randn(20,10))
▒  ℝ^10 -> ℝ^20

julia> MatrixOp(randn(20,10))
▒  ℝ^10 -> ℝ^20

julia> MatrixOp(Float64,(10,20),randn(20,10))
▒  ℝ^(10, 20) -> ℝ^(20, 20)

julia> MatrixOp(randn(20,10),4)
▒  ℝ^(10, 4) -> ℝ^(20, 4)

```

"""
struct MatrixOp{D, T, M <: AbstractMatrix{T}} <: LinearOperator
	A::M
	n_col_in::Integer
end

# Constructors

##TODO decide what to do when domainType is given, with conversion one loses pointer to data...
###standard constructor Operator{N}(DomainType::Type, DomainDim::NTuple{N,Int})
function MatrixOp(DomainType::Type, DomainDim::NTuple{N,Int}, A::M) where {N, T, M <: AbstractMatrix{T}}
	N > 2 && error("cannot multiply a Matrix by a n-dimensional Variable with n > 2")
	size(A,2) != DomainDim[1] && error("wrong input dimensions")
	if N == 1
		MatrixOp{DomainType, T, M}(A, 1)
	else
		MatrixOp{DomainType, T, M}(A, DomainDim[2])
	end
end
###

MatrixOp(A::M)                      where {M <: AbstractMatrix} = MatrixOp(eltype(A), (size(A,2),), A)
MatrixOp(D::Type, A::M)             where {M <: AbstractMatrix} = MatrixOp(D, (size(A,2),), A)
MatrixOp(A::M, n::Integer)          where {M <: AbstractMatrix} = MatrixOp(eltype(A), (size(A,2), n), A)
MatrixOp(D::Type, A::M, n::Integer) where {M <: AbstractMatrix} = MatrixOp(D, (size(A,2), n), A)

import Base: convert
convert(::Type{LinearOperator}, L::M) where {T, M<:AbstractMatrix{T}} = MatrixOp{T, T, M}(L,1)
convert(::Type{LinearOperator}, L::M, n::Integer) where {T, M<:AbstractMatrix{T}} = MatrixOp{T, T, M}(L, n)
convert(::Type{LinearOperator}, dom::Type, dim_in::Tuple, L::AbstractMatrix) = MatrixOp(dom, dim_in, L)

# Mappings

mul!(y::AbstractArray, L::MatrixOp{D, T, M}, b::AbstractArray) where {D, T, M} = mul!(y, L.A, b)
mul!(y::AbstractArray, L::AdjointOperator{MatrixOp{D, T, M}}, b::AbstractArray) where {D, T, M} = mul!(y, L.A.A', b)

# Special Case, real b, complex matrix
function mul!(y::AbstractArray, L::AdjointOperator{MatrixOp{D, T, M}}, b::AbstractArray) where {D <: Real , T <: Complex, M}
	yc = zeros(T,size(y))
	mul!(yc, L.A.A', b)
	y .= real.(yc)
end

# Properties

  domainType(L::MatrixOp{D, T}) where {D, T} = D
codomainType(L::MatrixOp{D, T}) where {D, T} = D <: Real && T <: Complex ? T : D

function size(L::MatrixOp)
	if L.n_col_in == 1
		( (size(L.A, 1),), (size(L.A, 2),) )
	else
		( (size(L.A, 1), L.n_col_in), (size(L.A, 2), L.n_col_in) )
	end
end

fun_name(L::MatrixOp) = "▒"

is_diagonal(L::MatrixOp) = isdiag(L.A)
is_full_row_rank(L::MatrixOp) = rank(L.A) == size(L.A, 1)
is_full_column_rank(L::MatrixOp) = rank(L.A) == size(L.A, 2)


text:
using QCMaterial
using QCMaterial.HartreeFock
using Test
using PyCall
using Requires
import PyCall: pyimport
import QCMaterial: uccgsd, convert_openfermion_op, up_index, down_index

function __init__()
    # Enable MPI tests only when MPI.jl is loaded
    @require MPI="da04e1cc-30fd-572f-bb4f-1f8673147195" include("mpi.jl")
end

include("util.jl")
include("core.jl")
include("hamiltonian.jl")
include("hartree_fock.jl")
include("ucc.jl")
include("computation.jl")


text:
using Dates

"""
    {{meta}}
Plug in specific meta information for a blog page. The `meta` local page
variable should be given as a list of tuples of pairs like so:
```
@def meta = [("property"=>"og:video", "content"=>"http://example.com/"),
             ("propery"=>"og:title", "content"=>"The Rock")]
```
"""
function hfun_meta()
    m = locvar(:meta)
    io = IOBuffer()
    for tuple in locvar(:meta)
        write(io, "<meta ")
        for (prop, val) in tuple
            write(io, "$prop=\"$val\" ")
        end
        write(io, ">")
    end
    return String(take!(io))
end

function two_digits(n::Int)::String
    pre = ifelse(div(n, 10) == 0, "0", "")
    return "$pre$n"
end

"""
    {{blogposts}}
Plug in the list of blog posts contained in the `/blog/` folder.
"""
@delay function hfun_blogposts()
    curyear = Dates.Year(Dates.today()).value
    io = IOBuffer()
    for year in curyear:-1:2016
        ys = "$year"
        year < curyear && write(io, "\n**$year**\n")
        for month in 12:-1:1
            ms = two_digits(month)
            # if you have an additional level of structure...
            for day in 31:-1:1
                ds = two_digits(day)
                # check the folder
                base = joinpath("news", ys, ms, ds)
                isdir(base) || continue
                # check the posts in the folder
                posts = filter!(p -> endswith(p, ".md"), readdir(base))
                for post in posts
                    ps       = splitext(post)[1]
                    url      = "/news/$ys/$ms/$ds/$ps/"
                    surl     = strip(url, '/')
                    title    = pagevar(surl, :title)
                    author   = pagevar(surl, :author)
                    isnothing(author) && (author = "")
                    pubdate  = "$ys-$ms-$ds"
                    write(io, "\n[$title]($url) $author ($pubdate) \n")
                end
            end
        end
    end
    # markdown conversion adds `<p>` beginning and end but
    # we want to  avoid this to avoid an empty separator
    r = Franklin.fd2html(String(take!(io)), internal=true)
    return r
end

function redirect(url)
    s = """
    <!-- REDIRECT - DO NOT MODIFY-->
    <!doctype html>
    <html>
      <head>
        <meta http-equiv="refresh" content="0; url=$url">
      </head>
    </html>
    """
    return s
end

function hfun_blog_post_redirects()
    basepath = Franklin.FOLDER_PATH[]
    for (root, _, files) in walkdir(joinpath(basepath, "news"))
        for file in files
            fullpath = joinpath(root, file)
            fullpath == joinpath(basepath, "news", "index.md") && continue
            # relative path: /news/...
            relpath = splitext(fullpath[length(basepath)+1:end])[1]
            dst = joinpath(basepath, "__site", strip(relpath, '/') * ".html")
            isfile(dst) && continue
            mkpath(splitdir(dst)[1])
            write(dst, redirect(relpath * "/"))
        end
    end
    return ""
end

function hfun_add_redirects()
  basepath = Franklin.FOLDER_PATH[]
  for n in readdir(basepath)
    endswith(n, ".md") || continue
    n in ("index.md", "config.md") && continue
    name = splitext(n)[1]
    dst = joinpath(basepath, "__site", name * ".html")
    isfile(dst) && continue
    prepath = globvar("prepath")
    pre = ""
    if prepath !== nothing && !isempty(prepath)
      pre = "/" * strip(prepath, '/')
    end
    s = redirect("$pre/$name/")
    write(dst, s)
  end
  return ""
end


text:
function deltasum(p0, ps, k, m2, T,δk=0.1)
    if p0 + δk > 2 * Epi(k, m2) > p0 - δk
        return (k * (-2 * k + ps) * coth(Epi(k, m2) / (2 * T))) /
               (8 * Epi(k, m2)^2 * (2 * Epi(k, m2) + p0) * pi^2)
    else
        return (k * (-2 * k + ps) * coth(sqrt(k^2 + m2) / (2 * T))) /
               (2 * sqrt(k^2 + m2) * (4 * (k^2 + m2) - p0^2) * pi^2)
    end
end

function F1tilde_delta(p0, p, k, m, T)
    (
        k *
        (2 * k - p) *
        csch(sqrt(k^2 + m) / (2 * T))^2 *
        (
            sqrt(k^2 + m) *
            (2 * k - p) *
            (4 * k + p) *
            (-4 * (k^2 + m) + p0^2) +
            (
                (2 * k - p) * (4 * k + p) * p0^2 +
                12 * (k^2 + m) * (8 * m + 2 * k * p + p^2 - 2 * p0^2)
            ) *
            T *
            sinh(sqrt(k^2 + m) / T)
        )
    ) / (192 * (k^2 + m)^(3 / 2) * (-4 * (k^2 + m) + p0^2)^2 * pi^2 * T)
end

function dkF1Tilde(p0, ps, k, m, T;kwarg...)
    if k > ps / 2
        hquadrature(
            x -> 2 * x * (-p0^2 + x^2)^-1 * ppfun(x, ps, k, Epi(k, m), T),
            2 * Epi(k, m),
            Epi(k + ps, m) + Epi(k, m);kwarg...
        )[1] +
        deltasum(p0, ps, k, m, T) +
        F1tilde_delta(p0, ps, k, m, T)
    elseif k <= ps / 2
        hquadrature(
            x -> 2 * x * (-p0^2 + x^2)^-1 * ppfun(x, ps, k, Epi(k, m), T),
            2 * Epi(k, m),
            Epi(k + ps, m) + Epi(k, m);kwarg...
        )[1]
    end
end

function dkF2Tilde(p0, ps, k, m, T;kwarg...)
    if k > ps / 2
        if Epi(k + ps, m) - Epi(k, m) < p0
            return hquadrature(
                x -> 2 * x * (-p0^2 + x^2)^-1 * pmfun(x, ps, k, Epi(k, m), T),
                0.0,
                Epi(k + ps, m) - Epi(k, m);kwarg...
            )[1]+pmfun_zerofix(p0, ps, k, m, T,9)
        else
            return quadgk_PV(
                x -> pmfun(x, ps, k, Epi(k, m), T),
                0.0,
                p0,
                Epi(ps + k, m) - Epi(k, m);kwarg...
            )[1]+pmfun_zerofix(p0, ps, k, m, T,9)
        end
    elseif k <= ps / 2
        if Epi(k + ps, m) - Epi(k, m) < p0 || p0 < Epi(k - ps, m) - Epi(k, m)
            return hquadrature(
                x -> 2 * x * (-p0^2 + x^2)^-1 * pmfun(x, ps, k, Epi(k, m), T),
                Epi(k - ps, m) - Epi(k, m),
                Epi(k + ps, m) - Epi(k, m);kwarg...
            )[1]
        else
            return quadgk_PV(
                x -> pmfun(x, ps, k, Epi(k, m), T),
                Epi(k - ps, m) - Epi(k, m),
                p0,
                Epi(ps + k, m) - Epi(k, m);kwarg...
            )[1]
        end
    end
end

dkF1TildeAll(q0, qs, k, m, T;kwarg...)=dkF1Tilde(abs(q0), qs, k, m, T;kwarg...)
dkF2TildeAll(q0, qs, k, m, T;kwarg...)=dkF2Tilde(abs(q0), qs, k, m, T;kwarg...)

function dklamdabar_Zero(q0, qs, k, m, T, lam4pik, Npi;kwarg...)
        lam4pik^2 *
        (
            (2 + Npi) * (dkF1TildeAll(k, m, T)+dkF2TildeAll(k, m, T)) +
             6 * dkF1TildeAll(q0, qs, k, m, T;kwarg...) +
             6 * dkF2TildeAll(q0, qs, k, m, T;kwarg...)
        )
end

function lamdabar_Zero(q0,qs,k, T, Npi, mfun,lamfun,UVScale;kwarg...)
    -hquadrature(
        kprim -> dklamdabar_Zero(q0, qs, kprim, mfun(kprim), T, lamfun(kprim), Npi;kwarg...),
        k,
        UVScale;kwarg...
        # initdiv=200,
    )[1] + lamfun(UVScale)
end


text:
using NBody
using JSON

nb = read_nbody()
# println(json(nb))
NBody.write_snapshot(nb)
NBody.write_snapshot(nb)


text:
function default_integer_type()
    return typeof(1)
end

function default_float_type()
    return typeof(1.0)
end

function reference_element_length()
    a,b = reference_interval_1d()
    return b-a
end

function in_reference_interval(x)
    a,b = reference_interval_1d()
    return a <= x <= b
end

function reference_cell(dim)
    if dim == 1
        xL = [-1.0]
        xR = [+1.0]
        return IntervalBox(xL,xR)
    elseif dim == 2
        xL = [-1.0,-1.0]
        xR = [+1.0,+1.0]
        return IntervalBox(xL,xR)
    else
        throw(ArgumentError("Expected dim ∈ {1,2}, got dim = $dim"))
    end
end

function number_of_faces(dim)
    if dim == 2
        return 4
    elseif dim == 3
        return 6
    else
        throw(ArgumentError("Expected dim ∈ {2,3}, got $dim == dim"))
    end
end

function reference_normals()
    return [[0.0,-1.0],[1.0,0.0],[0.0,1.0],[-1.0,0.0]]
end

function reference_cell_volume(dim)
    return reference_element_length()^dim
end

function neighbor_faceid(faceid)
    if faceid == 1
        return 3
    elseif faceid == 2
        return 4
    elseif faceid == 3
        return 1
    elseif faceid == 4
        return 2
    else
        throw(ArgumentError("Expected faceid ∈ {1,2,3,4}, got faceid = $faceid"))
    end
end

function jacobian(map::InterpolatingPolynomial,p)
    return gradient(map,p)
end

function determinant_jacobian(map::InterpolatingPolynomial{2},p)
    return norm(jacobian(map,p))
end

function nodal_coordinates(mesh,basis)

    dim = dimension(mesh)
    @assert dimension(basis) == dim
    NF = number_of_basis_functions(basis)

    ncells = mesh.total_number_of_elements
    coords = zeros(dim,NF*ncells)
    start = 1
    stop = NF
    for idx in 1:ncells
        xL,xR = CartesianMesh.element(mesh,idx)
        map = CellMap(xL,xR)
        coords[:,start:stop] = map(basis.points)
        start = stop+1
        stop += NF
    end
    return coords
end

function restrict_on_faces(func,box::IntervalBox{2})

    fb(x) = func(extend(x,2,box[2].lo))
    fr(y) = func(extend(y,1,box[1].hi))
    ft(x) = func(extend(x,2,box[2].hi))
    fl(y) = func(extend(y,1,box[1].lo))

    return [fb,fr,ft,fl]
end

function dimension(basis::TensorProductBasis{dim}) where {dim}
    return dim
end

function dimension(poly::InterpolatingPolynomial)
    return dimension(poly.basis)
end

function dimension(quad::QuadratureRule{dim}) where dim
    return dim
end

function dimension(mesh::UniformMesh{dim}) where dim
    return dim
end

function number_of_basis_functions(basis::T) where
    {T<:PolynomialBasis.AbstractBasis{dim,NF}} where {dim,NF}

    return NF
end

function number_of_quadrature_points(quad::QuadratureRule{D,Q}) where {D,Q}
    return Q
end

function Base.length(quad::QuadratureRule)
    return number_of_quadrature_points(quad)
end

function levelset_normal(poly,p::V,cellmap) where {V<:AbstractVector}
    g = vec(gradient(poly,p))
    invjac = inverse_jacobian(cellmap)
    n = diagm(invjac)*g
    return n/norm(n)
end

function levelset_normal(poly,p::M,cellmap) where {M<:AbstractMatrix}
    g = hcat([gradient(poly,p[:,i])' for i in 1:size(p)[2]]...)
    invjac = inverse_jacobian(cellmap)
    normals = diagm(invjac)*g
    for i in 1:size(normals)[2]
        n = normals[:,i]
        normals[:,i] = n/norm(n)
    end
    return normals
end

function tangents(normal::V) where {V<:AbstractVector}
    @assert length(normal) == 2
    return [normal[2],-normal[1]]
end

function tangents(normals::M) where {M<:AbstractMatrix}
    m,n = size(normals)
    @assert m == 2

    t = similar(normals)
    t[1,:] = normals[2,:]
    t[2,:] = -normals[1,:]
    return t
end

function scale_area(cellmap,normals)
    t = tangents(normals)
    invjac = inverse_jacobian(cellmap)
    den = sqrt.((t.^2)'*(invjac.^2))
    return 1.0 ./ den
end

function active_faces(facequads)
    isactiveface = [length(fq) > 0 ? true : false for fq in facequads]
end

function face_midpoints(cellmap::CellMap{2})
    xl = cellmap.xL
    xr = cellmap.xR
    xm = 0.5*(xl+xr)

    m1 = [xm[1],xl[2]]
    m2 = [xr[1],xm[2]]
    m3 = [xm[1],xr[2]]
    m4 = [xl[1],xm[2]]

    return [m1,m2,m3,m4]
end

function levelset_coefficients(distance_function,mesh,basis)
    coords = nodal_coordinates(mesh,basis)
    NF = number_of_basis_functions(basis)
    return reshape(distance_function(coords),NF,:)
end

function curve_normals(imap,points)
    npts = length(points)
    normals = zeros(2,npts)
    for (idx,x) in enumerate(points)
        t = gradient(imap,x)
        n = [-t[2],t[1]]
        normals[:,idx] .= n/norm(n)
    end
    return normals
end

function assign_uniform_cell_hybrid_coordinates!(coordinates,visited,
    phaseid,cellid,dgmesh,facetohelid,interfacehelid,refcoordinates,icoeffs)

    fhelid = facetohelid[phaseid,cellid]
    cellmap = CellMap(dgmesh.domain[cellid])
    for (faceid,helid) in enumerate(fhelid)
        if helid != 0 && !visited[helid]
            coordinates[helid] = [Array(cellmap(p)) for p in refcoordinates[faceid]]
            visited[helid] = true
        end
    end
    ihelid = interfacehelid[phaseid,cellid]
    if ihelid != 0 && !visited[ihelid]
        icoords = reshape(icoeffs[cellid],2,:)
        coordinates[ihelid] = [cellmap(icoords[:,i]) for i = 1:size(icoords)[2]]
        visited[ihelid] = true
    end
end

function hybrid_element_coordinates(dgmesh,ufs,hybrid_element_numbering)
    nhelid = hybrid_element_numbering.number_of_hybrid_elements

    coordinates = Vector{Any}(undef,nhelid)
    visited = zeros(Bool,nhelid)
    facetohelid = hybrid_element_numbering.facetohelid
    interfacehelid = hybrid_element_numbering.interfacehelid
    refcoordinates = [[fm(p) for p in Array(vec(ufs.sbasis.points))] for fm in dgmesh.facemaps]

    for cellid in 1:dgmesh.ncells
        for phaseid in 1:2
            assign_uniform_cell_hybrid_coordinates!(coordinates,visited,
                phaseid,cellid,dgmesh,facetohelid,interfacehelid,
                refcoordinates,ufs.icoeffs)
        end
    end
    coordinates = hcat([hcat(p...) for p in coordinates]...)
    return coordinates
end


text:
# interface proposal by Tom Breloff (@tbreloff)... comments welcome
# This uses the more basic formatting based on FormatSpec and the pyfmt method
# (formerly called fmt, which I repurposed)

# TODO: swap out FormatSpec for something that is able to use the "format" method, which has more options for units, prefixes, etc
# TODO: support rational numbers, autoscale, etc as in "format"

# --------------------------------------------------------------------------------------------------

# the DefaultSpec object is just something to hold onto the current FormatSpec.
# we keep the typechar around specically for the reset! function, to go back to the starting state

type DefaultSpec
    typechar::Char
    fspec::FormatSpec
    DefaultSpec(c::Char) = new(c, FormatSpec(c))
end

const DEFAULT_FORMATTERS = Dict{DataType, DefaultSpec}()

# adds a new default formatter for this type
defaultSpec!{T}(::Type{T}, c::Char) = (DEFAULT_FORMATTERS[T] = DefaultSpec(c); nothing)

# note: types T and K will now both share K's default
defaultSpec!{T,K}(::Type{T}, ::Type{K}) = (DEFAULT_FORMATTERS[T] = DEFAULT_FORMATTERS[K]; nothing)

# seed it with some basic default formatters
for (t, c) in [(Integer,'d'), (AbstractFloat,'f'), (Char,'c'), (AbstractString,'s')]
    defaultSpec!(t, c)
end

reset!{T}(::Type{T}) = (dspec = defaultSpec(T); dspec.fspec = FormatSpec(dspec.typechar); nothing)

# --------------------------------------------------------------------------------------------------

function addKWArgsFromSymbols(kwargs, syms::Symbol...)
    d = Dict(kwargs)
    for s in syms
        if s == :ljust || s == :left
            d[:align] = '<'
        elseif s == :rjust || s == :right
            d[:align] = '>'
        elseif s == :commas
            d[:tsep] = true
        elseif s == :zpad || s == :zeropad
            d[:zpad] = true
        elseif s == :ipre || s == :prefix
            d[:ipre] = true
        end
    end
    d
end

# --------------------------------------------------------------------------------------------------

# methods to get the current default objects
# note: if you want to set a default for an abstract type (i.e. AbstractFloat) you'll need to extend this method like here:
defaultSpec{T<:Integer}(::Type{T}) = DEFAULT_FORMATTERS[Integer]
defaultSpec{T<:AbstractFloat}(::Type{T}) = DEFAULT_FORMATTERS[AbstractFloat]
defaultSpec{T<:AbstractString}(::Type{T}) = DEFAULT_FORMATTERS[AbstractString]
function defaultSpec{T}(::Type{T})
    get(DEFAULT_FORMATTERS, T) do
        error("Missing default spec for type $T... call default!(T, c): $DEFAULT_FORMATTERS")
    end
end
defaultSpec(x) = defaultSpec(typeof(x))

fmt_default{T}(::Type{T}) = defaultSpec(T).fspec
fmt_default(x) = defaultSpec(x).fspec

# first resets the fmt_default spec to the given arg, then continue by updating with args and kwargs
fmt_default!{T}(::Type{T}, c::Char, args...; kwargs...) = (defaultSpec!(T,c); fmt_default!(T, args...; kwargs...))
fmt_default!{T,K}(::Type{T}, ::Type{K}, args...; kwargs...) = (defaultSpec!(T,K); fmt_default!(T, args...; kwargs...))

# update the fmt_default for a specific type
function fmt_default!{T}(::Type{T}, syms::Symbol...; kwargs...)
    if isempty(syms)

        # if there are no arguments, reset to initial defaults
        if isempty(kwargs)
            reset!(T)
            return
        end

        # otherwise update the spec
        dspec = defaultSpec(T)
        dspec.fspec = FormatSpec(dspec.fspec; kwargs...)

    else
        d = addKWArgsFromSymbols(kwargs, syms...)
        fmt_default!(T; d...)
    end
    nothing
end

# update the fmt_default for all types
function fmt_default!(syms::Symbol...; kwargs...)
    if isempty(syms)
        for k in keys(DEFAULT_FORMATTERS)
            fmt_default!(k; kwargs...)
        end
    else
        d = addKWArgsFromSymbols(kwargs, syms...)
        fmt_default!(; d...)
    end
    nothing
end

# --------------------------------------------------------------------------------------------------

# TODO: get rid of this entire hack by moving commas into pyfmt

function optionalCommas(x::Real, s::AbstractString, fspec::FormatSpec)
    dpos = findfirst(s, '.')
    prevwidth = length(s)

    if dpos == 0
        s = addcommas(s)
    else
        s = string(addcommas(s[1:dpos-1]), '.', s[dpos+1:end])
    end

    # check for excess width from commas
    w = length(s)
    if fspec.width > 0 && w > fspec.width && w > prevwidth
        # we may have made the string too wide with those commas... gotta fix it
        s = strip(s)
        n = fspec.width - length(s)
        if fspec.align == '<' # left alignment
            s = string(s, " "^n)
        else
            s = string(" "^n, s)
        end
    end

    s
end
optionalCommas(x, s::AbstractString, fspec::FormatSpec) = s

# --------------------------------------------------------------------------------------------------

# TODO: do more caching to optimize repeated calls

# creates a new FormatSpec by overriding the defaults and passes it to pyfmt
# note: adding kwargs is only appropriate for one-off formatting.
#       normally it will be much faster to change the fmt_default formatting as needed
function fmt(x; kwargs...)
    fspec = isempty(kwargs) ? fmt_default(x) : FormatSpec(fmt_default(x); kwargs...)
    s = pyfmt(fspec, x)

    # add the commas now... I was confused as to when this is done currently
    if fspec.tsep
        return optionalCommas(x, s, fspec)
    end
    s
end

# some helper method calls, which just convert to kwargs
fmt(x, width::Int, args...; kwargs...) = fmt(x, args...; width=width, kwargs...)

fmt(x, width::Int, prec::Int, args...; kwargs...) =
    fmt(x, args...; width=width, prec=prec, kwargs...)

# integrate some symbol shorthands into the keyword args
# note: as above, this will generate relavent kwargs, so to format in a tight loop, you should probably update the fmt_default
function fmt(x, syms::Symbol...; kwargs...)
    d = addKWArgsFromSymbols(kwargs, syms...)
    fmt(x; d...)
end


text:
using OrdinaryDiffEq, DiffEqDevTools, Test
using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear, prob_ode_2Dlinear

@testset "Initial Dt Tests" begin
  prob = prob_ode_linear
  sol =solve(prob,Rosenbrock32())
  dt₀ = sol.t[2]

  prob = prob_ode_2Dlinear
  sol =solve(prob,ExplicitRK(tableau=constructBogakiShampine3()))
  dt₀ = sol.t[2]

  @test  1e-7 < dt₀ < .1
  @test_throws ErrorException sol = solve(prob,Euler())
  #dt₀ = sol.t[2]

  sol3 =solve(prob,ExplicitRK(tableau=constructDormandPrince8_64bit()))
  dt₀ = sol3.t[2]

  @test 1e-7 < dt₀ < .3
end