text:
include("utils.jl")

filename = "data/earth_snapshots.jld2"
# filename = "data/small_earth_snapshots.jld2"
jl_file = jldopen(filename, "r+")

lat_grid = jl_file["grid"]["latitude"]
lon_grid = jl_file["grid"]["longitude"]
rad_grid = jl_file["grid"]["radius"]
λ = lon_grid
ϕ = lat_grid
r = rad_grid

fig = Figure(resolution = (1700, 1000))
s_string = "T" # state string
state_file =  jl_file[s_string]["1"][:,:,1,1]
clims = quantile.(Ref(state_file[:]), (0.0,1.0))
colormap =  :afmhot # :thermal # :binary
titlestring =  "Temperature"
ax1 = fig[1,1] = Axis(fig, title = titlestring, titlesize = 30)
# fig[3,2+3*2] = Label(fig, "ρv: latlon ", textsize = 30)
hm = heatmap!(ax1, λ, ϕ, state_file, colormap = colormap, colorrange = clims, interpolate = true, shading = false, show_axis=false)
cbar1 = Colorbar(fig, hm, label = L" $T$ [K]",width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize=25, tickalign = 1, height = Relative(3/4)
)
fig[1,2] = cbar1

s_string = "p" # state string
state_file =  jl_file[s_string]["1"][:,:,1,1] ./ 1000
clims = quantile.(Ref(state_file[:]), (0.0,1.0))
colormap = :bone_1 # :acton # :bamako # :balance # :thermal
titlestring =  "Pressure"
ax2 = fig[2,1] = Axis(fig, title = titlestring, titlesize = 30)
# fig[3,2+3*2] = Label(fig, "ρv: latlon ", textsize = 30)
hm2 = heatmap!(ax2, λ, ϕ, state_file, colormap = colormap, colorrange = clims, interpolate = true, shading = false, show_axis=false)
cbar2 = Colorbar(fig, hm2, label = L" $P$ [kpa]",width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize=25, tickalign = 1, height = Relative(3/4)
)
fig[2,2] = cbar2

s_string = "ρu" # state string
state_file =  jl_file[s_string]["1"][:,:,1,1] ./ jl_file["ρ"]["1"][:,:,1,1]
clims = quantile.(Ref(state_file[:]), (0.0,1.0))
clims = (-clims[2], clims[2])
colormap = :balance # :thermal
titlestring =  "Zonal Velocity"
ax3 = fig[1,3] = Axis(fig, title = titlestring, titlesize = 30)
# fig[3,2+3*2] = Label(fig, "ρv: latlon ", textsize = 30)
hm3 = heatmap!(ax3, λ, ϕ, state_file, colormap = colormap, colorrange = clims, interpolate = true, shading = false, show_axis=false)
cbar3 = Colorbar(fig, hm3, label = L" $u$ [m/s]",width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize=25, tickalign = 1, height = Relative(3/4)
)
fig[1,4] = cbar3

s_string = "ρv" # state string
state_file =  jl_file[s_string]["1"][:,:,1,1] ./ jl_file["ρ"]["1"][:,:,1,1]
clims = quantile.(Ref(state_file[:]), (0.0,1.0))
clims = (-clims[2], clims[2])
colormap = :balance # :thermal
titlestring =  "Meriodonal Velocity"
ax4 = fig[2,3] = Axis(fig, title = titlestring, titlesize = 30)
# fig[3,2+3*2] = Label(fig, "ρv: latlon ", textsize = 30)
hm4 = heatmap!(ax4, λ, ϕ, state_file, colormap = colormap, colorrange = clims, interpolate = true, shading = false, show_axis=false)
cbar4 = Colorbar(fig, hm4, label = L" $v$ [m/s]",width = 25, ticklabelsize = 30,
    labelsize = 30, ticksize=25, tickalign = 1, height = Relative(3/4)
)
fig[2,4] = cbar4

for ax in [ax1, ax2, ax3, ax4]
    ax.limits = (extrema(λ)..., extrema(ϕ)...)

    ax.xlabel = "Longitude [ᵒ]"
    ax.ylabel = "Latitude [ᵒ]"
    ax.xlabelsize = 25
    ax.ylabelsize = 25
    ax.yticks = ([-80, -60, -30, 0, 30, 60, 80], ["80S", "60S", "30S", "0", "30N", "60N", "80N"])
    ax.xticks = ([-160, -120, -80, -40, 0, 40, 80, 120, 160], ["160W", "120W", "80W", "40W", "0", "40E", "80E", "120E", "160E"])
end
update!(fig.scene)
display(fig)

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
@testset "Special Functions" begin
    vexpq = :(for i ∈ eachindex(a)
              b[i] = exp(a[i])
              end)
    lsvexp = LoopVectorization.loopset(vexpq);
    @test LoopVectorization.choose_order(lsvexp) == (Symbol[:i], :i, Symbol("##undefined##"), :i, 2, -1)

    function myvexp!(b, a)
        @inbounds for i ∈ eachindex(a)
            b[i] = exp(a[i])
        end
    end
    function myvexpavx!(b, a)
        @avx for i ∈ eachindex(a)
            b[i] = Base.exp(a[i])
        end
    end
    function myvexp_avx!(b, a)
        @_avx for i ∈ eachindex(a)
            b[i] = exp(a[i])
        end
    end
    function offset_exp!(A, B)
        @avx for i=1:size(A,1), j=1:size(B,2)
            A[i,j] = exp(B[i,j])
        end
    end
    function offset_expavx!(A, B)
        @avx for i=1:size(A,1), j=1:size(B,2)
            A[i,j] = exp(B[i,j])
        end
    end
    function offset_exp_avx!(A, B)
        @_avx for i=1:size(A,1), j=1:size(B,2)
            A[i,j] = exp(B[i,j])
        end
    end

    vexpsq = :(for i ∈ eachindex(a)
               s += exp(a[i])
               end)
    lsvexps = LoopVectorization.loopset(vexpsq);
    @test LoopVectorization.choose_order(lsvexps) == (Symbol[:i], :i, Symbol("##undefined##"), :i, 2, -1)

    function myvexp(a)
        s = zero(eltype(a))
        @inbounds for i ∈ eachindex(a)
            s += exp(a[i])
        end
        s
    end
    function myvexpavx(a)
        s = zero(eltype(a))
        @avx for i ∈ eachindex(a)
            s += exp(a[i])
        end
        s
    end

    function myvexp_avx(a)
        s = zero(eltype(a))
        @_avx for i ∈ eachindex(a)
            s += exp(a[i])
        end
        s
    end
    function trianglelogdetavx(L)
        ld = zero(eltype(L))
        @avx for i ∈ 1:size(L,1)
            ld += log(L[i,i])
        end
        ld
    end
    function trianglelogdet_avx(L)
        ld = zero(eltype(L))
        @_avx for i ∈ 1:size(L,1)
            ld += log(L[i,i])
        end
        ld
    end
    function testrepindshigherdim_avx(L)
        ld = zero(eltype(L))
        @_avx for i ∈ axes(L,1), j ∈ axes(L,2)
            ld += log(L[i,j,i])
        end
        ld
    end
    function testrepindshigherdimavx(L)
        ld = zero(eltype(L))
        @avx for i ∈ axes(L,1), j ∈ axes(L,2)
            ld += log(L[i,j,i])
        end
        ld
    end
    ldq = :(for i ∈ 1:size(L,1)
            ld += log(L[i,i])
            end)
    lsld = LoopVectorization.loopset(ldq);
    @test LoopVectorization.choose_order(lsld) == (Symbol[:i], :i, Symbol("##undefined##"), :i, 2, -1)

    function calc_sins!(res::AbstractArray{T}) where {T}
        code_phase_delta = T(0.01)
        @inbounds for i ∈ eachindex(res)
            res[i] = sin(i * code_phase_delta)
        end
    end
    function calc_sinsavx!(res::AbstractArray{T}) where {T}
        code_phase_delta = T(0.01)
        @avx for i ∈ eachindex(res)
            res[i] = sin(i * code_phase_delta)
        end
    end
    function calc_sins_avx!(res::AbstractArray{T}) where {T}
        code_phase_delta = T(0.01)
        @_avx for i ∈ eachindex(res)
            res[i] = sin(i * code_phase_delta)
        end
    end

    function logsumexp!(r::AbstractArray{T}, x::AbstractArray{T}) where {T}
        n = length(x)
        length(r) == n || throw(DimensionMismatch())
        isempty(x) && return -T(Inf)
        1 == LinearAlgebra.stride1(r) == LinearAlgebra.stride1(x) || throw(error("Arrays not strided"))

        u = maximum(x)                                       # max value used to re-center
        abs(u) == Inf && return any(isnan, x) ? T(NaN) : u   # check for non-finite values
        s = zero(T)
        @inbounds for i = 1:n
            tmp = exp(x[i] - u)
            r[i] = tmp
            s += tmp
        end
        invs = inv(s)
        r .*= invs

        return log1p(s-1) + u
    end
    function logsumexpavx!(r::AbstractArray{T}, x::AbstractArray{T}) where {T}
        n = length(x)
        length(r) == n || throw(DimensionMismatch())
        isempty(x) && return -T(Inf)
        1 == LinearAlgebra.stride1(r) == LinearAlgebra.stride1(x) || throw(error("Arrays not strided"))

        u = maximum(x)                                       # max value used to re-center
        abs(u) == Inf && return any(isnan, x) ? T(NaN) : u   # check for non-finite values
        s = zero(T)
        @avx for i = 1:n
            tmp = exp(x[i] - u)
            r[i] = tmp
            s += tmp
        end

        invs = inv(s)
        r .*= invs

        return log1p(s-1) + u
    end
    function logsumexp_avx!(r::AbstractArray{T}, x::AbstractArray{T}) where {T}
        n = length(x)
        length(r) == n || throw(DimensionMismatch())
        isempty(x) && return -T(Inf)
        1 == LinearAlgebra.stride1(r) == LinearAlgebra.stride1(x) || throw(error("Arrays not strided"))

        u = maximum(x)                                       # max value used to re-center
        abs(u) == Inf && return any(isnan, x) ? T(NaN) : u   # check for non-finite values
        s = zero(T)
        @_avx for i = 1:n
            tmp = exp(x[i] - u)
            r[i] = tmp
            s += tmp
        end

        invs = inv(s)
        r .*= invs

        return log1p(s-1) + u
    end
    # feq = :(for i = 1:n
    #         tmp = exp(x[i] - u)
    #         r[i] = tmp
    #         s += tmp
    #         end)
    # lsfeq = LoopVectorization.loopset(feq);
    # # lsfeq.operations

    function vpow0!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 0
        end; y
    end
    function vpown1!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ -1
        end; y
    end
    function vpow1!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 1
        end; y
    end
    function vpown2!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ -2
        end; y
    end
    function vpow2!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 2
        end; y
    end
    function vpown3!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ -3
        end; y
    end
    function vpow3!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 3
        end; y
    end
    function vpown4!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ -4
        end; y
    end
    function vpow4!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 4
        end; y
    end
    function vpown5!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ -5
        end; y
    end
    q = :(for i ∈ eachindex(y, x)
            y[i] = x[i] ^ -5
          end);
    ls = LoopVectorization.loopset(q);
    q2 = :(for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 5
          end);
    ls2 = LoopVectorization.loopset(q2);

    function vpow5!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 5
        end; y
    end
    function vpowf!(y, x)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ 2.3
        end; y
    end
    function vpowf!(y, x, p::Number)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ p
        end; y
    end
    function vpowf!(y, x, p::AbstractArray)
        @avx for i ∈ eachindex(y, x)
            y[i] = x[i] ^ p[i]
        end; y
    end

    function csetanh!(y, z, x)
        for j in axes(x, 2)
            for i = axes(x, 1)
                t2 = inv(tanh(x[i, j]))
                t1 = tanh(x[i, j])
                y[i, j] = z[i, j] * (-(1 - t1 ^ 2) * t2)
            end
        end
        y
    end
    function csetanhavx!(y, z, x)
        @avx for j in axes(x, 2)
            for i = axes(x, 1)
                t2 = inv(tanh(x[i, j]))
                t1 = tanh(x[i, j])
                y[i, j] = z[i, j] * (-(1 - t1 ^ 2) * t2)
            end
        end
        y
    end

    function transposedvectoraccess(x::AbstractVector{T}) where T
        N = length(x)
        ent = zeros(T, N)
        x isa AbstractVector && (x = x')
        for i = 1:N, j = 1:size(x,1)
            ent[i] += x[j,i]*log(x[j,i])
        end
        ent
    end
    function transposedvectoraccessavx(x::AbstractVector{T}) where T
        N = length(x)
        ent = zeros(T, N)
        x isa AbstractVector && (x = x')
        @avx for i = 1:N, j = 1:size(x,1)
            ent[i] += x[j,i]*log(x[j,i])
        end
        ent
    end

    function vsincos(x)
        y = Matrix{eltype(x)}(undef, length(x), 2)
        for i = eachindex(x)
            y[i,1], y[i,2] = sincos(x[i])
        end
        return y
    end
    function vsincosavx(x)
        y = Matrix{eltype(x)}(undef, length(x), 2)
        @avx for i = eachindex(x)
            y[i,1], y[i,2] = sincos(x[i])
        end
        return y
    end
    function sincosdot(x)
        a = zero(eltype(x))
        for i ∈ eachindex(x)
            s, c = sincos(x[i])
            a += s * c
        end
        a
    end
    function sincosdotavx(x)
        a = zero(eltype(x))
        @avx for i ∈ eachindex(x)
            s, c = sincos(x[i])
            a += s * c
        end
        a
    end

    for T ∈ (Float32, Float64)
        @show T, @__LINE__
        a = randn(T, 127);
        b1 = similar(a);
        b2 = similar(a);

        myvexp!(b1, a)
        myvexpavx!(b2, a)
        @test b1 ≈ b2
        fill!(b2, -999.9); myvexp_avx!(b2, a)
        @test b1 ≈ b2
        s = myvexp(a)
        @test s ≈ myvexpavx(a)
        @test s ≈ myvexp_avx(a)
        @test b1 ≈ @avx exp.(a)

        @test vsincos(a) ≈ vsincosavx(a)
        @test sincosdot(a) ≈ sincosdotavx(a)

        A = rand(T, 73, 73);
        ld = logdet(UpperTriangular(A))
        @test ld ≈ trianglelogdetavx(A)
        @test ld ≈ trianglelogdet_avx(A)
        Adim3 = rand(T, 37, 13, 37);
        ld = sum(i -> logdet(UpperTriangular(@view(Adim3[:,i,:]))), axes(Adim3,2))
        @test ld ≈ testrepindshigherdimavx(Adim3)
        @test ld ≈ testrepindshigherdim_avx(Adim3)

        x = rand(T, 999);
        r1 = similar(x);
        r2 = similar(x);
        lse = logsumexp!(r1, x);
        @test logsumexpavx!(r2, x) ≈ lse
        @test r1 ≈ r2
        fill!(r2, T(NaN));
        @test logsumexp_avx!(r2, x) ≈ lse
        @test r1 ≈ r2

        calc_sins!(r1)
        calc_sinsavx!(r2)
        @test r1 ≈ r2
        fill!(r2, NaN); calc_sins_avx!(r2)
        @test r1 ≈ r2

        N,M = 47,53
        B = reshape(cumsum(ones(T, 3N)),N,:)
        A1 = zeros(T, N, M)
        A2 = zeros(T, N, M)
        offset_exp!(A1, B)
        offset_expavx!(A2, B)
        @test A1 ≈ A2
        fill!(A2, 0); offset_exp_avx!(A2, B)
        @test A1 ≈ A2

        @test all(isone, vpow0!(r1, x))
        @test vpown1!(r1, x) ≈ map!(inv, r2, x)
        @test vpow1!(r1, x) == x
        @test vpown2!(r1, x) ≈ map!(abs2 ∘ inv, r2, x)
        @test vpow2!(r1, x) ≈ map!(abs2, r2, x)
        @test vpown3!(r1, x) ≈ (r2 .= x .^ -3)
        @test vpow3!(r1, x) ≈ (r2 .= x .^ 3)
        @test vpown4!(r1, x) ≈ (r2 .= x .^ -4)
        @test vpow4!(r1, x) ≈ (r2 .= x .^ 4)
        @test vpown5!(r1, x) ≈ (r2 .= x .^ -5)
        @test vpow5!(r1, x) ≈ (r2 .= x .^ 5)
        @test vpowf!(r1, x) ≈ (r2 .= x .^ 2.3)
        @test vpowf!(r1, x, -1.7) ≈ (r2 .= x .^ -1.7)
        p = randn(length(x));
        @test vpowf!(r1, x, x) ≈ (r2 .= x .^ x)

        X = rand(T, N, M); Z = rand(T, N, M);
        Y1 = similar(X); Y2 = similar(Y1);
        @test csetanh!(Y1, X, Z) ≈ csetanhavx!(Y2, X, Z)

        x = rand(T, 97);
        @test transposedvectoraccessavx(x) ≈ transposedvectoraccess(x)
    end
end

ext: jl


text:
const debug_callback_c = Ref{Ptr{Cvoid}}(C_NULL)

function init(;
    instance_layers = [],
    instance_extensions = [],
    device_extensions = [],
    enabled_features = PhysicalDeviceFeatures(),
    nqueues = 1,
    queue_flags = QUEUE_COMPUTE_BIT | QUEUE_GRAPHICS_BIT,
    with_validation = true,
)

    if with_validation && "VK_LAYER_KHRONOS_validation" ∉ instance_layers
        push!(instance_layers, "VK_LAYER_KHRONOS_validation")
    end
    if "VK_EXT_debug_utils" ∉ instance_extensions
        push!(instance_extensions, "VK_EXT_debug_utils")
    end

    available_layers = unwrap(enumerate_instance_layer_properties())
    unsupported_layers = filter(!in(getproperty.(available_layers, :layer_name)), instance_layers)
    if !isempty(unsupported_layers)
        error("Requesting unsupported instance layers: $unsupported_layers")
    end

    available_extensions = unwrap(enumerate_instance_extension_properties())
    unsupported_extensions = filter(!in(getproperty.(available_extensions, :extension_name)), instance_extensions)
    if !isempty(unsupported_extensions)
        error("Requesting unsupported instance extensions: $unsupported_extensions")
    end

    instance = Instance(instance_layers, instance_extensions; application_info = ApplicationInfo(v"1", v"1", v"1.2"))
    messenger = DebugUtilsMessengerEXT(
        instance,
        |(
            DEBUG_UTILS_MESSAGE_SEVERITY_VERBOSE_BIT_EXT,
            DEBUG_UTILS_MESSAGE_SEVERITY_INFO_BIT_EXT,
            DEBUG_UTILS_MESSAGE_SEVERITY_WARNING_BIT_EXT,
            DEBUG_UTILS_MESSAGE_SEVERITY_ERROR_BIT_EXT,
        ),
        |(
            DEBUG_UTILS_MESSAGE_TYPE_GENERAL_BIT_EXT,
            DEBUG_UTILS_MESSAGE_TYPE_PERFORMANCE_BIT_EXT,
            DEBUG_UTILS_MESSAGE_TYPE_VALIDATION_BIT_EXT,
        ),
        debug_callback_c[],
        function_pointer(instance, "vkCreateDebugUtilsMessengerEXT"),
        function_pointer(instance, "vkDestroyDebugUtilsMessengerEXT"),
    )
    physical_device = first(unwrap(enumerate_physical_devices(instance)))

    # TODO: check for supported device features
    available_extensions = unwrap(enumerate_device_extension_properties(physical_device))
    unsupported_extensions = filter(!in(getproperty.(available_extensions, :extension_name)), device_extensions)
    if !isempty(unsupported_extensions)
        error("Requesting unsupported device extensions: $unsupported_extensions")
    end

    queue_index = find_queue_family(physical_device, queue_flags)
    device = Device(
        physical_device,
        [DeviceQueueCreateInfo(queue_index, ones(Float32, nqueues))],
        [],
        device_extensions;
        enabled_features,
    )
    device, get_device_queue(device, 0, queue_index), messenger
end


text:
# LibFTD2XX.jl
#
# By Reuben Hill 2019, Gowerlabs Ltd, reuben@gowerlabs.co.uk
#
# Copyright (c) Gowerlabs Ltd.

using LibFTD2XX

include("util.jl")

numdevs = LibFTD2XX.createdeviceinfolist()
if numdevs > 0
  @info "found $numdevs devices. Running hardware tests..."
  include("hardware/wrapper.jl")
  include("hardware/LibFTD2XX.jl")
else
  @info "found $numdevs devices. Running nohardware tests..."
  include("nohardware/wrapper.jl")
  include("nohardware/LibFTD2XX.jl")
end


text:
module PackageName

# Write your package code here.

print("Hello Julia, Git and Kraken!")

end


text:
"""
    blocklossfn(predblock, yblock)

Construct a loss function that compares a batch of model outputs
`ŷs` and encoded targets `ys` and returns a scalar loss.

For example for `block = OneHotTensor{1}(classes)` (i.e. an encoded
`Label(classes)`), we have
`blocklossfn(block, block) == Flux.Losses.logitcrossentropy`.
"""
function blocklossfn end

blocklossfn(predblock) = blocklossfn(predblock, predblock)

function blocklossfn(outblock::OneHotTensor{0}, yblock::OneHotTensor{0})
    outblock.classes == yblock.classes || error("Classes of $outblock and $yblock differ!")
    return Flux.Losses.logitcrossentropy
end

function blocklossfn(outblock::OneHotTensorMulti{0}, yblock::OneHotTensorMulti{0})
    outblock.classes == yblock.classes || error("Classes of $outblock and $yblock differ!")
    return Flux.Losses.logitbinarycrossentropy
end

function blocklossfn(outblock::OneHotTensor{N}, yblock::OneHotTensor{N}) where N
    outblock.classes == yblock.classes || error("Classes of $outblock and $yblock differ!")
    return segmentationloss
end

function segmentationloss(ypreds, ys; kwargs...)
    # Has to be reshaped to 3D array since `logitcrossentropy(...; dims = 3)`
    # doesn't work on GPU
    sz = size(ypreds)
    ypreds = reshape(ypreds, :, sz[end-1], sz[end])
    ys = reshape(ys, :, size(ys, 3), size(ys, 4))
    Flux.Losses.logitcrossentropy(ypreds, ys; dims = 2, kwargs...)
end

function blocklossfn(outblock::KeypointTensor{N}, yblock::KeypointTensor{N}) where N
    outblock.sz == yblock.sz || error("Sizes of $outblock and $yblock differ!")
    return Flux.Losses.mse
end


text:
using LPWriter
using Base.Test

include("writer.jl")
include("reader.jl")


text:
import Pkg
Pkg.activate(joinpath(@__DIR__, "../../"))  # environment: CytofResearch

using BSON
using CSV
using CytofResearch
import DataFrames; const DF = DataFrames
using Distributions
import Random

max_stars_repo_name: luiarthur/CytofResearch

max_issues_repo_name: luiarthur/CytofResearch

max_forks_repo_name: luiarthur/CytofResearch


text:
export Hankel

struct Hankel{T} <: AbstractArray{T, 2}
    c :: Vector{T}
end
#Hankel{T}(c::Vector{T}) = length(c) % 2 == 1 ? Hankel{T}(c) : throw(ArgumentError(""))

getindex(H::Hankel, i::Int, j::Int) = H.c[i+j-1]
isassigned(H::Hankel, i::Int, j::Int) = isassigned(H.c, i+j-1)

size(H::Hankel, r::Int) = (r==1 || r==2) ? 1 + div(length(H.c),2) :
    throw(ArgumentError("Invalid dimension $r"))
size(H::Hankel) = size(H,1), size(H,2)

# Fast matrix x vector via fft()
function *(A::Hankel{T},x::Vector{T}) where T
    Toeplitz(A.c)*reverse(x)
end

function mul!(y::StridedVector{T},A::Hankel{T},x::StridedVector{T}) where T
    xx=reverse(x)
    return mul!(y,Toeplitz(A.c),xx)
end

function Matrix(H::Hankel{T}) where T
    n=size(H, 1)
    M=Array{T}(undef, n, n)
    for i=1:n
        M[:,i] = H.c[i:i+n-1]
    end
    M
end


text:
import Random

@inline function vote(Θ,x,c,k)

   s = 0
   for j=1:k
       s += c[j]*sign(Θ[j]'*x) # voting (+1 or -1 * c[j] weight)
   end
   s
end

function trainer(model::VotedPerceptron{T},
                 X::AbstractArray{T},
					  Y::Vector{T}) where T<:AbstractFloat

   shuffle_epoch = model.shuffle_epoch
   random_state  = model.random_state
   max_epochs    = model.max_epochs

   if random_state!=-1
      Random.seed!(random_state)
   end

   n,m         = size(X)
   X           = hcat(X,ones(n,1)) # adding bias
   history     = []
   nerrors,nlast_errors = Inf,0
   epochs      = 0
   k,Θ,c,α     = 1,Dict(1=>Random.rand(m+1)),Dict(1=>0),model.α
   #while  nerrors>0 && epochs < max_epochs
   while  epochs < max_epochs
   # stops when error is equal to zero or grater than last_error or reached max iterations
       # shuffle dataset
       if shuffle_epoch
          sind = Random.shuffle(1:n)
          x = X[sind,:]
          y = Y[sind]
       end
       nerrors = 0
       # weight updates for all samples
       for i=1:n
          xi = x[i,:]
          ξ   = sinal(Θ[k]'*xi) - y[i]
          if ξ==0
             c[k] += 1
          else
             nerrors+=1
             c[k+1] = 1
			    Θ[k+1] = Θ[k] - α * ξ * xi
             k     += 1
          end
		 end
       nlast_errors   = nerrors
       epochs+=1
       push!(history,nerrors)
   end
   if nerrors > 0
      warn("Perceptron: Not converged. Max epochs $(max_epochs) reached. Error history: $(history) \n Try to increase max_epochs or may be you have a non linear problem.")
   end
   model.Θ = Θ
   model.c = c
   model.k = k
   model.history = history
end

function predictor(model::VotedPerceptron{T},
	                X::AbstractArray{T}) where T<:AbstractFloat

   Θ = model.Θ
   α = model.α
   k = model.k
   c = model.c

   n   = size(X,1)
   y   = zeros(Real,n)
   X   = hcat(X,ones(n,1)) # adding bias
   for i=1:n
      y[i] = sinal(vote(Θ,X[i,:],c,k))
   end
   y

end


text:
module RTESolver

using StaticArrays
using CUDA
using Adapt

using ..Device: array_type, array_device, CUDADevice, CPU
using ..AngularDiscretizations
using ..Vmrs
using ..AtmosphericStates
using ..Fluxes
using ..BCs
using ..Sources
using ..RTE
using ..Optics
using ..LookUpTables

export solve_lw!, solve_sw!

include("RTESolverKernels.jl")

"""
    solve_lw!(
        slv::Solver{FT,I},
        max_threads::I,
        lkp_args...,
    ) where {I<:Int,FT<:AbstractFloat}

Solver for the longwave radiation problem
"""
function solve_lw!(
    slv::Solver{FT,I},
    max_threads::I,
    lkp_args...,
) where {I<:Int,FT<:AbstractFloat}
    (; as, op, bcs_lw, src_lw, flux_lw, fluxb_lw) = slv
    (; nlay, ncol) = as
    DA = array_type()
    nargs = length(lkp_args)
    @assert nargs < 3
    if nargs > 0
        @assert lkp_args[1] isa LookUpLW
        if nargs == 2
            @assert lkp_args[2] isa LookUpCld
        end
        (; n_gpt, major_gpt2bnd) = lkp_args[1]
        flux = fluxb_lw
    else
        n_gpt = 1
        major_gpt2bnd = DA(UInt8[1])
        flux = flux_lw
    end

    if op isa OneScalar
        rte_lw_noscat_solve!(
            flux,
            flux_lw,
            src_lw,
            bcs_lw,
            op,
            nlay,
            ncol,
            major_gpt2bnd,
            max_threads,
            as,
            lkp_args...,
        ) # no-scattering solver
    else
        rte_lw_2stream_solve!(
            flux,
            flux_lw,
            src_lw,
            bcs_lw,
            op,
            nlay,
            ncol,
            major_gpt2bnd,
            max_threads,
            as,
            lkp_args...,
        ) # 2-stream solver
    end
    return nothing
end

"""
    rte_lw_noscat_solve!(
        flux::FluxLW{FT},
        flux_lw::FluxLW{FT},
        src_lw::SourceLWNoScat{FT},
        bcs_lw::LwBCs{FT},
        op::OneScalar{FT},
        nlay,
        ncol,
        major_gpt2bnd::AbstractArray{UInt8,1},
        max_threads,
        as::AbstractAtmosphericState{FT},
        lkp_args...,
    ) where {FT<:AbstractFloat}

No scattering solver for the longwave problem.
LW fluxes, no scattering, mu (cosine of integration angle) specified by column
Does radiation calculation at user-supplied angles; converts radiances to flux
using user-supplied weights. Currently, only n_gauss_angles = 1 supported.
"""
function rte_lw_noscat_solve!(
    flux::FluxLW{FT},
    flux_lw::FluxLW{FT},
    src_lw::SourceLWNoScat{FT},
    bcs_lw::LwBCs{FT},
    op::OneScalar{FT},
    nlay,
    ncol,
    major_gpt2bnd::AbstractArray{UInt8,1},
    max_threads,
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    n_gpt = length(major_gpt2bnd)
    device = array_device(op.τ)
    if device === CUDADevice() # launch CUDA kernel
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (
            flux,
            flux_lw,
            src_lw,
            bcs_lw,
            op,
            nlay,
            ncol,
            major_gpt2bnd,
            as,
            lkp_args...,
        )
        @cuda threads = (tx) blocks = (bx) rte_lw_noscat_solve_CUDA!(args...)
    else # use Julia native multithreading
        for igpt = 1:n_gpt
            @inbounds Threads.@threads for gcol = 1:ncol
                igpt == 1 && set_flux_to_zero!(flux_lw, gcol)
                compute_optical_props!(op, as, src_lw, gcol, igpt, lkp_args...)
                rte_lw_noscat_source!(src_lw, op, gcol, nlay)
                rte_lw_noscat_transport!(
                    src_lw,
                    bcs_lw,
                    op,
                    gcol,
                    flux,
                    igpt,
                    major_gpt2bnd,
                    nlay,
                    nlev,
                )
                n_gpt > 1 && add_to_flux!(flux_lw, flux, gcol)
            end
        end
    end
    #------------------------------------------------
    return nothing
end
# ---------------------------------------------------------------
function rte_lw_noscat_solve_CUDA!(
    flux::FluxLW{FT},
    flux_lw::FluxLW{FT},
    src_lw::SourceLWNoScat{FT},
    bcs_lw::LwBCs{FT},
    op::OneScalar{FT},
    nlay,
    ncol,
    major_gpt2bnd,
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    n_gpt = length(major_gpt2bnd)
    if gcol ≤ ncol
        for igpt = 1:n_gpt
            igpt == 1 && set_flux_to_zero!(flux_lw, gcol)
            compute_optical_props!(op, as, src_lw, gcol, igpt, lkp_args...)
            rte_lw_noscat_source!(src_lw, op, gcol, nlay)
            rte_lw_noscat_transport!(
                src_lw,
                bcs_lw,
                op,
                gcol,
                flux,
                igpt,
                major_gpt2bnd,
                nlay,
                nlev,
            )
            n_gpt > 1 && add_to_flux!(flux_lw, flux, gcol)
        end
    end
    return nothing
end

"""
    rte_lw_2stream_solve!(
        flux::FluxLW{FT},
        flux_lw::FluxLW{FT},
        src_lw::SourceLW2Str{FT},
        bcs_lw::LwBCs{FT},
        op::TwoStream{FT},
        nlay,
        ncol,
        major_gpt2bnd::AbstractArray{UInt8,1},
        max_threads,
        as::AbstractAtmosphericState{FT},
        lkp_args...,
    ) where {FT<:AbstractFloat}

Two stream solver for the longwave problem.

RRTMGP provides source functions at each level using the spectral mapping
of each adjacent layer. Combine these for two-stream calculations
lw combine sources (rte_lw_2stream_combine_sources!)
combine RRTMGP-specific sources at levels
compute layer reflectance, transmittance
compute total source function at levels using linear-in-tau (rte_lw_2stream_source!)
transport (adding_lw!)
"""
function rte_lw_2stream_solve!(
    flux::FluxLW{FT},
    flux_lw::FluxLW{FT},
    src_lw::SourceLW2Str{FT},
    bcs_lw::LwBCs{FT},
    op::TwoStream{FT},
    nlay,
    ncol,
    major_gpt2bnd::AbstractArray{UInt8,1},
    max_threads,
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    n_gpt = length(major_gpt2bnd)
    device = array_device(op.τ)
    if device === CUDADevice() # launch CUDA kernel
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (
            flux,
            flux_lw,
            src_lw,
            bcs_lw,
            op,
            nlay,
            ncol,
            major_gpt2bnd,
            as,
            lkp_args...,
        )
        @cuda threads = (tx) blocks = (bx) rte_lw_2stream_solve_CUDA!(args...)
    else # use Julia native multithreading
        for igpt = 1:n_gpt
            Threads.@threads for gcol = 1:ncol
                igpt == 1 && set_flux_to_zero!(flux_lw, gcol)
                compute_optical_props!(op, as, src_lw, gcol, igpt, lkp_args...)
                rte_lw_2stream_combine_sources!(src_lw, gcol, nlev, ncol)
                rte_lw_2stream_source!(op, src_lw, gcol, nlay, ncol)
                adding_lw!(
                    flux,
                    src_lw,
                    bcs_lw,
                    gcol,
                    igpt,
                    major_gpt2bnd,
                    nlev,
                    ncol,
                )
                n_gpt > 1 && add_to_flux!(flux_lw, flux, gcol)
            end
        end
    end
    return nothing
end
# -------------------------------------------------------------------------------------------------
function rte_lw_2stream_solve_CUDA!(
    flux::FluxLW{FT},
    flux_lw::FluxLW{FT},
    src_lw::SourceLW2Str{FT},
    bcs_lw::LwBCs{FT},
    op::TwoStream{FT},
    nlay,
    ncol,
    major_gpt2bnd::AbstractArray{UInt8,1},
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    n_gpt = length(major_gpt2bnd)
    if gcol ≤ ncol
        for igpt = 1:n_gpt
            igpt == 1 && set_flux_to_zero!(flux_lw, gcol)
            compute_optical_props!(op, as, src_lw, gcol, igpt, lkp_args...)
            rte_lw_2stream_combine_sources!(src_lw, gcol, nlev, ncol)
            rte_lw_2stream_source!(op, src_lw, gcol, nlay, ncol)
            adding_lw!(
                flux,
                src_lw,
                bcs_lw,
                gcol,
                igpt,
                major_gpt2bnd,
                nlev,
                ncol,
            )
            n_gpt > 1 && add_to_flux!(flux_lw, flux, gcol)
        end
    end
    return nothing
end

"""
    solve_sw!(
        slv::Solver{FT,I},
        max_threads::I,
        lkp_args...,
    ) where {I<:Int,FT<:AbstractFloat}

Solver for the shortwave radiation problem
"""
function solve_sw!(
    slv::Solver{FT,I},
    max_threads::I,
    lkp_args...,
) where {I<:Int,FT<:AbstractFloat}
    (; as, op, bcs_sw, src_sw, flux_sw, fluxb_sw) = slv
    (; nlay, ncol) = as
    DA = array_type()
    nargs = length(lkp_args)
    @assert nargs < 3
    if nargs > 0
        @assert lkp_args[1] isa LookUpSW
        if nargs == 2
            @assert lkp_args[2] isa LookUpCld
        end
        (; n_gpt, major_gpt2bnd, solar_src_scaled) = lkp_args[1]
        flux = fluxb_sw
    else
        n_gpt = 1
        major_gpt2bnd = DA(UInt8[1])
        solar_src_scaled = DA(FT[1])
        flux = flux_sw
    end
    # solving radiative transfer equation
    if slv.op isa OneScalar
        rte_sw_noscat_solve!(
            flux,
            flux_sw,
            op,
            bcs_sw,
            nlay,
            ncol,
            solar_src_scaled,
            max_threads,
            as,
            lkp_args...,
        ) # no-scattering solver
    else
        rte_sw_2stream_solve!(
            flux,
            flux_sw,
            op,
            bcs_sw,
            src_sw,
            nlay,
            ncol,
            major_gpt2bnd,
            solar_src_scaled,
            max_threads,
            as,
            lkp_args...,
        ) # 2-stream solver
    end
    return nothing
end

"""
    rte_sw_noscat_solve!(
        flux::FluxSW{FT},
        flux_sw::FluxSW{FT},
        op::OneScalar{FT},
        bcs_sw::SwBCs{FT},
        nlay,
        ncol,
        solar_src_scaled::AbstractArray{FT,1},
        max_threads,
        as::AbstractAtmosphericState{FT},
        lkp_args...,
    ) where {FT<:AbstractFloat}

No-scattering solver for the shortwave problem.
(Extinction-only i.e. solar direct beam)
"""
function rte_sw_noscat_solve!(
    flux::FluxSW{FT},
    flux_sw::FluxSW{FT},
    op::OneScalar{FT},
    bcs_sw::SwBCs{FT},
    nlay,
    ncol,
    solar_src_scaled::AbstractArray{FT,1},
    max_threads,
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    n_gpt = length(solar_src_scaled)
    device = array_device(op.τ)
    # setting references for flux_sw
    if device === CUDADevice() # launching CUDA kernel
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (
            flux,
            flux_sw,
            op,
            bcs_sw,
            nlay,
            ncol,
            solar_src_scaled,
            as,
            lkp_args...,
        )
        @cuda threads = (tx) blocks = (bx) rte_sw_noscat_solve_CUDA!(args...)
    else # launch juila native multithreading
        for igpt = 1:n_gpt
            @inbounds Threads.@threads for gcol = 1:ncol
                igpt == 1 && set_flux_to_zero!(flux_sw, gcol)
                compute_optical_props!(op, as, gcol, igpt, lkp_args...)
                rte_sw_noscat_solve_kernel!(
                    flux,
                    op,
                    bcs_sw,
                    igpt,
                    solar_src_scaled,
                    gcol,
                    nlev,
                )
                n_gpt > 1 && add_to_flux!(flux_sw, flux, gcol)
            end
        end
    end
    #-----------------------------------------------------------------------------------
    return nothing
end
#--------------------------------------------------------------------------------------------------
function rte_sw_noscat_solve_CUDA!(
    flux::FluxSW{FT},
    flux_sw::FluxSW{FT},
    op::OneScalar{FT},
    bcs_sw::SwBCs{FT},
    nlay,
    ncol,
    solar_src_scaled::AbstractArray{FT,1},
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    n_gpt = length(solar_src_scaled)
    # setting references for flux_sw
    if gcol ≤ ncol
        for igpt = 1:n_gpt
            igpt == 1 && set_flux_to_zero!(flux_sw, gcol)
            compute_optical_props!(op, as, gcol, igpt, lkp_args...)
            rte_sw_noscat_solve_kernel!(
                flux,
                op,
                bcs_sw,
                igpt,
                solar_src_scaled,
                gcol,
                nlev,
            )
            n_gpt > 1 && add_to_flux!(flux_sw, flux, gcol)
        end
    end
    return nothing
end
#--------------------------------------------------------------------------------------------------
"""
    rte_sw_2stream_solve!(
        flux::FluxSW{FT},
        flux_sw::FluxSW{FT},
        op::TwoStream{FT},
        bcs_sw::SwBCs{FT},
        src_sw::SourceSW2Str{FT},
        nlay,
        ncol,
        major_gpt2bnd::AbstractArray{UInt8,1},
        solar_src_scaled::AbstractArray{FT,1},
        max_threads,
        as::AbstractAtmosphericState{FT},
        lkp_args...,
    ) where {FT<:AbstractFloat}

Two stream solver for the shortwave problem.
"""
function rte_sw_2stream_solve!(
    flux::FluxSW{FT},
    flux_sw::FluxSW{FT},
    op::TwoStream{FT},
    bcs_sw::SwBCs{FT},
    src_sw::SourceSW2Str{FT},
    nlay,
    ncol,
    major_gpt2bnd::AbstractArray{UInt8,1},
    solar_src_scaled::AbstractArray{FT,1},
    max_threads,
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    nlev = nlay + 1
    n_gpt = length(major_gpt2bnd)
    device = array_device(op.τ)
    if device === CUDADevice()
        tx = min(ncol, max_threads)
        bx = cld(ncol, tx)
        args = (
            flux,
            flux_sw,
            op,
            bcs_sw,
            src_sw,
            nlay,
            ncol,
            major_gpt2bnd,
            solar_src_scaled,
            as,
            lkp_args...,
        )
        @cuda threads = (tx) blocks = (bx) rte_sw_2stream_solve_CUDA!(args...)
    else # launch julia native multithreading
        # setting references for flux_sw
        for igpt = 1:n_gpt
            @inbounds Threads.@threads for gcol = 1:ncol
                igpt == 1 && set_flux_to_zero!(flux_sw, gcol)
                compute_optical_props!(op, as, gcol, igpt, lkp_args...)
                sw_two_stream!(op, src_sw, bcs_sw, gcol, nlay) # Cell properties: transmittance and reflectance for direct and diffuse radiation
                sw_source_2str!(
                    src_sw,
                    bcs_sw,
                    gcol,
                    flux,
                    igpt,
                    solar_src_scaled,
                    major_gpt2bnd,
                    nlay,
                ) # Direct-beam and source for diffuse radiation
                adding_sw!(
                    src_sw,
                    bcs_sw,
                    gcol,
                    flux,
                    igpt,
                    major_gpt2bnd,
                    nlev,
                ) # Transport
                n_gpt > 1 && add_to_flux!(flux_sw, flux, gcol)
            end
        end
    end
    return nothing
end
#--------------------------------------------------------------------------------------------------
function rte_sw_2stream_solve_CUDA!(
    flux::FluxSW{FT},
    flux_sw::FluxSW{FT},
    op::TwoStream{FT},
    bcs_sw::SwBCs{FT},
    src_sw::SourceSW2Str{FT},
    nlay,
    ncol,
    major_gpt2bnd::AbstractArray{UInt8,1},
    solar_src_scaled::AbstractArray{FT,1},
    as::AbstractAtmosphericState{FT},
    lkp_args...,
) where {FT<:AbstractFloat}
    gcol = threadIdx().x + (blockIdx().x - 1) * blockDim().x # global id
    nlev = nlay + 1
    n_gpt = length(major_gpt2bnd)
    if gcol ≤ ncol
        for igpt = 1:n_gpt
            igpt == 1 && set_flux_to_zero!(flux_sw, gcol)
            compute_optical_props!(op, as, gcol, igpt, lkp_args...)
            sw_two_stream!(op, src_sw, bcs_sw, gcol, nlay) # Cell properties: transmittance and reflectance for direct and diffuse radiation
            sw_source_2str!(
                src_sw,
                bcs_sw,
                gcol,
                flux,
                igpt,
                solar_src_scaled,
                major_gpt2bnd,
                nlay,
            ) # Direct-beam and source for diffuse radiation
            adding_sw!(src_sw, bcs_sw, gcol, flux, igpt, major_gpt2bnd, nlev) # Transport
            n_gpt > 1 && add_to_flux!(flux_sw, flux, gcol)
        end
    end
    return nothing
end
#--------------------------------------------------------------------------------------------------
end


text:
export PhspVector
using Mmap: Mmap

"""
PhspVector

A read only memory mapped Phase Space vector.
The user is responsible for closing the `io`.
"""
struct PhspVector{P,H} <: AbstractVector{P}
    io::IOStream
    header::H
    array::Vector{P}
    # owns_io::Bool
    function PhspVector{P,H}(io::IO, header::H) where {P,H}
        offset = if header isa EGSHeader
            sizeof(P)
        else
            0
        end
        dims = (get_nparticles(io, header),)
        seekstart(io)
        array = Mmap.mmap(io, Vector{P}, dims, offset)
        new{P, H}(io, header, array)
    end
end

function get_nparticles_header(h::EGSHeader)
    h.particlecount
end
function get_nparticles_header(h::IAEAHeader)
    parse(Int64, h.attributes[:PARTICLES])
end

function get_nparticles_iosize(P, io::IO)::Float64
    total_size  = bytelength(io)
    ret = (total_size / sizeof(P))
    if P <: EGSParticle
        ret -= 1 # first particle encodes the header
    end
    ret
end

function get_nparticles(io::IO, header)
    P = compressed_ptype(header)
    nheader = get_nparticles_header(header)
    ffilesize = get_nparticles_iosize(P, io)
    nfilesize = floor(Int64, ffilesize)
    if ffilesize != nfilesize
        @warn "File size indicates $(ffilesize) particles, which is not an integer. Assuming $(nfilesize) particles instead."
    end
    if nfilesize != nheader
        @warn "Particle count according to the header is $(nheader), while there are actually $(nfilesize) particles stored in the file."
    end
    nfilesize
end

function PhspVector(io::IO, h::Union{EGSHeader, IAEAHeader})
    P = compressed_ptype(h)
    H = typeof(h)
    PhspVector{P,H}(io, h)
end

function compressed_ptype(h::EGSHeader)
    ptype(h)
end
function compressed_ptype(h::IAEAHeader{Nf, Ni}) where {Nf, Ni}
    ptype(StaticIAEAHeader(h.record_contents, Nf, Ni))
end

function _PhspVector(path::AbstractString, fmt::FormatIAEA)
    _PhspVector(IAEAPath(path), fmt)
end

function _PhspVector(path::IAEAPath, ::FormatIAEA)
    io = open(path.phsp)
    h = load(path.header, IAEAHeader)
    PhspVector(io, h)
end

function _PhspVector(path::AbstractString, ::FormatEGS)
    io = open(path)
    h = consume_egs_header(io)
    seekstart(io)
    PhspVector(io, h)
end

function PhspVector(path::AbstractString)
    fmt = guess_format_from_path(path)
    _PhspVector(path, fmt)
end

function PhspVector(path::IAEAPath)
    _PhspVector(path, FormatIAEA())
end

for (f, args) in [
        [:(Base.size), ()],
        [:(Base.getindex), (:(i::Union{Integer, CartesianIndex}),)],
    ]
    @eval $f(o::PhspVector, $(args...)) = $f(o.array, $(args...))
end

Base.getindex(o::PhspVector, inds...) = view(o, inds...)


text:
module LinearSolve

using ArrayInterface: lu_instance
using Base: cache_dependencies, Bool
using Krylov
using LinearAlgebra
using Reexport
using SciMLBase: AbstractDiffEqOperator, AbstractLinearAlgorithm
using Setfield
using UnPack

@reexport using SciMLBase

abstract type SciMLLinearSolveAlgorithm end

include("common.jl")
include("factorization.jl")
include("krylov.jl")

export LUFactorization, SVDFactorization, QRFactorization
export KrylovJL

end


text:
module FixedArrayBench
# TODO, make everything work automatically and save history!

using BenchmarkTools
import FixedSizeArrays
import SIMD
import StaticArrays

include("missing_ops.jl")
const suite = BenchmarkGroup()
include("matrix.jl")
include("vector.jl")
# Define a parent BenchmarkGroup to contain our suite

function run(tunecache=true)
    vector_bench()
    if tunecache && isfile("params.jld")
        loadparams!(suite, JLD.load("params.jld", "suite"), :evals, :samples)
    else
        println("tuning...")
        tune!(suite)
        println("done tuning")
    end
    result = BenchmarkTools.run(suite)
    #JLD.save("result$(Dates.format(now(), "yyudH")).jld", "suite", result)
    result
end

end # module


text:
using Distributions
using LinearAlgebra
using ProgressMeter

abstract type MonteCarloMethod end

export DataAssimilation

"""
    DataAssimilation( forecasting, method, np, xt, sigma2)

parameters of the filtering method
 - method :chosen method (:AnEnKF, :AnEnKS, :AnPF)
 - N      : number of members (AnEnKF/AnEnKS) or particles (AnPF)
"""
struct DataAssimilation

    xb::Vector{Float64}
    B::Array{Float64,2}
    H::Array{Bool,2}
    R::Array{Float64,2}
    m::AbstractForecasting

    function DataAssimilation(m::AbstractForecasting, xt::TimeSeries, sigma2::Float64)

        xb = xt.u[1]
        B = 0.1 * Matrix(I, xt.nv, xt.nv)
        H = Matrix(I, xt.nv, xt.nv)
        R = sigma2 .* H

        new(xb, B, H, R, m)

    end

    function DataAssimilation(m::StateSpaceModel, xt::TimeSeries)

        xb = xt.u[1]
        B = 0.1 * Matrix(I, xt.nv, xt.nv)
        H = Matrix(I, xt.nv, xt.nv)
        R = m.sigma2_obs .* H

        new(xb, B, H, R, m)

    end

end


text:
function reverseandcombinetext(str::String)::String
    words = split(str,' ')
    # standard word sequence
    index = collect(1:lastindex(words))
    # calculate the reverseandcombine word sequence
    blocksize = 1
    while blocksize < lastindex(words)
        for bindex = 1:blocksize:lastindex(words)
            block = bindex:min(lastindex(words),bindex+blocksize-1)
            index[block] = -reverse(index[block])
        end
        blocksize *= 2
    end
    # combine the words
    buffer = IOBuffer(maxsize=length(str))
    for i in index
        write( buffer, i < 0 ? reverse(words[-i]) : words[i] )
    end
    String(take!(buffer))
end

max_stars_repo_name: UlrichBerntien/Codewars-Katas

max_issues_repo_name: UlrichBerntien/Codewars-Katas

max_forks_repo_name: UlrichBerntien/Codewars-Katas


text:
mutable struct Swipl
    initiated::Bool
end

Swipl() = Swipl(false)


text:
function cuITensor(::Type{T}, inds::IndexSet) where {T<:Number}
  return ITensor(Dense{float(T)}(CUDA.zeros(float(T), dim(inds))), inds)
end
cuITensor(::Type{T}, inds::Index...) where {T<:Number} = cuITensor(T, IndexSet(inds...))

cuITensor(is::IndexSet) = cuITensor(Float64, is)
cuITensor(inds::Index...) = cuITensor(IndexSet(inds...))

cuITensor() = ITensor()
function cuITensor(x::S, inds::IndexSet{N}) where {S<:Number,N}
  dat = CuVector{float(S)}(undef, dim(inds))
  fill!(dat, float(x))
  return ITensor(Dense{S}(dat), inds)
end
cuITensor(x::S, inds::Index...) where {S<:Number} = cuITensor(x, IndexSet(inds...))

#TODO: check that the size of the Array matches the Index dimensions
function cuITensor(A::Array{S}, inds) where {S<:Number}
  return ITensor(Dense(CuArray{S}(A)), inds)
end
function cuITensor(A::CuArray{S}, inds::IndexSet) where {S<:Number}
  return ITensor(Dense(A), inds)
end
cuITensor(A::Array{S}, inds::Index...) where {S<:Number} = cuITensor(A, IndexSet(inds...))
cuITensor(A::CuArray{S}, inds::Index...) where {S<:Number} = cuITensor(A, IndexSet(inds...))
function cuITensor(A::ITensor)
  return if storage(tensor(A)) isa ITensors.EmptyStorage
    cuITensor(zero(eltype(storage(tensor(A)))), inds(A)...)
  else
    cuITensor(data(tensor(A)), inds(A)...)
  end
end

cu(A::ITensor) = cuITensor(A)

function cpu(A::ITensor)
  typeof(data(storage(A))) <: CuArray && return ITensor(cpu(storage(A)), inds(A))
  return A
end

function randomCuITensor(::Type{S}, inds::Indices) where {S<:Real}
  T = cuITensor(S, inds)
  randn!(T)
  return T
end
function randomCuITensor(::Type{S}, inds::Indices) where {S<:Complex}
  Tr = cuITensor(real(S), inds)
  Ti = cuITensor(real(S), inds)
  randn!(Tr)
  randn!(Ti)
  return complex(Tr) + im * Ti
end
function randomCuITensor(::Type{S}, inds::Index...) where {S<:Number}
  return randomCuITensor(S, IndexSet(inds...))
end
randomCuITensor(inds::IndexSet) = randomCuITensor(Float64, inds)
randomCuITensor(inds::Index...) = randomCuITensor(Float64, IndexSet(inds...))

CuArray(T::ITensor) = CuArray(tensor(T))

function CuArray{ElT,N}(T::ITensor, is::Vararg{Index,N}) where {ElT,N}
  ndims(T) != N && throw(
    DimensionMismatch(
      "cannot convert an $(ndims(T)) dimensional ITensor to an $N-dimensional CuArray."
    ),
  )
  TT = tensor(permute(T, is...; allow_alias=true))
  return CuArray{ElT,N}(TT)::CuArray{ElT,N}
end

function CuArray{ElT}(T::ITensor, is::Vararg{Index,N}) where {ElT,N}
  return CuArray{ElT,N}(T, is...)
end

function CuArray(T::ITensor, is::Vararg{Index,N}) where {N}
  return CuArray{eltype(T),N}(T, is...)::CuArray{<:Number,N}
end

CUDA.CuMatrix(A::ITensor) = CuArray(A)

function CuVector(A::ITensor)
  if ndims(A) != 1
    throw(DimensionMismatch("Vector() expected a 1-index ITensor"))
  end
  return CuArray(A)
end

function CuMatrix(T::ITensor, i1::Index, i2::Index)
  ndims(T) != 2 &&
    throw(DimensionMismatch("ITensor must be order 2 to convert to a Matrix"))
  return CuArray(T, i1, i2)
end


text:
## auxiliary

const lowest = v"0"
const highest = v"999"

# some version comparisons need to ignore part of the version number
strip_patch(ver) = VersionNumber(ver.major, ver.minor)
strip_minor(ver) = VersionNumber(ver.major)

## devices supported by the CUDA toolkit

# Source:
# - https://en.wikipedia.org/wiki/CUDA#GPUs_supported
# - ptxas |& grep -A 10 '\--gpu-name'
const dev_cuda_db = Dict(
    v"1.0" => lowest:v"6.5",
    v"1.1" => lowest:v"6.5",
    v"1.2" => lowest:v"6.5",
    v"1.3" => lowest:v"6.5",
    v"2.0" => lowest:v"8.0",
    v"2.1" => lowest:v"8.0",
    v"3.0" => v"4.2":highest,
    v"3.2" => v"6.0":highest,
    v"3.5" => v"5.0":highest,
    v"3.7" => v"6.5":highest,
    v"5.0" => v"6.0":highest,
    v"5.2" => v"7.0":highest,
    v"5.3" => v"7.5":highest,
    v"6.0" => v"8.0":highest,
    v"6.1" => v"8.0":highest,
    v"6.2" => v"8.0":highest,
    v"7.0" => v"9.0":highest,
    v"7.2" => v"9.2":highest,
    v"7.5" => v"10.0":highest,
)

function devices_for_cuda(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in dev_cuda_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## PTX ISAs supported by the CUDA toolkit

# Source:
# - PTX ISA document, Release History table
const isa_cuda_db = Dict(
    v"1.0" => v"1.0":highest,
    v"1.1" => v"1.1":highest,
    v"1.2" => v"2.0":highest,
    v"1.3" => v"2.1":highest,
    v"1.4" => v"2.2":highest,
    v"1.5" => v"2.2":highest,
    v"2.0" => v"3.0":highest,
    v"2.1" => v"3.1":highest,
    v"2.2" => v"3.2":highest,
    v"2.3" => v"4.2":highest,
    v"3.0" => v"4.1":highest,
    v"3.1" => v"5.0":highest,
    v"3.2" => v"5.5":highest,
    v"4.0" => v"6.0":highest,
    v"4.1" => v"6.5":highest,
    v"4.2" => v"7.0":highest,
    v"4.3" => v"7.5":highest,
    v"5.0" => v"8.0":highest,
    v"6.0" => v"9.0":highest,
    v"6.1" => v"9.1":highest,
    v"6.2" => v"9.2":highest,
    v"6.3" => v"10.0":highest,
    v"6.4" => v"10.1":highest,
)

function isas_for_cuda(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in isa_cuda_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## devices supported by the LLVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const dev_llvm_db = Dict(
    v"2.0" => v"3.2":highest,
    v"2.1" => v"3.2":highest,
    v"3.0" => v"3.2":highest,
    v"3.2" => v"3.7":highest,
    v"3.5" => v"3.2":highest,
    v"3.7" => v"3.7":highest,
    v"5.0" => v"3.5":highest,
    v"5.2" => v"3.7":highest,
    v"5.3" => v"3.7":highest,
    v"6.0" => v"3.9":highest,
    v"6.1" => v"3.9":highest,
    v"6.2" => v"3.9":highest,
    v"7.0" => v"6.0":highest,
    v"7.2" => v"7.0":highest,
    v"7.5" => v"8.0":highest,
)

function devices_for_llvm(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in dev_llvm_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## PTX ISAs supported by the LVM NVPTX back-end

# Source: LLVM/lib/Target/NVPTX/NVPTX.td
const isa_llvm_db = Dict(
    v"3.0" => v"3.2":v"3.5",
    v"3.1" => v"3.2":v"3.5",
    v"3.2" => v"3.5":highest,
    v"4.0" => v"3.5":highest,
    v"4.1" => v"3.7":highest,
    v"4.2" => v"3.7":highest,
    v"4.3" => v"3.9":highest,
    v"5.0" => v"3.9":highest,
    v"6.0" => v"6.0":highest,
    v"6.1" => v"7.0":highest,
    v"6.3" => v"8.0":highest,
    v"6.4" => v"9.0":highest,
)

function isas_for_llvm(ver::VersionNumber)
    caps = Set{VersionNumber}()
    for (cap,r) in isa_llvm_db
        if strip_patch(ver) in r
            push!(caps, cap)
        end
    end
    return caps
end

## other

shader(cap::VersionNumber) = "sm_$(cap.major)$(cap.minor)"


text:
struct Label <: Drawable
    text::String
    size::Int8
    color::Color
    offset::Offset
end

function Label(text::String; color="#333", size=22)
    Label(
        text,
        size,
        Color(color),
        Offset(0,0),
    )
end

Label(::Nothing) = nothing
Label(c::Char) = Label(string(c))
Label(l::Label) = l
Label(x::Number) = Label(string(x))

function draw(l::Label, pos, cell_width; halign=:left, valign=:middle)
    fontsize(l.size)
    sethue(l.color.name)
    text(
        l.text,
        pos,
        halign=halign,
        valign=valign
    )
end


text:
using PartialP, Test
using PartialP: gradient

if PartialP.usetyped
  @info "Testing PartialP in type-hacks mode."
else
  @info "Testing PartialP in normal mode."
end

@testset "PartialP" begin

@testset "Features" begin
  include("features.jl")
end

@testset "Gradients" begin
  include("gradcheck.jl")
end

@testset "Complex" begin
  include("complex.jl")
end

@testset "Compiler" begin
  include("compiler.jl")
end

end


text:
module Green_Utils

using Cubature: hcubature, hquadrature
using QuadGK

# Defining several global variables
mutable struct HFVec # <--------------------------------------------------------------- Order of declaration IS IMPORTANT
    Niωn_::Int64 ## Discrtization of the Matsubara_Grid
    beta_::Int64 ## Temperature for the summation
    N_it_::Int64 ## Number of iterations
    dict_v_::Dict{String,Float64} ## Dict containing U and V
    Matsubara_G_::Vector{Complex{Float64}} ## Matsubara grid

    function HFVec(Niωn_::Int64, beta_::Int64, N_it_::Int64, dict_v_::Dict{String,Float64})
        Matsubara_G_ = Matsubara_Grid(beta_,Niωn_)

        return new(Niωn_, beta_, N_it_, dict_v_, Matsubara_G_)
    end
end

function Matsubara_Grid(beta_::Int64, Niωn_::Int64)
    Matsubara_G_ = Vector{Complex{Float64}}(undef,Niωn_)
    Niωn_ = mod(Niωn_,2) == 0 ? Niωn_ : Niωn_ + 1
    ## The information in the first Matsubara frequencies are more important
    Matsubara_G_[1:end] = 1.0im*[(2*n + 1)*pi/(beta_) for n in 1:Niωn_]

    return Matsubara_G_
end

struct Integral
    kx::Float64
    ky::Float64
    iωn::Complex{Float64}
end

struct Integral1D
    kx::Float64
    iωn::Complex{Float64}
end

# Defining some functions
function Potential(qx::Float64, qy::Float64, U::Float64, V::Float64)
    UV = U + V*(cos(qx) + cos(qy))

    return UV
end

function epsilonk(kx::Float64, ky::Float64; t::Float64=1.0, tp::Float64=-0.3, tpp::Float64=0.2)
    epsilonk1 = -2.0*t*(cos(kx) + cos(ky))
    epsilonk2 = -2.0*tp*(cos(kx+ky)+cos(kx-ky))
    epsilonk3 = -2.0*tpp*(cos(2.0*kx)+cos(2.0*ky))

    return epsilonk1
end

function Potential1D(qx::Float64, U::Float64, V::Float64)
    UV = U + V*(cos(qx))

    return UV
end

function epsilonk1D(kx::Float64; t::Float64=1.0)
    epsilonk1 = -2.0*t*(cos(kx))

    return epsilonk1
end

function Gk(qx::Float64, qy::Float64, kx::Float64, ky::Float64, iωn::Complex{Float64}, SE_func::Function)
    integ = Integral(kx + qx, ky + qy, iωn)
    gk = 1.0/(iωn - epsilonk(kx + qx, ky + qy) - SE_func(integ))

    return gk

end

function Init_Gk(qx::Float64, qy::Float64, kx::Float64, ky::Float64, iωn::Complex{Float64})
    gk_init = 1.0/(iωn - epsilonk(kx + qx, ky + qy) - (0.1-0.1im))

    return gk_init
end

function Gk1D(qx::Float64, kx::Float64, iωn::Complex{Float64}, SE_func::Function)
    integ = Integral1D(kx + qx, iωn)
    gk = 1.0/(iωn - epsilonk1D(kx + qx) - SE_func(integ))

    return gk

end

function Init_Gk1D(qx::Float64, kx::Float64, iωn::Complex{Float64})
    gk_init = 1.0/(iωn - epsilonk1D(kx + qx) - (0.1-0.1im))

    return gk_init
end

function Funct_Wrapper_one(funct::Function)
    function Inner_funct(k::Array{Float64,1})
        return funct(k[1], k[2])
    end

    return Inner_funct
end

function Funct_Wrapper_two(funct::Function, other::Integral)
    function Inner_funct(k::Array{Float64,1})
        return funct(k[1], k[2], other)
    end

    return Inner_funct
end

function Funct_Wrapper_two1D(funct::Function, other::Integral1D)
    function Inner_funct(k::Float64)
        return funct(k, other)
    end

    return Inner_funct
end

function integrate_complex(funct::Function, SE_funct::Function, ii::Int64, st::HFVec, BoundArr::Array{Array{Float64,1},1})
    dict = st.dict_v_
    if ii <= 1
        println(ii, " integrate_complex")
        dressed_funct = (qx::Float64, qy::Float64, kx::Float64, ky::Float64, iωn::Complex{Float64})->1.0*Potential(qx, qy, values(dict)...)*funct(qx, qy, kx ,ky, iωn)
    elseif ii > 1
        println(ii, " integrate_complex")
        dressed_funct = (qx::Float64, qy::Float64, kx::Float64, ky::Float64, iωn::Complex{Float64})->1.0*Potential(qx, qy, values(dict)...)*funct(qx, qy, kx, ky, iωn, SE_funct)
    end

    function real_funct(qx::Float64, qy::Float64, rest::Integral)
        return real(dressed_funct(qx, qy, rest.kx, rest.ky, rest.iωn))
    end

    function imag_funct(qx::Float64, qy::Float64, rest::Integral)
        return imag(dressed_funct(qx, qy, rest.kx, rest.ky, rest.iωn))
    end

    result_real = (remaining_var::Integral) -> hcubature(Funct_Wrapper_two(real_funct,remaining_var), BoundArr[1], BoundArr[2]; reltol=1.5e-2, abstol=1.5e-2, maxevals=100_000)[1]
    result_imag = (remaining_var::Integral) -> hcubature(Funct_Wrapper_two(imag_funct,remaining_var), BoundArr[1], BoundArr[2]; reltol=1.5e-2, abstol=1.5e-2, maxevals=100_000)[1]

    temp_func = (KK::Integral) -> 2.0*(2.0*pi)^(-2.0)*(result_real(KK) + 1.0im*result_imag(KK))

    return temp_func
end

function interation_process(st::HFVec, BoundArr::Array{Array{Float64,1},1}, SubLast::Int64)
    N_it = st.N_it_
    SE_funct = nothing
    Dict_array = Dict{Int64,Array{Function,1}}()
    for it in 1:N_it
        Funct_array = Array{Function,1}([])
        if it <= 1
            for sub in 1:SubLast
                println(it, " interation_process")
                #SE_funct = integrate_complex(Init_Gk, x -> x, it, st, BoundArr) ## SE_funct takes in Integral-type values
                #push!(Funct_array,SE_funct)
                to_add_lower = (convert(Float64,sub)-1.)*2.0*pi/SubLast
                to_add_upper = (convert(Float64,sub))*2.0*pi/SubLast
                BoundArr = Array{Array{Float64,1},1}([[-pi+to_add_lower,-pi+to_add_lower],[-pi+to_add_upper,-pi+to_add_upper]])
                println("it: ", it, "BoundArr: ", BoundArr)
                push!(Funct_array, integrate_complex(Init_Gk, x -> x, it, st, BoundArr))
            end
        elseif it > 1
            for sub in 1:SubLast
                println(it, " interation_process")
                #SE_temp = Array{Function,1}(undef,SubLast)
                to_add_lower = (convert(Float64,sub)-1.)*2.0*pi/SubLast
                to_add_upper = (convert(Float64,sub))*2.0*pi/SubLast
                BoundArr = Array{Array{Float64,1},1}([[-pi+to_add_lower,-pi+to_add_lower],[-pi+to_add_upper,-pi+to_add_upper]])
                println("BoundArr: ", BoundArr)
                push!(Funct_array, integrate_complex(Gk, Dict_array[it-1][sub], it, st, BoundArr))
            end
        end
        Dict_array[it] = Funct_array
    end

    return Dict_array
end

# function funct_wrapper_vector(funct::Function, integ::Integral1D, vec::Vector{Float64})
#     function wrapped_function(k::Float64,vec)
#         return(funct(k[1], integ, vec))
#     end
#     return wrapped_function
# end

function integrate_complex1D(funct::Function, SE_funct::Function, ii::Int64, st::HFVec, BoundArr::Array{Float64,1}, Gridk::Int64; opt::String="sum")
    dict = st.dict_v_
    if ii <= 1
        println(ii, " integrate_complex 1D")
        dressed_funct = (qx::Float64, kx::Float64, iωn::Complex{Float64})->1.0*Potential1D(qx,values(dict)...)*funct(qx, kx, iωn)
    elseif ii > 1
        println(ii, " integrate_complex 1D")
        dressed_funct = (qx::Float64, kx::Float64, iωn::Complex{Float64})->1.0*Potential1D(qx, values(dict)...)*funct(qx, kx, iωn, SE_funct)
    end

    function real_funct(qx::Float64, rest::Integral1D)
        return real(dressed_funct(qx, rest.kx, rest.iωn))
    end

    function imag_funct(qx::Float64, rest::Integral1D)
        return imag(dressed_funct(qx, rest.kx, rest.iωn))
    end

    # vec = Vector{Float64}(undef,2)
    # function vectorized_version(qx::Float64, rest::Integral1D, vec::Vector{Float64})
    #     length(vec) != 2 && throw(ErrorException("Size must be two!!"))
    #     realf = real(dressed_funct(qx, rest.kx, rest.iωn))
    #     imagf = imag(dressed_funct(qx, rest.kx, rest.iωn))
    #     vec[1] = realf; vec[2] = imagf

    #     return nothing
    # end

    # vectorized_funct = (KK::Integral1D) -> hquadrature(2, funct_wrapper_vector(vectorized_version,KK,vec), BoundArr[1], BoundArr[2]; reltol=1.5e-8, abstol=1.5e-8, maxevals=100_000)[1]

    # temp_func = (KK::Integral1D) -> 2.0*(2.0*pi)^(-1.0)*(vectorized_funct(KK)[1]+1.0im*vectorized_funct(KK)[2])

    if opt == "integral"
        result_real = (remaining_var::Integral1D) -> 2.0*(2.0*pi)^(-1.0)*quadgk(Funct_Wrapper_two1D(real_funct,remaining_var), BoundArr[1], BoundArr[2]; rtol=1.5e-4, atol=1.5e-4, maxevals=100_000)[1]
        result_imag = (remaining_var::Integral1D) -> 2.0*(2.0*pi)^(-1.0)*quadgk(Funct_Wrapper_two1D(imag_funct,remaining_var), BoundArr[1], BoundArr[2]; rtol=1.5e-4, atol=1.5e-4, maxevals=100_000)[1]
    elseif opt == "sum"
        result_real = (remaining_var::Integral1D) -> 2.0*(Gridk)^(-1.0)*sum(Funct_Wrapper_two1D(real_funct,remaining_var), range(BoundArr[1], stop=BoundArr[2], length=Gridk))
        result_imag = (remaining_var::Integral1D) -> 2.0*(Gridk)^(-1.0)*sum(Funct_Wrapper_two1D(imag_funct,remaining_var), range(BoundArr[1], stop=BoundArr[2], length=Gridk))
    else
        throw(ErrorException("Oups! opt parameter takes as input only: \"sum\" or \"integral\"."))
    end

    temp_func = (KK::Integral1D) -> 1.0*(result_real(KK) + 1.0im*result_imag(KK))

    return temp_func
end

function interation_process1D(st::HFVec, BoundArr::Array{Float64,1}, SubLast::Int64, Gridk::Int64)
    N_it = st.N_it_
    SE_funct = nothing
    Dict_array = Dict{Int64,Array{Function,1}}()
    for it in 1:N_it
        Funct_array = Array{Function,1}([])
        if it <= 1
            for sub in 1:SubLast
                println(it, " interation_process 1D")
                to_add_lower = (convert(Float64,sub)-1.)*2.0*pi/SubLast
                to_add_upper = (convert(Float64,sub))*2.0*pi/SubLast
                BoundArr = Array{Float64,1}([-pi+to_add_lower,-pi+to_add_upper])
                println("it: ", it, "BoundArr1D: ", BoundArr)
                push!(Funct_array, integrate_complex1D(Init_Gk1D, x -> x, it, st, BoundArr, Gridk))
            end
        elseif it > 1
            for sub in 1:SubLast
                println(it, " interation_process 1D")
                to_add_lower = (convert(Float64,sub)-1.)*2.0*pi/SubLast
                to_add_upper = (convert(Float64,sub))*2.0*pi/SubLast
                BoundArr = Array{Float64,1}([-pi+to_add_lower,-pi+to_add_upper])
                println("BoundArr1D: ", BoundArr)
                push!(Funct_array, integrate_complex1D(Gk1D, Dict_array[it-1][sub], it, st, BoundArr, Gridk))
            end
        end
        Dict_array[it] = Funct_array
    end

    return Dict_array
end

end ## End of module

max_stars_repo_name: oliviersimard/Hartree-Fock

max_issues_repo_name: oliviersimard/Hartree-Fock

max_forks_repo_name: oliviersimard/Hartree-Fock


text:
module TestExamples2DCurved

using Test
using Trixi

include("test_trixi.jl")

# pathof(Trixi) returns /path/to/Trixi/src/Trixi.jl, dirname gives the parent directory
EXAMPLES_DIR = joinpath(pathof(Trixi) |> dirname |> dirname, "examples", "2d")

@testset "Curved Mesh" begin
  @trixi_testset "elixir_advection_basic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_basic_curved.jl"),
      l2   = [9.14468177884088e-6],
      linf = [6.437440532947036e-5])
  end

  @trixi_testset "elixir_advection_extended_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended_curved.jl"),
      l2   = [4.842990962468553e-6],
      linf = [3.47372094784415e-5])
  end

  @trixi_testset "elixir_advection_extended_curved.jl with polydeg=4" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_extended_curved.jl"),
      l2   = [2.0610527374594057e-6],
      linf = [7.21793425673134e-6],
      atol = 1e-12, # required to make CI tests pass on macOS
      cells_per_dimension = (16, 23),
      polydeg = 4,
      cfl = 1.4)
  end

  @testset "elixir_advection_rotated.jl" begin
    @trixi_testset "elixir_advection_rotated.jl with α = 0.0" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_rotated.jl"),
        l2   = [7.013143474176369e-6],
        linf = [4.906526503622999e-5],
        alpha = 0.0)
    end

    @trixi_testset "elixir_advection_rotated.jl with α = 0.1" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_rotated.jl"),
        l2   = [7.013143474176369e-6],
        linf = [4.906526503622999e-5],
        alpha = 0.1)
    end

    @trixi_testset "elixir_advection_rotated.jl with α = 0.5 * pi" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_rotated.jl"),
        l2   = [7.013143474176369e-6],
        linf = [4.906526503622999e-5],
        alpha = 0.5 * pi)
    end
  end

  @trixi_testset "elixir_advection_parallelogram.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_parallelogram.jl"),
      l2   = [0.0005165995033861579],
      linf = [0.002506176163321605])
  end

  @trixi_testset "elixir_advection_waving_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_waving_flag.jl"),
      l2   = [0.00017007823031108628],
      linf = [0.0015264963372674245])
  end

  @trixi_testset "elixir_advection_free_stream_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_free_stream_curved.jl"),
      l2   = [5.984863383701255e-15],
      linf = [1.8207657603852567e-13])
  end

  @trixi_testset "elixir_advection_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_nonperiodic_curved.jl"),
      l2   = [0.00023766972629056245],
      linf = [0.004142508319267935])
  end

  @trixi_testset "elixir_advection_restart_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_curved.jl"),
      l2   = [6.398955192910044e-6],
      linf = [3.474337336717426e-5])
  end

  @trixi_testset "elixir_advection_restart_curved.jl with waving flag mesh" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_curved.jl"),
      l2   = [0.00017274040834067234],
      linf = [0.0015435741643734513],
      elixir_file="elixir_advection_waving_flag.jl",
      restart_file="restart_000041.h5")
  end

  @trixi_testset "elixir_advection_restart_curved.jl with free stream mesh" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_advection_restart_curved.jl"),
      l2   = [6.639325026542281e-15],
      linf = [1.829647544582258e-13],
      elixir_file="elixir_advection_free_stream_curved.jl",
      restart_file="restart_000068.h5")
  end

  @trixi_testset "elixir_euler_source_terms_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_curved.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5])
  end

  @testset "elixir_euler_source_terms_rotated.jl" begin
    @trixi_testset "elixir_euler_source_terms_rotated.jl with α = 0.0" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5],
      alpha = 0.0)
    end

    @trixi_testset "elixir_euler_source_terms_rotated.jl with α = 0.1" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517816067144339e-7, 1.1545335192659009e-6, 1.3105743360575469e-6, 4.2778880887284874e-6],
      linf = [8.357837601113971e-6, 9.243793987812055e-6, 1.1305611891110345e-5, 4.496194865932779e-5],
      alpha = 0.1)
    end

    @trixi_testset "elixir_euler_source_terms_rotated.jl with α = 0.2 * pi" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517916398386897e-7, 7.775434972297097e-7, 1.5639649498758719e-6, 4.277897693664215e-6],
      linf = [8.357487707666422e-6, 4.287882448716918e-6, 1.4423290043641401e-5, 4.496211536153538e-5],
      alpha = 0.2 * pi)
    end

    @trixi_testset "elixir_euler_source_terms_rotated.jl with α = 0.5 * pi" begin
      @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_rotated.jl"),
      l2   = [8.517808508019351e-7, 1.2350203856098537e-6, 1.2350203856728076e-6, 4.277886946638239e-6],
      linf = [8.357848139128876e-6, 1.0326302096741458e-5, 1.0326302101404394e-5, 4.496194024383726e-5],
      alpha = 0.5 * pi)
    end
  end

  @trixi_testset "elixir_euler_source_terms_parallelogram.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_parallelogram.jl"),
      l2   = [1.0893247337495366e-5, 1.0425688625462763e-5, 1.755222105014883e-5, 5.136512290929069e-5],
      linf = [7.449791413693951e-5, 7.621073881791673e-5, 0.00011093303834863733, 0.00039625209916493986])
  end

  @trixi_testset "elixir_euler_source_terms_waving_flag.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_source_terms_waving_flag.jl"),
      l2   = [3.2280277059090205e-5, 3.481614479752735e-5, 2.8784017747658748e-5, 0.00011549476000734391],
      linf = [0.00025339087459608223, 0.0003425481056145152, 0.0002454647901921625, 0.0012806891514367535])
  end

  @trixi_testset "elixir_euler_free_stream_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_curved.jl"),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 6.502853810985698e-13, 2.7000623958883807e-13],
      atol = 7.0e-13)
  end

  @trixi_testset "elixir_euler_free_stream_curved.jl with FluxRotated(flux_lax_friedrichs)" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_free_stream_curved.jl"),
      surface_flux=FluxRotated(flux_lax_friedrichs),
      l2   = [2.063350241405049e-15, 1.8571016296925367e-14, 3.1769447886391905e-14, 1.4104095258528071e-14],
      linf = [1.9539925233402755e-14, 2.9791447087035294e-13, 6.502853810985698e-13, 2.7000623958883807e-13],
      atol = 7.0e-13)
  end

  @trixi_testset "elixir_euler_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_nonperiodic_curved.jl"),
      l2   = [2.3653424742684444e-6, 2.1388875095440695e-6, 2.1388875095548492e-6, 6.010896863397195e-6],
      linf = [1.4080465931654018e-5, 1.7579850587257084e-5, 1.7579850592586155e-5, 5.956893531156027e-5])
  end

  @trixi_testset "elixir_euler_ec_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_euler_ec_curved.jl"),
      l2   = [0.03774170146315357, 0.028439692043402822, 0.028270724308772303, 0.13784229328899064],
      linf = [0.33571598859050134, 0.27831002252333553, 0.295995432205428, 1.2178776811996832],
      tspan = (0.0, 0.3))
  end

  @trixi_testset "elixir_hypdiff_nonperiodic_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_hypdiff_nonperiodic_curved.jl"),
      l2   = [0.8799744480157664, 0.8535008397034816, 0.7851383019164209],
      linf = [1.0771947577311836, 1.9143913544309838, 2.149549109115789],
      tspan = (0.0, 0.1))
  end

  @trixi_testset "elixir_mhd_ec_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_ec_curved.jl"),
      l2   = [0.049292513442569026, 0.06118472341688633, 0.06100981168439492, 0.03154775908310433,
              0.23198302161785492, 0.0247519941440438, 0.024448195911212046, 0.03543873189291887,
              0.0015617736595970992],
      linf = [0.23046364859986535, 0.30519572609771284, 0.39552198293921503, 0.2277249293363471,
              0.9953367157247124, 0.11969759086572018, 0.11429783654481473, 0.1761519810380463,
              0.037502059739661815],
      tspan = (0.0, 0.3))
  end

  @trixi_testset "elixir_mhd_alfven_wave_curved.jl" begin
    @test_trixi_include(joinpath(EXAMPLES_DIR, "elixir_mhd_alfven_wave_curved.jl"),
      l2   = [0.028937533507168618, 0.0062432054639327645, 0.005633373806188361, 0.007334018943006818,
              0.004877728963620626, 0.007192028498071723, 0.006999219481342092, 0.006729459133884995,
              0.004256125425790215],
      linf = [0.17577416304965765, 0.06225741511291068, 0.03880568731344702, 0.05278192154508075,
              0.03811117752054194, 0.0427542385663775, 0.037517543473435055, 0.047902503040458,
              0.04112917360970982],
      tspan = (0.0, 1.0))
  end
end

end # module


text:
immutable IborIndex{S <: AbstractString, TP <: TenorPeriod, CUR <: AbstractCurrency, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure} <: InterestRateIndex
  familyName::S
  tenor::TP
  fixingDays::Int
  currency::CUR
  fixingCalendar::B
  convention::C
  endOfMonth::Bool
  dc::DC
  ts::T
  pastFixings::Dict{Date, Float64}

  # call{S, I, B, C, DC}(::Type{IborIndex}, familyName::S, tenor::TenorPeriod, fixingDays::I, currency::AbstractCurrency, fixingCalendar::B,
  #                               convention::C, endOfMonth::Bool, dc::DC) =
  #   new{S, I, B, C, DC, TermStructure}(familyName, tenor, fixingDays, currency, fixingCalendar, convention, endOfMonth, dc)
  #
  # call{S, I, B, C, DC, T}(::Type{IborIndex}, familyName::S, tenor::TenorPeriod, fixingDays::I, currency::AbstractCurrency, fixingCalendar::B,
  #                               convention::C, endOfMonth::Bool, dc::DC, ts::T) =
  #   new{S, I, B, C, DC, T}(familyName, tenor, fixingDays, currency, fixingCalendar, convention, endOfMonth, dc, ts)
  # IborIndex(familyName::S, tenor::TenorPeriod, fixingDays::I, currency::AbstractCurrency, fixingCalendar::B, convention::C,
  #                   endOfMonth::Bool, dc::DC, ts::T = NullTermStructure()) = IborIndex(familyName, tenor, fixingDays, currency, fixingCalendar, convention,
  #                                                                                endOfMonth, dc, ts)
end

IborIndex{S <: AbstractString, TP <: TenorPeriod, CUR <: AbstractCurrency, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure}(familyName::S, tenor::TP, fixingDays::Int, currency::CUR,
          fixingCalendar::B, convention::C, endOfMonth::Bool, dc::DC, ts::T = NullTermStructure()) =
          IborIndex{S, TP, CUR, B, C, DC, T}(familyName, tenor, fixingDays, currency, fixingCalendar, convention, endOfMonth, dc, ts, Dict{Date, Float64}())

immutable LiborIndex{S <: AbstractString, TP <: TenorPeriod, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure} <: InterestRateIndex
  familyName::S
  tenor::TP
  fixingDays::Int
  currency::Currency
  fixingCalendar::B
  jointCalendar::JointCalendar
  convention::C
  endOfMonth::Bool
  dc::DC
  ts::T
  pastFixings::Dict{Date, Float64}

  # call{S, I, B, C, DC}(::Type{LiborIndex}, familyName::S, tenor::TenorPeriod, fixingDays::I, currency::Currency, fixingCalendar::B,
  #                               jointCalendar::JointCalendar, convention::C, endOfMonth::Bool, dc::DC) =
  #   new{S, I, B, C, DC, TermStructure}(familyName, tenor, fixingDays, currency, fixingCalendar, jointCalendar, convention, endOfMonth, dc)

end

# catch all constructor #
LiborIndex{S <: AbstractString, TP <: TenorPeriod, B <: BusinessCalendar, C <: BusinessDayConvention, DC <: DayCount, T <: TermStructure}(familyName::S, tenor::TP, fixingDays::Int, currency::Currency, fixingCalendar::B,
                  jointCalendar::JointCalendar, convention::C, endOfMonth::Bool, dc::DC, ts::T = NullTermStructure()) = LiborIndex{S, TP, B, C, DC, T}(familyName, tenor, fixingDays, currency, fixingCalendar, jointCalendar, convention,
                                                                                endOfMonth, dc, ts, Dict{Date, Float64}())

function LiborIndex(familyName::AbstractString, tenor::TenorPeriod, fixingDays::Int, currency::Currency, fixingCalendar::BusinessCalendar, dc::DayCount, yts::YieldTermStructure)
  endOfMonth = libor_eom(tenor.period)
  conv = libor_conv(tenor.period)
  jc = JointCalendar(QuantLib.Time.UKLSECalendar(), fixingCalendar)

  return LiborIndex(familyName, tenor, fixingDays, currency, fixingCalendar, jc, conv, endOfMonth, dc, yts)
end

fixing_date(idx::InterestRateIndex, d::Date) = advance(Dates.Day(-idx.fixingDays), idx.fixingCalendar, d, idx.convention)
maturity_date(idx::IborIndex, d::Date) = advance(idx.tenor.period, idx.fixingCalendar, d, idx.convention)
value_date(idx::InterestRateIndex, d::Date) = advance(Dates.Day(idx.fixingDays), idx.fixingCalendar, d)

function fixing(idx::InterestRateIndex, ts::TermStructure, _fixing_date::Date, forecast_todays_fixing::Bool=true)
  today = settings.evaluation_date
  if _fixing_date > today || (_fixing_date == today && forecast_todays_fixing)
    return forecast_fixing(idx, ts, _fixing_date)
  end

  # this is a path fixing
  pastFix = get(idx.pastFixings, _fixing_date, -1.0)

  if pastFix == -1.0
    # none found, default to forecast
    return forecast_fixing(idx, ts, _fixing_date)
  else
    return pastFix
  end
  # error("Not yet implemented for older dates than eval date")
end

function forecast_fixing(idx::InterestRateIndex, ts::TermStructure, _fixing_date::Date)
  d1 = value_date(idx, _fixing_date)
  d2 = maturity_date(idx, d1)
  t = year_fraction(idx.dc, d1, d2)
  return forecast_fixing(idx, ts, d1, d2, t)
end

function forecast_fixing(idx::InterestRateIndex, ts::TermStructure, d1::Date, d2::Date, t::Float64)
  disc1 = discount(ts, d1)
  disc2 = discount(ts, d2)

  return (disc1 / disc2 - 1.0) / t
end

is_valid_fixing_date(idx::InterestRateIndex, d::Date) = is_business_day(idx.fixingCalendar, d)

function add_fixing!(idx::InterestRateIndex, d::Date, fixingVal::Float64)
  idx.pastFixings[d] = fixingVal

  return idx
end

# Libor methods
function value_date(idx::LiborIndex, d::Date)
  new_d = advance(Dates.Day(idx.fixingDays), idx.fixingCalendar, d, idx.convention)
  return adjust(idx.jointCalendar, idx.convention, new_d)
end

maturity_date(idx::LiborIndex, d::Date) = advance(idx.tenor.period, idx.jointCalendar, d, idx.convention)

# types of indexes
euribor_index(tenor::TenorPeriod) = IborIndex("Euribor", tenor, 2, EURCurrency(), QuantLib.Time.TargetCalendar(), euribor_conv(tenor.period), euribor_eom(tenor.period), QuantLib.Time.Actual360())
euribor_index(tenor::TenorPeriod, ts::TermStructure) = IborIndex("Euribor", tenor, 2, EURCurrency(), QuantLib.Time.TargetCalendar(), euribor_conv(tenor.period), euribor_eom(tenor.period), QuantLib.Time.Actual360(), ts)

function usd_libor_index(tenor::TenorPeriod, yts::YieldTermStructure)
  return LiborIndex("USDLibor", tenor, 2, USDCurrency(), QuantLib.Time.USSettlementCalendar(), QuantLib.Time.Actual360(), yts)
end

euribor_conv(::Union{Base.Dates.Day, Base.Dates.Week}) = QuantLib.Time.Following()
euribor_conv(::Union{Base.Dates.Month, Base.Dates.Year}) = QuantLib.Time.ModifiedFollowing()

euribor_eom(::Union{Base.Dates.Day, Base.Dates.Week}) = false
euribor_eom(::Union{Base.Dates.Month, Base.Dates.Year}) = true

libor_conv(::Union{Base.Dates.Day, Base.Dates.Week}) = QuantLib.Time.Following()
libor_conv(::Union{Base.Dates.Month, Base.Dates.Year}) = QuantLib.Time.ModifiedFollowing()

libor_eom(::Union{Base.Dates.Day, Base.Dates.Week}) = false
libor_eom(::Union{Base.Dates.Month, Base.Dates.Year}) = true

# clone methods #
clone(idx::IborIndex, ts::TermStructure = idx.ts) = IborIndex(idx.familyName, idx.tenor, idx.fixingDays, idx.currency, idx.fixingCalendar, idx.convention, idx.endOfMonth, idx.dc, ts)
clone(idx::LiborIndex, ts::TermStructure = idx.ts) = LiborIndex(idx.familyName, idx.tenor, idx.fixingDays, idx.currency, idx.fixingCalendar, idx.jointCalendar, idx.convention, idx.endOfMonth, idx.dc, ts)

update_termstructure(idx::InterestRateIndex, ts::TermStructure) = clone(idx, ts)


text:
using ExcelReaders
using Dates
using PyCall
using DataValues
using Test

@testset "ExcelReaders" begin

# TODO Throw julia specific exceptions for these errors
@test_throws PyCall.PyError openxl("FileThatDoesNotExist.xlsx")
@test_throws PyCall.PyError openxl("runtests.jl")

filename = normpath(@__DIR__, "TestData.xlsx")
file = openxl(filename)
@test file.filename == "TestData.xlsx"

buffer = IOBuffer()
show(buffer, file)
@test String(take!(buffer)) == "ExcelFile <TestData.xlsx>"

for (k,v) in Dict(0=>"#NULL!",7=>"#DIV/0!",23 => "#REF!",42=>"#N/A",29=>"#NAME?",36=>"#NUM!",15=>"#VALUE!")
    errorcell = ExcelErrorCell(k)
    buffer = IOBuffer()
    show(buffer, errorcell)
    @test String(take!(buffer)) == v
end

# Read into DataValueArray
for f in [file, filename]
    @test_throws ErrorException readxl(f, "Sheet1!C4:G3")
    @test_throws ErrorException readxl(f, "Sheet1!G2:B5")
    @test_throws ErrorException readxl(f, "Sheet1!G5:B2")

    data = readxl(f, "Sheet1!C3:N7")
    @test size(data) == (5,12)
    @test data[4,1] == 2.0
    @test data[2,2] == "A"
    @test data[2,3] == true
    @test DataValues.isna(data[4,5])
    @test data[2,9] == Date(2015,3,3)
    @test data[3,9] == DateTime(2015,2,4,10,14)
    @test data[4,9] == DateTime(1988,4,9,0,0)
    @test data[5,9] == Time(15,2,0)
    @test data[3,10] == DateTime(1950,8,9,18,40)
    @test DataValues.isna(data[5,10])
    @test isa(data[2,11], ExcelErrorCell)
    @test isa(data[3,11], ExcelErrorCell)
    @test isa(data[4,12], ExcelErrorCell)
    @test DataValues.isna(data[5,12])

    # Test readxlsheet function
    @test_throws ErrorException readxlsheet(f, "Empty Sheet")
    for sheetinfo=["Second Sheet", 2]
        @test_throws ErrorException readxlsheet(f, sheetinfo, skipstartrows=-1)
        @test_throws ErrorException readxlsheet(f, sheetinfo, skipstartrows=:nonsense)

        @test_throws ErrorException readxlsheet(f, sheetinfo, skipstartcols=-1)
        @test_throws ErrorException readxlsheet(f, sheetinfo, skipstartcols=:nonsense)

        @test_throws ErrorException readxlsheet(f, sheetinfo, nrows=-1)
        @test_throws ErrorException readxlsheet(f, sheetinfo, nrows=:nonsense)

        @test_throws ErrorException readxlsheet(f, sheetinfo, ncols=-1)
        @test_throws ErrorException readxlsheet(f, sheetinfo, ncols=:nonsense)

        data = readxlsheet(f, sheetinfo)
        @test size(data) == (6, 6)
        @test data[2,1] == 1.
        @test data[5,2] == "CCC"
        @test data[3,3] == false
        @test data[6,6] == Time(15,2,00)
        @test DataValues.isna(data[4,3])
        @test DataValues.isna(data[4,6])

        data = readxlsheet(f, sheetinfo, skipstartrows=:blanks, skipstartcols=:blanks)
        @test size(data) == (6, 6)
        @test data[2,1] == 1.
        @test data[5,2] == "CCC"
        @test data[3,3] == false
        @test data[6,6] == Time(15,2,00)
        @test DataValues.isna(data[4,3])
        @test DataValues.isna(data[4,6])

        data = readxlsheet(f, sheetinfo, skipstartrows=0, skipstartcols=0)
        @test size(data) == (6+7, 6+3)
        @test data[2+7,1+3] == 1.
        @test data[5+7,2+3] == "CCC"
        @test data[3+7,3+3] == false
        @test data[6+7,6+3] == Time(15,2,00)
        @test DataValues.isna(data[4+7,3+3])
        @test DataValues.isna(data[4+7,6+3])

        data = readxlsheet(f, sheetinfo, skipstartrows=0, )
        @test size(data) == (6+7, 6)
        @test data[2+7,1] == 1.
        @test data[5+7,2] == "CCC"
        @test data[3+7,3] == false
        @test data[6+7,6] == Time(15,2,00)
        @test DataValues.isna(data[4+7,3])
        @test DataValues.isna(data[4+7,6])

        data = readxlsheet(f, sheetinfo, skipstartcols=0)
        @test size(data) == (6, 6+3)
        @test data[2,1+3] == 1.
        @test data[5,2+3] == "CCC"
        @test data[3,3+3] == false
        @test data[6,6+3] == Time(15,2,00)
        @test DataValues.isna(data[4,3+3])
        @test DataValues.isna(data[4,6+3])

        data = readxlsheet(f, sheetinfo, skipstartrows=1, skipstartcols=1, nrows=11, ncols=7)
        @test size(data) == (11, 7)
        @test data[2+6,1+2] == 1.
        @test data[5+6,2+2] == "CCC"
        @test data[3+6,3+2] == false
        @test_throws BoundsError data[6+6,6+2] == Time(15,2,00)
        @test DataValues.isna(data[4+6,2+2])
    end
end

end


text:
module TestRanges

using Test
using RLEVectors

@testset begin

    # disjoin_length
    @test disjoin_length([], []) == 0
    @test disjoin_length([4, 5], []) == 2
    @test disjoin_length([], [9, 10, 12]) == 3
    @test disjoin_length([1, 2, 3], [4, 5, 6]) == 6
    @test disjoin_length([1, 2, 3], [3, 5, 6]) == 5
    @test disjoin_length([1, 3], [2, 5, 6]) == 5
    @test disjoin_length(RLEVector([1, 2], [1, 3]), RLEVector([3, 2, 1], [2, 5, 6])) == 5

    # disjoin
    @test disjoin([1, 3, 8, 10], [3, 4, 9, 10]) == [1, 3, 4, 8, 9, 10]
    @test disjoin([3, 4, 9, 10], [1, 3, 8, 10]) == [1, 3, 4, 8, 9, 10]

    @test disjoin(RLEVector([1, 2, 3, 4], [1, 3, 8, 10]), RLEVector([4, 5, 6, 7], [3, 4, 9, 10])) ==
          ([1, 3, 4, 8, 9, 10], [1, 2, 3, 3, 4, 4], [4, 4, 5, 6, 6, 7])
    @test disjoin(RLEVector([4, 5, 6, 7], [3, 4, 9, 10]), RLEVector([1, 2, 3, 4], [1, 3, 8, 10])) ==
          ([1, 3, 4, 8, 9, 10], [4, 4, 5, 6, 6, 7], [1, 2, 3, 3, 4, 4])

    # rangeMeans

    x = RLEVector(collect(2:2:20))
    y = RLEVector(collect(1:1:10))

    ranges = [1:1, 1:2, 2:3, 4:8, 1:10, 9:10, 10:10]
    @test [mean(x[r]) for r in ranges] == rangeMeans(ranges, x)

end # testset

end # module