text:
# Tests for Radial DEA Models
@testset "AdditiveDEAModel" begin

    ## Test Additive DEA Models with FLS Book data
    # Test against results with R
    X = [5 13; 16 12; 16 26; 17 15; 18 14; 23 6; 25 10; 27 22; 37 14; 42 25; 5 17]
    Y = [12; 14; 25; 26; 8; 9; 27; 30; 31; 26; 12]

    # Additive CRS
    deaaddcrs1 = deaadd(X, Y, :Custom, wX = ones(size(X)), wY = ones(size(Y)), rts = :CRS)

    @test nobs(deaaddcrs1) == 11
    @test ninputs(deaaddcrs1) == 2
    @test noutputs(deaaddcrs1) == 1
    @test efficiency(deaaddcrs1) ≈ [0.0000000000;
                               10.76923077;
                               10.83783784;
                               0.0000000000;
                               22.15384615;
                               17.92307692;
                               0.0000000000;
                               12.07692308;
                               11.61379310;
                               35.000000000;
                               4.000000000]

    @test convert(Matrix, peers(deaaddcrs1)) ≈
    [ 1.0000000000  0  0 0.0000000000  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.5384615385  0  0 0.0000000000  0  0   0   0
      0.1216216216  0  0 0.9054054054  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 1.0000000000  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.3076923077  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.3461538462  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.0000000000  0  0 1.0000000000  0  0   0   0
      0.0000000000  0  0 1.1538461538  0  0 0.0000000000  0  0   0   0
      0.0000000000  0  0 0.4689655172  0  0 0.6965517241  0  0   0   0
      0.0000000000  0  0 1.0000000000  0  0 0.0000000000  0  0   0   0
      1.0000000000  0  0 0.0000000000  0  0 0.0000000000  0  0   0   0]

    # Additive VRS model
    deaaddvrs1 = deaadd(X, Y, :Custom, wX = ones(size(X)), wY = ones(size(Y)), rts = :VRS)

    @test nobs(deaaddvrs1) == 11
    @test ninputs(deaaddvrs1) == 2
    @test noutputs(deaaddvrs1) == 1
    @test efficiency(deaaddvrs1) ≈ [0.0000000000;
                               7.333333333;
                               0.0000000000;
                               0.0000000000;
                               18.000000000;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               35.000000000;
                               4.000000000]

    @test convert(Matrix, peers(deaaddvrs1)) ≈
    [1.000000000  0  0  0  0  0 0.00000000000  0  0   0   0;
     0.6666666667 0  0  0  0  0 0.3333333333   0  0   0   0;
     0.000000000  0  1  0  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0  1  0  0 0.00000000000  0  0   0   0;
     1.000000000  0  0  0  0  0 0.00000000000  0  0   0   0;
     0.000000000  0  0  0  0  1 0.00000000000  0  0   0   0;
     0.000000000  0  0  0  0  0 1.00000000000  0  0   0   0;
     0.000000000  0  0  0  0  0 0.00000000000  1  0   0   0;
     0.000000000  0  0  0  0  0 0.00000000000  0  1   0   0;
     0.000000000  0  0  1  0  0 0.00000000000  0  0   0   0;
     1.000000000  0  0  0  0  0 0.00000000000  0  0   0   0]

    # Check weights equal to one equals default model
    deaadddefaultcrs = deaadd(X, Y, rts = :CRS)
    @test efficiency(deaaddcrs1) ≈ efficiency(deaadddefaultcrs)

    deaadddefaultvrs = deaadd(X, Y)
    @test efficiency(deaaddvrs1) ≈ efficiency(deaadddefaultvrs)

    # Test model :Ones equals model with all weights equal to 1
    deaaddonescrs = deaadd(X, Y, :Ones, rts = :CRS)
    @test deaaddonescrs.weights == :Ones
    @test efficiency(deaaddonescrs) ≈ efficiency(deaadddefaultcrs)

    deaaddonesvrs = deaadd(X, Y, :Ones)
    @test deaaddonesvrs.weights == :Ones
    @test efficiency(deaaddonesvrs) ≈ efficiency(deaadddefaultvrs)

    # MIP CRS
    deaaddmipcrs = deaadd(X, Y, :MIP, rts = :CRS)
    @test deaaddmipcrs.weights == :MIP
    @test efficiency(deaaddmipcrs) ≈ [0.0000000000;
                               0.7577160494;
                               0.4168399168;
                               0.0000000000;
                               2.2219512195;
                               1.1478260870;
                               0.0000000000;
                               0.4867909868;
                               0.4041184041;
                               1.0726153846;
                               0.2352941176]
    @test deaaddmipcrs.slackX ≈ [0 0;
                                 3.037037037 6.814814815;
                                 0 10.837837838;
                                 0 0;
                                 0 0;
                                 8.000000000  0;
                                 0 0;
                                 7.384615385 4.692307692;
                                 8.296296296 2.518518519;
                                 0.000000000  8.200000000;
                                 0 4.000000000]

    @test deaaddmipcrs.slackY ≈ [0;
                                  0;
                                  0;
                                  0;
                                  17.77560976 ;
                                  7.20000000;
                                  0;
                                  0;
                                  0;
                                  19.36000000;
                                  0]

    # MIP VRS
    deaaddmipvrs = deaadd(X, Y, :MIP, rts = :VRS)
    @test deaaddmipvrs.weights == :MIP
    @test efficiency(deaaddmipvrs) ≈ [0.0000000000;
                               0.5075187970;
                               0.0000000000;
                               0.0000000000;
                               2.2039473684;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               0.0000000000;
                               1.0432234432;
                               0.2352941176]
    @test deaaddmipvrs.slackX ≈ [0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 0 0;
                                 17 15;
                                 0 4]
    @test deaaddmipvrs.slackY ≈ [0;
                                 7.105263158;
                                 0;
                                 0;
                                 17.631578947;
                                 0;
                                 0;
                                 0;
                                 0;
                                 1.000000000;
                                 0]

    # Normalized CRS
    deaaddnormalizedcrs = deaadd(X, Y, :Normalized, rts = :CRS)
    @test deaaddnormalizedcrs.weights == :Normalized
    @test efficiency(deaaddnormalizedcrs) ≈ [0.0000000000;
                                 1.3569256615;
                                 1.7407259078;
                                 0.0000000000;
                                 2.6877810368;
                                 1.6953153107;
                                 0.0000000000;
                                 1.6540884810;
                                 1.1212042237;
                                 4.0172855145;
                                 0.6424624298]

    @test deaaddnormalizedcrs.slackX ≈ [0.000000000  0.000000000;
                                 3.037037037  6.814814815;
                                 0.000000000 10.837837838;
                                 0.000000000  0.000000000;
                                 10.592592593 11.037037037;
                                 14.666666667  2.666666667;
                                 0.000000000  0.000000000;
                                 0.000000000 10.298429319;
                                 8.296296296  2.518518519;
                                 17.925925926 15.370370370;
                                 0.000000000  4.000000000]

    @test deaaddnormalizedcrs.slackY ≈ [0;
                                  0;
                                  0;
                                  0;
                                  0 ;
                                  0;
                                  0;
                                  0;
                                  0;
                                  0;
                                  0]

    # Normalized VRS
    deaaddnormalizedvrs = deaadd(X, Y, :Normalized, rts = :VRS)
    @test deaaddnormalizedvrs.weights == :Normalized
    @test efficiency(deaaddnormalizedvrs) ≈ [0.0000000000;
                                 0.8049248943;
                                 0.0000000000;
                                 0.0000000000;
                                 2.0149704908;
                                 0.0000000000;
                                 0.0000000000;
                                 0.0000000000;
                                 0.0000000000;
                                 3.9898943952;
                                 0.6424624298]

    @test deaaddnormalizedvrs.slackX ≈ [0  0.00;
                                 0  0.65;
                                 0  0.00;
                                 0  0.00;
                                 0  2.95;
                                 0  0.00;
                                 0  0.00;
                                 0  0.00;
                                 0  0.00;
                                 17 15.00;
                                 0  4.00]

    @test deaaddnormalizedvrs.slackY ≈ [0.00;
                                6.25;
                                0.00;
                                0.00;
                                13.75;
                                0.00;
                                0.00;
                                0.00;
                                0.00;
                                1.00;
                                0.00]

    # RAM CRS
    deaaddramcrs = deaadd(X, Y, :RAM, rts = :CRS)
    @test deaaddramcrs.weights == :RAM
    @test efficiency(deaaddramcrs) ≈ [0.0000000000;
                                0.14094094094;
                                0.18063063063;
                                0.0000000000;
                                0.27937937938;
                                0.17657657658;
                                0.0000000000;
                                0.17164048866;
                                0.11671671672;
                                0.41766766767;
                                0.06666666667]

    @test deaaddramcrs.slackX ≈ [ 0.000000000  0.000000000;
                                3.037037037  6.814814815;
                                0.000000000  10.837837838;
                                0.000000000  0.000000000;
                                10.592592593 11.037037037;
                                14.666666667 2.666666667;
                                0.000000000  0.000000000;
                                0.000000000  10.298429319;
                                8.296296296  2.518518519;
                                17.925925926 15.370370370;
                                0.000000000  4.000000000]

    @test deaaddramcrs.slackY ≈ [0;
                                   0;
                                   0;
                                   0;
                                   0 ;
                                   0;
                                   0;
                                   0;
                                   0;
                                   0;
                                   0]

    # RAM VRS
    deaaddramvrs = deaadd(X, Y, :RAM, rts = :VRS)
    @test deaaddramvrs.weights == :RAM
    @test efficiency(deaaddramvrs) ≈ [0.0000000000;
                              0.10297482838;
                              0.0000000000;
                              0.0000000000;
                              0.25553012967;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.41764590678 ;
                              0.06666666667]

    @test deaaddramvrs.slackX ≈ [0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                17 15;
                                0 4]
    @test deaaddramvrs.slackY ≈ [0;
                                7.105263158;
                                0;
                                0;
                                17.631578947;
                                0;
                                0;
                                0;
                                0;
                                1.000000000;
                                0]

    # BAM CRS
    deaaddbamcrs = deaadd(X, Y, :BAM, rts = :CRS)
    @test deaaddbamcrs.weights == :BAM
    @test efficiency(deaaddbamcrs) ≈ [0.0000000000;
                                0.4578892372;
                                0.3051750381;
                                0.0000000000;
                                0.6732181854;
                                0.3239568683;
                                0.0000000000;
                                0.5255235602;
                                0.1913580247;
                                0.6902867780;
                                0.1212121212]

    @test deaaddbamcrs.slackX ≈ [0.000000000  0.000000000;
                                4.110344828  6.000000000;
                                0.000000000  0.000000000;
                                0.000000000  0.000000000;
                                13.000000000  8.000000000;
                                17.493670886  0.000000000;
                                0.000000000  0.000000000;
                                0.000000000  9.225130890;
                                8.296296296  2.518518519;
                                13.296296296 13.518518519;
                                0.000000000  4.000000000]
    @test deaaddbamcrs.slackY ≈ [0;
                                0;
                                5.4931506849;
                                0;
                                0.4520547945 ;
                                0;
                                0;
                                1;
                                0;
                                5;
                                0]

    # BAM VRS
    deaaddbamvrs = deaadd(X, Y, :BAM, rts = :VRS)
    @test deaaddbamvrs.weights == :BAM
    @test efficiency(deaaddbamvrs) ≈ [0.0000000000;
                              0.1998936736;
                              0.0000000000;
                              0.0000000000;
                              0.4329710145;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.0000000000;
                              0.5713608345;
                              0.1212121212]

    @test deaaddbamvrs.slackX ≈ [0 0;
                                6.596491228 0;
                                0 0;
                                0 0;
                                13.000000000 1;
                                0 0;
                                0 0;
                                0 0;
                                0 0;
                                5 11;
                                0 4]
    @test deaaddbamvrs.slackY ≈ [0;
                                0;
                                0;
                                0;
                                4;
                                0;
                                0;
                                0;
                                0;
                                5;
                                0]

    # Test model with Custom weights
    deaddcustomcrs = deaadd(X, Y, :Custom, wX = 1 ./ X, wY = 1 ./ Y, rts = :CRS)
    @test deaddcustomcrs.weights == :Custom
    @test efficiency(deaddcustomcrs) ≈ efficiency(deaaddmipcrs)

    deaddcustomvrs = deaadd(X, Y, :Custom, wX = 1 ./ X, wY = 1 ./ Y, rts = :VRS)
    @test deaddcustomvrs.weights == :Custom
    @test efficiency(deaddcustomvrs) ≈ efficiency(deaaddmipvrs)

    ## Test if one-by-one DEA using evaluation and reference sets match initial results
    deaadddefaultcrs_ref_eff = zeros(size(X, 1))
    deaadddefaultvrs_ref_eff = zeros(size(X, 1))

    deaaddcustomcrs_ref_eff = zeros(size(X, 1))
    deaaddcustomvrs_ref_eff = zeros(size(X, 1))

    deaaddonescrs_ref_eff = zeros(size(X, 1))
    deaaddonesvrs_ref_eff = zeros(size(X, 1))

    deaaddmipcrs_ref_eff = zeros(size(X, 1))
    deaaddmipvrs_ref_eff = zeros(size(X, 1))

    Xref = X[:,:]
    Yref = Y[:,:]

    for i = 1:size(X, 1)
        Xeval = X[i:i,:]
        Xeval = Xeval[:,:]
        Yeval = Y[i:i,:]
        Yeval = Yeval[:,:]

        deaadddefaultcrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaadddefaultvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaaddcustomcrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Custom, rts = :CRS, wX = 1 ./ Xeval, wY = 1 ./ Yeval, Xref = Xref, Yref = Yref))[1]
        deaaddcustomvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Custom, rts = :VRS, wX = 1 ./ Xeval, wY = 1 ./ Yeval, Xref = Xref, Yref = Yref))[1]

        deaaddonescrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Ones, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaaddonesvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :Ones, rts = :VRS, Xref = Xref, Yref = Yref))[1]

        deaaddmipcrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :MIP, rts = :CRS, Xref = Xref, Yref = Yref))[1]
        deaaddmipvrs_ref_eff[i] = efficiency(deaadd(Xeval, Yeval, :MIP, rts = :VRS, Xref = Xref, Yref = Yref))[1]

    end

    @test deaadddefaultcrs_ref_eff ≈ efficiency(deaadddefaultcrs)
    @test deaadddefaultvrs_ref_eff ≈ efficiency(deaadddefaultvrs)

    @test deaaddcustomcrs_ref_eff ≈ efficiency(deaaddmipcrs)
    @test deaaddcustomvrs_ref_eff ≈ efficiency(deaaddmipvrs)

    @test deaaddonescrs_ref_eff ≈ efficiency(deaaddcrs1)
    @test deaaddonesvrs_ref_eff ≈ efficiency(deaaddvrs1)

    @test deaaddmipcrs_ref_eff ≈ efficiency(deaaddmipcrs)
    @test deaaddmipvrs_ref_eff ≈ efficiency(deaaddmipvrs)

    # Print
    show(IOBuffer(), deaaddcrs1)

    # Test errors
    @test_throws ErrorException deaadd([1; 2 ; 3], [4 ; 5], :Ones) #  Different number of observations
    @test_throws ErrorException deaadd([1; 2], [4 ; 5], Xref = [1; 2; 3; 4], :Ones) # Different number of observations in reference sets
    @test_throws ErrorException deaadd([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2], :Ones) # Different number of inputs
    @test_throws ErrorException deaadd([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5], :Ones) # Different number of inputs
    @test_throws ErrorException deaadd([1; 2; 3], [4; 5; 6], :Ones, rts = :Error) # Invalid returns to scale
    @test_throws ErrorException deaadd([1; 2; 3], [4; 5; 6], :Error) # Invalid model

    @test_throws ErrorException deaadd([1; 2 ; 3], [4 ; 5]) #  Different number of observations
    @test_throws ErrorException deaadd([1; 2], [4 ; 5], Xref = [1; 2; 3; 4]) # Different number of observations in reference sets
    @test_throws ErrorException deaadd([1 1; 2 2], [4 4; 5 5], Xref = [1 1 1; 2 2 2]) # Different number of inputs
    @test_throws ErrorException deaadd([1 1; 2 2], [4 4; 5 5], Yref = [4 4 4; 5 5 5]) # Different number of inputs
    @test_throws ErrorException deaadd([1; 2; 3], [4; 5; 6], rts = :Error) # Invalid returns to scale
    @test_throws ErrorException deaadd([1; 2; 3], [4; 5; 6], :Custom, wX = [1; 2; 3; 4]) # Different size of weights
    @test_throws ErrorException deaadd([1; 2; 3], [4; 5; 6], :Custom, wY = [4; 5; 6; 7]) # Different size of weights

end


text:
#==============================================================================#
#
#    Declarations of the main `types` used to define Ito diffusions
#    !! IMPORTANT !!
#    to define a diffusion one needs to define a struct and provide a drift
#    and a volatility coefficient functions:
#    `b`: denotes a drift
#    `σ`: denotes a volatility coefficient
#    `a`: (=:σσ') denotes a diffusivity coefficient
#    for linear diffusions
#    `B` and `β` jointly define the drift via `b(t,x):=Bₜx+βₜ`.
#
#==============================================================================#

"""
    DiffusionProcess{T,DP,DW,SS,EI}

Types inheriting from `DiffusionProcess` define Ito diffusions. `T` denotes the
datatype of each coordinate, `DP` the dimension of the stochastic process,
`DW` the dimension of the Wiener process, `SS` lists the state space
restrictions
"""
abstract type DiffusionProcess{T,DP,DW,SS} end

"""
    LinearDiffusion{T,DP,DW,SS} <: DiffusionProcess{T,DP,DW,SS}

Types inheriting from `LinearDiffusion` define a linear Ito-type diffusion, i.e.
solutions to stochastic differential equations of the form:
dXₜ = (BₜXₜ + βₜ)dt + σₜdWₜ, t∈[0,T], X₀=x₀.
"""
abstract type LinearDiffusion{T,DP,DW,SS} <: DiffusionProcess{T,DP,DW,SS} end

"""
    DiffusionDomain

Types inheriting from `DiffusionStateSpace` define the types of restrictions put
on the state-space of the stochastic process.
"""
abstract type DiffusionStateSpace end

const IndexedTime = Tuple{T,S} where {T<:Number,S<:Integer}

#------------------------------------------------------------------------------#
#       Names for functions defining drifts and volatility coefficients
#------------------------------------------------------------------------------#

"""
    b

Compute the drift of a diffusion out-of-place (should use StaticArrays).
"""
function b end

"""
    B

Compute matrix `B` of a linear diffusion out-of-place (should use StaticArrays).
"""
function B end

"""
    β

Compute vector `β` of a linear diffusion out-of-place (should use StaticArrays).
"""
function β end

"""
    σ

Compute the volatility coefficient of a diffusion out-of-place (should use
StaticArrays).
"""
function σ end

"""
    a

Compute the diffusion function σσ' of a diffusion out-of-place (should use
StaticArrays).
"""
function a end

"""
    b!

Compute the drift of a diffusion in-place (uses buffers to store temporary
results).
"""
function b! end

"""
    B!

Compute matrix `B` of a linear diffusion in-place (uses buffers to store
temporary results).
"""
function B! end

"""
    β!

Compute vector `β!` of a linear diffusion in-place (uses buffers to store
temporary results).
"""
function β! end

"""
    σ!

Compute the volatility coefficient of a diffusion in-place (uses buffers to
store temporary results).
"""
function σ! end

"""
    a!

Compute the diffusion function σσ' of a diffusion in-place (uses buffers to
store temporary results).
"""
function a! end

"""
    constdiff

Returns true if the diffusion coefficient does not depend on the state of the
process or time.
"""
function constdiff end

ext: jl


text:
#
# cg-tests.jl -
#
# Tests for linear conjugate gradients methods.
#

#isdefined(:LazyAlgebra) || include("../src/LazyAlgebra.jl")

module TestingLazyAlgebraConjGrad

using LazyAlgebra
using Test

const DEBUG = true

@testset "ConjGrad" begin
    types = (Float32, Float64)
    rows = (7,8)
    cols = (2,3,4)
    @testset "least-square fit ($T)" for T in types
        H = GeneralMatrix(randn(T, rows..., cols...))
        x = randn(T, cols)
        y = H*x + 0.01*randn(T, rows)

        # Least-squares solution using Julia linear algebra.
        yj = reshape(y, length(y))
        Hj = reshape(coefficients(H), length(y), length(x))
        Aj = Hj'*Hj
        bj = Hj'*yj
        x0 = reshape(Aj\bj, cols)

        # Tolerances.
        rtol = 1e-4
        atol = rtol*vnorm2(x0)

        # LHS matrix and RHS vector of the normal equations.
        A = H'*H
        b = H'*y
        n = length(b)
        x1 = conjgrad(A, b; maxiter=10n, restart=n,
                      quiet=false, verb=true,
                      gtol=(0,0), ftol=0, xtol=rtol/10)
        DEBUG && println("x1: ", vnorm2(x1 - x0)/vnorm2(x0))
        @test vnorm2(x1 - x0) ≤ atol

        # Exercise random starting point.
        x2 = conjgrad(A, b, randn(T, cols); maxiter=2n, restart=n,
                      quiet=true, verb=false,
                      gtol=(0,0), ftol=0, xtol=rtol/10)
        DEBUG && println("x2: ", vnorm2(x2 - x0)/vnorm2(x0))
        @test vnorm2(x2 - x0) ≤ atol

        # Exercise other convergence tests.
        x3 = conjgrad(A, b; maxiter=n, restart=n,
                      quiet=false, verb=true,
                      gtol=(0,0), ftol=0, xtol=0)
        DEBUG && println("x3: ", vnorm2(x3 - x0)/vnorm2(x0))
        @test vnorm2(x3 - x0) ≤ atol
        x4 = conjgrad(A, b; maxiter=10n, restart=n,
                      quiet=false, verb=true,
                      gtol=(0,0), ftol=1e-9, xtol=0)
        DEBUG && println("x4: ", vnorm2(x4 - x0)/vnorm2(x0))
        @test vnorm2(x4 - x0) ≤ atol
        x5 = conjgrad(A, b; maxiter=10n, restart=n,
                      quiet=false, verb=true,
                      gtol=(0,1e-6), ftol=0, xtol=0)
        DEBUG && println("x5: ", vnorm2(x5 - x0)/vnorm2(x0))
        @test vnorm2(x5 - x0) ≤ atol

        # Directly use an array.
        x6 = conjgrad(reshape(Aj, cols..., cols...), b;
                      maxiter=2n, restart=n,
                      quiet=true, verb=false,
                      gtol=(0,0), ftol=0, xtol=rtol/10)
        DEBUG && println("x6: ", vnorm2(x6 - x0)/vnorm2(x0))
        @test vnorm2(x6 - x0) ≤ atol

        # Force no iterations (should yield exactly the initial solution).
        x7 = conjgrad(A, b, x0;
                      maxiter=0, restart=n,
                      quiet=true, verb=false,
                      gtol=(0,0), ftol=0, xtol=0)
        DEBUG && println("x7: ", vnorm2(x7 - x0)/vnorm2(x0))
        @test vnorm2(x7 - x0) ≤ 0

        # Force restarts.
        x8 = conjgrad(A, b;
                      maxiter=4n, restart=(n>>1),
                      quiet=true, verb=false,
                      gtol=(0,0), ftol=0, xtol=0)
        DEBUG && println("x8: ", vnorm2(x8 - x0)/vnorm2(x0))
        @test vnorm2(x8 - x0) ≤ atol

        # Non-positive definite.
        @test vnorm2(conjgrad(-A, b; strict=false)) ≤ 0
        @test_throws LazyAlgebra.NonPositiveDefinite conjgrad(-A, b; verb=true,
                                                              quiet=false)

    end
end

nothing

end # module


text:
module DynamicAxisWarping

import ProgressMeter
using Statistics
using StatsBase
using Distances
using FillArrays
using ProgressMeter
using DataStructures
using SlidingDistancesBase
import SlidingDistancesBase: floattype, lastlength, setup_normalizer
using LoopVectorization
using Plots, Plots.PlotMeasures # Plots required both for @layout and for the margins
using Requires
using UnPack

export dtw,
       dtw_cost,
       soft_dtw_cost,
       DTWDistance,
       DTW,
       SoftDTW,
       FastDTW,
       distpath,
       dba,
       dbaclust,
       dtw_cost_matrix,
       soft_dtw_cost_matrix,
       DBAResult,
       fastdtw,
       radiuslimits,
       align_signals,
       dtwnn,
       DTWWorkspace,
       DTWSearchResult,
       GDTWWorkspace,
       sparse_distmat,
       gdtw,
       prepare_gdtw,
       iterative_gdtw!,
       gdtw_warpings,
       LinearInterpolation

export ZNormalizer,
       DiagonalZNormalizer,
       normalize

export stomp

include("utils.jl")

include("distance_interface.jl")
include("filters.jl")
include("dtw.jl")
include("gdtw.jl")
include("dtwnn.jl")
include("dba.jl")
include("dbaclust.jl")
include("windowed_matrix.jl")
include("fastdtw.jl")
include("datasets/datasets.jl")
include("plots.jl")

function __init__()
    @require MatrixProfile = "24e37439-14ec-4097-bda3-6a65822e2305" include("matrix_profile.jl")
end

end # module


text:
@testset "KPLS Pediction Tests (in sample)" begin

    @testset "Test KPLS Single Non Linear Target" begin

        Random.seed!(1)

        z(x)     = 4.26 * (exp.(-x) - 4 * exp.(-2.0*x) + 3 * exp.(-3.0*x))
        x_values = Array(range(0.0,step=3.5,length=100))
        z_pure   = z(x_values)
        noise    = Random.randn(100)
        z_noisy  = z_pure + noise
        X        = MLJBase.table(collect(x_values)[:,:])
        Y        = z_noisy #[:,:] #z_pure

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test MLJBase.mae(yhat, Y[train]) |> mean < 1e-2

    end

    @testset "Test KPLS Single Target (Linear Target)" begin

        X        = MLJBase.table([1 2; 2 4; 4.0 6])
        Y        = [-2; -4; -6.0] #[:,:]

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test MLJBase.mae(yhat, Y[train]) |> mean < 1e-2

        X        = MLJBase.table([1 2; 2 4; 4.0 6])
        Y        = [2; 4; 6.0]#[:,:]

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test MLJBase.mae(yhat, Y[train]) |> mean < 1e-2

    end

    @testset "Test KPLS Multiple Target (Linear Target)" begin

        X        = MLJBase.table([1; 2; 3.0][:,:])
        Y        = MLJBase.table([1 1; 2 2; 3 3.0]) #[:,:]

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test abs.(MLJBase.matrix(yhat) .- MLJBase.matrix(Y)[train,:]) |> mean < 1e-6

        X        = MLJBase.table([1; 2; 3.0][:,:])
        Y        = MLJBase.table([1 -1; 2 -2; 3 -3.0])#[:,:]

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test abs.( MLJBase.matrix(yhat) .- MLJBase.matrix(Y)[train,:]) |> mean < 1e-6

        @testset "Linear Prediction Tests " begin

        X        = MLJBase.table([1 2; 2 4; 4 6.0])
        Y        = MLJBase.table([4 2;6 4;8 6.0])#[:,:]

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test abs.( MLJBase.matrix(yhat) .- MLJBase.matrix(Y)[train,:]) |> mean < 1e-6

        X           = MLJBase.table([1 -2; 2 -4; 4 -6.0])
        Y           = MLJBase.table([-4 -2;-6 -4;-8 -6.0])#[:,:]

        pls_pipe       = MLJBase.@pipeline prediction_type=:deterministic Stand.Standardizer() PartialLeastSquaresRegressor.KPLSRegressor(n_factors=1,kernel="rbf",width=0.01)  target = Stand.Standardizer()

		pls_machine    = MLJBase.machine(pls_pipe, X, Y)

        train = range(1,stop=length(X))
        MLJBase.fit!(pls_machine, rows=train,force=true)
        yhat = MLJBase.predict(pls_machine, rows=train);
        @test abs.(MLJBase.matrix(yhat) .- MLJBase.matrix(Y)[train,:]) |> mean < 1e-6

        end

    end

end;


text:
@testset "util" begin
    x = [1, 2, 3]
    X = [1 0 0; 0 1 0; 0 0 1]

    @test Dolang._unpack_var(x, 1) == 1
    @test Dolang._unpack_var(x, 2) == 2
    @test Dolang._unpack_var(x, 3) == 3

    @test Dolang._unpack_var(X, 1) == [1, 0 ,0]
    @test Dolang._unpack_var(X, 2) == [0, 1 ,0]
    @test Dolang._unpack_var(X, 3) == [0, 0 ,1]

    Dolang._assign_var(x, 3, 1)
    Dolang._assign_var(x, 1, 2)
    Dolang._assign_var(x, 2, 3)
    @test x == [3, 1, 2]

    Dolang._assign_var(X, [0, 1 ,0], 3)
    Dolang._assign_var(X, [0, 0 ,1], 1)
    Dolang._assign_var(X, [1, 0 ,0], 2)
    @test X == [0 1 0; 0 0 1; 1 0 0]

    @test Dolang._output_size(4, x) == (4,)
    @test Dolang._output_size(4, x, x, x, x) == (4,)

    @test Dolang._output_size(4, X) == (3, 4)
    @test Dolang._output_size(4, X, x) == (3, 4)
    @test Dolang._output_size(4, x, X) == (3, 4)
    @test Dolang._output_size(4, X, x, X) == (3, 4)
    @test Dolang._output_size(4, x, X, x) == (3, 4)

    @test_throws DimensionMismatch Dolang._output_size(4, X, rand(4, 4))

    foo = Dolang._allocate_out(Int, 4, x)
    @test isa(foo, Vector{Int})
    @test size(foo) == (4,)

    foo = Dolang._allocate_out(Int, 4, x, x, x, x)
    @test isa(foo, Vector{Int})
    @test size(foo) == (4,)

    foo = Dolang._allocate_out(Int, 4, X)
    @test isa(foo, Matrix{Int})
    @test size(foo) == (3,4)

    foo = Dolang._allocate_out(Int, 4, X, x)
    @test isa(foo, Matrix{Int})
    @test size(foo) == (3,4)

    foo = Dolang._allocate_out(Int, 4, x, X)
    @test isa(foo, Matrix{Int})
    @test size(foo) == (3,4)

    foo = Dolang._allocate_out(Int, 4, x, X, x)
    @test isa(foo, Matrix{Int})
    @test size(foo) == (3,4)

    foo = Dolang._allocate_out(Int, 4, X, x, X)
    @test isa(foo, Matrix{Int})
    @test size(foo) == (3,4)

    @testset "_to_expr" begin
        @test Dolang._to_expr("foo") == Expr(:block, :foo)
        @test Dolang._to_expr(100) == Expr(:block, 100)
        @test Dolang._to_expr(:bar) == Expr(:block, :bar)
        @test Dolang._to_expr(:(x+y)) == :(x+y)
    end

    @testset "inf_to_Inf" begin
        x = rand()
        @test Dolang.inf_to_Inf(x) == x
        @test Dolang.inf_to_Inf(Inf) == Inf
        @test Dolang.inf_to_Inf(:inf) == Inf
        @test Dolang.inf_to_Inf(:(-inf)) == :(-$(Inf))
        @test Dolang.inf_to_Inf(:(x-inf)) == :(x-$(Inf))
    end

    @testset "solve_triangular_system" begin
        # very simple case
        d1 = Dict(:x => 1.0, :y => :(x+1))
        sol1 = OrderedDict{Symbol,Number}(:x => 1.0, :y => 2.0)
        @test Dolang.solve_triangular_system(d1) == sol1

        # fully specified numerical dict
        d2 = Dict(:x => 1.0, :y => 2.0)
        sol2 = OrderedDict{Symbol,Number}(:x => 1.0, :y => 2.0)
        @test Dolang.solve_triangular_system(d2) == sol2

        # unknown variable w
        d3 = Dict(:x => 1.0, :y => :(x+1), :z=> :(100*w))
        @test_throws ErrorException Dolang.solve_triangular_system(d3)

        # non-triangular system
        d4 = Dict(:x => :(y+1), :y => :(x-1))
        @test_throws ErrorException Dolang.solve_triangular_system(d4)
    end

end  # @testset "util"

@testset "triangular systems" begin

    @testset "solve_dependencies" begin
        system_ints = Dict{Int,Set{Int}}(
            1=>Set([2,3,4]),
            2=>Set([3,4]),
            3=>Set([4]),
            4=>Set([])
        )
        sol = Dolang.solve_dependencies(system_ints)
        @test sol==[4,3,2,1]

        system_syms = Dict{Any,Set{Any}}(
            :d=>Set([(:a,0),:b,:c]),
            :c=>Set([(:a,0),:b]),
            :b=>Set([(:a,0)]),
            (:a,0)=>Set([])
        )
        sol = Dolang.solve_dependencies(system_syms)
        @test sol==[(:a,0),:b,:c,:d]

        system_error = Dict{Any,Set{Any}}(
            :d=>Set([(:a,0),:b,:c]),
            :c=>Set([(:a,0),:b]),
            :b=>Set([(:a,0)]),
            (:a,0)=>Set([:d])
        )
        @test_throws Dolang.TriangularSystemException Dolang.solve_dependencies(system_error)
    end

    @testset "Reorder triangular block" begin
        eqs = Dict(
            (:k,0)=>:((1-delta)*k(-1)+i(-1)),
            (:y,0)=>:(A(0)*k(-1)^theta),
            (:A,0)=>:(rho*A(-1)+epsilon),
        )
        eqs_reordered = OrderedDict(
            (:A, 0) => :(rho * A(-1) + epsilon),
            (:y, 0) => :(A(0) * k(-1) ^ theta),
            (:k, 0) => :((1 - delta) * k(-1) + i(-1))
        )
        @test eqs_reordered == Dolang.reorder_triangular_block(eqs)
    end

    @testset "Solve definitions" begin

        defs = Dict(
            (:rho, 0) => :(c[t] / c[t-1]),
            (:V, 0)   => :(c[t] ^ (1 - gamma) / (1 - gamma)),
            (:c, 0)   => :(y[t] ^ theta - i[t]),
        )
        solved_defs = OrderedDict(
            (:V, 0)   => :(c[t] ^ (1 - gamma) / (1 - gamma)),
            (:c, 0)   => :(y[t] ^ theta - i[t]),
            (:c, -1)  => :(y[t-1] ^ theta - i[t-1]),
            (:rho, 0) => :(c[t] / c[t-1])
        )
        @test solved_defs == Dolang.solve_definitions(defs)

        solved_defs_reduced = OrderedDict(
            (:c, 1)   => :(y[t+1] ^ theta - i[t+1]),
            (:c, 0)  => :(y[t] ^ theta - i[t]),
            (:rho, 1) => :(c[t+1] / c[t])
        )

        @test solved_defs_reduced == Dolang.solve_definitions(defs, [(:rho,1)])

    end
end

max_issues_repo_name: EconForge/Dolang


text:
function sumstats(v, f = (x) -> @sprintf("%.2f", x))
  try
    v = convert(Array{Float64,1}, v)
  catch
  end
  "median = $(f(median(v))), mean = $(f(mean(v))) +/- $(f(std(v))), range = [$(f(minimum(v))), $(f(maximum(v)))]"
end

function format_time(t)
  if t < 5e-1
    @sprintf("%.2f ms", t*1e3)
  elseif t < 60.0
    @sprintf("%.2f s", t)
  elseif t < 30*60
    @sprintf("%.2f min", t/60.0)
  else
    @sprintf("%.2f hours", t/3600.0)
  end
end

function unique_filename(prefix = "result", suffix = ".txt")
  join([prefix, strftime("_%Y%m%d_%H%M%S_$(rand(1:int(1e6)))", time()), suffix])
end

function csvfile(header; fileprefix = "experiment",
  filepath = unique_filename(fileprefix, ".csv"), sep = ",")
  fh = open(filepath, "w")
  println(fh, join(header, sep))
  fh
end

# Run repeated experimental runs of a search function for each problem
# in a problem set.
function repeated_runs(searchf, problem_set, num_runs = 10; experiment = "exp")
  problems = collect(problem_set)
  num_problems = length(problems)
  fevals = zeros(Int64, num_runs, num_problems)
  fbests = zeros(num_runs, num_problems)
  times = zeros(num_runs, num_problems)
  reason_counts = [Dict{ASCIIString, Int64}() for i in 1:num_problems]

  file_prefix = strftime("$(experiment)_%Y%m%d_%H%M%S", time())

  run_csvfile = join([file_prefix, "_runs.csv"])

  summary_csvfh = 1 # Dummy value just so it exists...
  include_run_csv_header = true # Only the first time...
  add_summary_csv_header = true # Only the first time...

  # We add one run per problem (with runid 0) to ensure everything has been
  # compiled. It is not saved in the data though.
  for(i in 0:num_runs)
    # Random order of running each problem
    problem_indices = shuffle(collect(1:num_problems))
    for(pi in problem_indices)
      num, prob = problems[pi]
      dims = BlackBoxOptim.numdims(prob)
      println("Run $(i) of problem $(name(prob))")
      start_time = time()
      tic()
      params, paramheader, x, fb, nf, reason, archive = searchf(prob)
      t = toq()
      if i == 0
        break
      end
      times[i, pi] = t
      fbests[i, pi] = fb
      fevals[i, pi] = nf
      reason_counts[pi][reason] = get(reason_counts[pi], reason, 0) + 1

      # Save fitness history to a csv file
      save_fitness_history_to_csv_file(archive, run_csvfile;
        header_prefix = join(["Problem,Dimensions,RunId", paramheader], ","),
        line_prefix = join(["\"$(name(prob))\",$(dims),$(i)", params], ","),
        include_header = include_run_csv_header)
      include_run_csv_header = false # Only the first round...
      println("Saved fitness history to file: $(run_csvfile)")

      # Create the summary_csv file on the first loop through since we now
      # know the number of parameters
      if add_summary_csv_header
        summary_csvfh = csvfile(["Experiment", "Date", "Time", "RunId",
          "Problem", "Dimension", paramheader,
          "ElapsedTime", "FuncEvals", "TerminationReason",
          "Fitness"]; filepath = join([file_prefix, "_summary.csv"]))
        add_summary_csv_header = false
      end

      # Print to summary csv file
      println(summary_csvfh, join([experiment, strftime("%Y-%m-%d", start_time),
        strftime("%T", start_time), i, "\"$(name(prob))\"", dims, params,
        times[i,pi], fevals[i,pi], "\"$(reason)\"", fbests[i,pi]], ","))
      flush(summary_csvfh)
    end
  end

  for(i in 1:num_problems)
    p = problems[i][2]
    println("\nProblem: ", name(p))
    println("Dimensions: ", numdims(p))
    println("Fitness: ", sumstats(fbests[:,i], (x) -> @sprintf("%.2e", x)))
    println("Time: ", sumstats(times[:,i], format_time))
    println("Num. evals: ", sumstats(fevals[:,i]))
    println("Evals / dim: ", sumstats(fevals[:,i] / BlackBoxOptim.numdims(p), round))
    print("Reasons: "); show(reason_counts[i])
    println("")
  end

  close(summary_csvfh)

  println("Finished with experiment $(experiment).")

  return fbests, fevals, times
end

function compare_params(params, searchf, num_runs = 10)
  num_configs = length(params)
  times = zeros(num_runs, num_configs)
  fbests = zeros(num_runs, num_configs)
  fevals = zeros(num_runs, num_configs)

  for(r in 1:num_runs)
    for(i in 1:num_configs)
      tic()
      xb, fb, fe, reason = searchf(params[i]...)
      times[r,i] = toq()
      fbests[r, i] = fb
      fevals[r, i] = float(fe)
    end
  end

  println("\n\nResults per parameter config")
  println("----------------------------")
  sp = sortperm(mean(fbests, 1)[:])
  for(i in 1:num_configs)
    print(i, ". "); show(params[i]); println("")
    println("Fitness: ", sumstats(fbests[:,i], (x) -> @sprintf("%.6e", x)))
    println("Time: ", sumstats(times[:,i], format_time))
    println("Num. evals: ", sumstats(fevals[:,i]))
    println("")
  end

  return times, fbests, fevals
end

function run_test(f, n; mu = false)
  of = xtransform(n, f)
  tic()
  x, f, fevals, reason = cmsa_es(n, of; covarMatrixSampler = EigenCovarSampler,
    mu = mu,
    trace = true, known_fmin = 0.0)
  t = toq()
  print("x = ")
  show(x)
  println("\ntermination reason: $(reason)")
  println("fitness = $(f)")
  println("time = $(t)")
  println("fevals = $(fevals)")
  println("fevals/dim = $(fevals/n)")
end


text:
__precompile__(true)
module VirtualArrays

export VirtualArray, getindex, setindex!, size, length, eachindex
export virtual_cat, virtual_vcat, virtual_hcat

import Base.size, Base.getindex, Base.length, Base.setindex!, Base.eachindex

include("virtual_cat.jl")

type VirtualArray{T, N} <: AbstractArray{T, N}
    expanded_dim::Int # This is the dimension we expand along
    parents::Array{AbstractArray{T},1}

    function VirtualArray{T, N}(expanded_dim::Int, parents::AbstractArray{T,N}...)
        check_parents_dimensions(expanded_dim, parents...)
        return new(expanded_dim, AbstractArray[parent for parent in parents])
    end
    function VirtualArray{T, N}(parents::AbstractArray{T,N}...)
        check_parents_dimensions(1, parents...)
        return new(1, AbstractArray[parent for parent in parents])
    end
end

eachindex(v::VirtualArray) = 1:length(v)

function size(v::VirtualArray)
    result = [0]
    if !isempty(v.parents)
        result = collect(size(v.parents[1]))
        if length(result) >= v.expanded_dim
            expanded_size = 0
            for parent in v.parents
                expanded_size += size(parent)[v.expanded_dim]
            end
            result[v.expanded_dim] = expanded_size
        else
            expanded_size = length(v.parents)
            push!(result, ones(Int, v.expanded_dim - length(result))  ...)
            result[v.expanded_dim] = expanded_size
        end
    end
    return tuple(result...)
end

function length(v::VirtualArray)
    total = 0
    for parent in v.parents
        total += length(parent)
    end
    return total
end

getindex(v::VirtualArray) = getindex(v::VirtualArray, 1)

function getindex(v::VirtualArray, i::Int...)
    checkbounds(v, i...)
    i = expand_index(v, i...)

    for parent in v.parents
        if i[v.expanded_dim] <= get_dimension_size(parent, v.expanded_dim)
            return parent[i...]
        end
        i[v.expanded_dim] -= get_dimension_size(parent, v.expanded_dim)
    end
end

function getindex{T, N}(v::VirtualArray{T, N}, i::UnitRange...)
    checkbounds(v, i...)
    i = collect(i)
    i_dim = length(i)

    if N == i_dim || v.expanded_dim < i_dim
        getindex_range(v, i...)
    else
        num_runs = 1
        for j in v.expanded_dim+1:N
            num_runs *= size(v)[j]
        end

        getindex_range(v, i..., num_runs=num_runs)
    end
end

function setindex!(v::VirtualArray, n, i::Int...)
    checkbounds(v, i...)
    i = expand_index(v, i...)

    for parent in v.parents
        if i[v.expanded_dim] <= get_dimension_size(parent, v.expanded_dim)
            return parent[i...] = n
        end
        i[v.expanded_dim] -= get_dimension_size(parent, v.expanded_dim)
    end
end

function expand_index(v::VirtualArray, i::Int...)
    result = collect(i)
    len_needed = length(size(v))
    len_have = length(i)

    for at in len_have:len_needed - 1
        last_value = result[end]
        result[end] = fix_zero_index(last_value, size(v)[at])
        push!(result, divide(last_value, size(v)[at]))
    end

    return result
end

function fix_zero_index(value::Int, s::Int)
    value = value % s
    if value == 0
        return s
    else
        return value
    end
end

function divide(value::Int, s::Int)
    if value % s == 0
        return div(value, s)
    else
        return div(value, s) + 1
    end
end

function get_dimension_size(v::AbstractArray, i::Int)
    if length(size(v)) >= i
        return size(v)[i]
    else
        return 1
    end
end

function check_parents_dimensions(ignore_dim::Int, parents::AbstractArray...)
    if !isempty(parents)
        default_parent = parents[1]
        default_size = size(parents[1])
        for parent in parents
            curr_size = size(parent)
            for i in max(length(default_size), length(curr_size))
                if i != ignore_dim && get_dimension_size(parent, i) != get_dimension_size(default_parent, i)
                    throw(DimensionMismatch("mismatch in dimension $i"))
                end
            end
        end
    end
end

function getindex_range(v::AbstractArray, i::UnitRange...; num_runs::Int=1)
    result = []
    i = collect(i)

    # For indexing when we have to go through the parents multiple times.
    # For example, a 3d Virtual Array expanded along the 2nd dimensions and 1d indexing
    for j in 1:num_runs
        for parent in v.parents

            # This tells us which index in i we need to change to get the
            # index the user is asking for
            shifting_index = min(length(i), v.expanded_dim)

            # We alter the shifting_index by this much to abstract the array of arrays
            # We need to calculate it like this in case our indexing is < the expansion dimension
            # We need to calcualte for each parent because parents can have different
            # expanded dimensions length
            shift_by = 1
            for k in shifting_index:v.expanded_dim
                shift_by *= get_dimension_size(parent, k)
            end

            # If we have found an index that should be included in the result
            if i[shifting_index].start <= shift_by

                # We need to index into the parent at the right spot
                start_index = max(1,i[shifting_index].start)
                end_index = min(shift_by,i[shifting_index].stop)
                shifting_range = start_index:end_index

                # If we are going through the parents multiple times, we need to adjust accordingly
                shifting_range += shift_by * (j - 1)

                index = (i[1:shifting_index - 1]..., shifting_range, i[shifting_index + 1:end]...)

                push!(result, parent[index...])
            end

            # If we are done adding indexes to our results
            if i[shifting_index].stop <= shift_by
                return cat(shifting_index, result...)
            end

            # Moving onto the next one
            i[shifting_index] -= shift_by
        end
    end
end

end # module


text:
(function (dotwice f)
  (f)
  (f))

(function (printspam)
  (println "spam"))

(dotwice printspam)


text:
# Implementation of new type Vec
struct Vec
    x::Float32
    y::Float32
    z::Float32
end

# Implementation of new type Point
struct Point
    x::Float32
    y::Float32
    z::Float32
end

# Implementation of new type Normal
struct Normal
    x::Float32
    y::Float32
    z::Float32
end

struct Vec2D
    u::Float32
    v::Float32
end

# Implementation of new type Transformation
"""
    Transformation(m, invm)
    Transformation()

It creates a **Transformation**.

## Arguments
- *m* -> matrix of type Array{Array{Float32}};
- *invm* -> inverse matrix.

If not specified, both of the matrixes are the Identity matrix (4x4).
"""
mutable struct Transformation
    m::Matrix{Float32}
    invm::Matrix{Float32}

    Transformation(m=ID4x4, invm=ID4x4) = new(m, invm)
end

# Supporting methods for Transformation
function is_consistent(T::Transformation)
    prod = T.m * T.invm
    return isapprox(prod, ID4x4)
end

ID4x4 = Matrix{Float32}(I, 4, 4)

## BASIC OPERATIONS between Vecs, Normals, Points and Scalars ###############################################################

Base.:isapprox(V1::Vec, V2::Vec) = Base.isapprox(V1.x,V2.x, atol = 5f-6) && Base.isapprox(V1.y,V2.y, atol = 5f-6) && Base.isapprox(V1.z,V2.z, atol = 5f-6)

Base.:isapprox(n1::Normal, n2::Normal) = Base.isapprox(n1.x,n2.x, atol = 5f-6) && Base.isapprox(n1.y,n2.y, atol = 5f-6) && Base.isapprox(n1.z,n2.z, atol = 5f-6)

Base.:+(V1::Vec , V2::Vec )  = Vec((V1.x+V2.x), (V1.y+V2.y), (V1.z+V2.z))

Base.:-(V1::Vec , V2::Vec )   = Vec((V1.x-V2.x), (V1.y-V2.y), (V1.z-V2.z))

Base.:*(V1::Vec , scalar::Real)   = Vec(V1.x*scalar, V1.y*scalar, V1.z*scalar)

Base.:*(scalar::Real, V1::Normal )   = Normal(V1.x*Float32(scalar), V1.y*Float32(scalar), V1.z*Float32(scalar))

Base.:*(V1::Normal , scalar::Real)   = Normal(V1.x*Float32(scalar), V1.y*Float32(scalar), V1.z*Float32(scalar))

Base.:*(scalar::Real, V1::Vec )   = Vec(V1.x*scalar, V1.y*scalar, V1.z*scalar)
Base.:*(V1::Vec , V2::Vec )   = (V1.x*V2.x)+(V1.y*V2.y)+(V1.z*V2.z)

Base.:≈(V1::Vec2D, V2::Vec2D) = Base.isapprox(V1.u, V2.u, atol = 5f-6)  && isapprox(V1.v, V2.v, atol = 5f-6)

# Cross product Vecs
cross(V1::Vec, V2::Vec) = Vec( (V1.y * V2.z - V1.z * V2.y), (V1.z * V2.x - V1.x * V2.z), (V1.x * V2.y - V1.y * V2.x) )

# Norm
squared_norm(V1::Vec) = V1*V1

norm(V1::Vec) = sqrt(squared_norm(V1))

normalize(V1::Vec) = Vec( (V1.x/norm(V1)), (V1.y/norm(V1)), (V1.z/norm(V1)) )

# Aprroximation for two Point
Base.:isapprox(P1::Point, P2::Point) = Base.isapprox(P1.x,P2.x, atol = 5f-6) && Base.isapprox(P1.y,P2.y, atol = 5f-6) && Base.isapprox(P1.z,P2.z, atol = 5f-6)

# Sum Point + Vector returning a Point
Base.:+(P::Point, V::Vec) = Point(P.x + V.x, P.y + V.y, P.z + V.z)

# Sum Point + Point returning a Point
Base.:+(P1::Point, P2::Point) = Point(P1.x + P2.x, P1.y + P2.y, P1.z + P2.z)

# Difference Point - Point returning a Vector
Base.:-(P1::Point, P2::Point) = Vec(P1.x - P2.x, P1.y - P2.y, P1.z - P2.z)

# Difference Point - Vector returning a Point
Base.:-(P::Point, V::Vec) = Point(P.x - V.x, P.y - V.y, P.z - V.z)

# Product Point * scalar
Base.:*(P::Point, a) = Point(P.x*a, P.y*a, P.z*a)

# Convert a Point into a Vec
function toVec(point::Union{Point,Normal,Vec})
    return Vec(point.x,point.y,point.z)
end

function toPoint(vec::Vec)
    return Point(vec.x,vec.y,vec.z)
end
## TRANSFORMATON METHODS #######################################################################################################à

function Base.isapprox(M1::Transformation, M2::Transformation)
    return M1.m ≈ M2.m && M1.invm ≈ M2.invm
end

Base.:*(M1::Transformation, M2::Transformation) = Transformation(M1.m * M2.m, M2.invm * M1.invm)

Base.:*(M::Transformation, P::Point) = Point( P.x * M.m[1,1] + P.y * M.m[1,2] + P.z * M.m[1,3] + M.m[1,4],
                                              P.x * M.m[2,1] + P.y * M.m[2,2] + P.z * M.m[2,3] + M.m[2,4],
                                              P.x * M.m[3,1] + P.y * M.m[3,2] + P.z * M.m[3,3] + M.m[3,4] )

Base.:*(M::Transformation, V::Vec) = Vec( V.x * M.m[1,1] + V.y * M.m[1,2] + V.z * M.m[1,3],
                                          V.x * M.m[2,1] + V.y * M.m[2,2] + V.z * M.m[2,3],
                                          V.x * M.m[3,1] + V.y * M.m[3,2] + V.z * M.m[3,3] )

Base.:*(M::Transformation, N::Normal) = Normal( N.x * M.invm[1,1] + N.y * M.invm[2,1] + N.z * M.invm[3,1],
                                                N.x * M.invm[1,2] + N.y * M.invm[2,2] + N.z * M.invm[3,2],
                                                N.x * M.invm[1,3] + N.y * M.invm[2,3] + N.z * M.invm[3,3] )

function inverse(M::Transformation)
    return Transformation(M.invm, M.m)
end

# Defining translation, scaling and rotation
function translation(vec)
    m = [1.0f0 0.0f0 0.0f0 vec.x;
         0.0f0 1.0f0 0.0f0 vec.y;
         0.0f0 0.0f0 1.0f0 vec.z;
         0.0f0 0.0f0 0.0f0 1.0f0]
    invm = [1.0f0 0.0f0 0.0f0 -vec.x;
            0.0f0 1.0f0 0.0f0 -vec.y;
            0.0f0 0.0f0 1.0f0 -vec.z;
            0.0f0 0.0f0 0.0f0 1.0f0]

    return Transformation(m, invm)
end

function scaling(vec)
    m = [vec.x 0.0f0 0.0f0 0.0f0;
         0.0f0 vec.y 0.0f0 0.0f0;
         0.0f0 0.0f0 vec.z 0.0f0;
         0.0f0 0.0f0 0.0f0 1.0f0]
    invm = [1 / vec.x 0.0f0 0.0f0 0.0f0;
            0.0f0 1 / vec.y 0.0f0 0.0f0;
            0.0f0 0.0f0 1 / vec.z 0.0f0;
            0.0f0 0.0f0 0.0f0 1.0f0]

    return Transformation(m, invm)
end

# Rotations
"""
    rotation_x(angle_rad)

It defines a rotation around the x axis of an angle α (in RADIANTS!!).
"""
function rotation_x(angle_rad)
    sinang, cosang = Float32(sin(angle_rad)), Float32(cos(angle_rad))
    m = [1.0f0 0.0f0 0.0f0 0.0f0;
         0.0f0 cosang -sinang 0.0f0;
         0.0f0 sinang cosang 0.0f0;
         0.0f0 0.0f0 0.0f0 1.0f0]
    invm = [1.0f0 0.0f0 0.0f0 0.0f0;
            0.0f0 cosang sinang 0.0f0;
            0.0f0 -sinang cosang 0.0f0;
            0.0f0 0.0f0 0.0f0 1.0f0]

    return Transformation(m, invm)
end

"""
    rotation_y(angle_rad)

It defines a rotation around the y axis of an angle α (in RADIANTS!!).
"""
function rotation_y(angle_rad)
    sinang, cosang = Float32(sin(angle_rad)), Float32(cos(angle_rad))
    m = [cosang 0.0f0 sinang 0.0f0;
        0.0f0 1.0f0 0.0f0 0.0f0;
        -sinang 0.0f0 cosang 0.0f0;
        0.0f0 0.0f0 0.0f0 1.0f0]
    invm = [cosang 0.0f0 -sinang 0.0f0;
            0.0f0 1.0f0 0.0f0 0.0f0;
            sinang 0.0f0 cosang 0.0f0;
            0.0f0 0.0f0 0.0f0 1.0f0]

    return Transformation(m, invm)
end

"""
    rotation_z(angle_rad)

It defines a rotation around the z axis of an angle α (in RADIANTS!!).
"""
function rotation_z(angle_rad)
    sinang, cosang = Float32(sin(angle_rad)), Float32(cos(angle_rad))
    m = [cosang -sinang 0.0f0 0.0f0;
         sinang cosang 0.0f0 0.0f0;
         0.0f0 0.0f0 1.0f0 0.0f0;
         0.0f0 0.0f0 0.0f0 1.0f0]
    invm = [cosang sinang 0.0f0 0.0f0;
            -sinang cosang 0.0f0 0.0f0;
            0.0f0 0.0f0 1.0f0 0.0f0;
            0.0f0 0.0f0 0.0f0 1.0f0]

    return Transformation(m, invm)
end

###############################################################################
function create_onb(normal_passed::Union{Normal,Vec})
    normal = normalize(toVec(normal_passed))
    sign = copysign(1., normal.z)
    a = (-1. / (sign + normal.z))
    b = (normal.x * normal.y * a)

    e1 = Vec((1. + sign * normal.x * normal.x * a), (sign * b), (-sign * normal.x))
    e2 = Vec((b), (sign + normal.y * normal.y * a), (-normal.y))

    return e1, e2, Vec(normal.x, normal.y, normal.z)
end


text:
if !("." in LOAD_PATH) # for easier local testing
    push!(LOAD_PATH, ".")
end
using Test
using SingleColumnModels
using OrderedCollections
using SingleColumnModels.StateVecs
const scm = SingleColumnModels
using CLIMAParameters
using CLIMAParameters.Planet
import CLIMAParameters.Planet

ENV["GKSwstype"] = "nul"

const scm_dir = dirname(dirname(pathof(SingleColumnModels)));
output_dir = joinpath(scm_dir, "output", "BOMEX")
mkpath(output_dir)
const pycles_comparison = joinpath(output_dir, "pycles_comparison")
const plot_dir = joinpath(output_dir, "plots")

struct EarthParameterSet <: AbstractEarthParameterSet end
CLIMAParameters.Planet.MSLP(::EarthParameterSet) = 100000.0
const param_set = EarthParameterSet()

best_mse = OrderedDict()
best_mse[:q_tot_gm] = 6.7598618147035705e-01
best_mse[:a_up] = 1.1151497561239513e+02
best_mse[:w_up] = 6.7164800402645284e+01
best_mse[:q_tot_up] = 1.4426382216507470e+01
best_mse[:θ_liq_up] = 1.0040573552118677e+02
best_mse[:tke_en] = 9.4610544693497943e+01

# Define `compute_mse` and retrieve data file:
include(joinpath(scm_dir, "integration_tests", "utils", "compute_mse.jl"))
include(joinpath(scm_dir, "integration_tests", "utils", "plot_profiles.jl"))
data_file = Dataset(joinpath(PyCLES_output_dataset_path, "Bomex.nc"), "r")

@testset "BOMEX EDMF Solution Quality Assurance (QA) tests" begin

    grid, q, aux, params = scm.EDMF.run(param_set, scm.EDMF.BOMEX(), output_dir)
    computed_mse = compute_mse(
        grid,
        q,
        aux,
        params,
        data_file,
        "Bomex",
        best_mse,
        21600.0,
        pycles_comparison,
    )

    # Export plots (must export before tests):
    skip_fields = ["a_gm", "u_en", "u_ud_1", "v_en", "v_gm", "v_ud_1", "w_gm"]
    prog_fn = joinpath(output_dir, "prog_vs_time.nc")
    plot_profiles(prog_fn, plot_dir; skip_fields = skip_fields)
    # aux_fn = joinpath(output_dir, "aux_vs_time.nc")
    # plot_profiles(aux_fn, plot_dir)

    test_mse(computed_mse, best_mse, :q_tot_gm)
    test_mse(computed_mse, best_mse, :a_up)
    test_mse(computed_mse, best_mse, :w_up)
    test_mse(computed_mse, best_mse, :q_tot_up)
    test_mse(computed_mse, best_mse, :θ_liq_up)
    test_mse(computed_mse, best_mse, :tke_en)
end


text:
abstract type MemoryLayout end
struct UnknownLayout <: MemoryLayout end
abstract type AbstractStridedLayout <: MemoryLayout end
abstract type AbstractIncreasingStrides <: AbstractStridedLayout end
abstract type AbstractColumnMajor <: AbstractIncreasingStrides end
struct DenseColumnMajor <: AbstractColumnMajor end
struct ColumnMajor <: AbstractColumnMajor end
struct IncreasingStrides <: AbstractIncreasingStrides end
abstract type AbstractDecreasingStrides <: AbstractStridedLayout end
abstract type AbstractRowMajor <: AbstractDecreasingStrides end
struct DenseRowMajor <: AbstractRowMajor end
struct RowMajor <: AbstractRowMajor end
struct DecreasingStrides <: AbstractIncreasingStrides end
struct StridedLayout <: AbstractStridedLayout end
struct ScalarLayout <: MemoryLayout end

"""
    UnknownLayout()

is returned by `MemoryLayout(A)` if it is unknown how the entries of an array `A`
are stored in memory.
"""
UnknownLayout

"""
    AbstractStridedLayout

is an abstract type whose subtypes are returned by `MemoryLayout(A)`
if an array `A` has storage laid out at regular offsets in memory,
and which can therefore be passed to external C and Fortran functions expecting
this memory layout.

Julia's internal linear algebra machinery will automatically (and invisibly)
dispatch to BLAS and LAPACK routines if the memory layout is BLAS compatible and
the element type is a `Float32`, `Float64`, `ComplexF32`, or `ComplexF64`.
In this case, one must implement the strided array interface, which requires
overrides of `strides(A::MyMatrix)` and `unknown_convert(::Type{Ptr{T}}, A::MyMatrix)`.
"""
AbstractStridedLayout

"""
    DenseColumnMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
equivalent to an `Array`, so that `stride(A,1) == 1` and
`stride(A,i) ≡ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`. In particular,
if `A` is a matrix then `strides(A) == `(1, size(A,1))`.

Arrays with `DenseColumnMajor` memory layout must conform to the `DenseArray` interface.
"""
DenseColumnMajor

"""
    ColumnMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a column major array, so that `stride(A,1) == 1` and
`stride(A,i) ≥ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`.

Arrays with `ColumnMajor` memory layout must conform to the `DenseArray` interface.
"""
ColumnMajor

"""
    IncreasingStrides()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a strided array with  increasing strides, so that `stride(A,1) ≥ 1` and
`stride(A,i) ≥ size(A,i-1) * stride(A,i-1)` for `2 ≤ i ≤ ndims(A)`.
"""
IncreasingStrides

"""
    DenseRowMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a row major array with dense entries, so that `stride(A,ndims(A)) == 1` and
`stride(A,i) ≡ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`. In particular,
if `A` is a matrix then `strides(A) == `(size(A,2), 1)`.
"""
DenseRowMajor

"""
    RowMajor()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a row major array, so that `stride(A,ndims(A)) == 1` and
stride(A,i) ≥ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`.

If `A` is a matrix  with `RowMajor` memory layout, then
`transpose(A)` should return a matrix whose layout is `ColumnMajor`.
"""
RowMajor

"""
    DecreasingStrides()

is returned by `MemoryLayout(A)` if an array `A` has storage in memory
as a strided array with decreasing strides, so that `stride(A,ndims(A)) ≥ 1` and
stride(A,i) ≥ size(A,i+1) * stride(A,i+1)` for `1 ≤ i ≤ ndims(A)-1`.
"""
DecreasingStrides

"""
    StridedLayout()

is returned by `MemoryLayout(A)` if an array `A` has storage laid out at regular
offsets in memory. `Array`s with `StridedLayout` must conform to the `DenseArray` interface.
"""
StridedLayout

"""
    ScalarLayout()

is returned by `MemoryLayout(A)` if A is a scalar, which does not live in memory
"""
ScalarLayout

"""
    MemoryLayout(A)

specifies the layout in memory for an array `A`. When
you define a new `AbstractArray` type, you can choose to override
`MemoryLayout` to indicate how an array is stored in memory.
For example, if your matrix is column major with `stride(A,2) == size(A,1)`,
then override as follows:

    MemoryLayout(::MyMatrix) = DenseColumnMajor()

The default is `UnknownLayout()` to indicate that the layout
in memory is unknown.

Julia's internal linear algebra machinery will automatically (and invisibly)
dispatch to BLAS and LAPACK routines if the memory layout is compatible.
"""
@inline MemoryLayout(_) = UnknownLayout()

@inline MemoryLayout(::Number) = ScalarLayout()
@inline MemoryLayout(::DenseArray) = DenseColumnMajor()

@inline MemoryLayout(A::ReinterpretArray) = reinterpretedmemorylayout(MemoryLayout(parent(A)))
@inline reinterpretedmemorylayout(::MemoryLayout) = UnknownLayout()
@inline reinterpretedmemorylayout(::DenseColumnMajor) = DenseColumnMajor()

@inline MemoryLayout(A::ReshapedArray) = reshapedmemorylayout(MemoryLayout(parent(A)))
@inline reshapedmemorylayout(::MemoryLayout) = UnknownLayout()
@inline reshapedmemorylayout(::DenseColumnMajor) = DenseColumnMajor()

@inline MemoryLayout(A::SubArray) = subarraylayout(MemoryLayout(parent(A)), parentindices(A))
subarraylayout(_1, _2) = UnknownLayout()
subarraylayout(_1, _2, _3)= UnknownLayout()
subarraylayout(::DenseColumnMajor, ::Tuple{I}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    DenseColumnMajor()  # A[:] is DenseColumnMajor if A is DenseColumnMajor
subarraylayout(ml::AbstractColumnMajor, inds) = _column_subarraylayout1(ml, inds)
subarraylayout(::AbstractRowMajor, ::Tuple{I}) where I =
    UnknownLayout()  # A[:] does not have any structure if A is AbstractRowMajor
subarraylayout(ml::AbstractRowMajor, inds) = _row_subarraylayout1(ml, reverse(inds))
subarraylayout(ml::AbstractStridedLayout, inds) = _strided_subarraylayout(ml, inds)

_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:Union{Int,AbstractCartesianIndex} =
    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:Slice =
    DenseColumnMajor() # view(A,:,1,2) is a DenseColumnMajor vector
_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:AbstractUnitRange{Int} =
    DenseColumnMajor() # view(A,1:3,1,2) is a DenseColumnMajor vector
_column_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:Union{Int,AbstractCartesianIndex} =
    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
_column_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:AbstractUnitRange{Int} =
    DenseColumnMajor() # view(A,1:3,1,2) is a DenseColumnMajor vector
_column_subarraylayout1(::DenseColumnMajor, inds::Tuple{I,Vararg{Any}}) where I<:Slice =
    _column_subarraylayout(DenseColumnMajor(), DenseColumnMajor(), tail(inds))
_column_subarraylayout1(par::DenseColumnMajor, inds::Tuple{I,Vararg{Any}}) where I<:AbstractUnitRange{Int} =
    _column_subarraylayout(par, ColumnMajor(), tail(inds))
_column_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:AbstractUnitRange{Int} =
    _column_subarraylayout(par, ColumnMajor(), tail(inds))
_column_subarraylayout1(par::DenseColumnMajor, inds::Tuple{I,Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _column_subarraylayout(par, StridedLayout(), tail(inds))
_column_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _column_subarraylayout(par, StridedLayout(), tail(inds))
_column_subarraylayout1(par, inds) = UnknownLayout()
_column_subarraylayout(par, ret, ::Tuple{}) = ret
_column_subarraylayout(par, ret, ::Tuple{I}) where I = UnknownLayout()
_column_subarraylayout(::DenseColumnMajor, ::DenseColumnMajor, inds::Tuple{I,Vararg{Int}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    DenseColumnMajor() # A[:,1:3,1,2] is DenseColumnMajor if A is DenseColumnMajor
_column_subarraylayout(par::DenseColumnMajor, ::DenseColumnMajor, inds::Tuple{I, Vararg{Int}}) where I<:Slice =
    DenseColumnMajor()
_column_subarraylayout(par::DenseColumnMajor, ::DenseColumnMajor, inds::Tuple{I, Vararg{Any}}) where I<:Slice =
    _column_subarraylayout(par, DenseColumnMajor(), tail(inds))
_column_subarraylayout(par, ::AbstractColumnMajor, inds::Tuple{I, Vararg{Any}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    _column_subarraylayout(par, ColumnMajor(), tail(inds))
_column_subarraylayout(par, ::AbstractStridedLayout, inds::Tuple{I, Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _column_subarraylayout(par, StridedLayout(), tail(inds))

_row_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:Union{Int,AbstractCartesianIndex} =
    DenseColumnMajor() # view(A,1,1,2) is a scalar, which we include in DenseColumnMajor
_row_subarraylayout1(::DenseRowMajor, inds::Tuple{I,Vararg{Int}}) where I<:Slice =
    DenseColumnMajor() # view(A,1,2,:) is a DenseColumnMajor vector
_row_subarraylayout1(par, inds::Tuple{I,Vararg{Int}}) where I<:AbstractUnitRange{Int} =
    DenseColumnMajor() # view(A,1,2,1:3) is a DenseColumnMajor vector
_row_subarraylayout1(::DenseRowMajor, inds::Tuple{I,Vararg{Any}}) where I<:Slice =
    _row_subarraylayout(DenseRowMajor(), DenseRowMajor(), tail(inds))
_row_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:AbstractUnitRange{Int} =
    _row_subarraylayout(par, RowMajor(), tail(inds))
_row_subarraylayout1(par, inds::Tuple{I,Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _row_subarraylayout(par, StridedLayout(), tail(inds))
_row_subarraylayout1(par, inds) = UnknownLayout()
_row_subarraylayout(par, ret, ::Tuple{}) = ret
_row_subarraylayout(par, ret, ::Tuple{I}) where I = UnknownLayout()
_row_subarraylayout(::DenseRowMajor, ::DenseRowMajor, inds::Tuple{I,Vararg{Int}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    DenseRowMajor() # A[1,2,1:3,:] is DenseRowMajor if A is DenseRowMajor
_row_subarraylayout(par::DenseRowMajor, ::DenseRowMajor, inds::Tuple{I, Vararg{Int}}) where I<:Slice =
    DenseRowMajor()
_row_subarraylayout(par::DenseRowMajor, ::DenseRowMajor, inds::Tuple{I, Vararg{Any}}) where I<:Slice =
    _row_subarraylayout(par, DenseRowMajor(), tail(inds))
_row_subarraylayout(par::AbstractRowMajor, ::AbstractRowMajor, inds::Tuple{I, Vararg{Any}}) where I<:Union{AbstractUnitRange{Int},Int,AbstractCartesianIndex} =
    _row_subarraylayout(par, RowMajor(), tail(inds))
_row_subarraylayout(par::AbstractRowMajor, ::AbstractStridedLayout, inds::Tuple{I, Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _row_subarraylayout(par, StridedLayout(), tail(inds))

_strided_subarraylayout(par, inds) = UnknownLayout()
_strided_subarraylayout(par, ::Tuple{}) = StridedLayout()
_strided_subarraylayout(par, inds::Tuple{I, Vararg{Any}}) where I<:Union{RangeIndex,AbstractCartesianIndex} =
    _strided_subarraylayout(par, tail(inds))

# MemoryLayout of transposed and adjoint matrices
struct ConjLayout{ML<:MemoryLayout} <: MemoryLayout
    layout::ML
end

conjlayout(_1, _2) = UnknownLayout()
conjlayout(::Type{<:Complex}, M::ConjLayout) = M.layout
conjlayout(::Type{<:Complex}, M::AbstractStridedLayout) = ConjLayout(M)
conjlayout(::Type{<:Real}, M::MemoryLayout) = M

subarraylayout(M::ConjLayout, t::Tuple) = ConjLayout(subarraylayout(M.layout, t))

MemoryLayout(A::Transpose) = transposelayout(MemoryLayout(parent(A)))
MemoryLayout(A::Adjoint) = adjointlayout(eltype(A), MemoryLayout(parent(A)))
transposelayout(_) = UnknownLayout()
transposelayout(::StridedLayout) = StridedLayout()
transposelayout(::ColumnMajor) = RowMajor()
transposelayout(::RowMajor) = ColumnMajor()
transposelayout(::DenseColumnMajor) = DenseRowMajor()
transposelayout(::DenseRowMajor) = DenseColumnMajor()
transposelayout(M::ConjLayout) = ConjLayout(transposelayout(M.layout))
adjointlayout(::Type{T}, M::MemoryLayout) where T = transposelayout(conjlayout(T, M))

# MemoryLayout of Symmetric/Hermitian
"""
    SymmetricLayout(layout, uplo)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
as a symmetrized version of `layout`, where the entries used are dictated by the
`uplo`, which can be `'U'` or `L'`.

A matrix that has memory layout `SymmetricLayout(layout, uplo)` must overrided
`symmetricdata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] == B[k,j]` for `j ≥ k` if `uplo == 'U'` (`j ≤ k` if `uplo == 'L'`) and
`A[k,j] == B[j,k]` for `j < k` if `uplo == 'U'` (`j > k` if `uplo == 'L'`).
"""
struct SymmetricLayout{ML<:MemoryLayout} <: MemoryLayout
    layout::ML
    uplo::Char
end
SymmetricLayout(layout::ML, uplo) where ML<:MemoryLayout = SymmetricLayout{ML}(layout, uplo)

"""
    HermitianLayout(layout, uplo)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
as a hermitianized version of `layout`, where the entries used are dictated by the
`uplo`, which can be `'U'` or `L'`.

A matrix that has memory layout `HermitianLayout(layout, uplo)` must overrided
`hermitiandata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] == B[k,j]` for `j ≥ k` if `uplo == 'U'` (`j ≤ k` if `uplo == 'L'`) and
`A[k,j] == conj(B[j,k])` for `j < k` if `uplo == 'U'` (`j > k` if `uplo == 'L'`).
"""
struct HermitianLayout{ML<:MemoryLayout} <: MemoryLayout
    layout::ML
    uplo::Char
end
HermitianLayout(layout::ML, uplo) where ML<:MemoryLayout = HermitianLayout{ML}(layout, uplo)

MemoryLayout(A::Hermitian) = hermitianlayout(eltype(A), MemoryLayout(parent(A)), A.uplo)
MemoryLayout(A::Symmetric) = symmetriclayout(MemoryLayout(parent(A)), A.uplo)
hermitianlayout(_1, _2, _3) = UnknownLayout()
hermitianlayout(::Type{<:Complex}, layout::AbstractColumnMajor, uplo) = HermitianLayout(layout,uplo)
hermitianlayout(::Type{<:Real}, layout::AbstractColumnMajor, uplo) = SymmetricLayout(layout,uplo)
hermitianlayout(::Type{<:Complex}, layout::AbstractRowMajor, uplo) = HermitianLayout(layout,uplo)
hermitianlayout(::Type{<:Real}, layout::AbstractRowMajor, uplo) = SymmetricLayout(layout,uplo)
symmetriclayout(_1, _2) = UnknownLayout()
symmetriclayout(layout::AbstractColumnMajor, uplo) = SymmetricLayout(layout,uplo)
symmetriclayout(layout::AbstractRowMajor, uplo) = SymmetricLayout(layout,uplo)
transposelayout(S::SymmetricLayout) = S
adjointlayout(::Type{T}, S::SymmetricLayout) where T<:Real = S
adjointlayout(::Type{T}, S::HermitianLayout) where T = S
subarraylayout(S::SymmetricLayout, ::Tuple{<:Slice,<:Slice}) = S
subarraylayout(S::HermitianLayout, ::Tuple{<:Slice,<:Slice}) = S

symmetricdata(A::Symmetric) = A.data
symmetricdata(A::Hermitian{<:Real}) = A.data
symmetricdata(V::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = symmetricdata(parent(V))
symmetricdata(V::Adjoint{<:Real}) = symmetricdata(parent(V))
symmetricdata(V::Transpose) = symmetricdata(parent(V))
hermitiandata(A::Hermitian) = A.data
hermitiandata(V::SubArray{<:Any, 2, <:Any, <:Tuple{<:Slice,<:Slice}}) = hermitiandata(parent(V))
hermitiandata(V::Adjoint) = hermitiandata(parent(V))
hermitiandata(V::Transpose{<:Real}) = hermitiandata(parent(V))

# MemoryLayout of triangular matrices
struct TriangularLayout{UPLO,UNIT,ML} <: MemoryLayout
    layout::ML
end

TriangularLayout{UPLO,UNIT}(lay) where {UPLO,UNIT} = TriangularLayout{UPLO,UNIT,typeof(lay)}(lay)

"""
    LowerTriangularLayout(layout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `LowerTriangular(B)` where `B` satisfies `MemoryLayout(B) == layout`.

A matrix that has memory layout `LowerTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ zero(eltype(A))` for `j > k` and
`A[k,j] ≡ B[k,j]` for `j ≤ k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `UpperTriangularLayout`.
"""
const LowerTriangularLayout{ML} = TriangularLayout{'L','N',ML}

"""
    UnitLowerTriangularLayout(ML::MemoryLayout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `UnitLowerTriangular(B)` where `B` satisfies `MemoryLayout(B) == layout`.

A matrix that has memory layout `UnitLowerTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ zero(eltype(A))` for `j > k`,
`A[k,j] ≡ one(eltype(A))` for `j == k`,
`A[k,j] ≡ B[k,j]` for `j < k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `UnitUpperTriangularLayout`.
"""
const UnitLowerTriangularLayout{ML} = TriangularLayout{'L','U',ML}

"""
    UpperTriangularLayout(ML::MemoryLayout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `UpperTriangularLayout(B)` where `B` satisfies `MemoryLayout(B) == ML`.

A matrix that has memory layout `UpperTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ B[k,j]` for `j ≥ k` and
`A[k,j] ≡ zero(eltype(A))` for `j < k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `LowerTriangularLayout`.
"""
const UpperTriangularLayout{ML} = TriangularLayout{'U','N',ML}

"""
    UnitUpperTriangularLayout(ML::MemoryLayout)

is returned by `MemoryLayout(A)` if a matrix `A` has storage in memory
equivalent to a `UpperTriangularLayout(B)` where `B` satisfies `MemoryLayout(B) == ML`.

A matrix that has memory layout `UnitUpperTriangularLayout(layout)` must overrided
`triangulardata(A)` to return a matrix `B` such that `MemoryLayout(B) == layout` and
`A[k,j] ≡ B[k,j]` for `j > k`,
`A[k,j] ≡ one(eltype(A))` for `j == k`,
`A[k,j] ≡ zero(eltype(A))` for `j < k`.

Moreover, `transpose(A)` and `adjoint(A)` must return a matrix that has memory
layout `UnitLowerTriangularLayout`.
"""
UnitUpperTriangularLayout{ML} = TriangularLayout{'U','U',ML}

MemoryLayout(A::UpperTriangular) = triangularlayout(UpperTriangularLayout, MemoryLayout(parent(A)))
MemoryLayout(A::UnitUpperTriangular) = triangularlayout(UnitUpperTriangularLayout, MemoryLayout(parent(A)))
MemoryLayout(A::LowerTriangular) = triangularlayout(LowerTriangularLayout, MemoryLayout(parent(A)))
MemoryLayout(A::UnitLowerTriangular) = triangularlayout(UnitLowerTriangularLayout, MemoryLayout(parent(A)))
triangularlayout(_, ::MemoryLayout) = UnknownLayout()
triangularlayout(::Type{Tri}, ML::AbstractColumnMajor) where {Tri} = Tri(ML)
triangularlayout(::Type{Tri}, ML::AbstractRowMajor) where {Tri} = Tri(ML)
triangularlayout(::Type{Tri}, ML::ConjLayout{<:AbstractRowMajor}) where {Tri} = Tri(ML)
subarraylayout(layout::TriangularLayout, ::Tuple{<:Union{Slice,Base.OneTo},<:Union{Slice,Base.OneTo}}) = layout
conjlayout(::Type{<:Complex}, ml::TriangularLayout{UPLO,UNIT}) where {UPLO,UNIT} =
    TriangularLayout{UPLO,UNIT}(ConjLayout(ml.layout))

for (TriLayout, TriLayoutTrans) in ((UpperTriangularLayout,     LowerTriangularLayout),
                                    (UnitUpperTriangularLayout, UnitLowerTriangularLayout),
                                    (LowerTriangularLayout,     UpperTriangularLayout),
                                    (UnitLowerTriangularLayout, UnitUpperTriangularLayout))
    @eval transposelayout(ml::$TriLayout) = $TriLayoutTrans(transposelayout(ml.layout))
end

triangulardata(A::AbstractTriangular) = parent(A)
triangulardata(A::Adjoint) = Adjoint(triangulardata(parent(A)))
triangulardata(A::Transpose) = Transpose(triangulardata(parent(A)))
triangulardata(A::SubArray{<:Any,2,<:Any,<:Tuple{<:Union{Slice,Base.OneTo},<:Union{Slice,Base.OneTo}}}) =
    view(triangulardata(parent(A)), parentindices(A)...)

abstract type AbstractBandedLayout <: MemoryLayout end

struct DiagonalLayout{ML} <: AbstractBandedLayout
    layout::ML
end

struct SymTridiagonalLayout{ML} <: AbstractBandedLayout
    layout::ML
end

MemoryLayout(D::Diagonal) = DiagonalLayout(MemoryLayout(parent(D)))
diagonaldata(D::Diagonal) = parent(D)

MemoryLayout(D::SymTridiagonal) = SymTridiagonalLayout(MemoryLayout(D.dv))
diagonaldata(D::SymTridiagonal) = D.dv
offdiagonaldata(D::SymTridiagonal) = D.ev

transposelayout(ml::DiagonalLayout) = ml
transposelayout(ml::SymTridiagonalLayout) = ml

adjointlayout(_, ml::DiagonalLayout) = ml
adjointlayout(::Type{<:Real}, ml::SymTridiagonal) = ml

###
# Fill
####
abstract type AbstractFillLayout <: MemoryLayout end
struct FillLayout <: AbstractFillLayout end
struct ZerosLayout <: AbstractFillLayout end

MemoryLayout(::AbstractFill) = FillLayout()
MemoryLayout(::Zeros) = ZerosLayout()


text:
using Pkg

Pkg.activate(joinpath(@__DIR__, "..")); Pkg.instantiate()
Pkg.activate(); Pkg.instantiate();

pushfirst!(LOAD_PATH, joinpath(@__DIR__, ".."))

using SkipLists
using BenchmarkTools

# Run all benchmarks for the module that are contained in this directory

include("benchmark_generation.jl")
display(benchmark_generate_vector())
display(benchmark_generate_skiplist(SkipList))


text:
using ClangCompiler
using Test

const CC = ClangCompiler

@testset "FileManager" begin
    fm = CC.FileManager()
    @test fm.ptr != C_NULL

    @testset "FileEntry" begin
        p = joinpath(@__DIR__, "..", "code", "main.cpp") |> normpath
        f = CC.getFileEntry(fm, p)
        @test f.ptr != C_NULL
        @test CC.get_name(f) == p == CC.real_path_name(f)
        @test CC.isValid(f)
        @test !CC.isNamedPipe(f)
    end

    dispose(fm)
end

@testset "SourceManager" begin
    fm = CC.FileManager()
    @test fm.ptr != C_NULL

    sm = CC.SourceManager(fm)
    @test sm.ptr != C_NULL

    code = """
    int x;
    int y;
    """
    buffer = get_buffer(code)

    fid = CC.FileID(sm, buffer)
    @test fid.ptr != C_NULL

    CC.setMainFileID(sm, fid)

    fid2 = CC.getMainFileID(sm)
    @test fid2.ptr != C_NULL

    @test CC.value(fid) == CC.value(fid2)

    dispose(fid)
    dispose(fid2)
    dispose(sm)
    dispose(fm)
end


text:
include("codec.jl")

using MsgPack
using CodecZlib

export decodedata,encodedata,parsemmtf,writemmtf,fetchmmtf

"""
Decodes the MsgPack unpacked dict of the MMTF file.
"""
function decodedata(input_data::Dict)
    # Initialize optional field's default values
    dict = Dict(
        "bFactorList" => Array{Float32,1}[],
        "occupancyList" => Array{Float32,1}[],
        "atomIdList" => Array{Int32,1}[],
        "altLocList" => Array{Char,1}[],
        "insCodeList" => Array{Char,1}[],
        "sequenceIndexList" => Array{Int32,1}[],
        "chainNameList" => Array{String,1}[],
        "spaceGroup" => "",
        "bondAtomList" => Array{Int32,1}[],
        "bondOrderList" => Array{Int8,1}[],
        "mmtfVersion" => "",
        "mmtfProducer" => "",
        "structureId" => "",
        "title" => "",
        "experimentalMethods" => Array{Any,1}[],
        "depositionDate" => "",
        "releaseDate" => "",
        "entityList" => Array{Any,1}[],
        "bioAssemblyList" => Array{Any,1}[],
        "rFree" => "",
        "rWork" => "",
        "resolution" => "",
        "unitCell" => Array{Any,1}[],
        "secStructList" => Array{Int8,1}[],
        "ncsOperatorList" => Array{Any,1}[]
    )
    # Check if all mandatory fields are present in the input data.
    mandateFields = ["groupTypeList","xCoordList","yCoordList","zCoordList","groupIdList","groupList","chainsPerModel","groupsPerChain","chainIdList","numBonds","numChains","numModels","numAtoms","numGroups"]
    for field in mandateFields
        if !haskey(input_data, field)
            throw(ArgumentError("Mandatory field \"$field\" not available in input data!"))
        end
    end
    # decode the necessary fields alone
    toDecodeFields = ["groupTypeList","xCoordList","yCoordList","zCoordList","bFactorList","occupancyList","atomIdList","altLocList","insCodeList","groupIdList","sequenceIndexList","chainNameList","chainIdList","bondAtomList","bondOrderList","secStructList"]
    for (key,value) in input_data
        if key in toDecodeFields
            dict[key] = decodearray(value)
        else
            dict[key] = value
        end
    end
    # return the final dict
    return dict
end

"""
Encodes the MMTF Dict.
"""
function encodedata(MMTFDict::Dict)
    encodedMMTFDict = Dict()
    encodedMMTFDict["groupTypeList"] = encodearray(MMTFDict["groupTypeList"], 4, 0)
    encodedMMTFDict["xCoordList"] = encodearray(MMTFDict["xCoordList"], 10, COORD_DIVIDER)
    encodedMMTFDict["yCoordList"] = encodearray(MMTFDict["yCoordList"], 10, COORD_DIVIDER)
    encodedMMTFDict["zCoordList"] = encodearray(MMTFDict["zCoordList"], 10, COORD_DIVIDER)
    encodedMMTFDict["bFactorList"] = encodearray(MMTFDict["bFactorList"], 10, OCC_B_FACTOR_DIVIDER)
    encodedMMTFDict["occupancyList"] = encodearray(MMTFDict["occupancyList"], 9, OCC_B_FACTOR_DIVIDER)
    encodedMMTFDict["atomIdList"] = encodearray(MMTFDict["atomIdList"], 8, 0)
    encodedMMTFDict["altLocList"] = encodearray(MMTFDict["altLocList"], 6, 0)
    encodedMMTFDict["insCodeList"] = encodearray(MMTFDict["insCodeList"], 6, 0)
    encodedMMTFDict["groupIdList"] = encodearray(MMTFDict["groupIdList"], 8, 0)
    encodedMMTFDict["sequenceIndexList"] = encodearray(MMTFDict["sequenceIndexList"], 8, 0)
    encodedMMTFDict["chainNameList"] = encodearray(MMTFDict["chainNameList"], 5, 4)
    encodedMMTFDict["chainIdList"] = encodearray(MMTFDict["chainIdList"], 5, 4)
    encodedMMTFDict["bondAtomList"] = encodearray(MMTFDict["bondAtomList"], 4, 0)
    encodedMMTFDict["bondOrderList"] = encodearray(MMTFDict["bondOrderList"], 2, 0)
    encodedMMTFDict["secStructList"] = encodearray(MMTFDict["secStructList"], 2, 0)
    encodedMMTFDict["groupList"] = MMTFDict["groupList"]
    encodedMMTFDict["chainsPerModel"] = MMTFDict["chainsPerModel"]
    encodedMMTFDict["groupsPerChain"] = MMTFDict["groupsPerChain"]
    encodedMMTFDict["spaceGroup"] = MMTFDict["spaceGroup"]
    encodedMMTFDict["mmtfVersion"] = MMTFDict["mmtfVersion"]
    encodedMMTFDict["mmtfProducer"] = MMTFDict["mmtfProducer"]
    encodedMMTFDict["structureId"] = MMTFDict["structureId"]
    encodedMMTFDict["entityList"] = MMTFDict["entityList"]
    encodedMMTFDict["bioAssemblyList"] = MMTFDict["bioAssemblyList"]
    encodedMMTFDict["rFree"] = MMTFDict["rFree"]
    encodedMMTFDict["rWork"] = MMTFDict["rWork"]
    encodedMMTFDict["resolution"] = MMTFDict["resolution"]
    encodedMMTFDict["title"] = MMTFDict["title"]
    encodedMMTFDict["experimentalMethods"] = MMTFDict["experimentalMethods"]
    encodedMMTFDict["depositionDate"] = MMTFDict["depositionDate"]
    encodedMMTFDict["releaseDate"] = MMTFDict["releaseDate"]
    encodedMMTFDict["unitCell"] = MMTFDict["unitCell"]
    encodedMMTFDict["numBonds"] = MMTFDict["numBonds"]
    encodedMMTFDict["numChains"] = MMTFDict["numChains"]
    encodedMMTFDict["numModels"] = MMTFDict["numModels"]
    encodedMMTFDict["numAtoms"] = MMTFDict["numAtoms"]
    encodedMMTFDict["numGroups"] = MMTFDict["numGroups"]
    encodedMMTFDict["ncsOperatorList"] = MMTFDict["ncsOperatorList"]
    return encodedMMTFDict
end

"""
Parses the MMTF file into a MMTFDict.
filepath: Path of the input file.
"""
function parsemmtf(filepath::AbstractString; gzip::Bool=false)
    return open(filepath) do input
        parsemmtf(input,gzip=gzip)
    end
end

function parsemmtf(input::IO; gzip::Bool=false)
    if gzip
        decompressed = IOBuffer()
        write(decompressed, GzipDecompressorStream(input))
        decodedata(unpack(take!(decompressed)))
    else
        decodedata(unpack(input))
    end
end

"""
Writes the MMTFDict back to a file.
MMTFDict: Dict containing the decoded MMTF data to be written to file.
filepath: Path to write the file
"""
function writemmtf(MMTFDict::Dict, filepath::AbstractString; gzip::Bool=false)
    open(filepath, "w") do output
        writemmtf(MMTFDict,output,gzip=gzip)
    end
end

function writemmtf(MMTFDict::Dict, output::IO; gzip::Bool=false)
    if gzip
        stream=IOBuffer(pack(encodedata(MMTFDict)))
        write(output,GzipCompressorStream(stream))
    else
        pack(output,encodedata(MMTFDict))
    end
end

"""
Fetches a file from RCSB server and decodes it into a MMTFDict.
pdbid: ID of the PDB file to be fetched and decoded.
"""
function fetchmmtf(pdbid)
    tempfile = tempname()
    download("$(BASE_URL)/$(pdbid).mmtf.gz",tempfile)
    parsemmtf(tempfile,gzip=true)
end


text:
"""
≪(μ,ν)

# Absolute continuity

The following are equivalent:
1. μ ≪ ν
2. μ is absolutely continuous wrt ν
3. There exists a function f such that μ = ∫f dν
"""
function ≪ end

"""
≃(μ,ν)

# Equivalence of Measure

Measures μ and ν on the same space X are equivalent, written `μ ≃ ν`, if `μ ≪ ν` and `ν ≪ μ`.

Note that this is often written `~` in the literature, but this is overloaded in probabilistic programming, so we use alternate notation.
"""
function ≃(μ,ν)
return (μ≪ν && ν≪μ)
end