text:
module TestExamplesThreaded

using Test
using Trixi

include("test_trixi.jl")

# Start with a clean environment: remove Trixi output directory if it exists
outdir = "out"
Trixi.mpi_isroot() && isdir(outdir) && rm(outdir, recursive=true)

@testset "Threaded tests" begin
  @testset "TreeMesh" begin
    @trixi_testset "elixir_advection_restart.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_restart.jl"),
        l2   = [1.2148032444677485e-5],
        linf = [6.495644794757283e-5])
    end

    @trixi_testset "elixir_advection_amr_refine_twice.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_amr_refine_twice.jl"),
        l2   = [0.00019847333806230843],
        linf = [0.005591345460895569])
    end

    @trixi_testset "elixir_advection_amr_coarsen_twice.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_advection_amr_coarsen_twice.jl"),
        l2   = [0.00519897841357112],
        linf = [0.06272325552264647])
    end

    @trixi_testset "elixir_euler_source_terms_nonperiodic.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_source_terms_nonperiodic.jl"),
        l2   = [2.3652137675654753e-6, 2.1386731303685556e-6, 2.138673130413185e-6, 6.009920290578574e-6],
        linf = [1.4080448659026246e-5, 1.7581818010814487e-5, 1.758181801525538e-5, 5.9568540361709665e-5],
        rtol = 0.001)
    end

    @trixi_testset "elixir_euler_ec.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "tree_2d_dgsem", "elixir_euler_ec.jl"),
        l2   = [0.061751715597716854, 0.05018223615408711, 0.05018989446443463, 0.225871559730513],
        linf = [0.29347582879608825, 0.31081249232844693, 0.3107380389947736, 1.0540358049885143])
    end
  end

  @testset "StructuredMesh" begin
    @trixi_testset "elixir_advection_restart.jl with waving flag mesh" begin
      @test_trixi_include(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_advection_restart.jl"),
        l2   = [0.00017274040834067234],
        linf = [0.0015435741643734513],
        elixir_file="elixir_advection_waving_flag.jl",
        restart_file="restart_000041.h5")
    end

    @trixi_testset "elixir_mhd_ec.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "structured_2d_dgsem", "elixir_mhd_ec.jl"),
        l2   = [0.04937480811868297, 0.06117033019988596, 0.060998028674664716, 0.03155145889799417,
                0.2319175391388658, 0.02476283192966346, 0.024483244374818587, 0.035439957899127385,
                0.0016022148194667542],
        linf = [0.24749024430983746, 0.2990608279625713, 0.3966937932860247, 0.22265033744519683,
                0.9757376320946505, 0.12123736788315098, 0.12837436699267113, 0.17793825293524734,
                0.03460761690059514],
        tspan = (0.0, 0.3))
    end
  end

  @testset "UnstructuredMesh" begin
    @trixi_testset "elixir_ape_gauss_wall.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "unstructured_2d_dgsem", "elixir_ape_gauss_wall.jl"),
        l2   = [0.029331247985489625, 0.02934616721732521, 0.03803253571320854, 0.0,
                7.465932985352019e-16, 1.4931865970704038e-15, 1.4931865970704038e-15],
        linf = [0.3626825396196784, 0.3684490307932018, 0.8477478712580901, 0.0,
                8.881784197001252e-16, 1.7763568394002505e-15, 1.7763568394002505e-15],
        tspan = (0.0, 5.0))
    end
  end

  @testset "P4estMesh" begin
    @trixi_testset "elixir_euler_source_terms_nonperiodic_unstructured.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_euler_source_terms_nonperiodic_unstructured.jl"),
        l2   = [0.005689496253354319, 0.004522481295923261, 0.004522481295922983, 0.009971628336802528],
        linf = [0.05125433503504517, 0.05343803272241532, 0.053438032722404216, 0.09032097668196482])
    end

    @trixi_testset "elixir_eulergravity_convergence.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "p4est_2d_dgsem", "elixir_eulergravity_convergence.jl"),
        l2   = [0.00024871265138964204, 0.0003370077102132591, 0.0003370077102131964, 0.0007231525513793697],
        linf = [0.0015813032944647087, 0.0020494288423820173, 0.0020494288423824614, 0.004793821195083758],
        tspan = (0.0, 0.1))
    end
  end

  @testset "DGMulti" begin
    @trixi_testset "elixir_euler_weakform.jl (SBP, EC)" begin
      @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_euler_weakform.jl"),
        cells_per_dimension = (4, 4),
        volume_integral = VolumeIntegralFluxDifferencing(flux_ranocha),
        surface_integral = SurfaceIntegralWeakForm(flux_ranocha),
        approximation_type = SBP(),
        l2 = [0.01285864081726596, 0.010650165503847099, 0.01065016550381281, 0.026286162111579015],
        linf = [0.037333313274372504, 0.05308320130762212, 0.05308320130841948, 0.13378665881805185],
      )
    end

    @trixi_testset "elixir_euler_triangulate_pkg_mesh.jl" begin
      @test_trixi_include(joinpath(examples_dir(), "dgmulti_2d", "elixir_euler_triangulate_pkg_mesh.jl"),
        l2 = [4.664661209491976e-6, 3.7033509525940745e-6, 4.794877426562555e-6, 1.2682723101532175e-5],
        linf = [2.5099852761334418e-5, 2.2683684021362893e-5, 2.6180448559287584e-5, 5.5752932611508044e-5]
      )
    end
  end
end

# Clean up afterwards: delete Trixi output directory
Trixi.mpi_isroot() && isdir(outdir) && @test_nowarn rm(outdir, recursive=true)

end # module


text:
julia> include("minesweeper.jl")

Welcome to Minesweeper

Enter the gridsize x y:
6 4
Found 0 of 5 mines.

     1  2  3  4  5  6
    __________________
 1 | .  .  .  .  .  .
 2 | .  .  .  .  .  .
 3 | .  .  .  .  .  .
 4 | .  .  .  .  .  .

What do you do? ("o x y" to reveal a field; "m x y" to toggle a mine; close)
o 1 1
________________________________________________________________________________

Found 0 of 5 mines.

     1  2  3  4  5  6
    __________________
 1 | 1  .  .  .  .  .
 2 | .  .  .  .  .  .
 3 | .  .  .  .  .  .
 4 | .  .  .  .  .  .

What do you do? ("o x y" to reveal a field; "m x y" to toggle a mine; close)
o 1 2
________________________________________________________________________________

Found 0 of 5 mines.

     1  2  3  4  5  6
    __________________
 1 | 1  .  .  .  .  .
 2 | 2  .  .  .  .  .
 3 | .  .  .  .  .  .
 4 | .  .  .  .  .  .

What do you do? ("o x y" to reveal a field; "m x y" to toggle a mine; close)
o 2 1
________________________________________________________________________________

Found 0 of 5 mines.

     1  2  3  4  5  6
    __________________
 1 | 1  2  .  .  .  .
 2 | 2  .  .  .  .  .
 3 | .  .  .  .  .  .
 4 | .  .  .  .  .  .

What do you do? ("o x y" to reveal a field; "m x y" to toggle a mine; close)
o 6 4
________________________________________________________________________________

Found 0 of 5 mines.

     1  2  3  4  5  6
    __________________
 1 | 1  2  .  .  .  .
 2 | 2  .  2  1  1  1
 3 | .  .  2
 4 | .  .  1

What do you do? ("o x y" to reveal a field; "m x y" to toggle a mine; close)
m 2 2
________________________________________________________________________________

Found 1 of 5 mines.

     1  2  3  4  5  6
    __________________
 1 | 1  2  .  .  .  .
 2 | 2  ?  2  1  1  1
 3 | .  .  2
 4 | .  .  1

What do you do? ("o x y" to reveal a field; "m x y" to toggle a mine; close)
close
ERROR: You closed the game.

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
# More test utilities. Can't be included in KernelFunctions because they introduce a number
# of additional deps that we don't want to have in the main package.

# Check parameters of kernels

function test_params(kernel, reference)
    params_kernel = params(kernel)
    params_reference = params(reference)

    @test length(params_kernel) == length(params_reference)
    @test all(p == q for (p, q) in zip(params_kernel, params_reference))
end

# AD utilities

const FDM = FiniteDifferences.central_fdm(5, 1)

gradient(f, s::Symbol, args) = gradient(f, Val(s), args)

function gradient(f, ::Val{:Zygote}, args)
    g = first(Zygote.gradient(f, args))
    if isnothing(g)
        if args isa AbstractArray{<:Real}
            return zeros(size(args)) # To respect the same output as other ADs
        else
            return zeros.(size.(args))
        end
    else
        return g
    end
end

function gradient(f, ::Val{:ForwardDiff}, args)
    ForwardDiff.gradient(f, args)
end

function gradient(f, ::Val{:ReverseDiff}, args)
    ReverseDiff.gradient(f, args)
end

function gradient(f, ::Val{:FiniteDiff}, args)
    first(FiniteDifferences.grad(FDM, f, args))
end

function compare_gradient(f, AD::Symbol, args)
    grad_AD = gradient(f, AD, args)
    grad_FD = gradient(f, :FiniteDiff, args)
    @test grad_AD ≈ grad_FD atol=1e-8 rtol=1e-5
end

testfunction(k, A, B, dim) = sum(kernelmatrix(k, A, B, obsdim = dim))
testfunction(k, A, dim) = sum(kernelmatrix(k, A, obsdim = dim))

function test_ADs(kernelfunction, args = nothing; ADs = [:Zygote, :ForwardDiff, :ReverseDiff], dims = [3, 3])
    test_fd = test_FiniteDiff(kernelfunction, args, dims)
    if !test_fd.anynonpass
        for AD in ADs
            test_AD(AD, kernelfunction, args, dims)
        end
    end
end

function test_FiniteDiff(kernelfunction, args = nothing, dims = [3, 3])
    # Init arguments :
    k = if args === nothing
        kernelfunction()
    else
        kernelfunction(args)
    end
    rng = MersenneTwister(42)
    @testset "FiniteDifferences" begin
        if k isa SimpleKernel
            for d in log.([eps(), rand(rng)])
                @test_nowarn gradient(:FiniteDiff, [d]) do x
                    kappa(k, exp(first(x)))
                end
            end
        end
        ## Testing Kernel Functions
        x = rand(rng, dims[1])
        y = rand(rng, dims[1])
        @test_nowarn gradient(:FiniteDiff, x) do x
                k(x, y)
            end
        if !(args === nothing)
            @test_nowarn gradient(:FiniteDiff, args) do p
                kernelfunction(p)(x, y)
            end
        end
        ## Testing Kernel Matrices
        A = rand(rng, dims...)
        B = rand(rng, dims...)
        for dim in 1:2
            @test_nowarn gradient(:FiniteDiff, A) do a
                testfunction(k, a, dim)
            end
            @test_nowarn gradient(:FiniteDiff , A) do a
                testfunction(k, a, B, dim)
            end
            @test_nowarn gradient(:FiniteDiff, B) do b
                testfunction(k, A, b, dim)
            end
            if !(args === nothing)
                @test_nowarn gradient(:FiniteDiff, args) do p
                    testfunction(kernelfunction(p), A, B, dim)
                end
            end
        end
    end
end

function test_AD(AD::Symbol, kernelfunction, args = nothing, dims = [3, 3])
    @testset "$(AD)" begin
        # Test kappa function
        k = if args === nothing
            kernelfunction()
        else
            kernelfunction(args)
        end
        rng = MersenneTwister(42)
        if k isa SimpleKernel
            for d in log.([eps(), rand(rng)])
                compare_gradient(AD, [d]) do x
                    kappa(k, exp(x[1]))
                end
            end
        end
        # Testing kernel evaluations
        x = rand(rng, dims[1])
        y = rand(rng, dims[1])
        compare_gradient(AD, x) do x
            k(x, y)
        end
        compare_gradient(AD, y) do y
            k(x, y)
        end
        if !(args === nothing)
            compare_gradient(AD, args) do p
                kernelfunction(p)(x,y)
            end
        end
        # Testing kernel matrices
        A = rand(rng, dims...)
        B = rand(rng, dims...)
        for dim in 1:2
            compare_gradient(AD, A) do a
                testfunction(k, a, dim)
            end
            compare_gradient(AD, A) do a
                testfunction(k, a, B, dim)
            end
            compare_gradient(AD, B) do b
                testfunction(k, A, b, dim)
            end
            if !(args === nothing)
                compare_gradient(AD, args) do p
                    testfunction(kernelfunction(p), A, dim)
                end
            end
        end
    end
end


text:
export vectorize

function vectorize(str)
  str=replace(str,".*","*")
  str=replace(str,"./","/")
  str=replace(str,".^","^")
  str=replace(str,"*",".*")
  str=replace(str,"/","./")
  str=replace(str,"^",".^")
  return str
end


text:
#
# Normal Sampler
#
# parameters:
# 	(1) mu
#	(2) sigma
#

type NormalSampler <: ContinuousDistributionSampler
	paramranges::Vector{Tuple{Float64,Float64}}
	distribution::Normal
	function NormalSampler(params=Float64[])
		# We restrict the range of the real-values parameters to less than the entire float range
		# in order to avoid overflow errors in Distributions.Normal
		# We use sqrt of the bound for sigma on the basis that it is the calculation of the variance
		# that can cause an overflow
		s = new( [(-realmax(Float64), realmax(Float64)), (0, realmax(Float64))] )
		setparams!(s, isempty(params) ? [0.0, 1.0] : params)
		s
	end
end

function setparams!(s::NormalSampler, params)
	checkparamranges(s, params)
	μ, σ = params[1], max(nextfloat(0.0), params[2]) # sigma cannot be zero
	s.distribution = Normal(μ, σ)
end

getparams(s::NormalSampler) = [s.distribution.μ, s.distribution.σ]
# TODO should sample actually sample Inf?

estimateparams!(s::NormalSampler, traces) = estimateparams!(s, traces, 2)


text:
compressed_model(model::MetNet; sparsity_th = 0.66) =
    MetNet(compressed_copy(struct_to_dict(model);
        sparsity_th = sparsity_th))
uncompressed_model(model::MetNet) =
    MetNet(uncompressed_copy(struct_to_dict(model)))


text:
export return_calculate,
       return_accumulate,
       return_annualized;

function return_calculate(prices::DataFrame;
                method="DISCRETE",
                dateColumn="date")

        vars = names(prices)
        nVars = length(vars)
        vars = Symbol.(vars[vars.!=dateColumn])
        if nVars == length(vars)
                throw(ArgumentError(string("dateColumn: ", dateColumn, " not in DataFrame: ",vars)))
        end
        nVars = nVars-1

        p = Matrix(prices[!,vars])
        n = size(p,1)
        m = size(p,2)
        p2 = Array{Float64,2}(undef,n-1,m)
        # p = p[1:(n-1),:]

        @inbounds @avx for i ∈ 1:(n-1), j ∈ 1:m
            p2[i,j] = p[i+1,j] / p[i,j]
        end

        if uppercase(method) == "DISCRETE"
                p2 = p2 .- 1.0
        elseif uppercase(method) == "LOG"
                p2 = log.(p2)
        else
                throw(ArgumentError(string("method: ", method, " must be in (\"LOG\",\"DISCRETE\")")))
        end

        dates = prices[2:n,Symbol(dateColumn)]
        out = DataFrame(Symbol(dateColumn)=>dates)
        for i in 1:nVars
                out[!,vars[i]] = p2[:,i]
        end
        return(out)
end

function return_accumulate(returns::DataFrame;
                method="DISCRETE",
                toFreq="MONTH",
                dateColumn="date"
        )

        if uppercase(toFreq) in ("DAY", "DAILY")
                # println("Daily")
                byVars = [:year, :month, :day]
                _selector = "YMD"
        elseif uppercase(toFreq) in ("MONTH","MONTHLY","MTH")
                # println("Monthly")
                byVars = [:year, :month]
                _selector = "YM"
        elseif uppercase(toFreq) in ("QTR","QUARTERLY","QUARTER")
                # println("Quarterly")
                byVars = [:year, :qtr]
                _selector = "YQ"
        elseif uppercase(toFreq) in ("YR","YEARLY","YEAR")
                # println("Yearly")
                byVars = [:year]
                _selector = "Y"
        else
                throw(ArgumentError(string("toFreq: ", toFreq, " must be one of ",
                        ("DAY", "DAILY")," / ",
                        ("MONTH","MONTHLY","MTH")," / ",
                        ("QTR","QUARTERLY","QUARTER")," / ",
                        ("YR","YEARLY","YEAR")
                        )))
        end

        # println(byVars)

        vars = names(returns)
        nVars = length(vars)
        vars = Symbol.(vars[vars.!=dateColumn])
        if nVars == length(vars)
                throw(ArgumentError(string("dateColumn: ", dateColumn, " not in DataFrame: ",vars)))
        end
        nVars = nVars-1
        n = size(returns,1)
        dc = Symbol(dateColumn)

        out = returns[:, vars]
        dates = returns[:,dc]
        out[!,:year] =  Dates.Year.(dates)
        if findfirst("M", _selector) != nothing
                out[!,:month] = Dates.Month.(dates)
        end
        if findfirst("Q", _selector) != nothing
                out[!,:qtr] = Dates.quarterofyear.(dates)
        end
        if findfirst("D", _selector) != nothing
                out[!,:day] = Dates.Day.(dates)
        end

        if uppercase(method) == "DISCRETE"
                out = rename(
                        combine(
                            groupby(out,byVars),
                            [vars[i] => aggDiscreteReturns for i in 1:size(vars,1)]
                        ),
                        Dict(Symbol(String(string(vars[i],"_aggDiscreteReturns"))) => vars[i] for i in 1:size(vars,1))
                    )
        elseif uppercase(method) == "LOG"
                out = rename(
                        combine(
                            groupby(out,byVars),
                            [vars[i] => sum for i in 1:size(vars)[1]]
                        ),
                        Dict(Symbol(String(string(vars[i],"_sum"))) => vars[i] for i in 1:size(vars,1))
                    )
        else
                throw(ArgumentError(string("method: ", method, " must be in (\"LOG\",\"DISCRETE\")")))
        end

        if "YM" == _selector
                out[!, dc] = Date.(out.year,out.month).+ Dates.Month(1) .- Dates.Day(1)
        elseif "YQ"  == _selector
                out[!, dc] = Date.(out.year) + Dates.Month.( out.qtr .* 3) .- Dates.Day(1)
        elseif "YMD" == _selector
                out[!, Symbol(dateColumn)] = Date.(out.year,out.month,out.day)
        else
                out[!, Symbol(dateColumn)] = Date.(out.year) .+ Dates.Year(1) .- Dates.Day(1)
        end

        return(out[!, vcat(dc, vars)])
end

function return_annualized(returns::DataFrame; scale::Int=1,
                            method="DISCRETE",
                            dateColumn="date")

    vars = names(returns)
    nVars = length(vars)
    vars = Symbol.(vars[vars.!=dateColumn])
    if nVars == length(vars)
        throw(ArgumentError(string("dateColumn: ", dateColumn, " not in DataFrame: ",vars)))
    end
    nVars = nVars-1

    out = DataFrame(:_stat_=>"Annualized Return")

    s = convert(Float64,scale)

    if uppercase(method) == "DISCRETE"
            temp = geo_mean(returns,dateColumn=dateColumn)
#   TODO: Test if Threads.@treads is faster for a large number of variables...
            for i in 1:nVars
                out[!,vars[i]] .= (1.0 .+ temp[!,vars[i]]).^s .- 1.0
            end
    elseif uppercase(method) == "LOG"
        for i in 1:nVars
            out[!,vars[i]] .= (mean(returns[!,vars[i]])).*s
        end
    else
            throw(ArgumentError(string("method: ", method, " must be in (\"LOG\",\"DISCRETE\")")))
    end

    return(out)
end


text:
export snapshot, loadsnapshot

# macro snapshot(a...)
#     quote
        # dictopen("/tmp/snapshot.dictfile") do df
        #     for i = 1:length(a)
        #         #@show a[1] $a[i]
        #         df[string(a[i])] = [i]
        #     end
        # end
    # end
# end

function snapshot(a...)
    dictopen("/tmp/snapshot.dictfile","w") do df
        if length(a) == 1
            a = (:Dictfile_snapshot, a[1])
        end
        for i = 1:2:length(a)
            df[a[i]] = a[i+1]
        end
    end
end

loadsnapshot(a...) = dictopen("/tmp/snapshot.dictfile","r") do df
    if length(a) == 0
        r = df[]
        if collect(keys(r)) == [:Dictfile_snapshot]
            return df[:Dictfile_snapshot]
        end
        r
    else
        if length(a) == 1
            df[a[1]]
        else
            Dict(zip(a, [df[x] for x in a]))
        end
    end
end


text:
using Plots, Plots.PlotMeasures
using CSV, DataFrames
using Statistics
# %% Read in the fitting results of the Tensor-to-Scalar Ratio (r) of SO-SAT only and SO-SAT+FYST/Prime-Cam
milca_so = CSV.read("sosat/milca_results.csv", DataFrame)
println("MILCA:")
println(
    "- SO-SAT only: r = ",
    mean(milca_so.r_wo_FGmarg),
    " ± ",
    std(milca_so.r_wo_FGmarg),
)
milca = CSV.read("sosatfyst/milca_results.csv", DataFrame)
println(
    "- SO-SAT+FYST: r = ",
    mean(milca.r_wo_FGmarg),
    " ± ",
    std(milca.r_wo_FGmarg),
)
# %% Scatter Plot of r from SO-SAT only and SO-SAT+FYST/Prime-Cam
# Save to "figure7a.pdf"
p = scatter(
    milca.r_wo_FGmarg,
    milca_so.r_wo_FGmarg,
    xlab = "Tensor-to-scalar ratio, SO-SAT + FYST",
    ylab = "Tensor-to-scalar ratio, SO-SAT only",
    lab = "",
    legend = :best,
    xlims = [-0.005, 0.012],
    ylims = [-0.005, 0.012],
    aspect_ratio = :equal,
    markersize = 3,
    legendfontsize = 10,
    tickfontsize = 15,
    labelfontsize = 14,
    #left_margin = -100px,
)
p = plot!(
    milca.r_wo_FGmarg,
    milca.r_wo_FGmarg,
    lw = 2,
    lab = "r(SO-SAT only) = r(SO-SAT+FYST)",
)
p = plot!(
    milca.r_wo_FGmarg,
    milca.r_wo_FGmarg .+ 1.1e-3,
    lw = 2,
    lab = "r(SO-SAT only) = r(SO-SAT+FYST) + 0.0011",
)
savefig("figure7a.pdf")
display(p)
# %% Histogram. Save to "figure7b.pdf"
p = histogram(
    milca_so.r_wo_FGmarg,
    yticks = nothing,
    xlims = [-0.005, 0.012],
    xlab = "Tensor-to-scalar ratio, r",
    lab = "SO-SAT only",
    α = 0.7,
    legend = :best,
    legendfontsize = 15,
    tickfontsize = 13,
    labelfontsize = 20,
)
p = histogram!(milca.r_wo_FGmarg, lab = "SO-SAT+FYST", α = 0.7)
p = vline!([0], c = :grey, lw = 10, lab = "")
savefig("figure7b.pdf")
display(p)


text:
import Base: copy!
import Base: splice!
import Base: append!
import Base: push!
import Base: resize!
import Base: setindex!
import Base: getindex
import Base: map
import Base: length
import Base: eltype
import Base: endof
import Base: ndims
import Base: size
import Base: start
import Base: next
import Base: done

import GeometryTypes.SimpleRectangle

abstract type GPUArray{T, NDim} <: AbstractArray{T, NDim} end

#=
immutable GPUArray{T, NDim, GPUBuff <: GPUBuffer} <: DenseArray{T, NDim}
    buff::GPUBuff{T, NDim}
    size::NTuple{Int, NDim}
end

immutable BufferedGPUArray{GPUArr <: GPUArray}
    buff::GPUBuff{T, NDim}
    ram::Array{T, NDim}
end
=#

length(A::GPUArray)                                     = prod(size(A))
eltype(b::GPUArray{T, NDim}) where {T, NDim} = T
endof(A::GPUArray)                                      = length(A)
ndims(A::GPUArray{T, NDim}) where {T, NDim} = NDim
size(A::GPUArray)                                       = A.size
size(A::GPUArray, i::Integer)                           = i <= ndims(A) ? A.size[i] : 1

function checkdimensions(value::Array, ranges::Union{Integer, UnitRange}...)
    array_size   = size(value)
    indexes_size = map(length, ranges)

    (array_size != indexes_size) && throw(DimensionMismatch("asigning a $array_size to a $(indexes_size) location"))
    true
end
function to_range(index)
    map(index) do val
        isa(val, Integer) && return val:val
        isa(val, Range) && return val
        error("Indexing only defined for integers or ranges. Found: $val")
    end
end
setindex!(A::GPUArray{T, N}, value::Union{T, Array{T, N}}) where {T, N} = (A[1] = value)

function setindex!(A::GPUArray{T, N}, value, indexes...) where {T, N}
    ranges = to_range(indexes)
    v = isa(value, T) ? [value] : convert(Array{T,N}, value)
    setindex!(A, v, ranges...)
    nothing
end

setindex!(A::GPUArray{T, 2}, value::Vector{T}, i::Integer, range::UnitRange) where {T} =
   (A[i, range] = reshape(value, (length(value),1)))

function setindex!(A::GPUArray{T, N}, value::Array{T, N}, ranges::UnitRange...) where {T, N}
    checkbounds(A, ranges...)
    checkdimensions(value, ranges...)
    gpu_setindex!(A, value, ranges...)
    nothing
end

function update!(A::GPUArray{T, N}, value::Array{T, N}) where {T, N}
    if length(A) != length(value)
        if isa(A, GLBuffer)
            resize!(A, length(value))
        elseif isa(A, Texture) && ndims(A) == 2
            resize_nocopy!(A, size(value))
        else
            error("Dynamic resizing not implemented for $(typeof(A))")
        end
    end
    dims = map(x->1:x, size(A))
    A[dims...] = value
    nothing
end

function getindex(A::GPUArray{T, N}, i::Int) where {T, N}
    checkbounds(A, i)
    gpu_getindex(A, i:i)[1] # not as bad as its looks, as so far gpu data must be loaded into an array anyways
end
function getindex(A::GPUArray{T, N}, ranges::UnitRange...) where {T, N}
    checkbounds(A, ranges...)
    gpu_getindex(A, ranges...)
end

function getindex(A::GPUArray{T, N}, rect::SimpleRectangle) where {T, N}
    A[rect.x+1:rect.x+rect.w, rect.y+1:rect.y+rect.h]
end
function setindex!(A::GPUArray{T, N}, value::Array{T, N}, rect::SimpleRectangle) where {T, N}
    A[rect.x+1:rect.x+rect.w, rect.y+1:rect.y+rect.h] = value
end

mutable struct GPUVector{T} <: GPUArray{T, 1}
    buffer
    size
    real_length
end
GPUVector(x::GPUArray) = GPUVector{eltype(x)}(x, size(x), length(x))

function update!(A::GPUVector{T}, value::Vector{T}) where T
    if isa(A, GLBuffer) && (length(A) != length(value))
        resize!(A, length(value))
    end
    dims = map(x->1:x, size(A))
    A.buffer[dims...] = value
    nothing
end

length(v::GPUVector)            = prod(size(v))
size(v::GPUVector)              = v.size
size(v::GPUVector, i::Integer)  = v.size[i]
ndims(::GPUVector)              = 1
eltype(::GPUVector{T}) where {T}       = T
endof(A::GPUVector)             = length(A)

start(b::GPUVector)             = start(b.buffer)
next(b::GPUVector, state)       = next(b.buffer, state)
done(b::GPUVector, state)       = done(b.buffer, state)

gpu_data(A::GPUVector)          = A.buffer[1:length(A)]

getindex(v::GPUVector, index::Int)       = v.buffer[index]
getindex(v::GPUVector, index::UnitRange) = v.buffer[index]
setindex!(v::GPUVector{T}, value::T, index::Int) where {T} = v.buffer[index] = value
setindex!(v::GPUVector{T}, value::T, index::UnitRange) where {T} = v.buffer[index] = value

function grow_dimensions(real_length::Int, _size::Int, additonal_size::Int, growfactor::Real=1.5)
    new_dim = round(Int, real_length*growfactor)
    return max(new_dim, additonal_size+_size)
end
function Base.push!(v::GPUVector{T}, x::Vector{T}) where T
    lv, lx = length(v), length(x)
    if (v.real_length < lv+lx)
        resize!(v.buffer, grow_dimensions(v.real_length, lv, lx))
    end
    v.buffer[lv+1:(lv+lx)] = x
    v.real_length          = length(v.buffer)
    v.size                 = (lv+lx,)
    v
end
push!(v::GPUVector{T}, x::T) where {T} = push!(v, [x])
push!(v::GPUVector{T}, x::T...) where {T} = push!(v, [x...])
append!(v::GPUVector{T}, x::Vector{T}) where {T} = push!(v, x)

resize!(A::GPUArray{T, NDim}, dims::Int...) where {T, NDim} = resize!(A, dims)
function resize!(A::GPUArray{T, NDim}, newdims::NTuple{NDim, Int}) where {T, NDim}
    newdims == size(A) && return A
    gpu_resize!(A, newdims)
    A
end

function resize!(v::GPUVector, newlength::Int)
    if v.real_length >= newlength # is still big enough
        v.size = (max(0, newlength),)
        return v
    end
    resize!(v.buffer, grow_dimensions(v.real_length, length(v),  newlength-length(v)))
    v.size        = (newlength,)
    v.real_length = length(v.buffer)
end
function grow_at(v::GPUVector, index::Int, amount::Int)
    resize!(v, length(v)+amount)
    copy!(v, index, v, index+amount, amount)
end

function splice!(v::GPUVector{T}, index::UnitRange, x::Vector=T[]) where T
    lenv = length(v)
    elements_to_grow = length(x)-length(index) # -1
    buffer           = similar(v.buffer, length(v)+elements_to_grow)
    copy!(v.buffer, 1, buffer, 1, first(index)-1) # copy first half
    copy!(v.buffer, last(index)+1, buffer, first(index)+length(x), lenv-last(index)) # shift second half
    v.buffer      = buffer
    v.real_length = length(buffer)
    v.size        = (v.real_length,)
    copy!(x, 1, buffer, first(index), length(x)) # copy contents of insertion vector
    nothing
end
splice!(v::GPUVector{T}, index::Int, x::T) where {T} = v[index] = x
splice!(v::GPUVector{T}, index::Int, x::Vector=T[]) where {T} = splice!(v, index:index, map(T, x))

copy!(a::GPUVector, a_offset::Int, b::Vector, b_offset::Int, amount::Int)   = copy!(a.buffer, a_offset, b,        b_offset, amount)
copy!(a::GPUVector, a_offset::Int, b::GPUVector, b_offset::Int, amount::Int)= copy!(a.buffer, a_offset, b.buffer, b_offset, amount)

copy!(a::GPUArray, a_offset::Int, b::Vector,   b_offset::Int, amount::Int) = _copy!(a, a_offset, b, b_offset, amount)
copy!(a::Vector,   a_offset::Int, b::GPUArray, b_offset::Int, amount::Int) = _copy!(a, a_offset, b, b_offset, amount)
copy!(a::GPUArray, a_offset::Int, b::GPUArray, b_offset::Int, amount::Int) = _copy!(a, a_offset, b, b_offset, amount)

#don't overwrite Base.copy! with a::Vector, b::Vector
function _copy!(a::Union{Vector, GPUArray}, a_offset::Int, b::Union{Vector, GPUArray}, b_offset::Int, amount::Int)
    (amount <= 0) && return nothing
    @assert a_offset > 0 && (a_offset-1) + amount <= length(a) "a_offset $a_offset, amount $amount, lengtha $(length(a))"
    @assert b_offset > 0 && (b_offset-1) + amount <= length(b) "b_offset $b_offset, amount $amount, lengthb $(length(b))"
    unsafe_copy!(a, a_offset, b, b_offset, amount)
    return nothing
end

# Interface:
gpu_data(t)      = error("gpu_data not implemented for: $(typeof(t)). This happens, when you call data on an array, without implementing the GPUArray interface")
gpu_resize!(t)   = error("gpu_resize! not implemented for: $(typeof(t)). This happens, when you call resize! on an array, without implementing the GPUArray interface")
gpu_getindex(t)  = error("gpu_getindex not implemented for: $(typeof(t)). This happens, when you call getindex on an array, without implementing the GPUArray interface")
gpu_setindex!(t) = error("gpu_setindex! not implemented for: $(typeof(t)). This happens, when you call setindex! on an array, without implementing the GPUArray interface")
max_dim(t)       = error("max_dim not implemented for: $(typeof(t)). This happens, when you call setindex! on an array, without implementing the GPUArray interface")

function (::Type{T})(x::Signal) where T <: GPUArray
    gpu_mem = T(Reactive.value(x))
    preserve(const_lift(update!, gpu_mem, x))
    gpu_mem
end

const BaseSerializer = if isdefined(Base, :AbstractSerializer)
    Base.AbstractSerializer
elseif isdefined(Base, :SerializationState)
    Base.SerializationState
else
    error("No Serialization type found. Probably unsupported Julia version")
end

function Base.serialize(s::BaseSerializer, t::T) where T<:GPUArray
    Base.serialize_type(s, T)
    serialize(s, Array(t))
end
function Base.deserialize(s::BaseSerializer, ::Type{T}) where T<:GPUArray
    A = deserialize(s)
    T(A)
end

export data
export resize
export GPUArray
export GPUVector

export update!

export gpu_data
export gpu_resize!
export gpu_getindex
export gpu_setindex!
export max_dim


text:
using Gtk.ShortNames
import GR

function paint(w)
    ctx = Gtk.getgc(w)

    h = Gtk.height(w)
    w = Gtk.width(w)

    ENV["GKS_WSTYPE"] = "142"
    ENV["GKSconid"] = @sprintf("%lu", UInt64(ctx.ptr))

    Gtk.select_font_face(ctx, "Sans",
                         Cairo.FONT_SLANT_NORMAL, Cairo.FONT_WEIGHT_NORMAL);
    Gtk.move_to(ctx, 15, 15)
    Gtk.set_font_size(ctx, 14)
    Gtk.show_text(ctx, "Contour Plot using Gtk ...")

    srand(0)
    xd = -2 + 4 * rand(100)
    yd = -2 + 4 * rand(100)
    zd = [xd[i] * exp(-xd[i]^2 - yd[i]^2) for i = 1:100]

    GR.setviewport(0.15, 0.95, 0.1, 0.9)
    GR.setwindow(-2, 2, -2, 2)
    GR.setspace(-0.5, 0.5, 0, 90)
    GR.setmarkersize(1)
    GR.setmarkertype(GR.MARKERTYPE_SOLID_CIRCLE)
    GR.setcharheight(0.024)
    GR.settextalign(2, 0)
    GR.settextfontprec(3, 0)

    x, y, z = GR.gridit(xd, yd, zd, 200, 200)
    levels = linspace(-0.5, 0.5, 20)
    GR.surface(x, y, z, 5)
    GR.contour(x, y, levels, z, 0)
    GR.polymarker(xd, yd)
    GR.axes(0.25, 0.25, -2, -2, 2, 2, 0.01)

    GR.updatews()
end

win = Gtk.GtkWindow("Gtk")

canvas = Gtk.GtkCanvas(500, 500)
Gtk.push!(win, canvas)

Gtk.draw(paint, canvas)
Gtk.showall(win)

signal_connect(win, :destroy) do widget
    exit(0)
end

Gtk.gtk_main()


text:
#=
  to_num in Julia + ConstraintSolver.jl

  to_num(List, Base, Num) converts a list of integers to a number for
  a base Base. It is bidirectional but it is really recommended that
  the length of List is fixed.

  See examples below.

  Model created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#
using ConstraintSolver, JuMP
using Cbc, GLPK
const CS = ConstraintSolver
include("constraints_utils.jl")

#
# to_num and increasing are defined in constraints_utils.jl
#
function to_num_test(n=5,base=10)

    cbc_optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)
    glpk_optimizer = optimizer_with_attributes(GLPK.Optimizer)

    model = Model(optimizer_with_attributes(CS.Optimizer,   "all_solutions"=> true,
                                                            # "all_optimal_solutions"=>true,
                                                            "logging"=>[],

                                                            # "traverse_strategy"=>:BFS,
                                                            "traverse_strategy"=>:DFS, # <-
                                                            # "traverse_strategy"=>:DBFS,

                                                            # "branch_split"=>:Smallest,
                                                            # "branch_split"=>:Biggest,
                                                            "branch_split"=>:InHalf, # <-

                                                            # "simplify"=>false,
                                                            # "simplify"=>true, # default

                                                            "time_limit"=>16,

                                                            # "lp_optimizer" => cbc_optimizer, # 45.8s
                                                            # "lp_optimizer" => glpk_optimizer, # 21.5s
                                        ))

    @variable(model, 0 <= x[1:n] <= base-1, Int)
    @variable(model, 0 <= v <= 10^n-1, Int)

    to_num(model, x, base, v)
    increasing(model, x)
    @constraint(model, x in CS.AllDifferent())

    # Solve the problem
    println("solve")
    optimize!(model)

    status = JuMP.termination_status(model)
    println("status:$status")
    if status == MOI.OPTIMAL
        xx = convert.(Integer,JuMP.value.(x))
        vv = convert.(Integer,JuMP.value.(v))

        println("x:$xx")
        println("v:$vv")

        num_sols = MOI.get(model, MOI.ResultCount())
        println("\nnum_sols:$num_sols\n")

        for sol in 1:num_sols
            println("solution #$sol")
            xx = convert.(Integer,JuMP.value.(x,result=sol))
            vv = convert.(Integer,JuMP.value.(v,result=sol))
            println("$xx  <=> $vv\n")

        end
    end
end

to_num_test()

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
@testset "test MILP relaxation API" begin
    PR.logger_config!("error")
    m = Model(cbc_optimizer)
    @variable(m, -1.0 <= x <= 1.0)
    @variable(m, y)
    formulation_info =
        construct_univariate_relaxation!(m, a -> a^3, x, y, collect(-1.0:1.0:1.0), true)

    @test size(formulation_info.variables[:delta_1])[1] == 2
    @test size(formulation_info.variables[:delta_2])[1] == 2
    @test size(formulation_info.variables[:z])[1] == 2
    @test haskey(formulation_info.constraints, :x)
    @test haskey(formulation_info.constraints, :y)
    @test haskey(formulation_info.constraints, :first_delta)
    @test haskey(formulation_info.constraints, :below_z)
    @test haskey(formulation_info.constraints, :above_z)
end

@testset "test LP relaxation API" begin
    PR.logger_config!("error")
    m = Model(cbc_optimizer)
    @variable(m, -1.0 <= x <= 1.0)
    @variable(m, y)
    formulation_info =
        construct_univariate_relaxation!(m, a -> a^3, x, y, collect(-1.0:1.0:1.0), false)

    @test size(formulation_info.variables[:lambda])[1] == 4
    @test haskey(formulation_info.constraints, :sum_lambda)
    @test haskey(formulation_info.constraints, :x)
    @test haskey(formulation_info.constraints, :y)
end

@testset "test bilinear relaxation API (McCormick)" begin
    replicates = 10
    tolerance = 1e-3
    for r = 1:replicates
        x_lb, x_ub = 10 * rand(2) .* [-1, 1]
        y_lb, y_ub = 10 * rand(2) .* [-1, 1]

        m = JuMP.Model(ipopt_optimizer)
        JuMP.@variable(m, x_lb <= x <= x_ub)
        JuMP.@variable(m, y_lb <= y <= y_ub)
        JuMP.@variable(m, z)
        JuMP.@objective(m, Min, z)
        JuMP.@constraint(m, x * y == z)
        status = JuMP.optimize!(m)

        rm = JuMP.Model(cbc_optimizer)
        JuMP.@variable(rm, x_lb <= x <= x_ub)
        JuMP.@variable(rm, y_lb <= y <= y_ub)
        JuMP.@variable(rm, z)
        JuMP.@objective(rm, Min, z)
        construct_bilinear_relaxation!(rm, x, y, z, [x_lb, x_ub], [y_lb, y_ub])
        rstatus = JuMP.optimize!(rm)

        @test(JuMP.objective_value(rm) <= JuMP.objective_value(m) + tolerance)
        @test(rstatus == status)

        JuMP.set_objective_sense(m, MOI.MAX_SENSE)
        JuMP.set_objective_sense(rm, MOI.MAX_SENSE)

        status = JuMP.optimize!(m)
        rstatus = JuMP.optimize!(rm)

        @test(JuMP.objective_value(rm) >= JuMP.objective_value(m) - tolerance)
        @test(rstatus == status)
    end
end

function refine!(partition, val)
    num_partitions = length(partition) - 1
    val_in_partition = 0
    for i = 1:num_partitions
        if val > partition[i] && val < partition[i+1]
            val_in_partition = i
        end
    end
    left = val - (val - partition[val_in_partition]) * 0.5
    right = val + (partition[val_in_partition+1] - val) * 0.5
    insert!(partition, val_in_partition + 1, left)
    insert!(partition, val_in_partition + 2, right)
    return (val - partition[val_in_partition]) * 0.5
end

@testset "test bilinear MIP relaxation API" begin
    gap_x = Float64[]
    gap_y = Float64[]
    point = 0.0

    x_p = [-10.0, 10.0]
    y_p = [-10.0, 10.0]
    for _ = 1:20
        x_lb, x_ub = x_p[1], x_p[end]
        y_lb, y_ub = y_p[1], y_p[end]
        epsilon_x = refine!(x_p, point)
        epsilon_y = 10.0
        m = JuMP.Model(cbc_optimizer)
        JuMP.@variable(m, x_lb <= x <= x_ub)
        JuMP.@variable(m, y_lb <= y <= y_ub)
        JuMP.@variable(m, z)
        construct_bilinear_relaxation!(m, x, y, z, x_p, [y_lb, y_ub])
        JuMP.@constraint(m, x == point)
        JuMP.@constraint(m, y == point)
        JuMP.@objective(m, Min, z)
        JuMP.optimize!(m)
        min_y = JuMP.objective_value(m)
        JuMP.set_objective_sense(m, MOI.MAX_SENSE)
        JuMP.optimize!(m)
        max_y = JuMP.objective_value(m)
        @test abs(max_y) ≈ epsilon_x * epsilon_y atol = 1e-5
        @test abs(min_y) ≈ epsilon_x * epsilon_y atol = 1e-5
    end
    for _ = 1:20
        x_lb, x_ub = x_p[1], x_p[end]
        y_lb, y_ub = y_p[1], y_p[end]
        epsilon_x = 10.0
        epsilon_y = refine!(y_p, point)
        m = JuMP.Model(cbc_optimizer)
        JuMP.@variable(m, x_lb <= x <= x_ub)
        JuMP.@variable(m, y_lb <= y <= y_ub)
        JuMP.@variable(m, z)
        construct_bilinear_relaxation!(m, x, y, z, [x_lb, x_ub], y_p)
        JuMP.@constraint(m, x == point)
        JuMP.@constraint(m, y == point)
        JuMP.@objective(m, Min, z)
        JuMP.optimize!(m)
        min_x = JuMP.objective_value(m)
        JuMP.set_objective_sense(m, MOI.MAX_SENSE)
        JuMP.optimize!(m)
        max_x = JuMP.objective_value(m)
        @test abs(max_x) ≈ epsilon_x * epsilon_y atol = 1e-5
        @test abs(min_x) ≈ epsilon_x * epsilon_y atol = 1e-5
    end
end


text:
function JuMP.set_objective(m::InnerBilevelModel, sense::MOI.OptimizationSense,
    f::JuMP.AbstractJuMPScalar)
    set_mylevel_obj_sense(m, sense)
    set_mylevel_obj_function(m, f)
    level_f = replace_variables(f, bilevel_model(m), mylevel_model(m), mylevel_var_list(m), level(m))
    JuMP.set_objective(mylevel_model(m), sense, level_f)
end
JuMP.objective_sense(m::InnerBilevelModel) = mylevel_obj_sense(m)
JuMP.objective_function_type(m::InnerBilevelModel) = typeof(mylevel_obj_function(m))
JuMP.objective_function(m::InnerBilevelModel) = mylevel_obj_function(m)
function JuMP.objective_function(m::InnerBilevelModel, FT::Type)
    mylevel_obj_function(m) isa FT || error("The objective function is not of type $FT")
    mylevel_obj_function(m)
end

function JuMP.relative_gap(bm::BilevelModel)::Float64
    _check_solver(bm)
    return MOI.get(bm.solver, MOI.RelativeGap())
end
function JuMP.relative_gap(bm::UpperModel)::Float64
    return JuMP.relative_gap(bm.m)
end
function JuMP.dual_objective_value(bm::BilevelModel; result::Int=1)::Float64
    _check_solver(bm)
    return MOI.get(bm.solver, MOI.DualObjectiveValue(result))
end
function JuMP.objective_bound(bm::BilevelModel)::Float64
    _check_solver(bm)
    return MOI.get(bm.solver, MOI.ObjectiveBound())
end
function JuMP.objective_bound(bm::UpperModel)::Float64
    return JuMP.objective_bound(bm.m)
end

function JuMP.set_objective(m::BilevelModel, sense::MOI.OptimizationSense,
    f::JuMP.AbstractJuMPScalar)
    bilevel_obj_error()
end
JuMP.objective_sense(m::BilevelModel) = JuMP.objective_sense(m.upper)# bilevel_obj_error()
JuMP.objective_function_type(model::BilevelModel) = bilevel_obj_error()
JuMP.objective_function(model::BilevelModel) = bilevel_obj_error()
function JuMP.objective_function(model::BilevelModel, FT::Type)
    bilevel_obj_error()
end

bilevel_obj_error() = error("There is no objective for BilevelModel use Upper(.) and Lower(.)")

function JuMP.objective_value(model::BilevelModel)
    _check_solver(model)
    return MOI.get(model.solver, MOI.ObjectiveValue())
end
function JuMP.objective_value(model::UpperModel)
    return JuMP.objective_value(model.m)
end
function JuMP.objective_value(model::LowerModel)
    return lower_objective_value(model.m)
end

function lower_objective_value(model::BilevelModel; result::Int = 1)
    f = JuMP.objective_function(Lower(model))
    # Evaluate the lower objective expression
    return JuMP.value(f, v -> JuMP.value(v, result = result))
    # return JuMP.value(v -> inner_ref_to_value(Lower(model), v), f)
end


text:
immutable AccelTurnrate
    a::Float64
    ω::Float64
end
Base.show(io::IO, a::AccelTurnrate) = @printf(io, "AccelTurnrate(%6.3f,%6.3f)", a.a, a.ω)
Base.length(::Type{AccelTurnrate}) = 2
Base.convert(::Type{AccelTurnrate}, v::Vector{Float64}) = AccelTurnrate(v[1], v[2])
function Base.copy!(v::Vector{Float64}, a::AccelTurnrate)
    v[1] = a.a
    v[2] = a.ω
    v
end
function propagate(veh::Vehicle, action::AccelTurnrate, roadway::Roadway, Δt::Float64; n_integration_steps::Int=4)

    a = action.a # accel
    ω = action.ω # turnrate

    x = veh.state.posG.x
    y = veh.state.posG.y
    θ = veh.state.posG.θ
    v = veh.state.v

    δt = Δt / n_integration_steps

    for i in 1 : n_integration_steps
        x += v*cos(θ)*δt
        y += v*sin(θ)*δt
        θ += ω*δt
        v += a*δt
    end

    posG = VecSE2(x, y, θ)
    VehicleState(posG, roadway, v)
end
function Base.get(::Type{AccelTurnrate}, rec::SceneRecord, roadway::Roadway, vehicle_index::Int, pastframe::Int=0)
    accel = get(ACC, rec, roadway, vehicle_index, pastframe)
    turnrate = get(TURNRATEG, rec, roadway, vehicle_index, pastframe)
    AccelTurnrate(accel, turnrate)
end


text:
export Harmonic

abstract type BdInter end
Base.broadcastable(x::BdInter) = Ref(x)
struct NoBond<:BdInter end

struct Harmonic<:BdInter
    k::Array{Float64,1}
    l0::Array{Float64,1}
    function Harmonic(;k::Array{Float64,1}=[0.0],
                        l0::Array{Float64,1}=[0.0])
        new(k::Array{Float64,1},
            l0::Array{Float64,1})
    end
end

function init(::Harmonic,input::Dict)
    k_matrix=zeros(length(keys(input["harmonic bond"]["k"])))
    l0_matrix=zeros(length(keys(input["harmonic bond"]["l0"])))

    for atomi in 1:length(keys(input["harmonic bond"]["k"]))
        k_matrix[atomi]=input["harmonic bond"]["k"][string(atomi)]
        l0_matrix[atomi]=input["harmonic bond"]["l0"][string(atomi)]
    end
    bdinter=Harmonic(k=k_matrix,l0=l0_matrix)
    return bdinter
end


text:
mutable struct PIMDP
    states
    actions
    Pmin
    Pmax
    imdp_label_dict
    accepting_labels
    sink_labels
    initial_state
    policy
    state_history
    trajectory
end

"
Construct a PIMDP from an MDP by directly specifying unsafe states.
"
# TODO: Need multiple ways to do this. This is a hackish way for now.
function construct_PIMDP_from_IMDP(imdp, unsafe_states)

    # Problem 1: What are the accepting states?
    # Answer: the accepting states are the unsafe states. Then, we perform verification over the unsafe states and take the complementary
    # TODO: This is very hackkkkky
    sink_labels = zeros(1, length(imdp.states))
    sink_labels[unsafe_states] .= 1
    pimdp = PIMDP(imdp.states, imdp.actions, imdp.Pmin, imdp.Pmax, nothing, sink_labels, [],
                  nothing, nothing, nothing, nothing)

    return pimdp
end

"""
Construct a product IMDP from a (multimode) IMDP and DFA.
"""
function construct_DFA_IMDP_product(dfa, imdp)
    # Intial state
    qinit = δ(dfa.transitions, dfa.initial_state, imdp.labels[imdp.initial_state])
    sizeQ = length(dfa.states)
    pinit = (imdp.initial_state, qinit)
    dfa_acc_state = dfa.accepting_state
    dfa_sink_state = dfa.sink_state

    # Transition matrix size will be Nx|Q| -- or the original transition matrix permutated with the number of states in DFA
    Pmin = imdp.Pmin
    Pmax = imdp.Pmax
    N, M = size(Pmin)
    Pmin_new = spzeros(N*sizeQ, M*sizeQ)
    Pmax_new = spzeros(N*sizeQ, M*sizeQ)

    pimdp_states = []
    pimdp_actions = imdp.actions
    sizeA = length(pimdp_actions)
    @debug "Actions: " pimdp_actionsA
    for s in imdp.states
        for q in dfa.states
            new_state = (s, q)
            push!(pimdp_states, new_state)
            # Check for transitions to other states
        end
    end

    for sq in pimdp_states
        for a in pimdp_actions
            for sqp in pimdp_states
                qp_test = δ(dfa.transitions, sq[2], imdp.labels[sq[1]])
                if qp_test == sqp[2]
                    # Get the corresponding entry of the transition interval matrices
                    row_idx = (sq[1]-1)*sizeQ*sizeA + (sq[2]-1)*sizeA + a
                    if (sq[2] == dfa_acc_state && sqp[2] == dfa_acc_state) || (sq[2] == dfa_sink_state && sqp[2] == dfa_sink_state)
                        # Flush out the old probabilities
                        col_idx = (sq[1]-1)*sizeQ + sq[2]
                        # Pmin_new[row_idx, :] .= 0.  # Redundant??
                        # Pmax_new[row_idx, :] .= 0.
                        Pmin_new[row_idx, col_idx] = 1.0
                        Pmax_new[row_idx, col_idx] = 1.0
                    else
                        col_idx = (sqp[1]-1)*sizeQ + sqp[2]
                        Pmin_new[row_idx, col_idx] = Pmin[(sq[1]-1)*sizeA + a, sqp[1]]
                        Pmax_new[row_idx, col_idx] = Pmax[(sq[1]-1)*sizeA + a, sqp[1]]
                    end
                end
            end
        end
    end

    labels = zeros(1, M*sizeQ)
    if !isnothing(dfa_acc_state)
        labels[dfa_acc_state:sizeQ:M*sizeQ] .= 1
    end

    sink_labels = zeros(1, M*sizeQ)
    if !isnothing(dfa_sink_state)
        sink_labels[dfa_sink_state:sizeQ:M*sizeQ] .= 1
    end

    pimdp = PIMDP(pimdp_states, imdp.actions, Pmin_new, Pmax_new, imdp.labels, labels, sink_labels, nothing, nothing, nothing, nothing)
    return pimdp
end

"""
Write a PIMDP object to file in a format consistent with Morteza's synthesis tool.
"""
function write_pimdp_to_file(pimdp, filename)
    open(filename, "w") do f
        state_num = length(pimdp.states)
        action_num =length(pimdp.actions)
        @printf(f, "%d \n", state_num)
        @debug "Length actions: " length(pimdp.actions)
        @printf(f, "%d \n", length(pimdp.actions))
        # Get number of accepting states from the labels vector
        @printf(f, "%d \n", sum(pimdp.accepting_labels))
        pimdp_acc_labels = pimdp.accepting_labels
        acc_states = findall(>(0.), pimdp_acc_labels[:])
        [@printf(f, "%d ", acc_state-1) for acc_state in acc_states]
        @printf(f, "\n")
        pimdp_sink_labels = pimdp.sink_labels
        sink_states = findall(>(0.), pimdp_sink_labels[:])

        for i=1:state_num
            if isnothing(sink_states) || !(i∈sink_states)
                for action in pimdp.actions
                    row_idx = (i-1)*action_num + action
                    ij = findall(>(0.), pimdp.Pmax[row_idx, :])
                    # Something about if the upper bound is less than one? Perhaps for numerical issues?
                    @debug action, i
                    psum = sum(pimdp.Pmax[row_idx, :])
                    psum >= 1 ? nothing : throw(AssertionError("Bad max sum: $psum"))
                    for j=ij
                        @printf(f, "%d %d %d %f %f", i-1, action-1, j-1, pimdp.Pmin[row_idx, j], pimdp.Pmax[row_idx, j])
                        if (i < state_num || j < ij[end] || action < action_num)
                            @printf(f, "\n")
                        end
                    end
                end
            else
                [@printf(f, "%d %d %d %f %f\n", i-1, j, i-1, 1.0, 1.0) for j=0:action_num-1]
            end
        end
    end
end

function create_dot_graph(pimdp::PIMDP, filename::String)
    open(filename, "w") do f
        println(f, "digraph G {")
        println(f, "  rankdir=LR\n  node [shape=\"circle\"]\n  fontname=\"Lato\"\n  node [fontname=\"Lato\"]\n  edge [fontname=\"Lato\"]")
        println(f, "  size=\"8.2,8.2\" node[style=filled,fillcolor=\"#FDEDD3\"] edge[arrowhead=vee, arrowsize=.7]")

        # Initial State
        # TODO: Don't hardcode intial state
        @printf(f, "  I [label=\"\", style=invis, width=0]\n  I -> %d\n", 1)

        i = 1
        for state in pimdp.states
            @printf(f, "  %d [label=<(s%d,q%d)>, xlabel=<%d>]\n", i, state[1], state[2], pimdp.accepting_labels[(state[1]-1)*3 + state[2]])

            for action in pimdp.actions
                q_size = 3
                row_idx = (state[1]-1)*length(pimdp.actions)*q_size + (state[2]-1)*length(pimdp.actions) + action

                for idx in findall(>(0.), pimdp.Pmax[row_idx, :])
                    state_p = pimdp.states[idx]
                    col_idx = (state_p[1]-1)*q_size + state_p[2]
                    # TODO: fix the second index to be correct
                    @printf(f, "  %d -> %d [label=<a%d: %.1f-%.1f >]\n", i, idx, action, pimdp.Pmin[row_idx,col_idx], pimdp.Pmax[row_idx,col_idx])
                end
            end
            i+=1
        end
        println(f, "}")
    end
end

"""
Read a PIMDP from file.
"""
function read_pimdp_from_file(filename)
    pimdp = nothing
    open(filename, "r") do f
        num_states = tryparse(Int, readline(f))
        num_actions = tryparse(Int, readline(f))
        num_sink_states = tryparse(Int, readline(f))
        accept_states = parse.(Int, split(readline(f))) .+ 1
        Pmin = spzeros(num_states*num_actions, num_states)
        Pmax = spzeros(num_states*num_actions, num_states)
        while !eof(f)
            row_split = split(readline(f))
            if length(row_split) > 0
                q0 = parse(Int, row_split[1])
                a = parse(Int, row_split[2]) + 1
                qt = parse(Int, row_split[3])
                plow = parse(Float64, row_split[4])
                phigh = parse(Float64, row_split[5])
                row_idx = q0*num_actions + a
                col_idx = qt + 1
                Pmin[row_idx, col_idx] = plow
                Pmax[row_idx, col_idx] = phigh
            end
        end

        pimdp = PIMDP(collect(1:num_states), collect(1:num_actions), Pmin, Pmax,
        nothing, accept_states, [], nothing, nothing, nothing, nothing)
    end

    return pimdp
end

"""
Quickly Validate a (P)IMDP object.
"""
function validate_pimdp(pimdp)
    for (minrow, maxrow) in zip(eachrow(pimdp.Pmin), eachrow(pimdp.Pmax))
        @assert sum(maxrow) >= 1
        @assert sum(minrow) <= 1
    end
end

"""
Determine the IMDP state from the current PIMDP index and # of DFA states
"""
function pimdp_col_to_imdp_state(col, num_dfa_states)
    return ceil(col/num_dfa_states)
end

"""
% computes the true distribution given a range of distributions

% input:    pmin: lower bound transition probabilities
%           pmax: upper bound transition probabilities
%           indSorted: index of states according to their sorted list
% output:
%           p: true transition probability
"""
function get_true_transition_probabilities(pmin, pmax, indSorted)
    p = zeros(1,length(pmax))
    used = sum(pmin)
    remain = 1 - used

    for i in indSorted
        if pmax[i] <= (remain + pmin[i])
            p[i] = pmax[i]
        else
            p[i] = pmin[i] + remain
        end
        remain = maximum([0., remain - (pmax[i] - pmin[i]) ])
    end

    return p
end

"""
Get the PIMDP rows with the same IMDP state (non-sink/accept)
"""
function get_pimdp_row_idx_same_imdp(pimdp_state, dfa, num_actions, mode)
    num_dfa_states = length(dfa.states)
    imdp_idx = pimdp_state[1]
    pimdp_row_idxs = []
    for i=1:num_dfa_states
        if i == dfa.sink_state || i == dfa.accepting_state
            continue
        end
        push!(pimdp_row_idxs, map_pimdp_state_to_Pmat_row((imdp_idx, i), num_dfa_states, num_actions, mode))
    end

    return pimdp_row_idxs
end

"""
Get the PIMDP states 1-ring around the input PIMDP state
"""
function get_pimdp_ring_states(pimdp_state)
    # Given the imdp, get the ring states of the PIMDP
    # For how ours is constucted, this is
    imdp_state_idx = pimdp_state[1]

    # TODO: pass in the # of imdp states per row
    ring_imdp_states = get_imdp_ring_states(imdp_state_idx, 32, 1024)

    pimdp_ring_states = [(state, pimdp_state[2]) for state in ring_imdp_states]

    # ring_imdp_states = [imdp_state_idx + 1, imdp_state_idx - 1]
    return pimdp_ring_states
end


text:
"""
    coupled_states(E[; i₀=1])

Find all states coupled by the energy expression `E`, starting from
the state with index `i₀`. This can be useful to reduce the necessary
basis or to generate invariant sets for split-operator propagation.
"""
function coupled_states(E::AbstractSparseMatrix; i₀=1)
    m = size(E,1)
    visited = falses(m)
    visited[i₀] = true

    rows = rowvals(E)

    istack = [i₀]
    while !isempty(istack)
        icur = pop!(istack)
        neighbours = rows[nzrange(E, icur)]
        for n in neighbours
            if !visited[n]
                push!(istack, n)
                visited[n] = true
            end
        end
    end

    visited
end

function invariant_sets(E::AbstractSparseMatrix)
    m = size(E,1)
    visited = falses(m)

    sets = Vector{Vector{Int}}()

    while !all(visited)
        icur = findfirst(.!visited)
        set = coupled_states(E; i₀=icur)
        visited[:] .|= set
        j = findall(set)
        # If the state icur is the only one in its set, it may be that
        # it is actually not coupled to anything.
        length(j) == 1 && iszero(E[icur,j]) && continue
        push!(sets, j)
    end
    sets
end

export coupled_states, invariant_sets


text:
#Copyright (c) 2015 Rex Computing and Isaac Yonemoto

#see LICENSE.txt

#this work was supported in part by DARPA Contract D15PC00135

#max-speed-2.jl

#testing the use of max versus ternary operator for uint16s for __minimum_data_width

const count = 1000000

function testy()

  test1 = [(rand() < 0.5) ? UInt16(0) : rand(UInt16) for i=1:count]

  res::UInt16 = 0
  println("using max()")
  tic()
  for (i=1:count)
    res += max(0, test1[i] - 1)
  end
  toc()
  println(res)

  println("using max() with potential conversion")
  tic()
  for (i=1:count)
    res += UInt16(max(0, test1[i] - 1))
  end
  toc()
  println(res)

  res = 0
  println("using ternary op")
  tic()
  for (i=1:count)
    res += (test1[i] == 0) ? 0 : test1[i] - 1
  end
  toc()
  println(res)

  res = 0
  println("using ternary op with assignment")
  tic()
  for (i=1:count)
    temp = test1[i]
    res += (temp == 0) ? 0 : temp - 1
  end
  toc()
  println(res)
end

testy()
testy()
testy()

#and let's double check that the max result is the same.
function testy2()  #variable type annotations only valid within function blocks
  test2 = [0, 1, 2] #result should be [0, 0, 1]
  r::UInt16 = 0
  for (i=1:3)
    r = max(0, test2[i] - 1)
    println(r)
    println(typeof(r))
  end
end
testy2()

#30 Aug 2015
#results, third run only, 'res printlns' edited out.
#using max()
#elapsed time: 0.001808343 seconds
#using max() with potential conversion
#elapsed time: 0.001754662 seconds
#using ternary op
#elapsed time: 0.007478182 seconds
#using ternary op with assignment
#elapsed time: 0.007068262 seconds
#0
#0
#1
#result:  Using max() is still faster, and produces the desired results.

max_stars_repo_name: REX-Computing/unumjl

max_issues_repo_name: REX-Computing/unumjl

max_forks_repo_name: REX-Computing/unumjl


text:
using Pkg
Pkg.activate(".")

using tdpe

tdpe.print_solutions()

max_stars_repo_name: mithodin/ttpe

max_issues_repo_name: mithodin/ttpe

max_forks_repo_name: mithodin/ttpe


text:
"""
    AbstractVectorFold{T, P}
An abstract type used as an interface for folded vectors such as `VectorFold`.
To implement the interface and inherit from it, a new structure must define three fields:
- `pattern::P`. Note that both `length(::P)` and `rand(::P)` methods must be available
- `gap::T`
- `folds::int`
"""
abstract type AbstractVectorFold{T} <: AbstractVector{T} end

"""
    PatternFold{T, P}
A `Union` type used as an interface for folded patterns such as `VectorFold`.
To implement the interface and inherit from it, a new structure `MyFold{T[,P]}` must define three fields:
- `pattern::P`. Note that both `length(::P)` and `rand(::P)` methods must be available
- `gap::T`S
- `folds::int`
Finally one can redefine PatternFold{T}
```julia
PatternFold{T} = Union{AbstractVectorFold{T}, IntervalsFold{T}, MyFold{T[,P]}}
```
"""
PatternFold = Union{AbstractVectorFold, IntervalsFold}

"""
    pattern(<:PatternFold)
Return the *pattern* of any `PatternFold`. The pattern defines the occurences of the first fold.
"""
pattern(pf) = pf.pattern

"""
    gap(<:PatternFold)
Return the *gap* between the starts of consecutive folds.
"""
gap(pf) = pf.gap

"""
    folds(<:PatternFold)
Return the number of *folds*. An infinite folded pattern returns `0`.
"""
folds(pf) = pf.folds

# Forwards isempty, ndims
@forward PatternFold.pattern Base.isempty, Base.ndims

"""
    pattern_length(pf<:PatternFold)
Return the length of the basic pattern of `pf`.
"""
pattern_length(pf) = length(pattern(pf))

"""
    length(pf<:PatternFold)
Return the length of `pf` if unfolded.
"""
Base.length(pf::PatternFold) = pattern_length(pf) * folds(pf)
Base.size(pf::PatternFold) = (length(pf),)

"""
    rand(pf<:PatternFold)
Returns a random value of `pf` as if it was unfolded.
"""
function Base.rand(pf::PF) where {PF <: PatternFold}
    return Base.rand(pattern(pf)) + Base.rand(0:(folds(pf) - 1)) * gap(pf)
end

"""
    reset_pattern!(<:PatternFold)
Reset the *unfolded* pattern to the first fold.
"""
reset_pattern!(mvf) = set_fold!(mvf, 1)

"""
    fold(v::V, depth = 0)
returns a suitable `VectorFold`, which when unfolded gives the Vector V.
"""
function fold(v::V, depth = 0; kind = :mutable) where {T <: Real, V <: AbstractVector{T}}
    l = length(v)
    for i in 1:(l ÷ 2)
        gap = v[i + 1] - v[1]
        fold, r = divrem(l, i)
        if  r == 0 && check_pattern(v, i, gap, fold)
            # return VectorFold(fold(v[1:i], depth + 1), gap, fold)
            return make_vector_fold(v[1:i], gap, fold, kind)
        end
    end
    if depth == 0
        @warn "No non-degenerate patterns have been found" v
        return make_vector_fold(v, zero(T), 1, kind)
    else
        return v
    end
end

"""
    make_vector_fold(pattern, gap, fold, kind = :mutable)
A dispatcher to construct a folded vector. The `kind` of vector can be set to either `:mutable` (default) or `:immutable`. The default is faster in most cases but it depends on the `pattern`, `gap`, and `fold` parameters. For critical code, it is recommended to benchmark both options.
"""
function make_vector_fold(pattern, gap, fold, kind = :mutable)
    return make_vector_fold(pattern, gap, fold, Val(kind))
end
function make_vector_fold(isf; kind = :mutable)
    vf = kind == :mutable ? "VectorFold" : "IVectorFold"
    str = "The pattern of the interval is not a point." *
        " The IntervalsFold cannot be converted to a $vf"
    @assert is_points(isf) str

    return make_vector_fold([value(pattern(isf), :a)], gap(isf), folds(isf), kind)
end

function Base.print_array(io::IO, X::AbstractVectorFold)
    print(io, "\tPattern: $(pattern(X))\n\tGap: $(gap(X))\n\tFolds: $(folds(X))")
end

function Base.getindex(vf::AbstractVectorFold, key)
    d, r = divrem(key + 1, pattern_length(vf))
    return (d - 1) * gap(vf) + pattern(vf, r + 1)
end
Base.getindex(vf::AbstractVectorFold, key...) = map(k -> getindex(vf, k), key)

"""
    Base.rand(::Vector{AbstractVectorFold})
    Extend the `Base.rand` function to `Vector{AbstractVectorFold}`.
"""
Base.rand(v::Set{<:AbstractVectorFold}) = rand(collect(Iterators.flatten(v)))


text:
# """ Clear the screen, print an error message, and kill the program """
# function die(msg)
#     write(stdout, "\x1b[2J")
#     write(stdout, "\x1b[H")
#     error(msg)
# end

##########
# CURSOR #
##########

mutable struct Cursor
    x::Int
    y::Int
    rx::Int
end

##########
# EDITOR #
##########
"Editor instance"
mutable struct Editor
    "row view offset"
    rowoff::Int

    "column view offset"
    coloff::Int

    "terminal width"
    width::Int

    "terminal height"
    height::Int

    "currently edited file"
    filename::String

    "current status message"
    statusmsg::String

    "time reaining to display status"
    status_time::Float64

    "true if the buffer has been edited"
    dirty::Bool

    "if true, the key input loop will exit"
    quit::Bool

    "the cursor position"
    csr::Cursor

    "the test buffer"
    rows::Rows

    "terminal hosting this editor"
    term::REPL.Terminals.TTYTerminal

    "used by commands to store variables"
    params::Dict{Symbol, Dict{Symbol, Any}}
end

"""create a new editor with default params"""
function Editor()
    rowoff = 0
    coloff = 0
    width = 0
    height = 0
    filename = ""
    statusmsg = ""
    status_time = 0
    dirty = false
    quit = false

    csr = Cursor(1,1,1)
    rows = Rows()
    term = REPL.Terminals.TTYTerminal(get(ENV, "TERM", @static Sys.iswindows() ? "" : "dumb"), stdin, stdout, stderr)

    params = Dict{Symbol, Dict{Symbol, Any}}()

    Editor(rowoff, coloff, width, height, filename,
        statusmsg, status_time, dirty, quit,
        csr, rows, term, params)
end

###################
# FILE OPERATIONS #
###################

"""custom expand path function to ensure consistency"""
function editorExpandPath(filename::String)
    if filename == ""
        return ""
    else
        return abspath(expanduser(filename))
    end
end

""" Open a file and read the lines into the ed.rows array """
function editorOpen(ed::Editor, filename::String)
    try
        filename = editorExpandPath(filename)

        # If no file exists, create it
        !isfile(filename) && open(filename, "w") do f end

        open(filename, "r") do file
            clearRows(ed.rows)
            for line in eachline(file)
                appendRow(ed.rows, line)
            end
            ed.filename = filename
            ed.csr.x = 1
            ed.csr.y = 1
            ed.dirty = false
        end
    catch Exception
        setStatusMessage(ed, "Cannot open file $filename")
    end
end

function editorOpen(ed::Editor)
    if ed.dirty
        confirm = editorPrompt(ed, "There are unsaved changes. Open another file? [y/n]: ")
        if confirm != "y"
            setStatusMessage(ed, "Open aborted")
            return
        end
    end

    filename = editorPrompt(ed, "Open file: ")
    filename = editorExpandPath(filename)

    if filename != ""
        editorOpen(ed, filename)
    else
        setStatusMessage(ed, "Open aborted")
    end
end

""" Save the buffer to the file using the current file name """
function editorSave(ed::Editor)
    editorSave(ed, "")
end

""" Save the contents of the file buffer to disc """
function editorSave(ed::Editor, path::String)
    prev_filename = ed.filename

    try
        if path == ""
            if ed.filename == ""
                ed.filename = editorPrompt(ed, "Save as: ")
                if ed.filename == ""
                    setStatusMessage(ed, "Save aborted")
                    return
                end
            end
        else
            ed.filename = editorExpandPath(path)
        end

        open(ed.filename, "w") do f
            write(f, rowsToString(ed.rows))
        end
        setStatusMessage(ed, "File saved: $(ed.filename)")
        ed.dirty = false
    catch Exception
        # There was an error saving, restore original filename
        setStatusMessage(ed, "Unable to save: $(ed.filename)")
        ed.filename = prev_filename
    end

end

##########
# VISUAL #
##########

""" Given an arrow input, move the cursor """
function moveCursor(ed::Editor, key::UInt32)
    if key == ARROW_LEFT
        if ed.csr.x > 1
            ed.csr.x -= 1
        elseif ed.csr.y > 1
            # At start of line, move to end of prev line
            ed.csr.y -= 1
            ed.csr.x = 1+length(ed.rows[ed.csr.y].chars)
        end

    elseif key == ARROW_RIGHT
        onrow = ed.csr.y <= length(ed.rows)
        if onrow && ed.csr.x <= length(ed.rows[ed.csr.y].chars)
            ed.csr.x += 1
        elseif ed.csr.y < length(ed.rows) && ed.csr.x == 1 + length(ed.rows[ed.csr.y].chars)
            # At end of line, move to next line
            ed.csr.y += 1
            ed.csr.x = 1
        end

    elseif key == ARROW_UP
        ed.csr.y > 1 && (ed.csr.y -= 1)

    elseif key == ARROW_DOWN
        ed.csr.y < length(ed.rows) && (ed.csr.y += 1)
    end

    # Snap to end of line if we are further out
    rowlen = ed.csr.y < length(ed.rows)+1 ? length(ed.rows[ed.csr.y].chars)+1 : 1
    ed.csr.x > rowlen && (ed.csr.x = rowlen)
end

moveCursor(ed, key::Key) = moveCursor(ed, UInt32(key))

""" Scroll the screen based on the cursor position """
function scroll(ed::Editor)
    ed.csr.rx = 1
    if ed.csr.y <= length(ed.rows)
        ed.csr.rx = renderX(ed.rows[ed.csr.y], ed.csr.x)
    end

    # Vertical scrolling
    if ed.csr.y < ed.rowoff+1
        ed.rowoff = ed.csr.y-1
    end
    if ed.csr.y >= ed.rowoff+1 + ed.height
        ed.rowoff = ed.csr.y - ed.height
    end

    # Horizontal scrolling
    if ed.csr.rx < ed.coloff+1
        ed.coloff = ed.csr.rx-1
    end
    if ed.csr.rx >= ed.coloff+1 + ed.width
        ed.coloff = ed.csr.rx - ed.width
    end
end

function drawRows(ed::Editor, buf::IOBuffer)
    for y = 1:ed.height
        filerow = y + ed.rowoff
        y != 1 && write(buf, "\r\n")

        write(buf, "\x1b[K") # Clear line
        if filerow > length(ed.rows)
            if y == div(ed.height, 3) && ed.width > 40 && length(ed.rows) == 0
                msg = "Acorn Editor"
                padding = div(ed.width - length(msg), 2)
                if padding > 0
                    write(buf, "~")
                    padding -= 1
                end
                while (padding -= 1) > 0
                    write(buf, " ")
                end
                write(buf, msg)
            else
                write(buf, "~");
            end
        else
            len = length(ed.rows[filerow].render) - ed.coloff
            len = clamp(len, 0, ed.width)
            #write(buf, ed.rows[filerow].render[1+ed.coloff : ed.coloff + len])
            get_color = get_color_function(ed.filename)
            printrow(buf,
                ed.rows[filerow].render[1+ed.coloff : ed.coloff + len],
                get_color
            )
        end
    end
    # Write a newline to prepare for status bar
    write(buf, "\r\n");
end

function drawStatusBar(ed::Editor, buf::IOBuffer)
    write(buf, "\x1b[7m") # invert colors
    col = 1

    # left padding
    write(buf, ' ')
    col += 1

    # filename
    filename = configGet(:status_fullpath) ? editorExpandPath(ed.filename) : splitdir(ed.filename)[2]
    filestatus = string(filename, ed.dirty ? " *" : "")

    for i = 1:min(div(ed.width,2), length(filestatus))
        write(buf, filestatus[i])
        col += 1
    end

    linenum = string(ed.csr.y)

    while col < ed.width - length(linenum)
        write(buf, ' ')
        col += 1
    end

    write(buf, linenum, ' ')

    write(buf, "\x1b[m") # uninvert colors

    # make line for message bar
    write(buf, "\r\n")
end

function drawStatusMessage(ed::Editor, buf::IOBuffer)
    write(buf, "\x1b[K")
    if time() - ed.status_time < 5.0
        write(buf, ed.statusmsg[1:min(ed.width, length(ed.statusmsg))])
    end
end

function setStatusMessage(ed::Editor, msg::String)
    ed.statusmsg = msg
    ed.status_time = time()
end

function refreshScreen(ed::Editor)

    # Update terminal size
    ed.height = REPL.Terminals.height(ed.term) - 2 # status + msg bar = 2
	@static Sys.iswindows() ? (ed.height -= 1) : ed.height
    ed.width = REPL.Terminals.width(ed.term)

    scroll(ed)

    buf = IOBuffer()

    write(buf, "\x1b[?25l") # ?25l: Hide cursor
    write(buf, "\x1b[H")    # H: Move cursor to top left

    drawRows(ed, buf)
    drawStatusBar(ed, buf)
    drawStatusMessage(ed, buf)

    @printf(buf, "\x1b[%d;%dH", ed.csr.y-ed.rowoff,
                                ed.csr.rx-ed.coloff)

    write(buf, "\x1b[?25h") # ?25h: Show cursor

    write(stdout, String(take!(buf)))
end

function editorPrompt(ed::Editor, prompt::String;
                      callback=nothing,
                      buf::String="",
                      showcursor::Bool=true) ::String
    while true
        statusmsg = string(prompt, buf)
        setStatusMessage(ed, string(prompt, buf))
        refreshScreen(ed)

        if showcursor
            # Position the cursor at the end of the line
            @printf(stdout, "\x1b[%d;%dH", 999, length(statusmsg)+1)
        end

        c = Char(readKey())

        if c == '\x1b'
            setStatusMessage(ed, "")
            callback != nothing && callback(ed, buf, c)
            return ""
        elseif c == '\r'
            if length(buf) != 0
                setStatusMessage(ed, "")
                callback != nothing && callback(ed, buf, c)
                return buf
            end
        elseif UInt32(c) == BACKSPACE && length(buf) > 0
            buf = buf[1:end-1]
        elseif !iscntrl(c) && UInt32(c) < 128
            buf = string(buf, c)
        end

        callback != nothing && callback(ed, buf, c)
    end
end

##############
## Keyboard ##
##############

function processKeypress(ed::Editor)
    c = readKey();

    if c == UInt32('\x1b') # Esc
        setStatusMessage(ed, "Press ctrl-q to quit")
    elseif c == ctrl_key('p')
        runCommand(ed)
    elseif (c == ARROW_LEFT
            || c == ARROW_UP
            || c == ARROW_RIGHT
            || c == ARROW_DOWN)
        moveCursor(ed, c)
    elseif c == PAGE_UP || c == PAGE_DOWN
        lines = ed.height
        while (lines-=1) > 0
            moveCursor(ed, c == PAGE_UP ? ARROW_UP : ARROW_DOWN)
        end
    elseif c == HOME_KEY
        ed.csr.x = 0
    elseif c == END_KEY
        ed.csr.y <= length(ed.rows) && (ed.csr.x = length(ed.rows[ed.csr.y].chars) + 1)
    elseif c == UInt32('\r')
        editorInsertNewline(ed)
    elseif c == BACKSPACE
        editorDelChar(ed)
    elseif c == DEL_KEY
        moveCursor(ed, ARROW_RIGHT)
        editorDelChar(ed)
    # elseif c == ctrl_key('l')
    #     # Refresh screen
    #     return
    elseif c == UInt32('\t')
        editorInsertTab(ed)
    elseif iscntrl(Char(c)) #&& isKeyBound(Char(c))
        if isKeyBound(Char(c))
            runCommand(ed, getKeyBinding(Char(c)))
        else
            setStatusMessage(ed, "Key is not bound...")
        end
    elseif !iscntrl(Char(c)) && c < 1000
        # Chars above 1000 are a ::Key, see src/terminal.jl
        editorInsertChar(ed, c)
    end

end

#####################
# Editor Operations #
#####################

function editorInsertChar(ed::Editor, c::UInt32)
    # The cursor is able to move beyond the last row
    ed.csr.y == length(ed.rows)+1 && appendRow(ed.rows, "")
    insertChar!(ed.rows[ed.csr.y], ed.csr.x, Char(c))
    ed.csr.x += 1
    ed.dirty = true
end

function editorInsertTab(ed::Editor)
    # The cursor is able to move beyond the last row
    ed.csr.y == length(ed.rows)+1 && appendRow(ed.rows, "")

    # Insert character(s) into the row data
    mv_fwd = insertTab!(ed.rows[ed.csr.y], ed.csr.x)

    ed.csr.x += mv_fwd

    ed.dirty = true
end

function editorDelChar(ed::Editor)
    ed.csr.y == length(ed.rows)+1 && return
    ed.csr.x == 1 && ed.csr.y == 1 && return

    if ed.csr.x > 1
        deleteChar!(ed.rows[ed.csr.y], ed.csr.x -1)
        ed.csr.x -= 1
        ed.dirty = true
    else
        # Move cursor to end of prev line
        ed.csr.x = 1+length(ed.rows[ed.csr.y-1].chars)
        appendRowString(ed.rows[ed.csr.y-1], ed.rows[ed.csr.y].chars)
        editorDelRow(ed, ed.csr.y)
        ed.csr.y -= 1
    end
end

function editorInsertRow(ed::Editor, i::Int, str::String)
    row = Row(str)
    update!(row)
    insert!(ed.rows, i, row)
end

function editorInsertNewline(ed::Editor)
    if ed.csr.x == 1
        editorInsertRow(ed, ed.csr.y, "")
    else
        row = ed.rows[ed.csr.y]
        before = row.chars[1:ed.csr.x-1]
        after = row.chars[ed.csr.x:end]
        editorInsertRow(ed, ed.csr.y + 1, after)
        row.chars = before
        update!(row)
    end
    ed.csr.y += 1
    ed.csr.x = 1
end

function editorDelRow(ed::Editor, i::Int)
    i < 1 || i > length(ed.rows) && return
    deleteat!(ed.rows, i)
    ed.dirty = true
end

function editorQuit(ed::Editor; force::Bool=false)
    if !force && ed.dirty
        setStatusMessage(ed,
            "File has unsaved changes. Save changes or use <ctrl-p>'quit !' to quit anyway.")
    else
        write(stdout, "\x1b[2J")
        write(stdout, "\x1b[H")
        ed.quit = true
        !isinteractive() && exit(0)
    end
end

############
# COMMANDS #
############

function runCommand(ed::Editor)
    cmd = editorPrompt(ed, "> ")
    runCommand(ed, strip(cmd))
end

function runCommand(ed::Editor, command_str::AbstractString)
    cmd_arr = split(command_str, ' ', limit=2)

    # Get the command
    cmd = strip(cmd_arr[1])

    # Blank, do nothing
    cmd == "" && return

    # Command must be a valid identifier
    if !Base.isidentifier(cmd)
        setStatusMessage(ed, "'$cmd' is not a valid command name")
        return
    end

    cmd_sym = Symbol(cmd)

    # Get arguments if there are any
    args = ""
    if length(cmd_arr) > 1
        args = cmd_arr[2]
    end

    # If the command exists, run it
    if cmd_sym in keys(COMMANDS)
        # join(args): convert Substring to String
        runCommand(COMMANDS[cmd_sym], ed, join(args))
    else
        setStatusMessage(ed, "'$cmd_sym' is not a valid command")
    end
end

function runCommand(c::Command, ed::Editor, args::AbstractString)
    c.cmd(ed, args)
end


text:
# FIXME: Remove once fma is implemented in julia
if !isdefined(:fma)
  for (T, t, i) in ((Float32, "float", "f32"), (Float64, "double", "f64"))
    @eval begin
      function fma(a::$T, b::$T, c::$T)
        Base.llvmcall($("""%4 = call $t @llvm.fma.$i($t %0, $t %1, $t %2)
                           ret $t %4"""), $T, ($T, $T, $T), a, b, c)
      end

      # HACK:
      # The next line is required due to the interplay of the current
      # implementation of llvmcall and LLVM intrinsics.
      try fma($T(1.0), $T(1.0), $T(1.0)); end
    end
  end
  @doc """
    Computes `a * b + c` using the fused multiply-add (FMA) operation. IEEE
    754-2008 states:

    > The operation `fusedMultiplyAdd(x, y, z)` computes `(x * y) + z` as if with
    > unbounded range and precision, rounding only once to the destination format.
    > No underflow, overflow, or inexact exception (see 7) can arise due to the
    > multiplication, but only due to the addition; and so fusedMultiplyAdd
    > differs from a multiplication operation followed by an addition operation.

    References:

      * IEEE 754-2008, IEEE Standard for Floating-Point Arithmetic, August 2008.
    """ -> fma
end

@doc """
  Computes the error-free fused multiply-add `a * b + c = x + y + z`, where
  `x = fl(fma(a, b, c))`. This algorithm requires 17 flops. Different from
  `err_fma_compr`, it is possible that `y` and `z` are overlapping.

  References:

    * S. Boldo, and J.M. Muller, Some functions computable with a fused-mac,
      Proceedings of the 17th Symposium on Computer Arithmetic, IEEE Computer
      Society Press, 2005.
    * S. Boldo, and J.M. Muller, Exact and approximated error of the FMA, IEEE
      Transactions on Computers 60, pp.157-164, 2011.
  """ ->
function err_fma{T<:IEEE754}(a::T, b::T, c::T)
  x = fma(a, b, c)
  p, q = err_fast_mul(a, b)
  s, z = err_add(c, q) # not necessary if |c| ≥ |a b| or |c| ≥ |p| (requires 6 flops)
  u, v = err_add(p, s)
  y = (u-x) + v
  x, y, z
end

@doc """
  Computes the error-free fused multiply-add `a * b + c = x + y + z`, where
  `x = fl(fma(a, b, c))` and `x`, `y`, `z` are non-overlapping. This algorithm
  requires 20 flops.

  References:

  * S. Boldo, and J.M. Muller, Some functions computable with a fused-mac,
  Proceedings of the 17th Symposium on Computer Arithmetic, IEEE Computer
  Society Press, 2005.
  * S. Boldo, and J.M. Muller, Exact and approximated error of the FMA, IEEE
  Transactions on Computers 60, pp.157-164, 2011.
  """ ->
function err_fma_compr{T<:IEEE754}(a::T, b::T, c::T)
  x = fma(a, b, c)
  p, q = err_fast_mul(a, b)
  s, t = err_add(c, q) # not necessary if |c| ≥ |a b| or |c| ≥ |p| (requires 6 flops)
  u, v = err_add(p, s)
  w = (u-x) + v
  y, z = err_fast_add(w, t) # compress floating-point errors (not necessary, requires 3 flops)
  x, y, z
end

@doc """
  Instead of accurately computing the error-free fused multiply-add as in
  `err_fma_compr`, return only one error term `y`: `|(a * b + c) - x - y| ≤
  ulp(y)/2`.

  References:

    * S. Boldo, and J.M. Muller, Exact and approximated error of the FMA, IEEE
      Transactions on Computers 60, pp.157-164, 2011.
  """ ->
function err_fma_nearest{T<:IEEE754}(a::T, b::T, c::T)
  x = fma(a, b, c)
  p, q = err_fast_mul(a, b)
  s, t = err_add(c, q) # not necessary if |c| ≥ |a b| or |c| ≥ |p| (requires 6 flops)
  u, v = err_add(p, s)
  y = ((u-x) + v) + t
  x, y
end

@doc """
  Instead of accurately computing the error-free fused multiply-add as in
  `err_fma`, return an approximated error term `y`: `|(a * b + c) - x - y| ≤
  14 * 2^-p * |x|`, where `p` is the precision of the floating-point type used.
  That is, following `x` there are at least `p - 2` correct bits in `x + y`.

  References:

    * S. Boldo, and J.M. Muller, Exact and approximated error of the FMA, IEEE
      Transactions on Computers 60, pp.157-164, 2011.
  """ ->
function err_fma_approx{T<:IEEE754}(a::T, b::T, c::T)
  x = fma(a, b, c)
  p, q = err_fast_mul(a, b)
  u, v = err_add(p, c)
  y = (u-x) + (q+v)
  x, y
end


text:
using Test, YaoSym, YaoBase, YaoArrayRegister

@testset "constructors" begin
    @test ket"011" isa SymReg
    @test (ket"011")' isa AdjointSymReg
end

@testset "partial_tr" begin
    @vars α β
    r = α * ket"101" + β * ket"111"
    @test partial_tr(r, 2:3) == (α + β) * ket"1"
end

@testset "conversion" begin
    r = zero_state(3)
    sr = SymReg(r)

    for i in eachindex(r.state)
        @test r.state[i] ≈ N(sr.state[i])
    end
end


text:
using Gillespie
import Random: seed!

function F_l(x,parms)
    (Y1,Y2) = x
    (c1,c2,c3) = parms
    [c1*Y1,c2*Y1*Y2,c3*Y2]
end

x0 = [1000,1000]
nu = [[1 0];[-1 1];[0 -1]]
parms = [10.0,0.01,10.0]
tf = 2.0
seed!(1234)

gillespie_result = ssa(x0,F_l,nu,parms,tf)
gillespie_data = ssa_data(gillespie_result)

jensen_result = ssa(x0,F_l,nu,parms,tf,algo=:jensen,max_rate=100000.0,thin=true)
jensen_data = ssa_data(jensen_result)

jensen_alljumps_result = ssa(x0,F_l,nu,parms,tf,algo=:jensen,max_rate=100000.0,thin=false)
jensen_alljumps_data = ssa_data(jensen_alljumps_result)


text:
module C1_2
function f(q::Vector{T}, q_dot::Vector{T}, xi::T) where T
    l1 = q[1]
    l2 = q[2]
    l3 = q[3]
    l1_dot = q_dot[1]
    l2_dot = q_dot[2]
    l3_dot = q_dot[3]

    l1_dot.*(1.38654137505663e-11*l1 + 5.72153164875459e-12*l2 - 1.96418867892903e-11*l3 + 8.32613864009457e-15) - l2_dot.*(1.2466152575692e-11*l1 - 2.06113736247302e-11*l2 + 8.05356093237194e-12*l3 + 1.40368182567312e-14) - l3_dot.*(1.48441253668614e-12*l1 + 2.6325066615507e-11*l2 - 2.77725822370262e-11*l3 - 5.56082225561306e-15)
end
end


text:
using BasisFunctions: FFreq, DiscreteWeight, DiscreteProductWeight,
    innerproduct_fourier_part, default_dict_innerproduct, Measure, Weight,
	OuterProductArray
import BasisFunctions: quadweights, stencilarray, name, weights,
	grid, subindices, isnormalized, supermeasure,
    support, unsafe_weightfun, strings, innerproduct_native, restrict, discretemeasure
export submeasure, DiscreteTensorSubWeight, SubWeight, DiscreteSubWeight

"""
    struct SubWeight{M,D,T} <: Weight{T}

Submeasure of a Weight.
"""
struct SubWeight{M<:Weight,D<:Domain,T} <: Weight{T}
    measure     ::  M
    domain      ::  D

    function SubWeight(measure::Weight, domain::Domain)
        @assert eltype(domain) == eltype(support(measure))
        new{typeof(measure), typeof(domain), eltype(domain)}(measure, domain)
    end
end

"""
    submeasure(measure::Weight, domain::Domain)

Restrict a measure to a subset of its support
"""
submeasure(measure::Weight, domain::Domain) = SubWeight(measure, domain)
name(m::SubWeight) = "Restriction of a measure"
supermeasure(measure::SubWeight) = measure.measure
support(measure::SubWeight) = measure.domain
unsafe_weightfun(m::SubWeight, x) = unsafe_weightfun(supermeasure(m), x)

"""
    restrict(measure::Weight, domain::Domain)

Restrict a measure to a subset of its support
"""
restrict(measure::Weight, domain::Domain) = submeasure(measure, domain)
strings(m::SubWeight) = (name(m), (string(support(m)),), strings(supermeasure(m)))
submeasure(measure::ProductWeight, domain::ProductDomain) = ProductWeight(map(submeasure, components(measure), components(domain))...)

quadweights(grid::SubGrid, measure::SubWeight) =
    quadweights(supergrid(grid), supermeasure(measure))[subindices(grid)]

function innerproduct_native(b1::Fourier, i::FFreq, b2::Fourier, j::FFreq, m::SubWeight{<:FourierWeight}; options...)
	d = support(m)
	if d isa AbstractInterval
		innerproduct_fourier_part(b1, i, b2, j, infimum(d), supremum(d))
	else
		default_dict_innerproduct(b1, i, b2, j, m; options...)
	end
end

"""
    struct DiscreteSubWeight{T,M<:DiscreteWeight{T},G<:AbstractGrid} <: DiscreteWeight{T}

Submeasure of a discrete measure.
"""
struct DiscreteSubWeight{T,M<:DiscreteWeight{T},G<:AbstractGrid} <: DiscreteWeight{T}
    supermeasure   :: M
    subgrid        :: G
    DiscreteSubWeight(measure::DiscreteWeight{T}, grid::AbstractGrid) where T = new{T,typeof(measure),typeof(grid)}(measure, grid)
    DiscreteSubWeight(measure::DiscreteSubWeight{T}, grid::AbstractGrid) where T = error()
end

submeasure(measure::DiscreteWeight, domain::Domain) = _discretesubmeasure(subgrid(points(measure), domain), weights(measure))
subindices(measure::DiscreteSubWeight) = subindices(measure.subgrid)
supermeasure(measure::DiscreteSubWeight) = measure.supermeasure
discretemeasure(grid::Union{SubGrid,ProductSubGrid}) = DiscreteSubWeight(discretemeasure(supergrid(grid)), grid)
_discretesubmeasure(grid::Union{SubGrid,ProductSubGrid},weights) = DiscreteSubWeight(discretemeasure(supergrid(grid),weights), grid)
restrict(measure::DiscreteWeight, domain::Domain) = DiscreteSubWeight(measure, subgrid(points(measure), domain))
weights(measure::DiscreteSubWeight) = subweights(measure, subindices(measure), weights(supermeasure(measure)))
subweights(_, subindices, w) = w[subindices]
points(measure::DiscreteSubWeight) = measure.subgrid
isnormalized(::DiscreteSubWeight) = false
name(m::DiscreteSubWeight) = "Restriction of a "*name(supermeasure(m))

function discretemeasure(grid::Union{SubGrid,ProductSubGrid}, weights::Ones)
    @assert size(grid) == size(weights)
    DiscreteSubWeight(discretemeasure(supergrid(grid)), grid)
end

"""
    const DiscreteTensorSubWeight{T,G,W} = DiscreteSubWeight{T,M,G} where {T,M<:DiscreteProductWeight,G<:ProductSubGrid}

A tensor product of discrete submeasures.
"""
const DiscreteTensorSubWeight{T,G,W} = DiscreteSubWeight{T,M,G} where {T,M<:DiscreteProductWeight,G<:ProductSubGrid}
name(m::DiscreteTensorSubWeight) = "Tensor of submeasures (supermeasure:"*name(supermeasure(m))
components(m::DiscreteTensorSubWeight) = map(_discretesubmeasure,components(points(m)),components(weights(supermeasure(m))))
component(m::DiscreteTensorSubWeight, i) = _discretesubmeasure(component(points(m),i),component(weights(supermeasure(m)),i))
weights(m::DiscreteTensorSubWeight) = OuterProductArray(map(weights, components(m))...)

function stencilarray(m::DiscreteTensorSubWeight)
    A = Any[]
    push!(A, component(m,1))
    for i = 2:length(components(m))
        push!(A," ⊗ ")
        push!(A, component(m,i))
    end
    A
end


text:
#=

# 201: 2D Laplace equation
([source code](SOURCE_URL))

=#

module Example201_Laplace2D

using VoronoiFVM,ExtendableGrids
using GridVisualize

## Flux function which describes the flux
## between neigboring control volumes
function g!(f,u0,edge)
    u=unknowns(edge,u0)
    f[1]=u[1,1]-u[1,2]
end

function main(;Plotter=nothing)
    nspecies=1
    ispec=1
    X=collect(0:0.2:1)
    grid=VoronoiFVM.Grid(X,X)
    physics=VoronoiFVM.Physics(num_species=nspecies,flux=g!)
    sys=VoronoiFVM.System(grid,physics)
    enable_species!(sys,ispec,[1])
    boundary_dirichlet!(sys,ispec,1,0.0)
    boundary_dirichlet!(sys,ispec,3,1.0)
    inival=unknowns(sys,inival=0)
    solution=unknowns(sys)
    solve!(solution,inival,sys)
    scalarplot(grid,solution[1,:],Plotter=Plotter)
    return solution[7]
end

## Called by unit test

function test()
    main() ≈ 0.2
end

end


text:
using Sads: CrossSection

"""
Rectangular river channel cross section.

- `W⃰`: bankfull width
- `Yₘ⃰`: bankfull depth
- `Yₘ`: flow depth
- `S0`: bed slope
- `n`: Manning roughness coefficient

See also: [`Dingman`](@ref)
"""
mutable struct Rectangular <: CrossSection
    W⃰:: Real
    Yₘ⃰:: Real
    Yₘ:: Real
    S0:: Real
    n:: Real
end

"""
Dingman river channel cross section.

Uses generalized expressions for cross-section geometry based on at-a-station hydraulic geometry relations.

- `W⃰`: bankfull width
- `Yₘ⃰`: bankfull depth
- `Yₘ`: flow depth
- `r`: geometry exponent
- `S0`: bed slope
- `n`: Manning roughness coefficient

See also: [`Rectangular`](@ref)
"""
mutable struct Dingman <: CrossSection
    W⃰:: Real
    Yₘ⃰:: Real
    Yₘ:: Real
    r:: Real
    S0:: Real
    n:: Real
end

"""
Calculate water surface width of cross section.
"""
function width(xs::Dingman)
    xs.W⃰ * (xs.Yₘ / xs.Yₘ⃰)^(1 / xs.r)
end

function width(xs::Rectangular)
    xs.W⃰
end

"""
Calculate average depth of cross section.
"""
function depth(xs::Dingman)
    (xs.r / (xs.r + 1)) * xs.Yₘ
end

function depth(xs::Rectangular)
    xs.Yₘ
end

export width, depth


text:
"""
Bridges `CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES, T}`
to `CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T`.
"""
struct GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T} <: MOIBC.AbstractBridge
    cons_domain::Vector{MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, CP.Domain{T}}}
    con_gc::MOI.ConstraintIndex{MOI.VectorAffineFunction{T}, CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}}
end

function MOIBC.bridge_constraint(
    ::Type{GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T}},
    model,
    f::MOI.VectorOfVariables,
    s::CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES, T},
) where {T}
    return MOIBC.bridge_constraint(
        GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T},
        model,
        MOI.VectorAffineFunction{T}(f),
        s,
    )
end

function MOIBC.bridge_constraint(
    ::Type{GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T}},
    model,
    f::MOI.VectorAffineFunction{T},
    s::CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES, T},
) where {T}
    f_scalars = MOIU.scalarize(f)

    cons_domain = MOI.ConstraintIndex{MOI.ScalarAffineFunction{T}, CP.Domain{T}}[
        MOI.add_constraint(
            model,
            f_scalars[i],
            CP.Domain(Set(s.values)),
            # Create one set object per constraint to avoid sharing any state.
        )
        for i in 1:s.dimension
    ]

    con_gc = MOI.add_constraint(
        model,
        f,
        CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}(s.dimension, copy(s.values)),
        # Copy the values for each constraint to avoid sharing any state.
    )

    return GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge(cons_domain, con_gc)
end

function MOI.supports_constraint(
    ::Type{GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T}},
    ::Union{Type{MOI.VectorOfVariables}, Type{MOI.VectorAffineFunction{T}}},
    ::Type{CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.CLOSED_COUNTED_VALUES, T}},
) where {T}
    return true
end

function MOIB.added_constrained_variable_types(::Type{GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T}}) where {T}
    return Tuple{Type}[]
end

function MOIB.added_constraint_types(::Type{GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T}}) where {T}
    return [
        (MOI.ScalarAffineFunction{T}, CP.Domain{T}),
        (MOI.VectorAffineFunction{T}, CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T}),
    ]
end

function MOI.get(
    b::GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.ScalarAffineFunction{T}, CP.Domain{T},
    },
) where {T}
    return length(b.cons_domain)
end

function MOI.get(
    ::GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T},
    ::MOI.NumberOfConstraints{
        MOI.VectorAffineFunction{T}, CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T},
    },
) where {T}
    return 1
end

function MOI.get(
    b::GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.ScalarAffineFunction{T}, CP.Domain{T},
    },
) where {T}
    return copy(b.cons_domain)
end

function MOI.get(
    b::GlobalCardinalityFixedClosed2GlobalCardinalityFixedOpenBridge{T},
    ::MOI.ListOfConstraintIndices{
        MOI.VectorAffineFunction{T}, CP.GlobalCardinality{CP.FIXED_COUNTED_VALUES, CP.OPEN_COUNTED_VALUES, T},
    },
) where {T}
    return [b.con_gc]
end


text:
module RocketNeverObservableTest

using Test
using Rocket

include("../test_helpers.jl")

@testset "NeverObservable" begin

    println("Testing: never")

    @testset begin
        @test never()    == NeverObservable{Any}()
        @test never(Int) == NeverObservable{Int}()
        @test never(Int) == never(Int)
        @test never(Int) != never(Float64)
    end

    @testset begin
        source = never(Int)
        io = IOBuffer()

        show(io, source)

        printed = String(take!(io))

        @test occursin("NeverObservable", printed)
        @test occursin(string(eltype(source)), printed)
    end

    run_testset([
        (
            source = never(),
            values = @ts()
        )
    ])

end

end


text:
# Unless noted otherwise, lengths and coordinates are in microns

# const variables w/ fixed types perform better than default (global) variables
"Simulation Controls: other setters in ~/src/walk_logic/runsim!.jl"
global NUMBER_OF_WALKS = 100
global MAX_STEPS_PER_WALK = 3*1485604
global FLOW_BIAS = true
global CATALASE_ON_WALLS = false

function set_FLOW_BIAS(trueorfalse::Bool=true)
    global FLOW_BIAS = trueorfalse
    return nothing
end

function set_CATALASE_ON_WALLS(trueorfalse::Bool=true)
    global CATALASE_ON_WALLS = trueorfalse
    return nothing
end

"Solvent Diffusion Step Sizes (microns per 2.7 microseconds)"
const SECONDS_PER_STEP = 0.0000027077
STEP_SIZE_DICT = Dict("water" => 0.1, "enz" => 0.005, "ppd" => 0.001)
const WATER_STEP_SIZE = 0.1
const ENZ_STEP_SIZE = 0.005
const PPD_STEP_SIZE = 0.001

"Sensor/Walls Coordinates"
const WALL_Y = 1.5
const SENSOR_HALF_WIDTH = 0.5 * 25
const SENSOR_SPACING = 20
const SENSOR_CENTER_MAX_X = SENSOR_HALF_WIDTH
const SENSOR_INNER_ADJ_MIN_X = SENSOR_HALF_WIDTH + SENSOR_SPACING
const SENSOR_INNER_ADJ_MAX_X = 3 * SENSOR_HALF_WIDTH + SENSOR_SPACING
const SENSOR_OUTER_ADJ_MIN_X = 3 * SENSOR_HALF_WIDTH + 2 * SENSOR_SPACING
const SENSOR_OUTER_ADJ_MAX_X = 5 * SENSOR_HALF_WIDTH + 2 * SENSOR_SPACING

"PPD Layer Dimensions"
const PPD_MAX_Y = 0.15
const PPD_MIN_Y = 0

"Enzyme Layer Dimensions"
const ENZYME_LEFT_X = -13.5
const ENZYME_RIGHT_X = 13.5
const ENZYME_MAX_Y_FROM_WALL = 0.15
const ENZYME_MAX_Y = WALL_Y + ENZYME_MAX_Y_FROM_WALL

"Spawn Inside Enzyme Layer"
const BORDER_CORRECTION = 0.001
const SPAWN_LEFT_X = ENZYME_LEFT_X + BORDER_CORRECTION
const SPAWN_RIGHT_X = ENZYME_RIGHT_X - BORDER_CORRECTION
const SPAWN_ENZYME_MAX_Y = ENZYME_MAX_Y - BORDER_CORRECTION

"Enzyme-Water Corners"
const CORNER_CUT_IN_ENZ = ENZ_STEP_SIZE * sqrt(2) / 2
const CORNER_CUT_IN_WATER = WATER_STEP_SIZE * sqrt(2) / 2
# test set: locationbools here
const WATER_TO_ENZ_NORTHEAST = [ENZYME_RIGHT_X - CORNER_CUT_IN_ENZ, ENZYME_MAX_Y - CORNER_CUT_IN_ENZ]
const WATER_TO_ENZ_NORTHWEST = [ENZYME_LEFT_X + CORNER_CUT_IN_ENZ, ENZYME_MAX_Y - CORNER_CUT_IN_ENZ]
const ENZ_TO_WATER_NORTHEAST = [ENZYME_RIGHT_X + CORNER_CUT_IN_WATER, ENZYME_MAX_Y + CORNER_CUT_IN_WATER]
const ENZ_TO_WATER_NORTHWEST = [ENZYME_LEFT_X - CORNER_CUT_IN_WATER, ENZYME_MAX_Y + CORNER_CUT_IN_WATER]

"Escape Bound Limits"
const ESCAPE_X = 5 * SENSOR_HALF_WIDTH + 3 * SENSOR_SPACING
const ESCAPE_Y = ESCAPE_X

"Safe Bound Limits"
const MARGIN_OF_COLLISION = 2 * WATER_STEP_SIZE
const SAFE_MIN_Y = ENZYME_MAX_Y + MARGIN_OF_COLLISION
const SAFE_MAX_Y = ESCAPE_Y - MARGIN_OF_COLLISION
const SAFE_MAX_X = ESCAPE_X - MARGIN_OF_COLLISION

max_stars_repo_name: s-y-wu/HMCResearchRandomWalks

max_issues_repo_name: s-y-wu/HMCResearchRandomWalks

max_forks_repo_name: s-y-wu/HMCResearchRandomWalks


text:
"""
    CartInd(i1, i2, is...)

Represents a "scalar" index of a multidimensional associative (`Assoc`/`Arr`).
"""
struct CartInd{N, T <: NTuple{N,Any}}
    inds::T
end
CartInd(inds::T) where {N, T <: NTuple{N, Any}} = CartInd{N, T}(inds)
CartInd{N}(inds::T) where {N, T <: NTuple{N, Any}} = CartInd{N, T}(inds)
CartInd(inds...) = CartInd(inds)

"""
    CartProd(indices1, indices2, ...)

Represents the Cartesian product of the index sets `indices1`, `indices2`, ect. The elements
are `CartInd(first(indices1), first(indices2), ...)` and so-on.
"""
struct CartProd{N, I <: CartInd{N}, T <: CartInd{N}, Inds <: NTuple{N,Any}} <: Associative{I, T}
    indices::Inds
end
# TODO constructors that `collect` on `KeyIterator`s and so-on?

# The indices of a Cartesian product are another Cartesian Product.
# e.g. if c = CartProd(([:a,:b,:c], 2:4)), then indices(c) = CartProd(OneTo(3), OneTo(3))
indices(c::CartProd) = CartProd(map(indices, c.indices))
IndexStyle(c::CartProd) = IndexDirect()
@propagate_inbounds function getindex(c::CartProd{N}, i::CartInd{N})
    CartInd(map(getindex, c.indices, i.inds))
end

# Eventually this falls back to a Cartesian range
start(c::CartProd{N, NTuple{N, OneTo}}) = map(start, c.indices)
function next(c::CartProd{N}, inds::NTuple{N, Any})

end
done(c::CartProd{N}, inds::NTuple{N, Any}) = all(map(done, c.indices, i))

# TODO CartSlice - may be a mixture of ranges and scalars, dimensions following APL rules