text:
#######################
# static DSL function #
#######################

using FunctionalCollections: PersistentVector

struct Params
    prob_outlier::Float64
    slope::Float64
    intercept::Float64
    inlier_std::Float64
    outlier_std::Float64
end

@gen (static) function datum(x::Float64, (grad)(params::Params))
    is_outlier = @trace(bernoulli(params.prob_outlier), :z)
    std = is_outlier ? params.inlier_std : params.outlier_std
    y = @trace(normal(x * params.slope + params.intercept, std), :y)
    return y
end

data_fn = Map(datum)

"""
my documentation
"""
@gen (static) function model(xs::Vector{Float64})
    n = length(xs)
    inlier_std = @trace(gamma(1, 1), :inlier_std)
    outlier_std = @trace(gamma(1, 1), :outlier_std)
    slope = @trace(normal(0, 2), :slope)
    intercept = @trace(normal(0, 2), :intercept)
    params::Params = Params(0.5, inlier_std, outlier_std, slope, intercept)
    ys = @trace(data_fn(xs, fill(params, n)), :data)
    return ys
end

@gen (static) function at_choice_example_1(i::Int)
    ret = @trace(bernoulli(0.5), :x => i)
end

# @trace(choice_at(bernoulli)(0.5, i), :x)

@gen (static) function at_choice_example_2(i::Int)
    ret = @trace(bernoulli(0.5), :x => i => :y)
end

# @trace(call_at(choice_at(bernoulli))(0.5, i, :y), :x)

@gen function foo(mu)
    @trace(normal(mu, 1), :y)
end

@gen (static) function at_call_example_1(i::Int)
    mu = 1.123
    ret = @trace(foo(mu), :x => i)
end

@gen (static) function at_call_example_2(i::Int)
    mu = 1.123
    ret = @trace(foo(mu), :x => i => :y)
end

@testset "static DSL" begin

function get_node_by_name(ir, name::Symbol)
    nodes = filter((node) -> (node.name == name), ir.nodes)
    @assert length(nodes) == 1
    nodes[1]
end

@testset "check IR" begin

#####################
# check IR of datum #
#####################

ir = Gen.get_ir(typeof(datum))

# argument nodes
@test length(ir.arg_nodes) == 2
x = ir.arg_nodes[1]
params = ir.arg_nodes[2]
@test x.name == :x
@test x.typ == :Float64
@test !x.compute_grad
@test params.name == :params
@test params.typ == :Params
@test params.compute_grad

# choice nodes and call nodes
@test length(ir.choice_nodes) == 2
@test length(ir.call_nodes) == 0

# is_outlier
is_outlier = ir.choice_nodes[1]
@test is_outlier.name == :is_outlier
@test is_outlier.addr == :z
@test is_outlier.typ == QuoteNode(Bool)
@test is_outlier.dist == bernoulli
@test length(is_outlier.inputs) == 1

# std
std = get_node_by_name(ir, :std)
@test isa(std, Gen.JuliaNode)
@test std.name == :std
@test std.typ == QuoteNode(Any)
@test length(std.inputs) == 2
in1 = std.inputs[1]
in2 = std.inputs[2]
@test (in1 === is_outlier && in2 === params) || (in2 === is_outlier && in1 === params)

# y
y = ir.choice_nodes[2]
@test y.name == :y
@test y.addr == :y
@test y.typ == QuoteNode(Float64)
@test y.dist == normal
@test length(y.inputs) == 2
@test y.inputs[2] === std

# y_mean
y_mean = y.inputs[1]
@test isa(y_mean, Gen.JuliaNode)
@test y_mean.typ == QuoteNode(Any)
@test length(y_mean.inputs) == 2
in1 = y_mean.inputs[1]
in2 = y_mean.inputs[2]
@test (in1 === x && in2 === params) || (in2 === x && in1 === params)

# prob outlier
prob_outlier = is_outlier.inputs[1]
@test isa(prob_outlier, Gen.JuliaNode)
@test length(prob_outlier.inputs) == 1
@test prob_outlier.inputs[1] === params
@test prob_outlier.typ == QuoteNode(Any)

@test ir.return_node === y

#####################
# check IR of model #
#####################

ir = Gen.get_ir(typeof(model))
@test length(ir.arg_nodes) == 1
xs = ir.arg_nodes[1]
@test xs.name == :xs
@test xs.typ == :(Vector{Float64})
@test !xs.compute_grad

# choice nodes and call nodes
@test length(ir.choice_nodes) == 4
@test length(ir.call_nodes) == 1

# inlier_std
inlier_std = ir.choice_nodes[1]
@test inlier_std.name == :inlier_std
@test inlier_std.addr == :inlier_std
@test inlier_std.typ == QuoteNode(Float64)
@test inlier_std.dist == gamma
@test length(inlier_std.inputs) == 2

# outlier_std
outlier_std = ir.choice_nodes[2]
@test outlier_std.name == :outlier_std
@test outlier_std.addr == :outlier_std
@test outlier_std.typ == QuoteNode(Float64)
@test outlier_std.dist == gamma
@test length(outlier_std.inputs) == 2

# slope
slope = ir.choice_nodes[3]
@test slope.name == :slope
@test slope.addr == :slope
@test slope.typ == QuoteNode(Float64)
@test slope.dist == normal
@test length(slope.inputs) == 2

# intercept
intercept = ir.choice_nodes[4]
@test intercept.name == :intercept
@test intercept.addr == :intercept
@test intercept.typ == QuoteNode(Float64)
@test intercept.dist == normal
@test length(intercept.inputs) == 2

# data
ys = ir.call_nodes[1]
@test ys.name == :ys
@test ys.addr == :data
@test ys.typ == QuoteNode(PersistentVector{Float64})
@test ys.generative_function == data_fn
@test length(ys.inputs) == 2
@test ys.inputs[1] == xs

# params
params = get_node_by_name(ir, :params)
@test isa(params, Gen.JuliaNode)
@test params.name == :params
@test params.typ == :Params
@test length(params.inputs) == 4
@test slope in params.inputs
@test intercept in params.inputs
@test inlier_std in params.inputs
@test outlier_std in params.inputs

# n
n = get_node_by_name(ir, :n)
@test isa(n, Gen.JuliaNode)
@test n.name == :n
@test n.typ == QuoteNode(Any)
@test length(n.inputs) == 1
@test n.inputs[1] === xs

# params_vec
params_vec = ys.inputs[2]
@test isa(params_vec, Gen.JuliaNode)
@test params_vec.typ == QuoteNode(Any)
@test length(params_vec.inputs) == 2
in1 = params_vec.inputs[1]
in2 = params_vec.inputs[2]
@test (in1 === params && in2 === n) || (in2 === params && in1 === n)

@test ir.return_node === ys

end

@testset "at_choice" begin

# at_choice_example_1
ir = Gen.get_ir(typeof(at_choice_example_1))
i = get_node_by_name(ir, :i)
ret = get_node_by_name(ir, :ret)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 2
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test ret.inputs[2] === i
at = ret.generative_function
@test isa(at, Gen.ChoiceAtCombinator)
@test at.dist == bernoulli

# at_choice_example_2
ir = Gen.get_ir(typeof(at_choice_example_2))
i = get_node_by_name(ir, :i)
ret = get_node_by_name(ir, :ret)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 3
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test isa(ret.inputs[2], Gen.JuliaNode) # () -> :y
@test ret.inputs[3] === i
at = ret.generative_function
@test isa(at, Gen.CallAtCombinator)
at2 = at.kernel
@test isa(at2, Gen.ChoiceAtCombinator)
@test at2.dist == bernoulli

end

@testset "at_call" begin

# at_call_example_1
ir = Gen.get_ir(typeof(at_call_example_1))
i = get_node_by_name(ir, :i)
ret = get_node_by_name(ir, :ret)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 2
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test ret.inputs[2] === i
at = ret.generative_function
@test isa(at, Gen.CallAtCombinator)
@test at.kernel == foo

#at_call_example_2
ir = Gen.get_ir(typeof(at_call_example_2))
i = get_node_by_name(ir, :i)
ret = get_node_by_name(ir, :ret)
@test isa(ret, Gen.GenerativeFunctionCallNode)
@test ret.addr == :x
@test length(ret.inputs) == 3
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test isa(ret.inputs[1], Gen.JuliaNode) # () -> 0.5
@test isa(ret.inputs[2], Gen.JuliaNode) # () -> :y
@test ret.inputs[3] === i
at = ret.generative_function
@test isa(at, Gen.CallAtCombinator)
at2 = at.kernel
@test isa(at2, Gen.CallAtCombinator)
@test at2.kernel == foo

end

@testset "tracked diff annotation" begin

@gen (static, diffs) function bar(x)
    return x
end

@test Gen.get_options(typeof(bar)).track_diffs

@gen (static) function bar(x)
    return x
end

@test !Gen.get_options(typeof(bar)).track_diffs

end

@testset "no julia cache annotation" begin

@gen (static, nojuliacache) function bar(x)
    return x
end

@test !Gen.get_options(typeof(bar)).cache_julia_nodes

@gen (static) function bar(x)
    return x
end

@test Gen.get_options(typeof(bar)).cache_julia_nodes

end

@testset "trainable parameters" begin

@gen (static) function foo()
    @param theta::Float64
    return theta
end

ir = Gen.get_ir(typeof(foo))
theta = get_node_by_name(ir, :theta)
@test isa(theta, Gen.TrainableParameterNode)
@test ir.return_node === theta

end

@testset "use of 'end'" begin

@gen (static) function foo()
    x = [1, 2, 3, 4, 5, 6]
    y = x[3:end]
    return y
end

load_generated_functions()

@test foo() == [3, 4, 5, 6]

end

@testset "docstrings" begin
    io = IOBuffer()
    print(io, @doc model)
    @test String(take!(io)) == "my documentation\n"
end

end # @testset "static DSL"

max_stars_repo_name: dendisuhubdy/Gen

max_issues_repo_name: kumsh/Gen

max_forks_repo_name: kumsh/Gen


text:
@system LeafInitiation(Stage, Planting, FloralInitiation) begin
    SD: storage_days => 100 ~ preserve(u"d", parameter)
    ST: storage_temperature => 5 ~ preserve(u"°C", parameter)
    #FIXME: ThermalTime only accepts override of track, not preserve
    STP(ST): storage_tempreature_proxy ~ track(u"°C")
    SBF(context, T=STP, To=T_opt', Tx=T_ceil'): storage_beta_function ~ ::BetaFunction
    ILN: initial_leaves_at_harvest => 4 ~ preserve::int(parameter)
    ILS(LIR_max, β=SBF.ΔT, SD): initial_leaves_during_storage => begin
        LIR_max * β * SD
    end ~ preserve::int(round=:floor)
    initial_leaves(ILN, ILS) => ILN + ILS ~ preserve::int

    LIR_max: maximum_leaf_initiation_rate => 0.20 ~ preserve(u"d^-1", parameter)

    leaf_initiation(r=LIR_max, β=BF.ΔT) => r*β ~ accumulate(when=leaf_initiating)

    #HACK original Radish model assumed leaves are being initiated when the seeds are sown
    #HACK maize model model assumed leaf initiation begins when germination is over
    leaf_initiateable(planted) ~ flag

    # for maize
    #leaf_initiated(pheno.tassel_initiation.over) ~ flag
    # for Radish
    leaf_initiated(floral_initiated) ~ flag

    leaf_initiating(leaf_initiateable & !leaf_initiated) ~ flag

    # no MAX_LEAF_NO implied unlike original model
    leaves_initiated(initial_leaves, leaf_initiation) => begin
        initial_leaves + leaf_initiation
    end ~ track::int(round=:floor)
end

ext: jl


text:
#####
##### CPU box pressure solver
#####

function BoxPressureSolver(::CPU, grid, pressure_bcs, planner_flag=FFTW.PATIENT)
    kx², ky², kz² = generate_discrete_eigenvalues(grid, pressure_bcs)
    wavenumbers = (kx² = kx², ky² = ky², kz² = kz²)

    # Storage for RHS and Fourier coefficients is hard-coded to be Float64
    # because of precision issues with Float32.
    # See https://github.com/climate-machine/Oceananigans.jl/issues/55
    storage = zeros(Complex{Float64}, grid.Nx, grid.Ny, grid.Nz)

    @debug "Planning transforms for PressureSolver{Box, CPU}..."
    x_bc, y_bc, z_bc = pressure_bcs.x.left, pressure_bcs.y.left, pressure_bcs.z.left
    DCTxyz!  =  plan_forward_transform(storage, z_bc, [1, 2, 3], planner_flag)
    IDCTxyz! = plan_backward_transform(storage, z_bc, [1, 2, 3], planner_flag)
    @debug "Planning transforms for PressureSolver{Box, CPU} done!"

    transforms = (DCTxyz! =  DCTxyz!, IDCTxyz! = IDCTxyz!)

    return PressureSolver(Box(), CPU(), wavenumbers, storage, transforms, nothing)
end

function solve_poisson_equation!(solver::PressureSolver{Box, CPU}, grid)
    Nx, Ny, Nz, _ = unpack_grid(grid)
    kx², ky², kz² = solver.wavenumbers.kx², solver.wavenumbers.ky², solver.wavenumbers.kz²

    # We can use the same storage for the RHS and the solution ϕ.
    RHS, ϕ = solver.storage, solver.storage

    solver.transforms.DCTxyz! * RHS  # Calculate DCTˣʸᶻ(f) in place.

    # Solve the discrete Poisson equation in spectral space. We are computing
    # the Fourier coefficients of the solution from the Fourier coefficients
    # of the RHS.
    @. ϕ = -RHS / (kx² + ky² + kz²)

    # Setting DC component of the solution (the mean) to be zero. This is also
    # necessary because the source term to the Poisson equation has zero mean
    # and so the DC component comes out to be ∞.
    ϕ[1, 1, 1] = 0

    solver.transforms.IDCTxyz! * ϕ  # Calculate IDCTˣʸᶻ(ϕ) in place.

    # Must normalize by 2N for each dimension transformed via FFTW.REDFT.
    @. ϕ = ϕ / (8Nx*Ny*Nz)

    return nothing
end

#####
##### GPU channel pressure solver
#####

function BoxPressureSolver(::GPU, grid, pressure_bcs, no_args...)
    throw(ArgumentError("Box pressure solver not implemented for GPUs yet :("))
end


text:
# ---
# jupyter:
#   jupytext:
#     formats: ipynb,jl:hydrogen
#     text_representation:
#       extension: .jl
#       format_name: hydrogen
#       format_version: '1.3'
#       jupytext_version: 1.11.2
#   kernelspec:
#     display_name: Julia 1.6.2
#     language: julia
#     name: julia-1.6
# ---

# %% [markdown]
# https://discourse.julialang.org/t/problem-with-distrubuted-package/65199/5

# %%
using Distributed
using DifferentialEquations
# using Plots

@show addprocs()

@everywhere using DifferentialEquations
# Linear ODE which starts at 0.5 and solves from t=0.0 to t=1.0
@everywhere f(u,p,t) = 1.01*u
prob = ODEProblem(f,0.5,(0.0,1.0))

@everywhere function prob_func(prob,i,repeat) # i is the unique id 1:trajectories for the problem
  remake(prob,u0=rand()*prob.u0)
end

ensemble_prob = EnsembleProblem(prob,prob_func=prob_func)
sim = solve(ensemble_prob,Tsit5(),EnsembleDistributed(),trajectories=100)

# %%
versioninfo()

# %%
using Pkg
Pkg.status("DifferentialEquations")
Pkg.status("DiffEqBase"; mode=PKGMODE_MANIFEST)
Pkg.status("SciMLBase"; mode=PKGMODE_MANIFEST)

# %%

max_stars_repo_path: 0013/EnsembleProblem example.jl

max_stars_repo_name: genkuroki/public

max_issues_repo_path: 0013/EnsembleProblem example.jl

max_issues_repo_name: genkuroki/public

max_forks_repo_path: 0013/EnsembleProblem example.jl

max_forks_repo_name: genkuroki/public


text:
"""
    pve(y, X, β; l = IdentityLink())

Estimates phenotype's Proportion of Variance Explained (PVE) by typed genotypes
(i.e. chip heritability or SNP heritability).

# Model
We compute `Var(ŷ) / Var(y)` where `y` is the raw phenotypes, `X` contains
all the genotypes, and `ŷ = Xβ` is the predicted (average) phenotype
values from the statistical model β. Intercept is NOT included.
"""
function pve(
    y::AbstractVecOrMat,
    X::AbstractVecOrMat,
    β::AbstractVecOrMat;
    l::Link = IdentityLink()
    )
    μ = linkinv.(l, X * β) # mean
    return _pve(y, μ)
end

function _pve(y::AbstractVector, μ::AbstractVector)
    return var(μ) / var(y)
end

# each column is a trait
function _pve(Y::AbstractMatrix, μ::AbstractMatrix)
    size(Y, 2) == size(μ, 2) || error("Number of columns not equal in Y and μ!")
    return [_pve(@view(Y[i, :]), @view(μ[i, :])) for i in 1:size(Y, 2)]
end

function pve(v::IHTVariable)
    return _pve(v.y, v.μ)
end

function pve(v::mIHTVariable)
    return [_pve(@view(v.Y[i, :]), @view(v.μ[i, :])) for i in 1:ntraits(v)]
end


text:
module Foo

using Preferences

set!(key, value) = @set_preferences!(key=>value)
get(key) = @load_preference(key)

end # module


text:
##############################################################################
#
# VortexParticle2D.jl
#
# Part of CVortex.jl
# Representation of a 2D vortex particle.
#
# Copyright 2019-2020 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
##############################################################################

"""
	Representation of a 2D vortex particle in CVortex

You do not need this to use the CVortex API.

coord is a particle's position.
vorticity is the particle's vorticity
volume is the volume of the particle. This is only important for
viscous vortex particle strength exchange methods (not included in this
wrapper)
"""
struct VortexParticle2D
    coord :: Vec2f
    vorticity :: Float32
    volume :: Float32
end

function VortexParticle2D(coord::Vec2f, vort::Real)
    return VortexParticle2D(coord, vort, 0.0)
end

function VortexParticle2D(coord::Vector{<:Real}, vort::Real, vol::Real)
	@assert(length(coord)==2)
    return VortexParticle2D(Vec2f(coord), vort, vol)
end

function VortexParticle2D(coord::Vector{<:Real}, vort::Real)
	@assert(length(coord)==3)
    return VortexParticle2D(coord, vort, 0.0)
end

function particle_induced_velocity(
    inducing_particle_position :: Vector{<:Real},
    inducing_particle_vorticity :: Real,
	measurement_point :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real
	) :: Vector{Float32}

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec2f_vect(measurement_point, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")

    inducing_particle = VortexParticle2D(
        inducing_particle_position, inducing_particle_vorticity, 0.0)
    mes_pnt = Vec2f(measurement_point)
	ret = Vec2f(0., 0.)
	#=
    CVTX_EXPORT bsv_V2f cvtx_P2D_S2S_vel(
        const cvtx_P2D *self,
        const bsv_V2f mes_point,
        const cvtx_VortFunc *kernel,
        float regularisation_radius);
	=#
	ret = ccall(
			("cvtx_P2D_S2S_vel", libcvortex),
			Vec2f,
			(Ref{VortexParticle2D}, Vec2f, Ref{RegularisationFunction}, Cfloat),
			inducing_particle, mes_pnt, kernel, regularisation_radius
			)
	return Vector{Float32}(ret)
end

function particle_induced_velocity(
    inducing_particle_position :: Vector{<:Real},
    inducing_particle_vorticity :: Real,
	measurement_points :: Matrix{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real
	) :: Matrix{Float32}

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec2f_vect(measurement_points, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	ni = size(measurement_points)[1]
    mes_pnt = map(i->Vec2f(measurement_points[i,:]), 1:ni)

    inducing_particle = VortexParticle2D(
        inducing_particle_position, inducing_particle_vorticity, 0.0)
	ret = Vector{Vec2f}(undef, ni)
	#=
	CVTX_EXPORT void cvtx_P2D_S2M_vel(
		const cvtx_P2D* self,
		const bsv_V2f* mes_start,
		const int num_mes,
		bsv_V2f* result_array,
		const cvtx_VortFunc* kernel,
		float regularisation_radius);
	=#
	ccall(
		("cvtx_P2D_S2M_vel", libcvortex),
		Cvoid,
		(Ref{VortexParticle2D}, Ptr{Vec2f}, Cint,
		Ref{Vec2f}, Ref{RegularisationFunction}, Cfloat),
		inducing_particle, mes_pnt, ni, ret, kernel, regularisation_radius
		)
	return Matrix{Float32}(ret)
end

function particle_induced_velocity(
    inducing_particle_position :: Matrix{<:Real},
    inducing_particle_vorticity :: Vector{<:Real},
	measurement_point :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real
	) :: Vector{Float32}

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec2f_vect(measurement_point, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")

	np = size(inducing_particle_position)[1]
    inducing_particles = map(
        i->VortexParticle2D(
            inducing_particle_position[i, :],
            inducing_particle_vorticity[i], 0.0),
        1:np)
    mes_pnt = Vec2f(measurement_point)

	pargarr = Vector{Ptr{VortexParticle2D}}(undef, length(inducing_particles))
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	ret =Vec2f(0., 0.)
	GC.@preserve pargarr inducing_particles begin
		#=
		CVTX_EXPORT bsv_V2f cvtx_P2D_M2S_vel(
			const cvtx_P2D **array_start,
			const int num_particles,
			const bsv_V2f mes_point,
			const cvtx_VortFunc *kernel,
			float regularisation_radius);
		=#
		ret = ccall(
				("cvtx_P2D_M2S_vel", libcvortex),
				Vec2f,
				(Ref{Ptr{VortexParticle2D}}, Cint, Vec2f,
					Ref{RegularisationFunction}, Cfloat),
				pargarr, np, mes_pnt, kernel,	regularisation_radius)
	end
	return [ret.x, ret.y]
end

function particle_induced_velocity(
    inducing_particle_position :: Matrix{<:Real},
    inducing_particle_vorticity :: Vector{<:Real},
	measurement_points :: Matrix{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real
	) :: Matrix{Float32}

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity)
	convertable_to_Vec2f_vect(measurement_points, "measurement_points")
	convertable_to_F32(regularisation_radius, "regularisation_radius")

	np = size(inducing_particle_position)[1]
	ni = size(measurement_points)[1]
    inducing_particles = map(
        i->VortexParticle2D(
            inducing_particle_position[i, :],
            inducing_particle_vorticity[i], 0.0),
        1:np)
    mes_pnt = map(i->Vec2f(measurement_points[i,:]), 1:ni)

	pargarr = Vector{Ptr{VortexParticle2D}}(undef, np)
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	ret = Vector{Vec2f}(undef, ni)
	GC.@preserve pargarr inducing_particles begin
		#=
		CVTX_EXPORT void cvtx_P2D_M2M_vel(
			const cvtx_P2D **array_start,
			const int num_particles,
			const bsv_V2f *mes_start,
			const int num_mes,
			bsv_V2f *result_array,
			const cvtx_VortFunc *kernel,
			float regularisation_radius);
		=#
		ccall(
			("cvtx_P2D_M2M_vel", libcvortex),
			Cvoid,
			(Ptr{Ptr{VortexParticle2D}}, Cint, Ptr{Vec2f},
				Cint, Ref{Vec2f}, Ref{RegularisationFunction}, Cfloat),
			pargarr, np, mes_pnt, ni, ret, kernel, regularisation_radius)
	end
	return Matrix{Float32}(ret)
end

function particle_visc_induced_dvort(
    inducing_particle_position :: Vector{<:Real},
	inducing_particle_vorticity :: Real,
	inducing_particle_area :: Real,
    induced_particle_position :: Vector{<:Real},
	induced_particle_vorticity :: Real,
	induced_particle_area :: Real,
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real
	) :: Float32

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_area)
	check_particle_definition_2D(induced_particle_position,
		induced_particle_vorticity, induced_particle_area)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_2D!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

	np = size(inducing_particle_position)[1]
	inducing_particle = VortexParticle2D(
			inducing_particle_position,
			inducing_particle_vorticity, inducing_particle_area)
	induced_particle = VortexParticle2D(
			induced_particle_position,
			induced_particle_vorticity, induced_particle_area)

	#=
	CVTX_EXPORT float cvtx_P2D_S2S_visc_dvort(
		const cvtx_P2D * self,
		const cvtx_P2D * induced_particle,
		const cvtx_VortFunc * kernel,
		float regularisation_radius,
		float kinematic_visc);
	=#
	ret = ccall(
		("cvtx_P2D_S2S_visc_dvort", libcvortex),
		Cfloat,
		(Ref{VortexParticle2D}, Ref{VortexParticle2D},
			Ref{RegularisationFunction}, Cfloat, Cfloat),
		inducing_particle, induced_particle, kernel,
			regularisation_radius, kinematic_visc
		)
	return ret
end

function particle_visc_induced_dvort(
    inducing_particle_position :: Vector{<:Real},
	inducing_particle_vorticity :: Real,
	inducing_particle_area :: Real,
    induced_particle_position :: Matrix{<:Real},
	induced_particle_vorticity :: Vector{<:Real},
	induced_particle_area :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real
	) :: Vector{Float32}

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_area)
	check_particle_definition_2D(induced_particle_position,
		induced_particle_vorticity, induced_particle_area)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_2D!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

	ni = size(induced_particle_position)[1]
	inducing_particle = VortexParticle2D(
			inducing_particle_position,
			inducing_particle_vorticity, inducing_particle_area)
	induced_particles = map(
		i->VortexParticle2D(
			induced_particle_position[i, :],
			induced_particle_vorticity[i], induced_particle_area[i]),
		1:ni)
	indarg = Vector{Ptr{VortexParticle2D}}(undef, ni)
	for i = 1 : length(indarg)
		indarg[i] = Base.pointer(induced_particles, i)
	end
	ret = Vector{Float32}(undef, ni)

	GC.@preserve inducing_particle indarg induced_particles begin
		#=
		CVTX_EXPORT void cvtx_P2D_S2M_visc_dvort(
			const cvtx_P2D* self,
			const cvtx_P2D** induced_start,
			const int num_induced,
			float* result_array,
			const cvtx_VortFunc* kernel,
			float regularisation_radius,
			float kinematic_visc);
		=#
		ccall(
			("cvtx_P2D_S2M_visc_dvort", libcvortex),
			Cvoid,
			(Ref{VortexParticle2D}, Ptr{Ptr{VortexParticle2D}}, Cint,
				Ptr{Float32}, Ref{RegularisationFunction}, Cfloat, Cfloat),
			inducing_particle, indarg, ni, ret,
				kernel, regularisation_radius, kinematic_visc
			)
	end
	return ret
end

function particle_visc_induced_dvort(
    inducing_particle_position :: Matrix{<:Real},
	inducing_particle_vorticity :: Vector{<:Real},
	inducing_particle_area :: Vector{<:Real},
    induced_particle_position :: Vector{<:Real},
	induced_particle_vorticity :: Real,
	induced_particle_area :: Real,
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real
	) :: Float32

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_area)
	check_particle_definition_2D(induced_particle_position,
		induced_particle_vorticity, induced_particle_area)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_2D!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

	np = size(inducing_particle_position)[1]
	inducing_particles = map(
		i->VortexParticle2D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i], inducing_particle_area[i]),
		1:np)
	induced_particle = VortexParticle2D(
			induced_particle_position,
			induced_particle_vorticity, induced_particle_area)

	pargarr = Vector{Ptr{VortexParticle2D}}(undef, length(inducing_particles))
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	GC.@preserve pargarr inducing_particles begin
		#=
		CVTX_EXPORT float cvtx_P2D_M2S_visc_dvort(
			const cvtx_P2D **array_start,
			const int num_particles,
			const cvtx_P2D *induced_particle,
			const cvtx_VortFunc *kernel,
			float regularisation_radius,
			float kinematic_visc);
		=#
		ret = ccall(
			("cvtx_P2D_M2S_visc_dvort", libcvortex),
			Cfloat,
			(Ptr{Ptr{VortexParticle2D}}, Cint, Ref{VortexParticle2D},
				Ref{RegularisationFunction}, Cfloat, Cfloat),
			pargarr, length(inducing_particles), induced_particle,
				kernel, regularisation_radius, kinematic_visc
			)
	end
	return ret
end

function particle_visc_induced_dvort(
    inducing_particle_position :: Matrix{<:Real},
	inducing_particle_vorticity :: Vector{<:Real},
	inducing_particle_area :: Vector{<:Real},
    induced_particle_position :: Matrix{<:Real},
	induced_particle_vorticity :: Vector{<:Real},
	induced_particle_area :: Vector{<:Real},
	kernel :: RegularisationFunction,
	regularisation_radius :: Real,
	kinematic_visc :: Real
	) :: Vector{Float32}

	check_particle_definition_2D(inducing_particle_position,
		inducing_particle_vorticity, inducing_particle_area)
	check_particle_definition_2D(induced_particle_position,
		induced_particle_vorticity, induced_particle_area)
	convertable_to_F32(regularisation_radius, "regularisation_radius")
	@assert(kernel.eta_2D!=Cvoid, "You cannot use this regularisation "*
		"function for viscous simulations. Consider Winckelmans or Gaussian.")

	np = size(inducing_particle_position)[1]
	ni = size(induced_particle_position)[1]
	inducing_particles = map(
		i->VortexParticle2D(
			inducing_particle_position[i, :],
			inducing_particle_vorticity[i], inducing_particle_area[i]),
		1:np)
	induced_particles = map(
		i->VortexParticle2D(
			induced_particle_position[i, :],
			induced_particle_vorticity[i], induced_particle_area[i]),
		1:ni)

	pargarr = Vector{Ptr{VortexParticle2D}}(undef, length(inducing_particles))
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	indarg = Vector{Ptr{VortexParticle2D}}(undef, ni)
	for i = 1 : length(indarg)
		indarg[i] = Base.pointer(induced_particles, i)
	end
	ret = Vector{Float32}(undef, ni)
	GC.@preserve pargarr inducing_particles indarg induced_particles begin
		#=
		CVTX_EXPORT void cvtx_P2D_M2M_visc_dvort(
			const cvtx_P2D **array_start,
			const int num_particles,
			const cvtx_P2D **induced_start,
			const int num_induced,
			float *result_array,
			const cvtx_VortFunc *kernel,
			float regularisation_radius,
			float kinematic_visc)
		=#
		ccall(
			("cvtx_P2D_M2M_visc_dvort", libcvortex),
			Cvoid,
			(Ptr{Ptr{VortexParticle2D}}, Cint, Ptr{Ptr{VortexParticle2D}}, Cint,
				Ptr{Float32}, Ref{RegularisationFunction}, Cfloat, Cfloat),
			pargarr, length(inducing_particles), indarg, length(induced_particles),
				ret, kernel, regularisation_radius, kinematic_visc
			)
	end
	return ret
end

function redistribute_particles_on_grid(
    particle_positions :: Matrix{<:Real},
	particle_vorticities :: Vector{<:Real},
	redistribution_function :: RedistributionFunction,
	grid_density :: Real;
	negligible_vort::Real=1e-4,
	max_new_particles::Integer=-1)

	@assert(0 <= negligible_vort < 1, "The negligible_vort vort must be"*
		" in [0, 1). Given "*string(negligible_vort)*"." )
	@assert((max_new_particles==-1)||(0<=max_new_particles),
		"max_new_particles must be -1 (indicating any number of new particles"*
		") or a positive integer.")
	@assert(0<grid_density, "Grid density must be positive. Was "*
		string(grid_density)*".")
	check_particle_definition_2D(particle_positions,
		particle_vorticities)
	convertable_to_F32(grid_density, "grid_density")
	convertable_to_F32(negligible_vort, "negligible_vort")

	np = size(particle_positions)[1]
	inducing_particles = map(
		i->VortexParticle2D(
			particle_positions[i, :],
			particle_vorticities[i], 0.0),
		1:np)

	pargarr = Vector{Ptr{VortexParticle2D}}(undef, np)
	for i = 1 : length(pargarr)
		pargarr[i] = Base.pointer(inducing_particles, i)
	end
	GC.@preserve pargarr inducing_particles begin
		#=
		CVTX_EXPORT int cvtx_P2D_redistribute_on_grid(
			const cvtx_P2D **input_array_start,
			const int n_input_particles,
			cvtx_P2D *output_particles,		/* input is &(*cvtx_P2D) to write to */
			int max_output_particles,		/* Set to resultant num particles.   */
			const cvtx_RedistFunc *redistributor,
			float grid_density,
			float negligible_vort)
		=#
		if max_new_particles == -1
			max_new_particles = ccall(
				("cvtx_P2D_redistribute_on_grid", libcvortex),
				Cint,
				(Ptr{Ptr{VortexParticle2D}}, Cint, Ptr{VortexParticle2D},
					Cint, Ref{RedistributionFunction}, Cfloat, Cfloat),
				pargarr, np, C_NULL, 1, redistribution_function,
				grid_density, negligible_vort)
		end

		ret = Vector{VortexParticle2D}(undef, max_new_particles)
		nnp = ccall(
			("cvtx_P2D_redistribute_on_grid", libcvortex),
			Cint,
			(Ptr{Ptr{VortexParticle2D}}, Cint, Ptr{VortexParticle2D},
				Cint, Ref{RedistributionFunction}, Cfloat, Cfloat),
			pargarr, np, ret, max_new_particles, redistribution_function,
			grid_density, negligible_vort)
	end

	# nnp is the number of new particles.
	nvorts = zeros(Float32, nnp)
	nposns = zeros(Float32, nnp, 2)
	nareas = zeros(Float32, nnp)
	for i = 1 : nnp
		nvorts[i] = ret[i].vorticity
		nposns[i, :] = Vector{Float32}(ret[i].coord)
		nareas[i] = ret[i].volume
	end
	return nposns, nvorts, nareas
end


text:
using Test
using ArrayInterface
using ResizingTools
using ResizingTools: MethodUndefineError, DimBoundsError, tailn
using LoopVectorization
using FunctionIndices
using Static
using Aqua

const V = rand(3)
const M = rand(3, 3)
const T = rand(3, 3, 3)

_params(T::UnionAll) = _params(T.body)
_params(T::DataType) = T.parameters

@testset "ResizingTools" begin
    @testset "QA" begin
        @testset "Ambiguity" begin
            ambiguities = Test.detect_ambiguities(ResizingTools)
            filter!(ambiguities) do (m1, m2)
                p1 = _params(m1.sig)
                p2 = _params(m2.sig)
                for (t1, t2) in zip(p1, p2)
                    typeintersect(t1, t2) === Union{} && return false
                end
                return true
            end
            if !isempty(ambiguities)
                for (i, amb) in enumerate(ambiguities)
                    println("Ambiguity #$i:")
                    println("  ", amb[1])
                    println("  ", amb[2])
                end
            end
            @test isempty(ambiguities)
        end
        Aqua.test_all(ResizingTools; ambiguities=false, deps_compat=false)
    end

    @testset "arraymath" begin
        include("arraymath.jl")
    end
    @testset "resize" begin
        include("resize.jl")
    end
    @testset "MISC" begin
        include("misc.jl")
    end
end


text:
# ---
# title: 703. Kth Largest Element in a Stream
# id: problem703
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Heap, Design
# link: <https://leetcode.com/problems/kth-largest-element-in-a-stream/description/>
# hidden: true
# ---
#
# Design a class to find the `kth` largest element in a stream. Note that it is
# the `kth` largest element in the sorted order, not the `kth` distinct element.
#
# Implement `KthLargest` class:
#
#   * `KthLargest(int k, int[] nums)` Initializes the object with the integer `k` and the stream of integers `nums`.
#   * `int add(int val)` Returns the element representing the `kth` largest element in the stream.
#
#
#
# **Example 1:**
#
#
#
#     **Input**
#     ["KthLargest", "add", "add", "add", "add", "add"]
#     [[3, [4, 5, 8, 2]], [3], [5], [10], [9], [4]]
#     **Output**
#     [null, 4, 5, 5, 8, 8]
#
#     **Explanation**
#     KthLargest kthLargest = new KthLargest(3, [4, 5, 8, 2]);
#     kthLargest.add(3);   // return 4
#     kthLargest.add(5);   // return 5
#     kthLargest.add(10);  // return 5
#     kthLargest.add(9);   // return 8
#     kthLargest.add(4);   // return 8
#
#
#
#
# **Constraints:**
#
#   * `1 <= k <= 104`
#   * `0 <= nums.length <= 104`
#   * `-104 <= nums[i] <= 104`
#   * `-104 <= val <= 104`
#   * At most `104` calls will be made to `add`.
#   * It is guaranteed that there will be at least `k` elements in the array when you search for the `kth` element.
#
#
## @lc code=start
using LeetCode

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


text:
module MangalTestAttributes
    using Mangal
    using Test
    include("setup.jl")

    @test typeof(attributes()) <: Vector{MangalAttribute}

    n1 = attribute(6)
    @test typeof(n1) <: MangalAttribute

end


text:
__precompile__()

module PhysConsts

using Compat
export DATA

macro __DATA__()
    joinpath(dirname(@__DIR__), "data")
end

mkpath(@__DATA__)
DATA = Dict{String, Any}()

include("DataFile.jl")

# General Type Hierachy
include("Core.jl")
include("NIST.jl")
include("Defined.jl")

# Load Data
DATAFILES = Dict{String, DataFile}()
DATAFILES["NIST"] = DataFile(
    "https://nist.gov/srd/srd_data//srd121_allascii_2014.json",
    @__DATA__,
    "srd121_allascii_2014.json"
)
obtain(DATAFILES["NIST"])

# Constants Alias
include("Universal.jl")
include("Pauli.jl")

end # module


text:
"""
    Parameters

Type containing the parameters of a magnetism simulation.

Attributes
* `L`         : The side length of the lattice, e.g. for a square lattice there will be
    `L^2` spins
* `β`         : `1 / kT`, i.e. the inverse temperature
* `n_sweeps`  : The number of update sweeps to run over the lattice
* `cutoff`    : The first percentage of observable measurements to discard, useful to get
    rid of values during thermalization period
* `start_type`: `:cold` starts the lattice with all spins aligned, `:hot` all spins random
* `seed`      : Seed value for the random number generator
"""
struct Parameters
    L::Int
    β::Float64
    n_sweeps::Int
    cutoff::Float64
    start_type::Symbol
    seed::Int
    q::Int
end

function show(io::IO, params::Parameters)
    println("""
            L          = $(params.L)
            β          = $(params.β)
            n_sweeps   = $(params.n_sweeps)
            start_type = $(params.start_type)
            seed       = $(params.seed)
            q          = $(params.q)
            """)
end


text:
# arithmetic
Base.promote_rule(::Type{S}, ::Type{T})  where {T<:Sym, S<:Number}= T
Base.promote_rule(::Type{T}, ::Type{S})  where {T<:Sym, S<:Number}= T
Base.promote_rule(::Type{S}, ::Type{T})  where {T<:Sym, S<:Irrational}= T
Base.promote_rule(::Type{T}, ::Type{S})  where {T<:Sym, S<:Irrational}= T
Base.promote_rule(::Type{Sym}, ::Type{Bool}) = Sym
Base.promote_rule(::Type{Bool}, ::Type{Sym}) = Sym
Base.convert(::Type{T}, o::Py) where {T <: Sym} = T(o)
Base.convert(::Type{T}, o::Number) where {T <: Sym} = Sym(o)

Base.:(-)(x::Sym) = Sym(-unSym(x))
Base.:(+)(x::Sym, y::Sym) = Sym(+(unSym(x), unSym(y)))
Base.:(-)(x::Sym, y::Sym) = Sym(-(unSym(x), unSym(y)))
Base.:(*)(x::Sym, y::Sym) = Sym(*(unSym(x), unSym(y)))
Base.:(/)(x::Sym, y::Sym) = Sym(/(unSym(x), unSym(y)))
Base.:(^)(x::Sym, y::Int) = ↑(↓(x)^y)
Base.:(^)(x::Sym, y::Rational) = ^(promote(x, y)...)
Base.:(^)(x::Sym, y::Real) = ↑(↓(x)^y)
Base.:(^)(x::Sym, y::Sym) = Sym(^(unSym(x), unSym(y)))

# Bools
Base.:(+)(x::Sym, y::Bool) = x + Int(y)
Base.:(+)(x::Bool, y::Sym) = y + x
Base.:(-)(x::Sym, y::Bool) = x - Int(y)
Base.:(-)(x::Bool, y::Sym) = Int(x) - y
Base.:(*)(x::Sym, y::Bool) = x * Int(y)
Base.:(*)(x::Bool, y::Sym) = Int(x) * y
Base.:(/)(x::Sym, y::Bool) = x * Int(y)
Base.:(/)(x::Bool, y::Sym) = Int(x) * y
Base.:(^)(x::Sym, y::Bool) = x^Int(y)
Base.:(^)(x::Bool, y::Sym) = Int(x)^y

max_stars_repo_name: jverzani/SymPyCall

max_issues_repo_name: jverzani/SymPyCall

max_forks_repo_name: jverzani/SymPyCall


text:
# myfunctions2.jl

using Distributed

using Serialization

#const PATH_DB = "./test_invoicing.sqlite"
#const PATH_CSV = "./bank.csv"
#const PATH_JOURNAL = "test_journal.txt"
#const PATH_LEDGER = "test_ledger.txt"

# =================================
# task_0 - get orders from Sales
# =================================
function task_0(rx)
    tx = Channel(32)
    @async while true
        if isready(tx)
            start = take!(tx)
            @info("task_0 (master): $(typeof(start)) - $start - $(start == "START")")
            if start == "START"
                @info("task_0 (master) will start the process")
                orders = @fetch AppliSales.process()
                @info("task_0 (master) will put $(length(orders)) the orders on rx channel")
                put!(rx, orders)
                @info("task_0 (master) has put $(length(orders)) the orders on rx channel")
            end
        else
            @info("task_0 (master) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task_0

# =================================
# task_1 - process the orders
# =================================
function task_1(rx)
    tx = Channel(32)
    @async while true
        if isready(tx)
            orders = take!(tx)
            @info("task 1 (process the orders): $(typeof(orders))")
            if orders isa Array{AppliSales.Order, 1}
                @info("task_1 (process the orders) will process $(length(orders)) orders remotely")
                result = @fetch AppliAR.process(orders)
                @info("task_1 (process the orders) will put $(length(result)) journal entries on rx channel")
                put!(rx, result)
            end
        else
            @info("task_1 (process the orders) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task_1

# =================================
# task_2 - process journal entries
# =================================
function task_2(rx)
    tx = Channel(32)
    @async while true
        if isready(tx)
            entries = take!(tx)
            @info("task_2 (process journal entries): $(typeof(entries))")
            if entries isa Array{AppliGeneralLedger.JournalEntry,1}
                @info("task_2 (process journal entries) will process $(length(entries)) journal entries remotely")
                result = @fetch AppliGeneralLedger.process(entries)
                #@info("task_general_ledger saved $(length(result)) journal entries")
                @info("task_2 (process journal entries) saved the journal entries")
                #put!(tx, result)
            end
        else
            @info("task_2 (process journal entries) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task2_2

# =================================
# task_3 - process payments
# =================================
function task_3(rx)
    tx = Channel(32)
    @async while true
        if isready(tx)
            stms = take!(tx)
            @info("task_3 (process payments): $(typeof(stms))")
            if stms isa Array{AppliAR.BankStatement,1}
                @info("task_3 (process payments) will match unpaid invoices with bank statements")
                result = @fetch begin
                    unpaid_invoices = retrieve_unpaid_invoices()
                    AppliAR.process(unpaid_invoices, stms)
                end
                @info("task_3 (process payments) will put $(length(result)) journal entries on rx channel")
                put!(rx, result)
            end
        else
            @info("task_3 (process payments) is waiting for data")
            wait(tx)
        end
    end
    return tx
end # task_3

# =================================
# task_4 - process payments
# =================================
function task_4(rx)
    @info("task_4 - Start socket listener")
    server = listen(IPv4(0), 8000)
    @async while true
        sock = accept(server)
        @async while isopen(sock)
            #put!(rx, readline(sock, keep=true))
            put!(rx, deserialize(sock))
        end
    end
end # testk_4

# =================================
# task dispatcher
# =================================
function dispatcher()
    rx = Channel(32)

    # instantiate tasks
    tx0 = task_0(rx) # get orders from Sales
    tx1 = task_1(rx) # process the orders
    tx2 = task_2(rx) # process the journal entries
    tx3 = task_3(rx) # process the unpaid invoices
    task_4(rx) # listen to port 8000

    # definition Holy traits pattern Dispatcher in domain.jl
    dispatch(x::T) where {T} = dispatch(Dispatcher(T), x)

    dispatch(::T0, x) = put!(tx0, x)
    dispatch(::T1, x) = put!(tx1, x)
    dispatch(::T2, x) = put!(tx2, x)
    dispatch(::T3, x) = put!(tx3, x)
    # dispatch(::T4, x) = put!(tx0, x)

    @async while true
        if isready(rx)
            value = take!(rx)
            @info("Dispatcher received $(typeof(value))")
            dispatch(value)
        else
            wait(rx)
        end
    end
    return rx
end # dispatcher

using Sockets


text:
# ============================================================================================================
# Main API

using MixedModels, GLM
@reexport using MixedModels
import MixedModels: FeMat, createAL, reweight!, getθ,
                     _iscomparable, _criterion,
                     deviance, dof, dof_residual, nobs
import StatsModels: width

anova(models::Vararg{<: MixedModel};
        test::Type{<: GoodnessOfFit} = length(models) > 1 ? LRT : FTest,
        kwargs...) =
    anova(test, models...; kwargs...)

# ==================================================================================================================
# ANOVA by F test
# Linear mixed-effect models

function anova(::Type{FTest},
        model::LinearMixedModel;
        type::Int = 1,
        adjust_sigma::Bool = true,
        kwargs...)

    type == 2           && throw(ArgumentError("Type 2 anova is not implemented"))
    type in [1, 2, 3]   || throw(ArgumentError("Invalid type"))

    varβ = vcov(model)
    β = fixef(model)

    assign = asgn(first(model.formula.rhs))

    # calculate degree of freedom for factors and residuals
    df, resdf = calcdof(model)

    # use MMatrix/SizedMatrix ?
    if type == 1
        fs = abs2.(cholesky(Hermitian(inv(varβ))).U  * β)
        # adjust σ like linear regression
        adjust = 1.0
        offset = first(assign) - 1
        model.optsum.REML || adjust_sigma && (adjust = (nobs(model) - length(β)) / nobs(model))
        fstat = ntuple(last(assign) - offset) do fix
            sum(fs[findall(==(fix + offset), assign)]) / df[fix] * adjust
        end
    else
        # calculate block by block
        adjust = 1.0
        model.optsum.REML || adjust_sigma && (adjust = (nobs(model) - length(β)) / nobs(model))
        offset = first(assign) - 1
        fstat = ntuple(last(assign) - offset) do fix
            select = findall(==(fix + offset), assign)
            β[select]' * (varβ[select, select] \ β[select]) / df[fix] * adjust
        end
    end

    pvalue = ntuple(lastindex(fstat)) do id
            ccdf(FDist(df[id], resdf[id]), abs(fstat[id]))
    end
    AnovaResult{FTest}(model, type, df, ntuple(x->NaN, length(fstat)), fstat, pvalue, (resdof = resdf,))
end

#=
function anova(::Type{FTest},
        model::GeneralizedLinearMixedModel;
        type::Int = 1,
        adjust_sigma::Bool = true,
        kwargs...)
        # use refit! and deviance as in GLM?
end
=#
# ==================================================================================================================
# ANOVA by Likehood-ratio test
# Linear mixed-effect models

function anova(::Type{LRT}, model::LinearMixedModel; kwargs...)
    # check if fitted by ML
    # nested random effects for REML ?
    model.optsum.REML && throw(
        ArgumentError("""Likelihood-ratio tests for REML-fitted models are only valid when the fixed-effects specifications are identical"""))
    @warn "Fit all submodels"
    models = nestedmodels(model; null = isnullable(model))
    anova(LRT, models...; check = false, isnested = true)
end

# =================================================================================================================
# Nested models

function anova(::Type{LikelihoodRatioTest},
                models::Vararg{<: MixedModel};
                check::Bool = true,
                isnested::Bool = false,
                kwargs...)

    check && (_iscomparable(models...) || throw(
        ArgumentError("""Models are not comparable: are the objectives, data and, where appropriate, the link and family the same?""")))
    # _isnested (by QR) or isnested (by formula) are not part of _iscomparable:
    # isnested = true
    df = dof.(models)
    ord = sortperm(collect(df))
    df = df[ord]
    models = models[ord]
    _lrt_nested(models, df, deviance.(models), 1; nestedwarn = isnested)
end

#=
function anova(::Type{LikelihoodRatioTest},
                models::Vararg{<: GeneralizedLinearMixedModel};
                check::Bool = true,
                isnested::Bool = false,
                kwargs...)
    # need new nestedmodels
end
=#

# Compare to GLM
function anova(m0::Union{TableRegressionModel{<: Union{LinearModel, GeneralizedLinearModel}}, LinearModel, GeneralizedLinearModel},
                m::T,
                ms::Vararg{T};
                check::Bool = true,
                isnested::Bool = false,
                kwargs...) where {T <: MixedModel}
    # Contain _isnested (by QR) and test on formula
    check && (_iscomparable(m0, m) || throw(
        ArgumentError("""Models are not comparable: are the objectives, data and, where appropriate, the link and family the same?""")))
    check && (_iscomparable(m, ms...) || throw(
        ArgumentError("""Models are not comparable: are the objectives, data and, where appropriate, the link and family the same?""")))
    m = [m, ms...]
    df = dof.(m)
    ord = sortperm(df)
    df = (dof(m0), df[ord]...)
    models = (m0, m[ord]...)
    # isnested is not part of _iscomparable:
    # isnested = true
    dev = (_criterion(m0), deviance.(models[2:end])...)
    _lrt_nested(models, df, dev, 1; nestedwarn = isnested)
end

# =================================================================================================================================
# Fit new models

"""
    anova_lme(f::FormulaTerm, tbl; test::Type{<: GoodnessOfFit} = FTest, <keyword arguments>)

    anova_lme(test::Type{<: GoodnessOfFit}, f::FormulaTerm, tbl; <keyword arguments>)

    anova(test::Type{<: GoodnessOfFit}, ::Type{<: LinearMixedModel}, f::FormulaTerm, tbl;
            type::Int = 1,
            adjust_sigma::Bool = true, <keyword arguments>)

ANOVA for linear mixed-effect models.

The arguments `f` and `tbl` are `Formula` and `DataFrame`.
* `test`: `GoodnessOfFit`. The default is `FTest`.
* `type`: type of anova. Only 1, 3 are valid.
* `adjust_sigma` determines whether adjust σ to match that of linear mixed-effect model fitted by REML. The result will be slightly deviated from that of model fitted by REML.

Other keyword arguments
* `wts = []`
* `contrasts = Dict{Symbol,Any}()`
* `progress::Bool = true`
* `REML::Bool = true`

`anova_lme` generate a `LinearMixedModel` fitted with REML if applying [`FTest`](@ref); otherwise, a model fitted with ML.
"""
anova_lme(f::FormulaTerm, tbl;
        test::Type{<: GoodnessOfFit} = FTest,
        kwargs...) =
    anova(test, LinearMixedModel, f, tbl; kwargs...)

anova_lme(test::Type{<: GoodnessOfFit}, f::FormulaTerm, tbl;
        kwargs...) =
    anova(test, LinearMixedModel, f, tbl; kwargs...)

function anova(test::Type{<: GoodnessOfFit}, ::Type{<: LinearMixedModel}, f::FormulaTerm, tbl;
        type::Int = 1,
        adjust_sigma::Bool = true,
        wts = [],
        contrasts = Dict{Symbol,Any}(),
        progress::Bool = true,
        REML::Bool = test == FTest ? true : false)
    model = lme(f, tbl; wts, contrasts, progress, REML)
    anova(test, model; type, adjust_sigma)
end

"""
    lme(f::FormulaTerm, tbl; wts, contrasts, progress, REML)

An alias for `fit(LinearMixedModel, f, tbl; wts, contrasts, progress, REML)`.
"""
lme(f::FormulaTerm, tbl;
    wts = [],
    contrasts = Dict{Symbol, Any}(),
    progress::Bool = true,
    REML::Bool = false) =
    fit(LinearMixedModel, f, tbl;
        wts, contrasts, progress, REML)

"""
    glme(f::FormulaTerm, tbl, d::Distribution, l::Link; kwargs...)

An alias for `fit(GeneralizedLinearMixedModel, f, tbl, d, l; kwargs...)`
"""
glme(f::FormulaTerm, tbl, d::Distribution = Normal(), l::Link = canonicallink(d);
    kwargs...) =
    fit(GeneralizedLinearMixedModel, f, tbl, d, l; kwargs...)


text:
# External modules
using SparseArrays
using LinearAlgebra
using Random

######################
# NOTE: This is not a long term solution.
# Once the package is set up, we can put this dependency in the .toml
force_package_install = true;
need_restart = false;
try
    using SymEngine
catch e
    if force_package_install
        println("SymEngine package is not yet installed. Installing now.");
        using Pkg
        Pkg.add("SymEngine")
        using SymEngine

    else
        println("SymEngine package is not yet installed. It is required.\nWould you like to install now? y/n.");
        response = readline();
        if response=="Y" || response=="y"
            using Pkg
            Pkg.add("SymEngine")
            using SymEngine
        else
            println("It is a required package. Exiting Finch.")
        end
    end

end
try
    using WriteVTK
catch e
    if force_package_install
        println("WriteVTK package is not yet installed. Installing now.");
        using Pkg
        Pkg.add("WriteVTK")
        using WriteVTK

    else
        println("WriteVTK package is not yet installed. It is optional.\nWould you like to install now? y/n.");
        response = readline();
        if response=="Y" || response=="y"
            using Pkg
            Pkg.add("WriteVTK")
            using WriteVTK
        else
            println("Continuing without WriteVTK. Note that VTK output will not be possible.")
        end
    end

end
try
    using MPI
catch e
    if force_package_install
        println("MPI package is not yet installed. Installing now.");
        using Pkg
        Pkg.add("MPI")
        using MPI
        println("Note that MPI may require some manual configuration to work properly.")
        println("Now attempting to build the MPI package using an MPI implementation installed on your system.");
        ENV["JULIA_MPI_BINARY"] = "system";
        Pkg.build("MPI"; verbose=true);
        global need_restart = true;

    else
        println("MPI package is not yet installed. It is optional.\nWould you like to install now? y/n.");
        response = readline();
        if response=="Y" || response=="y"
            using Pkg
            Pkg.add("MPI")
            using MPI
            println("Note that MPI may require some manual configuration to work properly.")
            println("Now attempting to build the MPI package using an MPI implementation installed on your system.");
            ENV["JULIA_MPI_BINARY"] = "system";
            Pkg.build("MPI"; verbose=true);
            global need_restart = true;
        else
            println("Continuing without MPI.")
        end
    end

end
if @isdefined(MPI)
    try
        using METIS_jll
    catch e
        if force_package_install
            println("METIS_jll package is not yet installed. Installing now.");
            using Pkg
            Pkg.add("METIS_jll")
            using METIS_jll: libmetis

        else
            println("METIS_jll package is not yet installed. It is optional, but needed for effective MPI parallelization.\nWould you like to install now? y/n.");
            response = readline();
            if response=="Y" || response=="y"
                using Pkg
                Pkg.add("METIS_jll")
                using METIS_jll: libmetis
            else
                println("Continuing without METIS_jll.")
            end
        end

    end
end

try
    using PETSc
catch e
    println("PETSc package is not yet installed. It is optional and may cause issues depending on your system.\nWould you like to install now? y/n.");
    response = readline();
    if response=="Y" || response=="y"
        using Pkg
        Pkg.add("PETSc")
        try
            using PETSc
        catch e
            println("There was an issue with using PETSc. Check the Julia PETSc.jl documentation to install correctly.")
            exit(0);
        end
    else
        println("Continuing without PETSc.")
    end

end

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

# include these first
include("finch_constants.jl");
include("finch_config.jl");
include("finch_prob.jl");
include("finch_import_symbols.jl");

# include these next
include("logging.jl");

include("refel.jl");
include("mesh_data.jl");
include("mesh_read.jl");
include("mesh_write.jl")
include("grid.jl");
include("simple_mesh.jl");
include("recursive_ordering.jl");
include("tiled_ordering.jl");
include("ef_ordering.jl");
include("bad_ordering.jl");
include("partition.jl");
include("ghost_exchange.jl");

include("general_utils.jl");
include("function_utils.jl");
include("symexpression.jl");
include("variables.jl");
include("coefficient.jl");
include("parameter.jl");
include("indexer.jl");
include("geometric_factors.jl");
include("dg_utils.jl");
include("time_steppers.jl");
include("fv_utils.jl");
include("fv_neighborhood.jl");
include("grid_parent_child.jl");
include("polyharmonic_interp.jl");

include("output_data.jl");
include("cachsim_solve.jl");

# Finch submodules
include("cachesim/CachesimOut.jl");
using .CachesimOut
include("SymbolicParser.jl")
using .SymbolicParser
include("CodeGenerator.jl");
using .CodeGenerator
include("DGSolver.jl");
using .DGSolver
include("CGSolver.jl");
using .CGSolver
include("FVSolver.jl");
using .FVSolver
include("MixedSolver.jl");
using .MixedSolver

max_stars_repo_name: paralab/Finch

max_issues_repo_name: paralab/Finch

max_forks_repo_name: paralab/Finch


text:
using EconPDEs, Distributions

struct AchdouHanLasryLionsMollModel
    # income process parameters
    κy::Float64
    ybar::Float64
    σy::Float64

    r::Float64

    # utility parameters
    ρ::Float64
    γ::Float64

    amin::Float64
    amax::Float64
end

function AchdouHanLasryLionsMollModel(;κy = 0.1, ybar = 1.0, σy = 0.07, r = 0.03, ρ = 0.05, γ = 2.0, amin = 0.0, amax = 500.0)
    AchdouHanLasryLionsMollModel(κy, ybar, σy, r, ρ, γ, amin, amax)
end

function (m::AchdouHanLasryLionsMollModel)(state::NamedTuple, value::NamedTuple)
    (; κy, σy, ybar, r, ρ, γ, amin, amax) = m
    (; y, a) = state
    (; v, vy_up, vy_down, va_up, va_down, vyy, vya, vaa) = value
    μy = κy * (ybar - y)
    vy = (μy >= 0) ? vy_up : vy_down

    va = va_up
    iter = 0
    @label start
    va_up = max(va_up, eps())
    c = va^(-1 / γ)
    μa = y + r * a - c
    if (iter == 0) & (μa <= 0)
        iter += 1
        va = va_down
        @goto start
    end
    if (a ≈ amin) && (μa <= 0.0)
        va = (y + r * amin)^(-γ)
        c = y + r * amin
        μa = 0.0
    end
    vt = - (c^(1 - γ) / (1 - γ) + μa * va + μy * vy + 0.5 * vyy * σy^2 - ρ * v)
    return (; vt)
end

m = AchdouHanLasryLionsMollModel()
distribution = Gamma(2 * m.κy * m.ybar / m.σy^2, m.σy^2 / (2 * m.κy))
stategrid = OrderedDict(:y => range(quantile(distribution, 0.001), quantile(distribution, 0.999), length = 10),
                        :a =>  range(m.amin, m.amax, length = 100)
                        )
yend = OrderedDict(:v => [log(y + max(a, 0.0)) for y in stategrid[:y], a in stategrid[:a]])
result = pdesolve(m, stategrid, yend)
@assert result.residual_norm <= 1e-5

# finite horizon over 20 years
yend = OrderedDict(:v => [max(a + y)^(1-m.γ)/(1-m.γ) for y in stategrid[:y], a in stategrid[:a]])
τs = range(0, stop = 100, step = 1)
result  = pdesolve(m, stategrid, yend, τs)
@assert maximum(result.residual_norm) <= 1e-5

# Check marginal value of wealth converges to 1.0 at infinity
#b = ((m.r + (m.ρ - m.r)/m.γ))^(1/(1 - 1/m.γ))
#pw = (result[:v] * (1-m.γ)).^(1/(1-m.γ)-1) .* result[:va] ./ b


text:
dictkeys(d::Dict) = (collect(Symbol.(keys(d)))...,)
dictvalues(d::Dict) = (collect(values(d))...,)
namedtuple(d::Dict{T1, T2}) where {T1, T2} = NamedTuple{dictkeys(d)}(dictvalues(d))

"""
```
makeInstanceDict(valType::DataType)
```
Make a dictionary of enumerated variable instances. Helps in translating JSON files.
"""
function makeInstanceDict(valType::DataType)
    insts = instances(valType)
    dict = Dict{String, valType}()
    for inst in insts
        push!(dict, string(inst) => inst)
    end
    return dict
end

"""
```
subTypeTree(t; dict = Dict(), level = 1, cutoff = 0)
```
Create subtype dictionary. Adapted from https://github.com/JuliaLang/julia/issues/24741
"""
function subTypeTree(t; dict = Dict(), level = 1, cutoff = 0)
    level > cutoff ? push!(dict, t => supertype(t)) : nothing
    for s in subtypes(t)
        subTypeTree(s; level = level + 1, dict = dict, cutoff = cutoff)
    end
    return dict
end

"""
```
makeTypeDict(valType::DataType)
```
Inputs contain strings that correspond to DDD data types. This function atuomatically creates a dictionary for all concrete subtypes of a given `valType`.

# Examples

```@jldoctest
julia> abstract type MyAbstractType end
julia> struct MyStruct1 <: MyAbstractType end
julia> struct MyStruct2 <: MyAbstractType end
julia> makeTypeDict(MyAbstractType)
Dict{String,Any} with 4 entries:
  "DDD.MyStruct1()" => MyStruct1()
  "DDD.MyStruct2()" => MyStruct2()
  "MyStruct1()"     => MyStruct1()
  "MyStruct2()"     => MyStruct2()
```
"""
function makeTypeDict(valType::DataType; cutoff = 1)
    primitive = subTypeTree(valType; cutoff = cutoff)
    dict = Dict{String, Any}()
    for (key, val) in primitive
        strSubType = string(key) * "()"
        try
            push!(dict, strSubType => key())
            if strSubType[1:4] == "DDD."
                push!(dict, strSubType[5:end] => key())
            else
                push!(dict, "DDD." * strSubType => key())
            end
        catch
        end
    end
    return dict
end

"""
```
inclusiveComparison(data, args...)::Bool
```
Compare data to a tuple, return `true` if it is equal to any arg, `false` if it is not equal to any.

# Examples

```@jldoctest
julia> inclusiveComparison("f", 1,4,5,"f")
true
julia> inclusiveComparison(23.246, 1.5, 4, 5, "f")
false
```
"""
function inclusiveComparison(data, args...)::Bool
    for i in eachindex(args)
        if data == args[i]
            return true
        end
    end
    return false
end

"""
```
compStruct(arg1, arg2; verbose::Bool = false)
```
Compares values of the fields of two variables `arg1` and `arg2` with the same structure. If `verbose = true`, it will print which fields are different from each other.

# Examples

```@jldoctest
julia> struct MyStruct1; x; end
julia> test1 = MyStruct1(1)
MyStruct1(1)
julia> test2 = MyStruct1(5)
MyStruct1(5)
julia> compStruct(test1, test2; verbose = true)
Structures differ in field: x.
false
julia> compStruct(1, 1; verbose = true)
true
julia> compStruct(1, [1]; verbose = true)
false
```
"""
function compStruct(arg1, arg2; verbose::Bool = false)
    if typeof(arg1) != typeof(arg2)
        !verbose ? nothing :
        @warn "compStruct: Variables have different types:\n\ttypeof(arg1) = $(typeof(arg1))\n\ttypeof(arg2) = $(typeof(arg2))"
        return false
    end
    names = fieldnames(typeof(arg1))
    flag::Bool = true
    for name in names
        result = getproperty(arg1, name) == getproperty(arg2, name)
        if result == false
            flag = false
            !verbose ? nothing : @warn "Structures differ in field: $(name)."
        end
    end
    return flag
end

"""
```
internalAngle(n::Int)
```
Compute the interior angle of a regular polygon with `n` sides.
"""
internalAngle(n::Int) = (n - 2) * π / n

"""
```
externalAngle(n::Int)
```
Compute the exterior angle of a regular polygon with `n` sides.
"""
externalAngle(n::Int) = π - internalAngle(n)

"""
```
rot3D(xyz, uvw, abc, θ)
```
Rotate point `xyz` about the vector `uvw` that crosses point `abc` by the angle `θ`. Further details found [here](https://sites.google.com/site/glennmurray/Home/rotation-matrices-and-formulas/rotation-about-an-arbitrary-axis-in-3-dimensions).

# Examples

```@jldoctest
julia> rot3D([1;1;1],[1;0;0],[0;0;0],π/2)
3-element Array{Float64,1}:
  1.0
 -0.9999999999999999
  1.0
julia> rot3D([1;1;1],[1;0;0],[0;0;0],-π/2)
3-element Array{Float64,1}:
1.0
1.0
-0.9999999999999999
julia> rot3D([1;1;1],[1;0;0],[0;0;0],π)
3-element Array{Float64,1}:
  1.0
 -1.0000000000000002
 -0.9999999999999999
```
"""
function rot3D(xyz, uvw, abc, θ)
    uvw = normalize(uvw)

    cosθ = cos(θ)
    onemcosθ = 1 - cosθ
    sinθ = sin(θ)
    xyzDOTuvw = dot(xyz, uvw)

    return SVector{3, typeof(xyzDOTuvw)}(
        (
            abc[1] * (uvw[2] * uvw[2] + uvw[3] * uvw[3]) -
            uvw[1] * (abc[2] * uvw[2] + abc[3] * uvw[3] - xyzDOTuvw)
        ) * onemcosθ +
        xyz[1] * cosθ +
        sinθ * (-abc[3] * uvw[2] + abc[2] * uvw[3] - uvw[3] * xyz[2] + uvw[2] * xyz[3]),
        (
            abc[2] * (uvw[1] * uvw[1] + uvw[3] * uvw[3]) -
            uvw[2] * (abc[1] * uvw[1] + abc[3] * uvw[3] - xyzDOTuvw)
        ) * onemcosθ +
        xyz[2] * cosθ +
        sinθ * (abc[3] * uvw[1] - abc[1] * uvw[3] + uvw[3] * xyz[1] - uvw[1] * xyz[3]),
        (
            abc[3] * (uvw[1] * uvw[1] + uvw[2] * uvw[2]) -
            uvw[3] * (abc[1] * uvw[1] + abc[2] * uvw[2] - xyzDOTuvw)
        ) * onemcosθ +
        xyz[3] * cosθ +
        sinθ * (-abc[2] * uvw[1] + abc[1] * uvw[2] - uvw[2] * xyz[1] + uvw[1] * xyz[2]),
    )
end

"""
```
⊗(x::AbstractVector, y::AbstractVector)
```
Tensor product.
"""
⊗(x::AbstractVector{T1}, y::AbstractVector{T2}) where {T1, T2} = x * y'

"""
```
linePlaneIntersect(n::T, p0::T, l::T, l0::T) where {T <: AbstractVector}
```
Finds the intersect between a line and a plane. `n` is the plane normal, `p0` is a point on the plane, `l` is the line vector, `l0` is a point on a line.
"""
function linePlaneIntersect(n, p0, l, l0)
    den = l ⋅ n
    num = (p0 - l0) ⋅ n

    if den ≈ 0
        num ≈ 0 && return Inf
        return nothing
    end

    d = num / den
    intersect = l0 + d * l
    return intersect
end

"""
```
gausslegendre(n::Integer, a, b)
```
Compute Gauss-Legendre quadrature points and weights for the interval `[a, b]`.
"""
function FastGaussQuadrature.gausslegendre(n::Integer, a, b)
    bma = (b - a) / 2
    bpa = (b + a) / 2
    x, w = gausslegendre(n)
    @inbounds @simd for i in 1:n
        x[i] = bma * x[i] + bpa
        w[i] = bma * w[i]
    end
    return x, w
end
function gausslegendre2D(n)
    x, w = gausslegendre(n)
    X = zeros(2, n^2)
    W = zeros(n^2)
    @inbounds @simd for i in 1:n
        idx = n * (i - 1)
        for j in 1:n
            X[:, j + idx] .= (x[j], x[i])
            W[j + idx] = w[j] * w[i]
        end
    end
    return X, W
end
function gausslegendre3D(n)
    x, w = gausslegendre(n)
    X = zeros(3, n^3)
    W = zeros(n^3)
    @inbounds @simd for i in 1:n
        idx = n^2 * (i - 1)
        for j in 1:n
            idx2 = idx + n * (j - 1)
            for k in 1:n
                X[:, k + idx2] .= (x[k], x[j], x[i])
                W[k + idx2] = w[k] * w[j] * w[i]
            end
        end
    end
    return X, W
end

"""
```
safeNorm(x)
```
Safe normalisation.
"""
function safeNorm(x)
    elemT = eltype(x)
    xNorm = norm(x)
    if xNorm > eps(elemT)
        x = normalize(x)
    else
        x = zero(x)
        xNorm = eps(elemT)
    end
    return xNorm, x
end

"""
```
minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
```
Calculates the minimum distance between two segments `seg1 = x0 → x1`, `seg2 = y0 → y1`.

# Returns

- `distSq`: minimum distance squared between `L1` and `L2`.
- `dDistSqDt`: rate of change with respect to time of the minimum distance squared between `L1` and `L2`.
- `L1`: normalised position on `seg1` that is closest to `seg2`.
- `L2`: normalised position on `seg2` that is closest to `seg1`.
"""
function minimumDistance(x0, x1, y0, y1, vx0, vx1, vy0, vy1)
    seg1 = x1 - x0
    seg2 = y1 - y0
    vseg1 = vx1 - vx0
    vseg2 = vy1 - vy0

    A = seg1 ⋅ seg1
    B = 2 * seg1 ⋅ (x0 - y0)
    C = 2 * seg1 ⋅ seg2
    D = 2 * seg2 ⋅ (y0 - x0)
    E = seg2 ⋅ seg2
    F = x0 ⋅ x0 + y0 ⋅ y0
    G = C^2 - 4 * A * E

    elemT = typeof(G)

    # seg1 is a point.
    if A < eps(elemT)
        L1 = 0
        # seg2 is a point.
        E < eps(elemT) ? L2 = 0 : L2 = -D / E / 2
        # seg2 is a point.
    elseif E < eps(elemT)
        L2 = 0
        # seg1 is a point.
        A < eps(elemT) ? L1 = 0 : L1 = -B / A / 2
        # seg1 and two are close to each other.
    elseif abs(G) < eps(elemT)
        dist = SVector{4, elemT}(
            (y0 - x0) ⋅ (y0 - x0),
            (y1 - x0) ⋅ (y1 - x0),
            (y0 - x1) ⋅ (y0 - x1),
            (y1 - x1) ⋅ (y1 - x1),
        )

        missing, idx = findmin(dist)

        L1 = floor(idx / 2)
        L2 = mod(idx - 1, 2)
        # seg1 and seg2 are not points and away from each other.
    else
        L2 = (2 * A * D + B * C) / G
        L1 = (C * L2 - B) / A / 2
    end

    L1 = clamp(L1, 0, 1)
    L2 = clamp(L2, 0, 1)

    distSq = (x0 + seg1 * L1 - y0 - seg2 * L2) ⋅ (x0 + seg1 * L1 - y0 - seg2 * L2)
    dDistSqDt =
        2 * ((vx0 + vseg1 * L1 - vy0 - vseg2 * L2) ⋅ (x0 + seg1 * L1 - y0 - seg2 * L2))

    return distSq, dDistSqDt, L1, L2
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
#
#	This file contains some LinearAlgebraicRepresentation function that
#	 have not been updated to the last version yet. The problem is raised
#	 the `cat` acronym, here renamed `pcat`.
#	This file will be removed as soon as the appropriate changes are made.
#

function chainbasis2polygons(V,copEV,copFE)
	FE = [findnz(copFE[k,:])[1] for k=1:copFE.m]
	EV = [findnz(copEV[k,:])[1] for k=1:copEV.m]

	FEs = Array{Int64,1}[]
	EVs = Array{Array{Int64,1},1}[]
	FVs = Array{Int64,1}[]
	for f=1:copFE.m
		push!( FEs, collect(Set(CAGD.pcat([e for e in FE[f]]))) )
		# edges in EVs are aggregated by face, in order to answer point-classifications
		push!( EVs, [EV[e] for e in FE[f]]  )
		push!( FVs, collect(Set(CAGD.pcat([EV[e] for e in FE[f]])))  )
	end
	polygons = collect(zip(EVs,FVs,FEs)) # (EV,FV,FFE)s for polygon
	W = convert(Lar.Points,V')
	return W,polygons,FE
end

function internalpoints2d(W,copEV,copFE) # W by rows
	# transform each 2-cell in a solid (via Lar model)
	#-------------------------------------------------------------------------------
	U,pols,FE = chainbasis2polygons(W,copEV,copFE) # V by rows
	# compute, for each `pol` (3-cell) in `pols`, one `internalpoint`.
	#-------------------------------------------------------------------------------
	internalpoints = []
	for f=1:length(pols)
		(EV,FV,FE) = pols[f]
		#GL.VIEW([ GL.GLFrame, GL.GLLines(V,EV) ]);
		internalpoint = Lar.getinternalpoint2d(W,EV,FV, f, copEV,copFE)
		push!(internalpoints,internalpoint)
	end
	return internalpoints
end


text:
function one_plus_lambda(layers::Array, num_generations::Int, f::Function, 𝜆::Int)

    fits = zeros(num_generations)
    expert = Individual(layers::Array)
    evaluate!(expert, f)

    for i in eachindex(fits)
        population = Array{Individual}(undef, 𝜆)
        best = 1
        for j in eachindex(population)
            population[j] = mutate(expert)
            evaluate!(population[j], f)
            if population[j].fitness > population[best].fitness
                best = j
            end
        end

        if population[best].fitness >= expert.fitness
            expert = population[best]
        end

        fits[i] = expert.fitness
        print(i, "/",  num_generations, " > ", expert.fitness, "\t\t\r")
    end
    fits, expert
end

max_stars_repo_name: TemplierPaul/evolution

max_issues_repo_name: TemplierPaul/evolution

max_forks_repo_name: TemplierPaul/evolution


text:
using ImageBinarization
using Test, TestImages
using Statistics
using ImageCore, ImageTransformations
using ReferenceTests

include("testutils.jl")

@testset "ImageBinarization.jl" begin
    include("util.jl")

    include("algorithms/adaptive_threshold.jl")
    include("algorithms/balanced.jl")
    include("algorithms/entropy.jl")
    include("algorithms/intermodes.jl")
    include("algorithms/minimum.jl")
    include("algorithms/minimum_error.jl")
    include("algorithms/moments.jl")
    include("algorithms/niblack.jl")
    include("algorithms/otsu.jl")
    include("algorithms/polysegment.jl")
    include("algorithms/sauvola.jl")
    include("algorithms/unimodal.jl")
    include("algorithms/yen.jl")
end

nothing


text:
using NFFT, FFTW

FFTW.set_num_threads(Threads.nthreads())

# This is a reduced version of the actual problem
function _nfft!(p, f)

  NFFT.deconvolve!(p, f, p.tmpVec)     # deconvolve uses Polyester-based threading or no threading depending on storeDeconvolutionIdx

  @info "Timing After First FFT "
  @time p.forwardFFT * p.tmpVec
  @info "Timing After Second FFT"
  @time p.forwardFFT * p.tmpVec

 return
end

function doTimings()
  N = 1024;
  p1 = plan_nfft(zeros(2,N), (N,N); fftflags=NFFT.FFTW.MEASURE, storeDeconvolutionIdx=false);
  f = zeros(ComplexF64,(N,N));

  # warmup FFT
  p1.forwardFFT * p1.tmpVec

  @info "\n\nWith Deconvization Variant 1"
  _nfft!(p1, f);

  return
end

doTimings() # compiling
doTimings()


text:
export
    parent_ranges,
    get_range,
    set_range!,
    set_ranges!,
    expander_range

import Scruff: get_range

#####################################################################
#                                                                   #
# Algorithm building blocks for setting and expanding the ranges of #
# all the current instances in a runtime.                           #
#                                                                   #
# order is a topological sort of the network in the runtime.        #
# size is a uniform target size for each instance.                  #
#                                                                   #
#####################################################################

"""
    parent_ranges(runtime::Runtime, var::Variable{I,J,O}, depth = typemax(Int)) where {I,J,O}

Returns the ranges of the parents of the given variable.

See [`Scruff.get_range`](@ref)
"""
function parent_ranges(runtime::Runtime, var::Variable{I,J,O}, depth = typemax(Int)) where {I,J,O}
    net = get_network(runtime)
    pars = get_parents(net, var)
    result = []
    for p in pars
        pinst = current_instance(runtime, p)
        push!(result, get_range(runtime, pinst, depth))
    end
    T = isempty(result) ? Vector{O} : typejoin([typeof(x) for x in result]...)
    return convert(Vector{T}, result)
end

"""
    get_range(runtime::DynamicRuntime, v::Variable{I,J,O}, depth = 1) where {I,J,O}

Returns the range of the most recent instance of the given variable.

See [`Scruff.get_range`](@ref)
"""
function get_range(runtime::DynamicRuntime, v::Variable{I,J,O}, depth = 1) where {I,J,O}
    inst = current_instance(runtime, v)
    range = get_range(runtime, inst, depth)
    if range !== nothing
        return range
    elseif has_previous_instance(runtime, v)
        prev = previous_instance(runtime, v)
        if has_range(runtime, prev)
            return get_range(runtime, prev, depth)
        else
            return O[]
        end
    else
        return O[]
    end
end

"""
    get_range(runtime::InstantRuntime, v::Node{O}, depth = 1) where O

Returns the range of the given node.

See [`Scruff.get_range`](@ref)
"""
function get_range(runtime::InstantRuntime, v::Node{O}, depth = 1) where O
    inst = runtime.instances[v]
    range = get_range(runtime, inst, depth)
    if range !== nothing
        return range
    else
        return O[]
    end
end

"""
    set_ranges!(runtime::InstantRuntime, evidence = Dict{Symbol, Score}(),
        size = 10, depth = 1,
        order = topsort(get_initial_graph(get_network(runtime))),
        placeholder_beliefs = get_placeholder_beliefs(runtime))

Set the ranges of all current instances in the runtime.

This method first checks whether ranges exist for the runtime at the desired depth,
with the desired range size, and with the same evidence. If so, it doesn't do anything.
If the depth is less than 1, it doesn't do anything.
Otherwise, it uses the support operator to compute ranges of variables in order.
Placeholders should have ranges set already in `placeholder_beliefs`.
For expanders, it recursively sets the ranges of the subnetworks at depth - 1.

Returns a flag indicating whether any instance has a changed range.
"""
function set_ranges!(runtime::InstantRuntime,
    evidence::Dict{Symbol,Score} = Dict{Symbol, Score}(),
    size :: Int = 10, depth :: Int = 1,
    order = topsort(get_initial_graph(get_network(runtime))),
    placeholder_beliefs = get_placeholder_beliefs(runtime))
    if depth < 1
        return false
    end
    nodes = get_nodes(get_network(runtime))
    # If we require greater depth or greater size or different evidence, we have to do it again
    if has_state(runtime, :nodes) && has_state(runtime, :range_size) && has_state(runtime, :range_depth) && has_state(runtime, :range_evidence) &&
        get_state(runtime, :range_size) >= size  && get_state(runtime, :range_depth) >= depth && get_state(runtime, :range_evidence) == evidence &&
        get_state(runtime, :nodes) == nodes
        return false
    else
        set_state!(runtime, :range_size, size)
        set_state!(runtime, :range_depth, depth)
        set_state!(runtime, :range_evidence, evidence)
        set_state!(runtime, :nodes, nodes)
    end
    changed = false
    for v in order
        O = output_type(v)
        rng::Vector{O} = O[]
        curr = get_range(runtime, v, depth)
        inst = current_instance(runtime, v)
        if v isa Placeholder
            pi = placeholder_beliefs[v.name]
            rng = support(pi, (), size, output_type(v)[])
            chng = rng != curr
        else
            sf = get_sfunc(inst)
            if length(curr) < size
                parranges = parent_ranges(runtime, v, depth)
                if isa(sf, Expander)
                    rc = expander_range(runtime, v, size, depth)
                    rng= rc[1]
                    chng = rc[2]
                else
                    rng = support(sf, tuple(parranges...), size, curr)
                    chng = rng != curr
                end
                if v.name in keys(evidence)
                    ev = evidence[v.name]
                    if ev isa HardScore
                        if !(ev.value in rng)
                            push!(rng, ev.value)
                            chng = true
                        end
                    elseif ev isa LogScore
                        for k in keys(ev.logscores)
                            if !(k in rng)
                                push!(rng, k)
                                chng = true
                            end
                        end
                    end
                end
                changed = changed || chng
            else
                rng = curr
            end
        end
        set_range!(runtime, inst, rng, depth)
    end
    if changed
        delete_state!(runtime, :solution)
    end
    return changed
end

"""
    function expander_range(runtime :: Runtime, v :: Variable,
        target_size :: Int, depth :: Int)

    Recursively compute the range of the expander and subnetworks up to the given depth.

Computing the range of the expander expands enough of the parent
range to reach the desired target size, or expands all the parents fully.
"""
function expander_range(runtime :: Runtime, v :: Variable,
                        target_size :: Int, depth :: Int)
    changed = false
    net = runtime.network
    parents = get_parents(net, v)
    parent = parents[1]
    parrange = get_range(runtime, parent, depth)
    rangeset = Set()
    subranges = Dict()
    for p in parrange
        if expanded(runtime, v, p)
            subnet = expansion(runtime, v, p)
            subrun = subruntime(runtime, v, subnet)
            subout = current_instance(subrun, get_node(subnet, :out))
            if has_range(subrun, subout, depth - 1)
                subrange = get_range(subrun, subout, depth - 1)
                subranges[p] = subrange
                union!(rangeset, subrange)
            else
                subranges[p] = []
            end
        else
            subranges[p] = []
        end
    end
    todo = Dict()
    for p in parrange
        todo[p] = length(subranges[p])
    end
    while length(rangeset) < target_size && !isempty(todo)
        p = argmin(todo)
        delete!(todo, p)
        (subnet, chng) = expand!(runtime, v, p)
        changed = chng
        subrun = subruntime(runtime, v, subnet)
        if depth > 1
            order = topsort(get_initial_graph(subnet))
            outvar = get_node(subrun, :out)
            # Cannot have evidence on subnet
            if set_ranges!(subrun, Dict{Symbol, Score}(), target_size - length(rangeset), depth - 1, order)
                changed = true
            end
            inst = current_instance(subrun, outvar)
            rng = get_range(subrun, inst, depth - 1)
            union!(rangeset, rng)
        end
    end
    result = output_type(v)[]
    for x in rangeset
        push!(result, x)
    end
    sort!(result)
    return (result, changed)
end


text:
module common

    export trace_buffer_parts
    """
    buf_sizes, trace_idx = trace_buffer_parts(nbufs::Integer,ntraces::Integer)

    """
    function trace_buffer_parts(nbufs::Integer,ntraces::Integer)
        part = Vector{Int}(undef,nbufs)
        idxs = Vector{Int}(undef,nbufs+1)
        r::Int = rem(ntraces,nbufs)
        c::Int = ceil(Int,ntraces/nbufs)
        f::Int = floor(Int,ntraces/nbufs)
        part[1:r]       .= c
        part[r+1:nbufs] .= f
        @assert sum(part)==ntraces "FATAL ERROR: failed to properly partition $ntraces to $nbufs workers"
        for i=0:nbufs idxs[i+1]=sum(part[1:i])+1 end
        return part,idxs
    end

end


text:
import GitHub
using Nanosoldier, Test, BenchmarkTools
using Nanosoldier: BuildRef, JobSubmission, Config, BenchmarkJob, PkgEvalJob, AbstractJob
using BenchmarkTools: TrialEstimate, Parameters

#########
# setup #
#########

vinfo = """
Julia Version 0.4.3-pre+6
Commit adffe19 (2015-12-11 00:38 UTC)
Platform Info:
  System: Darwin (x86_64-apple-darwin13.4.0)
  CPU: Intel(R) Core(TM) i5-4288U CPU @ 2.60GHz
  WORD_SIZE: 64
  uname: Darwin 13.4.0 Darwin Kernel Version 13.4.0: Wed Dec 17 19:05:52 PST 2014; root:xnu-2422.115.10~1/RELEASE_X86_64 x86_64 i386
Memory: 8.0 GB (1025.07421875 MB free)
Uptime: 646189.0 sec
Load Avg:  1.38427734375  1.416015625  1.41455078125
Intel(R) Core(TM) i5-4288U CPU @ 2.60GHz:
       speed         user         nice          sys         idle          irq
#1  2600 MHz     186848 s          0 s     114241 s    1462955 s          0 s
#2  2600 MHz     114716 s          0 s      43882 s    1605408 s          0 s
#3  2600 MHz     189864 s          0 s      77629 s    1496513 s          0 s
#4  2600 MHz     114652 s          0 s      42497 s    1606856 s          0 s

  BLAS: libopenblas (USE64BITINT DYNAMIC_ARCH NO_AFFINITY Haswell)
  LAPACK: libopenblas64_
  LIBM: libopenlibm
  LLVM: libLLVM-3.3
"""

primary = BuildRef("ararslan/julia", "25c3659d6cec2ebf6e6c7d16b03adac76a47b42a", vinfo)
against = BuildRef("JuliaLang/julia", "bb73f3489d837e3339fce2c1aab283d3b2e97a4c", vinfo*"_against")
config = Config("user", Dict(Any => [getpid()]), GitHub.AnonymousAuth(), "test");
tagpred = "ALL && !(\"tag1\" || \"tag2\")"
pkgsel = "[\"Example\"]"

#####################################
# submission parsing and validation #
#####################################

function build_test_submission(jobtyp, submission_string)
    func, args, kwargs = Nanosoldier.parse_submission_string(submission_string)
    submission = JobSubmission(config, primary, primary.sha, "https://www.test.com", :commit, nothing, func, args, kwargs)
    @test Nanosoldier.isvalid(submission, jobtyp)
    return jobtyp(submission)
end

build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL)`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred)`")

build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, vs = \"JuliaLang/julia:master\")`")

build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, isdaily = true, vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", isdaily = true, vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, isdaily = true, vs = \"JuliaLang/julia:master\")`")

build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL; isdaily = true, vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\"; isdaily = true, vs = \"JuliaLang/julia:master\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred; isdaily = true, vs = \"JuliaLang/julia:master\")`")

build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, vs = \"JuliaLang/julia#v1.0.0\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(\"tag\", vs = \"JuliaLang/julia#v1.0.0\")`")
build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred, vs = \"JuliaLang/julia#v1.0.0\")`")

build_test_submission(PkgEvalJob, "@nanosoldier `runtests(ALL)`")
build_test_submission(PkgEvalJob, "@nanosoldier `runtests(\"pkg\")`")
build_test_submission(PkgEvalJob, "@nanosoldier `runtests($pkgsel)`")

#############################
# retrieval from job queue  #
#############################

non_daily_job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL)`")
daily_job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks(ALL, isdaily = true)`")

queue = [daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && job.isdaily
@test length(queue) == 1

queue = [non_daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && job.isdaily
@test length(queue) == 1

queue = [daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, false)
@test job === nothing
@test length(queue) == 2

queue = [non_daily_job, daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, false)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, false)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, PkgEvalJob, true)
@test job === nothing
@test length(queue) == 2

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, BenchmarkJob, true)
@test job !== nothing && job.isdaily
@test length(queue) == 1

queue = [daily_job, non_daily_job]
job = Nanosoldier.retrieve_job!(queue, Any, false)
@test job !== nothing && !job.isdaily
@test length(queue) == 1

#########################
# job report generation #
#########################

job = build_test_submission(BenchmarkJob, "@nanosoldier `runbenchmarks($tagpred)`")
@test job.tagpred == tagpred
@test job.against == nothing
job.against = against

results = Dict(
    "primary" => BenchmarkGroup([],
        "g" => BenchmarkGroup([],
            "h" => BenchmarkGroup([],
                "x"      => TrialEstimate(Parameters(), 1.0, 3.5, 1.0, 1.0),  # invariant
                ("y", 1) => TrialEstimate(Parameters(memory_tolerance = 0.03), 2.0, 1.0, 0.0, 1.0),  # regression/improvement
                ("y", 2) => TrialEstimate(Parameters(time_tolerance = 0.04), 0.5, 1.0, 1.0, 1.0),  # improvement
                "z"      => TrialEstimate(Parameters(memory_tolerance = 0.27, time_tolerance = 0.6), 1.0, 1.0, 5.0, 1.0),  # regression
                "∅"      => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0) # not in "against" group
            )
        )
    ),
    "against" => BenchmarkGroup([],
        "g" => BenchmarkGroup([],
            "h" => BenchmarkGroup([],
                "x"      => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0),
                ("y", 1) => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0),
                ("y", 2) => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0),
                "z"      => TrialEstimate(Parameters(), 1.0, 1.0, 1.0, 1.0)
            )
        )
    )
)

results["judged"] = BenchmarkTools.judge(results["primary"], results["against"])

@test begin
    mdpath = joinpath(@__DIR__, "report.md")
    open(mdpath, "r") do file
        read(file, String) == sprint(io -> Nanosoldier.printreport(io, job, results))
    end
end


text:
function test_user_mean()
    println("-- Testing UserMean recommender")

    data = DataAccessor(sparse([1 2 3; 4 5 6]))
    recommender = UserMean(data)
    build!(recommender)
    actual = predict(recommender, 1, 1)

    @test actual == 2.0
end

test_user_mean()


text:
global pcollect_opt_methods = [MultiobjOptimizerParams()]

@testset "pcollect(opt_params=$opt_params)" for opt_params in pcollect_opt_methods
    mosaics = Dict("A" => SetMosaic([Set([:a]), Set([:b]), Set([:c]), Set([:a, :b, :c])],
                                    Set([:a, :b, :c, :d, :e])))
    masks = Dict(1 => Set([:a, :b]),
                 2 => Set([:d, :e]),
                 3 => Set([:a, :c, :d]))
    res = OptEnrichedSetCover.pcollect(mosaics, masks, mode=:sequential,
                                       opt_params=opt_params,
                                       cover_params=CoverParams(sel_tax=-log(0.75)),
                                       enum_params=CoverEnumerationParams(max_set_score=10.0))
    @test length(res) == 2
    @test isa(res, Dict{Tuple{String, Int}, CoverCollection})
end


text:
module ContinuumIceModels

export ContinuumIceModel, time_step!, ViscoElasticRheology

using KernelAbstractions

using Oceananigans.AbstractOperations: ∂x, ∂y
using Oceananigans.Advection: CenteredSecondOrder
using Oceananigans.Architectures: device_event, device
using Oceananigans.BoundaryConditions: fill_halo_regions!
using Oceananigans.Fields: XFaceField, YFaceField, CenterField, Field, Center, Face
using Oceananigans.Models: AbstractModel
using Oceananigans.TimeSteppers: Clock, tick!
using Oceananigans.Utils: launch!

using Oceananigans.Operators

# Simulations interface
import Oceananigans: fields
import Oceananigans.TimeSteppers: time_step!, update_state!
import Oceananigans.Simulations: reset!

mutable struct ContinuumIceModel{Grid,
                                 Tim,
                                 Clo,
                                 Adv,
                                 Rheo,
                                 Drag,
                                 Vel,
                                 Tracers,
                                 Stresses,
                                 Ocean,
                                 Tend} <: AbstractModel{Nothing}
    grid :: Grid
    timestepper :: Tim # silly Oceananigans
    clock :: Clo
    tracer_advection :: Adv
    rheology :: Rheo
    drag :: Drag
    velocities :: Vel
    tracers :: Tracers # concentration, thickness, snow
    stresses :: Stresses
    ocean :: Ocean # ocean.u, ocean.v, ocean.ρ
    tendencies :: Tend
end

# Oceananigans.Simulations interface
fields(m::ContinuumIceModel) = merge(m.velocities, m.stresses)
update_state!(m::ContinuumIceModel) = fill_halo_regions!(fields(m), m.grid.architecture)
reset!(::Nothing) = nothing # silly Oceananigans # ⟨⟨ eyes ⟩⟩

"""
    ContinuumIceModel(; grid, ocean, rheology, kw...)

Return a continuum model for sea ice on an `ocean` state with `rheology` and `grid`.
"""
function ContinuumIceModel(; grid, ocean,
                           clock = Clock{eltype(grid)}(0, 0, 1),
                           tracer_advection = CenteredSecondOrder(),
                           rheology = ViscoElasticRheology(0.25, 1e3))

    # Check that grid is 2D with z=0 OR use 3D grid but place fields at the surface.

    velocities = (u = XFaceField(grid), v = YFaceField(grid))

    # One tracer for now
    tracers = (concentration = CenterField(grid), thickness=1)

    stresses = (σ₁₁ = CenterField(grid),
                σ₁₂ = Field{Face, Face, Center}(grid),
                σ₂₂ = CenterField(grid))

    drag = (; Cᴰ_ocean = 1e-3)

    tendencies = (
        u = XFaceField(grid),
        v = YFaceField(grid),
        σ₁₁ = CenterField(grid),
        σ₁₂ = Field{Face, Face, Center}(grid),
        σ₂₂ = CenterField(grid)
    )

    return ContinuumIceModel(grid,
                             nothing,
                             clock,
                             tracer_advection,
                             rheology,
                             drag,
                             velocities,
                             tracers,
                             stresses,
                             ocean,
                             tendencies)
end

"""
    ViscoElasticRheology(modulus=0.25, viscosity=1e3)

Return `ViscoElasticRheology` with elastic `modulus`
and dynamic `viscosity` representing viscoelastic ice model in which
the ice stress tensor obeys a prognostic viscoelastic
evolution equation.
"""
Base.@kwdef struct ViscoElasticRheology{T}
    modulus :: T = 0.25 # not old but Young
    viscosity :: T = 1e3
end

# Has free parameters...
Base.@kwdef struct Hibler97Pressure{T}
    p★ :: T = 2.75e5
    c★ :: T = 20.0
    e :: T = 2.0 # yield curve axis ratio
end

#####
##### Time-stepping
#####

@inline ϵ₁₁ᶜᶜᶜ(i, j, k, grid, U) = ∂xᶜᶜᶜ(i, j, k, grid, U.u)
@inline ϵ₁₂ᶠᶠᶜ(i, j, k, grid, U) = (∂yᶠᶠᶜ(i, j, k, grid, U.u) + ∂xᶠᶠᶜ(i, j, k, grid, U.v)) / 2
@inline ϵ₂₂ᶜᶜᶜ(i, j, k, grid, U) = ∂yᶜᶜᶜ(i, j, k, grid, U.v)

""" Calculate the right-hand-side of the free surface displacement (η) equation. """
@kernel function _calculate_tendencies!(R, grid, Uᵢ, Uₒ, σ, E, μ, Cᴰₒ, ρₒ)
    σ₁₁, σ₁₂, σ₂₂ = σ
    i, j = @index(Global, NTuple)
    k = grid.Nz

    @inbounds begin
        # Stress
        # calculate_stress_tendency!(R, i, j, k, grid, rheology, Uᵢ)
        R.σ₁₁[i, j, k] = E * (ϵ₁₁ᶜᶜᶜ(i, j, k, grid, Uᵢ) - σ₁₁[i, j, k] / 2μ)
        R.σ₁₂[i, j, k] = E * (ϵ₁₂ᶠᶠᶜ(i, j, k, grid, Uᵢ) - σ₁₂[i, j, k] / 2μ)
        R.σ₂₂[i, j, k] = E * (ϵ₂₂ᶜᶜᶜ(i, j, k, grid, Uᵢ) - σ₂₂[i, j, k] / 2μ)

        # Momentum
        # calculate_momentum_tendency!(R, i, j, k, grid, rheology, Uᵢ, Uₒ, Uₐ)
        x_stress = ∂xᶠᶜᶜ(i, j, k, grid, σ₁₁) + ∂yᶠᶜᶜ(i, j, k, grid, σ₁₂)
        y_stress = ∂yᶜᶠᶜ(i, j, k, grid, σ₂₂) + ∂xᶜᶠᶜ(i, j, k, grid, σ₁₂)

        ρᵢ = 910.0 # ice density TODO: make settable
        uᵣ = Uₒ.u[i, j, k] - Uᵢ.u[i, j, k]
        vᵣ = Uₒ.v[i, j, k] - Uᵢ.v[i, j, k]
        x_drag = Cᴰₒ * ρₒ / ρᵢ * sqrt(uᵣ^2 + vᵣ^2) * uᵣ
        y_drag = Cᴰₒ * ρₒ / ρᵢ * sqrt(uᵣ^2 + vᵣ^2) * vᵣ

        R.u[i, j, k] = x_stress + x_drag
        R.v[i, j, k] = y_stress + y_drag

        # Tracers...
    end
end

function calculate_tendencies!(R, grid, Uᵢ, Uₒ, σ, E, μ, Cᴰₒ, ρₒ)
    arch = grid.architecture
    event = launch!(arch, grid, :xy,
                    _calculate_tendencies!,
                    R, grid, Uᵢ, Uₒ, σ, E, μ, Cᴰₒ, ρₒ,
                    dependencies = device_event(arch))

    wait(device(arch), event)
    return nothing
end

function calculate_tendencies!(model)
    # Unpack parameters
    E   = model.rheology.modulus
    μ   = model.rheology.viscosity
    Cᴰₒ = model.drag.Cᴰ_ocean
    ρₒ  = model.ocean.ρ

    calculate_tendencies!(model.tendencies,
                          model.grid,
                          model.velocities,
                          model.ocean,
                          model.stresses,
                          E, μ, Cᴰₒ, ρₒ)

    #=
    # Rate of strain tensor
    ϵ = (∂x(uᵢ), (∂y(uᵢ) + ∂x(vᵢ)) / 2, ∂y(vᵢ))

    # Relative velocities
    uᵣ = uₒ - uᵢ
    vᵣ = vₒ - vᵢ

    # σ₁₁, σ₁₂, σ₂₂ = σ
    # ϵ₁₁, ϵ₁₂, ϵ₂₂ = ϵ
    # R.σ₁₁ .= E * (ϵ₁₁ - σ₁₁ / 2μ)
    # R.σ₁₂ .= E * (ϵ₁₂ - σ₁₂ / 2μ)
    # R.σ₂₂ .= E * (ϵ₂₂ - σ₂₂ / 2μ)
    =#

    return nothing
end

@kernel function step_fields!(F, grid, dF, Δt)
    i, j = @index(Global, NTuple)
    k = grid.Nz
    @inbounds begin
        F.u[i, j, k]   += Δt * dF.u[i, j, k]
        F.v[i, j, k]   += Δt * dF.v[i, j, k]
        F.σ₁₁[i, j, k] += Δt * dF.σ₁₁[i, j, k]
        F.σ₁₂[i, j, k] += Δt * dF.σ₁₂[i, j, k]
        F.σ₂₂[i, j, k] += Δt * dF.σ₂₂[i, j, k]
    end
end

function time_step!(model, Δt)
    calculate_tendencies!(model)

    arch = model.grid.architecture
    event = launch!(arch, model.grid, :xy,
                    step_fields!,
                    fields(model), model.grid, model.tendencies, Δt,
                    dependencies = device_event(arch))

    wait(device(arch), event)

    update_state!(model)
    tick!(model.clock, Δt)

    return nothing
end

end # module