text:
function upload_tripdata_to_s3(bucket_name)
    verify_file_in_s3(
        bucket_name,
        "tripdata.csv",
        "https://s3.amazonaws.com/nyc-tlc/trip+data/yellow_tripdata_2016-01.csv",
    )
end

@testset "Complex usage of BanyanDataFrames on medium dataset" begin
    run_with_job("Filter groupby") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        upload_tripdata_to_s3(bucket)
        tripdata = read_csv("s3://$(bucket)/tripdata.csv")

        # Read data and get properties
        tripdata = read_csv("s3://$(bucket)/tripdata.csv")
        tripdata = dropmissing(tripdata)
        @test nrow(tripdata) == 10906858
        @test ncol(tripdata) == 19
        @test size(tripdata) == (10906858, 19)
        @test collect(names(tripdata))[1] == "VendorID"
        @test collect(propertynames(tripdata))[1] == :VendorID

        # Get all trips with distance longer than 1.0, group by passenger count,
        # and get the average trip distance for each group
        distances = collect(
            sort(
                combine(
                    groupby(
                        filter(row -> row.trip_distance > 1.0, tripdata),
                        :passenger_count,
                    ),
                    :trip_distance => mean,
                ),
                :trip_distance_mean,
            ),
        )
        @test distances[:, :passenger_count] == [6, 3, 5, 2, 1, 7, 0, 9, 8, 4]
        @test round.(distances[:, :trip_distance_mean], digits = 3) ==
              [3.631, 3.725, 3.729, 5.296, 6.328, 7.15, 8.195, 8.276, 12.091, 19.087]
    end

    run_with_job("Operations on DateTime") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        upload_tripdata_to_s3(bucket)
        tripdata = read_csv("s3://$(bucket)/tripdata.csv")

        # Compute the hour of the day which has the greatest average trip distance.
        setindex!(
            tripdata,
            :start_time,
            map(
                t -> DateTime(t, "yyyy-mm-dd HH:MM:SS"),
                tripdata[:, :tpep_pickup_datetime],
            ),
        )
        setindex!(
            tripdata,
            :end_time,
            map(
                t -> DateTime(t, "yyyy-mm-dd HH:MM:SS"),
                tripdata[:, :tpep_dropoff_datetime],
            ),
        )
        setindex!(tripdata, :start_hour, map(t -> hour(t), tripdata[:, :start_time]))
        tripdata_grouped = groupby(tripdata, :start_hour)
        means = combine(tripdata_grouped, :trip_distance => mean)
        means_sorted = collect(sort(means, :trip_distance_mean))
        @test round.(means_sorted[:, :trip_distance_mean], digits = 3) == [
            2.523,
            2.529,
            2.577,
            2.597,
            2.627,
            2.753,
            2.831,
            2.875,
            2.893,
            2.904,
            3.048,
            3.247,
            3.248,
            3.329,
            3.381,
            3.469,
            3.527,
            3.824,
            4.109,
            4.603,
            5.872,
            6.081,
            12.34,
            21.456,
        ]
        @test means_sorted[:, :start_hour] == [
            9,
            8,
            18,
            12,
            19,
            13,
            7,
            15,
            16,
            20,
            21,
            1,
            2,
            0,
            23,
            3,
            6,
            17,
            4,
            5,
            11,
            22,
            10,
            14,
        ]
        @test means_sorted[:, :start_hour][end] == 14

        # Get the first trip for each pickup time, sort by drop off time, and filter for trips with start times on/after Jan. 30.
        tripdata_unique_start = sort(unique(tripdata, :start_time), :end_time)
        tripdata_unique_start_end_of_month = tripdata_unique_start[
            map(
                st -> st > DateTime(2016, 1, 30, 00, 00, 00, 00),
                tripdata_unique_start[:, :start_time],
            ),
            :,
        ]
        @test nrow(tripdata_unique_start) == 2368616
        @test nrow(tripdata_unique_start_end_of_month) == 163601
        unique_start_sum = collect(round(sum(tripdata_unique_start_end_of_month)))
        @test unique_start_sum == 490661
        pickup_res = collect(tripdata_unique_start_end_of_month[:, :tpep_pickup_datetime])
        @test pickup_res[1] == "2016-01-30 00:00:01"
        @test pickup_res[2] == "2016-01-30 00:00:02"
        @test pickup_res[163593] == "2016-01-31 23:57:26"
    end

    run_with_job("Groups and joins") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        upload_tripdata_to_s3(bucket)
        tripdata = read_csv("s3://$(bucket)/tripdata.csv")

        # Group by passenger count and compute min-max normalized trip distance. Get average for each group.
        # Sort by the average.
        gdf = groupby(tripdata, :passenger_count)

        # Method 1
        result_1 = transform(
            gdf,
            :trip_distance => x -> (x .- minimum(x)) ./ (maximum(x) - minimum(x)),
        )
        result_1 = rename(result_1, :trip_distance_function => :trip_distance_normalized)
        result_1 = sort(
            combine(groupby(result_1, :passenger_count), :trip_distance_normalized => mean),
            :trip_distance_normalized_mean,
        )

        # Method 2
        min_max = combine(gdf, :trip_distance => minimum, :trip_distance => maximum)
        result_2 = innerjoin(tripdata, min_max, on = :passenger_count)
        result_2[:, :trip_distance_normalized] = map(
            (d, min, max) -> (d - min) / (max - min),
            result_2[:, :trip_distance],
            result_2[:, :trip_distance_minimum],
            result_2[:, :trip_distance_maximum],
        )
        result_2 = sort(
            combine(groupby(result_2, :passenger_count), :trip_distance_normalized => mean),
            :trip_distance_normalized_mean,
        )

        result_1 = collect(result_1)
        result_2 = collect(result_2)
        @test round.(result_1[:, :trip_distance_normalized_mean], digits = 3) ==
              round.(result_2[:, :trip_distance_normalized_mean], digits = 3) ==
              [0.0, 0.0, 0.0, 0.013, 0.021, 0.031, 0.042, 0.112, 0.146, 0.181]
    end

    run_with_job("Select and subset") do job
        bucket = get_cluster_s3_bucket_name(get_cluster_name())
        upload_tripdata_to_s3(bucket)
        tripdata = read_csv("s3://$(bucket)/tripdata.csv")

        # Get trips which are longer than the average trip length for the day.
        setindex!(
            tripdata,
            :day,
            map(
                t -> day(DateTime(t, "yyyy-mm-dd HH:MM:SS")),
                tripdata[:, :tpep_pickup_datetime],
            ),
        )
        gdf = groupby(tripdata, :day)
        res = subset(gdf, :trip_distance => d -> d .> mean(d))

        # Save these trips to s3, with the pickup time, dropoff time, average trip distance for each day.
        res_2 = select(
            groupby(res, :day),
            :tpep_pickup_datetime,
            :tpep_dropoff_datetime,
            :trip_distance .=> mean,
        )
        write_csv(res_2, "s3://$(bucket)/tripdata_new.csv")

        @test sort(collect(combine(groupby(res, :passenger_count), nrow))[:, :nrow]) ==
              [3, 6, 11, 69, 53536, 88698, 109032, 149283, 403749, 1813376]
        @test sort(collect(combine(groupby(res_2, :day)))[:, :day]) == [
            120,
            2477,
            25868,
            26725,
            27639,
            40623,
            71419,
            80937,
            83528,
            85354,
            86604,
            88022,
            91998,
            95915,
            96279,
            97655,
            97973,
            98359,
            100627,
            100896,
            103117,
            105021,
            105257,
            106580,
            106687,
            106757,
            110904,
            112558,
            116446,
            117609,
            127809,
        ]
    end
end

max_stars_repo_name: banyan-team/banyan-julia

max_issues_repo_name: banyan-team/banyan-julia

max_forks_repo_name: banyan-team/banyan-julia


text:
@testset "labelmap" begin
    @test_throws ArgumentError labelmap(rand(2,2))
    @testset "Symbol" begin
        lm = @inferred labelmap([:yes,:no,:no,:yes,:yes])
        @test typeof(labelenc(lm)) <: LabelEnc.NativeLabels{Symbol,2}
        @test typeof(lm) <: Dict{Symbol,Vector{Int}}
        @test sort(@inferred(label(lm))) == sort([:yes, :no])
        @test typeof(label(lm)) <: Vector{Symbol}
        @test nlabel(lm) === 2
        @test lm[:yes] == [1,4,5]
        @test lm[:no] == [2,3]
        @test @inferred(labelmap!(lm, 6:7, [:yes,:maybe])) === lm
        @test nlabel(lm) === 3
        @test lm[:yes] == [1,4,5,6]
        @test lm[:no] == [2,3]
        @test lm[:maybe] == [7]
        @test @inferred(labelmap!(lm, 8, :no)) === lm
        @test nlabel(lm) === 3
        @test lm[:yes] == [1,4,5,6]
        @test lm[:no] == [2,3,8]
        @test lm[:maybe] == [7]
        @test @inferred(labelmap!(lm, 9, "no")) === lm
        @test nlabel(lm) === 3
        @test lm[:yes] == [1,4,5,6]
        @test lm[:no] == [2,3,8,9]
        @test lm[:maybe] == [7]
    end
    @testset "Float64" begin
        lm = @inferred labelmap([1.,-1,-1,1,1])
        @test typeof(labelenc(lm)) <: LabelEnc.MarginBased{Float64}
        @test typeof(lm) <: Dict{Float64,Vector{Int}}
        @test @inferred(label(lm)) == [1,-1]
        @test typeof(label(lm)) <: Vector{Float64}
        @test nlabel(lm) === 2
        @test lm[1] == [1,4,5]
        @test lm[-1] == [2,3]
        @test @inferred(labelmap!(lm, 6:7, [1,2])) === lm
        @test nlabel(lm) === 3
        @test lm[1] == [1,4,5,6]
        @test lm[-1] == [2,3]
        @test lm[2] == [7]
        @test @inferred(labelmap!(lm, 8, UInt8(2))) === lm
        @test lm[1] == [1,4,5,6]
        @test lm[-1] == [2,3]
        @test lm[2] == [7,8]
        @test_throws MethodError labelmap!(lm, 8, :no)
        @test nlabel(lm) === 3
        @test lm[1] == [1,4,5,6]
        @test lm[-1] == [2,3]
        @test lm[2] == [7,8]
        @test_throws MethodError labelmap!(lm, 9, "no")
        @test nlabel(lm) === 3
        @test lm[1] == [1,4,5,6]
        @test lm[-1] == [2,3]
        @test lm[2] == [7,8]
        @test_throws MethodError labelmap!(lm, 9:10, ["no","yes"])
        @test nlabel(lm) === 3
        @test lm[1] == [1,4,5,6]
        @test lm[-1] == [2,3]
        @test lm[2] == [7,8]
    end
end

@testset "labelfreq" begin
    @test_throws ArgumentError labelfreq(rand(2,2))
    @testset "from labelmap" begin
        lm = @inferred labelfreq(labelmap([:yes,:no,:no,:yes,:yes]))
        @test typeof(labelenc(lm)) <: LabelEnc.NativeLabels{Symbol,2}
        @test typeof(lm) <: Dict{Symbol,Int}
        @test sort(@inferred(label(lm))) == sort([:yes, :no])
        @test typeof(label(lm)) <: Vector{Symbol}
        @test nlabel(lm) === 2
        @test lm[:yes] == 3
        @test lm[:no] == 2
    end
    @testset "Symbol" begin
        lm = @inferred labelfreq([:yes,:no,:no,:yes,:yes])
        @test typeof(labelenc(lm)) <: LabelEnc.NativeLabels{Symbol,2}
        @test typeof(lm) <: Dict{Symbol,Int}
        @test sort(@inferred(label(lm))) == sort([:yes, :no])
        @test typeof(label(lm)) <: Vector{Symbol}
        @test nlabel(lm) === 2
        @test lm[:yes] == 3
        @test lm[:no] == 2
        @test @inferred(labelfreq!(lm, [:yes,:maybe])) === lm
        @test nlabel(lm) === 3
        @test lm[:yes] == 4
        @test lm[:no] == 2
        @test lm[:maybe] == 1
        @test @inferred(labelfreq!(lm, :no)) === lm
        @test nlabel(lm) === 3
        @test lm[:yes] == 4
        @test lm[:no] == 3
        @test lm[:maybe] == 1
        @test @inferred(labelfreq!(lm, "no")) === lm
        @test nlabel(lm) === 3
        @test lm[:yes] == 4
        @test lm[:no] == 4
        @test lm[:maybe] == 1
    end
    @testset "Float64" begin
        lm = @inferred labelfreq([1.,-1,-1,1,1])
        @test typeof(labelenc(lm)) <: LabelEnc.MarginBased{Float64}
        @test typeof(lm) <: Dict{Float64,Int}
        @test @inferred(label(lm)) == [1,-1]
        @test typeof(label(lm)) <: Vector{Float64}
        @test nlabel(lm) === 2
        @test lm[1] == 3
        @test lm[-1] == 2
        @test @inferred(labelfreq!(lm, [1,2])) === lm
        @test nlabel(lm) === 3
        @test lm[1] == 4
        @test lm[-1] == 2
        @test lm[2] == 1
        @test @inferred(labelfreq!(lm, UInt8(2))) === lm
        @test lm[1] == 4
        @test lm[-1] == 2
        @test lm[2] == 2
        @test_throws MethodError labelfreq!(lm, :no)
        @test nlabel(lm) === 3
        @test lm[1] == 4
        @test lm[-1] == 2
        @test lm[2] == 2
        @test_throws MethodError labelfreq!(lm, "no")
        @test nlabel(lm) === 3
        @test lm[1] == 4
        @test lm[-1] == 2
        @test lm[2] == 2
        @test_throws MethodError labelfreq!(lm, ["no","yes"])
        @test nlabel(lm) === 3
        @test lm[1] == 4
        @test lm[-1] == 2
        @test lm[2] == 2
    end
end

ext: jl


text:
using LinearAlgebra
using JuMP

function test(n::Int, verbose::Bool)
    model = Model()
    a = rand(n, n)
    @variable(model, x[1:n, 1:n])

    # Compile it
    x * a
    a * x * a

    elapsed = @elapsed x * a
    if verbose
        println("2D Matrix product `x * a`: $(elapsed)")
    end

    elapsed = @elapsed a * x * a
    if verbose
        println("2D Matrix product `a * x * a`: $(elapsed)")
    end

    return
end

# Compile the methods needed.
test(2, false)

for n in [10, 20, 50, 100]
    println("n = $n")
    test(n, true)
end

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
using Test

using StaticArrays

using Polyhedra

include("solvers.jl")
include("dummy_optimizer.jl")

@testset "Redundancy removal" begin
    @testset "LP-based" begin
        @testset "HRepresentation with $(typeof(x))" for (x, y, z, w) in [
            ([1, 0], [0, 1], [1, 1], [0, 0]),
            ([1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [0.0, 0.0]),
            ((@SVector [1, 0]), (@SVector [0, 1]), (@SVector [1, 1]), (@SVector [0, 0]))]
            T = eltype(x)
            AT = typeof(x)
            d = Polyhedra.FullDim(x)
            @testset "HAffineSpace" begin
                for hr in [hrep(typeof(HyperPlane(x, one(T)))[]; d=d),
                           intersect(HyperPlane(w, zero(T)))]
                    @test dim(hr) == 2
                    for h in [detecthlinearity(hr, lp_solver),
                              removehredundancy(hr, lp_solver)]
                        @test dim(h) == 2
                        @test h isa HyperPlanesIntersection{T, AT}
                        @test !hashyperplanes(h)
                        @test !hashalfspaces(h)
                    end
                end
                hr = intersect(HyperPlane(x, one(T)))
                @test dim(hr) == 1
                for h in [detecthlinearity(hr, lp_solver),
                          removehredundancy(hr, lp_solver)]
                    @test dim(h) == 1
                    @test h isa HyperPlanesIntersection{T, AT}
                    @test nhyperplanes(h) == 1
                    @test !hashalfspaces(h)
                end
                hr = intersect(HyperPlane(w, one(T)))
                @test dim(hr) == -1
                for h in [detecthlinearity(hr, lp_solver),
                          removehredundancy(hr, lp_solver)]
                    @test dim(h) == -1
                    @test h isa HyperPlanesIntersection{T, AT}
                    @test nhyperplanes(h) == 3
                    @test !hashalfspaces(h)
                end
            end
            @testset "Intersection" begin
                for hr in [HalfSpace(z + x, 3) ∩ HalfSpace(-x, -one(T)) ∩ HalfSpace(-2y, -2),
                           HalfSpace(x, one(T)) ∩ HalfSpace(y, one(T)) ∩ HalfSpace(-z, -2)]
                    h = detecthlinearity(hr, lp_solver, verbose=2)
                    @test h isa Polyhedra.Intersection{T, AT}
                    @test nhyperplanes(h) == 2
                    @test !hashalfspaces(h)
                end
                for hr in [HalfSpace(z, zero(T)) ∩ HalfSpace(-z, -zero(T)),
                           HalfSpace(x, one(T)) ∩ HalfSpace(-x, -one(T))]
                    h = detecthlinearity(hr, lp_solver, verbose=2)
                    @test h isa Polyhedra.Intersection{T, AT}
                    @test nhyperplanes(h) == 1
                    @test !hashalfspaces(h)
                end
                for hr in [HalfSpace(-x, one(T)) ∩ HalfSpace(-y, one(T)) ∩ HalfSpace(x, one(T)) ∩ HalfSpace(y, one(T)) ∩ HalfSpace(-x, -9),
                           HalfSpace(x, zero(T)) ∩ HalfSpace(-x, -one(T))]
                    h = detecthlinearity(hr, lp_solver, verbose=2)
                    @test h isa Polyhedra.Intersection{T, AT}
                    @test nhyperplanes(h) == 3
                    @test !hashalfspaces(h)
                end
            end
        end
        @testset "VRepresentation with $(typeof(x))" for (x, y, z, w) in [
            ([1, 0], [0, 1], [1, 1], [0, 0]),
            ([1.0, 0.0], [0.0, 1.0], [1.0, 1.0], [0.0, 0.0]),
            ((@SVector [1, 0]), (@SVector [0, 1]), (@SVector [1, 1]), (@SVector [0, 0]))]
            T = eltype(x)
            AT = typeof(x)
            d = Polyhedra.FullDim(x)
            @testset "VLinearSpace" begin
                for vr in [vrep(typeof(Line(x))[]; d=d),
                           convexhull(Line(w))]
                    for v in [detectvlinearity(vr, lp_solver),
                              removevredundancy(vr, lp_solver)]
                        @test v isa LinesHull{T, AT}
                        @test !haspoints(v)
                        @test !hasrays(v)
                        @test !haslines(v)
                    end
                end
                vr = convexhull(Line(x))
                for v in [detectvlinearity(vr, lp_solver),
                          removevredundancy(vr, lp_solver)]
                    @test v isa LinesHull{T, AT}
                    @test npoints(v) == 1
                    @test !hasrays(v)
                    @test nlines(v) == 1
                end
                for vr in [convexhull(Line(x), Line(y)),
                           convexhull(Line(x), Line(y), Line(z)),
                           convexhull(Line(x), Line(z), Line(y)),
                           convexhull(Line(x), Line(z), Line(z))]
                    for v in [detectvlinearity(vr, lp_solver),
                              removevredundancy(vr, lp_solver)]
                        @test v isa LinesHull{T, AT}
                        @test npoints(v) == 1
                        @test !hasrays(v)
                        @test nlines(v) == 2
                    end
                end
                vrepm1 = Polyhedra.emptyspace(vr)
                @test vrepm1 === @inferred removevredundancy(vrepm1, lp_solver)
            end
            @testset "Polyhedra.RaysHull" begin
                vr = vrep(typeof(Line(x))[], typeof(Ray(x))[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(Ray(x))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(Ray(x), Ray(y))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 2
                @test !haslines(rm)
                for vr in [convexhull(Ray(x), Ray(y), Ray(z)),
                           convexhull(Ray(x), Ray(z), Ray(y)),
                           convexhull(Ray(z), Ray(x), Ray(y))]
                    rm = removevredundancy(vr, lp_solver)
                    @test rm isa Polyhedra.RaysHull{T, AT}
                    @test npoints(rm) == 1
                    @test nrays(rm) == 2
                    rs = collect(rays(rm))
                    @test rs[1] == Ray(x)
                    @test rs[2] == Ray(y)
                    @test !haslines(rm)
                end
                vr = convexhull(Ray(x), Ray(y), Ray(-z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 3
                @test !haslines(rm)
                vr = convexhull(Ray(x), Ray(y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 2
                @test nlines(rm) == 1
                vr = convexhull(Ray(x), Ray(-y), Line(z))
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 1
                @test nlines(rm) == 1
                vr = conichull(Ray(x), Ray(y), Ray(-z))
                rm = detectvlinearity(vr, lp_solver)
                @test rm isa Polyhedra.RaysHull{T, AT}
                @test npoints(rm) == 1
                @test nrays(rm) == 0
                @test nlines(rm) == 2
            end
            @testset "Polyhedra.PointsHull" begin
                vr = vrep(typeof(x)[]; d=d)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test !haspoints(rm)
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 2
                @test !hasrays(rm)
                @test !haslines(rm)
                for vr in [convexhull(2x, 2y, z),
                           convexhull(2x, z, 2y),
                           convexhull(z, 2x, 2y)]
                    rm = removevredundancy(vr, lp_solver)
                    @test rm isa Polyhedra.PointsHull{T, AT}
                    @test npoints(rm) == 2
                    ps = collect(points(rm))
                    @test ps[1] == 2x
                    @test ps[2] == 2y
                    @test !hasrays(rm)
                    @test !haslines(rm)
                end
                vr = convexhull(x, y, z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 3
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y, -z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.PointsHull{T, AT}
                @test npoints(rm) == 3
                @test !hasrays(rm)
                @test !haslines(rm)
                vr = convexhull(x, y) + Ray(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 2
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(x, y) + Line(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 2
                @test !hasrays(rm)
                @test nlines(rm) == 1
                vr = convexhull(x, -y) + Ray(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 1
                @test first(points(rm)) == -y
                @test nrays(rm) == 1
                @test !haslines(rm)
                vr = convexhull(x, -y) + Line(z)
                rm = removevredundancy(vr, lp_solver)
                @test rm isa Polyhedra.Hull{T, AT}
                @test npoints(rm) == 1
                @test !hasrays(rm)
                @test nlines(rm) == 1
            end
        end
    end
    @testset "Representation-based" begin
        @testset "isredundant" begin
            @test isredundant(HyperPlane([0, 1], 0) ∩ HalfSpace([1, 1], 1), [0, 0], nl=0)
            @test !isredundant(HalfSpace([1, 1], 1) ∩ HyperPlane([0, 1], 0) ∩ HalfSpace([-1, -1], 1), [1, 0], nl=0)
            #@test !ishredundant(convexhull(SymPoint([1, 0])), HyperPlane([0, 1], 0), d=1)
            #@test isredundant(convexhull(SymPoint([1, 0])), HalfSpace([0, 1], 1), d=1)
            @test isredundant(HalfSpace([1, 1], 2) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0), [1, 1], nl=0)
            @test !isredundant(HalfSpace([1, 1], 2) ∩ HalfSpace([-1, 0], 0) ∩ HalfSpace([0, -1], 0), [1, 1], nl=0, strongly=true)
        end
#    @testset "Non-redundant SymPoint" begin
#        vr = convexhull(SymPoint([-1, 0]), SymPoint([0, 1]))
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1) ∩ HalfSpace([1, 1], 1) ∩ HalfSpace([-1, 1], 1)
#        vr = @inferred Polyhedra.removevredundancy(vr, hr)
#        @test vr isa Polyhedra.PointsHull{2,Int,Vector{Int}}
#        @test collect(sympoints(vr)) == [SymPoint([-1, 0]), SymPoint([0, 1])]
#        @test !haspoints(vr)
#        @test !hasallrays(vr)
#    end
        @testset "Point" begin
            # x - y ≤  1
            # x + y ≥ -1
            hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
            # [-1, 0] is weakly redundant as it belongs to the facet x + y = -1
            @test isredundant(hr, [-1, 0], nl=0)
            @test !isredundant(hr, [-1, 0], strongly=true, nl=0)
            # [0, -1] is not redundant as it is an extreme point
            @test !isredundant(hr, [0, -1], nl=0)
            @test !isredundant(hr, [0, -1], strongly=true, nl=0)
            # [0, 0] is strongly redundant as it is in the relative interior
            @test isredundant(hr, [0, 0], nl=0)
            @test isredundant(hr, [0, 0], strongly=true, nl=0)
            vr = convexhull([-1, 0], [0, -1], [0, 0]) + conichull([1, 1], [-1, 1])
            @test collect(removevredundancy(points(vr), hr, nl=0)) == [[0, -1]]
            @test collect(removevredundancy(points(vr), hr, strongly=true, nl=0)) == [[-1, 0], [0, -1]]
            vrr = Polyhedra.removevredundancy(vr, hr)
            @test collect(points(vrr)) == [[0, -1]]
            @test !haslines(vrr)
            @test collect(rays(vrr)) == [Ray([1, 1]), Ray([-1, 1])]
            vrr = Polyhedra.removevredundancy(vr, hr, strongly=true)
            @test collect(points(vrr)) == [[-1, 0], [0, -1]]
            @test !haslines(vrr)
            @test collect(rays(vrr)) == [Ray([1, 1]), Ray([-1, 1])]
            p = polyhedron(vr, DefaultLibrary{Rational{BigInt}}(Polyhedra.OppositeMockOptimizer))
            Polyhedra.computehrep!(p)
            removevredundancy!(p, strongly=true)
            @test collect(points(p)) == [[-1, 0], [0, -1]]
            @test !haslines(p)
            @test collect(rays(p)) == [Ray([1, 1]), Ray([-1, 1])]
            removevredundancy!(p, planar=false)
            @test collect(points(p)) == [[0, -1]]
            @test !haslines(p)
            @test collect(rays(p)) == [Ray([1, 1]), Ray([-1, 1])]
        end
#    @testset "Split SymPoint" begin
#        vr = convexhull(SymPoint([-1, 0]), SymPoint([0, 1])) + conichull([1, 1], [-1, 1])
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1)
#        vr = Polyhedra.removevredundancy(vr, hr)
#        @test !hassympoints(vr)
#        @test collect(points(vr)) == [[0, -1]]
#        @test !haslines(vr)
#        @test collect(rays(vr)) == [Ray([1, 1]), Ray([-1, 1])]
#
#        vr = convexhull(SymPoint([1, 0]), SymPoint([0, 1])) + conichull([-1, 1])
#        hr = HalfSpace([-1, -1], 1) ∩ HalfSpace([1, -1], 1) ∩ HalfSpace([1, 1], 1)
#        vr = Polyhedra.removevredundancy(vr, hr)
#        @test !hassympoints(vr)
#        @test collect(points(vr)) == [[1, 0], [0, -1]]
#        @test !haslines(vr)
#        @test collect(rays(vr)) == [Ray([-1, 1])]
#    end
        @testset "Solver error in redundancy" begin
            vr = convexhull([1.0, 0.0], [0.0, 1.0])
            err = ErrorException("Solver returned `OTHER_ERROR` when attempting to determine whether `[1.0, 0.0]` of index `Polyhedra.Index{Float64, Vector{Float64}}(1)` is redundant. Solver specific status: Dummy status")
            @test_throws err removevredundancy(vr, DummyOptimizer)
        end
    end
end
@testset "Duplicate removal" begin
    h = removeduplicates(hrep([1 1; -1 -1; 1 0; 0 -1], [1, -1, 1, 0]))
    @test nhyperplanes(h) == 1
    @test first(hyperplanes(h)) == HyperPlane([1, 1], 1)
    @test first(hyperplanes(h)) == HyperPlane([-1, -1], -1)
    @test nhalfspaces(h) == 1
    h = removeduplicates(hrep([1 1; -1 -1], [2, -1]))
    @test !hashyperplanes(h)
    @test nhalfspaces(h) == 2
#    v = removeduplicates(convexhull([1, 2], SymPoint([1, 2]), [-1, 1]))
#    @test typeof(v) == Polyhedra.PointsHull{2,Int,Vector{Int}}
#    @test collect(sympoints(v)) == [SymPoint([1, 2])]
#    @test collect(points(v)) == [[-1, 1]]
#    @test !hasallrays(v)
    for v in (#removeduplicates(convexhull([1, 2], [2, 1], SymPoint([1, 2]), [-1, 1]) + Line([0, 1])),
              #removeduplicates(Line([0, 1]) + convexhull([1, 2], [2, 1], SymPoint([1, 2]), [-1, 1])),
              removeduplicates(Line([0, 1]) + convexhull([-1, -2], [2, 1], [1, 2], [-1, 1], [-1, -2])),)
        @test typeof(v) == Polyhedra.Hull{Int, Vector{Int}, Int}
        #@test collect(sympoints(v)) == [SymPoint([1, 2])]
        @test collect(points(v)) == [[-1, -2], [2, 1], [1, 2]]
        @test !hasrays(v)
        @test collect(lines(v)) == [Line([0, 1])]
    end
end

@testset "Numerically delicate #270" begin
    A = [
        -1.0 0.0 0.0 0.0 0.0 0.0 0.0 0.0
        -0.00106131548486069 0.016634209829441304 -0.0006666769554416961 -1.0 -0.0 -0.0 -0.0 -0.0
        -0.0004116155738065077 0.01693983068952557 -2.6779479610930705e-5 -0.0 -1.0 -0.0 -0.0 -0.0
        -0.0024283354284051484 0.08467041915859919 -0.0005000035065344641 -0.0 -0.0 -1.0 -0.0 -0.0
        -0.002024394302735674 0.08367543404800237 -0.00013321128253368072 -0.0 -0.0 -0.0 -1.0 -0.0
        -0.0009309923368121566 0.009059293773203692 -0.0006675849264721348 -0.0 -0.0 -0.0 -0.0 -1.0
        0.00106131548486069 -0.016634209829441304 0.0006666769554416961 1.0 0.0 0.0 0.0 0.0
        0.0004116155738065077 -0.01693983068952557 2.6779479610930705e-5 0.0 1.0 0.0 0.0 0.0
        0.0024283354284051484 -0.08467041915859919 0.0005000035065344641 0.0 0.0 1.0 0.0 0.0
        0.002024394302735674 -0.08367543404800237 0.00013321128253368072 0.0 0.0 0.0 1.0 0.0
        0.0009309923368121566 -0.009059293773203692 0.0006675849264721348 0.0 0.0 0.0 0.0 1.0
        -0.08320472072273968 2.548630897981053 0.001134866706464998 0.0 0.0 0.0 0.0 0.0
        1210.55555992692 -37126.858723322744 -13.95377106842899 -0.0 -0.0 -0.0 -0.0 -0.0
        -1860.676468633397 57155.38024488611 18.21165578673308 -0.0 -0.0 -0.0 -0.0 -0.0
    ]
    # The last row is almost considered an equality, depends on `ztol`.
    # If it's detected as an equality then the polyhedron is detected as empty.
    # See https://github.com/JuliaPolyhedra/Polyhedra.jl/issues/270
    b = [-20.995, -1.007487377909021, -0.9877340129555137, -0.9705208833525414, -0.9548402631790559, -0.025711479466835108, 1.007487377909021, 0.9877340129555137, 0.9705208833525414, 0.9548402631790559, 0.025711479466835108, 0.8046759919340509, -11751.33121007926, 18148.933864538412]
    lib = DefaultLibrary{Float64}(lp_solver)
    P = polyhedron(hrep(A, b), lib)
    removevredundancy!(P, ztol=1e-4)
    @test nhyperplanes(P) >= 4
    @test nhalfspaces(P) <= 4
    @test npoints(P) > 0
end

import GLPK
@testset "Noise linearity #259" begin
    h = HalfSpace([-1.3333333333333333, -2.220446049250313e-16], 0.0) ∩ HalfSpace([0.0, 1.333333333333333], 0.0) ∩ HalfSpace([-2.7755575615628914e-17, -0.6666666666666667], 0.0) ∩ HalfSpace([-0.33333333333333337, 0.0], 0.0)
    # Presolve off on purpose to test that `λ_l∞` fixes it.
    hh = detecthlinearity(h, GLPK.Optimizer)
    @test nhyperplanes(hh) == 1
    hhh = removehredundancy(hh, GLPK.Optimizer)
    @test nhyperplanes(hhh) == 1
    @test nhalfspaces(hhh) == 1
end

max_stars_repo_name: blegat/Polyhedra

max_issues_repo_name: blegat/Polyhedra

max_forks_repo_name: blegat/Polyhedra


text:
#-----------------------------------------------------------------------
"""
        rms( x::Vector )

Computes the `rms` value for a vectors.
"""
rms=x->sqrt(dot(x,x)/length(x))            # Root Mean Square calculator

#-----------------------------------------------------------------------
"""
        PrintErr(  h::RealVec , v::RealVec  , e::RealVec , thresh::Tuple , pf::String , units::String )

This function prints the errors, their rms and 95% CL values, and the
maximum error.

"""
function PrintErr(  h::RealVec , v::RealVec  , e::RealVec ,
                    thresh::Tuple , pf::String , units::String )
    hrms = rms(h)
    vrms = rms(v)
    erms = rms(e)

    n=length(h)

    h95 = h[Int(0.05*n)]
    v95 = v[Int(0.05*n)]
    e95 = e[Int(0.05*n)]

    hcthr = count(x->x>thresh[1], h)
    vcthr = count(x->x>thresh[2], v)
    ecthr = count(x->x>thresh[3], e)

    open(pf*"-res.txt","a") do io
        println(io,"hrms =","\t",hrms,"\t",
                "# Horizontal error rms (",units,")")
        println(io,"vrms =","\t",vrms,"\t",
                "# Vertical error rms (",units,")")
        println(io,"erms =","\t",erms,"\t",
                "# Full error rms (",units,")")
        println(io,"h95 =","\t",h95,"\t",
                "# Horizontal error 95% CL (",units,")")
        println(io,"v95 =","\t",v95,"\t",
                "# Vertical error 95% CL (",units,")")
        println(io,"e95 =","\t",e95,"\t",
                "# Full error 95% CL (",units,")")
        println(io,"hcthr =","\t",hcthr,"\t",
                "# Horizontal errors above threshold of ",thresh[1]," ",
                units)
        println(io,"vcthr =","\t",vcthr,"\t",
                "# Vertical errors above threshold of ",thresh[2]," ",
                units)
        println(io,"ecthr =","\t",ecthr,"\t",
                "# Full errors above threshold of ",thresh[3]," ",units)
        println(io,"hmax =","\t",h[1],"\t",
                "# Largest horizontal error (",units,")")
        println(io,"vmax =","\t",v[1],"\t",
                "# Largest vertical error (",units,")")
        println(io,"emax =","\t",e[1],"\t",
                "# Largest total error (",units,")")
    end

    return (hrms,vrms,erms,h95,v95,e95,hcthr,vcthr,ecthr)
end     #---------------------------------------------------------------

#-----------------------------------------------------------------------
"""
PlotErr(   err::RealVec  , lines::Tuple , histpar::Tuple , xr::Tuple , yr::Tuple , labels::Tuple , plotname::String )

This function plots the errors.

"""
function PlotErr(   err::RealVec  , lines::Tuple , histpar::Tuple ,
                    xr::Tuple , yr::Tuple , labels::Tuple ,
                    plotname::String )
    plot!(vline([lines[1],lines[2]]))
    histogram!(
        err ,
        alpha=0.50 ,
        bins=histpar[2]:histpar[1]:histpar[3],
        xlims=xr ,
        ylims=yr ,
        xlabel=labels[1] ,
        ylabel=labels[2] ,
        leg=false ,
        fontfamily="Times"
    )
    savefig(plotname)
end     #---------------------------------------------------------------


text:
import .ExpressionExplorer
import .ExpressionExplorer: join_funcname_parts, SymbolsState, FunctionNameSignaturePair

"Return a copy of `old_topology`, but with recomputed results from `cells` taken into account."
function updated_topology(old_topology::NotebookTopology, notebook::Notebook, cells)

	updated_codes = Dict{Cell,ExprAnalysisCache}()
	updated_nodes = Dict{Cell,ReactiveNode}()
	unresolved_cells = copy(old_topology.unresolved_cells.c)
	for cell in cells
		old_code = old_topology.codes[cell]
		if old_code.code !== cell.code
			new_code = updated_codes[cell] = ExprAnalysisCache(notebook, cell)
			new_symstate = new_code.parsedcode |>
				ExpressionExplorer.try_compute_symbolreferences
			new_reactive_node = ReactiveNode(new_symstate)

			updated_nodes[cell] = new_reactive_node
		elseif old_code.forced_expr_id !== nothing
			# reset computer code
			updated_codes[cell] = ExprAnalysisCache(old_code; forced_expr_id=nothing, function_wrapped=false)
		end

		new_reactive_node = get(updated_nodes, cell, old_topology.nodes[cell])
		if !isempty(new_reactive_node.macrocalls)
			# The unresolved cells are the cells for wich we cannot create
			# a ReactiveNode yet, because they contains macrocalls.
			push!(unresolved_cells, cell)
		else
			pop!(unresolved_cells, cell, nothing)
		end
	end
	new_codes = merge(old_topology.codes, updated_codes)
	new_nodes = merge(old_topology.nodes, updated_nodes)

	for removed_cell in setdiff(keys(old_topology.nodes), notebook.cells)
		delete_unsafe!(new_nodes, removed_cell)
		delete_unsafe!(new_codes, removed_cell)
		delete!(unresolved_cells, removed_cell)
	end

	# TODO: this could be cached in a better way if it matters
	cell_order = notebook.cells == old_topology.cell_order ?
		old_topology.cell_order : # if the order is the same, we can reuse the old one
		ImmutableVector(notebook.cells) # (this will do a shallow copy)

	NotebookTopology(;
		nodes=new_nodes,
		codes=new_codes,
		unresolved_cells=ImmutableSet(unresolved_cells; skip_copy=true),
		cell_order
	)
end


text:
using Documenter, WGLMakie, JSServe

makedocs(; modules=[WGLMakie], format=Documenter.HTML(), pages=["Home" => "index.md"],
         repo="https://github.com/kool7d/WGLMakie.jl",
         sitename="WGLMakie.jl", authors="Dan Kool")

deploydocs(; repo="github.com/kool7d/WGLMakie.jl", push_preview=true)


text:
@testset "field" begin

N = [5,6,7]
ncmp = 4

# Test the size of create_field_array().
F = create_field_array(N, ncmp=ncmp)
@test size(F) == (N..., ncmp)

# Test if field_arr2vec() returns a view.
F .= rand(eltype(F), size(F))
f = field_arr2vec(F, order_cmpfirst=false)
f_cmpfirst = field_arr2vec(F, order_cmpfirst=true)

# Test if order_cmpfirst is working as intended.
@test reshape(transpose(reshape(f_cmpfirst, ncmp, :)), :) == f

end  # @testset "field"


text:
using Images, MXNet, JLD

function get_labels_from_file(path)
  return map(x -> parse(Int8, split(x)[2]), readlines(path))
end

PROJECT_IMG_FOLDER = joinpath("data", "book-covers");
DUMP_IMG_TRAIN = joinpath(PROJECT_IMG_FOLDER, "dump-train-inception.jld")
DUMP_IMG_TEST =  joinpath(PROJECT_IMG_FOLDER, "dump-test-inception.jld")
LABELS_TRAIN = joinpath(PROJECT_IMG_FOLDER, "bookcover30-labels-train.txt")
LABELS_TEST = joinpath(PROJECT_IMG_FOLDER, "bookcover30-labels-test.txt")

# get the data
train_img_data = JLD.load(DUMP_IMG_TRAIN, "data");
test_img_data = JLD.load(DUMP_IMG_TEST, "data");

train_labels = get_labels_from_file(LABELS_TRAIN);
test_labels = get_labels_from_file(LABELS_TEST);

# shuffle the dataset
shuffle_train_indices = shuffle(1:size(train_labels, 1));
shuffle_test_indices = shuffle(1:size(test_labels, 1));

train_img_data = train_img_data[:, shuffle_train_indices];
train_labels = train_labels[shuffle_train_indices];

test_img_data = test_img_data[:, shuffle_test_indices];
test_labels = test_labels[shuffle_test_indices];

# create data providers
valid_record_count = 5000
train_data_provider = mx.ArrayDataProvider(:data => train_img_data[:, 1:end-valid_record_count], :label => train_labels[1:end-valid_record_count], shuffle = true, batch_size = 500);
valid_data_provider = mx.ArrayDataProvider(:data => train_img_data[:, end-valid_record_count+1:end], :label => train_labels[end-valid_record_count+1:end], shuffle = true, batch_size = 500);

arch = @mx.chain mx.Variable(:data) =>
  mx.FullyConnected(num_hidden=512) =>
  mx.Activation(name=:relu1, act_type=:relu) =>
  mx.FullyConnected(num_hidden=30) =>
  mx.SoftmaxOutput(mx.Variable(:label))

nnet = mx.FeedForward(arch, context = mx.cpu())
mx.fit(nnet, mx.ADAM(), train_data_provider, eval_data = valid_data_provider, n_epoch = 10, eval_metric = mx.MultiMetric([mx.Accuracy(), mx.ACE()]));
mx.save_checkpoint(nnet, "weights/book-cover-classifier", mx.OptimizationState(1, 0, 0, 0))

# testing the model
test_data_provider = mx.ArrayDataProvider(:data => test_img_data, batch_size = 500);
predictions = mx.predict(nnet, test_data_provider);
results = map(x -> findmax(predictions[:, x])[2], 1:size(test_labels, 1));

println(round(mean((results - 1) .== test_labels), 2))


text:
get_integral(integrals::Vector, integral_map::Dict, symbolic_integral) =
    integrals[integral_map[symbolic_integral]]

function get_or_create_integral!(integrals, integral_map, symbolic_integral, atom; kwargs...)
    if symbolic_integral ∈ keys(integral_map)
        get_integral(integrals, integral_map, symbolic_integral)
    else
        i = create_integral(symbolic_integral, atom, integrals, integral_map; kwargs...)
        isnothing(i) && return
        push!(integrals, i)
        integral_map[symbolic_integral] = length(integrals)
        i
    end
end

function create_integral(symbolic_integral::OrbitalOverlap, atom::Atom, integrals, integral_map; kwargs...)
    a,b = symbolic_integral.a,symbolic_integral.b
    OrbitalOverlapIntegral(a, b, atom)
end

function create_integral(symbolic_integral::CoulombPotentialMultipole, atom::Atom, integrals, integral_map; kwargs...)
    a,b = symbolic_integral.a[1],symbolic_integral.b[1]
    k,g = symbolic_integral.o.k,symbolic_integral.o.g
    # All CoulombPotentialMultipoles that occur as common
    # integrals in an equation are direct potentials acting on
    # either the orbital, or on a source orbital coupled
    # through configuration interaction.
    HFPotential(:direct, k, a, b, atom, g; kwargs...)
end

function create_integral(ome::OrbitalMatrixElement{1,aO,<:OneBodyOperator,bO},
                         atom::Atom{T}, integrals, integral_map; kwargs...) where {aO<:SpinOrbital,bO<:SpinOrbital,T}
    a,o,b = ome.a,ome.o,ome.b
    op = Atoms.get_operator(o, atom, a[1], b[1]; kwargs...)
    iszero(op) && return ZeroIntegral{aO,bO,real(T)}()

    OperatorMatrixElement(a[1], b[1], op, atom, 1)
end

function create_integral(symbolic_integral::OrbitalMatrixElement{2,<:SpinOrbital,<:CoulombInteractionMultipole,<:SpinOrbital},
                         atom::Atom, integrals, integral_map; kwargs...)
    a,b = symbolic_integral.a[1],symbolic_integral.b[1]
    ∂I = diff(symbolic_integral, conj(a))
    # The operator in a two-electron integral is a HFPotential (that
    # may be shared with many other two-electron integrals). The
    # potential must be updated before the two-electron integral is
    # computed, and we achieve this by making sure it appears before
    # the two-electron integral in the list of integrals (which is
    # updated sequentially).
    hfpotential = get_or_create_integral!(integrals, integral_map, ∂I.operator, atom; kwargs...)
    # We know that the coefficient is purely numeric (and likely
    # unity), since we varied a single OrbitalMatrixElement with
    # respect to one of the orbitals.
    coeff = ∂I.factor.coeff
    OperatorMatrixElement(a, b, hfpotential.V̂, atom, coeff)
end

create_integral(symbolic_integral, atom::Atom, integrals, integral_map; kwargs...) =
    throw(ArgumentError("$(symbolic_integral) of type $(typeof(symbolic_integral)) not yet supported"))

function append_common_integrals!(integrals::Vector, integral_map::Dict, atom::Atom, equation_integrals; kwargs...)
    for symbolic_integral in equation_integrals
        get_or_create_integral!(integrals, integral_map, symbolic_integral, atom; kwargs...)
    end
end

"""
    pushterms!(terms, operator, equation_terms,
               integrals, integral_map, symbolic_integrals)

For each term in `equation_terms`, push a term, located at CI
coordinates `i,j`, of the overall orbital Hamiltonian to `terms`,
constructed from `operator` and a product of orbital `integrals`,
multiplied by an overall factor given by expression and multipole
expansions. `integrals` contain common [`OrbitalIntegral`](@ref)s
and `integral_map` maps from `symbolic_integrals` to `integrals`.
"""
function pushterms!(terms::Vector{<:OrbitalHamiltonianTerm{aO,bO,T}},
                    operator::QO,
                    equation_terms::Vector,
                    integrals::Vector,
                    integral_map::Dict{Any,Int},
                    symbolic_integrals) where {aO,bO,T,QO}
    for eq_term in equation_terms
        integrals = OrbitalIntegral{0}[
            get_integral(integrals, integral_map, symbolic_integrals[i])
            for i in eq_term.integrals
        ]
        push!(terms,
              OrbitalHamiltonianTerm{aO,bO,T,QO,typeof(integrals)}(
                  eq_term.i, eq_term.j, T(eq_term.coeff),
                  operator, integrals))
    end
end


text:
"""
    _substitute_views!(vals::Array{T, N}, comp_def) where {T, N}

For each value in `vals`, if the value is a `TimestepArray` swap in a new
TimestepArray with the same type parameterization but with its `data` field
holding a view of the original value's `data` defined by the first and last
indices of Component `comp_def`.
"""
# helper function to substitute views for data
function _substitute_views!(vals::Array{T, N}, comp_def) where {T, N}
    times = [keys(comp_def.dim_dict[:time])...]
    first_idx = findfirst(times .== comp_def.first)
    last_idx = findfirst(times .== comp_def.last)
    for (i, val) in enumerate(vals)
        if val isa TimestepArray
            vals[i] = _get_view(val, first_idx, last_idx)
        end
    end
end

"""
    _get_view(val::TimestepArray{T_TS, T, N, ti, S}, first_idx, last_idx) where {T_TS, T, N, ti, S}

Return a TimestepArray with the same type parameterization as the `val` TimestepArray,
but with its `data` field holding a view of the `val.data` based on the entered
`first-idx` and `last_idx`.
"""
function _get_view(val::TimestepArray{T_TS, T, N, ti, S}, first_idx, last_idx) where {T_TS, T, N, ti, S}

    idxs  = Array{Any}(fill(:, N))
    idxs[ti] = first_idx:last_idx
    # if we are making a connection, the val.data may already be a view, in which case
    # we need to return to the parent of that view before taking our new view
    return TimestepArray{T_TS, T, N, ti}(val.data isa SubArray ? view(val.data.parent, idxs...) : view(val.data, idxs...))
end

"""
    _instance_datatype(md::ModelDef, def::AbstractDatumDef)

Return the datatype of the AbstractDataumDef `def` in ModelDef `md`, which will
be used to create ModelInstance instance variables and parameters.
"""
function _instance_datatype(md::ModelDef, def::AbstractDatumDef)
    dtype = def.datatype == Number ? number_type(md) : def.datatype
    dims = dim_names(def)
    num_dims = dim_count(def)

    ti = get_time_index_position(def)

    if num_dims == 0
        T = ScalarModelParameter{dtype}

    elseif ti === nothing     # there's no time dimension
        T = Array{dtype, num_dims}

    else
        if isuniform(md)
            first, stepsize = first_and_step(md)
            last = last_period(md)
            first === nothing && @warn "_instance_datatype: first === nothing"
            T = TimestepArray{FixedTimestep{first, stepsize, last}, Union{dtype, Missing}, num_dims, ti}
        else
            times = time_labels(md)
            T = TimestepArray{VariableTimestep{(times...,)}, Union{dtype, Missing}, num_dims, ti}
        end
    end

    # @info "_instance_datatype returning $T"
    return T
end

"""
    _instantiate_datum(md::ModelDef, def::AbstractDatumDef)

Return the parameterized datum, broadly either Scalar or Array, pertaining to
AbstractDatumDef `def` in the Model Def `md`, that will support instantiate of parameters
and variables.
"""
# Create the Ref or Array that will hold the value(s) for a Parameter or Variable
function _instantiate_datum(md::ModelDef, def::AbstractDatumDef)
    dtype = _instance_datatype(md, def)
    dims = dim_names(def)
    num_dims = length(dims)

    # Scalar datum
    if num_dims == 0
        value = dtype(0)

    # Array datum, with :time dimension
    elseif dims[1] == :time

        if num_dims == 1
            value = dtype(dim_count(md, :time))
        else
            counts = dim_counts(md, Vector{Symbol}(dims))
            value = dtype <: AbstractArray ? dtype(undef, counts...) : dtype(counts...)
        end

    # Array datum, without :time dimension
    else
        # TBD: Handle unnamed indices properly
        counts = dim_counts(md, Vector{Symbol}(dims))
        value = dtype <: AbstractArray ? dtype(undef, counts...) : dtype(counts...)
    end

    return value
end

"""
    _instantiate_component_vars(md::ModelDef, comp_def::ComponentDef)

Instantiate a component `comp_def` in the model `md` and its variables (but not its
parameters). Return the resulting ComponentInstanceVariables.
"""
function _instantiate_component_vars(md::ModelDef, comp_def::ComponentDef)
    var_defs = variables(comp_def)

    names  = Symbol[nameof(def) for def in var_defs]
    values = Any[_instantiate_datum(md, def) for def in var_defs]
    _substitute_views!(values, comp_def)

    # this line was replaced with the one below because calling _instance_datatype
    # does not concretely type the S type parameter (the type of the TimestepArray's
    # and thus DataType[typeof(val) for val in values] errored when trying to
    # convert typeof(val<:TimestepArray) to a DataType

    # types  = DataType[_instance_datatype(md, def) for def in var_defs]
    types = DataType[typeof(val) for val in values]
    paths  = repeat(Any[comp_def.comp_path], length(names))

    return ComponentInstanceVariables(names, types, values, paths)
end

"""
    function _instantiate_vars(md::ModelDef)

Create the top-level variables for the Model Def `md` and return the dictionary
of the resulting ComponentInstanceVariables.
"""
function _instantiate_vars(md::ModelDef)
    vdict = Dict{ComponentPath, Any}()
    recurse(md, cd -> vdict[cd.comp_path] = _instantiate_component_vars(md, cd); leaf_only=true)
    return vdict
end

"""
    _find_paths_and_names(obj::AbstractComponentDef, datum_name::Symbol)

Recurses through sub components and finds the full path(s) to desired datum, and their
names at the leaf level. Returns a tuple (paths::Vector{ComponentPath}, datum_names::Vector{Symbol})
"""
function _find_paths_and_names(obj::AbstractComponentDef, datum_name::Symbol)

    # Base case-- leaf component
    if obj isa ComponentDef
        return ([nothing], [datum_name])
    end

    datumdef = obj[datum_name]
    if datumdef isa CompositeVariableDef
        refs = [datumdef.ref]   # CompositeVariableDef's can only point to one subcomponent
    else
        refs = datumdef.refs    # ComposteParameterDef's can have multiple refs
    end

    paths = []
    datum_names = []

    for ref in refs
        # Get the comp and datum's for the current ref
        next_obj = obj[ref.comp_name]
        next_datum_name = ref.datum_name

        # Recurse
        sub_paths, sub_datum_names = _find_paths_and_names(next_obj, next_datum_name)

        # Append the paths, and save with datum_names
        for (sp, dn) in zip(sub_paths, sub_datum_names)
            push!(paths, ComponentPath(next_obj.name, sp))
            push!(datum_names, dn)
        end
    end

    return (paths, datum_names)
end

"""
    _get_leaf_level_ipcs(md::ModelDef, conn::InternalParameterConnection)

Returns a vector of InternalParameterConnections that represent all of the connections at the leaf level
that need to be made under the hood as specified by `conn`.
"""
function _get_leaf_level_ipcs(md::ModelDef, conn::InternalParameterConnection)

    top_dst_path = conn.dst_comp_path
    comp = find_comp(md, top_dst_path)
    comp !== nothing || error("Cannot find $(top_dst_path) from $(md.comp_id)")
    par_sub_paths, param_names = _find_paths_and_names(comp, conn.dst_par_name)
    param_paths = [ComponentPath(top_dst_path, sub_path) for sub_path in par_sub_paths]

    top_src_path = conn.src_comp_path
    comp = find_comp(md, top_src_path)
    comp !== nothing || error("Cannot find $(top_src_path) from $(md.comp_id)")
    var_sub_path, var_name = _find_paths_and_names(comp, conn.src_var_name)
    var_path = ComponentPath(top_src_path, var_sub_path[1])

    ipcs = [InternalParameterConnection(var_path, var_name[1], param_path, param_name,
        conn.ignoreunits, conn.backup; backup_offset=conn.backup_offset) for (param_path, param_name) in
        zip(param_paths, param_names)]
    return ipcs
end

"""
    _get_leaf_level_epcs(md::AbstractCompositeComponentDef, epc::ExternalParameterConnection)

Returns a vector that has a new ExternalParameterConnections that represent all of the connections at the leaf level
that need to be made under the hood as specified by `epc`.
"""
function _get_leaf_level_epcs(md::ModelDef, epc::ExternalParameterConnection)

    comp = find_comp(md, epc.comp_path)
    comp !== nothing || error("Cannot find $(epc.comp_path) from $(md.comp_id)")
    par_sub_paths, param_names = _find_paths_and_names(comp, epc.param_name)

    leaf_epcs = ExternalParameterConnection[]
    model_param_name = epc.model_param_name

    top_path = epc.comp_path

    for (par_sub_path, param_name) in zip(par_sub_paths, param_names)
        param_path = ComponentPath(top_path, par_sub_path)
        epc = ExternalParameterConnection(param_path, param_name, model_param_name)
        push!(leaf_epcs, epc)
    end

    return leaf_epcs
end

# generic helper function to get connector component name
connector_comp_name(i::Int) = Symbol("ConnectorComp$i")

"""
    _collect_params(md::ModelDef, var_dict::Dict{ComponentPath, Any})

Collect all parameters in ModelDef `md` with connections to allocated variable
storage in `var_dict` and return a dictionary of (comp_path, par_name) => ModelParameter
elements.
"""
function _collect_params(md::ModelDef, var_dict::Dict{ComponentPath, Any})

    # @info "Collecting params for $(comp_def.comp_id)"

    # Iterate over connections to create parameters, referencing storage in vars
    conns = []
    recurse(md, cd -> append!(conns, internal_param_conns(cd)); composite_only=true)

    pdict = Dict{Tuple{ComponentPath, Symbol}, Any}()

    for conn in conns
        ipcs = _get_leaf_level_ipcs(md, conn)
        src_vars = var_dict[ipcs[1].src_comp_path]
        var_value_obj = get_property_obj(src_vars, ipcs[1].src_var_name)
        for ipc in ipcs
            _check_attributes(md, ipc)
            pdict[(ipc.dst_comp_path, ipc.dst_par_name)] = var_value_obj
        end
    end

    for epc in external_param_conns(md)
        param = model_param(md, epc.model_param_name)
        leaf_level_epcs = _get_leaf_level_epcs(md, epc)
        for leaf_epc in leaf_level_epcs
            pdict[(leaf_epc.comp_path, leaf_epc.param_name)] = (param isa ScalarModelParameter ? param : value(param))
        end
    end

    # Make the external parameter connections for the hidden ConnectorComps.
    # Connect each :input2 to its associated backup value.
    backups = []
    recurse(md, cd -> append!(backups, cd.backups); composite_only=true)

    for (i, backup) in enumerate(backups)
        conn_comp = compdef(md, connector_comp_name(i))
        conn_path = conn_comp.comp_path

        param = model_param(md, backup)
        pdict[(conn_path, :input2)] = (param isa ScalarModelParameter ? param : value(param))
    end

    return pdict
end

"""
    _instantiate_params(comp_def::ComponentDef, par_dict::Dict{Tuple{ComponentPath, Symbol}, Any})

Create the top-level parameters for the Model Def `md` using the parameter dictionary
`par_dict` and return the resulting ComponentInstanceParameters.
"""
function _instantiate_params(comp_def::ComponentDef, par_dict::Dict{Tuple{ComponentPath, Symbol}, Any})
    # @info "Instantiating params for $(comp_def.comp_path)"
    comp_path = comp_def.comp_path
    names = parameter_names(comp_def)
    vals  = Any[par_dict[(comp_path, name)] for name in names]
    _substitute_views!(vals, comp_def)
    types = DataType[typeof(val) for val in vals]
    paths = repeat([comp_def.comp_path], length(names))

    return ComponentInstanceParameters(names, types, vals, paths)
end

"""
    _build(comp_def::ComponentDef,
                    var_dict::Dict{ComponentPath, Any},
                    par_dict::Dict{Tuple{ComponentPath, Symbol}, Any},
                    time_bounds::Tuple{Int, Int})

Return a built leaf or composite LeafComponentInstance created using ComponentDef
`comp_def`, variables and parameters from `var_dict` and `par_dict` and the time
bounds set by `time_bounds`.
"""
function _build(comp_def::ComponentDef,
                var_dict::Dict{ComponentPath, Any},
                par_dict::Dict{Tuple{ComponentPath, Symbol}, Any},
                time_bounds::Tuple{Int, Int})
    # @info "_build leaf $(comp_def.comp_id)"
    # @info "  var_dict $(var_dict)"
    # @info "  par_dict $(par_dict)"

    pars = _instantiate_params(comp_def, par_dict)
    vars = var_dict[comp_def.comp_path]

    return LeafComponentInstance(comp_def, vars, pars, time_bounds)
end

"""
    _build(comp_def::AbstractCompositeComponentDef,
                    var_dict::Dict{ComponentPath, Any},
                    par_dict::Dict{Tuple{ComponentPath, Symbol}, Any},
                    time_bounds::Tuple{Int, Int})

Return a built CompositeComponentInstance created using AbstractCompositeComponentDef
`comp_def`, variables and parameters from `var_dict` and `par_dict` and the time
bounds set by `time_bounds`.
"""
function _build(comp_def::AbstractCompositeComponentDef,
                var_dict::Dict{ComponentPath, Any},
                par_dict::Dict{Tuple{ComponentPath, Symbol}, Any},
                time_bounds::Tuple{Int, Int})
    # @info "_build composite $(comp_def.comp_id)"
    # @info "  var_dict $(var_dict)"
    # @info "  par_dict $(par_dict)"

    comps = [_build(cd, var_dict, par_dict, time_bounds) for cd in compdefs(comp_def)]

    variables = _get_variables(comp_def)
    parameters = _get_parameters(comp_def)

    return CompositeComponentInstance(comps, comp_def, time_bounds, variables, parameters)
end

"""
    _get_variables(comp_def::AbstractCompositeComponentDef)

Return a vector of NamedTuples for all variables in the CompositeComponentInstance
`comp_def`.
"""
# helper functions for to create the variables and parameters NamedTuples for a
# CompositeComponentInstance
function _get_variables(comp_def::AbstractCompositeComponentDef)

    namespace = comp_def.namespace
    var_defs = filter(namespace -> isa(namespace.second, CompositeVariableDef), namespace)
    names = [k for (k,v) in var_defs]
    vals = [v.ref for (k,v) in var_defs]
    variables = (; zip(names, vals)...)

    return variables
end

"""
    _get_parameters(comp_def::AbstractCompositeComponentDef)

Return a vector of NamedTuples for all parameters in the CompositeComponentInstance
`comp_def`.
"""
function _get_parameters(comp_def::AbstractCompositeComponentDef)

    namespace = comp_def.namespace
    par_defs = filter(namespace -> isa(namespace.second, CompositeParameterDef), namespace)
    names = [k for (k,v) in par_defs]
    vals = [v.refs[1] for (k,v) in par_defs]
    parameters = (; zip(names, vals)...)

    return parameters
end

"""
    _build(md::ModelDef)

Build ModelDef `md` (lowest build function called by `build(md::ModelDef)`) and return the ModelInstance..
"""
function _build(md::ModelDef)

    # @info "_build(md)"
    add_connector_comps!(md)

    # check if any of the parameters initialized with the value(s) of nothing
    # are  still nothing
    nothingparams = nothing_params(md)
    if ! isempty(nothingparams)
        params = join([p.datum_name for p in nothingparams], "\n  ")
        error("Cannot build model; the following parameters still have values of `nothing` and need to be updated:\n  $params")
    end

    vdict = _instantiate_vars(md)
    pdict = _collect_params(md, vdict)

    # @info "vdict: $vdict"
    # @info "pdict: $pdict"

    t = dimension(md, :time)
    time_bounds = (firstindex(t), lastindex(t))

    _propagate_time_dim!(md, t) # this might not be needed, but is a final propagation to double check everything

    ci = _build(md, vdict, pdict, time_bounds)
    mi = ModelInstance(ci, md)
    return mi
end

"""
    build(m::Model)

Build Model `m` and return the ModelInstance.
"""
function build(m::Model)
    # Reference a copy in the ModelInstance to avoid changes underfoot
    md = deepcopy(m.md)
    mi = _build(md)
    return mi
end

"""
    build!(m::Model)

Build Model `m` for and dirty `m`'s ModelDef.
"""
function build!(m::Model)
    m.mi = build(m)
    m.md.dirty = false
    return nothing
end

"""
    create_marginal_model(base::Model, delta::Float64=1.0)

Create a `MarginalModel` where `base` is the baseline model and `delta` is the
difference used to create the `marginal` model.  Return the resulting `MarginaModel`
which shares the internal `ModelDef` between the `base` and `marginal`.
"""
function create_marginal_model(base::Model, delta::Float64=1.0)
    # Make sure the base has a ModelInstance before we copy since this
    # copies the ModelDef to avoid being affected by later changes.
    if ! is_built(base)
        build!(base)
    end

    # Create a marginal model, which shares the internal ModelDef between base and marginal
    mm = MarginalModel(base, delta)
end

"""
    Base.run(mm::MarginalModel; ntimesteps::Int=typemax(Int))

Run the marginal model `mm` once with `ntimesteps`.
"""
function Base.run(mm::MarginalModel; ntimesteps::Int=typemax(Int))
    run(mm.base, ntimesteps=ntimesteps)
    run(mm.modified, ntimesteps=ntimesteps)
end

"""
    build!(mm::MarginalModel)

Build MarginalModel `mm` by building both its `base` and `modified models`.
"""
function build!(mm::MarginalModel)
    build!(mm.base)
    build!(mm.modified)
end


text:
using Distributions, Plots, LaTeXStrings; pyplot()

alpha, beta = 8, 2
data = [2,1,0,0,1,0,2,2,5,2,4,0,3,2,5,0]

newAlpha, newBeta = alpha + sum(data), beta + length(data)
post = Gamma(newAlpha, 1/newBeta)

xGrid = quantile(post,0.01):0.001:quantile(post,0.99)
significance = 0.9; halfAlpha = (1-significance)/2

coverage(l,u) = cdf(post,u) - cdf(post,l)

function classicalCI(dist)
    l, u = mode(dist),mode(dist)
    bestl, bestu = l, u
    while  coverage(l,u) < significance
        l -= 0.00001; u += 0.00001
    end
    (l,u)
end
equalTailCI(dist) = (quantile(dist,halfAlpha), quantile(dist,1-halfAlpha))
function highestDensityCI(dist)
    height = 0.999 * maximum(pdf.(dist,xGrid))
    l,u = mode(dist),mode(dist)
    while coverage(l,u) <= significance
        range = filter(theta -> pdf(dist,theta) > height, xGrid)
        l,u = minimum(range), maximum(range)
        height -= 0.00001
    end
    (l,u)
end

l1, u1 = classicalCI(post)
l2, u2 = equalTailCI(post)
l3, u3 = highestDensityCI(post)
println("Classical: ", (l1,u1), "\tWidth: ",u1-l1,
	"\tCoverage: ", coverage(l1,u1))
println("Equal tails: ", (l2,u2), "\tWidth: ",u2-l2,
	"\tCoverage: ", coverage(l2,u2))
println("Highest density: ", (l3,u3), "\tWidth: ",u3-l3,
	"\tCoverage: ", coverage(l3,u3))

plot(xGrid,pdf.(post,xGrid),  yticks=(0:0.25:1.25),
	c=:black, label="Gamma Posterior Distribution",
    	xlims=(1.4, 2.9), ylims=(-0.4,1.25))
plot!([l1,u1],[-0.1,-0.1], label="Classic CI",
    	c=:blue, shape=:vline, ms=16)
plot!([l2,u2],[-0.2,-0.2], label="Equal Tail CI",
    	c=:red, shape=:vline, ms=16)
plot!([l3,u3],[-0.3,-0.3], label="Highest Density CI",
	c=:green, shape=:vline, ms=16, xlabel=L"\lambda", ylabel="Density")

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Plots, Printf, Statistics, LinearAlgebra
Dat = Float64  # Precision (double=Float64 or single=Float32)
# Helper fun
@views     av(A) = 0.25*(A[1:end-1,1:end-1].+A[2:end,1:end-1].+A[1:end-1,2:end].+A[2:end,2:end])
@views  av_xa(A) =  0.5*(A[1:end-1,:].+A[2:end,:])
@views  av_ya(A) =  0.5*(A[:,1:end-1].+A[:,2:end])
@views maxloc(A) = max.(A[1:end-2,1:end-2],A[1:end-2,2:end-1],A[1:end-2,3:end],
                        A[2:end-1,1:end-2],A[2:end-1,2:end-1],A[2:end-1,3:end],
                        A[3:end  ,1:end-2],A[3:end  ,2:end-1],A[3:end  ,3:end])
# 2D Stokes routine
@views function Stokes2D_vep()
    do_DP   = true               # do_DP=false: Von Mises, do_DP=true: Drucker-Prager (friction angle)
    η_reg   = 1.2e-2             # regularisation "viscosity"
    # Physics
    Lx, Ly  = 1.0, 1.0           # domain size
    radi    = 0.01               # inclusion radius
    τ_y     = 1.6                # yield stress. If do_DP=true, τ_y stand for the cohesion: c*cos(ϕ)
    sinϕ    = sind(30)*do_DP     # sinus of the friction angle
    μ0      = 1.0                # viscous viscosity
    G0      = 1.0                # elastic shear modulus
    Gi      = G0/(6.0-4.0*do_DP) # elastic shear modulus perturbation
    εbg     = 1.0                # background strain-rate
    dt      = μ0/G0/4.0 # assumes Maxwell time of 4
    # Numerics
    nt      = 10                 # number of time steps
    nx, ny  = 63, 63             # numerical grid resolution
    # Vdmp    = 4.0                # convergence acceleration (damping)
    # Vsc     = 2.0                # iterative time step limiter
    # Ptsc    = 6.0                # iterative time step limiter
    ε       = 1e-6               # nonlinear tolerence
    iterMax = 3e4                # max number of iters
    nout    = 200                # check frequency
    Re      = 5π                 # Reynolds number
    r       = 1.0                # Bulk to shear elastic modulus ratio
    CFL     = 0.06/sqrt(2)        # CFL number # DEBUG was 0.9
    # Preprocessing
    dx, dy  = Lx/nx, Ly/ny
    max_lxy = max(Lx,Ly)
    Vpdτ    = min(dx,dy)*CFL
    # Array initialisation
    Pt      = zeros(Dat, nx  ,ny  )
    dτ_Rho  = zeros(Dat, nx  ,ny  )
    Gdτ     = zeros(Dat, nx  ,ny  )
    ∇V      = zeros(Dat, nx  ,ny  )
    Vx      = zeros(Dat, nx+1,ny  )
    Vy      = zeros(Dat, nx  ,ny+1)
    Exx     = zeros(Dat, nx  ,ny  )
    Eyy     = zeros(Dat, nx  ,ny  )
    Exyv    = zeros(Dat, nx+1,ny+1)
    Exx1    = zeros(Dat, nx  ,ny  )
    Eyy1    = zeros(Dat, nx  ,ny  )
    Exy1    = zeros(Dat, nx  ,ny  )
    Exyv1   = zeros(Dat, nx+1,ny+1)
    Txx     = zeros(Dat, nx  ,ny  )
    Tyy     = zeros(Dat, nx  ,ny  )
    Txy     = zeros(Dat, nx  ,ny  )
    Txyv    = zeros(Dat, nx+1,ny+1)
    Txx_o   = zeros(Dat, nx  ,ny  )
    Tyy_o   = zeros(Dat, nx  ,ny  )
    Txy_o   = zeros(Dat, nx  ,ny  )
    Txyv_o  = zeros(Dat, nx+1,ny+1)
    Tii     = zeros(Dat, nx  ,ny  )
    Eii     = zeros(Dat, nx  ,ny  )
    F       = zeros(Dat, nx  ,ny  )
    Fchk    = zeros(Dat, nx  ,ny  )
    Pla     = zeros(Dat, nx  ,ny  )
    λ       = zeros(Dat, nx  ,ny  )
    dQdTxx  = zeros(Dat, nx  ,ny  )
    dQdTyy  = zeros(Dat, nx  ,ny  )
    dQdTxy  = zeros(Dat, nx  ,ny  )
    Rx      = zeros(Dat, nx-1,ny  )
    Ry      = zeros(Dat, nx  ,ny-1)
    dVx     = zeros(Dat, nx-1,ny  )
    dVy     = zeros(Dat, nx  ,ny-1)
    Rog     = zeros(Dat, nx  ,ny  )
    η_v     =    μ0*ones(Dat, nx, ny)
    η_e     = dt*G0*ones(Dat, nx, ny)
    η_ev    = dt*G0*ones(Dat, nx+1, ny+1)
    η_ve    =       ones(Dat, nx, ny)
    η_vep   =       ones(Dat, nx, ny)
    η_vepv  =       ones(Dat, nx+1, ny+1)
    # Initial condition
    xc, yc  = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny)
    xc, yc  = LinRange(dx/2, Lx-dx/2, nx), LinRange(dy/2, Ly-dy/2, ny)
    xv, yv  = LinRange(0.0, Lx, nx+1), LinRange(0.0, Ly, ny+1)
    (Xvx,Yvx) = ([x for x=xv,y=yc], [y for x=xv,y=yc])
    (Xvy,Yvy) = ([x for x=xc,y=yv], [y for x=xc,y=yv])
    radc      = (xc.-Lx./2).^2 .+ (yc'.-Ly./2).^2
    radv      = (xv.-Lx./2).^2 .+ (yv'.-Ly./2).^2
    η_e[radc.<radi] .= dt*Gi
    η_ev[radv.<radi].= dt*Gi
    η_ve   .= (1.0./η_e + 1.0./η_v).^-1
    Vx     .=   εbg.*Xvx
    Vy     .= .-εbg.*Yvy
    η_veτ   = copy(η_ve)
    η_vepτ  = copy(η_ve)
    # Time loop
    t=0.0; evo_t=[]; evo_Txx=[]
    for it = 1:nt
        iter=1; err=2*ε; err_evo1=[]; err_evo2=[]
        Txx_o.=Txx; Tyy_o.=Tyy; Txy_o.=av(Txyv); Txyv_o.=Txyv; λ.=0.0
        local itg
        while (err>ε && iter<=iterMax)

            # η_veτ[2:end-1,2:end-1] .= maxloc(η_ve); η_veτ[1,:] .= η_veτ[2,:]; η_veτ[end,:] .= η_veτ[end-1,:]; η_veτ[:,1] .= η_veτ[:,2]; η_veτ[:,end] .= η_veτ[:,end-1]
            # dτ_Rho .= Vpdτ .* max_lxy ./ Re ./ η_veτ  # could be η_veτ here with maxloc...

            η_vepτ[2:end-1,2:end-1] .= maxloc(η_vep); η_vepτ[1,:] .= η_vepτ[2,:]; η_vepτ[end,:] .= η_vepτ[end-1,:]; η_vepτ[:,1] .= η_vepτ[:,2]; η_vepτ[:,end] .= η_vepτ[:,end-1]
            dτ_Rho .= Vpdτ .* max_lxy ./ Re ./ η_vepτ # could be η_veτ here with maxloc...

            Gdτ    .= Vpdτ .^ 2 ./ dτ_Rho ./ (r .+ 2)
            # divergence - pressure
            ∇V     .= diff(Vx, dims=1)./dx .+ diff(Vy, dims=2)./dy
            Pt     .= Pt .- r .* Gdτ .* ∇V
            # strain rates
            Exx    .= diff(Vx, dims=1)./dx .- 1.0/3.0*∇V
            Eyy    .= diff(Vy, dims=2)./dy .- 1.0/3.0*∇V
            Exyv[2:end-1,2:end-1] .= 0.5.*(diff(Vx[2:end-1,:], dims=2)./dy .+ diff(Vy[:,2:end-1], dims=1)./dx)
            # visco-elastic strain rates
            Exx1   .=    Exx   .+ Txx_o ./2.0./η_e
            Eyy1   .=    Eyy   .+ Tyy_o ./2.0./η_e
            Exyv1  .=    Exyv  .+ Txyv_o./2.0./η_ev
            Exy1   .= av(Exyv) .+ Txy_o ./2.0./η_e
            Eii    .= sqrt.(0.5*(Exx1.^2 .+ Eyy1.^2) .+ Exy1.^2)
            # trial stress
            # Txx    .= 2.0.*η_ve.*Exx1
            # Tyy    .= 2.0.*η_ve.*Eyy1
            # Txy    .= 2.0.*η_ve.*Exy1
            Txx    .= (Txx .+ 2.0 .* Gdτ .* Exx1) ./ (Gdτ ./ η_ve .+ 1.0)
            Tyy    .= (Tyy .+ 2.0 .* Gdτ .* Eyy1) ./ (Gdτ ./ η_ve .+ 1.0)
            Txy    .= (Txy .+ 2.0 .* Gdτ .* Exy1) ./ (Gdτ ./ η_ve .+ 1.0)
            Tii    .= sqrt.(0.5*(Txx.^2 .+ Tyy.^2) .+ Txy.^2)
            # yield function
            F      .= Tii .- τ_y .- Pt.*sinϕ
            Pla    .= 0.0
            Pla    .= F .> 0.0
            λ      .= Pla.*F./(η_ve .+ η_reg)
            dQdTxx .= 0.5.*Txx./Tii
            dQdTyy .= 0.5.*Tyy./Tii
            dQdTxy .=      Txy./Tii
            # plastic corrections
            # Txx    .= 2.0.*η_ve.*(Exx1 .-      λ.*dQdTxx)
            # Tyy    .= 2.0.*η_ve.*(Eyy1 .-      λ.*dQdTyy)
            # Txy    .= 2.0.*η_ve.*(Exy1 .- 0.5.*λ.*dQdTxy)
            Txx    .= (Txx .+ 2.0 .* Gdτ .* (Exx1 .-        λ .* dQdTxx)) ./ (Gdτ ./ η_ve .+ 1.0)
            Tyy    .= (Tyy .+ 2.0 .* Gdτ .* (Eyy1 .-        λ .* dQdTyy)) ./ (Gdτ ./ η_ve .+ 1.0)
            Txy    .= (Txy .+ 2.0 .* Gdτ .* (Exy1 .- 0.5 .* λ .* dQdTxy)) ./ (Gdτ ./ η_ve .+ 1.0)
            Tii    .= sqrt.(0.5*(Txx.^2 .+ Tyy.^2) .+ Txy.^2)
            Fchk   .= Tii .- τ_y .- Pt.*sinϕ .- λ.*η_reg
            η_vep  .= Tii ./ 2.0 ./ Eii
            η_vepv[2:end-1,2:end-1] .= av(η_vep); η_vepv[1,:].=η_vepv[2,:]; η_vepv[end,:].=η_vepv[end-1,:]; η_vepv[:,1].=η_vepv[:,2]; η_vepv[:,end].=η_vepv[:,end-1]
            Txyv   .= 2.0 .* η_vepv .* Exyv1
            # velocities
            dVx    .= av_xa(dτ_Rho) .* (.-diff(Pt, dims=1)./dx .+ diff(Txx, dims=1)./dx .+ diff(Txyv[2:end-1,:], dims=2)./dy)
            dVy    .= av_ya(dτ_Rho) .* (.-diff(Pt, dims=2)./dy .+ diff(Tyy, dims=2)./dy .+ diff(Txyv[:,2:end-1], dims=1)./dx .+ av_ya(Rog))
            Vx[2:end-1,:] .= Vx[2:end-1,:] .+ dVx
            Vy[:,2:end-1] .= Vy[:,2:end-1] .+ dVy
            # convergence check
            if mod(iter, nout)==0
                Rx .= .-diff(Pt, dims=1)./dx .+ diff(Txx, dims=1)./dx .+ diff(Txyv[2:end-1,:], dims=2)./dy
                Ry .= .-diff(Pt, dims=2)./dy .+ diff(Tyy, dims=2)./dy .+ diff(Txyv[:,2:end-1], dims=1)./dx .+ av_ya(Rog)
                norm_Rx = norm(Rx)/length(Rx); norm_Ry = norm(Ry)/length(Ry); norm_∇V = norm(∇V)/length(∇V)
                err = maximum([norm_Rx, norm_Ry, norm_∇V])
                push!(err_evo1, err); push!(err_evo2, itg)
                @printf("it = %d, iter = %d, err = %1.2e norm[Rx=%1.2e, Ry=%1.2e, ∇V=%1.2e] (Fchk=%1.2e) \n", it, itg, err, norm_Rx, norm_Ry, norm_∇V, maximum(Fchk))
            end
            iter+=1; itg=iter
        end
        t = t + dt
        push!(evo_t, t); push!(evo_Txx, maximum(Txx))
        # Plotting
        p1 = heatmap(xv, yc, Vx' , aspect_ratio=1, xlims=(0, Lx), ylims=(dy/2, Ly-dy/2), c=:inferno, title="Vx")
        # p2 = heatmap(xc, yv, Vy' , aspect_ratio=1, xlims=(dx/2, Lx-dx/2), ylims=(0, Ly), c=:inferno, title="Vy")
        p2 = heatmap(xc, yc, η_vep' , aspect_ratio=1, xlims=(dx/2, Lx-dx/2), ylims=(0, Ly), c=:inferno, title="η_vep")
        p3 = heatmap(xc, yc, Tii' , aspect_ratio=1, xlims=(dx/2, Lx-dx/2), ylims=(0, Ly), c=:inferno, title="τii")
        p4 = plot(evo_t, evo_Txx , legend=false, xlabel="time", ylabel="max(τxx)", linewidth=0, markershape=:circle, framestyle=:box, markersize=3)
            plot!(evo_t, 2.0.*εbg.*μ0.*(1.0.-exp.(.-evo_t.*G0./μ0)), linewidth=2.0) # analytical solution for VE loading
            plot!(evo_t, 2.0.*εbg.*μ0.*ones(size(evo_t)), linewidth=2.0)            # viscous flow stress
            if !do_DP plot!(evo_t, τ_y*ones(size(evo_t)), linewidth=2.0) end        # von Mises yield stress
        display(plot(p1, p2, p3, p4))
    end
    return
end

Stokes2D_vep()


text:
using Krotov:
    SquareParametrization,
    TanhParametrization,
    TanhSqParametrization,
    LogisticParametrization,
    LogisticSqParametrization

using Plots
Plots.default(
    linewidth               = 3,
    size                    = (550, 300),
    legend                  = :top,
    foreground_color_legend = nothing,
    background_color_legend = RGBA(1, 1, 1, 0.8),
)

function plot_symmetric_parametrization_comparison()

    u_vals = collect(range(-3, 3, length=101))
    ϵ_vals = collect(range(-1, 1, length=101))

    ϵ_min = -1.0
    ϵ_max = 1.0

    pnl1 = plot(
        u_vals,
        u_vals;
        linestyle=:dash,
        color="black",
        label="",
        xlabel="u",
        ylabel="ϵ",
        legend=false
    )
    plot!(
        pnl1,
        u_vals,
        TanhParametrization(ϵ_min, ϵ_max).epsilon_of_u.(u_vals),
        label="Tanh"
    )
    plot!(
        pnl1,
        u_vals,
        LogisticParametrization(ϵ_min, ϵ_max).epsilon_of_u.(u_vals),
        label="Logistic(k=1)"
    )
    plot!(
        pnl1,
        u_vals,
        LogisticParametrization(ϵ_min, ϵ_max, k=4).epsilon_of_u.(u_vals),
        label="Logistic(k=4)"
    )
    ylims!(pnl1, (-1.2, 1.2))

    pnl2 = plot(
        ϵ_vals,
        ϵ_vals;
        linestyle=:dash,
        color="black",
        label="",
        xlabel="ϵ",
        ylabel="u"
    )
    plot!(
        pnl2,
        ϵ_vals,
        TanhParametrization(ϵ_min, ϵ_max).u_of_epsilon.(ϵ_vals),
        label="Tanh"
    )
    plot!(
        pnl2,
        ϵ_vals,
        LogisticParametrization(ϵ_min, ϵ_max).u_of_epsilon.(ϵ_vals),
        label="Logistic(k=1)"
    )
    plot!(
        pnl2,
        ϵ_vals,
        LogisticParametrization(ϵ_min, ϵ_max, k=4).u_of_epsilon.(ϵ_vals),
        label="Logistic(k=4)"
    )
    ylims!(pnl2, (-3, 3))

    pnl3 = plot(
        u_vals,
        [1.0 for _ in u_vals];
        linestyle=:dash,
        color="black",
        label="",
        xlabel="u",
        ylabel="∂ϵ/∂u",
        legend=false
    )
    plot!(
        pnl3,
        u_vals,
        TanhParametrization(ϵ_min, ϵ_max).de_du_derivative.(u_vals),
        label="Tanh"
    )
    plot!(
        pnl3,
        u_vals,
        LogisticParametrization(ϵ_min, ϵ_max).de_du_derivative.(u_vals),
        label="Logistic(k=1)"
    )
    plot!(
        pnl3,
        u_vals,
        LogisticParametrization(ϵ_min, ϵ_max, k=4).de_du_derivative.(u_vals),
        label="Logistic(k=4)"
    )
    ylims!(pnl3, (0, 2))

    plot(
        pnl1,
        pnl2,
        pnl3,
        layout=(1, 3),
        size=(1000, 300),
        left_margin=20Plots.px,
        bottom_margin=20Plots.px
    )

end

if abspath(PROGRAM_FILE) == @__FILE__
    gui(plot_symmetric_parametrization_comparison())
    readline()
end


text:
function dataset_preprocess(dataset::Type{Planetoid})
    return function preprocess(local_path)
        for subds in subdatasets(dataset)
            graph_file = @datadep_str "Planetoid/ind.$(subds).graph"
            trainX_file = @datadep_str "Planetoid/ind.$(subds).x"
            trainy_file = @datadep_str "Planetoid/ind.$(subds).y"
            testX_file = @datadep_str "Planetoid/ind.$(subds).tx"
            testy_file = @datadep_str "Planetoid/ind.$(subds).ty"
            allX_file = @datadep_str "Planetoid/ind.$(subds).allx"
            ally_file = @datadep_str "Planetoid/ind.$(subds).ally"

            train_X = read_data(trainX_file)
            train_y = read_data(trainy_file)
            test_X = read_data(testX_file)
            test_y = read_data(testy_file)
            all_X = read_data(allX_file)
            all_y = read_data(ally_file)
            graph = read_graph(graph_file)

            num_V = length(graph)
            sg = to_simplegraph(graph, num_V)
            num_E = ne(sg)
            feat_dim = size(all_X, 2)
            label_dim = size(all_y, 2)
            train_X, train_y = sparse(train_X'), sparse(train_y')
            test_X, test_y = sparse(test_X'), sparse(test_y')
            all_X, all_y = sparse(all_X'), sparse(all_y')
            raw = Dict(:graph=>graph, :train_X=>train_X, :train_y=>train_y,
                       :test_X=>test_X, :test_y=>test_y, :all_X=>all_X, :all_y=>all_y)
            meta = (graph=(num_V=num_V, num_E=num_E),
                    train=(features_dim=(feat_dim, size(train_X, 2)), labels_dim=(label_dim, size(train_y, 2))),
                    test=(features_dim=(feat_dim, size(test_X, 2)), labels_dim=(label_dim, size(test_y, 2))),
                    all=(features_dim=(feat_dim, size(all_X, 2)), labels_dim=(label_dim, size(all_y, 2)))
                    )

            graphfile = replace(graph_file, "ind.$(subds).graph"=>"$(subds).graph.jld2")
            trainfile = replace(graph_file, "ind.$(subds).graph"=>"$(subds).train.jld2")
            testfile = replace(graph_file, "ind.$(subds).graph"=>"$(subds).test.jld2")
            allfile = replace(graph_file, "ind.$(subds).graph"=>"$(subds).all.jld2")
            rawfile = replace(graph_file, "ind.$(subds).graph"=>"$(subds).raw.jld2")
            metadatafile = replace(graph_file, "ind.$(subds).graph"=>"$(subds).metadata.jld2")

            @save graphfile sg
            @save trainfile train_X train_y
            @save testfile test_X test_y
            @save allfile all_X all_y
            @save rawfile raw
            @save metadatafile meta
        end
    end
end

function read_data(filename)
    py"""
    import pickle
    from scipy.sparse import csr_matrix

    with open($filename,"rb") as f:
        u = pickle._Unpickler(f)
        u.encoding = "latin1"
        data = u.load()

    if type(data) is csr_matrix:
        data = data.toarray()
    """
    return SparseMatrixCSC(Array(py"data"))
end

read_index(filename) = map(x -> parse(Int64, x), readlines(filename))

function read_graph(filename)
    py"""
    import pickle

    with open($filename,"rb") as f:
        u = pickle._Unpickler(f)
        u.encoding = "latin1"
        data = u.load()
    """
    return Dict(py"data")
end

## Cora dataset

function dataset_preprocess(dataset::Type{Cora})
    return function preprocess(local_path)
        py"""
        import numpy as np
        import scipy.sparse as sp
        data = np.load($local_path, allow_pickle=True)
        A = sp.csr_matrix((data['adj_data'], data['adj_indices'], data['adj_indptr']), shape=data['adj_shape'])
        X = sp.csr_matrix((data['attr_data'], data['attr_indices'], data['attr_indptr']), shape=data['attr_shape'])
        """

        graph = SparseMatrixCSC(Array(py"A.toarray()"))
        X = SparseMatrixCSC(Array(py"X.toarray()"))
        y = py"data['labels']"
        raw = Dict(:graph=>graph, :all_X=>X, :all_y=>y)

        all_X = sparse(X')
        all_y = Matrix{UInt16}(y')
        sg = to_simpledigraph(graph)
        meta = (graph=(num_V=nv(sg), num_E=ne(sg)),
                all=(features_dim=size(all_X), labels_dim=size(all_y))
                )

        graphfile = replace(local_path, "cora.npz"=>"cora.graph.jld2")
        # trainfile = replace(local_path, "cora.npz"=>"cora.train.jld2")
        # testfile = replace(local_path, "cora.npz"=>"cora.test.jld2")
        allfile = replace(local_path, "cora.npz"=>"cora.all.jld2")
        rawfile = replace(local_path, "cora.npz"=>"cora.raw.jld2")
        metadatafile = replace(local_path, "cora.npz"=>"cora.metadata.jld2")

        @save graphfile sg
        # @save trainfile train_X train_y
        # @save testfile test_X test_y
        @save allfile all_X all_y
        @save rawfile raw
        @save metadatafile meta
    end
end

## PPI dataset

function dataset_preprocess(dataset::Type{PPI})
    return function preprocess(local_path)
        unzip_zip(local_path)

        for phase in ["train", "test", "valid"]
            graph_file = @datadep_str "PPI/$(phase)_graph.json"
            id_file = @datadep_str "PPI/$(phase)_graph_id.npy"
            X_file = @datadep_str "PPI/$(phase)_feats.npy"
            y_file = @datadep_str "PPI/$(phase)_labels.npy"

            py"""
            import numpy as np
            ids = np.load($id_file)
            X = np.load($X_file)
            y = np.load($y_file)
            """

            X = Matrix{Float32}(py"X")
            y = SparseMatrixCSC{Int32,Int64}(Array(py"y"))
            ids = Array(py"ids")
            graph = read_ppi_graph(graph_file)

            jld2file = replace(local_path, "ppi.zip"=>"ppi.$(phase).jld2")
            @save jld2file graph X y ids
        end
    end
end

function read_ppi_graph(filename::String)
    d = JSON.Parser.parsefile(filename)
    g = SimpleDiGraph{Int32}(length(d["nodes"]))

    for pair in d["links"]
        add_edge!(g, pair["source"], pair["target"])
    end
    g
end

## Reddit dataset

function dataset_preprocess(dataset::Type{Reddit})
    return function preprocess(local_path)
        unzip_zip(local_path)

        graph_file = datadep"Reddit/reddit_graph.npz"
        data_file = datadep"Reddit/reddit_data.npz"

        graphfile = replace(local_path, "reddit.zip"=>"reddit.graph.jld2")
        allfile = replace(local_path, "reddit.zip"=>"reddit.all.jld2")
        rawfile = replace(local_path, "reddit.zip"=>"reddit.raw.jld2")
        metadatafile = replace(local_path, "reddit.zip"=>"reddit.metadata.jld2")

        graph, X, y = to_reddit_rawfile(graph_file, data_file, rawfile)

        sg = to_simplegraph(graph)
        all_X = Matrix(X')
        all_y = Matrix{UInt16}(y')
        meta = (graph=(num_V=nv(sg), num_E=ne(sg)),
                all=(features_dim=size(all_X), labels_dim=size(all_y))
                )

        @save graphfile sg
        @save allfile all_X all_y
        @save metadatafile meta
    end
end

function to_reddit_rawfile(graph_file, data_file, rawfile)
    py"""
    import numpy as np
    import scipy.sparse as sp
    graph = np.load($graph_file, allow_pickle=True)
    data = np.load($data_file, allow_pickle=True)
    """

    graph = sparse(Vector(py"graph['row']") .+ 1,
                   Vector(py"graph['col']") .+ 1,
                   Vector{Int32}(py"graph['data']"))
    X = Matrix{Float32}(py"data['feature']")
    y = Vector{Int32}(py"data['label']")
    ids = Vector{Int32}(py"data['node_ids']")
    types = Vector{UInt8}(py"data['node_types']")
    raw = Dict(:graph=>graph, :X=>X, :y=>y, :ids=>ids, :types=>types)

    @save rawfile raw

    graph, X, y
end

## Reddit dataset

function dataset_preprocess(dataset::Type{QM7b})
    return function preprocess(local_path)
        vars = matread(local_path)
        names = vars["names"]
        X = vars["X"]
        T = Matrix{Float32}(vars["T"])
        raw = Dict(:names=>names, :X=>X, :T=>T)

        rawfile = replace(local_path, "qm7b.mat"=>"qm7b.raw.jld2")
        @save rawfile raw
    end
end

## Entities dataset

function dataset_preprocess(dataset::Type{Entities})
    return function preprocess(local_path)
        for subds in subdatasets(dataset)
            tgz_file = @datadep_str "Entities/$(subds).tgz"
            dir = joinpath(dirname(tgz_file), "$(subds)")
            unzip_tgz(tgz_file, dir)

            # nt_file = joinpath(dir, "$(dataset)_stripped.nt.gz")
            # train_file = joinpath(dir, "trainingSet.tsv")
            # test_file = joinpath(dir, "testSet.tsv")

            # py"""
            # import rdflib as rdf
            # import gzip
            # g = rdf.Graph()
            # with gzip.open($nt_file, 'rb') as file:
            #     g.parse(file, format='nt')
            # """

            # train = CSV.read(train_file, delim='\t')
            # test = CSV.read(test_file, delim='\t')
            # uri = HTTP.URI("http://data.bgs.ac.uk/id/Geochronology/Division/CN")
        end
    end
end


text:
"""
    ro = IdOffsetRange(r::AbstractUnitRange, offset=0)

Construct an "identity offset range". Numerically, `collect(ro) == collect(r) .+ offset`,
with the additional property that `axes(ro, 1) = axes(r, 1) .+ offset`.
When `r` starts at 1, then `ro[i] == i` and even `ro[ro] == ro`,
i.e., it's the "identity," which is the origin of the "Id" in `IdOffsetRange`.

# Examples

The most common case is shifting a range that starts at 1 (either `1:n` or `Base.OneTo(n)`):
```jldoctest ior
julia> using OffsetArrays: IdOffsetRange

julia> ro = IdOffsetRange(1:3, -2)
IdOffsetRange(values=-1:1, indices=-1:1)

julia> axes(ro, 1)
IdOffsetRange(values=-1:1, indices=-1:1)

julia> ro[-1]
-1

julia> ro[3]
ERROR: BoundsError: attempt to access 3-element $(IdOffsetRange{Int,UnitRange{Int}}) with indices -1:1 at index [3]
```

If the range doesn't start at 1, the values may be different from the indices:
```jldoctest ior
julia> ro = IdOffsetRange(11:13, -2)
IdOffsetRange(values=9:11, indices=-1:1)

julia> axes(ro, 1)     # 11:13 is indexed by 1:3, and the offset is also applied to the axes
IdOffsetRange(values=-1:1, indices=-1:1)

julia> ro[-1]
9

julia> ro[3]
ERROR: BoundsError: attempt to access 3-element $(IdOffsetRange{Int,UnitRange{Int}}) with indices -1:1 at index [3]
```

# Extended help

Construction/coercion preserves the (shifted) values of the input range, but may modify
the indices if required by the specified types. For example,

    r = OffsetArrays.IdOffsetRange{Int,UnitRange{Int}}(3:4)

has `r[1] == 3` and `r[2] == 4`, whereas

    r = OffsetArrays.IdOffsetRange{Int,Base.OneTo{Int}}(3:4)

has `r[3] == 3` and `r[4] == 4`, and `r[1]` would throw a `BoundsError`.
In this latter case, a shift in the axes was needed because `Base.OneTo` ranges
must start with value 1.

!!! warning

    In the future, *conversion* will preserve both the values and
    the indices, throwing an error when this is not achievable. For instance,

        r = convert(OffsetArrays.IdOffsetRange{Int,UnitRange{Int}}, 3:4)

    has `r[1] == 3` and `r[2] == 4` and would satisfy `r == 3:4`, whereas

    ```julia
    julia> convert(OffsetArrays.IdOffsetRange{Int,Base.OneTo{Int}}, 3:4)    # future behavior, not present behavior
    ERROR: ArgumentError: first element must be 1, got 3
    ```

    where the error will arise because the result could not have the same axes as the input.

    An important corollary is that `typeof(r1)(r2)` and `oftype(r1, r2)` will behave differently:
    the first coerces `r2` to be of the type of `r1`, whereas the second converts.
    Developers are urged to future-proof their code by choosing the behavior appropriate for each usage.
"""
struct IdOffsetRange{T<:Integer,I<:AbstractUnitRange{T}} <: AbstractUnitRange{T}
    parent::I
    offset::T

    function IdOffsetRange{T,I}(r::I, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}}
        _bool_check(T, r, offset)
        new{T,I}(r, offset)
    end

    #= This method is necessary to avoid a StackOverflowError in IdOffsetRange{T,I}(r::IdOffsetRange, offset::Integer).
    The type signature in that method is more specific than IdOffsetRange{T,I}(r::I, offset::T),
    so it ends up calling itself if I <: IdOffsetRange.
    =#
    function IdOffsetRange{T,IdOffsetRange{T,I}}(r::IdOffsetRange{T,I}, offset::T) where {T<:Integer,I<:AbstractUnitRange{T}}
        _bool_check(T, r, offset)
        new{T,IdOffsetRange{T,I}}(r, offset)
    end
end

function _bool_check(::Type{Bool}, r, offset)
    # disallow the construction of IdOffsetRange{Bool, UnitRange{Bool}}(true:true, true)
    if offset && (first(r) || last(r))
        throw(ArgumentError("values = $r and offset = $offset can not produce a boolean range"))
    end
    return nothing
end
_bool_check(::Type, r, offset) = nothing

# Construction/coercion from arbitrary AbstractUnitRanges
function IdOffsetRange{T,I}(r::AbstractUnitRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, o = offset_coerce(I, r)
    return IdOffsetRange{T,I}(rc, convert(T, o+offset)::T)
end
function IdOffsetRange{T}(r::AbstractUnitRange, offset::Integer = 0) where T<:Integer
    rc = convert(AbstractUnitRange{T}, r)::AbstractUnitRange{T}
    return IdOffsetRange{T,typeof(rc)}(rc, convert(T, offset)::T)
end
IdOffsetRange(r::AbstractUnitRange{T}, offset::Integer = 0) where T<:Integer =
    IdOffsetRange{T,typeof(r)}(r, convert(T, offset)::T)

# Coercion from other IdOffsetRanges
IdOffsetRange{T,I}(r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
function IdOffsetRange{T,I}(r::IdOffsetRange, offset::Integer = 0) where {T<:Integer,I<:AbstractUnitRange{T}}
    rc, offset_rc = offset_coerce(I, r.parent)
    return IdOffsetRange{T,I}(rc, convert(T, r.offset + offset + offset_rc)::T)
end
IdOffsetRange{T}(r::IdOffsetRange{T}) where {T<:Integer} = r
function IdOffsetRange{T}(r::IdOffsetRange, offset::Integer = 0) where T<:Integer
    return IdOffsetRange{T}(r.parent, r.offset + offset)
end
IdOffsetRange(r::IdOffsetRange) = r

# Constructor to make `show` round-trippable
function IdOffsetRange(; values::AbstractUnitRange{<:Integer}, indices::AbstractUnitRange{<:Integer})
    length(values) == length(indices) || throw(ArgumentError("values and indices must have the same length"))
    offset = first(indices) - 1
    return IdOffsetRange(values .- offset, offset)
end

# Conversions to an AbstractUnitRange{Int} (and to an OrdinalRange{Int,Int} on Julia v"1.6") are necessary
# to evaluate CartesianIndices for BigInt ranges, as their axes are also BigInt ranges
Base.AbstractUnitRange{T}(r::IdOffsetRange) where {T<:Integer} = IdOffsetRange{T}(r)

# A version upper bound on this may be set after https://github.com/JuliaLang/julia/pull/40038 is merged
if v"1.6" <= VERSION
    Base.OrdinalRange{T,T}(r::IdOffsetRange) where {T<:Integer} = IdOffsetRange{T}(r)
end

# TODO: uncomment these when Julia is ready
# # Conversion preserves both the values and the indices, throwing an InexactError if this
# # is not possible.
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange{T,I}) where {T<:Integer,I<:AbstractUnitRange{T}} = r
# Base.convert(::Type{IdOffsetRange{T,I}}, r::IdOffsetRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r.parent), r.offset)
# Base.convert(::Type{IdOffsetRange{T,I}}, r::AbstractUnitRange) where {T<:Integer,I<:AbstractUnitRange{T}} =
#     IdOffsetRange{T,I}(convert(I, r), 0)

offset_coerce(::Type{Base.OneTo{T}}, r::Base.OneTo) where T<:Integer = convert(Base.OneTo{T}, r), 0
function offset_coerce(::Type{Base.OneTo{T}}, r::AbstractUnitRange) where T<:Integer
    o = first(r) - 1
    return Base.OneTo{T}(last(r) - o), o
end
# function offset_coerce(::Type{Base.OneTo{T}}, r::IdOffsetRange) where T<:Integer
#     rc, o = offset_coerce(Base.OneTo{T}, r.parent)

# Fallback, specialze this method if `convert(I, r)` doesn't do what you need
offset_coerce(::Type{I}, r::AbstractUnitRange) where I<:AbstractUnitRange =
    convert(I, r)::I, 0

@inline Base.parent(r::IdOffsetRange) = r.parent
@inline Base.axes(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.axes1(r::IdOffsetRange) = IdOffsetRange(Base.axes1(r.parent), r.offset)
@inline Base.unsafe_indices(r::IdOffsetRange) = (Base.axes1(r),)
@inline Base.length(r::IdOffsetRange) = length(r.parent)
@inline Base.isempty(r::IdOffsetRange) = isempty(r.parent)
Base.reduced_index(i::IdOffsetRange) = typeof(i)(first(i):first(i))
# Workaround for #92 on Julia < 1.4
Base.reduced_index(i::IdentityUnitRange{<:IdOffsetRange}) = typeof(i)(first(i):first(i))
for f in [:firstindex, :lastindex]
    @eval @inline Base.$f(r::IdOffsetRange) = $f(r.parent) + r.offset
end
for f in [:first, :last]
    # coerce the type to deal with values that get promoted on addition (eg. Bool)
    @eval @inline Base.$f(r::IdOffsetRange) = eltype(r)($f(r.parent) + r.offset)
end

# Iteration for an IdOffsetRange
@inline Base.iterate(r::IdOffsetRange, i...) = _iterate(r, i...)
# In general we iterate over the parent term by term and add the offset.
# This might have some performance degradation when coupled with bounds-checking
# See https://github.com/JuliaArrays/OffsetArrays.jl/issues/214
@inline function _iterate(r::IdOffsetRange, i...)
    ret = iterate(r.parent, i...)
    ret === nothing && return nothing
    return (eltype(r)(ret[1] + r.offset), ret[2])
end
# Base.OneTo(n) is known to be exactly equivalent to the range 1:n,
# and has no specialized iteration defined for it,
# so we may add the offset to the range directly and iterate over the result
# This gets around the performance issue described in issue #214
# We use the helper function _addoffset to evaluate the range instead of broadcasting
# just in case this makes it easy for the compiler.
@inline _iterate(r::IdOffsetRange{<:Integer, <:Base.OneTo}, i...) = iterate(_addoffset(r.parent, r.offset), i...)

@inline function Base.getindex(r::IdOffsetRange, i::Integer)
    i isa Bool && throw(ArgumentError("invalid index: $i of type Bool"))
    @boundscheck checkbounds(r, i)
    @inbounds eltype(r)(r.parent[i - r.offset] + r.offset)
end

# Logical indexing following https://github.com/JuliaLang/julia/pull/31829
#= Helper function to perform logical indxeing for boolean ranges
The code implemented is a branch-free version of the following:

    range(first(s) ? first(r) : last(r), length=Int(last(s)))

See https://github.com/JuliaArrays/OffsetArrays.jl/pull/224#discussion_r595635143

Logical indexing does not preserve indices, unlike other forms of vector indexing
=#
@inline function _getindex(r, s::AbstractUnitRange{Bool})
    range(first(r) * first(s) + last(r) * !first(s), length=Int(last(s)))
end
@inline function _getindex(r, s::StepRange{Bool})
    range(first(r) * first(s) + last(r) * !first(s), step = oneunit(step(s)), length=Int(last(s)))
end
@inline function _getindex(r, s::AbstractUnitRange)
    @inbounds rs = r.parent[_subtractoffset(s, r.offset)] .+ r.offset
    _indexedby(rs, axes(s))
end
@inline function _getindex(r, s::StepRange)
    rs = @inbounds r.parent[s .- r.offset] .+ r.offset
    _indexedby(rs, axes(s))
end

for T in [:AbstractUnitRange, :StepRange]
    @eval @inline function Base.getindex(r::IdOffsetRange, s::$T{<:Integer})
        @boundscheck checkbounds(r, s)
        return _getindex(r, s)
    end
end

# These methods are necessary to avoid ambiguity
for R in [:IIUR, :IdOffsetRange]
    @eval @inline function Base.getindex(r::IdOffsetRange, s::$R)
        @boundscheck checkbounds(r, s)
        return _getindex(r, s)
    end
end

# offset-preserve broadcasting
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(-), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .- x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), r::IdOffsetRange{T}, x::Integer) where T =
    IdOffsetRange{T}(r.parent .+ x, r.offset)
Broadcast.broadcasted(::Base.Broadcast.DefaultArrayStyle{1}, ::typeof(+), x::Integer, r::IdOffsetRange{T}) where T =
    IdOffsetRange{T}(x .+ r.parent, r.offset)

Base.show(io::IO, r::IdOffsetRange) = print(io, IdOffsetRange, "(values=",first(r), ':', last(r),", indices=",first(eachindex(r)),':',last(eachindex(r)), ")")

# Optimizations
@inline Base.checkindex(::Type{Bool}, inds::IdOffsetRange, i::Real) = Base.checkindex(Bool, inds.parent, i - inds.offset)

if VERSION < v"1.5.2"
    # issue 100, 133: IdOffsetRange as another index-preserving case shouldn't comtribute offsets
    # fixed by https://github.com/JuliaLang/julia/pull/37204
    @inline Base.compute_offset1(parent, stride1::Integer, dims::Tuple{Int}, inds::Tuple{IdOffsetRange}, I::Tuple) =
        Base.compute_linindex(parent, I) - stride1*first(Base.axes1(inds[1]))
end

# This was deemed "too private" to extend: see issue #184
# # Fixes an inference failure in Base.mapfirst!
# # Test: A = OffsetArray(rand(4,4), (-3,5)); R = similar(A, (1:1, 6:9)); maximum!(R, A)
# if isdefined(Base, :_firstslice)
#     Base._firstslice(i::IdOffsetRange) = IdOffsetRange(Base._firstslice(i.parent), i.offset)
# end


text:
using Distributed
@everywhere import PWS

system_fn = () -> PWS.gene_expression_system(dtimes=0:0.1:10)
smc = PWS.SMCEstimate(128)
result = PWS.run_parallel(system_fn, smc, 32)
println(result)