text:
const OPTIMIZER = ()->MadNLP.Optimizer(
    linear_solver=MadNLPLapackCPU,
    print_level=MadNLP.ERROR
)

@testset "MINLPTests" begin
    ###
    ### src/nlp tests.
    ###

    MINLPTests.test_nlp(
        OPTIMIZER,
        exclude = [
            "005_011",  # Uses the function `\`
            "006_010",  # User-defined function without Hessian (autodiff only provides 1st order)
        ],
        objective_tol = 1e-5,
        primal_tol = 1e-5,
        dual_tol = NaN,
    )

    ###
    ### src/nlp-cvx tests.
    ###

    MINLPTests.test_nlp_cvx(
        OPTIMIZER,
    )
end


text:
import NLsolve: OnceDifferentiable, NewtonCache, newton_
import LineSearches

struct NLsolveNewton{T, AT, FT, DT, CT, ST, LT} <: AbstractNewtonSolver{T,AT}
    x::AT
    f::AT
    J::Matrix{T}

    F!::FT
    DF::DT

    line_search::ST
    linear_solver::LT

    cache::CT
    config::Options{T}
    status::NonlinearSolverStatus{T}

    function NLsolveNewton(x::AbstractVector{T}, f::AbstractVector{T}, J::AbstractMatrix{T},
                    F!::FT, DF::DT, cache::CT, line_search::ST, linear_solver::LT, config = Options()) where {T,FT,DT,CT,ST,LT}

        status = NonlinearSolverStatus{T}(length(x))

        new{T,typeof(x),FT,DT,CT,ST,LT}(x, f, J, F!, DF, line_search, linear_solver, cache, config, status)
    end
end

function NLsolveNewton(x::AbstractVector{T}, f::AbstractVector{T}, F!::Function; J!::Union{Function,Nothing}=nothing, mode = :autodiff, diff_type = :forward) where {T}
    linear_solver = LinearSolver(x)

    df = linear_solver.A

    if J! === nothing || mode == :autodiff
        DF = OnceDifferentiable(F!, x, f, df; autodiff=diff_type, inplace=true)
    else
        DF = OnceDifferentiable(F!, J!, x, f, df; inplace=true)
    end

    NLsolveNewton(x, f, df, F!, DF, NewtonCache(DF), LineSearches.Static(), linear_solver)
end

function linsolve!(s::NLsolveNewton, x, A, b)
    factorize!(s.linear_solver, A)
    ldiv!(x, s.linear_solver, b)
end

function solve!(x, s::NLsolveNewton)
    res = newton_(s.DF, x, s.config.x_abstol, s.config.f_abstol, s.config.max_iterations, false, false, false,
                  s.line_search, (x, A, b) -> linsolve!(s, x, A, b), s.cache)

    copyto!(x, res.zero)

    s.status.i = res.iterations
    s.status.rfₐ = res.residual_norm

    return x
end

ext: jl


text:
#=
---------------------------------------------------------------------
This function advances the solution in time using the LSRK45 Method by Carpenter-Kennedy 1994.

Written by F.X. Giraldo on July 12, 2021
           Department of Applied Mathematics
           Naval Postgraduate School
           Monterey; CA 93943-5216
---------------------------------------------------------------------
=#

include("create_rhs.jl")
include("plot_solution.jl")

function ti_LSRK!(q0,qb,coord,M,Dwe,intma,periodicity,time,ntime,dt,gravity,Δ_NL,space_method,plot_movie,case,DFloat)

    #Initialize RK coefficients
    RKA = (DFloat(0),
           DFloat(-567301805773)  / DFloat(1357537059087),
           DFloat(-2404267990393) / DFloat(2016746695238),
           DFloat(-3550918686646) / DFloat(2091501179385),
           DFloat(-1275806237668) / DFloat(842570457699 ))

    RKB = (DFloat(1432997174477) / DFloat(9575080441755 ),
           DFloat(5161836677717) / DFloat(13612068292357),
           DFloat(1720146321549) / DFloat(2090206949498 ),
           DFloat(3134564353537) / DFloat(4481467310338 ),
           DFloat(2277821191437) / DFloat(14882151754819))

    RKC = (DFloat(0),
           DFloat(1432997174477) / DFloat(9575080441755),
           DFloat(2526269341429) / DFloat(6820363962896),
           DFloat(2006345519317) / DFloat(3224310063776),
           DFloat(2802321613138) / DFloat(2924317926251))

    Npoin=size(q0,1)
    dq=zeros(DFloat,Npoin,2)
    qp=copy(q0)
    stages=length(RKA)

    #Time Integration
    for itime=1:ntime
        time=time + dt
        #RK Stages
        for s = 1:stages
            #Create RHS Matrix
            #= -----------Students add your CREATE_RHS routine here-----------=#
            #= -----------Students add your CREATE_RHS routine here-----------=#
            #R=create_rhs(qp,qb,M,Dwe,intma,periodicity,gravity,Δ_NL,DFloat)
            #= -----------Students add your CREATE_RHS routine here-----------=#
            #= -----------Students add your CREATE_RHS routine here-----------=#

            #Solve System
            @inbounds for I=1:Npoin
                dq[I,:] = RKA[s]*dq[I,:] + dt*R[I,:]
                qp[I,:] = qp[I,:] + RKB[s]*dq[I,:]
            end
            if periodicity[Npoin] == periodicity[1]
                qp[Npoin,:]=qp[1,:] #periodicity
            end
        end #s
        if mod(itime,10) == 0
            println(" itime = ",itime," time = ",time," qmax = ",maximum(qp)," qmin = ",minimum(qp))
        end

        #Update Q
        q0 .= qp

        #Plot Solution
        if (plot_movie)
            plot_solution(q0,qb,coord,space_method,time,case)
        end

    end #itime

    return(q0,time)

end


text:
@testset "Tree decomposition tests" begin
    # A simple graph to test flow cutter on.
    N = 10
    G = lg.SimpleGraph(N)
    for i = 1:N
        for j = i+1:N
            lg.add_edge!(G, i, j)
        end
    end

    # Check if the treewidth and number of bags is correct and check that the tree
    # decomposition has no edges.
    tree_decomp = flow_cutter(G, 30; seed=42)
    @test tree_decomp[:treewidth] == 9
    @test tree_decomp[:num_bags] == 1
    @test length(tree_decomp[:edges]) == 0

    # A simple graph with two disconnected cliques to test flow cutter on.
    N = 10; n = 5
    G = lg.SimpleGraph(N+n)
    # Clique 1
    for i = 1:N-1
        for j = i+1:N
            lg.add_edge!(G, i, j)
        end
    end
    # Clique 2
    for i = N+1:N+n-1
        for j = i+1:N+n
            lg.add_edge!(G, i, j)
        end
    end

    # Check if the treewidth, number of bags and number of edges is correct.
    tree_decomp = flow_cutter(G, 30; seed=42)
    @test tree_decomp[:treewidth] == N-1
    @test tree_decomp[:num_bags] == 2
    @test tree_decomp[:num_vertices] == N + n
    @test length(tree_decomp[:edges]) == 1

    # Test constructing a tree decomposition from an elimination order
    # A square lattice graph to test the min fill heuristic and conversion of a tree
    # decomposition to an elimination order.
    N = 5
    G = LabeledGraph(N*N)
    for i = 1:N-1
        for j = 1:N
            add_edge!(G, i + N*(j-1), i + N*(j-1) + 1)
            add_edge!(G, i + N*(j-1), i + N*(j-1) + N)
        end
    end
    treewidth_upperbound, min_fill_order = min_fill(G)
    B, T = build_clique_tree(G, min_fill_order)
    @test length(B) == lg.nv(T)
    @test treewidth_upperbound == maximum(length.(B)) - 1
end


text:
"""
GenX: An Configurable Capacity Expansion Model
Copyright (C) 2021,  Massachusetts Institute of Technology
This program is free software; you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation; either version 2 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.
A complete copy of the GNU General Public License v2 (GPLv2) is available
in LICENSE.txt.  Users uncompressing this from an archive may not have
received this license file.  If not, see <http://www.gnu.org/licenses/>.
"""

@doc raw"""
	investment_discharge(EP::Model, inputs::Dict, MinCapReq::Int)
This function defines the expressions and constraints keeping track of total available power generation/discharge capacity across all resources as well as constraints on capacity retirements.
The total capacity of each resource is defined as the sum of the existing capacity plus the newly invested capacity minus any retired capacity. Note for storage resources, additional energy and charge power capacity decisions and constraints are defined in the storage module.
```math
\begin{aligned}
& \Delta^{total}_{y,z} =(\overline{\Delta_{y,z}}+\Omega_{y,z}-\Delta_{y,z}) \forall y \in \mathcal{G}, z \in \mathcal{Z}
\end{aligned}
```
One cannot retire more capacity than existing capacity.
```math
\begin{aligned}
&\Delta_{y,z} \leq \overline{\Delta_{y,z}}
	\hspace{4 cm}  \forall y \in \mathcal{G}, z \in \mathcal{Z}
\end{aligned}
```
For resources where $\overline{\Omega_{y,z}}$ and $\underline{\Omega_{y,z}}$ is defined, then we impose constraints on minimum and maximum power capacity.
```math
\begin{aligned}
& \Delta^{total}_{y,z} \leq \overline{\Omega}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{G}, z \in \mathcal{Z} \\
& \Delta^{total}_{y,z}  \geq \underline{\Omega}_{y,z}
	\hspace{4 cm}  \forall y \in \mathcal{G}, z \in \mathcal{Z}
\end{aligned}
```
In addition, this function adds investment and fixed O\&M related costs related to discharge/generation capacity to the objective function:
```math
\begin{aligned}
& 	\sum_{y \in \mathcal{G} } \sum_{z \in \mathcal{Z}}
	\left( (\pi^{INVEST}_{y,z} \times \overline{\Omega}^{size}_{y,z} \times  \Omega_{y,z})
	+ (\pi^{FOM}_{y,z} \times \overline{\Omega}^{size}_{y,z} \times  \Delta^{total}_{y,z})\right)
\end{aligned}
```
"""
function investment_discharge(EP::Model, inputs::Dict, MinCapReq::Int)

	println("Investment Discharge Module")

	dfGen = inputs["dfGen"]

	G = inputs["G"] # Number of resources (generators, storage, DR, and DERs)

	NEW_CAP = inputs["NEW_CAP"] # Set of all resources eligible for new capacity
	RET_CAP = inputs["RET_CAP"] # Set of all resources eligible for capacity retirements
	COMMIT = inputs["COMMIT"] # Set of all resources eligible for unit commitment

	### Variables ###

	# Retired capacity of resource "y" from existing capacity
	@variable(EP, vRETCAP[y in RET_CAP] >= 0);

    # New installed capacity of resource "y"
	@variable(EP, vCAP[y in NEW_CAP] >= 0);

	### Expressions ###

	# Cap_Size is set to 1 for all variables when unit UCommit == 0
	# When UCommit > 0, Cap_Size is set to 1 for all variables except those where THERM == 1
	@expression(EP, eTotalCap[y in 1:G],
		if y in intersect(NEW_CAP, RET_CAP) # Resources eligible for new capacity and retirements
			if y in COMMIT
				dfGen[!,:Existing_Cap_MW][y] + dfGen[!,:Cap_Size][y]*(EP[:vCAP][y] - EP[:vRETCAP][y])
			else
				dfGen[!,:Existing_Cap_MW][y] + EP[:vCAP][y] - EP[:vRETCAP][y]
			end
		elseif y in setdiff(NEW_CAP, RET_CAP) # Resources eligible for only new capacity
			if y in COMMIT
				dfGen[!,:Existing_Cap_MW][y] + dfGen[!,:Cap_Size][y]*EP[:vCAP][y]
			else
				dfGen[!,:Existing_Cap_MW][y] + EP[:vCAP][y]
			end
		elseif y in setdiff(RET_CAP, NEW_CAP) # Resources eligible for only capacity retirements
			if y in COMMIT
				dfGen[!,:Existing_Cap_MW][y] - dfGen[!,:Cap_Size][y]*EP[:vRETCAP][y]
			else
				dfGen[!,:Existing_Cap_MW][y] - EP[:vRETCAP][y]
			end
		else # Resources not eligible for new capacity or retirements
			dfGen[!,:Existing_Cap_MW][y] + EP[:vZERO]
		end
	)

	## Objective Function Expressions ##

	# Fixed costs for resource "y" = annuitized investment cost plus fixed O&M costs
	# If resource is not eligible for new capacity, fixed costs are only O&M costs
	@expression(EP, eCFix[y in 1:G],
		if y in NEW_CAP # Resources eligible for new capacity
			if y in COMMIT
				dfGen[!,:Inv_Cost_per_MWyr][y]*dfGen[!,:Cap_Size][y]*vCAP[y] + dfGen[!,:Fixed_OM_Cost_per_MWyr][y]*eTotalCap[y]
			else
				dfGen[!,:Inv_Cost_per_MWyr][y]*vCAP[y] + dfGen[!,:Fixed_OM_Cost_per_MWyr][y]*eTotalCap[y]
			end
		else
			dfGen[!,:Fixed_OM_Cost_per_MWyr][y]*eTotalCap[y]
		end
	)

	# Sum individual resource contributions to fixed costs to get total fixed costs
	@expression(EP, eTotalCFix, sum(EP[:eCFix][y] for y in 1:G))

	# Add term to objective function expression
	EP[:eObj] += eTotalCFix

	### Constratints ###

	## Constraints on retirements and capacity additions
	# Cannot retire more capacity than existing capacity
	@constraint(EP, cMaxRetNoCommit[y in setdiff(RET_CAP,COMMIT)], vRETCAP[y] <= dfGen[!,:Existing_Cap_MW][y])
	@constraint(EP, cMaxRetCommit[y in intersect(RET_CAP,COMMIT)], dfGen[!,:Cap_Size][y]*vRETCAP[y] <= dfGen[!,:Existing_Cap_MW][y])

	## Constraints on new built capacity
	# Constraint on maximum capacity (if applicable) [set input to -1 if no constraint on maximum capacity]
	# DEV NOTE: This constraint may be violated in some cases where Existing_Cap_MW is >= Max_Cap_MW and lead to infeasabilty
	@constraint(EP, cMaxCap[y in intersect(dfGen[dfGen.Max_Cap_MW.>0,:R_ID], 1:G)], eTotalCap[y] <= dfGen[!,:Max_Cap_MW][y])

	# Constraint on minimum capacity (if applicable) [set input to -1 if no constraint on minimum capacity]
	# DEV NOTE: This constraint may be violated in some cases where Existing_Cap_MW is <= Min_Cap_MW and lead to infeasabilty
	@constraint(EP, cMinCap[y in intersect(dfGen[dfGen.Min_Cap_MW.>0,:R_ID], 1:G)], eTotalCap[y] >= dfGen[!,:Min_Cap_MW][y])

	if (MinCapReq == 1)
		@expression(EP, eMinCapResInvest[mincap = 1:inputs["NumberOfMinCapReqs"]], sum(EP[:eTotalCap][y] for y in dfGen[(dfGen[!,Symbol("MinCapTag_$mincap")].== 1) ,:][!,:R_ID]))
		EP[:eMinCapRes] += eMinCapResInvest
	end

	return EP
end

max_stars_repo_name: Betristor/LTESOM

max_issues_repo_name: Betristor/LTESOM

max_forks_repo_name: Betristor/LTESOM


text:
module ShaleDrillingLikelihood_OverallData_Test

# using Revise
# using Juno
# using Profile
# using ProfileView

using ShaleDrillingLikelihood
using Test
using StatsFuns
using Calculus
using Optim
using Random
using InteractiveUtils
using LinearAlgebra
using BenchmarkTools

using Base: OneTo

using ShaleDrillingLikelihood: AbstractDataSet,
    DataRoyalty, DataProduce, DataDrill,
    SimulationDraws,
    _nparm,
    TestDrillModel,
    theta_drill_ρ,
    theta_drill_ψ,
    theta_royalty_ρ,
    theta_produce_ψ,
    idx_drill, theta_drill,
    idx_royalty, theta_royalty,
    idx_produce, theta_produce,
    idx_drill_ρ, idx_royalty_ρ,
    idx_drill_ψ, idx_produce_ψ,
    split_thetas,
    merge_thetas,
    thetas,
    reward

println("testing overall royalty")

@testset "Access parameter vectors of joint model(s)" begin

    num_i = 500
    L = 3
    M = 500
    u = randn(num_i)
    v = randn(num_i)
    psi2 = u

    sim = SimulationDraws(M, num_i)

    # set up coefs
    θρ = 0.5
    αψ = 0.6

    # set up drilling coefs
    model_drill = TestDrillModel()
    θ_drill   = [αψ, 2.0, -2.0, -0.75, θρ]
    @test theta_drill_ρ(reward(model_drill), θ_drill) == θρ
    @test theta_drill_ψ(reward(model_drill), θ_drill) == αψ

    # set up royalty coefs
    θ_royalty = [θρ, 1.0,    -1.0, 1.0, 1.0,    -0.6, 0.6]  # dψdρ, ψ, β, κ
    @test theta_royalty_ρ(RoyaltyModel(), θ_royalty) == θρ

    # set up pdxn coefs
    θ_produce = vcat(αψ, rand(3), 0.3, 0.4)
    @test theta_produce_ψ(ProductionModel(), θ_produce) == αψ

    # big theta
    θ = vcat(θ_drill, θ_royalty[2:end], θ_produce[2:end])
    @test theta_drill(reward(model_drill), θ) == θ_drill

    # make data
    data_drill_opt = (num_zt=200, minmaxleases=1:2, nper_initial=10:20, nper_development=0:10, tstart=1:50,)
    data_drill   = DataDrill(u, v, TestDrillModel(), θ_drill; data_drill_opt...)
    data_produce = DataProduce(u, 10, 10:20, θ_produce)
    data_royalty = DataRoyalty(u,v,θ_royalty,L)

    @test idx_produce_ψ isa Function
    @test idx_drill_ψ isa Function
    coef_links = [(idx_produce_ψ, idx_drill_ψ,),]
    @test coef_links isa Vector{<:NTuple{2,Function}}

    data = DataSetofSets(data_drill, data_royalty, data_produce, coef_links)
    data_no_coeflinks = DataSetofSets(data_drill, data_royalty, data_produce)
    dataroypdxn = DataSetofSets(EmptyDataSet(), data_royalty, data_produce)
    datadrillonly = DataSetofSets(data_drill, EmptyDataSet(), EmptyDataSet() )

    @test length(data_drill) == num_i
    @test length(data_royalty) == num_i
    @test length(data_produce) == num_i
    @test all(length.(data) .== num_i)
    @test all(_nparm.(data) .== (_nparm(data_drill), _nparm(data_royalty), _nparm(data_produce)))

    @test _nparm(data) == length(vcat(θ_drill, θ_royalty[2:end], θ_produce[2:end]))
    @test _nparm(dataroypdxn) == length(vcat(θ_royalty, θ_produce))
    @test _nparm(datadrillonly) == length(θ_drill)

    @test (θ_drill, θ_royalty, θ_produce) == split_thetas(data, vcat(θ_drill, θ_royalty[2:end], θ_produce[2:end]))
    @test (θ_drill, θ_royalty, θ_produce) == split_thetas(data_no_coeflinks, vcat(θ_drill, θ_royalty[2:end], θ_produce))
    @test ([], θ_royalty, θ_produce,)     == split_thetas(dataroypdxn, vcat(θ_royalty, θ_produce))
    @test (θ_drill, [], []) == split_thetas(datadrillonly, θ_drill)

    @test theta_drill(  data, θ) == θ_drill
    @test theta_royalty(data, θ) == θ_royalty
    @test theta_produce(data, θ) == θ_produce
end

end # module


text:
# NOTE:
# This test file does end to end tests, while we don't want to check the details in the generated documents.
# Rather, we just assert we can `weave` all the supported formats without errors here.

# TODO:
# - more complex example
# - integration with other libraries, like Plots

@testset "end2end simple" begin
    include("test_simple.jl")
end  # @testset "end2end"


text:
###############################################################################
# TASK ########################################################################
###### The `task` field should not be accessed outside this section. ##########

function set_primal_start(task::Mosek.MSKtask, col::ColumnIndex, value::Float64)
    xx = Float64[value]
    for sol in [MSK_SOL_BAS, MSK_SOL_ITG]
        putxxslice(task, sol, col.value, col.value + Int32(1), xx)
    end
end
function set_primal_start(m::Optimizer, vi::MOI.VariableIndex, value::Float64)
    set_primal_start(m.task, mosek_index(m, vi), value)
end
function set_primal_start(task::Mosek.MSKtask, cols::ColumnIndices,
                          values::Vector{Float64})
    for sol in [MSK_SOL_BAS, MSK_SOL_ITG]
        if solutiondef(task, sol)
            xx = getxx(task, sol)
        else
            xx = zeros(Float64, getnumvar(task))
        end
        xx[cols.values] = values
        putxx(task, sol, xx)
    end
end
function set_primal_start(m::Optimizer, vis::Vector{MOI.VariableIndex},
                          values::Vector{Float64})
    if all(vi -> is_scalar(m, vi), vis)
        set_primal_start(m.task, columns(m, vis), values)
    else
        for (vi, value) in zip(vis, values)
            set_primal_start(m.task, mosek_index(m, vi), value)
        end
    end
end

###############################################################################
## INDEXING ###################################################################
###############################################################################

function variable_primal(m::Optimizer, N, col::ColumnIndex)
    return m.solutions[N].xx[col.value]
end
function variable_primal(m::Optimizer, N, mat::MatrixIndex)
    d = m.sd_dim[mat.matrix]
    r = d - mat.column + 1
    #   #entries full Δ       #entries right Δ      #entries above in lower Δ
    k = div((d + 1) * d, 2) - div((r + 1) * r, 2) + (mat.row - mat.column + 1)
    return m.solutions[N].barxj[mat.matrix][k]
end
function variable_primal(m::Optimizer, N, vi::MOI.VariableIndex)
    return variable_primal(m, N, mosek_index(m, vi))
end

###############################################################################
# MOI #########################################################################
###############################################################################

#### objective
function MOI.get(m::Optimizer, attr::MOI.ObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return getprimalobj(m.task, m.solutions[attr.result_index].whichsol)
end
function MOI.get(m::Optimizer, attr::MOI.DualObjectiveValue)
    MOI.check_result_index_bounds(m, attr)
    return getdualobj(m.task, m.solutions[attr.result_index].whichsol)
end

MOI.get(m::Optimizer,attr::MOI.ObjectiveBound) = getdouinf(m.task,MSK_DINF_MIO_OBJ_BOUND)

MOI.get(m::Optimizer,attr::MOI.RelativeGap) = getdouinf(m.task,MSK_DINF_MIO_OBJ_REL_GAP)

MOI.get(m::Optimizer,attr::MOI.SolveTime) = getdouinf(m.task,MSK_DINF_OPTIMIZER_TIME)

# NOTE: The MOSEK interface currently only supports Min and Max
function MOI.get(model::Optimizer, ::MOI.ObjectiveSense)
    if model.feasibility
        return MOI.FEASIBILITY_SENSE
    else
        sense = getobjsense(model.task)
        if sense == MSK_OBJECTIVE_SENSE_MINIMIZE
            MOI.MIN_SENSE
        else
            MOI.MAX_SENSE
        end
    end
end

function MOI.set(model::Optimizer,
                 attr::MOI.ObjectiveSense,
                 sense::MOI.OptimizationSense)
    if sense == MOI.MIN_SENSE
        model.feasibility = false
        putobjsense(model.task,MSK_OBJECTIVE_SENSE_MINIMIZE)
    elseif sense == MOI.MAX_SENSE
        model.feasibility = false
        putobjsense(model.task,MSK_OBJECTIVE_SENSE_MAXIMIZE)
    else
        @assert sense == MOI.FEASIBILITY_SENSE
        model.feasibility = true
        putobjsense(model.task,MSK_OBJECTIVE_SENSE_MINIMIZE)
        MOI.set(model,
                MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
                MOI.ScalarAffineFunction(MOI.ScalarAffineTerm{Float64}[], 0.0))
    end
end

#### Solver/Solution information

function MOI.get(m::Optimizer,attr::MOI.SimplexIterations)
    miosimiter = getlintinf(m.task,MSK_LIINF_MIO_SIMPLEX_ITER)
    if miosimiter > 0
        Int(miosimiter)
    else
        Int(getintinf(m.task,MSK_IINF_SIM_PRIMAL_ITER) + getintinf(m.task,MSK_IINF_SIM_DUAL_ITER))
    end
end

function MOI.get(m::Optimizer,attr::MOI.BarrierIterations)
    miosimiter = getlintinf(m.task,MSK_LIINF_MIO_INTPNT_ITER)
    if miosimiter > 0
        Int(miosimiter)
    else
        Int(getintinf(m.task,MSK_IINF_INTPNT_ITER))
    end
end

function MOI.get(m::Optimizer,attr::MOI.NodeCount)
        Int(getintinf(m.task,MSK_IINF_MIO_NUM_BRANCH))
end

MOI.get(m::Optimizer,attr::MOI.RawSolver) = m.task

MOI.get(m::Optimizer,attr::MOI.ResultCount) = length(m.solutions)

#### Problem information

function MOI.get(model::Optimizer,
                 ::MOI.NumberOfConstraints{MOI.ScalarAffineFunction{Float64},
                                           S}) where S <: ScalarLinearDomain
    F = MOI.ScalarAffineFunction{Float64}
    return count(id -> MOI.is_valid(model, MOI.ConstraintIndex{F, S}(id)),
                 allocatedlist(model.c_block))
end
function MOI.get(model::Optimizer,
                 ::MOI.ListOfConstraintIndices{MOI.ScalarAffineFunction{Float64},
                                               S}) where S <: ScalarLinearDomain
    F = MOI.ScalarAffineFunction{Float64}
    ids = filter(id -> MOI.is_valid(model, MOI.ConstraintIndex{F, S}(id)),
                 allocatedlist(model.c_block))
    return [MOI.ConstraintIndex{F, S}(id) for id in ids]
end
function MOI.get(model::Optimizer,
                 ::MOI.NumberOfConstraints{MOI.SingleVariable, S}) where S<:Union{ScalarLinearDomain,
                                                                                MOI.Integer}
    F = MOI.SingleVariable
    return count(id -> MOI.is_valid(model, MOI.ConstraintIndex{F, S}(id)),
                 allocatedlist(model.x_block))
end
function MOI.get(model::Optimizer,
                 ::MOI.ListOfConstraintIndices{MOI.SingleVariable, S}) where S<:Union{ScalarLinearDomain,
                                                                                MOI.Integer}
    F = MOI.SingleVariable
    ids = filter(id -> MOI.is_valid(model, MOI.ConstraintIndex{F, S}(id)),
                 allocatedlist(model.x_block))
    return [MOI.ConstraintIndex{F, S}(id) for id in ids]
end
function MOI.get(model::Optimizer,
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables, S}) where S<:VectorCone
    F = MOI.VectorOfVariables
    return count(i -> MOI.is_valid(model, MOI.ConstraintIndex{F, S}(i)),
                 eachindex(model.variable_to_vector_constraint_id))
end
function MOI.get(model::Optimizer,
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables, S}) where S<:VectorCone
    F = MOI.VectorOfVariables
    ids = filter(i -> MOI.is_valid(model, MOI.ConstraintIndex{F, S}(i)),
                 eachindex(model.variable_to_vector_constraint_id))
    return [MOI.ConstraintIndex{F, S}(id) for id in ids]
end
function MOI.get(model::Optimizer,
                 ::MOI.NumberOfConstraints{MOI.VectorOfVariables,
                                           MOI.PositiveSemidefiniteConeTriangle})
    # TODO this only works because deletion of PSD constraints is not supported yet
    return length(model.sd_dim)
end
function MOI.get(model::Optimizer,
                 ::MOI.ListOfConstraintIndices{MOI.VectorOfVariables,
                                               MOI.PositiveSemidefiniteConeTriangle})
    # TODO this only works because deletion of PSD constraints is not supported yet
    return [MOI.ConstraintIndex{MOI.VectorOfVariables,
                                MOI.PositiveSemidefiniteConeTriangle}(id) for id in 1:length(model.sd_dim)]
end

function MOI.get(model::Optimizer,
                 ::MOI.ListOfConstraints)
    list = Tuple{DataType, DataType}[]
    F = MOI.SingleVariable
    for D in [MOI.LessThan{Float64}, MOI.GreaterThan{Float64},
              MOI.EqualTo{Float64}, MOI.Interval{Float64},
              MOI.Integer]
        if !iszero(MOI.get(model, MOI.NumberOfConstraints{F, D}()))
            push!(list, (F, D))
        end
    end
    F = MOI.ScalarAffineFunction{Float64}
    for D in [MOI.LessThan{Float64}, MOI.GreaterThan{Float64},
              MOI.EqualTo{Float64}, MOI.Interval{Float64}]
        if !iszero(MOI.get(model, MOI.NumberOfConstraints{F, D}()))
            push!(list, (F, D))
        end
    end
    F = MOI.VectorOfVariables
    for D in [MOI.SecondOrderCone, MOI.RotatedSecondOrderCone,
              # TODO reenable for Mosek 9
              #MOI.PowerCone{Float64}, MOI.DualPowerCone{Float64},
              #MOI.ExponentialCone, MOI.DualExponentialCone,
              MOI.PositiveSemidefiniteConeTriangle]
        if !iszero(MOI.get(model, MOI.NumberOfConstraints{F, D}()))
            push!(list, (F, D))
        end
    end
    return list
end

#### Warm start values

MOI.supports(::Optimizer, ::MOI.VariablePrimalStart, ::Type{MOI.VariableIndex}) = true

function MOI.set(m::Optimizer, ::MOI.VariablePrimalStart,
                 v::MOI.VariableIndex, ::Nothing)
    set_primal_start(m, v, 0.0)
end
function MOI.set(m::Optimizer, ::MOI.VariablePrimalStart,
                 v::MOI.VariableIndex, val::Float64)
    set_primal_start(m, v, val)
end

function MOI.set(m::Optimizer, ::MOI.VariablePrimalStart,
                 vis::Vector{MOI.VariableIndex}, values::Vector{Float64})
    set_primal_start(m, vis, values)
end

# function MOI.set(m::Optimizer,attr::MOI.ConstraintDualStart, vs::Vector{MOI.ConstraintIndex}, vals::Vector{Float64})
#     subj = columns(vs)

#     for sol in [ MSK_SOL_BAS, MSK_SOL_ITG ]
#         if solutiondef(m.task,sol)
#             xx = getxx(m.task,sol)
#             xx[subj] = vals
#             putxx(m.task,sol,xx)
#         else
#             xx = zeros(Float64,getnumvar(m.task))
#             xx[subj] = vals
#             putxx(m.task,sol,xx)
#         end
#     end
# end

#### Variable solution values

function MOI.get(m::Optimizer, attr::MOI.VariablePrimal, vi::MOI.VariableIndex)
    MOI.check_result_index_bounds(m, attr)
    return variable_primal(m, attr.N, vi)
end
function MOI.get!(output::Vector{Float64}, m::Optimizer,
                  attr::MOI.VariablePrimal, vs::Vector{MOI.VariableIndex})
    MOI.check_result_index_bounds(m, attr)
    @assert eachindex(output) == eachindex(vs)
    for i in eachindex(output)
        output[i] = MOI.get(m, attr, vs[i])
    end
end
function MOI.get(m::Optimizer, attr::MOI.VariablePrimal,
                 vs::Vector{MOI.VariableIndex})
    MOI.check_result_index_bounds(m, attr)
    output = Vector{Float64}(undef, length(vs))
    MOI.get!(output, m, attr, vs)
    return output
end

#### Constraint solution values

function MOI.get(
    m     ::Optimizer,
    attr  ::MOI.ConstraintPrimal,
    ci  ::MOI.ConstraintIndex{MOI.SingleVariable,D}) where D
    MOI.check_result_index_bounds(m, attr)
    col = column(m, _variable(ci))
    return m.solutions[attr.N].xx[col.value]
end

# Semidefinite domain for a variable
function MOI.get!(
    output::Vector{Float64},
    m     ::Optimizer,
    attr  ::MOI.ConstraintPrimal,
    ci  ::MOI.ConstraintIndex{MOI.VectorOfVariables,
                              MOI.PositiveSemidefiniteConeTriangle})
    MOI.check_result_index_bounds(m, attr)
    whichsol = getsolcode(m,attr.N)
    output[1:length(output)] = reorder(getbarxj(m.task, whichsol, ci.value),
                                       MOI.PositiveSemidefiniteConeTriangle)
end

# Any other domain for variable vector
function MOI.get!(
    output::Vector{Float64},
    m     ::Optimizer,
    attr  ::MOI.ConstraintPrimal,
    ci  ::MOI.ConstraintIndex{MOI.VectorOfVariables,D}) where D
    MOI.check_result_index_bounds(m, attr)
    cols = columns(m, ci)
    output[1:length(output)] = reorder(m.solutions[attr.N].xx[cols.values], D)
end

function MOI.get(m     ::Optimizer,
                 attr  ::MOI.ConstraintPrimal,
                 ci  ::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D}) where D
    MOI.check_result_index_bounds(m, attr)
    cid = ref2id(ci)
    subi = getindex(m.c_block,cid)
    return m.solutions[attr.N].xc[subi]
end

function _variable_constraint_dual(sol::MosekSolution, col::ColumnIndex,
                                   ::Type{<:Union{MOI.Interval{Float64}, MOI.EqualTo{Float64}}})
    return sol.slx[col.value] - sol.sux[col.value]
end
function _variable_constraint_dual(sol::MosekSolution, col::ColumnIndex, ::Type{MOI.GreaterThan{Float64}})
    return sol.slx[col.value]
end
function _variable_constraint_dual(sol::MosekSolution, col::ColumnIndex, ::Type{MOI.LessThan{Float64}})
    return -sol.sux[col.value]
end
function MOI.get(m::Optimizer, attr::MOI.ConstraintDual,
                 ci::MOI.ConstraintIndex{MOI.SingleVariable, S}) where S <: ScalarLinearDomain
    MOI.check_result_index_bounds(m, attr)
    col = column(m, _variable(ci))
    dual = _variable_constraint_dual(m.solutions[attr.N], col, S)
    if getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE
        return dual
    else
        return -dual
    end
end

getsolcode(m::Optimizer, N) = m.solutions[N].whichsol

# The dual or primal of an SDP variable block is returned in lower triangular
# form but the constraint is in upper triangular form.
function reorder(x, ::Type{MOI.PositiveSemidefiniteConeTriangle})
    n = div(isqrt(1 + 8length(x)) - 1, 2)
    @assert length(x) == MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(n))
    y = similar(x)
    k = 0
    for j in 1:n, i in j:n
        k += 1
        y[div((i - 1) * i, 2) + j] = x[k]

    end
    @assert k == length(x)
    return y
end

function reorder(x, ::Type{<:Union{MOI.ExponentialCone,
                                   MOI.DualExponentialCone}})
    return [x[3], x[2], x[1]]
end
reorder(x, ::Type{<:VectorCone}) = x

# Semidefinite domain for a variable
function MOI.get!(
    output::Vector{Float64},
    m     ::Optimizer,
    attr  ::MOI.ConstraintDual,
    ci  ::MOI.ConstraintIndex{MOI.VectorOfVariables, MOI.PositiveSemidefiniteConeTriangle})
    MOI.check_result_index_bounds(m, attr)
    whichsol = getsolcode(m,attr.N)
    # It is in fact a real constraint and cid is the id of an ordinary constraint
    dual = reorder(getbarsj(m.task, whichsol, ci.value),
                   MOI.PositiveSemidefiniteConeTriangle)
    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        output[1:length(output)] .= dual
    else
        output[1:length(output)] .= -dual
    end
end

# Any other domain for variable vector
function MOI.get!(
    output::Vector{Float64},
    m     ::Optimizer,
    attr  ::MOI.ConstraintDual,
    ci  ::MOI.ConstraintIndex{MOI.VectorOfVariables,D}) where D
    MOI.check_result_index_bounds(m, attr)

    xcid = ref2id(ci)
    @assert(xcid > 0)

    cols = columns(m, ci)

    idx = reorder(1:length(output), D)

    if (getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE)
        output[idx] = m.solutions[attr.N].snx[cols.values]
    else
        output[idx] = -m.solutions[attr.N].snx[cols.values]
    end
end

function MOI.get(m     ::Optimizer,
                 attr  ::MOI.ConstraintDual,
                 ci  ::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64},D}) where D
    MOI.check_result_index_bounds(m, attr)

    cid = ref2id(ci)
    subi = getindex(m.c_block, cid)

    if getobjsense(m.task) == MSK_OBJECTIVE_SENSE_MINIMIZE
        m.solutions[attr.N].y[subi]
    else
        - m.solutions[attr.N].y[subi]
    end
end

solsize(m::Optimizer, ::MOI.ConstraintIndex{<:MOI.AbstractScalarFunction}) = 1
function solsize(m::Optimizer, ci::MOI.ConstraintIndex{MOI.VectorOfVariables})
    return getconeinfo(m.task, cone_id(m, ci))[3]
end
function solsize(m::Optimizer,
                 ci::MOI.ConstraintIndex{MOI.VectorOfVariables,
                                         MOI.PositiveSemidefiniteConeTriangle})
    d = m.sd_dim[ci.value]
    return MOI.dimension(MOI.PositiveSemidefiniteConeTriangle(d))
end

function MOI.get(m::Optimizer,
                 attr::Union{MOI.ConstraintPrimal, MOI.ConstraintDual},
                 ci::MOI.ConstraintIndex{<:MOI.AbstractVectorFunction})
    MOI.check_result_index_bounds(m, attr)
    cid = ref2id(ci)
    output = Vector{Float64}(undef, solsize(m, ci))
    MOI.get!(output, m, attr, ci)
    return output
end

#### Status codes
function MOI.get(m::Optimizer, attr::MOI.RawStatusString)
    if     m.trm === nothing
        return "MOI.OPTIMIZE_NOT_CALLED"
    elseif m.trm == MSK_RES_OK
        return join([string(sol.solsta) for sol in m.solutions], ", ")
    else
        return string(m.trm)
    end

end
function MOI.get(m::Optimizer, attr::MOI.TerminationStatus)
    if     m.trm === nothing
        MOI.OPTIMIZE_NOT_CALLED
    elseif m.trm == MSK_RES_OK
        # checking `any(sol -> sol.solsta == MSK_SOL_STA_PRIM_INFEAS_CER, m.solutions)`
        # doesn't work for MIP as there is not certificate, i.e. the solutions status is
        # `UNKNOWN`, only the problem status is `INFEAS`.
        if any(sol -> sol.prosta == MSK_PRO_STA_PRIM_INFEAS, m.solutions)
            MOI.INFEASIBLE
        elseif any(sol -> sol.prosta == MSK_PRO_STA_DUAL_INFEAS, m.solutions)
            MOI.DUAL_INFEASIBLE
        else
            MOI.OPTIMAL
        end
    elseif m.trm == MSK_RES_TRM_MAX_ITERATIONS
        MOI.ITERATION_LIMIT
    elseif m.trm == MSK_RES_TRM_MAX_TIME
        MOI.TIME_LIMIT
    elseif m.trm == MSK_RES_TRM_OBJECTIVE_RANGE
        MOI.OBJECTIVE_LIMIT
    elseif m.trm == MSK_RES_TRM_MIO_NUM_RELAXS
        MOI.OTHER_LIMIT
    elseif m.trm == MSK_RES_TRM_MIO_NUM_BRANCHES
        MOI.NODE_LIMIT
    elseif m.trm == MSK_RES_TRM_NUM_MAX_NUM_INT_SOLUTIONS
        MOI.SOLUTION_LIMIT
    elseif m.trm == MSK_RES_TRM_STALL
        MOI.SLOW_PROGRESS
    elseif m.trm == MSK_RES_TRM_USER_CALLBACK
        MOI.INTERRUPTED
    elseif m.trm == MSK_RES_TRM_MAX_NUM_SETBACKS
        MOI.OTHER_LIMIT
    elseif m.trm == MSK_RES_TRM_NUMERICAL_PROBLEM
        MOI.SLOW_PROGRESS
    elseif m.trm == MSK_RES_TRM_INTERNAL
        MOI.OTHER_ERROR
    elseif m.trm == MSK_RES_TRM_INTERNAL_STOP
        MOI.OTHER_ERROR
    else
        MOI.OTHER_ERROR
    end
end

function MOI.get(m::Optimizer, attr::MOI.PrimalStatus)
    if attr.N > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    solsta = m.solutions[attr.N].solsta
    if     solsta == MSK_SOL_STA_UNKNOWN
        MOI.UNKNOWN_RESULT_STATUS
    elseif solsta == MSK_SOL_STA_OPTIMAL
        MOI.FEASIBLE_POINT
    elseif solsta == MSK_SOL_STA_PRIM_FEAS
        MOI.FEASIBLE_POINT
    elseif solsta == MSK_SOL_STA_DUAL_FEAS
        MOI.UNKNOWN_RESULT_STATUS
    elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS
        MOI.FEASIBLE_POINT
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        MOI.NO_SOLUTION
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        MOI.INFEASIBILITY_CERTIFICATE
    elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER
        MOI.NO_SOLUTION
    elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER
        MOI.REDUCTION_CERTIFICATE
    elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL
        MOI.FEASIBLE_POINT
    else
        MOI.UNKNOWN_RESULT_STATUS
    end
end

function MOI.get(m::Optimizer,attr::MOI.DualStatus)
    if attr.N > MOI.get(m, MOI.ResultCount())
        return MOI.NO_SOLUTION
    end
    solsta = m.solutions[attr.N].solsta
    if     solsta == MSK_SOL_STA_UNKNOWN
        MOI.UNKNOWN_RESULT_STATUS
    elseif solsta == MSK_SOL_STA_OPTIMAL
        MOI.FEASIBLE_POINT
    elseif solsta == MSK_SOL_STA_PRIM_FEAS
        MOI.UNKNOWN_RESULT_STATUS
    elseif solsta == MSK_SOL_STA_DUAL_FEAS
        MOI.FEASIBLE_POINT
    elseif solsta == MSK_SOL_STA_PRIM_AND_DUAL_FEAS
        MOI.FEASIBLE_POINT
    elseif solsta == MSK_SOL_STA_PRIM_INFEAS_CER
        MOI.INFEASIBILITY_CERTIFICATE
    elseif solsta == MSK_SOL_STA_DUAL_INFEAS_CER
        MOI.NO_SOLUTION
    elseif solsta == MSK_SOL_STA_PRIM_ILLPOSED_CER
        MOI.REDUCTION_CERTIFICATE
    elseif solsta == MSK_SOL_STA_DUAL_ILLPOSED_CER
        MOI.NO_SOLUTION
    elseif solsta == MSK_SOL_STA_INTEGER_OPTIMAL
        MOI.NO_SOLUTION
    else
        MOI.UNKNOWN_RESULT_STATUS
    end
end


text:
export CircuitComponent, SeriesComponent, ParallelComponent, TransmissionLine

"""
    abstract type CircuitComponent

An abstract representation of a circuit component, including lumped circuit
elements like a capacitor, inductor, or resistor, but also distributed structures
such as a transmission line.
"""
abstract type CircuitComponent end

# the empty string is used for the ground net
ground = ""

"""
    Blackbox(ω::AbstractVector{<:Number}, comp::CircuitComponent)

Create a `Blackbox` model for a `CircuitComponent`.
"""
function Blackbox(ω::AbstractVector{<:Number}, comp::CircuitComponent)
    @assert isreal(ω)
    return canonical_gauge(Blackbox(ω, PSOModel(comp)))
end

"""
    SeriesComponent(p1::AbstractString, p2::AbstractString,
        k::Number, g::Number, c::Number)

Constructs an inductance (`1/k`), resistance (`1/g`), and capacitance (`c`)
in parallel connecting two ports `p1` and `p2`:

             ┌─ 1/k ─┐
    ─────┬───┼─ 1/g ─┼───┬─────
      p1 ╎   └─  c  ─┘   ╎ p2
         ╎               ╎
    ─────┴───────────────┴─────
"""
struct SeriesComponent <: CircuitComponent
    p1::String
    p2::String
    k::Float64
    g::Float64
    c::Float64
    function SeriesComponent(p1::AbstractString, p2::AbstractString, k::Number,
            g::Number, c::Number)
        @assert p1 != ground
        @assert p2 != ground
        @assert isreal(k) && isreal(g) && isreal(c)
        return new(p1, p2, k, g, c)
    end
end

"""
    Circuit(comp::SeriesComponent)

Create a `Circuit` model for a `SeriesComponent`.
"""
function Circuit(comp::SeriesComponent)
    c = Circuit([ground, comp.p1, comp.p2])
    set_inv_inductance!(c, comp.p1, comp.p2, comp.k)
    set_conductance!(c, comp.p1, comp.p2, comp.g)
    set_capacitance!(c, comp.p1, comp.p2, comp.c)
    return c
end

"""
    PSOModel(comp::SeriesComponent)

Create a `PSOModel` for a `SeriesComponent`.
"""
function PSOModel(comp::SeriesComponent)
    return PSOModel(Circuit(comp), [(comp.p1, ground), (comp.p2, ground)], [comp.p1, comp.p2])
end

"""
    ParallelComponent(p::AbstractString, k::Number, g::Number, c::Number)

Constructs an inductance (`1/k`), resistance (`1/g`), and capacitance (`c`) in
parallel across a port `p`:

    ──────────┬──────────
              │ port p
          ┌───┼───┐
          │   │   │
         1/k  c  1/g
          │   │   │
          └───┼───┘
              │
    ──────────┴──────────
"""
struct ParallelComponent <: CircuitComponent
    p::String
    k::Float64
    g::Float64
    c::Float64
    function ParallelComponent(p::AbstractString, k::Number, g::Number, c::Number)
        @assert p != ground
        @assert isreal(k) && isreal(g) && isreal(c)
        return new(p, k, g, c)
    end
end

"""
    Circuit(comp::ParallelComponent)

Create a `Circuit` model for a `ParallelComponent`.
"""
function Circuit(comp::ParallelComponent)
    c = Circuit([ground, comp.p])
    set_inv_inductance!(c, comp.p, ground, comp.k)
    set_conductance!(c, comp.p, ground, comp.g)
    set_capacitance!(c, comp.p, ground, comp.c)
    return c
end

"""
    PSOModel(comp::ParallelComponent)

Create a `PSOModel` for a `ParallelComponent`.
"""
PSOModel(comp::ParallelComponent) = PSOModel(Circuit(comp), [(comp.p, ground)], [comp.p])

struct TransmissionLine <: CircuitComponent
    ports::Vector{String}
    propagation_speed::Float64
    characteristic_impedance::Float64
    len::Float64 # length, but that's the name of a built in function
    port_locations::UniqueVector{Float64}
    δ::Float64 # ignored when converted directly to Blackbox
    function TransmissionLine(ports::AbstractVector{<:AbstractString},
            propagation_speed::Number, characteristic_impedance::Number,
            len::Number, locations::AbstractVector{<:Number}, δ::Number)
        @assert isreal(propagation_speed) && isreal(characteristic_impedance)
        @assert isreal(len) && isreal(δ)
        @assert !(ground in ports)
        @assert propagation_speed > zero(propagation_speed)
        @assert characteristic_impedance > zero(characteristic_impedance)
        @assert len > zero(len)
        @assert all(locations .> zero(eltype(locations)))
        @assert all(locations .< len)
        port_locations = UniqueVector([zero(len); locations; len])
        @assert length(ports) == length(port_locations)
        @assert δ >= zero(δ)
        return new(ports, propagation_speed, characteristic_impedance,
            len, port_locations, δ)
    end
end

"""
    TransmissionLine(ports::AbstractVector{<:AbstractString},
        propagation_speed::Number, characteristic_impedance::Number, len::Number;
        locations::AbstractVector{<:Number} = Float64[], δ::Number = len/100)

Create an abstract representation of a transmission line with given propagation
speed, characteristic impedance, and length. Ports will be present at the two
ends of the transmission line in addition to the `locations`. `δ` is the maximum
length of an LC stage in the circuit approximation to a transmission line. It is
ignored when constructing a Blackbox from a TransmissionLine.
"""
TransmissionLine(ports::AbstractVector{<:AbstractString}, propagation_speed::Number,
        characteristic_impedance::Number, len::Number;
        locations::AbstractVector{<:Number} = Float64[],
        δ::Number = len/100) =
    TransmissionLine(ports, propagation_speed, characteristic_impedance,
        len, locations, δ)

function Circuit_and_ports(comp::TransmissionLine,
        stages_per_segment::Union{Nothing, AbstractVector{Int}} = nothing)
    sort_inds = sortperm(comp.port_locations)
    port_locations, ports = comp.port_locations[sort_inds], comp.ports[sort_inds]
    ν, Z0, δ = comp.propagation_speed, comp.characteristic_impedance, comp.δ
    segment_lengths = port_locations[2:end] - port_locations[1:end-1]
    if isnothing(stages_per_segment)
        stages_per_segment = [ceil(Int, len / δ) for len in segment_lengths]
    end
    num_vertices = sum(stages_per_segment) + 2 # ground is 0
    c = Circuit(0:(num_vertices-1))
    v0 = 0
    port_edges = Tuple{Int, Int}[]
    for (segment_length, num_stages) in zip(segment_lengths, stages_per_segment)
        stage_length = float(segment_length / num_stages)
        inductance = stage_length * Z0 / ν
        capacitance = stage_length / (Z0 * ν)
        for v in (v0 + 1):(v0 + num_stages)
            set_inductance!(c, v, v + 1, inductance)
        end
        for v in [v0 + 1, v0 + num_stages + 1]
            set_capacitance!(c, v, 0, 0.5 * capacitance + get_capacitance(c, v, 0))
        end
        for v in (v0 + 2):(v0 + num_stages)
            set_capacitance!(c, v, 0, capacitance)
        end
        push!(port_edges, (v0 + 1, 0))
        v0 += num_stages
    end
    @assert v0 + 1 == num_vertices-1
    push!(port_edges, (num_vertices-1, 0))
    return c, port_edges, ports
end

"""
    Circuit(comp::TransmissionLine, vertex_prefix::AbstractString,
        stages_per_segment::Union{Nothing, AbstractVector{Int}} = nothing)

Create a `Circuit` model for a `TransmissonLine`.
"""
function Circuit(comp::TransmissionLine, vertex_prefix::AbstractString,
        stages_per_segment::Union{Nothing, AbstractVector{Int}} = nothing)
    c, port_edges, ports = Circuit_and_ports(comp, stages_per_segment)
    vertices = [ground; ["$(vertex_prefix)$v" for v in c.vertices[2:end]]]
    for (e, p) in zip(port_edges, ports)
        vertices[e[1]+1] = p
    end
    return Circuit(matrices(c)..., vertices)
end

"""
    PSOModel(comp::TransmissionLine)

Create a `PSOModel` for a `TransmissionLine`.
"""
PSOModel(comp::TransmissionLine) = PSOModel(Circuit_and_ports(comp)...)

"""
    Blackbox(ω::Vector{<:Number}, comp::TransmissionLine)

Create a `Blackbox` model from a `TransmissionLine`.
"""
function Blackbox(ω::Vector{<:Number}, comp::TransmissionLine)
    @assert isreal(ω)
    ν = comp.propagation_speed
    Z0 = comp.characteristic_impedance
    function Blackbox_from_params(ports::Vector, len::Real)
        u = exp.(-im * ω * len/ν)
        # S matrix is [[0, u] [u, 0]]
        # turn this into Y matrix using closed form expression
        denom = (1 .- u.^2) * Z0
        Y11_Y22 = (1 .+ u.^2) ./ denom
        Y12_Y21 = -2 * u ./ denom
        Y = [[[a, b] [b, a]] for (a, b) in zip(Y11_Y22, Y12_Y21)]
        U = eltype(eltype(Y))
        P = Matrix{U}(I, 2, 2)
        return Blackbox(ω, Y, P, ports)
    end
    sort_inds = sortperm(comp.port_locations)
    sort_port_locations = comp.port_locations[sort_inds]
    sort_ports = comp.ports[sort_inds]
    params = [(ports = [sort_ports[i], sort_ports[i + 1]],
              len = sort_port_locations[i+1] - sort_port_locations[i])
              for i in 1:(length(comp.ports)-1)]
    return connect([Blackbox_from_params(p...) for p in params])
end


text:
"""
    semidiscretise(u0func, semidisc::AbstractSemidiscretisation, tspan)

Apply the semidiscretisation `semidisc` to the initial data given by `u0func`
and return an `ODEProblem` with time span `tspan`.
"""
function semidiscretise(u0func, semidisc::AbstractSemidiscretisation, tspan)
    u0 = compute_coefficients(u0func, semidisc.derivative)
    ode = ODEProblem(semidisc, u0, tspan)
end

function evaluate_coefficients(u, semidisc::AbstractSemidiscretisation)
    evaluate_coefficients(u, semidisc.derivative)
end

function integrate(func, u, semidisc::AbstractSemidiscretisation)
    integrate(func, u, semidisc.derivative)
end

struct ScalarIntegralQuantities{T} <: FieldVector{2,T}
    mass::T
    energy::T
end

function DiffEqCallbacks.SavingCallback(semidisc::AbstractSemidiscretisation; kwargs...)
    T = eltype(semidisc.derivative)

    save_func = (u,t,integrator) -> integrate(u->ScalarIntegralQuantities(u,u^2),
                                                u, integrator.f)
    saved_values = SavedValues(T, ScalarIntegralQuantities{T})
    SavingCallback(save_func, saved_values; kwargs...)
end

function FilterCallback(filter::AbstractFilter)
    condition = (t, u, integrator) -> true
    affect!(integrator) = filter(integrator.u)

    DiscreteCallback(condition, affect!, save_positions=(false,false))
end


text:
"""
    function LPHGLET_Synthesis(dvec::Vector{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig; gltype::Symbol = :L, ϵ::Float64 = 0.3)

Perform Lapped-HGLET Synthesis transform

### Input Arguments
* `dvec`: the expansion coefficients corresponding to the chosen basis
* `GP`: a GraphPart object
* `BS`: a BasisSpec object
* `G`: a GraphSig object
* `gltype`: :L or :Lsym, indicating which eigenvectors are used
* `ϵ`: relative action bandwidth (default: 0.3)

### Output Argument
* `f`: the reconstructed signal
* `GS`: the reconstructed GraphSig object
"""
function LPHGLET_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig; gltype::Symbol = :L, ϵ::Float64 = 0.3)
    # Preliminaries
    W = G.W
    inds = GP.inds
    rs = GP.rs
    N = size(W, 1)
    jmax = size(rs, 2)
    Uf = Matrix{Float64}(I, N, N)
    used_node = Set()

    # fill in the appropriate entries of dmatrix
    dmatrix = dvec2dmatrix(dvec, GP, BS)

    f = zeros(size(dmatrix[:, jmax, :]))

    # Perform the synthesis transform
    for j = 1:jmax
        regioncount = count(!iszero, rs[:,j]) - 1
        # assemble orthogonal folding operator at level j - 1
        keep_folding!(Uf, used_node, W, GP; ϵ = ϵ, j = j - 1)
        for r = 1:regioncount
            # indices of current region
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            # indices of current region's nodes
            indrs = inds[indr, j]
            # number of nodes in current region
            n = length(indrs)

            # only proceed forward if coefficients do not exist
            if (j == jmax || count(!iszero, dmatrix[indr, j + 1, :]) == 0) && count(!iszero, dmatrix[indr, j, :]) > 0
                # compute the eigenvectors
                W_temp = W[indrs,indrs]
                D_temp = Diagonal(vec(sum(W_temp, dims = 1)))
                if gltype == :L
                    # compute the eigenvectors of L ==> svd(L)
                    v = svd(Matrix(D_temp - W_temp)).U
                elseif gltype == :Lsym
                    # check if one can assemble the Lsym
                    if minimum(sum(W[indrs, indrs], dims = 1)) > 10^3 * eps()
                        ### eigenvectors of L_sym ==> svd(L_sym)
                        D_temp_p = Diagonal(vec(sum(W_temp, dims = 1)).^(-0.5))
                        v = svd(Matrix(D_temp_p * (D_temp - W_temp) * D_temp_p)).U
                    else
                        ### eigenvectors of L ==> svd(L)
                        v = svd(Matrix(D_temp - W_temp)).U
                    end
                end
                v = v[:,end:-1:1]

                # standardize the eigenvector signs
                standardize_eigenvector_signs!(v)

                # construct unfolder operator custom to current region
                P = Uf[indrs, :]'

                # reconstruct the signal
                f += (P * v) * dmatrix[indr, j, :]

            end
        end
    end

    # creat a GraphSig object with the reconstructed data
    GS = deepcopy(G)
    replace_data!(GS, f)

    return f, GS
end

"""
    function LPHGLET_Analysis_All(G::GraphSig, GP::GraphPart; ϵ::Float64 = 0.3)

For a GraphSig object 'G', generate the 2 matrices of Lapped-HGLET expansion coefficients
corresponding to the eigenvectors of L and Lsym

### Input Arguments
* `G`:  a GraphSig object
* `GP`: a GraphPart object
* `ϵ`: relative action bandwidth (default: 0.3)

### Output Argument
* `dmatrixlH`:        the matrix of expansion coefficients for L
* `dmatrixlHsym`:     the matrix of expansion coefficients for Lsym
* `GP`:              a GraphPart object
"""
function LPHGLET_Analysis_All(G::GraphSig, GP::GraphPart; ϵ::Float64 = 0.3)
    # Preliminaries
    W = G.W
    inds = GP.inds
    rs = GP.rs
    N = size(W, 1)
    jmax = size(rs, 2)
    fcols = size(G.f, 2)
    Uf = Matrix{Float64}(I, N, N)
    used_node = Set()
    dmatrixlH = zeros(N, jmax, fcols)
    dmatrixlHsym = deepcopy(dmatrixlH)

    for j = 1:jmax
        regioncount = count(!iszero, rs[:,j]) - 1
        # assemble orthogonal folding operator at level j - 1
        keep_folding!(Uf, used_node, W, GP; ϵ = ϵ, j = j - 1)
        for r = 1:regioncount
            # indices of current region
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            # indices of current region's nodes
            indrs = inds[indr, j]
            # number of nodes in current region
            n = length(indrs)

            # compute the eigenvectors
            W_temp = W[indrs,indrs]
            D_temp = Diagonal(vec(sum(W_temp, dims = 1)))
            ## eigenvectors of L ==> svd(L)
            v = svd(Matrix(D_temp - W_temp)).U
            ## eigenvectors of L_sym ==> svd(L_sym)
            if minimum(sum(W[indrs, indrs], dims = 1)) > 10^3 * eps()
                ### eigenvectors of L_sym ==> svd(L_sym)
                D_temp_p = Diagonal(vec(sum(W_temp, dims = 1)).^(-0.5))
                v_sym = svd(Matrix(D_temp_p * (D_temp - W_temp) * D_temp_p)).U
            else
                ### eigenvectors of L ==> svd(L)
                v_sym = deepcopy(v)
            end

            # standardize the eigenvector signs
            v = v[:,end:-1:1]
            standardize_eigenvector_signs!(v)
            v_sym = v_sym[:,end:-1:1]
            standardize_eigenvector_signs!(v_sym)

            # construct unfolding operator custom to current region
            P = Uf[indrs, :]'
            # obtain the expansion coefficients
            dmatrixlH[indr, j, :] = (P * v)' * G.f
            dmatrixlHsym[indr, j, :] = (P * v_sym)' * G.f
        end
    end

    return dmatrixlH, dmatrixlHsym

end

function standardize_eigenvector_signs!(v)
    # standardize the eigenvector signs for HGLET (different with NGWPs)
    for col = 1:size(v, 2)
        row = 1
        standardized = false
        while !standardized
            if v[row, col] > 10^3 * eps()
                standardized = true
            elseif v[row,col] < -10^3 * eps()
                v[:, col] = -v[:, col]
            else
                row += 1
            end
        end
    end
end

"""
    HGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L)

assemble the whole HGLET dictionary

### Input Arguments
* `GP`: a GraphPart object
* `G`:  a GraphSig object
* `gltype`: `:L` or `:Lsym`

### Output Argument
* `dictionary`: the HGLET dictionary

"""
function HGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L)
    N = size(G.W, 1)
    jmax = size(GP.rs, 2)
    dictionary = zeros(N, jmax, N)
    for j = 1:jmax
        BS = BasisSpec(collect(enumerate(j * ones(Int, N))))
        dictionary[:, j, :] = HGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; gltype = gltype)[1]'
    end
    return dictionary
end

"""
    LPHGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L, ϵ::Float64 = 0.3)

assemble the whole LP-HGLET dictionary

### Input Arguments
* `GP`: a GraphPart object
* `G`:  a GraphSig object
* `gltype`: `:L` or `:Lsym`
* `ϵ`: relative action bandwidth (default: 0.3)

### Output Argument
* `dictionary`: the LP-HGLET dictionary

"""
function LPHGLET_dictionary(GP::GraphPart, G::GraphSig; gltype::Symbol = :L, ϵ::Float64 = 0.3)
    N = size(G.W, 1)
    jmax = size(GP.rs, 2)
    dictionary = zeros(N, jmax, N)
    for j = 1:jmax
        BS = BasisSpec(collect(enumerate(j * ones(Int, N))))
        dictionary[:, j, :] = LPHGLET_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; gltype = gltype, ϵ = ϵ)[1]'
    end
    return dictionary
end

function HGLET_DST4_Analysis(G::GraphSig, GP::GraphPart)
    # Preliminaries
    W = G.W
    ind = GP.ind
    rs = GP.rs
    N = size(G.W,1)
    jmax = size(rs,2)
    fcols = size(G.f,2)
    dmatrix = zeros(N,jmax,fcols)
    dmatrix[:,jmax,:] = G.f[ind,:]

    # Perform the HGLET analysis, i.e., generating the HGLET coefficients
    for j = jmax-1:-1:1
        regioncount = count(!iszero, rs[:,j]) - 1
        for r = 1:regioncount
            # the index that marks the start of the region
            rs1 = rs[r,j]

            # the index that is one after the end of the region
            rs3 = rs[r+1,j]

            # the number of points in the current region
            n = rs3 - rs1

            if n == 1
                dmatrix[rs1,j,:] = G.f[ind[rs1],:]
            elseif n > 1
                indrs = ind[rs1:rs3-1]
                W_temp = W[indrs,indrs]
                D_temp = Diagonal(vec(sum(W_temp, dims = 1)))
                L_temp = D_temp - W_temp
                L_temp[1, 1] = 3  # DST type-IV (left: Dirichlet, right: Neumann)
                v = svd(Matrix(L_temp)).U
                v = v[:,end:-1:1] # reorder the ev's in the decreasing ew's

                # standardize the eigenvector signs
                for col = 1:n
                    row = 1
                    standardized = false
                    while !standardized
                        if v[row,col] > 10^3*eps()
                            standardized = true
                        elseif v[row,col] < -10^3*eps()
                            v[:,col] = - v[:,col]
                            standardized = true
                        else
                            row += 1
                        end
                    end
                end

                # obtain the expansion coeffcients
                if gltype == :Lrw && normalizep
                    dmatrix[rs1:rs3-1,j,:] = v'*(D.^0.5)*G.f[indrs,:]
                else
                    dmatrix[rs1:rs3-1,j,:] = v'*G.f[indrs,:]
                end
            end
        end
    end
    return dmatrix
end

function HGLET_DST4_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec,
                         G::GraphSig)
    # Preliminaries
    W = G.W
    jmax = size(GP.rs,2)

    # Fill in the appropriate entries of dmatrix
    dmatrix = dvec2dmatrix(dvec,GP,BS)
    f = dmatrix[:,jmax,:]

    # Perform the signal synthesis from the given coefficients
    for j = jmax:-1:1
        regioncount = count(!iszero, GP.rs[:,j]) - 1
        for r = 1:regioncount
            # the index that marks the start of the region
            rs1 = GP.rs[r,j]

            # the index that is one after the end of the region
            rs3 = GP.rs[r+1,j]

            # the number of points in the current region
            n = rs3 - rs1

            # only proceed forward if coefficients do not exist
            if (j == jmax || count(!iszero, dmatrix[rs1:rs3-1,j+1,:]) == 0) && count(!iszero, dmatrix[rs1:rs3-1,j,:]) > 0

                if n == 1
                    f[rs1,:] = dmatrix[rs1,j,:]
                elseif n > 1
                    indrs = GP.ind[rs1:rs3-1]
                    W_temp = W[indrs,indrs]
                    D_temp = Diagonal(vec(sum(W_temp, dims = 1)))
                    L_temp = D_temp - W_temp
                    L_temp[1, 1] = 3  # DST type-IV (left: Dirichlet, right: Neumann)
                    v = svd(Matrix(L_temp)).U
                    v = v[:,end:-1:1] # reorder the ev's in the decreasing ew's

                    # standardize the eigenvector signs
                    standardize_eigenvector_signs!(v)

                    # reconstruct the signal
                    f[rs1:rs3-1,:] = v*dmatrix[rs1:rs3-1,j,:]
                end
            end
        end
    end

    # put the reconstructed values in the correct order
    f[GP.ind,:] = f

    # creat a GraphSig object with the reconstructed data
    GS = deepcopy(G)
    replace_data!(GS,f)

    return f, GS
end

function HGLET_DST4_dictionary(GP::GraphPart, G::GraphSig)
    N = size(G.W, 1)
    jmax = size(GP.rs, 2)
    dictionary = zeros(N, jmax, N)
    for j = 1:jmax
        BS = BasisSpec(collect(enumerate(j * ones(Int, N))))
        dictionary[:, j, :] = HGLET_DST4_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G)[1]'
    end
    return dictionary
end

function LPHGLET_DST4_Analysis(G::GraphSig, GP::GraphPart; ϵ::Float64 = 0.3)
    # Preliminaries
    W = G.W
    inds = GP.inds
    rs = GP.rs
    N = size(W, 1)
    jmax = size(rs, 2)
    fcols = size(G.f, 2)
    Uf = Matrix{Float64}(I, N, N)
    used_node = Set()
    dmatrixlH = zeros(N, jmax, fcols)

    for j = 1:jmax
        regioncount = count(!iszero, rs[:,j]) - 1
        # assemble orthogonal folding operator at level j - 1
        keep_folding!(Uf, used_node, W, GP; ϵ = ϵ, j = j - 1)
        for r = 1:regioncount
            # indices of current region
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            # indices of current region's nodes
            indrs = inds[indr, j]
            # number of nodes in current region
            n = length(indrs)

            # compute the eigenvectors
            W_temp = W[indrs,indrs]
            D_temp = Diagonal(vec(sum(W_temp, dims = 1)))
            L_temp = D_temp - W_temp
            L_temp[1, 1] = 3  # DST type-IV (left: Dirichlet, right: Neumann)
            v = svd(Matrix(L_temp)).U

            # standardize the eigenvector signs
            v = v[:,end:-1:1]
            standardize_eigenvector_signs!(v)

            # construct unfolding operator custom to current region
            P = Uf[indrs, :]'
            # obtain the expansion coefficients
            dmatrixlH[indr, j, :] = (P * v)' * G.f
        end
    end

    return dmatrixlH

end

function LPHGLET_DST4_Synthesis(dvec::Matrix{Float64}, GP::GraphPart, BS::BasisSpec, G::GraphSig; ϵ::Float64 = 0.3)
    # Preliminaries
    W = G.W
    inds = GP.inds
    rs = GP.rs
    N = size(W, 1)
    jmax = size(rs, 2)
    Uf = Matrix{Float64}(I, N, N)
    used_node = Set()

    # fill in the appropriate entries of dmatrix
    dmatrix = dvec2dmatrix(dvec, GP, BS)
    f = zeros(size(dmatrix[:, jmax, :]))

    # Perform the synthesis transform
    for j = 1:jmax
        regioncount = count(!iszero, rs[:,j]) - 1
        # assemble orthogonal folding operator at level j - 1
        MultiscaleGraphSignalTransforms.keep_folding!(Uf, used_node, W, GP; ϵ = ϵ, j = j - 1)
        for r = 1:regioncount
            # indices of current region
            indr = rs[r, j]:(rs[r + 1, j] - 1)
            # indices of current region's nodes
            indrs = inds[indr, j]
            # number of nodes in current region
            n = length(indrs)

            # only proceed forward if coefficients do not exist
            if (j == jmax || count(!iszero, dmatrix[indr, j + 1, :]) == 0) && count(!iszero, dmatrix[indr, j, :]) > 0
                # compute the eigenvectors
                W_temp = W[indrs,indrs]
                D_temp = Diagonal(vec(sum(W_temp, dims = 1)))
                L_temp = D_temp - W_temp
                L_temp[1, 1] = 3  # DST type-IV (left: Dirichlet, right: Neumann)
                v = svd(Matrix(L_temp)).U
                v = v[:,end:-1:1]

                # standardize the eigenvector signs
                standardize_eigenvector_signs!(v)

                # construct unfolder operator custom to current region
                P = Uf[indrs, :]'

                # reconstruct the signal
                f += (P * v) * dmatrix[indr, j, :]

            end
        end
    end

    # creat a GraphSig object with the reconstructed data
    GS = deepcopy(G)
    replace_data!(GS, f)

    return f, GS
end

function LPHGLET_DST4_dictionary(GP::GraphPart, G::GraphSig; ϵ::Float64 = 0.3)
    N = size(G.W, 1)
    jmax = size(GP.rs, 2)
    dictionary = zeros(N, jmax, N)
    for j = 1:jmax
        BS = BasisSpec(collect(enumerate(j * ones(Int, N))))
        dictionary[:, j, :] = LPHGLET_DST4_Synthesis(Matrix{Float64}(I, N, N), GP, BS, G; ϵ = ϵ)[1]'
    end
    return dictionary
end


text:
"""
    getcolourrange(N)

Return N colours for plotting.

  # Example
  ```julia-repl
julia> using PyPlot
julia> clrs = getcolourrange(4)
julia> x = LinRange(0, 2π, 100);
julia> for i in 1:4
          plot(x, sin.(x).+i, color=clrs[i])
       end
```
"""
function getcolourrange(N::Int)

  c = get_cmap("gist_rainbow", N)

  aux = c(LinRange(0,1,N))[:,1:3]

  return [aux[i,:] for i=1:size(aux,1)]

end


text:
using Test

include("game.jl")

# Julia 1.0 compat
# The function definition of eachrow is taken from Julia Base,
# released under the MIT license: https://julialang.org/license
if VERSION < v"1.1"
    @eval eachrow(A) = (view(A, i, :) for i in axes(A, 1))
end

@testset "fast attack" begin
    @test !can_do_fast_attack(true)
    @test  can_do_fast_attack(false)
end

@testset "spying" begin
    character_state_combinations = Bool[
        0 0 0 0;
        0 0 1 1;
        0 1 0 1;
        0 1 1 1;
        1 0 0 1;
        1 0 1 1;
        1 1 1 1;
    ]

    for state in eachrow(character_state_combinations)
        @test can_spy(state[1:3]...) == state[4]
    end
end

@testset "signaling prisoner" begin
    character_state_combinations = Bool[
        0 0 0;
        0 1 1;
        1 0 0;
        1 1 0;
    ]

    for state in eachrow(character_state_combinations)
        @test can_signal_prisoner(state[1:2]...) == state[3]
    end
end

@testset "freeing prisoner" begin
    character_state_combinations = Bool[
        0 0 0 0 0;
        0 0 0 1 1;
        0 0 1 0 1;
        0 0 1 1 1;
        0 1 0 0 0;
        0 1 0 1 0;
        0 1 1 0 0;
        0 1 1 1 0;
        1 0 0 0 0;
        1 0 0 1 1;
        1 0 1 0 0;
        1 0 1 1 1;
        1 1 0 0 0;
        1 1 0 1 0;
        1 1 1 0 0;
        1 1 1 1 0;
    ]

    for state in eachrow(character_state_combinations)
        @test can_free_prisoner(state[1:4]...) == state[5]
    end
end


text:
 @printf("\nTesting linear operators\n")

######## Conv ############
n,m = 5, 6
h = randn(m)
op = Conv(Float64,(n,),h)
x1 = randn(n)
y1 = test_op(op, x1, randn(n+m-1), verb)
y2 = conv(x1,h)

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
op = Conv(x1,h)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

########## DCT ############
n = 4
op = DCT(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, randn(n), verb)
y2 = dct(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
op = DCT((n,))
op = DCT(n,n)
op = DCT(Complex{Float64}, n,n)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == true
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

m = 10
op = DCT(n,m)
x1 = randn(n,m)

@test norm(op'*(op*x1) - x1) <= 1e-12
@test diag_AAc(op) == 1.
@test diag_AcA(op) == 1.

######### IDCT ############
n = 4
op = IDCT(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, randn(n), verb)
y2 = idct(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
op = IDCT((n,))
op = IDCT(n,n)
op = IDCT(Complex{Float64}, n,n)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == true
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

m = 10
op = IDCT(n,m)
x1 = randn(n,m)

@test norm(op'*(op*x1) - x1) <= 1e-12
@test diag_AAc(op) == 1.
@test diag_AcA(op) == 1.

######## DFT ############
# seems like there is an object called DFT in Base julia 0.7 (however in 1.0 was rm)
n,m = 4,7

op = AbstractOperators.DFT(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = fft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Complex{Float64},(n,))
x1 = randn(n)+im*randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = fft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = fft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Complex{Float64},(n,))
x1 = randn(n)+im*randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = fft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Float64,(n,),1)
x1 = randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = fft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Complex{Float64},(n,),1)
x1 = randn(n)+im*randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = fft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Float64,(n,m))
x1 = randn(n,m)
y1 = test_op(op, x1, fft(randn(n,m)), verb)
y2 = fft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Complex{Float64},(n,m))
x1 = randn(n,m)+im*randn(n,m)
y1 = test_op(op, x1, fft(randn(n,m)), verb)
y2 = fft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Float64,(m,n),1)
x1 = randn(m,n)
y1 = test_op(op, x1, fft(randn(m,n)), verb)
y2 = fft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.DFT(Complex{Float64},(n,m),2)
x1 = randn(n,m)+im*randn(n,m)
y1 = test_op(op, x1, fft(randn(n,m)), verb)
y2 = fft(x1,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
op = AbstractOperators.DFT((n,))
op = AbstractOperators.DFT(n,n)
op = AbstractOperators.DFT(Complex{Float64}, n,n)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == false
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

m = 10
op = AbstractOperators.DFT(n,m)
x1 = randn(n,m)
y1 = op*x1
@test norm(op'*(op*x1) - diag_AcA(op)*x1) <= 1e-12
@test norm(op*(op'*y1) - diag_AAc(op)*y1) <= 1e-12

######### IDFT ############
n,m = 5,6

op = IDFT(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = ifft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Complex{Float64},(n,))
x1 = randn(n)+im*randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = ifft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Float64,(n,),1)
x1 = randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = ifft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Complex{Float64},(n,),1)
x1 = randn(n)+im*randn(n)
y1 = test_op(op, x1, fft(randn(n)), verb)
y2 = ifft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = AbstractOperators.IDFT(Float64,(n,m))
x1 = randn(n,m)
y1 = test_op(op, x1, fft(randn(n,m)), verb)
y2 = ifft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Complex{Float64},(n,m))
x1 = randn(n,m)+im*randn(n,m)
y1 = test_op(op, x1, fft(randn(n,m)), verb)
y2 = ifft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Float64,(m,n),1)
x1 = randn(m,n)
y1 = test_op(op, x1, fft(randn(m,n)), verb)
y2 = ifft(x1,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Complex{Float64},(m,n),2)
x1 = randn(m,n)+im*randn(m,n)
y1 = test_op(op, x1, fft(randn(m,n)), verb)
y2 = ifft(x1,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 4,19,5
op = IDFT(Complex{Float64},(n,m,l),2)
x1 = fft(randn(n,m,l),2)
y1 = test_op(op, x1, ifft(x1,2), verb)
y2 = ifft(x1,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 4,18,5
op = IDFT(Complex{Float64},(n,m,l),(1,2))
x1 = fft(randn(n,m,l),(1,2))
y1 = test_op(op, x1, ifft(x1,(1,2)), verb)
y2 = ifft(x1,(1,2))

@test all(norm.(y1 .- y2) .<= 1e-12)

op = IDFT(Complex{Float64},(n,m,l),(3,2))
x1 = fft(randn(n,m,l),(3,2))
y1 = test_op(op, x1, ifft(x1,(3,2)), verb)
y2 = ifft(x1,(3,2))

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
op = IDFT((n,))
op = IDFT(n,n)
op = IDFT(Complex{Float64}, n,n)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == false
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

m = 10
op = IDFT(n,m)
x1 = randn(n,m)
y1 = op*x1
@test norm(op'*(op*x1) - diag_AcA(op)*x1) <= 1e-12
@test norm(op*(op'*y1) - diag_AAc(op)*y1) <= 1e-12

####### RDFT ############
n = 4
op = RDFT(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, rfft(x1), verb)
y2 = rfft(x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 4,8,5
op = RDFT(Float64,(n,m,l),2)
x1 = randn(n,m,l)
y1 = test_op(op, x1, rfft(x1,2), verb)
y2 = rfft(x1,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
op = RDFT((n,))
op = RDFT(n,n)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == false

####### IRDFT ############
n = 12
op = IRDFT(Complex{Float64},(div(n,2)+1,),n)
x1 = rfft(randn(n))
y1 = test_op(op, x1,irfft(randn(div(n,2)+1),n), verb)
y2 = irfft(x1,n)

@test all(norm.(y1 .- y2) .<= 1e-12)

n = 11
op = IRDFT(Complex{Float64},(div(n,2)+1,),n)
x1 = rfft(randn(n))
y1 = test_op(op, x1,irfft(randn(div(n,2)+1),n), verb)
y2 = irfft(x1,n)

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 4,19,5
op = IRDFT(Complex{Float64},(n,div(m,2)+1,l),m,2)
x1 = rfft(randn(n,m,l),2)
y1 = test_op(op, x1, irfft(x1,m,2), verb)
y2 = irfft(x1,m,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 4,18,5
op = IRDFT(Complex{Float64},(n,div(m,2)+1,l),m,2)
x1 = rfft(randn(n,m,l),2)
y1 = test_op(op, x1, irfft(x1,m,2), verb)
y2 = irfft(x1,m,2)

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 5,18,5
op = IRDFT(Complex{Float64},(div(n,2)+1,m,l),n,1)
x1 = rfft(randn(n,m,l),1)
y1 = test_op(op, x1, irfft(x1,n,1), verb)
y2 = irfft(x1,n,1)

@test all(norm.(y1 .- y2) .<= 1e-12)

## other constructors
op = IRDFT((10,),19)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == false

########## DiagOp ############
n = 4
d = randn(n)
op = DiagOp(Float64,(n,),d)
x1 = randn(n)
y1 = test_op(op, x1, randn(n), verb)
y2 = d.*x1

@test all(norm.(y1 .- y2) .<= 1e-12)

n = 4
d = randn(n)+im*randn(n)
op = DiagOp(Float64,(n,),d)
x1 = randn(n)
y1 = test_op(op, x1, randn(n).+im*randn(n), verb)
y2 = d.*x1

@test all(norm.(y1 .- y2) .<= 1e-12)

n = 4
d = pi
op = DiagOp(Float64,(n,),d)
x1 = randn(n)
y1 = test_op(op, x1, randn(n), verb)
y2 = d.*x1

@test all(norm.(y1 .- y2) .<= 1e-12)

n = 4
d = im
op = DiagOp(Float64,(n,),d)
x1 = randn(n)
y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
y2 = d.*x1

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
d = randn(4)
op = DiagOp(d)

d = randn(4).+im
op = DiagOp(d)

n = 4
d = pi
op = DiagOp((n,), d)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == true
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == false
@test is_invertible(DiagOp(ones(10)))       == true
@test is_invertible(DiagOp([ones(5);0]))    == false
@test is_full_row_rank(op)    == true
@test is_full_row_rank(DiagOp([ones(5);0]))    == false
@test is_full_column_rank(op) == true
@test is_full_column_rank(DiagOp([ones(5);0]))    == false

@test diag(op) == d
@test norm(op'*(op*x1) .- diag_AcA(op).*x1) <= 1e-12
@test norm(op*(op'*x1) .- diag_AAc(op).*x1) <= 1e-12

n = 4
d = pi
op = DiagOp((n,), d)
x1 = randn(n)

@test diag(op) == d
@test norm(op'*(op*x1) .- diag_AcA(op).*x1) <= 1e-12
@test norm(op*(op'*x1) .- diag_AAc(op).*x1) <= 1e-12

########## Eye ############
n = 4
op = Eye(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, randn(n), verb)

@test all(norm.(y1 .- x1) .<= 1e-12)

# other constructors
op = Eye(Float64, (n,))
op = Eye((n,))
op = Eye(n)
op = Eye(x1)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == true
@test is_diagonal(op)         == true
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == true
@test is_invertible(op)       == true
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

@test diag(op) == 1
@test diag_AcA(op) == 1
@test diag_AAc(op) == 1

######## Filt ############
n,m = 15,2
b,a = [1.;0.;1.;0.;0.], [1.;1.;1.]
op = Filt(Float64,(n,),b,a)
x1 = randn(n)
y1 = test_op(op, x1, randn(n), verb)
y2 = filt(b, a, x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

h = randn(10)
op = Filt(Float64,(n,m),h)
x1 = randn(n,m)
y1 = test_op(op, x1, randn(n,m), verb)
y2 = filt(h, [1.], x1)

@test all(norm.(y1 .- y2) .<= 1e-12)

# other constructors
Filt(n,  b, a)
Filt((n,m),  b, a)
Filt(n,  h)
Filt((n,),  h)
Filt(x1, b, a)
Filt(x1, b)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

###### FiniteDiff ############
n= 10
op = FiniteDiff(Float64,(n,))
x1 = randn(n)
y1 = test_op(op, x1, randn(n-1), verb)
y1 = op*collect(range(0,stop=1,length=n))
@test all(norm.(y1 .- 1/9) .<= 1e-12)

I1, J1, V1 = SparseArrays.spdiagm_internal(0 => ones(n-1))
I2, J2, V2 = SparseArrays.spdiagm_internal(1 => ones(n-1))
B = -sparse(I1,J1,V1,n-1,n)+sparse(I2,J2,V2,n-1,n)

@test norm(B*x1-op*x1) <= 1e-8

n,m= 10,5
op = FiniteDiff(Float64,(n,m))
x1 = randn(n,m)
y1 = test_op(op, x1, randn(n-1,m), verb)
y1 = op*repeat(collect(range(0,stop=1,length=n)),1,m)
@test all(norm.(y1 .- 1/9) .<= 1e-12)

n,m= 10,5
op = FiniteDiff(Float64,(n,m),2)
x1 = randn(n,m)
y1 = test_op(op, x1, randn(n,m-1), verb)
y1 = op*repeat(collect(range(0,stop=1,length=n)),1,m)
@test all(norm.(y1) .<= 1e-12)

n,m,l= 10,5,7
op = FiniteDiff(Float64,(n,m,l))
x1 = randn(n,m,l)
y1 = test_op(op, x1, randn(n-1,m,l), verb)
y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l)
@test all(norm.(y1 .- 1/9) .<= 1e-12)

n,m,l= 10,5,7
op = FiniteDiff(Float64,(n,m,l),2)
x1 = randn(n,m,l)
y1 = test_op(op, x1, randn(n,m-1,l), verb)
y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l)
@test all(norm.(y1) .<= 1e-12)

n,m,l= 10,5,7
op = FiniteDiff(Float64,(n,m,l),3)
x1 = randn(n,m,l)
y1 = test_op(op, x1, randn(n,m,l-1), verb)
y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l)
@test all(norm.(y1) .<= 1e-12)

n,m,l,i = 5,6,2,3
op = FiniteDiff(Float64,(n,m,l,i),1)
x1 = randn(n,m,l,i)
y1 = test_op(op, x1, randn(n-1,m,l,i), verb)
y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l*i),n,m,l,i)
@test all(norm.(y1 .- 1/(n-1)) .<= 1e-12)

n,m,l,i = 5,6,2,3
op = FiniteDiff(Float64,(n,m,l,i),4)
x1 = randn(n,m,l,i)
y1 = test_op(op, x1, randn(n,m,l,i-1), verb)
y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l*i),n,m,l,i)
@test norm(y1) <= 1e-12

@test_throws ErrorException FiniteDiff(Float64,(n,m,l), 4)

## other constructors
FiniteDiff((n,m))
FiniteDiff(x1)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == false

########## GetIndex ############
n,m = 5,4
k = 3
op = GetIndex(Float64,(n,),(1:k,))
x1 = randn(n)
y1 = test_op(op, x1, randn(k), verb)

@test all(norm.(y1 .- x1[1:k]) .<= 1e-12)

n,m = 5,4
k = 3
op = GetIndex(Float64,(n,m),(1:k,:))
x1 = randn(n,m)
y1 = test_op(op, x1, randn(k,m), verb)

@test all(norm.(y1 .- x1[1:k,:]) .<= 1e-12)

n,m = 5,4
op = GetIndex(Float64,(n,m),(:,2))
x1 = randn(n,m)
y1 = test_op(op, x1, randn(n), verb)

@test all(norm.(y1 .- x1[:,2]) .<= 1e-12)

n,m,l = 5,4,3
op = GetIndex(Float64,(n,m,l),(1:3,2,:))
x1 = randn(n,m,l)
y1 = test_op(op, x1, randn(3,3), verb)

@test all(norm.(y1 .- x1[1:3,2,:]) .<= 1e-12)

# other constructors
GetIndex((n,m), (1:k,:))
GetIndex(x1, (1:k,:))

@test_throws ErrorException GetIndex(Float64,(n,m),(1:k,:,:))
op = GetIndex(Float64,(n,m),(1:n,1:m))
@test typeof(op) <: Eye

op = GetIndex(Float64,(n,),(1:k,))

##properties
@test is_sliced(op)           == true
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == true
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == false

@test diag_AAc(op) == 1

####### MatrixOp ############

# real matrix, real input
n,m = 5,4
A = randn(n,m)
op = MatrixOp(Float64,(m,),A)
x1 = randn(m)
y1 = test_op(op, x1, randn(n), verb)
y2 = A*x1

# real matrix, complex input
n,m = 5,4
A = randn(n,m)
op = MatrixOp(Complex{Float64},(m,),A)
x1 = randn(m)+im.*randn(m)
y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
y2 = A*x1

# complex matrix, complex input
n,m = 5,4
A = randn(n,m)+im*randn(n,m)
op = MatrixOp(Complex{Float64},(m,),A)
x1 = randn(m)+im.*randn(m)
y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
y2 = A*x1

# complex matrix, real input
n,m = 5,4
A = randn(n,m)+im*randn(n,m)
op = MatrixOp(Float64,(m,),A)
x1 = randn(m)
y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
y2 = A*x1

@test all(norm.(y1 .- y2) .<= 1e-12)

# complex matrix, real matrix input
c = 3
op = MatrixOp(Float64,(m,c),A)
@test_throws ErrorException MatrixOp(Float64,(m,c,3),A)
@test_throws MethodError MatrixOp(Float64,(m,c),randn(n,m,2))
x1 = randn(m,c)
y1 = test_op(op, x1, randn(n,c).+randn(n,c), verb)
y2 = A*x1

# other constructors
op = MatrixOp(A)
op = MatrixOp(Float64, A)
op = MatrixOp(A, c)
op = MatrixOp(Float64, A, c)

op = convert(LinearOperator,A)
op = convert(LinearOperator,A,c)
op = convert(LinearOperator, Complex{Float64}, size(x1), A)

##properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(MatrixOp(randn(Random.seed!(0),3,4)))    == true
@test is_full_column_rank(MatrixOp(randn(Random.seed!(0),3,4))) == false

##### LMatrixOp ############

n,m = 5,6
b = randn(m)
op = LMatrixOp(Float64,(n,m),b)
x1 = randn(n,m)
y1 = test_op(op, x1, randn(n), verb)
y2 = x1*b

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m = 5,6
b = randn(m)+im*randn(m)
op = LMatrixOp(Complex{Float64},(n,m),b)
x1 = randn(n,m)+im*randn(n,m)
y1 = test_op(op, x1, randn(n)+im*randn(n), verb)
y2 = x1*b

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 5,6,7
b = randn(m,l)
op = LMatrixOp(Float64,(n,m),b)
x1 = randn(n,m)
y1 = test_op(op, x1, randn(n,l), verb)
y2 = x1*b

@test all(norm.(y1 .- y2) .<= 1e-12)

n,m,l = 5,6,7
b = randn(m,l)+im*randn(m,l)
op = LMatrixOp(Complex{Float64},(n,m),b)
x1 = randn(n,m)+im*randn(n,m)
y1 = test_op(op, x1, randn(n,l)+im*randn(n,l), verb)
y2 = x1*b

@test all(norm.(y1 .- y2) .<= 1e-12)

## other constructors
op = LMatrixOp(b,n)

##properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == false
@test is_full_column_rank(op) == false

######### MyLinOp ############

n,m = 5,4;
A = randn(n,m);
op = MyLinOp(Float64, (m,),(n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x))
x1 = randn(m)
y1 = test_op(op, x1, randn(n), verb)
y2 = A*x1

# other constructors
op = MyLinOp(Float64, (m,), Float64, (n,), (y,x) -> mul!(y,A,x), (y,x) -> mul!(y,A',x))

######### MIMOFilt ############
m,n = 10,2
b = [[1.;0.;1.;0.;0.],[1.;0.;1.;0.;0.]]
a = [[1.;1.;1.],[2.;2.;2.]]
op = MIMOFilt(Float64, (m,n), b, a)

x1 = randn(m,n)
y1 = test_op(op, x1, randn(m,1), verb)
y2 = filt(b[1],a[1],x1[:,1])+filt(b[2],a[2],x1[:,2])

@test all(norm.(y1 .- y2) .<= 1e-12)

m,n = 10,3; #time samples, number of inputs
b = [[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.],[1.;0.;1.], ];
a = [[1.;1.;1.],[2.;2.;2.],[      3.],[      4.],[      5.],[      6.], ];
op = MIMOFilt(Float64, (m,n), b, a)

x1 = randn(m,n)
y1 = test_op(op, x1, randn(m,2), verb)
y2 = [filt(b[1],a[1],x1[:,1])+filt(b[2],a[2],x1[:,2])+filt(b[3],a[3],x1[:,3]) filt(b[4],a[4],x1[:,1])+filt(b[5],a[5],x1[:,2])+filt(b[6],a[6],x1[:,3])]

@test all(norm.(y1 .- y2) .<= 1e-12)

m,n = 10,3
b = [randn(10),randn(5),randn(10),randn(2),randn(10),randn(10)]
a = [[1.],[1.],[1.],[1.],[1.],[1.]]
op = MIMOFilt(Float64, (m,n), b, a)

x1 = randn(m,n)
y1 = test_op(op, x1, randn(m,2), verb)
y2 = [filt(b[1],a[1],x1[:,1])+filt(b[2],a[2],x1[:,2])+filt(b[3],a[3],x1[:,3]) filt(b[4],a[4],x1[:,1])+filt(b[5],a[5],x1[:,2])+filt(b[6],a[6],x1[:,3])]

@test all(norm.(y1 .- y2) .<= 1e-12)

## other constructors
MIMOFilt((10,3),  b, a)
MIMOFilt((10,3),  b)
MIMOFilt(x1,  b, a)
MIMOFilt(x1,  b)

#errors
@test_throws ErrorException MIMOFilt(Float64, (10,3,2) ,b,a)
a2 = [[1.0f0],[1.0f0],[1.0f0],[1.0f0],[1.0f0],[1.0f0]]
b2 = convert.(Array{Float32,1},b)
@test_throws ErrorException MIMOFilt(Float64, (m,n),b2,a2)
@test_throws ErrorException MIMOFilt(Float64, (m,n),b,a[1:end-1])
push!(a2,[1.0f0])
push!(b2,randn(Float32,10))
@test_throws ErrorException MIMOFilt(Float32, (m,n),b2,a2)
a[1][1] = 0.
@test_throws ErrorException MIMOFilt(Float64, (m,n) ,b,a)

b = [randn(10),randn(5),randn(10),randn(2),randn(10),randn(10)]
a = [[1.],[1.],[1.],[1.],[1.],[1.]]
op = MIMOFilt(Float64, (m,n), b, a)

##properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

######## Variation ############

n,m = 10,5
op = Variation(Float64,(n,m))
x1 = randn(n,m)
y1 = test_op(op, x1, randn(m*n,2), verb)

y1 = op*repeat(collect(range(0,stop=1,length=n)),1,m)
@test all(norm.(y1[:,1] .- 1/(n-1) ) .<= 1e-12)
@test all(norm.(y1[:,2] ) .<= 1e-12)

Dx = spdiagm(0 => ones(n), -1 => -ones(n-1))
Dx[1,1],Dx[1,2] = -1,1
Dy = spdiagm(0 => ones(m), -1 => -ones(m-1))
Dy[1,1],Dy[1,2] = -1,1

Dxx = kron(sparse(I,m,m),Dx)
Dyy = kron(Dy,sparse(I,n,n))
TV = [Dxx;Dyy]

x1 = randn(n,m)
@test norm(op*x1-reshape(TV*(x1[:]),n*m,2))<1e-12

n,m,l = 10,5,3
op = Variation(Float64,(n,m,l))
x1 = randn(n,m,l)
y1 = test_op(op, x1, randn(m*n*l,3), verb)

y1 = op*reshape(repeat(collect(range(0,stop=1,length=n)),1,m*l),n,m,l)
@test all(norm.(y1[:,1] .- 1/(n-1) ) .<= 1e-12)
@test all(norm.(y1[:,2] ) .<= 1e-12)
@test all(norm.(y1[:,3] ) .<= 1e-12)

### other constructors
Variation(Float64, n,m)
Variation((n,m))
Variation(n,m)
Variation(x1)

##errors
@test_throws ErrorException Variation(Float64,(n,))

###properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == false
@test is_full_column_rank(op) == false

######## Xcorr ############
n,m = 5, 6
h = randn(m)
op = Xcorr(Float64,(n,),h)
x1 = randn(n)
y1 = test_op(op, x1, randn(n+m), verb)
y2 = xcorr(x1, h; padmode=:longest)

@test all(norm.(y1 .- y2) .<= 1e-12)
# other constructors
op = Xcorr(x1,h)

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == true
@test is_full_column_rank(op) == true

######## ZeroPad ############
n = (3,)
z = (5,)
op = ZeroPad(Float64,n,z)
x1 = randn(n)
y1 = test_op(op, x1, randn(n.+z), verb)
@test all(norm.(y1 .- [x1;zeros(5)] ) .<= 1e-12)

n = (3,2)
z = (5,3)
op = ZeroPad(Float64,n,z)
x1 = randn(n)
y1 = test_op(op, x1, randn(n.+z), verb)
y2 = zeros(n.+z)
y2[1:n[1],1:n[2]] = x1
@test all(norm.(y1 .- y2 ) .<= 1e-12)

n = (3,2,2)
z = (5,3,1)
op = ZeroPad(Float64,n,z)
x1 = randn(n)
y1 = test_op(op, x1, randn(n.+z), verb)
y2 = zeros(n.+z)
y2[1:n[1],1:n[2],1:n[3]] = x1
@test all(norm.(y1 .- y2 ) .<= 1e-12)

# other constructors
ZeroPad(n, z...)
ZeroPad(Float64, n, z...)
ZeroPad(n, z...)
ZeroPad(x1, z)
ZeroPad(x1, z...)

#errors
@test_throws ErrorException ZeroPad(Float64,n,(1,2))
@test_throws ErrorException ZeroPad(Float64,n,(1,-2,3))

#properties
@test is_linear(op)           == true
@test is_null(op)             == false
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == true
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == false
@test is_full_column_rank(op) == true

diag_AcA(op) == 1

########## Zeros ############
n = (3,4)
D = Float64
m = (5,2)
C = Complex{Float64}
op = Zeros(D,n,C,m)
x1 = randn(n)
y1 = test_op(op, x1, randn(m)+im*randn(m), verb)

#properties
@test is_linear(op)           == true
@test is_null(op)             == true
@test is_eye(op)              == false
@test is_diagonal(op)         == false
@test is_AcA_diagonal(op)     == false
@test is_AAc_diagonal(op)     == false
@test is_orthogonal(op)       == false
@test is_invertible(op)       == false
@test is_full_row_rank(op)    == false
@test is_full_column_rank(op) == false


text:
"""
	IterableWrapper(iter)

Wraps a Iterable object
"""

import Base: iterate, eltype, length, size, IteratorEltype, IteratorSize
export IterableTraitWrapper#, IteratorEltype, IteratorSize, iterate, eltype, length, size

struct IterableTraitWrapper{T} <: AbstractTraitWrapper
	object::T
	IterableTraitWrapper(t::T) where T = new{T}(t)
end

IteratorEltype(t::IterableTraitWrapper) = IteratorEltype(object(t))
IteratorSize(t::IterableTraitWrapper) = IteratorSize(object(t))
iterate(t::IterableTraitWrapper) = iterate(object(t))
iterate(t::IterableTraitWrapper, state) = iterate(object(t), state)
eltype(t::IterableTraitWrapper) = eltype(object(t))
length(t::IterableTraitWrapper) = length(object(t))
size(t::IterableTraitWrapper) = size(object(t))


text:
module TransportationMetaheuristics

end # module


text:
@testset "Clustering metrics" begin
    x = rand(50)

    for ps_type in (FearnheadParticles, ChenLiuParticles)
        @testset "$ps_type" begin
            ps = ps_type(100, NormalInverseChisq(0., 1, 2, 2), ChineseRestaurantProcess(1.))
            filter!(ps, x, false)

            asgn = assignments(ps)
            @test size(asgn) == (50, 100)
            @test eltype(asgn) == Int
            @test vec(maximum(asgn, dims=1)) == ncomponents.(ps.particles)
            @test all(minimum(asgn, dims=1) .== 1)

            asgn_sim = assignment_similarity(ps)
            @test size(asgn_sim) == (50, 50)
            @test issymmetric(asgn_sim)
            @test all(0 .≤ asgn_sim .≤ 1)
            @test all(diag(asgn_sim) .== 1)
        end
    end
end


text:
function trace!(tr, d, state, iteration, method::ZerothOrderOptimizer, options)
    dt = Dict()
    if options.extended_trace
        dt["x"] = copy(state.x)
    end
    update!(tr,
            state.iteration,
            d.F,
            NaN,
            dt,
            options.store_trace,
            options.show_trace,
            options.show_every,
            options.callback)
end

function assess_convergence(state::ZerothOrderState, d, options)
    false, false, false, false, false
end

f_abschange(d::AbstractObjective, state::ZerothOrderState) = convert(typeof(value(d)), NaN)
x_abschange(state::ZerothOrderState) = convert(real(eltype(state.x)), NaN)


text:
# Included into SamplingMethods.jl
#
# using SamplingMethods: CuSVTuple, SamplingMethod, SamplingMethodState, Snapshot

"""    ABF

Adaptive Biasing Force
"""
Base.@kwdef struct ABF{CV, G} <: SamplingMethod
    ξ::CV
    grid::G
    n::Int = 200
end

(self::ABF)(snapshot) = ABFState(self, snapshot)

struct ABFState{S, CV, G, B, H, U, V} <: SamplingMethodState
    snapshot::S
    ξ::CV
    grid::G
    bias::B
    hist::H
    ΣF::U
    F::V
    Wp::V
    Wp₋::V
    n::Int

    function ABFState(sampling_method::ABF, snapshot::S) where {S}
        grid = sampling_method.grid
        N = ndims(grid)
        shape = size(grid)

        ξ = DynamicCV(sampling_method.ξ, snapshot.positions, snapshot.ids)
        bias = fill!(similar(snapshot.forces), 0)
        hist = fill!(similar(snapshot.ids, shape), 0)
        ΣF = fill!(similar(snapshot.forces, (N, shape...)), 0)
        F = fill!(similar(snapshot.forces, N), 0)
        Wp = fill!(similar(snapshot.forces, N), 0)
        Wp₋ = fill!(similar(snapshot.forces, N), 0)

        return new{
            S, typeof(ξ), typeof(grid), typeof(bias), typeof(hist), typeof(ΣF), typeof(F)
        }(
            snapshot, ξ, grid, bias, hist, ΣF, F, Wp, Wp₋, sampling_method.n
        )
    end
end

(sm::ABFState)(sync_backend = nothing) = update!(sm; sync_backend = sync_backend)
#
function (sm::ABFState{<: Snapshot{<: CuSVTuple}})(; sync_backend = nothing)
    # When on GPU, make sure that forces are synchronized after each update.
    return CUDA.@sync update!(sm; sync_backend = sync_backend)
end

function update!(sm::ABFState{<: Snapshot}; sync_backend = nothing)
    ss = sm.snapshot
    rs = ss.positions
    vs = ss.velocities
    ms = ss.masses
    dt = ss.dt

    # Compute the collective variable and its Jacobian
    ξ, Jξ = sm.ξ(rs)
    # Compute momenta
    p = vec(ms' .* vs)
    #Wp₊ = Jξ \ p
    Wp₊ = ldiv(Jξ, p)
    # Second order backward finite difference
    Ẇp = (3 .* Wp₊ .- 4 .* sm.Wp .+ sm.Wp₋) ./ 2dt

    I_ξ = sm.grid[ξ]
    N_ξ = view(sm.hist, I_ξ)
    F_ξ = view(sm.ΣF, :, I_ξ)

    N_ξ .+= 1
    F_ξ .+= Ẇp .+ sm.F
    sm.F .= F_ξ ./ max.(N_ξ, sm.n)
    sm.Wp₋ .= sm.Wp
    sm.Wp .= Wp₊

    # Compute bias
    # TODO: Optimize with `mul!`
    sm.bias[:] .= reshape(Jξ, :, ndims(sm.grid)) * sm.F

    # TODO: Synchronization and biasing should be done outside
    # `update!(::SamplingMethodState)` so that if there are many sampling methods,
    # they can be run independently.

    # Forces may be computed asynchronously on the GPU, so we need to
    # synchronize them before applying the bias.
    if sync_backend !== nothing
        sync_backend()
    end
    # Add bias
    ss.forces .-= sm.bias

    return nothing
end


text:
function mirror(image, mask)
    ixs = map((ix, m) -> m ? reverse(ix) : ix, axes(image), mask)
    view(image, ixs...)
end

"""
    Params_map

Store parameters and inner pre-allocated matrices.
"""
struct Params_map{ImgArray,ResultArray}
    mirror_img::ImgArray
    mirror_result::ResultArray
end

"""
    Params_map(img, rmask)

Return parameters that can be used sequentially for all similar imgs.
"""
function Params_map(img)
    mirror_img = similar(img)
    mirror_result = similar(img, Float32)

    Params_map(
        mirror_img,
        mirror_result
    )
end

function corr_function_map!(
    result::CFMap{T1, N},
    img::AbstractArray{T2,N},
    map_params::Params_map,
    cf_params
) where T1 where T2 where N
    mirror_img = map_params.mirror_img
    mirror_result = map_params.mirror_result

    for mask in BoolMask(N)
        if result.cf_type == :central_symmetry && mask[end]
            continue
        elseif result.cf_type == :periodic_point_point && any(mask)
            continue
        end

        mirror_img .= mirror(img, mask)
        mirror_result .= 0

        correllation_function!(mirror_result, mirror_img, cf_params)

        v_result = cut_cfmap(result, mask)
        v_result .= mirror_result
    end
    result
end

"""
    corr_function_map(img, CF_constructor; parameters...)

Compute correlation function map in all meaningfull directions.

CF_constructor ∈ {Map.s2, Map.l2, Map.c2, Map.ss, Map.sv}
"""
function corr_function_map(
    img::AbstractArray{T,N},
    CF_constructor;
    parameters...
) where T where N
    cf_params, cf_type = CF_constructor(img; parameters...)

    map_params = Params_map(img)
    result = CFMap(img, cf_type)

    corr_function_map!(result, img, map_params, cf_params)
end