text:
using Test, JuMP

using Distributed
addprocs(2)
@everywhere using CharibdeOptim
@everywhere using IntervalArithmetic

@testset "Using Charibde for Constrained Optimsation" begin
      @everywhere using ModelingToolkit

      @everywhere vars = ModelingToolkit.@variables x y
      @everywhere C1 = constraint(vars, x+y, -Inf..4)
      @everywhere C2 = constraint(vars, x+3y, -Inf..9)
      @everywhere constraints = [C1, C2]
      (maxima, maximisers, info) = charibde_max(X->((x,y)=X;-(x-4)^2-(y-4)^2), IntervalBox(-4..4, -4..4), constraints)
      @test maxima ⊆ -8.01 .. -7.99
      @test maximisers[1] ⊆ (1.99 .. 2.01) × (1.99 .. 2.01)
end

@testset "Using JuMP syntax for Constrained Optimisation" begin
      model = Model(with_optimizer(CharibdeOptim.Optimizer))
      @variable(model, -4 <= x <= 4)
      @variable(model, -4 <= y <= 4)
      @NLconstraint(model, x+y<=4)
      @NLconstraint(model, 5<=x+3y<=9)
      @NLobjective(model, Max, -(x-4)^2-(y-4)^2)
      optimize!(model)

      @test JuMP.termination_status(model) == MOI.OPTIMAL
      @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
      @test JuMP.objective_value(model) ⊆ -8.01 .. -7.99
      @test JuMP.value(x) ⊆ (1.99 .. 2.01)
      @test JuMP.value(y) ⊆ (1.99 .. 2.01)
end

@testset "Using Interval bound and contract algorithm for Constrained Optimisation" begin
      vars = ModelingToolkit.@variables x y
      C1 = constraint(vars, x+y, -Inf..4)
      C2 = constraint(vars, x+3y, -Inf..9)

      (maxima, maximisers, info) = ibc_maximise(X->((x,y)=X;-(x-4)^2-(y-4)^2), IntervalBox(-4..4, -4..4),[C1, C2])
      @test maxima ⊆ -8.01 .. -7.99
      @test maximisers[1] ⊆ (1.99 .. 2.01) × (1.99 .. 2.01)
end

@testset "Optimising by Interval Branch and Contract Algorithm" begin
      (global_min, minimisers)= ibc_minimise(X->((x,y)= X;x^2 + y^2), IntervalBox(2..3, 3..4))
      @test global_min ⊆ 13 .. 13.01
      @test minimisers[1] ⊆ (2.0 .. 2.001) × (3.0 .. 3.001)
end

@testset "Optimising by Differencial Evolution" begin
      (global_min, minimiser)= diffevol_minimise(X->((x,y)= X;x^2 + y^2), IntervalBox(2..3, 3..4), iterations = 40)
      @test global_min ⊆ 13 .. 13.01
      @test minimiser ⊆ (2.0 .. 2.001) × (3.0 .. 3.001)
end

@testset "Optimising by Charibde (A hybrid approach) using only one worker" begin

      (global_min, minimisers, info) = charibde_min(X->((x,y)=X;x^2+y+1), IntervalBox(1..2, 2..3), workers = 1)
      @test global_min ⊆ 4.0 .. 4.01
      @test minimisers[1] ⊆ (1..1.001) × (2..2.001)

      (global_min, minimisers, info)= charibde_min(X->((x,y)=X;x^2 + y^2), IntervalBox(2..3, 3..4), workers = 1)
      @test global_min ⊆ 13 .. 13.01
      @test minimisers[1] ⊆ (2.0 .. 2.001) × (3..3.001)
end

@testset "Using JuMP syntax by using only one worker " begin           #for using two workers just dont pass 'workers' arguments as its value is set to 2
      model = Model(with_optimizer(CharibdeOptim.Optimizer, workers = 1))
      @variable(model, 1<=x<=2)
      @variable(model, 1<=y<=2)
      @NLobjective(model, Min, x^2+y^2)
      optimize!(model)

      @test JuMP.termination_status(model) == MOI.OPTIMAL
      @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
      @test JuMP.objective_value(model) ≈ 2.0
      @test JuMP.value(x) ≈ 1.0
      @test JuMP.value(y) ≈ 1.0
end

 # No need to add worker because a worker is already added while running testset "Using Charibde for Constrained Optimsation".
 # Otherwise we have to add a worker by using 'Distributed.addprocs(1)' and load the package on each worker
 # by '@everywhere using CharibdeOptim'.

@testset "Optimising by Charibde (A hybrid approach) using 2 workers" begin
      (global_min, minimisers, info) = charibde_min(X->((x,y)=X;x^3 + 2y + 5), IntervalBox(2..4, 2..3))
      @test global_min ⊆ 17.0 .. 17.01
      @test minimisers[1] ⊆ (2..2.001) × (2..2.001)

      (global_min, minimisers, info)= charibde_min(X->((x,y)=X;x^2 + y^2), IntervalBox(2..3, 3..4))
      @test global_min ⊆ 13 .. 13.01
      @test minimisers[1] ⊆ (2.0 .. 2.001) × (3..3.001)

end

@testset "Optimising difficult problem using JuMP" begin
      model = Model(with_optimizer(CharibdeOptim.Optimizer))

      @variable(model, 2 <= x1 <= 3)
      @variable(model, 3 <= x2 <= 4)
      @variable(model, 9 <= x3 <= 14)
      @variable(model, 2 <= x4 <= 3)
      @variable(model, 3 <= x5 <= 4)
      @variable(model, 9 <= x6 <= 14)
      @variable(model, 2 <= x7 <= 3)
      @variable(model, 3 <= x8 <= 4)
      @variable(model, 9 <= x9 <= 14)

      @NLobjective(model, Min, x1^2 + x2^2 + x3^4 - x4^7 - 200x5 - x6^5 - x7^9 + x8^5 - 8x9^3)

      optimize!(model)

      @test JuMP.termination_status(model) == MOI.OPTIMAL
      @test JuMP.primal_status(model) == MOI.FEASIBLE_POINT
      @test JuMP.objective_value(model) ≈ -575629.0
      @test JuMP.value(x1) ⊆ 1.99 .. 2.01
      @test JuMP.value(x2) ⊆ 2.99 .. 3.01
      @test JuMP.value(x3) ⊆ 8.99 .. 9.01
      @test JuMP.value(x4) ⊆ 2.99 .. 3.01
      @test JuMP.value(x5) ⊆ 3.99 .. 4.01
      @test JuMP.value(x6) ⊆ 13.99 .. 14.01
      @test JuMP.value(x7) ⊆ 2.99 .. 3.01
      @test JuMP.value(x8) ⊆ 2.99 .. 3.01
      @test JuMP.value(x9) ⊆ 13.99 .. 14.01
end

@testset "Optimising difficult problems using Charibde" begin
      f = X->((x1,x2,x3,x4,x5,x6,x7,x8,x9)=X;x1^2 + x2^2 + x3^4 - x4^7 - 200x5 - x6^5 - x7^9 + x8^5 - 8x9^3)
      X = IntervalBox(2..3, 3..4, 9..14, 2..3, 3..4, 9..14, 2..3, 3..4, 9..14)
      (global_min, minimisers, info) = charibde_min(f, X)

      @test global_min ⊆ -575630 .. -575628
      @test minimisers[1] ⊆ (1.99 .. 2.01) × (2.99 .. 3.01) × (8.99 .. 9.01) × (2.99 .. 3.01) × (3.99 .. 4.01) × (13.99 .. 14.01) × (2.99 .. 3.01) × (2.99 .. 3.01) × (13.99 .. 14.01)

end

max_issues_repo_name: yashcodes/Charibde

max_forks_repo_name: yashcodes/Charibde


text:
function feed(obj::DataProvider, arg0::Pattern, arg1::DataLoader, arg2::DataProvidersManager)
    return jcall(obj, "feed", jboolean, (Pattern, DataLoader, DataProvidersManager), arg0, arg1, arg2)
end

ext: jl


text:
"""
# Constants List:
- Physical Constants
   - `mu_0`: Permeability of vacuum, Unit: N/A^2;
   - `epsilon_0`: Vacuum dielectric constant, Unit: F/m;
   - `e`: elementary charge, Unit: C;
   - `k`: Electrostatic constant, Unit: N*m^2/C^2;
   - `k_b`: Boltzmann constant, Unit: m^2 kg / (s^2 K);
   - `n_a`: Avogadro constant, Unit: 1/mol;
   - `density_wat`: density of water, Unit: kg/m^3;
- Energy Converters:
   - `kcal2j`;
   - `kcalm2j`;
   - `kcalm2t`: kcal/mol to temparture in unit of K;
- Mass Converters:
   - `g2kg`;
   - `kg2g`;
   - `gm2g`;
   - `gm2kg`;
- Length Converters
   - `an2m`;
   - `an2nm`;
   - `nm2m`;
   - `cm2an`;
   - `dm2m`;
   - `m2dm`;
- Time Converters
   - `fs2s`;
   - `ps2s`;
   - `ns2s`;
"""
function constants()
    res = Dict()
    # Physical Constants
    res["mu_0"] = 4e-7 * π          # Permeability of vacuum
    res["epsilon_0"] = 8.854188e-12 # Vacuum dielectric constant
    res["e"] = 1.6e-19;             # Elementary charge
    res["k"] = 9e9;                 # Unit N*m^2/C^2
    res["k_b"] = 1.38065e-23;
    res["n_a"] = 6.02214e23;
    res["density_wat"] = 1e3;       # Unit: kg/m^3
    # Energy Converters
    res["kcal2j"] = 4.184e3;
    res["kcalm2j"] = Const_kcal2j/Const_n_a;
    res["kcalm2t"] = Const_kcalm2j/Const_k_b;
    # Mass Converters
    res["g2kg"] = 1e-3;
    res["kg2g"] = 1e3;
    res["gm2g"] = 1/Const_n_a;
    res["gm2kg"] = Const_g2kg/Const_n_a;
    # Length Converters
    res["an2m"] = 1e-10;
    res["an2nm"] = 1e-1;
    res["nm2m"] = 1e-9;
    res["cm2an"] = 1e8;
    res["dm2m"] = 1e-1;
    res["m2dm"] = 1/Const_dm2m;
    # Time Converters
    res["fs2s"] = 1e-15;
    res["ps2s"] = 1e-12;
    res["ns2s"] = 1e-9;

    return res
end

# Physical Constants
Const_mu_0 = 4e-7 * π;
Const_epsilong_0 = 8.854188e-12;
Const_e = 1.6e-19;
Const_k = 9e9;
Const_k_b = 1.38065e-23;
Const_n_a = 6.02214e23;
Const_density_wat = 1e3;      # Unit: kg/m^3
# Energy Converters
Const_kcal2j = 4.184e3;
Const_kcalm2j = Const_kcal2j/Const_n_a;
Const_kcalm2t = Const_kcalm2j/Const_k_b;
# Mass Converters
Const_g2kg = 1e-3;
Const_kg2g = 1e3;
Const_gm2g = 1/Const_n_a;
Const_gm2kg = Const_g2kg/Const_n_a;
# Length Converters
Const_an2m = 1e-10;
Const_an2nm = 1e-1;
Const_nm2m = 1e-9;
Const_cm2an = 1e8;
Const_dm2m = 1e-1;
Const_m2dm = 1/Const_dm2m;
# Time Converters
Const_fs2s = 1e-15;
Const_ps2s = 1e-12;
Const_ns2s = 1e-9;


text:
using jlmie
using Plots

# # xz plot range
# x = (-150:5:150)
# y = 0.0
# z = (-150:5:150)
# # calculation using broadcast
# r, theta, phi = cart2sph(x',y,z)
# # heatmap(x,z,r)
# # heatmap(x,z,theta)
# # heatmap(x,z,phi)

# meshgrid type calculation
# mat_x = x' .* ones(length(z))
# mat_z = ones(length(x))' .* z
# mat_r,mat_theta,mat_phi = cart2sph(mat_x,y,mat_z)
# calculation using broadcast
# r, theta, phi = cart2sph(x',y,z)
# heatmap(x,z,r)
# heatmap(x,z,theta)
# heatmap(x,z,phi)

x = (-150:1:150)
y = (-150:1:150)
z = (-150:1:150)

xyz = 2  # 1:x, 2:y, 3:z
if xyz == 1
    # xz plot range
    mat_x = x' .* ones(length(z))
    mat_z = ones(length(x))' .* z
    mat_r,mat_theta,mat_phi = cart2sph(mat_x,0,mat_z)
    p1 = heatmap(x,z,mat_r)
    p2 = heatmap(x,z,mat_theta)
    p3 = heatmap(x,z,mat_phi)
elseif xyz == 2
    # yz plot range
    mat_y = y' .* ones(length(y))
    mat_z = ones(length(y))' .* z
    mat_r,mat_theta,mat_phi = cart2sph(0,mat_y,mat_z)
    p1 = heatmap(y,z,mat_r)
    p2 = heatmap(y,z,mat_theta)
    p3 = heatmap(y,z,mat_phi)
elseif xyz == 3
    # xy plot range
    mat_x = x' .* ones(length(y))
    mat_y = ones(length(x))' .* y
    mat_r,mat_theta,mat_phi = cart2sph(mat_x,mat_y,0)
    p1 = heatmap(x,y,mat_r)
    p2 = heatmap(x,y,mat_theta)
    p3 = heatmap(x,y,mat_phi)
end
plot(p1,p2,p3)

max_stars_repo_name: Hinamoooon/jlmie

max_issues_repo_name: Hinamoooon/jlmie

max_forks_repo_name: Hinamoooon/jlmie


text:
module test_julia_sort

using Test

@test sortperm([5,6,7]) == 1:3
@test sortperm([6,5,7]) == [2,1,3]
@test sort([6,5,7]) == [5,6,7]

end # module test_julia_sort


text:
# Pentagonal numbers are generated by the formula, Pn=n(3n−1)/2. The first ten
# pentagonal numbers are:
#
# 1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ...
#
# It can be seen that P4 + P7 = 22 + 70 = 92 = P8. However, their difference,
# 70 − 22 = 48, is not pentagonal.
#
# Find the pair of pentagonal numbers, Pj and Pk, for which their sum and
# difference are pentagonal and D = |Pk − Pj| is minimised; what is the value
# of D?

using ProjectEulerSolutions

# Calculate pentagonal numbers, store in set, then check for whether
# differences and additions are also in the set.  Could probably be faster if
# searchsorted was used and the minimum bound for j was a function of i.
function p044solution(n::Integer=10)::Integer

    vals = collect(1:n)
    pentagonals = fld.(vals .* (3 .* vals .- 1), 2)
    pentagonalset = Set(pentagonals)

    for i in 4:n
        for j in 3:(i-1)
            if pentagonals[i] - pentagonals[j] in pentagonalset && pentagonals[i] + pentagonals[j] in pentagonalset
                return pentagonals[i] -  pentagonals[j]
            end
        end
    end

    return -1
end

p044 = Problems.Problem(p044solution)

Problems.benchmark(p044, 3000)

max_stars_repo_name: gnujosh/julia-euler-project

max_issues_repo_name: gnujosh/julia-euler-project

max_forks_repo_name: gnujosh/julia-euler-project


text:
# write your own tests here

include("testbits.jl")
include("testresults.jl")
include("testdistances.jl")
include("knr/testkvp.jl")
include("knr/testknr.jl")
include("graph/testvectorindex.jl")


text:
"""
Abstract type for time_series that are stored in a system.
Users never create them or get access to them.
Stores references to TimeSeriesData.
"""
abstract type TimeSeriesMetadata <: InfrastructureSystemsType end

abstract type ForecastMetadata <: TimeSeriesMetadata end

abstract type StaticTimeSeriesMetadata <: TimeSeriesMetadata end

get_count(ts::StaticTimeSeriesMetadata) = 1
get_initial_timestamp(ts::StaticTimeSeriesMetadata) = get_initial_timestamp(ts)
Base.length(ts::StaticTimeSeriesMetadata) = get_length(ts)
Base.length(ts::ForecastMetadata) = get_horizon(ts)

"""
Abstract type for time series stored in the system.
Components store references to these through TimeSeriesMetadata values so that data can
reside on storage media instead of memory.
"""
abstract type TimeSeriesData <: InfrastructureSystemsComponent end

# Subtypes must implement
# - Base.length
# - get_resolution
# - make_time_array
# - eltype_data


text:
# Tools
function extract_pkgdeps(proj_file)
    raw_deps = get(TOML.parsefile(proj_file), "deps", Dict())
    deps = Dict()
    for (name, uuid) in raw_deps
        deps[uuid] = name
    end
    deps
end

function find_projects(root = ".")
    projs = String[]
    c = 1
    root = abspath(root)

    # Look in root
    for projname in Base.project_names
        projfile = joinpath(root, projname)
        !isfile(projfile) && continue
        push!(projs, projfile)
        println(c, ": ", relpath(projfile, root))
        c += 1
    end

    # walk root
    for (dir, subdirs, _) in walkdir(root)
        for subdir in subdirs, projname in Base.project_names
            try
                projfile = joinpath(dir, subdir, projname)
                !isfile(projfile) && continue

                # ignores
                occursin("/.git/", projfile) && continue

                push!(projs, projfile)
                println(c, ": ", relpath(projfile, root))
                c += 1
            catch err; end
        end
    end
    projs
end

# Find package in registers
function findin_regs(name, uuidpkg)
    firstletter = string(first(name))
    founds = []
    for depot in DEPOT_PATH
        regs_root = joinpath(depot, "registries")
        !isdir(regs_root) && continue
        for reg in readdir(regs_root)
            reg_pkgdir = joinpath(regs_root, reg, firstletter, name)
            !isdir(reg_pkgdir) && continue

            # Check uuid
            pkg_toml = joinpath(reg_pkgdir, "Package.toml")
            !isfile(pkg_toml) && continue
            uuidfile = get(TOML.parsefile(pkg_toml), "uuid", "NOT_A_UUID")
            uuidfile != uuidpkg && continue

            push!(founds, reg_pkgdir)
        end
    end
    founds
end

# Extract existing versions from register
function extract_reg_pkgversions(regpath)
    ver_file = joinpath(regpath, "Versions.toml")
    versions = TOML.parsefile(ver_file)
    versions |> keys .|> VersionNumber
end

function extract_reg_pkguuid(regpath)
    pkg_file = joinpath(regpath, "Package.toml")
    dat = TOML.parsefile(pkg_file)
    return dat["uuid"]
end

function tempenv()
    temp_proj_dir = joinpath(tempdir(), string(rand(UInt)))
    mkpath(temp_proj_dir)
    return temp_proj_dir
end

function tempenv(proj_dir)
    temp_proj_dir = tempenv()
    for file_names in [Base.project_names, Base.manifest_names]
        for file_name in file_names
            file = joinpath(file_name, proj_dir)
            !isfile(file) && continue
            cp(file, temp_proj_dir; force = true, follow_symlinks = true)
        end
    end
    return temp_proj_dir
end


text:
using NetworkModels
using Documenter

DocMeta.setdocmeta!(NetworkModels, :DocTestSetup, :(using NetworkModels); recursive=true)

makedocs(;
    modules=[NetworkModels],
    authors="Cédric Simal, University of Namur",
    repo="https://github.com/csimal/NetworkModels.jl/blob/{commit}{path}#{line}",
    sitename="NetworkModels.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://csimal.github.io/NetworkModels.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/csimal/NetworkModels.jl",
    devbranch="main",
)


text:
using Test
using MuttsDicts

@testset "config_for_size" begin
    prev_config = MuttsDicts._config_for_size(1)
    while prev_config.next_config_size < 2^48
        c = MuttsDicts._config_for_size(prev_config.next_config_size-1)
        @test "$(c)" == "$(prev_config)"
        new_config = MuttsDicts._config_for_size(prev_config.next_config_size)
        @assert new_config.next_config_size > prev_config.next_config_size
        prev_config = new_config
    end
end

@testset "small_dict" begin
    dict = MuttsDict{Int,Int}()
    for i=1:100
        key = i
        value = i*317
        dict[key] = value
        @test length(dict) == i

        check = Tuple{Int,Int}[]
        for (ki,vi) in dict
            push!(check, (ki,vi))
        end
        sort!(check)
        @test length(check) == i
        for k=1:length(check)
            @test check[k] == (k,k*317)
        end
    end
end

@testset "dict" begin
    dict = MuttsDict{Int,Int}()
    N=2^20
    for i=1:N
        key = i
        value = i*317
        insert!(dict, key, value)
        @test getindex(dict,key) == value
        if (i == 7) || (i == 40) || (i == 120) || (i == 1000) || (mod(i,10000) == 0)
            prev_dict = dict
            dict = branch!(dict)
            @test dict !== prev_dict
            @test branch!(prev_dict) !== dict

            check = Vector{Tuple{Int,Int}}()
            for (k,v) in dict
                push!(check, (k,v))
            end
            @test length(check) == i
            sort!(check)
            for k=1:length(check)
                @test check[k] == (k,k*317)
            end
        end
    end
    for i=1:N
        @test getindex(dict, i) == i*317
    end
end

@testset "deletes" begin
    dict = MuttsDict{Int,Int}()
    N=2^20
    for i=1:N
        insert!(dict,i,i*317)
    end

    dict = branch!(dict)
    for i=1:N
        @test in((i,i*317),dict)
        @test in((i=>i*317),dict)
        delete!(dict,i)
        @test !in((i,i*317),dict)
    end

    empty_dict = MuttsDict{Int,Int}()
    diff = setdiff(dict,empty_dict)
    @test length(diff) == 0
end

@testset "setdiff" begin
    prev_dict = MuttsDict{Int,Int}()
    dict = branch!(prev_dict)
    N=2^20
    check_m = 100
    for i=1:N
        insert!(dict, i, i*317)
        if mod(i,check_m)==0
            diff = setdiff(dict,prev_dict)
            sort!(diff)
            @test length(diff) == check_m
            for k=1:check_m
                pk = i-check_m+k
                @test diff[k] == (pk, pk*317)
                if diff[k] != (pk,pk*317)
                    abort()
                end
            end
            diff2 = setdiff(prev_dict,dict)
            @test length(diff2) == 0
            prev_dict = dict
            dict = branch!(prev_dict)
        end
    end
end


text:
using Rasters, ArchGDAL, Dates, Pkg

function load_data(basedir, carrycap)
    # The directory for this project
    basedir = dirname(Pkg.project().path)
    # Load growthrates from tif files saved in growthmaps.jl
    datapath = joinpath(basedir, "data", "growthrates")
    # Set the timespan of the data. This will be looped by DynamicGrids
    # for years before or after this period.
    growthtimespan = DateTime(2020, 1):Month(1):DateTime(2020, 12)
    # Load, cut out Australia, and reorient the growth rate data
    aus = X(Between(113.0, 154.0)), Y(Between(-44.0, -10.0))
    rate_h, rate_p = map((:host, :parasitoid)) do species
        slices = map(readdir(joinpath(datapath, string(species)); join=true)) do path
            Raster(path; mappedcrs=EPSG(4326))[aus..., Band(1)] |>
                a -> replace_missing(a, 0.0f0) |>
                a -> permutedims(a, (Y, X))
        end
        A = cat(slices...; dims=Ti(growthtimespan))
    end

    # Take an array slice to use as a template
    A = rate_h[Ti(1)]

    # Make a masking layer for ocean areas
    mask = boolmask(A)

    # Create an initial population array for the host
    init_h = zero(A)
    # Add a host population around Cairns in Queensland
    cairns = Y(Between(-17, -20)), X(Between(144.0, 145.0))
    init_h[cairns...] .= carrycap

    # Initialise scattered parasitoid populations
    init_p = (_ -> (r = rand(); r < 0.1 ? carrycap * 1f-4 : 0.0f0)).(A)

    return (; init_h, init_p, rate_h, rate_p, mask)
end


text:
"""
see https://arxiv.org/abs/1701.07035
"""
@with_kw struct Vumps <: Algorithm
    tol_galerkin::Float64 = Defaults.tol
    tol_gauge::Float64 = Defaults.tolgauge
    maxiter::Int = Defaults.maxiter
    orthmaxiter::Int = Defaults.maxiter
    finalize::Function = (iter,state,ham,pars) -> (state,pars, true);
    verbose::Bool = Defaults.verbose
end

"
    find_groundstate(state,ham,alg,pars=params(state,ham))

    find the groundstate for ham using algorithm alg
"
function find_groundstate(state::InfiniteMPS, H::Hamiltonian,alg::Vumps,pars=params(state,H))
    galerkin  = 1+alg.tol_galerkin
    iter      = 1

    external_conv = false
    while true
        eigalg = Arnoldi(tol=galerkin/(4*sqrt(iter)))

        acjobs = map(enumerate(state.AC)) do (loc,ac)
            # @maarten, als ik julia start met 1 thread heb ik geen memory leak. Kan het zijn dat dit de oorzaak is ?
            # threads blijven mss openstaan ofzo met meer en meer ram
	    # delete dit maar als dit zever is ;)
            @Threads.spawn eigsolve(ac, 1, :SR, eigalg) do x
                ac_prime(x, loc, state, pars)
            end
        end

        cjobs = map(enumerate(state.CR)) do (loc,cr)
            @Threads.spawn eigsolve(cr, 1, :SR, eigalg) do x
                c_prime(x, loc, state, pars)
            end
        end

        newAs = map(zip(acjobs,cjobs)) do (acj,cj)
            (e,vac,ch) = fetch(acj)
            (e,vc,ch) = fetch(cj)

            QAc,_ = TensorKit.leftorth!(vac[1], alg=QRpos())
            Qc,_  = TensorKit.leftorth!(vc[1], alg=QRpos())

            QAc*adjoint(Qc)
        end

        state = InfiniteMPS(newAs; tol = alg.tol_gauge, maxiter = alg.orthmaxiter,leftgauged=true)
        galerkin   = calc_galerkin(state, pars)
        alg.verbose && @info "vumps @iteration $(iter) galerkin = $(galerkin)"

        (state,pars, external_conv) = alg.finalize(iter,state,H,pars);
        if (galerkin <= alg.tol_galerkin && external_conv ) || iter>=alg.maxiter
            iter>=alg.maxiter && println("vumps didn't converge $(galerkin)")
            return state, pars, galerkin
        end

        iter += 1
    end
end

"calculates the galerkin error"
calc_galerkin(state::InfiniteMPS, pars) = maximum([norm(leftnull(state.AC[loc])'*ac_prime(state.AC[loc], loc, state, pars)) for loc in 1:length(state)])


text:
#
# Date created: 2022-03-09
# Author: aradclif
#
#
############################################################################################

function sumsqdiff_body(N::Int, D)
    body = Expr(:block)
    a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
    params = D.parameters
    b = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
    δ = Expr(:(=), :Δ, Expr(:call, :-, a, b))
    e = Expr(:(=), c, Expr(:call, :+, c, Expr(:call, :*, :Δ, :Δ)))
    push!(body.args, δ)
    push!(body.args, e)
    body
end

function sumsqdiff_quote(N::Int, D)
    ls = loopgen(N)
    body = sumsqdiff_body(N, D)
    push!(ls.args, body)
    return quote
        @turbo $ls
        return C
    end
end
@generated function sumsqdiff!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                               B::AbstractArray{Tₘ, N}, dims::D) where {Tₒ, T, Tₘ, N, D}
    sumsqdiff_quote(N, D)
end

function lvvar(A::AbstractArray{T, N}, dims::NTuple{M, Int}, corrected::Bool) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = hvncat(ntuple(i -> 1, Val(N)), true, lvvar1(A))
    else
        B = lvmean(A, dims=dims)
        Dᴮ = size(B)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        Tₒ = Base.promote_op(/, T, Int)
        C = zeros(Tₒ, Dᴮ)
        sumsqdiff!(C, A, B, Dᴮ′)
        Dᴬ = size(A)
        denom = one(Tₒ)
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        denom = corrected ? denom - one(Tₒ) : denom
        x = inv(denom)
        smul!(C, x)
    end
    return C
end
lvvar(A::AbstractArray{T, N}, dims::Int, corrected) where {T, N} = lvvar(A, (dims,), corrected)
lvvar(A::AbstractArray{T, N}; dims=:, corrected=true) where {T, N} = lvvar(A, dims, corrected)
lvvar(A::AbstractArray{T, N}, ::Colon, corrected) where {T, N} = lvvar1(A, corrected)

function lvvar1(A::AbstractArray{T, N}, corrected::Bool=true) where {T, N}
    s = zero(Base.promote_op(+, T, Int))
    @turbo for i ∈ eachindex(A)
        s += A[i]
    end
    μ = s / length(A)
    ss = zero(Base.promote_op(/, T, Int))
    @turbo for i ∈ eachindex(A)
        Δ = A[i] - μ
        ss += Δ * Δ
    end
    return corrected ? ss / (length(A) - 1) : ss / length(A)
end

################ threaded version

function tsumsqdiff_quote(N::Int, D)
    ls = loopgen(N)
    body = sumsqdiff_body(N, D)
    push!(ls.args, body)
    return quote
        @tturbo $ls
        return C
    end
end
@generated function tsumsqdiff!(C::AbstractArray{Tₒ, N}, A::AbstractArray{T, N},
                                B::AbstractArray{Tₘ, N}, dims::D) where {Tₒ, T, Tₘ, N, D}
    tsumsqdiff_quote(N, D)
end

function lvtvar(A::AbstractArray{T, N}, dims::NTuple{M, Int}, corrected::Bool) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = hvncat(ntuple(i -> 1, Val(N)), true, lvtvar1(A))
    else
        B = lvmean(A, dims=dims)
        Dᴮ = size(B)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        Tₒ = Base.promote_op(/, T, Int)
        C = zeros(Tₒ, Dᴮ)
        tsumsqdiff!(C, A, B, Dᴮ′)
        Dᴬ = size(A)
        denom = one(Tₒ)
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        denom = corrected ? denom - one(Tₒ) : denom
        x = inv(denom)
        tsmul!(C, x)
    end
    return C
end
lvtvar(A::AbstractArray{T, N}, dims::Int, corrected) where {T, N} = lvtvar(A, (dims,), corrected)
lvtvar(A::AbstractArray{T, N}; dims=:, corrected=true) where {T, N} = lvtvar(A, dims, corrected)
lvtvar(A::AbstractArray{T, N}, ::Colon, corrected) where {T, N} = lvtvar1(A, corrected)

function lvtvar1(A::AbstractArray{T, N}, corrected::Bool=true) where {T, N}
    s = zero(Base.promote_op(+, T, Int))
    @tturbo for i ∈ eachindex(A)
        s += A[i]
    end
    μ = s / length(A)
    ss = zero(Base.promote_op(/, T, Int))
    @tturbo for i ∈ eachindex(A)
        Δ = A[i] - μ
        ss += Δ * Δ
    end
    return corrected ? ss / (length(A) - 1) : ss / length(A)
end

################
# function smul_sqrt_inplacebody(N::Int)
#     body = Expr(:block)
#     a = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     e = Expr(:(=), a, Expr(:call, :sqrt, Expr(:call, :*, a, :x)))
#     push!(body.args, e)
#     body
# end

# function smul_sqrt_inplace_quote(N::Int)
#     ls = loopgen(N)
#     body = smul_sqrt_inplacebody(N)
#     push!(ls.args, body)
#     return quote
#         @turbo $ls
#         return A
#     end
# end

# @generated function smul_sqrt!(A::AbstractArray{T, N}, x::T) where {T, N}
#     smul_sqrt_inplace_quote(N)
# end

function smul_sqrt!(A::AbstractArray{T, N}, x::T) where {T, N}
    @turbo for i ∈ eachindex(A)
        A[i] = sqrt(A[i] * x)
    end
    A
end

function lvstd(A::AbstractArray{T, N}, dims::NTuple{M, Int}, corrected::Bool) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = hvncat(ntuple(i -> 1, Val(N)), true, lvstd1(A))
    else
        B = lvmean(A, dims=dims)
        Dᴮ = size(B)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        Tₒ = Base.promote_op(/, T, Int)
        C = zeros(Tₒ, Dᴮ)
        sumsqdiff!(C, A, B, Dᴮ′)
        Dᴬ = size(A)
        denom = one(Tₒ)
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        denom = corrected ? denom - one(Tₒ) : denom
        x = inv(denom)
        smul_sqrt!(C, x)
    end
    return C
end
lvstd(A::AbstractArray{T, N}, dims::Int, corrected) where {T, N} = lvstd(A, (dims,), corrected)
lvstd(A::AbstractArray{T, N}; dims=:, corrected=true) where {T, N} = lvstd(A, dims, corrected)
lvstd(A::AbstractArray{T, N}, ::Colon, corrected) where {T, N} = lvstd1(A, corrected)

lvstd1(A::AbstractArray{T, N}, corrected::Bool=true) where {T, N} = √(lvvar1(A, corrected))

################ threaded version
function tsmul_sqrt!(A::AbstractArray{T, N}, x::T) where {T, N}
    @tturbo for i ∈ eachindex(A)
        A[i] = sqrt(A[i] * x)
    end
    A
end

function lvtstd(A::AbstractArray{T, N}, dims::NTuple{M, Int}, corrected::Bool) where {T, N, M}
    if ntuple(identity, Val(N)) ⊆ dims
        C = hvncat(ntuple(i -> 1, Val(N)), true, lvtstd1(A))
    else
        B = lvmean(A, dims=dims)
        Dᴮ = size(B)
        Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
        Tₒ = Base.promote_op(/, T, Int)
        C = zeros(Tₒ, Dᴮ)
        tsumsqdiff!(C, A, B, Dᴮ′)
        Dᴬ = size(A)
        denom = one(Tₒ)
        for d ∈ eachindex(Dᴬ)
            denom = d ∈ dims ? denom * Dᴬ[d] : denom
        end
        denom = corrected ? denom - one(Tₒ) : denom
        x = inv(denom)
        tsmul_sqrt!(C, x)
    end
    return C
end
lvtstd(A::AbstractArray{T, N}, dims::Int, corrected) where {T, N} = lvtstd(A, (dims,), corrected)
lvtstd(A::AbstractArray{T, N}; dims=:, corrected=true) where {T, N} = lvtstd(A, dims, corrected)
lvtstd(A::AbstractArray{T, N}, ::Colon, corrected) where {T, N} = lvtstd1(A, corrected)

lvtstd1(A::AbstractArray{T, N}, corrected::Bool=true) where {T, N} = √(lvtvar1(A, corrected))

# ################################################################
# # covariance in multiple dimensions

# function outerloopgen(N::Int)
#     loops = Expr(:for)
#     block = Expr(:block)
#     for d = N:-1:1
#         ex = Expr(:(=), Symbol(:j_, d), Expr(:call, :axes, :B, d))
#         push!(block.args, ex)
#     end
#     push!(loops.args, block)
#     loops
# end
# function outerloopgen(N::Int, D)
#     loops = Expr(:for)
#     block = Expr(:block)
#     params = D.parameters
#     for d = N:-1:1
#         if params[d] != Static.One
#             ex = Expr(:(=), Symbol(:j_, d), Expr(:call, :axes, :B, d))
#             push!(block.args, ex)
#         end
#     end
#     push!(loops.args, block)
#     loops
# end

# function sumsqdiffij_body(N::Int, D)
#     body = Expr(:block)
#     params = D.parameters
#     # # As it was originally conceived, but the more efficient form (iᵢ, iⱼ)
#     # # actually confers enhanced clarity.
#     # aᵢ = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     # bᵢ = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...)
#     # aⱼ = Expr(:ref, :A, ntuple(d -> Symbol(:j_, d), N)...)
#     # bⱼ = Expr(:ref, :B, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:j_, d), N)...)
#     # c = Expr(:ref, :C, ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)...,
#     #          ntuple(d -> params[d] == Static.One ? 1 : Symbol(:j_, d), N)...)
#     iᵢ = ntuple(d -> params[d] == Static.One ? 1 : Symbol(:i_, d), N)
#     iⱼ = ntuple(d -> params[d] == Static.One ? 1 : Symbol(:j_, d), N)
#     aᵢ = Expr(:ref, :A, ntuple(d -> Symbol(:i_, d), N)...)
#     bᵢ = Expr(:ref, :B, iᵢ...)
#     aⱼ = Expr(:ref, :A, ntuple(d -> params[d] == Static.One ? Symbol(:i_, d) : Symbol(:j_, d), N)...)
#     bⱼ = Expr(:ref, :B, iⱼ...)
#     c = Expr(:ref, :C, iᵢ..., iⱼ...)
#     δᵢ = Expr(:(=), :Δᵢ, Expr(:call, :-, aᵢ, bᵢ))
#     δⱼ = Expr(:(=), :Δⱼ, Expr(:call, :-, aⱼ, bⱼ))
#     e = Expr(:(=), c, Expr(:call, :+, c, Expr(:call, :*, :Δᵢ, :Δⱼ)))
#     push!(body.args, δᵢ)
#     push!(body.args, δⱼ)
#     push!(body.args, e)
#     body
# end

# function sumsqdiffij_quote(N::Int, D)
#     outer = outerloopgen(N)
#     # outer = outerloopgen(N, D)
#     inner = loopgen(N)
#     body = sumsqdiffij_body(N, D)
#     push!(inner.args, body)
#     push!(outer.args, inner)
#     outer
#     # return quote
#     #     @turbo $outer
#     #     return C
#     # end
# end

# @generated function sumsqdiffij!(C::AbstractArray{Tₒ, Nₒ}, A::AbstractArray{T, N},
#                                  B::AbstractArray{T, N}, dims::D) where {Tₒ, Nₒ, T, N, D}
#     sumsqdiffij_quote(N, D)
# end

# function lvcov(A::AbstractArray{T, N}, dims::NTuple{M, Int}, corrected::Bool) where {T, N, M}
#     B = lvmean(A, dims=dims)
#     Dᴮ = size(B)
#     Dᴮ′ = ntuple(d -> d ∈ dims ? StaticInt(1) : Dᴮ[d], N)
#     Tₒ = Base.promote_op(/, T, Int)
#     C = zeros(Tₒ, Dᴮ..., Dᴮ...)
#     sumsqdiffij!(C, A, B, Dᴮ′)
#     Dᴬ = size(A)
#     denom = one(Tₒ)
#     for d ∈ eachindex(Dᴬ)
#         denom = d ∈ dims ? denom * Dᴬ[d] : denom
#     end
#     denom = corrected ? denom - one(Tₒ) : denom
#     x = inv(denom)
#     smul!(C, x)
#     C
# end


text:
using LinearAlgebra, Statistics, Plots

x_vals = [0 0 0; 2 -3 -4];
y_vals = [0 0 0; 4 3 -3.5];

plot(x_vals, y_vals, arrow = true, color = :blue,
     legend = :none, xlims = (-5, 5), ylims = (-5,5),
     annotations = [(2.2, 4.4, "[2, 4]"),
                    (-3.3, 3.3, "[-3, 3]"),
                    (-4.4, -3.85, "[-4, -3.5]")],
    xticks = -5:1:5, yticks = -5:1:5,
    framestyle = :origin)

 # Illustrate scalar multiplication
 x = [2];
 scalars = [-2 1 2];
 vals = [0 0 0; x * scalars]
 labels = [(-3.6, -4.2, "-2x"), (2.4, 1.8, "x"), (4.4, 3.8, "2x")]

 plot(vals, vals, arrow = true, color = [:red :red :blue],
     legend = :none, xlims = (-5, 5), ylims = (-5, 5),
     annotations = labels, xticks = -5:1:5, yticks = -5:1:5,
     framestyle = :origin)

x = ones(3);
y = [2, 4, 6];
x + y;
4x; # equivalent to 4 * x and 4 .* x

dot(x, y);
norm(x); # Norm of x
sqrt(sum(abs2, x)); # Give the same result

# fixed linear function, to generate a plane
f(x, y) = 0.2x + 0.1y;

# lines to vectors
x_vec = [0 0; 3 3]
y_vec = [0 0; 4 -4]
z_vec = [0 0; f(3, 4) f(3, -4)]

# draw the plane
n = 20
grid = range(-5, 5, length = n)
z2 = [ f(grid[row], grid[col]) for row in 1:n, col in 1:n ]
wireframe(grid, grid, z2, fill = :blues, gridalpha =1 )
plot!(x_vec, y_vec, z_vec, color = [:blue :green], linewidth = 3, labels = "",
      colorbar = false)

A = [1 2
     3 4]
A'
A = ones(3,3);
2I
A + I

# Solve Linear Equations
A = [1.0 2.0; 3.0 4.0];
y = ones(2,1);
det(A);
A_inv = inv(A)
x = A_inv * y # solution
A * x # should equal y (a vector of ones)
A \ y # produces the same solution

A = [1.0 2.0; 2.0 1.0];
evals, evecs = eigen(A);


text:
using Gnuplot
ecycl_x(r,k,θ) = r*(k .+ 1).*cos.(θ) .- r*cos.((k .+ 1) .* θ)
ecycl_y(r,k,θ) = r*(k .+ 1).*sin.(θ) .- r*sin.((k .+ 1) .* θ)
θ = LinRange(0,6.2π,1000)
@gp ecycl_x(1,1,θ) ecycl_y(1,1,θ) "w l lw 2 t '1'"
@gp :- "set size square" title = "Epicycloid"
@gp :- "set key outside title 'k, r=2k' box opaque"
@gp :- xlabel = "x(θ) = r(k+1)cos(θ) -rcos((k+1)θ)"
@gp :- ylabel = "y(θ) = r(k+1)cos(θ) -rcos((k+1)θ) "
for k in 2:0.5:5.
    @gp :- ecycl_x(2k,k,θ) ecycl_y(2k,k,θ) "w l lw 2 t '$(k)'" :-
end
@gp
save(term="pngcairo font 'Consolas, 12' size 600,600", output="./code/plt2_ex6dot1.png") # hide

max_stars_repo_name: lazarusA/gnuplot-examples

max_issues_repo_name: lazarusA/gnuplot-examples

max_forks_repo_name: lazarusA/gnuplot-examples


text:
using DataFrames, CSV
data = CSV.read("../data/purchaseData.csv", copycols = true)

println("Grade of person 1: ", data[1, 3],
        ", ", data[1,:Grade],
        ", ", data.Grade[1], "\n")
println(data[[1,2,4], :], "\n")
println(data[13:15, :Name], "\n")
println(data.Name[13:15], "\n")
println(data[13:15, [:Name]])

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
code_llvm(unsafe_load, Tuple{Ptr{Float32}})

using TinyCassette
struct GPUctx end

code_llvm(TinyCassette.execute,
              Tuple{GPUctx, typeof(unsafe_load), Ptr{Float32}})


text:
# Readme example

x = -2:0.1:2
linear_f(x) = 2x ; quad_f(x) = x^2
y_linear = linear_f.(x)
y_quad = quad_f.(x)

using Plots
scatter(x, y_linear, label="linear")
scatter!(x, y_quad, label="quadratic", legend=:bottomright)
savefig("docs/linquad.png")

using StatsBase

@show(StatsBase.cor(x, y_linear))
@show(StatsBase.cor(x, y_quad))
@show(StatsBase.corspearman(x, y_linear))
@show(StatsBase.corspearman(x, y_quad))

@show(HoeffdingD.hoeffdingd(x, y_linear))
@show(HoeffdingD.hoeffdingd(x, y_quad))
@show(HoeffdingD.hoeffdingd(x, y_linear, 0.05))
@show(HoeffdingD.hoeffdingd(x, y_quad, 0.05))

max_stars_repo_name: ericqu/HoeffdingD

max_issues_repo_name: ericqu/HoeffdingD

max_forks_repo_name: ericqu/HoeffdingD


text:
###
### An abstract biological sequence type.
###
### This file is a part of BioJulia.
### License is MIT: https://github.com/BioJulia/BioSequences.jl/blob/master/LICENSE.md

abstract type BioSequence{A<:Alphabet} end

# Aliases and shorthands for describing subsets of the BioSequence type...
const NucleotideSeq = BioSequence{<:NucleicAcidAlphabet}
const AminoAcidSeq = BioSequence{AminoAcidAlphabet}

###
### Required traits and methods
###

# Base.length must be defined for every T<:BioSequence.

# As must the following...

"Get the length of a biological sequence."
@inline function Base.length(seq::BioSequence)
    error(
        string(
            "Base.length has not been defined for BioSequence type: ",
            typeof(seq),
            ". Any compatible concrete BioSequence subtype must have this method implemented."
        )
    )
end

"""
Return the data member of `seq` that stores the encoded sequence data.
"""
@inline function encoded_data(seq::BioSequence)
    error(
        string(
            "encoded_data has not been defined for BioSequence type: ",
            typeof(seq),
            ". Any compatible concrete BioSequence subtype must have this method implemented."
        )
    )
end

###
### Provided traits and methods
###

# These traits and methods are defined automatically for any subtype of BioSequence{A}.
# They may be overloaded for your concrete BioSequence sub-type if it is nessecery.

"Get the vector of bits storing a sequences packed encoded elements."
@inline encoded_data_type(seq::BioSequence) = typeof(encoded_data(seq))

"Get the element type of the vector of bits storing a sequences packed encoded elements."
@inline encoded_data_eltype(seq::BioSequence) = eltype(encoded_data_type(seq))

"""
Return the `Alpahbet` defining the possible biological symbols
and their encoding for a given biological sequence.
"""
@inline function Alphabet(::Type{<:BioSequence{A}}) where {A <: Alphabet}
    return A()
end

"Return the `Alpahbet` type that defines the biological symbols allowed for `seq`."
@inline function Alphabet(seq::BioSequence)
    return Alphabet(typeof(seq))
end

BioSymbols.alphabet(::Type{BioSequence{A}}) where {A<:Alphabet} = alphabet(A)

BitsPerSymbol(seq::BioSequence) = BitsPerSymbol(Alphabet(seq))

"Get the number of bits each symbol packed into a BioSequence uses, as an integer value."
bits_per_symbol(seq::BioSequence) = bits_per_symbol(Alphabet(seq))

"""
    join(::Type{T <: BioSequence}, seqs)

Concatenate all the `seqs` to a biosequence of type `T`.

# Examples
```julia> join(LongDNASeq, [dna"TAG", mer"AAC"])
"""
Base.join(::Type{T}, seqs) where {T <: BioSequence} = join!(T(), seqs)

function _join!(result::BioSequence, seqs, ::Val{resize}) where resize
    offset = 0
    for seq in seqs
        seqlen = length(seq)
        resize && length(result) < (offset + seqlen) && resize!(result, offset + seqlen)
        result[offset+1:offset+seqlen] = seq
        offset += seqlen
    end
    resize && resize!(result, offset)
    result
end

# This Union is the closest I can get to specifying a stateless iterator,
# such that we can iterate over it twice.
function join!(result::BioSequence, seqs::Union{AbstractArray, Tuple, AbstractSet, AbstractDict}
)
    # This early return is needed, because sum for empty iterables fail
    isempty(seqs) && return resize!(result, 0)
    resize!(result, sum(length, seqs))
    _join!(result, seqs, Val(false))
end

join!(result::BioSequence, seqs) = _join!(result, seqs, Val(true))

# The generic functions for any BioSequence...
include("indexing.jl")
include("conversion.jl")
include("predicates.jl")
include("find.jl")
include("printing.jl")
include("transformations.jl")
include("counting.jl")


text:
# Define a serialization format for the symbolic equations:
mutable struct Node
    #Holds operators, variables, constants in a tree
    degree::Int #0 for constant/variable, 1 for cos/sin, 2 for +/* etc.
    constant::Bool #false if variable
    val::CONST_TYPE
    # ------------------- (possibly undefined below)
    feature::Int #Either const value, or enumerates variable.
    op::Int #enumerates operator (separately for degree=1,2)
    l::Node
    r::Node

    Node(val::CONST_TYPE) =                                               new(0, true,                       val                                     ) #Leave other values undefined
    Node(feature::Int) =                                                  new(0, false, convert(CONST_TYPE, 0f0), feature                            )
    Node(op::Int, l::Node) =                                              new(1, false, convert(CONST_TYPE, 0f0),       0,      op,        l         )
    Node(op::Int, l::Union{CONST_TYPE, Int}) =                            new(1, false, convert(CONST_TYPE, 0f0),       0,      op,  Node(l)         )
    Node(op::Int, l::Node, r::Node) =                                     new(2, false, convert(CONST_TYPE, 0f0),       0,      op,        l,       r)
    Node(op::Int, l::Union{CONST_TYPE, Int}, r::Node) =                   new(2, false, convert(CONST_TYPE, 0f0),       0,      op,  Node(l),       r)
    Node(op::Int, l::Node, r::Union{CONST_TYPE, Int}) =                   new(2, false, convert(CONST_TYPE, 0f0),       0,      op,        l, Node(r))
    Node(op::Int, l::Union{CONST_TYPE, Int}, r::Union{CONST_TYPE, Int}) = new(2, false, convert(CONST_TYPE, 0f0),       0,      op,  Node(l), Node(r))
end

# Copy an equation (faster than deepcopy)
function copyNode(tree::Node)::Node
   if tree.degree == 0
       if tree.constant
           return Node(tree.val)
        else
           return Node(tree.feature)
        end
   elseif tree.degree == 1
       return Node(tree.op, copyNode(tree.l))
    else
        return Node(tree.op, copyNode(tree.l), copyNode(tree.r))
   end
end

# Count the operators, constants, variables in an equation
function countNodes(tree::Node)::Int
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + countNodes(tree.l)
    else
        return 1 + countNodes(tree.l) + countNodes(tree.r)
    end
end

# Count the max depth of a tree
function countDepth(tree::Node)::Int
    if tree.degree == 0
        return 1
    elseif tree.degree == 1
        return 1 + countDepth(tree.l)
    else
        return 1 + max(countDepth(tree.l), countDepth(tree.r))
    end
end

function stringOp(op::F, tree::Node, options::Options;
                  bracketed::Bool=false,
                  varMap::Union{Array{String, 1}, Nothing}=nothing)::String where {F}
    if op in [+, -, *, /, ^]
        l = stringTree(tree.l, options, bracketed=false, varMap=varMap)
        r = stringTree(tree.r, options, bracketed=false, varMap=varMap)
        if bracketed
            return "$l $(string(op)) $r"
        else
            return "($l $(string(op)) $r)"
        end
    else
        l = stringTree(tree.l, options, bracketed=true, varMap=varMap)
        r = stringTree(tree.r, options, bracketed=true, varMap=varMap)
        return "$(string(op))($l, $r)"
    end
end

# Convert an equation to a string
function stringTree(tree::Node, options::Options;
                    bracketed::Bool=false,
                    varMap::Union{Array{String, 1}, Nothing}=nothing)::String
    if tree.degree == 0
        if tree.constant
            return string(tree.val)
        else
            if varMap == nothing
                return "x$(tree.feature)"
            else
                return varMap[tree.feature::Int]
            end
        end
    elseif tree.degree == 1
        return "$(options.unaops[tree.op])($(stringTree(tree.l, options, bracketed=true, varMap=varMap)))"
    else
        return stringOp(options.binops[tree.op], tree, options, bracketed=bracketed, varMap=varMap)
    end
end

# Print an equation
function printTree(tree::Node, options::Options; varMap::Union{Array{String, 1}, Nothing}=nothing)
    println(stringTree(tree, options, varMap=varMap))
end

# Return a random node from the tree
function randomNode(tree::Node)::Node
    if tree.degree == 0
        return tree
    end
    a = countNodes(tree)
    b = 0
    c = 0
    if tree.degree >= 1
        b = countNodes(tree.l)
    end
    if tree.degree == 2
        c = countNodes(tree.r)
    end

    i = rand(1:1+b+c)
    if i <= b
        return randomNode(tree.l)
    elseif i == b + 1
        return tree
    end

    return randomNode(tree.r)
end

# Count the number of unary operators in the equation
function countUnaryOperators(tree::Node)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 1 + countUnaryOperators(tree.l)
    else
        return 0 + countUnaryOperators(tree.l) + countUnaryOperators(tree.r)
    end
end

# Count the number of binary operators in the equation
function countBinaryOperators(tree::Node)::Int
    if tree.degree == 0
        return 0
    elseif tree.degree == 1
        return 0 + countBinaryOperators(tree.l)
    else
        return 1 + countBinaryOperators(tree.l) + countBinaryOperators(tree.r)
    end
end

# Count the number of operators in the equation
function countOperators(tree::Node)::Int
    return countUnaryOperators(tree) + countBinaryOperators(tree)
end

# Count the number of constants in an equation
function countConstants(tree::Node)::Int
    if tree.degree == 0
        if tree.constant
            return 1
        else
            return 0
        end
    elseif tree.degree == 1
        return 0 + countConstants(tree.l)
    else
        return 0 + countConstants(tree.l) + countConstants(tree.r)
    end
end

# Get all the constants from a tree
function getConstants(tree::Node)::AbstractVector{CONST_TYPE}
    if tree.degree == 0
        if tree.constant
            return [tree.val]
        else
            return CONST_TYPE[]
        end
    elseif tree.degree == 1
        return getConstants(tree.l)
    else
        both = [getConstants(tree.l), getConstants(tree.r)]
        return [constant for subtree in both for constant in subtree]
    end
end

# Set all the constants inside a tree
function setConstants(tree::Node, constants::AbstractVector{T}) where {T<:Real}
    if tree.degree == 0
        if tree.constant
            tree.val = convert(CONST_TYPE, constants[1])
        end
    elseif tree.degree == 1
        setConstants(tree.l, constants)
    else
        numberLeft = countConstants(tree.l)
        setConstants(tree.l, constants)
        setConstants(tree.r, constants[numberLeft+1:end])
    end
end


text:
using Gadfly

plot([sin, cos], 0, 25, Scale.color_discrete_manual("red", "purple"))


text:
n = 2 .^ [6 6 7 8 9 10 11 12 13 14 15 16 17 18]

K = 8

MaxGB = 2
MaxGF = 4 * 1.6

m = K .* n
assoc_gbytes = zeros(1,length(n))
assoc_flops = zeros(1,length(n))
assoc_gflops = zeros(1,length(n))
assoc_time = zeros(1,length(n))
assoc_string_gbytes = zeros(1,length(n))
assoc_string_flops = zeros(1,length(n))
assoc_string_gflops = zeros(1,length(n))
assoc_string_time = zeros(1,length(n))

# Integer Labels
println("With Integer Labels")
for i = 1:length(n)

    global assoc_gflops, assoc_time, assoc_gbytes, assoc_flops

    ii = round.(Int, floor.(rand(m[i]) .* n[i]) .+ 1)
    jj = round.(Int, floor.(rand(m[i]) .* n[i]) .+ 1)
    A = Assoc(ii,jj,1.0)

    ii = round.(Int, floor.(rand(m[i]) .* n[i]) .+ 1)
    jj = round.(Int, floor.(rand(m[i]) .* n[i]) .+ 1)
    B = Assoc(ii,jj,1.0)

    assoc_time[i] = @elapsed C = A*B
    assoc_flops[i] = 2*sum(C)
    ii, jj, vv = find(C)
    assoc_gbytes[i] = assoc_gbytes[i] + (length(ii) + length(jj)) + 8 .* m[i] ./ 1e9
    assoc_gflops[i] = assoc_flops[i] ./ assoc_time[i] ./ 1e9
    print("Time: ", assoc_time[i])
    print(", GFlops: ", assoc_gflops[i])
    println(", GBytes: ", assoc_gbytes[i])
end

# String Labels
println("With String Labels")
for i = 1:length(n)

    global assoc_string_gflops, assoc_string_time, assoc_string_gbytes, assoc_string_flops

    ii = join(string.(rand(1:n[i],m[i])),",")*","
    jj = join(string.(rand(1:n[i],m[i])),",")*","
    A = Assoc(ii,jj,1.0)

    ii = join(string.(rand(1:n[i],m[i])),",")*","
    jj = join(string.(rand(1:n[i],m[i])),",")*","
    B = Assoc(ii,jj,1.0)

    assoc_string_time[i] = @elapsed C = A*B

    assoc_string_flops[i] = 2*sum(C)
    ii, jj, vv = find(C)
    assoc_string_gbytes[i] = assoc_string_gbytes[i] + (length(ii) + length(jj)) + 8 .* m[i] ./ 1e9
    assoc_string_gflops[i] = assoc_string_flops[i] ./ assoc_string_time[i] ./ 1e9
    print("Time: ", assoc_string_time[i])
    print(", GFlops: ", assoc_string_gflops[i])
    println(", GBytes: ", assoc_string_gbytes[i])
end


text:
function copy(model::RSDSGEModel)
    return RSDSGEModel([deepcopy(getfield(model, k)) for k = 1:length(fieldnames(model))]...)
end

function copy(meta::Meta)
    return Meta([deepcopy(getfield(meta, k)) for k = 1:length(fieldnames(meta))]...)
end

function copy(parameters::Parameters)
    return Parameters([deepcopy(getfield(parameters, k)) for k = 1:length(fieldnames(parameters))]...)
end

function copy(vars::Variables)
    return Variables([deepcopy(getfield(vars, k)) for k = 1:length(fieldnames(vars))]...)
end

function copy(eqs::Equations)
    return Equations([deepcopy(getfield(eqs, k)) for k = 1:length(fieldnames(eqs))]...)
end

function copy(ss::SteadyState)
    return SteadyState([deepcopy(getfield(ss, k)) for k = 1:length(fieldnames(ss))]...)
end

function copy(trans::TransitionMatrix)
    return TransitionMatrix([deepcopy(getfield(trans, k)) for k = 1:length(fieldnames(trans))]...)
end

function copy(pert::PerturbationSystem)
    return PerturbationSystem([deepcopy(getfield(pert, k)) for k = 1:length(fieldnames(pert))]...)
end

function createpairedtuplelist(keyarray::Array,valuearray::Array)
    out = ()
    for (key,val) in zip(keyarray,valuearray)
        out = (out...,(key,val))
    end
    return out
end

function checknames(names::Array{String})
    # This function throws an error for names that are unacceptable/reserved in sympy (i.e. pi, inf etc)
    badnames = Array{String}(0)
    push!(badnames,"Inf")
    push!(badnames,"inf")
    push!(badnames,"exp")
    push!(badnames,"pi")
    push!(badnames,"log")
    push!(badnames,"beta")

    for name in names
        if any(badnames.==name)
            error("Name: $(name) is a reserved name in SymPy. Please change it.")
        end
    end
    return
end

function sympify_catch(exp::String,eqno::Int)
    try
        return Expr = Sym(exp)
    catch err
        if string(err.val) == "PyObject SympifyError()"
            # This is a terrible way to do this, but I can't figure a better way...
            println("Syntax in equation $(eqno) is invalid.")
            rethrow(err)
        else
            println("Error sympifying equation $(eqno):")
            rethrow(err)
        end
    end
end

function checksymbols(equations::Equations,vars::Variables,shocks::Shocks,parameters::Parameters)
    # This function checks all the expressions for symbols that are unrecognised.
    # Catching this early is helpful for avoiding errors later.
    validsymbols = vcat(vars.leads_sym,shocks.syms,vars.contemps_sym,vars.states_sym,parameters.generic_sym[:])
    for (i,expression) in enumerate(equations.syms)
        symbollist = SymPy.free_symbols(expression)
        for symbol in symbollist
            if !(symbol in validsymbols)
                error("Unknown symbol $(symbol) in equation $(i).")
            end
        end
    end
    return
end


text:
#!/usr/bin/env julia
# -*- coding: UTF-8 -*-
# __julia-version__ = 1.7.0
# __author__        = "Max Kannenberg, Martin Scheidt"
# __copyright__     = "2021"
# __license__       = "ISC"

using TrainRuns, Test

trains = Dict()
paths = Dict()
settings = Dict()

@testset "load data" begin

  println("testing load train data")
  push!(trains, :freight      => @time Train("data/trains/freight.yaml"))
  push!(trains, :local        => @time Train("data/trains/local.yaml"))
  push!(trains, :longdistance => @time Train("data/trains/longdistance.yaml"))

  println("testing load path data")
  push!(paths, :const     => @time Path("data/paths/const.yaml"))
  push!(paths, :slope     => @time Path("data/paths/slope.yaml"))
  push!(paths, :speed     => @time Path("data/paths/speed.yaml"))
  push!(paths, :realworld => @time Path("data/paths/realworld.yaml"))

  println("testing load settings data")
  push!(settings, "default"       => @time Settings())
  push!(settings, "poi"           => @time Settings("data/settings/points_of_interest.yaml"))
  push!(settings, "drivingcourse" => @time Settings("data/settings/driving_course.yaml"))
  push!(settings, "everything"    => @time Settings("data/settings/everything.yaml"))
  push!(settings, "strip"         => @time Settings("data/settings/strip.yaml"))
  push!(settings, "time"          => @time Settings("data/settings/time.yaml"))
  push!(settings, "timestrip"     => @time Settings("data/settings/time_strip.yaml"))
  push!(settings, "velocity"      => @time Settings("data/settings/velocity.yaml"))

  @test typeof(first(paths)[2]) == Path
  @test typeof(first(settings)[2]) == Settings

end

println("====================")

tests = Base.Iterators.product(trains,paths)

## routine to generate the anticipated Dict()
# anticipated = Dict()
# for test in tests
#   println(test[1][1],"-",test[2][1])
#   result = @time trainrun(test[1][2],test[2][2])
#   push!(anticipated, Symbol(String(test[1][1]) * "_" * String(test[2][1])) => result )
# end

anticipated = Dict(
  :default => Dict(
    :longdistance_speed => 499.96109564970516,
    :freight_slope => 831.4768274141168,
    :local_slope => 396.99313307033276,
    :longdistance_const => 328.83479381353095,
    :freight_realworld => 8971.50124080998,
    :longdistance_slope => 329.22915822053164,
    :freight_const => 727.7969403041934,
    :longdistance_realworld => 2900.1198723158523,
    :local_speed => 524.3948201513945,
    :local_realworld => 3443.917823618831,
    :freight_speed => 733.2610572579886,
    :local_const => 392.7234008268302
  )
)

@testset "function trainrun()" begin

  @testset "Default settings" begin

    for test in tests
      test_name = String(test[1][1]) * "_" * String(test[2][1])
      println("testing $test_name")
      @time result = trainrun(test[1][2],test[2][2])
      expected = anticipated[:default][Symbol(test_name)]
      # compare result to test data set
      @test isapprox(result, expected, rtol=0.01)
      println("--------------------")
    end

  end
  println("====================")

end


text:
# type that represent an Edge or a Vertex and contains a reference to another object
abstract Container

# Type used to represent Unspecified values
immutable Unspecified end
const UnspecifiedValue = Unspecified()

function setpropertyvalue!(c::Container, propertykey::String, value::Any, withtracking::Bool = true)
	# set a single property value on the supplied container
	c.attachedproperties[propertykey] = value
	if withtracking
		trackchange(c)
	end
end

function setpropertyvalues!(c::Container, properties::Dict{String, Any} , withtracking::Bool = true)
	# set multiple property values on the supplied container based on a supplied property dictionary
	merge!(c.attachedproperties, properties)
	if withtracking
		trackchange(c)
	end
end

function clearproperties!(c::Container , withtracking::Bool = true)
	c.attachedproperties = Dict{String, Any}()
	if withtracking
		trackchange(c)
	end
end

function get(c::Container, name::String, default::Any=UnspecifiedValue)
	# get the value of a specified property, returning a default value if there is no other property value associated
	return Base.get(c.attachedproperties,name,default)
end

function gettracker(c::Container)
	associatedtracker = UnspecifiedValue

	if isdefined(c, :graph)
		associatedtracker = gettracker(c.graph)
	end

	return associatedtracker
end

function trackchange(c::Container)
	associatedtracker = gettracker(c)

	if associatedtracker != UnspecifiedValue
		trackchange!(associatedtracker,c)
	end
end

function trackadd(c::Container)
	associatedtracker = gettracker(c)

	if associatedtracker != UnspecifiedValue
		trackadd!(associatedtracker,c)
	end
end

function trackremove(c::Container)
	associatedtracker = gettracker(c)

	if associatedtracker != UnspecifiedValue
		trackremove!(associatedtracker,c)
	end
end


text:
a = 2
b = 3
c = 5

mf = BellMF(a, b, c)

# Evaluation tests
@assert Fuzzy.eval(mf, c) == 1
@assert Fuzzy.eval(mf, c - a) == Fuzzy.eval(mf, c + a) == 0.5

# Mean finding tests
@assert Fuzzy.mean_at(mf, 1) == Fuzzy.mean_at(mf, 0.6) == Fuzzy.mean_at(mf, 0.3) == Fuzzy.mean_at(mf, 0) == c


text:
"""
calculates main overlap metrics
based on https://github.com/Visceral-Project/EvaluateSegmentation/blob/master/source/DiceCoefficientMetric.h
https://github.com/Visceral-Project/EvaluateSegmentation/blob/master/source/JaccardCoefficientMetric.h
https://github.com/Visceral-Project/EvaluateSegmentation/blob/master/source/GlobalConsistencyError.h
"""
module MainOverlap
using ..BasicStructs, Parameters

export calculateBAsicOverlap, dice, jaccard, gce

"""
calculate dice coefficient based on precalulated constants
TnTpFpFns - list of basic metrics - in the order [tn,tp,fp,fn ]
res - ResultMetrics struct holding result of all metrics calculated for given run
return modified ResultMetrics with added Dice, Jaccard and global consistency error coefficients
"""
function calculateBAsicOverlap(TnTpFpFns::Vector{Float64} , res::ResultMetrics) ::ResultMetrics
    tn = TnTpFpFns[1]
    tp = TnTpFpFns[2]
    fp = TnTpFpFns[3]
    fn = TnTpFpFns[4]

    return setproperties(res, (dice=  dice(tn,tp,fp) ,jaccard=  jaccard(tn,tp,fp) ,gce =  gce(tn,tp,fp,fn)  ))

end #calculateDice

```@doc
Calculates Dice Coefficient  on the basis of true negative, true positive, false positive and false negatives
 ```
function dice(tp,fp, fn)::Float64
 return  2*tp/(2*tp + fp + fn)
end#dice

```@doc
Calculates Jaccard Coefficient  on the basis of true negative, true positive, false positive and false negatives
```
function jaccard(tp,fp, fn)::Float64
 return  tp/(tp + fp + fn)
end#jaccard

```@doc
Calculates Global Consistency error on the basis of true negative, true positive, false positive and false negatives
```
function gce(tn,tp,fp, fn)::Float64
    n = tn+fp+fn+tp;
    e1 = ( fn*(fn+ 2*tp)/(tp+fn) + fp*(fp + 2*tn)/(tn+fp) )/n;
    e2 = ( fp*(fp+2*tp)/(tp+fp) + fn*(fn + 2*tn)/(tn+fn) )/n;
    return min(e1, e2);
end#dice

end #MainOverlap


text:
export FT, RT, FieldType, RingType

if !@isdefined FieldType
    const FieldType = Float64
end

# only true for cartesian kernels, see homogeneous kernels
const RingType = FieldType
const FT = FieldType
const RT = RingType

if iscxxtype(FT) # define a couple more constructors, conversions, promotions
    # upscaling
    FT(x::AbstractFloat) = FT(Float64(x))
    FT(x::Integer) = FT(Int64(x))
    # have one for Rational
    function FT(x::Rational)
        isinf(x) && return FT(x*Inf) # preserve sign
        convert(FT, x.num)/convert(FT, x.den)
    end

    Base.promote_rule(::Type{<:Union{FT,CxxBaseRef{FT}}}, ::Type{<:Real}) = FT
    Base.promote_rule(::Type{<:Union{FT,CxxBaseRef{FT}}}, ::Type{Real}) = Real

    @cxxdereference (::Type{T})(x::FT) where {T<:Real} = T(float(x))

    @cxxdereference Base.oftype(x::FT, y) = convert(FT, y)
    @cxxdereference Base.isinteger(x::FT) = isinteger(float(x))
else
    for op ∈ (:<, :<=, :>, :>=, :/, :+, :-, :*, :^)
        @eval begin
            Base.$op(x::CxxBaseRef{FT}, y::CxxBaseRef{FT}) = $op(x[], y[]) # solves ambiguity
            Base.$op(x::CxxBaseRef{FT}, y) = $op(x[], y)
            Base.$op(x, y::CxxBaseRef{FT}) = $op(x, y[])
        end
    end
    Base.:+(x::CxxBaseRef{FT}) = +x[]
    Base.:-(x::CxxBaseRef{FT}) = -x[]
end

Base.convert(::Type{FT}, x::CxxBaseRef{FT}) = x[]
Base.convert(::Type{T}, x::CxxBaseRef{FT}) where {T<:Real} = convert(T, x[])

(::Type{T})(x::CxxBaseRef{FT}) where {T<:Real} = T(float(x))
FT(x::CxxBaseRef{FT}) = FT(x[])
Base.isnan(x::CxxBaseRef{FT}) = isnan(x[])

export AffTransformation2,
       AffTransformation3,
       Bbox2,
       Bbox3,
       Circle2,
       Circle3,
       CircularArc2,
       CircularArc3,
       Direction2,
       Direction3,
       IsoRectangle2,
       IsoCuboid3,
       Line2,
       Line3,
       Plane3,
       Point2,
       Point3,
       Ray2,
       Ray3,
       Segment2,
       Segment3,
       Sphere3,
       Tetrahedron3,
       Triangle2,
       Triangle3,
       Vector2,
       Vector3,
       WeightedPoint2,
       WeightedPoint3

const AffTransformation23 = Union{AffTransformation2,AffTransformation3}
const Bbox23 = Union{Bbox2,Bbox3}
const Circle23 = Union{Circle2,Circle3}
const CircularArc23 = Union{CircularArc2,CircularArc3}
const Direction23 = Union{Direction2,Direction3}
const IsoBox23 = Union{IsoRectangle2,IsoCuboid3}
const Line23 = Union{Line2,Line3}
const Point23  = Union{Point2,Point3}
const Ray23 = Union{Ray2,Ray3}
const Segment23 = Union{Segment2,Segment3}
const Triangle23 = Union{Triangle2,Triangle3}
const Vector23 = Union{Vector2,Vector3}
const WeightedPoint23 = Union{WeightedPoint2,WeightedPoint3}

# min and max specialize Base's respective functions
export # Operations
       a,
       approximate_area,
       approximate_squared_length,
       area_divided_by_pi,
       b,
       base1, base2,
       c,
       counterclockwise_in_between,
       d,
       delta,
       diametral_sphere,
       dilate,
       dimension,
       direction,
       dx, dy, dz,
       min_coord, max_coord,
       opposite,
       point,
       projection,
       source,
       squared_length,
       squared_length_divided_by_pi_square,
       supporting_line,
       supporting_plane,
       target,
       to_vector,
       vertex,
       volume,
       x_at_y,
       xmax,
       xmin,
       y_at_x,
       ymax,
       ymin,
       zmax,
       zmin,
       # 2D Conversion
       to_2d,
       to_3d,
       # Predicates
       bounded_side,
       collinear_has_on,
       has_on,
       has_on_boundary,
       has_on_bounded_side,
       has_on_negative_side,
       has_on_positive_side,
       has_on_unbounded_side,
       inverse,
       is_degenerate,
       is_even,
       is_horizontal,
       is_odd,
       is_reflection,
       is_rotation,
       is_scaling,
       is_translation,
       is_vertical,
       is_x_monotone,
       is_y_monotone,
       orientation,
       oriented_side,
       orthogonal_direction,
       orthogonal_transform,
       orthogonal_vector,
       # Access Functions
       cartesian,
       center,
       dimension,
       hm,
       homogeneous,
       hw,
       hx,
       hy,
       hz,
       left,
       m,
       orientation,
       right,
       squared_radius,
       supporting_circle,
       weight,
       x,
       y,
       z,
       # Miscellaneous
       area,
       bbox,
       direction,
       opposite,
       perpendicular,
       perpendicular_line,
       perpendicular_plane,
       to_vector,
       transform,
       vector

_pointfor(x) = _pointfor(typeof(x))
_pointfor(::Type) = Any
for T ∈ (:AffTransformation2, :AffTransformation3
       , :Bbox2, :Bbox3
       , :Circle2, :Circle3
       , :CircularArc2, :CircularArc3
       , :Direction2, :Direction3
       , :IsoRectangle2, :IsoCuboid3
       , :Line2, :Line3
       , :Plane3
       , :Point2, :Point3
       , :Ray2, :Ray3
       , :Segment2, :Segment3
       , :Sphere3
       , :Tetrahedron3
       , :Triangle2, :Triangle3
       , :Vector2, :Vector3
       , :WeightedPoint2, :WeightedPoint3)
    D = Symbol(last(string(T)))
    @eval _pointfor(::Type{<:reference_type_union($T)}) = $(Symbol(:Point, D))
end

## Origin & NullVector
export ORIGIN, NULL_VECTOR

"""
    ORIGIN

A symbolic constant which denotes the point at the origin.

This constant is used in the conversion between points and vectors.

# Example

```jldoctest
julia> p = Point2(1.0, 1.0)
PointC2(1, 1)

julia> v = p - ORIGIN
VectorC2(1, 1)

julia> q = ORIGIN + v
PointC2(1, 1)

julia> p == q
true
```

See also: [`Point2`](@ref), [`Point3`](@ref)
"""
ORIGIN

"""
    NULL_VECTOR

A symbolic constant used to construct zero length vectors.

See also: [`Vector2`](@ref), [`Vector3`](@ref)
"""
NULL_VECTOR

## Enumerations and Related Functions
export Sign,
       NEGATIVE,
       ZERO,
       POSITIVE

export Orientation,
       LEFT_TURN, RIGHT_TURN,
       CLOCKWISE, COUNTERCLOCKWISE,
       COLLINEAR, COPLANAR, DEGENERATE

export OrientedSide
       ON_NEGATIVE_SIDE,
       ON_ORIENTED_BOUNDARY,
       ON_POSITIVE_SIDE

export ComparisonResult,
       SMALLER,
       EQUAL,
       LARGER

export BoundedSide,
       ON_UNBOUNDED_SIDE,
       ON_BOUNDARY,
       ON_BOUNDED_SIDE

export Angle,
       OBTUSE,
       RIGHT,
       ACUTE

export BoxParameterSpace2,
       LEFT_BOUNDARY, RIGHT_BOUNDARY,
       BOTTOM_BOUNDARY, TOP_BOUNDARY,
       INTERIOR, EXTERIOR

"""
    Sign

See also: [`Orientation`](@ref)

|     Enumerator     |
|:-------------------|
| [`NEGATIVE`](@ref) |
| [`ZERO`](@ref)     |
| [`POSITIVE`](@ref) |
"""
Sign

# make Sign look more like a Number
Base.promote_rule(::Type{Sign}, ::Type{T}) where T <: Number = T
Base.promote_rule(::Type{Sign}, ::Type{Bool}) = Sign
(::Type{T})(x::Sign) where T <: AbstractFloat = T(convert(Integer, x))

# reimplement these in julia
Base.:-(x::Sign) = opposite(x)
Base.:*(a::Sign, b::Sign) = reinterpret(Sign, convert(Integer, a) * convert(Integer, b))

"""
    Orientation

See also: [`LEFT_TURN`](@ref), [`RIGHT_TURN`](@ref), [`COLLINEAR`](@ref),
[`CLOCKWISE`](@ref), [`COUNTERCLOCKWISE`](@ref), [`COPLANAR`](@ref)
"""
const Orientation = Sign

"""
    CLOCKWISE = NEGATIVE

See also: [`COUNTERCLOCKWISE`](@ref)
"""
CLOCKWISE

"""
    COLLINEAR = ZERO

See also: [`LEFT_TURN`](@ref), [`RIGHT_TURN`](@ref)
"""
COLLINEAR

"""
    COUNTERCLOCKWISE = POSITIVE

See also: [`CLOCKWISE`]
"""
COUNTERCLOCKWISE

"""
    LEFT_TURN = POSITIVE

See also: [`COLLINEAR`](@ref), [`RIGHT_TURN`](@ref)
"""
LEFT_TURN

"""
    RIGHT_TURN = NEGATIVE

See also: [`COLLINEAR`](@ref), [`LEFT_TURN`](@ref)
"""
RIGHT_TURN

"`COPLANAR = ZERO`" COPLANAR
"`DEGENERATE = ZERO`" DEGENERATE

"""
    OrientedSide

|           Enumerator           |
|:-------------------------------|
| [`ON_NEGATIVE_SIDE`](@ref)     |
| [`ON_ORIENTED_BOUNDARY`](@ref) |
| [`ON_POSITIVE_SIDE`](@ref)     |
"""
const OrientedSide = Sign

"""
    ON_NEGATIVE_SIDE

See also: [`ON_POSITIVE_SIDE`](@ref)
"""
ON_NEGATIVE_SIDE

"""
    ON_ORIENTED__BOUNDARY

See also: [`ON_NEGATIVE_SIDE`](@ref), [`ON_POSITIVE_SIDE`](@ref)
"""
ON_ORIENTED_BOUNDARY

"""
    ON_POSITIVE_SIDE

See also: [`ON_NEGATIVE_SIDE`](@ref)
"""
ON_POSITIVE_SIDE

"""
    opposite(o::OrientedSide)

Returns the opposite side (for example [`ON_POSITIVE_SIDE`](@ref) if `o ==
[ON_NEGATIVE_SIDE](@ref)`, or `[ON_ORIENTED_BOUNDARY](@ref) if o ==
[ON_ORIENTED_BOUNDARY](@ref)`)
"""
opposite(o::OrientedSide)

"""
    ComparisonResult

|    Enumerator     |
|:------------------|
| [`SMALLER`](@ref) |
| [`EQUAL`](@ref)   |
| [`LARGER`](@ref)  |
"""
const ComparisonResult = Sign

"""
    SMALLER = NEGATIVE

See also: [`LARGER`](@ref)
"""
SMALLER

"""
    EQUAL = ZERO

See also: [`SMALLER`](@ref), [`LARGER`](@ref)
"""
EQUAL

"""
    LARGER = POSITIVE

See also: [`SMALLER`](@ref)
"""
LARGER

"""
    BoundedSide

See also: [`opposite(o::BoundedSide)`](@ref)

|         Enumerator          |
|:----------------------------|
| [`ON_UNBOUNDED_SIDE`](@ref) |
| [`ON_BOUNDARY`](@ref)       |
| [`ON_BOUNDED_SIDE`](@ref)   |
"""
BoundedSide

"""
    ON_UNBOUNDED_SIDE = NEGATIVE

See also: [`ON_BOUNDED_SIDE`](@ref)
"""
ON_UNBOUNDED_SIDE

"""
    ON_BOUNDARY = ZERO

See also: [`ON_UNBOUNDED_SIDE`](@ref), [`ON_BOUNDED_SIDE`](@ref)
"""
ON_BOUNDARY

"""
    ON_BOUNDED_SIDE = POSITIVE

See also: [`ON_UNBOUNDED_SIDE`](@ref)
"""
ON_BOUNDED_SIDE

"""
    opposite(o::BoundedSide)

Returns the opposite side (for example [`ON_BOUNDED_SIDE`](@ref) if `o ==
[ON_UNBOUNDED_SIDE](@ref)`, or `[ON_BOUNDARY](@ref) if o ==
[ON_BOUNDARY](@ref)`)
"""
opposite(o::BoundedSide)

"""
    Angle

See also: [`angle`](@ref)

|    Enumerator    |
|:-----------------|
| [`OBTUSE`](@ref) |
| [`RIGHT`](@ref)  |
| [`ACUTE`](@ref)  |
"""
Angle

"""
    BoxParameterSpace2

|        Enumerator         |
|:--------------------------|
| [`LEFT_BOUNDARY`](@ref)   |
| [`RIGHT_BOUNDARY`](@ref)  |
| [`BOTTOM_BOUNDARY`](@ref) |
| [`TOP_BOUNDARY`](@ref)    |
| [`INTERIOR`](@ref)        |
| [`EXTERIOR`](@ref)        |
"""
BoxParameterSpace2

# despite undocumented, might as well just define `opposite` for all of them
for E ∈ (:Sign
       , :Angle
       , :BoundedSide
       , :BoxParameterSpace2)
    @eval opposite(x::$E) = reinterpret($E, -convert(Integer, x))
end

## 2D Kernel Objects
include("kernel/aff_transformation_2.jl")
include("kernel/bbox_2.jl")
include("kernel/circle_2.jl")
include("kernel/circular_arc_2.jl")
include("kernel/direction_2.jl")
include("kernel/iso_rectangle_2.jl")
include("kernel/line_2.jl")
include("kernel/point_2.jl")
include("kernel/ray_2.jl")
include("kernel/segment_2.jl")
include("kernel/triangle_2.jl")
include("kernel/vector_2.jl")
include("kernel/weighted_point_2.jl")

## 3D Kernel Objects
include("kernel/aff_transformation_3.jl")
include("kernel/bbox_3.jl")
include("kernel/circle_3.jl")
include("kernel/circular_arc_3.jl")
include("kernel/direction_3.jl")
include("kernel/iso_cuboid_3.jl")
include("kernel/line_3.jl")
include("kernel/plane_3.jl")
include("kernel/point_3.jl")
include("kernel/ray_3.jl")
include("kernel/segment_3.jl")
include("kernel/sphere_3.jl")
include("kernel/tetrahedron_3.jl")
include("kernel/triangle_3.jl")
include("kernel/vector_3.jl")
include("kernel/weighted_point_3.jl")


text:
### trunction

Base.trunc(dt::NanoDate, p::Type{Year}) = NanoDate(trunc(Date(dt), Year))
Base.trunc(dt::NanoDate, p::Type{Quarter}) = NanoDate(trunc(Date(dt), Quarter))
Base.trunc(dt::NanoDate, p::Type{Month}) = NanoDate(trunc(Date(dt), Month))
Base.trunc(dt::NanoDate, p::Type{Day}) = NanoDate(Date(dt))
Base.trunc(dt::NanoDate, p::Type{Hour}) = dt - Minute(dt) - Second(dt) - Millisecond(dt)
Base.trunc(dt::NanoDate, p::Type{Minute}) = dt - Second(dt) - Millisecond(dt)
Base.trunc(dt::NanoDate, p::Type{Second}) = dt - Millisecond(dt)
Base.trunc(dt::NanoDate, p::Type{Millisecond}) = dt - div(dt.nanosecs, NanosecondsPerMillisecond)
Base.trunc(dt::NanoDate, p::Type{Microsecond}) = dt - rem(dt.nanosecs, NanosecondsPerMicrosecond)
Base.trunc(dt::NanoDate, p::Type{Nanosecond}) = dt


text:
function JacobiP(m::Integer,α::Real,β::Real,z::Real) :: Float64
    sqrt((2*m+α+β+1)/(2^(α+β+1)))*
    sqrt(gamma(m+1)/(gamma(m+α+1)*gamma(m+β+1)*gamma(m+α+β+1)))*
    gamma(m+α+β+1)*jacobip(m, α, β, Float64(z))
end

function GradnJacobiP(m::Integer,α::Real,β::Real,z::Real,n::Integer=1) :: Float64
    if m<n
        return zero(Float64)
    else
        return sqrt((2*m+α+β+1)/(2^(α+β+2*n+1)))*
        sqrt(gamma(m+1)/(gamma(m+α+1)*gamma(m+β+1)*gamma(m+α+β+1)))*
        gamma(m+α+β+n+1)*jacobip(m-n, α+n, β+n, Float64(z))
    end
end


text:
module ExtremumPlanes

const SizeType = UInt16
const VolumeType = UInt32

abstract type RectangularCuboid end

function volume(a :: RectangularCuboid) :: VolumeType
  VolumeType(a.l) * VolumeType(a.w) * VolumeType(a.h)
end

# let us start with no rotation to keep things simple
struct Box <: RectangularCuboid
  l :: SizeType
  w :: SizeType
  h :: SizeType

  function Box(l, w, h)
    @assert l > 0
    @assert w > 0
    @assert h > 0

    new(l, w, h)
  end
end

struct Point
  x :: SizeType
  y :: SizeType
  z :: SizeType

  Point(x, y, z) = new(x, y, z)
end

struct Container <: RectangularCuboid
  l :: SizeType
  w :: SizeType
  h :: SizeType

  boxes :: Vector{Box}
  corners :: Vector{Point}

  Container(l, w, h, boxes = [], corners = [Point(0, 0, 0)]) = new(l, w, h, boxes, corners)
end

end # module ExtremumPlanes


text:
using Flight
using ComponentArrays

#this works conceptually, but to be done properly it requires a tree. a simple
#Node type with label, data, parent and children with a getproperty overload

function define_x0()
    aaa = nothing
    aab = ones(3)
    aba = ones(2)
    abb = ones(4)
    abc = nothing
    baa = ones(3)
    bab = nothing
    bba = nothing
    bbb = nothing
    bca = fill(2, 3)
    bcb = nothing
    bcc = fill(-5, 2)
    aa = (a = aaa, b = aab)
    ab = (a = aba, b = abb, c = abc)
    ba = (a = baa, b = bab)
    bb = (a = bba, b = bbb)
    bc = (a = bca, b = bcb, c = bcc)
    a = (a = aa, b=ab)
    b = (a = ba, b=bb, c=bc)
    root = (a=a, b=b)
    return root
end

function build_dict(node::Union{NamedTuple, Dict})
    # d = Dict{Symbol, Union{Dict, AbstractVector{<:Real}, Nothing}}()
    d = Dict{Symbol, Any}()
    for k in keys(node)
        d[k] = build_dict(node[k])
    end
    return d
end
build_dict(node::Union{AbstractVector{<:Real}, Nothing}) = node

function assemble_x0(node::Dict)
    blocks = Dict{Symbol, AbstractVector{<:Real}}()
    for label in keys(node)
        x_child = assemble_x0(node[label])
        if x_child !== nothing
            blocks[label] = x_child
        end
    end
    return (!isempty(blocks) ? ComponentVector(blocks) : nothing)
end
assemble_x0(node::AbstractVector) = node
assemble_x0(::Nothing) = nothing

function assign_x0!(node::Dict, x::ComponentVector)
    for label in keys(node)
            println("Called, label $label")
        if label in keys(x)
            assign_x0!(node[label], view(x,label))
            node[label] = view(x, label)
        end
    end
end
assign_x0!(::AbstractVector, ::AbstractVector) = nothing

x0_tree = define_x0() |> build_dict
x0 = assemble_x0(x0_tree)
assign_x0!(x0_tree, x0)


text:
"""
    luminance_spec(λ,r)`

Luminance spectrum of type `RSpec`
"""
function luminance_spec(λ,l)
    return LSpec(λ,l)
end

"""
`reflectance_spec(λ,r)`

Reflectance spectrum of type `RSpec`
"""
function reflectance_spec(λ,r)
    return RSpec(λ,r)
end

"""
    transmittance_spec(λ,t)

Tansmittance spectrum at unit thickness (x = 1.0).
"""
function transmittance_spec(λ,t)
    return TSpec(λ,t,1.0)
end

"""
    transmittance_spec(λ, t, x)

Transmittance spectrum with thickness `x` given relative to unit thickness, according to Bougouer’s Law/Lambert’s Law.

See: "Color Measurement" by David L. MacAdam, chapter 3.2: "Bouguer’s Law"
"""
function transmittance_spec(λ::Vector{Real},t::Vector{Real},x::Real)
    return TSpec(λ, t, x)
end

"""
    transmittance_spec(λ,t,d,d0)

Transmittance spectrum with thickness term `d`. Unit thickness is `d0`, according to Bouguer’s Law/Lambert’s Law.

See: "Color Measurement" by David L. MacAdam, chapter 3.2: "Bougouer’s Law"
"""
function transmittance_spec(λ::Vector{Real},t::Vector{Real},d::Real,d0::Real)
    transmittance_spec(λ,t,(d/d0))
end

"""
`reinterpolate(spectrum,λmin,Δλ,λmax,style=:natural)`
"""
function reinterpolate(spectrum,λmin,Δλ,λmax,istyle=:natural,estyle=:linear)
    if typeof(spectrum) == RSpec
        spl = cubicspline(spectrum.λ,spectrum.r,istyle)
        epsl = extrap(spl,)
        extrapolate(spl,xrange,extr::Symbol)
        l,s = interp(spl,(λmin:Δλ:λmax))
        return RSpec(l,s)
    elseif typeof(spectrum) == LSpec
        spl = cubicspline(spectrum.λ,spectrum.l,istyle)
        l,s = interp(spl,(λmin:Δλ,λmax))
        return LSpec(l,s)
    end
end

"""
    adaptated_color(spectrum)
"""
function adapted_color(illuminant1,λrange)
    lum=normalize_spec(blackbody_illuminant(T))
    wp=blackbody_whitepoint(SPECENV, T)
end


text:
"""
    diagview(M) = view(M, diagind(M))

Like `diag(M)` but makes a view.
"""
diagview(A::AbstractMatrix) = view(A, diagind(A))

diagview(A::LinearAlgebra.Diagonal) = A.diag

"""
    storage_type(A)

Return the type of the underlying matrix for PermutedDimsArray, Transpose, etc...,
(e.g., `CuArray{Float64,2}` or `Array{Float64,1}`).
"""
function storage_type(A::AbstractArray)
    P = parent(A)
    typeof(A) === typeof(P) ? typeof(A) : storage_type(P)
end
storage_type(A) = typeof(A)

"""
    B = orient(A, code)

Usually this calls `_orient(A, code)`, which reshapes `A` such that its
nontrivial axes lie in the directions where `code` contains a `:`,
by inserting axes on which `size(B, d) == 1` as needed.

When acting on `A::Transpose`, `A::PermutedDimsArray` etc, it will `collect(A)` first,
because reshaping these is very slow. However, this is only done when the underlying array
is a "normal" CPU `Array`, since e.g., for GPU arrays, `collect(A)` copies the array to
the CPU. Fortunately, the speed penalty for reshaping transposed GPU arrays is lower than
on the CPU.
"""
orient(A::AbstractArray, code::Tuple) = _orient(A::AbstractArray, code::Tuple)

@generated function _orient(A::AbstractArray, code::Tuple)
    list = []
    pretty = [] # just for error
    d = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(A,$d) ))
            push!(pretty, ":")
            d += 1
        else
            push!(list, 1)
            push!(pretty, "*")
        end
    end
    str = join(pretty, ", ")
    d-1 == ndims(A) || throw(ArgumentError(
        "orient(A, ($str)) expeceted $(d-1) dimensions, got ndims(A) = $(ndims(A))"))
    :(reshape(A, ($(list...),)))
end

# performance: transpose faster than reshape? https://github.com/JuliaArrays/LazyArrays.jl/issues/16
orient(A::AbstractVector{T}, ::Tuple{typeof(*),Colon}) where {T <: Number} = PermuteDims(A)
#=
V = rand(500);
@time @reduce W[i] := sum(j,k) V[i]*V[j]*V[k] lazy;
# was 0.140318 seconds seconds, now 0.030957 seconds, factor 4
=#

# performance: avoid reshaping lazy transposes (of CPU arrays), as this is very slow in broadcasting
const LazyPerm = Union{PermutedDimsArray, LinearAlgebra.Transpose, LinearAlgebra.Adjoint}
orient(A::Union{LazyPerm, SubArray{<:Any,<:Any,<:LazyPerm}}, code::Tuple) = begin
    if storage_type(A) <: Union{Array, SArray, Tuple}
        # for "normal" CPU arrays, collect before reshaping
        _orient(collect(A), code)
        # _orient(identity.(A), code) # alternative
    else
        # for other storage (e.g., GPU arrays), call the original algorithm with reshape,
        # since collect copies to CPU. Fortunately, reshape is not slow on GPU
        _orient(A,code)
    end
end
#=
A = rand(10,100); B = rand(100,100); C = rand(100,10);
@time @reduce D[a,d] := sum(b,c) A[a,b] * B[b,c] * C[c,d] lazy;
# was 0.142307 seconds 76.296 MiB, now 0.001855 seconds 10.516 KiB, factor 75
=#

# performance: avoid that if you can just extract A.parent. TODO make @generated perhaps?
const LazyRowVec = Union{
    LinearAlgebra.Transpose{<:Any,<:AbstractVector},
    LinearAlgebra.Adjoint{<:Real,<:AbstractVector} } # NB only Adjoint{Real}
orient(A::LazyRowVec, ::Tuple{typeof(*),Colon,Colon}) = orient(A.parent, (*,*,:))
orient(A::LazyRowVec, ::Tuple{Colon,typeof(*),Colon}) = orient(A.parent, (*,*,:))
orient(A::LazyRowVec, ::Tuple{typeof(*),typeof(*),Colon,Colon}) = orient(A.parent, (*,*,*,:))
orient(A::LazyRowVec, ::Tuple{typeof(*),Colon,typeof(*),Colon}) = orient(A.parent, (*,*,*,:))
orient(A::LazyRowVec, ::Tuple{Colon,typeof(*),typeof(*),Colon}) = orient(A.parent, (*,*,*,:))

using OffsetArrays
# reshaping causes offsets to be forgotten, https://github.com/mcabbott/TensorCast.jl/issues/11
# so explicitly copy them around:
@generated function _orient(A::OffsetArray, code::Tuple)
    off = []
    d = 1
    for s in code.parameters
        if s == Colon
            push!(off, :( A.offsets[$d] ))
            d += 1
        else
            push!(off, 0)
        end
    end
    :( OffsetArray(_orient(parent(A), code), ($(off...),)) )
end

"""
    rview(A, :,1,:) ≈ (@assert size(A,2)==1; view(A, :,1,:))

This simply reshapes `A` so as to remove a trivial dimension where indicated.
Throws an error if size of `A` is not 1 in the indicated dimension.

Will fail silently if given a `(:,3,:)` or `(:,\$c,:)`,
for which `needview!(code) = true`, so the macro should catch such cases.
"""
rview(A::AbstractArray, code...) = rview(A, code)

@generated function rview(A::AbstractArray, code::Tuple)
    list = []
    pretty = [] # just for error
    ex = quote end
    d = 1
    for s in code.parameters
        if s == Colon
            push!(list, :( size(A,$d) ))
            push!(pretty, ":")
            d += 1
        elseif s == Int
            push!(pretty, "1")
            str = "rview(A) expected size(A,$d) == 1"
            push!(ex.args, :(size(A,$d)==1|| throw(DimensionMismatch($str)) ) )
            d += 1
        end
    end
    str = join(pretty, ", ")
    d-1 == ndims(A) || throw(ArgumentError(
        "rview(A, ($str)) expeceted $(d-1) dimensions, got ndims(A) = $(ndims(A))"))
    push!(ex.args, :(reshape(A, ($(list...),))) )
    ex
end

rview(A::LazyRowVec, ::Tuple{Int,Colon}) = A.parent

rview(A::LazyPerm, code::Tuple) = view(A, code)
# cannot use rview(collect(A), code) as this may be on LHS, e.g. sum!(rview(),...)

"""
    rvec(x) = [x]
    rvec(A) = vec(A)

Really-vec... extends `LinearAlgebra.vec` to work on scalars too.
"""
rvec(x::Number) = [x]
rvec(A) = LinearAlgebra.vec(A)

"""
    mul!(Z,A,B)

Exactly `LinearAlgebra.mul!` except that it can write into a zero-array.
"""
mul!(Z,A,B) = LinearAlgebra.mul!(Z,A,B)
mul!(Z::AbstractArray{T,0}, A,B) where {T} = copyto!(Z, A * B)

"""
    star(x,y,...)

Like `*` but intended for multiplying sizes, and understands that `:` is a wildcard.
"""
star(x,y) = *(x,y)
star(::Colon,y) = Colon()
star(x,::Colon) = Colon()
star(x,y,zs...) = star(star(x,y), zs...)

"""
    PermuteDims(A::Matrix)
    PermuteDims(A::Vector)

Lazy like `transpose`, but not recursive:
calls `PermutedDimsArray` unless `eltype(A) <: Number`.
"""
PermuteDims(A::AbstractMatrix) = PermutedDimsArray(A, (2,1))
PermuteDims(A::AbstractMatrix{T}) where {T<:Number} = transpose(A)

PermuteDims(A::AbstractVector) = reshape(A,1,:)
PermuteDims(A::AbstractVector{T}) where {T<:Number} = transpose(A)

struct Reverse{D} end

"""
    Reverse{d}(A)

Lazy version of `reverse(A; dims=d)`. Also constructed by `Reverse(A; dims=d)`.
"""
function Reverse{D}(A::AbstractArray{T,N}) where {D,T,N}
    if D isa Tuple
        length(D)==0 && return A
        return Reverse{Base.tail(D)}(Reverse{first(D)}(A))
    end
    1 ≤ D ≤ N || throw(ArgumentError("dimension $D is not 1 ≤ $D ≤ $N"))
    tup = ntuple(i -> i==D ? reverse(axes(A,D)) : Colon(), N)
    view(A, tup...)
end

Reverse(A::AbstractArray; dims=ndims(A)) = _Reverse(A, dims...)

_Reverse(A::AbstractArray) = A
_Reverse(A::AbstractArray, d::Int, dims::Int...) = _Reverse(Reverse{d}(A), dims...)

Reverse!(A::AbstractArray; dims=ndims(A)) = A .= _Reverse(A, dims...)

# The reverse() function in Base is twice as quick as copying this view,
# somewhere I had a real type in order to exploit this?

using Random

struct Shuffle{D} end

"""
    Shuffle{d}(A)

For a vector this is a lazy version of `shuffle(A)`,
for a tensor this shuffles along one axis only.
"""
function Shuffle{D}(A::AbstractArray{T,N}) where {D,T,N}
    D isa Tuple && return _Shuffle(A, D)
    1 ≤ D ≤ N || throw(ArgumentError("dimension $D is not 1 ≤ $D ≤ $N"))
    tup = ntuple(i -> i==D ? shuffle(axes(A,D)) : Colon(), N)
    view(A, tup...)
end

Shuffle(A::AbstractArray; dims=ndims(A)) = _Shuffle(A, Tuple(dims))

_Shuffle(A::AbstractArray, dims::Tuple{}) = A
_Shuffle(A::AbstractArray, dims::Tuple) = _Shuffle(Shuffle{first(dims)}(A), Base.tail(dims))

Shuffle!(A::AbstractArray; dims=ndims(A)) = A .= _Shuffle(A, dims)


text:
using BinDeps

@BinDeps.setup

libhilbert = library_dependency("libhilbert")

builddir = joinpath(@__DIR__,"..","src","lib")

provides(SimpleBuild,
    (@build_steps begin
        ChangeDirectory( builddir )
        @build_steps begin
            `make`
        end
    end),
libhilbert, installed_libpath=builddir)

@BinDeps.install Dict(
    "libhilbert" => "libhilbert",
)


text:
# unit test the BFunc and bFunc

include("../src/DMD_Util.jl");
include("../src/DMD_Type.jl");
include("../src/DMD_Func.jl");

# Generate DMD Synthetic Data
#------------------------------------------------------------------------------
# dimensions
m = 100;    # temporal dimension
n = 100;    # spatial dimension
k = 3;      # number of modes
T = Float64;
# generate data
X, t, at, Bt = simDMD(m, n, k, T; seed=123456);
# specify the loss
function lossFunc(r)
	T = typeof(real(r[1]));
	return T(0.5)*sum(abs2, r)
end
#
function lossGrad(g, r)
end
#
# Create object
#------------------------------------------------------------------------------
params = DMDParams(k, X, t, lossFunc, lossGrad);
#
# Initialize with true parameter
#------------------------------------------------------------------------------
copyto!(params.a, at);
copyto!(params.B, Bt);
#
# Evaluate BFunc
#------------------------------------------------------------------------------
update_P!(params);
B_obj = BFunc(params);
#
if abs(B_obj) < 1e-6
	println("BFunc: OK")
else
	println("BFunc: Wrong, obj: $B_obj")
end
#
# Evaluate bFunc
#------------------------------------------------------------------------------
Random.randn!(params.R);
id = rand(1:n);
#
b_obj = bFunc(params, id);
#
if abs(b_obj) < 1e-6
	println("bFunc: OK")
else
	println("bFunc: Wrong, obj: $b_obj")
end


text:
# a!/b!
factorial_ratio(a::I, b::I) where {I<:Integer} = gamma(a+1)/gamma(b+1) #binomial(a,b)*factorial(a-b)

function coulomb_analytical(k::I, γ::R, ℓ::I, r̃::R) where {I<:Integer, R<:Real}
    e⁻ᵞʳ = exp(-γ*r̃)

    # S = -gamma(k+ℓ+1)/r̃^(ℓ+1)*(e⁻ᵞʳ-1)/γ^(k+ℓ+1)
    S = 2gamma(k+ℓ+1)*exp(-γ*r̃/2 -(ℓ+1)*log(r̃))*sinh(γ*r̃/2)/γ^(k+ℓ+1)

    for i ∈ 0:k+ℓ-1
        S -= factorial_ratio(k+ℓ,k+ℓ-i)*r̃^(k-i-1)*e⁻ᵞʳ/γ^(i+1)
    end

    for i ∈ 0:k-ℓ-1
        S += factorial_ratio(k-ℓ-1,k-ℓ-i-1)*r̃^(k-i-1)*e⁻ᵞʳ/γ^(i+1)
    end

    S
end


text:
"""
**The Logistic-Softmax likelihood**

The multiclass likelihood with a logistic-softmax mapping: : ``p(y=i|\\{f_k\\}) = \\sigma(f_i)/ \\sum_k \\sigma(f_k)``
where σ is the logistic function has the same properties as softmax.

---

For the analytical version, the likelihood is augmented multiple times to obtain :
```math
#TODO
```
Paper with details under submission
"""
struct LogisticSoftMaxLikelihood{T<:Real} <: MultiClassLikelihood{T}
    Y::AbstractVector{BitVector} #Mapping from instances to classes
    class_mapping::AbstractVector{Any} # Classes labels mapping
    ind_mapping::Dict{Any,Int} # Mapping from label to index
    y_class::AbstractVector{Int64} # GP Index for each sample
    c::AbstractVector{AbstractVector{T}} # Second moment of fₖ
    α::AbstractVector{T} # Variational parameter of Gamma distribution
    β::AbstractVector{T} # Variational parameter of Gamma distribution
    θ::AbstractVector{AbstractVector{T}} # Variational parameter of Polya-Gamma distribution
    γ::AbstractVector{AbstractVector{T}} # Variational parameter of Poisson distribution
    function LogisticSoftMaxLikelihood{T}() where {T<:Real}
        new{T}()
    end
    function LogisticSoftMaxLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class)
    end
    function LogisticSoftMaxLikelihood{T}(Y::AbstractVector{<:BitVector},
    class_mapping::AbstractVector, ind_mapping::Dict{<:Any,<:Int},y_class::AbstractVector{<:Int},
    c::AbstractVector{<:AbstractVector{<:Real}}, α::AbstractVector{<:Real},
    β::AbstractVector, θ::AbstractVector{<:AbstractVector},γ::AbstractVector{<:AbstractVector{<:Real}}) where {T<:Real}
        new{T}(Y,class_mapping,ind_mapping,y_class,c,α,β,θ,γ)
    end
end

function LogisticSoftMaxLikelihood()
    LogisticSoftMaxLikelihood{Float64}()
end

function pdf(l::LogisticSoftMaxLikelihood,f::AbstractVector)
    logisticsoftmax(f)
end

function pdf(l::LogisticSoftMaxLikelihood,y::Integer,f::AbstractVector)
    logisticsoftmax(f)[y]
end

function Base.show(io::IO,model::LogisticSoftMaxLikelihood{T}) where T
    print(io,"Logistic-Softmax Likelihood")
end

function init_likelihood(likelihood::LogisticSoftMaxLikelihood{T},inference::Inference{T},nLatent::Integer,nSamplesUsed::Integer,nFeatures::Integer) where T
    if inference isa AnalyticVI || inference isa GibbsSampling
        c = [ones(T,nSamplesUsed) for i in 1:nLatent]
        α = nLatent*ones(T,nSamplesUsed)
        β = nLatent*ones(T,nSamplesUsed)
        θ = [abs.(rand(T,nSamplesUsed))*2 for i in 1:nLatent]
        γ = [abs.(rand(T,nSamplesUsed)) for i in 1:nLatent]
        LogisticSoftMaxLikelihood{T}(likelihood.Y,likelihood.class_mapping,likelihood.ind_mapping,likelihood.y_class,c,α,β,θ,γ)
    else
        return likelihood
    end
end

## Local Updates Section ##

function local_updates!(model::VGP{T,<:LogisticSoftMaxLikelihood,<:AnalyticVI,V}) where {T,V}
    model.likelihood.c .= broadcast((Σ::V,μ::V)->sqrt.(Σ.+abs2.(μ)),diag.(model.Σ),model.μ)
    for _ in 1:2
        model.likelihood.γ .= broadcast((c::V,μ::V,ψα::V)->0.5/(model.likelihood.β[1])*exp.(ψα).*safe_expcosh.(-0.5*μ,0.5*c),
                                    model.likelihood.c,model.μ,[digamma.(model.likelihood.α)])
        model.likelihood.α .= 1.0.+(model.likelihood.γ...)
    end
    model.likelihood.θ .= broadcast((y::BitVector,γ::V,c::V)->0.5*(y.+γ)./c.*tanh.(0.5*c),model.likelihood.Y,model.likelihood.γ,model.likelihood.c)
end

function local_updates!(model::SVGP{T,<:LogisticSoftMaxLikelihood,<:AnalyticVI,V}) where {T,V}
    model.likelihood.c .= broadcast((μ::V,Σ::Symmetric{T,Matrix{T}},κ::Matrix{T},K̃::V)->sqrt.(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)),
                                    model.μ,model.Σ,model.κ,model.K̃)
    for _ in 1:5
        model.likelihood.γ .= broadcast((c::V,κμ::V,ψα::V)->(0.5/(model.likelihood.β[1]))*exp.(ψα).*safe_expcosh.(-0.5*κμ,0.5*c),
                                    model.likelihood.c,model.κ.*model.μ,[digamma.(model.likelihood.α)])
        model.likelihood.α .= 1.0.+(model.likelihood.γ...)
    end
    model.likelihood.θ .= broadcast((y::BitVector,γ::V,c::V)->0.5*(y[model.inference.MBIndices]+γ)./c.*tanh.(0.5.*c),
                                    model.likelihood.Y,model.likelihood.γ,model.likelihood.c)
    return nothing
end

function sample_local!(model::VGP{T,<:LogisticSoftMaxLikelihood,<:GibbsSampling}) where {T}
    model.likelihood.γ .= broadcast(μ::AbstractVector{<:Real}->rand.(Poisson.(0.5*model.likelihood.α.*safe_expcosh.(-0.5*μ,0.5*μ))), model.μ)
    model.likelihood.α .= rand.(Gamma.(1.0.+(model.likelihood.γ...),1.0./model.likelihood.β))
    pg = PolyaGammaDist()
    model.likelihood.θ .= broadcast((y::BitVector,γ::AbstractVector{<:Real},μ::AbstractVector{<:Real})->draw.([pg],y.+γ,μ),model.likelihood.Y,model.likelihood.γ,model.μ)
    return nothing
end

## Global Gradient Section ##

@inline ∇E_μ(model::VGP{T,<:LogisticSoftMaxLikelihood,<:GibbsorVI}) where {T} = 0.5.*(model.likelihood.Y.-model.likelihood.γ)
@inline ∇E_μ(model::VGP{T,<:LogisticSoftMaxLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5.*(model.likelihood.Y.-model.likelihood.γ)
@inline ∇E_μ(model::SVGP{T,<:LogisticSoftMaxLikelihood,<:GibbsorVI}) where {T} = 0.5.*(getindex.(model.likelihood.Y,[model.inference.MBIndices]).-model.likelihood.γ)
@inline ∇E_μ(model::SVGP{T,<:LogisticSoftMaxLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5.*(model.likelihood.Y[i][model.inference.MBIndices]-model.likelihood.γ[i])
@inline ∇E_Σ(model::AbstractGP{T,<:LogisticSoftMaxLikelihood,<:GibbsorVI}) where {T} = 0.5.*model.likelihood.θ
@inline ∇E_Σ(model::AbstractGP{T,<:LogisticSoftMaxLikelihood,<:GibbsorVI},i::Int) where {T} = 0.5.*model.likelihood.θ[i]

## ELBO Section ##

function ELBO(model::AbstractGP{T,<:LogisticSoftMaxLikelihood,<:AnalyticVI}) where {T}
    return expecLogLikelihood(model) - GaussianKL(model) - GammaEntropy(model) - PoissonKL(model) - PolyaGammaKL(model)
end

function expecLogLikelihood(model::VGP{T,<:LogisticSoftMaxLikelihood,<:AnalyticVI}) where {T}
    # model.likelihood.c .= broadcast((Σ,μ)->sqrt.(Σ.+abs2.(μ)),diag.(model.Σ),model.μ)
    tot = -model.nSample*logtwo
    tot += -sum(sum(model.likelihood.γ.+model.likelihood.Y))*logtwo
    tot +=  0.5*sum(broadcast((y,μ,γ,θ,c)->sum(μ.*(y-γ)-θ.*abs2.(c)),
                    model.likelihood.Y,model.μ,model.likelihood.γ,model.likelihood.θ,model.likelihood.c))
    return tot
end

function expecLogLikelihood(model::SVGP{T,<:LogisticSoftMaxLikelihood,<:AnalyticVI}) where {T}
    model.likelihood.c .= broadcast((μ::AbstractVector,Σ::AbstractMatrix,κ::AbstractMatrix,K̃::AbstractVector)->sqrt.(K̃+opt_diag(κ*Σ,κ)+abs2.(κ*μ)),
                                    model.μ,model.Σ,model.κ,model.K̃)
    tot = -model.inference.nSamplesUsed*logtwo
    tot += -sum(sum.(model.likelihood.γ.+getindex.(model.likelihood.Y,[model.inference.MBIndices])))*logtwo
    tot += 0.5*sum(broadcast((y,κμ,γ,θ,c)->sum((κμ).*(y[model.inference.MBIndices]-γ)-θ.*abs2.(c)),
                    model.likelihood.Y,model.κ.*model.μ,model.likelihood.γ,model.likelihood.θ,model.likelihood.c))
    return model.inference.ρ*tot
end

function PolyaGammaKL(model::VGP{T,<:LogisticSoftMaxLikelihood}) where {T}
    sum(broadcast(PolyaGammaKL,model.likelihood.Y.+model.likelihood.γ,model.likelihood.c,model.likelihood.θ))
end

function PolyaGammaKL(model::SVGP{T,<:LogisticSoftMaxLikelihood}) where {T}
    model.inference.ρ*sum(broadcast(PolyaGammaKL,getindex.(model.likelihood.Y,[model.inference.MBIndices]).+model.likelihood.γ,model.likelihood.c,model.likelihood.θ))
end

function PoissonKL(model::AbstractGP{T,<:LogisticSoftMaxLikelihood}) where {T}
    return model.inference.ρ*sum(broadcast(PoissonKL,model.likelihood.γ,[model.likelihood.α./model.likelihood.β],[digamma.(model.likelihood.α).-log.(model.likelihood.β)]))
end

## Numerical Gradient Section ##

function grad_samples(model::AbstractGP{T,<:LogisticSoftMaxLikelihood,<:NumericalVI},samples::AbstractMatrix{T},index::Integer) where {T}
    class = model.likelihood.y_class[index]::Int64
    grad_μ = zeros(T,model.nLatent)
    grad_Σ = zeros(T,model.nLatent)
    g_μ = similar(grad_μ)
    nSamples = size(samples,1)
    @views @inbounds for i in 1:nSamples
        σ = logistic.(samples[i,:])
        samples[i,:]  .= logisticsoftmax(samples[i,:])
        s = samples[i,class]
        g_μ .= grad_logisticsoftmax(samples[i,:],σ,class)/s
        grad_μ += g_μ
        grad_Σ += diaghessian_logisticsoftmax(samples[i,:],σ,class)/s - abs2.(g_μ)
    end
    for k in 1:model.nLatent
        model.inference.ν[k][index] = -grad_μ[k]/nSamples
        model.inference.λ[k][index] = grad_Σ[k]/nSamples
    end
end

function log_like_samples(model::AbstractGP{T,<:LogisticSoftMaxLikelihood},samples::AbstractMatrix,index::Integer) where {T}
    class = model.likelihood.y_class[index]
    nSamples = size(samples,1)
    loglike = zero(T)
    for i in 1:nSamples
        σ = logistic.(samples[i,:])
        loglike += log(σ[class])-log(sum(σ))
    end
    return loglike/nSamples
end

function grad_logisticsoftmax(s::AbstractVector{<:Real},σ::AbstractVector{<:Real},i::Integer)
    s[i]*(δ.(i,eachindex(σ)).-s).*(1.0.-σ)
end

function diaghessian_logisticsoftmax(s::AbstractVector{<:Real},σ::AbstractVector{<:Real},i::Integer)
    s[i]*(1.0.-σ).*(
    abs2.(δ.(i,eachindex(σ))-s).*(1.0.-σ)
    -s.*(1.0.-s).*(1.0.-σ)
    -σ.*(δ.(i,eachindex(σ))-s))
end

function hessian_logisticsoftmax(s::AbstractVector{T},σ::AbstractVector{T},i::Integer) where {T<:Real}
    m = length(s)
    hessian = zeros(T,m,m)
    @inbounds for j in 1:m
        for k in 1:m
            hessian[j,k] = (1-σ[j])*s[i]*(
            (δ(i,k)-s[k])*(1.0-σ[k])*(δ(i,j)-s[j])
            -s[j]*(δ(j,k)-s[k])*(1.0-σ[k])
            -δ(k,j)*σ[j]*(δ(i,j)-s[j]))
        end
    end
    return hessian
end