text:
using FTCTests
const FTC = FaultTolerantControl
using LinearAlgebra
using Plots

"""
# Refs
[1] K. G. Vamvoudakis and F. L. Lewis, “Online actor-critic algorithm to solve the continuous-time infinite horizon optimal control problem,” Automatica, vol. 46, no. 5, pp. 878–888, 2010, doi: 10.1016/j.automatica.2010.02.018.

# Notes
- x = state vector, [v, alpha, q]^T
- u = control input vector, [dele]
"""
function compute_minHSV(eff)
    # F16 aircraft longitudinal linear system model [1]
    A = [-1.01887 0.90506 -0.00215;
    0.82225 -1.07741 -0.17555;
    0 0 -1]
    B = [0 0 1]'
    C = Matrix(I, 3, 3)

    f(x, u, p, t) = A*x + B*(eff * u)
    g(x, u, p, t) = C*x

	# System dimension
	n = size(A)[1]
	m = 1
	l = size(C)[1]

    # Initial settings
    x0 = zeros(n,)
    u0 = zeros(m,)
    dt = 0.001
    tf = 1.0
    pr = zeros(3, 1)

    Wc = empirical_gramian(f, g, m, n, l; opt=:c, dt=dt, tf=tf, pr=pr, xs=x0, us=u0)
    Wo = empirical_gramian(f, g, m, n, l; opt=:o, dt=dt, tf=tf, pr=pr, xs=x0, us=u0)
	minHSV = min_HSV(Wc, Wo)
end

function plotting()
    lambda = range(0, 1, 20)
    HSVs = []
    for i = 1:length(lambda)
        HSVs = push!(HSVs, compute_minHSV(collect(lambda)[i]))
    end
    plot(lambda,
        HSVs,
        xlabel = "Actuator effectiveness",
        ylabel = "Minimum HSVs",
        label = nothing,
    )
end
nothing


text:
"""
    apply_Ω(δψ, ψ, H::Hamiltonian, Λ)

Compute the application of Ω defined at ψ to δψ. H is the Hamiltonian computed
from ψ and Λ is the set of Rayleigh coefficients ψk' * Hk * ψk at each k-point.
"""
function apply_Ω(δψ, ψ, H::Hamiltonian, Λ)
    δψ = proj_tangent(δψ, ψ)
    Ωδψ = [H.blocks[ik] * δψk - δψk * Λ[ik] for (ik, δψk) in enumerate(δψ)]
    proj_tangent!(Ωδψ, ψ)
end

"""
    apply_K(basis::PlaneWaveBasis, δψ, ψ, ρ, occupation)

Compute the application of K defined at ψ to δψ. ρ is the density issued from ψ.
δψ also generates a δρ, computed with `compute_δρ`.
"""
@views function apply_K(basis::PlaneWaveBasis, δψ, ψ, ρ, occupation)
    δψ = proj_tangent(δψ, ψ)
    δρ = compute_δρ(basis, ψ, δψ, occupation)
    δV = apply_kernel(basis, δρ; ρ=ρ)

    Kδψ = map(enumerate(ψ)) do (ik, ψk)
        kpt = basis.kpoints[ik]
        δVψk = similar(ψk)

        for n = 1:size(ψk, 2)
            ψnk_real = G_to_r(basis, kpt, ψk[:, n])
            δVψnk_real = δV[:, :, :, kpt.spin] .* ψnk_real
            δVψk[:, n] = r_to_G(basis, kpt, δVψnk_real)
        end
        δVψk
    end
    # ensure projection onto the tangent space
    proj_tangent!(Kδψ, ψ)
end

"""
    solve_ΩplusK(basis::PlaneWaveBasis{T}, ψ, res, occupation;
                 tol=1e-10, verbose=false) where T

Return δψ where (Ω+K) δψ = rhs
"""
function solve_ΩplusK(basis::PlaneWaveBasis{T}, ψ, rhs, occupation;
                      tol=1e-10, verbose=false) where T
    @assert mpi_nprocs() == 1  # Distributed implementation not yet available
    filled_occ = filled_occupation(basis.model)
    # for now, all orbitals have to be fully occupied -> need to strip them beforehand
    @assert all(all(occ_k .== filled_occ) for occ_k in occupation)

    # compute quantites at the point which define the tangent space
    ρ = compute_density(basis, ψ, occupation)
    _, H = energy_hamiltonian(basis, ψ, occupation; ρ=ρ)

    pack(ψ) = reinterpret_real(pack_ψ(ψ))
    unpack(x) = unpack_ψ(reinterpret_complex(x), size.(ψ))
    unsafe_unpack(x) = unsafe_unpack_ψ(reinterpret_complex(x), size.(ψ))

    # project rhs on the tangent space before starting
    proj_tangent!(rhs, ψ)
    rhs_pack = pack(rhs)

    # preconditioner
    Pks = [PreconditionerTPA(basis, kpt) for kpt in basis.kpoints]
    for ik = 1:length(Pks)
        precondprep!(Pks[ik], ψ[ik])
    end
    function f_ldiv!(x, y)
        δψ = unpack(y)
        proj_tangent!(δψ, ψ)
        Pδψ = [ Pks[ik] \ δψk for (ik, δψk) in enumerate(δψ)]
        proj_tangent!(Pδψ, ψ)
        x .= pack(Pδψ)
    end

    # Rayleigh-coefficients
    Λ = [ψk'Hψk for (ψk, Hψk) in zip(ψ, H * ψ)]

    # mapping of the linear system on the tangent space
    function ΩpK(x)
        δψ = unsafe_unpack(x)
        Kδψ = apply_K(basis, δψ, ψ, ρ, occupation)
        Ωδψ = apply_Ω(δψ, ψ, H, Λ)
        pack(Ωδψ + Kδψ)
    end
    J = LinearMap{T}(ΩpK, size(rhs_pack, 1))

    # solve (Ω+K) δψ = rhs on the tangent space with CG
    δψ, history = cg(J, rhs_pack, Pl=FunctionPreconditioner(f_ldiv!),
                  reltol=0, abstol=tol, verbose=verbose, log=true)

    (; δψ=unpack(δψ), history)
end

# Solves Ω+K using a split algorithm
# With χ04P = -Ω^-1,
# (Ω+K)^-1 = Ω^-1 (1 - K(1+Ω^-1 K)^-1 Ω^-1)
# (Ω+K)^-1 = -χ04P (1 + K(1 - χ04P K)^-1 χ04P)
# (Ω+K)^-1 = -χ04P (1 + E K2P (1 - χ02P K2P)^-1 R χ04P)
# where χ02P = R χ04P E and K2P = R K E
function solve_ΩplusK_split(ham::Hamiltonian, ρ::AbstractArray{T}, ψ, occupation, εF,
                            eigenvalues, rhs; tol=1e-8, tol_sternheimer=tol/10,
                            verbose=false, kwargs...) where T
    basis = ham.basis

    # compute δρ0 (ignoring interactions)
    δψ0 = apply_χ0_4P(ham, ψ, occupation, εF, eigenvalues, rhs;
                      reltol=0, abstol=tol_sternheimer, kwargs...)
    δρ0 = compute_δρ(basis, ψ, δψ0, occupation)

    pack(δρ)   = vec(δρ)
    unpack(δρ) = reshape(δρ, size(ρ))
    # compute total δρ
    function eps_fun(δρ)
        δρ = unpack(δρ)
        δV = apply_kernel(basis, δρ; ρ)
        # Would be nice to play with abstol / reltol etc. to avoid over-solving
        # for the initial GMRES steps.
        χ0δV = apply_χ0(ham, ψ, εF, eigenvalues, δV;
                        abstol=tol_sternheimer, reltol=0, kwargs...)
        pack(δρ - χ0δV)
    end
    J = LinearMap{T}(eps_fun, length(pack(δρ0)))
    δρ, history = gmres(J, pack(δρ0); reltol=0, abstol=tol, verbose, log=true)
    δV = apply_kernel(basis, unpack(δρ); ρ)

    δVψ = [DFTK.RealSpaceMultiplication(basis, kpt, @views δV[:, :, :, kpt.spin]) * ψ[ik]
           for (ik, kpt) in enumerate(basis.kpoints)]
    δψ1 = apply_χ0_4P(ham, ψ, occupation, εF, eigenvalues, δVψ;
                      abstol=tol_sternheimer, reltol=0, kwargs...)
    δψ  = .- (δψ0 .+ δψ1)
    (; δψ, history)
end

function solve_ΩplusK_split(basis::PlaneWaveBasis, ψ, rhs, occupation; kwargs...)
    ρ = compute_density(basis, ψ, occupation)
    _, H = energy_hamiltonian(basis, ψ, occupation; ρ)

    eigenvalues = [real.(eigvals(ψk'Hψk)) for (ψk, Hψk) in zip(ψ, H * ψ)]
    occupation, εF = compute_occupation(basis, eigenvalues)

    solve_ΩplusK_split(H, ρ, ψ, occupation, εF, eigenvalues, rhs; kwargs...)
end

function solve_ΩplusK_split(scfres::NamedTuple, rhs; kwargs...)
    n_ep_extra = scfres.n_ep_extra

    ψ       = [@view ψk[:, 1:end-n_ep_extra]     for ψk in scfres.ψ]
    ψ_extra = [@view ψk[:, end-n_ep_extra+1:end] for ψk in scfres.ψ]

    evals   = [εk[1:end-n_ep_extra]     for εk in scfres.eigenvalues]
    ε_extra = [εk[end-n_ep_extra+1:end] for εk in scfres.eigenvalues]

    occupation = [occk[1:end-n_ep_extra]     for occk in scfres.occupation]
    occ_extra  = [occk[end-n_ep_extra+1:end] for occk in scfres.occupation]

    @assert length(rhs) == length(ψ)
    rhs_capped = map(zip(ψ, rhs)) do (ψk, rhsk)
        n_bands = size(ψk, 2)
        @assert size(rhsk, 2) >= n_bands
        @view rhsk[:, 1:n_bands]
    end
    solve_ΩplusK_split(scfres.ham, scfres.ρ, ψ, occupation, scfres.εF, evals, rhs_capped;
                       ψ_extra, ε_extra, kwargs...)
end

ext: jl


text:
@inline function ExaTronKernels.nrm2!(wa, A::ROCDeviceArray{Float64,2}, n::Int)
    tx = workitemIdx().x

    v = 0.0
    if tx <= n
        @inbounds for j=1:n
            v += A[j,tx]^2
        end
        @inbounds wa[tx] = sqrt(v)
    end
    #=
    v = A[tx,ty]^2

    if tx > n || ty > n
        v = 0.0
    end

    # Sum over the x-dimension.
    offset = div(blockDim().x, 2)
    while offset > 0
        v += CUDA.shfl_down_sync(0xffffffff, v, offset)
        offset = div(offset, 2)
    end

    if tx == 1
        wa[ty] = sqrt(v)
    end
    =#
    AMDGPU.sync_workgroup()

    return
end

@inline function ExaTronKernels.dssyax(n::Int, A::ROCDeviceArray{Float64,2},
                        z::ROCDeviceArray{Float64,1},
                        q::ROCDeviceArray{Float64,1})
    tx = workitemIdx().x

    v = 0.0
    if tx <= n
        @inbounds for j=1:n
            v += A[tx,j]*z[j]
        end
        @inbounds q[tx] = v
    end
    #=
    v = 0.0
    if tx <= n && ty <= n
        v = A[ty,tx]*z[tx]
    end

    # Sum over the x-dimension: v = sum_tx A[ty,tx]*z[tx].
    # The thread with tx=1 will have the sum in v.

    offset = div(blockDim().x, 2)
    while offset > 0
        v += CUDA.shfl_down_sync(0xffffffff, v, offset)
        offset = div(offset, 2)
    end

    if tx == 1
        q[ty] = v
    end
    =#
    AMDGPU.sync_workgroup()

    return
end

@inline function ExaTronKernels.reorder!(n::Int, nfree::Int, B::ROCDeviceArray{Float64,2},
                          A::ROCDeviceArray{Float64,2}, indfree::ROCDeviceArray{Int,1},
                          iwa::ROCDeviceArray{Int,1})
    tx = workitemIdx().x

    #=
    if tx == 1 && ty == 1
        @inbounds for j=1:nfree
            jfree = indfree[j]
            B[j,j] = A[jfree,jfree]
            for i=jfree+1:n
                if iwa[i] > 0
                    B[iwa[i],j] = A[i,jfree]
                    B[j,iwa[i]] = B[iwa[i],j]
                end
            end
        end
    end
    =#
    if tx <= nfree
        @inbounds begin
            jfree = indfree[tx]
            B[tx,tx] = A[jfree,jfree]
            for i=jfree+1:n
                if iwa[i] > 0
                    B[iwa[i],tx] = A[i,jfree]
                    B[tx,iwa[i]] = B[iwa[i],tx]
                end
            end
        end
    end

    AMDGPU.sync_workgroup()

    return
end


text:
function variable_as_node!(ex, node::Type{T} where {T <: AbstractNode})
    ex
end

function variable_as_node!(ex::Symbol, node::Type{T} where {T <: AbstractNode})
    Expr(:call, node, QuoteNode(ex))
end

function variable_as_node!(ex::Expr, node::Type{T} where {T <: AbstractNode})
    to_quote = eachindex(ex.args)

    if ex.head == :call
        if ex.args[1] == :_
            length(ex.args) > 2 ? error("Unqote only a single argument. Right: `_(x)`, Wrong: `_(x, y)`.") : nothing
            return ex.args[2]
        end
        to_quote = to_quote[2:end]
    end
    for i in to_quote
        ex.args[i] = variable_as_node!(ex.args[i], node)
    end
    ex
end

variable_as_node!(ex) = variable_as_node!(ex, SimpleNode)

macro variable_as_node(ex)
    esc(variable_as_node!(ex, SimpleNode))
end


text:
macro component!(x)
    SM.terms!(SM.sort_terms!(SM.parse!(x)))
end

function component(x::Expr)

    (x.args[1] != :(=>)) && throw(ArgumentError("invalid expression encountered: $x"))

    y = @eval @component! $(x.args[3])

    if (typeof(y) <: Tuple) && (typeof(y[1]) <: SM.ConstantTerm)
        y = (y[2:end]..., y[1])
    end

    return y
end

macro micromodel(x...)
    Dict{Symbol, SM.TermOrTerms}(xi.args[2] => component(xi) for xi in x)
end

terms(y::Dict{Symbol, SM.TermOrTerms})    = mapreduce(terms, union, values(y))
termvars(y::Dict{Symbol, SM.TermOrTerms}) = mapreduce(termvars, union, values(y))
eterms(y::Dict{Symbol, SM.TermOrTerms})   = Tuple(mapreduce(SM.vectorize, union, values(y)))

function schema(y::Dict{Symbol, SM.TermOrTerms}, df::AbstractDataFrame, hints::Dict)
    schema(filter(SM.needs_schema, terms(y)), df, hints)
end

function apply_schema(y::Dict{Symbol, SM.TermOrTerms}, schema::SM.Schema)
    Dict{Symbol, SM.TermOrTerms}(i => apply_schema(v, schema) for (i, v) in y)
end

has_schema(y::Dict{Symbol, SM.TermOrTerms}) = all(has_schema, values(y))

function modelcols(y::Dict{Symbol, SM.TermOrTerms}, df::AbstractDataFrame)
    modelcols(SM.collect_matrix_terms(eterms(y)), df)
end

function asgn(y::Dict{Symbol, SM.TermOrTerms})

    Y = eterms(y)
    Y = Y[SM.asgn(Y)]

    return Dict{Symbol, Vector{Int}}(i => findall((in)(SM.vectorize(v)), Y) for (i, v) in y)
end

function asgn(y::Dict{Symbol, SM.TermOrTerms}, t::SM.TermOrTerms)

    Y = eterms(y)
    Y = Y[SM.asgn(Y)]

    return [findfirst(isequal(ti), Y) for ti in t]
end

max_stars_repo_name: lbittarello/Microeconometrics

max_issues_repo_name: lbittarello/Microeconometrics

max_forks_repo_name: lbittarello/Microeconometrics


text:
abstract type ACall <: AOperator end
function (f::Node)(x::Node)
    return Node("()", ACall, (f, x))
end
function (f::Node)(Vararg...)
    return Node("()", f, Vararg)
end


text:
using Sobol

export LandauDamping

"""
    Landau( α, kx)

Test structure to initialize a particles distribtion for
Landau damping test case in 1D1V and 1D2V

"""
struct LandauDamping

    alpha::Float64
    kx::Float64

end

function newton(r, alpha, kx)
    x0, x1 = 0.0, 1.0
    r *= 2π / kx
    while (abs(x1 - x0) > 1e-12)
        p = x0 + alpha * sin(kx * x0) / kx
        f = 1 + alpha * cos(kx * x0)
        x0, x1 = x1, x0 - (p - r) / f
    end
    x1
end

"""
    sample!(d, pg)

Sampling from a probability distribution to initialize a Landau damping in
1D2V space.

```math
f_0(x,v,t) = \\frac{n_0}{2π v_{th}^2} ( 1 + \\alpha cos(k_x x))
 exp( - \\frac{v_x^2+v_y^2}{2 v_{th}^2})
```
The newton function solves the equation ``P(x)-r=0`` with Newton’s method
```math
x^{n+1} = x^n – (P(x)-(2\\pi r / k)/f(x)
```
with
```math
P(x) = \\int_0^x (1 + \\alpha cos(k_x y)) dy = x + \\frac{\\alpha}{k_x} sin(k_x x)
```
"""
function sample!(pg::ParticleGroup{1,2}, d::LandauDamping)

    alpha, kx = d.alpha, d.kx

    s = SobolSeq(1)
    nbpart = pg.n_particles

    for i = 1:nbpart
        v = sqrt(-2 * log((i - 0.5) / nbpart))
        x = Sobol.next!(s)
        set_x!(pg, i, newton(x, alpha, kx))
        set_v!(pg, i, v)
        set_w!(pg, i, 2 * pi / kx / nbpart)
    end

end

"""
    sample!(d, pg)

Sampling from a probability distribution to initialize a Landau damping in
2D2V space.

```math
f_0(x,v,t) = \\frac{n_0}{2π v_{th}^2} ( 1 + \\alpha cos(k_x x))
 exp( - \\frac{v_x^2+v_y^2}{2 v_{th}^2})
```
"""
function sample!(pg::ParticleGroup{2,2}, mesh :: TwoDGrid, d::LandauDamping)

    alpha, kx = d.alpha, d.kx
    @assert mesh.dimx ≈ 2π / kx

    s = Sobol.SobolSeq(3)
    nbpart = pg.n_particles

    for i = 1:nbpart
        v = sqrt(-2 * log((i - 0.5) / nbpart))
        x, y, θ = Sobol.next!(s)
        θ = θ * 2π
        pg.array[1, i] = newton(x, alpha, kx)
        pg.array[2, i] = mesh.ymin + y * mesh.dimy
        pg.array[3, i] = v * cos(θ)
        pg.array[4, i] = v * sin(θ)
        pg.array[5, i] = (mesh.dimx * mesh.dimy) / nbpart
    end

end

"""
    sample!(d, pg)

Sampling from a probability distribution to initialize a Landau damping in
1D1V space.

```math
f_0(x,v,t) = \\frac{n_0}{2π v_{th}^2} ( 1 + \\alpha cos(k_x x)) exp( - \\frac{v^2}{2 v_{th}^2})
```
"""
function sample!( pg::ParticleGroup{1,1}, mesh :: OneDGrid, d::LandauDamping)

    s = Sobol.SobolSeq(2)
    alpha, kx = d.alpha, d.kx
    @assert mesh.dimx ≈ 2π / kx

    nbpart = pg.n_particles

    for i = 1:nbpart
        v = sqrt(-2 * log((i - 0.5) / nbpart))
        r1, r2 = Sobol.next!(s)
        θ = r1 * 2π
        pg.array[1,i] = newton(r2, alpha, kx)
        pg.array[2,i] = v * sin(θ)
        pg.array[3,i] = mesh.dimx / nbpart
    end

end


text:
using FEMSolids
using Documenter

DocMeta.setdocmeta!(FEMSolids, :DocTestSetup, :(using FEMSolids); recursive=true)

makedocs(;
    modules=[FEMSolids],
    authors="kimauth <auth@chalmers.se> and contributors",
    repo="https://github.com/kimauth/FEMSolids.jl/blob/{commit}{path}#{line}",
    sitename="FEMSolids.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kimauth.github.io/FEMSolids.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/kimauth/FEMSolids.jl.git",
    devbranch="main",
    push_preview=true,
)


text:
using SymEngine
using SymEngine: BasicType

sconj = SymFunction("conj")
Base.conj(x::Basic) = Basic(conj(SymEngine.BasicType(x)))
Base.conj(x::BasicType) = real(x) - im * imag(x)
Base.imag(x::BasicType{Val{:Constant}}) = Basic(0)
Base.imag(x::BasicType{Val{:Symbol}}) = Basic(0)

pretty_print_number(x; lengthonly=false) = pretty_print_number(stdout, x; lengthonly=lengthonly)
function pretty_print_number(io::IO, x; lengthonly=false)
    sx = string(x)
    lengthonly || print(io, sx)
    return length(sx)
end

function pretty_print_number(io::IO, x::AbstractFloat; lengthonly=false)
    closest_int = round(Int, x)
    if isapprox(x, closest_int, atol=1e-12)
        si = string(closest_int)
        lengthonly || print(io, si)
        return length(si)
    else
        sx = string(x)
        lengthonly || print(io, sx)
        return length(sx)
    end
end

function pretty_print_number(io::IO, x::Complex; atol::Real = 1e-12, lengthonly=false)
    l = 0
    if !isapprox(real(x), 0, atol=atol)
        l += pretty_print_number(io, real(x), lengthonly=lengthonly)
    end
    if !isapprox(imag(x), 0, atol=atol)
        if !isapprox(real(x), 0, atol=atol)
            lengthonly || print(imag(x) > 0 ? "+" : "")
            l += 1
        end
        l += pretty_print_number(io, imag(x), lengthonly=lengthonly)
        lengthonly || print(io, "I")
        l += 1
    else
        if isapprox(real(x), 0, atol=atol)
            lengthonly || print(io, "0")
            l += 1
        end
    end
    return l
end

pretty_print_matrix(m) = pretty_print_matrix(stdout, m)
function pretty_print_matrix(io::IO, m)
    minlen = maximum(pretty_print_number.(m, lengthonly=true))+1
    for i in 1:size(m,1)
        print(io, "[")
        for j in 1:size(m,2)
            l = pretty_print_number(m[i,j])
            print(" "^(minlen-l-(j==size(m,1))))
        end
        println(io, "]")
    end
end


text:
using Test
using DataFrames
using CSV

@testset "GraphNodeFeaturization" begin

    # custom lookup table...
    test_df = CSV.read(abspath(@__DIR__, "..", "test_data", "lookup_table.csv"), DataFrame)

    @testset "Construct, encode, decode" begin
        # make sure both constructors give the same results
        local fnames = ["X", "Block", "Atomic mass"]
        features = ElementFeatureDescriptor.(fnames)
        fzn1 = GraphNodeFeaturization(features)
        fzn2 = GraphNodeFeaturization(fnames)

        triangle_C_1 = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])
        triangle_C_2 = AtomGraph(Float32.([0 1 1; 1 0 1; 1 1 0]), ["C", "C", "C"])
        # (this one also tests broadcasting of featurize!)
        featurize!.([triangle_C_1, triangle_C_2], [fzn1, fzn2])
        @test triangle_C_1.encoded_features == triangle_C_2.encoded_features

        # test that other options work properly
        fzn3 = GraphNodeFeaturization(fnames, nbins = 2)
        F2 = AtomGraph(Float32.([0 1; 1 0]), ["F", "F"])
        featurize!(F2, fzn3)
        decoded_matrix = decode(fzn3, F2.encoded_features)
        decoded_ag = decode(F2)
        @test all(
            map(d -> d[1]["Block"] == d[2]["Block"] == "p", [decoded_matrix, decoded_ag]),
        )
    end

    # encodable_elements
    @testset "Encodable Elements" begin
        fzn3 = GraphNodeFeaturization(ElementFeatureDescriptor.(["Boiling point", "6d"]))
        fzn4 = GraphNodeFeaturization(ElementFeatureDescriptor.(["7s", "Valence"]))
        @test encodable_elements(fzn3) == ["Ac", "Th", "U"]
        @test encodable_elements(fzn4) == ["Fr", "Ra", "Ac", "Th"]

        # Custom lookup_table
        feature_1 = ElementFeatureDescriptor("MeaningOfLife", test_df)   # zero-value case - `As` has a value = 0
        feature_2 = ElementFeatureDescriptor("neg_nums", test_df)
        feature_3 = ElementFeatureDescriptor("first_letter", test_df)
        feature_4 = ElementFeatureDescriptor("noarsenic", test_df) # missing value case - `As` has a missing value

        @test encodable_elements(GraphNodeFeaturization([feature_1, feature_2])) ==
              ["C", "As", "Tc"]
        @test encodable_elements(GraphNodeFeaturization([feature_1, feature_4])) ==
              ["C", "Tc"]
        @test encodable_elements(GraphNodeFeaturization([feature_2, feature_3])) ==
              ["C", "As", "Tc"]
        @test encodable_elements(GraphNodeFeaturization([feature_2, feature_4])) ==
              ["C", "Tc"]
    end

    # chunk_vec helper fcn
    @testset "chunk_vec" begin
        vec = [1, 1, 0, 1, 0, 1, 0]
        @test_throws AssertionError chunk_vec(vec, [3, 3])
        @test chunk_vec(vec, [4, 1, 2]) == [[1, 1, 0, 1], [0], [1, 0]]
    end

end


text:
using Test, PseudoArcLengthContinuation, KrylovKit, RecursiveArrayTools, LinearAlgebra
const PALC = PseudoArcLengthContinuation
N = 10
M = 5
δ = 1e-8
####################################################################################################
# test the jacobian of the multiple shooting functional using Linear flow
# TODO do example with A matrix and exp(At)
vf(x) = x
flow(x, t) = exp(t) .* x
dflow(x, dx, t) = (flow(x, t), exp(t) .* dx)
section(x) = dot(x[1], ones(length(x[1])))
section(x::BorderedArray) = section(x.u)

fl = PALC.Flow(vf, flow, dflow)

probSh = PALC.ShootingProblem(fl,
	LinRange(0, 1, M+1) |> diff ,
	x -> section(x))

poguess = VectorOfArray([rand(N) for ii=1:M])
	po = BorderedArray(poguess, 1.)

res = probSh(po)

dpoguess = VectorOfArray([rand(N) for ii=1:M])
	dpo = BorderedArray(dpoguess, 2.)

dres = probSh(po, dpo)

# computation of finite differences
poguess2 = VectorOfArray([po.u[ii] .+ δ .* dpo.u[ii] for ii=1:M])
po2 = BorderedArray(poguess2, po.p + δ * dpo.p)
res2 = probSh(po2)

dresfd = BorderedArray((res2.u - res.u) ./δ, (res2.p - res.p)/δ)

@test dresfd.p - dres.p ≈ 0.0
@test norm(dresfd.u - dres.u, Inf) < 10δ

# use of AbstractArray structure
pov = vcat(vec(po.u), po.p)
dpov = vcat(vec(dpo.u), dpo.p)
resv = probSh(pov)
@test norm(res.u[1] - resv[1:N] , Inf) ≈ 0
dres = probSh(pov, dpov)
resfd = (probSh(pov .+ δ .* dpov) .- resv) ./ δ
@test norm(resfd - dres, Inf) < 10δ
####################################################################################################
# test the jacobian of the multiple shooting functional using nonLinear flow
vf(x) = x.^2
flow(x, t) = x ./ (1 .- t .* x)
dflow(x, dx, t) = (flow(x, t), dx ./ (1 .- t .* x).^2)
section(x) = dot(x[1], ones(length(x[1])))

fl = PALC.Flow(vf, flow, dflow)

probSh = PALC.ShootingProblem(fl,
	LinRange(0,1,M+1) |> diff ,
	x -> section(x))

poguess = VectorOfArray([rand(N) for ii=1:M])
	po = BorderedArray(poguess, 1.)

res = probSh(po)

dpoguess = VectorOfArray([rand(N) for ii=1:M])
	dpo = BorderedArray(dpoguess, 2.)

dres = probSh(po, dpo)

# computation of finite differences
poguess2 = VectorOfArray([po.u[ii] .+ δ .* dpo.u[ii] for ii=1:M])
po2 = BorderedArray(poguess2, po.p + δ * dpo.p)
res2 = probSh(po2)

dresfd = BorderedArray((res2.u - res.u) ./δ, (res2.p - res.p)/δ)

@test dresfd.p - dres.p ≈ 0.0
@test norm(dresfd.u - dres.u, Inf) < 10δ

# use of AbstractArray structure
pov = vcat(vec(po.u), po.p)
dpov = vcat(vec(dpo.u), dpo.p)
resv = probSh(pov)
@test norm(res.u[1] - resv[1:N] , Inf) ≈ 0
dres = probSh(pov, dpov)
resfd = (probSh(pov .+ δ .* dpov) .- resv) ./ δ
@test norm(resfd - dres, Inf) < 10δ
####################################################################################################
# test the hyperplane projections for Poincare Shooting
M = 3
normals = [rand(10) for ii=1:M]
centers = [rand(10) for ii=1:M]

hyper = PALC.HyperplaneSections(normals, centers)

x = 1:10 |>collect .|> Float64 |> vec
xb = PALC.R(hyper, x, 1)
# test
x2 = PALC.E(hyper, xb, 1)
# test that x2 in Sigma2
@test dot(x2 - centers[1], normals[1]) < 1e-14
# that we have Rk∘Ek = Id and Ek∘Rk = IdΣ
@test PALC.R(hyper, x2, 1) ≈ xb


text:
module StaticPolynomials

import Base: *, +, ^,
    convert, promote_rule, convert, show, isless, literal_pow,
    one, copy

export @polyvar,
       Variable,
       Monomial,
       Term,
       Polynomial

abstract type PolynomialLike end
abstract type TermLike <: PolynomialLike end
abstract type MonomialLike <: TermLike end
abstract type VariableLike <: MonomialLike end

immutable Variable{Name} <: VariableLike
end

macro polyvar(names...)
    exprs = [
        quote
            $(esc(name)) = Variable{$(esc(Expr(:quote, name)))}()
        end
        for name in names
    ]
    Expr(:block, exprs...)
end

one(v::Variable) = Monomial{tuple()}()
copy(v::Variable) = v

isless(v1::Variable{N1}, v2::Variable{N2}) where {N1, N2} = isless(N1, N2)
show(io::IO, v::Variable{Name}) where Name = print(io, Name)

immutable Power{Var}
    exponent::Rational{Int}
end

variable(::Power{V}) where {V} = V
exponent(p::Power) = p.exponent
isless(p1::Power{V1}, p2::Power{V2}) where {V1, V2} = (V1, p1.exponent) < (V2, p2.exponent)

function show(io::IO, p::Power)
    print(io, variable(p))
    d = denominator(exponent(p))
    n = numerator(exponent(p))
    if d == 1
        if n != 1
            print(io, "^", n)
        end
    else
        print(io, "^(", exponent(p), ")")
    end
end

immutable Monomial{Powers} <: MonomialLike
end

@generated function convert(::Type{<:Monomial}, v::V) where {V <: Variable}
    quote
        Monomial{$((Power{V()}(1),))}()
    end
end

@generated function variables(::Type{Monomial{P}}) where {P}
    :($(variable.(P)))
end

variables(m::Monomial) = variables(typeof(m))

function show(io::IO, m::Monomial{P}) where {P}
    for power in P
        show(io, power)
    end
end

immutable Term{Mono, T} <: TermLike
    coefficient::T
    Term{Mono}(c::T) where{Mono, T} = new{Mono, T}(c)
    Term{Mono, T}(c::T) where{Mono, T} = new{Mono, T}(c)
end

convert(::Type{<:Term}, mono::M) where {M <: Monomial} = Term{mono}(1)
convert(T::Type{<:Term}, v::Variable) = convert(T, Monomial(v))
convert(::Type{<:Term}, x::Number) = Term{Monomial{tuple()}()}(x)

function show(io::IO, t::Term{Mono, T}) where {Mono, T}
    print(io, t.coefficient, Mono)
end

immutable Polynomial{Terms <: Tuple} <: PolynomialLike
    terms::Terms
end

convert(::Type{<:Polynomial}, t::Term) = Polynomial((t,))
convert(T::Type{<:Polynomial}, m::MonomialLike) = convert(T, convert(Term, m))

function show(io::IO, p::Polynomial)
    if length(p.terms) == 0
        print(io, "0")
    else
        print(io, p.terms[1])
        for t in p.terms[2:end]
            print(io, " + ", t)
        end
    end
end

@generated function literal_pow(^, v::V, ::Type{Val{x}}) where {V <: Variable, x}
    quote
        Monomial{$((Power{V()}(x),))}()
    end
end

@generated function (*)(v1::V, v2::V) where {V <: Variable}
    quote
        Monomial{$((Power{V()}(2),))}()
    end
end

@generated function (*)(v1::V1, v2::V2) where {V1 <: Variable, V2 <: Variable}
    @assert V1 != V2
    p = Tuple(sort([Power{V1()}(1), Power{V2()}(1)]))
    quote
        Monomial{$(p)}()
    end
end

# @generated function promote_rule(::Type{Term{Monomial{P1}, T}}, ::Type{Term{Monomial{P2}, T}}) where {T, P1, P2}
#     vars = Set{Any}()
#     for v in P1
#         push!(vars, v)
#     end
#     for v in P2
#         push!(vars, v)
#     end
#     vars = Tuple(sort(collect(vars)))
#     quote
#         Term{Monomial{$(vars)}(), T}
#     end
# end
#
powers(::Monomial{P}) where {P} = P
powers(::Type{Term{M, T}}) where {M, T} = powers(M)

function add_impl(::Type{Polynomial{Terms1}}, ::Type{Polynomial{Terms2}}) where {Terms1, Terms2}
    # TODO: accidentally quadratic
    p1 = powers.(collect(Terms1.parameters))
    p2 = powers.(collect(Terms2.parameters))
    allpowers = sort(collect(Set([p1..., p2...])))
    args = []
    for p in allpowers
        i1 = indexin([p], p1)[1]
        i2 = indexin([p], p2)[1]
        if i1 != 0 && i2 != 0
            push!(args, :(p1.terms[$i1] + p2.terms[$i2]))
        elseif i2 != 0
            push!(args, :(p2.terms[$i2]))
        else
            @assert i1 != 0
            push!(args, :(p1.terms[$i1]))
        end
    end
    args
end

@generated function +(p1::Polynomial{Terms1}, p2::Polynomial{Terms2}) where {Terms1, Terms2}
    Expr(:call, :Polynomial, Expr(:tuple, add_impl(p1, p2)...))
end

@generated function *(t1::Term{M1, T1}, t2::Term{M2, T2}) where {M1, M2, T1, T2}
    p = Any[powers(M1)...]
    for p2 in powers(M2)
        matched = false
        for (i, p1) in enumerate(powers(M1))
            if variable(p2) == variable(p1)
                p[i] = Power{variable(p[i])}(exponent(p[i]) + exponent(p2))
                matched = true
                break
            end
        end
        if !matched
            push!(p, p2)
        end
    end
    quote
        Term{Monomial{$(Tuple(sort(p)))}()}(t1.coefficient * t2.coefficient)
    end
end

(+)(t1::PolynomialLike, t2::PolynomialLike) = convert(Polynomial, t1) + convert(Polynomial, t2)
(+)(t1::Term{P, T1}, t2::Term{P, T2}) where {P, T1, T2} = Term{P}(t1.coefficient + t2.coefficient)
(+)(v1::MonomialLike, v2::MonomialLike) = Term(v1) + Term(v2)
(+)(m::TermLike, x::Number) = m + Term(x)
(+)(x::Number, m::TermLike) = Term(x) + m
(*)(x::Number, t::T) where {T <: Term} = T(x * t.coefficient)
(*)(t::T, x::Number) where {T <: Term} = T(t.coefficient * x)
(*)(x::Number, m::MonomialLike) = x * Term(m)
(*)(m::MonomialLike, x::Number) = Term(m) * x
(*)(t1::TermLike, t2::TermLike) = Term(t1) * Term(t2)
(*)(x::Number, m::Monomial) = x * Term(m)

@generated function subs(m::Term{Mono, T}, s::Pair{Variable{N}, T2}) where {Mono, T, N, T2}
    coeff = :(m.coefficient)
    remaining = []
    for p in powers(Mono)
        if variable(p) == Variable{N}()
            coeff = Expr(:call, :*, coeff, :(s.second ^ $(exponent(p))))
        else
            push!(remaining, p)
        end
    end
    if isempty(remaining)
        coeff
    else
        quote
            Term{$(Monomial{Tuple(remaining)}())}($coeff)
        end
    end
end

@generated function subs(m::Term, s::Vararg{<: Pair, N}) where {N}
    expr = :(m)
    for j in 1:N
        expr = Expr(:call, :subs, expr, :(s[$j]))
    end
    expr
end

@generated function subs(p::Polynomial{Terms}, s::Vararg{<:Pair, N}) where {N, Terms}
    Expr(:call, :+, [:(subs(p.terms[$i], s...)) for i in 1:length(Terms.parameters)]...)
end

subs(x::Number, s...) = x

(m::Polynomial)(args...) = subs(m, args...)
(m::Term)(args...) = subs(m, args...)

end


text:
#-------------------------------------------------------------------
#* EMSO Model Library (EML) Copyright (C) 2004 - 2007 ALSOC.
#*
#* This LIBRARY is free software; you can distribute it and/or modify
#* it under the therms of the ALSOC FREE LICENSE as available at
#* http://www.enq.ufrgs.br/alsoc.
#*
#* EMSO Copyright (C) 2004 - 2007 ALSOC, original code
#* from http://www.rps.eng.br Copyright (C) 2002-2004.
#* All rights reserved.
#*
#* EMSO is distributed under the therms of the ALSOC LICENSE as
#* available at http://www.enq.ufrgs.br/alsoc.
#*
#*----------------------------------------------------------------------
#* Author: Maur좩cio Carvalho Maciel
#* $Id$
#*--------------------------------------------------------------------
type Diff_Dist
	Diff_Dist()=begin
		PP=outers.PP
		NComp=outers.NComp
		new(
			DanaPlugin(Dict{Symbol,Any}(
				:Brief=>"External Physical Properties",
				:Type=>"PP"
			)),
			DanaInteger(Dict{Symbol,Any}(
				:Brief=>"Number of chemical components",
				:Lower=>1
			)),
			area(Dict{Symbol,Any}(
				:Brief=>"Cross Section Area"
			)),
			volume(Dict{Symbol,Any}(
				:Brief=>"Total volume"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Feed stream",
				:PosX=>0,
				:PosY=>0.9385,
				:Symbol=>"_{in}"
			)),
			stream(Dict{Symbol,Any}(
				:Brief=>"Liquid inlet stream",
				:PosX=>0.5,
				:PosY=>0.1984,
				:Symbol=>"_{inL}"
			)),
			vapour_stream(Dict{Symbol,Any}(
				:Brief=>"Vapour outlet stream",
				:PosX=>1,
				:PosY=>0.1984,
				:Symbol=>"_{outV}"
			)),
			power(Dict{Symbol,Any}(
				:Brief=>"Heat supplied",
				:PosX=>1,
				:PosY=>0.9578,
				:Symbol=>"_{in}"
			)),
			fill(mol(Dict{Symbol,Any}(
				:Brief=>"Molar Holdup in the distillator"
			)),(NComp)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar liquid holdup"
			)),
			mol(Dict{Symbol,Any}(
				:Brief=>"Molar vapour holdup"
			)),
			energy(Dict{Symbol,Any}(
				:Brief=>"Total Energy holdup on distillator"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Liquid Molar Volume"
			)),
			volume_mol(Dict{Symbol,Any}(
				:Brief=>"Vapour Molar volume"
			)),
			length(Dict{Symbol,Any}(
				:Brief=>"Level of liquid phase",
				:Default=>1,
				:Lower=>0
			)),
			temperature(Dict{Symbol,Any}(
				:Brief=>"Temperature on distillator"
			)),
			pressure(Dict{Symbol,Any}(
				:Brief=>"Pressure on distillator"
			)),
			fill(fraction(Dict{Symbol,Any}(
				:Brief=>"Molar Fraction of the Liquid of the distillator"
			)),(NComp)),
			enth_mol(Dict{Symbol,Any}(
				:Brief=>"Molar Enthalpy of the liquid of the distillator"
			)),
			[
				:(diff(M)= Inlet.F*Inlet.z + InletLiquid.F*InletLiquid.z - OutletVapour.F*OutletVapour.z),
				:(diff(E) = Inlet.F*Inlet.h + InletLiquid.F*InletLiquid.h - OutletVapour.F*OutletVapour.h + InletQ),
				:(M = ML*x + MV*OutletVapour.z),
				:(E = ML*h + MV*OutletVapour.h - P*V),
				:(sum(x)=1.0),
				:(sum(x)=sum(OutletVapour.z)),
				:(volL = PP.LiquidVolume(T, P, x)),
				:(volV = PP.VapourVolume(OutletVapour.T, OutletVapour.P, OutletVapour.z)),
				:(PP.LiquidFugacityCoefficient(T, P, x)*x = PP.VapourFugacityCoefficient(OutletVapour.T, OutletVapour.P, OutletVapour.z)*OutletVapour.z),
				:(P = OutletVapour.P),
				:(T = OutletVapour.T),
				:(V = ML*volL + MV*volV),
				:(Level = ML*volL/Across),
				:(h = PP.LiquidEnthalpy(T, P, x)),
			],
			[
				"Component Molar Balance","Energy Balance","Molar Holdup","Energy Holdup","Mol fraction normalisation","","Liquid Volume","Vapour Volume","Chemical Equilibrium","Mechanical Equilibrium","Thermal Equilibrium","Geometry Constraint","Level of liquid phase","Enthalpy",
			],
			[:PP,:NComp,:Across,:V,],
			[:Inlet,:InletLiquid,:OutletVapour,:InletQ,:M,:ML,:MV,:E,:volL,:volV,:Level,:T,:P,:x,:h,]
		)
	end
	PP::DanaPlugin
	NComp::DanaInteger
	Across::area
	V::volume
	Inlet::stream
	InletLiquid::stream
	OutletVapour::vapour_stream
	InletQ::power
	M::Array{mol}
	ML::mol
	MV::mol
	E::energy
	volL::volume_mol
	volV::volume_mol
	Level::length
	T::temperature
	P::pressure
	x::Array{fraction}
	h::enth_mol
	equations::Array{Expr,1}
	equationNames::Array{String,1}
	parameters::Array{Symbol,1}
	variables::Array{Symbol,1}
	attributes::Dict{Symbol,Any}
end
export Diff_Dist
function setEquationFlow(in::Diff_Dist)
	addEquation(1)
	addEquation(2)
	addEquation(3)
	addEquation(4)
	addEquation(5)
	addEquation(6)
	addEquation(7)
	addEquation(8)
	addEquation(9)
	addEquation(10)
	addEquation(11)
	addEquation(12)
	addEquation(13)
	addEquation(14)
end
function atributes(in::Diff_Dist,_::Dict{Symbol,Any})
	fields::Dict{Symbol,Any}=Dict{Symbol,Any}()
	fields[:Pallete]=true
	fields[:Icon]="icon/BatchDist"
	fields[:Brief]="Model of a Batch Differential Distillation."
	fields[:Info]="== Assumptions ==
* perfect mixing of both phases;
* thermodynamics equilibrium;
* no liquid entrainment in the vapour stream.

== Specify ==
* the inlet stream;
* the liquid inlet stream;
* the molar flow of the vapour outlet stream.

== Initial Conditions ==
* the distillator temperature (T);
* the distillator liquid level (Level);
* (NoComps - 1) compositions in the distillator or in the OutletV.
"
	drive!(fields,_)
	return fields
end
Diff_Dist(_::Dict{Symbol,Any})=begin
	newModel=Diff_Dist()
	newModel.attributes=atributes(newModel,_)
	newModel
end


text:
load("kernels.jl")

dest1 = Array(Int, (2,3))
dest2 = Array(Int, (2,3))
B = [1 2 3
     4 5 6]
C = [1 0 1
     0 1 0]
D = [ 0  0  0
     10 10 10]

# code = quote
#     A = B.*C + D
#     dest1[...] = A
#     dest2[...] = A + C
# end
# ex = code_kernel(code)
# eval(ex)

@kernel let
    A = B.*C + D
    dest1[...] = A
    dest2[...] = A + C
end

function f(B,C,D)
    A = B.*C + D
    dest1 = A
    dest2 = A + C
    (dest1, dest2)
end
dest1f, dest2f = f(B,C,D)

println("kernel:     dest1 =\n$dest1")
println("Array ops:  dest1 =\n$dest1f")
println()
println("kernel:     dest2 =\n$dest2")
println("Array ops:  dest2 =\n$dest2f")
println()

max_stars_repo_name: toivoh/julia-kernels

max_issues_repo_name: toivoh/julia-kernels

max_forks_repo_name: toivoh/julia-kernels


text:
using ArgParse

function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin
        "SEED"
            help = "Seed for randomness"
            arg_type = Int
            default = 1
        "num_paths"
            help = "Number of price paths to generate"
            arg_type = Int
            default = 1
        "T"
            help = "Number of time periods"
            arg_type = Int
            default = 100
        "p₀"
            arg_type = Float64
            default = 100.0
        "M"
            help = "Number of price points"
            arg_type = Int
            default = 100
        "L"
            arg_type = Float64
            default = 50.0
        "D"
            arg_type = Float64
            default = 4.0
        "σ"
            arg_type = Float64
            default = 1.0
        "nu"
            arg_type = Float64
            default = 0.1
        "α"
            arg_type = Float64
            default = 20.0
        "λ"
            arg_type = Float64
            default = 1.0
        "μ"
            arg_type = Float64
            default = 0.5
    end

    return parse_args(s)
end


text:
module NetworkProtocols

using Printf
using Sockets
using UnsafeArrays

export MACAddress
export isunicast, ismulticast, isuniversal, islocal

export EthernetHeader, EthernetPacket
export ETHERTYPE_IPV4, ETHERTYPE_ARP
export decode_ethernet, ethertype_string

export IPv4Header, IPv4Packet
export IPPROTOCOL_IGMP, IPPROTOCOL_TCP, IPPROTOCOL_UDP
export decode_ipv4

export TCPHeader, TCPPacket
export UDPHeader, UDPPacket
export decode_tcp, decode_udp

export dispatch_ethernet

include("mac_address.jl")
include("ethernet.jl")
include("ip.jl")
include("tcp.jl")
include("udp.jl")
include("dispatch.jl")

module Experimental
    import FunctionWrappers.FunctionWrapper
    using Sockets
    using ..NetworkProtocols

    export TCPStream, tcpstreamid
    export TCPEvent, TCPE_OPEN, TCPE_CONNECTED, TCPE_CLOSE, TCPE_DISCONNECTED, TCPE_ERROR, TCPE_RETRANSMIT

    include("experimental/tcp_stream.jl")
end

end

max_stars_repo_name: ancapdev/NetworkProtocols

max_issues_repo_name: ancapdev/NetworkProtocols

max_forks_repo_name: ancapdev/NetworkProtocols


text:
using DelimitedFiles
include("GB11.jl")
import GLM: fit, predict, LinearModel
using Distributions,SpecialFunctions
using LinearAlgebra, MultivariateStats
using StatsBase, Random, Distributed
using DelimitedFiles, DataStructures
using CSV
using DataFrames

rng = Random.seed!(2019)

gbdt = GB11.GBDT(;
  sampling_rate = 0.75,                # Sampling rate
  learning_rate = 0.1,                # Learning rate
  num_iterations = 100,               # NuGBer of iterations
  tree_options = Dict(
      :nsubfeatures => 0,
      :min_samples_leaf => 1,
      :max_depth => 5))
gbl = GB11.GBLearner(gbdt, :regression)

######## σ² = 1 #########
temp_storedim = size(collect(10:1000),1)
fitval = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rssfitted = []
predicted = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted = []
Ytestval = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
μ_testval = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted2 = []
HyperS = Array{Union{Missing, Nothing, Int, Float64}, 2}(nothing, temp_storedim, 1000)
s
for pval in collect(10:11)
    # reproducibility
    Random.seed!(2019)
    rng = 2019

    ### simulate the training and testing data
    X, Y = GB11.test_data(500; p = pval, sigma = 1.0)
    X_test, Y_test = GB11.test_data(500; p = pval)
    μ = GB11.fried(X)
    μ_test = GB11.fried(X_test)

    println("The last simulated data:
        \nX:$(X[500, pval]) \nY:$(Y[500])
        \nX_test: $(X_test[500, pval]) \nY_test: $(Y_test[500])
        \nμ: $(μ[500]) \nμ_test: $(μ_test[500])\n\n")

    μ_testval[:, pval - 9] = μ_test
    Ytestval[:, pval - 9 ] = Y_test

    hypers = GB11.Hypers(X, Y, 1, 1, 0.95)
    hypers.σ = sqrt(1)  ### TODO: I am using
    hypers.τ = 3 / 10 * var(Y) / gbdt.num_iterations
    hypers.α = 0.95
    hypers.s =  hypers.α .* ones(size(X, 2)) / size(X, 2)

    # train the model
    mymodel=GB11.fit!(gbl, X, Y, hypers)
    # get the fitted value
    myfit  = GB11.predict!(gbl, X)

    # store the hyper.s for the ensembel tree
    HyperS[pval - 9, 1:size(hypers.s, 1)] .= hypers.s

    # store the fitted value
    fitval[:, pval - 9] = myfit

    # get the residual sum of squares
    rssfit = sqrt(mean(Y .- myfit) .^2)
    push!(rssfitted, rssfit)
    # get prediction
    mypred = GB11.predict!(gbl, X_test)

    println("The last pred is $(mypred[500])\n\n ")

    predicted[:, pval - 9] = mypred

    # get mse and root of mean squared error
    rmse = sqrt(mean((Y_test .- mypred ) .^2))

    println("The rmse b.w. Y_test and pred for p = $(pval) is $(rmse) ")

    push!(rmsefitted, rmse)

    # get the rmse(μ_test, mypred)
    rmse2 = sqrt(mean((μ_test .- mypred).^2))

    println("The rmse2 b.w. μ_test and pred for p = $(pval)  is $(rmse2).")

    push!(rmsefitted2, rmse2)
end

rssfitted
rmsefitted
rmsefitted2
Ytestval
μ_testval
HyperS

HyperS_copy = copy(HyperS)
for i in 1 : size(HyperS, 1)
    for j in 1 : size(HyperS, 2)
        if HyperS_copy[i, j] == nothing
            HyperS_copy[i, j] = 0.0
        end
    end
end
HyperS_copy

fitval = DataFrame(fitval)
CSV.write("_fitval_sigma1.csv", fitval)
rssfitted = DataFrame(rssfitted = rssfitted)
CSV.write("_fss_sigma1.csv", rssfitted)
rmsefitted = DataFrame(rmsefitted = rmsefitted)
CSV.write("_rmse_sigma1.csv", rmsefitted)
predicted = DataFrame(predicted)
CSV.write("_pred_sigma1.csv", predicted)
rmsefitted2 = DataFrame(rmsefitted2 = rmsefitted2)
CSV.write("_rmse2_sigma1.csv", rmsefitted2)
Ytestval = DataFrame(Ytestval)
CSV.write("_Ytest_sigma1.csv", Ytestval)
μ_testval = DataFrame(μ_testval)
CSV.write("_mu_test_sigma1.csv", μ_testval)
HyperS = DataFrame(HyperS)
HyperS_copy = DataFrame(HyperS_copy)
CSV.write("_HyperS_copy_sigma1.csv", HyperS_copy)
##### not runned yet

######## σ² = 10 #########
temp_storedim = size(collect(10:1000),1)
fitval_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rssfitted_sigma10 = []
predicted_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted_sigma10 = []
Ytestval_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
μ_testval_sigma10 = Array{Union{Missing, Nothing, Int, Float64}}(nothing, 500, temp_storedim)
rmsefitted2_sigma10 = []
HyperS_sigma10 = Array{Union{Missing, Nothing, Int, Float64}, 2}(nothing, temp_storedim, 1000)
s
for pval in collect(10:1000)
    # reproducibility
    Random.seed!(2019)
    rng = 2019

    ### simulate the training and testing data
    X, Y = GB11.test_data(500; p = pval, sigma = sqrt(10.0))
    X_test, Y_test = GB11.test_data(500; p = pval, sigma = sqrt(10.0))
    μ = GB11.fried(X)
    μ_test = GB11.fried(X_test)

    println("Simulated data:
    \nX: $(X[500, pval])
    \nY: $(Y[500])
    \nX_test: $(X_test[500, pval])
    \nY_test: $(Y_test[500])
    \nμ: $(μ[500])
    \nμ_test: $(μ_test[500])")

    μ_testval_sigma10[:, pval - 9] = μ_test
    Ytestval_sigma10[:, pval - 9] = Y_test

    hypers = GB11.Hypers(X, Y, 1, 1, 0.95)
    hypers.σ = sqrt(1)  # TODO: maybe change to 10 after checking the plot
    hypers.τ = 3 / 10 * var(Y) / gbdt.num_iterations
    hypers.α = 0.95
    hypers.s =  hypers.α .* ones(size(X, 2)) / size(X, 2)

    # train the model
    mymodel=GB11.fit!(gbl, X, Y, hypers)

    # store the hyper.s for the ensembel tree
    HyperS_sigma10[ pval - 9, 1:size(hypers.s, 1)] .= hypers.s

    # get the fitted value
    myfit  = GB11.predict!(gbl, X)
    # store the fitted value
    fitval_sigma10[:, pval - 9] = myfit

    # get the residual sum of squares
    rssfit = sqrt(mean((Y .- myfit) .^2))
    push!(rssfitted_sigma10, rssfit)
    # get prediction
    mypred = GB11.predict!(gbl, X_test)

    println("The last pred is $(mypred[500])")

    predicted_sigma10[:, pval - 9] = mypred

    # get mse and root of mean squared error
    rmse = sqrt(mean((Y_test .- mypred ) .^2))
    println("The p = $(pval), and the rmse is $rmse.")
    push!(rmsefitted_sigma10, rmse)

    # get the rmse(Y_test, μ_test)
    rmse2 = sqrt(mean((mypred .- μ_test) .^2))
    println("The p = $(pval) dim, rmse2 is $rmse2.")
    push!(rmsefitted2_sigma10, rmse2)
end

fitval_sigma10
rssfitted_sigma10
predicted_sigma10
rmsefitted_sigma10
Ytestval_sigma10
μ_testval_sigma10
rmsefitted2_sigma10
HyperS_sigma10

#predicted_sigma10 = DataFrame(predicted_sigma10)
rmsefitted_sigma10 = DataFrame(rmsefitted_sigma10=rmsefitted_sigma10)

#CSV.write("fitval_sigma10.csv", fitval_sigma10)
#CSV.write("rssfitted_sigma10.csv", rssfitted_sigma10)
#CSV.write("predicted_sigma10.csv", predicted_sigma10)
#CSV.write("rmsefitted_sigma10.csv", rmsefitted_sigma10)

fitval_sigma10 = DataFrame(fitval_sigma10)
CSV.write("_fitval_sigma10.csv", fitval_sigma10)
rssfitted_sigma10 = DataFrame(rssfitted_sigma10 = rssfitted_sigma10)
CSV.write("_fss_sigma10.csv", rssfitted_sigma10)

CSV.write("_rmse_sigma10.csv", rmsefitted_sigma10)
predicted_sigma10 = DataFrame(predicted_sigma10)
CSV.write("_pred_sigma10.csv", predicted_sigma10)
rmsefitted2_sigma10 = DataFrame(rmsefitted2_sigma10 = rmsefitted2_sigma10)
CSV.write("_rmse2_sigma10.csv", rmsefitted2_sigma10)
Ytestval_sigma10 = DataFrame(Ytestval_sigma10)
CSV.write("_Ytest_sigma10.csv", Ytestval_sigma10)
μ_testval_sigma10 = DataFrame(μ_testval_sigma10)
CSV.write("_mu_test_sigma10.csv", μ_testval_sigma10)

HyperS_sigma10 = DataFrame(HyperS_sigma10)

HyperS_copy_sigma10 = copy(HyperS_sigma10)

for i in 1 : size(HyperS_sigma10, 1)

    for j in 1 : size(HyperS_sigma10, 2)

        if HyperS_copy_sigma10[i, j] == nothing

            HyperS_copy_sigma10[i, j] = 0.0

        end

    end

end

HyperS_copy_sigma10
CSV.write("_HyperS_sigma10.csv", HyperS_sigma10)

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
"""
    PyList{T=PyObject}([o])

Wrap the Python list `o` (or anything satisfying the sequence interface) as a Julia vector with elements of type `T`.

If `o` is not given, an empty list is created.
"""
struct PyList{T} <: AbstractVector{T}
    ref::PyRef
    PyList{T}(o) where {T} = new{T}(ispyref(o) ? PyRef(o) : pylist(PyRef, o))
    PyList{T}() where {T} = new{T}(PyRef())
end
PyList(o) = PyList{PyObject}(o)
PyList() = PyList{PyObject}()
export PyList

ispyreftype(::Type{<:PyList}) = true
pyptr(x::PyList) = begin
    ptr = x.ref.ptr
    if isnull(ptr)
        ptr = x.ref.ptr = C.PyList_New(0)
    end
    ptr
end
Base.unsafe_convert(::Type{CPyPtr}, x::PyList) = checknull(pyptr(x))
C.PyObject_TryConvert__initial(o, ::Type{T}) where {T<:PyList} =
    C.putresult(T(pyborrowedref(o)))

# Base.length(x::PyList) = @pyv `len($x)`::Int
Base.length(x::PyList) = Int(pylen(x))

Base.size(x::PyList) = (length(x),)

Base.getindex(x::PyList{T}, i::Integer) where {T} = begin
    checkbounds(x, i)
    # The following line implements this function, but is typically 3x slower.
    # @pyv `$x[$(i-1)]`::T
    p = checknull(C.PySequence_GetItem(x, i - 1))
    r = C.PyObject_Convert(p, T)
    C.Py_DecRef(p)
    checkm1(r)
    C.takeresult(T)
end

Base.setindex!(x::PyList{T}, _v, i::Integer) where {T} = begin
    checkbounds(x, i)
    v = convertref(T, _v)
    # The following line implements this function, but is typically 10x slower
    # @py `$x[$(i-1)] = $v`
    vp = checknull(C.PyObject_From(v))
    err = C.PySequence_SetItem(x, i - 1, vp)
    C.Py_DecRef(vp)
    checkm1(err)
    x
end

Base.insert!(x::PyList{T}, i::Integer, v) where {T} = (
    i == length(x) + 1 || checkbounds(x, i); @py `$x.insert($(i-1), $(convertref(T, v)))`; x
)

Base.push!(x::PyList{T}, v) where {T} = (@py `$x.append($(convertref(T, v)))`; x)

Base.pushfirst!(x::PyList, v) = insert!(x, 1, v)

Base.pop!(x::PyList{T}) where {T} = @pyv `$x.pop()`::T

Base.delete!(x::PyList{T}, i::Integer) where {T} =
    (checkbounds(x, i); @py `$x.pop($(i-1))`; x)

Base.popfirst!(x::PyList) = pop!(x, 1)

Base.reverse!(x::PyList) = (@py `$x.reverse()`; x)

# TODO: support kwarg `by` (becomes python kwarg `key`)
# Base.sort!(x::PyList; rev::Bool=false) = (@py `$x.sort(reverse=$rev)`; x)

Base.empty!(x::PyList) = (@py `$x.clear()`; x)

Base.copy(x::PyList) = @pyv `$x.copy()`::typeof(x)


text:
### ScRNAseq.jl
###
### A julia package for scRNA-seq exploratory data analysis
##
### This file is NOT YET a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

module ScRNAseq

export
    FileIO,
    QualityControl,
    Transformation,
    Normalization,
    Embedding,
    Clustering,
    DifferentialExpression,
    GeneRegulatoryNetwork

include("helper.jl")
include("fileio/fileio.jl")
include("qualitycontrol/qualitycontrol.jl")
include("transformation/transformation.jl")
include("normalization/normalization.jl")
include("embedding/embedding.jl")
include("clustering/clustering.jl")
include("differentialexpression/differentialexpression.jl")
include("generegulatorynetwork/generegulatorynetwork.jl")

using .FileIO
using .QualityControl
using .Transformation
using .Normalization
using .Embedding
using .Clustering
using .DifferentialExpression
using .GeneRegulatoryNetwork

end  # module ScRNAseq


text:
#   Copyright (c) 2020-2022 Velexi Corporation
#
#   Licensed under the Apache License, Version 2.0 (the "License");
#   you may not use this file except in compliance with the License.
#   You may obtain a copy of the License at
#
#       http://www.apache.org/licenses/LICENSE-2.0
#
#   Unless required by applicable law or agreed to in writing, software
#   distributed under the License is distributed on an "AS IS" BASIS,
#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
#   See the License for the specific language governing permissions and
#   limitations under the License.

"""
add.jl defines methods for the +(x, y) function
"""

# --- Exports

import Base.:(+)

# --- Method definitions

# ------ Specializations involving an AbstractMultivector instance

# M::AbstractMultivector, B::One
# B::One,M::AbstractMultivector
@inline +(M::AbstractMultivector, B::One) = Multivector(vcat([B], blades(M)))
@inline +(B::One, M::AbstractMultivector) = B + M

# M::AbstractMultivector, B::Zero
# B::Zero, M::AbstractMultivector
@inline +(M::AbstractMultivector, B::Zero) = M
@inline +(B::Zero, M::AbstractMultivector) = M

# ------ Specializations involving a Blade instance

# B::Blade, C::Blade
function +(B::Blade, C::Blade)
    assert_dim_equal(B, C)

    if grade(B) == grade(C) == 1
        return Blade(volume(B) * basis(B) + volume(C) * basis(C))
    end

    # TODO: general blade sum
end

# ------ Specializations involving a Pseudoscalar instance

# B::Pseudoscalar, C::Pseudoscalar
@inline +(B::Pseudoscalar, C::Pseudoscalar) =
    Pseudoscalar(B, value=value(B) + value(C))

# ------ Specializations involving an AbstractScalar instance

# B::AbstractScalar, C::AbstractScalar
@inline +(B::AbstractScalar, C::AbstractScalar) =
    Scalar{typeof(value(B))}(value(B) + value(C))

# B::AbstractScalar, C::One
# B::One, C::AbstractScalar
@inline +(B::AbstractScalar, C::One) = Scalar{typeof(value(B))}(value(B) + 1)
@inline +(B::One, C::AbstractScalar) = C + B

# B::AbstractScalar, C::Zero
# B::Zero, C::AbstractScalar
@inline +(B::AbstractScalar, C::Zero) = B
@inline +(B::Zero, C::AbstractScalar) = C

# B::AbstractScalar, x::Real
# x::Real, B::AbstractScalar
@inline +(B::AbstractScalar, x::Real) = Scalar{typeof(value(B))}(value(B) + x)
@inline +(x::Real, B::AbstractScalar) = B + x

# ------ Specializations involving a One instance

# B::One, C::One
@inline +(B::One, C::One) = Scalar{typeof(value(B))}(2)

# B::One, C::Zero
# B::Zero, C::One
@inline +(B::One, C::Zero) = B
@inline +(B::Zero, C::One) = C

# ------ Specializations involving a Zero instance

# B::Zero, C::Zero
@inline +(B::Zero, C::Zero) = B

# B::Zero, v::Vector
# v::Vector, B::Zero
@inline +(B::Zero, v::Vector{<:Real}) = v
@inline +(v::Vector{<:Real}, B::Zero) = v


text:
isvalid(sides) = all(sides .> 0) && sides[1] + sides[2] >= sides[3]

function is_equilateral(sides)
    sort!(sides)
    isvalid(sides) && (sides[1] == sides[2] == sides[3])
end

function is_isosceles(sides)
    sort!(sides)
    isvalid(sides) && (sides[1] == [2] || sides[2] == sides[3])
end

function is_scalene(sides)
    sort!(sides)
    isvalid(sides) && (sides[1] != sides[2] != sides[3])
end

max_stars_repo_name: twolodzko/Learning

max_issues_repo_name: twolodzko/Learning

max_forks_repo_name: twolodzko/Learning


text:
# ---
# title: 377. Combination Sum IV
# id: problem377
# author: Tian Jun
# date: 2020-10-31
# difficulty: Medium
# categories: Dynamic Programming
# link: <https://leetcode.com/problems/combination-sum-iv/description/>
# hidden: true
# ---
#
# Given an integer array with all positive numbers and no duplicates, find the
# number of possible combinations that add up to a positive integer target.
#
# **Example:**
#
#
#
#     _**nums**_ = [1, 2, 3]
#     _**target**_ = 4
#
#     The possible combination ways are:
#     (1, 1, 1, 1)
#     (1, 1, 2)
#     (1, 2, 1)
#     (1, 3)
#     (2, 1, 1)
#     (2, 2)
#     (3, 1)
#
#     Note that different sequences are counted as different combinations.
#
#     Therefore the output is _**7**_.
#
#
#
#
# **Follow up:**
# What if negative numbers are allowed in the given array?
# How does it change the problem?
# What limitation we need to add to the question to allow negative numbers?
#
# **Credits:**
# Special thanks to [@pbrother](https://leetcode.com/pbrother/) for adding this
# problem and creating all test cases.
#
#
## @lc code=start
using LeetCode

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


text:
# This file is a part of JuliaFEM.
# License is MIT: see https://github.com/JuliaFEM/JuliaFEM.jl/blob/master/LICENSE.md

using JuliaFEM
using JuliaFEM.Preprocess
using JuliaFEM.Postprocess
using JuliaFEM.Testing

@testset "3d rod" begin
    mesh = aster_read_mesh(Pkg.dir("JuliaFEM")*"/test/testdata/primitives.med", "CYLINDER_20_TET4")
    problem = Problem(Heat, "rod of length 20", 1)
    problem.elements = create_elements(mesh, "CYLINDER")
    update!(problem, "temperature thermal conductivity", 200.0)
    outer = Problem(Heat, "outer surface", 1)
    outer.elements = create_elements(mesh, "FACE2", "FACE3")
    update!(outer, "temperature external temperature", 20.0)
    update!(outer, "temperature heat transfer coefficient", 1.0)
    #midline = Problem(Heat, "midline of rod", 1)
    #midline.elements = create_elements(mesh, "INNER_LINE")
    boundary = Problem(Dirichlet, "homogeneous dirichlet boundary", 1, "temperature")
    boundary.elements = create_elements(mesh, "FACE1")
    update!(boundary, "temperature 1", 100.0)
    #solver = LinearSolver(problem, outer, boundary, midline)
    solver = LinearSolver(problem, outer, boundary)
    solver()

    L = 20
    k = 200.0
    Tu = 20.0
    h = 1.0
    P = 2*pi
    A = pi
    α = h
    β = sqrt((h*P)/(k*A))
    T0 = 100.0
    C = [1.0 1.0; (α+k*β)*exp(β*L) (α-k*β)*exp(-β*L)] \ [T0-Tu, 0]
    T(x) = dot(C, [exp(β*x), exp(-β*x)]) + Tu

    T_diff = []
    for x in linspace(0, 20)
        T_FEM = problem("temperature", [x, 0.0, 0.0])[1]
        T_ACC = T(x)
        push!(T_diff, norm(T_FEM - T_ACC))
        info("x = $x, T_FEM = $T_FEM, T_ACC = $T_ACC")
    end
    info("mean diff = ", mean(T_diff))
    # mean diff = 1.14
    @test mean(T_diff) < 1.2
end


text:
"""
    PushBot
"""
mutable struct PushBot{T} <: Model{T}
    nq::Int
	nu::Int
	nw::Int
	nc::Int

    mb::T # mass
	ma::T
	l::T # length

    μ_world::T  # coefficient of friction
    μ_joint::T  # gravity
	g::T

	base::BaseMethods
	dyn::DynamicsMethods

	joint_friction::SVector
end

# Kinematics
function _kinematics(model::PushBot, q; mode = :com)
	θ, d  = q
	if mode == :d
		return [-model.l * sin(θ) + d * cos(θ);
				 model.l * cos(θ) + d * sin(θ)]
	elseif mode == :com
		return [-1.0 * model.l * sin(θ);
				 1.0 * model.l * cos(θ)]
	elseif mode == :ee
		return [-model.l * sin(θ);
				 model.l * cos(θ)]
	else
		@error "incorrect mode"
	 	return zeros(2)
	end
end

function _jacobian(model::PushBot, q; mode = :com)
	θ, d = q

	if mode == :d
 		return [(-model.l * cos(θ) - d * sin(θ)) cos(θ);
		        (-model.l * sin(θ) + d * cos(θ)) sin(θ)]
 	elseif mode == :com
 		return [-1.0 * model.l * cos(θ) 0.0;
		        -1.0 * model.l * sin(θ) 0.0]
	elseif mode == :ee
		return [-model.l * cos(θ) 0.0;
		        -model.l * sin(θ) 0.0]
 	else
 		@error "incorrect mode"
 	 	return
 	end
end

function kinematics(model::PushBot, q)
	[_kinematics(model, q, mode = :d);
	 _kinematics(model, q, mode = :d)]
end

function lagrangian(model::PushBot, q, q̇)
	L = 0.0

	vθ = _jacobian(model, q, mode = :com) * q̇
	L += 0.5 * model.mb * transpose(vθ) * vθ
	L -= model.mb * model.g * _kinematics(model, q, mode = :com)[2]

	vd1 = _jacobian(model, q, mode = :d) * q̇
	L += 0.5 * model.ma * transpose(vd1) * vd1
	L -= model.ma * model.g * _kinematics(model, q, mode = :d)[2]

	return L
end

function M_func(model::PushBot, q)
	Jθ = _jacobian(model, q, mode = :com)
	Jd = _jacobian(model, q, mode = :d)

	return model.mb * transpose(Jθ) * Jθ + model.ma * transpose(Jd) * Jd
end

function ϕ_func(model::PushBot, env::Environment, q)
	# walls at x = -0.5, x = 0.5
    SVector{2}([_kinematics(model, q, mode = :d)[1] + 0.5;
	            0.5 - _kinematics(model, q, mode = :d)[1]])
end

function rot2D(x)
	[cos(x) sin(x);
    -sin(x) cos(x)]
end

function J_func(model::PushBot, env::Environment, q)
	r1 = [0.0 -1.0; 1.0 0.0]
	r2 = [0.0 1.0; -1.0 0.0]

    SMatrix{4, 2}([r1 * _jacobian(model, q, mode = :d);
		           r2 * _jacobian(model, q, mode = :d)])
end

function B_func(model::PushBot, q)
	@SMatrix [1.0 * model.l 1.0;
	          1.0 1.0 / model.l]
end

function A_func(::PushBot, q)
	@SMatrix [1.0 0.0;
			  0.0 1.0]
end

# Parameters
g = 9.81 # gravity
μ_world = 0.5  # coefficient of friction
μ_joint = 10.0

mb = 1.0 # body mass
ma = 0.01 # mass of ee # CHANGED
l = 1.0  # leg mass

# Dimensions
nq = 2
nu = 2
nw = nq
nc = 2
nquat = 0

pushbot = PushBot(nq,nu,nw,nc,
			   mb, ma, l,
			   μ_world, μ_joint, g,
			   BaseMethods(), DynamicsMethods(),
			   SVector{2}(μ_joint * [1.0; 1.0]))

function friction_coefficients(model::PushBot)
	return [model.μ_world]
end

function initialize_z!(z, model::PushBot, idx::RoboDojo.IndicesZ, q)
    z .= 1.0
    z[idx.q] .= q
end


text:
using Test
using LinearAlgebra
using QuantumPropagators: initpropwrk, propstep!
using QuantumPropagators.Newton
using QuantumPropagators.SpectralRange: specrange
using QuantumControlBase:
    TimeDependentGradGenerator, GradGenerator, GradVector, resetgradvec!, evalcontrols
using QuantumControlBase.TestUtils
using Zygote

@testset "GradGenerator" begin

    N = 10  # size of Hilbert space
    ρ = 1.0  # spectral radius
    # We'll test with non-Hermitian Hamiltonians
    Ĥ₀ = random_complex_matrix(N, ρ)
    Ĥ₁ = random_complex_matrix(N, ρ)
    Ĥ₂ = random_complex_matrix(N, ρ)
    Zero = zeros(ComplexF64, N, N)
    ϵ₁ = t -> 1.0
    ϵ₂ = t -> 1.0
    Ĥ_of_t = (Ĥ₀, (Ĥ₁, ϵ₁), (Ĥ₂, ϵ₂))
    Ψ = random_state_vector(N)
    ψ_max = maximum(abs.(Ψ))
    Ψtgt = random_state_vector(N)
    𝕚 = 1im
    dt = 1.25

    vals_dict = IdDict(ϵ₁ => 1.0, ϵ₂ => 1.0)

    G̃_of_t = TimeDependentGradGenerator(Ĥ_of_t)
    G̃ = evalcontrols(G̃_of_t, vals_dict)
    Ĥ = evalcontrols(Ĥ_of_t, vals_dict)

    Û_Ψ = exp(-𝕚 * Ĥ * dt) * Ψ

    num_controls = length(Ĥ_of_t) - 1
    Ψ̃ = GradVector(Ψ, num_controls)
    # did the initialization work?
    @test norm(Ψ̃.state - Ψ) < 1e-14
    @test norm(Ψ̃.grad_states[1]) == 0.0
    @test norm(Ψ̃.grad_states[2]) == 0.0
    @test length(Ψ̃.grad_states) == 2

    @test maximum(abs.(Ψ)) == ψ_max  # is Ψ still exactly the same state?
    wrk = NewtonWrk(Ψ̃)
    newton!(Ψ̃, G̃, dt, wrk)
    Ψ̃_out = copy(Ψ̃)
    @test maximum(abs.(Ψ)) == ψ_max  # is Ψ still exactly the same state?
    @test norm(Ψ̃_out.state - Û_Ψ) < 1e-12  # Ψ̃_out.state correct?
    @test norm(Ψ̃_out.grad_states[1]) > 0
    @test norm(Ψ̃_out.grad_states[2]) > 0
    # did the copy work?
    @test norm(Ψ̃_out.state - Ψ̃.state) == 0.0
    @test norm(Ψ̃_out.grad_states[1] - Ψ̃.grad_states[1]) == 0.0
    @test norm(Ψ̃_out.grad_states[2] - Ψ̃.grad_states[2]) == 0.0

    # reset Ψ̃ to the original input (also to test the resetgradvec! method)
    resetgradvec!(Ψ̃, Ψ)
    # state should be the same as after previous initialization
    @test norm(Ψ̃.state - Ψ) < 1e-14
    @test norm(Ψ̃.grad_states[1]) == 0.0
    @test norm(Ψ̃.grad_states[2]) == 0.0

    @test norm(Ψ̃_out.state - Û_Ψ) < 1e-12  # Ψ̃_out.state still correct?

    ###########################################################################
    # Compare against explicit Grad-Gen
    # This checks whether the application of a GradGenerator to a GradVector an
    # all the linear-algebra methods are implemented correctly

    G̃_full = [
         Ĥ    Zero  Ĥ₁
        Zero   Ĥ    Ĥ₂
        Zero  Zero  Ĥ
    ]

    Ψ̃_full = [
        Ψ̃.grad_states[1]
        Ψ̃.grad_states[2]
        Ψ̃.state
    ]

    # proper initialization? grad_states should be zero
    @test norm(Ψ̃_full) == norm(Ψ̃.state) == norm(Ψ)

    Ψ̃_out_full = exp(-𝕚 * G̃_full * dt) * Ψ̃_full
    # propagation correct?
    @test norm(Ψ̃_out_full[2N+1:3N] - Û_Ψ) < 1e-12

    # do we get the same results as from newton?
    @test norm(Ψ̃_out_full[2N+1:3N] - Ψ̃_out.state) < 1e-12
    @test norm(Ψ̃_out_full[1:N] - Ψ̃_out.grad_states[1]) < 1e-12
    @test norm(Ψ̃_out_full[N+1:2N] - Ψ̃_out.grad_states[2]) < 1e-12

    ###########################################################################
    # Test custom expprop

    Ψ̃_full = [
        Ψ̃.grad_states[1]
        Ψ̃.grad_states[2]
        Ψ̃.state
    ]

    # proper initialization? grad_states should be zero
    @test norm(Ψ̃_full) == norm(Ψ̃.state) == norm(Ψ)

    Ψ̃_out_full = exp(-𝕚 * G̃_full * dt) * Ψ̃_full
    # propagation correct?
    @test norm(Ψ̃_out_full[2N+1:3N] - Û_Ψ) < 1e-12

    # do we get the same results as from newton?
    @test norm(Ψ̃_out_full[2N+1:3N] - Ψ̃_out.state) < 1e-12
    @test norm(Ψ̃_out_full[1:N] - Ψ̃_out.grad_states[1]) < 1e-12
    @test norm(Ψ̃_out_full[N+1:2N] - Ψ̃_out.grad_states[2]) < 1e-12

    ###########################################################################
    # Test standard expprop

    wrk_exp = initpropwrk(Ψ̃, [0, dt], Val(:expprop), G̃)
    propstep!(Ψ̃, G̃, dt, wrk_exp)
    Ψ̃_out_exp = copy(Ψ̃)
    @test norm(wrk_exp.Ψ_full - Ψ̃_full) ≈ 0.0
    @test norm(wrk_exp.G_full - G̃_full) ≈ 0.0
    @test norm(Ψ̃_out_exp - Ψ̃_out) < 1e-11
    resetgradvec!(Ψ̃, Ψ)

    ###########################################################################
    # Compare against explicit split-up Grad-Gens
    # This checks whether the extension of the Gradient-Vector to multiple
    # controls is correct (the literature generally only gives the expression
    # for a single control)

    G̃_full1 = [
         Ĥ   Ĥ₁
        Zero Ĥ
    ]
    G̃_full2 = [
         Ĥ   Ĥ₂
        Zero Ĥ
    ]
    Ψ̃_full1 = [
        Ψ̃.grad_states[1]
        Ψ̃.state
    ]

    @test maximum(abs.(Ψ)) == ψ_max  # is Ψ still exactly the same state?
    @test norm(Ψ̃_full1) == norm(Ψ)  # initialization correct?
    Ψ̃_full2 = vcat(Ψ̃.grad_states[2], Ψ̃.state)
    @test norm(Ψ̃_full2) == norm(Ψ)  # initialization correct?
    Ψ̃_out_full1 = exp(-𝕚 * G̃_full1 * dt) * Ψ̃_full1
    Ψ̃_out_full2 = exp(-𝕚 * G̃_full2 * dt) * Ψ̃_full2

    # propagation correct?
    @test norm(Ψ̃_out_full1[N+1:2N] - Û_Ψ) < 1e-12
    @test norm(Ψ̃_out_full2[N+1:2N] - Û_Ψ) < 1e-12

    # do we get the same results as with the combined grad-gen?
    @test norm(Ψ̃_out_full1[1:N] - Ψ̃_out_full[1:N]) < 1e-12
    @test norm(Ψ̃_out_full2[1:N] - Ψ̃_out_full[N+1:2N]) < 1e-12
    @test norm(Ψ̃_out_full1[N+1:2N] - Ψ̃_out_full[2N+1:3N]) < 1e-12
    @test norm(Ψ̃_out_full2[N+1:2N] - Ψ̃_out_full[2N+1:3N]) < 1e-12

    ###########################################################################
    # Compare against Zygote
    # This checks whether the gradients are correct, albeit a bit indirectly
    # (since Zygote can only calculate the gradient for a scalar function, in
    # this case the square-modulus of the overlap with a target)

    function F_sm(ϵ₁, ϵ₂)
        Û = exp(-𝕚 * (Ĥ₀ + ϵ₁ * Ĥ₁ + ϵ₂ * Ĥ₂) * dt)
        Ψ_T = Û * Ψ
        return abs2(Ψtgt ⋅ Ψ_T)
    end
    grad_zygote = collect(gradient(F_sm, 1.0, 1.0))

    @test norm(Ψ̃_out.state - Û_Ψ) < 1e-12  # still correct?
    τ = Ψtgt ⋅ Û_Ψ
    # `grad` is gradient of F_sm based on Newton-prop of GradGenerator
    # For ∂F/∂τ see Eq. (3.47) of Phd Thesis of Michael Goerz
    grad = [
        2 * real(conj(τ) * (Ψtgt ⋅ Ψ̃_out.grad_states[1])),
        2 * real(conj(τ) * (Ψtgt ⋅ Ψ̃_out.grad_states[2]))
    ]
    @test abs(grad_zygote[1] - grad[1]) < 1e-9
    @test abs(grad_zygote[2] - grad[2]) < 1e-9

    ###########################################################################
    # Backward-propagation of GradVector
    # This is a change of perspective from τ = ⟨Ψtgt | ⋅ (Û |Ψ⟩)
    # to τ = (⟨Ψtgt | Û) ⋅ |Ψ⟩) = ⟨χ|Ψ⟩ with |χ⟩ = exp[i Ĥ⁺ dt] |Ψtgt⟩. That
    # is, instead of forward-propagating the initial state, we
    # backward-propagate the target state

    χ̃ = GradVector(Ψtgt, num_controls)
    wrk_bw = NewtonWrk(χ̃)
    Ĥ_adj_of_t = (copy(Ĥ₀'), (copy(Ĥ₁'), ϵ₁), (copy(Ĥ₂'), ϵ₂))
    G̃_adj_of_t = TimeDependentGradGenerator(Ĥ_adj_of_t)
    G̃_adj = evalcontrols(G̃_adj_of_t, vals_dict)
    newton!(χ̃, G̃_adj, -dt, wrk)
    χ̃_out = copy(χ̃)
    Ψtgt_Û = χ̃_out.state
    τ_bw = Ψtgt_Û ⋅ Ψ
    @test τ_bw ≈ τ
    grad_bw = [
        2 * real(conj(τ_bw) * (χ̃_out.grad_states[1] ⋅ Ψ)),
        2 * real(conj(τ_bw) * (χ̃_out.grad_states[2] ⋅ Ψ))
    ]
    @test grad_bw[1] ≈ grad[1]
    @test grad_bw[2] ≈ grad[2]

    ###########################################################################
    # Compare against Taylor series
    # This checks the gradients more directly against an alternative method

    commutator(A, B) = A * B - B * A

    """Evaluate ∂/∂ϵ exp(-𝕚 Ĥ dt) via a Taylor expansion."""
    function U_grad(Ĥ, μ̂, dt)
        # See Eq. (14) in de Fouquieres et. al, JMR 212, 412 (2011)
        Û = exp(-𝕚 * Ĥ * dt)
        converged = false
        Ĉ = μ̂
        terms = [(-𝕚 * dt) * Ĉ]
        n = 2
        while !converged
            Ĉ = commutator(Ĥ, Ĉ)
            term = -((𝕚 * dt)^n / factorial(big(n))) * Ĉ
            push!(terms, term)
            converged = (norm(term) < 1e-12)
            n += 1
        end
        return Û * sum(terms)
    end

    grad_taylor = [U_grad(Ĥ, Ĥ₁, dt) * Ψ, U_grad(Ĥ, Ĥ₂, dt) * Ψ]

    @test norm(Ψ̃_out.grad_states[1] - grad_taylor[1]) < 1e-10
    @test norm(Ψ̃_out.grad_states[2] - grad_taylor[2]) < 1e-10

end

@testset "Gradgen specrad" begin

    N = 10  # size of Hilbert space
    ρ = 1.0  # spectral radius
    Ĥ₀ = random_hermitian_matrix(N, ρ)
    Ĥ₁ = random_hermitian_matrix(N, ρ)
    Ĥ₂ = random_hermitian_matrix(N, ρ)
    Zero = zeros(ComplexF64, N, N)
    ϵ₁ = t -> 1.0
    ϵ₂ = t -> 1.0
    Ĥ_of_t = (Ĥ₀, (Ĥ₁, ϵ₁), (Ĥ₂, ϵ₂))
    vals_dict = IdDict(ϵ₁ => 1.0, ϵ₂ => 1.0)
    G̃_of_t = TimeDependentGradGenerator(Ĥ_of_t)
    Ĥ = evalcontrols(Ĥ_of_t, vals_dict)
    G̃ = evalcontrols(G̃_of_t, vals_dict)

    G_expected = [
         Ĥ    Zero  Ĥ₁
        Zero   Ĥ    Ĥ₂
        Zero  Zero  Ĥ
    ]

    G = Array(G̃)

    @test norm(G_expected - G) < 1e-14

    @test maximum(imag.(eigvals(G))) < 1e-15
    @test maximum(imag.(eigvals(Ĥ))) < 1e-15

    H_E_min = minimum(real.(eigvals(Ĥ)))
    H_E_max = maximum(real.(eigvals(Ĥ)))

    G_E_min = minimum(real.(eigvals(G)))
    G_E_max = maximum(real.(eigvals(G)))
    @test abs(H_E_min - G_E_min) < 1e-12
    @test abs(H_E_max - G_E_max) < 1e-12

    G_range_diag = collect(specrange(G̃, method=:diag))
    @test eltype(G_range_diag) ≡ Float64
    @test norm(G_range_diag - [G_E_min, G_E_max]) < 1e-12

    H_range_diag = collect(specrange(Ĥ, method=:diag))
    @test norm(H_range_diag - G_range_diag) < 1e-12

    G_range_arnoldi = collect(specrange(G̃, method=:arnoldi, m_max=100))
    @test eltype(G_range_arnoldi) ≡ Float64
    @test norm(G_range_arnoldi - [H_E_min, H_E_max]) < 1e-2
    # `specrange(Ĥ, method=:arnoldi)` isn't very exact, so we don't
    # compare against that.

end

@testset "Gradgen-sparse specrad" begin
    N = 100  # size of Hilbert space
    ρ = 1.0  # spectral radius
    sparsity = 0.1
    Ĥ₀ = random_hermitian_sparse_matrix(N, ρ, sparsity)
    Ĥ₁ = random_hermitian_sparse_matrix(N, ρ, sparsity)
    Ĥ₂ = random_hermitian_sparse_matrix(N, ρ, sparsity)
    Zero = zeros(ComplexF64, N, N)
    ϵ₁ = t -> 1.0
    ϵ₂ = t -> 1.0
    Ĥ_of_t = (Ĥ₀, (Ĥ₁, ϵ₁), (Ĥ₂, ϵ₂))
    vals_dict = IdDict(ϵ₁ => 1.0, ϵ₂ => 1.0)
    G̃_of_t = TimeDependentGradGenerator(Ĥ_of_t)
    Ĥ = evalcontrols(Ĥ_of_t, vals_dict)
    G̃ = evalcontrols(G̃_of_t, vals_dict)

    G_expected = Array([
         Ĥ    Zero  Ĥ₁
        Zero   Ĥ    Ĥ₂
        Zero  Zero  Ĥ
    ])

    G = Array(G̃)
    @test G isa Matrix{ComplexF64}

    @test norm(G_expected - G) < 1e-12

    @test maximum(imag.(eigvals(G))) < 1e-15
    @test maximum(imag.(eigvals(Array(Ĥ)))) < 1e-15

    H_E_min = minimum(real.(eigvals(Array(Ĥ))))
    H_E_max = maximum(real.(eigvals(Array(Ĥ))))

    G_E_min = minimum(real.(eigvals(G)))
    G_E_max = maximum(real.(eigvals(G)))
    @test abs(H_E_min - G_E_min) < 1e-12
    @test abs(H_E_max - G_E_max) < 1e-12

    G_range_diag = collect(specrange(G̃, method=:diag))
    @test eltype(G_range_diag) ≡ Float64
    @test norm(G_range_diag - [G_E_min, G_E_max]) < 1e-12

    G_range_arnoldi = collect(specrange(G̃, method=:arnoldi, m_max=100))
    @test eltype(G_range_arnoldi) ≡ Float64
    @test norm(G_range_arnoldi - [G_E_min, G_E_max]) < 1e-1

end