text:
using TaylorIntegration, Test, DiffEqBase
using LinearAlgebra: norm
using StaticArrays
using Logging
import Logging: Warn

@testset "Testing `common.jl`" begin

    local TI = TaylorIntegration
    max_iters_reached() = "Maximum number of integration steps reached; exiting.\n"

    f(u,p,t) = u
    g(u,p,t) = cos(t)
    f!(du, u, p, t) = (du .= u)

    @testset "Test integration of ODE with numbers in common interface" begin
        u0 = 0.5
        tspan = (0.0, 1.0)
        prob = ODEProblem(f, u0, tspan)
        sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(50), abstol=1e-20))
        @test TI.alg_order(TaylorMethod(50)) == 50
        @test TI.alg_order(sol.alg) == TI.alg_order(TaylorMethod(50))
        @test abs(sol[end] - u0*exp(1)) < 1e-12
        u0 = 0.0
        tspan = (0.0, 11pi)
        prob = ODEProblem(g, u0, tspan)
        sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(50), abstol=1e-20))
        @test abs(sol[end] - sin(sol.t[end])) < 1e-12
    end

    # @testset "Test integration of ODE with abstract arrays in common interface" begin
    #     u0 = rand(4, 2)
    #     tspan = (0.0, 1.0)
    #     prob = ODEProblem(f!, u0, tspan)
    #     sol = solve(prob, TaylorMethod(50), abstol=1e-20)
    #     @test norm(sol[end] - u0.*exp(1)) < 1e-12
    #     f_oop(u, p, t) = u
    #     prob_oop = ODEProblem(f_oop, u0, tspan)
    #     sol_oop = solve(prob_oop, TaylorMethod(50), abstol=1e-20)
    #     @test norm(sol_oop[end] - u0.*exp(1)) < 1e-12
    #     @test sol.t == sol_oop.t
    #     @test sol.u == sol_oop.u
    # end

    local tspan = (0.0,5.0)
    local saveat_inputs = ([], 0:1:(tspan[2]+5), 0:1:tspan[2], 3:1:tspan[2], collect(0:1:tspan[2]))
    @testset "Test saveat behavior with numbers in common interface" begin
        u0 = 1.0
        prob = ODEProblem(f, u0, tspan)
        sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20))
        s = saveat_inputs[1]
        sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s))
        @test all(sol.t .== sol_taylor.t)
        @test all(sol_taylor.u .== sol.u)
        @test length(sol_taylor.t) == length(sol_taylor.u)
        s = saveat_inputs[2]
        sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s))
        @test sol_taylor.t == saveat_inputs[3]
        for s ∈ saveat_inputs[3:end]
            sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s))
            @test all(s .== sol_taylor.t)
            @test length(sol_taylor.t) == length(sol_taylor.u)
        end
    end

    @testset "Test saveat behavior with abstract arrays in common interface" begin
        u0 = rand(2)
        prob = ODEProblem(f!, u0, tspan)
        sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20))
        s = saveat_inputs[1]
        sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s))
        @test all(sol.t .== sol_taylor.t)
        @test all(sol_taylor.u .== sol.u)
        @test length(sol_taylor.t) == length(sol_taylor.u)
        s = saveat_inputs[2]
        sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s))
        @test sol_taylor.t == saveat_inputs[3]
        for s ∈ saveat_inputs[3:end]
            sol_taylor = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(20), abstol=1e-20, saveat = s))
            @test all(s .== sol_taylor.t)
            @test length(sol_taylor.t) == length(sol_taylor.u)
        end
    end

    function harmosc!(dx, x, p, t)
        dx[1] = x[2]
        dx[2] = - x[1]
        return nothing
    end
    @testset "Test consistency with taylorinteg" begin
        tspan = (0.0, 10pi)
        abstol = 1e-20
        order = 25 # Taylor expansion order wrt time
        u0 = [1.0; 0.0]
        prob = ODEProblem(harmosc!, u0, tspan)
        sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(order), abstol=abstol))
        tv1, xv1 = (@test_logs min_level=Logging.Warn taylorinteg(
            harmosc!, u0, tspan[1], tspan[2], order, abstol))
        @test sol.t == tv1
        @test xv1[end,:] == sol[end]
    end

    @testset "Test discrete callback in common interface" begin
        # discrete callback: example taken from DifferentialEquations.jl docs:
        # https://diffeq.sciml.ai/dev/features/callback_functions/#Using-Callbacks
        tspan = (0.0, 10pi)
        abstol = 1e-20
        order = 25 # Taylor expansion order wrt time
        u0 = [1.0; 0.0]
        prob = ODEProblem(harmosc!, u0, tspan)
        t_cb = 1.0pi
        condition(u,t,integrator) = t == t_cb
        affect!(integrator) = integrator.u[1] += 0.1
        cb = DiscreteCallback(condition,affect!)
        sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(order), abstol=abstol, tstops=[t_cb], callback=cb))
        @test sol.t[4] == t_cb
        @test sol.t[4] == sol.t[5]
        @test sol[4][1] + 0.1 == sol[5][1]
    end

    @testset "Test continuous callback in common interface" begin
        # continuous callback: example taken from DifferentialEquations.jl docs:
        # https://diffeq.sciml.ai/stable/features/callback_functions/#Example-1:-Bouncing-Ball
        @taylorize function f(du,u,p,t)
            local g_acc = p
            du[1] = u[2]
            du[2] = -g_acc + zero(u[2])
        end
        function condition(u,t,integrator)
            u[1]
        end
        function affect!(integrator)
            integrator.u[2] = -integrator.u[2]
            return nothing
        end
        cb = ContinuousCallback(condition,affect!)
        u0 = [50.0,0.0]
        tspan = (0.0,15.0)
        p = 9.8
        prob = ODEProblem(f,u0,tspan,p)
        # Avoid log-checks for Julia versions before v1.6
        if VERSION < v"1.6"
            sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb)
        else
            sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb))
        end
        tb = sqrt(2*50/9.8) # bounce time
        @test abs(tb - sol.t[9]) < 1e-14
        @test sol.t[9] == sol.t[10]
        @test sol[9][1] == sol[10][1]
        @test sol[9][2] == -sol[10][2] # check that callback was applied correctly (1st bounce)
        @test abs(3tb - sol.t[19]) < 1e-14
        @test sol.t[19] == sol.t[20]
        @test sol[19][1] == sol[20][1]
        @test sol[19][2] == -sol[20][2] # check that callback was applied correctly (2nd bounce)
    end

    @testset "Test vector continuous callback in common interface" begin
        # vector continuous callback: example taken from DifferentialEquations.jl docs:
        # https://diffeq.sciml.ai/dev/features/callback_functions/#VectorContinuousCallback-Example
        @taylorize function ff(du,u,p,t)
            local g_acc = p
            du[1] = u[2]
            du[2] = -g_acc + zero(u[2])
            du[3] = u[4]
            du[4] = zero(u[4])
        end

        function condition(out,u,t,integrator) # Event when event_f(u,t) == 0
            out[1] = u[1]
            out[2] = (u[3] - 10.0)u[3]
        end

        function affect!(integrator, idx)
            if idx == 1
                integrator.u[2] = -0.9integrator.u[2]
            elseif idx == 2
                integrator.u[4] = -0.9integrator.u[4]
            end
        end

        cb = VectorContinuousCallback(condition, affect!, 2)

        u0 = [50.0, 0.0, 0.0, 2.0]
        tspan = (0.0, 15.0)
        p = 9.8
        prob = ODEProblem(ff, u0, tspan, p)
        # Avoid log-checks for Julia versions before v1.6
        if VERSION < v"1.6"
            sol = solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb)
        else
            sol = (@test_logs min_level=Logging.Warn solve(prob, TaylorMethod(25), abstol=1e-16, callback=cb))
        end
        tb = sqrt(2*50/9.8) # bounce time
        @test abs(tb - sol.t[8]) < 1e-14
        @test sol.t[8] == sol.t[9]
        @test sol[9][1] == sol[8][1]
        @test sol[9][2] == -0.9sol[8][2]
        @test sol[9][3] == sol[8][3]
        @test sol[9][4] == sol[8][4]
        @test sol.t[13] == sol.t[14]
        @test sol[14][1] == sol[13][1]
        @test sol[14][2] == sol[13][2]
        @test sol[14][3] == sol[13][3]
        @test sol[14][4] == -0.9sol[13][4]
    end

    @testset "Test parsed jetcoeffs! method in common interface" begin
        b1 = 3.0
        order = 25 # Taylor expansion order wrt time
        abstol = 1.0e-20
        @taylorize xdot1(x, p, t) = b1-x^2
        @test (@isdefined xdot1)
        x0 = 1.0
        tspan = (0.0, 1000.0)
        prob = ODEProblem(xdot1, x0, tspan)
        sol1 = (@test_logs min_level=Logging.Warn solve(
            prob, TaylorMethod(order), abstol=abstol, parse_eqs=false))
        @test sol1.alg.parse_eqs == false
        sol2 = (@test_logs min_level=Logging.Warn solve(
            prob, TaylorMethod(order), abstol=abstol))
        @test sol2.alg.parse_eqs == true
        @test sol1.t == sol2.t
        @test sol1.u == sol2.u

        @taylorize function integ_vec(dx, x, p, t)
            local λ = p[1]
            dx[1] = cos(t)
            dx[2] = -λ*sin(t)
            return dx
        end
        @test (@isdefined integ_vec)
        x0 = [0.0, 1.0]
        tspan = (0.0, 11pi)
        prob = ODEProblem(integ_vec, x0, tspan, [1.0])
        sol1 = (@test_logs min_level=Logging.Warn solve(
            prob, TaylorMethod(order), abstol=abstol, parse_eqs=false))
        sol2 = (@test_logs min_level=Logging.Warn solve(
            prob, TaylorMethod(order), abstol=abstol)) # parse_eqs=true
        @test sol2.alg.parse_eqs == true
        @test length(sol1.t) == length(sol2.t)
        @test sol1.t == sol2.t
        @test sol1.u == sol2.u
        tv, xv = (@test_logs min_level=Logging.Warn taylorinteg(
            integ_vec, x0, tspan[1], tspan[2], order, abstol, [1.0]))
        @test sol1.t == tv
        @test sol1[1,:] == xv[:,1]
        @test sol1[2,:] == xv[:,2]
        @test norm(sol1[end][1] - sin(sol1.t[end]), Inf) < 1.0e-14
        @test norm(sol1[end][2] - cos(sol1.t[end]), Inf) < 1.0e-14

        V(x, y) = - (1-μ)/sqrt((x-μ)^2+y^2) - μ/sqrt((x+1-μ)^2+y^2)
        H_pcr3bp(x, y, px, py) = (px^2+py^2)/2 - (x*py-y*px) + V(x, y)
        H_pcr3bp(x) = H_pcr3bp(x...)
        J0 = -1.58 # Jacobi constant
        q0 = [-0.8, 0.0, 0.0, -0.6276410653920694] # initial condition
        μ = 0.01 # mass parameter
        tspan = (0.0, 1000.0)
        p = [μ]
        @taylorize function pcr3bp!(dq, q, param, t)
            local μ = param[1]
            local onemμ = 1 - μ
            x1 = q[1]-μ
            x1sq = x1^2
            y = q[2]
            ysq = y^2
            r1_1p5 = (x1sq+ysq)^1.5
            x2 = q[1]+onemμ
            x2sq = x2^2
            r2_1p5 = (x2sq+ysq)^1.5
            dq[1] = q[3] + q[2]
            dq[2] = q[4] - q[1]
            dq[3] = (-((onemμ*x1)/r1_1p5) - ((μ*x2)/r2_1p5)) + q[4]
            dq[4] = (-((onemμ*y )/r1_1p5) - ((μ*y )/r2_1p5)) - q[3]
            return nothing
        end
        prob = ODEProblem(pcr3bp!, q0, tspan, p)
        sol1 = (@test_logs min_level=Logging.Warn solve(
            prob, TaylorMethod(order), abstol=abstol))
        @test sol1.alg.parse_eqs == true
        sol2 = (@test_logs min_level=Logging.Warn solve(
            prob, TaylorMethod(order), abstol=abstol, parse_eqs=false))
        @test sol2.alg.parse_eqs == false
        @test norm( H_pcr3bp(sol1.u[end]) - J0 ) < 1e-10
        @test norm( H_pcr3bp(sol2.u[end]) - J0 ) < 1e-10
        @test sol1.u == sol2.u
        @test sol1.t == sol2.t
    end

    @testset "Test throwing errors in common interface" begin
        # NOTE: The original test defines u0 as a matrix (rand(4,2)) but that doesn't work
        u0 = rand(4)#rand(4, 2)
        tspan = (0.0, 1.0)
        prob1 = ODEProblem(f!, u0, tspan)
        sol = (@test_logs min_level=Logging.Warn solve(prob1, TaylorMethod(50), abstol=1e-20))

        # `order` is not specified
        @test_throws ErrorException solve(prob, TaylorMethod(), abstol=1e-20)
    end

    ### DynamicalODEProblem tests (see #108, #109)
    @testset "Test integration of DynamicalODEPoblem" begin
        function iip_q̇(dq,p,q,params,t)
            dq[1] = p[1]
            dq[2] = p[2]
        end

        function iip_ṗ(dp,p,q,params,t)
            dp[1] = -q[1] * (1 + 2q[2])
            dp[2] = -q[2] - (q[1]^2 - q[2]^2)
        end

        iip_q0 = [0.1, 0.]
        iip_p0 = [0., 0.5]

        function oop_q̇(p, q, params, t)
            p
        end

        function oop_ṗ(p, q, params, t)
            dp1 = -q[1] * (1 + 2q[2])
            dp2 = -q[2] - (q[1]^2 - q[2]^2)
            @SVector [dp1, dp2]
        end

        oop_q0 = @SVector [0.1, 0.]
        oop_p0 = @SVector [0., 0.5]

        T(p) = 1//2 * (p[1]^2 + p[2]^2)
        V(q) = 1//2 * (q[1]^2 + q[2]^2 + 2q[1]^2 * q[2]- 2//3 * q[2]^3)
        H(p,q, params) = T(p) + V(q)

        E = H(iip_p0, iip_q0, nothing)

        energy_err(sol) = maximum(i->H([sol[1,i], sol[2,i]], [sol[3,i], sol[4,i]], nothing)-E, 1:length(sol.u))

        iip_prob = DynamicalODEProblem(iip_ṗ, iip_q̇, iip_p0, iip_q0, (0., 100.))
        oop_prob = DynamicalODEProblem(oop_ṗ, oop_q̇, oop_p0, oop_q0, (0., 100.))

        sol1 = (@test_logs min_level=Logging.Warn solve(iip_prob, TaylorMethod(50), abstol=1e-20))
        @test energy_err(sol1) < 1e-10

        sol2 = (@test_logs min_level=Logging.Warn solve(oop_prob, TaylorMethod(50), abstol=1e-20))
        @test energy_err(sol2) < 1e-10

        @test sol1.t == sol2.t
        @test sol1.u[:] == sol2.u[:]
    end

end


text:
using Pkg
Pkg.activate(".")

using Images, CausticsEngineering

image = Images.load("./examples/cat_posing.jpg") # Check current working directory with pwd()
engineer_caustics(image);

ext: jl

max_stars_repo_name: MattFerraro/causticsEngineering


text:
module CachedArrays

#####
##### Exports
#####

export CachedArray,
    ReadableCachedArray, UnreadableCachedArray, WritableCachedArray, UnwritableCachedArray
export @annotate, tocached
export onobjects, @blockobjects, findobjects, @noescape

# state transitions
export readable, writable, release

#####
##### Deps
#####

# base
import Base: @lock

# stdlib
import Dates
import LinearAlgebra
import Mmap
import Random

# Dependencies
# Import "constructorof" to allow external modules to extend "constructorof" without
# directly invoking "ConstructablesBase'.
# Of course, that's what's ultimately happening, but this is just cleaner.
import ConstructionBase: ConstructionBase, constructorof
import DataStructures
import JSON
import Polyester
import SIMD
import MacroTools
import TimerOutputs

using DocStringExtensions

const DEBUG = true              # Enable asserts
const VERBOSE = false           # Diagnostic printing
const ENABLETIMING = true       # Debug timing
# const ALLOW_UNSAFE_FREE = true  # Enable the `unsafe_free` function.
include("options.jl")

#####
##### includes
#####

@enum AllocationPriority ForceLocal PreferLocal ForceRemote

# Bootstrap Utilities
include("utils/utils.jl")
include("utils/findnexttree.jl")
include("utils/freebuffer.jl")
include("utils/lru.jl")
include("object.jl")
include("allocators.jl")

# Heap implementations
include("memory/memory.jl")

# Cache eviction policies
include("policy.jl")
include("policies/local.jl")

# Implementation of the arrays and cache manager
include("manager.jl")
include("validation.jl")

# Array Implementstions
include("llvm.jl")
using .LoadStore: LoadStore

include("array.jl")
include("telemetry.jl")

# Fast "memcpy"
include("memcpy.jl")
include("lib.jl")
include("api.jl")

#####
##### Keep track of Managers
#####

# When Managers get created - hold them in a global array to keep them from getting
# GC'd before the arrays they track.
#
# We can periodically GC the mangers to see if they are no longer holding onto anything,
# at which point we're free to reclaim their resources.
const GlobalManagers = CacheManager[]
function gc_managers()
    # Find all the managers that can be garbage collected.
    #
    # Trigger a full garbage collection before to clean up anything that may be holding
    # onto a manager
    GC.gc(true)

    # Now, find all the managers that have been completely cleaned up.
    i = findall(cangc, GlobalManagers)
    i === nothing && return 0

    # Remove the managers from this list and then run a full garbage collection to make sure
    # they're well and completely gone.
    deleteat!(GlobalManagers, i)
    GC.gc(true)

    # Return the number of managers we removed
    return length(i)
end

end # module


text:
#=
# Based on the file JacobiGQ.m from the book:
# Nodal Discontinuous Galerkin Method - Hesthaven, Warburton
#  https://link.springer.com/book/10.1007/978-0-387-72067-8
=#

# Compute the Nth order Gauss quadrature points, x, and weights, w,
# associated with the Jacobi polynomial of type (alpha,beta).
function jacobi_gauss_quad(alpha::Int, beta::Int, N::Int)
    if N==0
        x = [-(alpha-beta)/(alpha+beta+2)];
        w = [2];
        return (x,w);
    end

    function sqrtvec(v::Array)
        w = zeros(size(v));
        for i=1:length(v)
            w[i] = sqrt(v[i]);
        end
        return w;
    end

    # Form symmetric matrix from recurrence.
    #J = zeros(N+1,N+1);
    h1 = 2*(0:N).+(alpha+beta);
    J = diagm(-1/2*(alpha^2-beta^2)./(h1.+2)./h1) +
        diagm(1 => 2 ./(h1[1:N].+2).*sqrtvec((1:N).*((1:N).+(alpha+beta)).*
        ((1:N).+alpha).*((1:N).+beta)./(h1[1:N].+1)./(h1[1:N].+3)));
    if (alpha+beta == 0)
        J[1,1]=0.0;
    end
    J = J + J';

    # Find eigenvalues and vectors
    E = eigen(J);
    #x = E.values;
    w = (E.vectors[1,:]').^2 .*2^(alpha+beta+1)/(alpha+beta+1)*factorial(alpha)*
        factorial(beta)/factorial(alpha+beta);
    w = w[:];
    return (E.values,w);
end

function jacobi_LGL_quad(N::Int)
    if N == 1
        lgl = [-1; 1];
        w = [1; 1];
    else
        # compute the points
        (r,w) = jacobi_gauss_quad(1,1,N-2);
        lgl = [-1; r ; 1];

        # compute the weights
        w = jacobi_polynomial(lgl, 0, 0, N);
        adgammaN = (2*N + 1) / (N * (N + 1));
        w = w.*w;
        w = adgammaN./w;
    end

    return (lgl, w);
end


text:
using LinearAlgebra
using StatsFuns
using Distances

function ∂ofK!(∂K, C, cK)
    ret = - dot(∂K * C, C) / 2
    ldiv!(cK, ∂K)
    ret -= tr(∂K) / 2
    return ret
end

function compute_cov!(X, K, ∂K, ℓ, ϵ, amp, dind)
    copyto!(X, X)
    rdiv!(X, Diagonal(ℓ))
    pairwise!(K, SqEuclidean(), X; dims=1)
    K .= exp.(.-K)
    if !isnothing(∂K)
        copyto!(∂K, K)
    end
    K .= amp .* K
    for i in dind
        K[i] += ϵ
    end
    cK = cholesky!(Symmetric(K))
end

function logprob_gen(X, Y)
    n = length(Y)
    L = similar(X)
    K = similar(X, n, n)
    Z = log.(Y)
    lnJ = -sum(Z)
    C = similar(Y)
    dind = diagind(K)
    f = function (ℓ, ϵ, amp)
        copyto!(L, X)
        cK = compute_cov!(L, K, nothing, ℓ, ϵ, amp, dind)
        lnd = logdet(cK)
        ldiv!(C, cK, Z)
        return lnJ - lnd / 2 - dot(Z, C) / 2 - n * log2π / 2
    end
    ∂_K = similar(K)
    ∂ℓ = similar(X, size(X, 2))
    fg = function (ℓ, ϵ, amp)
        copyto!(L, X)
        cK = compute_cov!(L, K, ∂_K, ℓ, ϵ, amp, dind)
        lnd = logdet(cK)
        ldiv!(C, cK, Z)
        ∂amp = ∂ofK!(∂_K, C, cK)
        for i in 1:length(ℓ)
            pairwise!(∂_K, SqEuclidean(), view(X, :, i:i); dims=1)
            ∂ℓ[i] = ∂ofK!(∂_K, C, cK)
        end
        copyto!(∂_K, I)
        ∂ϵ = ∂ofK!(∂_K, C, cK)
        ∂_K = similar(K)
        ∂ℓ = similar(X, size(X, 2))
        ret = lnJ - lnd / 2 - dot(Z, C) / 2 - n * log2π / 2
        return ret, ∂ℓ, ∂ϵ, ∂amp
    end
    return f, fg
end

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

# BEGIN 0.7 deprecations

# PR #22062
function set_remote_url(repo::LibGit2.GitRepo, url::AbstractString; remote::AbstractString="origin")
    Base.depwarn(string(
        "`LibGit2.set_remote_url(repo, url; remote=remote)` is deprecated, use ",
        "`LibGit2.set_remote_url(repo, remote, url)` instead."), :set_remote_url)
    set_remote_url(repo, remote, url)
end

function set_remote_url(path::AbstractString, url::AbstractString; remote::AbstractString="origin")
    Base.depwarn(string(
        "`LibGit2.set_remote_url(path, url; remote=remote)` is deprecated, use ",
        "`LibGit2.set_remote_url(path, remote, url)` instead."), :set_remote_url)
    set_remote_url(path, remote, url)
end

function prompt(msg::AbstractString; default::AbstractString="", password::Bool=false)
    Base.depwarn(string(
        "`LibGit2.prompt(msg::AbstractString; default::AbstractString=\"\", password::Bool=false)` is deprecated, use ",
        "`result = Base.prompt(msg, default=default, password=password); result === nothing ? \"\" : result` instead."), :prompt)
    coalesce(Base.prompt(msg, default=default, password=password), "")
end

# PR #26437
# when this deprecation is deleted, remove all calls to it, and remove the keyword of:
# `payload` from "src/LibGit2.jl"
function deprecate_payload_keyword(f, sig, payload)
    if payload !== nothing
        Base.depwarn(string(
            "`LibGit2.$f($sig; payload=cred)` is deprecated, use ",
            "`LibGit2.$f($sig; credentials=cred)` instead."), f)
    end
end

@deprecate get_creds!(cache::CachedCredentials, credid, default) get!(cache, credid, default)

@eval Base @deprecate merge!(repo::$(GitRepo), args...; kwargs...) $(LibGit2.merge!)(repo, args...; kwargs...)
@eval Base @deprecate push!(w::$(GitRevWalker), arg) $(LibGit2.push!)(w, arg)
@eval Base @deprecate count(diff::$(GitDiff)) $(LibGit2.count)(diff)
@eval Base @deprecate count(idx::$(GitIndex)) $(LibGit2.count)(idx)
@eval Base @deprecate count(rb::$(GitRebase)) $(LibGit2.count)(rb)
@eval Base @deprecate count(tree::$(GitTree)) $(LibGit2.count)(tree)
@eval Base @deprecate count(f::Function, walker::$(GitRevWalker); kwargs...) $(LibGit2.count)(f, walker; kwargs...)
@eval Base @deprecate map(f::Function, walker::$(GitRevWalker); kwargs...) $(LibGit2.map)(f, walker; kwargs...)

# PR #24594
@deprecate AbstractCredentials AbstractCredential false
@deprecate UserPasswordCredentials UserPasswordCredential false
@deprecate SSHCredentials SSHCredential false

@deprecate hex(id::LibGit2.GitHash)      string(id)
@deprecate hex(id::LibGit2.GitShortHash) string(id)

max_stars_repo_name: ninjin/julia

max_issues_repo_name: ninjin/julia

max_forks_repo_name: ninjin/julia


text:
__precompile__()
#
module sdelab_gq
#
export  sde_gq_solution, sde_gq_solution_extrap, set_fcn_mf
#
import NLsolve
import ForwardDiff
import Roots
#
using gauss_quadrature
#import gauss_quadrature
#
##########################
function init(fcn, opt)
  #=
  =#
end
############################
#

type Fcnmf
  #
  drift # u -> f(u) for u (n+1)-dimensions, f(u) scalar
  drift_ext # (z,q)-> q ( f(z,r(.))) integrate out 2nd arg
  Jdrift # Jacobian of drift (auto computed)
  Jdrift_ext # (z,q) -> q (Jdrift(z,r))
  Hdrift # Hessian of drift (auto computed)
  #
  diff # u-> g(u) diffusion for u (n+1)-dimensions, g(u) scalar
  diff_ext # (z,q)-> q ( g(z,r(.))) integrate out 2nd arg
  Jdiff # Jacobian of diffuion (auto computed)
  Hdiff # Hessian of diffusion (auto computed)
  #
  mf # mean-field fn x->r(x) for x scalar, r(mf) n-dimensional
  #
  la # scipt(L)a term in notes
  #
end
#
function set_fcn_mf(drift,diff,r)
  #=
  Define mean-field sde via Fcnmf data-type
  =#
  drift_ext=(z,q)->gauss_quadrature.do_quad(y->drift([z;r(y)]),q)
  diff_ext =(z,q)->gauss_quadrature.do_quad(y->diff( [z;r(y)]),q)
  # drift (n+1)-dim to 1-dim
  Jdrift=ForwardDiff.gradient(drift) # Jacobian (n+1)-dim -> (n+1)-dim
  # only need (1,1) derivative
  drift1=u->ForwardDiff.derivative(x->drift([x;u[2:end]]),u[1])
  Jdrift_ext=(z,q)->gauss_quadrature.do_quad(y->drift1([z;r(y)]),q)
  drift2=u->ForwardDiff.derivative(x->drift1([x;u[2:end]]),u[1])
  Hdrift=drift2
  #Hdrift=ForwardDiff.hessian(drift)
  #
  Jdiff=ForwardDiff.gradient(diff) # drift (n+1)-dim -> (n+1)-dim
  Hdiff=ForwardDiff.hessian(diff)
  diff1=u->ForwardDiff.derivative(x->diff([x;u[2:end]]),u[1])
  diff2=u->ForwardDiff.derivative(x->diff1([x;u[2:end]]),u[1])
  Hdiff=diff2
  #
  dr=ForwardDiff.derivative(r)
  ddr=ForwardDiff.derivative(dr)
  #
  la= (x,mf)->(dr(x)*drift([x;mf])
              + 0.5*ddr(x)*diff([x;mf])^2)
  #
  return Fcnmf(drift,drift_ext,Jdrift,Jdrift_ext,Hdrift,
               diff, diff_ext, Jdiff, Hdiff,
               r,la)
end
#######################
function sde_gq_solution(fcn::Fcnmf, tspan, quad0::QUAD, params)
  #
  # Apply QG integration.
  #
  # Inputs:
  # fcn   = definition of MF-SDE
  # tspan = vector of times to report
  # quad0 = initial quadrature
  # params= dict of parameters
  #
  # Outputs: quadrature and
  # stats=see below (to analyse behaviour of method)
  #
  T=maximum(tspan) # get final time
  noSteps=round(Int32,T/params["MaxStepSize"])
  t=(0:noSteps)/noSteps # grid of time points
  dt=t[2]-t[1] # fixed timestep
  #
  quad=deepcopy(quad0)
  # initialise
  solver=params["Solver"]
  if solver=="MFEM0"
    params["order"]=1
    quad_dw=gauss_quadrature.Gauss_Hermite([0,dt], 2)
    get= q->get_1st(fcn, q, dt, quad_dw, params)
  elseif solver=="MFEM0e"
    params["order"]=1
    quad_dw=gauss_quadrature.Gauss_Hermite([0,dt], 2)
    get= q->get_1st_ext(fcn, q, dt, quad_dw, params)
  elseif solver=="MFEM2"
    params["order"]=2
    quad_dw=gauss_quadrature.Gauss_Hermite([0,dt], 3)
    dtmat=dt_order2(dt,quad_dw)
    get= q->get_2nd(fcn, q, dtmat, quad_dw, params)
  elseif solver=="MFEM"
    params["order"]=1
    quad_dw=gauss_quadrature.Gauss_Hermite([0,dt], 2)
    function nonlin!(un_,up,out,r)
      out[:] =(un_-up-params["Alpha"]*fcn.drift(vec([un_[1];r]))*dt) # must avoid recreating out!
    end
    function jnonlin!(un_,out,r)
      tmp=fcn.Jdrift(vec([un_[1];r]) )[1]
      out[:,:]=1-params["Alpha"]*tmp
    end
    get= q->get_1st_imp(fcn, q, dt, quad_dw, params,nonlin!,jnonlin!)
  elseif solver=="MFEMe"
    params["order"]=1
    quad_dw=gauss_quadrature.Gauss_Hermite([0,dt], 2)
    function nonlin!(un_,up,out,q::QUAD)
      out[:] =(un_-up-params["Alpha"]*fcn.drift_ext(un_[1],q)*dt) # must avoid recreating out!
    end
    function jnonlin!(un_,out,q::QUAD)
      tmp=fcn.Jdrift_ext(un_[1],q)
      out[:,:]=1-params["Alpha"]*tmp
    end
    get= q->get_1st_imp_ext(fcn, q, dt, quad_dw, params,nonlin!,jnonlin!)
  end
  # initialise
  no_pts=zeros(Int32,noSteps+1); X_min=zeros(noSteps+1); X_max=zeros(noSteps+1)
  X_min[1]=minimum(quad.X); X_max[1]=maximum(quad.X); no_pts[1]=quad.length
  m=quad.length
  factor=params["order"]^3
  for i=1:noSteps
    quad=get(quad)
    # update degree and recompute Gauss rule/threshold
    if quad.length>factor*m && i<noSteps-2
      m=no_Gauss_points(m, dt, T-t[i+1], params)
      if params["periodic"]==false # non-periodic
        centre_mass=gauss_quadrature.do_quad(x->x,quad)
        r=abs(log(dt))*params["cutoff"]
        cutoff= [centre_mass-r,centre_mass+r]
        if params["no_partition"]>1
          gauss_quadrature.partition_simplify!(quad, m, params["no_partition"], cutoff)
        else
          gauss_quadrature.simplify!(quad, m, cutoff)
        end
      else # periodic
        if params["no_partition"]>1
          gauss_quadrature.partition_periodic_simplify!(quad, m, params["no_partition"], params["length_scale"])
        else
          gauss_quadrature.periodic_simplify!(quad, m, params["length_scale"])
        end
      end # end if periodic
    end # reduce quad.length
    # stats
    X_min[i+1]=minimum(quad.X); X_max[i+1]=maximum(quad.X);
    no_pts[i+1]=quad.length;
  end # time loop
  gauss_quadrature.sortquad!(quad)
  stats=Dict()
  stats["min"]=X_min;  stats["max"]=X_max
  stats["no_pts"]=no_pts;  stats["t"]=t
  return quad, stats
end
##########################
function sde_gq_solution_extrap(fcn::Fcnmf, tspan, quad0::QUAD, params)
  # Computes the Talay-Tubaro extrapolation
  # Returns a rule with negative weights
  # Only effective for first-order integrator
  #
  # Same calling syntax as sde_pdf_solution
  quad,stats=sde_gq_solution(fcn, tspan, quad0, params)
  params1=deepcopy(params)
  params1["MaxStepSize"]=params["MaxStepSize"]/2
  quad1,stats1=sde_gq_solution(fcn, tspan, quad0, params1)
  x=[quad1.X; quad.X]; w=[2*quad1.w; -quad.w]
  new_quad=QUAD(length(x), x ,w)
  gauss_quadrature.sortquad!(new_quad)
  return new_quad, stats1
end
#############################################
function get_1st_ext(fcnmf::Fcnmf, quad_in::QUAD, dt, quad_dw::QUAD, params)
  # Inputs: fcnmf =  drift and diffusion functions
  # quad_in=current quadrature approx
  # dt=time step
  #
  # Outputs: quad_out the next timestep quad approx by Euler-Maruyama (explicit)
  #
  # same as get_1st, except mean-field calculation is completed OUTSIDE the function eval
  Y=zeros(quad_in.length,quad_dw.length)
  W=deepcopy(Y)
  drift1=z->fcnmf.drift_ext(z,quad_in)
  diff1 =z->fcnmf.diff_ext(z,quad_in)
  for j=1:quad_in.length
    x=quad_in.X[j]
    xpfx_dt = x + drift1(x)*dt
    gx      =     diff1(x)
    #
    Y[j,:]=xpfx_dt +   quad_dw.X * gx
    W[j,:]=quad_in.w[j]*quad_dw.w
  end
  if params["periodic"]==true # periodic
     Y=mod(Y,params["length_scale"])
   end
  quad_out=gauss_quadrature.QUAD(length(Y), vec(Y), vec(W) )
  return quad_out
end

#############################################
function get_1st(fcnmf::Fcnmf, quad_in::QUAD, dt, quad_dw::QUAD, params)
  # Inputs: fcnmf =  drift and diffusion functions
  # quad_in=current quadrature approx
  # dt=time step
  #
  # Outputs: quad_out the next timestep quad approx by Euler-Maruyama (explicit)
  #
  Y=zeros(quad_in.length,quad_dw.length)
  W=deepcopy(Y)
  mf=gauss_quadrature.do_quad(y->fcnmf.mf(y),quad_in)
  for j=1:quad_in.length
    x=quad_in.X[j]
    xpfx_dt = x + fcnmf.drift([x;mf])*dt
    gx=fcnmf.diff([x;mf])
    #
    Y[j,:]=xpfx_dt +   quad_dw.X * gx
    W[j,:]=quad_in.w[j]*quad_dw.w
  end
  if params["periodic"]==true # periodic
     Y=mod(Y,params["length_scale"])
   end
  quad_out=gauss_quadrature.QUAD(length(Y), vec(Y), vec(W) )
  return quad_out
end
##############################################

function get_2nd(fcnmf::Fcnmf, quad_in::QUAD, dtmat, quad_dw::QUAD, params)
  # Inputs: fcn =  drift and diffusion functions
  # quad_in=current quadrature approx by 2nd order
  #
  Y=zeros(quad_in.length, quad_dw.length)
  W=deepcopy(Y)
  mf=gauss_quadrature.do_quad(fcnmf.mf, quad_in)
  La=gauss_quadrature.do_quad(x->fcnmf.la(x,mf), quad_in)
  #
  for j=1:quad_in.length
    Y[j,:]=quad_in.X[j]+ coeffs_order2(quad_in.X[j], mf, La, fcnmf)*dtmat
    W[j,:]=quad_in.w[j]*quad_dw.w
  end
  #
  if params["periodic"]==true
     Y=mod(Y,params["length_scale"])
   end
  quad_out=gauss_quadrature.QUAD(length(Y), vec(Y), vec(W) )
  return quad_out
end
##########################################
#(fcn, q, dt, quad_dw, params,nonlin!,jnonlin!))
function get_1st_imp(fcnmf::Fcnmf, quad_in::QUAD, dt, quad_dw::QUAD, params, nonlin!, jnonlin!)
  # Inputs: fcn =  drift and diffusion functions
  # quad_in=current quadrature approx
  # dt=time step
  #
  # Outputs: quad_out the next timestep quad approx
  P=quad_in.length
  sqrt_dt=sqrt(dt)
  Y=zeros(P,1)
  Z=zeros(P,1)
  mf=gauss_quadrature.do_quad(fcnmf.mf, quad_in)
  helper=(y1,y2)->NLsolve.nlsolve(
    (u_,out)-> nonlin!(u_,y1, out,mf),                       (u_,out)->jnonlin!(u_,    out,mf),
    [y2])
  for j=1:P
    x=quad_in.X[j]
    xpfx_dt  = x + (1-params["Alpha"])*fcnmf.drift([x;mf])*dt
    gx_sqrth = fcnmf.diff([x;mf])*sqrt_dt
    #
    un=xpfx_dt+gx_sqrth
    solve=helper(un,x)
    Y[j]=solve.zero[1]
    #
    un=xpfx_dt-gx_sqrth

    solve=helper(un,x)
    Z[j]=solve.zero[1]
  end
  quad_out=QUAD(2*P, vec([Y Z]), [quad_in.w;quad_in.w]/2)
  return quad_out
end
##########################################
#(fcn, q, dt, quad_dw, params,nonlin!,jnonlin!))
function get_1st_imp_ext(fcnmf::Fcnmf, quad_in::QUAD, dt, quad_dw::QUAD, params, nonlin!, jnonlin!)
  # Inputs: fcn =  drift and diffusion functions
  # quad_in=current quadrature approx
  # dt=time step
  #
  # Outputs: quad_out the next timestep quad approx
  P=quad_in.length
  sqrt_dt=sqrt(dt)
  Y=zeros(P,1)
  Z=zeros(P,1)
  drift1=z->fcnmf.drift_ext(z,quad_in)
  diff1 =z->fcnmf.diff_ext(z, quad_in)
  helper=(y1,y2)->NLsolve.nlsolve(
    (u_,out)-> nonlin!(u_,y1, out,quad_in),                       (u_,out)->jnonlin!(u_,    out,quad_in),
    [y2])
  for j=1:P
    x=quad_in.X[j]
    xpfx_dt  = x + (1-params["Alpha"])*drift1(x)*dt
    gx_sqrth = diff1(x)*sqrt_dt
    #
    un=xpfx_dt+gx_sqrth
    solve=helper(un,x)
    Y[j]=solve.zero[1]
    #
    un=xpfx_dt-gx_sqrth
    solve=helper(un,x)
    Z[j]=solve.zero[1]
  end
  quad_out=QUAD(2*P, vec([Y Z]), [quad_in.w;quad_in.w]/2)
  return quad_out
end
# #####################################################
# function getYZA_implicit(mf::Fcnmf, quad_in::QUAD, dt)
#   # Inputs: fcn =  drift and diffusion functions
#   # quad_in=current quadrature approx
#   # dt=time step
#   #
#   # Outputs: quad_out the next timestep quad approx
#   P=quad_in.length
#   sqrt_dt=sqrt(dt)
#   Y=zeros(P,1)
#   Z=zeros(P,1)
#   A=zeros(P,1)
#   for j=1:P
#     x=quad_in.X[j]
#     xpfx_dt  = x + (1-mf.Alpha)*do_quad(y->mf.f([x,y]),quad_in)*dt
#     gx_sqrth = sqrt(3)*do_quad(y->mf.g([x,y],sqrt_dt),quad_in)
#     #
#     un=xpfx_dt+gx_sqrth
#     solve=NLsolve.nlsolve((u_,out)->mf.nonlin!(u_,un,out,dt,quad_in),
#                           (u_,out)->mf.jnonlin!(u_,out,dt,quad_in),
#                           [x])
#     #
#     Y[j]=solve.zero[1]
#     #
#     un=xpfx_dt-gx_sqrth
#     solve=NLsolve.nlsolve((u_,out)->mf.nonlin!(u_,un,out,dt,quad_in),
#                           (u_,out)->mf.jnonlin!(u_,out,dt,quad_in),
#                           [x])
#     #
#     Z[j]=solve.zero[1]
#     #
#     un=xpfx_dt
#     solve=NLsolve.nlsolve((u_,out)->mf.nonlin!(u_,un,out,dt,quad_in),
#                           (u_,out)->mf.jnonlin!(u_,out,dt,quad_in),
#                           [x])
#     #
#     A[j]=solve.zero[1]
#   end
#   quad_out=QUAD(3*P, vec([Y Z A]),
#               [quad_in.w/6;quad_in.w/6;quad_in.w*(2/3)])
#   return quad_out
# end
######################################################
function no_Gauss_points(m_old,dt,t_remaining,params)
  # function m=no_Gauss_points(m_old, dt, t_remaining, params)
  # return number of Gauss points (at least m_old),
  # to preserve dt^(order) accuracy
  #
  m_new=m_old
  if t_remaining>0
    ldt=abs(log(dt))
    order=params["order"]
    if params["periodic"]==true
      lldt_term=params["length_scale"]/params["no_partition"]
    else
      lldt_term=log((8*params["length_scale"]^2)*ldt)
    end
    lt_r=(0.0)*abs(log(t_remaining))
    if params["skiphalf"]==true
      order1=order+0.5
    else
      order1=order+1
    end
    constant=order1*ldt - 0.5*lldt_term -(order-1)*lt_r+params["const"]
    linear=lldt_term+lt_r
    g=y->f(y,constant,linear)
    m_large=1e4
    low  =g(m_old)
    upper=g(m_large)
    if (low*upper<0)
      m_new=Roots.fzero(g, [m_old, m_large])
    end
  end
  return round(Int16,m_new)
end
#####################################################
function f(m, constant, linear)
  # helper function for no_Gauss_points (roots gives a value for number of Gauss points)
  # Inputs: m=number of Gauss points
  #
  y= constant + m * linear - lgamma(2*m-1)
  return y
end
#####################################################
function coeffs_order2(x,mf,La,fcn::Fcnmf)
  # Coefficients for second-order method
  #
  a    = fcn.drift([x;mf])
  grada=fcn.Jdrift([x;mf])
  b     = fcn.diff([x;mf])
  gradb =fcn.Jdiff([x;mf])
  #
  Lan=[a;La]
  #
  term1=gradb[1]*b
  term2=grada[1]*b +dot(gradb,Lan) + 0.5*fcn.Hdiff([x;mf]) * b^2
  term3=dot(grada,Lan) + 0.5 * fcn.Hdrift([x;mf])*b^2
  #
  return [a  b  term1  term2  term3] # row vector
end
#
function dt_order2(dt,quad_dw::QUAD)
  # multiplier for second-order method
  # want to generate increment by coeffs_order2 * dtmat
  dtmat=zeros(5,quad_dw.length)
  dtmat[1,:]=dt
  dtmat[2,:]=quad_dw.X
  dtmat[3,:]=0.5* (quad_dw.X.^2 - dt)
  dtmat[4,:]=0.5* quad_dw.X * dt
  dtmat[5,:]=0.5* dt^2
  return dtmat
end
#
end # module definition


text:
##################################################################
#
# Simplified interfaces that only uses Julia vectors and functions.
#
##################################################################

"""
    Insert a check that the given function call returns 0,
    throw an error otherwise. Only apply directly to function calls.
"""
macro checkflag(ex, throw_error = false)
    @assert Base.Meta.isexpr(ex, :call)
    fname = ex.args[1]
    quote
        flag = $(esc(ex))
        if flag < 0
            if $(esc(throw_error))
                @error($(string(fname, " failed with error code = ")), flag)
            else
                @warn($(string(fname, " failed with error code = ")), flag)
            end
        end
        flag
    end
end

mutable struct UserFunctionAndData{F}
    func::F
    data::Any

    UserFunctionAndData{F}(func, data) where {F} = new{F}(func, data)
end

UserFunctionAndData(func) = func
UserFunctionAndData(func, data::Nothing) = func

function kinsolfun(y::N_Vector, fy::N_Vector, userfun::UserFunctionAndData)
    userfun[].func(convert(Vector, fy), convert(Vector, y), userfun[].data)
    return KIN_SUCCESS
end

function kinsolfun(y::N_Vector, fy::N_Vector, userfun)
    userfun(convert(Vector, fy), convert(Vector, y))
    return KIN_SUCCESS
end

function kinsol(
    f,
    y0::Vector{Float64};
    userdata::Any = nothing,
    linear_solver = :Dense,
    jac_upper = 0,
    jac_lower = 0,
)
    # f, Function to be optimized of the form f(y::Vector{Float64}, fy::Vector{Float64})
    #    where `y` is the input vector, and `fy` is the result of the function
    # y0, Vector of initial values
    # return: the solution vector
    mem_ptr = KINCreate()
    (mem_ptr == C_NULL) && error("Failed to allocate KINSOL solver object")
    kmem = Handle(mem_ptr)

    y = copy(y0)

    # use the user_data field to pass a function
    #   see: https://github.com/JuliaLang/julia/issues/2554
    userfun = UserFunctionAndData(f, userdata)
    function getcfun(userfun::T) where {T}
        @cfunction(kinsolfun, Cint, (N_Vector, N_Vector, Ref{T}))
    end
    flag = @checkflag KINInit(kmem, getcfun(userfun), NVector(y0)) true
    if linear_solver == :Dense
        A = Sundials.SUNDenseMatrix(length(y0), length(y0))
        LS = Sundials.SUNLinSol_Dense(y0, A)
    elseif linear_solver == :Band
        A = Sundials.SUNBandMatrix(length(y0), jac_upper, jac_lower)
        LS = Sundials.SUNLinSol_Band(y0, A)
    end
    flag = @checkflag Sundials.KINDlsSetLinearSolver(kmem, LS, A) true
    flag = @checkflag KINSetUserData(kmem, userfun) true
    ## Solve problem
    scale = ones(length(y0))
    strategy = KIN_NONE
    flag = @checkflag KINSol(kmem, y, strategy, scale, scale) true

    empty!(kmem)

    return y
end

function cvodefun(t::Float64, y::N_Vector, yp::N_Vector, userfun::UserFunctionAndData)
    userfun.func(t, convert(Vector, y), convert(Vector, yp), userfun.data)
    return CV_SUCCESS
end

function cvodefun(t::Float64, y::N_Vector, yp::N_Vector, userfun)
    userfun(t, convert(Vector, y), convert(Vector, yp))
    return CV_SUCCESS
end

"""
`cvode(f, y0::Vector{Float64}, t::Vector{Float64}, userdata::Any=nothing;
       integrator=:BDF, reltol::Float64=1e-3, abstol::Float64=1e-6, callback=(mem,t,y)->true)`

* `f`, Function of the form
  `f(t, y::Vector{Float64}, yp::Vector{Float64})`
  where `y` is the input state vector, and `yp` is the output vector
  of time derivatives for the states `y`
* `y0`, Vector of initial values
* `t`, Vector of time values at which to record integration results
* `integrator`, the chosen integration algorithm. Default is `:BDF`
  , other option is `:Adams`
* `reltol`, Relative Tolerance to be used (default=1e-3)
* `abstol`, Absolute Tolerance to be used (default=1e-6)
* `callback`, Callback function of the form
  `callback(mem, t_k::Float64, y_k::Vector{Float64})::Bool`
  where `mem` is the integrator and `t_k`/`y_k` the time/state at timestep `k`.
  A return value of `false` exits the integrator at timestep `k` and returns the all the timesteps preceding `k`.
  (default=(mem,t_k,y_k)->true)

return: a solution matrix with time steps in `t` along rows and
        state variable `y` along columns
"""
function cvode(
    f::Function,
    y0::Vector{Float64},
    t::AbstractVector,
    userdata::Any = nothing;
    kwargs...,
)
    y = zeros(length(t), length(y0))
    n = cvode!(f, y, y0, t, userdata; kwargs...)
    return y[1:n, :]
end

function cvode!(
    f::Function,
    y::Matrix{Float64},
    y0::Vector{Float64},
    t::AbstractVector,
    userdata::Any = nothing;
    integrator = :BDF,
    reltol::Float64 = 1e-3,
    abstol::Float64 = 1e-6,
    callback = (x, y, z) -> true,
)
    if integrator == :BDF
        mem_ptr = CVodeCreate(CV_BDF)
    elseif integrator == :Adams
        mem_ptr = CVodeCreate(CV_ADAMS)
    end

    (mem_ptr == C_NULL) && error("Failed to allocate CVODE solver object")
    mem = Handle(mem_ptr)

    c = 1

    userfun = UserFunctionAndData(f, userdata)
    y0nv = NVector(y0)

    function getcfun(userfun::T) where {T}
        @cfunction(cvodefun, Cint, (realtype, N_Vector, N_Vector, Ref{T}))
    end
    flag = @checkflag CVodeInit(mem, getcfun(userfun), t[1], convert(N_Vector, y0nv)) true

    flag = @checkflag CVodeSetUserData(mem, userfun) true
    flag = @checkflag CVodeSStolerances(mem, reltol, abstol) true
    A = Sundials.SUNDenseMatrix(length(y0), length(y0))
    LS = Sundials.SUNLinSol_Dense(y0nv, A)
    flag = Sundials.@checkflag Sundials.CVDlsSetLinearSolver(mem, LS, A) true

    y[1, :] = y0
    ynv = NVector(copy(y0))
    tout = [0.0]
    for k in 2:length(t)
        flag = @checkflag CVode(mem, t[k], ynv, tout, CV_NORMAL) true
        if !callback(mem, t[k], ynv)
            break
        end
        y[k, :] = convert(Vector, ynv)
        c = c + 1
    end

    empty!(mem)
    Sundials.SUNLinSolFree_Dense(LS)
    Sundials.SUNMatDestroy_Dense(A)

    return c
end

function idasolfun(
    t::Float64,
    y::N_Vector,
    yp::N_Vector,
    r::N_Vector,
    userfun::UserFunctionAndData,
)
    userfun.func(
        t,
        convert(Vector, y),
        convert(Vector, yp),
        convert(Vector, r),
        userfun.data,
    )
    return IDA_SUCCESS
end

function idasolfun(t::Float64, y::N_Vector, yp::N_Vector, r::N_Vector, userfun)
    userfun(t, convert(Vector, y), convert(Vector, yp), convert(Vector, r))
    return IDA_SUCCESS
end

"""
`idasol(f, y0::Vector{Float64}, yp0::Vector{Float64}, t::Vector{Float64}, userdata::Any=nothing;
        reltol::Float64=1e-3, abstol::Float64=1e-6, diffstates::Union{Vector{Bool},Void}=nothing)`

* `f`, Function of the form
  `f(t, y::Vector{Float64}, yp::Vector{Float64}, r::Vector{Float64})``
  where `y` and `yp` are the input state and derivative vectors,
  and `r` is the output residual vector
* `y0`, Vector of initial values
* `yp0`, Vector of initial values of the derivatives
* `reltol`, Relative Tolerance to be used (default=1e-3)
* `abstol`, Absolute Tolerance to be used (default=1e-6)
* `diffstates`, Boolean vector, true for the positions such that `r` depends on `yp[k]`

return: (y,yp) two solution matrices representing the states and state derivatives
         with time steps in `t` along rows and state variable `y` or `yp` along columns
"""
function idasol(
    f,
    y0::Vector{Float64},
    yp0::Vector{Float64},
    t::Vector{Float64},
    userdata::Any = nothing;
    reltol::Float64 = 1e-3,
    abstol::Float64 = 1e-6,
    diffstates::Union{Vector{Bool}, Nothing} = nothing,
)
    mem_ptr = IDACreate()
    (mem_ptr == C_NULL) && error("Failed to allocate IDA solver object")
    mem = Handle(mem_ptr)

    yres = zeros(length(t), length(y0))
    ypres = zeros(length(t), length(y0))

    userfun = UserFunctionAndData(f, userdata)

    function getcfun(userfun::T) where {T}
        @cfunction(idasolfun, Cint, (realtype, N_Vector, N_Vector, N_Vector, Ref{T}))
    end
    flag = @checkflag IDAInit(mem, getcfun(userfun), t[1], y0, yp0) true
    flag = @checkflag IDASetUserData(mem, userfun) true
    flag = @checkflag IDASStolerances(mem, reltol, abstol) true

    A = Sundials.SUNDenseMatrix(length(y0), length(y0))
    LS = Sundials.SUNLinSol_Dense(y0, A)
    flag = Sundials.@checkflag Sundials.IDADlsSetLinearSolver(mem, LS, A) true

    rtest = zeros(length(y0))
    f(t[1], y0, yp0, rtest)
    if any(abs.(rtest) .>= reltol)
        if diffstates === nothing
            error("Must supply diffstates argument to use IDA initial value solver.")
        end
        flag = @checkflag IDASetId(mem, collect(Float64, diffstates)) true
        flag = @checkflag IDACalcIC(mem, IDA_YA_YDP_INIT, t[2]) true
    end
    yres[1, :] = y0
    ypres[1, :] = yp0
    y = copy(y0)
    yp = copy(yp0)
    tout = [0.0]
    for k in 2:length(t)
        retval = @checkflag IDASolve(mem, t[k], tout, y, yp, IDA_NORMAL) true
        yres[k, :] = y
        ypres[k, :] = yp
    end

    empty!(mem)
    Sundials.SUNLinSolFree_Dense(LS)
    Sundials.SUNMatDestroy_Dense(A)

    return yres, ypres
end


text:
struct IAPWS95Consts <: EoSParam
    Tc::Float64
    Pc::Float64
    Vc::Float64
    Mw::Float64
    n::Vector{Float64}
    t::Vector{Float64}
    d::Vector{Int}
    c::Vector{Int}
    function IAPWS95Consts()
        n = [0.12533547935523e-1, 0.78957634722828e1, -0.87803203303561e1,
        0.31802509345418, -0.26145533859358, -0.78199751687981e-2,
        0.88089493102134e-2,

        -0.66856572307965, 0.20433810950965, -0.66212605039687e-4,
        -0.19232721156002, -0.25709043003438, 0.16074868486251,
        -0.040092828925807, 0.39343422603254e-6, -0.75941377088144e-5,
        0.56250979351888e-3, -0.15608652257135e-4, 0.11537996422951e-8,
        .36582165144204e-6, -.13251180074668e-11, -.62639586912454e-9,
        -0.10793600908932, 0.17611491008752e-1, 0.22132295167546,
        -0.40247669763528, 0.58083399985759, 0.49969146990806e-2,
        -0.31358700712549e-1, -0.74315929710341, 0.47807329915480,
        0.20527940895948e-1, -0.13636435110343, 0.14180634400617e-1,
        0.83326504880713e-2, -0.29052336009585e-1, 0.38615085574206e-1,
        -0.20393486513704e-1, -0.16554050063734e-2, .19955571979541e-2,
        0.15870308324157e-3, -0.16388568342530e-4, 0.43613615723811e-1,
        0.34994005463765e-1, -0.76788197844621e-1, 0.22446277332006e-1,
        -0.62689710414685e-4, -0.55711118565645e-9, -0.19905718354408,
        0.31777497330738, -0.11841182425981]

        d = [1, 1, 1, 2, 2, 3, 4,

        1, 1, 1, 2, 2, 3, 4, 4, 5, 7, 9, 10, 11, 13, 15, 1, 2, 2, 2, 3,
        4, 4, 4, 5, 6, 6, 7, 9, 9, 9, 9, 9, 10, 10, 12, 3, 4, 4, 5, 14, 3, 6, 6, 6]

        t = [-0.5, 0.875, 1, 0.5, 0.75, 0.375, 1,
        4, 6, 12, 1, 5, 4, 2, 13, 9, 3, 4, 11, 4, 13, 1, 7, 1, 9, 10, 10, 3, 7,
        10, 10, 6, 10, 10, 1, 2, 3, 4, 8, 6, 9, 8, 16, 22, 23,23, 10, 50, 44, 46, 50]

        c = [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 2, 2, 2, 2, 2, 2,
        2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 2, 3, 3, 3, 3, 4, 6, 6,6, 6]
        return new(647.096,2.2064e7,5.594803726708074e-5,18.015268,
        n,t,d,c)
    end
end

struct IAPWS95 <: EmpiricHelmholtzModel
    components::Array{String,1}
    consts::IAPWS95Consts
    references::Vector{String}
end

function crit_pure(model::IAPWS95)
    return (model.consts.Tc,model.consts.Pc,model.consts.Vc)
end

IAPWS95() = IAPWS95(["water"],IAPWS95Consts(),["IAPWS R6-95(2018)"])

function _f0(_model::IAPWS95,δ,τ)
    n= (-8.3204464837497, 6.6832105275932, 3.00632,0.012436, 0.97315, 1.2795, 0.96956, 0.24873)
    γ = (0.0, 0.0, 0.0, 1.28728967, 3.53734222, 7.74073708, 9.24437796,27.5075105)
    res = log(δ)+n[1] + n[2]*τ + n[3]*log(τ)
    @inbounds for i = 4:8
        res += n[i]*log(-expm1(-γ[i]*τ))
    end
    return res
end

function _fr(model::IAPWS95,δ,τ)
    n = model.consts.n::Vector{Float64}
    d = model.consts.d::Vector{Int64}
    t = model.consts.t::Vector{Float64}
    c = model.consts.c::Vector{Int64}
    res=zero(δ+τ)
    for i = 1:7
        res += n[i]* (δ^d[i]) * (τ^t[i])
    end

    for i = 8:51
        ic = i-7
        res += n[i]* (δ^d[i]) * (τ^t[i]) * exp(-δ^c[ic])
    end

    nτt1,nτt2,nτt3 = (-0.31306260323435e2, 0.31546140237781e2*τ, -0.25213154341695e4*τ^4)
    δd = δ^3
    #for i = 52:54
    #(δ-ε)^2
    _δ = (δ-1.0)^2#
    gauss_d = 20*_δ
    gauss_1 = exp(-150*((τ-1.21)^2)  -gauss_d)
    gauss_2 = gauss_1
    gauss_3 = exp(-250*((τ-1.25)^2)  -gauss_d)

    res += δd*(nτt1*gauss_1 + nτt2*gauss_2 + nτt3*gauss_3)

    #for i = 55:56

    nδ1,nδ2 = (-0.14874640856724*δ, 0.31806110878444*δ)
    _τ = (τ-1.0)^2
    Θ = (1.0-τ) + 0.32*_δ^(1.6666666666666667)
    Δ = Θ^2 + 0.2*_δ^3.5
    Ψ1 = exp(- 28*_δ - 700*_τ)
    Ψ2 = exp(- 32*_δ - 800*_τ)
    Δb1 = Δ^0.85
    Δb2 = Δ^0.95
    res += nδ1*Δb1*Ψ1 + nδ2*Δb2*Ψ2

    return res

end

function a_ideal(model::IAPWS95,V,T,z=SA[1.0])
    Σz = only(z) #single component
    v = V/Σz
     mass_v =  v*1000.0*0.055508472036052976
     rho = one(mass_v)/mass_v
     δ = rho*0.003105590062111801 #/322
     τ = 647.096/T
     return 0.9999890238768239*_f0(model,δ,τ)
end

function a_res(model::IAPWS95,V,T,z=SA[1.0])
    Σz = only(z) #single component
    v = V/Σz
    #R value calculated from molecular weight and specific gas constant
     #return 8.3143713575874*T*_f(model, molar_to_weight(1/v,[model.molecularWeight],[1.0]),T)
     #println(molar_to_weight(1/v,[model.molecularWeight],[1.0]))'
     mass_v =  v*1000.0*0.055508472036052976
     rho = one(mass_v)/mass_v
     δ = rho*0.003105590062111801 #/322
     τ = 647.096/T
     return 0.9999890238768239*_fr(model,δ,τ)
end

function eos(model::IAPWS95, V, T, z=SA[1.0])
    Σz = only(z) #single component
    v = V/Σz
    #R value calculated from molecular weight and specific gas constant
     #return 8.3143713575874*T*_f(model, molar_to_weight(1/v,[model.molecularWeight],[1.0]),T)
     #println(molar_to_weight(1/v,[model.molecularWeight],[1.0]))'
     mass_v =  v*1000.0*0.055508472036052976
     rho = one(mass_v)/mass_v
     δ = rho*0.003105590062111801 #/322
     τ = 647.096/T
    return R̄*Σz*T*0.9999890238768239*(_fr(model,δ,τ)+_f0(model,δ,τ))
end

function eos_res(model::IAPWS95, V, T, z=SA[1.0])
    Σz = only(z) #single component
    v = V/Σz
    #R value calculated from molecular weight and specific gas constant
     #return 8.3143713575874*T*_f(model, molar_to_weight(1/v,[model.molecularWeight],[1.0]),T)
     #println(molar_to_weight(1/v,[model.molecularWeight],[1.0]))'
     mass_v =  v*1000.0*0.055508472036052976
     rho = one(mass_v)/mass_v
     δ = rho*0.003105590062111801 #/322
     τ = 647.096/T
    return R̄*Σz*T*0.9999890238768239*_fr(model,δ,τ)
end

const watersat_data = (;n = [0.116_705_214_527_67E4,
    -0.724_213_167_032_06E6,
    -0.170_738_469_400_92E2,
    0.120_208_247_024_70E5,
    -0.323_255_503_223_33E7,
    0.149_151_086_135_30E2,
    -0.482_326_573_615_91E4,
    0.405_113_405_420_57E6,
    -0.238_555_575_678_49,
    0.650_175_348_447_98E3]
)

#psat and tsat are modified from SteamTables.jl
function water_p_sat(t)
    t > 647.096 && return zero(t)/zero(t)
    n = watersat_data.n
    Θ = t + n[9]/(t - n[10])
    A =      Θ^2 + n[1]*Θ + n[2]
    B = n[3]*Θ^2 + n[4]*Θ + n[5]
    C = n[6]*Θ^2 + n[7]*Θ + n[8]
    P = (2C / (-B + √(B^2 - 4*A*C)))^4
    return P*1000000
end

function water_t_sat(p)
    p > 2.2064e7 && return zero(p)/zero(p)
    n = watersat_data.n
    P = p/1000000
    β = P^0.25
    E =      β^2 + n[3]*β + n[6]
    F = n[1]*β^2 + n[4]*β + n[7]
    G = n[2]*β^2 + n[5]*β + n[8]
    D = 2G / (-F - √(F^2 - 4*E*G))
    T = (n[10]+D-√((n[10]+D)^2 - 4(n[9]+n[10]*D)))/2.0
    return T
end

#from MoistAir.jl, liquid
function saturated_water_liquid(Tk)
    ρm= ( -0.2403360201e4 +
        Tk*(-0.140758895e1 +
        Tk * (0.1068287657e0 +
        Tk*(-0.2914492351e-3 + Tk*(0.373497936e-6 - Tk*0.21203787e-9))))) /
        ( -0.3424442728e1 + 0.1619785e-1*Tk )
    #kg/m3 -> m3/kg *(kg/mol)
    mw = 18.015268*0.001
    v = mw/ρm
end

function saturated_water_vapor(Tk)
    Blin = 0.70e-8 - 0.147184e-8 * exp(1734.29/Tk)
    Clin = 0.104e-14 - 0.335297e-17*exp(3645.09/Tk)
    p = water_p_sat(Tk)
    return R̄*Tk/p * (1.0 + p*(Blin  + p*Clin))
end

mw(model::IAPWS95) = model.consts.Mw
molecular_weight(model::IAPWS95,z=SA[1.]) = model.consts.Mw*1e-3

function x0_volume(model::IAPWS95,p,T,z=[1.0];phase = :unknown)
    if phase == :unknown || is_liquid(phase)
        if model.consts.Pc > p
            return sat_v = saturated_water_liquid(T)
        else
            return volume_virial(model,p,T,z) #must look for better initial point here
        end
    elseif is_vapour(phase)
        x0val = 1.1*saturated_water_vapor(T)
    elseif is_supercritical(phase)
        x0val = model.consts.Vc

    end
    return x0val
end

function x0_sat_pure(model::IAPWS95,T)
    vl = saturated_water_liquid(T)
    vg = saturated_water_vapor(T)
    x0  = (vl,vg)
    return log10.(x0)
end

function T_scale(model::IAPWS95,z=SA[1.0])
    return model.consts.Tc
end

function p_scale(model::IAPWS95,z=SA[1.0])
    return model.consts.Pc
end

Base.length(::IAPWS95) = 1

function Base.show(io::IO,model::IAPWS95)
    return eosshow(io,model)
end

function Base.show(io::IO,mime::MIME"text/plain",model::IAPWS95)
    return eosshow(io,mime,model)
end

lb_volume(model::IAPWS95, z=SA[1.0]; phase = :unknown) = 1.4393788065379039e-5

export IAPWS95,IAPWS95Ideal


text:
f = (x, y) -> x^3 – y + x * y
(x, y) -> x^3 – y + x * y

function (x)
    x + 2
end

(x) -> x + 2
x -> x + 2
() -> x * 5

(x, y, z) -> 3x + 2y - z

max_stars_repo_path: Module 1/Chapter03/Old/anonymous.jl

max_issues_repo_path: Module 1/Chapter03/Old/anonymous.jl

max_forks_repo_path: Module 1/Chapter03/Old/anonymous.jl


text:
include("lib_for_robot.jl")
function task20(r)
    go_to_west_south_corner_and_return_path!(r; go_around_barriers = true)

    my_ans = 0
    border_now = false
    side = Ost
    while !isborder(r, Nord)

        while !isborder(r, side)
            my_ans, border_now = find_special!(r, my_ans, border_now, Nord)
            move!(r,side)
        end
        my_ans, border_now = find_special!(r, my_ans, border_now, Nord)

        side = inverse_side(side)
        move!(r, Nord)
    end
    return my_ans
end


text:
mutable struct DynamicGenerator{
    M <: Machine,
    S <: Shaft,
    A <: AVR,
    TG <: TurbineGov,
    P <: PSS,
} <: DynamicInjection
    static_injector::Union{Nothing, Generator}
    ω_ref::Float64
    machine::M
    shaft::S
    avr::A
    prime_mover::TG
    pss::P
    n_states::Int64
    states::Vector{Symbol}
    ext::Dict{String, Any}
    internal::InfrastructureSystemsInternal
end

function DynamicGenerator(
    static_injector::Generator,
    ω_ref::Float64,
    machine::M,
    shaft::S,
    avr::A,
    prime_mover::TG,
    pss::P,
    ext::Dict{String, Any} = Dict{String, Any}(),
) where {M <: Machine, S <: Shaft, A <: AVR, TG <: TurbineGov, P <: PSS}

    n_states = (
        machine.n_states +
        shaft.n_states +
        avr.n_states +
        prime_mover.n_states +
        pss.n_states
    )
    states = vcat(machine.states, shaft.states, avr.states, prime_mover.states, pss.states)

    return DynamicGenerator{M, S, A, TG, P}(
        static_injector,
        ω_ref,
        machine,
        shaft,
        avr,
        prime_mover,
        pss,
        n_states,
        states,
        ext,
        InfrastructureSystemsInternal(),
    )
end

function DynamicGenerator(;
    static_injector::Generator,
    ω_ref::Float64,
    machine::M,
    shaft::S,
    avr::A,
    prime_mover::TG,
    pss::P,
    ext::Dict{String, Any} = Dict{String, Any}(),
) where {M <: Machine, S <: Shaft, A <: AVR, TG <: TurbineGov, P <: PSS}
    DynamicGenerator(static_injector, ω_ref, machine, shaft, avr, prime_mover, pss, ext)
end

IS.get_name(device::DynamicGenerator) = get_name(device.static_injector)
get_bus(device::DynamicGenerator) = get_bus(device.static_injector)
get_Sbase(device::DynamicGenerator) = device.machine.s_rated
get_states(device::DynamicGenerator) = device.states
get_n_states(device::DynamicGenerator) = device.n_states
get_ω_ref(device::DynamicGenerator) = device.ω_ref
get_V_ref(device::DynamicGenerator) = get_voltage(get_bus(device))
get_P_ref(device::DynamicGenerator) = get_activepower(device.static_injector)
get_Q_ref(device::DynamicGenerator) = get_reactivepower(device.static_injector)
get_basepower(device::DynamicGenerator) = get_basepower(device.static_injector)
get_machine(device::DynamicGenerator) = device.machine
get_shaft(device::DynamicGenerator) = device.shaft
get_avr(device::DynamicGenerator) = device.avr
get_prime_mover(device::DynamicGenerator) = device.prime_mover
get_pss(device::DynamicGenerator) = device.pss
get_static_injector(device::DynamicGenerator) = device.static_injector
get_ext(device::DynamicGenerator) = device.ext


text:
mutable struct GPRealisations{KEY}
    groupKeys::Vector{KEY}
    mgp::Dict{KEY,GPE}
    nobs::Int
    logNoise::Float64
    mean::Mean
    kernel::Kernel
    mll::Float64            # Marginal log-likelihood
    dmll::Vector{Float64}   # Gradient marginal log-likelihood
    function GPRealisations{KEY}(
            groupKeys::Vector{KEY}, mgp::Dict{KEY,GPE}, nobs::Int, logNoise::Float64,
            # groupKeys::Vector{KEY}, mgp::Dict{KEY,GPE}, nobs::Int, logNoise::Scalar,
            mean::Mean, kernel::Kernel) where {KEY}
        length(groupKeys) == length(mgp) || throw("groupKeys and mgp should have same length")
        gpreals = new(groupKeys, mgp, nobs, logNoise, mean, kernel)
        propagate_params!(gpreals)
        update_mll!(gpreals)
        return gpreals
    end
end
function GPRealisations(gpList::Vector{GPE}, groupKeys::Vector{KEY}) where {KEY}
    total_nobs = sum([gp.nobs for gp in gpList])
    # get kernel and parameters from first GP in list
    first_gp = gpList[1]
    dim = first_gp.dim
    logNoise = first_gp.logNoise
    kern = first_gp.kernel
    mean = first_gp.mean
    ngroups = length(groupKeys)
    @assert ngroups == length(gpList)
    gpDict = Dict{KEY, GPE}()
    istart = 1
    for j in 1:ngroups
        gp = gpList[j]
        key = groupKeys[j]
        nobs = gp.nobs
        gp.nobs > 0 || throw("empty group")
        gpDict[key] = gp
        istart += nobs
    end
    gpreals = GPRealisations{KEY}(
                groupKeys, gpDict,
                total_nobs, logNoise, mean, kern)
    return gpreals
end
function GPRealisations{KEY}(groupKeys, mgp, nobs,
        logNoise::Scalar, mean, kernel) where {KEY}
    logNoise_float = convert(Float64, logNoise)
    return GPRealisations{KEY}(groupKeys, mgp, nobs,
            logNoise_float, mean, kernel)
end

Base.getindex(gpreals::GPRealisations{KEY}, key::KEY) where {KEY} = gpreals.mgp[key]
Base.keys(gpreals::GPRealisations) = gpreals.groupKeys
Base.values(gpreals::GPRealisations) = values(gpreals.mgp)

function get_params(gpreals::GPRealisations;
                    noise::Bool=true, domean::Bool=true, kern::Bool=true)
    params = Float64[]
    if noise; push!(params, gpreals.logNoise); end
    if domean
        for gp in values(gpreals.mgp)
            append!(params, get_params(gp.mean))
        end
    end
    if kern; append!(params, get_params(gpreals.kernel)); end
    return params
end
function propagate_params!(gpreals::GPRealisations; noise::Bool=true, kern::Bool=true)
    for gp in values(gpreals.mgp)
        # harmonize parameters
        if kern
            gp.kernel = gpreals.kernel
        end
        if noise
            gp.logNoise = gpreals.logNoise
        end
    end
end
function set_params!(gpreals::GPRealisations, hyp::Vector{Float64};
                     noise::Bool=true, domean::Bool=true, kern::Bool=true)
    # println("mean=$(mean)")
    istart=1
    if noise
        gpreals.logNoise = hyp[istart]
        istart += 1
    end
    if domean
        for key in gpreals.groupKeys
            gp = gpreals.mgp[key]
            set_params!(gp.mean, hyp[istart:istart+num_params(gp.mean)-1])
            istart += num_params(gp.mean)
        end
    end
    if kern
        set_params!(gpreals.kernel, hyp[istart:end])
    end
    propagate_params!(gpreals, noise=noise, kern=kern)
end

function update_mll!(gpreals::GPRealisations)
    mll = 0.0
    for gp in values(gpreals.mgp)
        update_mll!(gp)
        mll += gp.mll
    end
    gpreals.mll = mll
    return mll
end
function update_mll_and_dmll!(gpreals::GPRealisations,
                              ααinvcKIs::Dict{Int,Matrix};
                              noise::Bool=true, domean::Bool=true, kern::Bool=true)
    gpreals.mll = 0.0
    gpreals.dmll = zeros(get_params(gpreals; noise=noise, domean=domean, kern=kern))
    imean=2
    ikern=collect((length(gpreals.dmll)-num_params(gpreals.kernel)+1):length(gpreals.dmll))
    for gp in values(gpreals.mgp)
        update_mll_and_dmll!(gp, ααinvcKIs[gp.nobs];
            noise=noise,  domean=domean, kern=kern)
        gpreals.mll += gp.mll
        dmll_indices=Int[]
        if noise
            push!(dmll_indices, 1)
        end
        if domean
            append!(dmll_indices, imean:imean+num_params(gp.mean)-1)
            imean+=num_params(gp.mean)
        end
        if kern
            append!(dmll_indices, ikern)
        end
        gpreals.dmll[dmll_indices] .+= gp.dmll
    end
    return gpreals.dmll
end
function update_mll_and_dmll!(gpreals::GPRealisations; kwargs...)
    ααinvcKIs = Dict{Int,Matrix}()
    for gp in values(gpreals.mgp)
        if haskey(ααinvcKIs, gp.nobs)
            continue
        end
        ααinvcKIs[gp.nobs] = Array{Float64}(gp.nobs, gp.nobs)
    end
    return update_mll_and_dmll!(gpreals, ααinvcKIs)
end

function get_optim_target(gpreals::GPRealisations;
                          noise::Bool=true, domean::Bool=true, kern::Bool=true)
    ααinvcKIs = Dict{Int,Matrix}()
    for gp in values(gpreals.mgp)
        if haskey(ααinvcKIs, gp.nobs)
            continue
        end
        ααinvcKIs[gp.nobs] = Array{Float64}(gp.nobs, gp.nobs)
    end
    function mll(hyp::Vector{Float64})
        try
            set_params!(gpreals, hyp; noise=noise, domean=domean, kern=kern)
            update_mll!(gpreals)
            return -gpreals.mll
        catch err
             if !all(isfinite(hyp))
                println(err)
                return Inf
            elseif isa(err, ArgumentError)
                println(err)
                return Inf
            elseif isa(err, Base.LinAlg.PosDefException)
                println(err)
                return Inf
            else
                throw(err)
            end
        end
    end

    function mll_and_dmll!(grad::Vector{Float64}, hyp::Vector{Float64})
        set_params!(gpreals, hyp; noise=noise, domean=domean, kern=kern)
        update_mll_and_dmll!(gpreals, ααinvcKIs; noise=noise, domean=domean, kern=kern)
        grad[:] = -gpreals.dmll
        return -gpreals.mll
    end
    function dmll!(grad::Vector{Float64}, hyp::Vector{Float64})
        mll_and_dmll!(grad, hyp)
    end

    func = OnceDifferentiable(mll, dmll!, mll_and_dmll!,
        get_params(gpreals;noise=noise,domean=domean,kern=kern))
    return func
end
function optimize!(gpreals::GPRealisations;
                   noise::Bool=true, domean::Bool=true, kern::Bool=true)
    func = get_optim_target(gpreals, noise=noise, domean=domean, kern=kern)
    init = get_params(gpreals;  noise=noise, domean=domean, kern=kern)  # Initial hyperparameter values
    nparams = length(init)
    lower = -Inf*ones(nparams)
    upper = Inf*ones(nparams)
    ikernstart=1
    if noise
        ikernstart+=1
        lower[1]=-10.0
    end
    upper[ikernstart:ikernstart+num_params(gpreals.kernel)-1]=20.0
    best_x = copy(init)
    best_y = Inf
    count = 0
    function myfunc(x::Vector, grad::Vector)
        count += 1
        if length(grad) > 0
            y = func.fg!(grad, x)
            if y < best_y
                best_x[:] = x
            end
            return y
        else
            try
                y = func.f(x)
                if y < best_y
                    best_x[:] = x
                end
                return y
            catch
                return Inf
            end
        end
    end

    opt = Opt(:LD_LBFGS, nparams)
    lower_bounds!(opt, lower)
    upper_bounds!(opt, upper)
    min_objective!(opt, myfunc)
    xtol_rel!(opt,1e-4)
    ftol_rel!(opt, 1e-20)

    (minf,minx,ret) = NLopt.optimize(opt, init)
    myfunc(best_x, [])
    return best_x, count
end


text:
cd(@__DIR__)
samples = readstring("all_samples.jl")

_, setup_rest = split(samples, "#setup")
setup_rest = split(setup_rest, "#cell")
setup = first(setup_rest)
cells = setup_rest[2:end]

open("visual_regression_funcs.jl", "w") do io
    println(io, setup)
    for (i, body) in enumerate(cells)
        body = chomp(body)
        println(io, "function testfunc_$i()")
        for line in split(body, "\n")
            println(io, " "^4, line)
        end
        println(io, "end")
    end
end

open("sample_funcs.jl", "w") do io
    println(io, setup)
    for (i, body) in enumerate(cells)
        body = chomp(body)
        println(io, "function testfunc_$i()")
        for line in split(body, "\n")
            println(io, " "^4, line)
        end
        println(io, "end")
    end
end


text:
import ReachabilityAnalysis: split

abstract type AbstractSplitter end

# ==================
# standard splitter
# ==================

struct Splitter{S} <: AbstractSplitter
    split_fun::S
end

split(s::Splitter, X) = s.split_fun(X)
Base.haskey(s::Splitter, k::Int) = k == 1
Base.getindex(s::Splitter, k::Int) = k == 1 ? s : error("key $k not found")

function NoSplitter()
    split_fun = X0 -> [X0]
    return Splitter(split_fun)
end

function BoxSplitter(partition=nothing)
    if partition == nothing
        # default: one split per dimension
        split_fun = X -> split(box_approximation(X), 2 * ones(Int, dim(X)))
    else
        split_fun = X -> split(box_approximation(X), partition)
    end
    return Splitter(split_fun)
end

function ZonotopeSplitter(generators=nothing, splits=nothing)
    if generators == nothing && splits == nothing
        # default: one split per generator
        function split_fun(Z)
            p = ngens(Z)
            generators = 1:p
            splits = ones(Int, p)
            return split(Z, generators, splits)
        end
    else
        split_fun = Z -> split(Z, generators, splits)
    end
    return Splitter(split_fun)
end

# ==================================
# splitter for different iterations
# ==================================

struct IndexedSplitter <: AbstractSplitter
    index2splitter::Dict{Int, Splitter}
end

Base.haskey(s::IndexedSplitter, k::Int) = haskey(s.index2splitter, k)
Base.getindex(s::IndexedSplitter, k::Int) = getindex(s.index2splitter, k)

# ==============================
# splitter based on state space
# ==============================

struct SignSplitter  <: AbstractSplitter
end

function split(s::SignSplitter, X::Interval{N}) where {N}
    l = low(X, 1)
    h = high(X, 1)
    if l < zero(N) && h > zero(N)
        return [Interval(l, zero(N)), Interval(zero(N), h)]
    else
        return [X]
    end
end


text:
@import ChargedParticleDynamics.GuidingCenter4d.TokamakDeeplyTrapped as problem

# run configuration
const Δt    = 800.
const ntime = 1250000

# plot configuration
const nplot = 1
const splot = true
const XLIM  = ( 1.00,  1.10)
const YLIM  = (-0.10, +0.10)

# output configuration
const title     = "Guiding Center 4D Tokamak Deeply Trapped"
const fig_dir   = run_id
const html_file = run_id * ".html"


text:
using EAGO, JuMP, Cthulhu

####### Example of interest
m = Model(EAGO.Optimizer)
set_optimizer_attribute(m, "output_iterations", 1)
set_optimizer_attribute(m, "iteration_limit", 5)
set_optimizer_attribute(m, "verbosity", 0)

# Define bounded variables
xL = [10.0; 0.0; 0.0; 0.0; 0.0; 85.0; 90.0; 3.0; 1.2; 145.0]
xU = [2000.0; 16000.0; 120.0; 5000.0; 2000.0; 93.0; 95.0; 12.0; 4.0; 162.0]
@variable(m, xL[1] <= x[1] <= xU[1])
@NLobjective(m, Max, -2.0 * x[1])

JuMP.optimize!(m)

####### Define code to reduce need to transverse tree to get to functions of interest
b = backend(m).optimizer.model.optimizer
nlr = b._working_problem._nonlinear_constr[1]
function g()
    EAGO.relax!(b, b._working_problem._nonlinear_constr[1], 1, true)
end
function q(b, nlr)
    EAGO.relax!(b, nlr, 1, true)
end
k(m) = JuMP.optimize!(m)

f() = k(m)

####### Define descend
@descend f()

max_stars_repo_name: PSORLab/RSActivationFunctions

max_issues_repo_name: PSORLab/RSActivationFunctions

max_forks_repo_name: PSORLab/RSActivationFunctions


text:
@testset "If" begin
s = """
wrap(pos::Int, len::Int) = true ? 1 : (pos > len ? len : pos)
"""
msgs = lintstr(s)
@test length(msgs) == 3
@test msgs[1].code == :W643
@test occursin("false branch is unreachable", msgs[1].message)
@test msgs[2].code == :I340
@test occursin("unused local variable", msgs[2].message)
@test msgs[3].code == :I340
@test occursin("unused local variable", msgs[3].message)

s = """
wrap(pos::Int, len::Int) = false ? 1 : (pos > len ? len : pos)
"""
msgs = lintstr(s)
@test length(msgs) == 1
@test msgs[1].code == :W642
@test occursin("true branch is unreachable", msgs[1].message)
@test Lint.line(msgs[1]) == 1

s = """
f(x) = (x=1) ? 1 : 2 # clearly not what we want
"""
msgs = lintstr(s)
@test msgs[1].code == :I472
@test occursin("assignment in the if-predicate clause", msgs[1].message)

s = """
f(x) = ifelse(length(x), 1 , 2) # clearly not what we want
"""
msgs = lintstr(s)
@test msgs[1].code == :E431
@test occursin("use of length(, msgs[1].message) in a Boolean context, use isempty()")

s = """
f(x,y) = (0 <= x < y = 6) ? 1 : 2 # clearly not what we want
"""
msgs = lintstr(s)
@test msgs[1].code == :I472
@test occursin("assignment in the if-predicate clause", msgs[1].message)

s = """
function f()
    if true
        println("hello")
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :W644
@test occursin("redundant if-true statement", msgs[1].message)

s = """
function f()
    v::Array{Int,1} = [1,2,3]
    if length(v)
        println("hello")
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E431
@test occursin("use of length(, msgs[1].message) in a Boolean context, use isempty()")

s = """
function f(b::Bool, x::Int, y::Int)
    a = b ? x : y
    @lintpragma("Info type a")
    a
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test occursin("typeof(a, msgs[1].message) == Int")

s = """
function f(b::Bool, x::Int, y::Any)
    a = b ? x : y
    @lintpragma("Info type a")
    a
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I271
@test occursin("typeof(a, msgs[1].message) == Any")

s = """
function f()
    if :a && !:b
        1
    else
        2
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :E512
@test msgs[1].variable == ":a"
@test occursin("Lint doesn't understand in a boolean context", msgs[1].message)
@test msgs[2].code == :E512
@test msgs[2].variable == ":b"
@test occursin("Lint doesn't understand in a boolean context", msgs[2].message)

s = """
function f(a, b)
    if a == 1 # MISSING && or ||
        b == 2
        1
    else
        2
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I571
@test occursin("the 1st statement under the true-branch is a boolean expression", msgs[1].message)

s = """
function f(a, b)
    if a == 1 # MISSING && or ||
        !(b < 2)
        1
    else
        2
    end
end
"""
msgs = lintstr(s)
@test msgs[1].code == :I571
@test occursin("the 1st statement under the true-branch is a boolean expression", msgs[1].message)

s = """
function f(a, b)
    if a == 1 (b < 2) || error("b needs to be < 2")
        1
    else
        2
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
function f(a, b)
    if a == 1 (b < 2) || return 3
        1
    else
        2
    end
end
"""
msgs = lintstr(s)
@test isempty(msgs)

s = """
1==1 && true
true && 1==1
1==1 || true || 1==1
1==1 || 1==1 || true
true || 1==1 || 1==1
"""
msgs = lintstr(s)
@test isempty(msgs)

# issue #149: it's not always safe to include modules
s = """
if 1 == 0  # suppress if false warning
    using FakeModule149
end
"""
msgs = lintstr(s)
@test isempty(msgs)
end


text:
module GaitSymmetry

export limits, reflect

export SymmetryFunc, Sel86, Rob87, Vag92, Plo05, Zif08, Roc14, Que20, Alv20, Alv20b

"""
    (::Type{<:SymmetryFunc})(x, y)

A functor-style type to calculate the symmetry between `x` and `y`.

# Examples
```julia-repl
julia> Zif08(1., 1.1)
3.029234437673622

julia> Rob87(1., 1.1)
9.52380952380953
```
"""
abstract type SymmetryFunc end

"""
    reflect(x)

Return an array where the first half is the reverse of `x` and inverted, and the last half
is `x`, unmodified. The first element of `x` will only occur once.

# Example
```julia-repl
julia> reflect(1//1:3//1)
5-element Array{Rational{Int64},1}:
 1//3
 1//2
 1//1
 2//1
 3//1
```
"""
function reflect(x)
    return [reverse(inv.(x[2:end])); x]
end

function (f::Type{<:SymmetryFunc})(x::T, y::U, args...) where {T, U}
    f(promote(x, y)..., args...)
end

"""
    inv(::Type{<:SymmetryFunc})(s, baseline=1)

Return the inverse of the given symmetry function. The inverse symmetry function calculates
two numbers that produce a similar symmetry, `s`, when used as inputs for the given symmetry
index. `baseline` is an assumption of the value of one of the two numbers.

# Examples
```julia-repl
julia> s = Rob87(2., 2.2)
9.52380952380953

julia> (x2, y2) = inv(Rob87)(s, 2.)
(2.0, 2.1999999999999997)

julia> s ≈ Rob87(x2, y2)
true

julia> (x1, y1) = inv(Rob87)(s1, 1.)
(1.0, 1.0999999999999999)

julia> s ≈ Rob87(x1, y1)
true
```
"""
Base.inv(::Type{<:SymmetryFunc})

"""
    convert(T <: SymmetryFunc, U <: SymmetryFunc) -> anonymous function

Return a function to convert from symmetry metric `U` to `T`

# Examples
```jldoctest
julia> s = Zif08(1., 1.1)
3.029234437673622

julia> convert(Rob87, Zif08)(s)
9.52380952380947

julia> convert(Rob87, Zif08)(s) ≈ Rob87(1,1.1)
true
```
"""
function Base.convert(::Type{T}, ::Type{U}) where {T <: SymmetryFunc, U <: SymmetryFunc}
    return (x) -> T(inv(U)(x)...)
end

function limits(S, T::Type{<:Number}=Float64)
    small_x = S(floatmin(T), floatmax(T))
    small_y = S(floatmax(T), floatmin(T))

    small_x = round(small_x; digits=2)
    small_y = round(small_y; digits=2)

    return minmax(small_x, small_y)
end

"""
    Sel86

Calculate the symmetry of `x` and `y` as the ratio ``y/x``

[1] R. Seliktar and J. Mizrahi, “Some Gait Characteristics of Below-Knee Amputees and Their
Reflection on the Ground Reaction Forces,” Engineering in Medicine, vol. 15, no. 1, pp.
27–34, Jan. 1986.
"""
struct Sel86 <: SymmetryFunc; end
const Ratio = Sel86

function (::Type{Sel86})(x::T, y::T) where T <: Number
    return y/x
end

Base.inv(::Type{Sel86}) = _inv_SymmetryRatio

function _inv_SymmetryRatio(s::T, baseline::T=one(T)) where T <: Number
    return (baseline, s / baseline)
end

"""
    Rob87

Calculate the symmetry of `x` and `y` as defined in [1]:
``
(y - x)/(0.5*(x+y))*100%
``

[1] R. O. Robinson, W. Herzog, and B. M. Nigg, “Use of force platform variables to quantify
the effects of chiropractic manipulation on gait symmetry,” J Manipulative Physiol Ther,
vol. 10, no. 4, pp. 172–176, Aug. 1987.
"""
struct Rob87 <: SymmetryFunc; end

function (::Type{Rob87})(x::T, y::T) where T <: Number
    (y - x)/(x+y)*200
end

Base.inv(::Type{Rob87}) = _inv_Rob87

function _inv_Rob87(s::T, baseline::T=one(T)) where T <: Number
    y = -baseline*(s+200)/(s-200)
    return (baseline, y)
end

"""
    Vag92

``
(y-x)/max(x,y)*100
``

[1] G. Vagenas and B. Hoshizaki, “A Multivariable Analysis of Lower Extremity Kinematic
Asymmetry in Running,” Journal of Applied Biomechanics, vol. 8, no. 1, pp.
11–29, Feb. 1992.

See also: [`Rob87`](@ref)
"""
struct Vag92 <: SymmetryFunc; end

function (::Type{Vag92})(x::T, y::T) where T <: Number
    return (y-x)/max(x,y)*100
end

Base.inv(::Type{Vag92}) = _inv_Vag92

function _inv_Vag92(s::T, baseline::T=one(T)) where T <: Number
    y_xmax = baseline*(1 + s/100)
    y_ymax = baseline/(1 - s/100)

    if Vag92(baseline, y_xmax) ≈ s
        y = y_xmax
    elseif Vag92(baseline, y_ymax) ≈ s
        y = y_ymax
    else
        # Logically should never happen
        throw(DomainError(s, "no possible x/y pair found"))
    end

    return (baseline, y)
end

"""
    Plo05

Calculate the symmetry of `x` and `y`, multiplied by 100% by default [1], as defined in [2]:
``
ln(y/x)
``

[1] M. Plotnik, N. Giladi, Y. Balash, C. Peretz, and J. M. Hausdorff, “Is freezing of gait
in Parkinson’s disease related to asymmetric motor function?,” Ann Neurol., vol. 57, no. 5,
pp. 656–663, May 2005, doi:10.1002/ana.20452.
[2] M. Plotnik, N. Giladi, and J. M. Hausdorff, “A new measure for quantifying the
bilateral coordination of human gait: effects of aging and Parkinson’s disease,” Exp Brain
Res, vol. 181, no. 4, pp. 561–570, Aug. 2007, doi:10.1007/s00221-006-0676-3.
"""
struct Plo05 <: SymmetryFunc; end

function (::Type{Plo05})(x::T, y::T; scaled::Bool=true) where T <: Number
    # !(signbit(x) ⊻ signbit(y)) || throw(DomainError((x,y), "Sign of `x` and `y` must match"))

    ga = log(y/x)
    if scaled
        ga *= 100
    end
    return ga
end

Base.abs(::Type{Plo05}) = _abs_Plo05
_abs_Plo05(x,y) = abs(Plo05(x,y))

Base.inv(::Type{Plo05}) = _inv_Plo05

function _inv_Plo05(s::T, baseline::T=one(T); scaled::Bool=true) where T <: Number
    !signbit(s) || throw(DomainError(s, "Plo05 only produces positive values"))

    if scaled
        y = baseline*exp(s/100)
    else
        y = baseline*exp(s)
    end
    return (baseline, y)
end

"""
    Zif08

Calculate the symmetry of `x` and `y`, originally representing the left and right sides,
respectively, as defined in [1]:
``
((45° - arctan(x/y))/90°)×100%
``

If `(45° - arctan(x/y)) > 90°`, the following equation is used:
``
((45° - arctan(x/y) - 180°)/90°)×100%
``

[1] R. A. Zifchock, I. Davis, J. Higginson, and T. Royer, “The symmetry angle: A novel,
robust method of quantifying asymmetry,” Gait & Posture, vol. 27, no. 4, pp. 622–627, May
2008, doi:10.1016/j.gaitpost.2007.08.006.
"""
struct Zif08 <: SymmetryFunc; end

function (::Type{Zif08})(x::T, y::T) where T <: Number
    a = 45 - atand(x, y)
    a -= (a > 90) ? 180 : 0
    return 100*a/90
end

Base.inv(::Type{Zif08}) = _inv_Zif08

function _inv_Zif08(s::T, baseline::T=one(T)) where T <: Number
    a′ = 90*s/100
    a′ += (a′ + 180 > 90) ? 180 : 0
    a′ = tand(45 - a′)
    return (baseline, baseline / a′)
end

"""
    Roc14(x, y)

``
y - x
``

[1] L. Rochester, B. Galna, S. Lord, and D. Burn, “The nature of dual-task interference
during gait in incident Parkinson’s disease,” Neuroscience, vol. 265, pp. 83–94, 2014, doi:
10.1016/j.neuroscience.2014.01.041.
"""
struct Roc14 <: SymmetryFunc; end

function (::Type{Roc14})(x::T, y::T) where T
    return y - x
end

Base.abs(::Type{Roc14}) = _abs_Roc14
_abs_Roc14(x,y) = abs(Roc14(x,y))

Base.inv(::Type{Roc14}) = _inv_Roc14

function _inv_Roc14(s::T, baseline::T=one(T)) where T
    return (baseline, s + baseline)
end

"""
    Que20

``
(x - y)/(max(0,x,y) - min(0,x,y))
``

[1] R. Queen, L. Dickerson, S. Ranganathan, and D. Schmitt, “A novel method for measuring
asymmetry in kinematic and kinetic variables: The normalized symmetry index,” J Biomech,
vol. 99, p. 109531, 2020, doi:10.1016/j.jbiomech.2019.109531.
"""
struct Que20 <: SymmetryFunc; end

function (::Type{Que20})(x::T, y::T) where T
    return (y - x)/(max(0,x,y) - min(0,x,y))
end

Base.inv(::Type{Que20}) = _inv_Que20

function _inv_Que20(s::T, baseline::T=one(T)) where T
    if s > 0
        y = -baseline/(s-1)
    else
        y = s*baseline + baseline
    end

    return (baseline, y)
end

"""
    Alv20(x, y)

``
(y-x)/sqrt(2*(x^2 + y^2))
``

[1] S. A. Alves, R. M. Ehrig, P. C. Raffalt, A. Bender, G. N. Duda, and A. N. Agres,

“Quantifying Asymmetry in Gait: The Weighted Universal Symmetry Index to Evaluate 3D Ground
Reaction Forces,” Frontiers Bioeng Biotechnology, vol. 8, p. 579511, 2020,
doi:10.3389/fbioe.2020.579511.
"""
struct Alv20 <: SymmetryFunc; end

function (::Type{Alv20})(x::T, y::T, σ::Union{Nothing,T}=nothing) where T
    if σ === nothing
        return (y - x)/hypot(√2*x, √2*y)
    else
        return (y - x)/hypot(√2*x, √2*y)*(1-√2*σ/√(2σ^2 + x^2 + y^2))
    end
end

Base.inv(::Type{Alv20}) = _inv_Alv20

function _inv_Alv20(s::T, baseline::T=one(T)) where T
    y = -baseline*(2s*√(1-s^2)+1)/(2s^2 -1)
    return (baseline, y)
end

struct Alv20b{T} <: SymmetryFunc
    σ::T
end

(alv::Alv20b{T})(x::T, y::T) where T = Alv20(x, y, alv.σ)
(::Type{Alv20b})(x::T, y::T) where T = Alv20(x, y)

end # module


text:
"""
Remarkable tagging system.
"""
module Tags

import Base: isless
using Base.Iterators
using AutoHashEquals
using DataStructures
using Distances
using Unicode
using ..Common

"""
An individual tag object representing a category that objects can be classified
under.
"""
@auto_hash_equals struct Tag
    name::String
end

"""
    tagname(::Tag)

Return the canonical, user-friendly name of this tag.
"""
tagname(t::Tag) = t.name
isless(t::Tag, u::Tag) = isless(tagname(t), tagname(u))

struct TagMatrix
    popularity::DefaultDict{Tag,Int,Int}
    correlation::DefaultDict{Tuple{Tag,Tag},Int,Int}
    uritable::Dict{String,Tag}

    TagMatrix() = new(
        DefaultDict{Tag,Int}(0),
        DefaultDict{Tuple{Tag,Tag}, Int}(0),
        Dict{String,Tag}())
end

"""
    tagobject(m::TagMatrix, t::AbstractString)

Obtain the `Tag` object corresponding to `t` from matrix `m`.
"""
function tagobject(m::TagMatrix, t::AbstractString)
    # Normalize the display name of the tag
    t = Unicode.normalize(t, :NFKC)

    # Normalize the tag name further to a URI-suitable variant
    uri = urinormalize(t)
    if haskey(m.uritable, uri)
        candidate = m.uritable[uri]
        if tagname(candidate) == t
            candidate
        else
            error("Cannot retrieve tag $(repr(t)) because it conflicts with tag $(candidate)")
        end
    else
        m.uritable[uri] = Tag(t)
    end
end
tagobject(::TagMatrix, t::Tag) = t

"""
    tags(m::TagMatrix)

Return an iterable over all tags seen in the TagMatrix `m`.
"""
tags(m::TagMatrix) = keys(m.popularity)

"""
    popularity(m::TagMatrix, t::AbstractString)

Return the popularity of tag `t` in `m`.
"""
popularity(m::TagMatrix, t::Tag) = m.popularity[t]
popularity(m::TagMatrix, t) = popularity(m, tagobject(m, t))

"""
    popular(m::TagMatrix)

Return a `Vector` of tags, ordered from most popular to least.
"""
popular(m::TagMatrix) = sort(collect(tags(m)), by=t -> -popularity(m, t))

"""
    joint(m::TagMatrix, t1, t2)

Return the total weight of entries which have both `t1` and `t2` as tags.
"""
joint(m::TagMatrix, t1::Tag, t2::Tag) =
    m.correlation[tuple(sort([t1, t2])...)]
joint(m::TagMatrix, t1, t2) = joint(m, tagobject(m, t1), tagobject(m, t2))

"""
    distance(m::TagMatrix, dist::SemiMetric, t::Tag, u::Tag)

Compute the total distance between tags `t` and `u` using the given metric.
Allowed metrics are Jaccard() and CosineDist() (which reduces to Ochiai
coefficient).
"""
function distance(m::TagMatrix, ::Jaccard, t::Tag, u::Tag)
    num = joint(m, t, u)
    den = popularity(m, t) + popularity(m, u) - num
    return 1.0 - num / den
end

function distance(m::TagMatrix, ::CosineDist, t::Tag, u::Tag)
    num = joint(m, t, u)
    den = √popularity(m, t) * √popularity(m, u)
    return 1.0 - num / den
end

"""
    relatedto(m::TagMatrix, tag, num=8)

Return up to the top `num` tags related to `tag`.
"""
function relatedto(m::TagMatrix, t::Tag, num=8)
    top = [(tag, distance(m, Jaccard(), t, tag))
           for tag in tags(m) if tag != t]
    filter!(x -> x[2] < 1, top)
    sort!(top, by=x -> x[2])
    take(top, num)
end
relatedto(m::TagMatrix, t, num=8) = relatedto(m, tagobject(m, t), num)

"""
    issubtag(m::TagMatrix, a, b)

Return `true` if `a` is a subtag of `b`. A tag `a` is defined to be a subtag of
`b` if at least 75% of items tagged with `a`, plus 0.75, are also tagged with
`b`, or if `a == b`.
"""
issubtag(m::TagMatrix, a::Tag, b::Tag) = a == b ||
    joint(m, a, b) >= 0.75 * popularity(m, a) + 0.75
issubtag(m::TagMatrix, a, b) = issubtag(m, tagobject(m, a), tagobject(m, b))

"""
    subtags(m::TagMatrix, tag)

Return all subtags of this tag, in order of size.
"""
function subtags(m::TagMatrix, t::Tag)
    result = [tag for tag in tags(m) if tag != t && issubtag(m, tag, t)]
    sort!(result; by=x -> popularity(m, x), rev=true)
    result
end
subtags(m::TagMatrix, t) = subtags(m, tagobject(m, t))

"""
    populate!(m::TagMatrix, tags, value=1)

Add an entry to the tag matrix, whose tags are given by `tags`, and whose
weight is given by `value`.
"""
function populate!(m::TagMatrix, tags::Vector{Tag}, value=1)
    for a in tags
        m.popularity[a] += value
        for b in tags
            if b > a
                m.correlation[(a, b)] += value
            end
        end
    end
end
populate!(m::TagMatrix, tags, value=1) =
    populate!(m, [tagobject(m, t) for t in tags], value)

struct TagTree
    root::Tag
    children::Vector{TagTree}
end
root(tr::TagTree) = tr.root
children(tr::TagTree) = tr.children

const TagForest = Vector{TagTree}

function forest(m::TagMatrix, rset=collect(tags(m)))
    rset = sort(rset; by=x -> popularity(m, x), rev=true)
    result = TagForest()
    while !isempty(rset)
        top = popfirst!(rset)
        lower = filter(x -> issubtag(m, x, top), rset)
        rset = filter(x -> !issubtag(m, x, top), rset)
        push!(result, TagTree(top, forest(m, lower)))
    end
    result
end

function tag! end

export TagMatrix, joint, relatedto, populate!, popular, subtags, tags, forest,
       root, children, tagname, Tag, tagobject, tag!

end