text:
printstyled("  units\n", color=:light_green)

# Cases from randSeisChannel
printstyled("    check that RandSeis uses valid UCUM units\n", color=:light_green)
for u in SeisIO.RandSeis.irregular_units
  @test(vucum(u))
  @test units2ucum(u) == u
end
S = randSeisData()
isv = validate_units(S)
@test isv == trues(S.n)

C = S[1]
@test validate_units(C)
@test_logs (:warn, "Error thrown for unit string: meters per second") SeisIO.vucum("meters per second")
@test vucum("m/s^2") == false

# Now some real units that show up constantly in web requests
printstyled("    test conversions to UCUM units\n", color=:light_green)

@test units2ucum("m/s^1")     == "m/s"
@test units2ucum("m/s**1")    == "m/s"
@test units2ucum("m s^-1")    == "m/s"
@test units2ucum("m*s^-1")    == "m/s"
@test units2ucum("m.s^-1")    == "m/s"
@test units2ucum("m⋅s^-1")    == "m/s"
@test units2ucum("m s**-1")   == "m/s"
@test units2ucum("m*s**-1")   == "m/s"
@test units2ucum("m.s**-1")   == "m/s"
@test units2ucum("m⋅s**-1")   == "m/s"
@test units2ucum("m × s^-1")  == "m/s"
@test units2ucum("m ⋅ s^-1")  == "m/s"

@test units2ucum("nm/s")      == "nm/s"
@test units2ucum("nm/s^1")    == "nm/s"
@test units2ucum("nm/s**1")   == "nm/s"
@test units2ucum("nm s^-1")   == "nm/s"
@test units2ucum("nm*s^-1")   == "nm/s"
@test units2ucum("nm.s^-1")   == "nm/s"
@test units2ucum("nm.s**-1")  == "nm/s"
@test units2ucum("nm*s**-1")  == "nm/s"
@test units2ucum("nm s**-1")  == "nm/s"

# Acceleration
@test units2ucum("m/s^2")   == "m/s2"
@test units2ucum("m/s**2")  == "m/s2"
@test units2ucum("m s^-2")  == "m/s2"
@test units2ucum("m*s^-2")  == "m/s2"
@test units2ucum("m.s^-2")  == "m/s2"
@test units2ucum("m⋅s^-2")  == "m/s2"
@test units2ucum("m s**-2") == "m/s2"
@test units2ucum("m*s**-2") == "m/s2"
@test units2ucum("m.s**-2") == "m/s2"
@test units2ucum("m⋅s**-2") == "m/s2"

# Fictitious cases; we're unlikely to see anything like these
@test units2ucum("s^43/s**3") == "s43/s3"
@test units2ucum("s**10/m^4") == "s10/m4"
@test units2ucum("s^43 m**-3") == "s43/m3"
@test units2ucum("s^43 cm^-3") == "s43/cm3"
@test units2ucum("s^43 cm^-3 g") == "s43/cm3.g"
@test units2ucum("s^43 cm**+3 g") == "s43.cm3.g"
@test units2ucum("m s^-1 K") == "m/s.K"
@test units2ucum("m s^-1*K") == "m/s.K"
@test units2ucum("s^43 cm^-3 V^-4") == "s43/cm3.V4"
@test units2ucum("s^43.cm^-3*V^-4") == "s43/cm3.V4"


text:
## ExoplanetsSysSim/examples/run_abc_parallel.jl
## (c) 2015 Eric B. Ford

if nworkers() == 1
  addprocs(8)
end

import ExoplanetsSysSim
import ABC
@everywhere include("abc_setup.jl")

import SysSimABC; @everywhere using SysSimABC

@everywhere abc_plan = setup_abc(3)

@time output = SysSimABC.run_abc(abc_plan)

@time output2 = SysSimABC.run_abc(abc_plan,output)
#@time output3 = SysSimABC.run_abc(abc_plan,output2)

ext: jl


text:
using PySyntax2, PyCall
using Test

@testset "PySyntax2.jl" begin
    # Write your own tests here.
    np = pyimport("numpy")
    np.random.seed(0)

    @test sprint(print, np.random.rand(10)) == "[0.548814, 0.715189, 0.602763, 0.544883, 0.423655, 0.645894, 0.437587, 0.891773, 0.963663, 0.383442]"
end


text:
@testset "Generic Evaluation" begin
    @testset "Generic Evaluation - Vector" begin
        in = 5
        x = ones(Int, in)
        y = ones(Int, in)
        ŷ = ones(Int, in)
        ri = abs(rand(Int)) % in + 1
        ŷ[ri] = 10

        # desired result
        # ŷ - y => [0, 9, 0, 0, ...]

        dataset = [(x, y)]
        # dummy funciton & variables
        model(_x) = ŷ
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [0.0, 1.0, 1.0, 0.0],))

        # find maximum per dims
        _maxfunc(y::AbstractVector, ŷ::AbstractVector) = maximum(abs.(y .- ŷ))
        _maxfunc(y::AbstractMatrix, ŷ::AbstractMatrix) = maximum(abs.(y .- ŷ), dims=[1])

        # find minimum per dims
        _minfunc(y::AbstractVector, ŷ::AbstractVector) = minimum(abs.(y .- ŷ))
        _minfunc(y::AbstractMatrix, ŷ::AbstractMatrix) = minimum(abs.(y .- ŷ), dims=[1])

        # check metric function is right
        @test _maxfunc(y, ŷ) == 9
        @test _minfunc(y, ŷ) == 0

        # single metric
        res = evaluation(dataset, model, statvals, _maxfunc)
        @test Int(res) == 9
        res = evaluation(dataset, model, statvals, _minfunc)
        @test Int(res) == 0

        # multiple metric
        res = evaluations(dataset, model, statvals, [_maxfunc, _minfunc])
        @test Int.(res) == (9, 0)
    end

     @testset "Generic Evaluation - Matrix" begin
        in = 5
        batch = 7
        x = ones(Int, in, batch)
        y = ones(Int, in, batch)
        ŷ = ones(Int, in, batch)
        for i in 1:batch
            ri = abs(rand(Int)) % in + 1
            ŷ[ri, i] = 10
        end

        # desired result
        # ŷ - y => [0, 9, 0, 0, ...]

        dataset = [(x, y)]
        # dummy funciton & variables
        model(_x) = ŷ
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [0.0, 1.0, 1.0, 0.0],))

        # find maximum per dims
        _maxfunc(y::AbstractVector, ŷ::AbstractVector) = maximum(abs.(y .- ŷ))
        _maxfunc(y::AbstractMatrix, ŷ::AbstractMatrix) = maximum(abs.(y .- ŷ), dims=[1])

        # find minimum per dims
        _minfunc(y::AbstractVector, ŷ::AbstractVector) = minimum(abs.(y .- ŷ))
        _minfunc(y::AbstractMatrix, ŷ::AbstractMatrix) = minimum(abs.(y .- ŷ), dims=[1])

        # check metric function is right
        @test _maxfunc(y, ŷ) == reshape(ones(Int, batch) * 9, 1, batch)
        @test _minfunc(y, ŷ) == reshape(ones(Int, batch) * 0, 1, batch)

        # single metric
        res = evaluation(dataset, model, statvals, _maxfunc)
        @test Int(res) == 9
        res = evaluation(dataset, model, statvals, _minfunc)
        @test Int(res) == 0

        # multiple metric
        res = evaluations(dataset, model, statvals, [_maxfunc, _minfunc])
        @test Int.(res) == (9, 0)
    end
end

@testset "RSR" begin
    @testset "best case (should be zero)" begin
        n = 100
        a = rand(n)
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = a
        @test isapprox(evaluation([(a, a)], model, statvals, :RSR), 0.0)
    end

    @testset "worst case (should be one)" begin
        n = 100
        a = ones(n)
        b = zeros(n)
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = b
        @test isapprox(evaluation([(a, b)], model, statvals, :RSR), 0.0)
    end
end

@testset "NSE" begin
    @testset "best case (should be one)" begin
        n = 100
        a = rand(n)
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = a
        @test isapprox(evaluation([(a, a)], model, statvals, :NSE), 1.0)
    end
end

@testset "IOA" begin
    @testset "best case (should be one)" begin
        n = 100
        a = rand(n)
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = a
        @test isapprox(evaluation([(a, a)], model, statvals, :IOA), 1.0)
    end

    @testset "worst case (should be zero)" begin
        n = 100
        a = ones(n)
        b = zeros(n)
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = b
        @test isapprox(evaluation([(a, b)], model, statvals, :IOA), 0.0)
    end
end

@testset "R2" begin
    @testset "best case (should be one)" begin
        n = 100
        a = rand(n)
         statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = a
        @test isapprox(evaluation([(a, a)], model, statvals, :R2), 1.0)
    end

    @testset "worst case (should be zero)" begin
        n = 100
        a = ones(n)
        b = zeros(n)
        statvals = ndsparse((
            dataset = ["total", "total", "total", "total"],
            type = ["μ", "σ", "maximum", "minimum"]),
            (value = [mean(a), std(a), 1.0, 0.0],))
        model(_x) = b
        @test isapprox(evaluation([(a, b)], model, statvals, :R2), 0.0)
    end
end


text:
struct Grenander end

struct FittedGrenanderDistribution{T,VT<:AbstractVector{T}} <:
       ContinuousUnivariateDistribution
    fs::VT
    Fs::VT
    locs::VT
end

Distributions.@distr_support FittedGrenanderDistribution 0.0 d.locs[end]

function fit(::Grenander, ps)
    low_level_grenander_fit = MultipleTesting.grenander(ps)
    density_rle = rle(low_level_grenander_fit[2])
    loc_idx = cumsum(density_rle[2])
    locs = [0.0; low_level_grenander_fit[1][loc_idx]]
    fs = density_rle[1]
    Fs = [0.0; cumsum(fs .* diff(locs))]
    FittedGrenanderDistribution(fs, Fs, locs)
end

function Distributions.pdf(grenander::FittedGrenanderDistribution, t::Real)
    if t > grenander.locs[end]
        return 0.0
    elseif t <= grenander.locs[1] #left-continuity
        idx = 1
    else
        idx = searchsortedfirst(grenander.locs, t) - 1
    end
    grenander.fs[idx]
end

function Distributions.cdf(grenander::FittedGrenanderDistribution, t::Real)
    if t > grenander.locs[end]
        F = 1.0
    elseif t <= grenander.locs[1] #locs[1] is 0.0, TODO: allow Dirac peak at 0?
        F = 0.0
    else
        idx = searchsortedfirst(grenander.locs, t) - 1
        loc_L = grenander.locs[idx]
        F_L = grenander.Fs[idx]
        loc_R = grenander.locs[idx+1]
        F_R = grenander.Fs[idx+1]
        λ = (t - loc_L) / (loc_R - loc_L)
        F = (1 - λ) * F_L + λ * F_R
    end
    F
end

function invert_subgradient(gr::FittedGrenanderDistribution, λ)
    if λ > gr.fs[1]
        loc_lb = 0.0#return Interval(0.0, 0.0)
        loc_ub = 0.0
    elseif λ < gr.fs[end]
        loc_lb = 1.0
        loc_ub = 1.0
    else
        idx = searchsortedfirst(gr.fs, λ, lt = >)
        loc_lb = gr.locs[idx]
        loc_lb_p1 = gr.locs[idx+1]
        loc_ub = pdf(gr, loc_lb_p1) == λ ? loc_lb_p1 : loc_lb
    end
    Interval(loc_lb, loc_ub)
end


text:
@testset "Operations with types" begin

    using CPIDataBase
    S = 10
    v = zeros(S, S)

    ## Test capitalize
    @test :capitalize in names(CPIDataBase)

    # Test capitalization with default base index
    cap = capitalize(v)
    @test all(cap .== 100)

    # Test capitalization with different base index
    @test all(capitalize(v, 110) .== 110)

    base_idx = rand(100:110, S)
    cap = capitalize(v, base_idx)
    @test all(100 .<= cap[1, :] .<= 110)

    # Test capitalize_addbase
    cap = CPIDataBase.capitalize_addbase(v, base_idx)
    @test all(100 .<= cap[1, :] .<= 110)
    @test size(cap, 1) == size(v, 1)+1

    ## Test varinterm
    # TODO

    ## Test varinteran
    # TODO

end


text:
param_ref  = vec(readdlm("../kepler62ef_planets.txt",',',Float64) )
include("../test_ttv.jl")

function calc_ttvs_given_param(param::Vector)
  output = test_ttv(5,40,20,param,WriteOutput=false,num_evals=1)
  return vcat( output[1], output[2] )
end

function set_target_data_obs(param::Vector)
  global const ttv_ref = calc_ttvs_given_param(param)
  global const sigma_ttv = 10/(60*24)*ones(length(ttv_ref))
  nothing
end

set_target_data_obs(param_ref)

param_pert_scale = 0.01
param_pert = sign(param_ref).*exp( log(abs(param_ref)).+param_pert_scale*randn(length(param_ref)) )

#=
ttv_ref    = calc_ttvs_given_param(param_ref)
sigma_ttv  = 10/(60*24)*ones(length(ttv_ref))
param_sgn  = sign(param_ref)
param_pert = param_sgn.*exp( log(abs(param_ref)).+0.001*randn(length(param_ref)) )
=#

function calc_chisq(theta::Vector)
   ttv_pred = calc_ttvs_given_param(theta)
   @assert( size(ttv_pred) == size(ttv_ref) )
   chisq = zero(eltype(theta))
   for i in 1:length(ttv_pred)
     delta = (ttv_pred[i]-ttv_ref[i])/sigma_ttv[i]
     chisq += delta*delta
   end
   return chisq
end

log_likelihood(theta::Vector) = -0.5*calc_chisq(theta)

using ForwardDiff
my_cache = ForwardDiff.ForwardDiffCache() # make new cache to pass in to our function

# (grad, allresults) = ForwardDiff.gradient(log_likelihood,param_ref, AllResults, cache=my_cache);
# log_target(theta::Vector) = log_likelihood(theta)
# grad_log_target(theta::Vector) = ForwardDiff.gradient(log_likelihood,theta, cache=my_cache);

function test_forw_diff(param_ref::Vector{Float64}; test_hessian::Bool = false)

println("# Testing evaluating function")
@time ttvs = calc_ttvs_given_param(param_ref);
@time ttvs_list = typeof(ttvs)[calc_ttvs_given_param(param_pert) for i in 1:1000];
println("# Ignore this (ttvs):", hash(ttvs_list))
return

println("# Testing evaluating function")
@time chisq = calc_chisq(param_ref);
@time chisq_list = typeof(chisq)[calc_chisq(param_pert) for i in 1:1000];
println("# Ignore this (Chisq):", hash(chisq_list))

println("# Testing gradient")
@time grad = ForwardDiff.gradient(log_likelihood, param_ref, cache=my_cache);
@time grad_list = typeof(grad)[ForwardDiff.gradient(log_likelihood, param_pert, cache=my_cache) for i in 1:100];
println("# Ignore this (gradient): ", hash(grad_list))

@time (grad, allresults) = ForwardDiff.gradient(log_likelihood,param_ref, AllResults, cache=my_cache);
@time grad_results_list = typeof((grad,allresults))[ForwardDiff.gradient(log_likelihood,param_pert, AllResults, cache=my_cache) for i in 1:100];
println("# Ignore this (eval & gradient): ", hash(grad_results_list))

println("# Testing Jacobbian")
@time jac = ForwardDiff.jacobian(calc_ttvs_given_param,param_ref, cache=my_cache);
@time jac_list = typeof(jac)[ForwardDiff.jacobian(calc_ttvs_given_param,param_ref, cache=my_cache) for i in 1:100];
println("# Ignore this (jacobian): ", hash(jac_list))

if test_hessian
println("# Testing Hessian")
@time hes_func = ForwardDiff.hessian(log_likelihood);
@time hes = hes_func(param_ref);
@time hes_list = typeof(hes)[hes_func(param_pert) for i in 1:100];
println("# Ignore this (hessian): ", hash(hes_list))

println("# Testing everything at once")
@time hes, allresults = ForwardDiff.hessian(log_likelihood,param_ref,AllResults);
@time hes, allresults = ForwardDiff.hessian(log_likelihood,param_pert,AllResults);
println("# Ignore this (hessian, all at once): ", hash(hesc_list))
end
return true
end # func test_forw_diff

test_forw_diff(param_ref)
# test_forw_diff(param_ref, test_hessian=true);  # Doesn't work (yet)

max_stars_repo_name: nwtuchow/TTVFaster

max_issues_repo_name: nwtuchow/TTVFaster

max_forks_repo_name: nwtuchow/TTVFaster


text:
"Return type from `rand(ω, args...)`.  Used for type stability."
function randrtype end

randrtype(ω, ::Type{T}) where T = T
randrtype(ω, ::UnitRange{T}) where T = T
randrtype(ω, ::Array{T}) where T = T


text:
"""
    ALTROSolver

Augmented Lagrangian Trajectory Optimizer (ALTRO) is a solver developed by the Robotic Exploration Lab at Stanford University.
    The solver is special-cased to solve Markov Decision Processes by leveraging the internal problem structure.

ALTRO consists of two "phases":
1) AL-iLQR: iLQR is used with an Augmented Lagrangian framework to solve the problem quickly to rough constraint satisfaction
2) Projected Newton: A collocation-flavored active-set solver projects the solution from AL-iLQR onto the feasible subspace to achieve machine-precision constraint satisfaction.

# Constructor

    ALTROSolver(prob::Problem, opts::SolverOptions; [infeasible, R_inf, kwarg_opts...])

The `infeasible` keyword is a boolean flag that specifies whether the solver should
augment the controls to make it artificially fully actuated, allowing state initialization.
The `R_inf` is the weight on these augmented controls. Any solver options can be passed
as additional keyword arguments and will be set in the solver.

# Getters
* `get_model`
* `get_objective`
* `get_constraints`
* `Altro.get_ilqr`
* `TO.get_initial_state`

# Other methods
* `Base.size`: returns `(n,m,N)`
* `TO.is_constrained`
"""
struct ALTROSolver{T,S,P} <: ConstrainedSolver{T}
    opts::SolverOptions{T}
    stats::SolverStats{T}
    solver_al::ALSolver{T,S}
    solver_pn::P
end

function ALTROSolver(prob::Problem{T}, opts::SolverOptions=SolverOptions();
        infeasible::Bool=false,
        R_inf::Real=1.0,
        use_static=usestaticdefault(get_model(prob)[1]),
        kwarg_opts...
    ) where {Q,T}
    if infeasible
        # Convert to an infeasible problem
        prob = InfeasibleProblem(prob, prob.Z, R_inf)

        # Set infeasible constraint parameters
        # con_inf = get_constraints(prob).constraints[end]
        # con_inf::ConstraintVals{T,Control,<:InfeasibleConstraint}
        # con_inf.params.μ0 = opts.penalty_initial_infeasible
        # con_inf.params.ϕ = opts.penalty_scaling_infeasible
    end
    set_options!(opts; kwarg_opts...)
    stats = SolverStats{T}(parent=solvername(ALTROSolver))
    solver_al = ALSolver(prob, opts, stats, use_static=use_static)
    solver_pn = ProjectedNewtonSolver(prob, opts, stats)
    S = typeof(solver_al.ilqr)
    solver = ALTROSolver{T,S,typeof(solver_pn)}(opts, stats, solver_al, solver_pn)
    reset!(solver)
    # set_options!(solver; opts...)
    solver
end

# Getters
get_ilqr(solver::ALTROSolver) = get_ilqr(solver.solver_al)
for method in (:(RD.dims), :(RD.state_dim), :(RD.errstate_dim), :(RD.control_dim),
              :(TO.get_trajectory), :(TO.get_objective), :(TO.get_model),
              :(TO.get_initial_state), :getlogger)
    @eval $method(solver::ALTROSolver) = $method(get_ilqr(solver))
end
TO.get_constraints(solver::ALTROSolver) = get_constraints(solver.solver_al)
stats(solver::ALTROSolver) = solver.stats
options(solver::ALTROSolver) = solver.opts
usestatic(solver::ALTROSolver) = usestatic(get_ilqr(solver))

is_constrained(solver::ALTROSolver) = !isempty(get_constraints(solver.solver_al))
solvername(::Type{<:ALTROSolver}) = :ALTRO

# Methods
function max_violation(solver::ALTROSolver)
    return max(max_violation(solver.solver_al), max_violation(solver.solver_pn))
end

function reset!(solver::ALTROSolver)
    # reset_solver!(solver)
    opts = options(solver)::SolverOptions
    reset!(stats(solver), opts.iterations, solvername(solver))
    reset!(solver.solver_al)
    # reset!(solver.solver_pn)

    # Reset constraints
    conset = get_constraints(solver)
    reset!(conset)
end


text:
# Helper functions for NEPTypes (imported from NEPTypes.jl)
export Mder_NEP;
export Mder_Mlincomb_NEP;

# Type representing function handle NEP
struct Mder_Mlincomb_NEP <: NEP
    n::Int;
    Mder_fun::Function #
    maxder_Mder::Int;
    Mlincomb_fun::Function
    maxder_Mlincomb::Int;
end
"""
    Mder_NEP(n,Mder_fun; maxder=max)

Creates a `NEP` from its `compute_Mder` function defined by the
function handle `Mder_fun`. The `Mder_fun(λ,der)` takes two parameters a
scalar `λ::Number`, derivative number  `der`. The size `n::Int` must also
be specified. The function `Mder_fun(λ,der)` should return the n x n matrix
corresponding to the  `der`th derivatve. If only a
limited number of derivatives are available, `maxder` should be set, e.g.,
if not derivatives are implemented, set `maxder=0`. The function
`compute_Mlicomb` is automatically available by (delegation) to
`compute_Mlincomb_from_Mder`.

Note: This is a convenience function it is not recommended for high performance computations, since, e.g., it will not maintain type stability.

# Example
 The following example defines a linear eigenvalue problem
`A0+λA1` defined in an external function.
```julia-repl
julia> using LinearAlgebra; # For the I function
julia> function my_Mder(s,der)
    A0=ones(Float64,3,3)-I; A0[1,1]=-1;
    A1=ones(Float64,3,3)*3; A1[2,3]=0;
    if (der==0)
       return A0+A1*s;
    elseif (der==1)
       return A1;
    else
       return zero(A0);
    end
end
julia> nep=Mder_NEP(3,my_Mder);
julia> (λ,v)=augnewton(nep,v=ones(3));
julia> norm(compute_Mder(nep,λ)*v)
5.551115123125783e-17
```
"""
function Mder_NEP(n,Mder_fun; maxder=typemax(Int))
    maxder_Mlincomb=-1; # This means we delegate every Mlincomb_call
    Mlincomb_fun=s->s; # Dummy function
    return Mder_Mlincomb_NEP(n,
                             Mder_fun,maxder,
                             Mlincomb_fun,maxder_Mlincomb);
end

"""
    Mder_Mlincomb_NEP(n,Mder_fun, [maxder_Mder,] Mlincomb_fun, [maxder_Mlincomb])

Creates a `NEP` from its `compute_Mder` and `compute_Mlincomb`functions
 defined by the function handles `Mder_fun` and `Mlincomb_fun`. The `Mlincomb_fun(λ,X)` takes two parameters a
scalar `λ::Number` and a matrix `X`.  The size `n::Int` must also
be specified. The function `Mlincomb_fun(λ,X)` should return a vector corresponding of the linear combination of derivatives multiplied with the vectors in X. If only a
limited number of derivatives are implemented, `maxder_Mder` or
`maxder_Mlincomb` should be set, e.g., if not derivatives are implemented,
set `maxder=0`.

See also `Mder_NEP`.

Note: This is a convenience function it is not recommended for high performance computations, since, e.g., it will not maintain type stability.

# Example

The following example defines a linear eigenvalue problem
`A0+λA1` defined in an external function.
```julia-repl
julia> using LinearAlgebra; # For the I function
julia> function my_Mder(s,der)
    A0=ones(Float64,3,3)-I; A0[1,1]=-1;
    A1=ones(Float64,3,3)*3; A1[2,3]=0;
    if (der==0)
       return A0+A1*s;
    elseif (der==1)
       return A1;
    else
       return zero(A0);
    end
end
julia> function my_Mlincomb(s,X) # Compute linear comb of derivatives
    A0=ones(Float64,3,3)-I; A0[1,1]=-1;
    A1=ones(Float64,3,3)*3; A1[2,3]=0;
    if (size(X,2) <= 1)
       return A0*X[:,1]+s*A1*X[:,1]
    else # This means: size(X,2) => 2
       return A0*X[:,1]+A1*(s*X[:,1]+X[:,2]);
    end
end
julia> nep=Mder_Mlincomb_NEP(3,my_Mder,my_Mlincomb);
julia> (λ,v)=augnewton(nep,v=[1.0; 2.3; 0.0])
julia> norm(compute_Mder(nep,λ)*v)
6.798699777552591e-17
```
"""
function Mder_Mlincomb_NEP(n,
                           Mder_fun, maxder_Mder,
                           Mlincomb_fun, maxder_Mlincomb::Int=typemax(Int))
    # Wrap the two functions to Function unless they are already Function
    if !isa(Mder_fun,Function)
        new_Mder_fun = (λ,der) -> Mder_fun(λ,der);
    else
        new_Mder_fun = Mder_fun;
    end
    if !isa(Mlincomb_fun,Function)
        new_Mlincomb_fun = (λ,X) -> Mlincomb_fun(λ,X);
    else
        new_Mlincomb_fun = Mlincomb_fun;
    end
    return Mder_Mlincomb_NEP(n,new_Mder_fun,maxder_Mder,
                             new_Mlincomb_fun,maxder_Mlincomb);
end
Mder_Mlincomb_NEP(n,Mder_fun,Mlincomb_fun) = Mder_Mlincomb_NEP(n,Mder_fun,typemax(Int),Mlincomb_fun);

function compute_Mder(nep::Mder_Mlincomb_NEP,λ::Number,der::Integer=0)
    if (der>nep.maxder_Mder)
        error("Derivatives higher than ",nep.maxder_Mder," are not available.");
    end
    return nep.Mder_fun(λ,der);
end

function compute_Mlincomb(nep::Mder_Mlincomb_NEP,λ::Number,V::AbstractVecOrMat)
    if ((size(V,2)-1)<=nep.maxder_Mlincomb)
        F=nep.Mlincomb_fun(λ,V); # Call external routine
        if (F isa AbstractMatrix) # implicitly assume it is a matrix (with one column or one row) cf #231
            return reshape(F,size(nep,1));
        else  # Implicitly assume it is a vector
            return F;
        end
    else
        # Delegate to if we do not have implementation of Mlincomb
        return compute_Mlincomb_from_Mder(nep,λ,V,ones(eltype(V),size(V,2)))
    end
end

function size(nep::Mder_Mlincomb_NEP)
    return (nep.n,nep.n)
end
function size(nep::Mder_Mlincomb_NEP,dim)
    return nep.n
end


text:
# ------------------------------------------------------------------
# Licensed under the MIT License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    Transform

A transform that takes a table as input and produces a new table.
Any transform implementing the `Transform` trait should implement the
[`apply`](@ref) function. If the transform [`isrevertible`](@ref),
then it should also implement the [`revert`](@ref) function.

A functor interface is automatically generated from the functions
above, which means that any transform implementing the `Transform`
trait can be evaluated directly at any table implementing the
[Tables.jl](https://github.com/JuliaData/Tables.jl) interface.
"""
abstract type Transform end

"""
    assertions(transform)

Returns a list of assertion functions for the `transform`. An assertion
function is a function that takes a table as input and checks if the table
is valid for the `transform`.
"""
function assertions end

"""
    isrevertible(transform)

Tells whether or not the `transform` is revertible, i.e. supports a
[`revert`](@ref) function. Defaults to `false` for new types.
"""
function isrevertible end

"""
    newtable, cache = apply(transform, table)

Apply the `transform` on the `table`. Return the `newtable` and a
`cache` to revert the transform later.
"""
function apply end

"""
    table = revert(transform, newtable, cache)

Revert the `transform` on the `newtable` using the `cache` from the
corresponding [`apply`](@ref) call and return the original `table`.
Only defined when the `transform` [`isrevertible`](@ref).
"""
function revert end

"""
    Stateless

This trait is useful to signal that we can [`reapply`](@ref) a transform
"fitted" with training data to "test" data without relying on the `cache`.
"""
abstract type Stateless <: Transform end

"""
    reapply(transform, table, cache)

Reapply the `transform` to (a possibly different) `table` using a `cache`
that was created with a previous [`apply`](@ref) call.
"""
function reapply end

"""
    Colwise

A transform that is applied column-by-column. In this case, the new type
only needs to implement [`colapply`](@ref), [`colrevert`](@ref) and
[`colcache`](@ref). Efficient fallbacks are provided that execute these
functions in parallel for all columns with multiple threads.
"""
abstract type Colwise <: Transform end

"""
    y = colapply(transform, x, c)

Apply `transform` to column `x` with cache `c` and return new column `y`.
"""
function colapply end

"""
    x = colrevert(transform, y, c)

Revert `transform` starting from column `y` with cache `c` and return
original column `x`. Only defined when the `transform` [`isrevertible`](@ref).
"""
function colrevert end

"""
    c = colcache(transform, x)

Produce cache `c` necessary to [`colapply`](@ref) the `transform` on `x`.
If the `transform` [`isrevertible`](@ref) then the cache `c` can also be
used in [`colrevert`](@ref).
"""
function colcache end

# --------------------
# TRANSFORM FALLBACKS
# --------------------

assertions(transform::Transform) = assertions(typeof(transform))
assertions(::Type{<:Transform}) = []

isrevertible(transform::Transform) = isrevertible(typeof(transform))
isrevertible(::Type{<:Transform}) = false

(transform::Transform)(table) = apply(transform, table) |> first

# --------------------
# STATELESS FALLBACKS
# --------------------

reapply(transform::Stateless, table, cache) = apply(transform, table) |> first

# ------------------
# COLWISE FALLBACKS
# ------------------

function apply(transform::Colwise, table)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # retrieve column names and values
  names = Tables.columnnames(table)
  cols  = Tables.columns(table)

  # function to transform a single column
  function colfunc(n)
    x = Tables.getcolumn(cols, n)
    c = colcache(transform, x)
    y = colapply(transform, x, c)
    (n => y), c
  end

  # parallel map with multiple threads
  vals = tcollect(colfunc(n) for n in names)

  # new table with transformed columns
  𝒯 = (; first.(vals)...) |> Tables.materializer(table)

  # cache values for each column
  𝒞 = last.(vals)

  # return new table and cache
  𝒯, 𝒞
end

function revert(transform::Colwise, newtable, cache)
  # basic checks
  @assert isrevertible(transform) "transform is not revertible"

  # transformed columns
  names = Tables.columnnames(newtable)
  cols  = Tables.columns(newtable)

  # function to transform a single column
  function colfunc(i)
    n = names[i]
    c = cache[i]
    y = Tables.getcolumn(cols, n)
    x = colrevert(transform, y, c)
    n => x
  end

  # parallel map with multiple threads
  vals = tcollect(colfunc(i) for i in 1:length(names))

  # new table with transformed columns
  (; vals...) |> Tables.materializer(newtable)
end

function reapply(transform::Colwise, table, cache)
  # basic checks
  for assertion in assertions(transform)
    assertion(table)
  end

  # retrieve column names and values
  names = Tables.columnnames(table)
  cols  = Tables.columns(table)

  # check that cache is valid
  @assert length(names) == length(cache) "invalid cache for table"

  # function to transform a single column
  function colfunc(i)
    n = names[i]
    c = cache[i]
    x = Tables.getcolumn(cols, n)
    y = colapply(transform, x, c)
    n => y
  end

  # parallel map with multiple threads
  vals = tcollect(colfunc(i) for i in 1:length(names))

  # new table with transformed columns
  (; vals...) |> Tables.materializer(table)
end

# ----------------
# IMPLEMENTATIONS
# ----------------

include("transforms/select.jl")
include("transforms/filter.jl")
include("transforms/rename.jl")
include("transforms/replace.jl")
include("transforms/coalesce.jl")
include("transforms/coerce.jl")
include("transforms/identity.jl")
include("transforms/center.jl")
include("transforms/scale.jl")
include("transforms/zscore.jl")
include("transforms/quantile.jl")
include("transforms/functional.jl")
include("transforms/eigenanalysis.jl")
include("transforms/sequential.jl")
include("transforms/parallel.jl")


text:
const YesDeploy = ["TravisCI", "GitHubActions", "GitLabCI"]
const GitHubPagesStyle = ["TravisCI", "GitHubActions"]

"""
    Logo(; light=nothing, dark=nothing)

Logo information for documentation.

## Keyword Arguments
- `light::AbstractString`: Path to a logo file for the light (default) theme.
- `dark::AbstractString`: Path to a logo file for the dark theme.
"""
@option struct Logo
    light::Union{String, Nothing} = nothing
    dark::Union{String, Nothing} = nothing
end

"""
    Documenter(;
        depoly="",
        make_jl=".ion/templates/package/docs/make.jl",
        index_md=".ion/templates/package/docs/src/index.md",
        assets=String[],
        logo=Logo(),
        canonical_url=make_canonical(T),
        makedocs_kwargs=Dict{Symbol, Any}(),
    )

Sets up documentation generation via [Documenter.jl](https://github.com/JuliaDocs/Documenter.jl).
Documentation deployment depends on `T`, where `T` is some supported CI plugin,
or `Nothing` to only support local documentation builds.

## Supported Type Parameters
- `GitHubActions`: Deploys documentation to [GitHub Pages](https://pages.github.com)
  with the help of [`GitHubActions`](@ref).
- `TravisCI`: Deploys documentation to [GitHub Pages](https://pages.github.com)
  with the help of [`TravisCI`](@ref).
- `GitLabCI`: Deploys documentation to [GitLab Pages](https://pages.gitlab.com)
  with the help of [`GitLabCI`](@ref).
- `NoDeploy` (default): Does not set up documentation deployment.

## Keyword Arguments
- `make_jl::AbstractString`: Template file for `make.jl`.
- `index_md::AbstractString`: Template file for `index.md`.
- `assets::Vector{<:AbstractString}`: Extra assets for the generated site.
- `logo::Logo`: A [`Logo`](@ref) containing documentation logo information.
- `canonical_url::Union{Function, Nothing}`: A function to generate the site's canonical URL.
  The default value will compute GitHub Pages and GitLab Pages URLs
  for [`TravisCI`](@ref) and [`GitLabCI`](@ref), respectively.
  If set to `nothing`, no canonical URL is set.
- `makedocs_kwargs::Dict{Symbol}`: Extra keyword arguments to be inserted into `makedocs`.
- `devbranch::Union{AbstractString, Nothing}`: Branch that will trigger docs deployment.

!!! note
    If deploying documentation with Travis CI, don't forget to complete
    [the required configuration](https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#SSH-Deploy-Keys-1).
"""
@plugin struct Documenter <: Plugin
    depoly::String = "" # TravisCI, GitHubActions, GitLabCI
    assets::Vector{String} = String[]
    logo::Logo = Logo()
    makedocs_kwargs::Dict{Symbol} = Dict{Symbol, Any}()
    canonical_url::Bool = !isempty(depoly)
    make_jl::String = default_file("docs", "make.jl")
    index_md::String = default_file("docs", "src", "index.md")
    devbranch::Union{String, Nothing} = nothing
end

gitignore(::Documenter) = ["/docs/build/"]
priority(::Documenter, ::Function) = DEFAULT_PRIORITY - 1  # We need SrcDir to go first.

function badges(p::Documenter)
    isempty(p.depoly) && return Badge[]

    if p.depoly in GitHubPagesStyle
        return [
            Badge(
                "Stable",
                "https://img.shields.io/badge/docs-stable-blue.svg",
                "https://{{{USER}}}.github.io/{{{PKG}}}.jl/stable",
            ),
            Badge(
                "Dev",
                "https://img.shields.io/badge/docs-dev-blue.svg",
                "https://{{{USER}}}.github.io/{{{PKG}}}.jl/dev",
            ),
        ]
    end

    if p.depoly == "GitLabCI"
        return Badge(
            "Dev",
            "https://img.shields.io/badge/docs-dev-blue.svg",
            # TODO: Support custom domain here.
            "https://{{{USER}}}.gitlab.io/{{{PKG}}}.jl/dev",
        )
    end
end

function view(p::Documenter, t::PackagePlan, pkg::AbstractString)
    base = Dict(
        "ASSETS" => map(basename, p.assets),
        "AUTHORS" => join(t.authors, ", "),
        "CANONICAL" => canonical_url(p, t, pkg),
        "HAS_ASSETS" => !isempty(p.assets),
        "MAKEDOCS_KWARGS" => map(((k, v),) -> k => repr(v), collect(p.makedocs_kwargs)),
        "PKG" => pkg,
        "REPO" => "$(t.host)/$(t.user)/$pkg.jl",
        "USER" => t.user,
        "BRANCH" => p.devbranch === nothing ? default_branch(t) : p.devbranch,
    )

    if p.depoly in GitHubPagesStyle
        base["HAS_DEPLOY"] = true
    end

    return base
end

function canonical_url(p::Documenter, t::PackagePlan, pkg::String)
    p.canonical_url || return
    if p.depoly in GitHubPagesStyle
        return github_pages_url(t, pkg)
    elseif p.depoly == "GitLabCI"
        return gitlab_pages_url(t, pkg)
    else
        return
    end
end

function validate(p::Documenter, t::PackagePlan)
    if isempty(p.depoly)
        foreach(p.assets) do a
            isfile(a) || throw(ArgumentError("Asset file $a does not exist"))
        end
        foreach((:light, :dark)) do k
            logo = getfield(p.logo, k)
            if logo !== nothing && !isfile(logo)
                throw(ArgumentError("Logo file $logo does not exist"))
            end
        end
    end

    invoke(validate, Tuple{Documenter, PackagePlan}, p, t)
    if !hasplugin(t, T)
        name = nameof(T)
        s = "Documenter: The $name plugin must be included for docs deployment to be set up"
        throw(ArgumentError(s))
    end
end

function hook(p::Documenter, t::PackagePlan, pkg_dir::AbstractString)
    pkg = basename(pkg_dir)
    docs_dir = joinpath(pkg_dir, "docs")

    # Generate files.
    make = render_file(p.make_jl, combined_view(p, t, pkg), tags(p))
    index = render_file(p.index_md, combined_view(p, t, pkg), tags(p))
    gen_file(joinpath(docs_dir, "make.jl"), make)
    gen_file(joinpath(docs_dir, "src", "index.md"), index)

    # Copy over any assets.
    assets_dir = joinpath(docs_dir, "src", "assets")
    mkpath(assets_dir)
    foreach(a -> cp(a, joinpath(assets_dir, basename(a))), p.assets)
    foreach((:light => "logo", :dark => "logo-dark")) do (k, f)
        logo = getfield(p.logo, k)
        if logo !== nothing
            _, ext = splitext(logo)
            cp(logo, joinpath(assets_dir, "$f$ext"))
        end
    end
    isempty(readdir(assets_dir)) && rm(assets_dir)

    # Create the documentation project.
    d = Dict{String, Any}(
        "deps" => Dict{String, Any}(
            "Documenter" => "e30172f5-a6a5-5a46-863b-614d45cd2de4"
        )
    )

    cd(docs_dir) do
        open("Project.toml", "w") do io
            TOML.print(io, d)
        end

        PkgCmd.withproject("Pkg.develop(PackageSpec(; path=\"..\"))", false, "dev $(t.name)")
    end
end

github_pages_url(t::PackagePlan, pkg::AbstractString) = "https://$(t.user).github.io/$pkg.jl"
gitlab_pages_url(t::PackagePlan, pkg::AbstractString) = "https://$(t.user).gitlab.io/$pkg.jl"

needs_username(::Documenter) = true

function customizable(::Type{<:Documenter})
    return (:canonical_url => NotCustomizable, :makedocs_kwargs => NotCustomizable)
end

function interactive(::Type{Documenter})
    styles = [NoDeploy, TravisCI, GitLabCI, GitHubActions]
    menu = RadioMenu(map(string, styles); pagesize=length(styles))
    println("Documenter deploy style:")
    idx = request(menu)
    return interactive(Documenter{styles[idx]})
end

function prompt(::Type{<:Documenter}, ::Type{Logo}, ::Val{:logo})
    light = Base.prompt("Enter value for 'logo.light' (String, default=nothing)")
    dark = Base.prompt("Enter value for 'logo.dark' (String, default=nothing)")
    return Logo(; light=light, dark=dark)
end


text:
using DataFrames, LinearAlgebra, MixedModels, Random, RData, SparseArrays, Statistics, StatsBase, Test

if !@isdefined(dat) || !isa(dat, Dict{Symbol, DataFrame})
    dat = Dict(Symbol(k) => v for (k, v) in load(joinpath(dirname(@__FILE__), "dat.rda")))
end

@testset "Dyestuff" begin
    fm1 = LinearMixedModel(@formula(Y ~ 1 + (1|G)), dat[:Dyestuff])

    @test nblocks(fm1.A) == (3, 3)
    @test size(fm1.trms) == (3, )
    @test nblocks(fm1.L.data) == (3, 3)
    @test lowerbd(fm1) == zeros(1)
    @test fm1.lowerbd == zeros(1)
    @test getθ(fm1) == ones(1)
    @test getΛ(fm1) == [1.0]
    fm1.θ = ones(1)
    @test getθ(fm1) == ones(1)

    @test objective(updateL!(setθ!(fm1, [0.713]))) ≈ 327.34216280955366
    MixedModels.describeblocks(IOBuffer(), fm1)

    fit!(fm1);
    @test :θ in propertynames(fm1)
    @test isapprox(objective(fm1), 327.3270598811428, atol=0.001)
    @test isapprox(getθ(fm1), [0.752580], atol=1.e-5)
    @test isapprox(fm1.λ, [0.752580], atol=1.e-5)
    @test isapprox(fm1.θ, [0.752580], atol=1.e-5)
    @test isapprox(deviance(fm1), 327.32705988, atol=0.001)
    @test isapprox(aic(fm1), 333.3270598811394, atol=0.001)
    @test isapprox(bic(fm1), 337.5306520261259, atol=0.001)
    @test fixef(fm1) ≈ [1527.5]
    @test fm1.β ≈ [1527.5]
    @test StatsBase.dof(fm1) == 3
    @test StatsBase.nobs(fm1) == 30
    @test MixedModels.fixef!(zeros(1),fm1) ≈ [1527.5]
    @test coef(fm1) ≈ [1527.5]
    @test fm1.σ == sdest(fm1)
    @test isapprox(fm1.σ, 49.510099986291145, atol=1.e-5)
    @test fm1.X == ones(30,1)
    @test fm1.y == dat[:Dyestuff][:Y]
    @test cond(fm1) == ones(1)
    cm = coeftable(fm1)
    @test length(cm.rownms) == 1
    @test length(cm.colnms) == 4
    @test MixedModels.fnames(fm1) == [:G]
    @test model_response(fm1) == dat[:Dyestuff][:Y]
    rfu = ranef(fm1, uscale = true)
    rfb = ranef(fm1)
    cv = condVar(fm1)
    @test abs(sum(rfu[1])) < 1.e-5
    @test length(cv) == 1
    @test size(cv[1]) == (1, 1, 6)
    show(IOBuffer(), fm1.optsum)

    @test isapprox(logdet(fm1), 8.06014522999825, atol=0.001)
    @test isapprox(varest(fm1), 2451.2501089607676, atol=0.001)
    @test isapprox(pwrss(fm1), 73537.49947885796, atol=0.001)
    @test isapprox(StatsBase.stderror(fm1), [17.69455188898009], atol=0.0001)

    vc = VarCorr(fm1)
    show(IOBuffer(), vc)
    @test vc.s == sdest(fm1)

    fit!(fm1, REML=true)
    @test isapprox(objective(fm1), 319.65427684225216, atol=0.0001)
    @test dof_residual(fm1) ≥ 0
    print(IOBuffer(), fm1)
end

@testset "Dyestuff2" begin
    fm = fit(LinearMixedModel, @formula(Y ~ 1 + (1 | G)), dat[:Dyestuff2])
    @test lowerbd(fm) == zeros(1)
    show(IOBuffer(), fm)
    @test getθ(fm)[1] < 1.0e-9
    @test objective(fm) ≈ 162.87303665382575
    @test abs(std(fm)[1][1]) < 1.0e-9
    @test std(fm)[2] ≈ [3.653231351374652]
    @test StatsBase.stderror(fm) ≈ [0.6669857396443261]
    @test coef(fm) ≈ [5.6656]
    @test logdet(fm) ≈ 0.0
    refit!(fm, dat[:Dyestuff][:Y])
    @test isapprox(objective(fm), 327.3270598811428, atol=0.001)
end

@testset "penicillin" begin
    fm = LinearMixedModel(@formula(Y ~ 1 + (1 | G) + (1 | H)), dat[:Penicillin]);
    @test size(fm) == (144, 1, 30, 2)
    @test getθ(fm) == ones(2)
    @test lowerbd(fm) == zeros(2)

    fit!(fm)

    @test isapprox(objective(fm), 332.18834867227616, atol=0.001)
    @test isapprox(coef(fm), [22.97222222222222], atol=0.001)
    @test isapprox(fixef(fm), [22.97222222222222], atol=0.001)
    @test coef(fm)[1] ≈ mean(dat[:Penicillin][:Y])
    @test isapprox(StatsBase.stderror(fm), [0.7445960346851368], atol=0.0001)
    @test isapprox(getθ(fm), [1.5375772376554968, 3.219751321180035], atol=0.001)
    @test isapprox(std(fm)[1], [0.8455645948223015], atol=0.0001)
    @test isapprox(std(fm)[2], [1.770647779277388], atol=0.0001)
    @test isapprox(varest(fm), 0.3024263987592062, atol=0.0001)
    @test isapprox(logdet(fm), 95.74614821367786, atol=0.001)
    rfu = ranef(fm, uscale=true)
    rfb = ranef(fm)
    @test length(rfb) == 2
end

@testset "pastes" begin
    fm = LinearMixedModel(@formula(Y ~ (1 | G) + (1 | H)), dat[:Pastes])
    @test size(fm) == (60, 1, 40, 2)
    @test getθ(fm) == ones(2)
    @test lowerbd(fm) == zeros(2)

    fit!(fm);

    @test isapprox(objective(fm), 247.99446586289676, atol=0.001)
    @test isapprox(coef(fm), [60.05333333333329], atol=0.001)
    @test isapprox(fixef(fm), [60.05333333333329], atol=0.001)
    @test isapprox(StatsBase.stderror(fm), [0.6421359883527029] , atol=0.0001)
    @test isapprox(getθ(fm), [3.5268858714382905, 1.3299230213750168], atol=0.001)
    @test isapprox(std(fm)[1], [2.904069002535747], atol=0.001)
    @test isapprox(std(fm)[2], [1.095070371687089], atol=0.0001)
    @test isapprox(std(fm)[3], [0.8234088395243269], atol=0.0001)
    @test isapprox(varest(fm), 0.6780020742644107, atol=0.0001)
    @test isapprox(logdet(fm), 101.0381339953986, atol=0.001)
end

@testset "InstEval" begin
    fm1 = LinearMixedModel(@formula(Y ~ 1 + A + (1 | G) + (1 | H) + (1 | I)), dat[:InstEval])
    @test size(fm1) == (73421, 2, 4114, 3)
    @test getθ(fm1) == ones(3)
    @test lowerbd(fm1) == zeros(3)

    fit!(fm1);

    @test isapprox(objective(fm1), 237721.7687745563, atol=0.001)
    ftd1 = fitted(fm1);
    @test size(ftd1) == (73421, )
    @test ftd1 == predict(fm1)
    @test isapprox(ftd1[1], 3.17876, atol=0.0001)
    resid1 = residuals(fm1);
    @test size(resid1) == (73421, )
    @test isapprox(resid1[1], 1.82124, atol=0.00001)

    fm2 = fit!(LinearMixedModel(@formula(Y ~ 1 + A*I + (1 | G) + (1 | H)), dat[:InstEval]))
    @test isapprox(objective(fm2), 237585.5534151694, atol=0.001)
    @test size(fm2) == (73421, 28, 4100, 2)
end

@testset "sleep" begin
    fm = LinearMixedModel(@formula(Y ~ 1 + U + (1 + U | G)), dat[:sleepstudy]);
    @test lowerbd(fm) == [0.0, -Inf, 0.0]
    A11 = fm.A[Block(1,1)]
    @test isa(A11, UniformBlockDiagonal{Float64})
    @test isa(fm.L.data[Block(1, 1)], UniformBlockDiagonal{Float64})
    @test size(A11) == (36, 36)
    @test A11.facevec[1] == [10. 45.; 45. 285.]
    @test length(A11.facevec) == 18
    updateL!(fm);
    b11 = LowerTriangular(fm.L.data[Block(1, 1)].facevec[1])
    @test b11 * b11' == fm.A[Block(1, 1)].facevec[1] + I
    @test count(!iszero, Matrix(fm.L.data[Block(1, 1)])) == 18 * 4

    fit!(fm)

    @test objective(fm) ≈ 1751.9393444647046
    @test isapprox(getθ(fm), [0.929221307, 0.01816838, 0.22264487096], atol=1.e-6)
    @test isapprox(pwrss(fm), 117889.46144025437)
    @test isapprox(logdet(fm), 73.90322021999222, atol=0.001)
    @test isapprox(StatsBase.stderror(fm), [6.632257721914501, 1.5022354739749826], atol=0.0001)
    @test coef(fm) ≈ [251.40510484848477,10.4672859595959]
    @test fixef(fm) ≈ [10.4672859595959, 251.40510484848477]
    @test isapprox(StatsBase.stderror(fm), [6.632246393963571, 1.502190605041084], atol=0.01)
    @test isapprox(std(fm)[1], [23.780468100188497, 5.716827903196682], atol=0.01)
    @test isapprox(logdet(fm), 73.90337187545992, atol=0.001)
    @test diag(cor(fm)[1]) ≈ ones(2)
    @test isapprox(cond(fm), [4.175251], atol=0.0001)
    @test loglikelihood(fm) ≈ -875.9696722323523
    show(IOBuffer(), fm)

    u3 = ranef(fm, uscale=true)
    @test length(u3) == 1
    @test size(u3[1]) == (2, 18)
    @test isapprox(u3[1][1, 1], 3.030300122575336, atol=0.001)
    u3n = ranef(fm, uscale=true, named=true)

    b3 = ranef(fm)
    @test length(b3) == 1
    @test size(b3[1]) == (2, 18)
    @test isapprox(b3[1][1, 1], 2.815819441982976, atol=0.001)

    simulate!(fm)  # to test one of the unscaledre methods

    fmnc = LinearMixedModel(@formula(Y ~ 1 + U + (1|G) + (0+U|G)), dat[:sleepstudy])
    @test size(fmnc) == (180,2,36,1)
    @test getθ(fmnc) == ones(2)
    @test lowerbd(fmnc) == zeros(2)

    fit!(fmnc)

    @test isapprox(deviance(fmnc), 1752.0032551398835, atol=0.001)
    @test isapprox(objective(fmnc), 1752.0032551398835, atol=0.001)
    @test coef(fmnc) ≈ [251.40510484848585, 10.467285959595715]
    @test fixef(fmnc) ≈ [10.467285959595715, 251.40510484848477]
    @test isapprox(StatsBase.stderror(fmnc), [6.707710260366577, 1.5193083237479683], atol=0.001)
    @test isapprox(getθ(fmnc), [0.9458106880922268, 0.22692826607677266], atol=0.0001)
    @test std(fmnc)[1] ≈ [24.171449463289047, 5.799379721123582]
    @test std(fmnc)[2] ≈ [25.556130034081047]
    @test isapprox(logdet(fmnc), 74.46952585564611, atol=0.001)
    cor(fmnc)

    MixedModels.lrt(fm, fmnc)

    fmrs = fit!(LinearMixedModel(@formula(Y ~ 1 + U + (0 + U|G)), dat[:sleepstudy]))
    @test isapprox(objective(fmrs), 1774.080315280528, rtol=0.00001)
    @test isapprox(getθ(fmrs), [0.24353985679033105], rtol=0.00001)
end

#= #takes too long on Travis
@testset "d3" begin
    fm = updateL!(LinearMixedModel(@formula(Y ~ 1 + U + (1+U|G) + (1+U|H) + (1+U|I)), dat[:d3]));
    @test isapprox(pwrss(fm), 5.1261847180180885e6, rtol = 1e-6)
    @test isapprox(objective(fm), 901641.2930413672, rtol = 1e-6)
    fit!(fm)
    @test isapprox(objective(fm), 884957.5540213, rtol = 1e-6)
    @test isapprox(coef(fm), [0.4991229873, 0.31130780953], atol = 1.e-4)
end
=#

@testset "simulate!" begin
    @test MixedModels.stddevcor(cholesky!(Matrix(I, 3, 3))) == (ones(3), Matrix(I, 3, 3))
    fm = fit!(LinearMixedModel(@formula(Y ~ 1 + (1 | G)), dat[:Dyestuff]))
    refit!(simulate!(Random.MersenneTwister(1234321), fm))
    @test isapprox(deviance(fm), 339.0218639362958, atol=0.001)
    refit!(fm, dat[:Dyestuff][:Y])
    Random.seed!(1234321)
    refit!(simulate!(fm))
    @test isapprox(deviance(fm), 339.0218639362958, atol=0.001)
    simulate!(fm, θ = getθ(fm))
    @test_throws DimensionMismatch refit!(fm, zeros(29))
    Random.seed!(1234321)
    dfr = bootstrap(10, fm)
    @test size(dfr) == (10, 5)
    @test names(dfr) == Symbol[:obj, :σ, :β₁, :θ₁, :σ₁]
end

@testset "Rank deficient" begin
    Random.seed!(0)
    data = DataFrame(x = rand(100), y = rand(100), z = repeat(1:20, 5))
    data.x2 = 1.5 * data.x
    model = fit(LinearMixedModel, @formula(y ~ x + x2 + (1|z)), data)
    @test length(fixef(model)) == 2
end


text:
using Plots

include(string(pwd(), "/src/interpreter.jl"))
include(string(pwd(), "/src/gensys.jl"))
include(string(pwd(), "/galis_example.jl"))

@polyvar pi pi_e y y_e i v v_l ep

#variables = [pi pi_e y y_e v v_l ep]

#Equations of the model

eq1 = pi - bet*pi_e - kappa*y
eq2 = y + 1/sig*(i - pi_e) - y_e
eq3 = i - phi_pi*pi - phi_y*y - v_l
eq4 = v - rho_v*v_l - ep

#Dictionary with all the equations

model = Dict("IS"=>eq2,"Phillips"=>eq1,"Taylor"=>eq3,"Ar Shock"=>eq4)

modelgensys = model2gensys(model,[ep])

sol = gensys(modelgensys)

fris = irf(sol,15,0.25)

plot(4*fris[:,1],w = 3)
plot!(4*irfs_true[2:15,1], line = :dash,w = 3)

plot(fris[:,2],w = 3)
plot!(irfs_true[2:15,2], line = :dash,w = 3)

plot(4*fris[:,3],w = 3)
plot!(4*irfs_true[2:15,4], line = :dash,w = 3)

plot(fris[:,3],w = 3)
plot!(irfs_true[2:15,4], line = :dash,w = 3)

max_stars_repo_name: danmrc/DSGE-HMC

max_issues_repo_name: danmrc/DSGE-HMC

max_forks_repo_name: danmrc/DSGE-HMC


text:
module OpenTrick

export opentrick, rawio, blockingtask, unsafe_clear

const tasks_pending=Dict{Condition, Task}()

mutable struct IOWrapper{T <: IO}
    io::T
    cond::Condition
    function IOWrapper(io::T, c) where T
        obj = new{T}(io, c)
        finalizer(obj) do obj
            notify(obj.cond)
        end
    end
end

"""
    opentrick(openfn[, args... [; <keyword arguments>]])

Call `openfn` with `(handlefn, args... ,kwargs ...)` as arguments,
return an `IOWrapper` instance. (NB:`handlefn` is provided by `opentrick`.)

# Arguments
- `openfn::Function` function actually called to obtain a `IO` instance.
  `openfn` must take a `Function(::IO)` instance as its first argument
- `args` optional arguments that will be passed to `openfn`
- `kwargs` optional keyword arguments that will be passed to `openfn`

# Examples
```jldoctest
julia> using OpenTrick

julia> filename = tempname();

julia> io = opentrick(open, filename, "w+");

julia> write(io, "hello world!")
12

julia> seek(io, 0);

julia> readline(io)
"hello world!"

```
"""
function opentrick(open::Function, args...; kwargs...)
    blockreturn() do cond
        open(args ...; kwargs...) do stream
            notifyreturn(cond, stream)
        end
    end
end

"""
    close(io; kwargs...)

Close `io` and unblock the corresponding blocking task.
"""
function Base.close(w::IOWrapper; kwargs...)
    try
        close(rawio(w); kwargs...)
    catch e
        rethrow(e)
    finally
        finalize(w)
        yield()
    end
end

for fname in (
    :read, :read!, :readbytes!, :unsafe_read, :readavailable,
    :readline, :readlines, :eachline, :readchomp, :readuntil, :bytesavailable,
    :write, :unsafe_write, :truncate, :flush,
    :print, :println, :printstyled, :showerror,
    :seek, :seekstart, :seekend, :skip, :skipchars, :position,
    :mark, :unmark, :reset, :ismarked,
    :isreadonly, :iswritable, :isreadable, :isopen, :eof,
    :countlines, :displaysize)
    eval(quote
        @assert isa(Base.$fname, Function)
        Base.$fname(w::IOWrapper, args...; kwargs...) = $fname(rawio(w), args...; kwargs...)
    end)
end

"""
    unsafe_clear()

Unblock all blocking tasks. All `io`s returned by `opentrick` will
be closed as a consequence.
"""
function unsafe_clear()
    for (cond, task) in tasks_pending
        notify(cond, InterruptException(), error=true)
        yield()
    end
end

"""
    rawio(io)

Return the actual `io` instance
"""
rawio(w::IOWrapper) = w.io

"""
    blockingtask(io)

Return the task blocking which prevents the `handlefn` passed to `openfn` from returning
"""
blockingtask(w::IOWrapper) = tasks_pending[w.cond]

"""
call blockreturn in other tasks like in @async block
"""
function notifyreturn(c, x)
    io = IOWrapper(x, c)
    @debug "notifing caller..." io
    notify(c, io, all=false) # caller ought to be waiting for it
    io = nothing # no longer keep reference to
    wait(c)
    @debug "notifyreturn finished"
end

"""
f must accept cond::Condition as its first argument
"""
function blockreturn(f::Function, args...; kwargs...)
    cond = Condition()
    task = @async begin
        try
            @debug "push!" cond current_task() tasks_pending
            push!(tasks_pending, cond => current_task() )
            @debug "calling ..." f args kwargs
            f(cond, args...; kwargs...)
            @debug "call returned" f
        catch e
            @debug "Caught Exception" e
            notify(cond, e, error=true)
        end
        @debug "delete!" current_task() tasks_pending
        delete!(tasks_pending, cond)
        @debug "deleted" tasks_pending

    end
    @debug "waiting" cond task
    return wait(cond)
end

end # module


text:
using Tables
using Recommenders: Movielens1M, Movielens100k, download, load_dataset
using Test

@testset "Movielens100k" begin
    ml100k = Movielens100k()
    download(ml100k, usecache = false, removezip = true, unzip = true)
    rating, user, movie = load_dataset(ml100k)

    @test size(rating) == (100000,)
    @test Tables.columnnames(rating) == [:userid, :movieid, :rating, :timestamp]

    @test size(user) == (943,)
    @test Tables.columnnames(user) == [:userid, :age, :gender, :occupation, :zipcode]

    @test size(movie) == (1682,)
    @test Tables.columnnames(movie) == [
        :movieid,
        :movie_title,
        :release_date,
        :video_release_date,
        :IMDbURL,
        :unknown,
        :Action,
        :Adventure,
        :Animation,
        :Childrens,
        :Comedy,
        :Crime,
        :Documentary,
        :Drama,
        :Fantasy,
        :FilmNoir,
        :Horror,
        :Musical,
        :Mystery,
        :Romance,
        :SciFi,
        :Thriller,
        :War,
        :Western,
    ]
end

@testset "Movielent1M" begin
    ml1m = Movielens1M()
    download(ml1m, usecache = false, removezip = true, unzip = true)
    rating, user, movie = load_dataset(ml1m)

    @test size(rating) == (1000209,)
    @test Tables.columnnames(rating) == [:userid, :movieid, :rating, :timestamp]

    @test size(user) == (6040,)
    @test Tables.columnnames(user) == [:userid, :gender, :age, :occupation, :zipcode]
    @test min(user.userid...) == 1
    @test max(user.userid...) == 6040

    @test size(movie) == (3883,)
    @test Tables.columnnames(movie) == [:movieid, :title, :genres]
    @test min(movie.movieid...) == 1
    @test max(movie.movieid...) == 3952
end


text:
# From node.

mutable struct FromNode <: TabularNode
    table::SQLTable

    FromNode(; table) =
        new(table)
end

FromNode(table) =
    FromNode(table = table)

"""
    From(; table)
    From(table)

`From` outputs the content of a database table.

```sql
SELECT ...
FROM \$table
```

# Examples

```jldoctest
julia> person = SQLTable(:person, columns = [:person_id, :year_of_birth]);

julia> q = From(person);

julia> print(render(q))
SELECT "person_1"."person_id", "person_1"."year_of_birth"
FROM "person" AS "person_1"
```
"""
From(args...; kws...) =
    FromNode(args...; kws...) |> SQLNode

dissect(scr::Symbol, ::typeof(From), pats::Vector{Any}) =
    dissect(scr, FromNode, pats)

Base.convert(::Type{AbstractSQLNode}, table::SQLTable) =
    FromNode(table)

function PrettyPrinting.quoteof(n::FromNode, ctx::QuoteContext)
    tex = get(ctx.vars, n.table, nothing)
    if tex === nothing
        tex = quoteof(n.table, limit = true)
    end
    Expr(:call, nameof(From), tex)
end

label(n::FromNode) =
    n.table.name


text:
using QuasiGeostrophy, LinearAlgebra, Test
# Polynomial Order
n = 16

# Domain [a b]
Ω = [-1.0 1.0]
b = Ω[end]
a = Ω[1]
L = b - a

D, cx = cheb(n)

x =  (b-a) .* (cx .+ 1) ./ 2 .+ a

∂ˣ = 2 / L .* D

# Equation Dy = y, y[x=-1] = 1
A = ∂ˣ - I

# Boundary Conditions
A[end, :] .= 0.0
A[end, end] = 1.0

# Solve
rhs = 0.0 * x
rhs[end] = 1.0
y = A \ rhs

# Exact solution
solution = exp.( x .+ 1)

@testset "ODE Test 1" begin
    bool = norm(y - solution) < eps(1.0) * 100
    @test bool
end

# Polynomial Order
n = 32

# Domain [a b]
Ω = [-1.0 0.0]
b = Ω[end]
a = Ω[1]
L = b - a

D, cx = cheb(n)

x =  (b-a) .* (cx .+ 1) ./ 2 .+ a

∂ˣ = 2 / L .* D

# Equation Dy = y, y[x=-1] = 1
A = ∂ˣ - I

# Boundary Conditions
A[end, :] .= 0.0
A[end, end] = 1.0

# Solve
rhs = 0.0 * x
rhs[end] = 1.0
y = A \ rhs

# Exact solution
solution = exp.( x .+ 1)

@testset "ODE Test 2" begin
    bool = norm(y - solution) < eps(1.0) * 100
    @test bool
end

max_stars_repo_name: upiterbarg/QuasiGeostrophy

max_issues_repo_name: sandreza/QuasiGeostrophy

max_forks_repo_name: upiterbarg/QuasiGeostrophy


text:
using Test
using NeuralAttentionlib

using Random
using Flux
using NNlib
using Static
using ChainRulesCore
using ChainRulesTestUtils

const tests = [
    "collapseddim",
    "matmul",
    "mask",
    "mha",
]

Random.seed!(0)

include("old_impl/old_impl.jl")
using .Old_Impl
using .Old_Impl: batched_triu!, batched_tril!

@testset "NeuralAttentionlib" begin
    for t in tests
        fp = joinpath(dirname(@__FILE__), "$t.jl")
        @info "Test $(uppercase(t))"
        include(fp)
    end
end


text:
# unambigously
Base.convert(::Type{CombinedParser},x::CombinedParser) =
    x

"""
    Base.convert(::Type{CombinedParser},x)

[`parser`](@ref)`(x)`.
"""
Base.convert(::Type{CombinedParser},x) =
    parser(x)

"""
    parser(x)

A [`ConstantParser`](@ref) matching `x`.
"""
parser(x) =
    ConstantParser(x)

"""
    parser(x::StepRange{Char,<:Integer})

[`ValueIn`](@ref) matching x.
"""
parser(x::StepRange) =
    ValueIn(x)
parser(x::CombinedParser) = x

"""
    parser(x::Pair{Symbol, P}) where P

A parser labelled with name `x.first`.
Labels are useful in printing and logging.

See also: [`@with_names`](@ref), [`with_name`](@ref), [`log_names`](@ref)
"""
parser(x::Pair{Symbol}) =
    NamedParser(x.first, parser(x.second))


text:
push!(LOAD_PATH, "../src")

using CASIM

render_option = Not_Render()
sim_horizon = 5000
num_run = 3
init_num_ac = 100
const_num_ac = [150, 200, 250, 300, 350, 400, 450]

"""
	VICAS + Nearest
	num_ac
"""
exp = "ConstNumAC_VICAS_10km_init_ac_" * string(init_num_ac) * "/" * VICAS_policyname # to make sub-folder

if !ispath("../simulation_results/ConstNumAC_VICAS_10km_init_ac_" * string(init_num_ac))
	mkdir("../simulation_results/ConstNumAC_VICAS_10km_init_ac_" * string(init_num_ac))
end
cas = VICAS()
cas_ctrl = Uniform_CAS(cas)
sensor = Nearest_Intruder_Tracker()
sensor_ctrl = Uniform_Sensor(sensor)

for run in 1:num_run
	for num_ac in const_num_ac
		scheme = "Nearest_num_ac_" *
			string(Int(num_ac)) * "_run_" * string(run)
		airspace_ctrl = Constant_Number_Controller(num_ac)
		simulate(exp, scheme, sim_horizon, airspace_ctrl, cas_ctrl, sensor_ctrl, render_option)
	end
end

max_stars_repo_name: parachutel/DeepCorrectionCAS


text:
"""
```
selection!(norm_weights, s_t1_temp, s_t_nontemp, ϵ_t;
    resampling_method = :multinomial)
```

Resample particles using `norm_weights`. This function modifies `s_t1_temp`,
`s_t_nontemp`, and `ϵ_t` in place.
"""
function selection!(norm_weights::Vector{Float64}, s_t1_temp::M, s_t_nontemp::M, ϵ_t::M;
                    resampling_method::Symbol = :multinomial) where M<:AbstractMatrix{Float64}
    # Resampling
    is = resample(norm_weights, method = resampling_method)

    # Update arrays using resampled indices
    s_t1_temp   .= s_t1_temp[:, is]
    s_t_nontemp .= s_t_nontemp[:, is]
    ϵ_t         .= ϵ_t[:, is]

    return nothing
end

"""
```
resample(weights; method = :systematic)
```

Return indices after resampling according to `method`, weighting by `weights`.
See https://xianblog.files.wordpress.com/2017/12/lawrenz.png for a helpful
visual comparison of multinomial and systematic resampling.
"""
function resample(weights::Vector{Float64}; method::Symbol = :systematic)
    n_particles = length(weights)

    # Normalize weights if necessary
    if sum(weights) != 1
        weights = weights ./ sum(weights)
    end

    if method == :systematic
        # Divide a circle of circumference 1 into n_particles segments, where
        # the jth segment has length weights[j] and lies on the half-closed
        # interval (cdf[j-1], cdf[j]]
        cdf = cumsum(weights)

        # Take a wheel with n_particles spokes, each spaced 1/n_particles
        # apart, and shift it from its initial position (where there is a spoke
        # at 0) by offset
        offset = rand() / n_particles

        # The ith spoke landed in the jth segment <=> set the ith new particle
        # equal to the jth old particle
        new_inds = zeros(Int, n_particles)
        j = 1
        for i in 1:n_particles
            spoke = (i-1)/n_particles + offset
            # Find j such that cdf[j-1] < spoke <= cdf[j]
            while spoke > cdf[j]
                j += 1
            end
            new_inds[i] = j
        end
        return new_inds

    elseif method == :multinomial
        return sample(1:n_particles, Weights(weights), n_particles, replace = true)

    else
        throw("Invalid resampling method: $method")
    end
end


text:
function tsat_check(t)
    if !(IF97_T3 ≤ T ≤ IF97_Tc)
        return throw(DomainError(T, "Temperature not between tripple and critical points"))
    else
        return nothing
    end
end

function is_liquid(x::Symbol)::Bool
    x in (:liquid,:l,:L,:LIQUID)
end

function is_gas(x::Symbol)::Bool
    x in (:gas,:g,:G,:GAS,:vapor,:VAPOR,:V,:v,:steam,:STEAM)
end

const satvl_b = [   1.992_740_64,
1.099_653_42,
-0.510_839_303,
-1.754_934_79,
-45.517_035_2,
-6.746_944_50E5]

satvv_c = [  -2.031_502_40,
-2.683_029_40,
-5.386_264_92,
-17.299_160_5,
-44.758_658_1,
-63.920_106_3
]
function massrho_l_impl(::SingleSatT,model::IndustrialWater,t::T) where T
    t =normalize_units(T)
    tsat_check(t)
    θ = T*IF97_Tc_inv
    τ = 1 - θ
    b = satvl_b
    return ρc*(T(1) + b[1]*τ^T(1/3) + b[2]*τ^T(2/3) + b[3]*τ^T(5/3) + b[4]*τ^T(16/3)
           + b[5]*τ^T(43/3) + b[6]*τ^T(110/3))
end

function massrho_v_impl(::SingleSatT,model::IndustrialWater,t::T) where T
    t =normalize_units(T)
    tsat_check(t)
    c = satvv_c
    θ = T*IF97_Tc_inv
    τ = 1 - θ

    return ρc*exp(c[1]*τ^T(2/6) + c[2]*τ^T(4/6) + c[3]*τ^T(8/6) + c[4]*τ^T(18/6)
    + c[5]*τ^T(37/6) + c[6]*τ^T(71/6))
end

function mass_density(mt::SingleSatT,model::IndustrialWater,st::ThermodynamicState,unit="kg/m3")
    sat = phase(FromState(),st)
    t = temperature(FromState(),t)
    if is_liquid(sat)
        res = massrho_l_impl(mt,model,t)
        return convert_unit("kg/m3",unit,res)
    elseif is_gas(sat)
        res =  massrho_v_impl(mt,model,t)
        return convert_unit("kg/m3",unit,res)
    else
        throw(error("unspecified required phase. Stablish phase = :gas or phase = :liquid"))
    end
end

function mass_volume(mt::SingleSatT,model::IndustrialWater,st::ThermodynamicState,unit)
    rho = mass_density(mt,model,st)
    res = one(res)/res
    return convert_unit(u"m^3/kg",unit,res)
end

#mol_volume does the following dispatch, in the case of SingleSatT

#mol_volume(model,st) - macro
#mol_volume(mt,model,st) - macro, multiplies mass volume with molar mass
#mass_volume(mt,model,st) - normally to mass_volume_impl, here to mass_density
#mass_density(mt,model,st)

#to add total and mol units, it suffices just to add the mass unit.

"""
    Psat2

    Do not call this function, but rather use Psat().
    Returns the saturated vapour pressure according to IAPWS SR1-86(1992)
    with T in K and P in MPa.
    This replaces the Region 4 equations for Psat and is used to calcualate
    the saturated enthalpy and entropy to match the values in IAPWS SR1-86(1992)
    It is not used in Tsat(), as the difference in vapour pressure is in the
    8th decimal at the boiling point.
"""

const dpsat_a = [ -7.85951783000,
1.84408259000,
-11.78664970000,
22.68074110000,
-15.96187190000,
1.80122502000
]

function ∂Psat∂T(T)
    a = dpsat_a
    θ = T*IF97_Tc_inv
    τ = 1 - θ
    ∂τ = - IF97_Tc_inv
    ∂T = one(T)
    τhalf = sqrt(τ)
    τ15 = τ*τhalf
    τ2 = τ*τ
    τ25 = τ2*τhalf
    τ3 = τ2*τ
    τ35=τ3*τhalf
    τ4=τ2*τ2
    τ65 = τ35*τ3
    τ75 = τ35*τ4

    x = IF97_Tc*(a[1]*τ + a[2]*τ15 + a[3]*τ3 + a[4]*τ35 + a[5]*τ4 + a[6]*τ75)
    ∂x = IF97_Tc*∂τ*(a[1] + 1.5*a[2]*τhalf + 3.0*a[3]*τ2 + 3.5*a[4]*τ25 + 4.0*a[5]*τ3 + 7.5*a[6]*τ65)
    xT = x/T
    #res =  Pc*exp(xT)
    ∂res = IF97_Pc*exp(xT)*(∂x - xT*∂T)/T
end


text:
using Test
using LinearAlgebra
using SummationByPartsOperators

for T in (Float32, Float64)
    xmin = zero(T)
    xmax = one(T)
    N = 51
    D₁ = periodic_derivative_operator(1, 4, xmin, xmax, N)
    D₂ = periodic_derivative_operator(2, 4, xmin, xmax, N)
    D₃ = periodic_derivative_operator(3, 4, xmin, xmax, N)

    x = grid(D₁)
    u = sinpi.(x)
    dest1 = fill(zero(eltype(u)), length(u))
    dest2 = fill(zero(eltype(u)), length(u))
    dest3 = fill(zero(eltype(u)), length(u))

    sum_12 = D₁ + D₂
    sum_123 = D₁ + D₂ + D₃
    sum_12_3 = sum_12 + D₃
    sum_3_12 = D₃ + sum_12

    for compact in (true, false)
        show(IOContext(devnull, :compact=>false), sum_12)
        show(IOContext(devnull, :compact=>false), sum_123)
        show(IOContext(devnull, :compact=>false), sum_12_3)
        show(IOContext(devnull, :compact=>false), sum_3_12)
    end

    # Compare mul! with β=0 and mul! without β.
    mul!(dest1, sum_12, u, one(T), zero(T))
    mul!(dest2, sum_12, u, one(T))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))

    # Test commutativity.
    mul!(dest1, sum_123, u, T(2))
    mul!(dest2, sum_12_3, u, T(2))
    mul!(dest3, sum_3_12, u, T(2))
    @test all(i->dest1[i] ≈ dest2[i], eachindex(u))
    @test all(i->dest1[i] ≈ dest3[i], eachindex(u))
end


text:
"""
sample!([rng,] a, k; exclude = ())

Sample `k` element from array `a` without repetition and eventually excluding elements in `exclude`.
Pay attention, it changes the order of the elements in `a`.
"""
function sample!(rng::AbstractRNG, a::AbstractArray, k::Integer; exclude = ())
    length(a) < k + length(exclude) && error("Array too short.")
    res = Vector{eltype(a)}()
    sizehint!(res, k)
    n = length(a)
    i = 1
    while length(res) < k
        r = rand(rng, 1:n-i+1)
        if !(a[r] in exclude)
            push!(res, a[r])
            a[r], a[n-i+1] = a[n-i+1], a[r]
            i += 1
        end
    end
    res
end

sample!(a::AbstractArray, k::Integer; exclude = ()) = sample!(getRNG(), a, k; exclude = exclude)

getRNG(seed::Integer = -1) = seed >= 0 ? MersenneTwister(seed) : Base.Random.GLOBAL_RNG


text:
"""
Metadata in a FLAC stream.

Used in both encoding and decoding of a stream.
Each type of metadata object contains an indicator
of its `typ`, an indicator of whether this is the
last metadata block and its length, in bytes.
"""
abstract type StreamMetaData end

@enum(MetaDataType,
      Info,
      Padding,
      Application,
      SeekTable,
      VorbisComment,
      CueSheet,
      Picture)

@enum(PictureType,
      Other,
      FileIconStandard,
      FileIcon,
      FrontCover,
      BackCover,
      LeafletPage,
      Media,
      LeadArtist,
      Artist,
      Conductor,
      Band,
      Composer,
      Lyricist,
      RecordingLocation,
      DuringRecording,
      DuringPerformance,
      VideoScreenCapture,
      Fish,                             # ?
      Illustration,
      BandLogo,
      PublisherLogo,
      Undefined)

"""
A block containing information on the stream including
`samplerate`, `channels`, `bitspersample`, `totalsamples`,
and `mdsum`.
"""
mutable struct InfoMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Int64
    minblocksize::Int32
    maxblocksize::Int32
    minframesize::Int32
    maxframesize::Int32
    samplerate::Int32
    channels::Int32
    bitspersample::Int32
    totalsamples::Int64
    md5sum::NTuple{16,UInt8}
    InfoMetaData() =
        new(Info, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, tuple(fill(0x00, 16)...))
end

"""
`len` bytes of padding in the FLAC stream.
"""
mutable struct PaddingMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Int64
    dummy::Cint
end

"""
Application metadata

Not sure what this is for.
"""
mutable struct ApplicationMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Int64
    id::NTuple{4, Char}
    data::Ptr{Void}
end

struct SeekPoint
    sample_number::Int64
    stream_offset::Int64
    frame_samples::UInt32
end

"""
An array of `SeekPoint`s
"""
mutable struct SeekTableMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Int64
    num_points::UInt32
    points::Ptr{SeekPoint}
end

"""
A single Vorbis comment.

Comments are usually key/value pairs of the form
`ARTIST=Miles Davis`, `YEAR=1965`, etc.
"""
struct VorbisCommentEntry
    len::UInt32
    entry::Ptr{UInt8}
end

"""
Vorbis comment metadata.  The vendor comment is always present.
"""
mutable struct VorbisCommentMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Cuint
    vendor::VorbisCommentEntry
    n::UInt32
    comments::Ptr{VorbisCommentEntry}
end

"""
A single cue sheet index
"""
struct CueSheetIndex
    offset::UInt64
    number::UInt8
end

"""
A single track annotation in a CueSheet.

I'm not sure about the offsets here.  In the C struct the `typ` and `pre_emphasis` fields
are single bits.
"""
struct CueSheetTrack
    offset::Int64
    number::UInt8
    isrc::NTuple{13, UInt8}
    typ::Cuint
#    pre_emphasis::Cuint
    num_indices::UInt8
    indices::Ptr{CueSheetIndex}
end

CueSheetTrack() = unsafe_load((:FLAC__StreamMetadata_CueSheet_Track, libflac),
                              Ptr{CueSheetTrack}, (Void,))

"""
Cue sheet meta data.

An array of `CueSheetTrack`s
"""
mutable struct CueSheetMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Int64
    media_catalog_number::NTuple{129,UInt8}
    lead_in::Int64
    is_cd::Cint
    num_tracks::Cuint
    tracks::Ptr{CueSheetTrack}
end

mutable struct PictureMetaData <: StreamMetaData
    typ::MetaDataType
    is_last::Cint
    len::Int64
    ptyp::PictureType
    mime_type::Ptr{UInt8}
    description::Ptr{UInt8}
    width::UInt32
    height::UInt32
    depth::UInt32
    color::UInt32
    data_length::UInt32
    data::Ptr{UInt8}
end

const MDTypes = DataType[InfoMetaData, PaddingMetaData, ApplicationMetaData,
                         SeekTableMetaData, VorbisCommentMetaData,
                         CueSheetMetaData, PictureMetaData];

# Zero-argument external constructors for the metadata types.
#
# The finalizer is added because the memory is allocated by libflac
# The first metadatatype, InfoMetaData, is skipped because it has an internal
# zero-argument constructor.  This has the happy side-effect of aligning the tags
# with the integer values of the MetaDataType.

for (ind,typ) in enumerate(MDTypes[2:end])
    @eval begin
        $(Symbol(typ))() =
            unsafe_load(ccall((:FLAC__metadata_object_new, libflac),
                              Ptr{$typ}, (MetaDataType,), MetaDataType($ind)))
    end
end

"""
Factory to construct a subtype of StreamMetaData from an opaque pointer.

Typically this is used in a callback function that is passed a `Ptr{Void}`.
"""
metadata(pt::Ptr{StreamMetaData}) =
    unsafe_load(convert(Ptr{MDTypes[unsafe_load(convert(Ptr{Int32}, pt)) + 1]}, pt))

## Specific constructors from file names
"""
Open the file, `fnm`, check that it is a flac stream and return any cue sheets, closing the file.
"""
function CueSheetMetaData(fnm::String)
    isfile(fnm) || throw(ArgumentError(string("\"",fnm, "\" is not a path to a readable file")))
    cue = Ptr{CueSheetMetaData}[C_NULL]
    ccall((:FLAC__metadata_get_cuesheet, libflac), Bool,
          (Ptr{UInt8}, Ptr{Ptr{CueSheetMetaData}}),
          fnm, cue) || error("call to get_cuesheet failed")
    unsafe_load(cue[1])
end

"""
Open the file, `fnm`, check that it is a flac stream and return the stream info, closing the file.
"""
function InfoMetaData(fnm::String)
    isfile(fnm) || throw(ArgumentError(string("\"",fnm, "\" is not a path to a readable file")))
    strinf = InfoMetaData()
    ccall((:FLAC__metadata_get_streaminfo,libflac), Bool,
          (Ptr{UInt8}, Ref{InfoMetaData}), fnm, strinf) || error("call to get_streaminfo failed")
    strinf
end

"""
Open the file, `fnm`, check that it is a flac stream and return any Vorbis tags, closing the file.
"""
function VorbisCommentMetaData(fnm::String)
    isfile(fnm) || throw(ArgumentError(string("\"",fnm, "\" is not a path to a readable file")))
    vorb = Ptr{VorbisCommentMetaData}[C_NULL]
    ccall((:FLAC__metadata_get_tags, libflac), Bool,
          (Ptr{UInt8}, Ptr{Ptr{VorbisCommentMetaData}}),
          fnm, vorb) || error("call to get_tags failed")
    unsafe_load(vorb[1])
end

"""
Create a Ptr{StreamMetaData} from a Dict of key/value pairs.

Both the key and the value are converted to strings.
"""
function Base.convert(::Type{Ptr{StreamMetaData}}, dd::Dict{K,V}) where {K,V}
    vcp = ccall((:FLAC__metadata_object_new, libflac),
                Ptr{StreamMetaData}, (MetaDataType, ), VorbisComment)
    for (k, v) in dd
        kk = string(k)
        vv = string(v)
        ccall((:FLAC__format_vorbiscomment_entry_name_is_legal, libflac), Bool,
              (Ptr{UInt8},Cuint), kk, length(kk)) ||
              error(string("\"",k, "\" is not a legal Vorbis comment name"))
        ccall((:FLAC__format_vorbiscomment_entry_value_is_legal, libflac), Bool,
              (Ptr{UInt8},Cuint), vv, -one(Cuint))||
              error(string("\"",v,"\" is not a legal Vorbis comment value"))
        ce = [VorbisCommentEntry(0, C_NULL)]
        ccall((:FLAC__metadata_object_vorbiscomment_entry_from_name_value_pair, libflac),
              Bool, (Ptr{VorbisCommentEntry}, Ptr{UInt8}, Ptr{UInt8}),
              ce, kk, vv) || error("failure to create vorbiscomment with key $k and value $v")
        ccall((:FLAC__metadata_object_vorbiscomment_append_comment, libflac), Bool,
              (Ptr{StreamMetaData}, VorbisCommentEntry, Bool),
              vcp, ce[1], true) || error("failure to append vorbiscomment with key $k and value $v")
    end
    vcp
end

function Base.convert(::Type{Dict}, vc::VorbisCommentMetaData)
    dd = Dict{String, String}()
    dd["vendor"] = string(vc.vendor)
    for e in unsafe_wrap(Vector{FLAC.VorbisCommentMetaData}, vc.comments, vc.n)
        k = Ptr{UInt8}[C_NULL]
        v = Ptr{UInt8}[C_NULL]
        ccall((:FLAC__metadata_object_vorbiscomment_entry_to_name_value_pair, libflac),
              Bool, (VorbisCommentEntry, Ptr{Ptr{UInt8}}, Ptr{Ptr{UInt8}}),
              e, k, v) || error("failure in entry_to_name_value_pair")
        dd[unsafe_string(k[1])] = unsafe_string(v[1])
    end
    dd
end

## show methods for various types of StreamMetaData

function Base.show(io::IO, ib::InfoMetaData)
    ib.typ == Info || error("StreamInfoMetaData must have typ == $InfoType")
    println(io,"StreamInfo metadata block")
    println(io," is_last ", ib.is_last)
    println(io," len ", ib.len)
    println(io," minblocksize: ", ib.minblocksize)
    println(io," maxblocksize: ", ib.maxblocksize)
    println(io," minframesize: ", ib.minframesize)
    println(io," maxframesize: ", ib.maxframesize)
    println(io," samplerate: ", ib.samplerate)
    println(io," channels: ", ib.channels)
    println(io," bitspersample: ", ib.bitspersample)
    println(io," totalsamples: ", ib.totalsamples)
    println(io," md5sum: ", ib.md5sum)
    println(io)
end

function Base.show(io::IO, pb::PaddingMetaData)
    pb.typ == Padding || error("PaddingMetaData object's typ must be $PaddingType")
    println(io, "Flac stream padding metadata of length ", pb.len, " bytes.")
    println(io)
end

function Base.show(io::IO, st::SeekTableMetaData)
    st.typ == SeekTable || error("SeekTableType object's typ must be $SeektableType")
    println(io, "Flac seek table metadata of length ", st.len, " bytes with ",
            st.num_points, " points.")
    pts = unsafe_wrap(convert(Ptr{SeekPoint}, st.points), st.num_points)
    println(io, pts)
end

Base.string(e::VorbisCommentEntry) = unsafe_string(e.entry, e.len)

Base.show(io::IO,e::VorbisCommentEntry) = println(io, string(e.entry))

function Base.show(io::IO,c::VorbisCommentMetaData)
    c.typ == VorbisComment || error("VorbisCommentMetaData's typ must be $VorbisComment")
    println(io,"vendor=", string(c.vendor))
    for cc in unsafe_wrap(Vector{VorbisCommentEntry}, c.comments, c.n)
        println(io, string(cc))
    end
end


text:
function each_cons(n, collection)
    [collection[i:(i + n - 1)] for i = 1:(length(collection) - n + 1)]
end

n = ["73167176531330624919225119674426574742355349194934",
     "96983520312774506326239578318016984801869478851843",
     "85861560789112949495459501737958331952853208805511",
     "12540698747158523863050715693290963295227443043557",
     "66896648950445244523161731856403098711121722383113",
     "62229893423380308135336276614282806444486645238749",
     "30358907296290491560440772390713810515859307960866",
     "70172427121883998797908792274921901699720888093776",
     "65727333001053367881220235421809751254540594752243",
     "52584907711670556013604839586446706324415722155397",
     "53697817977846174064955149290862569321978468622482",
     "83972241375657056057490261407972968652414535100474",
     "82166370484403199890008895243450658541227588666881",
     "16427171479924442928230863465674813919123162824586",
     "17866458359124566529476545682848912883142607690042",
     "24219022671055626321111109370544217506941658960408",
     "07198403850962455444362981230987879927244284909188",
     "84580156166097919133875499200524063689912560717606",
     "05886116467109405077541002256983155200055935729725",
     "71636269561882670428252483600823257530420752963450"]

println(maximum(map(digits -> prod(map(digit -> parse(Int, digit), split(digits, ""))), each_cons(5, reduce(*, n)))))

max_stars_repo_name: seaneshbaugh/rosetta-euler

max_issues_repo_name: seaneshbaugh/rosetta-euler

max_forks_repo_name: seaneshbaugh/rosetta-euler


text:
function checkdiff(f, δf, x0...)
    x = [x0...]
    y0 = f(x...)
    @assert length(y0) == 1 "Scalar functions only"
    y, ∇f = δf(x...)
    @assert isapprox(y0, y) "Return values do not match"
    #    @assert typeof(y0) === typeof(y) "Return type doesn't match"
    ∂x = ∇f()
    #    @assert typeof(x0) === typeof(∂x) "Gradient type doesn't match"
    checkgrad(f, x, ∂x)
end

function checkgrad(f, x, ∂x, h = 1e-8, atol = 1e-8)
    for k = 1:length(x)
        ∂xₖ = length(x) == 1 ? ∂x : ∂x[k]
        if isa(x[k], AbstractFloat)
            x1 = deepcopy(x)
            x1[k] -= h
            y1 = f(x1...)
            x2 = deepcopy(x)
            x2[k] += h
            y2 = f(x2...)
            isapprox(2h * ∂xₖ, y2-y1, atol=atol) || error("gradient for argument #$k doesn't match by $((y2-y1)/(2h*∂xₖ))")
        elseif isa(x[k], AbstractArray)
            for l = eachindex(x[k])
                x1 = deepcopy(x)
                x1[k][l] -= h
                y1 = f(x1...)
                x2 = deepcopy(x)
                x2[k][l] += h
                y2 = f(x2...)
                isapprox(2h * ∂xₖ[l], y2-y1, atol=atol) || error("gradient for argument #$k element $l doesn't match by $((y2-y1)/(2h*∂xₖ[l]))")
            end
        else error("not supported argument #$k type: $(typeof(x[k]))")
        end
    end
    true
end