text:
################################################################################
# For tractability, a target parameter must have the form:
#   \sum_{z, \ell, d}
#       int_constant(ℓ, d, z)
#       \times
#       \int_{int_limits(z)[1]}^{int_limits(z)[2]}
#       m(u_{ℓ} | Z_{-ℓ}) du
#
#   So defining a target parameter requires a pair of functions:
#   - int_limits, which takes z and returns a lower and upper limit
#   of integration
#   - int_constant, which takes (ℓ, d, z) and returns a scalar.
#
#   Note that expectation over Z is incorporated into int_constant(ℓ, d, z)
################################################################################
struct TargetParameter
    name::String
    int_limits::Function
    int_constant::Function

    function TargetParameter(name, int_limits, int_constant)
        new(name, int_limits, int_constant)
    end
end
export TargetParameter

function eval_tp(tp::TargetParameter, mtrs::Vector{Tuple{MTR, MTR}}, dgp)
    Γ⭑ = compute_Γ⭑(tp, [(mtr[1].basis, mtr[2].basis) for mtr in mtrs], dgp)
    total = 0
    for ℓ in 1:length(mtrs), d in 0:1
        total += sum(Γ⭑[ℓ, d + 1] .* mtrs[ℓ][d + 1].θ)
    end
    return total
end
export eval_tp

function compute_Γ⭑(tp::TargetParameter,
                    bases::Vector{Tuple{MTRBasis, MTRBasis}}, dgp::DGP)
    [compute_Γ⭑(tp, bases[ℓ][d + 1], d, ℓ, dgp)
     for ℓ in 1:length(bases), d in 0:1]
end

function compute_Γ⭑(tp::TargetParameter, basis::MTRBasis,
                    d::Integer, ℓ::Integer, dgp::DGP)
    Γ⭑ = zeros(length(basis.a), length(basis.b))
    for z in eachrow(dgp.suppZ)
        il = tp.int_limits(z)
        Γ⭑ += [(aj(z) * ibk(il[1], il[2]) * tp.int_constant(ℓ, d, z))
               for aj in basis.a, ibk in basis.ib]
    end
    return Γ⭑
end

function ey1(dgp::DGP; ℓ = 1)
    TargetParameter(
        "EY1",
        z -> (0,1),
        (l,d,z) -> d * (l == ℓ) * find_density(z, dgp)
    )
end
export ey1

function att(dgp::DGP; ℓ = 1)
    prd1 = dot(dgp.pscore, dgp.densZ)
    TargetParameter(
        "ATT",
        z -> (0, find_pscore(z, dgp)),
        (l,d,z) -> ((l == ℓ) * (2*d - 1)/prd1) * find_density(z, dgp)
    )
end
export att

function prte_plusδpercent(dgp::DGP, δ; ℓ = 1)
    @assert ℓ == 1 # haven't coded the other case
    limits = z -> (          find_pscore([0, z[2]], dgp),
                   (1 + δ) * find_pscore([1, z[2]], dgp))
    TargetParameter(
        "PRTE_plus$(δ)_ℓ=$(ℓ)",
        limits,
        (l,d,z) -> ((l == ℓ) * (2*d - 1) *
                    find_density(z, dgp) / (limits(z)[2] - limits(z)[1]))
    )
end
export prte_plusδpercent

################################################################################
# Special type of PRTE
#   dgpnew and dgpold should have the same "support" points for Z
#   note that you can set their probability to 0 in the old distribution
#   and that's totally fine in the program
#
# All that's shifting here is the marginal distribution of Z
################################################################################
function prte_newz(dgpold::DGP, dgpnew::DGP; ℓ = 1)
    @assert dgpold.suppZ == dgpnew.suppZ
    @assert dgpold.pscore == dgpnew.pscore
    TargetParameter(
        "PRTE New Z",
        z -> (0, find_pscore(z, dgpold)),
        (l, d, z) -> ((l == ℓ) * (2*d - 1) *
                     (find_density(z, dgpnew) - find_density(z, dgpold)))
    )
end
export prte_newz


text:
#=
Copyright (c) 2015, Intel Corporation

Redistribution and use in source and binary forms, with or without
modification, are permitted provided that the following conditions are met:

    * Redistributions of source code must retain the above copyright notice,
      this list of conditions and the following disclaimer.
    * Redistributions in binary form must reproduce the above copyright
      notice, this list of conditions and the following disclaimer in the
      documentation and/or other materials provided with the distribution.
    * Neither the name of Intel Corporation nor the names of its contributors
      may be used to endorse or promote products derived from this software
      without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE
FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL
DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR
SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER
CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
=#

using FactCheck
using Latte

net = Net(8)
data,  data_value   = MemoryDataLayer(net, :data, (20, ))
label, label_value = MemoryDataLayer(net, :label, (1,))
data_value[:]  = rand(Float32, size(data_value)...) * 256
label_value[:] = map(floor, rand(Float32, size(label_value)...) * 10)
fc1 = InnerProductLayer(:fc1, net, data, 20)
relu1    = ReLULayer(:relu1, net, fc1)
ratio = 0.5f0
drop = DropoutLayer(:drop, net, fc1, ratio)

init(net)
params = SolverParameters(
    lr_policy    = LRPolicy.Decay(.01f0, 5.0f-7),
    mom_policy   = MomPolicy.Fixed(0.9),
    max_epoch    = 300,
    regu_coef    = .0005)
sgd = SGD(params)

input = get_buffer(net, :fc1value)

facts("Testing Dropout Layer") do
    context("Forward") do
        forward(net; solver=sgd)
        randvals = get_buffer(net, :droprandval)
        weights = get_buffer(net, :fc1weights)
        expected = weights' * data_value
        expected = (expected .> 0.0f0) .* expected
        expected = (randvals .> ratio) .* expected .* (1.0 / ratio)
        value = get_buffer(net, :dropvalue)
        @fact value --> roughly(expected)
    end
    context("Backward") do
        top_diff = get_buffer(net, :drop∇)
        rand!(top_diff)
        randvals = get_buffer(net, :droprandval)
        expected = (top_diff .* (randvals .> ratio)) .* (1.0 / ratio)

        value = get_buffer(net, :fc1value)
        expected = (value .> 0.0f0) .* expected
        backward(net)

        @fact expected --> roughly(get_buffer(net, :fc1∇))
    end
end

FactCheck.exitstatus()

ext: jl


text:
import Base: convert

export
  MultiplyDistributionsRequest

struct MultiplyDistributionsRequest
    weights::Vector{Packed_BallTreeDensity}
end


text:
"""
Algorithm for solving Physics Informed Neural Networks problem.

Arguments:
* `dx` is discretization of grid
* `chain` is a Flux.jl chain with d dimensional input and 1 dimensional output,
* `init_params` is the initial parameter of the neural network,
* `phi` is trial solution,
* `autodiff` is a boolean variable that determines whether to use automatic, differentiation(not supported while) or numerical,
* `derivative` is method that calculate derivative,
* `strategy` is determines which training strategy will be used.
"""

struct PhysicsInformedNN{D,C,P,PH,DER,T,K}
  dx::D
  chain::C
  initθ::P
  phi::PH
  autodiff::Bool
  derivative::DER
  strategy::T
  kwargs::K
end

function PhysicsInformedNN(dx,
                           chain,
                           init_params = nothing;
                           _phi = nothing,
                           autodiff=false,
                           _derivative = nothing,
                           strategy = GridTraining(),
                           kwargs...)
    if init_params === nothing
        if chain isa FastChain
            initθ = DiffEqFlux.initial_params(chain)
        else
            initθ,re  = Flux.destructure(chain)
        end
    else
        initθ = init_params
    end
    isuinplace = chain.layers[end].out == 1
    dim = chain.layers[1].in

    if _phi == nothing
        phi = get_phi(chain,isuinplace)
    else
        phi = _phi
    end

    if _derivative == nothing
        derivative = get_derivative(dim,phi,autodiff,isuinplace)
    else
        derivative = _derivative
    end

    PhysicsInformedNN(dx,chain,initθ,phi,autodiff,derivative, strategy, kwargs)
end

abstract type TrainingStrategies  end

struct GridTraining <: TrainingStrategies end

struct StochasticTraining <:TrainingStrategies
    include_frac::Float64
end
function StochasticTraining(;include_frac=0.75)
    StochasticTraining(include_frac)
end

"""
Create dictionary: variable => unique number for variable

# Example 1

Dict{Symbol,Int64} with 3 entries:
  :y => 2
  :t => 3
  :x => 1

# Example 2

 Dict{Symbol,Int64} with 2 entries:
  :u1 => 1
  :u2 => 2
"""
get_dict_vars(vars) = Dict( [Symbol(v) .=> i for (i,v) in enumerate(vars)])

# Calculate an order of derivative
function count_order(_args)
    n = 0
    while (_args[1] == :derivative)
        n += 1
        _args = _args[2].args
    end
    return n
end

function get_depvars(_args)
    while (_args[1] == :derivative)
        _args = _args[2].args
    end
    return _args[1]
end
function get_indpvars(_args)
    while (_args[1] == :derivative)
        _args = _args[2].args
    end
    return _args[2:end-1]
end

# Wrapper for _simplified_derivative
function simplified_derivative(ex,indvars,depvars,dict_indvars,dict_depvars)
    if ex isa Expr
        ex.args = _simplified_derivative(ex.args,indvars,depvars,dict_indvars,dict_depvars)
    end
    return ex
end

"""
Simplify the derivative expression

# Examples

1. First derivative of function 'u(x,y)' of variables: 'x, y'
1. First derivative of function 'u(x,y)' with respect to x

Take expressions in the form: `derivative(u(x,y,θ), x)` to `derivative(unn, x, y, xnn, order, θ)`,

where
 unn - unique number for the function
 x,y - coordinates of point
 xnn - unique number for the variable
 order - order of derivative
 θ - parameters of neural network
"""
function _simplified_derivative(_args,indvars,depvars,dict_indvars,dict_depvars)
    for (i,e) in enumerate(_args)
        if !(e isa Expr)
            if e == :derivative && _args[end] != :θ
                order = count_order(_args)
                depvars = get_depvars(_args)
                indvars = get_indpvars(_args)
                depvars_num = dict_depvars[depvars]
                _args = [:derivative,depvars_num, indvars...,dict_indvars[_args[3]], order, :θ]
                break
            end
        else
            _args[i].args = _simplified_derivative(_args[i].args,indvars,depvars,dict_indvars,dict_depvars)
        end
    end
    return _args
end

"""
Parse ModelingToolkit equation form to the inner representation.

Example:

1)  1d ODE: Dt(u(t,θ)) ~ t +1

    Take expressions in the form: 'Equation(derivative(u(t, θ), t), t + 1)' to 'derivative(u, t, 1, 1, θ) - (t + 1)'

2)  2d PDE: Dxx(u(x,y,θ)) + Dyy(u(x,y,θ)) ~ -sin(pi*x)*sin(pi*y)

    Take expressions in the form:
     Equation(derivative(derivative(u(x, y, θ), x), x) + derivative(derivative(u(x, y, θ), y), y), -(sin(πx)) * sin(πy))
    to
     (derivative(1, x, y, 1, 2, θ) + derivative(1, x, y, 2, 2, θ)) - -(sin(πx)) * sin(πy)

3)  System of PDE: [Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0,
                    Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0]

    Take expressions in the form:
    2-element Array{Equation,1}:
        Equation(derivative(u1(x, y, θ), x) + 4 * derivative(u2(x, y, θ), y), ModelingToolkit.Constant(0))
        Equation(derivative(u2(x, y, θ), x) + 9 * derivative(u1(x, y, θ), y), ModelingToolkit.Constant(0))
    to
      [(derivative(1, x, y, 1, 1, θ) + 4 * derivative(2, x, y, 2, 1, θ)) - 0,
       (derivative(2, x, y, 1, 1, θ) + 9 * derivative(1, x, y, 2, 1, θ)) - 0]
"""
function parse_equation(eq,indvars,depvars,dict_indvars,dict_depvars)
    left_expr = simplified_derivative((ModelingToolkit.simplified_expr(eq.lhs)),
                                      indvars,depvars,dict_indvars,dict_depvars)
    right_expr = simplified_derivative((ModelingToolkit.simplified_expr(eq.rhs)),
                                       indvars,depvars,dict_indvars,dict_depvars)
    loss_func = :($left_expr - $right_expr)
end

"""
Build loss function for pde or boundary condition

# Examples: System of pde:

Take expressions in the form:

[Dx(u1(x,y,θ)) + 4*Dy(u2(x,y,θ)) ~ 0,
 Dx(u2(x,y,θ)) + 9*Dy(u1(x,y,θ)) ~ 0]

to

:((phi, vec, θ, derivative)->begin
          #= none:33 =#
          #= none:33 =#
          begin
              begin
                  u1 = ((x, y, θ)->begin
                              #= none:18 =#
                              (phi([x, y], θ))[1]
                          end)
                  u2 = ((x, y, θ)->begin
                              #= none:18 =#
                              (phi([x, y], θ))[2]
                          end)
              end
              let (x, y) = (vec[1], vec[2])
                  [(derivative(1, x, y, 1, 1, θ) + 4 * derivative(2, x, y, 2, 1, θ)) - 0,
                   (derivative(2, x, y, 1, 1, θ) + 9 * derivative(1, x, y, 2, 1, θ)) - 0]
              end
          end
      end)
"""
function build_loss_function(eqs,_indvars,_depvars)
    # dictionaries: variable -> unique number
    indvars = Expr.(_indvars)
    depvars = [d.name for d in _depvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return build_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars)
end

function build_loss_function(eqs,indvars,depvars,dict_indvars,dict_depvars)
    if !(eqs isa Array)
        eqs = [eqs]
    end
    loss_functions= []
    for eq in eqs
        loss_function = parse_equation(eq,depvars,indvars,dict_indvars,dict_depvars)
        push!(loss_functions,loss_function)
    end

    vars = :(phi, cord, θ, derivative)
    ex = Expr(:block)
    us = []
    for v in depvars
        var_num = dict_depvars[v]
        push!(us,:($v = ($(indvars...), θ) -> phi([$(indvars...)],θ)[$var_num]))
    end

    u_ex = ModelingToolkit.build_expr(:block, us)
    push!(ex.args,  u_ex)

    indvars_ex = [:($:cord[$i]) for (i, u) ∈ enumerate(indvars)]
    arg_pairs_indvars = indvars,indvars_ex

    left_arg_pairs, right_arg_pairs = arg_pairs_indvars
    vars_eq = Expr(:(=), ModelingToolkit.build_expr(:tuple, left_arg_pairs),
                            ModelingToolkit.build_expr(:tuple, right_arg_pairs))
    let_ex = Expr(:let, vars_eq, ModelingToolkit.build_expr(:vect, loss_functions))
    push!(ex.args,  let_ex)

    return :(($vars) -> begin $ex end)
end

# get agruments from bondary condition functions
function get_bc_argument(bcs,indvars,depvars,dict_indvars,dict_depvars)
    bc_args = []
    for _bc in bcs
        if !(_bc isa Array)
            _bc = [_bc]
        end
        left_expr = simplified_derivative(ModelingToolkit.simplified_expr(_bc[1].lhs),
                                          indvars,depvars,dict_indvars,dict_depvars)
        bc_arg = if (left_expr.args[1] == :derivative)
            left_expr.args[3:end-3]
        else
            left_expr.args[2:end-1]
        end

        push!(bc_args, bc_arg)
    end
    return bc_args
end

function generate_training_sets(domains,dx,bcs,_indvars,_depvars)
    indvars = Expr.(_indvars)
    depvars = [d.name for d in _depvars]
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)
    return generate_training_sets(domains,dx,bcs,indvars,depvars,dict_indvars,dict_depvars)
end
# Generate training set in the domain and on the boundary
function generate_training_sets(domains,dx,bcs,indvars,depvars,dict_indvars,dict_depvars)
    if dx isa Array
        dxs = dx
    else
        dxs = fill(dx,length(domains))
    end

    bound_args = get_bc_argument(bcs,indvars,depvars,dict_indvars,dict_depvars)

    spans = [d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)]
    dict_var_span = Dict([Symbol(d.variables) => d.domain.lower:dx:d.domain.upper for (d,dx) in zip(domains,dxs)])

    train_set = []
    for points in Iterators.product(spans...)
        push!(train_set, [points...])
    end

    train_bound_set = []
    for bt in bound_args
        _set = []
        span = [get(dict_var_span, b, b) for b in bt]
        for points in Iterators.product(span...)
            push!(_set, [points...])
        end
        push!(train_bound_set, _set)
    end

    flat_train_bound_set = collect(Iterators.flatten(train_bound_set))
    train_domain_set =  setdiff(train_set, flat_train_bound_set)

    [train_domain_set,train_bound_set,train_set]
end

function get_phi(chain,isuinplace)
    # The phi trial solution
    if chain isa FastChain
        if isuinplace
            phi = (x,θ) -> first(chain(adapt(typeof(θ),x),θ))
        else
            phi = (x,θ) -> chain(adapt(typeof(θ),x),θ)
        end
    else
        _,re  = Flux.destructure(chain)
        if isuinplace
            phi = (x,θ) -> first(re(θ)(adapt(typeof(θ),x)))
        else
            phi = (x,θ) -> re(θ)(adapt(typeof(θ),x))
        end
    end
    phi
end

function get_derivative(dim,phi,autodiff,isuinplace)
    if autodiff # automatic differentiation (not implemented yet)
        # derivative = (x,θ,n) -> ForwardDiff.gradient(x->phi(x,θ),x)[n]
    else # numerical differentiation
        epsilon = cbrt(eps(Float32))
        function get_ε(dim, der_num)
            ε = zeros(dim)
            ε[der_num] = epsilon
            ε
        end

        εs = [get_ε(dim,d) for d in 1:dim]

        derivative = (_x...) ->
        begin
            var_num = _x[1]
            x = collect(_x[2:end-3])
            der_num =_x[end-2]
            order = _x[end-1]
            θ = _x[end]
            ε = εs[der_num]
            return _derivative(x,θ,order,ε,der_num,var_num)
        end
        _derivative = (x,θ,order,ε,der_num,var_num) ->
        begin
            if order == 1
                #Five-point stencil
                # return (-phi(x+2ε,θ) + 8phi(x+ε,θ) - 8phi(x-ε,θ) + phi(x-2ε,θ))/(12*epsilon)
                if isuinplace
                    return (phi(x+ε,θ) - phi(x-ε,θ))/(2*epsilon)
                else
                    return (phi(x+ε,θ)[var_num] - phi(x-ε,θ)[var_num])/(2*epsilon)
                end
            else
                return (_derivative(x+ε,θ,order-1,ε,der_num,var_num)
                      - _derivative(x-ε,θ,order-1,ε,der_num,var_num))/(2*epsilon)
            end
        end
    end
    derivative
end

function get_loss_function(loss_function, train_set, phi, derivative, strategy)

    # norm coefficient for loss function
    τ = sum(length(set) for set in train_set)

    function inner_loss(loss_function,phi,x,θ,derivative)
        sum(loss_function(phi, x, θ, derivative))
    end

    if !(loss_function isa Array)
        loss_function = [loss_function]
        train_set = [train_set]
    end

    if strategy isa StochasticTraining
        include_frac = strategy.include_frac
        count_elements = []
        sets_size = []
        for j in 1:length(train_set)
            size_set = size(train_set[j])[1]
            count_element = convert(Int64,round(include_frac*size_set, digits=0))
            if count_element <= 2
                count_element = size_set
            end
            push!(sets_size,size_set)
            push!(count_elements,count_element)
        end
        loss = (θ) -> begin
            total = 0.
            for (j,l) in enumerate(loss_function)
                size_set = sets_size[j]
                for i in 1:count_elements[j]
                    index = convert(Int64, round(size_set*rand(1)[1] + 0.5, digits=0))
                    total += inner_loss(l,phi,train_set[j][index],θ,derivative)^2
                end
            end
            return (1.0f0/τ) * total
        end

    elseif strategy isa GridTraining
        loss = (θ) -> begin
            return (1.0f0/τ) *sum(sum(abs2,inner_loss(l,phi,x,θ,derivative)
                    for x in set) for (l,set) in zip(loss_function,train_set))
        end

    end
    return loss
end

# Convert a PDE problem into OptimizationProblem
function DiffEqBase.discretize(pde_system::PDESystem, discretization::PhysicsInformedNN)
    eqs = pde_system.eq
    bcs = pde_system.bcs

    domains = pde_system.domain
    # dimensionality of equation
    dim = length(domains)
    dim > 3 && error("While only dimensionality no more than 3")

    depvars = [d.name for d in pde_system.depvars]
    indvars = Expr.(pde_system.indvars)
    dict_indvars = get_dict_vars(indvars)
    dict_depvars = get_dict_vars(depvars)

    dx = discretization.dx
    chain = discretization.chain
    initθ = discretization.initθ
    phi = discretization.phi
    autodiff = discretization.autodiff
    derivative = discretization.derivative
    autodiff == true && error("Automatic differentiation is not support yet")
    strategy = discretization.strategy

    train_sets = generate_training_sets(domains,dx,bcs,
                                        indvars,depvars,
                                        dict_indvars,dict_depvars)

    # the points in the domain and on the boundary
    train_domain_set,train_bound_set,train_set = train_sets

    expr_pde_loss_function = build_loss_function(eqs,indvars,depvars,
                                                 dict_indvars,dict_depvars)
    expr_bc_loss_functions = [build_loss_function(bc,indvars,depvars,
                                                  dict_indvars,dict_depvars) for bc in bcs]

    pde_loss_function = get_loss_function(eval(expr_pde_loss_function),
                                          train_domain_set,
                                          phi,
                                          derivative,
                                          strategy)
    bc_loss_function = get_loss_function(eval.(expr_bc_loss_functions),
                                         train_bound_set,
                                         phi,
                                         derivative,
                                         strategy)

    function loss_function(θ,p)
        return pde_loss_function(θ) + bc_loss_function(θ)
    end

    f = OptimizationFunction(loss_function, initθ, GalacticOptim.AutoZygote())
    GalacticOptim.OptimizationProblem(f, initθ)
end


text:
radiality_centrality(g::AbstractGraph; parallel=:distributed) =
parallel == :distributed ? distr_radiality_centrality(g) : threaded_radiality_centrality(g)

function distr_radiality_centrality(g::AbstractGraph)::Vector{Float64}
    n_v = nv(g)
    vs = vertices(g)
    n = ne(g)
    meandists = SharedVector{Float64}(Int(n_v))
    maxdists = SharedVector{Float64}(Int(n_v))

    @sync @distributed for i = 1:n_v
        d = LightGraphs.dijkstra_shortest_paths(g, vs[i])
        maxdists[i] = maximum(d.dists)
        meandists[i] = sum(d.dists) / (n_v - 1)
        nothing
    end
    dmtr = maximum(maxdists)
    radialities = collect(meandists)
    return ((dmtr + 1) .- radialities) ./ dmtr
end

function threaded_radiality_centrality(g::AbstractGraph)::Vector{Float64}
    n_v = nv(g)
    vs = vertices(g)
    n = ne(g)
    meandists = Vector{Float64}(undef, n_v)
    maxdists = Vector{Float64}(undef, n_v)

    Base.Threads.@threads for i in vertices(g)
        d = LightGraphs.dijkstra_shortest_paths(g, vs[i])
        maxdists[i] = maximum(d.dists)
        meandists[i] = sum(d.dists) / (n_v - 1)
    end
    dmtr = maximum(maxdists)
    radialities = collect(meandists)
    return ((dmtr + 1) .- radialities) ./ dmtr
end


text:
import LibGit2: checkout!, clone, peel, tag_list, CloneOptions, GitCommit, GitHash, GitObject, GitTree
import Pkg: TOML

Packages = Dict()
PkgsDir = "packages/"
ActivePkgs = String[]

for url in readlines("tools/packages.txt")
    startswith(url, "#") && continue

    dir = mktempdir()
    @info "Cloning $(url) to $(dir)..."
    repo = clone(url, dir)

    Deps = Dict()
    Versions = Dict()

    project = TOML.parsefile(joinpath(dir, "Project.toml"))

    Package = Dict(
        "name" => project["name"],
        "uuid" => project["uuid"],
        "repo" => url
    )

    pkgpath = joinpath(PkgsDir, Package["name"])
    push!(ActivePkgs, Package["name"])

    Packages[Package["uuid"]] = Dict(
        "name" => Package["name"],
        "path" => pkgpath
    )

    for tag in tag_list(repo)
        obj = GitObject(repo, tag)
        obj_hash = string(GitHash(obj))

        tree = peel(obj)
        tree_hash = string(GitHash(tree))

        checkout!(repo, obj_hash)

        project = TOML.parsefile(joinpath(dir, "Project.toml"))
        version = project["version"]
	if match(r"^[\d\.]+$", version) === nothing
	    @info "Skipping (invalid) version $(version)"
            continue
        end

        @info "$(project["name"]) v$(version) => $(tree_hash)"

        Deps[version] = project["deps"]

        Versions[version] = Dict(
            "git-tree-sha1" => string(GitHash(tree_hash))
        )
    end

    rm(dir, recursive = true)

    mkpath(pkgpath)

    TOML.print(open(joinpath(pkgpath, "Deps.toml"), "w"), Deps, sorted = true)
    TOML.print(open(joinpath(pkgpath, "Package.toml"), "w"), Package, sorted = true)
    TOML.print(open(joinpath(pkgpath, "Versions.toml"), "w"), Versions, sorted = true)
end

# Cleanup packages not present in packages.txt
for pkg in readdir(PkgsDir)
    if !(pkg in ActivePkgs)
        @info "Removing $pkg"
        rm(joinpath(PkgsDir, pkg), recursive = true)
    end
end

registry = TOML.parsefile("Registry.toml")
registry["packages"] = Packages
TOML.print(open("Registry.toml", "w"), registry)

max_stars_repo_name: maxmouchet/JuliaRegistry

max_issues_repo_name: maxmouchet/JuliaRegistry

max_forks_repo_name: maxmouchet/JuliaRegistry


text:
export Term,     @term
export Symbolic, @syms
export Variable, @vars

using MacroTools

macro term(ex)
     :(convert(Term, $(esc(_term(ex)))))
end
function _term(ex::Expr)
    ex.head === :$ && return ex.args[1]
    ex.head === :. && return ex
    Expr(:call, Expr, Meta.quot(ex.head), _term.(ex.args)...)
end

_term(x) = x
_term(x::Symbol) = x === :im ? :literal_im : x

struct Term
    ex
end
Base.convert(::Type{Term}, ex::Term) = ex
Base.convert(::Type{Term}, ex) = Term(traverse(ex))
traverse(t::Term) = get(t)
traverse(ex::Expr) = Expr(ex.head, traverse.(ex.args)...)
traverse(x) = x

Base.:(==)(a::Term, b::Term) = a.ex == b.ex
Base.get(t::Term) = t.ex

Base.map(f, t::Term) = convert(Term, _map(f, get(t)))
_map(f, ex::Expr) = Expr(ex.head, map(f ∘ Term, ex.args)...)
_map(f, x) = x

Base.issubset(a::Term, b::Term) = !isempty(match(b, a))
Base.replace(t::Term, σ) = haskey(σ, get(t)) ? Term(σ[get(t)]) : map(x -> replace(x, σ), t)

struct Symbolic
    name::Symbol
end
(x::Symbolic)(xs...) = Expr(:call, x, xs...)
Base.show(io::IO, x::Symbolic) = print(io, x.name)
macro syms(xs::Symbol...)
    syms = (:($x = $(Symbolic(x))) for x ∈ xs)
    results = Expr(:tuple, xs...)
    esc(Expr(:block, syms..., results))
end

struct Im end
Base.show(io::IO, ::Im) = print(io, "im")
Base.convert(::Type{Term}, x::Complex) = Term(Expr(:call, +, real(x), Expr(:call, *, imag(x), Im())))

export literal_im
const literal_im = Im()

struct Variable
    name::Symbol
end
Base.show(io::IO, x::Variable) = print(io, x.name)
macro vars(xs::Symbol...)
    vars = (:($x = $(Variable(x))) for x ∈ xs)
    results = Expr(:tuple, xs...)
    esc(Expr(:block, vars..., results))
end

function _show_term(f::Function)
    # Inspired by: `show(::IO, ::Function)`
    ft = typeof(f)
    mt = ft.name.mt
    if isdefined(mt, :module) && isdefined(mt.module, mt.name) && getfield(mt.module, mt.name) === f
        Base.is_exported_from_stdlib(mt.name, mt.module) && return mt.name
        return :($(nameof(mt.module)).$(mt.name))
    else
        return f
    end
end

_show_term(::Irrational{sym}) where sym = sym
_show_term(x::Symbolic) = x.name
_show_term(x::Variable) = x.name
_show_term(x::Im) = :im
_show_term(x::Symbol) = Meta.quot(x)
_show_term(x) = x

function _show_term(x::Expr)
    x.head === :call || return x

    if x.args[1] === :sqrt
        return Expr(:call, :√, x.args[2:end]...)
    else
        return x
    end
end

function Base.show(io::IO, t::Term)
    ex = MacroTools.postwalk(_show_term, get(t))
    macro_call = Expr(:macrocall, Symbol("@term"), nothing, ex)
    repr = sprint(show, macro_call)[9:end-1]
    print(io, "@term(", repr, ")")
end


text:
# Delete a file
rm("input.txt")

# Delete a directory
rm("docs", recursive = true)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
export MapSliceOp

function mapSliceForeward(A, x::Vector, size1::Tuple, dim::Int)
  return vec( mapslices(u->A*u, reshape(x, size1), dims=dim) )
end

function mapSliceBackward(A, y::Vector, size2::Tuple, dim::Int)
  return vec( mapslices(v->adjoint(A)*v, reshape(y, size2), dims=dim) )
end

"""
    MapSliceOp(trafo, dim::Int64, size1::Tuple, size2::Tuple; T=ComplexF64)

generates an operator that applies `trafo` to each slice of dimension `dim` of an array

# Arguments
* `trafo`           - transformation to apply
* `size1`           - size of the Array to be transformed
* `size2`           - size of the resulting Array of applying trafo to all slices
* (`T=ComplexF64`)  - type of the transformation
"""
function MapSliceOp(trafo, dim::Int64, size1::Tuple, size2::Tuple; T=ComplexF64)
  return LinearOperator(prod(size2), prod(size1), false, false
            , x->mapSliceForeward(trafo, x, size1, dim)
            , nothing
            , y->mapSliceBackward(trafo, y, size2, dim) )
end


text:
using LazyArrays, Test
import LazyArrays: Add, AddArray, MulAdd, materialize!, MemoryLayout, ApplyLayout

@testset "Add/Subtract" begin
    @testset "Add" begin
        @testset "Mul-Add" begin
            A = Add(randn(5,5), randn(5,5))
            b = randn(5)
            c = similar(b)
            fill!(c,NaN)
            @test (c .= @~ A*b) ≈ A.args[1]*b + A.args[2]*b
        end

        @testset "gemv Float64" begin
            for A in (AddArray(randn(5,5), randn(5,5)),
                    AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))),
                b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))

                Ã = Array(A)
                c = similar(b)
                fill!(c,NaN)

                materialize!(MulAdd(1.0, A, b, 0.0, c))
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                c = similar(b)
                fill!(c,NaN)
                c .= @~ A*b
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                copyto!(c, @~ A*b)
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                b̃ = Array(b)
                copyto!(b̃, Mul(A,b̃))
                @test c ≈ b̃

                c .= @~ 2.0 * A * b
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 0.0, similar(c))

                c = Array(b)
                c .= @~ A*b + c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 2.0, Array(b))

                c = Array(b)
                c .= @~ 2.0 * A*b + c
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ 3.0 * A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))

                d = similar(c)
                c = Array(b)
                d .= @~ 3.0 * A*b + 2.0 * c
                @test d ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))
            end
        end

        @testset "gemm" begin
            for A in (AddArray(randn(5,5), randn(5,5)),
                    AddArray(randn(5,5), view(randn(9, 5), 1:2:9, :))),
                B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                    view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))

                Ã = Array(A)
                C = similar(B)

                C .= @~ A*B
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 0.0, similar(C))

                B .= @~ A*B
                @test C ≈ B

                C .= @~ 2.0 * A*B
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 0.0, similar(C))

                C = Array(B)
                C .= @~ A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 2.0, Array(B))

                C = Array(B)
                C .= @~ 2.0 * A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ 3.0 * A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))

                d = similar(C)
                C = Array(B)
                d .= @~ 3.0 * A*B + 2.0 * C
                @test d ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))
            end
        end
    end

    @testset "Subtract" begin
        @testset "Mul-Subtract" begin
            A = applied(-, randn(5,5), randn(5,5))
            b = randn(5)
            c = similar(b)
            fill!(c,NaN)
            @test (c .= @~ A*b) ≈ A.args[1]*b - A.args[2]*b
        end

        @testset "gemv Float64" begin
            for A in (ApplyArray(-, randn(5,5), randn(5,5)),
                    ApplyArray(-, randn(5,5), view(randn(9, 5), 1:2:9, :))),
                b in (randn(5), view(randn(5),:), view(randn(5),1:5), view(randn(9),1:2:9))

                Ã = Array(A)
                c = similar(b)
                fill!(c,NaN)

                materialize!(MulAdd(1.0, A, b, 0.0, c))
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                c = similar(b)
                fill!(c,NaN)
                c .= @~ A*b
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                copyto!(c, @~ A*b)
                @test c ≈ Ã*b ≈ BLAS.gemv!('N', 1.0, Ã, b, 0.0, similar(c))

                b̃ = Array(b)
                copyto!(b̃, Mul(A,b̃))
                @test c ≈ b̃

                c .= @~ 2.0 * A * b
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 0.0, similar(c))

                c = Array(b)
                c .= @~ A*b + c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 1.0, Ã, b, 2.0, Array(b))

                c = Array(b)
                c .= @~ 2.0 * A*b + c
                @test c ≈ BLAS.gemv!('N', 2.0, Ã, b, 1.0, Array(b))

                c = Array(b)
                c .= @~ 3.0 * A*b + 2.0 * c
                @test c ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))

                d = similar(c)
                c = Array(b)
                d .= @~ 3.0 * A*b + 2.0 * c
                @test d ≈ BLAS.gemv!('N', 3.0, Ã, b, 2.0, Array(b))
            end
        end

        @testset "gemm" begin
            for A in (ApplyArray(-, randn(5,5), randn(5,5)),
                    ApplyArray(-, randn(5,5), view(randn(9, 5), 1:2:9, :))),
                B in (randn(5,5), view(randn(5,5),:,:), view(randn(5,5),1:5,:),
                    view(randn(5,5),1:5,1:5), view(randn(5,5),:,1:5))

                Ã = Array(A)
                C = similar(B)

                C .= @~ A*B
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 0.0, similar(C))

                B .= @~ A*B
                @test C ≈ B

                C .= @~ 2.0 * A*B
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 0.0, similar(C))

                C = Array(B)
                C .= @~ A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 1.0, Ã, B, 2.0, Array(B))

                C = Array(B)
                C .= @~ 2.0 * A*B + C
                @test C ≈ BLAS.gemm!('N', 'N', 2.0, Ã, B, 1.0, Array(B))

                C = Array(B)
                C .= @~ 3.0 * A*B + 2.0 * C
                @test C ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))

                d = similar(C)
                C = Array(B)
                d .= @~ 3.0 * A*B + 2.0 * C
                @test d ≈ BLAS.gemm!('N', 'N', 3.0, Ã, B, 2.0, Array(B))
            end
        end
    end

    @testset "view" begin
        A = ApplyArray(+, randn(5,5), randn(5,5))
        Ã = Array(A)
        V = view(A,2:3,2:4)
        @test MemoryLayout(typeof(V)) == ApplyLayout{typeof(+)}()
        b = randn(3)
        c = randn(2)
        materialize!(MulAdd(1.0, V, b, 0.0, c))
        @test c ≈ view(Ã,2:3,2:4)*b

        A = ApplyArray(-, randn(5,5), randn(5,5))
        Ã = Array(A)
        V = view(A,2:3,2:4)
        @test MemoryLayout(typeof(V)) == ApplyLayout{typeof(-)}()
        b = randn(3)
        c = randn(2)
        materialize!(MulAdd(1.0, V, b, 0.0, c))
        @test c ≈ view(Ã,2:3,2:4)*b
    end

    @testset "Broadcast add" begin
        a = randn(5)
        ã = reshape(a,5,1)
        b = randn(5)
        A = randn(5,5)
        c = [2.3]

        for op in (+, -)
            @testset "$op" begin
                @testset "vec .$op mat" begin
                    B = BroadcastArray(op, a, A)
                    @test @inferred(B * b) isa Vector
                    @test @inferred(B * A) isa Matrix
                    @test @inferred(A * B) isa Matrix
                    @test B * B isa Matrix
                    @test B * b ≈ op.(a, A) * b
                    @test B * A ≈ op.(a, A) * A
                    @test A * B ≈ A * op.(a, A)
                    @test B * B ≈ op.(a, A) * op.(a, A)
                end
                @testset "mat .$op vec" begin
                    B = BroadcastArray(op, A, a)
                    @test B * b ≈ op.(A, a) * b
                    @test B * A ≈ op.(A, a) * A
                    @test A * B ≈ A * op.(A, a)
                end
                @testset "mat .$op mat" begin
                    B = BroadcastArray(op, A, 2A)
                    @test B * b ≈ op.(A, 2A) * b
                    @test B * A ≈ op.(A, 2A) * A
                    @test A * B ≈ A * op.(A, 2A)
                end
                @testset "vecmat .$op mat" begin
                    B = BroadcastArray(op, ã, A)
                    @test B * b ≈ op.(ã, A) * b
                    @test B * A ≈ op.(ã, A) * A
                    @test A * B ≈ A * op.(ã, A)
                end
                @testset "mat .$op vecmat" begin
                    B = BroadcastArray(op, A, ã)
                    @test B * b ≈ op.(A, ã) * b
                    @test B * A ≈ op.(A, ã) * A
                    @test A * B ≈ A * op.(A, ã)
                end
                @testset "rowvec .$op mat" begin
                    B = BroadcastArray(op, a', A)
                    @test B * b ≈ op.(a', A) * b
                    @test B * A ≈ op.(a', A) * A
                    @test A * B ≈ A * op.(a', A)
                end
                @testset "constvec .$op mat" begin
                    B = BroadcastArray(op, c, A)
                    B̃ = BroadcastArray(op, c[1], A)
                    @test B * b ≈ B̃ * b ≈ op.(c, A) * b
                    @test B * A ≈ B̃ * A ≈ op.(c, A) * A
                    @test A * B ≈ A * B̃ ≈ A * op.(c, A)
                end
                @testset "constvec .$op vec" begin
                    B = BroadcastArray(op, c, b)
                    @test B * reshape(b,1,5) ≈ op.(c, b) * reshape(b,1,5)
                end
            end
        end

        @testset "Mixed" begin
            B = BroadcastArray(+, A, 2A)
            C = BroadcastArray(-, A, 2A)
            @test B*C ≈ 3A * (-A)
            @test C*B ≈ (-A) * 3A
        end
    end
end


text:
export DryModel, EquilMoist

#### Moisture component in atmosphere model
abstract type MoistureModel end

vars_state_conservative(::MoistureModel, FT) = @vars()
vars_state_gradient(::MoistureModel, FT) = @vars()
vars_state_gradient_flux(::MoistureModel, FT) = @vars()
vars_state_auxiliary(::MoistureModel, FT) = @vars()

function atmos_nodal_update_auxiliary_state!(
    ::MoistureModel,
    m::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function flux_moisture!(
    ::MoistureModel,
    ::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_flux!(
    ::MoistureModel,
    diffusive,
    ∇transform,
    state,
    aux,
    t,
) end
function flux_second_order!(
    ::MoistureModel,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_t,
) end
function flux_first_order!(
    ::MoistureModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
) end
function compute_gradient_argument!(
    ::MoistureModel,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
) end

internal_energy(atmos::AtmosModel, state::Vars, aux::Vars) =
    internal_energy(atmos.moisture, atmos.orientation, state, aux)

@inline function internal_energy(
    moist::MoistureModel,
    orientation::Orientation,
    state::Vars,
    aux::Vars,
)
    MoistThermodynamics.internal_energy(
        state.ρ,
        state.ρe,
        state.ρu,
        gravitational_potential(orientation, aux),
    )
end

temperature(atmos::AtmosModel, ::MoistureModel, state::Vars, aux::Vars) =
    air_temperature(thermo_state(atmos, state, aux))
pressure(atmos::AtmosModel, ::MoistureModel, state::Vars, aux::Vars) =
    air_pressure(thermo_state(atmos, state, aux))
soundspeed(atmos::AtmosModel, ::MoistureModel, state::Vars, aux::Vars) =
    soundspeed_air(thermo_state(atmos, state, aux))

@inline function total_specific_enthalpy(
    atmos::AtmosModel,
    moist::MoistureModel,
    state::Vars,
    aux::Vars,
)
    phase = thermo_state(atmos, state, aux)
    R_m = gas_constant_air(phase)
    T = air_temperature(phase)
    e_tot = state.ρe * (1 / state.ρ)
    return e_tot + R_m * T
end

"""
    DryModel

Assumes the moisture components is in the dry limit.
"""
struct DryModel <: MoistureModel end

vars_state_auxiliary(::DryModel, FT) = @vars(θ_v::FT, air_T::FT)
@inline function atmos_nodal_update_auxiliary_state!(
    moist::DryModel,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    e_int = internal_energy(atmos, state, aux)
    ts = PhaseDry(atmos.param_set, e_int, state.ρ)
    aux.moisture.θ_v = virtual_pottemp(ts)
    aux.moisture.air_T = air_temperature(ts)
    nothing
end

thermo_state(atmos::AtmosModel, state::Vars, aux::Vars) =
    thermo_state(atmos, atmos.moisture, state, aux)

function thermo_state(
    atmos::AtmosModel,
    moist::DryModel,
    state::Vars,
    aux::Vars,
)
    return PhaseDry(
        atmos.param_set,
        internal_energy(atmos, state, aux),
        state.ρ,
    )
end

"""
    EquilMoist

Assumes the moisture components are computed via thermodynamic equilibrium.
"""
struct EquilMoist{FT} <: MoistureModel
    maxiter::Int
    tolerance::FT
end
EquilMoist{FT}(;
    maxiter::IT = 3,
    tolerance::FT = FT(1e-1),
) where {FT <: AbstractFloat, IT <: Int} = EquilMoist{FT}(maxiter, tolerance)

vars_state_conservative(::EquilMoist, FT) = @vars(ρq_tot::FT)
vars_state_gradient(::EquilMoist, FT) = @vars(q_tot::FT)
vars_state_gradient_flux(::EquilMoist, FT) = @vars(∇q_tot::SVector{3, FT})
vars_state_auxiliary(::EquilMoist, FT) =
    @vars(temperature::FT, θ_v::FT, q_liq::FT)

@inline function atmos_nodal_update_auxiliary_state!(
    moist::EquilMoist,
    atmos::AtmosModel,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ps = atmos.param_set
    e_int = internal_energy(atmos, state, aux)
    ts = PhaseEquil(
        ps,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        moist.maxiter,
        moist.tolerance,
    )
    aux.moisture.temperature = air_temperature(ts)
    aux.moisture.θ_v = virtual_pottemp(ts)
    aux.moisture.q_liq = PhasePartition(ts).liq
    nothing
end

function thermo_state(
    atmos::AtmosModel,
    moist::EquilMoist,
    state::Vars,
    aux::Vars,
)
    e_int = internal_energy(atmos, state, aux)
    ps = atmos.param_set
    PS = typeof(ps)
    FT = eltype(state)
    return PhaseEquil{FT, PS}(
        ps,
        e_int,
        state.ρ,
        state.moisture.ρq_tot / state.ρ,
        aux.moisture.temperature,
    )
end

function compute_gradient_argument!(
    moist::EquilMoist,
    transform::Vars,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρinv = 1 / state.ρ
    transform.moisture.q_tot = state.moisture.ρq_tot * ρinv
end

function compute_gradient_flux!(
    moist::EquilMoist,
    diffusive::Vars,
    ∇transform::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    # diffusive flux of q_tot
    diffusive.moisture.∇q_tot = ∇transform.moisture.q_tot
end

function flux_moisture!(
    moist::EquilMoist,
    atmos::AtmosModel,
    flux::Grad,
    state::Vars,
    aux::Vars,
    t::Real,
)
    ρ = state.ρ
    u = state.ρu / ρ
    flux.moisture.ρq_tot += u * state.moisture.ρq_tot
end

function flux_second_order!(
    moist::EquilMoist,
    flux::Grad,
    state::Vars,
    diffusive::Vars,
    aux::Vars,
    t::Real,
    D_t,
)
    d_q_tot = (-D_t) .* diffusive.moisture.∇q_tot
    flux_second_order!(moist, flux, state, d_q_tot)
end
#TODO: Consider whether to not pass ρ and ρu (not state), foc BCs reasons
function flux_second_order!(moist::EquilMoist, flux::Grad, state::Vars, d_q_tot)
    flux.ρ += d_q_tot * state.ρ
    flux.ρu += d_q_tot .* state.ρu'
    flux.moisture.ρq_tot += d_q_tot * state.ρ
end


text:
include("header.jl")

@testset "nonlinear" begin
     x = randn(10,2)
     s = Sigm()
     @test s(x) == sigm.(x)
     r = ReLU()
     @test r(x) == relu.(x)
     lgp = LogSoftMax()
     @test lgp(x) == logp(x)
     sft = SoftMax()
     @test sft(x) == softmax(x)
     elu = ELU()
     elu(x)
     @test true
end


text:
#==
This is the modified ZDT example taken from
Esmaile Khorram, K Khaledian, and M Khaledyan.
A numerical method for
constructing the pareto front of multi-objective optimization problems.
Journal of Computational and Applied Mathematics, 261:158–171, 2014.

This is an example of a disconnected Pareto front.
The standard equality-based NBI method fails miserably here, and
will hit an infeasibility error somewhere along the way.
It will be interesting to see how it works with
the inequality-extended NBI method
==#

using JuMP, MultiJuMP, Ipopt

m = multi_model(solver=IpoptSolver())
n = 30

l = -ones(n); l[1] = 0
u = ones(n)
@variable(m, l[i] <= x[i=1:n] <= u[i])
@NLexpression(m, f1, x[1])
@NLexpression(m, g, 1 + 9*sum{x[j]^2, j=2:n}/(n-1))
@NLexpression(m, h, 1 - sqrt((f1/g)) - (f1/g)*sin(10*π*f1))
@NLexpression(m, f2, g*h)

iv1 = zeros(n)
obj1 = SingleObjective(f1, sense = :Min,
                       iv = Dict{Symbol, Any}(:x => iv1))

# As the problem is nonconvex, we have to supply the
# initial value to get the global minimum of f2
iv2 = zeros(n)
iv2[1] = 0.85
obj2 = SingleObjective(f2, sense = :Min,
                       iv = Dict{Symbol, Any}(:x => iv2))

multim = get_multidata(m)
multim.objectives = [obj1, obj2]
multim.pointsperdim = 60
# NB: solve fails with an infeasibility error 65% through the algorithm,
# You can still plot the results by pasting in the plot_comparison() function below
solve(m, method = NBI(true))

function pltfun(x1)
    setvalue(x[1], x1)
    setvalue(x[2:n], zeros(n-1))
    getvalue(f2)
end

using Plots
plt = plot(multim, label="NBI points",
           legend = :bottomleft)
xvals = linspace(0,1,200)
plot!(plt, xvals, pltfun, markercolor="orange",
      label="Image boundary")


text:
include("solver_preamble.jl")
import SCS
factory = with_optimizer(SCS.Optimizer, verbose=0)
config = MOI.Test.TestConfig(atol=1e-5, rtol=1e-5)
@testset "Linear" begin
    Tests.linear_test(factory, config)
end
@testset "SOC" begin
    Tests.soc_test(factory, config)
end
@testset "SDP" begin
    Tests.sd_test(factory, config)
end


text:
module ViewHelper

using Genie, Validation, Helpers, App

export output_errors, output_flash

function output_errors(m, field::Symbol) :: String
  Validation.has_errors_for(m, field) ? """$(Validation.errors_to_string(m, field, "
\n", upper_case_first = true))""" : ""
end

function output_flash(params::Dict{Symbol,Any}) :: String
  ! isempty( flash(params) ) ? """
$(flash(params))
""" : ""
end

end

max_stars_repo_name: essenciary/genie-todo-mvc

max_issues_repo_name: essenciary/genie-todo-mvc

max_forks_repo_name: essenciary/genie-todo-mvc


text:
if haskey(ENV, "GIZA_JL_DLL") && ENV["GIZA_JL_DLL"] != ""
    _GIZALIB = ENV["GIZA_JL_DLL"]
else
    _GIZALIB = "libgiza.so"
end
open(joinpath(@__DIR__, "deps.jl"), "w") do io
    write(io, "const _GIZALIB = \"$_GIZALIB\"\n")
end
nothing


text:
module Writer

export write_pddl, write_domain, write_problem
export save_domain, save_problem

using Julog
using ..PDDL: Domain, Problem, Action, Event
using ..PDDL: DEFAULT_REQUIREMENTS, IMPLIED_REQUIREMENTS

"Write list of typed formulae in PDDL syntax."
function write_typed_list(formulae::Vector{<:Term}, types::Vector{Symbol})
    if all(x -> x == :object, types)
        return join(write_subformula.(formulae), " ")
    end
    str = ""
    for (i, (f, t)) in enumerate(zip(formulae, types))
        str *= "$(write_subformula(f)) "
        if i == length(formulae) || t != types[i+1]
            str *= "- $t "
        end
    end
    return str[1:end-1]
end

function write_typed_list(formulae::Vector{<:Term})
    formulae, types = zip([(t.args[1], t.name) for t in formulae]...)
    return write_typed_list(collect(formulae), collect(types))
end

function write_typed_list(formulae::Vector{<:Term}, types::Dict{<:Term,Symbol})
    types = Symbol[types[f] for f in formulae]
    return write_typed_list(formulae, types)
end

"Indent list of PDDL constants."
function indent_const_list(str::String, indent::Int, maxchars::Int=80)
    if length(str) < maxchars - indent return str end
    lines = String[]
    while length(str) > 0
        l_end = min(maxchars - indent, length(str))
        if l_end != length(str) l_end = findlast(' ', str[1:l_end]) end
        if isnothing(l_end) l_end = findfirst(' ', str) end
        if isnothing(l_end) l_end = length(str) end
        push!(lines, str[1:l_end])
        str = str[l_end+1:end]
    end
    return join(lines, "\n" * ' '^indent)
end

"Indent list of typed PDDL constants."
function indent_typed_list(str::String, indent::Int, maxchars::Int=80)
    if length(str) < maxchars - indent return str end
    if !occursin(" - ", str) return indent_const_list(strs[1]) end
    substrs = String[]
    while length(str) > 0
        idxs = findnext(r" - (\S+)", str, 1)
        if isnothing(idxs) idxs = [length(str)] end
        push!(substrs, str[1:idxs[end]])
        str = strip(str[idxs[end]+1:end])
    end
    substrs = indent_const_list.(substrs, indent, maxchars)
    return join(substrs, "\n" * ' '^indent)
end

"Write formula in PDDL syntax."
function write_formula(f::Compound)
    if f.name in [:exists, :forall]
        typecond, body = f.args
        var_str = write_typed_list(flatten_conjs(typecond))
        body_str = write_formula(body)
        return "($(f.name) ($var_str) $body_str)"
    elseif f.name in [:and, :or, :not, :when, :imply,
                      :(==), :>, :<, :!=, :>=, :<=, :+, :-, :*, :/,
                      :assign, :increase, :decrease,
                      Symbol("scale-up"), Symbol("scale-down")]
        name = f.name == :(==) ? "=" : string(f.name)
        args = join([write_formula(a) for a in f.args], " ")
        return "($name $args)"
    else
        args = join([write_subformula(a) for a in f.args], " ")
        return "($(f.name) $args)"
    end
end
write_formula(f::Var) = "?" * lowercasefirst(repr(f))
write_formula(f::Const) = "(" * repr(f) * ")"

write_subformula(f::Compound) = write_formula(f)
write_subformula(f::Var) = "?" * lowercasefirst(repr(f))
write_subformula(f::Const) = repr(f)

"Write to string in PDDL syntax."
write_pddl(f::Term) = write_formula(f)

"Write domain in PDDL syntax."
function write_domain(domain::Domain, indent::Int=2)
    strs = Dict{Symbol,String}()
    fields = [:requirements, :types, :constants, :predicates, :functions]
    strs[:requirements] = write_requirements(domain.requirements)
    strs[:types] = write_types(domain.types)
    strs[:constants] = indent_typed_list(
        write_typed_list(domain.constants, domain.constypes), 14)
    strs[:predicates] = write_predicates(domain.predicates, domain.predtypes)
    strs[:functions] = write_predicates(domain.functions, domain.functypes)
    strs = ["($(repr(k)) $(strs[k]))" for k in fields
            if haskey(strs, k) && length(strs[k]) > 0]
    append!(strs, write_axiom.(values(domain.axioms)))
    append!(strs, write_action.(values(domain.actions), 3))
    append!(strs, write_event.(values(domain.events), 3))
    pushfirst!(strs, "(define (domain $(domain.name))")
    return join(strs, "\n" * ' '^indent) * "\n)"
end
write_pddl(domain::Domain) = write_domain(domain)

"Write domain requirements."
function write_requirements(requirements::Dict{Symbol,Bool})
    return join([":$k" for (k, v) in requirements if v], " ")
end

"Write domain types."
function write_types(types::Dict{Symbol,Vector{Symbol}})
    strs = Dict{Symbol,String}()
    for (type, subtypes) in types
        if isempty(subtypes) || type == :object continue end
        subtype_str = join(subtypes, " ")
        strs[type] = "$subtype_str - $type"
    end
    maxtypes = collect(setdiff(types[:object], keys(strs)))
    strs[:object] = join(maxtypes, " ")
    return strip(join(values(strs), " "))
end

"Write domain predicates or functions."
function write_predicates(predicates::Dict{Symbol,Term},
                          predtypes::Dict{Symbol,Vector{Symbol}})
    strs = String[]
    for (name, pred) in predicates
        types = predtypes[name]
        args_str = write_typed_list(pred.args, types)
        push!(strs, "($name $args_str)")
    end
    return join(strs, " ")
end

"Write PDDL axiom / derived predicate."
function write_axiom(c::Clause, key=":derived")
    head_str = write_formula(c.head)
    body_str = length(c.body) == 1 ? write_formula(c.body[1]) :
        "(and " * join(write_formula.(c.body), " ") * ")"
    return "($key $head_str $body_str)"
end

"Write action in PDDL syntax."
function write_action(action::Action, indent::Int=1)
    strs = Dict{Symbol,String}()
    fields = [:action, :parameters, :precondition, :effect]
    strs[:action] = string(action.name)
    strs[:parameters] = "(" * write_typed_list(action.args, action.types) * ")"
    strs[:precondition] = write_formula(action.precond)
    strs[:effect] = write_formula(action.effect)
    strs = ["$(repr(k)) $(strs[k])" for k in fields]
    return "(" * join(strs, "\n" * ' '^indent) * ")"
end
write_pddl(action::Action) = write_action(action)

"Write event in PDDL syntax."
function write_event(event::Event, indent::Int=1)
    strs = Dict{Symbol,String}()
    fields = [:event, :precondition, :effect]
    strs[:event] = string(action.name)
    strs[:precondition] = write_formula(action.precond)
    strs[:effect] = write_formula(action.effect)
    strs = ["$(repr(k)) $(strs[k])" for k in fields]
    return "(" * join(strs, "\n" * ' '^indent) * ")"
end
write_pddl(event::Event) = write_event(event)

"Write problem in PDDL syntax."
function write_problem(problem::Problem, indent::Int=2)
    strs = Dict{Symbol,String}()
    fields = [:domain, :objects, :init, :goal, :metric]
    strs[:domain] = string(problem.domain)
    strs[:objects] = indent_typed_list(
        write_typed_list(problem.objects, problem.objtypes), 12)
    strs[:init] = write_init(problem.init, 9)
    strs[:goal] = write_formula(problem.goal)
    strs[:metric] = write_metric(problem.metric)
    strs = ["($(repr(k)) $(strs[k]))" for k in fields
            if haskey(strs, k) && length(strs[k]) > 0]
    pushfirst!(strs, "(define (problem $(problem.name))")
    return join(strs, "\n" * ' '^indent) * "\n)"
end
write_pddl(problem::Problem) = write_problem(problem)

"Write initial problem formulae in PDDL syntax."
function write_init(init::Vector{<:Term}, indent::Int=2, maxchars::Int=80)
    strs = write_formula.(init)
    if sum(length.(strs)) + length("(:init )") < maxchars
        return join(strs, " ") end
    return join(strs, "\n" * ' '^indent)
end

"Write metric for PDDL problem."
function write_metric(metric::Tuple{Int,Term})
    sign, formula = metric
    direction = sign > 0 ? "maximize" : "minimize"
    return direction * " " * write_formula(formula)
end
write_metric(::Nothing) = ""

"Save PDDL domain to specified path."
function save_domain(path::String, domain::Domain)
    open(f->write(f, write_domain(domain)), path, "w")
    return path
end

"Save PDDL problem to specified path."
function save_problem(path::String, problem::Problem)
    open(f->write(f, write_problem(problem)), path, "w")
    return path
end

end


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

mutable struct ListAccountSasResponse <: SwaggerModel
    accountSasToken::Any # spec type: Union{ Nothing, String } # spec name: accountSasToken

    function ListAccountSasResponse(;accountSasToken=nothing)
        o = new()
        validate_property(ListAccountSasResponse, Symbol("accountSasToken"), accountSasToken)
        setfield!(o, Symbol("accountSasToken"), accountSasToken)
        o
    end
end # type ListAccountSasResponse

const _property_map_ListAccountSasResponse = Dict{Symbol,Symbol}(Symbol("accountSasToken")=>Symbol("accountSasToken"))
const _property_types_ListAccountSasResponse = Dict{Symbol,String}(Symbol("accountSasToken")=>"String")
Base.propertynames(::Type{ ListAccountSasResponse }) = collect(keys(_property_map_ListAccountSasResponse))
Swagger.property_type(::Type{ ListAccountSasResponse }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_ListAccountSasResponse[name]))}
Swagger.field_name(::Type{ ListAccountSasResponse }, property_name::Symbol) =  _property_map_ListAccountSasResponse[property_name]

function check_required(o::ListAccountSasResponse)
    true
end

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


text:
# Stocks/Profiles

"""
    profile(conn::Connection, symbol::String)

Get company profile
"""
function company(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "company", path_params; symbol = symbol)
end

"""
    insider_roster(conn::Connection, symbol::String)

Returns the top 10 insiders, with the most recent information.
"""
function insider_roster(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "insider-roster", path_params; symbol = symbol)
end

"""
    insider_summary(conn::Connection, symbol::String)

Returns aggregated insiders summary data for the last 6 months.
"""
function insider_summary(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "insider-summary", path_params; symbol = symbol)
end

"""
    insider_transactions(conn::Connection, symbol::String)

Returns insider transactions.
"""
function insider_transactions(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "insider_transactions", path_params; symbol = symbol)
end

"""
    logo(conn::Connection, symbol::String)

Get company logo
"""
function logo(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "logo", path_params; symbol = symbol)
end

"""
    peers(conn::Connection, symbol::String)

Get company peer group
"""
function peers(conn::Connection, symbol::String)
    path_params = [:symbol]
    return get_stocks(conn, "peers", path_params; symbol = symbol)
end

max_stars_repo_name: GodolphinCapitalManagement/Iex

max_issues_repo_name: GodolphinCapitalManagement/Iex

max_forks_repo_name: GodolphinCapitalManagement/Iex


text:
using Distributed

@everywhere begin
    using DataStructures
    using Test
    using ParallelUtilities
    using ParallelUtilities.ClusterQueryUtils
    using ParallelUtilities.ClusterQueryUtils: oneworkerpernode
    using OffsetArrays
    import ParallelUtilities: BinaryTreeNode, BranchChannel,
    OrderedBinaryTree, SegmentedOrderedBinaryTree,
    parentnoderank, nchildren,
    createbranchchannels,
    workersactive,
    leafrankfoldedtree,
    TopTreeNode, SubTreeNode,
    NoSplat, reducedvalue

    function parentnoderank(tree::SegmentedOrderedBinaryTree, i::Integer)
        1 <= i <= length(tree) || throw(BoundsError(tree, i))

        Nmasternodes = length(keys(ParallelUtilities.workersonhosts(tree)))
        toptreenonleafnodes = length(tree.toptree) - Nmasternodes

        if toptreenonleafnodes == 0
            pr = parentnoderank(ParallelUtilities.unsegmentedtree(tree),i)

        elseif i <= toptreenonleafnodes
            #= In a SegmentedSequentialBinaryTree the leading indices
            are the parent nodes of the top tree, so ind = i
            In a SegmentedOrderedBinaryTree, the leaves are removed
            from the top tree, so only even numbers are left.
            In this case, index i of the full tree refers to index 2i of the
            top tree, so ind = 2i
            =#
            ind = ParallelUtilities.fulltree_to_toptree_index(tree.toptree, i)
            p = tree.toptree[ind].p
            #  Compute the parent of the node with rank ind on the top tree.
            # In a SegmentedSequentialBinaryTree this is what we want.
            # In a SegmentedOrderedBinaryTree, we need to convert this back to
            # the index of the full tree, that is div(pr, 2)
            pr_top = parentnoderank(tree.toptree, ind)
            pr = ParallelUtilities.toptree_to_fulltree_index(tree.toptree, pr_top)
        else
            subtree, rankinsubtree, nptotalprevhosts = ParallelUtilities.subtree_rank(tree, i)

            if rankinsubtree == ParallelUtilities.topnoderank(subtree)
                # masternode
                # parent will be on the top - tree
                p = subtree[rankinsubtree].p
                leafno = ParallelUtilities.masternodeindex(tree, p)
                Nmasternodes = length(keys(ParallelUtilities.workersonhosts(tree)))
                leafrank = ParallelUtilities.leafrankfoldedtree(tree.toptree, Nmasternodes, leafno)
                pr_top = parentnoderank(tree.toptree, leafrank)
                # Convert back to the rank on the full tree where the
                # leaves of the top tree aren't stored.
                pr = ParallelUtilities.toptree_to_fulltree_index(tree.toptree, pr_top)
            else
                # node on a sub - tree
                pr = parentnoderank(subtree, rankinsubtree)
                pr += nptotalprevhosts + toptreenonleafnodes
            end
        end

        return pr
    end
end

macro testsetwithinfo(str, ex)
    quote
        @info "Testing "*$str
        @testset $str begin $(esc(ex)); end;
    end
end

fmap_local(x) = x^2
fred_local(x) = x
fred_local(x, y) = x + y

function showworkernumber(ind, nw)
    # Cursor starts off at the beginning of the line
    print("\u1b[K") # clear till end of line
    print("Testing on worker $ind of $nw")
    # return the cursor to the beginning of the line
    endchar = ind == nw ? "\n" : "\r"
    print(endchar)
end

@testsetwithinfo "utilities" begin
    @testset "hostnames" begin
        hosts = hostnames()
        nodes = unique(hosts)
        @test nodenames() == nodes
        @test nodenames(hosts) == nodes
        np1 = nprocs_node(hosts, nodes)
        np2 = nprocs_node(hosts)
        np3 = nprocs_node()
        @test np1 == np2 == np3
        for node in nodes
            npnode = count(isequal(node), hosts)
            @test np1[node] == npnode
        end
        p1 = procs_node(workers(), hosts, nodes)
        for node in nodes
            pnode = workers()[findall(isequal(node), hosts)]
            @test p1[node] == pnode
        end
        np4 = nprocs_node(p1)
        @test np1 == np4
    end
end;

@testset "BinaryTree" begin
    @testsetwithinfo "BinaryTreeNode" begin
        @testset "Constructor" begin
            p = workers()[1]
            b = BinaryTreeNode(p, p, 0)
            @test nchildren(b) == 0
            b = BinaryTreeNode(p, p, 1)
            @test nchildren(b) == 1
            b = BinaryTreeNode(p, p, 2)
            @test nchildren(b) == 2

            @test_throws DomainError BinaryTreeNode(p, p, 3)
            @test_throws DomainError BinaryTreeNode(p, p,-1)
        end
    end

    @testsetwithinfo "BinaryTree" begin
        @testsetwithinfo "OrderedBinaryTree" begin
            @testset "pid and parent" begin
                for imax = 1:100
                    procs = 1:imax
                    tree = OrderedBinaryTree(procs)
                    @test length(tree) == length(procs)

                    topnoderank = ParallelUtilities.topnoderank(tree)
                    @test tree[topnoderank].parent == topnoderank
                    for rank in 1:length(tree)
                        node = tree[rank]
                        @test node.p == procs[rank]
                        @test node.parent == procs[parentnoderank(tree, rank)]
                    end
                    @test_throws BoundsError(tree, 0) parentnoderank(tree, 0)
                    @test_throws BoundsError(tree, imax + 1) parentnoderank(tree, imax + 1)
                end
            end

            @testset "nchildren" begin
                tree = OrderedBinaryTree(1:1)
                @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                @test_throws BoundsError(tree, 2) nchildren(tree, 2)
                @test ParallelUtilities.topnoderank(tree) == 1

                tree = OrderedBinaryTree(1:2)
                @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 1
                @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                @test_throws BoundsError(tree, 3) nchildren(tree, 3)
                @test ParallelUtilities.topnoderank(tree) == 2

                tree = OrderedBinaryTree(1:8)
                @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 2
                @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 0
                @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 2
                @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 1
                @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                @test_throws BoundsError(tree, 9) nchildren(tree, 9)
                @test ParallelUtilities.topnoderank(tree) == 8

                tree = OrderedBinaryTree(1:11)
                @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 2
                @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 0
                @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 2
                @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 2
                @test nchildren(tree, 9) == nchildren(tree[9]) == tree[9].nchildren == 0
                @test nchildren(tree, 10) == nchildren(tree[10]) == tree[10].nchildren == 2
                @test nchildren(tree, 11) == nchildren(tree[11]) == tree[11].nchildren == 0
                @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                @test_throws BoundsError(tree, 12) nchildren(tree, 12)
                @test ParallelUtilities.topnoderank(tree) == 8

                tree = OrderedBinaryTree(1:13)
                @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 2
                @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 0
                @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 2
                @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 2
                @test nchildren(tree, 9) == nchildren(tree[9]) == tree[9].nchildren == 0
                @test nchildren(tree, 10) == nchildren(tree[10]) == tree[10].nchildren == 2
                @test nchildren(tree, 11) == nchildren(tree[11]) == tree[11].nchildren == 0
                @test nchildren(tree, 12) == nchildren(tree[12]) == tree[12].nchildren == 2
                @test nchildren(tree, 13) == nchildren(tree[13]) == tree[13].nchildren == 0
                @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                @test_throws BoundsError(tree, 14) nchildren(tree, 14)
                @test ParallelUtilities.topnoderank(tree) == 8
            end

            @testset "level" begin
                tree = OrderedBinaryTree(1:15)
                @test ParallelUtilities.levels(tree) == 4

                @test ParallelUtilities.levelfromtop.((tree,), 1:2:15) == ones(Int, 8).*4
                @test ParallelUtilities.levelfromtop.((tree,), (2, 6, 10, 14)) == (3, 3, 3, 3)
                @test ParallelUtilities.levelfromtop.((tree,), (4, 12)) == (2, 2)
                @test ParallelUtilities.levelfromtop(tree, 8) == 1
                for p in [0, length(tree) + 1]
                    @test_throws BoundsError(tree, p) ParallelUtilities.levelfromtop(tree, p)
                end

                tree = OrderedBinaryTree(1:13)
                @test ParallelUtilities.levels(tree) == 4
                @test ParallelUtilities.levelfromtop.((tree,), 1:2:11) == ones(Int, 6).*4
                @test ParallelUtilities.levelfromtop.((tree,), (2, 6, 10, 13)) == (3, 3, 3, 3)
                @test ParallelUtilities.levelfromtop.((tree,), (4, 12)) == (2, 2)
                @test ParallelUtilities.levelfromtop(tree, 8) == 1
                for p in [0, length(tree) + 1]
                    @test_throws BoundsError(tree, p) ParallelUtilities.levelfromtop(tree, p)
                end
            end
        end

        @testsetwithinfo "SegmentedOrderedBinaryTree" begin
            @testsetwithinfo "single host" begin
                @testset "pid and parent" begin
                    for imax = 1:100
                        procs = 1:imax
                        workersonhosts = Dict("host" => procs)
                        tree = SegmentedOrderedBinaryTree(procs, workersonhosts)
                        treeOBT = OrderedBinaryTree(procs)
                        @test length(tree) == length(procs) == length(treeOBT)

                        topnoderank = ParallelUtilities.topnoderank(tree)
                        # The top node is its own parent
                        @test tree[topnoderank].parent == topnoderank
                        @test tree[topnoderank] == ParallelUtilities.topnode(tree)
                        for rank in 1:length(tree)
                            node = tree[rank]
                            parentnode = tree[parentnoderank(tree, rank)]
                            @test length(procs) > 1 ? nchildren(parentnode) > 0 : nchildren(parentnode) == 0
                            @test node.p == procs[rank]
                            @test node.parent == procs[parentnoderank(treeOBT, rank)]
                            @test parentnode.p == node.parent
                        end
                    end
                end;

                @testset "nchildren" begin
                    procs = 1:1
                    tree = SegmentedOrderedBinaryTree(procs, Dict("host" => procs))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 2) nchildren(tree, 2)
                    @test ParallelUtilities.topnoderank(tree) == 1

                    procs = 1:2
                    tree = SegmentedOrderedBinaryTree(procs, Dict("host" => procs))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                    @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 1
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 3) nchildren(tree, 3)
                    @test ParallelUtilities.topnoderank(tree) == 2

                    procs = 1:8
                    tree = SegmentedOrderedBinaryTree(procs, Dict("host" => procs))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                    @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                    @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                    @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 2
                    @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 0
                    @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 2
                    @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                    @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 1
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 9) nchildren(tree, 9)
                    @test ParallelUtilities.topnoderank(tree) == 8

                    procs = 1:11
                    tree = SegmentedOrderedBinaryTree(procs, Dict("host" => procs))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                    @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                    @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                    @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 2
                    @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 0
                    @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 2
                    @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                    @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 2
                    @test nchildren(tree, 9) == nchildren(tree[9]) == tree[9].nchildren == 0
                    @test nchildren(tree, 10) == nchildren(tree[10]) == tree[10].nchildren == 2
                    @test nchildren(tree, 11) == nchildren(tree[11]) == tree[11].nchildren == 0
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 12) nchildren(tree, 12)
                    @test ParallelUtilities.topnoderank(tree) == 8

                    procs = 1:13
                    tree = SegmentedOrderedBinaryTree(procs, Dict("host" => procs))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 0
                    @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                    @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                    @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 2
                    @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 0
                    @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 2
                    @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                    @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 2
                    @test nchildren(tree, 9) == nchildren(tree[9]) == tree[9].nchildren == 0
                    @test nchildren(tree, 10) == nchildren(tree[10]) == tree[10].nchildren == 2
                    @test nchildren(tree, 11) == nchildren(tree[11]) == tree[11].nchildren == 0
                    @test nchildren(tree, 12) == nchildren(tree[12]) == tree[12].nchildren == 2
                    @test nchildren(tree, 13) == nchildren(tree[13]) == tree[13].nchildren == 0
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 14) nchildren(tree, 14)
                    @test ParallelUtilities.topnoderank(tree) == 8
                end;
            end;

            @testsetwithinfo "multiple hosts" begin
                @testset "length" begin
                    procs = 1:2
                    tree = SegmentedOrderedBinaryTree(procs,
                        OrderedDict("host1" => 1:1,"host2" => 2:2))
                    @test length(tree) == 2 + 1

                    procs = 1:4
                    tree = SegmentedOrderedBinaryTree(procs,
                        OrderedDict("host1" => 1:2,"host2" => 3:4))

                    @test length(tree) == 4 + 1

                    procs = 1:12
                    tree = SegmentedOrderedBinaryTree(procs,
                        OrderedDict(
                            "host1" => 1:3,"host2" => 4:6,
                            "host3" => 7:9,"host4" => 10:12))

                    @test length(tree) == 12 + 3
                end;

                @testset "leafrankfoldedtree" begin
                    treeflag = OrderedBinaryTree(1:1)
                    @test leafrankfoldedtree(treeflag, 5, 1) == 1
                    @test leafrankfoldedtree(treeflag, 5, 2) == 3
                    @test leafrankfoldedtree(treeflag, 5, 3) == 5
                    @test leafrankfoldedtree(treeflag, 5, 4) == 7
                    @test leafrankfoldedtree(treeflag, 5, 5) == 9
                end;

                @testset "pid and parent" begin
                    for imax = 2:100
                        procs = 1:imax
                        mid = div(imax, 2)
                        workersonhosts = OrderedDict{String, Vector{Int}}()
                        workersonhosts["host1"] = procs[1:mid]
                        workersonhosts["host2"] = procs[mid + 1:end]
                        tree = SegmentedOrderedBinaryTree(procs, workersonhosts)

                        top = ParallelUtilities.topnoderank(tree)
                        @test tree[top] == ParallelUtilities.topnode(tree)
                        for (ind, rank) in enumerate(1:mid)
                            node = tree[rank + 1]
                            parentnode = tree[parentnoderank(tree, rank + 1)]
                            @test parentnode.p == node.parent
                            pnodes = workersonhosts["host1"]
                            @test node.p == pnodes[ind]
                            OBT = OrderedBinaryTree(pnodes)
                            if ind == ParallelUtilities.topnoderank(OBT)
                                # Special check for 2 hosts as
                                # there's only one node in the top tree
                                @test node.parent == ParallelUtilities.topnode(tree.toptree).p
                            else
                                @test node.parent == pnodes[parentnoderank(OBT, ind)]
                            end
                        end
                        for (ind, rank) in enumerate(mid + 1:imax)
                            node = tree[rank + 1]
                            parentnode = tree[parentnoderank(tree, rank + 1)]
                            @test parentnode.p == node.parent
                            pnodes = workersonhosts["host2"]
                            @test node.p == pnodes[ind]
                            OBT = OrderedBinaryTree(pnodes)
                            if ind == ParallelUtilities.topnoderank(OBT)
                                # Special check for 2 hosts as
                                # there's only one node in the top tree
                                @test node.parent == ParallelUtilities.topnode(tree.toptree).p
                            else
                                @test node.parent == pnodes[parentnoderank(OBT, ind)]
                            end
                        end
                    end
                end;

                @testset "nchildren" begin
                    procs = 1:2
                    tree = SegmentedOrderedBinaryTree(procs,
                        OrderedDict("host1" => 1:1,"host2" => 2:2))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 2
                    @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 0
                    @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 0
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 4) nchildren(tree, 4)

                    procs = 1:12
                    tree = SegmentedOrderedBinaryTree(procs,
                        OrderedDict(
                            "host1" => 1:3,"host2" => 4:6,
                            "host3" => 7:9,"host4" => 10:12))
                    @test nchildren(tree, 1) == nchildren(tree[1]) == tree[1].nchildren == 2
                    @test nchildren(tree, 2) == nchildren(tree[2]) == tree[2].nchildren == 2
                    @test nchildren(tree, 3) == nchildren(tree[3]) == tree[3].nchildren == 2
                    @test nchildren(tree, 4) == nchildren(tree[4]) == tree[4].nchildren == 0
                    @test nchildren(tree, 5) == nchildren(tree[5]) == tree[5].nchildren == 2
                    @test nchildren(tree, 6) == nchildren(tree[6]) == tree[6].nchildren == 0
                    @test nchildren(tree, 7) == nchildren(tree[7]) == tree[7].nchildren == 0
                    @test nchildren(tree, 8) == nchildren(tree[8]) == tree[8].nchildren == 2
                    @test nchildren(tree, 9) == nchildren(tree[9]) == tree[9].nchildren == 0
                    @test nchildren(tree, 10) == nchildren(tree[10]) == tree[10].nchildren == 0
                    @test nchildren(tree, 11) == nchildren(tree[11]) == tree[11].nchildren == 2
                    @test nchildren(tree, 12) == nchildren(tree[12]) == tree[12].nchildren == 0
                    @test nchildren(tree, 13) == nchildren(tree[13]) == tree[13].nchildren == 0
                    @test nchildren(tree, 14) == nchildren(tree[14]) == tree[14].nchildren == 2
                    @test nchildren(tree, 15) == nchildren(tree[15]) == tree[15].nchildren == 0
                    @test_throws BoundsError(tree, 0) nchildren(tree, 0)
                    @test_throws BoundsError(tree, 16) nchildren(tree, 16)
                end;
            end;
        end
    end
end;

@testsetwithinfo "reduction" begin

    @testset "BranchChannel" begin
        @test_throws DomainError BranchChannel(1, 3)
        parentchannel = RemoteChannel(() -> Channel(1))
        @test_throws DomainError BranchChannel(1, parentchannel, 3)
    end

    @testset "TopTreeNode" begin
        # Special test for this as this is usually not called when tests are carried out on the same machine
        parentchannel = RemoteChannel(() -> Channel(1))
        childrenchannel = RemoteChannel(() -> Channel(2))
        pipe = ParallelUtilities.BranchChannel(1, parentchannel, childrenchannel, 2)

        put!(childrenchannel, ParallelUtilities.pval(1, false, 1))
        put!(childrenchannel, ParallelUtilities.pval(2, false, 2))

        redval = reducedvalue(+, ParallelUtilities.TopTreeNode(1), pipe, nothing)
        @test redval === ParallelUtilities.pval(1, false, 3)

        put!(childrenchannel, ParallelUtilities.pval(1, false, 1))
        put!(childrenchannel, ParallelUtilities.pval(2, true, nothing))

        redval = reducedvalue(+, ParallelUtilities.TopTreeNode(1), pipe, nothing)
        @test redval === ParallelUtilities.pval(1, true, nothing)

        put!(childrenchannel, ParallelUtilities.pval(1, false, 1))
        put!(childrenchannel, ParallelUtilities.pval(2, false, 2))
        @test_throws Exception reducedvalue(x -> error(""), ParallelUtilities.TopTreeNode(1), pipe, nothing)
    end

    @testset "fake multiple hosts" begin
        tree = ParallelUtilities.SegmentedOrderedBinaryTree([1,1], OrderedDict("host1" => 1:1, "host2" => 1:1))
        branches = ParallelUtilities.createbranchchannels(tree)
        @test ParallelUtilities.pmapreduceworkers(x -> 1, +, (tree, branches), (1:4,)) == 4

        if nworkers() > 1
            p = procs_node()
            # Choose workers on the same node to avoid communication bottlenecks in testing
            w = first(values(p))
            tree = ParallelUtilities.SegmentedOrderedBinaryTree(w, OrderedDict("host1" => w[1]:w[1], "host2" => w[2]:w[end]))
            branches = ParallelUtilities.createbranchchannels(tree)
            @test ParallelUtilities.pmapreduceworkers(x -> 1, +, (tree, branches), (1:length(w),)) == length(w)
        end
    end
end

@testset "pmapreduce" begin
    @testsetwithinfo "pmapreduce" begin
        @testsetwithinfo "sum" begin
            @testsetwithinfo "comparison with mapreduce" begin
                for iterators in Any[(1:1,), (ones(2,2),), (1:10,)]
                    res_exp = mapreduce(x -> x^2, +, iterators...)
                    res = pmapreduce(x -> x^2, +, iterators...)
                    @test res_exp == res

                    res_exp = mapreduce(x -> x^2, +, iterators..., init = 100)
                    res = pmapreduce(x -> x^2, +, iterators..., init = 100)
                    @test res_exp == res
                end

                @testset "dictionary" begin
                    res = pmapreduce(x -> Dict(x => x), merge, 1:1)
                    res_exp = mapreduce(x -> Dict(x => x), merge, 1:1)
                    @test res == res_exp

                    res = pmapreduce(x -> Dict(x => x), merge, 1:200)
                    res_exp = mapreduce(x -> Dict(x => x), merge, 1:200)
                    @test res == res_exp

                    res = pmapreduce(x -> OrderedDict(x => x), merge, 1:20)
                    res_exp = mapreduce(x -> OrderedDict(x => x), merge, 1:20)
                    @test res == res_exp
                end

                iterators = (1:10, 2:2:20)
                res_exp = mapreduce((x, y) -> x*y, +, iterators...)
                res = pmapreduce((x, y) -> x*y, +, iterators...)
                @test res_exp == res

                res_exp = mapreduce((x, y) -> x*y, +, iterators..., init = 100)
                res = pmapreduce((x, y) -> x*y, +, iterators..., init = 100)
                @test res_exp == res

                iterators = (1:10, 2:2:20)
                iterators_product = Iterators.product(iterators...)
                res_exp = mapreduce(((x, y),) -> x*y, +, iterators_product)
                res = pmapreduce(((x, y),) -> x*y, +, iterators_product)
                @test res_exp == res

                res_exp_2itp = mapreduce(((x, y), (a, b)) -> x*a + y*b, +, iterators_product, iterators_product)
                res_2itp = pmapreduce(((x, y), (a, b)) -> x*a + y*b, +, iterators_product, iterators_product)
                @test res_2itp == res_exp_2itp

                iterators_product_putil = ParallelUtilities.product(iterators...)
                res_exp2 = mapreduce(((x, y),) -> x*y, +, iterators_product_putil)
                res2 = pmapreduce(((x, y),) -> x*y, +, iterators_product_putil)
                @test res_exp2 == res2
                @test res_exp2 == res_exp

                res_exp_2pup = mapreduce(((x, y), (a, b)) -> x*a + y*b, +, iterators_product_putil, iterators_product_putil)
                res_2pup = pmapreduce(((x, y), (a, b)) -> x*a + y*b, +, iterators_product_putil, iterators_product_putil)
                @test res_2pup == res_exp_2pup
                @test res_2pup == res_2itp
            end

            @testsetwithinfo "pmapreduce_productsplit" begin
                res_exp = sum(workers())
                @test pmapreduce_productsplit(x -> myid(), +, 1:nworkers()) == res_exp
                @test pmapreduce_productsplit(NoSplat(x -> myid()), +, 1:nworkers()) == res_exp
                @test pmapreduce_productsplit(x -> myid(), +, 1:nworkers(), 1:1) == res_exp
            end
        end;

        @testsetwithinfo "inplace assignment" begin
            res = pmapreduce_productsplit(x -> ones(2), ParallelUtilities.elementwisesum!, 1:10)
            resexp = mapreduce(x -> ones(2), +, 1:min(10, nworkers()))
            @test res == resexp

            res = pmapreduce_productsplit(x -> ones(2), ParallelUtilities.elementwiseproduct!, 1:4)
            resexp = mapreduce(x -> ones(2), (x,y) -> x .* y, 1:min(4, nworkers()))
            @test res == resexp

            res = pmapreduce_productsplit(x -> ones(2), ParallelUtilities.elementwisemin!, 1:4)
            resexp = mapreduce(x -> ones(2), (x,y) -> min.(x,y), 1:min(4, nworkers()))
            @test res == resexp

            res = pmapreduce_productsplit(x -> ones(2), ParallelUtilities.elementwisemax!, 1:4)
            resexp = mapreduce(x -> ones(2), (x,y) -> max.(x,y), 1:min(4, nworkers()))
            @test res == resexp
        end

        @testsetwithinfo "concatenation" begin
            @testsetwithinfo "comparison with mapreduce" begin
                resexp_vcat = mapreduce(identity, vcat, 1:nworkers())
                resexp_hcat = mapreduce(identity, hcat, 1:nworkers())
                res_vcat = pmapreduce(identity, vcat, 1:nworkers())
                res_hcat = pmapreduce(identity, hcat, 1:nworkers())
                @test res_vcat == resexp_vcat
                @test res_hcat == resexp_hcat
            end

            @testsetwithinfo "pmapreduce_productsplit" begin
                res_vcat = mapreduce(identity, vcat, ones(2) for i in 1:nworkers())
                res_hcat = mapreduce(identity, hcat, ones(2) for i in 1:nworkers())

                @test pmapreduce_productsplit(x -> ones(2), vcat, 1:nworkers()) == res_vcat
                @test pmapreduce_productsplit(x -> ones(2), hcat, 1:nworkers()) == res_hcat
            end
        end;

        @testsetwithinfo "run elsewhere" begin
            @testsetwithinfo "sum" begin
                res_exp = sum(workers())
                c = Channel(nworkers())
                tasks = Vector{Task}(undef, nworkers())
                @sync begin
                    for (ind, p) in enumerate(workers())
                        tasks[ind] = @async begin
                            try
                                res = @fetchfrom p pmapreduce_productsplit(x -> myid(), +, 1:nworkers())
                                put!(c,(ind, res, false))
                            catch
                                put!(c,(ind, 0, true))
                                rethrow()
                            end
                        end
                    end
                    for i = 1:nworkers()
                        ind, res, err = take!(c)
                        err && wait(tasks[ind])
                        @test res == res_exp
                        showworkernumber(i, nworkers())
                    end
                end
            end
            # concatenation where the rank is used in the mapping function
            # Preserves order of the iterators
            @testsetwithinfo "concatenation using rank" begin
                c = Channel(nworkers())
                tasks = Vector{Task}(undef, nworkers())
                @sync begin
                    for (ind, p) in enumerate(workers())
                        tasks[ind] = @async begin
                            try
                                res = @fetchfrom p (pmapreduce_productsplit(x -> x[1][1], vcat, 1:nworkers()) == mapreduce(identity, vcat, 1:nworkers()))
                                put!(c,(ind, res, false))
                            catch
                                put!(c,(ind, false, true))
                                rethrow()
                            end
                        end
                    end
                    for i = 1:nworkers()
                        ind, res, err = take!(c)
                        err && wait(tasks[ind])
                        @test res
                        showworkernumber(i, nworkers())
                    end
                end
            end
        end;

        @testsetwithinfo "errors" begin
            @test_throws Exception pmapreduce(x -> error("map"), +, 1:10)
            @test_throws Exception pmapreduce(identity, x -> error("reduce"), 1:10)
            @test_throws Exception pmapreduce(x -> error("map"), x -> error("reduce"), 1:10)

            @test_throws Exception pmapreduce(fmap, +, 1:10)
            @test_throws Exception pmapreduce(identity, fred, 1:10)
            @test_throws Exception pmapreduce(fmap, fred, 1:10)

            if nworkers() != nprocs()
                @test_throws Exception pmapreduce(fmap_local, +, 1:10)
                @test_throws Exception pmapreduce(identity, fred_local, 1:10)
                @test_throws Exception pmapreduce(fmap_local, fred, 1:10)
                @test_throws Exception pmapreduce(fmap_local, fred_local, 1:10)
            end
        end;
    end;
    @testsetwithinfo "pmapbatch" begin
        for (iterators, fmap) in Any[
            ((1:1,), x -> 1),
            ((1:10,), x -> 1),
            ((1:5,), x -> ones(1) * x),
            ((1:10, 1:10), (x,y) -> ones(3) * (x+y))]

            res = pmapbatch(fmap, iterators...)
            res_exp = pmap(fmap, iterators...)
            @test res == res_exp
        end

        v = pmapbatch_productsplit(x -> sum(sum(i) for i in x) * ones(2), 1:1, 1:1)
        @test v == [[2.0, 2.0]]
        v = pmapbatch_productsplit(x -> ParallelUtilities.workerrank(x), 1:nworkers(), 1:nworkers())
        @test v == [1:nworkers();]
    end
end;

@testset "ClusterQueryUtils" begin
    # These tests assume that all the workers are on the same node
    p = procs_node()
    myhost = Libc.gethostname()
    if myhost in keys(p)
        w_myhost = p[myhost]
        @test sort(workers_myhost(w_myhost)) == sort(w_myhost)
        @test sort(workers_myhost(WorkerPool(w_myhost))) == sort(w_myhost)

        w = oneworkerpernode(w_myhost)
        @test length(w) == 1
        @test (@fetchfrom w[1] Libc.gethostname()) == myhost

        pool = workerpool_nodes(WorkerPool(w_myhost))
        w_pool = workers(pool)
        @test length(w_pool) == 1
        @test (@fetchfrom w_pool[1] Libc.gethostname()) == myhost
    end

    w = workers(workerpool_nodes())
    @test !isempty(w)
    host_w = @fetchfrom w[1] Libc.gethostname()
    @test w[1] in p[host_w]
end


text:
struct ObsData{T, K, V} <: AbstractIdData
    obs::T
    var::Symbol
    id::Dict{K, V}
end

struct Descriptives{T <: ObsData} <: AbstractIDResult{T}
    data::T
    result::Dict
end

################################################################################


text:
using CairoMakie
using Clustering
using Distances
using Distributions
using KernelFunctions
using LinearAlgebra
using PrototypesCriticisms
using Random

"""
    output(D, protoids, critids, axis, title; color=nothing)

Print basic text output, and produce plot for the data, prototypes, and criticisms.
"""
function output(D, protoids, critids, title; color=nothing)
    println("$title:")
    println("  Prototypes: ", protoids)
    println("  Criticisms: ", critids)
    protos = D[:, vcat(protoids...)]
    crits = D[:, vcat(critids...)]
    fig = Figure()
    Axis(fig[1, 1], xgridvisible=false, ygridvisible=false, limits=(0, 6, 0, 4), xlabel=L"x", ylabel=L"y", title=title)
    if isnothing(color)
        scatter!(D[1, :], D[2, :], markersize=4)
    else
        scatter!(D[1, :], D[2, :], color=color, markersize=4)
    end
    scatter!(protos[1, :], protos[2, :], color=:red, marker=:cross)
    scatter!(crits[1, :], crits[2, :], color=:green, marker=:xcross)
    save("out/$(replace(title, ' ' => '_')).pdf", fig)
end

"""
    main()

Run the example.
"""
function main()
    Random.seed!(42)
    update_theme!(font="Libertinus Serif")
    mkpath("out")

    n = 40
    D1 = [rand(Normal(1, 0.1), n) rand(Normal(1, 0.3), n)]
    D2 = [rand(Normal(3, 0.3), n) rand(Normal(3, 0.3), n)]
    D3 = [rand(Normal(4, 0.5), n) rand(Normal(2, 0.3), n)]
    D4 = [rand(Normal(1.5, 0.1), n) rand(Normal(3, 0.1), n)]
    D = [D1; D2; D3; D4]'

    k = 4
    p = 1
    c = 1

    # k-medoids
    clustering = kmedoids(pairwise(Euclidean(), D), k)
    protoids = prototypes(clustering, p)
    critids = criticisms(clustering, c)
    output(D, protoids, critids, "k-medoids", color=assignments(clustering))

    # k-means
    clustering = kmeans(D, k)
    protoids = prototypes(clustering, p)
    critids = criticisms(clustering, c)
    output(D, protoids, critids, "k-means", color=assignments(clustering))

    # fuzzy c-means
    clustering = fuzzy_cmeans(D, k, 2)
    protoids = prototypes(clustering, p)
    critids = criticisms(clustering, c)
    output(D, protoids, critids, "fuzzy c-means", color=map(i -> i[2], argmax(clustering.weights, dims=2)))

    # affinity propagation
    S = -pairwise(Euclidean(), D, dims=2)
    S = S - diagm(0 => diag(S)) + median(S) * I
    clustering = affinityprop(S)
    protoids = prototypes(clustering, D, p)
    critids = criticisms(clustering, D, c)
    output(D, protoids, critids, "affinity propagation", color=assignments(clustering))

    # MMD-critic
    kernel = with_lengthscale(RBFKernel(), sqrt(size(D, 1)))
    protoids = prototypes(D, kernel, k)
    critids = criticisms(D, kernel, protoids, k)
    output(D, protoids, critids, "MMD-critic", color=:grey)
end

if abspath(PROGRAM_FILE) == @__FILE__
    main()
end