text:
module RocketEnumerateOperatorTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "operator: enumerate()" begin

    println("Testing: operator enumerate()")

    run_proxyshowcheck("Enumerate", enumerate())

    run_testset([
        (
            source      = from([ 3, 2, 1 ]) |> enumerate(),
            values      = @ts([ (1, 3), (2, 2), (3, 1), c ]),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = completed(Int) |> enumerate(),
            values      = @ts(c),
            source_type = Tuple{Int, Int}
        ),
        (
            source      = faulted(Float64, "e") |> enumerate(),
            values      = @ts(e("e")),
            source_type = Tuple{Int, Float64}
        ),
        (
            source      = never(String) |> enumerate(),
            values      = @ts(),
            source_type = Tuple{Int, String}
        )
    ])

end

end


text:
# Reading of each struct type has 2 methods
# - s has same type as sink
#   - read is performed at the element level and doesn't require any additional permutations
# - s has different type from sink
#   - read in as raw eltype and permutedims to fill in sink. This requires appropriately axisnames

########
# Read #
########
Base.read!(s::AxisStream, sink::ImageMeta) = read!(s, getfield(sink, :data))
Base.read!(s::AxisStream{S,<:BitTypes}, sink::Array{<:BitTypes}) where {S,N} = read!(getfield(s,:data), sink)
Base.read!(s::AxisStream{S,T,N}, sink::AxisArray{T,N}) where {S,T,N} = read!(getfield(s, :data), getfield(s,:data))
Base.read(s::AxisStream{S,T,N,I,Ax}) where {S,T,N,I,Ax} = AxisArray(read(s), axisparams(s))
# raw eltype → struct view can use Axis name permutation to semiautomate this
Base.read!(s::AxisStream{S,<:BitTypes,N}, sink::AxisArray{T,N}) where {S,T,N} = read!(s, getfield(sink, :data))

# Triangle
function Base.read!(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Triangle}) where {S,T,N,I,Ax}
    perm = permutation((:triangledim, setdiff(axisnames(s), :triangledim)...,), axisnames(s))
    return fill!(sink, triangleview(permutedims(read(getfield(s,:data)), perm)))
end

# Point
function Base.read!(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Point}) where {S,N,I,Ax}
    perm = permutation((:pointdim, setdiff(AxisArrays.axisnames(s), :pointdim)...,), AxisArrays.axisnames(s))
    return fill!(sink, pointview(permutedims(read(s), perm)))
end

# Quat
function Base.read!(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Quat}) where {S,N,I,Ax}
    perm = permutation((:quatdim, setdiff(axisnames(s), :quatdim)...,), axisnames(s))
    return fill!(sink, quatview(permutedims(read(getfield(s,:data)), perm)))
end

# SMatrix
function Base.read!(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{SMatrix}) where {S,N,I,Ax}
    perm = permutation((:matdim1, :matdim2, setdiff(axisnames(s), (:matdim1,:matdim2))...,), axisnames(s))
    return fill!(sink, matview(permutedims(read(getfield(s,:data)), perm)))
end

# SVector
function Base.read!(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{SVector}) where {S,N,I,Ax}
    perm = permutation((:vecdim, setdiff(axisnames(s), :vecdim)...,), axisnames(s))
    return fill!(sink, vecview(permutedims(read(getfield(s, :data)), perm)))
end

# ColorTypes
function Base.read!(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{T,N}) where {S,T<:Union{Color3,ColorAlpha},N,I,Ax}
    perm = permutation((:colordim, setdiff(axisnames(s), :colordim)...,), axisnames(s))
    return fill!(sink, colorview(T, PermutedDimsArray(read!(getfield(s, :data)), perm)))
end

# Distributions
function Base.read!(s::AxisStream{S,<:BitTypes}, sink::Array{T}) where {S,T<:Distribution}
    perm = permutation((:statdim, setdiff(axisnames(s), :statdim)...,), axisnames(s))
    return fill!(sink, statview(T, permutedims(read(s), perm)))
end

#########
# Write #
#########
Base.write(s::AxisStream{S,T,N,I,Ax}, sink::ImageMeta) where {S,T,N,I,D,Ax} = write(s, getfield(sink, :data))
Base.write(s::AxisStream{S,T,N,I,Ax}, sink::AxisArray) where {S,T,N,I,D,Ax} = write(s, getfield(sink, :data))

# Triangle
function Base.write(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Triangle}) where {S,T,N,I,Ax}
    perm = permutation(axisnames(s), (:triangledim, setdiff(axisnames(s), :triangledim)...,))
    write(getfield(s,:data), PermutedDimsArray(mappedview(extracttriangle, sink), perm))
end

# Point
function Base.write(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Point}) where {S,N,I,Ax}
    perm = permutation(axisnames(s), (:pointdim, setdiff(axisnames(s), :pointdim)...,))
    write(getfield(s,:data), PermutedDimsArray(mappedview(extractpoint, A), perm))
end

# Quat
function Base.write(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Quat}) where {S,N,I,Ax}
    perm = permutation(axisnames(s), (:quatdim, setdiff(axisnames(s), :quatdim)...,))
    write(getfield(s,:data), PermutedDimsArray(mappedview(extractquat, A), perm))
end

# SMatrix
function Base.write(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{Quat}) where {S,N,I,Ax}
    perm = permutation(axisnames(s), (:matdim1, :matdim2, setdiff(axisnames(s), :matdim1, :matdim2)...,))
    write(getfield(s,:data), PermutedDimsArray(extractsmat(sink), perm))
end

# SVector
function Base.write(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{SVector}) where {S,N,I,Ax}
    perm = permutation(axisnames(s), (:vecdim, setdiff(axisnames(s), :vecdim)...,))
    write(getfield(s,:data), PermutedDimsArray(extractsvec(sink), perm))
end

# ColorTypes
function Base.write(s::AxisStream{S,<:BitTypes,N,I,Ax}, sink::Array{<:Union{Color3,ColorAlpha}}) where {S,N,I,Ax}
    perm = permutation(axisnames(s), (:colordim, setdiff(axisnames(s), :colordim)...,))
    write(getfield(s,:data), PermutedDimsArray(channelview(sink), perm))
end

# Distributions
function Base.write(s::AxisStream{S,<:BitTypes}, sink::Array{T}) where {S,T<:Distribution}
    perm = permutation(axisnames(s), (:statdim, setdiff(axisnames(s), :statdim)...,))
    write(getfield(s,:data), PermutedDimsArray(extractstat(sink), perm))
end

ext: jl


text:
@enum QueryDetail begin
  LABEL
  SKELETON
  SUMMARY
  FULL
end

const assetGraphVizImg = "http://www.navability.io/wp-content/uploads/2022/03/factor_graph.png"
const assetGeomVizImg = "http://www.navability.io/wp-content/uploads/2022/03/geometric_map.png"

"""
    $TYPEDEF

Helper type for linking to App visualization of a factor graph for user:robot:session.
"""
struct GraphVizApp
  url::String
end

"""
    $TYPEDEF

Helper type for linking to App visualization of geometric map for user:robot:session.
"""
struct MapVizApp
  url::String
end


text:
function calculate_field(obj::GeoMagneticField, arg0::jdouble, arg1::jdouble, arg2::jdouble)
    return jcall(obj, "calculateField", GeoMagneticElements, (jdouble, jdouble, jdouble), arg0, arg1, arg2)
end

function get_decimal_year(::Type{GeoMagneticField}, arg0::jint, arg1::jint, arg2::jint)
    return jcall(GeoMagneticField, "getDecimalYear", jdouble, (jint, jint, jint), arg0, arg1, arg2)
end

function get_epoch(obj::GeoMagneticField)
    return jcall(obj, "getEpoch", jdouble, ())
end

function get_model_name(obj::GeoMagneticField)
    return jcall(obj, "getModelName", JString, ())
end

function supports_time_transform(obj::GeoMagneticField)
    return jcall(obj, "supportsTimeTransform", jboolean, ())
end

function transform_model(obj::GeoMagneticField, arg0::GeoMagneticField, arg1::jdouble)
    return jcall(obj, "transformModel", GeoMagneticField, (GeoMagneticField, jdouble), arg0, arg1)
end

function transform_model(obj::GeoMagneticField, arg0::jdouble)
    return jcall(obj, "transformModel", GeoMagneticField, (jdouble,), arg0)
end

function valid_from(obj::GeoMagneticField)
    return jcall(obj, "validFrom", jdouble, ())
end

function valid_to(obj::GeoMagneticField)
    return jcall(obj, "validTo", jdouble, ())
end


text:
abstract type AbstractTracerInfo end
abstract type AbstractTracerSource end

mutable struct Tracer
    xpos::Float64
    ypos::Float64
end

mutable struct TracerInfo
    Δt::Float64                     # time step
    T::Float64                      # final time
    nT::Int                         # number of time steps
    TracerArr::Vector{Tracer}
end

struct TracerLineSource <: AbstractTracerSource
    LineStart::Vector{Float64}
    LineEnd::Vector{Float64}
    nTracers::Int
end

struct TracerBoxSource <: AbstractTracerSource

end

function GenerateTracers(mesh::FluidMesh,u::Vector{T},v::Vector{T},TI::TracerInfo,file_name::String) where T<:Real
    NTracers = length(TI.TracerArr)

    # Generate time array
    timeArr = 0.0:TI.Δt:TI.T

    # needed to check whether point is inside domain
    xyMesh = NodesToArray(mesh.xy)
    cmMesh = ElementsToArray(mesh.cm)

    xy = zeros(NTracers,2)
    removeArrayVals = []

    # initialize xy array
    for tj=1:NTracers
        xy[tj,1] = TI.TracerArr[tj].xpos
        xy[tj,2] = TI.TracerArr[tj].ypos
    end

    # export vtk save file
    PARTICLES_TO_VTK(xy,1.0,string(file_name,"_0");time=0.0)

    # for each time pointer
    for ti=1:TI.nT
        for tj=1:NTracers
            # check whether any nodes fell outside of domain on last entry
            if !(isInsideDomain(xyMesh,cmMesh,[TI.TracerArr[tj].xpos,TI.TracerArr[tj].ypos]))
                push!(removeArrayVals,tj)
            end
        end

        # check to generate new tracers -- not yet

        # interpolate velocity at tracer positions
        for tj=1:NTracers
            # checks that tracer is still in domain
            if !(tj in removeArrayVals)
                # move tracers by flow field (requires interpolation)
                uval = pointTransform(mesh,u,[TI.TracerArr[tj].xpos,TI.TracerArr[tj].ypos])
                vval = pointTransform(mesh,v,[TI.TracerArr[tj].xpos,TI.TracerArr[tj].ypos])

                TI.TracerArr[tj].xpos += uval*TI.Δt
                TI.TracerArr[tj].ypos += vval*TI.Δt

                xy[tj,1] = TI.TracerArr[tj].xpos
                xy[tj,2] = TI.TracerArr[tj].ypos
            end
        end

        # export vtk save file
        PARTICLES_TO_VTK(xy,1.0,string(file_name,"_",ti);time=ti*TI.Δt)
    end
    nothing
end


text:
#
# A-B Call Center problem
#
using DiscreteEvents, Distributions, Random

mutable struct Caller
    id::Int
    t₁::Float64  # arrival time
    t₂::Float64  # beginning of service time
    t₃::Float64  # end of servive time
end

mutable struct Server
    id::Int
    S::Distribution  # service time distribution
    tbusy::Float64   # cumulative service time
end

function serve(c::Clock, s::Server, input::Channel, output::Vector{Caller}, limit::Int)
    call = take!(input)           # take a call
    call.t₂ = c.time              # record the beginning of service time
    @delay c s.S                  # delay for service time
    call.t₃ = c.time              # record the end of service time
    s.tbusy += call.t₃ - call.t₂  # log service time
    push!(output, call)           # hang up
    call.id ≥ limit && stop!(c)
end

function arrive(c::Clock, input::Channel, count::Vector{Int})
    count[1] += 1
    put!(input, Caller(count[1], c.time, 0.0, 0.0))
end

Random.seed!(123)
const N = 1000
const M_arr = Exponential(2.5)
const M_a   = Exponential(3)
const M_b   = Exponential(4)
count = [0]

clock = Clock()
input = Channel{Caller}(Inf)
output = Caller[]
s1 = Server(1, M_a, 0.0)
s2 = Server(2, M_b, 0.0)
@process serve(clock, s1, input, output, N)
@process serve(clock, s2, input, output, N)
@event arrive(clock,input,count) every M_arr
@run! clock 5000

# using Plots
# wt = [c.t₂ - c.t₁ for c in output]
# plot(wt, title="A-B-Call center ", xlabel="callers", ylabel="waiting time [min]", legend=false)
# savefig("ccenter.png")


text:
function alpha(p::Int)
    a = zeros(p, 1);
    for i = 0:p-1
        for j = 0:i
            for k = i:p-1
                a[i+1] = a[i+1] + (-1.0)^(k-j)/((k+j+1)*
                        factorial(j)*
                        factorial(p-1)*
                        factorial(p-1-k)*
                        factorial(i-j)*
                        factorial(k-i));
            end
        end
    end
    return a;
end

function spline_kernel(t::AbstractArray, p::Int)

    fp = factorial.(p-1:-1:0)
    a = alpha(p).*fp;
    Ut = (repeat(t,p,1).^Vector(p-1:-1:0))./fp;
    Vt = (repeat(t,p,1).^Vector(p:2*p-1)).*a;

    return Ut,Vt;
end

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
abstract type AbstractCostFunction end

struct SigmoidCrossEntropy <: AbstractCostFunction
    name::String
    function SigmoidCrossEntropy(name::String = "sigmimd cross entropy")
        new(name)
    end
end

struct SoftmaxCrossEntropy <: AbstractCostFunction
    name::String
    function SoftmaxCrossEntropy(name::String = "softmax cross entropy")
        new(name)
    end
end

struct SumSquareError <: AbstractCostFunction
    name::String
    function SumSquareError(name::String = "sum square error")
        new(name)
    end
end

function delta(
    cost_fn::AbstractCostFunction,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    # output is the output of final layer
    # [output_1 output_2 ,,, output_N]
    # [target_1 target_2 ,,, target_N]
    if typeof(cost_fn) != GaussianMixtureError
        size(output) == size(targets) || error("delta: size does not match")
    end
    return delta_impl(cost_fn, output, targets)
end

function delta_impl(
    cost_fn::SigmoidCrossEntropy,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    # `output` denotes "logits"
    probs = 1.0 ./ (1.0 .+ exp.(-output))
    return probs - targets
end

function delta_impl(
    cost_fn::SoftmaxCrossEntropy,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    output_max = [maximum(output[:, i]) for i = 1:size(output)[2]]
    output_ = output .- transpose(output_max)
    output_ = exp.(output_)
    probs = output_ ./ sum(output_, dims = 1)

    typeof(probs) == typeof(targets) || error("delta_impl: size does not match")
    return probs - targets
end

function delta_impl(
    cost_fn::SumSquareError,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    return output - targets
end

function cost(
    cost_fn::AbstractCostFunction,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    size(output) == size(targets) || error("cost: size does not match")
    return cost_impl(cost_fn, output, targets)
end

function cost_impl(
    cost_fn::SigmoidCrossEntropy,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    probs = 1.0 ./ (1.0 .+ exp.(-output))
    p = clamp.(probs, 1e-10, 1.0 - 1e-10)

    tmp = -(targets .* log.(p)) - ((1.0 .- targets) .* log.(1.0 .- p))
    return sum(tmp)
end

function cost_impl(
    cost_fn::SoftmaxCrossEntropy,
    ouput::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    output_max = [maximum(output[:, i]) for i = 1:size(output)[2]]
    output_ = output .- transpose(output_max)
    output_ = exp.(output_)
    probs = output_ ./ sum(output_, dims = 1)

    p = clamp.(probs, 1e-10, 1.0 - 1e-10)
    return -sum(targets .* log.(p))
end

function cost_impl(
    cost_fn::SumSquareError,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    diff = output - targets
    return sum(diff .* diff)
end

# not working yet
struct GaussianMixtureError <: AbstractCostFunction
    name::String
    _n_components::Int64
    function GaussianMixtureError(
        n_components::Int64,
        name::String = "gaussian mixture error",
    )
        new(name, n_components)
    end
end

function gaussian(x::Float64, mu::Float64, sigma2::Float64)
    return exp(-(x - mu)^2 / (2 * sigma2)) / sqrt(2.0 * pi * sigma2)
end

function delta_impl(
    cost_fn::GaussianMixtureError,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    # output is the output of final layer
    # [output_1 output_2 ,,, output_N] of size(n_feature, n_samples)
    # [target_1 target_2 ,,, target_N] of size(1, n_samples)
    n_samples = size(targets)[2]
    n_components = cost_fn._n_components
    @assert size(output)[1] == 3 * n_components

    ret = zeros(3 * n_components, n_samples)
    for n = 1:n_samples
        pis = view(output, 1:n_components, n)
        mus = view(output, n_components+1:2*n_components, n)
        sigma2s = view(output, 2*n_components+1:3*n_components, n)

        # (1) convert output to [pis, mus, sigma2s]
        # pis
        max_pi = maximum(pis)
        for k = 1:n_components
            pis[k] = pis[k] - max_pi
        end
        for k = 1:n_components
            pis[k] = exp(pis[k])
        end
        sum_pis = sum(pis)
        for k = 1:n_components
            pis[k] /= sum_pis
        end
        # sigma2s
        for k = 1:n_components
            sigma2s[k] = exp(sigma2s[k])^2
        end

        # (2) compute the graident
        # likelihoods of t_n for each gaussian pi_k N(t_n | mu_k, sigma_k)
        t_n = targets[1, n]
        probs = [pis[k] * gaussian(t_n, mus[k], sigma2s[k]) for k = 1:n_components]
        gammas = probs ./ sum(probs)
        ret[1:n_components, n] = pis - gammas
        ret[n_components+1:2*n_components, n] = gammas .* (mus .- t_n) ./ sigma2s
        ret[2*n_components+1:3*n_components, n] =
            gammas .* (1.0 .- (t_n .- mus) .* (t_n .- mus) ./ sigma2s)
    end

    return ret
end

function cost_impl(
    cost_fn::GaussianMixtureError,
    output::AbstractArray{Float64,2},
    targets::AbstractArray{Float64,2},
)
    n_samples = size(targets)[1]
    n_components = cost_fn._n_components
    @assert size(output)[1] == 3 * n_components

    cost = 0.0
    for n = 1:n_samples
        pis = view(output, 1:n_components, i)
        mus = view(output, n_components+1:2*n_components, i)
        sigma2s = view(output, 2*n_components+1:3*n_components, i)

        # (1) convert output to [pis, mus, sigma2s]
        # pis
        max_pi = maximum(pis)
        for k = 1:n_components
            pis[k] = pis[k] - max_pi
        end
        for k = 1:n_components
            pis[k] = exp(pis[k])
        end
        sum_pis = sum(pis)
        for k = 1:n_components
            pis[k] /= sum_pis
        end
        # sigma2s
        for k = 1:n_components
            sigma2s[k] = exp(sigma2s[k])
        end

        ln_sum = 0.0
        for k = 1:n_components
            ln_sum += pis[k] * gaussian(targets[n], mus[k], sigma2s[k])
        end
        cost -= log(ln_sum)
    end

    return cost
end

max_stars_repo_name: soblin/prml-julia

max_issues_repo_name: soblin/prml-julia

max_forks_repo_name: soblin/prml-julia


text:
function utility_test()
	@testset "Lifetime utility" begin
		w = CollegeStratWorker.make_test_worker(3);
		hV = [1.8, 2.3];
		assetV = [-0.7, 0.9, 3.2];
		workStartAge = 3;

		# On a grid
		utilM = lifetime_utility_grid(w, workStartAge, assetV, hV);
		@test size(utilM) == (length(assetV), length(hV))
		# Increasing in h
		@test all(utilM[:,2] .> utilM[:,1])
		# Increasing in assets
		util2M = lifetime_utility_grid(w, workStartAge, assetV .+ 0.001, hV);
		@test all(util2M .> utilM)

		# Test against state by state
		util3M = similar(utilM);
		for (ih, h) in enumerate(hV)
            for (ia, a) in enumerate(assetV)
                # Must be qualified b/c exported by two packages
				util3M[ia, ih] = CollegeStratWorker.lifetime_utility(w, workStartAge, assetV[ia], hV[ih]);
			end
		end
		@test all(util3M .≈ utilM)

		# Vector
		h2V = fill(hV[1], length(assetV));
		util3V = lifetime_utility_vector(w, workStartAge, assetV, h2V);
		@test all(util3V .≈ utilM[:, 1])

		# MU wealth
		ltIncome = [1.2 2.3; 3.4 4.5];
		mu = CollegeStratWorker.mu_wealth(w, workStartAge, ltIncome);
		@test size(ltIncome) == size(mu)
		@test all(mu .> 0.0)
		# Test by perturbation
		dy = 1e-5;
		util1M = [csw.lifetime_utility(w, workStartAge, lty)  for lty in ltIncome];
		util2M =
			[csw.lifetime_utility(w, workStartAge, lty .+ dy)  for lty in ltIncome];
		muM = (util2M .- util1M) ./ dy;
		@test all(muM .> 0.0)
		@test isapprox(muM, mu, rtol = 1e-3)
	end
end

function lt_util_fct_test()
	@testset "Lifetime Utility Function" begin
		wk = CollegeStratWorker.make_test_worker(3);
		workStartAge = 3;
		nk = 30;
		nh = 20;
		kGridV = LinRange(-10.0, 25.0, nk);
		hGridV = LinRange(1.2, 4.7, nh);
		f = lifetime_utility_function(wk, workStartAge;
			kMin = first(kGridV), kMax = last(kGridV),
			hMin = first(hGridV), hMax = last(hGridV));

		for ik = 2 : 3 : (nk-1)
			for ih = 2 : 4 : (nh-1)
				k = kGridV[ik] + 0.3;
				h = hGridV[ih] - 0.1;
				u1 = f(k, h);
				u2 = csw.lifetime_utility(wk, workStartAge, k, h);
				@test isapprox(u1, u2, rtol = 1e-3);
			end
		end

        dev = interpolation_deviation(wk, workStartAge, f);
        @test dev > 0
        @test dev < 1e-3

        uGrid, kGrid, hGrid = lifetime_utility_function_grid(f);
        uGrid2 = lifetime_utility_grid(wk, workStartAge, kGrid, hGrid);
        dev2 = maximum(abs.(uGrid2 .- uGrid) ./ max.(1.0, abs.(uGrid)));
        @test dev2 < 1e-3
	end
end

@testset "Lifetime utility" begin
	utility_test();
	lt_util_fct_test();
end

# -------------


text:
# A truncated version of Stiefel’s Conjugate Residual method
# cr(A, b, M, atol, rtol, γ, itmax, radius, verbose, linesearch) solves the linear system 'A * x = b' or the least-squares problem :
# 'min ‖b - A * x‖²' within a region of fixed radius.
#
# Marie-Ange Dahito, <marie-ange.dahito@polymtl.ca>
# Montreal, QC, June 2017

export cr

"""A truncated version of Stiefel’s Conjugate Residual method to solve the symmetric linear system Ax=b.
The matrix A must be positive semi-definite.

A preconditioner M may be provided in the form of a linear operator and is
assumed to be symmetric and positive definite.
In a linesearch context, 'linesearch' must be set to 'true'.
"""
function cr(A :: AbstractLinearOperator, b :: AbstractVector{T};
            M :: AbstractLinearOperator=opEye(size(A,1)), atol :: Float64=1.0e-8,
            rtol :: Float64=1.0e-6, γ :: Float64=1.0e-6, itmax :: Int=0,
            radius :: Float64=0.0, verbose :: Bool=false, linesearch :: Bool=false) where T <: Number

  if linesearch && (radius > 0.0)
    error("'linesearch' set to 'true' but radius > 0")
  end
  n = size(b, 1) # size of the problem
  (size(A, 1) == n & size(A, 2) == n) || error("Inconsistent problem size")
  verbose && @printf("CR: system of %d equations in %d variables\n", n, n)

  # Initial state.
  x = zeros(T, n) # initial estimation x = 0
  xNorm = 0.0
  r = M * b # initial residual r = M * (b - Ax) = M * b
  Ar = A * r
  ρ = @kdot(n, r, Ar)
  ρ == 0.0 && return (x, SimpleStats(true, false, [0.0], [], "x = 0 is a zero-residual solution"))
  p = copy(r)
  q = copy(Ar)
  if verbose
    m = 0.0 # quadratic model
  end

  iter = 0
  itmax == 0 && (itmax = 2 * n)

  rNorm = @knrm2(n, r) # ‖r‖
  rNorms = [rNorm] # Values of ‖r‖
  rNorm² = rNorm * rNorm
  pNorm = rNorm
  pNorm² = rNorm²
  pr = rNorm²
  abspr = pr
  pAp = ρ
  abspAp = abs(pAp)
  ArNorm = @knrm2(n, Ar) # ‖Ar‖
  ArNorms = [ArNorm]
  ε = atol + rtol * rNorm
  verbose && @printf("%5s %8s %8s %8s\n", "Iter", "‖x‖", "‖r‖", "quad")
  verbose && @printf("    %d  %8.1e %8.1e %8.1e\n", iter, xNorm, rNorm, m)

  descent = pr > 0.0 # pᵀr > 0 means p is a descent direction
  solved = rNorm ≤ ε
  tired = iter ≥ itmax
  on_boundary = false
  npcurv = false
  status = "unknown"

  while ! (solved || tired)
    if linesearch
      if (pAp ≤ γ * pNorm²) || (ρ ≤ γ * rNorm²)
        npcurv = true
        verbose && @printf("nonpositive curvature detected: pᵀAp = %8.1e and rᵀAr = %8.1e\n", pAp, ρ)
        stats = SimpleStats(solved, false, rNorms, ArNorms, "nonpositive curvature")
        iter == 0 && return (b, stats)
        return (x, stats)
      end
    elseif pAp ≤ 0.0 && radius == 0.0
      error("Indefinite system and no trust region")
    end
    Mq = M * q

    if radius > 0.0
      verbose && @printf("radius = %8.1e > 0 and ‖x‖ = %8.1e\n", radius, xNorm)
      # find t1 > 0 and t2 < 0 such that ‖x + ti * p‖² = radius²  (i = 1, 2)
      xNorm² = xNorm * xNorm
      t = to_boundary(x, p, radius; flip = false, xNorm2 = xNorm², dNorm2 = pNorm²)
      t1 = maximum(t) # > 0
      t2 = minimum(t) # < 0
      tr = maximum(to_boundary(x, r, radius; flip = false, xNorm2 = xNorm², dNorm2 = rNorm²))
      verbose && @printf("t1 = %8.1e, t2 = %8.1e and tr = %8.1e\n", t1, t2, tr)

      if abspAp ≤ γ * pNorm * @knrm2(n, q) # pᵀAp ≃ 0
        npcurv = true # nonpositive curvature
        verbose && @printf("pᵀAp = %8.1e ≃ 0\n", pAp)
        if abspr ≤ γ * pNorm * rNorm # pᵀr ≃ 0
          verbose && @printf("pᵀr = %8.1e ≃ 0, redefining p := r\n", pr)
          p = r # - ∇q(x)
          q = Ar
          # q(x + αr) = q(x) - α ‖r‖² + ½ α² rᵀAr
          # 1) if rᵀAr > 0, the quadratic decreases from α = 0 to α = ‖r‖² / rᵀAr
          # 2) if rᵀAr ≤ 0, the quadratic decreases to -∞ in the direction r
          if ρ > 0.0 # case 1
            verbose && @printf("quadratic is convex in direction r, curv = %8.1e\n", ρ)
            α = min(tr, rNorm² / ρ)
          else # case 2
            verbose && @printf("r is a direction of nonpositive curvature: %8.1e\n", ρ)
            α = tr
          end
        else
          # q_p = q(x + α_p * p) - q(x) = -α_p * rᵀp + ½ (α_p)² * pᵀAp
          # q_r = q(x + α_r * r) - q(x) = -α_r * ‖r‖² + ½ (α_r)² * rᵀAr
          # Δ = q_p - q_r. If Δ > 0, r is followed, else p is followed
          α = descent ? t1 : t2
          ρ > 0.0 && (tr = min(tr, rNorm² / ρ))
          Δ = -α * pr + tr * rNorm² - 0.5 * (tr)^2 * ρ # as pᵀAp = 0
          if Δ > 0.0 # direction r engenders a better decrease
            verbose && @printf("direction r engenders a bigger decrease. q_p - q_r = %8.1e > 0\n", Δ)
            verbose && @printf("redefining p := r\n")
            p = r
            q = Ar
            α = tr
          else
            verbose && @printf("direction p engenders an equal or a bigger decrease. q_p - q_r = %8.1e ≤ 0\n", Δ)
          end
        end

      elseif pAp > 0.0 && ρ > 0.0 # no negative curvature
        verbose && @printf("positive curvatures along p and r. pᵀAp = %8.1e and rᵀAr = %8.1e\n", pAp, ρ)
        α = ρ / @kdot(n, q, Mq)
        if α ≥ t1
          α = t1
          on_boundary = true
        end

      elseif pAp > 0.0 && ρ < 0.0
        npcurv = true
        verbose && @printf("pᵀAp = %8.1e > 0 and rᵀAr = %8.1e < 0\n", pAp, ρ)
        # q_p is minimal for α_p = rᵀp / pᵀAp
        α = descent ?  min(t1, pr / pAp) : max(t2, pr / pAp)
        Δ = -α * pr + tr * rNorm² + 0.5 * (α^2 * pAp - (tr)^2 * ρ)
        if Δ > 0.0
          verbose && @printf("direction r engenders a bigger decrease. q_p - q_r = %8.1e > 0\n", Δ)
          verbose && @printf("redefining p := r\n")
          p = r
          q = Ar
          α = tr
        else
          verbose && @printf("direction p engenders an equal or a bigger decrease. q_p - q_r = %8.1e ≤ 0\n", Δ)
        end

      elseif pAp < 0.0 && ρ > 0.0
        npcurv = true
        verbose && @printf("pᵀAp = %8.1e < 0 and rᵀAr = %8.1e > 0\n", pAp, ρ)
        α = descent ? t1 : t2
        tr = min(tr, rNorm² / ρ)
        Δ = -α * pr + tr * rNorm² + 0.5 * (α^2 * pAp - (tr)^2 * ρ)
        if Δ > 0.0
          verbose && @printf("direction r engenders a bigger decrease. q_p - q_r = %8.1e > 0\n", Δ)
          verbose && @printf("redefining p := r\n")
          p = r
          q = Ar
          α = tr
        else
          verbose && @printf("direction p engenders an equal or a bigger decrease. q_p - q_r = %8.1e ≤ 0\n", Δ)
        end

      elseif pAp < 0.0 && ρ < 0.0
        npcurv = true
        verbose && @printf("negative curvatures along p and r. pᵀAp = %8.1e and rᵀAr = %8.1e\n", pAp, ρ)
        α = descent ? t1 : t2
        Δ = -α * pr + tr * rNorm² + 0.5 * (α^2 * pAp - (tr)^2 * ρ)
        if Δ > 0.0
          verbose && @printf("direction r engenders a bigger decrease. q_p - q_r = %8.1e > 0\n", Δ)
          verbose && @printf("redefining p := r\n")
          p = r
          q = Ar
          α = tr
        else
          verbose && @printf("direction p engenders an equal or a bigger decrease. q_p - q_r = %8.1e ≤ 0\n", Δ)
        end
      end

    elseif radius == 0.0
      α = ρ / @kdot(n, q, Mq) # step
    end

    @kaxpy!(n,  α, p, x)
    xNorm = @knrm2(n, x)
    xNorm ≈ radius && (on_boundary = true)
    @kaxpy!(n, -α, Mq, r) # residual
    rNorm² = abs(rNorm² - α * ρ)
    rNorm = sqrt(rNorm²)
    push!(rNorms, rNorm)
    Ar = A * r
    ArNorm = @knrm2(n, Ar)
    push!(ArNorms, ArNorm)

    iter = iter + 1
    if verbose
      m = m - α * pr + 0.5 * α^2 * pAp
      @printf("    %d  %8.1e %8.1e %8.1e\n", iter, xNorm, rNorm, m)
    end

    solved = (rNorm <= ε) || npcurv || on_boundary
    tired = iter >= itmax

    (solved || tired) && continue
    ρbar = ρ
    ρ = @kdot(n, r, Ar)
    β = ρ / ρbar # step for the direction computation
    @kaxpby!(n, 1.0, r, β, p)
    @kaxpby!(n, 1.0, Ar, β, q)

    pNorm² = rNorm² + 2 * β * pr - 2 * β * α * pAp + β^2 * pNorm²
    pNorm = sqrt(pNorm²)
    pr = rNorm² + β * pr - β * α * pAp # pᵀr
    abspr = abs(pr)
    pAp = ρ + β^2 * pAp # pᵀq
    abspAp = abs(pAp)
    descent = pr > 0.0

  end
  verbose && @printf("\n")

  status = npcurv ? "nonpositive curvature" : (on_boundary ? "on trust-region boundary" : (tired ? "maximum number of iterations exceeded" : "solution good enough given atol and rtol"))
  stats = SimpleStats(solved, false, rNorms, ArNorms, status)
  return (x, stats)
end


text:
const World = Dict{Symbol,Object}

"""
Simulate a world of objects for some number of iterations and fixed timestep.

This repeatedly applies the hydroelastic force model to each pair of objects in
the scene, returning a vector of how the objects change over time.

The kinematics are implemented using a symplectic Euler method, which is more
numerically stable than forward Euler integration.
"""
function simulate(
    world::World;
    dt::Float64 = 1e-3,
    num_steps::Int64 = 1000,
    forces::Function = (a, α) -> (),
)::Dict{Symbol,Vector{Object}}
    ids = collect(keys(world))
    objects = Dict(id => [world[id]] for id in ids)

    for step = 1:num_steps
        # Hydroelastic computations for each pair of objects
        a = Dict(id => @SVector zeros(3) for id in ids)
        α = Dict(id => @SVector zeros(3) for id in ids)
        for i = 1:length(ids)
            for j = i+1:length(ids)
                u, v = ids[i], ids[j]
                results = compute_force(last(objects[u]), last(objects[v]))
                a[u] += 100.0 * results.F_AB
                a[v] += 100.0 * results.F_BA
                α[u] += 100.0 * results.τ_AB
                α[v] += 100.0 * results.τ_BA
            end
        end

        # Apply custom user-specified forces
        forces(a, α)

        # Semi-implicit Euler integration
        for id in ids
            obj = last(objects[id])
            obj = Object(obj.mesh, obj.pose, obj.v .+ dt .* a[id], obj.ω .+ dt .* α[id])
            obj = translate(obj, dt * obj.v)
            obj = changeRotation(obj, dt * obj.ω)
            push!(objects[id], obj)
        end
    end

    objects
end

export World, simulate

max_stars_repo_name: ekzhang/hydroelastics

max_issues_repo_name: ekzhang/hydroelastics

max_forks_repo_name: ekzhang/hydroelastics


text:
# #
# Study of the spectral stability of the Serre-Green-Naghdi equations.
# #
export Spectrum,figspecCW,figspecSW,Experiments
using WaterWaves1D,LinearAlgebra,FFTW,ProgressMeter,Plots;

"""
    Spectrum(param,η,u,c,ν;kwargs)

Compute the spectrum about a specified travelling solution of the
the linearized Green-Naghdi model.

# Argument
- `param` should be a NamedTuple
	e.g. `param = ( μ  = 0.1, ϵ  = 1, N  = 2^9, L  = 10*π)`

- `(η,u)` are real vectors of the aforementioned solution where
where

    - `η` is the surface deformation;
    - `u` corresponds to the layer-averaged velocity.
- 'c' is the velocity of the solution,
- 'ν' is a Floquet coefficient.

## Keywords
- `SGN :: Bool`: if `true` computes the Serre-Green-Naghdi (SGN) instead of Whitham-Green-Naghdi (WGN) system (default is `false`);

# Return values
an Eigen factorization object `L` which
contains the eigenvalues in `L.values` and the
eigenvectors in the columns of the matrix `L.vectors`.
(The `k`th eigenvector can be obtained from the slice `F.vectors[:, k]``.)

"""
function Spectrum(param,η,u,c,ν)

	μ 	= param.μ
	ϵ 	= param.ϵ
	mesh = Mesh(param)
	k = mesh.k
	x 	= mesh.x
	x₀ = mesh.x[1]
	∂ₓ	=  1im * mesh.k
	F₁ = 1 ./(1 .+ μ/3*k.^2)
	F₀ = sqrt(μ)*∂ₓ
	FFT = exp.(-1im*k*(x.-x₀)');
	IFFT = exp.(1im*k*(x.-x₀)')/length(x);
	Id = Diagonal(ones(size(x)));

	h = 1 .+ ϵ*η
	DxF(v) = ifft(F₀ .* fft(v))
	v= u - 1/3 ./h .* (DxF(h.^3 .*DxF(u)))

	Dx = IFFT*Diagonal(∂ₓ)*FFT + ν*Id#Diagonal(m.x./sqrt.(1 .+ m.x.^2))
	DxF₀ = IFFT*Diagonal(F₀)*FFT + ν*Id#Diagonal(m.x./sqrt.(1 .+ m.x.^2))

	A11 = -Dx*Diagonal(ϵ*u)
	A12 = -Dx*Diagonal(h)
	A13 = zeros(size(A11))
	A21 = -Dx* ( Id - Diagonal( h.*(DxF(ϵ*u).^2) ) )
	A22 = -Dx*( ϵ*Diagonal( v - u ) - Diagonal( (h.^2).*DxF(ϵ*u)) *DxF₀ )
	A23 = -Dx*Diagonal(ϵ*u)

	B1 = -Diagonal(1 ./h)*DxF₀*Diagonal( (h.^2).*DxF(ϵ*u) ) + 1/3*Diagonal(DxF(h.^3 .* DxF(ϵ*u))./h.^2)
	B2 = Id - 1/3*Diagonal(1 ./h)*DxF₀*Diagonal(h.^3)*DxF₀

	C2 = inv(B2)
	C1 = -C2*B1

	M11 = A11+A12*C1 + c*Dx
	M12 = A13+A12*C2
	M21 = A21+A22*C1
	M22 = A23+A22*C2 + c*Dx

	return eigen([ M11 M12  ; M21 M22 ])
end

function ell(a₀,a₁,k)
	h₀ = a₀
	h₂ = a₁ + a₀
	h₁ = (1-k^2)*a₁ + a₀
	return (h₀=h₀,h₁=h₁,h₂=h₂)
end

"""
	figspecCW(p;N,P,SGN)

Spectrum of the linearized SGN equations about a cnoidal wave.

- Argument `p` is a NamedTuple and should contain `h₀<h₁<h₂` definig the cnoidal wave.
- `N` is the number of collocation points.
- `P` is the number of periods in the mesh (default is `P=1`).
- `M` is the number of (equally spaced) Floquet points conputed.
- Solves SGN if `SGN` is `true`, WGN otherwise (not implemented yet).

"""
function figspecCW(p;N=2^6,P=1,M=2^6,SGN=true::Bool)
	param = merge( p,( μ  = 1, ϵ  = 1, N  = N, P = P ))

	(η,u,v,mesh,parc)=CnoidalWaveSerreGreenNaghdi(param;P=P)
	@show parc

	plt = plot(layout=(2,1))
	plot!(plt[1,1], mesh.x, [η u],label=["η" "u"])
	plot!(plt[2,1], fftshift(mesh.k),
	  [abs.(fftshift(fft(η))).+eps() abs.(fftshift(fft(u))).+eps()],
	  label=["|η̂|" "|û|"],yscale=:log10)
	display(plt)

	param2=merge(param,(L=P*parc.λ,))

	σ = Spectrum(param2,η,u,parc.c,1im*0.).values
	val = σ#[Int(7*end/8):end]
	@showprogress 1 for ν=-π/parc.λ/P/2:π/parc.λ/P/M:π/parc.λ/P/2
		σ = Spectrum(param2,η,u,parc.c,1im*ν).values
		val=[val;σ]
		# uncomment the following lines to stop when a unstable mode has been found
		# if maximum(real.(σ))>0.1
		# 	return Spectrum(param,η,u,parc.c,1im*ν),merge(parc,(ν=ν,)),η,u,v,mesh
		# 	break
		# end
	end
	#plt = scatter(values,label="")
	#display(plt)
	return val,parc,η,u,v,mesh

end

"""
	figspecSW(c;N,L,SGN)

Spectrum of the linearized SGN equations about a cnoidal wave.

- Argument `c` is the velocity definig the solitary wave.
- `N` is the number of collocation points.
- `L` is half-length of the mesh.
- `M` is the number of (equally spaced) Floquet points conputed.
- Solves SGN if `SGN` is `true`, WGN otherwise.

"""
function figspecSW(c;N=2^6,L=10*π,M=2^6,SGN=true::Bool)
	param = ( c=c, μ  = 1, ϵ  = 1, N  = N, L = L )

	if SGN == true
		(η,u,v,mesh)=SolitaryWaveSerreGreenNaghdi(param)
	else
		(η,u,v,mesh)=SolitaryWaveWhithamGreenNaghdi(param;method=3,verbose=true)
	end
	plt = plot(layout=(2,1))
	plot!(plt[1,1], mesh.x, [η u],label=["η" "u"])
	plot!(plt[2,1], fftshift(mesh.k),
	  [abs.(fftshift(fft(η))).+eps() abs.(fftshift(fft(u))).+eps()],
	  label=["|η̂|" "|û|"],yscale=:log10)
	display(plt)

	σ = Spectrum(param,η,u,c,1im*0.).values
	val = σ#[Int(7*end/8):end]
	@showprogress 1 for ν=-π/L/2:π/L/M:π/L/2
		σ = Spectrum(param,η,u,c,1im*ν).values
		val=[val;σ]
		# uncomment the following lines to stop when a unstable mode has been found
		# if maximum(real.(σ))>0.1
		# 	return Spectrum(param,η,u,c,1im*ν),η,u,v,mesh
		# 	break
		# end
	end
	#plt = scatter(values,label="")
	#display(plt)
	return val,η,u,v,mesh

end

"""
	Experiments(exp)

Some numerical experiments, plotting the spectrum of the linearized SGN equations about various travelling waves
- `exp==1`: a cnoidal wave with very small amplitude
- `exp==2`: a cnoidal wave with very large period
- `exp==3`: several cnoidal waves with increasing periods
- `exp==4`: a solitary wave

Return the plot.
"""
function Experiments(exp)

if exp == 1 # experiment 1: very small amplitude
	σ,param,η,u,v,mesh=figspecCW(ell(0.3,0.005,0.75);P=1,N=2^5)
	plt=scatter(σ)
	ylims!(0,0.1)
	σ,param,η,u,v,mesh=figspecCW(ell(0.3,0.005,0.75);P=1,N=2^5,M=2^19) # on peut limiter les Floquet dans (-π/parc.λ/P/20,π/parc.λ/P/20)
	scatter(σ[abs.(real.(σ)) .> 1e-10])
	#ylims!(0.03158,0.0316)
	ylims!(0.0113725,0.011375)

elseif exp == 2 # experiment 2: very large period
	σ,param,η,u,v,mesh=figspecCW(ell(0.3,0.1,0.999999999999999);P=1,N=2^8,M=2^8)
	plt=scatter(σ)
	ylims!(0.2,0.25)
	ylims!(0.23,0.235)

elseif exp == 3 # experiment 3 : growing period
	v1,param,η,u,v,mesh=figspecCW((h₀=0.99,h₁=1,h₂=11);P=1,N=2^7)
	v2,param,η,u,v,mesh=figspecCW((h₀=0.999999,h₁=1,h₂=11);P=1,N=2^7)
	v3,param,η,u,v,mesh=figspecCW((h₀=0.9999999999,h₁=1,h₂=11);P=1,N=2^7)
	v4,param,η,u,v,mesh=figspecCW((h₀=0.99999999999999,h₁=1,h₂=11);P=1,N=2^7)

	plt=scatter(v1[abs.(v1).<15],label="a=10,L=6")
	scatter!(plt,v2[abs.(v2).<15],label="a=10,L=12")
	scatter!(plt,v3[abs.(v3).<15],label="a=10,L=18")
	scatter!(plt,v4[abs.(v4).<15],label="a=10,L=22")

elseif exp == 4# experiment 4 : solitary wave
	σ,η,u,v,mesh=figspecSW(3;L=10*π,N=2^8)
	plt=scatter(σ)
end
display(plt)
return plt
end


text:
module TriangleIntersect

import Base: -,*,+,/, intersect

export Point, Triangle, Ray, Intersection, no_intersection

struct Point
    x::Float64
    y::Float64
    z::Float64
end

*(p::Point, n::Number) = Point(p.x*n, p.y*n, p.z*n)
*(n::Number, p::Point) = p*n
*(p1::Point, p2::Point) = p1.x*p2.x+p1.y*p2.y+p1.z*p2.z
-(p1::Point, p2::Point) = Point(p1.x-p2.x, p1.y-p2.y, p1.z-p2.z)
+(p1::Point, p2::Point) = Point(p1.x+p2.x, p1.y+p2.y, p1.z+p2.z)
cross(p1::Point, p2::Point) = Point(p1.y*p2.z-p1.z*p2.y, -p1.x*p2.z+p1.z*p2.x, p1.x*p2.y-p1.y*p2.x)

/(p::Point, n::Number) = Point(p.x/n, p.y/n, p.z/n)
unitize(p::Point) = p/(p*p)

struct Triangle
    a::Point
    v1::Point
    v2::Point
    normal::Point
    v1v1::Float64
    v2v2::Float64
    v1v2::Float64
    denom::Float64
    function Triangle(a::Point, b::Point, c::Point)
        v1 = b-a
        v2 = c-a
        normal = cross(v1, v2)
        v1v1 = v1*v1
        v2v2 = v2*v2
        v1v2 = v1*v2
        denom = v1v2*v1v2 - v1v1*v2v2
        new(a, v1, v2, normal, v1v1, v2v2, v1v2, denom)
    end
end

struct Ray
    origin::Point
    direction::Point
    Ray(a,b) = new(a,unitize(b))
end

struct Intersection
    ip::Point # intersecting point
    id::Float64 # intersecting distance
    is_intersection::Bool
end

const no_intersection = Intersection(Point(0,0,0), 0.0, false)

function intersect(r::Ray, t::Triangle)
    denom = t.normal*r.direction
    denom == 0 && return no_intersection
    ri = t.normal*(t.a - r.origin) / denom
    ri <= 0 && return no_intersection
    plane_intersection =  ri * r.direction + r.origin
    w = plane_intersection - t.a
    wv1 = w*t.v1
    wv2 = w*t.v2
    s_intersection = (t.v1v2*wv2 - t.v2v2*wv1) / t.denom
    s_intersection <= 0 && return no_intersection
    s_intersection >= 1 && return no_intersection
    t_intersection = (t.v1v2*wv1 - t.v1v1*wv2) / t.denom
    t_intersection <= 0 && return no_intersection
    t_intersection >= 1 && return no_intersection
    s_intersection + t_intersection >= 1 && return no_intersection
    Intersection(t.a + s_intersection*t.v1+t_intersection*t.v2, ri, true)
end

end # module


text:
#=

Implementation of commutator-free exponential time propagators, as
described in

  Alvermann, A., Fehske, H., & Littlewood, P. B. (2012). Numerical
  time propagation of quantum systems in radiation fields. New
  J. Phys., 14(10),
  105008. http://dx.doi.org/10.1088/1367-2630/14/10/105008

=#
struct CFET4BfCPropagator{E<:Exponentiator} <: MagnusPropagator
    B::LinearMap
    f::Function
    C::LinearMap
    a::Number
    Exp::E
end
CFET4BfCPropagator(B::LinearMap,
                   f::Function,
                   C::LinearMap,
                   Exp::E) where E<:Exponentiator =
                       CFET4BfCPropagator(B, f, C, 1, Exp)

function (p::CFET4BfCPropagator{E})(t::T, τ::T,
                                    v::AbstractVector, w::AbstractVector) where {E<:Exponentiator,
                                                                                T<:AbstractFloat}
    h1 = T(37/66 - 400/957*sqrt(5/3))
    h2 = T(-4/33)
    h3 = T(37/66 + 400/957*sqrt(5/3))
    h4 = T(-11/162)
    h5 = T(92/81)

    ff1 = p.f(t + T(1/2-sqrt(3/20))*τ)
    ff2 = p.f(t + τ/2)
    ff3 = p.f(t + T(1/2+sqrt(3/20))*τ)

    f1 = h1*ff1 + h2*ff2 + h3*ff3
    f2 = h4*ff1 + h5*ff2 + h4*ff3
    f3 = h3*ff1 + h2*ff2 + h1*ff3

    τ1 = 11τ/40
    τ2 = 9τ/20

    p.Exp(p.B + f3*p.C, p.a*τ1, v, w)
    p.Exp(p.B + f2*p.C, p.a*τ2, w, w)
    p.Exp(p.B + f1*p.C, p.a*τ1, w, w)
end

export CFET4BfCPropagator


text:
"""
**Laplace likelihood**

Laplace likelihood for regression: ``\\frac{1}{2\\beta}\\exp\\left(-\\frac{|y-f|}{\\beta}\\right)``
see [wiki page](https://en.wikipedia.org/wiki/Laplace_distribution)

```julia
LaplaceLikelihood(β::T=1.0)  #  Laplace likelihood with scale β
```

---

For the analytical solution, it is augmented via:
```math
p(y|f,\\omega) = \\mathcal{N}(y|f,\\omega^{-1})
```
where ``\\omega \\sim \\text{Exp}\\left(\\omega \\mid \\frac{1}{2 \\beta^2}\\right)``, and Exp is the [Exponential distribution](https://en.wikipedia.org/wiki/Exponential_distribution)
We approximate ``q(\\omega) = \\mathcal{GIG}\\left(\\omega \\mid a,b,p\\right)
"""
struct LaplaceLikelihood{T<:Real} <: RegressionLikelihood{T}
    β::LatentArray{T}
    a::LatentArray{T}
    p::LatentArray{T}
    b::LatentArray{Vector{T}}
    θ::LatentArray{Vector{T}} #Expected value of ω
    function LaplaceLikelihood{T}(β::AbstractVector{T}) where {T<:Real}
        new{T}(β,β.^-2,0.5*ones(size(β)))
    end
    function LaplaceLikelihood{T}(β::AbstractVector{T},b::AbstractVector{<:AbstractVector{T}},θ::AbstractVector{<:AbstractVector{T}}) where {T<:Real}
        new{T}(β,β.^(-2),0.5*ones(size(β)),b,θ)
    end
end

function LaplaceLikelihood(β::T=1.0) where {T<:Real}
    LaplaceLikelihood{T}([β])
end

function init_likelihood(likelihood::LaplaceLikelihood{T},inference::Inference{T},nLatent::Int,nSamplesUsed::Int) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        LaplaceLikelihood{T}(
        likelihood.β,
        [abs2.(T.(rand(T,nSamplesUsed))) for _ in 1:nLatent],
        [zeros(T,nSamplesUsed) for _ in 1:nLatent])
    else
        LaplaceLikelihood{T}([likelihood.β for _ in 1:nLatent])
    end
end

function pdf(l::LaplaceLikelihood,y::Real,f::Real)
    pdf(Laplace(f,l.β[1]),y) #WARNING multioutput invalid
end

function Base.show(io::IO,model::LaplaceLikelihood{T}) where T
    print(io,"Laplace likelihood")
end

function compute_proba(l::LaplaceLikelihood{T},μ::AbstractVector{T},σ²::AbstractVector{T}) where {T<:Real}
    N = length(μ)
    nSamples = 2000
    μ_pred = zeros(T,N)
    σ²_pred = zeros(T,N)
    temp_array = zeros(T,nSamples)
    for i in 1:N
        # e = expectation(Normal(μ[i],sqrt(β²[i])))
        # μ_pred[i] = μ[i]
        #
        # β²_pred[i] = e(x->pdf(LocationScale(x,1.0,st))^2) - e(x->pdf(LocationScale(x,1.0,st)))^2
        if σ²[i] <= 1e-3
            for j in 1:nSamples
                temp_array[j] = rand(Laplace(μ[i],l.β[1])) #WARNING Multiouput invalid
            end
        else
            d = Normal(μ[i],sqrt(σ²[i]))
            for j in 1:nSamples
                temp_array[j] = rand(Laplace(rand(d),l.β[1])) #WARNING multioutput invalid
            end
        end
        μ_pred[i] = μ[i];
        σ²_pred[i] = cov(temp_array)
    end
    return μ_pred,σ²_pred
end

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

function local_updates!(model::VGP{<:LaplaceLikelihood,<:AnalyticVI})
    model.likelihood.b .= broadcast((Σ,μ,y)->(Σ+abs2.(μ-y)),diag.(model.Σ),model.μ,model.y)
    model.likelihood.θ .= broadcast((a,b)->sqrt(a)./sqrt.(b),model.likelihood.a,model.likelihood.b)
end

function local_updates!(model::SVGP{<:LaplaceLikelihood,<:AnalyticVI})
    model.likelihood.b .= broadcast((K̃,κ,Σ,μ,y)->(K̃ + opt_diag(κ*Σ,κ) + abs2.(κ*μ-y[model.inference.MBIndices])),model.K̃,model.κ,model.Σ,model.μ,model.y)
    model.likelihood.θ .= broadcast((a,b)->sqrt(a)./sqrt.(b),model.likelihood.a,model.likelihood.b)
end

function sample_local!(model::VGP{<:LaplaceLikelihood,<:GibbsSampling})
    model.likelihood.ω .= NaN
    return nothing
end

""" Return the gradient of the expectation for latent GP `index` """
function cond_mean(model::VGP{<:LaplaceLikelihood,<:AnalyticVI},index::Integer)
    return model.likelihood.θ[index].*model.y[index]
end

function ∇μ(model::VGP{<:LaplaceLikelihood,<:AnalyticVI})
    return hadamard.(model.likelihood.θ,model.y)
end

""" Return the gradient of the expectation for latent GP `index` """
function cond_mean(model::SVGP{<:LaplaceLikelihood,<:AnalyticVI},index::Integer)
    return model.likelihood.θ[index].*model.y[index][model.inference.MBIndices]
end

function ∇μ(model::SVGP{<:LaplaceLikelihood,<:AnalyticVI})
    return hadamard.(model.likelihood.θ,getindex.(model.y,[model.inference.MBIndices]))
end

function ∇Σ(model::AbstractGP{<:LaplaceLikelihood,<:AnalyticVI})
    return model.likelihood.θ
end

function ELBO(model::AbstractGP{<:LaplaceLikelihood,<:AnalyticVI})
    return expecLogLikelihood(model) - GIGExpKL(model) - GaussianKL(model)
end

function expecLogLikelihood(model::VGP{LaplaceLikelihood{T},AnalyticVI{T}}) where T
    tot = -0.5*model.nLatent*model.nSample*log(twoπ)
    tot += 0.5.*sum(broadcast(θ->sum(log.(θ)),model.likelihood.θ))
    tot += -0.5.*sum(broadcast((θ,Σ,μ,y)->dot(θ,(Σ+abs2.(μ)-2.0*μ.*y-abs2.(y))),model.likelihood.θ,diag.(model.Σ),model.μ,model.y))
    return tot
end

function expecLogLikelihood(model::SVGP{LaplaceLikelihood{T},AnalyticVI{T}}) where T
    tot = -0.5*model.nLatent*model.inference.nSamplesUsed*log(twoπ)
    tot += 0.5.*sum(broadcast(θ->sum(log.(θ)),model.likelihood.θ))
    tot += -0.5.*sum(broadcast((θ,K̃,κ,Σ,μ,y)->dot(θ,(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)-2.0*(κ*μ).*y[model.inference.MBIndices]-abs2.(y[model.inference.MBIndices]))),model.likelihood.θ,model.K̃,model.κ,model.Σ,model.μ,model.y))
    return model.inference.ρ*tot
end

function GIGExpKL(model::AbstractGP{<:LaplaceLikelihood})
    GIGEntropy(model)-expecExponentialGIG(model)
end

function GIGEntropy(model::AbstractGP{<:LaplaceLikelihood})
    model.inference.ρ*sum(broadcast(GIGEntropy,model.likelihood.a,model.likelihood.b,model.likelihood.p))
end

function expecExponentialGIG(model::AbstractGP{<:LaplaceLikelihood})
    sum(broadcast((β,a,b)->sum(-log(2*β^2).-0.5*(a.*sqrt.(b)+b.*sqrt(a))./(a.*b*β^2)),model.likelihood.β,model.likelihood.a,model.likelihood.b))
end

function gradpdf(::LaplaceLikelihood,y::Int,f::T) where {T<:Real}
    @error "Not implemented yet"
end

function hessiandiagpdf(::LaplaceLikelihood,y::Int,f::T) where {T<:Real}
    @error "Not implemented yet"
end


text:
#= Julia syntax highlighting test.

This file derives from https://gist.github.com/Wilfred/f1aca44c61ed6e1df603
whose author is [@Wilfred](https://github.com/Wilfred). @Wilfred has put it in
the public domain:
  https://gist.github.com/Wilfred/f1aca44c61ed6e1df603#gistcomment-2948526

Changes from the original are governed by the license of the repository in
which this file is found.

This file is designed to test various corner cases of Julia
syntax highlighting.
=#

## Simple function definitions.
# Expected: `function` should be highlighted, as should `foo_bar!`.
function foo_bar!(x,y)
    x + y + 1
end

# Expected: foo_bar!` should be highlighted.
foo_bar!(x,y) = x + y

# Expected: `foo` should be highlighted.
Baz.foo(x) = 1

# Expected: `foo` should be highlighted.
foo(x::(Int,)) = 1

# Expected: `foo` should be highlighted.
foo(x, y=length(x))

## Function definitions in namespaces.
# Expected: `bar` should be highlighted.
function Foo.bar(x, y)
    x + 1
end

## Function definitions with type variables.
# Expected: `elsize` should be highlighted.
elsize{T}(::AbstractArray{T}) = sizeof(T)

function elsize{T}(::AbstractArray{T})
    sizeof(T)
end

## Nested brackets in function definitions.
# Expected: `cell` should be highlighted.
cell(dims::(Integer...)) = Array(Any, convert((Int...), dims))

# TODO: find an example with a nested type expression.

## Macro usage
# Expected: `@hello_world!` should be highlighted.
@hello_world! foo

## Builtin functions.
# Expected: `throw`, `error` and `super` should not be highlighted. There are
# too many built-in functions for this to be useful.
# https://github.com/JuliaLang/julia/commit/134867c69096fcb52afa5d5a7433892b5127e981
# https://github.com/JuliaLang/julia/pull/7963#issuecomment-52586261
throw(foo)
error("foo", bar, "baz")
super(Int)

## Strings
# Expected: highlight the string.
x = "foo \"bar\" baz"

# Expected: highlight the whole string.
x = "foo
bar"

# Expected: highlight the whole triple-quoted string.
x = """hello "world" foobar"""
y = """foo\\"""
z = """bar\""""

# Expected: highlight `$user`
x = "hello $user"

# Expected: don't highlight `$user`
x = "hello \$user"

# Expected: treat r as part of the string, so `r"a"` is highlighted.
x = r"0.1"

# Expected: treat ismx as part of the string, so `r"a"ismx` is highlighted.
x = r"a"ismx

# Expected: highlight `r"""a "b" c"""`
x = r"""a "b" c"""

# Expected: treat v as part of the string, so `v"0.1"` is highlighted.
x = v"0.1"

# Expected: treat b as part of the string, so `b"a"` is highlighted.
x = b"a"

# Bonus points:
# Expected: highlight the interpolation brackets `$(` and `)`
x = "hello $(user * user)"

# Bonus points:
# Expected: highlight regexp metacharacters `[` and `]`
x = r"[abc]"

# Bonus points:
# Expected: highlight escape sequences `\xff` and `\u2200`
x = b"DATA\xff\u2200"

## Characters
# Expected: highlight the character.
x = 'a'
y = '\u0'
z = '\U10ffff'
a = ' '
b = '"'

# Expected: don't highlight, as ' is an operator here, not a character delimiter.
a = b' + c'

# Bonus points:
# Expected: don't highlight the character
x = 'way too long so not a character'

## Comments
# Expected: highlight `# foo`
# foo

# Expected: highlight `#= foo\n bar =#`
#= foo
bar =#

# Expected: highlight `#= #= =# =#` (comments can nest).
#= #= =# =#

## Type declarations

# Expected highlight `Foo` and `Bar`
type Foo
    x::Bar
end

# Expected highlight `Foo` and `Bar`
immutable Foo
    x::Bar
end

# Expected: highlight `Foo` and `Bar`
abstract Foo <: Bar

# Expected: don't highlight x or y
x <: y

## Type annotations

# Expected: highlight `FooBar`
f(x::FooBar) = x

# Expected: highlight `Int8`
function foo()
    local x::Int8 = 5
    x
end

# Expected: highlight `T` and `Number`
same_type_numeric{T<:Number}(x::T, y::T) = true

## Variable delcarations

# Expected: highlight `x` and `y`
global x = "foo, bar = 2", y = 3

# Expected: highlight `x` and `y`
global x = foo(a, b), y = 3

# Expected: highlight `y`
const y = "hello world"

# Expected: highlight `x` and `y`
function foo()
    local x = f(1, 2), y = f(3, 4)
    x + y
end

# Expected: highlight `x` and `y`
let x = f(1, 2), y = f(3, 4)
    x + y
end

## Colons and end

# Expected: highlight `:foo`, `:end` and `:function`
:foo
x = :foo
y = :function
z = :end

# Expected: highlight index `[end]` differently to block delimiter `end`
if foo[end]
end

# Expected: highlight as index `end`
foo[bar:end]

# Expected: don't highlight `:123`
x = :123

# Expected: don't highlight `:baz`
foo[bar:baz]

# Expected: highlight `:baz`
foo[:baz]

# Expected: highlight `:baz`
foo(:baz,:baz)

# Note that `: foo` is currently a valid quoted symbol, this will hopefully
# change in 0.4: https://github.com/JuliaLang/julia/issues/5997

# Expected: highlight `:foo`
[1 :foo]

# Expected: highlight `:end`
[1 :end]

# Expected: don't highlight `:foo`
for x=1:foo
    print(x)
end

## Number highlighting

# Expected: highlight all these as numbers
x = 123
x = 1.1
x = .5
x = -1
x = 0x123abcdef
x = 0o7
x = 0b1011
x = 2.5e-4
x = 2.5E-4
x = 1e+00
x = 2.5f-4
x = 0x.4p-1

# Expected: highlight these as numbers or built-ins
x = Inf
x = NaN

# Expected: highlight `123`
y = 123x

# Debatable: highlight `π` as a number or built-in
# (note that `πx` is a single symbol, not `π * x`)
x = π

# Note that `x.3` is currently equivalent to `.3x` or `.3 * x`, but this
# may change: https://github.com/JuliaLang/julia/issues/6770

max_stars_repo_name: jarrison/julia-syntax-test-cases

max_issues_repo_name: jarrison/julia-syntax-test-cases

max_forks_repo_name: jarrison/julia-syntax-test-cases


text:
@testset "PointClouds tests" begin

typealias V3d SVector{3,Float64}

# Unit cube of points in 3D
position = [0.0  1.0  1.0  1.0  1.0  0.0  0.0  0.0
            0.0  0.0  1.0  0.0  1.0  1.0  1.0  0.0
            0.0  0.0  0.0  1.0  1.0  0.0  1.0  1.0
           ]
cloud = PointCloud(position)
cloud[:intensity] = collect(1:8)

@testset "Column access" begin
    @test haskey(cloud, :intensity)
    @test sort(collect(keys(cloud))) == [:intensity]

    @test cloud[:intensity][1]   == 1
    @test cloud[:intensity][end] == 8

    @test positions(cloud) === cloud.positions
end

@testset "Row access" begin
    @test length(cloud) == 8

    subset = cloud[3:4]
    @test positions(subset)[1] == V3d(1,1,0)
    @test subset[:intensity][2] == 4

    dupcloud = [cloud; cloud]
    @test positions(dupcloud[1:8]) == positions(dupcloud[9:end])
    @test dupcloud[1:8][:intensity] == dupcloud[9:end][:intensity]

    splitres = split_cloud(cloud, cloud[:intensity] .% 2)
    @test length(splitres[1]) == 4
    @test splitres[1][:intensity] == collect(1:2:8)
    @test length(splitres[0]) == 4
    @test splitres[0][:intensity] == collect(2:2:8)
end

@testset "Spatial indexing" begin
    @test knn(cloud, [0,0,0], 1)[1] == [1]
    @test sort(knn(cloud, V3d(0,0,0), 4)[1]) == [1,2,6,8]
    @test inrange(cloud, [0,0,0], 0.1) == [1]
    @test sort(inrange(cloud, V3d(0,0,0), 1.1)) == [1,2,6,8]

    # multiple points
    multi = [V3d(0,0,0), V3d(1,1,1)]
    @test knn(cloud, multi, 1)[1] == [[1], [5]]
    @test inrange(cloud, multi, 1.1) == [[1,2,6,8], [3,4,5,7]]
end

@testset "PCA" begin
    centres = [
        V3d(0,0,0),
        V3d(10,0,0),
        V3d(0,10,0),
    ]
    desired_normals = [
        V3d(1,0,0),
        V3d(1,1,0)/sqrt(2),
        V3d(1,0,1)/sqrt(2),
    ]

    position = V3d[]
    for j=1:3
        Nj = desired_normals[j]
        for i=1:10
            r = V3d(randn(),randn(),randn())
            r = r - dot(Nj, r)*Nj
            push!(position, r + centres[j])
        end
    end

    cloud = PointCloud(position)

    add_normals!(cloud)

    N = cloud[:normals]
    N1 = N[knn(cloud, centres[1], 1)[1][1]]
    N2 = N[knn(cloud, centres[2], 1)[1][1]]
    N3 = N[knn(cloud, centres[3], 1)[1][1]]
    @test abs(abs(dot(V3d(N1), desired_normals[1])) - 1) < 10*eps(Float32)
    @test abs(abs(dot(V3d(N2), desired_normals[2])) - 1) < 10*eps(Float32)
    @test abs(abs(dot(V3d(N3), desired_normals[3])) - 1) < 10*eps(Float32)
end

end # testset


text:
# aggregate2: find minimum including NAs
@testset "Aggregate data" begin
	dfi = DataFrame(Temp=[10,11,14,1,2,missing,4],
	   		datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,6),
	         		  DateTime(2010,1,1,18),
	           		  DateTime(2010,1,2,0),DateTime(2010,1,2,6),
	           		  DateTime(2010,1,2,12),DateTime(2010,1,2,18),
					  ]);
	dfa = aggregate2(dfi,resol=Dates.Day(1),fce=minimum);
	@test dfa[1,1] == 10
	@test ismissing(dfa[2,1])
	# aggregate2: compute sum removing NAs
	dfa = aggregate2(dfi,resol=Dates.Day(1),fce=x->sum(collect(skipmissing(x))))
	@test dfa[1,1] == 10+11+14
	@test dfa[2,1] == 1+2+4

	dfc = aggregate2(dfi,resol=Dates.Day(1),fce=maximum,rename=true);
	@test names(dfc) == [:Temp,:datetime]
end

@testset "Time regularization" begin
	dfi = DataFrame(Temp=[10,11,14,1,2,missing,4],
	   		datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,6),
	         		  DateTime(2010,1,1,18),
	           		  DateTime(2010,1,2,0),DateTime(2010,1,2,6),
	           		  DateTime(2010,1,2,12),DateTime(2010,1,2,18),
					  ]);
	reg_sample = time2regular(dfi,timecol=:datetime,resol=Dates.Hour(6))
	@test size(reg_sample,1) == 8
	@test ismissing(reg_sample[3,1])

	# Test irregular
	dfi = DataFrame(Temp=[10,11,14,1,2,missing,4],
	   		datetime=[DateTime(2010,1,1,0),DateTime(2010,1,1,6),
	         		  DateTime(2010,1,1,18),
	           		  DateTime(2010,1,2,0),DateTime(2010,1,2,6),
	           		  DateTime(2010,1,2,12),DateTime(2010,1,2,18),
					  ]);
	@test isregular(dfi[:datetime]) == false

	# Test regular
	timevec = collect(DateTime(2000,1,1):Dates.Hour(1):DateTime(2001,1,2,3,0,0));
	@test isregular(timevec) == true
end

@testset "Cut 2 time interval" begin
	datain = DataFrame(datetime=collect(DateTime(2000,1,1):Dates.Day(1):DateTime(2000,1,12)),
						 somevalues=collect(1:1:12))
	datatest = deepcopy(datain);
	# do not cut anything as the interval is outside given time vector
	starttime,endtime = DateTime(2000,2,2),DateTime(2000,2,11);
	cut2interval!(datatest,starttime,endtime,keepedges=(false,true))
	@test datatest == datain
	# cut, but keep starting point
	starttime,endtime = DateTime(2000,1,1),DateTime(2000,1,11);
	datatest = deepcopy(datain);
	cut2interval!(datatest,starttime,endtime,keepedges=(true,false))
	@test datatest[:datetime] == datain[:datetime][1:10]
	# cut, removing starting point and keeping end
	starttime,endtime = DateTime(2000,1,1),DateTime(2000,1,9);
	datatest = deepcopy(datain);
	cut2interval!(datatest,starttime,endtime,keepedges=(false,true))
	@test datatest[:datetime] == datain[:datetime][2:9]
	# Do not cut anything because of the specific setting (end time outside time
	# interval & keep first point)
	starttime,endtime = DateTime(2000,1,1),DateTime(2001,12,9);
	datatest = deepcopy(datain);
	cut2interval!(datatest,starttime,endtime,keepedges=(true,true))
	@test datatest == datain
end


text:
module IndexedTables

using PooledArrays, SparseArrays, Statistics, WeakRefStrings

using OnlineStatsBase: OnlineStat, fit!

using StructArrays: StructVector, StructArray, fieldarrays,
    refine_perm!, collect_structarray,
    append!!, replace_storage, GroupPerm,
    roweq, index_type

import Tables, TableTraits, IteratorInterfaceExtensions, TableTraitsUtils

import DataValues: DataValue, DataValueArray, isna

import DataAPI: Between, All

using InvertedIndices: Not
import StatsBase: transform

import Base:
    show, eltype, length, getindex, setindex!, ndims, map, convert, keys, values,
    ==, broadcast, empty!, copy, similar, sum, merge, merge!, mapslices,
    permutedims, sort, sort!, iterate, pairs, reduce, push!, size, permute!, issorted,
    sortperm, summary, resize!, vcat, append!, copyto!, view, tail,
    tuple_type_cons, tuple_type_head, tuple_type_tail, in, convert

#-----------------------------------------------------------------------# exports
export
    # macros
    @cols,
    # types
    AbstractNDSparse, All, ApplyColwise, Between, ColDict, Columns, IndexedTable,
    Keys, NDSparse, NextTable, Not,
    # functions
    aggregate!, antijoin, asofjoin, collect_columns, colnames,
    column, columns, convertdim, dimlabels, flatten, flush!, groupby, groupjoin,
    groupreduce, innerjoin, insertafter!, insertbefore!, insertcols, insertcolsafter,
    insertcolsbefore, leftgroupjoin, leftjoin, map_rows, naturalgroupjoin, naturaljoin,
    ncols, ndsparse, outergroupjoin, outerjoin, pkeynames, pkeys,
    reducedim_vec, reindex, rename, rows, select, selectkeys, selectvalues,
    stack, summarize, table, transform, unstack, update!, where, dropmissing, dropna

const Tup = Union{Tuple,NamedTuple}
const DimName = Union{Int,Symbol}

include("utils.jl")
include("columns.jl")
include("indexedtable.jl")
include("ndsparse.jl")
include("collect.jl")

#=
# Poor man's traits

# These support `colnames` and `columns`
const TableTrait = Union{AbstractVector, IndexedTable, NDSparse}

# These support `colnames`, `columns`,
# `pkeynames`, `permcache`, `cacheperm!`
=#

const Dataset = Union{IndexedTable, NDSparse}

# no-copy convert
_convert(::Type{IndexedTable}, x::IndexedTable) = x
function _convert(::Type{NDSparse}, t::IndexedTable)
    NDSparse(rows(t, pkeynames(t)), rows(t, excludecols(t, pkeynames(t))),
             copy=false, presorted=true)
end

function _convert(::Type{IndexedTable}, x::NDSparse)
    convert(IndexedTable, x.index, x.data;
            perms=x._table.perms,
            presorted=true, copy=false)
end

ndsparse(t::IndexedTable; kwargs...) = _convert(NDSparse, t; kwargs...)
table(t::NDSparse; kwargs...) = _convert(IndexedTable, t; kwargs...)

include("sortperm.jl")
include("indexing.jl") # x[y]
include("selection.jl")
include("reduce.jl")
include("flatten.jl")
include("join.jl")
include("reshape.jl")
include("tables.jl")
include("tabletraits.jl")

end # module


text:
module TurbulenceClosures

export
  IsotropicViscosity,
  ConstantIsotropicDiffusivity,
  ConstantAnisotropicDiffusivity,
  AnisotropicBiharmonicDiffusivity,
  TwoDimensionalLeith,
  ConstantSmagorinsky,
  SmagorinskyLilly,
  BlasiusSmagorinsky,
  AnisotropicMinimumDissipation,
  RozemaAnisotropicMinimumDissipation,
  VerstappenAnisotropicMinimumDissipation,

  TurbulentDiffusivities,
  calculate_diffusivities!,

  ∇_κ_∇c,
  ∇_κ_∇T,
  ∇_κ_∇S,
  ∂ⱼ_2ν_Σ₁ⱼ,
  ∂ⱼ_2ν_Σ₂ⱼ,
  ∂ⱼ_2ν_Σ₃ⱼ,

  cell_diffusion_timescale

using Oceananigans: @hascuda

@hascuda using CUDAdrv, CUDAnative

using
  Oceananigans,
  Oceananigans.Grids,
  Oceananigans.Operators,
  GPUifyLoops

import Oceananigans: with_tracers

using Oceananigans: AbstractArchitecture, AbstractGrid,
                    TracerFields, device, launch_config,
                    buoyancy_perturbation, buoyancy_frequency_squared, ∂x_b, ∂y_b, ∂z_b

####
#### Molecular viscosity and thermal diffusivity definitions
#### Approximate viscosities and thermal diffusivities for seawater at 20ᵒC and 35 psu,
#### according to Sharqawy et al., "Thermophysical properties of seawater: A review of
#### existing correlations and data" (2010).
####

const ν₀ = 1.05e-6
const κ₀ = 1.46e-7

####
#### Abstract types
####

"""
    TurbulenceClosure{FT}

Abstract supertype for turbulence closures with model parameters stored as properties of
type `FT`.
"""
abstract type TurbulenceClosure{FT} end

"""
    IsotropicViscosity{FT} <: TurbulenceClosure{FT}

Abstract supertype for turbulence closures that are defined by an isotropic viscosity
and isotropic diffusivities with model parameters stored as properties of type `FT`.
"""
abstract type IsotropicViscosity{FT} <: TurbulenceClosure{FT} end

"""
    TensorDiffusivity{FT} <: TurbulenceClosure{FT}

Abstract supertype for turbulence closures that are defined by a tensor viscosity and
tensor diffusivities with model parameters stored as properties of type `FT`.
"""
abstract type TensorDiffusivity{FT} <: TurbulenceClosure{FT} end

"""
    AbstractSmagorinsky{FT}

Abstract supertype for large eddy simulation models based off the model described
by Smagorinsky with model parameters stored as properties of type `FT`.
"""
abstract type AbstractSmagorinsky{FT} <: IsotropicViscosity{FT} end

"""
    AbstractAnisotropicMinimumDissipation{FT}

Abstract supertype for large eddy simulation models based on the anisotropic minimum
dissipation principle with model parameters stored as properties of type `FT`.
"""
abstract type AbstractAnisotropicMinimumDissipation{FT} <: IsotropicViscosity{FT} end

"""
    AbstractLeith{FT}

Abstract supertype for large eddy simulation models based on the Leith viscosity
principle with model parameters stored as properties of type `FT`.
"""
abstract type AbstractLeith{FT} <: IsotropicViscosity{FT} end

####
#### 'Tupled closure' implementation
####

for stress_div in (:∂ⱼ_2ν_Σ₁ⱼ, :∂ⱼ_2ν_Σ₂ⱼ, :∂ⱼ_2ν_Σ₃ⱼ)
    @eval begin
        @inline function $stress_div(i, j, k, grid::AbstractGrid{FT}, closure_tuple::Tuple, U,
                                     K_tuple, args...) where FT

            stress_div_ijk = zero(FT)

            ntuple(Val(length(closure_tuple))) do α
                @inbounds closure = closure_tuple[α]
                @inbounds K = K_tuple[α]
                stress_div_ijk += $stress_div(i, j, k, grid, closure, U, K, args...)
            end

            return stress_div_ijk
        end
    end
end

@inline function ∇_κ_∇c(i, j, k, grid::AbstractGrid{FT}, closure_tuple::Tuple,
                        c, tracer_index, K_tuple, args...) where FT
    flux_div_ijk = zero(FT)

    ntuple(Val(length(closure_tuple))) do α
        @inbounds closure = closure_tuple[α]
        @inbounds K = K_tuple[α]
        flux_div_ijk +=  ∇_κ_∇c(i, j, k, grid, closure, c, tracer_index, K, args...)
    end

    return flux_div_ijk
end

function calculate_diffusivities!(K_tuple::Tuple, arch, grid, closure_tuple::Tuple, args...)
    ntuple(Val(length(closure_tuple))) do α
        @inbounds closure = closure_tuple[α]
        @inbounds K = K_tuple[α]
        calculate_diffusivities!(K, arch, grid, closure, args...)
    end

    return nothing
end

TurbulentDiffusivities(arch::AbstractArchitecture, grid::AbstractGrid, tracers, closure_tuple::Tuple) =
    Tuple(TurbulentDiffusivities(arch, grid, tracers, closure) for closure in closure_tuple)

with_tracers(tracers, closure_tuple::Tuple) =
    Tuple(with_tracers(tracers, closure) for closure in closure_tuple)

####
#### Include module code
####

# Fallback constructor for diffusivity types without precomputed diffusivities:
TurbulentDiffusivities(arch::AbstractArchitecture, grid::AbstractGrid, args...) = nothing

include("turbulence_closure_utils.jl")
include("closure_operators.jl")
include("velocity_tracer_gradients.jl")

include("closure_tuples.jl")

include("turbulence_closure_implementations/constant_isotropic_diffusivity.jl")
include("turbulence_closure_implementations/constant_anisotropic_diffusivity.jl")
include("turbulence_closure_implementations/anisotropic_biharmonic_diffusivity.jl")
include("turbulence_closure_implementations/leith_enstrophy_diffusivity.jl")
include("turbulence_closure_implementations/smagorinsky_lilly.jl")
include("turbulence_closure_implementations/blasius_smagorinsky.jl")
include("turbulence_closure_implementations/verstappen_anisotropic_minimum_dissipation.jl")
include("turbulence_closure_implementations/rozema_anisotropic_minimum_dissipation.jl")

include("turbulence_closure_diagnostics.jl")

####
#### Some value judgements here
####

"""
    AnisotropicMinimumDissipation

An alias for `VerstappenAnisotropicMinimumDissipation`.
"""
const AnisotropicMinimumDissipation = VerstappenAnisotropicMinimumDissipation

"""
    ConstantSmagorinsky

An alias for `SmagorinskyLilly`.
"""
const ConstantSmagorinsky = SmagorinskyLilly

end # module


text:
# const urlCrest = "https://public-crest.eveonline.com/"
const urlCrest = "https://crest-tq.eveonline.com/"

const urlReplace = Dict{Char, Char}('/'=>'-', '?'=>'.')
url_to_filename(url::AbstractString) = replace(replace(url, urlCrest[1:end-1], ""), collect(keys(urlReplace)), c->urlReplace[c[1]]) * ".json"

int_keys(assoc::Associative) = Dict(parse(Int, k)=>v for (k,v) in assoc)

function json_read(fileName::AbstractString)
	data = nothing
	try
		open(fileName, "r") do s
			data = JSON.parse(s)
		end
	catch e
		if !isa(e, SystemError)
			rethrow()
		end
	end
	if isa(data, Associative)
		intKeys = all(keys(data)) do k
			try
				parse(Int, k)
				return true
			catch
				return false
			end
		end
		if intKeys
			data = int_keys(data)
		end
	end
	return data
end

function json_write(fileName::AbstractString, data; indent::Int = 2)
	open(fileName, "w+") do s
		JSON.print(s, data, indent)
	end
end

function read_from_cache(url::AbstractString, timeout::Float64)
	fileName = joinpath(cacheDir, url_to_filename(url))
	data = json_read(fileName)
	if time() - mtime(fileName) >= timeout
		data = nothing
	end
	return data
end

function write_to_cache(url::AbstractString, data)
	if !isdir(cacheDir)
		mkdir(cacheDir)
	end
	fileName = joinpath(cacheDir, url_to_filename(url))
	json_write(fileName, data)
	nothing
end

function clear_cache(ext::AbstractString)
	if isdir(cacheDir)
		map(f->endswith(f, ext) && rm(joinpath(cacheDir, f)), readdir(cacheDir))
	end
end

function get_crest(url::AbstractString, auth, timeoutHours::Float64 = convert(Float64, Inf))
	res = read_from_cache(url, timeoutHours * 60 * 60)
	if res != nothing
		return res
	end
	local queryRes
	orgUrl = url
	headers = Dict{AbstractString, AbstractString}()
	if auth != nothing
		headers["Authorization"] = auth["token_type"]*" "*auth["access_token"]
	end
	while true
		info("GET from $url")
		resp = Requests.get(url; headers=headers)
		if resp.status != 200
			info("Response status: $(resp.status)")
			return nothing
		end
		queryRes = JSON.parse(String(resp.data))
		if res == nothing
			res = queryRes
		else
			append!(res["items"], queryRes["items"])
		end
		if !haskey(queryRes, "next")
			break
		end
		url = queryRes["next"]["href"]
	end
	@assert !haskey(queryRes, "totalcount") || queryRes["totalcount"] == length(res["items"])
	if haskey(res, "items")
		res = res["items"]
	end
	write_to_cache(orgUrl, res)
	return res
end


text:
module MPOPIS

using Reexport
using CSV
using Printf
using LinearAlgebra
using IntervalSets
using Distributions
using StatsBase
using Random
@reexport import Random.seed!
using CovarianceEstimation
import CovarianceEstimation.LinearShrinkageTarget
import CovarianceEstimation.SimpleCovariance
@reexport using ReinforcementLearning
import ReinforcementLearning.AbstractEnv
import ReinforcementLearning.RLBase
using Plots
@reexport import Plots.plot
using ProgressMeter
using Dates

export
    MPPI_Policy,
    GMPPI_Policy,
    CEMPPI_Policy,
    CMAMPPI_Policy,
    μAISMPPI_Policy,
    μΣAISMPPI_Policy,
    PMCMPPI_Policy,
    NESMPPI_Policy,
    Track,
    CarRacingEnv,
    MultiCarRacingEnv,
    DroneEnv,
    within_track,
    calculate_β,
    exceed_β,
    block_diagm,
    _update_states_envs2env,
    _update_states_env2envs,
    simulate_car_racing,
    simulate_mountaincar,
    simulate_cartpole

abstract type AbstractWeightMethod end
abstract type AbstractPathIntegralPolicy end
Random.seed!(pol::AbstractPathIntegralPolicy, seed) = Random.seed!(pol.rng, seed)

abstract type AbstractGMPPI_Policy <: AbstractPathIntegralPolicy end

struct Cross_Entropy <: AbstractWeightMethod
    elite_threshold::Float64
    num_elite_samples::Int
end

struct Information_Theoretic <: AbstractWeightMethod
    λ::Float64
end

include("envs/car_racing_tracks/car_racing_tracks.jl")
include("envs/car_racing.jl")
include("envs/multi-car_racing.jl")
include("envs/drone_env.jl")
include("utils.jl")
include("mppi_mpopi_policies.jl")
include("envs/plots.jl")

include("examples/example_utils.jl")
include("examples/car_example.jl")
include("examples/mountaincar_example.jl")
include("examples/cartpole_example.jl")

end # module

max_stars_repo_name: sisl/MPOPIS

max_issues_repo_name: sisl/MPOPIS

max_forks_repo_name: sisl/MPOPIS


text:
"""
    JopMix(s::JetAbstractSpace{T,N}, nmix::NTuple{N,Int})

2D spatial mix, and is used to add smoothness to a domain.  JopMix is ported
from the CVX frequency domain FWI tools in SeisSpace.

# Notes:
* This appears to over-lap in funtionality with JopRoughness.

# Examples:
## 2D
```julia
nz, nx = 11,11
spc = JetSpace(Float32, nz, nx)
A = JopMix(spc, (5,5))
y = A * rand(dom)
x = A' * rand(rng)
```
"""
function JopMix(spc::JetAbstractSpace{T,N}, nmix::NTuple{N,Int}) where {T,N}
    for i in nmix
        i < 0 && error("mix lengths must be >= 0 -- $(nmix)")
    end
    JopLn(dom = spc, rng = spc, df! = JopMix_df!, df′! = JopMix_df′!, s= (nmix=nmix,))
end

export JopMix

function JopMix_df!(rngvec::AbstractArray{T,2}, domvec::AbstractArray{T,2}; nmix, kwargs...) where {T}
    rngvec .= 0
    nz,nx = size(domvec)
    lenz = div(nmix[1], 2)
    lenx = div(nmix[2], 2)
    wgt = (2 * lenz + 1) * (2 * lenx + 1)

    for kz = 1:nz
        jz1 = clamp(kz - lenz, 1, nz)
        jz2 = clamp(kz + lenz, 1, nz)
        for kx = 1:nx
            jx1 = clamp(kx - lenx, 1, nx)
            jx2 = clamp(kx + lenx, 1, nx)
            for jx = jx1:jx2
                for jz = jz1:jz2
                    rngvec[kz,kx] += domvec[jz,jx]
                end
            end
        end
    end
    rngvec ./= wgt
    rngvec
end

function JopMix_df′!(domvec::AbstractArray{T,2}, rngvec::AbstractArray{T,2}; nmix, kwargs...) where {T}
    domvec .= 0
    nz,nx = size(domvec)
    lenz = div(nmix[1], 2)
    lenx = div(nmix[2], 2)
    wgt = (2 * lenz + 1) * (2 * lenx + 1)

    for kz = 1:nz
        jz1 = clamp(kz - lenz, 1, nz)
        jz2 = clamp(kz + lenz, 1, nz)
        for kx = 1:nx
            jx1 = clamp(kx - lenx, 1, nx)
            jx2 = clamp(kx + lenx, 1, nx)
            for jx = jx1:jx2
                for jz = jz1:jz2
                    domvec[kz,kx] += rngvec[jz,jx]
                end
            end
        end
    end
    domvec ./= wgt
    domvec
end


text:
"""
    second-order cone
"""
function second_order_cone_projection(z)
    n = length(z)

    z0 = z[1]
    z1 = view(z, 2:n)

    if norm(z1) <= z0
        return z, true
    elseif norm(z1) <= -z0
        return zero(z), false
    else
        a = 0.5 * (1.0 + z0 / norm(z1))
        z_proj = zero(z)
        z_proj[1:end-1] = a * z1
        z_proj[end] = a * norm(z1)
        return z_proj, false
    end
end

function second_order_cone_product(z, s)
    n = length(z)
    SVector{n}([z' * s; z[1] * view(s, 2:n) + s[1] * view(z, 2:n)])
end

# check that second-order cone constraints are satisfied
function second_order_cone_check(x, idx_soc)
    for idx in idx_soc
        !second_order_cone_projection(view(x, idx))[2] && return false
    end
    return true
end

# check that inequality (non-negative orthant) constraints are satisfied
function inequality_check(x, idx_ineq)
    for i in idx_ineq
        if x[i] <= 0.0
            return false
        end
    end
    return true
end

function cone_check(x, idx_ineq, idx_soc)
    !inequality_check(x, idx_ineq) && return true
    !second_order_cone_check(x, idx_soc) && return true
end

function cone_plot(λ::AbstractVector{T}, Δ::AbstractVector{T};
		ϵ = 1e-20, N::Int = 100, verbose::Bool = false, show_α::Bool = false) where {T}
	plt = plot(aspect_ratio = 1.0, fmt = :svg)

	# unit circle
	θ = Vector(0:2π/N:(1.0+1/N)*2π)
	circ = [cos.(θ) , sin.(θ)]
	plot!(plt, circ[1], circ[2], label = false)

	# Keypoints
	λp = project(λ)
	λΔp = project(λ + Δ)
	verbose && @show soc_value(λ )
	verbose && @show soc_value(λ + Δ)
	scatter!(plt, λΔp[1:1],  λΔp[2:2],  markersize = 4.0, color = :orange, label = "λ+Δ")
	scatter!(plt, λp[1:1],   λp[2:2],   markersize = 8.0, color = :red,    label = "λ")

	if show_α
		α = soc_step_length(λ, Δ, τ = 1.0, ϵ = ϵ)
		λαΔp = project(λ + α * Δ)
		verbose && @show α
		verbose && @show soc_value(λ + α * Δ)
		scatter!(plt, λαΔp[1:1], λαΔp[2:2], markersize = 2.0, color = :black,   label = "λ+αΔ")
	end

	# Line
	A = Vector(0:1/N:1.0)
	L = [[project(λ + a * Δ)[1] for a in A], [project(λ + a * Δ)[2] for a in A]]
	plot!(plt, L[1], L[2], label = false)

	display(plt)
	return nothing
end

function project(x::AbstractVector{T}; ϵ::T = 1e-15) where {T}
	if length(x) == 3
		project_3D(x, ϵ = ϵ)
	elseif length(x) == 2
		project_2D(x, ϵ = ϵ)
	end
end

function project_3D(x::AbstractVector{T}; ϵ::T = 1e-15) where {T}
	@assert length(x) == 3
	xs = x[1]
	xv = x[2:end]
	r = norm(xv) / (xs + ϵ)
	v = xv ./ norm(xv)
	p = r * v
	return p
end

function project_2D(x::AbstractVector{T}; ϵ::T = 1e-15) where {T}
	@assert length(x) == 2
	xs = x[1]
	xv = x[2:end]
	r = norm(xv) / (xs + ϵ)
	v = xv ./ norm(xv)
	p = r * v
	return [p[1], 0.0]
end


text:
export read_cq_bvecs, read_cq_fvecs, CQ_parameters, dump_CQ_parameters

# Read float binary file produced by CQ, like the codebooks D
function read_cq_fvecs(fname::AbstractString)

  dim, count = zero(Cint), zero(Cint)
  vectors = Matrix{Cfloat}

  open(fname, "r") do fid
    count = read(fid, Cint, 1)[1]
    dim   = read(fid, Cint, 1)[1]
    vectors = read(fid, Cfloat, dim, count)
  end

  return vectors
end

# Read int binary file produced by CQ, like the codes B
function read_cq_bvecs(fname::AbstractString)

  dim, count = zero(Cint), zero(Cint)
  codes = Matrix{Cint}

  open(fname, "r") do fid
    count = read(fid, Cint, 1)[1]
    dim   = read(fid, Cint, 1)[1]
    codes = read(fid, Cint, dim, count)
  end

  return codes
end

# Parameters that CQ receives
@with_kw mutable struct CQ_parameters
  PQ::Bool=false
  NCQ::Bool=false
  CQ::Bool=true
  Search::Bool=true

  # global parameters
  # points_count=1000000
  points_count::Int=100000
  dictionaries_count::Int=8
  words_count::Int=256
  space_dimension::Int=128
  # points_file=/home/julieta/Desktop/local-search-quantization/data/sift/sift_base.fvecs
  points_file::String="/home/julieta/Desktop/local-search-quantization/data/sift/sift_learn.fvecs"
  output_file_prefix::String="/home/julieta/Desktop/CQ/build/temp/"
  max_iter::Int=30

  # PQ parameters
  distortion_tol::Float32=0.0001
  read_partition::Int=0
  partition_file::String=""
  # if 101 then using closure cluster, else lloyd kmeans
  kmeans_method::Int=101

  # NCQ and CQ parameters
  num_sep::Int=20
  # initial from outside, if 1 then set the file name of dictinary and codes
  initial_from_outside=0
  dictionary_file::String=""
  binary_codes_file::String=""

  # CQ parameters
  mu::Float32=0.0004f0

  # Search parameters
  queries_count::Int=10000
  groundtruth_length::Int=100
  result_length::Int=1000
  queries_file::String="/home/julieta/Desktop/local-search-quantization/data/sift/sift_query.fvecs"
  groundtruth_file::String="/home/julieta/Desktop/local-search-quantization/data/sift/sift_groundtruth.ivecs"
  trained_dictionary_file::String="/home/julieta/Desktop/CQ/build/temp/D"
  trained_binary_codes_file::String="/home/julieta/Desktop/CQ/build/temp/B"
  output_retrieved_results_file::String="/home/julieta/Desktop/CQ/build/temp/recall"
end

# Print CQ parameters to a config file that the CQ binary can consume
function dump_CQ_parameters(p::CQ_parameters, fname::String)
  open(fname, "w") do fid
    for (name, typ) in zip(fieldnames(CQ_parameters), CQ_parameters.types)
      if typ == Bool || typ == Int
        println(fid, name, "=", Int(getfield(p, name)))
      else
        println(fid, name, "=", getfield(p, name))
      end
    end
  end
end


text:
# PyTFNeuralPose2

ENV["PYTHON"] = "/usr/bin/python3.6"

using Pkg
Pkg.build("PyCall")

using PyCall

tf = pyimport("tensorflow")

tf.keras.models.load_model

Base.cd(ENV["USER"]*"/learning-odometry/")

#


text:
# Copyright (c) 2018: Matthew Wilhelm & Matthew Stuber.
# This code is licensed under MIT license (see LICENSE.md for full details)
#############################################################################
# EAGO
# A development environment for robust and global optimization
# See https://github.com/PSORLab/EAGO.jl
#############################################################################
# src/eago_optimizer/logging/log.jl
# Defines all type used to store solution information at given iterations
#############################################################################

"""
$(TYPEDEF)

A structure used to store information on the history of the solution procedure
for generating convergence plots and other analysis.

$(TYPEDFIELDS)
"""
Base.@kwdef mutable struct Log
    "Storage for lower bound calculated for current node."
    current_lower_bound::Vector{Float64} = Float64[]
    "Storage for upper bound calculated for current node."
    current_upper_bound::Vector{Float64} = Float64[]
    "Storage for preprocessing time of each iteration."
    preprocessing_time::Vector{Float64} = Float64[]
    "Storage for lower bounding time of each iteration."
    lower_problem_time::Vector{Float64} = Float64[]
    "Storage for upper bounding time of each iteration."
    upper_problem_time::Vector{Float64} = Float64[]
    "Storage for postprocessing time of each iteration."
    postprocessing_time::Vector{Float64} = Float64[]
    "Storage for preprocessing feasibility of each iteration."
    preprocessing_feas::Vector{Bool} = Bool[]
    "Storage for lower bounding feasibility of each iteration."
    lower_problem_feas::Vector{Bool} = Bool[]
    "Storage for upper bounding feasibility of each iteration."
    upper_problem_feas::Vector{Bool} = Bool[]
    "Storage for postprocessing feasibility of each iteration."
    postprocessing_feas::Vector{Bool} = Bool[]
    "Storage for best (global) lower bound at each iteration."
    global_lower_bound::Vector{Float64} = Float64[]
    "Storage for best (global) upper bound at each iteration."
    global_upper_bound::Vector{Float64} = Float64[]
    "Number of nodes at each iteration."
    node_count::Vector{Int64} = Int64[]
    "Run time at each iteration."
    run_time::Vector{Float64} = Float64[]
end

function Base.isempty(x::Log)
    for field in fieldnames(Log)
        if !isempty(getfield(x, field))
            return false
        end
    end
    return true
end

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
#
# This file is part of the SlowActors.jl Julia package,
# MIT license, part of https://github.com/JuliaActors
#

function _start_on_send(lk::Link{Mailbox})
    if !isnothing(lk.chn.A)
        if isnothing(lk.chn.A.t) || istaskdone(lk.chn.A.t)
            startActor(lk.chn)
        end
    end
end

"""
```
send!(lk::Link{Mailbox}, msg::Msg)
send!(lk::Link{Mailbox}, msg...)
```

Send a message `msg` to a slow actor. Start it if it doesn't run.

**Note:** You must reeimplement those as methods of `Actors.send!`.
"""
Actors.send!(lk::Link{Mailbox}, msg::Msg) = (push!(lk.chn, msg); _start_on_send(lk))
Actors.send!(lk::Link{Mailbox}, msg...) = (push!(lk.chn, msg); _start_on_send(lk))

# ------------------------------------------
# In order to make Actors.jl receive! work,
# you must reimplement the following three
# functions with your Link parameter:
# ------------------------------------------
Base.isready(mbx::Mailbox) = !isempty(mbx)
function Base.fetch(mbx::Mailbox) # this waits for a message
    timedwait(60, pollint=0.01) do
        !isempty(mbx)
    end == :ok ?
        lock(()->first(mbx.box), mbx.lock) :
        :timed_out
end
Base.take!(mbx::Mailbox) = popfirst!(mbx)


text:
function solve1(numbers, counts)
    gamma = 0
    for i in counts
        gamma = gamma*2 + (i > length(numbers)/2)
    end
    epsilon = (2^length(counts) - 1) ⊻ gamma
    return gamma * epsilon
end

function solve2(nums, pow)
    function descend_tree(numbers, power, inverse=false)
        if length(numbers) == 1
            return numbers[1]
        end
        matching = 2^power .& numbers
        sort!(matching)
        sort!(numbers, by=(x) -> x & (2^power))

        index = length(numbers) ÷ 2 + 1
        mid = matching[index]
        if (mid == 0) != inverse
            func = x -> (x & 2^power) == 0
        else
            func = x -> !((x & 2^power) == 0)
        end
        filter!(func, numbers)
        return descend_tree(numbers, power-1, inverse)
    end
    return descend_tree(copy(nums), pow, false) * descend_tree(copy(nums), pow, true)
end

function setup()
    counts, numbers = [], []
    for line in readlines()
        for i in eachindex(line)
            if length(counts) < i
                push!(counts, 0)
            end
            counts[i] += '1' == line[i]
        end
        push!(numbers, parse(Int, line, base=2))
    end
    println(solve1(numbers, counts))
    println(solve2(numbers, length(counts) - 1))
end

setup()

max_stars_repo_name: LiquidFun/adventofcode

max_issues_repo_name: LiquidFun/adventofcode

max_forks_repo_name: LiquidFun/adventofcode


text:
struct Solution{T} <: AbstractSolution
	iset::Vector{Hyperrectangle{T,Vector{T},Vector{T}}}
	options::Options
end

# constructor with no options
Solution(iset::Vector{Hyperrectangle}) = Solution(iset,Options())

function save_data(iter::Int, sol::Solution)
	save("solution.jld2", "$iter", sol)
end

function load_data(name::String)
	data = load(name)
	key = collect(Base.keys(data))[1]
	iter = parse(Int,key)
	sol = data[key]
	return iter, sol.iset
end


text:
import Oceananigans.Grids: required_halo_size

"""
    AnisotropicBiharmonicDiffusivity{FT, KH, KZ}

Parameters for anisotropic biharmonic diffusivity models.
"""
struct AnisotropicBiharmonicDiffusivity{FT, KX, KY, KZ} <: AbstractTensorDiffusivity
    νx :: FT
    νy :: FT
    νz :: FT
    κx :: KX
    κy :: KY
    κz :: KZ
end

"""
    AnisotropicBiharmonicDiffusivity(FT=Float64; νx=0, νy=0, νz=0, νh=nothing, κx=0, κy=0, κz=0, κh=nothing)

Returns parameters for a fourth-order, anisotropic biharmonic diffusivity closure with
constant x-, y, and z-direction biharmonic viscosities `νx`, `νy`, and `νz`,
and constant x-, y, and z-direction biharmonic diffusivities `κx`, `κy`, and `κz`,
`κx`, `κy`, and `κz` may be `NamedTuple`s with fields corresponding
to each tracer, or a single number to be a applied to all tracers.

If `νh` or `κh` are provided, then `νx = νy = νh` or `κx = κy = κh`.

The tracer flux divergence associated with an anisotropic biharmonic diffusivity is, for example

```math
    ∂ᵢ κᵢⱼ ∂ⱼc = - [κx ∂⁴x + κy ∂⁴y + κz ∂⁴z] c
```

"""
function AnisotropicBiharmonicDiffusivity(FT=Float64; νx=0, νy=0, νz=0, κx=0, κy=0, κz=0, νh=nothing, κh=nothing)
    if !isnothing(νh)
        νx = νh
        νy = νh
    end

    if !isnothing(κh)
        κx = κh
        κy = κh
    end

    κx = convert_diffusivity(FT, κx)
    κy = convert_diffusivity(FT, κy)
    κz = convert_diffusivity(FT, κz)

    return AnisotropicBiharmonicDiffusivity(FT(νx), FT(νy), FT(νz), κx, κy, κz)
end

required_halo_size(closure::AnisotropicBiharmonicDiffusivity) = 2

function with_tracers(tracers, closure::AnisotropicBiharmonicDiffusivity)
    κx = tracer_diffusivities(tracers, closure.κx)
    κy = tracer_diffusivities(tracers, closure.κy)
    κz = tracer_diffusivities(tracers, closure.κz)
    return AnisotropicBiharmonicDiffusivity(closure.νx, closure.νy, closure.νz, κx, κy, κz)
end

calculate_diffusivities!(K, arch, grid, closure::AnisotropicBiharmonicDiffusivity, args...) = nothing

@inline ∂ⱼ_2ν_Σ₁ⱼ(i, j, k, grid, clock, closure::AnisotropicBiharmonicDiffusivity, U, args...) = (
    - closure.νx * ∂⁴xᶠᵃᵃ(i, j, k, grid, U.u)
    - closure.νy * ∂⁴yᵃᶜᵃ(i, j, k, grid, U.u)
    - closure.νz * ∂⁴zᵃᵃᶜ(i, j, k, grid, U.u)
    )

@inline ∂ⱼ_2ν_Σ₂ⱼ(i, j, k, grid, clock, closure::AnisotropicBiharmonicDiffusivity, U, args...) = (
    - closure.νx * ∂⁴xᶜᵃᵃ(i, j, k, grid, U.v)
    - closure.νy * ∂⁴yᵃᶠᵃ(i, j, k, grid, U.v)
    - closure.νz * ∂⁴zᵃᵃᶜ(i, j, k, grid, U.v)
    )

@inline ∂ⱼ_2ν_Σ₃ⱼ(i, j, k, grid, clock, closure::AnisotropicBiharmonicDiffusivity, U, args...) = (
    - closure.νx * ∂⁴xᶜᵃᵃ(i, j, k, grid, U.w)
    - closure.νy * ∂⁴yᵃᶜᵃ(i, j, k, grid, U.w)
    - closure.νz * ∂⁴zᵃᵃᶠ(i, j, k, grid, U.w)
    )

@inline function ∇_κ_∇c(i, j, k, grid, clock, closure::AnisotropicBiharmonicDiffusivity,
                        c, ::Val{tracer_index}, args...) where tracer_index

    @inbounds κx = closure.κx[tracer_index]
    @inbounds κy = closure.κy[tracer_index]
    @inbounds κz = closure.κz[tracer_index]

    return (- κx * ∂⁴xᶜᵃᵃ(i, j, k, grid, c)
            - κy * ∂⁴yᵃᶜᵃ(i, j, k, grid, c)
            - κz * ∂⁴zᵃᵃᶜ(i, j, k, grid, c)
           )
end


text:
#This is for 2D domainwall Fermion
abstract type AbstractFermionfields_3D{NC} <: AbstractFermionfields{NC,3}
end

function Base.setindex!(x::AbstractFermionfields_3D{NC},v,i1,i2, i5,i6,i7)  where {NC}
    @inbounds x.w[i7][i1,i2, i5,i6] = v
end

function Base.getindex(x::AbstractFermionfields_3D{NC},i1,i2, i5,i6,i7) where {NC}
    @inbounds return x.w[i7][i1,i2, i5,i6]
end

Base.length(x::T) where T <: AbstractFermionfields_3D = x.L5*length(x.w[1])

function Base.setindex!(x::T,v,i)  where T <: AbstractFermionfields_3D
    #i = (iL-1)*NN+index
    NN = length(x.w[1])
    i4D = (i-1) % NN + 1
    iL = (i-i4D) ÷ NN + 1
    @inbounds x.w[iL][i4D] = v
end

function Base.getindex(x::T,i) where T <: AbstractFermionfields_3D
    #i = (iL-1)*NN+index
    NN = length(x.w[1])
    i4D = (i-1) % NN + 1
    iL = (i-i4D) ÷ NN + 1
    return @inbounds x.w[iL][i4D]
end

function substitute_fermion!(a::AbstractFermionfields_3D{NC},b::AbstractFermionfields_3D{NC}) where NC
    L5 = length(a.w)
    n1,n2, n5,n6 = size(a.w[1].f)
    for iL = 1:L5
        @inbounds for i6=1:n6
            for i5=1:n5
                #for i4=1:n4
                    #for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:NC
                                a.w[iL].f[i1,i2, i5,i6]= b.w[iL].f[i1,i2, i5,i6]
                            end
                        end
                    #end
                #end
            end
        end
    end
end

function clear_fermion!(a::Vector{<: AbstractFermionfields_3D{NC}}) where NC
    for μ=1:length(a)
        clear_fermion!(a[μ])
    end
end

function clear_fermion!(a::AbstractFermionfields_3D{NC}) where NC
    L5 = length(a.w)
    n1,n2, n5,n6 = size(a.w[1].f)
    for iL = 1:L5
        @inbounds for i6=1:n6
            for i5=1:n5
                #for i4=1:n4
                    #for i3=1:n3
                        for i2=1:n2
                            @simd for i1=1:NC
                                a.w[iL].f[i1,i2, i5,i6]= 0
                            end
                        end
                    #end
                #end
            end
        end
    end
end

function add_fermion!(c::AbstractFermionfields_3D{NC},α::Number,a::T1,β::Number,b::T2) where {NC,T1 <: Abstractfermion,T2 <: Abstractfermion}#c += alpha*a + beta*b
    L5 = length(a.w)
    for iL=1:L5
        add_fermion!(c.w[iL],α,a.w[iL],β,b.w[iL])
    end
end

function set_wing_fermion!(F::AbstractFermionfields_3D{NC})  where NC
    set_wing_fermion!(F,default_boundaryconditions[4])
end

function set_wing_fermion!(F::AbstractFermionfields_3D{NC},boundarycondition)  where NC
    L5 = length(F.w)
    for iL = 1:L5
        set_wing_fermion!(F.w[iL],boundarycondition)
    end
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::AbstractFermionfields_3D{NC}) where NC
    NX = x.NX
    #NY = x.NY
    #NZ = x.NZ
    NT = x.NT
    n6 = size(x.w[1].f)[end]
    σ = sqrt(1/2)

    L5 = length(x.w)
    for iL = 1:L5
        for ialpha = 1:n6
            for it=1:NT
                #for iz=1:NZ
                    #for iy=1:NY
                        for ix=1:NX
                            for ic=1:NC
                                x.w[iL][ic,ix, it,ialpha] = σ*randn()+im*σ*randn()
                            end
                        end
                    #end
                #end
            end
        end
        set_wing_fermion!(x.w[iL])
    end
    return
end

"""
c-------------------------------------------------c
c     Random number function for Gaussian  Noise
    with σ^2 = 1/2
c-------------------------------------------------c
    """
function gauss_distribution_fermion!(x::AbstractFermionfields_3D{NC},randomfunc,σ) where NC

    NX = x.NX
    #NY = x.NY
    #NZ = x.NZ
    NT = x.NT
    n6 = size(x.w[1].f)[end]
    #σ = sqrt(1/2)

    L5 = length(x.w)
    @inbounds for iL = 1:L5

        for mu = 1:n6
            for ic=1:NC
                for it=1:NT
                    #for iz=1:NZ
                        #for iy=1:NY
                            for ix=1:NX
                                v1 = sqrt(-log(randomfunc()+1e-10))
                                v2 = 2pi*randomfunc()

                                xr = v1*cos(v2)
                                xi = v1 * sin(v2)

                                x.w[iL][ic,ix, it,mu] = σ*xr + σ*im*xi
                            end
                        #end
                    #end
                end
            end
        end

        set_wing_fermion!(x.w[iL])
    end

    return
end

function gauss_distribution_fermion!(x::AbstractFermionfields_3D{NC},randomfunc) where NC
    σ = 1
    gauss_distribution_fermion!(x,randomfunc,σ)
end

#Overwrite Y with X*a + Y*b, where a and b are scalars. Return Y.
function LinearAlgebra.axpby!(a::Number, X::T, b::Number, Y::AbstractFermionfields_3D{NC}) where {NC,T <: AbstractFermionfields_3D}
    L5 = length(X.w)
    #n1,n2, n5,n6 = size(X.w[1].f)
    for iL = 1:L5
        Y.w[iL] = axpby!(a,X.w[iL],b,Y.w[iL])
    end
    return Y
end

function LinearAlgebra.axpy!(a::Number, X::T, Y::AbstractFermionfields_3D{NC}) where {NC,T <: AbstractFermionfields_3D}
    LinearAlgebra.axpby!(a,X,1,Y)
    return Y
end

function LinearAlgebra.dot(a::AbstractFermionfields_3D{NC},b::AbstractFermionfields_3D{NC}) where {NC}
    L5 = length(a.w)
    c = 0.0im
    for iL = 1:L5
        c += dot(a.w[iL],b.w[iL])
    end
    return c
end


text:
#---------------------------------------------------------------------#
#This routine advances the solution in time using Explicit RK methods
#written in Butcher tableau form.
#Written by F.X. Giraldo / P.R. Mugg on 9/22/19
#           Department of Applied Mathematics
#           Naval Postgraduate School
#           Monterey, CA 93943-5216
#---------------------------------------------------------------------#
function erk_butcher(q0,c,αI,βI,I,Δt,time)

   Q = zeros(I,1)
   R = zeros(I,1)

   c𝑖 = zeros(I)
   for i = 1:I
      c𝑖[i] = sum(αI[i,:])
   end

   Q[1] = q0
   t = time
   (rhs,f,g) = rhs_function(Q[1],c,t)
   R[1] = rhs

   for i = 2:I
      R_sum = 0
      for j = 1:i-1
         R_sum += αI[i,j]*R[j]
      end
      Q[i] = q0 + Δt*R_sum
      t = time + Δt*c𝑖[i]
      (rhs,f,g) = rhs_function(Q[i],c,t)
      R[i] = rhs
   end

   R_sum = 0
   for i = 1:I
      R_sum += βI[i]*R[i]
   end

   return qp = q0 + Δt*R_sum
end #end function


text:
while 42
	println("[Deep Thought] Answer to the Ultimate Question of Life, the Universe, and Everything")
end

max_stars_repo_name: gabriel-doriath-dohler/pjulia

max_issues_repo_name: gabriel-doriath-dohler/pjulia

max_forks_repo_name: gabriel-doriath-dohler/pjulia


text:
@testset "Smoothing" begin
  @testset "TaubinSmoothing" begin
    mesh  = readply(T, joinpath(datadir,"beethoven.ply"))
    smesh = smooth(mesh, TaubinSmoothing(30))
    # smoothing doesn't change the topology
    @test nvertices(smesh) == nvertices(mesh)
    @test nelements(smesh) == nelements(mesh)
    @test topology(smesh) == topology(mesh)
  end
end


text:
@inline @inbounds function henon(x, p, n) # In this case, p is a tuple of functions
    dx1 = 1.0 - p(n)[1]*x[1]^2 + x[2]
    dx2 = p(n)[2]*x[1]
    return SVector{2}(dx1, dx2)
end
# Jacobian:
@inline @inbounds function henon_J(x, p, n)
    J = @SMatrix [-2*p(n)[1]*x[1]      1.0;
                p(n)[2]            0.0]
    return J
end

"""
Henon map
"""
function henon(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters), jac=henon_J)
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

henonSim = Process(
    process = henon,
    X0 = [0.0, 0.0],
    parameter_profile = (constantParameter, constantParameter),
    parameter_profile_parameters = (1.4, 0.3),
    transient_t0 = 0,
    t0 = 1000,
    dt = 1,
    savedt = 1,
    tmax = 6000,
    alg = FunctionMap()) # The only discrete solver, pretty much
export henonSim

@inline @inbounds function logistic(x, 𝑟::Function, n)
    dx = 𝑟(n)*x*(1-x[1])
    return SVector{1}(dx)
end

"""Logistic map"""
function logistic(P::Process)
    seed(P.solver_rng)
    prob = ODEProblem(P.process, P.X0, (P.transient_t0, P.tmax), tuplef2ftuple(P.parameter_profile, P.parameter_profile_parameters))
    sol = dsolve(prob, P.alg; dt = P.dt, saveat=P.savedt, P.solver_opts...)
end

logisticSim = Process(
    process = logistic,
    X0 = [0.1],
    parameter_profile = constantParameter,
    parameter_profile_parameters = (4.0,), # (threshold, baseline, stepHeight)
    transient_t0 = 0,
    t0 = 1000,
    dt = 1,
    savedt = 1,
    tmax = 6000,
    alg = FunctionMap()) # The only discrete solver, pretty much
export logisticSim

"""
Ikeda map
"""
function ikeda(x, p, n)
    X, Y = x
    𝛼, 𝛽, 𝛾, 𝜇 = p(n)
    𝜙 = 𝛽 - 𝛼/(1 + X^2 + Y^2)
    dX = 𝛾 + 𝜇*(X*cos(𝜙) - Y*sin(𝜙))
    dY = 𝜇*(X*sin(𝜙) + Y*cos(𝜙))
    return SVector{2}(dX, dY)
end

ikedaSim = Process(
    process = ikeda(P) = process2solution(P),
    X0 = [0.0, 0.0],
    parameter_profile = ntuple(x->constantParameter, 4),
    parameter_profile_parameters = (17, 0.2, 0.9, 0.85), # [6.0, 0.4, 1.0, 0.9], # Sprott
    transient_t0 = 0,
    t0 = 1000,
    dt = 1,
    savedt = 1,
    tmax = 11000,
    alg = FunctionMap())
export ikedaSim

ikedaVis = ikedaSim(tmax = 200000)
export ikedaVis