text:
# ===========
# Indentation
# ===========

mutable struct Indentation
    ind::String
    step::String

    # default indentation: two spaces
    function Indentation(ind="", step="  ")
        return new(ind, step)
    end
end

function indent!(ind::Indentation)
    ind.ind = ind.ind * ind.step
end

function dedent!(ind::Indentation)
    ind.ind = ind.ind[1:end-(length(ind.step))]
end

function _write_indented(io, string, indentation)
    write(io, "$(indentation.ind)$string")
end

function _dedent(indentation)
    return
end

function _variable_name(id, dictionary)
    get(dictionary, id, "x" * string(id))
end

function _input_name(id, dictionary)
    get(dictionary, id, "u" * string(id))
end

# ======================
# Creation of model file
# ======================

function writesxmodel(filename, system; dictionary=Dict())
    open(filename, "w") do io
        write(io,
"<?xml version=\"1.0\" encoding=\"iso-8859-1\"?>\n" *
"<sspaceex xmlns=\"http://www-verimag.imag.fr/xml-namespaces/sspaceex\" version=\"0.2\" math=\"SpaceEx\">\n"
        )
        indentation = Indentation()
        indent!(indentation)
        _write_system(io, system, dictionary, indentation)
        dedent!(indentation)
        write(io, "</sspaceex>\n")
    end
    nothing
end

function _write_system(io, system::AbstractContinuousSystem, dictionary,
                       indentation)
    _write_component(io, system, dictionary, indentation)
end

function _write_component(io, system, dictionary, indentation)
    _write_indented(io, "<component id=\"system\">\n", indentation)
    indent!(indentation)
    _write_parameters(io, system, dictionary, indentation)
    _write_locations(io, system, dictionary, indentation)
    _write_transitions(io, system, dictionary, indentation)
    dedent!(indentation)
    _write_indented(io, "</component>\n", indentation)
end

function _write_parameters(io,
                           system::Union{LinearContinuousSystem,
                                         ConstrainedLinearContinuousSystem},
                           dictionary, indentation)
    _write_linear_system(io, system, dictionary, indentation)
end

function _write_parameters(io,
                           system::Union{LinearControlContinuousSystem,
                                         ConstrainedLinearControlContinuousSystem},
                           dictionary, indentation)
    _write_linear_system(io, system, dictionary, indentation)

    m = inputdim(system)
    controlled = false
    for i in 1:m
        name = _input_name(i, dictionary)
        _write_parameter(io, name, controlled, dictionary, indentation)
    end
end

function _write_linear_system(io, system::AbstractContinuousSystem,
                              dictionary, indentation)
    n = statedim(system)
    controlled = true
    for i in 1:n
        name = _variable_name(i, dictionary)
        _write_parameter(io, name, controlled, dictionary, indentation)
    end
end

function _write_parameter(io, name, controlled, dictionary, indentation)
    _write_indented(io, "<param name=\"$name\" type=\"real\" d1=\"1\" d2=\"1\" " *
        "local=\"false\" dynamics=\"any\" controlled=\"$controlled\" />\n",
        indentation)
end

function _write_locations(io, system::AbstractContinuousSystem, dictionary,
                          indentation)
    _write_location(io, system, 1, dictionary, indentation)
end

function _write_location(io, system, id, dictionary, indentation)
    name = "loc" * string(id)
    _write_indented(io, "<location id=\"$id\" name=\"$name\">\n", indentation)
    indent!(indentation)
    _write_invariant(io, system, dictionary, indentation)
    _write_flow(io, system, dictionary, indentation)
    dedent!(indentation)
    _write_indented(io, "</location>\n", indentation)
end

function _write_invariant(io, system, dictionary, indentation)
    if !isconstrained(system)
        # nothing to write
        return
    end

    _write_indented(io, "<invariant>", indentation)
    indent!(indentation)

    X = stateset(system)
    if X != nothing
        _write_invariant_states_specific(io, system, X, dictionary, indentation)
    end

    U = inputset(system)
    if U != nothing
        _write_invariant_inputs_specific(io, system, U, dictionary, indentation)
    end

    dedent!(indentation)
    write(io, "</invariant>\n")
end

function _write_invariant_states_specific(io, system, X, dictionary, indentation)
    println("WARNING: state constraints of type $(typeof(X)) are currenctly " *
        "not supported and will be ignored")
end

function _write_invariant_states_specific(io, system, X::Universe, dictionary,
                                          indentation)
    # nothing to write
end

function _write_invariant_states_specific(io, system, X::AbstractHyperrectangle,
                                          dictionary, indentation)
    for i in 1:dim(X)
        xi = _variable_name(i, dictionary)
        l = low(X, i)
        u = high(X, i)
        if i > 1
            write(io, " & ")
        end
        write(io, "$l <= $xi & $xi <= $u")
    end
end

function _write_invariant_inputs_specific(io, system, U, dictionary, indentation)
    println("WARNING: input constraints of type $(typeof(U)) are currenctly " *
        "not supported and will be ignored")
end

function _write_invariant_inputs_specific(io, system, U::Universe, dictionary,
                                          indentation)
    # nothing to write
end

function _write_invariant_inputs_specific(io, system, U::AbstractHyperrectangle,
                                          dictionary, indentation)
    for i in 1:dim(U)
        ui = _input_name(i, dictionary)
        l = low(U, i)
        u = high(U, i)
        if i > 1
            write(io, " & ")
        end
        write(io, "$l <= $ui & $ui <= $u")
    end
end

function _write_flow(io, system, dictionary, indentation)
    _write_indented(io, "<flow>", indentation)
    _write_flow_specific(io, system, dictionary)
    write(io, "</flow>\n")
end

function _write_flow_specific(io,
                              system::Union{LinearContinuousSystem,
                                            LinearControlContinuousSystem,
                                            ConstrainedLinearContinuousSystem,
                                            ConstrainedLinearControlContinuousSystem},
                              dictionary)
    A = state_matrix(system)
    n = statedim(system)

    B = input_matrix(system)  # returns `nothing` for systems with no inputs
    m = inputdim(system)

    for i in 1:n
        xi = _variable_name(i, dictionary)
        write(io, "$xi' == ")
        first = true
        for j in 1:n
            Aij = A[i, j]
            if Aij == 0
                continue
            end
            xj = _variable_name(j, dictionary)
            if Aij < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(Aij)) * $xj")
            first = false
        end
        for j in 1:m
            Bij = B[i, j]
            if Bij == 0
                continue
            end
            uj = _input_name(j, dictionary)
            if Bij < 0
                if first
                    prefix = "-"
                else
                    prefix = " - "
                end
            else
                if first
                    prefix = ""
                else
                    prefix = " + "
                end
            end
            write(io, "$prefix$(abs(Bij)) * $uj")
            first = false
        end
        if first
            # rhs was zero everywhere
            write(io, "0")
        end
        if i < n
            write(io, " & ")
        end
    end
end

function _write_transitions(io, system::AbstractContinuousSystem, dictionary,
                            indentation)
    # nothing to write
end


text:
# The PUG XML interface is much more painful, but it can run longer queries

"""
    cids = query_substructure_pug(;cid=nothing, smiles=nothing, smarts=nothing,        # specifier for the substructure to search for
                                   maxhits=200_000, poll_interval=10)

Retrieve a list of compounds containing a substructure specified via its `cid`, the SMILES string, or a SMARTS string.

# Example

```julia
julia> using PubChemCrawler

julia> cids = query_substructure_pug(smarts="[r13]Br")   # query brominated 13-atom rings
66-element Vector{$Int}:
  54533707
 153064026
 152829033
...
```

PUG searches can take a while to run (they poll for completion), but conversely they allow more complex, long-running searches
to succeed. See also [`query_substructure`](@ref).
"""
function query_substructure_pug(; kwargs...)
    xdoc = create_substructure_query(; kwargs...)
    cids = submit_substructure_query(xdoc; kwargs...)
    free(xdoc)
    return cids
end

# Implement a common pattern in the XML structure of PUG requests:
# <Tag1>
#   <Tag1_modifierA>
#     <Tag2>
#       <Tag2_modifierB>
# ...
function nest2(parent, tag, mod)
    nodetag = new_child(parent, tag)
    child = new_child(nodetag, tag*'_'*mod)
    return nodetag, child
end
function nest2(parent, tag, mods...)
    nodetag = new_child(parent, tag)
    children = [new_child(nodetag, tag*'_'*mod) for mod in mods]
    return nodetag, children
end

function create_substructure_query(;cid=nothing, smiles=nothing, smarts=nothing,          # inputs
    stereo="ignore",isotopes::Bool=false,charges::Bool=false,tautomers::Bool=false,rings::Bool=false,bonds::Bool=true,chains::Bool=true,hydrogen::Bool=false,
    maxhits::Int=2_000_000,
    kwargs...)

    query_data = if cid !== nothing
        (smiles === nothing && smarts === nothing) || throw(ArgumentError("only one of cid, smiles, or smarts can be specified"))
        string(cid)
    elseif smiles !== nothing
        smarts === nothing || throw(ArgumentError("only one of cid, smiles, or smarts can be specified"))
        smiles
    elseif smarts !== nothing
        smarts
    else
        throw(ArgumentError("one of cid, smiles, or smarts must be specified"))
    end
    stereoidx = 0
    if stereo != "ignore"
        stereoidx = findfirst(isequal(stereo), ["exact", "relative", "non-conflicting"])
        stereoidx === nothing && throw(ArgumentError("""stereo must be one of "ignore", "exact", "relative", or "non-conflicting", got $stereo"""))
    end

    xdoc = XMLDocument()
    xroot = create_root(xdoc, "PCT-Data")
    rootinput = new_child(xroot, "PCT-Data_input")
    input, inputquery = nest2(rootinput, "PCT-InputData", "query")
    query, querytype = nest2(inputquery, "PCT-Query", "type")
    querytype2, querytype2css = nest2(querytype, "PCT-QueryType", "css")
    querycomp, fields = nest2(querytype2css, "PCT-QueryCompoundCS", "query", "type", "results")
    add_text(new_child(fields[1], "PCT-QueryCompoundCS_query_data"), query_data)  # under "query"
    subss = new_child(fields[2], "PCT-QueryCompoundCS_type_subss")                # under "type"

    ssstr = "PCT-CSStructure"
    substruct = new_child(subss, ssstr)
    if stereoidx != 0
        cstereo = new_child(substruct, ssstr*"_stereo")
        set_attribute(cstereo, "value", stereo)
        add_text(cstereo, string(stereoidx))
    end
    for (name, val) in (("isotopes", isotopes), ("charges", charges), ("tautomers", tautomers), ("rings", rings), ("bonds", bonds), ("chains", chains), ("hydrogen", hydrogen))
        child = new_child(substruct, ssstr*'_'*name)
        set_attribute(child, "value", repr(val))
    end

    add_text(fields[3], string(maxhits))               # under "results"

    return xdoc
end

function create_poll(id)
    xdoc = XMLDocument()
    xroot = create_root(xdoc, "PCT-Data")
    rootinput = new_child(xroot, "PCT-Data_input")
    input, inputreq = nest2(rootinput, "PCT-InputData", "request")
    req, reqdata = nest2(inputreq, "PCT-Request", "reqid", "type")
    add_text(reqdata[1], id)
    set_attribute(reqdata[2], "value", "status")
    return xdoc
end

# function create_download(entrez; format="CSV", compression="gzip")
#     xdoc = XMLDocument()
#     xroot = create_root(xdoc, "PCT-Data")
#     rootinput = new_child(xroot, "PCT-Data_input")
#     input, inputreq = nest2(rootinput, "PCT-InputData", "download")
#     req, reqdata = nest2(inputreq, "PCT-Download", "uids", "format", "compression")
#     query, centrez = nest2(reqdata[1], "PCT-QueryUids", "entrez")
#     _, fields = nest2(centrez, "PCT-Entrez", "db", "query-key", "webenv")
#     add_text(fields[1], content(entrez["PCT-Entrez_db"][1]))
#     add_text(fields[2], content(entrez["PCT-Entrez_query-key"][1]))
#     add_text(fields[3], content(entrez["PCT-Entrez_webenv"][1]))
#     set_attribute(reqdata[2], "value", format)
#     set_attribute(reqdata[3], "value", compression)
#     return xdoc
# end

function get_top(xdoc)
    xdoc_root = root(xdoc)
    return xdoc_root["PCT-Data_output"][1]["PCT-OutputData"][1]
end
function get_status(top)
    elstatus = top["PCT-OutputData_status"][1]["PCT-Status-Message"][1]["PCT-Status-Message_status"][1]["PCT-Status"][1]
    return attribute(elstatus, "value")
end
function get_id(top)
    return content(top["PCT-OutputData_output"][1]["PCT-OutputData_output_waiting"][1]["PCT-Waiting"][1]["PCT-Waiting_reqid"][1])
end
function get_entrez(top)
    return top["PCT-OutputData_output"][1]["PCT-OutputData_output_entrez"][1]["PCT-Entrez"][1]
end

function submit_substructure_query(xdoc; poll_interval=10, kwargs...)
    r = HTTP.request("POST", "https://pubchem.ncbi.nlm.nih.gov/pug/pug.cgi", [], string(xdoc); kwargs...)
    xresp = parse_string(String(r.body))
    top = get_top(xresp)
    status = get_status(top)
    if status == "queued"
        # Wait for the results
        id = get_id(top)
        pollreq = create_poll(id)
        while true
            free(xresp)
            sleep(poll_interval)
            # poll for completion
            r = HTTP.request("POST", "https://pubchem.ncbi.nlm.nih.gov/pug/pug.cgi", [], string(pollreq); kwargs...)
            xresp = parse_string(String(r.body))
            top = get_top(xresp)
            status = get_status(top)
            status != "queued" && status != "running" && break
        end
        free(pollreq)
    end
    if status == "hit-limit"
        @warn "maxhits was hit, results are partial"
    elseif status != "success"
        error(string(xresp))
    end
    entrez = get_entrez(top)
    # Now retrieve the results
    url = "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/efetch.fcgi?retmode=text&rettype=uilist&WebEnvRq=1&db=pccompound"
    url *= "&query_key=" * content(entrez["PCT-Entrez_query-key"][1])
    url *= "&WebEnv=" * content(entrez["PCT-Entrez_webenv"][1])
    free(xresp)
    r = HTTP.request("GET", url)
    return parse.(Int, split(chomp(String(r.body)), '\n'))
    # xdl = create_download(entrez; format=format)
    # r = HTTP.request("POST", "https://pubchem.ncbi.nlm.nih.gov/pug/pug.cgi", [], string(xdl); kwargs...)
    # result = String(r.body)
    # free(xdl)
end

ext: jl


text:
using Zomato
using Compat.Test

# More comprehensive tests needed!
@test typeof(Zomato.authenticate("apikey")) == Zomato.Auth


text:
InputFolder = './Images/CellsPhaseContrast/';
OutputFolder = './Results/Images/CellsPhaseContrast/';
Rescale = 2;

@iA = '*.tif';

@fxg_gRestorePhaseContrast [iA] > [O];
params.w_smooth_spatio = 1;
params.w_sparsity = 0.001;
params.gamma = 0.1;
params.maxiter = 100;
params.tol = 1;
params.Rwid = 4000;		% Settings for "a" phase contrast
params.Wwid = 800;		% microscope and 1.3 um pixel size
params.MRadius = 3;		%
/endf

/show iA > ;
/show O > ;
/keep O > tif;

max_stars_repo_name: SebastienTs/LOBSTER

max_issues_repo_name: SebastienTs/LOBSTER

max_forks_repo_name: SebastienTs/LOBSTER


text:
# PeriodicTable.Elemen
using Unitful
using Pkg
using NumberIntervals
using Markdown

"""
    configuration(elm::Element)

The configuration of the shell occupancy in a specific ground-state element.
"""
configuration(elm::Element) = elm.el_config

"""
    element(z::Int)::PeriodicTable.Element

Covert an atomic number into a PeriodicTable.Element
"""
element(z::Integer) = PeriodicTable.elements[z]

"""
    a(elm::Element)

Return the mean atomic weight of the Element in amu
"""
a(elm::Element) = ustrip(elm.atomic_mass |> u"u")

"""
    z(elm::Element)

Return the atomic number of the Element.
"""
z(elm::Element) = elm.number

"""
    symbol(elm::Element)

Return the symbol like "H", "He", "Li", .... for the Element.
"""
symbol(elm::Element) = elm.symbol

"""
    name(elm::Element)

Return the name like "Hydrogen", "Helium",... for the Element.
"""
name(elm::Element) = elm.name

"""
    density(elm::Element)

Return the nominal density for the element in g/cm³.
"""
density(elm::Element) = ustrip(elm.density |> u"g/cm^3")

Base.:(==)(elm1::Element, elm2::Element) = z(elm1) == z(elm2)

function _pp(ss::String)
    p = findfirst('.', ss)
    o = findfirst('(',ss)
    c = findfirst(')',ss)
    v = parse(Float64, ss[1:o-1])
    dv = parse(Int, ss[o+1:c-1])*0.1^(o - p - 1)
    return uv(v,dv)
end

"""
    atomic_weight[elm::Element]::Union{NumberInterval, UncertainValue}

Atomic weights from the 2020 tabulation at https://ciaaw.org/atomic-weights.htm.  Not all elements are represented because not all
elements have a nominal isotopic distribution.  Some, like Tc, don't exist in nature.  Others, like Pm, are instable.  Most atomic weights
are given as `UncertainValue` instances while a few are `NumberInterval` instances.  For example, the atomic weight of Pb is highly variable
and is thus given as a range.  (See the website for more details.)
"""
const atomic_weight = Dict{Element, AbstractFloat}(
    elements[1] => NumberInterval(1.00784,1.00811),
    elements[2] => _pp("4.002602(2))"),
    elements[3] => NumberInterval(6.938,6.997),
    elements[4] =>_pp("9.0121831(5)"),
    elements[5] => NumberInterval(10.806,10.821),
    elements[6] => NumberInterval(12.0096,12.0116),
    elements[7] => NumberInterval(14.00643,14.00728),
    elements[8] => NumberInterval(15.99903,15.99977),
    elements[9] =>_pp("18.998403163(6)"),
    elements[10] =>_pp("20.1797(6)"),
    elements[11] =>_pp("22.98976928(2)"),
    elements[12] => NumberInterval(24.304,24.307),
    elements[13] =>_pp("26.9815384(3)"),
    elements[14] => NumberInterval(28.084,28.086),
    elements[15] =>_pp("30.973761998(5)"),
    elements[16] => NumberInterval(32.059,32.076),
    elements[17] => NumberInterval(35.446,35.457),
    elements[18] => NumberInterval(39.792,39.963),
    elements[19] =>_pp("39.0983(1)"),
    elements[20] =>_pp("40.078(4)"),
    elements[21] =>_pp("44.955908(5)"),
    elements[22] => _pp("47.867(1)"),
    elements[23] => _pp("50.9415(1)"),
    elements[24] => _pp("51.9961(6)"),
    elements[25] => _pp("54.938043(2)"),
    elements[26] => _pp("55.845(2)"),
    elements[27] => _pp("58.933194(3)"),
    elements[28] => _pp("58.6934(4)"),
    elements[29] => _pp("63.546(3)"),
    elements[30] => _pp("65.38(2)"),
    elements[31] => _pp("69.723(1)"),
    elements[32] => _pp("72.630(8)"),
    elements[33] => _pp("74.921595(6)"),
    elements[34] => _pp("78.971(8)"),
    elements[35] => NumberInterval(79.901,79.907),
    elements[36] => _pp("83.798(2)"),
    elements[37] => _pp("85.4678(3)"),
    elements[38] => _pp("87.62(1)"),
    elements[39] => _pp("88.90584(1)"),
    elements[40] => _pp("91.224(2)"),
    elements[41] => _pp("92.90637(1)"),
    elements[42] => _pp("95.95(1)"),
    elements[44] => _pp("101.07(2)"),
    elements[45] => _pp("102.90549(2)"),
    elements[46] => _pp("106.42(1)"),
    elements[47] => _pp("107.8682(2)"),
    elements[48] => _pp("112.414(4)"),
    elements[49] => _pp("114.818(1)"),
    elements[50] => _pp("118.710(7)"),
    elements[51] => _pp("121.760(1)"),
    elements[52] => _pp("127.60(3)"),
    elements[53] => _pp("126.90447(3)"),
    elements[54] => _pp("131.293(6)"),
    elements[55] => _pp("132.90545196(6)"),
    elements[56] => _pp("137.327(7)"),
    elements[57] => _pp("138.90547(7)"),
    elements[58] => _pp("140.116(1)"),
    elements[59] => _pp("140.90766(1)"),
    elements[60] => _pp("144.242(3)"),
    elements[62] => _pp("150.36(2)"),
    elements[63] => _pp("151.964(1)"),
    elements[64] => _pp("157.25(3)"),
    elements[65] => _pp("158.925354(8)"),
    elements[66] => _pp("162.500(1)"),
    elements[67] => _pp("164.930328(7)"),
    elements[68] => _pp("167.259(3)"),
    elements[69] => _pp("168.934218(6)"),
    elements[70] => _pp("173.045(10)"),
    elements[71] => _pp("174.9668(1)"),
    elements[72] => _pp("178.486(6)"),
    elements[73] => _pp("180.94788(2)"),
    elements[74] => _pp("183.84(1)"),
    elements[75] => _pp("186.207(1)"),
    elements[76] => _pp("190.23(3)"),
    elements[77] => _pp("192.217(2)"),
    elements[78] => _pp("195.084(9)"),
    elements[79] => _pp("196.966570(4)"),
    elements[80] => _pp("200.592(3)"),
    elements[81] => NumberInterval(204.382,204.385),
    elements[82] => NumberInterval(206.14,207.94),
    elements[83] => _pp("208.98040(1)"),
    elements[90] => _pp("232.0377(4)"),
    elements[91] => _pp("231.03588(1)"),
    elements[92] => _pp("238.02891(3)")
)

NeXLUncertainties.value(r::NumberInterval) = mid(r)

"""
    asa(::Type{DataFrame}, links::Dict{Element,String})

Create a DataFrame which contains a periodic table with links to URLs.
This doesn't work so well at the REPL when represented as text but works
nicely in HTML.
"""
function NeXLUncertainties.asa(::Type{DataFrame}, links::Dict{Element,String})
    df=DataFrame(
        IA=Markdown.MD[],
        IIA=Markdown.MD[],
        IIIB=Markdown.MD[],
        IVB=Markdown.MD[],
        VB=Markdown.MD[],
        VIB=Markdown.MD[],
        VIIB=Markdown.MD[],
        VIII_1=Markdown.MD[],
        VIII_2=Markdown.MD[],
        VIII_3=Markdown.MD[],
        IB=Markdown.MD[],
        IIB=Markdown.MD[],
        IIIA=Markdown.MD[],
        IVA=Markdown.MD[],
        VA=Markdown.MD[],
        VIA=Markdown.MD[],
        VIIA=Markdown.MD[],
        VIIIA=Markdown.MD[]
    )
    for row in 1:10
        r = fill(Markdown.parse(""), 18)
        for el in filter(e->e.ypos==row, collect(elements))
            r[el.xpos]= Markdown.parse(haskey(links,el) ? "[$(el.symbol)]($(links[el]))" : el.symbol)
        end
        push!(df, r)
    end
    return df
end


text:
using Test, XPORTA

@testset "chsh violation" begin

using BellScenario: QBase, QuantumBehavior

chsh_ineqs = XPORTA.read_ieq("./test/BellComm/integration/files/22-22_no_signaling.ieq").inequalities
bell_states = (
    1/sqrt(2)*[1;0;0;1],
    1/sqrt(2)*[1;0;0;-1],
    1/sqrt(2)*[0;1;1;0],
    1/sqrt(2)*[0;1;-1;0]
)

@testset "maximal violation with bell states" begin

    povm_a0 = [[1 0; 0 0], [0 0; 0 1]] # z measurement
    povm_a1 = [[0.5 0.5;0.5 0.5], [0.5 -0.5; -0.5 0.5]] # x measurement
    povm_b0 = [
        0.5*([1 0; 0 1] + 1/sqrt(2)*([1 0 ; 0 -1] + [0 1; 1 0])),
        0.5*([1 0; 0 1] - 1/sqrt(2)*([1 0 ; 0 -1] + [0 1; 1 0]))
    ] # -z - x measurement
    povm_b1 = [
        0.5*([1 0; 0 1] + 1/sqrt(2)*([1 0 ; 0 -1] - [0 1; 1 0])),
        0.5*([1 0; 0 1] - 1/sqrt(2)*([1 0 ; 0 -1] - [0 1; 1 0]))
    ] # z - x measuremnt

    @test QBase.is_povm(povm_a0)
    @test QBase.is_povm(povm_a1)
    @test QBase.is_povm(povm_b0)
    @test QBase.is_povm(povm_b1)

    for ψ in bell_states
        behavior = QuantumBehavior.bipartite_scenario([povm_a0,povm_a1],[povm_b0,povm_b1],ψ)

        violation_count = 0
        for ineq in eachrow(chsh_ineqs)

            exp_val = (cat(-1*ineq[end],ineq[1:end-1]',dims=2)*behavior)[1]
            if (exp_val > 0) & !(exp_val + 1 ≈ 1)
                violation_count += 1
            end
        end

        @test violation_count == 1
    end
end
end


text:
include("Solution.jl")
include("Params.jl")

""" Recombination Operations for the Genetic Algorithm

Implemented recombinations:
- Order 1 Crossover - order1cx
- Partially Mapped Crossover (PMX) - PMXcx
- Cycle Crossover - cyclecx

TODO:
- Edge Recombination
- Order Multiple Crossover
- Direct Insertion Crossover
- Cohesive Merging Procedure
- Scrambled Merging Procedure
"""

##############################################################################
# Recombination operations

function hasUniqueValues(y::Array{Int64})::Bool
    u = unique(y)
    d = Dict([(i,count(x->x==i,y)) for i in u])
    return all(x->x==true, [v for (k,v) in d] .== 1)
end

function order1op(p1::Array{Int64},p2::Array{Int64},s=nothing)
    n = length(p1)
    if ~((sort(p1) == 1:n) && (sort(p2) == 1:n))
        throw("Input parents must be a permutation arrays")
    end
    c_perm = zeros(Int, n)
    if s == nothing
        s = rand(1:n)
    end
    cp_idx = s:min(n,s+Int(floor(n/2))-1)
    if s+Int(floor(n/2))-1 > n
        cp_idx = vcat(cp_idx, 1:(s+Int(floor(n/2))-1) % n)
    end
    c_perm[cp_idx] = p1[cp_idx]

    # Child 1
    filled = 0
    i = s+Int(floor(n/2))
    if i > n
        i = (s+Int(floor(n/2))) % n
    end
    j = i
    while filled < Int(ceil(n/2))
        if i > n
            i = 1
        end
        if j > n
            j = 1
        end
        tmp = p2[i]
        if ~(tmp in c_perm)
            c_perm[j] = tmp
            filled += 1
            j += 1
        end
        i += 1
    end
    if ~hasUniqueValues(c_perm)
        throw("Child permutation must have unique values")
    end
    if sort(c_perm) != 1:n
        throw("Child must be a permutation")
    end
    return c_perm, s
end

function order1cx(p1::Solution,p2::Solution)
    # Order 1 crossover: http://www.rubicite.com/Tutorials/GeneticAlgorithms/CrossoverOperators/Order1CrossoverOperator.aspx
    c1_perm,s = order1op(p1.permutation,p2.permutation) # Child 1
    c2_perm,_ = order1op(p2.permutation,p1.permutation,s) # Child 2
    params = p1.params
    c1 = Solution(params,c1_perm)
    c2 = Solution(params,c2_perm)
    return (c1,c2)
end

function PMXop(p1::Array{Int64},p2::Array{Int64},s=nothing)
    n = length(p1)
    if ~((sort(p1) == 1:n) && (sort(p2) == 1:n))
        throw("Input parents must be permutation arrays")
    end
    c_perm = zeros(Int, n)
    if s == nothing
        s = rand(1:n)
    end
    cp_idx = s:min(n,s+Int(floor(n/2))-1)
    if s+Int(floor(n/2))-1 > n
        cp_idx = vcat(cp_idx, 1:(s+Int(floor(n/2))-1) % n)
    end
    c_perm[cp_idx] = p1[cp_idx]

    filled_idx = cp_idx
    for i in cp_idx
        p2_idx = i
        p2_val = p2[p2_idx]
        if ~(p2_val in c_perm[cp_idx])
            include_val = p2_val
            while (p2_idx in filled_idx)
                p1_idx = p2_idx
                p1_val = p1[p1_idx]
                p2_idx = findfirst(p2 .== p1_val)
                p2_val = p2[p2_idx]
            end
            c_perm[p2_idx] = include_val
        end
    end
    c_perm[c_perm .== 0] = p2[c_perm .== 0]

    if ~hasUniqueValues(c_perm)
        throw("Child permutation must have unique values")
    end
    if sort(c_perm) != 1:n
        throw("Child must be a permutation")
    end
    return c_perm, s
end

function PMXcx(p1::Solution,p2::Solution)
    # PMX crossover: http://www.rubicite.com/Tutorials/GeneticAlgorithms/CrossoverOperators/PMXCrossoverOperator.aspx
    c1_perm,s = PMXop(p1.permutation,p2.permutation) # Child 1
    c2_perm,_ = PMXop(p2.permutation,p1.permutation,s) # Child 2
    params = p1.params # TODO: p1.params and p2.params should be the same. Explicitly check if they are the same
    c1 = Solution(params,c1_perm)
    c2 = Solution(params,c2_perm)
    return (c1,c2)
end

function cycleop(p1::Array{Int64},p2::Array{Int64})
    n = length(p1)
    if ~((sort(p1) == 1:n) && (sort(p2) == 1:n))
        throw("Input parents must be permutation arrays")
    end
    c_perm = zeros(Int, n)
    allcycles = []
    visited_idx = []

    p2_idx = 1
    p2_val = p2[1]

    # Identify all cycles
    while length(visited_idx) < n
        cycle = []
        p2_idx = findfirst(x -> ~(x in visited_idx), 1:n)
        p2_val = p2[p2_idx]
        start_val = p2_val
        p1_val = 0
        while (p1_val != start_val)
            append!(cycle,p2_idx)
            p1_idx = p2_idx
            p1_val = p1[p1_idx]
            p2_idx = findfirst(p2 .== p1_val)
            p2_val = p2[p2_idx]
        end
        append!(visited_idx,cycle)
        if length(cycle) > 1
            push!(allcycles,cycle)
        end
    end

    # Copy alternate cycles to child
    for (i,cycle) in enumerate(allcycles)
        if i % 2 == 0
            c_perm[cycle] = p2[cycle]
        else
            c_perm[cycle] = p1[cycle]
        end
    end

    c_perm[c_perm .== 0] = p1[c_perm .== 0]

    if ~hasUniqueValues(c_perm)
        throw("Child permutation must have unique values")
    end
    if sort(c_perm) != 1:n
        throw("Child must be a permutation")
    end
    return c_perm
end

function cyclecx(p1::Solution,p2::Solution)
    # Cycle crossover: http://www.rubicite.com/Tutorials/GeneticAlgorithms/CrossoverOperators/CycleCrossoverOperator.aspx
    c1_perm = cycleop(p1.permutation,p2.permutation) # Child 1
    c2_perm = cycleop(p2.permutation,p1.permutation) # Child 2
    params = p1.params # TODO: p1.params and p2.params should be the same. Explicitly check if they are the same
    c1 = Solution(params,c1_perm)
    c2 = Solution(params,c2_perm)
    return (c1,c2)
end

function cohesiveMergecx(p1::Solution,p2::Solution)
    # Cohesive Merging crossover: Zvi Drezner, "A new Genetic Algorithm for the QAP"
    # TODO
    return (c1,c2)
end

function scrambledMergecx(p1::Solution,p2::Solution)
    # Scrambled Merging crossover: Zvi Drezner, "A new Genetic Algorithm for the QAP"
    # TODO
    return (c1,c2)
end

function recombinationOp(parents,recombine::Function=order1cx)::Array{Solution}
    offspring = Array{Solution}(undef, 0)
    for (p1,p2) in parents
        (c1,c2) = recombine(p1,p2)
        offspring = vcat(offspring, c1, c2)
    end
    return offspring
end

max_stars_repo_name: BSAraujo/QAP-Metaheuristics

max_issues_repo_name: BSAraujo/QAP-Metaheuristics

max_forks_repo_name: BSAraujo/QAP-Metaheuristics