text:
export Dense,
       DenseTensor,
       array,
       matrix,
       vector,
       contract,
       outer,
       read,
       permutedims!!,
       write,
       ⊗

#
# Dense storage
#

struct Dense{ElT, VecT<:AbstractVector} <: TensorStorage{ElT}
  data::VecT
  function Dense{ElT, VecT}(data::AbstractVector) where {ElT,
                                                         VecT<:AbstractVector{ElT}}
    return new{ElT, VecT}(data)
  end
end

function Dense(data::VecT) where {VecT <: AbstractVector{ElT}} where {ElT}
  return Dense{ElT,VecT}(data)
end

function Dense{ElR}(data::AbstractVector{ElT}) where {ElR,ElT}
  ElT == ElR ? Dense(data) : Dense(ElR.(data))
end

# Construct from a set of indices
function Dense{ElT, VecT}(inds) where {ElT,
                                       VecT <: AbstractVector{ElT}}
  return Dense(VecT(dim(inds)))
end

Dense{ElT}(dim::Integer) where {ElT} = Dense(zeros(ElT, dim))

Dense{ElT}(::UndefInitializer,
           dim::Integer) where {ElT} = Dense(Vector{ElT}(undef,dim))

Dense(::Type{ElT},
      dim::Integer) where {ElT} = Dense{ElT}(dim)

Dense(x::ElT,
      dim::Integer) where {ElT<:Number} = Dense(fill(x,dim))

Dense(dim::Integer) = Dense(Float64,dim)

Dense(::Type{ElT},
      ::UndefInitializer,
      dim::Integer) where {ElT} = Dense{ElT}(undef,dim)

Dense(::UndefInitializer,dim::Integer) = Dense(Float64,undef,dim)

Dense{ElT}() where {ElT} = Dense(ElT[])
Dense(::Type{ElT}) where {ElT} = Dense{ElT}()

#
# Random constructors
#

function Base.randn(::Type{StoreT},
                    dim::Integer) where {StoreT<:Dense}
  return Dense(randn(eltype(StoreT), dim))
end

Base.copy(D::Dense) = Dense(copy(data(D)))

function Base.complex(::Type{Dense{ElT,
                             Vector{ElT}}}) where {ElT}
  return Dense{complex(ElT),Vector{complex(ElT)}}
end

Base.similar(D::Dense) = Dense(similar(data(D)))

Base.similar(D::Dense,
             length::Int) = Dense(similar(data(D), length))

function Base.similar(::Type{<: Dense{ElT, VecT}},
                      length::Int) where {ElT,
                                          VecT <: AbstractVector{ElT}}
  return Dense(similar(Vector{ElT}, length))
end

Base.similar(D::Dense,
             ::Type{T}) where {T<:Number} = Dense(similar(data(D),T))

Base.zeros(DenseT::Type{<: Dense{ElT}},
           inds) where {ElT} = zeros(DenseT, dim(inds))

# TODO: should this do something different for SubArray?
Base.zeros(::Type{<: Dense{ElT}},
           dim::Int) where {ElT} = Dense{ElT}(dim)

function Base.promote_rule(::Type{<:Dense{ElT1,VecT1}},
                           ::Type{<:Dense{ElT2,VecT2}}) where {ElT1,VecT1,
                                                               ElT2,VecT2}
  ElR = promote_type(ElT1,ElT2)
  VecR = promote_type(VecT1,VecT2)
  return Dense{ElR,VecR}
end

# This is to get around the issue in Julia that:
# promote_type(Vector{ComplexF32},Vector{Float64}) == Vector{T} where T
function Base.promote_rule(::Type{<:Dense{ElT1,Vector{ElT1}}},
                           ::Type{<:Dense{ElT2,Vector{ElT2}}}) where {ElT1,ElT2}
  ElR = promote_type(ElT1,ElT2)
  VecR = Vector{ElR}
  return Dense{ElR,VecR}
end

# This is for type promotion for Scalar*Dense
function Base.promote_rule(::Type{<:Dense{ElT1,Vector{ElT1}}},
                           ::Type{ElT2}) where {ElT1,
                                                ElT2<:Number}
  ElR = promote_type(ElT1,ElT2)
  VecR = Vector{ElR}
  return Dense{ElR,VecR}
end

Base.convert(::Type{<:Dense{ElR,VecR}},
             D::Dense) where {ElR,VecR} = Dense(convert(VecR,data(D)))

# Make generic to all storage types?
Base.:*(D::Dense,x::Number) = Dense(x*data(D))
Base.:*(x::Number,D::Dense) = D*x

#
# DenseTensor (Tensor using Dense storage)
#

const DenseTensor{ElT,N,StoreT,IndsT} = Tensor{ElT,N,StoreT,IndsT} where {StoreT<:Dense}

DenseTensor(::Type{ElT},
            inds::Dims) where {ElT} = tensor(Dense(ElT,dim(inds)),inds)

DenseTensor(::Type{ElT},
            inds::Int...) where {ElT} = DenseTensor(ElT,inds)

DenseTensor(inds::Dims) = tensor(Dense(dim(inds)),inds)

DenseTensor(inds::Int...) = DenseTensor(inds)

DenseTensor(::Type{ElT},
            ::UndefInitializer,
            inds) where {ElT} = tensor(Dense(ElT,undef,dim(inds)),inds)

DenseTensor(::Type{ElT},
            ::UndefInitializer,
            inds::Int...) where {ElT} = DenseTensor(ElT,undef,inds)

DenseTensor(::UndefInitializer,
            inds::Dims) = tensor(Dense(undef,dim(inds)),inds)

DenseTensor(::UndefInitializer,
            inds::Int...) = DenseTensor(undef,inds)

# For convenience, direct Tensor constructors default to Dense
Tensor(::Type{ElT},inds::Dims) where {ElT} = DenseTensor(ElT,inds)
Tensor(::Type{ElT},inds::Int...) where {ElT} = DenseTensor(ElT,inds...)
Tensor(inds::Dims) = DenseTensor(inds)
Tensor(inds::Int...) = DenseTensor(inds...)

Tensor(::Type{ElT},
       ::UndefInitializer,
       inds::Dims) where {ElT} = DenseTensor(ElT,undef,inds)
Tensor(::Type{ElT},
       ::UndefInitializer,
       inds::Int...) where {ElT} = DenseTensor(ElT,undef,inds...)
Tensor(::UndefInitializer,
       inds::Dims) = DenseTensor(undef,inds)
Tensor(::UndefInitializer,
       inds::Int...) = DenseTensor(undef,inds...)

Tensor(A::Array{<:Number,N},
       inds::Dims{N}) where {N} = tensor(Dense(vec(A)),inds)

#
# Random constructors
#

function randomDenseTensor(::Type{ElT},
                           inds::Dims) where {ElT}
  return tensor(randn(Dense{ElT},dim(inds)),inds)
end

function randomDenseTensor(::Type{ElT},
                           inds::Int...) where {ElT}
  return randomDenseTensor(ElT, inds)
end

randomDenseTensor(inds::Dims) = randomDenseTensor(Float64, inds)

randomDenseTensor(inds::Int...) = randomDenseTensor(Float64, inds)

function randomTensor(::Type{ElT},
                      inds::Dims) where {ElT}
  return randomDenseTensor(ElT,inds)
end

function randomTensor(::Type{ElT},
                      inds::Int...) where {ElT}
  return randomDenseTensor(ElT, inds...)
end

randomTensor(inds::Dims) = randomDenseTensor(Float64, inds)

randomTensor(inds::Int...) = randomDenseTensor(Float64, inds)

# Basic functionality for AbstractArray interface
Base.IndexStyle(::Type{<:DenseTensor}) = IndexLinear()

# Override CartesianIndices iteration to iterate
# linearly through the Dense storage (faster)
Base.iterate(T::DenseTensor,args...) = iterate(store(T),args...)

function _zeros(TensorT::Type{<: DenseTensor},
                inds)
  return tensor(zeros(storetype(TensorT), dim(inds)), inds)
end

function Base.zeros(TensorT::Type{<: DenseTensor},
                    inds)
  return _zeros(TensorT, inds)
end

# To fix method ambiguity with zeros(::Type, ::Tuple)
function Base.zeros(TensorT::Type{<: DenseTensor},
                    inds::Dims)
  return _zeros(TensorT, inds)
end

function Base.zeros(TensorT::Type{<: DenseTensor},
                    inds::Tuple{})
  return _zeros(TensorT, inds)
end

function _similar(TensorT::Type{<: DenseTensor},
                  inds)
  return tensor(similar(storetype(TensorT), dim(inds)), inds)
end

function Base.similar(TensorT::Type{<: DenseTensor},
                      inds)
  return _similar(TensorT, inds)
end

# To fix method ambiguity with similar(::AbstractArray,::Tuple)
function Base.similar(TensorT::Type{<: DenseTensor},
                      inds::Dims)
  return _similar(TensorT, inds)
end

# To fix method ambiguity with similar(::AbstractArray,::Type)
function Base.similar(T::DenseTensor,
                      ::Type{ElT}) where {ElT}
  return tensor(similar(store(T), ElT), inds(T))
end

_similar(T::DenseTensor, inds) = similar(typeof(T), inds)

Base.similar(T::DenseTensor, inds) = _similar(T, inds)

# To fix method ambiguity with similar(::AbstractArray,::Tuple)
Base.similar(T::DenseTensor,inds::Dims) = _similar(T, inds)

#
# Single index
#

Base.@propagate_inbounds function _getindex(T::DenseTensor{<:Number,N},
                                            I::CartesianIndex{N}) where {N}
  return store(T)[@inbounds LinearIndices(T)[CartesianIndex(I)]]
end

Base.@propagate_inbounds function Base.getindex(T::DenseTensor{<:Number,N},
                                                I::Vararg{Int,N}) where {N}
  return _getindex(T,CartesianIndex(I))
end

Base.@propagate_inbounds function Base.getindex(T::DenseTensor{<:Number,N},
                                                I::CartesianIndex{N}) where {N}
  return _getindex(T,I)
end

Base.@propagate_inbounds function _setindex!(T::DenseTensor{<:Number,N},
                                             x::Number,
                                             I::CartesianIndex{N}) where {N}
  store(T)[@inbounds LinearIndices(T)[CartesianIndex(I)]] = x
  return T
end

Base.@propagate_inbounds function Base.setindex!(T::DenseTensor{<:Number,N},
                                                 x::Number,
                                                 I::Vararg{Int,N}) where {N}
  _setindex!(T,x,CartesianIndex(I))
  return T
end

Base.@propagate_inbounds function Base.setindex!(T::DenseTensor{<:Number,N},
                                                 x::Number,
                                                 I::CartesianIndex{N}) where {N}
  _setindex!(T,x,I)
  return T
end

#
# Linear indexing
#

Base.@propagate_inbounds Base.getindex(T::DenseTensor,i::Int) = store(T)[i]

Base.@propagate_inbounds Base.setindex!(T::DenseTensor,v,i::Int) = (store(T)[i] = v; T)

#
# Slicing
# TODO: this doesn't allow colon right now
# Create a DenseView that stores a Dense and an offset?
#

Base.@propagate_inbounds function _getindex(T::DenseTensor{ElT,N},
                                            I::CartesianIndices{N}) where {ElT,N}
  storeR = Dense(vec(@view array(T)[I]))
  indsR = Tuple(I[end]-I[1]+CartesianIndex(ntuple(_->1,Val(N))))
  return tensor(storeR, indsR)
end

Base.@propagate_inbounds function Base.getindex(T::DenseTensor{ElT,N},
                                                I...) where {ElT,N}
  return _getindex(T,CartesianIndices(I))
end

# Reshape a DenseTensor using the specified dimensions
# This returns a view into the same Tensor data
function Base.reshape(T::DenseTensor, dims)
  dim(T)==dim(dims) || error("Total new dimension must be the same as the old dimension")
  return tensor(store(T), dims)
end

# This version fixes method ambiguity with AbstractArray reshape
function Base.reshape(T::DenseTensor, dims::Dims)
  dim(T)==dim(dims) || error("Total new dimension must be the same as the old dimension")
  return tensor(store(T), dims)
end

function Base.reshape(T::DenseTensor, dims::Int...)
  return tensor(store(T), tuple(dims...))
end

# Create an Array that is a view of the Dense Tensor
# Useful for using Base Array functions
array(T::DenseTensor) = reshape(data(store(T)),dims(inds(T)))

function Base.Array{ElT,N}(T::DenseTensor{ElT,N}) where {ElT,N}
  return copy(array(T))
end

function Base.Array(T::DenseTensor{ElT,N}) where {ElT,N}
  return Array{ElT,N}(T)
end

# TODO: call permutedims!(R,T,perm,(r,t)->t)?
function Base.permutedims!(R::DenseTensor{<:Number,N},
                           T::DenseTensor{<:Number,N},
                           perm::NTuple{N,Int}) where {N}
  RA = array(R)
  TA = array(T)
  @strided RA .= permutedims(TA, perm)
  return R
end

function apply!(R::DenseTensor,
                T::DenseTensor,
                f::Function=(r,t)->t)
  RA = array(R)
  TA = array(T)
  @strided RA .= f.(RA, TA)
  return R
end

# Version that may overwrite the result or promote
# and return the result
# TODO: move to tensor.jl?
function permutedims!!(R::Tensor,
                       T::Tensor,
                       perm::NTuple{N,Int},
                       f::Function=(r,t)->t) where {N}
  RA = array(R)
  TA = array(T)
  if !is_trivial_permutation(perm)
    @strided RA .= f.(RA, permutedims(TA, perm))
  else
    @strided RA .= f.(RA, TA)
  end
  return R
end

# TODO: move to tensor.jl?
function Base.permutedims(T::Tensor{<:Number,N},
                          perm::NTuple{N,Int}) where {N}
  Tp = similar(T,permute(inds(T),perm))
  Tp = permutedims!!(Tp, T, perm)
  return Tp
end

Base.permutedims(T::Tensor,
                 perm::Tuple{Vararg{Int}}) = error("Permutation size must match tensor order")

# TODO: move to tensor.jl?
function Base.:*(x::Number,
                 T::Tensor)
  return tensor(x*store(T),inds(T))
end
Base.:*(T::Tensor, x::Number) = x*T

function Base.permutedims!(R::DenseTensor{<:Number,N},
                           T::DenseTensor{<:Number,N},
                           perm,
                           f::Function) where {N}
  RA = array(R)
  TA = array(T)
  @strided RA .= f.(RA, permutedims(TA, perm))
  return R
end

function outer!(R::DenseTensor{ElR},
                T1::DenseTensor{ElT1},
                T2::DenseTensor{ElT2}) where {ElR, ElT1, ElT2}
  if ElT1 != ElT2
    # TODO: use promote instead
    # T1,T2 = promote(T1,T2)

    ElT1T2 = promote_type(ElT1,ElT2)
    if ElT1 != ElT1T2
      # TODO: get this working
      # T1 = ElR.(T1)
      T1 = one(ElT1T2) * T1
    end
    if ElT2 != ElT1T2
      # TODO: get this working
      # T2 = ElR.(T2)
      T2 = one(ElT1T2) * T2
    end
  end

  v1 = data(T1)
  v2 = data(T2)
  RM = reshape(R, length(v1), length(v2))
  #RM .= v1 .* transpose(v2)
  #mul!(RM, v1, transpose(v2))
  BLAS.gemm!('N', 'T', one(ElR), v1, v2, zero(ElR), RM)
  return R
end

function outer!!(R::Tensor,
                 T1::Tensor,
                 T2::Tensor)
  outer!(R,T1,T2)
  return R
end

# TODO: call outer!!, make this generic
function outer(T1::DenseTensor{ElT1},
               T2::DenseTensor{ElT2}) where {ElT1,ElT2}
  array_outer = vec(array(T1))*transpose(vec(array(T2)))
  inds_outer = unioninds(inds(T1),inds(T2))
  return tensor(Dense{promote_type(ElT1,ElT2)}(vec(array_outer)),
                inds_outer)
end
const ⊗ = outer

# TODO: move to tensor.jl?
function contraction_output_type(::Type{TensorT1},
                                 ::Type{TensorT2},
                                 ::Type{IndsR}) where {TensorT1<:Tensor,
                                                       TensorT2<:Tensor,
                                                       IndsR}
  return similar_type(promote_type(TensorT1,TensorT2),IndsR)
end

function contraction_output(::TensorT1,
                            ::TensorT2,
                            indsR::IndsR) where {TensorT1<:DenseTensor,
                                                 TensorT2<:DenseTensor,
                                                 IndsR}
  TensorR = contraction_output_type(TensorT1,TensorT2,IndsR)
  return similar(TensorR,indsR)
end

# TODO: move to tensor.jl?
function contraction_output(T1::Tensor,
                            labelsT1,
                            T2::Tensor,
                            labelsT2,
                            labelsR)
  indsR = contract_inds(inds(T1),labelsT1,inds(T2),labelsT2,labelsR)
  R = contraction_output(T1,T2,indsR)
  return R
end

# TODO: move to tensor.jl?
function contract(T1::Tensor{<:Any, N1},
                  labelsT1,
                  T2::Tensor{<:Any, N2},
                  labelsT2,
                  labelsR = contract_labels(labelsT1,
                                            labelsT2)) where {N1, N2}
  # TODO: put the contract_inds logic into contraction_output,
  # call like R = contraction_ouput(T1,labelsT1,T2,labelsT2)
  #indsR = contract_inds(inds(T1),labelsT1,inds(T2),labelsT2,labelsR)
  R = contraction_output(T1, labelsT1,
                         T2, labelsT2,
                         labelsR)
  # contract!! version here since the output R may not
  # be mutable (like UniformDiag)
  R = contract!!(R, labelsR,
                 T1, labelsT1,
                 T2, labelsT2)
  return R
end

# Move to tensor.jl? Is this generic for all storage types?
function contract!!(R::Tensor{<:Number,NR},
                    labelsR::NTuple{NR},
                    T1::Tensor{<:Number,N1},
                    labelsT1::NTuple{N1},
                    T2::Tensor{<:Number,N2},
                    labelsT2::NTuple{N2},
                    α::Number=1, β::Number=0) where {NR,N1,N2}
  if N1==0
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for scalar ITensor with non-trivial α and β")
    # TODO: replace with an add! function?
    # What about doing `R .= T1[] .* PermutedDimsArray(T2,perm)`?
    perm = getperm(labelsR,labelsT2)
    R = permutedims!!(R,T2,perm,(r,t2)->T1[]*t2)
  elseif N2==0
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for scalar ITensor with non-trivial α and β")
    perm = getperm(labelsR,labelsT1)
    R = permutedims!!(R,T1,perm,(r,t1)->T2[]*t1)
  elseif N1+N2==NR
    (α ≠ 1 || β ≠ 0) &&
      error("contract!! not yet implemented for outer product tensor contraction with non-trivial α and β")
    # TODO: permute T1 and T2 appropriately first (can be more efficient
    # then permuting the result of T1⊗T2)
    # TODO: implement the in-place version directly
    R = outer!!(R, T1, T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR, labelsRp)
    if !is_trivial_permutation(perm)
      Rp = reshape(R, (inds(T1)..., inds(T2)...))
      R = permutedims!!(R, copy(Rp), perm)
    end
  else
    if α ≠ 1 || β ≠ 0
      R = _contract!!(R, labelsR,
                      T1, labelsT1,
                      T2, labelsT2,
                      α, β)
    else
      R = _contract!!(R, labelsR,
                      T1, labelsT1,
                      T2, labelsT2)
    end
  end
  return R
end

# Move to tensor.jl? Overload this function
# for immutable storage types
function _contract!!(R::Tensor, labelsR,
                     T1::Tensor, labelsT1,
                     T2::Tensor, labelsT2,
                     α::Number=1, β::Number=0)
  if α ≠ 1 || β ≠ 0
    contract!(R, labelsR,
              T1, labelsT1,
              T2, labelsT2,
              α, β)
  else
    contract!(R, labelsR,
              T1, labelsT1,
              T2, labelsT2)
  end
  return R
end

function contract!(R::DenseTensor{<:Number,NR},
                   labelsR,
                   T1::DenseTensor{ElT1,N1},
                   labelsT1,
                   T2::DenseTensor{ElT2,N2},
                   labelsT2,
                   α::Number=1,
                   β::Number=0) where {ElT1,ElT2,N1,N2,NR}
  if N1+N2==NR
    outer!(R,T1,T2)
    labelsRp = (labelsT1..., labelsT2...)
    perm = getperm(labelsR,labelsRp)
    if !is_trivial_permutation(perm)
      permutedims!(R, copy(R), perm)
    end
    return R
  end

  props = ContractionProperties(labelsT1,
                                labelsT2,
                                labelsR)
  compute_contraction_properties!(props,
                                  T1, T2, R)

  if ElT1 != ElT2
    # TODO: use promote instead
    # T1, T2 = promote(T1, T2)

    ElR = promote_type(ElT1, ElT2)
    if ElT1 != ElR
      # TODO: get this working
      # T1 = ElR.(T1)
      T1 = one(ElR) * T1
    end
    if ElT2 != ElR
      # TODO: get this working
      # T2 = ElR.(T2)
      T2 = one(ElR) * T2
    end
  end

  _contract!(R,T1,T2,props,α,β)
  return R
end

function _contract!(CT::DenseTensor{El, NC},
                    AT::DenseTensor{El, NA},
                    BT::DenseTensor{El, NB},
                    props::ContractionProperties,
                    α::Number = one(El),
                    β::Number = zero(El)) where {El, NC, NA, NB}
  # TODO: directly use Tensor instead of Array
  C = array(CT)
  A = array(AT)
  B = array(BT)

  tA = 'N'
  if props.permuteA
    pA = NTuple{NA, Int}(props.PA)
    @strided Ap = permutedims(A, pA)
    AM = reshape(Ap, props.dmid, props.dleft)
    tA = 'T'
  else
    #A doesn't have to be permuted
    if Atrans(props)
      AM = reshape(A, props.dmid, props.dleft)
      tA = 'T'
    else
      AM = reshape(A, props.dleft, props.dmid)
    end
  end

  tB = 'N'
  if props.permuteB
    pB = NTuple{NB, Int}(props.PB)
    @strided Bp = permutedims(B, pB)
    BM = reshape(Bp, props.dmid, props.dright)
  else
    if Btrans(props)
      BM = reshape(B, props.dright, props.dmid)
      tB = 'T'
    else
      BM = reshape(B, props.dmid, props.dright)
    end
  end

  # TODO: this logic may be wrong
  if props.permuteC
    # Need to copy here since we will be permuting
    # into C later
    CM = reshape(copy(C), props.dleft, props.dright)
  else
    if Ctrans(props)
      CM = reshape(C, props.dright, props.dleft)
      (AM, BM) = (BM, AM)
      if tA == tB
        tA = tB = (tA == 'T' ? 'N' : 'T')
      end
    else
      CM = reshape(C, props.dleft, props.dright)
    end
  end

  BLAS.gemm!(tA, tB, El(α), AM, BM, El(β), CM)

  if props.permuteC
    pC = NTuple{NC, Int}(props.PC)
    Cr = reshape(CM, props.newCrange...)
    # TODO: use invperm(pC) here?
    @strided C .= permutedims(Cr, pC)
  end
  return C
end

"""
permute_reshape(T::Tensor,pos...)

Takes a permutation that is split up into tuples. Index positions
within the tuples are combined.

For example:

permute_reshape(T,(3,2),1)

First T is permuted as `permutedims(3,2,1)`, then reshaped such
that the original indices 3 and 2 are combined.
"""
function permute_reshape(T::DenseTensor{ElT,NT,IndsT},
                         pos::Vararg{<:Any,N}) where {ElT,NT,IndsT,N}
  perm = flatten(pos...)

  length(perm)≠NT && error("Index positions must add up to order of Tensor ($N)")
  isperm(perm) || error("Index positions must be a permutation")

  dimsT = dims(T)
  indsT = inds(T)
  if !is_trivial_permutation(perm)
    T = permutedims(T, perm)
  end
  N==NT && return T
  newdims = MVector(ntuple(_->eltype(IndsT)(1),Val(N)))
  for i ∈ 1:N
    if length(pos[i])==1
      # No reshape needed, just use the
      # original index
      newdims[i] = indsT[pos[i][1]]
    else
      newdim_i = 1
      for p ∈ pos[i]
        newdim_i *= dimsT[p]
      end
      newdims[i] = eltype(IndsT)(newdim_i)
    end
  end
  newinds = similar_type(IndsT,Val{N})(Tuple(newdims))
  return reshape(T,newinds)
end

# svd of an order-n tensor according to positions Lpos
# and Rpos
function LinearAlgebra.svd(T::DenseTensor{<:Number,N,IndsT},
                           Lpos::NTuple{NL,Int},
                           Rpos::NTuple{NR,Int};
                           kwargs...) where {N,IndsT,NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  UM,S,VM,spec = svd(M;kwargs...)
  u = ind(UM,2)
  v = ind(VM,2)

  Linds = similar_type(IndsT,Val{NL})(ntuple(i->inds(T)[Lpos[i]],Val(NL)))
  Uinds = push(Linds,u)

  # TODO: do these positions need to be reversed?
  Rinds = similar_type(IndsT,Val{NR})(ntuple(i->inds(T)[Rpos[i]],Val(NR)))
  Vinds = push(Rinds,v)

  U = reshape(UM,Uinds)
  V = reshape(VM,Vinds)

  return U,S,V,spec
end

# qr decomposition of an order-n tensor according to
# positions Lpos and Rpos
function LinearAlgebra.qr(T::DenseTensor{<:Number,N,IndsT},
                          Lpos::NTuple{NL,Int},
                          Rpos::NTuple{NR,Int};kwargs...) where {N,IndsT,NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  QM,RM = qr(M;kwargs...)
  q = ind(QM,2)
  r = ind(RM,1)
  # TODO: simplify this by permuting inds(T) by (Lpos,Rpos)
  # then grab Linds,Rinds
  Linds = similar_type(IndsT,Val{NL})(ntuple(i->inds(T)[Lpos[i]],Val(NL)))
  Qinds = push(Linds,r)
  Q = reshape(QM,Qinds)
  Rinds = similar_type(IndsT,Val{NR})(ntuple(i->inds(T)[Rpos[i]],Val(NR)))
  Rinds = pushfirst(Rinds,r)
  R = reshape(RM,Rinds)
  return Q,R
end

# polar decomposition of an order-n tensor according to positions Lpos
# and Rpos
function polar(T::DenseTensor{<:Number,N,IndsT},
               Lpos::NTuple{NL,Int},
               Rpos::NTuple{NR,Int}) where {N,IndsT,NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  UM,PM = polar(M)

  # TODO: turn these into functions
  Linds = similar_type(IndsT,Val{NL})(ntuple(i->inds(T)[Lpos[i]],Val(NL)))
  Rinds = similar_type(IndsT,Val{NR})(ntuple(i->inds(T)[Rpos[i]],Val(NR)))

  # Use sim to create "similar" indices, in case
  # the indices have identifiers. If not this should
  # act as an identity operator
  simRinds = sim(Rinds)
  Uinds = (Linds..., simRinds...)
  Pinds = (simRinds..., Rinds...)

  U = reshape(UM,Uinds)
  P = reshape(PM,Pinds)
  return U,P
end

function LinearAlgebra.exp(T::DenseTensor{ElT,N},
                           Lpos::NTuple{NL,Int},
                           Rpos::NTuple{NR,Int};
                           ishermitian::Bool=false) where {ElT,N,
                                                           NL,NR}
  M = permute_reshape(T,Lpos,Rpos)
  indsTp = permute(inds(T), (Lpos...,Rpos...))
  if ishermitian
    expM = parent(exp(Hermitian(matrix(M))))
    return tensor(Dense{ElT}(vec(expM)), indsTp)
  else
    expM = exp(M)
    return reshape(expM,indsTp)
  end
end

function HDF5.write(parent::Union{HDF5File, HDF5Group},
                    name::String,
                    D::Store) where {Store <: Dense}
  g = g_create(parent,name)
  attrs(g)["type"] = "Dense{$(eltype(Store))}"
  attrs(g)["version"] = 1
  if eltype(D) != Nothing
    write(g,"data",D.data)
  end
end

function HDF5.read(parent::Union{HDF5File,HDF5Group},
                   name::AbstractString,
                   ::Type{Store}) where {Store <: Dense}
  g = g_open(parent,name)
  ElT = eltype(Store)
  typestr = "Dense{$ElT}"
  if read(attrs(g)["type"]) != typestr
    error("HDF5 group or file does not contain $typestr data")
  end
  if ElT == Nothing
    return Dense{Nothing}()
  end
  data = read(g,"data")
  return Dense{ElT}(data)
end

function Base.show(io::IO,
                   mime::MIME"text/plain",
                   T::DenseTensor)
  summary(io,T)
  print_tensor(io,T)
end


text:
module CMultipleBoxes

include("./../../constraints/geometric/constants.jl")
using .Constants

using JuMP

export cons_ss_mult_boxes

function cons_ss_mult_boxes(m)

    @variable(m, does_box_provide_var_support[i=1:length(boxes)-1, j=1:length(boxes)-1; i!=j], Bin)
    @constraint(m, c_how_many_boxes_provide_var_support1[i=1:(length(boxes)-1), j=1:length(boxes)-1; i!=j],
        (1 - m[:does_box_provide_var_support][i, j]) * area_of_box[i] >=
        mult_boxes_ss_var * area_of_box[i] - m[:pbs_box_ij][i, j]
    )

    @constraint(m, c_how_many_boxes_provide_var_support2[i=1:(length(boxes)-1), j=1:length(boxes)-1; i!=j],
        m[:does_box_provide_var_support][i, j] * area_of_box[i] >=
        m[:pbs_box_ij][i, j] - mult_boxes_ss_var * area_of_box[i]
    )

    @expression(m, sum_of_how_many_boxes_provide_var_support[i=1:length(boxes)-1],
        sum([m[:does_box_provide_var_support][i, j] for j in 1:length(boxes)-1 if i!=j])
    )

    @variable(m, is_box_above_ground[i=1:length(boxes)-1], Bin)
    @constraint(m, c_is_box_above_ground1[i=1:length(boxes)-1],
        (1 - is_box_above_ground[i]) * H >=
        0.002 - m[:flbz_front_left_bot_z][i]
    )

    @constraint(m, c_is_box_above_ground2[i=1:length(boxes)-1],
        is_box_above_ground[i] * H >=
        m[:flbz_front_left_bot_z][i] - 0.002
    )

    @constraint(m, c_ensure_two_or_more_support1[i=1:length(boxes)-1],
        number_of_boxes * (1 - is_box_above_ground[i]) >= number_of_boxes - m[:sum_of_how_many_boxes_provide_var_support][i]
    )

    @expression(m, mult_boxes_objective[i=1:(length(boxes)-1), j=1:length(boxes)-1; i!=j],
        (1 - m[:does_box_provide_var_support][i, j])
    )

    return m
end

end

ext: jl


text:
# Giacomo Zanella (2020) Informed Proposals for Local MCMC in Discrete Spaces
# Journal of the American Statistical Association

struct ZanellaMH{T3, V, CS} <: SpinSampler{T3}
    A::T3           # q × N × M
    A₁::T3          # 1 × N × M
    A₁₂::T3         # 1 × 1 × M
    α::T3           # 1 × 1 × M
    I::V            # M element vector
    batchshift::V   # M element vector
    𝕡::T3           # q × N × M
    q::Int          # number of spin states
    N::Int          # number of spins
    M::Int          # sample size
    catsampler::CS
    function ZanellaMH(X::AbstractArray{<:AbstractFloat,3})
        q, N, M = size(X)
        _eltype = promote_type(eltype(X), Float32)
        A  = similar(X, _eltype, q, N, M)
        A₁ = similar(X, _eltype, 1, N, M)
        A₁₂= similar(X, _eltype, 1, 1, M)
        α  = similar(X, _eltype, 1, 1, M)
        I  = similar(X, Int, M)
        batchshift = oftype(I, range(0; step = q*N, length = M))
        𝕡 = similar(X, _eltype, q, N, M)
        catsampler = InverseCDF(𝕡, q*N, M)
        return new{typeof(A), typeof(I), typeof(catsampler)}(
            A, A₁, A₁₂, α, I, batchshift, 𝕡, q, N, M, catsampler
        )
    end
end

(f::ZanellaMH)(X, Δenergy) = zanella_mh!(X, Δenergy, f)

function zanella_mh!(X::AbstractArray{T,3}, Δenergy, buffer = ZanellaMH(X)) where T
    (; A, A₁, A₁₂, α, I, batchshift, 𝕡, q, N, catsampler) = buffer
    _onehalf = convert(T, 1//2)

    ΔE = Δenergy(X)                 # E(xᵠᵏ)-E(x)
    _zanella_proposal!!(𝕡, A₁₂, ΔE) # proposal π(xᵠᵏ|x) and partition function Z(x)
    unitrange_rand!(I, 1:q*N, 𝕡, catsampler)
                                    # flipped spin and proposed spin state for each sample
    I .+= batchshift                # shift index with sample index

    α .= A₁₂                        # α = Z(x) = ∑ᵩₖ exp(-1//2 * [E(xᵠᵏ)-E(x)])

    fill!(A, zero(T))               # initialize
    fill!(@view(A[I]), one(T))      # onehot encode flipped spin and proposed spin state
    sum!(A₁, A)                     # onehot encode site of flipped spin
    @. A = X + (A - A₁ * X)         # proposed state
    @. A = A > _onehalf             # ensure the entries are in (0,1)

    ΔẼ = Δenergy(A)                 # E(x̃ᵠᵏ)-E(x̃)
    _zanella_proposal!!(𝕡, A₁₂, ΔẼ, onlypartionfunc=true)
                                    # Z(x̃)
    α ./= A₁₂                       # α = Z(x)/Z(x̃)

    @. α = min(α, 1)                # α = min(1, Z(x)/Z(x̃))
    @. α = α ≥ rand()               # accept proposal

    @. X = (1-α) * X + α * A        # update X
    @. X = X > _onehalf             # ensure the entries are in (0,1)
    return X
end

function _zanella_proposal!!(
        𝕡::AbstractArray{T,3}, A₁₂::AbstractArray{T,3}, ΔE::AbstractArray{T,3};
        onlypartionfunc = false
    ) where T<:Real
    _typemax = typemax(T)
    _onehalf = convert(T, 1//2)
    @. 𝕡 = -_onehalf * ΔE               # log probability of q(zᵠᵏ|z)
    @. 𝕡-= (X > _onehalf) * _typemax    # disallow transition to current state
    sum!(exp, A₁₂, 𝕡)                   # partition function
    onlypartionfunc || softmax!(𝕡, dims = 1:2)  # normalized proposal q(zᵠᵏ|z)
end


text:
# *********************************************************************************
# REopt, Copyright (c) 2019-2020, Alliance for Sustainable Energy, LLC.
# All rights reserved.
#
# Redistribution and use in source and binary forms, with or without modification,
# are permitted provided that the following conditions are met:
#
# Redistributions of source code must retain the above copyright notice, this list
# of conditions and the following disclaimer.
#
# Redistributions in binary form must reproduce the above copyright notice, this
# list of conditions and the following disclaimer in the documentation and/or other
# materials provided with the distribution.
#
# Neither the name of the copyright holder nor the names of its contributors may be
# used to endorse or promote products derived from this software without specific
# prior written permission.
#
# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
# ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
# WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED.
# IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT,
# INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING,
# BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
# LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE
# OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED
# OF THE POSSIBILITY OF SUCH DAMAGE.
# *********************************************************************************
module REopt

export
    Scenario,
    BAUScenario,
    REoptInputs,
    run_reopt,
    build_reopt!,
    reopt_results,
    simulate_outages,
    add_variables!,
    add_objective!,
    LinDistFlow,
    MPCScenario,
    MPCInputs,
    run_mpc,
    build_mpc!

import HTTP
import JSON
using LinDistFlow  # required to export LinDistFlow
import LinDistFlow
const LDF = LinDistFlow
using JuMP
using JuMP.Containers: DenseAxisArray
using Logging
using DelimitedFiles
using Dates
import MathOptInterface
import Dates: daysinmonth, Date, isleapyear
import DelimitedFiles: readdlm
const MOI = MathOptInterface
using Shapefile
using ArchGDAL
using Roots: fzero  # for IRR
global hdl = nothing

const EXISTING_BOILER_EFFICIENCY = 0.8

const GAL_PER_M3 = 264.172  # [gal/m^3]
const KWH_PER_GAL_DIESEL = 40.7  # [kWh/gal_diesel]
const KWH_PER_MMBTU = 293.07107  # [kWh/mmbtu]
const KWH_THERMAL_PER_TONHOUR = 3.51685
const TONNE_PER_LB = 1/2204.62  # [tonne/lb]

include("keys.jl")
include("core/types.jl")
include("core/utils.jl")

include("core/settings.jl")
include("core/site.jl")
include("core/financial.jl")
include("core/pv.jl")
include("core/wind.jl")
include("core/energy_storage/storage.jl")
include("core/energy_storage/electric_storage.jl")
include("core/energy_storage/thermal_storage.jl")
include("core/generator.jl")
include("core/doe_commercial_reference_building_loads.jl")
include("core/electric_load.jl")
include("core/existing_boiler.jl")
include("core/existing_chiller.jl")
include("core/absorption_chiller.jl")
include("core/flexible_hvac.jl")
include("core/heating_cooling_loads.jl")
include("core/electric_utility.jl")
include("core/prodfactor.jl")
include("core/urdb.jl")
include("core/electric_tariff.jl")
include("core/chp.jl")
include("core/scenario.jl")
include("core/bau_scenario.jl")
include("core/reopt_inputs.jl")
include("core/bau_inputs.jl")
include("core/cost_curve.jl")

include("constraints/outage_constraints.jl")
include("constraints/storage_constraints.jl")
include("constraints/flexible_hvac.jl")
include("constraints/load_balance.jl")
include("constraints/tech_constraints.jl")
include("constraints/electric_utility_constraints.jl")
include("constraints/generator_constraints.jl")
include("constraints/cost_curve_constraints.jl")
include("constraints/production_incentive_constraints.jl")
include("constraints/thermal_tech_constraints.jl")
include("constraints/chp_constraints.jl")
include("constraints/operating_reserve_constraints.jl")
include("constraints/battery_degradation.jl")

include("mpc/structs.jl")
include("mpc/scenario.jl")
include("mpc/inputs.jl")
include("mpc/constraints.jl")

include("core/techs.jl")

include("results/results.jl")
include("results/electric_tariff.jl")
include("results/electric_utility.jl")
include("results/proforma.jl")
include("results/financial.jl")
include("results/generator.jl")
include("results/pv.jl")
include("results/electric_storage.jl")
include("results/thermal_storage.jl")
include("results/outages.jl")
include("results/wind.jl")
include("results/electric_load.jl")
include("results/existing_boiler.jl")
include("results/existing_chiller.jl")
include("results/absorption_chiller.jl")
include("results/chp.jl")
include("results/flexible_hvac.jl")

include("core/reopt.jl")
include("core/reopt_multinode.jl")

include("outagesim/outage_simulator.jl")

include("lindistflow/extend.jl")

include("mpc/results.jl")
include("mpc/model.jl")

end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
module CurrenciesBase

import ..Currencies
using ..CurrencyData

import Base: +, -, *, /, ==
using Compat

# Exports
export AbstractMonetary, Monetary
export currency, decimals, majorunit, @usingcurrencies
export currencyinfo, iso4217num, iso4217alpha, shortsymbol, longsymbol
export newcurrency!, @usingcustomcurrency

# Monetary type, currencies, and arithmetic
include("monetary.jl")
include("currency.jl")
include("arithmetic.jl")
include("mixed.jl")

# Custom currencies and macros
include("usingcurrencies.jl")
include("custom.jl")

end  # module CurrenciesBase


text:
function MPT_SVD(arrayaux::Array{Float64},middlesizes::Matrix{Int64},acc::Float64)
    tensor = Array{Float64};
    array  = copy(arrayaux);
    if size(middlesizes,1) == 1
        tensor = reshape(array,Tuple(middlesizes));
        return TT_SVD(tensor,acc);
    else
        sizes   = Tuple(reshape(middlesizes',(length(middlesizes),1)));
        array   = reshape(array,sizes);
        permind = Tuple([ (i-1)*size(middlesizes,2)+j for j in 1:size(middlesizes,2) for i in 1:size(middlesizes,1) ]);
        array   = permutedims(array,permind);
        resind  = Tuple([prod(col) for col in eachcol(middlesizes)]);
        tensor  = reshape(array,resind)
        mps,err = TT_SVD(tensor,acc);
        rnks    = rank(mps);
        return MPT( [reshape(mps[i],(rnks[i][1], middlesizes[:,i]..., rnks[i][2])) for i = 1:order(mps)] ), err
    end
end

# TT-ALS for diagonal Matrix
function MPT_SVD(diagonal::Diagonal,middlesizes::Matrix,acc::Float64)
    # initial ttm with diagonal core
    N      = length(middlesizes);
    tensor = reshape(diag(diagonal),Tuple(middlesizes));
    tt,err = TT_SVD(tensor,acc);
    cores  = Vector{Array{Float64,4}}(undef,N);
    rnks   = rank(tt);
    for i = 1:N
        cores[i] = zeros(rnks[i][1],middlesizes[i],middlesizes[i],rnks[i][2])
        for ri = 1:rnks[i][1]
            for rii = 1:rnks[i][2]
                cores[i][ri,:,:,rii] = Matrix(Diagonal(tt[i][ri,:,rii]))
            end
        end
    end
    return MPT(cores)
end

# TT-ALS for symmetric Matrix
function MPT_SVD(sym::Symmetric,middlesizes::Matrix,acc::Float64)
    # initial ttm with diagonal core
    d          = Int(length(middlesizes)/2)
    symM       = Matrix(sym)
    tensor     = reshape(symM,Tuple(middlesizes));
    Uttm,S,err = TT_SVD(tensor,d,acc);
    r          = rank(Uttm)

    cores   = Vector{Array{Float64,3}}(undef,2*d)
    coresL  = Vector{Array{Float64,4}}(undef,d)
    coresU  = Vector{Array{Float64,4}}(undef,d)

    di     = collect(d:-1:1)
    for i = 1:d
        cores[i]   = Uttm[i]
        cores[d+i] = permutedims(Uttm[di[i]],[3,2,1])

        if i<d
            coresL[i]  = reshape(Uttm[i],(size(Uttm[i],1),size(Uttm[i],2),1,size(Uttm[i],3)))
            coresU[i]  = coresL[i]
        end
    end

    coresU[d]  = reshape(cores[d],(size(cores[d],1),size(cores[d],2),size(cores[d],3),1))
    cores[d]   = nmodeproduct(Matrix(sqrt.(Diagonal(S))),cores[d],3)
    cores[d+1] = nmodeproduct(Matrix(sqrt.(Diagonal(S))),cores[d+1],1)
    coresL[d]  = reshape(cores[d],(size(cores[d],1),size(cores[d],2),size(cores[d],3),1))

    return MPT(cores),MPT(coresL),MPT(coresU),S,err
end

max_stars_repo_name: clarazen/BigMat

max_issues_repo_name: clarazen/BigMat

max_forks_repo_name: clarazen/BigMat


text:
using Oceananigans

include("utils.jl")

Nx, Ny, Nz = 256, 1, 256
Lx, Ly, Lz = 2000, 1, 1000
Nt, Δt = 5000, 10

model = Model(N=(Nx, Ny, Nz), L=(Lx, Ly, Lz), ν=4e-2, κ=4e-2)

# Set a temperature perturbation with a Gaussian profile located at the center
# of the vertical slice. It roughly corresponds to a background temperature of
# T = 282.99 K and a bubble temperature of T = 283.01 K.
xC, zC = model.grid.xC, model.grid.zC
hot_bubble_perturbation = 0.01 * exp.(-100 * ((xC .- Lx/2).^2 .+ (zC .+ Lz/2)'.^2) / (Lx^2 + Lz^2))
model.tracers.T.data .= 282.99 .+ 2 .* reshape(hot_bubble_perturbation, (Nx, Ny, Nz))

# Add a NetCDF output writer that saves NetCDF files to the current directory
# "." with a filename prefix of "thermal_bubble_" every 10 iterations.
nc_writer = NetCDFOutputWriter(dir=".", prefix="thermal_bubble_", frequency=10)
push!(model.output_writers, nc_writer)

time_step!(model, Nt, Δt)

make_vertical_slice_movie(model, nc_writer, "T", Nt, Δt, model.eos.T₀)


text:
module Plotting

export countplot
export countplot!
export confidenceinterval
export confidenceinterval!

include("countplot.jl")
include("confidenceinterval.jl")

end


text:
module Hilbert
include("common.jl")
abstract type Basis end
abstract type Fock <: Basis end

"""
Fock Basis for fermion and hard-core boson, each site has two states with occupation number 0 and 1
e.g. |101>↓|100>↑ means a 3-site fock state, site 1: one up one down, site 2: no particle, site 3: one down
"""
struct BinaryFock <: Fock
    idx::Vector{Int} # index of basises
    sites::Int # site number
    # N::Int # particle number, only meaningful for conserved states
    # Sz::Int # magnetization, only meaningful for conserved states
    """
    # Arguments
    - N: Specify the total particle number by a number, an array of number or the symbol :all
    - Sz: Specify the total z-magnetization number by a number, an array of number or the symbol :all
    """
    function BinaryFock(sites, _N=:all, _Sz=:all)
        idx = Vector{Int}(undef, 0)
        for i in 0:2^(2sites) - 1
            Nup, Ndown = occupationTotal(i, sites)
            if _N == :all || Nup + Ndown in collect(_N)
                if _Sz == :all || Nup - Ndown in collect(_Sz)
                    append!(idx, i)
                end
            end
        end
        return new(idx, sites)
end
end

function state2idx(nup::Vector{Int}, ndown::Vector{Int})
    @assert length(nup) == length(ndown)
    sites = length(nup)
    idxup, idxdown = 0, 0
    for i in 1:sites
        idxup += 2^(sites - i) * nup[i]
        idxdown += 2^(sites - i) * ndown[i]
    end
    idx = idxdown * 2^sites + idxup
    return idx
end

function idx2state(idx::Int, sites::Int)
    idxdown = idx >> sites
    idxup = idx % 2^sites
    nup, ndown = zeros(Int, sites), zeros(Int, sites)
    for i in 1:sites
        nup[sites - i + 1] = idxup % 2
        ndown[sites - i + 1] = idxdown % 2
        idxup = idxup >> 1
        idxdown = idxdown >> 1
    end
    return nup, ndown
end

function occupationTotal(idx::Int, sites::Int)
    idxdown = idx >> sites
    idxup = idx % 2^sites
    nup, ndown = 0, 0
    for i in 1:sites
        nup += idxup % 2
        ndown += idxdown % 2
        idxup = idxup >> 1
        idxdown = idxdown >> 1
    end
    return nup, ndown
end

function occupation(idx, sites, site, spin)
    rsite = sites - site + 1
    if spin == UP
        idxup = idx % 2^sites
        idxup = idxup >> (rsite - 1)
        return idxup % 2
    else
        idxdown = idx >> sites
        idxdown = idxdown >> (rsite - 1)
        return idxdown % 2
end
end

function creation(basis::BinaryFock, site, spin)
    Dim = 4^basis.sites
    # Dim = length(basis.idx)
    c⁺ = spzeros(Float, Dim, Dim)
    for (ki, ket) in enumerate(basis.idx)
        nup, ndown = idx2state(ket, basis.sites)
        if spin == UP
            (nup[site] == 1) && continue
            sign = (-1.0)^(sum(ndown) + sum(nup[1:site])) # should be done before nup[site] is set to one
            nup[site] = 1
        else
            (ndown[site] == 1) && continue
            sign = (-1.0)^(sum(ndown[1:site])) # should be done before nup[site] is set to one
            ndown[site] = 1
        end
        idxbra = state2idx(nup, ndown)
        bi = indexin(idxbra, basis.idx)[1]
        @assert isnothing(bi) == false
        # if isnothing(bi)
        #     continue
        # end
        sbu, sbd = idx2state(idxbra, 2)
        sku, skd = idx2state(ket, 2)
        # if spin == UP
        #     println("<$(sbd[1])$(sbd[2])$(sbu[1])$(sbu[2])|cu⁺$site|$(skd[1])$(skd[2])$(sku[1])$(sku[2])> = ", sign)
        # else
        #     println("<$(sbd[1])$(sbd[2])$(sbu[1])$(sbu[2])|cd⁺$site|$(skd[1])$(skd[2])$(sku[1])$(sku[2])> = ", sign)
        # end
        c⁺[bi, ki] = sign
    end
    return c⁺
end

# state=BinaryFock(:fermi, 1)
end


text:
# * Level

"""
    Level(csf, J)

Given a [`CSF`](@ref) with a final [`Term`](@ref), a `Level`
additionally specifies a total angular momentum ``J``. By further
specifying a projection quantum number ``M_J``, we get a
[`State`](@ref).
"""
struct Level{O,IT,T}
    csf::CSF{O,IT,T}
    J::HalfInt
    function Level(csf::CSF{O,IT,T}, J::HalfInteger) where {O,IT,T}
        J ∈ J_range(last(csf.terms)) ||
            throw(ArgumentError("Invalid J = $(J) for term $(last(csf.terms))"))
        new{O,IT,T}(csf, J)
    end
end

Level(csf::CSF, J::Integer) = Level(csf, convert(HalfInteger, J))

function Base.show(io::IO, level::Level)
    write(io, "|")
    show(io, level.csf)
    write(io, ", J = $(level.J)⟩")
end

"""
    weight(l::Level)

Returns the statistical weight of the [`Level`](@ref) `l`, i.e. the
number of possible microstates: ``2J+1``.

# Example

```jldoctest
julia> l = Level(first(csfs(c"1s 2p")), 1)
|1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1⟩

julia> weight(l)
3
```
"""
weight(l::Level) = convert(Int, 2l.J + 1)

Base.:(==)(l1::Level, l2::Level) = ((l1.csf == l2.csf) && (l1.J == l2.J))

# It makes no sense to sort levels with different electron configurations
Base.isless(l1::Level, l2::Level) = ((l1.csf < l2.csf)) ||
    ((l1.csf == l2.csf)) && (l1.J < l2.J)

"""
    J_range(term::Term)

List the permissible values of the total angular momentum ``J`` for
`term`.

# Examples

```jldoctest
julia> J_range(T"¹S")
0:0

julia> J_range(T"²S")
1/2:1/2

julia> J_range(T"³P")
0:2

julia> J_range(T"²D")
3/2:5/2
```
"""
J_range(term::Term) = abs(term.L-term.S):(term.L+term.S)

"""
    J_range(J)

The permissible range of the total angular momentum ``J`` in ``jj``
coupling is simply the value of ``J`` for the final term.
"""
J_range(J::HalfInteger) = J:J
J_range(J::Integer) = J_range(convert(HalfInteger, J))

"""
    levels(csf)

Generate all permissible [`Level`](@ref)s given `csf`.

# Examples

```jldoctest
julia> levels.(csfs(c"1s 2p"))
2-element Array{Array{Level{Orbital{Int64},Term,Seniority},1},1}:
 [|1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1⟩]
 [|1s(₁²S|²S) 2p(₁²Pᵒ|³Pᵒ)-, J = 0⟩, |1s(₁²S|²S) 2p(₁²Pᵒ|³Pᵒ)-, J = 1⟩, |1s(₁²S|²S) 2p(₁²Pᵒ|³Pᵒ)-, J = 2⟩]

julia> levels.(csfs(rc"1s 2p"))
2-element Array{Array{Level{RelativisticOrbital{Int64},HalfIntegers.Half{Int64},Seniority},1},1}:
 [|1s(₁1/2|1/2) 2p(₁3/2|1)-, J = 1⟩]
 [|1s(₁1/2|1/2) 2p(₁3/2|2)-, J = 2⟩]

julia> levels.(csfs(rc"1s 2p-"))
2-element Array{Array{Level{RelativisticOrbital{Int64},HalfIntegers.Half{Int64},Seniority},1},1}:
 [|1s(₁1/2|1/2) 2p-(₁1/2|0)-, J = 0⟩]
 [|1s(₁1/2|1/2) 2p-(₁1/2|1)-, J = 1⟩]
```
"""
levels(csf::CSF) = sort([Level(csf,J) for J in J_range(last(csf.terms))])

# * State

@doc raw"""
    State(level, M_J)

A states defines, in addition to the total angular momentum ``J`` of
`level`, its projection ``M_J\in -J..J``.
"""
struct State{O,IT,T}
    level::Level{O,IT,T}
    M_J::HalfInt
    function State(level::Level{O,IT,T}, M_J::HalfInteger) where {O,IT,T}
        abs(M_J) ≤ level.J ||
            throw(ArgumentError("Invalid M_J = $(M_J) for level with J = $(level.J)"))
        new{O,IT,T}(level, M_J)
    end
end
State(level::Level, M_J::Integer) = State(level, convert(HalfInteger, M_J))

function Base.show(io::IO, state::State)
    write(io, "|")
    show(io, state.level.csf)
    write(io, ", J = $(state.level.J), M_J = $(state.M_J)⟩")
end

Base.:(==)(a::State,b::State) = a.level == b.level && a.M_J == b.M_J
Base.isless(a::State,b::State) = a.level < b.level || a.level == b.level && a.M_J < b.M_J

"""
   states(level)

Generate all permissible [`State`](@ref) given `level`.

# Example

```jldoctest
julia> l = Level(first(csfs(c"1s 2p")), 1)
|1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1⟩

julia> states(l)
3-element Array{State{Orbital{Int64},Term,Seniority},1}:
 |1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = -1⟩
 |1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = 0⟩
 |1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = 1⟩
```
"""
function states(level::Level{O,IT,T}) where {O,IT,T}
    J = level.J
    [State(level, M_J) for M_J ∈ -J:J]
end

"""
    states(csf)

Directly generate all permissible [`State`](@ref)s for `csf`.

# Example

```jldoctest
julia> c = first(csfs(c"1s 2p"))
1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-

julia> states(c)
1-element Array{Array{State{Orbital{Int64},Term,Seniority},1},1}:
 [|1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = -1⟩, |1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = 0⟩, |1s(₁²S|²S) 2p(₁²Pᵒ|¹Pᵒ)-, J = 1, M_J = 1⟩]
```
"""
states(csf::CSF) = states.(levels(csf))

export Level, weight, J_range, levels, State, states


text:
# example 7: argument groups

using ArgParse

function main(args)

    s = ArgParseSettings("Example 7 for argparse.jl: " *
                         "argument groups.")

    add_arg_group(s, "stack options") # add a group and sets it as the default
    @add_arg_table s begin            # all options (and arguments) in this table
                                      # will be assigned to the newly added group
        "--opt1"
            action = :append_const
            arg_type = ByteString
            constant = "O1"
            dest_name = "O_stack"
            help = "append O1 to the stack"
        "--opt2"
            action = :append_const
            arg_type = ByteString
            constant = "O2"
            dest_name = "O_stack"
            help = "append O2 to the stack"
    end

    add_arg_group(s, "weird options", "weird") # another group, this time with a tag which allows
                                               # to refer to it

    set_default_arg_group(s, "weird") # set the default group (useless here, since we just added it)

    @add_arg_table s begin
        "--awkward-option"
            nargs = '+'
            action = :append_arg
            dest_name = "awk"
            arg_type = ByteString
            range_tester = (x->x=="X"||x=="Y")
            metavar = "XY"
            help = "either X or Y; all XY's are " *
                   "stored in chunks"
    end

    set_default_arg_group(s) # reset the default arg group (which means arguments
                             # are automatically assigned to commands/options/pos.args
                             # groups)
    @add_arg_table s begin
        "-k"
            action = :store_const
            default = 0
            constant = 42
            help = "provide the answer"
        "--şİłłÿ"
            action = :store_true
            help = "an option with a silly name"
            group = "weird"   # this overrides the default group: this option
                              # will be grouped together with --awkward-option
    end

    parsed_args = parse_args(args, s)
    println("Parsed args:")
    for (key,val) in parsed_args
        println("  $key  =>  $(repr(val))")
    end
end

main(ARGS)


text:
using Select
using Test

function select_block_test(t1, t2, t3, t4)
    c1 = Channel{Symbol}(1)
    c2 = Channel{Int}(1)
    c3 = Channel(1)

    put!(c3,1)

    @async begin
        sleep(t1)
        put!(c1,:a)
    end

    @async begin
        sleep(t2)
        put!(c2,1)
    end

    @async begin
        sleep(t3)
        take!(c3)
    end

    task = @async begin
        sleep(t4)
        :task_done
    end

    @select begin
        c1 |> x           => "Got $x from c1"
        c2                =>  "Got a message from c2"
        c3 <| :write_test => "Wrote to c3"
        task |> z         => begin
            "Task finished with $z"
        end
    end
end

@testset "@select blocking" begin
    @test select_block_test(.5, 1, 1, 1) == "Got a from c1"
    @test select_block_test(1, .5, 1, 1) == "Got a message from c2"
    @test select_block_test(1, 1, .5, 1) == "Wrote to c3"
    @test select_block_test(1, 1, 1, .5) == "Task finished with task_done"
end

@testset "@select blocking, already ready" begin
    ch = Channel(1)
    put!(ch, 1)
    @test (@select begin
        ch |> x             => x
        @async(sleep(0.3))  => "timeout"
    end) == 1

    # Can immediately put into a buffered channel
    ch = Channel(1)
    @test (@select begin
        ch <| 1             => take!(ch)
        @async(sleep(0.3))  => "timeout"
    end) == 1
end

function select_nonblock_test(test)
    c = Channel(1)
    c2 = Channel(1)
    put!(c2, 1)
    if test == :take
        put!(c, 1)
    elseif test == :put
        take!(c2)
    elseif test == :default
    end

    @select begin
        c |> x => "Got $x from c"
        c2 <| 1 => "Wrote to c2"
        _ => "Default case"
    end
end

@testset "@select non-blocking" begin
    @test select_nonblock_test(:take) == "Got 1 from c"
    @test select_nonblock_test(:put) == "Wrote to c2"
    @test select_nonblock_test(:default) == "Default case"
end

@testset "self-references" begin
    # Test multiple times because the deadlock was only triggered if Task 1 is run before
    # Task 2 or Task 3. Running it multiple times will hopefully cover all the cases.
    for _ in 1:10
        # Test that the two clauses in `@select` can't trigger eachother (Here, the problem
        # would be if Task 1 puts into ch, then Task 2 sees Task 1 waiting, and thinks ch is
        # ready to take! and attempts to incorrectly proceed w/ the take!.) Instead, this should
        # timeout, since no one else is putting or taking on `ch`.
        ch = Channel()
        @test @select(begin
            ch <| "hi"          => "put"
            # This take!(ch) should not be triggered by the put
            ch |> x             => "take! |> $x"
            # This wait(ch) should also not be triggered by the put
            ch                  => "waiting on ch"
            @async(sleep(0.3))  => "timeout"
        end) == "timeout"
    end
end

# Other waitable things
# Conditions
    # Note that we may not support Base.Condition() because it cannot be used in multithreaded code.
    ## Simple condition test: wait until a condition is notified
    #c = Condition()
    #@async while notify(c) != 1 sleep(0.5) end
    #@test @select(begin
    #    c => "c"
    #end) == "c"

@sync begin
    c = Threads.Condition()
    t = @async begin
        success = false
        while true
            lock(c)
            success = notify(c) > 0
            unlock(c)
            if success break; end
            sleep(0.5)
        end
    end
    lock(c)
    @test @select(begin
        c => (unlock(c); "c")
    end) == "c"
end

# Multiple waiters
@sync begin
    coordinator = Channel()

    c1 = Threads.Condition()
    c2 = Threads.Condition()

    wait_on_conds() = begin
        lock(c1)
        lock(c2)
        put!(coordinator, 0)
        @select begin
            c1  => begin
                put!(coordinator, 0)
                unlock(c2); unlock(c1);
            end
            c2  => begin
                put!(coordinator, 0)
                unlock(c2); unlock(c1);
            end
        end
    end
    t1 = @async wait_on_conds()
    t2 = @async wait_on_conds()

    take!(coordinator); take!(coordinator)  # Wait for both tasks to start waiting

    # Now, everyone is waiting, so notifying c1 will wake up _both_ select tasks (and kill both c2 siblings)
    lock(c1)
    @test notify(c1) == 2
    unlock(c1)
    yield()
    yield()

    take!(coordinator); take!(coordinator)  # Wait for both tasks to finish running

    # Both select tasks will have killed their sibling clauses, so no one is listening on c2
    lock(c2)
    @test notify(c2) == 0
    unlock(c2)
end


text:
struct ConditionalTrace{T, S}
    FT::T
    plane::S
end
Base.IteratorSize(::ConditionalTrace) = Iterators.SizeUnknown()

"""
    conditional_trace(trace::Trace, (p, n))

ConditionalTrace `trace` on hyperplane `H` through p with norml n. Returns iterable object
iterating pairs `t => x` such that `x ∈ H`.

Iteration changes the vector `x` inplace,
`collect` creates necessary copies.
"""
conditional_trace(FT, f) = ConditionalTrace(FT, f)
function hittingtime(x, θ, (p, n), FT::Union{ZigZag, BouncyParticle})
    dot(p - x, n)/dot(θ, n)
end

function Base.iterate(D::ConditionalTrace{<:FactTrace})
    FT = D.FT
    t, x, θ = FT.t0, copy(FT.x0), copy(FT.θ0)
    iterate(D, (t, x, θ, 1))
end

function Base.iterate(D::ConditionalTrace{<:PDMPTrace})
    FT = D.FT
    t, x, θ, f = FT.t0, copy(FT.x0), copy(FT.θ0), copy(FT.f)
    iterate(D, (t, x, θ, f, 1))
end

function Base.iterate(D::ConditionalTrace{<:FactTrace{T}}, (t, x, θ, k)) where {T<:Union{FactBoomerang,ZigZag,JointFlow}}
    (p, n) = D.plane
    FT = D.FT
    while true
        k > length(FT.events) && return nothing
        τ = hittingtime(x, θ, (p, n), FT.F)
        ti, i, xi, θi = FT.events[k]
        if τ > 0 && t + τ <= ti
            t, x, θ = move_forward!(τ, t, x, θ, FT.F)
            return t => x, (t, x, θ, k)
        else # move to ti to change direction
            t, x, θ = move_forward!(ti - t, t, x, θ, FT.F)
            t = ti
            x[i] = xi
            θ[i] = θi
            k = k + 1
        end
    end
end

function Base.iterate(D::ConditionalTrace{<:PDMPTrace{T}}, (t, x, θ, f, k)) where {T<:Union{Boomerang,BouncyParticle}}
    (p,n) = D.plane
    FT = D.FT
    while true
        k > length(FT.events) && return nothing
        tn = first(FT.events[k])
        τ = hittingtime(x, θ, (p, n), FT.F)
        if τ > 0 && t + τ < tn
            t, x, θ = smove_forward!(τ, t, x, θ, f, FT.F)
            return t => x, (t, x, θ, f, k)
        else # move not more than to ti to change direction
            Δt = tn - t
            t, x, θ = smove_forward!(Δt, t, x, θ, f, FT.F)
            t, x, θ, f_ = FT.events[k]
            if !isnothing(f_)
                f = f_
            end
            k = k + 1
        end
    end
end

function Base.collect(D::ConditionalTrace{<:FactTrace})
    collect(t=>copy(x) for (t, x) in D)
end
function Base.collect(D::ConditionalTrace{<:PDMPTrace})
    collect(t=>copy(x) for (t, x) in D)
end

#####

conditional_trace(FS::NotFactSampler, f) = ConditionalTrace(FS, f)

function Base.iterate(D::ConditionalTrace{<:NotFactSampler})
    FS = D.FT
    ϕ = iterate(FS)
    ϕ === nothing && return nothing
    t, (x, θ) = deepcopy(FS.u0)
    ev, state = ϕ
    _, _, _, f_ = ev # CHECKME?
    if isnothing(f_)
        f = ones(Bool, length(x)) # simplified all moving
    else
        f = f_
    end

    iterate(D, ((t, x, θ, f), ev, state))
end

function Base.iterate(D::ConditionalTrace{<:NotFactSampler}, ((t, x, θ, f), ev, state))
    (p,n) = D.plane
    FS = D.FT
    while true
        tn = first(ev)
        τ = hittingtime(x, θ, (p, n), FS.F)
#        println(t, " ", τ, " ", tn, " ", x[end])
        if τ > 0 && t + τ < tn
            t, x, θ = smove_forward!(τ, t, x, θ, f, FS.F)
            return (t => copy(x)), ((t, x, θ, f), ev, state)
        else
            Δt = tn - t
            t, x, θ = smove_forward!(Δt, t, x, θ, f, FS.F)
        end
        t, x_, θ_, f_ = ev
        x .= x_
        θ .= θ_
        ϕ = iterate(FS, state)
        ϕ === nothing && return nothing
        ev, state = ϕ
        if !isnothing(f_)
            f = f_
        end
    end
end


text:
using Base: splat

test = """
    v...>>.vv>
    .vv>>.vv..
    >>.>v>...v
    >>v>>.>.v.
    v>v.vv.v..
    >.>>..v...
    .vv..>.>v.
    v.v..>>v.v
    ....v..v.>
    """

function parse_input(input::AbstractString)
    lines = split.(split(input, '\n'; keepempty = false), "") |>
        splat(hcat) |> permutedims
    return map(lines) do c
        return c == ">" ? 1 : c == "v" ? 2 : 0
    end
end

wrap_around(i, s) = (i - 1) % s + 1
right(C, i) = CartesianIndex(i[1], wrap_around(i[2] + 1, size(C, 2)))
down(C, i) = CartesianIndex(wrap_around(i[1] + 1, size(C, 1)), i[2])

function move!(C, dir)::Bool
    m = [right, down][dir]
    ri = findall(keys(C)) do i
        C[i] == dir && C[m(C, i)] == 0
    end
    isempty(ri) && return false
    C[ri] .= 0
    C[m.(Ref(C), ri)] .= dir
    return true
end

function step!(C)::Bool
    rm = move!(C, 1)
    dm = move!(C, 2)
    return rm || dm
end

function run(C)
    C = copy(C)
    s = 1
    while step!(C)
        s += 1
    end
    return s
end

q1(C) = run(C)

let C = parse_input(test)
    @assert q1(C) == 58
end

C = parse_input(read("day25.in", String))
println("Q1: ", q1(C))

max_stars_repo_name: dhanak/competitive-coding

max_issues_repo_name: dhanak/competitive-coding

max_forks_repo_name: dhanak/competitive-coding


text:
# For the system of ODEs du/dx = df(x,u), take one step with eighth order
# (ninth order locally) Runge Kutta.  y contains initial values at the
# start point x=t.  The final values at x=t+h are given in ys.  work
# is scratch space of size 7n, where n is the dimension of y and ys.
# Method of G. J. Cooper and J. H. Verner, SIAM Journal of Numerical
# Analysis 9, pp. 389-405 (1972).
function rk8!(df, t, h, y, ys, work)
    n = length(y)
    rk8aux(df, t, h, y, ys, view(work,1:n),
           view(work,  n+1:2*n), view(work,2*n+1:3*n), view(work,3*n+1:4*n),
           view(work,4*n+1:5*n), view(work,5*n+1:6*n), view(work,6*n+1:7*n))
end
function rk8aux(df, t, h, y, ys, yp, f1,f2,f3,f4,f5,f6)
    n = length(y)
    df(t, y, yp)
    for k = 1:n
        ys[k] = y[k] + h * (1/2) * yp[k]
    end
    df(t+h*(1/2), ys, f1)
    for k = 1:n
        ys[k] = y[k] + h * (1/4) * (yp[k] + f1[k])
    end
    df(t+h*(1/2), ys, f2)
    for k = 1:n
        ys[k] = y[k] + h * ((1/7) * yp[k] +
                            (-7+3*sqrt(21))/98 * f1[k] +
                            (21-5*sqrt(21))/49 * f2[k] )
    end
    df(t+h*(7-sqrt(21))/14, ys, f1)
    for k = 1:n
        ys[k] = y[k] + h * ((11-sqrt(21))/84 * yp[k] +
                            (18-4*sqrt(21))/63 * f2[k] +
                            (21+sqrt(21))/252 * f1[k])
    end
    df(t+h*(7-sqrt(21))/14, ys, f3)
    for k = 1:n
        ys[k] = y[k] + h * ((5-sqrt(21))/48 * yp[k] +
                            (9-sqrt(21))/36 * f2[k] +
                            (-231-14*sqrt(21))/360 * f1[k] +
                            (63+7*sqrt(21))/80 * f3[k])
    end
    df(t+h*(1/2), ys, f4)
    for k = 1:n
        ys[k] = y[k] + h * ((10+sqrt(21))/42 * yp[k] +
                            (-432-92*sqrt(21))/315 * f2[k] +
                            (633+145*sqrt(21))/90 * f1[k] +
                            (-504-115*sqrt(21))/70 * f3[k] +
                            (63+13*sqrt(21))/35 * f4[k])
    end
    df(t+h*(7+sqrt(21))/14, ys, f2)
    for k = 1:n
        ys[k] = y[k] + h * ((1/14) * yp[k] +
                            (14+3*sqrt(21))/126 * f3[k] +
                            (13+3*sqrt(21))/63 * f4[k] +
                            (1/9) * f2[k])
    end
    df(t+h*(7+sqrt(21))/14, ys, f1)
    for k = 1:n
        ys[k] = y[k] + h * ((1/32) * yp[k] +
                            (91+21*sqrt(21))/576 * f3[k] +
                            (11/72) * f4[k] +
                            (-385+75*sqrt(21))/1152 * f2[k] +
                            (63-13*sqrt(21))/128 * f1[k])
    end
    df(t+h*(1/2), ys, f5)
    for k = 1:n
        ys[k] = y[k] + h * ((1/14) * yp[k] +
                            (1/9) * f3[k] +
                            (-733+147*sqrt(21))/2205 * f4[k] +
                            (515-111*sqrt(21))/504 * f2[k] +
                            (-51+11*sqrt(21))/56 * f1[k] +
                            (132-28*sqrt(21))/245 * f5[k])
    end
    df(t+h*(7-sqrt(21))/14, ys, f6)
    for k = 1:n
        ys[k] = y[k] + h * ((-42-7*sqrt(21))/18 * f3[k] +
                            (-18-28*sqrt(21))/45 * f4[k] +
                            (-273+53*sqrt(21))/72 * f2[k] +
                            (301-53*sqrt(21))/72 * f1[k] +
                            (28+28*sqrt(21))/45 * f5[k] +
                            (49+7*sqrt(21))/18 * f6[k])
    end
    df(t+h, ys, f2)
    # approximate solution at t+h
    for k = 1:n
        ys[k] = y[k] + h * ((1/20) * yp[k] +
                            (49/180) * f1[k] +
                            (16/45) * f5[k] +
                            (49/180) * f6[k] +
                            (1/20) * f2[k])
    end
end


text:
# Test get/set objective coefficients in LP

using Gurobi
using Base.Test

env = Gurobi.Env()

# original model
#
#   maximize  2x + 2y
#
#	s.t. 0.2 <= x, y <= 1
#

model = gurobi_model(env;
	name="lp_03",
	sense=:maximize,
	f=[2.0, 2.0],
	lb=[0.2, 0.2],
	ub=[1.0, 1.0])

println(model)

lb_ = lowerbounds(model)
ub_ = upperbounds(model)
c_ = objcoeffs(model)

@test lb_ == [0.2, 0.2]
@test ub_ == [1.0, 1.0]
@test c_ == [2.0, 2.0]

optimize(model)

println()
println("soln = $(get_solution(model))")
println("objv = $(get_objval(model))")

# change objective (warm start)
#
#	maximize x - y
#
#	s.t. 0.2 <= x, y <= 1
#

set_objcoeffs!(model, [1, -1])
update_model!(model)

c_ = objcoeffs(model)
@test c_ == [1.0, -1.0]

optimize(model)

println()
println("soln = $(get_solution(model))")
println("objv = $(get_objval(model))")

gc()


text:
using Ramnet.Encoders
using Ramnet

using ..Buffers: MultiStepDynamicBuffer, add!

mutable struct QLearningDiscountedDiscriminatorAgent{O <: AbstractVector,A <: Real,E <: AbstractEncoder} <: AbstractAgent{O,A}
    actions::UnitRange{A}
    n::Int
    size::Int
    regressor_discount::Float64
    rl_discount::Float64
    ϵ::Float64
    encoder::E
    Q̂::Dict{A,RegressionDiscriminator}
    buffer::MultiStepDynamicBuffer{O,A}
    observation::Union{Nothing,O}
    done::Bool
    action::Union{Nothing,A}
    rng::MersenneTwister

    function QLearningDiscountedDiscriminatorAgent{O,A,E}(actions, n, size, regressor_discount, rl_discount, ϵ, encoder::E; seed::Union{Nothing,Int}=nothing) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
        !isnothing(seed) && seed < 0 && throw(DomainError(seed, "Seed must be non-negative"))
        rng = isnothing(seed) ? MersenneTwister() : MersenneTwister(seed)

        n < 1 && throw(DomainError(n, "Tuple size must be at least 1"))
        size < n && throw(DomainError(size, "Input size may not be smaller then tuple size"))
        !(0 ≤ regressor_discount ≤ 1) && throw(DomainError(regressor_discount, "Discount must lie in the [0, 1] interval"))

        Q̂ = Dict(a => RegressionDiscriminator(size, n; γ=regressor_discount) for a in actions)
        buffer = MultiStepDynamicBuffer{O,A}()

        new(
            actions,
            n,
            size,
            regressor_discount,
            rl_discount,
            ϵ,
            encoder,
            Q̂,
            buffer,
            nothing,
            false,
            nothing,
            rng
        )
    end
end

function observe!(agent::QLearningDiscountedDiscriminatorAgent{O,A,E}, observation::O) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    agent.observation = observation
    agent.done = false

    nothing
end

function observe!(agent::QLearningDiscountedDiscriminatorAgent{O,A,E}, action::A, reward::Float64, observation::O, done::Bool) where {O <: AbstractVector,A <: Real,E <: AbstractEncoder}
    add!(agent.buffer, agent.observation, agent.action, reward)

    agent.observation = observation
    agent.done = done

    nothing
end

function update!(agent::QLearningDiscountedDiscriminatorAgent)
    values = similar(agent.actions, Float64)
    q_max = q_values!(agent, agent.observation, values)

    t = popfirst!(agent.buffer)

    train!(
        agent.Q̂[t.action],
        encode(agent.encoder, t.observation),
        t.reward + agent.rl_discount * q_max
    )

    nothing
end

function Agents.reset!(agent::QLearningDiscountedDiscriminatorAgent)
    for action in agent.actions
        agent.Q̂[action] = RegressionDiscriminator(agent.size, agent.n; γ=agent.regressor_discount)
    end

    nothing
end


text:
"""
    pinv_M, rnull_optional, lnull_optional = pinv_qr(M::AbstractArray)

Computes Moore-Penrose pseudoinverse of `M` using a rank-revealing QR
factorization and optionally returns an orthonormal basis for right
and left null-space. This can be significantly faster than using
SVD for computing the pseudoinverse.  Moreover, `A` can be sparse or
dense.

# Examples
```julia-repl
julia> pinv_qr([3 8 17; 2 1 4; 8 3 21])
([-0.10344827586206914 1.344827586206895 -0.17241379310344793;0.1149425287356321 0.8390804597701144 -0.2528735632183907;0.022988505747126506 -0.6321839080459762 0.1494252873563217], Matrix{Float64}(undef, 3, 0), Matrix{Float64}(undef, 3, 0))
```
"""
function pinv_qr(M::AbstractArray{T}) where {T}
    m = size(M, 1)
    n = size(M, 2)

    F = qr(M)
    Q = F.Q*I  # so Q is allways m x m, with Matrix(F.Q) we get a m x n Matrix for m >= n
    R = Matrix(F.R)

    # compute rank
    abs_diag_R = abs.(diag(R))
    tol = 100 * minimum(size(M)) * maximum(abs_diag_R) * eps(T)  # MAGIC CONSTANT
    rnk = sum(tol .< abs_diag_R)

    # compute pseudoinverse
    M_pinv = Matrix{T}(undef, n, m)

    S = R[1:rnk, 1:rnk] \ R[1:rnk, rnk+1:end]
    W = R[1:rnk, 1:rnk] \ Q[:, 1:rnk]'
    X = (S' * S + I) \ (S' * W)
    M_pinv[rnk+1:end, :] = X
    M_pinv[1:rnk, :] = W - S*X

    # computing null spaces
    right_null_size = n - rnk
    left_null_size = m - rnk

    # right null space
    if right_null_size > 0
        rnull, = qr(Matrix([S; -I(right_null_size)]))
        rnull = Matrix(rnull)
    else
        rnull = zeros(T, n, 0)
    end

    # left null space
    if left_null_size > 0
        lnull = Q[: ,end - (left_null_size - 1): end]
    else
        lnull = zeros(T, m, 0)
    end

    return M_pinv, rnull, lnull
end

function pinv_qr(M::AbstractArray{<:Int})
    pinv_qr(float(M))
end


text:
module test_20_SeveralSignalsInOneDiagram

using ModiaPlot

include("test_base/test_20_SeveralSignalsInOneDiagram.jl")

end


text:
function DA_algo(m_pref::Array{Int64}, n_pref::Array{Int64}) #引数は選考表
    m = size(m_pref,2)
    n = size(n_pref,2)
    m_match = zeros(Int64, m)
    n_match = zeros(Int64, n)   #マッチング結果vectorの作成
    propose_number = zeros(Int64,m)  #告白回数の記憶vecotr
    no_propose = Int[]  #それぞれが「もう告白しなくなる」告白回数vector
    for man in 1:m
        push!(no_propose, find(m_pref[:, man] .== 0)[1] - 1)
    end
    while any(m_match .== 0) == true
        mikon = find(m_match .== 0)  #未婚menVector
        if all(propose_number[mikon] .>= no_propose[mikon]) == true
            break  #未婚menが「告白しない」状態なら脱loop
        else
            for each_mikon in mikon
                proposed = m_pref[propose_number[each_mikon] + 1, each_mikon]  #次のプロポーズ相手
                if proposed != 0
                    propose_number[each_mikon] += 1 #告白回数に+1
                    if find(n_pref[:,proposed] .== each_mikon)[1] < find(n_pref[:,proposed] .== n_match[proposed])[1]  #自分が相手の婚約者（いないなら、０）より、マシであれば成婚
                        if n_match[proposed] != 0           #婚約者がいれば、離婚します
                            m_match[n_match[proposed]] = 0
                        end
                    m_match[each_mikon] = proposed
                    n_match[proposed] = each_mikon
                    end
                end
            end
        end
    end
    return m_match, n_match
end

    function pref(m,n) #ランダム選考表の作成(今回は使わない)
        m_pref = shuffle!(collect(0:n))
        n_pref = shuffle!(collect(0:m))
        #1列目だけ作成しておく
        for i in　2:m
            each = shuffle!(collect(0:n))
            m_pref = hcat(m_pref, each)
        end
        for i in 2:n
            each = shuffle!(collect(0:m))
            n_pref = hcat(n_pref,each)
        end
        #あとはそれぞれloop
        return m_pref, n_pref
    end

max_stars_repo_name: M-okb/DA-Matching

max_issues_repo_name: M-okb/DA-Matching

max_forks_repo_name: M-okb/DA-Matching


text:
# Non Differentiable Functions
ChainRulesCore.@non_differentiable replicate(::Any)
ChainRulesCore.@non_differentiable update_statistics(::Any, ::Any, ::Any, ::Any, ::Any,
                                                     ::Any, ::Any)
ChainRulesCore.@non_differentiable generate_dropout_mask(::Any, ::Any, ::Any)
ChainRulesCore.@non_differentiable compute_adaptive_pooling_dims(::Any, ::Any)
ChainRulesCore.@non_differentiable glorot_normal(::Any...)
ChainRulesCore.@non_differentiable glorot_uniform(::Any...)
ChainRulesCore.@non_differentiable check_use_cuda()
ChainRulesCore.@non_differentiable istraining(::Any)

ChainRulesCore.Tangent{P}(; kwargs...) where {P <: AbstractExplicitLayer} = NoTangent()

ChainRulesCore.rrule(::typeof(istraining)) = true, _ -> (NoTangent(),)

function ChainRulesCore.rrule(::typeof(Base.broadcasted), ::typeof(identity), x)
    x, Δ -> (NoTangent(), NoTangent(), Δ)
end

# NNlib Functions
function ChainRulesCore.rrule(::typeof(batchnorm),
                              g::CuArray{T},
                              b::CuArray{T},
                              x::Union{CuArray{T, 4}, CuArray{T, 5}},
                              running_mean,
                              running_var,
                              momentum;
                              kwargs...) where {T <: CUDNNFloat}
    y = batchnorm(g, b, x, running_mean, running_var, momentum; kwargs...)
    function batchnorm_pullback(dy)
        dg, db, dx = ∇batchnorm(g, b, x, dy, running_mean, running_var, momentum; kwargs...)
        return NoTangent(), dg, db, dx, NoTangent(), NoTangent(), NoTangent()
    end
    return y, batchnorm_pullback
end

# Activation Rrules
function ChainRulesCore.rrule(::typeof(applyactivation), f::cudnnValidActivationTypes,
                              x::CuArray{T}) where {T <: CUDNNFloat}
    mode = getCUDNNActivationMode(f)
    y = CUDNN.cudnnActivationForward(x; mode)
    function applyactivation_pullback(Δ)
        return NoTangent(), NoTangent(), cudnnActivationBackward(y, Δ, x; mode), NoTangent()
    end
    return y, applyactivation_pullback
end

# Elementwise Functions
function ChainRulesCore.rrule(::typeof(elementwise_add), x, y) where {T}
    z = elementwise_add(x, y)
    _elementwise_add_pullback(Δ) = (NoTangent(), elementwise_add_pullback(x, y, Δ)...)
    return z, _elementwise_add_pullback
end

function ChainRulesCore.rrule(::typeof(elementwise_mul), x, y) where {T}
    z = elementwise_mul(x, y)
    _elementwise_mul_pullback(Δ) = (NoTangent(), elementwise_mul_pullback(x, y, Δ)...)
    return z, _elementwise_mul_pullback
end

# Zygote Fixes
function Zygote.accum(x::ComponentArray, ys::ComponentArray...)
    return ComponentArray(Zygote.accum(getdata(x), getdata.(ys)...), getaxes(x))
end

# Adapt Interface
function ChainRulesCore.rrule(::typeof(Array), x::CUDA.CuArray)
    return Array(x), d -> (NoTangent(), CUDA.cu(d))
end

function ChainRulesCore.rrule(::typeof(adapt_storage), to::LuxCPUAdaptor,
                              x::CUDA.AbstractGPUArray)
    return adapt_storage(to, x),
           d -> (NoTangent(), NoTangent(), adapt_storage(LuxCUDAAdaptor(), d))
end

function ChainRulesCore.rrule(::typeof(adapt_storage), to::LuxCUDAAdaptor, x::Array)
    return adapt_storage(to, x),
           d -> (NoTangent(), NoTangent(), adapt_storage(LuxCPUAdaptor(), d))
end

# RNN Helpers
## Taken from https://github.com/FluxML/Flux.jl/blob/1f82da4bfa051c809f7f3ce7dd7aeb43be515b14/src/layers/recurrent.jl#L9
function ChainRulesCore.rrule(::typeof(multigate), x::AbstractArray, c::Val{N}) where {N}
    function multigate_pullback(dy)
        dx = map!(zero, similar(x, float(eltype(x)), axes(x)), x)
        foreach(multigate(dx, c), dy) do dxᵢ, dyᵢ
            dyᵢ isa AbstractZero && return
            @. dxᵢ += dyᵢ
        end
        return (NoTangent(), dx, NoTangent(), NoTangent())
    end
    return multigate(x, c), multigate_pullback
end


text:
module Utils
using WAV

function window(n::Int, sr::Int, fn::String)
    w = wavread(fn)[1]
    c = w[1:20*sr, 1]
    len = size(c)[1] - (2n + 1)
    m = div(len, n)
    idxs = 1:2n:m*n
    xs = hcat([vec(c[i:i+n-1]) for i in idxs]...)
    ys = hcat([vec(c[i+n:i+2n-1]) for i in idxs]...)
    return xs, ys
end

function window_multiple(n::Int, sr::Int, fns::Vector{String})
    all_xs = Matrix(undef, n, 0)
    all_ys = Matrix(undef, n, 0)
    for f in fns
        xs, ys = window(n, sr, f)
        all_xs = hcat(all_xs, xs)
        all_ys = hcat(all_ys, ys)
    end
    return all_xs, all_ys
end
end

max_stars_repo_name: anandijain/audio

max_issues_repo_name: anandijain/audio

max_forks_repo_name: anandijain/audio


text:
function _test_cov_vec(ce::CovarianceEstimator, X::AbstractMatrix{<:Real})
    c = cov(ce, X)
    @inbounds for i in 1:size(X, 2)
        @test cov(ce, @view(X[:, i])) ≈ c[i, i] atol = 1e-12 rtol = 1e-16
    end
end

function _test_cov_vecvec(ce::CovarianceEstimator, X::AbstractMatrix{<:Real})
    c = cov(ce, X)
    @inbounds for i in 1:size(X, 2)
        for j in 1:size(X, 2)
            cij = cov(ce, @view(X[:, i]), @view(X[:, j]))
            @test cij ≈ c[i, j] atol = 1e-12 rtol = 1e-16
        end
    end
end

@testset "BiweightMidcovariance" begin
    for c in (9, 5)
        for modify_sample_size in (true, false)
            @testset "basic properties c=$c, modify_sample_size=$modify_sample_size" begin
                ce = BiweightMidcovariance(; c=c, modify_sample_size=modify_sample_size)
                for X in test_matrices
                    testTransposition(ce, X)
                    testDims(ce, X)
                    testUncorrelated(ce)
                    testTranslation(ce, X)

                    # Test that other dispatches are consistent with the one tested above.
                    _test_cov_vec(ce, X)
                    _test_cov_vecvec(ce, X)
                end
            end
        end
    end

    @testset "Z should give zeros, not NaNs" begin
        # This covers a difference discovered against astropy.
        # It was due to not implementing the special case that covariance is defined to be 0
        # when MAD = 0.
        @test cov(BiweightMidcovariance(), Z) == zeros(Float64, 3, 3)
    end

    @testset "astropy reference (c=9.0)" begin
        # These tests are against references generated with astropy.
        #
        # References generated with:
        #   * System: Linux 5.4.0-90-generic, x86_64, Intel 8700K
        #   * Software: Python 3.7.6, astropy 4.3.1, numpy 1.21.2
        #
        # import os
        # import numpy
        # from astropy.stats import biweight_midcovariance
        # prefix = os.path.join(
        #     os.path.expanduser("~"), ".julia", "dev", "CovarianceEstimation", "test",
        #     "test_matrices"
        # )
        # for sample in ["20x100", "50x50", "100x20", "50x50tdist3"]:
        #     for modify_sample_size in [False, True]:
        #         path = os.path.join(prefix, f"{sample}.csv")
        #         X = numpy.genfromtxt(path).T
        #         c = biweight_midcovariance(X, modify_sample_size=modify_sample_size)
        #         mss = "_mss" if modify_sample_size else ""
        #         out_path = os.path.join(prefix, f"{sample}_bwcov{mss}.csv")
        #         with open(out_path, "w") as file_:
        #             for i in range(c.shape[0]):
        #                 file_.write(" ".join(str(el) for el in c[i, :]))
        #                 file_.write("\n")
        names = ["20x100", "50x50", "100x20", "50x50tdist3"]
        _test_refs(BiweightMidcovariance(), names, "bwcov")
        _test_refs(BiweightMidcovariance(; modify_sample_size=true), names, "bwcov_mss")
    end
end


text:
# f(x)=A*X'
function FAXtranspose(A,X)
    if size(A)==(1,1)
        return (A[1].*X',nothing)
    elseif size(X)==(1,1)
        return (A.*X[1],nothing)
    else
        return (A*X',nothing)
    end
end

function FAXtranspose_inplace(value,auxvalue,A,X)
    if size(A)==(1,1)
        copy!(value,A[1]*X')
    elseif size(X)==(1,1)
        copy!(value,A*X[1])
    else
        BLAS.gemm!('N','T',1.0,A,X,0.0,value)
    end
end

function DAXtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,A,X)
    if derivativeIDX==1
        if size(A)==(1,1)
            axpy!(1.0,[sum(X'.*grad_c)],grad_n)
        elseif size(X)==(1,1)
            axpy!(X[1],grad_c,grad_n)
        else
            BLAS.gemm!('N','N',1.0,grad_c,X,1.0,grad_n) # grad_c*X
        end
    elseif derivativeIDX==2
        if size(A)==(1,1)
            axpy!(A[1],grad_c',grad_n)
        elseif size(X)==(1,1)
            axpy!(1.0,[sum(A.*grad_c)],grad_n)
        else
            BLAS.gemm!('T','N',1.0,grad_c,A,1.0,grad_n) # grad_c'*A
        end
    end
end

if PROC=="GPU"
#    function FAXtranspose(A::CudaArray,X::CudaArray)
#        if size(A)==(1,1)
#            tmp=CudaArray(Float64,size(X)); gfill!(tmp,A)
#            return (CUBLAS.gemm('N','T',tmp,X),nothing)
#        elseif size(X)==(1,1)
#            tmp=CudaArray(Float64,size(A)); gfill!(tmp,X)
#            return (CUBLAS.gemm('N','N',A,tmp),nothing)
#        else
#            return (CUBLAS.gemm('N','T',A,X),[])
#        end
#    end

    function FAXtranspose_inplace(value::CudaArray,auxvalue,A::CudaArray,X::CudaArray)
        if size(A)==(1,1)
            CUBLAS.geam!('T','N',1.0,X,0.0,value,value) # nb argument converse of Base.scale!
            scale!(A,value) # nb argument converse of Base.scale!
        elseif size(X)==(1,1)
            copy!(value,A); scale!(X,value) # nb argument converse of Base.scale!
        else
            CUBLAS.gemm!('N','T',1.0,A,X,0.0,value)
        end
    end

    function DAXtranspose(derivativeIDX,f_c,faux_c,grad_c,grad_n,A::CudaArray,X::CudaArray)
        if derivativeIDX==1
            if size(A)==(1,1)
                tmp=CudaArray(Float64,size(X))
                CUBLAS.geam!('T','N',1.0,X,0.0,tmp,tmp) # tmp=X'
                vmult!(1.0,tmp,grad_c,tmp) # tmp=X'.*grad_c
                sum_update!(1.0,tmp,1.0,grad_n) # grad_n+=sum(X'.*grad_c)
                free(tmp)
            elseif size(X)==(1,1)
                alphaaxpy!(1.0,X,grad_c,grad_n)
            else
                CUBLAS.gemm!('N','N',1.0,grad_c,X,1.0,grad_n) # grad_c*X
            end
        elseif derivativeIDX==2
            if size(A)==(1,1)
                tmp=CudaArray(Float64,size(grad_c))
                CUBLAS.geam!('T','N',1.0,grad_c,0.0,tmp,tmp) # tmp=grad_c'
                alphaaxpy!(1.0,A,tmp,grad_n); free(tmp)
            elseif size(X)==(1,1)
                tmp=CudaArray(Float64,size(A))
                vmult!(1.0,A,grad_c,tmp)
                sum_update!(1.0,tmp,1.0,grad_n) # sum(A.*grad_c)
                free(tmp)
            else
                CUBLAS.gemm!('T','N',1.0,grad_c,A,1.0,grad_n) # grad_c'*A
            end
        end
    end
end

Derivative[FAXtranspose]=DAXtranspose
export FAXtranspose

Inplace[FAXtranspose]=FAXtranspose_inplace

max_stars_repo_name: davidbarber/AutoDiff

max_issues_repo_name: davidbarber/AutoDiff

max_forks_repo_name: davidbarber/AutoDiff


text:
using Test
#using Arpack: eigs
using ExactDiagonalization.FEDW: NCBasis,bkron,boplus
using QuantumLattices
using LinearAlgebra: eigen,diagm
using ExactDiagonalization.Perturbation: ProjState,GSpinGenerator,jval,hamiltonianeff,SpinHeff,projstate

lattice=Lattice("L2P",[Point(PID(1),(0.0,)),Point(PID(2),(1.0,))])
config=Config{Fock{:f}}(pid->Fock{:f}(atom=pid.site%2, norbital=1, nspin=2, nnambu=2), lattice.pids)
table=Table(config,usualfockindextotuple)
t=Hopping(:t,1.0,1)
U=Hubbard(:U,4.0)
μ=Onsite(:μ,-0.5)
h0=Generator((U,μ),Bonds(lattice),config,half=false,table=table)
h1 =Generator((t,),Bonds(lattice),config,half=false,table=table)

@testset "projstate" begin
    basislowi=NCBasis{UInt}([2],[[(1,)]])
    pjslowi=projstate(expand(h0),basislowi,h0.table,pick=1:2)
    @test isa(pjslowi,ProjState)
    @test pjslowi.vals == [-0.50;-0.50]
    #show( pjslowi.vecs )
    @test pjslowi.basis.table ==  basislowi.table
    pjslowij = bkron(pjslowi,pjslowi)

    @test pjslowij.vals == [-1.0;-1.0;-1.0;-1.0]
    @test pjslowij.vecs ==  kron(pjslowi.vecs,pjslowi.vecs)
    #show(pjslowij.basis)

    basis_high_i = NCBasis{UInt}([2],[[(2,)]])
    basis_high_j = NCBasis{UInt}([2],[[(0,)]])

    pjs_high_i = projstate(expand(h0),basis_high_i,h0.table)
    pjs_high_j = projstate(expand(h0),basis_high_j,h0.table)
    pjs_high_ij = bkron(pjs_high_i,pjs_high_j)
    pjs_high_ji = bkron(pjs_high_j,pjs_high_i)
    pjs_high = boplus(pjs_high_ij,pjs_high_ji)

    basis_high = boplus( bkron(basis_high_i,basis_high_j),bkron(basis_high_j,basis_high_i) )
   # show(basis_high)
    pjs_high2=projstate(expand(h0),basis_high,h0.table)
    @test pjs_high.vals ==pjs_high2.vals
    #@test pjs_high.vecs == pjs_high2.vecs
    @test pjs_high.basis.table == pjs_high2.basis.table
    #show(pjs_high.basis)
end

@testset "GSpinGenerator and SpinHeff" begin
    namestri=["σᵢ⁰","σᵢˣ","σᵢʸ","σᵢᶻ"]
    namsi = tuple([Symbol(nm) for nm in namestri]...)
    σ0=[1.0 0.0;0.0 1.0]
    σx=[0.0 1.0; 1.0 0.0]
    σy=[0.0 -im;im 0.0]
    σz=[1.0 0.0; 0.0 -1.0]
    gsgi=GSpinGenerator(namsi,(σ0,σx,σy,σz))
    @test gsgi[:σᵢˣ]==σx
    @test gsgi[3]==(:σᵢʸ,σy)
    @test keys(gsgi) == namsi
    @test eltype(gsgi) == ComplexF64
    @test length(gsgi)  == 4
    @test dimension(gsgi) ==(2,2)
    namestrj=["σⱼ⁰","σⱼˣ","σⱼʸ","σⱼᶻ"]
    namsj = tuple([Symbol(nm) for nm in namestrj]...)
    gsgj = GSpinGenerator(namsj,(σ0,σx,σy,σz))
    gsgij=kron(gsgi,gsgj)
    @test gsgij.names==tuple([Symbol(i,",",j) for i in namestri for j in namestrj]...)
    @test gsgij[Symbol("σᵢˣ,σⱼʸ")]==kron(σx,σy)

    basislowij=NCBasis{UInt}([2,2],[[(1,),(1,)]])
    pjslowij=projstate(expand(h0),basislowij,h0.table,pick=1:4)
    basishighij = NCBasis{UInt}([2,2],[[(2,),(0,)],[(0,),(2,)]])
    pjshighij = projstate(expand(h0),basishighij,h0.table)
    heffterm = hamiltonianeff(pjslowij, pjshighij,expand(h1),h1.table)
    jex=jval(heffterm,gsgij)
    @test isa(jex,SpinHeff)
    @test jex[1] == diagm([-1.0+0im,-1.0,-1.0,-1.0])
    @test jex.jeff[[1,6,11,16]]==[-1.25,0.25,0.25,0.25]
    @test keys(jex) == (Symbol("σᵢ⁰,σⱼ⁰"), Symbol("σᵢ⁰,σⱼˣ"), Symbol("σᵢ⁰,σⱼʸ"), Symbol("σᵢ⁰,σⱼᶻ"), Symbol("σᵢˣ,σⱼ⁰"), Symbol("σᵢˣ,σⱼˣ"), Symbol("σᵢˣ,σⱼʸ"), Symbol("σᵢˣ,σⱼᶻ"), Symbol("σᵢʸ,σⱼ⁰"), Symbol("σᵢʸ,σⱼˣ"), Symbol("σᵢʸ,σⱼʸ"), Symbol("σᵢʸ,σⱼᶻ"), Symbol("σᵢᶻ,σⱼ⁰"), Symbol("σᵢᶻ,σⱼˣ"), Symbol("σᵢᶻ,σⱼʸ"), Symbol("σᵢᶻ,σⱼᶻ"))
    @test values(jex) == jex.jeff
    @test jex[Symbol("σᵢ⁰,σⱼ⁰")] == -1.25
    #show(jex)
    #===============================================================#
    #second method to get ProjState

    pjslowi2 = projstate(expand(h0),NCBasis{UInt}([2],[[(1,)]]),h0.table,pick=1:2)
    pjshighi2 = projstate(expand(h0),NCBasis{UInt}([2],[[(2,)]]),h0.table)
    pjshighi0 = projstate(expand(h0),NCBasis{UInt}([2],[[(0,)]]),h0.table)

    pjslowj2 = projstate(expand(h0),NCBasis{UInt}([2],[[(1,)]]),h0.table,pick=1:2)
    pjshighj2 = projstate(expand(h0),NCBasis{UInt}([2],[[(2,)]]),h0.table)
    pjshighj0 = projstate(expand(h0),NCBasis{UInt}([2],[[(0,)]]),h0.table)

    pjslowij2 = bkron(pjslowi2,pjslowj2)
    pjshighij2 = boplus(bkron(pjshighi2,pjshighj0),bkron(pjshighi0,pjshighj2))

    heffterm2 = hamiltonianeff(pjslowij2, pjshighij2,expand(h1),h1.table)
    jex2=jval(heffterm2,gsgij)
    @test jex2[1] == diagm([-1.0+0im,-1.0,-1.0,-1.0])
    @test jex2.jeff[[1,6,11,16]]==[-1.25,0.25,0.25,0.25]
end


text:
__precompile__()

module ImagineInterface

using JSON, Unitful, UnitAliases, Random, Mmap
using MappedArrays, IntervalSets, DataStructures, DSP
using ImagineFormat, ImagineHardware
import ImagineHardware:samprate
using AxisArrays
const axes = Base.axes

import Base: convert, show, length, size, isempty, ==, append!, pop!, empty!, replace!#, scale

using Unitful: μm, s, Hz, V

include("metadata_constants.jl")
#Load hardware parameters for all rigs
rig_dir = joinpath(dirname(@__DIR__), "rigs")
dictmap(f, d::T) where T = T((f(p) for p = pairs(d))) #used in rig files
for rig_file in readdir(rig_dir)
    include(joinpath(rig_dir, rig_file))
end
include("samplemapper.jl")
include("imaginesignal.jl")
include("hardware_templates.jl")
include("convenience.jl")
include("parse.jl")
include("sequence_analysis.jl")
include("validate_group.jl")
include("window_validate.jl")
include("validate_single.jl")
include("write.jl")
include("stack.jl")

#metadata_constants.jl
export chip_size,
        max_framerate,
        max_roi

#imaginecommand.jl
export RepeatedValue,
        RLEVector,
        full_length,
        ImagineSignal,
        name,
        rename!,
        duration,
        daq_channel,
        rig_name,
	rawtype,
        worldtype,
        samprate,
        intervals,
        interval_raw,
        interval_volts,
        interval_world,
        sequences,
	sequence_names,
	sequence_lookup,
        mapper,
        get_samples,
        add_sequence!,
        replace!,
        replicate!

#hardware_templates.jl
export rigtemplate

#convenience.jl
export  isfree,
        isdigital,
        isanalog,
        isoutput,
        isinput,
        ispos,
        isposmonitor,
        iscam,
        iscammonitor,
        islas,
        isstim,
        getname,
        getdigital,
        getanalog,
        getinputs,
        getoutputs,
        getfree,
        getfixed,
        getpositioners,
        getpositionermonitors,
        getcameras,
        getcameramonitors,
        getlasers,
        getstimuli,
        getgalvos,
        getgalvomonitors,
        hasmonitor,
        hasactuator,
        monitor_name,
        actuator_name,
        finddigital,
        findanalog,
        findinputs,
        findoutputs,
        findfree,
        findfixed,
        findpositioners,
        findcameras,
        findlasers,
        findstimuli

#parse.jl
export parse_command,
        parse_commands,
        parse_ai,
        parse_di,
        load_signals

#sequence_analysis.jl
export find_pulse_starts,
        find_pulse_stops,
        count_pulses

#validate_singles.jl
export validate_all

#write.jl
export write_commands

#stack.jl
export  apply_lp,
        gen_sweep,
        gen_sawtooth,
        gen_bidi_pos,
        spaced_intervals,
        gen_pulses,
        scale,
        gen_bidirectional_stack,
        gen_unidirectional_stack,
        gen_stepped_stack,
        gen_2d_timeseries
#deprecations
@deprecate ImagineCommand ImagineSignal
@deprecate decompress(args...;kwargs...) get_samples(args...;kwargs...)

end #module

max_forks_repo_name: HolyLab/ImagineInterface


text:
"""
$(TYPEDSIGNATURES)

Takes a Nth order ODESystem and returns a new ODESystem written in first order
form by defining new variables which represent the N-1 derivatives.
"""
function ode_order_lowering(sys::ODESystem)
    iv = get_iv(sys)
    eqs_lowered, new_vars = ode_order_lowering(equations(sys), iv, states(sys))
    @set! sys.eqs = eqs_lowered
    @set! sys.states = new_vars
    return sys
end

function dae_order_lowering(sys::ODESystem)
    iv = get_iv(sys)
    eqs_lowered, new_vars = dae_order_lowering(equations(sys), iv, states(sys))
    @set! sys.eqs = eqs_lowered
    @set! sys.states = new_vars
    return sys
end

function ode_order_lowering(eqs, iv, states)
    var_order = OrderedDict{Any, Int}()
    D = Differential(iv)
    diff_eqs = Equation[]
    diff_vars = []
    alge_eqs = Equation[]

    for (i, eq) in enumerate(eqs)
        if !isdiffeq(eq)
            push!(alge_eqs, eq)
        else
            var, maxorder = var_from_nested_derivative(eq.lhs)
            maxorder > get(var_order, var, 1) && (var_order[var] = maxorder)
            var′ = lower_varname(var, iv, maxorder - 1)
            rhs′ = diff2term(eq.rhs)
            push!(diff_vars, var′)
            push!(diff_eqs, D(var′) ~ rhs′)
        end
    end

    for (var, order) in var_order
        for o in (order - 1):-1:1
            lvar = lower_varname(var, iv, o - 1)
            rvar = lower_varname(var, iv, o)
            push!(diff_vars, lvar)

            rhs = rvar
            eq = Differential(iv)(lvar) ~ rhs
            push!(diff_eqs, eq)
        end
    end

    # we want to order the equations and variables to be `(diff, alge)`
    return (vcat(diff_eqs, alge_eqs), vcat(diff_vars, setdiff(states, diff_vars)))
end

function dae_order_lowering(eqs, iv, states)
    var_order = OrderedDict{Any, Int}()
    D = Differential(iv)
    diff_eqs = Equation[]
    diff_vars = OrderedSet()
    alge_eqs = Equation[]
    vars = Set()
    subs = Dict()

    for (i, eq) in enumerate(eqs)
        vars!(vars, eq)
        n_diffvars = 0
        for vv in vars
            isdifferential(vv) || continue
            var, maxorder = var_from_nested_derivative(vv)
            isparameter(var) && continue
            n_diffvars += 1
            order = get(var_order, var, nothing)
            seen = order !== nothing
            if !seen
                order = 1
            end
            maxorder > order && (var_order[var] = maxorder)
            var′ = lower_varname(var, iv, maxorder - 1)
            subs[vv] = D(var′)
            if !seen
                push!(diff_vars, var′)
            end
        end
        n_diffvars == 0 && push!(alge_eqs, eq)
        empty!(vars)
    end

    for (var, order) in var_order
        for o in (order - 1):-1:1
            lvar = lower_varname(var, iv, o - 1)
            rvar = lower_varname(var, iv, o)
            push!(diff_vars, lvar)

            rhs = rvar
            eq = Differential(iv)(lvar) ~ rhs
            push!(diff_eqs, eq)
        end
    end

    return ([diff_eqs; substitute.(eqs, (subs,))],
            vcat(collect(diff_vars), setdiff(states, diff_vars)))
end


text:
using DataAPI, PooledArrays

function compress_then_write(b::PooledVector, io)
    compress_refs = compress_then_write(b.refs, io)
    compress_poolindex = compress_then_write(b.pool, io)

    (type = PooledVector, refs = compress_refs, poolindex = compress_poolindex)
end

function column_loader(b::Type{PooledVector}, io, metadata)
    refs_meta = metadata.refs
    pi_meta = metadata.poolindex
    ref = column_loader(refs_meta.type, io, refs_meta)
    poolindex = column_loader(pi_meta.type, io, pi_meta)

    # TODO more efficient pooledArray construction
    # PooledVector{pi_meta.type}(ref, CategoricalPool{eltype(poolindex), eltype(ref)}(Array(poolindex)))

    ref = max.(1, min.(length(poolindex), ref))
    PooledArray(poolindex[ref])
end

max_stars_repo_name: xiaodaigh/JDF

max_issues_repo_name: xiaodaigh/JDF

max_forks_repo_name: xiaodaigh/JDF