text:
#!/usr/bin/env julia

#compute energy of a single sequence
#h=(q,N) and J=(q,N)
function compute_energy_single_sequence(h::Array{Float64,2},
                                        J::Array{Float64,4},
                                        S::Vector)

    N = size(h)[2]
    q = size(h)[1]
    E = 0.0
    for i = 1:N
        E -= h[S[i],i]
        for j = (i+1):N
			E -= J[S[i],S[j],i,j]
		end
	end
return E
end

function make_random_mutations_compute_energy(h, J, num_muts, ref_seq, length_ref_seq, freqs_pseudo)

		unsorted_columns = randperm(length_ref_seq)
		j = 1
		muts_done = 0
		mutated_sequence = copy(ref_seq)
		while muts_done <= num_muts
			column = unsorted_columns[j]
			if ref_seq[column] != 21 # If it is gap, we not consider the position
				unsorted_residues = randperm(20)
				if ref_seq[column] == unsorted_residues[1] # If the first random residue is the same as in the wild type, take the second. Take the first otherwise
					mut = unsorted_residues[2]
				else
					mut = unsorted_residues[1]
				end
				mutated_sequence[column] = mut
				muts_done = muts_done + 1
				# println("muts done "*string(muts_done))
			end
			j = j + 1
			# println(j)
		end
	return compute_energy_single_sequence(h,J,mutated_sequence)
end

function make_random_mutations_single_sequence(num_muts, ref_seq)

	unsorted_columns = randperm(length(ref_seq))
	j = 1
	muts_done = 0
	mutated_sequence = copy(ref_seq)
	while muts_done <= num_muts
		column = unsorted_columns[j]
		if ref_seq[column] != 21 # If it is gap, we not consider the position
			unsorted_residues = randperm(20)
			if ref_seq[column] == unsorted_residues[1] # If the first random residue is the same as in the wild type, take the second. Take the first otherwise
				mut = unsorted_residues[2]
			else
				mut = unsorted_residues[1]
			end
			mutated_sequence[column] = mut
			muts_done = muts_done + 1
		end
			j = j + 1
	end
	return mutated_sequence
end

function make_random_mutations_alignment(mutations, ref_seq)

	num_seqs = length(mutations[:,1])
	tam = length(ref_seq)
	alignment = Array{Int8}(tam, num_seqs)

	for seq in 1:num_seqs
		num_muts = length(mutations[seq,1])

		mutated_sequence = make_random_mutations_single_sequence(num_muts, ref_seq)
		for j in 1:tam
			alignment[j, seq] = mutated_sequence[j]
		end
	end
	return alignment
end

function compute_energy_independent_model(freqs, seq)
	energy_ind = 0
	for col in 1:length(freqs[1,:])
		aa = seq[col]
		energy_ind += log(freqs[aa,col])
	end
	return -energy_ind
end

function compute_diff_independent_only_muts(freqs, ref_seq, seq)
	diff_ind = 0
	for col in 1:length(seq)
		aa = seq[col]
		aa_ref_seq = ref_seq[col]
		if aa != aa_ref_seq
			diff_ind += log(freqs[aa,col]) - log(freqs[aa_ref_seq,col])
		end
	end
	return -diff_ind
end

function compute_diff_independent_only_muts_sequences(freqs, ref_seq, sequences)
	num_seqs = length(sequences[1,:])
	diff_ind = Vector(num_seqs)
	for i in 1:num_seqs
		seq = sequences[:,i]
		diff_ind[i] = compute_diff_independent_only_muts(freqs, ref_seq, seq)
	end
	return diff_ind
end

function compute_energy_independent_model_sequences(freqs, sequences)
	num_seqs = length(sequences[1,:])
	energy_ind = Vector(num_seqs)
	for i in 1:num_seqs
		seq = sequences[:,i]
		energy_ind[i] = compute_energy_independent_model(freqs, seq)
	end
	return energy_ind
end

function add_pseudocounts(freqs, alpha)
	tam_i = length(freqs[1,:])
	tam_j = length(freqs[:,1])
# 	freqs_pseudo = Array{Float64,2}(tam_j,tam_i)
    freqs_pseudo = zeros(Float64, tam_j,tam_i)
	for i in 1:tam_i
		for j in 1:tam_j
			freqs_pseudo[j,i] = (1-alpha)*freqs[j,i] + alpha/21
		end
	end

	return freqs_pseudo
end

# function make_mutation_according_to_IND(ref_seq, pos, freqs)
# 	aa = ref_seq[pos]
# 	mut_aa = aa
#
# 	if aa != 21 # If it is a gap, we do not change anything, just return the gap
# 		random_num = rand()
# 		found = 0
# 		while(!found)
#
# 			found = 1
# 		end
# 	end
#
# 	return mut_aa
# end

function normalized_cummulative_without_gaps(freqs)
	num_cols = length(freqs[1,:])
	normalized_freqs = Array{Float64,2}(20, num_cols)
	cummulative_normalized_freqs = Array{Float64,2}(20, num_cols)
	for col in 1:num_cols
		sum_freq = 0
		for aa in 1:20 # to avoid gap = 21
			sum_freq = sum_freq + freqs[aa,col]
		end
		# Normalize without gap
		for aa in 1:20 # to avoid gap = 21
			normalized_freqs[aa,col] = freqs[aa,col]/sum_freq
		end
		cummulative_normalized_freqs[1,col] = normalized_freqs[1,col]
		for aa in 2:20 # to avoid gap = 21
			cummulative_normalized_freqs[aa,col] = cummulative_normalized_freqs[aa-1,col]+normalized_freqs[aa,col]
		end
	end

	return cummulative_normalized_freqs
end

function get_mutations(ali, ref_seq)
	num_seqs = length(ali[1,:])
	num_pos = length(ali[:,1])

	mutations = []
	for seq in 1:num_seqs
		mutations_this_seq = []
		for pos in 1:num_pos
			if ali[pos, seq] != ref_seq[pos]
				push!(mutations_this_seq, [pos, ali[pos, seq]])
			end
		end
		push!(mutations, mutations_this_seq)
	end

	return mutations
end

function compute_delta_energy_mut_separately(mutations, ref_seq)
	num_seqs = length(mutations[:,1])
	diff_energy = zeros(num_seqs)
	for seq in 1:num_seqs
		num_muts = length(mutations[seq,1])
		sum_delta = 0
		for mut in 1:num_muts
			mutation = mutations[seq,1][mut]
			mut_seq = copy(ref_seq)
			mut_pos = mutation[1]
			mut_aa = mutation[2]
			mut_seq[mut_pos] = mut_aa
			energy_mut = compute_energy_single_sequence(h,J,mut_seq)
			delta = energy_mut - energy_ref
			sum_delta = sum_delta + delta
		end
		diff_energy[seq] = sum_delta
	end
	return diff_energy
end

function compute_delta_energy_mut_return_separately(mutations, ref_seq)
	num_seqs = length(mutations[:,1])
	diff_energy = zeros(num_seqs)

	effects = []
	for seq in 1:num_seqs
		num_muts = length(mutations[seq,1])
		sum_delta = 0
		mutations_this_seq = []
		for mut in 1:num_muts
			mutation = mutations[seq,1][mut]
			mut_seq = copy(ref_seq)
			mut_pos = mutation[1]
			mut_aa = mutation[2]
			mut_seq[mut_pos] = mut_aa
			energy_mut = compute_energy_single_sequence(h,J,mut_seq)
			delta = energy_mut - energy_ref
			sum_delta = sum_delta + delta
			push!(mutations_this_seq, delta)
		end
		diff_energy[seq] = sum_delta
		push!(effects, mutations_this_seq)
	end
	return effects
end

max_stars_repo_name: juan-rodriguez-rivas/covmut

max_issues_repo_name: juan-rodriguez-rivas/covmut

max_forks_repo_name: juan-rodriguez-rivas/covmut


text:
module LsoBase

using DataFrames

"""
    includedir(dir)

    Include all non-hidden .jl files from the given directory
"""
function includedir(dir::AbstractString)
    for filename in readdir(dir)
        if isfile(dir*"/"*filename) && endswith(filename, ".jl") && !startswith(filename, ".")
            include(dir*"/"*filename)
        end
    end
end

"""
    new_inf()

    Create new DataFrame to hold iteration info.
    Will have columns f, w, opt, iter, lsiter and time.
"""
new_inf() = DataFrame(
                w      = Array{Float64, 1}[],
                f      = Float64[],
                opt    = Float64[],
                iter   = Int[],
                lsiter = Int[],
                time   = Float64[]
            )

"""
    new_acc()

    Create new DataFrame to hold classifier info.
"""
new_acc() = DataFrame(
                gdopt      = ASCIIString[],
                ls         = ASCIIString[],
                batchsize  = Int[],
                seed       = Int[],
                data       = ASCIIString[],
                classifier = ASCIIString[],
                trainsize  = Int[],
                trainacc   = Float64[],
                testacc    = Float64[],
                time       = Float64[],
                iterrate   = Float64[]
            )

"""
    push_inf!(inf, w, opt, iter, lsiter)

    Will push the given iteration info to DataFrame inf.
"""
function push_inf!(inf::DataFrame, w::Array{Float64, 1},
                   f::Float64, opt::Float64, iter::Int, lsiter::Int, time::Float64=0.0)
    push!(inf, (w, f, opt, iter, lsiter, time))
end

function push_truth!(inf::DataFrame, f::Function, g::Function)
    println("\nComputing true f and opt values...")
    inf[:f_true]   = [ f(w)               for w in inf[:w] ]
    inf[:opt_true] = [ vecnorm(g(w), Inf) for w in inf[:w] ]
end

"""
    push_acc!(acc, gdopt, ls, batchsize, seed, data, trainsize, trainacc, testacc, time, iterrate)

    Will push the given classifier info to DataFrame acc.
"""
function push_acc!(acc::DataFrame, gdopt::ASCIIString, ls::ASCIIString, batchsize::Int, seed::Int,
                   data::ASCIIString, classifier::ASCIIString,
                   trainsize::Int, trainacc::Float64, testacc::Float64, time::Float64, iterrate::Float64)
    push!(acc, (gdopt, ls, batchsize, seed, data, classifier, trainsize, trainacc, testacc, time, iterrate))
end

"""
    store(filename, df)

    Store the given DataFrame, e.g, inf or acc, to the given location.
"""
store(df::DataFrame, filename::ASCIIString) = writetable(filename, df)

acc(y_true::Array{Float64, 1}, y_pred::Array{Float64, 1}) = mean(y_true .== y_pred)

end

ext: jl

max_stars_repo_name: mirkobunse/lso

max_issues_repo_name: mirkobunse/lso

max_forks_repo_name: mirkobunse/lso


text:
#=
    Vector valued constraints
=#
constrdict(m::LinQuadOptimizer, ::VLCI{MOI.Nonnegatives}) = cmap(m).nonnegatives
constrdict(m::LinQuadOptimizer, ::VLCI{MOI.Nonpositives}) = cmap(m).nonpositives
constrdict(m::LinQuadOptimizer, ::VLCI{MOI.Zeros})        = cmap(m).zeros

function MOI.addconstraint!(m::LinQuadOptimizer, func::VecLin, set::S) where S <: VecLinSets
    @assert MOI.dimension(set) == length(func.constants)

    nrows = get_number_linear_constraints(m)
    addlinearconstraint!(m, func, backend_type(m,set))
    nrows2 = get_number_linear_constraints(m)

    m.last_constraint_reference += 1
    ref = VLCI{S}(m.last_constraint_reference)

    dict = constrdict(m, ref)
    dict[ref] = collect(nrows+1:nrows2)
    for i in 1:MOI.dimension(set)
        push!(m.constraint_primal_solution, NaN)
        push!(m.constraint_dual_solution, NaN)
        push!(m.constraint_constant, func.constants[i])
    end
    ref
end

function addlinearconstraint!(m::LinQuadOptimizer, func::VecLin, sense::Cchar)
    outputindex   = [term.output_index for term in func.terms]
    columns       = [getcol(m, term.scalar_term.variable_index) for term in func.terms]
    coefficients  = [term.scalar_term.coefficient for term in func.terms]
    # sort into row order
    pidx = sortperm(outputindex)
    permute!(columns, pidx)
    permute!(coefficients, pidx)

    # check that there is at least a RHS for each row
    @assert maximum(outputindex) <= length(func.constants)
    # loop through to get starting position of each row
    row_pointers = Vector{Int}(length(func.constants))
    row_pointers[1] = 1
    row = 1
    for i in 2:length(pidx)
        if outputindex[pidx[i]] != outputindex[pidx[i-1]]
            row += 1
            row_pointers[row] = i
        end
    end
    A = CSRMatrix{Float64}(row_pointers, columns, coefficients)
    add_linear_constraints!(m, A, fill(sense, length(func.constants)), -func.constants)
end

MOI.canmodify(m::LinQuadOptimizer, ::Type{VLCI{S}}, ::Type{MOI.VectorConstantChange{Float64}}) where S <: VecLinSets = true
function MOI.modify!(m::LinQuadOptimizer, ref::VLCI{<: VecLinSets}, chg::MOI.VectorConstantChange{Float64})
    @assert length(chg.new_constant) == length(m[ref])
    for (r, v) in zip(m[ref], chg.new_constant)
        change_rhs_coefficient!(m, r, -v)
        m.constraint_constant[r] = v
    end
end

MOI.canmodify(m::LinQuadOptimizer, ::Type{VLCI{S}}, ::Type{MOI.MultirowChange{Float64}}) where S <: VecLinSets = true
function MOI.modify!(m::LinQuadOptimizer, ref::VLCI{<: VecLinSets}, chg::MOI.MultirowChange{Float64})
    col = m.variable_mapping[chg.variable]
    for (row, coef) in chg.new_coefficients
        change_matrix_coefficient!(m, row, col, coef)
    end
end

MOI.candelete(m::LinQuadOptimizer, c::VLCI{<:VecLinSets}) = MOI.isvalid(m, c)
function MOI.delete!(m::LinQuadOptimizer, c::VLCI{<:VecLinSets})
    deleteconstraintname!(m, c)
    dict = constrdict(m, c)
    # we delete rows from largest to smallest here so that we don't have
    # to worry about updating references in a greater numbered row, only to
    # modify it later.
    for row in sort(dict[c], rev=true)
        delete_linear_constraints!(m, row, row)
        deleteat!(m.constraint_primal_solution, row)
        deleteat!(m.constraint_dual_solution, row)
        deleteat!(m.constraint_constant, row)
        # shift all the other references
        shift_references_after_delete_affine!(m, row)
    end
    delete!(dict, c)
end

#=
    Constraint set of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintSet, ::Type{VLCI{S}}) where S <: VecLinSets = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintSet, c::VLCI{S}) where S <: VecLinSets
    constraint_indices = m[c]
    S(length(constraint_indices))
end

#=
    Constraint function of Linear function
=#

MOI.canget(m::LinQuadOptimizer, ::MOI.ConstraintFunction, ::Type{VLCI{S}}) where S <: VecLinSets = true
function MOI.get(m::LinQuadOptimizer, ::MOI.ConstraintFunction, c::VLCI{<: VecLinSets})
    ctrs = m[c]
    n = length(ctrs)
    constants = Float64[]
    terms = MOI.VectorAffineTerm{Float64}[]
    for i in 1:n
        rhs = get_rhs(m, ctrs[i])
        push!(constants, -rhs)
        # TODO more efficiently
        colidx, coefs = get_linear_constraint(m, ctrs[i])
        for (column, coefficient) in zip(colidx, coefs)
            push!(terms, MOI.VectorAffineTerm{Float64}(
                    i,
                    MOI.ScalarAffineTerm{Float64}(
                        coefficient, m.variable_references[column]
                    )
                )
            )
        end
    end
    return MOI.VectorAffineFunction(terms, constants)
end


text:
using Surrogates

#1D
lb = 0.0
ub = 5.0
obj_1D = x -> log(x)*exp(x)
x = sample(5,lb,ub,SobolSample())
y = obj_1D.(x)
my_second_order_poly = SecondOrderPolynomialSurrogate(x,y,lb,ub)
val = my_second_order_poly(5.0)
add_point!(my_second_order_poly,5.0,238.86)
add_point!(my_second_order_poly,[6.0,7.0],[722.84,2133.94])

#ND
lb = [0.0,0.0]
ub = [10.0,10.0]
obj_ND = x -> log(x[1])*exp(x[2])
x = sample(5,lb,ub,UniformSample())
y = obj_ND.(x)
my_second_order_poly = SecondOrderPolynomialSurrogate(x,y,lb,ub)
val = my_second_order_poly((5.0,7.0))
add_point!(my_second_order_poly,(5.0,7.0),1764.96)
add_point!(my_second_order_poly,[(1.5,1.5),(3.4,5.4)],[1.817,270.95])


text:
#___  ____       ___ _____   _   _       _ _
#|  \/  | |     |_  /  __ \ | | | |     (_| |
#| .  . | |       | | /  \/ | | | |_ __  _| |_ ___
#| |\/| | |       | | |     | | | | '_ \| | __/ _ \
#| |  | | |___/\__/ | \__/\ | |_| | | | | | || (_) |
#_______\_____\____/ \____/  _____|_|___|_____\_____ _____ _____
#| ___ \        (_)         | | \ \ / / / __  |  _  / __  |  _  |
#| |_/ _ __ ___  _  ___  ___| |_ \ V /  `' / /| |/' `' / /| |/' |
#|  __| '__/ _ \| |/ _ \/ __| __|/   \    / / |  /| | / / |  /| |
#| |  | | | (_) | |  __| (__| |_/ /^\ \ ./ /__\ |_/ ./ /__\ |_/ /
#\_|  |_|  \___/| |\___|\___|\__\/   \/ \_____/\___/\_____/\___/
#              _/ |
#             |__/
#
# This code is part of the proposal of the team "MLJC UniTo" - University of Turin
# for "ProjectX 2020" Climate Change for AI.
# The code is licensed under MIT 3.0
# Please read readme or comments for credits and further information.

# Compiler: Julia 1.5.3

# Short description of this file: 2D REACTION DIFFUSION EQUATION SOLVED WITH A CNN

#Warning: Reset Julia Kernel (pre-compilation required) to avoid errors in
#redefinitions of constants

using Pkg#; Pkg.activate("."); Pkg.instantiate()
Pkg.add("BSON")
Pkg.add("DifferentialEquations")

using PyPlot, Printf
using LinearAlgebra
using Flux, DiffEqFlux, Optim, DiffEqSensitivity
using BSON: @save, @load
using Flux: @epochs
using DifferentialEquations

#Domain definition
X = 1.0; T = 5.0;
dx = 0.04; dt = T/10;
Y = 1.0;
dy = 0.04;
x = collect(0:dx:X);
y = collect(0:dy:Y);
t = collect(0:dt:T);
Nx = 128 # Int64(X/dx+1);
Ny = 128 #Int64(Y/dy+1);
Nt = Int64(T/dt+1);

########################
# Generate training data
########################

const α₂ = 1.0
const α₃ = 1.0
const β₁ = 1.0
const β₂ = 1.0
const β₃ = 1.0
const r₁ = 1.0
const r₂ = 1.0
const DD = 100.0
const γ₁ = 0.1
const γ₂ = 0.1
const γ₃ = 0.1
const N = 128
const XX = reshape([i for i in 1:N for j in 1:N],N,N)
const YY = reshape([j for i in 1:N for j in 1:N],N,N)
const aaa = 1.0.*(XX.>=4*N/5)

const Mx = Tridiagonal([1.0 for i in 1:N-1],[-2.0 for i in 1:N],[1.0 for i in 1:N-1])
const My = copy(Mx)
Mx[2,1] = 2.0
Mx[end-1,end] = 2.0
My[1,2] = 2.0
My[end,end-1] = 2.0

# Define the discretized PDE as an ODE function
const MyA = zeros(N,N)
const AMx = zeros(N,N)
const DA = zeros(N,N)

function f(du,u,p,t)
   A = @view  u[:,:,1]
   B = @view  u[:,:,2]
   C = @view  u[:,:,3]
  dA = @view du[:,:,1]
  dB = @view du[:,:,2]
  dC = @view du[:,:,3]
  mul!(MyA,My,A)
  mul!(AMx,A,Mx)
  @. DA = DD*(MyA + AMx)
  @. dA = DA + aaa - β₁*A - r₁*A*B + r₂*C
  @. dB = α₂ - β₂*B - r₁*A*B + r₂*C
  @. dC = α₃ - β₃*C + r₁*A*B - r₂*C
end

u0 = zeros(N,N,3)

using BenchmarkTools
prob = ODEProblem(f,u0,(0.0,T))
#@btime sol = solve(prob, ROCK4(), reltol = 1e-8, abstol=1e-8, saveat=dt);
sol = solve(prob, ROCK4(), reltol = 1e-8, abstol=1e-8, saveat=dt); # 457.078 ms (3195 allocations: 238.33 MiB)

ode_data = Array(sol);

########################
# Define the neural PDE
########################
n_weights = 10

#for the reaction term
rx_nn = Chain(Dense(2, n_weights, tanh),
                Dense(n_weights, 2*n_weights, tanh),
                Dense(2*n_weights, n_weights, tanh),
                Dense(n_weights, 2))

#conv with bias with initial values as 1/dx^2
w_err = 0.0
init_w = reshape([1.1 -2.5 1.0], (3, 1, 1, 1))
diff_cnn_ = Conv(init_w, [0.], pad=(0,0,0,0))

#initialize D0 close to D/dx^2
D0 = [6.5]

p1,re1 = Flux.destructure(rx_nn)
p2,re2 = Flux.destructure(diff_cnn_)
p1
p2
p = [p1;p2;D0]
full_restructure(p) = re1(p[1:length(p1)]), re2(p[(length(p1)+1):end-1]), p[end]
p = p .*p'

function nn_ode(u,p,t)
    rx_nn = re1(p[1:length(p1)])

    u_cnn_1   = [p[end-4,end-4] * u[end,end,k] + p[end-3,end-3] * u[1,1,k] + p[end-2,end-2] * u[2,2,k] for k in 1:3]
    u_cnn     = [p[end-4,end-4] * u[i-1,j-1,k] + p[end-3,end-3] * u[i,j,k] + p[end-2,end-2] * u[i+1,j+1,k] for i in 2:Nx-1, j in 2:Ny-1, k in 1:3]
    u_cnn_end = [p[end-4,end-4] * u[end-1,end-1,k] + p[end-3,end-3] * u[end,end,k] + p[end-2,end-2] * u[1,1,k] for k in 1:3]

    [rx_nn([u[i,j,k], u[j,i,k]])[1] for i in 1:Nx, j in 1:Ny, k in 1:3] + p[end][end] * cat(reshape(u_cnn_1, (1, 1, 3)), u_cnn, reshape(u_cnn_end, (1, 1, 3)); dims = (1,2))
end

########################
# Solving the neural PDE and setting up loss function
########################

prob_nn = ODEProblem(nn_ode, u0, (0.0, T), p)

sol_nn = concrete_solve(prob_nn,VCABM(), u0, p, sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP()) )

function predict_rd(θ)
  # No ReverseDiff if using Flux
    Array(concrete_solve(prob_nn,VCABM(),u0,θ,saveat=dt,sensealg=InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end

#match data and force the weights of the CNN to add up to zero
function loss_rd(p)
    pred = predict_rd(p)
    sum(abs2, ode_data .- pred) + 10^2 * abs(sum(p[end-4 : end-2, end-4:end-2])), pred
end

########################
# Training
########################

#Optimizer
opt = ADAM(0.05)

global count = 0
global save_count = 0
save_freq = 50

train_arr = Float64[]
diff_arr = Float64[]
w1_arr = Float64[]
w2_arr = Float64[]
w3_arr = Float64[]

res1 = DiffEqFlux.sciml_train(loss_rd, p, opt, maxiters = 1, progress = true)

res2 = DiffEqFlux.sciml_train(loss_rd, res1.minimizer, opt, cb=cb, maxiters = 300)
res3 = DiffEqFlux.sciml_train(loss_rd, res2.minimizer, BFGS(), cb=cb, maxiters = 1000)

pstar = res3.minimizer


text:
export countmap, cucountmap, cucountmap!

using CUDA
using StatsBase
import StatsBase: countmap

cucountmap(v) = cu(v) |> countmap

# The winner of bencharmks/bencharmks-countmap
function cucountmap!(buffer, v::CuDeviceArray{T, 1, NN2}) where {T, NN2} #(buffer::CuArray{UInt32, 1, NN1}, v::CuArray{T, 1, NN2}) where {NN1, T, NN2}
    i = (blockIdx().x - 1) * blockDim().x + threadIdx().x
    stride = gridDim().x * blockDim().x

    addnum = 1 - Int(typemin(T))
    for j = i:stride:length(v)
        b = Int(v[j]) + addnum
        CUDA.@atomic buffer[b] = buffer[b] + 1
    end
    return
end

function countmap(v::CuArray{T}; threads = 256, blocks = 1024) where {T}
    st = sizeof(T)
    @assert st <= 2
    buffer = CUDA.zeros(Int, 2^(8st))
    CUDA.@sync @cuda threads = threads blocks = blocks cucountmap!(buffer, v)
    values = typemin(T) : typemax(T)
    Dict(zip(values, collect(buffer)))
end


text:
struct Environments
    I::MPS
    H::MPS
    InProgress::Matrix{ITensor}
end

function buildEdgeEnvironment(A::PEPS, H, left_H_terms, next_combiners, side::Symbol, col::Int; kwargs...)::Environments
    Ny, Nx              = size(A)
    is_cu::Bool         = is_gpu(A)
    up_combiners        = Vector{ITensor}(undef, Ny-1)
    fake_next_combiners = Vector{ITensor}(undef, Ny)
    fake_prev_combiners = is_cu ? fill(cuITensor(1.0), Ny) : fill(ITensor(1.0), Ny)
    I_mpo, fake_next_combiners, up_combiners = buildNewI(A, col, fake_prev_combiners, side)
    I_mps          = MPS(Ny, tensors(I_mpo), 0, Ny+1)
    @debug "Built new I"
    copyto!(next_combiners, fake_next_combiners)
    field_H_terms  = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms   = getDirectional(vcat(H[:, col]...), Vertical)
    vHs            = [buildNewVerticals(A, fake_prev_combiners, fake_next_combiners, up_combiners, vert_H_terms[vert_op], col) for vert_op in 1:length(vert_H_terms)]
    @debug "Built new Vs"
    fHs            = [buildNewFields(A, fake_prev_combiners, fake_next_combiners, up_combiners, field_H_terms[field_op], col) for field_op in 1:length(field_H_terms)]
    Hs             = MPS[MPS(Ny, tensors(H_term), 0, Ny+1) for H_term in vcat(vHs, fHs)]
    @inbounds for row in 1:Ny-1
        ci = linkindex(I_mps, row)
        ni = Index(ITensors.dim(ci), "u,Link,c$col,r$row")
        replaceindex!(I_mps[row], ci, ni)
        replaceindex!(I_mps[row+1], ci, ni)
    end
    maxdim::Int     = get(kwargs, :maxdim, 1)
    cutoff::Float64 = get(kwargs, :cutoff, 0.0)
    env_maxdim::Int = get(kwargs, :env_maxdim, maxdim)
    H_overall       = sum(Hs; cutoff=cutoff, maxdim=env_maxdim)
    @debug "Summed Hs, maxdim=$maxdim"
    side_H       = side == :left ? H[:, col] : H[:, col - 1]
    side_H_terms = getDirectional(vcat(side_H...), Horizontal)
    @debug "Trying to alloc $(Ny*length(side_H_terms))"
    in_progress  = Matrix{ITensor}(undef, Ny, length(side_H_terms))
    @inbounds for side_term in 1:length(side_H_terms)
        @debug "Generating edge bonds for term $side_term"
        in_progress[1:Ny, side_term] = generateEdgeDanglingBonds(A, up_combiners, side_H_terms[side_term], side, col)
    end
    @debug "Generated edge bonds"
    return Environments(I_mps, H_overall, in_progress)
end

function buildNextEnvironment(A::PEPS, prev_Env::Environments, H, previous_combiners::Vector{ITensor}, next_combiners::Vector{ITensor}, side::Symbol, col::Int; kwargs...)
    Ny, Nx = size(A)
    working_combiner = Vector{ITensor}(undef, Ny)
    up_combiners     = Vector{ITensor}(undef, Ny-1)
    @timeit "build I" begin
        I_mpo, working_combiner, up_combiners = buildNewI(A, col, previous_combiners, side)
    end
    @debug "Built I_mpo"
    copyto!(next_combiners, working_combiner)
    @inbounds for row in 1:Ny-1
        ci = linkindex(I_mpo, row)
        ni = Index(ITensors.dim(ci), "u,Link,c$col,r$row")
        replaceindex!(I_mpo[row], ci, ni)
        replaceindex!(I_mpo[row+1], ci, ni)
    end
    cutoff::Float64 = get(kwargs, :cutoff, 1e-13)
    maxdim::Int     = get(kwargs, :maxdim, 1)
    env_maxdim::Int = get(kwargs, :env_maxdim, maxdim)
    @timeit "build new_I and new_H" begin
        new_I     = applyMPO(I_mpo, prev_Env.I; cutoff=cutoff, maxdim=env_maxdim)
        new_H     = applyMPO(I_mpo, prev_Env.H; cutoff=cutoff, maxdim=env_maxdim)
    end
    @debug "Built new I and H"
    field_H_terms = getDirectional(vcat(H[:, col]...), Field)
    vert_H_terms  = getDirectional(vcat(H[:, col]...), Vertical)
    hori_H_terms  = getDirectional(vcat(H[:, col]...), Horizontal)
    side_H        = side == :left ? H[:, col] : H[:, col - 1]
    side_H_terms  = getDirectional(vcat(side_H...), Horizontal)
    H_term_count  = 1 + length(field_H_terms) + length(vert_H_terms)
    H_term_count += (side == :left ? length(side_H_terms) : length(hori_H_terms))
    final_H       = deepcopy(new_H)
    new_H_mps     = Vector{MPS}(undef, H_term_count)
    new_H_mps[1]  = deepcopy(new_H)
    @timeit "build new verts" begin
        vHs = [buildNewVerticals(A, previous_combiners, next_combiners, up_combiners, vert_H_terms[vert_op], col) for vert_op in 1:length(vert_H_terms)]
    end
    @debug "Built new Vs"
    @timeit "build new fields" begin
        fHs = [buildNewFields(A, previous_combiners, next_combiners, up_combiners, field_H_terms[field_op], col) for field_op in 1:length(field_H_terms)]
    end
    @timeit "build new H array" begin
        for (vv, vH) in enumerate(vHs)
            new_H_mps[1+vv] = applyMPO(vH, prev_Env.I; cutoff=cutoff, maxdim=env_maxdim)
        end
        for (ff, fH) in enumerate(fHs)
            new_H_mps[1+length(vHs)+ff] = applyMPO(fH, prev_Env.I; cutoff=cutoff, maxdim=env_maxdim)
        end
    end
    connect_H    = side == :left ? side_H_terms : hori_H_terms
    @timeit "connect dangling bonds" begin
        for (cc, cH) in enumerate(connect_H)
            new_H = connectDanglingBonds(A, next_combiners, up_combiners, cH, prev_Env.InProgress[:, cc], side, -1, col; kwargs...)
            new_H_mps[length(vert_H_terms) + length(field_H_terms) + 1 + cc] = MPS(Ny, new_H, 0, Ny+1)
        end
    end
    @debug "Connected dangling bonds"

    @timeit "sum H mps" begin
        H_overall    = sum(new_H_mps; cutoff=cutoff, maxdim=env_maxdim)
    end
    @debug "Summed Hs"
    gen_H_terms  = side == :left ? hori_H_terms : side_H_terms
    in_progress  = Matrix{ITensor}(undef, Ny, length(side_H_terms))
    @timeit "gen dangling bonds" begin
        @inbounds for side_term in 1:length(gen_H_terms)
            in_progress[1:Ny, side_term] = generateNextDanglingBonds(A, previous_combiners, next_combiners, up_combiners, gen_H_terms[side_term], prev_Env.I, side, col; kwargs...)
        end
    end
    @debug "Generated next dangling bonds"
    return Environments(new_I, H_overall, in_progress)
end

function buildNewVerticals(A::PEPS, previous_combiners::Vector, next_combiners::Vector{ITensor}, up_combiners::Vector{ITensor}, H, col::Int)::MPO
    Ny, Nx = size(A)
    is_cu               = is_gpu(A)
    col_site_inds       = [findindex(A[row, col], "Site") for row in 1:Ny]
    ops                 = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    vertical_row_a      = H.sites[1][1]
    vertical_row_b      = H.sites[2][1]
    ops[vertical_row_a] = replaceindex!(copy(H.ops[1]), H.site_ind, col_site_inds[vertical_row_a])
    ops[vertical_row_a] = replaceindex!(ops[vertical_row_a], H.site_ind', col_site_inds[vertical_row_a]')
    ops[vertical_row_b] = replaceindex!(copy(H.ops[2]), H.site_ind, col_site_inds[vertical_row_b])
    ops[vertical_row_b] = replaceindex!(ops[vertical_row_b], H.site_ind', col_site_inds[vertical_row_b]')
    internal_cmb_u      = is_cu ? vcat(cuITensor(1.0), up_combiners, cuITensor(1.0)) : vcat(ITensor(1.0), up_combiners, ITensor(1.0))
    AAs            = Vector{ITensor}(undef, Ny)
    AAs[1] = A[1, col] * ops[1] * prime(dag(A[1, col])) * next_combiners[1] * up_combiners[1] * previous_combiners[1]
    @inbounds for row in 2:Ny-1
        AAs[row] = A[row, col] * ops[row] * prime(dag(A[row, col])) * next_combiners[row] * previous_combiners[row] * up_combiners[row-1] * up_combiners[row]
    end
    AAs[Ny] = A[Ny, col] * ops[Ny] * prime(dag(A[Ny, col])) * next_combiners[Ny] * previous_combiners[Ny] * up_combiners[Ny-1]
    return MPO(Ny, AAs, 0, Ny+1)
end

function buildNewFields(A::PEPS, previous_combiners::Vector, next_combiners::Vector{ITensor}, up_combiners::Vector{ITensor}, H, col::Int)::MPO
    Ny, Nx = size(A)
    col_site_inds  = [findindex(A[row, col], "Site") for row in 1:Ny]
    is_cu         = is_gpu(A)
    ops            = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    field_row      = H.sites[1][1]
    ops[field_row] = replaceindex!(copy(H.ops[1]), H.site_ind, col_site_inds[field_row])
    ops[field_row] = replaceindex!(ops[field_row], H.site_ind', col_site_inds[field_row]')
    #internal_cmb_u = is_cu ? vcat(cuITensor(1.0), up_combiners, cuITensor(1.0)) : vcat(ITensor(1.0), up_combiners, ITensor(1.0))
    AAs            = Vector{ITensor}(undef, Ny)
    AAs[1] = A[1, col] * ops[1] * prime(dag(A[1, col])) * next_combiners[1] * up_combiners[1] * previous_combiners[1]
    @inbounds for row in 2:Ny-1
        AAs[row] = A[row, col] * ops[row] * prime(dag(A[row, col])) * next_combiners[row] * previous_combiners[row] * up_combiners[row-1] * up_combiners[row]
    end
    AAs[Ny] = A[Ny, col] * ops[Ny] * prime(dag(A[Ny, col])) * next_combiners[Ny] * previous_combiners[Ny] * up_combiners[Ny-1]
    #AAs            = [ A[row, col] * ops[row] * prime(dag(A[row, col])) * next_combiners[row] * previous_combiners[row] * internal_cmb_u[row] * internal_cmb_u[row+1] for row in 1:Ny ]
    return MPO(Ny, AAs, 0, Ny+1)
end

function buildNewI(A::PEPS, col::Int, previous_combiners::Vector, side::Symbol)::Tuple{MPO, Vector{ITensor}, Vector{ITensor}}
    Ny, Nx         = size(A)
    next_col       = side == :left ? col + 1 : col - 1 # side is handedness of environment
    AA             = [A[row, col] * prime(dag(A[row, col]), "Link") for row in 1:Ny]
    next_combiners = [combine(A[row, col], A[row, next_col], "Site,r$row,c$col") for row in 1:Ny]
    up_combiners   = [combine(A[row, col], A[row+1, col], "Link,CMB,c$col,r$row") for row in 1:Ny-1]
    @inbounds for row in 1:Ny
        AA[row] *= previous_combiners[row]
        AA[row] *= next_combiners[row]
        if row > 1
            AA[row] *= up_combiners[row-1]
        end
        if row < Ny
            AA[row] *= up_combiners[row]
        end
    end
    Iapp   = MPO(Ny, AA, 0, Ny+1)
    return Iapp, next_combiners, up_combiners
end

function generateEdgeDanglingBonds(A::PEPS, up_combiners::Vector{ITensor}, H, side::Symbol, col::Int)::Vector{ITensor}
    Ny, Nx         = size(A)
    is_cu         = is_gpu(A)
    op_row         = side == :left ? H.sites[1][1] : H.sites[2][1];
    H_op           = side == :left ? H.ops[1] : H.ops[2];
    col_site_inds  = [findindex(A[row, col], "Site") for row in 1:Ny]
    ops            = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    ops[op_row]    = replaceindex!(copy(H_op), H.site_ind, col_site_inds[op_row])
    ops[op_row]    = replaceindex!(ops[op_row], H.site_ind', col_site_inds[op_row]')
    #internal_cmb_u = is_cu ? vcat(cuITensor(1.0), up_combiners, cuITensor(1.0)) : vcat(ITensor(1.0), up_combiners, ITensor(1.0))
    #this_IP        = [A[row, col] * ops[row] * prime(dag(A[row, col])) * internal_cmb_u[row] * internal_cmb_u[row+1] for row in 1:Ny]
    this_IP        = Vector{ITensor}(undef, Ny)
    this_IP[1] = A[1, col] * ops[1] * prime(dag(A[1, col])) * up_combiners[1]
    for row in 2:Ny-1
        this_IP[row] = A[row, col] * ops[row] * prime(dag(A[row, col])) * up_combiners[row-1] * up_combiners[row]
    end
    this_IP[Ny] = A[Ny, col] * ops[Ny] * prime(dag(A[Ny, col])) * up_combiners[Ny-1]
    return this_IP
end

function generateNextDanglingBonds(A::PEPS, previous_combiners::Vector{ITensor}, next_combiners::Vector{ITensor}, up_combiners::Vector{ITensor}, H, Ident::MPS, side::Symbol, col::Int; kwargs...)::Vector{ITensor}
    Ny, Nx          = size(A)
    is_cu          = is_gpu(A)
    op_row          = side == :left ? H.sites[1][1] : H.sites[2][1]
    H_op            = side == :left ? H.ops[1] : H.ops[2]
    col_site_inds   = [findindex(A[row, col], "Site") for row in 1:Ny]
    ops             = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    ops[op_row]     = replaceindex!(copy(H_op), H.site_ind, col_site_inds[op_row])
    ops[op_row]     = replaceindex!(ops[op_row], H.site_ind', col_site_inds[op_row]')
    #internal_cmb_u  = is_cu ? vcat(cuITensor(1.0), up_combiners, cuITensor(1.0)) : vcat(ITensor(1.0), up_combiners, ITensor(1.0))
    #this_IP         = [A[row, col] * ops[row] * prime(dag(A[row, col])) * previous_combiners[row] * next_combiners[row] * internal_cmb_u[row] * internal_cmb_u[row+1] for row in 1:Ny]
    this_IP         = Vector{ITensor}(undef, Ny)
    this_IP[1] = A[1, col] * ops[1] * prime(dag(A[1, col])) * previous_combiners[1] * next_combiners[1] * up_combiners[1]
    for row in 2:Ny-1
        this_IP[row] = A[row, col] * ops[row] * prime(dag(A[row, col])) * previous_combiners[row] * next_combiners[row] * up_combiners[row-1] * up_combiners[row]
    end
    this_IP[Ny] = A[Ny, col] * ops[Ny] * prime(dag(A[Ny, col])) * previous_combiners[Ny] * next_combiners[Ny] * up_combiners[Ny-1]
    in_progress_MPO = MPO(Ny, this_IP, 0, Ny+1)
    cutoff::Float64 = get(kwargs, :cutoff, 1e-13)
    maxdim::Int     = get(kwargs, :maxdim, 1)
    env_maxdim::Int = get(kwargs, :env_maxdim, maxdim)
    result          = applyMPO(in_progress_MPO, Ident; cutoff=cutoff, maxdim=env_maxdim)
    return ITensor[result[row]*next_combiners[row] for row in 1:Ny]
end

function connectDanglingBonds(A::PEPS, next_combiners::Vector{ITensor}, up_combiners::Vector{ITensor}, oldH, in_progress::Vector{ITensor}, side::Symbol, work_row::Int, col::Int; kwargs...)::Vector{ITensor}
    Ny, Nx   = size(A)
    is_cu    = is_gpu(A)
    op_row_a = oldH.sites[1][1]
    op_row_b = oldH.sites[2][1]
    op       = side == :left ? oldH.ops[2] : oldH.ops[1]
    application_row = side == :left ? op_row_b : op_row_a
    col_site_inds  = [findindex(A[row, col], "Site") for row in 1:Ny]
    ops            = ITensor[spinI(spin_ind; is_gpu=is_cu) for spin_ind in col_site_inds]
    ops[application_row] = replaceindex!(copy(op), oldH.site_ind, col_site_inds[application_row])
    ops[application_row] = replaceindex!(ops[application_row], oldH.site_ind', col_site_inds[application_row]')
    #internal_cmb_u = is_cu ? vcat(cuITensor(1.0), up_combiners, cuITensor(1.0)) : vcat(ITensor(1.0), up_combiners, ITensor(1.0))
    in_prog_mps    = MPS(Ny, in_progress, 0, Ny + 1)
    this_IP        = Vector{ITensor}(undef, Ny)
    this_IP[1] = A[1, col] * ops[1] * prime(dag(A[1, col])) * next_combiners[1] * up_combiners[1]
    for row in 2:Ny-1
        this_IP[row] = A[row, col] * ops[row] * prime(dag(A[row, col])) * next_combiners[row] * up_combiners[row-1] * up_combiners[row]
    end
    this_IP[Ny] = A[Ny, col] * ops[Ny] * prime(dag(A[Ny, col])) * next_combiners[Ny] * up_combiners[Ny-1]
    #this_IP        = [A[row, col] * ops[row] * prime(dag(A[row, col])) * next_combiners[row] * internal_cmb_u[row] * internal_cmb_u[row+1] for row in 1:Ny]
    if 0 < work_row < Ny + 1
        this_IP[work_row] = ops[work_row]
    end
    completed_H = MPO(Ny, this_IP, 0, Ny+1)
    if work_row == -1
        dummy_cmbs      = [combiner(commoninds(completed_H[row], in_prog_mps[row]), tags="Site,r$row")[1] for row in 1:Ny]
        completed_H.A_  = dummy_cmbs .* tensors(completed_H)
        in_prog_mps.A_  = dummy_cmbs .* tensors(in_prog_mps)
        cutoff::Float64 = get(kwargs, :cutoff, 1e-13)
        maxdim::Int     = get(kwargs, :maxdim, 1)
        env_maxdim::Int = get(kwargs, :env_maxdim, maxdim)
        result          = applyMPO(completed_H, in_prog_mps; cutoff=cutoff, maxdim=env_maxdim)
        return tensors(result)
    else
        @inbounds for row in 1:Ny-1
            ci = linkIndex(completed_H, row)
            ni = Index(ITensors.dim(ci), "u,Link,c$col,r$row")
            replaceIndex!(completed_H[row],   ci, ni)
            replaceIndex!(completed_H[row+1], ci, ni)
        end
        return tensors(completed_H) .* tensors(in_prog_mps)
    end
end

function buildLs(A::PEPS, H; kwargs...)
    Ny, Nx = size(A)
    previous_combiners = Vector{ITensor}(undef, Ny)
    next_combiners     = Vector{ITensor}(undef, Ny)
    Ls                 = Vector{Environments}(undef, Nx)
    start_col::Int     = get(kwargs, :start_col, 1)
    if start_col == 1
        left_H_terms = getDirectional(H[1], Horizontal)
        @debug "Building left col $start_col"
        Ls[1] = buildEdgeEnvironment(A, H, left_H_terms, previous_combiners, :left, 1; kwargs...)
    elseif start_col - 1 > 1
        previous_combiners = [reconnect(commonindex(A[row, start_col], A[row, start_col - 1]), Ls[start_col-1].I[row]) for row in 1:Ny]
    end
    loop_col = start_col == 1 ? 2 : start_col
    @inbounds for col in loop_col:(Nx-1)
        @debug "Building left col $col"
        Ls[col] = buildNextEnvironment(A, Ls[col-1], H, previous_combiners, next_combiners, :left, col; kwargs...)
        previous_combiners = deepcopy(next_combiners)
    end
    return Ls
end

function buildRs(A::PEPS, H; kwargs...)
    Ny, Nx = size(A)
    previous_combiners = Vector{ITensor}(undef, Ny)
    next_combiners     = Vector{ITensor}(undef, Ny)
    start_col::Int     = get(kwargs, :start_col, Nx)
    Rs                 = Vector{Environments}(undef, Nx)
    if start_col == Nx
        right_H_terms = getDirectional(H[Nx-1], Horizontal)
        Rs[Nx]        = buildEdgeEnvironment(A, H, right_H_terms, previous_combiners, :right, Nx; kwargs...)
    elseif start_col + 1 < Nx
        previous_combiners = [reconnect(commonindex(A[row, start_col], A[row, start_col + 1]), Rs[start_col+1].I[row]) for row in 1:Ny]
    end
    loop_col = start_col == Nx ? Nx - 1 : start_col
    @inbounds for col in reverse(2:loop_col)
        @debug "Building right col $col"
        Rs[col] = buildNextEnvironment(A, Rs[col+1], H, previous_combiners, next_combiners, :right, col; kwargs...)
        previous_combiners = deepcopy(next_combiners)
    end
    return Rs
end


text:
using LinearAlgebra
using BenchmarkTools
using Keldysh

suite = BenchmarkGroup()

function setup_indexing_benchmark()
  suite["indexing"] = BenchmarkGroup()

  c = FullContour(tmax=2.0, β=5.0)
  grid = FullTimeGrid(c, 41, 101)
  N = length(grid)

  X = zeros(1,1,N,N)
  Y = Keldysh.GenericStorage(N,N,1,true)

  G1 = GenericTimeGF(grid, 1, true)
  G2 = FullTimeGF(grid, 1, fermionic, true)

  b = @benchmarkable $X[1,1,i,j] setup=(i=rand(1:$N); j=rand(1:$N))
  suite["indexing"]["Array"] = b

  b = @benchmarkable $Y[i,j] setup=(i=rand(1:$N); j=rand(1:$N))
  suite["indexing"]["GenericStorage"] = b

  b = @benchmarkable $G1[t1,t2] setup=(t1=$grid[rand(1:$N)]; t2=$grid[rand(1:$N)])
  suite["indexing"]["GenericTimeGF"] = b

  b = @benchmarkable $G2[t1,t2] setup=(t1=$grid[rand(1:$N)]; t2=$grid[rand(1:$N)])
  suite["indexing"]["FullTimeGF"] = b

  tune!(suite["indexing"])
end

function setup_interp_benchmark()
  suite["interp"] = BenchmarkGroup()

  c = FullContour(tmax=2.0, β=5.0)
  l = length(c)
  grid = FullTimeGrid(c, 41, 101)

  G1 = GenericTimeGF(grid, 1, true)

  b = @benchmarkable $G1(t1,t2) setup=(t1=Keldysh.get_point($c, $l*rand()); t2=Keldysh.get_point($c, $l*rand()))
  suite["interp"]["GenericTimeGF"] = b

  tune!(suite["interp"])
end

function main()
  setup_indexing_benchmark()
  setup_interp_benchmark()

  compare = false
  if length(ARGS) > 0

    if length(ARGS) != 2
      error("usage: julia test/benchmark.jl test/params.ref.json test/results.ref.json")
    end

    params_file = ARGS[1]
    p = BenchmarkTools.load(params_file)[1]
    loadparams!(suite, p, :evals, :samples)

    ref_file = ARGS[2]
    ref = BenchmarkTools.load(ref_file)[1]
    compare = true
  end

  results = run(suite, verbose=true)
  display(median(results))
  println()

  BenchmarkTools.save("params.json", params(suite))
  BenchmarkTools.save("results.json", results)

  if compare
    cmp = judge(median(results), median(ref))
    display(cmp)
  end
end

main()


text:
"""
    load_problem!(lp::AbstractProblem,
                  A::AbstractMatrix{T},
                  b::Vector{T},
                  c::Vector{T},
                  lo::Vector{T},
                  hi::Vector{T}) where T <: Real

Load problem data into a `Problem{T}`.
"""
function load_problem!(
    lp::AbstractProblem{T},
    A::AbstractMatrix{T},
    b::Vector{T},
    c::Vector{T},
    lo::Vector{T},
    hi::Vector{T}
) where T <: Real
    nc = size(A, 1)
    nc == length(b) || throw(DimensionMismatch(
        "constraint matrix has $nc rows, but rhs has only $(length(b))"
    ))

    nv = length(c)
    (nv == length(lo) && nv == length(hi)) || throw(DimensionMismatch(
        "solution vector has expected length $nv"
    ))

    nc <= nv || error("must have fewer constraints than variables")

    # store only the non-zero values of A
    cols = Vector{Column{T}}(undef, size(A, 2))
    for (j, col) in enumerate(eachcol(A))
        cols[j] = Column{T}()
        for (i, v) in enumerate(col)
            if v != 0
                push!(cols[j].ind, i)
                push!(cols[j].nzval, v)
            end
        end
    end

    lp.nc   = nc
    lp.nv   = nv
    lp.cols = cols
    lp.b    = b
    lp.c    = c
    lp.lo   = lo
    lp.hi   = hi

    nothing
end

"""
    solve(lp::AbstractProblem, maxiter::Int = 100, tol::T) where T

Solve the linear program. Returns a `Solution`.
"""
function solve(
    lp::AbstractProblem{T},
    maxiter::Int = 100,
    tol::T = sqrt(eps(T))
) where T <: Real
    # convert problem to standard form
    slp = reformulate(lp)

    # starting point
    iter = Iterate{T}(slp.nc, slp.nv, slp.nu)
    # FIXME: does this need its own func call?
    starting_pt!(iter)
    duality!(iter)

    tols = Tolerances{T}(tol)
    solv = Solver{T}(slp, iter, tols)

    # call internal HSD algorithm
    solve!(solv, maxiter)

    # extract solution from solver
    get_solution(solv, lp)
end


text:
## compose_matrix.jl : methods to compose example matrices

"""
anisotropic2d(ϵ, n, m)

2d anisotropic Laplacian on [0,1]^2 using n x m points and anisotropy number ϵ
"""
anisotropic2d(ϵ, n, m) = stencil2mat(SArray{Tuple{3,3}}([0 -m^2 0 -ϵ*n^2 2(ϵ*n^2+m^2) -ϵ*n^2 0 -m^2 0]), n, m)

"""
anisotropic2d(ϵ, β, n, m)

2d rotated anisotropic Laplacian on [0,1]^2 using n x m points with anisotropy number ϵ and angle β
"""
anisotropic2d(ϵ, β, n, m) = stencil2mat(SArray{Tuple{3,3}}([0 0 0 -n^2*(ϵ*cos(β)^2+sin(β)^2) 2n^2*(ϵ*cos(β)^2+sin(β)^2) -n^2*(ϵ*cos(β)^2+sin(β)^2) 0 0 0]), n, m) + stencil2mat(SArray{Tuple{3,3}}([n*m/2*((ϵ-1)*cos(β)*sin(β)) 0 -n*m/2*((ϵ-1)*cos(β)*sin(β)) 0 0 0 -n*m/2*((ϵ-1)*cos(β)*sin(β)) 0 n*m/2*((ϵ-1)*cos(β)*sin(β))]), n, m) + stencil2mat(SArray{Tuple{3,3}}([0 -m^2*(ϵ*sin(β)^2+cos(β)^2) 0 0 2m^2*(ϵ*sin(β)^2+cos(β)^2) 0 0 -m^2*(ϵ*sin(β)^2+cos(β)^2) 0]), n, m)

"""
anisotropic3d(ϵ₁, ϵ₂, n, m, l)

3d anisotropic Laplacian on [0,1]^3 using n x m x l points and anisotropy numbers ϵ₁ and ϵ₂
"""
anisotropic3d(ϵ₁, ϵ₂, n, m, l) = stencil2mat(SArray{Tuple{3,3,3}}([0 0 0 0 -l^2 0 0 0 0 0 -ϵ₂*m^2 0 -ϵ₁*n^2 2(ϵ₁*n^2+ϵ₂*m^2+l^2) -ϵ₁*n^2 0 -ϵ₂*m^2 0 0 0 0 0 -l^2 0 0 0 0]), n, m, l)


text:
include("../imports.jl")

@testset ExtendedTestSet "StopEarly" begin
    learner = testlearner(Recorder(), Metrics(), EarlyStopping(1), ProgressPrinter(), coeff = 3, opt = Descent(0.1))
    @test_throws CancelFittingException begin
        @suppress fit!(learner, 100)
    end
    @test_nowarn print(Base.DevNull(), EarlyStopping(1))
    @test_nowarn show(Base.DevNull(), EarlyStopping(FluxTraining.ES.NumberLimit(2)))
end


text:
module MaskArrays

using ConcreteStructs

EqualTo = Base.Fix2{typeof(isequal)}
@concrete terse struct MaskArray{T,N} <: AbstractArray{T,N}
    base
    indices
    imputed
    buffer
end

export imputed
imputed(ma::MaskArray) = ma.imputed

export buffer
buffer(ma::MaskArray) = ma.buffer

Base.length(x::MaskArray) = Base.length(x.base)
Base.size(x::MaskArray) = Base.size(x.base)

Base.getindex(ma::MaskArray, inds...) = getindex(ma.base, inds...)

export maskarray

function maskarray(x::AbstractArray{Union{T,Missing},N}) where {T,N}
    inds = CartesianIndices(x)

    # Get the locations of missing values
    missinginds = inds[findall(ismissing, x)]
    nummissing = length(missinginds)

    # # Build base array, fill with known data
    base = similar(x, T)
    for j in inds
        j ∈ missinginds && continue
        base[j] = x[j]
    end

    indices = Dict(zip(missinginds, 1:nummissing))
    vu = view(base, missinginds)

    return MaskArray{T,N}(base, indices, vu, vu)
end

export replace_buffer

function replace_buffer(ma::MaskArray{T,N}, x::AbstractVector) where {T,N}
    @assert eltype(x) == eltype(ma.base)
    @assert length(x) == length(ma.imputed)
    MaskArray{T,N}(ma.base, ma.indices, ma.imputed, x)
end

export sync!

function sync!(ma::MaskArray)
    ma.imputed .= ma.buffer
    return ma
end

end


text:
using BrainWave, Test

x = [1 2 3 4;
     5 6 7 8;
     9 10 11 12;
     13 14 15 16
     ]

@test isequal(deleteSlice2D(x, 1, 1), x[2:end,:])
@test isequal(deleteSlice2D(x, [1,4], 1), x[2:3,:])
@test isequal(deleteSlice2D(x, [4,1], 1), x[2:3,:])

@test isequal(deleteSlice2D(x, 1, 2), x[:,2:end])
@test isequal(deleteSlice2D(x, [1,4], 2), x[:,2:3])
@test isequal(deleteSlice2D(x, [4,1], 2), x[:,2:3])


text:
function print_results(fileflag, seqpa, seqpo, J, H, L, ctype, lambda_o, lambda_e, μext, μint, seqtest, seqins, init, fin, start, fname, satmf, beta, T0, el_time)

	enmf = compute_potts_en(J, H, seqpa, L,ctype)
	enhmm = compute_potts_en(J,H, seqtest, L ,ctype)
	costmf = compute_cost_function(J,H,seqpo,L,ctype,lambda_o,lambda_e,μext,μint)
	costhmm = compute_cost_function(J,H,seqins,L,ctype,lambda_o,lambda_e,μext,μint)
	distmf, mgmf, lgmf, dmmf = hammingdist(seqpa, seqtest)

	@printf(fileflag, ">%s/%d-%d\t", fname, init + parse(Int64, start) - 1, fin + parse(Int64, start) - 1)
	@printf(fileflag, "sat: %s beta: %.2f muext: %.1f muint: %.1f T0: %s Dist: %d %d %d %d Time: %.2f\n", satmf, beta, μext, μint, T0, distmf,mgmf,lgmf,dmmf,el_time)
	@printf(fileflag, "test     %s %.3f %.3f\n", seqtest, enhmm, costhmm)
	@printf(fileflag, "dcalign  %s %.3f %.3f\n", seqpa, enmf, costmf)
	@printf(fileflag, "ins test     %s\n", seqins)
	@printf(fileflag, "ins dcalign  %s\n", seqpo)
	flush(fileflag)
end

function align_all(
    q::Int,
    L::Int,
    filename_par::String,
    filename_full::String,
    inspen_file::String,
    μext::Real,
    μint::Real;
    delta::Int = 20,
    typel::Symbol = :bm,
    filename_align::String = "",
    filename_ins::String = "",
    β::Real = 1.0,
    filename_out::String = "out.fasta",
    filename_flag::String = "flag.dat",
    seedrng::Int = 0,
    mindec::Int = 50,
    nprint::Int = 100,
    maxiter::Int = 1000,
    epsconv::Real = 1e-5,
    verbose::Bool = false,
)

    if q == 5
        #ctype = Symbol("nbase")
        ctype = :nbase
        delta = 1000
    else
        #ctype = Symbol("amino")
        ctype = :amino
    end

    if typel == :bm
        gap = 0
    else
        gap = q
    end
    println("### Reading parameters from file ###")
    J, H = read_parameters(filename_par, q, L, gap = gap, typel = typel)
    file_out = open(filename_out, "a")
    println("### Reading sequences ###")
    typedic = true
    if filename_align == "" || filename_ins == ""
        println("No delta")
        println("File full: ", filename_full)
        al = enveloptoalign(filename_full, ctype = ctype)
        names = collect(keys(al))
        fake = '-'^L
    else
        println("### Cutting the full length sequence using delta =", delta)
        file_flag = open(filename_flag, "a")
        typedic = false
        println("File full: ", filename_full)
        println("File aligned seqs: ", filename_align)
        println("File ins seqs: ", filename_ins)
        al = enveloptoalign(
            filename_full,
            filename_align,
            filename_ins,
            delta = delta,
            ctype = ctype,
        )
    end

    M = length(al)

    Lambda_all = readdlm(inspen_file)

    lambda_o = Lambda_all[:, 1]
    lambda_e = Lambda_all[:, 2]
    damp = 0.0
    println("Tot number of sequences to be aligned ", M)

    tenp = ceil(Int64, M / 10)
    println("### Starting loop ###")
    for k = 1:M
        if !typedic
            (name, hit) = split(al[k][1], "/")
            (start, finish) = split(hit, "-")
            seqtest = al[k][3]
            seqins = al[k][4]
            fullseq = al[k][2]
            seq = Seq(seqtest, fullseq, ctype)
        else
            fname = names[k]
            fullseq = al[fname]
            if occursin("/", fname)
                (name, hit) = split(fname, "/")
                (start, finish) = split(hit, "-")
            else
                name = fname
                start = "1"
                finish = @sprintf("%s", L)
            end
            seq = Seq(fake, fullseq, ctype)
        end
        if parse(Int64, start) > parse(Int64, finish) # it often happens in RNA seq
            aux = finish
            finish = start
            start = aux
        end

        N = length(fullseq)
        T0 = true
        stime = time()
        res = palign(
            seq,
            β .* J,
            β .* H,
            β .* lambda_o,
            β .* lambda_e,
            ctype = ctype,
            mindec = mindec,
            nprint = nprint,
            maxiter = maxiter,
            T0 = T0,
            μint = β * μint,
            seed = seedrng,
            damp = damp,
            μext = β * μext,
            epsconv = epsconv,
            verbose = verbose,
        )
        el_time = time() -stime
        P = res[3].pbf.P
        out = decodeposterior(P, seq.strseq)
        seqpa0 = out.pa
        seqpo0 = out.po
        sat = check_assignment(P, false, N)
        if sat == false
            (seqpa0, P) = decimate_post(res[3], T0)
            out = decodeposterior(P, seq.strseq)
            seqpo0 = out.po
        end
        b0 = out.start
        e0 = out.finish
        enT0 = compute_potts_en(J, H, seqpa0, L, ctype)
        costT0 = compute_cost_function(J,H,seqpo0,L,ctype,lambda_o,lambda_e,μext,μint)

        if !typedic
            print_results(
                file_flag,
                seqpa0,
                seqpo0,
                β .* J,
                β .* H,
                L,
                ctype,
                β * lambda_o,
                β * lambda_e,
                β * μext,
                β * μint,
                seqtest,
                seqins,
                b0,
                e0,
                start,
                name,
                sat,
                β,
                T0,
                el_time,
            )
        end

        T0 = false
        stime = time()
        res = palign(
            seq,
            β .* J,
            β .* H,
            β .* lambda_o,
            β .* lambda_e,
            ctype = ctype,
            mindec = mindec,
            nprint = nprint,
            maxiter = maxiter,
            T0 = T0,
            μint = β * μint,
            seed = seedrng,
            damp = damp,
            μext = β * μext,
            epsconv = epsconv,
            verbose = verbose,
        )
        el_time = time() - stime
        P = res[3].pbf.P
        out = decodeposterior(P, seq.strseq)
        seqpa = out.pa
        seqpo = out.po
        sat = check_assignment(P, false, N)
        if sat == false
            (seqpa, P) = decimate_post(res[3], T0)
            out = decodeposterior(P, seq.strseq)
            seqpo = out.po
        end
        enT = compute_potts_en(J, H, seqpa, L, ctype)
        costT = compute_cost_function(J,H,seqpo,L,ctype,lambda_o,lambda_e,μext,μint)
        b = out.start
        e = out.finish

        #if enT < enT0
        if costT < costT0
            @printf(
                file_out,
                ">%s/%i-%i\n",
                name,
                b + parse(Int64, start) - 1,
                e + parse(Int64, start) - 1
            )
            @printf(file_out, "%s\n", seqpa)
        else
            @printf(
                file_out,
                ">%s/%i-%i\n",
                name,
                b0 + parse(Int64, start) - 1,
                e0 + parse(Int64, start) - 1
            )
            @printf(file_out, "%s\n", seqpa0)
        end
        flush(file_out)

        if !typedic
            print_results(
                file_flag,
                seqpa,
                seqpo,
                β .* J,
                β .* H,
                L,
                ctype,
                β * lambda_o,
                β * lambda_e,
                β * μext,
                β * μint,
                seqtest,
                seqins,
                b,
                e,
                start,
                name,
                sat,
                β,
                T0,
                el_time,
            )
        end

        if mod(k, tenp) == 0
            println("Done...", k, "/", M)
        end

    end
    println("### End ###")
end


text:
using Test
using AtomicGraphNets
using ChemistryFeaturization
using SimpleWeightedGraphs

@testset "CGCNN" begin
    # build a model with deterministic initialization of weights
    in_fea_len = 40
    conv_fea_len = 20
    pool_fea_len = 10
    model = Xie_model(in_fea_len, atom_conv_feature_length=conv_fea_len, pooled_feature_length=pool_fea_len, num_hidden_layers=2, initW=ones)

    # check that everything is the size it should be
    @test length(model)==5
    @test size(model[1].convweight) == size(model[1].selfweight) == (conv_fea_len, in_fea_len)
    @test size(model[2].convweight) == size(model[2].selfweight) == (conv_fea_len, conv_fea_len)
    @test size(model[4].W) == (pool_fea_len, pool_fea_len)
    @test size(model[5].W) == (1, pool_fea_len)

    # check that it evaluates to the right things, basically
    dummyfzn = [AtomFeat(:feat, [0]) for i in 1:40]
    input = AtomGraph(SimpleWeightedGraph{Int32,Float32}([0 1; 1 0]), ["H", "O"], ones(Float32, 40, 2), dummyfzn)
    output1 = model[1](input)
    #println(output1.features)
    int_mat = model[2].convweight * output1.features * output1.lapl + model[2].selfweight * output1.features + hcat([model[2].bias for i in 1:size(output1.features, 2)]...)
    #println(int_mat)
    #println(model[2].σ.(int_mat))
    #println(AtomicGraphNets.reg_norm(model[2].σ.(int_mat)))
    # TODO: figure out why the reg_norm step gives different results in REPL than in testing
    @test all(isapprox.(model[1:2](input).features, zeros(Float32, 20,2), atol=2e-3))
    @test isapprox(model(input)[1], 6.9, atol=3e-2)
end


text:
using ITensors

examples_dir = joinpath(dirname(pathof(ITensors)),
                        "..", "examples", "src")
include(joinpath(examples_dir, "electronk.jl"))
include(joinpath(examples_dir, "hubbard.jl"))

function run(; maxdim::Int,
               nsweeps::Int = 10,
               outputlevel::Int = 0,
               conserve_ky::Bool = true,
               cutoff::Float64 = 0.0,
               Nx::Int = 6,
               Ny::Int = 3,
               U::Float64 = 4.0,
               t::Float64 = 1.0)
  N = Nx * Ny
  sweeps = Sweeps(nsweeps)
  maxdims = min.(maxdim, [100, 200, 400, 800, 2000, 3000, maxdim])
  maxdim!(sweeps, maxdims...)
  cutoff!(sweeps, cutoff)
  noise!(sweeps, 1e-6, 1e-7, 1e-8, 0.0)
  sites = siteinds("ElecK", N;
                   conserve_qns = true,
                   conserve_ky = conserve_ky,
                   modulus_ky = Ny)
  ampo = hubbard(Nx = Nx, Ny = Ny, t = t, U = U, ky = true)
  H = MPO(ampo, sites)
  # Create start state
  state = Vector{String}(undef, N)
  for i in 1:N
    x = (i - 1) ÷ Ny
    y = (i - 1) % Ny
    if x % 2 == 0
      if y % 2 == 0
        state[i] = "Up"
      else
        state[i] = "Dn"
      end
    else
      if y % 2 == 0
        state[i] = "Dn"
      else
        state[i] = "Up"
      end
    end
  end
  ψ0 = productMPS(sites, state)
  energy, ψ = dmrg(H, ψ0, sweeps; outputlevel = outputlevel)
  if outputlevel > 0
    @show nsweeps
    @show cutoff
    @show energy
    @show flux(ψ)
  end
  return maxlinkdim(ψ)
end


text:
import Base.+, Base.-

export +, -

# Unary Negation

function -(x::Constant)
  return Constant(-x.value)
end

function -(x::AffineExpr)
  vars_to_coeffs_map = Dict{UInt64, Constant}()
  for (v, c) in x.vars_to_coeffs_map
    vars_to_coeffs_map[v] = -c
  end
  this = AffineExpr(:-, (x,), vars_to_coeffs_map, -x.constant, x.size)
  this.evaluate = ()->(-x.evaluate())
  return this
end

# Binary Addition

function +(x::Constant, y::Constant)
  if x.size == (1, 1) || y.size == (1, 1) || x.size == y.size
    return convert(Constant, x.value + y.value)
  else
    error("Cannot add two expressions with sizes $(x.size) and $(y.size)")
  end
end

function +(x::AffineExpr, y::Constant)
  if x.size != (1, 1) && y.size != (1, 1) && x.size != y.size
    error("Cannot add two expressions with sizes $(x.size) and $(y.size)")
  end

  if x.size == (1, 1) && y.size != (1, 1)
    x = repmat(x, y.size[1], y.size[2])
  end

  vars_to_coeffs_map = copy(x.vars_to_coeffs_map)
  constant = x.constant + vec(y)
  this = AffineExpr(:+, (x, y), vars_to_coeffs_map, constant, x.size)
  this.evaluate = ()->x.evaluate() + y.evaluate()
  return this
end

function +(x::Constant, y::AffineExpr)
  return y + x
end

function +(x::AffineExpr, y::AffineExpr)
  if x.size != (1, 1) && y.size != (1, 1) && x.size != y.size
    error("Cannot add two expressions with sizes $(x.size) and $(y.size)")
  end

  if x.size == (1, 1) && y.size != (1, 1)
    x = repmat(x, y.size[1], y.size[2])
  elseif y.size == (1, 1) && x.size != (1, 1)
    y = repmat(y, x.size[1], x.size[2])
  end

  vars_to_coeffs_map = copy(x.vars_to_coeffs_map)
  for (v, c) in y.vars_to_coeffs_map
    if haskey(vars_to_coeffs_map, v)
      vars_to_coeffs_map[v] = vars_to_coeffs_map[v] + c
    else
      vars_to_coeffs_map[v] = c
    end
  end
  constant = x.constant + y.constant
  this = AffineExpr(:+, (x, y), vars_to_coeffs_map, constant, x.size)
  this.evaluate = ()->x.evaluate() + y.evaluate()
  return this
end

+(x::AffineExpr, y::Numeric) = +(x, convert(Constant, y))
+(x::Numeric, y::AffineExpr) = +(y, convert(Constant, x))

function +(x::SumSquaresExpr, y::SumSquaresExpr)
  affines = copy(x.affines)
  multipliers = copy(x.multipliers)
  append!(affines, y.affines)
  append!(multipliers, y.multipliers)
  this = SumSquaresExpr(:+, affines, multipliers, x.scalar + y.scalar)
  return this
end

function +(x::Number, y::SumSquaresExpr)
  try
    x = convert(Float64, x)
  catch
    error("Only real scalars can be added to Sum Squares expressions")
  end
  if x < 0
    error("Only nonnegative scalars can be added to Sum Squares expressions")
  end
  this = SumSquaresExpr(:+, y.affines, y.multipliers, y.scalar + x)
  return this
end

+(x::SumSquaresExpr, y::Number) = +(y, x)

# Binary Subtraction
-(x::AffineExpr, y::AffineExpr) = +(x, -y)
-(x::AffineExpr, y::Constant) = +(x, -y)
-(x::Constant, y::AffineExpr) = +(-y, x)
-(x::AffineExpr, y::Numeric) = -(x, convert(Constant, y))
-(x::Numeric, y::AffineExpr) = -(convert(Constant, x), y)


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""Endpoints is a collection of endpoints that implement the actual service. Example:   Name: \"mysvc\",   Subsets: [     {       Addresses: [{\"ip\": \"10.10.1.1\"}, {\"ip\": \"10.10.2.2\"}],       Ports: [{\"name\": \"a\", \"port\": 8675}, {\"name\": \"b\", \"port\": 309}]     },     {       Addresses: [{\"ip\": \"10.10.3.3\"}],       Ports: [{\"name\": \"a\", \"port\": 93}, {\"name\": \"b\", \"port\": 76}]     },  ]

    IoK8sApiCoreV1Endpoints(;
        apiVersion=nothing,
        kind=nothing,
        metadata=nothing,
        subsets=nothing,
    )

    - apiVersion::String : APIVersion defines the versioned schema of this representation of an object. Servers should convert recognized schemas to the latest internal value, and may reject unrecognized values. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#resources
    - kind::String : Kind is a string value representing the REST resource this object represents. Servers may infer this from the endpoint the client submits requests to. Cannot be updated. In CamelCase. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#types-kinds
    - metadata::IoK8sApimachineryPkgApisMetaV1ObjectMeta : Standard object's metadata. More info: https://git.k8s.io/community/contributors/devel/sig-architecture/api-conventions.md#metadata
    - subsets::Vector{IoK8sApiCoreV1EndpointSubset} : The set of all endpoints is the union of all subsets. Addresses are placed into subsets according to the IPs they share. A single address with multiple ports, some of which are ready and some of which are not (because they come from different containers) will result in the address being displayed in different subsets for the different ports. No address will appear in both Addresses and NotReadyAddresses in the same subset. Sets of addresses and ports that comprise a service.
"""
mutable struct IoK8sApiCoreV1Endpoints <: SwaggerModel
    apiVersion::Any # spec type: Union{ Nothing, String } # spec name: apiVersion
    kind::Any # spec type: Union{ Nothing, String } # spec name: kind
    metadata::Any # spec type: Union{ Nothing, IoK8sApimachineryPkgApisMetaV1ObjectMeta } # spec name: metadata
    subsets::Any # spec type: Union{ Nothing, Vector{IoK8sApiCoreV1EndpointSubset} } # spec name: subsets

    function IoK8sApiCoreV1Endpoints(;apiVersion=nothing, kind=nothing, metadata=nothing, subsets=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1Endpoints, Symbol("apiVersion"), apiVersion)
        setfield!(o, Symbol("apiVersion"), apiVersion)
        validate_property(IoK8sApiCoreV1Endpoints, Symbol("kind"), kind)
        setfield!(o, Symbol("kind"), kind)
        validate_property(IoK8sApiCoreV1Endpoints, Symbol("metadata"), metadata)
        setfield!(o, Symbol("metadata"), metadata)
        validate_property(IoK8sApiCoreV1Endpoints, Symbol("subsets"), subsets)
        setfield!(o, Symbol("subsets"), subsets)
        o
    end
end # type IoK8sApiCoreV1Endpoints

const _property_map_IoK8sApiCoreV1Endpoints = Dict{Symbol,Symbol}(Symbol("apiVersion")=>Symbol("apiVersion"), Symbol("kind")=>Symbol("kind"), Symbol("metadata")=>Symbol("metadata"), Symbol("subsets")=>Symbol("subsets"))
const _property_types_IoK8sApiCoreV1Endpoints = Dict{Symbol,String}(Symbol("apiVersion")=>"String", Symbol("kind")=>"String", Symbol("metadata")=>"IoK8sApimachineryPkgApisMetaV1ObjectMeta", Symbol("subsets")=>"Vector{IoK8sApiCoreV1EndpointSubset}")
Base.propertynames(::Type{ IoK8sApiCoreV1Endpoints }) = collect(keys(_property_map_IoK8sApiCoreV1Endpoints))
Swagger.property_type(::Type{ IoK8sApiCoreV1Endpoints }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1Endpoints[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1Endpoints }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1Endpoints[property_name]

function check_required(o::IoK8sApiCoreV1Endpoints)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1Endpoints }, name::Symbol, val)
end


text:
using Test
using Latexify
using ModelingToolkit

### Tips for generating latex tests:
### Latexify has an unexported macro:
###
### Latexify.@generate_test
###
### which generates a test using a given latexify function.
### For example:
###
### Latexify.@generate_test latexify([1, 2, 3], [4, 5, 6]; env=:mdtable)
###
### This puts a ready-made test in your clipboard which you can paste into the
### test file.
###
### Just be sure to remove all such macros before you commit a change since it
### will cause issues with Travis.

@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x)*D(x-y)/D(z),
       0 ~ σ*x*(ρ-z)/10-y,
       D(z) ~ x*y^(2//3) - β*z]

@test latexify(eqs) ==
raw"\begin{align}
\frac{dx(t)}{dt} =& \frac{\sigma \left( \mathrm{y}\left( t \right) - \mathrm{x}\left( t \right) \right) \frac{d\left(\mathrm{x}\left( t \right) - \mathrm{y}\left( t \right)\right)}{dt}}{\frac{dz(t)}{dt}} \\
0 =& \frac{\sigma \mathrm{x}\left( t \right) \left( \rho - \mathrm{z}\left( t \right) \right)}{10} - \mathrm{y}\left( t \right) \\
\frac{dz(t)}{dt} =& \mathrm{x}\left( t \right) \left( \mathrm{y}\left( t \right) \right)^{\frac{2}{3}} - \beta \mathrm{z}\left( t \right)
\end{align}
"

@parameters t p[1:3]
@variables u[1:3](t)
@derivatives D'~t

eqs = [D(u[1]) ~ p[3]*(u[2]-u[1]),
       0 ~ p[2]*p[3]*u[1]*(p[1]-u[1])/10-u[2],
       D(u[3]) ~ u[1]*u[2]^(2//3) - p[3]*u[3]]

@test latexify(eqs) ==
raw"\begin{align}
\frac{du_1(t)}{dt} =& p_3 \left( \mathrm{u_2}\left( t \right) - \mathrm{u_1}\left( t \right) \right) \\
0 =& \frac{p_2 p_3 \mathrm{u_1}\left( t \right) \left( p_1 - \mathrm{u_1}\left( t \right) \right)}{10} - \mathrm{u_2}\left( t \right) \\
\frac{du_3(t)}{dt} =& \mathrm{u_1}\left( t \right) \left( \mathrm{u_2}\left( t \right) \right)^{\frac{2}{3}} - p_3 \mathrm{u_3}\left( t \right)
\end{align}
"

eqs = [D(u[1]) ~ p[3]*(u[2]-u[1]),
       D(u[2]) ~ p[2]*p[3]*u[1]*(p[1]-u[1])/10-u[2],
       D(u[3]) ~ u[1]*u[2]^(2//3) - p[3]*u[3]]

sys = ODESystem(eqs)

@test latexify(sys.eqs) ==
raw"\begin{align}
\frac{du_1(t)}{dt} =& p_3 \left( \mathrm{u_2}\left( t \right) - \mathrm{u_1}\left( t \right) \right) \\
\frac{du_2(t)}{dt} =& \frac{p_2 p_3 \mathrm{u_1}\left( t \right) \left( p_1 - \mathrm{u_1}\left( t \right) \right)}{10} - \mathrm{u_2}\left( t \right) \\
\frac{du_3(t)}{dt} =& \mathrm{u_1}\left( t \right) \left( \mathrm{u_2}\left( t \right) \right)^{\frac{2}{3}} - p_3 \mathrm{u_3}\left( t \right)
\end{align}
"

@test latexify(sys) ==
raw"$\begin{align}
\frac{du_1(t)}{dt} =& p_3 \left( \mathrm{u_2}\left( t \right) - \mathrm{u_1}\left( t \right) \right) \\
\frac{du_2(t)}{dt} =& \frac{p_2 p_3 \mathrm{u_1}\left( t \right) \left( p_1 - \mathrm{u_1}\left( t \right) \right)}{10} - \mathrm{u_2}\left( t \right) \\
\frac{du_3(t)}{dt} =& \mathrm{u_1}\left( t \right) \left( \mathrm{u_2}\left( t \right) \right)^{\frac{2}{3}} - p_3 \mathrm{u_3}\left( t \right)
\end{align}
$"


text:
##################
# BenchmarkGroup #
##################

immutable BenchmarkGroup
    tags::Vector{Any}
    data::Dict{Any,Any}
end

BenchmarkGroup(tags::Vector, args::Pair...) = BenchmarkGroup(tags, Dict(args...))
BenchmarkGroup(args::Pair...) = BenchmarkGroup([], args...)

function addgroup!(suite::BenchmarkGroup, id, args...)
    g = BenchmarkGroup(args...)
    suite[id] = g
    return g
end

# Dict-like methods #
#-------------------#

@compat Base.:(==)(a::BenchmarkGroup, b::BenchmarkGroup) = a.tags == b.tags && a.data == b.data
Base.copy(group::BenchmarkGroup) = BenchmarkGroup(copy(group.tags), copy(group.data))
Base.similar(group::BenchmarkGroup) = BenchmarkGroup(copy(group.tags), similar(group.data))
Base.isempty(group::BenchmarkGroup) = isempty(group.data)
Base.length(group::BenchmarkGroup) = length(group.data)
Base.getindex(group::BenchmarkGroup, i...) = getindex(group.data, i...)
Base.setindex!(group::BenchmarkGroup, i...) = setindex!(group.data, i...)
Base.delete!(group::BenchmarkGroup, k...) = delete!(group.data, k...)
Base.haskey(group::BenchmarkGroup, k) = haskey(group.data, k)
Base.keys(group::BenchmarkGroup) = keys(group.data)
Base.values(group::BenchmarkGroup) = values(group.data)
Base.start(group::BenchmarkGroup) = start(group.data)
Base.next(group::BenchmarkGroup, state) = next(group.data, state)
Base.done(group::BenchmarkGroup, state) = done(group.data, state)

# mapping/filtering #
#-------------------#

andexpr(a, b) = :($a && $b)
andreduce(preds) = reduce(andexpr, preds)

@generated function mapvals!(f, dest::BenchmarkGroup, srcs::BenchmarkGroup...)
    haskeys = andreduce([:(haskey(srcs[$i], k)) for i in 1:length(srcs)])
    getinds = [:(srcs[$i][k]) for i in 1:length(srcs)]
    return quote
        for k in keys(first(srcs))
            if $(haskeys)
                dest[k] = f($(getinds...))
            end
        end
        return dest
    end
end

mapvals!(f, group::BenchmarkGroup) = mapvals!(f, similar(group), group)
mapvals(f, groups::BenchmarkGroup...) = mapvals!(f, similar(first(groups)), groups...)

filtervals!(f, group::BenchmarkGroup) = (filter!((k, v) -> f(v), group.data); return group)
filtervals(f, group::BenchmarkGroup) = filtervals!(f, copy(group))

Base.filter!(f, group::BenchmarkGroup) = (filter!(f, group.data); return group)
Base.filter(f, group::BenchmarkGroup) = filter!(f, copy(group))

# benchmark-related methods #
#---------------------------#

Base.minimum(group::BenchmarkGroup) = mapvals(minimum, group)
Base.maximum(group::BenchmarkGroup) = mapvals(maximum, group)
Base.mean(group::BenchmarkGroup) = mapvals(mean, group)
Base.median(group::BenchmarkGroup) = mapvals(median, group)
Base.min(groups::BenchmarkGroup...) = mapvals(min, groups...)
Base.max(groups::BenchmarkGroup...) = mapvals(max, groups...)

Base.time(group::BenchmarkGroup) = mapvals(time, group)
gctime(group::BenchmarkGroup) = mapvals(gctime, group)
memory(group::BenchmarkGroup) = mapvals(memory, group)
allocs(group::BenchmarkGroup) = mapvals(allocs, group)
params(group::BenchmarkGroup) = mapvals(params, group)

ratio(groups::BenchmarkGroup...) = mapvals(ratio, groups...)
judge(groups::BenchmarkGroup...; kwargs...) = mapvals((x...) -> judge(x...; kwargs...), groups...)

rmskew!(group::BenchmarkGroup) = mapvals!(rmskew!, group)
rmskew(group::BenchmarkGroup) = mapvals(rmskew, group)

isregression(group::BenchmarkGroup) = any(isregression, values(group))
isimprovement(group::BenchmarkGroup) = any(isimprovement, values(group))
isinvariant(group::BenchmarkGroup) = all(isinvariant, values(group))

invariants(x) = x
regressions(x) = x
improvements(x) = x
invariants(group::BenchmarkGroup) = mapvals!(invariants, filtervals(isinvariant, group))
regressions(group::BenchmarkGroup) = mapvals!(regressions, filtervals(isregression, group))
improvements(group::BenchmarkGroup) = mapvals!(improvements, filtervals(isimprovement, group))

function loadparams!(group::BenchmarkGroup, paramsgroup::BenchmarkGroup, fields...)
    for (k, v) in group
        haskey(paramsgroup, k) && loadparams!(v, paramsgroup[k], fields...)
    end
    return group
end

# leaf iteration/indexing #
#-------------------------#

leaves(group::BenchmarkGroup) = leaves!(Any[], Any[], group)

function leaves!(results, parents, group::BenchmarkGroup)
    for (k, v) in group
        keys = vcat(parents, k)
        if isa(v, BenchmarkGroup)
            leaves!(results, keys, v)
        else
            push!(results, (keys, v))
        end
    end
    return results
end

function Base.getindex(group::BenchmarkGroup, keys::Vector)
    k = first(keys)
    v = length(keys) == 1 ? group[k] : group[k][keys[2:end]]
    return v
end

function Base.setindex!(group::BenchmarkGroup, x, keys::Vector)
    k = first(keys)
    if length(keys) == 1
        group[k] = x
        return x
    else
        return setindex!(group[k], x, keys[2:end])
    end
end

# tagging #
#---------#

immutable TagFilter{P}
    predicate::P
end

macro tagged(expr)
    return esc(:(BenchmarkTools.TagFilter(tags -> $(tagpredicate!(expr)))))
end

tagpredicate!(tag) = :(in($tag, tags))

function tagpredicate!(sym::Symbol)
    sym == :! && return sym
    sym == :ALL && return true
    return :(in($sym, tags))
end

# build the body of the tag predicate in place
function tagpredicate!(expr::Expr)
    expr.head == :quote && return :(in($expr, tags))
    for i in eachindex(expr.args)
        expr.args[i] = tagpredicate!(expr.args[i])
    end
    return expr
end

function Base.getindex(src::BenchmarkGroup, f::TagFilter)
    dest = similar(src)
    loadtagged!(f, dest, src, src, Any[], src.tags)
    return dest
end

# normal union doesn't have the behavior we want
# (e.g. union(["1"], "2") === ["1", '2'])
vcatunion(args...) = unique(vcat(args...))

function loadtagged!(f::TagFilter, dest::BenchmarkGroup, src::BenchmarkGroup,
                     group::BenchmarkGroup, keys::Vector, tags::Vector)
    ismatch = f.predicate(tags)
    if ismatch
        child_dest = createchild!(dest, src, keys)
        for (k, v) in group
            if isa(v, BenchmarkGroup)
                loadtagged!(f, dest, src, v, vcatunion(keys, k), vcatunion(tags, k, v.tags))
            else
                child_dest[k] = v
            end
        end
    else
        for (k, v) in group
            if isa(v, BenchmarkGroup)
                loadtagged!(f, dest, src, v, vcatunion(keys, k), vcatunion(tags, k, v.tags))
            end
        end
    end
    return dest
end

function createchild!(dest, src, keys)
    if isempty(keys)
        return dest
    else
        k = first(keys)
        src_child = src[k]
        if !(haskey(dest, k))
            dest_child = similar(src_child)
            dest[k] = dest_child
        else
            dest_child = dest[k]
        end
        return createchild!(dest_child, src_child, keys[2:end])
    end
end

# indexing by BenchmarkGroup #
#----------------------------#

function Base.getindex(group::BenchmarkGroup, x::BenchmarkGroup)
    result = BenchmarkGroup()
    for (k, v) in x
        result[k] = isa(v, BenchmarkGroup) ? group[k][v] : group[k]
    end
    return result
end

Base.setindex!(group::BenchmarkGroup, v, k::BenchmarkGroup) = error("A BenchmarkGroup cannot be a key in a BenchmarkGroup")

# pretty printing #
#-----------------#

tagrepr(tags) = string("[", join(map(repr, tags), ", "), "]")

Base.showall(io::IO, group::BenchmarkGroup) = show(io, group; verbose = true, limit = Inf)

# written this way for v0.5/v0.4 compatibility
_showcompact(io::IO, group::BenchmarkGroup) = print(io, "$(length(group))-element BenchmarkGroup($(tagrepr(group.tags)))")

if VERSION < v"0.5.0-dev+4305"
    Base.showcompact(io::IO, group::BenchmarkGroup) = _showcompact(io, group)
    function Base.show(io::IO, group::BenchmarkGroup, pad = ""; verbose = false, limit = 10)
        println(io, "$(length(group))-element BenchmarkTools.BenchmarkGroup:")
        print(io, pad, "  tags: ", tagrepr(group.tags))
        count = 1
        for (k, v) in group
            println(io)
            print(io, pad, "  ", repr(k), " => ")
            if verbose
                if isa(v, BenchmarkGroup)
                    show(io, v, "\t"*pad; verbose = verbose, limit = limit)
                else
                    show(io, v)
                end
            else
                showcompact(io, v)
            end
            count > limit && (println(io); print(io, pad, "  ⋮"); break)
            count += 1
        end
    end
else
    function Base.show(io::IO, group::BenchmarkGroup, pad = ""; verbose = false, limit = 10)
        if get(io, :multiline, true)
            println(io, "$(length(group))-element BenchmarkTools.BenchmarkGroup:")
            print(io, pad, "  tags: ", tagrepr(group.tags))
            count = 1
            element_io = verbose ? io : IOContext(io, :multiline => false)
            for (k, v) in group
                println(io)
                print(io, pad, "  ", repr(k), " => ")
                if verbose && isa(v, BenchmarkGroup)
                    show(element_io, v, "\t"*pad; verbose = verbose, limit = limit)
                else
                    show(element_io, v)
                end
                count > limit && (println(io); print(io, pad, "  ⋮"); break)
                count += 1
            end
        else
            _showcompact(io, group)
        end
    end
end


text:
using Test
using Metal

dims = (3,4)
a = round.(rand(Float32, dims) * 100)
b = round.(rand(Float32, dims) * 100)
c = similar(a)

d_a = MtlArray(a)
d_b = MtlArray(b)
d_c = MtlArray(c)

len = prod(dims)

# Get the Metal device
dev = device()
src = read(joinpath(@__DIR__, "vadd.metal"), String)
# Create a Metal library from a Metal Shading Language File
lib = MtlLibrary(dev, src, MtlCompileOptions())
# Parse the compiled kernel from the library
vadd = Metal.MtlKernel(dev, lib, "add_vectors")
# Call the kernel passing the Metal arrays by their buffers
Metal.mtlcall(vadd, Tuple{Core.LLVMPtr{Float32, 1},Core.LLVMPtr{Float32, 1},Core.LLVMPtr{Float32, 1}}, d_a.buffer, d_b.buffer, d_c.buffer; threads=len)

@test a+b ≈ Array(d_c)

# Reset for next test
lib  = nothing
vadd = nothing
d_c  = MtlArray(c)

# Create a Metal library from a Metal library file
lib = MtlLibraryFromFile(dev, joinpath(@__DIR__, "vadd.metallib"))
vadd = Metal.MtlKernel(dev, lib, "add_vectors")
Metal.mtlcall(vadd, Tuple{Core.LLVMPtr{Float32, 1},Core.LLVMPtr{Float32, 1},Core.LLVMPtr{Float32, 1}}, d_a.buffer, d_b.buffer, d_c.buffer; threads=len)

@test a+b ≈ Array(d_c)


text: summation(num) = (num*(num+1)) ÷ 2

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
module ClassifierMetrics

using RecipesBase # for creating a Plots.jl recipe
using SpecialFunctions

include("utils.jl")
include("binarymetrics.jl")
include("roc.jl")
include("det.jl")
include("pr.jl")
include("plotting.jl")

end # module


text:
using GraphQLGen
using Test
using Expronicon

@testset "GraphQLGen" begin
    @testset "get_leaf_type" begin
        ex = :(Union{Vector{Union{Person,Missing,Nothing}},Missing,Nothing})
        @test GraphQLGen.get_leaf_type(ex) == :Person

        ex = :(Union{Vector{Union{Person}},Missing,Nothing})
        @test GraphQLGen.get_leaf_type(ex) == :Person

        ex = :(Vector{Union{Person,Missing,Nothing}})
        @test GraphQLGen.get_leaf_type(ex) == :Person

        ex = :(Union{Person,Missing,Nothing})
        @test GraphQLGen.get_leaf_type(ex) == :Person

        ex = :Person
        @test GraphQLGen.get_leaf_type(ex) == :Person
    end

    @testset "scalar" begin
        str = """
        scalar A
        scalar B
        scalar C
        """
        scalar_type_map = Dict(:A => :Point, :B => Vector{UInt8})
        types, _ = GraphQLGen.tojl(GraphQLGen.parse(str), scalar_type_map)
        exprs = map(GraphQLGen.ExprPrettify.prettify, types)
        @test exprs[1] == :(const A = Point)
        @test exprs[2].args[1].head == :(=)
        @test exprs[2].args[1].args[1] == :(B)
        @test exprs[2].args[1].args[2] == Vector{UInt8}
        @test exprs[3].args[1].head == :(=)
        @test exprs[3].args[1].args[1] == :(C)
        @test exprs[3].args[1].args[2] == Any
    end

    @testset "enum" begin
        str = """
        enum Episode {
          NEWHOPE
          EMPIRE
          JEDI
        }
        """
        types, _ = GraphQLGen.tojl(GraphQLGen.parse(str))
        exprs = map(GraphQLGen.ExprPrettify.prettify, types)

        ret = :(@enum Episode begin
            NEWHOPE
            EMPIRE
            JEDI
        end)

        @test exprs[1] == GraphQLGen.ExprPrettify.prettify(ret)
    end

    @testset "union" begin
        str = """
        union U = U1 | U2
        """
        types, _ = GraphQLGen.tojl(GraphQLGen.parse(str))
        exprs = map(GraphQLGen.ExprPrettify.prettify, types)
        @test exprs[1] == :(const U = Union{U1,U2})
    end

    @testset "type fields" begin
        for t in [:type, :input]
            str = """
            $t A {
              field1: Int
              field2: ID!
              field3: [ID!]!
            }
            """

            types, _ = GraphQLGen.tojl(GraphQLGen.parse(str))
            ex = GraphQLGen.ExprPrettify.prettify(types[1])
            fields = ex.args[1].args[3].args

            @test fields[1].args[1] == :field1
            @test fields[1].args[2] == :(Union{Int,Missing,Nothing})

            @test fields[2].args[1] == :field2
            @test fields[2].args[2] == :String

            @test fields[3].args[1] == :field3
            @test fields[3].args[2] == :(Vector{String})

            # Union type
        end
    end

    # need to reorder expressions when writing Julia structs
    @testset "reorder types" begin
        str = """
        union U = U1 | A

        type C {
           field1: D
        }
        input A {
          field1: [B]
        }
        type D {
           field1: ID
        }
        input B {
           field1: C
        }

        scalar S

        enum E {
           E1
           E2
        }
        """

        # DAG is A -> B -> C -> D so the output order of types
        # should be D, C, B, A.
        #
        # scalar and enum types are moved to the top
        types, _ = GraphQLGen.tojl(GraphQLGen.parse(str))
        exprs = map(GraphQLGen.ExprPrettify.prettify, types)

        @test exprs[1].args[1].args[1] == :S
        @test exprs[2].args[3] == :E
        @test exprs[3].args[1].args[2] == :D
        @test exprs[4].args[1].args[2] == :C
        @test exprs[5].args[1].args[2] == :B
        @test exprs[6].args[1].args[2] == :A
        @test exprs[7] == :(const U = Union{U1,A})
    end

    @testset "cycles" begin
        str = """
        type C {
           field1: B
        }
        type B {
           field1: A
        }
        type A {
           field1: C
           field2: A
           field3: B!
        }
        """
        types, _ = GraphQLGen.tojl(GraphQLGen.parse(str))
        exprs = map(GraphQLGen.ExprPrettify.prettify, types)
        @test exprs[1].args[1].args[2] == :A
        @test exprs[2].args[1].args[2] == :B
        @test exprs[3].args[1].args[2] == :C

        @test exprs[1].args[2].head == :function
        f = JLFunction(exprs[1].args[2])
        @test f.kwargs[1] == :($(Expr(:kw, :field1, :nothing)))
        @test f.kwargs[2] == :($(Expr(:kw, :field2, :nothing)))
        @test f.kwargs[3] == :field3

        st = exprs[1].args[1]
        fields = st.args[3].args
        @test fields[1] == :field1
        @test fields[2] == :(field2::Union{A,Missing,Nothing})
        @test fields[3] == :field3

        @test exprs[1].args[end].head == :function
        f = JLFunction(exprs[1].args[end])
        funcdef = :(Base.getproperty(t::A, sym::Symbol))
        @test f.name == :(Base.getproperty)

        body = :(
            if s === Symbol("field1")
                getfield(t, Symbol("field1"))::Union{C,Missing,Nothing}
            elseif s === Symbol("field3")
                getfield(t, Symbol("field3"))::B
            else
                getfield(t, s)
            end
        )
        @test f.body == GraphQLGen.ExprPrettify.prettify(body)
    end

    @testset "functions" begin
        @testset "query" begin
            str = """
            schema {
              query: Query
            }

            type Query {
               "this returns some books"
               books(ids: [ID!]!): Book!
               authorBooks(authorName: String!, genre: Genre, limit: Int = 10): [Book!]
            }
            """
            _, functions = GraphQLGen.tojl(GraphQLGen.parse(str))
            exprs = map(GraphQLGen.ExprPrettify.prettify, functions)

            @test exprs[1].args[1].args[1] == false # not mutable
            @test exprs[1].args[1].args[2] == :books
            @test exprs[1].args[1].args[3] == :(query::String)
            @test exprs[1].args[2] == "\"\"\"\nthis returns some books\n\"\"\""
            f = JLFunction(exprs[1].args[3])
            @test f.name == :(f::books)
            @test f.args == Any[:(ids::Vector{String})]
            @test f.kwargs == Any[]

            f = JLFunction(exprs[2].args[2])
            @test f.name == :(f::authorBooks)
            @test f.args == Any[:(authorName::String)]
            @test f.kwargs[1] ==
                  :($(Expr(:kw, :(genre::Union{Genre,Missing,Nothing}), :nothing)))
            @test f.kwargs[2] == :($(Expr(:kw, :(limit::Union{Int,Missing,Nothing}), 10)))
        end

        @testset "mutation" begin
            str = """
            schema {
              mutation: Mutation
            }

            type Mutation {
              addBooks(input: [BookInput!]!): [ID!]!
            }

            input BookInput {
                author: String!
                name: String!
                year: Int
            }
            """

            types, functions = GraphQLGen.tojl(GraphQLGen.parse(str))
            @test length(types) == 1
            @test length(functions) == 1
            exprs = map(GraphQLGen.ExprPrettify.prettify, functions)

            f = JLFunction(exprs[1].args[2])
            @test f.name == :(f::addBooks)
            @test f.args == Any[:(input::Vector{BookInput})]
            @test f.kwargs == Any[]
        end
    end

    @testset "docstrings" begin
        @testset "types" begin
            str = """

            \"""The date a film was released.\"""
            scalar ReleaseDate

            \"""A single film.\"""
            type Film implements Node {
              \"""Title of this film.\"""
              title: String

              \"""Genre of the film.\"""
              genre: Genre
            }

            type TelevisionSeries implements Node {
              \"""Title of this tv show.\"""
              title: String

              \"""Genre of the tv show.\"""
              genre: Genre
            }

            \"""Input for a film.\"""
            input FilmInput {
              \"""The title of this film.\"""
              title: String!
              \"""Genre of the film.\"""
              genre: Genre!
            }

            \"""Genre of a film.\"""
            enum Genre {
              \"""An action film.\"""
              ACTION
              \"""A comedy film.\"""
              COMEDY
            }

            \"""Filmed on camera.\"""
            union CameraUsed = Film | TelevisionSeries
            """
            types, _ = GraphQLGen.tojl(GraphQLGen.parse(str))
            exprs = map(GraphQLGen.ExprPrettify.prettify, types)

            @test exprs[1].args[3] == "The date a film was released."

            @test exprs[2].args[3] == "Genre of a film."
            @test exprs[2].args[4].args[1] == Symbol("@enum")
            @test exprs[2].args[4].args[3] == :Genre
            fields = exprs[2].args[4].args[4]
            @test fields.args[1].args[3] == "An action film."
            @test fields.args[1].args[4] == :ACTION
            @test fields.args[2].args[3] == "A comedy film."
            @test fields.args[2].args[4] == :COMEDY

            st = exprs[3].args[1]
            @test st.args[2] == :TelevisionSeries
            fields = st.args[3]
            @test fields.args[1] == "Title of this tv show."
            @test fields.args[2] == :(title::Union{String,Missing,Nothing})
            @test fields.args[3] == "Genre of the tv show."
            @test fields.args[4] == :(genre::Union{Genre,Missing,Nothing})

            @test exprs[4].args[1] == "\"\"\"\nA single film.\n\"\"\""
            st = exprs[4].args[2]
            @test st.args[2] == :Film
            fields = st.args[3]
            @test fields.args[1] == "Title of this film."
            @test fields.args[2] == :(title::Union{String,Missing,Nothing})
            @test fields.args[3] == "Genre of the film."
            @test fields.args[4] == :(genre::Union{Genre,Missing,Nothing})

            @test exprs[5].args[3] == "Filmed on camera."

            @test exprs[6].args[1] == "\"\"\"\nInput for a film.\n\"\"\""
            st = exprs[6].args[2]
            @test st.args[2] == :FilmInput
            fields = st.args[3]
            @test fields.args[1] == "The title of this film."
            @test fields.args[2] == :(title::String)
            @test fields.args[3] == "Genre of the film."
            @test fields.args[4] == :(genre::Genre)
        end

        @testset "functions" begin
            str = """
            schema {
              query: Root
            }

            type Root {
              "Fetches an object given its ID."
              node(id: ID!): Node
            }
            """
            _, functions = GraphQLGen.tojl(GraphQLGen.parse(str))
            exprs = map(GraphQLGen.ExprPrettify.prettify, functions)

            @test exprs[1].args[2] == "\"\"\"\nFetches an object given its ID.\n\"\"\""
        end
    end
end


text:
"""
    grad(fdm, f, x::AbstractVector)

Approximate the gradient of `f` at `x` using `fdm`. Assumes that `f(x)` is scalar.
"""
function grad(fdm, f, x::Vector{T}) where T<:Real
    v, dx, tmp = fill(zero(T), size(x)), similar(x), similar(x)
    for n in eachindex(x)
        v[n] = one(T)
        dx[n] = fdm(function(ϵ)
                tmp .= x .+ ϵ .* v
                return f(tmp)
            end,
            zero(T),
        )
        v[n] = zero(T)
    end
    return dx
end

"""
    jacobian(fdm, f, x::AbstractVector{<:Real}, D::Int)
    jacobian(fdm, f, x::AbstractVector{<:Real})

Approximate the Jacobian of `f` at `x` using `fdm`. `f(x)` must be a length `D` vector. If
`D` is not provided, then `f(x)` is computed once to determine the output size.
"""
function jacobian(fdm, f, x::Vector{T}, D::Int) where {T<:Real}
    J = Matrix{T}(undef, D, length(x))
    for d in 1:D
        J[d, :] = grad(fdm, x->f(x)[d], x)
    end
    return J
end
jacobian(fdm, f, x::Vector{<:Real}) = jacobian(fdm, f, x, length(f(x)))

"""
    _jvp(fdm, f, x::Vector{<:Real}, ẋ::AbstractVector{<:Real})

Convenience function to compute `jacobian(f, x) * ẋ`.
"""
_jvp(fdm, f, x::Vector{<:Real}, ẋ::AV{<:Real}) = jacobian(fdm, f, x) * ẋ

"""
    _j′vp(fdm, f, ȳ::AbstractVector{<:Real}, x::Vector{<:Real})

Convenience function to compute `jacobian(f, x)' * ȳ`.
"""
_j′vp(fdm, f, ȳ::AV{<:Real}, x::Vector{<:Real}) = jacobian(fdm, f, x, length(ȳ))' * ȳ

"""
    jvp(fdm, f, x, ẋ)

Compute a Jacobian-vector product with any types of arguments for which `to_vec` is defined.
"""
function jvp(fdm, f, (x, ẋ)::Tuple{Any, Any})
    x_vec, vec_to_x = to_vec(x)
    _, vec_to_y = to_vec(f(x))
    return vec_to_y(_jvp(fdm, x_vec->to_vec(f(vec_to_x(x_vec)))[1], x_vec, to_vec(ẋ)[1]))
end
function jvp(fdm, f, xẋs::Tuple{Any, Any}...)
    x, ẋ = collect(zip(xẋs...))
    return jvp(fdm, xs->f(xs...), (x, ẋ))
end

"""
    j′vp(fdm, f, ȳ, x...)

Compute an adjoint with any types of arguments for which `to_vec` is defined.
"""
function j′vp(fdm, f, ȳ, x)
    x_vec, vec_to_x = to_vec(x)
    ȳ_vec, _ = to_vec(ȳ)
    return vec_to_x(_j′vp(fdm, x_vec->to_vec(f(vec_to_x(x_vec)))[1], ȳ_vec, x_vec))
end
j′vp(fdm, f, ȳ, xs...) = j′vp(fdm, xs->f(xs...), ȳ, xs)

"""
    to_vec(x)

Transform `x` into a `Vector`, and return a closure which inverts the transformation.
"""
to_vec(x::Real) = ([x], first)

# Arrays.
to_vec(x::Vector{<:Real}) = (x, identity)
to_vec(x::Array) = vec(x), x_vec->reshape(x_vec, size(x))

# AbstractArrays.
function to_vec(x::T) where {T<:LinearAlgebra.AbstractTriangular}
    x_vec, back = to_vec(Matrix(x))
    return x_vec, x_vec->T(reshape(back(x_vec), size(x)))
end
to_vec(x::Symmetric) = vec(Matrix(x)), x_vec->Symmetric(reshape(x_vec, size(x)))
to_vec(X::Diagonal) = vec(Matrix(X)), x_vec->Diagonal(reshape(x_vec, size(X)...))

function to_vec(X::T) where T<:Union{Adjoint,Transpose}
    U = T.name.wrapper
    return vec(Matrix(X)), x_vec->U(permutedims(reshape(x_vec, size(X))))
end

# Non-array data structures.
function to_vec(x::Tuple)
    x_vecs, x_backs = zip(map(to_vec, x)...)
    sz = cumsum([map(length, x_vecs)...])
    return vcat(x_vecs...), function(v)
        return ntuple(n->x_backs[n](v[sz[n]-length(x_vecs[n])+1:sz[n]]), length(x))
    end
end


text:
using Documenter
using AtomEnergyLevels

DocMeta.setdocmeta!(AtomEnergyLevels, :DocTestSetup, :(using AtomEnergyLevels); recursive=true)

makedocs(
    sitename = "AtomEnergyLevels",
    format = Documenter.HTML(),
    modules = [AtomEnergyLevels]
)

if "deploy" in ARGS
  include("../../faketravis.jl")
end

deploydocs(
    repo = "github.com/malykhin-sergei/AtomEnergyLevels.jl.git",
    target = "build",
    branch = "gh-pages",
    devbranch = "master",
)