text:
"""
@article{zitzler2000comparison,
  title={Comparison of multiobjective evolutionary algorithms: Empirical results},
  author={Zitzler, Eckart and Deb, Kalyanmoy and Thiele, Lothar},
  journal={Evolutionary computation},
  volume={8},
  number={2},
  pages={173--195},
  year={2000},
  publisher={MIT Press}
}

ZDT 测试函数的特点:
- 有固定的 2 个目标
- 有确定的 True Pareto Front
"""

abstract type ZDT <: TestProblem end

abstract type ZDT1 <: ZDT end

function ZDT1(m = 30, n_solutions = 100)
    f(x) = begin
        g = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1) )
        h = 1 - sqrt(x[1] / g)
        return [ x[1], g*h ]
    end

    bounds = Array([zeros(m) ones(m)]')

    x = range(0, 1, length=n_solutions)

    X = [vcat(x[i], zeros(m - 1)) for i in 1:n_solutions]
    pareto_set = [ Individual(x, f(x)) for x in X ]

    # @infiltrate

    return f, bounds, pareto_set
end

function truepf(x)
    "truepf"
end

thename(::Type{ZDT1}) = "ZDT1"

abstract type ZDT2 <: ZDT end

function ZDT2(D = 30, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1) )
        return [x[1], gx*(1 - (x[1] / gx)^2) ]
    end
    bounds = Array([zeros(D) ones(D)]')

    x = range(0, 1, length=n_solutions)

    x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ Individual(xx, f(xx)) for xx in x ]

    return f, bounds, pareto_set
end

abstract type ZDT3 <: ZDT end

function ZDT3(D = 30, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1) )
        a = x[1] / gx
        return [x[1], gx*(1 - sqrt(a) - a*sin(10π*x[1])) ]
    end
    bounds = Array([zeros(D) ones(D)]')

    if n_solutions < 6
        n_solutions = 6
    end

    regions = [ 0 0.0830015349
                0.182228780 0.2577623634
                0.4093136748 0.4538821041
                0.6183967944 0.6525117038
                0.8233317983 0.8518328654]

    n = Int(ceil(n_solutions / size(regions, 1)))
    x = Float64[]
    for i in 1:size(regions, 1)
        x = vcat(x, range(regions[i,1], regions[i,2], length=n))
    end

    x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ Individual(xx, f(xx)) for xx in x ]

    return f, bounds, get_non_dominated_solutions(pareto_set)
end

abstract type ZDT4 <: ZDT end

function ZDT4(D = 10, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 10*(length(x)-1) + sum( x[2:end].^2 - 10cos.(4π*x[2:end]))
        return [x[1], gx*(1 - sqrt(x[1] / gx)) ]
    end
    bounds = Array([-5zeros(D) 5ones(D)]')
    bounds[:,1] = [0, 1.0]

    x = range(0, 1, length=n_solutions)

    x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    pareto_set = [ Individual(xx, f(xx)) for xx in x ]

    return f, bounds, pareto_set
end

abstract type ZDT6 <: ZDT end

function ZDT6(D = 10, n_solutions = 100)
    f(x) = begin
        gx = 1.0 + 9.0 * ( sum(x[2:end]) / (length(x)-1.0) )^(0.25)
        ff1 = 1.0 - exp(-4.0x[1])*sin(6.0π*x[1])^6
        return [ ff1 , gx*(1.0 - (ff1 / gx)^2) ]
    end

    bounds = Array([zeros(D) ones(D)]')

    #x = range(0, 1, length=n_solutions)

    #x = [vcat(x[i], zeros(D - 1)) for i in 1:n_solutions]
    xx = range(0.2807753191, 1, length=100)
    yy = 1 .- (xx).^2
    pareto_set = [ Individual(zeros(0), [xx[i], yy[i]]) for i in 1:length(xx) ]

    return f, bounds, pareto_set
end


text:
module TestDay07

import AOC2021.Day07.part1, AOC2021.Day07.part2
using Test

test_input_raw = """
16,1,2,0,4,2,7,1,2,14"""

test_input = test_input_raw # no transformation needed here

@testset "Day 07" begin
    @testset "part 1" begin
        @test part1(test_input) == 37
    end
    @testset "part 2" begin
        @test part2(test_input) == 168
    end
end

end # module

ext: jl


text:
function problemstructure(v, draw_vertex; cell_width=10, cell_height = 10, fname="pstruct.svg")
    n = sum(v)
    d = Drawing((n+1)*cell_width, (n+1)*cell_height, fname)
    t = Table(n+1, n, cell_width, cell_height)
    sethue("black")
    line(t[1:1], t[n:1], :stroke)
    settext("\\mathbb{N}", t[1,1], valign = "bottom", halign="right")
    for (i,val) in enumerate(v)
        for s in 1:val
            draw_vertex(t[s, i])
        end
    end
    finish()
    preview()
end

max_stars_repo_name: Michiel-VL/PolyaViz

max_issues_repo_name: Michiel-VL/PolyaViz


text:
# Automatically generated using Clang.jl wrap_c, version 0.0.0

#using Compat

# begin enum ANONYMOUS_1
typealias ANONYMOUS_1 UInt32
const NVRTC_SUCCESS = (UInt32)(0)
const NVRTC_ERROR_OUT_OF_MEMORY = (UInt32)(1)
const NVRTC_ERROR_PROGRAM_CREATION_FAILURE = (UInt32)(2)
const NVRTC_ERROR_INVALID_INPUT = (UInt32)(3)
const NVRTC_ERROR_INVALID_PROGRAM = (UInt32)(4)
const NVRTC_ERROR_INVALID_OPTION = (UInt32)(5)
const NVRTC_ERROR_COMPILATION = (UInt32)(6)
const NVRTC_ERROR_BUILTIN_OPERATION_FAILURE = (UInt32)(7)
const NVRTC_ERROR_NO_NAME_EXPRESSIONS_AFTER_COMPILATION = (UInt32)(8)
const NVRTC_ERROR_NO_LOWERED_NAMES_BEFORE_COMPILATION = (UInt32)(9)
const NVRTC_ERROR_NAME_EXPRESSION_NOT_VALID = (UInt32)(10)
const NVRTC_ERROR_INTERNAL_ERROR = (UInt32)(11)
# end enum ANONYMOUS_1

# begin enum nvrtcResult
typealias nvrtcResult UInt32
const NVRTC_SUCCESS = (UInt32)(0)
const NVRTC_ERROR_OUT_OF_MEMORY = (UInt32)(1)
const NVRTC_ERROR_PROGRAM_CREATION_FAILURE = (UInt32)(2)
const NVRTC_ERROR_INVALID_INPUT = (UInt32)(3)
const NVRTC_ERROR_INVALID_PROGRAM = (UInt32)(4)
const NVRTC_ERROR_INVALID_OPTION = (UInt32)(5)
const NVRTC_ERROR_COMPILATION = (UInt32)(6)
const NVRTC_ERROR_BUILTIN_OPERATION_FAILURE = (UInt32)(7)
const NVRTC_ERROR_NO_NAME_EXPRESSIONS_AFTER_COMPILATION = (UInt32)(8)
const NVRTC_ERROR_NO_LOWERED_NAMES_BEFORE_COMPILATION = (UInt32)(9)
const NVRTC_ERROR_NAME_EXPRESSION_NOT_VALID = (UInt32)(10)
const NVRTC_ERROR_INTERNAL_ERROR = (UInt32)(11)
# end enum nvrtcResult

typealias _nvrtcProgram Void
typealias nvrtcProgram Ptr{_nvrtcProgram}


text:
function GDelta(CallPutFlag, S, X, T, r, b, v)

  d1 = (log(S / X) + (b + v^2 / 2) * T) / (v * sqrt(T))

  if CallPutFlag == "c"
    delta = exp((b - r) * T) * CND(d1)

  elseif CallPutFlag == "p"
    delta = -exp((b - r) * T) * CND(-d1)
  end

end


text:
using Documenter, MixedModels

makedocs()

deploydocs(
    deps   = Deps.pip("pygments", "mkdocs", "mkdocs-material", "python-markdown-math"),
    repo = "github.com/dmbates/MixedModels.jl.git",
    julia = "0.4",
    osname = "linux"
)


text:
if haskey(ENV, "CI") && ENV["CI"] == "true"
    ENV["PYTHON"] = ""
    using Pkg
    Pkg.build("PyCall")
end

using PyPlot
using Glob
using JLD2
using Oceananigans

using ConvergenceTests
using ConvergenceTests: compute_errors, extract_sizes

arch = CUDA.has_cuda() ? GPU() : CPU()

defaultcolors = plt.rcParams["axes.prop_cycle"].by_key()["color"]
removespine(side) = gca().spines[side].set_visible(false)
removespines(sides...) = [removespine(side) for side in sides]

filenames = glob("taylor_green*.jld2", joinpath(@__DIR__, "data"))

errors = compute_errors(
            (x, y, z, t) -> ConvergenceTests.DoublyPeriodicTaylorGreen.u(x, y, t),
            filenames...)

sizes = extract_sizes(filenames...)

Nx = map(sz -> sz[1], sizes)
L₁ = map(err -> err.L₁, errors)
L∞ = map(err -> err.L∞, errors)

close("all")
fig, ax = subplots()

ax.loglog(Nx, L₁, linestyle="None", marker="o", label="error, \$L_1\$-norm")
ax.loglog(Nx, L∞, linestyle="None", marker="^", label="error, \$L_\\infty\$-norm")
ax.loglog(Nx, L₁[end] * (Nx[end] ./ Nx).^2, "k-", linewidth=1, alpha=0.6,  label=L"\sim N_x^{-2}")

ax.set_xscale("log", base=2)
ax.set_yscale("log", base=10)

legend()

title("Convergence for freely-decaying Taylor-Green vortex")
removespines("top", "right")
ylabel("Norms of the absolute error, \$ | u_\\mathrm{simulation} - u_\\mathrm{analytical} | \$")
xlabel(L"N_x")

filename = "taylor_green_convergence_$(typeof(arch)).png"
filepath = joinpath(@__DIR__, "figs", filename)
mkpath(dirname(filepath))
savefig(filepath, dpi=480)

p = sortperm(Nx)
test_rate_of_convergence(L₁[p], Nx[p], expected=-2.0, atol=0.05, name="Taylor-Green L₁")
test_rate_of_convergence(L∞[p], Nx[p], expected=-2.0, atol=0.05, name="Taylor-Green L∞")


text:
using SpecialFunctions

"""
Standard deviation of the standard deviation of a Normal sample with size nObs

Source: http://stats.stackexchange.com/questions/631/standard-deviation-of-standard-deviation

IN
  xStd
     std dev of x
  nObs
     no of observations
"""
function std_std(xStdV :: Vector{T1}, nObsV :: Vector{T2}, dbg :: Bool) where
	{T1 <: AbstractFloat, T2 <: Integer}

	n = length(xStdV);

	stdStdV = zeros(n);
	for ix = 1 : n
	   nObs = Float64(nObsV[ix]);
	   if nObs < 300
		   # Gamma function from SpecialFunctions
		  gamma1 = gamma((nObs - 1.0) / 2.0)
		  gamma2 = gamma(nObs / 2.0)
	   else
		  # Approximation
		  #  See https://en.wikipedia.org/wiki/Gamma_function#Approximations
		  z1 = (nObs - 1.0) / 2.0 - 1.0;
		  z2 = nObs / 2.0 - 1.0;
		  gamma2 = 1.0;
		  gamma1 = exp(0.5 * log(z1) + z1 * (log(z1) - 1) - 0.5 * log(z2) - z2 * (log(z2) - 1));
	   end

	   # Ratio of gamma2 / gamma1
	   g21 = gamma2 / gamma1;

	   stdStdV[ix] = xStdV[ix] / g21 * sqrt((nObs-1) / 2 - g21 ^ 2);
	end

	##  Output check
	if dbg
	   @assert (all(stdStdV .>= 0));
	end

	return stdStdV
end


text:
"""
	arg_from_trees(t1, t2, MCCs)

Construct an `ARG` from two trees and a set of MCCs.
"""
function arg_from_trees(it1, it2, MCCs)
	if !ismissing(TreeTools.branch_length(it1.root)) || !ismissing(TreeTools.branch_length(it2.root))
		@warn "Root node of one of the input trees has a non missing branch length.
		This may cause issues when setting the branch length of the ARG"
	end
	#
	t1 = copy(it1)
	t2 = copy(it2)
	# Resolve trees
	resolve!(t1, t2, MCCs)
	# Find shared nodes in both trees
	X1, X2 = shared_nodes(t1, t2, MCCs)
	# Introduce shared singletons of one into the other
	fix_shared_singletons!(t1, t2, X1, X2, MCCs)
	# Make an initial ARG with the first tree
	arg, lm1 = arg_from_tree(t1, X1, 1)
	# Add the second tree
	lm2 = add_tree_to_arg!(arg, t2, X2, lm1, 2)
	# Point their root to the same node
	# set_global_root!(arg)
	# Reverse label map
	rlm = reverse_label_map(arg, lm1, lm2)
	# Deal with branch lengths of the ARG
	set_branch_length!(arg, t1, t2, rlm)

	return arg, rlm, lm1, lm2
end

"""
	arg_from_tree(t::Tree, SN::Dict, clr=1)

Build an initial `ARG` using one tree.
"""
function arg_from_tree(t::Tree, SN::Dict, clr=1)
	label_map = Dict{String, Any}()
	arg = ARG(Dict(), Dict(), Dict(), Dict())

	# Intialize root
	a_R = ARGNode()
	add_color!(a_R, clr)
	set_root!(a_R, clr)
	t.root.isleaf && set_label!(a_R, t.root.label)
	label_map[t.root.label] = a_R.label
	arg.nodes[a_R.label] = a_R
	arg.roots[clr] = a_R

	# Go down in tree
	for t_c in t.root.child
		grow_arg_from_treenode!(arg, a_R, t_c, SN, label_map, clr)
	end

	return arg, label_map
end

"""
	grow_arg_from_treenode!(
		arg::ARG,
		a_a::ARGNode,
		t_n::TreeTools.TreeNode,
		SN::Dict,
		label_map::Dict,
		clr,
	)

Build an initial `ARG` using one tree.
Create a new `a_n::ARGNode` corresponding to `t_n::TreeNode`, and graft it onto the current
tip of the ARG `a_a`.
If a reassortment occurs above `t_n`, a hybrid node `a_h` is created and the grafting is done
such that we have `a_a` --> `a_h` --> `a_n`.

## Inputs
- `SN::Dict` is a dictionary giving the shared/non-shared status of nodes in the trees.
  It is obtained with the function `shared_nodes`.
- `label_map::Dict`: Pointing from `TreeNode` to `ARGNode`.
"""
function grow_arg_from_treenode!(
	arg::ARG,
	a_a::AbstractARGNode,
	t_n::TreeTools.TreeNode,
	SN::Dict,
	label_map::Dict,
	clr,
)
	treenodetype = SN[t_n.label][1] # Is t_n shared / mcc_root / non shared
	root_in_other_tree = SN[t_n.label][3] # Is t_n the root in the other tree

	# Creating an ARGNode to represent t_n
	## If t_n is the root of an MCC (and not the root of the other tree)
	## Create a hybrid node a_h, a node a_n for t_n, and do a_a > a_h > a_n
	## Otherwise, create a node a_n for t_n, and do a_a > a_n
	a_n, a_h = if treenodetype == :mcc_root && !root_in_other_tree
		a_h = ARGNode(; hybrid=true)
		a_n = ARGNode()
		add_color!(a_h, clr)
		add_color!(a_n, clr)
		graft!(a_a, a_h, clr)
		graft!(a_h, a_n, clr)

		(a_n, a_h)
	else
		a_n = ARGNode()
		add_color!(a_n, clr)
		graft!(a_a, a_n, clr)

		(a_n, nothing)
	end
	t_n.isleaf && set_label!(a_n, t_n.label)

	# Adding new nodes to `arg`
	arg.nodes[label(a_n)] = a_n
	if !isnothing(a_h)
		arg.nodes[label(a_h)] = a_h
		arg.hybrids[label(a_h)] = a_h
	end
	t_n.isleaf && (arg.leaves[label(a_n)] = a_n)

	# Adding new node to label_map
	label_map[t_n.label] = label(a_n)

	# Recursing down the tree
	for t_c in t_n.child
		grow_arg_from_treenode!(
			arg, a_n, t_c, SN, label_map, clr
		)
	end

	return nothing
end

"""
	add_tree_to_arg!(arg, t, SN, lmref, clr = 2)

Add tree to an existing `ARG`, using color `clr`.
This is only designed for an ARG of two trees.
"""
function add_tree_to_arg!(arg, t, SN, lmref, clr = 2)
	lm = Dict()

	# Deal with root of t
	a_R = if SN[t.root.label][1] == :non_shared
		ARGNode()
	else
		n2 = SN[t.root.label][2] # Label of t.root in ref tree
		arg.nodes[lmref[n2]]
	end
	add_color!(a_R, clr)
	set_root!(a_R, clr)

	arg.roots[clr] = a_R
	arg.nodes[label(a_R)] = a_R
	lm[t.root.label] = a_R.label

	# Go down t
	for t_c in t.root.child
		add_treenode_to_arg!(arg, a_R, t_c, SN, lm, lmref, clr)
	end

	return lm
end

function add_treenode_to_arg!(
	arg::ARG,
	a_a::AbstractARGNode,
	t_n::TreeTools.TreeNode,
	SN::Dict,
	lm::Dict,
	lmref::Dict,
	clr
)
	treenodetype = SN[t_n.label][1]
	root_in_other_tree = SN[t_n.label][3]

	a_n = if treenodetype == :non_shared
		# New node : we create a new argnode a_n
		a_n = ARGNode()
		add_color!(a_n, clr)
		graft!(a_a, a_n, clr)
		arg.nodes[label(a_n)] = a_n
		a_n
	elseif treenodetype == :mcc_root
		# t_n is either (i) the root of the other tree, or (ii) the point of a reassortment
		# In either case, there is already a node a_n representing it
		n2 = SN[t_n.label][2]
		a_n = arg.nodes[lmref[n2]]
		add_color!(a_n, clr)
		if root_in_other_tree
			# (i)
			graft!(a_a, a_n, clr)
		else
			# (ii): Find the corresponding hybrid node
			a_h = ancestor(a_n, othercolor(clr))
			@assert ishybrid(a_h) "$label(a_h)"
			add_color!(a_h, clr)
			graft!(a_a, a_h, clr)
			graft!(a_h, a_n, clr)
		end
		a_n
	else
		# t_n is shared and not the root of an MCC
		# --> there is already an ARGNode corresponding to it
		n2 = SN[t_n.label][2]
		a_n = arg.nodes[lmref[n2]]
		add_color!(a_n, clr)
		@assert ancestor(a_n, othercolor(clr)) == a_a
		graft!(a_a, a_n, clr)
		a_n
	end

	# Adding to label map
	lm[t_n.label] = label(a_n)

	# Recursing down the tree
	for t_c in t_n.child
		add_treenode_to_arg!(
			arg, a_n, t_c, SN, lm, lmref, clr
		)
	end

	return nothing
end

"""
Map from `ARGNode` to `(TreeNode, TreeNode)`.
"""
function reverse_label_map(arg, lm1, lm2)
	rlm = Dict()
	for an in nodes(arg)
		n1 = findfirst(==(label(an)), lm1)
		n2 = findfirst(==(label(an)), lm2)
		rlm[label(an)] = (n1, n2)
		# Some checks
		if isnothing(n1) && isnothing(n2)
			@assert ishybrid(an)
		elseif !isnothing(n1) && isnothing(n2)
			@assert color(an) == [true, false]
		elseif isnothing(n1) && !isnothing(n2)
			@assert color(an) == [false, true]
		else
			@assert color(an) == [true, true]
		end
	end

	return rlm
end
"""
	set_branch_length!(arg, t1, t2, X1, X2)

Set branch lengths of `arg` following a simple heuristic.
For a given node `a` of the ARG:
- If `a` and the branch above it are fully shared, average of tree branch lengths
- If `a` belongs only to one tree, then use this tree for the branch length
- If `a` is the root of an MCC and `a_h` the reassortment above it, then let `τ` be the minimum branch length above `a` in either of the trees.
  Place `a_h` at `τ/2` above `a`, and use individual trees for the length above `a_h`.

### Note:
Due to missing values propagating, we have to deal with them carefuly.
"""
function set_branch_length!(arg::ARG, t1, t2, lm)
	v1 = Dict()
	v2 = Dict()
	for a in nodes(arg)
		if !ishybrid(a)
			n1, n2  = lm[label(a)]
			if isnothing(n1)
				# Use n2 only
				set_branch_length!(a, t2.lnodes[n2].tau, 2)
			elseif isnothing(n2)
				# Use n1 only
				set_branch_length!(a, t1.lnodes[n1].tau, 1)
			else
				# Shared by n1 and n2
				# Is the node above a hybrid?
				if ishybrid(ancestor(a, 1))
					# MCC root
					a_h = ancestor(a, 1)
					τ = if ismissing(t1.lnodes[n1].tau)
						t2.lnodes[n2].tau
					elseif ismissing(t2.lnodes[n2].tau)
						t1.lnodes[n1].tau
					else
						min(t1.lnodes[n1].tau, t2.lnodes[n2].tau)
					end
					aτ = τ/2
					hτ1 = t1.lnodes[n1].tau - aτ
					hτ2 = t2.lnodes[n2].tau - aτ
					set_branch_length!(a, aτ, 1, 2)
					set_branch_length!(a_h, hτ1, 1)
					set_branch_length!(a_h, hτ2, 2)
				else
					# Fully shared node - if root of one of the trees, deal with missing tau
					τ1, τ2 = if ismissing(t1.lnodes[n1].tau)
						if t1.lnodes[n1].isroot
							(missing, t2.lnodes[n2].tau)
						else
							(t2.lnodes[n2].tau, t2.lnodes[n2].tau)
						end
					elseif ismissing(t2.lnodes[n2].tau)
						if t2.lnodes[n2].isroot
							(t1.lnodes[n1].tau, missing)
						else
							(t1.lnodes[n1].tau, t1.lnodes[n1].tau)
						end
					else
						τ = (t1.lnodes[n1].tau + t2.lnodes[n2].tau)/2
						(τ, τ)
					end
					set_branch_length!(a, τ1, 1)
					set_branch_length!(a, τ2, 2)
				end
			end
		end
	end
end

"""
	shared_nodes(t1::Tree, t2::Tree, MCCs)

Return two dictionaries containing the status of nodes for each tree.
*E.g.* `S::Dict` for `t1` such that for `n1` a node in `t1`, `S[n1] = (status, n2, isroot, id_mcc)` with:
- `status ∈ {:non_shared, :shared, :mcc_root, :shared_singleton}`
- `n2`: if `status==:shared` or `status==:mcc_root`, node in `t2`. Otherwise, `nothing`.
- `isroot` if `n2` is the root of `t2` (i.e. root of other tree)
- `id_mcc`: if `status!=:non_shared`, identifier for MCC
"""
function shared_nodes(t1::Tree, t2::Tree, MCCs)
	X1 = Dict()
	X2 = Dict()
	for (i,mcc) in enumerate(MCCs)
		shared_nodes!(X1, X2, t1, t2, mcc, i)
	end
	for n1 in internals(t1)
		if !haskey(X1, n1.label)
			X1[n1.label] = (:non_shared, nothing, false, nothing)
		end
	end
	for n2 in internals(t2)
		if !haskey(X2, n2.label)
			X2[n2.label] = (:non_shared, nothing, false, nothing)
		end
	end

	# Quick test of reflectivity
	for (k1, v1) in X1
		if v1[1] == :shared && (X2[v1[2]][1] != :shared || X2[v1[2]][2] != k1)
			error("Inconsistent shared nodes")
		end
	end

	return X1, X2
end
#=
Shared nodes `X` of tree `t`, using `t2` as a second tree, for a given mcc.
=#
function shared_nodes!(X1::Dict, X2::Dict, t1::Tree, t2::Tree, mcc, id_mcc)
	# To identify singletons, we'll need the splits restricted to the mcc
	Sm = TreeKnit.splits_in_mcc(mcc, t1, t2)

	# Deal with the root first
	A1 = lca(t1, mcc)
	A2 = lca(t2, mcc)
	X1[A1.label] = (:mcc_root, A2.label, A2.isroot, id_mcc)
	X2[A2.label] = (:mcc_root, A1.label, A1.isroot, id_mcc)

	# Now adding internal nodes
	for label in mcc
		# Going up the two trees at the same time starting from this leaf
		n1 = t1.lnodes[label]
		n2 = t2.lnodes[label]
		#  A leaf is always shared
		!haskey(X1, n1.label) && (X1[n1.label] = (:shared, n2.label, n2.isroot, id_mcc))
		!haskey(X2, n2.label) && (X2[n2.label] = (:shared, n1.label, n1.isroot, id_mcc))

		a1 = n1.anc
		a2 = n2.anc
		while (n1 != A1 || n2 != A2) && (!haskey(X1, a1.label) || !haskey(X2, a2.label))
			# Are a1/a2 shared / shared singleton?
			f1 = is_shared_singleton(a1, n1, Sm[1]) ? :shared_singleton : :shared
			f2 = is_shared_singleton(a2, n2, Sm[2]) ? :shared_singleton : :shared
			if f1 == :shared && f2 == :shared
				# Both are shared, they are equivalent. Go up in the trees
				!haskey(X1, a1.label) && (X1[a1.label] = (:shared, a2.label, a2.isroot, id_mcc))
				!haskey(X2, a2.label) && (X2[a2.label] = (:shared, a1.label, a1.isroot, id_mcc))
				n1 = a1
				n2 = a2
				a1 = n1.anc
				a2 = n2.anc
				continue
			end
			if f1 == :shared_singleton
				# a1 does not have an equivalent in t2, but is shared. Go up in t1
				X1[a1.label] = (:shared_singleton, nothing, false, id_mcc)
				n1 = a1
				a1 = n1.anc
			end
			if f2 == :shared_singleton
				# a2 does not have an equivalent in t1, but is shared. Go up in t2
				X2[a2.label] = (:shared_singleton, nothing, false, id_mcc)
				n2 = a2
				a2 = n2.anc
			end
		end
	end
end

"""
	is_shared_singleton(n::TreeNode, c::TreeNode, Sm::SplitList)

Check whether `n` is a shared node that would be a singleton in one of the trees.
This is the case if
	- the splits defined by `n` and `c` are the same when restricted to an mcc
	- the split defined by `n` is a leaf split when restricted to an mcc

It's necessary to check this second condition since leaves are not stored in `SplitList`
  objects.

## Input
  - node `n`
  - node `c`, child of `n`
  - `Sm::TreeTools.SplitList`: splits defined by the MCC `c` belongs to.
"""
function is_shared_singleton(n::TreeNode, c::TreeNode, Sm::SplitList)
	if c.anc != n
		error("$(c.label) is not a child of $(n.label)")
	end
	if c.isleaf
		return TreeTools.is_leaf_split(Sm.splitmap[n.label], Sm.mask)
	else
		return Base.isequal(Sm.splitmap[n.label], Sm.splitmap[c.label], Sm.mask)
	end
end

"""
	fix_shared_singletons!(t1, t2, X1, X2, MCCs)

Introduce shared singletons of one tree in the other, to make reconstructing the ARG easier.
"""
function fix_shared_singletons!(t1, t2, X1, X2, MCCs)
	for (i,mcc) in enumerate(MCCs)
		fix_shared_singletons!(t1, t2, X1, X2, mcc)
	end
end
"""
	fix_shared_singletons!(t, Xref, tref)

Introduce shared singletons of `tref` into `t`.
"""
function fix_shared_singletons!(t1, t2, X1, X2, mcc::Vector{String})
	A1 = lca(t1, mcc...)
	A2 = lca(t2, mcc...)
	id_mcc = X1[A1.label][4]
	@assert X1[A1.label][1] == :mcc_root "$(A1)"
	@assert X2[A2.label][1] == :mcc_root "$(A2)"
	@assert X1[A1.label][2] == A2.label
	@assert X2[A2.label][2] == A1.label

	v1 = Dict()
	v2 = Dict()
	# Go up from each leaf in the two trees
	for m in mcc
		n1 = t1.lnodes[m]
		n2 = t2.lnodes[m]
		while n1 != A1 && !haskey(v1, n1.label) && n2 != A2 && !haskey(v2, n2.label)
			# look up in the two trees
			a1 = n1.anc
			a2 = n2.anc
			f1 = (X1[a1.label][1] == :shared_singleton)
			f2 = (X2[a2.label][1] == :shared_singleton)
			# println("$(n1.label) - $(n1.tau) - $(a1.label)")
			# println("$(n2.label) - $(n2.tau) - $(a2.label)")
			if f1
				if f2
					# Both are shared.
					# Choose one and fix it
					if !ismissing(n1.tau) && !ismissing(n2.tau)
						if n1.tau < n2.tau
							# a1 in t2
							introduce_singleton!(n2, a2, a1, n1.tau, X2, X1)
						else
							# a2 in t1
							introduce_singleton!(n1, a1, a2, n2.tau, X1, X2)
						end
					else
						# a1 in t2 by default
						introduce_singleton!(n2, a2, a1, n1.tau, X2, X1)
					end

				else
					# Introduce a1 in t2
					introduce_singleton!(n2, a2, a1, n1.tau, X2, X1)
				end
			elseif f2
				# Introduce a2 in a1
				introduce_singleton!(n1, a1, a2, n2.tau, X1, X2)
			end
			# At this point, n1 and n2 should have the same ancestor
			v1[n1.label] = true
			v2[n2.label] = true
			n1 = n1.anc
			n2 = n2.anc
			@assert X1[n1.label][2] == n2.label "1 - $(n1.label) - $(X1[n1.label]) - $(n2.label)"
			@assert X2[n2.label][2] == n1.label "2 - $(n2.label) - $(X2[n2.label])"
		end
	end
	node2tree!(t1, t1.root)
	node2tree!(t2, t2.root)
end

function introduce_singleton!(n, a, sref, τ, X, Xref)
	# println("Trying $(n.label) < $(sref.label) < $(a.label)")
	s = TreeTools.add_internal_singleton!(n, a, τ, make_random_label("Singleton"))
	X[s.label] = (:shared, sref.label, false, X[n.label][4])
	Xref[sref.label] = (:shared, s.label, false, Xref[sref.label][4])

	return nothing
end

max_issues_repo_name: PierreBarrat/TreeKnit

max_forks_repo_name: PierreBarrat/TreeKnit


text:
module SearchablePDFs

using Pkg
using Pkg.Artifacts
using Random

using ProgressMeter
using Scratch
using CSV
using Comonicon

using Poppler_jll
using unpaper_jll
using Tesseract_jll

export ocr

#####
##### Utilities
#####

function argument_error(msg; exception=isinteractive())
    if exception
        throw(ArgumentError(msg))
    else
        printstyled("ERROR:"; bold=true, color=:red)
        printstyled(" ", msg, "\n"; color=:red)
        exit(1)
    end
end

function require_extension(path, ext; exception=isinteractive())
    _ext = splitext(path)[2]
    _ext == ext ||
        argument_error("Expected $path to have file extension `$ext`; got `$(_ext)`";
                       exception)
    return nothing
end

function require_no_file(path; exception=isinteractive())
    isfile(path) && argument_error("File already exists at `$(path)`!"; exception)
    return nothing
end

# For now we will hardcode this choice of training data
function get_data_path()
    return artifact"tessdata_fast"
end

# a place to store intermediate files; we could use temporary directories,
# but I've occasionally run into permissions issues there
# so I'd prefer to use a more local location. This also means that we are
# more in charge of the cleanup, which can be good for debugging.
function get_scratch_dir(pdf)
    return joinpath(@get_scratch!("pdf_tmps"),
                    splitext(basename(pdf))[1] * "_" * string(randstring(10)))
end

# https://discourse.julialang.org/t/collecting-all-output-from-shell-commands/15592/7
function run_and_collect_logs(cmd::Cmd)
    out = Pipe()
    err = Pipe()
    process = run(pipeline(cmd; stdout=out, stderr=err); wait=false)
    close(out.in)
    close(err.in)

    stdout = @async String(read(out))
    stderr = @async String(read(err))
    wait(process)
    return (stdout=fetch(stdout), stderr=fetch(stderr), code=process.exitcode)
end

# There's gotta be a better way...
function num_pages(pdf)
    result = read(`$(Poppler_jll.pdfinfo()) $pdf`, String)
    m = match(r"Pages\:\s*([0-9]*)", result)
    return parse(Int, m.captures[1])
end

#####
##### Step 1: Extract images
#####

# Use Poppler to extract the image
function get_images(pdf, page_range::UnitRange{Int}, tmp, total_pages)
    logs = run_and_collect_logs(`$(Poppler_jll.pdftoppm()) -f $(first(page_range)) -l $(last(page_range)) $pdf -tiff -forcenum $(tmp)/page`)
    @debug "`pdftoppm`" logs
    paths = [joinpath(tmp, string("page-", lpad(page, ndigits(total_pages), '0'), ".tif"))
             for page in page_range]
    return paths, (; binary="pdftoppm", logs...)
end

# Clean up an image with unpaper
function unpaper(img)
    img_base, img_ext = splitext(img)
    img_unpaper = img_base * "_unpaper" * img_ext
    logs = run_and_collect_logs(`$(unpaper_jll.unpaper()) $img $img_unpaper`)
    return (; img_unpaper, logs=(; binary="unpaper", logs...))
end

#####
##### Step 2: Use tesseract to generate a one-page searchable PDF from an image
#####

function make_pdf(img; tesseract_nthreads)
    data_path = get_data_path() * "/"
    img_base, img_ext = splitext(img)
    output = img_base
    tesseract = addenv(Tesseract_jll.tesseract(), "OMP_THREAD_LIMIT" => tesseract_nthreads)
    cmd = `$tesseract -l eng+equ --tessdata-dir $data_path $img $output -c tessedit_create_pdf=1`
    @debug "Tesseracting!" img
    logs = run_and_collect_logs(cmd)
    @debug logs
    return (; pdf=output * ".pdf", logs=(; binary="tesseract", logs...))
end

#####
##### Step 3: collect all the PDFs into one with `pdfunite`
#####

function unite_pdfs(pdfs, output)
    logs = run_and_collect_logs(`$(Poppler_jll.pdfunite()) $pdfs $output`)
    return (; binary="pdfunite", logs...)
end

# unites all the pdfs in `pdfs`, `max_files` at a time.
# I ran into "too many open files" errors otherwise
# (which seems weird... maybe <https://github.com/JuliaLang/julia/issues/31126>? It was on MacOS)
function unite_many_pdfs!(unite_progress_meter, all_logs, tmp, pdfs, output_path;
                          max_files_per_unite=100)
    isdir(tmp) || mkdir(tmp)

    output_paths = map(enumerate(Iterators.partition(pdfs, max_files_per_unite))) do (i,
                                                                                current_pdfs)
        current_output_path = joinpath(tmp, string("section_", i, ".pdf"))
        unite_logs = unite_pdfs(current_pdfs, current_output_path)
        put!(all_logs, (; page=missing, unite_logs...))
        next!(unite_progress_meter; step=length(current_pdfs))
        return current_output_path
    end

    unite_logs = unite_pdfs(output_paths, output_path)
    put!(all_logs, (; page=missing, unite_logs...))

    next!(unite_progress_meter; step=length(output_paths))
    return nothing
end

#####
##### Apply steps 1 -- 3
#####

"""
    ocr(pdf, output_path=string(splitext(pdf)[1], "_OCR", ".pdf"); apply_unpaper=false,
             ntasks=Sys.CPU_THREADS - 1, tesseract_nthreads=1, pages=num_pages(pdf),
             cleanup_after=true, cleanup_at_exit=true, tmp=get_scratch_dir(pdf),
             verbose=true)

Reads in a PDF located at `pdf`, uses Tesseract to OCR each page and combines the results into a pdf located `output_path`.

Keyword arguments:

* `ntasks`: how many parallel tasks to use for launching `tesseract` and `pdftoppm`.
* `tesseract_nthreads`: how many threads to direct Tesseract to use
* `apply_unpaper`: whether or not to apply `unpaper` to try to improve the image quality
* `tmp`: a directory to store intermediate files. This directory is deleted at the end of the function if `cleanup_after` is set to `true`, and when the Julia session is ended if `cleanup_at_exit` is set to `true`.
* `pages=nothing`: the number of pages of the PDF to process; the default of `nothing` indicates all pages in the PDF. It can help in debugging to set this to something small.
* `verbose`: show a progress bar for each step of the process.

Set `ENV["JULIA_DEBUG"] = SearchablePDFs` to see (many) debug messages.
"""
function ocr(pdf, output_path=string(splitext(pdf)[1], "_OCR", ".pdf"); apply_unpaper=false,
             ntasks=Sys.CPU_THREADS - 1, tesseract_nthreads=1, pages=nothing,
             cleanup_after=true, cleanup_at_exit=true, tmp=get_scratch_dir(pdf),
             verbose=true, force=false, max_files_per_unite=100)
    isfile(pdf) || argument_error("Input file not found at `$pdf`"; exception=true)
    force || require_no_file(output_path; exception=true)
    require_extension(pdf, ".pdf"; exception=true)
    require_extension(output_path, ".pdf"; exception=true)

    total_pages = num_pages(pdf)

    if pages === nothing
        pages = total_pages
    elseif pages > total_pages
        argument_error("`pages` must be less than the total number of pages ($(total_pages))";
                       exception=true)
    end

    mkpath(tmp)
    if cleanup_at_exit
        atexit(() -> rm(tmp; force=true, recursive=true))
    end

    @debug "Found file" pdf pages tmp

    all_logs = Channel{@NamedTuple{page::Union{Int,UnitRange{Int},Missing},binary::String,
                                   stdout::String,stderr::String,code::Int}}(Inf)

    @debug "Generating images..."
    imag_prog = Progress(pages; desc="(1/3) Extracting images: ", enabled=verbose)

    img_paths_grps = asyncmap(Iterators.partition(1:pages, 20); ntasks) do page_range
        paths, pdftoppm_logs = get_images(pdf, page_range, tmp, total_pages)
        put!(all_logs, (; page=page_range, pdftoppm_logs...))
        next!(imag_prog; step=length(page_range))
        return paths
    end

    img_paths = reduce(vcat, img_paths_grps)

    @debug "Finished generating images. Starting tesseracting..."
    ocr_prog = Progress(pages; desc="(2/3) OCRing: ", enabled=verbose)
    pdfs = asyncmap(enumerate(img_paths); ntasks) do (page, img)
        @debug "img" page img
        if apply_unpaper
            img, unpaper_logs = unpaper(img)
            put!(all_logs, (; page, unpaper_logs...))
        end
        pdf, tesseract_logs = make_pdf(img; tesseract_nthreads)
        put!(all_logs, (; page, tesseract_logs...))
        next!(ocr_prog)
        return pdf
    end
    @debug "Finished processing pages. Uniting..."
    unite_dir = joinpath(tmp, "unite")
    unite_progress_meter = Progress(pages + cld(pages, max_files_per_unite) + 1;
                                    desc="(3/3) Collecting pages: ", enabled=verbose)
    unite_many_pdfs!(unite_progress_meter, all_logs, unite_dir, pdfs, output_path;
                     max_files_per_unite)
    @debug "Done uniting pdfs"
    if cleanup_after
        @debug "Cleaning up"
        rm(tmp; recursive=true, force=true)
    end
    @debug "Done"
    if verbose
        isfile(output_path) || @error "File was not generated, check the logs!"
    end
    close(all_logs)
    return (; output_path, logs=collect(all_logs), tmp)
end

#####
##### CLI interface
#####

"""
Create a searchable version of a PDF.
"""
@main function searchable(input_pdf::String,
                          output_path::String=string(splitext(input_pdf)[1], "_OCR",
                                                     ".pdf"); apply_unpaper::Bool=false,
                          ntasks::Int=Sys.CPU_THREADS - 1, tesseract_nthreads::Int=1,
                          keep_intermediates::Bool=false,
                          tmp::String=get_scratch_dir(input_pdf), quiet::Bool=false,
                          logfile::Union{Nothing,String}=nothing, force::Bool=false)
    # some of these are redundant with checks inside `ocr`; that's because we want to do them before the "Starting to ocr" message,
    # and we want them to exit if they fail in a non-interactive context, instead of printing a stacktracee.
    isfile(input_pdf) || argument_error("Input file not found at `$(input_pdf)`")
    force || require_no_file(output_path)
    require_extension(input_pdf, ".pdf")
    require_extension(output_path, ".pdf")
    if logfile !== nothing
        force || require_no_file(logfile)
        require_extension(logfile, ".csv")
    end
    verbose = !quiet
    verbose &&
        println("Starting to ocr `$(input_pdf)`; result will be located at `$(output_path)`.")
    result = ocr(input_pdf, output_path; apply_unpaper, ntasks, tesseract_nthreads,
                 cleanup_after=!keep_intermediates, cleanup_at_exit=!keep_intermediates,
                 tmp, verbose)
    verbose && println("\nOutput is located at `$(output_path)`.")
    if keep_intermediates && verbose
        println("Intermediate files located at `$tmp`.")
    end
    if logfile !== nothing
        verbose && println("Writing logs...")
        CSV.write(logfile, result.logs)
        verbose && println("Logs written to `$(logfile)`.")
    end

    return result
end

end # module


text:
# Julia Process example
# key: cmd, command, external, sh, shell

# test_exit.bat
# @echo off
# echo "Test Exit"
# exit 1

using Logging

bat_file = "text_exit.bat"
cmd = `cmd /c $bat_file`

#
# Method 1
#
println("Method 1")
try
    # IO redirected so we will see error from cmd by default
    run(cmd)
catch e
    #@show e
    proc = e.procs[1]  # Process
    @info "Exit Code = $(proc.exitcode)"
end

#
# Method 2
#
# IO not redirected by default
println("Method 2")
proc = run(pipeline(cmd, stderr=stderr), wait=false)
wait(proc)
@info "Exit Code = $(proc.exitcode)"

#
# Method 2.1 with IO redirect
#
println("Method 2.1")
cerr = IOBuffer()
cout = IOBuffer()
proc = run(pipeline(cmd, stdout=cout, stderr=cerr), wait=false)
wait(proc)
err = String(take!(cerr))
out = String(take!(cout))
@info "stderr:$err"
@info "stdout:$out"
@info "Exit Code = $(proc.exitcode)"


text:
function getkernel(parsed_args)

    afun = if parsed_args["activation"] == "triangle"
        TriangleWave()
    elseif parsed_args["activation"] == "prelu"
        PeriodicReLU()
    elseif parsed_args["activation"] == "sincos"
        SinCosActivation()
    else
        SinActivation()
    end

    if parsed_args["kernel"] == "Matern"
        return Matern(afun,parsed_args["nu"])
    elseif parsed_args["kernel"] == "RelU"
        return RELU()
    elseif parsed_args["kernel"] == "MaternLS"
        return MaternLS(parsed_args["nu"], parsed_args["ell"])
    elseif parsed_args["kernel"] == "RBFLS"
        return RBFLS()
    else
        return RBF(afun)
    end
end

max_stars_repo_name: AaltoML/PeriodicBNN

max_issues_repo_name: AaltoML/PeriodicBNN

max_forks_repo_name: AaltoML/PeriodicBNN


text:
# One dimensional wave equation
using ConservationLawsDiffEq
using OrdinaryDiffEq
const CL = ConservationLawsDiffEq

const CFL = 0.45
const Tend = 1.0
const cc = 1.0

Jf(u::AbstractVector) = [0.0 cc;cc 0.0]
f(u::AbstractVector) = [0.0 cc;cc 0.0]*u
f0(x) = [sin(4*π*x),0.0]

Nflux(ϕl::AbstractVector, ϕr::AbstractVector) = 0.5*(f(ϕl)+f(ϕr))
exact_sol(x, t::Float64) = [0.5*(sin(4*π*(-t+x))+sin(4*π*(t+x))),
0.5*(sin(4*π*(-t+x))-sin(4*π*(t+x)))]

# Now discretizate the domain
mesh = Uniform1DFVMesh(50, [-1.0, 1.0])

function get_problem(alg, mesh)
  #Compute discrete initial data
  u0 = getInitialState(mesh,f0,use_threads = true)

  # Now get a explicit semidiscretization (discrete in space) du_h(t)/dt = f_h(u_h(t))
  f_h = getSemiDiscretization(f,alg,mesh,[Periodic()]; Df = Jf, use_threads = false,numvars = 2)

  #Setup ODE problem for a time interval = [0.0,1.0]
  ode_prob = ODEProblem(f_h,u0,(0.0,Tend))

  #Setup callback in order to fix CFL constant value
  cb = getCFLCallback(f_h, CFL)

  #Estimate an initial dt
  dt = update_dt!(u0, f_h, CFL)
  return ode_prob, cb, dt
end

ode_prob, cb, dt = get_problem(FVSKTScheme(), mesh)
sol = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

ode_prob, cb, dt = get_problem(FVTecnoScheme(Nflux;order=3), mesh)
sol2 = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

ode_prob, cb, dt = get_problem(FVCompWENOScheme(), mesh)
sol3 = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

ode_prob, cb, dt = get_problem(FVCompMWENOScheme(), mesh)
sol4 = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

ode_prob, cb, dt = get_problem(FVSpecMWENOScheme(), mesh)
sol5 = solve(ode_prob,SSPRK22(); dt = dt, callback = cb)

#Estimate errors
u1_h = fv_solution(sol, mesh; vars = 2)
u2_h = fv_solution(sol2, mesh; vars = 2)
u3_h = fv_solution(sol3, mesh; vars = 2)
u4_h = fv_solution(sol4, mesh; vars = 2)
u5_h = fv_solution(sol5, mesh; vars = 2)
get_L1_error(exact_sol,u1_h)
get_L1_error(exact_sol,u2_h)
get_L1_error(exact_sol,u3_h)
get_L1_error(exact_sol,u4_h)
get_L1_error(exact_sol,u5_h)

#Plot
using Plots
plot(u1_h, vars=1, lab="KT y",line = (:dot,2))
plot!(u2_h, vars=1,lab="Tecno y",line=(:dot,3))
plot!(u3_h, vars=1,lab="Comp WENO y",line=(:dot,3))
plot!(u4_h, vars=1,lab="Comp MWENO y",line=(:dot,3))
plot!(u5_h, vars=1,lab="Spec MWENO y",line=(:dot,3))
plot!(CL.cell_centers(mesh), [exact_sol(x,Tend)[1] for x in CL.cell_centers(mesh)],lab="Ref y")


text:
module ArtifactInstallation
using Pkg.Artifacts, Test, Libdl
export do_test

function do_test()
    # First, check that `"c_simple"` is installed automatically
    artifacts_toml = joinpath(@__DIR__, "..", "Artifacts.toml")
    @test artifact_exists(artifact_hash("c_simple", artifacts_toml))

    # Test that we can use `artifact""` to get at c_simple
    c_simple_exe = joinpath(artifact"c_simple", "bin", "c_simple")
    if Sys.iswindows()
        c_simple_exe = "$(c_simple_exe).exe"
    end
    @test isfile(c_simple_exe)

    # Test that we can use a variable, not just a literal:
    c_simple = "c_simple"
    c_simple_exe = joinpath(@artifact_str(c_simple), "bin", "c_simple")
    if Sys.iswindows()
        c_simple_exe = "$(c_simple_exe).exe"
    end
    @test isfile(c_simple_exe)

    # Test that we can dlopen and ccall libc_simple
    libc_simple_path = if Sys.iswindows()
        joinpath(artifact"c_simple", "bin", "libc_simple.dll")
    elseif Sys.isapple()
        joinpath(artifact"c_simple", "lib", "libc_simple.dylib")
    else
        joinpath(artifact"c_simple", "lib", "libc_simple.so")
    end

    libc_simple = dlopen(libc_simple_path)
    @test libc_simple != nothing
    @test ccall(dlsym(libc_simple, :my_add), Cint, (Cint, Cint), 2, 3) == 5
    dlclose(libc_simple)
end

end


text:
ENV["GKSwstype"]="100" #src

#md # # Who am I ?
#md #
#md #  - My name is *Pierre Navaro*
#md #
#md #  - **Fortran 77 + PVM** : during my PhD 1998-2002 (Université du Havre)
#md #
#md #  - **Fortran 90-2003 + OpenMP-MPI** : Engineer in Strasbourg (2003-2015) at IRMA
#md #
#md #  - **Numpy + Cython, R + Rcpp** : Engineer in Rennes (2015-now) at IRMAR
#md #
#md #  - **Julia v1.0** since July 2018
#md #
#md #  ## Instructions to open the notebooks
#md #
#md #  https://github.com/JuliaVlasov/Numkin2019
#md #

#md # ---

#md #
#md # # Why Julia?
#md #
#md # - Born in 2009 and version 1.0 was released in August 2018.
#md # - High-level languages like Python and R let one explore and experiment rapidly, but can run slow.
#md # - Low-level languages like Fortran/C++ tend to take longer to develop, but run fast.
#md # - This is sometimes called the "two language problem" and is something the Julia developers set out to eliminate.
#md # - Julia's promise is to provide a "best of both worlds" experience for programmers who need to develop novel algorithms and bring them into production environments with minimal effort.
#md #

#md # ---

#md # # Julia's Engineering and Design Tradoffs
#md #
#md # - Type structures cannot be changed after being created (less dynamism but memory layout can be optimized for)
#md # - All functions are JIT compiled via LLVM (interactive lags but massive runtime improvements)
#md # - All functions specialize on types of arguments (more compilation but give generic programming structures)
#md # - Julia is interactive (use it like Python and R, but makes it harder to get binaries)
#md # - Julia has great methods for handling mutation (more optimization opportunities like C/Fortran, but more cognitive burden)
#md # - Julia's Base library and most packages are written in Julia, (you can understand the source, but learn a new package)
#md # - Julia has expensive tooling for code generation and metaprogramming (concise and more optimizations, but some codes can be for experienced users)
#md #
#md # To me, this gives me a language with a lot of depth which works well for computationally-expensive scientific applications.
#md #
#md # [© ChrisRackaukas](https://www.youtube.com/watch?v=zJ3R6vOhibA&feature=em-uploademail)

#md # ---

#md # # Type-Dispatch Programming
#md #
#md # - Centered around implementing the generic template of the algorithm not around building representations of data.
#md # - The data type choose how to efficiently implement the algorithm.
#md # - With this feature share and reuse code is very easy
#md #
#md # [JuliaCon 2019 | The Unreasonable Effectiveness of Multiple Dispatch | Stefan Karpinski](https://youtu.be/kc9HwsxE1OY)

#md # ---

using DifferentialEquations, Plots

g = 9.79 # Gravitational constants
L = 1.00 # Length of the pendulum

#Initial Conditions
u₀ = [0, π / 60] # Initial speed and initial angle
tspan = (0.0, 6.3) # time domain

#Define the problem
function simplependulum(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*θ
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-6)
#md nothing # hide

#md # ---

# Analytic solution
u = u₀[2] .* cos.(sqrt(g / L) .* sol.t)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")
plot!(sol.t, u, label = "Analytic")
#md savefig("pendulum1.svg"); nothing # hide

#md # ![](pendulum1.svg)

#md # ---

# [Numbers with Uncertainties](http://tutorials.juliadiffeq.org/html/type_handling/02-uncertainties.html)

using Measurements

g = 9.79 ± 0.02; # Gravitational constants
L = 1.00 ± 0.01; # Length of the pendulum

#Initial Conditions
u₀ = [0 ± 0, π / 60 ± 0.01] # Initial speed and initial angle

#Define the problem
function simplependulum(du, u, p, t)
    θ = u[1]
    dθ = u[2]
    du[1] = dθ
    du[2] = -(g/L)*θ
end

#Pass to solvers
prob = ODEProblem(simplependulum, u₀, tspan)
sol = solve(prob, Tsit5(), reltol = 1e-6);
nothing # hide

#md # ---

# Analytic solution
u = u₀[2] .* cos.(sqrt(g / L) .* sol.t)

plot(sol.t, getindex.(sol.u, 2), label = "Numerical")
plot!(sol.t, u, label = "Analytic")
#md savefig("pendulum2.svg"); nothing # hide

#md # ![](pendulum2.svg)

@test true #src

#md #  Next: [Vlasov-Ampere with FFT](/Numkin2019/02/build/index.html)


text:
export unfold, fold

"""
    unfold(X, W; stride=1, padding=0, dilation=1)
Extracts sliding local blocks from a batched input tensor. X is the input 5d vector of size
`(spatial_dims... , channels, batch_size)`. W is the size of kernel, in format
`(spatial_dims... , channels)`. Output has the size of `(L, channels*kernel_w*kernel_h*kernel_d, batch_size)`,
where L is the total number of blocks.

"""
function unfold(X::AbstractArray{T,M} where T, w_dim::NTuple{K}; stride=1, padding=0, dilation=1) where M where K
    x_dim = size(X)
    if ndims(X) > 5 || ndims(X) < 3
        throw(DimensionMismatch("X and W must be 3d/4d/5d for 1d/2d/3d image. got $(ndims(X))d input"))
    end

    if ndims(X)-2 != length(w_dim)-1
        throw(DimensionMismatch("spatial dimentions of image and kernel must be equal, got $(ndims(X)-2),$(length(w_dim)-1)"))
    end

    # reassign x_dim after converting it to a 3d image type input
    x_dim = ( x_dim[1:end-2]... , fill(1,5-ndims(X))... , x_dim[end-1:end]... )
    # w_dim must be in following format: (spatial_dims..., channels_in, channels_out)
    w_dim = ( w_dim[1:end-1]... , fill(1,4-length(w_dim))... , w_dim[end], w_dim[end] )
    X = reshape(X, x_dim)

    # Make DenseConvDims object
    cdims = DenseConvDims(x_dim, w_dim; stride=stride, padding=padding, dilation=dilation)

    # Calculate the total number of sliding blocks
    col_dim = (im2col_dims(cdims))[1:2] # im2col_dims() returns (col_dim_x, col_dim_y, thread_num)
    col = fill(0., col_dim[1],col_dim[2],x_dim[end]) # x_dim[end] is number of batches

    # Iterate through all batchs
    @views for i = 1:x_dim[end]
        im2col!(col[:,:,i], X[:,:,:,:,i], cdims)
    end
    return col
end

"""
    fold(col, out_dim, W; stride=1, padding=0, dilation=1)
Does the opposite of `unfold()`, Combines an array of sliding local blocks into a large containing
tensor. `col` is a 3d array of shape `(L, channels*kernel_w*kernel_h*kernel_d, batch_size)`, where,
L is the total number of blocks. out_dim is the spatial dimention of the required image. W is the
spatial dimentions of the kernel.

"""
function fold(col::AbstractArray{T,3} where T, out_dim::NTuple{M}, w_dim::NTuple{M}; stride=1, padding=0, dilation=1) where M
    # Validate input
    if length(out_dim) > 3
        throw(DimensionMismatch("output dimentions cannot be greater than 3, got $(ndims(out_dim))"))
    end

    # Create DenseConvDims object
    col_dim = size(col)
    channels = col_dim[2]÷prod(w_dim)
    x_dim = (out_dim... , fill(3-length(out_dim))... , channels,col_dim[3])
    w_dim = (w_dim... , fill(3-length(w_dim))... , channels,channels)
    cdims = DenseConvDims(x_dim,w_dim; stride=stride, padding=padding, dilation=dilation)

    img = fill(0., x_dim)

    # Iterate through all batchs
    @views for i = 1:x_dim[end]
        col2im!(img[:,:,:,:,i], col[:,:,i], cdims)
    end

    return reshape(img, (out_dim... , channels,col_dim[3]))
end


text:
#TODO:
  # schema resolution
  # parsing canonical form
"""
The Avro.jl package provides a pure Julia implementation for reading
writing data in the [avro format](http://avro.apache.org/docs/current/spec.html).

### Implementation status

It currently supports:

  * All primitive types
  * All nested/complex types
  * Logical types listed in the spec (Decimal, UUID, Date, Time, Timestamps, Duration)
  * Binary encoding/decoding
  * Reading/writing object container files via the Tables.jl interface
  * Supports the xz, zstd, deflate, and bzip2 compression codecs for object container files

Currently not supported are:

  * JSON encoding/decoding of objects
  * Single object encoding or schema fingerprints
  * Schema resolution
  * Protocol messages, calls, handshakes
  * Snappy compression

### Package motivation

Why use the avro format vs. other data formats? Some benefits include:
  * Very concise binary encoding, especially object container files with compression
  * Very fast reading/writing
  * Objects/data must have well-defined schema
  * One of the few "row-oriented" binary data formats

### Getting started

The Avro.jl package supports two main APIs to interact with avro data.
The first is similar to the JSON3.jl struct API for interacting with json data,
largely in part due to the similarities between the avro format and json.
This looks like:

```julia
buf = Avro.write(obj)
obj = Avro.read(buf, typeof(obj))
```

In short, we use [`Avro.write`](@ref) and provide an object `obj` to write
out in the avro format. We can optionally provide a filename or `IO` as a
first argument to write the data to.

We can then read the data back in using [`Avro.read`](@ref), where the first argument
must be a filename, `IO`, or any `AbstractVector{UInt8}` byte buffer containing
avro data. The 2nd argument is _required_, and is the type of data to be read.
This type can be provided as a simple Julia type (like `Avro.read(buf, Vector{String})`),
or as a parsed avro schema, like `Avro.read(buf, Avro.parseschema("schema.avsc"))`.
[`Avro.parseschema`](@ref) takes a filename or json string representing the avro schema of the
data to read and returns a "schema type" that can be passed to `Avro.read`.

The second alternative API allows "packaging" the data's schema with the data
in what the avro spec calls an "object container" file. While `Avro.read`/`Avro.write`
require the user to already know or pass the schema externally, [`Avro.writetable`](@ref)
and [`Avro.readtable`](@ref) can write/read avro object container files, and will take care
of any schema writing/reading, compression, etc. automatically. These table
functions unsurprisingly utilize the ubiquitous Tables.jl interface to facilitate
integrations with other formats.

```julia
# write our input_table out to a file named "data.avro" using the zstd compression codec
# input_table can be any Tables.jl-compatible source, like CSV.File, Arrow.Table, DataFrame, etc.
Avro.writetable("data.avro", input_table; compress=:zstd)

# we can also read avro data from object container files
# if file uses compression, it will be decompressed automatically
# the schema of the data is packaged in the object container file itself
# and will be parsed before constructing the file table
tbl = Avro.readtable("data.avro")
# the returned type is `Avro.Table`, which satisfies the Tables.jl interface
# which means it can be sent to any valid sink function, like
# Arrow.write("data.arrow", tbl), CSV.write("data.csv", tbl), or DataFrame(tbl)
```

"""
module Avro

using Mmap, UUIDs, Dates
using StructTypes, JSON3, Tables, SentinelArrays
using CodecBzip2, CodecXz, CodecZlib, CodecZstd#, CodecSnappy

import Base.open
import Base.close
import Base.iterate
import Base: ==

const DEBUG_LEVEL = Ref(0)

function setdebug!(level::Int)
    DEBUG_LEVEL[] = level
    return
end

function withdebug(f, level)
    lvl = DEBUG_LEVEL[]
    try
        setdebug!(level)
        f()
    finally
        setdebug!(lvl)
    end
end

macro debug(level, msg)
    esc(quote
        if DEBUG_LEVEL[] >= $level
            println(Base.string("DEBUG: ", $(QuoteNode(__source__.file)), ":", $(QuoteNode(__source__.line)), " ", $msg))
        end
    end)
end

include("utils.jl")
include("types.jl")
include("tables.jl")
include("rowwise.jl")

const COMPRESSORS = Dict{Symbol, Any}(
    :bzip2 => Bzip2Compressor[],
    :xz => XzCompressor[],
    :deflate => DeflateCompressor[],
    :zstd => ZstdCompressor[],
    # :snappy => SnappyCompressor[],
)

const DECOMPRESSORS = Dict{Symbol, Any}(
    :bzip2 => Bzip2Decompressor,
    :xz => XzDecompressor,
    :deflate => DeflateDecompressor,
    :zstd => ZstdDecompressor,
    # :snappy => SnappyDecompressor,
)

function __init__()
    for _ = 1:1 # replace 1 w/ Threads.nthreads() if we do multithreaded writing
        bzip2 = Bzip2Compressor()
        CodecBzip2.TranscodingStreams.initialize(bzip2)
        push!(COMPRESSORS[:bzip2], bzip2)

        xz = XzCompressor()
        CodecXz.TranscodingStreams.initialize(xz)
        push!(COMPRESSORS[:xz], xz)

        deflate = DeflateCompressor()
        CodecZlib.TranscodingStreams.initialize(deflate)
        push!(COMPRESSORS[:deflate], deflate)

        zstd = ZstdCompressor(; level=3)
        CodecZstd.TranscodingStreams.initialize(zstd)
        push!(COMPRESSORS[:zstd], zstd)

        # snappy = SnappyCompressor()
        # CodecSnappy.TranscodingStreams.initialize(snappy)
        # push!(COMPRESSORS[:snappy], snappy)
    end
    return
end

end


text:
using AutoHashEquals
using NamedArrays

import DataStructures
import Combinatorics
import Base.show
import Base.getindex
import Base.copy
import Base.==
import Base.isless
import Base.hash
import StatsBase

export AbstractTech, Tech, TechCombined, Product, System
export build_all_systems
export add_loop_techs!
export writedotfile
export prefilterTechList

# -----------
# define types

@auto_hash_equals struct Product
name::Symbol
end

Product(name::String) = Product(Symbol(name))
show(io::Base.IO, p::Product) = print(io, "$(p.name)")

isless(p1::Product, p2::Product) = isless(p1.name, p2.name)

abstract type AbstractTech end

@auto_hash_equals struct Tech <: AbstractTech
    inputs::Array{Product}
    outputs::Array{Product}
    name::String
    functional_group::Symbol
    appscore::Array{Float64}
    n_inputs::Int
    transC::Dict{String, Dict{Product, Float64}}
    transC_reliability::Dict{String, Float64}
end

const Source = Tech
const Sink = Tech

function Tech(inputs::Array{Product}, outputs::Array{Product},
              name::String, functional_group::Symbol,
              appscore::Float64, n_inputs::Int,
              transC::Dict{String, Dict{Product, Float64}},
              transC_reliability::Dict{String, Float64})
    Tech(inputs, outputs,
         name, functional_group,
         Float64[appscore],
         n_inputs,
         transC, transC_reliability)
end

"""
The `Tech` type represents Technolgies.
It consist of `inputs`, `outputs`, a `name`, a `functional_group`, and a transfer coefficients `transC`.
"""
function Tech{T<:String}(inputs::Array{T}, outputs::Array{T}, name::T, functional_group::T,
                         appscore::Float64,
                         transC::Dict{String, Dict{Product, Float64}},
                         transC_reliability::Dict{String, Float64})

    # sanity check
    for sub in keys(transC)
        s = 0.0
        for val = values(transC[sub])
            s += val
        end
        if !isapprox(s, 1.0)
            error("Tech $(name): The transfer coefficients for $sub do not sum to 1!")
        end
    end

    Tech([Product(x) for x in inputs],
         [Product(x) for x in outputs],
         name,
         Symbol(functional_group),
         appscore,
         size(inputs,1),
         transC,
         transC_reliability)
end

# Functions for pretty printing
function show(io::Base.IO, t::AbstractTech)
    instr = ["$(ii.name)" for ii in t.inputs]
    outstr = ["$(ii.name)" for ii in t.outputs]
    instr = length(instr) >0 ? instr : ["Source"]
    outstr = length(outstr) >0 ? outstr : ["Sink"]
    print(io, "$(t.name): ($(join(instr, ", "))) -> ($(join(outstr, ", ")))")
end

"""
A `Connection`is a Tuples{Product, sourceTech, sinkTech}.
"""
const Connection = Tuple{Product, AbstractTech, AbstractTech}

"""
Type of combined (looped) Techs.
"""
@auto_hash_equals struct TechCombined <: AbstractTech
    inputs::Array{Product}
    outputs::Array{Product}
    name::String
    functional_group::Symbol
    appscore::Array{Float64}
    n_inputs::Int
    internal_techs::Set{Tech}
    internal_connections::Set{Connection}
end

"""
A `System` consists of `techs` and `conncetions`.
"""
mutable struct System
    techs::Set{AbstractTech}
    connections::Set{Connection}
    complete::Bool
    properties::Dict
end

# ignore properties for isequal
Base.hash(sys::System, h::UInt) = hash(sys.techs, hash(sys.connections, hash(sys.complete, hash(:System, h))))
Base.:(==)(a::System, b::System) = isequal(a.techs, b.techs) && isequal(a.connections, b.connections) && isequal(a.complete, b.complete) && true

System{T <: AbstractTech}(techs::Set{T}, con::Set{Connection}, complete::Bool) = System(techs, con, complete, Dict())
System{T <: AbstractTech}(techs::Array{T}, con::Array{Connection}, complete::Bool) = System(Set(techs), Set(con), complete)
System{T <: AbstractTech}(techs::Array{T}, con::Array{Connection}) = System(Set(techs), Set(con), false)
System{T <: AbstractTech}(techs::Array{T}) = System(Set(techs), Set(Connection[]), false)

# Function to copy a System
function copy(sys::System)
    System(copy(sys.techs), copy(sys.connections), copy(sys.complete))
end

function show(io::Base.IO, s::System)
    !s.complete ? print(io, "Incomplete ") : nothing
    println(io, "System with $(length(s.techs)) technologies and $(length(s.connections)) connections: ")
    for i in s.connections
        println(io, " $(i[2].name) → ($(i[1].name)) → $(i[3].name)")
    end
end

# -----------
# helper functions

function get_outputs{T2<:Union{Array{T1}, Set{T1}} where T1<:AbstractTech}(techs::T2)
    outs = Product[]
    for t in techs
        append!(outs, t.outputs)
    end
    return outs
end

"""
returns all "open" outputs of a system

A dic
"""
function get_outputs(sys::System)
    # all outs
    outs = get_outputs(sys.techs)

    for c in sys.connections
        idx = find(outs .== c[1]) # get index of appl products mathing c
        deleteat!(outs, idx[1])   # remove one product
    end
    return outs
end

function get_inputs{T2<:Union{Array{T1}, Set{T1}} where T1<:AbstractTech}(techs::T2)
    ins = Product[]
    for t in techs
        append!(ins, t.inputs)
    end
    return ins
end

"""
return all "open" inputs of a system
"""
function get_inputs(sys::System)
    # all ins
    ins = DataStructures.counter(get_inputs(sys.techs))
    for c in sys.connections
        num = push!(ins, c[1], -1) # not very elegant...
        if num <= 0
            DataStructures.reset!(ins, c[1])
        end
    end
    return collect(keys(ins))
end

function is_complete(sys::System)
    length(get_outputs(sys)) == 0 &&
        length(get_inputs(sys)) == 0 # &&
    # any(length(t.outputs) == 0 for t in sys.techs) # check if at least one sink exist
end

"""
Return all technologies of the system that have an open `prod` output
"""
function get_openout_techs(sys::System, prod::Product)

    function is_connected(tech, sys)
        for c in sys.connections
            if c[2] == tech && c[1] == prod
                return true
            end
        end
        return false
    end

    matching_techs = filter(t -> prod in t.outputs, sys.techs) # Techs with matching outputs			# SLOW
    filter(t -> !is_connected(t, sys), matching_techs) # Techs open outputs				        # SLOW
end

"""
Return all technologies of the system that have an open `prod` input
"""
function get_openin_techs(sys::System, prod::Product)

    function is_connected(tech, sys)
        for c in sys.connections
            if c[3] == tech && c[1] == prod
                return true
            end
        end
        return false
    end

    matching_techs = filter(t -> prod in t.inputs, sys.techs) # Techs with matching inputs
    # filter(t -> !is_connected(t, sys), matching_techs) # Techs open inputs
    matching_techs
end

# pre filter the tech list
function get_candidates{T <: AbstractTech}(s::Array{T}, outs, k)
    # Limitation!
    # This function guaranties, that we never add a techs, which leads to new open inputs.
    # However, this banns not only loops, but also "triangles". E.g if only "A is given"
    #   A ---> B
    #    \    ^
    #     \  /
    #      >C
    # this system cannot be found! Because B, C are on the same stage
    # of the algorithm. (Exception: both inputs of B have the same product)
    n_out = length(outs)
    n_in_max = n_out - k + 1
    function condi(t::AbstractTech)
        issubset(t.inputs, outs) && (t.n_inputs <= n_in_max)

    end

    filter(condi, s)
end

# -----------
# functions to find all systems

# Return a vector of Systems
function build_system!{T <: AbstractTech}(sys::System, completesystems::Array{System},
                                          techs::Array{T}, islegal::Function,
                                          logfile::IO, hashset::Set{UInt64})

    # get Array of matching Techs Arrays
    candidates = get_candidates(sys, techs)

    for candidate in candidates
        # extend systems
        sys_exts = extend_system(sys, candidate)
        for sys_ext in sys_exts
            if sys_ext.complete && !(sys_ext in completesystems)
                push!(completesystems, sys_ext)
                println(logfile, now())
                println(logfile, sys_ext)
                flush(logfile)
            elseif !sys_ext.complete && islegal(sys_ext) && !(hash(sys_ext) in hashset)
                push!(hashset, hash(sys_ext))
                build_system!(sys_ext, completesystems,
                              techs, islegal, logfile, hashset)
            end
        end
    end
end

"""
Returns an Array of all possible `System`s starting with `source`. A source can be any technology with a least one output.
"""
function build_all_systems{T1 <: AbstractTech,
                           T2 <: AbstractTech}(source::Array{T1}, techs::Array{T2};
                                               islegal::Function=x -> true,
                                               logfile::IO=STDOUT)

    completesystems = System[]
    build_system!(System(source), completesystems, techs, islegal,
                  logfile, Set{UInt64}())

    # split TechCombineds
    split_techcombined!.(completesystems)
    completesystems = unique(completesystems) # before split techs system are different but once splitted they are identical (e.g. due to permuations of some techs in two combitechs)
    return completesystems
end

# replace combined techs with seperate internal techs
function split_techcombined!(sys)

    combtechs = filter(t -> typeof(t) == TechCombined, sys.techs)

    # remove TechCombineds
    filter!(t -> typeof(t) != TechCombined, sys.techs)

    # remove connection to/from TechCombineds
    allConn_TechCombs = filter(c -> any(contains.([c[2].name, c[3].name], "::")),
                               sys.connections)
    filter!(c -> !(c in allConn_TechCombs), sys.connections)

    # add internal techs/connections
    for t in combtechs
            union!(sys.techs, t.internal_techs)
            union!(sys.connections, t.internal_connections)
    end

    # add connection from/to internal techs
    for c in allConn_TechCombs
        prod, from_tech, to_tech = c
        if contains(from_tech.name, "::")

            int_con_prod = collect(filter(c -> c[1] == prod, from_tech.internal_connections))
            ffilter = function (t::AbstractTech)
                prod in t.outputs &&
                    !any(getindex.(int_con_prod, 2) == t)
            end

            from_tech = collect(filter(ffilter, from_tech.internal_techs))
        end
        if contains(to_tech.name, "::")

            int_con_prod = collect(filter(c -> c[1] == prod, to_tech.internal_connections))
            ffilter = function (t::AbstractTech)
                prod in t.inputs &&
                    !any(getindex.(int_con_prod, 3) == t)
            end

            to_tech = collect(filter(ffilter, to_tech.internal_techs))
        end

       if ! (typeof(from_tech) <: AbstractArray)
           from_tech = [from_tech]
       end
       if ! (typeof(to_tech) <: AbstractArray)
           to_tech = [to_tech]
       end
        for ft in [from_tech...]
            for tt in [to_tech...]
               push!(sys.connections, (prod, ft, tt))
            end
        end
    end
end

# converts a vecptr of TEchs and techsCombined to a purly Tech array
function meltTechs{T <: AbstractTech}(techs::Union{Array{T}, Set{T}})
    flattechs = Tech[]
    for t in techs
        if typeof(t)==Tech
            push!(flattechs, t)
        else
            append!(flattechs, t.internal_techs)
        end
    end
    return flattechs
end

# Returns techs that fit to an open system
function get_candidates{T <: AbstractTech}(sys::System, techs::Array{T})

    # filter out already used techs (also in TechCombined)
    systechs = meltTechs(sys.techs)
    ffilter(t::Tech) = !(t in systechs)
    ffilter(t::TechCombined) = length(intersect(t.internal_techs, systechs)) == 0
    techssub = filter(t -> ffilter(t), techs)

    outs = get_outputs(sys)

    matching_techs = Array{Array{AbstractTech},1}()
    n_out = length(outs)

    ## get matching combinations
    for k in 1:n_out
        for c in Combinatorics.combinations(get_candidates(techssub, outs, k), k) # try all combination of lenght k (in R: combn())
            inputs = get_inputs(c)
            if is_compatible(outs, inputs)
                push!(matching_techs, c)
            end

        end
    end

    ## check that no matching combination has dublicated Techs (e.g. not ok: [A, A::B])
    function filterdublicates(m)
        m_flat = meltTechs(m)
        length(m_flat) == length(unique(m_flat))
    end
    filter!(filterdublicates, matching_techs)

    return matching_techs
end

# test if outputs are compatible with next inputs.
function is_compatible(outputs, inputs)
    outputs = DataStructures.counter(outputs)
    inputs = DataStructures.counter(inputs)

    # check if outputs and inputs are compatible
    if length(setdiff(keys(outputs), keys(inputs))) != 0
        return false
    end

    # exclude spliting of one output to multiple inputs
    for prod in keys(outputs)
        if outputs[prod] < inputs[prod]
            return false
        end
    end
    return true
end

#-----------------------------------------------------
# helper for extend_system

"""
return all variation how sys can be extended with tech_comb
"""
function make_connections!(part_sys::Array{System}, extended_sys::Array{System})
    for s in part_sys
        if length(get_inputs(s)) == 0 # check if more connections must be added
            if is_complete(s) # check if the whole system is complete
                s.complete = true
            end
            push!(extended_sys, s)
        else
            part_sys_new = add_next_connection(s)
            make_connections!(part_sys_new, extended_sys)
        end
    end
end

"""
helper to add connections for the next product
"""
function add_next_connection(sys::System) # sys is the combination of the initial sys and tech_comb
    new_part_sys = System[]

    # select first open input product
    p = get_inputs(sys)[1]
    in_techs = get_openin_techs(sys, p)
    out_techs = get_openout_techs(sys, p)
    setdiff!(out_techs, in_techs) # to avoid looping back to the same tech
    # has additional benefit of finding loops between techs of tech_comb in case
    # of an output of one tech from tech_comb is part of the set of "open inputs"
    n = length(out_techs)
    if n > 0
        for connection_combs in Base.product(Base.Iterators.repeated(in_techs, n)...)
            if length(unique(connection_combs)) == length(in_techs) # no empty input techs allowed / all in_techs must be used
                sysi = copy(sys)
                for (i, out_tech) in enumerate(out_techs)
                    push!(sysi.connections, (p, out_tech, connection_combs[i])) # add new connection between out_techs and the combinations of in_techs for product p
                end
                push!(new_part_sys, sysi)
            end
        end
    end
    return new_part_sys
end

"""
Return an array of all possible extension of `sys` with the candidate technology
"""
function extend_system{T <: AbstractTech}(sys::System, tech_comb::Array{T})

    sysi = copy(sys)
    union!(sysi.techs, tech_comb)
    extended_sys = System[]
    make_connections!([sysi], extended_sys)

    return extended_sys
end

"""
Return an array possible Technologies (subset of techlist) for the given Sources.
The number of technologies is reduced by removing Techs that require an input that is not available with
the sources provided.
"""
function prefilterTechList{T1 <: AbstractTech, T2 <: AbstractTech}(currentSources::Array{T1}, sources::Array{T2},
                           sourcesAdd::Array{T2}, tech_list::Array{T2})

    # All Products that can be created by available sources
    otherSourcesProduct = vcat([t.outputs for t in sources]...)
    otherSourcesAddProduct = vcat([t.outputs for t in sourcesAdd]...)
    append!(otherSourcesProduct, otherSourcesAddProduct)

    # Outputs of all the used Products by current Source
    output_list = Product[]
    for ss_c in currentSources
        append!(output_list, ss_c.outputs)
    end

    otherSourcesProduct = filter(x -> !(x in output_list), otherSourcesProduct)
    otherSourcesProduct = map(x -> "$(x.name)", otherSourcesProduct) # convert to Strings
    push!(otherSourcesProduct, "^[.]") # this is a dummy patter to ensure that filter works if otherSourcesProduct is empty

    # # Additional Filter. REMOVE HARD CODING OF PRODUCT NAMES
    # for tproduct in otherSourcesProduct
    #     if tproduct == "excreta"
    #         append!(otherSourcesProduct, ["pithumus"])
    #     end
    # end

    # check of any the String in otherSourcesProduct is part of an input products name
    function ffilter(x)
        inputs = map(x -> "$(x.name)", x.inputs)
        match_inputs = filter(Regex(join(otherSourcesProduct, '|')), inputs)
        length(match_inputs) == 0
    end

    sub_tech_list = filter(ffilter, tech_list)

    return sub_tech_list
end

"""
add tech combinations with internal loops. Only for techs in group :S and :T
"""
function add_loop_techs!(tech_list::Array{AbstractTech}; groups = [:S, :T])
    for fgroup in groups
        tech_list_group = filter(t -> t.functional_group == fgroup, tech_list)
        for tech1 in tech_list_group
            for tech2 in tech_list_group
                if tech1 != tech2

                    outs1 = tech1.outputs
                    outs2 = tech2.outputs
                    ins1 = tech1.inputs
                    ins2 = tech2.inputs
                    exchange_1_2 = intersect(outs1, ins2)
                    exchange_2_1 = intersect(outs2, ins1)
                    if (length(exchange_1_2) >= 1 | length(exchange_2_1) >= 1) && length(intersect(exchange_1_2, exchange_2_1))==0 &&
                        length(intersect(ins1, ins2)) == 0
                        # matching tech are partners!
                        tech_news = make_looped_techs(tech1, tech2)
                        for tech_new in tech_news
                            if !(tech_new in tech_list) && length(tech_new.outputs)>0 && length(tech_new.inputs)>0
                                push!(tech_list, tech_new)
                            end
                        end
                    end
                end
            end
        end
    end
end

function make_looped_techs(tech1::Tech, tech2::Tech)

    # input_new = union(in1, in2) - intersect(out1, in2) - intersect(out2, in1)
    # output_new = union(out1, out2) - intersect(out1, in2) - intersect(out2, in1)

    inputs = DataStructures.counter(vcat(tech1.inputs, tech2.inputs))
    outputs = DataStructures.counter(vcat(tech1.outputs, tech2.outputs))

    internal_12 = DataStructures.counter(vcat(intersect(tech1.outputs, tech2.inputs)))
    internal_21 = DataStructures.counter(vcat(intersect(tech2.outputs, tech1.inputs)))

    internal_products = [collect(keys(internal_12)), collect(keys(internal_21))]

    # --- remove internal connections for new tech ins and outs
    # reduce counter
    for p in keys(merge(internal_12, internal_21))
        push!(inputs, p, -1)
        push!(outputs, p, -1)
    end
    # remove all products that occure 0 times
    ins = Product[p[1] for p in inputs if p[2]>0]
    outs = Product[p[1] for p in outputs if p[2]>0]

    # --- internal connections

    internal_connections = Set{Connection}()
    for p in internal_12
        if p[2]>0
            push!(internal_connections, (p[1], tech1, tech2))
        end
    end
    for p in internal_21
        if p[2]>0
            push!(internal_connections, (p[1], tech2, tech1))
        end
    end

    # name and apps score
    name = join(sort([tech1.name, tech2.name]), " :: ")
    appscore = sort(vcat(tech1.appscore, tech2.appscore))

    # make input variations
    newtechs = TechCombined[]
    tech = TechCombined(sort(ins), sort(outs), name, tech1.functional_group,
                        appscore, length(ins),
                        Set([tech1, tech2]), internal_connections)
    push!(newtechs, tech)
    for additional_inputs = Combinatorics.combinations(internal_products)
         ins2 = vcat([ins..., additional_inputs...]...)
        tech = TechCombined(sort(ins2), sort(outs), name, tech1.functional_group,
                             appscore, length(ins),
                             Set([tech1, tech2]), internal_connections)
         push!(newtechs, tech)
    end

    return newtechs
end


text:
using NLsolve

poly(x, p) = sum(p[n] * x^(n-1) for n in 1:length(p))
∂poly(x, p) = sum(n * p[n+1] * x^(n-1) for n in 1:length(p)-1)

"""
        fit_cubic(p1, p2, s1, s2)

Return the coeficients c such that f(x) = c₁ + c₂x + c₃x² + c₄x³ is a cubic polynomial that passes through the points p1 = (x1, y1) and p2 = (x2, y2) with slope s1 at p1 and slope s2 at p2.
"""
function fit_cubic(p1, p2, s1, s2)
    x1, y1 = p1
    x2, y2 = p2

    function f!(F, c)
        F[1] = poly(x1, c) - y1
        F[2] = poly(x2, c) - y2
        F[3] = ∂poly(x1, c) - s1
        F[4] = ∂poly(x2, c) - s2
    end

    results = nlsolve(f!, zeros(4))
    return results.zero
end


text:
# wrappers for functionality provided by the CUDA toolkit

# extensions to the C language
include("cuda/memory_shared.jl")
include("cuda/indexing.jl")
include("cuda/synchronization.jl")
include("cuda/warp_vote.jl")
include("cuda/warp_shuffle.jl")
include("cuda/output.jl")
include("cuda/assertion.jl")
include("cuda/memory_dynamic.jl")
include("cuda/atomics.jl")
include("cuda/misc.jl")
if VERSION >= v"1.5-"
include("cuda/wmma.jl")
end

# functionality from libdevice
#
# > The libdevice library is a collection of NVVM bitcode functions that implement common
# > functions for NVIDIA GPU devices, including math primitives and bit-manipulation
# > functions. These functions are optimized for particular GPU architectures, and are
# > intended to be linked with an NVVM IR module during compilation to PTX.
include("cuda/math.jl")
# TODO: native mathematical functions, CUDA C programming guide" > "C language extensions"
#       https://docs.nvidia.com/cuda/cuda-math-api/group__CUDA__MATH__INTRINSIC__DOUBLE.html
#       see /path/to/cuda/include/sm_20_intrinsics.h

# functionality from libcudadevrt
#
# The libcudadevrt library is a collection of PTX bitcode functions that implement parts of
# the CUDA API for execution on the device, such as device synchronization primitives,
# dynamic kernel APIs, etc.
using CEnum
include("cuda/libcudadevrt_common.jl")
include("cuda/libcudadevrt.jl")
include("cuda/cooperative_groups.jl")
include("cuda/dynamic_parallelism.jl")


text:
##### Sum wrapper

export Σfluxes, Σsources

"""
    flux

An individual flux.
See [`BalanceLaw`](@ref) for more info.
"""
function flux end

"""
    source

An individual source.
See [`BalanceLaw`](@ref) for more info.
"""
function source end

"""
    ntuple_sum(nt::NTuple{N,T}) where {N, T}

sum of `NTuple`, which requires more strict
type input than `sum`. This is added to better
synchronize the success/failure between CPU/GPU
runs to help improve debugging.
"""
ntuple_sum(nt::NTuple{N, T}) where {N, T} = sum(nt)

"""
    Σfluxes(fluxes::NTuple, bl, args)

Sum of the fluxes where
 - `fluxes` is an `NTuple{N, TendencyDef{Flux{O}}} where {N, O}`
 - `bl` is the balance law
 - `args` are the arguments passed to the individual `flux` functions
"""
function Σfluxes(
    pv::PV,
    fluxes::NTuple{N, TendencyDef{Flux{O}}},
    bl,
    args,
) where {N, O, PV}
    return ntuple_sum(
        ntuple(Val(N)) do i
            projection(pv, bl, fluxes[i], args, flux(pv, fluxes[i], bl, args))
        end,
    )
end

# Emptry scalar case:
function Σfluxes(
    pv::PV,
    fluxes::NTuple{0, TendencyDef{Flux{O}}},
    args...,
) where {O, PV}
    return SVector(0, 0, 0)
end

# Emptry vector case:
function Σfluxes(
    pv::PV,
    fluxes::NTuple{0, TendencyDef{Flux{O}}},
    args...,
) where {O, PV <: AbstractMomentumVariable}
    return SArray{Tuple{3, 3}}(ntuple(i -> 0, 9))
end

# Emptry tracer case:
function Σfluxes(
    pv::PV,
    fluxes::NTuple{0, TendencyDef{Flux{O}}},
    args...,
) where {O, N, PV <: AbstractTracersVariable{N}}
    return SArray{Tuple{3, N}}(ntuple(i -> 0, 3 * N))
end

"""
    Σsources(sources::NTuple, bl, args)

Sum of the sources where
 - `sources` is an `NTuple{N, TendencyDef{Source}} where {N}`
 - `bl` is the balance law
 - `args` are the arguments passed to the individual `source` functions
"""
function Σsources(
    pv::PV,
    sources::NTuple{N, TendencyDef{Source}},
    bl,
    args,
) where {N, PV}
    return ntuple_sum(
        ntuple(Val(N)) do i
            projection(
                pv,
                bl,
                sources[i],
                args,
                source(pv, sources[i], bl, args),
            )
        end,
    )
end

# Emptry scalar case:
function Σsources(
    pv::PV,
    sources::NTuple{0, TendencyDef{Source}},
    args...,
) where {PV}
    return 0
end

# Emptry vector case:
function Σsources(
    pv::AbstractMomentumVariable,
    sources::NTuple{0, TendencyDef{Source}},
    args...,
)
    return SVector(0, 0, 0)
end

# Emptry tracer case:
function Σsources(
    pv::AbstractTracersVariable{N},
    sources::NTuple{0, TendencyDef{Source}},
    args...,
) where {N}
    return SArray{Tuple{N}}(ntuple(i -> 0, N))
end


text:
# ---
# title: 951. Flip Equivalent Binary Trees
# id: problem951
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Tree
# link: <https://leetcode.com/problems/flip-equivalent-binary-trees/description/>
# hidden: true
# ---
#
# For a binary tree **T** , we can define a **flip operation** as follows:
# choose any node, and swap the left and right child subtrees.
#
# A binary tree **X**  is _flip equivalent_ to a binary tree **Y** if and only
# if we can make **X** equal to **Y** after some number of flip operations.
#
# Given the roots of two binary trees `root1` and `root2`, return `true` if the
# two trees are flip equivelent or `false` otherwise.
#
#
#
# **Example 1:**
#
# ![Flipped Trees
# Diagram](https://assets.leetcode.com/uploads/2018/11/29/tree_ex.png)
#
#
#
#     Input: root1 = [1,2,3,4,5,6,null,null,null,7,8], root2 = [1,3,2,null,6,4,5,null,null,null,null,8,7]
#     Output: true
#     Explanation: We flipped at nodes with values 1, 3, and 5.
#
#
# **Example 2:**
#
#
#
#     Input: root1 = [], root2 = []
#     Output: true
#
#
# **Example 3:**
#
#
#
#     Input: root1 = [], root2 = [1]
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: root1 = [0,null,1], root2 = []
#     Output: false
#
#
# **Example 5:**
#
#
#
#     Input: root1 = [0,null,1], root2 = [0,1]
#     Output: true
#
#
#
#
# **Constraints:**
#
#   * The number of nodes in each tree is in the range `[0, 100]`.
#   * Each tree will have **unique node values** in the range `[0, 99]`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using DiffEqGPU, StochasticDiffEq, Test, DiffEqNoiseProcess
using CUDA
using Random

function brusselator_f!(du,u,p,t)
 @inbounds begin
     du[1] = (p[1]-1)*u[1]+p[1]*u[1]^2+(u[1]+1)^2*u[2]
     du[2] = -p[1]*u[1]-p[1]*u[1]^2-(u[1]+1)^2*u[2]
 end
 nothing
end

function scalar_noise!(du,u,p,t)
 @inbounds begin
     du[1] = p[2]*u[1]*(1+u[1])
     du[2] = -p[2]*u[1]*(1+u[1])
 end
 nothing
end

function prob_func(prob, i, repeat)
    Random.seed!(seeds[i])
    W = WienerProcess(0.0,0.0,0.0)
    remake(prob,noise=W)
end

# fix seeds
seed = 100
Random.seed!(seed)
numtraj= 100
seeds = rand(UInt, numtraj)
W = WienerProcess(0.0,0.0,0.0)

CUDA.allowscalar(false)
u0 = [-0.1f0,0.0f0]
tspan = (0.0f0,100.0f0)
p = [1.9f0,0.1f0]

prob = SDEProblem(brusselator_f!,scalar_noise!,u0,tspan,p, noise=W)
ensembleprob = EnsembleProblem(prob, prob_func = prob_func)

@info "Brusselator"

#Performance check with nvvp
# CUDAnative.CUDAdrv.@profile
# check either on CPU with EnsembleCPUArray() or on GPU with EnsembleGPUArray()
@test_nowarn sol = @time solve(ensembleprob,DRI1(),EnsembleCPUArray(),trajectories=numtraj)
#sol = @time solve(ensembleprob,DRI1(),EnsembleGPUArray(),trajectories=numtraj)

# using Plots; plotly()  # or gr()
# using Plots; plot(sol,linealpha=0.6,color=:blue,vars=(0,1),title="Phase Space Plot")
# plot!(sol,linealpha=0.6,color=:red,vars=(0,2),title="Phase Space Plot")
# plot(sol,linealpha=0.6,color=:grey, vars=(1,2),title="Phase Space Plot")
#
# scatter(sol,linealpha=0.6,color=:blue,vars=(0,1),title="Phase Space Plot")
#
# summ = EnsembleSummary(sol,0.0f0:0.5f0:100f0)
# plot(summ,fillalpha=0.5)
#
#
# using DifferentialEquations.EnsembleAnalysis
# meansol = timeseries_steps_mean(sol)
# x1 = [x[1] for x in meansol.u]
# x2 = [x[2] for x in meansol.u]
# dts = []
# tmp1 = tspan[1]
# for tmp2 in meansol.t
#   global tmp1
#   push!(dts,tmp2-tmp1)
#   tmp1 = tmp2
# end
# #
# plot(x1,x2)
# plt = plot(dts)
# savefig(plt, "chosen_timesteps.png")

###
# oop not yet supported by DiffEqGPU
###
#
# function brusselator_f(u,p,t)
#  @inbounds begin
#    dx = (p[1]-1)*u[1]+p[1]*u[1]^2+(u[1]+1)^2*u[2]
#    dy = -p[1]*u[1]-p[1]*u[1]^2-(u[1]+1)^2*u[2]
#  end
#  return [dx,dy]
# end
#
# function scalar_noise(u,p,t)
#  @inbounds begin
#      dx = p[2]*u[1]*(1+u[1])
#      dy = -p[2]*u[1]*(1+u[1])
#  end
#  return [dx,dy]
# end
#
# function prob_func(prob, i, repeat)
#     Random.seed!(seeds[i])
#     W = WienerProcess(0.0,0.0,0.0)
#     remake(prob,noise=W)
# end
#
#
# # fix seeds
# seed = 100
# Random.seed!(seed)
# numtraj= 10
# seeds = rand(UInt, numtraj)
# W = WienerProcess(0.0,0.0,0.0)
#
# #CUDA.allowscalar(false)
# u0 = [-0.1f0,0.0f0]
# tspan = (0.0f0,100.0f0)
# p = [1.9f0,0.1f0]
#
# prob = SDEProblem(brusselator_f,scalar_noise,u0,tspan,p, noise=W)
#
#
# ensembleprob = EnsembleProblem(prob, prob_func = prob_func)
# sol = solve(ensembleprob,DRI1(),EnsembleCPUArray(),trajectories=numtraj)
#


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

"""
Problem u(X) = u₀ in Γ(d)
"""
type Dirichlet <: BoundaryProblem
    formulation :: Symbol
    variational :: Bool
    dual_basis :: Bool
    order :: Int
end

function Dirichlet()
    Dirichlet(:incremental, false, false, 1)
end

function get_unknown_field_name(::Type{Dirichlet})
    return "reaction force"
end

function get_formulation_type(problem::Problem{Dirichlet})
    return problem.properties.formulation
end

function assemble!(problem::Problem{Dirichlet}, time::Float64=0.0;
                   auto_initialize=true)
    # FIXME: boilerplate
    if !isempty(problem.assembly)
        warn("Assemble problem $(problem.name): problem.assembly is not empty and assembling, are you sure you know what are you doing?")
    end
    if isempty(problem.elements)
        warn("Assemble problem $(problem.name): problem.elements is empty, no elements in problem?")
    else
        first_element = first(problem.elements)
        unknown_field_name = get_unknown_field_name(problem)
        if !haskey(first_element, unknown_field_name)
            warn("Assemble problem $(problem.name): seems that problem is uninitialized.")
            if auto_initialize
                info("Initializing problem $(problem.name) at time $time automatically.")
                initialize!(problem, time)
            end
        end
    end

    if method_exists(assemble_prehook!, Tuple{typeof(problem), Float64})
        assemble_prehook!(problem, time)
    end

    if problem.properties.variational
        for element in get_elements(problem)
            assemble!(problem.assembly, problem, element, time)
        end
    else # nodal collocation
        field_vals = Dict{Int64, Float64}()
        field_name = get_parent_field_name(problem)
        field_dim = get_unknown_field_dimension(problem)
        for element in get_elements(problem)
            gdofs = get_gdofs(problem, element)
            for i=1:field_dim
                haskey(element, field_name*" $i") || continue
                ldofs = gdofs[i:field_dim:end]
                xis = get_reference_coordinates(typeof(element.properties))
                vals = Float64[]
                for xi in xis
                    g = element(field_name*" $i", xi, time)
                    # u = u_prev + Δu ⇒ Δu = u - u_prev
                    if haskey(element, field_name)
                        g_prev = element(field_name, xi, time)
                        g -= g_prev[i]
                    end
                    push!(vals, g)
                end
                for (dof, g) in zip(ldofs, vals)
                    field_vals[dof] = g
                end
            end
        end
        for (k, v) in field_vals
            push!(problem.assembly.C1, k, k, 1.0)
            push!(problem.assembly.C2, k, k, 1.0)
            push!(problem.assembly.g, k, 1, v)
        end
    end

    if method_exists(assemble_posthook!, Tuple{typeof(problem), Float64})
        assemble_posthook!(problem, time)
    end
end

function assemble!(assembly::Assembly, problem::Problem{Dirichlet},
                   element::Element, time::Float64)

    # get dimension and name of PARENT field
    nnodes = length(element)
    field_dim = get_unknown_field_dimension(problem)
    field_name = get_parent_field_name(problem)
    gdofs = get_gdofs(element, field_dim)
    props = problem.properties

    if problem.properties.dual_basis
        De, Me, Ae = get_dualbasis(element, time)
    else
        Ae = eye(nnodes)
        De = zeros(nnodes, nnodes)
        for ip in get_integration_points(element, props.order)
            N = element(ip, time)
            detJ = element(ip, time, Val{:detJ})
            De += ip.weight*N'*N*detJ
        end
    end

    # left hand side
    for i=1:field_dim
        ldofs = gdofs[i:field_dim:end]
        if haskey(element, field_name*" $i")
            add!(assembly.C1, ldofs, ldofs, De)
            add!(assembly.C2, ldofs, ldofs, De)
        end
    end

    # right hand side
    for ip in get_integration_points(element, props.order)
        detJ = element(ip, time, Val{:detJ})
        w = ip.weight*detJ
        N = element(ip, time)

        for i=1:field_dim
            ldofs = gdofs[i:field_dim:end]
            if haskey(element, field_name*" $i")
                g = element(field_name*" $i", ip, time)
                # u = u_prev + Δu ⇒ Δu = u - u_prev
                if haskey(element, field_name)
                    g_prev = element(field_name, ip, time)
                    g -= g_prev[i]
                end
                add!(assembly.g, ldofs, w*g*Ae*N')
            end
        end

    end

end

#=

function assemble!(assembly::Assembly, problem::Problem{DirichletProblem}, element::Element, time::Real)

    # get dimension and name of PARENT field
    field_dim = problem.parent_field_dim
    field_name = problem.parent_field_name

    gdofs = get_gdofs(element, field_dim)
    for ip in get_integration_points(element, Val{2})
        w = ip.weight
        J = get_jacobian(element, ip, time)
        JT = transpose(J)
        if size(JT, 2) == 1  # plane problem
            w *= norm(JT)
        else
            w *= norm(cross(JT[:,1], JT[:,2]))
        end
        N = element(ip, time)
        A = w*N'*N

        if haskey(element, field_name)
            # add all dimensions at once if defined
            # element["blaa"] = 0.0
            # or
            # element["blaa"] = Vector{Float64}[[0.1, 0.2], [0.3, 0.4]]
            g = element(field_name, ip, time)
            if length(g) != length(N)
                g = g*ones(length(N))
            end
            for i=1:field_dim
                ldofs = gdofs[i:field_dim:end]
                add!(assembly.C1, ldofs, ldofs, A)
                add!(assembly.C2, ldofs, ldofs, A)
            end
            add!(assembly.g, gdofs, w*g*N)
        end

        for i=1:field_dim
            if haskey(element, field_name*" $i")
                g = element(field_name*" $i", ip, time)
                ldofs = gdofs[i:field_dim:end]
                add!(assembly.C1, ldofs, ldofs, A)
                add!(assembly.C2, ldofs, ldofs, A)
                add!(assembly.g, ldofs, w*g*N)
            end
        end
    end
end
=#


text:
function divide_at_mean_seq(xs)
    m = mean(xs)
    larger = Vector{Int}(undef, length(xs))
    smaller = Vector{Int}(undef, length(xs))
    n_larger = 0
    n_smaller = 0
    for (i, x) in enumerate(xs)
        is_larger = x > m
        n_larger += is_larger
        n_smaller += !is_larger
        dest = ifelse(is_larger, larger, smaller)
        n = ifelse(is_larger, n_larger, n_smaller)
        @inbounds dest[n] = i
    end
    return resize!(smaller, n_smaller), resize!(larger, n_larger)
end

function divide_at_mean_threads(xs)
    m = ThreadsFolds.mean(xs)
    larger = Vector{Int}(undef, length(xs))
    smaller = Vector{Int}(undef, length(xs))
    n_larger = 0
    n_smaller = 0
    for (i, x) in enumerate(xs)
        is_larger = x > m
        n_larger += is_larger
        n_smaller += !is_larger
        dest = ifelse(is_larger, larger, smaller)
        n = ifelse(is_larger, n_larger, n_smaller)
        @inbounds dest[n] = i
    end
    return resize!(smaller, n_smaller), resize!(larger, n_larger)
end

function divide_at_mean_tapir(xs)
    m = TapirFolds.mean(xs)
    larger = Vector{Int}(undef, length(xs))
    smaller = Vector{Int}(undef, length(xs))
    n_larger = 0
    n_smaller = 0
    for (i, x) in enumerate(xs)
        is_larger = x > m
        n_larger += is_larger
        n_smaller += !is_larger
        dest = ifelse(is_larger, larger, smaller)
        n = ifelse(is_larger, n_larger, n_smaller)
        @inbounds dest[n] = i
    end
    return resize!(smaller, n_smaller), resize!(larger, n_larger)
end


text:
using Cuba

denom(Σinv::Matrix{Float64}) :: BigFloat  =  sqrt( (2*π)^3 / det(Σinv) )
function gaussint(Σinv::Matrix{Float64}, x::Vector{Float64}) :: BigFloat
    t = - BigFloat( (x'*Σinv*x)[1] )/2
    return exp(t)
    ;
end

function integrate(Σinv::Matrix{Float64}, a::Vector{Float64}, b::Vector{Float64}, alg, n::Int=100) :: Float64
    local d::Vector{Float64}
    function fun(x::Vector{Float64}, f::Vector{Float64})
        tmp = BigFloat(0)
        y = a + x.*d
        for di in -1:+1
            for dj in -1:+1
                for dk in -1:+1
                    v = gaussint(Σinv, y+[di/1000000,dj/1000000,dk/1000000].*d)
                    di!=0 || dj!=0 || dk!=0 || ( v*=6 )
                    tmp += v
                end
            end
        end
        f[1] = Float64( tmp )
    end

    divi = denom(Σinv)
    d = b-a

    if alg==:vegas
        integral,error,probability,neval,fail,nregions = vegas(fun,3,1)
        integral[1] *= d[1]*d[2]*d[3]/divi/32
        error[1]    *= d[1]*d[2]*d[3]/divi/32
    elseif alg==:suave
        integral,error,probability,neval,fail,nregions = suave(fun,3,1)
        integral[1] *= d[1]*d[2]*d[3]/divi/32
        error[1]    *= d[1]*d[2]*d[3]/divi/32
    elseif alg==:divonne
        integral,error,probability,neval,fail,nregions = divonne(fun,3,1)
        integral[1] *= d[1]*d[2]*d[3]/divi/32
        error[1]    *= d[1]*d[2]*d[3]/divi/32
    elseif alg==:cuhre
        integral,error,probability,neval,fail,nregions = cuhre(fun,3,1)
        integral[1] *= d[1]*d[2]*d[3]/divi/32
        error[1]    *= d[1]*d[2]*d[3]/divi/32
    elseif alg==:grid
        x = Vector{Float64}([0.,0.,0.])
        avg = BigFloat(0)
        for i = 0:(n-1)
            x[1] = ( i*b[1] + (n-i)*a[1] )/n
            for j = 0:(n-1)
                x[2] = ( j*b[2] + (n-j)*a[2] )/n
                for k = 0:(n-1)
                    x[3] = ( k*b[3] + (n-k)*a[3] )/n
                    p = gaussint(Σinv, x) / d
                    avg +=  p / n^3
                end
            end
        end
        integral,error,probability,neval,fail,nregions = avg,0,0,n,0,0
    else
        @assert false "Wrong Algorithm symbol"
    end
    @show integral error probability neval fail nregions
    return integral[1]
end

function make_apxgaussian_pdf(Σ::Matrix{Float64}, blocklength::Float64, zero::Float64=1.e-50, alg=:divonne, n::Int = 100) :: Dict{Tuple{Float64,Float64,Float64},Float64}
    @assert n>=4 "make_apxgaussian_pdf(): n>=4 needed"
    Σinv = inv(Σ)

    pdf = Dict{Tuple{Float64,Float64,Float64},Float64}()

    prc = precision(BigFloat)
    setprecision(BigFloat,512)

    boundary = Array{Tuple{Int,Int,Int}}(1)
    empty!(boundary)
    # Do central box, initialize boundary.
    i_,j_,k_ = (0,0,0)
    a = [ i_   *blocklength,   j_   *blocklength,   k_   *blocklength]
    b = [(i_+1)*blocklength,  (j_+1)*blocklength,  (k_+1)*blocklength]
    c = (a.+b)./2
    p = Float64( integrate(Σinv,a,b,alg,n) )
    if p > zero
        pdf[ (c[1],c[2],c[3]) ] = p
        push!(boundary, (i_,j_,k_) )
        # @show (i_,j_,k_) p
    end

    while !isempty(boundary)
        i,j,k = pop!(boundary)
        # @show (i,j,k) length(boundary)
        for di in -1:+1
            for dj in -1:+1
                for dk in -1:+1
                    i_,j_,k_ = (i+di,j+dj,k+dk)
                    a = [ i_   *blocklength,   j_   *blocklength,   k_   *blocklength]
                    b = [(i_+1)*blocklength,  (j_+1)*blocklength,  (k_+1)*blocklength]
                    c = (a.+b)./2
                    if (c[1],c[2],c[3]) ∉ keys(pdf)
                        p = Float64( integrate(Σinv,a,b,alg,n) )
                        # @show c p
                        if p > zero
                            pdf[ (c[1],c[2],c[3]) ] = p
                            push!(boundary, (i_,j_,k_) )
                            # @show (i_,j_,k_) p
                        end
                    end
                end
            end
        end
    end

    setprecision(BigFloat,prc)
    return pdf
    ;
end

function random_Σ() :: Matrix{Float64}
    A=ones(3,3)-2*rand(3,3);
    for i = 1:3
        A[:,i] ./= norm(A[:,i],2)
    end
    return A'*A
    ;
end

max_stars_repo_name: dot-at/newpid

max_issues_repo_name: dot-at/newpid

max_forks_repo_name: dot-at/newpid