text:
# ============================================================================ #
# SymbolicExpression
# ============================================================================ #

mutable struct SymbolicExpression <: AbstractSymbolicExpression
    value::Float64
    head::Symbol
    args::Vector{Any}

    function SymbolicExpression(head, args...)
        return new(default_variable_value(), head, collect(args)) # [args...])
    end
end

function SymbolicExpression(x::Constant)
    return SymbolicExpression(:+, x)
end

function SymbolicExpression(x::Variable)
    return SymbolicExpression(:+, x)
end

for op in (:+, :-, :/, :*, :^)
    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Number,
            Ty <: Union{Constant, Variable, SymbolicExpression}}

        c = Constant(x)
        ex = SymbolicExpression(Symbol(Base.$op), c, y)
        return ex
    end

    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Union{Constant, Variable, SymbolicExpression},
            Ty <: Number}

        c = Constant(y)
        ex = SymbolicExpression(Symbol(Base.$op), x, c)
        return ex
    end
end

for op in (:-, :/, :^) # (:+, :-, :/, :*, :^)
    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Union{Constant, Variable, SymbolicExpression},
            Ty <: Union{Constant, Variable, SymbolicExpression}}

        ex = SymbolicExpression(Symbol(Base.$op), x, y)
        return ex
    end
end

for op in (:+, :*)
    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Union{Variable, Constant},
            Ty <: Union{Variable, Constant}}

        ex = SymbolicExpression(Symbol(Base.$op), x, y)
        return ex
    end
end

for op in (:+, :*, :-, :/, :^)
    @eval function Base.$(op)(
        x::Tx, y::Ty) where {
            Tx <: Constant,
            Ty <: Constant}

        ee = Expr(:call, $(op), x.value, y.value)
        val = eval(ee)
        ex = Constant( val )
        # ex = SymbolicExpression(Symbol(Base.$op), x, y)
        return ex
    end
end

function Base.:(+)(
    x::Tx, y::Ty) where {
        Tx <: Union{Constant, Variable},
        Ty <: SymbolicExpression}
    if is_plus_operator(y.head)
        ex = SymbolicExpression(:+, x, y.args...)
    else
        ex = SymbolicExpression(:+, x, y)
    end
    return ex
end
function Base.:(+)(
    x::Tx, y::Ty) where {
        Tx <: SymbolicExpression,
        Ty <: Union{Constant, Variable}}
    if is_plus_operator(x.head)
        ex = SymbolicExpression(:+, x.args..., y)
    else
        ex = SymbolicExpression(:+, x, y)
    end
    return ex
end
function Base.:(+)(
    x::Tx, y::Ty) where {
        Tx <: SymbolicExpression,
        Ty <: SymbolicExpression}
    if is_plus_operator(x.head) && is_plus_operator(y.head)
        ex = SymbolicExpression(:+, x.args..., y.args...)
    else
        ex = SymbolicExpression(:+, x, y)
    end
    return ex
end

function Base.:(*)(
    x::Tx, y::Ty) where {
        Tx <: Union{Constant, Variable},
        Ty <: SymbolicExpression}
    if is_product_operator(y.head)
        ex = SymbolicExpression(:*, x, y.args...)
    else
        ex = SymbolicExpression(:*, x, y)
    end
    return ex
end
function Base.:(*)(
    x::Tx, y::Ty) where {
        Tx <: SymbolicExpression,
        Ty <: Union{Constant, Variable}}
    if is_product_operator(x.head)
        ex = SymbolicExpression(:*, x.args..., y)
    else
        ex = SymbolicExpression(:*, x, y)
    end
    return ex
end
function Base.:(*)(
    x::Tx, y::Ty) where {
        Tx <: SymbolicExpression,
        Ty <: SymbolicExpression}
    if is_product_operator(x.head) && is_product_operator(y.head)
        ex = SymbolicExpression(:*, x.args..., y.args...)
    else
        ex = SymbolicExpression(:*, x, y)
    end
    return ex
end

function Base.:(+)(
    x::Tx) where {Tx <: Union{Constant, Variable, SymbolicExpression}}
    return x
end

function Base.:(-)(
    x::Tx) where {Tx <: Union{Constant, Variable, SymbolicExpression}}
    ex = SymbolicExpression(:(-), x)
    return ex
end

for op in (:sqrt, :cbrt, :exp, :log, :log2, :log10, :sin, :cos, :tan, :abs, :sign)
    @eval function Base.$(op)(
        x::Tx) where {Tx <: Union{Constant, Variable, SymbolicExpression}}

        ex = SymbolicExpression(Symbol(Base.$op), x)
        return ex
    end
end

function evaluate(ex::SymbolicExpression)

    values = Vector{Any}(undef, length(ex.args))
    for (ii, arg) in enumerate(ex.args)
        values[ii] = evaluate(arg)
    end

    ex2call  = Expr(:call, ex.head, values...)
    val      = eval(ex2call)
    ex.value = val
    return val
end

function evaluate(ex::A) where {A <: AbstractArray}
    s = similar(ex)
    for ii in eachindex(ex)
        s[ii] = evaluate(ex[ii])
    end

    return s
end

function substitute!(
    ex::SymbolicExpression,
    old_var::T,
    new_ex::U) where {
        T  <: Union{Constant, Variable, SymbolicExpression},
        U <: Union{Constant, Variable, SymbolicExpression}}

    for (ii, ee) in enumerate(ex.args)
        if ee isa Variable || ee isa Constant
            if ee == old_var
                ex.args[ii] = new_ex
            end
        elseif ee isa SymbolicExpression
            substitute!(ee, old_var, new_ex)
        else
            error("something went wrong")
        end
    end

    return ex
end

function clone(ex::SymbolicExpression)
    head = ex.head
    L    = length(ex.args)

    cloned_args = Vector{Any}(undef, L)

    for ii = 1:L
        cloned_args[ii] = clone(ex.args[ii])
    end
    y = SymbolicExpression(head, cloned_args...)

    # if L == 1
    #     cloned_arg = clone(ex.args[1])
    #     y = SymbolicExpression(head, cloned_arg)
    # elseif L == 2
    #     cloned_lhs = clone(ex.args[1])
    #     cloned_rhs = clone(ex.args[2])
    #     y = SymbolicExpression(head, cloned_lhs, cloned_rhs)
    # else
    #     error("something went wrong")
    # end

    return y
end

max_stars_repo_name: HomoModelicus/julia

max_issues_repo_name: HomoModelicus/julia

max_forks_repo_name: HomoModelicus/julia


text:
# Script to run the four formulations -- Bertsimas and Mišić, Belloni et al., McBride and Zufryden, and the Utility formulation -- in their integer and relaxed forms,
# for the Toubia data set.
# This script produces the results needed for Tables EC.1 and EC.2 in Section EC.3.1 of the electronic companion of
# Bertsimas and Mišić.
#
# For each value of n and numPermutations, the script randomly selects a set of n products out of the universe
# of 3584, and numPermutations customers out of the set of 330, and solves the corresponding PLD problem with
# those parameters.
# The script takes two command line arguments that specify the "rep range". To run 20 repetitions for
# each n and numPermutations value (as per the paper), type
#
# > julia Toubia_Comparison.jl 1 20
#
# in the terminal.

using JuMP, Gurobi, Stats, Distributions #, HDF5, JLD

include("../optimalPLD_code/optimalPLD_createPath.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_constraints.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_BFSS_constraints.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_McBZ_constraints.jl")
include("../optimalPLD_code/optimalPLD_solveSAO_Utility1_constraints.jl")

srand(1204)
machine = 2;

n_vec = [10, 20, 50, 100]
numPerm_vec = [10, 100] #, 330]

numReps = 100;

counter = 0;

LR = parse(ARGS[1]);
UR = parse(ARGS[2]);

expdirpath = optimalPLD_createPath(machine, "exp_dir", "Toubia_Comparison_reps$(LR)to$(UR)");
mkdir(expdirpath)

b_vec = [5; 10; 20; 100];

MIP_timeLimit = 1*3600;
isRelaxation = false;

data_prefix = "Toubia2003"

outcsvfilepath = string(expdirpath, "out.csv");
outcsvhandle = open(outcsvfilepath, "w")
row_string = string("n", ",",
					"numPermutations", ",",
					"rep", ",",
					"b_RHS", ",",
					"sumxi", ",",
					"BM_expRevenue", ",",
					"BM_upperBound", ",",
					"BM_solutiontime", ",",
					"BM_MIPgap", ",",
					"BM_LP_expRevenue", ",",
					"BM_LP_upperBound", ",",
					"BM_LP_solutiontime", ",",
					"BM_LP_MIPgap", ",",
					"BM_LP_isIntegral", ",",
					"BFSS_expRevenue", ",",
					"BFSS_upperBound", ",",
					"BFSS_solutiontime", ",",
					"BFSS_MIPgap", ",",
					"BFSS_LP_expRevenue", ",",
					"BFSS_LP_upperBound", ",",
					"BFSS_LP_solutiontime", ",",
					"BFSS_LP_MIPgap", ",",
					"BFSS_LP_isIntegral", ",",
					"Utility1_expRevenue", ",",
					"Utility1_upperBound", ",",
					"Utility1_solutiontime", ",",
					"Utility1_MIPgap", ",",
					"Utility1_LP_expRevenue", ",",
					"Utility1_LP_upperBound", ",",
					"Utility1_LP_solutiontime", ",",
					"Utility1_LP_MIPgap", ",",
					"Utility1_LP_isIntegral", ",",
					"McBZ_expRevenue", ",",
					"McBZ_upperBound", ",",
					"McBZ_solutiontime", ",",
					"McBZ_MIPgap", ",",
					"McBZ_LP_expRevenue", ",",
					"McBZ_LP_upperBound", ",",
					"McBZ_LP_solutiontime", ",",
					"McBZ_LP_MIPgap", ",",
					"McBZ_LP_isIntegral" )
row_string = string(row_string, "\n");

print(outcsvhandle, row_string);
flush(outcsvhandle);

rep_range = LR:UR

counter = 0;
for nind = 1:length(n_vec)
	n = n_vec[nind]
	for nPind = 1:length(numPerm_vec)
		numPermutations = numPerm_vec[nPind];

		println("**** Reading from data directory...")
		data_name = string(data_prefix,"_neq3584_Keq330_v2")
		data_dir = optimalPLD_createPath(machine, "data_dir", data_name)
		revenues_mat = readcsv(string(data_dir, "revenues_mat.csv"))
		lambda_mat = readcsv(string(data_dir, "lambda_mat.csv"))
		orderings_mat = readcsv(string(data_dir, "orderings_mat.csv"))
		utilities_mat = readcsv(string(data_dir, "utilities_mat.csv"))

		orderings_mat = convert(Array{Int64}, orderings_mat);

		for rep in rep_range

			prod_universe = randperm(3584)[1:n];
			prod_universe = sort(prod_universe);
			push!(prod_universe, n+1);
			cust_universe = randperm(330)[1:numPermutations];
			cust_universe = sort(cust_universe);

			revenues = revenues_mat[1,:][prod_universe]
			lambda = 1/numPermutations *ones(numPermutations); #[1:numPermutations]
			orderings_raw = orderings_mat[cust_universe,:]

			orderings = zeros(Int64, numPermutations, n+1);
			for k=1:numPermutations
				p2r = sortperm(orderings_raw[k,:]);
				p2r = p2r[prod_universe];
				orderings[k,:] = sortperm(p2r);
			end

			utilities = utilities_mat[cust_universe, prod_universe];
			for k=1:numPermutations
				utilities[k,:] -= minimum(utilities[k,:]);
			end

			writecsv(string(expdirpath, "prod_universe_neq", n, "_Keq", numPermutations, "rep_", rep, ".csv"), prod_universe)
			writecsv(string(expdirpath, "cust_universe_neq", n, "_Keq", numPermutations, "rep_", rep, ".csv"), cust_universe)

			b_vec_2 = filter(i -> i < n, b_vec)
			push!(b_vec_2, n);

			for b_feas_2 in b_vec_2

				b_feas = b_feas_2;
				if (b_feas_2 >= n)
					b_feas = Int64[];
				end

				A_feas = ones(1,n);
				x_initial = [];

				println("**** Solving n = $n, nP = $numPermutations, rep = $rep ****")
				println("**** Solving Bertsimas and Misic MIO ****")
				counter = counter + 1;
				isRelaxation = false;
				x_val, y_val, BM_expRevenue, BM_upperBound, BM_solutiontime, BM_MIPgap = optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas)

				sumxi = sum(x_val);

				println("**** Solving Bertsimas and Misic LP ****")
				isRelaxation = true;
				BM_LP_x_val, BM_LP_y_val, BM_LP_expRevenue, BM_LP_upperBound, BM_LP_solutiontime, BM_LP_MIPgap = optimalPLD_solveSAO_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				BM_LP_isIntegral = isinteger(BM_LP_x_val) ? 1 : 0;

				println("**** Solving Beloni MIO ****")
				isRelaxation = false;
				x_val, y_val, BFSS_expRevenue, BFSS_upperBound, BFSS_solutiontime, BFSS_MIPgap = optimalPLD_solveSAO_BFSS_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

				println("**** Solving Beloni LP ****")
				isRelaxation = true;
				BFSS_LP_x_val, BFSS_y_val, BFSS_LP_expRevenue, BFSS_LP_upperBound, BFSS_LP_solutiontime, BFSS_LP_MIPgap = optimalPLD_solveSAO_BFSS_constraints( n, revenues, numPermutations, lambda, orderings, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				BFSS_LP_isIntegral = isinteger(BFSS_LP_x_val) ? 1 : 0;

				println("**** Solving McBride Zufryden MIO ****")
				isRelaxation = false;
				x_val, y_val, McBZ_expRevenue, McBZ_upperBound, McBZ_solutiontime, McBZ_MIPgap = optimalPLD_solveSAO_McBZ_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

				println("**** Solving McBride Zufryden LP ****")
				isRelaxation = true;
				McBZ_LP_x_val, McBZ_LP_y_val, McBZ_LP_expRevenue, McBZ_LP_upperBound, McBZ_LP_solutiontime, McBZ_LP_MIPgap = optimalPLD_solveSAO_McBZ_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				McBZ_LP_isIntegral = isinteger(McBZ_LP_x_val) ? 1 : 0;

				println("**** Solving Utility1 MIO ****")
				isRelaxation = false;
				x_val, y_val, Utility1_expRevenue, Utility1_upperBound, Utility1_solutiontime, Utility1_MIPgap = optimalPLD_solveSAO_Utility1_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )

				println("**** Solving Utility1 LP ****")
				isRelaxation = true;
				Utility1_LP_x_val, Utility1_LP_y_val, Utility1_LP_expRevenue, Utility1_LP_upperBound, Utility1_LP_solutiontime, Utility1_LP_MIPgap = optimalPLD_solveSAO_Utility1_constraints( n, revenues, numPermutations, lambda, utilities, MIP_timeLimit, isRelaxation, x_initial, A_feas, b_feas )
				Utility1_LP_isIntegral = isinteger(Utility1_LP_x_val) ? 1 : 0;

				row_string = string(n, ",",
						numPermutations, ",",
						rep, ",",
						b_feas_2, ",",
						sumxi, ",",
						BM_expRevenue, ",",
						BM_upperBound, ",",
						BM_solutiontime, ",",
						BM_MIPgap, ",",
						BM_LP_expRevenue, ",",
						BM_LP_upperBound, ",",
						BM_LP_solutiontime, ",",
						BM_LP_MIPgap, ",",
						BM_LP_isIntegral, ",",
						BFSS_expRevenue, ",",
						BFSS_upperBound, ",",
						BFSS_solutiontime, ",",
						BFSS_MIPgap, ",",
						BFSS_LP_expRevenue, ",",
						BFSS_LP_upperBound, ",",
						BFSS_LP_solutiontime, ",",
						BFSS_LP_MIPgap, ",",
						BFSS_LP_isIntegral, ",",
						Utility1_expRevenue, ",",
						Utility1_upperBound, ",",
						Utility1_solutiontime, ",",
						Utility1_MIPgap, ",",
						Utility1_LP_expRevenue, ",",
						Utility1_LP_upperBound, ",",
						Utility1_LP_solutiontime, ",",
						Utility1_LP_MIPgap, ",",
						Utility1_LP_isIntegral, ",",
						McBZ_expRevenue, ",",
						McBZ_upperBound, ",",
						McBZ_solutiontime, ",",
						McBZ_MIPgap, ",",
						McBZ_LP_expRevenue, ",",
						McBZ_LP_upperBound, ",",
						McBZ_LP_solutiontime, ",",
						McBZ_LP_MIPgap, ",",
						McBZ_LP_isIntegral )
				row_string = string(row_string, "\n");

				print(outcsvhandle, row_string);
				flush(outcsvhandle);
			end
		end
	end
end

ext: jl

max_stars_repo_name: vvmisic/EFCPLO

max_stars_repo_licenses: MIT Unlicense

max_issues_repo_name: vvmisic/EFCPLO

max_issues_repo_licenses: MIT Unlicense

max_forks_repo_name: vvmisic/EFCPLO

max_forks_repo_licenses: MIT Unlicense


text:
"""
    AbstractLearner

An abstract learning mechanism. Defined by required interface functions and whatever
data structures your learning mechanism needs.

Required interface functions:
- `prespike!(learner::AbstractLearner, w::Real, t::Integer, src_id::Integer, dest_id::Integer; dt::Real = 1.0)`:
    processes a pre-synaptic spike at time `t` along the synapse from `src_id` to `dest_id`
- `postspike!(learner::AbstractLearner, w::Real, t::Integer, src_id::Integer, dest_id::Integer; dt::Real = 1.0)`:
    processes a post-synaptic spike at time `t` along the synapse from `src_id` to `dest_id`
- `update!(learner::AbstractLearner, w::Real, t::Integer, src_id::Integer, dest_id::Integer; dt::Real = 1.0)`:
    return the updated weight along the synapse from `src_id` to `dest_id` since the last call to `update!()`

This could be, for example, Hebbian learning. The functions `prespike!()` and `postspike!()`
are called whenever the pre-synaptic neuron or post-synaptic neuron fires, respectively.
When each function is called, it processes the spike according to the learning rule, and
updates an internal weight change state, `Δw`.
When `update!()` is called, the current `Δw` is returned then reset to zero.
"""
abstract type AbstractLearner end

"""
    George

A dumb learner that does nothing when processing spikes.
Aptly named after my dog, George (https://darsnack.github.io/website/about)
"""
struct George <: AbstractLearner end
prespike!(learner::George, w, spikes; dt::Real = 1.0) = return
postspike!(learner::George, w, spikes; dt::Real = 1.0) = return
prespike!(learner::George, src, dst, w, spikes; dt::Real = 1.0) = return
postspike!(learner::George, src, dst, w, spikes; dt::Real = 1.0) = return
update!(learner::George, w, t::Integer; dt::Real = 1.0) = w

"""
    STDP

A simple spike-timing-dependent plasticity mechanism.

Fields:
- `A₊::Real`: positive weight change amplitude
- `A₋::Real`: negative weight change amplitude
- `τ₊::Real`: positive decay parameter
- `τ₋::Real`: negative decay parameter
- `lastpre::Array{Float64, 2}`: the last occurence of a pre-synaptic spike
    (matrix row = src, col = dest)
- `lastpost::Array{Float64, 2}`: the last occurence of a post-synaptic spike
    (matrix row = src, col = dest)
- `Δw::Array{Float64, 2}`: the current synaptic weight change
    (matrix row = src, col = dest)
"""
struct STDP{T<:Real, VT<:AbstractArray{<:Real}} <: AbstractLearner
    A₊::T
    A₋::T
    τ₊::T
    τ₋::T
    lastpre::VT
    lastpost::VT
end

"""
    STDP(A₀::Real, τ::Real, n::Integer)

Create an STDP learner for `n` neuron population with weight change amplitude `A₀` and decay `τ`.
"""
STDP(A₀::Real, τ::Real, n::Integer) = STDP{Float32, Matrix{Float32}}(A₀, -A₀, τ, τ, zeros(n, n), zeros(n, n))

function prespike!(learner::STDP, w, spikes; dt::Real = 1.0)
    f(x, y, w) = (w != 0) && (y > 0) ? y : x
    @cast learner.lastpre[i, j] = f(learner.lastpre[i, j], spikes[i], w[i, j])
end

function postspike!(learner::STDP, w, spikes; dt::Real = 1.0)
    f(x, y, w) = (w != 0) && (y > 0) ? y : x
    @cast learner.lastpost[i, j] = f(learner.lastpost[i, j], spikes[j], w[i, j])
end

_stdpfpos(A, τ, t, dt, x, y) = (x > y && x == t) ? A * exp(-abs(x - y) * dt / τ) : zero(A)
_stdpfneg(A, τ, t, dt, x, y) = (x < y && y == t) ? A * exp(-abs(x - y) * dt / τ) : zero(A)

function update!(learner::STDP, w, t::Integer; dt::Real = 1.0)
    A₊, A₋, τ₊, τ₋ = learner.A₊, learner.A₋, learner.τ₊, learner.τ₋
    map!((x, y, w) -> w + _stdpfpos(A₊, τ₊, t, dt, x, y) + _stdpfneg(A₋, τ₋, t, dt, x, y),
         w, learner.lastpost, learner.lastpre, w)
end


text:
"""
Index-based Graphs
==================

Nodes
-----

A node (sometimes also *vertex* called) in a graph is always represented as an integer.
To be able to use a node as an index, the nodes should be within `1:n` whereby `n` is the number of nodes.
Furtheron, we always use **signed** integers because we might include extra information within the first (sign) bit.

Edges
-----

By convention, a graph is considered *undirected* (unless specified else).
Hence we only need to store one direction of every edge `(u, v)`; by default we assume `u < v`.

TODO:
- differences
- incidence matrix

If any *flow* `α``for an edge `(u, v)` has to be computed, the orientation is defined as going from `v` to `u` (big to small), i.e. we compute

```julia
x[u] += α
x[v] -= α
```

Submodules
----------

The most important submodules are
- [`Tree`](@ref): Rooted Trees (including minimum spanning trees).
- [`Grid`](@ref): Specialized code for grid graphs.
- [`LinA`](@ref): Linear algebra (e.g. incidence matrix)
"""
module GraphIdx

include("weights.jl")

include("incidence.jl")
include("biadjacent.jl")
include("cluster.jl")
include("Utils.jl")
include("Tree.jl")

include("graph.jl")
include("weightedgraph.jl")
include("edgegraph.jl")
include("prim_mst.jl")
include("kruskal_mst.jl")

include("bits.jl")
include("lina.jl")
include("grid.jl")
include("grididx.jl")
include("stats.jl")
include("Io.jl")

include("graphviz.jl")

include("precompile.jl")
_precompile_()

end

max_stars_repo_name: EQt/graphidx

max_issues_repo_name: EQt/graphidx

max_forks_repo_name: EQt/graphidx


text:
# - return: data with the dimension of `[lon, lat, ntime]`

"""
    nc_read(file; band = nothing, period = nothing, raw = false)

# Parameters
- `raw`: boolean. It `true`, not replace na values.
- `period`: `[year_start, year_end]`
"""
function nc_read(file, band = nothing; type = nothing, period = nothing, raw = false)

    ds = Dataset(path_mnt(file))
    if band === nothing
        band = nc_bands(file)[1]
    end

    # @time data = ds[band].var[:] # not replace na values at here
    data = raw ? ds[band].var : ds[band]
    @time if period === nothing
        data = data[:]
    else
        if length(period) == 1; period = repeat([period], 2); end
        # change NC into Raster
        dates = ds["time"]
        years = Dates.year.(dates)
        ind = (years .>= period[1] .&& years .<= period[2]) |> findall
        # `ind` is continuous, but reading speed is faster when converting to `unitRange`
        # https://alexander-barth.github.io/NCDatasets.jl/stable/performance/
        ind = ind[1]:ind[end]
        # @show ind
        data = data[:, :, ind]
        # data = length(dates) != length(ind) ?  data[:, :, ind] : data
    end
    close(ds)

    if type !== nothing && eltype(data) != type
        data = @.(type(data))
    end
    data
    # data = SharedArray(data)
end


text:
using Base: @propagate_inbounds, OneTo, RefValue

struct KeyedArray{T,N,AT,KT} <: AbstractArray{T,N}
    data::AT
    keys::KT
end

const KeyedVector{T,AT,RT} = KeyedArray{T,1,AT,RT}
const KeyedMatrix{T,AT,RT} = KeyedArray{T,2,AT,RT}
const KeyedVecOrMat{T,AT,RT} = Union{KeyedVector{T,AT,RT}, KeyedMatrix{T,AT,RT}}

function KeyedArray(data::AbstractArray{T,N}, keys::Tuple) where {T,N}
    construction_check(data, keys)
    KeyedArray{T, N, typeof(data), typeof(keys)}(data, keys)
end
function KeyedArray(data::AbstractVector{T}, keys::RefValue{<:AbstractVector}) where {T}
    construction_check(data, (keys[],))
    KeyedArray{T, 1, typeof(data), typeof(keys)}(data, keys)
end
KeyedArray(data::AbstractVector, tup::Tuple{AbstractVector}) =
    KeyedArray(data, Ref(first(tup)))
KeyedArray(data::AbstractVector, arr::AbstractVector) =
    KeyedArray(data, tuple(arr))

function construction_check(data::AbstractArray, keys::Tuple)
    length(keys) == ndims(data) || throw(ArgumentError(
        "wrong number of key vectors, got $(length(keys)) with ndims(A) == $(ndims(data))"))
    keys isa Tuple{Vararg{AbstractVector}} || throw(ArgumentError(
        "key vectors must all be AbstractVectors"))
    map(v -> axes(v,1), keys) == axes(data) || throw(ArgumentError(
        "lengths of key vectors must match size of array (and their axes too)"))
end

KeyedArray(A::KeyedArray, k2::Tuple) =
    KeyedArray(parent(A), unify_keys(axiskeys(A), k2))
KeyedArray(A::KeyedVector, k2::Tuple{AbstractVector}) =
    KeyedArray(parent(A), unify_keys(axiskeys(A), k2)) # avoids an ambiguity

Base.size(x::KeyedArray) = size(parent(x))

Base.axes(x::KeyedArray) = axes(parent(x))

Base.parent(x::KeyedArray) = getfield(x, :data)
keyless(x::KeyedArray) = parent(x)
keyless(x) = x

axiskeys(x::KeyedArray) = getfield(x, :keys)
axiskeys(x::KeyedVector) = tuple(getindex(getfield(x, :keys)))
# avoiding Tuple(Ref()) as it's slow, https://github.com/JuliaLang/julia/pull/33674
axiskeys(x::KeyedArray, d::Int) = d<=ndims(x) ? getindex(axiskeys(x), d) : OneTo(1)
axiskeys(x::KeyedVector, d::Int) = d==1 ? getindex(getfield(x, :keys)) : OneTo(1)

Base.IndexStyle(A::KeyedArray) = IndexStyle(parent(A))

Base.eachindex(A::KeyedArray) = eachindex(parent(A))

# Base.keys(A::KeyedArray) = error("Base.keys(::KeyedArray) not defined, please open an issue if this happens unexpectedly.") # https://github.com/mcabbott/AxisKeys.jl/issues/20
Base.keys(A::KeyedArray) = Base.keys(parent(A))

Base.pointer(A::KeyedArray, i::Integer) = Base.pointer(parent(A), i)

Base.unsafe_convert(::Type{Ptr{T}}, A::KeyedArray{T}) where {T} =
    Base.unsafe_convert(Ptr{T}, parent(A))

Base.strides(A::KeyedArray) = Base.strides(parent(A))
Base.stride(A::KeyedArray, d::Integer) = Base.stride(parent(A), d)

for (get_or_view, key_get, maybe_copy) in [
        (:getindex, :keys_getindex, :copy),
        (:view, :keys_view, :identity)
    ]
    @eval begin

        @inline function Base.$get_or_view(A::KeyedArray, raw_inds...)
            inds = to_indices(A, raw_inds)
            @boundscheck checkbounds(parent(A), inds...)
            data = @inbounds $get_or_view(parent(A), inds...)
            data isa AbstractArray || return data # scalar output

            raw_keys = $key_get(axiskeys(A), inds)
            raw_keys === () && return data # things like A[A .> 0]

            new_keys = ntuple(ndims(data)) do d
                raw_keys === nothing && return axes(data, d)
                raw_keys[d]
            end
            KeyedArray(data, new_keys)
        end

        # drop all, for A[:] and A[A .> 0] with ndims>=2
        @inline $key_get(keys::Tuple{Any, Any, Vararg{Any}}, inds::Tuple{Base.LogicalIndex}) = ()
        @inline $key_get(keys::Tuple{Any, Any, Vararg{Any}}, inds::Tuple{Base.Slice}) = ()

        # drop one, for integer index
        @inline $key_get(keys::Tuple, inds::Tuple{Integer, Vararg{Any}}) =
            $key_get(tail(keys), tail(inds))

        # from a Colon, k[:] would copy too, but this avoids view([1,2,3], :)
        @inline $key_get(keys::Tuple, inds::Tuple{Base.Slice, Vararg{Any}}) =
            ($maybe_copy(first(keys)), $key_get(tail(keys), tail(inds))...)

        # this avoids making views of 1:10 etc, they are immutable anyway
        @inline function $key_get(keys::Tuple, inds::Tuple{AbstractVector, Vararg{Any}})
            got = if first(keys) isa AbstractRange
                @inbounds getindex(first(keys), first(inds))
            else
                @inbounds $get_or_view(first(keys), first(inds))
            end
            (got, $key_get(tail(keys), tail(inds))...)
        end

        # newindex=[CartesianIndex{0}()], uses up one ind, sets N keys to default i.e. axes
        @inline function $key_get(keys::Tuple, inds::Tuple{AbstractVector{CartesianIndex{0}}, Vararg{Any}})
            (OneTo(1), $key_get(keys, tail(inds))...)
        end
        @inline function $key_get(keys::Tuple, inds::Tuple{AbstractVector{CartesianIndex{N}}, Vararg{Any}}) where {N}
            _, keys_left = Base.IteratorsMD.split(keys, Val(N))
            (ntuple(_->nothing, N)..., $key_get(keys_left, tail(inds))...)
        end

        # terminating case, trailing 1s (already checked) could be left over
        @inline $key_get(keys::Tuple{}, inds::Tuple{}) = ()
        @inline $key_get(keys::Tuple{}, inds::Tuple{Integer, Vararg{Any}}) = ()

    end
end

@inline function Base.setindex!(A::KeyedArray, val, raw_inds...)
    I = Base.to_indices(A, raw_inds)
    @boundscheck checkbounds(A, I...)
    @inbounds setindex!(parent(A), val, I...)
    val
end

# saves a few allocations in A[1,:,Not(3)] .= 0 compared to view()
@inline function Base.dotview(A::KeyedArray, raw_inds...)
    I = Base.to_indices(A, raw_inds)
    @boundscheck checkbounds(A, I...)
    @inbounds Base.dotview(parent(A), I...)
end


text:
using WavefrontObj, FileIO, MeshIO, GeometryTypes, GLVisualize, GLAbstraction
using Base.Test

obj  = read(file"cat.obj", GLNormalMesh)

function normals(verts, faces)
	normals_result = fill(Point3{Float32}(0), length(verts))
	for face in faces
		i1 = Int(face[1])
		i2 = Int(face[2])
		i3 = Int(face[3])

		v1 = verts[i1]
		v2 = verts[i2]
		v3 = verts[i3]
		a = v2 - v1
		b = v3 - v1
		n = cross(a,b)
		normals_result[i1] = n+normals_result[i1]
		normals_result[i2] = n+normals_result[i2]
		normals_result[i3] = n+normals_result[i3]
	end
	map(normalize, normals_result)
	convert(Normal3{Float32}, normals_result)
end
empty!(obj.attributes.normal)
append!(obj.attributes.normal, normals(obj.vertices, obj.faces))

obj.faces[:] = obj.faces - Uint32(1)
function funcy(x,y,z)
    Vec3(sin(x),cos(y),tan(z))
end

N = 15
directions  = Vec3[funcy(4x/N*3f0,4y/N,4z/N) for x=1:N,y=1:N, z=1:N]
robj 		= visualize(directions, primitive=obj)

push!(GLVisualize.ROOT_SCREEN.renderlist, robj)

renderloop()


text:
module MaBLAS

# general design philosophy:
#   1: keep the number of defined types in check
#   2: don't split code into another file/module if the API and tests are not standalone

include("loopinfo.jl")
include("gemm.jl")

end # module


text:
@testset "conversion and promotion" begin
    x1 = convert(RevisedAffineForm{3}, 1)
    @test x1 isa RevisedAffineForm{3,Float64}
    @test x1.c == 1
    @test iszero(x1.ϵ)
    @test iszero(x1.Δ)

    @test convert(RevisedAffineForm{3}, big(1)) isa RevisedAffineForm{3,BigFloat}

    x1 = convert(RevisedAffineForm{3,Float32}, big(1))
    @test x1 isa RevisedAffineForm{3,Float32}
    @test x1.c == 1
    @test iszero(x1.ϵ)
    @test iszero(x1.Δ)

    x1 = convert(RevisedAffineForm{3}, 1 .. 2)
    @test x1 isa RevisedAffineForm{3,Float64}
    @test x1.c == mid(1 .. 2)
    @test iszero(x1.ϵ)
    @test x1.Δ == radius(1 .. 2)

    @test convert(RevisedAffineForm{3}, big(1 .. 2)) isa RevisedAffineForm{3,BigFloat}

    x1 = convert(RevisedAffineForm{3,Float32}, big(1 .. 2))
    @test x1 isa RevisedAffineForm{3,Float32}
    @test x1.c == 1.5
    @test iszero(x1.ϵ)
    @test x1.Δ == 0.5

    @test promote_type(RevisedAffineForm{2,Float64}, RevisedAffineForm{2,BigFloat}) ==
        RevisedAffineForm{2,BigFloat}

    @test promote_type(RevisedAffineForm{2,Float64}, BigInt) ==
        RevisedAffineForm{2,BigFloat}

    @test promote_type(RevisedAffineForm{2,Float64}, Interval{Rational{Int}}) ==
        RevisedAffineForm{2,Float64}
end


text:
using LinearAlgebra

A = [1 2 3; 4 5 6; 7 8 9]

for t in [
Symmetric
Hermitian
UpperTriangular
LowerTriangular
Tridiagonal
Diagonal
]
    show(stdout, MIME("text/plain"), t(A))
    println()
end

A = [1, 2, 3]
B = [5, 6]
for t in [
SymTridiagonal
]
    show(stdout, MIME("text/plain"), t(A, B))
    println()
end

A = [1, 2, 3]
B = [5, 6]
for t in [
Bidiagonal
]
    show(stdout, MIME("text/plain"), t(A, B, :U))
    println()
    show(stdout, MIME("text/plain"), t(A, B, :L))
    println()
end

for t in [
UniformScaling
]
    println(t(2))
end

#=
3×3 Symmetric{Int64,Array{Int64,2}}:
 1  2  3
 2  5  6
 3  6  9
3×3 Hermitian{Int64,Array{Int64,2}}:
 1  2  3
 2  5  6
 3  6  9
3×3 UpperTriangular{Int64,Array{Int64,2}}:
 1  2  3
 ⋅  5  6
 ⋅  ⋅  9
3×3 LowerTriangular{Int64,Array{Int64,2}}:
 1  ⋅  ⋅
 4  5  ⋅
 7  8  9
3×3 Tridiagonal{Int64,Array{Int64,1}}:
 1  2  ⋅
 4  5  6
 ⋅  8  9
3×3 Diagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 ⋅  5  ⋅
 ⋅  ⋅  9
3×3 SymTridiagonal{Int64,Array{Int64,1}}:
 1  5  ⋅
 5  2  6
 ⋅  6  3
3×3 Bidiagonal{Int64,Array{Int64,1}}:
 1  5  ⋅
 ⋅  2  6
 ⋅  ⋅  3
3×3 Bidiagonal{Int64,Array{Int64,1}}:
 1  ⋅  ⋅
 5  2  ⋅
 ⋅  6  3
UniformScaling{Int64}
2*I
=#


text:
using WoodburyMatrices
using Base.Test

srand(123)
n = 5

d = 1 .+ rand(n)
dl = -rand(n-1)
du = -rand(n-1)
v = randn(n)
B = randn(n,2)

U = randn(n,2)
V = randn(2,n)
C = randn(2,2)

for elty in (Float32, Float64, Complex64, Complex128, Int)
    if elty == Int
        srand(61516384)
        d = rand(1:100, n)
        dl = -rand(0:10, n-1)
        du = -rand(0:10, n-1)
        v = rand(1:100, n)
        B = rand(1:100, n, 2)

        # Woodbury
        U = rand(1:100, n, 2)
        V = rand(1:100, 2, n)
        C = rand(1:100, 2, 2)
    else
        d = convert(Vector{elty}, d)
        dl = convert(Vector{elty}, dl)
        du = convert(Vector{elty}, du)
        v = convert(Vector{elty}, v)
        B = convert(Matrix{elty}, B)
        U = convert(Matrix{elty}, U)
        V = convert(Matrix{elty}, V)
        C = convert(Matrix{elty}, C)
    end
    ε = eps(abs2(float(one(elty))))
    T = Tridiagonal(dl, d, du)
    @test size(T, 1) == n
    @test size(T) == (n, n)

    # Woodbury
    W = Woodbury(T, U, C, V)
    F = full(W)
    @test_approx_eq W*v F*v
    iFv = F\v
    @test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε # Revisit. Condition number is wrong
    @test abs((det(W) - det(F))/det(F)) <= n*cond(F)*ε # Revisit. Condition number is wrong
    iWv = similar(iFv)
    if elty != Int
        iWv = A_ldiv_B!(W, copy(v))
        @test_approx_eq iWv iFv
    end
end

# Sparse U and V
n = 5
d = 1 .+ rand(n)
dl = -rand(n-1)
du = -rand(n-1)
T = Tridiagonal(dl, d, du)
U = sprandn(n,2,0.2)
V = sprandn(2,n,0.2)
C = randn(2,2)
W = Woodbury(T, U, C, V)
F = full(W)

v = randn(n)
ε = eps()
iFv = F\v
@test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε

# Display
iob = IOBuffer()
show(iob, W)

# Vector-valued U and scalar-valued C
U = rand(n)
C = 2.3
V = rand(1,n)
W = Woodbury(T, U, C, V)
F = full(W)

v = randn(n)
ε = eps()
iFv = F\v
@test norm(W\v - iFv)/norm(iFv) <= n*cond(F)*ε

show(iob, W)

Wc = copy(W)

# Mismatched sizes
@test_throws DimensionMismatch Woodbury(rand(5,5),rand(5,2),rand(2,3),rand(3,5))
@test_throws DimensionMismatch Woodbury(rand(5,5),rand(5,2),rand(3,3),rand(2,5))

# spec builder
n = 10
num_nz = 3
for i in 1:5  # repeat 5 times
    args = [(rand(1:n), rand(1:n), rand()) for j in 1:num_nz]
    r, v, c = WoodburyMatrices.sparse_factors(Float64, n, args...)
    out = r*v*c

    by_hand = zeros(n, n)
    for (i, j, v) in args
        by_hand[i, j] = v
    end

    @test maxabs(out - by_hand) == 0.0
end

include("runtests_sym.jl")


text:
using Revise
using PyPlot, DelimitedFiles, Random, Statistics, transD_GP

Random.seed!(23)

zfixed = [-1e5]
ρfixed = [1e12]
nmax = 100

zstart = 0.0
extendfrac, dz = 1.06, 2.
zall, znall, zboundaries = transD_GP.setupz(zstart, extendfrac, dz=dz, n=40)
z, ρ, nfixed = transD_GP.makezρ(zboundaries; zfixed=zfixed, ρfixed=ρfixed)
@info z

## geometry parameters for tempest
zTx = -120
zRx = -82
x_rx = -113.0
y_rx = 0.
rx_roll = 0.
rx_pitch = 2.
rx_yaw = 0.
tx_roll = 0.
tx_pitch = 0.
tx_yaw = 0.
# electronics and stuff
include("electronics_halt.jl")
## fill in detail in ohm-m
ρ[(z.>=zstart) .& (z.<50)] .= 20.
ρ[(z.>=50) .&(z.<80)] .= 1.
ρ[(z.>=80) .&(z.<100)] .= 20
ρ[(z.>=100) .&(z.<200)] .= 50
ρ[(z.>=200) .&(z.<250)] .= 80
ρ[(z.>=250)]            .= 150

## create total field operator (required for nuisance inversion)
tempest = transD_GP.TEMPEST1DInversion.Bfield(
    zTx = zTx, zRx = zRx, x_rx = x_rx, y_rx = y_rx,
    rx_roll = rx_roll, rx_pitch = rx_pitch, rx_yaw = rx_yaw,
    tx_roll = tx_roll, tx_pitch = tx_pitch, tx_yaw = tx_yaw,
	ramp = ramp, times = times,
	z=z,
	ρ=ρ,
	addprimary = true #this ensures that the geometry update actually changes everything that needs to be
)
# plot before adding noise
transD_GP.TEMPEST1DInversion.plotmodelfield!(tempest,z,ρ)
## compute noisy data to invert
# remember noise in electronics_halt.jl are in fT!!
transD_GP.TEMPEST1DInversion.set_noisy_data!(tempest, z, ρ,
                        noisefracx = 0.02, noisefracz = 0.02,
                        halt_X = Hx_add_noise*1e-15, halt_Z = Hz_add_noise*1e-15)


text:
# More accurately, this should be called `BoundedUnivariateDensity`
struct UnivariateDensity{F<:Function, T<:Real}
    f::F
    lb::T
    ub::T
    function UnivariateDensity{F,T}(f::F, lb::T, ub::T) where {F<:Function, T<:Real}
        area = quadgk(f, lb, ub)[1]
        minf = Optim.minimum(Optim.optimize(f, lb, ub, Brent()))
        if area ≈ 1.0 && lb < ub && (minf ≥ 0 || minf ≈ 0.0)
            g(t) = lb ≤ t ≤ ub ? f(t) : zero(t)
            return new{typeof(g), T}(g, lb, ub)
        else
            throw(ArgumentError("Not a valid density."))
        end
    end
end

UnivariateDensity(f::F, lb::T, ub::T) where {F<:Function, T <: Real} = UnivariateDensity{F, T}(f, lb, ub)
UnivariateDensity(f::Function, lb::Real, ub::Real) = UnivariateDensity(f, promote(lb, ub)...)
UnivariateDensity(f::Real, lb, ub) = UnivariateDensity(_ -> 1 / (ub - lb), lb, ub)

function integrate(g::Function, d::UnivariateDensity, bounds::Vararg{T}) where {T <: Real}
    integrand(t) = g(t) * d.f(t)
    if length(bounds) == 0
        return quadgk(integrand, d.lb, d.ub)
    elseif length(bounds) == 1
        return quadgk(integrand, bounds[1], d.ub)
    else
        return quadgk(integrand, bounds[1], bounds[2])
    end
end

max_stars_repo_name: carlocab/BPPrivInfo

max_issues_repo_name: carlocab/BPPrivInfo

max_forks_repo_name: carlocab/BPPrivInfo


text:
#using Libdl
#using Unitful
using DifferentialEquations
using .Scattering: Method, Potential

struct VPA <: Method
    ρ::Float64  # Radius for where to truncate the potential
    rspan::Tuple{Float64, Float64}  # Range of the interaction
end
VPA() = VPA(1.0, (1e-4, 15.0))
VPA(ρ::Real) = VPA(ρ, (1e-4, 15.0))
VPA(rspan::Tuple{Float64, Float64}) = VPA(1.0, rspan)

function (vpa::VPA)(k, m, V::Potential;kw...)
    #Vρ(r) = if r ≤ vpa.ρ V(r) else 0 end
    #m = m/2

    function dδdr(δ, p, r)
        -1.0/k * 2m * V(r) * sin(k*r + δ)^2
    end

    δ0 = 0.0

    prob = ODEProblem(dδdr, δ0, vpa.rspan)
    sol = solve(prob, reltol=1e-8; kw...)

    δ = sol.u[end]
    #kcot = k/tan(δ)
    #δ/ pi * 180
    δ
end

max_stars_repo_name: ErlendLima/Scattering

max_issues_repo_name: ErlendLima/Scattering

max_forks_repo_name: ErlendLima/Scattering


text:
using OSQP
using Compat
using Compat.Test, Compat.SparseArrays, Compat.LinearAlgebra, Compat.Random

tests = [
    "basic.jl",
    "dual_infeasibility.jl",
    "feasibility.jl",
    "polishing.jl",
    "primal_infeasibility.jl",
    "unconstrained.jl",
    "warm_start.jl",
    "update_matrices.jl",
    "mpbinterface.jl",
    "MathOptInterfaceOSQP.jl"
    ]

println("Running tests:")

for curtest in tests
    println(" Test: $(curtest)")
    include(curtest)
end


text:
module gRPCClient

using LibCURL
using Downloads
using ProtoBuf

import Downloads: Curl
import ProtoBuf: call_method

export gRPCController, gRPCChannel, gRPCException, gRPCServiceCallException, gRPCMessageTooLargeException, gRPCStatus, gRPCCheck, StatusCode

abstract type gRPCException <: Exception end

include("limitio.jl")
include("curl.jl")
include("grpc.jl")
include("generate.jl")

function __init__()
    GRPC_STATIC_HEADERS[] = grpc_headers()
end

end # module


text:
@doc doc"""Represents a NeoVim tabpage.""" ->
type Tabpage
    nvim::Nvim
    ext::MsgPack.Ext
end


text:
# Autogenerated wrapper script for RAPID_jll for x86_64-linux-gnu
export rapid

using NetCDF_jll
using NetCDFF_jll
using PETSc_jll
JLLWrappers.@generate_wrapper_header("RAPID")
JLLWrappers.@declare_executable_product(rapid)
function __init__()
    JLLWrappers.@generate_init_header(NetCDF_jll, NetCDFF_jll, PETSc_jll)
    JLLWrappers.@init_executable_product(
        rapid,
        "bin/rapid",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
using Test, Statistics

pdir = joinpath(@__DIR__, "..", "examples", "mcmcchains")
cd(pdir) do

  isdir("tmp") && rm("tmp", recursive=true);

  include(joinpath(pdir, "mcmcchains.jl"))

  @test 18.0 <  mean(chn.value) < 22.0

  isdir("tmp") && rm("tmp", recursive=true);

end


text:
using SIIPExamples

pkgpath = dirname(dirname(pathof(SIIPExamples)))
include(joinpath(pkgpath, "test", "2_PowerSystems_examples", "02_parse_matpower.jl"))

folder = mktempdir()
path = joinpath(folder, "system.json")
println("Serializing to $path")
to_json(sys, path)

filesize(path) / (1024 * 1024) #MiB

sys2 = System(path)

# This file was generated using Literate.jl, https://github.com/fredrikekre/Literate.jl


text:
# Autogenerated wrapper script for NLopt_jll for i686-linux-musl-cxx11
export libnlopt

JLLWrappers.@generate_wrapper_header("NLopt")
JLLWrappers.@declare_library_product(libnlopt, "libnlopt.so.0")
function __init__()
    JLLWrappers.@generate_init_header()
    JLLWrappers.@init_library_product(
        libnlopt,
        "lib/libnlopt.so",
        RTLD_LAZY | RTLD_DEEPBIND,
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
@testset "arithmetic.jl" begin
    @testset "projection" begin
        for i in 1:5, j in 1:5, k in 1:5
            @test p1(i, j, k) == i
            @test p2(i, j, k) == j
            @test p3(i, j, k) == k
        end
    end
    @testset "addition" begin
        for i in 1:10, j in 1:10
            @test add(i, j) == i + j
        end
    end
    @testset "multi" begin
        for i in 1:10, j in 1:10
            @test mult(i, j) == i * j
        end
    end
    @testset "pred" begin
        for i in 1:10
            @test pred(i) == max(0, i - 1)
        end
    end
    @testset "cosub" begin
        for i in 1:10, j in rand(1:10, 10)
            @test cosub(i, j) == max(0, i - j)
        end
    end
    @testset "diff" begin
        for i in 1:10, j in rand(1:10, 10)
            @test dff(i, j) == abs(i - j)
        end
    end
    @testset "const1" begin
        for i in 1:10
            @test const1(i) == 1
        end
    end
    @testset "reminder" begin
        for i in 0:10, j in 1:5
            @test remainder(i, j) == i % j
        end
    end
    @testset "square" begin
        for i in 0:10
            @test square(i) == i ^ 2
        end
    end
    @testset "square_root" begin
        for i in 0:10
            @test floor_sqrt(i) == floor(Int, sqrt(i))
        end
    end
    @testset "traingle" begin
        for i in 1:10
            @test traingle(i) == binomial(i, 2)
        end
    end
    @testset "pair" begin
        for i in 0:8, j in 0:8
            @test pair(i, j) == binomial(i + j + 1, 2) + i
        end
    end
    @testset "inv-pair" begin
        for i in 0:6, j in 0:6-i
            @test fst(binomial(i + j + 1, 2) + i) == i
            @test scd(binomial(i + j + 1, 2) + i) == j
        end
    end

    @testset "fib" begin
        fibs = [0, 1, 1, 2, 3, 5, 8]
        for i in 0:5
            @test fibs[i + 1] == fib(i)
        end
    end
end


text:
module MaxHelpingHandSeekerBarrierColorControllerDisabler

using ..Ahorn, Maple

@mapdef Entity "MaxHelpingHand/SeekerBarrierColorControllerDisabler" SeekerBarrierColorControllerDisabler(x::Integer, y::Integer)

const placements = Ahorn.PlacementDict(
    "Seeker Barrier Color Controller (Disable)\n(max480's Helping Hand)" => Ahorn.EntityPlacement(
        SeekerBarrierColorControllerDisabler
    )
)

function Ahorn.selection(entity::SeekerBarrierColorControllerDisabler)
    x, y = Ahorn.position(entity)

    return Ahorn.Rectangle(x - 12, y - 12, 24, 24)
end

Ahorn.render(ctx::Ahorn.Cairo.CairoContext, entity::SeekerBarrierColorControllerDisabler, room::Maple.Room) = Ahorn.drawSprite(ctx, "ahorn/MaxHelpingHand/rainbowSpinnerColorControllerDisable", 0, 0)

end

max_stars_repo_name: Kitty-Cats/MaxHelpingHand

max_issues_repo_name: Kitty-Cats/MaxHelpingHand

max_forks_repo_name: Kitty-Cats/MaxHelpingHand


text:
function OrientFaceCart(n,m)
Orient=n[3];
return Orient
end


text:
#######################################################################
# Tools
#######################################################################
function meandist(h::Hamiltonian)
    distsum = 0.0
    num = 0
    ss = Elsa.sites(h.lattice)
    br = h.lattice.bravais.matrix
    for (row, col, dn) in Elsa.eachindex_nz(h)
        if row != col
            num += 1
            rsrc = ss[col]
            rdst = ss[row] + br * dn
            distsum += norm(rsrc - rdst)
        end
    end
    return iszero(num) ? 0.0 : distsum / num
end

function matrixidx(h::AbstractSparseMatrix, row, col)
    for ptr in nzrange(h, col)
        rowvals(h)[ptr] == row && return ptr
    end
    return 0
end

matrixidx(h::DenseMatrix, row, col) = LinearIndices(h)[row, col]

transparent(rgba::RGBAf0, v = 0.5) = RGBAf0(rgba.r, rgba.g, rgba.b, rgba.alpha * v)

function darken(rgba::RGBAf0, v = 0.66)
    r = max(0, min(rgba.r * (1 - v), 1))
    g = max(0, min(rgba.g * (1 - v), 1))
    b = max(0, min(rgba.b * (1 - v), 1))
    RGBAf0(r,g,b,rgba.alpha)
end
function lighten(rgba, v = 0.66)
    darken(rgba, -v)
end


text:
using Gadfly, DataFrames, Cairo

function read_bm_transforms(fname, framework)
    isfile(fname) || error("given file location $(fname) is invalid")
    data = DataFrame(
        category = String[],
        device = String[],
        transforms = String[],
        npoints = Int[],
        time_ms = Float64[],
    )
    for line in eachline(fname)
        raw = split(line)
        push!(data, (raw[1], raw[2], raw[3], parse(Int, raw[4]), parse(Float64, raw[5])))
    end
    data[!, :framework] .= framework
    return data
end

function read_bm_metrics(fname)
    isfile(fname) || error("given file location $(fname) is invalid")
    data = DataFrame(
        framework = String[],
        category = String[],
        device = String[],
        transforms = String[],
        npoints = Int[],
        time_ms = Float64[],
    )
    for line in eachline(fname)
        raw = split(line)
        push!(
            data,
            (raw[1], raw[2], raw[3], raw[4], parse(Int, raw[5]), parse(Float64, raw[6])),
        )
    end
    return data
end

function save_benchmarks(fname, bm, variant, xlabel)
    # using Theme(background_color = colorant"white")
    # as arg in plot to force white background
    p = plot(
        bm,
        xgroup = "transforms",
        ygroup = "device",
        color = "framework",
        x = "npoints",
        y = "time_ms",
        Guide.ylabel("Time (milliseconds)"),
        Guide.xlabel(xlabel),
        Theme(background_color = colorant"white"),
        Geom.subplot_grid(Geom.point, Geom.line, Scale.x_log2, Scale.y_log10),
    )
    draw(PNG(fname, 40cm, 20cm), p)
end

bm_flux3d = read_bm_transforms(joinpath(@__DIR__, "bm_flux3d.txt"), "Flux3D.jl")
bm_kaolin = read_bm_transforms(joinpath(@__DIR__, "bm_kaolin.txt"), "Kaolin")
bm = vcat(bm_flux3d, bm_kaolin)

bm_flux3d_metrics = read_bm_metrics(joinpath(@__DIR__, "bm_flux3d_metrics.txt"))
bm_kaolin_metrics = read_bm_metrics(joinpath(@__DIR__, "bm_kaolin_metrics.txt"))
bm_metrics = vcat(bm_flux3d_metrics, bm_kaolin_metrics)

ispath(joinpath(@__DIR__, "plots")) || mkdir(joinpath(@__DIR__, "plots"))

save_benchmarks(
    joinpath(@__DIR__, "plots/bm_pcloud.png"),
    bm[bm[!, :category].=="PointCloud", :],
    "PointCloud",
    "No. of points in PointCloud",
)
save_benchmarks(
    joinpath(@__DIR__, "plots/bm_trimesh.png"),
    bm[bm[!, :category].=="TriMesh", :],
    "TriMesh",
    "No. of verts in TriMesh",
)
save_benchmarks(
    joinpath(@__DIR__, "plots/bm_metrics.png"),
    bm_metrics,
    "Metrics",
    "No. of verts in TriMesh",
)


text:
# fitting

"""
    struct FittedSphere{A<:AbstractArray, R<:Real} <: FittedShape

Sphere primitive, defined by its center and radius.
Also stored, if the normals point outwards of the shape.
"""
struct FittedSphere{A<:AbstractArray, R<:Real} <: FittedShape
    center::A
    radius::R
    outwards::Bool
end

Base.show(io::IO, x::FittedSphere) =
    print(io, """sphere, R: $(x.radius)""")

Base.show(io::IO, ::MIME"text/plain", x::FittedSphere{A, R}) where {A, R} =
    print(io, "FittedSphere{$A, $R}\n", """center: $(x.center), R: $(x.radius), $(x.outwards ? "outwards" : "inwards" )""")

strt(x::FittedSphere) = "sphere"

function defaultshapeparameters(::Type{FittedSphere})
    # `sphere_par`: parameter in sphere fitting
    sp = (ϵ=0.3, α=deg2rad(5), sphere_par=0.02)
    return (sphere=sp,)
end

function fit2pointsphere(v, n, params)
    #@unpack parallelthrdeg, sphere_par = params
    @extract params : params_sphere=sphere
    @extract params_sphere : sphere_par
    @extract params.common : parallelthrdeg
    n1n = normalize(n[1])
    n2n = normalize(n[2])

    b = false
    if abs(dot(n1n, n2n)) > cosd(parallelthrdeg)
        # parallel normals
        # fit it, if bad, will fall out later
        centerp = (v[1]+v[2])/2
        return FittedSphere(centerp, norm(centerp-v[1]), b)
    else
        # not parallel normals
        # first check if they intersect
        g = v[2]-v[1]
        h = cross(n2n, g)
        k = cross(n2n, n1n)
        nk = norm(k)
        nh = norm(h)

        if nk < sphere_par || nh < sphere_par
            # if nk == 0 || nh == 0, but it's numeric
            # no intersection
            n2 = n2n × (n1n × n2n)
            n1 = n1n × (n2n × n1n)
            c1 = v[1] + dot( (v[2]-v[1]), n2 )/dot( n[1], n2 ) * n[1]
            c2 = v[2] + dot( (v[1]-v[2]), n1 )/dot( n[2], n1 ) * n[2]
            centerp = (c1+c2)/2
            r = (norm(v[1]-centerp) + norm(v[1]-centerp))/2
            return FittedSphere(centerp, r, b)
        else
            # intersection
            if dot(h, k) > 0
                # point the same direction -> +
                M = v[1] + nh/nk * n1n
                return FittedSphere(M, norm(M-v[1]), b)
            else
                # point in different direction -> -
                M = v[1] - nh/nk * n1n
                return FittedSphere(M, norm(M-v[1]), b)
            end
        end
    end
end

function setsphereOuterity(sp, b)
    FittedSphere(sp.center, sp.radius, b)
end

"""
    fit(::Type{FittedSphere}, p, n, pc, params)

Fit a sphere to 2 points. Additional points and their normals are used to validate the fit.
Return `nothing` if points do not fit to a sphere.
"""
function fit(::Type{FittedSphere}, p, n, pc, params)
    #@unpack ϵ_sphere, α_sphere = params
    @extract params : params_sphere=sphere
    @extract params_sphere : α_sphere=α ϵ_sphere=ϵ
    pl = length(p)
    @assert pl == length(n) "Size must be the same."
    @assert pl > 2 "Size must be at least 3."
    # "forcefit" a sphere
    sp = fit2pointsphere(p, n, params)
    # check if real sphere
    sp === nothing && return nothing
    thr = cos(α_sphere)
    vert_ok = falses(pl)
    norm_ok = falses(pl)
    invnorm_ok = falses(pl)
    for i in 1:pl
        # vertice check
        vert_ok[i] = abs(norm(p[i]-sp.center)-sp.radius) < ϵ_sphere
        # normal check
        dotp = dot( normalize(p[i]-sp.center), normalize(n[i]) )
        norm_ok[i] = dotp > thr
        invnorm_ok[i] = dotp < -thr
    end
    vert_ok == trues(pl) || return nothing
    norm_ok == trues(pl) && return setsphereOuterity(sp, true)
    invnorm_ok == trues(pl) && return setsphereOuterity(sp, false)
    return nothing
end

# bitmapping

function scorecandidate(pc, candidate::FittedSphere, subsetID, params)
    ps = @view pc.vertices[pc.subsets[subsetID]]
    ns = @view pc.normals[pc.subsets[subsetID]]
    ens = @view pc.isenabled[pc.subsets[subsetID]]

    cpl = compatiblesSphere(candidate, ps, ns, params)
    # verti: összes pont indexe, ami enabled és kompatibilis
    # lenne, ha működne, de inkább a boolean indexelést machináljuk
    verti = pc.subsets[subsetID]
    #underEn = uo.under .& cpl
    #overEn = uo.over .& cpl

    #inpoints = count(underEn) >= count(overEn) ? verti[underEn] : verti[overEn]
    inpoints = verti[cpl]
    score = estimatescore(length(pc.subsets[subsetID]), pc.size, length(inpoints))
    return (score, inpoints)
end

"""
    compatiblesSphere(plane, points, normals, eps, alpharad)

Create a bool-indexer array for those points that are compatible to the sphere.

Compatibility is measured with an `eps` distance to the sphere
and an `alpharad` angle to it's normal.
"""
function compatiblesSphere(sphere, points, normals, params)
    #@unpack ϵ_sphere, α_sphere = params
    @extract params : params_sphere=sphere
    @extract params_sphere : α ϵ
    @assert length(points) == length(normals) "Size must be the same."
    # eps check
    o = sphere.center
    R = sphere.radius
    # c1 = [abs(norm(a-o)-R) < ϵ_sphere for a in points]
    # alpha check
    #if sphere.outwards
    #    c2=[isparallel(normalize(points[i]-o), normals[i], α_sphere) && c1[i] for i in eachindex(points)]
    #else
    #    c2=[isparallel(normalize(o-points[i]), normals[i], α_sphere) && c1[i] for i in eachindex(points)]
    #end

    zpn = zip(points, normals)

    if sphere.outwards
        c2=[isparallel(normalize(p-o), n, α) && (abs(norm(p-o)-R) < ϵ) for (p,n) in zpn]
    else
        c2=[isparallel(normalize(o-p), n, α) && (abs(norm(p-o)-R) < ϵ) for (p,n) in zpn]
    end

    # TODO
    # normalize by it's own length or maximum length??? r+3ϵ
    # param_points = [unitdisk2square(normalize(a[1:2]-o[1:2])) for a in points]
    return c2
end

"""
    refit(s::T, pc, params) where {T<:FittedSphere}

Refit sphere.
"""
function refit(s::T, pc, params) where {T<:FittedSphere}
    # TODO: use octree for that
    pcv = @view pc.vertices[pc.isenabled]
    pcn = @view pc.normals[pc.isenabled]
    cpl = compatiblesSphere(s, pcv, pcn, params)
    # verti: összes pont indexe, ami enabled és kompatibilis
    verti = (1:pc.size)[pc.isenabled]
    #underEn = uo.under .& cpl
    #overEn = uo.over .& cpl
    #s.inpoints = count(underEn) >= count(overEn) ? verti[underEn] : verti[overEn]
    return ExtractedShape(s, verti[cpl])
end


text:
using Pidoh
using Graphs, SimpleWeightedGraphs
using Plots, DataFrames, PlotThemes, Statistics, DataFrames, LaTeXStrings
using StatsPlots

function create_experiment_TG(workspace, rep)
    algorithms :: Array{Any, 1} = []
    ins::Array{Instance, 1} = []
    for n in repeat([ 24, 36, 48,60],rep)
    # for n in repeat([ 10, 20],rep)
        mst = MST(TG(n))
        m = ne(mst.g)
        append!(algorithms, [
                RLSSDstar(stop=FixedBudget(10^9) , R=m^5, name=L"SD-RLS$^r$"),
                # RLSSDstarstar(stop=FixedBudget(10^9), R=m^4),
                ea1p1(stop=FixedBudget(10^9), mutation=UniformlyIndependentMutation(1/m), name = LaTeXString("(1+1)EA")),
                ea1p1(stop=FixedBudget(10^9), mutation=HeavyTailedMutation(1.5,m), name=L"(1+1)FEA$_{\beta}, \beta=1.5$"),
                RLS12(stop=FixedBudget(10^9), name=L"RLS$^{1,2}$"),
                ]
                )
        append!(ins, [Instance(mst, generator=RandBitStringIP(m), name=L"TG") for i in 1:4])
    end

    Experiment(workspace, Array{AbstractAlgorithm}(algorithms), ins, repeat=1)
end

function create_experiment_ER(workspace, rep)
    algorithms :: Array{Any, 1} = []
    ins::Array{Instance, 1} = []
    for n in repeat([40, 60,80,100], rep)
    # for n in repeat([10, 20, 30, 40], rep)
        mst = MST(ER(n))
        m = ne(mst.g)
        append!(algorithms, [
                RLSSDstar(stop=FixedBudget(10^9) , R=m^5, name=L"SD-RLS$^r$"),
                # RLSSDstarstar(stop=FixedBudget(10^9), R=m^4),
                ea1p1(stop=FixedBudget(10^9), mutation=UniformlyIndependentMutation(1/m), name = LaTeXString("(1+1)EA")),
                ea1p1(stop=FixedBudget(10^9), mutation=HeavyTailedMutation(1.5,m), name=L"(1+1)FEA$_{\beta}, \beta=1.5$"),
                RLS12(stop=FixedBudget(10^9), name=L"RLS$^{1,2}$"),
                ]
                )
        append!(ins, [Instance(mst, generator=RandBitStringIP(m), name=L"ER") for i in 1:4])
    end

    Experiment(workspace, Array{AbstractAlgorithm}(algorithms), ins, repeat=1)
end

exp_er = create_experiment_ER("ecj/er_ecj", 400)
exp_tg = create_experiment_TG("ecj/tg_ecj", 400)

run(exp_er)

run(exp_tg)

data_er = runtimes(exp_er)
data_tg = runtimes(exp_tg)

show(data, allcols=true)
data.ProblemParam_g

categories = :algorithm
xdata = :ProblemParam_g_weights_n
ydata = :runtime
data = data_tg

Plots.PyPlotBackend()

Plots.font("DejaVu")

begin
    sort!(data, :algorithm, rev=false)
    plt = plot()
    for cat in keys(groupby(data, categories))
        filtered_data = filter(row -> row[categories]==cat[1], data)
        aggregated_data = combine(groupby(filtered_data, xdata), ydata .=> [mean,std])
        label = L"^"*cat[1]
        plot!(plt, aggregated_data[!,1], aggregated_data[!,2], fillalpha=0.5, dpi=500,
            linewidth = 2,
            markersize=3, markershape=:circle,linealpha = 0.5, size=(320,320), yscale=:log10,
            label = label, xlabel=L"Number of vertices $(n)$", ylabel=string("Fitness function calls"),
            legend= :bottomright,
            xticks = unique(data.ProblemParam_g_weights_n)
        )
    end
    # savefig(plt, "er_line.png")
    savefig(plt, "tg_line.png")
    plt
end

begin
    sort!(data, :runtime, rev=false)
    boxplt = plot()
    groupedboxplot!(boxplt,  data.ProblemParam_g_weights_n , data.runtime, group = L"^".*data.AlgorithmParams_name_s, bar_width = 0.8,
    yscale=:log10,xlabel=L"Number of vertices $(n)$", ylabel=string("Fitness function calls"),
    fillalpha=1, dpi=500, linewidth = 1.2,
     size=(320,320),outliers=false, legend= :bottomright)
    #  savefig(boxplt, "er_box.png")
    savefig(boxplt, "tg_box.png")
    boxplt
end

string.(data.AlgorithmParams_name_s)

max_stars_repo_name: DTUComputeTONIA/SDRLS

max_issues_repo_name: DTUComputeTONIA/SDRLS

max_forks_repo_name: DTUComputeTONIA/SDRLS


text:
# Copyright 2018 IBM Corp.
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#     http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.

module TestAnnotationRDF
using Test

using Serd, Serd.RDF
using Catlab

using SemanticFlowGraphs

const R = Resource

A, B, C, D = Ob(Monocl, "A", "B", "C", "D")
f = Hom("f", A, B)
g = Hom("g", B, C)
h = Hom("h", D, D)

# Object annotations

prefix = RDF.Prefix("ex", "http://www.example.org/#")
annotation = ObAnnotation(
  AnnotationID("python", "mypkg", "a"),
  Dict(
    :class => "ClassA",
    :slots => [Dict("slot" => "attrB")],
  ),
  A, [f]
)
stmts = annotation_to_rdf(annotation, prefix)
node = R("ex", "python:mypkg:a")
@test Triple(node, R("monocl", "annotatedLanguage"), Literal("python")) in stmts
@test Triple(node, R("monocl", "annotatedPackage"), Literal("mypkg")) in stmts
@test Triple(node, R("monocl", "annotatedClass"), Literal("ClassA")) in stmts
@test Triple(node, R("monocl", "codeDefinition"), R("ex","A")) in stmts

slot_node = R("ex", "python:mypkg:a:slot1")
@test Triple(node, R("monocl", "annotatedSlot"), slot_node) in stmts
@test Triple(slot_node, R("monocl", "codeDefinition"), R("ex","f")) in stmts
@test Triple(slot_node, R("monocl", "codeSlot"), Literal("attrB")) in stmts

# Morphism annotations

annotation = HomAnnotation(
  AnnotationID("python", "mypkg", "a-do-composition"),
  Dict(
    :class => ["ClassA", "MixinB"],
    :method => "do_composition",
    :inputs => [Dict("slot" => 1)],
    :outputs => [Dict("slot" => "return")],
  ),
  compose(f,g)
)
node = R("ex", "python:mypkg:a-do-composition")
root_node = R("ex", "python:mypkg:a-do-composition:diagram:root")
stmts = annotation_to_rdf(annotation, prefix)
@test Triple(node, R("monocl", "annotatedLanguage"), Literal("python")) in stmts
@test Triple(node, R("monocl", "annotatedPackage"), Literal("mypkg")) in stmts
@test Triple(node, R("monocl", "annotatedClass"), Literal("ClassA")) in stmts
@test Triple(node, R("monocl", "annotatedClass"), Literal("MixinB")) in stmts
@test Triple(node, R("monocl", "annotatedMethod"), Literal("do_composition")) in stmts
@test Triple(node, R("monocl", "codeDefinition"), root_node) in stmts

end

max_stars_repo_name: cclauss/semanticflowgraph

max_issues_repo_name: cclauss/semanticflowgraph

max_forks_repo_name: cclauss/semanticflowgraph


text:
module MinFinVol

    include("Mesh.jl")
    include("Operations.jl")

end # module


text:
"""
    wrap_show_debug(expr)

Wrap `expr` with a `@show` macro call if `DEBUG_TENSOR_DSL` exists as
an environment flag, otherwise return `expr` as-is. This is applied at
compile time, so no overhead is incurred.
"""
function wrap_show_debug(expr)
    if "DEBUG_TENSOR_DSL" in keys(ENV)
        s = :(@show)
        push!(s.args, expr)
        s
    else
        expr
    end
end

remove_line_numbers(exprs) = filter(n -> !(n isa LineNumberNode), exprs)

generate_signature(f::Function, sym, TensorType::Symbol) =
    Expr(:call, sym, f(Expr(:(::), :tensor, TensorType))...)

generate_signature(f::Function, sym, TensorType::Expr) =
    Expr(:where,
         Expr(:call, sym, f(Expr(:(::), :tensor, TensorType.args[1]))...),
         TensorType.args[2])

"""
    strip′(s::Symbol)

Strip `s` of a trailing `′`, error if it is not present.
"""
function strip′(s::Symbol)
    s = string(s)
    s[end] == '′' || throw(ArgumentError("$s does not end in `′`"))
    Symbol(s[1:prevind(s, end, 1)])
end

"""
    identify_quantum_numbers(selection_rules)

Given the block of selection rules, identify the quantum numbers
pertaining to the tensor under consideration. It is assumed that the
left-hand side consists of a primed quantum number only, which is
specified to be equal to an expression or belonging to a set/interval.
"""
function identify_quantum_numbers(selection_rules::Expr)
    if selection_rules.head == :call
        γj′ = selection_rules.args[2]
        γj′, strip′(γj′)
    elseif selection_rules.head == :block
        γj′ = Symbol[]
        for selection_rule in remove_line_numbers(selection_rules.args)
            push!(γj′, selection_rule.args[2])
        end
        Expr(:tuple, γj′...),Expr(:tuple, map(strip′,γj′)...)
    end
end

recursepm(s) = [s]

"""
    recursepm(e::Expr)

Recursively expand the expression `e` for instances of `±` and `∓`,
generating `+` and `-` minus branches.
"""
function recursepm(e::Expr)
    h = e.head
    if first(e.args) == :(±)
        a = map(recursepm, e.args[2:end])
        ee = map([:(+),:(-)]) do op
            map(Iterators.product(a...)) do a
                Expr(h, op, a...)
            end
        end
        collect(Iterators.flatten(ee))
    elseif first(e.args) == :(∓)
        a = map(recursepm, e.args[2:end])
        ee = map([:(-),:(+)]) do op
            map(Iterators.product(a...)) do a
                Expr(h, op, a...)
            end
        end
        collect(Iterators.flatten(ee))
    else
        a = map(recursepm, e.args)
        ee = map(Iterators.product(a...)) do a
            Expr(h, a...)
        end
        vec(ee)
    end
end

"""
    recursepm(lhs, rhs)

Recursively expand `lhs` and `rhs` for instances of `±` and `∓` and
join all resulting cases into a short-circuiting `||` test.
"""
function recursepm(lhs, rhs)
    ls = recursepm(lhs)
    rs = recursepm(rhs)
    cases = vec(map(((l,r),) -> :($l == $r), Iterators.product(ls, rs)))
    if length(cases) == 1
        first(cases)
    else
        a = pop!(cases)
        b = pop!(cases)
        e = Expr(:(||), a, b)
        while !isempty(cases)
            a = pop!(cases)
            e = Expr(:(||), a, e)
        end
        e
    end
end

"""
    parse_selection_rule(rule)

Parse a selection rule in the DSL of [`@tensor`](@ref); if the
selection rule is an equality, its left- and right-hand sides are
recursively expanded for possible cases of `±` and `∓` with
[`recursepm`](@ref).
"""
function parse_selection_rule(rule)
    rule = if first(rule.args) == :(==)
        recursepm(rule.args[2], rule.args[3])
    else
        rule
    end
    :($rule || return true)
end

function push_selection_rules!(f, args, selection_rules)
    for selection_rule in (selection_rules.head == :call ?
                           [selection_rules] :
                           selection_rules.args)
        if selection_rule isa Expr &&
            selection_rule.head == :call &&
            # ~ implies this particular quantum number does not
            # influence whether this matrix element vanishes or not;
            # this is typically the principal quantum number n.
            first(selection_rule.args) ≠ :(~)
            push!(args, f(selection_rule))
        elseif selection_rule isa LineNumberNode
            push!(args, selection_rule)
        end
    end
end

"""
    generate_iszero(TensorType, selection_rules)

Generate a function for testing if the matrix element of `TensorType`
vanishes, given a set a quantum number deduced from `selection_rules`.
"""
function generate_iszero(TensorType, selection_rules)
    γj′,γj = identify_quantum_numbers(selection_rules)

    signature = generate_signature(t -> (γj′, t, γj), :(Base.iszero), TensorType)
    body = Expr(:block)
    push_selection_rules!(parse_selection_rule, body.args, selection_rules)
    push!(body.args, false)
    # TODO: Should auto-generate docstring for generated iszero
    # function.
    Expr(:function, signature, body)
end

"""
    generate_rme(TensorType, selection_rules, doc, rme)

Generate a function for computing the reduced matrix element of
`TensorType`, given a set a quantum number deduced from
`selection_rules`, along with the docstring and the definition `rme`.
"""
function generate_rme(TensorType, selection_rules, doc, rme)
    γj′,γj = identify_quantum_numbers(selection_rules)
    signature = generate_signature(t -> (γj′, t, γj), :rme, TensorType)
    body = Expr(:block, :(iszero($(γj′), tensor, $(γj)) && return 0), wrap_show_debug(rme))
    fun = Expr(:function, signature, body)
    :(@doc($doc, $fun))
end

"""
    generate_copulings(TensorType, selection_rules)

Generate a function that given the quantum numbers `γj` generates
lists of all permissible `γj′` for which the reduced matrix element
`⟨γj′||::TensorType||γj⟩` does not vanish. This is deduced from the
`selection_rules`.
"""
function generate_couplings(TensorType, selection_rules)
    indeterminates = filter(r -> r.head == :call && first(r.args) == :(~),
                            remove_line_numbers(selection_rules.head == :call ?
                                                [selection_rules] :
                                                selection_rules.args))
    γj′,γj = identify_quantum_numbers(selection_rules)
    signature = generate_signature(t -> (t, γj), :couplings, TensorType)
    body = Expr(:block)
    if !isempty(indeterminates)
        push!(body.args, Expr(:call, :error, "Cannot generate all couplings for "*
                              string(TensorType)*
                              " due to the following indeterminate quantum numbers: "*
                              join(string.(indeterminates), ", ")))
    else
        push_selection_rules!(body.args, selection_rules) do rule
            op = rule.args[1]
            lhs, rhs = rule.args[2], rule.args[3]
            rhs = if op == :(∈)
                rhs
            elseif op == :(==)
                if rhs isa Symbol
                    Expr(:vect, rhs)
                else
                    throw(ArgumentError("Don't know how to handle $(rhs)"))
                end
            else
                throw(ArgumentError("Unknown operator $op"))
            end
            Expr(:(=), lhs, rhs)
        end
        push!(body.args, γj′)
    end
    γj′s = γj′ isa Symbol ? string(γj′) : join(string.(γj′.args), "")
    γjs = γj isa Symbol ? string(γj) : join(string.(γj.args), "")
    TS = string(TensorType isa Expr ? TensorType.args[1] : TensorType)
    doc = """
    $(signature)

Generate all quantum numbers `$(γj′s)` for which
`⟨$(γj′s)||::$(TS)||$(γjs)⟩` does not vanish.
"""
    fun = Expr(:function, signature, body)
    :(@doc($doc, $fun))
end

"""
    @tensor(exprs, TensorType)

Macro to generate `Base.iszero`, [`rme`](@ref) and [`couplings`](@ref)
for `TensorType`, given a set of selection rules and an expression for
the reduced matrix element.
"""
macro tensor(exprs, TensorType)
    args = remove_line_numbers(exprs.args[2].args)
    selection_rules = args[1]
    iszero = generate_iszero(TensorType, selection_rules)
    rme = generate_rme(TensorType, selection_rules, args[2], args[3])
    couplings = generate_couplings(TensorType, selection_rules)
    esc(Expr(:block, iszero, rme, couplings))
end


text:
function start_gui(;psu=PST3201(0), dmm=GDM8246(0), fgen=GFG3015(0), scope=TDS2002B(0), relays=Relays(0)) #use keyword arguments
	window, ctx = init_gui()
    # Main while loop
    try
        clear_color = Cfloat[0.45, 0.55, 0.60, 1.00]
		fps = 200
		refresh_cnt = 0
		refresh_cnt_max = 52
		base = 2
		rev_state_dict = Dict("on"=>"off", "off"=>"on")
		# instantiate gui conf objects
		psu_conf = PST3201Conf()
		dmm_conf = GDM8246Conf()
		fgen_conf = GFG3015Conf()
		scope_conf = TDS2002BConf()            #Daniel
		relays_conf = RelaysConf()            #Daniel
        while !GLFW.WindowShouldClose(window)
			refresh_cnt == refresh_cnt_max && (refresh_cnt=0)
			refresh_cnt = refresh_cnt + 1

			GLFW.PollEvents()
			# start the Dear ImGui frame
			ImGui_ImplOpenGL3_NewFrame()
			ImGui_ImplGlfw_NewFrame()
			CImGui.NewFrame()

			# Display checkboxes for activating instrument front panels
			psu_conf, dmm_conf, fgen_conf, scope_conf, relays_conf = ShowMenuWindow(psu_conf, dmm_conf, fgen_conf, scope_conf, relays_conf)		#Daniel

			# DISPLAY INSTRUMENTS FRONT PANELS
			psu_conf.active && (psu_conf = ShowPSUWindow(psu, psu_conf, rev_state_dict, refresh_cnt, base))
			dmm_conf.active && (dmm_conf = ShowDMMWindow(dmm, dmm_conf, rev_state_dict, refresh_cnt, base))
			fgen_conf.active && (fgen_conf = ShowFGENWindow(fgen, fgen_conf, rev_state_dict, refresh_cnt, base))
			scope_conf.active && (scope_conf = ShowSCOPEWindow(scope, scope_conf, rev_state_dict, refresh_cnt, base))            #Daniel
			relays_conf.active && (relays_conf = ShowRelaysWindow(relays, relays_conf, refresh_cnt, base))
			# RENDER GUI
			render_gui(window, clear_color; fps=fps)
		end
    catch e
        @error "Error in renderloop!" exception=e
        Base.show_backtrace(stderr, catch_backtrace())
    finally
        destroy_gui(ctx,window)
    end
end

max_stars_repo_name: iuliancioarca/TIVM

max_issues_repo_name: iuliancioarca/TIVM

max_forks_repo_name: iuliancioarca/TIVM


text:
#unitary test GenericStatemod
x0 = ones(6)
state0 = GenericState(x0)

io = IOBuffer()
show(io, state0)

@test scoretype(state0) == Float64
@test xtype(state0) == Array{Float64, 1}

@test isnan(state0.current_time) #Default value of start_time is void
@test isnan(state0.current_score)
x1 = [1.0]
update!(state0, x = x1) #Check the update of state0
@test state0.x == x1
@test isnan(state0.current_time) #start_time must be unchanged
@test isnan(state0.current_score)

update!(state0, current_time = 1.0)
@test state0.x == x1 #must be unchanged
@test state0.current_time == 1.0
@test isnan(state0.current_score)

reinit!(state0, x0)
@test state0.x == x0
@test isnan(state0.current_time)
@test isnan(state0.current_score)

update!(state0, x = x1)
reinit!(state0, current_time = 0.5)
@test state0.x == x1
@test state0.current_time == 0.5
@test isnan(state0.current_score)

#Test _init_field
@test _init_field(typeof(zeros(2, 2))) == zeros(0, 0)
@test _init_field(SparseMatrixCSC{Float64, Int64}) == spzeros(0, 0)
@test _init_field(typeof(zeros(2))) == zeros(0)
@test _init_field(typeof(sparse(zeros(2)))) == spzeros(0)
@test isnan(_init_field(BigFloat))
@test isnan(_init_field(typeof(1.0)))
@test isnan(_init_field(Float32))
@test isnan(_init_field(Float16))
@test _init_field(Nothing) == nothing
@test ismissing(Main.Stopping._init_field(Missing))
@test !_init_field(typeof(true))
@test _init_field(typeof(1)) == -9223372036854775808

#_check_nan_miss
@test !Stopping._check_nan_miss(nothing)
@test !Stopping._check_nan_miss(Counters())
@test !Stopping._check_nan_miss(spzeros(0))
@test !Stopping._check_nan_miss(zeros(0))
@test !Stopping._check_nan_miss(missing)
@test !Stopping._check_nan_miss(spzeros(0, 0))


text:
@testset "Operators" begin
    g3 = PathGraph(5)
    g4 = PathDiGraph(5)

    for g in testlargegraphs(g3)
        T = eltype(g)
        c = @inferred(complement(g))
        @test nv(c) == 5
        @test ne(c) == 6

        gb = @inferred(blkdiag(g, g))
        @test nv(gb) == 10
        @test ne(gb) == 8

        hp = PathGraph(2)
        h = Graph{T}(hp)
        @test @inferred(intersect(g, h)) == h

        hp = PathGraph(4)
        h = Graph{T}(hp)

        z = @inferred(difference(g, h))
        @test nv(z) == 5
        @test ne(z) == 1
        z = @inferred(difference(h, g))
        @test nv(z) == 4
        @test ne(z) == 0
        z = @inferred(symmetric_difference(h, g))
        @test z == symmetric_difference(g, h)
        @test nv(z) == 5
        @test ne(z) == 1

        h = Graph{T}(6)
        add_edge!(h, 5, 6)
        e = SimpleEdge(5, 6)

        z = @inferred(union(g, h))
        @test has_edge(z, e)
        @test z == PathGraph(6)

        # Check merge_vertices function.
        h = Graph{T}(7)
        add_edge!(h, 2, 5)
        add_edge!(h, 3, 6)
        add_edge!(h, 1, 7)
        add_edge!(h, 6, 5)

        vs = [2, 3, 7, 3, 3, 2]
        hmerged = merge_vertices(h, vs)
        @test neighbors(hmerged, 1) == [2]
        @test neighbors(hmerged, 2) == [1, 4, 5]
        @test neighbors(hmerged, 3) == []
        @test neighbors(hmerged, 4) == [2, 5]

        new_map = @inferred(merge_vertices!(h, vs))
        @test new_map == [1, 2, 2, 3, 4, 5, 2]
        @test neighbors(h, 1) == [2]
        @test neighbors(h, 2) == [1, 4, 5]
        @test neighbors(h, 3) == []
        @test neighbors(hmerged, 4) == [2, 5]
        @test hmerged == h

        h = Graph{T}(7)
        add_edge!(h, 1, 2)
        add_edge!(h, 2, 3)
        add_edge!(h, 2, 4)
        add_edge!(h, 3, 4)
        add_edge!(h, 3, 7)
        new_map = @inferred(merge_vertices!(h, [2, 3, 2, 2]))
        @test new_map == [1, 2, 2, 3, 4, 5, 6]
        @test neighbors(h, 2) == [1, 3, 6]
        @test neighbors(h, 1) == [2]
        @test neighbors(h, 3) == [2]
        @test neighbors(h, 4) == Int[]
        @test neighbors(h, 6) == [2]
        @test ne(h) == 3
        @test nv(h) == 6

        h2 = Graph{T}(7)
        add_edge!(h2, 1, 2)
        add_edge!(h2, 2, 3)
        add_edge!(h2, 2, 4)
        add_edge!(h2, 3, 4)
        add_edge!(h2, 3, 7)
        add_edge!(h2, 6, 7)
        new_map = @inferred(merge_vertices!(h2, [2, 7, 3, 2]))
        @test new_map == [1, 2, 2, 3, 4, 5, 2]
        @test neighbors(h2, 2) == [1, 3, 5]
        @test neighbors(h2, 1) == [2]
        @test neighbors(h2, 3) == [2]
        @test neighbors(h2, 4) == Int[]
        @test neighbors(h2, 5) == [2]
        @test ne(h2) == 3
        @test nv(h2) == 5

    end
    for g in testlargedigraphs(g4)
        T = eltype(g)
        c = @inferred(complement(g))
        @test nv(c) == 5
        @test ne(c) == 16

        h = DiGraph{T}(6)
        add_edge!(h, 5, 6)
        e = SimpleEdge(5, 6)

        z = @inferred(union(g, h))
        @test has_edge(z, e)
        @test z == PathDiGraph(6)

        # Check merge_vertices function.
        # h = DiGraph{T}(7)
        # add_edge!(h, 2, 5)
        # add_edge!(h, 3, 6)
        # add_edge!(h, 1, 7)
        # add_edge!(h, 4, 3)
        # add_edge!(h, 6, 5)
        # new_map = @inferred(merge_vertices!(h, [2, 3, 7, 3, 3, 2]))
        # @test new_map == [1, 2, 2, 4, 5, 3, 2]
        # @test in_neighbors(h, 2) == [1, 4]
        # @test out_neighbors(h, 2) == [3, 5]
        # @test out_neighbors(h, 3) == [5]

    end

    re1 = Edge(2, 1)
    gr = @inferred(reverse(g4))
    for g in testdigraphs(gr)
        T = eltype(g)
        @test re1 in edges(g)
        @inferred(reverse!(g))
        @test g == DiGraph{T}(g4)
    end

    gx = CompleteGraph(2)

    for g in testlargegraphs(gx)
        T = eltype(g)
        hc = CompleteGraph(2)
        h = Graph{T}(hc)
        z = @inferred(blkdiag(g, h))
        @test nv(z) == nv(g) + nv(h)
        @test ne(z) == ne(g) + ne(h)
        @test has_edge(z, 1, 2)
        @test has_edge(z, 3, 4)
        @test !has_edge(z, 1, 3)
        @test !has_edge(z, 1, 4)
        @test !has_edge(z, 2, 3)
        @test !has_edge(z, 2, 4)
    end

    gx = Graph(2)
    for g in testgraphs(gx)
        T = eltype(g)
        h = Graph{T}(2)
        z = @inferred(join(g, h))
        @test nv(z) == nv(g) + nv(h)
        @test ne(z) == 4
        @test !has_edge(z, 1, 2)
        @test !has_edge(z, 3, 4)
        @test has_edge(z, 1, 3)
        @test has_edge(z, 1, 4)
        @test has_edge(z, 2, 3)
        @test has_edge(z, 2, 4)
    end

    px = PathGraph(10)
    for p in testgraphs(px)
        x = @inferred(p * ones(10))
        @test  x[1] == 1.0 && all(x[2:(end - 1)] .== 2.0) && x[end] == 1.0
        @test size(p) == (10, 10)
        @test size(p, 1) == size(p, 2) == 10
        @test size(p, 3) == 1
        @test sum(p, 1) == sum(p, 2)
        @test_throws ErrorException sum(p, 3)
        @test sparse(p) == adjacency_matrix(p)
        @test length(p) == 100
        @test ndims(p) == 2
        @test issymmetric(p)
    end

    gx = DiGraph(4)
    add_edge!(gx, 1, 2); add_edge!(gx, 2, 3); add_edge!(gx, 1, 3); add_edge!(gx, 3, 4)
    for g in testdigraphs(gx)
        @test @inferred(g * ones(nv(g))) == [2.0, 1.0, 1.0, 0.0]
        @test sum(g, 1) ==  [0, 1, 2, 1]
        @test sum(g, 2) ==  [2, 1, 1, 0]
        @test sum(g) == 4
        @test @inferred(!issymmetric(g))
    end

    nx = 20; ny = 21
    gx = PathGraph(ny)
    for g in testlargegraphs(gx)
        T = eltype(g)
        hp = PathGraph(nx)
        h = Graph{T}(hp)
        c = @inferred(cartesian_product(g, h))
        gz = @inferred(crosspath(ny, PathGraph(nx)))
        @test gz == c
    end
    function crosspath_slow(len, h)
        g = h
        m = nv(h)
        for i in 1:(len - 1)
            k = nv(g)
            g = blkdiag(g, h)
            for v in 1:m
                add_edge!(g, v + (k - m), v + k)
            end
        end
        return g
    end

    gx = CompleteGraph(2)
    for g in testgraphs(gx)
        h = @inferred(cartesian_product(g, g))
        @test nv(h) == 4
        @test ne(h) == 4

        h = @inferred(tensor_product(g, g))
        @test nv(h) == 4
        @test ne(h) == 1
    end
    g2 = CompleteGraph(2)
    for g in testgraphs(g2)
        @test crosspath_slow(2, g) == crosspath(2, g)
    end

    ## test subgraphs ##

    gb = smallgraph(:bull)
    for g in testgraphs(gb)
        n = 3
        h = @inferred(g[1:n])
        @test nv(h) == n
        @test ne(h) == 3

        h = @inferred(g[[1, 2, 4]])
        @test nv(h) == n
        @test ne(h) == 2

        h = @inferred(g[[1, 5]])
        @test nv(h) == 2
        @test ne(h) == 0
        @test typeof(h) == typeof(g)
    end

    gx = DiGraph(100, 200)
    for g in testdigraphs(gx)
        h = @inferred(g[5:26])
        @test nv(h) == 22
        @test typeof(h) == typeof(g)
        @test_throws ErrorException g[[1, 1]]

        r = 5:26
        h2, vm = @inferred(induced_subgraph(g, r))
        @test h2 == h
        @test vm == collect(r)
        @test h2 == g[r]
    end

    g10 = CompleteGraph(10)
    for g in testgraphs(g10)
        sg, vm = @inferred(induced_subgraph(g, 5:8))
        @test nv(sg) == 4
        @test ne(sg) == 6

        sg2, vm = @inferred(induced_subgraph(g, [5, 6, 7, 8]))
        @test sg2 == sg
        @test vm[4] == 8

        elist = [
          SimpleEdge(1, 2), SimpleEdge(2, 3), SimpleEdge(3, 4),
          SimpleEdge(4, 5), SimpleEdge(5, 1)
        ]
        sg, vm = @inferred(induced_subgraph(g, elist))
        @test sg == CycleGraph(5)
        @test sort(vm) == [1:5;]
    end

    gs = StarGraph(10)
    for g in testgraphs(gs)
        T = eltype(g)
        @test @inferred(egonet(g, 1, 0)) == Graph{T}(1)
        @test @inferred(egonet(g, 1, 1)) == g
        @test @inferred(ndims(g)) == 2
    end
end


text:
"""
Julia expression wrapper
"""
struct Expression
    expr::Expr
    syms::Vector{Symbol}
end
Expression(ex::Expr) = Expression(ex, sort!(symbols(ex)))
show(io::IO, e::Expression) = infix(io, e.expr)
(e::Expression)(val) = evaluate(val, e.expr, e.syms)

function symbols(ex::Expr)
    syms = Symbol[]
    for e in ex.args
        isa(e, Symbol) && push!(syms, e)
        isa(e, Expr) && append!(syms, symbols(e))
    end
    unique!(syms)
end

height(ex) = isa(ex, Expr) ? maximum( height(e) for e in ex.args )+1 : 0
nodes(ex) = !isa(ex, Expr) ? 1 : length(ex.args) > 0 ? sum( nodes(e) for e in ex.args ) : 0
length(ex) = nodes(ex)

function depth(root, ex; d=0)
    return if root == ex
        d
    elseif isa(root, Expr)
        maximum( depth(e, ex, d=d+1) for e in root.args )
    else
        -1
    end
end

function copyto!(ex1::Expr, ex2::Expr)
    ex1.head = ex2.head
    ex1.args = deepcopy(ex2.args)
    ex1
end

function randsubexpr(ex)
    !isa(ex, Expr) && return ex
    csize = map(nodes, ex.args[2:end])
    cidx = rand(1:sum(csize)+1)-1
    if cidx == 0
        ex
    else
        nidx = findfirst(i->cidx<=i, cumsum(csize))
        randsubexpr(ex.args[nidx+1])
    end
end

function getindex(ex::Expr, idx::Int)
    if idx == 0
        ex
    else
        csize = cumsum(map(nodes, ex.args[2:end]))
        nidx = findfirst(i->i>=idx, csize)
        nex = ex.args[nidx+1]
        !isa(nex, Expr) ? nex : nex[csize[nidx]-idx]
    end
end

function setindex!(ex::Expr, subex, idx::Int)
    csize = cumsum(map(nodes, ex.args[2:end]))
    nidx = findfirst(i->i>=idx, csize)
    nex = ex.args[nidx+1]
    if !isa(nex, Expr) || csize[nidx]-idx == 0
        ex.args[nidx+1] = subex
    else
        nex[csize[nidx]-idx] = subex
    end
end

function evaluate(val, ex::Expr, psyms::Vector{Symbol})
    exprm = ex.args
    exvals = (isa(nex, Expr) || isa(nex, Symbol) ? evaluate(val, nex, psyms) : nex for nex in exprm[2:end])
    exprm[1](exvals...)
end

function evaluate(val, ex::Symbol, psyms::Vector{Symbol})
    pidx = findfirst(isequal(ex), psyms)
    val[pidx]
end

iszeronum(root) = isa(root, Number) && iszero(root)
isonenum(root) = isa(root, Number) && isone(root)

"""
    simplify!(expr)

Simplify an algebraic expression
"""
function simplify!(root)
    # return if not an expression
    !isa(root, Expr) && return root
    # return if not a function
    root.head != :call && return root

    # recursively simplify arguments
    for i in 2:length(root.args)
        !isa(root.args[i], Expr) && continue
        root.args[i] = simplify!(root.args[i])
    end

    # some elementary simplification rules
    if root.args[1] == (-) && root.args[2] == root.args[3] # additive inverse: x-x = 0
        root = 0
    elseif root.args[1] == (/) && (root.args[2] == root.args[3] || iszeronum(root.args[3])) # look for multiplicative inverse x/x = 1
        root = 1
    elseif (root.args[1] == (*) || root.args[1] == (/)) &&  (iszeronum(root.args[2]) || iszeronum(root.args[3])) # look for 0: x*0 = 0*x = 0
        root = 0
    elseif (root.args[1] == (+) || root.args[1] == (-)) && iszeronum(root.args[3]) # look for 0: x ± 0 = 0
        root = root.args[2]
    elseif root.args[1] == (+) && iszeronum(root.args[2]) # look for 0: 0 ± x = 0
        root = root.args[3]
    elseif (root.args[1] == (*) || root.args[1] == (/)) && isonenum(root.args[3]) # x*1 = x || x/1 = x
        root = root.args[2]
    elseif root.args[1] == (*) && isonenum(root.args[2]) # 1*x = x
        root = root.args[3]
    else
        # x+x = 2x
        if root.args[1] == (+) && root.args[2] == root.args[3]
            root.args[1] = (*)
            root.args[2] = 2
        end
        # evaluate numerical arithmetic
        if all(e->!(isa(e, Symbol) || isa(e, Expr)), root.args[2:end])
            root = root.args[1](root.args[2:end]...)
        end
    end
    return root
end

function infix(io::IO, root; digits=3)
    if isa(root, Number)
        print(io, round(root, digits=digits))
    elseif isa(root, Expr) && root.head == :call
        if root.args[1] ∈ [+, -, *, /, ^]
            print(io, "(")
            infix(io, root.args[2])
            infix(io, root.args[1])
            infix(io, root.args[3])
            print(io, ")")
        else
            infix(io, root.args[1])
            print(io, "(")
            for ex in root.args[2:end]
                infix(io, ex)
            end
            print(io, ")")
        end
    else
        print(io, root)
    end
end
show(io::IO, ::MIME"text/html", e::Expression) = infix(io, e.expr)

latex(io::IO, root::Number; digits=3) = print(io, round(root, digits=digits))
latex(io::IO, root::Symbol; kwargs...) = print(io, root)
function latex(io::IO, root::Expr; digits=3)
    root.head != :call && print(io, expr)
    if root.args[1] == (/)
        print(io, "\\frac{")
        latex(io, root.args[2], digits=digits)
        print(io, "}{")
        latex(io, root.args[3], digits=digits)
        print(io, "}")
    elseif root.args[1] ∈ [+, -, *]
        print(io, "\\left(")
        latex(io, root.args[2], digits=digits)
        print(io, root.args[1])
        latex(io, root.args[3], digits=digits)
        print(io, "\\right)")
    elseif root.args[1] == (^)
        print(io, "{")
        latex(io, root.args[2], digits=digits)
        print(io, "}^{")
        latex(io, root.args[3], digits=digits)
        print(io, "}")
    else
        print(io, "\\")
        print(io, root.args[1])
        print(io, "\\left(")
        for ex in root.args[2:end]
            latex(io, ex, digits=digits)
        end
        print(io, "\\right)")
    end
end
show(io::IO, ::MIME"text/latex", e::Expression) = latex(io, e.expr)


text:
using Flux

#-----------
# Tools
#-----------

"""
    flatten_model(model)

Flattens the model weights into a 1D array.
"""
flatten_model(model::Chain) = reduce(vcat, [layer[:] for layer in Flux.params(model)])

"""
    load_flatten_model!(model, flatten_model)

Reshapes the `flatten_model` and load it into `model`.
"""
function load_flatten_model!(model::Chain, flatten_model::Vector{Float32})
    weights = [similar(layer) for layer in Flux.params(model)]

    start = 1
    for i in 1:length(weights)
        shape = size(weights[i])
        stop = reduce(*, shape)
        weights[i] = reshape(flatten_model[start:(start - 1 + stop)], shape)
        start += stop
    end

    Flux.loadparams!(model, weights)
end

# ---------
# MLP
# ---------

"""
    initialize_mlp(input, hidden, output)

Create a 2-hidden layers MLP of hidden params each.
"""
function initialize_mlp(input::Int, hidden::Int, output::Int)
    return Chain(
        Dense(input, hidden, relu),
        Dense(hidden, hidden, relu),
        Dense(hidden, output),
    )
end

# ---------
# CNN
# ---------

"""
    initialize_cnn(imgsize = (28,28,1), numclasses=10)

Create a 2-conv layers CNN to work with the MNIST dataset.
"""
function initialize_cnn(imgsize = (28,28,1), numclasses::Int=10)
    return Chain(
        # 5x5x1x32 + 32 = 832
        Conv((5, 5), imgsize[3]=>32, pad=(1,1), relu),
        MaxPool((2,2)),

        # 5x5x32x64 + 64 = 51264
        Conv((5, 5), 32=>64, pad=(1,1), relu),
        MaxPool((2,2)),

        flatten,

        # 5*5*64 * 512 + 512 = 819712
        Dense(1600, 512, relu),

        # 512 * 10 + 10 = 5130
        Dense(512, numclasses),
    )
end


text:
"""
  meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T<:Real

Replicates matlab function for generating structured meshes
"""
function meshgrid(vx::AbstractVector{T}, vy::AbstractVector{T}) where T<:Real
	m, n = length(vy), length(vx)
  vx = reshape(vx, 1, n)
  vy = reshape(vy, m, 1)

  return (repeat(vx, m, 1)', repeat(vy, 1, n)')
end

# square [x0, x1, y0, y1]
function squareMesh(square::Vector{T},N::Int,order::Int) where T<:Real
  x = range(square[1],stop=square[2],length=N)
  y = range(square[3],stop=square[4],length=N)

  # create nodes
  X,Y = meshgrid(x,y)
  xy = [X[:] Y[:]]

  # create element connectivity
  ni = size(X,1); # number of rows
  nNodes = size(xy,1);
  t2nidxMap = 1:nNodes-ni;
  topNode = ni:ni:nNodes-ni;
  t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode))
  k = t2nidxMap;
  cm = [k k.+1 k.+ni.+1 k.+ni]

  if order==2
    xy,cm = Order1ToOrder2(xy,cm)
    xy,cm = deleteDuplicates(xy,cm)
  end

  # create boundary nodes
  left=Int[]; top=Int[]; right=Int[]; bottom=Int[]
  for i=1:size(xy,1)
    # left
    if isMachZero(xy[i,1],square[1])
      push!(left,i)
    end
    # right
    if isMachZero(xy[i,1],square[2])
      push!(right,i)
    end
    # bottom
    if isMachZero(xy[i,2],square[3])
      push!(bottom,i)
    end
    # top
    if isMachZero(xy[i,2],square[4])
      push!(top,i)
    end
  end

  # generate bdry object
  bndNames = Index([:left,:top,:right,:bottom])
  bdry = DictList{Int}([left,top,right,bottom],bndNames)

  (order==1 ? o=:Linear : o=:Quadratic)

  # Generate Mesh!
  return ScalarMesh(ArrayToNodes(xy),ArrayToElements(cm),bdry,o)
end

function squareMeshFluid(square,N)
  mesh = squareMesh(square,N,2)

  xy = NodesToArray(mesh.xy)
  cm = ElementsToArray(mesh.cm)

  xyptemp,cmptemp = pMesh(xy,cm)

  xyp = ArrayToNodes(xyptemp)
  cmp = ArrayToElements(cmptemp)

  return FluidMesh(mesh.xy,mesh.cm,mesh.bdry,:Quadratic,xyp,cmp)
end

function axisymMesh(Nr,Nt,Rmax,order;phi=(x -> x))
  xi = linspace(0,1,Nr)
  t = linspace(0,pi,Nt)

  # create nodes
  XI,THETA = meshgrid(xi,t)
  xy = [XI[:] THETA[:]]

  # create element connectivity
  ni = size(XI,1); # number of rows
  nNodes = size(xy,1)
  t2nidxMap = 1:nNodes-ni;
  topNode = ni:ni:nNodes-ni;
  t2nidxMap = deleteat!(collect(t2nidxMap),collect(topNode))
  k = t2nidxMap;
  cm = [k k+1 k+ni+1 k+ni]

  if order==2
    xy,cm = Order1ToOrder2(xy,cm)
    xy,cm = deleteDuplicates(xy,cm)
  end

  # create boundary nodes
  inflow=Int[]; outflow=Int[]; axis=Int[]; sphere=Int[]
  for i=1:size(xy,1)
    # axis
    if isMachZero(xy[i,2],0.0) || isMachZero(xy[i,2],Float64(pi))
      push!(axis,i)
    end
    # inflow
    if isMachZero(xy[i,1],1) && xy[i,2]<=Float64(pi)/2
      push!(inflow,i)
    end
    # outflow
    if isMachZero(xy[i,1],1) && xy[i,2]>Float64(pi)/2
      push!(outflow,i)
    end
    # sphere
    if isMachZero(xy[i,1],0.0)
      push!(sphere,i)
    end
  end

  # transform xy nodes
  transform(r) = exp(log(Rmax)*phi(r))

  XI2    = xy[:,1]
  THETA2 = xy[:,2]

  R = transform.(XI2)
  X = R.*cos.(THETA2)
  Y = R.*sin.(THETA2)
  xy = [X[:] Y[:]]

  # generate bdry object
  bndNames = Index([:inflow,:outflow,:axis,:sphere])
  bdry = DictList{Int}([inflow,outflow,axis,sphere],bndNames)

  (order==1 ? o=:Linear : o=:Quadratic)

  # Generate Mesh!
  return ScalarMesh(ArrayToNodes(xy),ArrayToElements(cm),bdry,o)
end

function axisymFluid(Nr,Nt,Rmax;phi=(x -> x))
  mesh = axisymMesh(Nr,Nt,Rmax,2;phi=phi)

  xy = NodesToArray(mesh.xy)
  cm = ElementsToArray(mesh.cm)

  xyptemp,cmptemp = pMesh(xy,cm)

  xyp = ArrayToNodes(xyptemp)
  cmp = ArrayToElements(cmptemp)

  return FluidMesh(mesh.xy,mesh.cm,mesh.bdry,:Quadratic,xyp,cmp)
end

# function to convert array to Nodes
function ArrayToNodes(xy)
  return [Node(xy[i,1],xy[i,2]) for i=1:size(xy,1)]
end

# function to convert Nodes to array form
function NodesToArray(xy)
    # get xy into correct form
  x = [xy[i].x for i=1:length(xy)]
  y = [xy[i].y for i=1:length(xy)]
  return hcat(x,y)
end

# function to convert array to elements
function ArrayToElements(cm)
  return [Element(cm[i,:]) for i in 1:size(cm,1)]
end

# function to convert Elements to array form
function ElementsToArray(cm)
  return [cm[i].NodeList[j] for i=1:length(cm), j=1:length(cm[1].NodeList)]
end

# converts order 1 meshes to order 2 meshes
function Order1ToOrder2(xy,cm)
  nNodesOld = size(xy,1)
  nElm      = size(cm,1)

  xNew = Float64[]
  yNew = Float64[]
  cmNew = Array{Int}(undef,nElm,9)

  cmNew[:,1:4] = cm   # will throw error if cm isn't quad

  # for each element
  newNodeInd = nNodesOld + 1
  for el=1:nElm
    # for nodes 5 - 9
    for i=5:9
      x=0.0; y=0.0
      if i==5
        x = (xy[cm[el,1],1] + xy[cm[el,2],1])/2
        y = (xy[cm[el,1],2] + xy[cm[el,2],2])/2
      elseif i==6
        x = (xy[cm[el,2],1] + xy[cm[el,3],1])/2
        y = (xy[cm[el,2],2] + xy[cm[el,3],2])/2
      elseif i==7
        x = (xy[cm[el,3],1] + xy[cm[el,4],1])/2
        y = (xy[cm[el,3],2] + xy[cm[el,4],2])/2
      elseif i==8
        x = (xy[cm[el,4],1] + xy[cm[el,1],1])/2
        y = (xy[cm[el,4],2] + xy[cm[el,1],2])/2
      elseif i==9
        x = (xy[cm[el,1],1] + xy[cm[el,2],1] +
             xy[cm[el,3],1] + xy[cm[el,4],1])/4
        y = (xy[cm[el,1],2] + xy[cm[el,2],2] +
             xy[cm[el,3],2] + xy[cm[el,4],2])/4
      end
      # add node to xy
      push!(xNew,x)
      push!(yNew,y)

      # add node index to cm
      cmNew[el,i] = newNodeInd
      newNodeInd += 1
    end
  end

  xyNew = Array{Float64}(undef,newNodeInd-1,2)

  for i=1:(newNodeInd-1)
    if i <= nNodesOld
      xyNew[i,1] = xy[i,1]; xyNew[i,2] = xy[i,2]
    else
      k = i-nNodesOld
      xyNew[i,1] = xNew[k]; xyNew[i,2] = yNew[k]
    end
  end

  return xyNew, cmNew
end

function deleteDuplicates!(mesh)
  # transition data to arrays
  xyold = NodesToArray(mesh.xy)
  cmold = ElementsToArray(mesh.cm)

  # delete duplicates and reassign cm
  xynew,cmnew = deleteDuplicates(xyold,cmold)

  # create new mesh out of this

end

# need to make function that replaces/deletes duplicate nodes
# for the love of all that is holy, do not change anything in here. the indices
# are a nightmare
function deleteDuplicates(xy,cm)
  nNodesOld = size(xy,1)
  nElm   = size(cm,1)

  nodeIndPos5 = view(cm,:,5)
  nodeIndPos6 = view(cm,:,6)
  nodeIndPos7 = view(cm,:,7)
  nodeIndPos8 = view(cm,:,8)

  xyNew = Float64[]
  skipList = Int[]
  cmNew = copy(cm)

  # for each node in 5 or 7 position
  for ti=1:length(nodeIndPos5)
    i = nodeIndPos5[ti]
    if !(i in skipList)
      for tj=1:length(nodeIndPos7)
        j = nodeIndPos7[tj]
        # check whether x AND y match within epsilon
        if isMachZero(xy[i,1],xy[j,1]) &&
           isMachZero(xy[i,2],xy[j,2]) &&
           i!=j
          # find other index for node j
          kjfind = something(findfirst(isequal(j), cm), 0)
          # 1) assign corresponding node index to cm of other node
          cmNew[kjfind] = i
          # 3) add other node to list so that it isn't crossed later on
          push!(skipList,j)
          break
        end
      end
    end
  end

  # for each node in 6 or 8 position
  for ti=1:length(nodeIndPos6)
    i = nodeIndPos6[ti]
    if !(i in skipList)
      for tj=1:length(nodeIndPos8)
        j = nodeIndPos8[tj]
        # check whether x AND y match within epsilon
        if isMachZero(xy[i,1],xy[j,1]) &&
           isMachZero(xy[i,2],xy[j,2]) &&
           i!=j
          # find other index for node j
          kjfind = something(findfirst(isequal(j), cmNew), 0)
          # 1) assign corresponding node index to cm of other node
          cmNew[kjfind] = i
          # 3) add other node to list so that it isn't crossed later on
          push!(skipList,j)
          break
        end
      end
    end
  end

  cmsubtract = zeros(Int,(length(cm)))
  # forward sweep to decide how many to subtract from each cm
  tempind=1
  for i in cmNew
    for k in skipList
      if i>=k
        cmsubtract[tempind]+=1
      end
    end
    tempind+=1
  end

  # backward sweep to actually take off those deleted indices
  for i=1:length(cm)
    cmNew[i] -= cmsubtract[i]
  end

  xyNew = xy[setdiff(1:end, skipList), :]

  return xyNew, cmNew
end

# function to check whether two numbers are within machine epsilon
function isMachZero(x,y)
  ep = eps()
  if abs(x-y) < ep
    return true
  else
    return false
  end
end

"""
  pMesh(xy,cm) -> xyp, cmp

Generates pressure mesh equivalents given a velocity mesh (xy,cm).

Might need to generate map between velocity and pressure nodes ....?
"""
function pMesh(xy,cm)
  Nel         = size(cm,1)
  NumPNodes   = length(unique(@view cm[:,1:4]))
  xyp         = zeros(NumPNodes,2)
  cmp         = zeros(Int,(Nel,4))

  nodeCheckarr::Array{Int,1} = []

  elCounter = 1
  nodeCounter = 1
  nodeCheck = 1
  for el=1:Nel
    for i=1:4
      c = cm[el,i]
      if c in nodeCheckarr
        cmp[elCounter,i] = something(findfirst(isequal(c),nodeCheckarr),0)
      else
        push!(nodeCheckarr,c)
        xyp[nodeCounter,:] = @view xy[c,:]
        cmp[elCounter,i] = nodeCounter
        nodeCounter += 1
      end

      nodeCheck += 1
    end

    elCounter += 1
  end

  return xyp,cmp
end

"""
  lowerUV(xy,cm,u,v) -> u,v

computes velocity on coarser mesh using p-mesh
"""
function lowerUV(xy,cm,u,v)
  Nel       = size(cm,1)
  NumPNodes = length(unique(cm[:,1:4]))
  q1u   = zeros(NumPNodes,1)
  q1v   = zeros(NumPNodes,1)
  cmp = 0*Array{Int}(undef,Nel,4)
  nodeCheckarr::Array{Int,1} = []

  elCounter = 1
  nodeCounter = 1
  nodeCheck = 1
  for el=1:Nel
    for i=1:4
      c = cm[el,i]
      if c in nodeCheckarr
        cmp[elCounter,i] = something(findfirst(isequal(c),nodeCheckarr),0)
      else
        push!(nodeCheckarr,c)
        q1u[nodeCounter] = u[c]
        q1v[nodeCounter] = v[c]
        nodeCounter += 1
      end

      nodeCheck += 1
    end

    elCounter += 1
  end

  return q1u,q1v
end

"""
  lowerScalar(xy,cm,q2s) -> q1s

computes scalar on coarser mesh using p-mesh
"""
function lowerScalar(xy,cm,q2s)
  Nel       = size(cm,1)
  NumQ1Nodes = length(unique(cm[:,1:4]))
  q1s   = zeros(NumQ1Nodes,1)
  #cmp = 0*Array{Int}(Nel,4)
  nodeCheckarr::Array{Int,1} = []

  elCounter = 1
  nodeCounter = 1
  nodeCheck = 1
  for el=1:Nel
    for i=1:4
      c = cm[el,i]
      if c in nodeCheckarr
        #cmp[elCounter,i] = findfirst(nodeCheckarr,c)
      else
        push!(nodeCheckarr,c)
        q1s[nodeCounter] = q2s[c]
        nodeCounter += 1
      end

      nodeCheck += 1
    end

    elCounter += 1
  end

  return q1s
end


text:
# activates the Julia environment associated with the working folder
using Pkg
Pkg.activate(".")


text:
# UnitCommitment.jl: Optimization Package for Security-Constrained Unit Commitment
# Copyright (C) 2020, UChicago Argonne, LLC. All rights reserved.
# Released under the modified BSD license. See COPYING.md for more details.

using UnitCommitment, Cbc, JuMP

@testset "Initial conditions" begin
    # Load instance
    instance = UnitCommitment.read("$(pwd())/fixtures/case118-initcond.json.gz")
    optimizer = optimizer_with_attributes(Cbc.Optimizer, "logLevel" => 0)

    # All units should have unknown initial conditions
    for g in instance.units
        @test g.initial_power === nothing
        @test g.initial_status === nothing
    end

    # Generate initial conditions
    UnitCommitment.generate_initial_conditions!(instance, optimizer)

    # All units should now have known initial conditions
    for g in instance.units
        @test g.initial_power !== nothing
        @test g.initial_status !== nothing
    end

    # TODO: Check that initial conditions are feasible
end