text:
struct Config
    n_nodes::Int64
    n0::Int64
    k::Int64
    agentcount::Int64
    transmission_prob::Float64
    test_prob::Float64
    quarantine_duration::Int64
    behavior_reduction_factor::Float64
    behavior_activated::Bool
    tick_authority_recommendation::Int64
    tick_authority_policy::Int64
    n_public_spaces::Int64
    initial_infectives::Int64
    n_iter::Int64
end

Config(;
    n_nodes::Int64, n0::Int64, k::Int64, agentcount::Int64,
    transmission_prob::Float64, test_prob::Float64, quarantine_duration::Int64,
    behavior_reduction_factor::Float64, behavior_activated::Bool,
    tick_authority_recommendation::Int64, tick_authority_policy::Int64,
    n_public_spaces::Int64, initial_infectives::Int64, n_iter::Int64
) = Config(
    n_nodes, n0, k, agentcount,
    transmission_prob, test_prob, quarantine_duration,
    behavior_reduction_factor, behavior_activated,
    tick_authority_recommendation, tick_authority_policy,
    n_public_spaces, initial_infectives, n_iter
)

Config() = Config(
    n_nodes=1_000, n0=3, k=3, agentcount=100,
    transmission_prob=0.1, test_prob=0.6, quarantine_duration=140,
    behavior_reduction_factor=0.6, behavior_activated=true,
    tick_authority_recommendation=10, tick_authority_policy=50,
    n_public_spaces=30, initial_infectives=1, n_iter=1_000
)


text:
include("theislike_piml_setup.jl")

# Create generator of (wind,distance) tuples
data_train_batch = [[sample() for i = 1:1] for v in 1:10]
data_test = [[sample()] for i = 1:10]

# Place to save stuff
#folder_name = "AD_results_new"
epochs = 1:10
for epoch in epochs
	println("hi")
	#tt = @elapsed Flux.train!(loss, θ, data_train_batch, opt)
	tt = @elapsed Flux.train!(loss, θ, data_train_batch, opt, cb = cb)
	push!(train_time, tt)
	loss_train = sum(map(x->loss(x,train=false),data_train_batch))
	loss_test = sum(map(x->loss(x,train=false),data_test))
	rmse_train = sqrt(loss_train/length(data_train_batch))
	rmse_test = sqrt(loss_test/length(data_test))
	diffs_train = map(x->mydiff(x,train=false),data_train_batch)
	diffs_test = map(x->mydiff(x,train=false),data_test)
	meandiff_train = mean(diffs_train)
	meandiff_test = mean(diffs_test)
	stddiff_train = std(diffs_train)
	stddiff_test = std(diffs_test)
	# Terminal output
	println(string("epoch: ", epoch, " train rmse: ", rmse_train, " test rmse: ", rmse_test))
	# Save convergence metrics
	push!(losses_train, loss_train)
	push!(losses_test, loss_test)
	push!(rmses_train, rmse_train)
	push!(rmses_test, rmse_test)
	push!(meandiffs_train, meandiff_train)
	push!(meandiffs_test, meandiff_test)
	push!(stddiffs_train, stddiff_train)
	push!(stddiffs_test, stddiff_test)
	#plot_maps(epoch,display=false,folder_name=folder_name);
	#run(`convert $[string(folder_name,"/images/iters",lpad(epoch,4,"0"),".png"), string(folder_name,"/images/convergence",lpad(epoch,4,"0"),".png")] -append $[string(folder_name,"/images/combine",lpad(epoch,4,"0"),".png")]`)
end

#plot_losses(folder_name=folder_name, combine=true)

#@save string(folder_name,"/data/loss_data.jld2") meandiffs_train stddiffs_train epochs losses_test losses_train meandiffs_test stddiffs_test rmses_train rmses_test train_time θ well_crds mon_well_crds

#run(`ffmpeg -y -r 5 -f image2 -s 1920x1080 -i $folder_name/images/combine%04d.png -vcodec libx264 -crf 25  -pix_fmt yuv420p $folder_name/multitheisdp_1well.mp4`)

ext: jl


text:
module StateTest
using Base.Test
using Quantum

xb = basis(:X,[1:5]);
yb = basis(:Y,[1:5]);
xyb = kron(xb,yb);

@test 0*xb[1]==dvec([0], basis(xb[1]))
@test 1*xb[1]==dvec([1], basis(xb[1]))
@test xb[1]'*xb[1]==1
@test isdual(xb[1]',xb[1])
@test sum(xb[:])==DiracVector([1,1,1,1,1], xb)
@test xb[1]'*xyb[1]==1*yb[1]
@test separate(xyb[1])==xyb[1][:]
@test reduce(+,[xb[i]*yb[j]' for i=1:5, j=1:5]) == DiracMatrix(ones(5,5), xb, yb')
@test inner(xyb[1]', xyb[1], 1)==1
@test inner(xyb[1]', xyb[1], 2)==xb[1]'*yb[1]*yb[1]'*xb[1]
@test inner(xyb[1]', xyb[2], 1)==0
@test inner(xyb[1]', xyb[2], 2)==xb[1]'*yb[2]*yb[1]'*xb[1]

k = ket(:K,1)
s = ket(:S,"a")
v = [tensor(s,s),tensor(s,k),tensor(s',s'), tensor(s',k')]

@test typeof(v[1]) == Tensor{Ket{:S,ASCIIString}}
@test typeof(v[2]) == Tensor{Ket}
@test typeof(v[3]) == Tensor{Bra{:S,ASCIIString}}
@test typeof(v[4]) == Tensor{Bra}
@test typeof(tensor(v[1], s)) == Tensor{Ket{:S,ASCIIString}}
@test typeof(tensor(s, v[1])) == Tensor{Ket{:S,ASCIIString}}
@test typeof(tensor(v[2], s)) == Tensor{Ket}
@test typeof(tensor(s, v[2])) == Tensor{Ket}
@test typeof(tensor(v[3], s')) == Tensor{Bra{:S,ASCIIString}}
@test typeof(tensor(s', v[3])) == Tensor{Bra{:S,ASCIIString}}
@test typeof(tensor(v[4], s')) == Tensor{Bra}
@test typeof(tensor(s', v[4])) == Tensor{Bra}

@test typeof(tensor(v[1], v[1])) == Tensor{Ket{:S,ASCIIString}}
@test typeof(tensor(v[2], v[1])) == Tensor{Ket}
@test typeof(tensor(v[3], v[3])) == Tensor{Bra{:S,ASCIIString}}
@test typeof(tensor(v[4], v[3])) == Tensor{Bra}

end


text:
using ImageCore, Colors, FixedPointNumbers, OffsetArrays
using Test, Random

@testset "permuteddimsview" begin
    a = [1 3; 2 4]
    v = permuteddimsview(a, (1,2))
    @test v == a
    v = permuteddimsview(a, (2,1))
    @test v == a'
    a = rand(3,7,5)
    v = permuteddimsview(a, (2,3,1))
    @test v == permutedims(a, (2,3,1))
end

@testset "convert (deprecations)" begin
    @info "Deprecation warnings are expected"
    a = [RGB(1,0,0) RGB(0,0,1);
         RGB(0,1,0) RGB(1,1,1)]
    c = @inferred(convert(Array{BGR}, a))
    @test eltype(c) == BGR{N0f8}
    c = @inferred(convert(Array{BGR{Float32}}, a))
    @test eltype(c) == BGR{Float32}
    c = @inferred(convert(Array{Lab}, a))
    @test eltype(c) == Lab{Float32}
    for a in (rand(Float32, (4,5)),
              bitrand(4,5))
        b = @inferred(convert(Array{Gray}, a))
        @test eltype(b) == Gray{eltype(a)}
        b = @inferred(convert(Array{Gray{N0f8}}, a))
        @test eltype(b) == Gray{N0f8}
    end

    # Gray images wrapped by an OffsetArray.
    A = rand(8,8)
    for img in ( Gray.(A),
                 Gray.(N0f8.(A)),
                 Gray.(N0f16.(A)) )
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{Gray{Float32},2,Array{Gray{Float32}}},imgo))
        @test eltype(s) == Gray{Float32}
        @test s isa OffsetArray{Gray{Float32},2,Array{Gray{Float32},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    for img in ( Gray.(A),
                 Gray.(N0f8.(A)),
                 Gray.(N0f16.(A)) )
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{Gray{N0f8},2,Array{Gray{N0f8}}},imgo))
        @test eltype(s) == Gray{N0f8}
        @test s isa OffsetArray{Gray{N0f8},2,Array{Gray{N0f8},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    for img in ( Gray.(A),
                 Gray.(N0f8.(A)),
                 Gray.(N0f16.(A)) )
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{Gray{N0f16},2,Array{Gray{N0f16}}},imgo))
        @test eltype(s) == Gray{N0f16}
        @test s isa OffsetArray{Gray{N0f16},2,Array{Gray{N0f16},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    # Color images wrapped by an OffsetArray.
    A = rand(RGB{Float32},8,8)
    for img in ( A,
                 n0f8.(A),
                 n6f10.(A),
                 n4f12.(A),
                 n2f14.(A),
                 n0f16.(A))
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{RGB{N0f8},2,Array{RGB{N0f8}}},imgo))
        @test eltype(s) == RGB{N0f8}
        @test s isa OffsetArray{RGB{N0f8},2,Array{RGB{N0f8},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end

    A = rand(RGB{Float32},8,8)
    for img in ( A,
                 n0f8.(A),
                 n6f10.(A),
                 n4f12.(A),
                 n2f14.(A),
                 n0f16.(A))
        imgo = OffsetArray(img, -2, -1)
        s = @inferred(convert(OffsetArray{RGB{Float32},2,Array{RGB{Float32}}},imgo))
        @test eltype(s) == RGB{Float32}
        @test s isa OffsetArray{RGB{Float32},2,Array{RGB{Float32},2}}
        @test permutedims(permutedims(s,(2,1)),(2,1)) == s
        @test axes(s) === axes(imgo)
    end
end


text:
#=
Given the root to a binary search tree, find the second largest node in the tree.
=#

using Test
include("Solutions/problem36_second_largest_node.jl")

test_node = Node(9, Node(8, nothing, Node(3, nothing, nothing)), Node(2, nothing, nothing))
#=
                9
             /     \
            8      2
             |
             3
=#
@test secondlargestnode(test_node) == Node(8, nothing, Node(3, nothing, nothing))

max_stars_repo_name: DominiqueCaron/daily-coding-problem

max_issues_repo_name: DominiqueCaron/daily-coding-problem

max_forks_repo_name: DominiqueCaron/daily-coding-problem


text:
export Instruction, unsafe_delete!,
       metadata, metadata!,
       opcode,
       predicate_int, predicate_real

# forward definition of Instruction in src/core/value/constant.jl
register(Instruction, API.LLVMInstructionValueKind)

const instruction_opcodes = Vector{Type}(fill(Nothing, typemax(API.LLVMOpcode)+1))
function identify(::Type{Instruction}, ref::API.LLVMValueRef)
    opcode = API.LLVMGetInstructionOpcode(ref)
    typ = @inbounds instruction_opcodes[opcode+1]
    typ === Nothing && error("Unknown type opcode $opcode")
    return typ
end
function register(T::Type{<:Instruction}, opcode::API.LLVMOpcode)
    instruction_opcodes[opcode+1] = T
end

function refcheck(::Type{T}, ref::API.LLVMValueRef) where T<:Instruction
    ref==C_NULL && throw(UndefRefError())
    if Base.JLOptions().debug_level >= 2
        T′ = identify(Instruction, ref)
        if T != T′
            error("invalid conversion of $T′ instruction reference to $T")
        end
    end
end

# Construct a concretely typed instruction object from an abstract value ref
function Instruction(ref::API.LLVMValueRef)
    ref == C_NULL && throw(UndefRefError())
    T = identify(Instruction, ref)
    return T(ref)
end

Instruction(inst::Instruction) =
    Instruction(API.LLVMInstructionClone(inst))

unsafe_delete!(::BasicBlock, inst::Instruction) =
    API.LLVMInstructionEraseFromParent(inst)
Base.delete!(::BasicBlock, inst::Instruction) =
    API.LLVMInstructionRemoveFromParent(inst)

parent(inst::Instruction) =
    BasicBlock(API.LLVMGetInstructionParent(inst))

opcode(inst::Instruction) = API.LLVMGetInstructionOpcode(inst)

predicate_int(inst::Instruction) = API.LLVMGetICmpPredicate(inst)
predicate_real(inst::Instruction) = API.LLVMGetFCmpPredicate(inst)

## metadata iteration
# TODO: doesn't actually iterate, since we can't list the available keys

@enum(MD, MD_dbg = 0,
          MD_tbaa = 1,
          MD_prof = 2,
          MD_fpmath = 3,
          MD_range = 4,
          MD_tbaa_struct = 5,
          MD_invariant_load = 6,
          MD_alias_scope = 7,
          MD_noalias = 8,
          MD_nontemporal = 9,
          MD_mem_parallel_loop_access = 10,
          MD_nonnull = 11,
          MD_dereferenceable = 12,
          MD_dereferenceable_or_null = 13,
          MD_make_implicit = 14,
          MD_unpredictable = 15,
          MD_invariant_group = 16,
          MD_align = 17,
          MD_loop = 18,
          MD_type = 19,
          MD_section_prefix = 20,
          MD_absolute_symbol = 21,
          MD_associated = 22)

export InstructionMetadataDict

struct InstructionMetadataDict <: AbstractDict{MD,MetadataAsValue}
    inst::Instruction
end

metadata(inst::Instruction) = InstructionMetadataDict(inst)

Base.isempty(md::InstructionMetadataDict) =
  !convert(Core.Bool, API.LLVMHasMetadata(md.inst))

Base.haskey(md::InstructionMetadataDict, kind::MD) =
  API.LLVMGetMetadata(md.inst, kind) != C_NULL

function Base.getindex(md::InstructionMetadataDict, kind::MD)
    objref = API.LLVMGetMetadata(md.inst, kind)
    objref == C_NULL && throw(KeyError(kind))
    return Metadata(MetadataAsValue(objref))
  end

Base.setindex!(md::InstructionMetadataDict, node::Metadata, kind::MD) =
    API.LLVMSetMetadata(md.inst, kind, Value(node; ctx=context(md.inst)))

Base.delete!(md::InstructionMetadataDict, kind::MD) =
    API.LLVMSetMetadata(md.inst, kind, C_NULL)

## instruction types

const opcodes = [:Ret, :Br, :Switch, :IndirectBr, :Invoke, :Unreachable, :CallBr, :FNeg,
                 :Add, :FAdd, :Sub, :FSub, :Mul, :FMul, :UDiv, :SDiv, :FDiv, :URem, :SRem,
                 :FRem, :Shl, :LShr, :AShr, :And, :Or, :Xor, :Alloca, :Load, :Store,
                 :GetElementPtr, :Trunc, :ZExt, :SExt, :FPToUI, :FPToSI, :UIToFP, :SIToFP,
                 :FPTrunc, :FPExt, :PtrToInt, :IntToPtr, :BitCast, :AddrSpaceCast, :ICmp,
                 :FCmp, :PHI, :Call, :Select, :UserOp1, :UserOp2, :VAArg, :ExtractElement,
                 :InsertElement, :ShuffleVector, :ExtractValue, :InsertValue, :Freeze,
                 :Fence, :AtomicCmpXchg, :AtomicRMW, :Resume, :LandingPad, :CleanupRet,
                 :CatchRet, :CatchPad, :CleanupPad, :CatchSwitch]

for op in opcodes
    typename = Symbol(op, :Inst)
    enum = Symbol(:LLVM, op)
    @eval begin
        @checked struct $typename <: Instruction
            ref::API.LLVMValueRef
        end
        register($typename, API.$enum)
    end
end

## call sites and invocations

export callconv, callconv!,
       istailcall, tailcall!,
       called_value

callconv(inst::Instruction) = API.LLVMGetInstructionCallConv(inst)
callconv!(inst::Instruction, cc) =
    API.LLVMSetInstructionCallConv(inst, cc)

istailcall(inst::Instruction) = convert(Core.Bool, API.LLVMIsTailCall(inst))
tailcall!(inst::Instruction, bool) = API.LLVMSetTailCall(inst, convert(Bool, bool))

called_value(inst::Instruction) = Value(API.LLVMGetCalledValue(inst))

## terminators

export isterminator, isconditional, condition, condition!, default_dest

isterminator(inst::Instruction) = API.LLVMIsATerminatorInst(inst) != C_NULL

isconditional(br::Instruction) = convert(Core.Bool, API.LLVMIsConditional(br))

condition(br::Instruction) =
    Value(API.LLVMGetCondition(br))
condition!(br::Instruction, cond::Value) =
    API.LLVMSetCondition(br, cond)

default_dest(switch::Instruction) =
    BasicBlock(API.LLVMGetSwitchDefaultDest(switch))

# successor iteration

export successors

struct TerminatorSuccessorSet <: AbstractVector{BasicBlock}
    term::Instruction
end

successors(term::Instruction) = TerminatorSuccessorSet(term)

Base.size(iter::TerminatorSuccessorSet) = (API.LLVMGetNumSuccessors(iter.term),)

Base.IndexStyle(::TerminatorSuccessorSet) = IndexLinear()

function Base.getindex(iter::TerminatorSuccessorSet, i::Int)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return BasicBlock(API.LLVMGetSuccessor(iter.term, i-1))
end

Base.setindex!(iter::TerminatorSuccessorSet, bb::BasicBlock, i::Int) =
    API.LLVMSetSuccessor(iter.term, i-1, bb)

## phi nodes

# incoming iteration

export PhiIncomingSet

struct PhiIncomingSet <: AbstractVector{Tuple{Value,BasicBlock}}
    phi::Instruction
end

incoming(phi::Instruction) = PhiIncomingSet(phi)

Base.size(iter::PhiIncomingSet) = (API.LLVMCountIncoming(iter.phi),)

Base.IndexStyle(::PhiIncomingSet) = IndexLinear()

function Base.getindex(iter::PhiIncomingSet, i::Int)
    @boundscheck 1 <= i <= length(iter) || throw(BoundsError(iter, i))
    return tuple(Value(API.LLVMGetIncomingValue(iter.phi, i-1)),
                       BasicBlock(API.LLVMGetIncomingBlock(iter.phi, i-1)))
end

function Base.append!(iter::PhiIncomingSet, args::Vector{Tuple{V, BasicBlock}} where V <: Value)
    vals, blocks = zip(args...)
    API.LLVMAddIncoming(iter.phi, collect(vals), collect(blocks), length(args))
end

Base.push!(iter::PhiIncomingSet, args::Tuple{<:Value, BasicBlock}) = append!(iter, [args])


text:
using Random
include("struct/tree.jl")

"""
Création de deux listes d'indices pour les jeux de données d'entrainement et de test
Entrées :\n
    - n le nombre de données dans le dataset
    - p la proportion que représente le jeu de données test
Sorties :\n
    - train, la liste des indices des données d'entrainement
    - test, la liste des indices des données de tesr
"""
function train_test_indexes(n::Int64,p::Float64=0.2)

    # Fixe la graine aléatoire pour assurer la reproductibilité
    Random.seed!(1)
    rd = randperm(n)

    test = rd[1:ceil(Int,n*p)]
    train = rd[ceil(Int,n*p)+1:n]

    return train,test
end

"""
Retourne le nombre d'erreurs de prédiction d'un arbre pour un ensemble de données

Entrées :
- T : l'arbre
- x : les données à prédire
- y : la classe des données

Sortie :
- class::Vector{Int64} : class prédites (class[i] est la classe de la donnée x[i, :])
"""
function prediction_errors(T::Tree, x::Matrix{Float64}, y::Vector{Int64})
    dataCount = length(x[:, 1])
    featuresCount = length(x[1, :])

    errors = 0

    for i in 1:dataCount
        t = 1
        for d in 1:(T.D+1)
            if T.c[t] != -1
                errors += T.c[t] != y[i]
                break
            else
                if sum(T.a[j, t]*x[i, j] for j in 1:featuresCount) - T.b[t] < 0
                    t = t*2
                else
                    t = t*2 + 1
                end
            end
        end
    end
    return errors
end

"""
Retourne la prédiction d'un arbre pour un ensemble de données

Entrées :
- T : l'arbre
- x : les données à prédire

Sortie :
- class::Vector{Int64} : class prédites (class[i] est la classe de la donnée x[i, :])
"""
function predict_class(T::Tree, x::Matrix{Float64})
    dataCount = length(x[:, 1])
    featuresCount = length(x[1, :])
    class = zeros(Int64, dataCount)

    for i in 1:dataCount
        t = 1
        for d in 1:(T.D+1)
            if T.c[t] != -1
                class[i] = T.c[t]
                break
            else
                if sum(T.a[j, t]*x[i, j] for j in 1:featuresCount) - T.b[t] < 0
                    t = t*2
                else
                    t = t*2 + 1
                end
            end
        end
    end
    return class
end

"""
Change l'échelle des caractéristiques d'un dataset pour les situer dans [0, 1]

Entrée :
- X: les caractéristiques du dataset d'origine

Sortie :
- caractéristiques entre 0 et 1
"""
function centerData(X)

    result = Matrix{Float64}(X)

    # Pour chaque caractéristique
    for j in 1:size(result, 2)

        m = minimum(result[:, j])
        M = maximum(result[:, j])
        result[:, j] .-= m
        result[:, j] ./= M
    end

    return result
end

function centerAndSaveDataSet(X, Y::Vector{Int64}, outputFile::String)

    centeredX = centerData(X)

    open(outputFile, "w") do fout
        println(fout, "X = ", centeredX)
        println(fout, "Y = ", Y)
    end
end

max_stars_repo_name: boumpteryx/RODM

max_issues_repo_name: boumpteryx/RODM

max_forks_repo_name: boumpteryx/RODM


text:
"""
    metropolis(rbm, v; β = 1, steps = 1)

Metropolis-Hastings sampling from `rbm` at inverse temperature β, starting from
configuration `v`.
"""
function metropolis(rbm::RBM, v::AbstractArray; β::Real = 1, steps::Int = 1)
    @assert size(v)[1:ndims(rbm.visible)] == size(rbm.visible)
    for _ in 1:steps
        v = oftype(v, metropolis_once(rbm, v; β))
    end
    return v
end

function metropolis_once(rbm::RBM, v::AbstractArray; β::Real = 1)
    @assert size(v)[1:ndims(rbm.visible)] == size(rbm.visible)
    v_new = sample_v_from_v_once(rbm, v)
    ΔE = (β - 1) * (free_energy(rbm, v_new) - free_energy(rbm, v))
    ℐ = CartesianIndices(size(rbm.visible)) # index span of visible layer
    for n in CartesianIndices(ΔE)
        if ΔE[n] ≤ 0 || randexp() > ΔE[n]
            continue # accept move
        else
            v_new[ℐ, n] .= v[ℐ, n] # do not accept move
        end
    end
    return v_new
end

"""
    metropolis!(v, rbm; β = 1)

Metropolis-Hastings sampling from `rbm` at inverse temperature β.
Uses `v[:,:,..,:,1]` as initial configurations, and writes the
Monte-Carlo chains in `v[:,:,..,:,2:end]`.
"""
function metropolis!(v::AbstractArray, rbm::RBM; β::Real = 1)
    @assert ndims(v) > ndims(rbm.visible) # last dim contains nsteps
    @assert size(v)[1:ndims(rbm.visible)] == size(rbm.visible)
    for t in 2:size(v, ndims(v))
        v[.., t] .= metropolis_once(rbm, v[.., t - 1]; β)
    end
    return v
end

"""
    cold_metropolis(rbm, v; steps = 1)

Samples the `rbm` at zero temperature, starting from configuration `v`.
"""
function cold_metropolis(rbm::RBM, v::AbstractArray; steps::Int = 1)
    @assert size(v)[1:ndims(rbm.visible)] == size(rbm.visible)
    for _ in 1:steps
        v = oftype(v, cold_metropolis_once(rbm, v))
    end
    return v
end

function cold_metropolis_once(rbm::RBM, v::AbstractArray)
    @assert size(v)[1:ndims(rbm.visible)] == size(rbm.visible)
    h = mean_h_from_v(rbm, v)
    return mode_v_from_h(rbm, h)
end


text:
using Images, ImageCore
args = ARGS
if length(args) == 0
  println("No filename given.")
  exit(-1)
end
filen = shift!(args)

topclip = 0.01 # allow this proportion of pixels to be clipped into the top bin
bottomclip = 0.01 # allow this proportion of pixels to be clipped into the bottom bin

print("Loading $filen... ")
img = load(filen)
if img == nothing
  println("Image failed to load.")
  exit(-1)
end
println("Done!")
println("Loading processing vectors and histograms...")
sz = size(img)
px = sz[1]
py = sz[2]
pixcount = px * py
println("\Image size: ", px, "x", py, " pixels")
img = reinterpret(UInt16, img)
imgr = reshape(img[1, :, :], pixcount)
imgg = reshape(img[2, :, :], pixcount)
imgb = reshape(img[3, :, :], pixcount)
imgmax = zeros(imgr)
pixmax = zeros(UInt32, 65536)
for i = 1:pixcount
  peak = imgr[i]
  peak = imgg[i] > peak ? imgg[i] : peak
  peak = imgb[i] > peak ? imgb[i] : peak
  imgmax[i] = peak
  pixmax[peak + 1] += 1
end
println("Image and histograms loaded")

mapto = zeros(UInt32, 1:65536)
topcount = floor(pixcount * topclip)
sofar = 0
ptr = 65536
while sofar < topcount
  sofar += pixmax[ptr]
  mapto[ptr] = 65536
  ptr -= 1
end

topptr = ptr
println("Top clip level = ", topptr)

bottomcount = floor(pixcount * bottomclip)
sofar = 0
ptr = 1
while sofar < bottomcount
  sofar += pixmax[ptr]
  mapto[ptr] = 0
  ptr += 1
end

bottomptr = ptr
println("Bottom clip level = ", bottomptr)

# Linear moving
for i = bottomptr:topptr
  mapto[i] = floor(65536 * (i - bottomptr) / (topptr - bottomptr))
end

ten = floor(pixcount / 10)
ctr = ten
for i = 1:pixcount
  factor = mapto[imgmax[i] + 1] / (imgmax[i] + 1)
  rf = imgr[i] < bottomptr ? imgr[i] : bottomptr
  gf = imgg[i] < bottomptr ? imgg[i] : bottomptr
  bf = imgb[i] < bottomptr ? imgb[i] : bottomptr
  ctr -= 1
  if ctr == 0
    println(i, " of ", pixcount)
    ctr = ten
  end
  imgr[i] = floor(UInt16, factor * (imgr[i] - rf) + rf)
  imgg[i] = floor(UInt16, factor * (imgg[i] - gf) + gf)
  imgb[i] = floor(UInt16, factor * (imgb[i] - bf) + bf)
end

final = zeros(UInt16, 3, px, py)
final[1, :, :] = imgr
final[2, :, :] = imgg
final[3, :, :] = imgb

final = colorview(RGB{N0f16}, final)
#final = reinterpret()
oname = "q$filen"
print("Saving $oname... ")
save(oname, final)
println("Done!")

max_stars_repo_name: FatherGeoffHorton/AstroTools

max_issues_repo_name: FatherGeoffHorton/AstroTools

max_forks_repo_name: FatherGeoffHorton/AstroTools


text:
using Distributions: Multivariate, Distribution
# # Multivariate
export Mv

"""
Multivariate distribution: Random array where each variable is ciid given
values for parameters.

`Mv(dist, shape)`

# Arguments
- `dist` a variable class, i.e. `dist(id, ω)` must b defined, e.g. `Normal(0, 1)`
- `shape` Dimensions of Multivariate

# Returns

# Example
```julia
x = 1 ~ Normal(0, 1)
function f(id, ω)
  x(ω) + Uniform(0, 1)(id, ω)
end
xs = 2 ~ Mv(f, (3, 3))
randsample((x, xs))
```
"""
struct Mv{T, SHAPE}
  dist::T
  shape::SHAPE
end

Mv(dist, N::Integer) = Mv(dist, (N,))

traitlift(::Type{<:Mv}) = Lift()

# Base.eltype(Mv{T}) where {T} =

prim(d::Normal) = StdNormal()

func(d::Normal, x) = x * d.σ + d.μ

@inline Var.recurse(mv::Mv{<:Distribution}, id, ω) =
  map(x -> func(mv.dist, x), resolve(Mv(prim(mv.dist), mv.shape), id, ω))

f(x::Dims) = map(i->1:i, x)
g(x::Dims) = Iterators.product(f(x)...)
@inline Var.recurse(mv::Mv{<:T}, id, ω) where T =
  map(id_ -> mv.dist(append(id, id_), ω), g(mv.shape))

Base.rand(rng::AbstractRNG, mv::Mv{<:PrimDist}) =
  rand(rng, mv.dist, mv.shape)


text:
#=

    Indian GPA model

    From cplint http://cplint.ml.unife.it/example/inference/indian_gpa.pl
    """
    The Indian GPA Problem. From
    http://www.robots.ox.ac.uk/~fwood/anglican/examples/viewer/?worksheet=indian-gpa
    "This example was inspired by Stuart Russell...the problem is: if you observe
    that a student GPA is exactly 4.0 in a model of transcripts of students
    from the USA (GPA's from 0.0 to 4.0 and India (GPA's from 0.0 to
    10.0) what is the probability that the student is from India?...
    As we know from statistics, given the mixture distribution and given the
    fact that his/her GPA is exactly 4.0, the probability that the student
    is American must be 1.0.
    (i.e. zero probability that the student is from India)."
    Probabilistic logic program from
    https://github.com/davidenitti/DC/blob/master/examples/indian-gpa.pl

    """

    Also see the SPPL model indian-gpa.py which I used as inspiration for
    this Turing.jl model.

    Some tests:
    * Plain model
        Summary Statistics
        parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

        nationality    1.5230    0.4997     0.0158    0.0194   604.5633    0.9997      458.6975
            perfect    0.1520    0.3592     0.0114    0.0143   690.4058    1.0015      523.8283
                gpa    4.0086    2.9342     0.0928    0.1346   502.5789    0.9991      381.3194

        Quantiles
        parameters      2.5%     25.0%     50.0%     75.0%     97.5%
            Symbol   Float64   Float64   Float64   Float64   Float64

        nationality    1.0000    1.0000    2.0000    2.0000    2.0000
            perfect    0.0000    0.0000    0.0000    0.0000    1.0000
                gpa    0.2925    1.6643    3.3376    5.8096   10.0000

        Distributions of variable nationality
        usa        =>     523  (0.523000)
        india      =>     477  (0.477000)

    * true ~ Dirac(gpa == 4)
        Summary Statistics
        parameters      mean       std   naive_se      mcse         ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64     Float64   Float64       Float64

        nationality    1.9990    0.0316     0.0010    0.0010   1000.0080    1.0000      646.8357
            perfect    0.9980    0.0447     0.0014    0.0020    400.4815    1.0010      259.0437
                gpa    3.9959    0.1020     0.0032    0.0041    543.5300    1.0006      351.5718

        Distributions of variable nationality
        usa        =>     999  (0.999000)
        india      =>       1  (0.001000)

    * true ~ Dirac(perfect)
        Summary Statistics
        parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

        nationality    1.5200    0.4998     0.0158    0.0475    82.7796    0.9990       56.3893
            perfect    0.9980    0.0447     0.0014    0.0020   400.4815    1.0010      272.8076
                gpa    6.8704    3.0038     0.0950    0.2832    83.1920    0.9990       56.6703

        Distributions of variable nationality
        usa        =>     520  (0.520000)
        india      =>     480  (0.480000)

    * true ~ Dirac( (nationality == usa && gpa > 3) || (gpa >= 8 && gpa < 10))
        Summary Statistics
        parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

        nationality    1.7100    0.4540     0.0144    0.0326   226.4175    1.0000      152.7784
            perfect    0.2510    0.4338     0.0137    0.0303   234.7464    1.0005      158.3984
                gpa    5.2539    2.4201     0.0765    0.1821   223.4711    1.0000      150.7902

        Quantiles
        parameters      2.5%     25.0%     50.0%     75.0%     97.5%
            Symbol   Float64   Float64   Float64   Float64   Float64

        nationality    1.0000    1.0000    2.0000    2.0000    2.0000
            perfect    0.0000    0.0000    0.0000    1.0000    1.0000
                gpa    3.1515    3.6150    4.0000    8.1999    9.7962

        Distributions of variable nationality
        usa        =>     710  (0.710000)
        india      =>     290  (0.290000)

    * t1 ~ Dirac( (nationality == usa && gpa > 3) || (gpa >=8 && gpa < 10))
        Summary Statistics
        parameters      mean       std   naive_se      mcse        ess      rhat   ess_per_sec
            Symbol   Float64   Float64    Float64   Float64    Float64   Float64       Float64

        nationality    1.5020    0.5002     0.0158    0.0246   600.6463    0.9996      375.8738
            perfect    0.1120    0.3155     0.0100    0.0068   657.5644    0.9990      411.4921
                gpa    3.7279    2.8090     0.0888    0.1057   670.9449    0.9990      419.8654
                t1    0.2720    0.4452     0.0141    0.0154   645.9312    0.9990      404.2123

        Quantiles
        parameters      2.5%     25.0%     50.0%     75.0%     97.5%
            Symbol   Float64   Float64   Float64   Float64   Float64

        nationality    1.0000    1.0000    2.0000    2.0000    2.0000
            perfect    0.0000    0.0000    0.0000    0.0000    1.0000
                gpa    0.1251    1.4697    3.1822    4.9678   10.0000
                t1    0.0000    0.0000    0.0000    1.0000    1.0000

        Distributions of variable nationality
        usa        =>     502  (0.502000)
        india      =>     498  (0.498000)

=#

using Turing, StatsPlots, Distributions, StatsBase
using DataFrames,RDatasets
include("jl_utils.jl")

@model function indian_gpa()
    india = 1
    usa   = 2
    nationality ~ Categorical([0.5,0.5])
    if nationality == india
        perfect ~ flip(0.10)
        if perfect
            gpa ~ Dirac(10.0)
        else
            gpa ~ Uniform(0,10)
        end
    else
        perfect ~ flip(0.15)
        if perfect
            gpa ~ Dirac(4.0)
        else
            gpa ~ Uniform(0,4)
        end

    end

    # true ~ Dirac(gpa == 4)
    # true ~ Dirac(gpa > 4)
    # true ~ Dirac(perfect)
    # true ~ Dirac( (nationality == usa && gpa > 3) || (gpa >= 8 && gpa < 10))
    # t1 ~ Dirac( (nationality == usa && gpa > 3) || (gpa >=8 && gpa < 10))

end

model = indian_gpa()

# chns = sample(model, Prior(), 10_000)
# chns = sample(model, MH(), 1_000)
chns = sample(model, PG(5), 1_000)
# chns = sample(model, SMC(), 1_000)
# chns = sample(model, IS(), 10_000)

# chns = sample(model, HMC(0.1,6), 1_000)
# chns = sample(model, NUTS(), 1_000)

display(chns)

show_var_dist_pct(chns,:nationality,["india","usa"])

max_stars_repo_name: tias/hakank

max_issues_repo_name: tias/hakank

max_forks_repo_name: tias/hakank


text:
# ------------------------------------------------------------------
# Licensed under the ISC License. See LICENSE in the project root.
# ------------------------------------------------------------------

"""
    EmpiricalDistribution(values)

An empirical distribution holding continuous values.
"""
struct EmpiricalDistribution{T<:Real} <: ContinuousUnivariateDistribution
  values::Vector{T}

  function EmpiricalDistribution{T}(values) where {T<:Real}
    @assert !isempty(values) "values must be provided"
    new(sort(values))
  end
end
EmpiricalDistribution(values) = EmpiricalDistribution{eltype(values)}(values)

quantile(d::EmpiricalDistribution, p::Real) = quantile(d.values, p, sorted=true)

function cdf(d::EmpiricalDistribution{T}, x::T) where {T<:Real}
  v = d.values
  N = length(v)

  head, mid, tail = 1, 1, N
  while tail - head > 1
    mid = (head + tail) ÷ 2
    if x < v[mid]
      tail = mid
    else
      head = mid
    end
  end

  l, u = v[head], v[tail]

  if x < l
    return 0.
  elseif x > u
    return 1.
  else
    if head == tail
      return head / N
    else
      pl, pu = head / N, tail / N
      return (pu - pl) * (x - l) / (u - l) + pl
    end
  end
end

"""
    transform!(samples, origin, target)

Transform `samples` from `origin` distribution to `target` distribution in place.
"""
function transform!(samples::AbstractVector,
                    origin::ContinuousUnivariateDistribution,
                    target::ContinuousUnivariateDistribution)
  for i in eachindex(samples)
    p = cdf(origin, samples[i]) - eps()
    samples[i] = quantile(target, p)
  end
end


text:
struct QuadCost{TQ<:SMatrix}
    Q::TQ
end
@inline (c::QuadCost)(x) = 1/2*x'*c.Q*x
@inline function quad!(Q::MMatrix, l::MVector, c::QuadCost,
                       x::SVector{n}, xi::SVector{n}) where {n}
    Q[xi, xi] += c.Q
    l[xi] += c.Q * x
    return nothing
end

# Interface to implement NPlayerNavigationCost
@with_kw struct SoftConstr
    "The index of the state/input element this applies to."
    id::Int
    "The soft lower bound."
    min::Float64
    "The soft upper bound."
    max::Float64
    "The weight of the soft constraint."
    w::Float64
end
@inline function softconstr(val::Real, min::Real, max::Real, w::Real)
    @assert min < max
    gap = val < min ? val - min : val > max ? val - max : zero(w)
    return w*gap^2
end
@inline function (constr::SoftConstr)(x, xi)
    @unpack id, min, max, w = constr
    val = x[xi[id]]
    return softconstr(val, min, max, w)
end

@inline function quad!(Q::MMatrix, l::MVector, constr::SoftConstr, x::SVector, xi)
    @unpack id, w, min, max = constr
    idx = xi[id]

    @assert min < max
    if x[idx] < min
        Q[idx, idx] += 2w
        l[idx] += 2w*(x[idx] - min)
    elseif x[idx] > max
        Q[idx, idx] += 2w
        l[idx] += 2w*(x[idx] - max)
    end
    return nothing
end

@with_kw struct ProximityCost{TR,TW}
    "A collection of avoidance raddii this player want's to keep from others."
    rs::TR
    "A collection of costs this player has for colliding with others."
    ws::TW
end
ProximityCost(r, w, np) = ProximityCost(SVector{np}([r for _ in 1:np]),
                                        SVector{np}([w for _ in 1:np]))

"Proxmity cost for player i to player j. j is the index of the *other* player."
@inline function (pc::ProximityCost)(xp_i::SVector, xp_j::SVector, j::Int)
    @unpack rs, ws = pc
    Δxp = xp_i - xp_j
    return softconstr(sqrt(Δxp'*Δxp), rs[j], Inf, ws[j])
end

@inline function quad!(Q::MMatrix, l::MVector, pc::ProximityCost, x::SVector,
                       xyi_1, xyi_2, j::Int)
    @unpack rs, ws = pc
    r_avoid = rs[j]
    w = ws[j]

    x1, y1 = xyi_1[1], xyi_1[2]
    x2, y2 = xyi_2[1], xyi_2[2]
    Δx = x[x1] - x[x2]
    Δy = x[y1] - x[y2]
    Δnormsq = Δx^2 + Δy^2
    Δnorm = sqrt(Δnormsq)
    if Δnorm < r_avoid
        # cost model: w*(Δnorm - min)^2
        δx = 2*w*Δx*(Δnorm - r_avoid)/Δnorm
        δy = 2*w*Δy*(Δnorm - r_avoid)/Δnorm
        l[x1] += δx
        l[y1] += δy
        l[x2] -= δx
        l[y2] -= δy

        w̃ = 2*w/Δnormsq

        δxδx = (w̃*Δx^2) + (w̃*Δx^2*(r_avoid - Δnorm))/Δnorm - (w̃*Δnorm*(r_avoid - Δnorm))
        δyδy = (w̃*Δy^2) + (w̃*Δy^2*(r_avoid - Δnorm))/Δnorm - (w̃*Δnorm*(r_avoid - Δnorm))
        δxδy = (w̃*Δx*Δy) + (w̃*Δx*Δy*(r_avoid - Δnorm))/Δnorm

        Q[x1,x1] += δxδx
        Q[x1,x2] -= δxδx; Q[x2,x1] -= δxδx
        Q[x2,x2] += δxδx

        Q[y1,y1] += δyδy
        Q[y1,y2] -= δyδy; Q[y2,y1] -= δyδy
        Q[y2,y2] += δyδy

        Q[x1,y1] += δxδy; Q[y1,x1] += δxδy
        Q[x1,y2] -= δxδy; Q[y2,x1] -= δxδy
        Q[y1,x2] -= δxδy; Q[x2,y1] -= δxδy
        Q[x2,y2] += δxδy; Q[y2,x2] += δxδy
    end
end

@with_kw struct GoalCost{TG<:SVector, TQg<:SMatrix}
    t_active::Float64
    xg::TG
    Qg::TQg
end

@inline function (gc::GoalCost)(x, t)
    @unpack t_active, xg, Qg = gc
    if t >= t_active
        Δx = x - xg
        return 1/2 * Δx'*Qg*Δx
    else
        return 0.0
    end
end

@inline function quad!(Q::MMatrix, l::MVector, gc::GoalCost, x::SVector,
                       xid::SVector{n}, t::AbstractFloat) where {n}
    @unpack t_active, xg, Qg = gc
    if t >= t_active
        Q[xid, xid] += Qg
        l[xid] += Qg*(x-xg)
    end
    return nothing
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.
# *********************************************************************************
"""
    Degradation

Inputs used when `ElectricStorage.model_degradation` is `true`:
```julia
Base.@kwdef mutable struct Degradation
    calendar_fade_coefficient::Float64 = 2.46E-03
    cycle_fade_coefficient::Float64 = 7.82E-05
    installed_cost_per_kwh_declination_rate::Float64 = 0.05
    maintenance_strategy::String = "augmentation"  # one of ["augmentation", "replacement"]
    maintenance_cost_per_kwh::Vector{<:Real} = Real[]
end
```

None of the above values are required. If `ElectricStorage.model_degradation` is `true` then the
defaults above are used.
If the `maintenance_cost_per_kwh` is not provided then it is determined using the `ElectricStorage.installed_cost_per_kwh`
and the `installed_cost_per_kwh_declination_rate` along with a present worth factor ``f`` to account for the present cost
of buying a battery in the future. The present worth factor for each day is:

``
f(day) = \\frac{ (1-r_g)^\\frac{day}{365} } { (1+r_d)^\\frac{day}{365} }
``

where ``r_g`` = `installed_cost_per_kwh_declination_rate` and ``r_d`` = `p.s.financial.owner_discount_pct`.

Note this day-specific calculation of the present-worth factor accumulates differently from the annually updated discount
rate for other net-present value calculations in REopt, and has a higher effective discount rate as a result.  The present
worth factor is used in two different ways, depending on the `maintenance_strategy`, which is described below.

!!! warn
    When modeling degradation the following ElectricStorage inputs are not used:
    - `replace_cost_per_kw`
    - `replace_cost_per_kwh`
    - `inverter_replacement_year`
    - `battery_replacement_year`
    The are replaced by the `maintenance_cost_per_kwh` vector.

!!! note
    When providing the `maintenance_cost_per_kwh` it must have a length equal to `Financial.analysis_years*365`.

# Battery State Of Health
The state of health [`SOH`] is a linear function of the daily average state of charge [`Eavg`] and
the daily equivalent full cycles [`EFC`]. The initial `SOH` is set to the optimal battery energy capacity
(in kWh). The evolution of the `SOH` beyond the first day is:

``
SOH[d] = SOH[d-1] - h\\left(
    \\frac{1}{2} k_{cal} Eavg[d-1] / \\sqrt{d} + k_{cyc} EFC[d-1] \\quad \\forall d \\in \\{2\\dots D\\}
\\right)
``

where:
- ``k_{cal}`` is the `calendar_fade_coefficient`
- ``k_{cyc}`` is the `cycle_fade_coefficient`
- ``h`` is the hours per time step
- ``D`` is the total number of days, 365 * `analysis_years`

The `SOH` is used to determine the maintence cost of the storage system, which depends on the `maintenance_strategy`.

# Augmentation Maintenance Strategy
The augmentation maintenance strategy assumes that the battery energy capacity is maintained by replacing
degraded cells daily in terms of cost. Using the definition of the `SOH` above the maintenance cost is:

``
C_{\\text{aug}} = \\sum_{d \\in \\{2\\dots D\\}} 0.8 C_{\\text{install}} f(day) \\left( SOH[d-1] - SOH[d] \\right)
``

where
- the ``0.8`` factor accounts for sunk costs that do not need to be paid;
- ``C_{\\text{install}}`` is the `ElectricStorage.installed_cost_per_kwh`; and
- ``SOH[d-1] - SOH[d]`` is the incremental amount of battery capacity lost in a day.

The ``C_{\\text{aug}}`` is added to the objective function to be minimized with all other costs.

# Replacement Maintenance Strategy
Modeling the replacement maintenance strategy is more complex than the augmentation strategy.
Effectively the replacement strategy says that the battery has to be replaced once the `SOH` hits 80%
of the optimal, purchased capacity. It is possible that multiple replacements could be required under
this strategy.

!!! warn
    The "replacement" maintenance strategy requires integer variables and indicator constraints.
    Not all solvers support indicator constraints and some solvers are slow with integer variables.

The replacement strategy cost is:

``
C_{\\text{repl}} = B_{\\text{kWh}} N_{\\text{repl}} f(d_{80}) C_{\\text{install}}
``

where:
- ``B_{\\text{kWh}}`` is the optimal battery capacity (`ElectricStorage.size_kwh` in the results dictionary);
- ``N_{\\text{repl}}`` is the number of battery replacments required (a function of the month in which the `SOH` reaches 80% of original capacity);
-  ``f(d_{80})`` is the present worth factor at approximately the 15th day of the month that the `SOH` reaches 80% of original capacity.

The ``C_{\\text{repl}}`` is added to the objective function to be minimized with all other costs.

# Example of inputs
The following shows how one would use the degradation model in REopt via the [Scenario](@ref) inputs:
```javascript
{
    ...
    "ElectricStorage": {
        "installed_cost_per_kwh": 390,
        ...
        "model_degradation": true,
        "degradation": {
            "calendar_fade_coefficient": 2.86E-03,
            "cycle_fade_coefficient": 6.22E-05,
            "installed_cost_per_kwh_declination_rate": 0.06,
            "maintenance_strategy": "replacement",
            ...
        }
    },
    ...
}
```
Note that not all of the above inputs are necessary. When not providing `calendar_fade_coefficient` for example the default value will be used.
"""
Base.@kwdef mutable struct Degradation
    calendar_fade_coefficient::Float64 = 2.46E-03
    cycle_fade_coefficient::Float64 = 7.82E-05
    installed_cost_per_kwh_declination_rate::Float64 = 0.05
    maintenance_strategy::String = "augmentation"  # one of ["augmentation", "replacement"]
    maintenance_cost_per_kwh::Vector{<:Real} = Real[]
end

"""
    ElectricStorageDefaults

Electric storage system defaults. Overridden by user inputs provided in `ElectricStorage`.

```julia
Base.@kwdef struct ElectricStorageDefaults
    off_grid_flag::Bool = false
    min_kw::Real = 0.0
    max_kw::Real = 1.0e4
    min_kwh::Real = 0.0
    max_kwh::Real = 1.0e6
    internal_efficiency_pct::Float64 = 0.975
    inverter_efficiency_pct::Float64 = 0.96
    rectifier_efficiency_pct::Float64 = 0.96
    soc_min_pct::Float64 = 0.2
    soc_init_pct::Float64 = off_grid_flag ? 1.0 : 0.5
    can_grid_charge::Bool = off_grid_flag ? false : true
    installed_cost_per_kw::Real = 775.0
    installed_cost_per_kwh::Real = 388.0
    replace_cost_per_kw::Real = 440.0
    replace_cost_per_kwh::Real = 220.0
    inverter_replacement_year::Int = 10
    battery_replacement_year::Int = 10
    macrs_option_years::Int = 7
    macrs_bonus_pct::Float64 = 1.0
    macrs_itc_reduction::Float64 = 0.5
    total_itc_pct::Float64 = 0.0
    total_rebate_per_kw::Real = 0.0
    total_rebate_per_kwh::Real = 0.0
    charge_efficiency::Float64 = rectifier_efficiency_pct * internal_efficiency_pct^0.5
    discharge_efficiency::Float64 = inverter_efficiency_pct * internal_efficiency_pct^0.5
    grid_charge_efficiency::Float64 = can_grid_charge ? charge_efficiency : 0.0
    model_degradation::Bool = false
    degradation::Dict = Dict()
    minimum_avg_soc_fraction::Float64 = 0.0
end
```
"""
Base.@kwdef struct ElectricStorageDefaults
    off_grid_flag::Bool = false
    min_kw::Real = 0.0
    max_kw::Real = 1.0e4
    min_kwh::Real = 0.0
    max_kwh::Real = 1.0e6
    internal_efficiency_pct::Float64 = 0.975
    inverter_efficiency_pct::Float64 = 0.96
    rectifier_efficiency_pct::Float64 = 0.96
    soc_min_pct::Float64 = 0.2
    soc_init_pct::Float64 = off_grid_flag ? 1.0 : 0.5
    can_grid_charge::Bool = off_grid_flag ? false : true
    installed_cost_per_kw::Real = 775.0
    installed_cost_per_kwh::Real = 388.0
    replace_cost_per_kw::Real = 440.0
    replace_cost_per_kwh::Real = 220.0
    inverter_replacement_year::Int = 10
    battery_replacement_year::Int = 10
    macrs_option_years::Int = 7
    macrs_bonus_pct::Float64 = 1.0
    macrs_itc_reduction::Float64 = 0.5
    total_itc_pct::Float64 = 0.0
    total_rebate_per_kw::Real = 0.0
    total_rebate_per_kwh::Real = 0.0
    charge_efficiency::Float64 = rectifier_efficiency_pct * internal_efficiency_pct^0.5
    discharge_efficiency::Float64 = inverter_efficiency_pct * internal_efficiency_pct^0.5
    grid_charge_efficiency::Float64 = can_grid_charge ? charge_efficiency : 0.0
    model_degradation::Bool = false
    degradation::Dict = Dict()
    minimum_avg_soc_fraction::Float64 = 0.0
end

"""
    function ElectricStorage(d::Dict, f::Financial, settings::Settings)

Construct ElectricStorage struct from Dict with keys-val pairs from the
REopt ElectricStorage and Financial inputs.
"""
struct ElectricStorage <: AbstractElectricStorage
    min_kw::Real
    max_kw::Real
    min_kwh::Real
    max_kwh::Real
    internal_efficiency_pct::Float64
    inverter_efficiency_pct::Float64
    rectifier_efficiency_pct::Float64
    soc_min_pct::Float64
    soc_init_pct::Float64
    can_grid_charge::Bool
    installed_cost_per_kw::Real
    installed_cost_per_kwh::Real
    replace_cost_per_kw::Real
    replace_cost_per_kwh::Real
    inverter_replacement_year::Int
    battery_replacement_year::Int
    macrs_option_years::Int
    macrs_bonus_pct::Float64
    macrs_itc_reduction::Float64
    total_itc_pct::Float64
    total_rebate_per_kw::Real
    total_rebate_per_kwh::Real
    charge_efficiency::Float64
    discharge_efficiency::Float64
    grid_charge_efficiency::Float64
    net_present_cost_per_kw::Real
    net_present_cost_per_kwh::Real
    model_degradation::Bool
    degradation::Degradation
    minimum_avg_soc_fraction::Float64

    function ElectricStorage(d::Dict, f::Financial)
        s = ElectricStorageDefaults(;d...)

        if s.inverter_replacement_year >= f.analysis_years
            @warn "Battery inverter replacement costs (per_kw) will not be considered because inverter_replacement_year >= analysis_years."
        end

        if s.battery_replacement_year >= f.analysis_years
            @warn "Battery replacement costs (per_kwh) will not be considered because battery_replacement_year >= analysis_years."
        end

        net_present_cost_per_kw = effective_cost(;
            itc_basis = s.installed_cost_per_kw,
            replacement_cost = s.inverter_replacement_year >= f.analysis_years ? 0.0 : s.replace_cost_per_kw,
            replacement_year = s.inverter_replacement_year,
            discount_rate = f.owner_discount_pct,
            tax_rate = f.owner_tax_pct,
            itc = s.total_itc_pct,
            macrs_schedule = s.macrs_option_years == 7 ? f.macrs_seven_year : f.macrs_five_year,
            macrs_bonus_pct = s.macrs_bonus_pct,
            macrs_itc_reduction = s.macrs_itc_reduction,
            rebate_per_kw = s.total_rebate_per_kw
        )
        net_present_cost_per_kwh = effective_cost(;
            itc_basis = s.installed_cost_per_kwh,
            replacement_cost = s.battery_replacement_year >= f.analysis_years ? 0.0 : s.replace_cost_per_kwh,
            replacement_year = s.battery_replacement_year,
            discount_rate = f.owner_discount_pct,
            tax_rate = f.owner_tax_pct,
            itc = s.total_itc_pct,
            macrs_schedule = s.macrs_option_years == 7 ? f.macrs_seven_year : f.macrs_five_year,
            macrs_bonus_pct = s.macrs_bonus_pct,
            macrs_itc_reduction = s.macrs_itc_reduction
        )

        net_present_cost_per_kwh -= s.total_rebate_per_kwh

        if haskey(d, :degradation)
            degr = Degradation(;dictkeys_tosymbols(d[:degradation])...)
        else
            degr = Degradation()
        end

        # copy the replace_costs in case we need to change them
        replace_cost_per_kw = s.replace_cost_per_kw
        replace_cost_per_kwh = s.replace_cost_per_kwh
        if s.model_degradation
            if haskey(d, :replace_cost_per_kw) && d[:replace_cost_per_kw] != 0.0 ||
                haskey(d, :replace_cost_per_kwh) && d[:replace_cost_per_kwh] != 0.0
                @warn "Setting ElectricStorage replacement costs to zero. \nUsing degradation.maintenance_cost_per_kwh instead."
            end
            replace_cost_per_kw = 0.0
            replace_cost_per_kwh = 0.0
        end

        return new(
            s.min_kw,
            s.max_kw,
            s.min_kwh,
            s.max_kwh,
            s.internal_efficiency_pct,
            s.inverter_efficiency_pct,
            s.rectifier_efficiency_pct,
            s.soc_min_pct,
            s.soc_init_pct,
            s.can_grid_charge,
            s.installed_cost_per_kw,
            s.installed_cost_per_kwh,
            replace_cost_per_kw,
            replace_cost_per_kwh,
            s.inverter_replacement_year,
            s.battery_replacement_year,
            s.macrs_option_years,
            s.macrs_bonus_pct,
            s.macrs_itc_reduction,
            s.total_itc_pct,
            s.total_rebate_per_kw,
            s.total_rebate_per_kwh,
            s.charge_efficiency,
            s.discharge_efficiency,
            s.grid_charge_efficiency,
            net_present_cost_per_kw,
            net_present_cost_per_kwh,
            s.model_degradation,
            degr,
            s.minimum_avg_soc_fraction
        )
    end
end

max_stars_repo_name: NREL/REopt

max_issues_repo_name: NREL/REopt

max_forks_repo_name: NREL/REopt


text:
import Base.==
import Base.isequal
import Base.hash

abstract type AbsEnvironment end
abstract type AbsAction end
abstract type AbsState end

#returns successor states with their probabilities and related rewards
getSuccessors(state::AbsState, env::AbsEnvironment) = error("getSuccessors is unimplemented")
getSuccessors(state::AbsState, action::AbsAction, env::AbsEnvironment) = error("getSuccessors is unimplemented")
isTerminal(state::AbsState, env::AbsEnvironment) = error("isTerminal is unimplemented")
getActions(state::AbsState, env::AbsEnvironment) = error("getActions is unimplemented")

#takes a state and an action, returns an environment info
transfer(env::AbsEnvironment, state::AbsState, action::AbsAction) = error("transfer is unimplemented")
getInitialState(env::AbsEnvironment) = error("getInitialState is unimplemented")
getAllStates(env::AbsEnvironment) = error("getAllStates unimplemented")

render_env(env::AbsEnvironment) = error("render_env is unimplemented")

==(l::AbsAction, r::AbsAction) = error("== is unimplemented")
isequal(l::AbsAction, r::AbsAction) = error("isequal is unimplemented")
hash(x::AbsAction) = error("hash is unimplemented")

==(l::AbsState, r::AbsState) = error("== is unimplemented")
isequal(l::AbsState, r::AbsState) = error("isequal is unimplemented")
hash(x::AbsState) = error("hash is unimplemented")

max_stars_repo_name: ozanarkancan/ReinforcementLearning

max_issues_repo_name: ozanarkancan/ReinforcementLearning

max_forks_repo_name: ozanarkancan/ReinforcementLearning


text:
## Code for CategoricalArray

import Base: Array, convert, collect, copy, getindex, setindex!, similar, size,
             unique, vcat, in, summary, float, complex, copyto!

# Used for keyword argument default value
_isordered(x::AbstractCategoricalArray) = isordered(x)
_isordered(x::Any) = false

function reftype(sz::Int)
    if sz <= typemax(UInt8)
        return UInt8
    elseif sz <= typemax(UInt16)
        return UInt16
    elseif sz <= typemax(UInt32)
        return UInt32
    else
        return UInt64
    end
end

# This check is only there to print a user-friendly warning before
# a TypeError is thrown due to restrictions in the type signature
function check_supported_eltype(::Type{T}, ::Type{U}) where {T, U}
    T === Symbol &&
    throw(ArgumentError("CategoricalArray no longer supports Symbol as element type "*
                        "as that forces recompiling too many Julia Base methods: " *
                        "use strings instead, e.g. via categorical(string.(x))"))
    T <: Union{SupportedTypes, Missing} ||
        throw(ArgumentError("CategoricalArray only supports " *
                            "AbstractString, AbstractChar and Number element types " *
                            "(got element type $U)"))
end

fixstringtype(T::Type) = T <: SubString || T === AbstractString ? String : T
fixstringtype(T::Union) = Union{fixstringtype(T.a), fixstringtype(T.b)}
fixstringtype(::Type{Union{}}) = Union{}

# Find a narrow type that is supported to hold all elements if possible
function fixtype(A::AbstractArray{T}) where T
    if T <: Union{SupportedTypes, Missing}
        return fixstringtype(T)
    else
        U = fixstringtype(mapreduce(typeof, Base.promote_typejoin, A))
        check_supported_eltype(U, T)
        return U
    end
end

"""
    CategoricalArray{T}(undef, dims::Dims; levels=nothing, ordered=false)
    CategoricalArray{T}(undef, dims::Int...; levels=nothing, ordered=false)

Construct an uninitialized `CategoricalArray` with levels of type
`T <: $SupportedTypes` and dimensions `dims`.

The `levels` keyword argument can be a vector specifying possible values for the data
(this is equivalent to but more efficient than calling [`levels!`](@ref)
on the resulting array).
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

    CategoricalArray{T, N, R}(undef, dims::Dims; levels=nothing, ordered=false)
    CategoricalArray{T, N, R}(undef, dims::Int...; levels=nothing, ordered=false)

Similar to definition above, but uses reference type `R` instead of the default type
(`$DefaultRefType`).

    CategoricalArray(A::AbstractArray; levels=nothing, ordered=false)

Construct a new `CategoricalArray` with the values from `A` and the same element type.

The `levels` keyword argument can be a vector specifying possible values for the data
(this is equivalent to but more efficient than calling [`levels!`](@ref)
on the resulting array).
If `levels` is omitted and the element type supports it, levels are sorted
in ascending order; else, they are kept in their order of appearance in `A`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

If `A` is already a `CategoricalArray`, its levels, orderedness and reference type
are preserved unless explicitly overriden.
"""
function CategoricalArray end

"""
    CategoricalVector{T}(undef, m::Int; levels=nothing, ordered=false)

Construct an uninitialized `CategoricalVector` with levels of type
`T <: $SupportedTypes` and dimensions `dim`.

The `levels` keyword argument can be a vector specifying possible values for the data
(this is equivalent to but more efficient than calling [`levels!`](@ref)
on the resulting array).
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

    CategoricalVector{T, R}(undef, m::Int; levels=nothing, ordered=false)

Similar to definition above, but uses reference type `R` instead of the default type
(`$DefaultRefType`).

    CategoricalVector(A::AbstractVector; levels=nothing, ordered=false)

Construct a `CategoricalVector` with the values from `A` and the same element type.

The `levels` keyword argument can be a vector specifying possible values for the data
(this is equivalent to but more efficient than calling [`levels!`](@ref)
on the resulting array).
If `levels` is omitted and the element type supports it, levels are sorted
in ascending order; else, they are kept in their order of appearance in `A`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

If `A` is already a `CategoricalVector`, its levels, orderedness and reference type
are preserved unless explicitly overriden.
"""
function CategoricalVector end

"""
    CategoricalMatrix{T}(undef, m::Int, n::Int; levels=nothing, ordered=false)

Construct an uninitialized `CategoricalMatrix` with levels of type
`T <: $SupportedTypes` and dimensions `dim`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

    CategoricalMatrix{T, R}(undef, m::Int, n::Int; levels=nothing, ordered=false)

Similar to definition above, but uses reference type `R` instead of the default type
(`$DefaultRefType`).

    CategoricalMatrix(A::AbstractMatrix; levels=nothing, ordered=false)

Construct a `CategoricalMatrix` with the values from `A` and the same element type.

The `levels` keyword argument can be a vector specifying possible values for the data
(this is equivalent to but more efficient than calling [`levels!`](@ref)
on the resulting array).
If `levels` is omitted and the element type supports it, levels are sorted
in ascending order; else, they are kept in their order of appearance in `A`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

If `A` is already a `CategoricalMatrix`, its levels, orderedness and reference type
are preserved unless explicitly overriden.
"""
function CategoricalMatrix end

# UndefInitializer array constructors

CategoricalArray(::UndefInitializer, dims::Int...;
                 levels::Union{AbstractVector, Nothing}=nothing,
                 ordered::Bool=false) =
    CategoricalArray{String}(undef, dims, levels=levels, ordered=ordered)

function CategoricalArray{T, N, R}(::UndefInitializer, dims::NTuple{N,Int};
                                   levels::Union{AbstractVector, Nothing}=nothing,
                                   ordered::Bool=false) where {T, N, R}
    U = leveltype(nonmissingtype(T))
    S = T >: Missing ? Union{U, Missing} : U
    check_supported_eltype(S, T)
    V = CategoricalValue{U, R}
    levs = levels === nothing ? U[] : collect(U, levels)
    CategoricalArray{S, N}(zeros(R, dims), CategoricalPool{U, R, V}(levs, ordered))
end

CategoricalArray{T, N}(::UndefInitializer, dims::NTuple{N,Int};
                       levels::Union{AbstractVector, Nothing}=nothing,
                       ordered::Bool=false) where {T, N} =
    CategoricalArray{T, N, DefaultRefType}(undef, dims, levels=levels, ordered=ordered)
CategoricalArray{T, N}(::UndefInitializer, dims::NTuple{N,Int};
                       levels::Union{AbstractVector, Nothing}=nothing,
                       ordered::Bool=false) where
    {R, T <: Union{Missing, CategoricalValue{<:Any, R}}, N} =
    CategoricalArray{T, N, R}(undef, dims, levels=levels, ordered=ordered)
CategoricalArray{T}(::UndefInitializer, dims::NTuple{N,Int};
                    levels::Union{AbstractVector, Nothing}=nothing,
                    ordered::Bool=false) where {T, N} =
    CategoricalArray{T, N}(undef, dims, levels=levels, ordered=ordered)
CategoricalArray{T, 1}(::UndefInitializer, m::Int;
                       levels::Union{AbstractVector, Nothing}=nothing,
                       ordered::Bool=false) where {T} =
    CategoricalArray{T, 1}(undef, (m,), levels=levels, ordered=ordered)
CategoricalArray{T, 2}(::UndefInitializer, m::Int, n::Int;
                       levels::Union{AbstractVector, Nothing}=nothing,
                       ordered::Bool=false) where {T} =
    CategoricalArray{T, 2}(undef, (m, n), levels=levels, ordered=ordered)
CategoricalArray{T, 1, R}(::UndefInitializer, m::Int;
                          levels::Union{AbstractVector, Nothing}=nothing,
                          ordered=false) where {T, R} =
    CategoricalArray{T, 1, R}(undef, (m,), levels=levels, ordered=ordered)
# R <: Integer is required to prevent default constructor from being called instead
CategoricalArray{T, 2, R}(::UndefInitializer, m::Int, n::Int;
                          levels::Union{AbstractVector, Nothing}=nothing,
                          ordered::Bool=false) where {T, R <: Integer} =
    CategoricalArray{T, 2, R}(undef, (m, n), levels=levels, ordered=ordered)
CategoricalArray{T, 3, R}(::UndefInitializer, m::Int, n::Int, o::Int;
                          levels::Union{AbstractVector, Nothing}=nothing,
                          ordered::Bool=false) where {T, R} =
    CategoricalArray{T, 3, R}(undef, (m, n, o), levels=levels, ordered=ordered)
CategoricalArray{T}(::UndefInitializer, m::Int;
                    levels::Union{AbstractVector, Nothing}=nothing,
                    ordered::Bool=false) where {T} =
    CategoricalArray{T}(undef, (m,), levels=levels, ordered=ordered)
CategoricalArray{T}(::UndefInitializer, m::Int, n::Int;
                    levels::Union{AbstractVector, Nothing}=nothing,
                    ordered::Bool=false) where {T} =
    CategoricalArray{T}(undef, (m, n), levels=levels, ordered=ordered)
CategoricalArray{T}(::UndefInitializer, m::Int, n::Int, o::Int;
                    levels::Union{AbstractVector, Nothing}=nothing,
                    ordered::Bool=false) where {T} =
    CategoricalArray{T}(undef, (m, n, o), levels=levels, ordered=ordered)

CategoricalVector(::UndefInitializer, m::Integer;
                  levels::Union{AbstractVector, Nothing}=nothing,
                  ordered::Bool=false) =
    CategoricalArray(undef, m, levels=levels, ordered=ordered)

CategoricalMatrix(::UndefInitializer, m::Int, n::Int;
                  levels::Union{AbstractVector, Nothing}=nothing,
                  ordered::Bool=false) =
    CategoricalArray(undef, m, n, levels=levels, ordered=ordered)

## Constructors from arrays

function CategoricalArray{T, N, R}(A::CategoricalArray{S, N, Q};
                                   levels::Union{AbstractVector, Nothing}=nothing,
                                   ordered::Bool=_isordered(A)) where {S, T, N, Q, R}
    V = unwrap_catvaluetype(T)
    res = convert(CategoricalArray{V, N, R}, A)
    refs = res.refs === A.refs ? copy(res.refs) : res.refs
    pool = res.pool === A.pool ? copy(res.pool) : res.pool
    ordered!(pool, ordered)
    res = CategoricalArray{V, N}(refs, pool)
    if levels !== nothing
        # Calling levels! is faster than checking beforehand which values are used
        try
            levels!(res, levels)
        catch err
            err isa LevelsException || rethrow(err)
            throw(ArgumentError("encountered value(s) not in specified `levels`: " *
                                "$(setdiff(CategoricalArrays.levels(res), levels))"))
        end
    end
    return res
end

function CategoricalArray{T, N, R}(A::AbstractArray;
                                   levels::Union{AbstractVector, Nothing}=nothing,
                                   ordered::Bool=_isordered(A)) where {T, N, R}
    V = unwrap_catvaluetype(T)
    ordered!(_convert(CategoricalArray{V, N, R}, A, levels=levels), ordered)
end

# From AbstractArray

CategoricalArray{T, N}(A::AbstractArray{S, N};
                       levels::Union{AbstractVector, Nothing}=nothing,
                       ordered::Bool=_isordered(A)) where {S, T, N} =
    CategoricalArray{T, N, DefaultRefType}(A, levels=levels, ordered=ordered)
CategoricalArray{T}(A::AbstractArray{S, N};
                    levels::Union{AbstractVector, Nothing}=nothing,
                    ordered::Bool=_isordered(A)) where {S, T, N} =
    CategoricalArray{T, N}(A, levels=levels, ordered=ordered)
CategoricalArray(A::AbstractArray{T, N};
                 levels::Union{AbstractVector, Nothing}=nothing,
                 ordered::Bool=_isordered(A)) where {T, N} =
    CategoricalArray{fixtype(A), N}(A, levels=levels, ordered=ordered)

CategoricalVector(A::AbstractVector{T};
                  levels::Union{AbstractVector, Nothing}=nothing,
                  ordered::Bool=_isordered(A)) where {T} =
    CategoricalArray{fixtype(A), 1}(A, levels=levels, ordered=ordered)

CategoricalMatrix(A::AbstractMatrix{T};
                  levels::Union{AbstractVector, Nothing}=nothing,
                  ordered::Bool=_isordered(A)) where {T} =
    CategoricalArray{fixtype(A), 2}(A, levels=levels, ordered=ordered)

# From CategoricalArray (preserve R)
CategoricalArray{T, N}(A::CategoricalArray{S, N, R};
                       levels::Union{AbstractVector, Nothing}=nothing,
                       ordered::Bool=_isordered(A)) where {S, T, N, R} =
    CategoricalArray{T, N, R}(A, levels=levels, ordered=ordered)
CategoricalArray{T}(A::CategoricalArray{S, N, R};
                    levels::Union{AbstractVector, Nothing}=nothing,
                    ordered::Bool=_isordered(A)) where {S, T, N, R} =
    CategoricalArray{T, N, R}(A, levels=levels, ordered=ordered)
CategoricalArray(A::CategoricalArray{T, N, R};
                 levels::Union{AbstractVector, Nothing}=nothing,
                 ordered::Bool=_isordered(A)) where {T, N, R} =
    CategoricalArray{T, N, R}(A, levels=levels, ordered=ordered)

CategoricalVector(A::CategoricalArray{T, 1, R};
                  levels::Union{AbstractVector, Nothing}=nothing,
                  ordered::Bool=_isordered(A)) where {T, R} =
    CategoricalArray{T, 1, R}(A, levels=levels, ordered=ordered)

CategoricalMatrix(A::CategoricalArray{T, 2, R};
                  levels::Union{AbstractVector, Nothing}=nothing,
                  ordered::Bool=_isordered(A)) where {T, R} =
    CategoricalArray{T, 2, R}(A, levels=levels, ordered=ordered)

## Conversion methods

# From AbstractArray
convert(::Type{CategoricalArray{T, N}}, A::AbstractArray{S, N}) where {S, T, N} =
    convert(CategoricalArray{T, N, DefaultRefType}, A)
convert(::Type{CategoricalArray{T}}, A::AbstractArray{S, N}) where {S, T, N} =
    convert(CategoricalArray{T, N}, A)
convert(::Type{CategoricalArray}, A::AbstractArray{T, N}) where {T, N} =
    convert(CategoricalArray{fixtype(A), N}, A)

convert(::Type{CategoricalVector{T}}, A::AbstractVector) where {T} =
    convert(CategoricalVector{T, DefaultRefType}, A)
convert(::Type{CategoricalVector}, A::AbstractVector{T}) where {T} =
    convert(CategoricalVector{fixtype(A)}, A)
convert(::Type{CategoricalVector{T}},
        A::CategoricalVector{S, R}) where {S, T, R <: Integer} =
    convert(CategoricalVector{T, R}, A)
convert(::Type{CategoricalVector{T}}, A::CategoricalVector{T}) where {T} = A
convert(::Type{CategoricalVector}, A::CategoricalVector) = A

convert(::Type{CategoricalMatrix{T}}, A::AbstractMatrix) where {T} =
    convert(CategoricalMatrix{T, DefaultRefType}, A)
convert(::Type{CategoricalMatrix}, A::AbstractMatrix{T}) where {T} =
    convert(CategoricalMatrix{fixtype(A)}, A)
convert(::Type{CategoricalMatrix{T}},
        A::CategoricalMatrix{S, R}) where {S, T, R <: Integer} =
    convert(CategoricalMatrix{T, R}, A)
convert(::Type{CategoricalMatrix{T}}, A::CategoricalMatrix{T}) where {T} = A
convert(::Type{CategoricalMatrix}, A::CategoricalMatrix) = A

convert(::Type{CategoricalArray{T, N, R}}, A::AbstractArray{S, N}) where {S, T, N, R} =
    _convert(CategoricalArray{T, N, R}, A)

function _convert(::Type{CategoricalArray{T, N, R}}, A::AbstractArray{S, N};
                  levels::Union{AbstractVector, Nothing}=nothing) where {S, T, N, R}
    check_supported_eltype(T, T)

    res = CategoricalArray{T, N, R}(undef, size(A), levels=levels)
    copyto!(res, A)

    if levels !== nothing
        CategoricalArrays.levels(res) == levels ||
            throw(ArgumentError("encountered value(s) not in specified `levels`: " *
                                "$(setdiff(CategoricalArrays.levels(res), levels))"))
    else
        # if order is defined for level type, automatically apply it
        L = leveltype(res)
        if Base.OrderStyle(L) isa Base.Ordered
            levels!(res, sort(CategoricalArrays.levels(res)))
        elseif hasmethod(isless, (L, L))
            # isless may throw an error, e.g. for AbstractArray{T} of unordered T
            try
                levels!(res, sort(CategoricalArrays.levels(res)))
            catch e
                 e isa MethodError || rethrow(e)
            end
        end
    end

    res
end

# From CategoricalArray (preserve levels, ordering and R)
function convert(::Type{CategoricalArray{T, N, R}}, A::CategoricalArray{S, N}) where {S, T, N, R}
    if length(A.pool) > typemax(R)
        throw(LevelsException{T, R}(levels(A)[typemax(R)+1:end]))
    end

    if !(T >: Missing) && S >: Missing && any(iszero, A.refs)
        throw(MissingException("cannot convert CategoricalArray with missing values to a CategoricalArray{$T}"))
    end

    pool = convert(CategoricalPool{unwrap_catvaluetype(nonmissingtype(T)), R}, A.pool)
    refs = convert(Array{R, N}, A.refs)
    CategoricalArray{unwrap_catvaluetype(T), N}(refs, pool)
end

convert(::Type{CategoricalArray{T, N}},
        A::CategoricalArray{S, N, R}) where {S, T, N, R <: Integer} =
    convert(CategoricalArray{T, N, R}, A)
convert(::Type{CategoricalArray{T}},
        A::CategoricalArray{S, N, R}) where {S, T, N, R <: Integer} =
    convert(CategoricalArray{T, N, R}, A)
convert(::Type{CategoricalArray}, A::CategoricalArray{T, N, R}) where {T, N, R} =
    convert(CategoricalArray{T, N, R}, A)

# R<:Integer is needed for this method to be considered more specific
# than the generic one above (JuliaLang/julia#18443)
convert(::Type{CategoricalArray{T, N, R}}, A::CategoricalArray{T, N, R}) where {T, N, R<:Integer} = A
convert(::Type{CategoricalArray{T, N}}, A::CategoricalArray{T, N}) where {T, N} = A
convert(::Type{CategoricalArray{T}}, A::CategoricalArray{T}) where {T} = A
convert(::Type{CategoricalArray}, A::CategoricalArray) = A

function Base.:(==)(A::CategoricalArray{S}, B::CategoricalArray{T}) where {S, T}
    if size(A) != size(B)
        return false
    end
    anymissing = false
    if A.pool === B.pool
        @inbounds for (a, b) in zip(A.refs, B.refs)
            if a == 0 || b == 0
                (S >: Missing || T >: Missing) && (anymissing = true)
            elseif a != b
                return false
            end
        end
    else
        @inbounds for (a, b) in zip(A, B)
            eq = (a == b)
            if eq === false
                return false
            elseif S >: Missing || T >: Missing
                anymissing |= ismissing(eq)
            end
        end
    end
    return anymissing ? missing : true
end

function Base.isequal(A::CategoricalArray, B::CategoricalArray)
    if size(A) != size(B)
        return false
    end
    if A.pool === B.pool
        @inbounds for (a, b) in zip(A.refs, B.refs)
            if a != b
                return false
            end
        end
    else
        @inbounds for (a, b) in zip(A, B)
            if !isequal(a, b)
                return false
            end
        end
    end
    return true
end

size(A::CategoricalArray) = size(A.refs)
Base.IndexStyle(::Type{<:CategoricalArray}) = IndexLinear()

function update_refs!(A::CategoricalArray, newlevels::AbstractVector)
    oldlevels = levels(A)
    levelsmap = similar(A.refs, length(oldlevels)+1)
    # 0 maps to a missing value
    levelsmap[1] = 0
    levelsmap[2:end] .= something.(indexin(oldlevels, newlevels), 0)

    refs = A.refs
    @inbounds for (i, x) in enumerate(refs)
        refs[i] = levelsmap[x+1]
    end
    A
end

function merge_pools!(A::CatArrOrSub,
                      B::Union{CategoricalValue, CatArrOrSub};
                      updaterefs::Bool=true,
                      updatepool::Bool=true)
    newlevels, ordered = merge_pools(pool(A), pool(B))
    oldlevels = levels(A)
    pA = A isa SubArray ? parent(A) : A
    ordered!(pA, ordered)
    # If A's levels are an ordered superset of new (merged) pool, no need to recompute refs
    if updaterefs &&
        (length(newlevels) < length(oldlevels) ||
         view(newlevels, 1:length(oldlevels)) != oldlevels)
        update_refs!(pA, newlevels)
    end
    if updatepool
        pA.pool = typeof(pA.pool)(newlevels, ordered)
    end
    newlevels, ordered
end

@inline function setindex!(A::CategoricalArray, v::Any, I::Real...)
    @boundscheck checkbounds(A, I...)
    if v isa CategoricalValue && pool(v) !== pool(A) && pool(v) ⊈ pool(A)
        merge_pools!(A, v)
    end
    @inbounds A.refs[I...] = get!(A.pool, v)
end

Base.fill(v::CategoricalValue{T}, dims::NTuple{N, Integer}) where {T, N} =
    CategoricalArray{T, N}(fill(refcode(v), dims), copy(pool(v)))

# to avoid ambiguity
Base.fill(v::CategoricalValue, dims::Tuple{}) =
    invoke(fill, Tuple{CategoricalValue{T}, NTuple{N, Integer}} where {T, N}, v, dims)

function Base.fill!(A::CategoricalArray, v::Any)
    if v isa CategoricalValue && pool(v) !== pool(A) && pool(v) ⊈ pool(A)
        merge_pools!(A, v, updaterefs=false)
    end
    fill!(A.refs, get!(A.pool, v))
    A
end

# Methods preserving levels and more efficient than AbstractArray fallbacks
copy(A::CategoricalArray{T, N}) where {T, N} =
    CategoricalArray{T, N}(copy(A.refs), copy(A.pool))

function copyto!(dest::CatArrOrSub{T, N, R}, dstart::Integer,
                 src::CatArrOrSub{<:Any, N}, sstart::Integer,
                 n::Integer) where {T, N, R}
    n < 0 && throw(ArgumentError("tried to copy n=$n elements, but n should be nonnegative"))
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    if n > 0
        (dstart ∈ destinds && dstart+n-1 ∈ destinds) || throw(BoundsError(dest, dstart:dstart+n-1))
        (sstart ∈ srcinds  && sstart+n-1 ∈ srcinds)  || throw(BoundsError(src,  sstart:sstart+n-1))
    end

    drefs = refs(dest)
    srefs = refs(src)
    dpool = pool(dest)
    spool = pool(src)

    # try converting src to dest type to avoid partial copy corruption of dest
    # in the event that the src cannot be copied into dest
    slevs = convert(Vector{T}, levels(src))
    dlevs = levels(dest)
    if eltype(src) >: Missing && !(eltype(dest) >: Missing) && !all(x -> x > 0, srefs)
        throw(MissingException("cannot copy array with missing values to an array with element type $T"))
    end

    destp = dest isa SubArray ? parent(dest) : dest

    # For partial copy, need to recompute existing refs
    # TODO: for performance, avoid ajusting refs which are going to be overwritten
    updaterefs = isa(dest, SubArray) || !(n == length(dest) == length(src))
    newlevels, ordered = merge_pools!(dest, src, updaterefs=updaterefs, updatepool=false)

    # If destination levels are an ordered superset of source, no need to recompute refs
    if view(newlevels, 1:length(slevs)) == slevs
        copyto!(drefs, dstart, srefs, sstart, n)
    else # Otherwise, recompute refs according to new levels
        # Then adjust refs from source
        levelsmap = similar(drefs, length(slevs)+1)
        # 0 maps to a missing value
        levelsmap[1] = 0
        levelsmap[2:end] = indexin(slevs, newlevels)

        @inbounds for i = 0:(n-1)
            x = srefs[sstart+i]
            drefs[dstart+i] = levelsmap[x+1]
        end
    end
    # Need to allocate a new pool only if reordering destination levels
    if view(newlevels, 1:length(dlevs)) == dlevs
        levels!(dpool, newlevels, checkunique=false)
    else
        destp.pool = CategoricalPool{nonmissingtype(T), R}(newlevels, ordered)
    end

    dest
end

function copyto!(dest::CatArrOrSub{T1, N, R}, dstart::Integer,
                 src::AbstractArray{T2, N}, sstart::Integer,
                 n::Integer) where {T1, T2, N, R}
    n == 0 && return dest
    n < 0 && throw(ArgumentError("tried to copy n=$n elements, but n should be nonnegative"))
    destinds, srcinds = LinearIndices(dest), LinearIndices(src)
    (checkbounds(Bool, destinds, dstart) && checkbounds(Bool, destinds, dstart+n-1)) || throw(BoundsError(dest, dstart:dstart+n-1))
    (checkbounds(Bool, srcinds, sstart)  && checkbounds(Bool, srcinds, sstart+n-1))  || throw(BoundsError(src,  sstart:sstart+n-1))
    srclevs = DataAPI.refpool(src)
    srcrefs = DataAPI.refarray(src)
    # Fast path only supports refs which can be used to index
    # into vector mapping from source levels to destination levels
    if !(srclevs isa AbstractVector) || !(srcrefs isa AbstractArray{<:Integer})
        return invoke(copyto!, Tuple{AbstractArray, Integer, AbstractArray, Integer, Integer},
                      dest, dstart, src, sstart, n)
    end
    newdestlevs = destlevs = copy(levels(dest)) # copy since we need original levels below
    srclevsnm = T2 >: Missing ? setdiff(srclevs, [missing]) : srclevs
    if !(srclevsnm ⊆ destlevs)
        # if order is defined for level type, automatically apply it
        L = nonmissingtype(eltype(srclevsnm))
        if hasmethod(isless, Tuple{L, L})
            srclevsnm = srclevsnm === srclevs ? sort(srclevsnm) : sort!(srclevsnm)
        end
        newdestlevs = union(destlevs, srclevsnm)
        levels!(pool(dest), newdestlevs, checkunique=false)
    end
    levelsmap = something.(indexin(srclevs, [missing; newdestlevs])) .- 1
    destrefs = refs(dest)
    seen = fill(false, length(newdestlevs)+1)
    firstind = firstindex(srclevs)
    @inbounds for i in 0:(n-1)
        j = srcrefs[sstart+i] - firstind + 1
        ref = levelsmap[j]
        seen[ref+1] = true
        if !(T1 >: Missing) && T2 >: Missing && ref == 0
            throw(MethodError(convert, (T1, missing)))
        end
        destrefs[dstart+i] = ref
    end
    seennm = @view seen[2:end]
    if !all(seennm)
        destlevsset = Set(destlevs)
        keptlevs = [l for (i, l) in enumerate(newdestlevs)
                    if seennm[i] || l in destlevsset]
        levels!(dest, keptlevs)
    end
    dest
end

# This uses linear indexing even for IndexCartesian src, but
# the performance impact should be modest compared to the dict lookup
copyto!(dest::CatArrOrSub, src::AbstractArray) =
    copyto!(dest, 1, src, 1, length(src))

copyto!(dest::CatArrOrSub, dstart::Integer, src::AbstractArray) =
    copyto!(dest, dstart, src, 1, length(src))

if VERSION >= v"1.1"
    import Base: copy!
else
    import Future: copy!
end
copy!(dest::CatArrOrSub, src::AbstractArray) = copyto!(dest, 1, src, 1, length(src))
# To fix ambiguities
copy!(dest::CatArrOrSub{<:Any, 1}, src::AbstractArray{<:Any, 1}) =
    copyto!(dest, 1, src, 1, length(src))
copy!(dest::CatArrOrSub{T, 1}, src::AbstractArray{T, 1}) where {T} =
    copyto!(dest, 1, src, 1, length(src))

similar(A::CategoricalArray{S, M, R}, ::Type{T},
        dims::NTuple{N, Int}) where {T, N, S, M, R} =
    Array{T, N}(undef, dims)
similar(A::CategoricalArray{S, M, R}, ::Type{Missing},
        dims::NTuple{N, Int}) where {N, S, M, R} =
    Array{Missing, N}(missing, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{CategoricalValue{T, R}},
        dims::NTuple{N, Int}) where {R, T, N, S, M, Q} =
    CategoricalArray{T, N, R}(undef, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{CategoricalValue{T}},
        dims::NTuple{N, Int}) where {T, N, S, M, Q} =
    CategoricalArray{T, N, Q}(undef, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{Union{CategoricalValue{T, R}, Missing}},
        dims::NTuple{N, Int}) where {R, T, N, S, M, Q} =
    CategoricalArray{Union{T, Missing}, N, R}(undef, dims)
similar(A::CategoricalArray{S, M, Q}, ::Type{Union{CategoricalValue{T}, Missing}},
        dims::NTuple{N, Int}) where {T, N, S, M, Q} =
    CategoricalArray{Union{T, Missing}, N, Q}(undef, dims)

# AbstractRange methods are needed since collect uses 1:1 as dummy array
for A in (:Array, :Vector, :Matrix, :AbstractRange)
    @eval begin
        similar(A::$A, ::Type{CategoricalValue{T, R}},
                dims::NTuple{N, Int}=size(A)) where {T, R, N} =
            CategoricalArray{T, N, R}(undef, dims)
        similar(A::$A, ::Type{CategoricalValue{T}},
                dims::NTuple{N, Int}=size(A)) where {T, N} =
            CategoricalArray{T, N}(undef, dims)
        similar(A::$A, ::Type{Union{CategoricalValue{T, R}, Missing}},
                dims::NTuple{N, Int}=size(A)) where {T, R, N} =
            CategoricalArray{Union{T, Missing}, N, R}(undef, dims)
        similar(A::$A, ::Type{Union{CategoricalValue{T}, Missing}},
                dims::NTuple{N, Int}=size(A)) where {T, N} =
            CategoricalArray{Union{T, Missing}, N}(undef, dims)
    end
end

similar(::Type{T}, dims::Dims) where {U, R, T<:Array{CategoricalValue{U, R}}} =
    CategoricalArray{eltype(T)}(undef, dims)
similar(::Type{T}, dims::Dims) where {U, T<:Array{CategoricalValue{U}}} =
    CategoricalArray{eltype(T)}(undef, dims)
similar(::Type{T}, dims::Dims) where {U, R, T<:Array{Union{CategoricalValue{U, R}, Missing}}} =
    CategoricalArray{eltype(T)}(undef, dims)
similar(::Type{T}, dims::Dims) where {U, T<:Array{Union{CategoricalValue{U}, Missing}}} =
    CategoricalArray{eltype(T)}(undef, dims)

"""
    compress(A::CategoricalArray)

Return a copy of categorical array `A` using the smallest reference type able to hold the
number of [`levels`](@ref DataAPI.levels) of `A`.

While this will reduce memory use, this function is type-unstable, which can affect
performance inside the function where the call is made. Therefore, use it with caution.
"""
function compress(A::CategoricalArray{T, N}) where {T, N}
    R = reftype(length(levels(A.pool)))
    convert(CategoricalArray{T, N, R}, A)
end

"""
    decompress(A::CategoricalArray)

Return a copy of categorical array `A` using the default reference type ($DefaultRefType).
If `A` is using a small reference type (such as `UInt8` or `UInt16`) the decompressed array
will have room for more levels.

To avoid the need to call decompress, ensure [`compress`](@ref) is not called when creating
the categorical array.
"""
decompress(A::CategoricalArray{T, N}) where {T, N} =
    convert(CategoricalArray{T, N, DefaultRefType}, A)

function vcat(A::CategoricalArray...)
    ordered = any(isordered, A) && all(a->isordered(a) || isempty(levels(a)), A)
    newlevels, ordered = mergelevels(ordered, map(levels, A)...)

    refsvec = map(A) do a
        ii = convert(Vector{Int}, indexin(levels(a.pool), newlevels))
        [x==0 ? 0 : ii[x] for x in a.refs]::Array{Int,ndims(a)}
    end

    T = cat_promote_eltype(A...) >: Missing ?
        Union{eltype(newlevels), Missing} : eltype(newlevels)
    refs = DefaultRefType[refsvec...;]
    pool = CategoricalPool(newlevels, ordered)
    CategoricalArray{T, ndims(refs)}(refs, pool)
end

@inline function getindex(A::CategoricalArray{T}, I...) where {T}
    @boundscheck checkbounds(A, I...)
    # Let Array indexing code handle everything
    @inbounds r = A.refs[I...]

    if isa(r, Array)
        res = CategoricalArray{T, ndims(r)}(r, copy(A.pool))
        return ordered!(res, isordered(A))
    else
        r > 0 || throw(UndefRefError())
        @inbounds res = A.pool[r]
        return res
    end
end

leveltype(::Type{T}) where {T <: CategoricalArray} = leveltype(nonmissingtype(eltype(T)))

"""
    levels(x::CategoricalArray)
    levels(x::CategoricalValue)

Return the levels of categorical array or value `x`.
This may include levels which do not actually appear in the data
(see [`droplevels!`](@ref)).
"""
DataAPI.levels(A::CategoricalArray) = levels(A.pool)

"""
    levels!(A::CategoricalArray, newlevels::Vector; allowmissing::Bool=false)

Set the levels categorical array `A`. The order of appearance of levels will be respected
by [`levels`](@ref DataAPI.levels), which may affect display of results in some operations; if `A` is
ordered (see [`isordered`](@ref)), it will also be used for order comparisons
using `<`, `>` and similar operators. Reordering levels will never affect the values
of entries in the array.

If `A` accepts missing values (i.e. `eltype(A) >: Missing`) and `allowmissing=true`,
entries corresponding to omitted levels will be set to `missing`.
Else, `newlevels` must include all levels which appear in the data.
"""
function levels!(A::CategoricalArray{T, N, R}, newlevels::Vector;
                 allowmissing::Bool=false,
                 allow_missing::Union{Bool, Nothing}=nothing) where {T, N, R}
    if allow_missing !== nothing
        Base.depwarn("allow_missing argument is deprecated, use allowmissing instead",
                     :levels!)
        allowmissing = allow_missing
    end
    if !allunique(newlevels)
        throw(ArgumentError(string("duplicated levels found: ",
                                   join(unique(filter(x->sum(newlevels.==x)>1, newlevels)), ", "))))
    end

    oldlevels = levels(A.pool)

    # first pass to check whether, if some levels are removed, changes can be applied without error
    # TODO: save original levels and undo changes in case of error to skip this step
    # equivalent to issubset but faster due to JuliaLang/julia#24624
    if !isempty(setdiff(oldlevels, newlevels))
        deleted = [!(l in newlevels) for l in oldlevels]
        @inbounds for (i, x) in enumerate(A.refs)
            if T >: Missing
                !allowmissing && x > 0 && deleted[x] &&
                    throw(ArgumentError("cannot remove level $(repr(oldlevels[x])) as it " *
                                        "is used at position $i and allowmissing=false."))
            else
                x > 0 && deleted[x] &&
                    throw(ArgumentError("cannot remove level $(repr(oldlevels[x])) as it " *
                                        "is used at position $i. Change the array element " *
                                        "type to Union{$T, Missing} using convert if you want " *
                                        "to transform some levels to missing values."))
            end
        end
    end

    # replace the pool and recode refs to reflect new pool
    if newlevels != oldlevels
        newpool = CategoricalPool{nonmissingtype(T), R}(copy(newlevels), isordered(A.pool))
        update_refs!(A, newlevels)
        A.pool = newpool
    end

    A
end

function _unique(::Type{S},
                 refs::AbstractArray{T},
                 pool::CategoricalPool) where {S, T<:Integer}
    nlevels = length(levels(pool)) + 1
    order = fill(0, nlevels) # 0 indicates not seen
    # If we don't track missings, short-circuit even if none has been seen
    count = S >: Missing ? 0 : 1
    @inbounds for i in refs
        if order[i + 1] == 0
            count += 1
            order[i + 1] = count
            count == nlevels && break
        end
    end
    S[i == 1 ? missing : levels(pool)[i - 1] for i in sortperm(order) if order[i] != 0]
end

"""
    unique(A::CategoricalArray)

Return levels which appear in `A` in their order of appearance.
This function is significantly slower than [`levels`](@ref DataAPI.levels)
since it needs to check whether levels are used or not.
"""
unique(A::CategoricalArray{T}) where {T} = _unique(T, A.refs, A.pool)

"""
    droplevels!(A::CategoricalArray)

Drop levels which do not appear in categorical array `A` (so that they will no longer be
returned by [`levels`](@ref DataAPI.levels)).
"""
droplevels!(A::CategoricalArray) = levels!(A, intersect(levels(A), unique(A)))

"""
    isordered(A::CategoricalArray)

Test whether entries in `A` can be compared using `<`, `>` and similar operators,
using the ordering of levels.
"""
isordered(A::CategoricalArray) = isordered(A.pool)

"""
    ordered!(A::CategoricalArray, ordered::Bool)

Set whether entries in `A` can be compared using `<`, `>` and similar operators,
using the ordering of levels. Return the modified `A`.
"""
ordered!(A::CategoricalArray, ordered) = (ordered!(A.pool, ordered); return A)

function Base.resize!(A::CategoricalVector, n::Integer)
    n_orig = length(A)
    resize!(A.refs, n)
    if n > n_orig
        A.refs[n_orig+1:end] .= 0
    end
    A
end

function Base.push!(A::CategoricalVector, v::Any)
    if v isa CategoricalValue && pool(v) !== pool(A) && pool(v) ⊈ pool(A)
        merge_pools!(A, v)
    end
    r = get!(A.pool, v)
    push!(A.refs, r)
    A
end

function Base.append!(A::CategoricalVector, B::CatArrOrSub)
    if pool(B) !== pool(A) && pool(B) ⊈ pool(A)
        merge_pools!(A, B)
    end
    # TODO: optimize recoding
    len = length(A)
    len2 = length(B)
    resize!(A.refs, len + len2)
    for i = 1:len2
        A[len + i] = B[i]
    end
    return A
end

Base.empty!(A::CategoricalArray) = (empty!(A.refs); return A)

function Base.reshape(A::CategoricalArray{T, N}, dims::Dims) where {T, N}
    x = reshape(A.refs, dims)
    res = CategoricalArray{T, ndims(x)}(x, A.pool)
    ordered!(res, isordered(res))
end

"""
    categorical(A::AbstractArray; levels=nothing, ordered=false, compress=false)

Construct a categorical array with the values from `A`.

The `levels` keyword argument can be a vector specifying possible values for the data
(this is equivalent to but more efficient than calling [`levels!`](@ref)
on the resulting array).
If `levels` is omitted and the element type supports it, levels are sorted
in ascending order; else, they are kept in their order of appearance in `A`.
The `ordered` keyword argument determines whether the array values can be compared
according to the ordering of levels or not (see [`isordered`](@ref)).

If `compress` is `true`, the smallest reference type able to hold the
number of unique values in `A` will be used. While this will reduce memory use, passing
this parameter will also introduce a type instability which can affect performance inside
the function where the call is made. Therefore, use this option with caution (the
one-argument version does not suffer from this problem).

    categorical(A::CategoricalArray; compress=false, levels=nothing, ordered=false)

If `A` is already a `CategoricalArray`, its levels, orderedness and reference type
are preserved unless explicitly overriden.
"""
@inline function categorical(A::AbstractArray{T, N};
                             levels::Union{AbstractVector, Nothing}=nothing,
                             ordered=_isordered(A),
                             compress::Bool=false) where {T, N}
    # @inline is needed so that return type is inferred when compress is not provided
    RefType = compress ? reftype(length(unique(A))) : DefaultRefType
    CategoricalArray{fixtype(A), N, RefType}(A, levels=levels, ordered=ordered)
end
@inline function categorical(A::CategoricalArray{T, N, R};
                             levels::Union{AbstractVector, Nothing}=nothing,
                             ordered=_isordered(A),
                             compress::Bool=false) where {T, N, R}
    # @inline is needed so that return type is inferred when compress is not provided
    RefType = compress ? reftype(length(CategoricalArrays.levels(A))) : R
    CategoricalArray{T, N, RefType}(A, levels=levels, ordered=ordered)
end

function in(x::Any, y::CategoricalArray{T, N, R}) where {T, N, R}
    ref = get(y.pool, x, zero(R))
    ref != 0 ? ref in y.refs : false
end

function in(x::CategoricalValue, y::CategoricalArray{T, N, R}) where {T, N, R}
    if x.pool === y.pool
        return refcode(x) in y.refs
    else
        ref = get(y.pool, levels(x.pool)[refcode(x)], zero(R))
        return ref != 0 ? ref in y.refs : false
    end
end

Array(A::CategoricalArray{T}) where {T} = Array{T}(A)
collect(A::CategoricalArray) = copy(A)

# Defined for performance
collect(x::Base.SkipMissing{<: CatArrOrSub{T}}) where {T} =
    CategoricalVector{nonmissingtype(T)}(filter(v -> v > 0, refs(x.x)), copy(pool(x.x)))

function float(A::CategoricalArray{T}) where T
    if !isconcretetype(T)
        error("`float` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(float(zero(T)))}, A)
end
function complex(A::CategoricalArray{T}) where T
    if !isconcretetype(T)
        error("`complex` not defined on abstractly-typed arrays; please convert to a more specific type")
    end
    convert(AbstractArray{typeof(complex(zero(T)))}, A)
end

# Override AbstractArray method to avoid printing useless type parameters
if VERSION >= v"0.7.0-DEV.2657"
    summary(io::IO, A::CategoricalArray{T, N, R}) where {T, N, R} =
        print(io, Base.dims2string(size(A)), " $CategoricalArray{$T,$N,$R}")
else
    summary(A::CategoricalArray{T, N, R}) where {T, N, R} =
        string(Base.dims2string(size(A)), " $CategoricalArray{$T,$N,$R}")
end

refs(A::CategoricalArray) = A.refs
pool(A::CategoricalArray) = A.pool

Base.deleteat!(A::CategoricalArray, inds) = (deleteat!(A.refs, inds); A)

Base.Broadcast.broadcasted(::typeof(ismissing), A::CategoricalArray{T}) where {T} =
    T >: Missing ? Base.Broadcast.broadcasted(==, A.refs, 0) :
                   Base.Broadcast.broadcasted(_ -> false, A.refs)

Base.Broadcast.broadcasted(::typeof(!ismissing), A::CategoricalArray{T}) where {T} =
    T >: Missing ? Base.Broadcast.broadcasted(>, A.refs, 0) :
                   Base.Broadcast.broadcasted(_ -> true, A.refs)

function Base.Broadcast.broadcasted(::typeof(levelcode), A::CategoricalArray{T}) where {T}
    if T >: Missing
        Base.Broadcast.broadcasted(r -> r > 0 ? Signed(widen(r)) : missing, A.refs)
    else
        Base.Broadcast.broadcasted(r -> Signed(widen(r)), A.refs)
    end
end

function Base.sort!(v::CategoricalVector;
                    # alg is ignored since counting sort is more efficient
                    alg::Base.Algorithm=Base.Sort.defalg(v),
                    lt=isless,
                    by=identity,
                    rev::Bool=false,
                    order::Base.Ordering=Base.Forward)
    counts = zeros(UInt, length(v.pool) + (eltype(v) >: Missing))

    # do a count/histogram of the references
    @inbounds for ref in v.refs
        counts[ref + (eltype(v) >: Missing)] += 1
    end

    # compute the order in which to read from counts
    ord = Base.Sort.ord(lt, by, rev, order)
    seen = counts .> 0
    anymissing = eltype(v) >: Missing && seen[1]
    levs = eltype(v) >: Missing ?
        eltype(v)[i == 0 ? missing : CategoricalValue(v.pool, i) for i in 0:length(v.pool)] :
        eltype(v)[CategoricalValue(v.pool, i) for i in 1:length(v.pool)]
    sortedlevs = sort!(Vector(view(levs, seen)), order=ord)
    levelsmap = something.(indexin(sortedlevs, levs))
    j = 0
    refs = v.refs
    @inbounds for i in eachindex(levelsmap)
        ref = levelsmap[i]
        tmpj = j + counts[ref]
        refs[(j+1):tmpj] .= ref - (eltype(v) >: Missing)
        j = tmpj
    end

    return v
end

Base.repeat(a::CatArrOrSub{T, N},
            counts::Integer...) where {T, N} =
    CategoricalArray{T, N}(repeat(refs(a), counts...), copy(pool(a)))

Base.repeat(a::CatArrOrSub{T, N};
            inner = nothing, outer = nothing) where {T, N} =
    CategoricalArray{T, N}(repeat(refs(a), inner=inner, outer=outer), copy(pool(a)))

# JSON3 writing/reading
StructTypes.StructType(::Type{<:CategoricalVector}) = StructTypes.ArrayType()

StructTypes.construct(::Type{<:CategoricalArray}, A::AbstractVector) =
    constructgeneral(A)
StructTypes.construct(::Type{<:CategoricalArray}, A::Vector) =
    constructgeneral(A)

function constructgeneral(A)
    if eltype(A) === Any
        # unlike `replace`, broadcast narrows the type, which allows us to return small
        # union eltypes (e.g. Union{String,Missing})
        categorical(ifelse.(A .=== nothing, missing, A))
    elseif eltype(A) >: Nothing
        categorical(replace(A, nothing=>missing))
    else
        categorical(A)
    end
end

StructTypes.construct(::Type{<:CategoricalArray{Union{Missing, T}}},
                      A::AbstractVector) where {T} =
    CategoricalArray{Union{Missing, T}}(replace(A, nothing=>missing))
StructTypes.construct(::Type{<:CategoricalArray{Union{Missing, T}}},
                      A::Vector) where {T} =
    CategoricalArray{Union{Missing, T}}(replace(A, nothing=>missing))

# DataAPI refarray/refvalue/refpool support
struct CategoricalRefPool{T, P} <: AbstractVector{T}
    pool::P
end

Base.IndexStyle(::Type{<: CategoricalRefPool}) = Base.IndexLinear()

@inline function Base.getindex(x::CategoricalRefPool, i::Int)
    @boundscheck checkbounds(x, i)
    i > 0 ? @inbounds(x.pool[i]) : missing
end

Base.size(x::CategoricalRefPool{T}) where {T} = (length(x.pool) + (T >: Missing),)
Base.axes(x::CategoricalRefPool{T}) where {T} =
    ((T >: Missing ? 0 : 1):length(x.pool),)
Base.LinearIndices(x::CategoricalRefPool) = axes(x, 1)

DataAPI.refarray(A::CatArrOrSub) = refs(A)
DataAPI.refpool(A::CatArrOrSub{T}) where {T} =
    CategoricalRefPool{eltype(A), typeof(pool(A))}(pool(A))
DataAPI.invrefpool(A::CatArrOrSub{T}) where {T} =
    CategoricalInvRefPool{eltype(A), typeof(pool(A).invindex)}(pool(A).invindex)

@inline function DataAPI.refvalue(A::CatArrOrSub{T}, i::Integer) where T
    @boundscheck checkindex(Bool, (T >: Missing ? 0 : 1):length(pool(A)), i) ||
        throw(BoundsError())
    i > 0 ? @inbounds(pool(A)[i]) : missing
end

struct CategoricalInvRefPool{T, P}
    invpool::P
end

@inline function Base.haskey(x::CategoricalInvRefPool{T}, v) where {T}
    if T >: Missing && ismissing(v)
        return true
    else
        return haskey(x.invpool, v)
    end
end

@inline function Base.getindex(x::CategoricalInvRefPool{T}, v) where {T}
    if T >: Missing && ismissing(v)
        return zero(valtype(x.invpool))
    else
        return x.invpool[v]
    end
end

@inline function Base.get(x::CategoricalInvRefPool{T}, v, default) where {T}
    if T >: Missing && ismissing(v)
        return zero(valtype(x.invpool))
    else
        return get(x.invpool, v, default)
    end
end