text:
using DDD
using Test, FileIO, JSON, SparseArrays, LinearAlgebra

cd(@__DIR__)
@testset "Compare loading input and output parameters" begin
    # Load and create.
    fileDislocationParameters = "./testData/sampleDislocationParameters.json"
    fileMaterialParameters = "./testData/sampleMaterialParameters.json"
    fileFEMParameters = "./testData/sampleFEMParameters.json"
    fileIntegrationParameters = "./testData/sampleIntegrationParameters.json"
    fileSlipSystem = "./testData/BCC.json"
    fileDislocationLoop = "./testData/samplePrismShear.json"
    dlnParams, matParams, femParams, intParams, slipSystems, dislocationLoop =
        loadParameters(
            fileDislocationParameters,
            fileMaterialParameters,
            fileFEMParameters,
            fileIntegrationParameters,
            fileSlipSystem,
            fileDislocationLoop,
        )

    missing, missing, missing, missing, missing, dislocationLoop2 = loadParameters(
        fileDislocationParameters,
        fileMaterialParameters,
        fileFEMParameters,
        fileIntegrationParameters,
        fileSlipSystem,
        "./testData/sampleDislocation.json",
    )

    @test typeof(dislocationLoop2) <: DislocationLoop

    regularCuboidMesh = buildMesh(matParams, femParams)
    surfNode = regularCuboidMesh.surfNode

    testGamma, testForceDisp = Boundaries(femParams, regularCuboidMesh)

    fileDislocationLoop = "./testData/samplePrismShear.json"
    fileIntegTime = "./testData/sampleIntegrationTime.json"
    integTime = loadIntegrationTime(fileIntegTime)
    network = DislocationNetwork(dislocationLoop, memBuffer = 1)
    # Dump simulationJSON.
    paramDumpJSON = "./testData/sampleDump.json"
    paramDumpJLD2 = "./testData/sampleDump.jld2"
    saveJSON(
        paramDumpJSON,
        dlnParams,
        matParams,
        femParams,
        intParams,
        slipSystems,
        dislocationLoop,
        integTime,
        testGamma,
        testForceDisp,
    )
    save(
        paramDumpJLD2,
        "dlnParams",
        dlnParams,
        "matParams",
        matParams,
        "femParams",
        femParams,
        "intParams",
        intParams,
        "slipSystems",
        slipSystems,
        "dislocationLoop",
        dislocationLoop,
        "integTime",
        integTime,
        "testForceDisp",
        testForceDisp,
    )
    networkDumpJSON = "./testData/sampleNetwork.json"
    networkDumpJLD2 = "./testData/sampleNetwork.jld2"
    saveJSON(networkDumpJSON, network)
    save(networkDumpJLD2, "network", network)
    # Reload simulationJSON.
    simulationJSON = loadJSON(paramDumpJSON)
    dlnParams2 = loadDislocationParameters(simulationJSON[1])
    matParams2 = loadMaterialParameters(simulationJSON[2])
    femParams2 = loadFEMParameters(simulationJSON[3])
    intParams2 = loadIntegrationParameters(simulationJSON[4])
    slipSystems2 = loadSlipSystem(simulationJSON[5])
    dislocationLoop2 = zeros(DislocationLoop, length(simulationJSON[6]))
    simulationJSON[6][1]
    for i in eachindex(dislocationLoop2)
        dislocationLoop2[i] = loadDislocationLoop(simulationJSON[6][i], slipSystems2)
    end
    testBoundaries2 = loadBoundaries(simulationJSON[8])
    testForceDisp2 = loadForceDisplacement(simulationJSON[9])
    network2 = loadNetwork(networkDumpJSON)
    integTime2 = loadIntegrationTime(simulationJSON[7])
    # Reload simulationJLD2.
    dlnParams3,
    matParams3,
    femParams3,
    intParams3,
    slipSystems3,
    dislocationLoop3,
    integTime3,
    testForceDisp3 = load(
        paramDumpJLD2,
        "dlnParams",
        "matParams",
        "femParams",
        "intParams",
        "slipSystems",
        "dislocationLoop",
        "integTime",
        "testForceDisp",
    )
    network3 = load(networkDumpJLD2, "network")
    # Ensure the data is all the same.
    @test compStruct(dlnParams, dlnParams2; verbose = true)
    @test compStruct(matParams, matParams2; verbose = true)
    @test compStruct(femParams, femParams2; verbose = true)
    @test compStruct(intParams, intParams2; verbose = true)
    @test compStruct(slipSystems, slipSystems2; verbose = true)
    @test compStruct(dislocationLoop, dislocationLoop2; verbose = true)
    @test compStruct(network, network2; verbose = true)
    @test compStruct(integTime, integTime2; verbose = true)
    @test compStruct(testForceDisp, testForceDisp2; verbose = true)
    for i in fieldnames(typeof(testBoundaries2))
        isnothing(getproperty(testBoundaries2, i)) ? continue : nothing
        @test isequal(getproperty(testBoundaries2, i), getproperty(testBoundaries2, i))
    end

    @test compStruct(dlnParams, dlnParams3; verbose = true)
    @test compStruct(matParams, matParams3; verbose = true)
    @test compStruct(femParams, femParams3; verbose = true)
    @test compStruct(intParams, intParams3; verbose = true)
    @test compStruct(slipSystems, slipSystems3; verbose = true)
    @test compStruct(dislocationLoop, dislocationLoop3; verbose = true)
    @test compStruct(network, network3; verbose = true)
    @test compStruct(integTime, integTime3; verbose = true)
    @test compStruct(testForceDisp, testForceDisp3; verbose = true)

    sampleRegCubMeshDump = "./testData/sampleRegCubMesh.jld2"
    regularCuboidMesh = buildMesh(matParams, femParams)
    save(sampleRegCubMeshDump, "mesh", regularCuboidMesh)
    regularCuboidMesh2 = load(sampleRegCubMeshDump, "mesh")
    @test compStruct(regularCuboidMesh, regularCuboidMesh2; verbose = true)
    @test JSON.lower(nodeTypeFE(1)) == 1

    rm(paramDumpJSON; force = true)
    rm(paramDumpJLD2; force = true)
    rm(networkDumpJSON; force = true)
    rm(networkDumpJLD2; force = true)
    rm(sampleRegCubMeshDump; force = true)
end

max_stars_repo_name: dcelisgarza/DDD

max_forks_repo_name: dcelisgarza/DDD


text:
using JuMP, Cbc
A = [ –1   1  3; 1   3 –7]
b = [–5;  10]
c = [ 1  2  5]
(m, n) = size(A)
modelo = Model(with_optimizer(Cbc.Optimizer))
@variable(modelo, x[1:n] >= 0)
@objective(modelo, Max, sum(c[j]*x[j] for j=1:n))
@constraints modelo begin
    c1[i=1:m], sum(A[i,j]*x[j] for j=1:n) <= b[i]
    c2, x[1] <= 10
end

print(modelo)
optimize!(modelo)
println(“Função objetivo: ”, objective_value(modelo))
for j = 1:n
    println(“x[$j] = ”, value(x[j]))
end

ext: jl

max_stars_repo_name: jvlsg/julia-tsp

max_issues_repo_name: jvlsg/julia-tsp

max_forks_repo_name: jvlsg/julia-tsp


text:
################################################################################
#
# Functional objectives for WRIdual (ε version)
#
################################################################################

export objWRIdual_ε!, objWRIdual_ε_yres!, objWRIdual_ε_yfixed!, objWRIdual_ε_mfixed!

function objWRIdual_ε!(F,
					   Gm::Union{Array{R, 1}, Nothing}, Gy::Union{Array{Array{C, 2}, 1}, Nothing},
					   m::Array{R, 2}, y::Union{Array{Array{C, 2}, 1}, Nothing},
					   h::NTuple{2, R},
					   freqs::Array{R, 1},
					   x_rcv::Array{R, 1}, z_rcv::R,
					   x_src::Array{R, 1}, z_src::R,
					   src_enc::Array{C, 2},
					   dat::Array{Array{C, 2}, 1},
					   ε::Union{R, Array{Array{R, 2}, 1}};
					   nabc::Union{Nothing, NTuple{4, Int64}} = nothing,
					   α::Union{Nothing, R, Array{C, 2}} = nothing,
					   weight_array::Array{R, 2},
					   gradmprec_fun = g->g,
	  				   bounds::Union{NTuple{2, R}, NTuple{2, Array{R, 2}}, Nothing} = nothing,
					   mode_y::String = "provided",
					   mode_ε::String = "provided")

    # Initializing objective value
    L = R(0)

    # Zero-ing gradient
    if Gm != nothing
        Gm .= R(0)
    end

	# Projecting m to box constrained set
	m = proj_bounds(m, bounds)

    # Updating obj value and gradient (loop over frequencies)
    n = size(m)
    Ns = size(src_enc, 2)
    Nf = length(freqs)
    for idxf = 1:Nf

		# Encoding data
		dat_enc = dat[idxf]*src_enc

         # Set absorbing boundary size (frequency-dependent)
        if nabc == nothing
            nabc = computeNabc(h, m, freqs[idxf])
        end

        # Setting source/receiver restriction/injection operators & source weigths
        Pr = restrOpMat(n, h, x_rcv, z_rcv, nabc)
        Ps = restrOpMat(n, h, x_src, z_src, nabc)
        Fsrc = Matrix(adjoint(Ps)*src_enc)

		# Normalization factors
		ηd = sqrt.(sum(abs.(dat_enc).^R(2), dims = 1))
		ηf = sqrt.(sum(abs.(Fsrc).^R(2), dims = 1))

        # Helmholtz operator
        H = lu(HelmholtzLinOp(h, m, freqs[idxf], nabc)) # LU factorization

		# Computing y, if not provided (y = data residual)
		if mode_y == "provided"
			y_ = y[idxf]
		elseif mode_y == "residual"
			Uinc = C.(H\Fsrc)
			y_ = dat_enc.-Pr*Uinc
		end

		# Solve for adjoint wavefield
		V = C.(adjoint(H)\(adjoint(Pr)*y_))

		# ||F^T*y||_W^2 = <W*F^T*y, F^T*y>
		norm_FTy2 = sum(abs.(V./weight_array).^R(2), dims = 1)

		# <y, r>
		y_dot_r = sum(y_.*conj.(dat_enc), dims = 1)-sum(V.*conj.(Fsrc), dims = 1)

		# ||y||
		norm_y = sqrt.(sum(abs.(y_).^R(2), dims = 1))

		# Computing ε
		if mode_ε == "adaptive"
			if isa(ε, R)
				ε_ = ε*norm_y./ηd
			elseif isa(ε, Array{Array{R, 2}, 1})
				ε_ = ε[idxf].*norm_y./ηd
			end
		elseif mode_ε == "provided"
			if isa(ε, R)
				ε_ = ε
			elseif isa(ε, Array{Array{R, 2}, 1})
				ε_ = ε[idxf]
			end
		end

		# Computing optimized α
		if α == nothing
			α_ = (abs.(y_dot_r) .>= ε_.*ηd.*norm_y).*(abs.(y_dot_r).-ε_.*ηd.*norm_y)./(ηf.^R(2).*norm_FTy2).*conj.(y_dot_r)./abs.(y_dot_r)
		else
			α_ = α
		end

		# Objective value update
        if F != nothing
			L += (sum(-R(0.5).*abs.(α_).^R(2).*ηf.^R(2).*norm_FTy2.+real.(α_.*y_dot_r).-ε_.*ηd.*abs.(α_).*norm_y, dims = 2)/(Nf*Ns))[1]
        end

		# Solve for wavefield (extended source)
        if Gm != nothing || Gy != nothing
			Q = Fsrc.+ηf.^R(2).*α_.*V./weight_array.^R(2)
			U = C.(H\Q)
		end

		# Gradient update (y)
        if Gy != nothing || (Gm != nothing && mode_y == "residual")
            Gy_ = (conj.(α_).*(dat_enc-Pr*U)-ε_.*abs.(α_).*ηd.*y_./norm_y)./(Nf*Ns)
        end
		if Gy != nothing
			Gy[idxf] = Gy_
		end

		# Gradient update (m)
		if Gm != nothing
            Gm_ = -(R(2)*π*freqs[idxf]).^R(2).*restrict(reshape(sum(real.(α_.*conj.(U).*V), dims = 2), sz_abc(n, nabc)), nabc)./(Nf*Ns)
			if mode_y == "residual"
				VGy = C.(adjoint(H)\(adjoint(Pr)*Gy_))
				Gm_ .+= -(R(2)*π*freqs[idxf]).^R(2).*restrict(reshape(sum(real.(conj.(Uinc).*VGy), dims = 2), sz_abc(n, nabc)), nabc)
			end
            Gm .+= gradmprec_fun(Gm_)
        end

    end

	if F != nothing
        return L
    end

end

function objWRIdual_ε_yres!(F, G, m, h, freqs, x_rcv, z_rcv, x_src, z_src, src_enc, dat, ε; nabc = nothing, α = nothing, weight_array, gradmprec_fun = g->g, bounds = nothing, mode_ε = "provided")

	return objWRIdual_ε!(F, G, nothing, m, nothing, h, freqs, x_rcv, z_rcv, x_src, z_src, src_enc, dat, ε; nabc = nabc, α = α, weight_array = weight_array, gradmprec_fun = gradmprec_fun, bounds = bounds, mode_y = "residual", mode_ε = mode_ε)

end

function objWRIdual_ε_yfixed!(F, Gm, m, y, h, freqs, x_rcv, z_rcv, x_src, z_src, src_enc, dat, ε; nabc = nothing, α = nothing, weight_array = weight_array, gradmprec_fun = g->g, bounds = nothing, mode_ε = "provided")

	return objWRIdual_ε!(F, Gm, nothing, m, y, h, freqs, x_rcv, z_rcv, x_src, z_src, src_enc, dat, ε; nabc = nabc, α = α, weight_array = weight_array, gradmprec_fun = gradmprec_fun, bounds = bounds, mode_y = "provided", mode_ε = mode_ε)

end

function objWRIdual_ε_mfixed!(F, Gy::Union{Array{C, 2}, Nothing}, m, y::Array{C, 2}, h, freq::R, x_rcv, z_rcv, x_src, z_src, src_enc, dat::Array{C, 2}, ε; nabc = nothing, α = nothing, weight_array = weight_array, bounds = nothing, mode_ε = "provided")

	if Gy == nothing
		return objWRIdual_ε!(F, nothing, nothing, m, [y], h, [freq], x_rcv, z_rcv, x_src, z_src, src_enc, [dat], ε; nabc = nabc, α = α, weight_array = weight_array, bounds = bounds, mode_y = "provided", mode_ε = mode_ε)
	else
		Gy_ = Array{Array{C, 2}, 1}(undef, 1)
		L = objWRIdual_ε!(F, nothing, Gy_, m, [y], h, [freq], x_rcv, z_rcv, x_src, z_src, src_enc, [dat], ε; nabc = nabc, α = α, weight_array = weight_array, bounds = bounds, mode_y = "provided", mode_ε = mode_ε)
		Gy .= -Gy_[1]
		return L
	end

end


text:
module MatrixNetworks

using Compat

if VERSION < v"0.4-"
    using Docile
end

"""
Module ``MatrixNetworks``: Documentation on the module

- Option 1: start with a sparse matrix A:
- example: ``M = MatrixNetwork(A)``
- Option 2: start with row and column indexes for the nonzeros in the matrix
- example: ``M = MatrixNetwork(ei,ej)``

Available functions: (use ?function_name to get more documentation)
- bfs
- bipartite_matching
- clustercoeffs
- corenums
- cosineknn
- csr_to_sparse
- dfs
- dijkstra
- dirclustercoeffs
- floydwarshall
- largest_component
- mst_prim
- scomponents
- sparse_to_csr

You can check the readme file here: \n
"https://github.com/nassarhuda/MatrixNetworks.jl/blob/master/README.md"
"""
MatrixNetworks

include("MatrixNetwork.jl")

function MatrixNetwork{T}(A::SparseMatrixCSC{T,Int64})
    At = A'
    return MatrixNetwork(size(At,2),At.colptr,At.rowval,At.nzval)
end

function MatrixNetwork(ei::Vector{Int64},ej::Vector{Int64})
    At = sparse(ej,ei,true);
    return MatrixNetwork(size(At,2),At.colptr,At.rowval,At.nzval)
end

include("scomponents.jl")
include("csr_to_sparse.jl")
include("sparse_to_csr.jl")
include("bipartite_matching.jl")
include("bfs.jl")
include("dfs.jl")
include("clustercoeffs.jl")
include("corenums.jl")
include("floydwarshall.jl")
include("manage_data.jl")
include("largest_component.jl")
include("cosineknn.jl")
include("dirclustercoeffs.jl")
include("dijkstra.jl")
include("mst_prim.jl")

# export everything to make them accessible as functions
export MatrixNetwork, bipartite_matching, edge_list, create_sparse,
bipartite_matching_setup, bipartite_matching_indicator, bfs,
dfs, clustercoeffs, corenums, scomponents, strong_components_map,
readSMAT, enrich, load_matrix_network, matrix_network_datasets,
csr_to_sparse, load_matrix_network_metadata, floydwarshall, largest_component,
sparse_to_csr, cosineknn, dirclustercoeffs, dijkstra, mst_prim, mst_prim_matrix,
csr_to_sparse_matrix, edge_indicator

end # end module


text:
"""
    read_file(file::String)
    read_file(path::Array{String,1}, file<:File; kwargs...)
    read_file(path::Array{String,1}, file::XLSXInput)
    read_file(path::String, x::T) where T<:File
    read_from(editor::T) where T<:Edit
This method to reads .csv mapping files required for editing. These files must be stored in
the `data/coremaps` directory. It returns a .csv file.

# Arguments
- `path::Array{String,1}` or `path::String`: Path to file *location*;
    does not include file name.
- `file::String`: Full path to file, including file name.
- `file<:File`: A SLiDE DataType used to store information about a file. Options include:
    - [`SLiDE.CSVInput`](@ref)
    - [`SLiDE.XLSXInput`](@ref)
- `editor<:Edit`: A SLiDE DataType used to store information about an edit to make in a
    DataFrame. Specifically, this function might be called for edit types that include the
    field `file` in reference to
    - [`SLiDE.Group`](@ref)
    - [`SLiDE.Map`](@ref)

# Keywords
- `shorten::Bool = false` or `shorten::Int`: if an integer length is specified, the
    DataFrame will be shortened to the input value. This is meant to aid troubleshooting
    during development.

# Returns
- `df::DataFrame`: If the input is a csv or xlsx file, this method will return a DataFrame.
- `yml::Dict{Any,Any}`: If the input file is a yaml file, this method will return a
    dictionary. All keys that correspond with SLiDE DataStream DataTypes will be converted
    to (lists of) those types.
"""
function read_file(path::Array{String,1}, file::GAMSInput; shorten = false)
    filepath = joinpath(SLIDE_DIR, path..., file.name)
    xf = readlines(filepath)
    df = gams_to_dataframe(xf; colnames = file.col)
    return df
end

function read_file(path::Array{String,1}, file::CSVInput; shorten = false)
    filepath = joinpath(SLIDE_DIR, path..., file.name)
    df = CSV.read(filepath, DataFrame; silencewarnings = true, ignoreemptylines = true, comment = "#", missingstrings = ["","\xc9","..."])
    NUMTEST = min(10, size(df,1))

    # A column name containing the word "Column" indicates that the input csv file was
    # missing a column header. Multiple (more than 2?) columns with missing header suggests
    # that the first row in the .csv was not the header, but rather, the file began with
    # comments. Here, we find the header and reread .csv into a DataFrame.
    # CSV.read() was used rather than converting the 2-D Array read using  dlmread()
    # because this was faster, and includes the option to ignore missing rows.
    if sum(Int.(occursin.("Column", string.(propertynames(df))))) > 2 || size(df)[2] == 1
        xf = DelimitedFiles.readdlm(filepath, ',', Any, '\n')
        xf = xf[1:NUMTEST,:]

        HEAD = findmax(sum.(collect(eachrow(Int.(length.(xf) .!= 0)))) .> 1)[2]
        df = CSV.read(filepath, DataFrame, silencewarnings = true, ignoreemptylines = true,
            missingstrings = ["","\xc9","..."]; header = HEAD)
    end

    # Remove footer rows. These are identified by rows at the bottom of the sheet that are
    # empty with the exception of the first column.
    if all(ismissing.(values(df[end,2:end])))
        x = sum.([Int.(ismissing.(values(row))) for row in eachrow(df[end-NUMTEST:end,:])])
        FOOT = size(df)[1] - (length(x) - (findmax(x)[2]-1))
        df = df[1:FOOT,:]
    end

    shorten != false && (df = df[1:shorten,:])
    return unique(df)
end

function read_file(path::Array{String,1}, file::XLSXInput; shorten = false)
    filepath = joinpath(SLIDE_DIR, path..., file.name)
    xf = XLSX.readdata(filepath, file.sheet, file.range)

    # Delete rows containing only missing values.
    xf = xf[[!all(row) for row in eachrow(ismissing.(xf))],:]
    df = DataFrame(xf[2:end,:], Symbol.(xf[1,:]), makeunique = true)

    shorten != false && (df = df[1:shorten,:])
    return df
end

function read_file(path::String, file::T; shorten = false) where T <: File
    return read_file([path], file; shorten = shorten)
end

# !!!! Is this method too niche?
function read_file(editor::T) where T <: Edit
    # !!!! Should we avoid including specific paths within functions?
    # !!!! Need to throw error if this is called when "file" is not a field.
    # DIR = abspath(joinpath(dirname(Base.find_package("SLiDE")), "..", "data", "coremaps"))
    DIR = joinpath("data", "coremaps")
    df = read_file(joinpath(DIR, editor.file))
    return df
end

function read_file(file::String; colnames = false)
    file = joinpath(SLIDE_DIR, file)

    if occursin(".map", file) | occursin(".set", file)
        return gams_to_dataframe(readlines(file); colnames = colnames)
    end

    if occursin(".yml", file) | occursin(".yaml", file)
        y = YAML.load(open(file))
        # Here, we first list all sub-subtypes of DataStream (DataTypes that are used in
        # editing datasource files). Then, we find where they overlap with keys in the
        # dictionary read from the YAML file.
        TYPES = string.([IU.subtypes.(IU.subtypes(DataStream))...;])
        KEYS = intersect(TYPES, collect(keys(y)))
        [y[k] = load_from(datatype(k), y[k]) for k in KEYS]
        return y

    elseif occursin(".csv", file)
        df = CSV.read(file, DataFrame, silencewarnings = true, ignoreemptylines=true, comment = "#",
            missingstrings = ["","\xc9","..."])
        return df
    end
end

"""
    load_from(::Type{T}, df::DataFrame) where T <: Any
Load a DataFrame `df` into a structure of type T.

!!! note

    This requires that all structure fieldnames are also DataFrame column names.
    Extra dataframe columns are acceptable, although that information will not be used.

# Arguments
- `::Type{T} where T <: Any`: Any DataType.
- `df::DataFrame`: The DataFrame storing the information to store as a DataType.

# Returns
- `x<:Any`: The DataType specified as an argument.
- `lst::Array{T} where T<:Any`: A list of elements of the DataType specified as an argument
    given a multi-row DataFrame.

# Example

```julia
df = DataFrame(from = ["State"], to = ["region"])
load_from(Rename, df)
```
"""
function load_from(::Type{T}, d::Array{Dict{Any,Any},1}) where T <: Any
    lst = if all(isarray.(T.types))
        ensurearray(load_from(T, convert_type(DataFrame, d)))
    else
        vcat(ensurearray(load_from.(T, d))...)
    end
    return size(lst)[1] == 1 ? lst[1] : lst
end

function load_from(::Type{T}, d::Dict{Any,Any}) where T <: Any
    # Fill the datatype with the values in the dictionary keys, ensuring correct t.
    (fields, types) = (string.(fieldnames(T)), T.types)
    d = _load_path(d)

    # Fill the datatype with the input.
    # if length(unique(isarray.(types))) == 2
    if any(isarray.(types)) && !all(isarray.(types))
        # Restructure data into a list of inputs in the order and type required when
        # creating the datatype. Ensure that all array entries should, in fact, be arrays.
        inps = [_load_as_type(T, d[f], t) for (f,t) in zip(fields, types)]
        inpscorrect = isarray.(inps) .== isarray.(types)

        # If all inputs are of the correct structure, fill the data type.
        if all(inpscorrect)
            lst = [T(inps...)]
        # If some inputs are arrays when they shouldn't be, expand these into a new list of
        # dictionaries to create a list of datatypes, including all array values.
        # First, create a dictionary determining whether the entry needs to be split.
        # Then, split the dictionary into a list of arrays where necessary.
        else
            LEN = length(inps[findmax(.!inpscorrect)[2]])
            splitarray = Dict(fields[ii] => !inpscorrect[ii] for ii in 1:length(inps))
            lst = [Dict{Any,Any}(k => splitarray[k] ? d[k][ii] : d[k] for k in keys(d))
                for ii in 1:LEN]
            lst = ensurearray(load_from(T, lst))
        end
    else
        lst = ensurearray(load_from(T, convert_type(DataFrame, d)))
    end
    return size(lst)[1] == 1 ? lst[1] : lst
end

function load_from(::Type{T}, df::DataFrame) where T <: Any
    (fields, types) = (fieldnames(T), T.types)

    # Print warning if DataFrame is missing required columns.
    missing_fields = setdiff(fields, propertynames(df))
    if length(missing_fields) > 0
        @warn(string("DataFrame columns missing required fields to fill DataType ", Rename),
            missing_fields)
    end

    # If one of the struct fields is an ARRAY, we here assume that it is the length of the
    # entire DataFrame, and all other fields are duplicates.
    if any(isarray.(T.types))
        inps = [_load_as_type(T, df[:, f], t) for (f,t) in zip(fields, types)]
        lst = [T(inps...)]
    # If each row in the input df fills one and only one struct,
    # create a list of structures from each DataFrame row.
    else
        lst = [T((_load_as_type(T, row[f], t) for (f,t) in zip(fields, types))...)
            for row in eachrow(df)]
    end
    return size(lst)[1] == 1 ? lst[1] : lst
end

"""
    _load_path(d::Dict)
Edits directories containing a list of directories ending in a file name as one path.
"""
function _load_path(d::Dict)
    FILES = [".csv", ".xlsx", ".txt", ".map", ".set"]
    for (k, lst) in d
        if typeof(lst) .== Array{String,1}
            ii_file = [any(occursin.(FILES, x)) for x in lst]
            (ii_file[end] && .!any(ii_file[1:end-1])) && (d[k] = joinpath(lst...))
        end
    end
    return d
end

"""
    _load_as_type(::Type{Any}, entry, type::DataType)
Converts an entry to the required DataType
"""
function _load_as_type(entry, type::DataType)
    entry = ensurearray(convert_type.(type, entry))
    (!isarray(type) && length(entry) == 1) && (entry = entry[1])
    return entry
end

_load_as_type(::Type{T}, entry, type::DataType) where T<:Any = _load_as_type(entry, type)
_load_as_type(::Type{Drop},    entry, type::Type{Any})    = _load_as_type(_load_case(entry), type)
_load_as_type(::Type{Rename},  entry, type::Type{Symbol}) = _load_as_type(_load_case(entry), type)
_load_as_type(::Type{Replace}, entry, type::Type{Any})    = _load_as_type(_load_case(entry), type)
_load_as_type(::Type{Operate}, entry, type::Type{Symbol}) = _load_as_type(_load_axis(entry), type)

"""
    _load_case(entry::AbstractString)
Standardizes string identifiers that indicate a case change (upper-to-lower or vice-versa)
for easier editing.
"""
function _load_case(entry::AbstractString)
    test = lowercase(entry)

    occursin("lower", test) && (entry = "lower")
    occursin("upper", test) && (entry = "uppercasefirst" == test ? "uppercasefirst" : "upper")
    occursin("titlecase", test) && (entry = "titlecase")

    "all" == test    && (entry = "all")
    "unique" == test && (entry = "unique")
    return entry
end

_load_case(entry::Any) = entry

"""
    _load_axis(entry::Any)
"""
function _load_axis(entry::AbstractString)
    entry = convert_type(String, entry)
    ("1" == entry || occursin("row", lowercase(entry))) && (entry = "row")
    ("2" == entry || occursin("col", lowercase(entry))) && (entry = "col")
    return entry
end

_load_axis(entry::Any) = _load_axis(convert_type.(String, entry))

"""
    write_yaml(path, file::XLSXInput)
This function reads an XLSX file and writes a new yaml file containing the information in
each spreadsheet column. Sheet propertynames in the XLSX file correspond to the directory where new
files will be printed (`path/file.sheet/`). Yaml files will be named after the text in the
column's first row.

# Arguments
- `path::String` or `path::Array{String,1}`: XLSX file location. New yaml files will be
    printed here, as well.
- `file::XLSXInput` or `files::Array{XLSXInput,1}`: XLSX file information (file name, sheet name, sheet range) or list of multiple sheets.

# Returns
- `filenames::Array{String,1}`: List of yaml files

"""
function write_yaml(path, file::XLSXInput)
    println("\nGENERATE YAML FILES FROM ", file.name)

    # List all key words in the yaml file and use to add (purely aesthetic) spacing.
    KEYS = string.([IU.subtypes.(IU.subtypes(DataStream))...; "PathIn"], ":")

    # Read the XLSX file, identify relevant (not "missing"), and generate list of resultant
    # yaml file propertynames.
    println("  Reading ", file.name)
    df_all = read_file(path, file)
    df_all = df_all[:, .!occursin.(:missing, propertynames(df_all))]

    # Make sure the path exists and save the propertynames of yaml propertynames to generaate.
    path = joinpath(SLIDE_DIR, path, file.sheet)
    !isdir(path) && mkpath(path)
    filenames = joinpath.(path, string.(propertynames(df_all), ".yml"))

    # Iterate through columns. For each column, create a new yaml file and fill it with
    # the column text. Mark the yaml file as "Autogenerated" and ensure one space only
    # between each input file element that corresponds with a SLiDE Datastream subtype.
    for COL in 1:size(df_all,2)
        println("  Generating ", filenames[COL])

        # Remove all lines that do not contain text.
        lines = dropmissing(df_all, COL)[:,COL]
        lines = lines[match.(r"\S.*", lines) .!= nothing]

        open(filenames[COL], "w") do f
            println(f, string("# Autogenerated from ", file.name))
            for line in lines
                any(occursin.(KEYS, line)) && (println(f, ""))
                println(f, line)
            end
        end
    end
    return filenames
end

write_yaml(path::Array{String,1}, file::XLSXInput) = write_yaml(joinpath(path...), file)
write_yaml(path, files::Array{XLSXInput,1}) = vcat([write_yaml(path, f) for f in files]...)

"""
    run_yaml(filename::String)
    run_yaml(filenames::Array{String,1})
This file runs (a) yaml file(s) if it includes the line `Editable: true`. If given a list of
input files, the function will print a list of yaml files that were not marked as editable.
A file might not be "editable" if SLiDE functionality cannot make all of the specified edits.

# Arguments
- `filename::String` or `filenames::Array{String,1}`: yaml file name (or list of propertynames) to run

# Return
- `filename::String` or `filenames::Array{String,1}`: yaml file name (or list of propertynames) that
    was/were not ran by the function because they were annotated with `Editable: false`
"""
function run_yaml(filename::String)

    println("\n", filename)
    println("  Reading yaml file...")

    y = read_file(filename)

    if haskey(y, "Editable") && y["Editable"]
        # Iteratively make all edits and update output data file in accordance with
        # operations defined in the input YAML file.
        println("  Parsing and standardizing...")
        df = unique(edit_with(y))

        # Create the path where the output file will go if it doesn't already exist.
        path = joinpath(SLIDE_DIR, ensurearray(y["PathOut"])[1:end-1]...)
        !isdir(path) && mkpath(path)

        println("  Writing to ", joinpath(ensurearray(y["PathOut"])...))
        CSV.write(joinpath(SLIDE_DIR, ensurearray(y["PathOut"])...), df)
        return nothing
    else
        println("  Skipping... (not editable)")
        return filename
    end
end

function run_yaml(filenames::Array{String,1})
    println("\nSTANDARDIZE DATA FILES:")
    filenames = [run_yaml(f) for f in filenames]
    filenames = filenames[filenames .!== nothing]
    if length(filenames) > 0
        @warn(string("run_yaml() generated no output for:", string.("\n  ", filenames)...,
            "\nAdd \"Editable: true\" to yaml file to run automatically."))
    end
    return filenames
end

"""
    gams_to_dataframe(xf::Array{String,1}; colpropertynames = false)
This function converts a GAMS map or set to a DataFrame, expanding sets into multiple rows.

# Arguments
- `xf::Array{String,1}`: A list of rows of text from a .map or a .set input file.

# Keywords
- `colpropertynames = false`: A user has the option to specify the column propertynames of the output
    DataFrame. If none are specified, a default of `[missing, missing_1, ...]` will be used,
    consistent with the default column headers for `CSV.read()` if column propertynames are missing.

# Returns
- `df::DataFrame`: A dataframe representation of the GAMS map or set.
"""
function gams_to_dataframe(xf::Array{String,1}; colpropertynames = false)
    # Convert the input array into a DataFrame and use SLiDE editing capabilities to split
    # each rows into columns, based on the syntax.
    df = DataFrame(missing = xf)
    df = edit_with(df, Match(Regex("^(?<missing>\\S+)\\.(?<missing_1>[\\S^,]*)\\s*\"*(?<missing_2>[^\"]*),?"),
        :missing, [:missing, :missing_1, :missing_2]))
    ROWS, COLS = size(df)

    # Does the DataFrame row contain a set (indicated by parentheses)?
    df_set = match.(r"^\((.*)\)", df)
    df_isset = df_set .!== nothing

    # If so, expand into multiple rows by converting the row into a dictionary -- with
    # column propertynames as keys and the set divided into a list -- and back into a DataFrame.
    df = [[DataFrame(Dict(k => df_isset[ii,k] ? string.(split(df_set[ii,k][1], ",")) : df[ii,k]
        for k in propertynames(df))) for ii in 1:ROWS]...;]

    # If the user specified column propertynames, apply those here and
    # return a DataFrame sorted based on the mapping values.
    colpropertynames != false && (df = edit_with(df, Rename.(propertynames(df), colpropertynames)))
    return COLS > 1 ? sort(df, reverse(propertynames(df)[1:2])) : sort(df)
end

max_stars_repo_name: matbesancon/SLiDE

max_issues_repo_name: matbesancon/SLiDE

max_forks_repo_name: matbesancon/SLiDE


text:
# testing flag demography with constant fertility

# BASIC PARAMETER
par = (
    S  = 20,
    Sr = 10
)

# a collection to save results

#


text:
using CompScienceMeshes
using LinearAlgebra

o = point(3,4,2)
a = point(4,2,4)
b = point(2,8,2)
c = point(0,4,6)

tet = simplex(o,a,b,c)
qw = quadpoints(tet,4)
print(volume(tet),"==",abs(dot(a-o,cross(b-o,c-o)))/6,"==",abs(sum(w for (p,w) in qw)),"\n")

o,x,y,z = CompScienceMeshes.euclidianbasis(3)

tet = simplex(o,x,2y,3z)
volume(tet)

qw = quadpoints(tet,4)
length(qw)

print(abs(sum(w for (p,w) in qw)-volume(tet)) < 0.001,"\n")
print(sum(w*sin(p[1])*p[2]^2*p[3]^3 for (p,w) in qw))
#https://www.wolframalpha.com/input/?i=int_0^1+int_0^(2-2x)+int_0^(3-3x-1.5y)+sin(x)y^2z^3+dz+dy+dx


text:
# Contiguous on-device arrays

export ROCDeviceArray, ROCDeviceVector, ROCDeviceMatrix, ROCBoundsError

# construction

"""
    ROCDeviceArray(dims, ptr)
    ROCDeviceArray{T}(dims, ptr)
    ROCDeviceArray{T,A}(dims, ptr)
    ROCDeviceArray{T,A,N}(dims, ptr)

Construct an `N`-dimensional dense ROC device array with element type `T`
wrapping a pointer, where `N` is determined from the length of `dims` and `T`
is determined from the type of `ptr`. `dims` may be a single scalar, or a
tuple of integers corresponding to the lengths in each dimension). If the rank
`N` is supplied explicitly as in `Array{T,N}(dims)`, then it must match the
length of `dims`. The same applies to the element type `T`, which should match
the type of the pointer `ptr`.
"""
ROCDeviceArray

# NOTE: we can't support the typical `tuple or series of integer` style
# construction, because we're currently requiring a trailing pointer argument.

struct ROCDeviceArray{T,N,A} <: AbstractArray{T,N}
    shape::Dims{N}
    ptr::DevicePtr{T,A}

    # inner constructors, fully parameterized, exact types (ie. Int not <:Integer)
    ROCDeviceArray{T,N,A}(shape::Dims{N}, ptr::DevicePtr{T,A}) where {T,A,N} = new(shape,ptr)
end

const ROCDeviceVector = ROCDeviceArray{T,1,A} where {T,A}
const ROCDeviceMatrix = ROCDeviceArray{T,2,A} where {T,A}

# outer constructors, non-parameterized
ROCDeviceArray(dims::NTuple{N,<:Integer}, p::DevicePtr{T,A}) where {T,A,N} = ROCDeviceArray{T,N,A}(dims, p)
ROCDeviceArray(len::Integer,              p::DevicePtr{T,A}) where {T,A}   = ROCDeviceVector{T,A}((len,), p)

# outer constructors, partially parameterized
ROCDeviceArray{T}(dims::NTuple{N,<:Integer},   p::DevicePtr{T,A}) where {T,A,N} = ROCDeviceArray{T,N,A}(dims, p)
ROCDeviceArray{T}(len::Integer,                p::DevicePtr{T,A}) where {T,A}   = ROCDeviceVector{T,A}((len,), p)
ROCDeviceArray{T,N}(dims::NTuple{N,<:Integer}, p::DevicePtr{T,A}) where {T,A,N} = ROCDeviceArray{T,N,A}(dims, p)
ROCDeviceVector{T}(len::Integer,               p::DevicePtr{T,A}) where {T,A}   = ROCDeviceVector{T,A}((len,), p)

# outer constructors, fully parameterized
ROCDeviceArray{T,N,A}(dims::NTuple{N,<:Integer}, p::DevicePtr{T,A}) where {T,A,N} = ROCDeviceArray{T,N,A}(Int.(dims), p)
ROCDeviceVector{T,A}(len::Integer,               p::DevicePtr{T,A}) where {T,A}   = ROCDeviceVector{T,A}((Int(len),), p)

# getters

Base.pointer(a::ROCDeviceArray) = a.ptr
Base.pointer(a::ROCDeviceArray, i::Integer) =
    pointer(a) + (i - 1) * Base.elsize(a)

Base.elsize(::Type{<:ROCDeviceArray{T}}) where {T} = sizeof(T)
Base.size(g::ROCDeviceArray) = g.shape
Base.length(g::ROCDeviceArray) = prod(g.shape)

# conversions

Base.unsafe_convert(::Type{DevicePtr{T,A}}, a::ROCDeviceArray{T,N,A}) where {T,A,N} = pointer(a)

# indexing
# FIXME: Boundschecking

@inline function Base.getindex(A::ROCDeviceArray{T}, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = datatype_align(T)
    Base.unsafe_load(pointer(A), index, Val(align))::T
end

@inline function Base.setindex!(A::ROCDeviceArray{T}, x, index::Integer) where {T}
    @boundscheck checkbounds(A, index)
    align = datatype_align(T)
    Base.unsafe_store!(pointer(A), x, index, Val(align))
    return A
end

# other

Base.show(io::IO, a::ROCDeviceVector) =
    print(io, "$(length(a))-element device array at $(pointer(a))")
Base.show(io::IO, a::ROCDeviceArray) =
    print(io, "$(join(a.shape, '×')) device array at $(pointer(a))")

Base.show(io::IO, mime::MIME"text/plain", a::ROCDeviceArray) = show(io, a)

@inline function Base.unsafe_view(A::ROCDeviceVector{T}, I::Vararg{Base.ViewIndex,1}) where {T}
    ptr = pointer(A) + (I[1].start-1)*sizeof(T)
    len = I[1].stop - I[1].start + 1

    return ROCDeviceArray(len, ptr)
end

@inline function Base.iterate(A::ROCDeviceArray, i=1)
    if (i % UInt) - 1 < length(A)
        (@inbounds A[i], i + 1)
    else
        nothing
    end
end


text:
@windows? (
           begin
               # location of cudart64_xx.dll or cudart32_xx.dll have to be in PATH env var
               # ex: C:\Program Files\NVIDIA GPU Computing Toolkit\CUDA\v6.5\bin
               # (by default, it is done by CUDA toolkit installer)

               const dllnames = (WORD_SIZE==64) ?
                 ["curand64_70", "curand64_65", "curand64_60", "curand64_55", "curand64_50", "curand64_50_35"] :
                 ["curand32_70", "curand32_65", "curand32_60", "curand32_55", "curand32_50", "curand32_50_35"]
               const libcurand = Libdl.find_library(dllnames, [""])
           end
           : # linux or mac
           begin
               const libcurand = Libdl.find_library(["libcurand"],
                                                    ["/usr/local/cuda/lib", "/usr/local/cuda/lib64"])
           end)

if isempty(libcurand)
    error("CURAND library cannot be found")
end


text:
module InfiniteArrays
using Base, Statistics, LinearAlgebra, FillArrays, LazyArrays, DSP

import Base: *, +, -, /, \, ==, isinf, isfinite, sign, angle, show, isless,
            fld, cld, div, min, max, minimum, maximum, mod,
            <, ≤, >, ≥, promote_rule, convert, unsafe_convert, copy,
            size, step, isempty, length, first, last,
            getindex, setindex!, intersect, @_inline_meta,
            sort, sort!, issorted, sortperm, sum, in, broadcast,
            eltype, elsize, parent, parentindices, reinterpret,
            unaliascopy, dataids,
            real, imag, conj, transpose,
            exp, log, sqrt, cos, sin, tan, csc, sec, cot,
            cosh, sinh, tanh, csch, sech, coth, acos, asin, atan, acsc, asec, acot,
            acosh, asinh, atanh, acsch, asech, acoth, (:),
            AbstractMatrix, AbstractArray, checkindex, unsafe_length, OneTo,
            to_shape, _sub2ind, print_matrix, print_matrix_row, print_matrix_vdots,
            checkindex, Slice, @propagate_inbounds, @_propagate_inbounds_meta,
         	_in_range, _range, _rangestyle, Ordered,
         	ArithmeticWraps, floatrange, reverse, unitrange_last,
         	AbstractArray, AbstractVector, Array, Vector, Matrix,
         	axes, (:), _sub2ind_recurse, broadcast, promote_eltypeof,
         	diff, cumsum, show_delim_array, show_circular, Int,
         	similar, _unsafe_getindex, string, zeros, fill, permutedims,
         	cat_similar, vcat,
		 	reshape, ReshapedIndex, ind2sub_rs, _unsafe_getindex_rs,
         	searchsorted, Ordering, lt

using Base.Broadcast
import Base.Broadcast: BroadcastStyle, AbstractArrayStyle, Broadcasted, broadcasted,
                        @nexprs, @ncall, combine_eltypes, DefaultArrayStyle, instantiate, axistype

import LinearAlgebra: BlasInt, BlasFloat, norm, diag, diagm, ishermitian, issymmetric,
                             det, logdet, istriu, istril, adjoint, tr, AbstractTriangular,
                             norm2, norm1, normp

import Statistics: mean, median

import FillArrays: AbstractFill, getindex_value, fill_reshape
import LazyArrays: LazyArrayStyle, AbstractBandedLayout, MemoryLayout, LazyLayout, UnknownLayout,
                    ZerosLayout, @lazymul, AbstractArrayApplyStyle, CachedArray, CachedVector,
                    reshapedlayout

import DSP: conv

export ∞, Hcat, Vcat, Zeros, Ones, Fill, Eye, BroadcastArray, cache

include("Infinity.jl")
include("infrange.jl")
include("infarrays.jl")
include("reshapedarray.jl")

##
# Fill FillArrays
##

@lazymul Ones{<:Any,1,Tuple{OneToInf{Int}}}
@lazymul Fill{<:Any,1,Tuple{OneToInf{Int}}}
@lazymul Zeros{<:Any,1,Tuple{OneToInf{Int}}}

@lazymul Ones{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}
@lazymul Ones{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}
@lazymul Ones{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}

@lazymul Fill{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}
@lazymul Fill{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}
@lazymul Fill{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}

@lazymul Zeros{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}
@lazymul Zeros{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}
@lazymul Zeros{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}

length(::Ones{<:Any,1,Tuple{OneToInf{Int}}}) = ∞
length(::Fill{<:Any,1,Tuple{OneToInf{Int}}}) = ∞
length(::Zeros{<:Any,1,Tuple{OneToInf{Int}}}) = ∞
length(::Ones{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}) = ∞
length(::Ones{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}) = ∞
length(::Ones{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}) = ∞
length(::Fill{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}) = ∞
length(::Fill{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}) = ∞
length(::Fill{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}) = ∞
length(::Zeros{<:Any,2,Tuple{OneToInf{Int},OneToInf{Int}}}) = ∞
length(::Zeros{<:Any,2,<:Tuple{OneToInf{Int},<:Any}}) = ∞
length(::Zeros{<:Any,2,<:Tuple{<:Any,OneToInf{Int}}}) = ∞

for op in (:norm2, :norm1)
   @eval $op(a::Zeros{T,N,NTuple{N,OneToInf{Int}}}) where {T,N} = norm(getindex_value(a))
end

normp(a::Zeros{T,N,NTuple{N,OneToInf{Int}}}, p) where {T,N} = norm(getindex_value(a))

for N=1:3
   for op in (:norm2, :norm1)
      @eval function $op(a::AbstractFill{T,$N,NTuple{$N,OneToInf{Int}}}) where {T,N}
         z = norm(getindex_value(a))
         iszero(z) && return z
         typeof(z)(Inf)
      end
   end
   @eval function normp(a::AbstractFill{T,$N,NTuple{$N,OneToInf{Int}}}, p) where {T,N }
      z = norm(getindex_value(a))
      iszero(z) && return z
      typeof(z)(Inf)
   end
end

for Typ in (:Number, :AbstractVector)
   @eval begin
      vcat(a::$Typ, b::AbstractFill{<:Any,1,<:Tuple{<:OneToInf}}) = Vcat(a, b)
      vcat(a::$Typ, c::CachedVector{<:Any,<:Any,<:AbstractFill{<:Any,1,<:Tuple{<:OneToInf}}}) =
         CachedArray(vcat(a, view(c.data,1:c.datasize[1])), c.array)
   end
end

# cat_similar(A, T, ::Tuple{Infinity}) = zeros(T, ∞)

##
# Temporary hacks for base support
##
OneTo(::Infinity) = OneToInf()
function OneTo(x::OrientedInfinity)
    iszero(x.angle) && return OneTo(∞)
    throw(ArgumentError("Cannot create infinite range with negative length"))
end
OneTo{T}(::Infinity) where T<:Integer = OneToInf{T}()
UnitRange(start::Integer, ::Infinity) = InfUnitRange(start)
UnitRange{T}(start::Integer, ::Infinity) where T<:Real = InfUnitRange{T}(start)
OneTo(a::OneToInf) = a
OneTo{T}(::OneToInf) where T<:Integer = OneToInf{T}()

Int(::Infinity) = ∞

axistype(::OneTo{T}, ::OneToInf{V}) where {T,V} = OneToInf{promote_type(T,V)}()
axistype(::OneToInf{V}, ::OneTo{T}) where {T,V} = OneToInf{promote_type(T,V)}()

# sort.jl
# returns the range of indices of v equal to x
# if v does not contain x, returns a 0-length range
# indicating the insertion point of x
function searchsorted(v::AbstractVector, x, ilo::Int, ::Infinity, o::Ordering)
    lo = ilo-1
    hi = ∞
    @inbounds while lo < hi-1
        m = isinf(hi) ? lo + 1000 : (lo+hi)>>>1
        if lt(o, v[m], x)
            lo = m
        elseif lt(o, x, v[m])
            hi = m
        else
            a = searchsortedfirst(v, x, max(lo,ilo), m, o)
            b = searchsortedlast(v, x, m, hi, o)
            return a : b
        end
    end
    return (lo + 1) : (hi - 1)
end

end # module


text:
using IndexedTables
using NamedTuples
using OnlineStats
using DataValues
import DataValues: NA
using Base.Test

@testset "IndexedTables" begin

include("test_core.jl")
include("test_utils.jl")
include("test_tabletraits.jl")
include("test_collect.jl")

end


text:
function distance(H::Array{Int64, 2}, H_sim::Array{Int64, 2})

	dH::Array{Float64, 2} = convert(Array{Float64, 2}, H - H_sim)
	cum_dH::Array{Float64, 2} = cumsum(cumsum(dH, 1), 2)
	d::Float64 = log10( sum( cum_dH.^2 ) )

	return d
end


text:
using Test, ICOADSDict

"Function that returns a valid file name of a non-existing file for all OSs"
function tempname2()
    f, h = mktemp()
    close(h)
    rm(f)
    f
end
tlist = [Int8, UInt8, Int16, UInt16, Int32,UInt32, Int64,UInt64, Float32, Float64]

@testset "lineparsing" begin
    include("oneliner.jl")
end
nothing


text:
################################################################################
# The following data structures are the basis of OpenMendel.
# Complex data types allow related variables to be accessed as a whole.
# The plural of a variable indicates the number of instances of the variable.
# Consider removing those marked with <remove>.
################################################################################

export Person, Pedigree, NuclearFamily, Locus, SnpDataStruct

mutable struct Person
  people :: Int # number of people in sample
  populations :: Int # number of populations in sample
  name :: Vector{AbstractString} # <remove>
  pedigree :: Vector{Int} # <remove>
  mother :: Vector{Int} # index of mother
  father :: Vector{Int} # index of father
  male :: Vector{Bool} # true for male; false for female
  next_twin :: Vector{Int} # pointer to next twin
  primary_twin :: Vector{Int} # pointer to primary twin
  admixture :: Matrix{Float64} # person.admixture[person, population]
  children :: Vector{BitSet} # set of children
  spouse :: Vector{BitSet} # set of spouses
  genotype :: Matrix{Set{Tuple{Int, Int}}} # person.genotype[person, locus]
  homozygotes :: Matrix{Int} # person.homozygotes[person, locus]
  disease_status :: Vector{AbstractString} # <remove>
  variable :: Matrix{Float64} # person.variable[person, variable] <remove>
  variable_name :: Vector{AbstractString} # <remove>
end # Person

mutable struct Pedigree
  pedigrees :: Int # number of pedigrees in sample
  name :: Vector{AbstractString} # <remove>
  pedstart :: Vector{Int} # start of pedigree
  twin_finish :: Vector{Int} # end of pedigree including co-twins
  pedfinish :: Vector{Int} # end of pedigree excluding co-twins
  individuals :: Vector{Int} # count per pedigree including co-twins
  founders :: Vector{Int} # count per pedigree
  females :: Vector{Int} # count per pedigree
  males :: Vector{Int} # count per pedigree
  twins :: Vector{Int} # count per pedigree
  families :: Vector{Int} # count of nuclear families per pedigree
  loglikelihood :: Matrix{Float64} # pedigree.loglikelihood[pedigree, 1:2]
end # Pedigree

mutable struct NuclearFamily
  families :: Int # number of nuclear families across all pedigrees
  pedigree :: Vector{Int} # pedigree of nuclear family
  mother:: Vector{Int} # mother of nuclear family
  father :: Vector{Int} # father of nuclear family
  sib :: Vector{BitSet} # sib set of nuclear family
end # NuclearFamily

mutable struct Locus
  loci :: Int # number of loci in locus file
  model_loci :: Int # number of loci in current model
  trait :: Int # position of the trait locus among the model loci
  free_recombination :: Bool # true when model loci unlinked
##
  mutable_locus :: Int # locus subject to mutation
  donor_allele :: Int # allele subject to mutation
  recipient_allele :: Int # allele mutated to
  mutation_rate :: Vector{Float64} # female and male mutation rates
  name :: Vector{AbstractString} # <remove>
  chromosome :: Vector{AbstractString} # <remove>
  base_pairs :: Vector{Int} # <remove>
  morgans :: Matrix{Float64} # female and male genetic map distance <remove ?>
  theta :: Matrix{Float64} # model loci recombination fractions
  xlinked :: Vector{Bool} # true for xlinked; false for autosomal
  alleles :: Vector{Int} # number of alleles per locus
  allele_name :: Vector{Vector{AbstractString}} # locus.allele_name[loc][allele]
  frequency :: Vector{Matrix{Float64}} # locus.frequency[loc][population,allele]
  lumped_frequency :: Array{Float64, 3} # indices: pedigree, population, locus
  model_locus :: Vector{Int} # indices of loci currently modeled
  locus_field_in_person_frame :: Vector{Int}
end # Locus

#
# Old-style structure for SNP and person information.
# (Should be better coordinated with SnpArrays.)
#
mutable struct SnpDataStruct
  people::Int                           # number of rows (individuals)
  snps::Int                             # number of columns (snps)
  personid::Vector{AbstractString}      # names of individuals
  snpid::Vector{AbstractString}         # SNP ids
  chromosome::Vector{AbstractString}    # SNP chromosome
  genetic_distance::Vector{Float64}     # genetic distance
  basepairs::Vector{Int}                # SNP base pair position
  allele1::Vector{AbstractString}       # A1 code
  allele2::Vector{AbstractString}       # A2 code
  maf::Vector{Float64}                  # minor allele frequencies
  minor_allele::BitVector               # bitvector designating the minor allele
  snpmatrix::AbstractSnpArray           # matrix of genotypes or haplotypes
  missings_per_person::Vector{Float64}  # number of missing genotypes per person
  missings_per_snp::Vector{Float64}     # number of missing genotypes per snp
end # end SnpDataStruct

max_stars_repo_name: OpenMendel/MendelBase

max_issues_repo_name: OpenMendel/MendelBase

max_forks_repo_name: OpenMendel/MendelBase


text:
using LinkedLists

function consistent(l::LinkedList)
    node=l.node
    while node.next!=l.node
        @assert(node.next.prev==node)
        @assert(node.prev.next==node)
        node=node.next
    end
end

function walkthrough_list()
    l=LinkedList{Int}()
    @assert(isempty(l))
    @assert(length(l)==0)
    vals=[3,2,7,7,9]

    prepend!(l, vals)
    consistent(l)

    println(l)
    println(vals)
    compare(vals, l)
    compare_iterator(vals, l)

    @assert(!isempty(l))
    show(l)

    pushfirst!(vals, 8)
    pushfirst!(l, 8)
    compare(vals, l)
    consistent(l)

    @assert(length(l)==6)

    @assert(indextoposition(findall(x -> x∈[2,3], l) ,l) == findall(x -> x∈[2,3], vals))
    @assert(indextoposition(indexin([2,3,12], l),l) == indexin([2,3,12], vals))

    @assert(unique(l)==unique(vals))

    @assert(getindex(l, lastindex(l))==9)

    @assert(!in(13, l))
    @assert(in(2, l))
    @assert(eltype(l)==Int)
    @assert(first(l)==8)
    @assert(last(l)==9)
    show(l)
    println()
    push!(l, 11)
    consistent(l)
    @assert(pop!(l)==11)
    compare(vals, l)
    consistent(l)
    @assert(popfirst!(l)==8)
    popfirst!(vals)
    compare(vals, l)
    consistent(l)

    two=findfirst(x -> x==2 ,l)
    insert!(l, two, 6)
    insert!(vals, 2, 6)
    consistent(l)
    println(l)
    println(vals)
    compare(vals, l)
    deleteat!(l, two)
    deleteat!(vals, 3)
    compare(vals, l)
    splice!(l, findfirst(x -> x==9, l))
    deleteat!(vals, 5)
    compare(vals, l)
    splice!(l, findfirst(x -> x==6, l), 12)
    vals[2]=12
    compare(vals, l)
    append!(l, [23, 24])
    append!(vals, [23, 24])
    compare(vals, l)
    prepend!(l, [17, 16])
    prepend!(vals, [17, 16])
    println(l)
    println(vals)
    compare(vals, l)
end

walkthrough_list()


text:
module HyperModularity

import Combinatorics #
import Distributions #
import LinearAlgebra #
import Parameters #
import Random #
import SparseArrays #
import StatsBase #
import Optim #
import SpecialFunctions #

include("omega.jl")     # uses Combinatorics
include("HSBM.jl")      # uses Parameters, Distributions
include("dict_ops.jl")
include("utils.jl")
include("objectives.jl")
include("inference.jl") # uses Optim
include("analysis_helpers.jl")
include("vol.jl")       # uses StatsBase
include("louvain_utils.jl")
include("graph_louvain.jl") # uses Special Functions
include("read_data.jl")
include("AON_hyperlouvain.jl") # uses Random
include("general_hypergraph_louvain.jl")
include("simple_synthetic.jl")

export partitionize

# from analysis_helpers.jl
export downSampleEdges!
export subhypergraph
export removeEdges!
export kcore
export projectedGraph
export mutualInformation

# from HSBM.jl
export sampleSBM
export hypergraph

# from omega.jl
export partitionsUpTo
export partitionIntensityFunction
export IntensityFunction
export allOrNothingIntensityFunction
export sumOfExteriorDegreesIntensityFunction
export empiricalIntensityFunction

# from objectives.jl
export first_term_eval
export modularity
export logLikelihood

# from inference.jl
export estimateΩEmpirically
export learnParameters
export formObjective
export formObjectives

# from vol.jl
export evalSums
export evalConstants
export increments
export addIncrements
export second_term_eval
export momentIncrements

# from read_data.jl
export read_hypergraph_label_names
export read_hypergraph_labels
export read_hypergraph_edges
export read_hypergraph_data
export hypermodularity_datasets

# from AON_hypergraph_louvain.jl; uses Random
export Simple_AON_Louvain
export AON_Louvain
export learn_omega_aon
export modularity_aon

# from general_hypegraph_louvain.jl
export HyperLouvain
export SuperNodeLouvain
export evalCuts

# from graph_louvain.jl
export CliqueExpansion
export CliqueExpansionModularity
export StarExpansionModularity
export computeDyadicResolutionParameter
export dyadicModularity
export dyadicLogLikelihood
export dyadicLogLikelihoodDirect

# from simple_synthetic.jl
export SimpleSyntheticHypergraph

# from louvain_utils
export incidence2elist
export elist2incidence
export NeighborList

# from analysis_helpers
export mutualInformation

end

max_stars_repo_name: nveldt/TestPkg

max_issues_repo_name: nveldt/TestPkg

max_forks_repo_name: nveldt/TestPkg


text:
__precompile__()
module ASTInterpreter2

using DebuggerFramework
using DebuggerFramework: FileLocInfo, BufferLocInfo, Suppressed
using Base.Meta
using REPL.LineEdit
using REPL
import Base: +, deepcopy_internal
using Core: CodeInfo, SSAValue, SlotNumber, TypeMapEntry, SimpleVector, LineInfoNode, GotoNode, Slot, GeneratedFunctionStub
using Markdown

export @enter, @make_stack

include("interpret.jl")

struct JuliaProgramCounter
    next_stmt::Int
end
+(x::JuliaProgramCounter, y::Integer) = JuliaProgramCounter(x.next_stmt+y)

struct JuliaStackFrame
    meth::Method
    code::CodeInfo
    locals::Vector{Any}
    ssavalues::Vector{Any}
    used::BitSet
    sparams::Vector{Any}
    exception_frames::Vector{Int}
    last_exception::Ref{Any}
    pc::JuliaProgramCounter
    # A vector from names to the slotnumber of that name
    # for which a reference was last encountered.
    last_reference::Dict{Symbol, Int}
    wrapper::Bool
    generator::Bool
    # Display options
    fullpath::Bool
end
function JuliaStackFrame(frame::JuliaStackFrame, pc::JuliaProgramCounter; wrapper = frame.wrapper, generator=frame.generator, fullpath=frame.fullpath)
    JuliaStackFrame(frame.meth, frame.code, frame.locals,
                    frame.ssavalues, frame.used, frame.sparams,
                    frame.exception_frames, frame.last_exception,
                    pc, frame.last_reference, wrapper, generator,
                    fullpath)
end

is_loc_meta(expr, kind) = isexpr(expr, :meta) && length(expr.args) >= 1 && expr.args[1] === kind

function DebuggerFramework.locdesc(frame::JuliaStackFrame, specslottypes = false)
    sprint() do io
        argnames = frame.code.slotnames[2:frame.meth.nargs]
        spectypes = Any[Any for i=1:length(argnames)]
        print(io, frame.meth.name,'(')
        first = true
        for (argname, argT) in zip(argnames, spectypes)
            first || print(io, ", ")
            first = false
            print(io, argname)
            !(argT === Any) && print(io, "::", argT)
        end
        print(io, ") at ",
            frame.fullpath ? frame.meth.file :
            basename(String(frame.meth.file)),
            ":",frame.meth.line)
    end
end

function DebuggerFramework.print_locals(io::IO, frame::JuliaStackFrame)
    for i = 1:length(frame.locals)
        if !isa(frame.locals[i], Nothing)
            # #self# is only interesting if it has values inside of it. We already know
            # which function we're in otherwise.
            val = something(frame.locals[i])
            if frame.code.slotnames[i] == Symbol("#self#") && (isa(val, Type) || sizeof(val) == 0)
                continue
            end
            DebuggerFramework.print_var(io, frame.code.slotnames[i], something(frame.locals[i]), nothing)
        end
    end
    for i = 1:length(frame.sparams)
        DebuggerFramework.print_var(io, frame.meth.sparam_syms[i], frame.sparams[i], nothing)
    end
end

const SEARCH_PATH = []
__init__() = append!(SEARCH_PATH,[joinpath(Sys.BINDIR,"../share/julia/base/"),
    joinpath(Sys.BINDIR,"../include/")])
function loc_for_fname(file, line, defline)
    if startswith(string(file),"REPL[")
        hist_idx = parse(Int,string(file)[6:end-1])
        isdefined(Base, :active_repl) || return nothing, ""
        hp = Base.active_repl.interface.modes[1].hist
        return BufferLocInfo(hp.history[hp.start_idx+hist_idx], line, 0, defline)
    else
        for path in SEARCH_PATH
            fullpath = joinpath(path,string(file))
            if isfile(fullpath)
                return FileLocInfo(fullpath, line, 0, defline)
            end
        end
    end
    return nothing
end

function DebuggerFramework.locinfo(frame::JuliaStackFrame)
    loc_for_fname(frame.meth.file, location(frame), frame.meth.line)
end

function lookup_var_if_var(frame, x)
    if isa(x, Union{SSAValue, GlobalRef, SlotNumber}) || isexpr(x, :static_parameter) || isexpr(x, :the_exception) || isexpr(x, :boundscheck)
        return lookup_var(frame, x)
    end
    x
end

function DebuggerFramework.eval_code(state, frame::JuliaStackFrame, command)
    expr = Base.parse_input_line(command)
    local_vars = Any[]
    local_vals = Any[]
    for i = 1:length(frame.locals)
        if !isa(frame.locals[i], Nothing)
            push!(local_vars, frame.code.slotnames[i])
            push!(local_vals, QuoteNode(something(frame.locals[i])))
        end
    end
    for i = 1:length(frame.sparams)
        push!(local_vars, frame.meth.sparam_syms[i])
        push!(local_vals, QuoteNode(frame.sparams[i]))
    end
    res = gensym()
    eval_expr = Expr(:let,
        Expr(:block, map(x->Expr(:(=), x...), zip(local_vars, local_vals))...),
        Expr(:block,
            Expr(:(=), res, expr),
            Expr(:tuple, res, Expr(:tuple, local_vars...))
        ))
    eval_res, res = Core.eval(frame.meth.module, eval_expr)
    j = 1
    for i = 1:length(frame.locals)
        if !isa(frame.locals[i], Nothing)
            frame.locals[i] = Some(res[j])
            j += 1
        end
    end
    for i = 1:length(frame.sparams)
        frame.sparams[i] = res[j]
        j += 1
    end
    eval_res
end

function maybe_quote(x)
    (isa(x, Expr) || isa(x, Symbol)) ? QuoteNode(x) : x
end

function DebuggerFramework.print_next_state(io::IO, state, frame::JuliaStackFrame)
    print(io, "About to run: ")
    expr = pc_expr(frame, frame.pc)
    isa(expr, Expr) && (expr = copy(expr))
    if isexpr(expr, :(=))
        expr = expr.args[2]
    end
    if isexpr(expr, :call) || isexpr(expr, :return)
        expr.args = map(var->maybe_quote(lookup_var_if_var(frame, var)), expr.args)
    end
    if isa(expr, Expr)
        for (i, arg) in enumerate(expr.args)
            try
                nbytes = length(repr(arg))
                if nbytes > max(40, div(200, length(expr.args)))
                    expr.args[i] = Suppressed("$nbytes bytes of output")
                end
            catch
                expr.args[i] = Suppressed("printing error")
            end
        end
    end
    print(io, expr)
    println(io)
end

const all_commands = ("q", "s", "si", "finish", "bt", "loc", "ind", "shadow",
    "up", "down", "nc", "n", "se")

function DebuggerFramework.language_specific_prompt(state, frame::JuliaStackFrame)
    if haskey(state.language_modes, :julia)
        return state.language_modes[:julia]
    end
    julia_prompt = LineEdit.Prompt(DebuggerFramework.promptname(state.level, "julia");
        # Copy colors from the prompt object
        prompt_prefix = state.repl.prompt_color,
        prompt_suffix = (state.repl.envcolors ? Base.input_color : state.repl.input_color),
        complete = REPL.REPLCompletionProvider(),
        on_enter = REPL.return_callback)
    # 0.7 compat
    if isdefined(state.main_mode, :repl)
        julia_prompt.repl = state.main_mode.repl
    end
    julia_prompt.hist = state.main_mode.hist
    julia_prompt.hist.mode_mapping[:julia] = julia_prompt

    julia_prompt.on_done = (s,buf,ok)->begin
        if !ok
            LineEdit.transition(s, :abort)
            return false
        end
        xbuf = copy(buf)
        command = String(take!(buf))
        ok, result = DebuggerFramework.eval_code(state, command)
        REPL.print_response(state.repl, ok ? result : result[1], ok ? nothing : result[2], true, true)
        println(state.repl.t)

        if !ok
            # Convenience hack. We'll see if this is more useful or annoying
            for c in all_commands
                !startswith(command, c) && continue
                LineEdit.transition(s, state.main_mode)
                LineEdit.state(s, state.main_mode).input_buffer = xbuf
                break
            end
        end
        LineEdit.reset_state(s)
    end
    julia_prompt.keymap_dict = LineEdit.keymap([REPL.mode_keymap(state.main_mode);state.standard_keymap])
    state.language_modes[:julia] = julia_prompt
    return julia_prompt
end

function JuliaStackFrame(meth::Method)
    JuliaStackFrame(meth, Vector{Any}(),
        Vector{Any}(), Vector{Any}(), Vector{Any}(),
        Dict{Symbol, Int}(), false, false, true)
end

function DebuggerFramework.debug(meth::Method, args...)
    stack = [JuliaStackFrame(meth)]
    DebuggerFramework.RunDebugger(stack)
end

function to_function(x)
    if isa(x, Function) || isa(x, Core.IntrinsicFunction)
        x
    elseif isa(x, GlobalRef)
        eval(x)
    else
        x
    end
end

_Typeof(x) = isa(x,Type) ? Type{x} : typeof(x)

function is_function_def(ex)
    (isa(ex,Expr) && ex.head == :(=) && isexpr(ex.args[1],:call)) ||
    isexpr(ex,:function)
end

function is_generated(meth)
    isdefined(meth, :generator)
end

function determine_method_for_expr(expr; enter_generated = false)
    f = to_function(expr.args[1])
    allargs = expr.args
    # Extract keyword args
    local kwargs = Expr(:parameters)
    if length(allargs) > 1 && isexpr(allargs[2], :parameters)
        kwargs = splice!(allargs, 2)
    end
    if !isempty(kwargs.args)
        of = f
        f = Core.kwfunc(f)
        allargs = [f,reduce(vcat,Any[[ex.args[1];ex.args[2]] for ex in kwargs.args]),of,
            allargs[2:end]...]
    elseif f === Core._apply
        f = to_function(allargs[2])
        allargs = Base.append_any((allargs[2],), allargs[3:end]...)
    end
    # Can happen for thunks created by generated functions
    if !isa(f, Core.Builtin) && !isa(f, Core.IntrinsicFunction)
        args = allargs[2:end]
        argtypes = Tuple{map(_Typeof,args)...}
        method = try
            which(f, argtypes)
        catch err
            @show typeof(f)
            println(f)
            println(argtypes)
            rethrow(err)
        end
        argtypes = Tuple{_Typeof(f), argtypes.parameters...}
        args = allargs
        sig = method.sig
        isa(method, TypeMapEntry) && (method = method.func)
        # Get static parameters
        (ti, lenv) = ccall(:jl_type_intersection_with_env, Any, (Any, Any),
                            argtypes, sig)::SimpleVector
        if is_generated(method) && !enter_generated
            # If we're stepping into a staged function, we need to use
            # the specialization, rather than stepping thorugh the
            # unspecialized method.
            code = Core.Compiler.get_staged(Core.Compiler.code_for_method(method, argtypes, lenv, typemax(UInt), false))
        else
            if is_generated(method)
                args = map(_Typeof, args)
                code = get_source(method.generator)
            else
                code = get_source(method)
            end
        end
        return code, method, args, lenv
    end
    nothing
end

function get_source(meth)
    if isa(meth.source, Array{UInt8,1})
        return ccall(:jl_uncompress_ast, Any, (Any, Any), meth, meth.source)
    else
        return meth.source
    end
end

function get_source(g::GeneratedFunctionStub)
    b = g(g.argnames...)
    b isa CodeInfo && return b
    return eval(b)
end

function Base.deepcopy_internal(x::LineInfoNode, stackdict::IdDict)
    if haskey(stackdict, x)
        return stackdict[x]
    end
    deeper(x) = deepcopy_internal(x, stackdict)
    stackdict[x] = LineInfoNode(x.mod, deeper(x.method),
        deeper(x.file), deeper(x.line), deeper(x.inlined_at))
end

function copy_codeinfo(code::CodeInfo)
    old_code = code.code
    code.code = UInt8[]
    new_codeinfo = deepcopy(code)
    new_codeinfo.code = old_code
    code.code = old_code
    new_codeinfo
end

function prepare_locals(meth, code, argvals = (), generator = false)
    code = copy_codeinfo(code)
    # Construct the environment from the arguments
    argnames = code.slotnames[1:meth.nargs]
    locals = Array{Any}(undef, length(code.slotflags))
    ng = isa(code.ssavaluetypes, Int) ? code.ssavaluetypes : length(code.ssavaluetypes)
    ssavalues = Array{Any}(undef, ng)
    sparams = Array{Any}(undef, length(meth.sparam_syms))
    for i = 1:meth.nargs
        if meth.isva && i == length(argnames)
            locals[i] = length(argvals) >= i ? Some(tuple(argvals[i:end]...)) : Some(())
            break
        end
        locals[i] = length(argvals) >= i ? Some(argvals[i]) : Some(())
    end
    # add local variables initially undefined
    for i = (meth.nargs+1):length(code.slotnames)
        locals[i] = nothing
    end
    used = find_used(code)
    JuliaStackFrame(meth, code, locals, ssavalues, used, sparams,  Int[], nothing,
        JuliaProgramCounter(1), Dict{Symbol,Int}(), false, generator,
        true)
end

function enter_call_expr(expr; enter_generated = false)
    r = determine_method_for_expr(expr; enter_generated = enter_generated)
    if r !== nothing
        code, method, args, lenv = r
        frame = prepare_locals(method, code, args, enter_generated)
        # Add static parameters to environment
        for i = 1:length(lenv)
            frame.sparams[i] = lenv[i]
        end
        return frame
    end
    nothing
end

function maybe_step_through_wrapper!(stack)
    last = stack[1].code.code[end-1]
    isexpr(last, :(=)) && (last = last.args[2])
    is_kw = startswith(String(Base.unwrap_unionall(stack[1].meth.sig).parameters[1].name.name), "#kw")
    if is_kw || isexpr(last, :call) && any(x->x==SlotNumber(1), last.args)
        # If the last expr calls #self# or passes it to an implemetnation method,
        # this is a wrapper function that we might want to step though
        frame = stack[1]
        pc = frame.pc
        while pc != JuliaProgramCounter(length(frame.code.code)-1)
            pc = next_call!(frame, pc)
        end
        stack[1] = JuliaStackFrame(frame, pc; wrapper=true)
        pushfirst!(stack, enter_call_expr(Expr(:call, map(x->lookup_var_if_var(frame, x), last.args)...)))
        return maybe_step_through_wrapper!(stack)
    end
    stack
end

lower(mod, arg) = false ? expand(arg) : Meta.lower(mod, arg)

# This is a version of gen_call_with_extracted_types, except that is passes back the call expression
# for further processing.
function extract_args(__module__, ex0)
    if isa(ex0, Expr)
        kws = collect(filter(x->isexpr(x,:kw),ex0.args))
        if !isempty(kws)
            names = []
            values = Tuple(map(x-> begin
                push!(names,x.args[1])
                x.args[2]
            end,kws))
            names = Tuple(names)
            return Expr(:tuple,:(Core.kwfunc($(ex0.args[1]))),
                Expr(:call, NamedTuple{names,typeof(values)}, values),
                map(x->isexpr(x, :parameters) ? QuoteNode(x) : x,
                filter(x->!isexpr(x, :kw),ex0.args))...)
        else
            return Expr(:tuple,
                map(x->isexpr(x,:parameters) ? QuoteNode(x) : x, ex0.args)...)
        end
    end
    if isa(ex0, Expr) && ex0.head == :macrocall # Make @edit @time 1+2 edit the macro by using the types of the *expressions*
        return error("Macros are not supported in @enter")
    end
    ex = lower(__module__, ex0)
    exret = Expr(:none)
    if !isa(ex, Expr)
        return error("expression is not a function call or symbol")
    elseif ex.head == :call
        return Expr(:tuple,
                map(x->isexpr(x,:parameters) ? QuoteNode(x) : x, ex.args)...)
    elseif ex.head == :body
        a1 = ex.args[1]
        if isa(a1, Expr) && a1.head == :call
            a11 = a1.args[1]
            if a11 == :setindex!
                return Expr(:tuple,
                map(x->isexpr(x,:parameters) ? QuoteNode(x) : x, arg.args)...)
            end
        end
    end
    return error("expression is not a function call, "
               * "or is too complex for @enter to analyze; "
               * "break it down to simpler parts if possible")
end

function _make_stack(mod, arg)
    args = try
        extract_args(mod, arg)
    catch e
        return :(throw($e))
    end
    quote
        theargs = $(esc(args))
        stack = [ASTInterpreter2.enter_call_expr(Expr(:call,theargs...))]
        ASTInterpreter2.maybe_step_through_wrapper!(stack)
        stack[1] = ASTInterpreter2.JuliaStackFrame(stack[1], ASTInterpreter2.maybe_next_call!(stack[1]))
        stack
    end
end

macro make_stack(arg)
    _make_stack(__module__, arg)
end

macro enter(arg)
    quote
        let stack = $(_make_stack(__module__,arg))
            DebuggerFramework.RunDebugger(stack)
        end
    end
end

include("commands.jl")

end # module


text:
using MCPhylo, Random, Test

Random.seed!(123)

function simulate_NDARMA(N::Int, p::Int, q::Int, prob::Vector{Float64},
                         phi::Vector{Float64})
  X = zeros(Int64, N)
  X[1:p] = rand(Categorical(prob), p)
  d1 = Multinomial(1, phi)
  d2 = Categorical(prob)
  for t in (p+1):N
    alphabeta = rand(d1)
    eps = rand(d2, q + 1)
    X[t] = sum([X[(t-p):(t-1)]; eps] .* alphabeta)
  end
  return X
end

n = 10000
x1 = simulate_NDARMA(n, 1, 0, [0.25, 0.5, 0.05, 0.2], [0.8, 0.2])
x2 = simulate_NDARMA(n, 1, 0, [0.25, 0.5, 0.05, 0.2], [0.8, 0.2])
x3 = simulate_NDARMA(n, 1, 0, [0.25, 0.5, 0.05, 0.2], [0.8, 0.2])

sim = Chains(cat(x1,x2,x3,dims=3))
discretediag(sim)

res = MCPhylo.weiss(hcat(x1,x2,x3))
@test res[3] > 0.05

x1 = simulate_NDARMA(n, 1, 0, [0.25, 0.5, 0.05, 0.2], [0.8, 0.2])
x2 = simulate_NDARMA(n, 1, 0, [0.75, 0.15, 0.05, 0.05], [0.8, 0.2])
x3 = simulate_NDARMA(n, 1, 0, [0.25, 0.5, 0.05, 0.2], [0.8, 0.2])
res = MCPhylo.weiss(hcat(x1,x2,x3))
@test res[3] < 0.05

#mark that we got to the end of the test file succesfully
@test true


text:
This is a program.

    Now I fixed a bug.

    Now I fixed a second bug.

    Third bug fixed.

    This is dy/bugfix change.

max_stars_repo_name: denizyuret/gittut

max_issues_repo_name: denizyuret/gittut

max_forks_repo_name: denizyuret/gittut


text:
function δ(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E}
    for s in species(N)
        @assert haskey(L, s)
    end

    this_l = filter(p -> p.first in species(N), L)
    tl = [this_l[s] for s in species(N; dims=1)]
    bl = [this_l[s] for s in species(N; dims=2)]
    D = tl .== bl'
    return D
end

"""
    Q(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E}

Modularity of a network and its partition. The second argument is a dictionary
where every species of `N` is associated to an `Int64` value representing the
identity of the module. This function returns the same value of bipartite
networks and their unipartite projection.

#### References

Barber, M.J., 2007. Modularity and community detection in bipartite networks.
Phys. Rev. E 76, 066102. https://doi.org/10.1103/PhysRevE.76.066102

Newman, M.E., 2006. Modularity and community structure in networks. Proceedings
of the national academy of sciences 103, 8577–8582.
"""
function Q(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E}
    for s in species(N)
        @assert haskey(L, s)
    end

    # Degrees
    dkin, dkout = degree_in(N), degree_out(N)
    kin = map(x -> dkin[x], species(N; dims=2))
    kout = map(x -> dkout[x], species(N; dims=1))

    # Value of m -- sum of weights, total number of int, ...
    m = links(N)

    # Null model
    kikj = (kout .* kin')
    Pij = kikj ./ m

    # Difference
    diff = N.A .- Pij

    # Diff × delta
    dd = diff .* δ(N,L)

    return sum(dd)/m
end

"""
    Qr(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E}

Realized modularity -- this function returns a value giving the proportion of
all links that are within the same module. Higher values reflect a more strongly
modular partition (whereas `Q` represents the deviation of modularity from the
random expectation).

#### References

Poisot, T., 2013. An a posteriori measure of network modularity. F1000Research
2. https://doi.org/10.12688/f1000research.2-130.v3
"""
function Qr(N::T, L::Dict{E,Int64}) where {T<:AbstractEcologicalNetwork,E}
    for s in species(N)
        @assert haskey(L, s)
    end
    W = sum(N.A .* δ(N, L))
    B = links(N)
    return 2.0 * (W/B) - 1.0
end

function tidy_modules!(L::Dict{E,Int64}) where {E}
  l_values = sort(unique(collect(values(L))))
  for i in 1:length(l_values)
    for (k,v) in L
      if v == l_values[i]
        L[k] = i
      end
    end
  end
end


text:
@testset "LineSegment-LineSegment" begin
    for T ∈ Floats
        # basic intersection
        l₁ = LineSegment(Point{2,T}(0,  1), Point{2,T}(2, -1))
        l₂ = LineSegment(Point{2,T}(0, -1), Point{2,T}(2,  1))
        p = intersect(l₁, l₂)
        @test p ≈ Point{2,T}(1, 0)
        @test typeof(p) == Point{2,T}

        # vertex intersection
        l₂ = LineSegment(Point{2,T}(0, -1), Point{2,T}(2, -1))
        p = l₁ ∩ l₂
        @test p ≈ Point{2,T}(2, -1)

        # vertical
        l₁ = LineSegment(Point{2,T}(0,  1), Point{2,T}(2,   1))
        l₂ = LineSegment(Point{2,T}(1, 10), Point{2,T}(1, -10))
        p = intersect(l₁, l₂)
        @test p ≈ Point{2,T}(1, 1)

        # nearly vertical
        l₁ = LineSegment(Point{2,T}(-1, -100000), Point{2,T}(1,  100000))
        l₂ = LineSegment(Point{2,T}(-1,   10000), Point{2,T}(1,  -10000))
        p = l₁ ∩ l₂
        @test p ≈ Point{2,T}(0, 0)

        # parallel
        l₁ = LineSegment(Point{2,T}(0, 1), Point{2,T}(1, 1))
        l₂ = LineSegment(Point{2,T}(0, 0), Point{2,T}(1, 0))
        p = intersect(l₁, l₂)
        @test p ≈ Point{2,T}(1e6,1e6)

        # collinear
        l₁ = LineSegment(Point{2,T}(0, 0), Point{2,T}(2, 0))
        l₂ = LineSegment(Point{2,T}(0, 0), Point{2,T}(1, 0))
        p = intersect(l₁, l₂)
        @test p ≈ Point{2,T}(1e6,1e6)

        # line intersects, not segment (invalid s)
        l₁ = LineSegment(Point{2,T}(0, 0), Point{2,T}(2, 0    ))
        l₂ = LineSegment(Point{2,T}(1, 2), Point{2,T}(1, 1//10))
        p = l₁ ∩ l₂
        @test p ≈ Point{2,T}(1e6,1e6)

        # line intersects, not segment (invalid r)
        l₂ = LineSegment(Point{2,T}(0, 0), Point{2,T}(2, 0    ))
        l₁ = LineSegment(Point{2,T}(1, 2), Point{2,T}(1, 1//10))
        p = intersect(l₁, l₂)
        @test p ≈ Point{2,T}(1e6,1e6)
    end
end

max_stars_repo_name: khurrumsaleem/MOCNeutronTransport

max_issues_repo_name: khurrumsaleem/MOCNeutronTransport

max_forks_repo_name: khurrumsaleem/MOCNeutronTransport


text:
using InvariantMeasures: root, nthpreimage!
using StaticArrays
using ValidatedNumerics

@testset "Contractors" begin

@test root(x -> x^2-2, 1..2, 1e-13) ≈ sqrt(2)
# to ensure quadratic convergence
@test root(x -> x^2-2, 1..2, 1e-13; max_iter = 6) ≈ sqrt(2)

fs = (x -> x^2, x -> x^3)
X = [0..1, 0..1]
y = 0.1
nthpreimage!(y, fs, X; max_iter = 100)
@test X[1]^6 ≈ y
@test X[2]^3 ≈ y
X = @MVector[0..1, 0..1]
nthpreimage!(y, fs, X; max_iter = 9)
@test X[1]^6 ≈ y
@test X[2]^3 ≈ y

fs = (x -> x^2, x -> x^3, x -> x^4)
X = @MVector[0..1, 0..1, 0..1]
nthpreimage!(y, fs, X)
@test X[1]^24 ≈ y

end #testset


text:
Base.@kwdef mutable struct FileDialogModel <: AbstractDialogModel
    confirmed_path::Path = Path(pwd(), "")
    unconfirmed_path::Path = Path(pwd(), "")
end

mutable struct FileDialogControl <: AbstractDialogControl
    isenabled::Bool
end

function get_path(model::FileDialogModel, status::ConfirmedStatus)
    model.confirmed_path
end

function get_path(model::FileDialogModel, status::UnconfirmedStatus)
    model.unconfirmed_path
end

function set_path!(model::FileDialogModel, path::Path, status::ConfirmedStatus)
    model.confirmed_path = path
end

function set_path!(model::FileDialogModel, path::Path, status::UnconfirmedStatus)
    model.unconfirmed_path = path
end

function get_caption(property::FileDialogDisplayProperties)
    property.caption
end

function get_action(property::FileDialogDisplayProperties)
    property.action
end

function get_cursor_position(property::FileDialogDisplayProperties)
    property.position
end

function get_width(property::FileDialogDisplayProperties)
    property.get_width
end

function get_height(property::FileDialogDisplayProperties)
    property.height
end

function get_directory(p::Path)
    p.directory
end

function get_filename(p::Path)
    p.filename
end

function (control::FileDialogControl)(model::FileDialogModel, properties::AbstractDisplayProperties, operation::AbstractOperation)
    @c CImGui.Begin(get_caption(properties), &control.isenabled)
        facilitate_path_selection!(model)
        facilitate_directory_file_selection!(model)
        handle_cancellation_confirmation!(control, model,  get_action(properties), operation)
        handle_file_error_messaging()
    CImGui.End()
end

function facilitate_path_selection!(model::FileDialogModel)
    path = get_path(model,  UnconfirmedStatus())
    path_directories = splitpath(get_directory(path))
    selected_directory = Cint(length(path_directories))
    # Draw a button for each directory that constitutes the current path.
    for (index, d) in enumerate(path_directories)
        CImGui.Button(d) && (selected_directory = Cint(index);)
        CImGui.SameLine()
    end
    # If a button is clicked then we keep only the directories up-to and including the clicked button.
    truncated_directories = selected_directory == 1 ? joinpath(first(path_directories)) : joinpath(path_directories[1:selected_directory]...)
    path₂ = Path(truncated_directories, get_filename(path))
    set_path!(model, path₂, UnconfirmedStatus())
end

function facilitate_directory_file_selection!(model::FileDialogModel)
    # Make a list of files and directories that are visibile from the current directory.
    CImGui.NewLine()
    CImGui.BeginChild("Directory and File Listing", CImGui.ImVec2(CImGui.GetWindowWidth() * 0.98, -CImGui.GetWindowHeight() * 0.2))
        CImGui.Columns(1)
        handle_directory_selection!(model)
        handle_file_selection!(model)
    CImGui.EndChild()
end

function handle_directory_selection!(model::FileDialogModel)
    path = get_path(model, UnconfirmedStatus())
    current_directory = get_directory(path)
    visible_directories = filter(p->is_readable_dir(joinpath(current_directory, p)), readdir(current_directory))
    for (n, foldername) in enumerate(visible_directories)
        # When the user clicks on a directory then change directory by appending the selected directory to the current path.
        if CImGui.Selectable("[Dir] " * "$foldername")
            selected_directory = joinpath(current_directory, foldername)
            path₂ = Path(selected_directory, "")
            set_path!(model, path₂, UnconfirmedStatus())
        end
    end
end

function handle_file_selection!(model::FileDialogModel)
    path = get_path(model, UnconfirmedStatus())
    current_directory = get_directory(path)
    visible_files = filter(p->is_queryable_file(joinpath(current_directory, p)), readdir(current_directory))
    selected_file = Cint(0)
    for (n, filename) in enumerate(visible_files)
        if CImGui.Selectable("[File] " * "$filename")
            path₂ = Path(current_directory, filename)
            set_path!(model, path₂, UnconfirmedStatus())
        end
    end
    handle_keyboard!(model)
end

function handle_keyboard!(model::FileDialogModel)
    CImGui.Text("File Name:")
    CImGui.SameLine()
    path = get_path(model, UnconfirmedStatus())
    current_directory = get_directory(path)
    filename₀ = get_filename(path)
    filename₁ = filename₀*"\0"^(1)
    # Allow up to 255 characters for the filename.
    pad_null = max(0, 255 - length(filename₁) + 1)
    buffer = filename₁*"\0"^(pad_null)
    CImGui.InputText("",buffer, length(buffer))
    filename₂ = extract_string(buffer)
    path₂ = Path(current_directory, filename₂)
    set_path!(model, path₂, UnconfirmedStatus())
end

function handle_cancellation_confirmation!(control::FileDialogControl, model::FileDialogModel, action::String, operation::AbstractOperation)
    CImGui.Button("Cancel") && disable!(control)
    CImGui.SameLine()
    CImGui.Button(action) && (handle_confirmation!(control, model, operation))
end

function handle_confirmation!(control::FileDialogControl, model::FileDialogModel, operation::AbstractOperation)
    path = get_path(model, UnconfirmedStatus())
    directory = get_directory(path)
    filename = get_filename(path)
    path_to_file = joinpath(directory, filename)
    if is_queryable_file(path_to_file)
        path₂ = Path(directory, filename)
        set_path!(model, path₂, ConfirmedStatus())
        enable!(operation)
        disable!(control)
    else
        CImGui.OpenPopup("Does the file exist?")
    end
end

function handle_confirmation!(control::FileDialogControl, model::FileDialogModel, operation::AbstractExporter)
    path = get_path(model, UnconfirmedStatus())
    directory = get_directory(path)
    filename = get_filename(path)
    path_to_file = joinpath(directory, filename)
    path₂ = Path(directory, filename)
    set_path!(model, path₂, ConfirmedStatus())
    enable!(operation)
    disable!(control)
end

function handle_file_error_messaging()
    if CImGui.BeginPopupModal("Does the file exist?", C_NULL, CImGui.ImGuiWindowFlags_AlwaysAutoResize)
        CImGui.Text("Unable to access the specified file.\nPlease verify that: \n   (1) the file exists; \n   (2) you have permission to access the file.\n\n")
        CImGui.Separator()
        CImGui.Button("OK", (120, 0)) && CImGui.CloseCurrentPopup()
        CImGui.SetItemDefaultFocus()
        CImGui.EndPopup()
    end
end

# The isdir function might not have permissions to query certan folders and
# will thus throw an ERROR: "IOError: stat: permission denied (EACCES)"
function is_readable_dir(path)
    flag = false
    try
        flag = isdir(path)
    catch x
        flag = false
    end
    return flag
end

# The isfile function might not have permissions to query certan files and
# will thus throw an ERROR: "IOError: stat: permission denied (EACCES)"
function is_queryable_file(path)
    flag = false
    try
        flag = isfile(path)
    catch x
        flag = false
    end
    return flag
end

function is_readable_file(path)
    if is_queryable_file(path)
        return (uperm(path) & 0x04 > 0) ? true :  false
    else
        return false
    end
end

function is_writeable_file(path)
    if is_queryable_file(path)
        return (uperm(path) & 0x02 > 0) ?  true :  false
    else
        return false
    end
end


text:
using Gtk, Images, ImageView, ImageMagick
include("code/appGUI.jl")

max_stars_repo_name: o-Mateo-o/ImageEditorApp

max_issues_repo_name: o-Mateo-o/ImageEditorApp

max_forks_repo_name: o-Mateo-o/ImageEditorApp


text:
include("berryphase.jl")

using BerryPhase

const L = 6
const Jz = 1.0
const Jxy = 1.0

const deltas = [0.22, 0.24, 0.26, 0.28]
const S2 = 2
const M = 10

open("res.dat", "w") do io
  for delta in deltas
    bp = berryphase(altchain(S2, L, Jz, Jxy, delta), M)
    println(io, delta, " ", bp/pi)
    println(delta, " ", bp/pi)
  end
end

max_stars_repo_name: yomichi/HPhi-samples

max_issues_repo_name: yomichi/HPhi-samples

max_forks_repo_name: yomichi/HPhi-samples


text:
using Random
using Pkg

Pkg.add("StatsPlots")
Pkg.add("KernelDensity")

using StatsPlots
using KernelDensity
#=
Small model of a taxi cab getting directions of a random generator

It has a positon and a set of directions. The degree of freedoms it has are
4. It can go:

{Left,Right,Up,Down}

These we can translate to vectors:

(-1 0) (1 0) (0 -1) (0 1)

Which are generators of the group of our taxi cab together with vector addition.

=#
@enum Direction L R U D
@enum Reset Rst
const DirectionVector = Vector{Int64}
const channelBufferSize = 100
const startCompressionRatio = 100

function defaultCabDirection(setSeedChannel :: Channel{Integer}, channel::Channel{Direction}, seed :: Integer)
    generator = MersenneTwister(seed)
    while true
        if isready(setSeedChannel)
            generator = MersenneTwister(take!(setSeedChannel))
        end
        put!(channel,rand(generator,(L,R,U,D)))
    end
end

function dirToVec(dir::Direction) :: DirectionVector
    if dir == L # Left
        return [-1;0]
    elseif dir == R # Right
        return [1;0]
    elseif dir == U # Up
        return [0;1]
    elseif dir == D # Down
        return [0;-1]
    end
end

function directionsMapper(input :: Channel{Direction}, output :: Channel{DirectionVector})
    while true
        dir = take!(input)
        put!(output, dirToVec(dir))
    end
end

#=
We are abusing the group structure of the direction vectors aand the
fact that we are only interested in the position to compress
them very effeciently.
=#
function compressor(compressionRatio :: Int64, setCompressionRatio :: Channel{Int64}, input :: Channel{DirectionVector}, output :: Channel{DirectionVector})
    p = Array{DirectionVector,1}(undef,compressionRatio)
    while true

        if isready(setCompressionRatio)
            compressionRatio = take!(setCompressionRatio)
        end

        for i in 1:compressionRatio
            element = take!(input)
            p[i] = element
        end
        res = sum(p)
        put!(output,res)
    end
end

function observe!(n :: Integer, channel :: Channel{T}) :: Array{T} where T <: Any
    p = Array{T,1}(undef,n)
    for i in 1:n
        e = take!(channel)
        p[i] = e
    end
    return p
end

function observeUnzip!(n :: Integer, channel :: Channel{Vector{T}}) :: Tuple{Array{T},Array{T}} where T <: Any
    x = Array{T}(undef,n)
    y = Array{T}(undef,n)
    for i in 1:n
        e = take!(channel)
        x[i] = e[1]
        y[i] = e[2]
    end
    return (x,y)
end

function destroy!(n :: Integer, channel :: Channel)
    for i in 1:n
        take!(channel)
    end
end

function taxiCab(startPostion :: DirectionVector, reset :: Channel{Reset},  input :: Channel{DirectionVector}, output :: Channel{DirectionVector})
        currentPosition = startPostion
        while true

            if isready(reset)
                rst = take!(reset)
                currentPosition = startPostion
            end

            put!(output,currentPosition)
            dir = take!(input)
            currentPosition = currentPosition + dir
        end
end

directionsChannel = Channel{Direction}(channelBufferSize)
vectorsChannel = Channel{DirectionVector}(channelBufferSize)
compressedChannel = Channel{DirectionVector}(channelBufferSize)
taxiCabPositionChannel = Channel{DirectionVector}(channelBufferSize)
resetTaxiChannel = Channel{Reset}(1)
setCompressionRatioChannel = Channel{Int64}(1)
setSeedChannel = Channel{Integer}(1)
@async defaultCabDirection(setSeedChannel, directionsChannel,1234)
@async directionsMapper(directionsChannel,vectorsChannel)
@async compressor(startCompresssionRatio,setCompressionRatioChannel, vectorsChannel, compressedChannel)
@async taxiCab([0;0],resetTaxiChannel, compressedChannel,taxiCabPositionChannel)

x,y = observeUnzip!(1000, taxiCabPositionChannel)

function resetChannel(channel :: Channel)
    for i in 1:(channelBufferSize)
        take!(channel)
    end
end

# Simple interface to the taxi
function resetWalk(compressionRatio :: Int64, seed :: Integer)
    resetChannel(taxiCabPositionChannel)
    put!(resetTaxiChannel,Rst)
    put!(setCompressionRatioChannel,compressionRatio)
    put!(setSeedChannel, seed)
end

resetWalk(compressionRatio = 100, seed = 1)

density = kde((x,y))
histogram2d(x,y)
heatmap(density)
plot(density)
plot(x,y)
scatter(x,y)


text:
module MotionModels

using Random
using StaticArrays
using Distributions

export MotionData, MotionParams, MotionModel, step!

abstract type MotionData end
abstract type MotionParams end
abstract type MotionModel <: Function end
abstract type StochasticMotionModel <: Function end

abstract type PredictModel <: Function end
abstract type StochasticPredictModel <: Function end

#(m::M)(state_t, ctrl_t) where M<:MotionModel = nothing
#(sm::SM)(state_t, ctrl_t, rng) where SM<:StochasticMotionModel = nothing

#struct BasicPredictModel{D, P} <: PredictModel
#    data::D
#    params::P
#end
#
#function predict!(state_t, ctrl_t, dt, m::BasicPredictModel)
#    m.data.pose = state_t
#    step!(m.data, m.params, dt)
#    return m.data.pose
#end

include("ackermann.jl")

struct BasicStochasticPredictModel{D, P, CD, PD, ID} <: StochasticPredictModel
    data::D
    params::P
    ctrl_dist::CD
    process_dist::PD
    init_dist::ID
end

function reset!(particles::AbstractVector{Pose2D{P}}, pose::Pose2D, rng, m::BasicStochasticPredictModel) where {P}
    statev = pose.statev
    @inbounds @simd for i in eachindex(particles)
        particles[i] = Pose2D{P}(statev + rand(rng, m.init_dist))
    end
end

function predict!(particles::AbstractVector{Pose2D{P}}, ctrl::SVector{N, C}, rng, dt, m::BasicStochasticPredictModel) where {P,N,C}
    @inbounds @simd for i in eachindex(particles)
        m.data.pose = particles[i]
        m.data.ctrl = ctrl + SVector{N, C}(rand(rng, m.ctrl_dist))
        step!(m.data, m.params, dt)
        particles[i] = Pose2D{P}(m.data.pose.statev + rand(rng, m.process_dist))
    end
end

end # module


text:
using OrderedCollections: OrderedDict

@testset "ComputeEnvironment" begin
    @testset "constructor" begin
        arn = "arn:aws:batch:us-east-1:000000000000:compute-environment/ce"
        @test ComputeEnvironment(arn).arn == arn

        patch = describe_compute_environments_patch(
            OrderedDict(
                "computeEnvironmentName" => "ce-name",
                "computeEnvironmentArn" => arn,
            )
        )
        apply(patch) do
            @test ComputeEnvironment("ce-name").arn == arn
        end

        patch = describe_compute_environments_patch()
        apply(patch) do
            @test_throws ErrorException ComputeEnvironment("ce-name")
        end
    end

    @testset "max_vcpus" begin
        ce = ComputeEnvironment("arn:aws:batch:us-east-1:000000000000:compute-environment/ce")
        patch = describe_compute_environments_patch(
            OrderedDict(
                "computeEnvironmentArn" => ce.arn,
                "computeResources" => OrderedDict("maxvCpus" => 5),
            ),
        )

        apply(patch) do
            @test AWSBatch.max_vcpus(ce) == 5
        end
    end
end


text:
#export DiscretizedNoiseKernel, UniformNoise

using FFTW, LinearAlgebra

abstract type NoiseKernel end
function Base.:*(M::NoiseKernel, v)
	@error "Not implemented"
end

opnormbound(N::NormKind, M::NoiseKernel) = @error "Not Implemented"
opradius(N::NormKind, M::NoiseKernel) = @error "Not Implemented"
nonzero_per_row(M::NoiseKernel) = @error "Not Implemented"

struct DiscretizedNoiseKernelFFT{S<:AbstractVector, T<:AbstractVector} <: NoiseKernel
	v::S
	Mfft::T
	rad
	P
end

import IntervalArithmetic: Interval, mid, radius, @interval
import InvariantMeasures: opnormbound, Linf
DiscretizedNoiseKernelFFT(v::Vector{Real}) = DiscretizedNoiseKernelFFT(v, fft(v), 0, plan_fft(mid.(v)))
DiscretizedNoiseKernelFFT(v::Vector{Interval{T}}) where {T} = DiscretizedNoiseKernelFFT(v, fft(mid.(v)), opnormbound(L2, radius.(v)), plan_fft(mid.(v)))
Mfft(Q::DiscretizedNoiseKernelFFT) = Q.Mfft

#this is the Ulam discretization of the uniform noise of size ξ, on a partition of size k
function UniformNoiseFFT(ξ, k, boundarycondition = :periodic)
	n = Int64(floor(ξ*k))
	if boundarycondition == :periodic
		v = Interval.([ones(n); zeros(k-n)])*(1/(2*ξ))
		v+= Interval.([zeros(k-n); ones(n)])*(1/(2*ξ))
		v[n+1]+= @interval (ξ-n*1/k)*k/(2*ξ)
		v[k-n-1]+= @interval (ξ-n*1/k)*k/(2*ξ)
	end
	return DiscretizedNoiseKernelFFT(v)
end

function Base.:*(M::DiscretizedNoiseKernelFFT, v)
	P = M.P
	w = P*v
	@info w
	w = M.Mfft.*w
	@info w
	return real.(P\w)/length(v)
end

struct DiscretizedNoiseKernelUlam{S<:AbstractVector} <: NoiseKernel
	ξ
	v::S
	rad
	boundarycondition::Symbol
	w::Vector
	z::Vector
end

function UniformNoiseUlam(ξ, B::Ulam, boundarycondition = :periodic)
	k = length(B)
	n = 2*Int64(ceil(ξ*k))
	v = zeros(Interval{Float64}, n)
	a = 1/(2*Interval(ξ))
	v[2:n-1] = a*ones(Interval{Float64}, n-2)
	v[1] = (k-sum(v))/2
	v[n] = v[1]
	if boundarycondition == :periodic
		return DiscretizedNoiseKernelUlam(Interval(ξ),
										  mid.(v),
										  opnormbound(L1, v)-k,
										  boundarycondition,
										  zeros(k+n),
										  zeros(k))
	elseif boundarycondition == :reflecting
		return DiscretizedNoiseKernelUlam(Interval(ξ),
										  mid.(v),
										  opnormbound(L1, v)-k,
										  boundarycondition,
										  zeros(k+n+2),
										  zeros(k))
	end
end

#TODO, but at the moment this is fine, it is a Markov operator
opnormbound(::Type{L1}, M::DiscretizedNoiseKernelUlam) = 1.0
opradius(::Type{L1}, M::DiscretizedNoiseKernelUlam) = M.rad
nonzero_per_row(M::DiscretizedNoiseKernelUlam) = length(M.v)
dfly(::Type{TotalVariation}, ::Type{L1}, N::DiscretizedNoiseKernelUlam) = return 0.0, (1/(2*N.ξ)).hi

function Base.:*(M::DiscretizedNoiseKernelUlam, v)
	mult(M, v, Val(M.boundarycondition))
end

function mult(M::DiscretizedNoiseKernelUlam, v, ::Val{:periodic})
	n = length(M.v)
	k = length(v)
	l =n÷2

	M.w .= 0

	M.w[l+1:l+k] = v
	M.w[1:l] = v[end-l+1:end]
	M.w[end-l+1:end] = v[1:l]

	for i in 1:k
		h = @view M.w[i:i+n-1]
		v[i] = sum( M.v.* h)/k
	end

	return v
end

function mult(M::DiscretizedNoiseKernelUlam, v::Vector{Interval{T}}, ::Val{:periodic}) where {T}
	n = length(M.v)
	k = length(v)
	l =n÷2

	nrmv = opnormbound(L1, v)
	midv = mid.(v)
	radv = radius.(v)
	nrmrad = opnormbound(L1, radv)

	M.w .= 0

	M.w[l+1:l+k] = midv
	M.w[1:l] = midv[end-l+1:end]
	M.w[end-l+1:end] = midv[1:l]

	for i in 1:k
		h = @view M.w[i:i+n-1]
		v[i] = sum( M.v.* h)/k
	end
	δₖ = opradius(L1, M)
	γₖ = gamma(T, nonzero_per_row(M))
    nrm_MK = opnormbound(L1, M)
    normMK = nrm_MK ⊕₊ δₖ

	ϵ = (γₖ ⊗₊ normMK) ⊗₊ nrmv ⊕₊ normMK ⊗₊ nrmrad

	return v + fill(Interval(-ϵ ,ϵ), length(v))
end

function mult(M::DiscretizedNoiseKernelUlam, v, ::Val{:reflecting})
	n = length(M.v)
	k = length(v)
	l = n÷2+1

	@info v
	M.w .= 0
	M.w[l+1:l+k] = v
	v .= 0

	for i in 1:l
		h = @view M.w[i:i+n-1]
		@info h, i
		val = sum( M.v.* h)/k
		v[l+1-i] += val
	end

	for i in 1:k
		h = @view M.w[i+l:i+l+n-1]
		@info h, i+l
		val = sum( M.v.* h)/k
		v[i]+=val
	end
	@info "stop"
	for i in 1:l
		h = @view M.w[k+l+i:k+l+i+n-1]
		@info h
		val = sum( M.v.* h)/k
		v[end-i+1] +=val
	end

	return v
end


text:
include("ManipulateTanks.jl")
include("ManipulatePumps.jl")
include("UpdateExtrema.jl")
include("MaxMinLevels.jl")
include("PipeCoefs.jl")
include("PumpCoefs.jl")
include("CheckValveCoefs.jl")
include("FlowDirections.jl")
#Pipe Coefs
function Parameterize(wn_dict::Dict{Any,Any}, n::Int64, Q_lb::Float64, logspace_ratio::Float64, dH_critical::Float64, dense_coef::Float64, tight_coef::Float64)
#Thought: should we make one q_base with pipe and valve names?
    #Pipes
    Parameters = Dict{String,Any}()
    Parameters["Qmin"], Parameters["Qmax"], Parameters["Hmin"], Parameters["Hmax"]  = Q_extrema(wn_dict, Q_lb, tight_coef)

    Parameters["Q_basePipe"] = Q_base_pipe(wn_dict, Parameters["Qmin"], Parameters["Qmax"], n, Q_lb, logspace_ratio, dH_critical, dense_coef)
    Parameters["aPipe_over"], Parameters["bPipe_over"] = PipeCoefsOver(wn_dict, Parameters["Q_basePipe"])
    Parameters["aPipe_under"], Parameters["bPipe_under"] = PipeCoefsUnder(wn_dict, Parameters["Q_basePipe"])
    # Pumps
    Parameters["flow"], Parameters["head"], Parameters["power"] = FlowHeadPower(wn_dict)
    Parameters["aPumpPower_flow"], Parameters["bPumpPower_flow"] = PumpPowerCoefs_Flow(wn_dict, Parameters["flow"], Parameters["power"])
    Parameters["aPumpPower_head"], Parameters["bPumpPower_head"] = PumpPowerCoefs_Head(wn_dict, Parameters["head"], Parameters["power"])
    #Check Valves
    Parameters["Q_baseCV"] = Q_base_checkvalve(wn_dict, Parameters["Qmax"], n, Q_lb, logspace_ratio, dH_critical, dense_coef)
    Parameters["aCheckValve_over"], Parameters["bCheckValve_over"] = CheckValveCoefsOver(Parameters["Q_baseCV"])
    Parameters["aCheckValve_under"], Parameters["bCheckValve_under"] = CheckValveCoefsUnder(Parameters["Q_baseCV"])

    #Flow Directions
    Parameters["PositiveFlowLinks"], Parameters["NegativeFlowLinks"], Parameters["ReversibleFlowLinks"] = FlowDirection(wn_dict, Parameters["Qmin"], Parameters["Qmax"])
    return Parameters
end