text:
###############################################################################
#
#  Ray Class Group: Ctx
#
###############################################################################

mutable struct ctx_rayclassgrp
  order::NfOrd
  class_group_map::MapClassGrp #The class group mod n map
  n::Int #the n for n_quo
  diffC::fmpz #exponent of the full class group, divided by n
  units::Vector{FacElem{nf_elem, AnticNumberField}}
  princ_gens::Vector{FacElem{nf_elem, AnticNumberField}}

  computed::Vector{Tuple{Dict{NfOrdIdl, Int}, Bool, MapRayClassGrp}}
  multiplicative_groups::Dict{NfOrdIdl, GrpAbFinGenToAbsOrdQuoRingMultMap}

  function ctx_rayclassgrp()
    z = new()
    return z
  end
end

function rayclassgrp_ctx(O::NfOrd, expo::Int)

  C1, mC1 = class_group(O, use_aut = true)
  valclass = ppio(exponent(C1), fmpz(expo))[1]
  C, mC = Hecke.n_part_class_group(mC1, expo)
  ctx = ctx_rayclassgrp()
  ctx.order = O
  ctx.class_group_map = mC
  ctx.n = expo
  ctx.diffC = Int(divexact(exponent(C1), exponent(C)))
  ctx.multiplicative_groups = Dict{NfOrdIdl, GrpAbFinGenToAbsOrdQuoRingMultMap}()
  return ctx

end

function assure_elements_to_be_eval(ctx::ctx_rayclassgrp)
  if isdefined(ctx, :units)
    return nothing
  end
  OK = ctx.order
  U, mU = unit_group_fac_elem(OK)
  mC = ctx.class_group_map
  _assure_princ_gen(mC)
  n = ctx.n
  units = FacElem{nf_elem, AnticNumberField}[_preproc(mU(U[i]), fmpz(n)) for i = 1:ngens(U)]
  princ_gens = FacElem{nf_elem, AnticNumberField}[_preproc(mC.princ_gens[i][2], fmpz(n)) for i = 1:length(mC.princ_gens)]
  ctx.units = units
  ctx.princ_gens = princ_gens
  return nothing
end

###############################################################################
#
#  Ray Class Group: Fields function
#
###############################################################################

function empty_ray_class(m::NfOrdIdl)
  O = order(parent(m))
  X = abelian_group(Int[])

  local exp
  let O = O
    function exp(a::GrpAbFinGenElem)
      return FacElem(Dict(ideal(O,1) => fmpz(1)))
    end
  end

  local disclog
  let X = X
    function disclog(J::Union{NfOrdIdl, FacElem{NfOrdIdl}})
      return id(X)
    end
  end

  mp = Hecke.MapRayClassGrp()
  mp.header = Hecke.MapHeader(X, FacElemMon(parent(m)) , exp, disclog)
  mp.defining_modulus = (m, InfPlc[])
  return X,mp

end

function ray_class_group_quo(m::NfOrdIdl, y1::Dict{NfOrdIdl,Int}, y2::Dict{NfOrdIdl,Int}, inf_plc::Vector{InfPlc}, ctx::ctx_rayclassgrp; check::Bool = true, GRH::Bool = true)

  mC = ctx.class_group_map
  C = domain(mC)
  diffC = ctx.diffC

  if isempty(y1) && isempty(y2) && isempty(inf_plc)
    local exp_c
    let mC = mC
      function exp_c(a::GrpAbFinGenElem)
        return FacElem(Dict(mC(a) => 1))
      end
    end
    return class_as_ray_class(C, mC, exp_c, m, ctx.n)
  end

  O = ctx.order
  K = nf(O)
  n_quo = ctx.n

  lp = merge(max, y1, y2)

  powers = Vector{Tuple{NfOrdIdl, NfOrdIdl}}()
  quo_rings = Tuple{NfOrdQuoRing, Hecke.AbsOrdQuoMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[]
  groups_and_maps = Tuple{GrpAbFinGen, Hecke.GrpAbFinGenToAbsOrdQuoRingMultMap{NfAbsOrd{AnticNumberField,nf_elem},NfAbsOrdIdl{AnticNumberField,nf_elem},NfAbsOrdElem{AnticNumberField,nf_elem}}}[]
  for (PP, ee) in lp
    if isone(ee)
      dtame = Dict{NfOrdIdl, Int}(PP => 1)
      dwild = Dict{NfOrdIdl, Int}()
    else
      dwild = Dict{NfOrdIdl, Int}(PP => ee)
      if haskey(y1, PP)
        dtame = Dict{NfOrdIdl, Int}(PP => 1)
      else
        dtame = Dict{NfOrdIdl, Int}()
      end
    end
    QQ = PP^ee

    push!(powers, (PP, QQ))
    if haskey(ctx.multiplicative_groups, QQ)
      cached_map = ctx.multiplicative_groups[QQ]
      QQQ = codomain(cached_map)
      mQQQ = AbsOrdQuoMap(O, QQQ)
      push!(quo_rings, (QQQ, mQQQ))
      new_map = GrpAbFinGenToAbsOrdQuoRingMultMap(domain(cached_map), QQQ, copy(cached_map.generators), cached_map.discrete_logarithm)
      new_map.tame = copy(cached_map.tame)
      new_map.wild = copy(cached_map.wild)
      push!(groups_and_maps, (domain(new_map), new_map))
    else
      QQQ, mQQQ = quo(O, QQ)
      if ee == 1
        QQQ.factor = dtame
      else
        QQQ.factor = dwild
      end
      push!(quo_rings, (QQQ, mQQQ))
      gandm = _mult_grp_mod_n(quo_rings[end][1], dtame, dwild, n_quo)
      map_to_cache = GrpAbFinGenToAbsOrdQuoRingMultMap(gandm[1], QQQ, copy(gandm[2].generators), gandm[2].discrete_logarithm)
      map_to_cache.tame = copy(gandm[2].tame)
      map_to_cache.wild = copy(gandm[2].wild)
      ctx.multiplicative_groups[QQ] = map_to_cache
      push!(groups_and_maps, gandm)
    end
  end

  if isempty(groups_and_maps)
    nG = 0
    expon = fmpz(1)
  else
    nG = sum(ngens(x[1]) for x in groups_and_maps)
    expon = lcm([exponent(x[1]) for x in groups_and_maps])
  end

  if iseven(n_quo)
    p = InfPlc[ x for x in inf_plc if isreal(x) ]
  else
    p = InfPlc[]
  end
  H, lH = log_infinite_primes(O, p)
  expon = lcm(expon, exponent(H))

  if exponent(C)*expon < n_quo && check
    return empty_ray_class(m)
  end

  assure_elements_to_be_eval(ctx)
  exp_class, Kel = find_coprime_representatives(mC, m, lp)

  nU = length(ctx.units)
  # We construct the relation matrix and evaluate units and relations with the class group in the quotient by m
  # Then we compute the discrete logarithms

  R = zero_matrix(FlintZZ, 2*(ngens(C)+nG+ngens(H))+nU, ngens(C)+ngens(H)+nG)
  for i = 1:ncols(R)
    R[i+ngens(C)+nG+ngens(H)+nU, i] = n_quo
  end
  for i=1:ngens(C)
    R[i, i] = C.snf[i]
  end
  ind = 1
  for s = 1:length(quo_rings)
    G = groups_and_maps[s][1]
    @assert issnf(G)
    for i = 1:ngens(G)
      R[i+ngens(C)+ind-1, i+ngens(C)+ind-1] = G.snf[i]
    end
    ind += ngens(G)
  end
  for i = 1:ngens(H)
    R[ngens(C)+nG+i, ngens(C)+nG+i] = 2
  end

  @vprint :RayFacElem 1 "Collecting elements to be evaluated; first, units \n"
  tobeeval = Vector{FacElem{nf_elem, AnticNumberField}}(undef, nU+ngens(C))
  for i = 1:nU
    tobeeval[i] = ctx.units[i]
  end
  for i = 1:ngens(C)
    tobeeval[i+nU] = _preproc(ctx.princ_gens[i]*FacElem(Dict(Kel[i] => C.snf[i])), expon)
  end

  ind = 1
  for i = 1:length(groups_and_maps)
    exp_q = gcd(expon, norm(powers[i][2])- divexact(norm(powers[i][2]), norm(powers[i][1])))
    evals = fac_elems_eval(powers[i][1], powers[i][2], tobeeval, exp_q)
    Qr = quo_rings[i][1]
    mG = groups_and_maps[i][2]
    for j = 1:nU
      a = (mG\Qr(evals[j])).coeff
      for s = 1:ncols(a)
        R[j+nG+ngens(H)+ngens(C), ngens(C)+s+ind-1] = a[1, s]
      end
    end

    for j = 1:ngens(C)
      a = (mG\Qr(evals[j+nU])).coeff
      for s = 1:ncols(a)
        R[j, ngens(C)+ind+s-1] = -a[1, s]
      end
    end
    ind += ngens(groups_and_maps[i][1])
  end

  if !isempty(p)
    for j = 1:nU
      aa = lH(tobeeval[j])::GrpAbFinGenElem
      for s = 1:ngens(H)
        R[j+nG+ngens(C)+ngens(H), ngens(C)+ind-1+s] = aa[s]
      end
    end
    for j = 1:ngens(C)
      aa = lH(tobeeval[j+nU])::GrpAbFinGenElem
      for s = 1:ngens(H)
        R[j, ngens(C)+ind-1+s] = -aa[s]
      end
    end
  end

  X = abelian_group(R)
  X.exponent = n_quo
  if isone(order(X))
    mp = MapRayClassGrp()
    mp.header = MapHeader(X, FacElemMon(parent(m)))
    return X, mp
  end

  invd = invmod(fmpz(diffC), expon)
  local disclog
  let X = X, mC = mC, invd = invd, C = C, exp_class = exp_class, powers = powers, groups_and_maps = groups_and_maps, quo_rings = quo_rings, lH = lH, diffC = diffC, n_quo = n_quo, m = m, p = p, expon = expon

    # Discrete logarithm
    function disclog(J::FacElem{NfOrdIdl, NfOrdIdlSet})
      @vprint :RayFacElem 1 "Disc log of element $J \n"
      a1 = id(X)
      for (f, k) in J
        a1 += k*disclog(f)
      end
      return a1
    end

    function disclog(J::NfOrdIdl)
      @hassert :RayFacElem 1 iscoprime(J, m)
      if isone(J)
        @vprint :RayFacElem 1 "J is one \n"
        return id(X)
      end
      coeffs = zero_matrix(FlintZZ, 1, ngens(X))
      if J.is_principal == 1 && isdefined(J, :princ_gen)
        z = FacElem(Dict(J.princ_gen.elem_in_nf => diffC))
      else
        L = mC\J
        for i = 1:ngens(C)
          coeffs[1, i] = L[i]
        end
        @vprint :RayFacElem 1 "Disc log of element J in the Class Group: $(L.coeff) \n"
        eL = exp_class(L)
        inv!(eL)
        add_to_key!(eL.fac, J, 1)
        pow!(eL, diffC)
        @vprint :RayFacElem 1 "This ideal is principal: $eL \n"
        z = principal_generator_fac_elem(eL)
      end
      ii = 1
      z1 = _preproc(z, expon)
      for i = 1:length(powers)
        P, Q = powers[i]
        exponq = gcd(expon, norm(Q)-divexact(norm(Q), norm(P)))
        el = fac_elems_eval(P, Q, FacElem{nf_elem, AnticNumberField}[z1], exponq)
        y = (invd*(groups_and_maps[i][2]\quo_rings[i][1](el[1]))).coeff
        for s = 1:ncols(y)
          coeffs[1, ii-1+ngens(C)+s] = y[1, s]
        end
        ii += ngens(groups_and_maps[i][1])
      end
      if !isempty(p)
        b = lH(z).coeff
        for s = 1:ncols(b)
          coeffs[1, ii-1+s+ngens(C)] = b[1, s]
        end
      end
      return GrpAbFinGenElem(X, coeffs)
    end
  end

  Dgens = Tuple{NfOrdElem, GrpAbFinGenElem}[]
  ind = 1
  #We need generators of the full multiplicative group
  #In particular, we need the idempotents...
  if !isone(length(powers))
    for i = 1:length(powers)
      P, QQ = powers[i]
      mG = groups_and_maps[i][2]
      J = ideal(O, 1)
      minJ = fmpz(1)
      mins = fmpz(1)
      for j = 1:length(powers)
        if minimum(powers[j][1], copy = false) != minimum(P, copy = false)
          mins = lcm(mins, minimum(powers[j][2], copy = false))
          continue
        end
        if j != i
          J *= powers[j][2]
          minJ = lcm(minJ, minimum(powers[j][2], copy = false))
        end
      end
      J.minimum = minJ
      i1, i2 = idempotents(QQ, J)
      if !isone(mins)
        d, u1, v1 = gcdx(minimum(QQ, copy = false), mins)
        i1 = i1*(u1*minimum(QQ, copy = false) + v1*mins) + u1*minimum(QQ, copy = false) *i2
        i2 = v1*mins*i2
      end
      @hassert :RayFacElem 1 isone(i1+i2)
      @hassert :RayFacElem 1 i1 in QQ
      @hassert :RayFacElem 1 i2 in J
      @hassert :RayFacElem 1 i2 in ideal(O, mins)
      gens = mG.generators
      if isempty(p)
        if haskey(mG.tame, P)
          gens_tame = mG.tame[P].generators
          for s = 1:length(gens_tame)
            gens_tame[s] = gens_tame[s]*i2 + i1
            @hassert :RayFacElem 1 isone(mod(gens_tame[s], J*ideal(O, mins)))
          end
          mG.tame[P].generators = gens_tame
        end
        if haskey(mG.wild, P)
          gens_wild = mG.wild[P].generators
          for s = 1:length(gens_wild)
            gens_wild[s] = gens_wild[s]*i2 + i1
            @hassert :RayFacElem 1 isone(mod(gens_wild[s], J*ideal(O, mins)))
          end
          mG.wild[P].generators = gens_wild
        end
        for s = 1:length(gens)
          el_to_push = gens[s].elem*i2+i1
          @hassert :RayFacElem 1 isone(mod(el_to_push, J*ideal(O, mins)))
          push!(Dgens, (el_to_push, X[ngens(C)+ind-1+s]))
        end
      else
        if haskey(mG.tame, P)
          gens_tame = mG.tame[P].generators
          for s = 1:length(gens_tame)
            gens_tame[s] = make_positive(gens_tame[s]*i2 + i1, minimum(m, copy = false))
          end
          mG.tame[P].generators = gens_tame
        end
        if haskey(mG.wild, P)
          gens_wild = mG.wild[P].generators
          for s = 1:length(gens_wild)
            gens_wild[s] = make_positive(gens_wild[s]*i2 + i1, minimum(m, copy = false))
          end
          mG.wild[P].generators = gens_wild
        end
        for s = 1:length(gens)
          elgen = make_positive(gens[s].elem*i2 + i1, minimum(m, copy = false))
          push!(Dgens, (elgen, X[ngens(C)+ind-1+s]))
        end
      end
      ind += length(gens)
    end
  else
    mG = groups_and_maps[1][2]
    gens = mG.generators
    if !isempty(p)
      for s = 1:length(gens)
        elgen = make_positive(gens[s].elem, minimum(m, copy = false))
        push!(Dgens, (elgen, X[ngens(C)+s]))
      end
    else
      for s = 1:length(gens)
        elgen = gens[s].elem
        push!(Dgens, (elgen, X[ngens(C)+s]))
      end
    end
  end

  ind = 1
  for i = 1:length(powers)
    mG = groups_and_maps[i][2]
    for (prim, mprim) in mG.tame
      new_mprim = GrpAbFinGenToAbsOrdMap(domain(mprim), codomain(mprim), copy(mprim.generators), mprim.discrete_logarithm)
      dis = zero_matrix(FlintZZ, 1, ngens(X))
      to_be_c = mprim.disc_log.coeff
      for i = 1:length(to_be_c)
        dis[1, ind-1+i+ngens(C)] = to_be_c[1, i]
      end
      new_mprim.disc_log = GrpAbFinGenElem(X, dis)
      mG.tame[prim] = new_mprim
    end
    ind += ngens(domain(mG))
  end

  disc_log_inf = Dict{InfPlc, GrpAbFinGenElem}()
  for i = 1:length(p)
    eldi = zero_matrix(FlintZZ, 1, ngens(X))
    eldi[1, ngens(X) - length(inf_plc) + i] = 1
    disc_log_inf[p[i]] = GrpAbFinGenElem(X, eldi)
  end

  mp = MapRayClassGrp()
  mp.header = MapHeader(X, FacElemMon(parent(m)))
  mp.header.preimage = disclog
  mp.fact_mod = lp
  mp.defining_modulus = (m, p)
  mp.powers = powers
  mp.groups_and_maps = groups_and_maps
  mp.disc_log_inf_plc = disc_log_inf
  mp.gens_mult_grp_disc_log = Dgens
  mp.clgrpmap = mC
  X.exponent = fmpz(n_quo)
  return X, mp

end

function log_infinite_primes(O::NfOrd, p::Array{InfPlc,1})
  if isempty(p)
    S = abelian_group(Int[])

    local log1
    let S = S
      function log1(B::T) where T <: Union{nf_elem ,FacElem{nf_elem, AnticNumberField}}
        return id(S)
      end
    end
    return S, log1
  end

  S = abelian_group(Int[2 for i=1:length(p)])
  local log
  let S = S, p = p
    function log(B::T) where T <: Union{nf_elem ,FacElem{nf_elem, AnticNumberField}}
      emb = signs(B, p)
      ar = zero_matrix(FlintZZ, 1, length(p))
      for i = 1:length(p)
        if emb[p[i]] == -1
          ar[1, i] = 1
        end
      end
      return GrpAbFinGenElem(S, ar)
    end
  end
  return S, log
end

function ray_class_group_quo(O::NfOrd, m::Int, wprimes::Dict{NfOrdIdl,Int}, inf_plc::Array{InfPlc,1}, ctx::ctx_rayclassgrp; GRH::Bool = true)

  d1 = Dict{NfOrdIdl, Int}()
  lp = factor(m)
  I = ideal(O, 1)
  minI = fmpz(1)
  for q in keys(lp.fac)
    lq = prime_decomposition(O, q)
    for (P, e) in lq
      if !haskey(wprimes, P)
        I *= P
        d1[P] = 1
      else
        d1[P] = 1
      end
    end
    minI = q*minI
  end
  if !isempty(wprimes)
    for (p, v) in wprimes
      qq = p^v
      I *= qq
      minI = lcm(minI, minimum(qq, copy = false))
    end
  end
  I.minimum = minI
  return ray_class_group_quo(I, d1, wprimes, inf_plc, ctx, GRH = GRH)

end

function ray_class_group_quo(O::NfOrd, y::Dict{NfOrdIdl, Int}, inf_plc::Array{InfPlc, 1}, ctx::ctx_rayclassgrp; GRH::Bool = true, check::Bool = true)

  y1=Dict{NfOrdIdl,Int}()
  y2=Dict{NfOrdIdl,Int}()
  n = ctx.n
  for (q, e) in y
    if gcd(norm(q)-1, n) != 1
      y1[q] = 1
      if gcd(norm(q), n) != 1 && e >= 2
        y2[q] = e
      end
    elseif gcd(norm(q), n) != 1 && e >= 2
      y2[q] = e
    end
  end
  I = ideal(O,1)
  for (q, vq) in y1
    I *= q
  end
  for (q, vq) in y2
    I *= q^vq
  end
  return ray_class_group_quo(I, y1, y2, inf_plc, ctx, GRH = GRH, check = check)

end


text:
"""
    parse_file(file; flex_load=true, <keyword arguments>)

Parse a Matpower .m `file` or PTI (PSS(R)E-v33) .raw `file` into a FlexPlan data structure,
including non-dispatchable generators, DC components, storage and flexible loads.

`flex_load` specifies whether to process flexible load data.
Other keyword arguments, if any, are forwarded to `PowerModels.parse_file`.

Mandatory tables: `bus`, `gen`, `branch` (and `load_extra` if `flex_load==true`).

Optional tables: `gencost`, `ndgen`, `branch_oltc`, `storage`, `storage_extra`,
`ne_storage`, and tables used by PowerModelsACDC.
Other tables can be added as well: they will be made available in the returned object.
"""
function parse_file(file::String; flex_load=true, kwargs...)
    data = _PM.parse_file(file; kwargs...)
    add_gen_data!(data)
    if !haskey(data, "ne_branch")
        data["ne_branch"] = Dict{String,Any}()
    end
    if haskey(data, "busdc") || haskey(data, "busdc_ne")
        _PMACDC.process_additional_data!(data)
    end
    add_storage_data!(data)
    if flex_load
        if !haskey(data, "load_extra")
            Memento.error(_LOGGER, "No `load_extra` table found in input file.")
        end
        add_flexible_demand_data!(data)
    end
    return data
end

"Add a `dispatchable` bool field to all generators; add non-dispatchable generators to `data[\"gen\"]`."
function add_gen_data!(data::Dict{String,Any})
    for dgen in values(data["gen"])
        dgen["dispatchable"] = true
    end

    if haskey(data, "ndgen")
        offset = length(data["gen"])
        rescale      = x -> x/data["baseMVA"]
        rescale_dual = x -> x*data["baseMVA"]
        for ndgen in values(data["ndgen"])
            ndgen["dispatchable"] = false

            # Convert to p.u.
            _PM._apply_func!(ndgen, "pref", rescale)
            _PM._apply_func!(ndgen, "qmax", rescale)
            _PM._apply_func!(ndgen, "qmin", rescale)
            _PM._apply_func!(ndgen, "cost_gen", rescale_dual)
            _PM._apply_func!(ndgen, "cost_curt", rescale_dual)

            # Define active power bounds using the same names used by dispatchable
            # generators.
            ndgen["pmin"] = 0.0
            ndgen["pmax"] = ndgen["pref"]
            delete!(ndgen, "pref")

            # Convert the cost of power produced by non-dispatchable generators into
            # polynomial form (the same used by dispatchable generators).
            ndgen["cost"] = [ndgen["cost_gen"], 0.0]
            delete!(ndgen, "cost_gen")

            # Assign to non-dispatchable generators ids contiguous to dispatchable
            # generators so that each generator has an unique id.
            new_id = ndgen["index"] += offset
            data["gen"]["$new_id"] = ndgen
        end
        delete!(data, "ndgen")
    end

    return data
end

function add_storage_data!(data)
    if haskey(data, "storage")
        for (s, storage) in data["storage"]
            rescale_power = x -> x/data["baseMVA"]
            _PM._apply_func!(storage, "max_energy_absorption", rescale_power)
            _PM._apply_func!(storage, "stationary_energy_outflow", rescale_power)
            _PM._apply_func!(storage, "stationary_energy_inflow", rescale_power)
        end
    else
        data["storage"] = Dict{String,Any}()
    end

    if haskey(data, "ne_storage")
        for (s, storage) in data["ne_storage"]
            rescale_power = x -> x/data["baseMVA"]
            _PM._apply_func!(storage, "energy_rating", rescale_power)
            _PM._apply_func!(storage, "thermal_rating", rescale_power)
            _PM._apply_func!(storage, "discharge_rating", rescale_power)
            _PM._apply_func!(storage, "charge_rating", rescale_power)
            _PM._apply_func!(storage, "energy", rescale_power)
            _PM._apply_func!(storage, "ps", rescale_power)
            _PM._apply_func!(storage, "qs", rescale_power)
            _PM._apply_func!(storage, "q_loss", rescale_power)
            _PM._apply_func!(storage, "p_loss", rescale_power)
            _PM._apply_func!(storage, "qmax", rescale_power)
            _PM._apply_func!(storage, "qmin", rescale_power)
            _PM._apply_func!(storage, "max_energy_absorption", rescale_power)
            _PM._apply_func!(storage, "stationary_energy_outflow", rescale_power)
            _PM._apply_func!(storage, "stationary_energy_inflow", rescale_power)
        end
    else
        data["ne_storage"] = Dict{String,Any}()
    end

    return data
end

function add_flexible_demand_data!(data)
    for (le, load_extra) in data["load_extra"]

        # ID of load point
        idx = load_extra["load_id"]

        # Superior bound on voluntary load reduction (not consumed power) as a fraction of the total reference demand (0 ≤ pred_rel_max ≤ 1)
        data["load"]["$idx"]["pred_rel_max"] = load_extra["pred_rel_max"]

        # Superior bound on upward demand shifted as a fraction of the total reference demand (0 ≤ pshift_up_rel_max ≤ 1)
        data["load"]["$idx"]["pshift_up_rel_max"] = load_extra["pshift_up_rel_max"]

        # Superior bound on downward demand shifted as a fraction of the total reference demand (0 ≤ pshift_down_rel_max ≤ 1)
        data["load"]["$idx"]["pshift_down_rel_max"] = load_extra["pshift_down_rel_max"]

        # Superior bound on shifted energy as a fraction of the total reference demand (0 ≤ eshift_rel_max ≤ 1)
        if haskey(load_extra, "eshift_rel_max")
            data["load"]["$idx"]["eshift_rel_max"] = load_extra["eshift_rel_max"]
        end

        # Compensation for consuming less (i.e. voluntary demand reduction) (€/MWh)
        data["load"]["$idx"]["cost_red"] = load_extra["cost_red"]

        # Recovery period for upward demand shifting (h)
        if haskey(load_extra, "tshift_up")
            data["load"]["$idx"]["tshift_up"] = load_extra["tshift_up"]
        end

        # Recovery period for downward demand shifting (h)
        if haskey(load_extra, "tshift_down")
            data["load"]["$idx"]["tshift_down"] = load_extra["tshift_down"]
        end

        # Compensation for demand shifting (€/MWh), applied only to the power shifted upward to avoid double counting
        data["load"]["$idx"]["cost_shift"] = load_extra["cost_shift"]

        # Compensation for load curtailment (i.e. involuntary demand reduction) (€/MWh)
        data["load"]["$idx"]["cost_curt"] = load_extra["cost_curt"]

        # Investment costs for enabling flexible demand (€)
        data["load"]["$idx"]["cost_inv"] = load_extra["cost_inv"]

        # Whether load is flexible (boolean)
        data["load"]["$idx"]["flex"] = load_extra["flex"]

        # Superior bound on voluntary energy reduction as a fraction of the total reference demand (0 ≤ ered_rel_max ≤ 1)
        if haskey(load_extra, "ered_rel_max")
            data["load"]["$idx"]["ered_rel_max"] = load_extra["ered_rel_max"]
        end

        # Expected lifetime of flexibility-enabling equipment (years)
        data["load"]["$idx"]["lifetime"] = load_extra["lifetime"]

        # CO2 costs for enabling flexible demand (€)
        if haskey(load_extra, "co2_cost")
            data["load"]["$idx"]["co2_cost"] = load_extra["co2_cost"]
        end

        # Power factor angle θ, giving the reactive power as Q = P ⨉ tan(θ)
        if haskey(load_extra, "pf_angle")
            data["load"]["$idx"]["pf_angle"] = load_extra["pf_angle"]
        end

        # Rescale cost and power input values to the p.u. values used internally in the model
        rescale_cost = x -> x*data["baseMVA"]
        rescale_power = x -> x/data["baseMVA"]
        _PM._apply_func!(data["load"]["$idx"], "cost_red", rescale_cost)
        _PM._apply_func!(data["load"]["$idx"], "cost_shift", rescale_cost)
        _PM._apply_func!(data["load"]["$idx"], "cost_curt", rescale_cost)
    end
    delete!(data, "load_extra")
    return data
end

function add_generation_emission_data!(data)
    rescale_emission = x -> x * data["baseMVA"]
    for (g, gen) in data["gen"]
        _PM._apply_func!(gen, "emission_factor", rescale_emission)
    end
    return data
end

ext: jl


text:
"constraint: `c^2 + d^2 <= a*b`"
function relaxation_complex_product_conic(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef)
    a_lb, a_ub = _IM.variable_domain(a)
    b_lb, b_ub = _IM.variable_domain(b)

    @assert (a_lb >= 0 && b_lb >= 0) || (a_ub <= 0 && b_ub <= 0)

    JuMP.@constraint(m, [a/sqrt(2), b/sqrt(2), c, 0] in JuMP.RotatedSecondOrderCone())
end

function relaxation_complex_product_conic_on_off(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef, d::JuMP.VariableRef, z::JuMP.VariableRef)
    a_lb, a_ub = _IM.variable_domain(a)
    b_lb, b_ub = _IM.variable_domain(b)
    c_lb, c_ub = _IM.variable_domain(c)
    d_lb, d_ub = _IM.variable_domain(d)
    z_lb, z_ub = _IM.variable_domain(z)

    @assert c_lb <= 0 && c_ub >= 0
    @assert d_lb <= 0 && d_ub >= 0

    JuMP.@constraint(m, [a/sqrt(2)*z_ub, b/sqrt(2), c, d] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(m, [a_ub/sqrt(2)*z, b/sqrt(2), c, d] in JuMP.RotatedSecondOrderCone())
    JuMP.@constraint(m, [a/sqrt(2), b_ub/sqrt(2)*z, c, d] in JuMP.RotatedSecondOrderCone())
end

function relaxation_complex_product_on_off(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef, d::JuMP.VariableRef) #to be moved to _IM
    a_lb, a_ub = _IM.variable_domain(a)
    b_lb, b_ub = _IM.variable_domain(b)
    c_lb, c_ub = _IM.variable_domain(c)

    @assert c_lb <= 0 && c_ub >= 0

    JuMP.@constraint(m, c^2 <= a*b*z_ub)
    JuMP.@constraint(m, c^2 <= a_ub*b*z)
    JuMP.@constraint(m, c^2 <= a*b_ub*z)
end

"constraint: `c^2 + d^2 <= a*b`"
function relaxation_complex_product(m::JuMP.Model, a::JuMP.VariableRef, b::JuMP.VariableRef, c::JuMP.VariableRef)
    a_lb, a_ub = _IM.variable_domain(a)
    b_lb, b_ub = _IM.variable_domain(b)

    @assert (a_lb >= 0 && b_lb >= 0) || (a_ub <= 0 && b_ub <= 0)

    JuMP.@constraint(m, c^2 <= a*b)
end

function relaxation_semicont_variable_on_off(m::JuMP.Model, a::JuMP.VariableRef, z::JuMP.VariableRef)
    a_lb, a_ub = _IM.variable_domain(a)

    JuMP.@constraint(m, a <= a_ub*z)
    JuMP.@constraint(m, a >= a_lb*z)
end

function relaxation_variable_on_off(m::JuMP.Model, x::JuMP.VariableRef, y::JuMP.VariableRef, z::JuMP.VariableRef)
    x_lb, x_ub = _IM.variable_domain(x)

    JuMP.@constraint(m, y <= x_ub*z)
    JuMP.@constraint(m, y >= x_lb*z)
end


text:
function areaplot(;x::AbstractVector = Int[],
                   y::AbstractVector = Int[],
                   group::AbstractVector = Int[],
                   stacked::Bool = false,
                   normalize::Bool = false)

    if stacked
    	v = barplot(x = x, y = y, group = group, stacked = true)

        #Get unique table name
        table = v.data[1].name

	    v.scales[1].points = true
		  v.scales[1]._type = nothing

        innermark = VegaMark(_type="area", properties=VegaMarkProperties(enter = VegaMarkPropertySet(x = VegaValueRef(scale = "x", field = "x"),
                                                                                y = VegaValueRef(scale = "y", field = "y"))))

        innermark.properties.enter.y = VegaValueRef(scale = "y", field = "layout_start")
        innermark.properties.enter.y2 = VegaValueRef(scale = "y", field = "layout_end")
        innermark.properties.enter.fill = VegaValueRef(scale = "group", field = "group")
        innermark.properties.enter.interpolate = VegaValueRef(value="monotone")

        mark = VegaMark(_type = "group",
                        from = VegaMarkFrom(data=table,
                                            transform=[VegaTransform(Dict{Any,Any}("type"=>"stack", "groupby" => ["x"], "sortby" => ["group"], "field" => "y", "offset" => normalize == true? "normalize" : "zero")),
                                                       VegaTransform(Dict{Any,Any}("type"=>"facet", "groupby" => ["group"]))]),
                        marks = [innermark])

        v.marks = [mark]

      if normalize
        ylim!(v, min = 0, max = 1)
      end

	else
        v = barplot(x = x, y = y, group = group)

        v.marks[1]._type = "area"
        v.marks[1].properties.enter.width = nothing
        v.marks[1].properties.enter.interpolate = VegaValueRef(value = "monotone")
        v.scales[1]._type = "linear"
        v.scales[1].zero = false
    end

    return v
end

max_stars_repo_name: jpthor/vega


text:
const Match = Tuple{Rule,Union{Nothing,Pattern},Sub,Int64}
const MatchesBuf = Vector{Match}

function cached_ids(g::EGraph, p::Pattern)# ::Vector{Int64}
    keys(g.classes)
end

# FIXME
function cached_ids(g::EGraph, p::PatTerm)
    # println("pattern $p, $(p.head)")
    # println("all ids")
    # keys(g.classes) |> println
    # println("cached symbols")
    # cached = get(g.symcache, p.head, Set{Int64}())
    # println("symbols where $(p.head) appears")
    # appears = Set{Int64}()
    # for (id, class) ∈ g.classes
    #     for n ∈ class
    #         if n.head == p.head
    #             push!(appears, id)
    #         end
    #     end
    # end
    # # println(appears)
    # if !(cached == appears)
    #     @show cached
    #     @show appears
    # end

    keys(g.classes)
    # cached
    # get(g.symcache, p.head, [])
end

# function cached_ids(g::EGraph, p::PatLiteral)
#     get(g.symcache, p.val, [])
# end

function search_rule!(g::EGraph, r::SymbolicRule, id, buf)
    append!(buf, [(r, r.right, sub, id) for sub in ematch(g, r.left, id)])
end

function search_rule!(g::EGraph, r::DynamicRule, id, buf)
    append!(buf, [(r, nothing, sub, id) for sub in ematch(g, r.left, id)])
end

function search_rule!(g::EGraph, r::BidirRule, id, buf)
    append!(buf, [(r, r.right, sub, id) for sub in ematch(g, r.left, id)])
    append!(buf, [(r, r.left, sub, id) for sub in ematch(g, r.right, id)])
end

function eqsat_search_threaded!(egraph::EGraph, theory::Vector{<:Rule},
        scheduler::AbstractScheduler)::MatchesBuf
    matches = MatchesBuf()
    mlock = ReentrantLock()

    inequalities = filter(theory) do rule
        rule isa UnequalRule
    end
    # never skip contradiction checks
    Threads.@threads for rule ∈ inequalities
        ids = cached_ids(egraph, rule.left)

        for id ∈ ids
            lock(mlock) do
                search_rule!(egraph, rule, id, matches)
            end
        end
    end

    other_rules = filter(theory) do rule
        !(rule isa UnequalRule)
    end
    Threads.@threads for rule ∈ other_rules
        # don't apply banned rules
        if !cansearch(scheduler, rule)
            # println("skipping banned rule $rule")
            continue
        end

        rule_matches = []
        ids = cached_ids(egraph, rule.left)

        for id ∈ ids
            search_rule!(egraph, rule, id, rule_matches)
        end
        can_yield = inform!(scheduler, rule, rule_matches)
        if can_yield
            lock(mlock) do
                append!(matches, rule_matches)
            end
        end
    end
    return matches
end

function eqsat_search!(egraph::EGraph, theory::Vector{<:Rule},
    scheduler::AbstractScheduler)::MatchesBuf
    matches = MatchesBuf()

    inequalities = filter(theory) do rule
        rule isa UnequalRule
    end
# never skip contradiction checks
    for rule ∈ inequalities
        ids = cached_ids(egraph, rule.left)

        for id ∈ ids
            search_rule!(egraph, rule, id, matches)
        end
    end

    other_rules = filter(theory) do rule
        !(rule isa UnequalRule)
    end
    for rule ∈ other_rules
    # don't apply banned rules
        if !cansearch(scheduler, rule)
        # println("skipping banned rule $rule")
            continue
        end

        rule_matches = []

        ids = cached_ids(egraph, rule.left)

        for id ∈ ids
            search_rule!(egraph, rule, id, rule_matches)
        end

        can_yield = inform!(scheduler, rule, rule_matches)
        if can_yield
            append!(matches, rule_matches)
        end
    end
    return matches
end


text:
using Pkg
using NumericalMethodsforEngineers
using Test

a = [10. 5. 6.; 5. 20. 4.; 6. 4. 30.]
x0 = [1., 1., 1.]

tol = 1.0e-5
limit = 100

(iters, eigval, eigvec) = nmex(a, x0; tol=tol, limit=limit, show=8)

println("\nNo of iters:                   $iters")
println("Largest eigen value:           $eigval")
println("Corresponding eigen vector:    $eigvec")
println()

jeigval, jeigvec = eigen(a);

println("\nEigen values:                   $jeigval")
println("Corresponding eigen vectors:")
jeigvec |> display

@test round(maximum(jeigval), digits=3) == round(eigval, digits=3)


text:
using DSGE, ModelConstructors, Dates, Test

@testset "sys1: no rule-switch, sys2 parameter-switch in 1st forecast quarter but b/c no rule switch, implied system is the same" begin
    m = Model1002("ss10")
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) # Parameter and/or rule switch in 2nd forecast quarter
    ModelConstructors.set_regime_val!(m[:α], 1, m[:α].value)
    ModelConstructors.set_regime_val!(m[:α], 2, m[:α].value)
    ModelConstructors.set_regime_val!(m[:α], 3, m[:α].value / 2.)
    ModelConstructors.set_regime_val!(m[:α], 4, m[:α].value)

    sys1 = compute_system(m)

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)
    m <= Setting(:gensys_regimes, [1:2])
    m <= Setting(:gensys2_regimes, [2:4])
    sys2 = compute_system(m)

    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] != sys2[3, :TTT]
    @test sys1[4, :TTT] == sys2[4, :TTT]
end

@testset "sys1: no rule-switch, no parameter-switch. sys2: rule- and parameter-switch in 1st forecast quarter (so system 3 should be different)" begin
    m = Model1002("ss10")
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:α], 1, m[:α].value) #2.523)
    ModelConstructors.set_regime_val!(m[:α], 2, m[:α].value)
    ModelConstructors.set_regime_val!(m[:α], 3, 0.00001) #
    ModelConstructors.set_regime_val!(m[:α], 4, m[:α].value)

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)
    m <= Setting(:gensys_regimes, [1:2])
    m <= Setting(:gensys2_regimes, [2:4])

    sys2 = compute_system(m)

    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] != sys2[3, :TTT]
    @test sys1[4, :TTT] == sys2[4, :TTT]
end

@testset "sys1: no rule-switch, no parameter-switch. sys2: rule-switch in 1st forecast quarter, parameter-switch in 2nd forecastquarter" begin
    m = Model1002("ss10")
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:α], 1, m[:α].value) #2.523)
    ModelConstructors.set_regime_val!(m[:α], 2, m[:α].value)
    ModelConstructors.set_regime_val!(m[:α], 3, m[:α].value)
    ModelConstructors.set_regime_val!(m[:α], 4, 0.00001) #

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)
    m <= Setting(:gensys_regimes, [1:2])
    m <= Setting(:gensys2_regimes, [2:4])
    sys2 = compute_system(m)

    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] != sys2[3, :TTT]
    @test sys1[4, :TTT] != sys2[4, :TTT]
end

@testset "sys1: no rule-switch, no parameter-switch. sys2: no rule-switch, yes parameter-switch (in 1st forecast quarter)" begin
    m = Model1002("ss10")
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:α], 1, m[:α].value) #2.523)
    ModelConstructors.set_regime_val!(m[:α], 2, m[:α].value)
    ModelConstructors.set_regime_val!(m[:α], 3, 0.00001) #
    ModelConstructors.set_regime_val!(m[:α], 4, m[:α].value)

    sys2 = compute_system(m)

    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] != sys2[3, :TTT]
    @test sys1[4, :TTT] == sys2[4, :TTT]

end

### Now with QQ
@testset "sys1: rule-switch (1st forecast quarter), no parameter-switch. sys2: rule-switch (1st forecast quarter), no parameter-switch (parameter switch with same values)" begin
    m = Model1002("ss10")

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:σ_g], 1, m[:σ_g].value) #2.523)
    ModelConstructors.set_regime_val!(m[:σ_g], 2, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 3, m[:σ_g].value) #
    ModelConstructors.set_regime_val!(m[:σ_g], 4, m[:σ_g].value)
    sys2 = compute_system(m)

    @test sys1[1, :QQ] == sys2[1, :QQ]
    @test sys1[2, :QQ] == sys2[2, :QQ]
    @test sys1[3, :QQ] == sys2[3, :QQ]
    @test sys1[4, :QQ] == sys2[4, :QQ]
end

@testset "sys1: no rule-switch, no parameter-switch. sys2: rule- and parameter-switch in 1st forecast quarter (so system 3 should be different)" begin
    m = Model1002("ss10")

    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:σ_g], 1, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 2, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 3, 0.00001)
    ModelConstructors.set_regime_val!(m[:σ_g], 4, m[:σ_g].value)

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)

    sys2 = compute_system(m)

    @test sys1[1, :QQ] == sys2[1, :QQ]
    @test sys1[2, :QQ] == sys2[2, :QQ]
    @test sys1[3, :QQ] != sys2[3, :QQ]
    @test sys1[4, :QQ] == sys2[4, :QQ]
end

@testset "sys1: no rule-switch, sys2 parameter-switch in shocks in 1st forecast quarter so only QQ changes" begin
    # rule switching doesn't affect QQ here
    m = Model1002("ss10")

    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:σ_g], 1, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 2, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 3, m[:σ_g].value / 2.)
    ModelConstructors.set_regime_val!(m[:σ_g], 4, m[:σ_g].value)

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)
    m <= Setting(:gensys_regimes, [1:2])
    m <= Setting(:gensys2_regimes, [2:4])
    sys2 = compute_system(m)

    @test sys1[1, :QQ] == sys2[1, :QQ]
    @test sys1[2, :QQ] == sys2[2, :QQ]
    @test sys1[3, :QQ] != sys2[3, :QQ]
    @test sys1[4, :QQ] == sys2[4, :QQ]
    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] ≈ sys2[3, :TTT] # populated by gensys2, so slight numerical differences
    @test sys1[4, :TTT] == sys2[4, :TTT]
end

@testset "sys1: no rule-switch, no parameter-switch. sys2: rule-switch in 1st forecast quarter, parameter-switch in 2nd forecastquarter" begin
    # rule switching doesn't affect QQ here
    m = Model1002("ss10")

    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:σ_g], 1, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 2, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 3, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 4, 0.00001)

    # Use gensys2 for anticipated policy rule change
    m <= Setting(:gensys2, true)
    m <= Setting(:gensys_regimes, [1:2])
    m <= Setting(:gensys2_regimes, [2:4])
    sys2 = compute_system(m)

    @test sys1[1, :QQ] == sys2[1, :QQ]
    @test sys1[2, :QQ] == sys2[2, :QQ]
    @test sys1[3, :QQ] == sys2[3, :QQ]
    @test sys1[4, :QQ] != sys2[4, :QQ]
    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] ≈ sys2[3, :TTT] # populated by gensys2, so slight numerical differences
    @test sys1[4, :TTT] == sys2[4, :TTT]
end

@testset "sys1: no rule-switch, no parameter-switch. sys2: no rule-switch, yes parameter-switch (in 1st forecast quarter)" begin
    m = Model1002("ss10")
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)
    # Number of regimes OVERALL
    m <= Setting(:n_regimes, 4)
    # Number of regimes in forcast
    m <= Setting(:n_fcast_regimes, 2)
    # Number of regimes in history
    m <= Setting(:n_hist_regimes, 2)
    # Dates of all regime switches (including both rule and parameter switches)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30))) #Parameter and/or rule switch in 2nd forecast quarter

    sys1 = compute_system(m)

    ModelConstructors.set_regime_val!(m[:σ_g], 1, m[:σ_g].value) #2.523)
    ModelConstructors.set_regime_val!(m[:σ_g], 2, m[:σ_g].value)
    ModelConstructors.set_regime_val!(m[:σ_g], 3, 0.00001) #
    ModelConstructors.set_regime_val!(m[:σ_g], 4, m[:σ_g].value)

    sys2 = compute_system(m)

    @test sys1[1, :QQ] == sys2[1, :QQ]
    @test sys1[2, :QQ] == sys2[2, :QQ]
    @test sys1[3, :QQ] != sys2[3, :QQ]
    @test sys1[4, :QQ] == sys2[4, :QQ]
    @test sys1[1, :TTT] == sys2[1, :TTT]
    @test sys1[2, :TTT] == sys2[2, :TTT]
    @test sys1[3, :TTT] == sys2[3, :TTT]
    @test sys1[4, :TTT] == sys2[4, :TTT]
end

@testset "Adding redundant regimes does not affect gensys2" begin
    # Add additional check for using replace_eqcond "falsely"
    m = Model1002("ss10")
    # Always need regime switching on, even if you're not switching parameters
    m <= Setting(:regime_switching, true)

    m <= Setting(:date_forecast_start, Date(2020, 9, 30))
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30),
                                     5 => Date(2020, 12, 31),
                                     6 => Date(2021, 3, 31))) #Parameter and/or rule switch in 2nd forecast quarter
    setup_regime_switching_inds!(m)

    m <= Setting(:gensys2, true)
    m <= Setting(:replace_eqcond, true)
    m <= Setting(:regime_eqcond_info, Dict(4 => DSGE.EqcondEntry(DSGE.zlb_rule(), [1., 0.])))
    sys1 = compute_system(m)
    m <= Setting(:regime_dates, Dict(1 => date_presample_start(m), #Start
                                     2 => Date(1990, 3, 31), # Parameter switch in history
                                     3 => Date(2020, 6, 30), # Parameter and/or rule switch in 1st forecast quarter
                                     4 => Date(2020, 9, 30),
                                     5 => Date(2020, 12, 31),
                                     6 => Date(2021, 3, 31),
                                     7 => Date(2021, 6, 30))) #Parameter and/or rule switch in 2nd forecast quarter
    setup_regime_switching_inds!(m)
    sys2 = compute_system(m)

    for i in 1:n_regimes(sys1)
        @test sys1[i, :TTT] ≈ sys2[i, :TTT]
        @test sys1[i, :RRR] ≈ sys2[i, :RRR]
        @test sys1[i, :CCC] ≈ sys2[i, :CCC]
    end
    @test sys2[7, :TTT] ≈ sys1[6, :TTT]
    @test sys2[7, :RRR] ≈ sys1[6, :RRR]
    @test sys2[7, :CCC] ≈ sys1[6, :CCC]
end


text:
"""
Modia module with block component models (inspired from Modelica Standard Library.

* Developer: Hilding Elmqvist, Mogram AB
* Copyright (c) 2016-2018: Hilding Elmqvist, Toivo Henningsson, Martin Otter
* License: MIT (expat)

"""
module Blocks

using ..Instantiation
using Modia
export FirstOrder, Feedback, PI, Step, Sine, Switch, MIMO

@model Block begin
end

@model SISO begin
    # Single Input Single Output continuous control block
  @extends Block()
  u=Float()
  y=Float()
end

@model FirstOrder begin
    # First order transfer function block (= 1 pole)
  k=1   # Gain
  T=1   # Time Constant

  @extends SISO()
  @inherits u, y
@equations begin
  der(y) = (k*u - y)/T
  end
end

@model Feedback begin
    # Output difference between commanded and feedback input
  u1=Float()
  u2=Float()
  y=Float()
@equations begin
  y = u1 - u2
  end
end

@model PI begin
    # Proportional-Integral controller
  k=1   # Gain
  T=Parameter(1, min=1E-10)   # Time Constant (T>0 required)

  @extends SISO()
  @inherits u, y
  x=Float(start=0)   # State of block
@equations begin
  der(x) = u/T
  y = k*(x + u)
  end
end

@model SO begin
    # Single Output continuous control block
  @extends Block()
  y=Variable()
end

@model SignalSource begin
    # Base class for continuous signal source
  @extends SO()
  offset=0      # Offset of output signal y
  startTime=0   # Output y = offset for time < startTime
end

@model Step begin
    # Generate step signal of type Real
  height=1      # Height of step
  @extends SignalSource()
  @inherits y, offset, startTime
	t = Float(start=0.0)
@equations begin
  y = offset + (t < startTime ? 0 : height)
	der(t) = 1
  end
end

@model Sine begin
    # Generate sine signal
  amplitude=1 # Amplitude of sine wave
  freqHz=Parameter() # Frequency of sine wave
  phase=0 # Phase of sine wave
  offset=0      # Offset of output signal y
  startTime=0   # Output y = offset for time < startTime
  @extends SO()
  @inherits y
	t = Float(start=0.0)
@equations begin
#  y = offset + if time < startTime;  0 else amplitude*sin(2*pi*freqHz*(time - startTime) + phase) end
  y = offset + if t < startTime;  0 else amplitude*sin(2*pi*freqHz*(t - startTime) + phase) end
	der(t) = 1
  end
end

@model Sine2 begin
    # Generate sine signal
  amplitude=1 # Amplitude of sine wave
  freqHz=Parameter() # Frequency of sine wave
  phase=0 # Phase of sine wave
  @extends SignalSource()
  @inherits y, offset, startTime
@equations begin
  y = offset + if time < startTime;  0 else amplitude*sin(2*pi*freqHz*(time - startTime) + phase) end
  end
end

@model Switch begin
  sw=Boolean()
  u1=Variable(); u2=Variable()
  y=Variable()
@equations begin
  y = if sw; u1 else u2 end
  end
end

@model ABCD begin
  A=-1; B=1; C=1; D=0
  u=Float(); y=Float()
  x=Float(start=0)
@equations begin
  der(x) = A*x + B*u
  y = C*x + D*u
  end
end

end


text:
@testset "xnor constraint violated or solved?" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, !((sum(x) <= 2) ⊻ (sum(x) > 3)))
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 2; check_feasibility = false)
    @test CS.is_constraint_solved(com, constraint, constraint.fct, constraint.set)

    #################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, !((sum(x) < 2) ⊻ (sum(x) > 4)))
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[1]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.is_constraint_solved(com, constraint, constraint.fct, constraint.set)

    #################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((sum(x) > 3) ⊻ (sum(x) < 2))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    xor_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[xor_constraint.indices[1]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[xor_constraint.indices[2]], 2; check_feasibility = false)
    @test CS.is_constraint_solved(com, xor_constraint, xor_constraint.fct, xor_constraint.set)

    ###############################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((sum(x) > 3) ⊻ (sum(x) < 2))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    xor_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[xor_constraint.indices[1]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[xor_constraint.indices[2]], 4; check_feasibility = false)
    @test !CS.is_constraint_solved(com, xor_constraint, xor_constraint.fct, xor_constraint.set)

    ##################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((sum(x) > 3) ⊻ (sum(x) < 2))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    ##################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((sum(x) > 3) ⊻ (x in CS.AllDifferent()))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)

    ##################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((sum(x) > 3) ⊻ (x in CS.TableSet([0 1; 1 0])))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end

@testset "Indicator fixed to 1 xnor constraint violated or solved?" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {!((sum(x) > 3) ⊻ (sum(x) < 2))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_solved(com, constraint, constraint.fct, constraint.set)

    #################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {!((sum(x) < 3) ⊻ (sum(x) >= 2))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 1; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 1; check_feasibility = false)
    @test CS.is_constraint_solved(com, constraint, constraint.fct, constraint.set)

    #################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {!((sum(x) < 3) ⊻ (sum(x) > 1))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    xor_constraint = com.constraints[1].inner_constraint
    @test CS.fix!(com, variables[xor_constraint.indices[1]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[xor_constraint.indices[2]], 2; check_feasibility = false)
    @test CS.is_constraint_solved(com, xor_constraint, xor_constraint.fct, xor_constraint.set)

    ##################

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b => {!((sum(x) < 3) ⊻ (sum(x) > 1))})
    optimize!(m)
    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    @test CS.fix!(com, variables[constraint.indices[2]], 0; check_feasibility = false)
    @test CS.fix!(com, variables[constraint.indices[3]], 0; check_feasibility = false)
    @test CS.is_constraint_violated(com, constraint, constraint.fct, constraint.set)
end

@testset "reified xnor constraint prune_constraint!" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((x[1] <= 2) ⊻ (x[2] <= 2))})
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]

    constr_indices = constraint.indices
    @test CS.fix!(com, variables[constraint.indices[3]], 0; check_feasibility = false)
    CS.changed!(com, constraint, constraint.fct, constraint.set)
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test sort(CS.values(m, x[1])) == [0,1,2]
end

@testset "reified xnor constraint still feasible" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, b >= 1, Bin)
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, b := {!((x[1] <= 2) ⊻ (x[2] <= 2))})
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    constraint = com.constraints[1]
    xor_constraint = com.constraints[1].inner_constraint

    constr_indices = xor_constraint.indices
    @test CS.prune_constraint!(com, constraint, constraint.fct, constraint.set)
    @test CS.still_feasible(com, xor_constraint, xor_constraint.fct, xor_constraint.set, xor_constraint.indices[2], 0)
    @test CS.still_feasible(com, xor_constraint, xor_constraint.fct, xor_constraint.set, xor_constraint.indices[2], 3)
end

@testset "xnor constraint still feasible" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1]+x[2] <= 2) ⊻ (x[2] <= 2)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    @test CS.still_feasible(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set, 1, 0)
    @test CS.fix!(com, variables[2], 1; check_feasibility = false)
    @test !CS.still_feasible(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set, 1, 3)

    ###

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1] <= 2) ⊻ (x[2]+x[1] <= 2)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    @test CS.still_feasible(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set, 2, 0)
    @test CS.fix!(com, variables[1], 1; check_feasibility = false)
    @test !CS.still_feasible(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set, 2, 3)
end

@testset "xnor constraint prune_constraint" begin
    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 2 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1]+x[2] < 2) ⊻ (x[2] < 2)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    # is already solved as both are violated no pruning possible
    @test CS.prune_constraint!(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set)

    ###

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 2 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1]+x[2] >= 2) ⊻ (x[2] >= 2)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    # is already solved as both are subconstraints are solved => no pruning possible
    @test CS.prune_constraint!(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set)

    ###

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 2 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1]+x[2] >= 2) ⊻ (x[2] > 4)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    # lhs solved => rhs must be solved
    @test CS.prune_constraint!(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set)
    @test CS.values(variables[2]) == [5]

    ###

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 2 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1]+x[2] <= 2) ⊻ (x[2] > 4)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    # lhs violated => rhs must be complement solved
    @test CS.prune_constraint!(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set)
    @test sort!(CS.values(variables[2])) == [2,3,4]

    ###

    m = Model(optimizer_with_attributes(CS.Optimizer, "no_prune" => true, "logging" => []))
    @variable(m, 0 <= x[1:2] <= 5, Int)
    @constraint(m, !((x[1]+x[2] <= 5) ⊻ (x[2] > 5)))
    optimize!(m)

    com = CS.get_inner_model(m)

    variables = com.search_space
    xnor_constraint = com.constraints[1]

    constr_indices = xnor_constraint.indices
    # rhs violated => lhs must be complement solved
    @test CS.prune_constraint!(com, xnor_constraint, xnor_constraint.fct, xnor_constraint.set)
    # the variables can't be 0 as otherwise we can't fullfil the complement x[1]+x[2] > 5
    @test !CS.has(variables[1], 0)
    @test !CS.has(variables[2], 0)
end


text:
module MitochondrialDynamics

using ModelingToolkit

export make_model, DEFAULT_U0

include("utils.jl")
include("mtk.jl")

end # module

max_stars_repo_name: ntumitolab/mitodyn-ode

max_issues_repo_name: ntumitolab/mitodyn-ode

max_forks_repo_name: NTUMitoLab/mitodyn-ode


text:
#=
ndgrid-avect.jl
ndgrid for a set of AbstractVector inputs.
=#

export ndgrid

"""
    GridAV{T,d,D} <: AbstractGrid{T,d,D}
The `d`th component of `D`-dimensional `ndgrid(x, y, ...)`
where `1 ≤ d ≤ D` and `x, y, ...` are each an `AbstractVector`.
"""
struct GridAV{T,d,D} <: AbstractGrid{T,d,D}
    dims::Dims{D}
    v::AbstractVector{T}

    function GridAV(dims::Dims{D}, v::AbstractVector{T}, d::Int) where {D, T}
        1 ≤ d ≤ D || throw(ArgumentError("$d for $dims"))
        new{T,d,D}(dims, v)
    end
end

Base.size(a::GridAV) = a.dims
Base.eltype(::GridAV{T}) where T = T

@inline Base.@propagate_inbounds function Base.getindex(
    a::GridAV{T,d,D},
    i::Vararg{Int,D},
) where {T, d, D}
    @boundscheck checkbounds(a, i...)
#   @boundscheck checkbounds(i, d)
    @boundscheck checkbounds(a.v, i[d])
    return @inbounds a.v[ @inbounds i[d] ]
end

_grid(d::Int, dims::Dims, v::Base.OneTo{T}) where T = GridOT(T, dims, d)
_grid(d::Int, dims::Dims, v::UnitRange)      = GridUR(dims, v, d)
_grid(d::Int, dims::Dims, v::AbstractRange)  = GridAR(dims, v, d)
_grid(d::Int, dims::Dims, v::AbstractVector) = GridAV(dims, v, d)

_grid_type(d::Int, D::Int, ::Base.OneTo{T}) where T = GridOT{T, d, D}
_grid_type(d::Int, D::Int, ::UnitRange{T})  where T = GridUR{T, d, D}
_grid_type(d::Int, D::Int, ::AbstractRange{T})  where T = GridAR{T, d, D}
_grid_type(d::Int, D::Int, ::AbstractVector{T}) where T = GridAV{T, d, D}

"""
    (xg, yg, ...) = ndgrid(v1, v2, ...)
Construct `ndgrid` tuple for `AbstractVector` inputs.
Each output has a lazy grid type (subtype of `AbstractGrid`)
according to the corresponding input vector type.

# Examples

```jldoctest
julia> xg, yg = ndgrid(1:3, [:a, :b])
([1 1; 2 2; 3 3], [:a :b; :a :b; :a :b])

julia> xg
3×2 LazyGrids.GridUR{Int64, 1, 2}:
 1  1
 2  2
 3  3

julia> yg
3×2 LazyGrids.GridAV{Symbol, 2, 2}:
 :a  :b
 :a  :b
 :a  :b
```
"""
function ndgrid(vs::Vararg{AbstractVector})
    D = length(vs)
    dims = ntuple(d -> length(vs[d]), D)
    T = Tuple{ntuple(d -> _grid_type(d, D, vs[d]), D)...} # return type
    return ntuple(d -> _grid(d, dims, vs[d]), D)::T
end


text:
function GHIJjsPolynomials(arg0::jdouble, arg1::jdouble, arg2::jdouble, arg3::jdouble)
    return GHIJjsPolynomials((jdouble, jdouble, jdouble, jdouble), arg0, arg1, arg2, arg3)
end

function get_gjs(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getGjs", jdouble, (jint, jint), arg0, arg1)
end

function get_hjs(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getHjs", jdouble, (jint, jint), arg0, arg1)
end

function get_ijs(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getIjs", jdouble, (jint, jint), arg0, arg1)
end

function get_jjs(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getJjs", jdouble, (jint, jint), arg0, arg1)
end

function getd_gjsd_alpha(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdGjsdAlpha", jdouble, (jint, jint), arg0, arg1)
end

function getd_gjsd_beta(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdGjsdBeta", jdouble, (jint, jint), arg0, arg1)
end

function getd_gjsdh(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdGjsdh", jdouble, (jint, jint), arg0, arg1)
end

function getd_gjsdk(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdGjsdk", jdouble, (jint, jint), arg0, arg1)
end

function getd_hjsd_alpha(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdHjsdAlpha", jdouble, (jint, jint), arg0, arg1)
end

function getd_hjsd_beta(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdHjsdBeta", jdouble, (jint, jint), arg0, arg1)
end

function getd_hjsdh(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdHjsdh", jdouble, (jint, jint), arg0, arg1)
end

function getd_hjsdk(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdHjsdk", jdouble, (jint, jint), arg0, arg1)
end

function getd_ijsd_alpha(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdIjsdAlpha", jdouble, (jint, jint), arg0, arg1)
end

function getd_ijsd_beta(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdIjsdBeta", jdouble, (jint, jint), arg0, arg1)
end

function getd_ijsdh(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdIjsdh", jdouble, (jint, jint), arg0, arg1)
end

function getd_ijsdk(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdIjsdk", jdouble, (jint, jint), arg0, arg1)
end

function getd_jjsd_alpha(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdJjsdAlpha", jdouble, (jint, jint), arg0, arg1)
end

function getd_jjsd_beta(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdJjsdBeta", jdouble, (jint, jint), arg0, arg1)
end

function getd_jjsdh(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdJjsdh", jdouble, (jint, jint), arg0, arg1)
end

function getd_jjsdk(obj::GHIJjsPolynomials, arg0::jint, arg1::jint)
    return jcall(obj, "getdJjsdk", jdouble, (jint, jint), arg0, arg1)
end


text:
# -*- coding: utf-8 -*-
import HELICS

const h = HELICS

initstring = "-f 2 --name=mainbroker"
fedinitstring = "--broker=mainbroker --federates=1"
deltat = 0.01

helicsversion = h.helicsGetVersion()

println("PI SENDER: Helics version = $helicsversion")

# Create broker #
println("Creating Broker")
broker = h.helicsCreateBroker("zmq", "", initstring)
println("Created Broker")

println("Checking if Broker is connected")
isconnected = h.helicsBrokerIsConnected(broker)
println("Checked if Broker is connected")

if isconnected == 1
    println("Broker created and connected")
end

# Create Federate Info object that describes the federate properties #
fedinfo = h.helicsCreateFederateInfo()

# Set Federate name #
h.helicsFederateInfoSetCoreName(fedinfo, "TestA Federate")

# Set core type from string #
h.helicsFederateInfoSetCoreTypeFromString(fedinfo, "zmq")

# Federate init string #
h.helicsFederateInfoSetCoreInitString(fedinfo, fedinitstring)

# Set the message interval (timedelta) for federate. Note th#
# HELICS minimum message time interval is 1 ns and by default
# it uses a time delta of 1 second. What is provided to the
# setTimedelta routine is a multiplier for the default timedelta.

# Set one second message interval #
h.helicsFederateInfoSetTimeProperty(fedinfo, h.HELICS_PROPERTY_TIME_DELTA, deltat)

# Create value federate #
vfed = h.helicsCreateValueFederate("TestA Federate", fedinfo)
println("PI SENDER: Value federate created")

# Register the publication #
pub = h.helicsFederateRegisterGlobalTypePublication(vfed, "testA", "double", "")
println("PI SENDER: Publication registered")

# Enter execution mode #
h.helicsFederateEnterExecutingMode(vfed)
println("PI SENDER: Entering execution mode")

# This federate will be publishing deltat*pi for numsteps steps #
this_time = 0.0
value = pi

for t in 5:10
    val = value

    currenttime = h.helicsFederateRequestTime(vfed, Float64(t))

    h.helicsPublicationPublishDouble(pub, Float64(val))
    println(
        "PI SENDER: Sending value pi = $val at time $currenttime to PI RECEIVER"
    )

    sleep(1)
end

h.helicsFederateFinalize(vfed)
println("PI SENDER: Federate finalized")

while h.helicsBrokerIsConnected(broker)
    sleep(1)
end

h.helicsFederateFree(vfed)
h.helicsCloseLibrary()

println("PI SENDER: Broker disconnected")


text:
# Run tests on functions in src/by_reference_doc/shows.jl

using Spotify.Shows

@testset verbose = true "GET-request endpoints for shows" begin

    # Input argument from composite types defined in src/types.jl
    show_id = SpShowId()

    # Cycle through different input keywords for testing
    markets = ["US", "NL", "DE"]
    offsets = [10, 35]

    @testset "For market = $(market_id)" for market_id in markets

        @test ~isempty(show_get_single(show_id, market = "$market_id")[1])

        @test ~isempty(show_get_multiple("$show_id, 4rOoJ6Egrf8K2IrywzwOMk", market = "$market_id")[1])

        @test ~isempty(show_get_episodes(show_id, market = "$market_id")[1])

    end

    @test ~isempty(show_get_saved()[1])

    # Currently shows 403 error, implemented in shows.jl
    @test_broken ~isempty(show_get_contains("2MAi0BvDc6GTFvKFPXnkCL")[1])

end


text:
export ForwardModel, Loss, phi, localDescent, solveFiniteDimProblem, lmo, loss
abstract type ForwardModel end
abstract type Loss end
# Computes the forward operator.
# Inputs:
#   parameters:: p by k matrix. You might want a generic list instead if the parameter space is more complex.
#   weights:: p-dimensional vector of weights.
#
#   output:: d-dimensional vector.
#   $o = \sum_{i=1}^k w_i \psi(\theta_i).$
phi(model :: ForwardModel, parameters :: Matrix{Float64}, weights :: Vector{Float64}) =
   error("phi not implemented for model $(typeof(model)).")

# Perform some local optimization of
#   $(theta_1, \ldots, \theta_k) \mapsto \ell \left(\sum_{i=1}^k w_i \psf(\theta_i) \right).$
# Hopefully uses differentiability of $\psi.$
#
# Inputs:
#   loss :: A loss function.∆˚†
#   weights :: k-dimensional vector of (fixed) weights.
#   y :: d-dimensional target vector.
#
localDescent(model :: ForwardModel, loss :: Loss, thetas :: Matrix{Float64}, weights :: Vector{Float64}, y :: Vector{Float64})  =
   error("localDescent not implemented for model $(typeof(model)) and loss $(typeof(loss)).")

# Solve the (convex) optimization problem:
# $$ \min_{0 \le w \in \mathbb{R}^k} \ell \left(\sum_i w_i \psi(\theta_i) - y \right) .$$
solveFiniteDimProblem(model :: ForwardModel, loss :: Loss, parameters :: Array{Float64, 2}, y :: Array{Float64,1}, tau :: Float64)  =
   error("phi not implemented for model $(typeof(model)) and loss $(typeof(loss)).")

# Get the next atom.
# i.e. Find $ \arg\max_{\theta \in \Theta} | \langle\psi(\theta), v \rangle |.$
#
# If the model requires weights to be nonnegative, we instead need to find
# $ \arg\min_{\theta \in \Theta}  \langle\psi(\theta), v \rangle .$
#
# Returns the tuple (theta, score), where score is $-|\langle \psi(\theta), v \rangle$
lmo(model :: ForwardModel, v :: Vector{Float64})  =
   error("lmo not implemented for model $(typeof(model)).")
#Return the tuple (v, g) where $v = \ell(r)$ and $g = \nabla \ell(r).$
loss(loss :: Loss, r :: Vector{Float64}) = error("loss not implemented for loss $(typeof(loss)).")


text:
"""
`DAE`: Differential Algebraic Equation

Defines a differential algebraic initial value problem
```math
\\begin{align*}
\\dot{q} (t) &= f(t, q(t)) + u(t, q(t), \\lambda(t)) , & q(t_{0}) &= q_{0} , \\\\
0 &= \\phi (t, q(t), \\lambda(t)) , & \\lambda(t_{0}) &= \\lambda_{0} ,
\\end{align*}
```
with vector field ``f``, projection ``u``, algebraic constraint ``\\phi=0``,
initial conditions ``q_{0}`` and ``\\lambda_{0}``, the dynamical variable ``q``
taking values in ``\\mathbb{R}^{m}`` and the algebraic variable ``\\lambda``
taking values in ``\\mathbb{R}^{n}``.

### Fields

* `m`: dimension of dynamical variable ``q`` and the vector field ``f``
* `n`: dimension of algebraic variable ``\\lambda`` and the constraint ``\\phi``
* `f`: function computing the vector field
* `u`: function computing the projection
* `ϕ`: algebraic constraint
* `t₀`: initial time
* `q₀`: initial condition for dynamical variable ``q``
* `λ₀`: initial condition for algebraic variable ``\\lambda``

The function `f`, providing the vector field, takes three arguments,
`f(t, q, fq)`, the functions `u` and `ϕ`, providing the projection and the
algebraic constraint take four arguments, `u(t, q, λ, fu)` and `ϕ(t, q, λ, fϕ)`,
where `t` is the current time, `q` and `λ` are the current solution vectors,
and `fq`, `fu` and `fϕ` are the vectors which hold the result of evaluating the
vector field ``f``, the projection ``u`` and the algebraic constraint ``\\phi``
on `t`, `q` and `λ`.

### Example

```julia
    function f(t, q, fq)
        fq[1] = q[1]
        fq[2] = q[2]
    end

    function u(t, q, λ, fu)
        fu[1] = +λ[1]
        fu[2] = -λ[1]
    end

    function ϕ(t, q, λ, fϕ)
        fϕ[1] = q[2] - q[1]
    end

    t₀ = 0.
    q₀ = [1., 1.]
    λ₀ = [0.]

    dae = DAE(f, u, ϕ, t₀, q₀, λ₀)

```
"""
immutable DAE{T} <: Equation{T}
    m::Int
    n::Int
    f::Function
    u::Function
    ϕ::Function
    t₀::T
    q₀::Array{T,1}
    λ₀::Array{T,1}

    function DAE(m, n, f, u, ϕ, t₀, q₀, λ₀)
        @assert m == length(q₀)
        @assert n == length(λ₀)
        @assert m ≥ n
        @assert T == eltype(q₀)
        @assert T == eltype(λ₀)

        new(m, n, f, u, ϕ, t₀, q₀, λ₀)
    end
end

function DAE{T}(f::Function, u::Function, ϕ::Function, t₀::Real, q₀::Vector{T}, λ₀::Vector{T})
    DAE{T}(length(q₀), length(λ₀), f, u, ϕ, t₀, q₀, λ₀)
end

function DAE{T}(f::Function, u::Function, ϕ::Function, q₀::Vector{T}, λ₀::Vector{T})
    DAE(f, u, ϕ, 0, q₀, λ₀)
end


text:
# Autogenerated wrapper script for TOPCOM_jll for i686-w64-mingw32-cxx11
export chiro2allfinetriangs, chiro2alltriangs, chiro2circuits, chiro2cocircuits, chiro2dual, chiro2finetriang, chiro2finetriangs, chiro2mintriang, chiro2nallfinetriangs, chiro2nalltriangs, chiro2nfinetriangs, chiro2ntriangs, chiro2placingtriang, chiro2triangs, cocircuits2facets, points2allfinetriangs, points2alltriangs, points2chiro, points2facets, points2finetriang, points2finetriangs, points2flips, points2nallfinetriangs, points2nalltriangs, points2nfinetriangs, points2nflips, points2ntriangs, points2placingtriang, points2triangs, points2volume, santos_22_triang, santos_dim4_triang, santos_triang

using GMP_jll
using cddlib_jll
JLLWrappers.@generate_wrapper_header("TOPCOM")
JLLWrappers.@declare_executable_product(chiro2allfinetriangs)
JLLWrappers.@declare_executable_product(chiro2alltriangs)
JLLWrappers.@declare_executable_product(chiro2circuits)
JLLWrappers.@declare_executable_product(chiro2cocircuits)
JLLWrappers.@declare_executable_product(chiro2dual)
JLLWrappers.@declare_executable_product(chiro2finetriang)
JLLWrappers.@declare_executable_product(chiro2finetriangs)
JLLWrappers.@declare_executable_product(chiro2mintriang)
JLLWrappers.@declare_executable_product(chiro2nallfinetriangs)
JLLWrappers.@declare_executable_product(chiro2nalltriangs)
JLLWrappers.@declare_executable_product(chiro2nfinetriangs)
JLLWrappers.@declare_executable_product(chiro2ntriangs)
JLLWrappers.@declare_executable_product(chiro2placingtriang)
JLLWrappers.@declare_executable_product(chiro2triangs)
JLLWrappers.@declare_executable_product(cocircuits2facets)
JLLWrappers.@declare_executable_product(points2allfinetriangs)
JLLWrappers.@declare_executable_product(points2alltriangs)
JLLWrappers.@declare_executable_product(points2chiro)
JLLWrappers.@declare_executable_product(points2facets)
JLLWrappers.@declare_executable_product(points2finetriang)
JLLWrappers.@declare_executable_product(points2finetriangs)
JLLWrappers.@declare_executable_product(points2flips)
JLLWrappers.@declare_executable_product(points2nallfinetriangs)
JLLWrappers.@declare_executable_product(points2nalltriangs)
JLLWrappers.@declare_executable_product(points2nfinetriangs)
JLLWrappers.@declare_executable_product(points2nflips)
JLLWrappers.@declare_executable_product(points2ntriangs)
JLLWrappers.@declare_executable_product(points2placingtriang)
JLLWrappers.@declare_executable_product(points2triangs)
JLLWrappers.@declare_executable_product(points2volume)
JLLWrappers.@declare_executable_product(santos_22_triang)
JLLWrappers.@declare_executable_product(santos_dim4_triang)
JLLWrappers.@declare_executable_product(santos_triang)
function __init__()
    JLLWrappers.@generate_init_header(GMP_jll, cddlib_jll)
    JLLWrappers.@init_executable_product(
        chiro2allfinetriangs,
        "bin\\chiro2allfinetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2alltriangs,
        "bin\\chiro2alltriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2circuits,
        "bin\\chiro2circuits.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2cocircuits,
        "bin\\chiro2cocircuits.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2dual,
        "bin\\chiro2dual.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2finetriang,
        "bin\\chiro2finetriang.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2finetriangs,
        "bin\\chiro2finetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2mintriang,
        "bin\\chiro2mintriang.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2nallfinetriangs,
        "bin\\chiro2nallfinetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2nalltriangs,
        "bin\\chiro2nalltriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2nfinetriangs,
        "bin\\chiro2nfinetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2ntriangs,
        "bin\\chiro2ntriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2placingtriang,
        "bin\\chiro2placingtriang.exe",
    )

    JLLWrappers.@init_executable_product(
        chiro2triangs,
        "bin\\chiro2triangs.exe",
    )

    JLLWrappers.@init_executable_product(
        cocircuits2facets,
        "bin\\cocircuits2facets.exe",
    )

    JLLWrappers.@init_executable_product(
        points2allfinetriangs,
        "bin\\points2allfinetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2alltriangs,
        "bin\\points2alltriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2chiro,
        "bin\\points2chiro.exe",
    )

    JLLWrappers.@init_executable_product(
        points2facets,
        "bin\\points2facets.exe",
    )

    JLLWrappers.@init_executable_product(
        points2finetriang,
        "bin\\points2finetriang.exe",
    )

    JLLWrappers.@init_executable_product(
        points2finetriangs,
        "bin\\points2finetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2flips,
        "bin\\points2flips.exe",
    )

    JLLWrappers.@init_executable_product(
        points2nallfinetriangs,
        "bin\\points2nallfinetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2nalltriangs,
        "bin\\points2nalltriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2nfinetriangs,
        "bin\\points2nfinetriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2nflips,
        "bin\\points2nflips.exe",
    )

    JLLWrappers.@init_executable_product(
        points2ntriangs,
        "bin\\points2ntriangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2placingtriang,
        "bin\\points2placingtriang.exe",
    )

    JLLWrappers.@init_executable_product(
        points2triangs,
        "bin\\points2triangs.exe",
    )

    JLLWrappers.@init_executable_product(
        points2volume,
        "bin\\points2volume.exe",
    )

    JLLWrappers.@init_executable_product(
        santos_22_triang,
        "bin\\santos_22_triang.exe",
    )

    JLLWrappers.@init_executable_product(
        santos_dim4_triang,
        "bin\\santos_dim4_triang.exe",
    )

    JLLWrappers.@init_executable_product(
        santos_triang,
        "bin\\santos_triang.exe",
    )

    JLLWrappers.@generate_init_footer()
end  # __init__()


text:
abstract type BinaryNode{T} <: AbstractNode{T} end

### Interface

# Return the node's value.
getvalue(n::BinaryNode) = nothing

# Set the node's value.
setvalue!(n::BinaryNode) = nothing

# Return node's left child.
left(n::BinaryNode) = nothing

# Set the node' left child.
setleft!(n::BinaryNode, node) = nothing

# Return node's left child.
right(n::BinaryNode) = nothing

# Set the node' left child.
setright!(n::BinaryNode, node) = nothing

# Check if node is a leaf.
isleaf(n::BinaryNode) = nothing

# generic traversal and others
# function traversal() end

# preorder, inorder, postorder, levelorder


text:
@reexport module BlackJack

export BlackJackEnv

using ReinforcementLearningBase
using Random

const ACTIONS = [:hit, :stick]

const INDS = LinearIndices((
    2,  # player has ace or not
    11, # player's hands sum, 12 ~ 21, plus other cases that the sum is greater than 21
    10,  # dealer's hands sum, 2~11
))

deal_card() = rand(1:13)
is_ace(x) = x == 1
value(x) = is_ace(x) ? 11 : min(x, 10)

mutable struct Hands
    sum::Int
    cards::Vector{Int}
    n_usable_ace::Int
end

Hands() = Hands(0, [], false)

is_busted(h::Hands) = h.sum > 21

function Base.push!(h::Hands, x)
    is_busted(h) && throw(ArgumentError("cards in hand are already busted!"))
    push!(h.cards, x)
    h.sum += value(x)
    h.n_usable_ace += is_ace(x)

    if h.n_usable_ace > 0 && is_busted(h)
        h.sum -= 10
        h.n_usable_ace -= 1
    end
end

mutable struct BlackJackEnv <: AbstractEnv
    player_hands::Hands
    dealer_hands::Hands
    is_end::Bool
    reward::Float64
    is_exploring_start::Bool
    init::Union{Nothing,Tuple{Hands,Hands}}
end

RLBase.get_observation_space(env::BlackJackEnv) = DiscreteSpace(length(INDS))
RLBase.get_action_space(env::BlackJackEnv) = DiscreteSpace(2)

function BlackJackEnv(; is_exploring_start = false, init = nothing)
    env = BlackJackEnv(
        Hands(),
        Hands(),
        false,
        0.0,
        is_exploring_start,
        init,
    )
    init_hands!(env)
    env
end

function init_hands!(env::BlackJackEnv)
    player_hands, dealer_hands = Hands(), Hands()

    if env.is_exploring_start
        player_hands.sum = rand(12:21)
        player_hands.n_usable_ace = rand(Bool)
        dealer_hands.sum = rand(2:11)
        dealer_hands.n_usable_ace = dealer_hands.sum == 11 ? true : false
    elseif env.init === nothing
        while player_hands.sum <= 11
            push!(player_hands, deal_card())
        end
        push!(dealer_hands, deal_card())
    else
        player_hands, dealer_hands = deepcopy(env.init)
    end

    env.player_hands, env.dealer_hands = player_hands, dealer_hands
end

function (env::BlackJackEnv)(a::Int)
    if ACTIONS[a] == :hit
        push!(env.player_hands, deal_card())
        if is_busted(env.player_hands)
            env.is_end = true
            env.reward = -1.0
        else
            env.is_end = false
            env.reward = 0.0
        end
    elseif ACTIONS[a] == :stick
        while env.dealer_hands.sum < 17
            push!(env.dealer_hands, deal_card())
        end
        if is_busted(env.dealer_hands)
            env.reward = 1.0
        else
            if env.player_hands.sum > env.dealer_hands.sum
                env.reward = 1.0
            elseif env.player_hands.sum < env.dealer_hands.sum
                env.reward = -1.0
            else
                env.reward = 0.0
            end
        end
        env.is_end = true
    end
    nothing
end

function RLBase.reset!(env::BlackJackEnv)
    env.is_end = false
    env.reward = 0.0

    init_hands!(env)

    nothing
end

encode(env) =
    INDS[
        env.player_hands.n_usable_ace > 0 ? 1 : 2,
        12 <= env.player_hands.sum <= 21 ? env.player_hands.sum - 10 : 1,
        2 <= env.dealer_hands.sum <= 10 ? env.dealer_hands.sum : 1,
    ]

RLBase.observe(env::BlackJackEnv) =
    (reward = env.reward, terminal = env.is_end, state = encode(env))

end


text:
using RandomNumbers.PCG
if !@isdefined RandomNumbers
    include("../common.jl")
end

@info "Testing PCG"
stdout_ = stdout
pwd_ = pwd()
cd(dirname(@__FILE__))
rm("./actual"; force=true, recursive=true)
mkpath("./actual")

numbers = ['A', '2', '3', '4', '5', '6', '7', '8', '9', 'T', 'J', 'Q', 'K']
suit = ['h', 'c', 'd', 's'];

PCGStateOneseq()
PCGStateOneseq(UInt32)
PCGStateMCG()
PCGStateMCG(UInt32)
PCGStateSetseq()
PCGStateSetseq(UInt32)
PCGStateUnique()
PCGStateUnique(UInt32)

for (state_type_t, uint_type, method_symbol, return_type) in PCG_LIST
    state_type = eval(Symbol("PCGState$state_type_t"))
    method = Val{method_symbol}

    state_type(method)
    state_type(return_type, method)

    if state_type_t == :Setseq
        r = state_type(return_type, method, (42, 54))
        @test seed_type(r) == NTuple{2, uint_type}
    else
        r = state_type(return_type, method, 42)
        @test seed_type(r) == uint_type
    end
    @test copyto!(copy(r), r) == r

    # Unique state won't produce the same sequence every time.
    if state_type_t == :Unique
        t = rand(r, return_type)
        advance!(r, -1)
        @test t == rand(r, return_type)
        bounded_rand(r, 200701281 % return_type)
        continue
    end
    outfile = open(string(
        "./actual/check-$(lowercase("$state_type_t"))-$(sizeof(uint_type) << 3)-",
        "$(replace(lowercase("$method_symbol"), "_" => "-"))-$(sizeof(return_type) << 3).out"
    ), "w")
    redirect_stdout(outfile)
    for round in 1:5
        @printf "Round %d:\n" round
        @printf "%4dbit:" (sizeof(return_type) << 3)
        values = 10
        wrap = 10
        printstr = " 0x%04x"
        if return_type == UInt8
            values, wrap, printstr = 14, 14, " 0x%02x"
        elseif return_type == UInt32
            values, wrap, printstr = 6, 6, " 0x%08x"
        elseif return_type == UInt64
            values, wrap, printstr = 6, 3, " 0x%016llx"
        elseif return_type == UInt128
            values, wrap, printstr = 6, 2, " 0x%032llx"
        end

        for i in 1:values
            if i > 1 && (i - 1) % wrap == 0
                @printf "\n\t"
            end
            value = rand(r, return_type)
            @eval @printf $printstr $value
        end
        @printf "\n"

        @printf "  Again:"
        advance!(r, -(values % uint_type))
        for i in 1:values
            if i > 1 && (i - 1) % wrap == 0
                @printf "\n\t"
            end
            value = rand(r, return_type)
            @eval @printf $printstr $value
        end
        @printf "\n"

        @printf "  Coins: "
        for i in 1:65
            @printf "%c" (bounded_rand(r, 2 % return_type) == 1 ? 'H' : 'T')
        end
        @printf "\n"

        @printf "  Rolls:"
        for i in 1:33
            @printf " %d" (UInt32(bounded_rand(r, 6 % return_type)) + 1);
        end
        @printf "\n"

        cards = collect(0:51)
        for i = 52:-1:2
            chosen = bounded_rand(r, i % return_type)
            card = cards[chosen+1]
            cards[chosen+1] = cards[i]
            cards[i] = card
        end

        @printf "  Cards:"
        for i = 0:51
            @printf " %c%c" numbers[(cards[i+1] ÷ 4)+1] suit[(cards[i+1] % 4)+1]
            if (i+1) % 22 == 0
                @printf "\n\t"
            end
        end
        @printf "\n"

        @printf "\n"
    end
    close(outfile)
end
redirect_stdout(stdout_)

compare_dirs("expected", "actual")
cd(pwd_)


text:
using Logging2
using Logging
using Test
using Test: Ignored

struct IOWriteRecord <: IO
    writes
end

IOWriteRecord() = IOWriteRecord(Vector{String}())

function Base.unsafe_write(io::IOWriteRecord, p::Ptr{UInt8}, n::UInt)
    push!(io.writes, unsafe_string(p, n))
end

@testset "LineBufferedIO" begin
    function buffer_lines(strs...)
        w = IOWriteRecord()
        lbuf = Logging2.LineBufferedIO(w)
        for s in strs
            write(lbuf, s)
        end
        w.writes
    end

    @test buffer_lines("a", "\n") == ["a\n"]
    @test buffer_lines('a', '\n') == ["a\n"]
    @test buffer_lines("a", "b", "\n") == ["ab\n"]
    @test buffer_lines("a", "b", "\n", "c") == ["ab\n"]
    @test buffer_lines("a", "b", "\n", "c", "d\n") == ["ab\n", "cd\n"]
    @test buffer_lines("\n", "\n") == ["\n", "\n"]

    @test buffer_lines("ab\nc", "d") == ["ab\n"]
    @test buffer_lines("ab\nc", "d", "\n") == ["ab\n", "cd\n"]
    @test buffer_lines("ab\ncd\ne", "f", "\n") == ["ab\n", "cd\n", "ef\n"]
    @test buffer_lines("a", "b\ncd\ne", "f", "\n") == ["ab\n", "cd\n", "ef\n"]
    @test buffer_lines("a", "b\ncd\n", "ef", "\n") == ["ab\n", "cd\n", "ef\n"]
end

@testset "LoggingStream" begin
    @test_logs (:info, Text("L1"),     Logging, Ignored(), :asdf) #=
            =# (:info, Text("L2"),     Logging, Ignored(), :asdf) #=
            =# (:info, Text("L3\nL4"), Logging, Ignored(), :asdf) #=
    =# begin
        ls = Logging2.LoggingStream(current_logger(), id=:asdf)
        write(ls, "L1")
        write(ls, "L2\n")
        write(ls, "L3\nL4") # No buffering!
    end

    # Check that log level early-out testing works
    @test_logs min_level=Logging.Warn begin
        ls = Logging2.LoggingStream(current_logger(), id=:asdf)
        write(ls, "L1")
    end
end

@testset "redirect_stdout and redirect_stderr" begin
    @test_logs (:info,Text("Hi"))  (:info,Text("Hi2")) #=
        =#     (:info,Text("Hi3")) (:info,Text("Hi4")) #=
        =# redirect_stdout(current_logger()) do
            println("Hi")
            println("Hi2\nHi3")
            print("Hi")
            print("4")
            print("\n")
        end

    # test return value from `redirect_stdout`
    @test (@test_logs (:info,Text("Hi")) redirect_stdout(current_logger()) do
        println("Hi")
        101
    end) == 101

    @test_logs (:info,Text("Hi")) #=
        =# redirect_stderr(current_logger()) do
            println(stderr, "Hi")
        end

    @test_logs (:warn,Text("Hi")) #=
        =# redirect_stderr(current_logger(); level=Logging.Warn) do
            println(stderr, "Hi")
        end
end


text:
import GPUArrays: allowscalar, @allowscalar

function _getindex(xs::CuArray{T}, i::Integer) where T
  buf = Array{T}(undef)
  copyto!(buf, 1, xs, i, 1)
  buf[]
end

function _setindex!(xs::CuArray{T}, v::T, i::Integer) where T
  copyto!(xs, i, T[v], 1, 1)
end

## logical indexing

Base.getindex(xs::CuArray, bools::AbstractArray{Bool}) = getindex(xs, CuArray(bools))

function Base.getindex(xs::CuArray{T}, bools::CuArray{Bool}) where {T}
  bools = reshape(bools, prod(size(bools)))
  indices = cumsum(bools)  # unique indices for elements that are true

  n = _getindex(indices, length(indices))  # number that are true
  ys = CuArray{T}(undef, n)

  if n > 0
    function kernel(ys::CuDeviceArray{T}, xs::CuDeviceArray{T}, bools, indices)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if i <= length(xs) && bools[i]
            b = indices[i]   # new position
            ys[b] = xs[i]
        end

        return
    end

    function configurator(kernel)
        config = launch_configuration(kernel.fun)

        threads = min(length(indices), config.threads)
        blocks = cld(length(indices), threads)

        return (threads=threads, blocks=blocks)
    end

    @cuda name="logical_getindex" config=configurator kernel(ys, xs, bools, indices)
  end

  unsafe_free!(indices)

  return ys
end

## find*

function Base.findall(bools::CuArray{Bool})
    I = if VERSION >= v"1.2"
        keytype(bools)
    elseif bools isa CuVector
        Int
    else
        CartesianIndex{ndims(bools)}
    end
    indices = cumsum(reshape(bools, prod(size(bools))))

    n = _getindex(indices, length(indices))
    ys = CuArray{I}(undef, n)

    if n > 0
        function kernel(ys::CuDeviceArray, bools, indices)
            i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

            @inbounds if i <= length(bools) && bools[i]
                i′ = CartesianIndices(bools)[i]
                b = indices[i]   # new position
                ys[b] = i′
            end

            return
        end

        function configurator(kernel)
            config = launch_configuration(kernel.fun)

            threads = min(length(indices), config.threads)
            blocks = cld(length(indices), threads)

            return (threads=threads, blocks=blocks)
        end

        @cuda name="findall" config=configurator kernel(ys, bools, indices)
    end

    unsafe_free!(indices)

    return ys
end

function Base.findall(f::Function, A::CuArray)
    bools = map(f, A)
    ys = findall(bools)
    unsafe_free!(bools)
    return ys
end

function Base.findfirst(testf::Function, xs::CuArray)
    I = if VERSION >= v"1.2"
        keytype(xs)
    else
        eltype(keys(xs))
    end

    y = CuArray([typemax(Int)])

    function kernel(y::CuDeviceArray, xs::CuDeviceArray)
        i = threadIdx().x + (blockIdx().x - 1) * blockDim().x

        @inbounds if i <= length(xs) && testf(xs[i])
            CUDAnative.@atomic y[1] = min(y[1], i)
        end

        return
    end

    function configurator(kernel)
        config = launch_configuration(kernel.fun)

        threads = min(length(xs), config.threads)
        blocks = cld(length(xs), threads)

        return (threads=threads, blocks=blocks)
    end

    @cuda name="findfirst" config=configurator kernel(y, xs)

    first_i = _getindex(y, 1)
    return keys(xs)[first_i]
end

Base.findfirst(xs::CuArray{Bool}) = findfirst(identity, xs)

function Base.findmin(a::CuArray; dims=:)
    if dims != Colon()
        error("Unsupported")
    end
    m = minimum(a)
    i = findfirst(x->x==m, a)
    return (m, i)
end

function Base.findmax(a::CuArray; dims=:)
    if dims != Colon()
        error("Unsupported")
    end
    m = maximum(a)
    i = findfirst(x->x==m, a)
    return (m, i)
end