text:
@doc raw"""
    SymmetricPositiveSemidefiniteFixedRank{n,k,𝔽} <: AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType}

The [`Manifold`](@ref) $ \operatorname{SPS}_k(n)$ consisting of the real- or complex-valued
symmetric positive semidefinite matrices of size $n × n$ and rank $k$, i.e. the set

````math
\operatorname{SPS}_k(n) = \bigl\{
p  ∈ 𝔽^{n × n}\ \big|\ p^{\mathrm{H}} = p,
apa^{\mathrm{H}} \geq 0 \text{ for all } a ∈ 𝔽
\text{ and } \operatorname{rank}(p) = k\bigr\},
````
where $\cdot^{\mathrm{H}}$ denotes the Hermitian, i.e. complex conjugate transpose,
and the field $𝔽 ∈ \{ ℝ, ℂ\}$.
We sometimes $\operatorname{SPS}_{k,𝔽}(n)$, when distinguishing the real- and complex-valued
manifold is important.

An element is represented by $q ∈ 𝔽^{n × k}$ from the factorization $p = qq^{\mathrm{H}}$.
Note that since for any unitary (orthogonal) $A ∈ 𝔽^{n × n}$ we have
$(Aq)(Aq)^{\mathrm{H}} = qq^{\mathrm{H}} = p$, the representation is not unique, or in
other words, the manifold is a quotient manifold of $𝔽^{n × k}$.

The tangent space at $p$, $T_p\operatorname{SPS}_k(n)$, is also represented
by matrices $Y ∈ 𝔽^{n × k}$ and reads as

````math
T_p\operatorname{SPS}_k(n) = \bigl\{
X ∈ 𝔽^{n × n}\,|\,X = qY^{\mathrm{H}} + Yq^{\mathrm{H}}
\text{ i.e. } X = X^{\mathrm{H}}
\bigr\}.
````

Note that the metric used yields a non-complete manifold.
The metric was used in[^JourneeBachAbsilSepulchre2010][^MassartAbsil2020].

# Constructor

    SymmetricPositiveSemidefiniteFixedRank(n::Int, k::Int, field::AbstractNumbers=ℝ)

Generate the manifold of $n × n$ symmetric positive semidefinite matrices of rank $k$
over the `field` of real numbers `ℝ` or complex numbers `ℂ`.

[^JourneeBachAbsilSepulchre2010]:
    > Journée, M., Bach, F., Absil, P.-A., and Sepulchre, R.:
    > “Low-Rank Optimization on the Cone of Positive Semidefinite Matrices”,
    > SIAM Journal on Optimization (20)5, pp. 2327–2351, 2010.
    > doi: [10.1137/080731359](https://doi.org/10.1137/080731359),
    > arXiv: [0807.4423](http://arxiv.org/abs/0807.4423).
[^MassartAbsil2020]:
    > Massart, E., Absil, P.-A.:
    > "Quotient Geometry with Simple Geodesics for the Manifold of Fixed-Rank Positive-Semidefinite Matrices",
    > SIAM Journal on Matrix Analysis and Applications (41)1, pp. 171–198, 2020.
    > doi: [10.1137/18m1231389](https://doi.org/10.1137/18m1231389),
    > preprint: [sites.uclouvain.be/absil/2018.06](https://sites.uclouvain.be/absil/2018.06).
"""
struct SymmetricPositiveSemidefiniteFixedRank{n,k,𝔽} <:
       AbstractEmbeddedManifold{𝔽,DefaultIsometricEmbeddingType} end

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

@doc raw"""
    check_manifold_point(M::SymmetricPositiveSemidefiniteFixedRank{n,𝔽}, q; kwargs...)

Check whether `q` is a valid manifold point on the [`SymmetricPositiveSemidefiniteFixedRank`](@ref) `M`, i.e.
whether `p=q*q'` is a symmetric matrix of size `(n,n)` with values from the corresponding
[`AbstractNumbers`](@ref) `𝔽`.
The symmetry of `p` is not explicitly checked since by using `q` p is symmetric by construction.
The tolerance for the symmetry of `p` can and the rank of `q*q'` be set using `kwargs...`.
"""
function check_manifold_point(
    M::SymmetricPositiveSemidefiniteFixedRank{n,k,𝔽},
    q;
    kwargs...,
) where {n,k,𝔽}
    mpv =
        invoke(check_manifold_point, Tuple{supertype(typeof(M)),typeof(q)}, M, q; kwargs...)
    mpv === nothing || return mpv
    p = q * q'
    r = (VERSION >= v"1.1") ? rank(p * p'; kwargs...) : rank(p * p')
    if r < k
        return DomainError(
            r,
            "The point $(p) does not lie on $M, since its rank is less than $(k).",
        )
    end
    return nothing
end

"""
    check_tangent_vector(M::SymmetricPositiveSemidefiniteFixedRank{n,k,𝔽}, p, X; check_base_point = true, kwargs... )

Check whether `X` is a tangent vector to manifold point `p` on the
[`SymmetricPositiveSemidefiniteFixedRank`](@ref) `M`, i.e. `X` has to be a symmetric matrix of size `(n,n)`
and its values have to be from the correct [`AbstractNumbers`](@ref).
The optional parameter `check_base_point` indicates, whether to call
 [`check_manifold_point`](@ref)  for `p`.
The tolerance for the symmetry of `p` and `X` can be set using `kwargs...`.
"""
function check_tangent_vector(
    M::SymmetricPositiveSemidefiniteFixedRank{n,k,𝔽},
    q,
    Y;
    check_base_point = true,
    kwargs...,
) where {n,k,𝔽}
    if check_base_point
        mpe = check_manifold_point(M, q; kwargs...)
        mpe === nothing || return mpe
    end
    mpv = invoke(
        check_tangent_vector,
        Tuple{supertype(typeof(M)),typeof(q),typeof(Y)},
        M,
        q,
        Y;
        check_base_point = false, # already checked above
        kwargs...,
    )
    return mpv
end

function decorated_manifold(::SymmetricPositiveSemidefiniteFixedRank{N,K,𝔽}) where {N,K,𝔽}
    return Euclidean(N, K; field = 𝔽)
end

@doc raw"""
    distance(M::SymmetricPositiveSemidefiniteFixedRank, p, q)

Compute the distance between two points `p`, `q` on the
[`SymmetricPositiveSemidefiniteFixedRank`](@ref), which is the Frobenius norm of $Y$ which
minimizes $\lVert p - qY\rVert$ with respect to $Y$.
"""
distance(M::SymmetricPositiveSemidefiniteFixedRank, p, q) = norm(M, p, log(M, p, q))

@doc raw"""
    exp(M::SymmetricPositiveSemidefiniteFixedRank, q, Y)

Compute the exponential map on the [`SymmetricPositiveSemidefiniteFixedRank`](@ref),
which just reads

````math
    \exp_q Y = q+Y.
````
!!! note

    Since the manifold is represented in the embedding and is a quotient manifold,
    the exponential and logarithmic map are a bijection only with respect to the equivalence
    classes. Computing
    ````math
        q_2 = \exp_p(\log_pq)
    ````
    might yield a matrix $q_2\neq q$, but they represent the same point on the quotient
    manifold, i.e. $d_{\operatorname{SPS}_k(n)}(q_2,q) = 0$.
"""
exp(::SymmetricPositiveSemidefiniteFixedRank, ::Any, ::Any)

function exp!(::SymmetricPositiveSemidefiniteFixedRank, q2, q, Y)
    q2 .= q .+ Y
    return q2
end

@doc raw"""
    isapprox(M::SymmetricPositiveSemidefiniteFixedRank, p, q; kwargs...)

test, whether two points `p`, `q` are (approximately) nearly the same.
Since this is a quotient manifold in the embedding, the test is performed by checking
their distance, if they are not the same, i.e. that $d_{\mathcal M}(p,q) \approx 0$, where
the comparison is performed with the classical `isapprox`.
The `kwargs...` are passed on to this accordingly.
"""
function isapprox(M::SymmetricPositiveSemidefiniteFixedRank, p, q; kwargs...)
    return isapprox(norm(p - q), 0.0; kwargs...) ||
           isapprox(distance(M, p, q), 0.0; kwargs...)
end

@doc raw"""
    log(M::SymmetricPositiveSemidefiniteFixedRank, q, p)

Compute the logarithmic map on the [`SymmetricPositiveSemidefiniteFixedRank`](@ref) manifold
by minimizing $\lVert p - qY\rVert$ with respect to $Y$.

!!! note

    Since the manifold is represented in the embedding and is a quotient manifold,
    the exponential and logarithmic map are a bijection only with respect to the equivalence
    classes. Computing
    ````math
        q_2 = \exp_p(\log_pq)
    ````
    might yield a matrix $q_2\neq q$, but they represent the same point on the quotient
    manifold, i.e. $d_{\operatorname{SPS}_k(n)}(q_2,q) = 0$.
"""
log(::SymmetricPositiveSemidefiniteFixedRank, q, p)

function log!(::SymmetricPositiveSemidefiniteFixedRank, Z, q, p)
    s = svd(q' * p)
    return Z .= p * (s.V * s.U') - q
end

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

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

````math
\begin{aligned}
\dim \operatorname{SPS}_{k,ℝ}(n) &= kn - \frac{k(k-1)}{2},\\
\dim \operatorname{SPS}_{k,ℂ}(n) &= 2kn - k^2,
\end{aligned}
````

where the last $k^2$ is due to the zero imaginary part for Hermitian matrices diagonal
"""
manifold_dimension(::SymmetricPositiveSemidefiniteFixedRank)

@generated function manifold_dimension(
    ::SymmetricPositiveSemidefiniteFixedRank{N,K,ℝ},
) where {N,K}
    return K * N - div(K * (K - 1), 2)
end
@generated function manifold_dimension(
    ::SymmetricPositiveSemidefiniteFixedRank{N,K,ℂ},
) where {N,K}
    return 2 * K * N - K * K
end

function project!(::SymmetricPositiveSemidefiniteFixedRank, Z, q, Y)
    Z .= Y - q * lyap(q' * q, q' * Y - Y' * q)
    return Z
end

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

"""
     vector_transport_to(M::SymmetricPositiveSemidefiniteFixedRank, p, X, q)

 transport the tangent vector `X` at `p` to `q` by projecting it onto the tangent space
 at `q`.
 """
vector_transport_to(
    ::SymmetricPositiveSemidefiniteFixedRank,
    ::Any,
    ::Any,
    ::Any,
    ::ProjectionTransport,
)

function vector_transport_to!(
    M::SymmetricPositiveSemidefiniteFixedRank,
    Y,
    p,
    X,
    q,
    ::ProjectionTransport,
)
    project!(M, Y, q, X)
    return Y
end

@doc raw"""
     zero_tangent_vector(M::SymmetricPositiveSemidefiniteFixedRank, p)

returns the zero tangent vector in the tangent space of the symmetric positive
definite matrix `p` on the [`SymmetricPositiveSemidefiniteFixedRank`](@ref) manifold `M`.
"""
zero_tangent_vector(::SymmetricPositiveSemidefiniteFixedRank, ::Any...)

function zero_tangent_vector!(
    ::SymmetricPositiveSemidefiniteFixedRank{N,K},
    v,
    ::Any,
) where {N,K}
    fill!(v, 0)
    return v
end


text:
@cdef :PyFloat_FromDouble PyPtr (Cdouble,)
@cdef :PyFloat_AsDouble Cdouble (PyPtr,)

const PyFloat_Type__ref = Ref(PyNULL)
PyFloat_Type() = pyglobal(PyFloat_Type__ref, :PyFloat_Type)

PyFloat_Check(o) = Py_TypeCheck(o, PyFloat_Type())

PyFloat_CheckExact(o) = Py_TypeCheckExact(o, PyFloat_Type())

PyFloat_From(o::Union{Float16,Float32,Float64}) = PyFloat_FromDouble(o)

# "Floatable" means a 'float' or anything with a '__float__' method
PyFloatable_TryConvertRule_convert(o, ::Type{S}) where {S} = begin
    x = PyFloat_AsDouble(o)
    ism1(x) && PyErr_IsSet() && return -1
    putresult(convert(S, x))
end

PyFloatable_TryConvertRule_tryconvert(o, ::Type{S}) where {S} = begin
    x = PyFloat_AsDouble(o)
    ism1(x) && PyErr_IsSet() && return -1
    putresult(tryconvert(S, x))
end

ext: jl


text:
using Test

using GLPK

@testset "farmer" begin
    include("../examples/farmer.jl")
    status, objval, soln = DLP(m, with_optimizer(GLPK.Optimizer))
    @test status == :Optimal
    @test objval ≈ -108390
    @test soln ≈ [170, 80, 250]
end


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

JLLWrappers.@generate_main_file_header("SDPA")
JLLWrappers.@generate_main_file("SDPA", UUID("7fc90fd6-dbef-5a6a-93f8-169f2a2e705b"))
end  # module SDPA_jll


text:
@testset "Symmetric matrices" begin
    m = [
        1 2 1
        2 3 2
        1 2 1
    ]

    @test e(m) == m
    @test r90(m) == m
    @test r180(m) == m
    @test r270(m) == m
    @test tr(m) == m
    @test trr(m) == m
    @test fx(m) == m
    @test fy(m) == m
end

@testset "Non-symmetric matrices" begin
    m = [
        1 2 3
        4 5 6
        7 8 9
    ]

    @test e(m) == m
    @test r90(m) == [
        7 4 1
        8 5 2
        9 6 3
    ]
    @test r180(m) == [
        9 8 7
        6 5 4
        3 2 1
    ]
    @test r180(m) == [
        9 8 7
        6 5 4
        3 2 1
    ]
    @test r270(m) == [
        3 6 9
        2 5 8
        1 4 7
    ]
    @test tr(m) == [
        1 4 7
        2 5 8
        3 6 9
    ]
    @test trr(m) == [
        9 6 3
        8 5 2
        7 4 1
    ]
    @test fx(m) == [
        7 8 9
        4 5 6
        1 2 3
    ]
    @test fy(m) == [
        3 2 1
        6 5 4
        9 8 7
    ]
end

@testset "Group operation" begin
    # the '∘' function acts as function composition and group operation at the same time
    # it can be used prevent unneeded computations, e.g., (r180 ∘ r180)(m) does not actually
    # rotate m twice by 180 degrees but directly returns e(m)
    m = rand(10, 10)

    # validate that taking these shortcuts returns the expected results
    @test (e ∘ e)(m) == e(e(m))
    @test (e ∘ fx)(m) == e(fx(m))
    @test (e ∘ fy)(m) == e(fy(m))
    @test (e ∘ tr)(m) == e(tr(m))
    @test (e ∘ trr)(m) == e(trr(m))
    @test (e ∘ r90)(m) == e(r90(m))
    @test (e ∘ r180)(m) == e(r180(m))
    @test (e ∘ r270)(m) == e(r270(m))

    @test (fx ∘ e)(m) == fx(e(m))
    @test (fx ∘ fx)(m) == fx(fx(m))
    @test (fx ∘ fy)(m) == fx(fy(m))
    @test (fx ∘ tr)(m) == fx(tr(m))
    @test (fx ∘ trr)(m) == fx(trr(m))
    @test (fx ∘ r90)(m) == fx(r90(m))
    @test (fx ∘ r180)(m) == fx(r180(m))
    @test (fx ∘ r270)(m) == fx(r270(m))

    @test (fy ∘ e)(m) == fy(e(m))
    @test (fy ∘ fx)(m) == fy(fx(m))
    @test (fy ∘ fy)(m) == fy(fy(m))
    @test (fy ∘ tr)(m) == fy(tr(m))
    @test (fy ∘ trr)(m) == fy(trr(m))
    @test (fy ∘ r90)(m) == fy(r90(m))
    @test (fy ∘ r180)(m) == fy(r180(m))
    @test (fy ∘ r270)(m) == fy(r270(m))

    @test (tr ∘ e)(m) == tr(e(m))
    @test (tr ∘ fx)(m) == tr(fx(m))
    @test (tr ∘ fy)(m) == tr(fy(m))
    @test (tr ∘ tr)(m) == tr(tr(m))
    @test (tr ∘ trr)(m) == tr(trr(m))
    @test (tr ∘ r90)(m) == tr(r90(m))
    @test (tr ∘ r180)(m) == tr(r180(m))
    @test (tr ∘ r270)(m) == tr(r270(m))

    @test (trr ∘ e)(m) == trr(e(m))
    @test (trr ∘ fx)(m) == trr(fx(m))
    @test (trr ∘ fy)(m) == trr(fy(m))
    @test (trr ∘ tr)(m) == trr(tr(m))
    @test (trr ∘ trr)(m) == trr(trr(m))
    @test (trr ∘ r90)(m) == trr(r90(m))
    @test (trr ∘ r180)(m) == trr(r180(m))
    @test (trr ∘ r270)(m) == trr(r270(m))

    @test (r90 ∘ e)(m) == r90(e(m))
    @test (r90 ∘ fx)(m) == r90(fx(m))
    @test (r90 ∘ fy)(m) == r90(fy(m))
    @test (r90 ∘ tr)(m) == r90(tr(m))
    @test (r90 ∘ trr)(m) == r90(trr(m))
    @test (r90 ∘ r90)(m) == r90(r90(m))
    @test (r90 ∘ r180)(m) == r90(r180(m))
    @test (r90 ∘ r270)(m) == r90(r270(m))

    @test (r180 ∘ e)(m) == r180(e(m))
    @test (r180 ∘ fx)(m) == r180(fx(m))
    @test (r180 ∘ fy)(m) == r180(fy(m))
    @test (r180 ∘ tr)(m) == r180(tr(m))
    @test (r180 ∘ trr)(m) == r180(trr(m))
    @test (r180 ∘ r90)(m) == r180(r90(m))
    @test (r180 ∘ r180)(m) == r180(r180(m))
    @test (r180 ∘ r270)(m) == r180(r270(m))

    @test (r270 ∘ e)(m) == r270(e(m))
    @test (r270 ∘ fx)(m) == r270(fx(m))
    @test (r270 ∘ fy)(m) == r270(fy(m))
    @test (r270 ∘ tr)(m) == r270(tr(m))
    @test (r270 ∘ trr)(m) == r270(trr(m))
    @test (r270 ∘ r90)(m) == r270(r90(m))
    @test (r270 ∘ r180)(m) == r270(r180(m))
    @test (r270 ∘ r270)(m) == r270(r270(m))
end


text:
using StochasticDiffEq, DiffEqDevTools, LinearAlgebra, Random, Test
Random.seed!(100)
dts = (1/2) .^ (10:-1:2) #14->7 good plot

const σ_const = 0.87
const μ = 1.01

u0 = rand(2)
A = [-2.0 1.0;1.0 -2.0]
B = Diagonal([σ_const for i in 1:2])

function f_commute(du,u,p,t)
  mul!(du,A,u)
  du .+= 1.01u
end

function f_commute_oop(u,p,t)
  du = A*u
  du += 1.01u
end

function f_commute_analytic(u0,p,t,W)
 tmp = (A+1.01I-2*(B^2))*t + B*sum(W)
 exp(tmp)*u0
end

function σ(du,u,p,t)
  du[1,1] = σ_const*u[1]
  du[2,1] = σ_const*u[2]
  du[1,2] = σ_const*u[1]
  du[2,2] = σ_const*u[2]
  du[1,3] = σ_const*u[1]
  du[2,3] = σ_const*u[2]
  du[1,4] = σ_const*u[1]
  du[2,4] = σ_const*u[2]
end

function σ_oop(u,p,t)
  σ_const*[
           u[1]  u[1]  u[1]  u[1]
           u[2]  u[2]  u[2]  u[2]
           ]
end

ff_commute = SDEFunction(f_commute,σ,analytic=f_commute_analytic)

prob = SDEProblem(ff_commute,σ,u0,(0.0,1.0),noise_rate_prototype=rand(2,4))

sol = solve(prob,RKMilCommute(),dt=1/2^(8))
sol = solve(prob,RKMilGeneral(ii_approx=IICommutative()),dt=1/2^(8))
sol = solve(prob,RKMilGeneral(p=2),dt=1/2^(10))
sol = solve(prob,EM(),dt=1/2^(10))

dts = (1/2) .^ (10:-1:3) #14->7 good plot
sim2 = test_convergence(dts,prob,EM(),trajectories=Int(1e2))
@test abs(sim2.𝒪est[:final] - 0.5) < 0.2

sim2 = test_convergence(dts,prob,RKMilCommute(),trajectories=Int(2e2))
@test abs(sim2.𝒪est[:final] - 1.0) < 0.2

sim2 = test_convergence(dts,prob,RKMilGeneral(ii_approx=IICommutative()),trajectories=Int(2e2))
@test abs(sim2.𝒪est[:final] - 1.0) < 0.2

sim2 = test_convergence(dts,prob,RKMilGeneral(p=2),trajectories=Int(2e2))
@test abs(sim2.𝒪est[:final] - 1.0) < 0.2

ff_commute_oop = SDEFunction(f_commute_oop,σ_oop,analytic=f_commute_analytic)

proboop = SDEProblem(ff_commute_oop,σ_oop,u0,(0.0,1.0),noise_rate_prototype=rand(2,4))

sol = solve(proboop,RKMilCommute(),dt=1/2^(8))
sol = solve(proboop,RKMilGeneral(ii_approx=IICommutative()),dt=1/2^(8))
sol = solve(proboop,RKMilGeneral(p=2),dt=1/2^(10))
sol = solve(proboop,EM(),dt=1/2^(10))

dts = (1/2) .^ (10:-1:3) #14->7 good plot
sim2 = test_convergence(dts,proboop,EM(),trajectories=Int(1e2))
@test abs(sim2.𝒪est[:final] - 0.5) < 0.2

sim2 = test_convergence(dts,proboop,RKMilCommute(),trajectories=Int(2e2))
@test abs(sim2.𝒪est[:final] - 1.0) < 0.2

sim2 = test_convergence(dts,proboop,RKMilGeneral(ii_approx=IICommutative()),trajectories=Int(2e2))
@test abs(sim2.𝒪est[:final] - 1.0) < 0.2

sim2 = test_convergence(dts,proboop,RKMilGeneral(p=2),trajectories=Int(2e2))
@test abs(sim2.𝒪est[:final] - 1.0) < 0.2


text:
include("../src/Rhapsodie.jl")
include("test_data_simulation.jl")
#-----------------------------------------------
# Loading the dataset parameters
par=Rhapsodie.readdlm("test_results/Parameters.txt")

DSIZE=Int64(par[1]);
NTOT=Int64(par[2]);
Nframe=Int64(par[3]);
Nrot=Int64(par[4]);
Nangle=NTOT÷(Nframe*4)
Center=par[5:6];
DerotAng=Vector{Float64}();
Epsilon=Vector{Tuple{Array{Float64,1},Array{Float64,1}}}();

for iter=1:NTOT
    ind=div(iter-1, NTOT/4)+1
    push!(Epsilon,([0. ,0. ],par[end-1:end]));
    push!(DerotAng, deg2rad(-par[7]));
end

psf_center=Rhapsodie.readdlm("../data/PSF_centers_Airy.txt");

Rhapsodie.load_parameters((DSIZE, 2*DSIZE, NTOT), Nframe, Nrot, Nangle, Center, (psf_center[1:2], psf_center[3:4]), Epsilon, DerotAng)

for tau in [0.01, 0.03, 0.07, 0.1, 0.15, 0.25, 0.5]
#----------------------------------------------
# Loading the pre-processed data (i.e. bad pixels have been interpolated
    Rhapsodie.Load_Data("test_results/DATA_$tau-$DSIZE.fits",
                         "test_results/WEIGHT_$tau-$DSIZE.fits")
	Sdim=length(Rhapsodie.dataset)

    DATA=zeros(Rhapsodie.get_par().cols[1], Rhapsodie.get_par().cols[2], Sdim,2);
    WEIGHT=zeros(Rhapsodie.get_par().cols[1], Rhapsodie.get_par().cols[2], Sdim,2);

    ker= Rhapsodie.CatmullRomSpline(Float64, Rhapsodie.Flat)
    input_size=(Rhapsodie.get_par().rows[1], Rhapsodie.get_par().rows[2]÷2);
    output_size= Rhapsodie.get_par().cols[1:2];

# Pre processing
    for i=1:Sdim
        T1=Rhapsodie.TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, inv(Rhapsodie.Trans_Table[i][1]))
        T2=Rhapsodie.TwoDimensionalTransformInterpolator(output_size, input_size, ker, ker, inv(Rhapsodie.Trans_Table[i][2]))

        I1=T1*Rhapsodie.dataset[i].data[:,1:end÷2]
        I2=T2*Rhapsodie.dataset[i].data[:,end÷2+1:end]
        W1=T1*Rhapsodie.dataset[i].weights[:,1:end÷2]
        W2=T2*Rhapsodie.dataset[i].weights[:,end÷2+1:end]

        DATA[:,:,i,1]=I1;
        DATA[:,:,i,2]=I2;
        WEIGHT[:,:,i,1]=W1;
        WEIGHT[:,:,i,2]=W2;
    end

#-----------------------------------------------------
# Reconstruction
# Non-Linear Separable inverse method
    #MI=Rhapsodie.NonLinear_Method(DATA, WEIGHT);
    #Rhapsodie.write(MI, "test_results/Results_Separable_NonLinear_$tau-$DSIZE.fits")

# Linear Separable inverse method
    ML=Rhapsodie.Linear_Method(DATA, WEIGHT);
    Rhapsodie.write(ML, "test_results/Results_Separable_Linear_$tau-$DSIZE.fits")

# Double Ratio
    DR=Rhapsodie.Double_Ratio(DATA);
    Rhapsodie.write(DR, "test_results/Results_Separable_DoubleRatio_$tau-$DSIZE.fits")

# Double Difference
    DD=Rhapsodie.Double_Difference(DATA);
    Rhapsodie.write(DD, "test_results/Results_Separable_DoubleDifference_$tau-$DSIZE.fits")

    empty!(Rhapsodie.dataset);
end


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

using Markdown
using InteractiveUtils

# ╔═╡ 6abe6400-476f-11ec-0423-afbc13ad4a16
md"""
# Genie Plugins

Genie plugins are special Julia packages which extend Genie apps with powerful  functionality by providing specific integration points. A Genie plugin is made of two parts:

1. the Julia package exposing the core functionality of the plugin, and
2. a files payload (controllers, modules, views, database migrations, initializers, etc) which are copied into the client app upon plugin installation.

## Using Genie Plugins

The plugins are created by third party Genie/Julia developers. Take this simple demo plugin as an example: <https://github.com/GenieFramework/HelloPlugin.jl>

In order to add the plugin to an existing Genie app you need to:

Add the `HelloPlugin` package to your Genie app, just like any other Julia Pkg dependency:

"""

# ╔═╡ 9b87679e-66c0-4551-b166-44e56dd75745
md"""

```julia
julia> using Pkg

julia> Pkg.add(url="https://github.com/GenieFramework/HelloPlugin.jl")
```
"""

# ╔═╡ 2c393262-ee80-4a99-bf56-ac98b021be7f
md"""
Bring the package into scope:

"""

# ╔═╡ f8b6cd4b-891c-493b-8ef2-e9a108d20633
md"""

```julia
julia> using HelloPlugin
```
"""

# ╔═╡ 1abffbef-0109-417a-8788-0894fc64f40a
md"""
Install the plugin (this is a one time operation, when the package is added):

```julia

julia> HelloPlugin.install(@__DIR__)
```
"""

# ╔═╡ 257ae400-476b-4c3c-b3f9-766e583900ea
md"""
### Running the Plugin

The installation will add a new `hello` resource in `app/resources/hello/` in the form of `HelloController.jl` and `views/greet.jl.html`. Also, in your Genie app's `plugins/` folder you fill find a new file, `helloplugin.jl` (which is the plugins' initializer and is automatically loaded by Genie early in the bootstrap process).

The `helloplugin.jl` initializer is defining a new route `route("/hello", HelloController.greet)`. If you restart your Genie app and navigate to `/hello` you will get the plugin's greeting.

## Walkthrough

Create a new Genie app:
"""

# ╔═╡ 351a0df7-0566-4824-89d0-9ab3af54fc56
md"""

```julia

julia> using Genie

julia> Genie.newapp("Greetings", autostart = false)
```
"""

# ╔═╡ de53d6e1-f3eb-4b9c-8e4e-76a9fd7d5e4e
md"""

Add the plugin as a dependency:

```julia

julia> using Pkg

julia> Pkg.add(url="https://github.com/GenieFramework/HelloPlugin.jl")
```
"""

# ╔═╡ 8ad0a588-4569-4e1c-8bb3-46c8a19881a1
md"""
Bring the plugin into scope and run the installer (the installer should be run only once, upon adding the plugin package)

```julia
julia> using HelloPlugin
julia> HelloPlugin.install(@__DIR__)
```

The installation might show a series of logging messages informing about failure to copy some files or create folders. Normally it's nothing to worry about: these are due to the fact that some of the files and folders already exist in the app so they are not overwritten by the installer.

Restart the app to load the plugin:

```julia
julia> exit()
```

```bash
$ cd Greetings/
$ bin/repl
```

Start the server:

```julia
julia> Genie.startup()
```

Navigate to `http://localhost:8000/hello` to get the greeting from the plugin.

---

## Developing Genie Plugins

Genie provides an efficient scaffold for bootstraping a new plugin package. All you need to do is run this code to create your plugin project:

"""

# ╔═╡ 6678ed9c-3dc3-4615-8b8c-4da9a4c6aee2
md"""

```julia
julia> using genie
```
"""

# ╔═╡ d9a46470-9aac-48d5-9c2b-800b1b5b23b7
md"""

```julia
julia> Genie.Plugins.scaffold("GenieHelloPlugin") # use the actual name of your plugin
Generating project file
Generating project GenieHelloPlugin:
    GenieHelloPlugin/Project.toml
    GenieHelloPlugin/src/GenieHelloPlugin.jl
Scaffolding file structure

Adding dependencies
  Updating registry at `~/.julia/registries/General`
  Updating git-repo `https://github.com/JuliaRegistries/General.git`
  Updating git-repo `https://github.com/genieframework/Genie.jl`
 Resolving package versions...
  Updating `~/GenieHelloPlugin/Project.toml`
  [c43c736e] + Genie v0.9.4 #master (https://github.com/genieframework/Genie.jl)
  Updating `~/GenieHelloPlugin/Manifest.toml`

Initialized empty Git repository in /Users/adrian/GenieHelloPlugin/.git/
[master (root-commit) 30533f9] initial commit
 11 files changed, 261 insertions(+)

Congratulations, your plugin is ready!
You can use this default installation function in your plugin's module:
  function install(dest::String; force = false)
    src = abspath(normpath(joinpath(@__DIR__, "..", Genie.Plugins.FILES_FOLDER)))

    for f in readdir(src)
      isdir(f) || continue
      Genie.Plugins.install(joinpath(src, f), dest, force = force)
    end
  end
```

The scaffold command will create the file structure of your plugin, including the Julia project, the `.git` repo, and the file structure for integrating with Genie apps:

```
.
├── Manifest.toml
├── Project.toml
├── files
│   ├── app
│   │   ├── assets
│   │   │   ├── css
│   │   │   ├── fonts
│   │   │   └── js
│   │   ├── helpers
│   │   ├── layouts
│   │   └── resources
│   ├── db
│   │   ├── migrations
│   │   └── seeds
│   ├── lib
│   ├── plugins
│   │   └── geniehelloplugin.jl
│   └── task
└── src
    └── GenieHelloPlugin.jl
```

The core of the functionality shoud go into the `src/GenieHelloPlugin.jl` module. While everything placed within the `files/` folder should be copied into the corresponding folders of the Genie apps installing the plugin. You can add resources, controllers, models, database migrations, views, assets and any other files inside the `files/` folder to be copied.

The scaffolding will also create a `plugins/geniehelloplugin.jl` file - this is the initializer of the plugin and is meant to bootstrap the functionality of the plugin. Here you can load dependencies, define routes, set up configuration, etc.

Because any Genie plugin is a Julia `Pkg` project, you can add any other Julia packages as dependencies.

### The Installation Function

The main module file, present in `src/GenieHelloPlugin.jl` should also expose an `install(path::String)` function, responsible for copying the files of your plugin into the user Genie app. The `path` param is the root of the Genie app where the installation will be performed.

As copying the plugin's files is a standard but tedious operation, Genie provides some helpers to get you started. The `Genie.Plugins` module provides an `install(path::String, dest::String; force = false)` which can be used for copying the plugin's files to their destination in the app.

The scaffolding function will also recommend a default `install(path::String)` that you can use in your module:

```julia
function install(dest::String; force = false)
  src = abspath(normpath(joinpath(@__DIR__, "..", Genie.Plugins.FILES_FOLDER)))

  for f in readdir(src)
    isdir(f) || continue
    Genie.Plugins.install(joinpath(src, f), dest, force = force)
  end
end
```

You can use it as a starting point - and add any other specific extra logic to it.

"""

# ╔═╡ Cell order:
# ╟─6abe6400-476f-11ec-0423-afbc13ad4a16
# ╟─9b87679e-66c0-4551-b166-44e56dd75745
# ╟─2c393262-ee80-4a99-bf56-ac98b021be7f
# ╟─f8b6cd4b-891c-493b-8ef2-e9a108d20633
# ╟─1abffbef-0109-417a-8788-0894fc64f40a
# ╟─257ae400-476b-4c3c-b3f9-766e583900ea
# ╟─351a0df7-0566-4824-89d0-9ab3af54fc56
# ╟─de53d6e1-f3eb-4b9c-8e4e-76a9fd7d5e4e
# ╟─8ad0a588-4569-4e1c-8bb3-46c8a19881a1
# ╟─6678ed9c-3dc3-4615-8b8c-4da9a4c6aee2
# ╟─d9a46470-9aac-48d5-9c2b-800b1b5b23b7


text:
##################################################################################
# This file is part of StateSpaceEcon.jl
# BSD 3-Clause License
# Copyright (c) 2020, Bank of Canada
# All rights reserved.
##################################################################################

"""
    sim_nr!(x, sd, maxiter, tol, verbose [, sparse_solver])

Solve the simulation problem.
  * `x` - the array of data for the simulation. All initial, final and exogenous
    conditions are already in place.
  * `sd::AbstractSolverData` - the solver data constructed for the simulation
    problem.
  * `maxiter` - maximum number of iterations.
  * `tol` - desired accuracy.
  * `verbose` - whether or not to print progress information.
  * `sparse_solver` (optional) - a function called to solve the linear system A
    x = b for x. Defaults to A\\b

"""
function sim_nr!(x::AbstractArray{Float64}, sd::StackedTimeSolverData,
                maxiter::Int64, tol::Float64, verbose::Bool,
                sparse_solver::Function=(A, b) -> A \ b)
    for it = 1:maxiter
        @timer Fx, Jx = global_RJ(x, x, sd)
        @timer nFx = norm(Fx, Inf)
        if nFx < tol
            if verbose
                @info "$it, || Fx || = $(nFx)"
            end
            break
        end
        @timer Δx = sparse_solver(Jx, Fx)
        @timer nΔx = norm(vec(Δx), Inf)
        assign_update_step!(x, -1.0, Δx, sd)
        if verbose
            @info "$it, || Fx || = $(nFx), || Δx || = $(nΔx)"
        end
        if nΔx < tol
            break
        end
    end
    return nothing
end

export simulate
"""
    simulate(model, plan, data; <options>)

Run a simulation for the given model, simulation plan and exogenous data.

### Arguments
  * `model` - the [`Model`](@ref) instance to simulate.
  * `plan` - the [`Plan`](@ref) for the simulation.
  * `data` - a 2D `Array` containing the exogenous data. This includes the
    initial and final conditions.

### Options as keyword arguments
  * `fctype::`[`FCType`](@ref) - set the desired final condition type for the
    simulation. The default value is `fcgiven`. Other possible values
    include `fclevel` and `fcslope`.
  * `initial_guess::AbstractMatrix{Float64}` - a 2D `Array` containing the
    initial guess for the solution. This is used to start the Newton-Raphson
    algorithm. The default value is an empty array (`zeros(0,0)`), in which case
    we use the exogenous data for the initial condition. You can use the steady
    state solution using [`steadystatearray`](@ref).
  * `linearize::Bool` - set to `true` to instruct the solver to use the
    liearized model. If the model is already linearized, this option has the
    effect that the model gets linearized about the current steady stat and with
    the value of `deviation` given here. Otherwise the model is linearized about
    the steady state. After the simulation is computed, the model is restored to
    its original state. Default value is `false`.
  * `deviation::Bool` - set to `true` if the `data` is in deviations from the
    steady state. This is only relevant if the `linearize` option is set to
    `true`. Default value is `false`.
  * `anticipate::Bool` - set to `false` to instruct the solver that all shocks
    are unanticilated by the agents. Default value is `true`.
  * `verbose::Bool` - control whether or not to print progress information.
    Default value is taken from `model.options`.
  * `tol::Float64` - set the desired accuracy. Default value is taken from
    `model.options`.
  * `maxiter::Int` - algorithm fails if the desired accuracy is not reached
    within this maximum number of iterations. Default value is taken from
    `model.options`.

### See also:

### Examples

"""
function simulate(m::Model, p::Plan, exog_data::AbstractArray{Float64,2};
                    initial_guess::AbstractArray{Float64,2}=zeros(0, 0),
                    linearize::Bool=false,
                    deviation::Bool=false,
                    anticipate::Bool=true,
                    verbose::Bool=m.options.verbose,
                    tol::Float64=m.options.tol,
                    maxiter::Int64=m.options.maxiter,
                    fctype=getoption(m, :fctype, fcgiven),
                    expectation_horizon::Union{Nothing,Int64}=nothing,
                    sparse_solver::Function=(A, b) -> A \ b
    )
    NT = length(p.range)
    nauxs = length(m.auxvars)
    if size(exog_data) != (NT, length(m.varshks))
        error("Incorrect dimensions of exog_data. Expected $((NT, length(m.varshks))), got $(size(exog_data)).")
    end
    if !isempty(initial_guess) && size(initial_guess) != (NT, length(m.varshks))
        error("Incorrect dimensions of initial_guess. Expected $((NT, length(m.varshks))), got $(size(exog_data)).")
    end

    if deviation
        exog_data = copy(exog_data)
        local logvars = islog.(m.varshks) .| lsneglog.(m.varshks)
        local ss_data = steadystatearray(m, p)
        exog_data[:, logvars] .*= ss_data[:, logvars]
        exog_data[:, .! logvars] .+= ss_data[:, .! logvars]
    end

    exog_data = @timer ModelBaseEcon.update_auxvars(transform(exog_data, m), m)

    if !isempty(initial_guess)
        x = @timer ModelBaseEcon.update_auxvars(transform(initial_guess, m), m)
    else
        x = copy(exog_data)
    end

    if linearize
        throw(ErrorException("Linearization is disabled."))
        org_med = m.evaldata
        linearize!(m; deviation=deviation)
    end

    if anticipate
        @timer gdata = StackedTimeSolverData(m, p, fctype)
        @timer assign_exog_data!(x, exog_data, gdata)
        if verbose
            @info "Simulating $(p.range[1 + m.maxlag:NT - m.maxlead])" # anticipate gdata.FC
        end
        @timer sim_nr!(x, gdata, maxiter, tol, verbose, sparse_solver)
    else # unanticipated shocks
        init = 1:m.maxlag
        term = NT .+ (1 - m.maxlead:0)
        sim = 1 + m.maxlag:NT - m.maxlead
        nvars = length(m.variables)
        nshks = length(m.shocks)
        nauxs = length(m.auxvars)
        shkinds = nvars .+ (1:nshks)
        auxinds = nvars .+ nshks .+ (1:nauxs)
        varshkinds = 1:(nvars + nshks)
        allvarinds = 1:(nvars + nshks + nauxs)
        x[init, allvarinds] .= exog_data[init, allvarinds]
        # x[term, allvarinds] .= exog_data[term, allvarinds]
        x[sim, shkinds] .= 0.0
        t0 = first(sim)
        T = last(sim)
        if expectation_horizon === nothing
            # the original code, where we always simulate until the end with the true final conditions
            for t in sim
                exog_inds = p[t, Val(:inds)]
                psim = Plan(m, t:T)
                if (t != t0) && (psim[t0, Val(:inds)] == exog_inds) && (maximum(abs, exog_data[t,shkinds]) < tol)
                    continue
                end
                setexog!(psim, t0, exog_inds)
                @timer gdata = StackedTimeSolverData(m, psim, fctype)
                x[t,exog_inds] = exog_data[t,exog_inds]
                # @timer assign_exog_data!(x[psim.range,:], exog_data[psim.range,:], gdata)
                sim_range = UnitRange{Int}(psim.range)
                xx = view(x, sim_range, :)
                assign_final_condition!(xx, exog_data[sim_range,:], gdata)
                if verbose
                    @info "Simulating $(p.range[t:T])" # anticipate expectation_horizon gdata.FC
                end
                @timer sim_nr!(xx, gdata, maxiter, tol, verbose, sparse_solver)
                # if verbose
                #     @info "Result at $t: " x[[t], :]
                # end
            end
        else
            # the new code, where the first and last simulations use the true
            # simulation range and final condition, while the intermediate
            # simulations use expectation_horizon steps with fcslope
            if expectation_horizon == 0
                expectation_horizon = length(sim)
            elseif expectation_horizon < 10 * m.maxlead
                @warn "Expectation horizon may be too short for this model. Consider setting it to at least $(10 * m.maxlead)."
            end
            x = [x; zeros(expectation_horizon, size(x, 2))]
            ninit = length(init)
            nterm = length(term)
            # first simulation
            let t = t0
                # first run is with the full range, the true fctype,
                # and only the first period is imposed
                exog_inds = p[t, Val(:inds)]
                psim = Plan(m, t:T)
                setexog!(psim, t0, exog_inds)
                @timer sdata = StackedTimeSolverData(m, psim, fctype)
                x[t,exog_inds] = exog_data[t,exog_inds]
                sim_range = UnitRange{Int}(psim.range)
                xx = view(x, sim_range, :)
                assign_final_condition!(xx, exog_data[sim_range,:], sdata)
                if verbose
                    @info "Simulating $(p.range[t:T])" # anticipate expectation_horizon sdata.FC
                end
                sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver)
                # if verbose
                #     @info "Result at $t: " x[[t], :]
                # end
            end
            # intermediate simulations
            last_t::Int64 = t0
            psim = Plan(m, 0:expectation_horizon - 1)
            sdata = StackedTimeSolverData(m, psim, fcnatural)
            for t in sim[2:end]
                exog_inds = p[t, Val(:inds)]
                # we need to run a simulation if a variable is exogenous, or if a shock value is not zero
                # these intermediate simulations are always with fcnatural,
                #       have length equal to expectation_horizon and
                #       only the first period is imposed
                if (exog_inds == shkinds) && (maximum(abs, exog_data[t, shkinds]) <= tol)
                    continue
                end
                setexog!(psim, t0, exog_inds)
                update_plan!(sdata, m, psim)
                # note that the range always goes from 0 to expectation_horizon-1,
                # so we need to add t in order to get the correct set of rows of x
                sim_range = t .+ UnitRange(psim.range)
                xx = view(x, sim_range, :)
                # The initial conditions are already set
                # The exogenous values are already set as well, except for the first period
                # In other words, we only need to impose the first period here
                xx[t0, exog_inds] = exog_data[t, exog_inds]
                # Update the final conditions
                @timer assign_final_condition!(xx, similar(xx), sdata)
                if verbose
                    @info "Simulating $(p.range[t] .+ (0:expectation_horizon - 1))" # anticipate expectation_horizon sdata.FC
                end
                @timer sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver)
                # if verbose
                #     @info "Result at $t: " x[[t], :]
                # end
                last_t = t  # keep track of last simulation time
            end
            # last simulation
            if last_t > t0
                # do we need to re-run the last simulation?
                # if it didn't reach T, then yes
                # if the final condition is not fcnatural, then yes
                if (last_t + expectation_horizon < T) || (fctype != fcnatural)
                    psim = Plan(m, min(last_t + 1, T):T)
                    # there are no unanticipated shocks in this simulation
                    sdata = StackedTimeSolverData(m, psim, fctype)
                    # the initial conditions and the exogenous data are already in x
                    # we only need the final conditions
                    sim_range = UnitRange(psim.range)
                    xx = view(x, sim_range, :)
                    assign_final_condition!(xx, exog_data[sim_range, :], sdata)
                    if verbose
                        @info "Simulating $(p.range[last_t + 1:T])" # anticipate expectation_horizon sdata.FC
                    end
                    @timer sim_nr!(xx, sdata, maxiter, tol, verbose, sparse_solver)
                    # if verbose
                    #     @info "Result at $(last_t + 1): " x[[last_t + 1], :]
                    # end
                end
            end
            x = x[1:end - expectation_horizon, :]
        end
    end
    if linearize
        m.evaldata = org_med
    end

    x = x[:,1:end - nauxs]
    x .= inverse_transform(x, m)
    if deviation
        x[:, logvars] ./= ss_data[:, logvars]
        x[:, .! logvars] .-= ss_data[:, .! logvars]
    end

    return x
end

# The versions of simulate with Dict

# Simulate command, IRIS style with a range
@inline simulate(m::Model, data, rng, p::Plan=Plan(m, rng); kwargs...) = simulate(m, data, p[rng, m]; kwargs...)

# Simulate command, IRIS style but without a range
function simulate(m::Model, D1::Dict{<:AbstractString,<:Any}, plan::Plan;
    overlay::Bool=false, kwargs...)::Dict{String,<:Any}
    # Convert dictionary to Array{Float64,2}
    data01 = dict2array(D1, m.varshks, range=plan.range)
    ig = get(kwargs, :initial_guess, zeros(0, 0))
    if ig isa Dict
        ig = dict2array(ig, m.varshks, range=plan.range)
    end
    # Call native simulate command with Array{Float64,2}
    data02 = simulate(m, plan, data01; kwargs..., initial_guess=ig)
    # Reconvert Array{Float64,2} to Dict{String,Any}
    D2 = array2dict(data02, m.varshks, plan.range[1])
    # Overlay D1 and D2
    if overlay
        D3 = dictoverlay(D1, D2)
    else
        D3 = D2
    end
    return D3
end

import ..SimData

function simulate(m::Model, D1::SimData, plan::Plan; overlay::Bool=false, kwargs...)::typeof(D1)
    ret = copy(D1)
    ig = get(kwargs, :initial_guess, zeros(0, 0))
    if ig isa SimData
        ig = ig[plan.range, m.varshks]
    end
    sim = simulate(m, plan, D1[plan.range, m.varshks]; kwargs..., initial_guess=ig)
    # overlay the sim data onto ret
    ret[plan.range, m.varshks] = sim
    return ret
end


text:
#=
# Contains info about all nodes on the domain
# Unlike MeshData struct, this accounts for interior nodes and corresponds to nodal DOFs.
=#

struct Grid
    allnodes::Array{Float64}        # All node coordinates size = (dim, nnodes)
    # boundaries
    bdry::Array{Array{Int,1},1}     # Indices of boundary nodes for each BID (bdry[bid][nodes])*note:array of arrays
    bdryface::Array{Array{Int,1},1} # Indices of faces touching each BID (bdryface[bid][faces])*note:array of arrays
    bdrynorm::Array{Array{Float64,2},1} # Normal vector for boundary nodes for each BID (bdrynorm[bid][dim, nodes])*note:array of arrays
    bids::Array{Int,1}              # BID corresponding to rows of bdrynodes
    # elements
    loc2glb::Array{Int,2}           # local to global map for each element's nodes (size is (Np, nel))
    glbvertex::Array{Int,2}         # global indices of each elements' vertices (size is (Nvertex, nel))
    # faces (For CG, G=1. For DG, G=2)
    face2glb::Array{Int,3}          # local to global map for faces (size is (Nfp, G, Nfaces))
    element2face::Array{Int,2}      # face indices for each element (size is (Nfaces, nel))
    face2element::Array{Int,2}      # elements on both sides of a face, 0=boundary (size is (2, Nfaces))
    facenormals::Array{Float64,2}   # normal vector for each face
    faceRefelInd::Array{Int,2}      # Index for face within the refel for each side

    facebid::Array{Int,1}           # BID of each face (0=interior face)
end

etypetonv = [2, 3, 4, 4, 8, 6, 5, 2, 3, 4, 4, 8, 6, 5, 1, 4, 8, 6, 5]; # number of vertices for each type
etypetodim= [1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3, 3, 3, 1, 2, 2, 3, 3]; # dimension of each type
etypetonf = [2, 3, 4, 4, 6, 5, 5, 2, 3, 4, 4, 6, 5, 5, 1, 4, 6, 5, 5]; # number of faces for element types
etypetoftype=[0,1, 1, 2, 3, 3, 3, 0, 1, 1, 2, 3, 3, 3, 0, 0, 0, 0, 0]; # type of faces for this element type

# Build a grid from a mesh
function grid_from_mesh(mesh)
    log_entry("Building full grid from mesh data", 2);
    t_grid_from_mesh = Base.Libc.time();
    dim = config.dimension;
    ord = config.basis_order_min;
    nfaces = etypetonf[mesh.etypes[1]];
    #totalfaces = nfaces*mesh.nel;
    totalfaces = size(mesh.normals,2);
    nel = mesh.nel;
    if dim == 1
        facenvtx = 1
    else
        facenvtx = etypetonv[etypetoftype[mesh.etypes[1]]]; # Assumes one element type
    end
    nvtx = etypetonv[mesh.etypes[1]]; # Assumes one element type

    log_entry("Building reference element: "*string(dim)*"D, order="*string(ord)*", nfaces="*string(nfaces), 3);
    refel = build_refel(dim, ord, nfaces, config.elemental_nodes);

    if config.solver_type == DG # || config.solver_type == FV ????
        Gness = 2;
    else
        Gness = 1;
    end

    tol = 1e-12; # tolerance for is_same_node

    Np = refel.Np;                      # number of nodes per element
    bdry = [];                          # index(in x) of boundary nodes for each BID
    bdryfc = [];                        # index of faces touching each BID
    bdrynorm = [];                      # normal at boundary nodes
    bids = collectBIDs(mesh);           # BID list
    nbids = length(bids);
    for i=1:nbids
        push!(bdry, zeros(Int, 0));
        push!(bdryfc, zeros(Int,0));
        push!(bdrynorm, zeros(config.dimension,0));
    end
    loc2glb = zeros(Int, Np, nel)       # local to global index map for each element's nodes
    glbvertex = zeros(Int, nvtx, nel);     # local to global for vertices
    f2glb = zeros(Int, refel.Nfp[1], Gness, totalfaces);  # face node local to global
    element2face = zeros(Int, nfaces, nel);  # element to face map
    face2element = zeros(Int, 2, size(mesh.face2element,2));  # face to element map
    facenormals = zeros(dim, totalfaces); # normal vectors for every face
    faceRefelInd = zeros(Int, 2, totalfaces); # Index in refel for this face for elements on both sides
    facebid = zeros(Int, totalfaces); # BID of each face

    tmpallnodes = zeros(dim, mesh.nel*refel.Np);
    t_nodes1 = Base.Libc.time();
    for ei=1:nel
        # Find this element's nodes
        n_vert = etypetonv[mesh.etypes[ei]];
        e_vert = mesh.nodes[1:dim, mesh.elements[1:n_vert, ei]];

        if dim == 1
            e_x = line_refel_to_x(refel.r[:,1], e_vert);
        elseif dim == 2
            if nfaces == 3 # triangles
                (e_x, e_y) = triangle_refel_to_xy_(refel.r[:,1], refel.r[:,2], e_vert);
            else # quads
                (e_x, e_y) = quad_refel_to_xy(refel.r[:,1], refel.r[:,2], e_vert);
            end
        elseif dim == 3
            if nvtx == 8 # hexes
                (e_x, e_y, e_z) = hex_refel_to_xyz(refel.r[:,1], refel.r[:,2], refel.r[:,3], e_vert);
            else # tets
                (e_x, e_y, e_z) = tetrahedron_refel_to_xyz(refel.r[:,1], refel.r[:,2], refel.r[:,3], e_vert);
            end
        end

        # Add them to the tmp global nodes
        tmpallnodes[1, ((ei-1)*Np+1):(ei*Np)] = e_x;
        if dim > 1
            tmpallnodes[2, ((ei-1)*Np+1):(ei*Np)] = e_y;
            if dim > 2
                tmpallnodes[3, ((ei-1)*Np+1):(ei*Np)] = e_z;
            end
        end

        # temporary mapping
        loc2glb[:,ei] = ((ei-1)*Np+1):(ei*Np);
    end
    t_nodes1 = Base.Libc.time() - t_nodes1;
    log_entry("Initial node setup: "*string(t_nodes1), 3);

    if Gness == 1 # CG
        # Go back and remove duplicate nodes. Adjust loc2glb.
        t_nodes1 = Base.Libc.time();
        (allnodes, loc2glb) = remove_duplicate_nodes(tmpallnodes, loc2glb, tol=tol);
        t_nodes1 = Base.Libc.time() - t_nodes1;
        log_entry("Remove duplicate nodes: "*string(t_nodes1), 3);

    else
        allnodes = tmpallnodes; # DG grid is already made
    end

    # vertices, faces and boundary
    t_faces1 = Base.Libc.time();
    for ei=1:nel
        n_vert = etypetonv[mesh.etypes[ei]];
        mfids = mesh.element2face[:,ei];
        normals = mesh.normals[:,mfids];
        el_center = zeros(dim);

        # vertices and center
        for ni=1:Np
            el_center += allnodes[:,loc2glb[ni,ei]];

            for vi=1:n_vert
                if is_same_node(mesh.nodes[:, mesh.elements[vi,ei]], allnodes[:,loc2glb[ni,ei]], tol)
                    glbvertex[vi, ei] = loc2glb[ni,ei];
                end
            end
        end
        el_center ./= Np;

        # f2glb has duplicates. Compare to mesh faces and keep same ordering as mesh.
        # Copy normals and bdry info.
        # Set element2face map.
        meshfaces = mesh.element2face[:,ei];    # index from mesh
        if dim == 1
            test_same_face = is_same_node;
        elseif dim == 2
            # test_same_face = is_same_line;
            test_same_face = is_same_face_center;
        elseif dim == 3
            # test_same_face = is_same_plane;
            test_same_face = is_same_face_center;
        end

        for gfi=1:nfaces
            tmpf2glb = loc2glb[refel.face2local[gfi], ei];

            for mfi=1:nfaces
                thisfaceind = meshfaces[mfi];
                if test_same_face(mesh.nodes[:,mesh.face2vertex[:,thisfaceind]], allnodes[:, tmpf2glb], tol)
                    # This mesh face corresponds to this tmpf2glb face
                    # Put the tmpf2glb map into f2glb at the mesh index(thisfaceind).
                    # Move the f2glb[:,1,ind] to f2glb[:,2,ind] first if DG (Gness==2)
                    # Set element2face according to gfi(not mfi)
                    # Move face2element[1] to face2element[2] and put this one in [1]
                    if (Gness == 2) f2glb[:, 2, thisfaceind] = f2glb[:, 1, thisfaceind]; end
                    f2glb[:, 1, thisfaceind] = tmpf2glb;
                    element2face[gfi, ei] = thisfaceind;
                    face2element[2, thisfaceind] = face2element[1, thisfaceind];
                    face2element[1, thisfaceind] = ei;

                    # Find the normal for every face. The normal points from e1 to e2 or outward for boundary.
                    # Note that the normal stored in mesh_data could be pointing either way.
                    thisnormal = normals[:, mfi];
                    f_center = zeros(dim);
                    for ni=1:length(tmpf2glb)
                        f_center += allnodes[:, tmpf2glb[ni]];
                    end
                    f_center ./= length(tmpf2glb)
                    d1 = norm(f_center + thisnormal - el_center);
                    d2 = norm(f_center - thisnormal - el_center);
                    fdotn = sum((f_center-el_center) .* thisnormal);
                    # if d1 < d2 # normal is pointing in the wrong direction
                    #     thisnormal = -thisnormal;
                    # end
                    if fdotn < 0 # normal is pointing in the wrong direction
                        thisnormal = -thisnormal;
                    end
                    facenormals[:, thisfaceind] = thisnormal;

                    # Copy boundary info: bdry, bdryface, bdrynorm
                    mbid = mesh.bdryID[thisfaceind];
                    gbid = indexin([mbid], bids)[1];
                    nfacenodes = length(tmpf2glb);
                    if !(gbid === nothing) # This is a boundary face
                        append!(bdry[gbid], tmpf2glb);
                        push!(bdryfc[gbid], thisfaceind);
                        facebid[thisfaceind] = gbid;
                        thisnormal = normals[:, mfi];
                        normchunk = zeros(config.dimension, nfacenodes);
                        for ni=1:nfacenodes
                            normchunk[:,ni] = thisnormal;
                        end
                        bdrynorm[gbid] = hcat(bdrynorm[gbid], normchunk);
                    end
                end
            end
        end

    end # element loop
    t_faces1 = Base.Libc.time() - t_faces1;
    log_entry("Vertex, face, bdry setup: "*string(t_faces1), 3);

    # There are duplicates in the bdry info. Remove them
    t_faces1 = Base.Libc.time();
    newbdry = similar(bdry);
    newbdrynorm = similar(bdrynorm);
    for i=1:length(bdry)
        newbdry[i] = zeros(length(bdry[i]));
        newbdrynorm[i] = zeros(size(bdrynorm[i]));
    end

    for bidi=1:length(bids)
        nextbdryind = 1;
        for bi=1:length(bdry[bidi])
            found = false;
            for bidj=1:bidi
                for bj=1:length(newbdry[bidj])
                    if bdry[bidi][bi] == newbdry[bidj][bj]
                        # bi already exists in newbdry
                        found = true;
                        break;
                    end
                end
                if found
                    break;
                end
            end
            if !found # it's a new bdry node
                newbdry[bidi][nextbdryind] = bdry[bidi][bi];
                newbdrynorm[bidi][:,nextbdryind] = bdrynorm[bidi][:,bi];
                nextbdryind += 1;
            end
        end
        newbdry[bidi] = newbdry[bidi][1:nextbdryind-1];
        newbdrynorm[bidi] = newbdrynorm[bidi][:, 1:nextbdryind-1];
    end
    bdry = newbdry;
    bdrynorm = newbdrynorm;

    t_faces1 = Base.Libc.time() - t_faces1;
    log_entry("Remove duplicate bdry nodes: "*string(t_faces1), 3);

    # Refel index for each face
    for fi=1:totalfaces
        faceRefelInd[1,fi] = which_refel_face(f2glb[:,1,fi], allnodes, refel, loc2glb[:,face2element[1,fi]]);
        if face2element[2,fi] > 0
            faceRefelInd[2,fi] = which_refel_face(f2glb[:,Gness,fi], allnodes, refel, loc2glb[:,face2element[2,fi]]);
        end
    end

    t_grid_from_mesh = Base.Libc.time() - t_grid_from_mesh;
    log_entry("Total grid building time: "*string(t_grid_from_mesh), 2);

    return (refel, Grid(allnodes, bdry, bdryfc, bdrynorm, bids, loc2glb, glbvertex, f2glb, element2face, face2element, facenormals, faceRefelInd, facebid));
end

function collectBIDs(mesh)
    bids = [];
    for i=1:length(mesh.bdryID)
        if mesh.bdryID[i] > 0
            already = false;
            for j=1:length(bids)
                if mesh.bdryID[i] == bids[j]
                    already = true;
                end
            end
            if !already
                push!(bids, mesh.bdryID[i]);
            end
        end
    end
    return bids;
end

# Removes duplicate nodes and updates local to global maps
function remove_duplicate_nodes(nodes, loc2glb; tol=1e-12, depth=5, mincount=50, other2glb=[])
    # defaults
    # depth = 5; # 32768 for 3D
    # mincount = 50; # don't subdivide if less than this

    # Strategy: divide nodes into bins compare against nodes in bin
    tmpNnodes = size(nodes,2);
    dim = size(nodes,1);
    replace_with = zeros(Int, tmpNnodes); # holds index of replacement node being kept. 0 if this is kept.
    abins = zeros(Int, tmpNnodes);
    bin_ends = [tmpNnodes]; # The last index of each bin

    # init bins and find extrema
    xlim = [1e10, -1e10];
    ylim = [1e10, -1e10];
    zlim = [1e10, -1e10];
    for i=1:tmpNnodes
        abins[i] = i;
        xyz=nodes[:,i];
        xlim[1] = min(xlim[1], xyz[1]);
        xlim[2] = max(xlim[2], xyz[1]);
        if dim>1
            ylim[1] = min(ylim[1], xyz[2]);
            ylim[2] = max(ylim[2], xyz[2]);
            if dim > 2
                zlim[1] = min(zlim[1], xyz[3]);
                zlim[2] = max(zlim[2], xyz[3]);
            end
        end
    end

    if dim == 1
        (abins, bin_ends) = partition_nodes_in_bins_1d(nodes, abins, xlim, depth, mincount);
    elseif dim == 2
        (abins, bin_ends) = partition_nodes_in_bins_2d(nodes, abins, xlim, ylim, depth, mincount);
    else # dim==3
        (abins, bin_ends) = partition_nodes_in_bins_3d(nodes, abins, xlim, ylim, zlim, depth, mincount);
    end

    remove_count = 0;
    startind = 1;
    # Loop over nodes in bins and put numbers in replace_with where duplicates will be removed.
    for bini = 1:length(bin_ends)
        for ni=startind:bin_ends[bini]
            for nj=startind:(ni-1)
                # If node[nj] == node[ni], keep nj, remove ni
                if is_same_node(nodes[:,abins[ni]], nodes[:,abins[nj]], tol)
                    remove_count += 1;
                    replace_with[abins[ni]] = abins[nj];
                    break;
                end
            end
        end
        startind = bin_ends[bini] + 1;
    end

    # println(string(length(bin_ends))*" bins: " )
    # println(bin_ends);
    # println("Nnodes before: "*string(tmpNnodes)*" after: "*string(tmpNnodes-remove_count));

    # Loop over nodes and place in new allnodes array while adjusting replace_with
    Nnodes = tmpNnodes-remove_count;
    next_ind = 1;
    newnodes = zeros(size(nodes,1), Nnodes);
    new_homes = zeros(Int, tmpNnodes);
    for i=1:tmpNnodes
        if replace_with[i] == 0
            # A node to keep
            newnodes[:,next_ind] = nodes[:,i];
            new_homes[i] = next_ind;
            next_ind += 1;
        end
    end
    for i=1:tmpNnodes
        if replace_with[i] > 0
            # A node that was removed. update replace_with with new_homes
            replace_with[i] = new_homes[replace_with[i]];
        end
    end

    # Loop over each element's nodes and adjust loc2glb
    for ei=1:size(loc2glb,2) # nel
        for ni=1:size(loc2glb,1) # Np
            ind = loc2glb[ni,ei];
            if replace_with[ind] > 0
                loc2glb[ni,ei] = replace_with[ind];
            else
                loc2glb[ni,ei] = new_homes[ind];
            end
        end
    end

    # If another 2glb map is present, do the same to it
    if length(other2glb) > 0
        for i=1:length(other2glb) #
            ind = other2glb[i];
            if replace_with[ind] > 0
                other2glb[i] = replace_with[ind];
            else
                other2glb[i] = new_homes[ind];
            end
        end

        return (newnodes, loc2glb, other2glb);
    end

    return (newnodes, loc2glb);
end

# recursively subdivide bins of nodes
function partition_nodes_in_bins_1d(nodes, abin, xlim, depth, mincount)
    halfx = (xlim[1] + xlim[2])/2;
    N = length(abin);

    bbin = similar(abin); # Temporary storage
    lcount = 0;
    rcount = 0;
    for i=1:N
        if nodes[1,abin[i]] <= halfx
            lcount += 1;
            bbin[lcount] = abin[i];
        else
            rcount += 1;
            bbin[N-rcount+1] = abin[i];
        end
    end

    # Only recurse if depth and mincount allow
    if depth > 0
        if lcount >= mincount
            (abin[1:lcount], lbinends) = partition_nodes_in_bins_1d(nodes, bbin[1:lcount], [xlim[1], halfx], depth-1, mincount);
        else
            abin[1:lcount] = bbin[1:lcount];
            lbinends = [lcount];
        end
        if rcount >= mincount
            (abin[(lcount+1):N], rbinends) = partition_nodes_in_bins_1d(nodes, bbin[(lcount+1):N], [halfx, xlim[2]], depth-1, mincount);
        else
            abin[(lcount+1):N] = bbin[(lcount+1):N];
            rbinends = [rcount];
        end
        binends = append!(lbinends, rbinends .+ lcount);

    else
        abin = bbin;
        binends = [lcount, N];
    end

    return (abin, binends);
end

# recursively subdivide bins of nodes
function partition_nodes_in_bins_2d(nodes, abin, xlim, ylim, depth, mincount)
    halfx = (xlim[1] + xlim[2])/2;
    halfy = (ylim[1] + ylim[2])/2;
    N = length(abin);

    which_bin = zeros(Int, N) # Which bin will the node go in
    bin_counts = zeros(Int, 4) # How many go in each
    for i=1:N
        if nodes[1,abin[i]] <= halfx
            if nodes[2,abin[i]] <= halfy
                bin_counts[1] += 1;
                which_bin[i] = 1;
            else
                bin_counts[2] += 1;
                which_bin[i] = 2;
            end

        else
            if nodes[2,abin[i]] <= halfy
                bin_counts[3] += 1;
                which_bin[i] = 3;
            else
                bin_counts[4] += 1;
                which_bin[i] = 4;
            end
        end
    end

    # Make bbin arrays and fill them.
    bbins = [zeros(Int, bin_counts[1]), zeros(Int, bin_counts[2]), zeros(Int, bin_counts[3]), zeros(Int, bin_counts[4])];
    next_inds = [1,1,1,1];
    for i=1:N
        ind = which_bin[i];
        bbins[ind][next_inds[ind]] = abin[i];
        next_inds[ind] += 1;
    end

    # convenient offsets and limits
    starts = [1, 0,0,0];
    ends   = [bin_counts[1], 0,0,0];
    for i=2:4
        starts[i] = ends[i-1] + 1;
        ends[i] = bin_counts[i] + ends[i-1];
    end
    xlims = [[xlim[1], halfx], [xlim[1], halfx], [halfx, xlim[2]], [halfx, xlim[2]]];
    ylims = [[ylim[1], halfy], [halfy, ylim[2]], [ylim[1], halfy], [halfy, ylim[2]]];

    # Only recurse if depth and mincount allow
    if depth > 0
        binends = [];
        for i=1:4
            if bin_counts[i] >= mincount
                (abin[starts[i]:ends[i]], tmpbinends) = partition_nodes_in_bins_2d(nodes, bbins[i], xlims[i], ylims[i], depth-1, mincount);
                if i>1
                    append!(binends, tmpbinends .+ ends[i-1]);
                else
                    binends = tmpbinends;
                end

            else
                abin[starts[i]:ends[i]] = bbins[i];
                append!(binends, [ends[i]]);
            end
        end

    else
        abin = [bbins[1]; bbins[2]; bbins[3]; bbins[4]];
        binends = ends;
    end

    return (abin, binends);
end

# recursively subdivide bins of nodes
function partition_nodes_in_bins_3d(nodes, abin, xlim, ylim, zlim, depth, mincount)
    halfx = (xlim[1] + xlim[2])/2;
    halfy = (ylim[1] + ylim[2])/2;
    halfz = (zlim[1] + zlim[2])/2;
    N = length(abin);

    which_bin = zeros(Int, N) # Which bin will the node go in
    bin_counts = zeros(Int, 8) # How many go in each
    for i=1:N
        if nodes[1,abin[i]] <= halfx
            if nodes[2,abin[i]] <= halfy
                if nodes[2,abin[i]] <= halfz
                    bin_counts[1] += 1;
                    which_bin[i] = 1;
                else
                    bin_counts[2] += 1;
                    which_bin[i] = 2;
                end
            else
                if nodes[2,abin[i]] <= halfy
                    bin_counts[3] += 1;
                    which_bin[i] = 3;
                else
                    bin_counts[4] += 1;
                    which_bin[i] = 4;
                end
            end

        else
            if nodes[2,abin[i]] <= halfy
                if nodes[2,abin[i]] <= halfz
                    bin_counts[5] += 1;
                    which_bin[i] = 5;
                else
                    bin_counts[6] += 1;
                    which_bin[i] = 6;
                end
            else
                if nodes[2,abin[i]] <= halfy
                    bin_counts[7] += 1;
                    which_bin[i] = 7;
                else
                    bin_counts[8] += 1;
                    which_bin[i] = 8;
                end
            end
        end
    end

    # Make bbin arrays and fill them.
    bbins = [zeros(Int, bin_counts[1]), zeros(Int, bin_counts[2]), zeros(Int, bin_counts[3]), zeros(Int, bin_counts[4]),
             zeros(Int, bin_counts[5]), zeros(Int, bin_counts[6]), zeros(Int, bin_counts[7]), zeros(Int, bin_counts[8])];
    next_inds = [1,1,1,1,1,1,1,1];
    for i=1:N
        ind = which_bin[i];
        bbins[ind][next_inds[ind]] = abin[i];
        next_inds[ind] += 1;
    end

    # convenient offsets and limits
    starts = [1, 0,0,0,0,0,0,0];
    ends   = [bin_counts[1], 0,0,0,0,0,0,0];
    for i=2:8
        starts[i] = ends[i-1] + 1;
        ends[i] = bin_counts[i] + ends[i-1];
    end
    xlims = [[xlim[1], halfx], [xlim[1], halfx], [xlim[1], halfx], [xlim[1], halfx],
             [halfx, xlim[2]], [halfx, xlim[2]], [halfx, xlim[2]], [halfx, xlim[2]]];
    ylims = [[ylim[1], halfy], [ylim[1], halfy], [halfy, ylim[2]], [halfy, ylim[2]],
             [ylim[1], halfy], [ylim[1], halfy], [halfy, ylim[2]], [halfy, ylim[2]]];
    zlims = [[zlim[1], halfz], [halfz, zlim[2]], [zlim[1], halfz], [halfz, zlim[2]],
             [zlim[1], halfz], [halfz, zlim[2]], [zlim[1], halfz], [halfz, zlim[2]]];

    # Only recurse if depth and mincount allow
    if depth > 0
        binends = [];
        for i=1:8
            if bin_counts[i] >= mincount
                (abin[starts[i]:ends[i]], tmpbinends) = partition_nodes_in_bins_3d(nodes, bbins[i], xlims[i], ylims[i], zlims[i], depth-1, mincount);
                if i>1
                    append!(binends, tmpbinends .+ ends[i-1]);
                else
                    binends = tmpbinends;
                end
            else
                abin[starts[i]:ends[i]] = bbins[i];
                append!(binends, [ends[i]]);
            end
        end

    else
        abin = [bbins[1]; bbins[2]; bbins[3]; bbins[4]; bbins[5]; bbins[6]; bbins[7]; bbins[8]];
        binends = ends;
    end

    return (abin, binends);
end

#Extra remove later

function triangle_element_nodes_(refel, v)
    return  triangle_refel_to_xy_(refel.r[:,1], refel.r[:,2], v);
end

function line_refel_to_x(r, v)
    x = v[1] .+ 0.5 .* (v[2]-v[1]) .* (1 .+ r);

    return x;
end

function triangle_refel_to_xy_(r, s, v)
    x = 0.5 * (-(r .+ s) * v[1,1] .+ (1 .+ r) * v[1,2] .+ (1 .+ s) * v[1,3]);
    y = 0.5 * (-(r .+ s) * v[2,1] .+ (1 .+ r) * v[2,2] .+ (1 .+ s) * v[2,3]);

    return (x, y);
end

function quad_refel_to_xy(r, s, v)
    vx = v[1,:];
    vy = v[2,:];

    rp = 1 .+ r;
    rm = 1 .- r;
    sp = 1 .+ s;
    sm = 1 .- s;
    x = 0.25 .* (rm .* sm .* vx[1] .+ rp .* sm .* vx[2] .+ rp .* sp .* vx[3] .+ rm .* sp .* vx[4]);
    y = 0.25 .* (rm .* sm .* vy[1] .+ rp .* sm .* vy[2] .+ rp .* sp .* vy[3] .+ rm .* sp .* vy[4]);

    return (x, y);
end

function hex_refel_to_xyz(r, s, t, v)
    vx = v[1,:];
    vy = v[2,:];
    vz = v[3,:];

    rp = 1 .+ r;
    rm = 1 .- r;
    sp = 1 .+ s;
    sm = 1 .- s;
    tp = 1 .+ t;
    tm = 1 .- t;
    x = 0.125 .* (rm .* sm .* tm .* vx[1] .+ rp .* sm .* tm .* vx[2] .+ rp .* sp .* tm .* vx[3] .+ rm .* sp .* tm .* vx[4]
                    .+ rm .* sm .* tp .* vx[5] .+ rp .* sm .* tp .* vx[6] .+ rp .* sp .* tp .* vx[7] .+ rm .* sp .* tp .* vx[8]);
    y = 0.125 .* (rm .* sm .* tm .* vy[1] .+ rp .* sm .* tm .* vy[2] .+ rp .* sp .* tm .* vy[3] .+ rm .* sp .* tm .* vy[4]
                    .+ rm .* sm .* tp .* vy[5] .+ rp .* sm .* tp .* vy[6] .+ rp .* sp .* tp .* vy[7] .+ rm .* sp .* tp .* vy[8]);
    z = 0.125 .* (rm .* sm .* tm .* vz[1] .+ rp .* sm .* tm .* vz[2] .+ rp .* sp .* tm .* vz[3] .+ rm .* sp .* tm .* vz[4]
                    .+ rm .* sm .* tp .* vz[5] .+ rp .* sm .* tp .* vz[6] .+ rp .* sp .* tp .* vz[7] .+ rm .* sp .* tp .* vz[8]);

    return (x, y, z);
end

function tetrahedron_refel_to_xyz(r, s, t, v)
    d = v[:,1];
    A = [v[:,2].-d v[:,3].-d v[:,4].-d];

    np = length(r);
    mv = zeros(3,np);
    for i=1:np
        tmp = [(r[i]+1)/2, (s[i]+1)/2, (t[i]+1)/2];
        mv[:,i] = A*tmp + d;
    end

    x = mv[1,:]
    y = mv[2,:]
    z = mv[3,:]

    return (x, y, z);
end

# Returns true if the nodes are within tol of each other.
function is_same_node(x1, x2, tol)
    return sum(abs.(x1 - x2)) < tol
end

# Returns true if the centroid of each line is close enough.
function is_same_face_center(l1, l2, tol)
    n1 = size(l1,2);
    n2 = size(l2,2);
    center1 = zeros(size(l1,1));
    center2 = zeros(size(l2,1));
    for i=1:n1
        center1 = center1 .+ l1[:,i];
    end
    center1 = center1 ./ n1;

    for i=1:n2
        center2 = center2 .+ l2[:,i];
    end
    center2 = center2 ./ n2;

    return is_same_node(center1, center2, tol);
end

# Returns true if the two node lists have at least two of the same nodes.
function is_same_line(l1, l2, tol)
    found = 0;
    n1 = size(l1,2);
    n2 = size(l2,2);
    for i=1:n1
        for j=1:n2
            if is_same_node(l1[:,i], l2[:,j], tol)
                found += 1;
            end
            if found >= 2
                return true;
            end
        end
    end

    return false;
end

# Returns true if the two node lists have at least three of the same nodes.
function is_same_plane(p1, p2, tol)
    found = 0;
    n1 = size(p1,2);
    n2 = size(p2,2);
    for i=1:n1
        for j=1:n2
            if is_same_node(p1[:,i], p2[:,j], tol)
                found += 1;
            end
            if found >= 3
                return true;
            end
        end
    end

    return false;
end

function which_refel_face(f2glb, nodes, ref, elem)
    # Check f2glb against the face2local in refel
    fnodes = nodes[:,f2glb];
    for fi=1:ref.Nfaces
        refnodes = nodes[:, elem[ref.face2local[fi]]];

        if is_same_face(fnodes, refnodes, ref.dim)
            return fi;
        end
    end

    printerr("Couldn't match face when building grid (see which_refel_face() in grid.jl)");
end

function is_same_face(f1, f2, dim)
    tol = 1e-12;
    if dim == 1 # one point
        return is_same_node(f1, f2, tol);
    elseif dim == 2 # same line(two same points)
        return is_same_line(f1, f2, tol);
    elseif dim == 3 # same plane(three same points)
        return is_same_plane(f1, f2, tol);
    end
end

# Adds a boundary ID to some region. Find boundary points satifying on_bdry and moves them to a new set for this bid.
function add_boundary_ID_to_grid(bid, on_bdry, grid)
    # Find if this bid exists. If so, just add points to it, removing from others.
    ind = indexin([bid], grid.bids)[1];
    nbids = length(grid.bids);
    if ind === nothing
        # This is a new bid, add to bids, bdry, bdryface, bdrynorm
        ind = nbids + 1;
        nbids += 1;
        push!(grid.bids, bid);
        push!(grid.bdry, zeros(Int, 0));
        push!(grid.bdryface, zeros(Int, 0));
        push!(grid.bdrynorm, zeros(config.dimension, 0));
    end

    # Search all other bids for nodes and faces on this segment. Remove them there and add them here.
    # First find indices and count them. Then move.
    move_nodes = Array{Array{Int,1},1}(undef,nbids);
    node_count = zeros(Int, nbids);
    move_faces = Array{Array{Int,1},1}(undef,nbids);
    face_count = zeros(Int, nbids);
    for i=1:nbids
        bi = grid.bids[i];
        move_nodes[i] = [];
        move_faces[i] = [];
        if bi != bid
            # First the nodes
            for j=1:length(grid.bdry[i])
                nj = grid.bdry[i][j];
                if config.dimension == 1
                    if on_bdry(grid.allnodes[1, nj])
                        push!(move_nodes[i], nj);
                        node_count[i] += 1;
                    end
                elseif config.dimension == 2
                    if on_bdry(grid.allnodes[1, nj], grid.allnodes[2, nj])
                        push!(move_nodes[i], nj);
                        node_count[i] += 1;
                    end
                elseif config.dimension == 3
                    if on_bdry(grid.allnodes[1, nj], grid.allnodes[2, nj], grid.allnodes[3, nj])
                        push!(move_nodes[i], nj);
                        node_count[i] += 1;
                    end
                end
            end
            # Then the faces
            # Base this decision on the center of the face
            for j=1:length(grid.bdryface[i])
                fj = grid.bdryface[i][j];
                nfp = size(grid.face2glb,1)
                isbdryface = true
                # find the center
                fcenter = zeros(size(grid.allnodes,1));
                for ni=1:nfp
                    fcenter = fcenter + grid.allnodes[:,grid.face2glb[ni,1,fj]];
                end
                fcenter = fcenter./nfp;

                if config.dimension == 1
                    isbdryface = on_bdry(fcenter[1]);
                elseif config.dimension == 2
                    isbdryface = on_bdry(fcenter[1], fcenter[2]);
                elseif config.dimension == 3
                    isbdryface = on_bdry(fcenter[1], fcenter[2], fcenter[3]);
                end

                if isbdryface
                    push!(move_faces[i], fj);
                    face_count[i] += 1;
                end
            end
        end
    end # find indices

    # Move things from other bids to this one
    for i=1:nbids
        if i != ind
            # Add to this bid
            append!(grid.bdry[ind], move_nodes[i]);
            append!(grid.bdryface[ind], move_faces[i]);
            grid.bdrynorm[ind] = hcat(grid.bdrynorm[ind], grid.bdrynorm[i][:,indexin(move_nodes[i], grid.bdry[i])]);

            # Make sure all of the norms correspond to the face on this bdry
            startnodeind = length(grid.bdry[ind]) - length(move_nodes[i]);
            for ni=1:length(move_nodes[i])
                for fi=1:length(move_faces[i])
                    # Does this node lie on this face?
                    facenodeindex = indexin(move_nodes[i][ni], grid.face2glb[:,1,move_faces[i][fi]]);
                    if length(facenodeindex) > 0
                        grid.bdrynorm[ind][:,startnodeind + ni] = grid.facenormals[:,move_faces[i][fi]];
                        break;
                    end
                end
            end

            # Remove things from other bids
            # Remove bdrynorm and bdryfacenorm
            numremove = length(move_nodes[i]);
            if numremove > 0
                newbdrynorm = zeros(config.dimension, size(grid.bdrynorm[i],2) - numremove);
                nextind = 1;
                for j=1:length(grid.bdry[i])
                    keepit = true;
                    for k=1:numremove
                        if grid.bdry[i][j] == move_nodes[i][k]
                            keepit = false;
                            break;
                        end
                    end
                    if keepit
                        newbdrynorm[:,nextind] = grid.bdrynorm[i][:,j];
                        nextind += 1;
                    end
                end
                grid.bdrynorm[i] = newbdrynorm;
            end

            # Remove nodes
            deleteat!(grid.bdry[i], indexin(move_nodes[i], grid.bdry[i]));

            # Remove bdryface
            deleteat!(grid.bdryface[i], indexin(move_faces[i], grid.bdryface[i]));

            # Change the bid of moved faces in facebid
            for fi=1:length(move_faces[i])
                grid.facebid[move_faces[i][fi]] = bid;
            end
        end
    end

    log_entry("Added boundary ID: "*string(bid)*" including "*string(node_count)*" nodes, "*string(face_count)*" faces.");

end


text:
using Pkg

isinstalled(pkg::String) = any(x -> x.name == pkg && x.is_direct_dep, values(Pkg.dependencies()))

# Install required packages
for P in ["HTTP", "PyPlot"]
  !isinstalled(P) && Pkg.add(P)
end

# Download data
include("downloadData.jl")

using PyPlot, MRIReco

##################
## Data Loading ##
##################
f = ISMRMRDFile("data/ksp_knee_3dfse.h5")
acqRaw = RawAcquisitionData(f)

# kspace-center is not specified in the raw data file.
# So we let MRIReco estimate it
acqData = AcquisitionData(acqRaw; estimateProfileCenter=true)

# convert to 2d
acqData2d = convert3dTo2d(acqData)

# extract slices
sl = [50,100,150,200]
acqData2d.kdata = acqData2d.kdata[:,sl,:]
acqData2d.encodingSize = [274,208,1]

####################################
## generate coil sensitivity maps ##
####################################
@info "Espirit"
smaps = espirit(acqData2d, (6,6), 30, eigThresh_1=0.04, eigThresh_2=0.98)

###############################
## Reconstruction Parameters ##
###############################
params = Dict{Symbol, Any}()
params[:reco] = "multiCoil"
params[:reconSize] = Tuple(acqData2d.encodingSize[1:2])
params[:solver] = "admm"
params[:regularization] = "L1"
params[:sparseTrafo] = "Wavelet"
params[:λ] = 2.e-1
params[:iterations] = 30
params[:ρ] = 0.1
params[:absTol] = 1.e-4
params[:relTol] = 1.e-3
params[:tolInner] = 1.e-2
params[:senseMaps] = smaps
params[:normalizeReg] = true

############################
## Perform reconstruction ##
############################

@info "Perform Reconstruction "
@time img = fftshift( reconstruction(acqData2d, params).data[:,:,:,1,1], 2)

###########################
## Visualize the results ##
###########################

figure(1, figsize=(6,1.9))

for s in [1,2,3,4]
  subplot(1,4,s)
  imshow(abs.(img[:,:,s,1,1]))
  title("slice $(sl[s])")
  axis("off")
end

subplots_adjust(left=0.01,bottom=0.01,wspace=0.3,hspace=0.2,right=0.99,top=0.92)


text:
using NFFT, DataFrames, LinearAlgebra, LaTeXStrings, DelimitedFiles
using BenchmarkTools
using Plots; pgfplotsx()

const benchmarkTime = [4, 4]

NFFT.FFTW.set_num_threads(Threads.nthreads())
NFFT._use_threads[] = (Threads.nthreads() > 1)

const threads = [1,2,4,8]

#const NBase = [65536, 256, 32]
const NBase = [512*512, 512, 64]
const Ds = 1:3
const blockSizeBase = [
  [nextpow(2,round(Int,x)) for x=2 .^range(1,19,length=7)], #[div(NBase[1]*2,2^j) for j=(0:6).*2],
  [nextpow(2,round(Int,x)) for x=2 .^range(1,10,length=7)], #[div(NBase[2]*2,2^j) for j=0:6],
  [nextpow(2,round(Int,x)) for x=2 .^range(1,7,length=7)] #,[div(NBase[3]*2,2^j) for j=0:6],
]

const fftflags = NFFT.FFTW.MEASURE

function nfft_block_size_comparison(Ds=1:3, m = 4, σ = 2.0)
  println("\n\n ##### nfft_block_size_comparison ##### \n\n")

  df = DataFrame(Package=String[], Threads=Int[], D=Int[], M=Int[], N=Int[], m = Int[], σ=Float64[],
                   TimeTrafo=Float64[], TimeAdjoint=Float64[], blocking=Bool[], blockSizeBase=Int[] )

  for D in Ds
    @info "### Dimension D=$D ###"
    N = ntuple(d->NBase[D], D)
    M = prod(N)

    x = rand(D,M) .- 0.5
    fHat = randn(ComplexF64, M)
    fApprox = randn(ComplexF64, N)
    gHatApprox = randn(ComplexF64, M)

    for bl in blockSizeBase[D]
        @info "b=$(bl) D=$D"

          blockSize = ntuple(d-> bl, D)
          p = NFFTPlan(x, N; m, σ, precompute=POLYNOMIAL, blocking=true, blockSize=blockSize, fftflags=fftflags)

          @info "Adjoint benchmark:"
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[1]
          b = @benchmark mul!($fApprox, $(adjoint(p)), $fHat)
          tadjoint = minimum(b).time / 1e9

          @info "Trafo benchmark:"
          BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[2]
          b = @benchmark mul!($gHatApprox, $p, $fApprox)
          ttrafo = minimum(b).time / 1e9

          push!(df, ("NFFT.jl", Threads.nthreads(), D, M, N[D], m, σ, ttrafo, tadjoint, true, bl))
    end

    # non blocking
    p = NFFTPlan(x, N; m, σ, precompute=POLYNOMIAL, blocking=false)

    @info "Adjoint benchmark:"
    BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[1]
    b = @benchmark mul!($fApprox, $(adjoint(p)), $fHat)
    tadjoint = minimum(b).time / 1e9

    @info "Trafo benchmark:"
    BenchmarkTools.DEFAULT_PARAMETERS.seconds = benchmarkTime[2]
    b = @benchmark mul!($gHatApprox, $p, $fApprox)
    ttrafo = minimum(b).time / 1e9

    push!(df, ("NFFT.jl", Threads.nthreads(), D, M, N[D], m, σ, ttrafo, tadjoint, false, 0))
  end
  return df
end

function plot_performance_block_size(df, Ds)

  Plots.scalefontsizes()
  Plots.scalefontsizes(1.5)

  colors = [:black, :orange, :blue, :green, :brown, :gray, :blue, :purple, :yellow ]
  ls = [:solid, :solid, :solid, :solid, :solid, :dash, :solid, :dash, :solid]
  shape = [:circle, :circle, :circle, :xcross, :circle, :xcross, :xcross, :circle]

  pl = Matrix{Any}(undef, 2, length(Ds))
  for (i,D) in enumerate(Ds)
    titleTrafo = L"\textrm{NFFT}, \textrm{%$(D)D}"
    titleAdjoint = L"\textrm{NFFT}^H, \textrm{%$(D)D}"

    df1_ = df[df.σ.==2.0 .&& df.D.==D .&& df.blocking .== true ,:]
    df2_ = df[df.σ.==2.0 .&& df.D.==D .&& df.blocking .== false ,:]
    maxTimeTrafo = max(maximum(df1_[:,:TimeTrafo]), maximum(df2_[:,:TimeTrafo]))*1.1
    maxTimeAdjoint = max(maximum(df1_[:,:TimeAdjoint]),maximum(df2_[:,:TimeAdjoint]))*1.1

    blockSizes = df1_[df1_.Threads .== threads[1], :blockSizeBase]

    p1 = plot(blockSizes, yscale = :log2,
              df1_[df1_.Threads .== threads[1], :TimeTrafo], #ylims=(0.0,maxTimeTrafo),
              label="1 thread", lw=2, ylabel="Runtime / s", xlabel = i==length(Ds) ? "Block Size" : "",
              legend = nothing, title=titleTrafo, shape=:circle, c=:black)

    for p=2:length(threads)
      plot!(p1, blockSizes,
            df1_[df1_.Threads .== threads[p],:TimeTrafo],
              label="$(threads[p]) threads", lw=2, shape=shape[p], ls=ls[p],
              c=colors[p], msc=colors[p], mc=colors[p], ms=4, msw=2)
    end

    plot!(p1, blockSizes,
              df2_[df2_.Threads .== threads[1], :TimeTrafo].*ones(length(blockSizes)), #ylims=(0.0,maxTimeTrafo),
              label="regular: 1 thread", lw=2,  ls=:dash,
              c=colors[1])

    for p=2:length(threads)
      plot!(p1, blockSizes,
              df2_[df2_.Threads .== threads[p], :TimeTrafo].*ones(length(blockSizes)),
              label="regular: $(threads[p]) threads", lw=2,  ls=:dash,
              c=colors[p])
    end

    p2 = plot(blockSizes, yscale = :log2,
              df1_[df1_.Threads .== threads[1],:TimeAdjoint], #ylims=(0.0,maxTimeAdjoint),
              lw=2, xlabel = i==length(Ds) ? "Block Size" : "", label="1 thread",
              legend = i==2 ? :topright : nothing, title=titleAdjoint, shape=:circle, c=:black)

    for p=2:length(threads)
      plot!(p2, blockSizes,
            df1_[df1_.Threads .== threads[p],:TimeAdjoint],
              label="$(threads[p]) threads", lw=2, shape=shape[p], ls=ls[p],
              c=colors[p], msc=colors[p], mc=colors[p], ms=4, msw=2)
    end

    plot!(p2, blockSizes,
              df2_[df2_.Threads .== threads[1], :TimeAdjoint].*ones(length(blockSizes)), #ylims=(0.0,maxTimeAdjoint),
              lw=2, ls=:dash, label=nothing,
              c=colors[1])
    pl[1,i] = p1; pl[2,i] = p2;
  end

  #p = plot(p1, p2, layout=(1,2), size=(800,300), dpi=200)
  p = plot(vec(pl)..., layout=(length(Ds),2), size=(800,700), dpi=200)

  mkpath("./img/")
  savefig(p, "./img/performanceBlockSize.pdf")
  return p
end

### run the code ###

if haskey(ENV, "NFFT_PERF")
  df = nfft_block_size_comparison(Ds)

  if isfile("./data/performanceBlockSize.csv")
    data, header = readdlm("./data/performanceBlockSize.csv", ',', header=true);
    df_ = DataFrame(data, vec(header))
    append!(df, df_)
  end

  writedlm("./data/performanceBlockSize.csv", Iterators.flatten(([names(df)], eachrow(df))), ',')

else
  if false
    rm("./data/performanceBlockSize.csv", force=true)
    ENV["NFFT_PERF"] = 1
    for t in threads
      cmd = `julia -t $t performanceBlockSize.jl`
      @info cmd
      run(cmd)
    end

    delete!(ENV, "NFFT_PERF")
  end

  data, header = readdlm("./data/performanceBlockSize.csv", ',', header=true);
  df = DataFrame(data, vec(header))

  plot_performance_block_size(df, Ds)
end


text:
# -----------------------------------------------------------
#
# Using boundary terms and integration (summation) by parts
# to construct more accurate derivative matrices for time
# evolution using Lax-Friedricks-type fluxes
#
# All derivatives are computed here in pos basis first
# Then conjugated into hier basis
#
# -----------------------------------------------------------

# Efficiency criticality: MEDIUM
# Computations performed once,
# but take very long in multidimensional setting

# Accuracy criticality: HIGH
# Critical for accurate PDE evolution

# After the integration by parts, this corresponds to
# the remaining integral, ignoring the boundary term
function D_matrix(k::Int, level::Int)
    i = 1
    I = Int[]
    J = Int[]
    V = Float64[]
    for cell1 in 1:(1<<level)
        for mode1 in 1:k
            j = 1
            for cell2 in 1:(1<<level)
                for mode2 in 1:k
                    val = legvDv(level, cell1, mode1, cell2, mode2)
                    if abs(val) > 1.0e-15
                        push!(I, i)
                        push!(J, j)
                        push!(V, val)
                    end
                    j += 1
                end
            end
            i += 1
        end
    end
    return sparse(I, J, V, k * (1<<level), k * (1<<level), +)
end

# 1D D_matrix for anisotropic grids
function D_matrix(k::Int, level::Vector{Int})
    i = 1
    I = Int[]
    J = Int[]
    V = Float64[]
    for cell1 in 1:(1<<level) # 1<<level = 1*2^level also Anzahl Zellen bei dem level
        for mode1 in 1:k
            j = 1
            for cell2 in 1:(1<<level)
                for mode2 in 1:k
                    val = legvDv(level, cell1, mode1, cell2, mode2) # if cell1=cell2 calculate else 0.0
                    if abs(val) > 1.0e-15
                        push!(I, i)
                        push!(J, j)
                        push!(V, val)
                    end
                    j += 1
                end
            end
            # for mode2 in 1:k
            #     val = legvDv(level,cell1,mode1,cell1,mode2)
            #     push!(I, i)
            #     push!(J, j)
            #     push!(V, val)
            #     j += 1
            # end
            i += 1
        end
    end
    return sparse(I, J, V, k * (1<<level), k * (1<<level), +)
end

# -----------------------------------------------------
# Lax-Friedrichs flux matrix element on Legendre basis
# For now, use alpha = 0 only
# -----------------------------------------------------

# Periodic boundary term:
function periodic_legvLFv(level::Int, cell1::Int, mode1::Int,
                                      cell2::Int, mode2::Int; alpha::Real = 0)
    point1 = (cell2-1)//(1<<level)
    point2 = (cell2)//(1<<level)
    tiny   = 5.0e-16

    left1  = basis(level, cell1, mode1, point1-tiny)
    right1 = basis(level, cell1, mode1, point1+tiny)
    left2  = basis(level, cell1, mode1, point2-tiny)
    right2 = basis(level, cell1, mode1, point2+tiny)

    (cell2 == (1<<level)) && (right2 = basis(level, cell1, mode1, 0.0+tiny))
    (cell2 == 1)          && (left1  = basis(level, cell1, mode1, 1.0-tiny))

    LF1 = 1//2 * (left1 + right1) + alpha * (right1 - left1)
    LF2 = 1//2 * (left2 + right2) + alpha * (right2 - left2)

    val1 = basis(level, cell2, mode2, point1+tiny)
    val2 = basis(level, cell2, mode2, point2-tiny)

    # Net flux:
    return LF2 * val2 - LF1 * val1
end

function periodic_LF_matrix(k::Int, level::Int; alpha::Real = 0)
    i = 1
    I = Int[]
    J = Int[]
    V = Float64[]
    for cell1 in 1:(1<<level)
        for mode1 in 1:k
            j=1
            for cell2 in 1:(1<<level)
                for mode2 in 1:k
                    val = periodic_legvLFv(level, cell1, mode1,
                                                  cell2, mode2; alpha=alpha)
                    if abs(val)>1.0e-15
                        push!(I, i)
                        push!(J, j)
                        push!(V, val)
                    end
                    j += 1
                end
            end
            i += 1
        end
    end
    return sparse(I, J, V, k * (1<<level), k * (1<<level), +)
end

# ---------------------------------------------------------
# We can now build the full derivative matrix
# in both the position and hierarchical basis
# with some help from the hier2pos change of basis matrix
# ---------------------------------------------------------

function periodic_pos_DLF_matrix(k::Int, max_level::Int; alpha::Real = 0)
    A = -D_matrix(k, max_level) + periodic_LF_matrix(k, max_level; alpha = alpha)
    return A'
end

function periodic_hier_DLF_matrix(k::Int, max_level::Int; alpha::Real = 0)
    Q = hier2pos(k, max_level)
    A = periodic_pos_DLF_matrix(k, max_level; alpha = alpha)
    return *(Q', *(A, Q))
end

function periodic_nodal_DLF_matrix(k::Int, max_level::Int; alpha::Real = 0)
    n2pos = nodal2pos_1D(k, max_level)
    m2n = points2nodal_1D(k, max_level)*nodal2points_1D(k, max_level) #tippfehler modal2points_1D
    pos2m = pos2hier(k, max_level)
    pos2n = m2n * pos2m
    A = periodic_pos_DLF_matrix(k, max_level; alpha = alpha)
    return pos2n * A * n2pos
end

function periodic_point_DLF_matrix(k::Int, max_level::Int; alpha::Real = 0)

    p2pos = nodal2pos_1D(k, max_level) * points2nodal_1D(k, max_level)
    pos2p = modal2points_1D(k, max_level) * pos2hier(k, max_level)
    A = periodic_pos_DLF_matrix(k, max_level; alpha = alpha)
    return pos2p * A * p2pos
end

function periodic_DLF_matrix(k::Int, max_level::Int; alpha::Real = 0, basis = "hier")
    if basis == "hier"
        return periodic_hier_DLF_matrix(k, max_level; alpha=alpha)
    elseif basis == "pos"
        return periodic_pos_DLF_matrix(k, max_level; alpha=alpha)
    elseif basis == "nodal"
        return periodic_nodal_DLF_matrix(k, max_level; alpha=alpha)
    elseif basis == "point"
        return return periodic_point_DLF_matrix(k, max_level; alpha=alpha)
    else
        throw(ArgumentError)
    end
end

# Same function but for anisotropic grids
function periodic_DLF_matrix(k::Int, max_level::Vector{Int}; alpha::Real = 0, basis = "hier")
    max_level = maximum(max_level) # Reicht das? Oder muss für jede Dimension einzeln gerechnet werden?
    if basis == "hier"
        return periodic_hier_DLF_matrix(k, max_level; alpha=alpha)
    elseif basis == "pos"
        return periodic_pos_DLF_matrix(k, max_level; alpha=alpha)
    elseif basis == "nodal"
        return periodic_nodal_DLF_matrix(k, max_level; alpha=alpha)
    elseif basis == "point"
        #return return periodic_point_DLF_matrix(k, max_level; alpha=alpha) #wieso return zweimal?
        return periodic_point_DLF_matrix(k, max_level; alpha=alpha)
    else
        throw(ArgumentError)
    end
end


text:
abstract type Mammal end
habitat(::Mammal) = "planet Earth"

struct Whale <: Mammal
    mass::Float64
    habitat::String
end
Base.show(io::IO, ::Whale) = print(io, "a whale")
habitat(w::Whale) = w.habitat

struct Wolf <: Mammal
    mass::Float64
end
Base.show(io::IO, ::Wolf) = print(io, "a wolf")

arr = [Whale(1000, "ocean"), Wolf(50)]
println("Type of $arr is ", typeof(arr))
for a in arr
    println("Habitat of $a: ", habitat(a))
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using ApproxFun, OrdinaryDiffEq, Sundials
using DiffEqDevTools
using LinearAlgebra
using Plots; gr()

# Define the linear and nonlinear terms
function lin_term(N)
    dx = 1/(N + 1)
    du = ones(N-1) # off diagonal
    du2 = -ones(N-1) # off diagonal
    d = (1/2) * ones(N-2)
    d2 = (-1/2) * ones(N-2)
    diag=-2 * ones(N)
    DiffEqArrayOperator(5e-5*(1/dx^3) * diagm(-2 => d2, -1 => du, 1 => du2, 2 => d))
end

function nl_term(N)
    dx = 1/(N + 1)
    du = ones(N - 1) # super diagonal
    dl = -ones(N - 1) # lower diagonal
    D = (0.2/dx) * diagm(-1 => dl, 1 => du)

    tmp = zeros(N)
    function (du,u,p,t)
        @. tmp = u^2
        mul!(du, D, tmp)
    end
end

# Construct the problem
function kdv(N)
    f1 = lin_term(N)
    f2 = nl_term(N)
    dx = 1 / (N + 1)
    xs = (1:N) * dx

    μ0 = 0.3; σ0 = 0.05
    f0 = x -> 0.9*exp(-(x - μ0)^2 / (2 * σ0^2))
    u0 = f0.(xs)
    prob = SplitODEProblem(f1, f2, u0, (0.0, 1.0))
    xs, prob
end;

xs, prob = kdv(200)
sol = solve(prob, Tsit5(); abstol=1e-11, reltol=1e-11, dt=1e-7, adptive=false)
test_sol = TestSolution(sol);

tslices = [0.0 0.25 0.50 0.75 1.00]
ys = hcat((sol(t) for t in tslices)...)
labels = ["t = $t" for t in tslices]
fn=plot(xs, ys, label=labels)

const LS_Dense = LinSolveFactorize(lu)

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => IMEXEuler(), :dts => 1e-3 * multipliers),
          Dict(:alg => CNAB2(), :dts => 1e-5 * multipliers),
          Dict(:alg => CNLF2(), :dts => 1e-4 * multipliers),
          Dict(:alg => SBDF2(), :dts => 1e-3 * multipliers)]
labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="IMEX methods, dense linsolve, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [#Dict(:alg => IMEXEuler(linsolve=LinSolveGMRES()), :dts => 1e-5 * multipliers),
          Dict(:alg => CNAB2(linsolve=LinSolveGMRES()), :dts => 1e-5 * multipliers),
          Dict(:alg => CNLF2(linsolve=LinSolveGMRES()), :dts => 1e-5 * multipliers),
          Dict(:alg => SBDF2(linsolve=LinSolveGMRES()), :dts => 1e-4 * multipliers)]
labels = ["IMEXEuler" "CNAB2" "CNLF2" "SBDF2"]
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="IMEX methods, Krylov linsolve, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => NorsettEuler(), :dts => 1e-6 * multipliers),
          Dict(:alg => NorsettEuler(krylov=true, m=5), :dts => 1e-3 * multipliers),
          #Dict(:alg => NorsettEuler(krylov=true, m=20), :dts => 1e-3 * multipliers), matrix contains Infs or NaNs
          Dict(:alg => ETDRK2(), :dts => 1e-3 * multipliers),
          Dict(:alg => ETDRK2(krylov=true, m=20), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK2(krylov=true, m=20), :dts => 1e-2 * multipliers)]
labels = hcat("NorsettEuler (caching)", "NorsettEuler (m=5)",# "NorsettEuler (m=20)",
              "ETDRK2 (caching)", "ETDRK2 (m=5)"), "ETDRK2 (m=20)")
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="ExpRK methods, low order")

abstols = 0.1 .^ (5:8) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (1:4)
multipliers = 0.5 .^ (0:3)
setups = [Dict(:alg => CNAB2(), :dts => 1e-5 * multipliers),
          Dict(:alg => CNAB2(linsolve=LinSolveGMRES()), :dts => 1e-5 * multipliers),
          Dict(:alg => ETDRK2(), :dts => 1e-3 * multipliers)]
labels = ["CNAB2 (dense linsolve)" "CNAB2 (Krylov linsolve)" "ETDRK2 (m=5)"]
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="Between family, low orders")

abstols = 0.1 .^ (7:13)
reltols = 0.1 .^ (4:10)
setups = [Dict(:alg => KenCarp3()),
          Dict(:alg => KenCarp4()),
          Dict(:alg => KenCarp5()),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:Dense)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:Dense)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Dense))]
labels = hcat("KenCarp3", "KenCarp4", "KenCarp5", "ARKODE3", "ARKODE4", "ARKODE5")
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="IMEX methods, dense linsolve, medium order")

abstols = 0.1 .^ (7:13)
reltols = 0.1 .^ (4:10)
setups = [Dict(:alg => KenCarp3(linsolve=LinSolveGMRES())),
          Dict(:alg => KenCarp4(linsolve=LinSolveGMRES())),
          Dict(:alg => KenCarp5(linsolve=LinSolveGMRES())),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=3, linear_solver=:GMRES)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=4, linear_solver=:GMRES)),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES))]
labels = ["KenCarp3" "KenCarp4" "KenCarp5" "ARKODE3" "ARKODE4" "ARKODE5"]
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="IMEX methods, medium order")

abstols = 0.1 .^ (7:11) # all fixed dt methods so these don't matter much
reltols = 0.1 .^ (4:8)
multipliers = 0.5 .^ (0:4)
setups = [Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
          #Dict(:alg => ETDRK3(krylov=true, m=5), :dts => 1e-2 * multipliers),matrix contains Infs or NaNs
          Dict(:alg => ETDRK4(), :dts => 1e-3 * multipliers),
          #Dict(:alg => ETDRK4(krylov=true, m=5), :dts => 1e-2 * multipliers),matrix contains Infs or NaNs
          Dict(:alg => HochOst4(), :dts => 1e-2 * multipliers)]
          #Dict(:alg => HochOst4(krylov=true, m=5), :dts => 1e-2 * multipliers)] matrix contains Infs or NaNs
labels = hcat("ETDRK3 (caching)", "ETDRK4 (caching)",# "ETDRK3 (m=5)", "ETDRK4 (m=5)"
              "HochOst4 (caching)")#, "HochOst4 (m=5)")
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="ExpRK methods, medium order")

abstols = 0.1 .^ (7:11)
reltols = 0.1 .^ (4:8)
multipliers = 0.5 .^ (0:4)
setups = [Dict(:alg => KenCarp5()),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:Dense)),
          Dict(:alg => KenCarp5(linsolve=LinSolveGMRES())),
          Dict(:alg => ARKODE(Sundials.Implicit(), order=5, linear_solver=:GMRES)),
          Dict(:alg => ETDRK3(), :dts => 1e-2 * multipliers),
          Dict(:alg => ETDRK4(), :dts => 1e-3 * multipliers)]
labels = hcat("KenCarp5 (dense linsolve)", "ARKODE (dense linsolve)", "KenCarp5 (Krylov linsolve)",
              "ARKODE (Krylov linsolve)", "ETDRK3 (m=5)", "ETDRK4 (m=5)")
@time wp = WorkPrecisionSet(prob,abstols,reltols,setups;
                            print_names=true, names=labels,
                            numruns=5, error_estimate=:l2,
                            save_everystep=false, appxsol=test_sol, maxiters=Int(1e5));

plot(wp, label=labels, markershape=:auto, title="Between family, medium order")

using DiffEqBenchmarks
DiffEqBenchmarks.bench_footer(WEAVE_ARGS[:folder],WEAVE_ARGS[:file])


text:
module test_triangulate2
using Triangulate

function test()
    #     *
    #     |
    #  *--*--*

    triin = Triangulate.TriangulateIO(
        pointlist = Float64[0.0 1.0 2.0 1.0
                            0.0 0.0 0.0 1.0],
        segmentlist = Int32[1 2 2
                            2 3 4],
        # don't give segmentmarkerlist !
    )

    (triout, vorout) = triangulate("pcDQ", triin)

    return all([
        size(triout.pointlist, 2) == 4,
        size(triout.trianglelist, 2) == 2,
        size(triout.segmentlist, 2) == 5,
    ])
end

end