text:
# Problem evaluation #
# ========================== #
function eval_first_stage(stochasticprogram::StochasticProgram, x::AbstractVector)
    first_stage = get_stage_one(stochasticprogram)
    return evaluate_objective(objective_function(first_stage), x)
end
# Evaluation API #
# ========================== #
"""
    evaluate_decision(stochasticprogram::TwoStageStochasticProgram, decision::AbstractVector)

Evaluate the first-stage `decision` in `stochasticprogram`.

In other words, evaluate the first-stage objective at `decision` and solve outcome models of `decision` for every available scenario. The supplied `decision` must match the defined decision variables in `stochasticprogram`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function evaluate_decision(stochasticprogram::TwoStageStochasticProgram, decision::AbstractVector)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Ensure stochastic program has been generated at this point
    if deferred(stochasticprogram)
        generate!(stochasticprogram)
    end
    # Sanity checks on given decision vector
    all(.!(isnan.(decision))) || error("Given decision vector has NaN elements")
    # Restore structure (if necessary)
    restore_structure!(stochasticprogram.optimizer)
    # Dispatch evaluation on stochastic structure
    return evaluate_decision(structure(stochasticprogram), decision)
end
"""
    statistically_evaluate_decision(stochasticprogram::TwoStageStochasticProgram, decision::AbstractVector)

Statistically evaluate the first-stage `decision` in `stochasticprogram`, returning the evaluated value and the spread over the scenarios.

The supplied `decision` must match the defined decision variables in `stochasticprogram`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function statistically_evaluate_decision(stochasticprogram::TwoStageStochasticProgram, decision::AbstractVector)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Ensure stochastic program has been generated at this point
    if deferred(stochasticprogram)
        generate!(stochasticprogram)
    end
    # Sanity checks on given decision vector
    all(.!(isnan.(decision))) || error("Given decision vector has NaN elements")
    # Dispatch evaluation on stochastic structure
    return statistically_evaluate_decision(structure(stochasticprogram), decision)
end
"""
    evaluate_decision(stochasticprogram::TwoStageStochasticProgram, decision::AbstractVector, scenario::AbstractScenario)

Evaluate the result of taking the first-stage `decision` if `scenario` is the actual outcome in `stochasticprogram`. The supplied `decision` must match the defined decision variables in `stochasticprogram`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function evaluate_decision(stochasticprogram::TwoStageStochasticProgram,
                           decision::AbstractVector,
                           scenario::AbstractScenario)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticprogram.optimizer)
    # Sanity checks on given decision vector
    length(decision) == num_decisions(stochasticprogram) || error("Incorrect length of given decision vector, has ", length(decision), " should be ", num_decisions(stochasticprogram))
    all(.!(isnan.(decision))) || error("Given decision vector has NaN elements")
    # Generate and solve outcome model
    outcome = outcome_model(stochasticprogram, decision, scenario; optimizer = subproblem_optimizer(stochasticprogram))
    optimize!(outcome)
    status = termination_status(outcome)
    if status in AcceptableTermination
        return objective_value(outcome)
    else
        if status == MOI.INFEASIBLE
            return objective_sense(outcome) == MOI.MAX_SENSE ? -Inf : Inf
        elseif status == MOI.DUAL_INFEASIBLE
            return objective_sense(outcome) == MOI.MAX_SENSE ? Inf : -Inf
        else
            error("Outcome model could not be solved, returned status: $status")
        end
    end
end
"""
    evaluate_decision(stochasticmodel::StochasticModel{2}, decision::AbstractVector, sampler::AbstractSampler; kw...)

Return a statistical estimate of the objective of the two-stage `stochasticmodel` at `decision` in the form of a confidence interval at the current confidence level, over the scenario distribution induced by `sampler`.

In other words, evaluate `decision` on a sampled model and generate an confidence interval using the sample variance of the evaluation. The confidence level can be set through the [`Confidence`](@ref) attribute and the sample size can be set through the [`NumEvalSamples`](@ref) attribute.

The supplied `decision` must match the defined decision variables in `stochasticmodel`. If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.

See also: [`confidence_interval`](@ref)
"""
function evaluate_decision(stochasticmodel::StochasticModel{2}, decision::AbstractVector, sampler::AbstractSampler; kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    N = MOI.get(stochasticmodel, NumEvalSamples())
    # Calculate confidence interval using provided optimizer
    CI = let eval_model = instantiate(stochasticmodel, sampler, N; optimizer = optimizer, kw...)
        # Silence output
        MOI.set(eval_model, MOI.Silent(), true)
        # Sanity checks on given decision vector
        length(decision) == num_decisions(eval_model) || error("Incorrect length of given decision vector, has ", length(decision), " should be ", num_decisions(eval_model))
        all(.!(isnan.(decision))) || error("Given decision vector has NaN elements")
        # Confidence level
        α = 1 - confidence
        𝔼Q, σ = statistically_evaluate_decision(eval_model, decision)
        z = quantile(Normal(0,1), 1 - α)
        L = 𝔼Q - z * σ / sqrt(N)
        U = 𝔼Q + z * σ / sqrt(N)
        # Clear memory from temporary model
        clear!(eval_model)
        return ConfidenceInterval(L, U, confidence)
    end
    return CI
end
"""
    lower_bound(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; kw...)

Generate a confidence interval around a lower bound on the true optimum of the two-stage `stochasticmodel` at the current confidence level, over the scenario distribution induced by `sampler`.

The attribute [`NumSamples`](@ref) is the size of the sampled models used to generate the interval and generally governs how tight it is. The attribute [`NumLowerTrials`](@ref) is the number of sampled models. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function lower_bound(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get the instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    α = 1 - confidence
    N = MOI.get(stochasticmodel, NumSamples())
    M = MOI.get(stochasticmodel, NumLowerTrials())
    log = MOI.get(stochasticmodel, MOI.RawParameter("log"))
    keep = MOI.get(stochasticmodel, MOI.RawParameter("keep"))
    offset = MOI.get(stochasticmodel, MOI.RawParameter("offset"))
    indent = MOI.get(stochasticmodel, MOI.RawParameter("indent"))
    # Lower bound
    Qs = Vector{Float64}(undef, M)
    progress = Progress(M, 0.0, "$(repeat(" ", indent))Lower CI    ")
    log && sleep(0.1)
    log && ProgressMeter.update!(progress, 0, keep = false, offset = offset)
    for i = 1:M
        let sampled_model = instantiate(stochasticmodel, sampler, N; optimizer = optimizer, kw...)
            Qs[i] = VRP(sampled_model)
            # Clear memory from temporary model
            clear!(sampled_model)
        end
        log && ProgressMeter.update!(progress, i, keep = keep, offset = offset)
    end
    Q̂ = mean(Qs)
    σ = std(Qs)
    t = quantile(TDist(M-1), 1-α)
    L = Q̂ - t * σ / sqrt(M)
    U = Q̂ + t * σ / sqrt(M)
    return ConfidenceInterval(L, U, 1-α)
end
"""
    upper_bound(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; kw...)

Generate a confidence interval around an upper bound of the true optimum of the two-stage `stochasticmodel` at the current confidence level, over the scenario distribution induced by `sampler`, by generating and evaluating a candidate decision.

The attribute [`NumSamples`](@ref) is the size of the sampled model used to generate a candidate decision. The attribute [`NumUpperTrials`](@ref) is the number of sampled models and the attribute [`NumEvalSamples`](@ref) is the size of the evaluation models. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function upper_bound(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get the instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    α = 1 - confidence
    num_samples = MOI.get(stochasticmodel, NumSamples())
    # decision generation
    sampled_model = instantiate(stochasticmodel,
                                sampler,
                                num_samples;
                                optimizer = optimizer,
                                kw...)
    # Optimize
    optimize!(sampled_model)
    x̂ = optimal_decision(sampled_model)
    return upper_bound(stochasticmodel, x̂, sampler; kw...)
end
"""
    upper_bound(stochasticmodel::StochasticModel{2}, decision::AbstractVector, sampler::AbstractSampler; kw...)

Generate a confidence interval around an upper bound of the expected value of `decision` in the two-stage `stochasticmodel` at the current confidence level, over the scenario distribution induced by `sampler`.

The attribute [`NumUpperTrials`](@ref) is the number of sampled models and the attribute [`NumEvalSamples`](@ref) is the size of the evaluation models. The confidence level can be set through the [`Confidence`](@ref) attribute.

The supplied `decision` must match the defined decision variables in `stochasticmodel`. If an optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function upper_bound(stochasticmodel::StochasticModel{2}, decision::AbstractVector, sampler::AbstractSampler, kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Get the instance optimizer
    optimizer = MOI.get(stochasticmodel, InstanceOptimizer())
    # Get parameters
    confidence = MOI.get(stochasticmodel, Confidence())
    α = 1 - confidence
    N = MOI.get(stochasticmodel, NumEvalSamples())
    N = max(N, MOI.get(stochasticmodel, NumSamples()))
    T = MOI.get(stochasticmodel, NumUpperTrials())
    log = MOI.get(stochasticmodel, MOI.RawParameter("log"))
    keep = MOI.get(stochasticmodel, MOI.RawParameter("keep"))
    offset = MOI.get(stochasticmodel, MOI.RawParameter("offset"))
    indent = MOI.get(stochasticmodel, MOI.RawParameter("indent"))
    # Generate upper bound
    Q = Vector{Float64}(undef, T)
    progress = Progress(T, 0.0, "$(repeat(" ", indent))Upper CI    ")
    log && sleep(0.1)
    log && ProgressMeter.update!(progress, 0, keep = false, offset = offset - 1)
    for i = 1:T
        let eval_model = instantiate(stochasticmodel, sampler, N; optimizer = optimizer, kw...)
            # Silence output
            MOI.set(eval_model, MOI.Silent(), true)
            # Sanity checks on given decision vector
            length(decision) == num_decisions(eval_model) || error("Incorrect length of given decision vector, has ", length(decision), " should be ", num_decisions(eval_model))
            all(.!(isnan.(decision))) || error("Given decision vector has NaN elements")
            # Evaluate on sampled model
            Q[i] = evaluate_decision(eval_model, decision)
            # Clear memory from temporary model
            clear!(eval_model)
        end
        log && ProgressMeter.update!(progress, i, keep = keep, offset = offset - 1)
    end
    Q̂ = mean(Q)
    σ = std(Q)
    t = quantile(TDist(T - 1), 1 - α)
    L = Q̂ - t * σ / sqrt(T)
    U = Q̂ + t * σ / sqrt(T)
    return ConfidenceInterval(L, U, 1 - α)
end
"""
    confidence_interval(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler)

Generate a confidence interval around the true optimum of the two-stage `stochasticmodel` at level `confidence` using SAA, over the scenario distribution induced by `sampler`.

The attribute [`NumSamples`](@ref) is the size of the sampled models used to generate the interval and generally governs how tight it is. The attribute [`NumLowerTrials`](@ref) is the number of sampled models used in the lower bound calculation and the attribute [`NumUpperTrials`](@ref) is the number of sampled models used in the upper bound calculation. The attribute [`NumEvalSamples`](@ref) is the size of the sampled models used in the upper bound calculation. The confidence level can be set through the [`Confidence`](@ref) attribute.

If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function confidence_interval(stochasticmodel::StochasticModel{2}, sampler::AbstractSampler; kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Confidence level
    confidence = MOI.get(stochasticmodel, Confidence())
    # Modify confidence for two-sided interval
    α = (1 - confidence)/2
    MOI.set(stochasticmodel, Confidence(), 1 - α)
    # Lower bound
    lower_CI = lower_bound(stochasticmodel, sampler; kw...)
    L = lower(lower_CI)
    # Upper bound
    upper_CI = upper_bound(stochasticmodel, sampler; kw...)
    U = upper(upper_CI)
    # Restore confidence level
    MOI.set(stochasticmodel, Confidence(), confidence)
    # Check if confidence interval is valid
    if U <= L
        @warn "Could not calculate confidence interval at current level of confidence and sample size"
        return ConfidenceInterval(-Inf, Inf, confidence)
    end
    # Return confidence interval
    return ConfidenceInterval(L, U, confidence)
end
"""
    gap(stochasticmodel::StochasticModel{2}, decision::UnionAbstractVector, sampler::AbstractSampler)

Generate a confidence interval around the gap between the result of using `decision` and the true optimum of the two-stage `stochasticmodel` at the current confidence level, over the scenario distribution induced by `sampler`.

The attribute [`NumSamples`](@ref) is the size of the sampled models used to generate the interval and generally governs how tight it is. The attribute [`NumLowerTrials`](@ref) is the number of sampled models used in the lower bound calculation and the attribute [`NumUpperTrials`](@ref) is the number of sampled models used in the upper bound calculation. The attribute [`NumEvalSamples`](@ref) is the size of the sampled models used in the upper bound calculation. The confidence level can be set through the [`Confidence`](@ref) attribute.

The supplied `decision` must match the defined decision variables in `stochasticmodel`. If a sample-based optimizer has not been set yet (see [`set_optimizer`](@ref)), a `NoOptimizer` error is thrown.
"""
function gap(stochasticmodel::StochasticModel{2}, decision::AbstractVector, sampler::AbstractSampler; kw...)
    # Throw NoOptimizer error if no recognized optimizer has been provided
    check_provided_optimizer(stochasticmodel.optimizer)
    # Confidence level
    confidence = MOI.get(stochasticmodel, Confidence())
    # Modify confidence for two-sided interval
    α = (1-confidence)/2
    MOI.set(stochasticmodel, Confidence(), 1 - α)
    # Lower bound
    lower_CI = lower_bound(stochasticmodel, sampler; kw...)
    L = lower(lower_CI)
    # Upper bound
    upper_CI = upper_bound(stochasticmodel, decision, sampler; kw...)
    U = upper(upper_CI)
    # Restore confidence level
    MOI.set(stochasticmodel, Confidence(), confidence)
    # Return confidence interval
    return ConfidenceInterval(0., U - L, confidence)
end


text:
"""
    classify!(array, weights[, classifyMask])

Classify an array in-place into proportions based upon a list of class weights.
"""
function classify!(array, weights, mask = nothing)
    quantiles = zeros(length(weights))
    cumsum!(quantiles, weights)
    quantiles ./= quantiles[end]
    boundaryvalues = isnothing(mask) ?
    quantile(filter(isfinite, array), quantiles) :
    quantile(array[mask .& isfinite.(array)], quantiles)
    for i in eachindex(array)
        array[i] = isnan(array[i]) ? NaN : searchsortedfirst(boundaryvalues, array[i])
    end
    array
end
classify!(array, weights::Real, mask = nothing) = classify!(array, ones(weights), mask)

"""
    classify(array, weights[, classifyMask])

Classify an array into proportions based upon a list of class weights.
"""
classify(array, weights, mask = nothing) = classify!(copy(array), weights, mask)
classify(array, weights::Real, mask = nothing) = classify(array, ones(weights), mask)

function _clusterMean(clusterArray, array)
    clusters = Dict{Float64, Float64}()
    clustersum = Dict{Float64, Float64}()
    labels, nlabels = label(clusterArray)
    for ind in eachindex(labels, array)
        temp = labels[ind]
        if !haskey(clusters, temp)
            clusters[temp] = clustersum[temp] = 0.0
        end
        clusters[temp] += 1.0
        clustersum[temp] += array[ind]
    end
    for cl in keys(clusters)
        clustersum[cl] /= clusters[cl]
    end
    clustersum[NaN] = NaN
    _rescale!(get.(Ref(clustersum), labels, NaN))
end

"""
    blend(arrays[, scaling])

Blend arrays weighted by scaling factors.
"""
function blend(arrays, scaling::AbstractVector{<:Number} = ones(length(arrays)))
    if length(scaling) != length(arrays)
        throw(DimensionMismatch("The array of landscapes (n = $(length(arrays))) and scaling (n = $(length(scaling))) must have the same length"))
    end
    ret = sum(arrays .* scaling)
    _rescale!(ret)
end

"""
    blend(clusterarray, arrays[, scaling])

Blend a primary cluster NLM with other arrays in which the mean value per
cluster is weighted by scaling factors.
"""
function blend(clusterarray, arrays::AbstractVector, scaling::AbstractVector{<:Number} = ones(length(arrays)))
    ret = sum(_clusterMean.(Ref(clusterarray), arrays) .* scaling)
    _rescale!(clusterarray + ret)
end
blend(clusterarray, array, scaling = 1) = blend(clusterarray, [array], [scaling])

const _neighborhoods = Dict(
    :rook     => ((1, 0), (0, 1)),
    :diagonal => ((1, 0), (0, 1), (1, 1)),
    :queen    => ((1, 0), (0, 1), (1, 1), (1, -1)),
)

"""
    label(mat[, neighborhood])

Assign an arbitrary label to all clusters of contiguous matrix elements with the same value.
Returns a matrix of values and the total number of final clusters.
The `neighborhood` structure can be
`:rook`     `:queen`    `:diagonal`
 0 1 0        1 1 1        0 1 1
 1 x 1        1 x 1        1 x 1
 0 1 0        1 1 1        1 1 0
`:rook` is the default
"""
function label(mat, neighborhood = :rook)
    neighbors = _neighborhoods[neighborhood]
    m, n = size(mat)
    (m >= 3 && n >= 3) || error("The label algorithm requires the landscape to be at least 3 cells in each direction")

    # initialize objects
    ret = similar(mat)
    clusters = IntDisjointSets(0)

    # run through the matrix and make clusters
    for j in axes(mat, 2), i in axes(mat, 1)
        if isfinite(mat[i, j])
            same = [i - n[1] > 0 && j - n[2] > 0 && mat[i - n[1], j - n[2]] == mat[i, j] for n in neighbors]
            if count(same) == 0
                push!(clusters)
                ret[i, j] = length(clusters)
            elseif count(same) == 1
                n1, n2 = only(neighbors[same])
                ret[i, j] = ret[i - n1, j - n2]
            else
                vals = unique([ret[i - n[1], j - n[2]] for n in neighbors[same]])
                if length(vals) == 1
                    ret[i, j] = only(vals)
                else
                    for v in vals[2:end]
                        ret[i, j] = union!(clusters, Int(vals[1]), Int(v))
                    end
                end
            end
        else
            ret[i, j] = NaN
        end
    end

    # merge adjacent clusters with same value
    finalclusters = Set{eltype(ret)}()
    for i in eachindex(ret)
        ret[i] = isnan(ret[i]) ? NaN : find_root(clusters, Int(ret[i]))
        push!(finalclusters, ret[i])
    end

    # assign each cluster a random number in steps of 1 (good for plotting)
    randomcode = Dict(i => j for (i, j) in zip(finalclusters, 1.0:length(finalclusters)))
    randomcode[NaN] = NaN
    for i in eachindex(ret)
        ret[i] = randomcode[ret[i]]
    end

    ret, length(finalclusters)
end

ext: jl


text:
using .lib

module analyze
	export argAnalyze
	import ..Data,..CMHelp,..CMVersion,..MMNone,..MMSerial,..MMSpawn,..MMThread,..eq,..exitWithError,..clean

	@enum Key AKNull AKStdout AKStderr AKResult AKMultiple
	function argAnalyze(d::Data)
		l=ARGS

		if length(l)==0 exitWithError("引数が不足しています") end
		if eq(l[1],"-h","help","-help","--help") d.mode=CMHelp end
		if eq(l[1],"-v","version","-version","--version") d.mode=CMVersion end

		key::Key=AKNull
		n::Int=0

		for a in l
			n+=1
			if a=="" continue end

			proceed::Bool=true
			if eq(a,"-m","-multiple")
				d.multiple=MMSerial
				key=AKMultiple
			elseif eq(a,"-o","-out","-stdout") key=AKStdout
			elseif eq(a,"-e","-err","-stderr") key=AKStderr
			elseif eq(a,"-r","-result") key=AKResult
			else proceed=false end
			if proceed continue end

			if startswith(a,"-")
				exitWithError("不正なオプションが指定されています")
			elseif key!=AKNull
				proceed=true
				if key==AKStdout d.out=a end
				if key==AKStderr d.err=a end
				if key==AKResult d.result=a end
				if key==AKMultiple
					if eq(a,"none") d.multiple=MMNone
					elseif eq(a,"serial") d.multiple=MMSerial
					elseif eq(a,"spawn","parallel") d.multiple=MMSpawn
					elseif eq(a,"thread") d.multiple=MMThread
					else proceed=false end
				end
				key=AKNull
			end
			if proceed continue end

			d.command=view(l,n:length(l))
			break
		end

		if length(d.command)==0 exitWithError("実行する内容が指定されていません") end
	end

end

max_stars_repo_path: Time Measurement/Julia/analyze.jl

max_issues_repo_path: Time Measurement/Julia/analyze.jl

max_forks_repo_path: Time Measurement/Julia/analyze.jl


text:
mutable struct Tree
    parent
    left
    right
end

function Tree(parent)
    Tree(parent, nothing, nothing)
end

function Base.show(io::IO, tree::Tree)
    tree isa Number && return print(tree)
    print("[")
    print(tree.left)
    print(",")
    print(tree.right)
    print("]")
end

function make_tree(list, parent=nothing)
    list isa Number && return list
    curr = Tree(parent)
    curr.left = make_tree(list[1], curr)
    curr.right = make_tree(list[2], curr)
    return curr
end

function explode(curr, left_leaf, right_leaf)
    if left_leaf != nothing
        left_leaf.right isa Number ? left_leaf.right += curr.left : left_leaf.left += curr.left
    end
    if right_leaf != nothing
        right_leaf.left isa Number ? right_leaf.left += curr.right : right_leaf.right += curr.right
    end
    if curr.parent.left === curr
        curr.parent.left = 0
    else
        curr.parent.right = 0
    end
    return true
end

function find_leaves(curr, depth=0)
    curr isa Number && return []
    mid = curr.left isa Number || curr.right isa Number ? [(depth, curr)] : []
    vcat(find_leaves(curr.left, depth+1), mid, find_leaves(curr.right, depth+1))
end

function find_next_explode(tree)
    leaves = vcat([(0, nothing)], find_leaves(tree), [(0, nothing)])
    for i in eachindex(leaves)
        if leaves[i][1] >= 4 && leaves[i][2].left isa Number && leaves[i][2].right isa Number
            return explode(leaves[i][2], leaves[i-1][2], leaves[i+1][2])
        end
    end
    return false
end

function find_next_split(tree)
    found_split = false
    should_split(num) = num isa Number && num > 9
    function split_descend(curr)
        split_up(num) = Tree(curr, num ÷ 2, num - num ÷ 2)
        curr isa Number && return false
        found_split && return true
        should_split(curr.left) && (curr.left = split_up(curr.left)) != 1 && return found_split = true
        split_descend(curr.left) && return true
        should_split(curr.right) && (curr.right = split_up(curr.right)) != 1 && return found_split = true
        split_descend(curr.right) && return true
        return false
    end
    return split_descend(tree)
end

function reduce_tree(tree)
    i = 0
    while true
        find_next_explode(tree) && continue
        find_next_split(tree) && continue
        return tree
    end
end

function merge_trees(left_tree, right_tree)
    new_root = Tree(nothing)
    left_tree.parent = new_root
    right_tree.parent = new_root
    new_root.left = left_tree
    new_root.right = right_tree
    return new_root
end

function add_trees(tree1, tree2)
    merge_trees(tree1, tree2) |> reduce_tree
end

function magnitude(tree)
    tree isa Number && return tree
    return 3magnitude(tree.left) + 2magnitude(tree.right)
end

function solve_part1(trees)
    curr = trees[1]
    for next in trees[2:end]
        merged = merge_trees(curr, next)
        curr = reduce_tree(merged)
    end
    println(magnitude(curr))
end

function solve_part2(trees)
    best = 0
    for t1 in trees, t2 in trees
        if !(t1 == t2)
            best = max(best, add_trees(deepcopy(t1), deepcopy(t2)) |> magnitude)
        end
    end
    println(best)
end

function main()
    trees = readlines() .|> Meta.parse .|> eval .|> make_tree
    solve_part1(deepcopy(trees))
    solve_part2(deepcopy(trees))
end

main()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
import Base: show, *, -, +, isless, ==, convert, conj, truncate, zero, one, promote_rule, A_mul_B!

export variables

immutable Symbols
    names :: Array{ASCIIString,1}
end

immutable Poly{T<:Number}
    n       :: Int             # number of symbols in monomials
    m       :: Int             # number of monomials
    syms    :: Symbols
    c       :: Array{T,1}      # coefficients in polynomial
    alpha   :: Array{Int, 2}   # powers of symbols for all monomials
    deg     :: Int

    function Poly(syms::Symbols, alpha::Array{Int,2}, c::Array{T,1})
        m, n = size(alpha)
        if length(m) == 0 error("Poly must include at least one monomial term") end
        if length(c) != m error("incompatible dimensions") end
        deg = size(alpha, 1) > 0 ? maximum(sum(alpha,2)) : 0
        new(n, m, syms, c, alpha, deg)
    end

    # we allow a constant polynomial without a prober symbol basis
    Poly(c::T) = new(0, 1, Symbols([""]), [c], Array(Int, 1, 0), 0)
end

Poly{T<:Number}(syms::Symbols, alpha::Array{Int,2}, c::Array{T,1}) = Poly{T}(syms, alpha, c)
Poly{T<:Number}(c::T) = Poly{T}(c)

convert{S,T}(::Type{Poly{S}}, p::Poly{T}) = Poly(p.syms, p.alpha, convert(Array{S}, p.c))

zero{T<:Number}(::Type{Poly{T}}) = Poly(zero(T))
zero{T<:Number}(::Poly{T}) = Poly(zero(T))
one{T<:Number}(::Type{Poly{T}}) = Poly(one(T))
one{T<:Number}(::Poly{T}) = Poly(one(T))

function promote_rule{S<:Number,T<:Number}(::Type{S}, ::Type{Poly{T}} )
    Poly{promote_type(S,T)}
end

function promote_rule{S<:Number,T<:Number}(::Type{Poly{S}}, ::Type{Poly{T}} )
    Poly{promote_type(S,T)}
end

promote_rule(::Type{Rational}, ::Type{Int}) = Rational

function show{T}(io::IO, p::Poly{T})
    if p.m == 0 print(io, zero(T)) end
    for i=1:p.m
        if i==1
            print(io, p.c[i]>=0 ? "" : "-")
        else
            print(io, p.c[i]>=0 ? "+" : "-")
        end

        ci = abs(p.c[i])
        if sum(p.alpha[i,:]) == 0
            print(io, ci)
        else
            if ci != one(ci) print(io, "$(ci)*") end
            for k=1:p.n
                if p.alpha[i,k] > 0
                    print(io, p.syms.names[k])
                    if p.alpha[i,k] > 1
                        print(io, "^$(p.alpha[i,k])")
                    end
                    if k < p.n && sum(p.alpha[i,k+1:end]) > 0
                        print(io, "*")
                    end
                end
            end
        end
    end
end

-(p::Poly) = Poly(p.syms, p.alpha, -p.c)
+(p::Poly, a::Number) = p + Poly(a)
+(a::Number, p::Poly) = Poly(a) + p
-(p::Poly, a::Number) = p - Poly(a)
-(a::Number, p::Poly) = Poly(a) - p
*(p::Poly, a::Number) = Poly(p.syms, p.alpha, p.c*a)
*(a::Number, p::Poly) = *(p::Poly, a::Number)

function +(p1::Poly, p2::Poly)
    p1, p2 = promote_poly(p1, p2)
    simplify(Poly(p1.syms, vcat(p1.alpha, p2.alpha), vcat(p1.c, p2.c)))
end

function -(p1::Poly, p2::Poly)
    p1, p2 = promote_poly(p1, p2)
    simplify(Poly(p1.syms, vcat(p1.alpha, p2.alpha), vcat(p1.c, -p2.c)))
end

function *{S<:Number,T<:Number}(p1::Poly{S}, p2::Poly{T})
    p1, p2 = promote_poly(p1, p2)
    if p2.m == 1
        Poly(p1.syms, p1.alpha .+ p2.alpha, p1.c*p2.c[1])
    else
        sum(Poly{promote_type(S, T)}[ p1*Poly(p2.syms, p2.alpha[k,:], [p2.c[k]]) for k=1:p2.m ])
    end
end

function ^{T}(p::Poly{T}, a::Int)
    if a<0 error("power must be nonnegative") end
    if a==0 return Poly(p.syms, zeros(Int, 1, p.n), [one(T)]) end
    if p.n == 1 return Poly(p.syms, p.alpha*a, p.c.^a) end
    prod([p for k=1:a])
end

typealias MatOrVec{T} Union(Array{T,1},Array{T,2})

*{T<:Number,S<:Number}(a::T, v::MatOrVec{Poly{S}}) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))
*{T<:Number,S<:Number}(v::MatOrVec{Poly{T}}, a::S) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))
*{T<:Number,S<:Number}(a::Poly{T}, v::MatOrVec{Poly{S}}) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))
*{T<:Number,S<:Number}(v::MatOrVec{Poly{S}}, a::Poly{T}) = reshape(Poly{promote_type(T,S)}[ a*vi for vi=v ], size(v))

==(p1::Poly, p2::Poly) = (p1-p2).m == 0

conj(p::Poly) = Poly(p.syms, p.alpha, conj(p.c))
convert{S<:Number,T<:Number}(::Type{Poly{S}}, a::T) = Poly{S}(convert(S,a))
isconst(p::Poly) = p.m == 1 && p.n == 0

function A_mul_B!{S<:Number,T<:Number,}(alpha::Poly{S}, A::SparseMatrixCSC{S,Int}, x::Array{Poly{T},1}, beta::Poly{S}, y::Array{Poly{S},1})
    for i=1:length(y)
        if beta == zero(y[i])
            y[i] = zero(y[i])
        else
            y[i] *= beta
        end
    end

    for j=1:size(A,2)
        for k=A.colptr[j]:A.colptr[j+1]-1
            y[A.rowval[k]] += alpha*A.nzval[k]*x[j]
        end
    end
    y
end

# promote polynomial to share the same Symbol basis
function promote_poly{S,T}(p1::Poly{S}, p2::Poly{T})
    if (p1.syms == p2.syms) || (isconst(p1) && isconst(p2)) return (p1, p2) end
    if isconst(p1) && ~isconst(p2)
        return (Poly(p2.syms, zeros(Int, 1, p2.n), p1.c), p2)
    elseif ~isconst(p1) && isconst(p2)
        return (p1, Poly(p1.syms, zeros(Int, 1, p1.n), p2.c))
    else
        error("incompatible monomial bases")
    end
end

function evalpoly{T<:Number,S<:Number}(p::Poly{T}, x::Vector{S})
    val = zero(S)
    for j=1:p.m
        val += p.c[j]*prod(vec(x) .^ vec(p.alpha[j,:]))
    end
    val
end

# combine identical terms and remove zero terms
function simplify{T<:Number}(p::Poly{T})
    lgt = 0
    labeled = falses(p.m)
    # labeled monomial terms are either zero, or have been merged with other terms.
    for k=1:p.m
        if ~labeled[k]
            for j=k+1:p.m
                if (~labeled[j] && p.alpha[k,:] == p.alpha[j,:])
                    p.c[k] += p.c[j]
                    labeled[j] = true
                end
            end
            if p.c[k] != zero(T)
                lgt += 1
            else
                labeled[k] = true
            end
        end
    end

    # the unlabeled terms form the reduced polynomial
    alpha = Array(Int, lgt, p.n)
    c = Array(T, lgt)
    lgt = 1
    for k=1:p.m
        if ~labeled[k]
            alpha[lgt,:] = p.alpha[k,:]
            c[lgt] = p.c[k]
            lgt += 1
        end
    end
    perm = sortperm([ vec(alpha[i,1:p.n]) for i=1:size(alpha,1)], lt=grilex_isless)
    Poly(p.syms, alpha[perm,:], c[perm])
end

function truncate{T<:Number}(p::Poly{T}, threshold=1e-10)
    lgt = 0
    labeled = falses(p.m)
    for k=1:p.m
        if abs(p.c[k]) < threshold
            labeled[k] = true
        else
            lgt += 1
        end
    end

    # the unlabeled terms form the reduced polynomial
    alpha = Array(Int, lgt, p.n)
    c = Array(T, lgt)
    lgt = 1
    for k=1:p.m
        if ~labeled[k]
            alpha[lgt,:] = p.alpha[k,:]
            c[lgt] = p.c[k]
            lgt += 1
        end
    end
    perm = sortperm([ vec(alpha[i,1:p.n]) for i=1:size(alpha,1)], lt=grilex_isless)
    Poly(p.syms, alpha[perm,:], c[perm])
end

# ordering for polynomials
function isless{T<:Number}(p1::Poly{T}, p2::Poly{T})
    # sort polynomials by degree
    if max(p1.m, p2.m) > 1 return p1.deg < p2.deg end

    # otherwise use graded inverse lex-order for monomials
    grilex_isless(vec(p1.alpha), vec(p2.alpha))
end

# graded inverse lexicographic order (lowest total order first, then reverse lex)
function grilex_isless(a::Vector{Int}, b::Vector{Int})
    i, j = sum(a), sum(b)
    i == j ? ~lexless(a, b) : i < j
end

variables(syms::Symbols) = [Poly{Int}(syms, [zeros(Int,1,k-1) 1 zeros(Int,1,length(syms.names)-k)], [1]) for k=1:length(syms.names)]
variables(syms::Vector{ASCIIString}) = variables(Symbols(syms))


text:
function sampleDiscrete(p)
  # Returns a sample from a discrete probability mass function indexed by p
  # (assumes that p is already normalized)
  r = rand()
  pcum = 0;
  for i = 1:length(p)
    pcum += p[i]
    if r <= pcum
      return i
    end
  end
  return length(p)
end


text:
using DCELs

@testset "Vertex" begin
    v = Vertex()
    @test v.id[1] == 'V'
    @test isnan(v.x)
    @test isnan(v.y)
    @test v.edge === nothing

    v = Vertex(1, -1)
    @test v.id[1] == 'V'
    @test v.x == 1
    @test v.y == -1
    @test v.edge === nothing
end

@testset "Edge" begin
    h = Edge()
    @test h.id[1] == 'E'
    @test h.next === nothing
    @test h.prev === nothing
    @test h.twin === nothing
    @test h.dest === nothing
    @test h.face === nothing
end

@testset "Face" begin
    f = Face()
    @test f.id[1] == 'F'
    @test f.edge === nothing
end

@testset "Tessellation" begin
    tess = Tessellation()
    @test tess.vertices isa AbstractSet{Vertex}
    @test tess.edges isa AbstractSet{Edge}
    @test tess.faces isa AbstractSet{Face}
    @test isempty(tess.vertices)
    @test isempty(tess.edges)
    @test isempty(tess.faces)
end

@testset "Arithmetic" begin
    u, v = Vertex(0, 0), Vertex(1, 1)
    a, b = u + v, u - v
    @test a.x == 1 && a.y == 1
    @test b.x == -1 && b.y == -1
end


text:
using BenchmarkTools, Plots

function plot_benchmark()
    res0 = BenchmarkTools.load(joinpath(@__DIR__, "OMEinsum_master.json"))[]
    res1 = BenchmarkTools.load(joinpath(@__DIR__, "OMEinsum_refactor.json"))[]
    t0 = Float64[]
    t1 = Float64[]
    xl = String[]
    for k in keys(res0)
        push!(t0, minimum(res0[k]["Float64"]["large"].times))
        push!(t1, minimum(res1[k]["Float64"]["large"].times))
        push!(xl, k)
    end
    xs = 1:length(res0)
    plot(xs, t0, marker = (:hexagon, 10, 0.6, :green, stroke(3, 0.2, :black, :dot)), yscale=:log10, lw=0, label="master")
    for i in xs
        annotate!(i+0.3, t1[i], text(xl[i], :black, :left, 12))
    end
    plot!(xs, t1, marker = (:circle, 10, 0.6, :red, stroke(3, 0.2, :black, :dot)), lw=0, label="refactor")
    #Plots.xticks!(xs, xl)
end

function plot_compile()
    pl = []
    for dynamic in [true, false]
        for uniqueindex in [true]
            times = readdlm(joinpath(@__DIR__, "compiletime_8_$(dynamic)_$(uniqueindex).dat"))
            push!(pl, plot(1:1000, times, label=(dynamic ? "dynamic" : "static") * " ($(round(sum(times); sigdigits=2)))", ylim=(0,2)))
        end
    end
    times = readdlm(joinpath(@__DIR__, "compiletime_8_false_true_master.dat"))
    push!(pl, plot(1:1000, times, label="master ($(round(sum(times); sigdigits=2)))", ylim=(0,2)))
    plt = plot(pl..., xlabel="iteration", ylabel="time", layout=(3,1))
end
#plot_benchmark()
plot_compile()


text:
export get_config

function get_base_config()
    (seed = 0,           # evolution seed
     d_fitness = 1,      # fitness dimension
     n_population = 2,   # population size
     n_elite = 1,        # elites to carry over each generation
     n_gen = 10,         # number of generations
     log_gen = 1,        # log every x generations
     save_gen = 1,       # save population every x generation
     id = Dates.format(Dates.now(), "yyyymmddHHMMSSsss") # run id
     )
end

"convert Dict to named tuple"
function get_config(cfg::Dict)
    (; (Symbol(k)=>v for (k, v) in cfg)...)
end

"combine YAML file and kwargs, make sure ID is specified"
function get_config(cfg_file::String; kwargs...)
    cfg = YAML.load_file(cfg_file)
    for (k, v) in kwargs
        cfg[String(k)] = v
    end
    # generate id, use date if no existing id
    if ~(:id in keys(cfg))
        cfg["id"] = Dates.format(Dates.now(), "yyyymmddHHMMSSsss")
    end
    get_config(cfg)
end


text:
# Additional tests following bug reports

@testset "i10" begin
    X, _ = MB.make_blobs(500, 3, rng=555)
    m = AgglomerativeClustering()
    f, = MB.fit(m, 1, X)
    @test f !== nothing
end


text:
using Toxcore, Reactive
using Base.Test
function load_client(path_to_client_file)
	if isfile(path_to_client_file)
		savefile = open(path_to_client_file, "r")
		savedata = readbytes(savefile)
		close(savefile)
		tox = tox_new(tox_options_default(), savedata)
		return tox
	end
	return tox_new()
end
immutable ToxcoreClient
	tox::Ptr{Tox}
	signals::Dict{Symbol, Signal}
	function ToxcoreClient(path_to_client_file=Pkg.dir("Toxcore", "client", "bot_savedata.binary"))
		tox = load_client(path_to_client_file)
		new(tox, Dict{Symbol, Signal}())
	end
end
client = ToxcoreClient()


text:
function _get_julia_url(config::Config)
    return _get_julia_url(config.julia_version)
end

function _get_julia_url(version::AbstractString)
    if version == "nightly"
        julia_url = "https://julialangnightlies-s3.julialang.org/bin/linux/x64/julia-latest-linux64.tar.gz"
        asc_url = "$(julia_url).asc"
        return julia_url, asc_url
    else
        return _get_julia_url(VersionNumber(version))
    end
end

function _get_julia_url(version::VersionNumber)
    major = version.major
    minor = version.minor
    patch = version.patch
    julia_url = "https://julialang-s3.julialang.org/bin/linux/x64/$(major).$(minor)/julia-$(major).$(minor).$(patch)-linux-x86_64.tar.gz"
    asc_url = "$(julia_url).asc"
    return julia_url, asc_url
end


text:
mutable struct Network
    A
    B
    N
    T
    window
    read_layer
    write_layer
    embed_layer
    qnetwork
    encoder
    decoder
    encoder_hidden
    decoder_hidden
end

function Network(A, B, N, T, encoder_dim, decoder_dim, noise_dim;
              embed=0, nclass=10, atype=_atype, read_attn=true, write_attn=true)
    imgsize = A*B

    window = nothing
    if read_attn || write_attn
        window = AttentionWindow(decoder_dim, A, B, N; atype=atype)
    end

    read_layer = read_attn ? ReadAttention(window) : ReadNoAttention()
    write_layer = nothing
    if write_attn
        write_layer = WriteAttention(window, decoder_dim, N; atype=atype)
    else
        write_layer = WriteNoAttention(decoder_dim, imgsize; atype=atype)
    end

    embed_layer = embed == 0 ? nothing : Embedding(nclass, embed; atype=atype)
    qnetwork = QNet(decoder_dim, noise_dim, atype)
    encoder = RNN(2*N*N+decoder_dim, encoder_dim; dataType=_etype) # FIXME: adapt to attn
    decoder = RNN(noise_dim+embed, decoder_dim; dataType=_etype)
    encoder_hidden = []
    decoder_hidden = []

    return Network(
        A,
        B,
        N,
        T,
        window,
        read_layer,
        write_layer,
        embed_layer,
        qnetwork,
        encoder,
        decoder,
        encoder_hidden,
        decoder_hidden
    )
end

function Network(N, T, encoder_dim, decoder_dim, noise_dim;
              embed=0, nclass=10, atype=_atype, read_attn=true, write_atnn=true)
    A = B = N
    Network(A, B, N, T, encoder_dim, decoder_dim, noise_dim;
         embed=embed, atype=_atype, read_attn=read_attn, write_attn=write_attn)
end

# reconstruct
function (model::Network)(x, y=nothing; cprev=_atype(zeros(size(x))))
    empty!(model.encoder_hidden)
    empty!(model.decoder_hidden)
    output = DRAWOutput()
    atype = typeof(value(model.qnetwork.mu_layer.w))

    hdec = get_hdec(model, x)
    hinit = reshape(hdec, size(hdec,1), size(hdec,2), 1)
    model.decoder.h, model.decoder.c = hinit, 0.0f0
    model.encoder.h = model.encoder.c = 0.0f0
    for t = 1:model.T
        # update xhat and then read
        xhat = x - sigm.(cprev)
        rt = model.read_layer(x, xhat, hdec)

        # encoder
        input = vcat(rt, hdec)
        input = reshape(input, size(input, 1), size(input, 2), 1)
        model.encoder(input)
        henc = model.encoder.h
        henc = reshape(henc, size(henc)[1:2])

        # qnetwork
        z, mu, logsigma, sigma = model.qnetwork(henc)

        input = z
        if model.embed_layer != nothing
            y == nothing && error("You should pass labels also.")
            input = vcat(z, model.embed_layer(y))
        end

        # decoder
        model.decoder(input)
        hdec = model.decoder.h
        hdec = reshape(hdec, size(hdec)[1:2])

        # write and update draw output
        c = cprev + model.write_layer(hdec)
        push!(output, mu, logsigma, sigma, c)
        cprev = output.cs[end]
    end
    return output
end

# generate
function (model::Network)(batchsize::Int, y=nothing)
    empty!(model.encoder_hidden)
    empty!(model.decoder_hidden)
    output = DRAWOutput()
    model.decoder.h = model.decoder.c = 0.0f0
    for t = 1:model.T
        z = sample_noise(model, batchsize)
        c = t == 1 ? 0.0f0 : output.cs[end]

        input = z
        if model.embed_layer != nothing
            y == nothing && error("You should pass labels also.")
            input = vcat(z, model.embed_layer(y))
        end

        model.decoder(input)
        hdec = model.decoder.h
        hdec = reshape(hdec, size(hdec)[1:2])
        wt = model.write_layer(hdec)
        c = c .+ wt
        push!(output, c)
    end

    for i = 1:length(output.cs)
        output.cs[i] = Array(sigm.(output.cs[i]))
    end
    return output
end

function get_hdec(model::Network, x)
    num_hidden, batchsize = model.decoder.hiddenSize, size(x, 2)
    h = _atype(zeros(num_hidden, batchsize))
end

function sample_noise(model::Network, batchsize::Int)
    return sample_noise(model.qnetwork, batchsize; generation=true)
end

function binary_cross_entropy(x, x̂)
    F = _etype
    s = @. x * log(x̂ + F(1f-8)) + (1-x) * log(1 - x̂ + F(1f-8))
    return -mean(s)
end

function loss(model::Network, x, y=nothing; loss_values=[])
    output = model(x, y)
    xhat = sigm.(output.cs[end])
    Lx = binary_cross_entropy(x, xhat) * model.A * model.B
    kl_terms = []
    for t = 1:model.T
        mu_2 = square(output.mus[t])
        sigma_2 = square(output.sigmas[t])
        logsigma = output.logsigmas[t]
        kl = 0.5 * sum((mu_2 + sigma_2-2logsigma), dims=1) .- 0.5f0 * model.T
        push!(kl_terms, kl)
    end
    kl_sum = reduce(+, kl_terms) # == sum(kl_terms)
    Lz = mean(kl_sum)
    push!(loss_values, value(Lx), value(Lz))
    return Lx + Lz
end

function init_opt!(model::Network, optimizer="Adam()")
    for par in params(model)
        par.opt = eval(Meta.parse(optimizer))
    end
end

max_stars_repo_name: ilkerkesen/draw

max_issues_repo_name: ilkerkesen/draw

max_forks_repo_name: ilkerkesen/draw


text:
module CubicHermiteSpline

using ArgCheck

export CubicHermiteSplineInterpolation
export interp, grad

struct CubicHermiteSplineInterpolation
    x
    y
    gradient
end

"""
    basis(t)

Compute the basis functions of cubic Hermite spline interpolation. They are:
    H_{00} = 2t^3 - 3t^2 + 1 = (1+2t)(t-1)^2
    H_{10} = t^3 - 2t^2 + t = t(t-1)^2
    H_{01} = -2t^3 + 3t^2 = t^2(3-2t)
    H_{11} = t^3 - t^2 = t^2(t - 1)
"""
function basis(t)
    t2 = t * t
    it = t - 1
    it2 = it * it
    tt = 2 * t
    h00 = (1 + tt) * it2
    h10 = t * it2
    h01 = t2 * (3 - tt)
    h11 = t2 * it
    return h00, h10, h01, h11
end

"""
    basis_derivative(t)

Compute the basis functions of cubic Hermite spline interpolation. They are:
    H_{00} = 6t^2 - 6t = 6t(t-1)
    H_{10} = 3t^2 - 4t + 1 = (3t-1)(t-1)
    H_{01} = -6t^2 + 6t = -6t(t-1)
    H_{11} = 3t^2 - 2t = t(3t - 2)
"""
function basis_derivative(t)
    t2 = t * t
    h00 = 6*t2 - 6*t
    h10 = 3*t2 - 4*t + 1
    h01 = -6*t2 + 6*t
    h11 = 3*t2 - 2*t
    return h00, h10, h01, h11
end

function findinterval(v, x)
    for i in eachindex(x)
        if x[i] == v
            return i, x[i], nothing
        elseif x[i] > v
            return i-1, x[i-1], x[i]
        end
    end
end

function _interp(spl::CubicHermiteSplineInterpolation, v; grad=false)
    x = spl.x
    y = spl.y
    gradient = spl.gradient
    @argcheck v >= x[1]
    @argcheck v <= x[end]

    idx, x1, x2 = findinterval(v, x)
    if x2 === nothing
       return grad ? gradient[idx] : y[idx]
    end

    # mapping (x1, x2) to (0, 1), h is the scaling constant
    t = (v - x1) / (x2 - x1)
    h = x2 - x1
    y1, y2 = y[idx], y[idx+1]
    k1, k2 = gradient[idx], gradient[idx+1]

    h00, h10, h01, h11 = grad ? basis_derivative(t) : basis(t)
    r = y1*h00 + h*k1*h10 + y2*h01 + h*k2*h11

    # Need to rescale the interpolated gradient before return
    return grad ? (r/h) : r
end

(spl::CubicHermiteSplineInterpolation)(v::Real; grad=false) = _interp(spl, v; grad=grad)
(spl::CubicHermiteSplineInterpolation)(x::AbstractVector; grad=false) = spl.(x; grad=grad)

# handy methods
interp(spl::CubicHermiteSplineInterpolation, p) = spl(p)
grad(spl::CubicHermiteSplineInterpolation, p) = spl(p; grad=true)

end # module


text:
"""
Transient vibration analysis of clamped thin plate

A homogeneous square plate from the NAFEMS Benchmark, Test No. FV12.

The plate is discretized with shell elements.

"""

using LinearAlgebra
using SparseArrays
using Arpack
using FinEtools
using FinEtoolsDeforLinear
using FinEtoolsFlexStructures.FESetShellT3Module: FESetShellT3
using FinEtoolsFlexStructures.FEMMShellT3FFModule
using FinEtoolsFlexStructures.RotUtilModule: initial_Rfield, update_rotation_field!

function solve()
    E = 200e3*phun("MPa")
    nu = 0.3;
    rho= 8000*phun("KG/M^3");
    thickness = 0.05*phun("m");
    L = 10.0*phun("m");

    n = 4*16
    @info "Mesh: $n elements per side"

    tolerance = L/n/1000
    fens, fes = T3block(L,L,n,n);
    fens.xyz[:, 1] .-= L/2
    fens.xyz[:, 2] .-= L/2
    fens.xyz = xyz3(fens)

    mater = MatDeforElastIso(DeforModelRed3D, rho, E, nu, 0.0)

    sfes = FESetShellT3()
    accepttodelegate(fes, sfes)
    femm = FEMMShellT3FFModule.make(IntegDomain(fes, TriRule(1), thickness), mater)

# Construct the requisite fields, geometry and displacement
# Initialize configuration variables
    geom0 = NodalField(fens.xyz)
    u0 = NodalField(zeros(size(fens.xyz,1), 3))
    Rfield0 = initial_Rfield(fens)
    dchi = NodalField(zeros(size(fens.xyz,1), 6))

# Apply EBC's
    l1 = connectednodes(meshboundary(fes))
    for i in 1:6
        setebc!(dchi, l1, true, i)
    end
    applyebc!(dchi)
    numberdofs!(dchi);

# Assemble the system matrix
    FEMMShellT3FFModule.associategeometry!(femm, geom0)
    K = FEMMShellT3FFModule.stiffness(femm, geom0, u0, Rfield0, dchi);
    M = FEMMShellT3FFModule.mass(femm, geom0, dchi);
    # Check that the mass matrix is diagonal
    I, J, V = findnz(M)
    @assert length(I) == dchi.nfreedofs

    evals, evecs, nconv = eigs(K, M; nev=1, which=:LM, explicittransform=:none)
    @show typeof(evals), evals
    @show omega_max = sqrt(evals[1])
    # omega_max =   2.764502757286571e+06
    @show dt = Float64(0.99* 2/omega_max)
    @show typeof(dt)

    U0 = gathersysvec(dchi)
    U1 = deepcopy(U0)
    V0 = deepcopy(U0)
    V1 = deepcopy(U0)
    A0 = deepcopy(U0)
    A1 = deepcopy(U0)
    F = deepcopy(U0)
    E = deepcopy(U0)
    invM = deepcopy(U0)
    invM .= 1.0 ./ vec(diag(M))

    nsteps = 10000
    nbtw = 1000
    t = 0.0
    @time for step in 1:nsteps
        @. U1 = U0 + dt*V0 + ((dt^2)/2)*A0;
    end

end

solve()

nothing


text:
function setdefault(d::Associative, key, value)
# set d[key]=value if not present; return d[key]
    if !has(d, key);  d[key] = value;  end
    return d[key]
end

function update(dest::Associative, source)
# add entries to dest, overwriting existing ones
    for (key, value) in source
        dest[key] = value
    end
end

function makedict{K,V}(::Type{K}, ::Type{V}, sources...)
# make a HashTable{K,V} out of an Associative/(key,value) iterable
    d = HashTable{K,V}()
    for source in sources
        update(d, source)
    end
    d
end

# copy HashTable entries, but don't deep copy the keys/values
copydict{K,V}(d::HashTable{K,V}) = makedict(K,V,d)

#convert{K,V}(::Type{HashTable{K,V}}, d::Associative) = copydict{K,V}(d)

quote_sym(sym::Symbol) = expr(:quote, sym)

macro dict(args...)
# Create a HashTable{Symbol}
# Usage:
#     d = @dict x=1 s="hello"
# is equivalent to
#     d = {:x => 1, :s => "hello"}
# except that @dict throws an error on duplicate keys

    names = Set{Symbol}()
    toentry(arg) = error("@dict: malformed argument $arg")
    function toentry(ex::Expr)
        if ex.head != :(=); error("@dict: not an assignment: $ex"); end

        # check name duplicates
        name = ex.args[1]
        if has(names, name); error("@dict: duplicate definition of $name"); end
        add(names, name)

        expr(:(=>), quote_sym(name), ex.args[2])
    end
    expr(:cell1d, { toentry(arg) | arg in args } )
end

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
import CLBLAS
using OpenCL: cl
const clblas = CLBLAS
clblas.setup()

device, ctx, queue = clblas.get_next_compute_context()

N = Unsigned(5)
data = Vector{Float32}(5)
data[1] = 1.1
data[2] = 2.22
data[3] = 3.333
data[4] = 4.4444
data[5] = 5.55555

dat = Vector{Float32}(5)
dat[1] = 6.1
dat[2] = 7.22
dat[3] = 8.333
dat[4] = 9.4444
dat[5] = 10.55555

bufX = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=data)
bufY = cl.Buffer(Float32, ctx, (:r, :copy),hostbuf=dat)

bufDotP = cl.Buffer(Float32, ctx, :w, 1)
scratchBuff = cl.Buffer(Float32, ctx, :rw, length(data))

offDP = offy = offx = Unsigned(0)
incx = cl.cl_int(1)
incy = cl.cl_int(1)
ncq = cl.cl_uint(1)
clq = queue.id
newl = cl.cl_uint(0)

event = cl.UserEvent(ctx, retain=true)
ptrEvent = [event.id]

clblas.clblasSdot( N, bufDotP.id, offDP, bufX.id, offx, incx, bufY.id, offy, incy, scratchBuff.id, ncq, [clq], newl, C_NULL, ptrEvent);

cl.api.clWaitForEvents(cl.cl_uint(1), ptrEvent)

result = Vector{Float32}(length(bufDotP))
cl.enqueue_read_buffer(queue, bufDotP, result, Unsigned(0), nothing, true)

expected = dot(data, dat)
println(result)
println("------------")
println(expected)

try
    # locally passes without atol
    @assert(isapprox(norm(expected - result[1]), zero(Float32), atol=0.00002))
    info("success!")
finally
    clblas.teardown()
end


text:
"""Type class which represents a tagged dictionary

Tagged dictionaries are used to represent complex Julia types. Using a struct
instead of an actual Dictionary requires less memory allocation and allows us
to use multiple dispatch on the resulting tree structure.

It inherits abstract dict just for show."""
abstract type Tagged <: AbstractDict{Symbol, Any} end

"Type class for types which can occupuy the 'type' field in a struct"
abstract type TaggedStructType <: Tagged end

# Needed for show
Base.length(tt::T) where {T <: Tagged} = length(fieldnames(T))
Base.getindex(tt::Tagged, k::Symbol) = getfield(tt, k)
Base.iterate(tt::T) where {T <: Tagged} = let first = fieldnames(T)[1]
  (first => tt[first], 2)
end
Base.iterate(tt::T, state) where {T <: Tagged} = let names = fieldnames(T)
  if length(names) >= state
    next = names[state]
    (next => tt[next], state + 1)
  else
    nothing
  end
end

function applychildren!(f::Function, tt::T)::T where {T <: Tagged}
  for fn in fieldnames(T)
    setfield!(tt, fn, f(getfield(tt, fn)))
  end
  tt
end

struct TaggedBackref <: TaggedStructType
  ref::Int64
end

Base.show(io::IO, br::TaggedBackref) = print(io, "Ref(", br.ref, ")")

mutable struct TaggedTuple <: Tagged
  data::BSONArray
end

mutable struct TaggedSvec <: Tagged
  data::BSONArray
end

const TaggedParam = Union{Tagged, TypeVar, Type}

mutable struct TaggedType <: TaggedStructType
  name::Vector{String}
  params::Vector{TaggedParam}
end

mutable struct TaggedStruct <: TaggedStructType
  ttype::TaggedStructType
  data::Union{BSONArray, Vector{UInt8}}
end

mutable struct TaggedArray <: Tagged
  ttype::TaggedStructType
  size::BSONArray
  data::Union{Vector{UInt8}, BSONArray}
end

mutable struct TaggedAnonymous <: TaggedStructType
  typename::Union{TaggedStruct, TaggedBackref}
  params::Vector{TaggedParam}
end

struct TaggedRef <: Tagged
  path::Vector{String}
end

mutable struct TaggedUnionall <: Tagged
  var::Union{TaggedStruct, TaggedBackref, TaggedUnionall}
  body::Union{TaggedType, TaggedBackref, TaggedUnionall, TaggedAnonymous}
end

struct BackRefsWrapper
  root::Union{BSONDict, TaggedStruct}
  refs::BSONArray
end

function Base.show(io::IO, brw::BackRefsWrapper)
  summary(io, brw)
  print(io, ".root => ")
  show(io, MIME("text/plain"), brw.root)
  summary(io, brw)
  print(io, ".refs => ")
  show(io, MIME("text/plain"), brw.refs)
end

function applychildren!(f::Function, tt::Union{TaggedTuple, TaggedSvec})
  tt.data = f(tt.data)
  tt
end
applychildren!(f::Function, v::Vector{TaggedParam}) = applyvec!(f, v)
function applychildren!(f::Function, tt::TaggedType)::TaggedType
  tt.params = f(tt.params)
  tt
end
function applychildren!(f::Function, ts::TaggedStruct)::TaggedStruct
  ts.ttype = f(ts.ttype)
  if ts.data != nothing
    ts.data = f(ts.data)
  end
  ts
end
applychildren!(::Function, tr::TaggedRef) = tr

const PassthroughTypes = Union{Vector{UInt8}, Type{Union{}}, Symbol, String}
const RefVectorTypes = Union{Vector{TaggedParam}, BSONArray}

_raise_recursive(x::T, cache::IdDict{Any, Any}) where T = if isbitstype(T)
  x
elseif haskey(cache, x)
  cache[x]
else
  cache[x] = raise_recursive(x, cache)
end

raise_recursive(x, ::IdDict{Any, Any}) = x
function raise_recursive(d::BSONDict, cache::IdDict{Any, Any})::BSONDict
  cache[d] = d
  applychildren!(x -> _raise_recursive(x, cache), d)
end
function raise_recursive(v::T, cache::IdDict{Any, Any}) where {T <: RefVectorTypes}
  cache[v] = v
  applyvec!(x -> _raise_recursive(x, cache), v)
end
raise_recursive(tt::TaggedTuple, cache::IdDict{Any, Any})::Tuple =
  (raise_recursive(tt.data, cache)...,)
raise_recursive(tt::TaggedSvec, cache::IdDict{Any, Any})::SimpleVector =
  Core.svec(raise_recursive(tt.data, cache)...)
raise_recursive(tt::TaggedType, cache::IdDict{Any, Any})::Type =
  constructtype(resolve(tt.name), raise_recursive(tt.params, cache))
function raise_recursive(ts::TaggedStruct, cache::IdDict{Any, Any})
  T::Type = _raise_recursive(ts.ttype, cache)

  if ismutable(T)
    return newstruct_raw(cache, T, ts)
  end

  data = raise_recursive(ts.data, cache)
  if isprimitive(T)
    newprimitive(T, data)
  else
    newstruct(T, data...)
  end
end
function raise_recursive(ta::TaggedArray, cache::IdDict{Any, Any})
  T::DataType = _raise_recursive(ta.ttype, cache)
  size = raise_recursive(ta.size, cache)
  data() = raise_recursive(ta.data, cache)

  if isbitstype(T)
    if sizeof(T) == 0
      fill(T(), size...)
    else
      reshape(reinterpret_(T, data()), size...)
    end
  else
    Array{T}(reshape(data(), size...))
  end
end
function raise_recursive(ts::TaggedAnonymous, cache::IdDict{Any, Any})
  tn = _raise_recursive(ts.typename::TaggedStruct, cache)
  pr = _raise_recursive(ts.params, cache)
  constructtype(tn.wrapper, pr)
end
raise_recursive(ts::TaggedRef, cache::IdDict{Any, Any})::Module =
  resolve(ts.path)
raise_recursive(tu::TaggedUnionall, cache::IdDict{Any, Any})::UnionAll =
  UnionAll(_raise_recursive(tu.var, cache),
           _raise_recursive(tu.body, cache))


text:
# Variable (with applied transformations following the CF convention)
mutable struct CFVariable{T,N,TV,TA,TSA}  <: AbstractArray{T, N}
    # this var is generally a `Variable` type
    var::TV
    # Dict-like object for all attributes read from disk
    attrib::TA
    # a named tuple with fill value, scale factor, offset,...
    # immutable for type-stability
    _storage_attrib::TSA
end

NCDataset(var::CFVariable) = NCDataset(var.var)

"""
    sz = size(var::CFVariable)

Return a tuple of integers with the size of the variable `var`.

!!! note

    Note that the size of a variable can change, i.e. for a variable with an
    unlimited dimension.
"""
Base.size(v::CFVariable) = size(v.var)

############################################################
# Creating variables
############################################################

"""
    defVar(ds::NCDataset,name,vtype,dimnames; kwargs...)
    defVar(ds::NCDataset,name,data,dimnames; kwargs...)

Define a variable with the name `name` in the dataset `ds`.  `vtype` can be
Julia types in the table below (with the corresponding NetCDF type). The
parameter `dimnames` is a tuple with the names of the dimension.  For scalar
this parameter is the empty tuple `()`.
The variable is returned (of the type CFVariable).

Instead of providing the variable type one can directly give also the data `data` which
will be used to fill the NetCDF variable. In this case, the dimensions with
the appropriate size will be created as required using the names in `dimnames`.

If `data` is a vector or array of `DateTime` objects, then the dates
are saved as double-precision floats and units
"$(CFTime.DEFAULT_TIME_UNITS)" (unless a time unit
is specifed with the `attrib` keyword as described below). Dates are
converted to the default calendar in the CF conversion which is the
mixed Julian/Gregorian calendar.

## Keyword arguments

* `fillvalue`: A value filled in the NetCDF file to indicate missing data.
   It will be stored in the _FillValue attribute.
* `chunksizes`: Vector integers setting the chunk size. The total size of a chunk must be less than 4 GiB.
* `deflatelevel`: Compression level: 0 (default) means no compression and 9 means maximum compression. Each chunk will be compressed individually.
* `shuffle`: If true, the shuffle filter is activated which can improve the compression ratio.
* `checksum`: The checksum method can be `:fletcher32` or `:nochecksum` (checksumming is disabled, which is the default)
* `attrib`: An iterable of attribute name and attribute value pairs, for example a `Dict`, `DataStructures.OrderedDict` or simply a vector of pairs (see example below)
* `typename` (string): The name of the NetCDF type required for [vlen arrays](https://web.archive.org/save/https://www.unidata.ucar.edu/software/netcdf/netcdf-4/newdocs/netcdf-c/nc_005fdef_005fvlen.html)

`chunksizes`, `deflatelevel`, `shuffle` and `checksum` can only be
set on NetCDF 4 files.

## NetCDF data types

| NetCDF Type | Julia Type |
|-------------|------------|
| NC_BYTE     | Int8 |
| NC_UBYTE    | UInt8 |
| NC_SHORT    | Int16 |
| NC_INT      | Int32 |
| NC_INT64    | Int64 |
| NC_FLOAT    | Float32 |
| NC_DOUBLE   | Float64 |
| NC_CHAR     | Char |
| NC_STRING   | String |

## Example:

In this example, `scale_factor` and `add_offset` are applied when the `data`
is saved.

```julia-repl
julia> using DataStructures
julia> data = randn(3,5)
julia> NCDataset("test_file.nc","c") do ds
          defVar(ds,"temp",data,("lon","lat"), attrib = OrderedDict(
             "units" => "degree_Celsius",
             "add_offset" => -273.15,
             "scale_factor" => 0.1,
             "long_name" => "Temperature"
          ))
       end;
```

!!! note

    If the attributes `_FillValue`, `add_offset`, `scale_factor`, `units` and
    `calendar` are used, they should be defined when calling `defVar` by using the
    parameter `attrib` as shown in the example above.
"""
function defVar(ds::NCDataset,name,vtype::DataType,dimnames; kwargs...)
    # all keyword arguments as dictionary
    kw = Dict(k => v for (k,v) in kwargs)

    defmode(ds) # make sure that the file is in define mode
    dimids = Cint[nc_inq_dimid(ds.ncid,dimname) for dimname in dimnames[end:-1:1]]

    typeid =
        if vtype <: Vector
            # variable-length type
            typeid = nc_def_vlen(ds.ncid, kw[:typename], ncType[eltype(vtype)])
        else
            # base-type
            ncType[vtype]
        end

    varid = nc_def_var(ds.ncid,name,typeid,dimids)

    if haskey(kw,:chunksizes)
        storage = :chunked
        chunksizes = kw[:chunksizes]

        # this will fail on NetCDF-3 files
        nc_def_var_chunking(ds.ncid,varid,storage,reverse(chunksizes))
    end

    if haskey(kw,:shuffle) || haskey(kw,:deflatelevel)
        shuffle = get(kw,:shuffle,false)
        deflate = haskey(kw,:deflatelevel)
        deflate_level = get(kw,:deflatelevel,0)

        # this will fail on NetCDF-3 files
        nc_def_var_deflate(ds.ncid,varid,shuffle,deflate,deflate_level)
    end

    if haskey(kw,:checksum)
        checksum = kw[:checksum]
        nc_def_var_fletcher32(ds.ncid,varid,checksum)
    end

    if haskey(kw,:fillvalue)
        fillvalue = kw[:fillvalue]
        nofill = get(kw,:nofill,false)
        nc_def_var_fill(ds.ncid, varid, nofill, vtype(fillvalue))
    end

    if haskey(kw,:attrib)
        v = ds[name]
        for (attname,attval) in kw[:attrib]
            v.attrib[attname] = attval
        end
    end

    return ds[name]
end

# data has the type e.g. Array{Union{Missing,Float64},3}
function defVar(ds::NCDataset,
                name,
                data::AbstractArray{Union{Missing,nctype},N},
                dimnames;
                kwargs...) where nctype <: Union{Int8,UInt8,Int16,Int32,Int64,Float32,Float64} where N
    _defVar(ds::NCDataset,name,data,nctype,dimnames; kwargs...)
end

# data has the type e.g. Vector{DateTime}, Array{Union{Missing,DateTime},3} or
# Vector{DateTime360Day}
# Data is always stored as Float64 in the NetCDF file
function defVar(ds::NCDataset,
                name,
                data::AbstractArray{<:Union{Missing,nctype},N},
                dimnames;
                kwargs...) where nctype <: Union{DateTime,AbstractCFDateTime} where N
    _defVar(ds::NCDataset,name,data,Float64,dimnames; kwargs...)
end

function defVar(ds::NCDataset,name,data,dimnames; kwargs...)
    nctype = eltype(data)
    _defVar(ds::NCDataset,name,data,nctype,dimnames; kwargs...)
end

function _defVar(ds::NCDataset,name,data,nctype,dimnames; attrib = [], kwargs...)
    # define the dimensions if necessary
    for i = 1:length(dimnames)
        if !(dimnames[i] in ds.dim)
            ds.dim[dimnames[i]] = size(data,i)
        elseif !(dimnames[i] in unlimited(ds.dim))
            dimlen = ds.dim[dimnames[i]]

            if (dimlen !== size(data,i))
                throw(NetCDFError(
                    -1,"dimension $(dimnames[i]) is already defined with the " *
                    "length $dimlen. It cannot be redefined with a length of $(size(data,i))."))
            end
        end
    end

    T = eltype(data)
    attrib = collect(attrib)

    if T <: Union{TimeType,Missing}
        dattrib = Dict(attrib)
        if !haskey(dattrib,"units")
            push!(attrib,"units" => CFTime.DEFAULT_TIME_UNITS)
        end
        if !haskey(dattrib,"calendar")
            # these dates cannot be converted to the standard calendar
            if T <: Union{DateTime360Day,Missing}
                push!(attrib,"calendar" => "360_day")
            elseif T <: Union{DateTimeNoLeap,Missing}
                push!(attrib,"calendar" => "365_day")
            elseif T <: Union{DateTimeAllLeap,Missing}
                push!(attrib,"calendar" => "366_day")
            end
        end
    end

    v =
        if Missing <: T
            # make sure a fill value is set (it might be overwritten by kwargs...)
            defVar(ds,name,nctype,dimnames;
                   fillvalue = fillvalue(nctype),
                   attrib = attrib,
                   kwargs...)
        else
            defVar(ds,name,nctype,dimnames;
                   attrib = attrib,
                   kwargs...)
        end

    v[:] = data
    return v
end

function defVar(ds::NCDataset,name,data::T; kwargs...) where T <: Number
    v = defVar(ds,name,T,(); kwargs...)
    v[:] = data
    return v
end
export defVar

function _boundsParentVar(ds,name)
    # get from cache is available
    if ds._boundsmap !== nothing
        return get(ds._boundsmap,name,nothing)
    else
        for vname in keys(ds)
            v = variable(ds,vname)
            bounds = get(v.attrib,"bounds",nothing)
            if bounds == name
                return vname
            end
        end

        return nothing
    end
end

"""
    calendar,time_origin,time_factor = calendar_time(
        attrib;
        calendar = nothing,
        time_origin = nothing,
        time_factor = nothing)
        )

"""
function _calendar_time(
    attrib;
    calendar = nothing,
    time_origin = nothing,
    time_factor = nothing,
    )

    if haskey(attrib,"units")
        units = attrib["units"]
        if (units isa String) && occursin(" since ",units)
            calendar = lowercase(get(attrib,"calendar","standard"))
            try
                time_origin,time_factor = CFTime.timeunits(units, calendar)
            catch
                # keep defaults
            end
        end
    end

    return calendar,time_origin,time_factor
end

"""
    v = getindex(ds::NCDataset,varname::AbstractString)

Return the NetCDF variable `varname` in the dataset `ds` as a
`NCDataset.CFVariable`. The CF convention are honored when the
variable is indexed:
* `_FillValue` will be returned as `missing`
* `scale_factor` and `add_offset` are applied (output = `scale_factor` * data_in_file +  `add_offset`)
* time variables (recognized by the units attribute and possibly the calendar attribute) are returned usually as
  `DateTime` object. Note that `DateTimeAllLeap`, `DateTimeNoLeap` and
  `DateTime360Day` cannot be converted to the proleptic gregorian calendar used in
  julia and are returned as such. If a calendar is defined but not among the
ones specified in the CF convention, then the data in the NetCDF file is not
converted into a date structure.

Note that the attribute `missing_value` is not used to determine which
element is `missing` only `_FillValue`.

A call `getindex(ds,varname)` is usually written as `ds[varname]`.

If variable represents a cell boundary, the attributes `calendar` and `units` of the related NetCDF variables are used, if they are not specified. For example:

```
dimensions:
  time = UNLIMITED; // (5 currently)
  nv = 2;
variables:
  double time(time);
    time:long_name = "time";
    time:units = "hours since 1998-04-019 06:00:00";
    time:bounds = "time_bnds";
  double time_bnds(time,nv);
```

In this case, the variable `time_bnds` uses the units and calendar of `time`
because both variables are related thought the bounds attribute following the CF conventions.
"""
function Base.getindex(ds::AbstractDataset,varname::SymbolOrString)
    v = variable(ds,varname)

    fillvalue = get(v.attrib,"_FillValue",nothing)
    scale_factor = get(v.attrib,"scale_factor",nothing)
    add_offset = get(v.attrib,"add_offset",nothing)

    # special case for bounds variable who inherit
    # units and calendar from parent variables
    parentname = _boundsParentVar(ds,varname)

    if parentname == nothing
        calendar = nothing
        time_origin = nothing
        time_factor = nothing
    else
        calendar,time_origin,time_factor = _calendar_time(variable(ds,parentname).attrib)
    end

    calendar,time_origin,time_factor = _calendar_time(
        v.attrib;
        calendar = calendar,
        time_origin = time_origin,
        time_factor = time_factor)

    scaledtype = eltype(v)

    if eltype(v) <: Number
        if scale_factor != nothing
            scaledtype = promote_type(scaledtype, typeof(scale_factor))
        end

        if add_offset != nothing
            scaledtype = promote_type(scaledtype, typeof(add_offset))
        end
    end

    rettype = scaledtype

    # rettype can be a date if calendar is different from nothing
    if calendar != nothing
        DT = nothing
        try
            DT = CFTime.timetype(calendar)
        catch
        end

        if DT !== nothing
            # this is the only supported option for NCDatasets
            prefer_datetime = true

            if prefer_datetime &&
                (DT in [DateTimeStandard,DateTimeProlepticGregorian,DateTimeJulian])
                rettype = DateTime
            else
                rettype = DT
            end
        else
            @warn("unsupported calendar `$calendar`. Time units are ignored.")
        end
    end

    if fillvalue != nothing
        rettype = Union{Missing,rettype}
    end

    storage_attrib = (
        fillvalue = fillvalue,
        scale_factor = scale_factor,
        add_offset = add_offset,
        calendar = calendar,
        time_origin = time_origin,
        time_factor = time_factor,
    )

    return CFVariable{rettype,ndims(v),typeof(v),typeof(v.attrib),typeof(storage_attrib)}(
        v,v.attrib,storage_attrib)
end

"""
    dimnames(v::CFVariable)

Return a tuple of strings with the dimension names of the variable `v`.
"""
dimnames(v::CFVariable)  = dimnames(v.var)

"""
    dimsize(v::CFVariable)
Get the size of a `CFVariable` as a named tuple of dimension → length.
"""
function dimsize(v::CFVariable)
    s = size(v)
    names = Symbol.(dimnames(v))
    return NamedTuple{names}(s)
end
export dimsize

name(v::CFVariable) = name(v.var)
chunking(v::CFVariable,storage,chunksize) = chunking(v.var,storage,chunksize)
chunking(v::CFVariable) = chunking(v.var)

deflate(v::CFVariable,shuffle,dodeflate,deflate_level) = deflate(v.var,shuffle,dodeflate,deflate_level)
deflate(v::CFVariable) = deflate(v.var)

checksum(v::CFVariable,checksummethod) = checksum(v.var,checksummethod)
checksum(v::CFVariable) = checksum(v.var)

fillmode(v::CFVariable) = fillmode(v.var)

fillvalue(v::CFVariable) = v._storage_attrib.fillvalue
scale_factor(v::CFVariable) = v._storage_attrib.scale_factor
add_offset(v::CFVariable) = v._storage_attrib.add_offset
time_origin(v::CFVariable) = v._storage_attrib.time_origin
calendar(v::CFVariable) = v._storage_attrib.calendar
""""
    tf = time_factor(v::CFVariable)

The time unit in milliseconds. E.g. seconds would be 1000., days would be 86400000.
The result can also be `nothing` if the variable has no time units.
"""
time_factor(v::CFVariable) = v._storage_attrib.time_factor

############################################################
# CFVariable
############################################################

# fillvalue can be NaN (unfortunately)
@inline isfillvalue(data,fillvalue) = data == fillvalue
@inline isfillvalue(data,fillvalue::AbstractFloat) = (isnan(fillvalue) ? isnan(data) : data == fillvalue)

@inline CFtransform_missing(data,fv) = (isfillvalue(data,fv) ? missing : data)
@inline CFtransform_missing(data,fv::Nothing) = data

@inline CFtransform_replace_missing(data,fv) = (ismissing(data) ? fv : data)
@inline CFtransform_replace_missing(data,fv::Nothing) = data

@inline CFtransform_scale(data,scale_factor) = data*scale_factor
@inline CFtransform_scale(data,scale_factor::Nothing) = data
@inline CFtransform_scale(data::T,scale_factor) where T <: Union{Char,String} = data
@inline CFtransform_scale(data::T,scale_factor::Nothing) where T <: Union{Char,String} = data

@inline CFtransform_offset(data,add_offset) = data + add_offset
@inline CFtransform_offset(data,add_offset::Nothing) = data
@inline CFtransform_offset(data::T,add_factor) where T <: Union{Char,String} = data
@inline CFtransform_offset(data::T,add_factor::Nothing) where T <: Union{Char,String} = data

@inline asdate(data::Missing,time_origin,time_factor,DTcast) = data
@inline asdate(data,time_origin::Nothing,time_factor,DTcast) = data
@inline asdate(data::Missing,time_origin::Nothing,time_factor,DTcast) = data
@inline asdate(data,time_origin,time_factor,DTcast) =
    convert(DTcast,time_origin + Dates.Millisecond(round(Int64,time_factor * data)))

@inline fromdate(data::TimeType,time_origin,inv_time_factor) =
    Dates.value(data - time_origin) * inv_time_factor
@inline fromdate(data,time_origin,time_factor) = data

@inline function CFtransform(data,fv,scale_factor,add_offset,time_origin,time_factor,DTcast)
    return asdate(
        CFtransform_offset(
            CFtransform_scale(
                CFtransform_missing(data,fv),
                scale_factor),
            add_offset),
        time_origin,time_factor,DTcast)
end

# round float to integers
_approximate(::Type{T},data) where T <: Integer = round(T,data)
_approximate(::Type,data) = data

@inline function CFinvtransform(data,fv,inv_scale_factor,minus_offset,time_origin,inv_time_factor,DT)
    return _approximate(
        DT,
        CFtransform_replace_missing(
            CFtransform_scale(
                CFtransform_offset(
                    fromdate(data,time_origin,inv_time_factor),
                    minus_offset),
                inv_scale_factor),
            fv))
end

# this is really slow
# https://github.com/JuliaLang/julia/issues/28126
#@inline CFtransformdata(data,fv,scale_factor,add_offset,time_origin,time_factor,DTcast) =
#    # in boardcasting we trust..., or not
#    CFtransform.(data,fv,scale_factor,add_offset,time_origin,time_factor,DTcast)

# for scalars
@inline CFtransformdata(data,fv,scale_factor,add_offset,time_origin,time_factor,DTcast) =
    CFtransform(data,fv,scale_factor,add_offset,time_origin,time_factor,DTcast)

# for arrays
@inline function CFtransformdata(data::AbstractArray{T,N},fv,scale_factor,add_offset,time_origin,time_factor,DTcast) where {T,N}
    out = Array{DTcast,N}(undef,size(data))
    @inbounds @simd for i in eachindex(data)
        out[i] = CFtransform(data[i],fv,scale_factor,add_offset,time_origin,time_factor,DTcast)
    end
    return out
end

@inline _inv(x::Nothing) = nothing
@inline _inv(x) = 1/x
@inline _minus(x::Nothing) = nothing
@inline _minus(x) = -x

# # so slow
# @inline function CFinvtransformdata(data,fv,scale_factor,add_offset,time_origin,time_factor,DT)
#     inv_scale_factor = _inv(scale_factor)
#     minus_offset = _minus(add_offset)
#     inv_time_factor = _inv(time_factor)
#     return CFinvtransform.(data,fv,inv_scale_factor,minus_offset,time_origin,inv_time_factor,DT)
# end

# for arrays
@inline function CFinvtransformdata(data::AbstractArray{T,N},fv,scale_factor,add_offset,time_origin,time_factor,DT) where {T,N}
    inv_scale_factor = _inv(scale_factor)
    minus_offset = _minus(add_offset)
    inv_time_factor = _inv(time_factor)

    out = Array{DT,N}(undef,size(data))
    @inbounds @simd for i in eachindex(data)
        out[i] = CFinvtransform(data[i],fv,inv_scale_factor,minus_offset,time_origin,inv_time_factor,DT)
    end
    return out
end

# for scalar
@inline function CFinvtransformdata(data,fv,scale_factor,add_offset,time_origin,time_factor,DT)
    inv_scale_factor = _inv(scale_factor)
    minus_offset = _minus(add_offset)
    inv_time_factor = _inv(time_factor)

    return CFinvtransform(data,fv,inv_scale_factor,minus_offset,time_origin,inv_time_factor,DT)
end

# this function is necessary to avoid "iterating" over a single character in Julia 1.0 (fixed Julia 1.3)
# https://discourse.julialang.org/t/broadcasting-and-single-characters/16836
@inline CFtransformdata(data::Char,fv,scale_factor,add_offset,time_origin,time_factor,DTcast) = CFtransform_missing(data,fv)
@inline CFinvtransformdata(data::Char,fv,scale_factor,add_offset,time_origin,time_factor,DT) = CFtransform_replace_missing(data,fv)

function Base.getindex(v::CFVariable,
                       indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    data = v.var[indexes...]
    return CFtransformdata(data,fillvalue(v),scale_factor(v),add_offset(v),
                           time_origin(v),time_factor(v),eltype(v))
end

function Base.setindex!(v::CFVariable,data::Array{Missing,N},indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...) where N
    v.var[indexes...] = fill(fillvalue(v),size(data))
end

function Base.setindex!(v::CFVariable,data::Missing,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    v.var[indexes...] = fillvalue(v)
end

function Base.setindex!(v::CFVariable,data::Union{T,Array{T,N}},indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...) where N where T <: Union{AbstractCFDateTime,DateTime,Union{Missing,DateTime,AbstractCFDateTime}}

    if calendar(v) !== nothing
        # can throw an convertion error if calendar attribute already exists and
        # is incompatible with the provided data
        v.var[indexes...] = CFinvtransformdata(
            data,fillvalue(v),scale_factor(v),add_offset(v),
            time_origin(v),time_factor(v),eltype(v.var))
        return data
    end

    throw(NetCDFError(-1, "Time units and calendar must be defined during defVar and cannot change"))
end

function Base.setindex!(v::CFVariable,data,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    v.var[indexes...] = CFinvtransformdata(
        data,fillvalue(v),scale_factor(v),add_offset(v),
        time_origin(v),time_factor(v),eltype(v.var))

    return data
end

############################################################
# Convertion to array
############################################################

Base.Array(v::Union{CFVariable{T,N},Variable{T,N}}) where {T,N} = v[ntuple(i -> :, Val(N))...]

Base.show(io::IO,v::CFVariable; indent="") = Base.show(io::IO,v.var; indent=indent)

# necessary for IJulia if showing a variable from a closed file
Base.show(io::IO,::MIME"text/plain",v::Union{Variable,CFVariable}) = show(io,v)

Base.display(v::Union{Variable,CFVariable}) = show(stdout,v)


text:
# This file is a part of EasyPkg, licensed under the MIT License (MIT).

__precompile__()

module EasyPkg

include("include_helpers.jl")

EasyPkg.include_sources(
    "src_handling",
)

EasyPkg.include_all_sources("testing")

end # module


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule cargo_license_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("cargo_license")
JLLWrappers.@generate_main_file("cargo_license", UUID("e85b2321-5bf6-5335-9209-779fed322386"))
end  # module cargo_license_jll


text:
using Distributed
@show myid() nworkers()
@everywhere function sampling()
    return rand()^2 + rand()^2 < 1
end

function main()
    N_SAMPLE = 1000000
    area = @time @distributed (+) for i in 1:N_SAMPLE
        sampling()
    end
    @show 4area/N_SAMPLE
end

main()

max_stars_repo_name: Soyukke/julia-paralell-examples

max_issues_repo_name: Soyukke/julia-paralell-examples

max_forks_repo_name: Soyukke/julia-paralell-examples


text:
module DataDrivenDiffEq

using LinearAlgebra
using ModelingToolkit
using QuadGK
using Statistics
using DSP
using Compat

abstract type abstractBasis end;
abstract type abstractKoopmanOperator end;

include("./optimisers/Optimise.jl")
using .Optimise
export set_threshold!, set_threshold
export STRRidge, ADMM, SR3
export ADM

include("./basis.jl")
export Basis
export variables, jacobian, dynamics
export free_parameters, parameters, variables

include("./exact_dmd.jl")
export ExactDMD
export eigen, eigvals, eigvecs
export modes, frequencies, isstable
export dynamics, update!

include("./extended_dmd.jl")
export ExtendedDMD
export dynamics, linear_dynamics
export reduce_basis, update!

include("./dmdc.jl")
export DMDc
export eigen, eigvals, eigvecs
export get_dynamics, get_input_map, dynamics

include("./sindy.jl")
export SInDy
export sparse_regression, sparse_regression!

include("./isindy.jl")
export ISInDy

include("./utils.jl")
export AIC, AICC, BIC
export hankel, optimal_shrinkage, optimal_shrinkage!
export savitzky_golay
export burst_sampling, subsample

end # module


text:
using OrdinaryDiffEq
using DiffEqCallbacks
using DataInterpolations
using Distributions
using DiffEqFlux, Flux
using Random
using Plots;

Random.seed!(123);

function sira_ode(u,p,t)
    (S,I,C) = u
    (β,γ,α) = p
    dS = -β*S*(I^α)
    dI = β*S*(I^α) - γ*I
    dC = β*S*(I^α)
    [dS,dI,dC]
end;

solver = ROS34PW3();

N = 1000.0
p = [0.5, 0.25, 0.9]
u0 = [0.99, 0.01, 0.0]
tspan = (0., 40.)
δt = 1;

sira_prob = ODEProblem(sira_ode, u0, tspan, p)
sira_sol = solve(sira_prob, solver, saveat = δt);

train_time = 30.0
tsdata = Array(sira_sol(0:δt:train_time))
cdata = diff(tsdata[3,:])
noisy_data = rand.(Poisson.(N .* cdata));

tt = 0:δt:train_time
plot(tt[2:end],
     N .* cdata,
     xlabel = "Time",
     ylabel = "Number of new infected",
     label = "Model")
scatter!(tt,
         noisy_data,
         label = "Simulated data")

function sir_ude(u,p_,t,foi)
    S,I,C = u
    β,γ,α = p
    λ = foi([I],p_)[1]
    dS = -λ*S
    dI = λ*S - γ*I
    dC = λ*S
    [dS, dI, dC]
end;

function foi(ivec,p)
    t = 0:0.1:1
    f = LinearInterpolation([0.0;exp.(p)],t)
    return [f(ivec[1])]
end
p_ = log.(0.6 .* collect(0.1:0.1:1));

sir_psm = (u,p_,t) -> sir_ude(u,p_,t,foi)
prob_psm = ODEProblem(sir_psm,
                      u0,
                      (0.0, train_time),
                      p_);

function predict(θ, prob)
    Array(solve(prob,
                solver;
                u0 = u0,
                p = θ,
                saveat = δt,
                sensealg = InterpolatingAdjoint(autojacvec=ReverseDiffVJP())))
end;

function loss(θ, prob)
    pred = predict(θ, prob)
    cpred = abs.(N*diff(pred[3,:]))
    Flux.poisson_loss(cpred, float.(noisy_data)), cpred
end;

const losses = []
callback = function (p, l, pred)
    push!(losses, l)
    numloss = length(losses)
    if numloss % 20 == 0
        display("Epoch: " * string(numloss) * " Loss: " * string(l))
    end
    return false
end;

res_psm = DiffEqFlux.sciml_train((θ)->loss(θ,prob_psm),
                                  p_,
                                  cb = callback);

prob_psm_fit = ODEProblem(sir_psm, u0, tspan, res_psm.minimizer)
sol_psm_fit = solve(prob_psm_fit, solver, saveat = δt)
scatter(sira_sol, label=["True Susceptible" "True Infected" "True Recovered"],title="Fitted partially specified model")
plot!(sol_psm_fit, label=["Estimated Susceptible" "Estimated Infected" "Estimated Recovered"])
Plots.vline!([train_time],label="Training time")

Imax = maximum(tsdata[2,:])
Igrid = 0:0.01:1.0 # create a fine grid
β,γ,α = p
λ_true = β .* Igrid.^α
λ = [foi([I], res_psm.minimizer)[1] for I in Igrid]
scatter(Igrid,
        λ,
        xlabel="Proportion of population infected, I",
        ylab="Force of infection, λ",
        label="Model prediction")
Plots.vline!([Imax], color=:orange, label="Upper bound of training data")
plot!(Igrid, λ_true, color=:red, label="True function")


text:
#=
works for the following formats:
- `@forward(f)(args...)`
- `@forward f(args...) = ...`
- `@forward f = (args...) -> ...`
=#
function annotate_func_expr(typesym, _module_, expr)
    if isa(expr, Expr) && (expr.head == :(=) || expr.head == :function)
        lhs = expr.args[1]
        if isa(lhs, Expr) && (lhs.head == :call || lhs.head == :where) # named function definition site
            given_name = lhs.head == :where ? lhs.args[1].args[1] : lhs.args[1]
            @assert isa(given_name, Symbol) "potentially malformed function signature for $typesym"
            hidden_name = Symbol("#hidden_$(given_name)")
            if lhs.head == :where
                lhs.args[1].args[1] = hidden_name
            else
                lhs.args[1] = hidden_name
            end
            return quote
                $expr
                if !(isdefined($(_module_), $(Expr(:quote, given_name))))
                    const $(given_name) = ReverseDiff.$(typesym)($(hidden_name))
                end
            end
        elseif isa(lhs, Symbol) # variable assignment site
            expr.args[2] = :(ReverseDiff.$(typesym)($(expr.args[2])))
            return expr
        else
            error("failed to apply $typesym to expression $expr")
        end
    else # call site
        return :(ReverseDiff.$(typesym)($expr))
    end
end

#############################
# Forward-Mode Optimization #
#############################

struct ForwardOptimize{F}
    f::F
end

@inline ForwardOptimize(f::ForwardOptimize) = f

"""
    ReverseDiff.@forward(f)(args::Real...)
    ReverseDiff.@forward f(args::Real...) = ...
    ReverseDiff.@forward f = (args::Real...) -> ...

Declare that the given function should be differentiated using forward mode automatic
differentiation. Note that the macro can be used at either the definition site or at the
call site of `f`. Currently, only `length(args) <= 2` is supported. **Note that, if `f` is
defined within another function `g`, `f` should not close over any differentiable input of
`g`.** By using this macro, you are providing a guarantee that this property holds true.

This macro can be very beneficial for performance when intermediate functions in your
computation are low dimensional scalar functions, because it minimizes the number of
instructions that must be recorded to the tape. For example, take the function
`sigmoid(n) = 1. / (1. + exp(-n))`. Normally, using ReverseDiff to differentiate this
function would require recording 4 instructions (`-`, `exp`, `+`, and `/`). However, if we
apply the `@forward` macro, only one instruction will be recorded (`sigmoid`). The `sigmoid`
function will then be differentiated using ForwardDiff's `Dual` number type.

This is also beneficial for higher-order elementwise function application. ReverseDiff
overloads `map`/`broadcast` to dispatch on `@forward`-applied functions. For example,
`map(@forward(f), x)` will usually be more performant than `map(f, x)`.

ReverseDiff overloads many Base scalar functions to behave as `@forward` functions by
default. A full list is given by `DiffRules.diffrules()`.
"""
macro forward(ex)
    return esc(annotate_func_expr(:ForwardOptimize, __module__, ex))
end

# fallback #
#----------#

@inline (self::ForwardOptimize{F})(args...) where {F} = self.f(args...)

# unary #
#-------#

@inline function (self::ForwardOptimize{F})(t::TrackedReal{V,D}) where {F,V,D}
    T = promote_type(V, D)
    result = DiffResult(zero(T), zero(T))
    result = ForwardDiff.derivative!(result, self.f, value(t))
    tp = tape(t)
    out = track(DiffResults.value(result), D, tp)
    cache = RefValue(DiffResults.derivative(result))
    record!(tp, ScalarInstruction, self.f, t, out, cache)
    return out
end

# binary #
#--------#

@inline function (self::ForwardOptimize{F})(a::TrackedReal{V1,D}, b::TrackedReal{V2,D}) where {F,V1,V2,D}
    T = promote_type(V1, V2, D)
    result = DiffResults.GradientResult(SVector(zero(T), zero(T)))
    result = ForwardDiff.gradient!(result, x -> self.f(x[1], x[2]), SVector(value(a), value(b)))
    tp = tape(a, b)
    out = track(DiffResults.value(result), D, tp)
    cache = RefValue(DiffResults.gradient(result))
    record!(tp, ScalarInstruction, self.f, (a, b), out, cache)
    return out
end

@inline function (self::ForwardOptimize{F})(x::Real, t::TrackedReal{V,D}) where {F,V,D}
    T = promote_type(typeof(x), V, D)
    result = DiffResult(zero(T), zero(T))
    result = ForwardDiff.derivative!(result, vt -> self.f(x, vt), value(t))
    tp = tape(t)
    out = track(DiffResults.value(result), D, tp)
    dt = DiffResults.derivative(result)
    cache = RefValue(SVector(dt, dt))
    record!(tp, ScalarInstruction, self.f, (x, t), out, cache)
    return out
end

@inline function (self::ForwardOptimize{F})(t::TrackedReal{V,D}, x::Real) where {F,V,D}
    T = promote_type(typeof(x), V, D)
    result = DiffResult(zero(T), zero(T))
    result = ForwardDiff.derivative!(result, vt -> self.f(vt, x), value(t))
    tp = tape(t)
    out = track(DiffResults.value(result), D, tp)
    dt = DiffResults.derivative(result)
    cache = RefValue(SVector(dt, dt))
    record!(tp, ScalarInstruction, self.f, (t, x), out, cache)
    return out
end

@inline (self::ForwardOptimize{F})(x::Dual, t::TrackedReal) where {F} = invoke(self.f, Tuple{Dual,Real}, x, t)
@inline (self::ForwardOptimize{F})(t::TrackedReal, x::Dual) where {F} = invoke(self.f, Tuple{Real,Dual}, t, x)

#################################
# Skip Instruction Optimization #
#################################

struct SkipOptimize{F}
    f::F
end

@inline SkipOptimize(f::SkipOptimize) = f

"""
    ReverseDiff.@skip(f)(args::Real...)
    ReverseDiff.@skip f(args::Real...) = ...
    ReverseDiff.@skip f = (args::Real...) -> ...

Declare that the given function should be skipped during the instruction-recording phase of
differentiation. Note that the macro can be used at either the definition site or at the
call site of `f`. **Note that, if `f` is defined within another function `g`, `f` should not
close over any differentiable input of `g`.** By using this macro, you are providing a
guarantee that this property holds true.

ReverseDiff overloads many Base scalar functions to behave as `@skip` functions by default.
A full list is given by `ReverseDiff.SKIPPED_UNARY_SCALAR_FUNCS` and
`ReverseDiff.SKIPPED_BINARY_SCALAR_FUNCS`.
"""
macro skip(ex)
    return esc(annotate_func_expr(:SkipOptimize, __module__, ex))
end

@inline (self::SkipOptimize{F})(args...) where {F} = self.f(map(value, args)...)
@inline (self::SkipOptimize{F})(a) where {F} = self.f(value(a))
@inline (self::SkipOptimize{F})(a, b) where {F} = self.f(value(a), value(b))
@inline (self::SkipOptimize{F})(a, b, c) where {F} = self.f(value(a), value(b), value(c))

"""
    f(x) = dot(x, x)
    f(x::ReverseDiff.TrackedVector) = ReverseDiff.track(f, x)
    ReverseDiff.@grad function f(x)
        xv = ReverseDiff.value(x)
        return dot(xv, xv), Δ -> (Δ * 2 * xv,)
    end
The `@grad` macro provides a way for the users to define custom adjoints for single-output functions wrt to their input numbers or arrays.
"""
macro grad(expr)
    d = MacroTools.splitdef(expr)
    f = d[:name]
    closure = gensym(f)
    d[:name] = closure
    closure_ex = MacroTools.combinedef(d)

    @gensym tp output_value output back args kwargs
    args_ex = getargs_expr(d[:args])
    kwargs_ex = getkwargs_expr(d[:kwargs])
    return quote
        function $ReverseDiff.track(::typeof($f), $(d[:args]...); $(d[:kwargs]...)) where {$(d[:whereparams]...),}
            $closure_ex
            $args = $args_ex
            $kwargs = $kwargs_ex
            $tp = $ReverseDiff.tape($args...)
            $output_value, $back = $closure($args...; $kwargs...)
            $output = $ReverseDiff.track($output_value, $tp)
            $ReverseDiff.record!(
                $tp,
                $ReverseDiff.SpecialInstruction,
                $f,
                $args,
                $output,
                ($back, $closure, $kwargs),
            )
            return $output
        end

        if !hasmethod(
            $ReverseDiff.special_reverse_exec!,
            Tuple{$ReverseDiff.SpecialInstruction{typeof($f)}},
        )
            @noinline function $ReverseDiff.special_reverse_exec!(instruction::$ReverseDiff.SpecialInstruction{typeof($f)})
                output = instruction.output
                input = instruction.input
                back = instruction.cache[1]
                input_derivs = back($ReverseDiff.deriv(output))
                @assert input_derivs isa Tuple
                $ReverseDiff._add_to_deriv!.(input, input_derivs)
                $ReverseDiff.unseed!(output)
                return nothing
            end
        end

        if !hasmethod(
            $ReverseDiff.special_forward_exec!,
            Tuple{$ReverseDiff.SpecialInstruction{typeof($f)}},
        )
            @noinline function $ReverseDiff.special_forward_exec!(instruction::$ReverseDiff.SpecialInstruction{typeof($f)})
                output, input = instruction.output, instruction.input
                $ReverseDiff.pull_value!.(input)
                pullback = instruction.cache[2]
                kwargs = instruction.cache[3]
                out_value = pullback(input...; kwargs...)[1]
                $ReverseDiff.value!(output, out_value)
                return nothing
            end
        end
    end |> esc
end
_add_to_deriv!(d1, d2) = nothing
function _add_to_deriv!(d1::Union{TrackedReal, TrackedArray}, d2)
    increment_deriv!(d1, d2)
end
function getargs_expr(args_with_types)
    expr = Expr(:tuple)
    for at in args_with_types
        x, tosplat = remove_tp(at)
        if tosplat
            push!(expr.args, :($x...))
        else
            push!(expr.args, x)
        end
    end
    return expr
end
function getkwargs_expr(kwargs_with_types)
    syms = []
    final = nothing
    for at in kwargs_with_types
        final isa Nothing || throw("Invalid kwargs.")
        x, tosplat = remove_tp(at)
        if tosplat
            final = x
        else
            push!(syms, x)
        end
    end
    expr = length(syms) == 0 ? :(NamedTuple()) : Expr(:tuple, [:($f = $f) for f in syms]...)
    final = final == nothing ? :(NamedTuple()) : final
    return :(Base.merge($expr, $final))
end
function remove_tp(t)
    if @capture(t, X_::T_...)
        return X, true
    elseif @capture(t, X_::T_)
        return X, false
    elseif @capture(t, X_::T_ = V_)
        return X, false
    elseif @capture(t, ::typeof(T_)...)
        return T, true
    elseif @capture(t, ::typeof(T_))
        return T, false
    elseif @capture(t, X_...)
        return X, true
    elseif @capture(t, X_ = V_)
        return X, false
    else
        return t, false
    end
end


text:
using Flux

function build_model(name="MNIST"; kw...)
    return @eval $(Symbol(:build_, name))(;$kw...)
end

# Build model
function build_MNIST(; imgsize = (28,28, 1), nclasses = 10)
    cnn_output_size = Int.(floor.([imgsize[1]/8,imgsize[2]/8,32]))

    return Chain(
    # First convolution, operating upon a 28x28 image
    Conv((3, 3), imgsize[3]=>16, pad=(1,1), relu),
    MaxPool((2,2)),

    # Second convolution, operating upon a 14x14 image
    Conv((3, 3), 16=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    # Third convolution, operating upon a 7x7 image
    Conv((3, 3), 32=>32, pad=(1,1), relu),
    MaxPool((2,2)),

    # Reshape 3d tensor into a 2d one using `Flux.flatten`, at this point it should be (3, 3, 32, N)
    flatten,
    Dense(prod(cnn_output_size), 10))
end

# VGG16
function build_CIFAR10()
    Chain(
        Conv((3, 3), 3 => 64, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(64),
        Conv((3, 3), 64 => 64, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(64),
        MaxPool((2,2)),
        Conv((3, 3), 64 => 128, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(128),
        Conv((3, 3), 128 => 128, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(128),
        MaxPool((2,2)),
        Conv((3, 3), 128 => 256, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(256),
        Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(256),
        Conv((3, 3), 256 => 256, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(256),
        MaxPool((2,2)),
        Conv((3, 3), 256 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        MaxPool((2,2)),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        Conv((3, 3), 512 => 512, relu, pad=(1, 1), stride=(1, 1)),
        BatchNorm(512),
        MaxPool((2,2)),
        flatten,
        Dense(512, 4096, relu),
        Dropout(0.5),
        Dense(4096, 4096, relu),
        Dropout(0.5),
        Dense(4096, 10)
    )
end

max_stars_repo_name: sheevy/XConv

max_issues_repo_name: sheevy/XConv

max_forks_repo_name: sheevy/XConv