text:
#Utility to convert between xyz and internal (Z-matrix) coordinates
module Converter
"""
Contains functions that convert between xyz and z-matrix formats to serve as
input in the LIIC routines.
"""

using LinearAlgebra
using Formatting

export xyz_to_zmat, zmat_to_xyz
export Molecule, ZMatrix

struct Molecule
    Atoms::Vector{String}
    Coord::Array{Float64,2}
    Number::Int64
end

struct ZMatrix
    Atoms::Vector{String}
    IntVars::Vector{Float64}
    VarNames::Vector{String}
    Number::Int64
end

function import_molecule(xyz_file::String)
    io = open(xyz_file)
    natoms = parse(Int,readline(io))
    readline(io);
    reg=r"([a-zA-Z]+)\s+(-?\d+.?\d*)\s+(-?\d+.?\d*)\s+(-?\d+.?\d*)"
    atoms=Vector{String}()
    atom_coords=Array{Float64,2}(undef,natoms,3)
    for num in 1:natoms
        m=match(reg,readline(io))
        push!(atoms,string(m.captures[1]::SubString))
        atom_coords[num,:]=[parse(Float64,m.captures[2]),parse(Float64,m.captures[3]),parse(Float64,m.captures[4])]
    end
    close(io)
    return Molecule(atoms,atom_coords,natoms)
end

function bond(a1::Vector{Float64},a2::Vector{Float64})
    return norm(a1-a2)::Float64
end

function angle(a1::Vector{Float64},a2::Vector{Float64},a3::Vector{Float64})
    bond1 = a3 - a2
    bond2 = a2 - a1
    sintheta = norm(cross(bond1,bond2))
    costheta = dot(bond1,bond2)
    theta = atan(sintheta,costheta)
    angle = 180.0 * theta/pi
    # angle = acosd(dot(bond1,bond2)/(norm(bond1)*norm(bond2)))
    return angle::Float64
end

function dihedral(a1::Vector{Float64},a2::Vector{Float64},a3::Vector{Float64},a4::Vector{Float64})
    bond1 = a4-a3
    bond2 = a3-a2
    bond3 = a2-a1
    plane1 = cross(bond1,bond2)
    plane2 = cross(bond2,bond3)
    dihedral = acosd(dot(plane1,plane2)/(norm(plane1)*norm(plane2)))
    if dot(cross(plane1,plane2),bond2) < 0
        dihedral = -dihedral
    end
    return dihedral::Float64
end

function write_zmat(zmat::ZMatrix,withvars::Bool=true)
    if withvars==true
        natoms = zmat.Number
        for atom in 1:natoms
            if atom < 2
                printfmtln("{:3s}",zmat.Atoms[atom])
            elseif atom < 3
                printfmtln("{:3s} {:>4d}  {:11s}",zmat.Atoms[atom],(atom-1),zmat.VarNames[atom-1])
            elseif atom < 4
                printfmtln("{:3s} {:>4d}  {:11s} {:>4d}  {:11s}",zmat.Atoms[atom],(atom-1),
                zmat.VarNames[atom-1],(atom-2),zmat.VarNames[atom])
            else
                printfmtln("{:3s} {:>4d}  {:11s} {:>4d}  {:11s} {:>4d}  {:11s}",zmat.Atoms[atom],
                (atom-1),zmat.VarNames[3*atom-8],(atom-2),zmat.VarNames[3*atom-7],(atom-3),zmat.VarNames[3*atom-6])
            end
        end
        println("  Variables:")
        for var in 1:length(zmat.IntVars)
            printfmtln("{:3s} {:>11.8f}",zmat.VarNames[var],zmat.IntVars[var])
        end
    else
        natoms = length(zmat.Atoms)
        for atom in 1:natoms
            if atom < 2
                printfmtln("{:3s}",zmat.Atoms[atom])
            elseif atom < 3
                printfmtln("{:3s} {:>4d}  {:>11.8f}",zmat.Atoms[atom],(atom-1),zmat.IntVars[atom-1])
            elseif atom < 4
                printfmtln("{:3s} {:>4d}  {:>11.8f} {:>4d}  {:>11.8f}",zmat.Atoms[atom],(atom-1),
                zmat.IntVars[atom-1],(atom-2),zmat.IntVars[atom])
            else
                printfmtln("{:3s} {:>4d}  {:>11.8f} {:>4d}  {:>11.8f} {:>4d}  {:>11.8f}",zmat.Atoms[atom],
                (atom-1),zmat.IntVars[3*atom-8],(atom-2),zmat.IntVars[3*atom-7],(atom-3),zmat.IntVars[3*atom-6])
            end
        end
    end
end

function write_xyz(mol::Molecule,io=nothing)
    natoms = mol.Number
    if io !== nothing
        printfmtln(io,"{:>4d}",natoms)
        println(io)
        for atom in 1:natoms
            printfmtln(io,"{:3s}   {:>11.8f}   {:>11.8f}    {:>11.8f}",mol.Atoms[atom],mol.Coord[atom,1],mol.Coord[atom,2],mol.Coord[atom,3])
        end
    elseif io === nothing
        printfmtln("{:>4d}",natoms)
        println()
        for atom in 1:natoms
            printfmtln("{:3s}   {:>11.8f}   {:>11.8f}    {:>11.8f}",mol.Atoms[atom],mol.Coord[atom,1],mol.Coord[atom,2],mol.Coord[atom,3])
        end
    else
        println("Error in write_xyz")
    end
end

function xyz_to_zmat(molecule)
    natoms     = molecule.Number
    atoms      = molecule.Atoms
    intvars      = Vector{Float64}()
    varnames   = Vector{String}()
    for atom in 1:natoms
        if atom > 1
            bl     = bond(molecule.Coord[atom-1,:],molecule.Coord[atom,:])
            push!(intvars,bl)
            push!(varnames,"R$atom")
        end
        if atom > 2
            agl  = angle(molecule.Coord[atom-2,:],molecule.Coord[atom-1,:],molecule.Coord[atom,:])
            push!(intvars,agl)
            push!(varnames,"A$atom")
        end
        if atom > 3
            dhl  = dihedral(molecule.Coord[atom-3,:],molecule.Coord[atom-2,:],molecule.Coord[atom-1,:],molecule.Coord[atom,:])
            push!(intvars,dhl)
            push!(varnames,"D$atom")
        end
    end
    zmatrix    = ZMatrix(atoms,intvars,varnames,natoms)
    return zmatrix
end

function rotation_matrix(axis,angle)
    # Rotation about the origin in 3 dimensions using the Euler-Rodrigues formula
    axis = -axis / norm(axis)
    a = cosd(angle*0.5)
    (b,c,d) = axis * sind(angle*0.5)
    rotM = [a*a+b*b-c*c-d*d 2*(b*c-a*d) 2*(b*d+a*c);
    2*(b*c+a*d) a*a+c*c-b*b-d*d 2*(c*d-a*b);
    2*(b*d-a*c) 2*(c*d+a*b) a*a+d*d-b*b-c*c]
    return rotM::Matrix{Float64}
end

function zmat_to_xyz(zmat::ZMatrix)
    natoms=zmat.Number
    atom_coords = Array{Float64,2}(undef,natoms,3)

    if natoms > 0
        atom_coords[1,:]=[0.0,0.0,0.0]

        if natoms > 1
            r2 = zmat.IntVars[1]
            atom_coords[2,:]=[r2,0.0,0.0]

            if natoms > 2
                vec = atom_coords[2,:] - atom_coords[1,:]
                r3 = zmat.IntVars[2]
                a3 = zmat.IntVars[3]
                vec = r3 * vec / norm(vec)
                atom_coords[3,:] = atom_coords[2,:] + (rotation_matrix([0, 0, 1], a3) * vec)

                if natoms > 3
                    for atom in 4:natoms
                        rn = zmat.IntVars[3*atom-8]
                        an = zmat.IntVars[3*atom-7]
                        dn = zmat.IntVars[3*atom-6]

                        # ox = rn * cosd(an)
                        # oy = rn * cosd(dn) * sind(an)
                        # oz = rn * sind(an) * sind(dn)

                        # Coordinates of the last three atoms
                        am1 = atom_coords[atom-1,:]
                        am2 = atom_coords[atom-2,:]
                        am3 = atom_coords[atom-3,:]

                        # Vectors between them
                        vec1 = am1 - am2
                        vec2 = am2 - am3
                        # Plane normal defined by these vectors
                        nm = cross(vec1,vec2)

                        # Vector of length rn pointing along vec1
                        vecrn = rn * vec1 / norm(vec1)

                        # Rotate vecrn by angle an around the normal
                        vecan = rotation_matrix(nm,an)*vecrn

                        # Rotate vecan around dihedral angle
                        vecdn = rotation_matrix(vec1,dn)*vecan

                        # Add this vector to previous coordinate
                        atom_coords[atom,:] = am1 + vecdn
                    end
                end
            end
        end
    end

    mol=Molecule(zmat.Atoms,atom_coords,natoms)

    return mol::Molecule
end

end

max_stars_repo_name: dasha-shchep/LIIC

max_issues_repo_name: dasha-shchep/LIIC

max_forks_repo_name: dasha-shchep/LIIC


text:
include("utils.jl")

filename = []
title_names = []
# 1
push!(filename, "data/avg_earth_hs_he_12_hp_4_ve_12_vp_4_roefanov.jld2")
push!(title_names, "truth")

# 2
push!(filename, "data/avg_earth_hs_he_15_hp_1_ve_15_vp_1_roefanov.jld2")
push!(title_names, "pₕ = 1 , pᵥ = 1")

# 3
push!(filename, "data/avg_earth_hs_he_10_hp_2_ve_10_vp_2_roefanov.jld2")
push!(title_names, "pₕ = 2 , pᵥ = 2")

# 4
push!(filename, "data/avg_earth_hs_he_6_hp_4_ve_15_vp_1_roefanov.jld2")
push!(title_names, "pₕ = 4 , pᵥ = 1")

# 5
push!(filename, "data/avg_earth_hs_he_6_hp_4_ve_6_vp_4_roefanov.jld2")
push!(title_names, "pₕ = 4 , pᵥ = 4")

# 6
push!(filename, "data/avg_earth_hs_he_5_hp_5_ve_5_vp_5_roefanov.jld2")
push!(title_names, "pₕ = 5, pᵥ = 5")

# 7
push!(filename, "data/avg_earth_hs_he_5_hp_5_ve_5_vp_5_roefanov.jld2")
push!(title_names, "pₕ = 5, pᵥ = 5")

colorrange = (0, 360) # modify to extrema of "truth"
contour_levels = [0, 40, 80, 160, 240, 320]
fig = Figure(resolution = (1700+600, 1000+400))
add_label = true

for i in 1:6
    ii = (i-1)%3 + 1 # +1 on why 1 based indexing is wrong
    jj = (i-1)÷3 + 1 # +1 on why 1 based indexing is wrong
    jl_file = jldopen(filename[i], "r+")

    s_string = "uu"
    slice_zonal1, s_string = eddy_variance(s_string, jl_file)

    s_string = "vv"
    slice_zonal2, s_string = eddy_variance(s_string, jl_file)

    slice_zonal = 0.5 .* (slice_zonal1 + slice_zonal2)

    s_string = title_names[i]
    λ, ϕ, r, p_coord = grab_grid(jl_file)
    ax = fig[jj,ii] = Axis(fig, title = s_string, titlesize = 40)
    if ii !==1
    contour_heatmap!(ax, ϕ, p_coord, slice_zonal, contour_levels, colorrange,
                     add_labels = add_label, colormap = :thermometer)
    else
        contour_heatmap!(ax, ϕ, p_coord, slice_zonal, contour_levels, colorrange,
                     add_labels = add_label, colormap = :thermometer,
                     random_seed = 10)
    end
    if ii !== 1
        hideydecorations!(ax, grid = false)
    end
end

ext: jl

max_stars_repo_name: sandreza/HeldSuarezVisualizationScripts

max_issues_repo_name: sandreza/HeldSuarezVisualizationScripts

max_forks_repo_name: sandreza/HeldSuarezVisualizationScripts


text:
using Test
using ConvOptDL.Utils

@testset "Utils Test" begin
    A = hcat([ones(4)*i for i in 1:5]...)
    B = hcat([ones(3)*i for i in 1:5]...)
    refC = ones(4, 3, 5)
    for i in 1:5
        refC[:, :, i] .*= i^2
    end
    @test batch_outer(A, B) == refC
end


text:
using StaticArrays
using LinearAlgebra
import Base.-, Base.+, Base.*, Base./
import Base.getindex
import PyPlot

const Scalar = Float64;
const Vector = SVector{2,Float64};
const Label  = Int;

const ScalarList = Array{Scalar,1};
const VectorList = Array{Vector,1};
const LabelList  = Array{Label,1};
const LabelListList = Array{LabelList,1};

×(u::Vector, v::Vector) = u[1]*v[2] - u[2]*v[1];

struct PatchInfo
    name  :: String
    range :: UnitRange{Label}
end;

const PatchInfoList = Array{PatchInfo,1};

struct Mesh
    point      :: VectorList     # Souradnice vrcholu site
    centre     :: VectorList     # Souradnice stredu bunek
    volume     :: ScalarList     # Velikosti (objemy) bunek
    surface    :: VectorList     # Normalovy vektor na stenu
    facecentre :: VectorList     # Souradnice stredu strany
    owner      :: LabelList      # Vlastnik steny
    neighbor   :: LabelList      # Soused steny
    patch      :: PatchInfoList  # Seznam okrajovych podminek

    cell2point :: LabelListList  # Indexy vrcholu bunek
    face2point :: LabelListList  # Indexy vrcholu sten
end

"""
    Mesh(point, owner, neighbor, patch, cell2point, face2point)

    Vytvori sit ze zadanych souradnic vrcholu a informaci o konektivite
"""
function Mesh(point::VectorList, owner::LabelList, neighbor::LabelList, patch::PatchInfoList, cell2point::LabelListList, face2point::LabelListList)
    m = Mesh(point, VectorList(), ScalarList(), VectorList(), VectorList(), owner, neighbor, patch, cell2point, face2point)
    update!(m)
    return m
end

"""
    update!(mesh::Mesh)

    Vypocte objemy bunek, souradnice stredu bunek, normalove vektory na steny a stredy sten.
"""
function update!(mesh::Mesh)

    # Recalculate cell centres & volumes
    ncells = size(mesh.cell2point,1)
    empty!(mesh.centre)
    empty!(mesh.volume)
    for c=1:ncells
        pts = mesh.point[mesh.cell2point[c]]

        v = Scalar(0)
        x = Vector(0,0)
        p1 = pts[end]
        for i=1:size(pts,1)-2
            p2 = pts[i]
            p3 = pts[i+1]
            vv = (p2 - p1) × (p3 - p1) / 2
            v += vv
            x += vv*(p1 + p2 + p3)/3
        end
        push!(mesh.volume, v)
        push!(mesh.centre, x/v)
    end

    # Recalculate face sizes & positions
    empty!(mesh.surface)
    empty!(mesh.facecentre)
    nfaces = size(mesh.face2point,1)
    for f=1:nfaces
        pts = mesh.point[mesh.face2point[f]]
        push!(mesh.surface, Vector(pts[2][2]-pts[1][2], pts[1][1]-pts[2][1]))
        push!(mesh.facecentre, (pts[1]+pts[2])/2)
    end
end

"""
    cartesian_mesh(nx, ny)

    Vytvori kartezskou sit o nx*ny bunkach v jednotkovem ctverci. Okrajove podminky pojmenuje "left", "right", "top" a "bottom"
"""
function cartesian_mesh(nx, ny)
    Δx, Δy = 1.0/nx, 1.0/ny

    pid(i,j) = i + j*(nx+1) + 1
    cid(i,j) = i + (j-1)*nx

    # Mesh points
    point = VectorList()
    for j=0:ny, i=0:nx
        push!(point, Vector(i*Δx, j*Δy))
    end

    # Mesh cells
    cell2point = LabelListList()
    for j=1:ny, i=1:nx
        pts = [pid(i-1,j-1), pid(i,j-1), pid(i,j), pid(i-1,j)] |> LabelList
        push!(cell2point, pts)
    end

    # Mesh faces
    owner = LabelList()
    neighbor = LabelList()
    face2point = LabelListList()

    # Internal faces
    for j=1:ny, i=1:nx-1
        push!(owner, cid(i,j))
        push!(neighbor, cid(i+1,j))
        push!(face2point, [pid(i,j-1), pid(i,j)])
    end

    for j=1:ny-1, i=1:nx
        push!(owner, cid(i,j))
        push!(neighbor, cid(i,j+1))
        push!(face2point, [pid(i,j), pid(i-1,j)])
    end

    # Boundary patches
    # - bottom
    patch = PatchInfoList()
    j=0
    start = size(owner,1)+1
    for i=1:nx
        push!(owner, cid(i,j+1))
        push!(face2point, [pid(i-1,j), pid(i,j)])
    end
    push!(patch, PatchInfo("bottom", UnitRange(start, size(owner,1))))

    # - right
    i=nx
    start = size(owner,1)+1
    for j=1:ny
        push!(owner, cid(i,j))
        push!(face2point, [pid(i,j-1), pid(i,j)])
    end
    push!(patch, PatchInfo("right", UnitRange(start, size(owner,1))))

    # - top
    j=ny
    start = size(owner,1)+1
    for i=1:nx
        push!(owner, cid(i,j))
        push!(face2point, [pid(i,j), pid(i-1,j)])
    end
    push!(patch, PatchInfo("top", UnitRange(start, size(owner,1))))

    # - left
    i=0
    start = size(owner,1)+1
    for j=1:ny
        push!(owner, cid(i+1,j))
        push!(face2point, [pid(i,j), pid(i,j-1)])
    end
    push!(patch, PatchInfo("left", UnitRange(start, size(owner,1))))

    return Mesh(point, owner, neighbor, patch, cell2point, face2point)
end

cells(m::Mesh) = UnitRange(1,length(m.volume))

internal_faces(m::Mesh) = UnitRange(1,length(m.neighbor))

all_faces(m::Mesh) = UnitRange(1,length(m.owner))

boundary_patches(m::Mesh) = UnitRange(1,length(m.patch))

patch_names(m::Mesh) = [ m.patch[p].name for p in boundary_patches(m)]

patch_by_name(m::Mesh, name::String) = findfirst(x->x==name, patch_names(m))

patch_faces(m::Mesh, patch::Label) = m.patch[patch].range

function plot_mesh(mesh::Mesh, style="-k")
    for f in all_faces(mesh)
        pts = mesh.face2point[f]
        x = [mesh.point[p][1] for p in pts]
        y = [mesh.point[p][2] for p in pts]
        PyPlot.plot(x, y, style)
    end
end

function triangles(m::Mesh)
    tri = Array{Array{Int,1},1}()
    for c in cells(m)
        pids = m.cell2point[c]
        for i=3:length(pids)
            push!(tri, [pids[1], pids[i-1], pids[i]])
        end
    end
    return tri
end

max_stars_repo_name: furstj/NSP

max_issues_repo_name: furstj/NSP

max_forks_repo_name: furstj/NSP


text:
"""
OptimizeMe is a demonstration module used in illustrating how to improve code and generate effective `precompile` directives.
It has deliberate weaknesses in its design, and the analysis of these weaknesses via `@snoopi_deep` is discussed
in the documentation.
"""
module OptimizeMe

struct Container{T}
    value::T
end

function lotsa_containers()
    list = [1, 0x01, 0xffff, 2.0f0, 'a', [0], ("key", 42)]
    cs = Container.(list)
    println("lotsa containers:")
    display(cs)
end

concat_string(c1::Container, c2::Container) = string(c1.value) * ' ' * string(c2.value)

function contain_concrete(item1, item2)
    c1 = Container(item1)
    c2 = Container(item2)
    return concat_string(c1, c2)
end

function contain_list(list)
    cs = Container.(list)
    return concat_string(cs...)
end

struct Object
    x::Int
end

function makeobjects()
    xs = [1:5; 7]
    return Object.(xs)
end

function main()
    lotsa_containers()
    println(contain_concrete(3.14, "is great"))
    list = [2.718, "is jealous"]
    println(contain_list(list))
    display(makeobjects())
end

end


text:
@testset "Test 2D" begin

    p = [0. 0. 1. 1. 0.5;  0. 1. 0. 1. 0.5] # nodes
    u = [0.0; 0.0; 0.0; 0.0; 1.0]   # function values in nodes
    u2 = [0.0; 0.0; 0.0; 0.0; 2.0]  # second function values in nodes
    t = [0.5 0.5 0.499999; 0.5 0.499999 0.5]  # evaluation points

    dp = [0.5 0.5; 0.5 0.5]
    es = [1.0 0.0; 0.0 1.0]
    du = [0.; 100000.0]

    @testset "Test 2D-RK_H0 kernel" begin
        rk = RK_H0(0.001)
        s = interpolate(p, u, rk)
        σ = evaluate(s, t)
        @test σ[1] ≈ u[5]
        @test isapprox(σ[2], u[5], atol = 1e-5)
        @test isapprox(σ[3], u[5], atol = 1e-5)

        σ1 = evaluate_one(s, t[:,1])
        @test σ1[1] ≈ u[5]

        rk = RK_H0()
        s = prepare(p, rk) # prepare spline
        s = construct(s, u) # construct spline
        σ1 = evaluate(s, t) # evaluate spline in points
        @test σ[1] ≈ u[5]
        @test isapprox(σ[2], u[5], atol = 1e-5)
        @test isapprox(σ[3], u[5], atol = 1e-5)

        s = construct(s, u2)
        σ2 = evaluate(s, t)
        @test σ2[1] ≈ u2[5]
        @test isapprox(σ2[2], u2[5], atol = 1e-5)
        @test isapprox(σ2[3], u2[5], atol = 1e-5)

        cond = estimate_cond(s)
        @test cond ≈ 10.0

        eps = get_epsilon(s)
        @test isapprox(eps, 0.93, atol = 1e-2)

        est_eps = estimate_epsilon(p)
        @test isapprox(est_eps, 0.93, atol = 1e-2)
    end

    @testset "Test 2D-RK_H1 kernel" begin
        rk = RK_H1(0.001)
        s = interpolate(p, u, rk)
        cond = estimate_cond(s)
        @test cond ≈ 1.e11

        σ = evaluate(s, t)
        @test isapprox(σ[1], u[5], atol = 1e-5)
        @test isapprox(σ[2], u[5], atol = 1e-5)
        @test isapprox(σ[3], u[5], atol = 1e-5)

        grad = evaluate_gradient(s, p[:,5])
        @test abs(grad[1]) < 1.0e-5 && abs(grad[2]) < 1.0e-5

        rk = RK_H1()
        s = prepare(p, rk) # prepare spline
        cond = estimate_cond(s)
        @test cond ≈ 100.0

        s = construct(s, u) # construct spline

        σ1 = evaluate(s, t) # evaluate spline in points
        @test isapprox(σ[1], u[5], atol = 1e-5)
        @test isapprox(σ[2], u[5], atol = 1e-5)
        @test isapprox(σ[3], u[5], atol = 1e-5)

        s = construct(s, u2)
        σ2 = evaluate(s, t)
        @test isapprox(σ[1], u[5], atol = 1e-5)
        @test isapprox(σ2[2], u2[5], atol = 1e-5)
        @test isapprox(σ2[3], u2[5], atol = 1e-5)

        est_eps = estimate_epsilon(p, dp)
        @test est_eps ≈ 1.94 atol = 1e-2
        est_eps = estimate_epsilon(p, dp, RK_H1())
        @test est_eps ≈ 1.94 atol = 1e-2
        eps = get_epsilon(s)
        @test est_eps ≈ 1.94 atol = 1e-2

###
        eps = 0.0001
        rk = RK_H1(eps)
        s = interpolate(p, u, dp, es, du, rk)
        cond = estimate_cond(s)
        @test cond == 1.0e14

        σ1 = evaluate_one(s, p[:,5])
        @test !isapprox(σ1[1], u[5], atol = 0.1)

# Same test with extended precision
        rk = RK_H1(Double64(eps))
        p = Double64.(p)
        dp = Double64.(dp)
        es = Double64.(es)
        u = Double64.(u)
        du = Double64.(du)
        t = Double64.(t)
        s = prepare(p, dp, es, rk)
        s = construct(s, u, du)
        σ = evaluate(s, p)
        @test all(isapprox.(σ, u, atol = 1e-5))

        q = estimate_accuracy(s)
        @test q ≈ 14
    end

    @testset "Test 2D-RK_H2 kernel" begin
        p = [0. 0. 1. 1. 0.5;  0. 1. 0. 1. 0.5] # nodes
        u = [0.0; 0.0; 0.0; 0.0; 1.0]   # function values in nodes
        u2 = [0.0; 0.0; 0.0; 0.0; 2.0]  # second function values in nodes
        t = [0.5 0.5 0.499999; 0.5 0.499999 0.5]  # evaluation points

        dp = [0.5 0.5; 0.5 0.5]
        es = [1.0 0.0; 0.0 1.0]
        du = [0.; 100000.0]

        rk = RK_H2(0.01)
        s = interpolate(p, u, rk)
        cond = estimate_cond(s)
        @test cond ≈ 1e11

        σ = evaluate(s, t)
        @test isapprox(σ[1], u[5], atol = 1e-5)
        @test isapprox(σ[2], u[5], atol = 1e-5)
        @test isapprox(σ[3], u[5], atol = 1e-5)

        grad = evaluate_gradient(s, p[:,5])
        @test abs(grad[1]) < 1.0e-5 && abs(grad[2]) < 1.0e-5

        rk = RK_H2()
        s = prepare(p, rk) # prepare spline
        cond = estimate_cond(s)
        @test cond ≈ 100.0

        s = construct(s, u) # construct spline

        σ1 = evaluate(s, t) # evaluate spline in points
        @test isapprox(σ[1], u[5], atol = 1e-5)
        @test isapprox(σ[2], u[5], atol = 1e-5)
        @test isapprox(σ[3], u[5], atol = 1e-5)

        s = construct(s, u2)
        σ2 = evaluate(s, t)
        @test isapprox(σ[1], u[5], atol = 1e-5)
        @test isapprox(σ2[2], u2[5], atol = 1e-5)
        @test isapprox(σ2[3], u2[5], atol = 1e-5)
###
        eps = 0.01
        rk = RK_H2(eps)
        s = interpolate(p, u, dp, es, du, rk)
        cond = estimate_cond(s)
        @test cond == 1.0e13

        σ1 = evaluate_one(s, p[:,5])
        @test !isapprox(σ1[1], u[5], atol = 0.1)

        q = estimate_accuracy(s)
        @test q == 0

# Same tests with extended precision
        rk = RK_H2(Double64(eps))
        p = Double64.(p)
        dp = Double64.(dp)
        es = Double64.(es)
        u = Double64.(u)
        du = Double64.(du)
        t = Double64.(t)
        s = prepare(p, dp, es, rk)
        s = construct(s, u, du)

        σ = evaluate(s, p)
        @test all(isapprox.(σ, u, atol = 1e-5))

        q = estimate_accuracy(s)
        @test q ≈ 15
    end
end

@testset "Test 2D-Bis" begin

    p = collect([-1.0 2.0; -1.0 4.0; 3.0 2.0; 3.0 4.0; 1.0 3.0]') # function nodes
    u = [0.0; 0.0; 0.0; 0.0; 1.0]                        # function values in nodes

    t = collect([-1.0 3.0; 0.0 3.0; 1.0 3.0; 2.0 3.0; 3.0 3.0]')  # evaluation points

    @testset "Test 2D-Bis-RK_H0 kernel" begin
        spl = prepare(p, RK_H0(0.001))                   # prepare spline
        c = estimate_cond(spl)                                # get estimation of the problem's Gram matrix condition number
        @test c ≈ 100000.0

        spl = construct(spl, u)                          # construct spline
        vt = [1.0, 3.0]
        σ = evaluate_one(spl, vt)                            # evaluate spline in the node
        @test σ ≈ 1.0

        wt = [0.0, 3.0]
        σ1 = evaluate_one(spl, wt)

        u2 = [0.0; 0.0; 0.0; 0.0; 2.0]
        spl = construct(spl, u2)
        σ2 = evaluate_one(spl, wt)
        @test σ2 ≈ 2.0 * σ1

        spl = interpolate(p, u, RK_H0(0.001))            # prepare and construct spline
        σ = evaluate_one(spl, vt)
        @test σ ≈ 1.0
    end

    @testset "Test 2D-Bis-RK_H1 kernel" begin
        spl = prepare(p, RK_H1(0.001))
        c = estimate_cond(spl)
        @test c ≈ 1.0e11

        spl = construct(spl, u)
        vt = [1.0, 3.0]
        σ = evaluate_one(spl, vt)
        @test σ ≈ 1.0

        wt = [0.0, 3.0]
        σ1 = evaluate_one(spl, wt)

        u2 = [0.0; 0.0; 0.0; 0.0; 2.0]
        spl = construct(spl, u2)
        σ2 = evaluate_one(spl, wt)
        @test σ2 ≈ 2.0 * σ1

        spl = interpolate(p, u, RK_H1(0.001))
        σ = evaluate_one(spl, vt)
        @test σ ≈ 1.0
    end

    @testset "Test 2D-Bis-RK_H2 kernel" begin
        spl = prepare(p, RK_H2(0.001))
        c = estimate_cond(spl)
        @test c ≈ 1.0e15

        spl = construct(spl, u)

        vt = [1.0, 3.0]
        σ = evaluate_one(spl, vt)
        @test σ ≈ 1.0

        wt = [0.0, 3.0]
        σ1 = evaluate_one(spl, wt)
        u2 = [0.0; 0.0; 0.0; 0.0; 2.0]
        spl = construct(spl, u2)
        σ2 = evaluate_one(spl, wt)
        @test σ2 ≈ 2.0 * σ1

        spl = interpolate(p, u, RK_H2(0.001))
        σ = evaluate_one(spl, vt)
        @test σ ≈ 1.0
    end
end

@testset "Test 2D-Grad" begin
        p = collect([0.0 0.0; 1.0 0.0; 0.0 1.0]') # function nodes
        u = [0.0; 0.0; 1.0]                       # function values in nodes
        t = [0.5; 0.5]                            # evaluation points
        p2 = collect([0.0 0.0; 2.0 0.0; 0.0 2.0]') # function nodes
        u2 = [0.0; 0.0; 2.0]                       # function values in nodes
        t2 = [1.0; 1.0]                            # evaluation points
        @testset "Test 2D-Grad-RK_H1 kernel" begin
            spl = interpolate(p, u, RK_H1(0.001))
            grad = evaluate_gradient(spl, t)
            @test abs(grad[1] + 1.0) ≈ 1.0 atol = 1e-2
            @test abs(grad[2]) ≈ 1.0 atol = 1e-2

            spl = interpolate(p2, u2, RK_H1(0.001))
            grad = evaluate_gradient(spl, t2)
            @test abs(grad[1] + 1.0) ≈ 1.0 atol = 1e-2
            @test abs(grad[2]) ≈ 1.0 atol = 1e-2
        end

        @testset "Test 2D-Grad-RK_H2 kernel" begin
            spl = interpolate(p, u, RK_H2(0.001))
            grad = evaluate_gradient(spl, t)
            @test abs(grad[1] + 1.0) ≈ 1.0 atol = 1e-4
            @test abs(grad[2]) ≈ 1.0 atol = 1e-4

            spl = interpolate(p2, u2, RK_H2(0.001))
            grad = evaluate_gradient(spl, t2)
            @test abs(grad[1] + 1.0) ≈ 1.0 atol = 1e-4
            @test abs(grad[2]) ≈ 1.0 atol = 1e-4
        end
end


text:
function compare(yt, wt, yv, Gi, Ai) # using weighted MME
    nt = length(yt)
    nv = length(yv)

    #Ri = diagm(wt)
    XR = X .* wt
    X  = ones(nt)
    Z  = sparse([I zeros(nt, nv)])
    ZR = Z .* wt
    lhs = [ XR'X XR'Z
            ZR'X ZR'Z + Gi ]    # as va = ve = 1, h2=0.5.
    rhs = [ XR'yt
            ZR'yt ]
    bvg = lhs \ rhs
    c1 = cor(yv, bvg[nt+2:end])

    lhs = [ XR'X XR'Z
            ZR'X ZR'Z + Ai ]    # as va = ve = 1, h2=0.5.
    bva = lhs \ rhs
    c2 = cor(yv, bva[nt+2:end])
    println("Weighted\n- using G: $c1;\n- using A: $c2")
    return bvg[nt+2:end], bva[nt+2:end]
end

function compare(yt, yv, Gi, Ai) # not weighed MME
    nt = length(yt)
    nv = length(yv)

    X = ones(nt)
    Z = sparse([I zeros(nt, nv)])
    lhs = [ X'X X'Z
            Z'X Z'Z + Gi ]      # as va = ve = 1, h2=0.5.
    rhs = [ X'yt
            Z'yt ]
    bvg = lhs \ rhs
    c1 = cor(yv, bvg[nt+2:end])

    lhs = [ X'X X'Z
            Z'X Z'Z + Ai ]      # as va = ve = 1, h2=0.5.
    bva = lhs \ rhs
    c2 = cor(yv, bva[nt+2:end])
    println("Not weighted\n- using G: $c1;\n- using A: $c2")
    return bvg[nt+2:end], bva[nt+2:end]
end

"""
    function ccwcv(yt, wt, yv, Gi)
---
This is a results of the meeting on 2021-03-24 for
`c`ross `c`ountry `w`eighted `c`ross `v`alidation.

The conclusion was that to include weigt is a must.
"""
function ccwcv(yt, wt, yv, Gi)
    nt = length(yt)
    nv = length(yv)

    Ri = diagm(wt)
    X = ones(nt)
    Z = sparse([I zeros(nt, nv)])
    lhs = [ X'Ri*X X'Ri*Z
            Z'Ri*X Z'Ri*Z + Gi ] # as va = ve = 1, h2=0.5.
    rhs = [ X'Ri*yt
            Z'Ri*yt ]
    bv = lhs \ rhs
    cc = cor(yv, bv[nt+2:end])
    return cc, bv[nt+2:end]
end

"""
    function cccv(yt, yv, Gi)
---
This is a results of the meeting on 2021-03-24 for
`c`ross `c`ountry `c`ross `v`alidation.

The conclusion was that to include weigt is a must.
"""
function cccv(yt, yv, Gi)
    nt = length(yt)
    nv = length(yv)

    X = ones(nt)
    Z = sparse([I zeros(nt, nv)])
    lhs = [ X'X X'Z
            Z'X Z'Z + Gi ]      # as va = ve = 1, h2=0.5.
    rhs = [ X'yt
            Z'yt ]
    bv = lhs \ rhs
    cc = cor(yv, bv[nt+2:end])
    return cc, bv[nt+2:end]
end

"""
    function make_tv(ts, cs, bv, drp, erc)
---
This function returns training set `:ix` and `:Milk` phenotype data in `ts`,
returns `:ix` and `:Milk` EBV in `cs`.
InterBull ID names, or other real names are not necessary,
as `ix` are indices in `G` and `A`.
"""
function make_tv(ts, cs, bv, drp, erc)
    t1 = filter(row ->row.st == 't', cs)
    df1 = select(innerjoin(t1, ts, on = :ID), :ix, drp, erc)
    dft = filter(row -> row.Milk < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, drp)
    return dft, dfv
end

"""
    function mk_milk_dt(ts, cs, bv)
---
This function returns training set `:ix` and `:Milk` phenotype data in `ts`,
returns `:ix` and `:Milk` EBV in `cs`.
InterBull ID names, or other real names are not necessary,
as `ix` are indices in `G` and `A`.
"""
function mk_milk_dt(ts, cs, bv)
    t1  = filter(row -> row.st == 't', cs)
    df1 = select(innerjoin(t1, ts, on = :ID), :ix, :Milk, :Milk_ERC)
    dft = filter(row -> row.Milk < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, :Milk)
    return dft, dfv
end

"""
    function mk_scs_dt(ts, cs, bv)
---
This function returns training set `:ix`, `:SCS`, `:SCS_ERC` phenotype data in `ts`,
returns `:ix` and `:SCS` EBV in `cs`.
InterBull ID names, or other real names are not necessary,
as `ix` are indices in `G` and `A`.
"""
function mk_scs_dt(ts, cs, bv)
    t1  = filter(row -> row.st == 't', cs)
    df1 = select(innerjoin(t1, ts, on = :ID), :ix, :SCS, :SCS_ERC)
    dft = filter(row -> row.SCS < 999, df1)
    t2  = filter(row -> row.st == 'v', cs)
    dfv = select(innerjoin(t2, bv, on = :ID), :ix, :SCS)
    return dft, dfv
end

"""
    function cv_G(ts, vs, G)
---
Cross validation with only G matrix.
Returns correlation of weighted and not weighted on ERC.
"""
function cv_G(ts, vs, G, trait, erc)
    id = ts[:, :ix]
    yt = ts[:, trait]
    wt = ts[:, erc]
    nt = length(id)

    yv = vs[:, trait]
    append!(id, vs.ix)
    nv = length(yv)

    gi = inv(G[id, id])
    c1, _ = ccwcv(yt, wt, yv, gi)
    c2, _ = cccv(yt, yv, gi)
    return c1, c2
end


text:
"""
    function _rippa(A,a)

Estimate the Leave-One-Out (LOO) errors using rippas method. Complexity of 𝛰(3)
compared to calculating the exact LOO at a cost of 𝛰(4). `A` is the RBF matrix
and `a` is the weights of the RBF.
"""
function _rippa(A, a)
    N = size(A, 1)           #Number of Leave-One-Out (LOO) errors
    e = Matrix(I, N, N)
    E = Array{Float64}(undef,N)

    for k = 1:N
        xₖ = A \ e[:,k]       #xₖ solution to Ax[k] = e[k]   - Solved N times
        E[k] = a[k] / xₖ[k]   #Estimated error for k-th subset
    end

    return E
end

"""
    function RMSErrorLOO(interp,samples::Array{Float64,2},plan;
    cond_max = 1/eps(Float64)/10000, rippa = false)

Calculate the Leave-One-Out RMS error for a interpolation method.
`rippa` can be used to calculate the approximation of the LOO error
for Radial Basis Functions at a cost of 𝛰(3)
(compared to 𝛰(4)). `cond_max` sets the maximum allowed condition number for
matrix `A` used in the RBF calculation.
"""
function RMSErrorLOO(interp, samples, plan, smooth;
    cond_max::Float64=1e6, rippa::Bool=false, rbf_dist_metric = Euclidean())

    #initiate arrays
    N = length(samples)
    E = Array{Float64}(undef,N)       #Error
    A = zeros(Float64, N, N)    #RBF matrix A
    ests = Array{Float64}(undef,N)    #Function estimate based on LOO
    LOOinds = Array{Int}(undef,N - 1, N)
    for i = 1:N
        LOOinds[:,i] = filter(x -> x != i, 1:N) # Get the leave one out sub indices
    end

    RMSErrorLOO!(E, A, ests, LOOinds, interp, samples, plan, smooth;
    cond_max=cond_max, rippa=rippa, rbf_dist_metric = rbf_dist_metric)
end

function RMSErrorLOO!(E, A, ests, LOOinds, interp::U,
 samples, plan::T, smooth; cond_max::Float64=1e6,
  rippa::Bool=false, rbf_dist_metric = Euclidean()) where T <: AbstractArray where U <: AbstractArray

    N = length(samples)

    #perform Leave-One-Out estimation and calculate error
    if rippa
        @assert typeof(interp) <: Vector{ScatteredInterpolation.RadialBasisFunction} "Rippas
              algorithm only available for Radial Basis Functions"

        #interpolation object trained on the entire dataset
        itp, A = interpolate(interp, plan,
            samples, returnRBFmatrix=true, metric = rbf_dist_metric, smooth=smooth)

        #RBF error estimation based on Rippas algorithm
        E = _rippa(A, itp.w)

    else
        for i = 1:N

            if typeof(smooth) <: AbstractVector
                loo_smooth = smooth[LOOinds[:,i]]
            else
                loo_smooth = smooth
            end

            #interpolation object trained on the LOO information
            itp, A = interpolate(interp[LOOinds[:,i]], plan[:,LOOinds[:,i]],
                samples[LOOinds[:,i]], returnRBFmatrix=true, metric = rbf_dist_metric, smooth=loo_smooth)

             #evaluate the interpolation object in the LOO position
            ests[i] = ScatteredInterpolation.evaluate(itp, plan[:,i])[1]

            E[i] = samples[i] - ests[i]
        end
    end

    #check the conditioning of matrix A for RBFs
    isrbf = typeof(interp) <: ScatteredInterpolation.RadialBasisFunction

    if isrbf && (cond(A) > cond_max)
        RMSE = Inf
    else
        RMSE = sqrt(mean(E.^2))
    end
    return RMSE
end

function RMSErrorLOO!(E, A, ests, LOOinds, interp, samples, plan::T, smooth;
 cond_max::Float64=1e6, rippa::Bool=false, rbf_dist_metric = Euclidean()) where T <: AbstractArray

    N = length(samples)

    #perform Leave-One-Out estimation and calculate error
    if rippa
        @assert typeof(interp) <: ScatteredInterpolation.RadialBasisFunction "Rippas algorithm only available for Radial Basis Functions"

        #interpolation object trained on the entire dataset
        itp, A = interpolate(interp, plan,
         samples, returnRBFmatrix=true, metric = rbf_dist_metric, smooth=smooth)

        #RBF error estimation based on Rippas algorithm
        E = _rippa(A, itp.w)

    else
        for i = 1:N

            if typeof(smooth) <: AbstractVector
                loo_smooth = smooth[LOOinds[:,i]]
            else
                loo_smooth = smooth
            end

            #interpolation object trained on the LOO information
            itp, A = interpolate(interp, plan[:,LOOinds[:,i]],
                samples[LOOinds[:,i]], returnRBFmatrix=true, metric = rbf_dist_metric, smooth=loo_smooth)

            #evaluate the interpolation object in the LOO position
            ests[i] = ScatteredInterpolation.evaluate(itp, plan[:,i])[1]

            E[i] = samples[i] - ests[i]
        end
    end

    #check the conditioning of matrix A for RBFs
    isrbf = typeof(interp) <: ScatteredInterpolation.RadialBasisFunction

    if isrbf && (cond(A) > cond_max)
        RMSE = Inf
    else
        RMSE = sqrt(mean(E.^2))
    end
    return RMSE
end

"""
    function scale(old_x::T,new_min,new_max;old_min=minimum(old_x),
    old_max=maximum(old_x)) where T <: Real

Scale a scalar to match a new range.
"""
function scale(old_x::T,new_min,new_max;old_min=minimum(old_x),old_max=maximum(old_x)) where T <: Real

    old_x = (((old_x - old_min) * (new_max - new_min)) / (old_max - old_min)) + new_min
end

"""
    function scale(old_x::Array{T,1},new_min,new_max;old_min=minimum(old_x),
    old_max=maximum(old_x)) where T <: Real

Scale a vector to match a new range.
"""
function scale(old_x::Array{T,1},new_min,new_max;old_min=minimum(old_x),old_max=maximum(old_x)) where T <: Real

    newX = (((old_x .- old_min) .* (new_max .- new_min)) ./ (old_max .- old_min)) .+ new_min
end

"""
    function scale(old_x::Array{T,2},direction::Int,new_min,new_max;
    old_min=minimum(old_x,direction)::Array{T,2},old_max=maximum(old_x,direction)::Array{T,2}) where T <: Real

Scale a 2D matrix to match a new range along the specified `direction`.
"""
function scale(old_x::Array{T,2},direction::Int,new_min,new_max;old_min=minimum(old_x,dims=direction)::Array{T,2},old_max=maximum(old_x,dims=direction)::Array{T,2}) where T <: Real

    newX = mapslices(x -> scale(x,new_min,new_max), old_x, dims=direction)
end

scale(x::Missing,min_val,max_val;old_min,old_max) = missing

"""
    preprocess_point(points,optres;base_scale::Array{Float64,2})

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point(points,optres;base_scale::Array{Float64,2})

    old_min = minimum(base_scale,dims=2)
    old_max = maximum(base_scale,dims=2)

    preprocessed_point = similar(points)

    for i = 1:size(preprocessed_point,1)
        preprocessed_point[i,:] = scale(points[i,:],-1.0*optres.scaling[i],1.0*optres.scaling[i],
        old_min = old_min[i], old_max = old_max[i])
    end
    return preprocessed_point
end

"""
    preprocess_point!(preprocessed_point,points,optres;old_min::Array{Float64,2},old_max::Array{Float64,2})

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point!(preprocessed_point,points,optres;old_min::Array{Float64,2},old_max::Array{Float64,2})

    for i = 1:size(preprocessed_point,1)
        preprocessed_point[i,:] = scale(points[i,:],-1.0*optres.scaling[i],1.0*optres.scaling[i],
        old_min = old_min[i], old_max = old_max[i])
    end
    return preprocessed_point
end

"""
    preprocess_point!(preprocessed_point,points,optres::SurrogateModelOptim.RBFHypers{Bool,U};old_min::Array{Float64,2},old_max::Array{Float64,2}) where U

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point!(preprocessed_point,points,optres::SurrogateModelOptim.RBFHypers{Bool,U};old_min::Array{Float64,2},old_max::Array{Float64,2}) where U
    preprocessed_point = points
end

"""
    preprocess_point(points,optres::SurrogateModelOptim.RBFHypers{Bool,U};base_scale::Array{Float64,2}) where U

Scales the input points in order to evaluate the adaptively scaled RBF correctly.
"""
function preprocess_point(points,optres::SurrogateModelOptim.RBFHypers{Bool,U};base_scale::Array{Float64,2}) where U
    points
end

"""
    surrogate_evaluate(preprocessed_est_point,estimation_point,itp,optres,old_min,old_max)

Processes the input to scale correctly based on the options used.
"""
function surrogate_evaluate(preprocessed_est_point,estimation_point,itp,optres,old_min,old_max)
    preprocessed_est_point = preprocess_point!( preprocessed_est_point,
                                                estimation_point, optres;
                                                old_min=old_min,
                                                old_max=old_max,
                                                )

    # Evaluate the interpolation object
    res = ScatteredInterpolation.evaluate(itp, preprocessed_est_point)[1]
end

"""
    construct_search_range(plan::Array{Float64,2}, variable_kernel_width,
    min_rbf_width, max_rbf_width, variable_dim_scaling,
    min_scale, max_scale, smooth, max_smooth)

Creates array of tuples for the optimization of the RBF hyperparameters.
"""
function construct_search_range(plan::Array{Float64,2}, variable_kernel_width,
    min_rbf_width, max_rbf_width, variable_dim_scaling,
    min_scale, max_scale, smooth, max_smooth)

    n_dims, n_samples = size(plan)
    sr = Array{Tuple{Float64,Float64},1}()

    # Add the kernel width and type search range
    variable_kernel_width ? n_kerns = n_samples : n_kerns = 1
    push!(sr,   create_sr(  (min_range=min_rbf_width, max_range=max_rbf_width, n_times=n_kerns),
    (min_range=0.0, max_range=1.0, n_times=n_kerns))...)

    # Add the dimensional scaling search range
    variable_dim_scaling ? n_dim_scales = n_dims : n_dim_scales = 0
    push!(sr,   create_sr(  (min_range=min_scale, max_range=max_scale, n_times=n_dim_scales))...)

    # Add the ridge regression smoothing search range
    (smooth == :variable) && (n_smooth = n_samples)
    (smooth == :single)   && (n_smooth = 1)
    (smooth == :single_user)   && (n_smooth = 0)
    (smooth == false)     && (n_smooth = 0)
    push!(sr,   create_sr(  (min_range=0.0, max_range=max_smooth, n_times=n_smooth))...)

    return sr
end

"""
    create_sr(vargs::NamedTuple{(:min_range, :max_range, :n_times),Tuple{Float64,Float64,Int64}}...)

Facilitate the creation of array of tuples for the optimization of the RBF hyperparameters.
"""
function create_sr(vargs::NamedTuple{(:min_range, :max_range, :n_times),Tuple{Float64,Float64,Int64}}...)
    sr = Array{Tuple{Float64,Float64},1}()
    for varg in vargs
        for i = 1:varg[3]
            push!(sr,(varg[1],varg[2]))
        end
    end
    return sr
end

"""
    extract_vector_range(vargs::Int64...)

Supply length of subsets contained in a vector and receive a tuple containing
all the ranges needed to extract each subset in the vector.
"""
function extract_vector_range(vargs::Int64...)
    output = Array{Any,1}()

    count = 0
    for (i,varg) in enumerate(vargs)
        if varg == 0
            push!(output, false)
        else
            push!(output, (count+1):(count+varg))
            count += varg
            (length(output[end]) == 1) && (output[end] = output[end][1])
        end
    end
    return Tuple(output)
end

"""
    extract_bboptim_hypers(bboptim_fcall_vector,plan,kerns,
    variable_kernel_width,variable_dim_scaling,
    smooth,smooth_user)

Get the optimized RBF hyperparameters in useful format for further use.
"""
function extract_bboptim_hypers(bboptim_fcall_vector,plan,kerns,
    variable_kernel_width,variable_dim_scaling,
    smooth,smooth_user)

    n_dims, n_samples = size(plan)

    # Kernel width and type length
    variable_kernel_width ? n_kerns = n_samples : n_kerns = 1

    # Dimensional scaling length
    variable_dim_scaling ? n_dim_scales = n_dims : n_dim_scales = 0

    # Ridge regression smoothing length
    (smooth == :variable) && (n_smooth = n_samples)
    (smooth == :single)   && (n_smooth = 1)
    (smooth == :single_user)   && (n_smooth = 0)
    (smooth == false)     && (n_smooth = 0)

    width_inds, kernel_float_inds, scaling_inds, smooth_inds = extract_vector_range(n_kerns,
                                                            n_kerns,
                                                            n_dim_scales,
                                                            n_smooth
                                                            )

    # Arrange the width and smoothing results
    width = bboptim_fcall_vector[width_inds]
    kernel_float = bboptim_fcall_vector[kernel_float_inds]
    !(scaling_inds == false) ? scaling = bboptim_fcall_vector[scaling_inds] : scaling = scaling_inds
    (smooth == false) &&            (smooth = smooth_inds)
    (smooth == :single) &&          (smooth = bboptim_fcall_vector[smooth_inds])
    (smooth == :variable) &&        (smooth = bboptim_fcall_vector[smooth_inds])
    (smooth == :single_user) &&     (smooth = smooth_user)

    # Arrange the RBF kernel result
    kern_ind = round.(Int,scale(kernel_float,1,length(kerns),old_min=0,old_max=1))
    if variable_kernel_width
        kern = Vector{ScatteredInterpolation.RadialBasisFunction}(undef,size(plan,2))
        for i = 1:size(plan,2)
            kern[i] = kerns[kern_ind[i]](width[i])
        end
    elseif !variable_kernel_width
        kern = kerns[kern_ind](width)
    end

    return kern, scaling, smooth
end

"""
    interp_obj(inpt::Vector{Float64}, kerns, samples,
    plan::Array{Float64,2}; rippa::Bool = false,
    variable_kernel_width::Bool = true, variable_dim_scaling::Bool = true,
    smooth = false, cond_max=cond_max, rbf_dist_metric = Distances.Euclidean(),
    smooth_user::Float64 = 0.0)

Objective function for optimization of interpolation kernel function and width.
"""
function interp_obj(inpt::Vector{Float64}, kerns, samples,
        plan::Array{Float64,2}; rippa::Bool = false,
        variable_kernel_width::Bool = true, variable_dim_scaling::Bool = true,
        smooth = false, cond_max=cond_max, rbf_dist_metric = Distances.Euclidean(),
        smooth_user::Float64 = 0.0)

    kern, scaling, smooth = extract_bboptim_hypers( inpt,plan,kerns,variable_kernel_width,
                                                    variable_dim_scaling,smooth,
                                                    smooth_user)

    optres = RBFHypers(kern, scaling, smooth)

    #Preprocess the plan based on the settings used.
    preprocessed_plan = preprocess_point(plan,optres,base_scale=plan)

    #Wrapped in try catch-block to catch failure to solve linear eq. system
    E = try
        E = RMSErrorLOO(kern, samples, preprocessed_plan, smooth; rippa = rippa,
        cond_max = cond_max, rbf_dist_metric = rbf_dist_metric)
    catch
        E = Inf
    end

    return E
end

"""
    rbf_hypers_opt(samples_org::Array{Float64,2}, plan::Array{Float64,2}, options::Options)

Optimization function of Radial Basis Function kernel and width.
"""
function rbf_hypers_opt(samples_org::Array{Float64,2}, plan::Array{Float64,2}, options::Options)

    @unpack rippa, variable_kernel_width, variable_dim_scaling, rbf_opt_method,
            min_rbf_width, max_rbf_width, min_scale, max_scale, cond_max,
            rbf_dist_metric, rbf_opt_gens, kerns, rbf_opt_pop,
            smooth, max_smooth, smooth_user = options

    # Sample order for ScatteredInterpolation
    samples = vec(samples_org)

    # Create the hyperparameter search range based on the input options
    sr = construct_search_range(plan, variable_kernel_width,
                                min_rbf_width, max_rbf_width, variable_dim_scaling,
                                min_scale, max_scale, smooth, max_smooth)

    # RBF hyperparameter objective function
    itp_obj = function (x)
        interp_obj(x,kerns,samples,plan;
                rippa=rippa, variable_kernel_width=variable_kernel_width,
                variable_dim_scaling=variable_dim_scaling, smooth=smooth,
                cond_max=cond_max,rbf_dist_metric=rbf_dist_metric,)
    end

    # Optimize the interpolant hyperparameters
    res = bboptimize(itp_obj;
            Method=rbf_opt_method,SearchRange=sr, MaxFuncEvals=rbf_opt_gens,
            TraceMode=:silent, rbf_dist_metric=rbf_dist_metric,
            TargetFitness = 1e-5, FitnessTolerance = 1e-6,
            PopulationSize = rbf_opt_pop,
            MaxStepsWithoutProgress=rbf_opt_gens,
            MaxNumStepsWithoutFuncEvals=rbf_opt_gens,
            );

    kern, scaling, smooth = extract_bboptim_hypers( res.archive_output.best_candidate,
                                                    plan,kerns,variable_kernel_width,
                                                    variable_dim_scaling,smooth,smooth_user)

    # Return the optimized hyperparameters in the correct type
    return RBFHypers(kern, scaling, smooth)
end


text:
module UKCarbonIntensityData

using DataFrames
using DataFramesMeta
using Dates
using HTTP
using JSON
using Measures
using OrderedCollections
using PlotUtils
using RecipesBase
using Statistics
using TimeZones
using Measures

export get_regional_data, get_carbon_intensity, get_todays_forecast
export get_generation_mix
export todays_plot, carbon_intensity
export AVAILABLE_REGIONS

include("fetch.jl")
include("plot.jl")
include("mappings.jl")
end


text:
using Random

function generate(num)
    numbers = []
    i = 1
    @time while size(numbers)[1] < num
        append!(numbers, i)
        i += 1
    end
    return numbers
end

function exp(nums, exp)
    exp_numbers = []
    @time for item in nums
        append!(exp_numbers, item^exp)
    end
    return exp_numbers
end

function dot_exp(nums, exp)
    return @time nums.^exp
end

function shuffle_nums(nums)
    shuffled_nums = Random.shuffle(nums)
    return shuffled_nums
end

function sort_nums(nums)
    @time sort!(nums, alg = QuickSort)
    return nums
end

function list_multiplication(a, b)
    @time res = a .* b
    return res
end

numbers = generate(5000000)
exp_gen = exp(numbers, 2)
# dot_exp_gen = dot_exp(numbers, 2) # Same result as above, different code
arr1 = shuffle_nums(numbers)
arr2 = shuffle_nums(numbers)
unsorted_nums = list_multiplication(arr1, arr2)
sorted_nums = sort_nums(unsorted_nums)
println(sorted_nums)


text:
struct ArgIterator <: AbstractArgIter
    dict::Dict
    static_args::Dict{String, Any}
    arg_order::Vector{String}
    done::Bool
end

function ArgIterator(dict, static_arg::Dict=Dict(); arg_order=nothing)
    ArgIterator(dict,
                Dict{String, Any}(static_arg),
                arg_order==nothing ? collect(keys(dict)) : arg_order,
                false)
end

set_save_dir!(iter::ArgIterator, path) = iter.static_args["save_dir"] = path

function make_arguments(iter::ArgIterator, state)
    d = Dict{String, Any}()
    for (arg_idx, arg) in enumerate(iter.arg_order)
        if occursin("+", arg) #contains(arg, "+")
            ks = split(arg, "+")
            for (idx, k) ∈ enumerate(ks)
                d[k] = iter.dict[arg][state[2][arg_idx]][idx]
            end
        else
            d[arg] = iter.dict[arg][state[2][arg_idx]]
        end
    end
    merge!(d, iter.static_args)
    d
end

function Base.iterate(iter::ArgIterator)
    state = (1, ones(Int64, length(iter.arg_order)))
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function next_state(iter::ArgIterator, _state)
    state = _state[2]
    n_state = _state[1]

    state[end] += 1

    for (arg_idx, arg) in Iterators.reverse(enumerate(iter.arg_order))
        if arg_idx == 1
            return (n_state+1, state)
        end

        if state[arg_idx] > length(iter.dict[arg])
            state[arg_idx] = 1
            state[arg_idx - 1] += 1
        end
    end
end

function Base.iterate(iter::ArgIterator, state)
    if state[2][1] > length(iter.dict[iter.arg_order[1]])
        return nothing
    end
    arg_list = make_arguments(iter, state)
    return (state[1], arg_list), next_state(iter, state)
end

function Base.length(iter::ArgIterator)
    return *([length(iter.dict[key]) for key in iter.arg_order]...)
end


text:
using PyPlot

function plot_histogram_fit(histo, λ0, λ1, bright_ratio)
    ns = collect(0:(length(histo) - 1))
    model = sum(histo) * bipoisson(ns, λ0, λ1, bright_ratio)

    figure()
    plot(0:length(histo) - 1, histo, "ro")
    plot(ns, model, "-k")
    axvline(count_threshold(λ0, λ1), color="r", linestyle="dotted")
    xlim(-length(histo) / 100.0, length(histo))

    figure()
    relative_histo = histo ./ sum(histo)
    semilogy(0:length(histo) - 1, relative_histo, "ro")
    semilogy(ns, model ./ sum(histo), "-k")
    axvline(count_threshold(λ0, λ1), color="r", linestyle="dotted")
    xlim(-length(histo) / 100.0, length(histo))
end

function plot_01_2d_raw(data, xs, ys; diverging=false, ax=nothing, min=0, max=1)
    if ax == nothing
        ax = gca()
    end

    x_offset = (xs[2] - xs[1]) / 2
    corner_xs = [xs - x_offset; xs[end] + x_offset]

    y_offset = (ys[2] - ys[1]) / 2
    corner_ys = [ys - y_offset; ys[end] + y_offset]

    img = ax[:pcolormesh](corner_xs, corner_ys, data', cmap=(diverging ? "coolwarm" : plasma()), vmin=min, vmax=max)
    ax[:set_xlim](corner_xs[1], corner_xs[end])
    ax[:set_ylim](corner_ys[1], corner_ys[end])
    return img
end

function plot_01_2d(data, xs, ys, label=""; diverging=false, ax=nothing, min=0, max=1)
    img = plot_01_2d_raw(data, xs, ys, diverging=diverging, ax=ax, min=min, max=max)

    cb = colorbar(img, aspect=30)
    cb[:set_label](label)
    return img, cb
end


text:
include("NumFieldOrd/NfOrd.jl")
include("NumFieldOrd/NumFieldOrdIdl.jl")
include("NumFieldOrd/NfRelOrd/NfRelOrd.jl")
include("NumFieldOrd/NfRelOrd/OrdElem.jl")
include("NumFieldOrd/NfRelOrd/Ideal.jl")
include("NumFieldOrd/NfRelOrd/FracIdeal.jl")
include("Map/NfRelOrd.jl")
include("NumFieldOrd/NfRelOrd/ResidueRing.jl")
include("NumFieldOrd/NumFieldOrd.jl")
include("NumFieldOrd/NumFieldOrdElem.jl")


text:
WHITE_BG = 0.05N0f16

function transfertotif(img, filename, r = 0, l = 15)

    @showprogress @sprintf("Saving %s:", filename) for f = 1:size(img, 4)
        img1 = img[:, :, :, f]
        if r == 1
            x, y, z = size(img1)
            for i = x-l*3-5:x-l*3
                for j = y-l*4:y-l*3
                    for h = 1:z
                        img1[i, j, h] = WHITE_BG
                    end
                end
            end
        end
        img_save(img1, pwd(), @sprintf("%s_%d.tif", filename, f))
    end
end


text:
function GD_ellispe(∇f::Function, b, x_k, y_k, ϵ)
    k = 1
    x_vals = [x_k]
    y_vals = [y_k]
    while norm(∇f(x_k,y_k)) > ϵ
        x_k = b*((b-1)/(b+1))^k
        y_k = ((1-b)/(b+1))^k
        k+=1

        push!(x_vals, x_k)
        push!(y_vals, y_k)
    end
    return x_vals, y_vals
end

function Quad_Backtracking(f::Function, ∇f::Function, P, x_k, y_k, α, β, ϵ)

    x_vals = []
    y_vals = []

    ∇fx = ∇f(x_k, y_k)

    for i in 1:200
        if norm(∇fx) < ϵ
            break
        end
        ∇fx = ∇f(x_k, y_k)
        fx = f(x_k, y_k)
        Δxk = P\∇fx
        ∇fxT_Δxk = ∇fx'Δxk
        t = 1
        while f(x_k + t*Δxk[1], y_k + t*Δxk[2]) > fx + α*t*∇fxT_Δxk
            t = β * t
        end
        x_k = x_k + t*Δxk[1]
        y_k = y_k + t*Δxk[2]

        push!(x_vals, x_k)
        push!(y_vals, y_k)
    end
    return x_vals, y_vals
end

function updating_Hessien_Descent(f::Function, ∇f::Function, H::Function, x_k, y_k, l, α, β, ϵ)

    x_vals = []
    y_vals = []

    P = -1 * H(x_k, y_k)
    for i in 1:200
        ∇fx = ∇f(x_k, y_k)
        if norm(∇fx) < ϵ
            break
        end

        if i % l == 0
            P = -1 * H(x_k, y_k)
        end

        fx = f(x_k, y_k)
        Δxk = P\∇fx
        ∇fxT_Δxk = ∇fx'Δxk
        t=1
        while f(x_k + t*Δxk[1], y_k + t*Δxk[2]) > fx + α*t*∇fxT_Δxk
            t = β * t
        end

        x_k = x_k + t*Δxk[1]
        y_k = y_k + t*Δxk[2]

        push!(x_vals, x_k)
        push!(y_vals, y_k)
    end
    return x_vals, y_vals
end

;


text:
using Gen: mh, choicemap, @gen, @trace, uniform, normal, get_args, update, ChoiceMap, NoChange, bernoulli
using LinearAlgebra: norm
using GenDirectionalStats: Rot3

###############################
# random walk on some address #
###############################

@gen function random_walk_proposal(trace, address, width)
    prev_value = trace[address]
    @trace(normal(prev_value, width), address)
end

function random_walk_mh(trace, address, width)
    mh(trace, random_walk_proposal, (address, width))
end

export random_walk_mh

######################################################
# rotate by uniformly random angle around fixed axis #
######################################################

function compose_rotations(
        prev_q::Rotations.Rotation{3,Float64},
        change_q::Rotations.Rotation{3,Float64}, egocentric::Bool)
    # TODO check
    if egocentric
        return prev_q * change_q
    else
        return change_q * prev_q
    end
end

no_argdiffs(args) = map((_) -> NoChange(), args)

@gen function uniform_angle_fixed_axis_proposal(trace, addr, axis, egocentric)
    angle ~ uniform(0, 2 * pi)
    return angle
end

function uniform_angle_fixed_axis_involution(trace, ::ChoiceMap, angle::Real, prop_args; check=false)
    addr, axis, egocentric = prop_args
    prev_rot::Rot3 = trace[addr]
    change = Rotations.AngleAxis(angle, axis[1], axis[2], axis[3])
    new_rot = compose_rotations(prev_rot, change, egocentric)
    args = get_args(trace)
    argdiffs = no_argdiffs(args)
    new_trace, w = update(trace, args, argdiffs, choicemap((addr, new_rot)))
    back_angle = 2 * pi - angle
    backward_choices = choicemap((:angle, back_angle))
    return (new_trace, backward_choices, w)
end

function uniform_angle_fixed_axis_mh(trace, addr, axis; check=false, egocentric=true)
    return mh(
        trace,
        uniform_angle_fixed_axis_proposal,
        (addr, axis, egocentric),
        uniform_angle_fixed_axis_involution;
        check=check)
end

export uniform_angle_fixed_axis_mh

###########################################
# rotate by 180 degrees around fixed axis #
###########################################

@gen function flip_proposal(trace, addr, axis, egocentric)
    return nothing
end

function flip_involution(trace, ::ChoiceMap, ::Nothing, prop_args; check=false)
    addr, axis, egocentric = prop_args
    prev_rot::Rot3 = trace[addr]
    angle = pi
    change = Rotations.AngleAxis(angle, axis[1], axis[2], axis[3])
    new_rot = compose_rotations(prev_rot, change, egocentric)
    args = get_args(trace)
    argdiffs = no_argdiffs(args)
    new_trace, w = update(trace, args, argdiffs, choicemap((addr, new_rot)))
    backward_choices = choicemap()
    return (new_trace, backward_choices, w)
end

function flip_around_fixed_axis_mh(trace, addr, axis; check=false, egocentric=true)
    return mh(trace, flip_proposal, (addr, axis, egocentric), flip_involution; check=check)
end

export flip_around_fixed_axis_mh

##########################################
# random-walk on angle around fixed axis #
##########################################

@gen function small_angle_fixed_axis_proposal(trace, addr, axis, width, egocentric)
    angle_magnitude ~ uniform(0, width)
    direction ~ bernoulli(0.5)
    return (angle_magnitude, direction)
end

function small_angle_fixed_axis_involution(trace, ::ChoiceMap, prop_retval, prop_args; check=false)
    addr, axis, width, egocentric = prop_args
    angle_magnitude, direction = prop_retval
    prev_rot::Rot3 = trace[addr]
    if direction
        angle = angle_magnitude
    else
        angle = -angle_magnitude
    end
    change = Rotations.AngleAxis(angle, axis[1], axis[2], axis[3])
    new_rot = compose_rotations(prev_rot, change, egocentric)
    args = get_args(trace)
    argdiffs = no_argdiffs(args)
    new_trace, w = update(trace, args, argdiffs, choicemap((addr, new_rot)))
    backward_choices = choicemap((:direction, !direction), (:angle_magnitude, angle_magnitude))
    return (new_trace, backward_choices, w)
end

function small_angle_fixed_axis_mh(trace, addr, axis, width; check=false, egocentric=true)
    return mh(
        trace,
        small_angle_fixed_axis_proposal,
        (addr, axis, width, egocentric),
        small_angle_fixed_axis_involution;
        check=check)
end

export small_angle_fixed_axis_mh

#############################################
# random-walk on angle around a random axis #
#############################################

@gen function small_angle_random_axis_proposal(trace, addr, width, egocentric)
    angle = @trace(uniform(-width, width), :angle)
    axis_x = @trace(uniform(-1, 1), :x)
    axis_y = @trace(uniform(-1, 1), :y)
    axis_z = @trace(uniform(-1, 1), :z)
    axis = [axis_x, axis_y, axis_z]
    axis_norm = norm(axis)
    axis = axis / axis_norm
    (angle, axis, axis_norm)
end

function small_angle_random_axis_involution(trace, ::ChoiceMap, prop_retval, prop_args; check=false)
    addr, width, egocentric = prop_args
    angle, axis, axis_norm = prop_retval
    prev_rot::Rot3 = trace[addr]
    change = Rotations.AngleAxis(angle, axis[1], axis[2], axis[3])
    new_rot = compose_rotations(prev_rot, change, egocentric)
    args = get_args(trace)
    argdiffs = no_argdiffs(args)
    new_trace, w = update(trace, args, argdiffs, choicemap((addr, new_rot)))
    backward_choices = choicemap(
        (:angle, angle),
        (:x, -axis[1] * axis_norm),
        (:y, -axis[2] * axis_norm),
        (:z, -axis[3] * axis_norm))
    return (new_trace, backward_choices, w)
end

function small_angle_random_axis_mh(trace, addr, width; check=false, egocentric=true)
    return mh(
        trace,
        small_angle_random_axis_proposal,
        (addr, width, egocentric),
        small_angle_random_axis_involution;
        check=check)
end

export small_angle_random_axis_mh


text:
Using Checkers
using Math
using Test

f(x) = x^2
 @test_formany -10<x<10, f(x) >= 0
 @test_forall x in -10:10, f(x) >= 0

#based from the testing folder in Checkers
	@test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2

#1000 test will be run
#Passes if 100 tests meet the crierea of holding true. Will fail becuase generating 100 pairs of x,y, some will bound to fail condition

#100,000 tests will performed until 100 passing test
	@test_formany ntests = 100 maxtests = 100_000 0 < x < 10,  0 < y < 10,  x < y --> x^2 < y^2
#Same as top but will pass because 100000test will be performed.

#another test the conditional property of Julia

	@test_formany ntests = 100 maxtests = 100_000 0 < x < 7,  8 < y < 80,  x < y --> x^2 < y^2
#Same as top but will pass because 100000test will be performed.

f(x) = sin(x)
	@test_formany -100<x<100, f(x) >= -1

f(x) = sin(x)*cos(x)
 @test_formany -100<x<100, -1< f(x) < 1

#test if tan is right
function tanTest(x)
	sin(x)/cos(x)
	end

	@test_forall -1000<x<1000, tanTest(x) == tan(x)

#sin^2(x) + cos^2(x) = 1
function pythagFormula(x)
       sin(x)^2 + cos(x)^2
       end

       @test pythagFormula(1) == 1
       @test_formany -10<x<20 , pythagFormula(x) == 1

#test sin(-x) = -sin(x)
function trigidentities(x)
       sin(-x)
       end

       @test_formany 0<x<100, trigidentities(x) == -sin(x)

#test cos(-x) = cos(x)
function trigidentities2(x)
       cos(-x)
       end

       @test_formany 0<x<100, trigidentities2(x) == cos(x)

 #test tan(-x) = -tan(x)
function trigidentities3(x)
       tan(-x)
       end

       @test_formany 0<x<100, trigidentities3(x) == -tan(x)

#calculate volume function
function sphere_vol(r)
return 4/3*pi*r^3
end

#testing for all numbers between 1 and 100
	@test_formany 1<r<1000,sphere_vol(r) > 0
#will pass

#calculate quadratic formula/checkers cannot handle passing more than 1 variable
function quadratic2(a::Float64, b::Float64, c::Float64)

    sqr_term = sqrt(b^2-4a*c)
    r1 = quadratic(a, sqr_term, b)
    r2 = quadratic(a, -sqr_term, b)
r1, r2
end

#Test Failed
#Expression: (:((1 < a < 20,1 < b < 20,1 < c < 20,quad(a,b,c) > 0)),:(mode = test_formany))

#function return input
function stays(x)
       return x
       end

	@test_formany 1<x<10000, stays(x) == x

function additiveProperty(x)
       return x + 10
       end

	@test_formany 1<x<10000, additiveProperty(x) == 10 + x

function subtratProperty(x)
       return x - 10
       end

	@test_formany 1<x<10000, subtractProperty(x) == 10 - x

function multiplyProperty(x)
       return x * 10
       end

	@test_formany 1<x<10000, multiplyProperty(x) == 10 * x

	#associative property
	#cannot do this with checkers since it is only able to generate tests for one variable at a time
	#using Base Test package
	@test (1+2)+3 == 1 + (2 + 3)

function inverse(x)
	x + (-x)
	end

	@test_formany -1000<x<1000, inverse(x) == (-x) + x

function identity(x)
	x + 0
	end

	@test_formany -1000<x<1000, identity(x) == 0 + x

f(x) = 3x^6 - 7x + 1
	@test_exists 0<= x <=1, f(x) == 1

	#this gives an erro for some reason. What I am trying to test is that this test passes if at some point f(x) = 1 through 0<= x <=1.

#testing Periodicity of trig functions

 @test sin(1 + 2pi) == sin(1)
		Test Failed
  		Expression: sin(1 + 2pi) == sin(1)
  		 Evaluated: 0.8414709848078964 == 0.8414709848078965
ERROR: There was an error during testing

@test tan(1 + pi) == tan(1)
	  Test Failed
  	  Expression: tan(1 + pi) == tan(1)
   Evaluated: 1.5574077246549018 == 1.5574077246549023
ERROR: There was an error during testing

@test cos(1+2pi) == cos(1)
	  Test Failed
  	  Expression: cos(1 + 2pi) == cos(1)
   Evaluated: 0.5403023058681399 == 0.5403023058681398

#testing the double angle formula for sin
 @test sin(2*1) == 2*sin(1)*cos(1)
Test Failed
  Expression: sin(2 * 1) == 2 * sin(1) * cos(1)
   Evaluated: 0.9092974268256817 == 0.9092974268256818


text:
@testset "simplify-ps" begin
    s = "
aa
"
    @test F.simplify_ps(s) == "aa"
    s = "
aa"
    @test F.simplify_ps(s) == s
    s = "
aa
b
"
    @test F.simplify_ps(s) == s
end

# issue #549
@testset "div-note" begin
    # no nesting
    s = raw"""
    \newcommand{\note}[1]{@@note #1 @@}
    \note{hello}
    """ |> fd2html
    @test s // """

hello

        """
    # nesting 1 block
    s = raw"""
    \newcommand{\note}[1]{@@note @@title Note @@ #1 @@}
    \note{hello}
    """ |> fd2html
    @test s // """

Note

hello

        """
    #
    s = raw"""
    \newcommand{\note}[1]{@@note @@title Note @@ @@content #1 @@ @@}
    \note{hello}
    """ |> fd2html
    @test s // """

Note

hello

        """
end


text:
function ptp(x::AbstractArray{<:AbstractFloat, 1})
  return maximum(x) - minimum(x)
end

function ptp(x::AbstractArray{<:Complex{<:AbstractFloat}, 1})
  x_ = reinterpret(Float64, x)
  @views re′ = x_[1:2:end]
  @views im′ = x_[2:2:end]

  return ptp(re′) + 1.0im * ptp(im′)
end

function ptp(x::AbstractArray{<:AbstractFloat, N}, dims) where N
  return maximum(x, dims=dims) - minimum(x, dims=dims)
end

function ptp(x::AbstractArray{<:Complex{<:AbstractFloat}, N}, dims) where N
  x_ = reinterpret(Float64, x)
  @views re′ = reshape(x_[1:2:end], size(x)...)
  @views im′ = reshape(x_[2:2:end], size(x)...)

  return ptp(re′, dims) + 1.0im * ptp(im′, dims)
end


text:
module AntsArticle

export main, savedata, plotquality

using Serialization, DungBase

using Format2DB, Glob, DungAnalyse

using DataStructures, CoordinateTransformations, Rotations, DataFrames, Missings, Distributions, AngleBetweenVectors, LinearAlgebra, StatsBase, OnlineStats, Colors, PrettyTables, Measurements, HypothesisTests, GLM, DelimitedFiles, Printf

using CairoMakie, MakieLayout, FileIO, AbstractPlotting
import AbstractPlotting:px
CairoMakie.activate!()

include("preparedata.jl")
include("stats.jl")
include("plot.jl")
include("quality.jl")

"""
main()
Create all the tables and figures included in the manuscript.
"""
function main()
    data = deserialize("data")
    df = getdf(data)
    speeds!(df)
    descriptive_stats(df)
    displaced_stats(df)
    save_figures(df)
end

function goodpath(pathx)
    if isdir(pathx) && !isempty(readdir(pathx))
        x = basename(pathx)
        return !startswith(x, ['.', '_'])
    end
    return false
end

"""
savedata(path)
Convert the raw data to the standard dataset and save it.
The variable `path` is the path to the directory that
contains all the folders of all the experiments.
"""
function savedata(path)
    foreach(readdir(glob"source_*", tempdir())) do d
        rm(d, force = true, recursive = true)
    end
    todo = readdir(path, join = true)
    filter!(goodpath, todo)
    sources = Format2DB.main.(todo)
    source = DungAnalyse.joinsources(sources)
    foreach(sources) do d
        rm(d, force = true, recursive = true)
    end
    data = DungAnalyse.main(source)
    serialize("data", data)
end

"""
plotquality()
Plot some basic quality control.
"""
function plotquality()
    data = deserialize("data")
    df = getdf(data)
    plotquality(df)
end

end


text:
##################################################
## evaluate binary operations of symbolic objects
## XXX -- this may prove too narryw with the use of ::Sym

+(x::SymbolicObject, y::SymbolicObject) = x.__add__(y)
*(x::SymbolicObject, y::SymbolicObject) = x.__mul__(y)
-(x::SymbolicObject, y::SymbolicObject) = x.__sub__(y)
-(x::SymbolicObject)                    = x.__neg__()
/(x::SymbolicObject, y::SymbolicObject) = x.__div__(y)
^(x::SymbolicObject, y::SymbolicObject) = x.__pow__(y)
^(x::SymbolicObject, y::Rational) = x^convert(Sym,y)
#^(x::SymbolicObject, y::Integer) = x^convert(Sym,y) # no Union{Integer, Rational}, as that has ambiguity
//(x::SymbolicObject, y::Int) = x / Sym(y)
//(x::SymbolicObject, y::Rational) = x / Sym(y)
//(x::SymbolicObject, y::SymbolicObject) = x / y

\(x::SymbolicObject, y::SymbolicObject) = (y'/x')' # ?

inv(x::Sym) = x.__pow__(Sym(-1))


text:
module VPTrees

export
    # construct.jl
    Node,
    VPTree,

    # neighborlist.jl
    NeighborList,
    length,
    indices,
    distances,

    # knn.jl
    knn!,
    knn_sq!,

    # rangesearch.jl
    rangesearch!,
    rangesearch,
    rangesearch_sq!,
    rangesearch_sq,

    # display.jl
    printnode,
    show,
    display

include("construct.jl")
include("neighborlist.jl")
include("knn.jl")
include("rangesearch.jl")
include("display.jl")

end


text:
## connect do ProximalAlgorithms.jl as backend
ProximalCore.gradient!(grad, model::AbstractSem, parameters) = objective_gradient!(grad, model::AbstractSem, parameters)

#= function SemFit(optimization_result::Optim.MultivariateOptimizationResults, model::AbstractSem, start_val)
    return SemFit(
        optimization_result.minimum,
        optimization_result.minimizer,
        start_val,
        model,
        optimization_result
    )
end =#

function sem_fit(model::Sem{O, I, L, D}; start_val = start_val, kwargs...) where {O, I, L, D <: SemDiffProximal}

    if !isa(start_val, Vector)
        start_val = start_val(model; kwargs...)
    end

    if isnothing(model.diff.operator_h)
        solution, iterations = model.diff.algorithm(x0 = start_val, f = model, g = model.diff.operator_g)
    else
        solution, iterations = model.diff.algorithm(x0=start_val, f=model, g=model.diff.operator_g, h=model.diff.operator_h)
    end

    return SemFit(objective!(model, solution), solution, start_val, model, Dict(:iterations => iterations))

end


text:
export Interval

struct Interval
    a::Float64
    b::Float64
    Interval(a, b) = a > b ? error("invalid input") : new(a, b)
    Interval(a::Real) = new(a, a)
    Interval(a::Char) = a == '\U2205' ? new(Inf, -Inf) : error("invalid input") # U+2205 is ∅
end

Base.minimum(A::Interval) = A.a
Base.maximum(A::Interval) = A.b

function Base.:(∈)(a::Real, A::Interval)
    if a >= A.a && a <= A.b
        return true
    else
        return false
    end
end

function Base.isempty(A::Interval)
    if A.a > A.b
        return true
    else
        return false
    end
end

function Base.issubset(A::Interval, B::Interval)
    if A.a >= B.a && A.b <= B.b
return true
    else
        return false
    end
end

function Base.intersect(A::Interval, B::Interval)
    if maximum(A) >= minimum(B) && maximum(B) >= minimum(A)
        return Interval(max(minimum(A), minimum(B)), min(maximum(A), maximum(B)))
    else
        return Interval('∅')
    end
end

function Base.show(io::IO, A::Interval)
    if isempty(A)
        print(io, '∅')
    else
        print(io, '\u301a')
        print(io, A.a)
        print(io, ", ")
        print(io, A.b)
        print(io, '\u301b')
    end
end