text:
using Test, AdvancedHMC, Random
using Statistics: mean
using LinearAlgebra: dot
include("common.jl")

ϵ = 0.01
lf = Leapfrog(ϵ)

θ_init = randn(D)
h = Hamiltonian(UnitEuclideanMetric(D), ℓπ, ∂ℓπ∂θ)
τ = Trajectory{MultinomialTS}(Leapfrog(find_good_stepsize(h, θ_init)), GeneralisedNoUTurn())
r_init = AdvancedHMC.rand(h.metric)

@testset "Passing random number generator" begin
    τ_with_jittered_lf = Trajectory{MultinomialTS}(JitteredLeapfrog(find_good_stepsize(h, θ_init), 1.0), GeneralisedNoUTurn())
    for τ_test in [τ, τ_with_jittered_lf],
        seed in [1234, 5678, 90]
        rng = MersenneTwister(seed)
        z = AdvancedHMC.phasepoint(h, θ_init, r_init)
        z1′ = AdvancedHMC.transition(rng, τ_test, h, z).z

        rng = MersenneTwister(seed)
        z = AdvancedHMC.phasepoint(h, θ_init, r_init)
        z2′ = AdvancedHMC.transition(rng, τ_test, h, z).z

        @test z1′.θ == z2′.θ
        @test z1′.r == z2′.r
    end
end

@testset "TreeSampler" begin
    n_samples = 10_000
    z1 = AdvancedHMC.phasepoint(h, zeros(D), r_init)
    z2 = AdvancedHMC.phasepoint(h, ones(D), r_init)

    rng = MersenneTwister(1234)

    ℓu = rand()
    n1 = 2
    s1 = AdvancedHMC.SliceTS(z1, ℓu, n1)
    n2 = 1
    s2 = AdvancedHMC.SliceTS(z2, ℓu, n2)
    s3 = AdvancedHMC.combine(rng, s1, s2)
    @test s3.ℓu == ℓu
    @test s3.n == n1 + n2

    s3_θ = Vector(undef, n_samples)
    for i = 1:n_samples
        s3_θ[i] = AdvancedHMC.combine(rng, s1, s2).zcand.θ
    end
    @test mean(s3_θ) ≈ ones(D) * n2 / (n1 + n2) rtol=0.01

    w1 = 100
    s1 = AdvancedHMC.MultinomialTS(z1, log(w1))
    w2 = 150
    s2 = AdvancedHMC.MultinomialTS(z2, log(w2))
    s3 = AdvancedHMC.combine(rng, s1, s2)
    @test s3.ℓw ≈ log(w1 + w2)

    s3_θ = Vector(undef, n_samples)
    for i = 1:n_samples
        s3_θ[i] = AdvancedHMC.combine(rng, s1, s2).zcand.θ
    end
    @test mean(s3_θ) ≈ ones(D) * w2 / (w1 + w2) rtol=0.01
end

@testset "TerminationCriterion" begin
    tc = AdvancedHMC.ClassicNoUTurn()
    z1 = AdvancedHMC.phasepoint(h, θ_init, randn(D))
    z2 = AdvancedHMC.phasepoint(h, θ_init, randn(D))
    ts1 = AdvancedHMC.TurnStatistic(tc, z1)
    ts2 = AdvancedHMC.TurnStatistic(tc, z2)
    ts3 = AdvancedHMC.combine(ts1, ts2)
    @test ts1 == ts2 == ts3

    tc = AdvancedHMC.GeneralisedNoUTurn()
    r1 = randn(D)
    z1 = AdvancedHMC.phasepoint(h, θ_init, r1)
    r2 = randn(D)
    z2 = AdvancedHMC.phasepoint(h, θ_init, r2)
    ts1 = AdvancedHMC.TurnStatistic(tc, z1)
    ts2 = AdvancedHMC.TurnStatistic(tc, z2)
    ts3 = AdvancedHMC.combine(ts1, ts2)
    @test ts3.rho == r1 + r2
end

@testset "Termination" begin
    t00 = AdvancedHMC.Termination(false, false)
    t01 = AdvancedHMC.Termination(false, true)
    t10 = AdvancedHMC.Termination(true, false)
    t11 = AdvancedHMC.Termination(true, true)

    @test AdvancedHMC.isterminated(t00) == false
    @test AdvancedHMC.isterminated(t01) == true
    @test AdvancedHMC.isterminated(t10) == true
    @test AdvancedHMC.isterminated(t11) == true

    @test AdvancedHMC.isterminated(t00 * t00) == false
    @test AdvancedHMC.isterminated(t00 * t01) == true
    @test AdvancedHMC.isterminated(t00 * t10) == true
    @test AdvancedHMC.isterminated(t00 * t11) == true

    @test AdvancedHMC.isterminated(t01 * t00) == true
    @test AdvancedHMC.isterminated(t01 * t01) == true
    @test AdvancedHMC.isterminated(t01 * t10) == true
    @test AdvancedHMC.isterminated(t01 * t11) == true

    @test AdvancedHMC.isterminated(t10 * t00) == true
    @test AdvancedHMC.isterminated(t10 * t01) == true
    @test AdvancedHMC.isterminated(t10 * t10) == true
    @test AdvancedHMC.isterminated(t10 * t11) == true

    @test AdvancedHMC.isterminated(t11 * t00) == true
    @test AdvancedHMC.isterminated(t11 * t01) == true
    @test AdvancedHMC.isterminated(t11 * t10) == true
    @test AdvancedHMC.isterminated(t11 * t11) == true
end

@testset "BinaryTree" begin
    z = AdvancedHMC.phasepoint(h, θ_init, randn(D))

    t1 = AdvancedHMC.BinaryTree(z, z, AdvancedHMC.TurnStatistic(), 0.1, 1, -2.0)
    t2 = AdvancedHMC.BinaryTree(z, z, AdvancedHMC.TurnStatistic(), 1.1, 2, 1.0)
    t3 = AdvancedHMC.combine(t1, t2)

    @test t3.sum_α ≈ 1.2 atol=1e-9
    @test t3.nα == 3
    @test t3.ΔH_max == -2.0

    t4 = AdvancedHMC.BinaryTree(z, z, AdvancedHMC.TurnStatistic(), 1.1, 2, 3.0)
    t5 = AdvancedHMC.combine(t1, t4)

    @test t5.ΔH_max == 3.0
end

### Test ClassicNoUTurn and GeneralisedNoUTurn

function makeplot(
    plt,
    traj_θ,
    ts_list...
)
    function plotturn!(traj_θ, ts)
        s = 9.0
        idcs_nodiv = ts .== false
        idcs_div = ts .== true
        idcs_nodiv[1] = idcs_div[1] = false # avoid plotting the first point
        plt.scatter(traj_θ[1,idcs_nodiv], traj_θ[2,idcs_nodiv], s=s, c="black",  label="¬div")
        plt.scatter(traj_θ[1,idcs_div], traj_θ[2,idcs_div],  s=s, c="red",    label="div")
        plt.scatter(traj_θ[1,1], traj_θ[2,1], s=s, c="yellow", label="init")
    end

    fig = plt.figure(figsize=(16, 3))

    for (i, ts, title) in zip(
        1:length(ts_list),
        ts_list,
        [
            "Hand original (v = 1)",
            "AHMC original (v = 1)",
            "Hand generalised (v = 1)",
            "AHMC generalised (v = 1)"
        ]
    )
        plt.subplot(1, 4, i)
        plotturn!(traj_θ, ts)
        plt.gca().set_title(title)
        plt.legend()
    end

    return fig
end

function gettraj(rng, ϵ=0.1, n_steps=50)
    lf = Leapfrog(ϵ)

    q_init = randn(rng, D)
    p_init = AdvancedHMC.rand(rng, h.metric)
    z = AdvancedHMC.phasepoint(h, q_init, p_init)

    traj_z = Vector(undef, n_steps)
    traj_z[1] = z
    for i = 2:n_steps
        traj_z[i] = AdvancedHMC.step(lf, h, traj_z[i-1])
    end

    return traj_z
end

function hand_isturn(z0, z1, rho, v=1)
    θ0minusθ1 = z0.θ - z1.θ
    s = (dot(-θ0minusθ1, -z0.r) >= 0) || (dot(θ0minusθ1, z1.r) >= 0)
    return s
end

ahmc_isturn(z0, z1, rho, v=1) =
    AdvancedHMC.isterminated(ClassicNoUTurn(), h, AdvancedHMC.BinaryTree(z0, z1, AdvancedHMC.TurnStatistic(), 0, 0, 0.0)).dynamic

function hand_isturn_generalised(z0, z1, rho, v=1)
    s = (dot(rho, -z0.r) >= 0) || (dot(-rho, z1.r) >= 0)
    return s
end

ahmc_isturn_generalised(z0, z1, rho, v=1) =
    AdvancedHMC.isterminated(GeneralisedNoUTurn(), h, AdvancedHMC.BinaryTree(z0, z1, AdvancedHMC.TurnStatistic(rho), 0, 0, 0.0)).dynamic

function ahmc_isturn_strictgeneralised(z0, z1, rho, v=1)
    t = AdvancedHMC.isterminated(
        StrictGeneralisedNoUTurn(),
        h,
        AdvancedHMC.BinaryTree(z0, z1, AdvancedHMC.TurnStatistic(rho), 0, 0, 0.0),
        AdvancedHMC.BinaryTree(z0, z0, AdvancedHMC.TurnStatistic(rho - z1.r), 0, 0, 0.0),
        AdvancedHMC.BinaryTree(z1, z1, AdvancedHMC.TurnStatistic(rho - z0.r), 0, 0, 0.0)
    )
    return t.dynamic
end

"""
Check whether the subtree checks adequately detect U-turns.
"""
function check_subtree_u_turns(z0, z1, rho)
    t = AdvancedHMC.BinaryTree(z0, z1, AdvancedHMC.TurnStatistic(rho), 0, 0, 0.0)
    # The left and right subtree are created in such a way that the
    # check_left_subtree and check_right_subtree checks should be equivalent
    # to the general no U-turn check.
    tleft = AdvancedHMC.BinaryTree(z0, z0, AdvancedHMC.TurnStatistic(rho - z1.r), 0, 0, 0.0)
    tright = AdvancedHMC.BinaryTree(z1, z1, AdvancedHMC.TurnStatistic(rho - z0.r), 0, 0, 0.0)

    s1 = AdvancedHMC.isterminated(GeneralisedNoUTurn(), h, t)
    s2 = AdvancedHMC.check_left_subtree(h, t, tleft, tright)
    s3 = AdvancedHMC.check_right_subtree(h, t, tleft, tright)
    @test s1 == s2 == s3
end

@testset "ClassicNoUTurn" begin
    n_tests = 4
    for _ = 1:n_tests
        seed = abs(rand(Int8) + 128)
        rng = MersenneTwister(seed)
        @testset "seed = $seed" begin
            traj_z = gettraj(rng)
            traj_θ = hcat(map(z -> z.θ, traj_z)...)
            traj_r = hcat(map(z -> z.r, traj_z)...)
            rho = cumsum(traj_r, dims=2)

            ts_hand_isturn_fwd = hand_isturn.(Ref(traj_z[1]), traj_z, [rho[:,i] for i = 1:length(traj_z)], Ref(1))
            ts_ahmc_isturn_fwd = ahmc_isturn.(Ref(traj_z[1]), traj_z, [rho[:,i] for i = 1:length(traj_z)], Ref(1))

            ts_hand_isturn_generalised_fwd = hand_isturn_generalised.(Ref(traj_z[1]), traj_z, [rho[:,i] for i = 1:length(traj_z)], Ref(1))
            ts_ahmc_isturn_generalised_fwd = ahmc_isturn_generalised.(Ref(traj_z[1]), traj_z, [rho[:,i] for i = 1:length(traj_z)], Ref(1))

            ts_ahmc_isturn_strictgeneralised_fwd = ahmc_isturn_strictgeneralised.(Ref(traj_z[1]), traj_z, [rho[:,i] for i = 1:length(traj_z)], Ref(1))

            check_subtree_u_turns.(Ref(traj_z[1]), traj_z, [rho[:,i] for i = 1:length(traj_z)])

            @test ts_hand_isturn_fwd[2:end] ==
                ts_ahmc_isturn_fwd[2:end] ==
                ts_hand_isturn_generalised_fwd[2:end] ==
                ts_ahmc_isturn_generalised_fwd[2:end] ==
                ts_ahmc_isturn_strictgeneralised_fwd[2:end]

            if length(ARGS) > 0 && ARGS[1] == "--plot"
                import PyPlot
                fig = makeplot(
                    PyPlot,
                    traj_θ,
                    ts_hand_isturn_fwd,
                    ts_ahmc_isturn_fwd,
                    ts_hand_isturn_generalised_fwd,
                    ts_ahmc_isturn_generalised_fwd
                )
                fig.savefig("seed=$seed.png")
            end
        end
    end
end


text:
module IntervalBasedForest

using DecisionTree: apply_tree_proba, build_tree, Node, prune_tree
using Statistics, StableRNGs

"""
   `TimeSeriesForestClassifier(m, Xmatrix::Array, yplain::Array)`

TimeSeriesForestClassifier builds a forest of trees from the training set
`(Xmatrix, yplain)` using random intervals and summary features.
"""
function TimeSeriesForestClassifier(m, Xmatrix::Array, yplain::Array)
    n_trees = m.n_trees
    Interval_features, Intervals = InvFeatureGen(Xmatrix, n_trees, m.min_interval,
                                                 m.random_state)
    forest = Array{Node,1}()
    for i in range(1, stop=n_trees)
        tree = build_tree(yplain, Interval_features[i],
                         m.n_subfeatures,
                         m.max_depth,
                         m.min_samples_leaf,
                         m.min_samples_split,
                         m.min_purity_increase)
        if m.post_prune
            tree = prune_tree(tree, m.merge_purity_threshold)
        end
        push!(forest, tree)
    end
    forest, Intervals
end

"""
   `predict_new(Xmatrix, forest, Intervals::Array, integers_seen)`

Find predictions for all cases in Xmatrix. Built on top of `apply_tree_proba`.
"""
function predict_new(Xmatrix::Array, forest, Intervals::Array, integers_seen)
    n_instance, s_length = size(Xmatrix)
    n_trees = length(forest)
    Interval_features = InvFeatureGen(Xmatrix, Intervals, n_trees)
    sum = zeros(n_instance, length(integers_seen))
    for i=1:n_trees
        sum += apply_tree_proba(forest[i], Interval_features[i], integers_seen)
    end
    return sum/n_trees
end

"""
   `InvFeatureGen(Xmatrix::Array, n_trees::Int, min_interval::Int, random_state)`

InvFeatureGen selects random intervals, generating the mean, standard deviation
and slope of the random intervals resulting 3*√m features.
"""
function InvFeatureGen(Xmatrix::Array, n_trees::Int, min_interval::Int, random_state)
    n_instance, series_length = size(Xmatrix)
    Interval_features = Array{Array{Float64,2},1}()
    n_intervals = floor(Int, sqrt(series_length))
    Intervals = zeros(Int, n_trees, n_intervals, 2)
    for i in range(1, stop = n_trees)
       interval_feature = Array{Float64,2}(undef, 3*n_intervals, n_instance)
       for j in range(1, stop = n_intervals)
           if typeof(random_state) == StableRNGs.LehmerRNG
               Intervals[i,j,1] = rand(random_state, 1:(series_length - min_interval))
               len = rand(random_state, 1:(series_length - Intervals[i,j,1]))
           else
               Intervals[i,j,1] = rand(1:(series_length - min_interval))
               len = rand(1:(series_length - Intervals[i,j,1]))
           end
           if len < min_interval
               len = min_interval
           end
           Intervals[i,j,2] = Intervals[i,j,1] + len
           Y = Xmatrix[:, Intervals[i,j,1]:Intervals[i,j,2]]
           x = Array(1:size(Y)[2])
           means = mean(Y, dims=2)
           stds =  std(Y, dims=2)
           slope = (mean(transpose(x).*Y, dims=2) -
                    mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
           interval_feature[3*j-2,:] =  means
           interval_feature[3*j-1,:] =  stds
           interval_feature[3*j,:]   =  slope
       end
           push!(Interval_features, transpose(interval_feature))
    end
    return Interval_features, Intervals
end

function InvFeatureGen(Xmatrix::Array, Intervals::Array, n_trees::Int)
    n_instance, series_length = size(Xmatrix)
    Interval_features = Array{Array{Float64,2},1}()
    n_intervals = floor(Int, sqrt(series_length))
    for i in range(1, stop = n_trees)
       interval_feature = Array{Float64,2}(undef, 3*n_intervals, n_instance)
       for j in range(1, stop = n_intervals)
           Y = Xmatrix[:, Intervals[i,j,1]:Intervals[i,j,2]]
           x = Array(1:size(Y)[2])
           means = mean(Y, dims=2)
           stds =  std(Y, dims=2)
           slope = (mean(transpose(x).*Y, dims=2) -
                    mean(x)*mean(Y, dims=2)) / (mean(x.*x) - mean(x)^2)
           interval_feature[3*j-2,:] =  means
           interval_feature[3*j-1,:] =  stds
           interval_feature[3*j,:]   =  slope
       end
           push!(Interval_features, transpose(interval_feature))
    end
    return Interval_features
end

end # TimeSeriesForestClassifier

ext: jl


text:
using Documenter
using Plots
using ClustForOpt

makedocs(sitename="ClustForOpt.jl",
    authors = "Holger Teichgraeber, and Elias Kuepper",
    pages = [
        "Introduction" => "index.md",
        "Workflow" => "workflow.md",
        "Load Data" => "load_data.md",
        "Clustering" => "clust.md",
        "Optimization" => "opt.md"
        ],
    assets = [
        "assets/clust_for_opt_text.svg",
        "assets/opt_cep.svg",
        "assets/workflow.svg"])

deploydocs(repo = "github.com/holgerteichgraeber/ClustForOpt.jl.git", devbranch = "dev")


text:
# Apply boundary conditions for FV

function apply_boundary_conditions_to_face_rhs(var, fid, facefluxvec, t)
    multivar = typeof(var) <: Array;
    maxvarindex = 0;
    if multivar
        # multiple variables being solved for simultaneously
        dofs_per_node = 0;
        var_to_dofs = [];
        for vi=1:length(var)
            tmp = dofs_per_node;
            dofs_per_node += length(var[vi].symvar);
            push!(var_to_dofs, (tmp+1):dofs_per_node);
            maxvarindex = max(maxvarindex,var[vi].index);
        end
    else
        # one variable
        dofs_per_node = length(var.symvar);
        maxvarindex = var.index;
    end

    # Boundary conditions are applied to flux
    fbid = grid_data.facebid[fid]; # BID of this face
    if fbid > 0
        facex = grid_data.allnodes[:, grid_data.face2glb[:,1,fid]];  # face node coordinates

        if multivar
            dofind = 0;
            for vi=1:length(var)
                for compo=1:length(var[vi].symvar)
                    dofind = dofind + 1;
                    if prob.bc_type[var[vi].index, fbid] == NO_BC
                        # do nothing
                    elseif prob.bc_type[var[vi].index, fbid] == FLUX
                        # compute the value and add it to the flux directly
                        # Qvec = (refel.surf_wg[grid_data.faceRefelInd[1,fid]] .* geo_factors.face_detJ[fid])' * (refel.surf_Q[grid_data.faceRefelInd[1,fid]])[:, refel.face2local[grid_data.faceRefelInd[1,fid]]]
                        # Qvec = Qvec ./ geo_factors.area[fid];
                        # bflux = FV_flux_bc_rhs_only(prob.bc_func[var[vi].index, fbid][compo], facex, Qvec, t, dofind, dofs_per_node) .* geo_factors.area[fid];

                        bflux = FV_flux_bc_rhs_only_simple(prob.bc_func[var[vi].index, fbid][compo], facex, t) .* geo_factors.area[fid];

                        fluxvec[(eid-1)*dofs_per_node + dofind] += (bflux - facefluxvec[(fid-1)*dofs_per_node + dofind]) ./ geo_factors.volume[eid];
                        facefluxvec[(fid-1)*dofs_per_node + dofind] = bflux;
                    else
                        printerr("Unsupported boundary condition type: "*prob.bc_type[var[vi].index, fbid]);
                    end
                end
            end
        else
            for d=1:dofs_per_node
                dofind = d;
                if prob.bc_type[var.index, fbid] == NO_BC
                    # do nothing
                elseif prob.bc_type[var.index, fbid] == FLUX
                    # compute the value and add it to the flux directly
                    # Qvec = (refel.surf_wg[grid_data.faceRefelInd[1,fid]] .* geo_factors.face_detJ[fid])' * (refel.surf_Q[grid_data.faceRefelInd[1,fid]])[:, refel.face2local[grid_data.faceRefelInd[1,fid]]]
                    # Qvec = Qvec ./ geo_factors.area[fid];
                    # bflux = FV_flux_bc_rhs_only(prob.bc_func[var.index, fbid][d], facex, Qvec, t, dofind, dofs_per_node) .* geo_factors.area[fid];
                    bflux = FV_flux_bc_rhs_only_simple(prob.bc_func[var.index, fbid][d], facex, t) .* geo_factors.area[fid];

                    fluxvec[(eid-1)*dofs_per_node + dofind] += (bflux - facefluxvec[(fid-1)*dofs_per_node + dofind]) ./ geo_factors.volume[eid];
                    facefluxvec[(fid-1)*dofs_per_node + dofind] = bflux;
                else
                    printerr("Unsupported boundary condition type: "*prob.bc_type[var.index, fbid]);
                end
            end
        end
    end# BCs

end

function FV_flux_bc_rhs_only(val, facex, Qvec, t=0, dofind = 1, totaldofs = 1)
    if typeof(val) <: Number
        return val;

    elseif typeof(val) == Coefficient && typeof(val.value[1]) == GenFunction
        bvals = zeros(size(facex,2));
        if config.dimension == 1
            for i=1:length(bvals)
                bvals[i]=val.value[1].func(facex[1,i],0,0,t);
            end
        elseif config.dimension == 2
            for i=1:length(bvals)
                bvals[i]=val.value[1].func(facex[1,i],facex[2,i],0,t);
            end
        else
            for i=1:length(bvals)
                bvals[i]=val.value[1].func(facex[1,i],facex[2,i],facex[3,i],t);
            end
        end

    elseif typeof(val) == GenFunction
        bvals = zeros(size(facex,2));
        if config.dimension == 1
            for i=1:length(bvals)
                bvals[i]=val.func(facex[1,i],0,0,t);
            end
        elseif config.dimension == 2
            for i=1:length(bvals)
                bvals[i]=val.func(facex[1,i],facex[2,i],0,t);
            end
        else
            for i=1:length(bvals)
                bvals[i]=val.func(facex[1,i],facex[2,i],facex[3,i],t);
            end
        end
    end

    # Do quadrature over the face
    b = Qvec * bvals;

    return b[1];
end

function FV_flux_bc_rhs_only_simple(val, bctype, facex, t=0)
    if typeof(val) <: Number
        return val;

    elseif typeof(val) == Coefficient && typeof(val.value[1]) == GenFunction
        if config.dimension == 1
            bval = val.value[1].func(facex[1],0,0,t);
        elseif config.dimension == 2
            bval = val.value[1].func(facex[1],facex[2],0,t);
        else
            bval = val.value[1].func(facex[1],facex[2],facex[3],t);
        end

    elseif typeof(val) == GenFunction
        if config.dimension == 1
            bval = val.func(facex[1],0,0,t);
        elseif config.dimension == 2
            bval = val.func(facex[1],facex[2],0,t);
        else
            bval = val.func(facex[1],facex[2],facex[3],t);
        end
    end

    return bval;
end


text:
using StructViews

using Test

struct Point
    x
    y::Int
end

function create_point_array()
    points = Array{Point}(undef, 20, 30)

    for i in 1:length(points)
        points[i] = Point(i, i+100)
    end

    points
end

function create_point_array_and_view()
    points = create_point_array()
    view = StructView(points)

    return points, view
end

function testarrayssame(a, b)
    @test a == b
    @test size(a) == size(b)
    @test length(a) == length(b)
    @test eltype(a) === eltype(b)
    @test axes(a) === axes(b)

    if(isempty(a))
        return
    end

    @test all(a .=== b)
    @test all([a[i] === b[i] for i in 1:length(a)])
    @test all(all(a[i, :] .=== b[i, :]) for i in 1:20)
    @test all(all(eltype.(a[i, :]) .=== eltype.(b[i, :])) for i in 1:20)
end

@testset "basic_wrapper" begin
    points, view = create_point_array_and_view()
    @test points[2, 2] === view[2, 2]

    testarrayssame(points, view)
    @test parent(view) === points
    @test :x in propertynames(view)
    @test :y in propertynames(view)
    @test :parent in propertynames(view)
end

@testset "basic_field" begin
    points, view = create_point_array_and_view()

    X = (point->point.x).(points)
    Y = (point->point.y).(points)

    # X may contain by Any type
    X = convert(Array{Any}, X)

    testarrayssame(X, view.x)
    testarrayssame(Y, view.y)

    @test parent(view.x) === points
    @test parent(view.y) === points
end

@testset "field_view" begin
    points = create_point_array()

    X = (point->point.x).(points)
    Y = (point->point.y).(points)
    xview = FieldView{Int, :x}(points)
    yview = FieldView{:y}(points)
    testarrayssame(X, xview)
    testarrayssame(Y, yview)
    @test X[2, 2] === xview[2, 2]
    @test parent(xview) === points
    @test parent(yview) === points
end

struct Point1
    x::Complex{Int}
    y::Int
end

function create_point1_array()
    points = Array{Point1}(undef, 20, 30)

    for i in 1:length(points)
        points[i] = Point1(i+2im, i+100)
    end

    points
end

@testset "recursive_view" begin
    points = create_point1_array()

    X = (point->point.x).(points)
    Y = (point->point.y).(points)

    view = StructView(points)

    xview = view.x
    yview = view.y

    @test eltype(xview) === eltype(X)
    @test eltype(yview) === eltype(Y)

    testarrayssame(X, xview)
    testarrayssame(Y, yview)
    @test X[2, 2] === xview[2, 2]

    Xim = (point->point.x.im).(points)
    Xre = (point->point.x.re).(points)

    @test !(:re in propertynames(yview))

    @test :im in propertynames(xview)
    @test :re in propertynames(xview)

    ximview = xview.im
    xreview = xview.re
    @test eltype(ximview) === eltype(Xim)
    @test eltype(xreview) === eltype(Xre)

    testarrayssame(Xim, ximview)
    testarrayssame(Xre, xreview)

    @test Xim[2, 2] === ximview[2, 2]
end

struct H{T}
    t::T
    f::Int
end

struct Point2
    x::Complex
    y::Int
    p::Point2
    h::H{Point2}
end

function create_point2_array()
    Array{Point2}(undef, 0, 0)
end

@testset "recursive_view_with_cycle" begin
    points = create_point2_array()

    X = (point->point.x).(points)
    Y = (point->point.y).(points)

    view = StructView(points)

    @test !(:p in propertynames(view))
    @test !(:t in propertynames(view.h))
    @test (:f in propertynames(view.h))

    xview = view.x
    yview = view.y

    @test eltype(xview) === eltype(X)
    @test eltype(yview) === eltype(Y)

    testarrayssame(X, xview)
    testarrayssame(Y, yview)

    Xim = (point->point.x.im).(points)
    Xre = (point->point.x.re).(points)

    @test !(:re in propertynames(yview))

    @test :im in propertynames(xview)
    @test :re in propertynames(xview)

    ximview = xview.im
    xreview = xview.re
    @test eltype(ximview) === eltype(Xim)
    @test eltype(xreview) === eltype(Xre)

    testarrayssame(Xim, ximview)
    testarrayssame(Xre, xreview)
end

mutable struct MitableComplex
    im
    re
end

mutable struct Point3
    x
    y::MitableComplex
end

function create_point3_array()
    points = Array{Point3}(undef, 20, 30)

    for i in 1:length(points)
        points[i] = Point3(i, MitableComplex(i+1, i+100))
    end

    points
end

@testset "update_mutable" begin
    points = create_point3_array()
    view = StructView(points)

    X = convert(Array{Any}, (point->point.x).(points))

    testarrayssame(X, view.x)

    view.x .+= 1
    X .+= 1

    testarrayssame(X, view.x)

    view.x[1] = 123

    @test points[1].x === 123

    c = MitableComplex(1, 1)
    view.y[10] = c

    @test points[10].y === c

    view.y.im[11] = -10

    @test points[11].y.im === -10

    view.y.im[2, 2] = -11
    @test points[2, 2].y.im === -11
end

@testset "indexstyle" begin
    a = [1, 2, 3]
    @test IndexStyle(a) === IndexStyle(StructView(a))
    points = create_point_array()
    view = StructView(points)
    @test IndexStyle(points) === IndexStyle(view)
    @test IndexStyle(points) === IndexStyle(view.x)
end

@testset "update_immutable" begin
    points = create_point1_array()
    view = StructView(points)

    Y = (point->point.y).(points)

    view.y .+= 1

    Y .+= 1

    testarrayssame(Y, view.y)

    point = points[1]

    view.y .+= 1

    @test point !== points[1]
    @test point != points[1]

    view.y .-= 1

    @test point === points[1]
    @test point == points[1]

    point2 = points[2, 2]

    view.x.im .+= 1

    @test point2 !== points[2, 2]
    @test point2.x.im + 1 === points[2, 2].x.im

    view.x.im .-= 1

    @test point2 === points[2, 2]

    view.x.im[2, 2] = 123
    @test points[2, 2].x.im === 123
end

@testset "view_updated" begin
    a = [1, 2, 3]
    s = StructView(a)

    push!(s, 4)

    @test a == [1, 2, 3, 4]

    @test pop!(s) === 4

    @test a == [1, 2, 3]

    append!(s, [1, 2, 3])

    @test a == [1, 2, 3, 1, 2, 3]
end


text:
module MyGA
path = @__DIR__
include(joinpath(path,"Turtle.jl"))
include(joinpath(path,"DesitionTree.jl"))
include(joinpath(path,"Cannonball.jl"))
end # module


text:
"""
    trim(model, xu_trim_guess, xu_xdot_trimmed, i_xu_xdot_trimmed, lxu, uxu;
        derivatives = SNOW.ComplexStep(),
        solver = SNOW.IPOPT(),
        objective = (g, xu) -> LA.norm(g, Inf)
    )

Obtains a trimmed state of the model.

# Arguments

* `model<:AbstractModel` - model to be trimmed
* `xu_trim_guess::Vector{Float64}` - guess as to the concatenated state and controls when trimmed
* `xdot_trimmed::Vector{Float64}` - values taken on by any number of the time derivative of the state when trimmed
* `i_xdot_trimmed::Vector{Int}` - indices of the time derivative of the state corresponding to `xdot_trimmed` (these must be the same length)
* `lxu::Vector{Float64}` - lower bounds of the concatenated state and control vector
* `uxu::Vector{Float64}` - upper bounds of the concatenated state and control vector

# Keyword Arguments

* `derivatives = SNOW.ComplexStep()`
* `solver = SNOW.IPOPT()`
* `objective = (g, xu) -> LA.norm(g, Inf)`

"""
function trim(model, xu_trim_guess, x_xdot_trimmed, i_x_xdot_trimmed, lxu, uxu;
    derivatives = SNOW.ComplexStep(),
    solver = SNOW.IPOPT(),
    objective = (x, xdot, u) -> xdot[7]
)
    @assert length(x_xdot_trimmed) == length(i_x_xdot_trimmed) "length of xu_xdot_trimmed and i_xu_xdot_trimmed inconsistent"
    @assert length(xu_trim_guess) == model.nx + model.nu "length of xu_trim_guess must be the sum of states and controls nx + nu"

    # prepare optimizer options
    options = SNOW.Options(derivatives = derivatives, solver = solver)
    ng = length(x_xdot_trimmed) * 2
    lg = fill(-Inf, ng)
    # lg = zeros(ng) # equality constraint
    ug = zeros(ng)
    xdot = zeros(model.nx)

    function objective!(g, xu)
        x = xu[1:model.nx]
        u = xu[model.nx+1:end]
        xdot .= dynamics(model, x, u)
        x_xdot = vcat(x,xdot)
        g[1:Int(ng/2)] .= x_xdot[i_x_xdot_trimmed] .- x_xdot_trimmed
        g[1+Int(ng/2):end] .= x_xdot_trimmed .- x_xdot[i_x_xdot_trimmed]

        return objective(x, xdot, u)
    end

    # println("Sherlock!\n\txu_trim_guess = $xu_trim_guess\n\tng = $ng\n\ttlxu = $lxu\n\tuxu = $uxu\n\tlg = $lg\n\tug = $ug\n\toptions = $options")
    # println("\nTest ===:")
    # for i=1:length(xu_trim_guess)
    #     println(xu_trim_guess[i] .=== xu_trim_guess)
    # end
    xu_opt, f_opt, info = SNOW.minimize(objective!, xu_trim_guess, ng, lxu, uxu, lg, ug, options)

    return xu_opt, f_opt, info
end

max_stars_repo_name: byuflowlab/MultifidelityTrajectories

max_issues_repo_name: byuflowlab/MultifidelityTrajectories

max_forks_repo_name: byuflowlab/MultifidelityTrajectories


text:
@testset "kerneltensorproduct" begin
    rng = MersenneTwister(123456)
    u1 = rand(rng, 10)
    u2 = rand(rng, 10)
    v1 = rand(rng, 5)
    v2 = rand(rng, 5)

    # kernels
    k1 = SqExponentialKernel()
    k2 = ExponentialKernel()
    kernel1 = KernelTensorProduct(k1, k2)
    kernel2 = KernelTensorProduct([k1, k2])

    @test kernel1 == kernel2
    @test kernel1.kernels === (k1, k2) === KernelTensorProduct((k1, k2)).kernels
    @test length(kernel1) == length(kernel2) == 2
    @test string(kernel1) == (
        "Tensor product of 2 kernels:\n\tSquared Exponential Kernel\n\tExponential Kernel"
    )
    @test_throws DimensionMismatch kernel1(rand(3), rand(3))

    @testset "val" begin
        for (x, y) in (((v1, u1), (v2, u2)), ([v1, u1], [v2, u2]))
            val = k1(x[1], y[1]) * k2(x[2], y[2])

            @test kernel1(x, y) == kernel2(x, y) == val
        end
    end

    # Deprecations
    @test TensorProduct === KernelTensorProduct
    @test TensorProduct(k1, k2) == k1 ⊗ k2
    @test TensorProduct((k1, k2)) == k1 ⊗ k2
    @test TensorProduct([k1, k2]) == k1 ⊗ k2

    # Standardised tests.
    TestUtils.test_interface(kernel1, ColVecs{Float64})
    TestUtils.test_interface(kernel1, RowVecs{Float64})
    test_ADs(
        x -> KernelTensorProduct(SqExponentialKernel(), LinearKernel(; c=exp(x[1]))),
        rand(1);
        dims=[2, 2],
    )
    test_params(KernelTensorProduct(k1, k2), (k1, k2))

    @testset "single kernel" begin
        kernel = KernelTensorProduct(k1)
        @test length(kernel) == 1

        @testset "eval" begin
            for (x, y) in (((v1,), (v2,)), ([v1], [v2]))
                val = k1(x[1], y[1])

                @test kernel(x, y) == val
            end
        end
    end
end


text:
using Random
using Makie
using BasicBSpline
using ExportNURBS
Random.seed!(2)

## 1-dim B-spline manifold
p = 2 # degree of polynomial
k = Knots(1:12) # knot vector
P = FastBSplineSpace(p, k) # B-spline space
a = [[i-5, 3*sin(i^2)] for i in 1:dim(P)] # control points
M = BSplineCurve([P], a) # Define B-spline manifold
save_png("docs/src/img/1dim.png", M, unitlength = 50)

## 2-dim B-spline manifold
p = 2 # degree of polynomial
k = Knots(1:8) # knot vector
P = BSplineSpace(p, k) # B-spline space
rand_a = [rand(2) for i in 1:dim(P), j in 1:dim(P)]
a = [[2 * i - 6.5, 2 * j - 6.5] for i in 1:dim(P), j in 1:dim(P)] + rand_a # random generated control points
M = BSplineManifold([P, P], a) # Define B-spline manifold
save_png("docs/src/img/2dim.png", M, unitlength = 50)

## h-Refinement
k₊ = [Knots(3.3, 4.2), Knots(3.8, 3.2, 5.3)]
M′ = refinement(M, k₊ = k₊)
save_png("docs/src/img/2dim_h-refinement.png", M′, unitlength = 50)

## p-Refinement
p₊ = [1,2]
M′ = refinement(M, p₊ = p₊)
save_png("docs/src/img/2dim_p-refinement.png", M′, unitlength = 50)

## Makie
points = [M([u,v]) for u in range(3.0,6.0,length=50), v in range(3.0,6.0,length=50)]
X = [point[1] for point in points]
Y = [point[2] for point in points]
scene = Scene(resolution=(1000,1000))
Makie.surface!(X,Y)
Makie.xlims!(-5,5)
Makie.ylims!(-5,5)
save("docs/src/img/2dim_makie.png", scene)

## Fitting
p1 = 2
p2 = 2
k1 = Knots(-10:10) + p1 * Knots(-10, 10)
k2 = Knots(-10:10) + p2 * Knots(-10, 10)
P1 = FastBSplineSpace(p1, k1)
P2 = FastBSplineSpace(p2, k2)

f(u1,u2) = [2u1 + sin(u1) + cos(u2) + u2 / 2, 3u2 + sin(u2) + sin(u1) / 2 + u1^2 / 6] / 5

a = fittingcontrolpoints(f, P1, P2)
M = BSplineManifold([P1, P2], a)
save_png("docs/src/img/fitting.png", M, unitlength = 50, up = 10, down = -10, left = -10, right = 10)

## Coarse fitting
p1 = 2
p2 = 2
k1 = Knots(-10:5:10) + p1 * Knots(-10, 10)
k2 = Knots(-10:5:10) + p2 * Knots(-10, 10)
P1 = FastBSplineSpace(p1, k1)
P2 = FastBSplineSpace(p2, k2)

f(u) = [2u[1] + sin(u[1]) + cos(u[2]) + u[2] / 2, 3u[2] + sin(u[2]) + sin(u[1]) / 2 + u[1]^2 / 6] / 5

a = fittingcontrolpoints(f, [P1, P2])
M = BSplineManifold([P1, P2], a)
save_png("docs/src/img/fitting_coarse.png", M, unitlength = 50, up = 10, down = -10, left = -10, right = 10)

## Sine curve
p = 3
k = Knots(range(-2π, 2π, length = 8)) + p * Knots(-2π, 2π)
P = FastBSplineSpace(p, k)

f(u) = [u, sin(u)]
a = fittingcontrolpoints(f, P)
M = BSplineManifold([P], a)
save_svg("docs/src/img/sine_curve.svg", M, unitlength = 50, up = 2, down = -2, left = -8, right = 8)


text:
# https://github.com/cvxgrp/scs
"""
   unvec_symm(x, dim)
Returns a dim-by-dim symmetric matrix corresponding to `x`.
`x` is a vector of length dim*(dim + 1)/2, corresponding to a symmetric matrix
```
X = [ X11     X12/√2 ... X1k/√2
      X21/√2  X22    ... X2k/√2
      ...
      Xk1/√2  Xk2/√2 ... Xkk ],
```
where
`vec(X) = (X11, X12, X22, X13, X23, ..., Xkk)`

Note that the factor √2 preserves inner products:
`x'*c = Tr(unvec_symm(c, dim) * unvec_symm(x, dim))`
"""
function unvec_symm_scs(x, dim)
    X = zeros(eltype(x), dim, dim)
    idx = 1
    for i in 1:dim
        for j in 1:i
            if i == j
                X[i,j] = x[idx]
            else
                X[j,i] = X[i,j] = x[idx] / sqrt(2)
            end
            idx += 1
        end
    end
    return X
end

"""
   vec_symm(X)
Returns a vectorized representation of a symmetric matrix `X`.
`vec(X) = (X11, √2*X12, X22, √2*X13, X23, ..., Xkk)`

Note that the factor √2 preserves inner products:
`x'*c = Tr(unvec_symm(c, dim) * unvec_symm(x, dim))`
"""
function vec_symm_scs(X)
    x_vec = sqrt(2).*X[LinearAlgebra.tril(trues(size(X)))']
    idx = 1
    for i in 1:size(X)[1]
        x_vec[idx] =  x_vec[idx]/sqrt(2)
        idx += i + 1
    end
    return x_vec
end

# Write and read CP from file
function cp_to_file(file, params_dict; dense=true)
    A, b, c = params_dict[:A], params_dict[:b], params_dict[:c]
    if isnothing(params_dict[:x_star])
         _cp_to_file(file, (A, b, c), dense=dense)
    else
        x_star = params_dict[:x_star]
        y_star = params_dict[:y_star]
        s_star = params_dict[:s_star]
        _cp_to_file(file, (A, b, c), opt_vals=(x_star, y_star, s_star), dense=dense)
    end
end

function _cp_to_file(file, params; opt_vals=(), dense=true)
    if dense
        writedlm(file, (params..., opt_vals...))
    else
        writedlm(file, (findnz(params[1])..., params[2:3]..., opt_vals...))
    end
end

function cp_from_file(file; dense=true)
    file_vec = readlines(file)
    nvars = length(file_vec)

    if dense
        ~(nvars in [3, 6]) && throw(ArgumentError("Invalid file"))
        offset = 0
        A_flat = parse.(Float64, split(file_vec[1], '\t'))
        b = parse.(Float64, split(file_vec[2+offset], '\t'))
        c = parse.(Float64, split(file_vec[3+offset], '\t'))
        m,n = length(b), length(c)
        A = reshape(A_flat, (m,n))
    else
        ~(nvars in [5, 8]) && throw(ArgumentError("Invalid file"))
        offset = 2
        A_row_inds = parse.(Int, split(file_vec[1], '\t'))
        A_col_inds = parse.(Int, split(file_vec[2], '\t'))
        A_vals = parse.(Float64, split(file_vec[3], '\t'))
        A = sparse(A_row_inds, A_col_inds, A_vals)
        b = parse.(Float64, split(file_vec[2+offset], '\t'))
        c = parse.(Float64, split(file_vec[3+offset], '\t'))
    end

    if nvars in [6, 8]
        x_star = parse.(Float64, split(file_vec[4+offset], '\t'))
        y_star = parse.(Float64, split(file_vec[5+offset], '\t'))
        s_star = parse.(Float64, split(file_vec[6+offset], '\t'))
    else
        x_star, y_star, s_star = nothing, nothing, nothing
    end
    return Dict(:A => A, :b => b, :c => c, :x_star => x_star,
                :y_star => y_star, :s_star => s_star)
end

function derivatives_to_file(file, deriv, adjoint)
    packed_deriv = (deriv[:dA], deriv[:db], deriv[:dc], deriv[:dx], deriv[:dy], deriv[:ds])
    packed_adjoint = (adjoint[:dA], adjoint[:db], adjoint[:dc], adjoint[:dx], adjoint[:dy], adjoint[:ds])
    writedlm(file, (packed_deriv..., packed_adjoint...))
end

function derivatives_from_file(file)
    file_vec = readlines(file)

    function from_file(file_vec, offset)
        ret = Dict()
        ret[:db] = parse.(Float64, split(file_vec[2+offset], '\t'))
        ret[:dc] = parse.(Float64, split(file_vec[3+offset], '\t'))
        m, n = length(ret[:db]), length(ret[:dc])
        ret[:dA] = reshape(parse.(Float64, split(file_vec[1+offset], '\t')), (m,n))

        ret[:dx] = parse.(Float64, split(file_vec[4+offset], '\t'))
        ret[:dy] = parse.(Float64, split(file_vec[5+offset], '\t'))
        ret[:ds] = parse.(Float64, split(file_vec[6+offset], '\t'))
        return ret
    end

    # derivative, adjoint
    return from_file(file_vec, 0), from_file(file_vec, 6)
end

function _vectorized_index_row_major(i, j)
    if i <= j
        k = div((j-1)*j, 2) + i
    else
        k = div((i-1)*i, 2) + j
    end
    return k
end

function index_map_to_row_major(dim)
    return sortperm(index_map_to_col_major(dim))
end

function index_map_to_col_major(dim)
    index_map = []
    for i=1:dim
        for j=i:dim
            k = _vectorized_index_row_major(i, j)
            push!(index_map, k)
        end
    end
    index_map = sortperm(index_map)
    return index_map
end

function rewrite_sdps_to_row_major!(A, b, cones)
    i = 1
    for cone in cones
        if typeof(cone) <: MOI.PositiveSemidefiniteConeTriangle
            d = MOI.dimension(cone)
            s = cone.side_dimension
            index_map = index_map_to_row_major(cone.side_dimension)
            A[i:(i+d-1), :] .= A[i:(i+d-1), :][index_map, :]
            b[i:(i+d-1)] .= b[i:(i+d-1)][index_map]
        end
        i += MOI.dimension(cone)
    end
end

function rewrite_sdps_to_row_major(A, b, cones)
    i = 1
    A_ = copy(A)
    b_ = copy(b)
    for cone in cones
        if typeof(cone) <: MOI.PositiveSemidefiniteConeTriangle
            d = MOI.dimension(cone)
            s = cone.side_dimension
            index_map = index_map_to_row_major(cone.side_dimension)
            A_[i:(i+d-1), :] .= A[i:(i+d-1), :][index_map, :]
            b_[i:(i+d-1)] .= b[i:(i+d-1)][index_map]
        end
        i += MOI.dimension(cone)
    end
    return A_, b_
end

function rewrite_sdps_to_col_major!(A, b, cones)
    i = 1
    for cone in cones
        if typeof(cone) <: MOI.PositiveSemidefiniteConeTriangle
            d = MOI.dimension(cone)
            s = cone.side_dimension
            index_map = index_map_to_col_major(cone.side_dimension)
            A[i:(i+d-1), :] .= A[i:(i+d-1), :][index_map, :]
            b[i:(i+d-1)] .= b[i:(i+d-1)][index_map]
        end
        i += MOI.dimension(cone)
    end
end

function rewrite_sdps_to_col_major(A, b, cones)
    i = 1
    A_ = copy(A)
    b_ = copy(b)
    for cone in cones
        if typeof(cone) <: MOI.PositiveSemidefiniteConeTriangle
            d = MOI.dimension(cone)
            s = cone.side_dimension
            index_map = index_map_to_col_major(cone.side_dimension)
            A_ .= A[i:(i+d-1), :][index_map, :]
            b_ .= b[i:(i+d-1)][index_map]
        end
        i += MOI.dimension(cone)
    end
    return A_, b_
end


text:
#
# Other callbacks
#    ScfPlotTrace()        see src/plotting.jl (imported with Plots.jl)
#    ScfSaveCheckpoints()  see src/jld2io.jl   (imported with JLD2.jl)
#

"""
Default callback function for `self_consistent_field`, which prints a convergence table
"""
function ScfDefaultCallback()
    prev_energies = nothing
    function callback(info)
        !mpi_master() && return info  # Printing only on master
        if info.stage == :finalize
            info.converged || @warn "SCF not converged."
            return info
        end
        collinear = info.basis.model.spin_polarization == :collinear
        dVol      = info.basis.model.unit_cell_volume / prod(info.basis.fft_size)

        if info.n_iter == 1
            E_label = haskey(info.energies, "Entropy") ? "Free energy" : "Energy"
            magn    = collinear ? ("   Magnet", "   ------") : ("", "")
            @printf "n     %-12s      Eₙ-Eₙ₋₁     ρout-ρin%s   Diag\n" E_label magn[1]
            @printf "---   ---------------   ---------   --------%s   ----\n" magn[2]
        end
        E    = isnothing(info.energies) ? Inf : info.energies.total
        Δρ   = norm(info.ρout.fourier - info.ρin.fourier)
        magn = isnothing(info.ρ_spin_out) ? NaN : sum(info.ρ_spin_out.real) * dVol

        Estr   = (@sprintf "%+15.12f" round(E, sigdigits=13))[1:15]
        prev_E = prev_energies === nothing ? Inf : prev_energies.total
        ΔE     = prev_E == Inf ? "      NaN" : @sprintf "% 3.2e" E - prev_E
        Mstr = collinear ? "   $((@sprintf "%6.3f" round(magn, sigdigits=4))[1:6])" : ""
        diagiter = sum(info.diagonalization.iterations) / length(info.diagonalization.iterations)
        @printf "% 3d   %s   %s   %2.2e%s   % 3.1f \n" info.n_iter Estr ΔE Δρ Mstr diagiter
        prev_energies = info.energies

        flush(stdout)
        info
    end
end

"""
Flag convergence as soon as total energy change drops below tolerance
"""
function ScfConvergenceEnergy(tolerance)
    energy_total = NaN

    function is_converged(info)
        info.energies === nothing && return false # first iteration

        # The ρ change should also be small, otherwise we converge if the SCF is just stuck
        norm(info.ρout.fourier - info.ρin.fourier) > 10sqrt(tolerance) && return false

        etot_old = energy_total
        energy_total = info.energies.total
        abs(energy_total - etot_old) < tolerance
    end
    return is_converged
end

"""
Flag convergence by using the L2Norm of the change between
input density and unpreconditioned output density (ρout)
"""
function ScfConvergenceDensity(tolerance)
    info -> norm(info.ρout.fourier - info.ρin.fourier) < tolerance
end

"""
Determine the tolerance used for the next diagonalization. This function takes
``|ρnext - ρin|`` and multiplies it with `ratio_ρdiff` to get the next `diagtol`,
ensuring additionally that the returned value is between `diagtol_min` and `diagtol_max`
and never increases.
"""
function ScfDiagtol(;ratio_ρdiff=0.2, diagtol_min=nothing, diagtol_max=0.03)
    function determine_diagtol(info)
        isnothing(diagtol_min) && (diagtol_min = 100eps(real(eltype(info.ρin))))
        info.n_iter ≤ 1 && return diagtol_max
        info.n_iter == 2 && (diagtol_max /= 5)  # Enforce more accurate Bloch wave

        diagtol = norm(info.ρnext.fourier - info.ρin.fourier) * ratio_ρdiff
        # TODO Quantum espresso divides diagtol by the number of electrons
        diagtol = min(diagtol_max, diagtol)  # Don't overshoot
        diagtol = max(diagtol_min, diagtol)  # Don't undershoot
        @assert isfinite(diagtol)

        # Adjust maximum to ensure diagtol may only shrink during an SCF
        diagtol_max = min(diagtol, diagtol_max)
        diagtol
    end
end


text:
const modelName = "portfolio"

const modelFile = modelName * "_sigmodel_template.jl"

# directory name for storing results
const workingDir = "C:/Users/rkannan/Desktop/SA for CCP/Sigmoidal approximation/"
const baseDirName = workingDir * "experiments/" * modelName * "_sigapprox/"
mkpath(baseDirName)

#*========= INCLUDE MODEL FILES ================
include(modelFile) # parameters of the model
#*==============================================

const details_file = baseDirName * modelName * ".txt"
initializeFiles(details_file)

# loop over the sample sizes
for iter = 1:size(NumSamplesFactor,1)

	riskLevel::Float64 = riskLevels[1]

	numScenarios::Int64 = round(NumSamplesFactor[iter]/riskLevel)

	println("")
	println("***************************************************")
	println("Iteration #: ",iter,"/",size(NumSamplesFactor,1),"  with riskLevels = ",riskLevel," and numScenarios = ",numScenarios)
	println("***************************************************")
	println("")

	# create folder for storing results
	dirName::String = baseDirName * string(iter) * "/"
	mkpath(dirName)

	# write num samples and risk level to text file
	risk_file::String = dirName * "specifiedRiskLevel.txt"
	open(risk_file, "w") do f
		write(f,"$riskLevel \n")
	end

	scenario_file::String = dirName * "numScenarios.txt"
	open(scenario_file, "w") do f
		write(f,"$numScenarios \n")
	end

	#####
	# run as many replicates as needed
	for rep = 1:NumReplicates

		tic()

		println("REPLICATE #",rep)

		solveSigmoidalApproximationModel(riskLevel,numScenarios,dirName)

		gc()

		println("")

		elapsed_time::Float64 = toq()

		sleep(min(30,1.0+0.1*elapsed_time))

	end # for replicates
	#####

	sleep(5.0)

end # for risk levels

max_stars_repo_name: rohitkannan/SA-for-CCP

max_issues_repo_name: rohitkannan/SA-for-CCP

max_forks_repo_name: rohitkannan/SA-for-CCP


text:
using Statistics
using DelimitedFiles

include("MatrixGenerator.jl")
include("Calculations.jl")
include("MeanFieldSolver.jl")
include("Test.jl")

# get inputs, initialization
lines = readlines("InputBoseHubbard")
# inverse temperature
temp_name, β_str = split(lines[1])
const β = parse(Float64, β_str)
# number of sites
temp_name, num_sites_str = split(lines[2])
const num_sites = parse(Int64, num_sites_str)
# number of bosons
temp_name, num_bosons_str = split(lines[3])
const num_bosons = parse(Int64, num_bosons_str)
# length of the imaginary time interval
temp_name, Δτ_str = split(lines[4])
const Δτ = parse(Float64, Δτ_str)
# hopping constant
temp_name, t_str = split(lines[5])
const t = parse(Float64, t_str)
# on-site repulsion constant
temp_name, U_str = split(lines[6])
const U = parse(Float64, U_str)
# number of iterations
temp_name, num_iterations_str = split(lines[7])
const num_iterations = convert(Int64, parse(Float64, num_iterations_str))
# flag for the background subtraction
temp_name, flag_background_subtraction_str = split(lines[8])
const flag_background_subtraction = parse(Int16, flag_background_subtraction_str)
# flag for calculating the density matrix
temp_name, flag_density_matrix_str = split(lines[9])
const flag_density_matrix = convert(Int64, parse(Int16, flag_density_matrix_str))

const background_dist = ones(num_sites) / num_sites #  Equally distributed among sites
const residue_term = sum(β * (U/2) * background_dist.^2)

# Constants for calculations
# Kinetic matrix in each imaginary time slice is fixed throughout the calculation.
const kinetic_matrix_imgtime = complex(exp(-kinetic_matrix(num_sites, t) * Δτ / 2))
# number of imaginary time steps
const L = convert(Int64, round(β / Δτ))
# coefficient of the H-S transformed potential term
const AF_coefficient = sqrt(complex(-Δτ * U))

#mean_value, error_bar = Z_iteration(num_iteration, L)

#write_to_file(mean_value, error_bar, Δτ)

#function write_to_file(mean_value, error_bar, Δτ)
    #output_file = open("result","a+")
    #write(output_file, string(β), " ", string(Δτ), " ", string(mean_value), " ", string(error_bar), "\n")
    #close(output_file)
#end

Ek_mean, Ek_error, Ep_mean, Ep_error, E_tot_mean, E_tot_error_bar = measure_TotEnergy(num_iterations)
println("U = ", U)
println("Beta = ", β)
println("Ns = ", num_sites)
println("NBosons = ", num_bosons)
println("Ek = ", Ek_mean, ", Error = ", Ek_error)
println("Ep = ", Ep_mean, ", Error = ", Ep_error)
println("Etot = ", E_tot_mean, ", Error = ", E_tot_error_bar)

max_stars_repo_name: TongSericus/CanonicalAFQMC

max_issues_repo_name: TongSericus/CanonicalAFQMC

max_forks_repo_name: TongSericus/CanonicalAFQMC


text:
using Test
using Clustering

@testset "fuzzy_cmeans()" begin

@testset "Argument checks" begin
    Random.seed!(34568)
    @test_throws ArgumentError fuzzy_cmeans(randn(2, 3), 1, 2.0)
    @test_throws ArgumentError fuzzy_cmeans(randn(2, 3), 4, 2.0)
    @test_throws ArgumentError fuzzy_cmeans(randn(2, 3), 2, 1.0)
    for disp in keys(Clustering.DisplayLevels)
        @test fuzzy_cmeans(randn(2, 3), 2, 2.0, tol=0.1, display=disp) isa FuzzyCMeansResult
    end
end

Random.seed!(34568)

d = 3
n = 1000
k = 5

x = rand(d, n)

@testset "fuzziness = 2.0" begin
    fuzziness = 2.0
    Random.seed!(34568)
    r = fuzzy_cmeans(x, k, fuzziness)
    @test isa(r, FuzzyCMeansResult{Float64})
    @test nclusters(r) == k
    @test size(r.centers) == (d, k)
    @test size(r.weights) == (n, k)
    @test all(0 .<= r.weights .<= 1)
    @test sum(r.weights, dims=2) ≈ fill(1.0, n)

    @test wcounts(r) isa Vector{Float64}
    @test length(wcounts(r)) == n
    @test all(0 .<= wcounts(r) .<= n)
    @test sum(wcounts(r)) ≈ n
end

@testset "fuzziness = 3.0" begin
    fuzziness = 3.0
    Random.seed!(34568)
    r = fuzzy_cmeans(x, k, fuzziness)
    @test isa(r, FuzzyCMeansResult{Float64})
    @test nclusters(r) == k
    @test size(r.centers) == (d, k)
    @test size(r.weights) == (n, k)
    @test sum(r.weights, dims=2) ≈ fill(1.0, n)
    @test all(0 .<= r.weights .<= 1)

    @test wcounts(r) isa Vector{Float64}
    @test length(wcounts(r)) == n
    @test all(0 .<= wcounts(r) .<= n)
    @test sum(wcounts(r)) ≈ n
end

@testset "Abstract data matrix" begin
    fuzziness = 2.0
    Random.seed!(34568)
    r = fuzzy_cmeans(view(x, :, :), k, fuzziness)
    @test isa(r, FuzzyCMeansResult{Float64})
    @test nclusters(r) == k
    @test size(r.centers) == (d, k)
    @test size(r.weights) == (n, k)
    @test sum(r.weights, dims=2) ≈ fill(1.0, n)
    @test all(0 .<= r.weights .<= 1)

    @test wcounts(r) isa Vector{Float64}
    @test length(wcounts(r)) == n
    @test all(0 .<= wcounts(r) .<= n)
    @test sum(wcounts(r)) ≈ n
end

end


text:
using NearestNeighbors
import Parameters.@with_kw_noshow
export KDTree
export temporalprediction
export PredictionParameters

@with_kw_noshow struct PredictionParameters{Φ,BC,X,
        LM <: AbstractLocalModel, NT<:AbstractNeighborhood}
    em::AbstractSpatialEmbedding{Φ,BC,X}
    method::LM = AverageLocalModel(ω_safe)
    ntype::NT = FixedMassNeighborhood(3)
    treetype = KDTree
end

###########################################################################################
#                        Iterated Time Series Prediction                                  #
###########################################################################################

"""
    temporalprediction(U, em::AbstractSpatialEmbedding, tsteps; kwargs...)
Perform a spatio-temporal time series prediction for `tsteps` iterations,
using local weighted modeling [1] give a time series of the form
`U::AbstractVector{<:AbstractArray{T, Φ}}`.

The returned data always contains the final state of `U` as starting point
(total returned length is `tsteps+1`).
The reconstruction process is defined by `em`.
For available methods and interfaces see [`AbstractSpatialEmbedding`](@ref).

## Keyword Arguments
  * `ttype = KDTree` : Type/Constructor of tree structure. So far only tested with `KDTree`.
  * `method = AverageLocalModel(ω_safe)` : Subtype of [`AbstractLocalModel`](@ref).
  * `ntype = FixedMassNeighborhood(3)` : Subtype of `AbstractNeighborhood`.
  * `initial_ts = U` : Initial states for prediction (same type as `U`).
    Must have at least as many states as the maximum delay time used. Defaults to the training set `U`.
  * `progress = true` : To print progress done.

## Description
This method works similarly to [`localmodel_tsp`](@ref), by expanding the concept
of delay embedding to spatially extended systems. Instead of reconstructing
complete states of the system, local states are used.
See [`AbstractSpatialEmbedding`](@ref) for details on the embedding.
Predictions are then performed frame by frame and point py point. Once all values for a new
frame are found, the frame is added to the end of the timeseries and used to generate
new prediction queries for the next time step.

## Performance Notes
Be careful when choosing embedding parameters as memory usage and computation time
depend strongly on the resulting embedding dimension.

## References
[1] : U. Parlitz & C. Merkwirth, [Phys. Rev. Lett. **84**, pp 1890 (2000)](https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.84.1890)
"""
function temporalprediction(s,
    em::AbstractSpatialEmbedding,
    tsteps;
    method = AverageLocalModel(ω_safe),
    ntype = FixedMassNeighborhood(3),
    ttype=KDTree,
    kwargs...)

    params = PredictionParameters(em, method, ntype, ttype)
    return temporalprediction(params, s, tsteps; kwargs...)
end

function temporalprediction(params, s, tsteps; progress=true, kwargs...)
    progress && println("Reconstructing")
    R = reconstruct(s,params.em)

    #Prepare tree but remove the last reconstructed states first
    progress && println("Creating Tree")
    L = length(R)
    M = get_num_pt(params.em)

    tree = params.treetype(R[1:L-M])
    temporalprediction(params, s, tsteps, R, tree; progress=progress, kwargs...)
end

function temporalprediction(params, s,tsteps, R, tree;
        initial_ts=s, #optional start for prediction
        progress=true) where {T, Φ, BC, X}
    em = params.em
    @assert outdim(em) == size(R,2)
    @assert length(initial_ts) > get_τmax(em)
    num_pt = get_num_pt(em)
    #New state that will be predicted, allocate once and reuse
    state = similar(s[1])

    #Prepare starting point of prediction timeseries
    spred = working_ts(initial_ts, params.em)

    for n=1:tsteps
        progress && println("Working on Frame $(n)/$(tsteps)")
        queries = gen_queries(spred, em)

        #Iterate over queries/ spatial points
        for m=1:num_pt
            q = queries[m]

            #Find neighbors
            idxs,dists = neighborhood_and_distances(q,R,tree,params.ntype)

            xnn = R[idxs]
            #Retrieve ynn
            ynn = map(idxs) do idx
                #Indices idxs are indices of R. Convert to indices of s
                t,α = convert_idx(idx,em)
                s[t+1][α]
            end
            state[m] = params.method(q,xnn,ynn,dists)[1]
        end
        push!(spred,copy(state))
    end

    cut_off_beginning!(spred,tsteps)
    return spred
end

function working_ts(s,em)
    L = length(s)
    τmax = get_τmax(em)
    return s[L-τmax : L]
end

function gen_queries(s,em)
    L = length(s)
    τmax = get_τmax(em)
    s_slice = view( s, L-τmax:L)
    return reconstruct(s_slice, em)
end

function convert_idx(idx, em)
    τmax = get_τmax(em)
    num_pt = get_num_pt(em)
    t = 1 + (idx-1) ÷ num_pt + get_τmax(em)
    α = 1 + (idx-1) % num_pt
    return t,α
end

cut_off_beginning!(s,tsteps) = (N=length(s); deleteat!(s, 1:N-tsteps-1))


text:
#=

      Temporal convolution

=#

export TemporalConvolution

type TemporalConvolution{D<:Device, F<:Float} <: Layer
  kW::UInt

  dW::UInt

  input_framesize::UInt
  output_framesize::UInt

  weight::matrixParameters(D,F)

  bias::vectorParameters(D,F)

  # State
  output::denseRealMatrix(D,F)
  grad_input::denseRealMatrix(D,F)

  TemporalConvolution(device, kw, dW, input_framesize, output_framesize, weight, bias) =
    new(kw, dW, input_framesize, output_framesize, weight, bias, array(device, F, 0, 0), array(device, F, 0, 0))
end

# @doc doc"Creates a new temporal convolution operator

# input_framesize The size of the inputs vectors
# output_framesize The size of the output vectors
# kW The kernel width - how many inputs are taken into account
# dW The move width - the step size
# " ->
function TemporalConvolution{D<:Device, F<:Float}(device::D, ::Type{F}, input_framesize::UInt, output_framesize::UInt, kW::UInt, dW::Int=1)
  weight = matrixParameters(device, F, output_framesize, input_framesize * kW)
  bias = vectorParameters(device, F, output_framesize)

  self = TemporalConvolution{D,F}(device, kW, dW, input_framesize, output_framesize, weight, bias)

  return self
end

function init!{D<:Device, F<:Float}(s::TemporalConvolution{D,F}, stdv=Nullable{F}())
   if isnull(stdv)
      stdv = 1/ sqrt(s.kW * s.input_framesize)
   else
      stdv = stdv * math.sqrt(3)
   end

    uniform!(s.weight.values, -stdv, stdv)
    uniform!(s.bias.values, -stdv, stdv)
end

# Case where the input is (input_size, nb_samples)
function forward!{D<:Device,F<:Float}(m::TemporalConvolution{D,F}, input::DenseArray{F, 2})
  # Prepare
  nInputFrame = size(input, 2)
  nOutputFrame = div(nInputFrame - m.kW, m.dW) + 1

  output = @ensuresize m.output, nOutputFrame, m.output_framesize

  weight = m.weight.values

  # Compute the convolution
  pos::Int = 1 # Position in the input

  for k = 1:nOutputFrame
    outputview = unsafe_view(output, :, k)
    copy!(outputview, m.bias.values)
    inputview = flatten_view(view(input, :, pos:(pos + m.kW - 1)))

    BLAS.gemv!('N', one(F), weight, inputview, one(F), outputview)

    # Advance the window
    pos += m.dW
  end

  output
end

function compute_inputgradient!{D<:Device, F<:Float}(m::TemporalConvolution{D,F}, input::DenseRealMatrix, gradOutput::DenseRealMatrix)
 # Prepare
  nInputFrame = size(input, 2)
  nOutputFrame = div(nInputFrame - m.kW, m.dW) + 1

  grad_input = @ensuresize m.grad_input, (size(input, 1), size(input, 2))

  weight = m.weight.values

  # Compute the convolution
  pos::Int = 1 # Position in the input

  for k = 1:nOutputFrame
    outputview = unsafe_view(output, :, k)
    copy!(outputview, m.bias.values)
    inputview = flatten_view(view(input, :, pos:(pos + m.kW - 1)))

    BLAS.gemv!('N', one(F), weight, inputview, one(F), grad_input)

    # Advance the window
    pos += m.dW
  end
end

function update_gradient!{D<:Device, F<:Float}(m::TemporalConvolution{D,F}, input::DenseRealMatrix, gradOutput::DenseRealMatrix, scale::F=1.)
  @assert false
end


text:
using HttpServer
using WebSockets
using JSON

decodeMessage(msg) = bytestring(msg)

global jsonwebsocketconnections = Dict{Int,WebSocket}()

function send(msg)
  global jsonwebsocketconnections
  jsonmsg = JSON.json(msg)
  if length(jsonwebsocketconnections) == 0
    println("No one to send to! ", jsonmsg)
  end
  for (id, client) in jsonwebsocketconnections
    println("Sending to ", id, ": ", jsonmsg)
    try
      write(client, jsonmsg)
    catch error
      println(error)
    end
  end
end

function start_websocket_handler_while_processing(processingFunc; mainpage = "index.html", port = 8084)

  onepage = readall(mainpage)
  httph = HttpHandler() do req::Request, res::Response
    Response(onepage)
  end

  @async begin
    while true
      processingFunc()
    end
  end

  wsh = WebSocketHandler() do req, client
    global jsonwebsocketconnections
    println("Got request from client: ", client)
    jsonwebsocketconnections[client.id] = client
    while true
      msg = decodeMessage(read(client))
      println("Received from ", client.id, ": ", msg)
      try
        jsonmsg = JSON.parse(msg)
        handle_json_message(req, client, jsonmsg)
      catch error
        error("Could not parse and/or handle the msg: ", msg)
      end
    end
  end
  println("Now starting server!")
  server = Server(httph, wsh)
  run(server, 8084)
end

global currentrow = 1

session_timestamp() = begin
  t = time()
  strftime("%Y-%m-%dT%H:%M:%S.", t) * string(int(round(t % 1, 3) * 1000)) * "Z"
end

function rand_test_outcome()
  pr = rand()
  if pr < 0.10
    -1 # Failure
  elseif pr < 0.95
    1 # Pass
  else
    0 # Missing info
  end
end

# For now the processing is just to send some random data every now and then...
function processing_func()
  numtestcases = 10
  session = session_timestamp()
  @sync @parallel for tc in 1:numtestcases
    sleep(1 + 5 * rand()) # Fake that processing takes time...
    msgtype = (rand() < 0.80) ? "dataupdate" : "unspecified"
    send({"type" => msgtype, "name" => "tom",
      "session" => session, "testcase" => tc,
      "value" => rand_test_outcome()})
  end
end

start_websocket_handler_while_processing(processing_func)


text:
####################
#
# Point
#
####################

# Union of Shapefile Point types
ShpPoint = Union{Shapefile.Point, Shapefile.PointZ, Shapefile.PointM}

#=
For an array of ESRIShape points ???
1. convert to a different coordinate system (optional)
2. convert to a Compose circle

Note:
Units in Compose can be either:
- Context - cx, cy - using the axis coordinates
- Relative - w (width), h (height) - using the size of the image
- Absolute - pt, mm, cm, inch
If you just provide numbers e.g. if x_pts below was an array of Float64
the units would default to context units.
=#
function ESRItoComposeCircle(points::AbstractArray{T, 1}, convertcoords,
                             radius) where {T<:ShpPoint}
    x_pts = Array{Measures.Length{:cx, Float64}}(undef, 0)
    y_pts = Array{Measures.Length{:cy, Float64}}(undef, 0)
    for pt in points
        x, y = convertcoords(pt.x, pt.y)
        push!(x_pts, x * cx)
        push!(y_pts, y * cy)
    end
    Compose.circle(x_pts, y_pts, radius)
end

# Plot an array of points
function draw_shp(shapes::AbstractArray{T, 1}, canvas, convertcoords,
                  line_width, line_color, fill_color, radius) where {T<:ShpPoint}
    compose(canvas, (context(), ESRItoComposeCircle(shapes, convertcoords, radius),
                     linewidth(line_width), stroke(line_color), fill(fill_color)))
end


text:
"""
    job_query_by_id(id::Int64)

Search job by `job.id` in the queue.

Return `job::Job` if found, `nothing` if not found.
"""
function job_query_by_id(id::Int64)
    res = nothing
    wait_for_lock()
    try
        for job in JOB_QUEUE
            if job.id == id
                res = job
                break
            end
        end
        if isnothing(res)
            for job in JOB_QUEUE_OK
                if job.id == id
                    res = job
                    break
                end
            end
        end
    catch e
        rethrow(e)
    finally
        release_lock()
    end
    return res
end # function

job_query = job_query_by_id

function job_query_by_id_no_lock(id::Int64)
    for job in JOB_QUEUE
        if job.id == id
            return job
        end
    end
    for job in JOB_QUEUE_OK
        if job.id == id
            return job
        end
    end
    return nothing # not found
end # function

job_query = job_query_by_id


text:
# ---
# title: 1585. Check If String Is Transformable With Substring Sort Operations
# id: problem1585
# author: Tian Jun
# date: 2020-10-31
# difficulty: Hard
# categories: String, Greedy
# link: <https://leetcode.com/problems/check-if-string-is-transformable-with-substring-sort-operations/description/>
# hidden: true
# ---
#
# Given two strings `s` and `t`, you want to transform string `s` into string
# `t` using the following operation any number of times:
#
#   * Choose a **non-empty** substring in `s` and sort it in-place so the characters are in  **ascending order**.
#
# For example, applying the operation on the underlined substring in `"1 _4234_
# "` results in `"1 _2344_ "`.
#
# Return `true` if _it is possible to transform string`s` into string `t`_.
# Otherwise, return `false`.
#
# A **substring**  is a contiguous sequence of characters within a string.
#
#
#
# **Example 1:**
#
#
#
#     Input: s = "84532", t = "34852"
#     Output: true
#     Explanation: You can transform s into t using the following sort operations:
#     "84 _53_ 2" (from index 2 to 3) -> "84 _35_ 2"
#     " _843_ 52" (from index 0 to 2) -> " _348_ 52"
#
#
# **Example 2:**
#
#
#
#     Input: s = "34521", t = "23415"
#     Output: true
#     Explanation: You can transform s into t using the following sort operations:
#     " _3452_ 1" -> " _2345_ 1"
#     "234 _51_ " -> "234 _15_ "
#
#
# **Example 3:**
#
#
#
#     Input: s = "12345", t = "12435"
#     Output: false
#
#
# **Example 4:**
#
#
#
#     Input: s = "1", t = "2"
#     Output: false
#
#
#
#
# **Constraints:**
#
#   * `s.length == t.length`
#   * `1 <= s.length <= 105`
#   * `s` and `t` only contain digits from `'0'` to `'9'`.
#
#
## @lc code=start
using LeetCode

## add your code here:
## @lc code=end


text:
using Intervals: canonicalize

@testset "AnchoredInterval" begin
    dt = DateTime(2016, 8, 11, 2)

    @testset "constructor" begin
        expected = AnchoredInterval{Hour(-1), DateTime}(dt, Inclusivity(false, true))
        @test AnchoredInterval{Hour(-1), DateTime}(dt) == expected
        @test AnchoredInterval{Hour(-1)}(dt) == expected
        @test HourEnding{DateTime}(dt) == expected
        @test HourEnding(dt) == expected
        @test HE(dt) == expected
        @test HE(dt - Minute(59)) == expected

        expected = AnchoredInterval{Hour(1), DateTime}(dt, Inclusivity(true, false))
        @test AnchoredInterval{Hour(1), DateTime}(dt) == expected
        @test AnchoredInterval{Hour(1)}(dt) == expected
        @test HourBeginning{DateTime}(dt) == expected
        @test HourBeginning(dt) == expected
        @test HB(dt) == expected
        @test HB(dt + Minute(59)) == expected

        # Lazy inclusivity constructor
        @test HourEnding{DateTime}(dt, true, false) ==
            HourEnding{DateTime}(dt, Inclusivity(true, false))
        @test HourEnding(dt, true, false) == HourEnding(dt, Inclusivity(true, false))
        @test AnchoredInterval{Day(1), DateTime}(dt, false, false) ==
            AnchoredInterval{Day(1), DateTime}(dt, Inclusivity(false, false))
        @test AnchoredInterval{Day(1)}(dt, false, false) ==
            AnchoredInterval{Day(1)}(dt, Inclusivity(false, false))
        @test HE(dt, true, true) == HourEnding(dt, Inclusivity(true, true))
        @test HB(dt, true, true) == HourBeginning(dt, Inclusivity(true, true))

        # Non-period AnchoredIntervals
        @test AnchoredInterval{-10}(10) isa AnchoredInterval
        @test AnchoredInterval{25}('a') isa AnchoredInterval
    end

    @testset "non-ordered" begin
        ending = AnchoredInterval{-Inf}(Inf)
        @test isequal(first(ending), NaN)
        @test isequal(last(ending), Inf)
        @test ending == ending
        @test isequal(ending, ending)

        beginning = AnchoredInterval{Inf}(-Inf)
        @test isequal(first(beginning), -Inf)
        @test isequal(last(beginning), NaN)
        @test beginning == beginning
        @test isequal(beginning, beginning)
    end

    @testset "hash" begin
        # Need a complicated enough element type for this test to ever fail
        zdt = now(tz"Europe/London")
        a = HE(zdt)
        b = deepcopy(a)
        @test hash(a) == hash(b)
    end

    @testset "conversion" begin
        he = HourEnding(dt)
        hb = HourBeginning(dt)
        @test convert(DateTime, he) == dt
        @test convert(DateTime, hb) == dt
        @test convert(Interval, he) == Interval(dt - Hour(1), dt, Inclusivity(false, true))
        @test convert(Interval, hb) == Interval(dt, dt + Hour(1), Inclusivity(true, false))
        @test convert(Interval, he) == Interval(dt - Hour(1), dt, Inclusivity(false, true))
        @test convert(Interval, hb) == Interval(dt, dt + Hour(1), Inclusivity(true, false))
        @test convert(Interval{DateTime}, he) == Interval(dt - Hour(1), dt, Inclusivity(false, true))
        @test convert(Interval{DateTime}, hb) == Interval(dt, dt + Hour(1), Inclusivity(true, false))
    end

    @testset "eltype" begin
        @test eltype(AnchoredInterval{-10}(10)) == Int
        @test eltype(AnchoredInterval{25}('a')) == Char
        @test eltype(AnchoredInterval{Day(1)}(today())) == Date
        @test eltype(AnchoredInterval{Day(1),DateTime}(today())) == DateTime
        @test eltype(HourEnding(now())) == DateTime
    end

    @testset "accessors" begin
        inc = Inclusivity(true, true)
        P = Minute(-15)
        interval = AnchoredInterval{P}(dt, inc)

        @test first(interval) == DateTime(2016, 8, 11, 1, 45)
        @test last(interval) == dt
        @test span(interval) == -P
        @test inclusivity(interval) == inc

        inc = Inclusivity(false, false)
        P = Day(1)
        interval = AnchoredInterval{P}(Date(dt), inc)

        @test first(interval) == Date(2016, 8, 11)
        @test last(interval) == Date(2016, 8, 12)
        @test span(interval) == P
        @test inclusivity(interval) == inc

        # DST transition
        endpoint = ZonedDateTime(2018, 3, 11, 3, tz"America/Winnipeg")
        interval = AnchoredInterval{Hour(-2)}(endpoint)
        @test span(interval) == Hour(2)

        startpoint = ZonedDateTime(2018, 3, 11, tz"America/Winnipeg")
        interval = AnchoredInterval{Day(1)}(startpoint)
        @test first(interval) == startpoint
        @test last(interval) == ZonedDateTime(2018, 3, 12, tz"America/Winnipeg")
        @test span(interval) == Day(1)

        endpoint = ZonedDateTime(2018, 11, 4, 2, tz"America/Winnipeg")
        interval = AnchoredInterval{Hour(-2)}(endpoint)
        @test span(interval) == Hour(2)

        startpoint = ZonedDateTime(2018, 11, 4, tz"America/Winnipeg")
        interval = AnchoredInterval{Day(1)}(startpoint)
        @test first(interval) == startpoint
        @test last(interval) == ZonedDateTime(2018, 11, 5, tz"America/Winnipeg")
        @test span(interval) == Day(1)

        # Non-period AnchoredIntervals
        interval = AnchoredInterval{-10}(10)
        @test first(interval) == 0
        @test last(interval) == 10
        @test inclusivity(interval) == Inclusivity(false, true)
        @test span(interval) == 10

        interval = AnchoredInterval{25}('a')
        @test first(interval) == 'a'
        @test last(interval) == 'z'
        @test inclusivity(interval) == Inclusivity(true, false)
        @test span(interval) == 25
    end

    @testset "display" begin
        # Notes on compatibility changes and when they can be updated:
        #
        # When the minimum version of TimeZones requires that `repr` is fixed
        # - `repr(ZonedDateTime(...))` and `$ZonedDateTime` can be changed to be hardcoded
        #
        # When dropping VERSION < v"1.2.0-DEV.29" (https://github.com/JuliaLang/julia/pull/30200)
        # - `repr(Date(...))` and `repr(DateTime(...))` can be converted to hardcode strings
        #
        # When dropping VERSION < v"1.2.0-DEV.223" (https://github.com/JuliaLang/julia/pull/30817)
        # - `repr(Period(...))`can be converted to hardcode strings

        @test sprint(show, AnchoredInterval{Hour(-1)}) ==
            mod_prefix * "AnchoredInterval{$(repr(Hour(-1))),T} where T"
        @test sprint(show, AnchoredInterval{Hour(1)}) ==
            mod_prefix * "AnchoredInterval{$(repr(Hour(1))),T} where T"

        @test sprint(show, AnchoredInterval{Day(-1)}) ==
            mod_prefix * "AnchoredInterval{$(repr(Day(-1))),T} where T"
        @test sprint(show, AnchoredInterval{Day(1)}) ==
            mod_prefix * "AnchoredInterval{$(repr(Day(1))),T} where T"
        @test sprint(show, AnchoredInterval{Day(-1), DateTime}) ==
            mod_prefix * "AnchoredInterval{$(repr(Day(-1))),DateTime}"
        @test sprint(show, AnchoredInterval{Day(1), DateTime}) ==
            mod_prefix * "AnchoredInterval{$(repr(Day(1))),DateTime}"

        # Tuples contain fields: interval, printed, shown
        tests = [
            (
                HourEnding(dt),
                "(2016-08-11 HE02]",
                string(
                    "AnchoredInterval{$(repr(Hour(-1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2))), Inclusivity(false, true))",
                ),
            ),
            (
                HourEnding(DateTime(2013, 2, 13), Inclusivity(true, false)),
                "[2013-02-12 HE24)",
                string(
                    "AnchoredInterval{$(repr(Hour(-1))),DateTime}",
                    "($(repr(DateTime(2013, 2, 13))), Inclusivity(true, false))",
                ),
            ),
            (
                HourEnding(dt + Minute(15) + Second(30)),
                "(2016-08-11 HE02:15:30]",
                string(
                    "AnchoredInterval{$(repr(Hour(-1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2, 15, 30))), Inclusivity(false, true))",
                ),
            ),
            (
                HourEnding(dt + Millisecond(2)),
                "(2016-08-11 HE02:00:00.002]",
                string(
                    "AnchoredInterval{$(repr(Hour(-1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2, 0, 0, 2))), Inclusivity(false, true))",
                ),
            ),
            (
                HourEnding(DateTime(2013, 2, 13, 0, 1), Inclusivity(true, false)),
                "[2013-02-13 HE00:01:00)",
                string(
                    "AnchoredInterval{$(repr(Hour(-1))),DateTime}",
                    "($(repr(DateTime(2013, 2, 13, 0, 1))), Inclusivity(true, false))",
                ),
            ),
            (
                HourBeginning(dt),
                "[2016-08-11 HB02)",
                string(
                    "AnchoredInterval{$(repr(Hour(1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2))), Inclusivity(true, false))",
                ),
            ),
            (
                HourBeginning(DateTime(2013, 2, 13), Inclusivity(false, true)),
                "(2013-02-13 HB00]",
                string(
                    "AnchoredInterval{$(repr(Hour(1))),DateTime}",
                    "($(repr(DateTime(2013, 2, 13))), Inclusivity(false, true))",
                ),
            ),
            (
                HourEnding(ZonedDateTime(dt, tz"America/Winnipeg")),
                "(2016-08-11 HE02-05:00]",
                string(
                    "AnchoredInterval{$(repr(Hour(-1))),$ZonedDateTime}",
                    "($(repr(ZonedDateTime(dt, tz"America/Winnipeg"))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Year(-1)}(Date(dt)),
                "(YE 2016-08-11]",
                string(
                    "AnchoredInterval{$(repr(Year(-1))),Date}",
                    "($(repr(Date(2016, 8, 11))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Year(-1)}(ceil(Date(dt), Year)),
                "(YE 2017-01-01]",
                string(
                    "AnchoredInterval{$(repr(Year(-1))),Date}",
                    "($(repr(Date(2017, 1, 1))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Month(-1)}(dt),
                "(MoE 2016-08-11 02:00:00]",
                string(
                    "AnchoredInterval{$(repr(Month(-1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2, 0, 0))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Month(-1)}(ceil(dt, Month)),
                "(MoE 2016-09-01]",
                string(
                    "AnchoredInterval{$(repr(Month(-1))),DateTime}",
                    "($(repr(DateTime(2016, 9, 1))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Day(-1)}(DateTime(dt)),
                "(DE 2016-08-11 02:00:00]",
                string(
                    "AnchoredInterval{$(repr(Day(-1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Day(-1)}(ceil(DateTime(dt), Day)),
                "(DE 2016-08-12]",
                string(
                    "AnchoredInterval{$(repr(Day(-1))),DateTime}",
                    "($(repr(DateTime(2016, 8, 12))), Inclusivity(false, true))",
                ),
            ),
            (
                # Date(dt) will truncate the DateTime to the nearest day
                AnchoredInterval{Day(-1)}(Date(dt)),
                "(DE 2016-08-11]",
                string(
                    "AnchoredInterval{$(repr(Day(-1))),Date}",
                    "($(repr(Date(2016, 8, 11))), Inclusivity(false, true))",
                ),
            ),
            (
                # Prevent confusion when dealing with time zones by ensuring
                # that the full date and time are displayed
                AnchoredInterval{Day(-1)}(
                    ceil(ZonedDateTime(dt, tz"America/Winnipeg"), Day)
                ),
                "(DE 2016-08-12 00:00:00-05:00]",
                string(
                    "AnchoredInterval{$(repr(Day(-1))),$ZonedDateTime}",
                    "($(repr(ZonedDateTime(2016, 8, 12, tz"America/Winnipeg"))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Minute(-5)}(dt),
                "(2016-08-11 5ME02:00]",
                string(
                    "AnchoredInterval{$(repr(Minute(-5))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Second(-30)}(dt),
                "(2016-08-11 30SE02:00:00]",
                string(
                    "AnchoredInterval{$(repr(Second(-30))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2))), Inclusivity(false, true))",
                ),
            ),
            (
                AnchoredInterval{Millisecond(-10)}(dt),
                "(2016-08-11 10msE02:00:00.000]",
                string(
                    "AnchoredInterval{$(repr(Millisecond(-10))),DateTime}",
                    "($(repr(DateTime(2016, 8, 11, 2))), Inclusivity(false, true))",
                ),
           ),
        ]

        for (interval, printed, shown) in tests
            shown = mod_prefix * shown

            @test sprint(print, interval) == printed
            @test string(interval) == printed
            @test sprint(show, interval) == shown
            @test sprint(show, interval, context=:compact=>true) == printed
            @test repr(interval) == shown
        end

        interval = AnchoredInterval{Second(-10)}(Time(1, 0, 0))
        @test_nowarn string(interval)

        # Non-period AnchoredIntervals
        interval = AnchoredInterval{-10}(10)
        @test string(interval) == "(0 .. 10]"
        @test sprint(show, interval, context=:compact=>true) == string(interval)
        @test sprint(show, interval) ==
            mod_prefix * "AnchoredInterval{-10,$Int}(10, Inclusivity(false, true))"

        interval = AnchoredInterval{25}('a')
        @test string(interval) == "[a .. z)"
        @test sprint(show, interval, context=:compact=>true) == string(interval)
        @test sprint(show, interval) ==
            mod_prefix * "AnchoredInterval{25,Char}('a', Inclusivity(true, false))"
    end

    @testset "equality" begin
        he = HourEnding(dt)
        hb = HourBeginning(dt)
        me = AnchoredInterval{Minute(-1)}(dt)
        cpe = copy(he)
        cpb = copy(hb)

        @test he == cpe
        @test isequal(he, cpe)
        @test hash(he) == hash(cpe)

        @test he != hb
        @test !isequal(he, hb)
        @test hash(he) != hash(hb)

        @test he != me
        @test !isequal(he, me)
        @test hash(he) != hash(me)

        diff_inc = HourEnding(dt, Inclusivity(true, true))

        @test he != diff_inc
        @test !isequal(he, diff_inc)
        @test hash(he) != hash(diff_inc)

        # Overlap for an instant, so not disjoint
        @test isless(he, hb)
        @test !isless(hb, he)
        @test he < hb
        @test !(hb < he)
        @test !(he > hb)
        @test hb > he
        @test !(he ≪ hb)
        @test !(hb ≪ he)
        @test !(he ≫ hb)
        @test !(hb ≫ he)

        @test he ≪ he + Hour(1)
        @test hb ≪ hb + Hour(1)

        @test isless(diff_inc, diff_inc + Hour(2))
        @test isless(diff_inc, diff_inc + Hour(1))
        @test diff_inc < diff_inc + Hour(2)
        @test diff_inc < diff_inc + Hour(1)
        @test diff_inc ≪ diff_inc + Hour(2)
        @test !(diff_inc ≪ diff_inc + Hour(1))  # Overlap for an instant

        # DST transition
        hour1 = HourEnding(ZonedDateTime(2018, 11, 4, 1, tz"America/Winnipeg", 1))
        hour2 = HourEnding(ZonedDateTime(2018, 11, 4, 1, tz"America/Winnipeg", 2))
        @test hour1 != hour2
        @test !isequal(hour1, hour2)
        @test hash(hour1) != hash(hour2)
        @test hour1 ≪ hour2

        # Comparisons between AnchoredInterval{P, T} and T
        @test dt - Hour(2) < HourEnding(dt)
        @test dt - Hour(1) < HourEnding(dt, Inclusivity(false, false))
        @test !(dt - Hour(1) < HourEnding(dt, Inclusivity(true, true)))
        @test !(dt - Minute(30) < HourEnding(dt))
        @test !(dt < HourEnding(dt, Inclusivity(false, false)))
        @test !(dt < HourEnding(dt, Inclusivity(true, true)))
        @test !(dt + Hour(1) < HourEnding(dt))

        @test !(HourEnding(dt) < dt - Hour(2))
        @test !(HourEnding(dt, Inclusivity(false, false)) < dt - Hour(1))
        @test !(HourEnding(dt, Inclusivity(true, true)) < dt - Hour(1))
        @test HourEnding(dt) < dt - Minute(30)
        @test !(HourEnding(dt) ≪ dt - Minute(30))
        @test HourEnding(dt, Inclusivity(false, false)) < dt
        @test HourEnding(dt, Inclusivity(true, true)) < dt
        @test !(HourEnding(dt, Inclusivity(true, true)) ≪ dt)
        @test HourEnding(dt) < dt + Hour(1)
    end

    @testset "broadcasting" begin
        # Validate that an AnchoredInterval is treated as a scalar during broadcasting
        interval = HourEnding(DateTime(2016, 8, 11, 18))
        @test size(interval .== interval) == ()
    end

    @testset "sort" begin
        hb1 = HourBeginning(dt)
        he1 = HourEnding(dt)
        he2 = HourEnding(dt + Hour(1))
        he3 = HourEnding(dt + Hour(2))

        @test sort([hb1, he1, he2, he3]) == [he1, hb1, he2, he3]
        @test sort([hb1, he1, he2, he3]; rev=true) == [he3, he2, hb1, he1]
    end

    @testset "arithmetic" begin
        he = HourEnding(dt)
        hb = HourBeginning(dt)

        @test he + Hour(7) == HourEnding(dt + Hour(7))
        @test hb + Hour(7) == HourBeginning(dt + Hour(7))

        @test Month(1) + he == HourEnding(dt + Month(1))
        @test Month(1) + hb == HourBeginning(dt + Month(1))

        @test he - Day(2) == HourEnding(dt - Day(2))
        @test hb - Day(2) == HourBeginning(dt - Day(2))

        @test he + Minute(30) == HourEnding(dt + Minute(30))
        @test he - Minute(30) == HourEnding(dt - Minute(30))

        ai = AnchoredInterval{Minute(-60)}(dt)
        @test ai + Minute(30) == AnchoredInterval{Minute(-60)}(dt + Minute(30))
        @test ai - Minute(30) == AnchoredInterval{Minute(-60)}(dt - Minute(30))

        # Subtracting AnchoredInterval{P, T} from T doesn't make sense if T is a TimeType
        @test_throws MethodError Hour(1) - he
        @test_throws MethodError Hour(1) - hb

        # DST transition
        endpoint = ZonedDateTime(2018, 3, 11, 3, tz"America/Winnipeg")
        interval = AnchoredInterval{Hour(-2)}(endpoint)
        @test span(interval) == Hour(2)

        endpoint = ZonedDateTime(2018, 11, 4, 2, tz"America/Winnipeg")
        interval = AnchoredInterval{Hour(-2)}(endpoint)
        @test span(interval) == Hour(2)

        # Non-period AnchoredIntervals
        @test AnchoredInterval{10}('a') + 2 == AnchoredInterval{10}('c')
        @test AnchoredInterval{10}('d') - 2 == AnchoredInterval{10}('b')

        @test AnchoredInterval{-1}(20) + 2 == AnchoredInterval{-1}(22)
        @test AnchoredInterval{-1}(20) - 2 == AnchoredInterval{-1}(18)

        @test -AnchoredInterval{2}(10) == AnchoredInterval{-2}(-10) # -[10,12)==(-12,-10]
        @test -AnchoredInterval{-2}(10) == AnchoredInterval{2}(-10) # -(8,10]==[-10,-8)

        @test 15 - AnchoredInterval{-2}(10) == AnchoredInterval{2}(5)   # 15-(8,10]==[5,8)
        @test 15 - AnchoredInterval{2}(10) == AnchoredInterval{-2}(5)   # 15-[10,12)==(3,5]

        @test_throws MethodError -AnchoredInterval{10}('a')
        @test_throws MethodError -AnchoredInterval{-10}('z')

        @test_throws MethodError 10 - AnchoredInterval{10}('a')
        @test_throws MethodError 10 - AnchoredInterval{-10}('z')
    end

    @testset "range" begin
        r1 = HourEnding(dt):HourEnding(dt + Day(1))
        r2 = HourEnding(dt):Hour(2):HourEnding(dt + Day(7))

        @test r1 isa StepRange
        @test length(r1) == 25
        @test collect(r1) == map(HourEnding, dt:Hour(1):dt + Day(1))

        @test r2 isa StepRange
        @test length(r2) == 12 * 7 + 1
        @test collect(r2) == map(HourEnding, dt:Hour(2):dt + Day(7))

        # DST transition
        spring = ZonedDateTime(2018, 3, 11, tz"America/Winnipeg")
        fall = ZonedDateTime(2018, 11, 4, tz"America/Winnipeg")
        r3 = HourEnding(spring):HourEnding(spring + Day(1))
        r4 = HourEnding(fall):Hour(1):HourEnding(fall + Day(1))

        @test r3 isa StepRange
        @test length(r3) == 24
        @test collect(r3) == map(HourEnding, spring:Hour(1):spring + Day(1))

        @test r4 isa StepRange
        @test length(r4) == 26
        @test collect(r4) == map(HourEnding, fall:Hour(1):fall + Day(1))

        r5 = AnchoredInterval{-1}(3, false, true):2:AnchoredInterval{-1}(7, true, true)
        @test length(r5) == 3
        @test collect(r5) == [
            AnchoredInterval{-1}(3, false, true),
            AnchoredInterval{-1}(5, false, true),
            AnchoredInterval{-1}(7, false, true),
        ]
    end

    @testset "isempty" begin
        for P in [Year(1), Month(1), Day(1), Hour(1), Minute(1), Second(1)]
            for sign in [+, -]
                @test !isempty(AnchoredInterval{sign(P)}(dt, Inclusivity(false, false)))
                @test !isempty(AnchoredInterval{sign(P)}(dt, Inclusivity(false, true)))
                @test !isempty(AnchoredInterval{sign(P)}(dt, Inclusivity(true, false)))
                @test !isempty(AnchoredInterval{sign(P)}(dt, Inclusivity(true, true)))
            end
        end

        for P in [Year(0), Month(0), Day(0), Hour(0), Minute(0), Second(0)]
            @test isempty(AnchoredInterval{P}(dt, Inclusivity(false, false)))
            @test isempty(AnchoredInterval{P}(dt, Inclusivity(false, true)))
            @test isempty(AnchoredInterval{P}(dt, Inclusivity(true, false)))
            @test !isempty(AnchoredInterval{P}(dt, Inclusivity(true, true)))
        end
    end

    @testset "in" begin
        @test !in(dt + Hour(1), HourEnding(dt))
        @test in(dt, HourEnding(dt))
        @test !in(dt, HourEnding(dt, Inclusivity(true, false)))
        @test in(dt - Minute(30), HourEnding(dt))
        @test !in(dt - Hour(1), HourEnding(dt))
        @test in(dt - Hour(1), HourEnding(dt, Inclusivity(true, false)))
        @test !in(dt - Hour(2), HourEnding(dt))

        @test !in(dt - Hour(1), HourBeginning(dt))
        @test in(dt, HourBeginning(dt))
        @test !in(dt, HourBeginning(dt, Inclusivity(false, true)))
        @test in(dt + Minute(30), HourBeginning(dt))
        @test !in(dt + Hour(1), HourBeginning(dt))
        @test in(dt + Hour(1), HourBeginning(dt, Inclusivity(false, true)))
        @test !in(dt + Hour(2), HourBeginning(dt))

        zdt = ZonedDateTime(dt, tz"America/Winnipeg")
        @test in(ZonedDateTime(dt - Minute(30), tz"America/Winnipeg"), HourEnding(zdt))
        @test !in(ZonedDateTime(dt + Minute(30), tz"America/Winnipeg"), HourEnding(zdt))
        @test !in(ZonedDateTime(dt - Hour(1), tz"America/Winnipeg"), HourEnding(zdt))
        @test !in(ZonedDateTime(dt - Minute(30), tz"UTC"), HourEnding(zdt))
        @test in(ZonedDateTime(dt + Minute(270), tz"UTC"), HourEnding(zdt))
    end

    @testset "intersect" begin
        # Adjacent
        @test isempty(intersect(HourEnding(dt), HourEnding(dt + Hour(1))))

        # Single point overlap
        expected = AnchoredInterval{Hour(0)}(dt, Inclusivity(true, true))
        @test intersect(
            HourEnding(dt, Inclusivity(true, true)),
            HourEnding(dt + Hour(1), Inclusivity(true, true)),
        ) == expected

        # Hour overlap
        he = HourEnding(dt)
        @test intersect(he, AnchoredInterval{Hour(-2)}(dt)) == he
        @test intersect(AnchoredInterval{Hour(-3)}(dt + Hour(1)), he) == he
        @test intersect(HourBeginning(dt - Hour(1)), he) ==
            HourBeginning(dt - Hour(1), Inclusivity(false, false))

        # Identical save for inclusivity
        expected = HourEnding(dt, Inclusivity(false, false))
        @test intersect(
            HourEnding(dt, Inclusivity(false, false)),
            HourEnding(dt, Inclusivity(true, true)),
        ) == expected
        @test intersect(
            HourEnding(dt, Inclusivity(false, true)),
            HourEnding(dt, Inclusivity(true, false)),
        ) == expected

        # This should probably be an AnchoredInterval{Hour(0)}, but it's not important
        @test intersect(HourEnding(dt), HourBeginning(dt)) ==
            AnchoredInterval{Hour(0)}(dt, Inclusivity(true, true))

        # Non-period AnchoredIntervals
        @test intersect(AnchoredInterval{-2}(3), AnchoredInterval{-2}(4)) ==
            AnchoredInterval{-1}(3)
    end

    @testset "canonicalize" begin
        for s in (1, -1)
            @test canonicalize(Day, Millisecond(s * 3600000)) == Hour(s * 1)
            @test canonicalize(Hour, Millisecond(s * 3600000)) == Hour(s * 1)
            @test canonicalize(Minute, Millisecond(s * 3600000)) == Minute(s * 60)
            @test canonicalize(Second, Millisecond(s * 3600000)) == Second(s * 3600)
            @test canonicalize(Millisecond, Millisecond(s * 3600000)) ==
                Millisecond(s * 3600000)
            @test canonicalize(Second, Millisecond(s * 3601000)) == Second(s * 3601)
            @test canonicalize(Millisecond, Millisecond(s * 3600001)) ==
                Millisecond(s * 3600001)
        end

        # Can't promote past 1 week, because who knows how many days/weeks are in a month?
        @test canonicalize(Year, Millisecond(2419200000)) == Week(4)
        @test canonicalize(Year, Millisecond(4233600000)) == Week(7)
    end

    @testset "astimezone" begin
        zdt = ZonedDateTime(2013, 2, 13, 0, 30, tz"America/Winnipeg")

        for inclusivity in Inclusivity.(0:3)
            for tz in (tz"America/Winnipeg", tz"America/Regina", tz"UTC")
                @test isequal(
                    astimezone(HE(zdt, inclusivity), tz),
                    HE(astimezone(zdt, tz), inclusivity),
                )

                @test isequal(
                    astimezone(AnchoredInterval{Day(1)}(zdt, inclusivity), tz),
                    AnchoredInterval{Day(1)}(astimezone(zdt, tz), inclusivity),
                )
            end
        end
    end
end


text:
@assert proj[1,1] == 10
@assert (simpk'*proj)[1] == 10 - 30im
@assert (proj*simpk)[1] == 10 + 30im
@assert simpk'*proj*simpk == (simpk'*simpk)*(simpk'*simpk) == 100

@assert b*(k*b)*k == (b*k)*(b*k)
@assert b*k == 56 - 28im
@assert b*op*k == 2352 - 3136im

@test_approx_eq qubits'*qubits 1
@test_approx_eq norm(ptrace(bitdens, 2)) 1
@test_approx_eq trace(ptrace(bitdens, 2)) 1
@test_approx_eq trace(ptrace(belldens, 1)^2) .5
@test_approx_eq purity(bell) 1

@rep_op " H | n > = 1/√2 * ( | 0 > + (-1)^n * | 1 > )" 0:1
@assert matrep(H, 0:1) == 1/√2 * [1 1 ; 1 -1]
m = 0.5 * [1 1;1 -1]
@test_approx_eq matrep(tensor(H, H), 0:1, 0:1) [m m ; m -m]

@assert ptrace(ptrace(bitdens, 2), 1) == ptrace(ptrace(bitdens, 1), 2)

@assert act_on(bra(0), bell_unbal, 2)[1] == bell_unbal[1,0]
@assert act_on(bra(0), bell_unbal, 1)[1] == bell_unbal[0,1]

@assert anticommute(op,op) == 2 * op^2
@assert commute(op,op) == op^2 - op^2

b = d" 1/√2 * (< 0 | + < 1 |) "
itest = act_on(b, bell_unbal, 2)

@assert itest[0] == b[1]*bell_unbal[0,1]
@assert itest[1] == b[0]*bell_unbal[1,0]


text:
if VERSION >= v"1.6.0"
    using Downloads
    getfile(x) = Downloads.download(x)
else
    getfile(x) = download(x)
end

@testset "inflate" begin
    test_sequence = UInt8.([128,0,0,4,151,144,57,244,6,125,10,27,2,133,176,29,1,1])

    inflated = inflate(Balloons.LZW(), test_sequence)

    @test inflated == [0,0,73,242,14,250,6,250,40,216,40,182,7,128]

    bali_fp = getfile("https://github.com/tlnagy/exampletiffs/blob/master/bali.tif?raw=true")
    io = open(bali_fp)
    seek(io, 8)
    input = read(io, 3666)

    inflated = Balloons.inflate(Balloons.LZW(), input)

    @test length(inflated) == 7975
end


text:
"""
`module  JAC.Tools`
    ... a submodel of JAC that contains all methods for providing a toolbox for -- more o less --
        simple tools for the JAC program.
"""
module Tools

    using Interact

    """
    `struct  Tools.Settings`
        ... defines a type for the settings in estimating double-Auger and autoionization rates.

        + printBeforeComputation  ::Bool   ... True, if all energies and lines are printed before their evaluation.
    """
    struct Settings
        printBeforeComputation    ::Bool
    end

    """
    `Tools.perform(obs::Observable{Any})`  ... performs different tasks due to the given Observable.
    """
    function perform(obs::Observable{Any})
        ##x println(obs[])
        #
        if      obs[] == "Start task"          println("Nothing is implemented here.")
        elseif  obs[] == "Grid calculator"     Tools.taskGridCalculator()
        else                                   error("stop a")
        end
        return( nothing )
    end

    """
    `Tools.select()`  ... select different tools from a menu for which no further input is required.
    """
    function select()
        t1 = "Simple tools JAC without input: "
        b1 = dropdown(["Start task", "Grid calculator"])
        update = button("Update")
        ui = vbox( hbox( pad(1em, t1) ),
                   hbox( pad(1em, b1), pad(1em, update) )
                 )
        Interact.display(ui)
        output = Interact.@map (&update;  observe(b1)[] )
        return( output )
    end

    """
    `Tools.select(dict::Dict)`
        ... select different tools from a menu if proper results::dict are given in terms of a dictionary as obtained
            usually from an Atomic.Computation.
    """
    function select(dict::Dict)
        println("Tools.select(): A menu if a dict comes in; $dict ")
    end

    """
    `Tools.taskGridCalculator()`  ... grid calculator to determine useful grid parameters
    """
    function taskGridCalculator()
        t1 = "Grid calculator: Determine 'one' of the grid parameters if all others are given explicitly by the user. " *
             "There are two forms of the radial grid: "
        t2 = "   (i)    logarithmic grid   r(i)  = rnt * (Base.MathConstants.e^{(i-1)*h} - 1)            ... for bound-state electrons  \n" *
             "   (ii)   log-lin grid       ln(r(i)/rnt+1) +  (h/hp)*r(i) = (i-1)*h    ... for free electrons.  \n\n"
        t3 = "The logarithmic grid is used always if hp = 0. \n" *
             "The use of this 'calculator' has been found useful in Auger and photoionization studies as well as if " *
             "Rydberg-type electrons need to be taken into account in the bound-state density."
        #
        b1 = spinbox(label="rnt  ";    value=2.0e-6)
        b2 = spinbox(label="h     ";   value=5.0e-2)
        b3 = spinbox(label="hp    ";   value=0.0)
        b4 = spinbox(label="r_max ";   value=5.6e+2)
        update = button("Update")
        ui = vbox( hbox( pad(0em, t1) ),
                   hbox( pad(0em, t2) ),
                   hbox( pad(0em, t3) ),
                   hbox( pad(1em, b1), pad(1em, b2), pad(1em, b3), pad(1em, b4), pad(1em, update) )
                 )
        Interact.display(ui)
        Interact.@map (&update; taskGridCalculatorResults(observe(b1)[], observe(b2)[], observe(b3)[], observe(b4)[]) )
        return( nothing )
    end

    """
    `Tools.taskGridCalculatorResults()`  ... prints the results of the grid computations
    """
    function taskGridCalculatorResults(rnt::Float64, h::Float64, hp::Float64, rmax::Float64)
        println("** $rnt  $h  $hp  $rmax")
    end

end # module


text:
using GLMakie
let
    set_theme!(theme_black())
    # https://www.jpl.nasa.gov/edu/pdfs/scaless_reference.pdf
    # 1au = 149597870700 # meters
    distances = [0, 0.39, 0.72, 1, 1.52, 5.2, 9.54, 19.2, 30.06] * 149597870700 # meters
    speed = 299792458 # meters / second
    diameters = [696340 * 2, 4879, 12104, 12756, 6792, 142984, 120536, 51118, 49528] # km
    colors = ["snow4", "navajowhite", "lightskyblue", "orangered", "olive",
        "burlywood", "cyan3", "dodgerblue"]
    names = ["Mercury", "Venus", "Earth", "Mars", "Jupiter", "Saturn", "Uranus", "Neptune"]

    point = Node(Point2f(0, 0))
    tail = Node(Point2f[(0, 0)])
    ms = Node(20.0)
    tempo = Node(0.0)
    msP = diameters / maximum(diameters) * 700

    fig = Figure(resolution = (1600, 800))
    axS = Axis(fig[1, 1])
    axP = Axis(fig[1, 2]; title = "The speed of light\nin the solar system", titlesize = 35)
    scatter!(axS, [-3.9], [0]; markersize = msP[1], color = :white)
    scatter!(axP, distances[2:end], fill(0, length(distances[2:end]));
        markersize = msP[2:end], color = colors)
    colsize!(fig.layout, 1, Auto(0.1))
    colgap!(fig.layout, 0)
    hidedecorations!(axP; grid = false)
    hidexdecorations!(axS)
    hideydecorations!(axS; grid = false)
    hidespines!(axP)
    hidespines!(axS)
    xlims!(axS, -1, 1)
    ylims!(axS, -1, 1)
    Label(fig[1, 2, TopLeft()], "so fast! 😃", padding = (0, -150, 0, 0), textsize = 20)
    Label(fig[1, 2, TopRight()], "so slow 😢", padding = (-150, 0, 0, 0), textsize = 20)
    Label(fig[2, :], "Twitter: @LazarusAlon", color = :white,
        tellwidth = false, textsize = 30, halign = :right,
        font = "noto-sands-bold")
    Label(fig[2, :], "using Makie", color = "#f7c530",
        tellwidth = false, textsize = 35, halign = :left,
        font = "noto-sands-bold")
    text!(axP, names, position = tuple.(distances[2:end], [1, -1, 0.5, -0.5, 0.5, 0.5, 0.5, 0.5]))
    stem!(axP, distances[2:end], [1, -1, 0.5, -0.5, 0.5, 0.5, 0.5, 0.5]; color = colors,
        stemcolor = colors, trunkcolor = :white, trunklinestyle = :dashdotdot,
        trunkwidth = 1, stemwidth = 0.85, markersize = 8)

    for i = 1:5
        lines!(axP, @lift(length($tail) < 5 ? $tail : $tail[end:-1:end-Int64(round(0.01 * i * length($tail)))]);
            linewidth = 5 - i * 0.7, color = (:white, 1 - 0.1 * i))
    end
    scatter!(axP, point; color = (:white, 0.85), markersize = 20, marker = '▶')
    text!(axP, @lift("$($tempo) hours"), position = point; color = :yellow,
        align = (:center, :bottom))

    record(fig, joinpath(@__DIR__, "output", "flying.mp4"), framerate = 30) do io
        for frame in LinRange(0.001, distances[end] + distances[5], 1000)
            new_point = Point2f(frame, 0.0)
            push!(tail[], new_point)
            point[] = new_point
            tempo[] = round((frame / speed) / 60 / 60, digits = 2)
            tail[] = tail[]
            xlims!(axP, 0, frame + frame * 0.08)
            recordframe!(io)
        end
    end
    set_theme!() # reset to default
end
using Pkg # HIDE
Pkg.status(["GLMakie"]) # HIDE

max_stars_repo_name: lazarusA/BeautifulMakie

max_issues_repo_name: lazarusA/BeautifulMakie

max_forks_repo_name: lazarusA/BeautifulMakie