text:
include("../examples/lorenz.jl")
using Test
#using PyPlot
function test_attractor()
	s = [10., 28., 8.0/3]
	u = rand(3)
	n_runup = 20000
	for i = 1:n_runup
		u = next(u,s)
	end
	@show u
	n = 500000
	u_trj = step(u, s, n)
	u_trj = u_trj[:,1:10:end]
	x, y, z = view(u_trj,1,:), view(u_trj,2,:), view(u_trj,3,:)
	fig, ax = subplots(1,1)
	ax.plot(x,y,".",ms=2)
	ax.set_xlabel("x",fontsize=36)
	ax.set_ylabel("y",fontsize=36)
	ax.tick_params(axis="both",labelsize=36)
	fig, ax = subplots(1,1)
	ax.plot(x,z,".",ms=2)
	ax.set_xlabel("x",fontsize=36)
	ax.set_ylabel("z",fontsize=36)
	ax.tick_params(axis="both",labelsize=36)
end
function test_dstep()
	eps = 1.e-8
	s = [10., 28., 8.0/3]
	u = rand(3)
	u1px = u .+ eps*[1,0,0]
	u1mx = u .- eps*[1,0,0]
	u1py = u .+ eps*[0,1,0]
	u1my = u .- eps*[0,1,0]
	u1pz = u .+ eps*[0,0,1]
	u1mz = u .- eps*[0,0,1]

	du_dx_fd = (next(u1px,s) - next(u1mx,s))/(2*eps)
	du_dy_fd = (next(u1py,s) - next(u1my,s))/(2*eps)
	du_dz_fd = (next(u1pz,s) - next(u1mz,s))/(2*eps)

	du = dstep(u, s)
	du_dx = du[:,1]
	du_dy = du[:,2]
	du_dz = du[:,3]

	@show abs.(du_dx_fd .- du_dx)
	@show abs.(du_dy_fd .- du_dy)
	@show abs.(du_dz_fd .- du_dz)
	@test du_dx_fd ≈ du_dx atol=1.e-8
	@test du_dy_fd ≈ du_dy atol=1.e-8
	@test du_dz_fd ≈ du_dz atol=1.e-8
end
function test_dstep_arr()
	eps = 1.e-8
	s = [10., 28., 8.0/3]
	u = rand(3,1)
	u1px = u[:,1] .+ eps*[1,0,0]
	u1mx = u[:,1] .- eps*[1,0,0]
	u1py = u[:,1] .+ eps*[0,1,0]
	u1my = u[:,1] .- eps*[0,1,0]
	u1pz = u[:,1] .+ eps*[0,0,1]
	u1mz = u[:,1] .- eps*[0,0,1]

	du_dx_fd = (next(u1px,s) - next(u1mx,s))/(2*eps)
	du_dy_fd = (next(u1py,s) - next(u1my,s))/(2*eps)
	du_dz_fd = (next(u1pz,s) - next(u1mz,s))/(2*eps)

	du = dstep(u, s)
	@show du
	du_dx = du[:,1,1]
	du_dy = du[:,2,1]
	du_dz = du[:,3,1]

	@test du_dx_fd ≈ du_dx atol=1.e-8
	@test du_dy_fd ≈ du_dy atol=1.e-8
	@test du_dz_fd ≈ du_dz atol=1.e-8
end
function test_d2step()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3)
	u1px = u .+ eps*[1,0,0]
	u1mx = u .- eps*[1,0,0]
	u1py = u .+ eps*[0,1,0]
	u1my = u .- eps*[0,1,0]
	u1pz = u .+ eps*[0,0,1]
	u1mz = u .- eps*[0,0,1]

	ddu_dx_fd = (dstep(u1px,s) - dstep(u1mx,s))/(2*eps)
	ddu_dy_fd = (dstep(u1py,s) - dstep(u1my,s))/(2*eps)
	ddu_dz_fd = (dstep(u1pz,s) - dstep(u1mz,s))/(2*eps)

	ddu = d2step(u, s)
	@show ddu[:,:,1] .- ddu_dx_fd
	@show ddu[:,:,2] .- ddu_dy_fd
	@show ddu[:,:,3] .- ddu_dz_fd

	@test ddu[:,:,1] ≈ ddu_dx_fd atol=1.e-8
	@test ddu[:,:,2] ≈ ddu_dy_fd atol=1.e-8
	@test ddu[:,:,3] ≈ ddu_dz_fd atol=1.e-8

end
function test_d2step_arr()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3,1)
	u1px = u[:,1] .+ eps*[1,0,0]
	u1mx = u[:,1] .- eps*[1,0,0]
	u1py = u[:,1] .+ eps*[0,1,0]
	u1my = u[:,1] .- eps*[0,1,0]
	u1pz = u[:,1] .+ eps*[0,0,1]
	u1mz = u[:,1] .- eps*[0,0,1]

	ddu_dx_fd = (dstep(u1px,s) - dstep(u1mx,s))/(2*eps)
	ddu_dy_fd = (dstep(u1py,s) - dstep(u1my,s))/(2*eps)
	ddu_dz_fd = (dstep(u1pz,s) - dstep(u1mz,s))/(2*eps)

	ddu = d2step(u, s)
	@show ddu[:,:,1,1] .- ddu_dx_fd
	@show ddu[:,:,2,1] .- ddu_dy_fd
	@show ddu[:,:,3,1] .- ddu_dz_fd

	@test ddu[:,:,1,1] ≈ ddu_dx_fd atol=1.e-8
	@test ddu[:,:,2,1] ≈ ddu_dy_fd atol=1.e-8
	@test ddu[:,:,3,1] ≈ ddu_dz_fd atol=1.e-8

end
function test_pert()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3)
	sp = s .+ [0., eps, 0.]
	sm = s .- [0., eps, 0.]

	pp_fd = (next(u, sp) - next(u, sm))/(2*eps)
	pp = pert(u, s)

	@show pp_fd, pp
	@test pp_fd ≈ pp atol=1.e-8
end
function test_pert_arr()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3)
	sp = s .+ [0., eps, 0.]
	sm = s .- [0., eps, 0.]

	pp_fd = (next(u,sp) - next(u, sm))/(2*eps)
	u = reshape(u, 3, 1)
	pp = pert(u, s)[:,1]

	@show pp_fd, pp
	@test pp_fd ≈ pp atol=1.e-8
end
function test_dpert_next()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3)
	u1px = u .+ eps*[1,0,0]
	u1mx = u .- eps*[1,0,0]
	u1py = u .+ eps*[0,1,0]
	u1my = u .- eps*[0,1,0]
	u1pz = u .+ eps*[0,0,1]
	u1mz = u .- eps*[0,0,1]

	dpert_dx_fd = (pert(u1px,s) - pert(u1mx,s))/(2*eps)
	dpert_dy_fd = (pert(u1py,s) - pert(u1my,s))/(2*eps)
	dpert_dz_fd = (pert(u1pz,s) - pert(u1mz,s))/(2*eps)
	dpp = dpert_next(u, s)

	@show dpp[:,1] .- dpert_dx_fd
	@show dpp[:,2] .- dpert_dy_fd
	@show dpp[:,3] .- dpert_dz_fd

	@test dpp[:,1] ≈ dpert_dx_fd atol=1.e-8
	@test dpp[:,2] ≈ dpert_dy_fd atol=1.e-8
	@test dpp[:,3] ≈ dpert_dz_fd atol=1.e-8
end
function test_dpert_next_arr()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3,1)
	u1px = u[:,1] .+ eps*[1,0,0]
	u1mx = u[:,1] .- eps*[1,0,0]
	u1py = u[:,1] .+ eps*[0,1,0]
	u1my = u[:,1] .- eps*[0,1,0]
	u1pz = u[:,1] .+ eps*[0,0,1]
	u1mz = u[:,1] .- eps*[0,0,1]

	dpert_dx_fd = (pert(u1px,s) - pert(u1mx,s))/(2*eps)
	dpert_dy_fd = (pert(u1py,s) - pert(u1my,s))/(2*eps)
	dpert_dz_fd = (pert(u1pz,s) - pert(u1mz,s))/(2*eps)
	dpp = dpert_next(u, s)

	@show dpp[:,1,1] .- dpert_dx_fd
	@show dpp[:,2,1] .- dpert_dy_fd
	@show dpp[:,3,1] .- dpert_dz_fd

	@test dpp[:,1,1] ≈ dpert_dx_fd atol=1.e-8
	@test dpp[:,2,1] ≈ dpert_dy_fd atol=1.e-8
	@test dpp[:,3,1] ≈ dpert_dz_fd atol=1.e-8
end
function test_dmag_flow()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3)
	dmag_dx = (mag_flow(u + eps*[1.0,0,0],s) -
			   mag_flow(u - eps*[1.0,0,0],s))/(2*eps)
	dmag_dy = (mag_flow(u + eps*[0,1.0,0],s) -
			   mag_flow(u - eps*[0,1.0,0],s))/(2*eps)
	dmag_dz = (mag_flow(u + eps*[0,0,1.0],s) -
			   mag_flow(u - eps*[0,0,1.0],s))/(2*eps)

	dmag = dmag_flow(u, s)
	@test dmag[1] ≈ dmag_dx atol=1.e-8
	@test dmag[2] ≈ dmag_dy atol=1.e-8
	@test dmag[3] ≈ dmag_dz atol=1.e-8
end
function test_dmag_flow_arr()
	eps = 1.e-6
	s = [10., 28., 8.0/3]
	u = rand(3,1)
	dmag_dx = (mag_flow(u .+ eps*reshape([1.0,0,0],3,1),s) -
			   mag_flow(u .- eps*reshape([1.0,0,0],3,1),s))/(2*eps)
	dmag_dy = (mag_flow(u .+ eps*reshape([0,1.0,0],3,1),s) -
			   mag_flow(u .- eps*reshape([0,1.0,0],3,1),s))/(2*eps)
	dmag_dz = (mag_flow(u .+ eps*reshape([0,0,1.0],3,1),s) -
			   mag_flow(u .- eps*reshape([0,0,1.0],3,1),s))/(2*eps)

	dmag = dmag_flow(u, s)
	@test dmag[1,:] ≈ dmag_dx atol=1.e-8
	@test dmag[2,:] ≈ dmag_dy atol=1.e-8
	@test dmag[3,:] ≈ dmag_dz atol=1.e-8
end
function test_dunit_flow()
	eps = 1.e-6
	s = [10., 28., 8/3]
	u = rand(3)
	upx = u .+ eps*[1.0,0,0]
	umx = u .- eps*[1.0,0,0]
	upy = u .+ eps*[0,1.0,0]
	umy = u .- eps*[0,1.0,0]
	upz = u .- eps*[0,0,1.0]
	umz = u .- eps*[0,0,1.0]
	dfdx = (flow(upx,s)/mag_flow(upx,s) .-
			flow(umx,s)/mag_flow(umx,s))/(2*eps)
	dfdy = (flow(upy,s)/mag_flow(upy,s) .-
			flow(umy,s)/mag_flow(umy,s))/(2*eps)
	dfdz = (flow(upz,s)/mag_flow(upz,s) .-
			flow(umz,s)/mag_flow(umz,s))/(2*eps)
	df = dunit_flow(u, s)
	@show dfdx, df[:,1]

	@test df[:,1] ≈ dfdx atol=1.e-8
	@test df[:,2] ≈ dfdy atol=1.e-8
	@test df[:,3] ≈ dfdz atol=1.e-8

end

max_stars_repo_name: nishaChandramoorthy/linearResponse

max_issues_repo_name: nishaChandramoorthy/linearResponse

max_forks_repo_name: nishaChandramoorthy/linearResponse


text:
module BenchCollectors

using BangBang: Empty, collector, finish!, push!!
using BenchmarkTools: @benchmarkable, BenchmarkGroup
using SplittablesBase: amount

function map_manual(f, src)
    dest = similar(src)
    for i in eachindex(src, dest)
        @inbounds dest[i] = f(src[i])
    end
    return dest
end

function map_collector(f, src, unsafe = Val(true))
    dest = similar(src)
    c = collector(dest, unsafe)
    for x in src
        c = push!!(c, f(x))
    end
    return finish!(c)
end

function map_vector(f, src)
    dest = eltype(src)[]
    for x in src
        dest = push!!(dest, f(x))
    end
    return dest
end

@inline function _foldl(op, acc::T, itr) where {T}
    y = iterate(itr)
    y === nothing && return acc
    (x, state) = y
    result = op(acc, x)
    if result isa T
        return _foldl(op, result, itr, state)
    else
        return _foldl(op, result, itr, state)
    end
end

@inline function _foldl(op, acc::T, itr, state) where {T}
    while true
        y = iterate(itr, state)
        y === nothing && return acc
        (x, state) = y
        result = op(acc, x)
        if result isa T
            acc = result
        else
            return _foldl(op, result, itr, state)
        end
    end
end

collect_collector(src) = finish!(_foldl(push!!, collector(), src))

collect_vector(src) = _foldl(push!!, Empty(Vector), src)

collect_collector_nonexpanding(src) =
    finish!(_foldl(push!!, collector(resize!(Union{}[], amount(src)), Val(true)), src))

double(x) = 2 * x

suite = BenchmarkGroup()

let s0 = suite[:comp=>"map"] = BenchmarkGroup()
    xs = 1:10_000
    @assert map_manual(double, xs) == map_collector(double, xs) == map_vector(double, xs)
    s0[:impl=>"man"] = @benchmarkable(map_manual(double, $xs))
    s0[:impl=>"coll"] = @benchmarkable(map_collector(double, $xs))
    s0[:impl=>"safe_coll"] = @benchmarkable(map_collector(double, $xs, Val(false)))
    s0[:impl=>"vec"] = @benchmarkable(map_vector(double, $xs))
end

let s0 = suite[:comp=>"filter"] = BenchmarkGroup()
    xs = (x for x in 1:10_000 if isodd(x))
    @assert collect(xs) ==
            collect_collector(xs) ==
            collect_vector(xs) ==
            collect_collector_nonexpanding(xs)
    s0[:impl=>"base"] = @benchmarkable(collect($xs))
    s0[:impl=>"coll"] = @benchmarkable(collect_collector($xs))
    s0[:impl=>"vec"] = @benchmarkable(collect_vector($xs))
    s0[:impl=>"nonexpanding"] = @benchmarkable(collect_collector_nonexpanding($xs))
end

end  # module
BenchCollectors.suite

ext: jl


text:
##############################################################################
#
# test_3d_vortex_particles.jl
#
# Check 3D vortex particles behave as expected.
#
# Copyright 2020 HJA Bird
#
# Permission is hereby granted, free of charge, to any person obtaining a copy
# of this software and associated documentation files (the "Software"), to
# deal in the Software without restriction, including without limitation the
# rights to use, copy, modify, merge, publish, distribute, sublicense, and/or
# sell copies of the Software, and to permit persons to whom the Software is
# furnished to do so, subject to the following conditions:
#
# The above copyright notice and this permission notice shall be included in
# all copies or substantial portions of the Software.
#
# THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
# IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
# FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
# AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
# LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
# FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS
# IN THE SOFTWARE.
##############################################################################

function two_particles_basic_case_3d()
	@testset "3D 2 particle interaction" begin
		pl0 = [0,0,0]
		plz = [0,0,1]
		ply = [0,1,0]
		plx = [1,0,0]
		pvx = [1,0,0]
		pvy = [0,1,0]
		pvz = [0,0,1]

		# Singular kernel ---------------------------------------------------------------------------------------
		kernel = singular_regularisation()
		@test isapprox([0,-1,0]./(4*pi), particle_induced_velocity(pl0, pvx, plz, kernel, 1.); rtol=1e-5)
		@test isapprox([-1,0,0]./(4*pi), particle_induced_velocity(plz, pvy, pl0, kernel, 1.); rtol=1e-5)
		@test isapprox([0,0,0]./(4*pi), particle_induced_velocity(pl0, pvx, pl0, kernel, 1.); rtol=1e-5)
		@test isapprox([0,0,0]./(4*pi), particle_induced_velocity(plz, pvy, plz, kernel, 1.); rtol=1e-5)
		@test isapprox([0,0,0]./(4*pi), particle_induced_velocity(pl0, pvz, plz, kernel, 1.); rtol=1e-5)
		@test isapprox([0,0,0]./(4*pi), particle_induced_velocity(plz, pvz, pl0, kernel, 1.); rtol=1e-5)
		@test isapprox([0,-0.25,0]./(4*pi), particle_induced_velocity(pl0, pvx, 2*plz, kernel, 1.); rtol=1e-5)
		@test isapprox([-0.25,0,0]./(4*pi), particle_induced_velocity(2*plz, pvy, pl0, kernel, 1.); rtol=1e-5)
		@test isapprox([0,-2,0]./(4*pi), particle_induced_velocity(pl0, 2*pvx, plz, kernel, 1.); rtol=1e-5)
		@test isapprox([-2,0,0]./(4*pi), particle_induced_velocity(plz, 2*pvy, pl0, kernel, 1.); rtol=1e-5)
		@test isapprox([0,0,1]./(4*pi), particle_induced_velocity(pl0, pvx, ply, kernel, 1.); rtol=1e-5)
		@test isapprox([0,0,-1]./(4*pi), particle_induced_velocity(ply, pvx, pl0, kernel, 1.); rtol=1e-5)
		# Check these don't fail.
		kernel = gaussian_regularisation()
		kernel = planetary_regularisation()
		kernel = winckelmans_regularisation()
		# Check visc dvort.
		kernel = gaussian_regularisation()
		gaussr1 = sqrt(2/pi) * exp(-1/2)
		gaussr2 = sqrt(2/pi) * exp(-4/2)
		@test isapprox((pvz-pvy)*2*gaussr1,
			particle_visc_induced_dvort(ply, pvz, 1, pl0, pvy, 1, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvx-pvy)*2*gaussr1,
			particle_visc_induced_dvort(plx, pvx, 1, pl0, pvy, 1, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvz-pvz)*2*gaussr1,
			particle_visc_induced_dvort(ply, pvz, 1, pl0, pvz, 1, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvz-pvy)*2*gaussr2,
			particle_visc_induced_dvort(2*ply, pvz, 1, pl0, pvy, 1, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvx-pvy)*2*gaussr2,
			particle_visc_induced_dvort(2*plx, pvx, 1, pl0, pvy, 1, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvz-pvz)*2*gaussr2,
			particle_visc_induced_dvort(2*ply, pvz, 1, pl0, pvz, 1, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvz-pvy)*gaussr1,
			particle_visc_induced_dvort(ply, pvz, 0.5, pl0, pvy, 0.5, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvz-pvy)*gaussr1,
			particle_visc_induced_dvort(ply, 0.5*pvz, 0.5, pl0, pvy, 1.0, kernel, 1., 1.); rtol=1e-5)
		@test isapprox((pvz-pvy)*4*gaussr1,
			particle_visc_induced_dvort(ply, pvz, 1, pl0, pvy, 1, kernel, 1., 2.); rtol=1e-5)
		@test isapprox((pvz-pvy)*8*gaussr2,
			particle_visc_induced_dvort(ply, pvz, 1, pl0, pvy, 1, kernel, 0.5, 1.); rtol=1e-5)
	end
end

function test_correct_output_matrix_dims_3d_vortex_particles()
	@testset "3D particle correct output dimensions" begin
		n1 = 31
		n2 = 47
		kernel = singular_regularisation()
		@test (3,) == size(particle_induced_velocity(rand(3), rand(3), rand(3), kernel, 1.))
		@test (3,) == size(particle_induced_velocity(rand(n1, 3), rand(n1, 3), rand(3), kernel, 1.))
		@test (n1, 3) == size(particle_induced_velocity(rand(n1, 3), rand(n1, 3), rand(n1, 3), kernel, 1.))
		@test (n2, 3) == size(particle_induced_velocity(rand(n1, 3), rand(n1, 3), rand(n2, 3), kernel, 1.))
		@test (n1, 3) == size(particle_induced_velocity(rand(3), rand(3), rand(n1, 3), kernel, 1.))
		@test (3,) == size(particle_induced_dvort(rand(3), rand(3), rand(3), rand(3), kernel, 1.))
		@test (3,) == size(particle_induced_dvort(rand(n1, 3), rand(n1, 3), rand(3), rand(3), kernel, 1.))
		@test (n1, 3) == size(particle_induced_dvort(rand(n1, 3), rand(n1, 3), rand(n1, 3), rand(n1, 3), kernel, 1.))
		@test (n2, 3) == size(particle_induced_dvort(rand(n1, 3), rand(n1, 3), rand(n2, 3), rand(n2, 3), kernel, 1.))
		@test (n1, 3) == size(particle_induced_dvort(rand(3), rand(3), rand(n1, 3), rand(n1, 3), kernel, 1.))
		kernel = gaussian_regularisation()
		@test (3,) == size(particle_visc_induced_dvort(rand(3), rand(3), rand(), rand(3), rand(3), rand(), kernel, 1., 1))
		@test (3,) == size(particle_visc_induced_dvort(rand(n1, 3), rand(n1, 3), rand(n1), rand(3), rand(3), rand(), kernel, 1., 1))
		@test (n1, 3) == size(particle_visc_induced_dvort(rand(n1, 3), rand(n1, 3), rand(n1), rand(n1, 3), rand(n1, 3), rand(n1), kernel, 1., 1))
		@test (n2, 3) == size(particle_visc_induced_dvort(rand(n1, 3), rand(n1, 3), rand(n1), rand(n2, 3), rand(n2, 3), rand(n2), kernel, 1., 1))
		@test (n1, 3) == size(particle_visc_induced_dvort(rand(3), rand(3), rand(), rand(n1, 3), rand(n1, 3), rand(n1), kernel, 1., 1))
	end
end

test_correct_output_matrix_dims_3d_vortex_particles()
two_particles_basic_case_3d()


text:
using Revise
using Test, ECCOtour
using MITgcmTools, MeshArrays, Statistics, Dierckx
using GoogleDrive

@testset "ECCOtour.jl" begin

    pth = MeshArrays.GRID_LLC90
    γ = GridSpec("LatLonCap",pth)
    Γ = GridLoad(γ;option="full")
    nf = length(γ.fSize)
    z = depthlevels(γ)
    pstdz = pressurelevels(z)

    projectdir = dirname(Base.active_project())
    datadir = joinpath(projectdir,"data")
    srcdir = joinpath(projectdir,"src")

    !isdir(datadir) && mkdir(datadir)

    # download sample data sets

    # state 3d
    url = "https://docs.google.com/uc?export=download&id=1Sst5Y9AUbef1-Vk2ocBgOOiI2kudYRPx"
    filegz = google_download(url,datadir)
    cd(datadir)
    run(`tar xvzf $filegz`)

    # transport 3d # too large doesn't work due to virus scan
    #url = "https://docs.google.com/uc?export=download&id=1KKk8d_1nQFbM9xQjTelCmWTMfK3SA7U5"
    #filegz = google_download(url,datadir)

    cd(datadir)
    # workaround: use a shell script
    run(`sh $srcdir/download_google_drive.sh`)
    run(`tar xvzf trsp_3d_set1.0000000732.tar.gz`)

    ## specific for state
    # the state_3d monthly-average diagnostic output
    stateroot = "state_3d_set1" # 1: θ, 2: S

    # first filter for state_3d_set1
    filelist = searchdir(datadir,stateroot)
    # second filter for "data"
    statefile  = filter(x -> occursin("data",x),filelist)

    # Read from filelist, map to sigma-1, write to file
    fileroots = Vector{String}()
    fileroot = rstrip(statefile[1],['.','d','a','t','a'])
    push!(fileroots,fileroot)

    ## specific for transport
    transportroot = "trsp_3d_set1" # 1: θ, 2: S

    # first filter for state_3d_set1
    filelist = searchdir(datadir,transportroot)
    # second filter for "data"
    transportfile  = filter(x -> occursin("data",x),filelist)

    # Read from filelist, map to sigma-1, write to file
    fileroot = rstrip(transportfile[1],['.','d','a','t','a'])
    push!(fileroots,fileroot)

    @testset "IsopycnalSurfaces" begin

        # DEFINE THE LIST OF SIGMA1 VALUES.
        sig1grid = sigma1grid()
        eos_mitgcm = "JMD95"

        @testset "spline_interpolation" begin
            varsσ = mdsio2sigma1(datadir,datadir,fileroots,γ,pstdz,sig1grid,splorder=3,eos=eos_mitgcm)
            for ss in eachindex(sig1grid)
                @test maximum(MeshArrays.mask(varsσ["SALT"][:,ss],-Inf)) < 45.0
                @test minimum(MeshArrays.mask(varsσ["SALT"][:,ss],Inf)) ≥ 0.0
                @test maximum(MeshArrays.mask(varsσ["THETA"][:,ss],-Inf)) < 35.0
                @test minimum(MeshArrays.mask(varsσ["THETA"][:,ss],Inf)) ≥ -3.0
                @test maximum(MeshArrays.mask(varsσ["NVELMASS"][:,ss],-Inf)) < 3.0
                @test minimum(MeshArrays.mask(varsσ["NVELMASS"][:,ss],Inf)) ≥ -3.0
            end

        end

        @testset "linear_interpolation" begin

            varsσ = mdsio2sigma1(datadir,datadir,fileroots,γ,pstdz,sig1grid,linearinterp=true,eos=eos_mitgcm)

            for ss in eachindex(sig1grid)
                @test maximum(MeshArrays.mask(varsσ["SALT"][:,ss],-Inf)) < 45.0
                @test minimum(MeshArrays.mask(varsσ["SALT"][:,ss],Inf)) ≥ 0.0
                @test maximum(MeshArrays.mask(varsσ["THETA"][:,ss],-Inf)) < 35.0
                @test minimum(MeshArrays.mask(varsσ["THETA"][:,ss],Inf)) ≥ -3.0
                @test maximum(MeshArrays.mask(varsσ["NVELMASS"][:,ss],-Inf)) < 3.0
                @test minimum(MeshArrays.mask(varsσ["NVELMASS"][:,ss],Inf)) ≥ -3.0
            end

        end

        @testset "regularpoles" begin
            ####################################
            # Test the mapping onto the regularpoles grid.
            lat,lon = latlon(γ)

            # Set up Cartesian grid for interpolation.
            # Time for a structure.
            λC,λG,ϕC,ϕG,nx,ny,nyarc,nyantarc,λarc,λantarc =
                factors4regularpoles(γ)

            @testset "regularpoles 3d state" begin

                filein = fileroots[1]
                pathin = datadir

                @time varsregpoles =  mdsio2regularpoles(pathin,filein,γ,nx,ny,nyarc,λarc,nyantarc,λantarc)

                @test maximum(filter(!isnan,varsregpoles["SALT"])) < 50.
                @test minimum(filter(!isnan,varsregpoles["SALT"])) > 0.

            end

            @testset "regularpoles 3d transport" begin

                filein = fileroots[2]
                pathin = datadir

                @time varsregpoles =  mdsio2regularpoles(pathin,filein,γ,nx,ny,nyarc,λarc,nyantarc,λantarc)

                @test maximum(filter(!isnan,varsregpoles["NVELMASS"])) < 2.
                @test minimum(filter(!isnan,varsregpoles["NVELMASS"])) > -2.

            end

            @testset "regularpoles dxc" begin
                # load centered longitude
                # dxc = γ.read(γ.path*"XC.data",MeshArray(γ,Float64))
                vars = Dict("XC" => γ.read(γ.path*"XC.data",MeshArray(γ,Float64)))
                dxc_regpoles = vars2regularpoles(vars,γ,nx,ny,nyarc,λarc,nyantarc,λantarc)

                yy = 100
                for xx = 1:nx
                    println(@test isapprox(dxc_regpoles["XC"][xx,yy],λC[xx], rtol=1e-6))
                end
            end

        end

        @testset "MeshArrays" begin
            ######################################
            # Test the statistics for MeshArrays.
            lat,lon = latlon(γ)
            @test  -360.0 ≤ mean(lat,NaN) ≤ 360.0

            # can it find a NaN?
            latnan = deepcopy(lat)
            latnan[1][30,30] = NaN
            @test sum(nancount(latnan)-nancount(lat)) == 1

            # nancount work with different array size?
            lat2 = MeshArrays.gcmarray(γ,Float64,2)
            lat2[:,1] = latnan;
            lat2[:,2] = latnan;
            @test sum(nancount(lat2)) == 2*sum(nancount(latnan))

            @test maximum(MeshArrays.mask(latnan,-Inf)) ≤ 90.0
            @test minimum(MeshArrays.mask(lat,Inf)) ≥ -90.0

            @test !isnan(std(lat,mean(lat,NaN),NaN))

            # replace NaN with zero
            @test iszero(nancount(MeshArrays.mask(latnan,0.0)))

            # does this run?
            faststats(lat)

        end
    end
end


text:
import Base: show, string

function string(a::RingList)
    if length(a) == 0
        return "[ ]"
    end
    result = "[ "
    arrow = " → "
    for x in a
        result = result * string(x) * arrow
    end
    result = result * string(first(a,true)) * " ]"
    return result
end

function show(io::IO, a::RingList)
    print(io,string(a))
end

# function show(io::IO, a::RingList{T}) where {T}
#     v = Vector(a)
#     result = "RingList{$T}("
#     n = length(a)
#     for i = 1:n
#         result *= "$(v[i])"
#         if i < n
#             result *= ","
#         end
#     end
#     result *= ")"
#     print(io, result)
# end


text:
using DataStructures
"""    dir_base_ext(ljhname::AbstractString)
Given an ljh filename like "somedir/a_chan1.ljh" returns ("somedir,")
"""
function dir_base_ext(ljhname::AbstractString)::Tuple{String,String,String}
    if isdir(ljhname)
        dname = ljhname
        if dname[end] == '/'
            dname = dname[1:end-1]
        end
        bname = last(split(dname,'/'))
        if bname == ""
            bname = dname
        end
        return dname, bname, ".ljh"
    end
    bname,ext = splitext(basename(ljhname))
    ext = isempty(ext) ? ".ljh" : ext
    m = match(r"_chan\d+", bname)
    dname = dirname(ljhname)
    outdirname = dname == "" ? "." : dname # return . for instead of empty string for local dir
    outdirname, String(m == nothing ? bname : bname[1:m.offset-1]), ext
end
"    channel(ljhname::AbstractString)
Return the Channel number determined from the filename of an ljh file, looks for
the _chan1 part. Returns -1 if not found."
function channel(ljhname::AbstractString)::Int
    m = match(r"_chan(\d+)", ljhname)
    m == nothing ? -1 : parse(Int,m.captures[1])
end
"    fnames(ljhname::AbstractString, chans)
Returns a `Vector{String}` of ljh filenames `ljhname`, but with different
channel numbers one for each channel number in `chans`."
function fnames(ljhname::AbstractString, chans)
    dname, bname, ext = dir_base_ext(ljhname)
    String[joinpath(dname, "$(bname)_chan$c$ext") for c in chans]
end
"   fnames(ljhname::AbstractString, c::Int)
Retuns a new ljh filename with channel number `c`."
function fnames(ljhname::AbstractString, c::Int)
    dname, bname, ext = dir_base_ext(ljhname)
    joinpath(dname, "$(bname)_chan$c$ext")
end
"    allchannelnumbers(ljhname::AbstractString)
Returns a sorted `Vector{Int}` containing all the channel numbers of existing
ljh files in the same set as `ljhname`."
function allchannelnumbers(ljhname::AbstractString)
    dname, bname, ext = dir_base_ext(ljhname)
    potential_ljh = filter!(s->(startswith(s,bname)), readdir(dname))
    channels = filter!(x->x>=0,[channel(p) for p in potential_ljh])
    sort!(unique(channels))
end
"    function allchannels(ljhname::AbstractString,maxchannels=typemax(Int))
Returns an `OrderedDict` mapping channel number to filename containing all file
names of existing ljh files in the same set as `ljhname` and all channel numbers.
If there are more than `maxchannels` such files, only the first `maxchannels`
are included."
function allchannels(ljhname::AbstractString,maxchannels=typemax(Int))
    channels = allchannelnumbers(ljhname)
    n = min(length(channels), maxchannels)
    channels = channels[1:n]
    OrderedDict{Int,String}(ch=>fname for (ch,fname) in zip(channels, fnames(ljhname, channels)))
end

"""    outputname(ljhname::AbstractString, annotation::AbstractString, ext=".hdf5")
Returns the filename based on `ljhname` with an _annotation_, and with the desired extention `ext`.
`ext` may be specified with or without a leading period with identical results.
For example `outputname("abc","model")` returns `"abc_model.hdf5"` and
`outputname("abc","model","pdf")` returns `"abc_model.pdf"`."""
function outputname(ljhname::AbstractString, annotation::AbstractString, ext=".hdf5")
    ext = lstrip(ext,'.') #allows .pdf or pdf for 3rd argument
    d,b,e=dir_base_ext(ljhname)
    joinpath(d,b*"_"*annotation*"."*ext)
end

const sentinel_file_path = joinpath(expanduser("~"),".daq","latest_ljh_pulse.cur")
"    matter_writing_status()
returns a tuple (String, Bool) representing (filename, currently_open)"
function matter_writing_status()
    isfile(sentinel_file_path) || error("$(sentinel_file_path) must be a file")
    open(sentinel_file_path,"r") do sentinel_file
    lines = map(chomp,collect(eachline(sentinel_file)))
    if length(lines)>=1
      return lines[1], length(lines)==1
    else
      return "empty sentinel file", false
    end
	# the sentinel file has a second line that says closed when it has closed a file
	# so one line means open, two lines means closed
    end
end

"    write_sentinel_file(filename, writingbool)
MATTER writes it's writing status to a sentinal file, this emulates matters output for testing purposes.
write_sentinel_file(filename, writingbool), if writingbool is true, it writes that the file is still open.
Will create the `.daq` folder for the sentinel file if it does not already exist."
function write_sentinel_file(filename, writingbool)
    dname = dirname(sentinel_file_path)
    !isdir(dname) && mkdir(dname)
    open(sentinel_file_path,"w") do sentinel_file
        println(sentinel_file,filename)
        if !writingbool
            println(sentinel_file, "closed")
        end
    end
end

"    change_writing_status()
Change the sentinel_file writing status, false->true or true->false."
function change_writing_status()
    fname, writingbool = matter_writing_status()
    write_sentinel_file(fname, !writingbool)
end


text:
# Tables shifted by 1

# tables which don't depend on x,q

"""
    tabIln

Table of values of `I_ln=I(m,l,n)` used for computing Legendre integrals.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
const tabIln = zeros(Float64,2*N+1,2*N+1)

"""
    tabIaln

Table of values of `Ia_ln=Ia(m-1;l,m-1)` used for computing Legendre integrals.
See `Aoki et al. (1979)` for a definition (a=7/4).
"""
const tabIaln = zeros(Float64,2*N+1)

"""
    tabJln

Table of values of `J_ln=J(m,l,n)` used for computing Legendre integrals.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
const tabJln = zeros(Float64,2*N+1,2*N+1)

"""
    tabLegendre

Table of values of normalized Legendre associated functions
See `computePlmcostheta` of the `SphericalHarmonics` to see the normalization.
"""
const tabLegendre = zeros(Float64,N+1,nbK)

# tables which depend on x,q

"""
    tabOmega_serial

Table of values of the angular velocity of the system.
"""
tabOmega_serial = zeros(Float64,nbK)

"""
    tabKappaSqOverTwoOmega_serial

Table of values of the ratio of `kappa^2/(2 Omega)`.
"""
tabKappaSqOverTwoOmega_serial = zeros(Float64,nbK)

"""
    tabAln_serial

Table of values of the matrix `Aln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabAln_serial = zeros(Float64,N+1,N+1)

"""
    tabBln_serial

Table of values of the matrix `Bln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabBln_serial = zeros(Float64,N+1,N+1)

"""
    tabCln_serial

Table of values of the matrix `Cln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabCln_serial = zeros(Float64,N+1,N+1)

"""
    tabDln_serial

Table of values of the matrix `Dln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabDln_serial = zeros(Float64,N+1,N+1)

"""
    tabFln_serial

Table of values of the matrix `Fln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabFln_serial = zeros(Float64,N+1,N+1)

"""
    tabGln_serial

Table of values of the matrix `Gln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabGln_serial = zeros(Float64,N+1,N+1)

"""
    tabHln_serial

Table of values of the matrix `Hln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
tabHln_serial = zeros(Float64,N+1,N+1)

###############################################
# Fill Iln table
###############################################

"""
    I_seed!()

Initializes `tabIln[1,1]` with the value `I(m,m,m)`.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
function I_seed!()
    prod = 1/(3/4+m+1)
    for k=1:m
        prod *= 2*(2*k+1)/(3/4+m+1+k)
    end
    tabIln[1,1] = prod
end

"""
    I_firstLine!()

Initializes the first column `tabIln[k,1]` with the values `I(m,l,m)` for k=2,...,2N+1.
Initializes the first lign `tabIln[1,k]` by symmetry.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
function I_firstLine!()
    current = tabIln[1,1]
    for l=m+1:m+2*N
        current *= (l-3/4-m-1)/(l+3/4+m+1)*sqrt(((l+m)*(2*l+1))/((l-m)*(2*l-1)))
        tabIln[l-m+1,1] = current
        tabIln[1,l-m+1] = current #symmetry
   end
end

"""
    I_fill!()

Computes the table `tabIln` by recursion.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
function I_fill!()
    for n=m+1:m+N
        for l=n:m+2*N-(n-m)
            Ip = (sqrt(((l+m+1)*(l-m+1))/((2*l+1)*(2*l+3)))*tabIln[l-m+1+1,n-m+1-1]
                 + sqrt(((l+m)*(l-m))/((2*l+1)*(2*l-1))) *tabIln[l-m+1-1,n-m+1-1])
            tabIln[l-m+1,n-m+1] = sqrt(((2*n+1)*(2*n-1))/((n+m)*(n-m)))*Ip
            if (n>=m+2)
                tabIln[l-m+1,n-m+1] -= sqrt(((n+m-1)*(n-m-1)*(2*n+1))/((n+m)*(n-m)*(2*n-3)))*tabIln[l-m+1,n-m+1-2]
            end
            tabIln[n-m+1,l-m+1] = tabIln[l-m+1,n-m+1] #symmetry
        end
    end
end

"""
    tabIln!()

Computes the whole table `tabIln`.
See `Aoki et al. (1979)` (a=3/4).
"""
function tabIln!()
    I_seed!()
    I_firstLine!()
    I_fill!()
end

"""
    tabIln_clear!()

Clears the whole table `tabIln`.
"""
function tabIln_clear!()
    for i=1:2*N+1
        for j=1:2*N+1
            tabIln[i,j] = 0.0
        end
    end
end

###############################################
# Fill Ialm table : Ia(m-1;l,m-1) with l=m-1,...,m+2N-1
###############################################

"""
    Ia_seed!()

Initializes `tabIln[1,1]` with the value `Ia(m-1,m-1,m-1)`.
See `Aoki et al. (1979)` for a definition (a=7/4).
"""
function Ia_seed!()
    prod = 1/(7/4+m-1+1)
    for k=1:m-1
        prod *= 2*(2*k+1)/(7/4+m-1+1+k)
    end
    tabIaln[1] = prod
end

"""
    Ia_firstLine!()

Initializes the table `tabIaln[k]` with the values `Ia(m-1,l,m-1)` for k=2,...,2N+1.
See `Aoki et al. (1979)` for a definition (a=7/4).
"""
function Ia_firstLine!()
    current = tabIaln[1]
    for l=m:m+2*N-1
        current *= (l-7/4-m+1-1)/(l+7/4+m-1+1)*sqrt(((l+m-1)*(2*l+1))/((l-m+1)*(2*l-1)))
        tabIaln[l-m+2] = current
   end
end

"""
    tabIaln!()

Computes the whole table `tabIaln`.
See `Aoki et al. (1979)` (a=7/4).
"""
function tabIaln!()
    Ia_seed!()
    Ia_firstLine!()
end

"""
    tabIaln_clear!()

Clears the whole table `tabIaln`.
"""
function tabIaln_clear!()
    for i=1:2*N+1
        tabIaln[i] = 0.0
    end
end

###############################################
# Fill Jln table
###############################################

"""
    J_seed!()

Initializes `tabJln[1,1]` with the value `J(m,m,m)`.
Initializes `tabJln[2,1]` with the value `J(m,m+1,m)` and `tabJln[1,2]` by symmetry.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
function J_seed!()
     tabJln[1,1] = tabIln[1,1]*(3/4+2*m+1)/(m)
     tabJln[2,1] = -(7/4)*sqrt(2*m+3)/(3/4+2*m+1)*tabJln[1,1]
     tabJln[1,2] = tabJln[2,1] #symmetry
end

"""
    I_firstLine!()

Initializes the first column `tabJln[k,1]` with the values `J(m,l,m)` for k=3,...,2N+1.
Initializes the first lign `tabJln[1,k]` by symmetry.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
function J_firstLine!()
    for l=m+2:m+2*N
        current = (sqrt(((l-m)*(l-m-1)*(2*l+1))/((l+m)*(l+m-1)*(2*l-3)))*tabJln[l-m+1-2,1]
                  + 4*sqrt(((2*m+1)*(2*l+1)*(2*l-1))/((2*m)*(l+m)*(l+m-1)))*tabIaln[l-m+2-1])
        tabJln[l-m+1,1] = current
        tabJln[1,l-m+1] = current #symmetry
   end
end

"""
    J_fill!()

Computes the table `tabJln` by recursion.
See `Aoki et al. (1979)` for a definition (a=3/4).
"""
function J_fill!()
    for n=m+1:m+N
        for l=n:m+2*N-(n-m)
            Jp = (sqrt(((l+m+1)*(l-m+1))/((2*l+1)*(2*l+3)))*tabJln[l-m+1+1,n-m+1-1]
                 + sqrt(((l+m)*(l-m))/((2*l+1)*(2*l-1))) *tabJln[l-m+1-1,n-m+1-1])
            tabJln[l-m+1,n-m+1] = sqrt(((2*n+1)*(2*n-1))/((n+m)*(n-m)))*Jp
            if (n>=m+2)
                tabJln[l-m+1,n-m+1] -= sqrt(((n+m-1)*(n-m-1)*(2*n+1))/((n+m)*(n-m)*(2*n-3)))*tabJln[l-m+1,n-m+1-2]
            end
            tabJln[n-m+1,l-m+1] = tabJln[l-m+1,n-m+1] #symmetry
        end
    end
end

"""
    tabJln!()

Computes the whole table `tabJln`.
See `Aoki et al. (1979)` (a=3/4).
"""
function tabJln!()
    J_seed!()
    J_firstLine!()
    J_fill!()
end

"""
    tabJln_clear!()

Clears the whole table `tabJln`.
"""
function tabJln_clear!()
    for i=1:2*N+1
        for j=1:2*N+1
            tabJln[i,j] = 0.0
        end
    end
end

###############################################
# Fill tabLegendre table
###############################################

"""
    tabLegendre!()

Fills the table `tabLegendre` with mid-point rule evaluation of the normalized Legendre associated functions.
See `computePlmcostheta` of the `SphericalHarmonics` for a definition.
"""
function tabLegendre!()
    for k=1:nbK
        xik = -1 +(2/nbK)*(k-0.5)
        thk = acos(xik)
        P = computePlmcostheta(thk,m+N)
        for l=m:m+N
            tabLegendre[l-m+1,k] = sqrt(pi)*P[(l,m)]
        end
    end
end

"""
    tabLegendre_clear!()

Clears the table `tabLegendre`.
"""
function tabLegendre_clear!()
    for i=1:N+1
        for k=1:nbK
            tabLegendre[i,k] = 0.0
        end
    end
end

###############################################
# Fill tabOmega table
###############################################

"""
    tabOmega!(x,q,[args])

Fills the table `tabOmega` with mid-point rule evaluation of the dimensionless angular velocity.
"""
function tabOmega!(x::Float64=0.0, q::Float64=1.0,tabOmega=tabOmega_serial)
    for k=1:nbK
        xik = -1 +(2/nbK)*(k-0.5)
        tabOmega[k] = Omega(xik,x,q)
    end
end

"""
    tabOmega_clear!([args])

Clears the table `tabOmega`.
"""
function tabOmega_clear!(tabOmega=tabOmega_serial)
    for k=1:nbK
        tabOmega[k] = 0.0
    end
end

###############################################
# Fill tabAlphaSqOverTwoOmega table
###############################################

"""
    tabKappaSqOverTwoOmega!()

Fills the table `tabOmega` with mid-point rule evaluation of the dimensionless ratio `kappa^2/(2Omega)`.
"""
function tabKappaSqOverTwoOmega!(x::Float64=0.0, q::Float64=1.0,tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial)
    for k=1:nbK
        xik = -1 +(2/nbK)*(k-0.5)
        tabKappaSqOverTwoOmega[k] = kappaSqOverTwoOmega(xik,x,q)
    end
end

"""
    tabKappaSqOverTwoOmega_clear!([args])

Clears the table `tabKappaSqOverTwoOmega`.
"""
function tabKappaSqOverTwoOmega_clear!(tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial)
    for k=1:nbK
        tabKappaSqOverTwoOmega[k] = 0.0
    end
end

###############################################
# Fill matrix elements table
###############################################

"""
    tabAlnFln!([args])

Fills the tables `tabAln` and `tabFln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
function tabAlnFln!(tabAln=tabAln_serial,tabFln=tabFln_serial,tabOmega=tabOmega_serial)
    for l=m:m+N
        for n=l:m+N
            integral = 0
            for k=1:nbK
                integral += tabLegendre[l-m+1,k]*tabOmega[k]*tabLegendre[n-m+1,k]
            end
            integral *= 2/nbK
            tabAln[l-m+1,n-m+1] = abs(m)*integral
            tabFln[l-m+1,n-m+1] = 2*integral
            if (n>l)
                tabAln[n-m+1,l-m+1] = tabAln[l-m+1,n-m+1]
                tabFln[n-m+1,l-m+1] = tabFln[l-m+1,n-m+1]
            end
        end
    end
end

"""
    tabBln!(x,[args])

Fills the tables `tabBln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
function tabBln!(x::Float64=0.0,tabBln=tabBln_serial)
    for l=m:m+N
        for n=m:m+N
            tabBln[l-m+1,n-m+1] = sqrt(1-x)*((1/2)*(sqrt(((2*l+1)*(l+m+1)*(l-m+1))/(2*l+3))*tabJln[l-m+1+1,n-m+1]
                                + tabJln[l-m+1,n-m+1] ))
            if (l>m)
                tabBln[l-m+1,n-m+1] -= sqrt(1-x)*(1/2)*sqrt(((2*l+1)*(l+m)*(l-m))/(2*l-1))*tabJln[l-m+1-1,n-m+1]
            end
        end
    end
end

"""
    tabCln!(x,[args])

Fills the tables `tabCln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
function tabCln!(x::Float64=0.0,tabCln=tabCln_serial)
    for l=m:m+N
        for n=m:m+N
            tabCln[l-m+1,n-m+1] = m*sqrt(1-x)*tabJln[l-m+1,n-m+1]
        end
    end
end

"""
    tabDln!(x,[args])

Fills the tables `tabDln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
function tabDln!(x::Float64=0.0, tabDln=tabDln_serial)
    for l=m:m+N
        for n=m:m+N
            tabDln[l-m+1,n-m+1] = sqrt(1-x)*((1/2)*(1/(2*n+1)-(eps/3)/((1-x)^(2/3)))*((-sqrt(((2*n+1)*(n+m+1)*(n-m+1))/(2*n+3))
                                *tabIln[l-m+1,n-m+1+1])- tabIln[l-m+1,n-m+1] ))
            if (n>m)
                tabDln[l-m+1,n-m+1] += sqrt(1-x)*((1/2)*(1/(2*n+1)-(eps/3)/((1-x)^(2/3)))* (sqrt(((2*n+1)*(n+m)*(n-m))/(2*n-1))
                                       *tabIln[l-m+1,n-m+1-1]))
            end
        end
    end
end

"""
    tabGln!(x,[args])

Fills the tables `tabGln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
function tabGln!(x::Float64=0.0,tabGln=tabGln_serial)
    for l=m:m+N
        for n=m:m+N
            tabGln[l-m+1,n-m+1] = -m*sqrt(1-x)*(1/(2*n+1)-(eps/3)/((1-x)^(2/3)))*tabIln[l-m+1,n-m+1]
        end
    end
end

"""
    tabHln!([args])

Fills the tables `tabHln`.
See `Aoki et al. (1979)` for the case `x=0` and `q=1`.
"""
function tabHln!(tabHln=tabHln_serial,tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial)
    for l=m:m+N
        for n=l:m+N
            integral = 0
            for k=1:nbK
                integral += tabLegendre[l-m+1,k]*tabKappaSqOverTwoOmega[k]*tabLegendre[n-m+1,k]
            end
            integral *= 2/nbK
            tabHln[l-m+1,n-m+1] = integral
            if (n>l)
                tabHln[n-m+1,l-m+1] = tabHln[l-m+1,n-m+1]
            end
        end
    end
end

###############################################
# Fill matrix elements table
###############################################

"""
    table_constant_fill!()

Fills the recursion tables of `I`, `Ia` and `J`.
"""
function table_constant_fill!()
    tabIln!()
    tabIaln!()
    tabJln!()
    tabLegendre!()
end

"""
    table_function_fill!(x,q,[args])

Fills the function tables (angular velocity and epicyclic frequency).
"""
function table_function_fill!(x::Float64=0.0, q::Float64=1.0,
    tabOmega=tabOmega_serial, tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial)
    tabOmega!(x,q,tabOmega)
    tabKappaSqOverTwoOmega!(x,q,tabKappaSqOverTwoOmega)
end

"""
    matrix_fill!(x,q,[args])

Fills the matrix elements.
"""
function matrix_fill!(x::Float64=0.0, q::Float64=1.0,
    tabAln=tabAln_serial, tabBln=tabBln_serial, tabCln=tabCln_serial,
    tabDln=tabDln_serial, tabFln=tabFln_serial, tabGln=tabGln_serial,
    tabHln=tabHln_serial, tabOmega=tabOmega_serial,
    tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial)
    tabAlnFln!(tabAln,tabFln,tabOmega)
    tabBln!(x,tabBln)
    tabCln!(x,tabCln)
    tabDln!(x,tabDln)
    tabGln!(x,tabGln)
    tabHln!(tabHln,tabKappaSqOverTwoOmega)
end

######################################
# Clear table
######################################

"""
    tabAln_clear!([args])

Clears the matrix elements of `Aln`.
"""
function tabAln_clear!(tabAln=tabAln_serial)
    for i=1:N+1
        for j=1:N+1
            tabAln[i,j] = 0.0
        end
    end
end

"""
    tabBln_clear!([args])

Clears the matrix elements of `Bln`.
"""
function tabBln_clear!(tabBln=tabBln_serial)
    for i=1:N+1
        for j=1:N+1
            tabBln[i,j] = 0.0
        end
    end
end

"""
    tabCln_clear!([args])

Clears the matrix elements of `Cln`.
"""
function tabCln_clear!(tabCln=tabCln_serial)
    for i=1:N+1
        for j=1:N+1
            tabCln[i,j] = 0.0
        end
    end
end

"""
    tabDln_clear!([args])

Clears the matrix elements of `Dln`.
"""
function tabDln_clear!(tabDln=tabDln_serial)
    for i=1:N+1
        for j=1:N+1
            tabDln[i,j] = 0.0
        end
    end
end

"""
    tabFln_clear!([args])

Clears the matrix elements of `Fln`.
"""
function tabFln_clear!(tabFln=tabFln_serial)
    for i=1:N+1
        for j=1:N+1
            tabFln[i,j] = 0.0
        end
    end
end

"""
    tabGln_clear!([args])

Clears the matrix elements of `Gln`.
"""
function tabGln_clear!(tabGln=tabGln_serial)
    for i=1:N+1
        for j=1:N+1
            tabGln[i,j] = 0.0
        end
    end
end

"""
    tabHln_clear!([args])

Clears the matrix elements of `Hln`.
"""
function tabHln_clear!(tabHln=tabHln_serial)
    for i=1:N+1
        for j=1:N+1
            tabHln[i,j] = 0.0
        end
    end
end

"""
    table_function_clear!([args])

Clears the function tables.
"""
function table_function_clear!(tabOmega=tabOmega_serial,
    tabKappaSqOverTwoOmega=tabKappaSqOverTwoOmega_serial)
    tabOmega_clear!(tabOmega)
    tabKappaSqOverTwoOmega_clear!(tabKappaSqOverTwoOmega)
end

"""
    matrix_clear!([args])

Clears all the matrix tables.
"""
function matrix_clear!(tabAln=tabAln_serial, tabBln=tabBln_serial,
    tabCln=tabCln_serial, tabDln=tabDln_serial, tabFln=tabFln_serial,
    tabGln=tabGln_serial, tabHln=tabHln_serial)
    tabAln_clear!(tabAln)
    tabBln_clear!(tabBln)
    tabCln_clear!(tabCln)
    tabDln_clear!(tabDln)
    tabFln_clear!(tabFln)
    tabGln_clear!(tabGln)
    tabHln_clear!(tabHln)
end

######################################
# Initialize constant Tables
######################################

table_constant_fill!()

max_stars_repo_name: KerwannTEP/LiRA

max_issues_repo_name: KerwannTEP/LiRA

max_forks_repo_name: KerwannTEP/LiRA


text:
@testset "chaintransform" begin
    rng = MersenneTwister(123546)

    P = rand(rng, 3, 2)
    tp = LinearTransform(P)

    f(x) = sin.(x)
    tf = FunctionTransform(f)

    t = ChainTransform([tp, tf])

    # Check composition constructors.
    @test (tf ∘ ChainTransform([tp])).transforms == [tp, tf]
    @test (ChainTransform([tf]) ∘ tp).transforms == [tp, tf]

    # Verify correctness.
    x = ColVecs(randn(rng, 2, 3))
    x′ = map(t, x)

    @test all([t(x[n]) ≈ f(P * x[n]) for n in eachindex(x)])
    @test all([t(x[n]) ≈ x′[n] for n in eachindex(x)])

    # Verify printing works as expected.
    @test repr(tp ∘ tf) == "Chain of 2 transforms:\n\t - $(tf) |> $(tp)"
    test_ADs(
        x -> SEKernel() ∘ (ScaleTransform(exp(x[1])) ∘ ARDTransform(exp.(x[2:4]))),
        randn(rng, 4);
        ADs=[:ForwardDiff, :ReverseDiff],  # explicitly pass ADs to exclude :Zygote
    )
    @test_broken "test_AD of chain transform is currently broken in Zygote, see GitHub issue #263"
end


text:
"""
    UniversalFallback

The `UniversalFallback` can be applied on a [`MathOptInterface.ModelLike`](@ref)
`model` to create the model `UniversalFallback(model)` supporting *any*
constraint and attribute. This allows to have a specialized implementation in
`model` for performance critical constraints and attributes while still
supporting other attributes with a small performance penalty. Note that `model`
is unaware of constraints and attributes stored by `UniversalFallback` so this
is not appropriate if `model` is an optimizer (for this reason,
[`MathOptInterface.optimize!`](@ref) has not been implemented). In that case,
optimizer bridges should be used instead.
"""
mutable struct UniversalFallback{MT} <: MOI.ModelLike
    model::MT
    objective::Union{MOI.AbstractScalarFunction, Nothing}
    constraints::OrderedDict{Tuple{DataType, DataType}, OrderedDict} # See https://github.com/jump-dev/JuMP.jl/issues/1152 and https://github.com/jump-dev/JuMP.jl/issues/2238
    nextconstraintid::Int64
    con_to_name::Dict{CI, String}
    name_to_con::Union{Dict{String, MOI.ConstraintIndex}, Nothing}
    optattr::Dict{MOI.AbstractOptimizerAttribute, Any}
    modattr::Dict{MOI.AbstractModelAttribute, Any}
    varattr::Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}
    conattr::Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}
    function UniversalFallback{MT}(model::MOI.ModelLike) where {MT}
        new{typeof(model)}(model,
                           nothing,
                           OrderedDict{Tuple{DataType, DataType}, OrderedDict}(),
                           0,
                           Dict{CI, String}(),
                           nothing,
                           Dict{MOI.AbstractOptimizerAttribute, Any}(),
                           Dict{MOI.AbstractModelAttribute, Any}(),
                           Dict{MOI.AbstractVariableAttribute, Dict{VI, Any}}(),
                           Dict{MOI.AbstractConstraintAttribute, Dict{CI, Any}}())
    end
end
UniversalFallback(model::MOI.ModelLike) = UniversalFallback{typeof(model)}(model)

function Base.show(io::IO, U::UniversalFallback)
    s(n) = n == 1 ? "" : "s"
    indent = " "^get(io, :indent, 0)
    MOIU.print_with_acronym(io, summary(U))
    !(U.objective === nothing) && print(io, "\n$(indent)with objective")
    for (attr, name) in (   (U.constraints, "constraint"),
                            (U.optattr, "optimizer attribute"),
                            (U.modattr, "model attribute"),
                            (U.varattr, "variable attribute"),
                            (U.conattr, "constraint attribute") )
        n = length(attr)
        if n > 0
            print(io, "\n$(indent)with $n $name$(s(n))")
        end
    end
    print(io, "\n$(indent)fallback for ")
    show(IOContext(io, :indent => get(io, :indent, 0)+2), U.model)
end

function MOI.is_empty(uf::UniversalFallback)
    return MOI.is_empty(uf.model) && uf.objective === nothing && isempty(uf.constraints) &&
        isempty(uf.modattr) && isempty(uf.varattr) && isempty(uf.conattr)
end
function MOI.empty!(uf::UniversalFallback)
    MOI.empty!(uf.model)
    uf.objective = nothing
    empty!(uf.constraints)
    uf.nextconstraintid = 0
    empty!(uf.con_to_name)
    uf.name_to_con = nothing
    empty!(uf.modattr)
    empty!(uf.varattr)
    empty!(uf.conattr)
end
function MOI.copy_to(uf::UniversalFallback, src::MOI.ModelLike; kws...)
    MOIU.automatic_copy_to(uf, src; kws...)
end
function supports_default_copy_to(uf::UniversalFallback, copy_names::Bool)
    return supports_default_copy_to(uf.model, copy_names)
end

# References
MOI.is_valid(uf::UniversalFallback, idx::VI) = MOI.is_valid(uf.model, idx)
function MOI.is_valid(uf::UniversalFallback, idx::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.is_valid(uf.model, idx)
    else
        haskey(uf.constraints, (F, S)) && haskey(uf.constraints[(F, S)], idx)
    end
end
function MOI.delete(uf::UniversalFallback, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.delete(uf.model, ci)
    else
        if !MOI.is_valid(uf, ci)
            throw(MOI.InvalidIndex(ci))
        end
        delete!(uf.constraints[(F, S)], ci)
        delete!(uf.con_to_name, ci)
        uf.name_to_con = nothing
    end
    for d in values(uf.conattr)
        delete!(d, ci)
    end
end
function _remove_variable(uf::UniversalFallback,
                          constraints::OrderedDict{<:CI{MOI.SingleVariable}}, vi::VI)
    to_delete = keytype(constraints)[]
    for (ci, constraint) in constraints
        f::MOI.SingleVariable = constraint[1]
        if f.variable == vi
            push!(to_delete, ci)
        end
    end
    MOI.delete(uf, to_delete)
end
function _remove_variable(uf::UniversalFallback,
                          constraints::OrderedDict{CI{MOI.VectorOfVariables, S}},
                          vi::VI) where S
    to_delete = keytype(constraints)[]
    for (ci, constraint) in constraints
        f::MOI.VectorOfVariables, s = constraint
        if vi in f.variables
            if length(f.variables) > 1
                if MOI.supports_dimension_update(S)
                    constraints[ci] = remove_variable(f, s, vi)
                else
                    throw_delete_variable_in_vov(vi)
                end
            else
                push!(to_delete, ci)
            end
        end
    end
    MOI.delete(uf, to_delete)
end
function _remove_variable(::UniversalFallback, constraints::OrderedDict{<:CI}, vi::VI)
    for (ci, constraint) in constraints
        f, s = constraint
        constraints[ci] = remove_variable(f, s, vi)
    end
end
function _remove_vector_of_variables(
    uf::UniversalFallback, constraints::OrderedDict{<:CI{MOI.VectorOfVariables}},
    vis::Vector{VI}
)
    to_delete = keytype(constraints)[]
    for (ci, constraint) in constraints
        f::MOI.VectorOfVariables = constraint[1]
        if vis == f.variables
            push!(to_delete, ci)
        end
    end
    MOI.delete(uf, to_delete)
end
function _remove_vector_of_variables(
    ::UniversalFallback, ::OrderedDict{<:CI}, ::Vector{VI})
end
function MOI.delete(uf::UniversalFallback, vi::VI)
    MOI.delete(uf.model, vi)
    for d in values(uf.varattr)
        delete!(d, vi)
    end
    if uf.objective !== nothing
        uf.objective = remove_variable(uf.objective, vi)
    end
    for (_, constraints) in uf.constraints
        _remove_variable(uf, constraints, vi)
    end
end
function MOI.delete(uf::UniversalFallback, vis::Vector{VI})
    MOI.delete(uf.model, vis)
    for d in values(uf.varattr)
        for vi in vis
            delete!(d, vi)
        end
    end
    if uf.objective !== nothing
        uf.objective = remove_variable(uf.objective, vis)
    end
    for (_, constraints) in uf.constraints
        _remove_vector_of_variables(uf, constraints, vis)
        for vi in vis
            _remove_variable(uf, constraints, vi)
        end
    end
end

# Attributes
_get(uf, attr::MOI.AbstractOptimizerAttribute)          = uf.optattr[attr]
_get(uf, attr::MOI.AbstractModelAttribute)              = uf.modattr[attr]
function _get(uf, attr::MOI.AbstractVariableAttribute, vi::VI)
    attribute_dict = get(uf.varattr, attr, nothing)
    if attribute_dict === nothing
        # It means the attribute is not set to any variable so in particular, it
        # is not set for `vi`
        return nothing
    end
    return get(attribute_dict, vi, nothing)
end
function _get(uf, attr::MOI.AbstractConstraintAttribute, ci::CI)
    attribute_dict = get(uf.conattr, attr, nothing)
    if attribute_dict === nothing
        # It means the attribute is not set to any constraint so in particular,
        # it is not set for `ci`
        return nothing
    end
    return get(attribute_dict, ci, nothing)
end
function _get(uf, attr::MOI.CanonicalConstraintFunction, ci::MOI.ConstraintIndex)
    return MOI.get_fallback(uf, attr, ci)
    func = MOI.get(uf, MOI.ConstraintFunction(), ci)
    if is_canonical(func)
        return func
    else
        return canonical(func)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::Union{MOI.AbstractOptimizerAttribute,
                             MOI.AbstractModelAttribute})
    if !MOI.is_copyable(attr) || MOI.supports(uf.model, attr)
        MOI.get(uf.model, attr)
    else
        _get(uf, attr)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::MOI.AbstractConstraintAttribute,
                 idx::MOI.ConstraintIndex{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S) &&
        (!MOI.is_copyable(attr) || MOI.supports(uf.model, attr, typeof(idx)))
        MOI.get(uf.model, attr, idx)
    else
        _get(uf, attr, idx)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::MOI.AbstractVariableAttribute, idx::MOI.VariableIndex)
    if !MOI.is_copyable(attr) || MOI.supports(uf.model, attr, typeof(idx))
        MOI.get(uf.model, attr, idx)
    else
        _get(uf, attr, idx)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::MOI.NumberOfConstraints{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.get(uf.model, attr)
    else
        return length(get(uf.constraints, (F, S), OrderedDict{CI{F, S}, Tuple{F, S}}()))
    end
end
function MOI.get(uf::UniversalFallback,
                 listattr::MOI.ListOfConstraintIndices{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, listattr)
    else
        collect(keys(get(uf.constraints, (F, S), OrderedDict{CI{F, S}, Tuple{F, S}}())))
    end
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfConstraints)
    list = MOI.get(uf.model, listattr)
    for (FS, constraints) in uf.constraints
        if !isempty(constraints)
            push!(list, FS)
        end
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfOptimizerAttributesSet)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.optattr)
        push!(list, attr)
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfModelAttributesSet)
    list = MOI.get(uf.model, listattr)
    if uf.objective !== nothing
        push!(list, MOI.ObjectiveFunction{typeof(uf.objective)}())
    end
    for attr in keys(uf.modattr)
        push!(list, attr)
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfVariableAttributesSet)
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.varattr)
        push!(list, attr)
    end
    list
end
function MOI.get(uf::UniversalFallback, listattr::MOI.ListOfConstraintAttributesSet{F, S}) where {F, S}
    list = MOI.get(uf.model, listattr)
    for attr in keys(uf.conattr)
        push!(list, attr)
    end
    return list
end

# Objective
function MOI.set(uf::UniversalFallback, attr::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense) where T
    if sense == MOI.FEASIBILITY_SENSE
        uf.objective = nothing
    end
    MOI.set(uf.model, attr, sense)
end
function MOI.get(uf::UniversalFallback,
                 attr::MOI.ObjectiveFunctionType)
    if uf.objective === nothing
        return MOI.get(uf.model, attr)
    else
        return typeof(uf.objective)
    end
end
function MOI.get(uf::UniversalFallback,
                 attr::MOI.ObjectiveFunction{F})::F where F
    if uf.objective === nothing
        return MOI.get(uf.model, attr)
    else
        return uf.objective
    end
end
function MOI.set(uf::UniversalFallback,
                 attr::MOI.ObjectiveFunction,
                 func::MOI.AbstractScalarFunction)
    if MOI.supports(uf.model, attr)
        MOI.set(uf.model, attr, func)
        # Clear any fallback objective
        uf.objective = nothing
    else
        uf.objective = copy(func)
        # Clear any `model` objective
        sense = MOI.get(uf.model, MOI.ObjectiveSense())
        MOI.set(uf.model, MOI.ObjectiveSense(), MOI.FEASIBILITY_SENSE)
        MOI.set(uf.model, MOI.ObjectiveSense(), sense)
    end
end

function MOI.modify(uf::UniversalFallback, obj::MOI.ObjectiveFunction, change::MOI.AbstractFunctionModification) where F
    if uf.objective === nothing
        MOI.modify(uf.model, obj, change)
    else
        uf.objective = modify_function(uf.objective, change)
    end
end

# Name
# The names of constraints not supported by `uf.model` need to be handled
function MOI.set(uf::UniversalFallback, attr::MOI.ConstraintName, ci::CI{F, S}, name::String) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, attr, ci, name)
    else
        uf.con_to_name[ci] = name
        uf.name_to_con = nothing # Invalidate the name map.
    end
    return
end
function MOI.get(uf::UniversalFallback, attr::MOI.ConstraintName, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.get(uf.model, attr, ci)
    else
        return get(uf.con_to_name, ci, EMPTYSTRING)
    end
end

MOI.get(uf::UniversalFallback, ::Type{VI}, name::String) = MOI.get(uf.model, VI, name)

check_type_and_multiple_names(::Type, ::Nothing, ::Nothing, name) = nothing
check_type_and_multiple_names(::Type{T}, value::T, ::Nothing, name) where T = value
check_type_and_multiple_names(::Type, ::Any, ::Nothing, name) where T = nothing
check_type_and_multiple_names(::Type{T}, ::Nothing, value::T, name) where T = value
check_type_and_multiple_names(::Type, ::Nothing, ::Any, name) where T = nothing
function check_type_and_multiple_names(T::Type, ::Any, ::Any, name)
    throw_multiple_name_error(T, name)
end
function MOI.get(uf::UniversalFallback, ::Type{CI{F, S}}, name::String) where {F, S}
    if uf.name_to_con === nothing
        uf.name_to_con = build_name_to_con_map(uf.con_to_name)
    end
    if MOI.supports_constraint(uf.model, F, S)
        ci = MOI.get(uf.model, CI{F, S}, name)
    else
        # There is no `F`-in-`S` constraint in `b.model`, `ci` is only queried
        # to check for duplicate names.
        ci = MOI.get(uf.model, CI, name)
    end
    ci_uf = get(uf.name_to_con, name, nothing)
    throw_if_multiple_with_name(ci_uf, name)
    return check_type_and_multiple_names(CI{F, S}, ci_uf, ci, name)
end
function MOI.get(uf::UniversalFallback, ::Type{CI}, name::String)
    if uf.name_to_con === nothing
        uf.name_to_con = build_name_to_con_map(uf.con_to_name)
    end
    ci_uf = get(uf.name_to_con, name, nothing)
    throw_if_multiple_with_name(ci_uf, name)
    return check_type_and_multiple_names(
        CI, ci_uf, MOI.get(uf.model, CI, name), name)
end

_set(uf, attr::MOI.AbstractOptimizerAttribute, value) = uf.optattr[attr] = value
_set(uf, attr::MOI.AbstractModelAttribute, value)     = uf.modattr[attr] = value
function _set(uf, attr::MOI.AbstractVariableAttribute, vi::VI, value)
    if !haskey(uf.varattr, attr)
        uf.varattr[attr] = Dict{VI, Any}()
    end
    uf.varattr[attr][vi] = value
end
function _set(uf, attr::MOI.AbstractConstraintAttribute, ci::CI, value)
    if !haskey(uf.conattr, attr)
        uf.conattr[attr] = Dict{CI, Any}()
    end
    uf.conattr[attr][ci] = value
end
MOI.supports(::UniversalFallback, ::Union{MOI.AbstractModelAttribute, MOI.AbstractOptimizerAttribute}) = true
function MOI.set(uf::UniversalFallback, attr::Union{MOI.AbstractOptimizerAttribute, MOI.AbstractModelAttribute}, value)
    if MOI.supports(uf.model, attr)
        return MOI.set(uf.model, attr, value)
    else
        return _set(uf, attr, value)
    end
end
MOI.supports(::UniversalFallback, ::Union{MOI.AbstractVariableAttribute, MOI.AbstractConstraintAttribute}, ::Type{<:MOI.Index}) = true
function MOI.set(uf::UniversalFallback, attr::MOI.AbstractVariableAttribute, idx::VI, value)
    if MOI.supports(uf.model, attr, typeof(idx))
        return MOI.set(uf.model, attr, idx, value)
    else
        return _set(uf, attr, idx, value)
    end
end
function MOI.set(uf::UniversalFallback, attr::MOI.AbstractConstraintAttribute, idx::CI{F, S}, value) where {F, S}
    if MOI.supports_constraint(uf.model, F, S) && MOI.supports(uf.model, attr, CI{F, S})
        return MOI.set(uf.model, attr, idx, value)
    else
        return _set(uf, attr, idx, value)
    end
end

# Constraints
MOI.supports_constraint(uf::UniversalFallback, ::Type{F}, ::Type{S}) where {F<:MOI.AbstractFunction, S<:MOI.AbstractSet} = true
function _new_constraint_index(uf, f::MOI.SingleVariable, s::MOI.AbstractScalarSet)
    return CI{MOI.SingleVariable, typeof(s)}(f.variable.value)
end
function _new_constraint_index(uf, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    uf.nextconstraintid += 1
    return CI{typeof(f), typeof(s)}(uf.nextconstraintid)
end
function MOI.add_constraint(uf::UniversalFallback, f::MOI.AbstractFunction, s::MOI.AbstractSet)
    F = typeof(f)
    S = typeof(s)
    if MOI.supports_constraint(uf.model, F, S)
        return MOI.add_constraint(uf.model, f, s)
    else
        constraints = get!(uf.constraints, (F, S)) do
            OrderedDict{CI{F, S}, Tuple{F, S}}()
        end::OrderedDict{CI{F, S}, Tuple{F, S}}
        ci = _new_constraint_index(uf, canonical(f), copy(s))
        constraints[ci] = (f, s)
        return ci
    end
end
function MOI.modify(uf::UniversalFallback, ci::CI{F, S}, change::MOI.AbstractFunctionModification) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.modify(uf.model, ci, change)
    else
        (f, s) = uf.constraints[(F, S)][ci]
        uf.constraints[(F, S)][ci] = (modify_function(f, change), s)
    end
end

function MOI.get(uf::UniversalFallback, attr::MOI.ConstraintFunction, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, attr, ci)
    else
        MOI.throw_if_not_valid(uf, ci)
        uf.constraints[(F, S)][ci][1]
    end
end
function MOI.get(uf::UniversalFallback, attr::MOI.ConstraintSet, ci::CI{F, S}) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.get(uf.model, attr, ci)
    else
        MOI.throw_if_not_valid(uf, ci)
        uf.constraints[(F, S)][ci][2]
    end
end
function MOI.set(uf::UniversalFallback, ::MOI.ConstraintFunction, ci::CI{F,S}, func::F) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, MOI.ConstraintFunction(), ci, func)
    else
        MOI.throw_if_not_valid(uf, ci)
        if F == MOI.SingleVariable
            throw(MOI.SettingSingleVariableFunctionNotAllowed())
        end
        (_, s) = uf.constraints[(F, S)][ci]
        uf.constraints[(F, S)][ci] = (func, s)
    end
end
function MOI.set(uf::UniversalFallback, ::MOI.ConstraintSet, ci::CI{F,S}, set::S) where {F, S}
    if MOI.supports_constraint(uf.model, F, S)
        MOI.set(uf.model, MOI.ConstraintSet(), ci, set)
    else
        MOI.throw_if_not_valid(uf, ci)
        (f, _) = uf.constraints[(F, S)][ci]
        uf.constraints[(F, S)][ci] = (f, set)
    end
end

# Variables
MOI.add_variable(uf::UniversalFallback) = MOI.add_variable(uf.model)
MOI.add_variables(uf::UniversalFallback, n) = MOI.add_variables(uf.model, n)


text:
using Plots
using CSV
using Printf
using DataFrames
using LaTeXStrings
using LsqFit            # Required for the statistics on the fitting process
using FrictionModels

# Read in a CSV file containing the data to fit.
#
# The file is structured as follows:
# Column 1: Time
# Column 2: Velocity
# Column 3: Normal force
# Column 4: Friction Force
println(pwd())
data = CSV.read("examples//test_data_1.csv", DataFrame)

# Construct an initial estimate for model parameters
init_guess = LuGreModel(
    0.25,
    0.15,
    0.01,
    1e6,
    1e3,
    0.1
)

# We can specify limits on each parameter.  If no limit is desired for A
# specific parameter, we can simply input -Inf or Inf for either a lower or
# upper constraint respectively.
lb = [0.2, 0.05, 0.0, 1e3, 0.0, 0.0]
ub = [1.0, 0.2, Inf, Inf, Inf, Inf]

# Fit the model
results = fit_model(
    init_guess,     # Initial Guess
    data[:,1],      # Time
    data[:,4],      # Friction Force
    data[:,3],      # Normal Force
    data[:,2],      # Velocity
    [0.0],          # Initial Condition Vector
    lower = lb,     # Lower Bounds (optional)
    upper = ub      # Upper Bounds (optional)
)

# Solve the model
rsp = friction(results.model, data[:,1], data[:,3], data[:,2], [0.0])

# Plot the data
plt = plot(
    data[:,2], rsp.f,
    xlabel = L"v(t)",
    ylabel = L"F(t)",
    label = "Fitted Model",
    lw = 2,
    legend = :bottomright
)
plot!(
    plt,
    data[:,2], data[:,4],
    st = :scatter,
    label = "Measured Data"
)
display(plt)

# Print out the coefficients of the fitted model
@printf("Fitted Model Coefficients:\n")

# Print out statistitics related to the fit
sigma = stderror(results.fit)
confidence = confidence_interval(results.fit)

# Print out the results
@printf(
    "static_coefficient: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.static_coefficient,
    sigma[1],
    confidence[1][1],
    confidence[1][2]
)
@printf(
    "coulomb_coefficient: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.coulomb_coefficient,
    sigma[2],
    confidence[2][1],
    confidence[2][2]
)
@printf(
    "stribeck_velocity: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.stribeck_velocity,
    sigma[3],
    confidence[3][1],
    confidence[3][2]
)
@printf(
    "bristle_stiffness: %e\n\tError: %e\n\tConfidence Interval: (%e, %e)\n",
    results.model.bristle_stiffness,
    sigma[4],
    confidence[4][1],
    confidence[4][2]
)
@printf(
    "bristle_damping: %e\n\tError: %e\n\tConfidence Interval: (%e, %e)\n",
    results.model.bristle_damping,
    sigma[5],
    confidence[5][1],
    confidence[5][2]
)
@printf(
    "viscous_damping: %f\n\tError: %f\n\tConfidence Interval: (%f, %f)\n",
    results.model.viscous_damping,
    sigma[6],
    confidence[6][1],
    confidence[6][2]
)


text:
#=
20 named functions ( 19 without counting the mit function)  (counting helper functions)
10 helper functions

Divide & Conquer
Does a Bisection Sort on an array where it measures the length of
our Objective Function (of the Interval vector) which can either be:=#
#------------------
#= A:
 if it is equal to 1,
    1. first compare the interval's  contents
    if left value (of a vector) is higher than the right one, swap vector contents
    2.then we have to stop (as no more middles are available to explore)
    congratulations, you've reached the end #HALT! =#
#------------------
#= B: However if distance between Interval bounds is higher than 1, hence, there is at least 1 new middle to compute & find out:
    1.Compute the middle: either there is a whole number, representing an actual index, with Intervals (a,mid) , (mid+1, b) [2 middles mid , mid+1]
     or we've got a fractional middle, that needs to be ceiled (as above) & floored (as below) with intervals:  (a, below),... ,(above, b) (2 middles below & above):

    -new: further care with Interval processing is required: say below+1 & above+1 exits & unique then ranges should be: (a, below), (below+1,above), (above+1, b) #TODO or
    -- Assuming: below-1 & above-1 exits & unique,
        then the Intervals: (a,below-1), (below, above-1), (above, b) [if each Interval bound exist, as a whole Int Number]
    --- need to introduce another helper metric, Captures & stores the speed (of algorithm). as both workflows are valid, it would all boil down to the number #1: speed of the Algorithm #better-benchmark

    [this bifurcates the workflow of the program, hence, try one function at once, [the validity of an Interval is the Priority]   ]
    at each Interval-building, we also check the contents of each bounds
    2. for each new Interval (we got ), compare it's contents, check distance (& we're back to square 1, to A ) #Done!
# Are the middles compared too (well, as we'd be comparing every Interval, that would be inlusive of generated middles) =#
#------------------
#=if middlle:
    - is Interval even?: Yes -> has 1 middle #now-here
    - if not: there's a fractionalMid, approximate index to get below & above
else
 go left middle(a,b,arr) # mid, mid+1 (isEven), OR # below & above (!isEven)

 go right middle(a,b,arr) # mid, mid+1 (isEven), OR # below & above (!isEven)

 #tip: (execution) time can be misleading

-
 -buildAboveSoBelow(below, above, a)

 go right
 go left

sumInt = sumInterval(a,b)
  isEven2() #true

=#

import Base: @propagate_inbounds, @inbounds # compiler inlines a function, while retaining the caller's inbounds context
import BenchmarkTools: @btime, @time, @benchmark
UnexpMsg = "ERROR: unexpected input:  please check input arguments , then try again  "
positiveMsg = "ERROR: Only Positive input arguments are allowed - please check input arguments, then try again"
outofBoundsMsg = "ERROR: input's length is larger than original vector- kindly check again "

#@benchmark  +(10000, 1000000) #simple + is WAY much Bloated, it isn't simple, anymore :S
#@benchmark sumInterval(10000, 1000000) =#

@propagate_inbounds isPositive(num) =  0 < num ? true : false

#Euclidean Distnace https://www.sciencedirect.com/topics/mathematics/euclidean-distance#:~:text=The%20weighted%20Euclidean%20distance%20measure,element%20of%20the%20average%20profile.

    #isnumeric(a) && isnumeric(b) && a > 0 && b > 0 ? abs(max(a, b) - min(a, b)) :
    #positiveMsg

#@benchmark sumInterval()
@propagate_inbounds function sumInterval(a, b) # subtracts 1 (complies with julia logic)
    cond = abs(a + b) - 1
    @inbounds if cond >= 1
        return cond
    end
end

#--------------
#3 if-statement-logic (if-elseif-else)

""" checks evinity of a bound between point a to point b   #TODO: Research  >= v.s. >
```input:
a: lower bound
b: upper bound

calculates eucledian distance between a & b
1. checks if condition is positive
2. checks evenity

```
```output:
true  if it's positive & even
false: if either #improve
```
"""
@propagate_inbounds function isEven2(a = 1, b = 10) # = readable #preferred #fast  #optimized  #2 if-statements
    try
        @inbounds if a >= 0 && b >= 0   #==#
            distance = sumInterval(a, b) #1.distance
            evenCond = distance % 2 == 0 ? true : false
            _isPositive = isPositive(distance)
            @inbounds if evenCond && _isPositive #ok
                return true
            elseif !evenCond || !_isPositive
                return false
            end
        else
            throw(error("Positive number error"))
        end
    catch notPositiveInputError
        #Exception
        @error positiveMsg exception = (notPositiveInputError, catch_backtrace())
        println(positiveMsg) #positive arguments error
    end #==#
end

@propagate_inbounds function isEven2(m = 10) # = readable #preferred #fast  #optimiz-able  #2 if-statements
    try
        @inbounds if m >= 0  # && b >= 0   #==#
            distance = m  #sumInterval(m) #1.distance
            evenCond = distance % 2 == 0 ? true : false
            _isPositive = isPositive(distance)
            @inbounds if evenCond && _isPositive #ok
                return true
            elseif !evenCond || !_isPositive
                return false
            end
        else
            throw(error("Positive number error"))
        end
    catch notPositiveInputError

        @error positiveMsg exception = (notPositiveInputError, catch_backtrace())
        println(positiveMsg) #positive arguments error
    end #==#
end
#--------------------

isEven2(1, 10)# not divisible by 2  -then got fractionalMid , estimate: above & below

#@benchmark isEven1()
@benchmark isEven2(1, 10) #slightly does better

@benchmark isEven2() #due to optimization max time drops from 667.437 to 211.268 = 30.483
#in return of a `slight` incremental degredation in the overall mean = 39,860 =  9.860

#=
mM(s1)=minMax(s1)= Range (min … max):(min … max)= 30.553 ns … 667.437 ns
mM(s2)=minMax(s2) =Range (min … max): = 30.483 ns … 211.268 ns
mM(s1)= minMax(s1)= (mean ± σ)= 38.030 ns ±  18.992 ns
nM(s2)=minMax(s2)=(mean ± σ) = 39.860 ns ±  19.883 ns

diffmM(s1,s2) = mM(s2) - mM(s1) = 30.553 ns … 667.437 ns - 30.483 ns … 211.268 ns
mM(s2) - mM(s1) = 30.483 ns … 211.268 ns = 180.785
mM(s2) / mM(s1) =
mM(s2) / mM(s1) =

=#

#1stOrder
#RangeRatios
#Low

diffLow = 30.553 - 30.483 # lower ranges of both trials are near i.e. how
RatioLow= min(30.483,30.553) / max(30.483,30.553) # low sides almost exact 0.9977088992897588

#High
diff2sampleHigh= diffHigh = 667.437 - 211.268 # difference high, 2 samples  # # <--------------------
RatioHigh = max(30.483,30.553) / min(30.483,30.553) # RatioHigh  1.0022963619066365

#2nd order  - probabilty High difference
#dispersuib diffHigh - diffLow

diffHighLow =N = uncertainRangeHighLow = diffHigh - diffLow # <--------------------
lengthRatio =  diffLow /diffHigh

probdiffLow =   diffLow/N # inverse of that ratio Lo i.e. diffLo/ rangeHiLo = ProbabilityRatioLo #almost not happening
probdiffHigh = diffHigh / N #  probabilty High difference   1

#to put things into a metric, Ratio:
uncertainRatio = N / diffHigh #almost 1 (almost equal)

#try Infer ------------------------------------------------------------------------------------------------------------------------
lendiff_diff =  diffHigh - diffHighLow #<-------------------run this
diffHighLow ≈  diffHigh

Difference2sampleHigh =  diff2sampleHigh - diffHighLow
#------------------------------------------------------------------------------------------------------------

#2nd order Ratio
diff2sampleHigh / diffHighLow

#infer
diff_Difference2sampleHigh = Difference2sampleHigh - lendiff_diff # diffHigh - diffHighLow
diff_Difference2sampleHigh == Difference2sampleHigh - lendiff_diff

 diff2sampleHigh  -   diffHigh   - (diffHighLow  - diffHighLow )
(diff2sampleHigh)  -   diffHigh    +(diffHighLow  - diffHighLow )
#----------------------------------------------------------------------
diff2sampleHigh  -   diffHigh
diff2sampleHigh == diffHigh
diff2sampleHigh === diffHigh
#---------------------------------------
N = uncertainRangeHighLow = diffHigh - diffLow ≈  diffHigh

N === diffHigh - diffLow

#--------------------------------------------------------------------------------
#=concludes
isEven2 is slighly better than isEven1 - besides it's more compact, readable
one obvious reason is: number of if statements: isEven2 has only 2 ifs, while is isEven1 has 3
=#

@propagate_inbounds isEven(a, b) = isEven2(a, b)
@propagate_inbounds isEven(m) = isEven2(m)
#either 1 of the conditions are true all time

#Great handy function:
@propagate_inbounds ϟ(a, b) =
    (a > b) ⊻ (b > a) ? max(a, b) - min(a, b) : println(positiveMsg) #  abs(a - b)

#--- buildInterval

@propagate_inbounds function buildInterval(a, b)

    return collect((a, b))
end

vector = buildInterval(1, 2)

@propagate_inbounds function buildInterval(tuple)

    return collect((tuple[1], tuple[2]))
end

vector =  buildInterval(((1,2)))

function oldschoolSwap!(x, y)
    tmp = x
    x = y
    y = tmp
    return x, y
end
#----------------------

middles = [] # need to track middles

#----------------------
""" compares vector a, it's element at first index α with second element at index β

```input:
a: first index of comparison
b: second index of comparison
arr: the original vector array
```

```output:
an ordered tuple of the corrected indicies (of the vector array)

```
""" #requires: of indexOf,oldschoolSwap!

function compareVector(a = 1, b = 2, arr = [2, 1, 3, 4])
    response = nothing
    try #1. we call this function when we'd like to compare index α with index β of a Vector array  # do your thing

        firstContent = Int(findfirst(isequal(a), arr)) #indexOf(first)
        lastContent = Int(findfirst(isequal(b), arr)) #indexOf(last)

        if firstContent > lastContent # correct
            response = @inbounds arr[a], arr[b] = oldschoolSwap!(arr[a], arr[b]) #plain content swap in julia  #swap array contents directly

        elseif firstContent < lastContent #only possible - correct situation (to deal with)
            #Intent: skip
            return response = 1
        else #2. throw frisbe error here
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end

    catch UnexpectedError # 3. catch `materialize` (UnexpectedError object )
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )
    return response
end

function compareVector(v=[1,2,3], arr = [2, 1, 3, 4])

    try #1. we call this function when we'd like to compare index α with index β of a Vector array  # do your thing
        response = nothing
        a = v[1]
        b = length(v)-1

        firstContent = Int(findfirst(isequal(a), arr)) #indexOf(first)
        lastContent = Int(findfirst(isequal(b), arr)) #indexOf(last)

        if firstContent > lastContent # correct
            response = @inbounds arr[a], arr[b] = oldschoolSwap!(arr[a], arr[b]) #plain content swap in julia  #swap array contents directly

        elseif firstContent < lastContent #only possible - correct situation (to deal with)
            #Intent: skip
            return response = 1
        else #2. throw frisbe error here
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end

    catch UnexpectedError # 3. catch `materialize` (UnexpectedError object )
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )
    return response
end
res = compareVector(1,2,[2, 1, 3, 4])
typeof(res)# no comparison is needed #ERROR

#if res
count=1
if res === nothing return
    compareVector(1)
end
typeof(res)
v = buildInterval(res)
#we're good to go!

#--------------------------------------------------------
@propagate_inbounds function indexOf(i, v::Vector)
    try
        res = findfirst(isequal(i), v)
        if isnumeric(res) #  isa Number
            return res
        else
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end
    catch UnexpectedError
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())
    end
end

a = [2, 1, 3, 4]
_first = Int(indexOf(1, a)) # copy(a[st])
_last = Int(indexOf(2, a)) # copy(a[ed])

#----------

"""
```
input:

v:  subset vector of indicies (index-space)
a: whole element vector (element-space)
i: the offset from the beginning
```
"""
@propagate_inbounds function replaceVector(v = [1, 2], a = [2, 3, 4, 5]; i = 1) #optimized
  try
        lenV = copy(length(v))
        lenA = copy(length(a))
        j=1
        subsetIndex= j+i
        @inbounds if lenV <= lenA # first assumption
            if subsetIndex <= lenA
            for j ∈ 1:lenV
                a[subsetIndex] = v[subsetIndex]
                #a[2] = v[2] #... lenV
            end

            #elseif lenV > lenA #v can't be larger than original vector a (throws an outofBoundsError )
            else
                throw(error("Unexpected error Occured"))
            end
        end
    catch UnexpectedError
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())
    end
    return a
end

v = [1, 2]
lenV = copy(length(v))
lenA = copy(length(a))

#------------------
a = [2, 3, 4, 5]
v = [1, 2]

tuple = compareVector( [1, 2], [2, 3, 4, 5])

v = buildInterval(tuple)# pass-in a tuple   # buildInterval(tuple[1],tuple[2])
v = buildInterval(tuple)# pass-in a tuple   # buildInterval(tuple[1],tuple[2])
#v = getindex(::, ::Int64)

a = replaceVector(v, a) #returns correct range
#--------------------------------------------
#Uncoment for Debugging
#=
@debug a = replaceVector(v,a)
@debug α = replaceVector2(v,a)
@debug a == α && a === α #true - doesn't matter which one to pick
=#

#-----------------

res = compareVector(1,2) #(1,2) tuple
typeof(res)#tuple

#if res == 1
#else if res is valid
v = buildInterval(res) #to vector

 if res == 1 #no need to proceed further
    return
 elseif typeof(res) == VecOrMat
    replaceVector(v,a;i)
 end

function responseHandling(res,v,a = [1, 2, 3, 4];i=1) #TODO: check value
    if res == +1 #no need to proceed further
        return
    elseif typeof(res) == VecOrMat
   return      replaceVector(v,a;i)
    end

end
#a = replaceVector2(v,a) #ERROR:

#----------- questionable #builds an arbitrary rational numbers fron bound a to b  #0ld-thinking #check-if-Working
@propagate_inbounds function buildRangeAroundPoint(a, mid, b) #tobe deleted
    #checked #works but unhelpful #building theoretical ranges  won't belp in  sorting
    if a >= 0 && mid >= 0 && b >= 0
        q = []
        @inbounds push!(q, buildInterval(a, mid))
        @inbounds push!(q, buildInterval(mid + 1, b))

        return q
    else
        println(positiveMsg)
    end
end

#---transpose (new) (2) possible functions
#buildInterval: note different type of functions used to map the transpose
@propagate_inbounds transpose(buildInterval(1, 5)) # transpose
@propagate_inbounds buildInterval(1, 5)' #adjoint
#-----
buildRangeAroundPoint(1, 5, 9) #errourneous # i wanna 1+10 -> 5 but careful: 1+10 = 11 (odd) & 10-1 = 9 (odd) 11-#TODO: should not be niether 9 nor 11 it should be 10 (10)

#=what-if scenario :
1.what-if the range is already overcalculated: need to fix, using the higher figure: 1+10 -1 =  10 (even)-> 1 middle (as anticipated)
2. generalize to any range by adding -1
#Either:
r(a = 1, b = 19) = abs(b + a) - 1 #Or
length((1:10)) #perfect! # length: Built-in function understands it, as well #use it #instead of a+B range calculation (of the middle )

# 9:11-+1:11-9:11 #:10 =#
#---------------------
""" Warning: Assumes there exits a rational series of numbers from point a to point b
this function applies above & below on orgirinal vecto a
```input:
α: current lower bound (of Interval)
below:
above:
β: current upper bound (of Interval)
interval: original vector interval
`

```output:

```
# Examples:
""" # There's a translation map for a vector space a as well as index space  <---------------
# at this stage: calculated: isEven, no-> calc fracMiddle ceil above, floor below  then we call this function
    @propagate_inbounds function buildAboveSoBelow(below, above, a=[1,2,3,4]) #Checked #requires arr  a as input #a.k.a divideConquer
    try
        if  below >= 0 && above >= 0 && b >= 0

            tuple = compareVector(below, above, a)
            interval = buildInterval(tuple)
            replaceVector(interval, a)

            #right()
            exploreInterval(1,2,[1,2])

            #left()
            exploreInterval(3,4,[3,4])

        else
            throw(error("input arguments must be positive"))
        end

    catch notPositiveInputError
        @error positiveMsg exception = (notPositiveInputError, catch_backtrace())
    end
    return a
end

function middle(a, b) # b  + a -1
    try
        _sum = copy(sumInterval(a, b)) #    b + a - 1
        isItEven = copy(iseven(_sum)) # TODO: surround by a copy
        mid = copy(_sum / 2)  #  -1  #Idea: precalculate mid here (_sum /2 )
        isWhole = nothing
        below = -1
        above = -1
        if isItEven
            # 1 middle calculate
            mid = Int(mid)
            isWhole = true
            return mid, mid + 1, isWhole
        elseif !isItEven
            # calculate fractionalMid
            below = floor(mid)
            above = ceil(mid)
            isWhole = false
            return below, above, isWhole # the differenece is still 1, only way to discriminate is by using isWhole
        else
            throw(error("Unexpected error occured"))
        end
    catch UnexpectedError
        @error "ERROR: UnexpectedError occured" exception =
            (UnexpectedError, catch_backtrace)
    end
end

function exploreInterval(a=3,b=4,arr=[3,4])
   response = nothing
    _first =  arr[a];
   _last =  arr[b];
    cond =  _first <last
     try #1. we call this function when we'd like to compare index α with index β of a Vector array  # do your thing

        firstContent = Int(findfirst(isequal(a), arr)) #indexOf(first)
        lastContent = Int(findfirst(isequal(b), arr)) #indexOf(last)

        if firstContent > lastContent  # && # correct
            response = @inbounds a[a], a[b] = oldschoolSwap!(a[a], a[b]) #plain content swap in julia  #swap array contents directly

        elseif firstContent < lastContent && #only possible - correct situation (to deal with)
            #Intent: skip
            return@error "UnexpectedValueError:Unexpected value for s" exception = (UnexpectedValueError, catch_backtrace())
        #elseif firstContent == lastContent

    else
        throw(error("Unexpected input arguments"))
    end
    catch UnexpectedError
        @error positiveMsg exception = (UnexpectedError, catch_backtrace())
    end
    #in the end, compare range itself, is there any more space to explore?
   checklastRange(3,4,[3,4])
end

res = exploreInterval(3,4,[3,4]) #bounderror # attempt to access n+1 while having n length vector

function  checklastRange(a=3,b=4,arr=[3,4])
try
    difference = ϟ(a,b)
    if difference == 1 #there's exactly 1 return 0 #halt!
        return 0 # congratulations # you've rearched the end #HALT!

    elseif difference >=1
        #calculate Interval difference difference
        sum = ϟ(a,b)
        sum =  Int( sum //2)
        response =  isEven2(sum)
         HandleResponse(a,b,arr)
      #TODO:   compareVector(a,b,arr)
    else
        throw(error("Unexpected input arguments"))
    end
    catch UnexpectedError
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())
    end #<---------------------------------

end

#function middle()

function divideConquer(a=[4,5,6,7,3]) # <----we are here
    #check length
   l =  copy(length(a)-1 )
   isItEven = copy(isEven2(l))
    evenHandling(isItEven) = isItEven ? middle() : buildAboveSoBelow()
   #call middle
     m = middle(a) #found first middle
    pushfirst!(m, middles)
    #Q. where we go next?
    goright(a) # <----
    goleft(a) #<---

end

#----------the End--------

#= in the same context, we can logmsg(), debug(), info(), warn()

macro logmsg(level, exs...) logmsg_code((@_sourceinfo)..., esc(level), exs...) end
macro debug(exs...) logmsg_code((@_sourceinfo)..., :Debug, exs...) end
macro  info(exs...) logmsg_code((@_sourceinfo)..., :Info,  exs...) end
macro  warn(exs...) logmsg_code((@_sourceinfo)..., :Warn,  exs...) end

=#

#-------------------------------
"""with a vector provided, making things easier with Comparing, on the spot

""" #Vital function  #here
@propagate_inbounds function buildAboveSoBelow(
    α,
    β,
    a = [1, 2, 3, 4, 5],
) #Checked
try
        if  α >= 0 && β >= 0 #at this point: the state variable of above & below is uncertain (cannot be sure )
            #Q. is there a need for comparison?  yes, inside compareVector if no need for change tuple would be empty
            tuple = compareVector(α, β, a) #correct
            interval = buildInterval(tuple)
            length(interval)
            replaceVector(interval, a)# TODO: replace at where ?

        else
            throw(error("Unexpected input arguments"))
        end
    catch UnexpectedError
        @error positiveMsg exception = (UnexpectedError, catch_backtrace())
    end

   # return q =#
        return a

end

#---------------------------------------
#=test: building above & below (should be automatic)=#
# (α, below, above, β, interval=[1,2,3,4,5] ) #TODO: #@code_llvm()
function compareVector((1:10;),a)
end
tuple = compareVector(1, 2, a) #correct
interval = buildInterval(tuple)
 getindex(::Nothing,::Int64)
length(interval)

#TODO: what's next checklastInterval
a = buildAboveSoBelow(1, 5, arr) #should return the array a #ERROR: check arguments  #exhausts mid ranges[2] = 6 , 3 element, left, mid, right # still: left, Right

#----test middle ---

# mid = middle(a, b) # function calls itself!
mid = sumInterval(a, b) #mid Pre-cursor
cond = isEven2(mid)#divisible by 2

if cond #isEven,  got mid, mid +1 #start of next range
    mid=Int(mid//2)
    q = buildRangeAroundPoint(α, mid, β)
    # push!(ranges, buildRangeAroundPoint(a, m
end
#---- middle -----------

#e.g. Example

arr = [1, 22, 34, 44, 55]
#=
for i in 1:length(arr)
    println(popfirst!(arr))
end
=#
actualSize = length(arr) + 1  #length = size - 1 #5

#=
 if actualSize == 1
#    return
elseif actualSize >= 1 # still values to choose from
#    evinity = isEven(actualSize)
    #sumInterval(arr[])
    if !evinity #!isEven(actualSize) # false

    #end
    elseif evinity

    end
=#
isEven(actualSize)
#---

q = []
#-------------------------
""" important function for divideConquer: the terminal condition where the distance of ranges become 1  """ # last (first) function of a divideConquer cycle
function isToNotDivide(range) #last function to execute # bounds error
    cond = abs(length(range)[2] - (range)[1])  # upkoppa #condition Not Divide
    # response = nothing
    if cond == 1 # check current range # if 1 (nothing more to explore )
        #To Not Divide this range, get out of program
        return true

    else#if cond > 1 # there are still ranges in-between  (to explore)
        # To Divide
        return false
        #else
    end

end
isToNotDivide

#---testing length properties
length(6:10) - length(3:7) === length(6:10) - length(5:9)
length(6:10) - length(6:10)
#they  not equal
#-------------- range op (-)
6:10-3:6 #6:10 - 3:10  = 6:7:6 #enigma #can it be go around 6 7 times (the only explination that makes sense)
[3]
[2]
(1:2)[1]

isToNotDivide(buildInterval(1:2;)) #bounds error
#isToNotDivide(1:4)

"""calls divideConquer - as  at the end there is only One""" # completes the infinite dragon
function calldivideConquer(a, Interval) # contains arr #passes
    divideConquer(a) #, Interval[1], length(Interval[1]))
    divideConquer(a) #, Interval[2], length(Interval[2]))

end
#the absurd call of the  divide (divine)
#finish !
#--------------------------------------------

#----------------------------------------------
#example
arr = [1, 2, 3, 4, 5]
Length = copy(length(arr))
actual = arr[1] + Length #first(arr) + length(arr) #emulates an actual errourneous input
a = indexOf(arr[1], arr) # 1
b = actualSize - 1; #fudging the number
fractionalMid = length(arr) / 2 # 5/2 = 2.5 #does not exist go to the nearest neighbor #mid = 2.5
above = Int(ceil(fractionalMid)) # 3
below = Int(floor(fractionalMid)) # 2
 buildAboveSoBelow(below,above,arr)
function middle(a=[1,3,4,5,6])
    ϟ(a[1],length(a)-1)
end
#------------------------------------------------

#----------------------
#middle(arr) # find middle #q old thinking detedted

compareVector(1,2,arr)
a = middle(arr)#1, 10) # true
typeof(a)
l = length(a) # if length == 3 : Right middle, Left situation #examine each on its own
#if subrange == 1 return
length(q[1])
#middle(q[1])
 q[2]
 q[3] #finished middle

#goright(q[3]) # right is the higher side  #TODO:

#---Datas fields -----
ranges = [] # all we want is middles

#--- test: buildRangeAroundPoint
@btime resMid = buildRangeAroundPoint(1, 5, 10) # 104.920 ns - 105.297 ns  (4 allocations: 320 bytes)

# ranges = buildMiddle(1, 2, 4) # error brings up numbers not in the list (as if its reading off a theoretical array found only in the head )

typeof(ranges)
#range = popfirst!(ranges)
#end buildMiddle -----

#--- test middle(a,b) ---

#--- test buildAboveSoBelow----
#--- test - easy :
a = [1, 2, 3, 4, 5]
res = buildAboveSoBelow(1, 3, a) # a <below < above < b

#fine
res = buildAboveSoBelow(1, -1, a) #TODO: check input arguments are only positive
typeof(res) #nothing, rather than an errourneous outcome
#TODO: responseHandling

#up until negative - 1 : Interval is omitted, returns an empty Int64[] Interval , else return other positive Interval # fine
#---end

#---
#check n/2 is a whole number
#--- Left Half
#-- right Half
#---------------------------------------------------

function goleft(a = [1, 2, 3, 4], α = 1, β = length(arr) - 1)
    res = compareVector(α, β ,a) # TODO: middle(arr, a, b) #or # TODO: middle(a, b, arr)
    if res == +1 #no need to proceed further
        return 0
    elseif typeof(res) == VecOrMat
        replaceVector()
    end
end

@benchmark @check_args 1+10^6
@benchmark 1+10^6
#= #UncomentMe
function goright(a = [1, 2, 3, 4], α = 1, β = length(arr) - 1)  #mid + 1, b = length(arr))
    response = nothing
    try
        if res == 1     # ?   return 0 # : replaceVector() #uncommentMe  #no need to proceed further
        #return 0
          response = HandleResponse()
          return response
        else
        throw(error("Unexpected input arguments"))
        end
        elseif typeof(res) == VecOrMat
        replaceVector()
        end
    catch UnexpectedError
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end  #finaly return
    return response
end =#

#=
function HandleResponse(α, β ,a) # todo the hardwork
    res = compareVector(α, β ,a)
    if res == 1     # ?   return 0 # : replaceVector() #uncommentMe  #no need to proceed further

    @check_args lenA = copy(length(arr))
    #size = lenA - 1 # ture last max possible reachable range, in an Array
    isEven(lenA-1) #either even or not
    #middle here:
    middle = Int(lenA // 2)

end =#
#= the end =#
#function definition: formula: b + a - 1  - sumInterval
#recap
# checkdifference(a,b) =  # if interval length equal to 1
sumInterval(a,b) = 0 < b && 0 < a ? b + a - 1 : -1 ; #TODO: '-1' Handling
isEven2(a,b) #TODO: evenHandling - now true ? middle(a,b) : buildAboveSoBelow(a,b)
isEven2(1,5) #true
function middle2(a,b)
 if 0 < b && 0 < a #positive argument @
    Int(sumInterval(a,b) //2)
 else
     -1 #it's safe to throw an error
 end
end
#---
@propogate_inbounds function buildAboveSoBelow(a=1,b=5) # if false
#get exact value
end
#---------start here
res = sumInterval(1,5)
function buildAboveSoBelow(num = 5) # if false
#mid gets exact value
    mid = num /2
    below  =  Int( floor(mid) )
    above = Int( ceil(mid) )
    return below, above
end
below, above =  buildAboveSoBelow(5)
a= 1; b = 5;#say #-now build corresponding 2 range(tuples )
range1 = (a,below)
typeof(range1)
range1 = collect(range1)
range2  = (below+1, above) #normal way of thinking havent covered this case,yet
range2 = collect(range2) #where range is range of inicies
#possible errourneous (without arr )
function compareVector(range=range2) #a = 1, b = 2, arr = [2, 1, 3, 4])
    response = nothing
    try #1. we call this function when we'd like to compare index α with index β of a Vector array  # do your thing
        firstContent = Int(findfirst(isequal(range2[1]), range)) #indexOf(first)
        lastContent = Int(findfirst(isequal(range2[2]), range)) #indexOf(last)

        if firstContent > lastContent # correct
            response = @inbounds arr[a], arr[b] = oldschoolSwap!(arr[a], arr[b]) #plain content swap in julia  #swap array contents directly

        elseif firstContent < lastContent #only possible - correct situation (to deal with)
            #Intent: skip
            return response = 1
        else #2. throw frisbe error here
            throw(error("Unexpected Error")) # 2. throw(error(ExceptionError))
        end

    catch UnexpectedError # 3. catch `materialize` (UnexpectedError object )
        @error UnexpMsg exception = (UnexpectedError, catch_backtrace())   # define Exception here, passing arguments 1. positiveError object, 2. call catch_backtrace() (to catch it)
    end #ends try - finally afterthat return whatever correct value you've been working on  (if not already )
    return response
end

#check distance, first  (we have 2 ranges range1, range2 )

if  range1[2] - range1[1]  == 1
    a = compareVector(range1)
end
a
a = compareVector(range1)

firstContent = Int(findfirst(isequal(range1[1]), range1)) #indexOf(first)
lastContent = Int(findfirst(isequal(range1[2]), range1)) #indexOf(last)

    if firstContent > lastContent # correct
        response = @inbounds arr[a], arr[b] = oldschoolSwap!(arr[a], arr[b]) #plain content swap in julia  #swap array contents directly
    end
response
#assume first range is 1, 2
#must be handling for it , alerting that it's tur , first < last , returns 0 for finish correctly
#special unhandeled case: distance == 0 where 3 - 3 = 0
range2[2]
range2[1]
#arr 1,2,3
#last range 3,3 returns 0  #singleton = 3 #scalar value - not a range ! #must insertat sth (or pushlast ) ?
a =[1,2,3]
 #cannot do that! \
#input can be a ramge if dist(range)>=1  or if 0 , return that scalar element, at its index i.e. a[3]
a[3]
# there must be a functionality keeps track of middles ,informs that oh, we're out of middles , we're finished, we can exit safely wihth no errors
#middle to raPnge functinality as well

#--------------------------------------------------------------------------

max_stars_repo_name: adamwillisMastery/DeepLearner

max_issues_repo_name: adamwillisMastery/DeepLearner

max_forks_repo_name: adamwillisXanax/DeepLearner