text:
###############################################################################
## Least squares problem parametrized by regularizer reg
#    min_x  0.5 * ||Ax-y||² + g(reg, x)
struct LeastsquaresPb{Tr, Tm, Txman} <: CompositeProblem
    A::Matrix{Float64}
    y::Vector{Float64}
    regularizer::Tr
    n::Int64
    x0::Txman
    M_x0::Tm
end

problem_dimension(pb::LeastsquaresPb) = pb.n

## f
# 0th order
f(pb::LeastsquaresPb, x) = 0.5 * norm(pb.A * x - pb.y)^2
# f(pb::LeastsquaresPb, x::AbstractMatrix) = f(pb, view(x, :))

# 1st order
∇f!(pb::LeastsquaresPb, res, x) = (res .= transpose(pb.A) * (pb.A * x - pb.y))
∇f(pb::LeastsquaresPb, x) = transpose(pb.A) * (pb.A * x - pb.y)

# # ∇f!(pb::LeastsquaresPb, res, x::AbstractMatrix) = ∇f!(pb, view(res, :), view(x, :))
# ∇f(pb::LeastsquaresPb, x::AbstractMatrix) = reshape(transpose(pb.A) * (pb.A * view(x, :) - pb.y), size(x))

# 2nd order
∇²f_h!(pb::LeastsquaresPb, res, x, h) = (res .= transpose(pb.A) * (pb.A * h))
∇²f_h(pb::LeastsquaresPb, x, h) = transpose(pb.A) * (pb.A * h)

# conditioning
get_gradlips(pb::LeastsquaresPb) = opnorm(pb.A)^2
get_μ_cvx(pb::LeastsquaresPb) = (svdvals(pb.A)[end])^2

## instance generation:
function build_original_signal(::regularizer_l1, n, sparsity, seed)
    Random.seed!(seed)
    x0 = rand(Normal(), n)

    Random.seed!(seed)
    inds_nz = Set(randperm(n)[1:sparsity])

    nnz_coords = BitVector(zeros(n))
    for ind in inds_nz
        nnz_coords[ind] = 1
    end

    x0 .*= nnz_coords

    return x0, l1Manifold(nnz_coords)
end

function build_original_signal(reg::regularizer_distball, n, sparsity, seed)
    Random.seed!(seed)
    x0 = rand(Normal(), n)
    x0 *= reg.r / norm(x0, reg.p)

    return x0, PSphere(reg.p, reg.r, n)
end

function build_original_signal(reg::regularizer_lnuclear, mn, sparsity, seed)
    m, n = mn

    Random.seed!(seed)
    x0 = rand(Normal(), m, n)

    F = svd(x0)
    x0 = F.U[:, 1:sparsity] * Diagonal(F.S[1:sparsity]) * F.Vt[1:sparsity, :]

    return vec(x0), FixedRankMatrices(m, n, sparsity)
end

function build_original_signal(reg::regularizer_group{Tr}, n, sparsity, seed) where {Tr}
    x0 = zeros(n)
    Ms = []
    for (i, group) in enumerate(reg.groups)
        x0i, M0 = build_original_signal(reg.regs[i], length(group), sparsity, seed)
        x0[group] .= x0i
        push!(Ms, M0)
    end

    return x0, ProductManifold(Ms...)
end

function update_regularizationstrength!(regularizer, delta)
    regularizer.λ = delta
    return
end
function update_regularizationstrength!(
    regularizer::regularizer_group{Tr},
    delta,
) where {Tr}
    for i in 1:length(regularizer.regs)
        update_regularizationstrength!(regularizer.regs[i], delta)
    end
    return
end

"""
    get_random_qualifiedleastsquares(n, m, regularizer, sparsity; seed=1234)

Build a least-squares problem that fits the signal recovery framework. The optimization
problem solution should be close to the original signal `x0` and lie exactly on the same
optimal manifold `M_x0`.
"""
function get_random_qualifiedleastsquares(n, m, regularizer, sparsity; seed = 1234, A=nothing, delta = 0.01)
    # A is drawn from the standard normal distribution
    # b = Ax0 + e where e is taken from the normal distribution with standard deviation 0.001.
    # We set λ1 so that the original sparsity is ultimately recovered.
    n_vectorialpb = isa(n, Tuple) ? n[1]*n[2] : n

    Random.seed!(seed)
    if isnothing(A)
        A = rand(Normal(), m, n_vectorialpb)
    end

    x0, M0 = build_original_signal(regularizer, n, sparsity, seed)

    Random.seed!(seed + 3)
    e = rand(Normal(0, delta^2), m)

    y = A * x0 + e

    update_regularizationstrength!(regularizer, delta)
    return LeastsquaresPb(A, y, regularizer, n_vectorialpb, x0, M0)
end

function get_random_qualifiedlasso(n, m, sparsity; seed = 1234)
    return get_random_qualifiedleastsquares(
        n,
        m,
        regularizer_l1(1.0),
        sparsity;
        seed = seed,
    )
end

# function get_randomlasso(n, m, sparsity; reg=regularizer_l1, seed=1234)
#     # A is drawn from the standard normal distribution
#     # b = Ax0 + e where e is taken from the normal distribution with standard deviation 0.001.
#     # We set λ1 so that the original sparsity is ultimately recovered.

#     elseif reg <: regularizer_l12
#         T = reg.parameters[1]
#         ngroups = Int(ceil(n/T))

#         Random.seed!(seed)
#         inds_nz = Set(randperm(ngroups)[1:sparsity])

#         for i in 1:ngroups-1
#             !(i in inds_nz) && (x0[T*i-3:T*i] .= 0)
#         end
#         !(ngroups in inds_nz) && (x0[T*ngroups-3:end] .= 0)

function get_lasso_MLE(;n=20, m=15, sparsity=0.5, δ=0.1, seed=1234)
    @assert 0 ≤ sparsity ≤ 1

    Random.seed!(seed)
    A = rand(Normal(), m, n)

    Random.seed!(seed+1)
    x0 = rand(Normal(), n)
    Random.seed!(seed+2)
    optstructure = l1Manifold(rand(Bernoulli(1-sparsity), n))
    x0 = project(optstructure, x0)

    # measures are linear model + (unbiaised) noise
    y = A*x0 + rand(Normal(0, δ^2), m)

    return LeastsquaresPb(A, y, regularizer_l1(δ), n, x0, optstructure)
end


text:
#Why the p in grid interfaces?
abstract type GridInterface{q <: QuadratureRule} end
abstract type GridVessel{q, B <: GridBuild, K, p}  <: GridInterface{q} end
abstract type CollapsedGrid{q <: QuadratureRule} end
struct FlatGrid{q,N,C} <: CollapsedGrid{q}
    nodes::N
    cache::C
    weights::Vector{Float64}
    density::Vector{Float64}
    baseline::Vector{Float64}
end
struct MatrixCache{p}
  d::Dict{Int,Array{Float64,2}}
end
struct DynamicGridVessel{q, B, K, p, N, C} <: GridVessel{q, B, K, p}
    grids::Dict{K, FlatGrid{q,N,C}}
    U::Matrix{Float64}
    mats::MatrixCache{p}
end
struct MatrixVecSVec{p, T <: Real}
    M::Matrix{T}
    S::Vector{SVector{p,T}}
end

ext: jl


text:
##########################################################################
# Copyright 2017 Samuel Ridler.
#
# 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.
##########################################################################

# Run a search algorithm to optimise ambulance deployment, given some objective function (objFn).
# Given an initial deployment, an ambulance is added where it provides the
# greatest benefit, then an ambulance that provides the least benefit is removed.
# The search stops when no further improvements can be made.
# Ambulances are assumed to be identical, and station capacities are ignored.

using JEMSS
using Random
using Statistics
using StatsFuns
using Base.Iterators: Stateful

const StationsNumAmbs = Vector{Int} # type alias

function deploymentSearch!(sim::Simulation, deployments::Vector{Deployment}; outputFolder::String = "")
	isdir(outputFolder) || mkdir(outputFolder)
	@assert(isdir(outputFolder))
	@assert(!isempty(deployments))

	# parameters
	global solFilename = "$outputFolder/solutions.csv" # final solutions from each search
	global deploymentsOutputFilename = "$outputFolder/deployments.csv" # solutions, stored as deployments
	global logFilename = "$outputFolder/log.csv" # log search progress to file
	global sense = :max # :min or :max; direction of optimisation for objective function
	global conf = 0.95 # statistical confidence level
	global doPrint = true
	# parallel (multi-threading):
	global parallel = false
	global numThreads = 1

	# some parameter checks
	@assert(sense == :min || sense == :max)
	@assert(0 <= conf < 1) # should be between 0 and 1, but not equal to 1 otherwise confidence interval is infinite
	@assert(isa(parallel, Bool))
	@assert(numThreads >= 1)

	global nullObjVal = MeanAndHalfWidth(NaN,NaN)

	# keep track of station ambulance counts tried, and their objective values
	global stationsNumAmbsObjVal = Dict{StationsNumAmbs, MeanAndHalfWidth}()
	global stationsNumAmbsStats = Dict{StationsNumAmbs, Dict{String, Any}}()

	# save the station ambulance counts found for each search
	global stationsNumAmbsSols = Vector{StationsNumAmbs}() # stationsNumAmbsSols[i] is solution for deployments[i]

	repeatedSearch!(sim, deployments)
end

function objValLookup(stationsNumAmbs::StationsNumAmbs)::MeanAndHalfWidth
	return get(stationsNumAmbsObjVal, stationsNumAmbs, nullObjVal) # return objective value if found, otherwise nullObjVal
end

function meanAndHalfWidth(x::Vector{T}; conf::Float = conf)::MeanAndHalfWidth where T <: Real
	x = convert(Vector{Float}, x)
	return MeanAndHalfWidth(mean(x), tDistrHalfWidth(x; conf = conf))
end

# Objective value for a single period
function objFn(period::SimPeriodStats)::Float
	return period.call.numResponsesInTime / period.call.numCalls # fraction of calls reached in time
end

# For completed simulation replications, calculate and return the objective function value mean and half-width of the mean.
function objFn(sim::Simulation)::MeanAndHalfWidth
	@assert(all(rep -> rep.complete, sim.reps))
	periods = getRepsPeriodStatsList(sim.reps)
	return meanAndHalfWidth([objFn(p) for p in periods])
end

# add ambulance to sim at given station (but not from sim.backup)
function addAmb!(sim::Simulation; stationIndex::Int = 1)
	@assert(!sim.used)
	amb = Ambulance() # deepcopy(sim.ambulances[1])
	amb.index = sim.numAmbs + 1
	amb.stationIndex = stationIndex
	JEMSS.initAmbulance!(sim, amb)

	push!(sim.ambulances, amb)
	sim.numAmbs += 1
end

# remove ambulance from sim (but not from sim.backup)
function removeAmb!(sim::Simulation; ambIndex::Int = sim.numAmbs)
	@assert(!sim.used)
	deleteat!(sim.ambulances, ambIndex)
	sim.numAmbs -= 1
	for event in filter(e -> e.ambIndex == ambIndex, sim.eventList)
		JEMSS.deleteEvent!(sim.eventList, event)
	end
end

# apply stationsNumAmbs to sim replications, and return objective value (and half-width) and whether lookup was used
# uses look-up if possible
# mutates: sim, stationsNumAmbsObjVal
function simObjVal!(sim::Simulation, stationsNumAmbs::StationsNumAmbs)::Tuple{Tuple{Float,Float},Bool}
	global stationsNumAmbsObjVal, parallel, numThreads
	objVal = objValLookup(stationsNumAmbs)
	if objVal != nullObjVal return ((objVal.mean, objVal.halfWidth), true) end
	numAmbsDiff = sum(stationsNumAmbs) - sim.numAmbs # number of ambulances to add to sim
	function runRep!(rep::Simulation)
		resetRep!(rep)
		for i = 1:numAmbsDiff addAmb!(rep) end
		for i = 1:-numAmbsDiff removeAmb!(rep) end
		applyStationsNumAmbs!(rep, stationsNumAmbs)
		simulateRep!(rep)
		@assert(getStationsNumAmbs(rep) == stationsNumAmbs) # check that rep had the current deployment applied
	end
	if parallel == true
		runParallel!(runRep!, sim.reps...; numThreads = numThreads)
	else
		for rep in sim.reps runRep!(rep) end
	end
	objVal = objFn(sim)
	stationsNumAmbsObjVal[stationsNumAmbs] = objVal
	return ((objVal.mean, objVal.halfWidth), false)
end

# print number of ambulances at each station
function printStationsNumAmbs(stationsNumAmbs::StationsNumAmbs)
	for i = 1:length(stationsNumAmbs)
		println("station ", i, ": ", stationsNumAmbs[i], " ambulance(s)")
	end
end

function getSimStats(sim::Simulation)
	# some sim statistics to write to log file
	global periodDuration, simStatsKeys, simStatsEmpty
	if all(r -> r.complete, sim.reps)
		periods = getRepsPeriodStatsList(sim.reps)
		@assert(all(p -> isapprox(p.duration, periodDuration), periods))
		stats = flatten(statsDictFromPeriodStatsList(periods))
		return merge(Dict(["$(key)_mean" => stats[key].mean for key in simStatsKeys]), Dict(["$(key)_halfWidth" => stats[key].halfWidth for key in simStatsKeys]))
	end
	return deepcopy(simStatsEmpty)
end

global simStatsEmpty = flatten(statsDictFromPeriodStatsList(SimPeriodStats[]))
global simStatsKeys = collect(keys(simStatsEmpty))
global simStatsEmpty = merge(Dict(["$(key)_mean" => NaN for key in simStatsKeys]), Dict(["$(key)_halfWidth" => NaN for key in simStatsKeys]))

global solFileHeader = vcat(["search", "objVal", "objValHalfWidth"], sort(collect(keys(simStatsEmpty)))) # ... and stationsNumAmbs
global logFileHeader = vcat(["search", "iter", "oper", "stationIndex", "usedLookup", "searchDurationSeconds", "objVal", "objValHalfWidth", "bestObjVal"], sort(collect(keys(simStatsEmpty)))) # ... and stationsNumAmbs
global logFileDict = Dict{String,Any}([s => "" for s in logFileHeader])

# mutates: file
function fileWriteDlmLine!(file::IOStream, fileHeader::Vector{String}, data::Dict{String,T}, stationsNumAmbs::StationsNumAmbs) where T <: Any
	line = []
	for header in fileHeader
		push!(line, data[header])
	end
	line = vcat(line, stationsNumAmbs)
	writeDlmLine!(file, line...)
	flush(file)
end

# perform search, starting at each of the deployments provided
# mutates: sim
function repeatedSearch!(sim::Simulation, deployments::Vector{Deployment})
	global doPrint

	# open files for writing solution
	global solFilename, logFilename
	solFile = open(solFilename, "w")
	logFile = open(logFilename, "w")

	# write misc data to files
	global warmUpDuration, periodDuration, conf
	numSearches = length(deployments)
	miscTable = Table("miscData", ["numAmbs", "numStations", "numSearches", "conf"];
		rows = [[sim.numAmbs, sim.numStations, numSearches, conf]]) # note that numCalls may vary between sim replications
	repsTable = Table("repsData", ["numReps", "warmUpDuration", "periodDuration"];
		rows = [[sim.numReps, warmUpDuration, periodDuration]])
	writeTablesToFile!(logFile, [miscTable, repsTable])
	writeTablesToFile!(solFile, [miscTable, repsTable])

	# write file headers
	writeDlmLine!(solFile, "solutions")
	writeDlmLine!(logFile, "log")
	writeDlmLine!(solFile, solFileHeader..., ["station_$i numAmbs" for i = 1:sim.numStations]...)
	writeDlmLine!(logFile, logFileHeader..., ["station_$i numAmbs" for i = 1:sim.numStations]...)
	flush(solFile)
	flush(logFile)

	for i = 1:numSearches
		logFileDict["search"] = i

		# initial deployment
		stationsNumAmbs = deploymentToStationsNumAmbs(deployments[i], sim.numStations)

		# perform search until no improvements can be made
		doPrint && println()
		doPrint && println("Search ", i, " (of ", numSearches, ")")
		doPrint && printStationsNumAmbs(stationsNumAmbs)
		stationsNumAmbs = search!(sim, stationsNumAmbs, logFile)
		doPrint && println("Finished search $i (of $numSearches); solution:")
		doPrint && printStationsNumAmbs(stationsNumAmbs)

		# save best station ambulance count found for this search
		global stationsNumAmbsSols
		push!(stationsNumAmbsSols, stationsNumAmbs)

		# write solution to file
		objValMeanAndHalfWidth = objValLookup(stationsNumAmbs)
		solFileDict = Dict("search" => i, "objVal" => objValMeanAndHalfWidth.mean, "objValHalfWidth" => objValMeanAndHalfWidth.halfWidth)
		merge!(solFileDict, stationsNumAmbsStats[stationsNumAmbs])
		fileWriteDlmLine!(solFile, solFileHeader, solFileDict, stationsNumAmbs)

		global deploymentsOutputFilename
		writeDeploymentsFile(deploymentsOutputFilename, map(stationsNumAmbsToDeployment, stationsNumAmbsSols), sim.numStations)
	end

	# print out all results from each finished search
	doPrint && println()
	doPrint && println("Station ambulance counts from completed searches:")
	for i = 1:numSearches
		doPrint && println()
		doPrint && println("Search: ", i)
		stationsNumAmbs = stationsNumAmbsSols[i]
		doPrint && printStationsNumAmbs(stationsNumAmbs)
		doPrint && println("objective value = ", objValLookup(stationsNumAmbs).mean)
	end

	close(solFile)
	close(logFile)
end

# Search of ambulance deployment, starting at stationsNumAmbs.
# Given an initial deployment, an ambulance is added where it provides the
# greatest benefit, then an ambulance that provides the least benefit is removed.
# Continue search until no improvement can be made.
# mutates: sim, logFile
function search!(sim::Simulation, stationsNumAmbs::StationsNumAmbs, logFile::IOStream)::StationsNumAmbs

	global logFileHeader, logFileDict, stationsNumAmbsStats, sense, doPrint

	# track iterations
	iter = 1 # count iterations performed
	startTime = time()
	getSearchDuration() = round(time() - startTime, digits = 2)

	# calculate objective value for starting point
	(objVal, objValHalfWidth), usedLookup = simObjVal!(sim, stationsNumAmbs)
	doPrint && println("starting objective value: ", objVal)
	bestObjVal = objVal # current best objective value
	bestStationsNumAmbs = copy(stationsNumAmbs)

	# write starting point
	for (key, value) in logFileDict
		if key != "search"
			logFileDict[key] = "" # reset value
		end
	end
	stats = getSimStats(sim)
	if !usedLookup stationsNumAmbsStats[stationsNumAmbs] = stats end
	merge!(logFileDict, stats)
	merge!(logFileDict, Dict("iter" => iter, "usedLookup" => Int(usedLookup), "searchDurationSeconds" => getSearchDuration(),
		"objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal))
	fileWriteDlmLine!(logFile, logFileHeader, logFileDict, stationsNumAmbs)

	# search
	doPrint && println("Starting search.")
	endPoint = nothing
	iter += 1
	numStations = sim.numStations # shorthand
	while true
		bestTempStationsNumAmbs = nothing # will hold new deployment for one more ambulance
		bestTempObjVal = 0
		addStationIndex = 0
		for i = 1:numStations
			# test addition of ambulance to station i
			stationsNumAmbs = copy(bestStationsNumAmbs)
			stationsNumAmbs[i] += 1
			doPrint && print("add: i = $i")

			(objVal, objValHalfWidth), usedLookup = simObjVal!(sim, stationsNumAmbs)
			if !usedLookup
				stats = stationsNumAmbsStats[stationsNumAmbs] = getSimStats(sim)
				merge!(logFileDict, stats)
			else
				merge!(logFileDict, stationsNumAmbsStats[stationsNumAmbs])
				doPrint && print("; done")
			end
			doPrint && println("; objective value: $objVal")

			if i == 1 || (sense == :max && bestTempObjVal < objVal) || (sense == :min && bestTempObjVal > objVal)
				bestTempObjVal = objVal
				bestTempStationsNumAmbs = stationsNumAmbs
				addStationIndex = i
			end

			merge!(logFileDict, Dict("iter" => iter, "oper" => "+", "stationIndex" => i, "usedLookup" => Int(usedLookup),
				"searchDurationSeconds" => getSearchDuration(), "objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal))
			fileWriteDlmLine!(logFile, logFileHeader, logFileDict, stationsNumAmbs)
		end
		doPrint && println("added ambulance to station: $addStationIndex")

		rmStationIndex = 0
		for i = 1:numStations
			# test removal of ambulance from station i
			stationsNumAmbs = copy(bestTempStationsNumAmbs)
			stationsNumAmbs[i] -= 1
			if stationsNumAmbs[i] < 0 continue end
			doPrint && print("rm: i = $i")

			(objVal, objValHalfWidth), usedLookup = simObjVal!(sim, stationsNumAmbs)
			if !usedLookup
				stats = stationsNumAmbsStats[stationsNumAmbs] = getSimStats(sim)
				merge!(logFileDict, stats)
			else
				merge!(logFileDict, stationsNumAmbsStats[stationsNumAmbs])
				doPrint && print("; done")
			end
			doPrint && println("; objective value: $objVal")

			if (sense == :max && bestObjVal < objVal) || (sense == :min && bestObjVal > objVal)
				bestObjVal = objVal
				bestStationsNumAmbs = stationsNumAmbs
				rmStationIndex = i
			end

			merge!(logFileDict, Dict("iter" => iter, "oper" => "-", "stationIndex" => i, "usedLookup" => Int(usedLookup),
				"searchDurationSeconds" => getSearchDuration(), "objVal" => objVal, "objValHalfWidth" => objValHalfWidth, "bestObjVal" => bestObjVal))
			fileWriteDlmLine!(logFile, logFileHeader, logFileDict, stationsNumAmbs)
		end

		if rmStationIndex == 0
			# no improvement to deployment, search is complete
			doPrint && println("Finished search.")
			return bestStationsNumAmbs
		else
			doPrint && println("removed ambulance from station: $rmStationIndex")
			iter += 1
		end
	end
end

# mutates: sim
function setSimStatsCapture!(sim::Simulation, periodDurationsIter::Stateful, warmUpDuration::Float)
	stats = sim.stats = deepcopy(sim.backup.stats)
	stats.doCapture = true
	(stats.periodDurationsIter, stats.warmUpDuration) = (periodDurationsIter, warmUpDuration)
	stats.nextCaptureTime = sim.startTime + (stats.warmUpDuration > 0 ? stats.warmUpDuration : first(stats.periodDurationsIter))
	sim.backup.stats = deepcopy(sim.stats)
end

# parameters
configFilename = "sim_config.xml"
outputFolder = "output"
numSearches = 1
deploymentRng = Random.MersenneTwister(0) # useful for reproducing results, if using random restarts
# replication parameters
warmUpDuration = 0.1
periodDuration = 1.0

# parameter checks
@assert(numSearches >= 1)
@assert(warmUpDuration >= 0)
@assert(periodDuration > 0)

# initialise sim
t = time()
@assert(isfile(configFilename))
println("Initialising simulation from config: ", configFilename)
sim = initSim(configFilename, createBackup = true)
@assert(sim.numReps >= 2, "Need at least two simulation replications (sim.reps) for statistics.")

# set sim stats capturing
setSimStatsCapture!(sim, Stateful(periodDuration), warmUpDuration)
periodEndTime = sim.startTime + warmUpDuration + periodDuration
for rep in sim.reps
	# statistics should stop collecting before last call arrives (if not earlier), to avoid cool-down period
	@assert(periodEndTime <= rep.calls[end].arrivalTime)
end

makeRepsRunnable!(sim) # copy changes to sim (stats) to sim.reps

# generate deployments
deployments = makeRandDeployments(sim, numSearches; rng = deploymentRng)

# run
deploymentSearch!(sim, deployments, outputFolder = outputFolder)
println("total runtime: ", round(time()-t, digits = 2), " seconds")


text:
using FixedSizeArrays

type Heightfield
    w::Float64
    h::Float64
    y::Matrix{Float64}
end

Heightfield(w::Float64,h::Float64, nx::Integer, nz::Integer) =
    Heightfield(w, h, zeros(Float64, nx, nz))

Heightfield(hf::Heightfield) = Heightfield(hf.w, hf.h, copy(hf.y))

function triangle_list(hf::Heightfield)
    nx,nz = size(hf.y)
    x = linspace(-1,1,nx)*hf.w/2
    z = linspace(-1,1,nz)*hf.h/2
    dx = hf.w/(nx-1)
    dz = hf.h/(nz-1)

    vertices = Vector{Point{3,Float64}}()
    for xi in 1:nx
        for zi in 1:nz
            push!(vertices, Point(x[xi], hf.y[xi,zi], z[zi]))
        end
    end

    Dx = Tridiagonal(-ones(nx-1),zeros(nx),ones(nx-1))
    Dz = Tridiagonal(-ones(nz-1),zeros(nz),ones(nz-1))

    dydx = Dx*hf.y
    dydz = (Dz*hf.y')'

    normals = Vector{Vec{3,Float64}}()
    for xi in 1:nx
        for zi in 1:nz
            a = Vec(2dx, dydx[xi,zi], 0.)
            b = Vec(0., dydz[xi,zi], 2dz)
            push!(normals, -normalize(cross(a,b)))
        end
    end

    faces = Vector{Tuple}()
    for xi in 0:nx-2
        for zi in 1:nz-1
            push!(faces, (zi + xi*nz, zi + (xi+1)*nz, zi + 1 + (xi+1)*nz))
            push!(faces, (zi + xi*nz, zi + 1 + (xi+1)*nz, zi + 1 + xi*nz))
        end
    end

    vertices, normals, faces
end

export Heightfield, view, triangle_list


text:
using Plots
using RecipesBase
import StatsBase: countmap

@userplot CountPlot

@recipe function f(cp::CountPlot)
    bar_width --> 0.75
    color --> :skyblue
    legend --> nothing
    seriestype := :bar

    data_points = cp.args[1]
    counted_data = sort(countmap(data_points))

    label_ticks = keys(counted_data) |> collect
    count_values = values(counted_data) |> collect

    given_orientation = get(plotattributes, :orientation, nothing)

    if given_orientation == :horizontal
        yticks := get(plotattributes, :yticks, nothing) ==
            nothing ? label_ticks : nothing
        y, x = label_ticks, count_values
    else
        xticks := get(plotattributes, :xticks, nothing) ==
            nothing ? label_ticks : nothing
        x, y = label_ticks, count_values
    end
end


text:
using .Networks

using JuMP
using Dates
using DataStructures
using Printf
using Parameters
"""
REF_SCHEDULE_TYPE_IN_TSO : Indicates which schedule to use as reference for pilotables state/levels needed
                            for sequencing constraints and TSO objective function.
"""
@with_kw mutable struct TSOBilevelConfigs
    CONSIDER_DELTAS::Bool = true
    CONSIDER_N_1_CSTRS::Bool = false
    TSO_LIMIT_PENALTY::Float64 = 1e-3
    TSO_LOL_PENALTY::Float64 = 1e5
    TSO_CAPPING_COST::Float64 = 1.
    TSO_PILOTABLE_BOUNDING_COST::Float64 = 1.
    USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST::Bool = true
    MARKET_LOL_PENALTY::Float64 = 1e5
    MARKET_CAPPING_COST::Float64 = 1.
    out_path::Union{Nothing,String} = nothing
    problem_name::String = "TSOBilevel"
    LINK_SCENARIOS_LIMIT::Bool = true
    LINK_SCENARIOS_PILOTABLE_LEVEL::Bool = false
    LINK_SCENARIOS_PILOTABLE_ON::Bool = false
    LINK_SCENARIOS_PILOTABLE_LEVEL_MARKET::Bool = false
    big_m = 1e6
    REF_SCHEDULE_TYPE_IN_TSO::Union{Market,TSO} = Market();
end

##########################################################
#        upper problem : TSO
##########################################################

@with_kw struct TSOBilevelTSOLimitableModel <: AbstractLimitableModel
    #gen,ts,s
    delta_p = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_injected = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_limit = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    b_is_limited = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_limit_x_is_limited = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_capping = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #ts,s
    p_global_capping = SortedDict{Tuple{DateTime,String},VariableRef}();
end

@with_kw struct TSOBilevelTSOPilotableModel <: AbstractPilotableModel
    #gen,ts,s
    p_imposition_min = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_imposition_max = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    b_start = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    b_on = SortedDict{Tuple{String,DateTime,String},VariableRef}();
end

@with_kw struct TSOBilevelTSOLoLModel <: AbstractLoLModel
    #bus,ts,s #Loss of Load
    p_loss_of_load = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #ts,s
    p_global_loss_of_load = SortedDict{Tuple{DateTime,String},VariableRef}();
end

@with_kw mutable struct TSOBilevelTSOObjectiveModel <: AbstractObjectiveModel
    deltas = GenericAffExpr{Float64, VariableRef}()

    limitable_cost = GenericAffExpr{Float64, VariableRef}()
    pilotable_cost = GenericAffExpr{Float64, VariableRef}()

    penalty = GenericAffExpr{Float64, VariableRef}()

    full_obj_1 = GenericAffExpr{Float64, VariableRef}()
    full_obj_2 = GenericAffExpr{Float64, VariableRef}()
end

##########################################################
#        lower problem : Market
##########################################################
@with_kw struct TSOBilevelMarketLimitableModel <: AbstractLimitableModel
    #ts,s #FIXME not sure if this should be a lower or an upper variable
    p_global_injected = SortedDict{Tuple{DateTime,String},VariableRef}();
    #ts,s
    p_global_capping = SortedDict{Tuple{DateTime,String},VariableRef}();
end

"""
# Note :
    deltas are added here to allow automatically computing them and their constraints
    but the constraints are not included in the reformulation (stationarity nor complementarity)
    cause they can be seen as part of the upper
    constraints discussed here are : delta >= +- (p_injected - P_ref)
"""
@with_kw struct TSOBilevelMarketPilotableModel <: AbstractPilotableModel
    #gen,ts,s
    delta_p = SortedDict{Tuple{String,DateTime,String},VariableRef}();
    #gen,ts,s
    p_injected = SortedDict{Tuple{String,DateTime,String},VariableRef}();
end

@with_kw struct TSOBilevelMarketLoLModel <: AbstractLoLModel
    #ts,s
    p_global_loss_of_load = SortedDict{Tuple{DateTime,String},VariableRef}();
end

@with_kw mutable struct TSOBilevelMarketObjectiveModel <: AbstractObjectiveModel
    limitable_cost = GenericAffExpr{Float64, VariableRef}()
    pilotable_cost = GenericAffExpr{Float64, VariableRef}()

    penalty = GenericAffExpr{Float64, VariableRef}()

    full_obj = GenericAffExpr{Float64, VariableRef}()
end

##########################################################
#        TSOBilevel Model
##########################################################

@with_kw struct TSOBilevelTSOModelContainer <: AbstractModelContainer
    model::Model
    limitable_model::TSOBilevelTSOLimitableModel = TSOBilevelTSOLimitableModel()
    pilotable_model::TSOBilevelTSOPilotableModel = TSOBilevelTSOPilotableModel()
    lol_model::TSOBilevelTSOLoLModel = TSOBilevelTSOLoLModel()
    objective_model::TSOBilevelTSOObjectiveModel = TSOBilevelTSOObjectiveModel()
    #branch,ts,s,ptdf_case
    flows::SortedDict{Tuple{String,DateTime,String,String},AffExpr} =
        SortedDict{Tuple{String,DateTime,String,String},AffExpr}()
    rso_constraint::SortedDict{Tuple{String,DateTime,String,String},ConstraintRef} =
        SortedDict{Tuple{String,DateTime,String,String},ConstraintRef}()
end
@with_kw struct TSOBilevelMarketModelContainer <: AbstractModelContainer
    model::Model
    limitable_model::TSOBilevelMarketLimitableModel = TSOBilevelMarketLimitableModel()
    pilotable_model::TSOBilevelMarketPilotableModel = TSOBilevelMarketPilotableModel()
    lol_model::TSOBilevelMarketLoLModel = TSOBilevelMarketLoLModel()
    objective_model::TSOBilevelMarketObjectiveModel = TSOBilevelMarketObjectiveModel()
    #ts,s
    eod_constraint::SortedDict{Tuple{Dates.DateTime,String}, ConstraintRef} =
        SortedDict{Tuple{Dates.DateTime,String}, ConstraintRef}()
end
@with_kw struct TSOBilevelKKTModelContainer <: AbstractModelContainer
    #TODO create adequate struct to link each dual kkt variable, indicator variable, and constraint to each other
    # e.g.
    # eod_model = ts,s -> Struct{dual, indicator, ConstraintRef}
    # pilotable_min = id,ts,s -> Struct{dual, indicator, ConstraintRef} or Struct{dual, Nothing, ConstraintRef}
    # than simply call reformulate_kkt(cstr) when building lower problem's constraints provided that objective is created first
    model::Model
    #ts,s
    eod_duals::SortedDict{Tuple{DateTime,String},VariableRef} =
        SortedDict{Tuple{DateTime,String},VariableRef}()
    capping_duals::SortedDict{Tuple{DateTime,String},VariableRef} =
        SortedDict{Tuple{DateTime,String},VariableRef}()
    capping_indicators::SortedDict{Tuple{DateTime,String},VariableRef} =
        SortedDict{Tuple{DateTime,String},VariableRef}()
    loss_of_load_duals::SortedDict{Tuple{DateTime,String},VariableRef} =
        SortedDict{Tuple{DateTime,String},VariableRef}()
    loss_of_load_indicators::SortedDict{Tuple{DateTime,String},VariableRef} =
        SortedDict{Tuple{DateTime,String},VariableRef}()
    #pilotable_id,ts,s
    firmness_duals::SortedDict{Tuple{String,DateTime,String},VariableRef} =
        SortedDict{Tuple{String,DateTime,String},VariableRef}()
    pmin_duals::SortedDict{Tuple{String,DateTime,String},VariableRef} =
        SortedDict{Tuple{String,DateTime,String},VariableRef}()
    pmin_indicators::SortedDict{Tuple{String,DateTime,String},VariableRef} =
        SortedDict{Tuple{String,DateTime,String},VariableRef}()
    pmax_duals::SortedDict{Tuple{String,DateTime,String},VariableRef} =
        SortedDict{Tuple{String,DateTime,String},VariableRef}()
    pmax_indicators::SortedDict{Tuple{String,DateTime,String},VariableRef} =
        SortedDict{Tuple{String,DateTime,String},VariableRef}()
end
TSOBilevelModel = BilevelModelContainer{TSOBilevelTSOModelContainer, TSOBilevelMarketModelContainer, TSOBilevelKKTModelContainer}
function BilevelModelContainer{TSOBilevelTSOModelContainer,TSOBilevelMarketModelContainer,TSOBilevelKKTModelContainer}()
    bilevel_model = Model()
    upper = TSOBilevelTSOModelContainer(model=bilevel_model)
    lower = TSOBilevelMarketModelContainer(model=bilevel_model)
    kkt_model = TSOBilevelKKTModelContainer(model=bilevel_model)
    return BilevelModelContainer(bilevel_model, upper, lower, kkt_model)
end

"""
Computes the actual capping in the Limitable model

    The global capping of the upper limitable model can be lower than the actual capping
    The actual capping is given by the localised capping in the upper problem or the global capping in the lower problem
"""
function sum_capping(limitable_model::TSOBilevelTSOLimitableModel, ts,s, network::Networks.Network)
    sum_l = 0.
    for gen in Networks.get_generators_of_type(network, Networks.LIMITABLE)
        gen_id = Networks.get_id(gen)
        sum_l += limitable_model.p_capping[gen_id,ts,s]
    end
    return sum_l
end

"""
Computes the actual loss of load (LoL) in the LoL model

    The global LoL of the upper limitable model can be lower than the actual LoL
    The actual LoL is given by the localised LoL in the upper problem or the global LoL in the lower problem
"""
function sum_lol(lol_model::TSOBilevelTSOLoLModel, ts, s, network::Networks.Network)
    sum_l = 0.
    for bus in Networks.get_buses(network)
        bus_id = Networks.get_id(bus)
        sum_l += lol_model.p_loss_of_load[bus_id,ts,s]
    end
    return sum_l
end

function has_positive_slack(model_container::TSOBilevelModel)::Bool
    return has_positive_value(model_container.lower.lol_model.p_global_loss_of_load) #If TSO cut => market did too
end

function has_positive_slack(model_container::TSOBilevelTSOModelContainer)::Bool
    return has_positive_value(model_container.lol_model.p_loss_of_load) #If market cut => TSO localised the LoL
end

function get_upper_obj_expr(bilevel_model::TSOBilevelModel)
    return bilevel_model.upper.objective_model.full_obj_2
end
function get_lower_obj_expr(bilevel_model::TSOBilevelModel)
    return bilevel_model.lower.objective_model.full_obj
end

##########################################################
#        upper problem : TSO functions : TSOBilevelTSOModelContainer
##########################################################
function create_tso_vars!( model_container::TSOBilevelTSOModelContainer,
                            network::Networks.Network,
                            target_timepoints::Vector{Dates.DateTime},
                            scenarios::Vector{String},
                            preceding_market_schedule::Schedule
                            )
    pilotables_list_l = Networks.get_generators_of_type(network, Networks.PILOTABLE)
    limitables_list_l = Networks.get_generators_of_type(network, Networks.LIMITABLE)
    buses_list_l = Networks.get_buses(network)
    add_pilotables_vars!(model_container,
                        pilotables_list_l, target_timepoints, scenarios,
                        imposition_vars=true, commitment_vars=true,
                        prefix="tso")
    add_limitables_vars!(model_container,
                        target_timepoints, scenarios,
                        limitables_list_l,
                        preceding_market_schedule,
                        delta_vars=true,
                        injection_vars=true, limit_vars=true, local_capping_vars=true, global_capping_vars=true,
                        prefix="tso"
                        )
    add_lol_vars!(model_container,
                target_timepoints, scenarios, buses_list_l,
                global_lol_vars=true, local_lol_vars=true,
                prefix="tso")
end

function add_tso_constraints!(bimodel_container::TSOBilevelModel,
                            target_timepoints, scenarios, network,
                            firmness, reference_schedule, generators_initial_state,
                            uncertainties_at_ech::UncertaintiesAtEch,
                            configs::TSOBilevelConfigs)
    model = bimodel_container.model
    tso_model_container::TSOBilevelTSOModelContainer = bimodel_container.upper
    market_model_container::TSOBilevelMarketModelContainer = bimodel_container.lower

    pilotables_list_l = Networks.get_generators_of_type(network, Networks.PILOTABLE)
    unit_commitment_constraints!(model,
                    tso_model_container.pilotable_model, pilotables_list_l,  target_timepoints, scenarios,
                    firmness, reference_schedule, generators_initial_state,
                    always_link_scenarios=configs.LINK_SCENARIOS_PILOTABLE_ON,
                    prefix="tso")
    power_imposition_constraints!(model, tso_model_container.pilotable_model,
                    pilotables_list_l, target_timepoints, scenarios,
                    firmness, reference_schedule;
                    always_link_scenarios=configs.LINK_SCENARIOS_PILOTABLE_LEVEL,
                    prefix="tso")

    limitables_list_l = Networks.get_generators_of_type(network, Networks.LIMITABLE)
    limitables_ids_l = map(lim_gen->Networks.get_id(lim_gen), limitables_list_l)
    limitable_model = tso_model_container.limitable_model
    limitable_power_constraints!(model,
                        limitable_model, limitables_list_l, target_timepoints, scenarios,
                        firmness, uncertainties_at_ech,
                        always_link_scenarios=configs.LINK_SCENARIOS_LIMIT,
                        prefix="tso")
    local_capping_constraints!(model,
                            limitable_model, limitables_list_l, target_timepoints, scenarios,
                            uncertainties_at_ech, prefix="tso")
    global_capping_constraints!(model,
                            limitable_model, limitables_list_l,
                            target_timepoints, scenarios,
                            uncertainties_at_ech,
                            prefix="tso")
    distribution_constraint!(model,
                            get_global_capping(market_model_container.limitable_model),
                            get_local_capping(tso_model_container.limitable_model),
                            limitables_ids_l, target_timepoints, scenarios,
                            cstr_prefix_name="tso_distribute_capping")
    distribution_constraint!(model,
                            market_model_container.limitable_model.p_global_injected, #TODO rename to avoid onfusion between localised injections and global ones
                            get_p_injected(tso_model_container.limitable_model),
                            limitables_ids_l, target_timepoints, scenarios,
                            cstr_prefix_name="tso_distribute_enr_injections")

    buses_list_l = Networks.get_buses(network)
    buses_ids_l = map(bus->Networks.get_id(bus), buses_list_l)
    lol_model = tso_model_container.lol_model
    local_lol_constraints!(model,
                        lol_model,
                        buses_list_l, target_timepoints, scenarios,
                        uncertainties_at_ech,
                        prefix="tso")
    global_lol_constraints!(model,
                        lol_model, buses_list_l, target_timepoints, scenarios,
                        uncertainties_at_ech,
                        prefix="tso")
    distribution_constraint!(model,
                            get_global_lol(market_model_container.lol_model),
                            get_local_lol(tso_model_container.lol_model),
                            buses_ids_l, target_timepoints, scenarios,
                            cstr_prefix_name="tso_distribute_lol")

    add_rso_flows_exprs(tso_model_container.flows,
                market_model_container.pilotable_model,
                tso_model_container.limitable_model,
                tso_model_container.lol_model,
                all_combinations(network, target_timepoints, scenarios),
                uncertainties_at_ech, network)
    combinations = (configs.CONSIDER_N_1_CSTRS) ? all_combinations(network, target_timepoints, scenarios) :
                                                    all_n_combinations(network, target_timepoints, scenarios)
    rso_constraints!(bimodel_container.model,
                    tso_model_container.rso_constraint,
                    tso_model_container.flows,
                    combinations,
                    network,
                    prefix="tso")

    return bimodel_container
end

function create_tso_objectives!(bimodel_container::TSOBilevelModel,
                                target_timepoints, scenarios, network,
                                preceding_market_schedule::Schedule,
                                capping_cost, loss_of_load_cost,
                                limit_penalty,
                                pilotable_bounding_cost,
                                use_prop_cost_for_bounding::Bool)
    tso_model_container = bimodel_container.upper
    market_model_container = bimodel_container.lower
    objective_model = tso_model_container.objective_model

    # cost for deviating from market schedule
    for (_, var_delta) in market_model_container.pilotable_model.delta_p
        objective_model.deltas += var_delta
    end
    for (_, var_delta) in tso_model_container.limitable_model.delta_p
        objective_model.deltas += var_delta
    end

    # objective_model.pilotable_cost
    add_tsobilevel_impositions_cost!(tso_model_container,
                                    target_timepoints, scenarios, network,
                                    preceding_market_schedule,
                                    pilotable_bounding_cost,
                                    use_prop_cost_for_bounding)

    # limitable_cost : capping (fr. ecretement)
    objective_model.limitable_cost += coeffxsum(tso_model_container.limitable_model.p_global_capping, capping_cost)

    # cost for cutting consumption (lol) and avoid limiting for no reason
    objective_model.penalty += coeffxsum(tso_model_container.lol_model.p_global_loss_of_load, loss_of_load_cost)
    objective_model.penalty += coeffxsum(tso_model_container.limitable_model.b_is_limited, limit_penalty)

    # cost for cutting load/consumption by the market
    objective_model.penalty += coeffxsum(market_model_container.lol_model.p_global_loss_of_load, loss_of_load_cost)

    #Objective 1 : add penalty ?
    objective_model.full_obj_1 = objective_model.deltas + objective_model.penalty
    #Objective 2 : add penalty ?
    objective_model.full_obj_2 = ( objective_model.pilotable_cost +
                                objective_model.limitable_cost +
                                objective_model.penalty )

    return tso_model_container
end

function add_tsobilevel_impositions_cost!(model_container::TSOBilevelTSOModelContainer,
                                        target_timepoints, scenarios, network,
                                        preceding_market_schedule,
                                        pilotable_bounding_cost,
                                        use_prop_cost_for_bounding::Bool)
    tso_pilotable_model = model_container.pilotable_model
    objective_expr = model_container.objective_model.pilotable_cost

    for gen in Networks.get_generators_of_type(network, Networks.PILOTABLE)
        gen_id = Networks.get_id(gen)
        if use_prop_cost_for_bounding
            cost = Networks.get_prop_cost(gen)
        else
            cost = pilotable_bounding_cost
        end

        for ts in target_timepoints
            for s in scenarios
                commitment = get_commitment_value(preceding_market_schedule, gen_id, ts, s)
                if  !Networks.needs_commitment(gen) || (!ismissing(commitment) && (commitment==ON))
                    add_tsobilevel_started_impositions_cost!(objective_expr, gen,
                                                            tso_pilotable_model.p_imposition_min[gen_id, ts, s],
                                                            tso_pilotable_model.p_imposition_max[gen_id, ts, s],
                                                            cost)
                else
                    add_tsobilevel_non_started_impositions_cost!(objective_expr, model_container,
                                                                 gen, ts, s, cost)
                end
            end
        end
    end

    return model_container
end
function add_tsobilevel_started_impositions_cost!(objective_expr::AffExpr,
                                                gen::Generator,
                                                pmin_var, pmax_var,
                                                pilotable_bounding_cost
                                                )
    p_max = Networks.get_p_max(gen)
    p_min = Networks.get_p_min(gen)

    add_to_expression!(objective_expr, pilotable_bounding_cost * (p_max - pmax_var))
    add_to_expression!(objective_expr, pilotable_bounding_cost * (pmin_var - p_min))

    return objective_expr
end
function add_tsobilevel_non_started_impositions_cost!(objective_expr::AffExpr,
                                                      tso_model_container::TSOBilevelTSOModelContainer,
                                                    gen, ts, s,
                                                    pilotable_bounding_cost
                                                    )
    #need a second expression for units that do not need a commitment (p_min=0 and no b_on)
    @assert Networks.needs_commitment(gen)

    tso_pilotable_model = tso_model_container.pilotable_model
    gen_id = Networks.get_id(gen)

    p_max = Networks.get_p_max(gen)
    # p_min = Networks.get_p_min(gen)

    #need to cost reducing pmax otherwise TSO may always limit pmax when starting a unit
    pmax_var = tso_pilotable_model.p_imposition_max[gen_id, ts, s]
    b_on_var = tso_pilotable_model.b_on[gen_id, ts, s]
    add_to_expression!(objective_expr, pilotable_bounding_cost * (p_max*b_on_var - pmax_var) )

    pmin_var = tso_pilotable_model.p_imposition_min[gen_id, ts, s]
    add_to_expression!(objective_expr, pilotable_bounding_cost * pmin_var)
    # add_to_expression!(objective_expr, pilotable_bounding_cost * p_min * b_on_var)

    return objective_expr
end

##########################################################
#        lower problem : Market functtons  : TSOBilevelMarketModelContainer
##########################################################

function create_market_vars!(model_container::TSOBilevelMarketModelContainer,
                            network::Networks.Network,
                            target_timepoints::Vector{Dates.DateTime},
                            scenarios::Vector{String},
                            preceding_market_schedule::Schedule
                            )

    pilotables_list_l = Networks.get_generators_of_type(network, Networks.PILOTABLE)
    add_pilotables_vars!(model_container,
                        pilotables_list_l, target_timepoints, scenarios,
                        preceding_market_schedule,
                        delta_vars=true,
                        injection_vars=true,
                        prefix="market_"
                        )
    add_limitables_vars!(model_container,
                        target_timepoints, scenarios,
                        global_capping_vars=true,
                        global_injection_vars=true,
                        prefix="market_"
                        )
    add_lol_vars!(model_container,
                target_timepoints, scenarios,
                global_lol_vars=true,
                prefix="market_")
end

function add_firmness_duals!(kkt_model_container::TSOBilevelKKTModelContainer,
                            pilotables_list, target_timepoints, scenarios,
                            firmness,
                            link_scenarios_pilotable_level_market;
                            prefix::String="")
    for gen_pilotable in pilotables_list
        gen_id = Networks.get_id(gen_pilotable)
        for ts in target_timepoints
            decision_firmness_l = get_power_level_firmness(firmness, gen_id, ts)
            for s in scenarios[2:end]
                if requires_linking(decision_firmness_l, link_scenarios_pilotable_level_market)
                    #create duals relative to firmness constraints of injected pilotable power in market
                    name = @sprintf("c_firmness[%s,%s,%s]",gen_id,ts,s)
                    add_dual!(kkt_model_container.model, kkt_model_container.firmness_duals,
                                (gen_id,ts,s), name, false)
                end
            end
        end
    end
    return kkt_model_container
end

function add_eod_constraints_duals!(kkt_model_container::TSOBilevelKKTModelContainer,
                            target_timepoints, scenarios)
    for ts in target_timepoints
        for s in scenarios
            #create duals relative to EOD constraint
            name = @sprintf("c_eod[%s,%s]",ts,s)
            add_dual!(kkt_model_container.model, kkt_model_container.eod_duals, (ts,s), name, false)
        end
    end
end

function add_min_global_capping_duals!(kkt_model_container::TSOBilevelKKTModelContainer,
                                    target_timepoints, scenarios)
    for ts in target_timepoints
        for s in scenarios
            #create duals and indicators relative to TSO min capping constraint
            name = @sprintf("c_min_e[%s,%s]",ts,s)
            add_dual_and_indicator!(kkt_model_container.model,
                                    kkt_model_container.capping_duals, kkt_model_container.capping_indicators, (ts,s),
                                    name, true)
        end
    end
end

function add_min_global_lol_duals!(market_model_container::TSOBilevelMarketModelContainer,
                                kkt_model_container::TSOBilevelKKTModelContainer,
                                target_timepoints, scenarios)
    for ts in target_timepoints
        for s in scenarios
            #create duals and indicators relative to TSO min LoL constraint
            name = @sprintf("c_min_lol[%s,%s]",ts,s)
            add_dual_and_indicator!(kkt_model_container.model,
                                    kkt_model_container.loss_of_load_duals, kkt_model_container.loss_of_load_indicators, (ts,s),
                                    name, true)
        end
    end
    return market_model_container
end
function add_impositions_duals!(kkt_model_container::TSOBilevelKKTModelContainer,
                            pilotables_list, target_timepoints, scenarios)
    for ts in target_timepoints
        for s in scenarios
            for gen in pilotables_list
                gen_id = Networks.get_id(gen)
                add_pilotable_pmin_duals!(kkt_model_container,
                                        gen_id, ts, s)
                add_pilotable_pmax_duals!(kkt_model_container,
                                        gen_id, ts, s)
            end
        end
    end

    return kkt_model_container
end
function add_pilotable_pmin_duals!(kkt_model_container::TSOBilevelKKTModelContainer,
                                gen_id, ts, s)
    #create duals and indicators relative to tso pmin constraint
    name = @sprintf("c_tso_pmin[%s,%s,%s]",gen_id,ts,s)
    add_dual_and_indicator!(kkt_model_container.model,
                            kkt_model_container.pmin_duals, kkt_model_container.pmin_indicators, (gen_id,ts,s),
                            name, true)

    return kkt_model_container
end
function add_pilotable_pmax_duals!(kkt_model_container::TSOBilevelKKTModelContainer,
                                gen_id, ts, s)
    #create duals and indicators relative to tso pmax constraint
    name = @sprintf("c_tso_pmax[%s,%s,%s]",gen_id,ts,s)
    add_dual_and_indicator!(kkt_model_container.model,
                            kkt_model_container.pmax_duals, kkt_model_container.pmax_indicators, (gen_id,ts,s),
                            name, true)

    return kkt_model_container
end

function add_market_constraints!(bimodel_container::TSOBilevelModel,
                            target_timepoints, scenarios, network,
                            firmness, #reference_schedule, generators_initial_state,
                            uncertainties_at_ech::UncertaintiesAtEch,
                            configs::TSOBilevelConfigs)
    model = bimodel_container.model
    tso_model_container::TSOBilevelTSOModelContainer = bimodel_container.upper
    market_model_container::TSOBilevelMarketModelContainer = bimodel_container.lower
    kkt_model_container::TSOBilevelKKTModelContainer = bimodel_container.kkt_model

    pilotables_list_l = Networks.get_generators_of_type(network, Networks.PILOTABLE)
    pilotable_power_constraints!(model,
                                market_model_container.pilotable_model,
                                pilotables_list_l, target_timepoints, scenarios,
                                firmness,
                                missing,
                                ignore_sequencing_cstrs=true, #Sequencing constraints are imposed in TSO upper problem
                                always_link_scenarios=configs.LINK_SCENARIOS_PILOTABLE_LEVEL_MARKET,
                                prefix="market_")
    add_firmness_duals!(kkt_model_container,
                        pilotables_list_l, target_timepoints, scenarios, firmness,
                        configs.LINK_SCENARIOS_PILOTABLE_LEVEL_MARKET,
                        prefix="market_")
    respect_impositions_constraints!(model,
                                market_model_container.pilotable_model, #Injections are in the lower problem
                                tso_model_container.pilotable_model, #Impositions are in the upper problem
                                pilotables_list_l,  target_timepoints, scenarios,
                                prefix="market_")
    add_impositions_duals!(kkt_model_container,
                        pilotables_list_l, target_timepoints, scenarios)

    limitables_list_l = Networks.get_generators_of_type(network, Networks.LIMITABLE)
    global_capping_constraints!(model,
                            market_model_container.limitable_model,
                            limitables_list_l, target_timepoints, scenarios,
                            uncertainties_at_ech,
                            min_cap=get_global_capping(tso_model_container.limitable_model)
                            )
    add_min_global_capping_duals!(kkt_model_container,
                            target_timepoints, scenarios)
    global_injected_constraints!(model,
                            market_model_container.limitable_model,
                            limitables_list_l, target_timepoints, scenarios,
                            uncertainties_at_ech, prefix="market_")

    buses_list_l = Networks.get_buses(network)
    lol_model = market_model_container.lol_model
    global_lol_constraints!(model,
                        lol_model, buses_list_l, target_timepoints, scenarios,
                        uncertainties_at_ech,
                        min_lol=get_global_lol(tso_model_container.lol_model),
                        prefix="market")
    add_min_global_lol_duals!(market_model_container, kkt_model_container,
                        target_timepoints, scenarios)

    eod_constraints!(market_model_container.model, market_model_container.eod_constraint,
                    market_model_container.pilotable_model,
                    tso_model_container.limitable_model,
                    market_model_container.lol_model,
                    target_timepoints, scenarios,
                    uncertainties_at_ech, network
                    )
    add_eod_constraints_duals!(kkt_model_container, target_timepoints, scenarios)

    return bimodel_container
end

function create_market_objectives!(model_container::TSOBilevelMarketModelContainer,
                                network,
                                capping_cost, loss_of_load_cost)
    objective_model = model_container.objective_model

    # model_container.pilotable_cost
    add_prop_cost!(model_container.objective_model.pilotable_cost,
                            model_container.pilotable_model.p_injected, network)

    # limitable_cost : capping (fr. ecretement)
    objective_model.limitable_cost += coeffxsum(model_container.limitable_model.p_global_capping, capping_cost)

    # cost for cutting load/consumption
    objective_model.penalty += coeffxsum(model_container.lol_model.p_global_loss_of_load, loss_of_load_cost)

    objective_model.full_obj = ( objective_model.pilotable_cost +
                                objective_model.limitable_cost +
                                objective_model.penalty )

    return model_container
end

##########################################################
#        kkt reformulation : TSOBilevelKKTModelContainer
##########################################################

function add_kkt_stationarity_constraints!(kkt_model::TSOBilevelKKTModelContainer,
                                            target_timepoints, scenarios, network,
                                            capping_cost, loss_of_load_cost)
    #FIXME can be generic by iterating on lower variables to construct each stationarity constraint
    # iterate on the objective and lower constraints to extract their coefficients, but need to link each cnstraint to its dual var
    add_capping_stationarity_constraints!(kkt_model, target_timepoints, scenarios, capping_cost)
    add_loss_of_load_stationarity_constraints!(kkt_model, target_timepoints, scenarios, loss_of_load_cost)
    add_pilotable_stationarity_constraints!(kkt_model, target_timepoints, scenarios, network)
end

function add_loss_of_load_stationarity_constraints!(kkt_model::TSOBilevelKKTModelContainer,
                                                target_timepoints, scenarios, loss_of_load_cost)
    for ts in target_timepoints
        for s in scenarios
            # @assert ( capping_cost ≈ coefficient(market_model.objective_model.full_obj,
            #                                                 market_model.lol_model.p_global_loss_of_load[ts,s]) )
            name = @sprintf("c_stationarity_lol[%s,%s]",ts,s)
            @constraint(kkt_model.model,
                        loss_of_load_cost + kkt_model.eod_duals[ts,s] - kkt_model.loss_of_load_duals[ts,s] == 0,
                        base_name=name)
        end
    end
end

function add_capping_stationarity_constraints!(kkt_model::TSOBilevelKKTModelContainer,
                                                target_timepoints, scenarios, capping_cost)
    for ts in target_timepoints
        for s in scenarios
            # @assert ( capping_cost ≈ coefficient(market_model.objective_model.full_obj,
            #                                                 market_model.limitable_model.p_global_capping[ts,s]) )
            name = @sprintf("c_stationarity_e[%s,%s]",ts,s)
            @constraint(kkt_model.model,
                        capping_cost - kkt_model.eod_duals[ts,s] - kkt_model.capping_duals[ts,s] == 0,
                        base_name=name)
        end
    end
end

function firmness_duals_sationarity_expr(kkt_model,
                                        gen_id, ts, scenario,
                                        scenarios)::AffExpr
    if length(scenarios) == 1
        result_l = 0.
    elseif scenario == scenarios[1]
        sum_l = sum( get(kkt_model.firmness_duals, (gen_id,ts,s), 0.) for s in scenarios[2:end])
        result_l = -sum_l
    else
        result_l = get(kkt_model.firmness_duals, (gen_id,ts,scenario), 0.)
    end
    return result_l
end
function add_pilotable_stationarity_constraints!(kkt_model::TSOBilevelKKTModelContainer,
                                                target_timepoints, scenarios, network)
    for pilotable_gen in Networks.get_generators_of_type(network, Networks.PILOTABLE)
        gen_id = Networks.get_id(pilotable_gen)
        gen_prop_cost = Networks.get_prop_cost(pilotable_gen)
        for ts in target_timepoints
            for s in scenarios
                # @assert ( pilotable_bounding_cost ≈ coefficient(market_model.objective_model.full_obj,
                #                                             market_model.pilotable_model.p_injected[gen_id,ts,s]) )
                name = @sprintf("c_stationarity_pilotable_p[%s,%s,%s]",gen_id,ts,s)
                firmness_sationarity_expr = firmness_duals_sationarity_expr(kkt_model,
                                                                                gen_id, ts, s,
                                                                                scenarios)
                @constraint(kkt_model.model,
                            0 == gen_prop_cost
                                    + kkt_model.eod_duals[ts,s]
                                    - kkt_model.pmin_duals[gen_id,ts,s]
                                    + kkt_model.pmax_duals[gen_id,ts,s]
                                    + firmness_sationarity_expr,
                            base_name=name)
            end
        end
    end
end

function add_kkt_complementarity_constraints!(model_container::TSOBilevelModel,
                                            big_m, target_timepoints, scenarios, network)
    #FIXME can be done iteratively and generically if we loop on constraint expressions and know their corresponding kkt vars
    add_emin_complementarity_constraints!(model_container, big_m, target_timepoints, scenarios)
    add_lolmin_complementarity_constraints!(model_container, big_m, target_timepoints, scenarios)
    add_pmin_complementarity_constraints!(model_container, big_m, target_timepoints, scenarios, network)
    add_pmax_complementarity_constraints!(model_container, big_m, target_timepoints, scenarios, network)
end

"""
    complementarity constraint linked to @ref PSCOPF.global_capping_constraints!
"""
function add_emin_complementarity_constraints!(model_container::TSOBilevelModel,
                                            big_m, target_timepoints, scenarios)
    tso_limitable_model = model_container.upper.limitable_model
    market_limitable_model = model_container.lower.limitable_model
    kkt_model = model_container.kkt_model

    for ts in target_timepoints
        for s in scenarios
            kkt_var = kkt_model.capping_duals[ts,s]
            cstr_expr = market_limitable_model.p_global_capping[ts,s] - tso_limitable_model.p_global_capping[ts,s]
            b_indicator = kkt_model.capping_indicators[ts,s]
            ub_cstr = compute_ub(cstr_expr, big_m)
            formulate_complementarity_constraints!(kkt_model.model, kkt_var, cstr_expr, b_indicator, big_m, ub_cstr)
        end
    end
    return model_container
end

"""
    complementarity constraint linked to @ref PSCOPF.global_lol_constraints!
"""
function add_lolmin_complementarity_constraints!(model_container::TSOBilevelModel,
                                                big_m, target_timepoints, scenarios)
    tso_lol_model = model_container.upper.lol_model
    market_lol_model = model_container.lower.lol_model
    kkt_model = model_container.kkt_model

    for ts in target_timepoints
        for s in scenarios
            kkt_var = kkt_model.loss_of_load_duals[ts,s]
            cstr_expr = market_lol_model.p_global_loss_of_load[ts,s] - tso_lol_model.p_global_loss_of_load[ts,s]
            b_indicator = kkt_model.loss_of_load_indicators[ts,s]
            ub_cstr = compute_ub(cstr_expr, big_m)
            formulate_complementarity_constraints!(kkt_model.model, kkt_var, cstr_expr, b_indicator, big_m, ub_cstr)
        end
    end
    return model_container
end

"""
    complementarity constraint linked to the p_min part of @ref PSCOPF.respect_impositions_constraints!
"""
function add_pmin_complementarity_constraints!(model_container::TSOBilevelModel,
                                                big_m, target_timepoints, scenarios, network)
    tso_pilotable_model = model_container.upper.pilotable_model
    market_pilotable_model = model_container.lower.pilotable_model
    kkt_model = model_container.kkt_model

    for ts in target_timepoints
        for s in scenarios
            for pilotable_gen in Networks.get_generators_of_type(network, Networks.PILOTABLE)
                gen_id = Networks.get_id(pilotable_gen)

                kkt_var = kkt_model.pmin_duals[gen_id,ts,s]
                cstr_expr = market_pilotable_model.p_injected[gen_id,ts,s] - tso_pilotable_model.p_imposition_min[gen_id,ts,s]
                b_indicator = kkt_model.pmin_indicators[gen_id,ts,s]
                ub_cstr = compute_ub(cstr_expr, big_m) #or get_p_max(pilotable_gen)
                formulate_complementarity_constraints!(kkt_model.model, kkt_var, cstr_expr, b_indicator, big_m, ub_cstr)
            end
        end
    end
    return model_container
end

"""
    complementarity constraint linked to the p_max part of @ref PSCOPF.respect_impositions_constraints!
"""
function add_pmax_complementarity_constraints!(model_container::TSOBilevelModel,
                                                big_m, target_timepoints, scenarios, network)
    tso_pilotable_model = model_container.upper.pilotable_model
    market_pilotable_model = model_container.lower.pilotable_model
    kkt_model = model_container.kkt_model

    for ts in target_timepoints
        for s in scenarios
            for pilotable_gen in Networks.get_generators_of_type(network, Networks.PILOTABLE)
                gen_id = Networks.get_id(pilotable_gen)

                kkt_var = kkt_model.pmax_duals[gen_id,ts,s]
                cstr_expr = tso_pilotable_model.p_imposition_max[gen_id,ts,s] - market_pilotable_model.p_injected[gen_id,ts,s]
                b_indicator = kkt_model.pmax_indicators[gen_id,ts,s]
                ub_cstr = compute_ub(cstr_expr, big_m) #or get_p_max(pilotable_gen)
                formulate_complementarity_constraints!(kkt_model.model, kkt_var, cstr_expr, b_indicator, big_m, ub_cstr)
            end
        end
    end
    return model_container
end

#############################
# Utils
#############################
function coeffxsum(vars_dict::AbstractDict{T,V}, coeff::Float64
                )::GenericAffExpr{Float64, VariableRef} where T <: Tuple where V <: AbstractVariableRef
    terms = [var_l=>coeff for (_, var_l) in vars_dict]
    result = GenericAffExpr{Float64, VariableRef}(0., terms)
    return result
end

##########################################################
#        TSOBilevel
##########################################################
function tso_bilevel(network::Networks.Network,
                    target_timepoints::Vector{Dates.DateTime},
                    generators_initial_state::SortedDict{String,GeneratorState},
                    scenarios::Vector{String},
                    uncertainties_at_ech::UncertaintiesAtEch,
                    firmness::Firmness,
                    preceding_market_schedule::Schedule,
                    preceding_tso_schedule::Schedule,
                    # gratis_starts::Set{Tuple{String,Dates.DateTime}},
                    configs::TSOBilevelConfigs
                    )

    @assert(configs.big_m >= configs.MARKET_LOL_PENALTY)
    @assert(configs.big_m >= configs.MARKET_CAPPING_COST)
    @assert(all( configs.big_m >= Networks.get_prop_cost(gen)
                for gen in Networks.get_generators_of_type(network, Networks.PILOTABLE) ))

    if is_market(configs.REF_SCHEDULE_TYPE_IN_TSO)
        reference_schedule = preceding_market_schedule
    elseif is_tso(configs.REF_SCHEDULE_TYPE_IN_TSO)
        reference_schedule = preceding_tso_schedule
    else
        throw( error("Invalid REF_SCHEDULE_TYPE_IN_TSO config.") )
    end

    bimodel_container_l = TSOBilevelModel()

    create_tso_vars!(bimodel_container_l.upper,
                    network, target_timepoints, scenarios,
                    preceding_market_schedule)
    create_market_vars!(bimodel_container_l.lower,
                        network, target_timepoints, scenarios,
                        preceding_market_schedule)

    #this is the expression no objective is added to the jump model
    create_market_objectives!(bimodel_container_l.lower, network,
                            configs.MARKET_CAPPING_COST, configs.MARKET_LOL_PENALTY)

    #constraints may use upper and lower vars at the same time
    add_tso_constraints!(bimodel_container_l, target_timepoints, scenarios, network,
                        firmness, reference_schedule, generators_initial_state,
                        uncertainties_at_ech, configs)
    add_market_constraints!(bimodel_container_l, target_timepoints, scenarios, network, firmness, uncertainties_at_ech, configs)
    #kkt stationarity
    add_kkt_stationarity_constraints!(bimodel_container_l.kkt_model,
                                    target_timepoints, scenarios, network,
                                    configs.MARKET_CAPPING_COST, configs.MARKET_LOL_PENALTY)
    #kkt complementarity
    add_kkt_complementarity_constraints!(bimodel_container_l, configs.big_m, target_timepoints, scenarios, network)

    create_tso_objectives!(bimodel_container_l,
                        target_timepoints, scenarios, network,
                        reference_schedule, #reference to see which units are currently on
                        configs.TSO_CAPPING_COST, configs.TSO_LOL_PENALTY,
                        configs.TSO_LIMIT_PENALTY,
                        configs.TSO_PILOTABLE_BOUNDING_COST, configs.USE_UNITS_PROP_COST_AS_TSO_BOUNDING_COST)

    launch_solve!(bimodel_container_l, configs)

    return bimodel_container_l
end

function launch_solve!(bimodel_container::TSOBilevelModel, configs::TSOBilevelConfigs)
    # first_iter = true
    # while(!isempty(to_add) || first_iter)
    #     first_iter = false
    #     solve!(bimodel_container, configs.problem_name, configs.out_path)

    #     to_add = verify_rso()
    #     add_constraints(to_add)
    # end

    #the upper problem's objective model holds the two objective expressions and the deltas expression
    if configs.CONSIDER_DELTAS
        solve_2steps_deltas!(bimodel_container.upper, configs)
    else
        obj = bimodel_container.upper.objective_model.full_obj_2
        @objective(get_model(bimodel_container), Min, obj)
        solve!(bimodel_container, configs.problem_name, configs.out_path)
    end
    @info("Lower Objective Value : $(value(bimodel_container.lower.objective_model.full_obj))")
end

max_stars_repo_name: klorel/pscopf

max_issues_repo_name: klorel/pscopf

max_forks_repo_name: klorel/pscopf


text:
using SIMD

"""
    chipsort_small!(data, Val(V), Val(J))

Sort a small array using a vectorized sorting network and bitonic merge networks. The initial sort is over `J` SIMD vectors of size `V`. The array should typically fit inside register memory, for instance 64 Int32 values on an AVX2 machine.

# Examples
```jldoctest
julia> chipsort_small!(Int32[1:16...] .* Int32(1729) .%0x100, Val(4), Val(4))'
1×16 LinearAlgebra.Adjoint{Int32,Array{Int32,1}}:
 4  8  12  16  67  71  75  79  130  134  138  142  193  197  201  205
```
"""
chipsort_small!(data::AbstractVector{T}, ::Val{V}, ::Val{J}) where {T,V,J} = sort_block!(data, Val(V), Val(J))

@generated function sort_block!(data::AbstractVector{T}, ::Val{V}, ::Val{J}) where {T,V,J}
    ex = [Expr(:meta, :inline)]

    for j in 1:J
        v = Symbol("input_", j)
        push!(ex, :($v = vload(Vec{V,T}, pointer(data, 1+($j-1)*V))))
    end
    vecs=[Symbol("input_", j) for j in 1:J]
    append!(ex, [
        :(output = merge_vecs(transpose_vecs(sort_net($(vecs...))...)...)),
        :(vstore(output, pointer(data, 1))),
        :(return data)
    ])
    quote $(ex...) end
end

function sort_blocks!(data::AbstractVector{T}, ::Val{L}, ::Val{N}) where {L,N,T}
    block_size = N*L
    num_blocks = div(length(data), block_size)

    for m in 1:num_blocks
        sort_block!((@view data[1 + (m-1)*block_size:m*block_size]), Val(N), Val(L))
    end
    data
end

"""Sort blocks from the input array, optionally transposing to generate sorted sequences of size V (default), or further
merging those into sequences of size V*L."""
function sort_vecs!(input::AbstractVector{T}, ::Val{J}, ::Val{V}, ::Val{Transpose}=Val(true), ::Val{Merge}=Val(false)) where {V,J,T,Transpose,Merge}
    block_size = V*J
    num_blocks = div(length(input), block_size)

    for m in 1:num_blocks
        block = ntuple(j->vload(Vec{V, T}, input, 1 + (m-1)*block_size + (j-1)*V), Val(J))
        sorted_vecs = if Merge
            merge_vecs(transpose_vecs(sort_net(block...)...)...)
        elseif Transpose
            transpose_vecs(sort_net(block...)...)
        else
            sort_net(block...)
        end

        if Merge
            vstore(sorted_vecs, input, 1 + (m-1)*(V*J))
        elseif Transpose
            for v in 1:V
                vstore(sorted_vecs[v], input, 1 + (m-1)*(V*J) + (v-1)*J)
            end
        else
            for j in 1:J
                vstore(sorted_vecs[j], input, 1 + (m-1)*(V*J) + (j-1)*V)
            end
        end
    end
    input
end


text:
using Base: CFunction, Float32, Float64

mutable struct PathBuilder
    controlsLevel::Array{String}
end

function pushLabel!(builder::PathBuilder, label::String)
    push!(builder.controlsLevel, label)
end

function popLabel!(builder::PathBuilder)
    deleteat!(builder.controlsLevel, lastindex(builder.controlsLevel))
end

function buildPath(builder::PathBuilder, label::String)
    path = join(builder.controlsLevel, "/")
    res = "/$path/$label"
    for c in [' ', '#', '*', ',', '?', '[', ']', '{', '}', '(', ')']
        res = replace(res, c => '_')
    end
    res
end

mutable struct GlueFns
    openTabBox::CFunction
    openHorizontalBox::CFunction
    openVerticalBox::CFunction
    closeBox::CFunction
    addButton::CFunction
    addCheckButton::CFunction
    addVerticalSlider::CFunction
    addHorizontalSlider::CFunction
    addNumEntry::CFunction
    addHorizontalBargraph::CFunction
    addVerticalBargraph::CFunction
    # addSoundfile::CFunction
    declare::CFunction
end

struct UIRange{T}
    init::T
    min::T
    max::T
    step::T
end

mutable struct UIGlue{T}
    uiInterface::Ptr{Cvoid}
    openTabBox::Ptr{Cvoid}
    openHorizontalBox::Ptr{Cvoid}
    openVerticalBox::Ptr{Cvoid}
    closeBox::Ptr{Cvoid}
    addButton::Ptr{Cvoid}
    addCheckButton::Ptr{Cvoid}
    addVerticalSlider::Ptr{Cvoid}
    addHorizontalSlider::Ptr{Cvoid}
    addNumEntry::Ptr{Cvoid}
    addHorizontalBargraph::Ptr{Cvoid}
    addVerticalBargraph::Ptr{Cvoid}
    addSoundfile::Ptr{Cvoid}
    declare::Ptr{Cvoid}

    paths::Dict{String,Ptr{T}}
    ranges::Dict{String,UIRange{T}}
    pathBuilder::PathBuilder

    gluefns::GlueFns

    function UIGlue{T}() where {T}
        uglue = new{T}()
        uglue.paths = Dict{String, Ptr{T}}()
        uglue.pathBuilder = PathBuilder([])
        uglue.ranges = Dict{String, UIRange{T}}()
        initGlue!(uglue)
    end
end

function initGlue!(uglue::UIGlue{T}) where {T}

    function _openTabBox(ui, label)::Cvoid
        pushLabel!(uglue.pathBuilder, unsafe_string(label))
        nothing
    end
    openTabBox = @cfunction($_openTabBox, Cvoid, (Ptr{Cvoid}, Cstring))

    function _openHorizontalBox(ui, label)::Cvoid
        pushLabel!(uglue.pathBuilder, unsafe_string(label))
        nothing
    end
    openHorizontalBox = @cfunction($_openHorizontalBox, Cvoid, (Ptr{Cvoid}, Cstring))

    function _openVerticalBox(ui, label)::Cvoid
        pushLabel!(uglue.pathBuilder, unsafe_string(label))
        nothing
    end
    openVerticalBox = @cfunction($_openVerticalBox, Cvoid, (Ptr{Cvoid}, Cstring))

    function _closeBox(ui)::Cvoid
        popLabel!(uglue.pathBuilder)
        nothing
    end
    closeBox = @cfunction($_closeBox, Cvoid, (Ptr{Cvoid},))

    function _addButton(ui, label, zone)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        nothing
    end
    addButton = @cfunction(
        $_addButton,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}))

    function _addCheckButton(ui, label, zone)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        nothing
    end
    addCheckButton = @cfunction(
        $_addCheckButton,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}))

    function _addVerticalSlider(ui, label, zone, init, fmin, fmax, step)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        uglue.ranges[path] = UIRange(init, fmin, fmax, step)
        nothing
    end
    addVerticalSlider = @cfunction(
        $_addVerticalSlider,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}, T, T, T, T))

    function _addHorizontalSlider(ui, label, zone, init, fmin, fmax, step)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        uglue.ranges[path] = UIRange(init, fmin, fmax, step)
        nothing
    end
    addHorizontalSlider = @cfunction(
        $_addHorizontalSlider,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}, T, T, T, T))

    function _addNumEntry(ui, label, zone, init, fmin, fmax, step)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        uglue.ranges[path] = UIRange(init, fmin, fmax, step)
        nothing
    end
    addNumEntry = @cfunction(
        $_addNumEntry,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}, T, T, T, T))

    function _addHorizontalBargraph(ui, label, zone, fmin, fmax)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        nothing
    end
    addHorizontalBargraph = @cfunction(
        $_addHorizontalBargraph,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}, T, T))

    function _addVerticalBargraph(ui, label, zone, fmin, fmax)::Cvoid
        path = buildPath(uglue.pathBuilder, unsafe_string(label))
        uglue.paths[path] = zone
        nothing
    end
    addVerticalBargraph = @cfunction(
        $_addVerticalBargraph,
        Cvoid, (Ptr{Cvoid}, Cstring, Ptr{T}, T, T))

    function _declare(ui, zone, key, value)::Cvoid
        nothing
    end
    declare = @cfunction(
        $_declare,
        Cvoid, (Ptr{Cvoid}, Ptr{T}, Cstring, Cstring))

# // -- soundfiles
# typedef void (* addSoundfileFun) (void* ui_interface, const char* label, const char* url, struct Soundfile** sf_zone);

    uglue.gluefns = GlueFns(
        openTabBox,
        openHorizontalBox,
        openVerticalBox,
        closeBox,
        addButton,
        addCheckButton,
        addVerticalSlider,
        addHorizontalSlider,
        addNumEntry,
        addHorizontalBargraph,
        addVerticalBargraph,
        declare,
    )
    uglue.openTabBox = uglue.gluefns.openTabBox.ptr
    uglue.openHorizontalBox = uglue.gluefns.openHorizontalBox.ptr
    uglue.openVerticalBox = uglue.gluefns.openVerticalBox.ptr
    uglue.closeBox = uglue.gluefns.closeBox.ptr
    uglue.addButton = uglue.gluefns.addButton.ptr
    uglue.addCheckButton = uglue.gluefns.addCheckButton.ptr
    uglue.addVerticalSlider = uglue.gluefns.addVerticalSlider.ptr
    uglue.addHorizontalSlider = uglue.gluefns.addHorizontalSlider.ptr
    uglue.addNumEntry = uglue.gluefns.addNumEntry.ptr
    uglue.addHorizontalBargraph = uglue.gluefns.addHorizontalBargraph.ptr
    uglue.addVerticalBargraph = uglue.gluefns.addVerticalBargraph.ptr
    uglue.declare = uglue.gluefns.declare.ptr

    uglue
end


text:
using Test, QMPS, Yao, Random, Statistics
using MSQR

# Un-comment if testing GPU acceleration compatibility.
# using CuYao

@testset "SWAPtest(SWAPtrain.jl)" begin
    seedNum = 1234
    m = 1
    n = 3
    Random.seed!(seedNum)
    cgen = dispatch!(DCbuilder(n,4).circuit, :random)
    reg1 = rand_state(n)
    reg2 = rand_state(n)
    reg3 = zero_state(n) |> cgen
    olp12m1 = SWAPtest(reg1, reg2, nMeasure=m)
    olp12m2 = SWAPtest(join(zero_state(1), reg1, reg2), nMeasure=m)
    olp23m1 = SWAPtest(reg2, cgen, nMeasure=m)
    olp23m2 = SWAPtest(cgen, regAll = join(zero_state(n+1), reg2), nMeasure=m)

    cST = chain(7, put(7, 7=>H),
                   put(7, 7=>shift(0)),
                   chain(7, control(7, 7, (6,3)=>SWAP),
                            control(7, 7, (5,2)=>SWAP),
                            control(7, 7, (4,1)=>SWAP)),
                   put(7, 7=>H))
    op = put(2n+1, 2n+1=>Z)
    reg12t = repeat(join(zero_state(1), reg1, reg2), m) |> cST
    reg23t = repeat(join(zero_state(1), reg3, reg2), m) |> cST
    olp12t = expect(op, reg12t) |> mean |> real
    olp23t = expect(op, reg23t) |> mean |> real

    @test olp12m1.reg ≈ olp12m2.reg
    @test olp23m1.reg ≈ olp23m2.reg
    @test olp12m1.reg ≈ reg12t
    @test olp23m1.reg ≈ reg23t
    @test olp12m2.witnessOp == op
    @test olp12m2.circuit == cST
    @test isapprox(olp12m1.overlap, olp12m2.overlap, atol=10e-12)
    @test isapprox(olp23m1.overlap, olp23m2.overlap, atol=10e-12)
    @test isapprox(olp12m1.overlap, olp12t, atol=10e-12)
    @test isapprox(olp23m1.overlap, olp23t, atol=10e-12)
end

@testset "MPSSwapTest.jl" begin
    seedNum = 1234
    # MScircuit(nBitT::Int64, vBit::Int64, rBit::Int64, MPSblocks::Array{CompositeBlock,1}; ϕ::Float64=0.0)
    ## MPSblcoks = MPSC("CS", 4, 1, 1).mpsBlocks
    n0 = 4
    v0 = 1
    r0 = 1
    Random.seed!(seedNum)
    mps0cs = MPSC("CS", n0, v0, r0)
    csMblocks = mps0cs.mpsBlocks
    CSc = MScircuit(n0, v0, r0, mps0cs.mpsBlocks)
    CSct = chain(7, put(7, 7=>H),
                   chain(7, subroutine(7, csMblocks[1], 5:6), control(7, 7, (4,6)=>SWAP), Measure(7, locs=(6), resetto=0)),
                   chain(7, subroutine(7, csMblocks[2], 5:6), control(7, 7, (3,6)=>SWAP), Measure(7, locs=(6), resetto=0)),
                   chain(7, subroutine(7, csMblocks[3], 5:6), control(7, 7, (2,6)=>SWAP), Measure(7, locs=(6), resetto=0)),
                   control(7, 7, (1,5)=>SWAP),
                   put(7, 7=>H))
    opx(n) = chain(n, [put(n, i=>X) for i=1:n])
    opy(n) = chain(n, [put(n, i=>Y) for i=1:n])
    opz(n) = chain(n, [put(n, i=>Z) for i=1:n])
    function cr_test(reg::ArrayReg, cr::CompositeBlock, crt::CompositeBlock)
        n = nqubits(reg)
        Random.seed!(seedNum)
        reg1x = expect(opx(n), copy(reg) |> cr) |> mean |> real
        Random.seed!(seedNum)
        reg1y = expect(opy(n), copy(reg) |> cr) |> mean |> real
        Random.seed!(seedNum)
        reg1z = expect(opz(n), copy(reg) |> cr) |> mean |> real
        Random.seed!(seedNum)
        reg2x = expect(opx(n), copy(reg) |> crt) |> mean |> real
        Random.seed!(seedNum)
        reg2y = expect(opy(n), copy(reg) |> crt) |> mean |> real
        Random.seed!(seedNum)
        reg2z = expect(opz(n), copy(reg) |> crt) |> mean |> real
        @test reg1x ≈ reg2x
        @test reg1y ≈ reg2y
        @test reg1z ≈ reg2z
    end
    cr_test(repeat(join(zero_state(3),rand_state(4)),5000), CSc, CSct)

    ## MPSblcoks = MPSC(("DC",2), 8, 2, 2).mpsBlocks
    n1 = 8
    v1 = 2
    r1 = 2
    d1 = 2
    Random.seed!(seedNum)
    mps0dc = MPSC(("DC",d1), n1, v1, r1)
    dcMblocks = mps0dc.mpsBlocks
    DCc = MScircuit(n1, v1, r1, mps0dc.mpsBlocks)
    DCct = chain(13, put(13, 13=>H),
                     chain(13, subroutine(13, dcMblocks[1], 9:12), control(13, 13, (12, 8)=>SWAP), control(13, 13, (11, 7)=>SWAP), Measure(13, locs=(11,12), resetto=0)),
                     chain(13, subroutine(13, dcMblocks[2], 9:12), control(13, 13, (12, 6)=>SWAP), control(13, 13, (11, 5)=>SWAP), Measure(13, locs=(11,12), resetto=0)),
                     chain(13, subroutine(13, dcMblocks[3], 9:12), control(13, 13, (12, 4)=>SWAP), control(13, 13, (11, 3)=>SWAP), Measure(13, locs=(11,12), resetto=0)),
                     chain(13, control(13, 13, (10, 2)=>SWAP), control(13, 13, (9, 1)=>SWAP)),
                     put(13, 13=>H))
    cr_test(repeat(join(zero_state(5),rand_state(8)),5000), DCc, DCct)

    # MStest(regT::ArrayReg, MSCircuit::ChainBlock; nMeasure::Int64=1, useCuYao::Bool=CUDA_ON)
    # MStest(MSCircuit::ChainBlock; regAll::ArrayReg)
    # MSTtest(regT::ArrayReg, MSCircuit::ChainBlock, cExtend::ChainBlock; nMeasure::Int64=1, useCuYao::Bool=CUDA_ON)
    ## MPSblcoks = MPSC("CS", 4, 1, 1).mpsBlocks
    m = 500
    n = 4
    v = 1
    r = 1
    Random.seed!(seedNum)
    regT = rand_state(n)
    mps1 = MPSC("CS", n, v, r)
    c1 = MScircuit(n, v, r, mps1.mpsBlocks)
    Random.seed!(seedNum)
    t1_0_m1 = MStest(regT, c1, nMeasure=m)
    Random.seed!(seedNum)
    t1_0_m2 = MStest(c1, regAll=repeat(join(zero_state(nqubits(c1)-n),regT),m))
    Random.seed!(seedNum)
    t1 = MSTtest(regT, c1, mps1.cExtend, nMeasure=m)
    @test t1_0_m1.reg ≈ t1_0_m2.reg
    @test t1_0_m1.witnessOp == t1_0_m2.witnessOp
    @test isapprox(t1_0_m1.overlap, t1_0_m2.overlap, atol=10e-12)
    @test isapprox(t1_0_m1.overlap, t1.Aoverlap, atol=10e-12)
    @test isapprox(t1.Eoverlap, t1.Aoverlap, atol=1.001abs(t1.error*t1.Eoverlap))
    Ext_overlap = SWAPtest(zero_state(n)|>mps1.cExtend, regT, nMeasure=m).overlap
    @test isapprox.(t1.Eoverlap, Ext_overlap, atol=1.001abs(t1.error*t1.Eoverlap))

    ## MPSblcoks = MPSC(("DC",2), 6, 2, 2).mpsBlocks
    n = 4
    v = 2
    r = 1
    d = 3
    Random.seed!(seedNum)
    regT = rand_state(n)
    mps2 = MPSC(("DC",d), n, v, r)
    c2 = MScircuit(n, v, r, mps2.mpsBlocks)
    p = MSCpar(c2)
    @test (p.nBitA, p.nBitT, p.vBit, p.rBit, p.depth) == (nqubits(c2), n, v, r, d)
    Random.seed!(seedNum)
    t2_0_m1 = MStest(regT, c2, nMeasure=m)
    Random.seed!(seedNum)
    t2_0_m2 = MStest(c2, regAll=repeat(join(zero_state(p.nBitA-n),regT),m))
    Random.seed!(seedNum)
    t2 = MSTtest(regT, c2, mps2.cExtend, nMeasure=m)
    @test t2_0_m1.reg ≈ t2_0_m2.reg
    @test t2_0_m1.witnessOp == t2_0_m2.witnessOp
    @test isapprox(t2_0_m1.overlap, t2_0_m2.overlap, atol=10e-12)
    @test isapprox(t2_0_m1.overlap, t2.Aoverlap, atol=10e-12)
    @test isapprox(t2.Eoverlap, t2.Aoverlap, atol=1.001abs(t2.error*t2.Eoverlap))
    Ext_overlap = SWAPtest(zero_state(n)|>mps2.cExtend, regT, nMeasure=m).overlap
    @test isapprox.(t2.Eoverlap, Ext_overlap, atol=1.001abs(t2.error*t2.Eoverlap))
end

@testset "MSQRtrain!(MSQRtrain.jl) + SWAPtrain!(SWAPtrain.jl)" begin
    seedNum = 1234
    n = 3
    v = 1
    r = 1
    d = 2
    Random.seed!(seedNum)
    regT = rand_state(n)
    Random.seed!(seedNum)
    mps = MPSC(("DC",d),n,v,r)
    c = MScircuit(n, v, r, mps.mpsBlocks)
    md1 = ("ADAM", 0.05)
    md2 = ("default", 0.2)
    m = 150
    ## Testing function exceptions.
    @test_throws ErrorException GDescent("default")

    mps_m1q = deepcopy(mps)
    c_m1q = MScircuit(n, v, r, mps_m1q.mpsBlocks)
    mps_s1q = deepcopy(mps)
    c_s1q = mps_s1q.cExtend
    mps_m2q = deepcopy(mps)
    c_m2q = MScircuit(n, v, r, mps_m2q.mpsBlocks)
    mps_s2q = deepcopy(mps)
    c_s2q = mps_s2q.cExtend
    mps_m1n = deepcopy(mps)
    c_m1n = MScircuit(n, v, r, mps_m1n.mpsBlocks)
    mps_s1n = deepcopy(mps)
    c_s1n = mps_s1n.cExtend
    mps_m2n = deepcopy(mps)
    c_m2n = MScircuit(n, v, r, mps_m2n.mpsBlocks)
    mps_s2n = deepcopy(mps)
    c_s2n = mps_s2n.cExtend

    showON = false
    Random.seed!(seedNum)
    mres1q = MSQRtrain!(regT, c_m1q, 40, nMeasure=m, GDmethod = md1, show=true) # Test show option.
    Random.seed!(seedNum)
    mres2q = MSQRtrain!(regT, c_m2q, 50, nMeasure=m, GDmethod = md2, show=showON)
    Random.seed!(seedNum)
    sres1q = SWAPtrain!(regT, c_s1q, 40, nMeasure=m, GDmethod = md1, show=true) # Test show option.
    Random.seed!(seedNum)
    sres2q = SWAPtrain!(regT, c_s2q, 50, nMeasure=m, GDmethod = md2, show=showON)
    Random.seed!(seedNum)
    mres1n = MSQRtrain!(regT, c_m1n, 40, nMeasure=m, GDmethod = md1, Gmethod=("Ndiff", 0.05), show=showON)
    Random.seed!(seedNum)
    mres2n = MSQRtrain!(regT, c_m2n, 50, nMeasure=m, GDmethod = md2, Gmethod=("Ndiff", 0.05), show=showON)
    Random.seed!(seedNum)
    sres1n = SWAPtrain!(regT, c_s1n, 40, nMeasure=m, GDmethod = md1, Gmethod=("Ndiff", 0.05), show=showON)
    Random.seed!(seedNum)
    sres2n = SWAPtrain!(regT, c_s2n, 50, nMeasure=m, GDmethod = md2, Gmethod=("Ndiff", 0.05), show=showON)
    ## Test different methods / optional arguments of the functions.
    c_m1n_2 = MScircuit(n, v, r, deepcopy(mps).mpsBlocks)
    Random.seed!(seedNum)
    a = MSQRtrain!(regT, c_m1n_2, 150, nMeasure=200, GDmethod = ("ADAM", 0.05, (0.9, 0.999)), Gmethod="Ndiff", show=showON)
    mres1n_2 = MSQRtrain!(regT, c_m1n_2, :auto, nMeasure=m, GDmethod = md1, show=showON, ConvTh=(5e-3,1e-2))
    c_s1n_2 = deepcopy(mps).cExtend
    Random.seed!(seedNum)
    b = SWAPtrain!(regT, c_s1n_2, 20, nMeasure=m, GDmethod = md1, show=showON)
    sres1n_2 = SWAPtrain!(regT, c_s1n_2, :auto, nMeasure=m, GDmethod = "ADAM", show=showON, ConvTh=(5e-3,1e-2))

    # If all the trainings converge in the end.
    cuti = 9
    tol = 0.06
    mres1qConv = mres1q[end-cuti:end]    # MSQR + ADAM + Qdiff
    mres2qConv = mres2q[end-cuti:end]    # MSQR + CONS + Qdiff
    sres1qConv = sres1q[end-cuti:end]    # SWAP + ADAM + Qdiff
    sres2qConv = sres2q[end-cuti:end]    # SWAP + CONS + Qdiff
    mres1nConv = mres1n[end-cuti:end]    # MSQR + ADAM + Ndiff
    mres2nConv = mres2n[end-cuti:end]    # MSQR + CONS + Ndiff
    sres1nConv = sres1n[end-cuti:end]    # SWAP + ADAM + Ndiff
    sres2nConv = sres2n[end-cuti:end]    # SWAP + CONS + Ndiff
    # @show mres1n_2Conv = mres1n_2[end-cuti:end]
    # @show sres1n_2Conv = sres1n_2[end-cuti:end]
    mres1n_2Conv = mres1n_2[end-cuti:end]
    sres1n_2Conv = sres1n_2[end-cuti:end]

    mres1qMean = mres1qConv |> mean
    mres2qMean = mres2qConv |> mean
    sres1qMean = sres1qConv |> mean
    sres2qMean = sres2qConv |> mean
    mres1nMean = mres1nConv |> mean
    mres2nMean = mres2nConv |> mean
    sres1nMean = sres1nConv |> mean
    sres2nMean = sres2nConv |> mean
    mres1n_2Mean = mres1n_2Conv |> mean
    sres1n_2Mean = sres1n_2Conv |> mean

    @test std(mres1qConv) ≤ tol*mres1qMean
    @test std(mres2qConv) ≤ tol*mres2qMean
    @test std(sres1qConv) ≤ tol*sres1qMean
    @test std(sres2qConv) ≤ tol*sres2qMean
    @test std(mres1nConv) ≤ tol*mres1nMean
    @test std(mres2nConv) ≤ tol*mres2nMean
    @test std(sres1nConv) ≤ tol*sres1nMean
    @test std(sres2nConv) ≤ tol*sres2nMean

    # If ADAM and default SGD converge to the approximately same value.
    ## (Including Ndiff and Qdiff situations)
    tol2 = 0.06
    @test isapprox(mres1qMean, mres2qMean, atol=tol2*max( mres1qMean, mres2qMean ))
    @test isapprox(sres1qMean, sres2qMean, atol=tol2*max( sres1qMean, sres2qMean ))
    @test isapprox(mres1nMean, mres2nMean, atol=tol2*max( mres1nMean, mres2nMean ))
    @test isapprox(sres1nMean, sres2nMean, atol=tol2*max( sres1nMean, sres2nMean ))

    # If MSQRtrain and SWAPtrain converge to the approximately same value.
    ## (Including Ndiff and Qdiff situations)
    tol3 = 0.06
    @test isapprox(mres1qMean, sres1qMean, atol=tol3*max( mres1qMean, sres1qMean ))
    @test isapprox(mres2qMean, sres2qMean, atol=tol3*max( mres2qMean, sres2qMean ))
    @test isapprox(mres1nMean, sres1nMean, atol=tol3*max( mres1nMean, sres1nMean ))
    @test isapprox(mres2nMean, sres2nMean, atol=tol3*max( mres2nMean, sres2nMean ))

    # If the auto-train actually converges to theoretical maximum.
    @test isapprox(mres1n_2Mean, 1.0, atol = 1.5e-2)
    @test isapprox(sres1n_2Mean, 1.0, atol = 1.5e-2)

    # If the training curves monotonically rise.
    function trendCompr(res::Array{Float64, 1})
        midv = middle(res)
        @test res[end] > midv > res[1]
        max = (2midv-res[1])
        @test res[end] / max > 0.9
    end
    trendCompr(mres1q)
    trendCompr(sres1q)
    trendCompr(mres2q)
    trendCompr(sres2q)
    trendCompr(mres1n)
    trendCompr(sres1n)
    trendCompr(mres2n)
    trendCompr(sres2n)

    # If the overlaps from MSQR have enough accuracy.
    function oltest(regT::ArrayReg, cGen::ChainBlock, nM::Int64, ol0::Float64)
        nq = nqubits(cGen)
        regG = zero_state(nq, nbatch=nM) |> cGen
        ol = ((regT.state'*regG.state)[1] |> abs)^2
        @test isapprox(ol, ol0, atol=0.005)
    end
    oltest(regT, mps_m1q.cExtend, m, mres1q[end])
    oltest(regT, c_s1q, m, sres1q[end])
    oltest(regT, mps_m2q.cExtend, m, mres2q[end])
    oltest(regT, c_s2q, m, sres2q[end])
    oltest(regT, mps_m1n.cExtend, m, mres1n[end])
    oltest(regT, c_s1n, m, sres1n[end])
    oltest(regT, mps_m2n.cExtend, m, mres2n[end])
    oltest(regT, c_s2n, m, sres2n[end])
end


text:
abstract type SparseDiffToolsColoringAlgorithm <: ArrayInterface.ColoringAlgorithm end
struct GreedyD1Color <: SparseDiffToolsColoringAlgorithm end
struct BacktrackingColor <: SparseDiffToolsColoringAlgorithm end
struct ContractionColor <: SparseDiffToolsColoringAlgorithm end
struct GreedyStar1Color <: SparseDiffToolsColoringAlgorithm end
struct GreedyStar2Color <: SparseDiffToolsColoringAlgorithm end
struct AcyclicColoring <: SparseDiffToolsColoringAlgorithm end

"""
    matrix_colors(A, alg::ColoringAlgorithm = GreedyD1Color())

Return the colorvec vector for the matrix A using the chosen coloring
algorithm. If a known analytical solution exists, that is used instead.
The coloring defaults to a greedy distance-1 coloring.

Note that if A isa SparseMatrixCSC, the sparsity pattern is defined by structural nonzeroes,
ie includes explicitly stored zeros.
"""
function ArrayInterface.matrix_colors(A::AbstractMatrix, alg::SparseDiffToolsColoringAlgorithm = GreedyD1Color(); partition_by_rows::Bool = false)
    _A = A isa SparseMatrixCSC ? A : sparse(A) # Avoid the copy
    A_graph = matrix2graph(_A, partition_by_rows)
    return color_graph(A_graph, alg)
end


text:
# ------------------------------------------------------------------------------------------
# # Conditionals
#
# #### with the `if` keyword
# In Julia, the syntax
#
# ```julia
# if *condition 1*
#     *option 1*
# elseif *condition 2*
#     *option 2*
# else
#     *option 3*
# end
# ```
#
# allows us to conditionally evaluate one of our options.
#

# For example, we might want to implement the FizzBuzz test: given a number, N, print "Fizz"
# if N is divisible by 3, "Buzz" if N is divisible by 5, and "FizzBuzz" if N is divisible by
# 3 and 5. Otherwise just print the number itself! Enter your choice for `N` here:
# ------------------------------------------------------------------------------------------

N =

if (N % 3 == 0) && (N % 5 == 0) # `&&` means "AND"; % computes the remainder after division
    println("FizzBuzz")
elseif N % 3 == 0
    println("Fizz")
elseif N % 5 == 0
    println("Buzz")
else
    println(N)
end

# ------------------------------------------------------------------------------------------
# #### with ternary operators
#
# For this last block, we could instead use the ternary operator with the syntax
#
# ```julia
# a ? b : c
# ```
#
# which equates to
#
# ```julia
# if a
#     b
# else
#     c
# end
# ```
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# Now let's say we want to return the larger of two numbers. Give `x` and `y` values here:
# ------------------------------------------------------------------------------------------

x =
y =

# ------------------------------------------------------------------------------------------
# Using the `if` and `else` keywords, we might write:
# ------------------------------------------------------------------------------------------

if x > y
    x
else
    y
end

# ------------------------------------------------------------------------------------------
# and as a ternary operator, the conditional looks like this:
# ------------------------------------------------------------------------------------------

(x > y) ? x : y

# ------------------------------------------------------------------------------------------
# #### with short-circuit evaluation
#
# We've already seen expressions with the syntax
# ```julia
# a && b
# ```
# to return true if both `a` and `b` are true. Of course, if `a` is false, Julia doesn't
# even need to know the value of `b` in order to determine that the overall result will be
# false. So Julia doesn't even need to check what `b` is; it can just "short-circuit" and
# immediately return `false`.  The second argument `b` might be a more complicated
# expression like a function call with a side-effect, in which case it won't even be called:
# ------------------------------------------------------------------------------------------

false && (println("hi"); true)

true && (println("hi"); true)

# ------------------------------------------------------------------------------------------
# On the other hand, if `a` is true, Julia knows it can just return the value of `b` as the
# overall expression. This means that `b` doesn't necessarily need evaluate to `true` or
# `false`!  `b` could even be an error:
# ------------------------------------------------------------------------------------------

(x > 0) && error("x cannot be greater than 0")

# ------------------------------------------------------------------------------------------
# Similarly, check out the `||` operator, which also uses short-circuit evaluation to
# perform the "or" operation.
# ------------------------------------------------------------------------------------------

true || println("hi")

# ------------------------------------------------------------------------------------------
# and
# ------------------------------------------------------------------------------------------

false || println("hi")

# ------------------------------------------------------------------------------------------
# ### Exercises
#
# #### 5.1
# Write a conditional statement that prints a number if the number is even and the string
# "odd" if the number is odd.
# ------------------------------------------------------------------------------------------

# ------------------------------------------------------------------------------------------
# #### 5.2
# Rewrite the code from 5.1 using a ternary operator.
# ------------------------------------------------------------------------------------------

max_stars_repo_name: grenkoca/JuliaTutorials

max_issues_repo_name: grenkoca/JuliaTutorials

max_forks_repo_name: grenkoca/JuliaTutorials


text:
"""
inspectQNITensor is a developer-level debugging tool
to look at internals or properties of QNITensors
"""
function inspectQNITensor(T::ITensor,is::QNIndexSet)
  #@show T.store.blockoffsets
  #@show T.store.data
  println("Block fluxes:")
  for b in nzblocks(T)
    @show flux(T,b)
  end
end
inspectQNITensor(T::ITensor,is::IndexSet) = nothing
inspectQNITensor(T::ITensor) = inspectQNITensor(T,inds(T))


text:
# This file is a part of ArraysOfArrays.jl, licensed under the MIT License (MIT).

__precompile__(true)

"""
    module ArraysOfArrays

Efficient storage and handling of nested arrays.

ArraysOfArrays provides two different types of nested arrays:
[`ArrayOfSimilarArrays`](@ref section_ArrayOfSimilarArrays) and
[`VectorOfArrays`](@ref section_VectorOfArrays).
"""
module ArraysOfArrays

using Adapt
using Requires
using Statistics
using UnsafeArrays
using ChainRulesCore

include("util.jl")
include("functions.jl")
include("array_of_similar_arrays.jl")
include("vector_of_arrays.jl")

function __init__()
    @require StaticArrays = "90137ffa-7385-5640-81b9-e52037218182" include("staticarrays_support.jl")
    @require StatsBase = "2913bbd2-ae8a-5f71-8c99-4fb6c76f3a91" include("statsbase_support.jl")
end

end # module