text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

@test length(ExponentialBackOff(n=10)) == 10
@test collect(ExponentialBackOff(n=10, first_delay=0.01))[1] == 0.01
@test maximum(ExponentialBackOff(n=10, max_delay=0.06)) == 0.06
ratio(x) = x[2:end]./x[1:end-1]
@test all(x->x ≈ 10.0, ratio(collect(ExponentialBackOff(n=10, max_delay=Inf, factor=10, jitter=0.0))))
srand(12345)
@test (mean(ratio(collect(ExponentialBackOff(n=100, max_delay=Inf, factor=1, jitter=0.1)))) - 1.0) < 1e-4

let
    function foo_error(c, n)
        c[1] += 1
        if c[1] <= n
            error("foo")
        end
        return 7
    end

    # Success on first attempt
    c = [0]
    @test retry(foo_error)(c,0) == 7
    @test c[1] == 1

    # Success on second attempt
    c = [0]
    @test retry(foo_error)(c,1) == 7
    @test c[1] == 2

    # 2 failed retry attempts, so exception is raised
    c = [0]
    ex = try retry(foo_error, delays=ExponentialBackOff(n=2))(c,3) catch e; e end
    @test ex.msg == "foo"
    @test c[1] == 3

    c = [0]
    ex = try retry(foo_error, check=(s,e)->(s,isa(e, ErrorException)))(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 2

    c = [0]
    ex = try retry(foo_error, check=(s,e)->(s,e.msg == "foo"))(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 2

    # No retry if condition does not match
    c = [0]
    ex = try retry(foo_error, check=(s,e)->(s,e.msg == "bar"))(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 1

    c = [0]
    ex = try retry(foo_error, check=(s,e)->(s,try e.http_status_code == "503" end != true))(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 2

    c = [0]
    ex = try retry(foo_error, check=(s,e)->(s,isa(e,SystemError)))(c,2) catch e; e end
    @test typeof(ex) == ErrorException
    @test ex.msg == "foo"
    @test c[1] == 1
end


text:
module Spring

using Plots
toto(x) = cos(x)

struct MyStruct
    myfield::Float32
end

function go()
    x=[Float64(i) for i ∈ 1:40]
    y=toto.(x)
    plot(x,y)

    mf=MyStruct(3)
    @show mf

end

end # module

ext: jl


text:
# Test the functions that move between constrained and unconstrained
# parameterizations.

using Celeste: Transform, SensitiveFloats
using DerivativeTestUtils

"""
Generate parameters within the given bounds.
"""
function generate_valid_parameters(
    NumType::DataType, bounds_vec::Vector{ParamBounds})

    @assert NumType <: Number
    S = length(bounds_vec)
    vp = convert(VariationalParams{NumType},
                                 [ zeros(NumType, length(ids)) for s = 1:S ])
        for s=1:S, (param, constraint_vec) in bounds_vec[s]
        is_box = isa(constraint_vec, Array{ParamBox})
        if is_box
            # Box parameters.
            for ind in 1:length(getfield(ids, param))
                constraint = constraint_vec[ind]
                constraint.ub == Inf ?
                    vp[s][getfield(ids, param)[ind]] = constraint.lb + 1.0:
                    vp[s][getfield(ids, param)[ind]] =
                        0.5 * (constraint.ub - constraint.lb) +
                        constraint.lb
            end
        else
            # Simplex parameters can ignore the bounds.
            param_size = size(getfield(ids, param))
            if length(param_size) == 2
                # matrix simplex
                for col in 1:param_size[2]
                    vp[s][getfield(ids, param)[:, col]] = 1 / param_size[1]
                end
            else
                # vector simplex
                vp[s][getfield(ids, param)] = 1 / length(getfield(ids, param))
            end
        end
        end

    vp
end

function test_parameter_conversion()
    blob, ea, body = gen_three_body_dataset();

    transform = get_mp_transform(ea.vp, ea.active_sources, loc_width=1.0);

    function check_transform(transform::DataTransform, ea::ElboArgs)
        original_vp = deepcopy(ea.vp);
        ea_check = deepcopy(ea);

        # Check that the constrain and unconstrain operations undo each other.
        vp_free = Transform.from_vp(transform, ea.vp)
        Transform.to_vp!(transform, vp_free, ea_check.vp)

        for id in fieldnames(ids), s in 1:ea.S
            @test_approx_eq_eps(original_vp[s][getfield(ids, id)],
                                ea_check.vp[s][getfield(ids, id)], 1e-6)
        end

        # Check conversion to and from a vector.
        omitted_ids = Array(Int, 0)
        kept_ids = setdiff(1:length(UnconstrainedParams), omitted_ids)
        vp = deepcopy(ea.vp)
        x = Transform.vp_to_array(transform, vp, omitted_ids)
        @test length(x) == length(vp_free[1]) * length(ea.active_sources)

        vp2 = generate_valid_parameters(Float64, transform.bounds)
        Transform.array_to_vp!(transform, x, vp2, kept_ids)

        for id in fieldnames(ids), si in eachindex(transform.active_sources)
            s = transform.active_sources[si]
            @test_approx_eq_eps(original_vp[s][getfield(ids, id)], vp2[si][getfield(ids, id)], 1e-6)
        end

    end

    transform = get_mp_transform(ea.vp, ea.active_sources, loc_width=1.0);
    check_transform(transform, ea)

    # Test transforming only active sources.
    ea1 = deepcopy(ea);
    ea1.active_sources = [1]
    transform1 = Transform.get_mp_transform(ea1.vp, ea1.active_sources)

    @assert transform1.S == ea.S
    @assert length(transform1.active_sources) == 1
    check_transform(transform1, ea1)

end

function test_transform_simplex_functions()
    function simplex_and_unsimplex{NumType <: Number}(
            param::Vector{NumType}, simplex_box::Transform.SimplexBox)

        param_free = Transform.unsimplexify_parameter(param, simplex_box)
        new_param = Transform.simplexify_parameter(param_free, simplex_box)
        @test_approx_eq param new_param
    end

    for this_scale = [ 1.0, 2.0 ], lb = [0.1, 0.0 ]
        param = Float64[ 0.2, 0.8 ]

        simplex_box = Transform.SimplexBox(lb, this_scale, length(param))
        simplex_and_unsimplex(param, simplex_box)

        # Test that the edges work.
        simplex_and_unsimplex(Float64[ lb, 1 - lb ], simplex_box)
        simplex_and_unsimplex(Float64[ 1 - lb, lb ], simplex_box)

        # Test the scaling
        unscaled_simplex_box = Transform.SimplexBox(lb, 1.0, length(param))
        @test_approx_eq(
            Transform.unsimplexify_parameter(param, simplex_box),
            this_scale * Transform.unsimplexify_parameter(param, unscaled_simplex_box))

        # Test the bound checking
        @test_throws(AssertionError,
            Transform.unsimplexify_parameter([ lb - 1e-6, 1 - lb + 1e-6 ], simplex_box))
        @test_throws(AssertionError,
            Transform.unsimplexify_parameter([ 0.3, 0.8 ], simplex_box))
        @test_throws(AssertionError,
            Transform.unsimplexify_parameter([ 0.2, 0.3, 0.5 ], simplex_box))
        @test_throws(AssertionError, Transform.simplexify_parameter([ 1., 2. ], simplex_box))
    end
end

function test_transform_box_functions()
    function box_and_unbox{NumType <: Number}(param::NumType, param_box::ParamBox)
        param_free = Transform.unbox_parameter(param, param_box)
        new_param = Transform.box_parameter(param_free, param_box)
        @test_approx_eq param new_param
    end

    for this_scale = [ 1.0, 2.0 ], lb = [-10.0, 0.1], ub = [0.5, Inf]
        #println(this_scale, " ", lb, " ", ub)
        param = 0.2
        param_box = Transform.ParamBox(lb, ub, this_scale)
        box_and_unbox(param, param_box)

        # Test that the edges work.
        box_and_unbox(lb, param_box)
        ub < Inf && box_and_unbox(ub, param_box)

        # Test the scaling
        unscaled_param_box = Transform.ParamBox(lb, ub, 1.0)
        @test_approx_eq(
            Transform.unbox_parameter(param, param_box),
            this_scale * Transform.unbox_parameter(param, unscaled_param_box))

        # Test the bound checking
        @test_throws AssertionError Transform.unbox_parameter(lb - 1.0, param_box)
        ub < Inf &&
            @test_throws AssertionError Transform.unbox_parameter(ub + 1.0, param_box)
    end
end

function test_basic_transforms()
    @test_approx_eq 0.99 Transform.logit(Transform.inv_logit(0.99))
    @test_approx_eq -6.0 Transform.inv_logit(Transform.logit(-6.0))

    @test_approx_eq(
        [ 0.99, 0.001 ], Transform.logit.(Transform.inv_logit.([ 0.99, 0.001 ])))
    @test_approx_eq(
        [ -6.0, 0.5 ], Transform.inv_logit.(Transform.logit.([ -6.0, 0.5 ])))

    z = Float64[ 2.0, 4.0 ]
    z ./= sum(z)
    x = Transform.unconstrain_simplex(z)

    @test_approx_eq Transform.constrain_to_simplex(x) z

    @test Transform.inv_logit(1.0) == Inf
    @test Transform.inv_logit(0.0) == -Inf

    @test Transform.logit(Inf) == 1.0
    @test Transform.logit(-Inf) == 0.0

    @test_approx_eq Transform.constrain_to_simplex([-Inf]) [0.0, 1.0]
    @test_approx_eq Transform.unconstrain_simplex([0.0, 1.0]) [-Inf]

    @test_approx_eq Transform.constrain_to_simplex([Inf]) [1.0, 0.0]
    @test_approx_eq Transform.unconstrain_simplex([1.0, 0.0]) [Inf]

    @test_approx_eq Transform.constrain_to_simplex([Inf, 5]) [1.0, 0.0, 0.0]
end

function test_enforce_bounds()
  blob, ea_original, three_bodies = gen_three_body_dataset();
  transform = get_mp_transform(ea_original.vp, ea_original.active_sources);

  ea = deepcopy(ea_original);
  ea.vp[1][ids.a[1, 1]] = transform.bounds[1][:a][1].lb - 0.00001
  ea.vp[2][ids.r1[1]] = transform.bounds[2][:r1][1].lb - 1.0
  ea.vp[2][ids.r1[2]] = transform.bounds[2][:r1][1].ub + 1.0
  ea.vp[3][ids.k[1, 1]] = transform.bounds[3][:k][1].lb - 0.00001

  @test_throws AssertionError Transform.from_vp(transform, ea.vp)
  Transform.enforce_bounds!(ea.vp, ea.active_sources, transform)

  # Check that it now works and all values are finite.
  x_trans = Transform.from_vp(transform, ea.vp)
  for s = 1:ea.S
      @test !any(Bool[ isinf(x) for x in x_trans[1] ])
  end

  # Test a corner case of the simplex bounds to make sure that normalizing
  # doesn't violate the minimization constraints
  ea = deepcopy(ea_original);
  constraint = transform.bounds[1][:a][1]
  ea.vp[1][ids.a[1, 1]] = transform.bounds[1][:a][1].lb - 0.00001
  ea.vp[1][ids.a[2, 1]] = 100
  @test_throws AssertionError Transform.from_vp(transform, ea.vp)
  Transform.enforce_bounds!(ea.vp, ea.active_sources, transform)
  # Check that it runs without an error now
  Transform.from_vp(transform, ea.vp)

  # Test with only one active source.
  ea = deepcopy(ea_original);
  sa = 2
  ea.active_sources = [sa]
  transform = get_mp_transform(ea.vp, ea.active_sources);

  ea.vp[sa][ids.a[1, 1]] = transform.bounds[1][:a][1].lb - 0.00001
  ea.vp[sa][ids.r1[1]] = transform.bounds[1][:r1][1].lb - 1.0
  ea.vp[sa][ids.r1[2]] = transform.bounds[1][:r1][1].ub + 1.0
  ea.vp[sa][ids.k[1, 1]] = transform.bounds[1][:k][1].lb - 0.00001

  @test_throws AssertionError Transform.from_vp(transform, ea.vp)
  Transform.enforce_bounds!(ea.vp, ea.active_sources, transform)

  # Check that it now works and all values are finite.
  x_trans = Transform.from_vp(transform, ea.vp)
  for s = 1:ea.S
      @test !any(Bool[ isinf(x) for x in x_trans[1] ])
  end
end

test_transform_box_functions()
test_parameter_conversion()
test_transform_simplex_functions()
test_basic_transforms()
test_enforce_bounds()


text:
using ComputationalHomology
using LinearAlgebra
using SparseArrays
using Distances
using Test
using Random

include("cells.jl")
include("chains.jl")
include("simplicial.jl")
include("bitmap.jl")
include("cw.jl")
include("boundary.jl")
include("filtration.jl")
include("constructions.jl")
include("homology.jl")
include("persistence.jl")
include("landscape.jl")
include("pimage.jl")
include("distances.jl")
include("tropic.jl")


text:
# test the SliceSampling module

include("SliceSampling.jl")
using Distributions
using Gadfly
using Random
#### I : test the Normal Target Distribution #####
function target(x)
    pdf(Normal(10, 2), x)  # density of given x value, destracted from N(10, 2²)
end

points = SliceSampling.sliceSample(10000, target, [0, 20], 0.1)
fromTarget = rand(Normal(10, 2), 10000)
# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.histogram, Theme(default_color=colorant"green"))
plot(p1, p2, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target N(μ = 10, σ = 2)"],["orange","green"]))

# Normal distribution is fine

#### II: test the Beta Target Distribution #####
function target(x)
    pdf(Beta(2,5), x)  # density of given x value, destracted from N(10, 2²)
end

points = SliceSampling.sliceSample(10000, target, [0, 1], 0.001)
fromTarget = rand(Beta(2,5), 10000)
# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.histogram, Theme(default_color=colorant"green"))
plot(p1, p2, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target N(μ = 10, σ = 2)"],["orange","green"]))

#plot(x = fromTarget, Geom.histogram)

# Beta is fine

#### III: test the Cauchy Target Distribution #####
#### TODO: half Cauchy
function target(x)
    pdf(Cauchy(2,0.5), x)  # density of given x value, destracted from N(10, 2²)
end

fromTarget = rand(Cauchy(2,0.5), 10000)
fromTarget = fromTarget[-5 .< fromTarget .< 10]
points = SliceSampling.sliceSample(length(fromTarget), target, [-5, 10], 0.001)

# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.density, Theme(default_color=colorant"green"))
plot(p1, p2, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target N(μ = 10, σ = 2)"],["orange","green"]))

#plot(x = fromTarget, Geom.histogram)

# Cauchy not OK? TODO
# Cauchy not OK? NOT that OK:
    # the slice sampling hist is covering the Target hist

#### III.2: test the Half Cauchy Target Distribution #####
# as Half Cauchy is not available raight now
# HalfCauchy of σ² = Inverse Beta of λ > 0
include("InvBeta.jl")
rng = Random.seed!(2019)
function target(x)
    InvBeta.pdf(InvBeta.InverseBeta(.5,.5), x)  # density of given x value, destracted from N(10, 2²)
end

fromTarget = InvBeta.rand(rng, InvBeta.InverseBeta(.5, .5), 10000)
fromTarget = fromTarget[0 .< fromTarget .< 50]
points = SliceSampling.sliceSample(length(fromTarget), target, [0.000001, 50], 0.1)

# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.histogram, Theme(default_color=colorant"green"))
plot(p2, p1, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target N(μ = 10, σ = 2)"],["orange","green"]))

#plot(x = fromTarget, Geom.histogram)
# HalfCauchy/ InverseBeta is fine!

#### IV: test the Gamma Target Distribution #####
function target(x)
    pdf(Gamma(3, 2), x)  # density of given x value, destracted from N(10, 2²)
end

fromTarget = rand(Cauchy(3,2), 10000)
fromTarget = fromTarget[0 .< fromTarget .<16]
points = SliceSampling.sliceSample(length(fromTarget), target, [0, 16], 0.1)

# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.histogram, Theme(default_color=colorant"green"))
plot(p1, p2, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target"],["orange","green"]))

#plot(x = fromTarget, Geom.histogram)
# Gamma not OK: TODO
    # The sampled data has heavier tails

#### V: test the InverseGamma Target Distribution #####
function target(x)
    pdf(InverseGamma(3,1), x)  # density of given x value, destracted from N(10, 2²)
end

fromTarget = rand(InverseGamma(3,1), 10000)
fromTarget = fromTarget[0.0001 .< fromTarget .< 5]
points = SliceSampling.sliceSample(length(fromTarget), target, [0.0001, 5], 0.01)

# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.histogram, Theme(default_color=colorant"green"))
plot(p1, p2, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target N(μ = 10, σ = 2)"],["orange","green"]))

#plot(x = fromTarget, Geom.histogram)
# InverseGamma is OK !

#### V: test the Exponential Target Distribution #####
function target(x)
    pdf(Exponential(.5), x)  # density of given x value, destracted from N(10, 2²)
end

fromTarget = rand(Exponential(.5), 10000)
fromTarget = fromTarget[0 .< fromTarget .< 5]
points = SliceSampling.sliceSample(length(fromTarget), target, [0, 5], 0.01)

# to plot the results
p1 = layer(x = points, Geom.histogram, Theme(default_color=colorant"orange"))
p2 = layer(x = fromTarget, Geom.histogram, Theme(default_color=colorant"green"))
plot(p2, p1, Theme(grid_color = colorant"white", grid_color_focused = colorant"white"), Guide.manual_color_key("Sampled vs Target",["SliceSampling","Target N(μ = 10, σ = 2)"],["orange","green"]))

#plot(x = fromTarget, Geom.histogram)
# Exponential is OK !

#plot(x = fromTarget, Geom.histogram)

max_stars_repo_name: YinpuLi/SparsityBoosting

max_issues_repo_name: YinpuLi/SparsityBoosting

max_forks_repo_name: YinpuLi/SparsityBoosting


text:
abstract type AbstractQueue{T} end

### Interface

# Insert an item at the rear of the queue.
enqueue!(q::AbstractQueue, item) = nothing

# Remove and return element at the front of the queue.
dequeue!(q::AbstractQueue) = nothing

# Return the front element.
frontvalue(q::AbstractQueue) = nothing

# Return the number of elements in the queue.
length(q::AbstractQueue) = nothing

# Check if queue is empty.
isempty(q::AbstractQueue) = nothing

### Generic
eltype(q::AbstractQueue{T}) where T = T


text:
if Base.HOME_PROJECT[] !== nothing
    Base.HOME_PROJECT[] = abspath(Base.HOME_PROJECT[])
end

push!(LOAD_PATH, joinpath("..", "src"))

using Documenter, Stripeline

makedocs(
    modules = [Stripeline],
    format = Documenter.HTML(
        prettyurls = get(ENV, "CI", nothing) == "true"
    ),
    sitename = "Stripeline.jl",
    pages = [
        "Introduction" => "index.md",
        "Basic functions" => "basic.md",
        "Instrument database" => "instrumentdb.md",
        "Scanning strategy" => "scanning.md",
        "Data acquisition" => "acquisition.md",
        "Map making" => "mapmaking.md",
	    "Simulation tutorial" => "simulation_tutorial.md",
    ])

deploydocs(
    repo = "github.com/lspestrip/Stripeline.jl.git",
)


text:
# Write the sys source cache in format readable by Base._read_dependency_src
cachefile = ARGS[1]
open(cachefile, "w") do io
    for (_, filename) in Base._included_files
        src = read(filename, String)
        write(io, hton(Int32(sizeof(filename))))
        write(io, filename)
        write(io, hton(UInt64(sizeof(src))))
        write(io, src)
    end
    write(io, Int32(0))
end

max_stars_repo_name: blegat/julia

max_issues_repo_name: JuliaLangEs/julieta

max_forks_repo_name: JuliaLangEs/julieta


text:
using OVERTVerify
using LazySets
using Dates

controller_name = ARGS[1] # e.g. "big"
controller = "nnet_files/jmlr/tora_$(controller_name)_controller.nnet"
println("controller is: $controller")

query = OvertQuery(
	Tora,      # problem
	controller, # network file
	Id(),    # last layer activation layer Id()=linear, or ReLU()=relu
	"MIP",     # query solver, "MIP" or "ReluPlex"
	15,        # ntime
	0.1,       # dt
	-1,        # N_overt
	)

input_set = Hyperrectangle(low=[0.6, -0.7, -0.4, 0.5], high=[0.7, -0.6, -0.3, 0.6])
concretization_intervals = [5, 5, 5]
t1 = Dates.time()
concrete_state_sets, symbolic_state_sets, concrete_meas_sets, symbolic_meas_sets = symbolic_reachability_with_concretization(query, input_set, concretization_intervals)
t2 = Dates.time()
dt = (t2-t1)
print("elapsed time= $(dt) seconds")

# clean up sets
init_set0, reachable_state_sets = clean_up_sets(concrete_state_sets, symbolic_state_sets, concretization_intervals)

# we want to check inclusion in the safe set:
constraint1 = HalfSpace([1., 0., 0., 0.], 2.) # x1 <= 2
constraint2 = HalfSpace([-1., 0., 0., 0.], 2.) # -x1 <= 2 aka x1 >= -2
safe_set = HPolyhedron([constraint1, constraint2])
t1 = time()
safe_steps = reachable_state_sets .⊆ Ref(safe_set)
violations = .!safe_steps
safe = all(safe_steps)
dt_check = time() - t1

using JLD2
JLD2.@save "src/examples/jmlr/data/tora_reachability_$(controller_name)_controller_data.jld2" query input_set safe_set concrete_state_sets symbolic_state_sets concrete_meas_sets symbolic_meas_sets reachable_state_sets dt safe violations dt_check concretization_intervals


text:
#-----------------------------------------------------------
# Woodbury inequality for efficient inverse calculation
#-----------------------------------------------------------

"""
    woodbury_A_plus_UCV(; A⁻¹=A⁻¹, U=U, C⁻¹=C⁻¹, V=V)

Efficient computation of `(A + UCV)⁻¹`` when `A⁻¹`` is available (and `C` can be easily inverted).
It does so by applying the Woodbury identity is used:

# Examples
```julia-repl
julia> using LinearAlgebra
julia> A = randn(5,5); A = A*A'; U = randn(5,3); C = Diagonal(rand(3)); V = randn(3,5);
julia> A⁻¹ = A\\I # assume inverse is available
julia> B1 = inv(A + U*C*V) # regular way of computing inverse
julia> B2 = woodbury_A_plus_UCV(; A⁻¹=A⁻¹, U=U, C⁻¹=C\\I, V=V)
julia> maximum(abs.(B1.-B2))<1e-10 # verify that they are numerically equal, should return true
```
"""
woodbury_A_plus_UCV(; A⁻¹=A⁻¹, U=U, C⁻¹=C⁻¹, V=V) =  A⁻¹ - A⁻¹*U*((C⁻¹ + V*A⁻¹*U)\V)*A⁻¹

#-----------------------------------------------------------
# Make matrix symmetric by adding its transpose
#-----------------------------------------------------------

function makematrixsymmetric!(A)
    I,J = size(A)
    for i in 1:I
        for j in i+1:J
            @inbounds tmp = (A[i,j]+A[j,i])/2
            @inbounds A[i,j] = A[j,i] = tmp
        end
    end
    nothing
end

function makematrixsymmetric(A)
    (A+A')/2
end

#-----------------------------------------------------------
# Add a small constant, i.e. jitter to the diagonal
#-----------------------------------------------------------

function addjitter!(A, JITTER = 1e-8)

    I, J = size(A)

    if I !== J
        error("Matrix passed to  addjitter! must be square")
    end

    for i in 1:I
        @inbounds A[i,i] += JITTER
    end

    nothing
end

addjitter(A, JITTER = 1e-8) = A + JITTER*I

#-----------------------------------------------------------
# Find closest positive definite matrix
#-----------------------------------------------------------

"""
    nearestposdef(A; minimumeigenvalue = 1e-6)

Nearest positive definite matrix in Frobenious norm.
Also called the projection of the matrix onto the cone of positive definite matrices.
See: https://nhigham.com/2021/01/26/what-is-the-nearest-positive-semidefinite-matrix/
```
"""
function nearestposdef(A; minimumeigenvalue = 1e-6)

    local Evalues, Evector = eigen(A)

    local newA = Evector * Diagonal(max.(Evalues, minimumeigenvalue)) * Evector'

    makematrixsymmetric!(newA)

    newA

end


text:
###############################################
#FEniCS tutorial demo program: Incompressible Navier-Stokes equations
#for channel flow (Poisseuille) on the unit square using the
#Incremental Pressure Correction Scheme (IPCS).
#  u' + u . nabla(u)) - div(σ(u, p)) = f
#                                 div(u) = 0
###############################################

module ft07

using FenicsPy

T = 10.0           # final time
num_steps = 500    # number of time steps
Δt = T / num_steps # time step size
μ = 1             # kinematic viscosity
ρ = 1            # density

# Create mesh and define function spaces
mesh = UnitSquareMesh(16, 16)
V = VectorFunctionSpace(mesh, "P", 2)
Q = FunctionSpace(mesh, "P", 1)

# Define boundaries
inflow  = "near(x[0], 0)"
outflow = "near(x[0], 1)"
walls   = "near(x[1], 0) || near(x[1], 1)"

# Define boundary conditions
bcu_noslip  = DirichletBC(V, Constant((0, 0)), walls)
bcp_inflow  = DirichletBC(Q, Constant(8), inflow)
bcp_outflow = DirichletBC(Q, Constant(0), outflow)
bcu = [bcu_noslip]
bcp = [bcp_inflow, bcp_outflow]

# Define trial and test functions
u = TrialFunction(V)
v = TestFunction(V)
p = TrialFunction(Q)
q = TestFunction(Q)

# Define functions for solutions at previous and current time steps
u_n = FeFunction(V)
u_  = FeFunction(V)
p_n = FeFunction(Q)
p_  = FeFunction(Q)

# Define expressions used in variational forms
U   = 0.5*(u_n + u)
n   = FacetNormal(mesh)
f   = Constant((0, 0))

# Define strain-rate tensor
ϵ(u) = sym(nabla_grad(u))

# Define stress tensor
σ(u, p) = 2*μ*ϵ(u) - p*Identity(len(u))

# Define variational problem for step 1
F1 = ρ*dot((u - u_n) / Δt, v)*dx +
        ρ*dot(dot(u_n, nabla_grad(u_n)), v)*dx +
        inner(σ(U, p_n), ϵ(v))*dx +
        dot(p_n*n, v)*ds - dot(μ*nabla_grad(U)*n, v)*ds -
        dot(f, v)*dx
a1 = lhs(F1)
L1 = rhs(F1)

# Define variational problem for step 2
a2 = dot(nabla_grad(p), nabla_grad(q))*dx
L2 = dot(nabla_grad(p_n), nabla_grad(q))*dx - (1/Δt)*div(u_)*q*dx

# Define variational problem for step 3
a3 = dot(u, v)*dx
L3 = dot(u_, v)*dx - Δt*dot(nabla_grad(p_ - p_n), v)*dx

# Assemble matrices
A1 = assemble(a1)
A2 = assemble(a2)
A3 = assemble(a3)

# Apply boundary conditions to matrices
[bc.apply(A1) for bc in bcu]
[bc.apply(A2) for bc in bcp]

# Time-stepping
t = 0
for i in 1:num_steps

    # Update current time
    global t += Δt

    # Step 1: Tentative velocity step
    b1 = assemble(L1)
    [bc.apply(b1) for bc in bcu]
    solve(A1, u_.vector(), b1)

    # Step 2: Pressure correction step
    b2 = assemble(L2)
    [bc.apply(b2) for bc in bcp]
    solve(A2, p_.vector(), b2)

    # Step 3: Velocity correction step
    b3 = assemble(L3)
    solve(A3, u_.vector(), b3)

    # Plot solution
    # plot(u_)

    # Compute error
    u_e = Expression(("4*x[1]*(1.0 - x[1])", "0"), degree=2)
    u_e = interpolate(u_e, V)
    error = max(abs.(array(u_e) - array(u_))...)
    println("t = ", t, ", \t error = ", error)
    println("max u:", max(array(u_)...))

    # Update previous solution
    u_n.assign(u_)
    p_n.assign(p_)

end

end # module ft07


text:
export gnuplot
export errors

function gnuplot( filename :: String, fields :: MeshFields )

    open(filename, "w") do f

        nx = fields.mesh.nx
        ny = fields.mesh.ny
        dx = fields.mesh.dx
        dy = fields.mesh.dy

        for i in 1:nx+1
            for j in 1:ny+1
                write(f,
	    	        string((i-1)*dx),        "  ",
	    	        string((j-1)*dy),        "  ",
	    	        string(fields.e[1,i,j]), "  ",
	    	        string(fields.e[2,i,j]), "  ",
	    	        string(fields.ρ[i,j]),  "\n")
            end
	        write(f,"\n")
        end

    end

end

function errors( computed :: MeshFields, reference :: MeshFields )

    err = maximum(abs.(computed.e .- reference.e))

    err

end


text:
export xtitle, ytitle, x2title, y2title,
       xticks, yticks, x2ticks, y2ticks

####
#### [x|y|x2|y2]title[!]
####

"""
    _title(...)

Internal function to set the title of axes (`el==:title`) or axis objects
(`el==:xtitle`...).
"""
function _title(
            el::Symbol,
            text::String = "";
            #
            axes::Option{Axes2D} = nothing,
            opts...
        )::PreviewFigure

    axes = check_axes(axes)
    # retrieve on what object to apply the title either the current axes or a sub-axis
    obj = (el == :axis) ? axes : getfield(axes, el)
    # create a new title object with the text and apply properties
    obj.title = Title(; text)
    set_properties!(obj.title;  opts...)
    return PreviewFigure(gcf())
end

# Generate xtitle, and associated for each axis
for axs ∈ ("", "x", "y", "x2", "y2")
    f   = Symbol(axs * "title")
    f2  = Symbol(axs * "label")
    ex = quote
        # mutate
        function $f(t::String = ""; o...)
            _title(
                Symbol($axs * "axis"), t;
                axes = gca(), o...
            )
        end
        # synonyms xlabel, ylabel, etc
        !isempty($axs) && ($f2 = $f)
    end
    eval(ex)
end

####
#### [x|y|x2|y2]ticks
####

function _ticks(
            axis_sym::Symbol,
            loc::Vector{Float64} = Float64[],
            lab::Vector{String}  = String[];
            #
            axes::Option{Axes2D} = nothing,
            opts...
        )::PreviewFigure

    axes = check_axes(axes)
    # retrieve the appropriate axis
    axis = getfield(axes, axis_sym)
    # create a new ticks object
    axis.ticks = Ticks()
    # if locs are empty, just pass options and return
    if isempty(loc)
        isempty(lab) || throw(
            ArgumentError("""
                Cannot pass ticks labels without specifying ticks locations.
                """
            )
        )
        set_properties!(axis.ticks;  opts...)
        return PreviewFigure(gcf())
    end
    axis.ticks.places = loc
    # process labels if any are passed
    if !isempty(lab)
        length(lab) == length(loc) || throw(
            ArgumentError("""
                Ticks locations and labels must have the same length.
                """
            )
        )
        axis.ticks.labels = TicksLabels(names=lab)
    end
    # set remaining properties and return
    set_properties!(axis.ticks;  opts...)
    return PreviewFigure(gcf())
end

# Generate xticks, and associated for each axis
for axs ∈ ("x", "y", "x2", "y2")
    f  = Symbol(axs * "ticks")
    a  = Symbol(axs * "axis")
    ex = quote
        # xticks([1, 2], ["a", "b"]; o...)
        function $f(loc::AVR, lab=String[]; o...)
            _ticks(
                Symbol($axs * "axis"), Float64.(loc), lab;
                o...
            )
        end
        # xticks("off"; o...)
        function $f(s::String=""; o...)
            isempty(s) && return _ticks(Symbol($axs * "axis"); o...)
            ax = getfield(gca(), Symbol($axs * "axis"))
            s_lc = lowercase(s)
            if s_lc == "off"
                reset!(ax.ticks)
                ax.ticks.off = true
                ax.ticks.labels.off = true
            elseif s_lc == "on"
                ax.ticks.off = false
                ax.ticks.labels.off = false
                set_properties!(ax.ticks;  o...)
            else
                throw(ArgumentError("Unrecognised shorthand: $s"))
            end
            return PreviewFigure(gcf())
        end
    end
    eval(ex)
end


text:
function pull_request_build(::NewVersion,
                            pr::GitHub.PullRequest,
                            current_pr_head_commit_sha::String,
                            registry::GitHub.Repo;
                            auth::GitHub.Authorization,
                            authorized_authors::Vector{String},
                            registry_head::String,
                            registry_master::String,
                            suggest_onepointzero::Bool,
                            whoami::String,
                            registry_deps::Vector{<:AbstractString} = String[])::Nothing
    # first check if the PR is open, and the author is authorized - if not, then quit
    # then, delete ALL reviews by me
    # then check rules 1-6. if fail, post comment.
    # if everything passed, add an approving review by me
    #
    # Rules:
    # 1. Only changes a subset of the following files:
    #     - `E/Example/Compat.toml`
    #     - `E/Example/Deps.toml`
    #     - `E/Example/Versions.toml`
    # 2. Sequential version number
    #     - if the last version was 1.2.3 then the next can be 1.2.4, 1.3.0 or 2.0.0
    #     - does not apply to JLL packages
    # 3. Compat for all dependencies
    #     - there should be a [compat] entry for Julia
    #     - all [deps] should also have [compat] entries
    #     - all [compat] entries should have upper bounds
    #     - dependencies that are standard libraries do not need [compat] entries
    #     - dependencies that are JLL packages do not need [compat] entries
    # 4. If it is a patch release, then it does not narrow the Julia compat range
    # 5. (only applies to JLL packages) The only dependencies of the package are:
    #     - Pkg
    #     - Libdl
    #     - other JLL packages
    # 6. Version can be installed
    #     - given the proposed changes to the registry, can we resolve and install the new version of the package?
    #     - i.e. can we run `Pkg.add("Foo")`
    # 7. Version can be loaded
    #     - once it's been installed (and built?), can we load the code?
    #     - i.e. can we run `import Foo`
    pkg, version = parse_pull_request_title(NewVersion(), pr)
    this_is_jll_package = is_jll_name(pkg)
    @info("This is a new package pull request", pkg, version, this_is_jll_package)
    pr_author_login = author_login(pr)
    if is_open(pr)
        if pr_author_login in authorized_authors
            description = "New version. Pending."
            params = Dict("state" => "pending",
                          "context" => "automerge/decision",
                          "description" => description)
            my_retry(() -> GitHub.create_status(registry,
                                                current_pr_head_commit_sha;
                                                auth = auth,
                                                params=params))
            g1, m1 = pr_only_changes_allowed_files(NewVersion(),
                                                   registry,
                                                   pr,
                                                   pkg;
                                                   auth = auth)
            if this_is_jll_package
                g2 = true
                m2 = ""
                release_type = :jll_release
            else
                g2, m2, release_type = meets_sequential_version_number(pkg,
                                                                       version;
                                                                       registry_head = registry_head,
                                                                       registry_master = registry_master)
            end

            g3, m3 = meets_compat_for_all_deps(registry_head,
                                               pkg,
                                               version)
            g4_if_patch, m4_if_patch = meets_patch_release_does_not_narrow_julia_compat(pkg,
                                                                                version;
                                                                                registry_head = registry_head,
                                                                                registry_master = registry_master)
            if release_type == :patch
                g4 = g4_if_patch
                m4 = m4_if_patch
            else
                g4 = true
                m4 = ""
            end
            g5_if_jll, m5_if_jll = meets_allowed_jll_nonrecursive_dependencies(registry_head,
                                                                               pkg,
                                                                               version)
            if this_is_jll_package
                g5 = g5_if_jll
                m5 = m5_if_jll
            else
                g5 = true
                m5 = ""
            end
            @info("Only modifies the files that it's allowed to modify",
                  meets_this_guideline = g1,
                  message = m1)
            @info("Sequential version number",
                  meets_this_guideline = g2,
                  message = m2)
            @info("Compat (with upper bound) for all dependencies",
                  meets_this_guideline = g3,
                  message = m3)
            @info("If it is a patch release, then it does not narrow the Julia compat range",
                  meets_this_guideline = g4,
                  message = m4)
            @info("If this is a JLL package, only deps are Pkg, Libdl, and other JLL packages",
                  meets_this_guideline = g5,
                  message = m5)
            g1through5 = Bool[g1, g2, g3, g4, g5]
            if !all(g1through5)
                description = "New version. Failed."
                params = Dict("state" => "failure",
                              "context" => "automerge/decision",
                              "description" => description)
                my_retry(() -> GitHub.create_status(registry,
                                                    current_pr_head_commit_sha;
                                                    auth = auth,
                                                    params = params))
            end
            g6, m6 = meets_version_can_be_pkg_added(registry_head,
                                                    pkg,
                                                    version;
                                                    registry_deps = registry_deps)
            @info("Version can be `Pkg.add`ed",
                  meets_this_guideline = g6,
                  message = m6)
            g7, m7 = meets_version_can_be_imported(registry_head,
                                                   pkg,
                                                   version;
                                                   registry_deps = registry_deps)
            @info("Version can be `import`ed",
                  meets_this_guideline = g7,
                  message = m7)
            g1through7 = Bool[g1, g2, g3, g4, g5, g6, g7]
            allmessages1through7 = String[m1, m2, m3, m4, m5, m6, m7]
            if all(g1through7) # success
                description = "New version. Approved. name=\"$(pkg)\". sha=\"$(current_pr_head_commit_sha)\""
                params = Dict("state" => "success",
                              "context" => "automerge/decision",
                              "description" => description)
                my_retry(() -> GitHub.create_status(registry,
                                                    current_pr_head_commit_sha;
                                                    auth = auth,
                                                    params = params))
                this_pr_comment_pass = comment_text_pass(NewVersion(),
                                                         suggest_onepointzero,
                                                         version)
                my_retry(() -> update_automerge_comment!(registry,
                                                         pr;
                                                         auth = auth,
                                                         body = this_pr_comment_pass,
                                                         whoami = whoami))
                return nothing
            else # failure
                description = "New version. Failed."
                params = Dict("state" => "failure",
                              "context" => "automerge/decision",
                              "description" => description)
                my_retry(() -> GitHub.create_status(registry,
                                                    current_pr_head_commit_sha;
                                                    auth=auth,
                                                    params=params))
                failingmessages1through7 = allmessages1through7[.!g1through7]
                this_pr_comment_fail = comment_text_fail(NewVersion(),
                                                         failingmessages1through7,
                                                         suggest_onepointzero,
                                                         version)
                my_retry(() -> update_automerge_comment!(registry,
                                                         pr;
                                                         body = this_pr_comment_fail,
                                                         auth = auth,
                                                         whoami = whoami))
                throw(AutoMergeGuidelinesNotMet("The automerge guidelines were not met."))
            end
        else
            throw(AutoMergeAuthorNotAuthorized("Author $(pr_author_login) is not authorized to automerge. Exiting..."))
        end
    else
        throw(AutoMergePullRequestNotOpen("The pull request is not open. Exiting..."))
    end
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
# https://discourse.julialang.org/t/vector-of-matrices-vs-multidimensional-arrays/9602/5
# 5d2360465457a3f77ddc131e has TOU demand
# 59bc22705457a3372642da67 has monthly tiered demand (no TOU demand)

"""
    Base.@kwdef struct URDBrate

Contains some of the data for ElectricTariff
"""
struct URDBrate
    energy_rates::Array{Float64,2}  # time X tier
    energy_tier_limits::Array{Real,1}
    n_energy_tiers::Int

    n_monthly_demand_tiers::Int
    monthly_demand_tier_limits::Array{Real,1}
    monthly_demand_rates::Array{Float64,2}  # month X tier

    n_tou_demand_tiers::Int
    tou_demand_tier_limits::Array{Real,1}
    tou_demand_rates::Array{Float64,2}  # ratchet X tier
    tou_demand_ratchet_time_steps::Array{Array{Int64,1},1}  # length = n_tou_demand_ratchets

    demand_lookback_months::AbstractArray{Int,1}
    demand_lookback_percent::Float64
    demand_lookback_range::Int

    fixed_monthly_charge::Float64
    annual_min_charge::Float64
    min_monthly_charge::Float64

    sell_rates::Array{Float64,2}  # time X tier
end

"""
    URDBrate(urdb_label::String, year::Int=2019; time_steps_per_hour=1)

download URDB dict, parse into reopt inputs, return URDBrate struct.
    year is required to align weekday/weekend schedules.
"""
function URDBrate(urdb_label::String, year::Int=2019; time_steps_per_hour=1)
    urdb_response = download_urdb(urdb_label)
    URDBrate(urdb_response, year; time_steps_per_hour=time_steps_per_hour)
end

"""
    URDBrate(util_name::String, rate_name::String, year::Int=2019; time_steps_per_hour=1)

download URDB dict, parse into reopt inputs, return URDBrate struct.
    year is required to align weekday/weekend schedules.
"""
function URDBrate(util_name::String, rate_name::String, year::Int=2019; time_steps_per_hour=1)
    urdb_response = download_urdb(util_name, rate_name)
    URDBrate(urdb_response, year; time_steps_per_hour=time_steps_per_hour)
end

"""
    URDBrate(urdb_response::Dict, year::Int)

process URDB response dict, parse into reopt inputs, return URDBrate struct.
    year is required to align weekday/weekend schedules.
"""
function URDBrate(urdb_response::Dict, year::Int=2019; time_steps_per_hour=1)

    demand_min = get(urdb_response, "peakkwcapacitymin", 0.0)  # TODO add check for site min demand against tariff?

    n_monthly_demand_tiers, monthly_demand_tier_limits, monthly_demand_rates,
      n_tou_demand_tiers, tou_demand_tier_limits, tou_demand_rates, tou_demand_ratchet_time_steps =
      parse_demand_rates(urdb_response, year)

    energy_rates, energy_tier_limits, n_energy_tiers, sell_rates =
        parse_urdb_energy_costs(urdb_response, year; time_steps_per_hour=time_steps_per_hour)

    fixed_monthly_charge, annual_min_charge, min_monthly_charge = parse_urdb_fixed_charges(urdb_response)

    demand_lookback_months, demand_lookback_percent, demand_lookback_range = parse_urdb_lookback_charges(urdb_response)

    URDBrate(
        energy_rates,
        energy_tier_limits,
        n_energy_tiers,

        n_monthly_demand_tiers,
        monthly_demand_tier_limits,
        monthly_demand_rates,

        n_tou_demand_tiers,
        tou_demand_tier_limits,
        tou_demand_rates,
        tou_demand_ratchet_time_steps,

        demand_lookback_months,
        demand_lookback_percent,
        demand_lookback_range,

        fixed_monthly_charge,
        annual_min_charge,
        min_monthly_charge,

        sell_rates
    )
end

function download_urdb(urdb_label::String; version::Int=8)
    url = string("https://api.openei.org/utility_rates", "?api_key=", urdb_key,
                "&version=", version , "&format=json", "&detail=full",
                "&getpage=", urdb_label
    )
    response = nothing
    try
        @info "Checking URDB for " urdb_label
        r = HTTP.get(url, require_ssl_verification=false)  # cannot verify on NREL VPN
        response = JSON.parse(String(r.body))
        if r.status != 200
            error("Bad response from URDB: $(response["errors"])")  # TODO URDB has "errors"?
        end
    catch e
        error("Error occurred :$(e)")
    end

    rates = response["items"]  # response['items'] contains a vector of dicts
    if length(rates) == 0
        error("Could not find $(urdb_label) in URDB.")
    end
    if rates[1]["label"] == urdb_label
        return rates[1]
    else
        error("Could not find $(urdb_label) in URDB.")
    end
end

function download_urdb(util_name::String, rate_name::String; version::Int=8)
    url = string("https://api.openei.org/utility_rates", "?api_key=", urdb_key,
                "&version=", version , "&format=json", "&detail=full",
                "&ratesforutility=", replace(util_name, "&" => "%26")
    )
    response = nothing
    try
        @info "Checking URDB for " rate_name
        r = HTTP.get(url, require_ssl_verification=false)  # cannot verify on NREL VPN
        response = JSON.parse(String(r.body))
        if r.status != 200
            error("Bad response from URDB: $(response["errors"])")  # TODO URDB has "errors"?
        end
    catch e
        error("Error occurred :$(e)")
    end

    rates = response["items"]  # response['items'] contains a vector of dicts
    if length(rates) == 0
        error("Could not find $(rate_name) in URDB.")
    end

    matched_rates = []
    start_dates = []

    for rate in rates
        if contains(rate_name, rate["name"])
            push!(matched_rates, rate)  # urdb can contain multiple rates of same name
            if "startdate" in keys(rate)
                push!(start_dates, rate["startdate"])
            end
        end
    end

    # find the newest rate of those that match the rate_name
    newest_index = 1  # covers the case where one rate is returned without a "startdate"
    if length(start_dates) > 1 && length(start_dates) == length(matched_rates)
        _, newest_index = findmax(start_dates)
    end

    if length(matched_rates) == 0
        error("Could not find $(rate_name) in URDB.")
    end

    return matched_rates[newest_index]
end

"""
    parse_urdb_energy_costs(d::Dict, year::Int; time_steps_per_hour=1, bigM = 1.0e8)

use URDB dict to return rates, energy_cost_vector, energy_tier_limits_kwh where:
    - rates is vector summary of rates within URDB (used for average rates when necessary)
    - energy_cost_vector is a vector of vectors with inner vectors for each energy rate tier,
        inner vectors are costs in each time step
    - energy_tier_limits_kwh is a vector of upper kWh limits for each energy tier
"""
function parse_urdb_energy_costs(d::Dict, year::Int; time_steps_per_hour=1, bigM = 1.0e8)
    if length(d["energyratestructure"]) == 0
        error("No energyratestructure in URDB response.")
    end
    # TODO check bigM (in multiple functions)
    energy_tiers = Float64[]
    for energy_rate in d["energyratestructure"]
        append!(energy_tiers, length(energy_rate))
    end
    energy_tier_set = Set(energy_tiers)
    if length(energy_tier_set) > 1
        @warn "energy periods contain different numbers of tiers, using limits of period with most tiers"
    end
    period_with_max_tiers = findall(energy_tiers .== maximum(energy_tiers))[1]
    n_energy_tiers = Int(maximum(energy_tier_set))

    rates = Float64[]
    energy_tier_limits_kwh = Float64[]
    non_kwh_units = false

    for energy_tier in d["energyratestructure"][period_with_max_tiers]
        # energy_tier is a dictionary, eg. {'max': 1000, 'rate': 0.07531, 'adj': 0.0119, 'unit': 'kWh'}
        energy_tier_max = get(energy_tier, "max", bigM)

        if "rate" in keys(energy_tier) || "adj" in keys(energy_tier)  || "sell" in keys(energy_tier)
            append!(energy_tier_limits_kwh, energy_tier_max)
        end

        if "unit" in keys(energy_tier) && string(energy_tier["unit"]) != "kWh"
            @warn "Using average rate in tier due to exotic units of " energy_tier["unit"]
            non_kwh_units = true
        end

        append!(rates, get(energy_tier, "rate", 0) + get(energy_tier, "adj", 0))
    end

    if non_kwh_units
        rate_average = sum(rates) / maximum([length(rates), 1])
        n_energy_tiers = 1
        energy_tier_limits_kwh = Float64[bigM]
    end

    energy_cost_vector = Float64[]
    sell_vector = Float64[]

    for tier in range(1, stop=n_energy_tiers)

        for month in range(1, stop=12)
            n_days = daysinmonth(Date(string(year) * "-" * string(month)))

            for day in range(1, stop=n_days)

                for hour in range(1, stop=24)

                    # NOTE: periods are zero indexed
                    if dayofweek(Date(year, month, day)) < 6  # Monday == 1
                        period = d["energyweekdayschedule"][month][hour] + 1
                    else
                        period = d["energyweekendschedule"][month][hour] + 1
                    end
                    # workaround for cases where there are different numbers of tiers in periods
                    n_tiers_in_period = length(d["energyratestructure"][period])
                    if n_tiers_in_period == 1
                        tier_use = 1
                    elseif tier > n_tiers_in_period
                        tier_use = n_tiers_in_period
                    else
                        tier_use = tier
                    end
                    if non_kwh_units
                        rate = rate_average
                    else
                        rate = get(d["energyratestructure"][period][tier_use], "rate", 0)
                    end
                    total_rate = rate + get(d["energyratestructure"][period][tier_use], "adj", 0)
                    sell = get(d["energyratestructure"][period][tier_use], "sell", 0)

                    for step in range(1, stop=time_steps_per_hour)  # repeat hourly rates intrahour
                        append!(energy_cost_vector, round(total_rate, digits=6))
                        append!(sell_vector, round(-sell, digits=6))
                    end
                end
            end
        end
    end
    energy_rates = reshape(energy_cost_vector, (:, n_energy_tiers))
    sell_rates = reshape(sell_vector, (:, n_energy_tiers))
    return energy_rates, energy_tier_limits_kwh, n_energy_tiers, sell_rates
end

"""
    parse_demand_rates(d::Dict, year::Int; bigM=1.0e8)

Parse monthly ("flat") and TOU demand rates
    can modify URDB dict when there is inconsistent numbers of tiers in rate structures
"""
function parse_demand_rates(d::Dict, year::Int; bigM=1.0e8)

    if haskey(d, "flatdemandstructure")
        scrub_urdb_demand_tiers!(d["flatdemandstructure"])
        monthly_demand_tier_limits = parse_urdb_demand_tiers(d["flatdemandstructure"])
        n_monthly_demand_tiers = length(monthly_demand_tier_limits)
        monthly_demand_rates = parse_urdb_monthly_demand(d, n_monthly_demand_tiers)
    else
        monthly_demand_tier_limits = []
        n_monthly_demand_tiers = 1
        monthly_demand_rates = Array{Float64,2}(undef, 0, 0)
    end

    if haskey(d, "demandratestructure")
        scrub_urdb_demand_tiers!(d["demandratestructure"])
        tou_demand_tier_limits = parse_urdb_demand_tiers(d["demandratestructure"])
        n_tou_demand_tiers = length(tou_demand_tier_limits)
        ratchet_time_steps, tou_demand_rates = parse_urdb_tou_demand(d, year=year, n_tiers=n_tou_demand_tiers)
    else
        tou_demand_tier_limits = []
        n_tou_demand_tiers = 0
        ratchet_time_steps = []
        tou_demand_rates = Array{Float64,2}(undef, 0, 0)
    end

    return n_monthly_demand_tiers, monthly_demand_tier_limits, monthly_demand_rates,
           n_tou_demand_tiers, tou_demand_tier_limits, tou_demand_rates, ratchet_time_steps

end

"""
    scrub_urdb_demand_tiers!(A::Array)

validate flatdemandstructure and demandratestructure have equal number of tiers across periods
"""
function scrub_urdb_demand_tiers!(A::Array)
    if length(A) == 0
        return
    end
    len_tiers = Int[length(r) for r in A]
    len_tiers_set = Set(len_tiers)
    n_tiers = maximum(len_tiers_set)

    if length(len_tiers_set) > 1
        @warn """Demand rate structure has varying number of tiers in periods.
                 Making the number of tiers the same across all periods by repeating the last tier."""
        for (i, rate) in enumerate(A)
            n_tiers_in_period = length(rate)
            if n_tiers_in_period != n_tiers
                rate_new = rate
                last_tier = rate[n_tiers_in_period]
                for j in range(1, stop=n_tiers - n_tiers_in_period)
                    append!(rate_new, last_tier)
                end
                A[i] = rate_new
            end
        end
    end
end

"""
    parse_urdb_demand_tiers(A::Array; bigM=1.0e8)

set up and validate demand tiers
    returns demand_tiers::Array{Float64, n_tiers}
"""
function parse_urdb_demand_tiers(A::Array; bigM=1.0e8)
    if length(A) == 0
        return []
    end
    len_tiers = Int[length(r) for r in A]
    n_tiers = maximum(len_tiers)
    period_with_max_tiers = findall(len_tiers .== maximum(len_tiers))[1]

    # set up tiers and validate that the highest tier has the same value across periods
    demand_tiers = Dict()
    demand_maxes = Float64[]
    for period in range(1, stop=length(A))
        demand_max = Float64[]
        for tier in A[period]
            append!(demand_max, get(tier, "max", bigM))
        end
        demand_tiers[period] = demand_max
        append!(demand_maxes, demand_max[end])  # TODO should this be maximum(demand_max)?
    end

    # test if the highest tier is the same across all periods
    if length(Set(demand_maxes)) > 1
        @warn "Highest demand tiers do not match across periods: using max tier from largest set of tiers."
    end
    return demand_tiers[period_with_max_tiers]
end

"""
    parse_urdb_monthly_demand(d::Dict)

return monthly demand rates as array{month, tier}
"""
function parse_urdb_monthly_demand(d::Dict, n_tiers)
    if !haskey(d, "flatdemandmonths")
        return []
    end
    if length(d["flatdemandmonths"]) == 0
        return []
    end

    demand_rates = zeros(12, n_tiers)  # array(month, tier)
    for month in range(1, stop=12)
        period = d["flatdemandmonths"][month] + 1  # URDB uses zero-indexing
        rates = d["flatdemandstructure"][period]  # vector of dicts

        for (t, tier) in enumerate(rates)
            demand_rates[month, t] = round(get(tier, "rate", 0.0) + get(tier, "adj", 0.0), digits=6)
        end
    end
    return demand_rates
end

"""
    parse_urdb_tou_demand(d::Dict; year::Int, n_tiers::Int)

return array of arrary for ratchet time steps, tou demand rates array{ratchet, tier}
"""
function parse_urdb_tou_demand(d::Dict; year::Int, n_tiers::Int)
    if !haskey(d, "demandratestructure")
        return [], []
    end
    n_periods = length(d["demandratestructure"])
    ratchet_time_steps = Array[]
    rates_vec = Float64[]  # array(ratchet_num, tier), reshape later
    n_ratchets = 0  # counter

    for month in range(1, stop=12)
        for period in range(0, stop=n_periods)
            time_steps = get_tou_demand_steps(d, year=year, month=month, period=period-1)
            if length(time_steps) > 0  # can be zero! not every month contains same number of periods
                n_ratchets += 1
                append!(ratchet_time_steps, [time_steps])
                for (t, tier) in enumerate(d["demandratestructure"][period])
                    append!(rates_vec, round(get(tier, "rate", 0.0) + get(tier, "adj", 0.0), digits=6))
                end
            end
        end
    end
    rates = reshape(rates_vec, (:, n_tiers))  # Array{Float64,2}
    ratchet_time_steps = convert(Array{Array{Int64,1},1}, ratchet_time_steps)
    return ratchet_time_steps, rates
end

"""
    get_tou_demand_steps(d::Dict; year::Int, month::Int, period::Int, time_steps_per_hour=1)

return Array{Int, 1} for time_steps in ratchet (aka period)
"""
function get_tou_demand_steps(d::Dict; year::Int, month::Int, period::Int, time_steps_per_hour=1)
    step_array = Int[]
    start_step = 1
    start_hour = 1

    if month > 1
        plus_days = 0
        for m in range(1, stop=month-1)
            plus_days += daysinmonth(Date(string(year) * "-" * string(m)))
            if m == 2 && isleapyear(year)
                plus_days -= 1
            end
        end
        start_hour += plus_days * 24
        start_step = start_hour * time_steps_per_hour
    end

    hour_of_year = start_hour
    step_of_year = start_step

    for day in range(1, stop=daysinmonth(Date(string(year) * "-" * string(month))))
        for hour in range(1, stop=24)
            if dayofweek(Date(year, month, day)) < 6 &&
               d["demandweekdayschedule"][month][hour] == period
                append!(step_array, step_of_year)
            elseif dayofweek(Date(year, month, day)) > 5 &&
               d["demandweekendschedule"][month][hour] == period
                append!(step_array, step_of_year)
            end
            step_of_year += 1
        end
        hour_of_year += 1
    end
    return step_array
end

"""
    parse_urdb_fixed_charges(d::Dict)

return fixed_monthly, annual_min, min_monthly :: Float64
"""
function parse_urdb_fixed_charges(d::Dict)
    fixed_monthly = 0.0
    annual_min = 0.0
    min_monthly = 0.0

    # first try $/month, then check if $/day exists, as of 1/28/2020 there were only $/day and $month entries in the URDB
    fixed_monthly = Float64(get(d, "fixedmonthlycharge", 0.0))
    if fixed_monthly == 0.0
        if get(d, "fixedchargeunits", "") == "\$/month"
            fixed_monthly = Float64(get(d, "fixedchargefirstmeter", 0.0))
        elseif get(d, "fixedchargeunits", "") == "\$/day"
            fixed_monthly = Float64(get(d, "fixedchargefirstmeter", 0.0) * 30.4375)
            # scalar intended to approximate annual charges over 12 month period, derived from 365.25/12
        elseif get(d, "fixedchargeunits", "") == "\$/year"
            fixed_monthly = Float64(get(d, "fixedchargefirstmeter", 0.0) / 12)
        elseif !isnothing(get(d, "fixedchargefirstmeter",  nothing))
            @warn "A valid value for fixedchargeunits (\$/month, \$/day, or \$/year) was not provided in urdb_response so the value provided for fixedchargefirstmeter will be ignored."
        end
    end

    if get(d, "minchargeunits", "") == "\$/month"
        min_monthly = Float64(get(d, "mincharge", 0.0))
        # first try $/month, then check if $/day or $/year exists, as of 1/28/2020 these were the only unit types in the urdb
    elseif get(d, "minchargeunits", "") == "\$/day"
        min_monthly = Float64(get(d, "mincharge", 0.0) * 30.4375 )
        # scalar intended to approximate annual charges over 12 month period, derived from 365.25/12
    elseif get(d, "minchargeunits", "") == "\$/year"
        annual_min = Float64(get(d, "mincharge", 0.0))
    elseif !isnothing(get(d, "minchargeunits",  nothing))
        @warn "A valid value for minchargeunits (\$/month, \$/day, or \$/year) was not provided in urdb_response so the value provided for mincharge will be ignored."
    end

    return fixed_monthly, annual_min, min_monthly
end

"""
    parse_urdb_lookback_charges(d::Dict)

URDB lookback fields:
- lookbackMonths
    - Type: array
    - Array of 12 booleans, true or false, indicating months in which lookbackPercent applies.
        If any of these is true, lookbackRange should be zero.
- lookbackPercent
    - Type: decimal
    - Lookback percentage. Applies to either lookbackMonths with value=1, or a lookbackRange.
- lookbackRange
    - Type: integer
    - Number of months for which lookbackPercent applies. If not 0, lookbackMonths values should all be 0.
"""
function parse_urdb_lookback_charges(d::Dict)
    lookback_months = get(d, "lookbackMonths", Int[])
    lookback_percent = Float64(get(d, "lookbackPercent", 0.0))
    lookback_range = Int64(get(d, "lookbackRange", 0.0))

    reopt_lookback_months = Int[]
    if lookback_range != 0 && length(lookback_months) == 12
        for mth in range(1, stop=12)
            if lookback_months[mth] == 1
                push!(reopt_lookback_months, mth)
            end
        end
    end
    return reopt_lookback_months, lookback_percent, lookback_range
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
"""
Authentication mechanisms
Only SASL-Plain supported for now
"""
mutable struct HiveAuth
    mechanism::AbstractString
    callback::Function

    HiveAuth(mechanism::AbstractString=SASL_MECH_PLAIN, callback::Function=Thrift.sasl_callback_default) = new(mechanism, callback)
end

function HiveAuthSASLPlain(uid::AbstractString, passwd::AbstractString; zid::AbstractString="")
    function callback(part::Symbol)
        (part == :authcid) && (return uid)
        (part == :passwd) && (return passwd)
        (part == :show) && (return uid)
        (part == :mechanism) && (return "SASL-Plain")
        return zid
    end
    HiveAuth(SASL_MECH_PLAIN, callback)
end

function show(io::IO, auth::HiveAuth)
    uid = auth.callback(:show)
    mech = auth.callback(:mechanism)
    println(io, "HiveAuth ($mech): $uid")
    nothing
end

"""
HiveConn holds the thrift connection and protocol objects.
It also holds the hive session handle for this connection.
"""
mutable struct HiveConn
    transport::TTransport
    protocol::TProtocol
    client::TCLIServiceClient
    handle::TOpenSessionResp
    connstr::AbstractString

    function HiveConn(host::AbstractString, port::Integer, auth::HiveAuth; tprotocol::Symbol=:binary)
        transport = TSASLClientTransport(TSocket(host, port), auth.mechanism, auth.callback)
        if tprotocol === :binary
            protocol = TBinaryProtocol(transport, true)
        elseif tprotocol === :compact
            protocol = TCompactProtocol(transport)
        else
            error("Unsupported protocol: $tprotocol")
        end
        client = TCLIServiceClient(protocol)
        uid = auth.callback(:show)
        connstr = "hive2://$(uid)@$(host):$port"
        new(transport, protocol, client, connect(transport, client), connstr)
    end

    function connect(transport::TTransport, client::TCLIServiceClient)
        open(transport)
        request = thriftbuild(TOpenSessionReq, Dict(:client_protocol => TProtocolVersion.HIVE_CLI_SERVICE_PROTOCOL_V8))
        response = OpenSession(client, request)
        check_status(response.status)
        response
    end
end

function show(io::IO, conn::HiveConn)
    println(io, conn.connstr)
    nothing
end

function close(conn::HiveConn)
    request = thriftbuild(TCloseSessionReq, Dict(:sessionHandle => conn.handle.sessionHandle))
    response = CloseSession(conn.client, request)
    check_status(response.status)
end

"""
HiveSession holds a connection and session status
"""
mutable struct HiveSession
    conn::HiveConn

    function HiveSession(host::AbstractString="localhost", port::Integer=10000, auth::HiveAuth=HiveAuth(); tprotocol::Symbol=:binary)
        new(HiveConn(host, port, auth; tprotocol=tprotocol))
    end
end

function show(io::IO, sess::HiveSession)
    print(io, "HiveSession: ")
    show(io, sess.conn)
end
close(session::HiveSession) = close(session.conn)

"""
Info Type
"""
const InfoType = TGetInfoType
function get_info(session::HiveSession, info_type::Int32)
    conn = session.conn
    request = thriftbuild(TGetInfoReq, Dict(:sessionHandle => conn.handle.sessionHandle, :infoType => info_type))
    response = GetInfo(conn.client, request)
    check_status(response.status)

    val = response.infoValue
    for fldname in fieldnames(TGetInfoValue)
        isfilled(val, fldname) && (return getfield(val, fldname))
    end
    nothing
end


text:
using Documenter, GeoDatasets

makedocs(modules = [GeoDatasets], sitename = "GeoDatasets.jl")

deploydocs(
    repo = "github.com/JuliaGeo/GeoDatasets.jl.git",
)

max_stars_repo_name: JuliaGeo/GeoDatasets

max_issues_repo_name: JuliaGeo/GeoDatasets

max_forks_repo_name: JuliaGeo/GeoDatasets


text:
const ALPHABET = [Char(x) for x in vcat(collect.([42:43, 48:57, 65:90, 97:122])...)]
const INV_ALPHABET = Dict(c => i for (i,c) in enumerate(ALPHABET))

_segment_width(l::Integer) = ceil(Int, log2(l+1))
encode(i::Integer, l::Integer) = string(i, base=2, pad=_segment_width(l))
function decode(c::AbstractString, l::Integer)
    k = min(_segment_width(l), length(c))
    parse(Int, c[1:k], base=2), c[k+1:end]
end

function stringify(c::AbstractString)
    if length(c) % 6 != 0
        c = c * '0' ^ mod(6-length(c), 6)
    end
    join(ALPHABET[parse(Int, x, base=2) + 1] for x in [c[i:i+5] for i in 1:6:(length(c)-1)])
end

function destringify(c::AbstractString)
    join(string(INV_ALPHABET[x] - 1, base=2, pad=6) for x in c)
end

descendants(m::BagModel) = [m.im]
descendants(n::AbstractBagNode) = [n.data]
descendants(m::ProductModel) = m.ms
descendants(n::AbstractTreeNode) = n.data
descendants_n(m) = length(descendants(m))

function ith_child(m::T, i::Integer) where T
    try
        return descendants(m)[i]
    catch
        @error "Invalid index!"
    end
end

function _walk(m::Union{ArrayModel, ArrayNode}, c::AbstractString)
    !isempty(c) || return m
    i, nc = decode(c, length(c))
    if i == 0 && Set(nc) ⊆ ['0']
        return m
    else
        @error "Invalid index!"
    end
end

function _walk(m::Union{BagModel, AbstractBagNode}, c::AbstractString)
    !isempty(c) || return m
    i, nc = decode(c, 1)
    if i == 0
        if Set(nc) ⊆ ['0']
            return m
        else
            @error "Invalid index!"
        end
    end
    _walk(descendants(m)[1], nc)
end

function _walk(m::ProductModel, c::AbstractString)
    !isempty(c) || return m
    i, nc = decode(c, length(m.ms))
    if i == 0
        if Set(nc) ⊆ ['0']
            return m
        else
            @error "Invalid index!"
        end
    end
    _walk(m.ms[i], nc)
end

function _walk(n::AbstractTreeNode, c::AbstractString)
    !isempty(c) || return n
    i, nc = decode(c, length(n.data))
    0 <= i <= length(n.data) || @error "Invalid index!"
    if i == 0
        if Set(nc) ⊆ ['0']
            return n
        else
            @error "Invalid index!"
        end
    end
    _walk(n.data[i], nc)
end

function list_traversal(m::Union{ArrayNode, ArrayModel}, s::String="")
    [stringify(s)]
end

function list_traversal(m::Union{AbstractNode, MillModel}, s::String="")
    d = descendants(m)
    n = length(d)
    vcat(stringify(s), [list_traversal(d[i], s * encode(i, n)) for i in 1:n]...)
end

show_traversal(n::AbstractNode) = dsprint(Base.stdout, n, tr=true)
show_traversal(m::MillModel) = modelprint(Base.stdout, m, tr=true)
Base.getindex(n::AbstractNode, i::AbstractString) = i == "" ? n : _walk(n, destringify(i))
Base.getindex(m::MillModel, i::AbstractString) = i == "" ? m : _walk(m, destringify(i))

tr_repr(s::AbstractString, traversal::Bool) = traversal ? " [$(stringify(s))]" : ""

encode_traversal(m::AbstractNode, idxs::Integer...) = stringify(_encode_traversal(m, idxs...))

function _encode_traversal(m, idxs...)
    !isempty(idxs) || return ""
    n = ith_child(m, idxs[1])
    return encode(idxs[1], descendants_n(m)) * _encode_traversal(n, idxs[2:end]...)
end


text:
using LinearAlgebraicRepresentation
Lar = LinearAlgebraicRepresentation
using SparseArrays;
using ViewerGL, LinearAlgebra
GL = ViewerGL

V = [0.6899599 0.502989 0.8326865 0.6457156 0.8586244 0.6716535 1.001351 0.8143801 0.5361978 0.3644763 0.5783999 0.4066785 0.5468975 0.375176 0.5890997 0.4173782;
 0.2417823 0.3845089 0.4884655 0.6311921 0.1912528 0.3339794 0.4379359 0.5806626 0.2867544 0.3289566 0.4588041 0.5010062 0.2854601 0.3276623 0.4575098 0.4997119;
 0.0757586 0.2444231 0.0252291 0.1938936 0.3054883 0.4741529 0.2549588 0.4236233 0.1554979 0.1661976 0.1542036 0.1649033 0.3323244 0.3430241 0.3310301 0.3417298]

CV = Array{Int64,1}[[1, 2, 3, 4, 5, 6, 7, 8], [9, 10, 11, 12, 13, 14, 15, 16]]

FV = Array{Int64,1}[[1, 2, 3, 4], [5, 6, 7, 8], [1, 2, 5, 6], [3, 4, 7, 8], [1, 3, 5, 7], [2, 4, 6, 8], [9, 10, 11, 12], [13, 14, 15, 16], [9, 10, 13, 14], [11, 12, 15, 16], [9, 11, 13, 15], [10, 12, 14, 16]]

EV = Array{Int64,1}[[1, 2], [3, 4], [5, 6], [7, 8], [1, 3], [2, 4], [5, 7], [6, 8], [1, 5], [2, 6], [3, 7], [4, 8], [9, 10], [11, 12], [13, 14], [15, 16], [9, 11], [10, 12], [13, 15], [14, 16], [9, 13], [10, 14], [11, 15], [12, 16]]

GL.VIEW([ GL.GLPol(V,CV, GL.COLORS[1]) ]);
cop_EV = Lar.coboundary_0(EV::Lar.Cells);
cop_FE = Lar.coboundary_1(V, FV::Lar.Cells, EV::Lar.Cells);
W = convert(Lar.Points, V');#cobordo 1
V, copEV, copFE, copCF = Lar.Arrangement.spatial_arrangement(W::Lar.Points, cop_EV::Lar.ChainOp, cop_FE::Lar.ChainOp);
triangulated_faces = Lar.triangulate2D(V, [copEV, copFE]);
FVs = convert(Array{Lar.Cells}, triangulated_faces);
VVV = convert(Lar.Points, V');
GL.VIEW( GL.GLExplode(VVV,FVs,1,1,1,99,1) );
EVs = Lar.FV2EVs(copEV, copFE); # polygonal face fragments
GL.VIEW( GL.GLExplode(VVV,EVs,1.5,1.5,1.5,99,1) );


text:
# Plotting examples
# =================
# ## Plotters
#  All plot functions in ExtendableGrids.jl have a `Plotter` keyword argument
#  which defaults to `nothing`.  This allows to pass a module as plotting backend
#  without creating a dependency. Fully supported are `PyPlot` and `GLMakie`.
#  `WGLMakie` and `CairoMakie` work in principle but in the moment don't deliver
#  all necessary functionality.
#
# Also supported is [`VTKView`](https://github.com/j-fu/VTKView.jl)  which is exeprimental and works only on linux
#
# Generally, grids and P1 FEM functions on grids can be plotted up to now.
#
# ## Grid plots
# Here, we define some sample grids for plotting purposes.
function grid1d(;n=50)
    X=collect(0:1/n:1)
    g=simplexgrid(X)
end

function grid2d(;n=20)
    X=collect(0:1/n:1)
    g=simplexgrid(X,X)
end

function grid3d(;n=15)
    X=collect(0:1/n:1)
    g=simplexgrid(X,X,X)
end
#
# Now, we can use the plot command of ExtendableGrids to plot grids
# Note the kwargs `xplane`, `yplane` and `zplane` which allow to control
# cutplanes which peel off some elements from the grid in 3d and allow to
# explore the inner triangulation.
#
# For Makie and VTKView, the cutplane values can be controlled interactively.
function plotting_grid3d(;Plotter=nothing)
    visualize(grid3d(), Plotter=Plotter, zplane=0.49)
end
# ![](plotting_grid3d.svg)

# ## Function plots
# Let us define some functions
function func1d(;n=50)
    g=grid1d(n=n)
    g,map(x->sinpi(2*x[1]),g)
end

function func2d(;n=20)
    g=grid2d(n=n)
    g,map((x,y)->sinpi(2*x)*sinpi(3.5*y),g)
end

function func3d(;n=15)
    g=grid3d(n=n)
    g, map((x,y,z)->sinpi(2*x)*sinpi(3.5*y)*sinpi(1.5*z),g)
end

#
# Plotting a function then goes as follows:
# `xplane`, `yplane` and `zplane` now define cut planes where
# the function projection is plotted as a heatmap.
# The additional `flevel` keyword argument allows
# to control an isolevel.
#
# For Makie and VTKView, the cutplane values and the flevel can be controlled interactively.
function plotting_func3d(;Plotter=nothing)
    g,f=func3d()
    visualize(g,f, Plotter=Plotter, zplane=0.49,xplane=0.49,flevel=0.25)
end
# ![](plotting_func3d.svg)

# ## Multiscene plots
# We can combine multiple plots into one scene according to
# some layout grid given by the layout parameter.
#
# The ',' key for GLMakie and the '*' key for VTKView allow to
# switch between gallery view (default) and focused view of only
# one subscene.
function plotting_multiscene(;Plotter=nothing)
    p=GridVisualizer(;Plotter=Plotter,layout=(2,3),clear=true,resolution=(800,500))

    # 3D must come first for Makie, see https://github.com/JuliaPlots/Makie.jl/issues/805
    visualize!(p[1,3],grid3d(),zplane=0.49,title="3D grid")
    visualize!(p[2,3],func3d()...,zplane=0.49,flevel=0.5,colormap=:bamako, title="3D grid function")
    visualize!(p[1,1],grid1d(), title="1D grid")
    visualize!(p[2,1],func1d()..., title="1D grid function")
    visualize!(p[1,2],grid2d(),title="2D grid")
    visualize!(p[2,2],func2d()...,colormap=:bamako,title="2D grid function")
    reveal(p)
end
# ![](plotting_multiscene.svg)


text:
int main() {
    double dA;
    double dB;
    dA=1.4e-3;
    dB=0.4e-3;
    if (dA-dB == 0.001) {
        printInt(99);
    }
    int iA;
    int iB;
    iA=342;
    iB=5123123;
    printInt(iA-iB);
    return 0;
}

max_stars_repo_name: carlostome/Javalette-Compiler

max_issues_repo_name: carlostome/Javalette-Compiler

max_forks_repo_name: carlostome/Javalette-Compiler


text:
using MAT
using jInv.Mesh
using EikonalInv
using jInv.LinearSolvers
using jInv.ForwardShare
using jInv.InverseSolve
using jInv.Utils
using ArgParse
using KrylovMethods

include("getDCResistivitySourcesAndRecAll.jl")

function parse_commandline()
	s = ArgParseSettings()

	@add_arg_table s begin

	"--n1"
		help = "number of cells in x1 direction "
		arg_type = Int
		default = 48
	"--n2"
		help = "number of cells in x2 direction "
		arg_type = Int
		default = 48
	"--n3"
		help = "number of cells in x3 direction "
		arg_type = Int
		default = 24

	"--nsources"
		help = "number of sources per worker"
		arg_type= Int
		default = 36
	"--out"
		help ="file to write into"
		default = "times.csv"
	end
	return parse_args(s)
end

function main()
	parsed_args = parse_commandline()

	nInv      = [128;128;64]
	m         = 1500 + 3000*rand(tuple(nInv...))
	# random conductivity model can be replaced by SEG model described in
	#
	# Aminzadeh, F., Brac, J., and Kunz, T., 1997. 3D Salt and Overthrust models. SEG/EAGE Modeling Series, No. 1: Distribution CD of Salt and
	# Overthrust models, SEG Book Series Tulsa, Oklahoma
	#
	# matfile   = matread("3Dseg12812864.mat")
	# m         = matfile["VELc"]

	n1 = parsed_args["n1"]
	n2 = parsed_args["n2"]
	n3 = parsed_args["n3"]
	out = parsed_args["out"]
	nsources = parsed_args["nsources"]

	# set number of threads for openbla to one for Eikonal
	@everywhere set_num_threads(1)

	domain  = [0;13.5;0;13.5;0;4.2]
	n       = [n1;n2;n3]
	M    = getRegularMesh(domain,n-1)

	# interpolate model to right resolution
	if any(M.n .!= nInv)
		MInv = getRegularMesh(domain,nInv)
		m   = getInterpolationMatrix(MInv,getRegularMesh(domain,n))*vec(m)
	end

	# get sources and receivers
	RCVfile = "DATA_SEG(64,64,32)_rcvMap.dat";
	SRCfile = "DATA_SEG(64,64,32)_srcMap.dat";
	srcNodeMap = readSrcRcvLocationFile(SRCfile,M);
	rcvNodeMap = readSrcRcvLocationFile(RCVfile,M);
	Q = generateSrcRcvProjOperators(M.n+1,srcNodeMap);
	Q = Q.*(1/(norm(M.h)^2));
	Q = Q[:,1:nsources]
	P = generateSrcRcvProjOperators(M.n+1,rcvNodeMap);
	nsrc = size(Q,2)

	@printf "problem size: n1=%d\tn2=%d\tn3=%d\tns=%d\n" n1 n2 n3 nsrc
	@printf "nworkers()=%d\n" nworkers()

	# setup pFor
	HO = false
	Qp = kron(Q,ones(1,nworkers())); # dulicate sources so that all workers solve same problems
	Qs = kron(Q[:,1],ones(1,nworkers()))
	(pForp,contDivEIK,SourcesSubIndEIK) = getEikonalInvParam(M,Qp,P,HO,nworkers())
	(pFors,contDivEIK,SourcesSubIndEIK) = getEikonalInvParam(M,Qs,P,HO,nworkers())

	# warm up
	wTime = @elapsed begin
	res = getData(vec(m),pFors,ones(length(pFors)),true)
	end
	@printf "%3.2f\t" wTime
	clear!(res[1])
	clear!(res[2])
	clear!(pFors)

	tTime = @elapsed begin
		m = vec(m);
	for k=1:5
		getData(m,pForp,ones(length(pForp)),true)
	end
	end
	tTime /= 5.0;
	@printf "%3.2f\n" tTime
	f = open(out,"a")
	str = @sprintf "%d,%d,%d,%d,%d,,%3.5f,%3.5f\n" nworkers() n1 n2 n3 nsrc wTime tTime
	write(f,str )
	close(f)
end
main()


text:
export sigmoid,∂sigmoid

function sigmoid(rᵢ)
  rₒ=1./(1+exp(-rᵢ))
end
function ∂sigmoid(rᵢ)
  y=sigmoid(rᵢ)
  rₒ=y.*(1-y)
end


text:
# f(x)=A*diagX
function FAmultdiagX(A::Array,X::Array)
    if size(A)==(1,1)
        return (A[1].*diagm(vec(X)),nothing)
    elseif size(X)==(1,1)
        return (A.*X[1],nothing)
    else
        return (A*diagm(vec(X)),nothing)
    end
end

function FAmultdiagX_inplace(value,auxvalue,A::Array,X::Array)
    if size(A)==(1,1)
        copy!(value,A[1].*diagm(vec(X)))
    elseif size(X)==(1,1)
        copy!(value,A.*X[1])
    else
        copy!(value,A*diagm(vec(X))) # BLAS has no dgmm routine
    end
end

function DAmultdiagX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X)
    if derivativeIDX==1
        if size(A)==(1,1)
            axpy!(1.0,[sum(X.*diag(grad_c))],grad_n) # sum_i x_i gradc_ii
        elseif size(X)==(1,1)
            axpy!(X[1],grad_c,grad_n) # x gradc_ii
        else
            axpy!(1.0,grad_c.*repmat(vec(X)',size(grad_c,1),1),grad_n)
#            axpy!(1.0,grad_c.*repmat(X,size(grad_c,1),1),grad_n) ## CHECK!
        end
    elseif derivativeIDX==2
        if size(A)==(1,1)
            axpy!(A[1],diag(grad_c),grad_n)
        elseif size(X)==(1,1)
            axpy!(1.0,[sum(A.*grad_c)],grad_n)
        else
            axpy!(1.0,sum(A.*grad_c,1),grad_n)
        end
    end
end

if PROC=="GPU"
    function FAmultdiagX_inplace(value::CudaArray,auxvalue,A::CudaArray,X::CudaArray)
        if size(A)==(1,1)
            diagm!(X,value)
            scale!(A,value)
        elseif size(X)==(1,1)
            copy!(value,A)
            scale!(X,value) # nb argument converse of Base.scale!
        else
            CUBLAS.dgmm!('R',A,vec(X),value)
        end
    end

    function DAmultdiagX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray)
        switched=false
        if size(X,2)>1
            switched=true
            X.dims=(size(X,2),1)
        end

        if derivativeIDX==1
            if size(A)==(1,1)
                tmp=CudaArray(Float64,(size(X,1),1))
                diag!(grad_c,tmp)
                vmult!(1.0,X,tmp,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
            elseif size(X)==(1,1)
                alphaaxpy!(1.0,X,grad_c,grad_n) # x gradc_ii
            else
                onc=CudaArray(Float64,(size(grad_c,1),1))
                fill!(onc,1.0);
                tmp=CudaArray(Float64,size(grad_c))
                gemm!('N','T',1.0,onc,X,0.0,tmp)
                vmultupdate!(1.0,grad_c,tmp,grad_n)
                free(tmp); free(onc)
                #                axpy!(1.0,grad_c.*repmat(X',size(grad_c,1),1),grad_n)
            end
        elseif derivativeIDX==2
            if size(A)==(1,1)
#                axpy!(A[1],diag(grad_c),grad_n)
                tmp=CudaArray(Float64,size(grad_n))
                diag!(grad_c,tmp)
                alphaaxpy!(1.0,A,tmp,grad_n)
                free(tmp)
            elseif size(X)==(1,1)
#                axpy!(1.0,[sum(A.*grad_c)],grad_n)
                tmp=CudaArray(Float64,size(A))
                vmult!(1.0,A,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
                free(tmp)
            else
#                axpy!(1.0,sum(A.*grad_c,1),grad_n)
                tmp=CudaArray(Float64,size(A))
                vmult!(1.0,A,grad_c,tmp)
                tmp2=CudaArray(Float64,(1,size(tmp,2)))
                sum!(tmp,1,0.0,tmp2)
                axpy!(1.0,vec(tmp2),grad_n)
                free(tmp); free(tmp2)
            end
        end

        if switched
            X.dims=(1,size(X,1))
        end
    end

    function DAmultdiagX(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray{Float32},X::CudaArray{Float32})
        switched=false
        if size(X,2)>1
            switched=true
            X.dims=(size(X,2),1)
        end
        if derivativeIDX==1
            if size(A)==(1,1)
                tmp=CudaArray(Float32,(size(X,1),1))
                diag!(grad_c,tmp)
                vmult!(1.0,X,tmp,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
            elseif size(X)==(1,1)
                alphaaxpy!(1.0,X,grad_c,grad_n) # x gradc_ii
            else
                onc=CudaArray(Float32,(size(grad_c,1),1))
                fill!(onc,1.0);
                tmp=CudaArray(Float32,size(grad_c))
                gemm!('N','T',1.0,vec(onc),X,0.0,tmp)
                vmultupdate!(1.0,grad_c,tmp,grad_n)
                free(tmp); free(onc)
                #                axpy!(1.0,grad_c.*repmat(X',size(grad_c,1),1),grad_n)
            end
        elseif derivativeIDX==2
            if size(A)==(1,1)
#                axpy!(A[1],diag(grad_c),grad_n)
                tmp=CudaArray(Float32,size(grad_n))
                diag!(grad_c,tmp)
                alphaaxpy!(1.0,A,tmp,grad_n)
                free(tmp)
            elseif size(X)==(1,1)
#                axpy!(1.0,[sum(A.*grad_c)],grad_n)
                tmp=CudaArray(Float32,size(A))
                vmult!(1.0,A,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n)
                free(tmp)
            else
#                axpy!(1.0,sum(A.*grad_c,1),grad_n)
                tmp=CudaArray(Float32,size(A))
                vmult!(1.0,A,grad_c,tmp)
                tmp2=CudaArray(Float32,(1,size(tmp,2)))
                sum!(tmp,1,0.0,tmp2)
                axpy!(1.0,vec(tmp2),grad_n)
                free(tmp); free(tmp2)
            end
        end

        if switched
            X.dims=(1,size(X,1))
        end

    end

end

Derivative[FAmultdiagX]=DAmultdiagX
export FAmultdiagX

Inplace[FAmultdiagX]=FAmultdiagX_inplace

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
"""
Return the longitudinal acceleration to come to a stop at the end of the current lane
follows idm style braking
"""
function stop_at_end(model::DriverModel, veh::Vehicle, roadway::Roadway)
    return stop_at_end(model, veh.state, roadway)
end
function stop_at_end(model::DriverModel, veh::VehicleState, roadway::Roadway)
    # run the IDM model, consider that the stop line is a fixed vehicle
    headway = get_dist_to_end(veh, roadway)
    return stop_at_dist(model, veh, headway)
end

"""
Return the longitudinal acceleration to come to a stop at the distance `headway` following idm
"""
function stop_at_dist(model::DriverModel, veh::VehicleState, headway::Float64)
    s_min = model.stop_delta
    acc = Inf
    v_des = 0.
    v = veh.v
    idm = model.navigator
    if headway > 0.0
        Δv = v_des - v
        s_des = s_min + v*idm.T - v*Δv / (2*sqrt(idm.a_max*idm.d_cmf))
        v_ratio = idm.v_des > 0.0 ? (v/idm.v_des) : 1.0
        acc = idm.a_max * (1.0 - v_ratio^idm.δ - (s_des/headway)^2)
    else # probably unnecessary
        acc = -model.navigator.d_max
    end

    acc = clamp(acc, -idm.d_max, idm.a_max)
    if abs(acc) < model.accel_tol
        acc = 0.
    end
    return acc
end

function stop_at_dist(model::DriverModel, veh::Vehicle, headway::Float64)
    return stop_at_dist(model, veh.state, headway)
end

"""
set the state stop! to true if veh is stopped at the intersection
"""
function update_stop!(model::DriverModel, veh::Vehicle, roadway::Roadway)
    dist_to_end = get_dist_to_end(veh, roadway)
    if veh.state.v ≈ 0. && isapprox(dist_to_end - model.stop_delta, 0, atol=0.5) # parameterize rtol?
        model.stop = true
    end
    model.stop = false
end

"""
Return the distance to the end of the lane
"""
function get_dist_to_end(veh::Vehicle, roadway::Roadway)
    return get_dist_to_end(veh.state, roadway)
end
function get_dist_to_end(veh::VehicleState, roadway::Roadway)
    lane = get_lane(roadway, veh)
    s_end = get_end(lane)
    dist_to_end = s_end - veh.posF.s
    # println("dist_to_end ", dist_to_end)
    # println("State ", veh.state.posF)
    return dist_to_end
end

max_stars_repo_name: MaximeBouton/AutomotivePOMDPs

max_issues_repo_name: MaximeBouton/AutomotivePOMDPs

max_forks_repo_name: MaximeBouton/AutomotivePOMDPs


text:
lines = readlines("input.txt")

function main()
    pairs = Dict('(' => ')', '[' => ']', '{' => '}', '<' => '>')
    revpairs = Dict(v => k for (k,v) = pairs)
    penalties = Dict(')' => 1, ']' => 2, '}' => 3, '>' => 4)
    Σs = []
    for line = lines
        Σ = 0
        stack = []
        for char = collect(line)
            if char in ['(', '[', '{', '<']
                pushfirst!(stack, char)
            else
                if popfirst!(stack) != revpairs[char]
                    stack = []
                    break
                end
            end
        end
        if length(stack) > 0
            for char = stack
                Σ = Σ * 5 + penalties[pairs[char]]
            end
            append!(Σs, Σ)
        end
    end
    println(sort(Σs)[1+length(Σs)÷2])
end

main()


text:
using Printf
using JLD2
using Plots
using Oceananigans
using Oceananigans.Grids
using Oceananigans.Advection
using Oceananigans.OutputWriters

# Useful for headless plotting with Plots.jl.
# See: https://github.com/jheinen/GR.jl/issues/278
ENV["GKSwstype"] = 100

function make_output_dirs()
    simulation_dir = joinpath(@__DIR__, "simulation")
    s = mkpath(simulation_dir)
    animation_path = joinpath(@__DIR__, "animation")
    a = mkpath(animation_path)
    return s, a
end
"""
    simulate_advection_diffusion(; Nx, Δt, stop_time, U, κ, initial_condition, output_interval, filename)

Run a 1D advection-diffusion simulation of a passive tracer with `Nx` grid points with a time step
`Δt` value until t = `stop_time` with constant velocity `U` and constant diffusivity `κ`.
The `inital_condition` should be a function of `x` describing the initial profile. Output will be saved to
JLD2 with `filename`.jld2 and animated every `output_interval` units of time.
"""
function simulate_advection_diffusion(; Nx, Δt, stop_time, U, κ, initial_condition, output_interval, filename)
    topo = (Periodic, Periodic, Periodic)
    domain = (x=(-1, 1), y=(0, 1), z=(0, 1))
    grid = RegularCartesianGrid(topology=topo, size=(Nx, 1, 1); domain...)
    simulation_dir, animation_dir = make_output_dirs()

    model = IncompressibleModel(
        architecture = CPU(),
                grid = grid,
         timestepper = :RungeKutta3,
           advection = WENO5(),
            buoyancy = nothing,
             tracers = :c,
             closure = IsotropicDiffusivity(κ=κ)
    )

    set!(model, u = U, c = (x, y, z) -> initial_condition(x))

    simulation = Simulation(model, Δt=Δt, stop_time=stop_time, iteration_interval=20)

    simulation.output_writers[:tracer] =
        JLD2OutputWriter(model, (c=model.tracers.c,), schedule=TimeInterval(output_interval),
                         dir=simulation_dir, prefix=filename, force=true)

    run!(simulation)

    file = jldopen(joinpath(simulation_dir, "$filename.jld2"))
    iterations = parse.(Int, keys(file["timeseries/t"]))

    anim = @animate for (i, iter) in enumerate(iterations)
        x = xnodes(model.tracers.c)
        t = file["timeseries/t/$iter"]
        c = file["timeseries/c/$iter"][:]

        plot(x, c, linewidth=2, title=@sprintf("t = %.3f", t),
             label="", xlabel="x", ylabel="Tracer", xlims=(-1, 1), ylims=(-1, 1))
    end

    gif(anim, joinpath(animation_dir, "$filename.gif"), fps=15)
    mp4(anim, joinpath(animation_dir, "$filename.mp4"), fps=15)
end

c₀_cosine(x) = cos(π*x)

simulate_advection_diffusion(Nx=32, Δt=0.01, stop_time=2, U=1, κ=0, initial_condition=c₀_cosine,
                             output_interval=0.05, filename="cosine_advection")

simulate_advection_diffusion(Nx=32, Δt=0.01, stop_time=2, U=0, κ=0.1, initial_condition=c₀_cosine,
                             output_interval=0.05, filename="cosine_diffusion")

simulate_advection_diffusion(Nx=32, Δt=0.01, stop_time=2, U=1, κ=0.1, initial_condition=c₀_cosine,
                             output_interval=0.05, filename="cosine_advection_diffusion")

max_stars_repo_name: drozzy/NeuralOceananigans

max_issues_repo_name: drozzy/NeuralOceananigans

max_forks_repo_name: drozzy/NeuralOceananigans


text:
# helpers useful for users
export deparameterize_type

@doc """Removes type parameters from types, e.g. Array{Int}->Array.

     It is often useful to make an associated type with this to match
     against methods which do not specialize on the type parameters.
     """ ->
deparameterize_type(A::DataType) = A.name.primary #eval(A.name.module, A.name.name)::DataType
deparameterize_type(A::TypeConstructor) = error("TypeConstructor not supported by deparameterize_type.")
  # could do A.body.name.primary but that would remove fixed parameters, like the 1 in Vector

# Internal helpers
##################
function eval_curmod(expr::Union{Symbol,Expr,QuoteNode})
    # evaluates a symbol or expression in the current module.
    # I.e. the one where the macro definition is.
    return eval(current_module(),expr)
end

# To iterate over code blocks dropping the line-number bits:
immutable Lines
    block::Expr
end
Base.start(lns::Lines) = 1
function Base.next(lns::Lines, nr)
    for i=nr:length(lns.block.args)
        if isa(lns.block.args[i], Expr) && !(lns.block.args[i].head==:line)
            return lns.block.args[i], i+1
        end
    end
    return -1
end
function Base.done(lns::Lines, nr)
    if next(lns::Lines, nr)==-1
        true
    else
        false
    end
end

@doc "Checks all elements of a collection are equal" ->
allequal(x) = reduce(&, [x[1]==xx for xx in x])

## Parsing
####
# translate a symbol in a array of expressions and/or symbols
function translate!(ex::Vector{Any}, di::Dict)
    for (i,e) in enumerate(ex)
        if isa(e, Symbol)
            ex[i] = get(di, e, e)
        else
            translate!(e.args, di)
        end
    end
    nothing
end

# expressions like :(I<:Int) are not parsed into TypeVar expressions
# but subtype comparisons.  This function translates this
function subt2tvar!(exs::Vector{Any})
    for (i,ex) in enumerate(exs)
        if isa(ex, Symbol)
            # do nothing
        elseif ex.head==:comparison
            exs[i] = Expr(:<:, ex.args[1], ex.args[3])
        else
            subt2tvar!(ex.args)
        end
    end
    nothing
end

function tvar2tvar!(exs::Vector{Any})
    # tranlates x<:Int -> TypeVar(:X, Int)
    for (i,ex) in enumerate(exs)
        if isa(ex, Symbol)
            # do nothing
        elseif ex.head==:<:
            var = ex.args[1]
            exs[i] = :(TypeVar(symbol($(string(var))), $(ex.args[2])))
        else
            tvar2tvar!(ex.args)
        end
    end
    nothing
end

# check whether a type is parameterized
isparameterized(t::DataType) = length(t.parameters)==0 ? false : true

function hasparameters(t::DataType)
    if isparameterized(t)
        if reduce(&, [map(typeof,t.parameters)...].==DataType)
            return true  # i.e. return true for Dict{Int,Int}
        else
            return false # but not for Dict{Int}
        end
    else
        return false
    end
end

# debugging
# import Traits: FakeMethod,  replace_concrete_tvars, get_tms_fms, isfitting
function get_tms_fms(Tr, gf::Function, ind=1)
    tm = collect(methods(Tr().methods[gf]))[1]
    fms = collect(methods(gf, NTuple{length(tm.sig),Any}))
    tm, fms
end


text:
"""
    Models for random hypergraphs
        * random model
        * k-uniform model
        * d-uniform model
        * preferential-attachment
"""

"""
    randomH(nVertices, nEdges)

Generate a *random* hypergraph without any structural constraints.

**The Algorithm**
Given two integer parameters *nVertices* and *nEdges* (the number of nodes and hyperedges, respectively),
the algorithm computes - for each hyperedge *he={1,...,m}* -
a random number *s ϵ [1, n]* (i.e. the hyperedge size).
Then, the algorithm selects uniformly at random *s* vertices from *V* to be added in *he*.
"""
function randomH(nVertices, nEdges)
    mx = Matrix{Union{Nothing,Bool}}(nothing, nVertices,nEdges)
    for e in 1:size(mx,2)
        nv = rand(1:size(mx,1))
        mx[sample(1:size(mx,1), nv;replace=false), e] .= true
    end

    h = Hypergraph(mx)
    if all(length.(h.v2he) .> 0)
        return h
    else
        return randomH(nVertices, nEdges)
    end
end

"""
    randomHkuniform(nVertices, nEdges, k)

Generates a *k*-uniform hypergraph, i.e. an hypergraph where each hyperedge has size *k*.

**The Algorithm**
The algorithm proceeds as the *randomH*, forcing the size of each hyperedge equal to *k*.
"""
function randomHkuniform(nVertices, nEdges, k)
    mx = Matrix{Union{Nothing,Bool}}(nothing, nVertices,nEdges)
    for e in 1:size(mx,2)
        nv = k#rand(1:size(mx,1))#rand(2:5)
        mx[sample(1:size(mx,1), nv;replace=false), e] .= true
    end

    h = Hypergraph(mx)
    if all(length.(h.v2he) .> 0)
        return h
    else
        return randomH(nVertices, nEdges)
    end
end

"""
    randomHduniform(nVertices, nEdges, d)

Generates a *d*-uniform hypergraph, where each node has degree *d*.

**The Algorithm**
The algorithm exploits the *k*-uniform approach described for the *randomHkuniform* method
to build a *d*-uniform hypergraph *H* having *nVertices* nodes and *nEdges* edges.
It returns the hypergraph H⃰ dual of *H*.
"""
function randomHduniform(nVertices, nEdges, d)
    mx = Matrix{Union{Nothing,Bool}}(nothing, nVertices,nEdges)

    for v in 1:size(mx,1)
        ne = d
        mx[v, sample(1:size(mx,2), ne;replace=false)] .= true
    end

    h = Hypergraph(mx)
    if all(length.(h.v2he) .> 0)
        return h
    else
        return randomH(nVertices, nEdges)
    end
end

"""
    randomHpreferential(nVertices, p)

Generate a hypergraph with a preferential attachment rule between nodes, as presented in
*Avin, C., Lotker, Z., and Peleg, D.Random preferential attachment hyper-graphs.Computer Science 23(2015).*

**The Algorithm**
The algorithm starts with a random graph with 5 nodes and 5 edges.
It iteratively adds a node or a edge, according to a given parameter *p*,
which defines the probability of creating a new node or a new hyperedge.

More in detail, the connections with the new node/hyperedge are generated according to
a preferential attachment policy.
"""
function randomHpreferential(nVertices, p)
    H₀ = randomH(10,10)
    H = H₀
    while nhv(H) < nVertices
        r = rand()
        y = rand(1:nhv(H))
        Y = Dict{Int,Bool}()
        if r < p
            #add a vertex
            v = SimpleHypergraphs.add_vertex!(H)
            push!(Y,v=>true)
            for v in nextNodes(H,y-1)
                push!(Y,v)
            end
        else
            for v in nextNodes(H,y)
                push!(Y,v)
            end
        end
            #add a hyperedge
            SimpleHypergraphs.add_hyperedge!(H, vertices=Y)
    end
    H
end

function nextNodes(h,size)
    nodes = Dict{Int, Bool}()

    ids = collect(1:nhv(h))
    degrees = length.(h.v2he)

    for s=1:size
        psum = collect(1:length(ids))

        psum[1] = degrees[ids[1]]
        for j=2:length(ids)
            psum[j] = psum[j-1] + degrees[ids[j]]
        end

        number = rand(1:psum[length(psum)])
        bucket = -1
        index=1
        for i=1:length(psum)
            if number <= psum[i]
                bucket = ids[i]
                index = i
                break
            end
        end

        push!(nodes, bucket=>true)
        deleteat!(ids, index)
    end
    nodes
end


text:
export assert_approx_equal

"""
    assert_approx_equal(x, y, ε_abs, ε_rel[, desc])

Assert that `x` is approximately equal to `y`.

Let `eps_z = eps_abs / eps_rel`. Call `x` and `y` small if
`abs(x) < eps_z` and `abs(y) < eps_z`, and call `x` and `y` large otherwise.  If this
function returns `True`, then it is guaranteed that
`abs(x - y) < 2 eps_rel max(abs(x), abs(y))` if `x` and `y` are large, and
`abs(x - y) < 2 eps_abs` if `x` and `y` are small.

# Arguments
- `x`: First object to compare.
- `y`: Second object to compare.
- `ε_abs`: Absolute tolerance.
- `ε_rel`: Relative tolerance.
- `desc`: Description of the comparison. Omit or set to `false` to have no description.
"""
function assert_approx_equal(x, y, ε_abs, ε_rel, desc)
    if abs(x - y) >= ε_abs + ε_rel * max(abs(x), abs(y))
        msg = "$(desc != false ? "\"$desc\": " : "")large deviation from reference:\n" *
              "  relative error:  $(@sprintf "%.3e" abs(x - y) / max(abs(x), abs(y)))\n" *
              "    tolerance:     $(@sprintf "%.3e" ε_rel)\n" *
              "  absolute error:  $(@sprintf "%.3e" abs(x - y))\n" *
              "    tolerance:     $(@sprintf "%.3e" ε_abs)\n"
        throw(ErrorException(msg))
    end
    return true
end
assert_approx_equal(x, y, ε_abs, ε_rel) = assert_approx_equal(x, y, ε_abs, ε_rel, false)