text:
# Base64 Encoder
# ==============

struct Base64Encoder <: Codec
    table::CodeTable64
    state::State
    buffer::Buffer
end

Base64Encoder(table::CodeTable64) = Base64Encoder(table, State(), Buffer(2))

"""
    Base64Encoder(;urlsafe::Bool=false)

Create a base64 encoding codec.

Arguments
---------
- `urlsafe`: use `-` and `_` as the last two values
"""
function Base64Encoder(;urlsafe::Bool=false)
    if urlsafe
        table = BASE64_URLSAFE
    else
        table = BASE64_STD
    end
    return Base64Encoder(table)
end

const Base64EncoderStream{S} = TranscodingStream{Base64Encoder,S} where S<:IO

"""
    Base64EncoderStream(stream::IO; kwargs...)

Create a base64 encoding stream (see `Base64Encoder` for `kwargs`).
"""
function Base64EncoderStream(stream::IO; kwargs...)
    return TranscodingStream(Base64Encoder(;kwargs...), stream)
end

function TranscodingStreams.startproc(
        codec :: Base64Encoder,
        state :: Symbol,
        error :: Error)
    start!(codec.state)
    empty!(codec.buffer)
    return :ok
end

function TranscodingStreams.process(
        codec  :: Base64Encoder,
        input  :: Memory,
        output :: Memory,
        error  :: Error)
    table = codec.table
    state = codec.state
    buffer = codec.buffer

    # Check if we can encode data.
    if !is_running(state)
        error[] = ArgumentError("encoding is already finished")
        return 0, 0, :error
    elseif output.size < 4
        # Need more output space.
        return 0, 0, :ok
    end

    # Load the first three bytes.
    i = j = 0
    while buffer.size < 2 && i < input.size
        buffer[buffer.size+=1] = input[i+=1]
    end
    b1 = b2 = b3 = 0x00
    npad = 0
    status = :ok
    if i < input.size
        b1 = buffer[1]
        b2 = buffer[2]
        b3 = input[i+=1]
    elseif input.size == 0
        # Found the end of the input.
        if buffer.size == 0
            finish!(state)
            return i, j, :end
        elseif buffer.size == 1
            b1 = buffer[1]
            npad = 2
        elseif buffer.size == 2
            b1 = buffer[1]
            b2 = buffer[2]
            npad = 1
        else
            @unreachable
        end
        status = :end
    else
        # Need more data to encode.
        return i, j, :ok
    end
    empty!(buffer)

    # Encode the body.
    @inbounds while true
        output[j+1] = encode(table, b1 >> 2          )
        output[j+2] = encode(table, b1 << 4 | b2 >> 4)
        output[j+3] = encode(table, b2 << 2 | b3 >> 6)
        output[j+4] = encode(table,           b3     )
        j += 4
        if i + 3 ≤ input.size && j + 4 ≤ output.size
            b1 = input[i+1]
            b2 = input[i+2]
            b3 = input[i+3]
            i += 3
        else
            break
        end
    end

    # Epilogue.
    while npad > 0
        output[j-npad+1] = table.padcode
        npad -= 1
    end
    if status == :end || status == :error
        finish!(state)
    end
    return i, j, status
end


text:
SensitivityAlg(args...;kwargs...) = @error("The SensitivtyAlg choice mechanism was completely overhauled. Please consult the local sensitivity documentation for more information")

abstract type AbstractForwardSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractAdjointSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end
abstract type AbstractSecondOrderSensitivityAlgorithm{CS,AD,FDT} <: DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT} end

struct ForwardSensitivity{CS,AD,FDT} <: AbstractForwardSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::Bool
end
Base.@pure function ForwardSensitivity(;
                                       chunk_size=0,autodiff=true,
                                       diff_type=Val{:central},
                                       autojacvec=autodiff)
  ForwardSensitivity{chunk_size,autodiff,diff_type}(autojacvec)
end

struct ForwardDiffSensitivity{CS,CTS} <: AbstractForwardSensitivityAlgorithm{CS,Nothing,Nothing}
end
Base.@pure function ForwardDiffSensitivity(;chunk_size=0,convert_tspan=nothing)
  ForwardDiffSensitivity{chunk_size,convert_tspan}()
end

struct BacksolveAdjoint{CS,AD,FDT,VJP} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  checkpointing::Bool
end
Base.@pure function BacksolveAdjoint(;chunk_size=0,autodiff=true,
                                      diff_type=Val{:central},
                                      autojacvec=autodiff,
                                      checkpointing=true)
  BacksolveAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec)}(autojacvec,checkpointing)
end

struct InterpolatingAdjoint{CS,AD,FDT,VJP} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  checkpointing::Bool
end
Base.@pure function InterpolatingAdjoint(;chunk_size=0,autodiff=true,
                                         diff_type=Val{:central},
                                         autojacvec=autodiff,
                                         checkpointing=false)
  InterpolatingAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec)}(autojacvec,checkpointing)
end

struct QuadratureAdjoint{CS,AD,FDT,VJP} <: AbstractAdjointSensitivityAlgorithm{CS,AD,FDT}
  autojacvec::VJP
  abstol::Float64
  reltol::Float64
  compile::Bool
end
Base.@pure function QuadratureAdjoint(;chunk_size=0,autodiff=true,
                                         diff_type=Val{:central},
                                         autojacvec=autodiff,abstol=1e-6,
                                         reltol=1e-3,compile=false)
  QuadratureAdjoint{chunk_size,autodiff,diff_type,typeof(autojacvec)}(autojacvec,abstol,reltol,compile)
end

struct TrackerAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end
struct ZygoteAdjoint <: AbstractAdjointSensitivityAlgorithm{nothing,true,nothing} end

abstract type VJPChoice end
struct ZygoteVJP <: VJPChoice end
struct TrackerVJP <: VJPChoice end
struct ReverseDiffVJP{compile} <: VJPChoice
  ReverseDiffVJP(compile=false) = new{compile}()
end

@inline convert_tspan(::ForwardDiffSensitivity{CS,CTS}) where {CS,CTS} = CTS
@inline alg_autodiff(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = AD
@inline get_chunksize(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = CS
@inline diff_type(alg::DiffEqBase.AbstractSensitivityAlgorithm{CS,AD,FDT}) where {CS,AD,FDT} = FDT
@inline function get_jacvec(alg::DiffEqBase.AbstractSensitivityAlgorithm)
  alg.autojacvec isa Bool ? alg.autojacvec : true
end
@inline ischeckpointing(alg::DiffEqBase.AbstractSensitivityAlgorithm, ::Vararg) = isdefined(alg, :checkpointing) ? alg.checkpointing : false
@inline ischeckpointing(alg::InterpolatingAdjoint, sol) = alg.checkpointing || !sol.dense
@inline compile_tape(vjp::ReverseDiffVJP{compile}) where compile = compile
@inline compile_tape(autojacvec::Bool) = false
@inline compile_tape(sensealg::QuadratureAdjoint) = sensealg.compile

struct ForwardDiffOverAdjoint{A} <: AbstractSecondOrderSensitivityAlgorithm{nothing,true,nothing}
  adjalg::A
end

ext: jl


text:
export getTetraMesh

function getTetraMesh(P::Array{Float64},T::Array{Int64})

	T = sort(T,2)

	(edges,faces) = getTetraMeshNumbering(T)

	nc = size(T,1)
  nf = size(faces,1)
  ne = size(edges,1)
	nn = size(P,1)

	Div   = spzeros(0,0)
	Curl  = spzeros(0,0)
	Grad  = spzeros(0,0)

  V = Float64[]
  F = Float64[]
  L = Float64[]

	return TetraMesh(P,T,nc,nf,ne,nn,Div,Curl,Grad,V,F,L,faces,edges)

end

function getTetraMesh(PFile::String,TFile::String)

  file  = open(PFile, "r")
  lines = readlines(file)
  close(file)

  n = length(lines)

  P = zeros(Float64, n, 3)

  for i = 1:n
    line = split(lines[i])
    P[i,1] = parse(Float64, line[1])
    P[i,2] = parse(Float64, line[2])
    P[i,3] = parse(Float64, line[3])
  end

  file  = open(TFile, "r")
  lines = readlines(file)
  close(file)

  n = length(lines)

  T = zeros(Int64, n, 4)

  for i = 1:n
    line = split(lines[i])
    T[i,1] = parse(Int64, line[1])
    T[i,2] = parse(Int64, line[2])
    T[i,3] = parse(Int64, line[3])
    T[i,4] = parse(Int64, line[4])
  end

  return getTetraMesh(P,T)

end

max_stars_repo_name: JuliaInv/JTetra

max_issues_repo_name: JuliaInv/JTetra

max_forks_repo_name: JuliaInv/JTetra


text:
# Basic structures

"""
    Orbital

Stores an orbital.

### Fields
- 'pos::SVector{3,Float64}': position of the orbital within the unit cell
- 'lm::NTuple{2, Int}': orbital quantum numbers: (total angular momentum, azimuthal angular momentum)
"""
struct Orbital
    pos::SVector{3,Float64}
    lm::NTuple{2, Int}
end

"""
    Orbital(pos::Vector; lm=(0, 0))

Initialized an orbital centered at `pos` of type (l, m)"
"""
function Orbital(pos::AbstractVector, lm=(0, 0))

    return Orbital(toSVector(pos, 3), lm)
end

const replacement_wf = Dict(
    :s => (0, 0),
    :px => (1, 1),
    :py => (1, -1),
    :pz => (1, 0),
    :dxy => (2, -2),
    :dx2y2 => (2, 2),
    :dxz => (2, 1),
    :dyz => (2, -1),
    :dz2 => (2, 0)
)

"""
    Orbital(pos::AbstractVector, wf::Symbol)

Initialized an orbital centered at `pos` of type `wf`.
Known types: `:s`, `:px`, `:py`, `:pz`, `:dxy`, `:dx2y2`, `:dxz`, `:dyz`, `:dz2`"
"""
function Orbital(pos::AbstractVector, wf::Symbol)

    @assert wf in keys(replacement_wf) "Shorthand for wavefunction of known type. Available shorthands: `:s`, `:px`, `:py`,`:pz`, `:dxy`, `:dx2y2`, `:dxz`, `:dyz`, `:dz2`."

    return Orbital(toSVector(pos, 3), replacement_wf[wf])
end

"""
    Hopping{D, T}

Stores a hopping with amplitude of type `T` in a `D`-dimensional lattice.

### Fields
- 'from::Int': start orbital
- 'to::Int': end orbital
- 'dir::SVector{D, Int}': direction of the hopping (in units of Bravais basis vectors.
- 'val::T': value of hopping
"""
struct Hopping{D, T}
    from::Int
    to::Int
    dir::SVector{D, Int}
    val::T
end

"""
    TBHamiltonian{D, T}

Stores a tight-binding system define in real space for a `D`-dimensional lattice (embedded in 3D) with hoppings of type `T`

### Fields
The tight-binsing system consists of:
- 'latbasis::SMatrix{D, D, Float64}': the basis for a `D`-dimensional Bravais lattice. The basis vectors are stored as columns of the matrix basis, i. e. a1 = basis[:, 1]
- 'orbitals::Vector{Orbital}': a list of orbitals within the unit cell
- 'hoppings::Vector{Hopping{D, T}': a list of hoppings (which includ intra- an intercell hoppings)
- 'sorted::Bool': whether the list of hoppings is been sorted by (from, to). Use of sparce matrices requires sorting.
"""
struct TBHamiltonian{D, T}
    latbasis::SMatrix{D, D, Float64}
    orbitals::Vector{Orbital}
    hoppings::Vector{Hopping{D, T}}
    sorted::Vector{Bool}
end


text:
module ConstraintSolver

using DataStructures
using Formatting
using JSON
using JuMP:
    @variable,
    @constraint,
    @objective,
    Model,
    optimizer_with_attributes,
    VariableRef,
    backend,
    set_optimizer,
    direct_model,
    optimize!,
    objective_value,
    set_lower_bound,
    set_upper_bound,
    termination_status
import JuMP.sense_to_set
import JuMP
import ConstraintProgrammingExtensions
using LightGraphs
using MathOptInterface
using MatrixNetworks
using Random
using Statistics
using StatsBase
using StatsFuns

const CPE = ConstraintProgrammingExtensions

const CS = ConstraintSolver
const CS_RNG = MersenneTwister(1)
const MOI = MathOptInterface
const MOIB = MathOptInterface.Bridges
const MOIBC = MathOptInterface.Bridges.Constraint
const MOIU = MOI.Utilities

const SVF = MOI.SingleVariable
const SAF = MOI.ScalarAffineFunction
const VAF = MOI.VectorAffineFunction

# indices
const VI = MOI.VariableIndex
const CI = MOI.ConstraintIndex

const VAR_TYPES = Union{MOI.ZeroOne,MOI.Integer}

include("types.jl")
const CoM = ConstraintSolverModel

include("tablelogger.jl")
include("options.jl")

include("type_inits.jl")

include("util.jl")
include("branching.jl")
include("traversing.jl")
include("lp_model.jl")
include("MOI_wrapper/MOI_wrapper.jl")
include("printing.jl")
include("logs.jl")
include("Variable.jl")
include("objective.jl")
include("constraints.jl")

include("constraints/all_different.jl")
include("constraints/complement.jl")
include("constraints/boolset.jl")
include("constraints/and.jl")
include("constraints/or.jl")
include("constraints/xor.jl")
include("constraints/xnor.jl")
include("constraints/linear_constraints.jl")
include("constraints/svc.jl")
include("constraints/all_equal.jl")
include("constraints/not_equal.jl")
include("constraints/table.jl")
include("constraints/activator_constraints.jl")
include("constraints/indicator.jl")
include("constraints/reified.jl")
include("constraints/geqset.jl")

include("pruning.jl")
include("simplify.jl")

"""
    get_inner_model(::Model) or get_inner_model(::MOIB.LazyBridgeOptimizer{Optimizer})

Return the ConstraintSolverModel for the Model or Optimizer
"""
function get_inner_model(m::Model)
    JuMP.backend(m).optimizer.model.model.inner
end

function get_inner_model(o::MOIB.LazyBridgeOptimizer{Optimizer})
    o.model.inner
end

"""
    fulfills_constraints(com::CS.CoM, vidx, value)

Return whether the model is still feasible after setting the variable at position `vidx` to `value`.
"""
function fulfills_constraints(com::CS.CoM, vidx, value)
    # variable doesn't have any constraint
    if vidx > length(com.subscription)
        return true
    end
    feasible = true
    constraints =  com.constraints
    for ci in com.subscription[vidx]
        constraint = constraints[ci]
        # only call if the function got initialized already
        if constraint.is_initialized
            feasible =
                still_feasible(com, constraint, constraint.fct, constraint.set, vidx, value)
            !feasible && break
        end
    end
    return feasible
end

"""
    set_pvals!(com::CS.CoM, constraint::Constraint)

Compute the possible values inside this constraint and set it as constraint.pvals
"""
function set_pvals!(com::CS.CoM, constraint::Constraint)
    indices = constraint.indices
    # nothing to do i.e x <= x will be x-x <= 0 => 0 <= 0 without variables
    length(indices) == 0 && return
    variables = Variable[v for v in com.search_space[indices]]
    pvals_intervals = Vector{NamedTuple}()
    push!(pvals_intervals, (from = variables[1].lower_bound, to = variables[1].upper_bound))
    for i in 1:length(indices)
        extra_from = variables[i].min
        extra_to = variables[i].max
        comp_inside = false
        for cpvals in pvals_intervals
            if extra_from >= cpvals.from && extra_to <= cpvals.to
                # completely inside the interval already
                comp_inside = true
                break
            elseif extra_from >= cpvals.from && extra_from <= cpvals.to
                extra_from = cpvals.to + 1
            elseif extra_to <= cpvals.to && extra_to >= cpvals.from
                extra_to = cpvals.from - 1
            end
        end
        if !comp_inside && extra_to >= extra_from
            push!(pvals_intervals, (from = extra_from, to = extra_to))
        end
    end
    pvals = collect((pvals_intervals[1].from):(pvals_intervals[1].to))
    for interval in pvals_intervals[2:end]
        pvals = vcat(pvals, collect((interval.from):(interval.to)))
    end
    constraint.pvals = pvals
    if constraint isa IndicatorConstraint || constraint isa ReifiedConstraint
        set_pvals!(com, constraint.inner_constraint)
    end
    if constraint isa BoolConstraint
        set_pvals!(com, constraint.lhs)
        set_pvals!(com, constraint.rhs)
    end
end

"""
    add_constraint!(com::CS.CoM, constraint::Constraint; set_pvals=true)

Add a constraint to the model and set pvals if `set_pvals=true` as well.
Pushes the new constraint to the subscription vector of the involved variables.
"""
function add_constraint!(com::CS.CoM, constraint::Constraint; set_pvals = true)
    @assert constraint.idx == length(com.constraints) + 1
    push!(com.constraints, constraint)
    set_pvals && set_pvals!(com, constraint)
    for vidx in constraint.indices
        push!(com.subscription[vidx], constraint.idx)
    end
end

"""
    get_best_bound(com::CS.CoM, backtrack_obj::BacktrackObj; vidx=0, lb=0, ub=0)

Return the best bound if setting the variable with idx: `vidx` to
    lb <= var[vidx] <= ub if vidx != 0
Without an objective function return 0.
"""
function get_best_bound(com::CS.CoM, backtrack_obj::BacktrackObj; vidx = 0, lb = 0, ub = 0)
    if com.sense == MOI.FEASIBILITY_SENSE
        return zero(com.best_bound)
    end
    best_bound = zero(com.best_bound)
    if com.options.lp_optimizer !== nothing
        best_bound = get_best_bound_lp(com, backtrack_obj, vidx, lb, ub)
    else
        best_bound = get_best_bound(com, backtrack_obj, com.objective, vidx, lb, ub)
    end
    # vidx != 0 can mean a temporary backtrack object that doesn't get added to backtrack_vec
    if vidx == 0 && backtrack_obj.status == :Open
        set_update_backtrack_pq!(com, backtrack_obj; best_bound = best_bound)
    end
    return best_bound
end

"""
    checkout_from_to!(com::CS.CoM, from_nidx::Int, to_nidx::Int)

Change the state of the search space given the current position in the tree (`from_nidx`) and the index we want
to change to (`to_nidx`)
"""
function checkout_from_to!(com::CS.CoM, from_nidx::Int, to_nidx::Int)
    backtrack_vec = com.backtrack_vec
    from = backtrack_vec[from_nidx]
    to = backtrack_vec[to_nidx]
    if to.parent_idx == from.idx
        return
    end
    reverse_pruning!(com, from.idx)

    prune_steps = Vector{Int}()
    # first go to same level if new is higher in the tree
    if to.depth < from.depth
        depth = from.depth
        parent_nidx = from.parent_idx
        parent = backtrack_vec[parent_nidx]
        while to.depth < depth
            reverse_pruning!(com, parent_nidx)
            parent = backtrack_vec[parent_nidx]
            parent_nidx = parent.parent_idx
            depth -= 1
        end
        if parent_nidx == to.parent_idx
            return
        else
            from = parent
        end
    elseif from.depth < to.depth
        depth = to.depth
        parent_nidx = to.parent_idx
        parent = backtrack_vec[parent_nidx]
        while from.depth < depth
            pushfirst!(prune_steps, parent_nidx)
            parent = backtrack_vec[parent_nidx]
            parent_nidx = parent.parent_idx
            depth -= 1
        end

        to = parent
        if backtrack_vec[prune_steps[1]].parent_idx == from.parent_idx
            !isempty(prune_steps) && restore_prune!(com, prune_steps)
            return
        end
    end
    @assert from.depth == to.depth
    # same depth but different parent
    # => level up until same parent
    while from.parent_idx != to.parent_idx
        reverse_pruning!(com, from.parent_idx)
        from = backtrack_vec[from.parent_idx]

        pushfirst!(prune_steps, to.parent_idx)
        to = backtrack_vec[to.parent_idx]
    end

    !isempty(prune_steps) && restore_prune!(com, prune_steps)
end

"""
    set_state_to_best_sol!(com::CS.CoM, last_backtrack_id::Int)

Set the state of the model to the best solution found
"""
function set_state_to_best_sol!(com::CS.CoM, last_backtrack_id::Int)
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    backtrack_vec = com.backtrack_vec
    # find one of the best solutions
    sol, sol_id = findmin([
        backtrack_vec[sol.backtrack_id].best_bound * obj_factor for sol in com.solutions
    ])
    backtrack_id = com.solutions[sol_id].backtrack_id
    checkout_from_to!(com, last_backtrack_id, backtrack_id)
    # prune the last step as checkout_from_to! excludes the to part
    restore_prune!(com, backtrack_id)
end

"""
    addBacktrackObj2Backtrack_vec!(backtrack_vec, backtrack_obj, com::CS.CoM)

Add a backtrack object to the backtrack vector and create necessary vectors and maybe include it in the logs
"""
function addBacktrackObj2Backtrack_vec!(backtrack_vec, backtrack_obj, com::CS.CoM)
    push!(backtrack_vec, backtrack_obj)
    @assert length(backtrack_vec) == backtrack_obj.idx
    add2priorityqueue(com, backtrack_obj)

    create_log_node(com)
end

"""
    backtrack_vec::Vector{BacktrackObj{T}}, com::CS.CoM{T}, parent_idx, depth, vidx;
     check_bound=false, only_one=false, compute_bound=true)

Create two branches with two additional `BacktrackObj`s and add them to `backtrack_vec`.
Normally bounds are not checked so it gets added even if it's worse.
One can use `only_one` if one wants to probe instead one strain instead of creating the two
    possible children.
If the bound is not important one can set `compute_bound = false` to improve running time.
"""
function add2backtrack_vec!(
    backtrack_vec::Vector{BacktrackObj{T}},
    com::CS.CoM{T},
    parent_idx,
    vidx;
    check_bound = false,
    only_one = false,
    compute_bound = true,
) where {T<:Real}
    @assert !check_bound || compute_bound
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    left_lb, left_ub, right_lb, right_ub =
        get_split_pvals(com, com.branch_split, com.search_space[vidx])

    #=
        Check whether the new node is needed which depends on
        - Is there a solution already?
            - no => Add
        - Do we want all solutions?
            - yes => Add
        - Do we want all optimal solutions?
            - yes => Add if better or same as previous optimal one
    =#

    # left branch
    backtrack_obj = new_BacktrackObj(com, parent_idx, vidx, left_lb, left_ub)
    if compute_bound
        backtrack_obj.best_bound =
            get_best_bound(com, backtrack_obj; vidx = vidx, lb = left_lb, ub = left_ub)
    end
    # only include nodes which have a better objective than the current best solution if one was found already
    if com.options.all_solutions ||
       !check_bound ||
       length(com.solutions) == 0 ||
       backtrack_obj.best_bound * obj_factor < com.best_sol * obj_factor ||
       com.options.all_optimal_solutions &&
       backtrack_obj.best_bound * obj_factor <= com.best_sol * obj_factor

        addBacktrackObj2Backtrack_vec!(backtrack_vec, backtrack_obj, com)
    end
    only_one && return
    @assert left_ub < right_lb
    # right branch
    backtrack_obj = new_BacktrackObj(com, parent_idx, vidx, right_lb, right_ub)
    if compute_bound
        backtrack_obj.best_bound =
            get_best_bound(com, backtrack_obj; vidx = vidx, lb = right_lb, ub = right_ub)
    end
    if com.options.all_solutions ||
       !check_bound ||
       length(com.solutions) == 0 ||
       backtrack_obj.best_bound * obj_factor < com.best_sol ||
       com.options.all_optimal_solutions &&
       backtrack_obj.best_bound * obj_factor <= com.best_sol * obj_factor

        addBacktrackObj2Backtrack_vec!(backtrack_vec, backtrack_obj, com)
    end
end

"""
    set_bounds!(com, backtrack_obj)

Set lower/upper bounds for the current variable index `backtrack_obj.vidx`.
Return if simple removable is still feasible
"""
function set_bounds!(com, backtrack_obj)
    vidx = backtrack_obj.vidx
    !remove_above!(com, com.search_space[vidx], backtrack_obj.ub) && return false
    !remove_below!(com, com.search_space[vidx], backtrack_obj.lb) && return false
    return true
end

"""
    add_new_solution!(com::CS.CoM, backtrack_vec::Vector{BacktrackObj{T}}, backtrack_obj::BacktrackObj{T}, log_table) where T <: Real

A new solution was found.
- Add it to the solutions objects
Return true if backtracking can be stopped
"""
function add_new_solution!(
    com::CS.CoM,
    backtrack_vec::Vector{BacktrackObj{T}},
    backtrack_obj::BacktrackObj{T},
    log_table,
) where {T<:Real}
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    find_more_solutions = com.options.all_solutions || com.options.all_optimal_solutions

    new_sol = get_best_bound(com, backtrack_obj)
    backtrack_obj.best_bound = new_sol
    if length(com.solutions) == 0 || obj_factor * new_sol <= obj_factor * com.best_sol
        # also push it to the solutions object
        new_sol_obj = Solution(new_sol, CS.value.(com.search_space), backtrack_obj.idx)
        push!(com.solutions, new_sol_obj)
        com.best_sol = new_sol
        log_table && (last_table_row = update_table_log(com, backtrack_vec; force = true))
        if com.best_sol == com.best_bound && !find_more_solutions
            return true
        end
        # set all nodes to :Worse if they can't achieve a better solution
        for bo in backtrack_vec
            if bo.status == :Open && obj_factor * bo.best_bound >= obj_factor * com.best_sol
                bo.status = :Worse
            end
        end
    else # if new solution was found but it's worse
        log_table && (last_table_row = update_table_log(com, backtrack_vec; force = true))
        if com.options.all_solutions
            new_sol_obj = Solution(new_sol, CS.value.(com.search_space), backtrack_obj.idx)
            push!(com.solutions, new_sol_obj)
        end
    end
    # change the traverse strategy for example if it was :DBFS and we found the first solution
    old_traverse_strategy = com.traverse_strategy
    com.traverse_strategy = get_traverse_strategy(com; options = com.options)
    if com.traverse_strategy != old_traverse_strategy
        changed_traverse_strategy!(com, old_traverse_strategy)
    end
    return false
end

"""
    checkout_new_node!(com::CS.CoM, last_id, new_id)

If last id is not 0 then changes from last_id to new_id and sets `com.c_backtrack_idx`
"""
function checkout_new_node!(com::CS.CoM, last_id, new_id)
    if last_id != 0
        com.c_backtrack_idx = 0
        checkout_from_to!(com, last_id, new_id)
        com.c_backtrack_idx = new_id
    end
end

"""
    found_best_node(com::CS.CoM)

Return whether a optimal solution was found
"""
function found_best_node(com::CS.CoM)
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    return length(com.solutions) > 0 &&
           obj_factor * com.best_bound >= obj_factor * com.best_sol
end

"""
    handle_infeasible!(com::CS.CoM; finish_pruning=false)

Handle infeasibility:
- finish pruning if `finish_pruning` is true
- log if desired
- increase `backtrack_reverses`

Return true to make calls like `!feasible && handle_infeasible!(com) && continue` possible
"""
function handle_infeasible!(com::CS.CoM; finish_pruning = false)
    # need to call as some function might have pruned something.
    # Just need to be sure that we save the latest states
    finish_pruning && call_finished_pruning!(com)
    last_backtrack_id = com.c_backtrack_idx
    com.backtrack_vec[last_backtrack_id].is_feasible = false
    update_log_node!(com, last_backtrack_id; feasible = false)
    com.info.backtrack_reverses += 1
    return true
end

function solve_with_backtrack!(com, max_bt_steps; sorting = true)
    com.info.backtrack_fixes = 1

    log_table = false
    if :Table in com.options.logging
        log_table = true
    end

    check_bounds = com.sense != MOI.FEASIBILITY_SENSE
    status, last_backtrack_id =
        backtrack!(com, max_bt_steps; sorting = sorting, log_table = log_table, check_bounds = check_bounds)

    status != :TBD && return status

    if length(com.solutions) > 0
        set_state_to_best_sol!(com, last_backtrack_id)
        com.best_bound = com.best_sol
        if time() - com.start_time > com.options.time_limit
            return :Time
        else
            return :Solved
        end
    end

    if time() - com.start_time > com.options.time_limit
        return :Time
    else
        return :Infeasible
    end
end

"""
    backtrack!(com::CS.CoM, max_bt_steps; sorting=true)

Start backtracking and stop after `max_bt_steps`.
If `sorting` is set to `false` the same ordering is used as when used without objective this has only an effect when an objective is used.
Return :Solved or :Infeasible if proven or `:NotSolved` if interrupted by `max_bt_steps`.
"""
function backtrack!(
    com::CS.CoM,
    max_bt_steps;
    sorting = true,
    log_table = true,
    first_parent_idx = 1,
    single_path = false,
    compute_bounds = true,
    check_bounds = true,
    cb_finished_pruning = (args...) -> nothing,
)

    branch_var = get_next_branch_variable(com)
    branch_var.is_solution && return :Solved, first_parent_idx
    !branch_var.is_feasible && return :Infeasible, first_parent_idx

    backtrack_vec = com.backtrack_vec

    find_more_solutions = com.options.all_solutions || com.options.all_optimal_solutions

    log_table && println(get_header(com.options.table))

    backtrack_vec = com.backtrack_vec

    add2backtrack_vec!(
        backtrack_vec,
        com,
        first_parent_idx, # parent_idx
        branch_var.vidx;
        only_one = single_path,
        compute_bound = compute_bounds,
        check_bound = check_bounds,
    )
    last_backtrack_id = first_parent_idx

    started = true
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1

    while length(backtrack_vec) > 0
        # get next open backtrack object
        if !started
            # close the previous backtrack object
            close_node!(com, last_backtrack_id)
            update_log_node!(com, last_backtrack_id)
        end
        # run at least once so that everything is well defined
        if !started && time() - com.start_time > com.options.time_limit
            break
        end

        if !started && compute_bounds
            update_best_bound!(com)
        end
        found, backtrack_obj = get_next_node(com, backtrack_vec, sorting)
        !found && break

        # there is no better node => return best solution
        !find_more_solutions && found_best_node(com) && break

        vidx = backtrack_obj.vidx

        com.c_backtrack_idx = backtrack_obj.idx
        checkout_new_node!(com, last_backtrack_id, backtrack_obj.idx)

        # if backtracking was started
        # => remove all values which are root infeasible
        # this is done here instead of at the root node to avoid removing variables
        # in nodes which already have children
        if started
            for root_infeasible in com.root_infeasible_vars
                for val in (root_infeasible.lb):(root_infeasible.ub)
                    if has(com.search_space[root_infeasible.vidx], val)
                        if !rm!(com, com.search_space[root_infeasible.vidx], val)
                            return :Infeasible, last_backtrack_id
                        end
                    end
                end
            end
        end

        started = false
        last_backtrack_id = backtrack_obj.idx

        # limit the variable bounds
        if !set_bounds!(com, backtrack_obj)
            update_log_node!(com, last_backtrack_id; feasible = false)
            continue
        end

        constraints = com.constraints[com.subscription[vidx]]
        com.info.backtrack_fixes += 1

        # prune completely start with all that changed by the fix or by updating best bound
        feasible = prune!(com)
        !feasible && handle_infeasible!(com; finish_pruning = true) && continue
        call_finished_pruning!(com)
        cb_finished_pruning(com)

        if log_table
            last_table_row = update_table_log(com, backtrack_vec)
        end

        branch_var = get_next_branch_variable(com)
        if branch_var.is_solution
            finished = add_new_solution!(com, backtrack_vec, backtrack_obj, log_table)
            update_log_node!(com, last_backtrack_id)
            if finished
                # close the previous backtrack object
                close_node!(com, last_backtrack_id)
                return :Solved, last_backtrack_id
            end
            continue
        end

        if com.info.backtrack_fixes > max_bt_steps
            close_node!(com, last_backtrack_id)
            update_log_node!(com, last_backtrack_id)
            return :NotSolved, last_backtrack_id
        end

        update_log_node!(com, backtrack_obj.idx)

        last_backtrack_obj = backtrack_vec[last_backtrack_id]
        add2backtrack_vec!(
            backtrack_vec,
            com,
            last_backtrack_obj.idx,
            branch_var.vidx;
            only_one = single_path,
            compute_bound = compute_bounds,
            check_bound = check_bounds,
        )
    end

    close_node!(com, last_backtrack_id)
    update_log_node!(com, last_backtrack_id)
    return :TBD, last_backtrack_id
end

"""
    set_in_all_different!(com::CS.CoM)

Set `constraint.in_all_different` if all variables in the constraint are part of the same `all_different` constraint.
"""
function set_in_all_different!(com::CS.CoM; constraints = com.constraints)
    for constraint in constraints
        if :in_all_different in fieldnames(typeof(constraint))
            if !constraint.in_all_different
                subscriptions_idxs =
                    [[i for i in com.subscription[v]] for v in constraint.indices]
                length(subscriptions_idxs) == 0 && continue
                intersects = intersect(subscriptions_idxs...)

                for i in intersects
                    if isa(com.constraints[i].set, CPE.AllDifferent)
                        constraint.in_all_different = true
                        push!(com.constraints[i].sub_constraint_idxs, constraint.idx)
                    end
                end
            end
        end
    end
end

"""
    sort_solutions!(com::CS.CoM)

Order com.solutions by objective
"""
function sort_solutions!(com::CS.CoM)
    obj_factor = com.sense == MOI.MIN_SENSE ? 1 : -1
    sort!(com.solutions, by = s -> s.incumbent * obj_factor)
end

function print_info(com::CS.CoM)
    println("# Variables: ", length(com.search_space))
    println("# Constraints: ", length(com.constraints))
    for field in fieldnames(CS.NumberConstraintTypes)
        field_str = uppercasefirst(String(field))
        val = getfield(com.info.n_constraint_types, field)
        val != 0 && println(" - # $field_str: $val")
    end
    println()
end

"""
    solve!(com::CS.CoM, options::SolverOptions)

Solve the constraint model based on the given settings.
"""
function solve!(com::CS.CoM)
    options = com.options
    Random.seed!(CS_RNG, options.seed)
    backtrack = options.backtrack
    max_bt_steps = options.max_bt_steps
    backtrack_sorting = options.backtrack_sorting
    keep_logs = options.keep_logs
    if options.traverse_strategy == :Auto
        options.traverse_strategy = get_auto_traverse_strategy(com)
    end
    if options.branch_strategy == :Auto
        options.branch_strategy = get_auto_branch_strategy(com)
    end
    com.traverse_strategy = get_traverse_strategy(com; options = options)
    com.branch_strategy = get_branch_strategy(; options = options)
    com.branch_split = get_branch_split(; options = options)

    if :Info in com.options.logging
        print_info(com)
    end
    com.start_time = time()

    !set_init_fixes!(com) && return :Infeasible
    set_in_all_different!(com)

    # initialize constraints if `init_constraint!` exists for the constraint
    !init_constraints!(com) && return :Infeasible

    com.input[:logs] = keep_logs
    if keep_logs
        com.init_search_space = deepcopy(com.search_space)
    end

    # check for better constraints
    if options.simplify
        added_con_idxs = simplify!(com)
        if length(added_con_idxs) > 0
            set_in_all_different!(com; constraints = com.constraints[added_con_idxs])
            !init_constraints!(com; constraints = com.constraints[added_con_idxs]) &&
                return :Infeasible
            !update_init_constraints!(com; constraints = com.constraints[added_con_idxs]) &&
                return :Infeasible
            recompute_subscriptions(com)
        end
    end

    options.no_prune && return :NotSolved

    # check if all feasible even if for example everything is fixed
    feasible = prune!(com; pre_backtrack = true, initial_check = true)
    # finished pruning will be called in second call a few lines down...

    if !feasible
        com.solve_time = time() - com.start_time
        return :Infeasible
    end
    if all(v -> isfixed(v), com.search_space)
        # root node is the first backtrack obj
        dummy_backtrack_obj = BacktrackObj(com)
        com.best_bound = get_best_bound(com, dummy_backtrack_obj)

        backtrack_vec = com.backtrack_vec
        addBacktrackObj2Backtrack_vec!(backtrack_vec, dummy_backtrack_obj, com)

        com.best_sol = com.best_bound
        com.solve_time = time() - com.start_time
        new_sol_obj = Solution(com.best_sol, CS.value.(com.search_space), 1)
        push!(com.solutions, new_sol_obj)
        return :Solved
    end
    feasible = prune!(com; pre_backtrack = true)
    call_finished_pruning!(com)

    # root node is the first backtrack obj
    dummy_backtrack_obj = BacktrackObj(com)
    dummy_backtrack_obj.step_nr = com.c_step_nr
    com.best_bound = get_best_bound(com, dummy_backtrack_obj)

    backtrack_vec = com.backtrack_vec
    addBacktrackObj2Backtrack_vec!(backtrack_vec, dummy_backtrack_obj, com)

    if keep_logs
        com.logs[1].step_nr = 1
    end

    if !feasible
        com.solve_time = time() - com.start_time
        return :Infeasible
    end

    if all(v -> isfixed(v), com.search_space)
        com.best_sol = com.best_bound
        com.solve_time = time() - com.start_time
        new_sol_obj = Solution(com.best_sol, CS.value.(com.search_space), 1)
        push!(com.solutions, new_sol_obj)
        return :Solved
    end
    if backtrack
        com.info.backtracked = true
        if time() - com.start_time > com.options.time_limit
            com.solve_time = time() - com.start_time
            return :Time
        end
        keep_logs && update_log_node!(com, 1)
        status = solve_with_backtrack!(com, max_bt_steps; sorting = backtrack_sorting)
        com.solve_time = time() - com.start_time
        return status
    else
        @info "Backtracking is turned off."
        com.solve_time = time() - com.start_time
        return :NotSolved
    end
end

@deprecate AllDifferentSet AllDifferent false
@deprecate EqualSet AllEqual false

end # module


text:
abstract type AbstractHydroFormulation <: AbstractDeviceFormulation end

abstract type AbstractHydroDispatchFormulation <: AbstractHydroFormulation end

abstract type AbstractHydroUnitCommitment <: AbstractHydroFormulation end

struct HydroFixed <: AbstractHydroFormulation end

struct HydroDispatchRunOfRiver <: AbstractHydroDispatchFormulation end

struct HydroDispatchSeasonalFlow <: AbstractHydroDispatchFormulation end

struct HydroCommitmentRunOfRiver <: AbstractHydroUnitCommitment end

struct HydroCommitmentSeasonalFlow <: AbstractHydroUnitCommitment end

########################### Hydro generation variables #################################

function activepower_variables!(canonical::Canonical,
                               devices::IS.FlattenIteratorWrapper{H}) where H<:PSY.HydroGen

    add_variable(canonical,
                 devices,
                 Symbol("P_$(H)"),
                 false,
                 :nodal_balance_active;
                 lb_value = d -> d.tech.activepowerlimits.min,
                 ub_value = d -> d.tech.activepowerlimits.max,
                 init_value = d -> PSY.get_activepower(PSY.get_tech(d)))

    return

end

function reactivepower_variables!(canonical::Canonical,
                                 devices::IS.FlattenIteratorWrapper{H}) where H<:PSY.HydroGen

    add_variable(canonical,
                 devices,
                 Symbol("Q_$(H)"),
                 false,
                 :nodal_balance_reactive,
                 ub_value = d -> d.tech.reactivepowerlimits.max,
                 lb_value = d -> d.tech.reactivepowerlimits.min,
                 init_value = d -> d.tech.reactivepower)

    return

end

"""
This function add the variables for power generation commitment to the model
"""
function commitment_variables!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{H}) where {H<:PSY.HydroGen}

    time_steps = model_time_steps(canonical)
    var_names = [Symbol("ON_$(H)"), Symbol("START_$(H)"), Symbol("STOP_$(H)")]

    for v in var_names
        add_variable(canonical, devices, v, true)
    end

    return

end

### Constraints for Thermal Generation without commitment variables ####
"""
This function adds the Commitment Status constraint when there are CommitmentVariables
"""
function commitment_constraints!(canonical::Canonical,
                                 devices::IS.FlattenIteratorWrapper{H},
                                 device_formulation::Type{D},
                                 system_formulation::Type{S}) where {H<:PSY.HydroGen,
                                                                     D<:AbstractHydroUnitCommitment,
                                                                     S<:PM.AbstractPowerModel}

    key = ICKey(DeviceStatus, H)

    if !(key in keys(canonical.initial_conditions))
        error("Initial status conditions not provided. This can lead to unwanted results")
    end

    device_commitment(canonical,
                     canonical.initial_conditions[key],
                     Symbol("commitment_$(H)"),
                     (Symbol("START_$(H)"),
                      Symbol("STOP_$(H)"),
                      Symbol("ON_$(H)"))
                      )

    return

end

####################################### Reactive Power Constraints #########################
function reactivepower_constraints!(canonical::Canonical,
                                    devices::IS.FlattenIteratorWrapper{H},
                                    device_formulation::Type{AbstractHydroDispatchFormulation},
                                    system_formulation::Type{<:PM.AbstractPowerModel}) where H<:PSY.HydroGen

    range_data = Vector{NamedMinMax}(undef, length(devices))

    for (ix, d) in enumerate(devices)
        tech = PSY.get_tech(d)
        name = PSY.get_name(d)
        if isnothing(PSY.get_reactivepowerlimits(tech))
            limits = (min = 0.0, max = 0.0)
            range_data[ix] = (PSY.get_name(d), limits)
            @warn("Reactive Power Limits of $(name) are nothing. Q_$(name) is set to 0.0")
        else
            range_data[ix] = (name, PSY.get_reactivepowerlimits(tech))
        end
    end

    device_range(canonical,
                range_data,
                Symbol("reactiverange_$(H)"),
                Symbol("Q_$(H)"))

    return

end

######################## output constraints without Time Series ############################
function _get_time_series(canonical::Canonical,
                          devices::IS.FlattenIteratorWrapper{<:PSY.HydroGen})

    initial_time = model_initial_time(canonical)
    use_forecast_data = model_uses_forecasts(canonical)
    parameters = model_has_parameters(canonical)
    time_steps = model_time_steps(canonical)
    device_total = length(devices)
    ts_data_active = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)
    ts_data_reactive = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)

    for (ix, device) in enumerate(devices)
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        tech = PSY.get_tech(device)
        # pf = sin(acos(PSY.get_powerfactor(PSY.get_tech(device))))
        active_power = use_forecast_data ? PSY.get_rating(tech) : PSY.get_activepower(device)
        if use_forecast_data
            ts_vector = TS.values(PSY.get_data(PSY.get_forecast(PSY.Deterministic,
                                                                device,
                                                                initial_time,
                                                                "rating")))
        else
            ts_vector = ones(time_steps[end])
        end
        ts_data_active[ix] = (name, bus_number, active_power, ts_vector)
        ts_data_reactive[ix] = ts_data_active[ix] # (name, bus_number, active_power * pf, ts_vector)
    end

    return ts_data_active, ts_data_reactive

end

function activepower_constraints!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{H},
                                device_formulation::Type{<:AbstractHydroDispatchFormulation},
                                system_formulation::Type{<:PM.AbstractPowerModel}) where H<:PSY.HydroGen

    parameters = model_has_parameters(canonical)
    use_forecast_data = model_uses_forecasts(canonical)

    if !parameters && !use_forecast_data
        range_data = [(PSY.get_name(d), (min = 0.0, max = PSY.get_rating(PSY.get_tech(d)))) for d in devices]
        device_range(canonical,
                    range_data,
                    Symbol("activerange_$(H)"),
                    Symbol("P_$(H)"))
        return
    end

    ts_data_active, _ = _get_time_series(canonical, devices)
    if parameters
        device_timeseries_param_ub(canonical,
                            ts_data_active,
                            Symbol("activerange_$(H)"),
                            UpdateRef{H}(:rating),
                            Symbol("P_$(H)"))
    else
        device_timeseries_ub(canonical,
                            ts_data_active,
                            Symbol("activerange_$(H)"),
                            Symbol("P_$(H)"))
    end

    return

end

function activepower_constraints!(canonical::Canonical,
                                devices::IS.FlattenIteratorWrapper{H},
                                device_formulation::Type{<:AbstractHydroUnitCommitment},
                                system_formulation::Type{<:PM.AbstractPowerModel}) where H<:PSY.HydroGen

    parameters = model_has_parameters(canonical)
    use_forecast_data = model_uses_forecasts(canonical)

    if !parameters && !use_forecast_data
        range_data = [(PSY.get_name(d), (min = 0.0, max = PSY.get_rating(PSY.get_tech(d)))) for d in devices]
        device_semicontinuousrange(canonical,
                                    range_data,
                                    Symbol("activerange_$(H)"),
                                    Symbol("P_$(H)"),
                                    Symbol("ON_$(H)"))
        return
    end

    ts_data_active, _ = _get_time_series(canonical, devices)
    if parameters
        device_timeseries_ub_bigM(canonical,
                            ts_data_active,
                            Symbol("activerange_$(H)"),
                            Symbol("P_$(H)"),
                            UpdateRef{H}(:rating),
                            Symbol("ON_$(H)"))
    else
        device_timeseries_ub_bin(canonical,
                            ts_data_active,
                            Symbol("activerange_$(H)"),
                            Symbol("P_$(H)"),
                            Symbol("ON_$(H)"))
    end

    return

end

########################## Make initial Conditions for a Model #############################
function initial_conditions!(canonical::Canonical,
                            devices::IS.FlattenIteratorWrapper{H},
                            device_formulation::Type{<:AbstractHydroUnitCommitment}) where {H<:PSY.HydroGen}

    status_init(canonical, devices)
    output_init(canonical, devices)
    duration_init(canonical, devices)

    return

end

function initial_conditions!(canonical::Canonical,
                            devices::IS.FlattenIteratorWrapper{H},
                            device_formulation::Type{D}) where {H<:PSY.HydroGen,
                                                                D<:AbstractHydroDispatchFormulation}

    output_init(canonical, devices)

    return

end

########################## Addition of to the nodal balances ###############################
function nodal_expression!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{H},
                           system_formulation::Type{<:PM.AbstractPowerModel}) where H<:PSY.HydroGen

    parameters = model_has_parameters(canonical)
    ts_data_active, ts_data_reactive = _get_time_series(canonical, devices)

    if parameters
        include_parameters(canonical,
                           ts_data_active,
                           UpdateRef{H}(:rating),
                           :nodal_balance_active)
        include_parameters(canonical,
                           ts_data_reactive,
                           UpdateRef{H}(:rating),
                           :nodal_balance_reactive)
        return
    end

    for t in model_time_steps(canonical)
        for device_value in ts_data_active
            _add_to_expression!(canonical.expressions[:nodal_balance_active],
                            device_value[2],
                            t,
                            device_value[3]*device_value[4][t])
        end
        for device_value in ts_data_reactive
            _add_to_expression!(canonical.expressions[:nodal_balance_reactive],
                            device_value[2],
                            t,
                            device_value[3]*device_value[4][t])
        end
    end

    return

end

function nodal_expression!(canonical::Canonical,
                           devices::IS.FlattenIteratorWrapper{H},
                           system_formulation::Type{<:PM.AbstractActivePowerModel}) where H<:PSY.HydroGen

    parameters = model_has_parameters(canonical)
    ts_data_active, _  = _get_time_series(canonical, devices)

    if parameters
        include_parameters(canonical,
                           ts_data_active,
                           UpdateRef{H}(:rating),
                           :nodal_balance_active)
        return
    end

    for t in model_time_steps(canonical)
        for device_value in ts_data_active
            _add_to_expression!(canonical.expressions[:nodal_balance_active],
                            device_value[2],
                            t,
                            device_value[3]*device_value[4][t])
        end
    end

    return

end

##################################### Hydro generation cost ############################
function cost_function(canonical::Canonical,
                       devices::IS.FlattenIteratorWrapper{PSY.HydroDispatch},
                       device_formulation::Type{D},
                       system_formulation::Type{<:PM.AbstractPowerModel}) where D<:AbstractHydroFormulation

    add_to_cost(canonical,
                devices,
                Symbol("P_HydroDispatch"),
                :fixed,
                -1.0)

    return

end

##################################### Water/Energy Budget Constraint ############################
function _get_budget(canonical::Canonical,
                    devices::IS.FlattenIteratorWrapper{H}) where H<:PSY.HydroGen

    initial_time = model_initial_time(canonical)
    use_forecast_data = model_uses_forecasts(canonical)
    parameters = model_has_parameters(canonical)
    time_steps = model_time_steps(canonical)
    device_total = length(devices)
    budget_data = Vector{Tuple{String, Int64, Float64, Vector{Float64}}}(undef, device_total)

    for (ix, device) in enumerate(devices)
        bus_number = PSY.get_number(PSY.get_bus(device))
        name = PSY.get_name(device)
        tech = PSY.get_tech(device)
        # This is where you would get the water/energy storage capacity
        # which is then multiplied by the forecast value to get you the energy budget
        energy_capacity = use_forecast_data ? PSY.get_storagecapacity(device) : PSY.get_activepower(device)
        if use_forecast_data
            ts_vector = TS.values(PSY.get_data(PSY.get_forecast(PSY.Deterministic,
                                                                device,
                                                                initial_time,
                                                                "storagecapacity")))
        else
            ts_vector = ones(time_steps[end])
        end
        budget_data[ix] = (name, bus_number, energy_capacity, ts_vector)
    end
    return budget_data
end

function budget_constraints!(canonical::Canonical,
                    devices::IS.FlattenIteratorWrapper{H},
                    device_formulation::Type{<:AbstractHydroDispatchFormulation},
                    system_formulation::Type{<:PM.AbstractPowerModel}) where H<:PSY.HydroGen

    parameters = model_has_parameters(canonical)
    budget_data  = _get_budget(canonical, devices)
    if parameters
        device_budget_param_ub(canonical,
                            budget_data,
                            Symbol("budget_$(H)"), # TODO: better name for this constraint
                            UpdateRef{H}(:storagecapacity),
                            Symbol("P_$(H)"))
    else
        device_budget_param_ub(canonical,
                            budget_data,
                            Symbol("budget_$(H)"), # TODO: better name for this constraint
                            Symbol("P_$(H)"))
    end
end

function device_budget_param_ub(canonical::Canonical,
                            budget_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                            cons_name::Symbol,
                            param_reference::UpdateRef,
                            var_name::Symbol)

    time_steps = model_time_steps(canonical)
    variable = get_variable(canonical, var_name)
    set_name = (r[1] for r in budget_data)
    no_of_budgets = length(budget_data[1][4])
    time_lengths = time_steps/length(budget_data[1][4])
    time_chunks = reshape(time_steps, (time_lengths, no_of_budgets))
    constraint = _add_cons_container!(canonical, cons_name, set_name, no_of_budgets)
    param = _add_param_container!(canonical, param_reference, names, no_of_budgets)

    for data in budget_data, i in 1:no_of_budgets
        name = data[1]
        forecast = data[4][i]
        multiplier = data[3]
        param[name] = PJ.add_parameter(canonical.JuMPmodel, forecast)
        constraint[name] = JuMP.@constraint(canonical.JuMPmodel,
                    sum([variable[name, t] for t in time_chunks[:, i]]) <= multiplier*param[name])
    end

    return

end

function device_budget_ub(canonical::Canonical,
                            budget_data::Vector{Tuple{String, Int64, Float64, Vector{Float64}}},
                            cons_name::Symbol,
                            var_name::Symbol)

    time_steps = model_time_steps(canonical)
    variable = get_variable(canonical, var_name)
    set_name = (r[1] for r in budget_data)
    no_of_budgets = length(budget_data[1][4])
    time_lengths = time_steps/length(budget_data[1][4])
    time_chunks = reshape(time_steps, (time_lengths, no_of_budgets))
    constraint = _add_cons_container!(canonical, cons_name, set_name, no_of_budgets)

    for data in budget_data, i in 1:no_of_budgets
        name = data[1]
        forecast = data[4][i]
        multiplier = data[3]
        constraint[name] = JuMP.@constraint(canonical.JuMPmodel,
                    sum([variable[name, t] for t in time_chunks[:, i]]) <= multiplier*forecast)
    end

    return

end


text:
module Initializations

import PyCall: PyObject, pycall

import ..Keras
import ..Keras: PyDoc

const keras_initializer_obj = [
    "Zeros",
    "Ones",
    "Constant",
    "RandomNormal",
    "RandomUniform",
    "TruncatedNormal",
    "VarianceScaling",
    "Orthogonal",
    "Identity",

]

for i in keras_initializer_obj
    init_name = Symbol(i)

    @eval begin
        struct $init_name
            obj::PyObject

            @doc PyDoc(Keras._initializers, Symbol($i)) function $init_name(args...; kwargs...)
                new(Keras._initializers[Symbol($i)](args...; kwargs...))
            end
        end

        PyObject(initializer::$(init_name)) = initializer.obj
    end
end

const keras_initializer_funcs = [
    "lecun_uniform",
    "glorot_normal",
    "glorot_uniform",
    "he_normal",
    "he_uniform",
]

for i in keras_initializer_funcs
    init_name = Symbol(i)

    @eval begin
        @doc PyDoc(Keras._initializers, Symbol($i)) function $init_name(args...; kwargs...)
            Keras._initializers[Symbol($i)](args...; kwargs...)
        end
    end
end

end


text:
using SolverPeeker
using Test
using JuMP
using MathOptInterface
const MOI = MathOptInterface

# Model with linear constraints and objective
model1 = Model()
@variable(model1, x >= 5, Int)
@variable(model1, y <= 5)
@variable(model1, 0 <= z <= 50)
@objective(model1, Max, 3x+5y-2z+2)
@constraint(model1, x+2y <= 45)
@constraint(model1, z - y >= 4)
@constraint(model1, 0 <= 2z <= 7)
@constraint(model1, x+y+2z <= 55)

function fill_small_test_model!(model::JuMP.Model)
    # The model does not need to make sense, just use many different features.
    @variable(model, a[1:5] >= 0, Int)
    @variable(model, b[6:10], Bin)
    @variable(model, c[1:3] == 0)
    @variable(model, 10 <= d[1:3] <= 20)
    @constraint(model, con1, sum(a) + sum(b) <= 5)
    @constraint(model, con2, sum(b) >= 3)
    @constraint(model, con3, sum(d[1:2]) >= 5)
    @constraint(model, con4, sum(d) <= (sum(c) + 10))
    @objective(model, Max, sum(a) - sum(b) + sum(d))
    return model
end

@testset "lprep.jl" begin
    #lpmod = lpmodel(model1)
    #show(lpmodel)
    lprep = LPRep(model1)
    show(lprep.A)
    print("\n\n")
    #show(typeof(model1))
    #show(lprep.obj_constant)
    print("\n\n")
    #show(typeof(lpmodel.sense))
    #cis_i = MOI.get(lpmodel, MOI.ListOfConstraintIndices{MathOptInterface.SingleVariable,MathOptInterface.Integer}())
end

max_stars_repo_name: henriean/solver-peeker

max_issues_repo_name: henriean/solver-peeker

max_forks_repo_name: henriean/solver-peeker


text:
rm(joinpath(@__DIR__, "..", "Manifest.toml"), force = true) # Remove local Manifest.toml

test_env = joinpath(@__DIR__, "..", "env", "test")
push!(LOAD_PATH, test_env)
push!(LOAD_PATH, joinpath(@__DIR__, ".."))

using Pkg
Pkg.activate(test_env)
Pkg.instantiate(; verbose = true)

using Test
using PkgTemplates

@testset "Test generate" begin
    julia_version = "1.4"
    org = "charleskawczynski"
    pkg_name = "MyGeneratedPkg.jl"
    directory = joinpath(@__DIR__, "..", "..", "GeneratedPackages")

    ci = ContinuousIntegration(
        azure=Azure(;org=org, pkg_name=pkg_name, julia_version=julia_version),
        travis=Travis(;julia_version=julia_version),
        appveyor=AppVeyor(;julia_version=julia_version),
    )

    t = Template(;
        julia_version=julia_version,
        directory=directory,
        pkg_name=pkg_name,
        ci=ci,
    )

    mkpath(directory)

    generate(t)

    rm(t.pkg_root, recursive=true, force=true) # Remove local generated packages
    rm(replace(t.pkg_root, ".jl" => ""), recursive=true, force=true) # Remove local generated packages
    generate(t)

    result = open(f->read(f, String), joinpath(t.pkg_root, ".gitignore"))
    expected = PkgTemplates.gitignore()
    @test result == expected

    result = open(f->read(f, String), joinpath(t.pkg_root, "test", "runtests.jl"))
    expected = PkgTemplates.runtests(t.pkg_name_no_ext)
    @test result == expected

    result = open(f->read(f, String), joinpath(t.pkg_root, "bors.toml"))
    expected = PkgTemplates.toml(t.bors)
    @test result == expected
end

max_stars_repo_name: charleskawczynski/PkgTemplates

max_issues_repo_name: charleskawczynski/PkgTemplates

max_forks_repo_name: charleskawczynski/PkgTemplates


text:
# Stream management

export
    CuStream, default_stream, synchronize, destroy

immutable CuStream
    handle::Ptr{Void}
    blocking::Bool
    priority::Int
end

default_stream() = CuStream(convert(Ptr{Void}, 0), true, 0)

synchronize(s::CuStream) = @cucall(:cuStreamSynchronize, (Ptr{Void},), s.handle)

destroy(s::CuStream) = @cucall(:cuStreamDestroy, (Ptr{Void},), s.handle)


text:
using TupleTools

export NotImplementedError, LocationConflictError, QubitMismatchError

# NOTE: kwargs do not involve in multiple dispatch
#       no need to store kwargs
struct NotImplementedError{ArgsT} <: Exception
    name::Symbol
    args::ArgsT
end

NotImplementedError(name::Symbol) = NotImplementedError(name, ())

function Base.show(io::IO, e::NotImplementedError)
    str = join(map(typeof, e.args), ", ::")
    str = "::" * str
    print(io, "NotImplementedError: $(e.name) is not implemented for (", str, "), ",
              "please implement this method for your custom type")
end

function Base.show(io::IO, e::NotImplementedError{Tuple{}})
    print(io, "$(e.name) is not implemented.")
end

struct LocationConflictError <: Exception
    msg::String
end

function Base.show(io::IO, e::LocationConflictError)
    print(io, "LocationConflictError: ", e.msg)
    # print(io, "locations of $(e.blk1) and $(e.blk2) is conflict.")
end

# NOTE: More detailed error msg?
"""
    QubitMismatchError <: Exception

Qubit number mismatch error when applying a Block to a Register or concatenating Blocks.
"""
struct QubitMismatchError <: Exception
    msg::String
end

function show(io::IO, e::QubitMismatchError)
    print(io, e.msg)
end

export islocs_inbounds, islocs_conflict

_sort(x::Vector; by=identity) = sort(x, by=by)
_sort(x::Tuple; by=identity) = TupleTools.sort(x, by=by)

# NOTE: this method assumes its input is not empty, it gets rid of errors
nonempty_minimum(x::UnitRange) = x.start
nonempty_minimum(x::Integer) = x
nonempty_maximum(x::UnitRange) = x.stop
nonempty_maximum(x::Integer) = x

const AddressVector{T} = Vector{T} where {T <: Union{Integer, UnitRange}}
const AddressNTuple{N, T} = NTuple{N, T} where {T <: Union{Integer, UnitRange}}
const AddressList{T} = Union{AddressVector{T}, AddressNTuple{N, T} where N}

"""
    islocs_inbounds(n, locs) -> Bool

Check if the input locations are inside given bounds `n`.
"""
function islocs_inbounds(n::Int, locs::AddressList)
    length(locs) == 0 && return true
    locs = _sort(locs, by=x->nonempty_minimum(x))
    (minimum(first(locs)) > 0 && maximum(last(locs)) <= n) || return false
    return true
end

"""
    islocs_conflict(locs) -> Bool

Check if the input locations has conflicts.
"""
function islocs_conflict(locs::AddressList)
    locs = _sort(locs)
    for (nxt, cur) in zip(locs[2:end], locs[1:end-1])
        nonempty_minimum(nxt) > nonempty_maximum(cur) || return true
    end
    return false
end

function process_msgs(msgs...; default="")
    msg = isempty(msgs) ? default : msgs[1]
    if isa(msg, AbstractString)
        msg = msg # pass-through
    elseif !isempty(msgs) && (isa(msg, Expr) || isa(msg, Symbol))
        # message is an expression needing evaluating
        msg = :(Main.Base.string($(esc(msg))))
    elseif isdefined(Main, :Base) && isdefined(Main.Base, :string) && applicable(Main.Base.string, msg)
        msg = Main.Base.string(msg)
    else
        # string() might not be defined during bootstrap
        msg = :(Main.Base.string($(Expr(:quote,msg))))
    end
    return msg
end

# NOTE: we may use @assert in the future
#       these macro will help us keep original APIs
export @assert_locs_safe, @assert_locs_inbounds

"""
    @assert_locs_inbounds <number of total qubits> <locations list> [<msg>]

Assert if all the locations are inbounds.
"""
macro assert_locs_inbounds(n, locs, msgs...)
    msg = process_msgs(msgs...; default="locations is out of bounds!")

    return quote
        islocs_inbounds($(esc(n)), $(esc(locs))) || error($msg)
        nothing
    end
end

"""
    @assert_locs_safe <number of total qubits> <locations list> [<msg>]

Assert if all the locations are:
    - inbounds.
    - do not have any conflict.
"""
macro assert_locs_safe(n, locs, msgs...)
    msg = process_msgs(msgs...; default="locations conflict.")
    return quote
        @assert_locs_inbounds $(esc(n)) $(esc(locs))
        islocs_conflict($(esc(locs))) && throw(LocationConflictError($msg))
        nothing
    end
end


text:
# ```math
# x1[a_, b_, c_, t_] :=  Sin[.5 t] - a Sin[b t]*Cos[t] - .1 c Sin[10 b t];
# y1[a_, b_, c_, t_] :=  Cos[.5 t] - a Sin[b t]*Sin[t] - .1 c Cos[10 b t];
#
# GraphicsGrid[Partition[ParallelTable[
#    With[{
#      a = RandomReal[{-1.5, 1.5}],
#      b = RandomInteger[{3, 15}],
#      c = RandomReal[{0, 1.5}],
#      clr1 = Black,
#      clr2 = RGBColor @@ RandomReal[1, 3],
#      clr3 = RGBColor @@ RandomReal[1, 3],
#      thick = RandomReal[{.04, .5}],
#      tm = 1},
#     ParametricPlot[
#      Evaluate[{{x1[a, b, c, t], y1[a, b, c, t]}, {x1[a, b, c, t],
#         y1[a, b, c, t]}}], {t, 0, tm 4 \[Pi]},
#      PlotStyle -> {{clr2, Thickness[0.001` + 0.05` thick]}, {clr3,
#         Thickness[0.001` + 0.01` thick]}}, Axes -> False,
#      PlotPoints -> 200, PlotRange -> All, Background -> clr1]]
#    , {n, 32}], 4], ImageSize -> 600]
# ```
using Distributions
using OdinSon
using PyPlot

function colorflake()
    clr1 = RGB(rand(3)...)
    clr2 = RGB(rand(3)...)
    thick = rand(Uniform(2, 7)) # need to get this better, MM linewidth is based on size of figure
    tm = 1 # this seems needless in the current code

    a = rand(Uniform(0, 1.5))
    b = rand(3:15)
    c = rand(Uniform(0, 1.5))
    # this creates a closure over the parameters, might not be the best way
    # this is also vectorized since sin/cos are, can't wait for f.() syntax!
    x1(t) = sin(0.5*t) - a*sin(b*t).*cos(t) - 0.1*c*sin(10*b*t)
    y1(t) = cos(0.5*t) - a*sin(b*t).*sin(t) - 0.1*c*cos(10*b*t)

    t = linspace(0, tm*4.0*π, 5000)
    plot(x1(t), y1(t), c=clr1, linewidth=thick)
    plot(x1(t), y1(t), c=clr2, linewidth=0.2*thick)
    gca()[:set_aspect]("equal")
    gca()[:set_axis_bgcolor](NC"black")
end

colorflake()

max_stars_repo_name: gabrielgellner/OdinSon

max_issues_repo_name: gabrielgellner/OdinSon

max_forks_repo_name: gabrielgellner/OdinSon


text:
function tokenize(tu::LibClang.CXTranslationUnit, cursor::LibClang.CXCursor)
	tokens = Ref{Ptr{LibClang.CXToken}}(C_NULL)
	count = Ref{Cuint}(0)
	try
		cxrange = LibClang.clang_getCursorExtent(cursor)
		LibClang.clang_tokenize(tu, cxrange, tokens, count)

		return map(1:count[]) do i
			token = unsafe_load(tokens[], i)
			kind = LibClang.clang_getTokenKind(token)
			text = string(tu, token)
			return kind => text
		end
	finally
		tokens[] == C_NULL || LibClang.clang_disposeTokens(tu, tokens[], count[])
	end
end

function convert_macro(tu::LibClang.CXTranslationUnit, cursor::LibClang.CXCursor, exports::Vector{String})
	name = convert_name(cursor)
	expr = nothing
	deps = Set{String}()
	tokens = []
	try
		Bool(LibClang.clang_Cursor_isMacroFunctionLike(cursor)) && error("Function-like macros are not yet supported")

		tokens = tokenize(tu, cursor)
		e = convert_macro_def(deps, tokens, 2:length(tokens))  # NOTE: first token is the macro identifier itself, so skip it
		if !isnothing(e)
			Meta.parse(e)
			expr = e
		end
	catch
		tokens = isempty(tokens) ? "" : "\n"*join(map(last, tokens), ' ')
		@warn "Unable to convert macro `$(name)` near $(string(CodeLocation(cursor)))$(tokens)"
	end
	isnothing(expr) && return nothing

	expr = "macro $(name)() return quote $(expr) end end"
	comments = Dict{String, Comment}(convert_comment(cursor, "@"*name))

	if isempty(deps) && !(name in exports)  # if macro has no deps and not yet exported, declare a const
		expr = "$(expr)\nconst $(name) = (@$(name))"
		push!(comments, convert_comment(cursor, name))
	end

	return name => Converted(expr, comments)
end

function convert_macro_def(deps::Set{String}, tokens::Vector{<:Pair}, range::UnitRange)
	if length(range) <= 0
		return nothing
	elseif length(range) == 1 && tokens[range.start].first == LibClang.CXToken_Identifier
		sym = tokens[range.start].second
		push!(deps, sym)
		return "(@$(sym))"
	elseif length(range) == 1 && tokens[range.start].first == LibClang.CXToken_Literal
		return convert_literal(tokens[range.start].second)
	elseif length(range) >= 1 && all(i -> tokens[i].first == LibClang.CXToken_Keyword, range)
		return convert_keyword(join(map(i -> tokens[i].second, range), ' '))
	elseif length(range) >= 2 && tokens[range.start].first == LibClang.CXToken_Punctuation && tokens[range.start].second in ("+", "-", "~")
		return convert_operator(tokens[range.start].second)*convert_macro_def(deps, tokens, (range.start+1):range.stop)
	elseif tokens[range.start] in (
			(LibClang.CXToken_Punctuation => "("),
			(LibClang.CXToken_Punctuation => "{"),
			(LibClang.CXToken_Punctuation => "["),
	)
		(pre, post) = if tokens[range.start] == (LibClang.CXToken_Punctuation => "(")
			("(", ")")
		elseif tokens[range.start] == (LibClang.CXToken_Punctuation => "{")
			("begin", "end")
		elseif tokens[range.start] == (LibClang.CXToken_Punctuation => "[")
			("[(", ") + 1]")
		end

		sub = find_subrange(tokens, range)
		cvt = pre*convert_macro_def(deps, tokens, (sub.start+1):(sub.stop-1))*post

		return sub.stop == range.stop ? cvt : cvt*" "*convert_macro_def(deps, tokens, (sub.stop+1):range.stop)
	end

	error("Unhandled macro conversion:  $(join(map(last, tokens[range]), ' '))")
end

function find_subrange(tokens::Vector{<:Pair}, range::UnitRange)
	inc = tokens[range.start]
	dec = if inc == (LibClang.CXToken_Punctuation => "(")
		LibClang.CXToken_Punctuation => ")"
	elseif inc == (LibClang.CXToken_Punctuation => "{")
		LibClang.CXToken_Punctuation => "}"
	elseif inc == (LibClang.CXToken_Punctuation => "[")
		LibClang.CXToken_Punctuation => "]"
	end

	depth = 0
	stop = 0
	for ind in range.start:range.stop
		stop = ind
		tokens[stop] == inc && (depth += 1)
		tokens[stop] == dec && (depth -= 1)
		depth == 0 && break
	end
	depth == 0 || error("Failed to find subrange $(join(map(last, tokens[range]), ' '))")
	return range.start:stop
end

function convert_keyword(keywords::String)
	if keywords in ("char", "signed char")
		return "\$(𝐣𝐥.Cchar)"
	elseif keywords in ("int",)
		return "\$(𝐣𝐥.Cint)"
	elseif keywords in ("short", "short int")
		return "\$(𝐣𝐥.Cshort)"
	elseif keywords in ("long", "long int",)
		return "\$(𝐣𝐥.Clong)"
	elseif keywords in ("long long", "long long int")
		return "\$(𝐣𝐥.Clonglong)"
	elseif keywords in ("unsigned char",)
		return "\$(𝐣𝐥.Cuchar)"
	elseif keywords in ("unsigned int",)
		return "\$(𝐣𝐥.Cuint)"
	elseif keywords in ("unsigned short", "unsigned short int")
		return "\$(𝐣𝐥.Cushort)"
	elseif keywords in ("unsigned long", "unsigned long int")
		return "\$(𝐣𝐥.Culong)"
	elseif keywords in ("unsigned long long", "unsigned long long int")
		return "\$(𝐣𝐥.Culonglong)"
	elseif keywords in ("extern",)
		return nothing
	end

	error("Unhandled keywords:  $(keywords)")
end

function convert_literal(literal::String)
	m = match(r"^'(.*)'$", literal)
	!isnothing(m) && return repr(Meta.parse(m.match))

	m = match(r"^\"(.*)\"$", literal)
	!isnothing(m) && return repr(Meta.parse(m.match))

	m = match(r"^(0|0x|0b)?([\da-f]+)([ul]*)$", lowercase(literal))
	if !isnothing(m)
		(pre, val, suf) = m.captures
		isUnsigned = occursin('u', suf)
		howLong = count(isequal('l'), suf)

		typ = howLong == 0 ? Cint : (howLong == 1 ? Clong : Clonglong)
		typ = isUnsigned || !isnothing(pre) ? unsigned(typ) : typ
		pre = isnothing(pre) ? "" : (pre == "0" ? "0o" : pre)

		return repr(parse(typ, pre*val))
	end

	m = match(r"^(\d*\.)?(\d*)(e[-+]?\d+)?([fl])?$", lowercase(literal))
	if !isnothing(m)
		(val1, val2, exp, suf) = m.captures
		val = val1*val2
		val = val[end] == '.' ? val*"0" : val
		val = val[1] == '.' ? "0"*val : val

		suf == "l" && error("Unable to convert `long double` literal:  $(literal)")
		typ = suf == "f" ? Cfloat : Cdouble
		exp = isnothing(exp) ? "" : exp

		return repr(parse(typ, val*exp))
	end

	error("Failed to convert literal:  $(literal)")
end

function convert_operator(op::String)
	return op == "^" ? "⊻" : op
end


text:
import ConvexHull: Point, grahamscan, isrightturn
using Test

@test Point(4.2, 42.0).x == 4.2
@test Point(4.2, 42.0).y == 42.0

@test (Point(1.0, 1.0) < Point(1.0, 1.0)) == false
@test (Point(1.0, 1.0) < Point(1.0, 2.0)) == true
@test (Point(1.0, 0.0) < Point(1.0, 1.0)) == true
@test (Point(1.0, 0.0) < Point(0.0, 1.0)) == false
@test (Point(1.0, 2.0) < Point(2.0, 1.0)) == true

@test isrightturn(Point(0, 0), Point(0, 1), Point(1, 1)) == true
@test isrightturn(Point(0, 0), Point(1, 1), Point(2, 1)) == true
@test isrightturn(Point(0, 0), Point(1, 1), Point(2, 2)) == false
@test isrightturn(Point(0, 0), Point(1, 1), Point(1, 2)) == false
@test isrightturn(Point(0, 0), Point(-1, -1), Point(-2, -2)) == false
@test isrightturn(Point(0, 0), Point(-1, -1), Point(-1, -2)) == false
@test isrightturn(Point(0, 0), Point(-1, -1), Point(-2, -1)) == true

@test grahamscan([Point(1, 1), Point(-1, -1), Point(-1, 1), Point(1, -1)]) ==
                 [Point(-1, -1), Point(-1, 1), Point(1, 1), Point(1, -1)]

@test grahamscan([Point(-1, 0), Point(1, 0), Point(0, 1)]) ==
                 [Point(-1, 0), Point(0, 1), Point(1, 0)]

@test grahamscan([Point(-1, 0), Point(-1, 0), Point(1, 0), Point(0, 1)]) ==
                 [Point(-1, 0), Point(0, 1), Point(1, 0)]


text:
include("../imports.jl")

@testset ExtendedTestSet "`ParamGroups`" begin
    model = Chain(Dense(3, 5), Dense(5, 3))
    paramgroups = ParamGroups(IndexGrouper([1, 2]), model)

    @test getgroup(paramgroups, model[1].weight) == 1
    @test getgroup(paramgroups, model[2].weight) == 2
    @test getgroup(paramgroups, rand(10)) === nothing
end


text:
# we write a function that uses Monte Carlo to estimate the volume of the unit hypersphere in d dimensions using n trials
function sphereVolMC(d, n)

	# we initialize the counter variable
	count = 0

	# we perform the Monte Carlo trial n times
	for i = 1:n

		# we generate a vector with d random floats from [0,1]
		x = rand(d)

		# we calculate the squared norm of x
		norm2 = x' * x

		# if x is inside the hypersphere we increment the count
		if norm2 < 1
			count += 1
		end
	end

	# we calculate the estimated volume given by
	vol = 2^d * count / n

	# we return the estimated volume
	return vol
end


text:
# Gaussian drift proposal
@gen function gaussian_drift(tr::Trace, params=Dict(:β => (0.3, 0.01, Inf)))
    for (addr, distargs) in params
        if distargs isa Tuple && length(distargs) == 3
            sigma, lb, ub = distargs
            {addr} ~ trunc_normal(tr[addr], sigma, lb, ub)()
        else
            sigma = distargs
            {addr} ~ normal(tr[addr], sigma)()
        end
    end
end

# Guess beta from the data from initial infection counts
@gen function init_beta(obs::ChoiceMap, t::Int)
    params = create_params()
    cities = 1:length(params[:Ns])
    infected_0 = sum(params[:Is])
    infected_t = sum(obs[:step => t => :obs => :infected => c] for c in cities)
    β_hat = (infected_t/infected_0) ^ (1/t) - 1
    β ~ trunc_normal(β_hat, 0.1, 0, Inf)()
end


text:
using GenomePermutations
using Documenter

cd("docs")
DocMeta.setdocmeta!(
    GenomePermutations, :DocTestSetup,
    :(using GenomePermutations);
    recursive=true)
# REVISE REGEX
makedocs(sitename="GenomePermutations.jl", modules = [GenomePermutations])


text:
# Quadrotor in Maze

# model
model = Dynamics.quadrotor
# model = Dynamics.quadrotor_euler
model_d = rk3(model)
n = model.n; m = model.m
q0 = [1.;0.;0.;0.] # unit quaternion

x0 = zeros(T,n)
x0[1:3] = [0.; 0.; 10.]
x0[4:7] = q0

xf = zero(x0)
xf[1:3] = [0.;60.; 10.]
xf[4:7] = q0;

# cost
Q = (1.0e-3)*Diagonal(I,n)
Q[4:7,4:7] = (1.0e-3)*Diagonal(I,4)
R = (1.0e-2)*Diagonal(I,m)
Qf = 1.0*Diagonal(I,n)

u_min = 0.
u_max = 50.
x_max = Inf*ones(model.n)
x_min = -Inf*ones(model.n)

x_max[1:3] = [25.0; Inf; 20]
x_min[1:3] = [-25.0; -Inf; 0.]
bnd_u = BoundConstraint(n,m,u_min=u_min, u_max=u_max)
bnd = BoundConstraint(n,m,u_min=u_min, u_max=u_max, x_min=x_min, x_max=x_max)

xf_no_quat_U = copy(xf)
xf_no_quat_L = copy(xf)
xf_no_quat_U[4:7] .= Inf
xf_no_quat_L[4:7] .= -Inf
xf_no_quat_U[8:10] .= 0.
xf_no_quat_L[8:10] .= 0.
bnd_xf = BoundConstraint(n,m,x_min=xf_no_quat_L,x_max=xf_no_quat_U)
# goal = goal_constraint(xf)

N = 101 # number of knot points
tf = 5.0
dt = tf/(N-1) # total time

U_hover = [0.5*9.81/4.0*ones(m) for k = 1:N-1] # initial hovering control trajectory
obj = LQRObjective(Q, R, Qf, xf, N) # objective with same stagewise costs

# Constraints
cyl_obs = [(0., 10., 3.),
          (10., 30., 3.),
          (-13., 25., 2.),
          (5.,50.,4.)]

sph_obs = [( 0., 40., 5., 2.),
           (-5., 15., 3., 1.),
           (10., 20., 7., 2.)]
r_quad = 2.0

cylinder_constraint = let cylinders=cyl_obs, spheres=sph_obs, r_quad=r_quad
    function cylinder_constraint(c,x,u)
        for (i,cyl) in enumerate(cylinders)
            c[i] = circle_constraint(x,cyl[1],cyl[2],cyl[3]+r_quad)
        end
    end
end

sphere_constraint = let cylinders=cyl_obs, spheres=sph_obs, r_quad=r_quad
    function sphere_constraint(c,x,u)
        for (i,sphere) in enumerate(spheres)
            c[i] = TrajectoryOptimization.sphere_constraint(x,sphere[1],sphere[2],sphere[3],sphere[3]+r_quad)
        end
    end
end

cyl_con = Constraint{Inequality}(cylinder_constraint,n,m,length(cyl_obs),:cylinders)
sph_con = Constraint{Inequality}(  sphere_constraint,n,m,length(sph_obs),:spheres)

constraints = Constraints(N)
constraints[1] += bnd_u
for k = 2:N-1
    constraints[k] += bnd + cyl_con + sph_con
end
constraints[N] += bnd_xf

quad_obs = Problem(model_d, obj, constraints=constraints, x0=x0, xf=xf, N=N, dt=dt)
initial_controls!(quad_obs,U_hover); # initialize problem with control


text:
@inline ArrayInterface.size(A::PtrArray) = A.size
@inline ArrayInterface.size(A::StrideArray) = A.ptr.size
@inline VectorizationBase.bytestrides(A::StridedPointer) = A.strd
@inline VectorizationBase.bytestrides(A::PtrArray) = A.ptr.strd
@inline VectorizationBase.bytestrides(A::StrideArray) = A.ptr.ptr.strd

@inline bytestride(A, n) = VectorizationBase.bytestrides(A)[n]

@generated function ArrayInterface.strides(A::PtrArray{S,D,T,N}) where {S,D,T,N}
    shifter = static_expr(VectorizationBase.intlog2(sizeof(T)))
    x = Expr(:tuple)
    for n in 1:N
        push!(x.args, Expr(:call, :(>>>), Expr(:ref, :x, n), shifter))
    end
    quote
        $(Expr(:meta,:inline))
        x = A.ptr.strd
        $x
    end
end
@inline ArrayInterface.strides(A::StrideArray) = strides(A.ptr)

@inline Base.size(A::AbstractStrideArray) = map(Int, ArrayInterface.size(A))
@inline Base.strides(A::AbstractStrideArray) = map(Int, ArrayInterface.strides(A))

@inline create_axis(s, ::Zero) = CloseOpen(s)
@inline create_axis(s, ::One) = One():s
@inline create_axis(s, o) = CloseOpen(s, s+o)

@inline ArrayInterface.axes(A::AbstractStrideArray) = map(create_axis, size(A), offsets(A))
@inline Base.axes(A::AbstractStrideArray) = axes(A)

@inline ArrayInterface.offsets(A::PtrArray) = A.ptr.offsets
@inline ArrayInterface.offsets(A::StrideArray) = A.ptr.ptr.offsets

@inline ArrayInterface.static_length(A::AbstractStrideArray) = prod(size(A))

# type stable, because index known at compile time
@inline type_stable_select(t::NTuple, ::StaticInt{N}) where {N} = t[N]
@inline type_stable_select(t::Tuple, ::StaticInt{N}) where {N} = t[N]
# type stable, because tuple is homogenous
@inline type_stable_select(t::NTuple, i::Integer) = t[i]
# make the tuple homogenous before indexing
@inline type_stable_select(t::Tuple, i::Integer) = map(Int, t)[i]

@inline function Base.axes(A::AbstractStrideVector, i::Integer)
    if i == 1
        o = type_stable_select(offsets(A), i)
        s = type_stable_select(size(A), i)
        return create_axis(s, o)
    else
        return One():1
    end
end
@inline function Base.axes(A::AbstractStrideArray, i::Integer)
    o = type_stable_select(offsets(A), i)
    s = type_stable_select(size(A), i)
    create_axis(s, o)
end

@inline zeroindex(r::ArrayInterface.OptionallyStaticUnitRange{One}) = CloseOpen(Zero(), last(r))
@inline zeroindex(r::Base.OneTo) = CloseOpen(Zero(), last(r))
@inline zeroindex(r::AbstractUnitRange) = Zero():(last(r)-first(r))

@inline zeroindex(r::CloseOpen{Zero}) = r
@inline zeroindex(r::ArrayInterface.OptionallyStaticUnitRange{Zero}) = r


text:
__precompile__()

module CircuitSimulator

using ResumableFunctions
using Base.Meta
using Base.CoreLogging: Debug, global_logger
using Logging: ConsoleLogger
global_logger(ConsoleLogger(stderr, Debug))

include("def_macro.jl")
include("types.jl")
include("spiceunits.jl")
include("parse_netlist.jl")
include("readcards.jl")
include("modified_nodal_analysis.jl")

end # module


text:
import Tables: rows, columns, istable, rowaccess, columnaccess, schema, Schema

import Base: propertynames, getproperty, getindex

export istable, columns

istable(::Type{JDFFile}) = true
istable(::JDFFile) = true

rowaccess(::JDFFile) = false
columnaccess(::JDFFile) = true

rowaccess(::Type{<:JDFFile}) = false
columnaccess(::Type{<:JDFFile}) = true

propertynames(jdf::JDFFile) = names(jdf)

getproperty(jdf::JDFFile, col::Symbol) =  JDF.load(jdf; cols = [col]).columns[col]

schema(jdf::JDFFile) = begin
    meta = metadata(jdf)
    Schema(meta.names, map(x -> x.type, meta.metadatas))
end

columns(jdf::JDFFile) = jdf

# this is the table type specific to JDF
struct Table
    columns::NamedTuple
end

nrow(t::Table) = length(t.columns[1])

ncol(t::Table) = length(t.columns)

Tables.columns(t::Table) = t.columns

Tables.istable(t::Table) = true

function Base.getindex(t::Table, col::Symbol)
    t.columns[col]
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF