text:
export set_relaxation, get_defaultparams

"""
    relax_ctx = set_relaxation(pb::Problem; ismultiordered=false, issparse=false, symmetries=Set(), hierarchykind=:Complex, renamevars=false, di=Dict{String, Int}(), d=-1)

    Build a `relax_ctx` object containing relaxation choices and problem features : order by constraint, relaxation order by constraint...
"""
function set_relaxation(pb::Problem; ismultiordered::Bool=false,
                                     issparse::Bool=false,
                                     symmetries::Array{DataType, 1}=DataType[],
                                     hierarchykind::Symbol=:Complex,
                                     renamevars::Bool=false,
                                     di::Dict{String, Int}=Dict{String, Int}(),
                                     d::Int=-1,
                                     params=Dict())

    # Check that all variables have a type fitting the hierarchy kind
    for (varname, vartype) in pb.variables
        (vartype<:Int) && error(LOGGER, "set_relaxation() : variable $varname,$vartype is integer, unfit for SOS relaxation.\nConsider relaxing it and adding a complementarity constraint.")
        (hierarchykind==:Complex) && !(vartype<:Complex) && error(LOGGER, "set_relaxation() : variable $varname,$vartype should be complex for complex hierarchy.")
        (hierarchykind==:Real) && !(vartype<:Real) && error(LOGGER, "set_relaxation() : variable $varname,$vartype should be real for real hierarchy.")
    end

    relax_ctx = RelaxationContext()

    relaxparams = relax_ctx.relaxparams

    # Compute each constraint degree
    relctx_setki!(relax_ctx, pb)

    # Store each SDP multiplier type
    relctx_setSDPmulttypes!(relax_ctx, pb, hierarchykind)

    # Relaxation order management
    relctx_setdi!(relax_ctx, pb, di, d)

    rel_ctx_setsymetries!(relax_ctx, pb, symmetries)

    log_POPcharact!(relax_ctx, pb)

    relaxparams[:opt_issparse] = issparse
    relaxparams[:opt_hierarchykind] = hierarchykind
    relaxparams[:opt_multiordered] = ismultiordered
    relaxparams[:opt_globalorder] = relax_ctx.di[get_momentcstrname()]
    (PhaseInvariance in symmetries) && (relaxparams[:opt_sym_phaseinv] = true)

    for (param, val) in params
        if haskey(relax_ctx.relaxparams, param)
            relax_ctx.relaxparams[param] = val
        else
            warn(LOGGER, "set_relaxation(): Unhandled parameter $param")
        end
    end

    init_output(relax_ctx)

    print_build_relctx(relax_ctx, pb)
    return relax_ctx
end

function relctx_setki!(relax_ctx, pb::Problem)
    ki = relax_ctx.ki
    ki[get_momentcstrname()] = 0
    for (cstrname, cstr) in pb.constraints
        cstrtype = get_cstrtype(cstr)
        if cstrtype == :ineqdouble
            cstrname_lo, cstrname_up = get_cstrname(cstrname, cstrtype)
            ki[cstrname_lo] = max(cstr.p.degree.explvar, cstr.p.degree.conjvar)
            ki[cstrname_up] = max(cstr.p.degree.explvar, cstr.p.degree.conjvar)
        else
            ki[get_cstrname(cstrname, cstrtype)] = max(cstr.p.degree.explvar, cstr.p.degree.conjvar)
        end
    end
    return
end

function relctx_setSDPmulttypes!(relax_ctx, pb::Problem, hierarchykind)
    ctrtypes = relax_ctx.cstrtypes
    for (cstrname, cstr) in pb.constraints
        cstrtype = get_cstrtype(cstr)
        if cstrtype == :ineqdouble
            cstrname_lo, cstrname_up = get_cstrname(cstrname, cstrtype)
            ctrtypes[cstrname_lo] = (hierarchykind==:Complex ? :SDPC : :SDP)
            ctrtypes[cstrname_up] = (hierarchykind==:Complex ? :SDPC : :SDP)
        elseif cstrtype == :eq
            ctrtypes[get_cstrname(cstrname, cstrtype)] = (hierarchykind==:Complex ? :SymC : :Sym)
        else
            ctrtypes[get_cstrname(cstrname, cstrtype)] = (hierarchykind==:Complex ? :SDPC : :SDP)
        end
    end
    ctrtypes[get_momentcstrname()] = (hierarchykind==:Complex ? :SDPC : :SDP)
    return
end

function relctx_setdi!(relax_ctx, pb::Problem, di, d)
    di_relax = relax_ctx.di
    !((di == Dict{String, Int}()) && (d==-1)) || error(LOGGER, "RelaxationContext(): Either di or d should be provided as input.")

    ## Checking order by constraint
    for (cstrname, cstr) in pb.constraints
        cur_order = haskey(di, cstrname) ? di[cstrname] : d

        # Check provided di is suitable wrt constraint degree, add
        cstrtype = get_cstrtype(cstr)
        if cstrtype == :ineqdouble
            cstrname_lo, cstrname_up = get_cstrname(cstrname, cstrtype)
            cur_ki = relax_ctx.ki[cstrname_lo]
            (0 ≤ cur_order-ceil(cur_ki/2)) || warn(LOGGER, "RelaxationContext(): Provided order ($cur_order) is lower than constraint $cstrname order ($cur_ki). \nUsing value ceil($cur_ki/2).")
            # (cur_ki <= cur_order) || warn(LOGGER, "RelaxationContext(): Provided order ($cur_order) is lower than constraint $cstrname order ($cur_ki). \nUsing value $cur_ki, hierarchy may be multiordered.")
            di_relax[cstrname_lo] = cur_order
            di_relax[cstrname_up] = cur_order
            # di_relax[cstrname_lo] = max(cur_order, ceil(cur_ki/2))
            # di_relax[cstrname_up] = max(cur_order, ceil(cur_ki/2))
        else # :eq, :ineqlo, :ineqhi
            cur_ki = relax_ctx.ki[get_cstrname(cstrname, cstrtype)]
            (0 ≤ cur_order-ceil(cur_ki/2)) || warn(LOGGER, "RelaxationContext(): Provided order ($cur_order) is lower than constraint $cstrname order ($cur_ki). \nUsing value ceil($cur_ki/2).")
            # (cur_ki <= cur_order) || warn(LOGGER, "RelaxationContext(): Provided order ($cur_order) is lower than constraint $cstrname order ($cur_ki). \nUsing value $cur_ki, hierarchy may be multiordered.")
            di_relax[get_cstrname(cstrname, cstrtype)] = cur_order
            # di_relax[get_cstrname(cstrname, cstrtype)] = max(cur_order, ceil(cur_ki/2))
        end
    end

    ## Setting order for moment contraint(s)
    if haskey(di, get_momentcstrname())
        di_relax[get_momentcstrname()] = di[get_momentcstrname()]
    elseif d!=-1
        di_relax[get_momentcstrname()] = d
    else
        di_relax[get_momentcstrname()] = maximum(values(di_relax))
    end

    # Checking order of objective
    obj_degree = max(pb.objective.degree.explvar, pb.objective.degree.conjvar)
    # if obj_degree > di_relax[get_momentcstrname()]
    #     warn(LOGGER, "RelaxationContext(): Moment matrix order $(di_relax[get_momentcstrname()]) is lower than objective degree ($obj_degree). \nUsing value $obj_degree, hierarchy may be multiordered.")
    #     di_relax[get_momentcstrname()] = ceil(obj_degree/2)
    # end

    return
end

function rel_ctx_setsymetries!(relax_ctx, pb, symmetries)
    pbsymmetries = relax_ctx.symmetries
    isa(symmetries, Array) || error(LOGGER, "set_relaxation(): symmetries should be an Array of types.")
    for symtype in symmetries
        if has_symmetry(relax_ctx, pb, symtype)
            push!(pbsymmetries, symtype)
        end
    end
    relax_ctx.symmetries = pbsymmetries

    if has_symmetry(relax_ctx, pb, PhaseInvariance)
        relax_ctx.relaxparams[:pb_isphaseinv] = true
    end
    return
end

function log_POPcharact!(relax_ctx::RelaxationContext, pb::Problem)
    params = relax_ctx.relaxparams

    params[:pb_nvar] = length(pb.variables)
    for (varname, vartype) in pb.variables
        if vartype <: Complex
            params[:pb_nvar_cplx] += 1
        elseif vartype <: Real
            params[:pb_nvar_real] += 1
        elseif vartype <: Bool
            params[:pb_nvar_bin] += 1
        end
    end

    params[:pb_nctr] = length(pb.constraints)
    params[:pb_nctr_eq] = -1
    params[:pb_nctr_ineqsimple] = -1
    params[:pb_nctr_ineqdouble] = -1

    params[:pb_maxpolydeg] = maximum(values(relax_ctx.ki))
end

function get_defaultparams()
    defparams = OrderedDict{Symbol, Any}(:pb_name=>"",          # POP parameters
                        :pb_nvar=>0,
                        :pb_nvar_cplx=>0,
                        :pb_nvar_real=>0,
                        :pb_nvar_bin=>0,
                        :pb_nctr=>0,
                        :pb_nctr_eq=>-1,
                        :pb_nctr_ineqsimple=>-1,
                        :pb_nctr_ineqdouble=>-1,
                        :pb_maxpolydeg=>-1,
                        :pb_isphaseinv=>false,
                        :slv_mmtrel_t=>-1.0,                    # Relaxation time and memory consumption
                        :slv_mmtrel_bytes=>-1,
                        :slv_sosrel_t=>-1.0,
                        :slv_sosrel_bytes=>-1,
                        :slv_mskstruct_t=>-1.0,
                        :slv_mskstruct_bytes=>-1,
                        :slv_prosta=>"",                        # SDP solution values
                        :slv_solsta=>"",
                        :slv_primobj=>-1.0,
                        :slv_dualobj=>-1.0,
                        :slv_primfeas=>-1.0,
                        :slv_dualfeas=>-1.0,
                        :slv_solvetime=>-1.0,
                        :opt_hierarchykind=>:Undef,
                        :opt_issparse=>false,                   # Relaxation parameters
                        :opt_multiordered=>false,
                        :opt_globalorder=>-1,
                        :opt_sym_phaseinv=>false,
                        :opt_nb_cliques=>-1,
                        :opt_msk_maxtime=>-1,                   # Default -1 is no time limit
                        :opt_outmode=>0,                        # 0: screen, 1: file, 2: both
                        :opt_outlev=>1,                         # 0: none, 1:summary at moment relaxation, sos relaxation, 2: detailled information, 3: full problems
                        :opt_outname=>"momentsos.log",
                        :opt_outcsv=>0,                         # 0: no csv is written, 1: csv is written
                        :opt_outcsvname=>"momentsos_solve.csv")

    return defparams
end


text:
using Test
using LinearAlgebra

include("../../lib/KTBC.jl")
using .KTBC: CreateKTBC

# 1
K₂ = CreateKTBC(2)[1]

Q = (1 / √2) * [
  1 1
  1 -1
]

Λ = [
  1 0
  0 3
]

@test K₂ == round.(Q * Λ * Q')

# 3, 4
K₅ = convert(Array{Int}, CreateKTBC(5)[1])
eK₅ = round.(eigvals(K₅), digits=3)

@test eK₅ == round.([
  2 - √3
  2 - 1
  2 - 0
  2 + 1
  2 + √3
], digits=3)

@test eK₅ == round.(2 * ones(5, 1) - 2 * cos.((1:5) * π / 6), digits=3)[1:5]

Q = round.(eigvecs(K₅), digits=3)
DST_C = round.(Q * diagm(0 => [1,-1,-1,-1,-1]), digits=3)

JK = (1:5) * (1:5)'

DST = sin.(JK * π / 6)/√3

@test DST_C == round.(DST, digits=3)
@test round.(DST', digits=3) == round.(DST^-1, digits=3)

# 6
T₆ = convert(Array{Int}, CreateKTBC(6)[2])
eT₆ = round.(eigvals(T₆), digits=3)

@test eT₆ == round.(2 * ones(6, 1) - 2 * cos.((.5:5.5) * π / 6.5), digits=3)[1:6]

# 7
C₄ = convert(Array{Int}, CreateKTBC(4)[4])
Λ, Q = eigen(C₄)

F = [
  1 1 1 1
  1 im im^2 im^3
  1 im^2 im^4 im^6
  1 im^3 im^6 im^9
]

A = Q \ F

# 16
A₁ = [
  .6 .4
  .4 .6
]

eigvals(A₁)

@test !all(1 .> eigvals(A₁))

# (A₂)ᵏ -> 0
A₂ = [
  .6 .9
  .1 .6
]

eigvals(A₂)

@test all(1 .> eigvals(A₂))

# 19
B = [
  3 1
  0 2
]

evalsB, evecB = eigen(B)

@test B^10 == evecB * (diagm(0 => evalsB)^10) * evecB^-1

ext: jl


text:
import Test
import Test: @test, @test_throws, @test_broken, @test_skip, @test_warn, @test_nowarn, @inferred

# Re-documenting the assertions from Test manually, because we need Documenter to pick them up.
# We could use their original docs instead of using links, but Documenter tries to
# run doctests on them, which fail.
# TODO: perhaps it is possible to create all these in a loop to avoid repetition.

@doc "See [`Test.@test`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@test)."
:(@test)

@doc "See [`Test.@test_throws`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@test_throws)."
:(@test_throws)

@doc "See [`Test.@test_broken`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@test_broken)."
:(@test_broken)

@doc "See [`Test.@test_skip`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@test_skip)."
:(@test_skip)

@doc "See [`Test.@test_warn`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@test_warn)."
:(@test_warn)

@doc "See [`Test.@test_nowarn`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@test_nowarn)."
:(@test_nowarn)

@doc "See [`Test.@inferred`](https://docs.julialang.org/en/latest/stdlib/Test/#Test.@inferred)."
:(@inferred)

struct ReturnValue{T} <: Test.Result
    value :: T
end

"""
    @test_result expr

Records a result from the test.
The result of `expr` will be displayed in the report by calling `string()` on it.
"""
macro test_result(expr)
    quote
        Test.record(Test.get_testset(), ReturnValue($(esc(expr))))
    end
end

struct FailExplanation <: Test.Result
    description :: String
end

"""
    @test_fail descr

Report a fail, providing an additional description (must be convertable to `String`).
The description will be displayed in the final report at the end of the test run.
"""
macro test_fail(descr)
    quote
        Test.record(Test.get_testset(), FailExplanation($(esc(descr))))
    end
end

struct CriticalFailException <: Exception
end

const _test_macro_symbols = Symbol.([
    "@test",
    "@test_fail",
    "@test_throws",
    "@test_broken",
    "@inferred",
    "@test_warn",
    "@test_nowarn",
    ])

"""
    @critical expr

Terminates the testcase on failure of an assertion `expr`.
`expr` must start from one of [`@test`](@ref), [`@test_fail`](@ref),
[`@test_throws`](@ref), [`@test_broken`](@ref), [`@inferred`](@ref),
[`@test_warn`](@ref), [`@test_nowarn`](@ref).
"""
macro critical(expr)
    if expr.head != :macrocall || !(expr.args[1] in _test_macro_symbols)
        error("@critical must precede an assertion macro that records a result")
    end

    quote
        $(esc(expr))
        if is_failed(Test.get_testset().results[end])
            throw(CriticalFailException())
        end
    end
end

struct TestcaseOutcome
    results :: Array{Test.Result, 1}
    elapsed_time :: Float64
    output :: String
end

is_failed(::Any) = false
is_failed(::Test.Fail) = true
is_failed(::Test.Error) = true
is_failed(::FailExplanation) = true
is_failed(outcome::TestcaseOutcome) = any(map(is_failed, outcome.results))

mutable struct JuteTestSet <: Test.AbstractTestSet
    results :: Array{Test.Result, 1}

    JuteTestSet(descr; results=[]) = new(results)
end

function Test.record(ts::JuteTestSet, res::Test.Result)
    push!(ts.results, res)
end

function Test.finish(ts::JuteTestSet) end

function run_testcase(tc::Testcase, args, capture_output::Bool=false)
    succeeded = true
    results = Test.Result[]

    elapsed_time, output = with_output_capture(!capture_output) do
        t = time_ns()
        Test.@testset JuteTestSet results=:($results) begin
            try
                Base.invokelatest(tc.func, args...)
            catch e
                if typeof(e) == CriticalFailException
                    # Used to terminate the testcase from inside of any number of nested calls
                else
                    rethrow(e)
                end
            end
        end
        (time_ns() - t) / 1e9
    end

    if length(results) == 0
        push!(results, Test.Pass(:test, nothing, nothing, nothing))
    end

    TestcaseOutcome(results, elapsed_time, output)
end


text:
# This file implements different ways to specify a degree

using Combinatorics
using JuLIP.Chemistry: atomic_number

import Base: ==

export SparseSHIP

abstract type BasisSpec{BO, NZ} end
abstract type AnalyticBasisSpec{BO, NZ} <: BasisSpec{BO, NZ} end

"""
`get_basisspec(<: BasisSpec) -> allKL, Nu`
"""
function get_basisspec end

VecOrTup = Union{AbstractVector, Tuple}

_convert_Zs(Zs::Symbol) = (atomic_number(Zs),)
_convert_Zs(Zs::Integer) = (Zs,)
_convert_Zs(Zs::VecOrTup) = _convert_Zs( tuple(Zs...) )
_convert_Zs(Zs::NTuple{NZ, Symbol}) where {NZ} = atomic_number.(Zs)
_convert_Zs(Zs::NTuple{NZ, <: Integer}) where {NZ} = Int16.(Zs)

z2i(spec::AnalyticBasisSpec, z::Integer) = spec.z2i[z]

i2z(spec::AnalyticBasisSpec, z::Integer) = spec.Zs[i]

nspecies(spec::AnalyticBasisSpec) = length(spec.Zs)

admissible(D::AnalyticBasisSpec, k, l) = deg(D, k, l) <= D.deg

"""
`SparseSHIP` : a general sparse-grid type degree definition,
```
deg({k}, {l}) = csp ∑ᵢ (kᵢ + wL * lᵢ)
              + chc ∏ᵢ max(ahc, bhc + kᵢ + wL * lᵢ)
```

Constructors:
```
SparseSHIP(bo, deg; kwargs...)
SparseSHIP(Zs, bo, deg; kwargs...)
```
where the keyword arguments are (with defaults)
* `wL = 1.5`
* `csp = 1.0`
* `chc = 0.0`
* `ahc = 0.0`
* `bhc = 0.0`
"""
struct SparseSHIP{BO, NZ} <: AnalyticBasisSpec{BO, NZ}
   deg::IntS
   wL::Float64
   csp::Float64
   chc::Float64
   ahc::Float64
   bhc::Float64
   # --------------------
   Zs::NTuple{NZ, Int16}
   valbo::Val{BO}
   z2i::Dict{Int16, Int16}
end

==(s1::SparseSHIP, s2::SparseSHIP) =
      all( getfield(s1, i) == getfield(s2, i)
           for i = 1:fieldcount(SparseSHIP) )

SparseSHIP(bo::Integer, deg::Integer; kwargs...) =
      SparseSHIP(:X, bo, deg; kwargs...)

function SparseSHIP(Zs, bo::Integer, deg::Integer;
                    wL = 1.5,
                    csp = 1.0,
                    chc = 0.0,
                    ahc = 0.0,
                    bhc = 0.0)
   @assert wL > 0
   @assert deg > 0
   @assert bo >= 0
   Zs = _convert_Zs(Zs)
   z2i = Dict([ Int16(z) => Int16(i) for (i, z) in enumerate(Zs) ]...)
   return SparseSHIP(IntS(deg), Float64(wL), Float64(csp),
                     Float64(chc), Float64(ahc), Float64(bhc),
                     Zs, Val(bo), z2i)
end

deg(D::SparseSHIP, k::Integer, l::Integer) =
      k + D.wL * l

deg(D::SparseSHIP, kk::VecOrTup, ll::VecOrTup) =
      D.csp * sum( deg(D, k, l) for (k, l) in zip(kk, ll) ) +
      D.chc * prod( max(D.ahc, D.bhc + deg(D, k, l)) for (k, l) in zip(kk, ll) )

# if one k is non-zero all other ks and ls are zero then we get
#    csp * k + chc * max(ahc, bhc)^{N-1} * max(ahc, bhc + k)
# take body-order = 1 (minimal) then we get just
#    csp * k + chc * max(ahc, bhc + k)
function maxK(D::SparseSHIP)
   allk = 0:ceil(Int, D.deg / D.csp)
   # degs = D.csp * allk + D.chc * max.(ahc, bhc .+ allk)
   degs = [ deg(D, (k,), (0,)) for k in allk ]
   admissible = findall(degs .<= D.deg)
   return maximum(allk[admissible])
end

# For a pure 2-body potential we don't need an angular component
maxL(D::SparseSHIP{1}, k::Integer =  0) = 0

# if just one k is non-zero (or 0), just one l non-zero and all other
# ks and ls are zero then we get
#    csp * (k+wL*l) + chc * max(ahc, bhc + k + wL*l)
function maxL(D::SparseSHIP, k::Integer = 0)
   alll = 0:ceil(Int, (D.deg / D.csp - k) / D.wL)
   degs = [ deg(D, (k,), (l,)) for l in alll ]
   admissible = findall(degs .<= D.deg)
   if isempty(admissible)
      return 0
   else
      return maximum(alll[admissible])
   end
end

Dict(D::SparseSHIP{BO}) where {BO} = Dict("__id__" => "SHIPs_SparseSHIP",
                                "deg" => D.deg,
                                "wL" => D.wL,
                                "csp" => D.csp,
                                "chc" => D.chc,
                                "ahc" => D.ahc,
                                "bhc" => D.bhc,
                                "Zs" => D.Zs,
                                "bo" => BO)

convert(::Val{:SHIPs_SparseSHIP}, D::Dict) =
      SparseSHIP(D["Zs"], D["bo"], D["deg"],
                 wL = D["wL"], csp = D["csp"],
                 chc = D["chc"], ahc = D["ahc"], bhc = D["bhc"] )

# ---------------------------------------------------------------
#    generating the basis specification
# ---------------------------------------------------------------

function generate_KL(spec::AnalyticBasisSpec, TI = IntS, TF=Float64)
   allKL = NamedTuple{(:k, :l), Tuple{TI,TI}}[]
   degs = TF[]
   # find k, l such that deg(spec, k, l) ≦ deg
   for k = 0:maxK(spec), l = 0:maxL(spec, k)
      @assert admissible(spec, k, l)
      push!(allKL, (k=k, l=l))
      push!(degs, deg(spec, k, l))
   end
   # sort allKL according to total degree
   I = sortperm(degs)
   return allKL[I], degs[I]
end

# NOTE: this is a very rudimentary generate_ZKL that could be
#       overloaded and so that different `allKL` collections are
#       created for each species.

function generate_ZKL(spec::AnalyticBasisSpec, TI = IntS, TF=Float64)
   allKL, degs = generate_KL(spec, TI, TF)
   allZKL = ntuple( _->copy(allKL), nspecies(spec) )
   return return allZKL
end

"""
returns kk, ll
where kk, ll is BO-tuples of k and l indices
"""
function _kl(ν::StaticVector{N}, KL) where {N}
   kk = SVector( ntuple(i -> KL[i][ν[i]].k, N) )
   ll = SVector( ntuple(i -> KL[i][ν[i]].l, N) )
   return kk, ll
end

function _kl(ν::StaticVector{N}, izz::StaticVector{N}, KLZ) where {N}
   kk = SVector( ntuple(α -> KLZ[izz[α]][ν[α]].k, N) )
   ll = SVector( ntuple(α -> KLZ[izz[α]][ν[α]].l, N) )
   return kk, ll
end

"""
create a vector of Nu arrays with the right type information
for each body-order
"""
function _init_NuZ(bo::Integer, nz::Integer, TI=IntS)
   Nu = Matrix{Vector}(undef, bo, nz)
   for n = 1:bo, iz = 1:nz
      Nu[n, iz] = NamedTuple{       (:izz,              :ν),
                               Tuple{SVector{n, Int16}, SVector{n, TI}} }[]
   end
   # convert into an SVector to make the length a type parameters
   return SMatrix{bo,nz,Vector}(Nu)
end

function generate_ZKL_tuples(spec::AnalyticBasisSpec{BO}, rotcoefs;
                            filter=true) where {BO}
   # all possible (k, l) pairs
   allZKL = generate_ZKL(spec)
   # separate arrays for all body-orders and species
   NuZ = _init_NuZ(BO, nspecies(spec))
   for N = 1:BO
      _generate_ZKL_tuples!(NuZ[N,1], spec, rotcoefs, allZKL, Val(N); filter=filter)
      for iz = 2:nspecies(spec)
         append!(NuZ[N, iz], NuZ[N, 1])
      end
   end
   return allZKL, NuZ
end

"""
takes two N-body basis functions each specified by a species vector `zzi` and
an index vector `νi` into the A-basis and decides whether they define the
same basis function. E.g.,
```
ν = SVector(1,2,2)
z1 = SVector(1,1,2)
z2 = SVector(1,2,1)
z3 = SVector(2,1,1)
@show _iseqB(z1, ν, z2, ν) # -> true
@show _iseqB(z1, ν, z3, ν) # -> false
```
"""
@generated function _iseqB(zz1::StaticVector{N}, ν1::StaticVector{N},
                           zz2::StaticVector{N}, ν2::StaticVector{N}
                           ) where {N}
   ex1 = "a1 = SVector(" * prod("(zz1[$i], ν1[$i])," for i = 1:N)[1:end-1] * ")"
   ex2 = "a2 = SVector(" * prod("(zz2[$i], ν2[$i])," for i = 1:N)[1:end-1] * ")"
   quote
      $(Meta.parse(ex1))
      $(Meta.parse(ex2))
      return sort(a1) == sort(a2)
   end
end

function _generate_ZKL_tuples!(NuZ, spec::AnalyticBasisSpec, rotcoefs, ZKL, ::Val{BO};
                               filter=true) where {BO}
   nz = nspecies(spec)
   izz = @MVector ones(Int16, BO)
   izz_tmp = SVector{BO, Int16}[]

   # temporary storage for all ν-tuples for a given zz combination
   Nu = SVector{BO, IntS}[]
   Nu_old = SVector{BO, IntS}[]

   for izz_ci in CartesianIndices(ntuple(_->(1:nz), BO))
      izz = Int16.(SVector(Tuple(izz_ci)...))
      # don't use this zz unless it is sorted; we will look at all permutations
      # of zz below; this means lots of skipping, but who cares, this part of
      # the loop is cheap. All the cost comes later.
      if !issorted(Tuple(izz))
         continue
      end

      # collect all possible ν tuples for the current species combination
      # specified by `zz`.  A (zz, ν) combination specifies a basis function
      #    ∏_a A[zₐ][νₐ]     (actually izₐ instead of zₐ and ignoring the m's)
      empty!(Nu)
      _generate_KL_tuples!(Nu, spec, rotcoefs, ZKL[izz]; filter=filter)

      # now loop through all the ν tuples we found to push them into NuZ
      for ν in Nu
         # for each such ν we also need to check whether any permutations
         # of `izz` give a new basis function as well. E.g.,
         #   A[1,n] A[1,m] == A[1,m] A[1,n]
         #   A[1,n] A[2,n] == A[2,n] A[1,n]
         # but
         #   A[1,n] A[2,m] != A[2,n] A[1,m]
         empty!(izz_tmp)
         push!(izz_tmp, izz)
         push!(NuZ, (izz = izz, ν = ν))
         for izzp in unique(permutations(izz))
            if !any( _iseqB(SVector(izzp...), ν, izz1, ν)  for izz1 in izz_tmp )
               push!(izz_tmp, izzp)
               push!(NuZ, (izz = izzp, ν = ν))
            end
         end
      end
   end
end

function _generate_KL_tuples!(Nu::Vector{<: SVector{BO}},
                              spec::AnalyticBasisSpec,
                              rotcoefs, ZKLs;
                              filter=true) where {BO}
   # the first item is just (1, ..., 1)
   # we can choose (k1, l1), (k2, l2) ... by indexing into allKL
   # then we start incrementing until we hit the maximum degree
   # while retaining the ordering ν₁ ≤ ν₂ ≤ …
   lastidx = 0
   ν = @MVector ones(IntS, BO)
   kk, ll = _kl(ν, ZKLs)
   while true
      # check whether the current ν tuple is admissible
      # the first condition is that its max index is small enough
      isadmissible = all(ν[i] <= length(ZKLs[i]) for i = 1:BO)
      if isadmissible
         # the second condition is that the multivariate degree it defines
         # is small enough => for that we first have to compute the corresponding
         # k and l vectors
         kk, ll = _kl(ν, ZKLs)
         isadmissible = admissible(spec, kk, ll)
      end

      # we want to increment `curindex`, but if we've reach the maximum degree
      # then we need to move to the next index down

      # if the current tuple ν has admissible degree ...
      if isadmissible
         # ... then we add it to the stack  ...
         #     (at least if it is an admissible basis function respecting
         #      all the symmetries - this is checked by filter_tuple)
         if !filter || filter_tuple(ll, rotcoefs)
            push!(Nu, SVector(ν))
         end
         # ... and increment it
         lastidx = BO
         ν[lastidx] += 1
      else
         # we have overshot, _deg(ν) > deg; we must go back down, by
         # decreasing the index at which we increment
         if lastidx == 1
            break
         end
         ν[lastidx-1:end] .= ν[lastidx-1] + 1
         lastidx -= 1
      end
   end
   return nothing
end

function filter_tuple(ll, rotcoefs)
   if isodd(sum(ll))
      return false
   end
   Bcoefs = SHIPs.Rotations.single_B(rotcoefs, ll)
   return norm(Bcoefs) > 1e-12
end

max_stars_repo_name: cortner/SHIPs

max_issues_repo_name: cortner/SHIPs

max_forks_repo_name: cortner/SHIPs


text:
using Base.Test
using Distributions
using ConjugatePriors

import ConjugatePriors: posterior, posterior_rand, posterior_mode, posterior_randmodel, fit_map

n = 100
w = rand(100)

# auxiliary tools

function ccount(K, x)
	r = zeros(K)
	for i = 1:length(x)
		r[x[i]] += 1.0
	end
	r
end

function ccount(K, x, w)
	r = zeros(K)
	for i = 1:length(x)
		r[x[i]] += w[i]
	end
	r
end

# Beta - Bernoulli

pri = Beta(1.0, 2.0)

x = rand(Bernoulli(0.3), n)
p = posterior(pri, Bernoulli, x)
@test isa(p, Beta)
@test_approx_eq p.α pri.α + sum(x)
@test_approx_eq p.β pri.β + (n - sum(x))

f = fit_map(pri, Bernoulli, x)
@test isa(f, Bernoulli)
@test_approx_eq succprob(f) mode(p)

p = posterior(pri, Bernoulli, x, w)
@test isa(p, Beta)
@test_approx_eq p.α pri.α + sum(x .* w)
@test_approx_eq p.β pri.β + (sum(w) - sum(x .* w))

f = fit_map(pri, Bernoulli, x, w)
@test isa(f, Bernoulli)
@test_approx_eq succprob(f) mode(p)

# posterior_rand & posterior_randmodel

pri = Beta(1.0, 2.0)
x = rand(Bernoulli(0.3), n)
post = posterior(pri, Bernoulli, x)

pv = posterior_rand(pri, Bernoulli, x)
@test isa(pv, Float64)
@test 0. <= pv <= 1.

pv = posterior_rand(pri, Bernoulli, x, w)
@test isa(pv, Float64)
@test 0. <= pv <= 1.

pm = posterior_randmodel(pri, Bernoulli, x)
@test isa(pm, Bernoulli)
@test 0. <= succprob(pm) <= 1.

pm = posterior_randmodel(pri, Bernoulli, x, w)
@test isa(pm, Bernoulli)
@test 0. <= succprob(pm) <= 1.

# Beta - Binomial

x = rand(Binomial(10, 0.3), n)
p = posterior(pri, Binomial, (10, x))
@test isa(p, Beta)
@test_approx_eq p.α pri.α + sum(x)
@test_approx_eq p.β pri.β + (10n - sum(x))

f = fit_map(pri, Binomial, (10, x))
@test isa(f, Binomial)
@test ntrials(f) == 10
@test_approx_eq succprob(f) mode(p)

p = posterior(pri, Binomial, (10, x), w)
@test isa(p, Beta)
@test_approx_eq p.α pri.α + sum(x .* w)
@test_approx_eq p.β pri.β + (10 * sum(w) - sum(x .* w))

f = fit_map(pri, Binomial, (10, x), w)
@test isa(f, Binomial)
@test ntrials(f) == 10
@test_approx_eq succprob(f) mode(p)

# Dirichlet - Categorical

pri = Dirichlet([1., 2., 3.])

x = rand(Categorical([0.2, 0.3, 0.5]), n)
p = posterior(pri, Categorical, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + ccount(3, x)

f = fit_map(pri, Categorical, x)
@test isa(f, Categorical)
@test_approx_eq probs(f) mode(p)

p = posterior(pri, Categorical, x, w)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + ccount(3, x, w)

f = fit_map(pri, Categorical, x, w)
@test isa(f, Categorical)
@test_approx_eq probs(f) mode(p)

# Dirichlet - Multinomial

x = rand(Multinomial(100, [0.2, 0.3, 0.5]), 1)
p = posterior(pri, Multinomial, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + x

r = posterior_mode(pri, Multinomial, x)
@test_approx_eq r mode(p)

x = rand(Multinomial(10, [0.2, 0.3, 0.5]), n)
p = posterior(pri, Multinomial, x)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + vec(sum(x, 2))

r = posterior_mode(pri, Multinomial, x)
@test_approx_eq r mode(p)

p = posterior(pri, Multinomial, x, w)
@test isa(p, Dirichlet)
@test_approx_eq p.alpha pri.alpha + vec(x * w)

r = posterior_mode(pri, Multinomial, x, w)
@test_approx_eq r mode(p)

# Gamma - Exponential

pri = Gamma(1.5, 2.0)

x = rand(Exponential(2.0), n)
p = posterior(pri, Exponential, x)
@test isa(p, Gamma)
@test_approx_eq shape(p) shape(pri) + n
@test_approx_eq rate(p) rate(pri) + sum(x)

f = fit_map(pri, Exponential, x)
@test isa(f, Exponential)
@test_approx_eq rate(f) mode(p)

p = posterior(pri, Exponential, x, w)
@test isa(p, Gamma)
@test_approx_eq shape(p) shape(pri) + sum(w)
@test_approx_eq rate(p) rate(pri) + sum(x .* w)

f = fit_map(pri, Exponential, x, w)
@test isa(f, Exponential)
@test_approx_eq rate(f) mode(p)


text:
# Getters and setters with microbiome-flavored names

# species -> features
const nfeatures = nspecies
const getfeature = getspecies
const featurenames = specnames
const featuretotals = speciestotals

# sites -> samples
const nsamples = nsites
const getsample = getsite
const samplenames = sitenames
const sampletotals = sitetotals


text:
module StridedTensorXD

import Base: parent, size, strides, tail, setindex
using Base: @propagate_inbounds, RangeIndex, Dims
using LinearAlgebra
using TupleLabXD
using TupleLabXD: StaticLength

export StridedView, @strided, @unsafe_strided, sreshape, sview

"""
    const StridedTensorXD.SliceIndex = Union{RangeIndex,Colon}

SliceIndex is the union name of standard `RangeIndex` and `Colon`, where
`RangeIndex = Union{Int, AbstractRange{Int}, AbstractUnitRange{Int}}`.
"""
const SliceIndex = Union{RangeIndex,Colon}

"""
    const StridedTensorXD._NTHREADS = Ref(1)

The global constant represents the number of threads in use, default to `Ref(1)`.
"""
const _NTHREADS = Ref(1)

"""
    StridedTensorXD.set_num_threads(n::Int)

Set the number of threads, i.e., set the value of the global constant `_NTHREADS`.
"""
function set_num_threads(n::Int)
    N = Threads.nthreads()
    if n > N
        n = N
        _set_num_threads_warn(n)
    end
    # Ref objects are dereferenced (loaded or stored) with [],
    # i.e., explicitly give a value for a bitstype reference.
    _NTHREADS[] = n
end

@noinline function _set_num_threads_warn(n)
    @warn "Maximal number of threads limited by number of Julia threads,
            setting number of threads equal to Threads.nthreads() = $n"
end

"""
    StridedTensorXD.get_num_threads() -> Int

Get the current number of threads in use now.
"""
get_num_threads() = _NTHREADS[]

"""
    StridedTensorXD.enable_threads()

Set the number of threads to the maximum value available.
"""
enable_threads() = set_num_threads(Threads.nthreads())

"""
    StridedTensorXD.disable_threads()

Set the number of threads to 1.
"""
disable_threads() = set_num_threads(1)

"""
    StridedTensorXD.__init__()

Initialize by setting the number of threads as `Threads.nthreads()`.
"""
function __init__()
    set_num_threads(Threads.nthreads())
end

"""
    const StridedTensorXD._use_threaded_mul = Ref(false)

The global constant represents whether we use multithreading for matrix multiplicatins,
default to `false`.
"""
const _use_threaded_mul = Ref(false)

"""
    StridedTensorXD.use_threaded_mul()

Check whether we are using multithreading for matrix multiplicatins.
"""
use_threaded_mul() = _use_threaded_mul[]

"""
    StridedTensorXD.enable_threaded_mul()

Enable the use of multithreading for matrix multiplicatins.
"""
function enable_threaded_mul()
    _use_threaded_mul[] = true
    return
end

"""
    StridedTensorXD.disable_threaded_mul()

Disable the use of multithreading for matrix multiplicatins.
"""
function disable_threaded_mul()
    _use_threaded_mul[] = false
    return
end

include("abstractstridedview.jl")
include("stridedview.jl")
include("unsafestridedview.jl")
include("mapreduce.jl")
include("broadcast.jl")
include("macros.jl")
end


text:
################################################################################
# Add Collision Avoidance
################################################################################

function add_collision_avoidance!(game_con::GameConstraintValues, i::Int, j::Int, radius::T) where {T}
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	px = probsize.px

	add_constraint!(game_con.state_conlist[i], CollisionConstraint(n,px[i],px[j],radius), 2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end

function add_collision_avoidance!(game_con::GameConstraintValues, radius::Vector{T}) where {T}
	probsize = game_con.probsize
	p = probsize.p
	@assert p == length(radius)
	for i = 1:p
		for j in setdiff(1:p,i)
			ri = radius[i]
			rj = radius[j]
			add_collision_avoidance!(game_con, i, j, ri+rj)
		end
	end
	return nothing
end

function add_collision_avoidance!(game_con::GameConstraintValues, radius::T) where {T}
	p = game_con.probsize.p
	add_collision_avoidance!(game_con, radius*ones(p))
	return nothing
end

################################################################################
# Add Spherical Collision Avoidance
################################################################################

function add_spherical_collision_avoidance!(game_con::GameConstraintValues, i::Int, j::Int, radius::T) where {T}
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	pz = probsize.pz

	pzi = pz[i][1:3] # assumes that the 3 first dimensions of the state of player i correspond to x,y,z.
	pzj = pz[j][1:3] # assumes that the 3 first dimensions of the state of player i correspond to x,y,z.
	add_constraint!(game_con.state_conlist[i], CollisionConstraint(n,pzi,pzj,radius), 2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end

function add_spherical_collision_avoidance!(game_con::GameConstraintValues, radius::Vector{T}) where {T}
	probsize = game_con.probsize
	p = probsize.p
	@assert p == length(radius)
	for i = 1:p
		for j in setdiff(1:p,i)
			ri = radius[i]
			rj = radius[j]
			add_spherical_collision_avoidance!(game_con, i, j, ri+rj)
		end
	end
	return nothing
end

function add_spherical_collision_avoidance!(game_con::GameConstraintValues, radius::T) where {T}
	p = game_con.probsize.p
	add_spherical_collision_avoidance!(game_con, radius*ones(p))
	return nothing
end

################################################################################
# Add State Bounds
################################################################################

function add_state_bound!(game_con::GameConstraintValues, i::Int, x_max::AbstractVector, x_min::AbstractVector)
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	add_constraint!(game_con.state_conlist[i], StateBoundConstraint(n,x_max=x_max,x_min=x_min), 2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end

################################################################################
# Add Control Bounds
################################################################################

function add_control_bound!(game_con::GameConstraintValues, u_max::AbstractVector, u_min::AbstractVector)
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	add_constraint!(game_con.control_conlist, ControlBoundConstraint(m,u_max=u_max,u_min=u_min), 1:N-1)
	con  = game_con.control_conlist.constraints[end]
	inds = game_con.control_conlist.inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.control_conval, conval)
	return nothing
end

################################################################################
# Add Circle Constraint
################################################################################

function add_circle_constraint!(game_con::GameConstraintValues, i::Int,
	xc::AbstractVector, yc::AbstractVector, radius::AbstractVector)
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	px = probsize.px

	add_constraint!(
		game_con.state_conlist[i],
		TrajectoryOptimization.CircleConstraint(n, xc, yc, radius, px[i][1], px[i][2]),
		2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end

function add_circle_constraint!(game_con::GameConstraintValues,
	xc::AbstractVector, yc::AbstractVector, radius::AbstractVector)
	p = game_con.probsize.p
	for i = 1:p
		add_circle_constraint!(game_con, i, xc, yc, radius)
	end
	return nothing
end

################################################################################
# Wall Constraint
################################################################################
abstract type AbstractWall end

mutable struct Wall <: AbstractWall
	p1::AbstractVector # initial point of the boundary
	p2::AbstractVector # final point of the boundary
	v::AbstractVector  # vector orthogonal to (p2 - p1) and indicating the forbiden halfspace
end

function add_wall_constraint!(game_con::GameConstraintValues, i::Int, walls::AbstractVector{Wall})
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	px = probsize.px

	n_wall = length(walls)
	T = eltype(walls[1].p1)
	x1 = SVector{n_wall,T}([wall.p1[1] for wall in walls])
	y1 = SVector{n_wall,T}([wall.p1[2] for wall in walls])
	x2 = SVector{n_wall,T}([wall.p2[1] for wall in walls])
	y2 = SVector{n_wall,T}([wall.p2[2] for wall in walls])
	xv = SVector{n_wall,T}([wall.v[1] for wall in walls])
	yv = SVector{n_wall,T}([wall.v[2] for wall in walls])

	add_constraint!(
		game_con.state_conlist[i],
		WallConstraint(n, x1, y1, x2, y2, xv, yv, px[i][1], px[i][2]),
		2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end

function add_wall_constraint!(game_con::GameConstraintValues, walls::AbstractVector{<:AbstractWall})
	p = game_con.probsize.p
	for i = 1:p
		add_wall_constraint!(game_con, i, walls)
	end
	return nothing
end

################################################################################
# Wall3D Constraint
################################################################################

mutable struct Wall3D <: AbstractWall
	p1::AbstractVector # initial point of the boundary
	p2::AbstractVector # middle point of the boundary
	p3::AbstractVector # final point of the boundary
	v::AbstractVector  # vector orthogonal to the plane defined by (p1,p2,p3) and indicating the forbiden halfspace
end

function add_wall_constraint!(game_con::GameConstraintValues, i::Int, walls::AbstractVector{Wall3D})
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	pz = probsize.pz

	n_wall = length(walls)
	T = eltype(walls[1].p1)
	x1 = SVector{n_wall,T}([wall.p1[1] for wall in walls])
	y1 = SVector{n_wall,T}([wall.p1[2] for wall in walls])
	z1 = SVector{n_wall,T}([wall.p1[3] for wall in walls])
	x2 = SVector{n_wall,T}([wall.p2[1] for wall in walls])
	y2 = SVector{n_wall,T}([wall.p2[2] for wall in walls])
	z2 = SVector{n_wall,T}([wall.p2[3] for wall in walls])
	x3 = SVector{n_wall,T}([wall.p3[1] for wall in walls])
	y3 = SVector{n_wall,T}([wall.p3[2] for wall in walls])
	z3 = SVector{n_wall,T}([wall.p3[3] for wall in walls])
	xv = SVector{n_wall,T}([wall.v[1] for wall in walls])
	yv = SVector{n_wall,T}([wall.v[2] for wall in walls])
	zv = SVector{n_wall,T}([wall.v[3] for wall in walls])

	add_constraint!(
		game_con.state_conlist[i],
		Wall3DConstraint(n,
			x1, y1, z1,
			x2, y2, z2,
			x3, y3, z3,
			xv, yv, zv, pz[i][1], pz[i][2], pz[i][3]),
		2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end
################################################################################
# CylinderWall Constraint
################################################################################

mutable struct CylinderWall <: AbstractWall
	p::AbstractVector # initial point of the boundary
	v::Symbol         # middle point of the boundary
	l::Real           # final point of the boundary
	r::Real           # vector orthogonal to the plane defined by (p1,p2,p3) and indicating the forbiden halfspace
end

function add_wall_constraint!(game_con::GameConstraintValues, i::Int, walls::AbstractVector{CylinderWall})
	probsize = game_con.probsize
	N = probsize.N
	n = probsize.n
	m = probsize.m
	p = probsize.p
	pz = probsize.pz

	n_wall = length(walls)
	T = eltype(walls[1].p)
	p1 = SVector{n_wall,T}([wall.p[1] for wall in walls])
	p2 = SVector{n_wall,T}([wall.p[2] for wall in walls])
	p3 = SVector{n_wall,T}([wall.p[3] for wall in walls])
	v = SVector{n_wall,Symbol}([wall.v for wall in walls])
	l = SVector{n_wall,T}([wall.l for wall in walls])
	r = SVector{n_wall,T}([wall.r for wall in walls])

	add_constraint!(
		game_con.state_conlist[i],
		CylinderConstraint(n, p1, p2, p3, v, l, r, pz[i][1], pz[i][2], pz[i][3]),
		2:N)
	con  = game_con.state_conlist[i].constraints[end]
	inds = game_con.state_conlist[i].inds[end]
	conval = Altro.ALConVal(n,m,con,inds)
	push!(game_con.state_conval[i], conval)
	return nothing
end

################################################################################
# Helpers
################################################################################
function TrajectoryOptimization.cost_expansion!(conval::TrajectoryOptimization.AbstractConstraintValues)
    s = TrajectoryOptimization.sense(conval)
    for (i,k) in enumerate(conval.inds)
        TrajectoryOptimization.cost_expansion!(s, conval, i)
    end
    return nothing
end

function reset!(game_con::GameConstraintValues)
    reset_duals!(game_con)
    reset_penalties!(game_con)
    return nothing
end

function reset_duals!(game_con::GameConstraintValues)
    p = game_con.probsize.p
    for i = 1:p
        for conval in game_con.state_conval[i]
            Altro.reset_duals!(conval)
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        Altro.reset_duals!(conval)
    end
    return nothing
end

function reset_penalties!(game_con::GameConstraintValues)
    p = game_con.probsize.p
    for i = 1:p
        for conval in game_con.state_conval[i]
            Altro.reset_penalties!(conval)
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        Altro.reset_penalties!(conval)
    end
    return nothing
end

function penalty_update!(game_con::GameConstraintValues)
    p = game_con.probsize.p
    for i = 1:p
        for conval in game_con.state_conval[i]
            Altro.penalty_update!(conval)
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        Altro.penalty_update!(conval)
    end
    return nothing
end

function penalty_update!(con::Altro.ALConVal)
    for i in eachindex(con.μ)
        con.μ[i] .*= con.params.ϕ
    end
end

function dual_update!(game_con::GameConstraintValues)
    p = game_con.probsize.p
	α_dual = game_con.α_dual
	αx_dual = game_con.αx_dual
    for i = 1:p
        for conval in game_con.state_conval[i]
            # Altro.dual_update!(conval)
			dual_update!(conval, αx_dual[i])
        end
    end
    # Control constraints
    for conval in game_con.control_conval
		# Altro.dual_update!(conval)
		dual_update!(conval, α_dual)
    end
    return nothing
end

function evaluate!(game_con::GameConstraintValues, traj::Traj)
    p = game_con.probsize.p
    for i = 1:p
        for conval in game_con.state_conval[i]
            TrajectoryOptimization.evaluate!(conval, traj)
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        TrajectoryOptimization.evaluate!(conval, traj)
    end
    return nothing
end

import RobotDynamics.jacobian!
function jacobian!(game_con::GameConstraintValues, traj::Traj)
    p = game_con.probsize.p
    for i = 1:p
        for conval in game_con.state_conval[i]
            TrajectoryOptimization.jacobian!(conval, traj)
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        TrajectoryOptimization.jacobian!(conval, traj)
    end
    return nothing
end

function update_active_set!(game_con::GameConstraintValues, traj::Traj)
	evaluate!(game_con, traj)
	update_active_set!(game_con)
    return nothing
end

function update_active_set!(game_con::GameConstraintValues)
    p = game_con.probsize.p
    for i = 1:p
        for conval in game_con.state_conval[i]
            Altro.update_active_set!(conval, Altro.Val(game_con.active_set_tolerance))
        end
    end
    # Control constraints
    for conval in game_con.control_conval
        Altro.update_active_set!(conval, Altro.Val(game_con.active_set_tolerance))
    end
    return nothing
end

################################################################################
# Dual Update
################################################################################

function dual_update!(conval::Altro.ALConVal, α_dual::T) where {T}
	c = conval.vals
	λ = conval.λ
	μ = conval.μ
	λ_max = conval.params.λ_max
	cone = TrajectoryOptimization.sense(conval.con)
	for i in eachindex(conval.inds)
		λ[i] .= dual_update(cone, SVector(λ[i]), SVector(c[i]), SVector(μ[i]), λ_max, α_dual)
	end
end

function dual_update(::Equality, λ, c, μ, λmax, α_dual)
	λbar = λ + α_dual*μ .* c
	return clamp.(λbar, -λmax, λmax)
end

function dual_update(::Inequality, λ, c, μ, λmax, α_dual)
 	λbar = λ + α_dual*μ .* c
	return clamp.(λbar, 0, λmax)  # project onto the dual cone via max(0,x)
end

function dual_update(cone::Altro.SecondOrderCone, λ, c, μ, λmax, α_dual)
	 λbar = λ - α_dual*μ .* c
	 return TrajectoryOptimization.projection(cone, λbar)  # project onto the dual cone
end


text:
function plotlybase_syncplot(plt::Plot)
    plt.o = ORCA.PlotlyBase.Plot()
    traces = ORCA.PlotlyBase.GenericTrace[]
    for series_dict in plotly_series(plt)
        plotly_type = pop!(series_dict, :type)
        push!(traces, ORCA.PlotlyBase.GenericTrace(plotly_type; series_dict...))
    end
    ORCA.PlotlyBase.addtraces!(plt.o, traces...)
    layout = plotly_layout(plt)
    w, h = plt[:size]
    ORCA.PlotlyBase.relayout!(plt.o, layout, width = w, height = h)
    return plt.o
end

_show(io::IO, ::MIME{Symbol("image/png")}, plt::Plot{PlotlyBackend}) = ORCA.PlotlyBase.savefig(io, plotlybase_syncplot(plt), format = "png")


text:
"""
L = lagrange_basis_poly(x0,x,N,j)
--------------------------------------------------------------------------------------
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/26/2016, Last Modified: 1/25/2016
Citations: This function was influenced by the lagrange() function [located here](https://github.com/pjabardo/Jacobi.jl/blob/master/src/gauss_quad.jl)
--------------------------------------------------------------------------------------
# Input Arguments
* `x0`: point to approximate function value at
* `x`: x data to used calculate basis polynomials
* `N`: order of Lagrange interpolating polynomial
* `j`: index of interest

# Output Arguments
* `L`: Lagrange basis polynomials

A basic description of Lagrange interpolating polynomials is provided [here](http://127.0.0.1:8000/lagrange_poly.html#lagrange-poly)

"""
function lagrange_basis_poly(x,x_data,Nc,j)
    L = 1;
    for idx in 1:Nc+1 # use all of the data
      if idx!=j
        L = L*(x - x_data[idx])/(x_data[j]-x_data[idx]);
      end
    end
  return L
end
"""
y = interpolate_lagrange(x::Vector{T}, x_data::Vector{T}, y_data::Vector{T}) where {T <: Number}
"""
function interpolate_lagrange(x, x_data, y_data)

    Nc = length(x_data) - 1

    if length(x_data)!=length(y_data)
        error(string("\n",
                      "-------------------------------------------------------", "\n",
                      "There is an error with the data vector lengths!!", "\n",
                      "-------------------------------------------------------", "\n",
                      "The following variables should be equal:", "\n",
                      "length(x_data) = ",length(x_data),"\n",
                      "length(y_data) = ",length(y_data),"\n"
                      )
              )
      end
      ns = length(x)
      L = zeros(Float64,Nc+1,ns)
      x = x[:]; x_data = x_data[:]; y_data = y_data[:]; # make sure data is in a column
      for idx in 1:Nc+1
        for j in 1:ns
            L[idx,j] = lagrange_basis_poly(x[j],x_data,Nc,idx);
        end
      end
      return y_data'*L
end

"""
scale_w(w,ta,tb)
--------------------------------------------------------------------------------------
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/23/2017, Last Modified: 9/18/2017
--------------------------------------------------------------------------------------
"""
function scale_w(w::Vector{T},ta::T,tb::T)::Vector{T} where { T <: Number }
    return w .* (tb - ta) ./ 2
end

"""
scale_tau(tau,ta,tb)
--------------------------------------------------------------------------------------
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/23/2017, Last Modified: 12/06/2019
--------------------------------------------------------------------------------------
"""
function scale_tau(tau::Vector{T}, ta::T, tb::T)::Vector{T} where { T <: Number }
    return ( tau .* (tb - ta) .+ (ta + tb) ) ./ 2
end

"""
L = lagrange_basis_poly!(x,x_data,L)
--------------------------------------------------------------------------------------
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 12/26/2016, Last Modified: 12/02/2019
Last Modified By: Alexander Buck
Citations: This function was influenced by the lagrange() function [located here](https://github.com/pjabardo/Jacobi.jl/blob/master/src/gauss_quad.jl)
--------------------------------------------------------------------------------------
# Input Arguments
* `x`: point to approximate function value at
* `x_data`: x data to used calculate basis polynomials
* `L`: array of interest

# Output Arguments
* `L`: Lagrange basis polynomials

A basic description of Lagrange interpolating polynomials is provided [here](http://127.0.0.1:8000/lagrange_poly.html#lagrange-poly)

# TODO get ride of Nc and replace it with Nck[int]
"""
function lagrange_basis_poly!(x::AbstractArray{T},x_data,L::AbstractArray{T}) where {T<:Number}
    Nc = length(x_data) - 1
    ns = length(x);
    L = zeros(Float64,Nc+1,ns);
    for idx in 1:Nc+1
      for j in 1:ns
        L[idx,j] = lagrange_basis_poly(x[j],x_data,Nc,idx);
      end
    end
    return L
end
lagrange_basis_poly(x::AbstractArray{T},x_data,Nc) where {T<:Number} = lagrange_basis_poly!(x::AbstractArray{T},x_data,Nc,zeros(x))

"""
D = polyDiff(x);
--------------------------------------------------------------------------
Last modifed for julia on 12/06/2019 by Mattias Fält
Original Author: JJ.A.C. Weideman, S.C. Reddy 1998
Original Function Name: poldif.m  |  Source: [matlabcentral](https://www.mathworks.com/matlabcentral/fileexchange/29-dmsuite)
https://pdfs.semanticscholar.org/bae2/1eb9458f194887bc8d7808383f56d7f4dca0.pdf
https://math.stackexchange.com/questions/1105160/evaluate-derivative-of-lagrange-polynomials-at-construction-points
--------------------------------------------------------------------------
# Input Arguments
* `x::Vector`: Vector of N distinct nodes.

# Output Arguments
* `D::Array{Float64,2}`: Differention Matrix

"""
function polyDiff(x)  #TODO get rid of B stuff

    N = length(x)

    x = x[:]

    M = 1

    alpha = ones(N, 1)

    B = zeros(M, N)

    L = Matrix{Bool}(LinearAlgebra.I, N, N)

    XX = repeat(x, 1, N)

    # DX contains entries x(k)-x(j).
    DX = XX - XX'

    # Put 1's one the main diagonal.
    DX[L] = ones(N, 1)

    # Quantities c(j)
    c = alpha .* prod(DX, dims=2)

    C = repeat(c, 1, N)

    # Matrix with entries c(k)/c(j).
    C = C ./ C'

    # Z contains entries 1/(x(k)-x(j))
    Z = 1 ./ DX

    # with zeros on the diagonal.
    Z[L] = zeros(N, 1)

    # X is same as Z', but with
    X = Z'

    # diagonal entries removed.
    flag = trues(size(X));
    flag[L] .= false;
    X = X[flag];

    X = reshape(X, N-1, N)

    # Initialize Y and D matrices.
    Y = ones(N-1, N)

    # Y is matrix of cumulative sums,
    D = Matrix{Float64}(LinearAlgebra.I, N, N)

    temp = reshape(B[1,:], 1, N)

    # Diagonals
    Y = cumsum([ temp ; Y[1:N-1,:] .* X ], dims=1)

    # Off-diagonals
    D = Z .* (C .* repeat(LinearAlgebra.diag(D),1,N) - D)

    # Correct the diagonal
    D[L] = Y[N, :]

    return D

end

"""
# t must be the time vector
# V is any vector that you are interpolating
--------------------------------------------------------------------------------------\n
Author: Huckleberry Febbo, Graduate Student, University of Michigan
Date Create: 2/06/2017, Last Modified: 12/06/2019 \n
--------------------------------------------------------------------------------------\n
"""
function linearSpline(t::Vector, V::Vector)
    println("running linearSpline")
    if !isequal(length(t),length(V))
        error("!isequal(length(t),length(V))")
    end
    # remove any repeating values
    M = Array{Bool}(undef, length(t)); M[1:length(t)] .= false;
    for i in 1:length(t)-1
        if t[i]==t[i+1]
            M[i]=true
        else
            M[i]=false
        end
    end
    rm_idx = findall(M)
    if (length(t)==length(rm_idx))
        error("No time has elapsed and there will be an issue with interpolation. \n
            Cannot simulate the vehicle.")
    end

    # initialize vetors
    t_new = Array{Float64}(undef, length(t)-length(rm_idx))
    V_new = Array{Float64}(undef, length(t)-length(rm_idx))
    q = 1
    for i in 1:length(V) #TODO put an error message here if V and t are different sizes
        if !M[i]
            t_new[q] = t[i]
            V_new[q] = V[i]
            q = q + 1
        end
    end

    # make interpolant using Dierckx.jl
    #Spline1D(t_new,V_new,k=1)    # linear spline

    # make interpolant using Interpolations.jl
    knots = (t_new,)
    it = interpolate(knots,V_new,Gridded(Linear()))
    et = extrapolate(it, Flat())

    return et
end


text:
using Revise, MCMCChains, LinearAlgebra, UnPack

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

#holds the main state vector
#unfinished and untested
abstract type AbstractCCState{T} <: AbstractVector{T} end
struct CCState{TΠ<:AbstractVector, Tparams<:Dict, Tparamsidx, T} <: AbstractCCState{T}
  Π::TΠ
  params::Tparams
  paramsidx::Tparamsidx

  CCState(Π::TΠ, params::Tparams, paramsidx::Tparamsidx) where {
    TΠ, Tparams, Tparamsidx}=new{TΠ, Tparams, Tparamsidx, eltype(Π)}(Π, params, paramsidx)
end

#forms a state vector given a parameter index
#main constructor, but internally facing only
#the paramsidx object should map all params to indices in a vector
function CCState(paramsidx::Dict, ::Type{T}) where T
  Nparams = sum(length.(values(paramsidx)))
  Π = rand(T, Nparams)

  #build high-performance array wraps for indexing into the state vector
  params::Dict{Symbol, Union{Ref,AbstractArray{T}}} = Dict()
  for (n,p) ∈ pairs(paramsidx)
    Tp = typeof(p)
    if Tp <: Int
      params[n] = Ref(Π, p)
    elseif Tp <: AbstractVector
      params[n] = less_unsafe_wrap(Vector{T}, source=Π, start=minimum(p), dims=size(p))
    elseif Tp <: AbstractMatrix
      params[n] = less_unsafe_wrap(Matrix{T}, source=Π, start=minimum(p), dims=size(p))
    else
      @assert false
    end
  end

  #integrity checks
  length(Π)==sum(length.(values(params)))
  coversindex = falses(length(Π)) #keep track of which indices params references
  for (n,p) ∈ pairs(paramsidx)
    @assert size(p) ≡ size(params[n]) #check dimensions
    indicesused = [p...;]
    @assert all(reshape(Π[indicesused], size(params[n])) .≡ params[n]) #checks for substantial identity
    @assert all((!).(coversindex[indicesused])) #makes sure we have no overlaps
    coversindex[indicesused] .= true
  end
  @assert all(coversindex) #makes sure params references all indices of Π

  return CCState(Π,params,paramsidx)
end

###The below provides 1d indexing of the parameters
Base.length(Ψ::AbstractCCState) = length(Ψ.Π)
Base.size(Ψ::AbstractCCState) = (Base.length(Ψ),)
Base.IndexStyle(::Type{<:AbstractCCState}) = IndexLinear()

Base.getindex(Ψ::AbstractCCState, inds) = Ψ.Π[inds]
Base.setindex!(Ψ::AbstractCCState, x::Real, i::Int) = (Ψ.Π[i] = x)

#shortcut for directly indexing into the capacity state
function Base.getproperty(Ψ::TΨ, Fparam::Symbol) where TΨ <: AbstractCCState
  hasfield(TΨ, Fparam) && return getfield(Ψ, Fparam)
  Ψ.params[Fparam]
end

function Base.show(io::IO, ::MIME"text/plain", Ψ::AbstractCCState{T}) where T
  for k ∈ keys(Ψ.params)
    println(io, "Showing contents of $k: ")
    printmln(Ψ.params[k])
  end
end

function Base.deepcopy(Ψ::AbstractCCState{T}) where T
  Ψcpy = CCState(Ψ.paramsidx, T)

  #integrity checks
  Ψcpy .= Ψ
  @assert all(Ψcpy.Π .≡ Ψ.Π)
  for s ∈ keys(Ψ.params)
    @assert all(Ψcpy.params[s] .≡ Ψ.params[s])
  end

  return Ψcpy
end

#contains the MCMC chain
mutable struct CapacityChains{
    Tchain<:Chains, Tvalue, Tparams<:Dict, Tparamsidx<:Dict,
      Thyper<:Dict, Tidxiter, T}

  chain::Tchain #holds an MCMC chain object
  value::Tvalue
  params::Tparams
  paramsidx::Tparamsidx
  hyper::Thyper

  i::Int #the iteration counter
  previdx::Int
  idxiter::Tidxiter #index iterator, relevant during burn-in

  Nburn::Int
  burncomplete::Bool
  function CapacityChains(;
      chain::Tchain,
      params::Tparams,
      paramsidx::Tparamsidx,
      hyper::Thyper,
      Nburn::Int) where {Tchain, Tparams,Tparamsidx,Thyper}

    value = chain.value
    previdx=0
    idxiter=Iterators.cycle(1:(size(chain.value,1))) |> Iterators.Stateful

    burncomplete=Nburn≤0

    return new{Tchain, typeof(value), Tparams, Tparamsidx, Thyper, typeof(idxiter), eltype(chain.value)}(
      chain, value, params, paramsidx, hyper, 1, previdx, idxiter, Nburn, burncomplete)
  end
end

#main constructor for the MCMC chain
function CapacityChains(::Type{T};
  Nburn::Int = error("Nburn is required"),
  Ndraws::Int = error("Ndraws is required"),
  paramstemplate::Union{NamedTuple,Dict} = error("paramstemplate is required"),
  hyper::Union{Dict,NamedTuple}) where T <: Real

  #total number of parameters in the chain
  Nparams::Int = sum([length(p) for p ∈ values(paramstemplate)])
  paramnames::Vector{String} = Vector{String}(undef, Nparams)
  col::Int = 1 #the column in the values matrix
  paramsidx::Dict{Symbol, Union{Int, Vector{Int}, Matrix{Int}}} = Dict()
  for (i,(n,p)) ∈ enumerate(pairs(paramstemplate))
    nparams = length(p)
    @assert nparams == prod(size(p))

    if typeof(p) <: Real #the parameter is a single value
      paramnames[col] = n|>string
      paramsidx[n] = col
      col += 1

    #a corresponding  vector of column indices for a parameter vector
    elseif typeof(p) <: AbstractVector
      paramsidx[n] = similar(p, Int)
      colend = col+length(p)-1 #the end of the range
      paramsidx[n] .= col:colend

      #assign a string name to each parameter
      paramnames[col:colend] .= (coord->"$n[$coord]").(col:colend)

      col = col + length(p)
      @assert col == colend + 1
    elseif typeof(p) <: AbstractMatrix
      paramsidx[n] = similar(p, Int)
      colend = col+length(p)-1 #the end of the range
      paramsidx[n] .= reshape(col:colend, size(p)) #creates the indices

      #assign a string name to each parameter. The string includes name and index e.g. a[2,3]
      paramnames[col:colend] .= (coord->"$n[$(coord[1]),$(coord[2])]").(
        Iterators.product(1:size(p,1), 1:size(p,2)) |> collect |> vec)
      col = col + length(p)
      @assert col == colend + 1
    else
      error("Unrecognized type of parameter $n: $(typeof(p))")
    end
  end
  @assert Nparams == col - 1

  valuetemplate = rand(T, Ndraws, Nparams, 1)
  chain = Chains(valuetemplate, paramnames)

  #the below views are carefully set to allow linear indexing
  params = Dict(n => view(chain.value,:,vec([p...;]),1) for (n,p) ∈ pairs(paramsidx))

  cleanedhyper::Dict = typeof(hyper) <: Dict ? hyper : (hyper |> Dict)

  return CapacityChains(hyper=cleanedhyper; chain, params, paramsidx, Nburn)
end

#gets the index within the value array
function Base.getproperty(cc::CapacityChains, sym::Symbol)
  (sym === :idx) && return peek(cc.idxiter)
  return getfield(cc, sym)
end

#increments the state of the chain
function increment!(cc::CapacityChains)
  if cc.i == cc.Nburn #if we are at the end of the burning period, reset
    @assert !cc.burncomplete
    cc.burncomplete=true
    resetidx!(cc)
  else
    cc.previdx = cc.idx
    popfirst!(cc.idxiter)
    (cc.burncomplete) && (cc.idx < cc.previdx) && error("Attempted to increment chain past chain length")
  end

  cc.i+=1
  return nothing
end

#updates the chain based on the state
function updateandincrement!(cc::CapacityChains, Ψ::AbstractCCState; verify::Bool=true)

  cc.value[cc.idx, :, 1] .= Ψ.Π
  increment!(cc)

  @unpack previdx = cc #check the integrity of the assignment
  if verify
    for n ∈ keys(cc.params)
      (all(cc.params[n][previdx, :, 1] .== [getproperty(Ψ, n)...;])) || error(
        "!all(cc.params[n][previdx, :, 1] .== [getproperty(Ψ, n)...;])
        n=$n, previdx=$previdx
        cc.params[n][previdx, :, 1]=$(cc.params[n][previdx, :, 1])
        [getproperty(Ψ, n)...;]=$([getproperty(Ψ, n)...;])
        "
      )
    end
  end

  return nothing
end

#a couple of convenience methods
Base.length(cc::CapacityChains) = size(cc.value,1)
Base.size(cc::CapacityChains) = size(cc.value)[1:2]
Base.size(cc::CapacityChains, i::Int) = size(cc.value)[i]

#resets the circular indexing to put values in chronological order
function resetidx!(cc::CapacityChains)
  @unpack value = cc
  value .= @view value[[(cc.idx+1):end...;1:cc.idx...;], :, :]

  #reset the index
  while cc.idx ≠ 1
    popfirst!(cc.idxiter)
  end
  cc.previdx = length(cc)

  return nothing
end

CCState(cc::CapacityChains, ::Type{T}=eltype(cc.value)) where T = CCState(cc.paramsidx, T)

#Base.Broadcast.broadcastable(yp::YearPeriod) = Ref(yp)
#helper function that assigns an initial value
function CCState(cc::CapacityChains, initial::Union{NamedTuple,Dict},
  ::Type{T}=eltype(cc.value)) where T

  Ψ = CCState(cc, T)

  for (n,p) ∈ pairs(initial)
    Ψ.params[n] .= initial[n]
  end

  #integrity check
  for (n,p) ∈ pairs(Ψ.params)
    @assert all(Ψ.params[n] .== initial[n])
  end

  return Ψ
end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
# Chapter 1.1

# Visualizing a geometric series

## Julia code to generate a geometric series
using Plots:bar
p = 0.5
n = 1:10
X = p .^ n
bar(n, X, legend=false)

## compute N choose K
# Julia code to (N choose K) and K!
n = 10
k = 2
binomial(n, k)
factorial(k)

# Chapter 1.4

# Inner product of two vectors

## Julia code to perform and inner product
x = [1, 0, -1]
y = [3, 2, 0]
x'y

## Norm of a vector
# Julia code to compute the norm
using LinearAlgebra:norm
x = [1, 0, -1]
norm(x)

## Weighted norm of a vector
# Julia code to compute the weighted norm
W = [1. 2. 3.; 4. 5. 6.; 7. 8. 9.]
x = [2, -1, 1]
z = x'W*x

## System of linear equations
# Julia code to solve Xβ = y
X = [1. 3.; -2. 7.; 0. 1.]
y = [2., 1., 0.]
β = X\y


text:
export run_tfopf, run_ac_tfopf, run_dc_tfopf

""
function run_ac_tfopf(file, solver; kwargs...)
    return run_tfopf(file, ACPPowerModel, solver; kwargs...)
end

""
function run_dc_tfopf(file, solver; kwargs...)
    return run_tfopf(file, DCPPowerModel, solver; kwargs...)
end

""
function run_tfopf(file, model_constructor, solver; kwargs...)
    return run_generic_model(file, model_constructor, solver, post_tfopf; solution_builder = get_solution_tf, kwargs...)
end

""
function post_tfopf(pm::GenericPowerModel)
    PowerModels.variable_voltage(pm)
    PowerModels.variable_generation(pm)
    PowerModels.variable_branch_flow(pm) #d
    PowerModels.variable_dcline_flow(pm)
    variable_transformation(pm)

    PowerModels.objective_min_fuel_cost(pm)

    PowerModels.constraint_voltage(pm)
    for i in PowerModels.ids(pm, :ref_buses)
        PowerModels.constraint_theta_ref(pm, i)
    end

    for i in PowerModels.ids(pm, :bus)
        PowerModels.constraint_kcl_shunt(pm, i)
    end

    for i in PowerModels.ids(pm, :branch)
        branch = PowerModels.ref(pm, :branch, i)
        if branch["shiftable"] == false && branch["tappable"] == false
            PowerModels.constraint_ohms_yt_from(pm, i)
            PowerModels.constraint_ohms_yt_to(pm, i)
            constraint_link_voltage_magnitudes(pm, i)
        else
            constraint_variable_transformer_y_from(pm, i)
            constraint_variable_transformer_y_to(pm, i)
        end
        PowerModels.constraint_voltage_angle_difference(pm, i)

        PowerModels.constraint_thermal_limit_from(pm, i)
        PowerModels.constraint_thermal_limit_to(pm, i)
    end
    for i in PowerModels.ids(pm, :dcline)
        PowerModels.constraint_dcline(pm, i)
    end
end

max_issues_repo_name: frederikgeth/PowerModelsReliability

max_forks_repo_name: frederikgeth/PowerModelsReliability


text:
using NeuralNetDiffEq
using Flux
using DifferentialEquations
using LinearAlgebra
d = 100 # number of dimensions
X0 = fill(0.0f0, d) # initial value of stochastic control process
tspan = (0.0f0, 1.0f0)
λ = 1.0f0

g(X) = log(0.5f0 + 0.5f0 * sum(X.^2))
f(X,u,σᵀ∇u,p,t) = -λ * sum(σᵀ∇u.^2)
μ_f(X,p,t) = zero(X)  # Vector d x 1 λ
σ_f(X,p,t) = Diagonal(sqrt(2.0f0) * ones(Float32, d)) # Matrix d x d
prob = TerminalPDEProblem(g, f, μ_f, σ_f, X0, tspan)
hls = 10 + d # hidden layer size
opt = Flux.ADAM(0.01)  # optimizer
# sub-neural network approximating solutions at the desired point
u0 = Flux.Chain(Dense(d, hls, relu),
                Dense(hls, hls, relu),
                Dense(hls, 1))
# sub-neural network approximating the spatial gradients at time point
σᵀ∇u = Flux.Chain(Dense(d + 1, hls, relu),
                  Dense(hls, hls, relu),
                  Dense(hls, hls, relu),
                  Dense(hls, d))
pdealg = NNPDENS(u0, σᵀ∇u, opt=opt)
@time ans = solve(prob, pdealg, verbose=true, maxiters=100, trajectories=100,
                            alg=EM(), dt=1.2, pabstol=1f-2)


text:
module Poisson1DTests

using Gridap
using Test

u(x) = x[1]^2 + 1
∇u(x) = ∇(u)(x)
f(x) = -Δ(u)(x)

domain = (0,3)
cells = (10,)
model = CartesianDiscreteModel(domain,cells)

order = 2
V = FESpace(
  model=model,
  reffe=:Lagrangian,
  order=order,
  valuetype=Float64,
  conformity=:H1,
  dirichlet_tags=2)

U = TrialFESpace(V,u)

degree = 2*order
trian = Triangulation(model)
quad = CellQuadrature(trian,degree)

btrian = BoundaryTriangulation(model,1)
bquad = CellQuadrature(btrian,degree)
nb = get_normal_vector(btrian)

a(u,v) = ∇(v)⋅∇(u)
l(v) = v*f
t_Ω = AffineFETerm(a,l,trian,quad)

l_b(v) = v*(nb⋅∇u)
t_b = FESource(l_b,btrian,bquad)

op = AffineFEOperator(U,V,t_Ω,t_b)
uh = solve(op)

e = u - uh

l2(u) = inner(u,u)
sh1(u) = a(u,u)
h1(u) = sh1(u) + l2(u)

el2 = sqrt(sum( integrate(l2(e),trian,quad) ))
eh1 = sqrt(sum( integrate(h1(e),trian,quad) ))
ul2 = sqrt(sum( integrate(l2(uh),trian,quad) ))
uh1 = sqrt(sum( integrate(h1(uh),trian,quad) ))

@test el2/ul2 < 1.e-8
@test eh1/uh1 < 1.e-7

end # module


text:
# recipe.jl

using RecipesBase

@recipe function g(f::TaylorModel1)
    ffp = fp_rpa(f)
    fT = polynomial(ffp)
    Δ = remainder(ffp)
    ξ0 = ffp.x0

    seriesalpha --> 0.3
    seriestype := :shape

    xs = range(f.dom.lo, stop=f.dom.hi, length=100)
    evals = fT.(xs .- ξ0) .+ Δ

    # make polygon:
    xs = [xs; reverse(xs); xs[1]]
    ys = [inf.(evals); reverse(sup.(evals)); inf(evals[1])]

    xs, ys
end

@recipe function g(f::RTaylorModel1)
    ffp = fp_rpa(f)
    fT = polynomial(ffp)
    Δ = remainder(ffp)
    ξ0 = ffp.x0
    order = get_order(f)+1

    seriesalpha --> 0.5
    seriestype := :shape

    xs = range(f.dom.lo, stop=f.dom.hi, length=100)
    evals = fT.(xs .- ξ0)

    corrs = (xs .- ξ0) .^ order
    Δrel = Δ .* corrs
    evalslo = inf.(evals + Δrel)
    evalshi = sup.(evals + Δrel)

    xs = [xs; reverse(xs); xs[1]]
    ys = [evalslo; reverse(evalshi); evalslo[1]]

    xs, ys
end

@recipe function g(sol::TMSol; vars=(0,1), timediv=1)
    @assert length(vars) == length(unique(vars)) == 2 "`vars` must have 2 distinct integer indices"

    return _plot_intvbox(sol, vars=vars, timediv=timediv)
end

function _plot_intvbox(sol::TMSol; vars=(0,1), timediv=1)
    domT = mince_in_time(sol, var=0, timediv=timediv)

    if 0 ∈ vars
        tup0 = findfirst(0 ∈ vars)
        var1 = vars[3-tup0]
        v1 = _mince_in_time(sol, domT, var1, timediv)
        if tup0 == 1
            return @. IntervalBox(domT, v1)
        else
            return @. IntervalBox(v1, domT)
        end
    end

    var1, var2 = vars
    v1 = _mince_in_time(sol, domT, var1, timediv)
    v2 = _mince_in_time(sol, domT, var2, timediv)
    return @. IntervalBox(v1, v2)
end

"""
    mince_in_time(sol::TMSol; var=0, timediv=1) --> ::Vector{Interval}

For `var=0`, this function divides the time domain of each entry of `sol` in
`timediv` parts (`timediv==1` is the initial domain), and returns the time
intervals where the solution holds. This is useful for plotting or finding
specific events.
For `var ≥ 1`, this function evaluates the flowpipes at the split domain
intervals, which is useful to decrease the overapproximations associated
to the whole time domain.
"""
function mince_in_time(sol::TMSol; var::Int=0, timediv::Int=1)
    @assert timediv > 0 "`timediv must be 1 or larger"
    @assert 0 ≤ var ≤ get_numvars(sol)

    domT = _mince_in_time(sol, Val(true), timediv)
    var == 0 && return domT
    return _mince_in_time(sol, domT, var, timediv)
end

# Mince in time var (var == 0)
function _mince_in_time(sol::TMSol, ::Val{true}, timediv::Int=1)
    # Case timediv == 1
    if timediv == 1
        return expansion_point(sol) .+ domain(sol)
    end

    # Case timediv > 1
    domT = Vector{typeof(domain(sol,1))}(undef, timediv*length(sol))
    i0 = 1
    @inbounds for indT in eachindex(sol)
        i1 = indT*timediv
        domT[i0:i1] .= expansion_point(sol,indT) .+ mince(domain(sol,indT), timediv)
        i0 = i1 + 1
    end

    return domT
end

# Mince other var (var > 0)
function _mince_in_time(sol::TMSol, domT::Vector{Interval{T}}, var::Int,
        timediv::Int=1) where {T}
    N = get_numvars(sol)
    @assert 1 ≤ var ≤ N

    # Case timediv == 1
    if timediv == 1
        return getindex.(flowpipe(sol), var)
    end

    # Case timediv > 1
    vv = similar(domT)
    normalized_box = symmetric_box(N, Float64)
    δt = mince(domain(sol,1), timediv)

    i0 = 1
    @inbounds for indT in eachindex(sol)
        i1 = indT*timediv
        δt .= mince(domain(sol,indT), timediv)
        @. vv[i0:i1] = evaluate(evaluate(sol[indT,var], δt), (normalized_box,))
        i0 = i1 + 1
    end

    return vv
end


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

"""
Formulation described in:

    Morales-España, G., Latorre, J. M., & Ramos, A. (2013). Tight and compact
    MILP formulation for the thermal unit commitment problem. IEEE Transactions
    on Power Systems, 28(4), 4897-4908. DOI: https://doi.org/10.1109/TPWRS.2013.2251373
"""
module MorLatRam2013

import ..RampingFormulation
import ..StartupCostsFormulation

struct Ramping <: RampingFormulation end
struct StartupCosts <: StartupCostsFormulation end

end


text:
"""
Default parameters for all convenience methods that are exported to the end user.

See `OptRunController` for the description.
"""
const DefaultParameters = ParamsDict(
    :NumDimensions  => :NotSpecified, # Dimension of problem to be optimized
    :SearchRange    => (-1.0, 1.0), # Default search range to use per dimension unless specified
    :SearchSpace    => false, # Search space can be directly specified and will then take precedence over NumDimensions and SearchRange.
    :FitnessScheme  => MinimizingFitnessScheme, # fitness scheme to be used
    :TargetFitness => nothing, # optimal (target) fitness, if known

    :Method => :adaptive_de_rand_1_bin_radiuslimited,

    :MaxTime        => 0.0,
    :MaxFuncEvals   => 0,
    :MaxSteps       => 10000,
    :MaxStepsWithoutProgress => 10000,
    :MinDeltaFitnessTolerance => 1e-50,
    :FitnessTolerance => 1e-8,

    :MaxNumStepsWithoutFuncEvals => 100,

    :NumRepetitions => 1,     # Number of repetitions to run for each optimizer for each problem

    :TraceMode      => :compact,  # Print tracing information during the optimization
    :TraceInterval  => 0.50,  # Minimum number of seconds between consecutive trace messages printed to STDOUT
    :SaveTrace      => false,
    :SaveFitnessTraceToCsv => false, # Save a csv file with information about the major fitness improvement events (only the first event in each fitness magnitude class is saved)
    :SaveParameters => false, # Save parameters to a json file for later scrutiny

    :CallbackFunction => x -> x, # Function to callback to, here just the identity function.
    :CallbackInterval  => -1.0,  # Minimum number of seconds between consecutive callbacks. If <0.0 we never callback (which is the default).

    :RandomizeRngSeed => true, # Randomize the RngSeed value before using any random numbers.
    :RngSeed        => 1234,   # The specific random seed to set before any random numbers are generated. The seed is randomly selected if RandomizeRngSeed is true, and this parameter is updated with its actual value.

    :PopulationSize => 50
)

function check_and_create_search_space(params::Parameters)
    # If an explicit SearchSpace has been given we use that. It takes precedence.
    if haskey(params, :SearchSpace)
        ss = params[:SearchSpace]
        if isa(ss, SearchSpace)
            return ss
        elseif isa(ss, typeof([(0.0, 1.0)]))
            return RectSearchSpace(ss)
        elseif ss == false
            # silently fallthrough to the other means of search space specification
        else
            throw(ArgumentError("Using $(typeof(ss)) for SearchSpace is not supported"))
        end
    end

    # If no explicit search space has been given we must create one from the
    # other related parameters.
    if haskey(params, :SearchRange)
        sr = params[:SearchRange]
        # Check that a valid search range has been stated and create the search_space
        # based on it, or bail out.
        if isa(sr, Tuple{<:Real, <:Real})
            ndim = params[:NumDimensions]
            if ndim == :NotSpecified
                throw(ArgumentError("You MUST specify NumDimensions= in a solution when giving a SearchRange=$(sr)"))
            end
            return RectSearchSpace(params[:NumDimensions], sr)
        elseif isa(sr, Vector{Tuple{<:Real, <:Real}})
            return RectSearchSpace(sr)
        else
            throw(ArgumentError("Using $(typeof(sr)) for SearchRange is not supported."))
        end
    end

    # No valid search space have been specified => bail.
    throw(ArgumentError("Invalid search space specification ("*
                        "SearchSpace = $(params[:SearchSpace]), "*
                        "SearchRange = $(params[:SearchRange]), "*
                        "NumDimensions = $(params[:NumDimensions]))"))
end

function check_valid!(params::Parameters)
    # Check that max_time is larger than zero if it has been specified.
    if haskey(params, :MaxTime)
        if !isa(params[:MaxTime], Number) || params[:MaxTime] < 0.0
            throw(ArgumentError("MaxTime parameter must be a non-negative number"))
        elseif params[:MaxTime] > 0.0
            params[:MaxTime] = convert(Float64, params[:MaxTime])
            params[:MaxFuncEvals] = 0
            params[:MaxSteps] = 0
        end
    end

    # Check that a valid number of fevals has been specified. Print warning if higher than 1e8.
    if haskey(params,:MaxFuncEvals)
        if !isa(params[:MaxFuncEvals], Integer) || params[:MaxFuncEvals] < 0.0
            throw(ArgumentError("MaxFuncEvals parameter MUST be a non-negative integer"))
        elseif params[:MaxFuncEvals] > 0.0
            if params[:MaxFuncEvals] >= 1e8
                @warn("Number of allowed function evals is $(params[:MaxFuncEvals]); this can take a LONG time")
            end
            params[:MaxFuncEvals] = convert(Int, params[:MaxFuncEvals])
            params[:MaxSteps] = 0
        end
    end

    # Check that a valid number of iterations has been specified. Print warning if higher than 1e8.
    if haskey(params, :MaxSteps)
        if !isa(params[:MaxSteps], Number) || params[:MaxSteps] < 0.0
            throw(ArgumentError("The number of iterations (MaxSteps) MUST be a non-negative number"))
        elseif params[:MaxSteps] > 0.0
            if params[:MaxSteps] >= 1e8
                @warn("Number of allowed iterations is $(params[:MaxSteps]); this can take a LONG time")
            end
            params[:MaxSteps] = convert(Int, params[:MaxSteps])
        end
    end

    # Check that a valid population size has been given.
    if params[:PopulationSize] < 2
        # FIXME why? What if we use popsize of 1 for optimizers that improve on one solution?
        throw(ArgumentError("The population size MUST be at least 2"))
    end

    method = params[:Method]
    # Check that a valid method has been specified and then set up the optimizer
    if !isa(method, Symbol) || method ∉ BlackBoxOptim.MethodNames
        throw(ArgumentError("The method specified, $(method), is NOT among the valid methods:   $(MethodNames)"))
    end
end


text:
"""
    Layer(transformation, data, positional::AbstractVector, named::AbstractDictionary)

Algebraic object encoding a single layer of a visualization. It is composed of a dataset,
positional and named arguments, as well as a transformation to be applied to those.
`Layer` objects can be multiplied, yielding a novel `Layer` object, or added,
yielding a [`AlgebraOfGraphics.Layers`](@ref) object.
"""
Base.@kwdef struct Layer
    transformation::Any=identity
    data::Any=nothing
    positional::Arguments=Arguments()
    named::NamedArguments=NamedArguments()
end

transformation(f) = Layer(transformation=f)

data(df) = Layer(data=columns(df))

mapping(args...; kwargs...) = Layer(positional=collect(Any, args), named=NamedArguments(kwargs))

⨟(f, g) = f === identity ? g : g === identity ? f : g ∘ f

function Base.:*(l1::Layer, l2::Layer)
    transformation = l1.transformation ⨟ l2.transformation
    data = isnothing(l2.data) ? l1.data : l2.data
    positional = vcat(l1.positional, l2.positional)
    named = merge(l1.named, l2.named)
    return Layer(; transformation, data, positional, named)
end


text:
import Base:show, size, summary

const show_os = 12
const unindexed_fields = (:c, :n)

summary(R::RawData) = string(typeof(R), " with ", length(R.t),
                              " window", size(R.x,1) == 1 ? "" : "s")
summary(F::FFTData) = string(typeof(F), " with ", length(F.t), " fft",
                             size(F.fft,2) == 1 ? "" : "s")
summary(C::CorrData) = string(typeof(C), " with ", length(C.t),
                              " correlation", size(C.corr,1) == 1 ? "" : "s")

function show(io::IO, R::RawData)
  W = max(80, displaysize(io)[2]) - show_os
  w = min(W, 35)
  nc = size(R.x,1) == 0 ? 0 : 1
  N = min(nc, div(W-1, w))
  M = min(N+1, nc)
  rawlen,numraw = size(R.x)
  println(io, "RawData with ", length(R.t), " windows")
  Field = fieldnames(RawData)
  for f in Field
    if (f in unindexed_fields) == false
      targ = getfield(R, f)
      t = typeof(targ)
      fstr = uppercase(String(f))
      print(io, lpad(fstr, show_os-2), ": ")
      if f == :notes || f == :misc
        SeisIO.show_str(io, String[string(length(targ), " entries")], w, N)
      elseif (t <: AbstractFloat || t <: InstrumentPosition || t<: InstrumentResponse)
        println(io, repr("text/plain", targ, context=:compact=>true))
      elseif f == :t
        if length(R.t) > 0
          SeisIO.show_str(io, String[timestamp(R.t[1]), " (", string(length(R.t)), " windows)"], w, N)
        else
          SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
        end
      elseif f == :x
        SeisIO.show_str(io, String[summary(targ)], w, N)
      else
        SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
      end
    end
  end
  return nothing
end
show(R::RawData) = show(stdout, R)

function show(io::IO, F::FFTData)
  W = max(80, displaysize(io)[2]) - show_os
  w = min(W, 35)
  nc = size(F.fft,1) == 0 ? 0 : 1
  N = min(nc, div(W-1, w))
  M = min(N+1, nc)
  fftlen,numfft = size(F.fft)
  println(io, "FFTData with ", length(F.t), " ffts")
  Field = fieldnames(FFTData)
  for f in Field
    if (f in unindexed_fields) == false
      targ = getfield(F, f)
      t = typeof(targ)
      fstr = uppercase(String(f))
      print(io, lpad(fstr, show_os-2), ": ")
      if f == :notes || f == :misc
        SeisIO.show_str(io, String[string(length(targ), " entries")], w, N)
      elseif (t <: AbstractFloat || t <: InstrumentPosition || t<: InstrumentResponse)
        println(io, repr("text/plain", targ, context=:compact=>true))
      elseif f == :t
        if length(F.t) > 0
          SeisIO.show_str(io, String[timestamp(F.t[1]), " (", string(length(F.t)), " FFTs)"], w, N)
        else
          SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
        end
      elseif f == :fft
        SeisIO.show_str(io, String[summary(targ)], w, N)
      else
        SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
      end
    end
  end
  return nothing
end
show(F::FFTData) = show(stdout, F)

function show(io::IO, C::CorrData)
  W = max(80, displaysize(io)[2]) - show_os
  w = min(W, 35)
  nc = size(C.corr,1) == 0 ? 0 : 1
  N = min(nc, div(W-1, w))
  M = min(N+1, nc)
  corrlen,numcorr = size(C.corr)
  println(io, "CorrData with ", length(C.t), " Corrs")
  Field = fieldnames(CorrData)
  for f in Field
    if (f in unindexed_fields) == false
      targ = getfield(C, f)
      t = typeof(targ)
      fstr = uppercase(String(f))
      print(io, lpad(fstr, show_os-2), ": ")
      if f == :notes || f == :misc
        SeisIO.show_str(io, String[string(length(targ), " entries")], w, N)
      elseif (t <: AbstractFloat || t <: InstrumentPosition || t<: InstrumentResponse)
        println(io, repr("text/plain", targ, context=:compact=>true))
      elseif f == :t
        if length(C.t) > 0
          SeisIO.show_str(io, String[timestamp(C.t[1]), " (", string(length(C.t)), " Corrs)"], w, N)
        else
          SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
        end
      elseif f == :corr
        SeisIO.show_str(io, String[summary(targ)], w, N)
      else
        SeisIO.show_str(io, String[repr("text/plain", targ, context=:compact=>true)], w, N)
      end
    end
  end
  return nothing
end
show(C::CorrData) = show(stdout, C)


text:
module FileOps

export find_files, find_directories

# =========================================================================== #
find_files(dir::AbstractString, re=r".*") = return do_match(dir, re, isfile)
# --------------------------------------------------------------------------- #
find_directories(dir::AbstractString, re=r".*") = return do_match(dir, re, isdir)

# =========================================================================== #
function do_match(dir::AbstractString, re::Regex, f::Function)
    if !isdir(dir)
        error("Input is not a vaild directory path")
    end
    files = [joinpath(dir, x) for x in readdir(dir)]
    return filter(x->occursin(re, x) && f(x), files)
end
# =========================================================================== #
end # END MODULE