text:
function solver_EPWE(solver::Solver, ky::Number, omega::Number, polarisation::Polarisation; Norm=false)
    epc = solver.epc.*mask(size(solver.epc)[1])
    muc = solver.muc.*mask(size(solver.muc)[1])
    global K = hcat(diag_R(solver.Ky),diag_R(solver.Kx))/2/pi
    M = length(K[:,1])
    LHS11 = zeros(ComplexF64, M, M)
    LHS12 = eye(M)
    Kp = K + ones(M) * [0 ky]

    if polarisation == TE
        LHS21 = epc * (omega^2*muc*eye(M) .- inv(epc).*(Kp*transpose(Kp)))
        LHS22 = -epc*(inv(epc).*(K[:,1]*ones(1,M)+ones(M)*K[:,1]'))
        LHS = [LHS11 LHS12; LHS21 LHS22]
        RHS = epc
        label = L"H_z"
    elseif polarisation == TM
        LHS21 = (omega^2*epc - diagm(K[:,1].^2 + (ky .+ K[:,2]).^2))
        LHS22 = diagm(-2*K[:,1])
        LHS = [LHS11 LHS12; LHS21 LHS22]
        RHS = muc
        label = L"E_z"
    end
    kT, VT = eigen(LHS,sortby=x -> -abs(x))
    VT = VT[1:M, :]
    # Normalization
    if Norm == true
    VT = VT * diagm(vec(1 ./sqrt.(sum(abs.(VT).^2,dims=1))))
    maxVT,idVT=findmax(abs.(VT), dims=1)
    VT=VT*diagm(vec(maxVT./VT[idVT]))
    end
    #direction
    dT=1/omega*diag_R(VT'*inv(RHS)*(diagm(K[:,1])*VT+VT*diagm(kT)))

    # find the solution in 1st BZ
    i0= (abs.(real.(dT)) .> abs.(imag.(dT))).*((real.(dT)+imag.(dT)) .< 0) + (abs.(real.(dT)) .< abs.(imag.(dT))).*(imag.(kT).< 0)
    i0 = Bool.(i0)

    i3 = []
    for i4 = 0:1
        i0 = .!i0
        i1= 1:2*M
        i1 = i1[i0]
        i2 = sortperm(abs.(real.(kT[i0])))
        buf = sort(abs.(real.(kT[i0])))
        i1=i1[i2]
        i2 = (abs.([imag.(kT[i1]); Inf] - [Inf; imag.(kT[i1])] ) .> 1e-10) .| (abs.([imag.(kT[i1]); Inf]) .<= 1e-10) .| (abs.(abs.([real.(kT[i1]); Inf]) .- abs.([Inf; real.(kT[i1])])) .> 1e-10) .| (abs.(abs.([real.(kT[i1]); Inf]) .+ abs.([Inf; real.(kT[i1])]) .- 1) .>= 0.1)
        i1=i1[i2[1:(length(i2)-1)]]
        i3=[i3; i1[1:Int(sqrt(M))]]
    end
    kT=kT[i3]
    dT=dT[i3]
    VT=VT[:,i3]
    return kT, dT, VT
end

function abcd(kT::VecOrMat,dT::VecOrMat,VT::VecOrMat,omega::Number,ky::Number;background=1)
    M = size(VT)[1]
    m = Int((sqrt(M)-1)/2)
    Ky=Int(round(ky))
    ky=ky-Ky
    msum=zeros((2*m+1),(2*m+1)^2)
    for i0=1:(2*m+1)
        msum[i0,((i0-1)*(2*m+1)+1):(i0*(2*m+1))] .= 1
    end

    ep = background
    abcd1=[-eye(2*m+1) msum*VT[:,1:(2*m+1)]; diagm(sqrt.(Complex.(omega^2 .- ((-m:m) .+ky).^2))) msum*(inv(ep)*(diagm(K[:,1])*VT[:,1:(2*m+1)]+VT[:,1:(2*m+1)]*diagm(kT[1:(2*m+1)])))]
    input = zeros(ComplexF64, 4*m+2)
    input[Ky+m+1] = 1
    input[Ky+3*m+2] = sqrt(Complex(omega^2-(ky+Ky)^2))
    output=abcd1\input
    rT=output[1:(2*m+1)]
    tT=output[(2*m+2):(4*m+2)]

    kxi = sqrt(Complex(omega^2 - ky^2))
    if kxi==0
        kxi=eps()
    end
    kxr=sqrt.(Complex.(omega^2 .- (ky - round(ky) .+ (-m:m)).^2))

    i0=abs.(imag.(kxr[1:(2*m+1)])) .<= 1e-15
    #    i0 = Bool.[i0]
    if length(rT[i0])==0
        RT = 0
    else
        RT = rT[i0]'*diagm(kxr[i0])*rT[i0]/kxi;
    end

    i0 = abs.(imag.(dT[1:(2*m+1)])) .<= abs.(real.(dT[1:(2*m+1)]))
    if length(tT[i0])==0
        TT=0
    else
        TT=tT[i0]'*diagm(dT[1:2m+1][i0])*tT[i0]/kxi*omega;
    end

#    RT = abs(RT); TT = abs(TT)
    return RT, TT, rT, tT
end


text:
using DataStructures
using Base.Meta

import Base.ImmutableDict

const EClassId = Int64

abstract type AbstractENode{T} end

mutable struct ENodeTerm{T} <: AbstractENode{T}
    exprhead::Union{Symbol, Nothing}
    operation::Any
    args::Vector{EClassId}
    hash::Ref{UInt} # hash cache
end

function ENodeTerm{T}(exprhead, operation, c_ids) where {T}
    ENodeTerm{T}(exprhead, operation, c_ids, Ref{UInt}(0))
end

function Base.:(==)(a::ENodeTerm, b::ENodeTerm)
    isequal(a.args, b.args) &&
    isequal(a.exprhead, b.exprhead) && isequal(a.operation, b.operation)
end

TermInterface.istree(n::ENodeTerm) = true
TermInterface.exprhead(n::ENodeTerm) = n.exprhead
TermInterface.operation(n::ENodeTerm) = n.operation
TermInterface.arguments(n::ENodeTerm) = n.args
TermInterface.arity(n::ENodeTerm) = length(n.args)

# This optimization comes from SymbolicUtils
# The hash of an enode is cached to avoid recomputing it.
# Shaves off a lot of time in accessing dictionaries with ENodes as keys.
function Base.hash(t::ENodeTerm{T}, salt::UInt) where {T}
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hash(t.args,  hash(t.exprhead, hash(t.operation, hash(T, salt))))
    t.hash[] = h′
    return h′
end

function toexpr(n::ENodeTerm)
    eh = exprhead(n)
    if isnothing(eh)
        return operation(n) # n is a constant enode
    end
    similarterm(Expr, operation(n), map(i -> Symbol(i, "ₑ"), arguments(n)); exprhead=exprhead(n))
end

# ==================================================
# ENode Literal
# ==================================================

mutable struct ENodeLiteral{T} <: AbstractENode{T}
    value::T
    hash::Ref{UInt}
end

TermInterface.istree(n::ENodeLiteral) = false
TermInterface.exprhead(n::ENodeLiteral) = nothing
TermInterface.operation(n::ENodeLiteral) = n.value
TermInterface.arity(n::ENodeLiteral) = 0

ENodeLiteral(a::T) where{T} = ENodeLiteral{T}(a, Ref{UInt}(0))

Base.:(==)(a::ENodeLiteral, b::ENodeLiteral) = isequal(a.value, b.value)

function Base.hash(t::ENodeLiteral{T}, salt::UInt) where {T}
    !iszero(salt) && return hash(hash(t, zero(UInt)), salt)
    h = t.hash[]
    !iszero(h) && return h
    h′ = hash(t.value, hash(T, salt))
    t.hash[] = h′
    return h′
end

termtype(x::AbstractENode{T}) where T = T

toexpr(n::ENodeLiteral) = operation(n)

function Base.show(io::IO, x::ENodeTerm{T}) where {T}
    print(io, "ENode{$T}(", toexpr(x), ")")
end

Base.show(io::IO, x::ENodeLiteral) = print(io, toexpr(x))

ext: jl


text:
using Jets, JetPack, JetPackDSP, Statistics, Test

@testset "convolve, 1D" begin
    m = sin.([0:0.004:511*0.004;] * 2 * pi * 2.0)

    A = JopConvolve(JetSpace(Float64, 512), JetSpace(Float64, 512), [0.5], x0=0.0, dx=0.004)
    d = A * m
    @test (0.5*m) ≈ (d)

    A = JopConvolve(JetSpace(Float64, 512), JetSpace(Float64, 512), [2.0], x0=-10*0.004, dx=0.004)
    d = A * m
    for i=1:502
        @test isapprox(2.0*m[10+i], d[i], atol=1e-12)
    end
    for i=1:9
        @test isapprox(2.0*m[i], d[502+i], atol=1e-12)
    end

    m = rand(domain(A))
    d = rand(range(A))
    lhs, rhs = dot_product_test(A, m, d)
    @test isapprox((lhs - rhs)/(lhs + rhs), 0.0, atol=1e-10)
end

@testset "convolve, 2D" begin
    A = JopConvolve(JetSpace(Float64,128,128), rand(Float64,32,32), x0=(-16.,-16.), dx=(1.,1.))
    lhs, rhs = dot_product_test(A, rand(domain(A)), rand(range(A)))
    @test lhs ≈ rhs
end

@testset "convolve smoothing, 2D-2D" begin
    # even
    P = JopPad(JetSpace(Float64,128,128), -10:128+11, -10:128+11, extend=true)
    S = JopConvolve(range(P), smoother=:rect, n=(1,1))
    R = JopPad(domain(P), -10:128+11, -10:128+11, extend=false)
    @test iseven(size(range(P),1))
    @test iseven(size(range(P),2))

    m = rand(domain(P))
    d = (R' ∘ S ∘ P)*m
    for i = 2:127, j = 2:127
        @test d[i,j] ≈ mean(m[i-1:i+1,j-1:j+1])
    end
    @test sum(d) ≈ sum(m)

    # odd
    P = JopPad(JetSpace(Float64,128,128), -7:128+9, -7:128+9, extend=true)
    S = JopConvolve(range(P), smoother=:rect, n=(1,1))
    R = JopPad(domain(P), -7:128+9, -7:128+9, extend=false)
    @test isodd(size(range(P),1))
    @test isodd(size(range(P),2))

    m = rand(domain(P))
    d = (R' ∘ S ∘ P)*m
    for i = 2:127, j = 2:127
        @test d[i,j] ≈ mean(m[i-1:i+1,j-1:j+1])
    end
    @test sum(d) ≈ sum(m)
end

@testset "smoothing, 2D-1D" begin
    # even
    P = JopPad(JetSpace(Float64,128,128), -10:128+11, -10:128+11, extend=true)
    S = JopConvolve(range(P), smoother=:rect, n=(1,0))
    R = JopPad(domain(P), -10:128+11, -10:128+11, extend=false)
    @test iseven(size(range(P),1))
    @test iseven(size(range(P),2))

    m = rand(domain(P))
    d = (R' ∘ S ∘ P)*m
    for i = 2:127, j = 2:127
        @test d[i,j] ≈ mean(m[i-1:i+1,j])
    end
    @test sum(d) ≈ sum(m)
    @test sum(d,dims=1) ≈ sum(m,dims=1)

    # odd
    P = JopPad(JetSpace(Float64,128,128), -7:128+9, -7:128+9, extend=true)
    S = JopConvolve(range(P), smoother=:rect, n=(1,0))
    R = JopPad(domain(P), -7:128+9, -7:128+9, extend=false)
    @test isodd(size(range(P),1))
    @test isodd(size(range(P),2))

    m = rand(domain(P))
    d = (R' ∘ S ∘ P)*m
    for i = 2:127, j = 2:127
        @test d[i,j] ≈ mean(m[i-1:i+1,j])
    end
    @test sum(d) ≈ sum(m)
    @test sum(d,dims=1) ≈ sum(m,dims=1)

    # even
    P = JopPad(JetSpace(Float64,128,128), -10:128+11, -10:128+11, extend=true)
    S = JopConvolve(range(P), smoother=:rect, n=(0,1))
    R = JopPad(domain(P), -10:128+11, -10:128+11, extend=false)
    @test iseven(size(range(P),1))
    @test iseven(size(range(P),2))

    m = rand(domain(P))
    d = (R' ∘ S ∘ P)*m
    for i = 2:127, j = 2:127
        @test d[i,j] ≈ mean(m[i,j-1:j+1])
    end
    @test sum(d) ≈ sum(m)
    @test sum(d,dims=2) ≈ sum(m,dims=2)

    # odd
    P = JopPad(JetSpace(Float64,128,128), -7:128+9, -7:128+9, extend=true)
    S = JopConvolve(range(P), smoother=:rect, n=(0,1))
    R = JopPad(domain(P), -7:128+9, -7:128+9, extend=false)
    @test isodd(size(range(P),1))
    @test isodd(size(range(P),2))

    m = rand(domain(P))
    d = (R' ∘ S ∘ P)*m
    for i = 2:127, j = 2:127
        @test d[i,j] ≈ mean(m[i,j-1:j+1])
    end
    @test sum(d) ≈ sum(m)
    @test sum(d,dims=2) ≈ sum(m,dims=2)
end


text:
@inline function save_noise!(W::SimpleNoiseProcess)
  if W.t != W.curt
    push!(W.W,copy(W.curW))
    push!(W.t,copy(W.curt))
    if W.Z != nothing
      push!(W.Z,copy(W.curZ))
    end
  end
end

@inline function accept_step!(W::SimpleNoiseProcess,dt,u,p,setup_next=true)

  W.curt += W.dt
  W.iter += 1

  if isinplace(W)
    @.. W.curW += W.dW
    if W.Z != nothing
      @.. W.curZ += W.dZ
    end
  else
    W.curW += W.dW
    if W.Z != nothing
      W.curZ += W.dZ
    end
  end

  if W.save_everystep
    push!(W.W,copy(W.curW))
    push!(W.t,copy(W.curt))
    if W.Z != nothing
      push!(W.Z,copy(W.curZ))
    end
  end

  W.dt = dt #dtpropose
  # Setup next step
  if setup_next
    setup_next_step!(W::SimpleNoiseProcess,u,p)
  end
end

@inline function setup_next_step!(W::SimpleNoiseProcess,u,p)
  calculate_step!(W,W.dt,u,p)
end

@inline function calculate_step!(W::SimpleNoiseProcess,dt,u,p)
  if isinplace(W)
    W.dist(W.dW,W,dt,u,p,W.curt,W.rng)
    if W.Z != nothing
      W.dist(W.dZ,W,dt,u,p,W.curt,W.rng)
    end
  else
    W.dW = W.dist(W,dt,u,p,W.curt,W.rng)
    if W.Z != nothing
      W.dZ = W.dist(W,dt,u,p,W.curt,W.rng)
    end
  end
  W.dt = dt
end

@inline function reject_step!(W::SimpleNoiseProcess,dtnew,u,p)
  error("SimpleNoiseProcess cannot be used with adaptivity rejections")
end

@inline function interpolate!(W::SimpleNoiseProcess,u,p,t)
  if sign(W.dt)*t > sign(W.dt)*W.t[end] # Steps past W
    dt = t - W.t[end]
    if isinplace(W)
      W.dist(W.dW,W,dt,u,p,t,W.rng)
      W.curW .+= W.dW
      if W.Z != nothing
        W.dist(W.dZ,W,dt,u,p,t,W.rng)
        W.curZ .+= W.dZ
      end
    else
      W.dW = W.dist(W,dt,u,p,t,W.rng)
      W.curW += W.dW
      if W.Z != nothing
        W.dZ = W.dist(W,dt,u,p,t,W.rng)
        W.curZ += W.dZ
      end
    end
    out1 = copy(W.curW)
    if W.save_everystep
      push!(W.t,t)
      push!(W.W,out1)
    end
    if W.Z != nothing
      out2= copy(W.curZ)
      if W.save_everystep
        push!(W.Z,out2)
      end
    else
      out2 = nothing
    end
    return out1,out2
  else # Bridge
    i = searchsortedfirst(W.t,t)
    if t == W.t[i]
      if isinplace(W)
        W.curW .= W.W[i]
      else
        W.curW = W.W[i]
      end
      if W.Z != nothing
        if isinplace(W)
          W.curZ .= W.Z[i]
        else
          W.curZ = W.Z[i]
        end
        return copy(W.curW),copy(W.curZ)
      else
        return copy(W.curW),nothing
      end
    else
      W0,Wh = W.W[i-1],W.W[i]
      if W.Z != nothing
        Z0,Zh = W.Z[i-1],W.Z[i]
      end
      h = W.t[i]-W.t[i-1]
      q = (t-W.t[i-1])/h
      if isinplace(W)
        new_curW = similar(W.dW)
        W.bridge(new_curW,W,W0,Wh,q,h,u,p,t,W.rng)
        if iscontinuous(W)
          @. new_curW += (1-q)*W0
        else
          @. new_curW += W0
        end
        if W.Z != nothing
          new_curZ = similar(W.dZ)
          W.bridge(new_curZ,W,Z0,Zh,q,h,u,p,t,W.rng)
          if iscontinuous(W)
            @. new_curZ += (1-q)*Z0
          else
            @. new_curZ += Z0
          end
        else
          new_curZ = nothing
        end
      else
        new_curW = W.bridge(W,W0,Wh,q,h,u,p,t,W.rng)
        if iscontinuous(W)
          # This should actually be based on the function for computing the mean
          # flow of the noise process, but for now we'll just handle Wiener and
          # Poisson
          new_curW += (1-q)*W0
        else
          new_curW += W0
        end
        if W.Z != nothing
          new_curZ = W.bridge(W,Z0,Zh,q,h,u,p,t,W.rng)
          if iscontinuous(W)
            new_curZ += (1-q)*Z0
          else
            new_curZ += Z0
          end
        else
          new_curZ = nothing
        end
      end
      W.curW = new_curW
      if W.save_everystep
        insert!(W.W,i,new_curW)
        insert!(W.t,i,t)
      end
      if W.Z != nothing
        W.curZ = new_curZ
        if W.save_everystep
          insert!(W.Z,i,new_curZ)
        end
      end
      return new_curW,new_curZ
    end
  end
end

@inline function interpolate!(out1,out2,W::SimpleNoiseProcess,u,p,t)
  if sign(W.dt)*t > sign(W.dt)*W.t[end] # Steps past W
    dt = t - W.t[end]
    W.dist(W.dW,W,dt,u,p,t,W.rng)
    out1 .+= W.dW
    if W.Z != nothing
      W.dist(W.dZ,W,dt,u,p,t,W.rng)
      out2 .+= W.dZ
    end
    if W.save_everystep
      push!(W.t,t)
      push!(W.W,copy(out1))
      if W.Z != nothing
        push!(W.Z,copy(out2))
      end
    end
  else # Bridge
    i = searchsortedfirst(W.t,t)
    if t == W.t[i]
      out1 .= W.W[i]
      if W.Z != nothing
        out2 .= W.Z[i]
      end
    else
      W0,Wh = W.W[i-1],W.W[i]
      if W.Z != nothing
        Z0,Zh = W.Z[i-1],W.Z[i]
      end
      h = W.t[i]-W.t[i-1]
      q = (t-W.t[i-1])/h
      W.bridge(out1,W,W0,Wh,q,h,u,p,t,W.rng)
      out1 .+= (1-q)*W0
      if W.Z != nothing
        W.bridge(out2,W,Z0,Zh,q,h,u,p,t,W.rng)
        out2 .+= (1-q)*Z0
      end
      W.curW .= out1
      if W.save_everystep
        insert!(W.W,i,copy(out1))
        insert!(W.t,i,t)
      end
      if W.Z != nothing
        W.curZ .= out2
        if W.save_everystep
          insert!(W.Z,i,copy(out2))
        end
      end
    end
  end
end


text:
function package_show(connection::CKANConnection,
                      id::String; querryres::Bool=true)
    url = connection.url
    querry = package_show(url, id)
    result = querry["result"]
    resources_dicts = result["resources"]
    resources_ids = [x["id"] for x in resources_dicts]

    if querryres
        return CKANDataset(connection, result["id"], result,
                           resource_show.([connection], resources_ids))
    else
        return CKANDataset(connection, result["id"], result, resources_ids)
    end

end

function package_show(url::String, id::String)
    parameters = Pair("id", id)
    return querry_api(url, "package_show", parameters)
end

function resource_show(connection::CKANConnection, id::String)
    url = connection.url
    querry = resource_show(url, id)
    result = querry["result"]
    return CKANResource(connection, result["id"], result["package_id"], result)
end

function resource_show(url::String, id::String)
    parameters = Pair("id", id)
    return querry_api(url, "resource_show", parameters)
end

function package_list(connection::CKANConnection;
                      querryres::Bool=true, limit::Int64=20)
    url = connection.url
    querry = package_list(url)
    package_ids = querry["result"]

    if length(package_ids) > limit
        @info "$(length(package_ids)) datasets, querying the first $limit"
        package_ids = package_ids[1:limit]
    end

    return package_show.([connection], package_ids, querryres=querryres)
end

function package_list(url::String)
    return querry_api(url, "package_list", ("" => ""))
end

function package_search(connection::CKANConnection, parameters::Pair...;
                        querrypkg::Bool=true, querryres::Bool=true)
    url = connection.url
    querry = package_search(url, parameters...)

    result = querry["result"]
    count = result["count"]

    dicts = result["results"]
    ids = [x["id"] for x in dicts]

    @info "Search returned $count results"
    if querrypkg
        return package_show.([connection], ids, querryres=querryres)
    else
        return result
    end
end

function package_search(url::String, parameters::Pair...)
    return querry_api(url, "package_search", parameters...)
end

# function resource_search(connection::CKANConnection, parameters::Pair)
# end

# function resource_search(url::String, parameters::Pair)
#     return querry_api(url, "resource_search", parameters)
# end


text:
__precompile__()

module ElectronicStructure

include("m_Utility.jl")
include("m_LagrangeFunction.jl")
include("m_PsPot.jl")

end


text:
abstract type AbstractBoundaryOperator end

struct BoundaryOperator <: AbstractBoundaryOperator
    operators::Any
    facetooperator::Any
    ncells::Any
    function BoundaryOperator(operators, facetooperator)
        nphase, nface, ncells = size(facetooperator)
        @assert nphase == 2
        @assert nface == 4
        @assert all(facetooperator .>= 0)
        @assert all(facetooperator .<= length(operators))
        new(operators, facetooperator, ncells)
    end
end

function boundary_mass_operator(basis, facequads, cutmesh, onboundary, penalty)
    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    cellsign = cell_sign(cutmesh)
    ncells = length(cellsign)
    facedetjac = face_determinant_jacobian(cutmesh)
    dim = dimension(basis)

    refquads = reference_face_quadrature(facequads)
    operators = [
        penalty * mass_matrix(basis, q, detjac, dim)
        for (q, detjac) in zip(refquads, facedetjac)
    ]
    facetooperator = zeros(Int, 2, nfaces, ncells)

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cellsign[cellid]
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 1
                        facetooperator[1, faceid, cellid] = faceid
                    elseif s == -1
                        facetooperator[2, faceid, cellid] = faceid
                    elseif s == 0
                        pquad = facequads[+1, faceid, cellid]
                        positiveop =
                            penalty * mass_matrix(basis, pquad, facedetjac[faceid], dim)
                        push!(operators, positiveop)
                        facetooperator[1, faceid, cellid] = length(operators)

                        nquad = facequads[-1, faceid, cellid]
                        negativeop =
                            penalty * mass_matrix(basis, nquad, facedetjac[faceid], dim)
                        push!(operators, negativeop)
                        facetooperator[2, faceid, cellid] = length(operators)
                    else
                        error("Expected cellsign ∈ {-1,0,1}, got cellsign = $s")
                    end
                end
            end
        end
    end
    BoundaryOperator(operators, facetooperator)
end

function boundary_traction_operator(basis, facequads, stiffness, cutmesh, onboundary)

    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    cellsign = cell_sign(cutmesh)
    ncells = length(cellsign)

    cellmap = cell_map(cutmesh, 1)
    facedetjac = face_determinant_jacobian(cellmap)
    dim = dimension(basis)

    refnormals = reference_face_normals()
    refquads = reference_face_quadrature(facequads)
    operators1 = [
        face_traction_operator(basis, q, n, stiffness[+1], detjac, cellmap)
        for (q, n, detjac) in zip(refquads, refnormals, facedetjac)
    ]
    operators2 = [
        face_traction_operator(basis, q, n, stiffness[-1], detjac, cellmap)
        for (q, n, detjac) in zip(refquads, refnormals, facedetjac)
    ]

    operators = vcat(operators1, operators2)
    facetooperator = zeros(Int, 2, nfaces, ncells)

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cellsign[cellid]
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 1
                        facetooperator[1, faceid, cellid] = faceid
                    elseif s == -1
                        facetooperator[2, faceid, cellid] = nfaces + faceid
                    elseif s == 0
                        pquad = facequads[+1, faceid, cellid]
                        positiveop = face_traction_operator(
                            basis,
                            pquad,
                            refnormals[faceid],
                            stiffness[1],
                            facedetjac[faceid],
                            cellmap,
                        )
                        push!(operators, positiveop)
                        facetooperator[1, faceid, cellid] = length(operators)

                        nquad = facequads[-1, faceid, cellid]
                        negativeop = face_traction_operator(
                            basis,
                            nquad,
                            refnormals[faceid],
                            stiffness[-1],
                            facedetjac[faceid],
                            cellmap,
                        )
                        push!(operators, negativeop)
                        facetooperator[2, faceid, cellid] = length(operators)
                    else
                        error("Expected cellsign ∈ {-1,0,+1}, got cellsign = $s")
                    end
                end
            end
        end
    end
    BoundaryOperator(operators, facetooperator)
end

function boundary_displacement_rhs(bcfunc, basis, facequads, cutmesh, onboundary, penalty)
    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    cellsign = cell_sign(cutmesh)
    facedetjac = face_determinant_jacobian(cell_map(cutmesh, 1))
    ncells = length(cellsign)

    vectors = []
    facetovectors = zeros(Int, 2, nfaces, ncells)

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cellsign[cellid]
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    @assert s == +1 || s == 0 || s == -1
                    if s == +1 || s == 0
                        cm = cell_map(cutmesh,+1,cellid)
                        pquad = facequads[+1, faceid, cellid]
                        rhs =
                            penalty *
                            linear_form(bcfunc, basis, pquad, cm, facedetjac[faceid])
                        push!(vectors, rhs)
                        facetovectors[1, faceid, cellid] = length(vectors)
                    end
                    if s == -1 || s == 0
                        cm = cell_map(cutmesh,-1,cellid)
                        nquad = facequads[-1, faceid, cellid]
                        rhs =
                            penalty *
                            linear_form(bcfunc, basis, nquad, cm, facedetjac[faceid])
                        push!(vectors, rhs)
                        facetovectors[2, faceid, cellid] = length(vectors)
                    end
                end
            end
        end
    end
    BoundaryOperator(vectors, facetovectors)
end

struct BoundaryComponentOperator <: AbstractBoundaryOperator
    operators::Any
    facetooperator::Any
    ncells::Any
    component::Any
    function BoundaryComponentOperator(operators, facetooperator, component)
        nphase, nface, ncells = size(facetooperator)
        @assert nphase == 2
        @assert nface == 4
        @assert all(facetooperator .>= 0)
        @assert all(facetooperator .<= length(operators))
        new(operators, facetooperator, ncells, component)
    end
end

function boundary_mass_component_operator(
    basis,
    facequads,
    cutmesh,
    onboundary,
    component,
    penalty,
)
    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    cellsign = cell_sign(cutmesh)
    ncells = length(cellsign)
    facedetjac = face_determinant_jacobian(cell_map(cutmesh, 1))
    dim = dimension(basis)

    refquads = reference_face_quadrature(facequads)
    operators = [
        penalty * component_mass_matrix(basis, q, component, detjac)
        for (q, detjac) in zip(refquads, facedetjac)
    ]
    facetooperator = zeros(Int, 2, nfaces, ncells)

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cellsign[cellid]
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 1
                        facetooperator[1, faceid, cellid] = faceid
                    elseif s == -1
                        facetooperator[2, faceid, cellid] = faceid
                    elseif s == 0
                        pquad = facequads[+1, faceid, cellid]
                        positiveop =
                            penalty * component_mass_matrix(
                                basis,
                                pquad,
                                component,
                                facedetjac[faceid],
                            )
                        push!(operators, positiveop)
                        facetooperator[1, faceid, cellid] = length(operators)

                        nquad = facequads[-1, faceid, cellid]
                        negativeop =
                            penalty * component_mass_matrix(
                                basis,
                                nquad,
                                component,
                                facedetjac[faceid],
                            )
                        push!(operators, negativeop)
                        facetooperator[2, faceid, cellid] = length(operators)
                    else
                        error("Expected cellsign ∈ {-1,0,1}, got cellsign = $s")
                    end
                end
            end
        end
    end
    BoundaryComponentOperator(operators, facetooperator, component)
end

function boundary_traction_component_operator(
    basis,
    facequads,
    stiffness,
    cutmesh,
    onboundary,
    component,
)

    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    cellsign = cell_sign(cutmesh)
    ncells = length(cellsign)

    cellmap = cell_map(cutmesh, 1)
    facedetjac = face_determinant_jacobian(cellmap)
    dim = dimension(basis)

    refnormals = reference_face_normals()
    refquads = reference_face_quadrature(facequads)
    operators1 = [
        face_traction_component_operator(
            basis,
            q,
            component,
            n,
            stiffness[+1],
            detjac,
            cellmap,
        ) for (q, n, detjac) in zip(refquads, refnormals, facedetjac)
    ]
    operators2 = [
        face_traction_component_operator(
            basis,
            q,
            component,
            n,
            stiffness[-1],
            detjac,
            cellmap,
        ) for (q, n, detjac) in zip(refquads, refnormals, facedetjac)
    ]

    operators = vcat(operators1, operators2)
    facetooperator = zeros(Int, 2, nfaces, ncells)

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cellsign[cellid]
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    if s == 1
                        facetooperator[1, faceid, cellid] = faceid
                    elseif s == -1
                        facetooperator[2, faceid, cellid] = nfaces + faceid
                    elseif s == 0
                        pquad = facequads[+1, faceid, cellid]
                        positiveop = face_traction_component_operator(
                            basis,
                            pquad,
                            component,
                            refnormals[faceid],
                            stiffness[+1],
                            facedetjac[faceid],
                            cellmap,
                        )
                        push!(operators, positiveop)
                        facetooperator[1, faceid, cellid] = length(operators)

                        nquad = facequads[-1, faceid, cellid]
                        negativeop = face_traction_component_operator(
                            basis,
                            nquad,
                            component,
                            refnormals[faceid],
                            stiffness[-1],
                            facedetjac[faceid],
                            cellmap,
                        )
                        push!(operators, negativeop)
                        facetooperator[2, faceid, cellid] = length(operators)
                    else
                        error("Expected cellsign ∈ {-1,0,+1}, got cellsign = $s")
                    end
                end
            end
        end
    end
    BoundaryComponentOperator(operators, facetooperator, component)
end

function boundary_displacement_component_rhs(
    bcfunc,
    basis,
    facequads,
    cutmesh,
    onboundary,
    component,
    penalty,
)

    cellids = findall(is_boundary_cell(cutmesh))
    cellconnectivity = cell_connectivity(cutmesh)
    facemidpoints = reference_face_midpoints()
    nfaces = length(facemidpoints)
    cellsign = cell_sign(cutmesh)
    facedetjac = face_determinant_jacobian(cell_map(cutmesh, 1))
    ncells = length(cellsign)

    vectors = []
    facetovectors = zeros(Int, 2, nfaces, ncells)

    for cellid in cellids
        cellmap = cell_map(cutmesh, cellid)
        s = cellsign[cellid]
        for faceid = 1:nfaces
            if cellconnectivity[faceid, cellid] == 0
                if onboundary(cellmap(facemidpoints[faceid]))
                    @assert s == +1 || s == 0 || s == -1
                    if s == +1 || s == 0
                        pquad = facequads[+1, faceid, cellid]
                        rhs =
                            penalty * component_linear_form(
                                bcfunc,
                                basis,
                                pquad,
                                component,
                                cellmap,
                                facedetjac[faceid],
                            )
                        push!(vectors, rhs)
                        facetovectors[1, faceid, cellid] = length(vectors)
                    end
                    if s == -1 || s == 0
                        nquad = facequads[-1, faceid, cellid]
                        rhs =
                            penalty * component_linear_form(
                                bcfunc,
                                basis,
                                nquad,
                                component,
                                cellmap,
                                facedetjac[faceid],
                            )
                        push!(vectors, rhs)
                        facetovectors[2, faceid, cellid] = length(vectors)
                    end
                end
            end
        end
    end
    BoundaryComponentOperator(vectors, facetovectors, component)
end

function Base.getindex(bm::BO, s, faceid, cellid) where {BO<:AbstractBoundaryOperator}
    phaseid = cell_sign_to_row(s)
    idx = bm.facetooperator[phaseid, faceid, cellid]
    idx > 0 ||
        error("Expected idx > 0, got idx = $idx. Check if [s,faceid,cellid] =[$s,$faceid,$cellid] has a valid operator")
    return bm.operators[idx]
end

function Base.show(io::IO, bm::BoundaryOperator)
    ncells = bm.ncells
    noperators = length(bm.operators)
    str = "BoundaryOperator\n\tNum. Cells: $ncells\n\tNum. Unique Operators: $noperators"
    print(io, str)
end

function Base.show(io::IO, bm::BoundaryComponentOperator)
    ncells = bm.ncells
    noperators = length(bm.operators)
    component = bm.component
    str = "BoundaryComponentOperator\n\tNum. Cells: $ncells\n\tNum. Unique Operators: $noperators\n\tComponent: $component"
    print(io, str)
end

function number_of_cells(bm::BO) where {BO<:AbstractBoundaryOperator}
    return bm.ncells
end

function has_operator(bm::BO, s, faceid, cellid) where {BO<:AbstractBoundaryOperator}
    @assert s == -1 || s == +1
    phaseid = s == +1 ? 1 : 2
    return bm.facetooperator[phaseid, faceid, cellid] > 0
end

abstract type AbstractDisplacementCondition end

struct DisplacementCondition <: AbstractDisplacementCondition
    tractionoperator::Any
    massoperator::Any
    displacementrhs::Any
    penalty::Any
    ncells::Any
    function DisplacementCondition(
        tractionoperator::BoundaryOperator,
        massoperator::BoundaryOperator,
        displacementrhs::BoundaryOperator,
        penalty,
    )
        ncells = number_of_cells(massoperator)
        @assert number_of_cells(tractionoperator) == ncells
        @assert number_of_cells(displacementrhs) == ncells
        new(tractionoperator, massoperator, displacementrhs, penalty, ncells)
    end
end

function DisplacementCondition(
    displacementfunc,
    basis,
    facequads,
    stiffness,
    cutmesh,
    onboundary,
    penalty,
)
    massoperator = boundary_mass_operator(basis, facequads, cutmesh, onboundary, penalty)
    tractionoperator =
        boundary_traction_operator(basis, facequads, stiffness, cutmesh, onboundary)
    displacementrhs = boundary_displacement_rhs(
        displacementfunc,
        basis,
        facequads,
        cutmesh,
        onboundary,
        penalty,
    )
    return DisplacementCondition(tractionoperator, massoperator, displacementrhs, penalty)
end

struct DisplacementComponentCondition <: AbstractDisplacementCondition
    tractionoperator::Any
    massoperator::Any
    displacementrhs::Any
    penalty::Any
    component::Any
    ncells::Any
    function DisplacementComponentCondition(
        tractionoperator::BoundaryComponentOperator,
        massoperator::BoundaryComponentOperator,
        displacementrhs::BoundaryComponentOperator,
        penalty,
        component,
    )
        ncells = number_of_cells(massoperator)
        @assert number_of_cells(tractionoperator) == ncells
        @assert number_of_cells(displacementrhs) == ncells
        new(tractionoperator, massoperator, displacementrhs, penalty, component, ncells)
    end
end

function DisplacementComponentCondition(
    displacementfunc,
    basis,
    facequads,
    stiffness,
    cutmesh,
    onboundary,
    component,
    penalty,
)

    massoperator = boundary_mass_component_operator(
        basis,
        facequads,
        cutmesh,
        onboundary,
        component,
        penalty,
    )
    tractionoperator = boundary_traction_component_operator(
        basis,
        facequads,
        stiffness,
        cutmesh,
        onboundary,
        component,
    )
    displacementrhs = boundary_displacement_component_rhs(
        displacementfunc,
        basis,
        facequads,
        cutmesh,
        onboundary,
        component,
        penalty,
    )

    return DisplacementComponentCondition(
        tractionoperator,
        massoperator,
        displacementrhs,
        penalty,
        component,
    )
end

function traction_operator(
    dispcondition::DC,
    s,
    faceid,
    cellid,
) where {DC<:AbstractDisplacementCondition}
    return dispcondition.tractionoperator[s, faceid, cellid]
end

function mass_operator(
    dispcondition::DC,
    s,
    faceid,
    cellid,
) where {DC<:AbstractDisplacementCondition}
    return dispcondition.massoperator[s, faceid, cellid]
end

function displacement_rhs(
    dispcondition::DC,
    s,
    faceid,
    cellid,
) where {DC<:AbstractDisplacementCondition}
    return dispcondition.displacementrhs[s, faceid, cellid]
end

function Base.show(io::IO, dispcondition::DisplacementCondition)
    ncells = dispcondition.ncells
    penalty = dispcondition.penalty
    str = "DisplacementCondition\n\tNum. Cells: $ncells\n\tPenalty: $penalty"
    print(io, str)
end

function Base.show(io::IO, dispcondition::DisplacementComponentCondition)
    ncells = dispcondition.ncells
    penalty = dispcondition.penalty
    component = dispcondition.component
    str = "DisplacementComponentCondition\n\tNum. Cells: $ncells\n\tPenalty: $penalty\n\tComponent: $component"
    print(io, str)
end

function has_operator(
    dispcondition::DC,
    s,
    faceid,
    cellid,
) where {DC<:AbstractDisplacementCondition}
    flag = has_operator(dispcondition.tractionoperator, s, faceid, cellid)
    @assert has_operator(dispcondition.massoperator, s, faceid, cellid) == flag
    @assert has_operator(dispcondition.displacementrhs, s, faceid, cellid) == flag
    return flag
end


text:
#=
BFSParentsMultiSourceTest:
- Julia version: 1.6.2
- Author: samuel
- Date: 2021-09-11
=#

using SuiteSparseGraphBLAS
using SparseArrays
using LinearAlgebra
using Test

include("/Users/samuel/Documents/Jupyter/GraphBlasAlgorithmsFinal/src/BFSParentsMultiSource.jl")
S = GBMatrix(sparse([[1, 0, 0] [0, 0, 0] [0, 3, 0] [0, 0, 4] [0, 0, 0] [0, 0, 0] [0, 0, 0]]))
matrix =  GBMatrix(sparse([[0, 0, 0, 1, 0, 0, 0] [1, 0, 0, 0, 0, 0, 0] [0, 0, 0, 1, 0, 1, 1] [1, 0, 0, 0, 0, 0, 1] [0, 1, 0, 0, 0, 0, 1] [0, 0, 1, 0, 1, 0, 0] [0, 1, 0, 0, 0, 0, 0]]))
@test ms_bfsp(matrix, S, 7, 3) == [[0, nothing, 4] [1, nothing, 1] [4, 0, 4] [1, nothing, 0] [2, nothing, 2] [3, 3, 3] [2, nothing, 2]]


text:
using IterativeSolvers
using Test

@testset "Hessenberg" begin

    # Some well-conditioned Hessenberg matrix
    H1 = [
        1.19789 1.42354 -0.0371401  0.0118481 -0.0362113  0.00269463;
        1.46142 4.01953  0.890729  -0.0157701 -0.0300656 -0.0191307;
        0.0     1.08456  3.35179    0.941966   0.0439339 -0.072888;
        0.0     0.0      1.29071    3.1746     0.853378   0.0202058;
        0.0     0.0      0.0        1.32227    3.06086    1.18129;
        0.0     0.0      0.0        0.0        1.58682    2.99037;
        0.0     0.0      0.0        0.0        0.0        1.45345
    ]

    H2 = [
        1.14661+1.00222im  1.38611+0.0100307im  -0.0280623+0.00494128im    0.0205642-0.00867472im;
        1.43066+0.0im      4.11145+1.0144im       0.883563+0.00583816im  -0.00332255-0.0274278im;
            0.0+0.0im       1.0865+0.0im            3.2653+0.990342im       0.942434+0.00586481im;
            0.0+0.0im          0.0+0.0im           1.31146+0.0im             3.11491+1.01587im;
            0.0+0.0im          0.0+0.0im               0.0+0.0im             1.42175+0.0im
    ]

    for H = (H1, H2)
        T = eltype(H)

        # Fist standard basis vector as rhs
        rhs = [i == 1 ? one(T) : zero(T) for i = 1 : size(H, 1)]

        # Compare \ against the optimized version.
        solution_with_residual = copy(rhs)
        ldiv!(IterativeSolvers.FastHessenberg(copy(H)), solution_with_residual)
        solution = H \ rhs

        # First part is the solution
        @test solution_with_residual[1 : size(H, 2)] ≈ H \ rhs

        # Last element is the residual
        @test abs(last(solution_with_residual)) ≈ norm(H * solution - rhs)
    end
end


text:
"""
```
checkbounds{VT,FT,FD,FO}(m::AbstractMesh{VT,Face{FD,FT,FO}})
```

Check the `Face` indices to ensure they are in the bounds of the vertex
array of the `AbstractMesh`.
"""
function Base.checkbounds{VT,FT,FD,FO}(m::AbstractMesh{VT,Face{FD,FT,FO}})
    isempty(faces(m)) && return true # nothing to worry about I guess

    # index max and min
    const i = one(FO) + FO # normalize face offset
    s = length(vertices(m)) + FO

    for face in faces(m)
        # I hope this is unrolled
        for elt in face
            i <= elt && elt <= s || return false
        end
    end
    return true
end


text:
# using Revise, Test
using PiecewiseDeterministicMarkovProcesses, LinearAlgebra, Random, DifferentialEquations, Sundials
	const PDMP = PiecewiseDeterministicMarkovProcesses

function F!(ẋ, xc, xd, parms, t)
	if mod(xd[1], 2) == 0
		ẋ[1] = 10xc[1]
	else
		ẋ[1] = -3xc[1]^2
	end
end

R(x) = 10.0 + x

function R!(rate, xc, xd, parms, t, issum::Bool)
	if issum == false
		rate[1] = R(xc[1])
		rate[2] = xd[2]
		return 0., 0.
	else
		r = R(xc[1]) + xd[2]
		return r, r
	end
end

function Rcst!(rate, xc, xd, parms, t, issum::Bool)
	if issum == false
		# rate[1] = R(xc[1])
		rate[2] = xd[2]
		return 0., 0.
	else
		r = xd[2]
		return r, r
	end
end

function Rvar!(rate, xc, xd, parms, t, issum::Bool)
	if issum == false
		rate[1] = R(xc[1])
		# rate[2] = xd[2]
		return 0., 0.
	else
		r = R(xc[1])
		return r, r
	end
end

xc0 = [1.0]
	xd0 = [0, 0]

	nu = [1 0;0 -1]
	parms = [.0]
	ti = 0.0
	tf = 100000.
	nj = 20

# test the different way to write rates
rate0 = zeros(2)
Rvar = PDMP.VariableRate(R!)
Rcmp = PDMP.CompositeRate(Rcst!, Rvar!)

# using BenchmarkTools
# @btime R!($rate0, $[10.0], $[3,0], $parms, 0., false)
# @btime Rvar!($rate0, $[10.0], $[3,0], $parms, 0., false)
# @btime Rcst!($rate0, $[10.0], $[3,0], $parms, 0., false)
# @btime $Rvar($rate0, $[10.0], $[3,0], $parms, 0., false)
# @btime $Rcmp($rate0, $[10.0], $[3,0], $parms, 0., false)
#
# @btime $Rvar($rate0, $[10.0], $[3,0], $parms, 0., true)
# @btime $Rcmp($rate0, $[10.0], $[3,0], $parms, 0., true)

out1 =   R!(rate0, [1.0], [2,0], parms, 0., true)
outv = Rvar(rate0, [1.0], [2,0], parms, 0., true)
@test out1 == outv
outc = Rcmp(rate0, [1.0], [2,0], parms, 0., true)
@test out1 == outc
outc = Rcmp(rate0, [10.0], [3,0], parms, 0., true)
out1 =   R!(rate0,  [1.0], [3,0], parms, 0., true)
@test out1 != outc
outc = Rcmp(rate0, [10.0], [3,0], parms, 0., false)
outc = Rcmp(rate0, [10.0], [3,0], parms, 0., true)
out1 =   R!(rate0, [10.0], [3,0], parms, 0., true)
@test out1 == outc

algo = CHV(Tsit5())

# using Plots

Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	res0 =  PDMP.solve(problem, algo; n_jumps = nj)
	# plot(res0.time, res0.xc[1,:], marker = :d)

# here the rate function is constant in between jumps
pbvar = PDMP.PDMPProblem(F!, Rvar, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	resvar = PDMP.solve(pbvar, algo; n_jumps = nj)
	@test resvar.time == res0.time
	# plot!(resvar.time, resvar.xc[1,:], label = "Variable")

# here the rate function is constant in between jumps
pbcmp = PDMP.PDMPProblem(F!, Rcmp, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	rescmp = PDMP.solve(pbcmp, algo; n_jumps = nj)
	@test rescmp.time == res0.time
	# plot!(rescmp.time, rescmp.xc[1,:], label = "Composite")

# using PrettyTables
# datat = hcat(res0.time, resvar.time, rescst.time)
# datax = hcat(res0.xc', resvar.xc', rescst.xc')
#
# PrettyTables.pretty_table(datat)
# 	PrettyTables.pretty_table(datax)

# test for the end time
tf = 0.6
Random.seed!(8)
	problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	res0 =  PDMP.solve(problem, algo; n_jumps = nj)
	# plot(res0.time, res0.xc[1,:], marker = :d)

pbvar = PDMP.PDMPProblem(F!, Rvar, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	resvar = PDMP.solve(pbvar, algo; n_jumps = nj)
	@test resvar.time == res0.time
	# plot!(rescmp.time, resvar.xc[1,:])

pbcmp = PDMP.PDMPProblem(F!, Rcmp, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	rescmp = PDMP.solve(pbcmp, algo; n_jumps = nj)
	@test rescmp.time == res0.time
	# plot!(rescmp.time, rescmp.xc[1,:])

# test for allocations
Random.seed!(8)
	res0 = PDMP.solve(problem, algo; n_jumps = 3nj, save_positions = (false, false))
	resvar = PDMP.solve(pbvar, algo; n_jumps = 3nj, save_positions = (false, false))
	rescmp = PDMP.solve(pbcmp, algo; n_jumps = 3nj, save_positions = (false, false))

	Random.seed!(8)
	res0 = @timed PDMP.solve(problem, algo; n_jumps = nj, save_positions = (false, false))
	Random.seed!(8)
	resvar = @timed PDMP.solve(pbvar, algo; n_jumps = nj, save_positions = (false, false))
	Random.seed!(8)
	rescmp = @timed PDMP.solve(pbcmp, algo; n_jumps = nj, save_positions = (false, false))
	alloc0 = res0[5].poolalloc
	allocvar = resvar[5].poolalloc
	alloccmp = rescmp[5].poolalloc
	Random.seed!(8)
	res0 = @timed PDMP.solve(problem, algo; n_jumps = 3nj, save_positions = (false, false))
	Random.seed!(8)
	resvar = @timed PDMP.solve(pbvar, algo; n_jumps = 3nj, save_positions = (false, false))
	Random.seed!(8)
	rescmp = @timed PDMP.solve(pbcmp, algo; n_jumps = 3nj, save_positions = (false, false))
	# @test res0[5].poolalloc != alloc0
	# @test resvar[5].poolalloc != allocvar
	# @test rescmp[5].poolalloc != alloccmp

println("--> Allocations with Composite struct")
	Random.seed!(8)
	pbcmp = PDMP.PDMPProblem(F!, Rcmp, nu, xc0, xd0, parms, (ti, 100000.))
	rescmp = @time PDMP.solve(pbcmp, algo; n_jumps = nj, save_positions = (false, false))
	# plot(rescmp.time, rescmp.xc[1,:], label = "Composite", marker=:d)
	Random.seed!(8)
	rescmp = @time PDMP.solve(pbcmp, algo; n_jumps = 2nj, save_positions = (false, false))
	# plot!(rescmp.time, rescmp.xc[1,:], label = "Composite", marker=:c)
	Random.seed!(8)
	rescmp = @time PDMP.solve(pbcmp, algo; n_jumps = 3nj, save_positions = (false, false))
	# plot!(rescmp.time, rescmp.xc[1,:], label = "Composite")

# test with different CompositeStruct made of all variables
tf = 100000.
problem = PDMP.PDMPProblem(F!, R!, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	res0 =  PDMP.solve(problem, algo; n_jumps = nj)

	Rcmp = CompositeRate(VariableRate(Rcst!), VariableRate(Rvar!))
	pbcmp = PDMP.PDMPProblem(F!, Rcmp, nu, xc0, xd0, parms, (ti, tf))
	Random.seed!(8)
	rescmp = PDMP.solve(pbcmp, algo; n_jumps = nj)
	@test rescmp.time == res0.time

println("--> Allocations with Composite struct made of VariableRates")
	Random.seed!(8)
	pbcmp = PDMP.PDMPProblem(F!, Rcmp, nu, xc0, xd0, parms, (ti, 100000.))
	rescmp = @time PDMP.solve(pbcmp, algo; n_jumps = nj, save_positions = (false, false))
	Random.seed!(8)
	rescmp = @time PDMP.solve(pbcmp, algo; n_jumps = 2nj, save_positions = (false, false))
	Random.seed!(8)
	rescmp = @time PDMP.solve(pbcmp, algo; n_jumps = 3nj, save_positions = (false, false))


text:
# include("AngularMomentum.jl")

abstract type Pseudopotentials end

mutable struct SphericalPseudopotentials <: Pseudopotentials
    ell::AngMom
    Q::AngMom

    L::AngMomVector # total angular momentum on sphere
    m::AngMomVector # == 2ell - L == relative angular momenta on sphere
    V::Vector{Float64} # the actual pseudopotentials
end

# use the expression in the thesis by Wooten: http://trace.tennessee.edu/cgi/viewcontent.cgi?article=2884&context=utk_graddiss
# have checked numerically that this is equivalent to the expression in Toke and Jain, PRL 96, 246805 (2006) [Eq. (6)]
# NB: this assumes chord distance separation on the sphere when defining the 2-body potential V(r1 - r2)
function SphericalPseudopotentials(ell::AngMom, Q::AngMom, Vk::Vector{Float64})
    @assert isangmom(ell) && isangmom(Q)
    @assert ell >= Q

    L = collect(0:2ell)
    V = zeros(Float64, length(L))
    for i in 1:length(L)
        Vtmp = 0.
        for k in 0:Int(2ell)
            Vtmp += Vk[k+1] * Wigner6j(L[i], ell, ell, k, ell, ell) * Wigner3j(ell, k, ell, -Q, 0, Q)^2
        end
        V[i] = Vtmp * (-1)^(2Q + L[i]) * (2ell + 1)^2
    end

    m = 2ell - L
    p = sortperm(m) # return L, m, and V according to increasing relative angular momentum m = 2ell - L

    return SphericalPseudopotentials(ell, Q, L[p], m[p], V[p])
end

# for pure Coulomb interaction, the parameters Vk = 1/(radius of sphere) = const. (in Coulomb units e^2/(ϵ * ell0))
function spherical_Coulomb_pseudopotentials(ell::AngMom, Q::AngMom, R::Float64)::SphericalPseudopotentials
    return SphericalPseudopotentials(ell, Q, fill(1/R, Int(2ell)+1))
end

# pseudopotentials for LL index n = 0, 1, ...
# work in units assuming that we are in an effective LLL problem with monopole strength ell (= Q)
# i.e., so that the radius of the sphere is sqrt(ell) (in units of the magnetic length ell0)
# using this function will put energies in Coulomb units e^2/(ϵ * ell0) for the effective LLL problem
function spherical_Coulomb_pseudopotentials(ell::AngMom, n::Int)::SphericalPseudopotentials
    return spherical_Coulomb_pseudopotentials(ell, ell - n, sqrt(ell))
end

function set_pseudopotentials!(pp::SphericalPseudopotentials, Vnew::Dict)
    for m in keys(Vnew)
        @assert isangmom(m)
        ind = find(pp.m .== m)
        if length(ind) == 0
            warn("m = $(m) not found in pp.m")
        end
        pp.V[ind] = convert(Float64, Vnew[m])
    end
    return pp
end

# dinosaur, but keep it in here for reference
function V_Coulomb(L::Int, Nphi::Int, R::Float64)
    @assert 0 <= L <= Nphi
    # the following closed-form expression blows up for large Q, so we use the 2nd expression (see notes)
    # return (2/R) * (binomial(Int(4Q - 2L), Int(2Q - L))
    #              *  binomial(Int(4Q + 2L + 2), Int(2Q + L + 1))
    #              /  binomial(Int(4Q + 2), Int(2Q + 1))^2)
    if L == 0
        return (1/R) * float((Nphi + 1)/(2Nphi + 1))
    else # L >= 1
        return (1/R) * (prod([(2*(Nphi + L + 1 - i) - 1)/(Nphi + L + 1 - i) for i in 0:L-1])
                     *  prod([(Nphi + 1 - i)/(2*(Nphi + 1 - i) - 1) for i in 0:L]))
    end
end

function build_Vmat(pp::SphericalPseudopotentials, CG_table::Array{Vector{Float64}, 3})
    Norb = length(pp.V)
    @assert size(CG_table)[1] == size(CG_table)[2] == size(CG_table)[3]
    @assert size(CG_table)[1] == Norb

    println("Building 2-body interaction matrix ...")
    tic()
    Vmat = zeros(Float64, (Norb, Norb, Norb, Norb))
    for alpha1 in 1:Norb
        println("Working on alpha1 = $(alpha1): $(alpha1)/$(Norb)")
        for alpha2 in 1:Norb
            for alpha3 in 1:Norb
                for alpha4 in 1:Norb
                    V = 0.
                    for L in 0:Int(2pp.ell)
                        for Mz in -L:L
                            @assert pp.L[end - L] == L
                            V += CG_table[alpha1, alpha2, L + 1][Mz + (L + 1)] * conj(CG_table[alpha3, alpha4, L + 1][Mz + (L + 1)]) * pp.V[end - L]
                        end
                    end
                    Vmat[alpha1, alpha2, alpha3, alpha4] = V
                end
            end
        end
    end
    toc()
    println()

    return Vmat
end

max_stars_repo_name: mishmash/QHED

max_issues_repo_name: mishmash/QHED

max_forks_repo_name: mishmash/QHED


text:
#-------------------------------------------------------------------------------
# Set of useful functions related to LinearModel struct.
#-------------------------------------------------------------------------------

export QuadraticModel

"""

Define the type for quadratic models.

    - 'n': model dimension.
    - 'm': number of interpolation points.
    - 'imin': index of the function interpolated.
    - 'kopt': index of the best point so far in 'Y'.
    - 'g': gradient of the model.
    - 'gopt': gradient of the model in 'xopt'.
    - 'hq': holds the explicit second derivatives of the quadratic model.
    - 'pq': holds parameters of the implicit second derivatives of the quadratic model.
    - 'xbase': origin of the sample set and center of the model.
    - 'xopt': best point so far (in terms of function values).
    - 'fval': set of the function values of the interpolation points.
    - 'dst': distances between 'xopt' and other interpolation points.
    - 'Y': set of interpolation points, shifted from the center of the sample set 'xbase'.
    - 'BMAT': holds the elements of 'Ξ', with the exception of its first column,
            and the elements of 'Υ', with the exception of its first row and column.
    - 'ZMAT': holds the elements of 'Z', from the factorization 'Ω = ZZ^T'.

"""
struct QuadraticModel <: AbstractModel

    n     :: Int64
    m     :: Int64
    imin  :: Ref{Int64}
    kopt  :: Ref{Int64}
    g     :: AbstractVector
    gopt  :: AbstractVector
    hq    :: AbstractVector
    pq    :: AbstractVector
    xbase :: AbstractVector
    xopt  :: AbstractVector
    fval  :: AbstractVector
    dst   :: AbstractVector
    Y     :: AbstractMatrix
    BMAT  :: AbstractMatrix
    ZMAT  :: AbstractMatrix

end

"""

    QuadraticModel(n::Int64, m::Int64)

Create a quadratic interpolation model for a function (not specified yet) of R^{n}, with 'm' interpolation points.

    - 'n': dimension of the model.
    - 'm': number of interpolation points.

Returns a model of QuadraticModel type.

"""
QuadraticModel(n::Int64, m::Int64) = QuadraticModel(
    n, m, Ref{Int64}(), Ref{Int64}(),
    zeros(Float64, n), zeros(Float64, n), zeros(Float64, convert(Int64, n * ( n + 1 ) / 2) ),
    zeros(Float64, m), zeros(Float64, n), zeros(Float64, n), zeros(Int64, m),
    zeros(Float64, m - 1), zeros(Float64, m - 1, n), zeros(Float64, n + m, n), zeros(Float64, m, m - n - 1) )

"""

    construct_model!(func_list::Array{Function, 1}, imin_idx::Int64,
                        δ::Float64, fbase::Float64, xbase::Vector{Float64},
                        ao::Vector{Float64}, bo::Vector{Float64},
                        model::QuadraticModel)

Constructs the model based in the given information.

    - 'func_list': list containing the functions that determine the objective
    function fmin.

    - 'imin_idx': index belonging to the set I_{min}('xbase').

    - 'δ': radius of the sample set.

    - 'fbase': objective function value in 'xbase'.

    - 'xbase': n-dimensional vector (origin of the sample set).

    - 'ao': n-dimensional vector with the shifted lower bounds.

    - 'bo': n-dimensional vector with the shifted upper bounds.

The function modifies the argument:

    - 'model': model of QuadraticModel type.

"""
function construct_model!(
                            func_list::Array{Function, 1},
                            imin_idx::Int64,
                            δ::Float64,
                            fbase::Float64,
                            xbase::Vector{Float64},
                            ao::Vector{Float64},
                            bo::Vector{Float64},
                            model::QuadraticModel
                            )

    # Saves information about the index i ∈ I_{min}('xbase') and f_{i}('xbase').
    model.imin[] = imin_idx
    model.fval[1] = fbase

    # Saves information about 'model.xbase'.
    @. model.xbase = xbase

    # Usefull constant and variables
    aux = 1.0 / δ ^ 2.0
    kopt = 0
    fopt = fbase

    for i = 1:( model.m - 1 )

        # Usefull variables
        j = i - model.n
        k = i + 1

        if i ≤ ( 2 * model.n )

            if i ≤ model.n

                if bo[i] == 0.0

                    α = - δ

                else

                    α = δ

                end

                model.Y[i, i] = α
                model.dst[i] = δ

            else

                α = model.Y[j, j]

                if ao[j] == 0.0

                    β = min( 2.0 * δ, bo[j] )

                elseif bo[j] == 0.0

                    β = max( - 2.0 * δ, ao[j] )

                else

                    β = - δ

                end

                model.Y[i, j] = β
                model.dst[i] = abs( β )

            end

        else

            tmp = div( i - model.n - 1, model.n )
            p_j = i - tmp * model.n - model.n
            q_j = p_j + tmp

            if q_j > model.n

                tmp = p_j
                p_j = q_j - model.n
                q_j = tmp

            end

            model.Y[i, p_j] = model.Y[p_j, p_j]
            model.Y[i, q_j] = model.Y[q_j, q_j]
            model.dst[i] = sqrt( model.Y[i, p_j] ^ 2.0 + model.Y[i, q_j] ^ 2.0 )

        end

        # Uses 'xbase' has workspace to reconstruct the sample points.
        reconstruct_original_point!( model, i, a, b, ao, bo, xbase )
        model.fval[k] = fi_eval( func_list, imin_idx, xbase )

        # Searches for the least function value to determine 'kopt'.
        if model.fval[k] < fopt

            kopt = i
            fopt = model.fval[k]

        end

        # Determines 'g', 'hq', 'BMAT' and 'ZMAT'.
        if i ≤ ( 2 * model.n )

            if i ≤ model.n

                model.g[i] = ( model.fval[k] - fbase ) / α

                if model.m < ( k + model.n )

                    model.BMAT[1, i] = - 1.0 / α
                    model.BMAT[k, i] = 1.0 / α
                    model.BMAT[model.m + i, i] = - 0.5 * δ ^ 2.0

                end

            else

                idx_h = convert( Int64, j * ( j + 1 ) / 2 )
                tmp = ( model.fval[k] - fbase ) / β
                diff = β - α
                model.hq[idx_h] = 2.0 * ( tmp - model.g[j] ) / diff
                model.g[j] = ( model.g[j] * β - tmp * α ) / diff

                # Reorders the sample points if necessary.
                if α * β < 0.0

                    if model.fval[k] < model.fval[k - model.n]

                        tmp = model.fval[k - model.n]
                        model.fval[k - model.n] = model.fval[k]
                        model.fval[k] = tmp

                        if kopt == i

                            kopt = i - model.n

                        end

                        model.Y[j, j] = β
                        model.Y[i, j] = α
                        model.dst[j] = abs( β )
                        model.dst[i] = abs( α )

                    end

                end

                model.BMAT[1, j] = - ( α + β ) / ( α * β )
                model.BMAT[k, j] = - 0.5 / model.Y[i - model.n, j]
                model.BMAT[k - model.n, j] = - model.BMAT[1, j] - model.BMAT[k, j]
                model.ZMAT[1, j] = sqrt( 2.0 ) / ( α * β )
                model.ZMAT[k, j] = sqrt( 0.5 ) / δ ^ 2.0
                model.ZMAT[k - model.n, j] = - model.ZMAT[1, j] - model.ZMAT[k, j]

            end

        else

            idx_h = convert( Int64, q_j * ( q_j - 1 ) / 2 + p_j )
            model.ZMAT[1, j] = aux
            model.ZMAT[k, j] = aux
            model.ZMAT[q_j + 1, j] = - aux
            model.ZMAT[p_j + 1, j] = - aux
            tmp = model.Y[i, q_j] * model.Y[i, p_j]
            model.hq[idx_h] = ( fbase - model.fval[q_j + 1] - model.fval[p_j + 1] + model.fval[k] ) / tmp

        end

    end

    # Reconstructs the point 'xopt'
    reconstruct_original_point!( model, kopt, a, b, ao, bo, xbase )
    @. model.xopt = xbase
    model.kopt[] = kopt

    # Computes 'gopt' and updates 'dst' if necessary
    if kopt == 0

        @. model.gopt = model.g

    else

        update_gopt!( model, true )

        for i=1:(model.m - 1)

            if i == kopt

                continue

            else

                model.dst[i] = norm( model.Y[i, :] - model.Y[kopt, :] )

            end

        end

    end

end

"""

    update_gopt!(model::QuadraticModel, first_call::Bool)

Computes the gradient of the model in 'xopt'.

    - 'model': model of QuadraticModel type.

    - 'first_call' boolean that indicates if the model was newly created.

"""
function update_gopt!(
                        model::QuadraticModel,
                        first_call::Bool
                        )

    # Computes H*d, where d = xopt - xbase, and store in 'gopt'.
    mul_hess_vec!( model, model.Y[model.kopt[], :], model.gopt )

    # Computes gopt = g + H * d
    @. model.gopt += model.g

    # Why? Search in BOBYQA's paper and implementation
    if !( first_call )

        for k=1:( model.m - 1 )

            tmp = model.pq[k] * dot( model.Y[k, :], model.Y[model.kopt[], :] )

            for i=1:model.n

                model.gopt[i] += tmp * Y[k, i]

            end

        end

    end

end

"""

    mul_hess_vec!(model::QuadraticModel, vec::Vector{Float64}, x::Vector{Float64})

Calculates the product between the model's Hessian and a vector, i.e., computes H * vec.

    - 'model': model of QuadraticModel type.

    - 'vec': n-dimensional vector.

The function modifies the argument:

    - 'x': n-dimensional vector (solution)

"""
function mul_hess_vec!(
                        model::QuadraticModel,
                        vec::Vector{Float64},
                        x::Vector{Float64}
                        )

    for j=1:model.n

        x[j] = 0.0

        for i=1:j

            idx_h += 1

            if i < j

                x[j] += model.hq[idx_h] * vec[i]

            end

            x[i] += model.hq[idx_h] * vec[j]

        end

    end

end

"""

    mul_vec_hess_vec!(model::QuadraticModel, v::Vector{Float64}, w::Vector{Float64})

Computes the product v^{T} * H * w, where H is the hessian of the model.

    - 'model': model of QuadraticModel type.

    - 'v': n-dimensional vector.

    - 'w': n-dimensional vector.

Return the value of the product.

"""
function mul_vec_hess_vec!(
                            model::QuadraticModel,
                            v::Vector{Float64},
                            w::Vector{Float64}
                            )

    aux = 0.0

    for j=1:model.n

        for i=1:j

            idx_h += 1

            if i < j

                aux += model.hq[idx_h] * w[i] * v[j]

            end

            aux += model.hq[idx_h] * w[j] * v[i]

        end

    end

    return aux

end

"""

    stationarity_measure(model::LinearModel, a::Vector{Float64},
                            b::Vector{Float64})

Calculates the stationarity measure π_{k} = || P_{Ω}( x_{k} - gopt_{k} ) - x_{k} ||.

    - 'model': model of QuadraticModel type.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

Returns the stationarity measure.

"""
function stationarity_measure(
                                model::QuadraticModel,
                                a::Vector{Float64},
                                b::Vector{Float64}
                                )

    aux = 0.0
    for i = 1:model.n

        aux += ( min( max( a[i], model.xopt[i] - model.gopt[i] ), b[i] ) - model.xopt[i] ) ^ 2.0

    end

    return sqrt(aux)

end

"""

    compute_active_set!(model::LinearModel, a::Vector{Float64},
                        b::Vector{Float64}, active_set::Vector{Bool})

Computes the set of active constraints.

    - 'model': model of QuadraticModel type.

    - 'a': n-dimensional vector with the lower bounds.

    - 'b': n-dimensional vector with the upper bounds.

The function modifies the argument:

    - 'active_set': boolean n-dimensional vector with the active constraints.

"""
function compute_active_set!(
                        model::QuadraticModel,
                        a::Vector{Float64},
                        b::Vector{Float64},
                        active_set::Vector{Bool}
                        )

    for i=1:model.n

        if ( model.xopt[i] == a[i] && model.gopt[i] ≥ 0.0 ) || ( model.xopt[i] == b[i] && model.gopt[i] ≤ 0.0 )

            active_set[i] = true

        else

            active_set[i] = false

        end

    end

end

max_stars_repo_name: aschwertner/LOWDER

max_issues_repo_name: aschwertner/LOWDER

max_forks_repo_name: aschwertner/LOWDER


text:
## -------------------------------------------------------------------
function ToyModel3D()
    net = Dict()

    # glyc
    # https://en.wikipedia.org/wiki/Glycolysis
    # Glucose + 2 NAD+ + 2 ADP + 2 Pi → 2 pyruvate + 2 NADH + 2 H+ + 2 ATP

    # ppp
    # https://en.wikipedia.org/wiki/Pentose_phosphate_pathway
    # The overall reaction for this process is: Glucose 6-phosphate + 2 NADP+ + H2O → ribulose 5-phosphate + 2 NADPH + 2 H+ + CO.

    # resp
    # https://en.wikipedia.org/wiki/Oxidative_phosphorylation
    # 1/2 O2 + NADH + H+ → H2O + NAD+
    # O2 + 2*NADH + 2*H+ → 2*H2O + 2*NAD+
    # When one NADH is oxidized through the electron transfer chain, three ATPs are produced
    # When one NAPDH is oxidized through the electron transfer chain, two ATPs are produced
    # we use 2.5 as average
    # O2 + 2*NADH + 2*H+ → 2*H2O + 2*NAD+ 5*ATP

    # tac
    # https://en.wikipedia.org/wiki/Citric_acid_cycle
    # Products of the first turn of the cycle are one GTP (or ATP), three NADH, one FADH2 and two CO2.
    # Because two acetyl-CoA molecules are produced from each glucose molecule, two cycles are required per glucose molecule.
    # Therefore, at the end of two cycles, the products are: two GTP, six NADH, two FADH2, and four CO2.

    # ferm
    # http://ecolistudentportal.org/article_fermentation#_
    # Many fermenting bacteria including E. coli can convert acetyl-CoA to
    # acetate to generate ATP. The two key enzymes are phosphate acetyltransferase and
    # acetate kinase. This mini-pathway generates 1 mole of acetate and
    # 1 mole of ATP (by SLP) per acetyl-CoA consumed.
    # acetyl-CoA + Pi → acetyl-P + CoA
    # acetyl-P + ADP → acetate + ATP
    # AcCoa -> Ac + Atp

    # z biomass
    # Data From:
    # Varma, (1993): 2465–73. https://doi.org/10.1128/AEM.59.8.2465-2473.1993.
    # FIG 2
    # GBY = 0.068 # mmol GLC/gDW
    GBY = 1 / 0.068 # gDW/ mmol

    # feistGenomescaleMetabolicReconstruction2007 DOI: 10.1038/msb4100155
    # Using chemostat data for E. coli growing aerobically on
    # glucose (see Supplementary information), we estimated the
    # GAM and NGAM costs (Figure 3D). We found that an NGAM
    # value of 8.39 mmol ATP gDW1 h1 and a GAM value of
    # 59.81 mmol ATP gDW1 best fit the experimental data.
    GAM = 59.81 # EColi GAM 59.81 mmol ATP gDW1
    NGAM = 8.4 # EColi NGAM 8.4 mmolATP/gDW/h

    # uptakes limits
    AB = 1000.0 # abs bound
    # This model is bounded by the maximum rates found for EColi.
    # Data From:
    # Varma, (1993): 2465–73. https://doi.org/10.1128/AEM.59.8.2465-2473.1993.
    # Extract max exchages from FIG 3 to form the maximum polytope
    Vg = 20      # (glucose uptake limit) mmol/gDW h
    Vz = 0.5     # D < 0.5 below the acetate switch [23] Nature, 528(7580):99–104, 2015.

    net[:S] = [
    #             glyc   ppp    resp   tac    ferm   z     atpm   ua     ug     uo
    #= Glc   =#  -1.0   -1.0    0.0    0.0    0.0   -GBY   0.0    0.0   +1.0    0.0 ;  #= Glc   =#
    #= Atp   =#  +2.0    0.0   +5.0   +1.0   +1.0   -GAM  -NGAM   0.0    0.0    0.0 ;  #= Atp   =#
    #= Ac    =#   0.0    0.0    0.0    0.0   +1.0    0.0   0.0   +1.0    0.0    0.0 ;  #= Ac    =#
    #= NADH  =#  +4.0   +2.0   -2.0   +4.0    0.0    0.0   0.0    0.0    0.0    0.0 ;  #= NADH  =#
    #= AcCoa =#  +1.0   +1.0    0.0   -1.0   -1.0    0.0   0.0    0.0    0.0    0.0 ;  #= AcCoa =#
    #= Oxy   =#   0.0    0.0   -1.0    0.0    0.0    0.0   0.0    0.0    0.0   +1.0 ;  #= Oxy   =#
    ]

    net[:mets] = ["Glc", "Atp", "Ac", "NADH", "AcCoa", "Oxy"]
    net[:b] =    [0.0  , 0.0  , 0.0 , 0.0   , 0.0    , 0.0  ] # const exchanges

    net[:rxns] = ["glyc"  , "ppp"  , "resp" , "tac" , "ferm" , "z"   , "atpm", "ua"  ,  "ug" , "uo" ]
    net[:lb]   = [ 0.0    , 0.0    , 0.0    , 0.0   , 0.0    ,  0.0  ,  1.0  , -AB   ,  0.0  , 0.0  ]
    net[:ub]   = [ AB     , AB     , AB     , AB    , AB     ,  Vz   ,  1.0  ,  0.0  ,  Vg   ,  AB  ]
    net[:c]    = [ 0.0    , 0.0    , 0.0    , 0.0   , 0.0    ,  -1.0 ,  0.0  ,  0.0  ,  0.0  , 0.0  ]

    return reducebox!(MetNet(;net...))
end


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdspace

"""
    'feffery'_antdspace(;kwargs...)
    'feffery'_antdspace(children::Any;kwargs...)
    'feffery'_antdspace(children_maker::Function;kwargs...)

An AntdSpace component.

Keyword arguments:
- `children` (a list of or a singular dash component, string or number; optional): The content of the tab - will only be displayed if this tab is selected
- `id` (String; optional)
- `addSplitLine` (Bool; optional)
- `align` (a value equal to: 'start', 'end', 'center', 'baseline'; optional)
- `className` (String; optional)
- `direction` (a value equal to: 'vertical', 'horizontal'; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `size` (a value equal to: 'small', 'middle', 'large' | Real; optional)
- `style` (Dict; optional)
- `wrap` (Bool; optional)
"""
function 'feffery'_antdspace(; kwargs...)
        available_props = Symbol[:children, :id, :addSplitLine, :align, :className, :direction, :loading_state, :size, :style, :wrap]
        wild_props = Symbol[]
        return Component("'feffery'_antdspace", "AntdSpace", "feffery_antd_components", available_props, wild_props; kwargs...)
end

'feffery'_antdspace(children::Any; kwargs...) = 'feffery'_antdspace(;kwargs..., children = children)
'feffery'_antdspace(children_maker::Function; kwargs...) = 'feffery'_antdspace(children_maker(); kwargs...)

max_stars_repo_name: limmyforget/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
using JuLIP, QMMM2, SHIPs, PrettyTables
using JuLIP.MLIPs

module SWqmmm

using JuLIP, QMMM2, NeighbourLists, LinearAlgebra, SHIPs
using JuLIP.MLIPs: IPSuperBasis

pot() = StillingerWeber()
at_train() = bulk(:Si, cubic=true) * 2

function sketch()
   sw =  pot()
   rcut = cutoff(sw)
   at = at_train()
   return QMMM2.data_djEs_sketch(at, 1e-3, [rcut, rcut, rcut])
end

function training_set()
   sw = pot()
   sw_sketch = sketch()
   QMMM2.eval_dataset!(sw_sketch, sw; key="sw")
end

function get_basis(bo, deg; rinfact = 0.85, wL = 1.5)
   sw = pot()
   r0 = rnn(:Si)
   rcut = cutoff(sw)
   rin = rinfact * r0
   trans = PolyTransform(2, rnn(:Si))
   fcut = PolyCutoff1s(2, rin, rcut)
   shipB = SHIPBasis(SparseSHIP(bo, :Si, deg, wL), trans, fcut)
   basis = IPSuperBasis(OneBody(1.0), shipB)
end

function train_ship(bo, deg;
            wL = 1.5,
            rinfact = 0.85,
            weights = Dict("Es" => 100.0, "dEs" => 10.0, "d2Esh" => 1.0),
            rtol = 1e-5
         )
   trainset = training_set()
   basis = get_basis(bo, deg; rinfact=rinfact, wL=wL)
   🚢, fitinfo = QMMM2.lsqfit(basis, trainset, weights; rtol=rtol, key="sw")
end

end

# construct a training set
train = SWqmmm.training_set()

for deg = 2:6
   🚢, fitinfo = SWqmmm.train_ship(3, deg)
end


text:
empty!(server._documents)

LanguageServer.textDocument_didOpen_notification(LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem(uri"untitled:none", "julia", 0, "")), server, server.jr_endpoint)
@test LanguageServer.hasdocument(server, uri"untitled:none")

LanguageServer.textDocument_didOpen_notification(LanguageServer.DidOpenTextDocumentParams(LanguageServer.TextDocumentItem(uri"untitled:none", "julia", 0, "")), server, server.jr_endpoint)
@test LanguageServer.hasdocument(server, uri"untitled:none")

LanguageServer.textDocument_didSave_notification(LanguageServer.DidSaveTextDocumentParams(LanguageServer.TextDocumentIdentifier(uri"untitled:none"), ""), server, server.jr_endpoint)

LanguageServer.textDocument_didChange_notification(LanguageServer.DidChangeTextDocumentParams(LanguageServer.VersionedTextDocumentIdentifier(uri"untitled:none", 0), [LanguageServer.TextDocumentContentChangeEvent(missing, missing, "ran")]), server, server.jr_endpoint)

LanguageServer.textDocument_didClose_notification(LanguageServer.DidCloseTextDocumentParams(LanguageServer.TextDocumentIdentifier(uri"untitled:none")), server, server.jr_endpoint)
@test !LanguageServer.hasdocument(server, uri"untitled:none")


text:
using MethodAnalysis

"""
    is_atrisk_type(tt)

Given a Tuple-type signature (e.g., `Tuple{typeof(sum),Vector{Int}}`), determine whether this signature
is "at risk" for invalidation. Essentially it returns `true` if one or more arguments are of abstract type,
although there are prominent exceptions:

- Constructor calls with arbitrary argument types
- `convert(X, x)` where `isa(x, X)`
- `setindex!` and `push!` methods where the valtype is a subtype of the eltype (for AbstractDicts, likewise for the keytype)
- `getindex`, `length`, `isempty`, and `iterate` on any tuple

All of these are "allowed," meaning that they return `false`.
Moreover, some specific non-concrete argument types---like `Union`s of concrete types and `Function`---
do not trigger a return of `true`, although other at-risk argument types can lead to an overall `true` return
for the signature.
"""
function is_atrisk_type(@nospecialize(typ))
    # signatures like `convert(Vector, a)`, `foo(::Vararg{Synbol,N}) where N` do not seem to pose a problem
    isa(typ, TypeVar) && return false
    # isbits parameters are not a problem
    isa(typ, Type) || return false
    if isa(typ, UnionAll)
        typ = Base.unwrap_unionall(typ)
    end
    # Exclude signatures with Union{}
    typ === Union{} && return false
    isa(typ, Union) && return is_atrisk_type(typ.a) | is_atrisk_type(typ.b)
    # Type{T}: signatures like `convert(::Type{AbstractString}, ::String)` are not problematic
    typ <: Type && return false
    if typ <: Tuple && length(typ.parameters) >= 1
        p1 = typ.parameters[1]
        # Constructor calls are not themselves a problem (any `convert`s they trigger might be, but those are covered)
        isa(p1, Type) && p1 <: Type && return false
        # convert(::Type{T}, ::S) where S<:T is not problematic
        if p1 === typeof(Base.convert) || p1 === typeof(Core.convert)
            p2, p3 = typ.parameters[2], typ.parameters[3]
            if isa(p2, Type)
                p2 = Base.unwrap_unionall(p2)
                if isa(p2, DataType) && length(p2.parameters) === 1
                    T = p2.parameters[1]
                    if isa(T, TypeVar)
                        T = T.ub
                    end
                    isa(p3, Type) && isa(T, Type) && p3 <: T && return false
                end
            end
        # `getindex`, `length`, etc are OK for various Tuple{T1,T2,...}
        elseif p1 === typeof(Base.getindex) ||
               p1 === typeof(Base.length)  ||
               p1 === typeof(Base.isempty) ||
               p1 === typeof(Base.iterate) || p1 === typeof(Core.iterate)
            p2 = typ.parameters[2]
            if isa(p2, Type)
                p2 = Base.unwrap_unionall(p2)
                p2 <: Tuple && return false
            end
        # show(io::IO, x) is OK as long as typeof(x) is safe
        elseif p1 === typeof(Base.show) || p1 === typeof(Base.print) || p1 === typeof(Base.println)
            # is_atrisk_type(typ.parameters[2]) && return true
            for i = 3:length(typ.parameters)
                is_atrisk_type(typ.parameters[i]) && return true
            end
            return false
        # setindex!(a, x, idx) and push!(a, x) are safe if typeof(x) <: eltype(a)
        elseif (p1 === typeof(Base.setindex!) || p1 === typeof(Base.push!)) && length(typ.parameters) >= 3
            p2, p3 = typ.parameters[2], typ.parameters[3]
            if isconcretetype(p2)
                if p2 <: AbstractDict && length(typ.parameters) >= 4
                    p4 = typ.parameters[4]
                    p3 <: valtype(p2) && p4 <: keytype(p2) && return false
                    p2 <: IdDict && return false   # these are annotated @nospecialize
                else
                    p3 <: eltype(p2) && return false
                end
            end
        # vararg `string` methods should not be specialized
        elseif p1 === typeof(Base.string) && length(typ.parameters) > 2
            return false
        end
    end
    # Standard DataTypes
    isconcretetype(typ) && return false
    # ::Function args are excluded
    typ === Function && return false
    !isempty(typ.parameters) && (any(is_atrisk_type, typ.parameters) || return false)
    return true
end

@assert  is_atrisk_type(Tuple{typeof(==),Any,Any})
@assert  is_atrisk_type(Tuple{typeof(==),Symbol,Any})
@assert  is_atrisk_type(Tuple{typeof(==),Any,Symbol})
@assert !is_atrisk_type(Tuple{typeof(==),Symbol,Symbol})
@assert !is_atrisk_type(Tuple{typeof(convert),Type{Any},Any})
@assert !is_atrisk_type(Tuple{typeof(convert),Type{AbstractString},AbstractString})
@assert !is_atrisk_type(Tuple{typeof(convert),Type{AbstractString},String})
@assert  is_atrisk_type(Tuple{typeof(convert),Type{String},AbstractString})
@assert !is_atrisk_type(Tuple{typeof(convert),Type{Union{Int,Float32}},Int})
@assert !is_atrisk_type(Tuple{typeof(convert),Type{Union{Int,Float32}},Int32})
@assert  is_atrisk_type(Tuple{typeof(convert),Type{Union{Int,Float32}},Integer})
@assert !is_atrisk_type(Tuple{typeof(convert),Type{T} where T<:Union{Int,Float32},Int})
@assert !is_atrisk_type(Tuple{typeof(map),Function,Vector{Any}})
@assert !is_atrisk_type(Tuple{typeof(getindex),Dict{Union{String,Int},Any},Union{String,Int}})
@assert  is_atrisk_type(Tuple{typeof(getindex),Dict{Union{String,Int},Any},Any})
@assert !is_atrisk_type(Tuple{Type{BoundsError},Any,Any})
@assert  is_atrisk_type(Tuple{typeof(sin),Any})
@assert !is_atrisk_type(Tuple{typeof(length),Tuple{Any,Any}})
@assert  is_atrisk_type(Tuple{typeof(setindex!),Vector{Int},Any,Int})
@assert !is_atrisk_type(Tuple{typeof(setindex!),Vector{Any},Any,Int})
@assert  is_atrisk_type(Tuple{typeof(push!),Vector{Int},Any})
@assert !is_atrisk_type(Tuple{typeof(push!),Vector{Any},Any})

# Get the name of a method as written in the code. This strips keyword-method mangling.
function codename(sym::Symbol)
    symstr = String(sym)
    # Body methods
    m = match(r"^#(.*?)#\d+$", symstr)
    m !== nothing && return Symbol(only(m.captures))
    # kw methods
    m = match(r"^(.*?)##kw$", symstr)
    m !== nothing && return Symbol(only(m.captures))
    return sym
end

isexported(mi::Core.MethodInstance) = isdefined(Main, codename(mi.def.name))
getfunc(mi::Core.MethodInstance) = getfunc(mi.def)
getfunc(m::Method) = getfield(m.module, m.name)
nmethods(mi::Core.MethodInstance) = length(methods(getfunc(mi)))

# Test whether a module is Core.Compiler or inside it
# Methods there are protected from invalidation by other means
function fromcc(mod::Module)
    fn = fullname(mod)
    return length(fn) >= 2 && fn[1] === :Core && fn[2] === :Compiler
end

const mis = Dict{Method,Vector{Core.MethodInstance}}()
visit() do item
    if item isa Method && !fromcc(item.module)
        m = item
        mis[m] = methodinstances(m)
        return false
    end
    return true
end

# Count # of backedges for MethodInstances with abstract types
const becounter = Dict{Core.MethodInstance,Int}()
visit() do item
    if item isa Core.MethodInstance && !fromcc(item.def.module)
        if is_atrisk_type(item.specTypes)
            becounter[item] = length(all_backedges(item))
        end
        return false
    end
    return true
end

prs = sort!(collect(becounter); by=last)
open("/tmp/methdata_$VERSION.log", "w") do io
    for (mi, c) in prs
        c == 0 && continue
        println(io, mi.specTypes=>c)
    end
end

# Split into exported & private functions
miexp = Pair{Core.MethodInstance,Int}[]
mipriv = similar(miexp)
for (mi, c) in prs
    if isexported(mi)
        push!(miexp, mi=>c)
    else
        push!(mipriv, mi=>c)
    end
end


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

module TestBroadcastInternals

using Base.Broadcast: broadcast_indices, check_broadcast_indices,
                      check_broadcast_shape, newindex, _bcs, _bcsm
using Base: Test, OneTo

@test @inferred(_bcs((), (3,5), (3,5))) == (3,5)
@test @inferred(_bcs((), (3,1), (3,5))) == (3,5)
@test @inferred(_bcs((), (3,),  (3,5))) == (3,5)
@test @inferred(_bcs((), (3,5), (3,)))  == (3,5)
@test_throws DimensionMismatch _bcs((), (3,5), (4,5))
@test_throws DimensionMismatch _bcs((), (3,5), (3,4))
@test @inferred(_bcs((), (-1:1, 2:5), (-1:1, 2:5))) == (-1:1, 2:5)
@test @inferred(_bcs((), (-1:1, 2:5), (1, 2:5)))    == (-1:1, 2:5)
@test @inferred(_bcs((), (-1:1, 1),   (1, 2:5)))    == (-1:1, 2:5)
@test @inferred(_bcs((), (-1:1,),     (-1:1, 2:5))) == (-1:1, 2:5)
@test_throws DimensionMismatch _bcs((), (-1:1, 2:6), (-1:1, 2:5))
@test_throws DimensionMismatch _bcs((), (-1:1, 2:5), (2, 2:5))

@test @inferred(broadcast_indices(zeros(3,4), zeros(3,4))) == (OneTo(3),OneTo(4))
@test @inferred(broadcast_indices(zeros(3,4), zeros(3)))   == (OneTo(3),OneTo(4))
@test @inferred(broadcast_indices(zeros(3),   zeros(3,4))) == (OneTo(3),OneTo(4))
@test @inferred(broadcast_indices(zeros(3), zeros(1,4), zeros(1))) == (OneTo(3),OneTo(4))

check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,5))
check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,1))
check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3))
check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,5), zeros(3))
check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,5), 1)
check_broadcast_indices((OneTo(3),OneTo(5)), 5, 2)
@test_throws DimensionMismatch check_broadcast_indices((OneTo(3),OneTo(5)), zeros(2,5))
@test_throws DimensionMismatch check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,4))
@test_throws DimensionMismatch check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,4,2))
@test_throws DimensionMismatch check_broadcast_indices((OneTo(3),OneTo(5)), zeros(3,5), zeros(2))
check_broadcast_indices((-1:1, 6:9), 1)

check_broadcast_shape((-1:1, 6:9), (-1:1, 6:9))
check_broadcast_shape((-1:1, 6:9), (-1:1, 1))
check_broadcast_shape((-1:1, 6:9), (1, 6:9))
@test_throws DimensionMismatch check_broadcast_shape((-1:1, 6:9), (-1, 6:9))
@test_throws DimensionMismatch check_broadcast_shape((-1:1, 6:9), (-1:1, 6))

ci(x) = CartesianIndex(x)
@test @inferred(newindex(ci((2,2)), (true, true), (-1,-1)))   == ci((2,2))
@test @inferred(newindex(ci((2,2)), (true, false), (-1,-1)))  == ci((2,-1))
@test @inferred(newindex(ci((2,2)), (false, true), (-1,-1)))  == ci((-1,2))
@test @inferred(newindex(ci((2,2)), (false, false), (-1,-1))) == ci((-1,-1))
@test @inferred(newindex(ci((2,2)), (true,), (-1,-1)))   == ci((2,))
@test @inferred(newindex(ci((2,2)), (true,), (-1,)))   == ci((2,))
@test @inferred(newindex(ci((2,2)), (false,), (-1,))) == ci((-1,))
@test @inferred(newindex(ci((2,2)), (), ())) == ci(())

end

function as_sub(x::AbstractVector)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = view(y, 2:2:length(y))
    y[:] = x[:]
    y
end
function as_sub(x::AbstractMatrix)
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = view(y, 2:2:size(y,1), 2:2:size(y,2))
    for j=1:size(x,2)
        for i=1:size(x,1)
            y[i,j] = x[i,j]
        end
    end
    y
end
function as_sub{T}(x::AbstractArray{T,3})
    y = similar(x, eltype(x), tuple(([size(x)...]*2)...))
    y = view(y, 2:2:size(y,1), 2:2:size(y,2), 2:2:size(y,3))
    for k=1:size(x,3)
        for j=1:size(x,2)
            for i=1:size(x,1)
                y[i,j,k] = x[i,j,k]
            end
        end
    end
    y
end

bittest(f::Function, ewf::Function, a...) = (@test ewf(a...) == BitArray(broadcast(f, a...)))
n1 = 21
n2 = 32
n3 = 17
rb = 1:5

for arr in (identity, as_sub)
    @test broadcast(+, arr(eye(2)), arr([1, 4])) == [2 1; 4 5]
    @test broadcast(+, arr(eye(2)), arr([1  4])) == [2 4; 1 5]
    @test broadcast(+, arr([1  0]), arr([1, 4])) == [2 1; 5 4]
    @test broadcast(+, arr([1, 0]), arr([1  4])) == [2 5; 1 4]
    @test broadcast(+, arr([1, 0]), arr([1, 4])) == [2, 4]
    @test broadcast(+, arr([1, 0]), 2) == [3, 2]

    @test @inferred(arr(eye(2)) .+ arr([1, 4])) == arr([2 1; 4 5])
    @test arr(eye(2)) .+ arr([1  4]) == arr([2 4; 1 5])
    @test arr([1  0]) .+ arr([1, 4]) == arr([2 1; 5 4])
    @test arr([1, 0]) .+ arr([1  4]) == arr([2 5; 1 4])
    @test arr([1, 0]) .+ arr([1, 4]) == arr([2, 4])
    @test arr([1]) .+ arr([]) == arr([])

    A = arr(eye(2)); @test broadcast!(+, A, A, arr([1, 4])) == arr([2 1; 4 5])
    A = arr(eye(2)); @test broadcast!(+, A, A, arr([1  4])) == arr([2 4; 1 5])
    A = arr([1  0]); @test_throws DimensionMismatch broadcast!(+, A, A, arr([1, 4]))
    A = arr([1  0]); @test broadcast!(+, A, A, arr([1  4])) == arr([2 4])
    A = arr([1  0]); @test broadcast!(+, A, A, 2) == arr([3 2])

    @test arr([ 1    2])   .* arr([3,   4])   == [ 3 6; 4 8]
    @test arr([24.0 12.0]) ./ arr([2.0, 3.0]) == [12 6; 8 4]
    @test arr([1 2]) ./ arr([3, 4]) == [1/3 2/3; 1/4 2/4]
    @test arr([1 2]) .\ arr([3, 4]) == [3 1.5; 4 2]
    @test arr([3 4]) .^ arr([1, 2]) == [3 4; 9 16]
    @test arr(BitArray([true false])) .* arr(BitArray([true, true])) == [true false; true false]
    @test arr(BitArray([true false])) .^ arr(BitArray([false, true])) == [true true; true false]
    @test arr(BitArray([true false])) .^ arr([0, 3]) == [true true; true false]

    M = arr([11 12; 21 22])
    @test broadcast_getindex(M, eye(Int, 2).+1,arr([1, 2])) == [21 11; 12 22]
    @test_throws BoundsError broadcast_getindex(M, eye(Int, 2).+1,arr([1, -1]))
    @test_throws BoundsError broadcast_getindex(M, eye(Int, 2).+1,arr([1, 2]), [2])
    @test broadcast_getindex(M, eye(Int, 2).+1,arr([2, 1]), [1]) == [22 12; 11 21]

    A = arr(zeros(2,2))
    broadcast_setindex!(A, arr([21 11; 12 22]), eye(Int, 2).+1,arr([1, 2]))
    @test A == M
    broadcast_setindex!(A, 5, [1,2], [2 2])
    @test A == [11 5; 21 5]
    broadcast_setindex!(A, 7, [1,2], [1 2])
    @test A == fill(7, 2, 2)
    A = arr(zeros(3,3))
    broadcast_setindex!(A, 10:12, 1:3, 1:3)
    @test A == diagm(10:12)
    @test_throws BoundsError broadcast_setindex!(A, 7, [1,-1], [1 2])

    for (f, ewf) in (((==), (.==)),
                     ((<) , (.<) ),
                     ((!=), (.!=)),
                     ((<=), (.<=)))
        bittest(f, ewf, arr(eye(2)), arr([1, 4]))
        bittest(f, ewf, arr(eye(2)), arr([1  4]))
        bittest(f, ewf, arr([0, 1]), arr([1  4]))
        bittest(f, ewf, arr([0  1]), arr([1, 4]))
        bittest(f, ewf, arr([1, 0]), arr([1, 4]))
        bittest(f, ewf, arr(rand(rb, n1, n2, n3)), arr(rand(rb, n1, n2, n3)))
        bittest(f, ewf, arr(rand(rb,  1, n2, n3)), arr(rand(rb, n1,  1, n3)))
        bittest(f, ewf, arr(rand(rb,  1, n2,  1)), arr(rand(rb, n1,  1, n3)))
        bittest(f, ewf, arr(bitrand(n1, n2, n3)), arr(bitrand(n1, n2, n3)))
    end
end

r1 = 1:1
r2 = 1:5
ratio = [1,1/2,1/3,1/4,1/5]
@test r1.*r2 == [1:5;]
@test r1./r2 == ratio
m = [1:2;]'
@test m.*r2 == [1:5 2:2:10]
@test m./r2 ≈ [ratio 2ratio]
@test m./[r2;] ≈ [ratio 2ratio]

@test @inferred([0,1.2].+reshape([0,-2],1,1,2)) == reshape([0 -2; 1.2 -0.8],2,1,2)
rt = Base.return_types(.+, Tuple{Array{Float64, 3}, Array{Int, 1}})
@test length(rt) == 1 && rt[1] == Array{Float64, 3}
rt = Base.return_types(broadcast, Tuple{typeof(.+), Array{Float64, 3}, Array{Int, 3}})
@test length(rt) == 1 && rt[1] == Array{Float64, 3}
rt = Base.return_types(broadcast!, Tuple{Function, Array{Float64, 3}, Array{Float64, 3}, Array{Int, 1}})
@test length(rt) == 1 && rt[1] == Array{Float64, 3}

# f.(args...) syntax (#15032)
let x = [1,3.2,4.7], y = [3.5, pi, 1e-4], α = 0.2342
    @test sin.(x) == broadcast(sin, x)
    @test sin.(α) == broadcast(sin, α)
    @test sin.(3.2) == broadcast(sin, 3.2) == sin(3.2)
    @test factorial.(3) == broadcast(factorial, 3)
    @test atan2.(x, y) == broadcast(atan2, x, y)
    @test atan2.(x, y') == broadcast(atan2, x, y')
    @test atan2.(x, α) == broadcast(atan2, x, α)
    @test atan2.(α, y') == broadcast(atan2, α, y')
end

# issue 14725
let a = Number[2, 2.0, 4//2, 2+0im] / 2
    @test eltype(a) == Number
end
let a = Real[2, 2.0, 4//2] / 2
    @test eltype(a) == Real
end
let a = Real[2, 2.0, 4//2] / 2.0
    @test eltype(a) == Real
end

# issue 16164
let a = broadcast(Float32, [3, 4, 5])
    @test eltype(a) == Float32
end

# broadcasting scalars:
@test sin.(1) === broadcast(sin, 1) === sin(1)
@test (()->1234).() === broadcast(()->1234) === 1234

# issue #4883
@test isa(broadcast(tuple, [1 2 3], ["a", "b", "c"]), Matrix{Tuple{Int,String}})
@test isa(broadcast((x,y)->(x==1?1.0:x,y), [1 2 3], ["a", "b", "c"]), Matrix{Tuple{Real,String}})
let a = length.(["foo", "bar"])
    @test isa(a, Vector{Int})
    @test a == [3, 3]
end
let a = sin.([1, 2])
    @test isa(a, Vector{Float64})
    @test a ≈ [0.8414709848078965, 0.9092974268256817]
end

# PR #17300: loop fusion
@test (x->x+1).((x->x+2).((x->x+3).(1:10))) == collect(7:16)
let A = [sqrt(i)+j for i = 1:3, j=1:4]
    @test atan2.(log.(A), sum(A,1)) == broadcast(atan2, broadcast(log, A), sum(A, 1))
end
let x = sin.(1:10)
    @test atan2.((x->x+1).(x), (x->x+2).(x)) == broadcast(atan2, x+1, x+2) == broadcast(atan2, x.+1, x.+2)
    @test sin.(atan2.([x+1,x+2]...)) == sin.(atan2.(x+1,x+2))
    @test sin.(atan2.(x, 3.7)) == broadcast(x -> sin(atan2(x,3.7)), x)
    @test atan2.(x, 3.7) == broadcast(x -> atan2(x,3.7), x) == broadcast(atan2, x, 3.7)
end
# Use side effects to check for loop fusion.
let g = Int[]
    f17300(x) = begin; push!(g, x); x+2; end
    f17300.(f17300.(f17300.(1:3)))
    @test g == [1,3,5, 2,4,6, 3,5,7]
end
# fusion with splatted args:
let x = sin.(1:10), a = [x]
    @test cos.(x) == cos.(a...)
    @test atan2.(x,x) == atan2.(a..., a...) == atan2.([x, x]...)
    @test atan2.(x, cos.(x)) == atan2.(a..., cos.(x)) == broadcast(atan2, x, cos.(a...)) == broadcast(atan2, a..., cos.(a...))
    @test ((args...)->cos(args[1])).(x) == cos.(x) == ((y,args...)->cos(y)).(x)
end
@test atan2.(3,4) == atan2(3,4) == (() -> atan2(3,4)).()
# fusion with keyword args:
let x = [1:4;]
    f17300kw(x; y=0) = x + y
    @test f17300kw.(x) == x
    @test f17300kw.(x, y=1) == f17300kw.(x; y=1) == f17300kw.(x; [(:y,1)]...) == x .+ 1
    @test f17300kw.(sin.(x), y=1) == f17300kw.(sin.(x); y=1) == sin.(x) .+ 1
    @test sin.(f17300kw.(x, y=1)) == sin.(f17300kw.(x; y=1)) == sin.(x .+ 1)
end

# PR #17510: Fused in-place assignment
let x = [1:4;], y = x
    y .= 2:5
    @test y === x == [2:5;]
    y .= factorial.(x)
    @test y === x == [2,6,24,120]
    y .= 7
    @test y === x == [7,7,7,7]
    y .= factorial.(3)
    @test y === x == [6,6,6,6]
    f17510() = 9
    y .= f17510.()
    @test y === x == [9,9,9,9]
    y .-= 1
    @test y === x == [8,8,8,8]
    y .-= 1:4
    @test y === x == [7,6,5,4]
    x[1:2] .= 1
    @test y === x == [1,1,5,4]
    x[1:2] .+= [2,3]
    @test y === x == [3,4,5,4]
    x[:] .= 0
    @test y === x == [0,0,0,0]
    x[2:end] .= 1:3
    @test y === x == [0,1,2,3]
end
let a = [[4, 5], [6, 7]]
    a[1] .= 3
    @test a == [[3, 3], [6, 7]]
end
let d = Dict(:foo => [1,3,7], (3,4) => [5,9])
    d[:foo] .+= 2
    @test d[:foo] == [3,5,9]
    d[3,4] .-= 1
    @test d[3,4] == [4,8]
end
let identity = error, x = [1,2,3]
    x .= 1 # make sure it goes to broadcast!(Base.identity, ...), not identity
    @test x == [1,1,1]
end

# make sure scalars are inlined, which causes f.(x,scalar) to lower to a "thunk"
import Base.Meta: isexpr
@test isexpr(expand(:(f.(x,y))), :call)
@test isexpr(expand(:(f.(x,1))), :thunk)
@test isexpr(expand(:(f.(x,1.0))), :thunk)
@test isexpr(expand(:(f.(x,$π))), :thunk)

# PR 16988
@test Base.promote_op(+, Bool) === Int
@test isa(broadcast(+, [true]), Array{Int,1})

# issue #17304
let foo = [[1,2,3],[4,5,6],[7,8,9]]
    @test max.(foo...) == broadcast(max, foo...) == [7,8,9]
end

# Issue 17314
@test broadcast(x->log(log(log(x))), [1000]) == [log(log(log(1000)))]
let f17314 = x -> x < 0 ? false : x
    @test eltype(broadcast(f17314, 1:3)) === Int
    @test eltype(broadcast(f17314, -1:1)) === Integer
    @test eltype(broadcast(f17314, Int[])) === Union{Bool,Int}
end
let io = IOBuffer()
    broadcast(x->print(io,x), 1:5) # broadcast with side effects
    @test take!(io) == [0x31,0x32,0x33,0x34,0x35]
end

# Issue 18176
let f18176(a, b, c) = a + b + c
    @test f18176.(1.0:2, 3, 4) == f18176.(3.0, 1.0:2, 4.0) == broadcast(f18176, 3, 4, 1.0:2)
end

# Issue #17984
let A17984 = []
    @test isa(abs.(A17984), Array{Any,1})
end

# Issue #16966
@test parse.(Int, "1") == 1
@test parse.(Int, ["1", "2"]) == [1, 2]
@test trunc.((Int,), [1.2, 3.4]) == [1, 3]
@test abs.((1, -2)) == (1, 2)
@test broadcast(+, 1.0, (0, -2.0)) == (1.0,-1.0)
@test broadcast(+, 1.0, (0, -2.0), [1]) == [2.0, 0.0]
@test broadcast(*, ["Hello"], ", ", ["World"], "!") == ["Hello, World!"]

# Ensure that even strange constructors that break `T(x)::T` work with broadcast
immutable StrangeType18623 end
StrangeType18623(x) = x
StrangeType18623(x,y) = (x,y)
@test @inferred broadcast(StrangeType18623, 1:3) == [1,2,3]
@test @inferred broadcast(StrangeType18623, 1:3, 4:6) == [(1,4),(2,5),(3,6)]

@test typeof(Int.(Number[1, 2, 3])) === typeof((x->Int(x)).(Number[1, 2, 3]))

@test @inferred broadcast(CartesianIndex, 1:2) == [CartesianIndex(1), CartesianIndex(2)]
@test @inferred broadcast(CartesianIndex, 1:2, 3:4) == [CartesianIndex(1,3), CartesianIndex(2,4)]

# Issue 18622
@test @inferred muladd.([1.0], [2.0], [3.0])::Vector{Float64} == [5.0]
@test @inferred tuple.(1:3, 4:6, 7:9)::Vector{Tuple{Int,Int,Int}} == [(1,4,7), (2,5,8), (3,6,9)]

max_stars_repo_name: habemus-papadum/julia

max_issues_repo_name: habemus-papadum/julia

max_forks_repo_name: habemus-papadum/julia


text:
module IndexingTests

using Test
using InplaceArrays.TensorValues
using StaticArrays

a = (3,4,5,1)

v = VectorValue{4}(a)

@test eltype(v) == Int
@test eltype(typeof(v)) == Int

@test size(v) == (4,)
@test length(v) == 4

for (k,i) in enumerate(eachindex(v))
  @test v[i] == a[k]
end

t = TensorValue{2}(a)

@test size(t) == (2,2)
@test length(t) == 4

for (k,i) in enumerate(eachindex(t))
  @test t[i] == a[k]
end

@test t[2,1] == 4

@test t[2] == 4

for (k,ti) in enumerate(t)
  @test ti == a[k]
end

v = @SMatrix zeros(2,3)
w = MultiValue(v)
@test CartesianIndices(w) == CartesianIndices(v)
@test LinearIndices(w) == LinearIndices(v)

end # module IndexingTests