text:
module TestMonoidalDirectedWiringDiagrams
using Test

using Catlab.Theories, Catlab.WiringDiagrams

# Categorical interface
#######################

# Category
#---------

# Generators
A, B = Ports([:A]), Ports([:B])
f = singleton_diagram(Box(:f,A,B))
g = singleton_diagram(Box(:g,B,A))
@test nboxes(f) == 1
@test boxes(f) == [ Box(:f,A,B) ]
@test nwires(f) == 2

# Composition
@test nboxes(compose(f,g)) == 2
@test boxes(compose(f,g)) == [ Box(:f,A,B), Box(:g,B,A) ]
@test nwires(compose(f,g)) == 3

# Domains and codomains
@test dom(f) == Ports([:A])
@test codom(f) == Ports([:B])
@test dom(compose(f,g)) == Ports([:A])
@test codom(compose(f,g)) == Ports([:A])

# Associativity
@test compose(compose(f,g),f) == compose(f,compose(g,f))

# Identity
@test compose(id(dom(f)), f) == f
@test compose(f, id(codom(f))) == f

# Symmetric monoidal category
#----------------------------

# Domains and codomains
@test dom(otimes(f,g)) == otimes(dom(f),dom(g))
@test codom(otimes(f,g)) == otimes(codom(f),codom(g))

# Associativity and unit
X, Y = Ports([:A,:B]), Ports([:C,:D])
I = munit(Ports)
@test otimes(X,I) == X
@test otimes(I,X) == X
@test otimes(otimes(X,Y),X) == otimes(X,otimes(Y,X))
@test otimes(otimes(f,g),f) == otimes(f,otimes(g,f))

# Braiding
@test compose(braid(X,Y),braid(Y,X)) == id(otimes(X,Y))

# Permutations
W = otimes(X,Y)
@test permute(W, [1,2,3,4]) == id(W)
@test permute(W, [1,2,3,4], inverse=true) == id(W)
@test permute(W, [3,4,1,2]) == braid(X,Y)
@test permute(W, [3,4,1,2], inverse=true) == braid(Y,X)
@test_throws AssertionError permute(W, [1,2])

# Diagonals
#----------

# Basic composition
d = WiringDiagram(dom(f), otimes(codom(f),codom(f)))
fv1 = add_box!(d, first(boxes(f)))
fv2 = add_box!(d, first(boxes(f)))
add_wires!(d, [
  (input_id(d),1) => (fv1,1),
  (input_id(d),1) => (fv2,1),
  (fv1,1) => (output_id(d),1),
  (fv2,1) => (output_id(d),2),
])
@test compose(mcopy(dom(f)), otimes(f,f)) == d

# Domains and codomains
@test dom(mcopy(Ports([:A]))) == Ports([:A])
@test codom(mcopy(Ports([:A]))) == Ports([:A,:A])
@test dom(mcopy(Ports([:A,:B]),3)) == Ports([:A,:B])
@test codom(mcopy(Ports([:A,:B]),3)) == Ports([:A,:B,:A,:B,:A,:B])

# Associativity
A = Ports([:A])
@test compose(mcopy(A), otimes(id(A),mcopy(A))) == mcopy(A,3)
@test compose(mcopy(A), otimes(mcopy(A),id(A))) == mcopy(A,3)

# Commutativity
@test compose(mcopy(A), braid(A,A)) == mcopy(A)

# Unitality
@test compose(mcopy(A), otimes(id(A),delete(A))) == id(A)

# Cartesian categories
A = Ports{CartesianCategory}([:A])
@test mcopy(A) == mcopy(Ports([:A]))
@test delete(A) == delete(Ports([:A]))
@test_throws MethodError mmerge(A)
@test_throws MethodError create(A)

# Codiagonals
#------------

# Domains and codomains
@test dom(mmerge(Ports([:A]))) == Ports([:A,:A])
@test codom(mmerge(Ports([:A]))) == Ports([:A])
@test dom(mmerge(Ports([:A,:B]),3)) == Ports([:A,:B,:A,:B,:A,:B])
@test codom(mmerge(Ports([:A,:B]),3)) == Ports([:A,:B])

# Associativity
A = Ports([:A])
@test compose(otimes(id(A),mmerge(A)), mmerge(A)) == mmerge(A,3)
@test compose(otimes(mmerge(A),id(A)), mmerge(A)) == mmerge(A,3)

# Commutativity
@test compose(braid(A,A), mmerge(A)) == mmerge(A)

# Unitality
@test compose(otimes(id(A),create(A)), mmerge(A)) == id(A)

# Cocartesian categories
A = Ports{CocartesianCategory}([:A])
@test mmerge(A) == mmerge(Ports([:A]))
@test create(A) == create(Ports([:A]))
@test_throws MethodError mcopy(A)
@test_throws MethodError delete(A)

# Bidiagonals
#------------

# Monoidal categories with bidiagonals, and non-naturality of explicit
# representation.
A = Ports{MonoidalCategoryWithBidiagonals}([:A])
@test boxes(mcopy(A)) == [ Junction(:A,1,2) ]
@test boxes(mcopy(otimes(A,A))) == repeat([ Junction(:A,1,2) ], 2)
@test compose(create(A), mcopy(A)) != create(otimes(A,A))
@test compose(mmerge(A), delete(A)) != delete(otimes(A,A))

# Biproduct categories, and naturality of implicit representation.
A = Ports{BiproductCategory}([:A])
@test compose(create(A), mcopy(A)) == create(otimes(A,A))
@test compose(mmerge(A), delete(A)) == delete(otimes(A,A))

# Dagger category
#----------------

f = singleton_diagram(DaggerSymmetricMonoidalCategory, Box(:f,[:A],[:B]))
g = singleton_diagram(DaggerSymmetricMonoidalCategory, Box(:g,[:B],[:A]))

@test boxes(dagger(f)) == [ BoxOp{:dagger}(Box(:f,[:A],[:B])) ]

# Domain and codomain
@test dom(dagger(f)) == codom(f)
@test codom(dagger(f)) == dom(f)

# Functoriality
@test is_permuted_equal(dagger(compose(f,g)), compose(dagger(g),dagger(f)), [2,1])
@test dagger(otimes(f,g)) == otimes(dagger(f),dagger(g))

# Involutivity
@test dagger(dagger(f)) == f
@test dagger(dagger(compose(f,g))) == compose(f,g)
@test dagger(dagger(otimes(f,g))) == otimes(f,g)

# Compact closed category
#------------------------

### Duals

A, B = [ Ports{CompactClosedCategory}([sym]) for sym in [:A, :B] ]
I = munit(typeof(A))

@test boxes(dunit(A)) == [ Junction(:A, [], [PortOp{:dual}(:A), :A]) ]
@test boxes(dcounit(A)) == [ Junction(:A, [:A, PortOp{:dual}(:A)], []) ]

# Domains and codomains
@test dom(dunit(A)) == I
@test codom(dunit(A)) == otimes(dual(A),A)
@test dom(dcounit(A)) == otimes(A,dual(A))
@test codom(dcounit(A)) == I

# Coherence relations
@test is_permuted_equal(dunit(otimes(A,B)),
  compose(dunit(B), otimes(id(dual(B)), dunit(A), id(B))), [2,1])
@test is_permuted_equal(dcounit(otimes(A,B)),
  compose(otimes(id(A), dcounit(B), id(dual(A))), dcounit(A)), [2,1])

### Adjoint mates

f = singleton_diagram(CompactClosedCategory, Box(:f,[:A],[:B]))
g = singleton_diagram(CompactClosedCategory, Box(:g,[:B],[:A]))

@test boxes(mate(f)) == [ BoxOp{:mate}(Box(:f,[:A],[:B])) ]

# Domain and codomain
@test dom(mate(f)) == dual(codom(f))
@test codom(mate(f)) == dual(dom(f))

# Functoriality
@test is_permuted_equal(mate(compose(f,g)), compose(mate(g),mate(f)), [2,1])
@test is_permuted_equal(mate(otimes(f,g)), otimes(mate(g),mate(f)), [2,1])

# Involutivity
@test mate(mate(f)) == f
@test mate(mate(compose(f,g))) == compose(f,g)
@test mate(mate(otimes(f,g))) == otimes(f,g)

# Traced monoidal category
#-------------------------

const TracedMon = TracedMonoidalCategory
A, B, X, Y = [ Ports{TracedMon}([sym]) for sym in [:A,:B,:X,:Y] ]
I = munit(typeof(A))
f = singleton_diagram(TracedMon, Box(:f, [:X,:A], [:X,:B]))

# Domain and codomain
@test dom(trace(X, A, B, f)) == A
@test codom(trace(X, A, B, f)) == B

# Naturality
g = singleton_diagram(TracedMon, Box(:g, [:A], [:A]))
h = singleton_diagram(TracedMon, Box(:h, [:B], [:B]))
@test trace(X, compose(id(X)⊗g, f, id(X)⊗h)) == compose(g, trace(X,f), h)

# Stength, aka superposing
g = singleton_diagram(TracedMon, Box(:g, [:C], [:C]))
@test trace(X, otimes(f,g)) == otimes(trace(X,f), g)

# Symmetry sliding
f = singleton_diagram(TracedMon, Box(:f, [:X,:Y,:A], [:Y,:X,:B]))
@test trace(X⊗Y, compose(f, braid(Y,X)⊗id(B))) ==
  trace(Y⊗X, compose(braid(Y,X)⊗id(A), f))

# Vanishing
f = singleton_diagram(TracedMon, Box(:f, [:X,:Y,:A], [:X,:Y,:B]))
@test trace(X⊗Y, f) == trace(Y, trace(X, f))
f = singleton_diagram(TracedMon, Box(:f, [:A], [:B]))
@test trace(I, f) == f

# Yanking
@test trace(X, braid(X,X)) == id(X)

# Bicategory of relations
#------------------------

A, B = [ Ports{BicategoryRelations}([sym]) for sym in [:A, :B] ]
R = singleton_diagram(BicategoryRelations, Box(:R,[:A],[:B]))
S = singleton_diagram(BicategoryRelations, Box(:S,[:A],[:B]))

# Domains and codomains
@test dom(meet(R,S)) == A
@test codom(meet(R,S)) == B
@test dom(top(A,B)) == A
@test codom(top(A,B)) == B

# Units and counits
@test dunit(A) == merge_junctions(compose(create(A), mcopy(A)))
@test dcounit(A) == merge_junctions(compose(mmerge(A), delete(A)))

# Abelian bicategory of relations
#--------------------------------

A, B = [ Ports{AbelianBicategoryRelations}([sym]) for sym in [:A, :B] ]
R = singleton_diagram(AbelianBicategoryRelations, Box(:R,[:A],[:B]))
S = singleton_diagram(AbelianBicategoryRelations, Box(:S,[:A],[:B]))

# Domains and codomains.
@test dom(plus(A)) == dom(mmerge(A)) && codom(plus(A)) == codom(mmerge(A))
@test dom(coplus(A)) == dom(mcopy(A)) && codom(coplus(A)) == codom(mcopy(A))
@test dom(zero(A)) == dom(create(A)) && codom(zero(A)) == codom(create(A))
@test dom(cozero(A)) == dom(delete(A)) && codom(cozero(A)) == codom(delete(A))

@test dom(join(R,S)) == A
@test codom(join(R,S)) == B
@test dom(bottom(A,B)) == A
@test codom(bottom(A,B)) == B

# Distinctness of "multiplicative" and "additive" (co)diagonals.
@test plus(A) != mmerge(A)
@test coplus(A) != mcopy(A)
@test zero(A) != create(A)
@test cozero(A) != delete(A)

# Junctions
###########

# Add and remove junctions
#-------------------------

A, B, C = [ Ports([sym]) for sym in [:A, :B, :C] ]
f = singleton_diagram(Box(:f, [:A], [:B]))
g = singleton_diagram(Box(:g, [:B], [:C]))

# Copies.
d = compose(f, mcopy(B))
junctioned = compose(f, junction_diagram(B,1,2))
@test add_junctions(d) == junctioned
@test rem_junctions(junctioned) == d

d = compose(mcopy(A), otimes(f,f))
junctioned = compose(junction_diagram(A,1,2), otimes(f,f))
@test is_permuted_equal(add_junctions(d), junctioned, [3,1,2])
@test rem_junctions(junctioned) == d

# Merges.
d = compose(mmerge(A), f)
junctioned = compose(junction_diagram(A,2,1), f)
@test is_permuted_equal(add_junctions(d), junctioned, [2,1])
@test rem_junctions(junctioned) == d

d = compose(otimes(f,f), mmerge(B))
junctioned = compose(otimes(f,f), junction_diagram(B,2,1))
@test add_junctions(d) == junctioned
@test rem_junctions(junctioned) == d

# Deletions.
d = compose(f, delete(B))
junctioned = compose(f, junction_diagram(B,1,0))
@test add_junctions(d) == junctioned
@test rem_junctions(junctioned) == d

# Creations.
d = compose(create(A), f)
junctioned = compose(junction_diagram(A,0,1), f)
@test is_permuted_equal(add_junctions(d), junctioned, [2,1])
@test rem_junctions(junctioned) == d

# Copies, merges, deletions, and creations, all at once.
d = compose(create(A), f, mcopy(B), mmerge(B), g, delete(C))
junctioned = compose(junction_diagram(A,0,1), f, junction_diagram(B,1,2),
                     junction_diagram(B,2,1), g, junction_diagram(C,1,0))
actual = add_junctions(d)
# XXX: An isomorphism test would be more convenient.
perm = [ findfirst([b] .== boxes(actual)) for b in boxes(junctioned) ]
@test is_permuted_equal(actual, junctioned, perm)
@test rem_junctions(junctioned) == d

# Simplify junctions
#-------------------

A = Ports([:A])
j = junction_diagram

# Comonoid laws.
@test merge_junctions(j(A,1,2)⋅(j(A,1,2)⊗id(A))) ==
  j(A,1,3)⋅permute(A⊗A⊗A,[3,1,2]) # FIXME: Shouldn't need permutation.
@test merge_junctions(j(A,1,2)⋅(id(A)⊗j(A,1,2))) == j(A,1,3)
@test merge_junctions(j(A,1,2)⋅(j(A,1,0)⊗id(A))) == j(A,1,1)
@test merge_junctions(j(A,1,2)⋅(id(A)⊗j(A,1,0))) == j(A,1,1)

# Caps and cups.
@test merge_junctions(j(A,0,1)⋅j(A,1,2)) == j(A,0,2)
@test merge_junctions(j(A,2,1)⋅j(A,1,0)) == j(A,2,0)

# Zigzag laws.
@test merge_junctions((id(A)⊗j(A,0,2))⋅(j(A,2,0)⊗id(A))) == j(A,1,1)
@test merge_junctions((j(A,0,2)⊗id(A))⋅(id(A)⊗j(A,2,0))) == j(A,1,1)

end


text:
using SweepContractor

L=3; d=2

LTN = LabelledTensorNetwork{Tuple{Int,Int,Int}}()
for i∈1:L, j∈1:L, k∈1:L
    adj=Tuple{Int,Int,Int}[];
    i>1 && push!(adj,(i-1,j,k))
    i<L && push!(adj,(i+1,j,k))
    j>1 && push!(adj,(i,j-1,k))
    j<L && push!(adj,(i,j+1,k))
    k>1 && push!(adj,(i,j,k-1))
    k<L && push!(adj,(i,j,k+1))
    LTN[i,j,k] = Tensor(
        adj,
        randn(d*ones(Int,length(adj))...),
        i+0.01*randn(),
        j+0.01*randn()
    )
end

println("example_3d")
for χ∈d.^(4:9)
    sweep=sweep_contract(LTN, χ, 2*χ)
    println("χ=$χ:\t",ldexp(sweep...))
end

ext: jl


text:
#==========================================================================================
                                Constructors
==========================================================================================#
function DiaSymSemiseparableMatrix(U::AbstractArray, V::AbstractArray, d::AbstractArray)
    if size(U,1) == size(V,1) && size(U,2) == size(V,2) && length(d) == size(U,2)
        return DiaSymSemiseparableMatrix(size(U,2),size(U,1),U,V,d);
    else
        error("Dimension mismatch between the generators U, V and d")
    end
end
function DiaSymSemiseparableMatrix(L::SymSemiseparableMatrix, d::AbstractArray)
    return DiaSymSemiseparableMatrix(L.n, L.p, L.Ut, L.Vt, d)
end
function DiaSymSemiseparableMatrix(L::DiaSymSemiseparableCholesky)
    V, d = dss_create_vd(L.Ut, L.Wt, L.ds);
    return DiaSymSemiseparableMatrix(L.n, L.p, L.Ut, V, d)
end
#==========================================================================================
                    Defining Matrix Properties / Overloading Base
==========================================================================================#
Base.Matrix(K::DiaSymSemiseparableMatrix) = tril(K.Ut'*K.Vt) + triu(K.Vt'*K.Ut,1) + Diagonal(K.d)
Base.size(K::DiaSymSemiseparableMatrix)   = (K.n,K.n)
function getindex(K::DiaSymSemiseparableMatrix{T}, i::Int, j::Int) where T
	i > j && return dot(K.Ut[:,i],K.Vt[:,j])
	j == i && return dot(K.Vt[:,i],K.Ut[:,j]) + K.d[i]
	return dot(K.Vt[:,i],K.Ut[:,j])
end
Base.propertynames(F::DiaSymSemiseparableMatrix, private::Bool=false) =
    (private ? fieldnames(typeof(F)) : ())
#==========================================================================================
                            Overloading LinearAlgebra routines
==========================================================================================#
LinearAlgebra.cholesky(K::DiaSymSemiseparableMatrix) = DiaSymSemiseparableCholesky(K)
LinearAlgebra.logdet(K::DiaSymSemiseparableMatrix) = 2.0*logdet(cholesky(K))
LinearAlgebra.det(K::DiaSymSemiseparableMatrix)    = det(cholesky(K))^2
LinearAlgebra.logdet(L::Adjoint{<:Any,<:DiaSymSemiseparableMatrix}) = logdet(L.parent)
LinearAlgebra.det(K::Adjoint{<:Any,<:DiaSymSemiseparableMatrix})    = det(cholesky(K.parent))^2
#==========================================================================================
                        Defining multiplication and inverse
==========================================================================================#
function LinearAlgebra.mul!(y::AbstractArray, L::DiaSymSemiseparableMatrix, b::AbstractArray)
    dss_mul_mat!(y, L.Ut, L.Vt, L.d, b)
    return y
end
function LinearAlgebra.mul!(y::AbstractArray, L::Adjoint{<:Any,<:DiaSymSemiseparableMatrix}, b::AbstractArray)
    dss_mul_mat!(y, L.parent.Ut, L.parent.Vt, L.parent.d, b)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray, K::DiaSymSemiseparableMatrix, b::AbstractArray)
	L = DiaSymSemiseparableCholesky(K)
	y .= L'\(L\b)
    return y
end
function LinearAlgebra.ldiv!(y::AbstractArray, K::Adjoint{<:Any,<:DiaSymSemiseparableMatrix}, b::AbstractArray)
	L = DiaSymSemiseparableCholesky(K.parent)
	y .= L'\(L\b)
    return y
end
function (Base.:\)(K::DiaSymSemiseparableMatrix, b::AbstractVecOrMat)
    y = similar(b)
	L = DiaSymSemiseparableCholesky(K)
	y .= L'\(L\b)
    return y
end
function (Base.:\)(K::Adjoint{<:Any,<:DiaSymSemiseparableMatrix}, b::AbstractVecOrMat)
    y = similar(b)
	L = DiaSymSemiseparableCholesky(K.parent)
	y .= L'\(L\b)
    return y
end

#===========================================================================================
                Cholesky factoriaztion of:  Higher-order quasiseparable matrices
===========================================================================================#
"""
    dss_mul_mat!(Y, U, V, d, X)

Computes the matrix-matrix product `Y = (tril(U*V') + triu(V*U',1) + diag(d))*X`.
"""
function dss_mul_mat!(Y, U, V, d, X)
    Ubar = U*X
    Vbar = zeros(eltype(X),size(Ubar))
    @inbounds for (u,v,y,x,i) in zip(eachcol(U),eachcol(V),eachrow(Y),eachrow(X),eachindex(d))
        add_inner_minus!(Ubar,u,x)
        add_inner_plus!(Vbar,v,x)
        add_Y_diag!(y,Ubar,Vbar,u,v,x,d[i]) # Y[i,:] = tmpU'*Vbar + tmpV'*Ubar + d[i]*tmpX
    end
	return Y
end

"""
    dss_create_vd(U, W, dbar)

Using `L = tril(U*W',-1) + Diagonal(dbar)`, compute `V` and `d` such that
`LL = tril(UV') + triu(V'*U,1) + Diagonal(d)`.
"""
function dss_create_vd(U, W, dbar)
    p, n = size(U)
    d = zeros(n)
    V = zeros(n,p)
    P = zeros(p,p)
    @inbounds for (u,w,v,i) in zip(eachcol(U),eachcol(W),eachcol(V),eachindex(dbar))
        d[i]  = dbar[i]^2 - dot(u,w)
        v    .= P*u
        add_outer_product!(P,w)
    end
    return V, d
end


text:
mutable struct fCTM <: TopicModel
	K::Int
	M::Int
	V::Int
	N::Vector{Int}
	C::Vector{Int}
	corp::Corpus
	topics::VectorList{Int}
	mu::Vector{Float64}
	sigma::Matrix{Float64}
	eta::Float64
	beta::Matrix{Float64}
	fbeta::Matrix{Float64}
	kappa::Vector{Float64}
	lambda::VectorList{Float64}
	vsq::VectorList{Float64}
	lzeta::Float64
	tau::VectorList{Float64}
	phi::Matrix{Float64}
	invsigma::Matrix{Float64}
	newbeta::Matrix{Float64}
	newkappa::Vector{Float64}
	elbo::Float64
	newelbo::Float64

	function fCTM(corp::Corpus, K::Integer)
		@assert ispositive(K)
		@assert !isempty(corp)
		checkcorp(corp)

		M, V, U = size(corp)
		N = [length(doc) for doc in corp]
		C = [size(doc) for doc in corp]

		topics = [collect(1:V) for _ in 1:K]

		mu = zeros(K)
		sigma = eye(K)
		invsigma = eye(K)
		eta = 0.95
		neweta = 0.0
		beta = rand(Dirichlet(V, 1.0), K)'
		newbeta = zeros(K, V)
		fbeta = beta
		kappa = rand(Dirichlet(V, 1.0))
		newkappa = zeros(V)
		lambda = [zeros(K) for _ in 1:M]
		vsq = [ones(K) for _ in 1:M]
		lzeta = 0.0
		tau = [fill(eta, N[d]) for d in 1:M]
		phi = ones(K, N[1]) / K

		model = new(K, M, V, N, C, copy(corp), topics, mu, sigma, eta, beta, fbeta, kappa, lambda, vsq, lzeta, tau, phi)
		fixmodel!(model, check=false)

		for d in 1:M
			model.phi = ones(K, N[d]) / K
			updateNewELBO!(model, d)
		end
		model.phi = ones(K, N[1]) / K
		updateELBO!(model)
		return model
	end
end

function Elogpeta(model::fCTM, d::Int)
	x = 0.5 * (logdet(model.invsigma) - model.K * log(2pi) - dot(diag(model.invsigma), model.vsq[d]) - dot(model.lambda[d] - model.mu, model.invsigma * (model.lambda[d] - model.mu)))
	return x
end

function Elogpc(model::fCTM, d::Int)
	counts = model.corp[d].counts
	y = dot(model.tau[d], counts)
	x = log(@boink model.eta^y * (1 - model.eta)^(model.C[d] - y))
	return x
end

function Elogpz(model::fCTM, d::Int)
	counts = model.corp[d].counts
	x = dot(model.phi' * model.lambda[d], counts) + model.C[d] * model.lzeta
	return x
end

function Elogpw(model::fCTM, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	x = sum(model.phi .* log.(@boink model.beta[:,terms]) * (model.tau[d] .* counts)) + dot(1 - model.tau[d], log.(@boink model.kappa[terms]))
	return x
end

function Elogqeta(model::fCTM, d::Int)
	x = -entropy(MvNormal(model.lambda[d], diagm(model.vsq[d])))
	return x
end

function Elogqc(model::fCTM, d::Int)
	counts = model.corp[d].counts
	x = -sum([c * entropy(Bernoulli(model.tau[d][n])) for (n, c) in enumerate(counts)])
	return x
end

function Elogqz(model::fCTM, d::Int)
	counts = model.corp[d].counts
	x = -sum([c * entropy(Categorical(model.phi[:,n])) for (n, c) in enumerate(counts)])
	return x
end

function updateELBO!(model::fCTM)
	model.elbo = model.newelbo
	model.newelbo = 0
	return model.elbo
end

function updateNewELBO!(model::fCTM, d::Int)
	model.newelbo += (Elogpeta(model, d)
					+ Elogpc(model, d)
					+ Elogpz(model, d)
					+ Elogpw(model, d)
					- Elogqeta(model, d)
					- Elogqc(model, d)
					- Elogqz(model, d))
end

function updateMu!(model::fCTM)
	model.mu = sum(model.lambda) / model.M
end

function updateSigma!(model::fCTM)
	model.sigma = diagm(sum(model.vsq)) / model.M + Base.covm(hcat(model.lambda...), model.mu, 2, false)
	model.invsigma = inv(model.sigma)
end

function updateEta!(model::fCTM)
	model.eta = sum([dot(model.tau[d], model.corp[d].counts) for d in 1:model.M]) / sum(model.C)
end

function updateBeta!(model::fCTM)
	model.beta = model.newbeta ./ sum(model.newbeta, 2)
	model.newbeta = zeros(model.K, model.V)
end

function updateNewBeta!(model::fCTM, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.newbeta[:,terms] += model.phi .* (model.tau[d] .* counts)'
end

function updateKappa!(model::fCTM)
	model.kappa = model.newkappa / sum(model.newkappa)
	model.newkappa = zeros(model.V)
end

function updateNewKappa!(model::fCTM, d::Int)
	terms, counts = model.corp[d].terms, model.corp[d].counts
	model.newkappa[terms] += (1 - model.tau[d]) .* counts
end

function updateLambda!(model::fCTM, d::Int, niter::Integer, ntol::Real)
	"Newton's method."

	counts = model.corp[d].counts
	for _ in 1:niter
		lambdaGrad = model.invsigma * (model.mu - model.lambda[d]) + model.phi * counts - model.C[d] * exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta)
		lambdaInvHess = -inv(eye(model.K) + model.C[d] * model.sigma * diagm(exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta))) * model.sigma
		model.lambda[d] -= lambdaInvHess * lambdaGrad
		if norm(lambdaGrad) < ntol
			break
		end
	end
end

function updateVsq!(model::fCTM, d::Int, niter::Integer, ntol::Real)
	"Newton's method."

	for _ in 1:niter
		rho = 1.0
		vsqGrad = -0.5 * (diag(model.invsigma) + model.C[d] * exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta) - 1 ./ model.vsq[d])
		vsqInvHess = -1 ./ (0.25 * model.C[d] * exp.(model.lambda[d] + 0.5 * model.vsq[d] - model.lzeta) + 0.5 ./ model.vsq[d].^2)
		p = vsqInvHess .* vsqGrad

		while minimum(model.vsq[d] - rho * p) <= 0
			rho *= 0.5
		end
		model.vsq[d] -= rho * p

		if norm(vsqGrad) < ntol
			break
		end
	end
	@bumper model.vsq[d]
end

function updateLzeta!(model::fCTM, d::Int)
	model.lzeta = logsumexp(model.lambda[d] + 0.5 * model.vsq[d])
end

function updateTau!(model::fCTM, d::Int)
	terms = model.corp[d].terms
	model.tau[d] = model.eta ./ (@boink model.eta + (1 - model.eta) * (model.kappa[terms] .* vec(prod(model.beta[:,terms].^-model.phi, 1))))
end

function updatePhi!(model::fCTM, d::Int)
	terms = model.corp[d].terms
	model.phi = addlogistic(model.tau[d]' .* log.(model.beta[:,terms]) .+ model.lambda[d], 1)
end

function train!(model::fCTM; iter::Integer=150, tol::Real=1.0, niter=1000, ntol::Real=1/model.K^2, viter::Integer=10, vtol::Real=1/model.K^2, chkelbo::Integer=1)
	@assert all(.!isnegative.([tol, ntol, vtol]))
	@assert all(ispositive.([iter, niter, viter, chkelbo]))

	for k in 1:iter
		chk = (k % chkelbo == 0)
		for d in 1:model.M
			for _ in 1:viter
				oldlambda = model.lambda[d]
				updatePhi!(model, d)
				updateLzeta!(model, d)
				updateLambda!(model, d, niter, ntol)
				updateVsq!(model, d, niter, ntol)
				if norm(oldlambda - model.lambda[d]) < vtol
					break
				end
			end
			updateTau!(model, d)
			updateNewBeta!(model, d)
			updateNewKappa!(model, d)
			chk && updateNewELBO!(model, d)
		end
		updateMu!(model)
		updateSigma!(model)
		updateEta!(model)
		updateBeta!(model)
		updateKappa!(model)
		if checkELBO!(model, k, chk, tol)
			break
		end
	end
	updatePhi!(model, 1)
	updateLzeta!(model, 1)
	@bumper model.fbeta = model.beta .* (model.kappa' .<= 0)
	model.fbeta ./= sum(model.fbeta, 2)
	model.topics = [reverse(sortperm(vec(model.fbeta[i,:]))) for i in 1:model.K]
	nothing
end


text:
module MLPlotsTests

using VisualRegressionTests
using FactCheck

using LearnBase
using OnlineAI
import ROCAnalysis
import ValueHistories

using MLPlots, Plots

# don't let pyplot use a gui... it'll crash
# note: Agg will set gui -> :none in PyPlot
ENV["MPLBACKEND"] = "Agg"
try
    @eval import PyPlot
    info("Matplotlib version: $(PyPlot.matplotlib[:__version__])")
end
pyplot(size=(400,300), reuse=true)

refdir = Pkg.dir("MLPlots", "test", "refimg")

# run a visual regression test comparing the output to the saved reference png
function dotest(testname, func)
    srand(1)
    reffn = joinpath(refdir, "$testname.png")
    vtest = VisualTest(func, reffn)
    result = test_images(vtest)
    @fact success(result) --> true
end

# builds a testable function which saves a png to the location `fn`
# use: VisualTest(@plotfunc(plot(rand(10))), "/tmp/tmp.png")
macro plottest(testname, expr)
    esc(quote
        dotest(string($testname), fn -> begin
            $expr
            png(fn)
        end)
    end)
end

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

facts("LearnBase") do
    with(size=(300,300), guidefont=font(9), tickfont=font(7)) do
        @plottest "learnbase1" begin
            plot(LossFunctions.HingeLoss())
        end
        @plottest "learnbase2" begin
            plot([ZeroOneLoss(), L1HingeLoss(), L2HingeLoss(), LogitMarginLoss()])
        end
        @plottest "learnbase3" begin
            subplot([ZeroOneLoss(), L1HingeLoss(), L2HingeLoss(), LogitMarginLoss()], size=(400,400))
        end
        @plottest "learnbase4" begin
            plot(LossFunctions.L2DistLoss())
        end
        @plottest "learnbase5" begin
            plot([L2DistLoss(), L1DistLoss(), EpsilonInsLoss(.4), LogitDistLoss()])
        end
        @plottest "learnbase6" begin
            subplot([L2DistLoss(), L1DistLoss(), EpsilonInsLoss(.4), LogitDistLoss()], size=(400,400))
        end
    end
end

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

facts("CorrPlot") do
    @plottest "corrplot" begin
        M = randn(1000, 4)
        M[:,2] += 0.8M[:,1]
        M[:,3] -= 0.7M[:,1]
        corrplot(M, size=(500,500))
    end
end

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

facts("OnlineAI") do

    # build and fit a neural net, then show it
    @plottest "onlineai1" begin
        net = buildClassificationNet(3,1,[15,10,5])
        for i in 1:1000
            x = rand(3)
            fit!(net, x, Float64[sum(x)])
        end
        plot(net)
    end

    @plottest "onlineai2" begin
        n = 20
        spikes = SpikeTrains(n, title = "Spike Trains", color = :darkblue)
        for t=1:100, i=1:n
            if rand() < 0.1
                push!(spikes, i, t)
            end
        end
        spikes.plt
    end
end

# ---------------------------------------------------------\

facts("ROCAnalysis") do
    @plottest "rocanalysis" begin
        plot(ROCAnalysis.roc(2+2randn(1000), -2+2randn(100000)))
    end
end

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

facts("ValueHistories") do
    @plottest "valuehistories1" begin
        history = ValueHistories.DynMultivalueHistory()
        for i=1:100
            x = 0.1i
            push!(history, :a, x, sin(x))
            push!(history, :wrongtype, x, "$(sin(x))")
            if i % 10 == 0
                push!(history, :b, x, cos(x))
            end
        end
        plot(history)
    end

    @plottest "valuehistories2" begin
        history = ValueHistories.QueueUnivalueHistory(Int)
        for i = 1:100
            push!(history, i, 2i)
        end
        plot(history)
    end

    @plottest "valuehistories3" begin
        history = ValueHistories.VectorUnivalueHistory(Int)
        for i = 1:100
            push!(history, i, 100-i)
        end
        plot(history)
    end
end

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

FactCheck.exitstatus()
end # module


text:
module FlatArraysBenchmarking

using BenchmarkTools, FlatArrays
using Base:
    has_offset_axes

function print_left_justified(str::AbstractString, len::Integer, c::Char=' ')
    print(str)
    (n = len - length(str)) > 0 && print(c^n)
end

function print_right_justified(str::AbstractString, len::Integer, c::Char=' ')
    (n = len - length(str)) > 0 && print(c^n)
    print(str)
end

# Use eachindex()
function add0!(dest::AbstractArray{T,N},
               A::AbstractArray{T,N},
               B::AbstractArray{T,N}) where {T,N}
    # Checks must be skipped because they take time.
    #@assert axes(dest) == axes(A) == axes(B)
    @inbounds @simd for i in eachindex(dest, A, B)
        dest[i] = A[i] + B[i]
    end
    return dest
end

# Use linear indexing.
function add1!(dest::AbstractArray{T,N},
               A::AbstractArray{T,N},
               B::AbstractArray{T,N}) where {T,N}
    # Checks must be skipped because they take time.
    #@assert !Base.has_offset_axes(dest, A, B)
    #@assert (n = length(dest)) == length(A) == length(B)
    n = length(dest)
    @inbounds @simd for i in Base.OneTo(n)
        dest[i] = A[i] + B[i]
    end
    return dest
end

# Use 2D indexing.
function add2!(dest::AbstractArray{T,2},
               A::AbstractArray{T,2},
               B::AbstractArray{T,2}) where {T}
    # Checks must be skipped because they take time.
    #@assert (m = size(dest,1)) == size(A,1) == size(B,1)
    #@assert (n = size(dest,2)) == size(A,2) == size(B,2)
    m = size(dest,1)
    n = size(dest,2)
    @inbounds for j in 1:n
        @simd for i in 1:m
            dest[i,j] = A[i,j] + B[i,j]
        end
    end
    return dest
end

# Use 2D indexing + faking.
function add2p!(dest::AbstractArray{T,2},
                A::AbstractArray{T,2},
                B::AbstractArray{T,2}) where {T}
    # Checks must be skipped because they take time.
    #@assert (m = size(dest,1)) == size(A,1) == size(B,1)
    #@assert (n = size(dest,2)) == size(A,2) == size(B,2)
    m = size(dest,1)
    n = size(dest,2)
    @inbounds for j in 1:n
        off = m*(j - 1)
        @simd for i in 1:m
            k = i + off
            dest[k] = A[k] + B[k]
        end
    end
    return dest
end

# Use CartesianIndex()
function addI!(dest::AbstractArray{T,N},
               A::AbstractArray{T,N},
               B::AbstractArray{T,N}) where {T,N}
    # Checks must be skipped because they take time.
    #@assert axes(dest) == axes(A) == axes(B)
    @inbounds @simd for i in CartesianIndices(dest)
        dest[i] = A[i] + B[i]
    end
    return dest
end

function runtests()
    dims = (3,4,5,6,7)
    A = randn(dims)
    B = view(A, :, :, :, 1:4, 3) # flat view
    C = view(A, :, :, :, 2:5, :) # non-flat view
    @assert isflat(A)
    @assert isflat(B)
    @assert !isflat(C)

    nelem = length(A)
    nrows = prod(dims[1:3])
    ncols = prod(dims[4:end])
    V = flatten(A, nelem)
    M = flatten(A, nrows, ncols)

    l = 30
    c = '.'
    if false
    println("In the following tests, A is a regular array, B and C are views of A,")
    println("B is flat, C is not flat, V is a flat vector from A and M is a flat")
    println("matrix from A.")
    println()
    print_left_justified("ndims(A) ",               l, c); @btime ndims($A);
    print_left_justified("ndims(B) ",               l, c); @btime ndims($B);
    print_left_justified("ndims(C) ",               l, c); @btime ndims($C);
    print_left_justified("ndims(V) ",               l, c); @btime ndims($V);
    print_left_justified("ndims(M) ",               l, c); @btime ndims($M);
    print_left_justified("firstindex(A) ",          l, c); @btime firstindex($A);
    print_left_justified("firstindex(B) ",          l, c); @btime firstindex($B);
    print_left_justified("firstindex(C) ",          l, c); @btime firstindex($C);
    print_left_justified("firstindex(V) ",          l, c); @btime firstindex($V);
    print_left_justified("firstindex(M) ",          l, c); @btime firstindex($M);
    print_left_justified("lastindex(A) ",           l, c); @btime lastindex($A);
    print_left_justified("lastindex(B) ",           l, c); @btime lastindex($B);
    print_left_justified("lastindex(C) ",           l, c); @btime lastindex($C);
    print_left_justified("lastindex(V) ",           l, c); @btime lastindex($V);
    print_left_justified("lastindex(M) ",           l, c); @btime lastindex($M);
    print_left_justified("iterate(A) ",             l, c); @btime iterate($A);
    print_left_justified("iterate(B) ",             l, c); @btime iterate($B);
    print_left_justified("iterate(C) ",             l, c); @btime iterate($C);
    print_left_justified("iterate(V) ",             l, c); @btime iterate($V);
    print_left_justified("iterate(M) ",             l, c); @btime iterate($M);
    print_left_justified("IndexStyle(A) ",          l, c); @btime IndexStyle($A);
    print_left_justified("IndexStyle(B) ",          l, c); @btime IndexStyle($B);
    print_left_justified("IndexStyle(C) ",          l, c); @btime IndexStyle($C);
    print_left_justified("IndexStyle(V) ",          l, c); @btime IndexStyle($V);
    print_left_justified("IndexStyle(M) ",          l, c); @btime IndexStyle($M);
    print_left_justified("has_offset_axes(A) ",     l, c); @btime has_offset_axes($A);
    print_left_justified("has_offset_axes(B) ",     l, c); @btime has_offset_axes($B);
    print_left_justified("has_offset_axes(C) ",     l, c); @btime has_offset_axes($C);
    print_left_justified("has_offset_axes(V) ",     l, c); @btime has_offset_axes($V);
    print_left_justified("has_offset_axes(M) ",     l, c); @btime has_offset_axes($M);
    print_left_justified("reshape(A,nelem) ",       l, c); @btime reshape($A,$nelem);
    print_left_justified("flatten(A,nelem) ",       l, c); @btime flatten($A,$nelem);
    print_left_justified("reshape(A,nrows,ncols) ", l, c); @btime reshape($A,$nrows,$ncols);
    print_left_justified("flatten(A,nrows,ncols) ", l, c); @btime flatten($A,$nrows,$ncols);
    print_left_justified("isflat(A) ",              l, c); @btime isflat($A);
    print_left_justified("isflat(B) ",              l, c); @btime isflat($B);
    print_left_justified("isflat(C) ",              l, c); @btime isflat($C);
    print_left_justified("isflat(V) ",              l, c); @btime isflat($V);
    print_left_justified("isflat(M) ",              l, c); @btime isflat($M);
    println()
    println("In the following tests, A is a regular array, R is a reshaped version of A,")
    println("M is a flat version of A and S is a random array of same size as R and M.")
    println()
    R = reshape(A, nrows, ncols)
    S = randn(size(R))
    print_left_justified("R'*R ",                   l, c); @btime $R'*$R;
    print_left_justified("S'*R ",                   l, c); @btime $S'*$R;
    print_left_justified("R'*S ",                   l, c); @btime $R'*$S;
    print_left_justified("M'*M ",                   l, c); @btime $M'*$M;
    print_left_justified("M'*R ",                   l, c); @btime $M'*$R;
    print_left_justified("R'*M ",                   l, c); @btime $R'*$M;
    println()
    end
    println("In the following tests:")
    println(" - `add0!` uses `eachindex`;")
    println(" - `add1!` uses 1D linear index;")
    println(" - `add2!` uses 2D linear indices;")
    println(" - `add2p!` uses 2D linear indices with a common stride for all arrays;")
    println(" - `addI!` uses CartesianIndex;")
    println()
    println("Tests with regular Array's of size $(dims):")
    A1 = randn(dims)
    A2 = randn(dims)
    A0 = similar(A1)
    print_left_justified("  add0!(A0, A1, A2) ", l, c); @btime add0!($A0, $A1, $A2);
    print_left_justified("  add1!(A0, A1, A2) ", l, c); @btime add1!($A0, $A1, $A2);
    print_left_justified("  addI!(A0, A1, A2) ", l, c); @btime addI!($A0, $A1, $A2);
    println()
    println("Tests with regular Vector's of length $(nelem):")
    B1 = randn(nelem) # true vector
    B2 = randn(nelem) # true vector
    B0 = similar(B1)  # true vector
    print_left_justified("  add0!(B0, B1, B2) ", l, c); @btime add0!($B0, $B1, $B2);
    print_left_justified("  add1!(B0, B1, B2) ", l, c); @btime add1!($B0, $B1, $B2);
    print_left_justified("  addI!(B0, B1, B2) ", l, c); @btime addI!($B0, $B1, $B2);
    println()
    println("Tests with FlatVector's of length $(nelem):")
    V1 = flatten(A1, nelem)
    V2 = flatten(A2, nelem)
    V0 = flatten(A0, nelem)
    print_left_justified("  add0!(V0, V1, V2) ", l, c); @btime add0!($V0, $V1, $V2);
    print_left_justified("  add1!(V0, V1, V2) ", l, c); @btime add1!($V0, $V1, $V2);
    print_left_justified("  addI!(V0, V1, V2) ", l, c); @btime addI!($V0, $V1, $V2);
    println()
    println("Tests with reshaped Vector's of length $(nelem):")
    R1 = reshape(A1, nelem)
    R2 = reshape(A2, nelem)
    R0 = reshape(A0, nelem)
    print_left_justified("  add0!(R0, R1, R2) ", l, c); @btime add0!($R0, $R1, $R2);
    print_left_justified("  add1!(R0, R1, R2) ", l, c); @btime add1!($R0, $R1, $R2);
    print_left_justified("  addI!(R0, R1, R2) ", l, c); @btime addI!($R0, $R1, $R2);
    println()
    println("Tests with regular Matrix's of size $((nrows, ncols)):")
    C1 = randn(nrows,ncols) # true matrix
    C2 = randn(nrows,ncols) # true matrix
    C0 = similar(C1)        # true matrix
    print_left_justified("  add0!(C0, C1, C2) ", l, c); @btime add0!($C0, $C1, $C2);
    print_left_justified("  add1!(C0, C1, C2) ", l, c); @btime add1!($C0, $C1, $C2);
    print_left_justified("  add2!(C0, C1, C2) ", l, c); @btime add2!($C0, $C1, $C2);
    print_left_justified("  add2p!(C0, C1, C2) ",l, c); @btime add2p!($C0, $C1, $C2);
    print_left_justified("  addI!(C0, C1, C2) ", l, c); @btime addI!($C0, $C1, $C2);
    println()
    println("Tests with FlatMatrix's of size $((nrows, ncols)):")
    M1 = flatten(A1, nrows, ncols)
    M2 = flatten(A2, nrows, ncols)
    M0 = flatten(A0, nrows, ncols)
    print_left_justified("  add0!(M0, M1, M2) ", l, c); @btime add0!($M0, $M1, $M2);
    print_left_justified("  add1!(M0, M1, M2) ", l, c); @btime add1!($M0, $M1, $M2);
    print_left_justified("  add2!(M0, M1, M2) ", l, c); @btime add2!($M0, $M1, $M2);
    print_left_justified("  add2p!(M0, M1, M2) ",l, c); @btime add2p!($M0, $M1, $M2);
    print_left_justified("  addI!(M0, M1, M2) ", l, c); @btime addI!($M0, $M1, $M2);
    println()
    println("Tests with reshaped Matrix's of size $((nrows, ncols)):")
    S1 = reshape(A1, nrows, ncols)
    S2 = reshape(A2, nrows, ncols)
    S0 = reshape(A0, nrows, ncols)
    print_left_justified("  add0!(S0, S1, S2) ", l, c); @btime add0!($S0, $S1, $S2);
    print_left_justified("  add1!(S0, S1, S2) ", l, c); @btime add1!($S0, $S1, $S2);
    print_left_justified("  add2!(S0, S1, S2) ", l, c); @btime add2!($S0, $S1, $S2);
    print_left_justified("  add2p!(S0, S1, S2) ",l, c); @btime add2p!($S0, $S1, $S2);
    print_left_justified("  addI!(S0, S1, S2) ", l, c); @btime addI!($S0, $S1, $S2);

end
runtests();
end # module

max_stars_repo_name: emmt/FlatArrays

max_issues_repo_name: emmt/FlatArrays

max_forks_repo_name: emmt/FlatArrays


text:
# independent constructor
function Source(fullpath::Union{AbstractString,IO};

              delim=COMMA,
              quotechar=QUOTE,
              escapechar=ESCAPE,
              missingstring::AbstractString="",
              null::Union{AbstractString,Nothing}=nothing,

              header::Union{Integer, UnitRange{Int}, Vector}=1, # header can be a row number, range of rows, or actual string vector
              datarow::Int=-1, # by default, data starts immediately after header or start of file
              types=Type[],
              allowmissing::Symbol=:all,
              nullable::Union{Bool,Missing,Nothing}=nothing,
              dateformat=nothing,
              decimal=PERIOD,
              truestring="true",
              falsestring="false",
              categorical::Bool=true,
              weakrefstrings::Union{Bool,Nothing}=nothing,
              strings::Symbol=:intern,

              footerskip::Int=0,
              rows_for_type_detect::Int=100,
              rows::Int=0,
              use_mmap::Bool=true)
    # make sure character args are UInt8
    isascii(delim) || throw(ArgumentError("non-ASCII characters not supported for delim argument: $delim"))
    isascii(quotechar) || throw(ArgumentError("non-ASCII characters not supported for quotechar argument: $quotechar"))
    isascii(escapechar) || throw(ArgumentError("non-ASCII characters not supported for escapechar argument: $escapechar"))
    return CSV.Source(fullpath=fullpath,
                        options=CSV.Options(delim=typeof(delim) <: String ? UInt8(first(delim)) : (delim % UInt8),
                                            quotechar=typeof(quotechar) <: String ? UInt8(first(quotechar)) : (quotechar % UInt8),
                                            escapechar=typeof(escapechar) <: String ? UInt8(first(escapechar)) : (escapechar % UInt8),
                                            missingstring=missingstring, null=null, dateformat=dateformat, decimal=decimal,
                                            truestring=truestring, falsestring=falsestring, internstrings=(strings === :intern)),
                        header=header, datarow=datarow, types=types, allowmissing=allowmissing, nullable=nullable, categorical=categorical,
                        weakrefstrings=weakrefstrings, strings=strings, footerskip=footerskip,
                        rows_for_type_detect=rows_for_type_detect, rows=rows, use_mmap=use_mmap)
end

function Source(;fullpath::Union{AbstractString,IO}="",
                options::CSV.Options{D}=CSV.Options(),

                header::Union{Integer,UnitRange{Int},Vector}=1, # header can be a row number, range of rows, or actual string vector
                datarow::Int=-1, # by default, data starts immediately after header or start of file
                types=Type[],
                allowmissing::Symbol=:all,
                nullable::Union{Bool,Missing,Nothing}=nothing,
                categorical::Bool=true,
                weakrefstrings::Union{Bool,Nothing}=nothing,
                strings::Symbol=:intern,

                footerskip::Int=0,
                rows_for_type_detect::Int=100,
                rows::Int=0,
                use_mmap::Bool=true) where {D}
    # argument checks
    isa(fullpath, AbstractString) && (isfile(fullpath) || throw(ArgumentError("\"$fullpath\" is not a valid file")))
    header = (isa(header, Integer) && header == 1 && datarow == 1) ? -1 : header
    isa(header, Integer) && datarow != -1 && (datarow > header || throw(ArgumentError("data row ($datarow) must come after header row ($header)")))
    if options.null !== nothing
        resize!(options.missingstring, length(options.null))
        copyto!(options.missingstring, options.null)
        Base.depwarn("null option is deprecated, use missingstring instead", :Source)
    end
    if nullable !== nothing
        allowmissing = ismissing(nullable) ? :auto :
                       nullable            ? :all  : :none
        Base.depwarn("nullable=$nullable argument is deprecated, use allowmissing=$(repr(allowmissing)) instead", :Source)
    end
    if weakrefstrings !== nothing
        strings = weakrefstrings ? :weakref : :raw
        Base.depwarn("weakrefstrings argument is deprecated, use strings=$(repr(strings)) instead", :Source)
    end

    # open the file for property detection
    if isa(fullpath, IOBuffer)
        source = fullpath
        fs = bytesavailable(fullpath)
        fullpath = "<IOBuffer>"
    elseif isa(fullpath, IO)
        source = IOBuffer(Base.read(fullpath))
        fs = bytesavailable(source)
        fullpath = isdefined(fullpath, :name) ? fullpath.name : "__IO__"
    else
        source = open(fullpath, "r") do f
            IOBuffer(use_mmap ? Mmap.mmap(f) : Base.read(f))
        end
        fs = filesize(fullpath)
    end
    options.datarow != -1 && (datarow = options.datarow)
    options.rows != 0 && (rows = options.rows)
    options.header != 1 && (header = options.header)
    !isempty(options.types) && (types = options.types)
    startpos = position(source)
    rows = rows == 0 ? CSV.countlines(source, options.quotechar, options.escapechar) : rows
    seek(source, startpos)
    # BOM character detection
    if fs > 0 && peekbyte(source) == 0xef
        readbyte(source)
        readbyte(source) == 0xbb || seek(source, startpos)
        readbyte(source) == 0xbf || seek(source, startpos)
    end
    datarow = datarow == -1 ? (isa(header, Vector) ? 0 : last(header)) + 1 : datarow # by default, data starts on line after header
    rows = fs == 0 ? -1 : max(-1, rows - datarow + 1 - footerskip) # rows now equals the actual number of rows in the dataset

    # figure out # of columns and header, either an Integer, AbstractRange, or Vector{String}
    # also ensure that `f` is positioned at the start of data
    row_vals = RawField[]
    if isa(header, Integer)
        # default header = 1
        if header <= 0
            CSV.skipto!(source,1,datarow,options.quotechar,options.escapechar)
            datapos = position(source)
            CSV.readsplitline!(row_vals, source,options.delim,options.quotechar,options.escapechar)
            seek(source, datapos)
            columnnames = ["Column$i" for i = eachindex(row_vals)]
        else
            CSV.skipto!(source,1,header,options.quotechar,options.escapechar)
            columnnames = [strip(x.value) for x in CSV.readsplitline!(row_vals, source,options.delim,options.quotechar,options.escapechar)]
            datarow != header+1 && CSV.skipto!(source,header+1,datarow,options.quotechar,options.escapechar)
            datapos = position(source)
        end
    elseif isa(header, AbstractRange)
        CSV.skipto!(source,1,first(header),options.quotechar,options.escapechar)
        columnnames = [x.value for x in readsplitline!(row_vals,source,options.delim,options.quotechar,options.escapechar)]
        for row = first(header):(last(header)-1)
            for (i,c) in enumerate([x.value for x in readsplitline!(row_vals,source,options.delim,options.quotechar,options.escapechar)])
                columnnames[i] *= "_" * c
            end
        end
        datarow != last(header)+1 && CSV.skipto!(source,last(header)+1,datarow,options.quotechar,options.escapechar)
        datapos = position(source)
    elseif fs == 0
        datapos = position(source)
        columnnames = header
        cols = length(columnnames)
    else
        CSV.skipto!(source,1,datarow,options.quotechar,options.escapechar)
        datapos = position(source)
        readsplitline!(row_vals,source,options.delim,options.quotechar,options.escapechar)
        seek(source,datapos)
        if isempty(header)
            columnnames = ["Column$i" for i in eachindex(row_vals)]
        else
            length(header) == length(row_vals) || throw(ArgumentError("The length of provided header ($(length(header))) doesn't match the number of columns at row $datarow ($(length(row_vals)))"))
            columnnames = header
        end
    end

    # Detect column types
    cols = length(columnnames)
    if isa(types, Vector) && length(types) == cols
        # types might be a Vector{DataType}, which will be a problem if Unions are needed
        columntypes = convert(Vector{Type}, types)
    elseif isa(types, Dict) || isempty(types)
        columntypes = fill!(Vector{Type}(undef, cols), Any)
        levels = [Dict{WeakRefString{UInt8}, Int}() for _ = 1:cols]
        lineschecked = 0
        while !eof(source) && lineschecked < min(rows < 0 ? rows_for_type_detect : rows, rows_for_type_detect)
            lineschecked += 1
            # println("type detecting on row = $lineschecked...")
            for i = 1:cols
                # print("\tdetecting col = $i...")
                typ = CSV.detecttype(source, options, columntypes[i], levels[i])::Type
                # print(typ)
                columntypes[i] = CSV.promote_type2(columntypes[i], typ)
                # println("...promoting to: ", columntypes[i])
            end
        end
        if options.dateformat === nothing && any(x->Missings.T(x) <: Dates.TimeType, columntypes)
            # FIXME: doesn't this override any custom options?!
            # auto-detected TimeType
            options = Options(delim=options.delim, quotechar=options.quotechar, escapechar=options.escapechar,
                              missingstring=options.missingstring, dateformat=Dates.ISODateTimeFormat, decimal=options.decimal,
                              datarow=options.datarow, rows=options.rows, header=options.header, types=options.types)
        end
        if categorical
            for i = 1:cols
                T = columntypes[i]
                if length(levels[i]) / sum(values(levels[i])) < .67 && T !== Missing && Missings.T(T) <: WeakRefString
                    columntypes[i] = substitute(T, CategoricalArrays.catvaluetype(Missings.T(T), UInt32))
                end
            end
        end
    else
        throw(ArgumentError("$cols number of columns detected; `types` argument has $(length(types)) entries"))
    end
    if isa(types, Dict)
        if isa(types, Dict{String})
            @static if VERSION >= v"0.7.0-DEV.3627"
                colinds = indexin(keys(types), columnnames)
            else
                colinds = indexin(collect(keys(types)), columnnames)
            end
        else
            colinds = keys(types)
        end
        for (col, typ) in zip(colinds, values(types))
            columntypes[col] = typ
        end
        autocols = setdiff(1:cols, colinds)
    elseif isempty(types)
        autocols = collect(1:cols)
    else
        autocols = Int[]
    end
    if strings !== :weakref
        columntypes = Type[(T !== Missing && Missings.T(T) <: WeakRefString) ? substitute(T, String) : T for T in columntypes]
    end
    if allowmissing != :auto
        if allowmissing == :all # allow missing values in all automatically detected columns
            for i = autocols
                T = columntypes[i]
                columntypes[i] = Union{Missings.T(T), Missing}
            end
        elseif allowmissing == :none # disallow missing values in all automatically detected columns
            for i = autocols
                T = columntypes[i]
                columntypes[i] = Missings.T(T)
            end
        else
            throw(ArgumentError("allowmissing must be either :all, :none or :auto"))
        end
    end
    seek(source, datapos)
    sch = Data.Schema(columntypes, columnnames, ifelse(rows < 0, missing, rows))
    return Source(sch, options, source, String(fullpath), datapos)
end

# construct a new Source from a Sink
Source(s::CSV.Sink) = CSV.Source(fullpath=s.fullpath, options=s.options)

# Data.Source interface
"reset a `CSV.Source` to its beginning to be ready to parse data from again"
Data.reset!(s::CSV.Source) = (seek(s.io, s.datapos); return nothing)
Data.schema(source::CSV.Source) = source.schema
Data.accesspattern(::Type{<:CSV.Source}) = Data.Sequential
@inline Data.isdone(io::CSV.Source, row, col, rows, cols) = eof(io.io) || (!ismissing(rows) && row > rows)
@inline Data.isdone(io::Source, row, col) = Data.isdone(io, row, col, size(io.schema)...)
Data.streamtype(::Type{<:CSV.Source}, ::Type{Data.Field}) = true
@inline Data.streamfrom(source::CSV.Source, ::Type{Data.Field}, ::Type{T}, row, col::Int) where {T} = CSV.parsefield(source.io, T, source.options, row, col)
Data.reference(source::CSV.Source) = source.io.data

"""
`CSV.read(fullpath::Union{AbstractString,IO}, sink::Type{T}=DataFrame, args...; kwargs...)` => `typeof(sink)`

`CSV.read(fullpath::Union{AbstractString,IO}, sink::Data.Sink; kwargs...)` => `Data.Sink`

parses a delimited file into a Julia structure (a DataFrame by default, but any valid `Data.Sink` may be requested).

Minimal error-reporting happens w/ `CSV.read` for performance reasons; for problematic csv files, try [`CSV.validate`](@ref) which takes exact same arguments as `CSV.read` and provides much more information for why reading the file failed.

Positional arguments:

* `fullpath`; can be a file name (string) or other `IO` instance
* `sink::Type{T}`; `DataFrame` by default, but may also be other `Data.Sink` types that support streaming via `Data.Field` interface; note that the method argument can be the *type* of `Data.Sink`, plus any required arguments the sink may need (`args...`).
                    or an already constructed `sink` may be passed (2nd method above)

Keyword Arguments:

* `delim::Union{Char,UInt8}`: how fields in the file are delimited; default `','`
* `quotechar::Union{Char,UInt8}`: the character that indicates a quoted field that may contain the `delim` or newlines; default `'"'`
* `escapechar::Union{Char,UInt8}`: the character that escapes a `quotechar` in a quoted field; default `'\\'`
* `missingstring::String`: indicates how missing values are represented in the dataset; default `""`
* `dateformat::Union{AbstractString,Dates.DateFormat}`: how dates/datetimes are represented in the dataset; default `Base.Dates.ISODateTimeFormat`
* `decimal::Union{Char,UInt8}`: character to recognize as the decimal point in a float number, e.g. `3.14` or `3,14`; default `'.'`
* `truestring`: string to represent `true::Bool` values in a csv file; default `"true"`. Note that `truestring` and `falsestring` cannot start with the same character.
* `falsestring`: string to represent `false::Bool` values in a csv file; default `"false"`
* `header`: column names can be provided manually as a complete Vector{String}, or as an Int/AbstractRange which indicates the row/rows that contain the column names
* `datarow::Int`: specifies the row on which the actual data starts in the file; by default, the data is expected on the next row after the header row(s); for a file without column names (header), specify `datarow=1`
* `types`: column types can be provided manually as a complete Vector{Type}, or in a Dict to reference individual columns by name or number
* `allowmissing::Symbol=:all`: indicates whether columns should allow for missing values or not, that is whether their element type should be `Union{T,Missing}`; by default, all columns are allowed to contain missing values. If set to `:none`, no column can contain missing values, and if set to `:auto`, only colums which contain missing values in the first `rows_for_type_detect` rows are allowed to contain missing values. Column types specified via `types` are not affected by this argument.
* `footerskip::Int`: indicates the number of rows to skip at the end of the file
* `rows_for_type_detect::Int=100`: indicates how many rows should be read to infer the types of columns
* `rows::Int`: indicates the total number of rows to read from the file; by default the file is pre-parsed to count the # of rows; `-1` can be passed to skip a full-file scan, but the `Data.Sink` must be set up to account for a potentially unknown # of rows
* `use_mmap::Bool=true`: whether the underlying file will be mmapped or not while parsing; note that on Windows machines, the underlying file will not be "deletable" until Julia GC has run (can be run manually via `gc()`) due to the use of a finalizer when reading the file.
* `append::Bool=false`: if the `sink` argument provided is an existing table, `append=true` will append the source's data to the existing data instead of doing a full replace
* `transforms::Dict{Union{String,Int},Function}`: a Dict of transforms to apply to values as they are parsed. Note that a column can be specified by either number or column name.
* `transpose::Bool=false`: when reading the underlying csv data, rows should be treated as columns and columns as rows, thus the resulting dataset will be the "transpose" of the actual csv data.
* `categorical::Bool=true`: read string column as a `CategoricalArray` ([ref](https://github.com/JuliaData/CategoricalArrays.jl)), as long as the % of unique values seen during type detection is less than 67%. This will dramatically reduce memory use in cases where the number of unique values is small.
* `strings::Symbol`=:intern: indicates how to treat strings. By default strings are interned, i.e. a global pool of already encountered strings is used to avoid allocating a new `String` object for each field. If `strings=:raw`, string interning is disabled, which can be faster when most strings are unique. If `strings=:weakref`, [`WeakRefStrings`](https://github.com/quinnj/WeakRefStrings.jl) package is used used to speed up file parsing by avoiding copies; can only be used for the `Sink` objects that support `WeakRefStringArray` columns (note that `WeakRefStringArray` still returns regular `String` elements).

Example usage:
```
julia> dt = CSV.read("bids.csv")
7656334×9 DataFrames.DataFrame
│ Row     │ bid_id  │ bidder_id                               │ auction │ merchandise      │ device      │
├─────────┼─────────┼─────────────────────────────────────────┼─────────┼──────────────────┼─────────────┤
│ 1       │ 0       │ "8dac2b259fd1c6d1120e519fb1ac14fbqvax8" │ "ewmzr" │ "jewelry"        │ "phone0"    │
│ 2       │ 1       │ "668d393e858e8126275433046bbd35c6tywop" │ "aeqok" │ "furniture"      │ "phone1"    │
│ 3       │ 2       │ "aa5f360084278b35d746fa6af3a7a1a5ra3xe" │ "wa00e" │ "home goods"     │ "phone2"    │
...
```

Other example invocations may include:
```julia
# read in a tab-delimited file
CSV.read(file; delim='\t')

# read in a comma-delimited file with missing values represented as '\\N', such as a MySQL export
CSV.read(file; missingstring="\\N")

# read a csv file that happens to have column names in the first column, and grouped data in rows instead of columns
CSV.read(file; transpose=true)

# manually provided column names; must match # of columns of data in file
# this assumes there is no header row in the file itself, so data parsing will start at the very beginning of the file
CSV.read(file; header=["col1", "col2", "col3"])

# manually provided column names, even though the file itself has column names on the first row
# `datarow` is specified to ensure data parsing occurs at correct location
CSV.read(file; header=["col1", "col2", "col3"], datarow=2)

# types provided manually; as a vector, must match length of columns in actual data
CSV.read(file; types=[Int, Int, Float64])

# types provided manually; as a Dict, can specify columns by # or column name
CSV.read(file; types=Dict(3=>Float64, 6=>String))
CSV.read(file; types=Dict("col3"=>Float64, "col6"=>String))

# manually provided # of rows; if known beforehand, this will improve parsing speed
# this is also a way to limit the # of rows to be read in a file if only a sample is needed
CSV.read(file; rows=10000)

# for data files, `file` and `file2`, with the same structure, read both into a single DataFrame
# note that `df` is used as a 2nd argument in the 2nd call to `CSV.read` and the keyword argument
# `append=true` is passed
df = CSV.read(file)
df = CSV.read(file2, df; append=true)

# manually construct a `CSV.Source` once, then stream its data to both a DataFrame
# and SQLite table `sqlite_table` in the SQLite database `db`
# note the use of `CSV.reset!` to ensure the `source` can be streamed from again
source = CSV.Source(file)
df1 = CSV.read(source, DataFrame)
CSV.reset!(source)
db = SQLite.DB()
sq1 = CSV.read(source, SQLite.Sink, db, "sqlite_table")
```
"""
function read end

function read(fullpath::Union{AbstractString,IO}, sink::Type=DataFrame, args...; append::Bool=false, transforms::AbstractDict=Dict{Int,Function}(), transpose::Bool=false, kwargs...)
    source = transpose ? TransposedSource(fullpath; kwargs...) : Source(fullpath; kwargs...)
    sink = Data.stream!(source, sink, args...; append=append, transforms=transforms)
    return Data.close!(sink)
end

function read(fullpath::Union{AbstractString,IO}, sink::T; append::Bool=false, transforms::AbstractDict=Dict{Int,Function}(), transpose::Bool=false, kwargs...) where {T}
    source = transpose ? TransposedSource(fullpath; kwargs...) : Source(fullpath; kwargs...)
    sink = Data.stream!(source, sink; append=append, transforms=transforms)
    return Data.close!(sink)
end

read(source::CSV.Source, sink=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink, args...; append=append, transforms=transforms); return Data.close!(sink))
read(source::CSV.Source, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T} = (sink = Data.stream!(source, sink; append=append, transforms=transforms); return Data.close!(sink))
read(source::CSV.TransposedSource, sink=DataFrame, args...; append::Bool=false, transforms::Dict=Dict{Int,Function}()) = (sink = Data.stream!(source, sink, args...; append=append, transforms=transforms); return Data.close!(sink))
read(source::CSV.TransposedSource, sink::T; append::Bool=false, transforms::Dict=Dict{Int,Function}()) where {T} = (sink = Data.stream!(source, sink; append=append, transforms=transforms); return Data.close!(sink))


text:
#
# This file is part of the DiscreteEvents.jl Julia package, MIT license
#
# Paul Bayer, 2020
#
# This is a Julia package for discrete event simulation
#

import Base.invokelatest

function event! end

# Function barrier for functions and expressions as well for arguments and
# keywords of a `fun`. It allows Expr, Symbol and `fun` as arguments for `fun`.
_evaluate(y) = y    # catchall, gives back the argument
function _evaluate(y::T) where {T<:Function}
    if (threadid() == 1) # || (parentmodule(y) != Main)
        y()
    else
        invokelatest(y)
    end
end
_evaluate(y::T) where {T<:Tuple{Vararg{<:Any}}} = _evaluate.(y)
_evaluate(kw::T) where {T<:Iterators.Pairs} = (; zip(keys(kw), map(i->_evaluate(i), values(kw)) )...)
function _evaluate(y::T) where {T<:Union{Symbol,Expr}}  # Symbol,Expr: eval
    try
        ret = Core.eval(Main,y)
        @warn "Evaluating expressions is slow, use functions instead" maxlog=1
        return ret
    catch
        return y
    end
end

# Function barrier for executing `fun`s.
_invoke(f::F, ::Nothing, ::Nothing) where {F<:Function} = f()
_invoke(f::F, arg::T, ::Nothing) where {F<:Function,T<:Tuple{Vararg{<:Any}}} = f(_evaluate.(arg)...)
_invoke(f::F, ::Nothing, kw) where {F<:Function} = f(; _evaluate(kw)...)
_invoke(f::F, arg::T, kw) where {F<:Function,T<:Tuple{Vararg{<:Any}}} = f(_evaluate.(arg)...; _evaluate(kw)...)
_invoke(f::typeof(event!), arg::T, ::Nothing) where T<:Tuple{Vararg{<:Any}} = f(arg...)
_invoke(f::typeof(event!), arg::T, kw) where T<:Tuple{Vararg{<:Any}} = f(arg...; kw...)

"""
    fun(f::Function, args...; kwargs...)

Save a function `f` and its arguments in a closure for later execution.

# Arguments
The arguments `args...` and keyword arguments `kwargs...` to fun are passed to `f` at
execution but may change their values between beeing captured in `fun` and `f`s later
execution. If `f` needs their current values at execution time there are two possibilities:
1. `fun` can take `fun`s, function closures, symbols or expressions at the place of values
    or variable arguments. They are evaluated at event time just before being passed to f.
    There is one exception: if `f` is an `event!`, its arguments are passed on unevaluated.
2.  A mutable type argument (Array, struct ...) is always current. You can
    also change its content from within a function.

If using `Symbol`s or `Expr` in `fun` you get a one time warning. They are
evaluated at global scope in Module `Main` only and therefore cannot
be used by other modules.
"""
@inline function fun(f::F, args::Vararg{Any, N}; kwargs...) where {F<:Function,N}
    args = ifelse(isempty(args), nothing, args)
    kwargs = ifelse(isempty(kwargs), nothing, kwargs)
    () -> _invoke(f, args, kwargs)
end


text:
import ..Curlcurls._bulk_derivative
import ..Curlcurls._domain_wall_derivative
import ..Curlcurls._surface_derivative

function Curl(lattice::Lattice{1},nnm,nnp,indices,interior,surface,domain_wall) where T
	dx = lattice.dx
	N = length(interior)

	∂ₓ_bulk = _bulk_derivative(lattice,N,interior,domain_wall,surface,nnm,nnp)
	∂ₓ_domain_wall = _domain_wall_derivative(lattice,N,interior,domain_wall,nnm,nnp,1,indices)
	∂ₓ_surface = _surface_derivative(lattice,N,interior,surface,nnm,nnp,indices)
	return Curl{1}(∂ₓ_bulk + ∂ₓ_domain_wall + ∂ₓ_surface)
end

(cc::Curl{1})() = kron(sparse([3,2],[2,3],[1,-1],3,3),cc.curl)

function (cc::Curl{1})(ky::Real,kz::Real)
	cc0 = cc()
	cc1 = kron(sparse([2,3,1,1],[1,1,2,3],1im*[kz,-ky,-kz,ky],3,3),sparse(I,size(cc0)...))
	return cc0+cc1
end


text:
# Adaptive Gauss-Kronrod quadrature routines (arbitrary precision),
# originally written and contributed to Julia by Steven G. Johnson, 2013.
#
# This file was formerly a part of Julia. License is MIT: http://julialang.org/license

"""
The `QuadGK` module implements 1d numerical integration by an adaptive Gauss–Kronrod algorithm.

The package provides three functions: `quadgk`, `gauss`, and `kronrod`.

* `quadgk` performs adaptive integration over arbitrary intervals
* `gauss` computes Gaussian quadrature points and weights for integrating over the interval [-1, 1]
* `kronrod` computes Kronrod points, weights, and embedded Gaussian quadrature weights for integrating over [-1, 1].

Typical usage looks like:
```
using QuadGK
integral, err = quadgk(x -> exp(-x^2), 0, 1, rtol=1e-8)
```
which computes the integral of exp(–x²) from x=0 to x=1 to a relative tolerance of 10⁻⁸,
and returns the approximate `integral = 0.746824132812427` and error estimate
`err = 7.887024366937112e-13`.
"""
module QuadGK

export quadgk, quadgk!, gauss, kronrod, alloc_segbuf

using DataStructures, LinearAlgebra
import Base.Order.Reverse

# an in-place integrand function f!(result, x) and
# temporary mutable data (e.g. arrays) of type R for evalrule
struct InplaceIntegrand{F,R,RI}
    f!::F

    # temporary arrays for evalrule
    fg::R
    fk::R
    Ig::R
    Ik::R
    fx::R
    Idiff::RI

    # final result array
    I::RI
end

InplaceIntegrand(f!::F, I::RI, fx::R) where {F,RI,R} =
    InplaceIntegrand{F,R,RI}(f!, similar(fx), similar(fx), similar(fx), similar(fx), fx, similar(I), I)

include("gausskronrod.jl")
include("evalrule.jl")
include("adapt.jl")
include("weightedgauss.jl")

end # module QuadGK


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

module GeoTables

using Meshes

import Tables
import Meshes

import GADM
import Shapefile as SHP
import ArchGDAL as AG
import GeoInterface as GI

include("conversion.jl")
include("geotable.jl")

"""
    load(fname, layer=0)

Load geospatial table from file `fname` and convert the
`geometry` column to Meshes.jl geometries. Optionally,
specify the layer of geometries to read in the file.

Currently supported file types are:

- `*.shp` via Shapefile.jl
- Other formats via ArchGDAL.jl
"""
function load(fname, layer=0)
  if endswith(fname, ".shp")
    table = SHP.Table(fname)
  else # fallback to GDAL
    file  = AG.read(fname)
    data  = AG.getlayer(file, layer)
    table = Tables.columns(data)
  end
  GeoTable(table)
end

"""
    gadm(country, subregions...; children=true)

(Down)load GADM table using `GADM.get` and convert
the `geometry` column to Meshes.jl geometries.

If `children` is `true`, return the table with the
geometries of the subregions, otherwise return a single
geometry for the final subregion in the specification.
"""
function gadm(country, subregions...; children=false)
  data  = GADM.get(country, subregions...; children=children)
  table = children ? data[2] : data
  GeoTable(table)
end

end


text:
using Base.Test
using PowerSeries

# Represent the series 1.0 + 1.0*x - 2.0*x^2
a = series(1.0, 1.0, -2.0)

# Represent the series 1.0 + 0.0*x + 1.0*x^2
b = series(1.0, 0.0, 1.0)

# Series add linearly
@test a+b == series(2.0, 1.0, -1.0)

# The output of series operations is truncated to match the input size.
# Represents (1+x-2x^2)(1 + x^2) = 1+x-x^2+o(x^3)
@test a*b == series(1.0,1.0,-1.0)

# Functions with known derivatives can easily be overloaded to operate on
# power series.
# You can generate the taylor series of a function about a point x up to
# e.g. 6th order by computing f(Series(x, 1.0, 0.0, 0.0, 0.0, 0.0))
x = series(0.0, 1.0, 0.0, 0.0, 0.0, 0.0)
@test sin(x) == series(0.0, 1.0, 0.0, -1.0/6.0, 0.0, 1.0/120)
@test log(1 + x) == series(0.0, 1.0, -1.0/2, 1.0/3, -1.0/4, 1.0/5)
@test 1/(1 - x) == series(1.0, 1.0, 1.0, 1.0, 1.0, 1.0)

@test polyval(a, 0.1) == 1.08
@test polyval(a, -0.1) == 0.88
@test polyder(a) == series(1.0, -4.0)
@test polyint(a) == series(0.0, 1.0, 1.0/2, -2.0/3)
@test restrict(a) == series(1.0, 1.0)

# Power series are one of the best ways to take higher order derivatives
# of generic functions.
#
# This example tests promotion rules in the constructor because of the
# 1 and 0, instead of 1.0 and 0.0
f(x) = exp(-x^2)
f2(x) = polyder(polyder(f(series(x, 1, 0))))
@test f2(2.0) == let x = 2.0; -2exp(-x^2)+4x^2*exp(-x^2); end

# macro-expansion-time conditional.
# This is probably not very good or general, but gets around the fact that
# the signture for @test_throws changed between versions 0.2 and 0.3
macro when(cond, ex)
  eval(cond) ? esc(ex) : nothing
end

# PowerSeries comes with types defined for series up to order 7. By default,
# trying to construct a higher order series is a type error.
@when(
  VERSION >= VersionNumber(0,2,typemax(Int64)),
  @test_throws MethodError series(0, 1, 2, 3, 4, 5, 6, 7, 8, 9)
)

# If you want to work with higher order series, you can generate types up
# to a given order with PowerSeries.generate(order)
PowerSeries.generate(9)
@test series(0, 1, 2, 3, 4, 5, 6, 7, 8, 9) != false

@test series(1, 2) - 1 == series(0, 2)

# Regression test for https://github.com/jwmerrill/PowerSeries.jl/issues/6
let x=series(0.0, 1.0, 0.0, 0.0)
  @test x*x == x^2
end

@test series(1.0, 2.0, 3.0) < 2.0
@test 2.0 > series(1.0, 2.0, 3.0)
@test 0.5 < series(1.0, 2.0, 3.0)
@test series(1.0, 2.0, 3.0) > 0.5
@test !(series(1.0, 2.0, 3.0) < 1.0)
@test !(series(1.0, 2.0, 3.0) > 1.0)
@test series(1.0, 1.0) < series(2.0, 0.0)
@test series(2.0, 0.0) > series(1.0, 1.0)
@test !(series(1.0, 2.0) < series(1.0, 0.0))
@test !(series(1.0, 0.0) > series(1.0, 2.0))

# Only 1 argument functions
fns = [
  sqrt,
  exp,
  log,
  sin,
  cos,
  tan,
  asin,
  acos,
  atan,
  sinh,
  cosh,
  tanh,
  asinh,
  acosh,
  atanh,
  csc,
  sec,
  cot,
  acsc,
  asec,
  acot,
  csch,
  sech,
  coth,
  acsch,
  asech,
  acoth,
  gamma,
  floor,
  ceil,
  round,
  sign,
  abs
]

for val in [-2.3, -1.2, -0.1, 0.7, 3.6]
  for fn in fns
    if (val < 1 && fn === acosh) continue end
    if (val < 0 && fn === sqrt) continue end
    if (val < 0 && fn === log) continue end
    if (abs(val) > 1 && fn === asin) continue end
    if (abs(val) > 1 && fn === acos) continue end
    if (abs(val) > 1 && fn === atanh) continue end
    if ((val < 0 || val > 1) && fn === asech) continue end
    if (abs(val) < 1 && fn === acsc) continue end
    if (abs(val) < 1 && fn === asec) continue end
    if (abs(val) < 1 && fn === acoth) continue end
    @test_approx_eq_eps constant(polyder(fn(series(val, 1.0)))) 5000*(fn(val + .0001) - fn(val- .0001)) .01
  end
end


text:
using PyCall
using PyPlot

@pyimport matplotlib.animation as anim  # After using PyPlot, so matplotlib installed via conda.
const wave_propagation_filename = "wave_propagation.mp4"

function saveall(output_path=".")
    fig = figure()
    ims = []
    for i_time in 1:10
        im = PyPlot.imshow(rand(4, 8), cmap="viridis")
        if 1 == i_time
            PyPlot.colorbar()
        end
        PyPlot.title("Wave amplitude [dB]")
        PyPlot.xlabel("Depth [m]")
        PyPlot.ylabel("Azimuth [m]")
        push!(ims, PyCall.PyObject[im])
    end
    #= close() =#

    # If matplotlib complains, ensure that
    # a) ffmpeg is installed with libx264 support, and
    # b) matplotlib is built with ffmpeg support enabled
    ani = anim.ArtistAnimation(fig, ims, interval=30, blit=true, repeat=false)
    anim_filename = wave_propagation_filename
    mywriter = anim.FFMpegWriter()
    ani[:save](joinpath(output_path, anim_filename), mywriter)
end

saveall()

max_stars_repo_name: cmey/matplotlib-animation-ci

max_issues_repo_name: cmey/matplotlib-animation-ci

max_forks_repo_name: cmey/matplotlib-animation-ci


text:
using JuMP
using Ipopt
using ForwardDiff

struct BoundaryConstraint # make a default setting if unbounded
    initialState::Array # make sure these are Nx1 size matricies
    finalState::Array # incorperate lower bound and ub
end

mutable struct TrajProblem
    objectiveFunc::Function
    dynamicsFunc::Function
    pathConstraintFunc::Function
    controlVectorGuess::Array
    stateVectorGuess::Array # each row represents a state, maybe use an add state guess function which adds each state guess
    timeStep::Union{Array, Float64} # if only one value input, turn this into an array of that one timestep, must be one less than state length
    boundaryConstraints::BoundaryConstraint # add checks that these are dimensionally correct
    numStates::Int   # can determine these from size of guess
    numControls::Int
    numCollocationPoints::Int # user / program specified
end

# work around to make current problem globally available
mutable struct CurrentProblem
    nullableProblem::Union{TrajProblem,Nothing}
end

const CURRENT_PROBLEM = CurrentProblem(nothing)
setCurrentProblem(problem::TrajProblem) = (CURRENT_PROBLEM.nullableProblem = problem)
getCurrentProblem() = CURRENT_PROBLEM.nullableProblem

function solve(problem::TrajProblem) # multiple dispatch on this function
    setCurrentProblem(problem)
    model = Model(with_optimizer(Ipopt.Optimizer, max_cpu_time=60.0)) # remove cpu time at some point

    # assign state and control varibles
    @variable(model, x[1:problem.numStates, 1:problem.numCollocationPoints])
    @variable(model, u[1:problem.numControls, 1:problem.numCollocationPoints])
    @constraint(model,x[1,:] >= 1.0)
    @constraint(model,x[3,:] >= 0.6)
    # set state and control variables to problem guess
    [set_start_value(x[i,j], problem.stateVectorGuess[i,j]) for i = 1:problem.numStates, j = 1:problem.numCollocationPoints]
    [set_start_value(u[i,j], problem.controlVectorGuess[i,j]) for i = 1:problem.numControls, j = 1:problem.numCollocationPoints]

    #register custom defined functions
    register(model, :collocateConstraint, (problem.numStates + problem.numControls) * problem.numCollocationPoints + 2, collocateConstraint, ΔcollocateConstraint)
    register(model, :pathConstraint, (problem.numStates + problem.numControls) * problem.numCollocationPoints + 1, pathConstraint, ΔpathConstraint)
    register(model, :objectiveFuncInterp, (problem.numStates + problem.numControls) * problem.numCollocationPoints, objectiveFuncInterp, ΔobjectiveFuncInterp)

    # state boundary constraints
    # initial state
    [fix(x[i,1], problem.boundaryConstraints.initialState[i]) for i = 1:length(problem.boundaryConstraints.initialState)]   # adjust this for lb and ub

    for i = 1:length(problem.boundaryConstraints.finalState)
        if !ismissing(problem.boundaryConstraints.finalState[i])
            fix(x[i,end], problem.boundaryConstraints.finalState[i])
        end
    end

    for i = 1:length(problem.pathConstraintFunc(problem.stateVectorGuess, problem.controlVectorGuess))
        @NLconstraint(model, pathConstraint(i,x...,u...) <=0)
    end

    for i = 1:problem.numStates
        for j = 1:problem.numCollocationPoints - 1
            @NLconstraint(model, collocateConstraint(i,j,x...,u...) == 0)
        end
    end

    @NLobjective(model, Min, objectiveFuncInterp(x...,u...))
    optimize!(model)
    return value.(x), value.(u)
end

function objectiveFuncInterp(stateControlVector...)
    problem = getCurrentProblem()
    # assemble state and control vector from tuple inputs
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)

    # return interpolated integral
    return sum(0.5 .* problem.timeStep[1] .* (problem.objectiveFunc(stateVector[2:end,:], controlVector[2:end]) .+ problem.objectiveFunc(stateVector[2:end,:], controlVector[2:end])) ) # should we use a map function here? test with a time, also needs to be rewritten for control vector with multiple dimensions
end

function objectiveFuncInterp(stateControlVector::Array)  # can we type union this? so we don't have to define two functions, only one
    problem = getCurrentProblem()
    # assemble state and control vector from tuple inputs
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)
    # return interpolated integral
    return [sum(0.5 .* problem.timeStep[1] .* (problem.objectiveFunc(stateVector[2:end,:], controlVector[2:end]) .+ problem.objectiveFunc(stateVector[1:end-1,:], controlVector[1:end-1])))]  # should we use a map function here? test with a time, also needs to be rewritten for control vector with multiple dimensions
end

function ΔobjectiveFuncInterp(g, stateControlVector...)
    problem = getCurrentProblem()
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)
    jacobian = ForwardDiff.jacobian(objectiveFuncInterp, [stateVector; controlVector])
    g[:] = getStateControlFromXU(problem, jacobian)
    return g
end

function collocateConstraint(ζr, ζc, stateControlVector...) # make sure timestep vector matches length of state vector
    ζr = convert(Int, ζr)
    ζc = convert(Int, ζc) # convert from jumps conversion to float64

    problem = getCurrentProblem()
    # assemble state and control vector from tuple inputs
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)

    ΔstateVector = 0.5 .* problem.timeStep[1] .* (problem.dynamicsFunc(stateVector[:,2:end], controlVector[:,2:end]) + problem.dynamicsFunc(stateVector[:,1:end-1], controlVector[:,1:end-1])) # how do we ensure dynamicsFunc has the right inputs and outputs if it is user defined? also have a look at making timestep dynamic for each value
    ζ = stateVector[:, 2:end] - stateVector[:,1:end-1] - ΔstateVector
    return ζ[ζr, ζc]
end

function collocateConstraint(stateControlVector::Array) # make sure timestep vector matches length of state vector
    problem = getCurrentProblem()
    # assemble state and control vector from tuple inputs
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)

    ΔstateVector = 0.5 .* problem.timeStep[1] .* (problem.dynamicsFunc(stateVector[:,2:end], controlVector[:,2:end]) + problem.dynamicsFunc(stateVector[:,1:end-1], controlVector[:,1:end-1])) # how do we ensure dynamicsFunc has the right inputs and outputs if it is user defined? also have a look at making timestep dynamic for each value
    ζ = stateVector[:, 2:end] - stateVector[:,1:end-1] - ΔstateVector
    return ζ
end

function ΔcollocateConstraint(g, ζr, ζc, stateControlVector...)
    problem = getCurrentProblem()
    g[1:2] .= 0 # gradient of xr and xc is zero as they're just indexes and will not be changed
    ζr = convert(Int, ζr)
    ζc = convert(Int, ζc) # convert from jumps conversion to float64
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)
    jacobian = ForwardDiff.jacobian(collocateConstraint, [stateVector; controlVector])
    jacobianRow = (ζc - 1) * problem.numStates + ζr # get the row of the jacobian that matches our output function
    g[3:end] = getStateControlFromXU(problem, jacobian[jacobianRow,:])
    return g
end

function pathConstraint(ζr, stateControlVector...) # make sure timestep vector matches length of state vector
    ζr = convert(Int, ζr)
    problem = getCurrentProblem()
    # assemble state and control vector from tuple inputs
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)

    ζ = problem.pathConstraintFunc(stateVector, controlVector)
    return ζ[ζr]
end

function pathConstraint(stateControlVector::Array) # make sure timestep vector matches length of state vector
    problem = getCurrentProblem()
    # assemble state and control vector from tuple inputs
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)
    ζ = problem.pathConstraintFunc(stateVector, controlVector)
    return ζ
end

function ΔpathConstraint(g, ζr,stateControlVector...)
    problem = getCurrentProblem()
    g[1] = 0 # gradient of xr and xc is zero as they're just indexes and will not be changed
    ζr = convert(Int, ζr)
    stateVector, controlVector = getXUFromStateControl(problem, stateControlVector)
    jacobian = ForwardDiff.jacobian(pathConstraint, [stateVector; controlVector])
    g[2:end] = getStateControlFromXU(problem, jacobian[ζr,:])
    return g
end

function getXUFromStateControl(problem, stateControlVector::Tuple) # get separate state and control vector matricies from input tuple, add type to problem
    stateControlVector = collect(stateControlVector)
    stateVector = zeros(problem.numStates, problem.numCollocationPoints)
    controlVector = zeros(problem.numControls, problem.numCollocationPoints)
    for i = 1:length(stateControlVector) ÷ (problem.numStates + problem.numControls)
        jx = (i-1) * (problem.numStates) + 1
        ju = (i-1) * (problem.numControls) + problem.numCollocationPoints * problem.numStates + 1
        stateVector[:,i] = stateControlVector[jx:jx+problem.numStates-1]
        controlVector[:,i] = stateControlVector[ju:ju+problem.numControls-1]
    end
    return stateVector, controlVector
end

function getXUFromStateControl(problem, stateControlVector::Array) # get separate state and control vector matricies from input tuple
    stateVector = Array{Union{Missing, typeof(stateControlVector[1]) }}(missing, problem.numStates, problem.numCollocationPoints)
    controlVector = Array{Union{Missing,typeof(stateControlVector[1]) }}(missing, problem.numControls, problem.numCollocationPoints)
    stateVector[:,:] = stateControlVector[1:problem.numStates, :]
    controlVector[:,:] = stateControlVector[problem.numStates + 1, :]
    return stateVector, controlVector
end

function getStateControlFromXU(problem, stateControlVector::Array) # get tuple of statecontrols from a 1D array with states and controls in a column row pattern
    x = ones(1,problem.numCollocationPoints*problem.numStates)
    u = ones(1,problem.numCollocationPoints*problem.numControls)
    for i = 1:problem.numCollocationPoints
        j = (i-1)*(problem.numStates + problem.numControls) + 1
        x[(i-1)*problem.numStates+1:i*problem.numStates] = stateControlVector[j:j+problem.numStates-1]
        u[(i-1)*problem.numControls+1:i*problem.numControls] = stateControlVector[j+problem.numStates:j+problem.numStates+problem.numControls-1]
    end
    return [x u]
end

####################################################################################################################################
### problem specific functions #####################################################################################################
####################################################################################################################################

function dynamicsFunc(stateVector, controlVector)
    h = stateVector[1,:]
    v = stateVector[2,:]
    m = stateVector[3,:]
    T = controlVector[1,:]
    D(h0,h,v) = 320 * v^2 * exp(-500*((h-h0)/h0))
    g(h0,h) = (h0/h)^2

    x1dot = h
    x2dot = map( (T,hi,v,m) -> (T-D(h[1],hi,v)) / m - g(h[1],hi), T,h,v,m)
    x3dot = -2T
    return [x1dot';x2dot';x3dot']
end
objectiveFunc(stateVector, controlVector) = -stateVector[1,end]

function pathConstraintFunc(stateVector, controlVector) # less than or equal to output value. Output must by 1 D vector of constraints
    T_max = 3.5
    T_min = 0
    h_0 = 1
    m_0 = 1
    m_f = 0.6

    bounds = zeros(typeof(stateVector[1,1]),6,size(stateVector,2))
    bounds[1,:] = -m_f .- stateVector[3,:] # weight constraints
    bounds[2,:] = -m_0 .+ stateVector[3,:]
    bounds[3,:] = -T_min .- controlVector[1,:] # thurst constraints
    bounds[4,:] = -T_max .+ controlVector[1,:]
    bounds[5,:] = -stateVector[2,:] # velocity constraint
    bounds[6,:] = h_0 .- stateVector[1,:] # velocity constraint
    return [bounds...]
end

####################################################################################################################################
### example code  ##################################################################################################################
####################################################################################################################################

numCP = 30
T = 0.2
controlVectorGuess = zeros(1,numCP)
timeStep = ones(1,numCP-1) * T/(numCP-1)
time = pushfirst!(cumsum(timeStep';dims=1)[:,1],0.0)
stateVectorGuess = ones((3,numCP))
stateVectorGuess[2,:] = time .* (1 .- time)
stateVectorGuess[3,:] = -0.4.*time.+1

boundaryConstraints = BoundaryConstraint([0;1;1],[missing;missing;0.6])
problem = TrajProblem(objectiveFunc, dynamicsFunc, pathConstraintFunc, controlVectorGuess, stateVectorGuess, timeStep, boundaryConstraints, 3, 1, numCP)
x,u = solve(problem)

using Plots
plotly()
plot(time, x[1,:])
plot(time, x[2,:])
plot(time, u')

max_stars_repo_name: sprockmonty/JuliaFYP

max_issues_repo_name: sprockmonty/JuliaFYP

max_forks_repo_name: sprockmonty/JuliaFYP


text: Dict{String,Any}()

max_stars_repo_name: vanillajonathan/julia

max_issues_repo_name: vanillajonathan/julia

max_forks_repo_name: vanillajonathan/julia


text:
a = [1 2 3; 4 5 6; 7 8 9]
b = [1 0 0; 0 2 0; 9 9 9]

@assert a + b == [2 2 3; 4 7 6; 16 17 18]
@assert a * 2 == [2 4 6; 8 10 12; 14 16 18]
@assert a * b == [28 31 27; 58 64 54; 88 97 81]  # Matrix multiply
@assert a .* b == [1 0 0; 0 10 0; 63 72 81]      # Elementwise mul

@assert -a == [-1 -2 -3; -4 -5 -6; -7 -8 -9]

max_stars_repo_name: idrougge/ple

max_issues_repo_name: idrougge/ple

max_forks_repo_name: idrougge/ple


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -9.75097970893739 <= q <= 8.436779502208893)

                     add_NL_constraint(m, :(swish(0.5456914178120993 + -0.994137029231358*swish(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.38451384631825425*swish(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + 0.9673084611668079*swish(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) + swish(0.6143810270772536 + 0.32603207961941205*swish(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.7660966505536861*swish(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + 0.30918160231413605*swish(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) + swish(-0.4079478457750749 + 0.06443443960646311*swish(0.4043679177252084 + -0.9159920740102909*$(x[1]) + 0.630569400327222*$(x[2]) + 0.15548796895771444*$(x[3])) + 0.09351478982242645*swish(-0.7218630242098185 + -0.7615826451198555*$(x[1]) + 0.556650647741431*$(x[2]) + -0.7074053991868854*$(x[3])) + -0.7258007123919756*swish(-0.4849043710003591 + 0.9895364493832828*$(x[1]) + -0.17212002157434236*$(x[2]) + -0.9444649970468406*$(x[3]))) - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "1432.max-difference-you-can-get-from-changing-an-integer.jl" begin
    @test max_score_1432([1,2,3,4,5,6,1], 3) == 12
    @test max_score_1432([2,2,2], 2) == 4
    @test max_score_1432([9,7,7,9,7,7,9], 7) == 55
    @test max_score_1432([1,79,80,1,1,1,200,1], 3) == 202
    @test max_score_1432([1,1000,1], 1) == 1
end


text:
function get_player_actions_dict(grid)
    dict = Dict{String,Array{Int32,1}}()
    player_actions = convert(Array{String,1},Griddly.get_player_available_actions(grid))
    for action_name in player_actions
        action_input = convert(Array{Int32,1},Griddly.get_input_ids(grid,action_name))
        dict[action_name] = action_input
    end
    return dict
end

function get_map_object_ids_char(grid)
    dict = Dict{Int32,String}()
    list_object_ids = convert(Array{Int32,1},Griddly.get_object_ids_list(grid))
    for object_id in list_object_ids
        object_char = convert(String,Griddly.get_object_char_from_id(grid,object_id))
        dict[object_id] = object_char
    end
    return dict
end

# Need the ObserverType to be ASCII
function pretty_print_ASCII(grid,mapping,observation;sep=" ")
    mapping[0] = "."
    _, width, height = size(observation)
    to_print = []
    for i in 0:(width-1)
        line = ""
        for j in 1:height
            line = string(line,sep,mapping[observation[i*width+j]])
        end
        println(line)
    end
    println("--------------------------")
end


text:
using Revise
using ZonalFlow
using JLD2

tspan   = (0.0,1000.0);
tsargs  = (
            dt=0.001,
            adaptive=false,
            progress=true,
            progress_steps=100000,
            save_everystep=false,
            dense=false,
            save_noise=false,
            saveat=5
           );

domain  = Domain(extent=(2π,2π),res=(8,8));
coeffs  = Coefficients(Ω=0.0,θ=0.0,μ=0.0,ν=0.02,ν₄=0.0); # β=0
forcing = Kolmogorov(A₁=-1.0,A₄=-2.0);
prob    = BetaPlane(domain,coeffs,forcing);

eqs     = [NL(),GQL(7),GCE2(7),GQL(0),CE2(),GCE2(Λ=0,poscheck=true,poscheckat=10),GQL(1),GCE2(1),GQL(3),GCE2(3)];
labels  = ["nl","gql_m","gce2_m","ql_ts5","ce2_ts5","gce2_0","gql_1","gce2_1","gql_3","gce2_3"];
# sols    = integrate(prob,eqs,tspan;tsargs...);

dn      = "validate/kolmogorov/8x8/"

nl      = integrate(prob,eqs[1],tspan;tsargs...)
write(prob,eqs[1],nl,dn=dn,fn=labels[1])
@save dn*labels[1]*".jld2" nl

gql_1      = integrate(prob,GQL(1),tspan;tsargs...)

gql_7      = integrate(prob,eqs[2],tspan;tsargs...)
write(prob,eqs[2],gql_7,dn=dn,fn=labels[2])
@save dn*labels[2]*".jld2" gql_7

gce2_7     = integrate(prob,eqs[3],tspan;tsargs...)
write(prob,eqs[3],gce2_7,dn=dn,fn=labels[3])
@save dn*labels[3]*".jld2" gce2_7

ql      = integrate(prob,eqs[4],tspan;u0=0.0,tsargs...)
write(prob,eqs[4],ql,dn=dn,fn=labels[4])
@save dn*labels[4]*".jld2" ql
@load dn*labels[4]*".jld2" ql

ce2     = integrate(prob,eqs[5],tspan;tsargs...)
write(prob,eqs[5],ce2,dn=dn,fn=labels[5])
@save dn*labels[5]*".jld2" ce2

gce2      = integrate(prob,eqs[6],tspan;tsargs...)
write(prob,eqs[6],gce2,dn=dn,fn=labels[6])
@save dn*labels[6]*".jld2" gce2

gql_1     = integrate(prob,eqs[7],tspan;tsargs...)
write(prob,eqs[7],gql_1,dn=dn,fn=labels[7])
@save dn*labels[7]*".jld2" gql_1

gce2_1    = integrate(prob,eqs[8],tspan;u0=gql_1.u[end],tsargs...)
write(prob,eqs[8],gce2_1,dn=dn,fn=labels[8])
@save dn*labels[8]*".jld2" gce2_1
# @load dn*labels[8]*".jld2" gce2_1

@info "Running:" eqs[9]
gql_3     = integrate(prob,eqs[9],tspan;tsargs...)
write(prob,eqs[9],gql_2,dn=dn,fn=labels[9])
@save dn*labels[9]*".jld2" gql_2

@info "Running:" eqs[10]
gce2_3     = integrate(prob,eqs[10],tspan;tsargs...)
write(prob,eqs[10],gce2_3,dn=dn,fn=labels[10])
@save dn*labels[10]*".jld2" gce2_3

# sol    = integrate(prob,NL(),tspan;tsargs...);
# write(prob,NL(),sol,dn="validate/kolmogorov/6x6/",fn="nl")
#
# sol    = integrate(prob,GQL(Λ=prob.d.nx-1),tspan;tsargs...);
# write(prob,GQL(Λ=prob.d.nx-1),sol,dn="validate/kolmogorov/6x6/",fn="gql")
#
# sol    = integrate(prob,GCE2(Λ=prob.d.nx-1),tspan;tsargs...);
# write(prob,GCE2(Λ=prob.d.nx-1),sol,dn="validate/kolmogorov/6x6/",fn="gce2")
#
# sol    = integrate(prob,GQL(0),tspan;tsargs...);
# write(prob,GQL(0),sol,dn="validate/kolmogorov/6x6/",fn="ql")
#
# sol    = integrate(prob,CE2(),tspan;u0=sol.u[end],tsargs...);
# write(prob,CE2(),sol,dn="validate/kolmogorov/6x6/",fn="ce2")
#
# sol    = integrate(prob,GCE2(0),tspan;tsargs...);
# write(prob,GCE2(0),sol,dn="validate/kolmogorov/6x6/",fn="gce2_0")
#
# sol    = integrate(prob,GQL(1),tspan;tsargs...);
# write(prob,GQL(0),sol,dn="validate/kolmogorov/6x6/",fn="gql_1")
#
# sol    = integrate(prob,GCE2(1),tspan;tsargs...);
# write(prob,GQL(0),sol,dn="validate/kolmogorov/6x6/",fn="gce2_1")


text:
using Bloomie
using Base.Test

# check blpapi version is higher than 3.8.0.0
@test (x->Int(1e6*x[1]+1e4*x[2]+1e2*x[3]+x[4]))(Bloomie.get_version())>=3080000


text:
#!/usr/bin/env julia

#= ID3v1.1 parser
   http://rubyquiz.com/quiz136.html

   Usage:
     ./script <mp3 file> [...]
=#

module ID3
using Printf

const tagsSize = 128

struct Tags
  song :: String
  artist :: String
  album  :: String
  year :: String
  comment :: String
  genre_index :: UInt8
  track_number :: Union{UInt8, Nothing}
end

function load_genres()
  open(joinpath(@__DIR__, "genres")) do f
    readlines(f)
  end
end

genres = load_genres()

function genre(t :: Tags) :: Union{String, Nothing}
  if 0 <= t.genre_index < length(genres)
    genres[t.genre_index+1]
  end
end

function read_tags_string(f) :: String
  seekend(f)
  offset = position(f)
  if offset >= tagsSize
    skip(f, -tagsSize)
    read(f, String)
  else
    ""
  end
end

function extract_field(tag_string :: String, interval) :: String
  string(strip(SubString(tag_string, interval), '\0'))
end

"""
Parse tags from a string and return either nothing or a Tags struct.
"""
function parse_tags(tag_string :: String) :: Union{Tags, Nothing}
  if length(tag_string) != tagsSize || SubString(tag_string, 1:3) != "TAG"
    return nothing
  end
  # 1       4     34       64      94     98       128     129
  # | "TAG" | song | artist | album | year | comment | genre |
  comment = SubString(tag_string, 98:127)
  track_number = nothing

  if comment[29] == '\0' && comment[30] != '\0'
    track_number = UInt8(comment[30])
    comment = SubString(comment, 1, 29)
  end

  Tags(
       extract_field(tag_string, 4:33),
       extract_field(tag_string, 34:63),
       extract_field(tag_string, 64:93),
       extract_field(tag_string, 94:97),
       string(strip(comment, '\0')),
       UInt8(tag_string[128]),
       track_number,
      )
end

function print_tags(filename :: String, tags :: Union{Tags, Nothing})
  @printf("%s\n", filename)
  if tags == nothing
    return println("  (no ID3 tags)")
  end
  @printf("  Song:   %s\n", tags.song)
  @printf("  Artist: %s\n", tags.artist)
  @printf("  Album:  %s\n", tags.album)
  if tags.track_number != nothing
    @printf("  Track number: %d\n", tags.track_number)
  end
  @printf("  Year:   %s\n", tags.year)
  @printf("  Genre:  %s\n", genre(tags))
  if !isempty(tags.comment)
    @printf("  Comment: %s\n", tags.comment)
  end
end

export read_tags_string, parse_tags, print_tags, load_genres, genre

end # module end

using Printf

if PROGRAM_FILE != "" && realpath(@__FILE__) == realpath(PROGRAM_FILE)
  using .ID3

  if isempty(ARGS)
    println("Usage:\n  $(PROGRAM_FILE) <mp3 file> [...]")
  else
    for filename in ARGS
      tag_string = try
        open(filename) do f
          read_tags_string(f)
        end
      catch e
        @printf(stderr, "Error reading file '%s': %s\n", filename, e)
        exit(1)
      end
      tags = parse_tags(tag_string)
      print_tags(filename, tags)
    end
  end
end

max_stars_repo_name: bfontaine/Katas

max_issues_repo_name: bfontaine/Katas

max_forks_repo_name: bfontaine/Katas


text:
@reexport module ModelKit

export @var,
    @unique_var,
    Variable,
    Expression,
    coefficients,
    degree,
    degrees,
    differentiate,
    dense_poly,
    evaluate,
    expand,
    exponents_coefficients,
    expressions,
    horner,
    is_homogeneous,
    multi_degrees,
    parameters,
    nparameters,
    nvariables,
    monomials,
    optimize,
    parameters,
    subs,
    support_coefficients,
    rand_poly,
    to_dict,
    to_number,
    variables,
    System,
    Homotopy,
    TaylorVector,
    variable_groups,
    vectors

import LinearAlgebra
import MultivariatePolynomials:
    MultivariatePolynomials,
    coefficients,
    degree,
    differentiate,
    nvariables,
    monomials,
    subs,
    variables
const MP = MultivariatePolynomials

include("model_kit/symengine.jl")
include("model_kit/symbolic.jl")
include("model_kit/instructions.jl")
include("model_kit/codegen.jl")

end # module


text:
#73
na = NamedArray([1, 2], ([1, missing],), ("A",))
let one = na[Name(1)], two = na[Name(missing)]
    @test one == 1
    @test two == 2
end

#39
include("init-namedarrays.jl")
v = n[1, :]
@test sin.(v).array == sin.(v.array)
@test namesanddim(sin.(v)) == namesanddim(v)

#110
n = NamedArray([1 2; 3 4], ([11,12], [13,14]))
s = sum(n, dims=1)
table = n ./ s
@test namesanddim(table) == namesanddim(n)

#105
a = [10, 20, 30]
n = NamedArray(a)
filter!(x -> x > 15, n)
@test n == a == [20, 30]
@test names(n, 1) == ["2", "3"]


text:
import DataFrames

function table(df::DataFrames.DataFrame; args...) :: String
  grh.table(; class="table", args...) do
    grh.thead() do
      grh.tr() do
        grh.collection(names(df)) do item
          grh.th(scope="col") do
            item |> string
          end
        end
      end
    end *
    grh.tbody() do
      grh.collection(eachrow(df) |> collect) do row
        grh.tr() do
          grh.collection(Array(row)) do item
            grh.td() do
              item |> string
            end
          end
        end
      end
    end
  end
end


text:
module TestVarArgs
using Test

f() = println("()")
f(a::Vararg{T,N}) where {T,N} = println("v", tuple(typeof(a).types...))
f(a::Tuple) =  println(tuple(typeof(a).types...))
f(a...) = println(tuple(typeof(a).types...), "...")

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx