text:
#Structs used for Forest-Photosynthesis-Hydraulic model

#Constants
mutable struct Constants{T<:Float64}
    M_H2O::T #Molar mass water (Kg mol⁻¹)
    M_C::T #Molar mass Carbon (Kg mol⁻¹
    ρ_H2O::T #Density water (Kg m⁻³)
    g::T #gravitational acceleration (m s⁻²)
    r::T #ratio water conductance:co2 conductance (-)
end
function Constants(;M_H2O::T=0.018,M_C::T=0.012,ρ_H2O::T=997.0,g::T=9.82,r::T=1.6) where {T<:Float64}
    Constants(M_H2O,M_C,ρ_H2O,g,r)
end

#Struct containing environmental variables
mutable struct EnvironmentStruct{T<:Float64}
    I₀::T #Above canopy irradiance (mol m⁻² s⁻¹)
    Cₐ::T #Ambient carbon dioxide partial pressure (Pa)
    P::T #Atmospheric pressure (Pa)
    Tₐ::T #Ambient air tempreture (°C)
    VPD::T #Vapour-pressure deficit
    θₛ::T #Soil volumetric water content (-)
end
#Standard vaules
EnvironmentStruct(;I₀::T=820.0*10^-6,Cₐ::T=400.0/10.0,P::T=1.0*10^5,
Tₐ::T=25.0,VPD::T=800.0,θₛ::T = 0.15) where {T<:Float64} =  EnvironmentStruct(I₀,Cₐ,P,Tₐ,VPD,θₛ)

mutable struct ArrheniusKineticRate{T<:Float64}
    K_ref::T #Kinetic rate at referance temperature
    T_ref::T #Reference temperature (°C)
    Eₐ::T #Activation energy (J mol⁻¹)
end
#Calcualte Kinetic rate at temperature Temp (°C)
function (K::ArrheniusKineticRate)(T::Float64)
    T += 273.15 #Temperature in Kelvin
    T_ref = K.T_ref+273.15 #Reference temperature in Kelvin
    R = 8.3145 #Gas constant (J K⁻¹ mol⁻¹)
    return K.K_ref*exp(K.Eₐ*(1/T_ref-1/T)/R)
end

mutable struct EnzymeKineticRate{T<:Float64}
    Kₒₚₜ::T #Kinetic rate at optimal temperature
    Tₒₚₜ::T #Optimal temperature (°C)
    E_a::T #Activation energy (J mol⁻¹)
    E_d::T #Deactivation energy (J mol⁻¹)
end
#Calcualte Kinetic rate at temperature Temp (°C)
function (K::EnzymeKineticRate)(T::Float64)
    T += 273.15 #Temperature in Kelvin
    Tₒₚₜ = K.Tₒₚₜ+273.15 #Optimal temperature in Kelvin
    R = 8.3145 #Gas constant (J K⁻¹ mol⁻¹)
    return K.Kₒₚₜ*K.E_d*exp(K.E_a*(T-Tₒₚₜ)/(T*R*Tₒₚₜ))/(K.E_d-K.E_a*(1-exp(K.E_d*(T-Tₒₚₜ)/(T*R*Tₒₚₜ))))
end

mutable struct PhotoKineticRates{T<:ArrheniusKineticRate}
    K_C::T #Carboxylation Michaelis-Menten constant
    K_O::T #Oxygenation Michaelis-Menten constant
    Γ::T #The CO2 compensation point in the absence of mitochondrial (day) respiration (Pa)
    b_r::T #The rate of dark respiration to Vmax ratio (Pa)
    Jmax::EnzymeKineticRate #
    O::Float64 #Intercellular partial pressure of oxygen (Pa)
end
function PhotoKineticRates(;Temp_ref::T=25.0, #Reference temperature (°C)
    K_C_ref::T = 40.40, K_C_E_a::T = 59.36e3, #Reference value (Pa) and activation energy (J mol⁻¹) for carboxylation Michaelis-Menten constant
    K_O_ref::T = 248.0e2, K_O_E_a::T = 35.94e3, #Reference value (Pa) and activation energy (J mol⁻¹) for oxygenation Michaelis-Menten constant
    Γ_ref::T = 4.17, Γ_E_a::T = 23.42e3, #Reference value (Pa) and activation energy (J mol⁻¹) for CO2 compensation point
    b_r_ref::T = 0.02, b_r_E_a::T = 7.88e3, #Reference value (-) and activation energy (J mol⁻¹) for dark respiration to Vmax ratio
    Jmaxₒₚₜ::T = 184.0e-6, Jmax_Tempₒₚₜ::T = 31.85, Jmax_E_a::T = 47.4e3, Jmax_E_d::T = 200.0e3,
    #Optimal value (mol m⁻² s⁻¹), optimal temperature (°C), activation (J mol⁻¹) and deactivaiton (J mol⁻¹) energy for Jamx
    O::T = 20.5e3) where {T<:Float64}

    K_C = ArrheniusKineticRate(K_C_ref,Temp_ref,K_C_E_a)
    K_O = ArrheniusKineticRate(K_O_ref,Temp_ref,K_O_E_a)
    Γ = ArrheniusKineticRate(Γ_ref,Temp_ref,Γ_E_a)
    b_r = ArrheniusKineticRate(b_r_ref,Temp_ref,b_r_E_a)
    Jmax = EnzymeKineticRate(Jmaxₒₚₜ,Jmax_Tempₒₚₜ,Jmax_E_a,Jmax_E_d)

    PhotoKineticRates(K_C,K_O,Γ,b_r,Jmax,O)
end

#Parameters values used in the Farquhar photosynthesis model
mutable struct PhotoPar{T<:Float64}
    K::T #Combined Michaelis-Menten constant (Pa)
    b_r::T #Ratio between rate of dark respiration and Vmax (Pa)
    Γ::T #The CO2 compensation point in the absence of mitochondrial (day) respiration (Pa)
    α::T #quantum yield of electron tansport
    θ::T #Curvature of the light response curve
    b_Jmax::T #The ratio between Jmax at tempreatur Temp and Jmax at optimal tempreture (Jmaxₒₚₜ)
end
#Reference values at temperature = 25°C
PhotoPar(;K::T=73.8,b_r::T=0.02,Γ::T=4.17,α::T=0.36,θ::T=0.7,b_Jmax::T=0.812) where {T<:Float64} =
PhotoPar(K,b_r,Γ,α,θ,b_Jmax)
#Function for initiate PhotoPar object at a tempreture of Temp °C
function PhotoPar(kinetic::PhotoKineticRates,Temp::T;α::T=0.36,θ::T=0.7) where {T<:Float64}
    K = kinetic.K_C(Temp)*(1+kinetic.O/kinetic.K_O(Temp))
    b_r = kinetic.b_r(Temp)
    Γ = kinetic.Γ(Temp)
    b_Jmax = kinetic.Jmax(Temp)/kinetic.Jmax.Kₒₚₜ

    PhotoPar(K,b_r,Γ,α,θ,b_Jmax)
end
#Calc PhotoPar values at temperature Temp (⁰C)
function PhotoPar!(photo::PhotoPar,kinetic::PhotoKineticRates,Temp::Float64)
    K = kinetic.K_C(Temp)*(1+kinetic.O/kinetic.K_O(Temp))
    b_r = kinetic.b_r(Temp)
    Γ = kinetic.Γ(Temp)
    b_Jmax = kinetic.Jmax(Temp)/kinetic.Jmax.Kₒₚₜ

    photo.K = K
    photo.b_r = b_r
    photo.Γ = Γ
    photo.b_Jmax = b_Jmax

    return nothing
end

#Struct containing variables which defines the size of a mean tree in the stand
mutable struct TreeSize{T<:Float64}
    Wf::T #Foliage weight (kg)
    Ww::T #Foliage try weight (kg)
    H::T #Tree height (m)
    Hs::T #Height to crown base (m)
    As::T #sapwood area at crown base (m²)
    B::T #Basal area (single tree) (m²)
    N::T #Stem density (# trees m⁻² ground area)
end

#Struct containing tree and stand parameters
mutable struct TreePar{T<:Float64}
    αf::T #Foliage weight to Sapwood area ratio (Kg m⁻²)
    ρw::T #Sapwood density (Kg m³)
    β₁::T #Parameter for estimating the average length of an active pipe (β₁H+β₂Hs)
    β₂::T #Parameter for estimating the average length of an active pipe (β₁H+β₂Hs)
    Tf::T #Average longevity of foliage (time⁻¹)
    Tr::T #Average longevity of fine root (time⁻¹)
    y::T #Conversion efficiency of C to biomass (includes growth respiration)
    z::T #Scaling exponent As∝(H-Hs)ᶻ
    Nₛ::T #Maximum N uptake per fine root mass (Kg N Kg⁻¹ Wr)
    rₘ::T #Specific Nitrogen maintenance respiration rate (Kg⁻¹ N year⁻¹)
    k::T #Light extinction coefficient (-)
    m::T #Average leaf transmittance (-)
    a_Jmax::T #Slope of the Nitrogen per leaf area (Nₐ)-Jmaxₒₚₜ line (mol m⁻² leaf s⁻¹ Nₐ⁻¹)
    b_Jmax::T #Intercept of the Nitrogen per leaf area (Nₐ)-Jmaxₒₚₜ line (mol m⁻² leaf s⁻¹)
    LMA::T #Leaf mass area (Kg foliage m⁻² leaf)
    Kr::T #Maximum stand fine root mass wich result in 1/2 of maximum N uptake (Kg C m⁻²)
    rW::T #Sapwood to foliage nitrogen concentration (mass) (-)
    rR::T #Fine root to foliage nitrogen concentration (mass) (-)
    r_gₛ::T #total leaf conductance (stomatal+mesophyll) to stomatal conductance (gₜ/gₛ)
    Xₜ::T #Factor [0,1] accounting for the delayed effect of temperature on gross primary production (-)
end
#Standard values
TreePar(;αf::T=460.0,ρw::T=400.0,β₁::T=1.27,β₂::T=-0.27,Tf::T=3.33,Tr::T=1.25,y::T=1.54,z::T=1.86,
Nₛ::T=0.04,rₘ::T=24.0,k::T=0.52,m::T=0.05,a_Jmax::T=0.033,b_Jmax::T=-1.1e-5,LMA::T=0.256
,Kr::T = 0.35,rW::T = 0.07,rR::T = 0.6,r_gₛ::T = 0.42,Xₜ::T = 1.0) where {T<:Float64}  =
TreePar(αf,ρw,β₁,β₂,Tf,Tr,y,z,Nₛ,rₘ,k,m,a_Jmax,b_Jmax,LMA,Kr,rW,rR,r_gₛ,Xₜ)

#Parameteres used in the Hydraulics model
mutable struct HydraulicsPar{T<:Float64}
    ψ₅₀::T #water potential which causes 50% loss in leaf conductance (MPa)
    b::T #Sensitivity of leaf conductance to water potential (-)
    i::T #Sensitivity to hydraulic failure
    Kₓₗ₀::T #Maximum xylem-leaf hydraulic conductance (mol m⁻² leaf s⁻¹ MPa⁻¹)
    Kₛᵣ₀::T #Maximum soil hydraulic conductance
end
#Standard values
HydraulicsPar(;ψ₅₀::T=-2.0,b::T=2.0,i::T=1.0,Kₓₗ₀::T=0.01,Kₛᵣ₀::T=0.1) where {T<:Float64} =
HydraulicsPar(ψ₅₀,b,i,Kₓₗ₀,Kₛᵣ₀)

#Collection of structs used for the Photosynthesis and Hydraulic model
mutable struct CCPHStruct
    cons::Constants
    env::EnvironmentStruct
    treepar::TreePar
    treesize::TreeSize
    photopar::PhotoPar
    hydPar::HydraulicsPar
end

#Struct collecting output from the Coupled Canopy Photosynthesis and Hydraulic model
mutable struct CCPHOutput
    P::Float64
    αr::Float64
    ψ_c::Float64
    Kₓₗ::Float64
    K_cost::Float64
    Gain::Float64
end

#Struct collecting times series results from growth simulation
mutable struct CCPHTS{T<:Float64}
    Wf::Array{T,1}
    Ww::Array{T,1}
    H::Array{T,1}
    Hs::Array{T,1}
    As::Array{T,1}
    B::Array{T,1}
    N::Array{T,1}
    P::Array{T,1}
    αr::Array{T,1}
    ψ_c::Array{T,1}
    Kₓₗ::Array{T,1}
    K_cost::Array{T,1}
    Gain::Array{T,1}
    gₛ::Array{T,1}
    Nₘ_f::Array{T,1}
    gₛ_crit::Array{T,1}
end
CCPHTS() = CCPHTS(Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],
Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],Float64[],Float64[])
function CCPHTS!(ccphts::CCPHTS,treesize::TreeSize)
    push!(ccphts.Wf,treesize.Wf)
    push!(ccphts.Ww,treesize.Ww)
    push!(ccphts.H,treesize.H)
    push!(ccphts.Hs,treesize.Hs)
    push!(ccphts.As,treesize.As)
    push!(ccphts.B,treesize.B)
    push!(ccphts.N,treesize.N)
end
function CCPHTS!(ccphts::CCPHTS,modeloutput::CCPHOutput,gₛ::T,Nₘ_f::T,gₛ_crit::T) where {T<:Float64}
    push!(ccphts.P,modeloutput.P)
    push!(ccphts.αr,modeloutput.αr)
    push!(ccphts.ψ_c,modeloutput.ψ_c)
    push!(ccphts.Kₓₗ,modeloutput.Kₓₗ)
    push!(ccphts.K_cost,modeloutput.K_cost)
    push!(ccphts.Gain,modeloutput.Gain)
    push!(ccphts.gₛ,gₛ)
    push!(ccphts.Nₘ_f,Nₘ_f)
    push!(ccphts.gₛ_crit,gₛ_crit)
end

#Containter for weather related time series
mutable struct WeatherTS{T<:Float64}
    date::Array{Dates.DateTime,1} #Dates
    daylight::Array{T,1} #Daylight lenght, for each timestep, given in seconds
    tot_annual_daylight::Array{T,1} #Total annual daylight lenght during growth period, given in seconds
    PAR::Array{T,1} #mol s⁻¹ m⁻²
    temp::Array{T,1} #Mean air temperature (⁰C)
    VPD::Array{T,1} #Pa
    θₛ::Array{T,1} #Volumetric soil water content
    acclimation_fac::Array{T,1} #Factor [0,1] accounting for annual acclimation cycle, see Mäkelä 2004 & 2008
end


text:
using Plots
using Random

# Number of observations
T = 50

# Variances
Q = 0.1
R = 0.5

# Define model with prior N(0, 1000)
model = LocalLevel(zeros(T), Q, R)

# Simulate
model.y, x = simulate(model)

# Filter
kf = kalmanfilter(model)

# Smoother
ks = kalmansmoother(model, kf.priors, kf.posteriors)

# Plot
scatter(y, label = "Measurements", legend = :topright, title = "Univariate local level model",
        markerstrokewidth = 0, markercolor = "black", foreground_color_legend = nothing)
plot!(kf.posteriors.μ, label = "Filter", linewidth = 2, color = "red")
plot!(ks.μ, label = "Smoother", linewidth = 2, color = "green")

ext: jl


text:
import OrdinaryDiffEq
import DiffEqBase: DynamicalODEProblem
import DiffEqSensitivity:
    solve,
    ODEProblem,
    ODEAdjointProblem,
    InterpolatingAdjoint,
    ZygoteVJP,
    ReverseDiffVJP
import RecursiveArrayTools: ArrayPartition

sol = solve(
    DynamicalODEProblem(
        (v, x, p, t) -> [0.0, 0.0],

        # ERROR: LoadError: type Nothing has no field x
        # (v, x, p, t) -> [0.0, 0.0],

        # ERROR: LoadError: MethodError: no method matching ndims(::Type{Nothing})
        (v, x, p, t) -> v,

        [0.0, 0.0],
        [0.0, 0.0],
        (0.0, 1.0),
    ),OrdinaryDiffEq.Tsit5()
)

solve(
    ODEAdjointProblem(
        sol,
        InterpolatingAdjoint(),
        (out, x, p, t, i) -> (out .= 0),
        [sol.t[end]],
    ),OrdinaryDiffEq.Tsit5()
)

dyn_v(v_ap, x_ap, p, t) = ArrayPartition(zeros(), [0.0])
# Originally, I imagined that this may be a bug in Zygote, and it still may be, but I tried doing a pullback on this
# function on its own and didn't have any trouble with that. So I'm led to believe that it has something to do with
# how DiffEqSensitivity is invoking Zygote. At least this was as far as I was able to simplify the reproduction.
dyn_x(v_ap, x_ap, p, t) = begin
    # ERROR: LoadError: MethodError: no method matching ndims(::Type{NamedTuple{(:x,),Tuple{Tuple{Nothing,Array{Float64,1}}}}})
    v = v_ap.x[2]

    # ERROR: LoadError: type Nothing has no field x
    # v = [0.0]
    ArrayPartition(zeros(), v)
end

v0 = [-1.0]
x0 = [0.75]

sol = solve(
    DynamicalODEProblem(
        dyn_v,
        dyn_x,
        ArrayPartition(zeros(), v0),
        ArrayPartition(zeros(), x0),
        (0.0, 1.0),
        zeros()
    ),OrdinaryDiffEq.Tsit5(),
    # Without setting parameters, we end up with https://github.com/SciML/DifferentialEquations.jl/issues/679 again.
    p = zeros()
)

g = ArrayPartition(ArrayPartition(zeros(), zero(v0)), ArrayPartition(zeros(), zero(x0)))
bwd_sol = solve(
    ODEAdjointProblem(
        sol,
        InterpolatingAdjoint(autojacvec=ZygoteVJP()),
        # Also fails, but due to a different bug:
        # InterpolatingAdjoint(autojacvec=ReverseDiffVJP()),
        (out, x, p, t, i) -> (out[:] = g),
        [sol.t[end]],
    ),OrdinaryDiffEq.Tsit5()
)


text:
using ITensors

#
# DMRG calculation of the extended Hubbard model
# ground state wavefunction, and spin densities
#

# Usage:
#
# Run:
#
# $ julia exthubbard.jl input.jl
#
# or with options, like
#
# $ julia exthubbard.jl input.jl N=100 Npart=20 nsweep=4 cutoff=1e-11 maxdim=20,40,80,100 noise=1e-12,1e-13
#
# In the REPL (if you just type `julia`), you can modify the command line options as follows (or just modify "input.jl" and reinclude it):
#
# julia> push!(ARGS, "input.jl");
#
# julia> push!(ARGS, "N=100", "Npart=20");
#
# julia> include("exthubbard.jl");
#

# Parse arguments to overide defaults
args = argsdict(; prefix="file")

# Include the file from the first positional
# argument. If none is specified, use
# input.jl by default.
include(get(args, "file1", "input.jl"))

# Create the sweeps object from the input
# file
sweeps = Sweeps(nsw, sweeps_args)
# Extract the original (default) values
maxdim = get_maxdims(sweeps)
mindim = get_mindims(sweeps)
cutoff = get_cutoffs(sweeps)
noise = get_noises(sweeps)

# Extract the argument values from
# the dictionary

N = get(args, "N", N)
Npart = get(args, "Npart", Npart)
t1 = get(args, "t1", t1)
t2 = get(args, "t2", t2)
U = get(args, "U", U)
V1 = get(args, "V1", V1)
nsw = get(args, "nsweep", nsw)
maxdim = get(args, "maxdim", maxdim)
mindim = get(args, "mindim", mindim)
cutoff = get(args, "cutoff", cutoff)
noise = get(args, "noise", noise)

#
# Alternatively, define all of the
# variables for all of the inputs
# programatically. You can use the
# keyword argument `as_symbols = true`
# in the `parse_args` function to
# automatically turn the keys of
# the dictionary into symbols.
#

#for (arg, val) in args
#  arg == "input_file" && continue
#  @eval $(Symbol(arg)) = $val
#end

sweeps = Sweeps(nsw)
maxdim!(sweeps, maxdim...)
mindim!(sweeps, mindim...)
cutoff!(sweeps, cutoff...)
noise!(sweeps, noise...)

@show sweeps
@show N, Npart

sites = siteinds("Electron", N; conserve_qns=true)

ampo = OpSum()
for i in 1:N
  ampo .+= U, "Nupdn", i
end
for b in 1:(N - 1)
  ampo .+= -t1, "Cdagup", b, "Cup", b + 1
  ampo .+= -t1, "Cdagup", b + 1, "Cup", b
  ampo .+= -t1, "Cdagdn", b, "Cdn", b + 1
  ampo .+= -t1, "Cdagdn", b + 1, "Cdn", b
  ampo .+= V1, "Ntot", b, "Ntot", b + 1
end
for b in 1:(N - 2)
  ampo .+= -t2, "Cdagup", b, "Cup", b + 2
  ampo .+= -t2, "Cdagup", b + 2, "Cup", b
  ampo .+= -t2, "Cdagdn", b, "Cdn", b + 2
  ampo .+= -t2, "Cdagdn", b + 2, "Cdn", b
end
H = MPO(ampo, sites)

state = ["Emp" for n in 1:N]
p = Ref(Npart)
for i in N:-1:1
  if p[] > i
    println("Doubly occupying site $i")
    state[i] = "UpDn"
    p[] -= 2
  elseif p[] > 0
    println("Singly occupying site $i")
    state[i] = (isodd(i) ? "Up" : "Dn")
    p[] -= 1
  end
end
# Initialize wavefunction to be bond
# dimension 10 random MPS with number
# of particles the same as `state`
psi0 = randomMPS(sites, state, 10)

# Check total number of particles:
@show flux(psi0)

# Start DMRG calculation:
energy, psi = dmrg(H, psi0, sweeps)

upd = fill(0.0, N)
dnd = fill(0.0, N)
for j in 1:N
  orthogonalize!(psi, j)
  psidag_j = dag(prime(psi[j], "Site"))
  upd[j] = scalar(psidag_j * op(sites, "Nup", j) * psi[j])
  dnd[j] = scalar(psidag_j * op(sites, "Ndn", j) * psi[j])
end

println("Up Density:")
for j in 1:N
  println("$j $(upd[j])")
end
println()

println("Dn Density:")
for j in 1:N
  println("$j $(dnd[j])")
end
println()

println("Total Density:")
for j in 1:N
  println("$j $(upd[j]+dnd[j])")
end
println()

println("\nGround State Energy = $energy")


text:
"""
    AssemblyModule

Module for assemblers  of system matrices and vectors.
"""
module AssemblyModule

__precompile__(true)

using ..FTypesModule: FInt, FFlt, FCplxFlt, FFltVec, FIntVec, FFltMat, FIntMat, FMat, FVec, FDataDict
using  SparseArrays: sparse, spzeros
using  LinearAlgebra: diag

"""
    AbstractSysmatAssembler

Abstract type of system-matrix assembler.
"""
abstract type AbstractSysmatAssembler end;

"""
    SysmatAssemblerSparse{T<:Number} <: AbstractSysmatAssembler

Type for assembling a sparse global matrix from elementwise matrices.

!!! note

    All fields of the datatype are private. No need to access them directly.
"""
mutable struct SysmatAssemblerSparse{T<:Number} <: AbstractSysmatAssembler
    buffer_length::FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer::FInt;
    ndofs_row::FInt; ndofs_col::FInt;
    nomatrixresult::Bool
end

"""
    SysmatAssemblerSparse(zero::T=0.0, nomatrixresult = false) where {T<:Number}

Construct blank system matrix assembler.

The matrix entries are of type `T`. The assembler either produces a sparse
matrix (when `nomatrixresult = true`), or does not (when `nomatrixresult =
false`). When the assembler does not produce the sparse matrix when
`makematrix!` is called, it still can be constructed from the buffers stored in
the assembler.

# Example

This is how a sparse matrix is assembled from two rectangular dense matrices.
```
    a = SysmatAssemblerSparse(0.0)
    startassembly!(a, 5, 5, 3, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
        0.786024  0.00206713  0.995379  0.780298
        0.845816  0.198459    0.355149  0.224996]
    assemble!(a, m, [1 7 5], [5 2 1 4])
    m = [0.146618  0.53471   0.614342    0.737833
         0.479719  0.41354   0.00760941  0.836455
         0.254868  0.476189  0.460794    0.00919633
         0.159064  0.261821  0.317078    0.77646
         0.643538  0.429817  0.59788     0.958909]
    assemble!(a, m, [2 3 1 7 5], [6 7 3 4])
    A = makematrix!(a)
```

When the `nomatrixresult` is set as true, no matrix is produced.
```
    a = SysmatAssemblerSparse(0.0, true)
    startassembly!(a, 5, 5, 3, 7, 7)
    m = [0.24406   0.599773    0.833404  0.0420141
        0.786024  0.00206713  0.995379  0.780298
        0.845816  0.198459    0.355149  0.224996]
    assemble!(a, m, [1 7 5], [5 2 1 4])
    m = [0.146618  0.53471   0.614342    0.737833
         0.479719  0.41354   0.00760941  0.836455
         0.254868  0.476189  0.460794    0.00919633
         0.159064  0.261821  0.317078    0.77646
         0.643538  0.429817  0.59788     0.958909]
    assemble!(a, m, [2 3 1 7 5], [6 7 3 4])
    A = makematrix!(a)
```
Here `A` is a sparse zero matrix. To construct the correct matrix is still
possible, for instance like this:
```
    a.nomatrixresult = false
    A = makematrix!(a)
```
At this point all the buffers of the assembler have been cleared, and
`makematrix!(a) ` is no longer possible.

"""
function SysmatAssemblerSparse(zero::T=0.0, nomatrixresult = false) where {T<:Number}
    return SysmatAssemblerSparse{T}(0,[zero],[0],[0],0,0,0,nomatrixresult)
end

"""
    startassembly!(self::SysmatAssemblerSparse{T},
      elem_mat_nrows::FInt, elem_mat_ncols::FInt, elem_mat_nmatrices::FInt,
      ndofs_row::FInt, ndofs_col::FInt) where {T<:Number}

Start the assembly of a global matrix.

The method makes buffers for matrix assembly. It must be called before
the first call to the method `assemble!`.
- `elem_mat_nrows`= number of rows in typical element matrix,
- `elem_mat_ncols`= number of columns in a typical element matrix,
- `elem_mat_nmatrices`= number of element matrices,
- `ndofs_row`= Total number of equations in the row direction,
- `ndofs_col`= Total number of equations in the column direction.
"""
function startassembly!(self::SysmatAssemblerSparse{T}, elem_mat_nrows::FInt, elem_mat_ncols::FInt, elem_mat_nmatrices::FInt, ndofs_row::FInt, ndofs_col::FInt) where {T<:Number}
    self.buffer_length = elem_mat_nmatrices*elem_mat_nrows*elem_mat_ncols;
    self.rowbuffer = zeros(FInt,self.buffer_length);
    self.colbuffer = zeros(FInt,self.buffer_length);
    self.matbuffer = zeros(T,self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs_row = ndofs_row;
    self.ndofs_col = ndofs_col;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparse{T}, mat::FMat{T},
      dofnums_row::FIntMat, dofnums_col::FIntMat) where {T<:Number}

Assemble a rectangular matrix.
"""
function assemble!(self::SysmatAssemblerSparse{T}, mat::FMat{T}, dofnums_row::FIntVec, dofnums_col::FIntVec) where {T<:Number}
    # Assembly of a rectangular matrix.
    # The method assembles a rectangular matrix using the two vectors of
    # equation numbers for the rows and columns.
    nrows=length(dofnums_row); ncolumns=length(dofnums_col);
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    @inbounds for j=1:ncolumns
        @inbounds for i=1:nrows
            self.matbuffer[p] = mat[i,j] # serialized matrix
            self.rowbuffer[p] = dofnums_row[i];
            self.colbuffer[p] = dofnums_col[j];
            p=p+1
        end
    end
    self.buffer_pointer=p;
    return self
end

function assemble!(self::SysmatAssemblerSparse{T}, mat::FMat{T}, dofnums_row::FIntMat, dofnums_col::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums_row), vec(dofnums_col))
end

"""
    makematrix!(self::SysmatAssemblerSparse)

Make a sparse matrix.
"""
function makematrix!(self::SysmatAssemblerSparse)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    # Here we will go through the rows and columns, and whenever the row or
    # the column refer to indexes outside of the limits of the matrix, the
    # corresponding value will be set to 0 and assembled to row and column 1.
    @inbounds for j=1:self.buffer_pointer-1
        if (self.rowbuffer[j] > self.ndofs_row) || (self.rowbuffer[j] <= 0)
            self.rowbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
        if (self.colbuffer[j] > self.ndofs_col) || (self.colbuffer[j] <= 0)
            self.colbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
    end
    if self.nomatrixresult
        # No actual sparse matrix is returned. The entire result of the assembly
        # is preserved in the assembler buffers.
        return spzeros(self.ndofs_row, self.ndofs_col);
    else
        # The sparse matrix is constructed and returned. The  buffers used for
        # the assembly are cleared.
        S = sparse(self.rowbuffer[1:self.buffer_pointer-1],
                     self.colbuffer[1:self.buffer_pointer-1],
                     self.matbuffer[1:self.buffer_pointer-1],
                     self.ndofs_row, self.ndofs_col);
        self = SysmatAssemblerSparse(zero(eltype(self.matbuffer)))# get rid of the buffers
        return S
    end
end

"""
    SysmatAssemblerSparseSymm{T<:Number} <: AbstractSysmatAssembler

Assembler for a **symmetric square** matrix  assembled from symmetric square
matrices.
"""
mutable struct SysmatAssemblerSparseSymm{T<:Number} <: AbstractSysmatAssembler
    # Type for assembling of a sparse global matrix from elementwise matrices.
    buffer_length:: FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer:: FInt;
    ndofs:: FInt;
end

"""
    SysmatAssemblerSparseSymm(zero::T=0.0) where {T<:Number}

Construct blank system matrix assembler for symmetric matrices. The matrix
entries are of type `T`.

# Example

This is how a symmetric sparse matrix is assembled from two square dense matrices.
```
	a = SysmatAssemblerSparseSymm(0.0)
	startassembly!(a, 5, 5, 3, 7, 7)
	m = [0.24406   0.599773    0.833404  0.0420141
		0.786024  0.00206713  0.995379  0.780298
		0.845816  0.198459    0.355149  0.224996]
	assemble!(a, m'*m, [5 2 1 4], [5 2 1 4])
	m = [0.146618  0.53471   0.614342    0.737833
		 0.479719  0.41354   0.00760941  0.836455
		 0.254868  0.476189  0.460794    0.00919633
		 0.159064  0.261821  0.317078    0.77646
		 0.643538  0.429817  0.59788     0.958909]
	assemble!(a, m'*m, [2 3 1 5], [2 3 1 5])
	A = makematrix!(a)
```
# See also
[]
"""
function SysmatAssemblerSparseSymm(zero::T=0.0) where {T<:Number}
    return SysmatAssemblerSparseSymm{T}(0,[zero],[0],[0],0,0)
end

"""
    startassembly!(self::SysmatAssemblerSparseSymm{T},
      elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt,
      ndofs::FInt, ignore2::FInt) where {T<:Number}

Start the assembly of a symmetric square global matrix.

The method makes buffers for matrix assembly. It must be called before
the first call to the method `assemble!`.
- `elem_mat_nrows`= number of rows in typical element matrix,
- `elem_mat_ncols`= number of columns in a typical element matrix,
- `elem_mat_nmatrices`= number of element matrices,
- `ndofs_row`= Total number of equations in the row direction,
- `ndofs_col`= Total number of equations in the column direction.
"""
function startassembly!(self::SysmatAssemblerSparseSymm{T}, elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt, ndofs::FInt, ignore2::FInt) where {T<:Number}
    self.buffer_length = elem_mat_nmatrices*elem_mat_dim^2;
    self.rowbuffer = zeros(FInt,self.buffer_length);
    self.colbuffer = zeros(FInt,self.buffer_length);
    self.matbuffer = zeros(T,self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs = ndofs;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}

Assemble a square symmetric matrix.

`dofnums` are the row degree of freedom numbers, the column degree of freedom
number input is ignored (the row and column numbers are the same).
"""
function assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}
    # Assembly of a square symmetric matrix.
    # The method assembles the lower triangle of the square symmetric matrix using the two vectors of
    # equation numbers for the rows and columns.
    nrows=length(dofnums); ncolumns=nrows;
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    @inbounds for j=1:ncolumns
        @inbounds for i=j:nrows
            self.matbuffer[p] = mat[i,j] # serialized matrix
            self.rowbuffer[p] = dofnums[i];
            self.colbuffer[p] = dofnums[j];
            p=p+1
        end
    end
    self.buffer_pointer=p;
    return self
end

"""
	assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T},
	  dofnums::FIntMat, ignore::FIntMat) where {T<:Number}

Assemble a square symmetric matrix.

`dofnums` are the row degree of freedom numbers, the column degree of freedom
number input is ignored (the row and column numbers are the same).
"""
function assemble!(self::SysmatAssemblerSparseSymm{T}, mat::FMat{T}, dofnums::FIntMat, ignore::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums), vec(ignore))
end

"""
    makematrix!(self::SysmatAssemblerSparseSymm)

Make a sparse symmetric square matrix.
"""
function makematrix!(self::SysmatAssemblerSparseSymm)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    # Here we will go through the rows and columns, and whenever the row or
    # the column refer to indexes outside of the limits of the matrix, the
    # corresponding value will be set to 0 and assembled to row and column 1.
    @inbounds for j=1:self.buffer_pointer-1
        if (self.rowbuffer[j] > self.ndofs) || (self.rowbuffer[j] <= 0)
            self.rowbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
        if (self.colbuffer[j] > self.ndofs) || (self.colbuffer[j] <= 0)
            self.colbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
    end
    S = sparse(self.rowbuffer[1:self.buffer_pointer-1],
               self.colbuffer[1:self.buffer_pointer-1],
               self.matbuffer[1:self.buffer_pointer-1],
               self.ndofs, self.ndofs);
    #  Now we need to construct the other triangle of the matrix. The diagonal
    #  will be duplicated.
    S = S + transpose(S);
    @inbounds for j=1:size(S,1)
        S[j,j]=S[j,j]/2.0;      # the diagonal is there twice; fix it;
    end
    self = SysmatAssemblerSparseSymm(zero(eltype(self.matbuffer))) # get rid of the buffers
    return S
end

"""
    SysmatAssemblerReduced{T<:Number} <: AbstractSysmatAssembler

Type for assembling a sparse global matrix from elementwise matrices.

!!! note

    All fields of the datatype are private. No need to access them directly.
"""
mutable struct SysmatAssemblerReduced{T<:Number} <: AbstractSysmatAssembler
    m::Matrix{T} # reduced system matrix
    ndofs_row::FInt; ndofs_col::FInt;
    red_ndofs_row::FInt; red_ndofs_col::FInt;
    t::Matrix{T} # transformation matrix
end

function SysmatAssemblerReduced(t::Matrix{T}) where {T<:Number}
    ndofs_row = ndofs_col = size(t, 1)
    red_ndofs_row = red_ndofs_col = size(t, 2)
    m = fill(zero(T), red_ndofs_row, red_ndofs_col)
    return SysmatAssemblerReduced{T}(m, ndofs_row, ndofs_col, red_ndofs_row, red_ndofs_col, t)
end

function startassembly!(self::SysmatAssemblerReduced{T}, elem_mat_nrows::FInt, elem_mat_ncols::FInt, elem_mat_nmatrices::FInt, ndofs_row::FInt, ndofs_col::FInt) where {T<:Number}
    @assert self.ndofs_row == ndofs_row
    @assert self.ndofs_col == ndofs_col
    self.m[:] .= zero(T)
    return self
end

function assemble!(self::SysmatAssemblerReduced{T}, mat::FMat{T}, dofnums_row::FIntVec, dofnums_col::FIntVec) where {T<:Number}
    R, C = size(mat)
    @assert R == C "The matrix must be square"
    @assert dofnums_row == dofnums_col "The degree of freedom numbers must be the same for rows and columns"
    for i in 1:R
        gi = dofnums_row[i]
        if gi < 1
            mat[i, :] .= zero(T)
            dofnums_row[i] = 1
        end
    end
    for i in 1:C
        gi = dofnums_col[i]
        if gi < 1
            mat[:, i] .= zero(T)
            dofnums_col[i] = 1
        end
    end
    lt = self.t[dofnums_row, :]
    self.m .+= lt' * mat * lt
    return self
end

function assemble!(self::SysmatAssemblerReduced{T}, mat::FMat{T}, dofnums_row::FIntMat, dofnums_col::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums_row), vec(dofnums_col))
end

"""
    makematrix!(self::SysmatAssemblerReduced)

Make a sparse matrix.
"""
function makematrix!(self::SysmatAssemblerReduced)
    return self.m
end

# =============================================================================
# =============================================================================
# =============================================================================
# =============================================================================

"""
    AbstractSysvecAssembler

Abstract type of system vector assembler.
"""
abstract type AbstractSysvecAssembler end;

"""
    startassembly!(self::SysvecAssembler{T},
      ndofs_row::FInt) where {T<:Number}

Start assembly.

The method makes the buffer for the vector assembly. It must be called before
the first call to the method assemble.

`ndofs_row`= Total number of degrees of freedom.
"""
function startassembly!(self::SV,  ndofs_row::FInt) where {SV<:AbstractSysvecAssembler, T<:Number}
end

"""
    assemble!(self::SysvecAssembler{T}, vec::MV,
      dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}

Assemble an elementwise vector.

The method assembles a column element vector using the vector of degree of
freedom numbers for the rows.
"""
function assemble!(self::SV, vec::MV, dofnums::D) where {SV<:AbstractSysvecAssembler, T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}
end

"""
    makevector!(self::SysvecAssembler)

Make the global vector.
"""
function makevector!(self::SV) where {SV<:AbstractSysvecAssembler}
end

"""
    SysvecAssembler

Assembler for the system vector.
"""
mutable struct SysvecAssembler{T<:Number} <: AbstractSysvecAssembler
    F_buffer::Vector{T};
    ndofs::FInt
end

"""
    SysvecAssembler(zero::T=0.0) where {T<:Number}

Construct blank system vector assembler. The vector entries are of type `T`.
"""
function SysvecAssembler(zero::T=0.0) where {T<:Number}
    return SysvecAssembler([zero], 1)
end

"""
    startassembly!(self::SysvecAssembler{T},
      ndofs_row::FInt) where {T<:Number}

Start assembly.

The method makes the buffer for the vector assembly. It must be called before
the first call to the method assemble.

`ndofs_row`= Total number of degrees of freedom.
"""
function startassembly!(self::SysvecAssembler{T},  ndofs_row::FInt) where {T<:Number}
    self.ndofs = ndofs_row
    self.F_buffer = zeros(T,self.ndofs);
end

"""
    assemble!(self::SysvecAssembler{T}, vec::MV,
      dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}

Assemble an elementwise vector.

The method assembles a column element vector using the vector of degree of
freedom numbers for the rows.
"""
function assemble!(self::SysvecAssembler{T}, vec::MV,
  dofnums::D) where {T<:Number, MV<:AbstractArray{T}, D<:AbstractArray{FInt}}
    for i = 1:length(dofnums)
        gi = dofnums[i];
        if (0 < gi <= self.ndofs)
            self.F_buffer[gi] = self.F_buffer[gi] + vec[i];
        end
    end
end

"""
    makevector!(self::SysvecAssembler)

Make the global vector.
"""
function makevector!(self::SysvecAssembler)
  return deepcopy(self.F_buffer);
end

"""
    SysmatAssemblerSparseHRZLumpingSymm{T<:Number} <: AbstractSysmatAssembler

Assembler for a **symmetric lumped square** matrix  assembled from  **symmetric square**
matrices.

Reference: A note on mass lumping and related processes in the finite element method,
E. Hinton, T. Rock, O. C. Zienkiewicz, Earthquake Engineering & Structural Dynamics,
volume 4, number 3, 245--249, 1976.
}

!!! note

    This assembler can compute and assemble diagonalized mass matrices.
    However, if the meaning of the entries of the mass matrix  differs
    (translation versus rotation), the mass matrices will not be computed
    correctly. Put bluntly: it can only be used for homogeneous mass matrices,
    all translation degrees of freedom, for instance.
"""
mutable struct SysmatAssemblerSparseHRZLumpingSymm{T<:Number} <: AbstractSysmatAssembler
    # Type for assembling of a sparse global matrix from elementwise matrices.
    buffer_length:: FInt;
    matbuffer::Vector{T};
    rowbuffer::Vector{FInt};
    colbuffer::Vector{FInt};
    buffer_pointer:: FInt;
    ndofs:: FInt;
end

"""
    SysmatAssemblerSparseHRZLumpingSymm(zer::T=0.0) where {T<:Number}

Construct blank system matrix assembler. The matrix entries are of type `T`.

"""
function SysmatAssemblerSparseHRZLumpingSymm(zer::T=0.0) where {T<:Number}
    return SysmatAssemblerSparseHRZLumpingSymm{T}(0,[zer],[0],[0],0,0)
end

"""
    startassembly!(self::SysmatAssemblerSparseHRZLumpingSymm{T},
      elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt,
      ndofs::FInt, ignore2::FInt) where {T<:Number}

Start the assembly of a symmetric lumped square global matrix.

The method makes buffers for matrix assembly. It must be called before
the first call to the method `assemble!`.
- `elem_mat_nrows`= number of rows in typical element matrix,
- `ignore1`= number of columns in a typical element matrix: equal to `elem_mat_nrows`,
- `elem_mat_nmatrices`= number of element matrices,
- `ndofs_row`= Total number of equations in the row direction,
- `ignore2`= Total number of equations in the column direction: equal to `ndofs_row`.
"""
function startassembly!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, elem_mat_dim::FInt, ignore1::FInt, elem_mat_nmatrices::FInt, ndofs::FInt,ignore2::FInt) where {T<:Number}
    self.buffer_length = elem_mat_nmatrices*elem_mat_dim^2;
    self.rowbuffer = zeros(FInt,self.buffer_length);
    self.colbuffer = zeros(FInt,self.buffer_length);
    self.matbuffer = zeros(T,self.buffer_length);
    self.buffer_pointer = 1;
    self.ndofs = ndofs;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T},
      dofnums::FIntMat, ignore::FIntMat) where {T<:Number}

Assemble a HRZ-lumped square symmetric matrix.

Assembly of a HRZ-lumped square symmetric matrix. The method assembles the
scaled diagonal of the square symmetric matrix using the two vectors of
equation numbers for the rows and columns.
"""
function assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T},  dofnums::FIntVec, ignore::FIntVec) where {T<:Number}

    nrows=length(dofnums); ncolumns=nrows;
    p = self.buffer_pointer
    @assert p+ncolumns*nrows <= self.buffer_length+1
    @assert size(mat) == (nrows, ncolumns)
    @assert nrows == ncolumns # only square matrices allowed
    # Now comes the lumping procedure
    em2 = sum(sum(mat, dims = 1)); # total mass times the number of space dimensions
    dem2 = zero(eltype(mat)); # total mass on the diagonal times the number of space dimensions
    for i in 1:nrows
        dem2 += mat[i, i]
    end
    ffactor = em2/dem2 # total-element-mass compensation factor
    @inbounds for j=1:ncolumns
        @inbounds for i=j:nrows
            if i == j # only the diagonal elements are assembled
                self.matbuffer[p] = mat[i,j]*ffactor # serialized matrix
                self.rowbuffer[p] = dofnums[i];
                self.colbuffer[p] = dofnums[j];
                p=p+1
            end
        end
    end
    self.buffer_pointer=p;
    return self
end

"""
    assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T},
        dofnums::FIntMat, ignore::FIntMat) where {T<:Number}

Assemble an HRZ-lumped square symmetric matrix.
"""
function assemble!(self::SysmatAssemblerSparseHRZLumpingSymm{T}, mat::FMat{T}, dofnums::FIntMat, ignore::FIntMat) where {T<:Number}
    return assemble!(self, mat, vec(dofnums), vec(ignore))
end

"""
    makematrix!(self::SysmatAssemblerSparseHRZLumpingSymm)

Make a sparse HRZ-lumped **symmetric square**  matrix.
"""
function makematrix!(self::SysmatAssemblerSparseHRZLumpingSymm)
    # Make a sparse matrix.
    # The method makes a sparse matrix from the assembly buffers.
    @assert length(self.rowbuffer) >= self.buffer_pointer-1
    @assert length(self.colbuffer) >= self.buffer_pointer-1
    # Here we will go through the rows and columns, and whenever the row or
    # the column refer to indexes outside of the limits of the matrix, the
    # corresponding value will be set to 0 and assembled to row and column 1.
    @inbounds for j=1:self.buffer_pointer-1
        if (self.rowbuffer[j] > self.ndofs) || (self.rowbuffer[j] <= 0)
            self.rowbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
        if (self.colbuffer[j] > self.ndofs) || (self.colbuffer[j] <= 0)
            self.colbuffer[j] = 1;
            self.matbuffer[j] = 0.0
        end
    end
    S = sparse(self.rowbuffer[1:self.buffer_pointer-1],
               self.colbuffer[1:self.buffer_pointer-1],
               self.matbuffer[1:self.buffer_pointer-1],
               self.ndofs, self.ndofs);
    # Construct the other half of the matrix. (Even though this one really should be diagonal!)
    S = S + transpose(S);    # construct the other triangle
    @inbounds for j=1:size(S,1)
        S[j,j]=S[j,j]/2.0;      # the diagonal is there twice; fix it;
    end
    self=SysmatAssemblerSparse(0.0*self.matbuffer[1])# get rid of the buffers
    return S
end

end


text:
function f_ipnlt(t_o, pnlt)

    # Number of time steps
    n_t = size(t_o)[1]

    # Interpolate time, pnlt and C
    dt = 0.5
    n_ip = Int64(ceil((t_o[end]-t_o[1])/dt))

    T = eltype(pnlt)
    t_ip = zeros(T, n_ip)
    for i in range(1, n_ip, step=1)
        t_ip[i] = t_o[1] + (i-1)*dt
    end

    # Interpolate the data
    f_ipnlterp = PCHIPInterpolation.Interpolator(t_o, reshape(pnlt, (n_t, )))
    pnlt_ip = f_ipnlterp.(t_ip)

    # Compute max. PNLT
    pnltm = maximum(pnlt_ip)

    # ICAO Annex 16 procedures (p132 - App. 2-18)
    f_int = 10 .^(pnlt_ip / 10.)

    # Compute integration bounds
    D = 10 * log10(sum(f_int)) - pnltm - 10 * log10(10. / dt)

    # Compute ipnlt
    ipnlt = pnltm + D

    return ipnlt
end

max_stars_repo_name: MIT-LAE/pyNA

max_issues_repo_name: MIT-LAE/pyNA

max_forks_repo_name: MIT-LAE/pyNA


text:
module ParXRec

using LightXML

export parxrec, parrec, xmlrec, SeriesInfo, ImageInfo, RecHdr

Base.@kwdef mutable struct SeriesInfo
    Patient_Name               :: String  = ""
    Examination_Name           :: String  = ""
    Protocol_Name              :: String  = ""
    Examination_Date           :: String  = ""      # TODO: Date
    Examination_Time           :: String  = ""      # TODO: Time
    Series_Data_Type           :: String  = ""
    Aquisition_Number          :: Int32   = 0
    Reconstruction_Number      :: Int32   = 0
    Scan_Duration              :: Float32 = 0
    Max_No_Phases              :: Int32   = 1
    Max_No_Echoes              :: Int32   = 1
    Max_No_Slices              :: Int32   = 1
    Max_No_Dynamics            :: Int32   = 1
    Max_No_Mixes               :: Int16   = 1
    Max_No_B_Values            :: Int32   = 1
    Max_No_Gradient_Orients    :: Int32   = 1
    No_Label_Types             :: Int32   = 0
    Patient_Position           :: String  = ""
    Preparation_Direction      :: String  = ""
    Technique                  :: String  = ""
    Scan_Resolution_X          :: Int16   = 0
    Scan_Resolution_Y          :: Int32   = 0
    Scan_Mode                  :: String  = ""
    Repetition_Times           :: NTuple{2,Float32} = (0,0)
    FOV_AP                     :: Float32 = 0
    FOV_FH                     :: Float32 = 0
    FOV_RL                     :: Float32 = 0
    Water_Fat_Shift            :: Float32 = 0
    Angulation_AP              :: Float32 = 0
    Angulation_FH              :: Float32 = 0
    Angulation_RL              :: Float32 = 0
    Off_Center_AP              :: Float32 = 0
    Off_Center_FH              :: Float32 = 0
    Off_Center_RL              :: Float32 = 0
    Flow_Compensation          :: Bool    = false
    Presaturation              :: Bool    = false
    Phase_Encoding_Velocity    :: NTuple{3,Float32} = (0,0,0)
    MTC                        :: Bool    = false
    SPIR                       :: Bool    = false
    EPI_factor                 :: Int32   = 0
    Dynamic_Scan               :: Bool    = false
    Diffusion                  :: Bool    = false
    Diffusion_Echo_Time        :: Float32 = 0
    PhotometricInterpretation  :: String  = ""
end

Base.@kwdef mutable struct ImageKey
    Slice       :: Int32  = 1
    Echo        :: Int32  = 1
    Dynamic     :: Int32  = 1
    Phase       :: Int32  = 1
    BValue      :: Int32  = 1
    Grad_Orient :: Int32  = 1
    Label_Type  :: String = ""      # Enumeration
    Type        :: String = ""      # Enumeration
    Sequence    :: String = ""      # Enumeration
    Index       :: Int32  = 0
end

Base.@kwdef mutable struct ImageInfo
    Key                        :: ImageKey
    Pixel_Size                 :: UInt16  = 0
    Scan_Percentage            :: Float64 = 0
    Resolution_X               :: UInt16  = 0
    Resolution_Y               :: UInt16  = 0
    Rescale_Intercept          :: Float64 = 0
    Rescale_Slope              :: Float64 = 0
    Scale_Slope                :: Float32 = 0
    Window_Center              :: Float64 = 0
    Window_Width               :: Float64 = 0
    Slice_Thickness            :: Float64 = 0
    Slice_Gap                  :: Float64 = 0
    Display_Orientation        :: String  = ""      # Enumeration
    fMRI_Status_Indication     :: Int16   = 0
    Image_Type_Ed_Es           :: String  = ""      # Enumeration
    Pixel_Spacing              :: NTuple{2,Float64} = (0,0)
    Echo_Time                  :: Float64 = 0
    Dyn_Scan_Begin_Time        :: Float32 = 0
    Trigger_Time               :: Float64 = 0
    Diffusion_B_Factor         :: Float32 = 0
    No_Averages                :: Float64 = 0
    Image_Flip_Angle           :: Float64 = 0
    Cardiac_Frequency          :: Int32   = 0
    Min_RR_Interval            :: Int32   = 0
    Max_RR_Interval            :: Int32   = 0
    TURBO_Factor               :: Int32   = 0
    Inversion_Delay            :: Float64 = 0
    Contrast_Type              :: String  = ""
    Diffusion_Anisotropy_Type  :: String  = ""
    Diffusion_AP               :: Float32 = 0
    Diffusion_FH               :: Float32 = 0
    Diffusion_RL               :: Float32 = 0
    Angulation_AP              :: Float64 = 0
    Angulation_FH              :: Float64 = 0
    Angulation_RL              :: Float64 = 0
    Offcenter_AP               :: Float64 = 0
    Offcenter_FH               :: Float64 = 0
    Offcenter_RL               :: Float64 = 0
    Slice_Orientation          :: String  = ""      # Enumeration
    Image_Planar_Configuration :: UInt16  = 0
    Samples_Per_Pixel          :: UInt16  = 0
end

mutable struct RecHdr
    series :: SeriesInfo
    images :: Vector{ImageInfo}
end

####
#### ParXRec
####

parxrec(path::AbstractString; load_data::Bool = true) =
    parxrec(Float64, path, load_data = load_data)

function parxrec(T::DataType, path::AbstractString; load_data::Bool = true)
    hdrpath, recpath = parxrec_paths(path, load_data)
    ext = splitext(hdrpath)[2]
    hdr = lowercase(ext) == ".par" ? readpar(hdrpath) : readxml(hdrpath)
    return load_data ? (hdr, readrec(T, recpath, hdr)) : hdr
end

parrec(path::AbstractString; load_data::Bool = true) =
    parrec(Float64, path, load_data = load_data)

function parrec(T::DataType, path::AbstractString; load_data::Bool = true)
    hdrpath, recpath = parxrec_paths(path, load_data, ".PAR")
    hdr = readpar(hdrpath)
    return load_data ? (hdr, readrec(T, recpath, hdr)) : hdr
end

xmlrec(path::AbstractString; load_data::Bool = true) =
    xmlrec(Float64, path, load_data = load_data)

function xmlrec(T::DataType, path::AbstractString; load_data::Bool = true)
    hdrpath, recpath = parxrec_paths(path, load_data, ".XML")
    hdr = readxml(hdrpath)
    return load_data ? (hdr, readrec(T, recpath, hdr)) : hdr
end

function parxrec_paths(path::AbstractString, load_data::Bool, ext = "")
    isfile(path) || error("file $path does not exist")
    path1, ext1 = splitext(path)

    if lowercase(ext1) == ".xml" || lowercase(ext1) == ".par"
        ext == "" || lowercase(ext1) == lowercase(ext) ||
            error("wrong file extension: $path (expected $ext)")

        hdrpath = path
        recpath = load_data ? (
            isfile(path1 * ".REC") ? path1 * ".REC" :
            isfile(path1 * ".rec") ? path1 * ".rec" :
            error(".REC file $path1 does not exist")
        ) : ""

    elseif lowercase(ext1) == ".rec"
        recpath = path
        hdrpath = ext == "" ? (
            isfile(path1 * ".PAR") ? path1 * ".PAR" :
            isfile(path1 * ".XML") ? path1 * ".XML" :
            isfile(path1 * ".par") ? path1 * ".par" :
            isfile(path1 * ".xml") ? path1 * ".xml" :
            error(".PAR/.XML file $path1 does not exist")

        ) : lowercase(ext) == ".par" ? (
            isfile(path1 * ".PAR") ? path1 * ".PAR" :
            isfile(path1 * ".par") ? path1 * ".par" :
            error(".PAR file $path1 does not exist")

        ) : lowercase(ext) == ".xml" ? (
            isfile(path1 * ".XML") ? path1 * ".XML" :
            isfile(path1 * ".xml") ? path1 * ".xml" :
            error(".XML file $path1 does not exist")

        ) : error("file extension $ext is not a valid Philips Rec extension")

    else
        error("file $path does not have valid Philips Rec extension")
    end

    return hdrpath, recpath
end

####
#### .REC
####

readrec(path::AbstractString, hdr::RecHdr) =
    readrec(Float64, path, hdr)

function readrec(T::DataType, path::AbstractString, hdr::RecHdr)
    validate(hdr)
    hdr = hdr.images

    prec = Dict{Int,DataType}(
        8  => UInt8,
        16 => UInt16,
        32 => Float32,
        64 => Float64,
    )

    nx       = Int(hdr[1].Resolution_X)
    ny       = Int(hdr[2].Resolution_X)
    nz       = maximum((i -> i.Key.Slice).(hdr))
    echoes   = unique((i -> i.Key.Echo).(hdr))
    dynamics = unique((i -> i.Key.Dynamic).(hdr))
    phases   = unique((i -> i.Key.Phase).(hdr))
    bvals    = unique((i -> i.Key.BValue).(hdr))
    gorients = unique((i -> i.Key.Grad_Orient).(hdr))
    types    = Dict(t => i for (i, t) in enumerate(unique((i -> i.Key.Type).(hdr))))
    labels   = Dict(l => i for (i, l) in enumerate(unique((i -> i.Key.Label_Type).(hdr))))
    seqs     = Dict(s => i for (i, s) in enumerate(unique((i -> i.Key.Sequence).(hdr))))

    if length(unique((i -> i.Pixel_Size).(hdr))) == 1
        sz = (nx * ny * hdr[1].Pixel_Size÷8 * length(hdr)) - filesize(path)
        sz < 0 && @warn "REC file larger than expected from hdr"
        sz > 0 && @warn "REC file smaller than expected from hdr"
    end

    data = zeros(T,
        nx, ny, nz, length.((echoes, dynamics, phases, bvals, gorients, labels, types, seqs))...
    )

    open(path, "r") do io
        for (i, img) in enumerate(hdr)
            ps = Int(img.Pixel_Size)

            seek(io, img.Key.Index * (nx * ny * ps÷8))
            slice = read!(io, Array{prec[ps]}(undef, nx, ny))

            ss  = inv(img.Scale_Slope)
            ri  = ss * img.Rescale_Intercept / img.Rescale_Slope

            i3  = img.Key.Slice
            i4  = img.Key.Echo
            i5  = img.Key.Dynamic
            i6  = img.Key.Phase
            i7  = img.Key.BValue
            i8  = img.Key.Grad_Orient
            i9  = labels[img.Key.Label_Type]
            i10 = types[img.Key.Type]
            i11 = seqs[img.Key.Sequence]

            @. data[:,:,i3,i4,i5,i6,i7,i8,i9,i10,i11] = ss*slice + ri
        end
    end

    return squeeze(data)
end

####
#### .XML
####

const XMLREC_HEADER           = "PRIDE_V5"
const XMLREC_SERIES_HEADER    = "Series_Info"
const XMLREC_IMAGE_ARR_HEADER = "Image_Array"
const XMLREC_IMAGE_HEADER     = "Image_Info"
const XMLREC_IMAGE_KEY_HEADER = "Key"
const XMLREC_ATTRIB_HEADER    = "Attribute"

function readxml(path::AbstractString)
    doc = parse_file(path)

    xroot = root(doc)
    name(xroot) == XMLREC_HEADER || @warn "unknown XML/REC header: $(name(xroot))"

    # Series Info
    cur = find_element(xroot, XMLREC_SERIES_HEADER)
    cur !== nothing || error("<Series_Info> tag not found")

    Ts = NamedTuple{(fieldnames(SeriesInfo)...,)}(SeriesInfo.types)
    series = SeriesInfo()

    for e in cur[XMLREC_ATTRIB_HEADER]
        key, val = _tryparse(Ts, e)
        val !== nothing ?
            setfield!(series, key, val) :
            @warn "unknown series attribute: $key. Skipping..."
    end

    # Image Info
    cur = find_element(xroot, XMLREC_IMAGE_ARR_HEADER)
    cur !== nothing || error("<Image_Array> tag not found")

    Tk = NamedTuple{(fieldnames(ImageKey)...,)}(ImageKey.types)
    Ti = NamedTuple{(fieldnames(ImageInfo)...,)}(ImageInfo.types)
    images = Vector{ImageInfo}()

    for cur in cur[XMLREC_IMAGE_HEADER]
        k = find_element(cur, XMLREC_IMAGE_KEY_HEADER)
        k !== nothing || error("<Key> tag not found")

        imagekey = ImageKey()
        for e in k[XMLREC_ATTRIB_HEADER]
            key, val = _tryparse(Tk, e)
            val !== nothing ?
                setfield!(imagekey, key, val) :
                @warn "unknown key attribute: $key. Skipping..."
        end

        image = ImageInfo(Key = imagekey)
        for e in cur[XMLREC_ATTRIB_HEADER]
            key, val = _tryparse(Ti, e)
            val !== nothing ?
                setfield!(image, key, val) :
                @warn "unknown image attribute: $key. Skipping..."
        end

        push!(images, image)
    end

    free(doc)

    return RecHdr(series, images)
end

####
#### .PAR
####

const PARSERIES = Dict{Symbol,String}(
    :Patient_Name              => "Patient name",
    :Examination_Name          => "Examination name",
    :Protocol_Name             => "Protocol name",
    :Examination_Date          => "Examination date/time",
    :Examination_Time          => "Examination date/time",
    :Series_Data_Type          => "Series Type",
    :Aquisition_Number         => "Acquisition nr",
    :Reconstruction_Number     => "Reconstruction nr",
    :Scan_Duration             => "Scan Duration",
    :Max_No_Phases             => "Max. number of cardiac phases",
    :Max_No_Echoes             => "Max. number of echoes",
    :Max_No_Slices             => "Max. number of slices",
    :Max_No_Dynamics           => "Max. number of dynamics",
    :Max_No_Mixes              => "Max. number of mixes",
    :Max_No_B_Values           => "Max. number of diffusion",
    :Max_No_Gradient_Orients   => "Max. number of gradient",
    :No_Label_Types            => "Number of label types",
    :Patient_Position          => "Patient position",
    :Preparation_Direction     => "Preparation direction",
    :Technique                 => "Technique",
    :Scan_Resolution_X         => "Scan resolution",
    :Scan_Resolution_Y         => "Scan resolution",
    :Scan_Mode                 => "Scan mode",
    :Repetition_Times          => "Repetition time",
    :FOV_AP                    => "FOV",
    :FOV_FH                    => "FOV",
    :FOV_RL                    => "FOV",
    :Water_Fat_Shift           => "Water Fat",
    :Angulation_AP             => "Angulation",
    :Angulation_FH             => "Angulation",
    :Angulation_RL             => "Angulation",
    :Off_Center_AP             => "Off Centre",
    :Off_Center_FH             => "Off Centre",
    :Off_Center_RL             => "Off Centre",
    :Flow_Compensation         => "Flow",
    :Presaturation             => "Presaturation",
    :Phase_Encoding_Velocity   => "Phase encoding",
    :MTC                       => "MTC",
    :SPIR                      => "SPIR",
    :EPI_factor                => "EPI factor",
    :Dynamic_Scan              => "Dynamic scan",
    :Diffusion                 => "Diffusion   ",
    :Diffusion_Echo_Time       => "Diffusion echo time",
)

const PARKEY = Dict{Symbol,Int}(
    :Slice       => 1,
    :Echo        => 2,
    :Dynamic     => 3,
    :Phase       => 4,
    :BValue      => 42,
    :Grad_Orient => 43,
    :Label_Type  => 49,
    :Type        => 5,
    :Sequence    => 6,
    :Index       => 7,
)

const PARIMAGE = Dict{Symbol,Int}(
    :Pixel_Size                 => 8,
    :Scan_Percentage            => 9,
    :Resolution_X               => 10,
    :Resolution_Y               => 11,
    :Rescale_Intercept          => 12,
    :Rescale_Slope              => 13,
    :Scale_Slope                => 14,
    :Window_Center              => 15,
    :Window_Width               => 16,
    :Slice_Thickness            => 23,
    :Slice_Gap                  => 24,
    :Display_Orientation        => 25,
    :fMRI_Status_Indication     => 27,
    :Image_Type_Ed_Es           => 28,
    :Pixel_Spacing              => 29,  # [29, 30]
    :Echo_Time                  => 31,
    :Dyn_Scan_Begin_Time        => 32,
    :Trigger_Time               => 33,
    :Diffusion_B_Factor         => 34,
    :No_Averages                => 35,
    :Image_Flip_Angle           => 36,
    :Cardiac_Frequency          => 37,
    :Min_RR_Interval            => 38,
    :Max_RR_Interval            => 39,
    :TURBO_Factor               => 40,
    :Inversion_Delay            => 41,
    :Contrast_Type              => 44,
    :Diffusion_Anisotropy_Type  => 45,
    :Diffusion_AP               => 46,
    :Diffusion_FH               => 47,
    :Diffusion_RL               => 48,
    :Angulation_AP              => 17,
    :Angulation_FH              => 18,
    :Angulation_RL              => 19,
    :Offcenter_AP               => 20,
    :Offcenter_FH               => 21,
    :Offcenter_RL               => 22,
    :Slice_Orientation          => 26,
)

function readpar(path::AbstractString)
    doc = read(path, String)

    # Series Info / General Information
    Ts = NamedTuple{(fieldnames(SeriesInfo)...,)}(SeriesInfo.types)
    series = SeriesInfo()

    # lines starting with `.`
    buf = [strip.(split(x.match, ": ")) for x in eachmatch(r"(?<=\n\.).*?\n", doc)];

    for (k, v) in PARSERIES
        str = buf[findfirst(a -> occursin(v, a[1]), buf)][2]
        val =
            k == :Examination_Date ?
                strip(split(str, "/")[1]) |> string :
            k == :Examination_Time ?
                strip(split(str, "/")[2]) |> string :
            k ∈ [:Scan_Resolution_X, :FOV_AP, :Angulation_AP, :Off_Center_AP] ?
                _parse(Ts[k], split(str)[1]) :
            k ∈ [:Scan_Resolution_Y, :FOV_FH, :Angulation_FH, :Off_Center_FH] ?
                _parse(Ts[k], split(str)[2]) :
            k ∈ [:FOV_RL, :Angulation_RL, :Off_Center_RL] ?
                _parse(Ts[k], split(str)[3]) :
            # else
                _parse(Ts[k], str)

        setfield!(series, k, val)
    end

    # Image Info
    Tk = NamedTuple{(fieldnames(ImageKey)...,)}(ImageKey.types)
    Ti = NamedTuple{(fieldnames(ImageInfo)...,)}(ImageInfo.types)
    images = Vector{ImageInfo}()

    # lines starting with ` ` or digit
    buf = [split(x.match) for x in eachmatch(r"(?<=\n)( |\d+).*?\n", doc)];

    for line in buf
        imagekey = ImageKey()
        for (k, v) in PARKEY
            setfield!(imagekey, k, _parse(Tk[k], line[v]))
        end
        image = ImageInfo(Key = imagekey)
        for (k, v) in PARIMAGE
            k == :Pixel_Spacing ?
                setfield!(image, k, _parse(Ti[k], line[v]*" "*line[v+1])) :
                setfield!(image, k, _parse(Ti[k], line[v]))
        end
        push!(images, image)
    end

    return RecHdr(series, images)
end

####
#### Utilities
####

function validate(hdr::RecHdr)
    series = hdr.series
    images = hdr.images

    keys = (i -> i.Key).(images)
    keymax = (v) -> maximum(getfield.(keys, v))
    keylength = (v) -> length(unique(getfield.(keys, v)))
    imglength = (v) -> length(unique(getfield.(images, v)))

    vals = [
        (series.Max_No_Slices, :Slice, "slices"),
        (series.Max_No_Echoes, :Echo, "echoes"),
        (series.Max_No_Dynamics, :Dynamic, "dynamics"),
        (series.Max_No_Phases, :Phase, "phases"),
        (series.Max_No_B_Values, :BValue, "bvalues"),
        (series.Max_No_Gradient_Orients, :Grad_Orient, "gradient orientations"),
    ]

    # Check whether data size in series info is consistent with image info
    wstr1 = (f, s, i) ->
        "maximum number of $f not matching: $s (series info) vs. $i (image info)"

    for (v1,v2,v3) in vals
        v1 == keymax(v2) || @warn wstr1(v3, v1, keymax(v2))
    end

    series.No_Label_Types + 1 == keylength(:Label_Type) ||
        @warn wstr1("label types", series.No_Label_Types, keylength(:Label_Type))

    # Non-exhaustive image info check
    wstr2 = (f, m, l) ->
        "inconsistent number of $f in image info: $m (maximum) vs. $l (length(unique))"

    for (_, v1, v2) in vals
        keymax(v1) == keylength(v1) ||
            @warn wstr2(v2, keymax(v1), keylength(v1))
    end

    keymax(:Index) + 1 == keylength(:Index) ||
        @warn wstr2("indices", keymax(:Index), keylength(:Index))

    imglength(:Pixel_Size) == 1 ||
        @warn "multiple values for `Pixel Size` found"

    imglength(:Slice_Thickness) == 1 ||
        @warn "multiple values for `Slice Thickness` found"

    imglength(:Slice_Gap) == 1 ||
        @warn "multiple values for `Slice Gap` found"

    length(unique(first.(getfield.(images, :Pixel_Spacing)))) == 1 ||
        @warn "multiple values for `Pixel Spacing` (x) found"

    length(unique(last.(getfield.(images, :Pixel_Spacing)))) == 1 ||
        @warn "multiple values for `Pixel Spacing` (y) found"

    imglength(:Resolution_X) == 1 ||
        error("multiple values for `Resolution X` found")

    imglength(:Resolution_Y) == 1 ||
        error("multiple values for `Resolution Y` found")

    return nothing
end

function squeeze(x::Array{T}) where {T}
    N = count(n -> n > 1, size(x))
    dims = ([d for d in 1:ndims(x) if size(x, d) == 1]...,)
    return dropdims(x, dims = dims) :: Array{T,N}
end

_parse(::Type{String}, val::AbstractString) = string(val)

_parse(::Type{T}, val::AbstractString) where {T<:Tuple} =
    (_parse.(T.types, split(val))...,)

_parse(::Type{Bool}, val::AbstractString) =
    (val == "1" || uppercase(val) == "Y") ? true :
    (val == "0" || uppercase(val) == "N") ? false : _parse(Bool, val)

_parse(type, str; kwargs...) = Base.parse(type, str; kwargs...) # fallback

function _tryparse(T::NamedTuple, e::XMLElement)
    key = attribute(e, "Name") |> x -> replace(strip(x), " " => "_") |> Symbol
    val = haskey(T, key) ? _parse(T[key], strip(content(e))) : nothing
    return key, val
end

_tryparse(type, str; kwargs...) = Base.tryparse(type, str; kwargs...) # fallback

end # module


text:
# First generate some random data and test functions in Base on it
const NVALS = 1000

input = Dict(
    t=>[
        [ (randindomain(t, NVALS, domain),) for (_, _, domain) in base_unary_real ];
        [ (randindomain(t, NVALS, domain1), randindomain(t, NVALS, domain2))
            for (_, _, domain1, domain2) in base_binary_real ]
    ]
    for t in (Float32, Float64)
)

fns = [[x[1:2] for x in base_unary_real]; [x[1:2] for x in base_binary_real]]

# output = Dict(t=>[fns[i](input[t][i]...) for i = 1:length(fns)] for t in (Float32, Float64))

@testset "Definitions and Comparison with Base for Reals" begin

  for t in (Float32, Float64), i = 1:length(fns)

    base_fn = eval(:($(fns[i][1]).$(fns[i][2])))
    vml_fn = eval(:(IntelVectorMath.$(fns[i][2])))
    # vml_fn! = eval(:(IntelVectorMath.$(fns[i][2])!))

    Test.@test which(vml_fn, typeof(input[t][i])).module == IntelVectorMath

    # Test.test_approx_eq(output[t][i], fn(input[t][i]...), "Base $t $fn", "IntelVectorMath $t $fn")
    Test.@test vml_fn(input[t][i]...) ≈ base_fn.(input[t][i]...)

  end

end

@testset "Error Handling and Settings" begin

  # Verify that we still throw DomainErrors
  Test.@test_throws DomainError IntelVectorMath.sqrt([-1.0])
  Test.@test_throws DomainError IntelVectorMath.log([-1.0])

  # Setting accuracy
  vml_set_accuracy(VML_LA)
  Test.@test vml_get_accuracy() == VML_LA

  vml_set_accuracy(VML_EP)
  Test.@test vml_get_accuracy() == VML_EP

end

@testset "@overload macro" begin

    @test IntelVectorMath.exp([1.0]) ≈ exp.([1.0])
    @test_throws MethodError Base.exp([1.0])
    @test (@overload log exp) isa String
    @test Base.exp([1.0]) ≈ exp.([1.0])

    @test_throws MethodError Base.atan([1.0], [2.0])
    @test (@overload atan) isa String
    @test Base.atan([1.0], [2.0]) ≈ atan.([1.0], [2.0])

end


text:
using Distributions
using PyPlot
srand(1)
# generate 100 observations from correlated normal variates
n = 100
dist = MvNormal([0.0; 0.0], [1.0 0.5; 0.5 1.0])
r = rand(dist, n)'
# create 100 000 bootstrap replications
# and fetch time and memory used
@time bootcor = Float64[cor(r[sample(1:n, n),:])[1, 2] for i in 1:10^5]
# calculate kernel density estimator
kdeboot = kde(bootcor)
# plot results
plt.hist(bootcor, 50, normed = 1)
plot(kdeboot.x, kdeboot.density, color = "y", linewidth = 3.0)
axvline(0.5, color = "r", linewidth = 3.0)
savefig("corboot.pdf", format = "pdf")
# save results to pd


text:
using EarthDataLab, NetCDF, YAXArrays
using Test

newcubedir = mktempdir()
YAXdir(newcubedir)
# Download Cube subset
c = esdc()

cgermany = c[
  region = "Germany",
  var = ["gross", "net_ecosystem", "air_temperature_2m", "terrestrial_ecosystem", "land_surface_temperature"],
  time = 2000:2010
]
savecube(cgermany,"germanycube",
  chunksize=Dict("lon"=>20,"lat"=>20,"time"=>92))
YAXArrays.YAXDefaults.cubedir[] = joinpath(newcubedir,"germanycube")
include("access.jl")
include("analysis.jl")
#include("artype.jl")
include("transform.jl")
include("remap.jl")
include("tabletocube.jl")


text:
"""
	gmtwich(cmd0::String="", arg1=[], kwargs...)

Time domain filtering of 1-D data tables.

Full option list at [`gmtwich`](http://gmt.soest.hawaii.edu/doc/latest/gmtwich.html)

Parameters
----------

- **A** : **with_permissions** : -- Bool or [] --

    Only consider files that the user has permission to read [Default consider all files found].
    [`-A`](http://gmt.soest.hawaii.edu/doc/latest/gmtwich.html#a)
- **C** : **confirm** : -- Bool or [] --

    Instead of reporting the paths, print the confirmation Y if the file is found and N if it is not.
    [`-C`](http://gmt.soest.hawaii.edu/doc/latest/gmtwich.html#c)
- **D** : **report_dir** : -- Bool or [] --

    Instead of reporting the paths, print the directories that contains the files.
    [`-D`](http://gmt.soest.hawaii.edu/doc/latest/gmtwich.html#d)
- **G** : **download** : -- Str or [] --      Flags = [c|l|u]

    If a file argument is a downloadable file (either a full URL, a @file for downloading from
    the GMT Site Cache, or @earth_relief_*.grd) we will try to download the file if it is not
    found in your local data or cache dirs.
    [`-G`](http://gmt.soest.hawaii.edu/doc/latest/gmtwich.html#g)
- $(GMT.opt_V)
"""
function gmtwich(cmd0::String="", arg1=[]; kwargs...)

	length(kwargs) == 0 && return monolitic("gmtwich", cmd0, arg1)	# Speedy mode

	d = KW(kwargs)
	cmd = parse_V_params("", d)

    cmd = add_opt(cmd, 'A', d, [:A :with_permissions])
    cmd = add_opt(cmd, 'C', d, [:C :confirm])
	cmd = add_opt(cmd, 'D', d, [:D :report_dir])
	cmd = add_opt(cmd, 'G', d, [:G :download])

	cmd, got_fname, arg1 = find_data(d, cmd0, cmd, 1, arg1)
	return common_grd(d, cmd, got_fname, 1, "gmtwich", arg1)		# Finish build cmd and run it
end

# ---------------------------------------------------------------------------------------------------
gmtwich(arg1=[], cmd0::String=""; kw...) = gmtwich(cmd0, arg1; kw...)


text:
# # Projects

"""
    struct Project

A `struct` that holds details of a "project", which is defined as

> a configuration file (in TOML format), along with a collection of associated
> source and support files.

"Projects" can be created in two ways, by `Module`, or configuration path.
"""
Base.@kwdef mutable struct Project
    path::Union{Nothing,AbstractPath} = nothing
    env::Dict{String,Any} = Dict()
    tree::FileTree = maketree("." => [])
    pages::Vector{AbstractPath} = []
    loaded::Dict{AbstractPath,Project} = Dict()
    globals::Dict{String,Any} = Dict()
end
Base.show(io::IO, p::Project) = print(io, "$Project($(p.path))")

# ## [`Project`](#) Constructors

"""
    Project(mod)

Create a new [`Project`](#) object from the given module `mod`.
"""
Project(mod::Module; kws...) = Project(findproject(mod); kws...)

"""
    Project(path)

Create a new [`Project`](#) object from the given configuration `path`. The
`path` must be a TOML file.
"""
Project(path::AbstractString; kws...) = Project(Path(path); kws...)

function Project(path::AbstractPath; loaded=Dict{AbstractPath,Project}(), globals=Dict{String,Any}())
    if haskey(loaded, path)
        return loaded[path]
    elseif isfile(path)
        path = abspath(path)
        cd(dirname(path)) do
            path = canonicalize(path)
            env = loadtoml(path, globals)
            tree = loadtree(env, path)
            tree, env = loadtheme(tree, env)
            tree, pages = loadpages(tree, env)
            tree, pages = loaddocs(tree, env, pages)
            loaded[path] = Project(
                path    = path,
                env     = env,
                tree    = tree,
                pages   = pages,
                loaded  = loaded,
                globals = globals,
            )
            return loaded[path]
        end
    else
        ## Raise error only when at 'toplevel' of project creation.
        return isempty(loaded) ? error("not a file: '$path'.") : Project()
    end
end
Project((name, uuid)::Pair; kws...) = Project(findmodule(name, uuid); kws...)
Project(project::Project; kws...) = project
Project(::Nothing; kws...) = nothing

# ## Helper Functions

"""
    findproject(mod)

Returns the path to the configuration file for the module `mod`.

When no path can be found, for example `Base` and it's modules, then `nothing`
is returned instead.
"""
function findproject(mod::Module)
    root = Base.moduleroot(mod)
    meth = first(methods(root.eval))
    file = string(meth.file)
    isabspath(file) || return nothing
    dir = dirname(dirname(realpath(file)))
    for toml in ("Project.toml", "JuliaProject.toml")
        f = joinpath(dir, toml)
        isfile(f) && return f
    end
    return nothing
end
findproject(id::Base.PkgId) = findproject(findmodule(id))

function findmodule(id::Base.PkgId)
    haskey(Base.loaded_modules, id) || Base.require(id)
    return Base.loaded_modules[id]
end
findmodule(name::AbstractString, uuid::AbstractString) = findmodule(Base.PkgId(Base.UUID(uuid), name))
findmodule(env::AbstractDict) = findmodule(env["name"], get(env, "uuid", nothing))
findmodule(name::AbstractString, ::Nothing) = nothing

function findmodules(env::AbstractDict)
    roots = env["publish"]["modules"]
    mods = Set{Module}()
    if isempty(roots)
        mod = findmodule(env)
        mod isa Module && push!(mods, mod)
    else
        for root in roots
            bind = binding(root)
            if Docs.defined(bind)
                push!(mods, Docs.resolve(bind))
            end
        end
    end
    return mods
end

"""
    update!(project)

Reload contents of the given `project`.
"""
function update!(p::Project)
    ## TODO: make this more efficient, only update parts of project that change.
    delete!(p.loaded, p.path)
    q = Project(string(p.path); loaded=p.loaded, globals=p.globals)
    p.path = q.path
    p.env = q.env
    p.tree = q.tree
    p.pages = q.pages
    p.loaded = q.loaded
    p.globals = q.globals
    p.loaded[p.path] = p
    return p
end

function loaddeps!(project::Project)
    @info "loading project dependencies."
    for each in project.env["deps"]
        Project(each; loaded=project.loaded, globals=project.globals)
    end
end


text:
module SCOPE
using Libdl # for NI dlls
include("../INSTR_TYPES.jl")
include("../lib_locations.jl")
include("../visa/VISA.jl")
#
include("SCOPE_common.jl")
# Lecroy
# specifics by manufacturer/model
include("./LecroyHDO6054A/LCR6xxx_Acquisition.jl")
include("./LecroyHDO6054A/LCR6xxx_Horizontal.jl")
include("./LecroyHDO6054A/LCR6xxx_Math.jl")
include("./LecroyHDO6054A/LCR6xxx_Measurements.jl")
include("./LecroyHDO6054A/LCR6xxx_Misc.jl")
include("./LecroyHDO6054A/LCR6xxx_Synchronization.jl")
include("./LecroyHDO6054A/LCR6xxx_Triggering.jl")
include("./LecroyHDO6054A/LCR6xxx_Vertical.jl")

# PXI
include("NIScope/NIScope.jl")

end #endmodule


text:
using BenchmarkTools
using StatsBase
# function lrIndices(trnidx::Vector{Int},f::DecisionTrees.PooledArraysDTM.PooledArray{String,T,1},subset::Array) where T<:Unsigned
include("c:\\temp\\2.jl") # PDA module
using .PooledArraysDTM

"""
function lrIndices!(trnidx::Vector{Int},f,subset::Array)
modifies trnidx to 'become' the left child index, returns r
l and r are defined by subset. l corresponds to the observations 'in' subset
"""
function lrIndices!(trnidx::Vector{Int}, f, subset::Array)
    if eltype(f.pool) <: Number
    # NOTE: subset may not necessarily start at one here (because after a few splits, certain 'middle' values of the pool may not exist in the current node)
    # still, the definition of the split should be of the form x<subset[end] -> go left (otherwise go right)
    # therefore we can consider this function for numerical variables
        return lrIndicesForContiguousSubsetStartingAtONE!(trnidx, f, subset)
    else
        if isContiguous(subset)
            if (length(subset) > 0) && (isone(subset[1]))
                return lrIndicesForContiguousSubsetStartingAtONE!(trnidx, f, subset)
            else
                return lrIndicesForContiguousSubsetNOTStartingAtONE!(trnidx, f, subset)
            end
        else
            return lrIndicesDefault!(trnidx, f, subset)
        end
    end
end

"""
function lrIndices!(trnidx::Vector{Int},f,subset::Array)
modifies trnidx to 'become' the left child index, returns r
l and r are defined by subset. l corresponds to the observations 'in' subset
"""
function lrIndicesDefault!(trnidx::Vector{Int}, f, subset::Array)
	r = Vector{Int}(undef, 0)
	sizehint!(r, length(trnidx))
	indexSize = length(trnidx)
    for jj = indexSize:-1:1
        @inbounds i = trnidx[jj]
		@inbounds thisref = f.refs[i]
		if !in(thisref, subset)
			pushfirst!(r, i)
			deleteat!(trnidx, jj)
		end
	end
	return r
end

"""
function lrIndicesForContiguousSubsetStartingAtONE!(trnidx::Vector{Int},f,subset::Array)
this is a similar function as lrIndices! (->see the docs of lrIndices!)
Here subset needs to be of the form collect(1:n)
"""
function lrIndicesForContiguousSubsetStartingAtONE!(trnidx::Vector{Int}, f, subset::Array)
	r = Vector{Int}(undef, 0)
	sizehint!(r, length(trnidx))
    subsetEnd = subset[end]
    indexSize = length(trnidx)
    for jj = indexSize:-1:1
        @inbounds i = trnidx[jj]
		@inbounds thisref = f.refs[i]
		if !(thisref <= subsetEnd)
			pushfirst!(r, i)
			deleteat!(trnidx, jj)
		end
	end
	return r
end

"""
this is a similar function as lrIndices! (->see the docs of lrIndices!)
Here subset needs to be of the form collect(m:n)
"""
function lrIndicesForContiguousSubsetNOTStartingAtONE!(trnidx::Vector{Int}, f, subset::Array)
	r = Vector{Int}(undef, 0)
	sizehint!(r, length(trnidx))
    subsetEnd = subset[end]
    indexSize = length(trnidx)
	subsetEnd = subset[end]
	subsetStart = subset[1]
    for jj = indexSize:-1:1
        @inbounds i = trnidx[jj]
		@inbounds thisref = f.refs[i]
		# if !((thisref<=subsetEnd)&&(thisref>=subsetStart))
        if (thisref > subsetEnd) || (thisref < subsetStart)
			pushfirst!(r, i)
			deleteat!(trnidx, jj)
		end
	end
	return r
end

"""
isContiguous(v) returns true if the vector v is of the form collect(n:m)
"""
function isContiguous(subset::Vector)
    mi = subset[1]
    ma = subset[end]
    return (length(subset) == (ma - mi + 1)) && issorted(subset)
end

import Random

ntest = 0
whereStartsAtONE = 0
nn = 15
trni = sample(1:nn, floor(Int, .4 * nn), replace=false)
sort!(trni)
strArr = [randstring(1) for i = 1:nn]
floatArr = float.(abs.(rand(Int, nn) .% 9))
pdStr = PooledArray(strArr)
pdFloat = PooledArray(floatArr)

for klU = 1:400000
    # @show klU
    pd = rand() < .5 ? pdStr : pdFloat

    subsetStrOrFloat = sample(pd.pool, 4, replace=false)
    subset = convert(Vector{eltype(pd.refs)}, findall((in)(subsetStrOrFloat), pd.pool))
    sort!(subset)

    # if isContiguous(subset)

    whereStartsAtONE += isone(subset[1])
    ntest += 1

    trniWillBeModified = deepcopy(trni)
    rrr = lrIndices!(trniWillBeModified, pd, subset)
    if eltype(pd.pool) == Float64
        ll2, rr2 = DecisionTrees.lrIndicesForContiguousSubsetStartingAtONE(trni, pd, subset)
    else
        ll2, rr2 = DecisionTrees.lrIndices(trni, pd, subset)
    end
        #=
        @show ll2,rr2
        @show subset
        @show pd
        @show pd.pool
        @show sort!(unique(pd.refs))
        @show ll2
        @show rrr
        @show trni
        @show trniWillBeModified =#
    @assert ll2 == trniWillBeModified
    @assert rr2 == rrr

        # @btime lrIndices($trni,$pd,$subset)
        # @btime lrIndices2($trni,$pd,$subset)
    # end
end


text:
using Test
include("../src/metric.jl")

@testset "Basic Loss Functions" begin

    @testset "MSE" begin
        @test mse(1, 1) == 0
        @test mse([1,1], [1,1]) == 0
        @test mse([1,1], [3,4]) == 6.5
    end

    @testset "MAE" begin
        @test mae(1, 1) == 0
        @test mae([1,1], [1,1]) == 0
        @test mae([1,1], [3,4]) == 2.5
    end

    @testset "RMSE" begin
        @test rmse(1, 1) == 0
        @test rmse([1,1], [1,1]) == 0
        @test rmse([1,1], [3,4]) == 2.5495097567963922
    end

    @testset "RMAE" begin
        @test rmae(1, 1) == 0
        @test rmae([1,1], [1,1]) == 0
        @test rmae([1,1], [3,4]) == 1.5811388300841898
    end

    @testset "RMSLE" begin
        @test rmsle(1, 1) == 0
        @test rmsle([1,1], [1,1]) == 0
        @test rmsle([100,101], [30,400]) == 1.2785158831395131
    end

    @testset "R2SCORE" begin
        @test r2score(1, 1) == 1
        @test r2score([1,1], [1,1]) == 1
        @test r2score([1,1], [3,4]) == -11.5
    end

    @testset "BINCROSSENTROPY" begin
        @test binarycrossentropy(1, 0.9) == 1
        @test binarycrossentropy([0,0,1], [0.2,0.3,0.8]) == 0.3861431206792365
    end

end


text:
using Random
using LinearAlgebra

"""
Basis Vector = ej
"""
function basisVec(l, i)
    return vcat(zeros(i-1), 1.0, zeros(l-i))
end

"""
Generating a random vector based on Rademacher Distribution
"""
function randomRademacherVector(x::Int64, T::Type)
    o = one(T)
    t = 2*o
    return rand(-o:t:o, x)
end

"""
Lanczos function for eigen vlaues
"""
function lanczos(A, x, m)
    q = (x/norm(x))
    r = similar(q)

    mul!(r,A,q)
    alpha = q' * r

    r = r - alpha .* q

    Alpha = Vector{Float64}(undef, (m+1))
    Beta = Vector{Float64}(undef, m)
    Alpha[1] = alpha

    for j = 2:(m+1)
        Beta[j-1] = norm(r)
        v = q
        q = r/(Beta[j-1])

        mul!(r, A, q)
        r .-= Beta[j-1] .* v
        alpha = q' * r

        r .= r .- alpha .* q

        Alpha[j] = alpha
    end
    T = SymTridiagonal(Alpha, Beta)
    return T
end

"""
Actual Algorithm Implementataion
"""
function slq(A::AbstractMatrix, f::Function, m::Int64, nv::Int64)
    if !(isposdef(A))
        throw("Matrix passed is not a Positive Semi-Definite Martix")
    else
        trace = 0
        for i = 1:nv
            vl = randomRademacherVector(size(A,2), Float64)
            T = lanczos(A, vl, m)

            Y = eigvecs(T)
            theta = eigvals(T)

            # for k = 1:size(theta,1)
            for k = 1:(m+1)
                trace = trace + Y[1,k]^2 * f(theta[k])
            end
        end
        return (size(A, 1)/nv) * trace
    end
end