text:
using Statistics,Serialization,StatsPlots, DataFrames, LaTeXStrings

filename_prefix = "..//data//main//xi06//"
plotname_prefix = "..//figures//main_xi06_comparison_" # store the plots here

show_capacity = false # for infected plot

# set intial values as in coivd par ini
datat = 7 # interval of data collection
T = 582 #582  # time periods

# define period of plotting for economic plots
start_week = 2 # cut of first week
end_week = 84
# set periods for covid plots
inf_start_day = 1
inf_end_day = 400 # cut of after period 400
inf_start_week = 1
inf_end_week = 58 # 58*7 = 406 and 42*7 = 294

datapoint = fld(T,datat)
include("$(filename_prefix)beta5//alpha1//covid_par_ini.jl")

### blue LINES
# parameter values for xi = 0.6
paras = ["beta5//alpha25","beta5//alpha50","beta5//alpha75","beta5//alpha1"]

# scatterplot
gdploss_means = []
totcas_means = []

for (i,p) in enumerate(paras)
    data_prefix = string("$(filename_prefix)", p , "//")
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    push!(gdploss_means,gdploss)
    push!(totcas_means,totcas)
    if p == "beta5//alpha1"
        global point_A_x = gdploss
        global point_A_y = totcas
    end
    if p == "beta5//alpha25"
        global point_B_x = gdploss
        global point_B_y = totcas
    end
end

ll = length(paras)
gdploss_means = reverse(gdploss_means)
totcas_means = reverse(totcas_means)
suppfig3a = quiver(gdploss_means[1:ll-1],totcas_means[1:ll-1],quiver=(gdploss_means[2:ll]-gdploss_means[1:ll-1],totcas_means[2:ll]-totcas_means[1:ll-1]), color = [:red], linewidth = 2, xlabel = "Average GDP loss [%]", ylabel = "Mortality [%]")
scatter!(suppfig3a, gdploss_means, totcas_means, markersize = 5, markercolor = [:red], legend=false)

fig2a = quiver(gdploss_means[1:ll-1],totcas_means[1:ll-1],quiver=(gdploss_means[2:ll]-gdploss_means[1:ll-1],totcas_means[2:ll]-totcas_means[1:ll-1]), color = [:red], linewidth = 2, xlabel = "Average GDP loss [%]", ylabel = "Mortality [%]")
scatter!(fig2a, gdploss_means, totcas_means, markersize = 5, markercolor = [:red], legend=false)

# set boundaries
ylims!(suppfig3a,0.008,0.07)
xlims!(suppfig3a,0,10)

ylims!(fig2a,0.008,0.055)
xlims!(fig2a,0,10)

### black LINES
paras = ["beta10//alpha25","beta10//alpha50","beta10//alpha75","beta10//alpha1"]

# scatterplot
gdploss_means2 = []
totcas_means2 = []

for (i,p) in enumerate(paras)
    data_prefix = string("$(filename_prefix)", p , "//")
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    push!(gdploss_means2,gdploss)
    push!(totcas_means2,totcas)
end

ll = length(paras)
gdploss_means2 = reverse(gdploss_means2)
totcas_means2 = reverse(totcas_means2)
quiver!(suppfig3a, gdploss_means2[1:ll-1],totcas_means2[1:ll-1],quiver=(gdploss_means2[2:ll]-gdploss_means2[1:ll-1],totcas_means2[2:ll]-totcas_means2[1:ll-1]), color = [:black], linewidth = 2)
scatter!(suppfig3a, gdploss_means2, totcas_means2, markersize = 5, markercolor = [:black], legend=false)

### Green LINES
paras = ["beta30//alpha25","beta30//alpha50","beta30//alpha75","beta30//alpha1"]

# scatterplot
gdploss_means3 = []
totcas_means3 = []

for (i,p) in enumerate(paras)
    data_prefix = string("$(filename_prefix)", p , "//")
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    push!(gdploss_means3,gdploss)
    push!(totcas_means3,totcas)
    if p == "beta30//alpha25"
        global point_D_x = gdploss
        global point_D_y = totcas
    end
end

ll = length(paras)
gdploss_means3 = reverse(gdploss_means3)
totcas_means3 = reverse(totcas_means3)
quiver!(suppfig3a, gdploss_means3[1:ll-1],totcas_means3[1:ll-1],quiver=(gdploss_means3[2:ll]-gdploss_means3[1:ll-1],totcas_means3[2:ll]-totcas_means3[1:ll-1]), color = [:green], linewidth = 2)
scatter!(suppfig3a, gdploss_means3, totcas_means3, markersize = 5, markercolor = [:green], legend=false)

### blue LINES
paras = ["beta5//alpha1","beta10//alpha1","beta30//alpha1"]

# scatterplot
gdploss_means4 = []
totcas_means4 = []

for (i,p) in enumerate(paras)
    data_prefix = string("$(filename_prefix)", p , "//")
    # store data scatterplot
    worker_results = deserialize(open("$(data_prefix)batchdata.dat"))
    gdploss = mean((results -> results[:togdploss]).(worker_results))
    totcas = mean(100*(results -> results[:totcas]).(worker_results))
    gdpmean = mean(hcat((results -> results[:gdppercaptraj]).(worker_results)...), dims=2)
    global nn = size(gdpmean)[1]
    push!(gdploss_means4,gdploss)
    push!(totcas_means4,totcas)
    if p == "beta30//alpha1"
        global point_C_x = gdploss
        global point_C_y = totcas
    end
end

ll = length(paras)
quiver!(suppfig3a, gdploss_means4[1:ll-1],totcas_means4[1:ll-1],quiver=(gdploss_means4[2:ll]-gdploss_means4[1:ll-1],totcas_means4[2:ll]-totcas_means4[1:ll-1]), color = [:blue], linewidth = 2)
scatter!(suppfig3a, gdploss_means4, totcas_means4, markersize = 5, markercolor = [:blue], legend=false)

quiver!(fig2a, gdploss_means4[1:ll-1],totcas_means4[1:ll-1],quiver=(gdploss_means4[2:ll]-gdploss_means4[1:ll-1],totcas_means4[2:ll]-totcas_means4[1:ll-1]), color = [:blue], linewidth = 2)
scatter!(fig2a, gdploss_means4, totcas_means4, markersize = 5, markercolor = [:blue], legend=false)

#### add specific points A B C
scatter!(suppfig3a, [point_A_x], [point_A_y], markersize = 6, markercolor = [:black])
annotate!(suppfig3a,[point_A_x+0.4], [point_A_y+0.0005], text.([L"A_1"]) , fontsize=60 )

scatter!(suppfig3a, [point_B_x], [point_B_y], markersize = 6, markercolor = [:black])
annotate!(suppfig3a,[point_B_x], [point_B_y-0.003], text.([L"B_1"]) , fontsize=160 )

scatter!(suppfig3a, [point_C_x], [point_C_y], markersize = 6, markercolor = [:black])
annotate!(suppfig3a,[point_C_x+0.4], [point_C_y], text.([L"C_1"]) , fontsize=160 )

scatter!(suppfig3a, [point_D_x], [point_D_y], markersize = 6, markercolor = [:black])
annotate!(suppfig3a,[point_D_x], [point_D_y-0.003], text.([L"D_1"]) , fontsize=160 )

plot!(suppfig3a, size=(600,400))
savefig(suppfig3a,"$(plotname_prefix)supp_fig_3a.pdf")

### add specific points A B C
scatter!(fig2a, [point_A_x], [point_A_y], markersize = 6, markercolor = [:black])
annotate!(fig2a,[point_A_x+0.4], [point_A_y+0.0005], text.([L"A_1"]) , fontsize=60 )
# scatter!(plarrows, [point_A_x-0.3], [point_A_y], markersize = 0, markercolor = [:white], series_annotations = text.([L"\textrm{\sffamily A_1}"], :top))

scatter!(fig2a, [point_B_x], [point_B_y], markersize = 6, markercolor = [:black])
annotate!(fig2a,[point_B_x+0.4], [point_B_y+0.0005], text.([L"B_1"]) , fontsize=160 )
#scatter!(plarrows, [point_B_x], [point_B_y-0.003], markersize = 0, markercolor = [:white], series_annotations = text.([L"B_1"], :top))

scatter!(fig2a, [point_C_x], [point_C_y], markersize = 6, markercolor = [:black])
annotate!(fig2a,[point_C_x+0.4], [point_C_y], text.([L"C_1"]) , fontsize=160 )

plot!(fig2a, size=(600,400))
savefig(fig2a,"$(plotname_prefix)fig_2a.pdf")


text:
module Regression

using LinearAlgebra
using Statistics
export LinearRegression, Lasso, Ridge, make_design_matrix

include("utils.jl")
include("Base.jl")
include("BFM.jl")
include("LR.jl")
include("RR.jl")

end

ext: jl


text:
# Indexing --------------------
function get{T, S}(tn::TreeNode{T, S})
    return tn
end

function getindex{T, S}(tn::TreeNode{T, S}, key::S)
    tn.children[key]
end

function getindex{T, S}(tn::TreeNode{T, S}, key::Vector{S})
    if isempty(key) error("Empty path/key") end

    keys = copy(key)

    current_key = shift!(keys)

    if haskey(tn.children, current_key)
        if length(key) > 0
            return(getindex(tn.children[current_key], keys))
        else
            return(tn.children[current_key])
        end
    else
        KeyError(current_key)
    end
end

"""
Verify that a node is a root node.
"""
function isroot{T, S}(tn::TreeNode{T, S})
    isnull(tn.parent)
end

function isroot{T, S}(tn::Nullable{TreeNode{T, S}})
    isnull(tn) ? Nullable{bool}() : isroot(get(tn))
end

"""
Verify if a node is a leaf (ie does not have children).
"""
function isleaf{T, S}(tn::TreeNode{T, S})
    return length(tn.children) == 0
end

function isleaf{T, S}(tn::Nullable{TreeNode{T, S}})
    return isnull(tn) ? Nullable{Bool}() : length(get(tn).children) == 0
end

"""
Find the root node of a tree.
"""
function getroot{T, S}(tn::TreeNode{T, S})
    if isroot(tn)
        return(tn)
    else
        getroot(get(tn.parent))
    end
end

"""
Compute depth of a node (its distance from root)
"""
function getdepth{T, S}(tn::TreeNode{T, S})
    if isroot(tn)
        return(0)
    else
        return(1 + getdepth(get(tn.parent)))
    end
end

function getdepth{T, S}(tn::Nullable{TreeNode{T, S}})
    if isnull(tn)
        return 0
    else
        return getdepth(get(tn))
    end
end

"""
Get sibilings of a node.
"""
function getsiblings{T, S}(tn::TreeNode{T, S})
    if (!isnull(tn.parent))
        allchildren = copy(get(tn.parent).children)
        delete!(allchildren, tn.id) # remove itself from the result
    else nothing end
end


text:
# test_bb_operations()
A = BoundingBox(10,10,30,20)
B = BoundingBox(5,20,20,10)
C = A+B
D = BoundingBox(5,10,35,20)
@test C == D

# @test_throws InexactError BoundingBox(0,0,-1,-1)

# test_find_mesh_bb()
  dst = [4.0 2.0;
          8.0 2.0;
          6.0 10.0]
  bounds = find_mesh_bb(dst)
  @test bounds == BoundingBox(4,2,4,8)

  dst = [0.0 0.0;
          0.0 0.0;
          0.0 0.0]
  bounds = find_mesh_bb(dst)
  @test bounds == BoundingBox(0,0,0,0)

# test_tform_bb()
tform = [1 0 0;
        0 1 0;
        0 0 1];
bb = BoundingBox(100,100,100,100)
tbb = tform_bb(bb, tform)
@test bb == tbb

tform = [1 0 0;
        0 1 0;
        100 100 1];
bb = BoundingBox(100,100,100,100)
tbb = tform_bb(bb, tform)
@test BoundingBox(200,200,100,100) == tbb

sz = (100, 100)
tform = [cos(pi/6) -sin(pi/6) 0;
        sin(pi/6) cos(pi/6) 0;
        0 0 1];
bb = BoundingBox(0, -50, 136.60254037844388, 136.60254037844388)
tbb = tform_bb(sz_to_bb(sz), tform)
@test_approx_eq bb_to_pts(bb) bb_to_pts(tbb)

# test_snap_bb()
bb = snap_bb(sz_to_bb((100.5, 200.75)))
tbb = sz_to_bb((101, 201))
@test bb_to_pts(bb) == bb_to_pts(tbb)

bb = snap_bb(sz_to_bb((100.9, 200.1)))
tbb = sz_to_bb((101, 201))
@test bb_to_pts(bb) == bb_to_pts(tbb)

bb = snap_bb(BoundingBox(100.9, 200.1, 100.5, 100.2))
tbb = sz_to_bb((101, 201))
tbb = BoundingBox(100, 200, 102, 101)
@test bb_to_pts(bb) == bb_to_pts(tbb)

a = BoundingBox(0,0,10,10)
b = BoundingBox(0,0,5,5)
@test a+b == BoundingBox(0,0,10,10)
@test a-b == BoundingBox(0,0,5,5)

a = BoundingBox(2,2,10,10)
b = BoundingBox(0,0,5,5)
@test a+b == BoundingBox(0,0,12,12)
@test a-b == BoundingBox(2,2,3,3)

a = BoundingBox(20,20,10,10)
b = BoundingBox(0,0,5,5)
@test a+b == BoundingBox(0,0,30,30)
@test isnan((a-b).i)


text:
@safetestset "Output """).body |> String ==
    """
Good morning
"""

  @test Html.html("""
Good morning
""").body |> String ==
    """
Good morning
"""

  @test Html.html("""
Good morning
""").body |> String ==
    """
Good morning
"""

  @test Html.html("""
Good morning
""").body |> String ==
    """
Good morning
"""

  @test Html.html("""
Good morning
""").body |> String ==
    """
Good morning
"""

  @test Html.html(filepath("views/outputscripttags.jl.html")).body |> String |> fws ==
    """<!DOCTYPE html>
Greetings

    """ |> fws

  @test Html.html(filepath("views/outputscripttags.jl.html"), layout=filepath("views/layoutscripttags.jl.html")).body |> String |> fws ==
    """<!DOCTYPE html>
Layout header

Greetings

Layout footer

    """ |> fws
end


text:
 @testset "LineSegment-QuadraticSegment" begin
     for T ∈ Floats
        # intersect
        X₁ = Point{2,T}(0, 0)
        X₂ = Point{2,T}(2, 0)
        X₃ = Point{2,T}(1, 1)
        X₄ = Point{2,T}(1, 0)
        X₅ = Point{2,T}(1, 2)

        # 1 intersection, straight
        q = QuadraticSegment(X₁, X₂, Point{2,T}(1//2, 0))
        l = LineSegment(Point{2,T}(1,-1), Point{2,T}(1,1))
        (point1, point2) = intersect(l, q)
        @test point1 ≈ Point{2,T}(1, 0)
        @test point2 ≈ Point{2,T}(1e6, 1e6)

        # 1 intersection
        q = QuadraticSegment(X₁, X₂, X₃)
        l = LineSegment(X₄, X₅)
        (point1, point2) = intersect(l, q)
        @test point1 ≈ Point{2,T}(1, 1)
        @test point2 ≈ Point{2,T}(1e6, 1e6)

        # 2 intersections
        X₄ = Point{2,T}(0, 3//4)
        X₅ = Point{2,T}(2, 3//4)
        l = LineSegment(X₄, X₅)
        (point1, point2) = l ∩ q
        @test point1 ≈ Point{2,T}(1//2, 3//4)
        @test point2 ≈ Point{2,T}(3//2, 3//4)

        # 0 intersections
        X₄ = Point{2,T}(0, 3)
        X₅ = Point{2,T}(2, 3)
        l = LineSegment(X₄, X₅)
        (point1, point2) = intersect(l, q)
        @test point1 ≈ Point{2,T}(1e6, 1e6)
        @test point2 ≈ Point{2,T}(1e6, 1e6)
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
       myFunction (a;
                   x::Float64=sin(1)) = 0


text:
# This file is a part of Julia. License is MIT: https://julialang.org/license

## Basic functions ##

"""
    AbstractQuasiArray{T,N}

Supertype for `N`-dimensional arrays (or array-like types) with elements of type `T`.
[`Array`](@ref) and other types are subtypes of this. See the manual section on the
[`AbstractQuasiArray` interface](@ref man-interface-array).
"""
AbstractQuasiArray

convert(::Type{T}, a::T) where {T<:AbstractQuasiArray} = a
convert(::Type{AbstractQuasiArray{T}}, a::AbstractQuasiArray) where {T} = AbstractQuasiArray{T}(a)
convert(::Type{AbstractQuasiArray{T,N}}, a::AbstractQuasiArray{<:Any,N}) where {T,N} = AbstractQuasiArray{T,N}(a)

Array{T}(a::AbstractQuasiArray) where T = T[a[k] for k in eachindex(a)]
Array{T,N}(a::AbstractQuasiArray{<:Any,N}) where {T,N} = Array{T}(a)
Array(a::AbstractQuasiArray{T}) where T = Array{T}(a)
Matrix(a::AbstractQuasiMatrix{T}) where T = Array{T}(a)
Vector(a::AbstractQuasiVector{T}) where T = Array{T}(a)
AbstractArray(a::AbstractQuasiArray) = Array(a)
AbstractMatrix(a::AbstractQuasiMatrix) = Matrix(a)
AbstractVector(a::AbstractQuasiVector) = Vector(a)
AbstractArray{T}(a::AbstractQuasiArray) where T = Array{T}(a)
AbstractMatrix{T}(a::AbstractQuasiMatrix) where T = Matrix{T}(a)
AbstractVector{T}(a::AbstractQuasiVector) where T = Vector{T}(a)

convert(::Type{Array{T}}, a::AbstractQuasiArray) where T = Array{T}(a)
convert(::Type{Array{T,N}}, a::AbstractQuasiArray{<:Any,N}) where {T,N} = Array{T}(a)
convert(::Type{Array}, a::AbstractQuasiArray{T}) where T = Array{T}(a)
convert(::Type{Matrix}, a::AbstractQuasiMatrix{T}) where T = Array{T}(a)
convert(::Type{Vector}, a::AbstractQuasiVector{T}) where T = Array{T}(a)

convert(::Type{AbstractArray{T}}, a::AbstractQuasiArray) where T = Array{T}(a)
convert(::Type{AbstractArray}, a::AbstractQuasiArray{T}) where T = convert(AbstractArray{T}, a)
convert(::Type{AbstractArray{T,N}}, a::AbstractQuasiArray{<:Any,N}) where {T,N} = convert(AbstractArray{T}, a)
convert(::Type{AbstractMatrix}, a::AbstractQuasiMatrix) = convert(AbstractArray, a)
convert(::Type{AbstractVector}, a::AbstractQuasiVector) = convert(AbstractArray, a)

"""
    size(A::AbstractQuasiArray, [dim])

Return a tuple containing the dimensions of `A`. Optionally you can specify a
dimension to just get the length of that dimension.

Note that `size` may not be defined for arrays with non-standard indices, in which case [`axes`](@ref)
may be useful. See the manual chapter on [arrays with custom indices](@ref man-custom-indices).

# Examples
```jldoctest
julia> A = fill(1, (2,3,4));

julia> size(A)
(2, 3, 4)

julia> size(A, 2)
3
```
"""
size(t::AbstractQuasiArray{T,N}, d) where {T,N} = d <= N ? size(t)[d] : 1

"""
    axes(A, d)

Return the valid range of indices for array `A` along dimension `d`.

See also [`size`](@ref), and the manual chapter on [arrays with custom indices](@ref man-custom-indices).

# Examples
```jldoctest
julia> A = fill(1, (5,6,7));

julia> axes(A, 2)
Base.OneTo(6)
```
"""
function axes(A::AbstractQuasiArray{T,N}, d) where {T,N}
    @_inline_meta
    d <= N ? axes(A)[d] : OneTo(1)
end

# Performance optimization: get rid of a branch on `d` in `axes(A, d)`
# for d=1. 1d arrays are heavily used, and the first dimension comes up
# in other applications.
axes1(A::AbstractQuasiArray{<:Any,0}) = OneTo(1)
axes1(A::AbstractQuasiArray) = (@_inline_meta; axes(A)[1])

keys(a::AbstractQuasiArray) = CartesianIndices(axes(a))
keys(a::AbstractQuasiVector) = LinearIndices(a)

eltype(::Type{<:AbstractQuasiArray{E}}) where {E} = @isdefined(E) ? E : Any
elsize(A::AbstractQuasiArray) = elsize(typeof(A))

"""
    ndims(A::AbstractQuasiArray) -> Integer

Return the number of dimensions of `A`.

# Examples
```jldoctest
julia> A = fill(1, (3,4,5));

julia> ndims(A)
3
```
"""
ndims(::AbstractQuasiArray{T,N}) where {T,N} = N
ndims(::Type{<:AbstractQuasiArray{T,N}}) where {T,N} = N

"""
    length(collection) -> Integer

Return the number of elements in the collection.

Use [`lastindex`](@ref) to get the last valid index of an indexable collection.

# Examples
```jldoctest
julia> length(1:5)
5

julia> length([1, 2, 3, 4])
4

julia> length([1 2; 3 4])
4
```
"""
length

"""
    length(A::AbstractQuasiArray)

Return the number of elements in the array, defaults to `prod(size(A))`.

# Examples
```jldoctest
julia> length([1, 2, 3, 4])
4

julia> length([1 2; 3 4])
4
```
"""
length(t::AbstractQuasiArray) = (@_inline_meta; prod(size(t)))

# eachindex iterates over all indices. IndexCartesian definitions are later.
eachindex(A::AbstractQuasiVector) = (@_inline_meta(); axes1(A))

"""
    eachindex(A...)

Create an iterable object for visiting each index of an `AbstractQuasiArray` `A` in an efficient
manner. For array types that have opted into fast linear indexing (like `Array`), this is
simply the range `1:length(A)`. For other array types, return a specialized Cartesian
range to efficiently index into the array with indices specified for every dimension. For
other iterables, including strings and dictionaries, return an iterator object
supporting arbitrary index types (e.g. unevenly spaced or non-integer indices).

If you supply more than one `AbstractQuasiArray` argument, `eachindex` will create an
iterable object that is fast for all arguments (a [`UnitRange`](@ref)
if all inputs have fast linear indexing, a [`CartesianIndices`](@ref)
otherwise).
If the arrays have different sizes and/or dimensionalities, `eachindex` will return an
iterable that spans the largest range along each dimension.

# Examples
```jldoctest
julia> A = [1 2; 3 4];

julia> for i in eachindex(A) # linear indexing
           println(i)
       end
1
2
3
4

julia> for i in eachindex(view(A, 1:2, 1:1)) # Cartesian indexing
           println(i)
       end
CartesianIndex(1, 1)
CartesianIndex(2, 1)
```
"""
eachindex(A::AbstractQuasiArray) = (@_inline_meta(); eachindex(IndexStyle(A), A))

function eachindex(A::AbstractQuasiArray, B::AbstractQuasiArray)
    @_inline_meta
    eachindex(IndexStyle(A,B), A, B)
end
function eachindex(A::AbstractQuasiArray, B::AbstractQuasiArray...)
    @_inline_meta
    eachindex(IndexStyle(A,B...), A, B...)
end
eachindex(::IndexLinear, A::AbstractQuasiArray) = (@_inline_meta; OneTo(length(A)))
eachindex(::IndexLinear, A::AbstractQuasiVector) = (@_inline_meta; axes1(A))
function eachindex(::IndexLinear, A::AbstractQuasiArray, B::AbstractQuasiArray...)
    @_inline_meta
    indsA = eachindex(IndexLinear(), A)
    _all_match_first(X->eachindex(IndexLinear(), X), indsA, B...) ||
        throw_eachindex_mismatch(IndexLinear(), A, B...)
    indsA
end

# keys with an IndexStyle
keys(s::IndexStyle, A::AbstractQuasiArray, B::AbstractQuasiArray...) = eachindex(s, A, B...)

"""
    lastindex(collection) -> index
    lastindex(collection, d) -> index

Return the last index of `collection`. If `d` is given, return the last index of `collection` along dimension `d`.

The syntaxes `A[end]` and `A[end, end]` lower to `A[lastindex(A)]` and
`A[lastindex(A, 1), lastindex(A, 2)]`, respectively.

# Examples
```jldoctest
julia> lastindex([1,2,4])
3

julia> lastindex(rand(3,4,5), 2)
4
```
"""
lastindex(a::AbstractQuasiArray) = (@_inline_meta; last(eachindex(IndexLinear(), a)))
lastindex(a::AbstractQuasiArray, d) = (@_inline_meta; last(axes(a, d)))

"""
    firstindex(collection) -> index
    firstindex(collection, d) -> index

Return the first index of `collection`. If `d` is given, return the first index of `collection` along dimension `d`.

# Examples
```jldoctest
julia> firstindex([1,2,4])
1

julia> firstindex(rand(3,4,5), 2)
1
```
"""
firstindex(a::AbstractQuasiArray) = (@_inline_meta; first(eachindex(IndexLinear(), a)))
firstindex(a::AbstractQuasiArray, d) = (@_inline_meta; first(axes(a, d)))

first(a::AbstractQuasiArray) = a[first(eachindex(a))]
stride(A::AbstractQuasiArray, k::Integer) = strides(A)[k]

function isassigned(a::AbstractQuasiArray, i...)
    try
        a[i...]
        true
    catch e
        if isa(e, BoundsError) || isa(e, UndefRefError)
            return false
        else
            rethrow(e)
        end
    end
end

function checkbounds(::Type{Bool}, A::AbstractQuasiArray, I...)
    @_inline_meta
    checkbounds_indices(Bool, axes(A), I)
end

# As a special extension, allow using logical arrays that match the source array exactly
function checkbounds(::Type{Bool}, A::AbstractQuasiArray{<:Any,N}, I::AbstractQuasiArray{Bool,N}) where N
    @_inline_meta
    axes(A) == axes(I)
end

function checkbounds(A::AbstractQuasiArray, I...)
    @_inline_meta
    checkbounds(Bool, A, to_indices(A,I)...) || throw_boundserror(A, I)
    nothing
end

const QuasiDimOrInd = Union{Number, AbstractQuasiOrVector{<:Number}}

similar(a::AbstractQuasiArray{T}) where {T}                             = similar(a, T)
similar(a::AbstractQuasiArray, ::Type{T}) where {T}                     = similar(a, T, axes(a))
similar(a::AbstractQuasiArray{T}, dims::Tuple) where {T}                = similar(a, T, dims)
similar(a::AbstractQuasiArray{T}, dims::QuasiDimOrInd...) where {T}          = similar(a, T, dims)
similar(a::AbstractQuasiArray, ::Type{T}, dims::QuasiDimOrInd...) where {T}  = similar(a, T, dims)
similar(::Type{<:AbstractQuasiArray{T}}, shape::NTuple{N,AbstractQuasiOrVector}) where {N,T} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, shape))
similar(a::AbstractQuasiArray, ::Type{T}, dims::NTuple{N,AbstractQuasiOrVector}) where {T,N} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, dims))
similar(a::AbstractQuasiArray, ::Type{T}, dims::Vararg{AbstractQuasiOrVector,N}) where {T,N} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, dims))

similar(a::AbstractQuasiArray{T}, m::Int) where {T}              = Vector{T}(undef, m)
similar(a::AbstractQuasiArray, T::Type, dims::Dims{N}) where {N} = Array{T,N}(undef, dims)
similar(a::AbstractQuasiArray{T}, dims::Dims{N}) where {T,N}     = Array{T,N}(undef, dims)

similar(::Type{T}, dims::QuasiDimOrInd...) where {T<:AbstractQuasiArray} = similar(T, dims)
similar(::Type{T}, shape::Tuple{AbstractQuasiVector,Vararg{QuasiDimOrInd}}) where {N,T<:AbstractArray{<:Any,N}} =
    QuasiArray{T,N}(undef, convert.(AbstractVector, shape))

empty(a::AbstractQuasiVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()

# like empty, but should return a mutable collection, a Vector by default
emptymutable(a::AbstractQuasiVector{T}, ::Type{U}=T) where {T,U} = Vector{U}()

## copy between abstract arrays - generally more efficient
## since a single index variable can be used.

copyto!(dest::AbstractQuasiArray, src::AbstractQuasiArray) =
    copyto!(IndexStyle(dest), dest, IndexStyle(src), src)

function copyto!(::IndexStyle, dest::AbstractQuasiArray, ::IndexStyle, src::AbstractQuasiArray)
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    isempty(srcinds) || (checkbounds(Bool, destinds, first(srcinds)) && checkbounds(Bool, destinds, last(srcinds))) ||
        throw(BoundsError(dest, srcinds))
    @inbounds for i in srcinds
        dest[i] = src[i]
    end
    return dest
end

function copyto!(::IndexStyle, dest::AbstractQuasiArray, ::IndexCartesian, src::AbstractQuasiArray)
    axes(dest) == axes(src) || throw(BoundsError(dest, axes(src)))
    i = 0
    @inbounds for i in eachindex(src)
        dest[i] = src[i]
    end
    return dest
end

function copy(a::AbstractQuasiArray)
    @_propagate_inbounds_meta
    copymutable(a)
end

function copymutable(a::AbstractQuasiArray)
    @_propagate_inbounds_meta
    copyto!(similar(a), a)
end

## iteration support for arrays by iterating over `eachindex` in the array ##
# Allows fast iteration by default for both IndexLinear and IndexCartesian arrays

# While the definitions for IndexLinear are all simple enough to inline on their
# own, IndexCartesian's CartesianIndices is more complicated and requires explicit
# inlining.
function iterate(A::AbstractQuasiArray, state=(eachindex(A),))
    y = iterate(state...)
    y === nothing && return nothing
    A[y[1]], (state[1], tail(y)...)
end

isempty(a::AbstractQuasiArray) = (length(a) == 0)

getindex(A::AbstractQuasiArray, I...) = _getindex(indextype(A), A, I)

function _getindex(::Type{IND}, A::AbstractQuasiArray, I) where IND
    @_propagate_inbounds_meta
    error_if_canonical_getindex(IndexStyle(A), A, I)
    _getindex(IND, IndexStyle(A), A, to_indices(A, I))
end
function unsafe_getindex(A::AbstractQuasiArray, I...)
    @_inline_meta
    @inbounds r = getindex(A, I...)
    r
end

## Internal definitions
_getindex(_, ::IndexStyle, A::AbstractQuasiArray, I) = layout_getindex(A, I...)

## IndexCartesian Scalar indexing: Canonical method is full dimensionality of indices
function _getindex(::Type{IND}, ::IndexCartesian, A::AbstractQuasiArray, I::IND) where {M,IND}
    @_inline_meta
    @boundscheck checkbounds(A, I...) # generally _to_subscript_indices requires bounds checking
    @inbounds r = getindex(A, _to_subscript_indices(A, I...)...)
    r
end

error_if_canonical_getindex(::IndexCartesian, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_getindex(indextype(A), A, I)

error_if_canonical_getindex(::IndexLinear, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_getindex(indextype(A), A, I)

_error_if_canonical_getindex(::Type{IND}, A::AbstractQuasiArray{T,N}, I::IND) where {T,N,IND} =
    error("getindex not defined for ", typeof(A))

_error_if_canonical_getindex(::Type, ::AbstractQuasiArray, ::Any...) = nothing

_to_subscript_indices(A::AbstractQuasiArray, i) = (@_inline_meta; _unsafe_ind2sub(A, i))
_to_subscript_indices(A::AbstractQuasiArray{T,N}) where {T,N} = (@_inline_meta; fill_to_length((), 1, Val(N)))
_to_subscript_indices(A::AbstractQuasiArray{T,0}) where {T} = ()
_to_subscript_indices(A::AbstractQuasiArray{T,0}, i) where {T} = ()
_to_subscript_indices(A::AbstractQuasiArray{T,0}, I...) where {T} = ()
function _to_subscript_indices(A::AbstractQuasiArray{T,N}, I...) where {T,N}
    @_inline_meta
    J, Jrem = IteratorsMD.split(I, Val(N))
    _to_subscript_indices(A, J, Jrem)
end
_to_subscript_indices(A::AbstractQuasiArray, J::Tuple, Jrem::Tuple{}) =
    __to_subscript_indices(A, axes(A), J, Jrem)
function __to_subscript_indices(A::AbstractQuasiArray,
        ::Tuple{AbstractUnitRange,Vararg{AbstractUnitRange}}, J::Tuple, Jrem::Tuple{})
    @_inline_meta
    (J..., map(first, tail(_remaining_size(J, axes(A))))...)
end
_to_subscript_indices(A::AbstractQuasiArray{T,N}, I::Vararg{Any,N}) where {T,N} = I

## Setindex! is defined similarly. We first dispatch to an internal _setindex!
# function that allows dispatch on array storage

setindex!(A::AbstractQuasiArray, v, I...) = _setindex!(indextype(A), A, v, I)

function _setindex!(::Type{IND}, A::AbstractQuasiArray, v, I) where IND
    @_propagate_inbounds_meta
    error_if_canonical_setindex(IndexStyle(A), A, I)
    _setindex!(IND, IndexStyle(A), A, v, to_indices(A, I))
end

error_if_canonical_setindex(::IndexCartesian, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_setindex(indextype(A), A, I)

error_if_canonical_setindex(::IndexLinear, A::AbstractQuasiArray{T,N}, I::Tuple) where {T,N} =
    _error_if_canonical_setindex(indextype(A), A, I)

_error_if_canonical_setindex(::Type{IND}, A::AbstractQuasiArray{T,N}, I::IND) where {T,N,IND} =
    error("setindex! not defined for ", typeof(A))

_error_if_canonical_setindex(::Type, ::AbstractQuasiArray, ::Any...) = nothing
## Internal definitions
_setindex!(::Type, ::IndexStyle, A::AbstractQuasiArray, v, I) =
    error("setindex! for $(typeof(A)) with types $(typeof(I)) is not supported")

# IndexCartesian Scalar indexing
function _setindex!(::Type{IND}, ::IndexCartesian, A::AbstractQuasiArray{T,N}, v, I::NTuple{N,Any}) where {T,N,IND}
    @_propagate_inbounds_meta
    setindex!(A, v, I...)
end
function _setindex!(::Type{IND}, ::IndexCartesian, A::AbstractQuasiArray, v, I::NTuple{M,Any}) where {M,IND}
    @_inline_meta
    @boundscheck checkbounds(A, I...)
    @inbounds r = setindex!(A, v, _to_subscript_indices(IND, A, I...)...)
    r
end

parent(a::AbstractQuasiArray) = a

unalias(dest, A::AbstractQuasiArray) = mightalias(dest, A) ? unaliascopy(A) : A
unaliascopy(A::AbstractQuasiArray)::typeof(A) = (@_noinline_meta; _unaliascopy(A, copy(A)))

"""
    Base.mightalias(A::AbstractQuasiArray, B::AbstractQuasiArray)

Perform a conservative test to check if quasi arrays `A` and `B` might share the same memory.

By default, this simply checks if either of the arrays reference the same memory
regions, as identified by their [`Base.dataids`](@ref).
"""
mightalias(A::AbstractQuasiArray, B::AbstractQuasiArray) = !_isdisjoint(dataids(A), dataids(B))

"""
    Base.dataids(A::AbstractQuasiArray)

Return a tuple of `UInt`s that represent the mutable data segments of an array.

Custom arrays that would like to opt-in to aliasing detection of their component
parts can specialize this method to return the concatenation of the `dataids` of
their component parts.  A typical definition for an array that wraps a parent is
`Base.dataids(C::CustomArray) = dataids(C.parent)`.
"""
dataids(A::AbstractQuasiArray) = (UInt(objectid(A)),)

## Reductions and accumulates ##

function isequal(A::AbstractQuasiArray, B::AbstractQuasiArray)
    if A === B return true end
    if axes(A) != axes(B)
        return false
    end
    for (a, b) in zip(A, B)
        if !isequal(a, b)
            return false
        end
    end
    return true
end

function _equals(_, _, A, B)
    if axes(A) != axes(B)
        return false
    end
    anymissing = false
    for (a, b) in zip(A, B)
        eq = (a == b)
        if ismissing(eq)
            anymissing = true
        elseif !eq
            return false
        end
    end
    return anymissing ? missing : true
end

(==)(A::AbstractQuasiArray, B::AbstractQuasiArray) = _equals(MemoryLayout(A), MemoryLayout(B), A, B)

##
# show
##

show(io::IO, A::AbstractQuasiArray) = summary(io, A)

struct QuasiArrayLayout <: MemoryLayout end
MemoryLayout(::Type{<:AbstractQuasiArray}) = QuasiArrayLayout()


text:
using DataValues
using Dates

a, b = DataValue("14:00:00"), DataValue("15:15:00")

Dates.value(DateTime(a,"HH:MM:SS") - DateTime(b,"HH:MM:SS"))

c = DateTime(a,"HH:MM:SS") - DateTime(b,"HH:MM:SS")


text:
struct SymEGRQSCholesky{T,UT<:AbstractMatrix,WT<:AbstractMatrix,dT<:AbstractArray} <: AbstractMatrix{T}
	Ut::UT
    Wt::WT
	d::dT
    n::Int
    p::Int
    function SymEGRQSCholesky{T,UT,WT,dT}(Ut,Wt,d,n,p) where
		{T,UT<:AbstractMatrix,WT<:AbstractMatrix,dT<:AbstractArray}
		Up, Un = size(Ut)
		Wp, Wn = size(Wt)
		(Up == Wp && Un == Wn && Un == length(d)) || throw(DimensionMismatch())
        new(Ut,Wt,d,n,p)
    end
end

#### Creating W and dbar s.t. L = tril(UW',-1) + diag(dbar) ####
function dss_create_wdbar(Ut::AbstractArray, Vt::AbstractArray, d::AbstractArray)
    p, n = size(Ut);
    P  = zeros(p, p);
    Wt  = similar(Vt);
    dbar = zeros(n);
    @inbounds for i = 1:n
        tmpU  = Ut[:,i]
        tmpW  = Vt[:,i] - P*tmpU;
		tmpds = sqrt(abs(tmpU'*tmpW + d[i]));
        #tmpds = sqrt(tmpU'*tmpW + d[i]);
        tmpW  = tmpW/tmpds
        Wt[:,i] = tmpW;
        dbar[i] = tmpds;
        P += tmpW*tmpW';
    end
    return Wt, dbar
end

function cholesky(K::SymEGRQSMatrix{T,UT,VT,dT}) where
	{T,UT<:AbstractMatrix,VT<:AbstractMatrix,dT<:AbstractArray}
	Wt, dbar = dss_create_wdbar(K.Ut,K.Vt,K.d)
	SymEGRQSCholesky{T,typeof(K.Ut),typeof(Wt),typeof(dbar)}(K.Ut,Wt,dbar,K.n,K.p)
end

function cholesky(K::SymEGRSSMatrix{T,UT,VT}, σ::Number) where
	{T,UT<:AbstractMatrix,VT<:AbstractMatrix}
	Wt, dbar = dss_create_wdbar(K.Ut,K.Vt,ones(K.n)*σ)
	SymEGRQSCholesky{T,typeof(K.Ut),typeof(Wt),typeof(dbar)}(K.Ut,Wt,dbar,K.n,K.p)
end

function SymEGRQSCholesky(K::SymEGRSSMatrix{T,UT,VT}, σ::Number) where
	{T,UT<:AbstractMatrix,VT<:AbstractMatrix}
	Wt, dbar = dss_create_wdbar(K.V,K.V,ones(K.n)*σ)
	SymEGRQSCholesky{T,typeof(K.Ut),typeof(Wt),typeof(dbar)}(K.Ut,Wt,dbar,K.n,K.p)
end

function SymEGRQSCholesky(Ut::UT, Wt::WT, d::AbstractArray) where
	{T,UT<:AbstractMatrix,WT<:AbstractMatrix}
	p, n = size(Ut);
	SymEGRQSCholesky{eltype(Ut),typeof(Ut),typeof(Wt),typeof(d)}(Ut,Wt,d,n,p)
end

########################################################################
#### Helpful properties. Not nessecarily computionally efficient    ####
########################################################################
Matrix(K::SymEGRQSCholesky) = getproperty(K,:L)
size(K::SymEGRQSCholesky) = (K.n, K.n)
size(K::SymEGRQSCholesky,d::Int) = (1 <= d && d <=2) ? size(K)[d] : throw(ArgumentError("Invalid dimension $d"))

function getindex(K::SymEGRQSCholesky, i::Int, j::Int)
	i > j && return dot(K.Ut[:,i], K.Wt[:,j])
	i == j && return K.d[i]
	return 0
end

function getproperty(K::SymEGRQSCholesky, d::Symbol)
    if d === :U
        return UpperTriangular(triu(K.Wt'*K.Ut,1) + Diagonal(K.d))
    elseif d === :L
        return LowerTriangular(tril(K.Ut'*K.Wt,-1) + Diagonal(K.d))
    else
        return getfield(K, d)
    end
end

Base.propertynames(F::SymEGRQSCholesky, private::Bool=false) =
    (:U, :L, (private ? fieldnames(typeof(F)) : ())...)

function Base.show(io::IO, mime::MIME{Symbol("text/plain")},
		 K::SymEGRQSCholesky{<:Any,<:AbstractArray,<:AbstractArray,<:AbstractArray})
    summary(io, K); println(io)
    show(io, mime, K.L)
end

########################################################################
#### Cholesky factorization of Higher-order quasiseparable matrices ####
########################################################################

#### Forward substitution ####
function dss_forward!(X::AbstractArray,Ut::AbstractArray,Wt::AbstractArray,
                     ds::AbstractArray,B::AbstractArray)
    p, n = size(Ut)
    mx = size(B,2)
    Wbar = zeros(p,mx);
    @inbounds for i = 1:n
        tmpU = Ut[:,i];
        tmpW = Wt[:,i];
        X[i:i,:] = (B[i:i,:] - tmpU'*Wbar)/ds[i];
        Wbar += tmpW .* X[i:i,:];
    end
end
#### Backward substitution ####
function dssa_backward!(X::AbstractArray,Ut::AbstractArray,Wt::AbstractArray,
                       ds::AbstractArray,B::AbstractArray)
    p, n = size(Ut)
    mx = size(B,2)
    Ubar = zeros(p,mx);
    @inbounds for i = n:-1:1
        tmpU = Ut[:,i];
        tmpW = Wt[:,i];
        X[i:i,:] = (B[i:i,:] - tmpW'*Ubar)/ds[i];
        Ubar += tmpU .* X[i:i,:];
    end
end

#### Multiplying with L ####
function dss_tri_mul!(Y::AbstractArray,K::SymEGRQSCholesky,X::AbstractArray)
     p, n = size(K.Ut)
     mx = size(X, 2)
     Wbar = zeros(p, mx)
     @inbounds for i = 1:n
         tmpW = K.Wt[:,i]
         tmpU = K.Ut[:,i]
         tmpX = X[i:i,:];
         Y[i,:] = tmpU'*Wbar + K.d[i]*tmpX;
         Wbar  +=  tmpW*tmpX;
     end
	 return Y
end

#### Multiplying with L' ####
function dssa_tri_mul!(Y::AbstractArray,K::SymEGRQSCholesky,X::AbstractArray)
     p, n = size(K.Ut)
     mx = size(X, 2)
     Ubar = zeros(p,mx)
     @inbounds for i = n:-1:1
         tmpW = K.Wt[:,i]
         tmpU = K.Ut[:,i]
         tmpX = X[i:i,:]
         Y[i,:] = tmpW'*Ubar + K.d[i]*tmpX;
         Ubar = Ubar + tmpU*tmpX;
     end
	 return Y
end

#### Squared norm of columns of L = tril(UW',-1) + diag(dbar) ####
function squared_norm_cols(Ut::AbstractArray,Wt::AbstractArray,
                         dbar::AbstractArray)
    p, n = size(Ut)
    P = zeros(p, p)
    c = zeros(n)
    @inbounds for i = n:-1:1
        tmpW = Wt[:,i]
        tmpU = Ut[:,i]
        c[i] = dbar[i]^2 + tmpW'*P*tmpW
        P += tmpU*tmpU'
    end
    return c
end
#### Implicit inverse of  L = tril(UW',-1) + diag(dbar) ####
function dss_create_yz(Ut::AbstractArray, Wt::AbstractArray,
                     dbar::AbstractArray)
    p, n = size(Ut)
    Y = zeros(n,p)
    Z = zeros(n,p)
    dss_forward!(Y, Ut, Wt, dbar, Ut')
    dssa_backward!(Z, Ut, Wt, dbar, Wt')
    # Probably best not to use inv
    return copy(Y'), copy((Z*inv(Ut*Z - Diagonal(ones(p))))')
	#return Y, Z*((U'*Z - Diagonal(ones(m)))\Diagonal(ones(m)))
end

function fro_norm_L(L::SymEGRQSCholesky)
	return sum(squared_norm_cols(L.Ut, L.Wt, L.d))
end

function trinv(L::SymEGRQSCholesky)
	Yt, Zt = dss_create_yz(L.Ut, L.Wt, L.d)
	return sum(squared_norm_cols(Yt, Zt, L.d.^(-1)))
end

function tr(Ky::SymEGRQSCholesky, K::SymEGRSSMatrix)
	n = Ky.n;
	p = Ky.p;
	c = Ky.d;
	Ut = K.Ut;
	Vt = K.Vt;
	Yt, Zt = dss_create_yz(Ky.Ut, Ky.Wt, Ky.d);
	b = 0;
	P = zeros(p,p);
	R = zeros(p,p);
	@inbounds for k = 1:Ky.n
		yk = Yt[:,k];
		zk = Zt[:,k];
		uk = Ut[:,k];
		vk = Vt[:,k];
		cki = c[k]^(-1);
		b += yk'*P*yk + 2*yk'*R*uk*cki + uk'*vk*(cki^2);
		P += ((uk'*vk)*zk)*zk' + zk*(R*uk)' + (R*uk)*zk';
		R += zk*vk';
	end
	return b
end

#### Log-determinant ####
function logdet(K::SymEGRQSCholesky)
	dd = sum(log.(K.d))
    return dd
end

#### Determinant ####
function det(K::SymEGRQSCholesky)
    return exp(logdet(K))
end

mul!(y::AbstractVecOrMat, K::SymEGRQSCholesky, x::AbstractVecOrMat) =
	dss_tri_mul!(y,K,x)
mul!(y::AbstractVecOrMat, K::Adjoint{<:Any,<:SymEGRQSCholesky}, x::AbstractVecOrMat) =
	dssa_tri_mul!(y,K.parent,x)
function (\)(F::SymEGRQSCholesky, B::AbstractVecOrMat)
	X = similar(B)
	dss_forward!(X,F.Ut,F.Wt,F.d,B)
	return X
end

function (\)(F::Adjoint{<:Any,<:SymEGRQSCholesky}, B::AbstractVecOrMat)
	Y = similar(B)
	dssa_backward!(Y,F.parent.Ut,F.parent.Wt,F.parent.d,B)
	return Y
end

max_stars_repo_name: mipals/SymEGRSSMatrices


text:
using Distributions

"""
General Interface for Computation of confidence interval of price

		IC=confinter(mcProcess,rfCurve,mcBaseData,payoff_,alpha)

Where:\n
		mcProcess          = Process to be simulated.
		rfCurve  = Zero Rate Data.
		mcBaseData = Basic properties of MonteCarlo simulation
		payoff_ = Payoff(s) to be priced
		alpha [Optional, default to 99%] = confidence level

		Price     = Price of the derivative

"""
function confinter(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoff::AbstractPayoff, alpha::Real = 0.99)
    set_seed(mcConfig)
    T = maturity(abstractPayoff)
    Nsim = mcConfig.Nsim
    S = simulate(mcProcess, rfCurve, mcConfig, T)
    Payoff = payoff(S, abstractPayoff, rfCurve)
    mean1 = mean(Payoff)
    var1 = var(Payoff)
    alpha_ = 1 - alpha
    dist1 = Distributions.TDist(Nsim)
    tstar = quantile(dist1, 1 - alpha_ / 2.0)
    IC = (mean1 - tstar * sqrt(var1) / Nsim, mean1 + tstar * sqrt(var1) / Nsim)
    return IC
end

function confinter(mcProcess::BaseProcess, rfCurve::AbstractZeroRateCurve, mcConfig::MonteCarloConfiguration, abstractPayoffs::Array{abstractPayoff_}, alpha::Real = 0.99) where {abstractPayoff_ <: AbstractPayoff}
    set_seed(mcConfig)
    maxT = maximum([maturity(abstractPayoff) for abstractPayoff in abstractPayoffs])
    Nsim = mcConfig.Nsim
    S = simulate(mcProcess, rfCurve, mcConfig, maxT)
    Means = [mean(payoff(S, abstractPayoff, rfCurve, maxT)) for abstractPayoff in abstractPayoffs]
    Vars = [var(payoff(S, abstractPayoff, rfCurve, maxT)) for abstractPayoff in abstractPayoffs]
    alpha_ = 1 - alpha
    dist1 = Distributions.TDist(Nsim)
    tstar = quantile(dist1, 1 - alpha_ / 2.0)
    IC = [(mean1 - tstar * sqrt(var1) / Nsim, mean1 + tstar * sqrt(var1) / Nsim) for (mean1, var1) in zip(Means, Vars)]
    return IC
end


text:
@doc raw"""
    SpecialEuclidean(n)

Special Euclidean group $\mathrm{SE}(n)$, the group of rigid motions.

``\mathrm{SE}(n)`` is the semidirect product of the [`TranslationGroup`](@ref) on $ℝ^n$ and
[`SpecialOrthogonal(n)`](@ref)

````math
\mathrm{SE}(n) ≐ \mathrm{T}(n) ⋊_θ \mathrm{SO}(n),
````

where $θ$ is the canonical action of $\mathrm{SO}(n)$ on $\mathrm{T}(n)$ by vector rotation.

This constructor is equivalent to calling

```julia
Tn = TranslationGroup(n)
SOn = SpecialOrthogonal(n)
SemidirectProductGroup(Tn, SOn, RotationAction(Tn, SOn))
```

Points on $\mathrm{SE}(n)$ may be represented as points on the underlying product manifold
$\mathrm{T}(n) × \mathrm{SO}(n)$. For group-specific functions, they may also be
represented as affine matrices with size `(n + 1, n + 1)` (see [`affine_matrix`](@ref)), for
which the group operation is [`MultiplicationOperation`](@ref).
"""
const SpecialEuclidean{N} = SemidirectProductGroup{
    ℝ,
    TranslationGroup{Tuple{N},ℝ},
    SpecialOrthogonal{N},
    RotationAction{TranslationGroup{Tuple{N},ℝ},SpecialOrthogonal{N},LeftAction},
}

function SpecialEuclidean(n)
    Tn = TranslationGroup(n)
    SOn = SpecialOrthogonal(n)
    A = RotationAction(Tn, SOn)
    return SemidirectProductGroup(Tn, SOn, A)
end

Base.show(io::IO, ::SpecialEuclidean{n}) where {n} = print(io, "SpecialEuclidean($(n))")

Base.@propagate_inbounds function submanifold_component(
    ::SpecialEuclidean{n},
    p::AbstractMatrix,
    ::Val{1},
) where {n}
    return view(p, 1:n, n + 1)
end
Base.@propagate_inbounds function submanifold_component(
    ::SpecialEuclidean{n},
    p::AbstractMatrix,
    ::Val{2},
) where {n}
    return view(p, 1:n, 1:n)
end

function submanifold_components(G::SpecialEuclidean{n}, p::AbstractMatrix) where {n}
    @assert size(p) == (n + 1, n + 1)
    @inbounds t = submanifold_component(G, p, Val(1))
    @inbounds R = submanifold_component(G, p, Val(2))
    return (t, R)
end

Base.@propagate_inbounds function _padpoint!(
    ::SpecialEuclidean{n},
    q::AbstractMatrix,
) where {n}
    for i in 1:n
        q[n + 1, i] = 0
    end
    q[n + 1, n + 1] = 1
    return q
end

Base.@propagate_inbounds function _padvector!(
    ::SpecialEuclidean{n},
    X::AbstractMatrix,
) where {n}
    for i in 1:(n + 1)
        X[n + 1, i] = 0
    end
    return X
end

@doc raw"""
    adjoint_action(::SpecialEuclidean{3}, p, fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}})

Adjoint action of the [`SpecialEuclidean`](@ref) group on the vector with coefficients `fX`
tangent at point `p`.

The formula for the coefficients reads ``t×(R⋅ω) + R⋅r`` for the translation part and
``R⋅ω`` for the rotation part, where `t` is the translation part of `p`, `R` is the rotation
matrix part of `p`, `r` is the translation part of `fX` and `ω` is the rotation part of `fX`,
``×`` is the cross product and ``⋅`` is the matrix product.
"""
function adjoint_action(::SpecialEuclidean{3}, p, fX::TFVector{<:Any,VeeOrthogonalBasis{ℝ}})
    t = p.parts[1]
    R = p.parts[2]
    r = fX.data[SA[1, 2, 3]]
    ω = fX.data[SA[4, 5, 6]]
    Rω = R * ω
    return TFVector([cross(t, Rω) + R * r; Rω], fX.basis)
end

@doc raw"""
    affine_matrix(G::SpecialEuclidean, p) -> AbstractMatrix

Represent the point $p ∈ \mathrm{SE}(n)$ as an affine matrix.
For $p = (t, R) ∈ \mathrm{SE}(n)$, where $t ∈ \mathrm{T}(n), R ∈ \mathrm{SO}(n)$, the
affine representation is the $n + 1 × n + 1$ matrix

````math
\begin{pmatrix}
R & t \\
0^\mathrm{T} & 1
\end{pmatrix}.
````

This function embeds $\mathrm{SE}(n)$ in the general linear group $\mathrm{GL}(n+1)$.
It is an isometric embedding and group homomorphism [^RicoMartinez1988].

See also [`screw_matrix`](@ref) for matrix representations of the Lie algebra.

[^RicoMartinez1988]:
    > Rico Martinez, J. M., “Representations of the Euclidean group and its applications
    > to the kinematics of spatial chains,” PhD Thesis, University of Florida, 1988.
"""
function affine_matrix(G::SpecialEuclidean{n}, p) where {n}
    pis = submanifold_components(G, p)
    pmat = allocate_result(G, affine_matrix, pis...)
    map(copyto!, submanifold_components(G, pmat), pis)
    @inbounds _padpoint!(G, pmat)
    return pmat
end
affine_matrix(::SpecialEuclidean{n}, p::AbstractMatrix) where {n} = p
function affine_matrix(::GT, ::Identity{GT}) where {n,GT<:SpecialEuclidean{n}}
    s = maybesize(Size(n, n))
    s isa Size && return SDiagonal{n,Float64}(I)
    return Diagonal{Float64}(I, n)
end

@doc raw"""
    screw_matrix(G::SpecialEuclidean, X) -> AbstractMatrix

Represent the Lie algebra element $X ∈ 𝔰𝔢(n) = T_e \mathrm{SE}(n)$ as a screw matrix.
For $X = (b, Ω) ∈ 𝔰𝔢(n)$, where $Ω ∈ 𝔰𝔬(n) = T_e \mathrm{SO}(n)$, the screw representation is
the $n + 1 × n + 1$ matrix

````math
\begin{pmatrix}
Ω & b \\
0^\mathrm{T} & 0
\end{pmatrix}.
````

This function embeds $𝔰𝔢(n)$ in the general linear Lie algebra $𝔤𝔩(n+1)$ but it's not
a homomorphic embedding (see [`SpecialEuclideanInGeneralLinear`](@ref) for a homomorphic one).

See also [`affine_matrix`](@ref) for matrix representations of the Lie group.
"""
function screw_matrix(G::SpecialEuclidean{n}, X) where {n}
    Xis = submanifold_components(G, X)
    Xmat = allocate_result(G, screw_matrix, Xis...)
    map(copyto!, submanifold_components(G, Xmat), Xis)
    @inbounds _padvector!(G, Xmat)
    return Xmat
end
screw_matrix(::SpecialEuclidean{n}, X::AbstractMatrix) where {n} = X

function allocate_result(G::SpecialEuclidean{n}, f::typeof(affine_matrix), p...) where {n}
    return allocate(p[1], Size(n + 1, n + 1))
end
function allocate_result(G::SpecialEuclidean{n}, f::typeof(screw_matrix), X...) where {n}
    return allocate(X[1], Size(n + 1, n + 1))
end

compose(::SpecialEuclidean, p::AbstractMatrix, q::AbstractMatrix) = p * q

function compose!(
    ::SpecialEuclidean,
    x::AbstractMatrix,
    p::AbstractMatrix,
    q::AbstractMatrix,
)
    return mul!(x, p, q)
end

@doc raw"""
    group_exp(G::SpecialEuclidean{n}, X)

Compute the group exponential of $X = (b, Ω) ∈ 𝔰𝔢(n)$, where $b ∈ 𝔱(n)$ and $Ω ∈ 𝔰𝔬(n)$:

````math
\exp X = (t, R),
````

where $t ∈ \mathrm{T}(n)$ and $R = \exp Ω$ is the group exponential on $\mathrm{SO}(n)$.

In the [`screw_matrix`](@ref) representation, the group exponential is the matrix
exponential (see [`group_exp`](@ref)).
"""
group_exp(::SpecialEuclidean, ::Any)

@doc raw"""
    group_exp(G::SpecialEuclidean{2}, X)

Compute the group exponential of $X = (b, Ω) ∈ 𝔰𝔢(2)$, where $b ∈ 𝔱(2)$ and $Ω ∈ 𝔰𝔬(2)$:

````math
\exp X = (t, R) = (U(θ) b, \exp Ω),
````

where $t ∈ \mathrm{T}(2)$, $R = \exp Ω$ is the group exponential on $\mathrm{SO}(2)$,

````math
U(θ) = \frac{\sin θ}{θ} I_2 + \frac{1 - \cos θ}{θ^2} Ω,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_exp(::SpecialEuclidean{2}, ::Any)

@doc raw"""
    group_exp(G::SpecialEuclidean{3}, X)

Compute the group exponential of $X = (b, Ω) ∈ 𝔰𝔢(3)$, where $b ∈ 𝔱(3)$ and $Ω ∈ 𝔰𝔬(3)$:

````math
\exp X = (t, R) = (U(θ) b, \exp Ω),
````

where $t ∈ \mathrm{T}(3)$, $R = \exp Ω$ is the group exponential on $\mathrm{SO}(3)$,

````math
U(θ) = I_3 + \frac{1 - \cos θ}{θ^2} Ω + \frac{θ - \sin θ}{θ^3} Ω^2,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_exp(::SpecialEuclidean{3}, ::Any)

function group_exp!(G::SpecialEuclidean, q, X)
    Xmat = screw_matrix(G, X)
    qmat = exp(Xmat)
    map(copyto!, submanifold_components(G, q), submanifold_components(G, qmat))
    _padpoint!(G, q)
    return q
end
function group_exp!(G::SpecialEuclidean{2}, q, X)
    SO2 = submanifold(G, 2)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(R) == (2, 2)
    @assert size(t) == (2,)
    @assert size(b) == (2,)

    θ = vee(SO2, Identity(SO2, Ω), Ω)[1]
    sinθ, cosθ = sincos(θ)
    if θ ≈ 0
        α = 1 - θ^2 / 6
        β = θ / 2
    else
        α = sinθ / θ
        β = (1 - cosθ) / θ
    end

    @inbounds begin
        R[1] = cosθ
        R[2] = sinθ
        R[3] = -sinθ
        R[4] = cosθ
        t[1] = α * b[1] - β * b[2]
        t[2] = α * b[2] + β * b[1]
        _padpoint!(G, q)
    end
    return q
end
function group_exp!(G::SpecialEuclidean{3}, q, X)
    SO3 = submanifold(G, 2)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(R) == (3, 3)
    @assert size(t) == (3,)

    θ = norm(SO3, Identity(SO3, Ω), Ω) / sqrt(2)
    θ² = θ^2
    if θ ≈ 0
        α = 1 - θ² / 6
        β = θ / 2
        γ = 1 / 6 - θ² / 120
    else
        sinθ, cosθ = sincos(θ)
        α = sinθ / θ
        β = (1 - cosθ) / θ²
        γ = (1 - α) / θ²
    end

    Ω² = Ω^2
    Jₗ = I + β .* Ω .+ γ .* Ω²
    R .= I + α .* Ω .+ β .* Ω²
    copyto!(t, Jₗ * b)
    @inbounds _padpoint!(G, q)
    return q
end

@doc raw"""
    group_log(G::SpecialEuclidean{n}, p) where {n}

Compute the group logarithm of $p = (t, R) ∈ \mathrm{SE}(n)$, where $t ∈ \mathrm{T}(n)$
and $R ∈ \mathrm{SO}(n)$:

````math
\log p = (b, Ω),
````

where $b ∈ 𝔱(n)$ and $Ω = \log R ∈ 𝔰𝔬(n)$ is the group logarithm on $\mathrm{SO}(n)$.

In the [`affine_matrix`](@ref) representation, the group logarithm is the matrix logarithm
(see [`group_log`](@ref)):
"""
group_log(::SpecialEuclidean, ::Any)

@doc raw"""
    group_log(G::SpecialEuclidean{2}, p)

Compute the group logarithm of $p = (t, R) ∈ \mathrm{SE}(2)$, where $t ∈ \mathrm{T}(2)$
and $R ∈ \mathrm{SO}(2)$:

````math
\log p = (b, Ω) = (U(θ)^{-1} t, \log R),
````

where $b ∈ 𝔱(2)$, $Ω = \log R ∈ 𝔰𝔬(2)$ is the group logarithm on $\mathrm{SO}(2)$,

````math
U(θ) = \frac{\sin θ}{θ} I_2 + \frac{1 - \cos θ}{θ^2} Ω,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_log(::SpecialEuclidean{2}, ::Any)

@doc raw"""
    group_log(G::SpecialEuclidean{3}, p)

Compute the group logarithm of $p = (t, R) ∈ \mathrm{SE}(3)$, where $t ∈ \mathrm{T}(3)$
and $R ∈ \mathrm{SO}(3)$:

````math
\log p = (b, Ω) = (U(θ)^{-1} t, \log R),
````

where $b ∈ 𝔱(3)$, $Ω = \log R ∈ 𝔰𝔬(3)$ is the group logarithm on $\mathrm{SO}(3)$,

````math
U(θ) = I_3 + \frac{1 - \cos θ}{θ^2} Ω + \frac{θ - \sin θ}{θ^3} Ω^2,
````

and $θ = \frac{1}{\sqrt{2}} \lVert Ω \rVert_e$
(see [`norm`](@ref norm(M::Rotations, p, X))) is the angle of the rotation.
"""
group_log(::SpecialEuclidean{3}, ::Any)

function group_log!(G::SpecialEuclidean, X, q)
    qmat = affine_matrix(G, q)
    Xmat = real(log_safe(qmat))
    map(copyto!, submanifold_components(G, X), submanifold_components(G, Xmat))
    _padvector!(G, X)
    return X
end
function group_log!(G::SpecialEuclidean{2}, X, q)
    SO2 = submanifold(G, 2)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(b) == (2,)

    group_log!(SO2, Ω, R)
    @inbounds θ = Ω[2]
    β = θ / 2
    α = θ ≈ 0 ? 1 - β^2 / 3 : β * cot(β)

    @inbounds begin
        b[1] = α * t[1] + β * t[2]
        b[2] = α * t[2] - β * t[1]
        _padvector!(G, X)
    end
    return X
end
function group_log!(G::SpecialEuclidean{3}, X, q)
    b, Ω = submanifold_components(G, X)
    t, R = submanifold_components(G, q)
    @assert size(Ω) == (3, 3)
    @assert size(b) == (3,)

    trR = tr(R)
    cosθ = (trR - 1) / 2
    θ = acos(clamp(cosθ, -1, 1))
    θ² = θ^2
    if θ ≈ 0
        α = 1 / 2 + θ² / 12
        β = 1 / 12 + θ² / 720
    else
        sinθ = sin(θ)
        α = θ / sinθ / 2
        β = 1 / θ² - (1 + cosθ) / 2 / θ / sinθ
    end

    Ω .= (R .- transpose(R)) .* α
    Jₗ⁻¹ = I - Ω ./ 2 .+ β .* Ω^2
    mul!(b, Jₗ⁻¹, t)
    @inbounds _padvector!(G, X)
    return X
end

"""
    lie_bracket(G::SpecialEuclidean, X::ProductRepr, Y::ProductRepr)
    lie_bracket(G::SpecialEuclidean, X::AbstractMatrix, Y::AbstractMatrix)

Calculate the Lie bracket between elements `X` and `Y` of the special Euclidean Lie
algebra. For the matrix representation (which can be obtained using [`screw_matrix`](@ref))
the formula is ``[X, Y] = XY-YX``, while in the [`ProductRepr`](@ref) representation the
formula reads ``[X, Y] = [(t_1, R_1), (t_2, R_2)] = (R_1 t_2 - R_2 t_1, R_1 R_2 - R_2 R_1)``.
"""
function lie_bracket(G::SpecialEuclidean, X::ProductRepr, Y::ProductRepr)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    return ProductRepr(hX * nY - hY * nX, lie_bracket(G.manifold.manifolds[2], hX, hY))
end
function lie_bracket(::SpecialEuclidean, X::AbstractMatrix, Y::AbstractMatrix)
    return X * Y - Y * X
end

function lie_bracket!(G::SpecialEuclidean, Z, X, Y)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    nZ, hZ = submanifold_components(G, Z)
    lie_bracket!(G.manifold.manifolds[2], hZ, hX, hY)
    nZ .= hX * nY .- hY * nX
    @inbounds _padvector!(G, Z)
    return Z
end

"""
    translate_diff(G::SpecialEuclidean, p, q, X, ::RightAction)

Differential of the right action of the [`SpecialEuclidean`](@ref) group on itself.
The formula for the rotation part is the differential of the right rotation action, while
the formula for the translation part reads
````math
R_q⋅X_R⋅t_p + X_t
````
where ``R_q`` is the rotation part of `q`, ``X_R`` is the rotation part of `X`, ``t_p``
is the translation part of `p` and ``X_t`` is the translation part of `X`.
"""
translate_diff(G::SpecialEuclidean, p, q, X, ::RightAction)

function translate_diff!(G::SpecialEuclidean, Y, p, q, X, ::RightAction)
    np, hp = submanifold_components(G, p)
    nq, hq = submanifold_components(G, q)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    hY .= hp' * hX * hp
    copyto!(nY, hq * (hX * np) + nX)
    @inbounds _padvector!(G, Y)
    return Y
end
function translate_diff!(G::SpecialEuclidean, Y, ::Identity, q, X, ::RightAction)
    copyto!(G, Y, X)
    @inbounds _padvector!(G, Y)
    return Y
end
function translate_diff!(G::SpecialEuclidean, Y, p, ::Identity, X, ::RightAction)
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    nY, hY = submanifold_components(G, Y)
    hY .= hp' * hX * hp
    copyto!(nY, hX * np + nX)
    @inbounds _padvector!(G, Y)
end
function translate_diff!(G::SpecialEuclidean, Y, ::Identity, ::Identity, X, ::RightAction)
    copyto!(G, Y, X)
    @inbounds _padvector!(G, Y)
    return Y
end

@doc raw"""
    SpecialEuclideanInGeneralLinear

An explicit isometric and homomorphic embedding of $\mathrm{SE}(n)$ in $\mathrm{GL}(n+1)$
and $𝔰𝔢(n)$ in $𝔤𝔩(n+1)$.
Note that this is *not* a transparently isometric embedding.

# Constructor

    SpecialEuclideanInGeneralLinear(n)
"""
const SpecialEuclideanInGeneralLinear =
    EmbeddedManifold{ℝ,<:SpecialEuclidean,<:GeneralLinear}

function SpecialEuclideanInGeneralLinear(n)
    return EmbeddedManifold(SpecialEuclidean(n), GeneralLinear(n + 1))
end

"""
    embed(M::SpecialEuclideanInGeneralLinear, p)

Embed the point `p` on [`SpecialEuclidean`](@ref) in the [`GeneralLinear`](@ref) group.
The embedding is calculated using [`affine_matrix`](@ref).
"""
function embed(M::SpecialEuclideanInGeneralLinear, p)
    G = M.manifold
    return affine_matrix(G, p)
end
"""
    embed(M::SpecialEuclideanInGeneralLinear, p, X)

Embed the tangent vector X at point `p` on [`SpecialEuclidean`](@ref) in the
[`GeneralLinear`](@ref) group. Point `p` can use any representation valid for
`SpecialEuclidean`. The embedding is similar from the one defined by [`screw_matrix`](@ref)
but the translation part is multiplied by inverse of the rotation part.
"""
function embed(M::SpecialEuclideanInGeneralLinear, p, X)
    G = M.manifold
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    Y = allocate_result(G, screw_matrix, nX, hX)
    nY, hY = submanifold_components(G, Y)
    copyto!(hY, hX)
    copyto!(nY, hp' * nX)
    @inbounds _padvector!(G, Y)
    return Y
end

function embed!(M::SpecialEuclideanInGeneralLinear, q, p)
    return copyto!(q, embed(M, p))
end
function embed!(M::SpecialEuclideanInGeneralLinear, Y, p, X)
    return copyto!(Y, embed(M, p, X))
end

"""
    project(M::SpecialEuclideanInGeneralLinear, p)

Project point `p` in [`GeneralLinear`](@ref) to the [`SpecialEuclidean`](@ref) group.
This is performed by extracting the rotation and translation part as in [`affine_matrix`](@ref).
"""
function project(M::SpecialEuclideanInGeneralLinear, p)
    G = M.manifold
    np, hp = submanifold_components(G, p)
    return ProductRepr(np, hp)
end
"""
    project(M::SpecialEuclideanInGeneralLinear, p, X)

Project tangent vector `X` at point `p` in [`GeneralLinear`](@ref) to the
[`SpecialEuclidean`](@ref) Lie algebra.
This reverses the transformation performed by [`embed`](@ref embed(M::SpecialEuclideanInGeneralLinear, p, X))
"""
function project(M::SpecialEuclideanInGeneralLinear, p, X)
    G = M.manifold
    np, hp = submanifold_components(G, p)
    nX, hX = submanifold_components(G, X)
    return ProductRepr(hp * nX, hX)
end

function project!(M::SpecialEuclideanInGeneralLinear, q, p)
    return copyto!(q, project(M, p))
end
function project!(M::SpecialEuclideanInGeneralLinear, Y, p, X)
    return copyto!(Y, project(M, p, X))
end


text:
function TriMesh_ElementAdjacency(eConn)
#= -------------------------------------------------------------------------78--

   Discussion:

     TriMesh_ElementAdjacency determines the adjacent elements for each
     triangular element and for each vertex in a triangular mesh.  This is based
     on the function TRIANGULARION_TRIANGLE_NEIGHBORS written by John Burkardt.

   Usage:

     eAdjacency = TriMesh_ElementAdjacency(eConn)
     # planned: eAdjacency,vConn = TriMesh_ElementAdjacency(eConn)

     where:

     eConn     is the element connectivity.  Standard orientation is assumed
               so the first three nodes in the connectivity are the vertices.
               Dimension: [ nElements, nLocalNodes ]

     eAdjacency  is a [ n_elements, 3 ] array which contains the adjacent
                 element to each triangle edge.  A value of -1 is given if that
                 edge has no adjacent element (e.g., on the boundary)

                                 3
                                | \
                         edge 2 |  \ edge 1
                                | e \
                                1----2
                                edge 3

     #vConn    is a structure of n_node arrays that list elements that are
               connected to each node.  vConn[1].e_list contains a list of the
               elements that include node 1, etc.  (planned)

   Licensing:
     This code is distributed under the MIT license.

   Modified:

     18 January 2021

   Authors:

     Jeff Borggaard and John Burkardt
     part of FEMfunctions.jl

##--------------------------------------------------------------------------78=#

  nElements, nelDOF = size(eConn)

#
#  Detect and correct 0-based indexing.
#
# minN = minimum(eConn)
# if (minN==0)
#   @printf('TriMesh_ElementAdjacency: element connectivity may be 0-based, ')
#   @printf('indexing\n will attempt to adjust and proceed\n')
#   eConn = eConn .+ 1
# end

  nNodes = maximum(eConn)

#
#  Create the triangle neighbor array.
#

#  Step 1.
#  From the list of nodes for triangle T, of the form: (I,J,K)
#  construct the three neighbor relations:
#
#    (I,J,3,T) or (J,I,3,T),
#    (J,K,1,T) or (K,J,1,T),
#    (K,I,2,T) or (I,K,2,T)
#
#  where we choose (I,J,3,T) if I < J, or else (J,I,3,T)

  faces = Array{Int64,2}(undef,4,3*nElements)
  for nElem = 1:nElements

    i = eConn[nElem,1];
    j = eConn[nElem,2];
    k = eConn[nElem,3];

    if i<j
      faces[:,1+3*(nElem-1)] = [ i, j, 3, nElem ]
    else
      faces[:,1+3*(nElem-1)] = [ j, i, 3, nElem ]
    end

    if j<k
      faces[:,2+3*(nElem-1)] = [ j, k, 1, nElem ]
    else
      faces[:,2+3*(nElem-1)] = [ k, j, 1, nElem ]
    end

    if k<i
      faces[:,3+3*(nElem-1)] = [ k, i, 2, nElem ]
    else
      faces[:,3+3*(nElem-1)] = [ i, k, 2, nElem ]
    end

  end

#  Step 2. Perform an ascending dictionary sort on the adjacency relations.
#  We only intend to sort on rows 1:2; the routine we call here
#  sorts on rows 1 through 4 but that won't hurt us.
#
#  What we need is to find cases where two triangles share a face.
#  By sorting the columns of the FACES array, we will put shared faces
#  next to each other.
#
#  this scrambles the element information.
#   faces = sort(faces,1);

#   faces = sortslices(faces,dims=1)?
   index2 = sortperm(faces[2,:])
   faces = faces[:,index2]

   index1 = sortperm(faces[1,:])
   faces = faces[:,index1]

#  Step 3. Neighboring triangles show up as consecutive columns with
#  identical first two entries.  Whenever you spot this happening,
#  make the appropriate entries in eAdjacency.

  eAdjacency = fill(-1,nElements,3)

  face = 1
  while true

    if 3nElements<=face
      break
    end

    if ( faces[1:2,face] == faces[1:2,face+1] )
      face1 = faces[3,face]
      elem1 = faces[4,face]
      face2 = faces[3,face+1]
      elem2 = faces[4,face+1]
      eAdjacency[elem1,face1] = elem2
      eAdjacency[elem2,face2] = elem1
      face = face + 2
    else # there is no neighbor, move on
      face = face + 1;
    end

  end

#  if nargout==2
#    for n=1:nNodes
#      vConn[n].eList = []
#    end
#    for n=1:nElements  % for each node, add the current element to e_list
#      for m=1:nelDOF
#        vConn[eConn[n,m]].eList = [ vConn[eConn[n,m]].eList, n ]
#      end
#    end
#  end

  return eAdjacency #,v_conn

end # function TriMesh_ElementAdjacency


text:
export Ricker
"""
    Ricker(epp::Union{ElasticPropagatorParams, AcousticPropagatorParams},
    a::Union{PyObject, <:Real},
    shift::Union{PyObject, <:Real},
    amp::Union{PyObject, <:Real}=1.0)

Returns a Ricker wavelet (a tensor).
- `epp`: a `ElasticPropagatorParams` or an `AcousticPropagatorParams`
- `a`: Width parameter
- `shift`: Center of the Ricker wavelet
- `amp`: Amplitude of the Ricker wavelet

```math
f(x) = \\mathrm{amp}A (1 - x^2/a^2) exp(-x^2/2 a^2)
```
where
```math
A = 2/sqrt(3a)pi^1/4
```
"""
function Ricker(epp,
        a,
        shift,
        amp=1.0)
    NT, T = epp.NSTEP, epp.NSTEP*epp.DELTAT
    A = @. 2 / (sqrt(3 * a) * (pi^0.25))
    wsq = @. a^2
    vec =  collect(1:NT) .-shift
    xsq = @. vec^2
    mod = @. (1 - xsq / wsq)
    gauss = @. exp(-xsq / (2 * wsq))
    total = @. amp * A * mod * gauss
    return total
end

struct SafeIndex{T}
    arg::T
end

function SafeIndex(args::Tuple)
    SafeIndex(CartesianIndex(args))
end
SafeIndex(args::Int...) = SafeIndex(args)

function Base.getindex(x::AbstractArray, si::SafeIndex)
    getindex(x, si.arg)
end

function Base.setindex!(x::AbstractArray, val, si::SafeIndex)
    setindex!(x, val, si.arg)
end


text:
@testset "customop" begin
    if isdir("temp")
        rm("temp", force=true, recursive=true)
    end
    mkdir("temp")
    cd("temp")
    customop()
    customop()
    customop(;julia=true)
    cp("$(@__DIR__)/../examples/while_loop/DirichletBD/DirichletBD.cpp", "DirichletBD.cpp", force=true)
    cp("$(@__DIR__)/../examples/while_loop/DirichletBD/DirichletBD.h", "DirichletBD.h", force=true)
    cp("$(@__DIR__)/../examples/while_loop/DirichletBD/CMakeLists.txt", "CMakeLists.txt", force=true)
    mkdir("build")
    cd("build")
    ADCME.cmake()
    ADCME.make()
    cd("..")
    rm("temp", force=true, recursive=true)
    @test true
end


text:
"""
    assigndist_betaprime(α, β; trunc_lower = -Inf, trunc_upper = Inf)

Assign parameters to a Beta prime distribution with parameters `α` and `β`, optionally
truncating the distribution.
"""
function assigndist_betaprime(α, β; trunc_lower = -Inf, trunc_upper = Inf)
    Truncated(BetaPrime(α, β), trunc_lower, trunc_upper)
end

export assigndist_betaprime


text:
export Euler

"""
    Euler()

This uses the explicit Euler method for prediction, also known as the
tangent predictor.
"""
struct Euler <: AbstractPredictor end
struct EulerCache{T} <: AbstractPredictorCache
    A::Matrix{T}
    b::Vector{T}
end

cache(::Euler, H, x, ẋ, t) = EulerCache(jacobian(H, x, t), dt(H, x, t))

function predict!(xnext, ::Euler, cache::EulerCache, H::HomotopyWithCache, x, t, Δt, ẋ, Jac)
    @inbounds for i=1:length(x)
        xnext[i] = x[i] + Δt * ẋ[i]
    end
    nothing
end

order(::Euler) = 2


text:
## tests of find_zero interface
using Roots
using Compat.Test
import SpecialFunctions.erf

meths = [Order0(), Order1(), Order2(), Order5(), Order8(), Order16()]

## basics
fns = [(x -> x^2 - exp(x) - 3x + 2,   0.257530285439860, [-1, 0.5]),
       (x -> x^3 + 4x^2 -10,          1.365230013414097, [0.5]), # not -1.0, which causes errors in most
       (x -> exp(x)*sin(x) - 2x - 5, -2.523245230732555, [-1.0, 0.0]),
       (x -> log(x) + sqrt(x) - 5,    8.309432694231571, [15, 20]),
       (x -> sqrt(x) - 1/x - 3,       9.633595562832696, [20, 25]),

       (x -> exp(-x) - cos(x), 1.292695719373398, [-1.0, .75, 1.5]),  # not 0.6 which has issues with the bump
       (x -> cos(x)^2 - x/5, 1.085982678007472, [0.5, 1.0]),
       (x -> x^10 - x^3 - x - 1, -0.674177935277052, [0.25, 0.5]),
       (x -> sin(x) - x + 1, 1.934563210752024, [2.0, 20.0]),
       (x -> exp(-x^2 + x + 2) - cos(x) + x^3 + 1, -1.0899423334391694, [-0.9, -0.6]),

       (x -> asin(x^2 -1) - x/2 + 1, 0.594810968398369, [0.3, 0.9]),
       (x -> tanh(x) - tan(x), 7.068582745628732, [5.5, 7.1])   # not 5.0, as some methods go off to vertical asymptote
]

for (i, (f, xstar, xs)) in enumerate(fns)
    for m in meths
        for x0_ in xs
            out = try
                xn = find_zero(f, x0_, m)
                @test abs(xn - xstar) < 1e-14 || abs(f(xn)) < 1e-13
                "."
            catch err
                "*"
            end
            out == "*" && println("Fn $i, method $m, x0=$x0_ failed")
        end
    end
end

## Some more extreme tests from http://ir.igsnrr.ac.cn/bitstream/311030/8840/1/%E4%BE%AF%E9%BA%9F%E7%A7%91(SCI)2.pdf
## All the methods do poorly here due to the multiplicities of the functions
## even though there should be convergnce of either |xn+1 - xn| or |f(xn)|
multiplicity_tests = [
         (x -> (x - sqrt(5))^4 / ((x-1)^2 + 2),    4.3,  2.236067977499790)
         (x -> (sin(x)^2 - 2x + 1)^5,              4.5,  0.71483582544138924)
         (x -> (8x*exp(-x^2) -2x - 3)^8,          -2.0, -1.7903531791589544)
         (x -> (2x*cos(x) + x^2 - 3)^10/(x^2 + 1), 3.0,  2.9806452794385368)
         (x -> (exp(-x^2 + x + 3) - x + 2)^9,      2.5,  2.4905398276083051)

         (x -> (exp(-x) + 2sin(x))^4,              3.0,  3.1627488709263654)
         (x -> (log(x^2 + 3x + 5) - 2x + 7)^8,     5.5,  5.4690123359101421)
         (x -> (sqrt(x^2 + 2x + 5) - 2sin(x) - x^2 + 3)^5,   2.3, 2.3319676558839640)
         (x -> (x-2)^4 / ( (x-1)^2 + 1),           1.8,  2.0000000000000000)
         (x -> abs(x - 2.5)^(15/4) * exp(x),       2.2,  2.5)

         (x -> (sqrt(x) - 1/x - 1)^7,              2.0,  2.147899035704787)
         (x -> (log(x) + sqrt(x) - 5)^3,           8.0,  8.309432694231572)
         (x -> (sin(x)*cos(x) - x^3 + 1)^9,        1.0,  1.117078770687451)
         (x -> ((x-3)*exp(x))^5,                   2.5,  3.0000000000000000)
         (x -> (log(x) + sqrt(x^4 + 1) -2)^7,      1.0,  1.222813963628973)
         ]

for (i, (fn_, x0_, xstar)) in enumerate(multiplicity_tests)
    for m in meths
        #println("$i: $m")
        @test  abs(find_zero(fn_, x0_, m, maxevals=100) - xstar) < 1e-1 # wow, not too ambitious here, 9th powers...
    end
end

## issues with starting near a maxima. Some bounce out of it, but
## one would expect all to have issues
fn, xstar = x -> x^3 + 4x^2 -10,  1.365230013414097
for m in [Order0(), Order1()]
    @test_throws Roots.ConvergenceFailed find_zero(fn, -1.0, m)
end
for m in [Order2(), Order5(), Order8(), Order16()]
    @test find_zero(fn, -1.0, m) ≈ xstar
end

## non-simple root
fn, xstar, x0 = x -> cos(x) - 1, 0.0, 0.1
for m in meths
    xn = find_zero(fn, x0, m)
    @test abs(fn(xn)) <= 1e-10
end

## issue with large steps
fn, x0 = x -> x^20 - 1, 0.5
for m in meths[2:end] # not 0, as it uses bracket
    @test_throws Roots.ConvergenceFailed find_zero(fn, x0, m)
end

## different types of initial values
for m in meths
    @test find_zero(sin, 3, m) ≈ pi
    @test find_zero(sin, 3.0, m) ≈ pi
    @test find_zero(sin, big(3), m) ≈ pi
    @test find_zero(sin, big(3.0), m) ≈ pi
end

## Methods guarded with a bracket
# fn, xstar, x0 = (x -> sin(x) - x - 1,  -1.9345632107520243, 2)
# @test_throws Roots.ConvergenceFailed find_zero(fn, x0, Order2())
# for m in meths
#     @test find_zero(fn, x0, m, bracket=[-2,3]) ≈ xstar
# end

# fn, xstar, x0 = (x -> x * exp( - x ), 0, 1.0)
# @test find_zero(fn, x0, Order0(), bracket=[-1,2]) ≈ xstar
# @test find_zero(fn, 7.0, Order0(), bracket=[-1,2]) ≈ xstar  # out of bracket

## bisection methods
@test find_zero(x -> cos(x) - x, [0, pi], Bisection()) ≈ 0.7390851332151607
@test (find_zero(x -> sin(x), [big(3), 4], Bisection()) |> sin) < 1e-70
@test find_zero(x -> sin(x), [4,3], Bisection()) ≈ pi
@test find_zero(x -> sin(x), [big(4),3], Bisection()) ≈ pi
@test find_zero(x -> Inf*sign(x-1.0), (-Inf, Inf)) ≈ 1
#@test_throws ArgumentError  find_zero(x -> sin(x), 3.0, Bisection())

@test find_zero(x -> cos(x) - x, [0, pi], FalsePosition()) ≈ 0.7390851332151607
@test (find_zero(x -> sin(x), [big(3), 4], FalsePosition(), maxevals=1000) |> sin) < 1e-70
@test find_zero(x -> sin(x), [4,3], FalsePosition()) ≈ pi

galadino_probs = [(x -> x^3 - 1, [.5, 1.5]),
                  (x -> x^2 * (x^2/3 + sqrt(2) * sin(x)) - sqrt(3)/18, [.1, 1]),
                  (x -> 11x^11 - 1, [0.1, 1]),
                  (x ->  x^3 + 1, [-1.8, 0]),
                  (x ->  x^3 - 2x - 5, [2.0, 3]),

                  ((x,n=5)  -> 2x * exp(-n) + 1 - 2exp(-n*x) , [0,1]),
                  ((x,n=10) -> 2x * exp(-n) + 1 - 2exp(-n*x) , [0,1]),
                  ((x,n=20) -> 2x * exp(-n) + 1 - 2exp(-n*x) , [0,1]),

                  ((x,n=5)  -> (1 + (1-n)^2) * x^2 - (1 - n*x)^2 , [0,1]),
                  ((x,n=10) -> (1 + (1-n)^2) * x^2 - (1 - n*x)^2 , [0,1]),
                  ((x,n=20) -> (1 + (1-n)^2) * x^2 - (1 - n*x)^2 , [0,1]),

                  ((x,n=5)  -> x^2 - (1-x)^n , [0,1]),
                  ((x,n=10) -> x^2 - (1-x)^n , [0,1]),
                  ((x,n=20) -> x^2 - (1-x)^n , [0,1]),

                  ((x,n=5)  -> (1 + (1-n)^4)*x - (1 - n*x)^4 , [0,1]),
                  ((x,n=10) -> (1 + (1-n)^4)*x - (1 - n*x)^4 , [0,1]),
                  ((x,n=20) -> (1 + (1-n)^4)*x - (1 - n*x)^4 , [0,1]),

                  ((x,n=5)  -> exp(-n*x)*(x-1) + x^n , [0,1]),
                  ((x,n=10) -> exp(-n*x)*(x-1) + x^n , [0,1]),
                  ((x,n=20) -> exp(-n*x)*(x-1) + x^n , [0,1]),

                  ((x,n=5)  -> x^2 + sin(x/n) - 1/4 , [0,1]),
                  ((x,n=10) -> x^2 + sin(x/n) - 1/4 , [0,1]),
                  ((x,n=10) -> x^2 + sin(x/n) - 1/4 , [0,1])
                  ]

for (fn_, ab) in galadino_probs
    for m in [Roots.A42(), Bisection(), (FalsePosition(i) for i in 1:12)...]
        global x0 = find_zero(fn_, ab, m, maxevals=120)
        @test abs(fn_(x0)) <= 1e-14
    end
end

## Can use (a,b) or [a,b] for bisection
fn = x -> x^5 - x - 1
for m in [Roots.A42(), Bisection(), (FalsePosition(i) for i in 1:12)...]
    global x0 = find_zero(fn, (1,2.0), m)
    @test abs(fn(x0)) <= 1e-14
end

## defaults for method argument
@test find_zero(x -> cbrt(x), 1) ≈ 0.0 # order0()
@test find_zero(sin, [3,4]) ≈ π   # Bisection()

## Callable objects
struct _SampleCallableObject end
_SampleCallableObject(x) = x^5 - x - 1
for m in meths
    @test find_zero(_SampleCallableObject, 1.0, m) ≈ 1.1673039782614187
end

### a wrapper to count function calls, say
mutable struct Cnt
    cnt::Int
    f
    Cnt(f) = new(0, f)
end
(f::Cnt)(x) = (f.cnt += 1; f.f(x))

g = Cnt(x -> x^5 - x - 1)
for m in meths
    @test find_zero(g, 1.0, m) ≈ 1.1673039782614187
end

## test tolerance arguments
fn, xstar = x -> sin(x) - x + 1, 1.9345632107520243
@test find_zero(fn, 20.0, Order2())  ≈ xstar   # needs 16 iterations, 33 fn evaluations
@test abs(fn(find_zero(fn, 20.0, Order2(), atol=1e-2)) - xstar) > 1e-12
@test abs(fn(find_zero(fn, 20.0, Order2(), rtol=1e-2)) - xstar) > 1e-12
@test_throws Roots.ConvergenceFailed find_zero(fn, 20.0, Order2(), maxevals=5)
@test_throws Roots.ConvergenceFailed find_zero(fn, 20.0, Order2(), maxfnevals=10)

## test robustness of Order0
## tests from: http://people.sc.fsu.edu/~jburkardt/cpp_src/test_zero/test_zero.html
function _newton_baffler(x)
    if ( x - 0.0 ) < -0.25
        0.75 * ( x - 0 ) - 0.3125
    elseif  ( x - 0 ) < 0.25
        2.0 * ( x - 0 )
    else
        0.75 * ( x - 0 ) + 0.3125
    end
end
pathological = [
                (x ->sin( x ) - x / 2, .1),
                (x-> 2 * x - exp( - x ), 1),
                (x -> x * exp( - x ), 0.99),
                (x -> exp( x ) - 1 / ( 10 * x )^2, .1),
                (x -> ( x + 3 ) * ( x - 1 )^2, 1),

                (x -> exp( x ) - 2 - 1 / ( 10 * x )^2 + 2 / ( 100 * x )^3, 1),
                (x -> cbrt(x), 1),
                (x -> cos( x ) - x, 1),
                (x-> _newton_baffler(x), 8),
                (x -> 20.0 * x / ( 100.0 * x^2 + 1.0), 0.095),

                (x ->  ( 4.0 + x^2) * ( 2.0 + x ) * ( 2.0 - x )  / ( 16.0 * x * x * x * x + 0.00001 ), 1),
                (x -> (x == 1.0) ? float(0) : sign(x-1.0) * exp(log(1e4) + log(abs(x - 1.0)) - 1.0/(x-1.0)^2), 1),
                (x -> 0.00000000001 * (x - 100.0), 1),
                (x -> 1.0 / ( ( x - 0.3 ) * ( x - 0.3 ) + 0.01 ) + 1.0 / ( ( x - 0.9 ) * ( x - 0.9 ) + 0.04 ) + 2.0 * x - 5.2, -1),
                (x -> ( 1 - 6x^2) * cbrt(x) * exp(-x^2) / (3*x), -0.25),

                (x -> ( pi * ( x - 5.0 ) / 180.0 ) - 0.8 * sin( pi * x / 180.0 ), 1),
                (x -> x^3 - 2*x - 5, 2),
(x -> 1e6 * (x^7 -7x^6 +21x^5 -35x^4 +35x^3-21x^2+7x-1),  0.990),
(x -> cos(100*x)-4*erf(30*x-10), 0.0)
                ]

for (fn_, x0_) in pathological
    find_zero(fn_, x0_, Order0())
end

## issue tests: put in tests to ensure closed issues don't reappear.

## issue #94; tolerances not matching documentation
test_94 = function(;kwargs...)
    g, T = 1.62850, 14.60000
    α, t1, tf = 0.00347, 40.91375, 131.86573
    y, ya, yf = 0.0, 9000.0, 10000.0
    vy = sqrt(2g*(ya-y))
    θ0, θ1 = atan(α*tf), atan(α*(tf-t1))
    I_sintan(x) = tan(x)/2cos(x) - atanh(tan(x/2))
    I_sintan(x, y) = I_sintan(y) - I_sintan(x)
    function lhs(θ)
        tRem = (vy - T/α*(sec(θ1) - sec(θ))) / g
        val = -yf + y + vy*tRem - 0.5g*tRem^2 - T/α^2*I_sintan(θ, θ1)
        val
    end

    meth = Roots.FalsePosition()
    f, x0 = lhs, [atan(α*tf), atan(α*(tf-t1))]
    F = Roots.DerivativeFree(lhs)
    state = Roots.init_state(meth, F, x0)
    options = Roots.init_options(meth, state)
    find_zero(meth, F, options, state)

    @test state.steps <= 15
end
test_94()

@testset "init_state!" begin
    ## add init_state! method
    g1(x) = x^5 - x - 1
    x0_, xstar_ = (1.0, 2.0), 1.1673039782614187
    M = Roots.A42()
    state = Roots.init_state(M, g1, x0_)
    options = Roots.init_options(M, state)
    for M in (Roots.A42(), Roots.Bisection(), Roots.FalsePosition())
        Roots.init_state!(state, M, g1, x0_)
        @test find_zero(M, g1, options, state) ≈ xstar_
    end

    ## hybrid
    g1(x) = exp(x) - x^4
    x0_, xstar_ = (5.0, 20.0), 8.613169456441398
    M = Roots.Bisection()
    state = Roots.init_state(M, g1, x0_)
    options = Roots.init_options(M, state, xatol=1/2)
    find_zero(M, g1, options, state)
    M = Roots.Order1()
    Roots.init_state!(state, M, g1, state.m[1])
    Roots.init_options!(options, M)
    @test find_zero(M, g1, options, state) ≈ xstar_

    # test conversion between states/options
    Ms = vcat([Roots.FalsePosition(i) for i in 1:12], Roots.A42(), Roots.AlefeldPotraShi(), Roots.Brent(), Roots.Bisection())
    Ns = [Roots.Order1(), Roots.Order2(), Roots.Order5(), Roots.Order8(), Roots.Order16()]

    g1(x) = exp(x) - x^4
    x0_, xstar_ = (5.0, 20.0), 8.613169456441398
    M = Ms[1]
    state = Roots.init_state(M, g1, x0_)
    options = Roots.init_options(M, state)

    for M in Ms
        @test Roots.init_state!(state, M, g1, x0_) == nothing
        @test Roots.init_options!(options, M) == nothing
        @test Roots.update_state(M, g1, state, options) == nothing
        @test Roots.assess_convergence(M, state, options) == false # none in just 1 step
    end

    x0_ = 8.0
    for M in Ns
        @test Roots.init_state!(state, M, g1, x0_) == nothing
        @test Roots.init_options!(options, M) == nothing
        @test Roots.update_state(M, g1, state, options)  == nothing
        @test Roots.assess_convergence(M, state, options) == false  # none in 1 step
    end

end


text:
using Lasso, GLM, Distributions, GLMNet, FactCheck

datapath = joinpath(dirname(@__FILE__), "data")

testpath(T::DataType, d::Normal, l::GLM.Link, nsamples::Int, nfeatures::Int) =
    joinpath(datapath, "$(T)_$(typeof(d).name.name)_$(typeof(l).name.name)_$(nsamples)_$(nfeatures).tsv")

function makeX(ρ, nsamples, nfeatures, sparse)
    Σ = fill(ρ, nfeatures, nfeatures)
    Σ[diagind(Σ)] = 1
    X = rand(MvNormal(Σ), nsamples)'
    sparse && (X[randperm(length(X))[1:round(Int, length(X)*0.95)]] = 0)
    β = [(-1)^j*exp(-2*(j-1)/20) for j = 1:nfeatures]
    (X, β)
end

randdist(::Normal, x) = rand(Normal(x))
randdist(::Binomial, x) = rand(Bernoulli(x))
randdist(::Poisson, x) = rand(Poisson(x))
function genrand(T::DataType, d::Distribution, l::GLM.Link, nsamples::Int, nfeatures::Int, sparse::Bool)
    X, coef = makeX(0.0, nsamples, nfeatures, sparse)
    y = X*coef
    for i = 1:length(y)
        y[i] = randdist(d, linkinv(l, y[i]))
    end
    (X, y)
end

# Test against GLMNet
facts("LassoPath") do
    for (dist, link) in ((Normal(), IdentityLink()), (Binomial(), LogitLink()), (Poisson(), LogLink()))
        context("$(typeof(dist).name.name) $(typeof(link).name.name)") do
            for sp in (false, true)
                srand(1337)
                context(sp ? "sparse" : "dense") do
                    (X, y) = genrand(Float64, dist, link, 1000, 10, sp)
                    yoff = randn(length(y))
                    for intercept = (false, true)
                        context("$(intercept ? "w/" : "w/o") intercept") do
                            for alpha = [1, 0.5]
                                context("alpha = $alpha") do
                                    for offset = Vector{Float64}[Float64[], yoff]
                                        context("$(isempty(offset) ? "w/o" : "w/") offset") do
                                            # First fit with GLMNet
                                            if isa(dist, Normal)
                                                yp = isempty(offset) ? y : y + offset
                                                ypstd = std(yp, corrected=false)
                                                # glmnet does this on entry, which changes λ mappings, but not
                                                # coefficients. Should we?
                                                yp ./= ypstd
                                                !isempty(offset) && (offset ./= ypstd)
                                                y ./= ypstd
                                                g = glmnet(X, yp, dist, intercept=intercept, alpha=alpha, tol=10*eps())
                                            elseif isa(dist, Binomial)
                                                yp = zeros(size(y, 1), 2)
                                                yp[:, 1] = y .== 0
                                                yp[:, 2] = y .== 1
                                                g = glmnet(X, yp, dist, intercept=intercept, alpha=alpha, tol=10*eps(),
                                                           offsets=isempty(offset) ? zeros(length(y)) : offset)
                                            else
                                                g = glmnet(X, y, dist, intercept=intercept, alpha=alpha, tol=10*eps(),
                                                           offsets=isempty(offset) ? zeros(length(y)) : offset)
                                            end
                                            gbeta = convert(Matrix{Float64}, g.betas)

                                            for randomize = VERSION >= v"0.4-dev+1915" ? [false, true] : [false]
                                                context(randomize ? "random" : "sequential") do
                                                    niter = 0
                                                    for naivealgorithm in (false, true)
                                                         context(naivealgorithm ? "naive" : "covariance") do
                                                            for spfit in (false, true)
                                                                context(spfit ? "as SparseMatrixCSC" : "as Matrix") do
                                                                    # Now fit with Lasso
                                                                    l = fit(LassoPath, spfit ? sparse(X) : X, y, dist, link,
                                                                            λ=g.lambda, naivealgorithm=naivealgorithm, intercept=intercept,
                                                                            cd_tol=10*eps(), irls_tol=10*eps(), criterion=:coef, randomize=randomize,
                                                                            α=alpha, offset=offset)

                                                                    # rd = (l.coefs - gbeta)./gbeta
                                                                    # rd[!isfinite(rd)] = 0
                                                                    # println("         coefs adiff = $(maxabs(l.coefs - gbeta)) rdiff = $(maxabs(rd))")
                                                                    # rd = (l.b0 - g.a0)./g.a0
                                                                    # rd[!isfinite(rd)] = 0
                                                                    # println("         b0    adiff = $(maxabs(l.b0 - g.a0)) rdiff = $(maxabs(rd))")
                                                                    @fact l.coefs --> roughly(gbeta, 5e-7)
                                                                    @fact l.b0 --> roughly(g.a0, 2e-5)

                                                                    # Ensure same number of iterations with all algorithms
                                                                    if niter == 0
                                                                        niter = l.niter
                                                                    else
                                                                        @fact l.niter --> niter
                                                                    end
                                                                end
                                                            end
                                                        end
                                                    end
                                                end
                                            end
                                        end
                                    end
                                end
                            end
                        end
                    end
                end
            end
        end
    end
end

# Test for sparse matrices


text:
# Autogenerated wrapper script for Libtask_jll for x86_64-w64-mingw32
export libtask_julia

JLLWrappers.@generate_wrapper_header("Libtask")
JLLWrappers.@declare_library_product(libtask_julia, "libtask_julia.dll")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libtask_julia,
        "bin\\libtask_julia.dll",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
function Base.:(==)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:CompoundPeriod}
    return Base.:(==)(canonical(x), canonical(y))
end
function Base.:(!=)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:CompoundPeriod}
    return Base.:(!=)(canonical(x), canonical(y))
end

function Base.:(<)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:CompoundPeriod}
    cx, cy = canonical(x), canonical(y)
    xmonths, ymonths = Months(cx), Months(cy)
    xsecs, ysecs = Seconds(cx), Seconds(cy)
    xnanos, ynanos = Nanoseconds(x-xsecs), Nanoseconds(y-ysecs)
    return xmonths < ymonths || (xmonths == ymonths && (xsecs < ysecs || (xsecs == ysecs && xnanos < ynanos)))
end

function Base.:(>)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:CompoundPeriod}
    cx, cy = canonical(x), canonical(y)
    xmonths, ymonths = Months(cx), Months(cy)
    xsecs, ysecs = Seconds(cx), Seconds(cy)
    xnanos, ynanos = Nanoseconds(x-xsecs), Nanoseconds(y-ysecs)
    return xmonths > ymonths || (xmonths == ymonths && (xsecs > ysecs || (xsecs == ysecs && xnanos > ynanos)))
end

function Base.:(<=)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:CompoundPeriod}
    cx, cy = canonical(x), canonical(y)
    xmonths, ymonths = Months(cx), Months(cy)
    xsecs, ysecs = Seconds(cx), Seconds(cy)
    xnanos, ynanos = Nanoseconds(x-xsecs), Nanoseconds(y-ysecs)
    return xmonths < ymonths || (xmonths == ymonths && (xsecs < ysecs || (xsecs == ysecs && xnanos <= ynanos)))
end

function Base.:(>=)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:CompoundPeriod}
    cx, cy = canonical(x), canonical(y)
    xmonths, ymonths = Months(cx), Months(cy)
    xsecs, ysecs = Seconds(cx), Seconds(cy)
    xnanos, ynanos = Nanoseconds(x-xsecs), Nanoseconds(y-ysecs)
    return xmonths > ymonths || (xmonths == ymonths && (xsecs > ysecs || (xsecs == ysecs && xnanos >= ynanos)))
end

function Base.:(==)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:Period}
  mx, my = Months(x), Months(y)
  mx !== my && return false

  mintyp = Dates.periodisless(mintype(x)(1), P2(1)) ? mintype(x) : P2
  typ = plural(mintyp)
  xx = typ(x)
  yy = typ(y)
  xx.value == yy.value
end

function Base.:(!=)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:Period}
  mx, my = Months(x), Months(y)
  mx === my && return false

  mintyp = Dates.periodisless(mintype(x)(1), P2(1)) ? mintype(x) : P2
  typ = plural(mintyp)
  xx = typ(x)
  yy = typ(y)
  xx.value != yy.value
end

function Base.:(<)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:Period}
  mx, my = Months(x), Months(y)
  mx < my && return true
  mx !== my && return false

  mintyp = Dates.periodisless(mintype(x)(1), P2(1)) ? mintype(x) : P2
  typ = plural(mintyp)
  xx = typ(x)
  yy = typ(y)
  xx.value < yy.value
end

function Base.:(<=)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:Period}
  mx, my = Months(x), Months(y)
  mx < my && return true
  mx > my && return false

  mintyp = Dates.periodisless(mintype(x)(1), P2(1)) ? mintype(x) : P2
  typ = plural(mintyp)
  xx = typ(x)
  yy = typ(y)
  xx.value <= yy.value
end

function Base.:(>)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:Period}
  mx, my = Months(x), Months(y)
  mx > my && return true
  mx !== my && return false

  mintyp = Dates.periodisless(mintype(x)(1), P2(1)) ? mintype(x) : P2
  typ = plural(mintyp)
  xx = typ(x)
  yy = typ(y)
  xx.value > yy.value
end

function Base.:(>=)(x::P1, y::P2) where {P1<:CompoundPeriod, P2<:Period}
  mx, my = Months(x), Months(y)
  mx > my && return true
  mx < my && return false

  mintyp = Dates.periodisless(mintype(x)(1), P2(1)) ? mintype(x) : P2
  typ = plural(mintyp)
  xx = typ(x)
  yy = typ(y)
  xx.value >= yy.value
end

const PeriodTypes = (:Year, :Month, :Week, :Day, :Hour, :Minute, :Second, :Millisecond, :Microsecond, :Nanosecond)

for P in PeriodTypes
    @eval Base.isless(::Type{$P}, ::Type{$P}) = false
    @eval Base.isequal(::Type{$P}, ::Type{$P}) = true
end
for idx1 in 1:(length(PeriodTypes)-1)
    P1 = PeriodTypes[idx1]
    for idx2 in (idx1+1):length(PeriodTypes)
        P2 = PeriodTypes[idx2]
        @eval Base.isless(::Type{$P2}, ::Type{$P1}) = true
        @eval Base.isless(::Type{$P1}, ::Type{$P2}) = false
        @eval Base.isequal(::Type{$P2}, ::Type{$P1}) = false
        @eval Base.isequal(::Type{$P1}, ::Type{$P2}) = false
    end
end


text:
module CubedSphereDiscreteModelsTestsSeq
  using Gridap
  using GridapGeosciences
  using FillArrays
  using Test

  include("../CubedSphereDiscreteModelsTests.jl")
  include("../ConvergenceAnalysisTools.jl")

  n_values=generate_n_values(2)
  hs=[2.0/n for n in n_values]

  model_args_series=zip(n_values,Fill(1,length(n_values)))
  _,_,s=convergence_discrete_cubed_sphere_surface(hs,model_args_series,0)
  @test round(s,digits=1) ≈ 2.0

  model_args_series=zip(n_values,Fill(2,length(n_values)))
  _,_,s=convergence_discrete_cubed_sphere_surface(hs,model_args_series,4)
  @test round(s,digits=1) ≈ 4.0

  model_args_series=[(n,) for n in n_values]
  _,_,s=convergence_discrete_cubed_sphere_surface(hs,model_args_series,4)
  @test round(s,digits=1) ≈ 5.7

end

max_stars_repo_name: santiagobadia/GridapGeosciences

max_issues_repo_name: santiagobadia/GridapGeosciences

max_forks_repo_name: santiagobadia/GridapGeosciences