text:
function proposal2(current, cholV, lb, ub)
    trial = copy(current)
    if rand() > 0.1
        trial += cholV'*randn(size(trial))
    else
        trial = lb + (ub - lb).*rand(size(lb,1))
    end
    return trial
end

max_stars_repo_name: mcreel/SV

max_issues_repo_name: mcreel/SV

max_forks_repo_name: mcreel/SV


text:
using GeometricIntegrators
using GeometricIntegrators.TestProblems.HarmonicOscillatorProblem
using Test

set_config(:nls_atol_break, Inf)
set_config(:nls_rtol_break, Inf)
set_config(:nls_stol_break, Inf)

Δt = 0.1
nt = 10000
ns = 10

nsave1 = 1
nwrte1 = 100

nsave2 = 10
nwrte2 = 100

h5file = "test.hdf5"

tab = getTableauImplicitMidpoint()
ode = harmonic_oscillator_ode(vcat(rand(1,ns), zeros(1,ns)))

# @testset "$(rpad("Serial Simulation",80))" begin

    sim1 = Simulation(ode, tab, Δt, "Harmonic Oscillator Test 1", h5file, nt; nsave=nsave1, nwrite=nwrte1)
    run!(sim1)
    @test isfile(h5file)
    ssol1 = SSolutionODE(h5file)
    @test ssol1.ntime == nt
    @test ssol1.nsave == nsave1
    rm(h5file)

    sim2 = Simulation(ode, tab, Δt, "Harmonic Oscillator Test 2", h5file, nt; nsave=nsave2, nwrite=nwrte2)
    run!(sim2)
    @test isfile(h5file)
    ssol2 = SSolutionODE(h5file)
    @test ssol2.ntime == nt
    @test ssol2.nsave == nsave2
    rm(h5file)

# end

# @testset "$(rpad("Parallel Simulation (" * string(Threads.nthreads()) * " Threads)",80))" begin

    sim1 = ParallelSimulation(ode, tab, Δt, "Harmonic Oscillator Test 1 (" * string(Threads.nthreads()) * " Threads)", h5file, nt; nsave=nsave1, nwrite=nwrte1)
    run!(sim1)
    @test isfile(h5file)
    psol1 = SSolutionODE(h5file)
    @test psol1.ntime == nt
    @test psol1.nsave == nsave1
    rm(h5file)

    sim2 = ParallelSimulation(ode, tab, Δt, "Harmonic Oscillator Test 2 (" * string(Threads.nthreads()) * " Threads)", h5file, nt; nsave=nsave2, nwrite=nwrte2)
    run!(sim2)
    @test isfile(h5file)
    psol2 = SSolutionODE(h5file)
    @test psol2.ntime == nt
    @test psol2.nsave == nsave2
    rm(h5file)

# end

@test psol1.q == ssol1.q
@test psol2.q == ssol2.q

ext: jl


text:
# ============================================================
# Transformation Methods between nodal and modal formulations
# ============================================================

include("transform_triangle.jl")

"""
    vandermonde_matrix(N, r)
    vandermonde_matrix(N, r, s)

Compute Vandermonde matrix for node-mode transform

Vû = u

- @arg N: polynomial degree
- @arg r: local x axis
- @arg s: local y axis
"""
function vandermonde_matrix(N, r)
    V1D = zeros(eltype(r), length(r), N + 1)

    for j = 1:N+1
        V1D[:, j] .= JacobiP(r, 0, 0, j - 1)
    end

    return V1D
end

vandermonde_matrix(::Type{Line}, N, r) = vandermonde_matrix(N, r)

function vandermonde_matrix(::Type{Tri}, N, r, s)
    Np = (N + 1) * (N + 2) ÷ 2
    V2D = zeros(eltype(r), length(r), Np)
    a, b = rs_ab(r, s)

    sk = 1
    for i = 0:N
        for j = 0:N-i
            V2D[:, sk] .= simplex_basis(a, b, i, j)
            sk += 1
        end
    end

    return V2D
end

function vandermonde_matrix(::Type{Quad}, N, r, s)
    Np = (N + 1)^2
    V = zeros(eltype(r), length(r), Np)

    sk = 1
    for i = 0:N
        for j = 0:N
            V[:, sk] = JacobiP(r, 0, 0, i) .* JacobiP(s, 0, 0, j)
            sk += 1
        end
    end

    return V
end

"""
    ∂vandermonde_matrix(N, r)
    ∂vandermonde_matrix(N, r, s)

gradient of the modal basis (i,j) at (r,s) at order N

"""
function ∂vandermonde_matrix(N::T, r) where {T<:Integer}
    Vr = zeros(length(r), N + 1)

    for i = 0:N
        Vr[:, i+1] .= ∂JacobiP(r, 0, 0, i)
    end

    return Vr
end

∂vandermonde_matrix(::Type{Line}, N, r) = ∂vandermonde_matrix(N, r)

function ∂vandermonde_matrix(::Type{Tri}, N::T, r, s) where {T<:Integer}
    V2Dr = zeros(eltype(r), length(r), (N + 1) * (N + 2) ÷ 2)
    V2Ds = zeros(eltype(r), length(r), (N + 1) * (N + 2) ÷ 2)

    # tensor-product coordinates
    a, b = rs_ab(r, s)

    sk = 1
    for i = 0:N
        for j = 0:N-i
            V2Dr[:, sk], V2Ds[:, sk] = ∂simplex_basis(a, b, i, j)
            sk += 1
        end
    end

    return V2Dr, V2Ds
end

function ∂vandermonde_matrix(::Type{Quad}, N::T, r, s) where {T<:Integer}
    V2Dr = zeros(eltype(r), length(r), (N + 1)^2)
    V2Ds = zeros(eltype(r), length(r), (N + 1)^2)

    sk = 1
    for i = 0:N
        for j = 0:N
            V2Dr[:, sk] .= ∂JacobiP(r, 0, 0, i) .* JacobiP(s, 0, 0, j)
            V2Ds[:, sk] .= JacobiP(r, 0, 0, i) .* ∂JacobiP(s, 0, 0, j)
            sk += 1
        end
    end

    return V2Dr, V2Ds
end


text:
"""
    VectorSpaceType

Abstract type for tangent spaces, cotangent spaces, their tensor products,
exterior products, etc.

Every vector space `VS` is supposed to provide:
* a method of constructing vectors,
* basic operations: addition, subtraction, multiplication by a scalar
  and negation (unary minus),
* [`zero_vector!(VS, v, x)`](@ref) to construct zero vectors at point `x`,
* `allocate(v)` and `allocate(v, T)` for vector `v` and type `T`,
* `copyto!(v, w)` for vectors `v` and `w`,
* `number_eltype(v)` for vector `v`,
* [`vector_space_dimension(::VectorBundleFibers{<:typeof(VS)}) where VS`](@ref).

Optionally:
* inner product via `inner` (used to provide Riemannian metric on vector
  bundles),
* [`flat`](@ref) and [`sharp`](@ref),
* `norm` (by default uses `inner`),
* [`project_vector`](@ref) (for embedded vector spaces),
* [`representation_size`](@ref) (if support for [`ProductArray`](@ref) is desired),
* broadcasting for basic operations.
"""
abstract type VectorSpaceType end

struct TangentSpaceType <: VectorSpaceType end

struct CotangentSpaceType <: VectorSpaceType end

TCoTSpaceType = Union{TangentSpaceType,CotangentSpaceType}

const TangentSpace = TangentSpaceType()
const CotangentSpace = CotangentSpaceType()

"""
    TensorProductType(spaces::VectorSpaceType...)

Vector space type corresponding to the tensor product of given vector space
types.
"""
struct TensorProductType{TS<:Tuple} <: VectorSpaceType
    spaces::TS
end

TensorProductType(spaces::VectorSpaceType...) = TensorProductType{typeof(spaces)}(spaces)

"""
    VectorBundleFibers(VS::VectorSpaceType, M::Manifold)

Type representing a family of vector spaces (fibers) of a vector bundle over `M`
with vector spaces of type `VS`. In contrast with `VectorBundle`, operations
on `VectorBundleFibers` expect point-like and vector-like parts to be
passed separately instead of being bundled together. It can be thought of
as a representation of vector spaces from a vector bundle but without
storing the point at which a vector space is attached (which is specified
separately in various functions).
"""
struct VectorBundleFibers{TVS<:VectorSpaceType,TM<:Manifold}
    VS::TVS
    M::TM
end

const TangentBundleFibers{M} = VectorBundleFibers{TangentSpaceType,M}

TangentBundleFibers(M::Manifold) = VectorBundleFibers(TangentSpace, M)

const CotangentBundleFibers{M} = VectorBundleFibers{CotangentSpaceType,M}

CotangentBundleFibers(M::Manifold) = VectorBundleFibers(CotangentSpace, M)

"""
    VectorSpaceAtPoint(fiber::VectorBundleFibers, x)

A vector space (fiber type `fiber` of a vector bundle) at point `x` from
the manifold `fiber.M`.
"""
struct VectorSpaceAtPoint{TFiber<:VectorBundleFibers,TX}
    fiber::TFiber
    x::TX
end

"""
    TangentSpaceAtPoint(M::Manifold, x)

Return an object of type [`VectorSpaceAtPoint`](@ref) representing tangent
space at `x`.
"""
TangentSpaceAtPoint(M::Manifold, x) = VectorSpaceAtPoint(TangentBundleFibers(M), x)

"""
    CotangentSpaceAtPoint(M::Manifold, x)

Return an object of type [`VectorSpaceAtPoint`](@ref) representing cotangent
space at `x`.
"""
CotangentSpaceAtPoint(M::Manifold, x) = VectorSpaceAtPoint(CotangentBundleFibers(M), x)

"""
    VectorBundle(M::Manifold, type::VectorSpaceType)

Vector bundle on manifold `M` of type `type`.
"""
struct VectorBundle{TVS<:VectorSpaceType,TM<:Manifold} <: Manifold
    type::TVS
    M::TM
    VS::VectorBundleFibers{TVS,TM}
end

function VectorBundle(VS::TVS, M::TM) where {TVS<:VectorSpaceType,TM<:Manifold}
    return VectorBundle{TVS,TM}(VS, M, VectorBundleFibers(VS, M))
end

const TangentBundle{M} = VectorBundle{TangentSpaceType,M}

TangentBundle(M::Manifold) = VectorBundle(TangentSpace, M)

const CotangentBundle{M} = VectorBundle{CotangentSpaceType,M}

CotangentBundle(M::Manifold) = VectorBundle(CotangentSpace, M)

"""
    FVector(type::VectorSpaceType, data)

Decorator indicating that the vector `data` is from a fiber of a vector bundle
of type `type`.
"""
struct FVector{TType<:VectorSpaceType,TData}
    type::TType
    data::TData
end

const TFVector = FVector{TangentSpaceType}
const CoTFVector = FVector{CotangentSpaceType}

struct PrecomputedVectorBundleOrthonormalBasis{
    F,
    TBase<:AbstractPrecomputedOrthonormalBasis{F},
    TVec<:AbstractPrecomputedOrthonormalBasis{F},
} <: AbstractPrecomputedOrthonormalBasis{F}
    base_basis::TBase
    vec_basis::TVec
end

(+)(v1::FVector, v2::FVector) = FVector(v1.type, v1.data + v2.data)

(-)(v1::FVector, v2::FVector) = FVector(v1.type, v1.data - v2.data)
(-)(v::FVector) = FVector(v.type, -v.data)

(*)(a::Number, v::FVector) = FVector(v.type, a * v.data)

function copyto!(y::FVector, x::FVector)
    copyto!(y.data, x.data)
    return y
end

base_manifold(B::VectorBundleFibers) = base_manifold(B.M)
base_manifold(B::VectorSpaceAtPoint) = base_manifold(B.fiber)
base_manifold(B::VectorBundle) = base_manifold(B.M)

"""
    bundle_projection(B::VectorBundle, x::ProductRepr)

Projection of point `x` from the bundle `M` to the base manifold.
Returns the point on the base manifold `B.M` at which the vector part
of `x` is attached.
"""
bundle_projection(B::VectorBundle, x) = submanifold_component(B.M, x, Val(1))

"""
    distance(B::VectorBundleFibers, x, v, w)

Distance between vectors `v` and `w` from the vector space at point `x`
from the manifold `M.M`, that is the base manifold of `M`.
"""
distance(B::VectorBundleFibers, x, v, w) = norm(B, x, v - w)
@doc doc"""
    distance(B::VectorBundle, x, y)

Distance between points $x$ and $y$ from the
vector bundle `B` over manifold `B.VS` (denoted $M$).

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
    Similarly, $y = (p_y, \xi_y)$.

The distance is calculated as

$d_B(x, y) = \sqrt{d_M(p_x, p_y)^2 + d_F(\xi_x, \xi_{y\to x})^2}$

where $d_M$ is the distance on manifold $M$, $d_F$ is the distance
between two vectors from the fiber $F$ and $\xi_{y\to x}$ is the result
of parallel transport of vector $\xi_y$ to point $p_x$. The default
behavior of [`vector_transport_to`](@ref) is used to compute the vector
transport.
"""
function distance(B::VectorBundle, x, y)
    px, ξx = submanifold_components(B.M, x)
    py, ξy = submanifold_components(B.M, y)
    dist_man = distance(B.M, px, py)
    vy_x = vector_transport_to(B.M, py, ξy, px)
    dist_vec = distance(B.VS, px, ξx, vy_x)
    return sqrt(dist_man^2 + dist_vec^2)
end

number_eltype(::Type{FVector{TType,TData}}) where {TType<:VectorSpaceType,TData} = number_eltype(TData)
number_eltype(v::FVector) = number_eltype(v.data)

@doc doc"""
    exp(B::VectorBundle, x, v)

Exponential map of tangent vector $v$ at point $x$ from
vector bundle `B` over manifold `B.VS` (denoted $M$).

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
  * The tangent vector $v = (\xi_{v,M}, \xi_{v,F}) \in T_{x}B$ where
    $\xi_{v,M}$ is a tangent vector from the tangent space $T_{p_x}M$ and
    $\xi_{v,F}$ is a tangent vector from the tangent space $T_{\xi_x}F$ (isomorphic to $F$).

The exponential map is calculated as

$\exp_{x}(v) = (\exp_{p_x}(\xi_{v,M}), \xi_{\exp})$

where $\xi_{\exp}$ is the result of vector transport of $\xi_x + \xi_{v,F}$
to the point $\exp_{p_x}(\xi_{v,M})$.
The sum $\xi_x + \xi_{v,F}$ corresponds to the exponential map in the vector space $F$.
"""
exp(::VectorBundle, ::Any)

function exp!(B::VectorBundle, y, x, v)
    px, ξx = submanifold_components(B.M, x)
    py, ξy = submanifold_components(B.M, y)
    ξvM, ξvF = submanifold_components(B.M, v)
    exp!(B.M, py, px, ξvM)
    vector_transport_to!(B.M, ξy, px, ξx + ξvF, py)
    return y
end

@doc doc"""
    flat(M::Manifold, x, w::FVector)

Compute the flat isomorphism (one of the musical isomorphisms) of tangent vector `w`
from the vector space of type `M` at point `x` from the underlying [`Manifold`](@ref).

The function can be used for example to transform vectors
from the tangent bundle to vectors from the cotangent bundle
$\flat \colon T\mathcal M \to T^{*}\mathcal M$
"""
function flat(M::Manifold, x, w::FVector)
    v = allocate_result(M, flat, w, x)
    return flat!(M, v, x, w)
end

function flat!(M::Manifold, v::FVector, x, w::FVector)
    error(
        "flat! not implemented for vector bundle fibers space " *
        "of type $(typeof(M)), vector of type $(typeof(v)), point of " *
        "type $(typeof(x)) and vector of type $(typeof(w)).",
    )
end

function get_basis(M::VectorBundle, x, B::DiagonalizingOrthonormalBasis)
    xp1 = submanifold_component(x, Val(1))
    bv1 = DiagonalizingOrthonormalBasis(submanifold_component(B.v, Val(1)))
    b1 = get_basis(M.M, xp1, bv1)
    bv2 = DiagonalizingOrthonormalBasis(submanifold_component(B.v, Val(2)))
    b2 = get_basis(M.VS, xp1, bv2)
    return PrecomputedVectorBundleOrthonormalBasis(b1, b2)
end
function get_basis(M::TangentBundleFibers, x, B::DiagonalizingOrthonormalBasis)
    return get_basis(M.M, x, B)
end

function get_coordinates(M::VectorBundle, x, v, B::ArbitraryOrthonormalBasis) where {N}
    px, ξx = submanifold_components(M.M, x)
    ξvM, ξvF = submanifold_components(M.M, v)
    coord1 = get_coordinates(M.M, px, ξvM, B)
    coord2 = get_coordinates(M.VS, px, ξvF, B)
    return vcat(coord1, coord2)
end
function get_coordinates(
    M::VectorBundle,
    x,
    v,
    B::PrecomputedVectorBundleOrthonormalBasis,
) where {N}
    px, ξx = submanifold_components(M.M, x)
    ξvM, ξvF = submanifold_components(M.M, v)
    coord1 = get_coordinates(M.M, px, ξvM, B.base_basis)
    coord2 = get_coordinates(M.VS, px, ξvF, B.vec_basis)
    return vcat(coord1, coord2)
end
function get_coordinates(M::TangentBundleFibers, x, v, B::AbstractBasis) where {N}
    return get_coordinates(M.M, x, v, B)
end

function get_vector(M::VectorBundle, x, v, B::ArbitraryOrthonormalBasis) where {N}
    mdim = manifold_dimension(M.M)
    xp1 = submanifold_component(x, Val(1))
    v1 = get_vector(M.M, xp1, v[1:mdim], B)
    v2 = get_vector(M.VS, xp1, v[mdim+1:end], B)
    return ProductRepr(v1, v2)
end
function get_vector(
    M::VectorBundle,
    x,
    v,
    B::PrecomputedVectorBundleOrthonormalBasis,
) where {N}
    mdim = manifold_dimension(M.M)
    xp1 = submanifold_component(x, Val(1))
    v1 = get_vector(M.M, xp1, v[1:mdim], B.base_basis)
    v2 = get_vector(M.VS, xp1, v[mdim+1:end], B.vec_basis)
    return ProductRepr(v1, v2)
end
function get_vector(M::TangentBundleFibers, x, v, B::AbstractBasis) where {N}
    return get_vector(M.M, x, v, B)
end
function get_vectors(M::VectorBundle, x, B::PrecomputedVectorBundleOrthonormalBasis)
    xp1 = submanifold_component(x, Val(1))
    zero_m = zero_tangent_vector(M.M, xp1)
    zero_f = zero_vector(M.VS, xp1)
    vs = typeof(ProductRepr(zero_m, zero_f))[]
    for bv in get_vectors(M.M, xp1, B.base_basis)
        push!(vs, ProductRepr(bv, zero_f))
    end
    for bv in get_vectors(M.VS, xp1, B.vec_basis)
        push!(vs, ProductRepr(zero_m, bv))
    end
    return vs
end

get_vectors(::VectorBundleFibers, x, B::PrecomputedOrthonormalBasis) = B.vectors
function get_vectors(::VectorBundleFibers, x, B::PrecomputedDiagonalizingOrthonormalBasis)
    return B.vectors
end

Base.@propagate_inbounds getindex(x::FVector, i) = getindex(x.data, i)

"""
    inner(B::VectorBundleFibers, x, v, w)

Inner product of vectors `v` and `w` from the vector space of type `B.VS`
at point `x` from manifold `B.M`.
"""
function inner(B::VectorBundleFibers, x, v, w)
    error(
        "inner not defined for vector space family of type $(typeof(B)), " *
        "point of type $(typeof(x)) and " *
        "vectors of types $(typeof(v)) and $(typeof(w)).",
    )
end
inner(B::VectorBundleFibers{<:TangentSpaceType}, x, v, w) = inner(B.M, x, v, w)
function inner(B::VectorBundleFibers{<:CotangentSpaceType}, x, v, w)
    return inner(
        B.M,
        x,
        sharp(B.M, x, FVector(CotangentSpace, v)).data,
        sharp(B.M, x, FVector(CotangentSpace, w)).data,
    )
end
@doc doc"""
    inner(B::VectorBundle, x, v, w)

Inner product of tangent vectors `v` and `w` at point `x` from the
vector bundle `B` over manifold `B.VS` (denoted $M$).

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
  * The tangent vector $v = (\xi_{v,M}, \xi_{v,F}) \in T_{x}B$ where
    $\xi_{v,M}$ is a tangent vector from the tangent space $T_{p_x}M$ and
    $\xi_{v,F}$ is a tangent vector from the tangent space $T_{\xi_x}F$ (isomorphic to $F$).
    Similarly for the other tangent vector $w = (\xi_{w,M}, \xi_{w,F}) \in T_{x}B$.

The inner product is calculated as

$\langle v, w \rangle_{B} = \langle \xi_{v,M}, \xi_{w,M} \rangle_{M} + \langle \xi_{v,F}, \xi_{w,F} \rangle_{F}.$
"""
function inner(B::VectorBundle, x, v, w)
    px, ξx = submanifold_components(B.M, x)
    ξvM, ξvF = submanifold_components(B.M, v)
    ξwM, ξwF = submanifold_components(B.M, w)
    return inner(B.M, px, ξvM, ξwM) + inner(B.VS, ξx, ξvF, ξwF)
end

function isapprox(B::VectorBundle, x, y; kwargs...)
    px, ξx = submanifold_components(B.M, x)
    py, ξy = submanifold_components(B.M, y)
    return isapprox(B.M, px, py; kwargs...) && isapprox(ξx, ξy; kwargs...)
end
function isapprox(B::VectorBundle, x, v, w; kwargs...)
    px, ξx = submanifold_components(B.M, x)
    ξvM, ξvF = submanifold_components(B.M, v)
    ξwM, ξwF = submanifold_components(B.M, w)
    return isapprox(B.M, ξvM, ξwM; kwargs...) && isapprox(B.M, px, ξvF, ξwF; kwargs...)
end

@doc doc"""
    log(B::VectorBundle, x, y)

Logarithmic map of the point $y$ at point $x$ from
vector bundle `B` over manifold `B.VS` (denoted $M$).

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
    Similarly, $y = (p_y, \xi_y)$.

The logarithmic map is calculated as

$\log_{x}(y) = (\log_{p_x}(p_y), \xi_{\log} - \xi_x)$

where $\xi_{\log}$ is the result of vector transport of $\xi_y$
to the point $p_x$.
The difference $\xi_{\log} - \xi_x$ corresponds to the logarithmic map in the vector space $F$.
"""
log(::VectorBundle, ::Any...)

function log!(B::VectorBundle, v, x, y)
    px, ξx = submanifold_components(B.M, x)
    py, ξy = submanifold_components(B.M, y)
    ξvM, ξvF = submanifold_components(B.M, v)
    log!(B.M, ξvM, px, py)
    vector_transport_to!(B.M, ξvF, py, ξy, px)
    copyto!(ξvF, ξvF - ξx)
    return v
end

manifold_dimension(B::VectorBundle) = manifold_dimension(B.M) + vector_space_dimension(B.VS)

"""
    norm(B::VectorBundleFibers, x, v)

Norm of the vector `v` from the vector space of type `B.VS`
at point `x` from manifold `B.M`.
"""
norm(B::VectorBundleFibers, x, v) = sqrt(inner(B, x, v, v))
norm(B::VectorBundleFibers{<:TangentSpaceType}, x, v) = norm(B.M, x, v)

@doc doc"""
    project_point(B::VectorBundle, x)

Project the point $x$ from the ambient space of the vector bundle `B`
over manifold `B.VS` (denoted $M$) to the vector bundle.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x$ belongs to the ambient space of $M$
    and $\xi_x$ belongs to the ambient space of the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.

The projection is calculated by projecting the point $p_x$ to the manifold $M$
and then projecting the vector $\xi_x$ to the tangent space $T_{p_x}M$.
"""
project_point(::VectorBundle, ::Any...)

function project_point!(B::VectorBundle, x)
    px, ξx = submanifold_components(B.M, x)
    project_point!(B.M, px)
    project_tangent!(B.M, ξx, px, ξx)
    return x
end

@doc doc"""
    project_tangent(B::VectorBundle, x, v)

Project the element $v$ of the ambient space of the tangent space $T_x B$
to the tangent space $T_x B$.

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.
  * The vector $x = (\xi_{v,M}, \xi_{v,F})$ where $p_x$ belongs to the ambient space of $T_{p_x}M$
    and $\xi_{v,F}$ belongs to the ambient space of the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.

The projection is calculated by projecting $\xi_{v,M}$ to tangent space $T_{p_x}M$
and then projecting the vector $\xi_{v,F}$ to the fiber $F$.
"""
project_tangent(::VectorBundle, ::Any...)

function project_tangent!(B::VectorBundle, w, x, v)
    px, ξx = submanifold_components(B.M, x)
    ξvM, ξvF = submanifold_components(B.M, v)
    ξwM, ξwF = submanifold_components(B.M, w)
    project_tangent!(B.M, ξwM, px, ξvM)
    project_tangent!(B.M, ξwF, px, ξvF)
    return w
end

"""
    project_vector(B::VectorBundleFibers, x, w)

Project vector `w` from the vector space of type `B.VS` at point `x`.
"""
function project_vector(B::VectorBundleFibers, x, w)
    v = allocate_result(B, project_vector, x, w)
    return project_vector!(B, v, x, w)
end

function project_vector!(B::VectorBundleFibers{<:TangentSpaceType}, v, x, w)
    return project_tangent!(B.M, v, x, w)
end
function project_vector!(B::VectorBundleFibers, v, x, w)
    error("project_vector! not implemented for vector space family of type $(typeof(B)), output vector of type $(typeof(v)) and input vector at point $(typeof(x)) with type of w $(typeof(w)).")
end

Base.@propagate_inbounds setindex!(x::FVector, val, i) = setindex!(x.data, val, i)

representation_size(B::VectorBundleFibers{<:TCoTSpaceType}) = representation_size(B.M)
function representation_size(B::VectorBundle)
    len_manifold = prod(representation_size(B.M))
    len_vs = prod(representation_size(B.VS))
    return (len_manifold + len_vs,)
end

@doc doc"""
    sharp(M::Manifold, x, w::FVector)

Compute the sharp isomorphism (one of the musical isomorphisms) of vector `w`
from the vector space `M` at point `x` from the underlying [`Manifold`](@ref).

The function can be used for example to transform vectors
from the cotangent bundle to vectors from the tangent bundle
$\sharp \colon T^{*}\mathcal M \to T\mathcal M$
"""
function sharp(M::Manifold, x, w::FVector)
    v = allocate_result(M, sharp, w, x)
    return sharp!(M, v, x, w)
end

function sharp!(M::Manifold, v::FVector, x, w::FVector)
    error(
        "sharp! not implemented for vector bundle fibers space " *
        "of type $(typeof(M)), vector of type $(typeof(v)), point of " *
        "type $(typeof(x)) and vector of type $(typeof(w)).",
    )
end

allocate(x::FVector) = FVector(x.type, allocate(x.data))
allocate(x::FVector, ::Type{T}) where {T} = FVector(x.type, allocate(x.data, T))

"""
    allocate_result(B::VectorBundleFibers, f, x...)

Allocates an array for the result of function `f` that is
an element of the vector space of type `B.VS` on manifold `B.M`
and arguments `x...` for implementing the non-modifying operation
using the modifying operation.
"""
function allocate_result(B::VectorBundleFibers, f, x...)
    T = allocate_result_type(B, f, x)
    return allocate(x[1], T)
end
function allocate_result(M::Manifold, ::typeof(flat), w::TFVector, x)
    return FVector(CotangentSpace, allocate(w.data))
end
function allocate_result(M::Manifold, ::typeof(sharp), w::CoTFVector, x)
    return FVector(TangentSpace, allocate(w.data))
end

"""
    allocate_result_type(B::VectorBundleFibers, f, args::NTuple{N,Any}) where N

Returns type of element of the array that will represent the result of
function `f` for representing an operation with result in the vector space `VS`
for manifold `M` on given arguments (passed at a tuple).
"""
function allocate_result_type(B::VectorBundleFibers, f, args::NTuple{N,Any}) where N
    T = typeof(reduce(+, one(number_eltype(eti)) for eti ∈ args))
    return T
end

size(x::FVector) = size(x.data)

function submanifold_component(M::Manifold, x::FVector, i::Val)
    return submanifold_component(M, x.data, i)
end
submanifold_component(x::FVector, i::Val) = submanifold_component(x.data, i)

submanifold_components(M::Manifold, x::FVector) = submanifold_components(M, x.data)
submanifold_components(x::FVector) = submanifold_components(x.data)

"""
    vector_space_dimension(B::VectorBundleFibers)

Dimension of the vector space of type `B`.
"""
function vector_space_dimension(B::VectorBundleFibers)
    error("vector_space_dimension not implemented for vector space family $(typeof(B)).")
end
vector_space_dimension(B::VectorBundleFibers{<:TCoTSpaceType}) = manifold_dimension(B.M)
function vector_space_dimension(B::VectorBundleFibers{<:TensorProductType})
    dim = 1
    for space in B.VS.spaces
        dim *= vector_space_dimension(VectorBundleFibers(space, B.M))
    end
    return dim
end

"""
    zero_vector(B::VectorBundleFibers, x)

Compute the zero vector from the vector space of type `B.VS` at point `x`
from manifold `B.M`.
"""
function zero_vector(B::VectorBundleFibers, x)
    v = allocate_result(B, zero_vector, x)
    return zero_vector!(B, v, x)
end

"""
    zero_vector!(B::VectorBundleFibers, v, x)

Save the zero vector from the vector space of type `B.VS` at point `x`
from manifold `B.M` to `v`.
"""
function zero_vector!(B::VectorBundleFibers, v, x)
    error("zero_vector! not implemented for vector space family of type $(typeof(B)).")
end
function zero_vector!(B::VectorBundleFibers{<:TangentSpaceType}, v, x)
    return zero_tangent_vector!(B.M, v, x)
end

@doc doc"""
    zero_tangent_vector(B::VectorBundle, x)

Zero tangent vector at point $x$ from the vector bundle `B`
over manifold `B.VS` (denoted $M$). The zero vector belongs to the space $T_{x}B$

Notation:
  * The point $x = (p_x, \xi_x)$ where $p_x \in M$ and $\xi_x$ belongs to the
    fiber $F=\pi^{-1}(\{p_x\})$ of the vector bundle $B$ where $\pi$ is the
    canonical projection of that vector bundle $B$.

The zero vector is calculated as

$\mathbf{0}_{x} = (\mathbf{0}_{p_x}, \mathbf{0}_F)$

where $\mathbf{0}_{p_x}$ is the zero tangent vector from $T_{p_x}M$ and
$\mathbf{0}_F$ is the zero element of the vector space $F$.
"""
zero_tangent_vector(::VectorBundle, ::Any...)

function zero_tangent_vector!(B::VectorBundle, v, x)
    px, ξx = submanifold_components(B.M, x)
    ξvM, ξvF = submanifold_components(B.M, v)
    zero_tangent_vector!(B.M, ξvM, px)
    zero_vector!(B.VS, ξvF, ξx)
    return v
end


text:
@with_kw struct AndersonCache{Tm, Tv, Tmcv, Tmgv, Tmrv, Tmqv, Tvγ, Tvmc}
    G::Tm
    Gv::Tmgv
    R::Tm
    Rv::Tmrv
    Rcv::Tmcv
    Q::Tm
    Qv::Tmqv
    fold::Tv
    fcur::Tv
    Δf::Tv
    gold::Tv
    gcur::Tv
    Δg::Tv
    m::Int
    γ::Tvγ
    m_eff_cache::Tvmc
    x_cache::Tv
end

function AndersonCache(x0, m)
    n = length(x0)
    m_corrected = min(n, m)
    G = similar(x0, n, m_corrected)
    Gv = [view(G, :, 1:i) for i = 1:m_corrected]
    R = zeros(eltype(x0), m_corrected, m_corrected)
    Rv = [view(R, 1:i, 1:i) for i = 1:m_corrected]
    Rcv = [similar(x0, i, i) for i = 1:m_corrected]
    Q = similar(x0, n, m_corrected)
    Qv = [view(Q, :, 1:i) for i = 1:m_corrected]
    fold = similar(x0)
    fcur = similar(x0)
    Δf = similar(fcur)
    gold = similar(x0)
    gcur = similar(x0)
    Δg = similar(gcur)
    γ = similar(x0, m)
    γv = [view(γ, 1:i) for i = 1:m_corrected]

    m_eff_cache = similar(x0, m_corrected)
    m_cache = [view(m_eff_cache, 1:i) for i = 1:m_corrected]
    x_cache = similar(x0)

    AndersonCache(G, Gv, R, Rv, Rcv, Q, Qv, fold, fcur, Δf, gold, gcur, Δg, m_corrected, γv, m_cache, x_cache)
end

# if some term converges exactly, you can get lapack errors (try to set an x to the true value)
function anderson(g, x0; itermax=1000, m=10, delay=0, mutate_x0=false, droptol=1e10, beta=nothing)
    if mutate_x0
        x = x0
    else
        x = copy(x0)
    end
    cache = AndersonCache(x0, m)
    anderson!(g, x, cache; itermax=itermax, delay=delay, droptol=droptol, beta=beta)
end

function anderson!(g, x::AbstractArray{T}, cache; itermax = 1000, delay=0, tol=sqrt(eps(T)), droptol=T(1e10), beta=nothing) where T
    @unpack G, Gv, R, Rv, Rcv, Q, Qv, fold, fcur, Δf, gold, gcur, Δg, m, γ, m_eff_cache, x_cache = cache
    function_iterations = delay+1 # always do one function iteration to set up gold and fold
    x, function_converged = function_iteration!(g, x, function_iterations, tol, gcur, fcur)
    if function_converged
        println("Converged before acceleration began")
        return x, true, 0
    end
    # k is now 0
    if m == 0
        x, function_converged = function_iteration!(g, x, itermax, tol, gcur, fcur)
        if function_converged
            println("Converged without acceleration (m=0).")
            return x, true, 0
        else
            return x, false, 0
        end
    end

    # Since we always do one function iteration we can store the current values
    # in gold and fold. We've already checked for convergence.
    gold .= gcur
    fold .= fcur
    # m_eff is the current history length (<= cache)
    m_eff = 0
    verbose = false
    # start acceleration
    for k = 1:itermax
        #gcur .= g(x)
        g(gcur, x)
        fcur .= gcur .- x
        verbose && println()
        verbose && println("$k) Residual: ", norm(fcur, Inf))
        if norm(fcur, Inf) < tol
            return x, true, k
        end
        Δf .= fcur .- fold
        Δg .= gcur .- gold

        # now that we have the difference, we can update the cache arrays
        # for f and g
        copyto!(fold, fcur)
        copyto!(gold, gcur)

        # if effective memory counter is below the maximum
        # memory length, use the iteration index. Otherwise, update from the
        # right.
        m_eff = m_eff + 1
        if m_eff <= m
            G[:, m_eff] = Δg
        else
            for i = 1:size(G, 1)
                for j = 2:m
                    G[i, j-1] = G[i, j]
                end
            end
            G[:, end] = Δg
        end
        # increment effective memory counter

        # if m_eff == 1 simply update R[1, 1] and Q[:, 1], else we need to
        # loop over rows in R[:, m_eff] and columns in Q. Potentially, we
        # need to delete a column from the right in Q and the last row/column
        # in R first, if the effective memory counter is larger than the
        # maximum memory length.
        if m_eff > 1
            if m_eff > m
                m_eff = m_eff - 1
                qrdelete!(Qv[m_eff], Rv[m_eff])
            end
            for i = 1:m_eff-1
                # R[i, m_eff] = dot(Q[:, i], Δf)
                R[i, m_eff] = T(0)
                for j = 1:length(Δf)
                    R[i, m_eff] += Q[j, i] * Δf[j]
                end
#                Δf .= Δf .- R[i, m_eff].*Q[:, i]
                for j = 1:length(Δf)
                    Δf[j] = Δf[j] - R[i, m_eff]*Q[j, i]
                end
            end
        end
       R[m_eff, m_eff] = norm(Δf, 2)
       Q[:, m_eff] .= Δf./R[m_eff, m_eff]

        if !isa(droptol, Nothing)
            condR = cond(Rv[m_eff])
            while condR > droptol && m_eff > 1
                qrdelete!(Qv[m_eff], Rv[m_eff])
                for j = 2:m_eff
                    for i = 1:size(G, 1)
                        G[i, j-1] = G[i, j]
                    end
                end
                m_eff = m_eff - 1
                condR = cond(Rv[m_eff])
            end
        end

        Rcv[m_eff] .= Rv[m_eff]
        #    γ[1:m_eff] .= Rv[m_eff]\Qv[m_eff]'*fcur
        F = qr!(Rcv[m_eff])
        ldiv!(γ[m_eff], F, mul!(m_eff_cache[m_eff], Qv[m_eff]', fcur))
        x .= gcur .- mul!(x_cache, Gv[m_eff], γ[m_eff])
        if !isa(beta, Nothing)
            x .= x .- (1 .- beta).*(fcur .- mul!(x_cache, Qv[m_eff], mul!(m_eff_cache[m_eff], Rv[m_eff], γ[m_eff])))
        end
    end
    printstyled("!!!!!!!!!\n"; color=9)
    printstyled("!Failure!\n"; color=9)
    printstyled("!!!!!!!!!\n"; color=9)
    printstyled("Exited with\n"; color=9)
    printstyled(" * residual inf-norm $(norm(fcur, Inf))\n"; color=9)
    printstyled(" * effective memory length of $m_eff\n"; color=9)
    printstyled("and\n"; color=9)
    printstyled(" - ran for $delay initial function iterations\n"; color=9)
    printstyled(" - ran for $itermax accelerated iterations\n"; color=9)

    x, false, itermax
end


text:
#transfer
transfer_left(v::MPSBondTensor, A::GenericMPSTensor, Ab::GenericMPSTensor=A) =
    _permute_as(_permute_front(Ab)' * _permute_front(_permute_front(v)*_permute_tail(A)), v)
transfer_right(v::MPSBondTensor, A::GenericMPSTensor, Ab::GenericMPSTensor=A) =
    _permute_as(_permute_tail(_permute_front(A)*_permute_front(v)) * _permute_tail(Ab)', v)

#transfer for 2 mpo tensors
transfer_left(v::MPSBondTensor,A::MPOTensor,B::MPOTensor) = @tensor t[-1;-2] := v[1,2]*A[2,3,-2,4]*conj(B[1,3,-1,4])
transfer_right(v::MPSBondTensor,A::MPOTensor,B::MPOTensor) = @tensor t[-1;-2] := A[-1,3,1,4]*conj(B[-2,3,2,4])*v[1,2]

#transfer, but there is a utility leg in the middle that is passed through
transfer_left(v::MPSTensor, A::MPSTensor, Ab::MPSTensor=A) = @tensor v[-1 -2;-3] := v[1,-2,2]*A[2,3,-3]*conj(Ab[1,3,-1])
transfer_right(v::MPSTensor, A::MPSTensor, Ab::MPSTensor=A) = @tensor v[-1 -2;-3] := A[-1,3,1]*v[1,-2,2]*conj(Ab[-3,3,2])

#mpo transfer
transfer_left(v::MPSTensor,O::MPOTensor,A::MPSTensor,Ab::MPSTensor) = @tensor v[-1 -2;-3] := v[4,2,1]*A[1,3,-3]*O[2,5,-2,3]*conj(Ab[4,5,-1])
transfer_right(v::MPSTensor,O::MPOTensor,A::MPSTensor,Ab::MPSTensor) = @tensor v[-1 -2;-3] := A[-1,4,5]*O[-2,2,3,4]*conj(Ab[-3,2,1])*v[5,3,1]

#utility, allowing transfering with arrays
function transfer_left(v,A::AbstractArray,Ab::AbstractArray=A;rvec=nothing,lvec=nothing)
    for (a,b) in zip(A,Ab)
        v = transfer_left(v,a,b)
    end

    if rvec != nothing && lvec != nothing
        if v isa MPSBondTensor #normal transfer
            @tensor v[-1;-2]-=rvec[1,2]*v[2,1]*lvec[-1,-2]
        elseif v isa MPSTensor #utiity leg in the middle
            @tensor v[-1 -2;-3]-=rvec[1,2]*v[2,-2,1]*lvec[-1,-3]
        else #what have you just given me?
            @assert false
        end
    end

    return v
end
function transfer_right(v,A::AbstractArray,Ab::AbstractArray=A;rvec=nothing,lvec=nothing)
    for (a,b) in Iterators.reverse(zip(A,Ab))
        v = transfer_right(v,a,b)
    end

    if rvec != nothing && lvec != nothing
        if v isa MPSBondTensor #normal transfer
            @tensor v[-1;-2]-=lvec[1,2]*v[2,1]*rvec[-1,-2]
        elseif v isa MPSTensor #utiity leg in the middle
            @tensor v[-1 -2;-3]-=lvec[1,2]*v[2,-2,1]*rvec[-1,-3]
        else
            @assert false
        end
    end

    return v
end
transfer_left(v,O::AbstractArray,A::AbstractArray,Ab::AbstractArray) = reduce((v,x)->transfer_left(v,x[1],x[2],x[3]),zip(O,A,Ab),init=v)
transfer_right(v,O::AbstractArray,A::AbstractArray,Ab::AbstractArray) = reduce((v,x)->transfer_right(v,x[1],x[2],x[3]),Iterators.reverse(zip(O,A,Ab)),init=v)


text:
struct CGMesh
    dim::Any
    cellmaps::Any
    nodalcoordinates::Any
    nodalconnectivity::Any
    cellconnectivity::Any
    numcells::Any
    numnodes::Any
    nodesperelement::Any
    x0::Any
    meshwidths::Any
    nelements::Any
    nfmside::Any
    mesh::Any
    elementsize::Any
end

function CGMesh(mesh, nodesperelement)
    dim = dimension(mesh)
    cellmaps = construct_cell_maps(mesh)
    x0 = reference_corner(mesh)
    meshwidths = mesh_widths(mesh)
    nelements = elements_per_mesh_side(mesh)
    elementsize = meshwidths ./ nelements
    nfeside = nodes_per_element_side(nodesperelement, dim)
    nfmside = nodes_per_mesh_side(nelements, nfeside)
    numnodes = prod(nfmside)

    nodalcoordinates =
        cg_nodal_coordinates(x0, mesh_widths(mesh), nelements, nfmside)
    nodalconnectivity =
        cg_nodal_connectivity(nfmside, nfeside, nodesperelement, nelements)
    cellconnectivity = cell_connectivity(mesh)
    numcells = number_of_elements(mesh)
    CGMesh(
        dim,
        cellmaps,
        nodalcoordinates,
        nodalconnectivity,
        cellconnectivity,
        numcells,
        numnodes,
        nodesperelement,
        x0,
        meshwidths,
        nelements,
        nfmside,
        mesh,
        elementsize,
    )
end

function CGMesh(x0, meshwidths, nelements, nodesperelement)
    mesh = UniformMesh(x0, meshwidths, nelements)
    return CGMesh(mesh, nodesperelement)
end

function Base.show(io::IO, mesh::CGMesh)
    dim = dimension(mesh)
    x0 = reference_corner(mesh)
    nelements = elements_per_mesh_side(mesh)
    meshwidths = mesh_widths(mesh)
    ncells = number_of_cells(mesh)
    nodesperelement = nodes_per_element(mesh)
    nnodes = number_of_nodes(mesh)
    nf = nodes_per_element(mesh)

    str =
        "CGMesh\n\tDimension : $dim\n\tCorner : $x0\n\tWidth : $meshwidths\n\t" *
        "Elements/Side : $nelements\n\tNum. Cells : $ncells\n\t" *
        "Nodes/Element : $nodesperelement\n\t" *
        "Num. Nodes : $nnodes"
    print(io, str)
end

function element_size(mesh::CGMesh)
    return mesh.elementsize
end

function dimension(mesh::CGMesh)
    return mesh.dim
end

function reference_corner(mesh::CGMesh)
    return mesh.x0
end

function mesh_widths(mesh::CGMesh)
    return mesh.meshwidths
end

function elements_per_mesh_side(mesh::CGMesh)
    return mesh.nelements
end

function number_of_cells(mesh::CGMesh)
    mesh.numcells
end

function number_of_nodes(mesh::CGMesh)
    return mesh.numnodes
end

function nodes_per_element(mesh::CGMesh)
    return mesh.nodesperelement
end

function nodal_connectivity(mesh::CGMesh)
    return mesh.nodalconnectivity
end

function nodal_connectivity(mesh::CGMesh, cellid::Int)
    return mesh.nodalconnectivity[:, cellid]
end

function nodes_per_mesh_side(mesh::CGMesh)
    return mesh.nfmside
end

function nodal_coordinates(mesh::CGMesh)
    return mesh.nodalcoordinates
end

function background_mesh(mesh::CGMesh)
    return mesh.mesh
end

function cell_map(mesh::CGMesh, cellid)
    return mesh.cellmaps[cellid]
end

function jacobian(mesh::CGMesh)
    return jacobian(cell_map(mesh, 1))
end

function inverse_jacobian(mesh::CGMesh)
    return inverse_jacobian(cell_map(mesh, 1))
end
################################################################################

function nodes_per_element_side(nodesperelement, dim)
    nfeside = (nodesperelement)^(1.0 / dim)
    @assert isinteger(nfeside)
    return round(Int, nfeside)
end

function nodes_per_mesh_side(nelements, nfeside)
    return (nfeside - 1) * nelements .+ 1
end

function cg_nodal_coordinates(x0, widths, nelements, nfmside)
    dim = length(x0)
    @assert length(x0) == length(widths) == length(nelements) == dim

    if dim == 1
        xrange = range(x0[1], stop = x0[1] + widths[1], length = nfmside[1])
        return Array(transpose(xrange))
    elseif dim == 2
        xrange = range(x0[1], stop = x0[1] + widths[1], length = nfmside[1])
        yrange = range(x0[2], stop = x0[2] + widths[2], length = nfmside[2])

        ycoords = repeat(yrange, outer = nfmside[1])
        xcoords = repeat(xrange, inner = nfmside[2])
        return vcat(xcoords', ycoords')
    else
        error("3D currently not supported")
    end
end

function cg_nodal_connectivity_2d(nfmside, nfeside, nodesperelement, nelements)
    numnodes = prod(nfmside)
    ncells = prod(nelements)

    nodeids = reshape(1:numnodes, reverse(nfmside)...)
    connectivity = zeros(Int, nodesperelement, ncells)
    colstart = 1
    colend = nfeside

    cellid = 1
    for col = 1:nelements[1]
        rowstart = 1
        rowend = nfeside
        for row = 1:nelements[2]
            connectivity[:, cellid] .=
                vec(nodeids[rowstart:rowend, colstart:colend])
            cellid += 1
            rowstart = rowend
            rowend += nfeside - 1
        end
        colstart = colend
        colend += nfeside - 1
    end
    return connectivity
end

function cg_nodal_connectivity_1d(nfmside, nodesperelement, nelements)
    numnodes = prod(nfmside)
    ncells = prod(nelements)

    connectivity = zeros(Int, nodesperelement, ncells)
    nodeids = 1:numnodes
    start = 1
    stop = nodesperelement

    for cellid = 1:ncells
        connectivity[:, cellid] = nodeids[start:stop]
        start = stop
        stop = start + nodesperelement - 1
    end
    return connectivity
end

function cg_nodal_connectivity(nfmside, nfeside, nodesperelement, nelements)
    dim = length(nfmside)
    @assert length(nfmside) == length(nelements)

    if dim == 1
        return cg_nodal_connectivity_1d(nfmside,nodesperelement,nelements)
    elseif dim == 2
        return cg_nodal_connectivity_2d(nfmside,nfeside,nodesperelement,nelements)
    else
        error("3D not supported")
    end
end


text:
export LocalCSRGraph, numRows, maxEntry, minEntry

"""
    LocalCSRGraph{EntriesType, IndexType}()
    LocalCSRGraph(entries::AbstractArray{EntriesType, 1}, rowMap::AbstractArray{IndexType, 1})

A compressed row storage array.  Used by CSRGraph to store local structure.
`EntriesType` is the type of the data being held
`IndexType` is the type used to represent the indices
"""
mutable struct LocalCSRGraph{EntriesType, IndexType <: Integer}
    entries::AbstractArray{EntriesType, 1}
    rowMap::AbstractArray{IndexType, 1}
end

function LocalCSRGraph{EntriesType, IndexType}() where{EntriesType, IndexType <: Integer}
    LocalCSRGraph(Array{EntriesType, 1}(undef, 0), Array{IndexType, 1}(undef, 0))
end

"""
    numRows(::LocalCSRGraph{EntriesType, IndexType})::IndexType

Gets the number of rows in the storage
"""
function numRows(graph::LocalCSRGraph{EntriesType, IndexType})::IndexType where {
        EntriesType, IndexType <: Integer}
    len = length(graph.rowMap)
    if len != 0
        len - 1
    else
        0
    end
end

"""
    maxEntry(::LocalCSRGraph{EntriesType})::EntriesType

Finds the entry with the maximum value.
"""
function maxEntry(graph::LocalCSRGraph{EntriesType})::EntriesType where {
        EntriesType}
    if length(graph.entries) != 0
        maximum(graph.entries)
    else
        throw(InvalidArgumentError("Cannot find the maximum of an empty graph"))
    end
end

"""
    minEntry(::LocalCSRGraph{EntriesType})::EntriesType

Finds the entry with the minimum value.
"""
function minEntry(graph::LocalCSRGraph{EntriesType})::EntriesType where {
        EntriesType}
    if length(graph.entries) != 0
        minimum(graph.entries)
    else
        throw(InvalidArgumentError("Cannot find the minimum of an empty graph"))
    end
end


text:
module Advection
  export AdvectionLaw

  import ..Atum
  using StaticArrays: SVector

  @Base.kwdef struct AdvectionLaw{FT, D} <: Atum.AbstractBalanceLaw{FT, D, 1}
    u⃗::SVector{D, FT} = ones(SVector{D, FT})
  end

  Atum.flux(law::AdvectionLaw, q, x⃗) = law.u⃗ * q'
  Atum.wavespeed(law::AdvectionLaw, n⃗, q, x⃗) = abs(n⃗' * law.u⃗)
end


text:
using JuMP
using Xpress

function reopt_model(MAXTIME)

    REopt = direct_model(Xpress.Optimizer(MAXTIME=-MAXTIME, OUTPUTLOG = 0))
	return REopt

end


text:
######################################################################################
#
# demo 6: Interactive brightness/contrast adjustment of live webcam with trackbars
#
######################################################################################

device_index = CAP_ANY

# C++ OpenCV code
cxx"""
void liveVideo_adjust(int device_index) {

    // open Video device
    cv::VideoCapture capture(device_index);

     //Check that it is indeed open
    if (!capture.isOpened())
    {
       std::cout<< "Can not open camera!";
       exit(0);
    }

    // Create a window to show the frames
    namedWindow("Welcome!", cv::WINDOW_AUTOSIZE);

    // Create trackbars for brightness and contrast
    int iValueForBrightness = 50;
    int iValueForContrast = 50;

    cv::createTrackbar("Brightness", "Welcome!", &iValueForBrightness, 100);
    cv::createTrackbar("Contrast", "Welcome!", &iValueForContrast, 100);

    // Set the frame size to the camera capture WIDTH and HEIGHT
    double dWidth = capture.get(cv::CAP_PROP_FRAME_WIDTH);
    double dHeight = capture.get(cv::CAP_PROP_FRAME_HEIGHT);

    int Height = static_cast<int>(dHeight);
    int Width = static_cast<int>(dWidth);

    // Initialize the cv::Mat structures
    cv::Mat frame(Width,Height, CV_8UC3);
    cv::Mat converted(Width,Height, CV_8UC3);

    // Loop until user presses ESC or frame is empty
    while(true)
    {
        bool Success = capture.read(frame);

        if (!Success) {
            break;
        }

        int iBrightness = iValueForBrightness -50;
        double dContrast = iValueForContrast/50.0;

        //cv::cvtColor(frame, frame, cv::COLOR_BGR2GRAY);
        //cv::equalizeHist(frame, converted);

        frame.convertTo(converted, -1, dContrast, iBrightness);

        cv::imshow("Welcome!", converted);

        if (cv::waitKey(30) == 27)
        {
            cv::destroyWindow("Welcome!");
            break;
        }
    }
}
"""

@cxx liveVideo_adjust(device_index)


text:
import FixedPointNumbers: N0f8

#defined in Plotting/palettes.jl
colormapdict = Dict{String,Vector{RGB{Float32}}}(

	 "Spectral"  =>  Spectral,
	 "RdYlGn"    =>  RdYlGn,
	 "GnYlRd"    =>  reverse(RdYlGn),
	 "RdBu"      =>  RdBu,
	 "PiYG"      =>  PiYG,
	 "PRGn"      =>  PRGn,
	 "RdYlBu"    =>  RdYlBu,
	 "BuYlRd"    =>  BuYlRd,#reverse(RdYlBu),
	 "BrBG"      =>  BrBG,
	 "RdGy"      =>  RdGy,
	 "GyRd"      =>  reverse(RdGy),
	 "PuOr"      =>  PuOr,
	 "Arnold"    =>  reverse(PiYG),
	 "BuRd"      =>  reverse(RdBu),
	 "GnYlRd"    =>  reverse(RdYlGn),
	 "Paired"    =>  Paired,
	 "Purples"   =>  Purples,
	 "BlPu"		 =>  BlPu,
	 "BkOr"   =>  BkOr,
	 "GrBl" => GrBl,
	 "distgclrs" => distgclrs,
	 "bone"		=>	bone,
	 "bone_rev"		=>	reverse(bone),
	 "YG" => YG)

abstract type AbstractColorMap
end

mutable struct ColorMap <: AbstractColorMap
	 map::Vector{RGB{Float32}}
	 function ColorMap(str::String)
		  clrmap = new()
		  clrmap.map = generate(colormapdict[str],1:256)
		  return clrmap
	 end
end

using Nullables

mutable struct BoundedColorMap
	 clrmap::ColorMap
	 mindta::Nullable{Float32}
	 maxdta::Nullable{Float32}
	 function BoundedColorMap(clrmap::ColorMap)
		  return new(clrmap,Nullable{Float32}(),Nullable{Float32}())
	 end
end

import Statistics.mean
function mean(clrs::Vector{RGB{Float32}})
	 r,g,b,a = zeros(Float32,4)
	 n = length(clrs)
	 for clr in clrs
		  r+=clr.r
		  g+=clr.g
		  b+=clr.b
	 end
	 r/=n;g/=n;b/=n
	 return RGB{Float32}(r,g,b)
end

function generate(colormap::Vector{RGB{Float32}},data::AbstractVector)
	 inds = map(x->Int(floor(x)),linspace(1,length(data),length(colormap)))
	 clrs = map(mean,map(x->colormap[find(y->y==x,inds)],unique(inds)))
	 inds[1] += sum(map(i->inds[i+1]-inds[i],1:(length(inds)-1)))-length(data)
	 return length(data) == length(clrs) ? clrs : vcat(map(i->(inds[i+1]-inds[i]) > 1 ? linspace(clrs[i],clrs[i+1],inds[i+1]-inds[i]) : clrs[i] ,1:(length(inds)-1))...)
end

mutable struct JetColorMap <: AbstractColorMap
	 map::Vector{RGBA{Float32}}
	 function JetColorMap()
		  cmap = Array(RGBA{Float32}, 256)
		  for i = 0:255
				n=4*i/256
				cmap[i+1] = RGBA(min(max(min(n-1.5,-n+4.5),0.),1.),
									 min(max(min(n-0.5,-n+3.5),0.),1.),
									 min(max(min(n+0.5,-n+2.5),0.),1.),1.)
		  end
		  return new(cmap)
	 end
end

mutable struct HSVColorMap <: AbstractColorMap
	 map::Vector{RGBA{Float32}}
	 function HSVColorMap()
		  return new(linspace(Colors.HSV(0,1,1),Colors.HSV(330,1,1),256))
	 end
end

function apply!(cm::BoundedColorMap, data::Union{Array{Float32,1},Array{Float64,1},AbstractVector{Int}})
	 cm.mindta = minimum(data)
	 cm.maxdta = maximum(data)
	 return apply(cm.clrmap,data)
end

function apply(cm::AbstractColorMap,data::Union{Array{Float64,1},Array{Float32,1}})
	 color = Array{RGBA{Float32}}(length(data))
	 mindta = minimum(data)
	 maxdta = maximum(data)
	 interv = maxdta-mindta
	 if isapprox(abs(interv),0.,atol=1e-10)
		  fill!(color,cm.map[floor(Int,length(cm.map)/2)])
	 else
		  pos = 0.0
		  for i = 1:length(data)
				pos = (data[i]-mindta)/interv * 255.0 + 1.0
				clr_1 = cm.map[floor(Int64,pos)]
				clr_2 = cm.map[ceil(Int64,pos)]
				weight = pos%1
				clr = weighted_color_mean(weight,clr_2,clr_1)
				#clr = (1-weight)*clr_1+weight*clr_2
				color[i] = clr
		  end
	 end
	 return color
end

function apply(cm::AbstractColorMap,data::AbstractVector{Int})
	 color = Array{RGBA{Float32}}(length(data))
	 mindta = minimum(data)
	 maxdta = maximum(data)
	 interv = maxdta-mindta
	 if interv < 1
		  fill!(color,cm.map[floor(Int,length(cm.map)/2)])
	 else
		  pos = 0.0
		  for i = 1:length(data)
				pos = (data[i]-mindta)/interv * 255.0 + 1.0
				clr_1 = cm.map[floor(Int64,pos)]
				clr_2 = cm.map[ceil(Int64,pos)]
				weight = pos%1
				clr = weighted_color_mean(weight,clr_2,clr_1)
				#clr = (1-weight)*clr_1+weight*clr_2
				color[i] = clr
		  end
	 end
	 return color
end


text:
#=

=#

"""
    Returns an array containing the degree for each vertex
    in the given graph.

    Args:
        graph:
            'Graph' structure where the vertices and edges are placed.
        mode:
            {'total', 'in', 'out'} string to signal which type of
            degree to return.
    Return:
        degree:
            Integer representing the value of the degree.
"""
function get_degree(graph::Graph, mode::String)
    degree = zeros(Int, length(graph.vertices))
    for j in 1:length(graph.vertices)
        vertex = graph.vertices[j]
        if mode=="total"
            degree[j] = length(vertex.edges_source) + length(vertex.edges_target)
        elseif mode=="in"
            degree[j] = length(vertex.edges_target)
        elseif mode=="out"
            degree[j] = length(vertex.edges_source)
        end
    end
    return degree
end

"""
    Returns an array containing the indexes of each incoming
    neighbor of 'node' in 'graph'.

    Args:
        node:
            Integer representing the index of the node in the
            parsed graph for which we want to extract the neighbors.
        graph:
            'Graph' structure where the vertices and edges are placed.
    Return:
        rtn_neighbors:
            Array of integers representing all the indexes of the incoming
            neighbors.
"""
function get_in_neighbors(node::Int, graph::Graph)
    node_obj = graph.vertices[node]

    # -- edges where 'node' is the target --
    incoming_edges = node_obj.edges_target
    incoming_neighbors = Int[]
    for edge in incoming_edges
        # fundamental check
        if edge.target!=node
            print("BUG. aborted.\n")
            return
        end
        push!(incoming_neighbors, edge.source)
    end
    rtn_neighbors = collect(Set(incoming_neighbors))
    return rtn_neighbors
end

"""
    Returns an array containing the indexes of each outgoing
    neighbor of 'node' in 'graph'.

    Args:
        node:
            Integer representing the index of the node in the
            parsed graph for which we want to extract the neighbors.
        graph:
            'Graph' structure where the vertices and edges are placed.
    Return:
        rtn_neighbors:
            Array of integers representing all the indexes of the outgoing
            neighbors.
"""
function get_out_neighbors(node::Int, graph::Graph)
    node_obj = graph.vertices[node]

    # -- edges where 'node' is the source --
    outgoing_edges = node_obj.edges_source
    outgoing_neighbors = Int[]
    for edge in outgoing_edges
        # fundamental check
        if edge.source!=node
            print("BUG. aborted.\n")
        end
        push!(outgoing_neighbors, edge.target)
    end
    rtn_neighbors = collect(Set(outgoing_neighbors))
    return rtn_neighbors
end

"""
    Returns an array containing the indexes of each neighbor (in and out)
    of 'node' in 'graph'.

    Args:
        node:
            Integer representing the index of the node in the
            parsed graph for which we want to extract the neighbors.
        graph:
            'Graph' structure where the vertices and edges are placed.
    Return:
        rtn_neighbors:
            Array of integers representing all indexes of all neighbors.
"""
function get_all_neighbors(node::Int, graph::Graph)
    incoming_neighbors = get_in_neighbors(node, graph)
    outcoming_neighbors = get_out_neighbors(node, graph)

    neighbors = Int[]
    append!(neighbors, incoming_neighbors)
    append!(neighbors, outcoming_neighbors)
    rtn_neighbors = collect(Set(neighbors))
    return rtn_neighbors
end

"""
    Breadth-first search from 'source' in 'graph'. Returns three arrays:
    'color', 'dist' and 'parent'.

    Args:
        source:
            Integer index of the source node of the BFS algorithm.
        graph:
            'Graph' structure where the vertices and edges are placed.
    Return:
        color:
            visited status of the node (1=visited)
        dist:
            shortest distance between 'source' and the current node.
        parent:
            parent of each node in the tree from 'source'.
"""
function bfs_search(source::Int, graph::Graph)
    N = length(graph.vertices)
    color = [-1 for j in 1:N]
    dist = [-1 for j in 1:N]
    parent = [-1 for j in 1:N]

    # -1/0/1 -> white/gray/black
    color[source] = 0
    dist[source] = 0
    parent[source] = -1

    queue = Int[]
    push!(queue, source)
    while length(queue)>0
        u = pop!(queue)
        u_adj = get_out_neighbors(u, graph)
        for w in u_adj
            if color[w]==-1
                color[w] = 0
                dist[w] = dist[u]+1
                parent[w] = u
                push!(queue, w)
            end
        end
        color[u] = 1
    end
    return color, dist, parent
end

function dfs_search(graph::Graph)
    N = length(graph.vertices)
    color = [-1 for j in 1:N ]
    dist = [-1 for j in 1:N ]
    parent = [-1 for j in 1:N ]
    finished = [-1 for j in 1:N ]

    time = [0]
    for u in 1:N
        if color[u]==-1
            dfs_visit(u, graph, time, color, dist, parent, finished)
        end
    end
    return color, parent, finished
end

"""
    Invert the direction of all edges of 'graph'.
"""
function transpose_graph(graph::Graph)
    edges = graph.edges
    for edge in edges
        src = edge.source
        tgt = edge.target
        edge.source = tgt
        edge.target = src
    end
    for node in graph.vertices
        aux_src = node.edges_source
        aux_tgt = node.edges_target
        node.edges_source = aux_tgt
        node.edges_target = aux_src
    end
end

function get_root(node::Int, parent::Array{Int,1})
    r = node
    while parent[r]!=-1
        r = parent[r]
    end
    return r
end

function dfs_visit(u::Int, graph::Graph, time::Array{Int},
                   color::Array{Int}, dist::Array{Int},
                   parent::Array{Int}, finished::Array{Int})
    time[1] += 1
    dist[u] = time[1]
    color[u] = 0
    u_adj = get_out_neighbors(u, graph)
    for v in u_adj
        if color[v]==-1
            parent[v] = u
            dfs_visit(v, graph, time, color, dist, parent, finished)
        end
    end
    color[u] = 1
    time[1] += 1
    finished[u] = time[1]
end

"""
    Extract the partition of strongly connected components of the graph.
    Obtained using the Kosajaru algorithm.

    Args:
        graph:
            'Graph' structure where the vertices and edges are placed.
        return_dict:
            {default=false} to signal whether the function should return
            a formatted dictionary where keys are root nodes of components
            and values are arrays containing the nodes belonging to this
            strongly connected component.
    Return:
        node_labels:
            Array of integers with length equal to the total number of nodes.
            The jth index of this array represent the index of the SCC in which
            the node j belongs to.
        unique_labels:
            Array containing all the unique labels for each SCC found.
        scc_trees(return_dict=true):
            Dictionary where keys are SCC labels and values are arrays containing
            all the nodes belonging to the current SCC label.

"""
function extract_strong(graph::Graph, return_dict=false)
    N = graph.N
    color, parent, finished = dfs_search(graph)

    # Create the tranpose graph from 'graph'.
    graph_t = copy_graph(graph)
    transpose_graph(graph_t)

    time = [0]
    color_t = [-1 for j in 1:N ]
    dist_t = [-1 for j in 1:N ]
    parent_t = [-1 for j in 1:N ]
    finished_t = [-1 for j in 1:N ]
    # Apply second DFS to 'graph_t' in decreasing order of 'finished'.
    node_ordering = sortperm(finished, rev=true)
    for u in node_ordering
        if color_t[u]==-1
            dfs_visit(u, graph_t, time, color_t, dist_t, parent_t, finished_t)
        end
    end

    # Now, each DFS tree in 'parent_t' represents an strongly connected component.
    scc_trees = Dict{Int, Array{Int}}()
    node_labels = [-1 for j in 1:N]
    for u in 1:N
        root = get_root(u, parent_t)
        if get(scc_trees, root, -1)==-1
            scc_trees[root] = Int[]
        end
        push!(scc_trees[root], u)
        node_labels[u] = root
    end
    unique_labels = collect(Int, Set(node_labels))
    if return_dict
        return node_labels, unique_labels, scc_trees
    end
    return node_labels, unique_labels
end


text:
using MarkovModelAndMMs.MarkovStuff
using MarkovModelAndMMs.MMs
using MarkovModelAndMMs.CostFunctions
using LinearAlgebra
using Test

let m = MarkovProcess([-1//2 1//4 1//4; 3//4 -1//1 1//4; 1//2 1//6 -2//3]),
	c = MarkovChain([0 1//2 1//2; 3//4 0 1//4; 3//4 1//4 0]),
	λ = [1//2, 1//1, 2//3]
	@test MarkovProcess(c, λ).G == m.G
	@test imbeddedchain(m).P == c.P
end


text:
### A Pluto.jl notebook ###
# v0.11.14

using Markdown
using InteractiveUtils

# ╔═╡ a31e1a3a-f637-11ea-1e0f-2d61b88049db
using Suppressor, InteractiveUtils

# ╔═╡ 80b2135e-f1ed-11ea-17ee-656b8062026a
md"""
# What are macros?

**Macros** change the syntax of parts of a program.

While functions work on values, a macro works on a **piece of Julia code**, which is represented as a data structure in Julia itself, and **rewrites it to a different piece of code**! This can be a very powerful thing to do.
"""

# ╔═╡ 7acec72a-f1fa-11ea-2753-efcf2cee8335
md"""As a simple example, let's look at the macro called `@elapsed`, which measures how much time it takes to run a piece of code.

It does this by taking the input expression and adding extra code around it.

First let's see it in action! We'll use a function that takes some time to run: the function `peakflops()` estimates how fast your computer runs. (It reports an estimate of the number of flops, i.e. floating-point operations per second, when doing matrix multiplications.)
"""

# ╔═╡ 0f7c908c-f1fb-11ea-18fb-c7e70278ec46
peakflops()

# ╔═╡ f43139a2-f1ed-11ea-0cdf-678ffac4374a
@elapsed peakflops()

# ╔═╡ 67b3fb78-f1ee-11ea-2d0a-a92e81580917
md"""
**NOTE:** a macro is always called using the `@` prefix character; this is a visual cue that the following piece of code is *not what is actually executed*, but rather is *transformed in some way*.
"""

# ╔═╡ 5677f2d6-f1ee-11ea-3f9b-b59c8a8d07ff
md"""
You can look at the transformed expression by using `@macroexpand`. Since this name begins with `@`, we can see that this is actually itself a macro!
"""

# ╔═╡ 4b2cbfc4-f1ee-11ea-17d6-a72fd22e6d89
@macroexpand @elapsed peakflops()

# ╔═╡ c4e4e402-f1fa-11ea-3e5e-4344dcc7e17b
md"Let's remove the line number information:"

# ╔═╡ 9639e822-f1ee-11ea-19b6-5f377bb4b809
Base.remove_linenums!( @macroexpand @elapsed peakflops() )

# ╔═╡ ea65d68c-f1fa-11ea-1cde-ebf4b7958b8a
md"We see that the macro has added more code around the original `peakflops()` call.
Julia *replaces* the original code by the new code; this new code is what will actually be used (compiled)."

# ╔═╡ a5816132-f1ee-11ea-268b-05086212324b
md"""
Here, the `remove_linenums!` function actually receives a Julia **expression** as its argument.
"""

# ╔═╡ e3bd8fbc-f1fa-11ea-331f-bbe8f038f46a
md"## Expressions in Julia"

# ╔═╡ cb5e6a8a-f1fb-11ea-2f15-e7486afdfa4a
md"One feature that makes Julia powerful is its ability (influenced by the Lisp language) to represent Julia code as Julia objects, which can be manipulated using the usual Julia features.
"

# ╔═╡ faa18946-f1fb-11ea-0f5b-290683f29ebf
md"Let's look at some simple expressions.

If we type an expression, Julia will **evaluate** it and return its value:
"

# ╔═╡ c5465fe0-f1ee-11ea-2f61-3172509f693d
x = 1 + 2

# ╔═╡ 6bb8d276-f1ff-11ea-2ada-1916761bc6a8
md"The variable `x` now contains the value `3`, the result of evaluating the expression `1 + 2`."

# ╔═╡ 10d909b4-f1fc-11ea-080a-f9f87b3512c5
md"We can prevent Julia from evaluating an expression by **quoting** it, which we can do in one of two ways:"

# ╔═╡ c88dfe88-f1ee-11ea-0327-8b18a25c511c
expr = :(1 + 2)

# ╔═╡ 80ec2c30-f1ff-11ea-30ea-ef5dd01e26ae
expr2 = quote
			1 + 2
		end

# ╔═╡ 60778b00-f1ff-11ea-011f-231f3ab18de1
md"This *leaves the expression alone, without evaluating it*, giving as result the unevaluated expression. This is a Julia object of type `Expr` that we have not seen before:"

# ╔═╡ bdb08e52-f1ff-11ea-1881-4bb2bb3f5c77
typeof(expr)

# ╔═╡ 6a0fdfc8-f637-11ea-04d8-a55a4c42df3e
md"As an example, we can create the following expression:"

# ╔═╡ cbce40ee-f1ee-11ea-2e05-1d523acbf8b9
expr3 = :(y + 1)

# ╔═╡ 7f1cad88-f637-11ea-1cf2-4d5b9a4fb50f
md"If we evaluate this expression, it gives an error, since there is no variable `y` in the current notebook:"

# ╔═╡ 906e3234-f637-11ea-1f77-679dc897ec69
y + 1

# ╔═╡ 940f888c-f637-11ea-20a0-1ddda8fb7908
md"Since expressions are Julia objects, we can examine them using Julia's usual mechanisms. For example, the `dump` function allows us to inspect the fields inside an object:"

# ╔═╡ 8dab58d4-f639-11ea-1b0b-536927639c6c
md"We can extract the fields inside the object:"

# ╔═╡ e46e6c96-f1ee-11ea-3e7f-4f613282ccb1
expr3.head

# ╔═╡ e631bede-f1ee-11ea-0fab-7516ab535cf5
expr3.args

# ╔═╡ 95cf0d76-f639-11ea-34dc-f32852ac32b8
md"and even modify them:"

# ╔═╡ 9a587ba4-f639-11ea-2bef-37fed4d45a7f
expr3.args[2]

# ╔═╡ 9cd4334e-f639-11ea-26c6-e9782451ef02
expr3.args[2] = :x

# ╔═╡ a037391e-f639-11ea-2723-455ccff4b024
md"(`:x` and `:y` here are `Symbol`s, a special type of string.)"

# ╔═╡ b6e67c1a-f639-11ea-3abe-f75671ed4ca2
expr3

# ╔═╡ cb0cb434-f639-11ea-1ea2-93009c372885
md"We see that we have been able to *modify the `Expr` object, which is a piece of Julia code, using Julia itself! This is the simplest example of **metaprogramming**, i.e. writing a program that modifies another *program* (instead of modifying data).

This is what macros do internally: transform Julia expressions into something else.
"

# ╔═╡ 0ff6e154-f1ef-11ea-2ee8-61eb98523642
md"""
## Conclusion

When we see something like `@foo`, we recognise that it's a **syntactic transformation** and may modify the behavior of the language in some way!

Being able to represent Julia code within Julia is very useful, and hence we frequently see macros and expression manipulation being used in Julia code.
"""

# ╔═╡ aa3a731c-f637-11ea-3afc-9144bdeb7d6f
md"## Function library"

# ╔═╡ b440b66c-f637-11ea-393d-9d444aa58501
function with_terminal(f)
	local spam_out, spam_err
	@color_output false begin
		spam_out = @capture_out begin
			spam_err = @capture_err begin
				f()
			end
		end
	end
	spam_out, spam_err

	HTML("""

$(Markdown.htmlesc(spam_out))

	""")
end

# ╔═╡ 77437eaa-f639-11ea-0c4d-3756cbd589da
with_terminal() do
	dump(expr3)
end

# ╔═╡ Cell order:
# ╟─80b2135e-f1ed-11ea-17ee-656b8062026a
# ╟─7acec72a-f1fa-11ea-2753-efcf2cee8335
# ╠═0f7c908c-f1fb-11ea-18fb-c7e70278ec46
# ╠═f43139a2-f1ed-11ea-0cdf-678ffac4374a
# ╟─67b3fb78-f1ee-11ea-2d0a-a92e81580917
# ╟─5677f2d6-f1ee-11ea-3f9b-b59c8a8d07ff
# ╠═4b2cbfc4-f1ee-11ea-17d6-a72fd22e6d89
# ╟─c4e4e402-f1fa-11ea-3e5e-4344dcc7e17b
# ╠═9639e822-f1ee-11ea-19b6-5f377bb4b809
# ╟─ea65d68c-f1fa-11ea-1cde-ebf4b7958b8a
# ╟─a5816132-f1ee-11ea-268b-05086212324b
# ╟─e3bd8fbc-f1fa-11ea-331f-bbe8f038f46a
# ╟─cb5e6a8a-f1fb-11ea-2f15-e7486afdfa4a
# ╟─faa18946-f1fb-11ea-0f5b-290683f29ebf
# ╠═c5465fe0-f1ee-11ea-2f61-3172509f693d
# ╟─6bb8d276-f1ff-11ea-2ada-1916761bc6a8
# ╟─10d909b4-f1fc-11ea-080a-f9f87b3512c5
# ╠═c88dfe88-f1ee-11ea-0327-8b18a25c511c
# ╠═80ec2c30-f1ff-11ea-30ea-ef5dd01e26ae
# ╟─60778b00-f1ff-11ea-011f-231f3ab18de1
# ╠═bdb08e52-f1ff-11ea-1881-4bb2bb3f5c77
# ╟─6a0fdfc8-f637-11ea-04d8-a55a4c42df3e
# ╠═cbce40ee-f1ee-11ea-2e05-1d523acbf8b9
# ╟─7f1cad88-f637-11ea-1cf2-4d5b9a4fb50f
# ╠═906e3234-f637-11ea-1f77-679dc897ec69
# ╟─940f888c-f637-11ea-20a0-1ddda8fb7908
# ╠═77437eaa-f639-11ea-0c4d-3756cbd589da
# ╠═8dab58d4-f639-11ea-1b0b-536927639c6c
# ╠═e46e6c96-f1ee-11ea-3e7f-4f613282ccb1
# ╠═e631bede-f1ee-11ea-0fab-7516ab535cf5
# ╠═95cf0d76-f639-11ea-34dc-f32852ac32b8
# ╠═9a587ba4-f639-11ea-2bef-37fed4d45a7f
# ╠═9cd4334e-f639-11ea-26c6-e9782451ef02
# ╠═a037391e-f639-11ea-2723-455ccff4b024
# ╠═b6e67c1a-f639-11ea-3abe-f75671ed4ca2
# ╟─cb0cb434-f639-11ea-1ea2-93009c372885
# ╟─0ff6e154-f1ef-11ea-2ee8-61eb98523642
# ╟─aa3a731c-f637-11ea-3afc-9144bdeb7d6f
# ╠═a31e1a3a-f637-11ea-1e0f-2d61b88049db
# ╠═b440b66c-f637-11ea-393d-9d444aa58501

max_stars_repo_name: mathijsvdv/ComputationalThinking

max_issues_repo_name: mathijsvdv/ComputationalThinking

max_forks_repo_name: mathijsvdv/ComputationalThinking


text:
using PowerAnalysis
using Documenter

makedocs(;
    modules=[PowerAnalysis],
    authors="John Myles White",
    repo="https://github.com/johnmyleswhite/PowerAnalysis.jl/blob/{commit}{path}#L{line}",
    sitename="PowerAnalysis.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://johnmyleswhite.github.io/PowerAnalysis.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/johnmyleswhite/PowerAnalysis.jl",
)


text:
using Test, LinearAlgebra, OrderedCollections
using WaveProp
using WaveProp.Nystrom
using WaveProp.Geometry
using WaveProp.Integration
using WaveProp.Mesh

@testset "Single element" begin
    pde  = op = Helmholtz(dim=2,k=1)
    el          = LagrangeLine((1,0),(1.1,0))
    q   = GaussLegendre(20)
    x̂,ŵ = q()
    yi  = map(el,x̂)
    νi  = map(x->normal(el,x),x̂)
    nb = 3*length(yi)
    G   = SingleLayerKernel(op)
    dG  = DoubleLayerKernel(op)
    ##
    ϕ  = y->cos(y[1])  # regular density
    x  = el(0.51)
    ν  = normal(el,0.51)
    w  = Nystrom.singular_weights_ldim(G,el,yi,νi,x,:onsurface)
    Ie  = integrate(y->G(x,y)*ϕ(y),el)
    Ia = w*ϕ.(yi)
    @test norm(Ie - Ia,Inf)/norm(Ie) < 1e-5
    y₀  = el(0.51)
    x   = x + 1e-2*ν # outside
    w  = Nystrom.singular_weights_ldim(G,el,yi,νi,x,:outside)
    Ie  = integrate(y->G(x,y)*ϕ(y),el)
    Ia = w*ϕ.(yi)
    @test norm(Ie - Ia,Inf)/norm(Ie) < 1e-5
    x   = x - 1e-2*ν # inside
    w  = Nystrom.singular_weights_ldim(G,el,yi,νi,x,:inside)
    Ie  = integrate(y->G(x,y)*ϕ(y),el)
    Ia = w*ϕ.(yi)
    @test norm(Ie - Ia,Inf)/norm(Ie) < 1e-5
end

@testset "Single obstacle" begin
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    # geo  = Circle()
    Geometry.clear!()
    Ω,M   = WaveProp.IO.gmsh_disk()
    Γ     = boundary(Ω)
    # generate a Nystrom mesh with Gauss-Legendre quadrature
    qrule   = GaussLegendre(3)
    e2qrule = OrderedDict(E=>qrule for E in etypes(M))
    mesh    = NystromMesh(M,Γ,e2qrule)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    Sldim = Nystrom.assemble_ldim(S)
    Dldim = Nystrom.assemble_ldim(D)
    e1    = WaveProp.Nystrom.error_interior_green_identity(Sldim,Dldim,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
    K     = AdjointDoubleLayerOperator(pde,mesh)
    H     = HyperSingularOperator(pde,mesh)
    e0    = WaveProp.Nystrom.error_interior_derivative_green_identity(K,H,γ₀u,γ₁u)
    Kldim  = Nystrom.assemble_ldim(K)
    Hldim  = Nystrom.assemble_ldim(H)
    e1    = WaveProp.Nystrom.error_interior_derivative_green_identity(Kldim,Hldim,γ₀u,γ₁u)
    @test 10*norm(e1,Inf) < norm(e0,Inf)
    @test norm(e1,Inf) < 1e-4
end

@testset "Single obstacle" begin
    pde  = Helmholtz(dim=2,k=1)
    xout = SVector(-10,0)
    u    = (x)   -> SingleLayerKernel(pde)(xout,x)
    dudn = (x,n) -> DoubleLayerKernel(pde)(xout,x,n)
    Geometry.clear!()
    geo  = Kite()
    Ω    = Domain(geo)
    M   = meshgen(Ω,h=0.05)
    Γ     = boundary(Ω)
    # generate a Nystrom mesh with Gauss-Legendre quadrature
    qrule   = GaussLegendre(5)
    e2qrule = OrderedDict(E=>qrule for E in etypes(M))
    mesh    = NystromMesh(M,Γ,e2qrule)
    γ₀u   = γ₀(u,mesh)
    γ₁u   = γ₁(dudn,mesh)
    S     = SingleLayerOperator(pde,mesh)
    D     = DoubleLayerOperator(pde,mesh)
    e0    = WaveProp.Nystrom.error_interior_green_identity(S,D,γ₀u,γ₁u)
    Sldim = Nystrom.assemble_ldim(S)
    Dldim = Nystrom.assemble_ldim(D)
    e1    = WaveProp.Nystrom.error_interior_green_identity(Sldim,Dldim,γ₀u,γ₁u)
    @test norm(e0,Inf) > 1e-5
    @test norm(e1,Inf) < 1e-5
end

max_stars_repo_name: WaveProp/WaveProp

max_issues_repo_name: WaveProp/WaveProp

max_forks_repo_name: WaveProp/WaveProp


text:
# This file includes testset for Simulation

@testset "SimulationTestSet" begin
    @info "Running SimulationTestSet ..."

    # Simulation construction
    model = Model()
    clock = Clock(0 : 0.01 : 1.)
    simname = string(uuid4())
    simdir = tempdir()
    sim = Simulation(model, clock, simdir=simdir, simname=simname, logger=SimpleLogger())
    @test sim.model === model
    @test startswith(basename(sim.path), "Simulation-")
    @test sim.path == joinpath(simdir, "Simulation-" * simname)
    @test sim.state == :idle
    @test sim.retcode == :unknown
    @test sim.name == "Simulation-" * simname

    # Check Writer files
    model = Model()
    clock = Clock(0. : 0.01 : 1.)
    addnode!(model, SinewaveGenerator(), label=:gen)
    addnode!(model, Writer(), label=:writer)
    addbranch!(model, :gen => :writer)
    dname1 = dirname(getnode(model, :writer).component.file.path)
    simname = string(uuid4())
    simdir = tempdir()
    sim = Simulation(model, clock, simdir=simdir, simname=simname)
    @test dirname(getnode(model, :writer).component.file.path) == sim.path
    @test dname1 != sim.path

    # Report Simulation
    sim = simulate!(model, 0., 0.01, 10.)
    report(sim)
    filename = joinpath(sim.path, "report.jld2")
    @test isfile(filename)
    data = load(filename)
    @test data["name"] == sim.name
    @test data["path"] == sim.path
    @test data["state"] == sim.state
    @test data["retcode"] == sim.retcode

    @info "Done SimulationTestSet."
end # testset


text:
"""
Get the conference schedule as a nested JSON object.
On first call, the schedule is downloaded from Pretalx and cached for further usage.
"""
function get_conference_json()
    if !isassigned(conf_json)
        data = urldownload(CONFERENCE_SCHEDULE_JSON_URL)
        try
            conf_json[] = data.schedule.conference
        catch err
            error("Invalid conference schedule JSON file.")
            println(err)
        end
    end

    return conf_json[]
end

"""
Given a room (i.e. a fixed day and track), it finds the talks that are running now (it only compares times).
"""
function _find_current_talk_in_room(room::JSON3.Array; now=Dates.now())
    for talk in room
        start_time = Time(talk.start)
        dur = Time(talk.duration)
        end_time = start_time + Hour(dur) + Minute(dur)
        if start_time <= Time(now) <= end_time
            return talk
        end
    end
    return nothing
end

"""
Given a fixed day, it finds the talks in all tracks / rooms that are running now (it only compares times).

Returns a vector of tuples of the type `(room::Symbol, talk::JSON3.Object)`.
"""
function _find_current_talks_on_day(d::JSON3.Object; now=Dates.now())
    query_result = [(r, _find_current_talk_in_room(d.rooms[r]; now=now)) for r in keys(d.rooms)]
    return filter(x->!isnothing(x[2]), query_result)
end

# now_fake = DateTime("2020-07-29T16:30:00.000")

function get_running_talks(; now=Dates.now())
    conf = get_conference_json()
    days = Date[Date(d.date) for d in conf.days]

    dayidx = findfirst(isequal(Date(now)), days)
    if isnothing(dayidx)
        @info "There is no JuliaCon program today!"
        return nothing
    end

    d = conf.days[dayidx]
    if !(DateTime(d.day_start[1:end-6]) <= now <= DateTime(d.day_end[1:end-6]))
        @info "There is no JuliaCon program now!"
        return nothing
    end

    current_talks = _find_current_talks_on_day(d; now=now)
    return current_talks
end

function _track2color(track::Symbol)
    if track == Symbol("Red Track")
        return :red
    elseif track == Symbol("Green Track")
        return :green
    elseif track == Symbol("Purple Track")
        return :magenta
    else
        return :default
    end
end

function _print_running_talks(current_talks; now=Dates.now())
    !isnothing(current_talks) || return nothing
    # println()
    # println(Dates.format(Dates.now(), "HH:MM dd-mm-YYYY"))
    for (track, talk) in current_talks
        println()
        printstyled(track, bold=true, color=_track2color(track))
        println()
        println("\t", talk.title, " (", talk.type,")")
        println("\t", "└─ ", talk.url)
    end
    println("\n")
    println("(Full schedule: https://pretalx.com/juliacon2021/schedule)")
    return nothing
end

function now(; now=Dates.now())
    current_talks = get_running_talks(; now=now)
    _print_running_talks(current_talks; now=now)
    return nothing
end

function get_today(; now=Dates.now())
    conf = get_conference_json()
    days = [Date(d.date) for d in conf.days]

    dayidx = findfirst(isequal(Date(now)), days)
    if isnothing(dayidx)
        @info "There is no JuliaCon program today!"
        return nothing
    end

    d = conf.days[dayidx]
    schedule = Vector{Tuple{String, Matrix{String}}}(undef, length(d.rooms))
    i = 1
    for (track, talks) in d.rooms
        timetable = Matrix{String}(undef, length(talks), 4)
        for (j, talk) in enumerate(talks)
            timetable[j,1] = talk.start
            timetable[j,2] = talk.title
            timetable[j,3] = talk.type
            timetable[j,4] = talk.duration
        end
        schedule[i] = (string(track), timetable)
        i+=1
    end
    return schedule
end

function today(; now=Dates.now(), track=nothing)
    track_schedules = get_today(; now=now)
    isnothing(track_schedules) && return nothing
    header = (["Time", "Title", "Type"],)
    header_crayon = crayon"dark_gray bold"
    border_crayon = crayon"dark_gray"
    h_times = Highlighter((data, i, j) -> j == 1, crayon"white bold")
    for (tr, sched) in track_schedules
        !isnothing(track) && tr != track && continue
        h_current = _get_current_talk_highlighter(sched; now=now)
        println()
        pretty_table(@view sched[:,1:3];
            title = tr,
            title_crayon = Crayon(foreground = _track2color(Symbol(tr)), bold = true),
            header = header,
            header_crayon = header_crayon,
            border_crayon = border_crayon,
            highlighters = (h_times, h_current),
            tf = tf_unicode_rounded,
            alignment = :l,
        )
    end
    println()
    println("Check out https://pretalx.com/juliacon2021/schedule for more information.")
    println()
    printstyled("(Currently running talks are highlighted in ")
    printstyled("yellow", color=:yellow)
    printstyled(".)")
    println()
    return nothing
end

function _get_current_talk_highlighter(sched; now=Dates.now())
    for (i, talk) in enumerate(eachrow(sched))
        start_time = Time(talk[1])
        dur = Time(talk[4])
        end_time = start_time + Hour(dur) + Minute(dur)
        if start_time <= Time(now) <= end_time
            return Highlighter((data, m, n) -> m == i, crayon"yellow")
        end
    end
    return nothing
end