text:
mutable struct Profile{T<:AbstractMatrix,S}
    data::T
    var::DataFrame
    obs::DataFrame
    layers::Dict{Symbol,S}
    pipeline::OrderedDict
end

function Profile(data::M, var::DataFrame, obs::DataFrame; T=float(eltype(data))) where {M<:AbstractMatrix}
    @assert (nrow(var), nrow(obs)) == size(data)
    S = Union{Matrix{T},SparseMatrixCSC{T,UInt32}}
    data = T.(data)
    layers = Dict{Symbol,S}()
    pipeline = OrderedDict{Symbol,Dict}()
    Profile{typeof(data),S}(data, copy(var), copy(obs), layers, pipeline)
end

obsnames(p::Profile) = names(p.obs)
varnames(p::Profile) = names(p.var)
layernames(p::Profile) = keys(p.layers)
nrow(p::Profile) = size(p.data, 1)
ncol(p::Profile) = size(p.data, 2)
nvar(p::Profile) = size(p.data, 1)
nobs(p::Profile) = size(p.data, 2)

function Base.show(io::IO, p::Profile)
    println(io, "Profile(n_var × n_obs = ", nrow(p), " × ", ncol(p), ")")
    isempty(p.obs) || println(io, "    obs: ", join(obsnames(p), ", "))
    isempty(p.var) || println(io, "    var: ", join(varnames(p), ", "))
    isempty(p.layers) || println(io, "    layers: ", join(layernames(p), ", "))
    isempty(p.pipeline) || println(io, "    pipeline: ", join(keys(p.pipeline), ", "))
end

Base.copy(p::Profile) = Profile(copy(p.data), copy(p.var), copy(p.obs), copy(p.layers), copy(p.pipeline))

Base.maximum(p::Profile) = maximum(p.data)
Base.minimum(p::Profile) = minimum(p.data)

Base.size(p::Profile) = size(p.data)
Base.axes(p::Profile) = axes(p.data)

function Base.getindex(p::Profile, inds...)
    p_ = Profile(getindex(p.data, inds[1], inds[2]),
                 getindex(p.var, inds[1], :),
                 getindex(p.obs, inds[2], :))
    for (k, v) in p.layers
        if size(v) == size(p_.data)
            p_.layers[k] = getindex(v, inds[1], inds[2])
        end
    end
    p_.pipeline = copy(p.pipeline)
    p_
end

Base.setproperty!(prof::Profile, name::Symbol, x) = setproperty!(prof, Val(name), x)
Base.setproperty!(prof::Profile, ::Val{S}, x) where S = setfield!(prof, S, x)

function Base.setproperty!(prof::Profile, ::Val{:obs}, x)
    @assert nrow(x) == size(prof.data, 2)
    setfield!(prof, :obs, x)
end

function Base.setproperty!(prof::Profile, ::Val{:var}, x)
    @assert nrow(x) == size(prof.data, 1)
    setfield!(prof, :var, x)
end

Base.filter(x::Pair{Symbol,T}, prof::Profile) where {T} = filter!(x, copy(prof))

function Base.filter!(x::Pair{Symbol,T}, prof::Profile) where {T}
    col, f = x
    sel = f.(prof.var[:,col])
    filter!(x, prof.var)
    filter_layers!(prof, var_idx=sel)
    prof.data = prof.data[sel, :]
    prof
end

function filter_layers!(prof::Profile; var_idx=(:), obs_idx=(:))
    for k in keys(prof.layers)
        if size(prof.layers[k]) == size(prof.data)
            prof.layers[k] = prof.layers[k][var_idx, obs_idx]
        end
    end
    prof
end

function get_gene_expr(prof::Profile, gene_name, ::Nothing=nothing)
    idx = collect(1:nrow(prof))[prof.var.index .== gene_name]
    return view(prof.data, idx, :)
end

function get_gene_expr(prof::Profile, gene_name, layer::Symbol)
    idx = collect(1:nrow(prof))[prof.var.index .== gene_name]
    return view(prof.layers[layer], idx, :)
end


text:
#Functions used to summarise GWAS outputs
function calc_summaries(submodel)
    calc_tstats!(submodel)
    calc_pvalues_norm!(submodel)
end

function calc_tstats!(submodel)
    idx = findall(x -> x != -9, submodel.output_stats[:,1]) # all where MAF != -9
    if submodel.type == "additive"
        submodel.output_stats[idx,5] = submodel.output_stats[idx,3] ./ sqrt.(submodel.output_stats[idx,4])
    else
        submodel.output_stats[idx,7] = submodel.output_stats[idx,3] ./ sqrt.(submodel.output_stats[idx,4])
        submodel.output_stats[idx,8] = submodel.output_stats[idx,5] ./ sqrt.(submodel.output_stats[idx,6])
    end
end

function calc_pvalues_norm!(submodel)
    dist = Normal()
    idx = findall(x -> x != -9, submodel.output_stats[:,1])
    if submodel.type == "additive"
        submodel.output_stats[idx,6] = 2*ccdf.(dist, abs.(submodel.output_stats[idx,5]))
    else
        submodel.output_stats[idx,9] = 2*ccdf.(dist, abs.(submodel.output_stats[idx,7]))
        submodel.output_stats[idx,10] = 2*ccdf.(dist, abs.(submodel.output_stats[idx,8]))
    end
end

function write_summaries(submodel)
    open(string(submodel.output_filename, ".part", submodel.id), "a") do io
        writedlm(io, [submodel.snpset.markerInfo submodel.output_stats], '\t'; header = false)
    end
end

function concatenate_summaries(model)
    fulloutput = Any["CHR" "BP" "SNP" "REF" "ALT" "MAF" "K" "BETA_A" "VAR_A" "BETA_D" "VAR_D" "T_A" "T_D" "P_A" "P_D"] # may need a check
    if model.type == "additive"
        fulloutput = Any["CHR" "BP" "SNP" "REF" "ALT" "MAF" "K" "BETA_A" "VAR_A" "T_A" "P_A"] # may
    end
    for future in model.futures
        wait(future)
    end
    for i in 1:model.submodelcount
        filename = string(model.outputfilename, ".part", (i-1))
	if filesize(filename) > 0
	        fulloutput = vcat(fulloutput, readdlm(filename, '\t', Any, '\n'; header = false))
	end
        run(`rm $filename`)
    end
    open(string(model.outputfilename, ".txt"), "a") do io
        writedlm(io, fulloutput, '\t'; header = false)
    end
end

ext: jl

max_stars_repo_name: egmreynolds/MANA

max_issues_repo_name: egmreynolds/MANA

max_forks_repo_name: egmreynolds/MANA


text:
using ACSF
using Test, LinearAlgebra, BenchmarkTools
using JuLIP, JuLIP.Testing
using ACSF: acsf

@testset "ACSF.jl" begin
    include("test.jl")
    #include("test2.jl")
end


text:
module DNN

using BenchmarkTools
using StructuredOptimization, AbstractOperators, ProximalOperators
using PyPlot

function dnn(D::Matrix, n::Int, lambda)

	dim_in, N = size(D)
	W1 = Variable(sqrt(2/(2+n)*N)*randn(n,dim_in)) # first layer
	S1 = Sigmoid((n,N))
	b1 = Variable(1)
	W2 = Variable(sqrt(2/(2*n))*randn(n,n))        # second layer
	S2 = Sigmoid((n,N))
	b2 = Variable(1)
	W3 = Variable(sqrt(2/(n+1))*randn(1,n))        # third layer
	S3 = Sigmoid((1,N))
	b3 = Variable(1)

	L1 = S1*(W1* D.+b1)
	L2 = S2*(W2*L1.+b2)
	y  = S3*(W3*L2.+b3)

	lambda1 = lambda/(2*n)
	lambda2 = lambda/(n+n)
	lambda3 = lambda/n

	reg = lambda1*norm(W1,1)+lambda2*norm(W2,1)+lambda3*norm(W3,1)
	return y, reg
end

function create_dataset()

	srand(11)
	#construct dataset
	Na, Nb = 800,800
	N = Na+Nb
	rhoa = 5.
	a = hcat([ rhoa.*[ cos(theta);  sin(theta)] for theta in linspace(0,2*pi,Na) ]...)
	a .+= 1.5.*randn(size(a))
	b = zeros(2, div(Nb,2))
	rhob = 8.
	b = [b hcat([ rhob.*[ cos(theta);  sin(theta)] for theta in linspace(0,2*pi,div(Nb,2)) ]...)]
	b .+= randn(size(b))
	yt = [ones(Bool,Na); zeros(Bool,Nb)] #labels

	A = [a b]
	#preprocess data
	A .-= mean(A)
	A ./= sqrt.(var(A,2))

	return Na, Nb, A, yt

end

function set_up()
	Na, Nb, A, yt = create_dataset()
	KCV = false
	# training
	K = 5  # KCV folds
	L = 5   # number of candidate lambdas
	#
	n = 7  # inner layers nodes
	if KCV == true  #K-fold cross validation
		lambda = cross_validation(Na, Nb, A, yt, n, K, L)
	else
		lambda = 0.0031622776601683794 # set lambda manually/already trained
	end
	#create deep neural network
	y, reg = dnn(A, n, lambda)
	return Na, n, y, yt, A, reg
end

function run_demo()
	slv = PANOC(tol = 1e-4, maxit = 50000)
	setup = set_up()
	solve_problem!(slv,setup...)
	return setup
end

function cross_validation(Na,Nb,A,yt, n,K,L)

	Nacv, Nbcv = div(Na,K), div(Nb,K)
	fa,fb = reshape(randperm(K*Nacv),Nacv,K),reshape(randperm(K*Nbcv),Nbcv,K) #fold indices

	ytcv = [ones(Bool,Nacv);zeros(Bool,Nbcv)]              #labels validation
	yttr = [ones(Bool,(K-1)*Nacv);zeros(Bool,(K-1)*Nbcv)]  #labels training

	#initialize data matrices
	Dcv = zeros(2,Nacv+Nbcv)
	Dtr = zeros(2,(K-1)*(Nacv+Nbcv))

	ycv,     = dnn(Dcv,n,1.)
	ytr, reg = dnn(Dtr,n,1.)

	lambdas = logspace(-4,-1,L)

	cv_err = zeros(K,L) # cross validation error
	ft_err = zeros(K,L) # fit error
	@views a, b = A[:,1:Na], A[:,Na+1:end]

	for k = 1:K, l = 1:L

		if l == 1
			@printf("|------------------------------------------------------------------------| \n")
			ytr, reg = dnn(Dtr,n,1.) #re-initialize weights
		end

		@views Dcv .= [a[:,fa[:,k]] b[:,fb[:,k]]]                 # validation set
		idx_tr = find(1:K .!= k)
		@views Dtr .= [a[:,fa[:,idx_tr][:]] b[:,fb[:,idx_tr][:]]] # training set

		slv = PANOC(verbose = 0, tol = 1e-4, maxit = 5000)
		@minimize crossentropy(ytr,yttr)+lambdas[l]*reg with slv

		ft_err[k,l] = CrossEntropy(yttr)(operator(ytr)*(~).(variables(ytr)))
		cv_err[k,l] = CrossEntropy(ytcv)(operator(ycv)*(~).(variables(ytr)))

		@printf("| fold %2.1d / %2.1d | λ = %4.4f | %2.1d / %2.1d | ϵ cv =  %4.4f  | ϵ ft =  %4.4f | \n",
			k, K, lambdas[l], l, L, cv_err[k,l], ft_err[k,l] )
	end

	println("trained lambda: $( lambdas[indmin(mean(cv_err,1))] )")
	#use lambda that gives minimum cross validation error
	return lambdas[indmin(mean(cv_err,1))]
end

function solve_problem!(slv, Na, n, y, yt, A, reg)
	it, = @minimize crossentropy(y,yt)+reg with slv
	return it
end

function benchmark(;verb = 0, samples = 5, seconds = 100, tol = 1e-4, maxit = 50000 )

	suite = BenchmarkGroup()

	solvers = ["ZeroFPR",
               "PANOC",
               "PG"]
	slv_opt = [
               "(verbose = $verb, tol = $tol, maxit = $maxit)",
               "(verbose = $verb, tol = $tol, maxit = $maxit)",
               "(verbose = $verb, tol = $tol, maxit = $maxit)"]

	its = Dict([(sol,0.) for sol in solvers])
	for i in eachindex(solvers)

		setup = set_up()
		solver = eval(parse(solvers[i]*slv_opt[i]))

		suite[solvers[i]] =
		@benchmarkable(it = solve_problem!(solver, setup...),
			       setup = (
					it = 0;
					setup = deepcopy($setup);
					solver = deepcopy($solver) ),
			       teardown = (
					  $its[$solvers[$i]] = it;
					  ),
			       evals = 1, samples = samples, seconds = seconds)
	end

	results = run(suite, verbose = (verb != 0))
	println("DNN its")
	println(its)
	return results
end

function show_results(Na, n, y, yt, A, reg)

	xx = linspace(-3,3,200)
	yy = linspace(-3,3,200)

	R = vcat([[xi yi] for yi in yy, xi in xx][:]...)'
	yr,        = dnn(R,n,1.)
	regions = operator(yr)*(~).(variables(y))
	regions = reshape(regions,200,200)

	figure()
	pcolormesh(xx,yy,regions)
	cnt = contour(xx,yy,regions, levels = 0.7*[maximum(regions)] )
	plot(A[1,Na+1:end],A[2,Na+1:end], "ko")
	plot(A[1,1:Na],A[2,1:Na], "r*")
end

end


text:
using DoWhile
using Base.Test

# write your own tests here
i = 0
@do begin
    i += 1
end :while i < 10

@test i == 10

@do begin
    i -= 1
end :while i > 5

@test i == 5


text:
export isleft_ideal, isright_ideal

@doc Markdown.doc"""
    algebra(a::AbsAlgAssIdl) -> AbsAlgAss

Returns the algebra containing $a$.
"""
algebra(a::AbsAlgAssIdl) = a.algebra

iszero(a::AbsAlgAssIdl) = (a.iszero == 1)

###############################################################################
#
#  String I/O
#
###############################################################################

function show(io::IO, a::AbsAlgAssIdl)
  print(io, "Ideal of ")
  print(io, algebra(a))
  println(io, " with basis matrix")
  print(io, basis_matrix(a, copy = false))
end

################################################################################
#
#  Deepcopy
#
################################################################################

function Base.deepcopy_internal(a::AbsAlgAssIdl, dict::IdDict)
  b = typeof(a)(algebra(a))
  for i in fieldnames(typeof(a))
    if isdefined(a, i)
      if i != :algebra
        setfield!(b, i, Base.deepcopy_internal(getfield(a, i), dict))
      end
    end
  end
  return b
end

################################################################################
#
#  Basis (matrices)
#
################################################################################

function assure_has_basis(a::AbsAlgAssIdl)
  if isdefined(a, :basis)
    return nothing
  end

  A = algebra(a)
  M = basis_matrix(a, copy = false)
  a.basis = Vector{elem_type(A)}(undef, nrows(M))
  for i = 1:nrows(M)
    a.basis[i] = elem_from_mat_row(A, M, i)
  end
  return nothing
end

@doc Markdown.doc"""
    basis(a::AbsAlgAssIdl; copy::Bool = true) -> Vector{AbsAlgAssElem}

Returns the basis of $a$.
"""
function basis(a::AbsAlgAssIdl; copy::Bool = true)
  assure_has_basis(a)
  if copy
    return deepcopy(a.basis)
  else
    return a.basis
  end
end

@doc Markdown.doc"""
    basis_matrix(a::AbsAlgAssIdl; copy::Bool = true) -> MatElem

Returns the basis matrix of $a$ with respect to the basis of the algebra.
"""
function basis_matrix(a::AbsAlgAssIdl; copy::Bool = true)
  if copy
    return deepcopy(a.basis_matrix)
  else
    return a.basis_matrix
  end
end

################################################################################
#
#  Inclusion of elements in ideals
#
################################################################################

@doc Markdown.doc"""
    in(x::AbsAlgAssElem, a::AbsAlgAssIdl) -> Bool

Returns `true` if $x$ is an element of $a$ and `false` otherwise.
"""
function in(x::T, a::AbsAlgAssIdl{S, T, U}) where {S, T, U}
  A = algebra(a)
  M = matrix(base_ring(A), 1, dim(A), coefficients(x, copy = false))
  return rank(vcat(basis_matrix(a, copy = false), M)) == nrows(basis_matrix(a, copy = false)) # so far we assume nrows(basis_matrix) == rank(basis_matrix)
end

################################################################################
#
#  Test right/left
#
################################################################################

function _test_ideal_sidedness(a::AbsAlgAssIdl, side::Symbol)
  A = algebra(a)
  ba = basis(a, copy = false)
  t = A()
  for i = 1:dim(A)
    for j = 1:length(ba)
      if side == :left
        t = mul!(t, A[i], ba[j])
      elseif side == :right
        t = mul!(t, ba[j], A[i])
      else
        error("side must be either :left or :right")
      end
      if !(t in a)
        return false
      end
    end
  end
  return true
end

@doc Markdown.doc"""
    isright_ideal(a::AbsAlgAssIdl) -> Bool
    isright_ideal(a::AlgAssAbsOrdIdl) -> Bool
    isright_ideal(a::AlgAssRelOrdIdl) -> Bool

Returns `true` if $a$ is a right ideal and `false` otherwise.
"""
function isright_ideal(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl, AlgAssRelOrdIdl })
  if a.isright == 1
    return true
  elseif a.isright == 2
    return false
  end

  if _test_ideal_sidedness(a, :right)
    a.isright = 1
    return true
  end

  a.isright = 2
  return false
end

@doc Markdown.doc"""
    isleft_ideal(a::AbsAlgAssIdl) -> Bool
    isleft_ideal(a::AlgAssAbsOrdIdl) -> Bool
    isleft_ideal(a::AlgAssRelOrdIdl) -> Bool

Returns `true` if $a$ is a left ideal and `false` otherwise.
"""
function isleft_ideal(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl, AlgAssRelOrdIdl })
  if a.isleft == 1
    return true
  elseif a.isleft == 2
    return false
  end

  if _test_ideal_sidedness(a, :left)
    a.isleft = 1
    return true
  end

  a.isleft = 2
  return false
end

################################################################################
#
#  Arithmetic
#
################################################################################

@doc Markdown.doc"""
    +(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> AbsAlgAssIdl

Returns $a + b$.
"""
function +(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
  if iszero(a)
    return deepcopy(b)
  elseif iszero(b)
    return deepcopy(a)
  end

  M = vcat(basis_matrix(a), basis_matrix(b))
  r = rref!(M)
  if r != nrows(M)
    M = sub(M, 1:r, 1:ncols(M))
  end
  return ideal(algebra(a), M, :nothing, true)
end

@doc Markdown.doc"""
    *(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> AbsAlgAssIdl

Returns $a \cdot b$.
"""
function *(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where {S, T, U}
  if iszero(a)
    return deepcopy(a)
  elseif iszero(b)
    return deepcopy(b)
  end

  A = algebra(a)
  ba = basis(a, copy = false)
  bb = basis(b, copy = false)
  M = zero_matrix(base_ring(A), length(ba)*length(bb), dim(A))
  for i = 1:length(ba)
    ii = (i - 1)*length(bb)
    for j = 1:length(bb)
      elem_to_mat_row!(M, ii + j, ba[i]*bb[j])
    end
  end
  return ideal(algebra(a), M, :nothing)
end

@doc Markdown.doc"""
    ^(a::AbsAlgAssIdl, e::Union{ Int, fmpz }) -> AbsAlgAssIdl

Returns $a^e$.
"""
^(A::AbsAlgAssIdl, e::Int) = Base.power_by_squaring(A, e)
^(A::AbsAlgAssIdl, e::fmpz) = Base.power_by_squaring(A, BigInt(e))

function one(a::AbsAlgAssIdl)
  A = algebra(a)
  return ideal(A, identity_matrix(base_ring(A), dim(A)), :twosided, true)
end

function Base.copy(a::AbsAlgAssIdl)
  return a
end

function *(x::AbsAlgAssElem, a::AbsAlgAssIdl)
  @assert isleft_ideal(a) "Not a left ideal"
  if iszero(a)
    return deepcopy(a)
  end

  basis_a = basis(a, copy = false)
  return ideal_from_gens(algebra(a), [ x*basis_a[i] for i = 1:length(basis_a) ])
end

function *(a::AbsAlgAssIdl, x::AbsAlgAssElem)
  @assert isright_ideal(a) "Not a right ideal"
  if iszero(a)
    return deepcopy(a)
  end

  basis_a = basis(a, copy = false)
  return ideal_from_gens(algebra(a), [ basis_a[i]*x for i = 1:length(basis_a) ])
end
################################################################################
#
#  Equality
#
################################################################################

@doc Markdown.doc"""
    ==(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> Bool

Returns `true` if $a$ and $b$ are equal and `false` otherwise.
"""
function ==(a::AbsAlgAssIdl, b::AbsAlgAssIdl)
  algebra(a) != algebra(b) && return false
  return basis_matrix(a, copy = false) == basis_matrix(b, copy = false)
end

################################################################################
#
#  Construction
#
################################################################################

@doc Markdown.doc"""
    ideal_from_gens(A::AbsAlgAss, b::Vector{ <: AbsAlgAssElem},
                    side::Symbol = :nothing)
      -> AbsAlgAssIdl

Returns the ideal of $A$ generated by the elements of `b` as a subspace of $A$.
"""
function ideal_from_gens(A::AbsAlgAss, b::Vector{T}, side::Symbol = :nothing) where { T <: AbsAlgAssElem }
  if length(b) == 0
    M = zero_matrix(base_ring(A), 0, dim(A))
    return ideal(A, M, side, true)
  end

  @assert parent(b[1]) == A

  M = zero_matrix(base_ring(A), length(b), dim(A))
  for i = 1:length(b)
    elem_to_mat_row!(M, i, b[i])
  end
  return ideal(A, M, side)
end

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, x::AbsAlgAssElem) -> AbsAlgAssIdl

Returns the twosided principal ideal of $A$ generated by $x$.
"""
function ideal(A::AbsAlgAss, x::AbsAlgAssElem)
  t1 = A()
  t2 = A()
  M = zero_matrix(base_ring(A), dim(A)^2, dim(A))
  for i = 1:dim(A)
    t1 = mul!(t1, A[i], x)
    ii = (i - 1)*dim(A)
    for j = 1:dim(A)
      t2 = mul!(t2, t1, A[j])
      elem_to_mat_row!(M, ii + j, t2)
    end
  end

  return ideal(A, M, :twosided)
end

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, x::AbsAlgAssElem, action::Symbol) -> AbsAlgAssIdl

Returns the ideal $x \cdot A$ if `action == :left`, and $A \cdot x$ if
`action == :right`.
"""
function ideal(A::AbsAlgAss, x::AbsAlgAssElem, action::Symbol)
  M = representation_matrix(x, action)
  a = ideal(A, M)

  if action == :left
    a.isright = 1
  elseif action == :right
    a.isleft = 1
  end

  return a
end

@doc Markdown.doc"""
    *(A::AbsAlgAss, x::AbsAlgAssElem) -> AbsAlgAssIdl
    *(x::AbsAlgAssElem, A::AbsAlgAss) -> AbsAlgAssIdl

Returns the ideal $A \cdot x$ or $x \cdot A$ respectively.
"""
*(A::AbsAlgAss, x::AbsAlgAssElem) = ideal(A, x, :right)
*(x::AbsAlgAssElem, A::AbsAlgAss) = ideal(A, x, :left)

@doc Markdown.doc"""
    ideal(A::AbsAlgAss, M::MatElem, side::Symbol = :nothing, M_in_rref::Bool = false)
      -> AbsAlgAssIdl

Returns the ideal of $A$ with basis matrix $M$.
If the ideal is known to be a right/left/twosided ideal of $A$, `side` may be
set to `:right`/`:left`/`:twosided` respectively.
If `M_in_rref == true`, it is assumed that $M$ is already in row reduced echelon
form.
"""
function ideal(A::AbsAlgAss, M::MatElem, side::Symbol = :nothing, M_in_rref::Bool = false)
  @assert base_ring(M) == base_ring(A)
  @assert ncols(M) == dim(A)
  if !M_in_rref
    r, N = rref(M)
    if r == 0
      a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, zero_matrix(base_ring(A), 0, dim(A)))
      a.iszero = 1
      return a
    end
    if r != nrows(N)
      M = sub(N, 1:r, 1:ncols(N))
    else
      M = N
    end
  end
  if M_in_rref && nrows(M) == 0
    a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, M)
    a.iszero = 1
    return a
  end

  a = AbsAlgAssIdl{typeof(A), typeof(M)}(A, M)
  _set_sidedness(a, side)
  a.iszero = 2
  return a
end

# Helper function to set the side-flags
# side can be :right, :left or :twosided
function _set_sidedness(a::Union{ AbsAlgAssIdl, AlgAssAbsOrdIdl, AlgAssRelOrdIdl }, side::Symbol)
  if side == :right
    a.isleft = 0
    a.isright = 1
  elseif side == :left
    a.isleft = 1
    a.isright = 0
  elseif side == :twosided
    a.isleft = 1
    a.isright = 1
  else
    a.isleft = 0
    a.isright = 0
  end
  return nothing
end

################################################################################
#
#  Quotient rings
#
################################################################################

@doc Markdown.doc"""
    quo(A::AbsAlgAss, a::AbsAlgAssIdl) -> AbsAlgAss, AbsAlgAssMor

Returns the quotient algebra $A/a$ and the projection map $A \to A/a$.
"""
function quo(A::S, a::AbsAlgAssIdl{S, T, U}) where { S, T, U }
  @assert A == algebra(a)
  K = base_ring(A)

  # First compute the vector space quotient
  Ma = basis_matrix(a, copy = false)
  M = hcat(deepcopy(transpose(Ma)), identity_matrix(K, dim(A)))
  r = rref!(M)
  pivot_cols = Vector{Int}()
  j = 1
  for i = 1:ncols(M)
    if !iszero(M[j, i])
      if i > nrows(Ma)
        push!(pivot_cols, i - nrows(Ma))
      end
      j += 1
      if j > nrows(M)
        break
      end
    end
  end

  # We now have the basis (basis of the quotient, basis of the ideal)
  n = dim(A) - nrows(Ma)
  M = vcat(zero_matrix(K, n, dim(A)), Ma)
  oneK = K(1)
  zeroK = K()
  for i = 1:n
    M[i, pivot_cols[i]] = oneK
  end
  iM = inv(M)

  N = sub(M, 1:n, 1:dim(A))
  NN = sub(iM, 1:dim(A), 1:n)

  # Lift a basis of the quotient to A
  quotient_basis = Vector{elem_type(A)}(undef, n)
  b = zero_matrix(K, 1, n)
  for i = 1:n
    b[1, i] = oneK
    bN = b*N
    quotient_basis[i] = A([ bN[1, i] for i = 1:dim(A) ])
    b[1, i] = zeroK
  end

  # Build the multiplication table
  t = A()
  s = zero_matrix(K, 1, dim(A))
  mult_table = Array{elem_type(K), 3}(undef, n, n, n)
  for i = 1:n
    for j = 1:n
      t = mul!(t, quotient_basis[i], quotient_basis[j])
      elem_to_mat_row!(s, 1, t)
      sNN = s*NN
      mult_table[i, j, :] = [ sNN[1, k] for k = 1:n ]
    end
  end

  B = AlgAss(K, mult_table)
  AtoB = hom(A, B, NN, N)
  return B, AtoB
end

# Assumes b \subseteq a
@doc Markdown.doc"""
    quo(a::AbsAlgAssIdl, b::AbsAlgAssIdl) -> AbsAlgAss, AbsAlgAssMor

Given ideals $b \subseteq a$, this function returns the quotient algebra $a/b$
and the projection map $a \to a/b$.
"""
function quo(a::AbsAlgAssIdl{S, T, U}, b::AbsAlgAssIdl{S, T, U}) where { S, T, U }
  @assert algebra(a) == algebra(b)
  A = algebra(a)
  K = base_ring(A)

  # First compute the vector space quotient
  Ma = basis_matrix(a, copy = false)
  Mb = basis_matrix(b, copy = false)
  M = hcat(deepcopy(transpose(Mb)), deepcopy(transpose(Ma)))
  r = rref!(M)
  pivot_cols = Vector{Int}()
  j = 1
  for i = 1:ncols(M)
    if !iszero(M[j, i])
      if i > nrows(Mb)
        push!(pivot_cols, i - nrows(Mb))
      end
      j += 1
      if j > nrows(M)
        break
      end
    end
  end

  # Build the basis matrix for the quotient
  M = zero_matrix(K, dim(A), dim(A))
  n = nrows(Ma) - nrows(Mb)
  for i = 1:n
    for j = 1:dim(A)
      M[i, j] = deepcopy(Ma[pivot_cols[i], j])
    end
  end

  mult_table = _build_subalgebra_mult_table!(A, M)
  # M is now in rref

  B = AlgAss(K, mult_table)
  MM = sub(M, 1:dim(B), 1:dim(A))

  AtoB = AbsAlgAssMor{typeof(A), typeof(B), typeof(MM)}(A, B)

  N = transpose(vcat(MM, Mb)) # Another basis matrix for a
  function _image(x::AbsAlgAssElem)
    t, y = can_solve_with_solution(N, matrix(K, dim(A), 1, coefficients(x, copy = false)))
    if t
      return B([ y[i, 1] for i = 1:dim(B) ])
    else
      error("Element is not in the domain")
    end
  end

  function _preimage(x::AbsAlgAssElem)
    t = zero_matrix(K, 1, dim(B))
    for i = 1:dim(B)
      t[1, i] = x.coeffs[i]
    end
    tt = t*MM
    return A([ tt[1, i] for i = 1:dim(A) ])
  end

  AtoB.header.image = _image
  AtoB.header.preimage = _preimage
  return B, AtoB
end

################################################################################
#
#  Random elements
#
################################################################################

# TODO: implement for ::Type{AbsAlgAssIdl}
Random.gentype(a::AbsAlgAssIdl) = elem_type(algebra(a))

function rand(rng::AbstractRNG, a_sp::Random.SamplerTrivial{<:AbsAlgAssIdl})
  a = a_sp[]
  A = algebra(a)
  x = A()
  for b in basis(a, copy = false)
    x += rand(rng, base_ring(A))*b
  end
  return x
end

function rand(a::AbsAlgAssIdl, rng::UnitRange{Int})
  A = algebra(a)
  x = A()
  for b in basis(a, copy = false)
    x += rand(base_ring(A), rng)*b
  end
  return x
end

################################################################################
#
#  Reduction of element modulo ideal
#
################################################################################

function mod(x::AbsAlgAssElem, a::AbsAlgAssIdl)
  if iszero(a)
    return deepcopy(x)
  end

  c = coefficients(x)
  M = basis_matrix(a, copy = false) # Assumed to be in upper right rref
  k = 1
  for i = 1:nrows(M)
    while iszero(M[i, k])
      k += 1
    end
    if iszero(c[k])
      continue
    end

    t = divexact(c[k], M[i, k])
    for j = k:dim(algebra(a))
      c[j] = c[j] - t*M[i, j]
    end
  end
  return algebra(a)(c)
end

################################################################################
#
#  Principal generators (in full matrix algebras)
#
################################################################################

function left_principal_generator(a::AbsAlgAssIdl{S, T, U}) where { S <: AlgMat, T, U }
  @assert isleft_ideal(a) "Not a left ideal"
  A = algebra(a)
  if dim(A) != degree(A)^2*dim_of_coefficient_ring(A)
    error("Only implemented for full matrix algebras")
  end

  if iscanonical(A)
    e11 = A[1]
  else
    t = zero_matrix(coefficient_ring(A), degree(A), degree(A))
    t[1, 1] = one(coefficient_ring(A))
    e11 = A(t)
    t[1, 1] = zero(coefficient_ring(A))
  end
  ea = e11*a

  x = A()
  for i = 1:length(basis(ea, copy = false))
    if iscanonical(A)
      ei1 = A[i]
    else
      t[i, 1] = one(coefficient_ring(A))
      ei1 = A(t)
      t[i, 1] = zero(coefficient_ring(A))
    end
    x += ei1*basis(ea, copy = false)[i]
  end
  return x
end

function right_principal_generator(a::AbsAlgAssIdl{S, T, U}) where { S <: AlgMat, T, U }
  @assert isright_ideal(a) "Not a right ideal"
  A = algebra(a)
  if dim(A) != degree(A)^2*dim_of_coefficient_ring(A)
    error("Only implemented for full matrix algebras")
  end

  if iscanonical(A)
    e11 = A[1]
  else
    t = zero_matrix(coefficient_ring(A), degree(A), degree(A))
    t[1, 1] = one(coefficient_ring(A))
    e11 = A(t)
    t[1, 1] = zero(coefficient_ring(A))
  end
  ae = a*e11

  x = A()
  for i = 1:length(basis(ae, copy = false))
    if iscanonical(A)
      e1i = A[(i - 1)*degree(A) + 1]
    else
      t[1, i] = one(coefficient_ring(A))
      e1i = A(t)
      t[1, i] = zero(coefficient_ring(A))
    end
    x += basis(ae, copy = false)[i]*e1i
  end
  return x
end


text:
using Dash
using DashHtmlComponents
using DashCoreComponents

app = dash()
app.layout = html_div() do
    dcc_input(id=(type = "input", index = 1), value="input-1"),
    html_div(id="container"),
    html_div(id="output-outer"),
    html_button("Show content", id="btn")
end

callback!(app, Output("container", "children"), Input("btn", "n_clicks")) do n
    if !isnothing(n) && n > 0
        return html_div(
            [
                dcc_input(id=(type = "input", index = 2), value="input-2"),
                html_div(id="output-inner"),
            ]
        )
    else
        return "No content initially"
    end
end

function trigger_info()
    triggered = callback_context().triggered
    trig_string = !isempty(triggered) ? "Truthy" : "Falsy"
    prop_ids = join(getproperty.(triggered, :prop_id), ", ")
    return "triggered is $trig_string with prop_ids $prop_ids"
end

callback!(app,
    Output("output-inner", "children"),
    Input((type = "input", index = ALL), "value"),
) do wc_inputs
    return trigger_info()
end

callback!(app,
    Output("output-outer", "children"),
    Input((type = "input", index = ALL), "value")
    )  do value
    return trigger_info()
end

run_server(app)


text:
module Shear_Modulus

using JuLIP
using Minimise_Energy
using Deform
using Pressure
using MultiHZ

export cal_Shear_Modulus

# abstract type AbstractShearInfo end

struct ShearInfo
    G_Modulus::Float64
    ShearType::String
    ShearStrain::Array{Float64,1}
    ShearStress::Array{Float64,1}
    Config_Gamma_Hist::Array{AbstractAtoms,1}
end

function cal_Shear_Modulus(at::AbstractAtoms, StrainRate::Float64, StrainStep::Integer)
    N = StrainStep
    γ = StrainRate
    ShearStress = Array{Float64,1}(N)
    ShearStrain = γ * Vector(1:N)
    Config_Gamma_Hist = Array{AbstractAtoms,1}(N)

    deformType = "strain"
    deformValue = γ

    Config_Temp = deepcopy(at)
    for i = 1:N
        Config_Gamma_Hist[i] = deform_Conig(Config_Temp, deformType, deformValue)
        ShearStress[i] = -get_data(Config_Temp, "Pressure")[1,2]
        Config_Temp = deepcopy(Config_Gamma_Hist[i])
    end

end

end


text:
using JuMP, EAGO

                     m = Model()

                     EAGO.register_eago_operators!(m)

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

                     @variable(m, -33.59531116254298 <= q <= 31.312008728451346)

                     add_NL_constraint(m, :((0.7464435003747418 + -0.7676428720221384*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + -0.5937285886485864*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + 0.70936536775501*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.9647291374973435*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.527810271934432*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)*(1 + erf((0.7464435003747418 + -0.7676428720221384*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + -0.5937285886485864*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + 0.70936536775501*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.9647291374973435*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.527810271934432*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.35636357579153266 + 0.8665784119884141*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + 0.03175667478129274*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + -0.7799384758999897*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.4104070349373172*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + -0.7416668360909755*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.35636357579153266 + 0.8665784119884141*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + 0.03175667478129274*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + -0.7799384758999897*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.4104070349373172*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + -0.7416668360909755*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + (0.5038932713674513 + 0.7223377417739489*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + -0.17110213406275054*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + 0.7726394930743403*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.6118790711423934*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.16250007841990932*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)*(1 + erf((0.5038932713674513 + 0.7223377417739489*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + -0.17110213406275054*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + 0.7726394930743403*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.6118790711423934*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.16250007841990932*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.7190092552996492 + 0.6952008118484776*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + 0.975536624348428*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + 0.03367235210941821*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.8059867473566515*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.5874603423029647*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.7190092552996492 + 0.6952008118484776*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + 0.975536624348428*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + 0.03367235210941821*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.8059867473566515*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.5874603423029647*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 + (-0.5715134646750415 + 0.6325466436773857*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + -0.2633742432221431*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + -0.3786600783719858*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.5923927477083897*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.9114166052727355*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)*(1 + erf((-0.5715134646750415 + 0.6325466436773857*(-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))*(1 + erf((-0.6140753345551935 + 0.31490583169576514*$(x[1]) + 0.936890701125813*$(x[2]) + 0.4336839739639933*$(x[3]) + 0.7163194904782597*$(x[4]))/sqrt(2)))/2 + -0.2633742432221431*(0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))*(1 + erf((0.6840055128551685 + -0.06479708000334483*$(x[1]) + -0.1318852460587414*$(x[2]) + -0.9505832967102239*$(x[3]) + -0.4743222114166801*$(x[4]))/sqrt(2)))/2 + -0.3786600783719858*(-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))*(1 + erf((-0.704811998387469 + -0.05810722048544559*$(x[1]) + 0.19858600700809026*$(x[2]) + 0.8369128555874457*$(x[3]) + -0.8633890952295191*$(x[4]))/sqrt(2)))/2 + 0.5923927477083897*(0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))*(1 + erf((0.32884738257732327 + -0.8422524081564249*$(x[1]) + -0.01136909567652955*$(x[2]) + -0.695962689032747*$(x[3]) + -0.6201182075752167*$(x[4]))/sqrt(2)))/2 + 0.9114166052727355*(-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))*(1 + erf((-0.18842030384345643 + 0.5908436325070316*$(x[1]) + -0.482219078235679*$(x[2]) + -0.678588104206248*$(x[3]) + -0.8841369791890994*$(x[4]))/sqrt(2)))/2)/sqrt(2)))/2 - $q <= 0.0))

                     @objective(m, Min, q)

                     return m

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
module TestMomentum

using Test
import GraphIdx
import TreeLas: Momentum
import TreeLas: MGT

include("demo3x7.jl")
n1, n2 = size(y)
lam = GraphIdx.Const(0.1)

@testset "Momentum Demo 3x7              " begin
    graph = GraphIdx.Grid.GridGraph(n1, n2)
    x = Momentum.gaplas(copy(y), graph, lam, max_iter=10, Momentum.MomMem, verbose=true)
    @test all(minimum(y) .<= x)
    @test all(maximum(y) .>= x)
    if false    # compare to old implementation?
        x2 = MGT.gaplas(copy(y), graph, lam, max_iter=10, learn=0.85, verbose=true)
    end
end

end

max_stars_repo_name: EQt/treelas

max_issues_repo_name: EQt/treelas

max_forks_repo_name: EQt/treelas


text:
export Sequence, Parallel

# --- Sequence

type Sequence
    list::Array{Layer}

    Sequence() = new([])
end

function push!(s::Sequence, a_module::Layer)
    push!(s.list, a_module)
end

function forward!(s::Sequence, input)
    for a_module = s.list
        forward!(a_module, input)
    end
end

# --- Parallel

type Parallel
    list::Array{Layer}

    Parallel() = new([])
end

function push!(p::Parallel, a_module::Layer)
    push!(p.list, a_module)
end

function forward!(p::Parallel, input::Array)
    @assert length(input) == length(p.list)
    for a_module = p.list
        forward!(a_module, input)
    end
end

ListBasedContainer = Union{Parallel, Sequence}

function init!(m::ListBasedContainer)
    for submodule in m.list
        init!(submodule)
    end
end

function init_gradient!(m::ListBasedContainer)
    for submodule in m.list
        init_gradient!(submodule)
    end
end


text:
function closestpair(P::Vector{Vector{T}}) where T <: Number
    N = length(P)
    if N < 2 return (Inf, ()) end
    mindst = norm(P[1] - P[2])
    minpts = (P[1], P[2])
    for i in 1:N-1, j in i+1:N
        tmpdst = norm(P[i] - P[j])
        if tmpdst < mindst
            mindst = tmpdst
            minpts = (P[i], P[j])
        end
    end
    return mindst, minpts
end

closestpair([[0, -0.3], [1., 1.], [1.5, 2], [2, 2], [3, 3]])


text:
struct ForestFires <: Tabular end
url(::ForestFires) = "https://archive.ics.uci.edu/ml/machine-learning-databases/forest-fires/forestfires.csv"
checksum(::ForestFires) = "0d6586a1fa52f55bef48578aef14eb97273f1e9330e1a53423df497a77065253"
prep(::ForestFires) = path -> preprocess(path, ForestFires())
target(::ForestFires) = 13
headers(::ForestFires) = 1
categorical(::ForestFires) = [3;4]
size(::ForestFires) = (517, 0, 0)
problem(::ForestFires) = Regression

max_stars_repo_name: GruAna/PreprocessData

max_issues_repo_name: GruAna/PreprocessData

max_forks_repo_name: GruAna/PreprocessData


text:
module CRC32

export crc32

using Compat

function maketable(poly::UInt32)
	tab = zeros(UInt32, 256)
	for i in 0:255
		crc = @compat UInt32(i)
		for _ in 1:8
			if (crc&1) == 1
				crc = (crc >> 1) $ poly
			else
				crc >>= 1
			end
		end
		tab[i+1] = crc
	end
	tab
end

const table = maketable(0xedb88320)

function crc32(data::Vector{UInt8}, crc::Integer=0)
	crc = ~@compat UInt32(crc)
	for b in data
		crc = table[(@compat UInt8(crc&0xff) $ b) + 1] $ (crc >> 8)
	end
	~crc
end

crc32(data::String, crc::Integer=0) = crc32(convert(Vector{UInt8}, data), crc)

end # module


text:
function bm()
	# set up git repo. that's a submodule at matlab/
	tdir = @__DIR__
	mldir = joinpath(tdir,"..","matlab")
	cd(mldir)
	br = chomp(read(`git rev-parse --abbrev-ref HEAD`,String))
	if br !="bm"
		run(`git checkout bm`)
	end

	# run fedor's matlab code
	run(`/Applications/MATLAB_R2019b.app/bin/matlab -batch "bench"`)

	cd(joinpath(tdir,".."))
	# set up param dict
    pd = Dict(:na => 500,
			   :beta => 0.95,
			   :sigma => 0.35,
			   :R => 1.05,
			   :lambda => 0.000002)

	# accuracy of solution has been tested in F_test unit test.
	println()
    DCEGM.runf(par = pd);  # run once to compile
	println("julia timing:")
	@benchmark DCEGM.runf(par = $pd)
	return nothing

end


text:
# abstract type Verb end
# struct Move <: Verb end
# struct Yank <: Verb end
# struct Delete <: Verb end
# struct Change <: Verb end
# struct Find <: Verb
#     forward :: Bool
# end
# Find(;forward::Bool=true) = Find(forward)

# struct Action
#     verb :: Verb
#     quantity :: Integer
#     motion :: Motion
# end

# Action(motion :: Motion) = Action(Move(), 1, motion)
# Action(verb :: Verb, motion :: Motion) = Action(verb, 1, motion)

abstract type VimMode end
abstract type VimCommand end

# InsertMode is the standard Julia REPL
struct InsertMode <: VimMode end
struct NormalMode <: VimMode end
# A mode in which the user selects the motion which will
# be used for `Action`
# abstract type AbstractSelectMode{T <: Action} <: VimMode end
# eltype(::AbstractSelectMode{T}) where {T} = T

"""
In this mode a user selects the motion
to apply the Action to. For example, this mode
will be enabled after `d` in the command `diw`
"""
# struct MotionMode{T <: Action} <: AbstractSelectMode{T} end

# const NormalMode = MotionMode{Move}

# """
# A mode which requires the char value
# of the pressed key
# """
# struct FindChar{T <: Action} <: AbstractSelectMode{T} end
# struct ToChar{T <: Action} <: AbstractSelectMode{T} end

# struct SelectRegister <: VimCommand end


text:
_with_vdw(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in keys(vanderwaalsradius)

_with_cov(a::PDBAtom, resname_a::ASCIIString) = a.element in keys(covalentradius)

ishydrophobic(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in _hydrophobic

isaromatic(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in _aromatic

iscationic(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in _cationic

isanionic(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in _anionic

ishbonddonor(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in keys(_hbond_donor)

ishbondacceptor(a::PDBAtom, resname_a::ASCIIString) = (resname_a, a.atom) in keys(_hbond_acceptor)

"""
`any(f::Function, a::PDBResidue, b::PDBResidue, criteria::Function)`

Test if the function `f` is true for any pair of atoms between the residues `a` and `b`.
This function only test atoms that returns `true` for the fuction `criteria`.
"""
function any(f::Function, a::PDBResidue, b::PDBResidue, criteria::Function)
    resname_a, resname_b = a.id.name, b.id.name
    indices_a = find(x -> criteria(x, resname_a), a.atoms)
    indices_b = find(x -> criteria(x, resname_b), b.atoms)
    if length(indices_a) != 0 && length(indices_b) != 0
        @inbounds for i in indices_a
            for j in indices_b
                if f(a.atoms[i], b.atoms[j], resname_a, resname_b)
                    return(true)
                end
            end
        end
    end
    return(false)
end

# Interaction types
# =================

# van der Waals
# -------------

"""
Test if two atoms or residues are in van der Waals contact using:
`distance(a,b) <= 0.5 + vanderwaalsradius[a] + vanderwaalsradius[b]`.
It returns distance `<= 0.5` if the atoms aren't in `vanderwaalsradius`.
"""
function vanderwaals(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    return( distance(a,b) <= 0.5 +
               get(vanderwaalsradius, (resname_a, a.atom), 0.0) +
               get(vanderwaalsradius, (resname_b, b.atom), 0.0) )
end

vanderwaals(a::PDBResidue, b::PDBResidue) = any(vanderwaals, a, b, _with_vdw)

# van der Waals clash
# -------------------

"""
Returns `true` if the distance between the atoms is less than the sum of the `vanderwaalsradius` of the atoms.
If the atoms aren't on the list (i.e. `OXT`), the `vanderwaalsradius` of the element is used.
If there is not data in the dict,  distance `0.0` is used.
"""
function vanderwaalsclash(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    return( distance(a,b) <= get(vanderwaalsradius, (resname_a, a.atom), get(vanderwaalsradius, (resname_a, a.element), 0.0)) +
               get(vanderwaalsradius, (resname_b, b.atom), get(vanderwaalsradius, (resname_b, b.element), 0.0)) )
end

vanderwaalsclash(a::PDBResidue, b::PDBResidue) = any(vanderwaalsclash, a, b, _with_vdw)

# Covalent
# --------

"Returns `true` if the distance between atoms is less than the sum of the `covalentradius` of each atom."
function covalent(a::PDBAtom, b::PDBAtom, resname_a, resname_b) # any(... calls it with the res names
    return( distance(a,b) <= get(covalentradius, a.element, 0.0) +
               get(covalentradius, b.element, 0.0) )
end

covalent(a::PDBAtom, b::PDBAtom) = covalent(a, b, "", "")

covalent(a::PDBResidue, b::PDBResidue) = any(covalent, a, b, _with_cov)

# Disulphide
# ----------

_issulphurcys(a::PDBAtom, resname_a) = resname_a == "CYS" && a.element == "S"

"Returns `true` if two `CYS`'s `S` are at 2.08 Å or less"
function disulphide(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if _issulphurcys(a, resname_a) && _issulphurcys(b, resname_b)
        return(distance(a,b) <= 2.08)
    end
    return(false)
end

disulphide(a::PDBResidue, b::PDBResidue) = any(disulphide, a, b, _issulphurcys)

# Aromatic-Sulphur
# ----------------

_issulphur(a::PDBAtom) = a.element == "S"

"""
Returns `true` if an sulphur and an aromatic atoms are 5.3 Å or less"
"""
function aromaticsulphur(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if ( _issulphur(a) && isaromatic(b, resname_b) ) || ( _issulphur(b) && isaromatic(a, resname_a) )
        return(distance(a,b) <= 5.3)
    end
    return(false)
end

_issulphuroraromatic(a::PDBAtom, resname_a) = _issulphur(a) || isaromatic(a, resname_a)

aromaticsulphur(a::PDBResidue, b::PDBResidue) = any(aromaticsulphur, a, b, _issulphuroraromatic)

# Π-Cation
# --------

"""
There's a Π-Cation interaction if a cationic and an aromatic atoms are at 6.0 Å or less
"""
function pication(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if ( iscationic(a, resname_a) && isaromatic(b, resname_b) ) || ( iscationic(b, resname_b) && isaromatic(a, resname_a) )
        return(distance(a,b) <= 6.0)
    end
    return(false)
end

_iscationicoraromatic(a::PDBAtom, resname_a) = iscationic(a, resname_a) || isaromatic(a, resname_a)

pication(a::PDBResidue, b::PDBResidue) = any(pication, a, b, _iscationicoraromatic)

# Aromatic
# --------

"""
There's an aromatic interaction if centriods are at 6.0 Å or less.
"""
function aromatic(a::PDBResidue, b::PDBResidue)
    if a.id.name in _aromatic_res && b.id.name in _aromatic_res && distance(a, b) <= 6.0
        centres_a = _centre(_get_plane(a))
        centres_b = _centre(_get_plane(b))
        return( any( Bool[ distance(centroid_a, centroid_b) <= 6.0 for centroid_a in centres_a, centroid_b in centres_b ] ) )
    end
    false
end

# Ionic
# -----

"""
There's an ionic interaction if a cationic and an anionic atoms are at 6.0 Å or less.
"""
function ionic(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if ( iscationic(a, resname_a) && isanionic(b, resname_b) ) || ( iscationic(b, resname_b) && isanionic(a, resname_a) )
        return(distance(a,b) <= 6.0)
    end
    return(false)
end

_iscationicoranionic(a::PDBAtom, resname_a) = iscationic(a, resname_a) || isanionic(a, resname_a)

ionic(a::PDBResidue, b::PDBResidue) = any(ionic, a, b, _iscationicoranionic)

# Hydrophobic contact
# -------------------

"""
There's an hydrophobic interaction if two hydrophobic atoms are at 5.0 Å or less.
"""
function hydrophobic(a::PDBAtom, b::PDBAtom, resname_a, resname_b)
    if ishydrophobic(a, resname_a) && ishydrophobic(b, resname_b)
        return(distance(a,b) <= 5.0)
    end
    return(false)
end

hydrophobic(a::PDBResidue, b::PDBResidue) = any(hydrophobic, a, b, ishydrophobic)

# Hydrogen bonds
# --------------

function _find_antecedent(res::PDBResidue, a::PDBAtom)
    ids = _hbond_acceptor[(res.id.name, a.atom)]
    N = length(res)
    indices = Array(Int,N)
    j = 0
    @inbounds for i in 1:N
        if res.atoms[i].atom in ids
            # && (res.atoms[i].element != "H") && (res.atoms[i] != a) && covalent(res.atoms[i], a)
            j += 1
            indices[j] = i
        end
    end
    resize!(indices, j)
end

function _find_h(res::PDBResidue, a::PDBAtom)
    ids = _hbond_donor[(res.id.name, a.atom)]
    N = length(res)
    indices = Array(Int,N)
    j = 0
    @inbounds for i in 1:N
        if res.atoms[i].atom in ids
            j += 1
            indices[j] = i
        end
    end
    resize!(indices, j)
end

function _hbond_kernel(donor, acceptor, indices_donor, indices_acceptor)
    @inbounds for i in indices_donor
        don = donor.atoms[i]
        indices_h = _find_h(donor, don)
        if length(indices_h) == 0
            continue
        end
        for j in indices_acceptor
            acc = acceptor.atoms[j]
            indices_ant = _find_antecedent(acceptor, acc)
            if distance(don, acc) <= 3.9 && length(indices_ant) != 0
                for k in indices_h
                    hyd = donor.atoms[k]
                    if distance(hyd, acc) <= 2.5 && angle(don, hyd, acc) >= 90.0
                        for ant in indices_ant
                            if angle(don, acc, acceptor.atoms[ant]) >= 90.0 && angle(hyd, acc, acceptor.atoms[ant]) >= 90.0
                                return(true)
                            end
                        end
                    end
                end
            end
        end
    end
    return(false)
end

function _hydrogenbond_don_acc(donor::PDBResidue, acceptor::PDBResidue)
    if donor != acceptor
        indices_donor = find(x -> ishbonddonor(x, donor.id.name), donor.atoms)
        indices_acceptor = find(x -> ishbondacceptor(x, acceptor.id.name), acceptor.atoms)
        if length(indices_donor) != 0 && length(indices_acceptor) != 0
            return(_hbond_kernel(donor, acceptor, indices_donor, indices_acceptor))
        end
    end
    return(false)
end

"""
This function only works if there are hydrogens in the structure.
The criteria for a hydrogen bond are:

- d(Ai, Aj) < 3.9Å
- d(Ah, Aacc) < 2.5Å
- θ(Adon, Ah, Aacc) > 90°
- θ(Adon, Aacc, Aacc-antecedent) > 90°
- θ(Ah, Aacc, Aacc-antecedent) > 90°

Where Ah is the donated hydrogen atom, Adon is the hydrogen bond donor atom,
Aacc is the hydrogen bond acceptor atom and Aacc-antecednt is the atom antecedent to the hydrogen bond acceptor atom.
"""
hydrogenbond(a::PDBResidue, b::PDBResidue) = _hydrogenbond_don_acc(a,b) || _hydrogenbond_don_acc(b,a)

# STRIDE Hydrogen bonds
# =====================

function stridehydrogenbond(filename::ASCIIString; model::ASCIIString="1", group::ASCIIString="ATOM", path::ASCIIString="stride")
    out = split(readall(`$path -h $filename`),'\n')
    pairs = Set{(PDBResidueIdentifier,PDBResidueIdentifier)}()
    for line in out
        if length(line) > 3 && ( line[1:3] == "DNR" || line[1:3] == "ACC" )
            push!(pairs, (PDBResidueIdentifier(Nullable{Int}(), replace(line[12:15],' ', ""), replace(line[6:8],' ', ""), group, model, line[10:10]),
                          PDBResidueIdentifier(Nullable{Int}(), replace(line[32:35],' ', ""), replace(line[26:28],' ', ""), group, model, line[30:30])))
        end
    end
    sizehint!(pairs,length(pairs))
end

# CHIMERA Hydrogen bonds
# ======================

function chimerahydrogenbond(filename::ASCIIString; chain::ASCIIString="A", model::ASCIIString="1", commands::ASCIIString="relax 0 intraRes 0")
    sel = "select :*.$chain; select ~ligand & ~solvent" #ATOM
    out =  split(readall( `echo "select #$model; $sel; hbonds log 1 selRestrict both intermodel 0 batch 1 namingStyle simple $commands"` |> `chimera -n $filename` ), '\n')
    parser = r"^([A-Z]{3})\s+(\S+)\.(\S)\s+\S+\s+([A-Z]{3})\s+(\S+)\.(\S)\s+\S+\s+[A-Z]{3}\s+\S+\.\S\s+\S+"
    pairs = Set{(PDBResidueIdentifier,PDBResidueIdentifier)}()
    for line in out
        m = match(parser, line)
        if m !== nothing && length(m.captures) == 6
            push!(pairs, (PDBResidueIdentifier(Nullable{Int}(), m.captures[2], m.captures[1], "ATOM", model, m.captures[3]),
                          PDBResidueIdentifier(Nullable{Int}(), m.captures[5], m.captures[4], "ATOM", model, m.captures[6])))
        end
    end
    sizehint!(pairs,length(pairs))
end


text:
macro exported_function(name::Symbol)
    quote
        function $name end
        export $name
    end |> esc
end


text:
using ChainRulesCore
export logdetjac
import ChainRulesCore: frule, rrule

## Tape types and utilities

"""
    I = InvertibleOperationsTape(Y, layer_blocks, counter_block, counter_layer, logdet)

Invertible global state type, it keeps track of invertible blocks of operations (each block being a sequence of contiguous invertible layers)
"""
mutable struct InvertibleOperationsTape
    Y::Array{Any,1}
    layer_blocks::Array{Any,1}
    counter_block::Int64
    counter_layer::Int64
    logdet::Union{Nothing,Float32}
end

"""
 I = InvertibleOperationsTape() = InvertibleOperationsTape([], [], 0, 0, nothing)

Default  Constructor
"""
InvertibleOperationsTape() = InvertibleOperationsTape([], [], 0, 0, nothing)

# Initialize global state
const GLOBAL_STATE_INVOPS = InvertibleOperationsTape()
export GLOBAL_STATE_INVOPS

"""
Get current state of the tape
"""
function current(state::InvertibleOperationsTape)
    state.counter_block == 0 && throw(ArgumentError("Please, run forward pass first, to rebuild global state."))
    return state.Y[state.counter_block]
end

"""
Reset the state of the tape
"""
function reset!(state::InvertibleOperationsTape)
    state.Y = []
    state.layer_blocks = []
    state.counter_block = 0
    state.counter_layer = 0
    state.logdet = nothing
end

"""
Determine if the input is related to a new block of invertible operations
"""
isa_newblock(state::InvertibleOperationsTape, X) = (state.counter_block == 0) || !(state.Y[end] == X)

"""
Error if mismatch between state and network
"""
function check_coherence(state::InvertibleOperationsTape, net::Union{NeuralNetLayer,InvertibleNetwork})
    if state.counter_block != 0 && state.counter_layer != 0 && state.layer_blocks[state.counter_block][state.counter_layer] != net
        reset!(state)
        throw(ArgumentError("Current state does not correspond to current layer, resetting state..."))
    end
end

"""
Update state in the forward pass.
"""
function forward_update!(state::InvertibleOperationsTape, X::Array{Float32,N}, Y::Array{Float32,N}, logdet::Union{Nothing,Float32}, net::Union{NeuralNetLayer,InvertibleNetwork}) where N

    if isa_newblock(state, X)
        push!(state.Y, Y)
        push!(state.layer_blocks, [net])
        state.counter_block += 1
        state.counter_layer = 1
    else
        state.Y[state.counter_block] = Y
        push!(state.layer_blocks[state.counter_block], net)
        state.counter_layer += 1
    end
    if logdet isa Float32
        state.logdet === nothing ? (state.logdet = logdet) : (state.logdet += logdet)
    end

end

"""
Update state in the backward pass
"""
function backward_update!(state::InvertibleOperationsTape, X::Array{Float32,N}) where N

    if state.counter_layer == 1 # Check if first layer of current block
        state.Y[state.counter_block] = nothing
        state.counter_block -= 1
        (state.counter_block != 0) && (state.counter_layer = length(state.layer_blocks[state.counter_block]))
    else
        state.Y[state.counter_block] = X
        state.counter_layer -= 1
    end

    state.counter_block == 0 && reset!(state) # reset state when first block/first layer is reached

end

## Chain rules for invertible networks
# General pullback function
function pullback(net::Union{NeuralNetLayer,InvertibleNetwork}, ΔY::Array{Float32,N};
                 state::InvertibleOperationsTape=GLOBAL_STATE_INVOPS) where N

    # Check state coherency
    check_coherence(state, net)

    # Backward pass
    ΔX, X_ = net.backward(ΔY, current(state))

    # Update state
    backward_update!(state, X_)

    return nothing, ΔX
end

# Reverse-mode AD rule
function ChainRulesCore.rrule(net::Union{NeuralNetLayer,InvertibleNetwork}, X;
                              state::InvertibleOperationsTape=GLOBAL_STATE_INVOPS)

    # Forward pass
    net.logdet ? ((Y, logdet) = net.forward(X)) : (Y = net.forward(X); logdet = nothing)

    # Update state
    forward_update!(state, X, Y, logdet, net)

    # Pullback
    ∂Y_T(ΔY) = pullback(net, ΔY; state=state)

    return Y, ∂Y_T
end

## Logdet utilities for Zygote pullback

logdetjac(; state::InvertibleOperationsTape=GLOBAL_STATE_INVOPS) = state.logdet


text:
function test_generic(::Type{T}) where T <: CasadiSymbolicObject
    @testset "$( string("Transform ", T, " to null, scalar, Vector, and Matrix ") )" begin
        null_symbol = T("null_symbol", 0)
        scalar = T("scalar", 1)
        V = T("V1", 3) - T("V2", 3)
        M = T("M1", 2, 4) - T("M2", 2, 4)

        @test Vector(null_symbol) == Vector{Any}()
        @test Matrix(null_symbol) == Vector{Any}(undef,0)
        @test casadi.is_equal(T( Vector(scalar) ), scalar)
        @test casadi.is_equal(T( Matrix(scalar) ), scalar)
        @test casadi.is_equal(T( Vector(V) ), V, 1)
        @test casadi.is_equal(T( Matrix(M) ), M, 1)
    end

    @testset "$( string("Get index for ", T, "                                 ") )" begin
        s = randn(); S = T(s)
        v = randn(4); V = T(v)
        m = randn(5,5); M = T(m)

        @test to_julia(S[1])   ≈ s
        @test to_julia(S[1:1]) ≈ s
        @test to_julia(S[:])   ≈ s
        @test to_julia(V[3])   ≈ v[3]
        @test to_julia(V[2:4]) ≈ v[2:4]
        @test to_julia(V[:])   ≈ v[:]
        @test to_julia(M[3])   ≈ m[3]
        @test to_julia(M[4:9]) ≈ m[4:9]
        @test to_julia(M[:])   ≈ m[:]

        @test to_julia(S[1,1])   ≈ s
        @test to_julia(S[1:1,1]) ≈ s
        @test to_julia(S[1,:])   ≈ s
        @test to_julia(V[2,1])   ≈ v[2]
        @test to_julia(V[2:4,1]) ≈ v[2:4]
        @test to_julia(V[3,:])   ≈ v[3]
        @test to_julia(M[3,5])   ≈ m[3,5]
        @test to_julia(M[3,2:4]) ≈ m[3,2:4]
        @test to_julia(M[1:3,:]) ≈ m[1:3,:]
    end

    @testset "$( string("Set index for ", T, "                                 ") )" begin
        s = randn(); S = T(s)
        v = randn(4); V = T(v)
        m = randn(5,5); M = T(m)

        S[1] = 1.
        @test to_julia(S) == 1.
        S[1:1] = 2.
        @test to_julia(S) == 2.
        S[:] = 3.
        @test to_julia(S) == 3.
        V[2] = 1.
        v[2] = 1.
        @test to_julia(V) ≈ v
        V[2:3] = T([2. ; 3])
        v[2:3] = [2. ; 3]
        @test to_julia(V) ≈ v
        V[:] = T([1. ; 2 ; 3 ; 4])
        v[:] = [1. ; 2 ; 3 ; 4]
        @test to_julia(V) ≈ v
        M[17] = m[17] = rand()
        @test to_julia(M) ≈ m
        M[3:8] = zeros(T,6)
        m[3:8] .= 0.
        @test to_julia(M) ≈ m
        M[:] = ones(T,25)
        m[:] .= 1.
        @test to_julia(M) ≈ m

        S[1,1] = 2.
        @test to_julia(S) == 2.
        S[1,:] = 3.
        @test to_julia(S) == 3.
        S[:,1:1] = 4.
        @test to_julia(S) == 4.
        S[:,:] = 5.
        @test to_julia(S) == 5.
        V[3,1] = v[3,1] = 1.
        @test to_julia(V) ≈ v
        V[3,:] = 3.
        v[3] = 3.
        @test to_julia(V) ≈ v
        V[2:3,1:1] = 4.
        v[2:3] .= 4.
        @test to_julia(V) ≈ v
        V[:,:] = 5.
        v .= 5.
        @test to_julia(V) ≈ v
        M[3,4] = m[3,4] = 8.
        @test to_julia(M) ≈ m
        M[3,:] = 9.
        m[3,:] .= 9.
        @test to_julia(M) ≈ m
        M[2:3,4:5] = 10.
        m[2:3,4:5] .= 10.
        @test to_julia(M) ≈ m
        M[:,:] = -5.
        m .= -5.
        @test to_julia(M) ≈ m
    end

    @testset "$( string("Last index for ", T, "                                ") )" begin
        s = randn(); S = T(s)
        v = randn(4); V = T(v)
        m = randn(5,5); M = T(m)

        @test to_julia(S[end]) ≈ s
        @test to_julia(S[:,end]) ≈ s
        @test to_julia(S[end,end]) ≈ s
        @test to_julia(V[end]) ≈ v[end]
        @test to_julia(V[end,:]) ≈ v[end]
        @test to_julia(V[:,end]) ≈ v
        @test to_julia(V[end,end]) ≈ v[end]
        @test to_julia(M[end]) ≈ m[end]
        @test to_julia(M[end,:]) ≈ m[end,:]
        @test to_julia(M[:,end]) ≈ m[:,end]
        @test to_julia(M[end,end]) ≈ m[end]
    end

    @testset "$( string("One, zero, ones, zeros for ", T, "                    ") )" begin
        @test to_julia( one(T) ) == 1.
        @test to_julia( one( T("x",5,5) ) ) == one( zeros(5,5) )
        @test_throws DimensionMismatch(
          "multiplicative identity defined only for square matrices") one( zeros(T,3,4) )

        @test to_julia( zero(T) ) == 0.
        @test to_julia( zero( T("x",3,4) ) ) == zeros(3,4)

        @test to_julia( ones(T,6) ) == ones(6)
        @test to_julia( ones(T,7,8) ) == ones(7,8)

        @test to_julia( zeros(T,2) ) == zeros(2)
        @test to_julia( zeros(T,1,9) ) == zeros(1,9)
    end

    @testset "$( string("Size related operations for ", T, "                   ") )" begin
        @test size(zeros(T,3), 1)  == 3
        @test size(zeros(T,3), 2)  == 1
        @test size(ones(T,4,5), 1) == 4
        @test size(ones(T,4,5), 2) == 5
        @test_throws DimensionMismatch(
          "arraysize: dimension out of range") size(zeros(T,3,3),3)
    end

    @testset "$( string("Concatenations for ", T, "                            ") )" begin
        M = rand(3,3)

        @test to_julia( hcat([T(M) ; T(M)]) ) ≈ hcat(M, M)
        @test to_julia( vcat([T(M) ; T(M)]) ) ≈ vcat(M, M)
    end

    @testset "$( string("Matrix operations for ", T, "                         ") )" begin
        m = rand(5,3)
        s = rand(6,6)
        M = T(m)
        S = T(s)

        @test to_julia( adjoint(M) ) ≈ m'
        @test to_julia( repeat(M,2,3) ) ≈ repeat(m,2,3)
    end

    @testset "$( string("Broadcasting for ", T, "                              ") )" begin
        m₁ = rand(3,6)
        M₁ = T(m₁)
        m₂ = rand(3,6)
        M₂ = T(m₂)

        @test to_julia(M₁.*M₂)  ≈ m₁.*m₂
        @test to_julia(M₁*M₂')  ≈ m₁*m₂'
        @test to_julia(M₁'*M₂)  ≈ m₁'*m₂
        @test to_julia(M₁.^M₂)  ≈ m₁.^m₂
        @test to_julia(M₁.^3)   ≈ m₁.^3
        @test to_julia(M₁.^2.1) ≈ m₁.^2.1
    end

    @testset "$( string("Return types for ", T, "                              ") )" begin
        M = T( rand(5,3) )
        S = T( rand(6,6) )

        @test eltype( transpose(M) ) == T
        @test eltype( transpose( Matrix(M) ) ) == T
        @test eltype( Symmetric( Matrix(S) ) ) == T

        @test eltype( Symmetric( Matrix( T(rand(5,5)) ) ) * Matrix( T(rand(5,5)) ) ) == T
    end

    @testset "$( string("Solve linear systems with ", T, "                     ") )" begin
        A = Matrix(2one(T(3,3)))
        b = Vector(T([2, 4, 8]))

        @test to_julia.(A\b) ≈ [1.0, 2.0, 4.0]
    end
end


text:
using FreeType, FreeTypeAbstraction
using FreeType: FT_Vector
using Test

library = Ref{FT_Library}()
error = FT_Init_FreeType(library)
@test error == 0

@testset "test FT_Outline_Decompose" begin

refface = Ref{FT_Face}()
@test FT_New_Face(library[], joinpath(dirname(pathof(FreeTypeAbstraction)), "..", "test", "hack_regular.ttf"), 0, refface) == 0

glyph_index = FT_Get_Char_Index(refface[], 'J')
@test glyph_index == 0x00000031
@test FT_Set_Char_Size(refface[], 0, 16*64, 3, 3) == 0
@test FT_Load_Glyph(refface[], glyph_index, FT_LOAD_NO_SCALE | FT_LOAD_NO_BITMAP) == 0

function pos(p::Ptr{FT_Vector})
    v = unsafe_load(p)
    (v.x, v.y)
end

paths = []
function move_to_func(to, user)
    push!(paths, (:move, pos(to)))
    Cint(0)
end
function line_to_func(to, user)
    push!(paths, (:line, pos(to)))
    Cint(0)
end
function conic_to_func(control, to, user)
    push!(paths, (:conic, pos.((control, to))...))
    Cint(0)
end
function cubic_to_func(control1, control2, to, user)
    push!(paths, (:cubic, pos.((control1, control2, to))...))
    Cint(0)
end

move_f = @cfunction $move_to_func Cint (Ptr{FT_Vector}, Ptr{Cvoid})
line_f = @cfunction $line_to_func Cint (Ptr{FT_Vector}, Ptr{Cvoid})
conic_f = @cfunction $conic_to_func Cint (Ptr{FT_Vector}, Ptr{FT_Vector}, Ptr{Cvoid})
cubic_f = @cfunction $cubic_to_func Cint (Ptr{FT_Vector}, Ptr{FT_Vector}, Ptr{FT_Vector}, Ptr{Cvoid})

GC.@preserve move_f line_f conic_f cubic_f begin
face = unsafe_load(refface[])
glyph = unsafe_load(face.glyph)
outline_funcs = FreeType.FT_Outline_Funcs(Base.unsafe_convert.(Ptr{Cvoid}, (move_f, line_f, conic_f, cubic_f))..., 0, 0)
FT_Outline_Decompose(pointer_from_objref.((Ref(glyph.outline), Ref(outline_funcs)))..., C_NULL)
end

@test paths == [(:move, (502, -29)), (:conic, (399, -29), (307, -7)), (:conic, (210, 16), (109, 61)), (:line, (109, 297)), (:conic, (200, 216), (297, 176)), (:conic, (396, 135), (499, 135)), (:conic, (566, 135), (617, 153)), (:conic, (669, 171), (698, 210)), (:conic, (754, 284), (754, 487)), (:line, (754, 1323)), (:line, (373, 1323)), (:line, (373, 1493)), (:line, (956, 1493)), (:line, (956, 487)), (:conic, (956, 343), (930, 246)), (:conic, (905, 150), (850, 89)), (:conic, (795, 28), (709, 0)), (:conic, (624, -29), (502, -29))]

@test FT_Done_FreeType(library[]) == 0

end

## since there are no meaningful tests, we simply run FreeTypeAbstraction's tests for now
# include(joinpath(dirname(pathof(FreeTypeAbstraction)), "..", "test", "runtests.jl"))


text:
using ORInterface
using Test

@testset "ORInterface.jl" begin
    # Write your tests here.
end


text:
"""
	countAAA(a::AbstractArray, left=1, right=length(a))

Counts how many time such pattern is repeated in the array `a`.
It works through recursion.
Here is an example:

```jldoctest
julia> countAAA([0, 1, 1, 1, 0])
1
```
"""
function countAAA(a::AbstractArray, left=1, right=length(a))
	if right-left ≤ 1
		return 0
	elseif right-left == 2
		if a[left] == a[left+1] == a[right]
			return 1
		else
			return 0
		end
	end
	mid = (right+left) ÷ 2
	return countAAA(a, left, mid+1) + countAAA(a, mid, right)
end


text:
"""
    ThreadSafeVarInfo

A `ThreadSafeVarInfo` object wraps an [`AbstractVarInfo`](@ref) object and an
array of log probabilities for thread-safe execution of a probabilistic model.
"""
struct ThreadSafeVarInfo{V<:AbstractVarInfo,L} <: AbstractVarInfo
    varinfo::V
    logps::L
end
function ThreadSafeVarInfo(vi::AbstractVarInfo)
    return ThreadSafeVarInfo(vi, [Ref(zero(getlogp(vi))) for _ in 1:Threads.nthreads()])
end
ThreadSafeVarInfo(vi::ThreadSafeVarInfo) = vi

# Instead of updating the log probability of the underlying variables we
# just update the array of log probabilities.
function acclogp!(vi::ThreadSafeVarInfo, logp)
    vi.logps[Threads.threadid()][] += logp
    return vi
end

# The current log probability of the variables has to be computed from
# both the wrapped variables and the thread-specific log probabilities.
getlogp(vi::ThreadSafeVarInfo) = getlogp(vi.varinfo) + sum(getindex, vi.logps)

# TODO: Make remaining methods thread-safe.

function resetlogp!(vi::ThreadSafeVarInfo)
    for x in vi.logps
        x[] = zero(x[])
    end
    return resetlogp!(vi.varinfo)
end
function setlogp!(vi::ThreadSafeVarInfo, logp)
    for x in vi.logps
        x[] = zero(x[])
    end
    return setlogp!(vi.varinfo, logp)
end

get_num_produce(vi::ThreadSafeVarInfo) = get_num_produce(vi.varinfo)
increment_num_produce!(vi::ThreadSafeVarInfo) = increment_num_produce!(vi.varinfo)
reset_num_produce!(vi::ThreadSafeVarInfo) = reset_num_produce!(vi.varinfo)
set_num_produce!(vi::ThreadSafeVarInfo, n::Int) = set_num_produce!(vi.varinfo, n)

syms(vi::ThreadSafeVarInfo) = syms(vi.varinfo)

function setgid!(vi::ThreadSafeVarInfo, gid::Selector, vn::VarName)
    setgid!(vi.varinfo, gid, vn)
end
setorder!(vi::ThreadSafeVarInfo, vn::VarName, index::Int) = setorder!(vi.varinfo, vn, index)
setval!(vi::ThreadSafeVarInfo, val, vn::VarName) = setval!(vi.varinfo, val, vn)

keys(vi::ThreadSafeVarInfo) = keys(vi.varinfo)
haskey(vi::ThreadSafeVarInfo, vn::VarName) = haskey(vi.varinfo, vn)

link!(vi::ThreadSafeVarInfo, spl::AbstractSampler) = link!(vi.varinfo, spl)
invlink!(vi::ThreadSafeVarInfo, spl::AbstractSampler) = invlink!(vi.varinfo, spl)
islinked(vi::ThreadSafeVarInfo, spl::AbstractSampler) = islinked(vi.varinfo, spl)

getindex(vi::ThreadSafeVarInfo, spl::AbstractSampler) = getindex(vi.varinfo, spl)
getindex(vi::ThreadSafeVarInfo, spl::SampleFromPrior) = getindex(vi.varinfo, spl)
getindex(vi::ThreadSafeVarInfo, spl::SampleFromUniform) = getindex(vi.varinfo, spl)
getindex(vi::ThreadSafeVarInfo, vn::VarName) = getindex(vi.varinfo, vn)
getindex(vi::ThreadSafeVarInfo, vns::Vector{<:VarName}) = getindex(vi.varinfo, vns)

function setindex!(vi::ThreadSafeVarInfo, val, spl::AbstractSampler)
    setindex!(vi.varinfo, val, spl)
end
function setindex!(vi::ThreadSafeVarInfo, val, spl::SampleFromPrior)
    setindex!(vi.varinfo, val, spl)
end
function setindex!(vi::ThreadSafeVarInfo, val, spl::SampleFromUniform)
    setindex!(vi.varinfo, val, spl)
end

function set_retained_vns_del_by_spl!(vi::ThreadSafeVarInfo, spl::Sampler)
    return set_retained_vns_del_by_spl!(vi.varinfo, spl)
end

isempty(vi::ThreadSafeVarInfo) = isempty(vi.varinfo)
function empty!(vi::ThreadSafeVarInfo)
    empty!(vi.varinfo)
    fill!(vi.logps, zero(getlogp(vi)))
    return vi
end

function push!(
    vi::ThreadSafeVarInfo,
    vn::VarName,
    r,
    dist::Distribution,
    gidset::Set{Selector}
)
    push!(vi.varinfo, vn, r, dist, gidset)
end

function unset_flag!(vi::ThreadSafeVarInfo, vn::VarName, flag::String)
    return unset_flag!(vi.varinfo, vn, flag)
end
function is_flagged(vi::ThreadSafeVarInfo, vn::VarName, flag::String)
    return is_flagged(vi.varinfo, vn, flag)
end

tonamedtuple(vi::ThreadSafeVarInfo) = tonamedtuple(vi.varinfo)


text:
export compute_energy

let n::Int64=0
	global increment_n() = (n+=1)
	global reset_n() = (n=0)
	global get_n() = n
end
let resolve::Bool = true
	global set_resolve(r) = (resolve = r)
	global get_resolve() = resolve
end

"""
"""
function compute_energy(conf::Array{Bool,1}, g::Graph)
	length(conf) != length(g.leaves) && error("`conf` and `g` do not have the same length")
	E = 0
	if length(g.leaves) + length(g.internals) == 1
		return E
	end
	@inbounds for (i,s) in enumerate(conf)
		if s
			for k1 in 1:g.K
				# Ancestors in tree k1
				a1 = g.leaves[i].anc[k1]
				# If ancestor is identical to leaf for given configuration, go up
				# ie go up to the first non trivial split
				while is_trivial_split(a1.conf, conf) && !a1.isroot
					a1 = a1.anc::SplitNode
				end
				for k2 in (k1+1):g.K
					# Same for
					a2 = g.leaves[i].anc[k2]
					while is_trivial_split(a2.conf, conf) && !a2.isroot
						a2 = a2.anc
					end
					if get_resolve() && !are_equal_with_resolution(g, a1.conf, a2.conf, conf, k1, k2)
						E += 1
					elseif !get_resolve() && !are_equal(a1.conf, a2.conf, conf)
						E += 1
					end
				end
			end
		end
	end
	return E
end

"""
	 are_equal(nconf1, nconf2, conf)

Are internal configurations `nconf1` and `nconf2` equal for leaf configuration `conf`?
"""
function are_equal(nconf1, nconf2, conf)
	_in = max(length(nconf1), length(nconf2)) > 25 ? insorted : in

	@inbounds for n1 in nconf1
		if conf[n1] && !_in(n1, nconf2)
			return false
		end
	end

	@inbounds for n2 in nconf2
		if conf[n2] && !_in(n2, nconf1)
			return false
		end
	end

	return true
end

"""
	are_equal_with_resolution(g, aconf1, aconf2, conf, k1, k2)
"""
function are_equal_with_resolution(g, aconf1, aconf2, conf, k1::Int, k2::Int)
	if are_equal(aconf1, aconf2, conf)
		return true
	elseif is_contained(aconf1, aconf2, conf) && are_equal_with_resolution(g, aconf1, aconf2, k2, conf)
		return true
	elseif is_contained(aconf2, aconf1, conf) && are_equal_with_resolution(g, aconf2, aconf1, k1, conf)
		return true
	end
	return false
end

"""
	are_equal_with_resolution(g, aconf1, aconf2, k2, conf)

For every leaf `n` in `a1.conf`, check if all ancestors of `n` in tree `k2` up to `a2`
  have a split that is contained in `a1.conf`, for `conf` as a leaves state.
  If so, the split `a1.conf` (for `conf`) can be transformed into a clade in the other
  tree `k2` by adding one internal node.

**Expects `is_contained(a1.conf, a2.conf, conf)` to return `true`.**
"""
function are_equal_with_resolution(g::SplitGraph.Graph, aconf1, aconf2, k2::Int, conf)
	@inbounds for i in aconf1
		if conf[i]
			a = g.leaves[i].anc[k2]
			while a.conf != aconf2
				if !is_contained(a.conf, aconf1, conf)
					return false
				end
				a = a.anc::SplitNode
			end
		end
	end
	return true
end

function is_trivial_split(nodeconf, conf)
	n = 0
	@inbounds for i in nodeconf
		if conf[i]
			n += 1
		end
	end
	return n < 2
end

"""
	is_contained(nconf1, nconf2, conf)

Is `nconf1` in `nconf2` for leaf configuration `conf`?
"""
function is_contained(nconf1, nconf2, conf)
	_in = max(length(nconf2)) > 25 ? insorted : in
	@inbounds for i in nconf1
		if conf[i] && !_in(i, nconf2)
			return false
		end
	end
	return true
end

"""
"""
function compute_F(conf::Array{Bool,1}, g::Graph, γ::Real)
	E = compute_energy(conf, g)
	return E + γ*(length(conf) - sum(conf))
end

"""
"""
function doMCMC(g::Graph, conf::Array{Bool,1}, M::Int64; T=1, γ=1)
	_conf = copy(conf)
	E = compute_energy(_conf, g)
	F = E + γ*(length(conf) - sum(conf))
	##
	# ee = zeros(Int,M+1)
	# ff = zeros(Float64,M+1)
	# ee[1] = E
	# ff[1] = F
	##
	Fmin = F
	oconf = [copy(_conf)]
	for m in 1:M
		E, F = mcmcstep!(_conf, g, F, T, γ)
		# ee[m+1] = E
		# ff[m+1] = F
		# If new minimum is found
		if F < Fmin
			Fmin = F
			# oconf = [copy(_conf)]
			oconf[1] .= _conf
			deleteat!(oconf, 2:length(oconf))
		end
		# If equal minimum is found
		if F == Fmin && !in(_conf, oconf) #mapreduce(x->x!=_conf, *, oconf)
			push!(oconf, copy(_conf))
		end
	end
	return oconf, _conf, Fmin
end

"""
"""
function mcmcstep!(conf, g, F, T, γ)
	i = rand(1:length(conf))
	conf[i] = !conf[i]
	Enew = compute_energy(conf, g)
	Fnew = Enew + γ*(length(conf) - sum(conf))
	if Fnew < F || exp(-(Fnew-F)/T) > rand()
		return Enew, Fnew
	else
		conf[i] = !conf[i]
		return (round(Int64, F-γ*(length(conf) - sum(conf))), F)
	end
end

"""
	sa_opt(g::Graph; Trange=1.:-0.01:0.1, γ=1.05, M=1000, rep=1, resolve=true)

Call `_sa_opt` repeatedly to find a set of optimal confs.
"""
function sa_opt(g::Graph; Trange=reverse(1e-3:1e-2:0.6), γ=2., M=10, rep=1, resolve=true)
	set_resolve(resolve)
	#
	oconf = Any[]
	F = Int64[]
	Fmin = Inf
	nfound = 0
	for r in 1:rep
		oconf_, F_ = _sa_opt(g, γ, Trange, M)
		Fm = minimum(F_)
		if Fm == Fmin
			append!(oconf, oconf_)
			append!(F, F_)
			nfound += 1
		elseif Fm < Fmin
			Fmin = Fm
			oconf = oconf_
			F = F_
			nfound = 1
		end
	end
	return unique(oconf), F, nfound
end

function _sa_opt(g::Graph, γ, Trange, M)
	reset_chance = 0.0
	conf = ones(Bool, length(g.leaves))
	oconf = [copy(conf)]
	F = Float64[Inf]
	Fmin = F[1]
	for T in Trange
		if rand() < reset_chance
			tmp_oconf, conf, fmin = doMCMC(g, oconf[rand(1:length(oconf))], M, T=T,γ=γ)
		else
			tmp_oconf, conf, fmin = doMCMC(g, conf, M, T=T,γ=γ)
		end
		append!(F,fmin)
		# If a better conf is found than all configurations in oconf
		# (which is the min of `f` from doMCMC), completely replace oconf
		if fmin < Fmin
			oconf = tmp_oconf
			Fmin = fmin
		# If equally good confs have been found
		elseif fmin == Fmin
			append!(oconf, tmp_oconf)
			oconf = unique(oconf)
		end
	end
	return oconf,F
end

"""
	count_mismatches(g::Graph)

Count the number of topological mismatches in `g`.
  Equivalent to `compute_energy(conf, g)` with `conf = ones(Bool)`.
"""
function count_mismatches(g::Graph)
	conf = ones(Bool, length(g.leaves))
	return compute_energy(conf, g)
end

"""
	count_mismatches(t::Vararg{Tree})
"""
function count_mismatches(trees::Vararg{Tree})
	treelist = [copy(t, TreeTools.EmptyData) for t in trees]
	mcc = naive_mccs(treelist)
	mcc_names = TreeKnit.name_mcc_clades!(treelist, mcc)
	for (i,t) in enumerate(treelist)
		TreeKnit.reduce_to_mcc!(t, mcc)
	end
	g = trees2graph(treelist)
	conf = ones(Bool, length(g.leaves))

	return compute_energy(conf, g)
end

max_stars_repo_name: PierreBarrat/OldTreeKnit

max_issues_repo_name: PierreBarrat/OldTreeKnit

max_forks_repo_name: PierreBarrat/OldTreeKnit


text:
__precompile__()

module GEOGloWS

using CSV
using HTTP
using Retry
using Dates
using JSON3
using Logging
using TimeZones
using DataFrames

const baseurl = "https://geoglows.ecmwf.int/api"

"""
    _make_response(url::AbstractString, endpoint::AbstractString, parameters::AbstractDict{Symbol,<:Any})

Private function for submitting get request to particular endpoint with parameters
"""
function _make_request(
    url::AbstractString,
    endpoint::AbstractString,
    parameters::AbstractDict{Symbol,<:Any},
)

    parameterstring = "?"

    # convert dict to string of key=value
    for (key, value) in parameters
        x = "$(string(key))=$(string(value))&"
        parameterstring *= x
    end

    # cutoff the trailing '&' from the last addition of parameter
    parameterstring = parameterstring[1:end-1]

    # create a url to get data
    requesturl = joinpath(baseurl, endpoint, parameterstring)

    # send request
    _make_request(requesturl)
end

"""
    _make_response(url::AbstractString, endpoint::AbstractString)

Private function for submitting get request to particular endpoint
"""
function _make_request(url::AbstractString, endpoint::AbstractString)

    # create a url to get data
    requesturl = joinpath(baseurl, endpoint)

    # send request
    _make_request(requesturl)
end

"""
    _make_response(requesturl::AbstractString)

Private function for submitting get request with retries and exponential backoff
"""
function _make_request(requesturl::AbstractString)

    @debug "making request to $requesturl"

    # use exponential backoff with 4 retries
    @repeat 4 try
        response = HTTP.get(requesturl)

    catch e
        @delay_retry if e.status < 200 && e.status >= 500
        end
    end
end

"""
    _read_response(r::HTTP.Messages.Response, fmt::Symbol; with_dates::Bool = false)

Private function for taking api response and reading to proper format
"""
function _read_response(
    r::HTTP.Messages.Response,
    fmt::Symbol;
    with_dates::Bool = false,
)
    if fmt == :csv
        result = CSV.read(r.body, DataFrame)

        if with_dates
            result[!, :datetime] =
                DateTime.(result[!, :datetime], DateFormat("Y-m-dTHH:MM:SSz"))
        end

    elseif fmt == :json
        result = JSON3.read(r.body, Dict)
    else
        @error "could not determine return format, please choose either ':csv' or ':json'"
    end

    return result
end

include("geoglowsutils.jl")
include("forecaststreamflow.jl")
include("historicalstreamflow.jl")

end # module


text:
# Use the access token to access the Spotify Web API
# https://developer.spotify.com/documentation/general/guides/authorization-guide/#client-credentials-flow

"Access the Spotify Web API"
function spotify_request(url_ext::String, method::String= "GET"; scope = "client-credentials")
    # Not so cool, but this is how we get rid of spaces in vectors of strings:
    url_ext = replace(url_ext, ' ' => "")
    url = "https://api.spotify.com/v1/$url_ext"
    headers = header_maker(;scope)
    isnothing(headers) && return()
    resp = HTTP.Messages.Response()
    try
        resp = HTTP.request(method, url, headers)
    catch e
        knowntype = e isa HTTP.ExceptionRequest.StatusError && e.status == 401
        if knowntype
            response_body = e.response.body |> String
            @error response_body
            return nothing
        else
            @warn "HTTP.request call: method = $method\n headers = $headers \n $url_ext"
            @error e
            return nothing
        end
    end
    response_body = resp.body |> String
    response_body |> JSON3.read
end

#=
"""
Access the Spotify web api through the 'Implicit grant flow'.
Implicit grant flow is for clients that are implemented entirely
using JavaScript and running in the resource owner’s browser.
Julia: No prob.!
..but we do launch a browser for login..
"""
function other_request()
    error("ouch - not used")
    url = OAUTH_AUTHORIZE_URL
    headers = header_maker("Bearer")
    resp = HTTP.Messages.Response()
    method = "GET"
    try
        resp = HTTP.request(method, url, headers)
    catch e
        request = "HTTP.request call: method = $method\n  url = $url\n  headers = $headers \n"
        @error request
        @error e
        return nothing
    end
    response_body = resp.body |> String
    response_body |> JSON3.read
end
=#