text:
@doc raw"""
    SkewHermitianMatrices{n,𝔽} <: AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding}

The [`AbstractManifold`](@ref) $ \operatorname{SkewHerm}(n)$ consisting of the real- or
complex-valued skew-hermitian matrices of size ``n × n``, i.e. the set

````math
\operatorname{SkewHerm}(n) = \bigl\{p  ∈ 𝔽^{n × n}\ \big|\ p^{\mathrm{H}} = -p \bigr\},
````
where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transpose,
and the field $𝔽 ∈ \{ ℝ, ℂ, ℍ\}$.

Though it is slightly redundant, usually the matrices are stored as ``n × n`` arrays.

Note that in this representation, the real-valued part of the diagonal must be zero,
which is also reflected in the
[`manifold_dimension`](@ref manifold_dimension(::SkewHermitianMatrices{N,𝔽}) where {N,𝔽}).

# Constructor

    SkewHermitianMatrices(n::Int, field::AbstractNumbers=ℝ)

Generate the manifold of ``n × n`` skew-hermitian matrices.
"""
struct SkewHermitianMatrices{n,𝔽} <:
       AbstractEmbeddedManifold{𝔽,TransparentIsometricEmbedding} end

function SkewHermitianMatrices(n::Int, field::AbstractNumbers=ℝ)
    return SkewHermitianMatrices{n,field}()
end

@doc raw"""
    SkewSymmetricMatrices{n}

Generate the manifold of ``n × n`` real skew-symmetric matrices.
This is equivalent to [`SkewHermitianMatrices(n, ℝ)`](@ref).

# Constructor

    SkewSymmetricMatrices(n::Int)
"""
const SkewSymmetricMatrices{n} = SkewHermitianMatrices{n,ℝ}

SkewSymmetricMatrices(n::Int) = SkewSymmetricMatrices{n}()
@deprecate SkewSymmetricMatrices(n::Int, 𝔽) SkewHermitianMatrices(n, 𝔽)

function allocation_promotion_function(
    ::SkewHermitianMatrices{<:Any,ℂ},
    ::typeof(get_vector),
    args::Tuple,
)
    return complex
end

@doc raw"""
    check_point(M::SkewHermitianMatrices{n,𝔽}, p; kwargs...)

Check whether `p` is a valid manifold point on the [`SkewHermitianMatrices`](@ref) `M`, i.e.
whether `p` is a skew-hermitian matrix of size `(n,n)` with values from the corresponding
[`AbstractNumbers`](@ref) `𝔽`.

The tolerance for the skew-symmetry of `p` can be set using `kwargs...`.
"""
function check_point(M::SkewHermitianMatrices{n,𝔽}, p; kwargs...) where {n,𝔽}
    mpv = check_point(decorated_manifold(M), p; kwargs...)
    mpv === nothing || return mpv
    if !isapprox(p, -p'; kwargs...)
        return DomainError(
            norm(p + p'),
            "The point $(p) does not lie on $M, since it is not skew-hermitian.",
        )
    end
    return nothing
end

"""
    check_vector(M::SkewHermitianMatrices{n}, p, X; kwargs... )

Check whether `X` is a tangent vector to manifold point `p` on the
[`SkewHermitianMatrices`](@ref) `M`, i.e. `X` must be a skew-hermitian matrix of size `(n,n)`
and its values have to be from the correct [`AbstractNumbers`](@ref).
The tolerance for the skew-symmetry of `p` and `X` can be set using `kwargs...`.
"""
function check_vector(M::SkewHermitianMatrices, p, X; kwargs...)
    return check_point(M, X; kwargs...)  # manifold is its own tangent space
end

decorated_manifold(M::SkewHermitianMatrices{N,𝔽}) where {N,𝔽} = Euclidean(N, N; field=𝔽)

function get_basis(M::SkewHermitianMatrices, p, B::DiagonalizingOrthonormalBasis)
    Ξ = get_basis(M, p, DefaultOrthonormalBasis()).data
    κ = zeros(real(eltype(p)), manifold_dimension(M))
    return CachedBasis(B, κ, Ξ)
end

function get_coordinates!(
    M::SkewSymmetricMatrices{N},
    Y,
    p,
    X,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
) where {N}
    dim = manifold_dimension(M)
    @assert size(Y) == (dim,)
    @assert size(X) == (N, N)
    @assert dim == div(N * (N - 1), 2)
    k = 1
    for i in 1:N, j in (i + 1):N
        @inbounds Y[k] = X[i, j] * sqrt(2)
        k += 1
    end
    return Y
end
function get_coordinates!(
    M::SkewHermitianMatrices{N,ℂ},
    Y,
    p,
    X,
    ::DefaultOrthonormalBasis{ℂ,TangentSpaceType},
) where {N}
    dim = manifold_dimension(M)
    @assert size(Y) == (dim,)
    @assert size(X) == (N, N)
    @assert dim == N^2
    k = 1
    for i in 1:N, j in i:N
        if i == j # real-part zero on the diagonal -> just one basis vector per diag entry
            @inbounds Y[k] = imag(X[i, j])
            k += 1
        else
            @inbounds Y[k] = real(X[i, j]) * sqrt(2)
            @inbounds Y[k + 1] = imag(X[i, j]) * sqrt(2)
            k += 2
        end
    end
    return Y
end

function get_vector!(
    M::SkewSymmetricMatrices{N},
    Y,
    p,
    X,
    ::DefaultOrthonormalBasis{ℝ,TangentSpaceType},
) where {N}
    dim = manifold_dimension(M)
    @assert size(X) == (dim,)
    @assert size(Y) == (N, N)
    k = 1
    for i in 1:N
        Y[i, i] = convert(eltype(p), 0.0)
    end
    for i in 1:N, j in (i + 1):N
        @inbounds Y[i, j] = X[k] / sqrt(2)
        @inbounds Y[j, i] = -X[k] / sqrt(2)
        k += 1
    end
    return Y
end
function get_vector!(
    M::SkewHermitianMatrices{N,ℂ},
    Y,
    p,
    X,
    ::DefaultOrthonormalBasis{ℂ,TangentSpaceType},
) where {N}
    dim = manifold_dimension(M)
    @assert size(X) == (dim,)
    @assert size(Y) == (N, N)
    k = 1
    for i in 1:N, j in i:N
        if i == j # real zero on the diag
            @inbounds Y[i, j] = X[k] * 1im
            k += 1
        else
            @inbounds Y[i, j] = (X[k] + X[k + 1] * 1im) / sqrt(2)
            k += 2
            @inbounds Y[j, i] = -conj(Y[i, j])
        end
    end
    return Y
end

@doc raw"""
    manifold_dimension(M::SkewHermitianMatrices{n,𝔽})

Return the dimension of the [`SkewHermitianMatrices`](@ref) matrix `M` over the number
system `𝔽`, i.e.

````math
\dim \mathrm{SkewHerm}(n,ℝ) = \frac{n(n+1)}{2} \dim_ℝ 𝔽 - n,
````

where ``\dim_ℝ 𝔽`` is the [`real_dimension`](@ref) of ``𝔽``. The first term corresponds to
only the upper triangular elements of the matrix being unique, and the second term
corresponds to the constraint that the real part of the diagonal be zero.
"""
function manifold_dimension(::SkewHermitianMatrices{N,𝔽}) where {N,𝔽}
    return div(N * (N + 1), 2) * real_dimension(𝔽) - N
end

@doc raw"""
    project(M::SkewHermitianMatrices, p)

Projects `p` from the embedding onto the [`SkewHermitianMatrices`](@ref) `M`, i.e.

````math
\operatorname{proj}_{\operatorname{SkewHerm}(n)}(p) = \frac{1}{2} \bigl( p - p^{\mathrm{H}} \bigr),
````

where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
"""
project(::SkewHermitianMatrices, ::Any)

function project!(M::SkewHermitianMatrices, q, p)
    q .= (p .- p') ./ 2
    return q
end

@doc raw"""
    project(M::SkewHermitianMatrices, p, X)

Project the matrix `X` onto the tangent space at `p` on the [`SkewHermitianMatrices`](@ref) `M`,

````math
\operatorname{proj}_p(X) = \frac{1}{2} \bigl( X - X^{\mathrm{H}} \bigr),
````

where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transposed.
"""
project(::SkewHermitianMatrices, ::Any, ::Any)

project!(M::SkewHermitianMatrices, Y, p, X) = project!(M, Y, X)

function Base.show(io::IO, ::SkewHermitianMatrices{n,F}) where {n,F}
    return print(io, "SkewHermitianMatrices($(n), $(F))")
end
function Base.show(io::IO, ::SkewSymmetricMatrices{n}) where {n}
    return print(io, "SkewSymmetricMatrices($(n))")
end


text:
# ---
# title: 101. Symmetric Tree
# id: problem101
# author: Tian Jun
# date: 2020-10-31
# difficulty: Easy
# categories: Tree, Depth-first Search, Breadth-first Search
# link: <https://leetcode.com/problems/symmetric-tree/description/>
# hidden: true
# ---
#
# Given a binary tree, check whether it is a mirror of itself (ie, symmetric
# around its center).
#
# For example, this binary tree `[1,2,2,3,4,4,3]` is symmetric:
#
#
#
#         1
#        / \
#       2   2
#      / \ / \
#     3  4 4  3
#
#
#
#
# But the following `[1,2,2,null,3,null,3]` is not:
#
#
#
#         1
#        / \
#       2   2
#        \   \
#        3    3
#
#
#
#
# **Follow up:** Solve it both recursively and iteratively.
#
#
## @lc code=start
using LeetCode

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

ext: jl


text:
#
# types.jl -
#
# Types defined for Nüvü Camēras.
#

mutable struct Camera
    handle::NC.Cam # opaque handle
    function Camera(handle::NC.Cam)
        obj = new(handle)
        if handle.ptr != C_NULL
            finalizer(obj, close)
        end
        return obj
    end
end

mutable struct FrameGrabber
    handle::NC.Grab # opaque handle
    function FrameGrabber(handle::NC.Grab)
        obj = new(handle)
        if handle.ptr != C_NULL
            finalizer(obj, close)
        end
        return obj
    end
end

mutable struct ControllerList
    handle::NC.CtrlList # opaque handle
    cnt::Int
    function ControllerList(handle::NC.CtrlList)
        obj = new(handle, 0)
        if handle.ptr != C_NULL
            finalizer(obj, close)
            obj.cnt = NC.getSize(handle)
        end
        return obj
    end
end


text:
module TestSpatialDiscretization

using Modia
using ModiaMath.plot

const n=5
@model SpatialDiscretization begin
  u = Var(start=zeros(n))
  c = 1
  @equations begin
#    der(u[2:n]) = ( u[2:n] - u[1:n-1] ) / c
    der(u) = [1; ( u[2:n] - u[1:n-1] ) / c]
  end
end

result = simulate(SpatialDiscretization, 1)
plot(result, "u")

end


text:
@testset "CellValues" begin
for (func_interpol, quad_rule) in  (
                                    (Lagrange{1, RefCube, 1}(), QuadratureRule{1, RefCube}(2)),
                                    (Lagrange{1, RefCube, 2}(), QuadratureRule{1, RefCube}(2)),
                                    (Lagrange{2, RefCube, 1}(), QuadratureRule{2, RefCube}(2)),
                                    (Lagrange{2, RefCube, 2}(), QuadratureRule{2, RefCube}(2)),
                                    (Lagrange{2, RefTetrahedron, 1}(), QuadratureRule{2, RefTetrahedron}(2)),
                                    (Lagrange{2, RefTetrahedron, 2}(), QuadratureRule{2, RefTetrahedron}(2)),
                                    (Lagrange{3, RefCube, 1}(), QuadratureRule{3, RefCube}(2)),
                                    (Serendipity{2, RefCube, 2}(), QuadratureRule{2, RefCube}(2)),
                                    (Lagrange{3, RefTetrahedron, 1}(), QuadratureRule{3, RefTetrahedron}(2)),
                                    (Lagrange{3, RefTetrahedron, 2}(), QuadratureRule{3, RefTetrahedron}(2))
                                   )

    for fe_valtype in (CellScalarValues, CellVectorValues)
        cv = fe_valtype(quad_rule, func_interpol)
        ndim = JuAFEM.getdim(func_interpol)
        n_basefuncs = getnbasefunctions(func_interpol)

        fe_valtype == CellScalarValues && @test getnbasefunctions(cv) == n_basefuncs
        fe_valtype == CellVectorValues && @test getnbasefunctions(cv) == n_basefuncs * JuAFEM.getdim(func_interpol)

        x, n = valid_coordinates_and_normals(func_interpol)
        reinit!(cv, x)

        # We test this by applying a given deformation gradient on all the nodes.
        # Since this is a linear deformation we should get back the exact values
        # from the interpolation.
        u = Vec{ndim, Float64}[zero(Tensor{1,ndim}) for i in 1:n_basefuncs]
        u_scal = zeros(n_basefuncs)
        H = rand(Tensor{2, ndim})
        V = rand(Tensor{1, ndim})
        for i in 1:n_basefuncs
            u[i] = H ⋅ x[i]
            u_scal[i] = V ⋅ x[i]
        end
        u_vector = reinterpret(Float64, u)

        for i in 1:length(getpoints(quad_rule))
            @test function_gradient(cv, i, u) ≈ H
            @test function_symmetric_gradient(cv, i, u) ≈ 0.5(H + H')
            @test function_divergence(cv, i, u) ≈ tr(H)
            ndim == 3 && @test function_curl(cv, i, u) ≈ JuAFEM.curl(H)
            function_value(cv, i, u)
            if isa(cv, CellScalarValues)
                @test function_gradient(cv, i, u_scal) ≈ V
                function_value(cv, i, u_scal)
            elseif isa(cv, CellVectorValues)
                @test function_gradient(cv, i, u_vector) ≈ function_gradient(cv, i, u) ≈ H
                @test function_value(cv, i, u_vector) ≈ function_value(cv, i, u)
                @test function_divergence(cv, i, u_vector) ≈ function_divergence(cv, i, u) ≈ tr(H)
                ndim == 3 && @test function_curl(cv, i, u_vector) ≈ JuAFEM.curl(H)
            end
        end

        # Test of volume
        vol = 0.0
        for i in 1:getnquadpoints(cv)
            vol += getdetJdV(cv,i)
        end
        @test vol ≈ calculate_volume(func_interpol, x)

        # Test quadrature rule after reinit! with ref. coords
        x = JuAFEM.reference_coordinates(func_interpol)
        reinit!(cv, x)
        vol = 0.0
        for i in 1:getnquadpoints(cv)
            vol += getdetJdV(cv,i)
        end
        @test vol ≈ reference_volume(func_interpol)

        # Test spatial coordinate (after reinit with ref.coords we should get back the quad_points)
        for (i, qp_x) in enumerate(getpoints(quad_rule))
            @test spatial_coordinate(cv, i, x) ≈ qp_x
        end

        # test copy
        cvc = copy(cv)
        for fname in fieldnames(typeof(cv))
            @test typeof(cv) == typeof(cvc)
            @test pointer(getfield(cv, fname)) != pointer(getfield(cvc, fname))
            @test getfield(cv, fname) == getfield(cvc, fname)
        end
    end
end

@testset "#265" begin
    dim = 1
    deg = 1
    grid = generate_grid(Line, (2,))
    ip_fe = Lagrange{dim, RefCube, deg}()
    dh = DofHandler(grid)
    push!(dh, :u, 1, ip_fe)
    close!(dh);
    cell = first(CellIterator(dh))
    ip_geo = Lagrange{dim, RefCube, 2}()
    qr = QuadratureRule{dim, RefCube}(deg+1)
    cv = CellScalarValues(qr, ip_fe, ip_geo)
    res = @test_throws ArgumentError reinit!(cv, cell)
    @test occursin("#265", res.value.msg)
    ip_geo = Lagrange{dim, RefCube, 1}()
    cv = CellScalarValues(qr, ip_fe, ip_geo)
    reinit!(cv, cell)
end

end # of testset


text:
export EPITrajectory, epiNodes, epiDensity

"""
    EPITrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , EPI_factor::Int64=1
                  , profileOffset= :equispaced
                  , kargs...)

returns a 2d cartesian trajectory.

# Arguments
* `numProfiles::Int64`            - number of profiles
* `numSamplingPerProfile::Int64`  - number of sampling points per profile
* (`TE::Float64=0.0`)             - echo time in s
* (`AQ::Float64=1.e-3`)           - readout duration in s (per profile)
* (`EPI_factor::Int64=1`)         - EPI factor, e.g. how many profiles to acquire per shot
* (`profileOffset= :equispaced`)  - equispaced or random ordering of the shots
"""
function EPITrajectory(numProfiles, numSamplingPerProfile
                  ; TE::Float64=0.0
                  , AQ::Float64=1.e-3
                  , EPI_factor::Int64=1
                  , profileOffset= :equispaced
                  , kargs...)
  nodes = epiNodes(numProfiles, numSamplingPerProfile; EPI_factor=EPI_factor, profileOffset=profileOffset)
  times = readoutTimes(numProfiles, numSamplingPerProfile; TE=TE, AQ=AQ)
  return Trajectory("EPI", nodes, times, TE, AQ, numProfiles, numSamplingPerProfile, 1, true, false)
end

function epiNodes(numProfiles::Int64
                    , numSamplingPerProfile::Int64
                    ; EPI_factor::Int64=1
                    , profileOffset= :equispaced
                    , kargs...)
  nodes = zeros(2,numSamplingPerProfile, numProfiles)

  if profileOffset == :equispaced
      samplesperline = numSamplingPerProfile/EPI_factor -1
      for l = 1:numProfiles
        for gamma1 = 1:EPI_factor
          for gamma2 = 1:numSamplingPerProfile/EPI_factor
            # Node index
            index = Int64( (gamma1-1) * (numSamplingPerProfile/EPI_factor) + gamma2 )
            # Nodes containing kx values, second index is the actual sample index, l profile index
            nodes[1,index,l] =  (-1)^(gamma1-1) * (gamma2-1) / samplesperline + ((gamma1-1) % 2 ) - 0.5
            # Nodes containing ky values, second index is the actual sample index, l profile index
            nodes[2,index,l] = ( (l-1) + (gamma1-1)*numProfiles)/(numProfiles*EPI_factor) - 0.5
          end
        end
      end
  elseif profileOffset == :random
      samplesperline = numSamplingPerProfile/EPI_factor -1
      for l = 1:numProfiles
        # Sorting the random samples
        angleOffsets = sort(rand(EPI_factor)-0.5)
        for gamma1 = 1:EPI_factor
          for gamma2 = 1:numSamplingPerProfile/EPI_factor
            # Node index
            index = Int64( (gamma1-1) * (numSamplingPerProfile/EPI_factor) + gamma2 )
            # Nodes containing kx values, second index is the actual sample index, l profile index
            nodes[1,index,l] =  (-1)^(gamma1-1) * (gamma2-1) / samplesperline + ((gamma1-1) % 2 ) - 0.5
            # Nodes containing ky values, second index is the actual sample index, l profile index
            nodes[2,index,l] = angleOffsets[gamma1]
          end
        end
      end
  end
  return reshape(nodes, 2, numSamplingPerProfile*numProfiles)
end

function epiDensity(numProfiles::Int64, numSamplingPerProfile::Int64)
  @error "not implement yet"
end


text:
using ApproxFun, Random, InfiniteArrays, Test
import ApproxFunBase: testbandedblockbandedoperator, blocklengths

@time include("ReadmeTest.jl")
@time include("ExtrasTest.jl")
@time include("NumberTypeTest.jl")
# @time include("FractionalTest.jl")

@testset "Chebyshev and Fourier" begin
    @test norm(Fun(x->Fun(cos,Fourier(-π .. π),20)(x),20)-Fun(cos,20)) <100eps()
    @test norm(Fun(x->Fun(cos,Fourier(-π .. π))(x))-Fun(cos)) <100eps()
    @test norm(Fun(x->Fun(cos,Laurent)(x))-Fun(cos)) <100eps()
end

@testset "Piecewise space definite integral" begin
    Γ=Segment(-im,1.0-im) ∪ Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6))) ∪ Circle(2.0,0.2)
        z=Fun(Γ)

    S=PiecewiseSpace(map(d->isa(d,Circle) ? Fourier(d) : JacobiWeight(0.5,0.5,Ultraspherical(1,d)),components(Γ)))

    B=DefiniteLineIntegral(S)

    Random.seed!(0)
    f=Fun(S,rand(20))
    @test B*f ≈ linesum(component(f,1)) + linesum(component(f,2)) + linesum(component(f,3))
end

@testset "Extending function" begin
    Γ = Segment(-im,1.0-im) ∪ Curve(Fun(x->exp(0.8im)*(x+x^2-1+im*(x-4x^3+x^4)/6))) ∪ Circle(2.0,0.2)

    @test isempty(component(Γ,1)\component(Γ,1))
    @test Γ \ component(Γ,1) == component(Γ,2) ∪ component(Γ,3)

    @test norm(Fun(ones(component(Γ,1)),Γ) - Fun(x->x ∈ component(Γ,1) ? 1.0 : 0.0,Γ)) == 0
end

@testset "periodic x interval" begin
    dθ=PeriodicSegment(-2.,2.);dt=Interval(0,1.)
    d=dθ×dt
    Dθ=Derivative(d,[1,0]);Dt=Derivative(d,[0,1])
    u0=Fun(θ->exp(-20θ^2),dθ)

    ε = 0.1
    @time u=\([I⊗ldirichlet(dt);Dt-ε*Dθ^2-Dθ],[u0;0.];tolerance=1E-4)
    @test ≈(u(0.1,0.2),0.3103472600253807;atol=1E-2)

    A=Dt+Dθ
    testbandedblockbandedoperator(A)

    @time u=\([I⊗ldirichlet(dt);Dt+Dθ],[u0;0.0];tolerance=1E-6)
    @test ≈(u(0.2,0.1),u0(0.1);atol=1E-6)
end

@testset "Laplace in a strip" begin
    d=PeriodicSegment() × ChebyshevInterval()
    g=Fun((x,y)->real(cos(x+im*y)),∂(d))  # boundary data

    @test g(0.1,1.0) ≈ real(cos(0.1+im))
    @test g(0.1,-1.0) ≈ real(cos(0.1-im))
    v=[g;0]
    @test v(0.1,-1) ≈ [real(cos(0.1-im));0]

    A=[Dirichlet(d);Laplacian(d)]
    a = space(v)
    b = rangespace(A)

    @test Fun(component(v[1],1), component(b[1],1))(0.1,-1.0) ≈ v(0.1,-1.0)[1]
    @test Fun(component(v[1],2), component(b[1],2))(0.1,-1.0) ≈ v(0.1,-1.0)[1]
    @test ApproxFun.default_Fun(v[1] , b[1])(0.1,1.0) ≈ v(0.1,1.0)[1]
end

@testset "Bug in Multiplication" begin
    dom = Interval(0.001, 1) × PeriodicSegment(-pi, pi)

    @test blocklengths(Space(dom)) == 2:2:∞

    r,r2 = Fun((r,t) -> [r;r^2], dom)

    @test r(0.1,0.2) ≈ 0.1
    @test r2(0.1,0.2) ≈ 0.1^2

    sp = Space(dom)
    Dr = Derivative(sp, [1,0])
    @test ApproxFun.blockbandwidths(Dr) == (-1,1)
    @test ApproxFun.subblockbandwidths(Dr)  == (1,3)

    Dθ = Derivative(sp, [0,1])
    Mr = Multiplication(Fun( (r, θ) -> r, sp ), sp)
    rDr = Mr * Dr

    testbandedblockbandedoperator(rDr)
end

@testset "Periodic x Interval" begin
    d=PeriodicSegment() × ChebyshevInterval()

    u_ex=Fun((x,y)->real(cos(x+im*y)),d)
    @test u_ex(1.0,0.1) ≈ real(cos(1.0+im*0.1)) atol=10eps()

    B=Dirichlet(Space(d))

    @test B.order == 0  # tests stupid bug
    g=Fun((x,y)->real(cos(x+im*y)),rangespace(B))  # boundary data

    @test norm((B*u_ex-g).coefficients) < 100eps()

    testbandedblockbandedoperator(Laplacian(d))

    @time u=[B;Laplacian(d)]\[g;0.]

    @test u(.1,.2) ≈ real(cos(.1+.2im))
end

@testset "Mix Fourier-Chebyshev (#602)" begin
    s = Chebyshev(-π..π)
    a = Fun(t-> 1+sin(cos(2t)), s)
    L = Derivative() + a
    f = Fun(t->exp(sin(10t)), s)
    B = periodic(s,0)
    @time uChebyshev = [B;L] \ [0.;f]

    s = Fourier(-π..π)
    a = Fun(t-> 1+sin(cos(2t)), s)
    L = Derivative() + a
    f = Fun(t->exp(sin(10t)), s)
    @time uFourier = L\f

    @test norm(uFourier-uChebyshev) ≤ 100eps()
end

@testset "Conversion" begin
    f=Fun(t->[cos(t) 0;sin(t) 1],-π..π)
    g=Fun(f,Space(PeriodicSegment(-π,π)))
    @test g(.1) ≈ f(.1)
end

@testset "definite integral" begin
    Σ = DefiniteIntegral()

    f1 = Fun(t->cos(cos(t)),-π..π)
    f = Fun(t->cos(cos(t)),Laurent(-π..π))

    @test sum(f1) ≈ Σ*f
end

@testset "Sampling" begin
    ff=(x,y)->(x-y)^2*exp(-x^2/2-y^2/2)

    f=Fun(ff,Domain(-4..4)^2)
    r=ApproxFun.sample(f,5000)

    #We can compare the histogram to the 1-point correlation
    g=sum(f,1)/sum(f)
    @test  g(0.1) ≈ 0.2004758624973169

    # check bisection inv
    f = Fun(x -> exp(-x^2/2),-5..5)
    g = cumsum(f)
    @test g(ApproxFun.bisectioninv(g,0.5)) ≈ 0.5
end

@testset "piecewise sample (#635)" begin
    f = abs(Fun(sin, -5..5))
    @test integrate(f)(-4.0) ≈ -(cos(-4.0) - cos(-5.0))
    @test -(cos(-π) - cos(-5.0)) + cos(-3.0) - cos(-π) ≈ integrate(f)(-3.0)
    r = ApproxFun.sample(f,10)
    @test maximum(r) ≤ 5
    @test minimum(r) ≥ -5
end

f=Fun(exp)
x=sample(f,100000)
x=sample(f,100000)
@time x=sample(f,100000)
println("Sample: Time should be ~0.13")

@testset "Null space" begin
    d=ChebyshevInterval()
    D=Derivative(d)
    A=D^2-I
    @time κ=nullspace(A)
    @test length(κ) == 2

    c=[κ(0.);κ'(0.)]\[exp(0.),exp(0.)]
    u=(κ*c)[1]

    @test u(1.0) ≈ ℯ

    d=(-50..5.)
    x=Fun(identity,d)
    D=Derivative(d)
    @time u=nullspace(D^2-x)
    c=[u(leftendpoint(d)); u(rightendpoint(d))]\[airyai.(endpoints(d))...]
    @test norm((u*c)[1]-Fun(airyai,d))<10000eps()
end


text:
@recipe function f(fc::FORECAST;
                   titlefont  = font(8, "Courier"),
                   xtickfont  = font(3, "Courier"),
                   ytickfont  = font(3, "Courier"),
                   yguidefont = font(5, "Courier"),
                   xguidefont = font(5, "Courier"),
                   legendfont = font(5, "Courier"),
                   xlims = nothing,
                   ylims = nothing,
                   markershape = :circle,
                   markersize = 2,
                   markerstrokewidth = 0,
                   title = fc.call,
                   xguide = "Time",
                   legend = nothing)

    # Default values
    titlefont  := titlefont
    xtickfont  := xtickfont
    ytickfont  := ytickfont
    yguidefont := yguidefont
    xguidefont := xguidefont
    legendfont := legendfont
    title      := title
    xguide     := xguide
    markershape := markershape
    markersize := markersize
    markerstrokewidth := markerstrokewidth

    # Making the forecast to occupy the last 1/4 of xlims
    sx = size(fc.model.x,1)
    sm = size(fc.mean,1)
    ts = vcat(fc.model.x[:,1],fc.mean[:,1])
    xr1 = max(1,sx-3*size(fc.mean,1))
    xr2 = sx+sm
    xr = eltype(ts) in [Date, DateTime, Time] ?
        (ts[xr1], ts[xr2]) : (Dates.value(ts[xr1]), Dates.value(ts[xr2]))
    xlims := isnothing(xlims) ? xr : xlims
    xend = Array(fc.model.x[xr1:end,2:end])
    my = minimum([minimum(xend),minimum(Array(fc.lower[:,2:end]))])
    My = maximum([maximum(xend),maximum(Array(fc.upper[:,2:end]))])
    ypad = 0.1 * maximum([abs(My),abs(my)])
    ylims := isnothing(ylims) ? (my-ypad,My+ypad) : ylims

    xts = fc.model.x
    names_x = names(xts)[2:end]
    x = Array(xts[:,2:end])
    ts = xts[:,1]
    name_ts = names(xts)[1]

    legend := isnothing(legend) ? (size(x,2) <= 4 ? :top : :outerright) : legend

    @series begin
        label := permutedims(fc.model.varnames)
        color_palette := :tab10
        seriescolor := reshape(collect(1:size(x,2)),1,size(x,2))
        seriestype := :line

        ts,x
    end

    alpha_dimming = 2
    if size(x,2) == 1
        @series begin
            label := nothing
            color_palette := :tab10
            seriescolor := reshape(collect(1:size(x,2)),1,size(x,2))
            linestyle := :dot
            seriesalpha := 0
            fillalpha := 1/(alpha_dimming*size(x,2))

            rmean = Array(fc.mean[:,2:end])
            lower1 = Array(fc.lower[:,2:div(end-1,2)+1])
            upper1 = Array(fc.upper[:,2:div(end-1,2)+1])

            ribbon :=  (rmean - lower1, upper1 - rmean)

            fc.mean[:,1], Array(fc.mean[:,2:end])

        end
    end

    @series begin
        label := nothing
        color_palette := :tab10
        seriescolor := reshape(collect(1:size(x,2)),1,size(x,2))
        seriestype := :line
        linestyle := :dot
        seriesalpha := 0
        fillalpha := 1/(alpha_dimming*size(x,2))

        rmean = Array(fc.mean[:,2:end])
        lower2 = Array(fc.lower[:,div(end-1,2)+2:end])
        upper2 = Array(fc.upper[:,div(end-1,2)+2:end])

        ribbon :=  (rmean - lower2, upper2 - rmean)

        fc.mean[:,1], Array(fc.mean[:,2:end])
    end

    @series begin

        label := nothing
        color_palette := :tab10
        seriescolor := reshape(collect(1:size(x,2)),1,size(x,2))
        seriestype := :line
        linestyle := :dot
        seriesalpha := 1

        fc.mean[:,1], Array(fc.mean[:,2:end])
    end

end


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function _eqpole{T<:AbstractFloat}(l::T, b::T, southpole::Bool)
    sgn = southpole ? -1.0 : 1.0
    l = deg2rad(sgn*l)
    b = deg2rad(sgn*b)
    sq = 2.0*(1.0 - sin(b))
    sq = sq >= 0.0 ? sq : 0.0
    r = 18.0*3.53553391*sqrt(sq)
    return r*sin(l), r*cos(l)
end

"""
    eqpole(l, b) -> x, y

### Purpose ###

Convert right ascension \$l\$ and declination \$b\$ to coordinate \$(x, y)\$
using an equal-area polar projection.

### Explanation ###

The output \$x\$ and \$y\$ coordinates are scaled to be in the range \$[-90, 90]\$
and to go from equator to pole to equator.  Output map points can be
centered on the north pole or south pole.

### Arguments ###

* `l`: longitude, scalar or vector, in degrees
* `b`: latitude, same number of elements as right ascension, in degrees
* `southpole` (optional boolean keyword): keyword to indicate that the plot is
  to be centered on the south pole instead of the north pole.  Default is
  `false`.

### Output ###

The 2-tuple \$(x, y)\$:

* \$x\$ coordinate, same number of elements as right ascension, normalized to be
  in the range \$[-90, 90]\$.
* \$y\$ coordinate, same number of elements as declination, normalized to be
  in the range \$[-90, 90]\$.

### Example ###

``` julia
eqpole(100, 35, southpole=true)
# => (-111.18287262822456,-19.604540237028665)
eqpole(80, 19)
# => (72.78853915267848,12.83458333897169)
```

### Notes ###

Code of this function is based on IDL Astronomy User's Library.
"""
eqpole(l::Real, b::Real; southpole::Bool=false) =
    _eqpole(promote(float(l), float(b))..., southpole)

function eqpole{L<:Real, B<:Real}(l::AbstractArray{L},
                                  b::AbstractArray{B};
                                  southpole::Bool=false)
    @assert length(l) == length(b)
    typel = typeof(float(one(L)))
    x = similar(l, typel)
    y = similar(l, typel)
    for i in eachindex(l)
        x[i], y[i] = eqpole(l[i], b[i], southpole=southpole)
    end
    return x, y
end


text:
## sample random effects given latent class membership, observed data and parameters
function sample_rand(g::Int, data::Data, param::Param)
    @unpack beta, mu, S, L, sigma, b = param
    @unpack Y, times, n, Xl, Z = data

    N, q = size(b)
    z = Matrix{Float64}(undef, N, q)
    invSigma = inv_cov_SL(S[g,:], L[g,:,:])

    tt = 1:n[1]
    Zi = Z[1,tt,:]
    ll_mean = Zi'*(Y[1,tt] - Xl[1,tt,:]*beta[g,:])
    ll_inv_cov = Zi'*Zi / sigma^2
    post_cov = inv(Symmetric(invSigma + ll_inv_cov))
    post_mean = post_cov*(invSigma*mu[g,:] + ll_mean / sigma^2)

    rand(MvNormal(post_mean, post_cov))
end

f1(t, eta=-1.5) = (1 .+ t).^eta .- 1
f2(t, nu=0.0) = t.^(nu+1) ./ (t.+1).^nu

## sample longitudinal trjactory given latent class membership, random effects, covariates and parameters
function sample_lmm(time::AbstractVector{T}, g::Int, b::Vector{T}, data::Data, param::Param) where T <: Real
    @unpack beta, sigma = param
    @unpack Xl, n = data

    X = Xl[1,1,:]
    y = zeros(T, length(time))

    for j in 1:length(time)
        Z = [1, f1(time[j]), f2(time[j])]
        Xl = vcat(X[1:3], X[4:6]*f1(time[j]), X[7:9]*f2(time[j]))
        mean_y = dot(Xl, beta[g,:]) + dot(Z, b)
        y[j] = mean_y ##rand(Normal(mean_y, sigma))
    end
    return y
end

surv_g(time, shape, scale, delta, Xe) = exp.(-exp(dot(Xe, delta))*weibull_cumhaz(time, shape, scale))

## prediction of conditional survival function and longitudinal trajectory conditional on new_data
## new_data is interpreted as data from a new subject if out_of_sample=true or as additional data from an existing
## subject in data if out_of_sample=false (for within sample prediction new data can be same as old data)
## for new subjects latent class indicator and random effects are sampled
## for existing subjects the existing posterior samples of latent class indicator and random effects are used
function predict(time::AbstractVector{T}, new_data::Data, x::Tf; out_of_sample=true) where {T <: Real, Tf <: Fit}
    param = get_samples(x)

    if length(new_data.T) > 1
        @warn "Only using first subject in new_data"
    end

    t0 = new_data.T[1]
    @assert minimum(time) >= t0

    obs_t = new_data.times[1, 1:new_data.n[1]]
    obs_y = new_data.Y[1,1:new_data.n[1]]

    if t0 > maximum(obs_t)
        time = vcat(t0, sort(time))
    end

    M = length(param)
    b = Matrix{T}(undef, M, size(param[1].b,2))
    y = Matrix{T}(undef, M, length(t))
    S = Matrix{T}(undef, M, length(t))
    g = Vector{Int}(undef, M)

    for m in 1:M
        theta = param[m]

        if out_of_sample
            g[m] = sample_cind(theta, new_data)[1] ##sample_class(new_data, theta)
            b[m,:] = sample_rand(g[m], new_data, theta)
        else
            g[m] = theta.cind[1]
            b[m,:] = theta.b[1,:]
        end

        y[m,:] = sample_lmm(time, g[m], b[m,:], new_data, theta)
        S0 = surv_g(new_data.T[1], theta.bl_shape[g[m]], theta.bl_scale[g[m]], theta.delta[g[m],:], new_data.Xe[1,:])
        S[m,:] = min.(surv_g(time, theta.bl_shape[g[m]], theta.bl_scale[g[m]], theta.delta[g[m],:], new_data.Xe[1,:]) ./ S0, 1.0)
    end

    summary_predict(obs_t, obs_y, time, y, S)
end

function summary_predict(obs_t, obs_y, t, y, S)
    S0 = repeat([1.0], inner=length(obs_t))

    S_ci_lo = vcat(S0, [quantile(S[:,i], 0.025) for i in 1:length(t)])
    S_ci_hi = vcat(S0, [quantile(S[:,i], 0.975) for i in 1:length(t)])

    y_ci_lo = vcat(obs_y, [quantile(y[:,i], 0.025) for i in 1:length(t)])
    y_ci_hi = vcat(obs_y, [quantile(y[:,i], 0.975) for i in 1:length(t)])

    DataFrame(time=vcat(obs_t, t), type=vcat(repeat(["observed"], inner=length(obs_t)), repeat(["predicted"], inner=length(t))),
              y=vcat(obs_y, mean(y, dims=1)[1,:]), sd_y=vcat(zeros(Float64, length(obs_t)), sqrt.(var(y, dims=1))[1,:]), ci_lo_y=y_ci_lo, ci_hi_y=y_ci_hi,
              S=vcat(S0, mean(S, dims=1)[1,:]), sd_S=vcat(zeros(Float64, length(obs_t)), sqrt.(var(S, dims=1)[1,:])), ci_lo_S=S_ci_lo, ci_hi_S=S_ci_hi)
end

function plot_predict(df::DataFrame; use_R=true)
    if use_R
        R"library(ggplot2)"
        @rput df
        R"ggplot(df, aes(x=df$time, y=df$y, ymin=df$ci_lo_y, ymax=df$ci_hi_y)) + geom_line() + geom_ribbon()"
        R"ggplot(df, aes(x=df$time, y=df$S, ymin=df$ci_lo_S, ymax=df$ci_hi_S)) + geom_line() + geom_ribbon()"
    else
        df_obs = df[df.type .== "observed",:]
        df_pred = df[df.type .== "predicted",:]

        p1 = plot(layer(x=df_pred.time, y=df_pred.y, ymin=df_pred.ci_lo_y, ymax=df_pred.ci_hi_y, Geom.line, Geom.ribbon),
                  layer(x=df_obs.time, y=df_obs.y, Geom.point), Guide.xlabel("Time"), Guide.ylabel("log(PSA)"),
                  Coord.cartesian(xmin=0, xmax=15, ymin=0, ymax=10))
        p2 = plot(x=df_pred.time, y=df_pred.S, ymin=df_pred.ci_lo_S, ymax=df_pred.ci_hi_S, Geom.line, Geom.ribbon, Guide.xlabel("Time"), Guide.ylabel("Survival Probability"), Coord.cartesian(xmin=0, xmax=15, ymin=0, ymax=1))

        hstack(p1, p2)
    end
end


text:
function dump_init(dgngrp, currtime, st::AbstractStateType)
    FT = eltype(Settings.Q)
    bl = Settings.dg.balance_law
    mpicomm = Settings.mpicomm
    mpirank = MPI.Comm_rank(mpicomm)

    if mpirank == 0
        # get dimensions for the interpolated grid
        dims = dimensions(dgngrp.interpol)

        # set up the variables we're going to be writing
        vars = OrderedDict()
        statenames = flattenednames(vars_state(bl, st, FT))
        for varname in statenames
            vars[varname] = (tuple(collect(keys(dims))...), FT, Dict())
        end

        dprefix = @sprintf(
            "%s_%s-%s",
            dgngrp.out_prefix,
            dgngrp.name,
            Settings.starttime,
        )
        dfilename = joinpath(Settings.output_dir, dprefix)
        init_data(dgngrp.writer, dfilename, dims, vars)
    end

    return nothing
end

max_stars_repo_name: mwarusz/CLIMA

max_issues_repo_name: mwarusz/CLIMA

max_forks_repo_name: mwarusz/CLIMA


text:
# Algorith integration rules
#
# References:
#
#  - J. Berntsen, T.O. Espelid and A. Genz. "An adaptive algorithm for the
#    approximate calculation of multiple integrals". ACM Trans. Math. Soft.,
#    17, 437-451, 1991.
#
#  - J. Berntsen and T.O. Espelid. "On the construction of higher degree
#    three-dimensional embedded integrations rules". SIAM J. Numer. Anal., 25,
#    222-234, 1988.
#
#  - A. C. Genz and A. A. Malik. "An imbedded family of fully symmetric
#    numerical integration rules". SIAM J. Numer. Anal., 20, 580-588, 1983.

function points_per_generator(n)
    p₁ = 2n
    p₂ = p₁ * (n - 1)
    p₃ = 2p₂
    p₅ = p₃ * (n - 2)
    p₄ = p₅ ÷ 3
    p₆ = 1 << n
    return (1, p₁, p₂, p₃, p₄, p₅, p₆)
end

### 11th degree rule
#
# (K₀, K₁, K₃, K₄, K₅, K₆) = (1, 5, 2, 0, 3, 2, 0); (n = 3)

# function evaluations per region
# sum(map(*, (K₀, K₁, K₃, K₄, K₅, K₆), points_per_generator(3)))
const L₁₁ = 127

const α₁ = 0.19
const α₂ = 0.50
const α₃ = 0.75
const α₄ = 0.80
const α₅ = 0.995

const β₁ = 0.9987344998351400
const β₂ = 0.7793703685672423

const ɛ₁ = 0.9999698993088767
const ɛ₂ = 0.7902637224771788
const ɛ₃ = 0.4403396687650737

const ζ₁ = 0.4378478459006862
const η₁ = 0.9549373822794593

const ζ₂ = 0.9661093133630747
const η₂ = 0.4577105877763134

const α = (α₁, α₂, α₃, α₄, α₅)
const β = (β₁, β₂)
const ɛ = (ɛ₁, ɛ₂, ɛ₃)
const ζ = (ζ₁, ζ₂)
const η = (η₁, η₂)

# weight for (0, 0, ..., 0)
const p₀ = (  0.007923078151097301  ,  1.7150062482246855    ,
              1.9360149789495251    ,  0.517082819560578     ,
              2.0544045038185197    )
# weights for α
const a₁ = (  0.06797177392788224   , -0.3755893815889209    ,
             -0.3673449403754268    ,  0.01445269144914044   ,
              0.013777599884901202  )
const a₂ = (  0.0010869865388070375 ,  0.1488632145140549    ,
              0.029297786578981855  , -0.3601489663995932    ,
             -0.576806291790441     )
const a₃ = (  0.18386336622127977   , -0.2497046640620823    ,
             -0.1151883520260315    ,  0.3628307003418485    ,
              0.03726835047700328   )
const a₄ = (  0.0336211977782927    ,  0.1792501419135204    ,
              0.050866582208722165  ,  0.007148802650872727  ,
              0.0068148789397772195 )
const a₅ = (  0.010137511233340512  ,  0.0034461267589738897 ,
              0.04453911087786472   , -0.09222852896022965   ,
              0.0572316973385185    )
# weights for β
const b₁ = (  0.0016876486839852163 , -0.005140483185555828  ,
             -0.02287828257125902   ,  0.01719339732471725   ,
             -0.044930187438112855  )
const b₂ = (  0.13464685645128072   ,  0.006536017839876424  ,
              0.029089262163458317  , -0.102141653746035     ,
              0.027292365738663484  )
# weights for ɛ
const d₁ = (  0.0017501458846003852 , -0.00065134549392297   ,
             -0.002898884350669207  , -0.007504397861080488  ,
              0.0003547473950556992 )
const d₂ = (  0.07752336383837458   , -0.006304672433547204  ,
             -0.028059634133074954  ,  0.01648362537726711   ,
              0.01571366799739551   )
const d₃ = (  0.2461864902770248    ,  0.012669593997882627  ,
              0.05638741361145884   ,  0.05234610158469334   ,
              0.04990099219278569   )
# weights for ζ and η
const e₁ = (  0.06797944868483038   , -0.005454241018647931  ,
             -0.02427469611942451   ,  0.014454323316130661  ,
              0.013779155526667704  )
const e₂ = (  0.014199628233007138  ,  0.004826995274768427  ,
              0.021483070341828822  ,  0.003019236275367777  ,
              0.0028782064230998723 )

const wa = (a₁, a₂, a₃, a₄, a₅)
const wb = (b₁, b₂)
const wd = (d₁, d₂, d₃)
const we = (e₁, e₂)

# unique(-(w[i+1] / w[i]) for w in (p₀, wa..., wb..., wd..., we...), i = 2:4)
const μ₁ = ( -1.1288675950618956    ,  -0.9780493229637744   ,
             -0.19681011641875891   ,  -0.46129836003941455  ,
             -0.28377429253731273   , -12.924397154539543    ,
             -4.45060935819115      )
const μ₂ = ( -0.2670861667822143    ,   0.0393436518667435   ,
             12.29270222952484      ,   3.1498905397991295   ,
             -0.14054025925191632   ,   2.0707312548995196   ,
              0.7515160839177919    ,   3.5113181342338917   ,
             -2.5887192979423617    ,   0.5874497614292568   ,
             -0.9283295372507698    ,   0.5954481673022622   )
const μ₃ = ( -3.9730666463921054    ,  -0.9532895608672675   ,
             -1.6015769739859858    ,  -0.10271553769262119  ,
              0.6205422333386417    ,   2.613223354846873    ,
              0.26720113428477693   ,   0.047271933288012324 )

# norms for each μ in μᵢ
#
# let ws = ((p₀,), wa, wb, wd, we)
#     ns = points_per_generator(3)[[1, 2, 3, 5, 6]]
#     μs = (μ₁, μ₂, μ₃)
#     s = Float64[]
#     for i in 1:length(μs), μ in μs[i]
#         sum = 0.0
#         for j in 1:length(n), w in ws[j]
#             sum += n[j] * abs(μ[i] * w[i+1] + w[i+2])
#         end
#         push!(s, 1/sum)
#     end
#     (s...)
# end
const s₁ = (  0.19485823831443944   ,   0.21677927747604941  ,
              0.155629450215002     ,   0.19365063964604642  ,
              0.16991906305920448   ,   0.010482710226773457 ,
              0.035771397525237954  )
const s₂ = (  0.11602998413860789   ,   0.12595256004163577  ,
              0.010803324347691502  ,   0.043690435761900195 ,
              0.12171190392783988   ,   0.0592433918191275   ,
              0.09766409013104478   ,   0.03936956866824965  ,
              0.037753982678781134  ,   0.10512211588403095  ,
              0.07529459306375538   ,   0.10482954809408154  )
const s₃ = (  0.038152895079104944  ,   0.1253631725603185   ,
              0.10318930777035908   ,   0.13160544230723548  ,
              0.08802194500708062   ,   0.03776436580400292  ,
              0.11024221575948996   ,   0.1221778013467861   )

# heuristic constants
const c = (4, 4, 0.5, 3, 0.5, 0.25)

### 9th degree rule
#
# (K₀, K₁, K₃, K₄, K₅, K₆) = (1, 4, 2, 1, 0, 0, 0); (n = 2)
# (K₀, K₁, K₃, K₄, K₅, K₆) = (1, 4, 1, 1, 1, 0, 1); (n ≥ 3)

# δ² = big"0.4707"
#const δ  = parse(BigFloat, "0.686075797561756291400285229032935720681" *
#                             "3963449254746025939952564970604511034176")

# λ₁² = 4 / (15 - 5 / δ²)
#const λ₁ = parse(BigFloat, "0.955907315804539012385720839434366012946" *
#                             "392550072657165493502911217854564875694")

# s = p = λ₁²
# d = (1 - λ₁² / δ²) / 27
# λ₂² = (1 // 7 - s / 5 - d) / (1 // 5 - s / 3 - d / δ²)
#const λ₂ = parse(BigFloat, "0.406057174738239735599606880474655826651" *
#                             "5592938007285994340533769219899878305832")

# s += λ₂²
# p *= λ₂²
# d *= (1 - λ₂² / δ²) / 3
# λ₃² = (1 // 9 - s / 7 + p / 5 - d) / (1 // 7 - s / 5 + p / 3 - d / δ²)
#const λ₃ = parse(BigFloat, "0.895254709252356257641540961306199707857" *
#                             "0302757591422795174722770397905378692193")


text:
using NonLinearReactionAdvectionDiffusionWithFrontData;

@testset "Дельта-функция                                   " begin
    N = 20;
    X = [i/N for i in 0:N];

    d = [ delta(x, X, 0.5) for x in X ]

    @testset "формула прямоугольников" begin
        I = 0.0;
        for n in 1:N
            δx = X[n+1] - X[n]
            I += d[n] * δx
        end
        @test isapprox(I, 1.0)
    end

    @testset "формула трапеций" begin
        I = 0.0;
        for n in 1:N
            δx = X[n+1] - X[n]
            I += ((d[n] + d[n+1]) / 2.0) * δx
        end
        @test isapprox(I, 1.0)
    end

    # Возвращает только один ненулевой элемент
    @test length( filter( x -> x != 0, d) ) == 1

    @testset "Проверка области определения дельта-функции" begin
        @test_throws DomainError delta(-2, X, 0.5)
        @test_throws DomainError delta(2, X, 0.5)
        @test_throws DomainError delta(0.5, X, -2)
        @test_throws DomainError delta(0.5, X, 2)
        @test_throws DomainError delta(0.5, X, 1.0)
    end

end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule VIGRA_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("VIGRA")
JLLWrappers.@generate_main_file("VIGRA", UUID("75c9d399-9f49-5d0c-80c8-ccece7b7b166"))
end  # module VIGRA_jll


text:
"""
    slider(id, label, min, max; default = Int(round((max + min) / 2))

Input element that allows users to select among a set of integers using a slider
"""
function slider(id::AbstractString, label::AbstractString, min::Integer, max::Integer; default::Integer = Int(round((max + min) / 2)))
    UIElement("""
    <v-slider
        v-model = "$id"
        min="$min"
        max="$max"
        thumb-label
        label="$label">
        </v-slider>"""),
    UIModel(id, "$default")
end

"""
    text_input(id, label; default = "")

Free text field for string input by users

Note that a text input always returns a String to the server. If you need users
to be able to freely enter a number, use `number_input`.
"""
function text_input(id::AbstractString, label::AbstractString; default::AbstractString = "")
    UIElement("""
    <v-text-field
        v-model = "$id"
        label="$label"
        filled>
    </v-text-field>
    """),
    UIModel(id, "\"$default\"")
end

"""
    number_input(id, label; default = 0)

Form input that only accepts numbers

Note that `number_input` always returns a float to the server. If you need a string
use `text_input`. If you need integers, you can use a `slider` or round the result.
"""
function number_input(id::AbstractString, label::AbstractString; default::Integer = 0)
    UIElement("""
    <v-text-field
        v-model.number = "$id"
        label="$label"
        type = "number"
        filled>
    </v-text-field>
    """),
    UIModel(id, "$default")
end

"""
    button(id, label; color, size)

Add a button to a UI. Buttons return `true` to the server when clicked, and false otherwise.

size can be one of `x-small`, `small`, `normal`, `large`, `x-large`
color can be any valid color (see docs on colors for a full list) - e.g. `primary`, `error`, `teal` etc
"""
function button(id::AbstractString, label::AbstractString; color::AbstractString = "primary", size::AbstractString = "normal")
    UIElement("""
    <v-btn $size color="$color" v-model="$id" @click="fetch_update_$id(true)">$label</v-btn>
    """),
    UIModel(id, "false")
end

"""
    button_group(id, buttons...)

Group together a set of buttons to provide `radio`-like exclusive selections.

The list of `buttons` given as input should only contain `button`s. The number (0-based index)
of the selected button is returned to the server. A selection is mandatory and exclusive.

Component buttons can still have their own `id`s and will still tell the server when they
are pushed.
"""
function button_group(id::AbstractString, buttons...)
    UIElement("""
    <v-btn-toggle
         v-model="$id"
         rounded
         mandatory
       >"""),
    buttons...,
    UIElement("</v-btn-toggle>"),
    UIModel(id, "null")
end

"""
    floating_action_button(id, label; location = "bottom right", color = "red")

Add a floating action button to your UI at `location`
"""
function floating_action_button(id::AbstractString, label::AbstractString; location::AbstractString = "bottom right", color::AbstractString = "red", size::AbstractString = "normal")
    UIElement("""<v-btn absolute fab $location color="$color" @click="fetch_update_$id(true)">$label</v-btn>"""),
    UIModel(id, "false")
end

"""
    date_picker(id::AbstractString; color::AbstractString = "primary")

Add a date picker to your UI. Date is returned as a String in ISO 8601 form (YYYY-MM-DD)
"""
function date_picker(id::AbstractString; color::AbstractString = "primary")
    UIElement("""<v-date-picker v-model="$id" color="$color" header-color="$color"></v-date-picker>"""),
    UIModel(id, "new Date().toISOString().substr(0, 10)")
end

"""
    time_picker(id; color = "primary", default = nothing)

Add a clock time picker to your UI. Time is returned as a String, in HH:MM 24 hour format.
Returns `nothing` if user has not yet picked a time.
"""
function time_picker(id::AbstractString; color::AbstractString = "primary", default::Union{AbstractString, Nothing} = nothing)
    UIElement("""<v-time-picker v-model="$id" color = "$color"></v-time-picker>"""),
    UIModel(id, isnothing(default) ? "null" : "\"$default\"")
end

"""
    function selector(id, label, items; multiple = false, autocomplete = false)

Create a selection box in your UI for users to choose among options. Allow multiple selections
with `multiple`. Let users type in the box to filter options by setting `autocomplete` to `true`.
`items` can either be a string for a static (javascript) array of options (i.e. "['a first option',
'the second option']") or an id that corresponds to a function in your Server that returns an
array of options (potentially dynamically).

Returns either an Array (if no or multiple elements are selected) or the type of the individual
element if only one is selected. The type that is returned depends on the types of the elements
in the select list.
"""
function selector(id::AbstractString, label::AbstractString, items::AbstractString; multiple::Bool = false, autocomplete::Bool = false)
    if occursin("[", items)
        items_model = ()
    else
        items_model = UIModel(items, "[]"), UIWatch(items, "this.$id = []")
    end
    tag = autocomplete ? "v-autocomplete" : "v-select"
    UIElement("""
    <$tag
            v-model="$id"
            :items="$items"
            attach
            label="$label"
            filled
            $(multiple ? "multiple\nchips" : "")></$tag>"""),
    UIModel(id, "[]"),
    items_model
end

"""
    checkbox(id, label; default = false)

Add a checkbox that returns true or false based on whether it is ticked.
"""
function checkbox(id, label; default = false)
    UIElement("""
    <v-checkbox
      v-model="$id"
      label="$label"></v-checkbox>"""),
    UIModel(id, "$default")
end

quote_wrap_if_string(x::AbstractString) = "'$x'"
quote_wrap_if_string(x) = x

"""
    radio(id, values, labels)

Give users a set of options to choose from (exclusively).

Unlike a select, all radio options are displayed (not a drop down list), and cannot be
dynamic.
"""
function radio(id::AbstractString, values::Array, labels::Array{<:AbstractString, 1})
    length(values) != length(labels) && error("`values` and `labels` must have the same length")

    buttons = map((value, label) -> UIElement("""<v-radio label="$label" :value="$(quote_wrap_if_string(value))"></v-radio>"""), values, labels)

    UIElement("""<v-radio-group v-model="$id" :mandatory="false">"""),
    buttons...,
    UIElement("</v-radio-group>"),
    UIModel(id, "null")
end

"""
    switch(id, label; default = false)

Add a switch that returns true or false based on whether it is turned on.
"""
function switch(id, label; default = false)
    UIElement("""
    <v-switch
      v-model="$id"
      label="$label"
    ></v-switch>
    """),
  UIModel(id, "$default")
end

"""
    function list(id, items...; title = nothing)

Create a list of items for users to choose from. `items...` should be populated with
`list_item`s.

`title` is optional; if specified it provides a grey title at the top of the list.

The selected list item is returned to the server as a zero-based index.
"""
function list(id, items...; title = nothing)
    if title != nothing
        title = UIElement("<v-subheader>REPORTS</v-subheader>")
    else
        title = ()
    end
    UIElement("""<v-list><v-list-item-group v-model = "$id" color="primary">"""),
    title,
    items...,
    UIElement("</v-list-item-group></v-list>"),
    UIModel(id, "null")
end

"""
    list_item(title, subtitle = nothing)

Should only be used as part of a `list`

`subtitle` is optional. It provides more detail in a smaller font beneath the main item.
"""
function list_item(title, subtitle = nothing)
    if subtitle != nothing
        UIElement("
        <v-list-item two-line>
            <v-list-item-content>
            <v-list-item-title>")
        title,
        UIElement("</v-list-item-title>
            <v-list-item-subtitle>"),
        subtitle,
        UIElement("</v-list-item-subtitle>
            </v-list-item-content>
            </v-list-item>")
    else
        UIElement("
        <v-list-item>
            <v-list-item-content>
            <v-list-item-title>"),
        title,
        UIElement("</v-list-item-title>
            </v-list-item-content>
            </v-list-item>")
    end
end


text:
using BenchmarkTools
using GraphAlgorithms
using LightGraphs,SimpleWeightedGraphs

DIGRAPHS = Dict{String,DiGraph}(
    "complete20" => complete_digraph(20),
    "path40" => path_digraph(40)
)
GRAPHS = Dict{String,Graph}(
    "petersen" => smallgraph(:petersen),
    "tutte" => smallgraph(:tutte),
    "complete20" => complete_graph(20)
)

suite = BenchmarkGroup()
include("operators.jl")
include("sets.jl")
include("diversity.jl")
include("fractional_numbers.jl")

tune!(suite);
results = run(suite, verbose=true, seconds=60, samples=10000)


text:
using ChainRulesTestUtils, FiniteDifferences, Zygote, Adapt, CUDA
CUDA.allowscalar(false)

# global GRAPH_T = :coo
# global TEST_GPU = true

const rule_config = Zygote.ZygoteRuleConfig()

# Using this until https://github.com/JuliaDiff/FiniteDifferences.jl/issues/188 is fixed
function FiniteDifferences.to_vec(x::Integer)
    Integer_from_vec(v) = x
    return Int[x], Integer_from_vec
end

# Test that forward pass on cpu and gpu are the same.
# Tests also gradient on cpu and gpu comparing with
# finite difference methods.
# Test gradients with respects to layer weights and to input.
# If `g` has edge features, it is assumed that the layer can be
# use them in the forward pass as `l(g, x, e)`.
# Test also gradient with repspect to `e`.
function test_layer(l, g::GNNGraph; atol = 1e-6, rtol = 1e-5,
                                 exclude_grad_fields = [],
                                 broken_grad_fields =[],
                                 verbose = false,
                                 test_gpu = TEST_GPU,
                                 outsize = nothing,
                                 outtype = :node,
                                )

    # TODO these give errors, probably some bugs in ChainRulesTestUtils
    # test_rrule(rule_config, x -> l(g, x), x; rrule_f=rrule_via_ad, check_inferred=false)
    # test_rrule(rule_config, l -> l(g, x), l; rrule_f=rrule_via_ad, check_inferred=false)

    isnothing(node_features(g)) && error("Plese add node data to the input graph")
    fdm = central_fdm(5, 1)

    x = node_features(g)
    e = edge_features(g)
    use_edge_feat = !isnothing(e)

    x64, e64, l64, g64 = to64.([x, e, l, g]) # needed for accurate FiniteDifferences' grad
    xgpu, egpu, lgpu, ggpu = gpu.([x, e, l, g])

    f(l, g::GNNGraph) = l(g)
    f(l, g::GNNGraph, x, e) = use_edge_feat ? l(g, x, e) : l(g, x)

    loss(l, g::GNNGraph) = if outtype == :node
                                sum(node_features(f(l, g)))
                            elseif outtype == :edge
                                sum(edge_features(f(l, g)))
                            elseif outtype == :graph
                                sum(graph_features(f(l, g)))
                            elseif outtype == :node_edge
                                gnew = f(l, g)
                                sum(node_features(gnew)) + sum(edge_features(gnew))
                            end

    function loss(l, g::GNNGraph, x, e)
        y = f(l, g, x, e)
        if outtype == :node_edge
            return sum(y[1]) + sum(y[2])
        else
            return sum(y)
        end
    end

    # TEST OUTPUT
    y = f(l, g, x, e)
    if outtype == :node_edge
        @assert y isa Tuple
        @test eltype(y[1]) == eltype(x)
        @test eltype(y[2]) == eltype(e)
        @test all(isfinite, y[1])
        @test all(isfinite, y[2])
        if !isnothing(outsize)
            @test size(y[1]) == outsize[1]
            @test size(y[2]) == outsize[2]
        end
    else
        @test eltype(y) == eltype(x)
        @test all(isfinite, y)
        if !isnothing(outsize)
            @test size(y) == outsize
        end
    end

    # test same output on different graph formats
    gcoo = GNNGraph(g, graph_type=:coo)
    ycoo = f(l, gcoo, x, e)
    if outtype == :node_edge
        @test ycoo[1] ≈ y[1]
        @test ycoo[2] ≈ y[2]
    else
        @test ycoo ≈ y
    end

    g′ = f(l, g)
    if outtype == :node
        @test g′.ndata.x ≈ y
    elseif outtype == :edge
        @test g′.edata.e ≈ y
    elseif outtype == :graph
        @test g′.gdata.u ≈ y
    elseif outtype == :node_edge
        @test g′.ndata.x ≈ y[1]
        @test g′.edata.e ≈ y[2]
    else
        @error "wrong outtype $outtype"
    end
    if test_gpu
        ygpu = f(lgpu, ggpu, xgpu, egpu)
        if outtype == :node_edge
            @test ygpu[1] isa CuArray
            @test eltype(ygpu[1]) == eltype(xgpu)
            @test Array(ygpu[1]) ≈ y[1]
            @test ygpu[2] isa CuArray
            @test eltype(ygpu[2]) == eltype(xgpu)
            @test Array(ygpu[2]) ≈ y[2]
        else
            @test ygpu isa CuArray
            @test eltype(ygpu) == eltype(xgpu)
            @test Array(ygpu) ≈ y
        end
    end

    # TEST x INPUT GRADIENT
    x̄  = gradient(x -> loss(l, g, x, e), x)[1]
    x̄_fd = FiniteDifferences.grad(fdm, x64 -> loss(l64, g64, x64, e64), x64)[1]
    @test eltype(x̄) == eltype(x)
    @test x̄ ≈ x̄_fd    atol=atol rtol=rtol

    if test_gpu
        x̄gpu  = gradient(xgpu -> loss(lgpu, ggpu, xgpu, egpu), xgpu)[1]
        @test x̄gpu isa CuArray
        @test eltype(x̄gpu) == eltype(x)
        @test Array(x̄gpu) ≈ x̄   atol=atol rtol=rtol
    end

    # TEST e INPUT GRADIENT
    if e !== nothing
        verbose && println("Test e gradient cpu")
        ē  = gradient(e -> loss(l, g, x, e), e)[1]
        ē_fd = FiniteDifferences.grad(fdm, e64 -> loss(l64, g64, x64, e64), e64)[1]
        @test eltype(ē) == eltype(e)
        @test ē ≈ ē_fd    atol=atol rtol=rtol

        if test_gpu
            verbose && println("Test e gradient gpu")
            ēgpu  = gradient(egpu -> loss(lgpu, ggpu, xgpu, egpu), egpu)[1]
            @test ēgpu isa CuArray
            @test eltype(ēgpu) == eltype(ē)
            @test Array(ēgpu) ≈ ē   atol=atol rtol=rtol
        end
    end

    # TEST LAYER GRADIENT - l(g, x, e)
    l̄ = gradient(l -> loss(l, g, x, e), l)[1]
    l̄_fd = FiniteDifferences.grad(fdm, l64 -> loss(l64, g64, x64, e64), l64)[1]
    test_approx_structs(l, l̄, l̄_fd; atol, rtol, broken_grad_fields, exclude_grad_fields, verbose)

    if test_gpu
        l̄gpu = gradient(lgpu -> loss(lgpu, ggpu, xgpu, egpu), lgpu)[1]
        test_approx_structs(lgpu, l̄gpu, l̄; atol, rtol, broken_grad_fields, exclude_grad_fields, verbose)
    end

    # TEST LAYER GRADIENT - l(g)
    l̄ = gradient(l -> loss(l, g), l)[1]
    test_approx_structs(l, l̄, l̄_fd; atol, rtol, broken_grad_fields, exclude_grad_fields, verbose)

    return true
end

function test_approx_structs(l, l̄, l̄2; atol=1e-5, rtol=1e-5,
            broken_grad_fields=[],
            exclude_grad_fields=[],
            verbose=false)

    l̄ = l̄ isa Base.RefValue ? l̄[] : l̄           # Zygote wraps gradient of mutables in RefValue
    l̄2 = l̄2 isa Base.RefValue ? l̄2[] : l̄2           # Zygote wraps gradient of mutables in RefValue

    for f in fieldnames(typeof(l))
        f ∈ exclude_grad_fields && continue
        f̄, f̄2 = getfield(l̄, f), getfield(l̄2, f)
        x = getfield(l, f)
        if verbose
            println()
            @show f x f̄ f̄2
        end
        if isnothing(f̄)
            verbose && println("A")
            @test !(f̄2 isa AbstractArray) || isapprox(f̄2, fill!(similar(f̄2), 0); atol=atol, rtol=rtol)
        elseif f̄ isa Union{AbstractArray, Number}
            verbose && println("B")
            @test eltype(f̄) == eltype(x)
            if x isa CuArray
                @test f̄ isa CuArray
                f̄ = Array(f̄)
            end
            if f ∈ broken_grad_fields
                @test_broken f̄ ≈ f̄2   atol=atol rtol=rtol
            else
                @test f̄ ≈ f̄2   atol=atol rtol=rtol
            end
        else
            verbose && println("C")
            test_approx_structs(x, f̄, f̄2; atol, rtol, exclude_grad_fields, broken_grad_fields, verbose)
        end
    end
    return true
end

"""
    to32(m)

Convert the `eltype` of model's parameters to `Float32` or `Int32`.
"""
function to32(m)
    f(x::AbstractArray) = eltype(x) <: Integer ? adapt(Int32, x) : adapt(Float32, x)
    f(x::Number) = typeof(x) <: Integer ? adapt(Int32, x) : adapt(Float32, x)
    f(x) = adapt(Float32, x)
    return fmap(f, m)
end

"""
    to64(m)

Convert the `eltype` of model's parameters to `Float64` or `Int64`.
"""
function to64(m)
    f(x::AbstractArray) = eltype(x) <: Integer ? adapt(Int64, x) : adapt(Float64, x)
    f(x::Number) = typeof(x) <: Integer ? adapt(Int64, x) : adapt(Float64, x)
    f(x) = adapt(Float64, x)
    return fmap(f, m)
end


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

ENV["LINES"] = 20

using IJulia
notebook(dir=".")


text:
module readwrite

export readschema, writeschema

using TOML
using ..types

"Returns either a Schema or a TableSchema, depending on the contents of the config file."
function readschema(filename::String)
    d = TOML.parsefile(filename)
    haskey(d, "columns") && return TableSchema(d)  # Config is for a TableSchema
    Schema(d)  # Config is for a Schema
end

writeschema(outfile::String, schema::Schema)           = toml_to_file(outfile, schema_to_dict(schema))
writeschema(outfile::String, tableschema::TableSchema) = toml_to_file(outfile, tableschema_to_dict(tableschema))

function toml_to_file(outfile::String, d)
    io = open(outfile, "w")
    TOML.print(io, d)
    close(io)
end

################################################################################
# Utils

function schema_to_dict(schema::Schema)
    result = Dict{String, Any}()
    result["name"]        = String(schema.name)
    result["description"] = schema.description
    result["tables"]      = Dict(String(tablename) => tableschema_to_dict(tableschema) for (tablename, tableschema) in schema.tables)
    result
end

function tableschema_to_dict(tableschema::TableSchema)
    result = Dict{String, Any}()
    result["name"]        = String(tableschema.name)
    result["description"] = tableschema.description
    result["primarykey"]  = String.(tableschema.primarykey)
    columns = Dict{String, Any}[]  # colname => colschema
    for colname in tableschema.columnorder
        push!(columns, colschema_to_dict(tableschema.colname2colschema[colname]))
    end
    result["columns"] = columns
    if !isempty(tableschema.intrarow_constraints)
        result["intrarow_constraints"] = Dict(msg => func_as_supplied for (func_as_supplied, f, msg) in tableschema.intrarow_constraints)
    end
    result
end

function colschema_to_dict(colschema::ColumnSchema)
    result = Dict{String, Any}()
    result["name"]          = string(colschema.name)
    result["description"]   = colschema.description
    result["datatype"]      = string(colschema.datatype)
    result["iscategorical"] = colschema.iscategorical
    result["isrequired"]    = colschema.isrequired
    result["isunique"]      = colschema.isunique
    result["validvalues"]   = format_validvalues(colschema.validvalues, colschema.valueorder)
    if !isnothing(colschema.parser_as_supplied)
        result["parser"] = colschema.parser_as_supplied
    end
    result
end

format_validvalues(vv::DataType, valueorder)      = string(vv)
format_validvalues(vv::AbstractRange, valueorder) = string(vv)
format_validvalues(vv::Set, valueorder::Nothing)  = sort!([x for x in vv])
format_validvalues(vv::Set, valueorder::Vector)   = valueorder

end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Word2Vec"
version = v"0.1.0"

# Collection of sources required to complete build
sources = [
    GitSource("https://github.com/tmikolov/word2vec.git", "20c129af10659f7c50e86e3be406df663beff438"),
    DirectorySource("./bundled")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir
for f in ${WORKSPACE}/srcdir/patches/*.patch; do
    atomic_patch -p1 ${f}
done
cd word2vec/
export CFLAGS="-lm -pthread -O3 -Wall -funroll-loops -Wno-unused-result"
make word2vec
make word2phrase
mkdir -p $bindir
cp word2vec$exeext $bindir
cp word2phrase$exeext $bindir
install_license LICENSE
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = supported_platforms()
# Too many posix assumptions
filter!(!Sys.iswindows, platforms)

# The products that we will ensure are always built
products = [
    ExecutableProduct("word2vec", :word2vec)
    ExecutableProduct("word2phrase", :word2phrase)
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies)

max_stars_repo_name: sharanry/Yggdrasil

max_issues_repo_name: sharanry/Yggdrasil

max_forks_repo_name: sharanry/Yggdrasil


text:
# This is it: the bottom of the ocean trench, the last place the sleigh keys
# could be. Your submarine's experimental antenna still isn't boosted enough to
# detect the keys, but they must be here. All you need to do is reach the
# seafloor and find them.
#
# At least, you'd touch down on the seafloor if you could; unfortunately, it's
# completely covered by two large herds of sea cucumbers, and there isn't an
# open space large enough for your submarine.
#
# You suspect that the Elves must have done this before, because just then you
# discover the phone number of a deep-sea marine biologist on a handwritten
# note taped to the wall of the submarine's cockpit.
#
# "Sea cucumbers? Yeah, they're probably hunting for food. But don't worry,
# they're predictable critters: they move in perfectly straight lines, only
# moving forward when there's space to do so. They're actually quite polite!"
#
# You explain that you'd like to predict when you could land your submarine.
#
# "Oh that's easy, they'll eventually pile up and leave enough space for--
# wait, did you say submarine? And the only place with that many sea cucumbers
# would be at the very bottom of the Mariana--" You hang up the phone.
#
# There are two herds of sea cucumbers sharing the same region; one always
# moves east (>), while the other always moves south (v). Each location can
# contain at most one sea cucumber; the remaining locations are empty (.). The
# submarine helpfully generates a map of the situation (your puzzle input). For
# example:
#
# v...>>.vv>
# .vv>>.vv..
# >>.>v>...v
# >>v>>.>.v.
# v>v.vv.v..
# >.>>..v...
# .vv..>.>v.
# v.v..>>v.v
# ....v..v.>
#
# Every step, the sea cucumbers in the east-facing herd attempt to move forward
# one location, then the sea cucumbers in the south-facing herd attempt to move
# forward one location. When a herd moves forward, every sea cucumber in the
# herd first simultaneously considers whether there is a sea cucumber in the
# adjacent location it's facing (even another sea cucumber facing the same
# direction), and then every sea cucumber facing an empty location
# simultaneously moves into that location.
#
# So, in a situation like this:
#
# ...>>>>>...
#
# After one step, only the rightmost sea cucumber would have moved:
#
# ...>>>>.>..
#
# After the next step, two sea cucumbers move:
#
# ...>>>.>.>.
#
# During a single step, the east-facing herd moves first, then the south-facing
# herd moves. So, given this situation:
#
# ..........
# .>v....v..
# .......>..
# ..........
#
# After a single step, of the sea cucumbers on the left, only the south-facing
# sea cucumber has moved (as it wasn't out of the way in time for the
# east-facing cucumber on the left to move), but both sea cucumbers on the
# right have moved (as the east-facing sea cucumber moved out of the way of the
# south-facing sea cucumber):
#
# ..........
# .>........
# ..v....v>.
# ..........
#
# Due to strong water currents in the area, sea cucumbers that move off the
# right edge of the map appear on the left edge, and sea cucumbers that move
# off the bottom edge of the map appear on the top edge. Sea cucumbers always
# check whether their destination location is empty before moving, even if that
# destination is on the opposite side of the map:
#
# Initial state:
# ...>...
# .......
# ......>
# v.....>
# ......>
# .......
# ..vvv..
#
# After 1 step:
# ..vv>..
# .......
# >......
# v.....>
# >......
# .......
# ....v..
#
# After 2 steps:
# ....v>.
# ..vv...
# .>.....
# ......>
# v>.....
# .......
# .......
#
# After 3 steps:
# ......>
# ..v.v..
# ..>v...
# >......
# ..>....
# v......
# .......
#
# After 4 steps:
# >......
# ..v....
# ..>.v..
# .>.v...
# ...>...
# .......
# v......
#
# To find a safe place to land your submarine, the sea cucumbers need to stop
# moving. Again consider the first example:
#
# Initial state:
# v...>>.vv>
# .vv>>.vv..
# >>.>v>...v
# >>v>>.>.v.
# v>v.vv.v..
# >.>>..v...
# .vv..>.>v.
# v.v..>>v.v
# ....v..v.>
#
# After 1 step:
# ....>.>v.>
# v.v>.>v.v.
# >v>>..>v..
# >>v>v>.>.v
# .>v.v...v.
# v>>.>vvv..
# ..v...>>..
# vv...>>vv.
# >.v.v..v.v
#
# After 2 steps:
# >.v.v>>..v
# v.v.>>vv..
# >v>.>.>.v.
# >>v>v.>v>.
# .>..v....v
# .>v>>.v.v.
# v....v>v>.
# .vv..>>v..
# v>.....vv.
#
# After 3 steps:
# v>v.v>.>v.
# v...>>.v.v
# >vv>.>v>..
# >>v>v.>.v>
# ..>....v..
# .>.>v>v..v
# ..v..v>vv>
# v.v..>>v..
# .v>....v..
#
# After 4 steps:
# v>..v.>>..
# v.v.>.>.v.
# >vv.>>.v>v
# >>.>..v>.>
# ..v>v...v.
# ..>>.>vv..
# >.v.vv>v.v
# .....>>vv.
# vvv>...v..
#
# After 5 steps:
# vv>...>v>.
# v.v.v>.>v.
# >.v.>.>.>v
# >v>.>..v>>
# ..v>v.v...
# ..>.>>vvv.
# .>...v>v..
# ..v.v>>v.v
# v.v.>...v.
#
# ...
#
# After 10 steps:
# ..>..>>vv.
# v.....>>.v
# ..v.v>>>v>
# v>.>v.>>>.
# ..v>v.vv.v
# .v.>>>.v..
# v.v..>v>..
# ..v...>v.>
# .vv..v>vv.
#
# ...
#
# After 20 steps:
# v>.....>>.
# >vv>.....v
# .>v>v.vv>>
# v>>>v.>v.>
# ....vv>v..
# .v.>>>vvv.
# ..v..>>vv.
# v.v...>>.v
# ..v.....v>
#
# ...
#
# After 30 steps:
# .vv.v..>>>
# v>...v...>
# >.v>.>vv.>
# >v>.>.>v.>
# .>..v.vv..
# ..v>..>>v.
# ....v>..>v
# v.v...>vv>
# v.v...>vvv
#
# ...
#
# After 40 steps:
# >>v>v..v..
# ..>>v..vv.
# ..>>>v.>.v
# ..>>>>vvv>
# v.....>...
# v.v...>v>>
# >vv.....v>
# .>v...v.>v
# vvv.v..v.>
#
# ...
#
# After 50 steps:
# ..>>v>vv.v
# ..v.>>vv..
# v.>>v>>v..
# ..>>>>>vv.
# vvv....>vv
# ..v....>>>
# v>.......>
# .vv>....v>
# .>v.vv.v..
#
# ...
#
# After 55 steps:
# ..>>v>vv..
# ..v.>>vv..
# ..>>v>>vv.
# ..>>>>>vv.
# v......>vv
# v>v....>>v
# vvv...>..>
# >vv.....>.
# .>v.vv.v..
#
# After 56 steps:
# ..>>v>vv..
# ..v.>>vv..
# ..>>v>>vv.
# ..>>>>>vv.
# v......>vv
# v>v....>>v
# vvv....>.>
# >vv......>
# .>v.vv.v..
#
# After 57 steps:
# ..>>v>vv..
# ..v.>>vv..
# ..>>v>>vv.
# ..>>>>>vv.
# v......>vv
# v>v....>>v
# vvv.....>>
# >vv......>
# .>v.vv.v..
#
# After 58 steps:
# ..>>v>vv..
# ..v.>>vv..
# ..>>v>>vv.
# ..>>>>>vv.
# v......>vv
# v>v....>>v
# vvv.....>>
# >vv......>
# .>v.vv.v..
#
# In this example, the sea cucumbers stop moving after 58 steps.
#
# Find somewhere safe to land your submarine. What is the first step on which
# no sea cucumbers move?

# read input
cucumbers = vcat(map(l -> reshape(Vector{Char}(l), 1, :), eachline("input.txt"))...)
(height, width) = size(cucumbers)

# do eet
for step in Iterators.countfrom(1)
  nextcucumbers = copy(cucumbers)

  # move cucumbers to the right
  shifted = @view nextcucumbers[:, [2:width..., 1]]
  selected = nextcucumbers .== '>' .&& shifted .== '.'
  shifted[selected] .= '>'
  nextcucumbers[selected] .= '.'

  # move cucumbers down
  shifted = @view nextcucumbers[[2:height..., 1], :]
  selected = nextcucumbers .== 'v' .&& shifted .== '.'
  shifted[selected] .= 'v'
  nextcucumbers[selected] .= '.'

  if cucumbers == nextcucumbers
    println("Step: ", step)
    break
  end
  global cucumbers = nextcucumbers
end


text:
using VortexLattice

# wing
xle = [0.0, 0.2]
yle = [0.0, 5.0]
zle = [0.0, 1.0]
chord = [1.0, 0.6]
theta = [-1.0*pi/180, -1.0*pi/180]
phi = [0.0, 0.0]
fc = fill((xc) -> 0, 2) # camberline function for each section
ns = 12
nc = 6
spacing_s = Uniform()
spacing_c = Uniform()
mirror = false

# horizontal stabilizer
xle_h = [0.0, 0.14]
yle_h = [0.0, 1.25]
zle_h = [0.0, 0.0]
chord_h = [0.7, 0.42]
theta_h = [0.0, 0.0]
phi_h = [0.0, 0.0]
fc_h = fill((xc) -> 0, 2) #camberline function for each section
ns_h = 6
nc_h = 3
spacing_s_h = Uniform()
spacing_c_h = Uniform()
mirror_h = false

# vertical stabilizer
xle_v = [0.0, 0.14]
yle_v = [0.0, 0.0]
zle_v = [0.0, 1.0]
chord_v = [0.7, 0.42]
theta_v = [0.0, 0.0]
phi_v = [0.0, 0.0]
fc_v = fill((xc) -> 0, 2) #camberline function for each section
ns_v = 5
nc_v = 3
spacing_s_v = Uniform()
spacing_c_v = Uniform()
mirror_v = false

Sref = 8.0
cref = 0.9
bref = 10.0
rref = [1.9437, 0.0, 0.3708]#TODO Always make this the correct CoM
Vinf =  21.0967#Was originally 1
ref = Reference(Sref, cref, bref, rref, Vinf)

alpha = 1.2129*pi/180
beta = 0.0
Omega = [0.0; 0.0; 0.0]
fs = Freestream(Vinf, alpha, beta, Omega)

symmetric = [false, false, false,false,false]

# generate surface panels for wing
wgrid, wing = wing_to_surface_panels(xle, yle, zle, chord, theta, phi, ns, nc;
    mirror=mirror, fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)
wgrid2, wing2 = wing_to_surface_panels(xle, -yle, zle, chord, theta, phi, ns, nc;
        mirror=mirror, fc = fc, spacing_s=spacing_s, spacing_c=spacing_c)

# generate surface panels for horizontal tail
hgrid, htail = wing_to_surface_panels(xle_h, yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
    mirror=mirror_h, fc=fc_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
translate!(hgrid, [4.0, 0.0, 0.0])
translate!(htail, [4.0, 0.0, 0.0])
hgrid2, htail2 = wing_to_surface_panels(xle_h, -yle_h, zle_h, chord_h, theta_h, phi_h, ns_h, nc_h;
    mirror=mirror_h, fc=fc_h, spacing_s=spacing_s_h, spacing_c=spacing_c_h)
translate!(hgrid2, [4.0, 0.0, 0.0])
translate!(htail2, [4.0, 0.0, 0.0])

# generate surface panels for vertical tail
vgrid, vtail = wing_to_surface_panels(xle_v, yle_v, zle_v, chord_v, theta_v, phi_v, ns_v, nc_v;
    mirror=mirror_v, fc=fc_v, spacing_s=spacing_s_v, spacing_c=spacing_c_v)
translate!(vgrid, [4.0, 0.0, 0.0])
translate!(vtail, [4.0, 0.0, 0.0])

grids = [wgrid, hgrid,wgrid2, hgrid2, vgrid] #TODO
surfaces = [wing, htail,wing2, htail2, vtail] #TODO
surface_id = [1, 2, 3,4,5]

system = steady_analysis(surfaces, ref, fs; symmetric=symmetric, surface_id=surface_id)

CF, CM = body_forces(system; frame=Wind())

CDiff = far_field_drag(system)

CD, CY, CL = CF
Cl, Cm, Cn = CM

dCF, dCM = stability_derivatives(system)

CDa, CYa, CLa = dCF.alpha
Cla, Cma, Cna = dCM.alpha
CDb, CYb, CLb = dCF.beta
Clb, Cmb, Cnb = dCM.beta
CDp, CYp, CLp = dCF.p
Clp, Cmp, Cnp = dCM.p
CDq, CYq, CLq = dCF.q
Clq, Cmq, Cnq = dCM.q
CDr, CYr, CLr = dCF.r
Clr, Cmr, Cnr = dCM.r

#Visualize the model
#properties= get_surface_properties(system)
#write_vtk("Anim/model",surfaces,properties; symmetric)

using SixDOF

CLM = 0.0 # Mach derivative
CLdf = 0.0  # flaps derivative
CLde = 0.0  # elevator derivative
alphas = 15*pi/180
U0=21.0967
exp_Re = -0.2  # exponent in Reynolds number calc
Mcc = 0.7  # crest critcal Mach number
e = 0.80  # Oswald efficiency#TODO Adjust
CDdf = 0.0  # flaps
CDde = 0.0  # elevator
CDda = 0.0  # aileron
CDdr = 0.0  # rudder
CYda = 0.0 # Roll control (aileron) derivative
CYdr = 0.0 # Yaw control (rudder) derivative
Clda = 0.0  # Roll (aileron) control derivative
Cldr = 0.0  #Yaw (rudder) control derivative
Cmdf = 0.0
Cmde = 0.0 # Pitch control derivative
CmM = 0.0
Cnda = 0.0  # Roll (aileron) control derivative
Cndr = 0.0  # Yaw (rudder) control derivative

m = 7.25755
Ixx = 38.05
Iyy = 26.46
Izz = 62.97
Ixz = -2.133
mp = MassProp(m, Ixx, Iyy, Izz, Ixz)

Wi = [0.0, 0.0, 0.0]
Wb = [0.0, 0.0, 0.0]
rho = 1.225
asound = 300.0
g = 9.81
atm = ConstantAtmosphere(Wi, Wb, rho, asound, g)

controller = ConstantController(0.0, 0.0, 0.0, 0.0, 0.0)

sd = StabilityDeriv(CL, CLa, CLq, CLM, CLdf, CLde, alphas,
    CD, U0, exp_Re, e, Mcc, CDdf, CDde, CDda, CDdr,
    CYb, -CYp, -CYr, CYda, CYdr,
    -Clb, Clp, Clr, Clda, Cldr,
    Cm, Cma, Cmq, CmM, Cmdf, Cmde,
    -Cnb, Cnp, Cnr, Cnda, Cndr) #some are negative to account for changes in refrence frame.
nothing # hide

CT0 = 0.11221
CT1 = -0.13803
CT2 = -0.047394
CQ0 = 0.0062
CQ1 = 0.00314
CQ2 = -0.015729
D = 10*0.0254
num = 0
type = COCOUNTER
R = 0.5
Kv = 2500.0 * pi/30
i0 = 0.3
voltage = 8.0
propulsion = MotorPropBatteryDataFit(CT2, CT1, CT0, CQ2, CQ1, CQ0, D, num, type, R, Kv, i0, voltage)
nothing # hide

inertial2 = UniformGravitationalField()
nothing # hide

Vinf = U0
alpha = 1.21290*pi/180
s0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, U0, 0.0, 0.0, 0.0, 0.0, 0.0]
tspan = (0.0, 4.0)
p = mp, ref, sd, propulsion, inertial2, atm, controller

#Tyson testing his stuff
JACOBIAN=getjacobian(s0,p)
s0 = [0.0, 0.0, 0.0, 0.0, 0.0, 0.0, U0, 0.0, 0.0, 0.0, 0.0, 0.0]
p = mp, ref, sd, propulsion, inertial2, atm, controller

leftWing = Wing([0.0 0.0 0.0],[-.1 -1 -.2], .3, .1, 0.0, 3.0)
rightWing = Wing([0.0 0.0 0.0],[-.1  1 -.2], .3, .1, 0.0, 3.0)
leftTail = Wing([-0.8 0.0 0.0],[-.05 -.33 -0.0], .3, .05, 0.0, 0.0)
rightTail = Wing([-0.8 0.0 0.0],[-.05  .33 -0.0], .3, .05, 0.0, 0.0)
rudder = Wing([-0.8 0.0 0.0],[-.05  0.0 -0.33], .3, .05, 0.0, 0.0)
planeModel = (leftWing,rightWing,leftTail,rightTail,rudder)
animvtk(JACOBIAN, planeModel,"Test")

using LinearAlgebra
Eigenvals, Eigenvecs =eigen(JACOBIAN)

longitudinal = [JACOBIAN[7,7] JACOBIAN[7,9] JACOBIAN[7,11] JACOBIAN[7,5]
                JACOBIAN[9,7] JACOBIAN[9,9] JACOBIAN[9,11] JACOBIAN[9,5]
                JACOBIAN[11,7] JACOBIAN[11,9] JACOBIAN[11,11] JACOBIAN[11,5]
                JACOBIAN[5,7] JACOBIAN[5,9] JACOBIAN[5,11] JACOBIAN[5,5]]
EigenvalsLo, EigenvecsLo =eigen(longitudinal)
lateral = [JACOBIAN[8,8] JACOBIAN[8,10] JACOBIAN[8,12] JACOBIAN[8,4]
                JACOBIAN[10,8] JACOBIAN[10,10] JACOBIAN[10,12] JACOBIAN[10,4]
                JACOBIAN[12,8] JACOBIAN[12,10] JACOBIAN[12,12] JACOBIAN[12,4]
                JACOBIAN[4,8] JACOBIAN[4,10] JACOBIAN[4,12] JACOBIAN[4,4]]
EigenvalsLa, EigenvecsLa =eigen(lateral)

#Hunting down errors--------------------------------
s1=s0
state = State(s0...)
Fa,Ma = aeroforces(sd, atm, state, controller, ref, mp)
Fp,Mp = propulsionforces(propulsion, atm, state, controller, ref, mp)
Fg,Mg = gravityforces(inertial2, atm, state, controller, ref, mp)
ds = zeros(length(s0))
ds[1:3]=s0[7:9]
#required for use of sixdof!() Not implemented.
time = 0.0
#Also must be reset each time sixdof!() is used.
state0=s0
params0=p
sixdof!(ds, state0, params0, time)#sixdof!() modifies ds so that it is now the defailt derivatives
#Result: values in ds give the derivatives about the center point.
dsc=ds

s1[10]=s1[10]+.001
state1 = State(s1...)
Fap,Map = aeroforces(sd, atm, state1, controller, ref, mp)
Fpp,Mpp = propulsionforces(propulsion, atm, state1, controller, ref, mp)
Fgp,Mgp = gravityforces(inertial2, atm, state1, controller, ref, mp)
#F = Fa + Fp + Fg
time = 0.0
state0=s0
params0=p
ds = zeros(length(s0))
ds[1:3]=s0[7:9]
#Perturb s at correct index
state0[10]=state0[10]+.001
#Run sixdof!() for perturbed value.
sixdof!(ds, state0, params0, time)
dsp=ds#perturbed ds

col11=(dsp-dsc)/.001


text:
# AUTO GENERATED FILE - DO NOT EDIT

export 'feffery'_antdcascader

"""
    'feffery'_antdcascader(;kwargs...)

An AntdCascader component.

Keyword arguments:
- `id` (String; optional)
- `bordered` (Bool; optional)
- `changeOnSelect` (Bool; optional)
- `className` (String; optional)
- `defaultValue` (Array of Strings; optional)
- `disabled` (Bool; optional)
- `expandTrigger` (a value equal to: 'click', 'hover'; optional)
- `loading_state` (optional): . loading_state has the following type: lists containing elements 'is_loading', 'prop_name', 'component_name'.
Those elements have the following types:
  - `is_loading` (Bool; optional): Determines if the component is loading or not
  - `prop_name` (String; optional): Holds which property is loading
  - `component_name` (String; optional): Holds the name of the component that is loading
- `maxTagCount` (Real | a value equal to: 'responsive'; optional)
- `multiple` (Bool; optional)
- `options` (optional)
- `placeholder` (String; optional)
- `placement` (a value equal to: 'bottomLeft', 'bottomRight', 'topLeft', 'topRight'; optional)
- `size` (a value equal to: 'small', 'middle', 'large'; optional)
- `style` (Dict; optional)
- `value` (Array of Bool | Real | String | Dict | Arrays; optional)
"""
function 'feffery'_antdcascader(; kwargs...)
        available_props = Symbol[:id, :bordered, :changeOnSelect, :className, :defaultValue, :disabled, :expandTrigger, :loading_state, :maxTagCount, :multiple, :options, :placeholder, :placement, :size, :style, :value]
        wild_props = Symbol[]
        return Component("'feffery'_antdcascader", "AntdCascader", "feffery_antd_components", available_props, wild_props; kwargs...)
end

max_stars_repo_name: limmyforget/feffery-antd-components

max_issues_repo_name: limmyforget/feffery-antd-components

max_forks_repo_name: limmyforget/feffery-antd-components


text:
# Before running this, please make sure to activate and instantiate the environment
# corresponding to [this `Project.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Project.toml) and [this `Manifest.toml`](https://raw.githubusercontent.com/alan-turing-institute/DataScienceTutorials.jl/master/Manifest.toml)
# so that you get an environment which matches the one used to generate the tutorials:
#
# ```julia
# cd("DataScienceTutorials") # cd to folder with the *.toml
# using Pkg; Pkg.activate("."); Pkg.instantiate()
# ```

# In stacking one blends the predictions of different regressors or# classifiers to gain, in some cases, better performance than naive# averaging or majority vote.
# Here we illustrate how to build a two-model stack as an MLJ learning# network, which we export as a new stand-alone composite model# type `MyTwoStack`. This will make the stack that we build completely# re-usable (new data, new models). This means we can apply# meta-algorithms, such as performance evaluation and tuning, to the# stack, exactly as we would for any other model.
# Our main purpose here is to demonstrate the flexibility of MLJ's# composite model interface. Eventually, MLJ will provide built-in# composite types or macros to achieve the same results in a few# lines, which will suffice for routine stacking tasks.
# ## Basic stacking using out-of-sample base learner predictions
# A rather general stacking protocol was first described in a [1992# paper](https://www.sciencedirect.com/science/article/abs/pii/S0893608005800231)# by David Wolpert. For a generic introduction to the basic two-layer# stack described here, see [this blog# post](https://burakhimmetoglu.com/2016/12/01/stacking-models-for-improved-predictions/)# of Burak Himmetoglu.
# A basic stack consists of a number of base learners (two, in this# illustration) and a single adjudicating model.
# When a stacked model is called to make a prediction, the individual# predictions of the base learners are made the columns of an *input*# table for the adjudicating model, which then outputs the final# prediction. However, it is crucial to understand that the flow of# data *during training* is not the same.
# The base model predictions used to train the adjudicating model are# *not* the predictions of the base learners fitted to all the# training data. Rather, to prevent the adjudicator giving too much# weight to the base learners with low *training* error, the input# data is first split into a number of folds (as in cross-validation),# a base learner is trained on each fold complement individually, and# corresponding predictions on the folds are spliced together to form# a full-length prediction called the *out-of-sample prediction*.
# For illustrative purposes we use just three folds. Each base learner# will get three separate machines, for training on each fold# complement, and a fourth machine, trained on all the supplied data,# for use in the prediction flow.
# We build the learning network with dummy data at the source nodes,# so the reader inspects the workings of the network as it is built (by# calling `fit!` on nodes, and by calling the nodes themselves). As# usual, this data is not seen by the exported composite model type,# and the component models we choose are just default values for the# hyperparameters of the composite model.
using MLJ
using PyPlot

using StableRNGs

# Some models we will use:
linear = @load LinearRegressor pkg=MLJLinearModels
knn = @load KNNRegressor; knn.K = 4
tree_booster = @load EvoTreeRegressor; tree_booster.nrounds = 100
forest = @load RandomForestRegressor pkg=DecisionTree; forest.n_trees = 500
svm = @load SVMRegressor;

# ### Warm-up exercise: Define a model type to average predictions
# Let's define a composite model type `MyAverageTwo` that# averages the predictions of two deterministic regressors. Here's the learning network:
X = source()
y = source()

model1 = linear
model2 = knn

m1 = machine(model1, X, y)
y1 = predict(m1, X)

m2 = machine(model2, X, y)
y2 = predict(m2, X)

yhat = 0.5*y1 + 0.5*y2

# In preparation for export, we wrap the learning network in a# learning network machine, which specifies what the source nodes are,# and which node is for prediction. As our exported model will make# point-predictions (as opposed to probabilistic ones), we use a# `Deterministic` ["surrogate" model](https://alan-turing-institute.github.io/MLJ.jl/dev/performance_measures/):
mach = machine(Deterministic(), X, y; predict=yhat)

# Note that we cannot actually fit this machine because we chose not# to wrap our source nodes `X` and `y` in data.
# Here's the macro call that "exports" the learning network as a new# composite model `MyAverageTwo`:
@from_network mach begin
    mutable struct MyAverageTwo
        regressor1=model1
        regressor2=model2
    end
end

# Note that, unlike a normal struct definition, the defaults `model1`# and `model2` must be specified, and they must refer to model# instances in the learning network.
# We can now create an instance of the new type:
average_two = MyAverageTwo()

# Evaluating this average model on the Boston data set, and comparing# with the base model predictions:
function print_performance(model, data...)
    e = evaluate(model, data...;
                 resampling=CV(rng=StableRNG(1234), nfolds=8),
                 measure=rms,
                 verbosity=0)
    μ = round(e.measurement[1], sigdigits=5)
    ste = round(std(e.per_fold[1])/sqrt(8), digits=5)
    println("$model = $μ ± $(2*ste)")
end;

X, y = @load_boston

print_performance(linear, X, y)
print_performance(knn, X, y)
print_performance(average_two, X, y)

# ## Stacking proper# ### Helper functions:
# To generate folds for generating out-of-sample predictions, we define
folds(data, nfolds) =
    partition(1:nrows(data), (1/nfolds for i in 1:(nfolds-1))...);

# For example, we have:
f = folds(1:10, 3)

# It will also be convenient to use the MLJ method `restrict(X, f, i)`# that restricts data `X` to the `i`th element (fold) of `f`, and# `corestrict(X, f, i)` that restricts to the corresponding fold# complement (the concatenation of all but the `i`th# fold).
# For example, we have:
corestrict(string.(1:10), f, 2)

# ### Choose some test data (optional) and some component models (defaults for the composite model):
figure(figsize=(8,6))
steps(x) = x < -3/2 ? -1 : (x < 3/2 ? 0 : 1)
x = Float64[-4, -1, 2, -3, 0, 3, -2, 1, 4]
Xraw = (x = x, )
yraw = steps.(x);
idxsort = sortperm(x)
xsort = x[idxsort]
ysort = yraw[idxsort]
step(xsort, ysort, label="truth", where="mid")
plot(x, yraw, ls="none", marker="o", label="data")
xlim(-4.5, 4.5)
legend()

# \fig{s1.svg}
# Some models to stack (which we can change later):
model1 = linear
model2 = knn

# The adjudicating model:
judge = linear

# ### Define the training nodes
# Let's instantiate some input and target source nodes for the# learning network, wrapping the play data defined above in source# nodes:
X = source(Xraw)
y = source(yraw)

# Our first internal node will represent the three folds (vectors of row# indices) for creating the out-of-sample predictions. We would like# to define `f = folds(X, 3)` but this will not work because `X` is# not a table, just a node representing a table. We could fix this# by using the @node macro:
f = @node folds(X, 3)

# Now `f` is itself a node, and so callable:
f()

# However, we can also just overload `folds` to work on nodes, using the# `node` *function*:
folds(X::AbstractNode, nfolds) = node(XX->folds(XX, nfolds), X)
f = folds(X, 3)
f()

# In the case of `restrict` and `corestrict`, which also don't# operate on nodes, method overloading will save us writing `@node`# all the time:
MLJ.restrict(X::AbstractNode, f::AbstractNode, i) =
    node((XX, ff) -> restrict(XX, ff, i), X, f);
MLJ.corestrict(X::AbstractNode, f::AbstractNode, i) =
    node((XX, ff) -> corestrict(XX, ff, i), X, f);

# We are now ready to define machines for training `model1` on each# fold-complement:
m11 = machine(model1, corestrict(X, f, 1), corestrict(y, f, 1))
m12 = machine(model1, corestrict(X, f, 2), corestrict(y, f, 2))
m13 = machine(model1, corestrict(X, f, 3), corestrict(y, f, 3))

# Define each out-of-sample prediction of `model1`:
y11 = predict(m11, restrict(X, f, 1));
y12 = predict(m12, restrict(X, f, 2));
y13 = predict(m13, restrict(X, f, 3));

# Splice together the out-of-sample predictions for model1:
y1_oos = vcat(y11, y12, y13);

# Note there is no need to overload the `vcat` function to work on# nodes; it does so out of the box, as does `hcat` and basic# arithmetic operations.
# Since our source nodes are wrapping data, we can optionally check# our network so far, by calling fitting and calling `y1_oos`:
fit!(y1_oos, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, y1_oos(), ls="none", marker="o", label="linear oos")
legend()

# \fig{s2.svg}
# We now repeat the procedure for the other model:
m21 = machine(model2, corestrict(X, f, 1), corestrict(y, f, 1))
m22 = machine(model2, corestrict(X, f, 2), corestrict(y, f, 2))
m23 = machine(model2, corestrict(X, f, 3), corestrict(y, f, 3))
y21 = predict(m21, restrict(X, f, 1));
y22 = predict(m22, restrict(X, f, 2));
y23 = predict(m23, restrict(X, f, 3));

# And testing the knn out-of-sample prediction:
y2_oos = vcat(y21, y22, y23);
fit!(y2_oos, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, y2_oos(), ls="none", marker="o", label="knn oos")
legend()

# \fig{s3.svg}
# Now that we have the out-of-sample base learner predictions, we are# ready to merge them into the adjudicator's input table and construct# the machine for training the adjudicator:
X_oos = MLJ.table(hcat(y1_oos, y2_oos))
m_judge = machine(judge, X_oos, y)

# Are we done with constructing machines? Well, not quite. Recall that# when we use the stack to make predictions on new data, we will be# feeding the adjudicator ordinary predictions of the base learners# (rather than out-of-sample predictions). But so far, we have only# defined machines to train the base learners on fold complements, not# on the full data, which we do now:
m1 = machine(model1, X, y)
m2 = machine(model2, X, y)

# ### Define nodes still needed for prediction
# To obtain the final prediction, `yhat`, we get the base learner# predictions, based on training with all data, and feed them to the# adjudicator:
y1 = predict(m1, X);
y2 = predict(m2, X);
X_judge = MLJ.table(hcat(y1, y2))
yhat = predict(m_judge, X_judge)

# Let's check the final prediction node can be fit and called:
fit!(yhat, verbosity=0)

figure(figsize=(8,6))
step(xsort, ysort, label="truth", where="mid")
plot(x, yhat(), ls="none", marker="o", label="yhat")
legend()

# \fig{s4}
# Although of little statistical significance here, we note that# stacking gives a lower *training* error than naive averaging:
e1 = rms(y1(), y())
e2 = rms(y2(), y())
emean = rms(0.5*y1() + 0.5*y2(), y())
estack = rms(yhat(), y())
@show e1 e2 emean estack;

# ## Export the learning network as a new model type
# The learning network (less the data wrapped in the source nodes)# amounts to a specification of a new composite model type for# two-model stacks, trained with three-fold resampling of base model# predictions. Let's create the new type `MyTwoModelStack`, in the# same way we exported the network for model averaging:
@from_network machine(Deterministic(), X, y; predict=yhat) begin
    mutable struct MyTwoModelStack
        regressor1=model1
        regressor2=model2
        judge=judge
    end
end

my_two_model_stack = MyTwoModelStack()

# And this completes the definition of our re-usable stacking model type.
# ## Applying `MyTwoModelStack` to some data
# Without undertaking any hyperparameter optimization, we evaluate the# performance of a tree boosting algorithm and a support vector# machine on a synthetic data set. As adjudicator, we'll use a random# forest.
# We use a synthetic set to give an example where stacking is# effective but the data is not too large. (As synthetic data is based# on perturbations to linear models, we are deliberately avoiding# linear models in stacking illustration.)
X, y = make_regression(1000, 20; sparse=0.75, noise=0.1, rng=123);

# #### Define the stack and compare performance
avg = MyAverageTwo(regressor1=tree_booster,
                   regressor2=svm)

stack = MyTwoModelStack(regressor1=tree_booster,
                        regressor2=svm,
                        judge=forest) # judge=linear

all_models = [tree_booster, svm, forest, avg, stack];

for model in all_models
    print_performance(model, X, y)
end

# #### Tuning a stack
# A standard abuse of good data hygiene is to optimize stack component# models *separately* and then tune the adjudicating model# hyperparameters (using the same resampling of the data) with the# base learners fixed. Although more computationally expensive, better# generalization might be expected by applying tuning to the stack as# a whole, either simultaneously, or in in cheaper sequential# steps. Since our stack is a stand-alone model, this is readily# implemented.
# As a proof of concept, let's see how to tune one of the base model# hyperparameters, based on performance of the stack as a whole:
r = range(stack, :(regressor2.C), lower = 0.01, upper = 10, scale=:log)
tuned_stack = TunedModel(model=stack,
                         ranges=r,
                         tuning=Grid(shuffle=false),
                         measure=rms,
                         resampling=Holdout())

mach = fit!(machine(tuned_stack,  X, y), verbosity=0)
best_stack = fitted_params(mach).best_model
best_stack.regressor2.C

# Let's evaluate the best stack using the same data resampling used to# the evaluate the various untuned models earlier (now we are neglecting# data hygiene!):
print_performance(best_stack, X, y)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
module WeightedStats

export weighted_mean, weighted_median

include("mean.jl")

include("median.jl")

end # module