text:
# Test variational information

using Base.Test
using Clustering

a1 = [1, 1, 1, 2, 2, 2, 3, 3, 3, 3]
a2 = [1, 1, 1, 1, 2, 2, 2, 2, 2, 2]
a3 = [3, 3, 3, 2, 2, 2, 1, 1, 1, 1]

(ARI, RI, MI, HI) = randindex(a1, a1)

@test_approx_eq_eps ARI 1.0 1.0e-12
@test_approx_eq_eps RI  1.0 1.0e-12
@test_approx_eq_eps MI  0.0 1.0e-12
@test_approx_eq_eps HI  1.0 1.0e-12

(ARI, RI, MI, HI) = randindex(a1, a3)

@test_approx_eq_eps ARI 1.0 1.0e-12
@test_approx_eq_eps RI  1.0 1.0e-12
@test_approx_eq_eps MI  0.0 1.0e-12
@test_approx_eq_eps HI  1.0 1.0e-12

(ARI, RI, MI, HI) = randindex(a1, a2)

@test_approx_eq_eps ARI 0.403669 1.0e-5
@test_approx_eq_eps RI  0.711111 1.0e-5
@test_approx_eq_eps MI  0.288888 1.0e-5
@test_approx_eq_eps HI  0.422222 1.0e-5

@test randindex(a1, a2) == randindex(a2, a1)


text:
module TestUnionArrays
using Test

env_test_cuda = lowercase(get(ENV, "UNIONARRAYS_JL_TEST_CUDA", "auto"))

const TEST_CUDA = if env_test_cuda == "true"
    import CUDA
    true
elseif env_test_cuda == "auto"
    try
        import CUDA
        CUDA.has_cuda_gpu()
    catch
        false
    end
else
    false
end

const TEST_GPU = TEST_CUDA

TEST_CUDA && CUDA.allowscalar(false)

find_test(subdir = "") = sort([
    joinpath(subdir, file) for file in readdir(joinpath(@__DIR__, subdir)) if
    match(r"^test_.*\.jl$", file) !== nothing
])

@testset "$file" for file in find_test()
    TEST_GPU || include(file)
end

@testset "$file" for file in find_test("cuda")
    if VERSION < v"1.6-"
        basename(file) == "test_kernels.jl" && continue
        basename(file) == "test_doctest.jl" && continue
    end
    TEST_CUDA && include(file)
end

end  # module

ext: jl


text:
meanVect = [27.155405405405407, 26.163835263835264]
covMat = [16.1253895583728 13.046961200891612; 13.046961200891612 12.367253570765161]

max_stars_repo_name: newptcai/StatsWithJuliaBook

max_issues_repo_name: newptcai/StatsWithJuliaBook

max_forks_repo_name: newptcai/StatsWithJuliaBook


text:
# Plotting results

using Plots
plotly();
print("Plotting results...")

ArraysToPlot = Dict{Int, Any}() # Priority , Arrays of RequirementFulfillment corresponding to the priority
ArraysLabels = Dict{Int, Any}() # Priority , Array of lablels
for i in 1:length(MPObjectives)
    try
        ArraysToPlot[MPObjectives[i].PlotNb] = [ArraysToPlot[MPObjectives[i].PlotNb] RequirementFulfillment[:,i]]
        ArraysLabels[MPObjectives[i].PlotNb] = [ArraysLabels[MPObjectives[i].PlotNb] reduce(*,MPObjectives[i].Objectives)]
    catch
        if MPObjectives[i].PlotNb > 0
            ArraysToPlot[MPObjectives[i].PlotNb] = RequirementFulfillment[:,i]
            ArraysLabels[MPObjectives[i].PlotNb] = [reduce(*,MPObjectives[i].Objectives)]
            if ArraysLabels[MPObjectives[i].PlotNb][1] == ""
                ArraysLabels[MPObjectives[i].PlotNb] = "Total Population"
            end
        end
    end
end
for p in ArraysToPlot
    gui(plot(1:NBYears, p[2], size=(800,600), label = ArraysLabels[p[1]], lw = 3, ylim = (0,1.1*maximum(p[2])), title = "$(SimulationName) NbPlot : $(p[1])"))
end

for i in 1:length(MPObjectives)
    if MPObjectives[i].PlotNb == 0
        titles = ((MPObjectives[i].Objectives[1] == "")? "Total population": reduce(*,MPObjectives[i].Objectives))
        p1 = plot(1:NBYears, RequirementFulfillment[:,i], size=(800,600), lw = 3, ylim = (0,1.1*maximum(RequirementFulfillment[:,i])), title = "$(SimulationName) $(titles)")
        plot!(p1, [1, NBYears], [(1 + MPObjectives[i].InitTolerance)*MPObjectives[i].Number, (1 + MPObjectives[i].EndTolerance)*MPObjectives[i].Number], color = :red, legend=false)
        plot!(p1, [1, NBYears], [(1 - MPObjectives[i].InitTolerance)*MPObjectives[i].Number, (1 - MPObjectives[i].EndTolerance)*MPObjectives[i].Number], color = :red, legend=false)
        gui(p1)
    end
end
println(" End.")

max_stars_repo_name: Omazaria/CareerPathBeDef

max_issues_repo_name: Omazaria/CareerPathBeDef

max_forks_repo_name: Omazaria/CareerPathBeDef


text:
@inline base_aggregate(f::BaseFunction, s::T, n, x::T,       y::T)       where T =
    base_aggregate(f, s, x, y), n + 1
@inline base_aggregate( ::BaseFunction, s::T, n, x::T,        ::Missing) where T = s, n
@inline base_aggregate( ::BaseFunction, s::T, n,  ::Missing, y::T)       where T = s, n
@inline base_aggregate( ::BaseFunction, s::T, n,  ::Missing,  ::Missing) where T = s, n

@inline function unsafe_base_evaluate(
    f::BaseFunction,
    x::AbstractArray{mT1},
    y::AbstractArray{mT2}
) where mT1 <: Union{T, Missing} where mT2 <: Union{T, Missing} where T <: AbstractFloat

    s, n = base_initiate(f, T), 0
    # s, n = base_initiate(f, promote_type(eltype(x), eltype(y))), 0
    @inbounds for i in eachindex(x, y)
        s, n = base_aggregate(f, s, n, x[i], y[i])
    end

    n > 0 ? base_return(f, s / n) : missing
end

# TODO: ask MLkernels.jl to  incorporate these changes
for orientation in (:row, :col)

    row_oriented = orientation == :row
    dim_obs      = row_oriented ? 1 : 2

    @eval begin

        @inline function allocate_basematrix(
             ::Val{$(Meta.quot(orientation))},
            X::AbstractMatrix{Tm}
        ) where Tm
            Array{Tm}(undef, size(X,$dim_obs), size(X,$dim_obs))
        end

        @inline function allocate_basematrix(
             ::Val{$(Meta.quot(orientation))},
            X::AbstractMatrix{mT1},
            Y::AbstractMatrix{mT2}
        ) where mT1 <: Union{T, Missing} where mT2 <: Union{T, Missing} where T <: AbstractFloat
            Array{Union{T, Missing}}(undef, size(X,$dim_obs), size(Y,$dim_obs))
        end

    end
end

function basematrix!(
    σ::Orientation,
    P::Matrix{Union{T, Missing}},
    f::BaseFunction,
    X::Matrix{mT},
    symmetrize::Bool
) where
    {mT <: Union{T, Missing}} where
    {T <: AbstractFloat}

    n = checkdimensions(σ, P, X)

    @inbounds for j in 1:n
        xj = subvector(σ, X, j)
        for i in 1:j
            xi = subvector(σ, X, i)
            P[i, j] = unsafe_base_evaluate(f, xi, xj)
        end
    end

    symmetrize && copytri!(P, 'U', false)

    return P
end

function basematrix!(
    σ::Orientation,
    P::Matrix{Union{T, Missing}},
    f::BaseFunction,
    X::AbstractMatrix{mT1},
    Y::AbstractMatrix{mT2},
) where
    mT1 <: Union{T, Missing} where
    mT2 <: Union{T, Missing} where
    T <: AbstractFloat

    n, m = checkdimensions(σ, P, X, Y)

    @inbounds for j = 1:m
        yj = subvector(σ, Y, j)
        for i = 1:n
            xi = subvector(σ, X, i)
            P[i,j] = unsafe_base_evaluate(f, xi, yj)
        end
    end

    return P
end


text:
function interpolate(obj::TrivariateGridInterpolator, arg0::Vector{jdouble}, arg1::Vector{jdouble}, arg2::Vector{jdouble}, arg3::Vector{Vector{jdouble}})
    return jcall(obj, "interpolate", TrivariateFunction, (Vector{jdouble}, Vector{jdouble}, Vector{jdouble}, Vector{Vector{jdouble}}), arg0, arg1, arg2, arg3)
end


text:
function isCoprime(a, b)
    if b == 0
        return a == 1
    end
    return isCoprime(b, a % b)
end

max_stars_repo_path: Utility Functions/Coprime.jl

max_issues_repo_path: Utility Functions/Coprime.jl

max_forks_repo_path: Utility Functions/Coprime.jl


text:
function Armijo(nlp, f, g, x, d, τ₀)
  # Simple Armijo backtracking
  hp0 = g' * d
  t = 1.0
  ft = obj(nlp, x + t * d)
  while ft > (f + τ₀ * t * hp0)
    t /= 2.0
    ft = obj(nlp, x + t * d)
  end
end


text:
module IntegratedOptics

using  LinearAlgebra
using  SparseArrays

include("types.jl")
include("Interface/Simulation.jl")
include("LinearAlgebra/LinearAlgebra.jl")
include("Solvers/FDFD/FDFDSolver.jl")
include("Solvers/ModeSolver/FullVectorialModeSolver.jl")
include("Optimization/Optimization.jl")
include("Plots/Plots.jl")

end # module


text:
"""
This part of the code was taken from https://github.com/ADThomas-astro/oxaf/blob/master/oxaf.py .
Credit to A.D. Thomas.
Code was adapted from Xspec for https://arxiv.org/pdf/1611.05165.pdf .
"""

using Qsosed, Test, PyCall

py"""
import numpy as np

def donthcomp_in_python(ear, param):
    param = np.array(param)
    param = np.insert(param,0,0)
    ne = ear.size  # Length of energy bin vector
    # Note that this model does not calculate errors.
    #c     xth is the energy array (units m_e c^2)
    #c     spnth is the nonthermal spectrum alone (E F_E)
    #c     sptot is the total spectrum array (E F_E), = spref if no reflection
    zfactor = 1.0 + param[5]
    #c  calculate internal source spectrum
    #                           blackbody temp,   plasma temp,      Gamma
    xth, nth, spt = _thcompton(param[3] / 511.0, param[2] / 511.0, param[1])
    # The temperatures are normalized by 511 keV, the electron rest energy
    # Calculate normfac:
    xninv = 511.0 / zfactor
    ih = 1
    xx = 1.0 / xninv
    while (ih < nth and xx > xth[ih]):
        ih = ih + 1
    il = ih - 1
    spp = spt[il] + (spt[ih] - spt[il]) * (xx - xth[il]) / (xth[ih] - xth[il])
    normfac = 1.0 / spp

    #c     zero arrays
    photar = np.zeros(ne)
    prim   = np.zeros(ne)
    #c     put primary into final array only if scale >= 0.
    j = 0
    for i in range(0, ne):
        while (j <= nth and 511.0 * xth[j] < ear[i] * zfactor):
            j = j + 1
        if (j <= nth):
            if (j > 0):
                jl = j - 1
                prim[i] = spt[jl] + ((ear[i] / 511.0 * zfactor - xth[jl]) *
                                     (spt[jl + 1] - spt[jl]) /
                                     (xth[jl + 1] - xth[jl])                 )
            else:
                prim[i] = spt[0]
    for i in range(1, ne):
        photar[i] = (0.5 * (prim[i] / ear[i]**2 + prim[i - 1] / ear[i - 1]**2)
                         * (ear[i] - ear[i - 1]) * normfac                    )

    return photar

def _thcompton(tempbb, theta, gamma):
    #c use internally Thomson optical depth
    tautom = np.sqrt(2.250 + 3.0 / (theta * ((gamma + .50)**2 - 2.250))) - 1.50

    # Initialise arrays
    dphdot = np.zeros(900); rel = np.zeros(900); c2 = np.zeros(900)
    sptot  = np.zeros(900); bet = np.zeros(900); x  = np.zeros(900)

    #c JMAX  -  # OF PHOTON ENERGIES
    #c delta is the 10 - log interval of the photon array.
    delta = 0.02
    deltal = delta * np.log(10.0)
    xmin = 1e-4 * tempbb
    xmax = 40.0 * theta
    jmax = min(899, int(np.log10(xmax / xmin) / delta) + 1)

    #c X  -  ARRAY FOR PHOTON ENERGIES
    # Energy array is normalized by 511 keV, the rest energy of an electron
    x[:(jmax + 1)] = xmin * 10.0**(np.arange(jmax + 1) * delta)

    #c compute c2(x), and rel(x) arrays
    #c c2(x) is the relativistic correction to Kompaneets equation
    #c rel(x) is the Klein - Nishina cross section divided by the
    #c Thomson crossection
    for j in range(0, jmax):
        w = x[j]
    #c c2 is the Cooper's coefficient calculated at w1
    #c w1 is x(j + 1 / 2) (x(i) defined up to jmax + 1)
        w1 = np.sqrt(x[j] * x[j + 1])
        c2[j] = (w1**4 / (1.0 + 4.60 * w1 + 1.1 * w1 * w1))
        if (w <= 0.05):
            #c use asymptotic limit for rel(x) for x less than 0.05
            rel[j] = (1.0 - 2.0 * w + 26.0 * w * w * 0.2)
        else:
            z1 = (1.0 + w) / w**3
            z2 = 1.0 + 2.0 * w
            z3 = np.log(z2)
            z4 = 2.0 * w * (1.0 + w) / z2
            z5 = z3 / 2.0 / w
            z6 = (1.0 + 3.0 * w) / z2 / z2
            rel[j] = (0.75 * (z1 * (z4 - z3) + z5 - z6))

    #c the thermal emission spectrum
    jmaxth = min(900, int(np.log10(50 * tempbb / xmin) / delta))
    if (jmaxth > jmax):
       jmaxth = jmax
    planck = 15.0 / (np.pi * tempbb)**4
    dphdot[:jmaxth] = planck * x[:jmaxth]**2 / (np.exp(x[:jmaxth] / tempbb)-1)

    #c compute beta array, the probability of escape per Thomson time.
    #c bet evaluated for spherical geometry and nearly uniform sources.
    #c Between x = 0.1 and 1.0, a function flz modifies beta to allow
    #c the increasingly large energy change per scattering to gradually
    #c eliminate spatial diffusion
    jnr  = int(np.log10(0.10 / xmin) / delta + 1)
    jnr  = min(jnr, jmax - 1)
    jrel = int(np.log10(1 / xmin) / delta + 1)
    jrel = min(jrel, jmax)
    xnr  = x[jnr - 1]
    xr   = x[jrel - 1]
    for j in range(0, jnr - 1):
        taukn = tautom * rel[j]
        bet[j] = 1.0 / tautom / (1.0 + taukn / 3.0)
    for j in range(jnr - 1, jrel):
        taukn = tautom * rel[j]
        arg = (x[j] - xnr) / (xr - xnr)
        flz = 1 - arg
        bet[j] = 1.0 / tautom / (1.0 + taukn / 3.0 * flz)
    for j in range(jrel, jmax):
        bet[j] = 1.0 / tautom

    dphesc = _thermlc(tautom, theta, deltal, x, jmax, dphdot, bet, c2)

    #c     the spectrum in E F_E
    for j in range(0, jmax - 1):
        sptot[j] = dphesc[j] * x[j]**2

    return x, jmax, sptot

def _thermlc(tautom, theta, deltal, x, jmax, dphdot, bet, c2):
    dphesc = np.zeros(900)  # Initialise the output
    a = np.zeros(900); b   = np.zeros(900); c = np.zeros(900)
    d = np.zeros(900); alp = np.zeros(900); u = np.zeros(900)
    g = np.zeros(900); gam = np.zeros(900)

    #c u(x) is the dimensionless photon occupation number
    c20 = tautom / deltal

    #c determine u
    #c define coefficients going into equation
    #c a(j) * u(j + 1) + b(j) * u(j) + c(j) * u(j - 1) = d(j)
    for j in range(1, jmax - 1):
        w1 = np.sqrt( x[j] * x[j + 1] )
        w2 = np.sqrt( x[j - 1] * x[j] )
        #c  w1 is x(j + 1 / 2)
        #c  w2 is x(j - 1 / 2)
        a[j] =  -c20 * c2[j] * (theta / deltal / w1 + 0.5)
        t1 =  -c20 * c2[j] * (0.5 - theta / deltal / w1)
        t2 = c20 * c2[j - 1] * (theta / deltal / w2 + 0.5)
        t3 = x[j]**3 * (tautom * bet[j])
        b[j] = t1 + t2 + t3
        c[j] = c20 * c2[j - 1] * (0.5 - theta / deltal / w2)
        d[j] = x[j] * dphdot[j]

    #c define constants going into boundary terms
    #c u(1) = aa * u(2) (zero flux at lowest energy)
    #c u(jx2) given from region 2 above
    x32 = np.sqrt(x[0] * x[1])
    aa = (theta / deltal / x32 + 0.5) / (theta / deltal / x32 - 0.5)

    #c zero flux at the highest energy
    u[jmax - 1] = 0.0

    #c invert tridiagonal matrix
    alp[1] = b[1] + c[1] * aa
    gam[1] = a[1] / alp[1]
    for j in range(2, jmax - 1):
        alp[j] = b[j] - c[j] * gam[j - 1]
        gam[j] = a[j] / alp[j]
    g[1] = d[1] / alp[1]
    for j in range(2, jmax - 2):
        g[j] = (d[j] - c[j] * g[j - 1]) / alp[j]
    g[jmax - 2] = (d[jmax - 2] - a[jmax - 2] * u[jmax - 1]
                               - c[jmax - 2] * g[jmax - 3]) / alp[jmax - 2]
    u[jmax - 2] = g[jmax - 2]
    for j in range(2, jmax + 1):
        jj = jmax - j
        u[jj] = g[jj] - gam[jj] * u[jj + 1]
    u[0] = aa * u[1]
    #c compute new value of dph(x) and new value of dphesc(x)
    dphesc[:jmax] = x[:jmax] * x[:jmax] * u[:jmax] * bet[:jmax] * tautom

    return dphesc
"""

@testset "Test that they agree with ADT implementation in python" begin
    ear = 10 .^ range(-4, 2, length=50)
    #parameters = [2.27963222215913, 100, 0.01676681194689287, 0, 0]
    for gamma in range(2, 10, length=10)
        for kt_e in range(10, 1000, length=10)
            for t_corona_kev in range(0.001, 0.1, length=10)
                parameters = [gamma, kt_e, t_corona_kev, 0, 0]
                expected = py"donthcomp_in_python"(ear, parameters)
                results = donthcomp(ear, parameters)
                @test expected ≈ results rtol = 1e-3
            end
        end
    end
end


text:
#-------------------------------------------------------------------------------------------------------#
#    The library contains material characteristics needed for heat and moisture transfer calculation    #
#                                                                                #
#    Characteristics are stored in material characteristics functions accesed with the material key     #
#-------------------------------------------------------------------------------------------------------#

#----------------------------------------------------------------------#
#    Moisture storage function (relative humidity to water content)    #

# RHwc(ϕ,a=0) -> "m" is relative humidity (RH) and "a" is a material key #

function RHwc(ϕ,a=0)

    if a==0
        println("ERROR! Please specify ϕaterial key a in function RHwc(ϕ,a)!")

    #Lime-cement mortar, fine
    elseif a==10
        if ϕ>=0&&ϕ<0.30
            12.79/0.3*ϕ
        elseif ϕ>=0.30&&ϕ<0.50
            12.79+3.09/0.2*(ϕ-0.3)
        elseif ϕ>=0.50&&ϕ<0.70
            15.88+5.09/0.2*(ϕ-0.5)
        elseif ϕ>=0.70&&ϕ<0.80
            20.98+4.68/0.1*(ϕ-0.7)
        elseif ϕ>=0.80&&ϕ<=0.90
            25.66+9.62/0.1*(ϕ-0.8)
        elseif ϕ>=0.90&&ϕ<=0.95
            35.28+11.65/0.05*(ϕ-0.9)
        else
            println("ERROR! Relative humidity greater than 95 %!")
        end

    #Mineral wool
    elseif a==11
        if ϕ>=0&&ϕ<0.50
            0.461/0.5*ϕ
        elseif ϕ>=0.50&&ϕ<0.60
            0.461+0.226/0.1*(ϕ-0.5)
        elseif ϕ>=0.60&&ϕ<0.70
            0.687+0.373/0.1*(ϕ-0.6)
        elseif ϕ>=0.70&&ϕ<0.80
            1.06+0.73/0.1*(ϕ-0.7)
        elseif ϕ>=0.80&&ϕ<0.85
            1.79+0.7/0.05*(ϕ-0.8)
        elseif ϕ>=0.85&&ϕ<0.90
            2.49+1.34/0.05*(ϕ-0.85)
        elseif ϕ>=0.90&&ϕ<0.91
            3.83+0.43/0.01*(ϕ-0.9)
        elseif ϕ>=0.91&&ϕ<0.92
            4.26+0.52/0.01*(ϕ-0.91)
        elseif ϕ>=0.92&&ϕ<0.93
            4.78+0.65/0.01*(ϕ-0.92)
        elseif ϕ>=0.93&&ϕ<0.94
            5.43+0.84/0.01*(ϕ-0.93)
        elseif ϕ>=0.94&&ϕ<0.95
            6.27+1.11/0.01*(ϕ-0.94)
        elseif ϕ>=0.95&&ϕ<0.96
            7.38+1.56/0.01*(ϕ-0.95)
        elseif ϕ>=0.96&&ϕ<0.97
            8.94+2.36/0.01*(ϕ-0.96)
        elseif ϕ>=0.97&&ϕ<0.98
            11.3+3.8/0.01*(ϕ-0.97)
        elseif ϕ>=0.98&&ϕ<0.99
            15.1+7.6/0.01*(ϕ-0.98)
        elseif ϕ>=0.99&&ϕ<0.995
            22.7+7.32/0.005*(ϕ-0.99)
        elseif ϕ>=0.995&&ϕ<=100
            30.2+14.6/0.005*(ϕ-0.995)
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Concrete
    elseif a==12
        if ϕ>=0&&ϕ<0.35
            27/0.35*ϕ
        elseif ϕ>=0.35&&ϕ<0.50
            27+7/0.15*(ϕ-0.35)
        elseif ϕ>=0.50&&ϕ<0.70
            34+15.5/0.2*(ϕ-0.5)
        elseif ϕ>=0.70&&ϕ<0.80
            49.5+14.5/0.1*(ϕ-0.7)
        elseif ϕ>=0.80&&ϕ<0.90
            64+24.5/0.1*(ϕ-0.8)
        elseif ϕ>=0.90&&ϕ<0.95
            88.5+20.5/0.05*(ϕ-0.9)
        elseif ϕ>=0.95&&ϕ<0.96
            109+5.5/0.01*(ϕ-0.95)
        elseif ϕ>=0.96&&ϕ<0.97
            114.5+10/0.01*(ϕ-0.96)
        elseif ϕ>=0.97&&ϕ<0.976
            124.5+5.5/0.006*(ϕ-0.97)
        elseif ϕ>=0.976&&ϕ<=1
            130+10/0.024*(ϕ-0.976)
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Reinforced concrete
    elseif a==13
        if ϕ>=0&&ϕ<0.35
            27/0.35*ϕ/100
        elseif ϕ>=0.35&&ϕ<0.50
            27+7/0.15*(ϕ-0.35)
        elseif ϕ>=0.50&&ϕ<0.70
            34+15.5/0.2*(ϕ-0.5)
        elseif ϕ>=0.70&&ϕ<0.80
            49.5+14.5/0.1*(ϕ-0.7)
        elseif ϕ>=0.80&&ϕ<0.90
            64+24.5/0.1*(ϕ-0.8)
        elseif ϕ>=0.90&&ϕ<0.95
            88.5+20.5/0.05*(ϕ-0.9)
        elseif ϕ>=0.95&&ϕ<0.96
            109+5.5/0.01*(ϕ-0.95)
        elseif ϕ>=0.96&&ϕ<0.97
            114.5+10/0.01*(ϕ-0.96)
        elseif ϕ>=0.97&&ϕ<0.976
            124.5+5.5/0.006*(ϕ-0.97)
        elseif ϕ>=0.976&&ϕ<=1
            130+10/0.024*(ϕ-0.976)
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Extruded polystyrene
    elseif a==14
            if ϕ>=0&&ϕ<0.10
                0.009/0.1*ϕ
            elseif ϕ>=0.10&&ϕ<0.50
                0.009+0.071/0.4*(ϕ-0.1)
            elseif ϕ>=0.50&&ϕ<0.80
                0.08+0.23/0.3*(ϕ-0.5)
            elseif ϕ>=0.80&&ϕ<0.90
                0.31+0.35/0.1*(ϕ-0.8)
            elseif ϕ>=0.90&&ϕ<0.97
                0.66+1.29/0.07*(ϕ-0.9)
            else
                println("ERROR! Relative humidity greater than 100 %!")
            end

    #PVC hydro membrane
    elseif a==15
            if ϕ>=0&&ϕ<0.5
                0.000485/0.5*ϕ
            elseif ϕ>=0.5&&ϕ<0.6
                0.000485+0.000239/0.1*(ϕ-0.5)
            elseif ϕ>=0.60&&ϕ<0.70
                0.000724+0.000396/0.1*(ϕ-0.6)
            elseif ϕ>=0.70&&ϕ<0.80
                0.00112+0.00076/0.1*(ϕ-0.7)
            elseif ϕ>=0.80&&ϕ<0.85
                0.00188+0.0074/0.05*(ϕ-0.8)
            elseif ϕ>=0.85&&ϕ<0.9
                0.00262+0.00141/0.05*(ϕ-0.85)
            elseif ϕ>=0.9&&ϕ<0.91
                0.00403+0.00045/0.01*(ϕ-0.9)
            elseif ϕ>=0.91&&ϕ<0.92
                0.00448+0.00055/0.01*(ϕ-0.91)
            elseif ϕ>=0.92&&ϕ<0.93
                0.00503+0.00069/0.01*(ϕ-0.92)
            elseif ϕ>=0.93&&ϕ<0.94
                0.00572+0.00088/0.01*(ϕ-0.93)
            elseif ϕ>=0.94&&ϕ<0.95
                0.0066+0.00117/0.01*(ϕ-0.94)
            elseif ϕ>=0.95&&ϕ<0.96
                0.00777+0.00164/0.01*(ϕ-0.95)
            elseif ϕ>=0.96&&ϕ<0.97
                0.00941+0.00249/0.01*(ϕ-0.96)
            elseif ϕ>=0.97&&ϕ<0.98
                0.0119+0.004/0.01*(ϕ-0.97)
            elseif ϕ>=0.98&&ϕ<0.99
                0.0159+0.008/0.01*(ϕ-0.98)
            elseif ϕ>=0.99&&ϕ<0.995
                0.0239+0.0079/0.005*(ϕ-0.99)
            elseif ϕ>=0.995&&ϕ<1
                0.0318+0.0153/0.005*(ϕ-0.995)
            else
                println("ERROR! Relative humidity greater than 100 %!")
            end

        #Cement screed
        elseif a==16
            if ϕ>=0&&ϕ<0.35
                82.5/0.35*ϕ
            elseif ϕ>=0.35&&ϕ<0.50
                82.5+25/0.15*(ϕ-0.35)
            elseif ϕ>=0.50&&ϕ<0.7
                107.5+47.5/0.2*(ϕ-0.5)
            elseif ϕ>=0.7&&ϕ<0.8
                155+45/0.1*(ϕ-0.7)
            elseif ϕ>=0.8&&ϕ<0.9
                200+77.5/0.1*(ϕ-0.8)
            elseif ϕ>=0.9&&ϕ<0.95
                277.5+62.5/0.05*(ϕ-0.9)
            elseif ϕ>=0.95&&ϕ<0.96
                340+17.5/0.01*(ϕ-0.95)
            elseif ϕ>=0.96&&ϕ<0.97
                357.5+30/0.01*(ϕ-0.96)
            elseif ϕ>=0.97&&ϕ<0.976
                387.5+22.5/0.006*(ϕ-0.97)
            elseif ϕ>=0.976&&ϕ<1
                410+70/0.024*(ϕ-0.976)
            else
                println("ERROR! Relative humidity greater than 100 %!")
            end

# add new material above! #
    else
        println("ERROR! Material is not in the database!")
    end
end

#----------------------------------------------------------------------#
#    Moisture storage function (water content to relative humidity)    #

# wcRH(m,a=0) -> "m" is water content in kg/m^3 and "a" is a material key #

function wcRH(m,a=0)

    if a==0
        println("ERROR! Please specify material key a in function wcRH(m,a)!")

    #Lime-cement mortar, fine
    elseif a==10
        if m>=0&&m<12.79
            0.3/12.79*m
        elseif m>=12.79&&m<15.88
            0.3+0.2/3.09*(m-12.79)
        elseif m>=15.88&&m<20.98
            0.5+0.2/5.1*(m-15.88)
        elseif m>=20.98&&m<25.66
            0.7+0.1/4.68*(m-20.98)
        elseif m>=25.66&&m<=35.28
            0.8+0.1/9.62*(m-25.66)
        elseif m>=35.28&&m<=46.93
            0.9+0.05/11.65*(m-35.28)
        else
            println("ERROR! Moisture content greater than 46.93 kg/m^3!")
        end

    #Mineral wool
    elseif a==11
        if m>=0&&m<0.461
            0.5/0.461*m
        elseif m>=0.461&&m<0.687
            0.5+0.1/0.226+(m-0.461)
        elseif m>=0.687&&m<1.06
            0.6+0.1/0.373*(m-0.687)
        elseif m>=1.06&&m<1.79
            0.7+0.1/0.73*(m-1.06)
        elseif m>=1.79&&m<2.49
            0.8+0.05/0.7*(m-1.79)
        elseif m>=2.49&&m<3.83
            0.85+0.05/1.34*(m-2.49)
        elseif m>=3.83&&m<4.26
            0.9+0.01/0.43*(m-3.83)
        elseif m>=4.26&&m<4.78
            0.91+0.01/0.52*(m-4.26)
        elseif m>=4.78&&m<5.43
            0.92+0.01/0.65*(m-4.78)
        elseif m>=5.43&&m<6.27
            0.93+0.01/0.84*(m-5.43)
        elseif m>=6.27&&m<7.38
            0.94+0.01/1.11*(m-6.27)
        elseif m>=7.38&&m<8.94
            0.95+0.01/1.56*(m-7.38)
        elseif m>=8.94&&m<11.3
            0.96+0.01/2.36*(m-8.94)
        elseif m>=11.3&&m<15.1
            0.97+0.01/3.8*(m-11.3)
        elseif m>=15.1&&m<22.7
            0.98+0.01/7.6*(m-15.1)
        elseif m>=22.7&&m<30.2
            0.99+0.005/7.5*(m-22.7)
        elseif m>=30.2&&m<=44.8
            0.995+0.005/14.6*(m-30.2)
        else
            println("ERROR! Moisture content greater than 44,8 kg/m^3!")
        end

    #Concrete
    elseif a==12
        if m>=0&&m<27
            0.35/27*m
        elseif m>=27&&m<34
            0.35+0.15/7*(m-27)
        elseif m>=34&&m<49.5
            0.5+0.2/15.5*(m-34)
        elseif m>=49.5&&m<64
            0.7+0.1/14.5*(m-49.5)
        elseif m>=64&&m<88.5
            0.8+0.1/24.5*(m-64)
        elseif m>=88.5&&m<109
            0.9+0.05/20.5*(m-88.5)
        elseif m>=109&&m<114.5
            0.95+0.01/5.5*(m-109)
        elseif m>=114.5&&m<124.5
            0.96+0.1/10*(m-114.5)
        elseif m>=124.5&&m<130
            0.97+0.006/5.5*(m-124.5)
        elseif m>=130&&m<=140
            0.976+0.024/10*(m-130)
        else
            println("ERROR! Moisture content greater than 140 kg/m^3!")
        end

    #Reinforced concrete
    elseif a==13
        if m>=0&&m<27
            0.35/27*m
        elseif m>=27&&m<34
            0.35+0.15/7*(m-27)
        elseif m>=34&&m<49.5
            0.5+0.2/15.5*(m-34)
        elseif m>=49.5&&m<64
            0.7+0.1/14.5*(m-49.5)
        elseif m>=64&&m<88.5
            0.8+0.1/24.5*(m-64)
        elseif m>=88.5&&m<109
            0.9+0.05/20.5*(m-88.5)
        elseif m>=109&&m<114.5
            0.95+0.01/5.5*(m-109)
        elseif m>=114.5&&m<124.5
            0.96+0.1/10*(m-114.5)
        elseif m>=124.5&&m<130
            0.97+0.006/5.5*(m-124.5)
        elseif m>=130&&m<=140
            0.976+0.024/10*(m-130)
        else
            println("ERROR! Moisture content greater than 140 kg/m^3!")
        end

    #Extruded polystyrene
    elseif a==14
        if m>=0&&m<0.009
            0.1/0.009*m
        elseif m>=0.009&&m<0.08
            0.1+0.4/0.071*(m-0.009)
        elseif m>=0.08&&m<0.31
            0.5+0.3/0.23*(m-0.08)
        elseif m>=0.31&&m<0.66
            0.8+0.1/0.35*(m-0.31)
        elseif m>=0.66&&m<1.95
            0.9+0.07/1.29*(m-0.66)
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #PVC hydro membrane
    elseif a==15
        if m>=0&&m<0.000485
            0.5/0.000485*m
        elseif m>=0.000485&&m<0.000724
            0.5+0.1/0.000239*(m-0.000485)
        elseif m>=0.000724&&m<0.00112
            0.6+0.1/0.000396*(m-0.000724)
        elseif m>=0.00112&&m<0.00188
            0.7+0.1/0.00076*(m-0.0011)
        elseif m>=0.00188&&m<0.00262
            0.8+0.05/0.0074*(m-0.00112)
        elseif m>=0.00262&&m<0.00403
            0.85+0.05/0.00141*(m-0.00262)
        elseif m>=0.00403&&m<0.00448
            0.9+0.01/0.00045*(m-0.00403)
        elseif m>=0.00448&&m<0.00503
            0.91+0.01/0.00055*(m-0.00448)
        elseif m>=0.00503&&m<0.00572
            0.92+0.01/0.00069*(m-0.00503)
        elseif m>=0.00572&&m<0.0066
            0.93+0.01/0.00088*(m-0.00572)
        elseif m>=0.0066&&m<0.00777
            0.94+0.01/0.00117*(m-0.0066)
        elseif m>=0.00777&&m<0.00941
            0.95+0.01/0.00164*(m-0.00777)
        elseif m>=0.00941&&m<0.0119
            0.96+0.01/0.00249*(m-0.00941)
        elseif m>=0.0119&&m<0.0159
            0.97+0.01/0.004*(m-0.0119)
        elseif m>=0.0159&&m<0.0239
            0.98+0.01/0.008*(m-0.0159)
        elseif m>=0.0239&&m<0.0318
            0.99+0.005/0.0079*(m-0.0239)
        elseif m>=0.0318&&m<0.0471
            0.995+0.005/0.0153*(m-0.0318)
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Cement screed
    elseif a==16
        if m>=0&&m<82.5
            0.35/82.5*m
        elseif m>=82.5&&m<107.5
            0.35+0.15/25*(m-82.5)
        elseif m>=107.5&&m<155
            0.5+0.2/47.5*(m-107.5)
        elseif m>=155&&m<200
            0.7+0.1/45*(m-155)
        elseif m>=200&&m<277.5
            0.8+0.1/77.5*(m-200)
        elseif m>=277.5&&m<340
            0.9+0.05/62.5*(m-277.5)
        elseif m>=340&&m<357.5
            0.95+0.01/17.5*(m-340)
        elseif m>=357.5&&m<387.5
            0.96+0.01/30*(m-357.5)
        elseif m>=387.5&&m<410
            0.97+0.006/22.5*(m-387.5)
        elseif m>=410&&m<480
            0.976+0.024/70*(m-410)
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

# add new material above! #
    else
        println("ERROR! Material is not in the database!")
    end
end

#---------------------------------------#
#    Moisture dependent conductivity    #

# λ_m(m,a=0) -> "m" is relative humidity [%] and "a" is a material key #

# function λ_m(m,a=0)
#
#     if a==0
#             println("ERROR! Please specify material key a in function λ_m(m,a)!")
#
#     #Lime-cement mortar, fine
#     elseif a==10
#         if m>=0&&RHwc(m,a)<=210
#             0.6+0.69/210*RHwc(m,a)
#         else
#             println("ERROR! Relative humidity greater than 100 %!")
#         end
#
#     #Mineral wool
#     elseif a==11
#         if m>=0&&RHwc(m,a)<10
#             0.04
#         elseif RHwc(m,a)>=10&&RHwc(m,a)<20
#             0.04+0.001/10*(RHwc(m,a)-10)
#         elseif RHwc(m,a)>=20&&RHwc(m,a)<=50
#             0.041+0.002/30*(RHwc(m,a)-20)
#         else
#             println("ERROR! Relative humidity greater than 100 %!")
#         end
#
#     #Concrete
#     elseif a==12
#         if m>=0&&m<=100
#             1.7
#         else
#             println("ERROR! Relative humidity greater than 100 %!")
#         end
#
#
# # add new material above! #
#     else
#         println("ERROR! Material is not in the database!")
#     end
# end

#---------------------------------------------------------------------#
#    Moisture dependent water vapour diffusion resistance factor μ    #

# μ_m(m,a=0) -> "m" is relative humidity [%] and "a" is a material key #

function μ_m(ϕ,a=0)

    if a==0
            println("ERROR! Please specify material key a in function μ_m(ϕ,a)!")

    #Lime-cement mortar, fine
    elseif a==10
        10
        # if ϕ>=0&&ϕ<0.50
        #     147
        # elseif ϕ>=0.50&&ϕ<0.70
        #     147-28/0.20*(ϕ-0.50)
        # elseif ϕ>=0.70&&ϕ<0.80
        #     119-59.5/0.10*(ϕ-0.70)
        # elseif ϕ>=0.80&&ϕ<0.90
        #     59.5-40.3/0.10*(ϕ-0.80)
        # elseif ϕ>=0.90&&ϕ<0.95
        #     19.2-15.03/0.5*(ϕ-0.90)
        # elseif ϕ>=0.95&&ϕ<0.96
        #     4.17-0.6/0.01*(ϕ-0.95)
        # elseif ϕ>=0.96&&ϕ<0.97
        #     3.57-1.65/0.01*(ϕ-0.96)
        # elseif ϕ>=0.97&&ϕ<0.976
        #     1.92-0.67/0.06*(ϕ-0.97)
        # elseif ϕ>=0.976#&&ϕ<=1
        #     1.25
        # else
        #     println("ERROR! Relative humidity greater than 100 %!")
        # end

    #Mineral wool
    elseif a==11
        #if ϕ>=0#&&ϕ<1
            1.2
        # else
        #     println("ERROR! Relative humidity greater than 100 %!")
        # end

    #Concrete
    elseif a==12
        130
        # if ϕ>=0&&ϕ<0.50
        #     50-8.38/0.50*ϕ
        # elseif ϕ>=0.50&&ϕ<0.70
        #     41.62-7.84/0.20*(ϕ-50)
        # elseif ϕ>=0.70&&ϕ<0.80
        #     33.78-8.83/0.10*(ϕ-0.70)
        # elseif ϕ>=0.80&&ϕ<0.90
        #     24.95-14.37/0.10*(ϕ-0.80)
        # elseif ϕ>=0.90&&ϕ<0.95
        #     10.58-7.31/0.5*(ϕ-0.90)
        # elseif ϕ>=0.95#&&ϕ<=1
        #     3.27
        # else
        #     println("ERROR! Relative humidity greater than 100 %!")
        # end

    #Reinforced concrete
    elseif a==13
        130
        # if ϕ>=0&&ϕ<0.50
        #     50-8.38/0.50*ϕ
        # elseif ϕ>=0.50&&ϕ<0.70
        #     41.62-7.84/0.20*(ϕ-0.50)
        # elseif ϕ>=0.70&&ϕ<0.80
        #     33.78-8.83/0.10*(ϕ-0.70)
        # elseif ϕ>=0.80&&ϕ<0.90
        #     24.95-14.37/0.10*(ϕ-0.80)
        # elseif ϕ>=0.90&&ϕ<0.95
        #     10.58-7.31/0.5*(ϕ-0.90)
        # elseif ϕ>=0.95#&&ϕ<=1
        #     3.27
        # else
        #     println("ERROR! Relative humidity greater than 100 %!")
        # end

    #Extruded polystyrene
    elseif a==14
        200
        # if ϕ>=0#&&ϕ<1
        #     170.56
        # end

    #PVC hydro membrane
    elseif a==15
        100000
        # if ϕ>=0#&&ϕ<1
        #     10761
        # end

    #Cement screed
    elseif a==16
        50
        # if ϕ>=0&&ϕ<0.7
        #     25/0.7*ϕ
        # elseif ϕ>=0.7&&ϕ<0.8
        #     25-5/0.1*(ϕ-0.7)
        # elseif ϕ>=0.8&&ϕ<0.9
        #     20-13.75/0.1*(ϕ-0.8)
        # elseif ϕ>=0.9&&ϕ<0.95
        #     6.25-3.12/0.05*(ϕ-0.9)
        # elseif ϕ>=0.95&&ϕ<0.96
        #     3.13-0.62/0.01*(ϕ-0.95)
        # elseif ϕ>=0.96&&ϕ<0.97
        #     2.5-0.71/0.01*(ϕ-0.96)
        # elseif ϕ>=0.97&&ϕ<0.976
        #     1.79-0.32/0.006*(ϕ-0.97)
        # elseif ϕ>=0.976#&&ϕ<1
        #     1.47-0.17/0.024*(ϕ-0.976)
        # end

# add new material above! #
    else
        println("ERROR! Material is not in the database!")
    end
end

### Moisture storage function -- capacity ###

# dw_dϕ(ϕ,a=0)

function dw_dϕ(ϕ,a=0)

    if a==0
        println("ERROR! Please specify ϕaterial key a in function RHwc(ϕ,a)!")

    #Lime-cement mortar, fine
    elseif a==10
        if ϕ>=0&&ϕ<0.30
            12.79/30
        elseif ϕ>=0.30&&ϕ<0.50
            3.09/20
        elseif ϕ>=0.50&&ϕ<0.70
            5.09/20
        elseif ϕ>=0.70&&ϕ<0.80
            4.68/10
        elseif ϕ>=0.80&&ϕ<=0.90
            9.62/10
        elseif ϕ>=0.90#&&ϕ<=1#&&ϕ<=0.95
            11.65/5
        else
            println("ERROR! Relative humidity greater than 95 %!")
        end

    #Mineral wool
    elseif a==11
        if ϕ>=0&&ϕ<0.50
            0.461/50
        elseif ϕ>=0.50&&ϕ<0.60
            0.226/10
        elseif ϕ>=0.60&&ϕ<0.70
            0.373/10
        elseif ϕ>=0.70&&ϕ<0.80
            0.73/10
        elseif ϕ>=0.80&&ϕ<0.85
            0.7/5
        elseif ϕ>=0.85&&ϕ<0.90
            1.34/5
        elseif ϕ>=0.90&&ϕ<0.91
            0.43/1
        elseif ϕ>=0.91&&ϕ<0.92
            0.52/1
        elseif ϕ>=0.92&&ϕ<0.93
            0.65/1
        elseif ϕ>=0.93&&ϕ<0.94
            0.84/1
        elseif ϕ>=0.94&&ϕ<0.95
            1.11/1
        elseif ϕ>=0.95&&ϕ<0.96
            1.56/1
        elseif ϕ>=0.96&&ϕ<0.97
            2.36/1
        elseif ϕ>=0.97&&ϕ<0.98
            3.8/1
        elseif ϕ>=0.98&&ϕ<0.99
            7.6/1
        elseif ϕ>=0.99&&ϕ<0.995
            7.32/0.5
        elseif ϕ>=0.995#&&ϕ<=1
            14.6/0.5
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Concrete
    elseif a==12
        if ϕ>=0&&ϕ<0.35
            27/35
        elseif ϕ>=0.35&&ϕ<0.50
            7/15
        elseif ϕ>=0.50&&ϕ<0.70
            15.5/20
        elseif ϕ>=0.70&&ϕ<0.80
            14.5/10
        elseif ϕ>=0.80&&ϕ<0.90
            24.5/10
        elseif ϕ>=0.90&&ϕ<0.95
            20.5/5
        elseif ϕ>=0.95&&ϕ<0.96
            5.5/1
        elseif ϕ>=0.96&&ϕ<0.97
            10/1
        elseif ϕ>=0.97&&ϕ<0.976
            5.5/0.6
        elseif ϕ>=0.976#&&ϕ<=1
            10/2.4
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Reinforced concrete
    elseif a==13
        if ϕ>=0&&ϕ<0.35
            27/35
        elseif ϕ>=0.35&&ϕ<0.50
            7/15
        elseif ϕ>=0.50&&ϕ<0.70
            15.5/20
        elseif ϕ>=0.70&&ϕ<0.80
            14.5/10
        elseif ϕ>=0.80&&ϕ<0.90
            24.5/10
        elseif ϕ>=0.90&&ϕ<0.95
            20.5/5
        elseif ϕ>=0.95&&ϕ<0.96
            5.5/1
        elseif ϕ>=0.96&&ϕ<0.97
            10/1
        elseif ϕ>=0.97&&ϕ<0.976
            5.5/0.6
        elseif ϕ>=0.976#&&ϕ<=1
            10/2.4
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Extruded polystyrene
    elseif a==14
        if ϕ>=0&&ϕ<0.10
            0.009/10
        elseif ϕ>=0.10&&ϕ<0.50
            0.071/40
        elseif ϕ>=0.50&&ϕ<0.80
            0.23/30
        elseif ϕ>=0.80&&ϕ<0.90
            0.35/10
        elseif ϕ>=0.90#&&ϕ<=1#ϕ<0.97
            1.29/7
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #PVC hydro membrane
    elseif a==15
        if ϕ>=0&&ϕ<0.5
            0.000485/50
        elseif ϕ>=0.5&&ϕ<0.6
            0.000239/10
        elseif ϕ>=0.60&&ϕ<0.70
            0.000396/10
        elseif ϕ>=0.70&&ϕ<0.80
            0.00076/10
        elseif ϕ>=0.80&&ϕ<0.85
            0.0074/5
        elseif ϕ>=0.85&&ϕ<0.9
            0.00141/5
        elseif ϕ>=0.9&&ϕ<0.91
            0.00045/1
        elseif ϕ>=0.91&&ϕ<0.92
            0.00055/1
        elseif ϕ>=0.92&&ϕ<0.93
            0.00069/1
        elseif ϕ>=0.93&&ϕ<0.94
            0.00088/1
        elseif ϕ>=0.94&&ϕ<0.95
            0.00117/1
        elseif ϕ>=0.95&&ϕ<0.96
            0.00164/1
        elseif ϕ>=0.96&&ϕ<0.97
            0.00249/1
        elseif ϕ>=0.97&&ϕ<0.98
            0.004/1
        elseif ϕ>=0.98&&ϕ<0.99
            0.008/1
        elseif ϕ>=0.99&&ϕ<0.995
            0.0079/0.5
        elseif ϕ>=0.995#&&ϕ<1
            0.0153/0.5
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

    #Cement screed
    elseif a==16
        if ϕ>=0&&ϕ<0.35
            82.5/35
        elseif ϕ>=0.35&&ϕ<0.50
            25/15
        elseif ϕ>=0.50&&ϕ<0.7
            47.5/20
        elseif ϕ>=0.7&&ϕ<0.8
            45/10
        elseif ϕ>=0.8&&ϕ<0.9
            77.5/10
        elseif ϕ>=0.9&&ϕ<0.95
            62.5/5
        elseif ϕ>=0.95&&ϕ<0.96
            17.5/1
        elseif ϕ>=0.96&&ϕ<0.97
            30/1
        elseif ϕ>=0.97&&ϕ<0.976
            22.5/0.6
        elseif ϕ>=0.976#&&ϕ<1
            70/2.4
        else
            println("ERROR! Relative humidity greater than 100 %!")
        end

# add new material above! #
    else
        println("ERROR! Material is not in the database!")
    end
end


text:
"""Serializes a InfrastructureSystemsType to a JSON file."""
function to_json(obj::T, filename::AbstractString) where {T <: InfrastructureSystemsType}
    result = open(filename, "w") do io
        return to_json(io, obj)
    end

    @info "Serialized $T to $filename"
    return result
end

"""Serializes a InfrastructureSystemsType to a JSON string."""
function to_json(obj::T)::String where {T <: InfrastructureSystemsType}
    return JSON2.write(obj)
end

"""JSON Serializes a InfrastructureSystemsType to an IO stream in JSON."""
function to_json(io::IO, obj::T) where {T <: InfrastructureSystemsType}
    return JSON2.write(io, obj)
end

"""Deserializes a InfrastructureSystemsType from a JSON filename."""
function from_json(::Type{T}, filename::String) where {T <: InfrastructureSystemsType}
    return open(filename) do io
        from_json(io, T)
    end
end

"""Deserializes a InfrastructureSystemsType from String or IO."""
function from_json(io::Union{IO, String}, ::Type{T}) where {T <: InfrastructureSystemsType}
    return JSON2.read(io, T)
end

"""Enables JSON deserialization of TimeSeries.TimeArray.
The default implementation fails because the data field is defined as an AbstractArray.
Deserialization can't determine the actual concrete type.
"""
function JSON2.read(io::IO, ::Type{T}) where {T <: TimeSeries.TimeArray}
    data = JSON2.read(io)
    timestamp = [Dates.DateTime(x) for x in data.timestamp]
    colnames = [Symbol(x) for x in data.colnames]
    dim2 = length(colnames)
    dim1 = Int(length(data.values) / dim2)

    for i in eachindex(data.values)
        data.values[i] = Float64(data.values[i])
    end

    if length(colnames) > 1
        vals = reshape(data.values, dim1, dim2)
    else
        vals = data.values
    end

    return TimeSeries.TimeArray(timestamp, vals, colnames)
end

"""Enables JSON deserialization of Dates.Period.
The default implementation fails because the field is defined as abstract.
Encode the type when serializing so that the correct value can be deserialized.
"""
function JSON2.write(resolution::Dates.Period)
    return JSON2.write(encode_for_json(resolution))
end

function JSON2.write(io::IO, resolution::Dates.Period)
    return JSON2.write(io, encode_for_json(resolution))
end

function encode_for_json(resolution::Dates.Period)
    return (value = resolution.value, unit = strip_module_name(string(typeof(resolution))))
end

function JSON2.read(io::IO, ::Type{T}) where {T <: Dates.Period}
    data = JSON2.read(io)
    return getfield(Dates, Symbol(data.unit))(data.value)
end

"""
The next few methods fix serialization of UUIDs. The underlying type of a UUID is a UInt128.
JSON2 tries to encode this as a number in JSON. Encoding integers greater than can
be stored in a signed 64-bit integer sometimes does not work - at least when using the
JSON2.@pretty option. The number gets converted to a float in scientific notation, and so
the UUID is truncated and essentially lost. These functions cause JSON2 to encode UUIDs as
strings and then convert them back during deserialization.
"""

function JSON2.write(uuid::Base.UUID)
    return JSON2.write(encode_for_json(uuid))
end

function JSON2.write(io::IO, uuid::Base.UUID)
    return JSON2.write(io, encode_for_json(uuid))
end

function JSON2.read(io::IO, ::Type{Base.UUID})
    data = JSON2.read(io)
    return Base.UUID(data.value)
end

function encode_for_json(uuid::Base.UUID)
    return (value = string(uuid),)
end

"""The next set of methods fix serialization for Complex numbers."""

function JSON2.write(value::Complex)
    return JSON2.write(encode_for_json(value))
end

function JSON2.write(io::IO, value::Complex)
    return JSON2.write(io, encode_for_json(value))
end

function JSON2.read(io::IO, ::Type{Complex})
    data = JSON2.read(io)
    return Complex(data.real, data.imag)
end

function encode_for_json(value::Complex)
    return (real = real(value), imag = imag(value))
end

# Refer to docstrings in services.jl.

function JSON2.write(io::IO, forecast::Forecast)
    return JSON2.write(io, encode_for_json(forecast))
end

function JSON2.write(forecast::Forecast)
    return JSON2.write(encode_for_json(forecast))
end

function encode_for_json(forecast::T) where {T <: Forecast}
    fields = [x for x in fieldnames(T) if x != :data]
    vals = []

    for name in fields
        if name == :data
            # The timeseries is stored within SystemForecasts.
            continue
        end
        val = getfield(forecast, name)
        if val isa InfrastructureSystemsType
            push!(vals, get_uuid(val))
        else
            push!(vals, val)
        end
    end

    push!(fields, :type)
    push!(vals, strip_module_name(string(T.name)))
    return NamedTuple{Tuple(fields)}(vals)
end


text:
function OEMWriter()
    return OEMWriter(())
end

function OEMWriter(arg0::OEMWriter_InterpolationMethod, arg1::JString, arg2::JString, arg3::JString)
    return OEMWriter((OEMWriter_InterpolationMethod, JString, JString, JString), arg0, arg1, arg2, arg3)
end

function write(obj::EphemerisFileWriter, arg0::JString, arg1::EphemerisFile)
    return jcall(obj, "write", void, (JString, EphemerisFile), arg0, arg1)
end

function write(obj::OEMWriter, arg0::Appendable, arg1::EphemerisFile)
    return jcall(obj, "write", void, (Appendable, EphemerisFile), arg0, arg1)
end


text:
"""
    read data, sample ID;s and observation ID's to DataFrame.
"""
function readCooccurrence(source::String, sort_on=nothing; abs_cutoff=nothing, rel_cutoff=nothing, uncommon_observation=-1)::DataFrame
    #= parse parameters =#
    if sort_on === nothing || startswith(lowercase(sort_on), "sample")
        x = "sample"
        y = "observation"
    elseif startswith(lowercase(sort_on), "observation")
        x = "observation"
        y = "sample"
    else
        throw(ArgumentError("Can't sort on $(sort_on). Choose sample [default] or observation"))
    end
    if abs_cutoff !== nothing && rel_cutoff !== nothing
        throw(ArgumentError("Can't have both a relative and absolute cutoff."))
    end
    if rel_cutoff !== nothing && (0 > rel_cutoff || rel_cutoff > 1)
        throw(ArgumentError("rel_cutoff is $(rel_cutoff), but should `nothing` or be between 0 and 1."))
    end
    #= read HDF5 file =#
    h5file = HDF5.h5open(source, "r") do h5
        read(h5)
    end
    #= set preprocess parameters =#
    n_obs = length(h5file["sample"]["matrix"]["indptr"])-1
    if abs_cutoff !== nothing && (0 > abs_cutoff || abs_cutoff > n_obs)
        throw(ArgumentError("abs_cutoff is $(abs_cutoff), but should be `nothing` or between 0 and $(n_obs)."))
    end
    if abs_cutoff !== nothing
        preprocess = true
        min_obs = abs_cutoff
    elseif rel_cutoff !== nothing
        preprocess = true
        min_obs = n_obs*rel_cutoff
    else
        preprocess = false
    end
    #= write data to dataframe =#
    indptr = h5file[x]["matrix"]["indptr"]
    # x_col  = [parse(Int, h5file[y]["ids"][i+1]) for i in h5file[x]["matrix"]["indices"]]
    x_col  = [i for i in h5file[x]["matrix"]["indices"]]
    y_col  = [Int32(i) for a  in 1:length(indptr)-1 for i in fill(a, indptr[a+1]-indptr[a])]
    data   = h5file[x]["matrix"]["data"]
    if x == "sample"
        df = DataFrame(data = data, observation = x_col, sample = y_col)
    else
        df = DataFrame(data = data, observation = y_col, sample = x_col)
    end

    # df.observation = parse.(Int, df.observation)
    #= preprocess rare/uncommon observations =#
    c = countmap(df.observation)
    if preprocess
        for i in 1:length(df.observation)
            if c[df.observation[i]] < min_obs
                df.observation[i] = uncommon_observation
            end
        end
    end
    return df
end

function readBIOM(path::String)::Dict{String, Any}
    dict = h5open(path, "r") do h5
        read(h5)
    end
end

function isBIOM(path::String)::Bool
    di = readBIOM(path)
    isBIOM(di)
end

function isBIOM(di::Dict{String, Any})::Bool
    k₁ = keys(di)
    if  length(k₁) != 2 && ["sample", "observation"] ⊆ k₁
        return false
    else
        k₂ₒ = keys(di["observation"])
        k₂ₛ = keys(di["sample"])
        if  length(k₂ₒ) != 4 && ["ids", "matrix", "metadata", "group-metadata"] ⊆ k₂ₒ &&
            length(k₂ₛ) != 4 && ["ids", "matrix", "metadata", "group-metadata"] ⊆ k₂ₛ
            return false
        else
            k₃ₒ = keys(di["observation"]["matrix"])
            k₃ₛ = keys(di["sample"]["matrix"])
            if  length(k₃ₒ) != 3 && ["data", "indptr", "indices"] ⊆ k₃ₒ &&
                length(k₃ₛ) != 3 && ["data", "indptr", "indices"] ⊆ k₃ₛ
                return false
            else
                return true
            end
        end
    end

end

function writeBIOM(path::String, df::DataFrame;
                    sample_meta::Union{Nothing,Dict{String, <:Any}}       = nothing,
                    sample_group_meta::Union{Nothing,Dict{String, <:Any}} = nothing,
                    obs_meta::Union{Nothing,Dict{String, <:Any}}          = nothing,
                    obs_group_meta::Union{Nothing,Dict{String, <:Any}}    = nothing)
    if  ["sample", "observation", "data"] ⊈ names(df)
        Throw(ArgumentError("Something bad happend."))
    end
    dfₒ   = sort(df, :observation)
    dfₛ   = sort(df, :sample)
	dataₒ = Array{Float64}(dfₒ.data)
	samₒ  = Array{Int32}(dfₒ.sample)
	obsₒ  = Array{Int32}([[findfirst(==(i), dfₒ.observation)-1 for i in unique(dfₒ.observation)]... , length(dfₒ.observation)])
	idₒ   = Array{String}(string.(sort(unique(dfₒ.observation))))
    dataₛ = Array{Float64}(dfₛ.data)
	obsₛ  = Array{Int32}(dfₛ.observation)
	samₛ  = Array{Int32}([[findfirst(==(i), dfₛ.sample)-1 for i in unique(dfₛ.sample)]... , length(dfₛ.sample)])
	idₛ   = Array{String}(string.(sort(unique(dfₛ.sample))))
    HDF5.h5write(path, "sample/ids", idₛ)
    HDF5.h5write(path, "sample/matrix/data", dataₛ)
    HDF5.h5write(path, "sample/matrix/indptr", samₛ)
    HDF5.h5write(path, "sample/matrix/indices",obsₛ)
    HDF5.h5write(path, "observation/ids", idₒ)
    HDF5.h5write(path, "observation/matrix/data", dataₒ)
    HDF5.h5write(path, "observation/matrix/indptr", obsₒ)
    HDF5.h5write(path, "observation/matrix/indices", samₒ)
    if sample_meta !== nothing
        for k in keys(sample_meta)
            HDF5.h5write(path, "sample/metadata/$(k)", sample_meta[k])
        end
    end
    if sample_group_meta !== nothing
        for k in keys(sample_group_meta)
            HDF5.h5write(path, "sample/group-metadata/$(k)", sample_group_meta[k])
        end
    end
    if obs_meta !== nothing
        for k in keys(obs_meta)
            HDF5.h5write(path, "observation/metadata/$(k)", obs_meta[k])
        end
    end
    if obs_group_meta !== nothing
        for k in keys(obs_group_meta)
            HDF5.h5write(path, "sample/group-metadata/$(k)", obs_group_meta[k])
        end
    end
end

function collapseBIOM(biomin::String, biomout::String, field::String; on::Int)
    if isfile(biomout)
        throw(ArgumentError("$(biomout) already exists. Can't overwrite files."))
    end
    if biomin == biomout
        throw(ArgumentError("Can't change/override/append to biom file."))
    end
    if field != "taxonomy"
        throw(ErrorException("Collapsing on $(field) is not implemented (yet)."))
    end
    indict = readBIOM(biomin)
    if !haskey(indict["observation"]["metadata"], field)
        throw(KeyError(field))
    end
    if size(indict["observation"]["metadata"][field], 1) < on
        throw(BoundsError())
    end
    df = readCooccurrence(biomin)
    tax = indict["observation"]["metadata"]["taxonomy"]
    dftax =  transform!(df, :observation => (x -> tax[on, x]) => :taxonomy)
    remove = Array{Int}([])
    sort!(dftax, [:taxonomy, :observation])
    for i in 2:size(dftax,1)
        if length(dftax[i, :taxonomy]) > 3 && dftax[i,:taxonomy] == dftax[i-1, :taxonomy]
            dftax[i,:observation] = dftax[i-1,:observation]
        end
    end
    sort!(dftax, [:sample, :observation])
    for i in 2:size(dftax,1)
        if dftax[i,:observation] == dftax[i-1,:observation] &&
            dftax[i, :sample] == dftax[i-1, :sample] &&
            length(dftax[i, :taxonomy]) > 3
            dftax[i,:data] += dftax[i-1,:data]
            push!(remove, i-1)
        end
    end
    delete!(dftax, remove)

    sort!(dftax, :observation)
    col_tax = [dftax[1, :taxonomy]]
    for i in 2:size(dftax,1)
        if dftax[i,:observation] != dftax[i-1,:observation]
            push!(col_tax, dftax[i,:taxonomy])
        end
    end
    outdict = writeBIOM(biomout, select(dftax, :sample, :data, :observation), sample_meta=Dict("taxonomy" => col_tax))
end

# function preprocessBIOM!(inpath::String, outpath::String, rel_cutoff::Float64; uncommon_observation=-1)
#     @warn "This function is not very well tested yet. Please don't overwrite files that shouldn't be potentially corrupted"
#     if rel_cutoff !== nothing && (0 > rel_cutoff || rel_cutoff > 1)
#         throw(ArgumentError("rel_cutoff is $(rel_cutoff), but should `nothing` or be between 0 and 1."))
#     end
#     #= read HDF5 file =#
#     h5in = HDF5.h5open(inpath, "r") do h5
#         read(h5)
#     end
#     #= set preprocess parameters =#
#     n_obs = length(h5in["sample"]["matrix"]["indptr"])-1
#     min_obs = n_obs*rel_cutoff
#     #= add synthetic observation =#
#     push!(hdin["observation"]["ids"], "synthetic")
#     syn_id = length(hdin["observation"]["ids"])
#     #= change observation under cut off to synthetic =#
#     observations = hdin["sample"]["matrix"]["indices"]
#     obs_count = countmap(observations)
#     removed_ids = Set()
#     if preprocess
#         for i in 1:length(observations)
#             if obs_count[i] < min_obs
#                 push!(observations[i])
#                 observations[i] = syn_id
#             end
#         end
#     end
#     #= adjust ids index =#
#     obs_indptr = hdin["observation"]["matrix"]["indptr"]
#     updated_obs_indptr = [i for i in obs_indptr if !contains(removed_ids, i)]

#     #= sort data for hdf5 =#
#     sample_indptr = h5in["sample"]["matrix"]["indptr"]
#     df = DataFrame(observation = observations,
#                    data        = h5in["sample"]["matrix"]["data"],
#                    sample      = [i-1 for a  in 1:length(sample_indptr)-1 for i in fill(a, sample_indptr[a+1]-sample_indptr[a])] )
#     sort!(df, :observation)
#     h5in["observation"]["matrix"]["data"] = df.data
#     h5in["observation"]["matrix"]["indices"] = df.sample

# end


text:
# # Random particles in a circle

# The code [particles_in_circle.jl](particles_in_circle.jl) compares the scattered wave from one big circle, with the scattered wave from a circle filled with small particles.

# ```julia
using MultipleScattering

#You can also pick your own shape, an generate random particles inside it
#with a certain radius ands volume fraction
radius = 0.3
volfrac = 0.45
centre = [0.,0.]
big_radius = 3.0

particle_medium = Acoustic(2; ρ=0.0, c=0.0) # 2D particle with density ρ = 0.0 and soundspeed c = 0.0
particle_shape = Circle(radius)

circle = Sphere(centre, big_radius)

particles = random_particles(particle_medium, particle_shape; region_shape = circle, volume_fraction = volfrac, seed=1)

x = [-10.,0.] # position to receive the reflected wave
host_medium = Acoustic(2; ρ=1.0, c=1.0)
source =  plane_source(host_medium; position = x, direction = [1.0,0.])
simulation = FrequencySimulation(particles, source)
# ```
# The particles chosen are impenetrable, i.e. the wave is 100\% reflected. So this circle filled with scatterers should act like one big particle.
# ```julia
big_particle = Particle(particle_medium, circle)
big_particle_simulation = FrequencySimulation([big_particle], source)

#define a bounding box for plot
    bottomleft = [-10, -2*big_radius]
    topright = [big_radius, 2*big_radius]
    box = Box([bottomleft, topright])

using Plots
height = 300
#gr(size=(1.4*height,height))
# pyplot(leg=false, size=(1.4*height,height))

ω = 0.5
#plot(big_particle_simulation, ω; res=15, bounds = box);
#plot!(big_particle)

#savefig("plot_field_big.png")

#plot(simulation, ω; res=15, bounds = box);
#plot!(particles, linecolor = :green)

#savefig("plot_field.png")
# ```
# Resulting in the figures:

# ![The field with big particle](plot_field_big.png)
# ![The field with particles](plot_field.png)

# If we compare the response measured at the listener `[-10., 0.]`, they should be very similar:
# ```julia
#define angular frequency range
ωs = collect(LinRange(0.1,1.0,10))
result = run(simulation, x, ωs)
big_result = run(big_particle_simulation, x, ωs)

#plot(result, lab = "scattering from particles")
#plot!(big_result,
#    lab = "scattering from big particle",
#    title="Compare scattered wave from one big particle, \n and a circle filled with small particles")
# ```
# ![The response comparison](plot_response_compare.png)


text:
i = 0
while true
    print(i,"\n")
    global i += 1
end


text:
module GriddingMachine

# include submodules
include("Collections.jl")
include("Indexer.jl")
include("Requestor.jl")
include("GriddingMachineOld.jl")

end # module


text:
function dftReg(imgRef::AbstractArray{Complex{T},N},imgF::AbstractArray{Complex{T},N},usfac) where {T, N}
    dftReg(CPU1(),imgRef,imgF,usfac)
end
"
`dftReg{N}(imgRef::AbstractArray{Complex{Float64},N},imgF::AbstractArray{Complex{Float64},N},usfac)`

Main internal function : takes the Fourier transforms of the arrays to register as inputs (`imgRef`/`imgF`) and returns a dictionary containing the residual error and the shift along the dimensions of the arrays, with the level of subpixel precision provided by `usfac`"
function dftReg(resource::CPU1,imgRef::AbstractArray{Complex{T},N},imgF::AbstractArray{Complex{T},N},usfac) where {T, N}
    if usfac==0
        ## Compute error for no pixel shift
        CCmax = sum(imgRef.*conj(imgF))
        rfzero = sum(abs2, imgRef)
        rgzero = sum(abs2, imgF)
        error = 1 - CCmax*conj(CCmax)/(rgzero*rfzero)
        diffphase = atan(imag(CCmax),real(CCmax))
        output = Dict("error" => error)
    elseif usfac==1
        ## Whole-pixel shift - Compute crosscorrelation by an IFFT and locate the peak
        L = length(imgRef)
        CC = ifft(imgRef.*conj(imgF))
        loc = argmax(abs.(CC))
        CCmax=CC[loc]
        rfzero = sum(abs2, imgRef)/L
        rgzero = sum(abs2, imgF)/L
        error = abs(1 - CCmax*conj(CCmax)/(rgzero*rfzero))
        diffphase = atan(imag(CCmax),real(CCmax))

        indi = size(imgRef)
        ind2 = tuple([div(x,2) for x in indi]...)

        #locI = [ind2sub(indi,loc)...]
        locI = [Tuple(loc)...]

        shift = zeros(size(locI))
        for i in eachindex(locI)
            if locI[i]>ind2[i]
                shift[i]=locI[i]-indi[i]-1
            else shift[i]=locI[i]-1
            end
        end

        output = Dict("error"=>error,"shift"=>shift,"diffphase"=>diffphase)
    else
        ## Partial pixel shift

        ##First upsample by a factor of 2 to obtain initial estimate
        ##Embed Fourier data in a 2x larger array
        dimR = [size(imgRef)...]
        ranges = [(x+1-div(x,2)):(x+1+div(x-1,2)) for x in dimR]
        dimRL = map(x -> x*2, dimR)
        CC = zeros(Complex{Float64},tuple(dimRL...))
        CC[ranges...] = fftshift(imgRef).*conj(fftshift(imgF))

        ##  Compute crosscorrelation and locate the peak
        CC = ifft(ifftshift(CC))
        loc = argmax(abs.(CC))

        indi = size(CC)
        #locI = [ind2sub(indi,loc)...]
        locI = [Tuple(loc)...]
        CCmax = CC[loc]
        ## Obtain shift in original pixel grid from the position of the crosscorrelation peak

        ind2 = tuple([div(x,2) for x in indi]...)

        shift = zeros(size(locI))
        for i in eachindex(locI)
            if locI[i]>ind2[i]
                shift[i]=locI[i]-indi[i]-1
            else shift[i]=locI[i]-1
            end
        end
        shift = shift/2

        ## If upsampling > 2, then refine estimate with matrix multiply DFT
        if usfac > 2
            ### DFT Computation ###
            # Initial shift estimate in upsampled grid
            shift = round.(Integer,shift*usfac)/usfac
            dftShift = div(ceil(usfac*1.5),2) ## center of output array at dftshift+1
            ## Matrix multiplies DFT around the current shift estimate
            CC = conj(dftups(imgF.*conj(imgRef),ceil(Integer,usfac*1.5),usfac,dftShift.-shift*usfac))/(prod(ind2)*usfac^N)
            ## Locate maximum and map back to original pixel grid
            loc = argmax(abs.(CC))
            # locI = ind2sub(size(CC),loc)
            locI = Tuple(loc)
            CCmax = CC[loc]
            rg00 = dftups(imgRef.*conj(imgRef),1,usfac)[1]/(prod(ind2)*usfac^N)
            rf00 = dftups(imgF.*conj(imgF),1,usfac)[1]/(prod(ind2)*usfac^N)
            locI = map((x) -> x - dftShift - 1,locI)

            for i in eachindex(shift)
                shift[i]=shift[i]+locI[i]/usfac
            end

        else
            rg00 = sum(imgRef.*conj(imgRef))/(prod(indi))
            rf00 = sum(imgF.*conj(imgF))/(prod(indi))
        end
        error = 1 - CCmax*conj(CCmax)/(rg00*rf00)
        error = sqrt(abs.(error))
        diffphase = atan(imag(CCmax),real(CCmax))
        ## If its only one row or column the shift along that dimension has no effect. Set to zero.
        shift[[div(x,2) for x in size(imgRef)].==1].=0

        output = Dict("error"=>error,"shift"=>shift,"diffphase"=>diffphase)

    end
    output
end

"
`dftups{T,N}(inp::AbstractArray{T,N},no,usfac::Int=1,off=zeros(N))`

Upsampled DFT by matrix multiplication, computes an upsampled DFT in just a small region. `no` is the size of the region in pixels, `offset` the position in the full array, `usfac` the upsampling parameter."
function dftups(inp::AbstractArray{T,N},no,usfac::Int=1,offset=zeros(N)) where {T,N}
    sz = [size(inp)...]
    permV = 1:N
    for i in permV
        inp = permutedims(inp,[i;deleteat!(collect(permV),i)])
        kern = exp.((-1im*2*pi/(sz[i]*usfac))*((0:(no-1)).-offset[i])*transpose(ifftshift(0:(sz[i]-1)).-floor(sz[i]/2)))
        d = size(inp)[2:N]
        inp = kern * reshape(inp, Val(2))
        inp = reshape(inp,(no,d...))
    end
    permutedims(inp,collect(ndims(inp):-1:1))
end

"
`subPixShift(imgft::AbstractArray{Complex{Float64}},shift::Array{Float64,1})`
    Shift the image `imgft` (in Fourier space) by the amount provided in the vector `shift`."
function subPixShift(imgft::AbstractArray{Complex{T},N},shift::Array{Float64,2},diffphase::Array{Float64,1}) where {T,N}
    sz = [size(imgft)[1:(N-1)]...]
    nim = size(imgft)[N]
    Z=0
    resh = [repeat([1],inner=N-1)...,nim]
    for i in eachindex(sz)
        shifti = ifftshift((-div(sz[i],2)):(ceil(Integer,sz[i]/2)-1)).*shift[:,i]'/sz[i]
        reshTemp = copy(resh)
        reshTemp[i] = size(shifti)[1]
        Z = Z .- reshape(shifti,reshTemp...)
    end

    Greg = imgft .* exp.(2im*pi*Z)
    Greg = Greg .* reshape(exp.(1im*diffphase),resh...)
    Greg
end

function subPixShift(imgft::AbstractArray{Complex{T}},shift::Array{Float64,1},diffphase) where {T}
    sz = [size(imgft)...]
    N=0
    for i in eachindex(sz)
        shifti = ifftshift((-div(sz[i],2)):(ceil(Integer,sz[i]/2)-1))*shift[i]/sz[i]
        resh = (repeat([1],inner=[i-1])...,length(shifti))
        N = N .- reshape(shifti,resh)
    end

    Greg = imgft .* exp.(2im*pi*N)
    Greg = Greg .* exp(1im*diffphase)
    Greg
end

function stackDftReg(imgser::AbstractArray{T,N};ref::AbstractArray{T,N1}=reshape(selectdim(imgser,N,1),size(imgser)[1:(N-1)]),ufac::Int=10) where {T,N,N1}
    stackDftReg(CPU1(),imgser,ref=ref,ufac=ufac)
end
"
`stackDftReg{T,N,N1}(imgser::AbstractArray{T,N};ref::AbstractArray{T,N1}=reshape(slicedim(imgser,N,1),size(imgser)[1:(N-1)]),ufac::Int=10)`

`dftReg` applied to a full array. Each array along the last dimension of `imgser` is aligned to `ref` (by default the first image of the series, with precision `ufac`. Returns an array of `Dict` containing the translation information."
function stackDftReg(resource::CPU1,imgser::AbstractArray{T,N};ref::AbstractArray{T,N1}=reshape(selectdim(imgser,N,1),size(imgser)[1:(N-1)]),ufac::Int=10) where {T,N,N1}
    if ((N1 != (N - 1)) & (N1 != N))
        error("Reference image has the wrong dimensionality")
    end
    ref = fft(ref)
    imgF = fft(imgser,(1:N1...,))

    if N1 == (N-1)
        imgF = [reshape(selectdim(imgF,N,i),size(ref)) for i = 1:size(imgF)[N]]
        results = pmap(imgF) do im
            dftReg(ref,im,ufac)
        end
    else
        results = dftReg(ref,imgF,ufac)
    end
    return results
end

"
`alignFromDft{T,N}(img2reg::AbstractArray{T,N},dftRegRes::Array{Any,1})`

Given an array and a `Dict` of translations as returned by `dftReg`, returns the aligned array."
function alignFromDict(img2reg::AbstractArray{T,N},dftRegRes::Array{Dict{String,Any},1}) where {T,N}
    alignFromDict(CPU1(),img2reg,dftRegRes)
end

function alignFromDict(resource::CPU1,img2reg::AbstractArray{T,N},dftRegRes::Array{Dict{String,Any},1}) where {T,N}
    if (length(dftRegRes) != size(img2reg)[N])
        error("Alignment results and image stack dimensionalities don't match.")
    end

    img2regF = fft(img2reg,(1:(N-1)...,))
    szF = size(img2reg)[1:(N-1)]

    imRes = subPixShift(img2regF,shifts_from_res(dftRegRes),diffphase_from_res(dftRegRes))
    imRes = real(ifft(imRes,1:(N-1)))

    imRes
end

## Only a single Dict, means we expect one image of the same size
function alignFromDict(img2reg::AbstractArray{T,N},dftRegRes::Dict) where {T,N}
    if (length(dftRegRes["shift"]) != N)
        error("Alignment results and image dimensionalities don't match.")
    end
    frameft = real(ifft(subPixShift(fft(img2reg),dftRegRes["shift"],dftRegRes["diffphase"])))
end

function shifts_from_res(dftRegRes::Array{Dict{String,Any},1})
    shifts = zeros(length(dftRegRes),length(dftRegRes[1]["shift"]))
    for i in eachindex(dftRegRes)
        shifts[i,:] = dftRegRes[i]["shift"]
    end
    shifts
end

function diffphase_from_res(dftRegRes::Array{Dict{String,Any},1})
    dP = zeros(length(dftRegRes))
    for i in eachindex(dftRegRes)
        dP[i] = dftRegRes[i]["diffphase"]
    end
    dP
end


text:
function onpatch!(x::AbstractVector, v̄::PVector{T, N}) where {T,N}
    ranges = ProjectiveVectors.dimension_indices(v̄)
    for range in ranges
        λ = zero(eltype(x))
        @inbounds for i in range
            λ += v̄[i] * x[i]
        end
        λ⁻¹ = @fastmath inv(λ)
        for i in range
            x[i] *= λ⁻¹
        end
    end
    x
end

function evaluate_patch!(u, v̄::PVector{S, N}, x::PVector{T,N}) where {S, T, N}
    ranges = ProjectiveVectors.dimension_indices(v̄)
    n = length(u) - N
    for (k, range) in enumerate(ranges)
        out = -one(eltype(x))
        for i in range
            out += v̄[i] * x[i]
        end
        u[n + k] = out
    end
    nothing
end

function jacobian_patch!(U, v̄::PVector{S, N}, x::PVector) where {S, T, N}
    ranges = ProjectiveVectors.dimension_indices(v̄)
    n = size(U, 1) - N
    for j in 1:size(U, 2), i in (n+1):size(U, 1)
        U[i, j] = zero(eltype(U))
    end
    for (k, range) in enumerate(ranges)
        for j in range
            U[n + k, j] = v̄[j]
        end
    end
    nothing
end


text:
using Test
using Quadmath

for T in (Float64, Int32, Int64)
    @test Float128(T(1)) + Float128(T(2)) == Float128(T(3))
    @test Float128(T(1)) + Float128(T(2)) <= Float128(T(3))
    @test Float128(T(1)) + Float128(T(2)) != Float128(T(4))
    @test Float128(T(1)) + Float128(T(2)) < Float128(T(4))
    @test T(Float128(T(1)) + Float128(T(2))) === T(3)
end

@test Base.exponent_one(Float128) == reinterpret(UInt128, Float128(1.0))


text:
using OrdinaryDiffEq
using Trixi
using Plots

######################Functions that need to be loaded ###########
# mapping O-mesh
function mapping(xi_, eta_)
    ξ = xi_
    η = eta_
    x = 2 * (2 + ξ ) * cos(π * (η + 1))
    y = 2 * (2 + ξ ) * sin(π * (η + 1))
    return SVector(x, y)
end
function mapping1zu1(xi_, eta_)
    x = xi_
    y = eta_
    return SVector(x, y)
end
function WR2_initial_condition_constant(x, t, equations::Trixi.AbstractEquations)
    rho = 1.0
    rho_v1 = 0.1
    rho_v2 = -0.2
    rho_e = 10.0
    return SVector(rho, rho_v1, rho_v2, rho_e)
end
function boundary_condition_constant( u_inner, orientation, direction, x, t,
    surface_flux_function,
    equations::Trixi.AbstractEquations)
    # Far Field Conditions
    u_boundary = initial_condition(x , t, equations)

    # Calculate boundary flux
    if direction in (2, 4) # u_inner is "left" of boundary, u_boundary is "right" of boundary
        flux = surface_flux_function(u_inner, u_boundary, orientation, equations)
    else # direction == 4 # u_boundary is "left" of boundary, u_inner is "right" of boundary
        flux = surface_flux_function(u_boundary, u_inner, orientation, equations)
    end
    return flux
  end

########################    Setting    ######################################
CFL = 0.01
tspan = (0.0, 0.5)
N = 3
c = 16
mu = 0.001

# semidiscretization of the compressible Euler equations

equations = CompressibleEulerEquations2D(1.4, viscous = true, mu = mu)

initial_condition = WR2_initial_condition_constant
boundary_conditions = boundary_condition_constant # boundary_condition_periodic
cells_per_dimension = (c, c)
# mesh = CurvedMesh(cells_per_dimension, mapping, periodicity = true)
mesh = CurvedMesh(cells_per_dimension, (-1.0, -1.0), (1.0, 1.0))

surface_flux = FluxPlusDissipation(flux_chandrashekar, DissipationLocalLaxFriedrichs(max_abs_speed_naive))
basis = LobattoLegendreBasis(N)
volume_flux  = flux_chandrashekar
volume_integral = VolumeIntegralFluxDifferencing(volume_flux)
solver = DGSEM(basis, surface_flux, volume_integral)
semi = SemidiscretizationHyperbolic(mesh, equations, initial_condition, solver)#,boundary_conditions=boundary_conditions)

###############################################################################
# ODE solvers, callbacks etc.

ode = semidiscretize(semi, tspan)

summary_callback = SummaryCallback()

analysis_interval = 100
analysis_callback = AnalysisCallback(semi, interval=analysis_interval, save_analysis=true,
                                      # extra_analysis_errors=(:conservation_error,),
                                      # extra_analysis_integrals=(entropy, energy_total,
                                      # energy_kinetic, energy_internal)
                                      )

alive_callback = AliveCallback(analysis_interval=analysis_interval)

save_solution = SaveSolutionCallback(interval=100,
                                     save_initial_solution=true,
                                     save_final_solution=true,
                                     solution_variables=cons2prim)

stepsize_callback = StepsizeCallback(cfl=CFL)

callbacks = CallbackSet(summary_callback,
                        analysis_callback,
                        #alive_callback,
                        #save_solution,
                        stepsize_callback)

###############################################################################
# run the simulation

sol = solve(ode, CarpenterKennedy2N54(williamson_condition=false),
            dt=1.0, # solve needs some value here but it will be overwritten by the stepsize_callback
            save_everystep=false, callback=callbacks);
summary_callback() # print the timer summary

plot(sol)

max_stars_repo_path: src/WR2/MS5 viscous Freestream.jl

max_issues_repo_path: src/WR2/MS5 viscous Freestream.jl

max_forks_repo_path: src/WR2/MS5 viscous Freestream.jl


text:
missed = Set();
const parameterswithnumval = union(trivalwithnumval, Set(["P_TRIPLE","pcrit","p_triple","MOLEMASS","rhomolar_critical",
"MOLARMASS","Pcrit","gas_constant","Tcrit","p_reducing","Tmin","rhocrit","molarmass","T_CRITICAL",
"P_min","pmax","T_min","T_reducing","T_MAX","P_MIN","T_triple","P_MAX","ptriple","Ttriple",
"acentric","p_critical","Tmax","T_TRIPLE","rhomolar_reducing","MOLAR_MASS",
"P_CRITICAL","T_max","molemass","T_MIN","rhomass_critical","T_critical",
"P_max","RHOMASS_CRITICAL","molar_mass","pmin"]));
global counter = 0;
global longunits = Set();
for p in coolpropparameters
    global longunit = get_parameter_information_string(p, "long") * " | " * get_parameter_information_string(p, "units");
    note = "";
    if (!in(longunit, longunits))
        push!(longunits, longunit);
    else
        note = " *Duplicated* "
    end
    for fluid in coolpropfluids
        try
            res = ("$(PropsSI(p, fluid))");
            (!in(p, parameterswithnumval)) && push!(missed, p);
            note *= " **Constant Property** "
            global counter+=1;
        break;
        catch err
        end
    end
end
length(missed) > 0 && @warn "missed parameters with numerical value:\n $missed)";
@test length(parameterswithnumval) == counter
@info "parameters with numerical values are as expected";


text:
using GeometricBase
using GeometricIntegrators.Solutions
using GeometricIntegrators.Solutions: save#load,
using GeometricProblems.HarmonicOscillator
using Test
import HDF5

ode = harmonic_oscillator_ode()
dae = harmonic_oscillator_dae()
pode = harmonic_oscillator_pode()
pdae = harmonic_oscillator_pdae()

nt = 10
Δt = 0.1
ni = 5

t0 = 0.0
x0 = rand(2)
y0 = rand(2)
q0 = rand(1)
p0 = rand(1)
z0 = rand(3)
λ0 = rand(1)
μ0 = rand(2)

t1 = 1.0
x1 = [rand(2) for j = 1:ni]
y1 = [rand(2) for j = 1:ni]
q1 = [rand(1) for j = 1:ni]
p1 = [rand(1) for j = 1:ni]
z1 = [rand(3) for j = 1:ni]
λ1 = [rand(1) for j = 1:ni]
μ1 = [rand(2) for j = 1:ni]

x100 = rand(ndims(ode), 100)
y100 = rand(ndims(ode), 100)
z100 = rand(ndims(dae), 100)
q100 = rand(ndims(pode), 100)
p100 = rand(ndims(pode), 100)
λ100 = rand(dae.m, 100)
μ100 = rand(pdae.m, 100)

h5file = "test.hdf5"

@testset "$(rpad("HDF5 I/O for ODE Solution",80))" begin
    # test general hdf5 functions
    sol1 = Solution(ode, Δt, nt)
    h5io = SolutionHDF5(h5file, sol1; overwrite = true)
    @test typeof(hdf5(h5io)) == HDF5.File
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)

    sol1 = Solution(ode, Δt, nt)
    h5io = SolutionHDF5(h5file; overwrite = false)
    @test typeof(hdf5(h5io)) == HDF5.File
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)

    rm(h5file)

    # test hdf5 in- and output
    sol1 = Solution(harmonic_oscillator_ode(x0), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    # h5io = SolutionHDF5(h5file) # TODO
    sol2 = SolutionODE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_ode(x1), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionODE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(ode, Δt, 100, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            set_solution!(sol1, x100[:, (i-1)*10+j], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol = 1E-14
    @test sol2.q.d == vcat(ode.q₀, [x100[:, j] for j in axes(x100, 2)])
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(ode, Δt, 100, nsave = 2, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            set_solution!(sol1, x100[:, (i-1)*10+j], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol = 1E-14
    @test sol2.q.d == vcat(ode.q₀, [x100[:, j] for j in axes(x100, 2)[2:2:end]])
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)

end

@testset "$(rpad("HDF5 I/O for PODE Solution",80))" begin
    sol1 = Solution(harmonic_oscillator_pode(q0, p0), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionPODE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_pode(q1, p1), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionPODE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(pode, Δt, 100, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            oj = (i - 1) * 10 + j
            set_solution!(sol1, q100[:, oj], p100[:, oj], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionPODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol = 1E-14
    @test sol2.q.d == vcat(pode.q₀, [q100[:, j] for j in axes(q100, 2)])
    @test sol2.p.d == vcat(pode.p₀, [p100[:, j] for j in axes(p100, 2)])
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(pode, Δt, 100, nsave = 2, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            oj = (i - 1) * 10 + j
            set_solution!(sol1, q100[:, oj], p100[:, oj], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionPODE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol = 1E-14
    @test sol2.q.d == vcat(pode.q₀, [q100[:, j] for j in axes(q100, 2)[2:2:end]])
    @test sol2.p.d == vcat(pode.p₀, [p100[:, j] for j in axes(p100, 2)[2:2:end]])
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("HDF5 I/O for DAE Solution",80))" begin
    sol1 = Solution(harmonic_oscillator_dae(z0, λ0), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionDAE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_dae(z1, λ1), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionDAE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(dae, Δt, 100, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            oj = (i - 1) * 10 + j
            set_solution!(sol1, z100[:, oj], λ100[:, oj], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol = 1E-14
    @test sol2.q.d == vcat(dae.q₀, [z100[:, j] for j in axes(z100, 2)])
    @test sol2.λ.d == vcat(dae.λ₀, [λ100[:, j] for j in axes(λ100, 2)])
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(dae, Δt, 100, nsave = 2, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            oj = (i - 1) * 10 + j
            set_solution!(sol1, z100[:, oj], λ100[:, oj], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol = 1E-14
    @test sol2.q.d == vcat(dae.q₀, [z100[:, j] for j in axes(z100, 2)[2:2:end]])
    @test sol2.λ.d == vcat(dae.λ₀, [λ100[:, j] for j in axes(λ100, 2)[2:2:end]])
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end

@testset "$(rpad("HDF5 I/O for PDAE Solution",80))" begin
    sol1 = Solution(harmonic_oscillator_pdae(x0, y0, μ0), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionPDAE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(harmonic_oscillator_pdae(x1, y1, μ1), Δt, nt)
    h5io = SolutionHDF5(h5file, sol1)
    save(h5io, sol1)
    close(h5io)
    @test isfile(h5file)
    sol2 = SolutionPDAE(h5file)
    @test sol1 != sol2
    @test sol1.t == sol2.t
    @test sol1.q == sol2.q
    @test sol1.p == sol2.p
    @test sol1.λ == sol2.λ
    @test sol1.ntime == sol2.ntime
    @test sol1.nsave == sol2.nsave
    rm(h5file)

    sol1 = Solution(pdae, Δt, 100, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            oj = (i - 1) * 10 + j
            set_solution!(sol1, x100[:, oj], y100[:, oj], μ100[:, oj], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionPDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:100) atol = 1E-14
    @test sol2.q.d == vcat(pdae.q₀, [x100[:, j] for j in axes(x100, 2)])
    @test sol2.p.d == vcat(pdae.p₀, [y100[:, j] for j in axes(y100, 2)])
    @test sol2.λ.d == vcat(pdae.λ₀, [μ100[:, j] for j in axes(μ100, 2)])
    # @test sol2.q.d == hcat(reshape(pdae.q₀, (ndims(pdae),1)), x100)
    # @test sol2.p.d == hcat(reshape(pdae.p₀, (ndims(pdae),1)), y100)
    # @test sol2.λ.d == hcat(reshape(pdae.λ₀, (pdae.m,1)), μ100)
    @test sol2.ntime == 100
    @test sol2.nsave == 1
    rm(h5file)

    sol1 = Solution(pdae, Δt, 100, nsave = 2, nwrite = 10)
    h5io = SolutionHDF5(h5file, sol1)
    for i = 1:10
        for j = 1:10
            oj = (i - 1) * 10 + j
            set_solution!(sol1, x100[:, oj], y100[:, oj], μ100[:, oj], j)
        end
        save(h5io, sol1)
        reset!(sol1)
    end
    close(h5io)
    sol2 = SolutionPDAE(h5file)
    @test sol2.t.t ≈ Δt .* collect(0:2:100) atol = 1E-14
    @test sol2.q.d == vcat(pdae.q₀, [x100[:, j] for j in axes(x100, 2)[2:2:end]])
    @test sol2.p.d == vcat(pdae.p₀, [y100[:, j] for j in axes(y100, 2)[2:2:end]])
    @test sol2.λ.d == vcat(pdae.λ₀, [μ100[:, j] for j in axes(μ100, 2)[2:2:end]])
    # @test sol2.q.d == hcat(reshape(pdae.q₀, (ndims(pdae),1)), x100)[:,1:2:end]
    # @test sol2.p.d == hcat(reshape(pdae.p₀, (ndims(pdae),1)), y100)[:,1:2:end]
    # @test sol2.λ.d == hcat(reshape(pdae.λ₀, (pdae.m,1)), μ100)[:,1:2:end]
    @test sol2.ntime == 100
    @test sol2.nsave == 2
    rm(h5file)
end


text:
import DataStructures:update!
export AbstractActionSelector

"""
    AbstractActionSelector

A subtype of `AbstractActionSelector` is used to generate an action
given the estimated valud of different actions.

| Required Methods| Brief Description |
|:----------------|:------------------|
| `selector(values)` | `selector`, an instance of `AbstractActionSelector`, must be a callable object which takes in an estimation and returns an action |
"""
abstract type AbstractActionSelector end

function update!(as::AbstractActionSelector, s, a) end


text:
barbell = QNetwork()
A = BasicNode("A")
B = BasicNode("B")
AB = BasicChannel(A, B)
AB.costs = Dict("loss"=>1, "Z"=>0.5)
for i in [A, B, AB]
    add(barbell, i)
end
refresh_graph!(barbell)

max_stars_repo_name: FalafelGood/QuNet

max_issues_repo_name: FalafelGood/QuNet

max_forks_repo_name: FalafelGood/QuNet


text:
module FastRounding

export add_round, sub_round, mul_round, square_round,
       inv_round, div_round, sqrt_round

using ErrorfreeArithmetic

const SysFloat = Union{Float64, Float32}  # fast iff fma is available in hardware

setrounding(Float64, RoundNearest)
setrounding(Float32, RoundNearest)

@inline function add_round(a::T, b::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_sum(a, b)
    return round_errorfree(hi, lo, rounding)
end
add_round(a::T, b::T) where {T<:SysFloat} = a + b

@inline function sub_round(a::T, b::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_diff(a, b)
    return round_errorfree(hi, lo, rounding)
end
sub_round(a::T, b::T) where {T<:SysFloat} = a - b

@inline function mul_round(a::T, b::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_prod(a, b)
    return round_errorfree(hi, lo, rounding)
end
mul_round(a::T, b::T) where {T<:SysFloat} = a * b

@inline function inv_round(a::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_inv(a)
    return round_errorfree(hi, lo, rounding)
end
inv_round(a::T) where {T<:SysFloat} = inv(a)

@inline function div_round(a::T, b::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_div(a, b)
    return round_errorfree(hi, lo, rounding)
end
div_round(a::T, b::T) where {T<:SysFloat} = a / b

@inline function square_round(a::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_square(a)
    return round_errorfree(hi, lo, rounding)
end
square_round(a::T) where {T<:SysFloat} = a^2

@inline function sqrt_round(a::T, rounding::R)::T where {T<:SysFloat, R<:RoundingMode}
    hi, lo = two_sqrt(a)
    return round_errorfree(hi, lo, rounding)
end
sqrt_round(a::T) where {T<:SysFloat} = sqrt(a)

#=
    To perform arithmetic with directed rounding more rapidly
      we use error-free transformations to control rounding.
=#

@inline function round_errorfree(hi::T, lo::T, ::RoundingMode{:Nearest})::T where {T<:SysFloat}
     !isinf(hi) && return hi
     return signbit(hi) ? T(-Inf) : T(Inf)
end

@inline function round_errorfree(hi::T, lo::T, ::RoundingMode{:ToZero})::T where {T<:SysFloat}
     !isinf(hi) && return signbit(hi) != signbit(lo) ? next_nearerto_zero(hi) : hi
     return signbit(hi) ? nextfloat(T(-Inf)) : prevfloat(T(Inf))
end

@inline function round_errorfree(hi::T, lo::T, ::RoundingMode{:FromZero})::T where {T<:SysFloat}
    !isinf(hi) && return signbit(hi) == signbit(lo) ? next_awayfrom_zero(hi) : hi
    return signbit(hi) ? T(-Inf) : T(Inf)
end

@inline function round_errorfree(hi::T, lo::T, ::RoundingMode{:Up})::T where {T<:SysFloat}
    !isinf(hi) && return (lo<=zero(T) || isnan(lo))  ? hi : nextfloat(hi)
    return signbit(hi) ? nextfloat(T(-Inf)) : T(Inf)
end

@inline function round_errorfree(hi::T, lo::T, ::RoundingMode{:Down})::T where {T<:SysFloat}
    !isinf(hi) && return (lo>=zero(T) || isnan(lo))  ? hi : prevfloat(hi)
    return signbit(hi) ? T(-Inf) : prevfloat(T(Inf))
end

@inline next_nearerto_zero(x::T) where {T<:SysFloat} = !signbit(x) ? prevfloat(x) : nextfloat(x)
@inline next_awayfrom_zero(x::T) where {T<:SysFloat} = !signbit(x) ? nextfloat(x) : prevfloat(x)

end # module


text:
struct QStructs{FT<:Real}
	A  :: FT
	Lx :: FT
	Qx :: FT
	Ly :: FT
	Qy :: FT
end

function QfieldProperties(
	FT = Float64;
	A  :: Real = 1.,
	Lx :: Real = 2.,
	Qx :: Real = 0.,
	Ly :: Real = 2.,
	Qy :: Real = 0.
)

	return QStructs{FT}(A,Lx,Qx,Ly,Qy)

end

function createQ(
	QParams::Array{QStructs{FT}},
	g::Grid{FT},
) where FT <: Real

	xc = g.xc
	yc = g.yc
	nQ = length(QParams)
	Q  = zeros(FT,g.nx,g.ny)

	for Qs in QParams, jj = 1 : g.ny, ii = 1 : g.nx
		if abs(xc[ii]-Qs.Qx) < Qs.Lx
			Q[ii,jj] += Qs.A * cos(pi*(xc[ii]-Qs.Qx)*0.5/Qs.Lx) *
							   exp(-((yc[jj]-Qs.Qy)/Qs.Ly)^2)
		end
	end

	return Q

end

function show(io::IO, Q::QStructs{FT}) where FT <: Real
    print(
		io,
		"The Q{$FT}-field forcing is defined as follows:\n",
		"           Amplitude : ", Q.A, '\n',
		"        Size (Lx,Ly) : ", (Q.Lx,Q.Ly), '\n',
		"    Position (Qx,Qy) : ", (Q.Qx,Q.Qy), '\n'
	)
end


text:
using FirstOrderSolvers: conjugategradient!, CGdata

Random.seed!(2)
# a = randn(1000,1000)
# A = a'a
A = rand(1000,1000)
A = A'*A
b = randn(1000)
x0 = randn(1000)

x = copy(x0)
#anrm = sqrt(sum(abs2, A, 2))[:]
#A = A./anrm
#b = b./anrm
function test(A, b, x0, x, cgdata)
    x .= x0
    conjugategradient!(x, A, b, cgdata.r, cgdata.p, cgdata.z)
end

cgdata = CGdata(similar(b), similar(b), similar(b))
conjugategradient!(x, A, b, cgdata.r, cgdata.p, cgdata.z, max_iters = 100)
n0 = norm(A*x-b)    # 48
conjugategradient!(x, A, b, cgdata.r, cgdata.p, cgdata.z, max_iters = 5000)
n1 = norm(A*x-b)    # 3e-6

@test n1 < 1e-5

xcopy = x .+ 1e-5.*randn(size(x))
n2 = norm(A*xcopy-b) #0.9
conjugategradient!(xcopy, A, b, cgdata.r, cgdata.p, cgdata.z, max_iters = 100)
n3 = norm(A*xcopy-b) #3e-5

@test n3 < 10*n2


text:
function correct_regions!(wrapped, visited, regions)
    for r in regions
        wrapped[visited .== r] .-= (2π * median(round.(wrapped[visited .== r] ./ 2π)))
    end
end

function merge_regions!(wrapped, visited, nregions, weights)
    mask = sum(weights; dims=1)
    region_size = countmap(visited) # TODO could use weight instead
    offsets = zeros(nregions, nregions)
    offset_counts = zeros(Int, nregions, nregions)
    stridelist = getdimoffsets(wrapped)
    for dim in 1:3
        neighbor = CartesianIndex(ntuple(i->i==dim ? 1 : 0, 3))
        for I in CartesianIndices(wrapped)
            J = I + neighbor
            if checkbounds(Bool, wrapped, J)# && weights[dim, I] < 255
                ri = visited[I]
                rj = visited[J]
                if ri != 0 && rj != 0 && ri != rj
                    w = 255 - weights[dim,I]
                    if w == 255 w = 0 end
                    offsets[ri, rj] += (wrapped[I] - wrapped[J]) * w
                    offset_counts[ri, rj] += w
                end
            end
        end
    end
    for i in 1:nregions, j in i:nregions
        offset_counts[i,j] += offset_counts[j,i]
        offset_counts[j,i] = offset_counts[i,j]
        offsets[i,j] -= offsets[j,i]
        offsets[j,i] = -offsets[i,j]
    end
    corrected = falses(nregions)
    remaining_regions = Int[]
    while !all(corrected)
        largest_uncorrected_region = try
            findmax(filter(p -> first(p) != 0 && !corrected[first(p)], region_size))[2]
        catch
            @show region_size size(corrected) nregions
            throw(error())
        end
        # TODO correct region? No, offsets are already calculated
        corrected[largest_uncorrected_region] = true
        push!(remaining_regions, largest_uncorrected_region)

        # TODO multiple rounds until no change?
        sorted_offsets = get_offset_count_sorted(offset_counts, corrected)
        for I in sorted_offsets
            (i,j) = Tuple(I)
            offset = round((offsets[I] / offset_counts[I]) / 2π)
            if offset != 0
                wrapped[visited .== j] .+= offset * 2π
                visited[visited .== j] .= i
                # TODO region merging offset_count and offset calculation
            end
            corrected[j] = true
            offset_counts[i,j] = offset_counts[j,i] = -1
        end
    end
    return remaining_regions
end

function get_offset_count_sorted(offset_counts, corrected)
    f(I) = corrected[I[1]] && !corrected[I[2]]
    return sort(filter(f, findall(offset_counts .> 0)), by=i->offset_counts[i], rev=true)
end


text:
export select_region_ft, select_region_rft
export rft_fix_after, rft_fix_before, ft_fix_after, ft_fix_before

"""
    select_region_ft(mat,new_size)

Performs the necessary Fourier-space operations of resampling
in the space of ft (meaning the already circshifted version of fft).

Note that in dependence of the size, fixes in Fourier-space are applied under the assumption
that the corresponding data in Real space is real-valued only. You should used this function
if you want to select a region (extract) in Fourier space and know that you are dealing with
a real-valued signal. For complex-valued data, the ordinary `select_region` function should be
used instead.
`new_size`.

The size of the array view after the operation finished. The Fourier-center
is always assumed to align before and after the padding operation.

 # Examples
```jldoctest
julia> x = [1 20 3; 4 500 6; -7 821 923]
3×3 Matrix{Int64}:
  1   20    3
  4  500    6
 -7  821  923

julia> ffts(x)
3×3 ShiftedArrays.CircShiftedArray{ComplexF64, 2, Matrix{ComplexF64}}:
   106.5+390.577im  -1099.5-1062.61im   1000.5+700.615im
 -1138.5+354.204im   2271.0+0.0im      -1138.5-354.204im
  1000.5-700.615im  -1099.5+1062.61im    106.5-390.577im

julia> select_region_ft(ffts(x), (4,4))
4×4 PaddedView(0.0 + 0.0im, OffsetArray(::ShiftedArrays.CircShiftedArray{ComplexF64, 2, Matrix{ComplexF64}}, 2:4, 2:4), (Base.OneTo(4), Base.OneTo(4))) with eltype ComplexF64:
 0.0+0.0im      0.0+0.0im          0.0+0.0im          0.0+0.0im
 0.0+0.0im    106.5+390.577im  -1099.5-1062.61im   1000.5+700.615im
 0.0+0.0im  -1138.5+354.204im   2271.0+0.0im      -1138.5-354.204im
 0.0+0.0im   1000.5-700.615im  -1099.5+1062.61im    106.5-390.577im

julia> x = [1 20; 4 500; -7 821; -2 2]
4×2 Matrix{Int64}:
  1   20
  4  500
 -7  821
 -2    2

julia> ffts(x)
4×2 ShiftedArrays.CircShiftedArray{ComplexF64, 2, Matrix{ComplexF64}}:
  -347.0+0.0im     331.0+0.0im
   809.0-492.0im  -793.0+504.0im
 -1347.0+0.0im    1339.0+0.0im
   809.0+492.0im  -793.0-504.0im

julia> select_region_ft(ffts(x), (5,3))
5×3 FourierTools.FourierSplit{ComplexF64, 2, FourierTools.FourierSplit{ComplexF64, 2, PaddedViews.PaddedView{ComplexF64, 2, Tuple{Base.OneTo{Int64}, Base.OneTo{Int64}}, OffsetArrays.OffsetMatrix{ComplexF64, ShiftedArrays.CircShiftedArray{ComplexF64, 2, Matrix{ComplexF64}}}}}}:
 -86.75+0.0im     165.5+0.0im    -86.75+0.0im
  404.5-246.0im  -793.0+504.0im   404.5-246.0im
 -673.5+0.0im    1339.0+0.0im    -673.5+0.0im
  404.5+246.0im  -793.0-504.0im   404.5+246.0im
 -86.75+0.0im     165.5+0.0im    -86.75+0.0im
```
"""
function select_region_ft(mat,new_size)
    old_size = size(mat)
    mat_fixed_before = ft_fix_before(mat,old_size,new_size)
    mat_pad = ft_pad(mat_fixed_before,new_size)
    # afterwards we add the highest pos. frequency to the highest lowest one
    return ft_fix_after(mat_pad ,old_size,new_size)
end

"""
    select_region_rft(mat,old_size, new_size)

performs the necessary Fourier-space operations of resampling
in the space of rft (meaning the already circshifted version of rfft).
Note that in dependence of the size, fixes in Fourier-space are applied, which is
why you should used this function if you want to select a region (extract) in rft space.
Since rfts always assume the corresponding Real-space data to be real, this rule should
have no exception for rft data.

`old_size`.
The size of the corresponding real-space array before it was rfted.

`new_size`.
The size of the corresponding real-space array view after the operation finished. The Fourier-center
is always assumed to align before and after the padding aperation.
"""
function select_region_rft(mat,old_size, new_size)
    rft_old_size = size(mat)
    rft_new_size = Base.setindex(new_size,new_size[1]÷2 +1, 1)
    return rft_fix_after(rft_pad(
        rft_fix_before(mat,old_size,new_size),
        rft_new_size),old_size,new_size)
end

"""
    select_region(mat,new_size)

performs the necessary Fourier-space operations of resampling
in the space of ft (meaning the already circshifted version of fft).

`new_size`.
The size of the array view after the operation finished.

`center`.
Specifies the center of the new view in coordinates of the old view. By default an alignment of the Fourier-centers is assumed.
# Examples
```jldoctest
julia> using FFTW, FourierTools

julia> select_region(ones(3,3),new_size=(7,7),center=(1,3))
7×7 PaddedView(0.0, OffsetArray(::Matrix{Float64}, 4:6, 2:4), (Base.OneTo(7), Base.OneTo(7))) with eltype Float64:
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  1.0  1.0  1.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0  0.0  0.0
```
"""
function select_region(mat; new_size=size(mat), center=ft_center_diff(size(mat)).+1, pad_value=zero(eltype(mat)))
    new_size = Tuple(expand_size(new_size, size(mat)))
    center = Tuple(expand_size(center, ft_center_diff(size(mat)).+1))
    oldcenter = ft_center_diff(new_size).+1
    PaddedView(pad_value, mat,new_size, oldcenter .- center.+1);
end

function ft_pad(mat, new_size)
    return select_region(mat;new_size=new_size)
end

function rft_pad(mat, new_size)
    c2 = rft_center_diff(size(mat))
    c2 = Base.setindex(c2, new_size[1].÷2, 1);
    return select_region(mat;new_size=new_size, center=c2.+1)
end

function ft_fix_before(mat, size_old, size_new; start_dim=1)
    for d = start_dim:ndims(mat)
        sn = size_new[d]
        so = size_old[d]
        if sn < so && iseven(sn)
            L1 = (size_old[d] -size_new[d] )÷2 +1
            mat = FourierJoin(mat, d, L1)
        end
    end
    return mat
end

function ft_fix_after(mat,size_old,size_new; start_dim=1)
    start_dim
    ndims(mat)
    for d=start_dim:ndims(mat)
        sn = size_new[d]
        so = size_old[d]
        if sn > so && iseven(so)
            L1 = (size_new[d]-size_old[d])÷2+1
            mat = FourierSplit(mat,d,L1)
        end
        # if equal do nothing
    end
    return mat
end

function rft_fix_first_dim_before(mat,size_old,size_new;dim=1)
    sn = size_new[dim] # Note that this dim is the corresponding real-space size
    so = size_old[dim] # Note that this dim is the corresponding real-space size
    if sn < so && iseven(sn) # result size is even upon cropping
        L1 = size_new[dim] ÷ 2 + 1
        mat = FourierJoin(mat, dim, L1, L1) # a hack to dublicate the value
    end
    return mat
end

function rft_fix_first_dim_after(mat,size_old,size_new;dim=1)
    sn = size_new[dim] # Note that this dim is the corresponding real-space size
    so = size_old[dim] # Note that this dim is the corresponding real-space size
    if sn > so && iseven(so) # source size is even upon padding
        L1 = size_old[dim] ÷ 2 + 1
        mat = FourierSplit(mat,dim,L1,-1) # This hack prevents a second position to be affected
    end
    # if equal do nothing
    return mat
end

function rft_fix_before(mat,size_old,size_new)
    mat=rft_fix_first_dim_before(mat,size_old,size_new;dim=1) # ignore the first dimension
    ft_fix_before(mat,size_old,size_new;start_dim=2) # ignore the first dimension
end

function rft_fix_after(mat,size_old,size_new)
    mat = rft_fix_first_dim_after(mat,size_old,size_new;dim=1) # ignore the first dimension
    ft_fix_after(mat,size_old,size_new;start_dim=2) # ignore the first dimension
end


text:
# Toggles for nGraph environmental variables
enable_codegen() = ENV["NGRAPH_CODEGEN"] = 1
disable_codegen() = delete!(ENV, "NGRAPH_CODEGEN")

"""
    codegen_debug()

Set environmental variables to allow internal codegeneration to display debug information
to `stdout`.
"""
function codegen_debug()
    # The actual values that are assigned don't actually matter. It just matters that these
    # variables are assigned to begin with.
    ENV["NGRAPH_COMPILER_DEBUGINFO_ENABLE"] = true
    ENV["NGRAPH_COMPILER_DIAG_ENABLE"] = true
    ENV["NGRAPH_COMPILER_REPORT_ENABLE"] = true
end

enable_timing() = ENV["NGRAPH_CPU_TRACING"] = true
disable_timing() = delete!(ENV, "NGRAPH_CPU_TRACING")

# Pass Enables
abstract type AbstractPassAttribute end
struct ReuseMemory <: AbstractPassAttribute end
express(::ReuseMemory) = "ReuseMemory=1"

set_pass_attributes(x::AbstractPassAttribute...) = (ENV["NGRAPH_PASS_ATTRIBUTES"] = join(express.(x), ";"))