text:
# This file is a part of Julia. License is MIT: http://julialang.org/license

### Data

spv_x1 = SparseVector(8, [2, 5, 6], [1.25, -0.75, 3.5])

@test isa(spv_x1, SparseVector{Float64,Int})

x1_full = zeros(length(spv_x1))
x1_full[SparseArrays.nonzeroinds(spv_x1)] = nonzeros(spv_x1)

### Basic Properties

let x = spv_x1
    @test eltype(x) == Float64
    @test ndims(x) == 1
    @test length(x) == 8
    @test size(x) == (8,)
    @test size(x,1) == 8
    @test size(x,2) == 1
    @test !isempty(x)

    @test countnz(x) == 3
    @test nnz(x) == 3
    @test SparseArrays.nonzeroinds(x) == [2, 5, 6]
    @test nonzeros(x) == [1.25, -0.75, 3.5]
end

@test count(SparseVector(8, [2, 5, 6], [true,false,true])) == 2

# full

for (x, xf) in [(spv_x1, x1_full)]
    @test isa(Array(x), Vector{Float64})
    @test Array(x) == xf
end

### Show

@test contains(string(spv_x1), "1.25")
@test contains(string(spv_x1), "-0.75")
@test contains(string(spv_x1), "3.5")

### Other Constructors

### Comparison helper to ensure exact equality with internal structure
function exact_equal(x::AbstractSparseVector, y::AbstractSparseVector)
    eltype(x) == eltype(y) &&
    eltype(SparseArrays.nonzeroinds(x)) == eltype(SparseArrays.nonzeroinds(y)) &&
    length(x) == length(y) &&
    SparseArrays.nonzeroinds(x) == SparseArrays.nonzeroinds(y) &&
    nonzeros(x) == nonzeros(y)
end

# construct empty sparse vector

@test exact_equal(spzeros(Float64, 8), SparseVector(8, Int[], Float64[]))

# from list of indices and values

@test exact_equal(
    sparsevec(Int[], Float64[], 8),
    SparseVector(8, Int[], Float64[]))

@test exact_equal(
    sparsevec(Int[], Float64[]),
    SparseVector(0, Int[], Float64[]))

@test exact_equal(
    sparsevec([3, 3], [5.0, -5.0], 8),
    SparseVector(8, [3], [0.0]))

@test exact_equal(
    sparsevec([2, 3, 6], [12.0, 18.0, 25.0]),
    SparseVector(6, [2, 3, 6], [12.0, 18.0, 25.0]))

let x0 = SparseVector(8, [2, 3, 6], [12.0, 18.0, 25.0])
    @test exact_equal(
        sparsevec([2, 3, 6], [12.0, 18.0, 25.0], 8), x0)

    @test exact_equal(
        sparsevec([3, 6, 2], [18.0, 25.0, 12.0], 8), x0)

    @test exact_equal(
        sparsevec([2, 3, 4, 4, 6], [12.0, 18.0, 5.0, -5.0, 25.0], 8),
        SparseVector(8, [2, 3, 4, 6], [12.0, 18.0, 0.0, 25.0]))

    @test exact_equal(
        sparsevec([1, 1, 1, 2, 3, 3, 6], [2.0, 3.0, -5.0, 12.0, 10.0, 8.0, 25.0], 8),
        SparseVector(8, [1, 2, 3, 6], [0.0, 12.0, 18.0, 25.0]))

    @test exact_equal(
        sparsevec([2, 3, 6, 7, 7], [12.0, 18.0, 25.0, 5.0, -5.0], 8),
        SparseVector(8, [2, 3, 6, 7], [12.0, 18.0, 25.0, 0.0]))
end

@test exact_equal(
    sparsevec(Any[1, 3], [1, 1]),
    sparsevec([1, 3], [1, 1]))

@test exact_equal(
    sparsevec(Any[1, 3], [1, 1], 5),
    sparsevec([1, 3], [1, 1], 5))

# from dictionary

function my_intmap(x)
    a = Dict{Int,eltype(x)}()
    for i in SparseArrays.nonzeroinds(x)
        a[i] = x[i]
    end
    return a
end

let x = spv_x1
    a = my_intmap(x)
    xc = sparsevec(a, 8)
    @test exact_equal(x, xc)

    xc = sparsevec(a)
    @test exact_equal(xc, SparseVector(6, [2, 5, 6], [1.25, -0.75, 3.5]))

    d = Dict{Int, Float64}((1 => 0.0, 2 => 1.0, 3 => 2.0))
    @test exact_equal(sparsevec(d), SparseVector(3, [1, 2, 3], [0.0, 1.0, 2.0]))
end

# spones - copies structure, but replaces nzvals with ones
let x = SparseVector(8, [2, 3, 6], [12.0, 18.0, 25.0])
    y = spones(x)
    @test (x .!= 0) == (y .!= 0)
    @test y == SparseVector(8, [2, 3, 6], [1.0, 1.0, 1.0])
end

# sprand & sprandn

let xr = sprand(1000, 0.9)
    @test isa(xr, SparseVector{Float64,Int})
    @test length(xr) == 1000
    @test all(nonzeros(xr) .>= 0.0)
end

let xr = sprand(Float32, 1000, 0.9)
    @test isa(xr, SparseVector{Float32,Int})
    @test length(xr) == 1000
    @test all(nonzeros(xr) .>= 0.0)
end

let xr = sprandn(1000, 0.9)
    @test isa(xr, SparseVector{Float64,Int})
    @test length(xr) == 1000
    if !isempty(nonzeros(xr))
        @test any(nonzeros(xr) .> 0.0) && any(nonzeros(xr) .< 0.0)
    end
end

let xr = sprand(Bool, 1000, 0.9)
    @test isa(xr, SparseVector{Bool,Int})
    @test length(xr) == 1000
    @test all(nonzeros(xr))
end

let r1 = MersenneTwister(), r2 = MersenneTwister()
    @test sprand(r1, 100, .9) == sprand(r2, 100, .9)
    @test sprandn(r1, 100, .9) == sprandn(r2, 100, .9)
    @test sprand(r1, Bool, 100, .9, ) == sprand(r2,  Bool, 100, .9)
end

### Element access

# getindex

# single integer index
for (x, xf) in [(spv_x1, x1_full)]
    for i = 1:length(x)
        @test x[i] == xf[i]
    end
end

# generic array index
let x = sprand(100, 0.5)
    I = rand(1:length(x), 20)
    r = x[I]
    @test isa(r, SparseVector{Float64,Int})
    @test all(nonzeros(r) .!= 0.0)
    @test Array(r) == Array(x)[I]
end

# boolean array index
let x = sprand(10, 10, 0.5)
    I = rand(1:size(x, 2), 10)
    bI = falses(size(x, 2))
    bI[I] = true
    r = x[1,bI]
    @test isa(r, SparseVector{Float64,Int})
    @test all(nonzeros(r) .!= 0.0)
    @test Array(r) == Array(x)[1,bI]
end

let x = sprand(10, 0.5)
    I = rand(1:length(x), 5)
    bI = falses(length(x))
    bI[I] = true
    r = x[bI]
    @test isa(r, SparseVector{Float64,Int})
    @test all(nonzeros(r) .!= 0.0)
    @test Array(r) == Array(x)[bI]
end

# setindex

let xc = spzeros(Float64, 8)
    xc[3] = 2.0
    @test exact_equal(xc, SparseVector(8, [3], [2.0]))
end

let xc = copy(spv_x1)
    xc[5] = 2.0
    @test exact_equal(xc, SparseVector(8, [2, 5, 6], [1.25, 2.0, 3.5]))
end

let xc = copy(spv_x1)
    xc[3] = 4.0
    @test exact_equal(xc, SparseVector(8, [2, 3, 5, 6], [1.25, 4.0, -0.75, 3.5]))

    xc[1] = 6.0
    @test exact_equal(xc, SparseVector(8, [1, 2, 3, 5, 6], [6.0, 1.25, 4.0, -0.75, 3.5]))

    xc[8] = -1.5
    @test exact_equal(xc, SparseVector(8, [1, 2, 3, 5, 6, 8], [6.0, 1.25, 4.0, -0.75, 3.5, -1.5]))
end

let xc = copy(spv_x1)
    xc[5] = 0.0
    @test exact_equal(xc, SparseVector(8, [2, 5, 6], [1.25, 0.0, 3.5]))

    xc[6] = 0.0
    @test exact_equal(xc, SparseVector(8, [2, 5, 6], [1.25, 0.0, 0.0]))

    xc[2] = 0.0
    @test exact_equal(xc, SparseVector(8, [2, 5, 6], [0.0, 0.0, 0.0]))

    xc[1] = 0.0
    @test exact_equal(xc, SparseVector(8, [2, 5, 6], [0.0, 0.0, 0.0]))
end

## dropstored! tests
let x = SparseVector(10, [2, 7, 9], [2.0, 7.0, 9.0])
    # Test argument bounds checking for dropstored!(x, i)
    @test_throws BoundsError Base.SparseArrays.dropstored!(x, 0)
    @test_throws BoundsError Base.SparseArrays.dropstored!(x, 11)
    # Test behavior of dropstored!(x, i)
    # --> Test dropping a single stored entry
    @test Base.SparseArrays.dropstored!(x, 2) == SparseVector(10, [7, 9], [7.0, 9.0])
    # --> Test dropping a single nonstored entry
    @test Base.SparseArrays.dropstored!(x, 5) == SparseVector(10, [7, 9], [7.0, 9.0])
end

# find and findnz tests
@test find(spv_x1) == find(x1_full)
@test findnz(spv_x1) == (find(x1_full), filter(x->x!=0, x1_full))
let xc = SparseVector(8, [2, 3, 5], [1.25, 0, -0.75]), fc = Array(xc)
    @test find(xc) == find(fc)
    @test findnz(xc) == ([2, 5], [1.25, -0.75])
end

### Array manipulation

# copy

let x = spv_x1
    xc = copy(x)
    @test isa(xc, SparseVector{Float64,Int})
    @test x.nzind !== xc.nzval
    @test x.nzval !== xc.nzval
    @test exact_equal(x, xc)
end

let a = SparseVector(8, [2, 5, 6], Int32[12, 35, 72])
    # vec
    @test vec(a) == a

    # reinterpret
    au = reinterpret(UInt32, a)
    @test isa(au, SparseVector{UInt32,Int})
    @test exact_equal(au, SparseVector(8, [2, 5, 6], UInt32[12, 35, 72]))

    # float
    af = float(a)
    @test float(af) == af
    @test isa(af, SparseVector{Float64,Int})
    @test exact_equal(af, SparseVector(8, [2, 5, 6], [12., 35., 72.]))
    @test sparsevec(transpose(transpose(af))) == af

    # complex
    acp = complex(af)
    @test complex(acp) == acp
    @test isa(acp, SparseVector{Complex128,Int})
    @test exact_equal(acp, SparseVector(8, [2, 5, 6], complex([12., 35., 72.])))
    @test sparsevec(ctranspose(ctranspose(acp))) == acp
end

let x1 = SparseVector(8, [2, 5, 6], [12.2, 1.4, 5.0])
    x2 = SparseVector(8, [3, 4], [1.2, 3.4])
    copy!(x2, x1)
    @test x2 == x1
    x2 = SparseVector(8, [2, 4, 8], [10.3, 7.4, 3.1])
    copy!(x2, x1)
    @test x2 == x1
    x2 = SparseVector(8, [1, 3, 4, 7], [0.3, 1.2, 3.4, 0.1])
    copy!(x2, x1)
    @test x2 == x1
    x2 = SparseVector(10, [3, 4], [1.2, 3.4])
    copy!(x2, x1)
    @test x2[1:8] == x1
    @test x2[9:10] == spzeros(2)
    x2 = SparseVector(10, [3, 4, 9], [1.2, 3.4, 17.8])
    copy!(x2, x1)
    @test x2[1:8] == x1
    @test x2[9] == 17.8
    @test x2[10] == 0
    x2 = SparseVector(10, [3, 4, 5, 6, 9], [8.3, 7.2, 1.2, 3.4, 17.8])
    copy!(x2, x1)
    @test x2[1:8] == x1
    @test x2[9] == 17.8
    @test x2[10] == 0
    x2 = SparseVector(6, [3, 4], [1.2, 3.4])
    @test_throws BoundsError copy!(x2, x1)
end

let x1 = sparse([2, 1, 2], [1, 3, 3], [12.2, 1.4, 5.0], 2, 4)
    x2 = SparseVector(8, [3, 4], [1.2, 3.4])
    copy!(x2, x1)
    @test x2[:] == x1[:]
    x2 = SparseVector(8, [2, 4, 8], [10.3, 7.4, 3.1])
    copy!(x2, x1)
    @test x2[:] == x1[:]
    x2 = SparseVector(8, [1, 3, 4, 7], [0.3, 1.2, 3.4, 0.1])
    copy!(x2, x1)
    @test x2[:] == x1[:]
    x2 = SparseVector(10, [3, 4], [1.2, 3.4])
    copy!(x2, x1)
    @test x2[1:8] == x1[:]
    @test x2[9:10] == spzeros(2)
    x2 = SparseVector(10, [3, 4, 9], [1.2, 3.4, 17.8])
    copy!(x2, x1)
    @test x2[1:8] == x1[:]
    @test x2[9] == 17.8
    @test x2[10] == 0
    x2 = SparseVector(10, [3, 4, 5, 6, 9], [8.3, 7.2, 1.2, 3.4, 17.8])
    copy!(x2, x1)
    @test x2[1:8] == x1[:]
    @test x2[9] == 17.8
    @test x2[10] == 0
    x2 = SparseVector(6, [3, 4], [1.2, 3.4])
    @test_throws BoundsError copy!(x2, x1)
end

let x1 = SparseVector(8, [2, 5, 6], [12.2, 1.4, 5.0])
    x2 = sparse([1, 2], [2, 2], [1.2, 3.4], 2, 4)
    copy!(x2, x1)
    @test x2[:] == x1[:]
    x2 = sparse([2, 2, 2], [1, 3, 4], [10.3, 7.4, 3.1], 2, 4)
    copy!(x2, x1)
    @test x2[:] == x1[:]
    x2 = sparse([1, 1, 2, 1], [1, 2, 2, 4], [0.3, 1.2, 3.4, 0.1], 2, 4)
    copy!(x2, x1)
    @test x2[:] == x1[:]
    x2 = sparse([1, 2], [2, 2], [1.2, 3.4], 2, 5)
    copy!(x2, x1)
    @test x2[1:8] == x1
    @test x2[9:10] == spzeros(2)
    x2 = sparse([1, 2, 1], [2, 2, 5], [1.2, 3.4, 17.8], 2, 5)
    copy!(x2, x1)
    @test x2[1:8] == x1
    @test x2[9] == 17.8
    @test x2[10] == 0
    x2 = sparse([1, 2, 1, 2, 1], [2, 2, 3, 3, 5], [8.3, 7.2, 1.2, 3.4, 17.8], 2, 5)
    copy!(x2, x1)
    @test x2[1:8] == x1
    @test x2[9] == 17.8
    @test x2[10] == 0
    x2 = sparse([1, 2], [2, 2], [1.2, 3.4], 2, 3)
    @test_throws BoundsError copy!(x2, x1)
end

### Type conversion

let x = convert(SparseVector, sparse([2, 5, 6], [1, 1, 1], [1.25, -0.75, 3.5], 8, 1))
    @test isa(x, SparseVector{Float64,Int})
    @test exact_equal(x, spv_x1)
end

let x = spv_x1, xf = x1_full
    xc = convert(SparseVector, xf)
    @test isa(xc, SparseVector{Float64,Int})
    @test exact_equal(xc, x)

    xc = convert(SparseVector{Float32,Int}, x)
    xf32 = SparseVector(8, [2, 5, 6], [1.25f0, -0.75f0, 3.5f0])
    @test isa(xc, SparseVector{Float32,Int})
    @test exact_equal(xc, xf32)

    xc = convert(SparseVector{Float32}, x)
    @test isa(xc, SparseVector{Float32,Int})
    @test exact_equal(xc, xf32)

    xm = convert(SparseMatrixCSC, x)
    @test isa(xm, SparseMatrixCSC{Float64,Int})
    @test Array(xm) == reshape(xf, 8, 1)

    xm = convert(SparseMatrixCSC{Float32}, x)
    @test isa(xm, SparseMatrixCSC{Float32,Int})
    @test Array(xm) == reshape(convert(Vector{Float32}, xf), 8, 1)
end

### Concatenation

let m = 80, n = 100
    A = Array{SparseVector{Float64,Int}}(n)
    tnnz = 0
    for i = 1:length(A)
        A[i] = sprand(m, 0.3)
        tnnz += nnz(A[i])
    end

    H = hcat(A...)
    @test isa(H, SparseMatrixCSC{Float64,Int})
    @test size(H) == (m, n)
    @test nnz(H) == tnnz
    Hr = zeros(m, n)
    for j = 1:n
        Hr[:,j] = Array(A[j])
    end
    @test Array(H) == Hr

    V = vcat(A...)
    @test isa(V, SparseVector{Float64,Int})
    @test length(V) == m * n
    Vr = vec(Hr)
    @test Array(V) == Vr
end

# Test that concatenations of combinations of sparse vectors with various other
# matrix/vector types yield sparse arrays
let N = 4
    spvec = spzeros(N)
    spmat = spzeros(N, 1)
    densevec = ones(N)
    densemat = ones(N, 1)
    diagmat = Diagonal(ones(4))
    # Test that concatenations of pairwise combinations of sparse vectors with dense
    # vectors/matrices, sparse matrices, or special matrices yield sparse arrays
    for othervecormat in (densevec, densemat, spmat)
        @test issparse(vcat(spvec, othervecormat))
        @test issparse(vcat(othervecormat, spvec))
    end
    for othervecormat in (densevec, densemat, spmat, diagmat)
        @test issparse(hcat(spvec, othervecormat))
        @test issparse(hcat(othervecormat, spvec))
        @test issparse(hvcat((2,), spvec, othervecormat))
        @test issparse(hvcat((2,), othervecormat, spvec))
        @test issparse(cat((1,2), spvec, othervecormat))
        @test issparse(cat((1,2), othervecormat, spvec))
    end
    # The preceding tests should cover multi-way combinations of those types, but for good
    # measure test a few multi-way combinations involving those types
    @test issparse(vcat(spvec, densevec, spmat, densemat))
    @test issparse(vcat(densevec, spvec, densemat, spmat))
    @test issparse(hcat(spvec, densemat, spmat, densevec, diagmat))
    @test issparse(hcat(densemat, spmat, spvec, densevec, diagmat))
    @test issparse(hvcat((5,), diagmat, densevec, spvec, densemat, spmat))
    @test issparse(hvcat((5,), spvec, densemat, diagmat, densevec, spmat))
    @test issparse(cat((1,2), densemat, diagmat, spmat, densevec, spvec))
    @test issparse(cat((1,2), spvec, diagmat, densevec, spmat, densemat))
end

## sparsemat: combinations with sparse matrix

let S = sprand(4, 8, 0.5)
    Sf = Array(S)
    @assert isa(Sf, Matrix{Float64})

    # get a single column
    for j = 1:size(S,2)
        col = S[:, j]
        @test isa(col, SparseVector{Float64,Int})
        @test length(col) == size(S,1)
        @test Array(col) == Sf[:,j]
    end

    # Get a reshaped vector
    v = S[:]
    @test isa(v, SparseVector{Float64,Int})
    @test length(v) == length(S)
    @test Array(v) == Sf[:]

    # Get a linear subset
    for i=0:length(S)
        v = S[1:i]
        @test isa(v, SparseVector{Float64,Int})
        @test length(v) == i
        @test Array(v) == Sf[1:i]
    end
    for i=1:length(S)+1
        v = S[i:end]
        @test isa(v, SparseVector{Float64,Int})
        @test length(v) == length(S) - i + 1
        @test Array(v) == Sf[i:end]
    end
    for i=0:div(length(S),2)
        v = S[1+i:end-i]
        @test isa(v, SparseVector{Float64,Int})
        @test length(v) == length(S) - 2i
        @test Array(v) == Sf[1+i:end-i]
    end
end

let r = [1,10], S = sparse(r, r, r)
    Sf = Array(S)
    @assert isa(Sf, Matrix{Int})

    inds = [1,1,1,1,1,1]
    v = S[inds]
    @test isa(v, SparseVector{Int,Int})
    @test length(v) == length(inds)
    @test Array(v) == Sf[inds]

    inds = [2,2,2,2,2,2]
    v = S[inds]
    @test isa(v, SparseVector{Int,Int})
    @test length(v) == length(inds)
    @test Array(v) == Sf[inds]

    # get a single column
    for j = 1:size(S,2)
        col = S[:, j]
        @test isa(col, SparseVector{Int,Int})
        @test length(col) == size(S,1)
        @test Array(col) == Sf[:,j]
    end

    # Get a reshaped vector
    v = S[:]
    @test isa(v, SparseVector{Int,Int})
    @test length(v) == length(S)
    @test Array(v) == Sf[:]

    # Get a linear subset
    for i=0:length(S)
        v = S[1:i]
        @test isa(v, SparseVector{Int,Int})
        @test length(v) == i
        @test Array(v) == Sf[1:i]
    end
    for i=1:length(S)+1
        v = S[i:end]
        @test isa(v, SparseVector{Int,Int})
        @test length(v) == length(S) - i + 1
        @test Array(v) == Sf[i:end]
    end
    for i=0:div(length(S),2)
        v = S[1+i:end-i]
        @test isa(v, SparseVector{Int,Int})
        @test length(v) == length(S) - 2i
        @test Array(v) == Sf[1+i:end-i]
    end
end

## math

### Data

rnd_x0 = sprand(50, 0.6)
rnd_x0f = Array(rnd_x0)

rnd_x1 = sprand(50, 0.7) * 4.0
rnd_x1f = Array(rnd_x1)

spv_x1 = SparseVector(8, [2, 5, 6], [1.25, -0.75, 3.5])
spv_x2 = SparseVector(8, [1, 2, 6, 7], [3.25, 4.0, -5.5, -6.0])

### Arithmetic operations

let x = spv_x1, x2 = spv_x2
    # negate
    @test exact_equal(-x, SparseVector(8, [2, 5, 6], [-1.25, 0.75, -3.5]))

    # abs and abs2
    @test exact_equal(abs.(x), SparseVector(8, [2, 5, 6], abs.([1.25, -0.75, 3.5])))
    @test exact_equal(abs2.(x), SparseVector(8, [2, 5, 6], abs2.([1.25, -0.75, 3.5])))

    # plus and minus
    xa = SparseVector(8, [1,2,5,6,7], [3.25,5.25,-0.75,-2.0,-6.0])

    @test exact_equal(x + x, x * 2)
    @test exact_equal(x + x2, xa)
    @test exact_equal(x2 + x, xa)

    xb = SparseVector(8, [1,2,5,6,7], [-3.25,-2.75,-0.75,9.0,6.0])
    @test exact_equal(x - x, SparseVector(8, Int[], Float64[]))
    @test exact_equal(x - x2, xb)
    @test exact_equal(x2 - x, -xb)

    @test Array(x) + x2 == Array(xa)
    @test Array(x) - x2 == Array(xb)
    @test x + Array(x2) == Array(xa)
    @test x - Array(x2) == Array(xb)

    # multiplies
    xm = SparseVector(8, [2, 6], [5.0, -19.25])
    @test exact_equal(x .* x, abs2.(x))
    @test exact_equal(x .* x2, xm)
    @test exact_equal(x2 .* x, xm)

    @test Array(x) .* x2 == Array(xm)
    @test x .* Array(x2) == Array(xm)

    # max & min
    @test exact_equal(max.(x, x), x)
    @test exact_equal(min.(x, x), x)
    @test exact_equal(max.(x, x2),
        SparseVector(8, Int[1, 2, 6], Float64[3.25, 4.0, 3.5]))
    @test exact_equal(min.(x, x2),
        SparseVector(8, Int[2, 5, 6, 7], Float64[1.25, -0.75, -5.5, -6.0]))
end

### Complex

let x = spv_x1, x2 = spv_x2
    # complex
    @test exact_equal(complex.(x, x),
        SparseVector(8, [2,5,6], [1.25+1.25im, -0.75-0.75im, 3.5+3.5im]))
    @test exact_equal(complex.(x, x2),
        SparseVector(8, [1,2,5,6,7], [3.25im, 1.25+4.0im, -0.75+0.0im, 3.5-5.5im, -6.0im]))
    @test exact_equal(complex.(x2, x),
        SparseVector(8, [1,2,5,6,7], [3.25+0.0im, 4.0+1.25im, -0.75im, -5.5+3.5im, -6.0+0.0im]))

    # real & imag

    @test real(x) === x
    @test exact_equal(imag(x), spzeros(Float64, length(x)))

    xcp = complex.(x, x2)
    @test exact_equal(real(xcp), x)
    @test exact_equal(imag(xcp), x2)
end

### Zero-preserving math functions: sparse -> sparse

function check_nz2z_z2z{T}(f::Function, x::SparseVector{T}, xf::Vector{T})
    R = typeof(f(zero(T)))
    r = f(x)
    isa(r, AbstractSparseVector) || error("$f(x) is not a sparse vector.")
    eltype(r) == R || error("$f(x) results in eltype = $(eltype(r)), expect $R")
    all(r.nzval .!= 0) || error("$f(x) contains zeros in nzval.")
    Array(r) == f.(xf) || error("Incorrect results found in $f(x).")
end

for f in [floor, ceil, trunc, round]
    check_nz2z_z2z(f, rnd_x1, rnd_x1f)
end

for f in [log1p, expm1,
          sin, tan, sinpi, sind, tand,
          asin, atan, asind, atand,
          sinh, tanh, asinh, atanh]
    check_nz2z_z2z(f, rnd_x0, rnd_x0f)
end

### Non-zero-preserving math functions: sparse -> dense

function check_z2nz{T}(f::Function, x::SparseVector{T}, xf::Vector{T})
    R = typeof(f(zero(T)))
    r = f(x)
    isa(r, Vector) || error("$f(x) is not a dense vector.")
    eltype(r) == R || error("$f(x) results in eltype = $(eltype(r)), expect $R")
    r == f.(xf) || error("Incorrect results found in $f(x).")
end

for f in [exp, exp2, exp10, log, log2, log10,
          cos, csc, cot, sec, cospi,
          cosd, cscd, cotd, secd,
          acos, acot, acosd, acotd,
          cosh, csch, coth, sech, acsch, asech]
    check_z2nz(f, rnd_x0, rnd_x0f)
end

### Reduction

# sum, vecnorm

let x = spv_x1
    @test sum(x) == 4.0
    @test sum(abs, x) == 5.5
    @test sum(abs2, x) == 14.375

    @test vecnorm(x) == sqrt(14.375)
    @test vecnorm(x, 1) == 5.5
    @test vecnorm(x, 2) == sqrt(14.375)
    @test vecnorm(x, Inf) == 3.5
end

# maximum, minimum

let x = spv_x1
    @test maximum(x) == 3.5
    @test minimum(x) == -0.75
    @test maximum(abs, x) == 3.5
    @test minimum(abs, x) == 0.0
end

let x = abs.(spv_x1)
    @test maximum(x) == 3.5
    @test minimum(x) == 0.0
end

let x = -abs.(spv_x1)
    @test maximum(x) == 0.0
    @test minimum(x) == -3.5
end

let x = SparseVector(3, [1, 2, 3], [-4.5, 2.5, 3.5])
    @test maximum(x) == 3.5
    @test minimum(x) == -4.5
    @test maximum(abs, x) == 4.5
    @test minimum(abs, x) == 2.5
end

let x = spzeros(Float64, 8)
    @test maximum(x) == 0.0
    @test minimum(x) == 0.0
    @test maximum(abs, x) == 0.0
    @test minimum(abs, x) == 0.0
end

### linalg

### BLAS Level-1

let x = sprand(16, 0.5), x2 = sprand(16, 0.4)
    xf = Array(x)
    xf2 = Array(x2)

    # axpy!
    for c in [1.0, -1.0, 2.0, -2.0]
        y = Array(x)
        @test Base.axpy!(c, x2, y) === y
        @test y == Array(x2 * c + x)
    end

    # scale
    let α = 2.5, sx = SparseVector(x.n, x.nzind, x.nzval * α)
        @test exact_equal(x * α, sx)
        @test exact_equal(x * (α + 0.0*im), complex(sx))
        @test exact_equal(α * x, sx)
        @test exact_equal((α + 0.0*im) * x, complex(sx))
        @test exact_equal(x * α, sx)
        @test exact_equal(α * x, sx)
        @test exact_equal(x .* α, sx)
        @test exact_equal(α .* x, sx)
        @test exact_equal(x / α, SparseVector(x.n, x.nzind, x.nzval / α))

        xc = copy(x)
        @test scale!(xc, α) === xc
        @test exact_equal(xc, sx)
        xc = copy(x)
        @test scale!(α, xc) === xc
        @test exact_equal(xc, sx)
        xc = copy(x)
        @test scale!(xc, complex(α, 0.0)) === xc
        @test exact_equal(xc, sx)
        xc = copy(x)
        @test scale!(complex(α, 0.0), xc) === xc
        @test exact_equal(xc, sx)
    end

    # dot
    let dv = dot(xf, xf2)
        @test dot(x, x) == sum(abs2, x)
        @test dot(x2, x2) == sum(abs2, x2)
        @test dot(x, x2) ≈ dv
        @test dot(x2, x) ≈ dv
        @test dot(Array(x), x2) ≈ dv
        @test dot(x, Array(x2)) ≈ dv
    end
end

let x = complex.(sprand(32, 0.6), sprand(32, 0.6)),
    y = complex.(sprand(32, 0.6), sprand(32, 0.6))
    xf = Array(x)::Vector{Complex128}
    yf = Array(y)::Vector{Complex128}
    @test dot(x, x) ≈ dot(xf, xf)
    @test dot(x, y) ≈ dot(xf, yf)
end

### BLAS Level-2:

## dense A * sparse x -> dense y

let A = randn(9, 16), x = sprand(16, 0.7)
    xf = Array(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α*A*xf + β*y
        @test A_mul_B!(α, A, x, β, y) === y
        @test y ≈ rr
    end
    y = A*x
    @test isa(y, Vector{Float64})
    @test A*x ≈ A*xf
end

let A = randn(16, 9), x = sprand(16, 0.7)
    xf = Array(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α*A'xf + β*y
        @test At_mul_B!(α, A, x, β, y) === y
        @test y ≈ rr
    end
    y = At_mul_B(A, x)
    @test isa(y, Vector{Float64})
    @test y ≈ At_mul_B(A, xf)
end

## sparse A * sparse x -> dense y

let A = sprandn(9, 16, 0.5), x = sprand(16, 0.7)
    Af = Array(A)
    xf = Array(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α*Af*xf + β*y
        @test A_mul_B!(α, A, x, β, y) === y
        @test y ≈ rr
    end
    y = SparseArrays.densemv(A, x)
    @test isa(y, Vector{Float64})
    @test y ≈ Af*xf
end

let A = sprandn(16, 9, 0.5), x = sprand(16, 0.7)
    Af = Array(A)
    xf = Array(x)
    for α in [0.0, 1.0, 2.0], β in [0.0, 0.5, 1.0]
        y = rand(9)
        rr = α*Af'xf + β*y
        @test At_mul_B!(α, A, x, β, y) === y
        @test y ≈ rr
    end
    y = SparseArrays.densemv(A, x; trans='T')
    @test isa(y, Vector{Float64})
    @test y ≈ At_mul_B(Af, xf)
end

let A = complex.(sprandn(7, 8, 0.5), sprandn(7, 8, 0.5)),
    x = complex.(sprandn(8, 0.6), sprandn(8, 0.6)),
    x2 = complex.(sprandn(7, 0.75), sprandn(7, 0.75))
    Af = Array(A)
    xf = Array(x)
    x2f = Array(x2)
    @test SparseArrays.densemv(A, x; trans='N') ≈ Af * xf
    @test SparseArrays.densemv(A, x2; trans='T') ≈ Af.' * x2f
    @test SparseArrays.densemv(A, x2; trans='C') ≈ Af'x2f
    @test_throws ArgumentError SparseArrays.densemv(A, x; trans='D')
end

## sparse A * sparse x -> sparse y

let A = sprandn(9, 16, 0.5), x = sprand(16, 0.7), x2 = sprand(9, 0.7)
    Af = Array(A)
    xf = Array(x)
    x2f = Array(x2)

    y = A*x
    @test isa(y, SparseVector{Float64,Int})
    @test all(nonzeros(y) .!= 0.0)
    @test Array(y) ≈ Af * xf

    y = At_mul_B(A, x2)
    @test isa(y, SparseVector{Float64,Int})
    @test all(nonzeros(y) .!= 0.0)
    @test Array(y) ≈ Af'x2f
end

let A = complex.(sprandn(7, 8, 0.5), sprandn(7, 8, 0.5)),
    x = complex.(sprandn(8, 0.6), sprandn(8, 0.6)),
    x2 = complex.(sprandn(7, 0.75), sprandn(7, 0.75))
    Af = Array(A)
    xf = Array(x)
    x2f = Array(x2)

    y = A*x
    @test isa(y, SparseVector{Complex128,Int})
    @test Array(y) ≈ Af * xf

    y = At_mul_B(A, x2)
    @test isa(y, SparseVector{Complex128,Int})
    @test Array(y) ≈ Af.' * x2f

    y = Ac_mul_B(A, x2)
    @test isa(y, SparseVector{Complex128,Int})
    @test Array(y) ≈ Af'x2f
end

# left-division operations involving triangular matrices and sparse vectors (#14005)
let m = 10
    sparsefloatvecs = SparseVector[sprand(m, 0.4) for k in 1:3]
    sparseintvecs = SparseVector[SparseVector(m, sprvec.nzind, round.(Int, sprvec.nzval*10)) for sprvec in sparsefloatvecs]
    sparsecomplexvecs = SparseVector[SparseVector(m, sprvec.nzind, complex.(sprvec.nzval, sprvec.nzval)) for sprvec in sparsefloatvecs]

    sprmat = sprand(m, m, 0.2)
    sparsefloatmat = speye(m) + sprmat/(2m)
    sparsecomplexmat = speye(m) + SparseMatrixCSC(m, m, sprmat.colptr, sprmat.rowval, complex.(sprmat.nzval, sprmat.nzval)/(4m))
    sparseintmat = speye(Int, m)*10m + SparseMatrixCSC(m, m, sprmat.colptr, sprmat.rowval, round.(Int, sprmat.nzval*10))

    denseintmat = eye(Int, m)*10m + rand(1:m, m, m)
    densefloatmat = eye(m) + randn(m, m)/(2m)
    densecomplexmat = eye(m) + complex.(randn(m, m), randn(m, m))/(4m)

    inttypes = (Int32, Int64, BigInt)
    floattypes = (Float32, Float64, BigFloat)
    complextypes = (Complex{Float32}, Complex{Float64})
    eltypes = (inttypes..., floattypes..., complextypes...)

    for eltypemat in eltypes
        (densemat, sparsemat) = eltypemat in inttypes ? (denseintmat, sparseintmat) :
                                eltypemat in floattypes ? (densefloatmat, sparsefloatmat) :
                                eltypemat in complextypes && (densecomplexmat, sparsecomplexmat)
        densemat = convert(Matrix{eltypemat}, densemat)
        sparsemat = convert(SparseMatrixCSC{eltypemat}, sparsemat)
        trimats = (LowerTriangular(densemat), UpperTriangular(densemat),
                   LowerTriangular(sparsemat), UpperTriangular(sparsemat) )
        unittrimats = (Base.LinAlg.UnitLowerTriangular(densemat), Base.LinAlg.UnitUpperTriangular(densemat),
                       Base.LinAlg.UnitLowerTriangular(sparsemat), Base.LinAlg.UnitUpperTriangular(sparsemat) )

        for eltypevec in eltypes
            spvecs = eltypevec in inttypes ? sparseintvecs :
                     eltypevec in floattypes ? sparsefloatvecs :
                     eltypevec in complextypes && sparsecomplexvecs
            spvecs = SparseVector[SparseVector(m, spvec.nzind, convert(Vector{eltypevec}, spvec.nzval)) for spvec in spvecs]

            for spvec in spvecs
                fspvec = convert(Array, spvec)
                # test out-of-place left-division methods
                for mat in (trimats..., unittrimats...), func in (\, At_ldiv_B, Ac_ldiv_B)
                    @test func(mat, spvec) ≈ func(mat, fspvec)
                end
                # test in-place left-division methods not involving quotients
                if eltypevec == typeof(zero(eltypemat)*zero(eltypevec) + zero(eltypemat)*zero(eltypevec))
                    for mat in unittrimats, func in (A_ldiv_B!, Base.LinAlg.At_ldiv_B!, Base.LinAlg.Ac_ldiv_B!)
                        @test func(mat, copy(spvec)) ≈ func(mat, copy(fspvec))
                    end
                end
                # test in-place left-division methods involving quotients
                if eltypevec == typeof((zero(eltypemat)*zero(eltypevec) + zero(eltypemat)*zero(eltypevec))/one(eltypemat))
                    for mat in trimats, func in (A_ldiv_B!, Base.LinAlg.At_ldiv_B!, Base.LinAlg.Ac_ldiv_B!)
                        @test func(mat, copy(spvec)) ≈ func(mat, copy(fspvec))
                    end
                end
            end
        end
    end
end
# The preceding tests miss the edge case where the sparse vector is empty (#16716)
let origmat = [-1.5 -0.7; 0.0 1.0]
    transmat = transpose(origmat)
    utmat = UpperTriangular(origmat)
    ltmat = LowerTriangular(transmat)
    uutmat = Base.LinAlg.UnitUpperTriangular(origmat)
    ultmat = Base.LinAlg.UnitLowerTriangular(transmat)

    zerospvec = spzeros(Float64, 2)
    zerodvec = zeros(Float64, 2)

    for mat in (utmat, ltmat, uutmat, ultmat)
        for func in (\, At_ldiv_B, Ac_ldiv_B)
            @test isequal((func)(mat, zerospvec), zerodvec)
        end
        for ipfunc in (A_ldiv_B!, Base.LinAlg.At_ldiv_B!, Base.LinAlg.Ac_ldiv_B!)
            @test isequal((ipfunc)(mat, copy(zerospvec)), zerospvec)
        end
    end
end

# fkeep!
let x = sparsevec(1:7, [3., 2., -1., 1., -2., -3., 3.], 7)
    # droptol
    xdrop = Base.droptol!(copy(x), 1.5)
    @test exact_equal(xdrop, SparseVector(7, [1, 2, 5, 6, 7], [3., 2., -2., -3., 3.]))
    Base.droptol!(xdrop, 2.5)
    @test exact_equal(xdrop, SparseVector(7, [1, 6, 7], [3., -3., 3.]))
    Base.droptol!(xdrop, 3.)
    @test exact_equal(xdrop, SparseVector(7, Int[], Float64[]))

    xdrop = copy(x)
    # This will keep index 1, 3, 4, 7 in xdrop
    f_drop(i, x) = (abs(x) == 1.) || (i in [1, 7])
    Base.SparseArrays.fkeep!(xdrop, f_drop)
    @test exact_equal(xdrop, SparseVector(7, [1, 3, 4, 7], [3., -1., 1., 3.]))
end

# dropzeros[!]
let testdims = (10, 20, 30), nzprob = 0.4, targetnumposzeros = 5, targetnumnegzeros = 5
    for m in testdims
        v = sprand(m, nzprob)
        struczerosv = find(x -> x == 0, v)
        poszerosinds = unique(rand(struczerosv, targetnumposzeros))
        negzerosinds = unique(rand(struczerosv, targetnumnegzeros))
        vposzeros = setindex!(copy(v), 2, poszerosinds)
        vnegzeros = setindex!(copy(v), -2, negzerosinds)
        vbothsigns = setindex!(copy(vposzeros), -2, negzerosinds)
        map!(x -> x == 2 ? 0.0 : x, vposzeros.nzval, vposzeros.nzval)
        map!(x -> x == -2 ? -0.0 : x, vnegzeros.nzval, vnegzeros.nzval)
        map!(x -> x == 2 ? 0.0 : x == -2 ? -0.0 : x, vbothsigns.nzval, vbothsigns.nzval)
        for vwithzeros in (vposzeros, vnegzeros, vbothsigns)
            # Basic functionality / dropzeros!
            @test dropzeros!(copy(vwithzeros)) == v
            @test dropzeros!(copy(vwithzeros), false) == v
            # Basic functionality / dropzeros
            @test dropzeros(vwithzeros) == v
            @test dropzeros(vwithzeros, false) == v
            # Check trimming works as expected
            @test length(dropzeros!(copy(vwithzeros)).nzval) == length(v.nzval)
            @test length(dropzeros!(copy(vwithzeros)).nzind) == length(v.nzind)
            @test length(dropzeros!(copy(vwithzeros), false).nzval) == length(vwithzeros.nzval)
            @test length(dropzeros!(copy(vwithzeros), false).nzind) == length(vwithzeros.nzind)
        end
    end
    # original dropzeros! test
    xdrop = sparsevec(1:7, [3., 2., -1., 1., -2., -3., 3.], 7)
    xdrop.nzval[[2, 4, 6]] = 0.0
    Base.SparseArrays.dropzeros!(xdrop)
    @test exact_equal(xdrop, SparseVector(7, [1, 3, 5, 7], [3, -1., -2., 3.]))
end

# It's tempting to share data between a SparseVector and a SparseMatrix,
# but if that's done, then modifications to one or the other will cause
# an inconsistent state:
sv = sparse(1:10)
sm = convert(SparseMatrixCSC, sv)
sv[1] = 0
@test Array(sm)[2:end] == collect(2:10)

# Ensure that sparsevec with all-zero values returns an array of zeros
@test sparsevec([1,2,3],[0,0,0]) == [0,0,0]

# Compare stored zero semantics between SparseVector and SparseMatrixCSC
let S = SparseMatrixCSC(10,1,[1,6],[1,3,5,6,7],[0,1,2,0,3]), x = SparseVector(10,[1,3,5,6,7],[0,1,2,0,3])
    @test nnz(S) == nnz(x) == 5
    for I = (:, 1:10, collect(1:10))
        @test S[I,1] == S[I] == x[I] == x
        @test nnz(S[I,1]) == nnz(S[I]) == nnz(x[I]) == nnz(x)
    end
    for I = (2:9, 1:2, 9:10, [3,6,1], [10,9,8], [])
        @test S[I,1] == S[I] == x[I]
        @test nnz(S[I,1]) == nnz(S[I]) == nnz(x[I])
    end
    @test S[[1 3 5; 2 4 6]] == x[[1 3 5; 2 4 6]]
    @test nnz(S[[1 3 5; 2 4 6]]) == nnz(x[[1 3 5; 2 4 6]])
end

# Issue 14013
s14013 = sparse([10.0 0.0 30.0; 0.0 1.0 0.0])
a14013 = [10.0 0.0 30.0; 0.0 1.0 0.0]
@test s14013 == a14013
@test vec(s14013) == s14013[:] == a14013[:]
@test Array(s14013)[1,:] == s14013[1,:] == a14013[1,:] == [10.0, 0.0, 30.0]
@test Array(s14013)[2,:] == s14013[2,:] == a14013[2,:] == [0.0, 1.0, 0.0]

# Issue 14046
s14046 = sprand(5, 1.0)
@test spzeros(5) + s14046 == s14046
@test 2*s14046 == s14046 + s14046

# Issue 14589
#test vectors with no zero elements
x = sparsevec(1:7, [3., 2., -1., 1., -2., -3., 3.], 7)
@test collect(sort(x)) == sort(collect(x))
#test vectors with all zero elements
x = sparsevec(Int64[], Float64[], 7)
@test collect(sort(x)) == sort(collect(x))
#test vector with sparsity approx 1/2
x = sparsevec(1:7, [3., 2., -1., 1., -2., -3., 3.], 15)
@test collect(sort(x)) == sort(collect(x))
# apply three distinct tranformations where zeros sort into start/middle/end
@test collect(sort(x, by=abs)) == sort(collect(x), by=abs)
@test collect(sort(x, by=sign)) == sort(collect(x), by=sign)
@test collect(sort(x, by=inv)) == sort(collect(x), by=inv)

#fill!
for Tv in [Float32, Float64, Int64, Int32, Complex128]
    for Ti in [Int16, Int32, Int64, BigInt]
        sptypes = (SparseMatrixCSC{Tv, Ti}, SparseVector{Tv, Ti})
        sizes = [(3, 4), (3,)]
        for (siz, Sp) in zip(sizes, sptypes)
            arr = rand(Tv, siz...)
            sparr = Sp(arr)
            fillval = rand(Tv)
            fill!(sparr, fillval)
            @test Array(sparr) == fillval * ones(arr)
            fill!(sparr, 0)
            @test Array(sparr) == zeros(arr)
        end
    end
end

# ref 13130 and 16661
@test issparse([sprand(10,10,.1) sprand(10,.1)])
@test issparse([sprand(10,1,.1); sprand(10,.1)])

@test issparse([sprand(10,10,.1) rand(10)])
@test issparse([sprand(10,1,.1)  rand(10)])
@test issparse([sprand(10,2,.1) sprand(10,1,.1) rand(10)])
@test issparse([sprand(10,1,.1); rand(10)])

@test issparse([sprand(10,.1)  rand(10)])
@test issparse([sprand(10,.1); rand(10)])


text:
__precompile__()

module Utils

    # export @define, @reexport
    #
    # include("utils/macro_utils.jl")

    export istriustrict, istrilstrict
    export simd_scale!, simd_copy!, simd_copy_scale!,
           simd_copy_xy_first!, simd_copy_yx_first!,
           simd_copy_yx_first_last!,
           simd_xpy!, simd_axpy!, simd_wxpy!, simd_waxpy!, simd_mult!

    include("utils/matrix_utils.jl")

end

ext: jl


text:
using Documenter, RasterDataSources

makedocs(
    sitename = "RasterDataSources.jl",
)

deploydocs(
    repo = "github.com/cesaraustralia/RasterDataSources.jl.git",
)


text:
using MAT, JLD2

@inline function string_as_varname(s::AbstractString)
    S = Symbol(s)
    @eval $S
end

#@inline function test_against_matlab(filename::String)
    vars = matread(filename)
    for v in vars
        cur_type = typeof(v[2])
        if (v[1] == "t0")
            continue
        elseif v[1] == "a_grid"
            @assert v[2] ≈ a_grid
        elseif v[1] == "a_g_grid"
            @assert v[2] ≈ a_g_grid
        elseif v[1] == "b_grid"
            @assert v[2] ≈ b_grid
        elseif v[1] == "y_grid"
            @assert v[2] ≈ y_grid
        elseif v[1] == "y_g_grid"
            @assert v[2] ≈ y_g_grid
        elseif v[1] == "r_a_grid"
            @assert v[2] ≈ r_a_grid
        elseif v[1] == "r_b_grid"
            @assert v[2] ≈ r_b_grid
        elseif v[1] == "r_a_g_grid"
            @assert v[2] ≈ r_a_g_grid
        elseif v[1] == "r_b_g_grid"
            @assert v[2] ≈ r_b_g_grid
        elseif v[1] == "daf_g_grid"
            @assert v[2] ≈ daf_g_grid
        elseif v[1] == "daf_grid"
            @assert v[2] ≈ daf_grid
        elseif v[1] == "dab_grid"
            @assert v[2] ≈ dab_grid
        elseif v[1] == "dab_g_grid"
            @assert v[2] ≈ dab_g_grid
        elseif v[1] == "dab_g_tilde_grid"
            @assert v[2] ≈ dab_g_tilde_grid
        elseif v[1] == "dbf_grid"
            @assert v[2] ≈ dbf_grid
        elseif v[1] == "dbf_g_grid"
            @assert v[2] ≈ dbf_g_grid
        elseif v[1] == "dbb_grid"
            @assert v[2] ≈ dbb_grid
        elseif v[1] == "dbb_g_grid"
            @assert v[2] ≈ dbb_g_grid
        elseif v[1] == "trans_grid"
            @assert v[2] ≈ trans_grid
        elseif v[1] == "l_g_grid"
            @assert v[2] ≈ l_g_grid
        elseif v[1] == "w_grid"
            @assert v[2] ≈ w_grid
        elseif v[1] == "c_0"
            @assert v[2] ≈ c_0
        elseif v[1] == "V_0"
            @assert v[2] ≈ V_0
        elseif v[1] == "gg0"
            @assert v[2] ≈ gg0
        elseif v[1] == "gg"
            @assert v[2] ≈ gg
        elseif v[1] == "r_b"
            @assert v[2] ≈ r_b
        elseif v[1] == "r_b_borr"
            @assert v[2] ≈ r_b_borr
        elseif v[1] == "y_dist"
            @assert v[2] ≈ y_dist
        elseif v[1] == "y_mean"
            @assert v[2] ≈ y_mean
        elseif v[1] == "aalpha"
            @assert v[2] ≈ aalpha
        elseif v[1] == "ddeath"
            @assert v[2] ≈ ddeath
        elseif v[1] == "ddelta"
            @assert v[2] ≈ ddelta
        elseif v[1] == "rrho"
            @assert v[2] ≈ rrho
        elseif v[1] == "chi0"
            @assert v[2] ≈ chi0
        elseif v[1] == "chi1"
            @assert v[2] ≈ chi1
        elseif v[1] == "chi2"
            @assert v[2] ≈ chi2
        elseif v[1] == "a_lb"
            @assert v[2] ≈ a_lb
        elseif v[1] == "pam"
            @assert v[2] ≈ pam
        elseif v[1] == "xxi"
            @assert v[2] ≈ xxi
        elseif v[1] == "ggamma"
            @assert v[2] ≈ ggamma
        elseif v[1] == "tau_I"
            @assert v[2] ≈ tau_I
        elseif v[1] == "trans"
            @assert v[2] ≈ trans
        elseif v[1] == "lambda"
            @assert v[2] ≈ lambda
        elseif v[1] == "K_liquid"
            @assert v[2] ≈ K_liquid
        elseif v[1] == "aggregate_variables"
            @assert v[2] ≈ aggregate_variables
        elseif v[1] == "distributional_variables"
            @assert v[2] ≈ distributional_variables
        elseif v[1] == "distributional_variables_1"
            #@assert v[2] ≈ distributional_variables_1
            continue
        else
            @show "NOT PRINTING", v[1]
            #=
            try
                if !(string_as_varname(v[1]) ≈ v[2])
                    @show "NEQ", v[2]
                else
                    @show "WORKSWORKSWORKS", v[2]
                end
                @assert string_as_varname(v[1]) ≈ v[2]
            catch err
                if isa(err, UndefVarError)
                    println(err)
                else
                    throw(err)
                end
            end
            =#
        end
    end
    println("All tests passed for: ", filename)
#end


text:
using BenchmarkTools

versioninfo(verbose=true)

# ## HWloc
using Pkg
Pkg.add("Hwloc")

using Hwloc
Hwloc.num_physical_cores()

# ##Starting Threads

using Base.Threads
nthreads()

threadid()

a=zeros(4)
for i in 1:nthreads()
     a[threadid()] = threadid()
end
a

# ## The @Threads macro
    @threads for i in 1:nthreads()
               a[threadid()] = threadid()
           end

# ### Prefix sum

function prefix_threads!(y)
    l=length(y)
    k=ceil(Int, log2(l))
    for j=1:k
        @threads for i=2^j:2^j:min(l, 2^k)
            @inbounds y[i] = y[i-2^(j-1)] + y[i]
        end
    end
    for j=(k-1):-1:1
        @threads for i=3*2^(j-1):2^j:min(l, 2^k)
            @inbounds y[i] = y[i-2^(j-1)] + y[i]
        end
    end
    y
end

# # Thread Safety

# ## Monte Carlo
using Random
function darts_in_circle(n, rng=Random.GLOBAL_RNG)
   inside = 0
   for i in 1:n
      if rand(rng)^2 + rand(rng)^2 < 1
         inside += 1
      end
   end
 return inside
end

function pi_serial(n)
    return 4 * darts_in_circle(n) / n
end

using Base.Threads

const rnglist = [MersenneTwister() for i in 1:nthreads()];

function pi_threads(n, loops)
   inside = zeros(Int, loops)
   @threads for i in 1:loops
      rng = rnglist[threadid()]
      inside[threadid()] = darts_in_circle(n, rng)
   end
   return 4 * sum(inside) / (n*loops)
end

pi_threads(2_500_000, 4)

 @btime pi_serial(10_000_000)

 @btime pi_threads(2_500_000, 4)

 # ## Atomics

 function sum_thread_base(x)
   r = zero(eltype(x))
   @threads for i in eachindex(x)
      @inbounds r += x[i]
   end
   return r
end

a=rand(10_000_000);

@btime sum($a)

@btime sum_thread_base($a)

function sum_thread_atomic(x)
   r = Atomic{eltype(x)}(zero(eltype(x)))
   @threads for i in eachindex(x)
      @inbounds atomic_add!(r, x[i])
   end
  return r[]
end

 @btime sum_thread_atomic($a)

 function sum_thread_split(A)
   r = Atomic{eltype(A)}(zero(eltype(A)))
   len, rem = divrem(length(A), nthreads())
   #Split the array equally among the threads
   @threads for t in 1:nthreads()
      r[] = zero(eltype(A))
      @simd for i in (1:len) .+ (t-1)*len
         @inbounds r[] += A[i]
      end
      atomic_add!(r, r[])
    end
   result = r[]
   #process up the remaining data
   @simd for i in length(A)-rem+1:length(A)
      @inbounds result += A[i]
   end
   return result
end

@btime sum_thread_split($a)

@btime sum_thread_split($a)

# ## Synchronisation primitives

const f = open(tempname(), "a+")

const mt = Base.Threads.Mutex();

@threads for i in 1:50
     r = pi_serial(10_000_000)
     lock(mt)
     write(f, "From $(threadid()), pi = $r\n")
     unlock(mt)
end

close(f)

const s = Base.Semaphore(2);
@threads for i in 1:nthreads()
   Base.acquire(s)
   r = pi_serial(10_000_000)
   Core.println("Calculated pi = $r in Thread $(threadid())")
   Base.release(s)
end

# # Threaded libraries

a = rand(1000, 1000);

b = rand(1000, 1000);

@btime $a*$b;

# ## Oversubscriptions

function matmul_serial(x)
   first_num = zeros(length(x))
   for i in eachindex(x)
      @inbounds first_num[i] = (x[i]'*x[i])[1]
   end
   return first_num
end

function matmul_thread(x)
   first_num = zeros(length(x))
   @threads for i in eachindex(x)
      @inbounds first_num[i] = (x[i]'*x[i])[1]
   end
   return first_num
end

m = [rand(1000, 1000) for _ in 1:10];

@btime matmul_serial(m);

@btime matmul_thread(m);

using LinearAlgebra
BLAS.set_num_threads(1)

@btime matmul_thread(m);


text:
using OrderedCollections: OrderedDict

function ordered_dict_show(dict, padchar)
    if length(dict) == 0
        return string(typeof(dict), " with no entries")
    elseif length(dict) == 1
        return string(typeof(dict), " with 1 entry:", '\n',
                      padchar, "   └── ", dict.keys[1], " => ", typeof(dict.vals[1]))
    else
        return string(typeof(dict), " with $(length(dict)) entries:", '\n',
                      Tuple(string(padchar,
                                   "   ├── ", name, " => ", typeof(dict[name]), '\n')
                            for name in dict.keys[1:end-1]
                           )...,
                           padchar, "   └── ", dict.keys[end], " => ", typeof(dict.vals[end])
                      )
    end
end


text:
"""
    struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D}
      # Private Fields
    end
"""
struct CartesianDiscreteModel{D,T,F} <: DiscreteModel{D,D}
  grid::CartesianGrid{D,T,F}
  grid_topology::UnstructuredGridTopology{D,D,T,Oriented}
  face_labeling::FaceLabeling
  @doc """
      CartesianDiscreteModel(desc::CartesianDescriptor)

  Inner constructor
  """
  function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F}) where {D,T,F}
    grid = CartesianGrid(desc)
    _grid = UnstructuredGrid(grid)
    if any(desc.isperiodic)
      topo = _cartesian_grid_topology_with_periodic_bcs(_grid, desc.isperiodic, desc.partition)
    else
      topo = UnstructuredGridTopology(_grid)
    end
    nfaces = [num_faces(topo,d) for d in 0:num_cell_dims(topo)]
    labels = FaceLabeling(nfaces)
    _fill_cartesian_face_labeling!(labels,topo)
    new{D,T,F}(grid,topo,labels)
  end

  @doc """
      CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F},
                             cmin::CartesianIndex,
                             cmax::CartesianIndex)

      Builds a CartesianDiscreteModel object which represents a subgrid of
      a (larger) grid represented by desc. This subgrid is described by its
      D-dimensional minimum (cmin) and maximum (cmax) CartesianIndex
      identifiers.

  Inner constructor
  """
  function CartesianDiscreteModel(desc::CartesianDescriptor{D,T,F},
                                  cmin::CartesianIndex,
                                  cmax::CartesianIndex) where {D,T,F}

     suborigin = Tuple(desc.origin) .+ (Tuple(cmin) .- 1) .* desc.sizes
     subpartition = Tuple(cmax) .- Tuple(cmin) .+ 1
     subsizes = desc.sizes
     subdesc =
       CartesianDescriptor(Point(suborigin), subsizes, subpartition; map=desc.map, isperiodic=desc.isperiodic)

     grid = CartesianGrid(subdesc)
     _grid = UnstructuredGrid(grid)
     topo = UnstructuredGridTopology(_grid)
     nfaces = [num_faces(topo, d) for d = 0:num_cell_dims(topo)]
     labels = FaceLabeling(nfaces)
     _fill_subgrid_cartesian_face_labeling!(labels,topo,subdesc,desc,cmin)
     new{D,T,F}(grid, topo, labels)
  end
end

"""
    CartesianDiscreteModel(args...)

Same args needed to construct a `CartesianDescriptor`
"""
function CartesianDiscreteModel(args...; kwargs...)
  desc = CartesianDescriptor(args...; kwargs...)
  CartesianDiscreteModel(desc)
end

"""
    get_cartesian_descriptor(model::CartesianDiscreteModel)
"""
function get_cartesian_descriptor(model::CartesianDiscreteModel)
  get_cartesian_descriptor(model.grid)
end

# Interface

get_grid(model::CartesianDiscreteModel) = model.grid

get_grid_topology(model::CartesianDiscreteModel) = model.grid_topology

get_face_labeling(model::CartesianDiscreteModel) = model.face_labeling

# These needed to be type stable

function get_face_nodes(model::CartesianDiscreteModel,d::Integer)
  face_nodes::Table{Int32,Vector{Int32},Vector{Int32}} = compute_face_nodes(model,d)
  face_nodes
end

function get_face_type(model::CartesianDiscreteModel,d::Integer)
  _, face_to_ftype::Vector{Int8} = compute_reffaces(ReferenceFE{d},model)
  face_to_ftype
end

function get_reffaces(::Type{ReferenceFE{d}},model::CartesianDiscreteModel) where d
  reffaces::Vector{LagrangianRefFE{d}},_ = compute_reffaces(ReferenceFE{d},model)
  reffaces
end

# Helpers

function _fill_cartesian_face_labeling!(labels,topo)
  _fill_cartesian_entities!(labels,topo)
  _add_cartesian_tags!(labels,topo)
end

function _fill_cartesian_entities!(labels,topo)
  D = num_cell_dims(topo)
  d_to_dface_to_entity = labels.d_to_dface_to_entity
  polytope = first(get_polytopes(topo))
  dim_to_offset = get_offsets(polytope)
  interior_id = num_faces(polytope)
  boundary_id = -1
  for d in 0:(D-1)
    face_to_cells = get_faces(topo,d,D)
    cell_to_faces = get_faces(topo,D,d)
    offset = dim_to_offset[d+1]
    _generate_pre_geolabel!(
      d_to_dface_to_entity[d+1],
      face_to_cells,
      cell_to_faces,
      offset,
      interior_id,boundary_id,d,D)
  end
  _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id)
  fill!(d_to_dface_to_entity[end],interior_id)
end

function _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id)
  for d = 0:(D-2)
    for j = (d+1):(D-1)
      dface_to_jfaces = get_faces(topo, d, j)
      dface_to_geolabel = d_to_dface_to_entity[d+1]
      jface_to_geolabel = d_to_dface_to_entity[j+1]
      _fix_dface_geolabels!(
        dface_to_geolabel,
        jface_to_geolabel,
        dface_to_jfaces.data,
        dface_to_jfaces.ptrs,
        interior_id,
        boundary_id,
      )
    end
  end
end

function _add_cartesian_tags!(labels,topo)
  D = num_cell_dims(topo)
  polytope = first(get_polytopes(topo))
  interior_id = num_faces(polytope)
  boundary_ids = collect(1:(interior_id-1))
  for i in boundary_ids
    name = lpad(i,ceil(Int,log10(interior_id)),'0')
    add_tag!(labels,"tag_$(name)",[i])
  end
  add_tag!(labels,"interior",[interior_id])
  add_tag!(labels,"boundary",boundary_ids)
end

function _generate_pre_geolabel!(
  face_to_geolabel,
  face_to_cells,
  cell_to_faces,
  offset,
  interior_id,
  boundary_id,d,D)

  nfaces = length(face_to_cells)
  fill!(face_to_geolabel,interior_id)

  max_ncells_around = 2^(D-d)

  _generate_pre_geolabel_kernel!(
    face_to_geolabel,
    face_to_cells.data,
    face_to_cells.ptrs,
    cell_to_faces.data,
    cell_to_faces.ptrs,
    offset,boundary_id,max_ncells_around)

  face_to_geolabel
end

function _generate_pre_geolabel_kernel!(
  face_to_geolabel,
  face_to_cells_data,
  face_to_cells_ptrs,
  cell_to_faces_data,
  cell_to_faces_ptrs,
  offset,boundary_id,max_ncells_around)
  nfaces = length(face_to_geolabel)
  for face in 1:nfaces
    a = face_to_cells_ptrs[face]-1
    ncells_around = face_to_cells_ptrs[face+1] - (a+1)
    if ncells_around == 1
      icell_around = 1
      cell = face_to_cells_data[a+icell_around]
      b = cell_to_faces_ptrs[cell]-1
      nlfaces = cell_to_faces_ptrs[cell+1] - (b+1)
      for lface in 1:nlfaces
        face2 = cell_to_faces_data[b+lface]
        if face == face2
          face_to_geolabel[face] = lface + offset
          break
        end
      end
    elseif ncells_around != max_ncells_around
      face_to_geolabel[face] = boundary_id
    end
  end
end

function _fix_dface_geolabels!(
  dface_to_geolabel,
  jface_to_geolabel,
  dface_to_jfaces_data,
  dface_to_jfaces_ptrs,
  interior_id,boundary_id)

  ndfaces = length(dface_to_jfaces_ptrs)-1
  for dface in 1:ndfaces
    if dface_to_geolabel[dface] != boundary_id
      continue
    end
    a = dface_to_jfaces_ptrs[dface]
    b = dface_to_jfaces_ptrs[dface+1]-1
    for p in a:b
      jface = dface_to_jfaces_data[p]
      geolabel = jface_to_geolabel[jface]
      if geolabel != interior_id && geolabel != boundary_id
        dface_to_geolabel[dface] = geolabel
        break
      end
    end
  end
end

function _fill_subgrid_cartesian_face_labeling!(labels,topo,subdesc,desc,cmin)
   _fill_subgrid_cartesian_entities!(labels,topo,subdesc,desc,cmin)
   _add_cartesian_tags!(labels,topo)
end

function _fill_subgrid_cartesian_entities!(labels, topo, subdesc, desc, cmin)
  D = num_cell_dims(topo)
  d_to_dface_to_entity = labels.d_to_dface_to_entity
  gcis = CartesianIndices(desc.partition)
  subcis = CartesianIndices(subdesc.partition)

  polytope = first(get_polytopes(topo))
  face_labeling = labels
  offsets = get_offsets(polytope)
  interior_id = num_faces(polytope)
  boundary_id = -1

  minus_one_ci = CartesianIndex(tfill(-1, Val{D}()))

  polytope_d_face_to_jfaces = Matrix{Vector{Vector{Int}}}(undef, (D, D))
  for d = 0:(D-1)
    for j = d+1:D-1
      polytope_d_face_to_jfaces[d+1, j+1] = get_faces(polytope, d, j)
    end
  end

  face_deltas = _find_ncube_face_neighbor_deltas(polytope)
  for d = 0:(D-1)
    face_to_cells = get_faces(topo, d, D)
    cell_to_faces = get_faces(topo, D, d)
    face_to_geolabel = face_labeling.d_to_dface_to_entity[d+1]
    _generate_subgrid_pregeo_label!(
      face_to_geolabel,
      d,
      D,
      offsets,
      polytope_d_face_to_jfaces,
      minus_one_ci,
      face_to_cells,
      cell_to_faces,
      subcis,
      gcis,
      face_deltas,
      cmin,
      interior_id,
      boundary_id,
    )
  end
  _fix_geolabels(D, topo, d_to_dface_to_entity, interior_id, boundary_id)
  fill!(d_to_dface_to_entity[end], interior_id)
end

function _generate_subgrid_pregeo_label!(
  face_to_geolabel,
  d,
  D,
  offsets,
  polytope_d_face_to_jfaces,
  minus_one_ci,
  face_to_cells,
  cell_to_faces,
  subcis,
  gcis,
  face_deltas,
  cmin,
  interior_id,
  boundary_id,
)
  for face_gid = 1:length(face_to_geolabel)
    cell_gid = face_to_cells.data[face_to_cells.ptrs[face_gid]]
    a = cell_to_faces.ptrs[cell_gid]
    b = cell_to_faces.ptrs[cell_gid+1] - 1
    gci = (subcis[cell_gid] + minus_one_ci) + cmin

    face_lid = -1
    for j = a:b
      if (cell_to_faces.data[j] == face_gid)
        face_lid = j - a + 1
        break
      end
    end
    @assert face_lid != -1

    face_lid += offsets[d+1]
    # Check whether cell neighbour across face face_lid belongs to the
    # global grid. If yes, the current face is actually at the interior
    if ((gci + face_deltas[face_lid]) in gcis)
      face_to_geolabel[face_gid] = interior_id
    else
      cell_found = false
      for j = d+1:D-1
        dface_to_jfaces =
          polytope_d_face_to_jfaces[d+1, j+1][face_lid-offsets[d+1]]
        cell_found = _is_there_interior_cell_across_higher_dim_faces(
          dface_to_jfaces,
          offsets[j+1],
          gcis,
          gci,
          face_deltas,
        )
        cell_found && break
      end
      if (cell_found)
        # The current face is at least in two subgrid cells
        face_to_geolabel[face_gid] = boundary_id
      else
        # The current face is only in one subgrid cell
        face_to_geolabel[face_gid] = face_lid
      end
    end
  end
end

function _is_there_interior_cell_across_higher_dim_faces(
  dface_to_jfaces,
  offset_j,
  gcis,
  gci,
  face_deltas,
)
  for k in dface_to_jfaces
    jface_lid = k + offset_j
    if (isassigned(face_deltas, jface_lid))
      if ((gci + face_deltas[jface_lid]) in gcis)
        return true
      end
    end
  end
  return false
end

"""
  _find_ncube_face_neighbor_deltas(p::ExtrusionPolytope{D}) -> Vector{CartesianIndex}

  Given an n-cube type ExtrusionPolytope{D}, returns V=Vector{CartesianIndex} with as many
  entries as the number of faces in the boundary of the Polytope. For an entry face_lid
  in this vector, V[face_lid] returns what has to be added to the CartesianIndex of a
  cell in order to obtain the CartesianIndex of the cell neighbour of K across the face F
  with local ID face_lid.
"""
function _find_ncube_face_neighbor_deltas(p::ExtrusionPolytope{D}) where {D}
  nfaces = num_faces(p)
  delta_faces = Vector{CartesianIndex{D}}(undef, nfaces - 1)
  for face_lid = 1:nfaces-1
    delta_faces[face_lid] = _find_ncube_face_neighbor_delta(p, face_lid)
  end
  delta_faces
end

function _find_ncube_face_neighbor_delta(p::ExtrusionPolytope{D}, face_lid) where {D}
  @assert is_n_cube(p)
  result = fill(0, D)
  face = p.dface.nfaces[face_lid]
  for d = 1:D
    if (face.extrusion[d] == 0)
      result[d] = (face.anchor[d] == 0) ? -1 : 1
    end
  end
  return CartesianIndex(Tuple(result))
end

# Cartesian grid topology with periodic BC

function _cartesian_grid_topology_with_periodic_bcs(grid::UnstructuredGrid,
  isperiodic::NTuple,
  partition)

  cell_to_vertices, vertex_to_node =
    _generate_cell_to_vertices_from_grid(grid, isperiodic, partition)
  _generate_grid_topology_from_grid(grid,cell_to_vertices,vertex_to_node)
end

function _generate_cell_to_vertices_from_grid(grid::UnstructuredGrid,
  isperiodic::NTuple, partition)

  if is_first_order(grid)
    nodes = get_cell_node_ids(grid)
    cell_to_vertices = copy(nodes)

    nnodes = num_nodes(grid)
    num_nodes_x_dir = [partition[i]+1 for i in 1:length(partition)]
    point_to_isperiodic, slave_point_to_point, slave_point_to_master_point =
      _generate_slave_to_master_point(num_nodes_x_dir,isperiodic, nnodes)

    vertex_to_point = findall( .! point_to_isperiodic)
    point_to_vertex = fill(-1,length(point_to_isperiodic))
    point_to_vertex[vertex_to_point] = 1:length(vertex_to_point)
    point_to_vertex[slave_point_to_point] = point_to_vertex[slave_point_to_master_point]

    cell_to_vertices = Table(lazy_map(Broadcasting(Reindex(point_to_vertex)),nodes))

    vertex_to_node = vertex_to_point
    node_to_vertex = point_to_vertex
  else
    @notimplemented
  end
  (cell_to_vertices,vertex_to_node, node_to_vertex)
end

function _generate_slave_to_master_point(num_nodes_x_dir::Vector{Int},
  isperiodic::NTuple, num_nodes::Int)

  point_to_isperiodic = fill(false,num_nodes)

  slave_ijk_bounds = Array{Any,1}(undef,length(isperiodic))
  master_ijk_bounds = Array{Any,1}(undef,length(isperiodic))

  linear_indices = LinearIndices(Tuple(num_nodes_x_dir))
  periodic_dirs = findall(x->x==true, isperiodic)
  for periodic_dir in periodic_dirs
    for i in 1:length(isperiodic)
      if i == periodic_dir
        slave_ijk_bounds[i] = num_nodes_x_dir[i]
      else
        slave_ijk_bounds[i] = 1:num_nodes_x_dir[i]
      end
    end
    slave_ijk_set = CartesianIndices(Tuple(slave_ijk_bounds))
    point_to_isperiodic[linear_indices[slave_ijk_set]] .= true
  end

  slave_point_to_point = findall( point_to_isperiodic)
  slave_point_to_master_point = Array{Int32,1}(undef,length(slave_point_to_point))

  cartesian_indices = CartesianIndices(Tuple(num_nodes_x_dir))
  for (i,point) in enumerate(slave_point_to_point)
    ijk = collect(cartesian_indices[point].I)
    for i in periodic_dirs
      if ijk[i] == num_nodes_x_dir[i]
        ijk[i] = 1
      end
    end

    master_point_ijk = CartesianIndex(Tuple(ijk))
    slave_point_to_master_point[i] = linear_indices[master_point_ijk]
  end

  point_to_isperiodic, slave_point_to_point, slave_point_to_master_point
end


text:
using RobustMCTS
using Base.Test

# write your own tests here
include("grid_world.jl")


text:
###########################
# Action Space Definition #
###########################

# Terminal States
function POMDPs.isterminal(problem::SatPlanningProblem, s::SatMDPState)
    done = false

    # Have remaining actions
    if length(problem.lt_feasible_actions[(s.last_cdo_action.id, s.last_action.id)]) == 0
        done = true
    end

    # if s.power <= 0
    #     done = true
    # end

    return done
end
POMDPs.discount(problem::SatPlanningProblem) = problem.solve_discount

# Generate next state from state and action
function POMDPs.gen(problem::SatPlanningProblem, state::SatMDPState, action::Opportunity, rng::AbstractRNG)

    # Current state time
    time0 = state.time
    time  = state.time

    # Store last viable action
    last_cdo_action = state.last_cdo_action

    # Resource
    power = state.power
    data  = state.data

    # Copy observed request_ids
    request_ids = copy(state.request_ids)

    # Finish state
    done = false

    # Advance time to next action
    time = action.t_start

    # Update data generate
    power_generated = 0.0
    data_generated  = problem.spacecraft[1].datagen_backorbit*(time - time0)

    if typeof(action) == Noop
        # Do nothing if Noop
    elseif typeof(action) == Sunpoint
        # Charge for duration of sunpoint
        power_generated += problem.spacecraft[1].powergen_sunpoint*(time - time0)

    elseif typeof(action) == Collect
        # Update last action
        last_cdo_action = action

        collect_generation = action.duration * problem.spacecraft[1].datagen_image

        if data + collect_generation < 1.0
            # Only perform collect if we have capacilty
            data_generated  += collect_generation
            power_generated += action.duration * problem.spacecraft[1].powergen_image

            push!(request_ids, action.location.id)
        end

    elseif typeof(action) == Contact
        # Update last action
        last_cdo_action = action

        data_generated  += action.duration * problem.spacecraft[1].datagen_downlink
        power_generated += action.duration * problem.spacecraft[1].powergen_downlink
    else
        throw(ErrorException("Unknown Action Type: $(string(action))."))
    end

    # Update power and data numbers
    power += power_generated
    data  += data_generated

    # Ensure resource limits are enforced
    if power > 1.0 power = 1.0 end

    if power < 0.0
        power = 0.0
    end

    if data > 1.0 data = 1.0 end
    if data < 0.0 data = 0.0 end

    # Compute next state and reward for next state
    sp = SatMDPState(time=time,
        last_action=action,
        last_cdo_action=last_cdo_action,
        request_ids=request_ids,
        power=power,
        data=data
    )

    # spr = reward(problem, state, action)

    return (sp=sp,)
end

# State reward function
function POMDPs.reward(problem::SatPlanningProblem, state::SatMDPState, action::Opportunity)
    r = 0.0

    # println("Reward Action: $action")

    # Update data generate
    power_generated = 0.0
    data_generated  = problem.spacecraft[1].datagen_backorbit*(action.t_start - state.time)

    if typeof(action) == Collect
        # Resources consumed by action
        data_generated  += action.duration * problem.spacecraft[1].datagen_image
        power_generated += action.duration * problem.spacecraft[1].powergen_image

        # Only reward if we have enough spare data space
        if (state.data + data_generated < 1.0)
            # println("Action Location ID: $(action.location.id) - Requests: $(state.request_ids)")
            if !(action.location.id in state.request_ids)
                # r += action.reward

                # Semi-markov reward update
                t_diff = abs(action.t_start - state.time)
                r += action.reward*POMDPs.discount(problem)^(t_diff)
            end
            # elseif problem.solve_allow_repeats == true
            #     r += action.reward*0.25
            # end
        else
            # Penalize collection when full up
            # r -= 1000
        end

    elseif typeof(action) == Contact
        # r += 0.01*action.duration

        # Update data generation
        data_generated  += action.duration * problem.spacecraft[1].datagen_downlink
        power_generated += action.duration * problem.spacecraft[1].powergen_downlink

        n_img_downlinked = action.duration * abs(problem.spacecraft[1].datagen_downlink / problem.spacecraft[1].datagen_image)
        r += n_img_downlinked*0.01
    elseif typeof(action) == Noop
        r += 0.0
    elseif typeof(action) == Sunpoint
        # Charge for duration of sunpoint
        power_generated += problem.spacecraft[1].powergen_sunpoint*(action.t_start - state.time)

        # Have tiny bit of reward for charging
        # r += (action.t_start - state.time)/state.last_cdo_action.duration
        r += 0.0001*(action.t_start - state.time)
    end

    # Update resources to see if violations occur
    power = state.power + power_generated
    data  = state.data  + data_generated

    # Penalize running out of power
    if power <= 0.3
        r -= 10000
    end

    # if 0.2 < power <= 0.4
    #     r -= 10
    # end

    # Penalize being full on data
    if data >= 0.75
        r -= 10000
    end

    # Don't penalize overcharge
    # if power > 1
    #     reward -= 10.0
    # end

    # Don't penalize full capacity
    # if data > 1
    #     reward -= 10.0
    # end

    # Don't penalize having everything down
    # if data < 0
    #     reward -= 10.0
    # end

    # println("Reward Reward: $r")
    return r
end

# Action Space
POMDPs.actions(problem::SatPlanningProblem) = problem.actions
function POMDPs.actions(problem::SatPlanningProblem, state::SatMDPState)
    return problem.lt_feasible_actions[(state.last_cdo_action.id, state.last_action.id)]
end

# Initial state funciton
function POMDPs.initialstate(problem::SatPlanningProblem, rng::AbstractRNG)

    # Get Initial Opportunity
    init_opp = problem.opportunities[1]

    # Create Initial State
    state = SatMDPState(time=init_opp.t_start, last_action=init_opp, last_cdo_action=init_opp)

    return state
end


text:
export SymMatrix, AbstractMomentMatrix, MomentMatrix, SparseMomentMatrix
export getmat, moment_matrix, symmetric_setindex!

using SemialgebraicSets

abstract type AbstractMomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMeasureLike{T} end

"""
    mutable struct MomentMatrix{T, B<:MultivariateBases.AbstractPolynomialBasis} <: AbstractMeasureLike{T}
        Q::SymMatrix{T}
        basis::B
        support::Union{Nothing, AlgebraicSet}
    end

Measure ``\\nu`` represented by the moments of the monomial matrix ``x x^\\top`` in the symmetric matrix `Q`.
The set of points that are zeros of all the polynomials ``p`` such that ``\\mathbb{E}_{\\nu}[p] = 0`` is stored in the field `support` when it is computed.
"""
mutable struct MomentMatrix{T, B<:MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B}
    Q::SymMatrix{T}
    basis::B
    support::Union{Nothing, AlgebraicSet}
end
MomentMatrix{T, B}(Q::SymMatrix{T}, basis::MB.AbstractPolynomialBasis) where {T, B} = MomentMatrix{T, B}(Q, basis, nothing)
function MomentMatrix(Q::SymMatrix{T}, basis::MB.AbstractPolynomialBasis) where T
    return MomentMatrix{T, typeof(basis)}(Q, basis)
end

MP.variables(μ::MomentMatrix) = variables(μ.basis)
MP.nvariables(μ::MomentMatrix) = nvariables(μ.basis)

function MomentMatrix{T}(f::Function, basis::MB.AbstractPolynomialBasis, σ=1:length(basis)) where T
    return MomentMatrix(trimat(T, f, length(basis), σ), basis)
end
function MomentMatrix{T}(f::Function, monos::AbstractVector) where T
    σ, sorted_monos = sortmonovec(monos)
    return MomentMatrix{T}(f, MB.MonomialBasis(sorted_monos), σ)
end

"""
    moment_matrix(μ::Measure, x)

Creates a matrix the moment matrix for the moment matrix  ``x x^\\top`` using the moments of `μ`.
"""
function moment_matrix(μ::Measure{T}, X) where T
    function getmom(i, j)
        x = X[i] * X[j]
        for m in moments(μ)
            if monomial(m) == x
                return moment_value(m)
            end
        end
        throw(ArgumentError("μ does not have the moment $(x)"))
    end
    return MomentMatrix{T}(getmom, X)
end

function MomentMatrix(Q::AbstractMatrix{T}, basis::MB.AbstractPolynomialBasis, σ) where T
    return MomentMatrix{T}((i, j) -> Q[σ[i], σ[j]], basis)
end
function MomentMatrix(Q::AbstractMatrix, monos::AbstractVector)
    σ, sorted_monos = MP.sortmonovec(monos)
    return MomentMatrix(Q, MB.MonomialBasis(sorted_monos), σ)
end
moment_matrix(Q::AbstractMatrix, monos) = MomentMatrix(Q, monos)

function getmat(μ::MomentMatrix)
    Matrix(μ.Q)
end

function measure(ν::MomentMatrix{T, <:MB.MonomialBasis}) where T
    n = length(ν.basis)
    monos = ν.basis.elements
    measure(ν.Q.Q, [monos[i] * monos[j] for i in 1:n for j in 1:i])
end

struct SparseMomentMatrix{T, B <: MB.AbstractPolynomialBasis} <: AbstractMomentMatrix{T, B}
    sub_moment_matrices::Vector{MomentMatrix{T, B}}
end


text:
module S3

# Etag should always be quoted - in XML as well as headers...
using AWS
using AWS.Crypto
using Requests
using LightXML
using Base.get

using URIParser

import AWS: xml, elements_by_tagname

include("s3_types.jl")

def_rto = 0.0
null_cb(curl) = nothing

type RequestOptions
    blocking::Bool
    query_params::Vector{Tuple}
    request_timeout::Float64
    callback::Union{Function,Bool}
    content_type::String
    headers::Vector{Tuple}
    ostream::Union{IO, String, Void}
    auto_content_type::Bool

    RequestOptions(; blocking=true, query_params=Array{Tuple}(0), request_timeout=def_rto, callback=null_cb, content_type="", headers=Array{Tuple}(0), ostream=nothing, auto_content_type=true) =
    new(blocking, query_params, request_timeout, callback, content_type, headers, ostream, auto_content_type)
end

type Response
    body
    headers :: Dict{String, Vector{String}}
    http_code
    total_time
    bytes_recd::Integer

    Response() = new(nothing, Dict{String, Vector{String}}(), 0, 0.0, 0)
end

macro req_n_process(env, ro, resp_obj_type)
    quote
        s3_resp = do_request($(esc(env)), $(esc(ro)))
        if (isa(s3_resp.obj, String) &&  (length(s3_resp.obj) > 0))
			s3_resp.pd = LightXML.root(LightXML.parse_string(s3_resp.obj))
            s3_resp.obj = $(esc(resp_obj_type))(s3_resp.pd)
        end
        s3_resp
    end
end

type S3Response
    content_length::Int
    #Connection (open or closed) - we will not use this
    date::String    #The date and time Amazon S3 responded, for example, Wed, 01 Mar 2009 12:00:00 GMT.
    server::String  #The name of the server that created the response.
    eTag::String
    http_code::Int

    # Common amz fields here
    delete_marker::Bool
    id_2::String
    request_id::String
    version_id::String

    # all header fields
    headers::Dict

    obj::Any
    ## pd::Union{ETree, Void}
    pd::Union{LightXML.XMLElement, Void}

    S3Response() = new(0, "", "", "", 0, false, "","","",Dict{String,String}(), nothing, nothing)
end
export S3Response

type RO # RequestOptions
    verb::Symbol                  # HTTP verb
    bkt::String                   # bucket
    key::String                   # object id typically
    sub_res::Vector{Tuple}
    http_hdrs::Vector{Tuple}
    amz_hdrs::Vector{Tuple}
    cont_typ::String
    body::String
    istream::Any
    ostream::Any

    RO() = RO(:GET, "", "")
    RO(verb, bkt, key) = new(verb, bkt, key, Tuple[], Tuple[], Tuple[], "", "", nothing, nothing)
end
export RO

function list_all_buckets(env::AWSEnv)
    ro = RO(:GET, "", "")
    @req_n_process(env, ro, ListAllMyBucketsResult)
end

function del_bkt(env::AWSEnv, bkt::String)
    ro = RO(:DELETE, bkt, "")

    s3_resp = do_request(env, ro)
    s3_resp
end

function del_bkt_cors(env::AWSEnv, bkt::String)
    ro = RO(:DELETE, bkt, "")
    ro.sub_res=[("cors", "")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function del_bkt_lifecycle(env::AWSEnv, bkt::String)
    ro = RO(:DELETE, bkt, "")
    ro.sub_res=[("lifecycle", "")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function del_bkt_policy(env::AWSEnv, bkt::String)
    ro = RO(:DELETE, bkt, "")
    ro.sub_res=[("policy", "")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function del_bkt_tagging(env::AWSEnv, bkt::String)
    ro = RO(:DELETE, bkt, "")
    ro.sub_res=[("tagging", "")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function del_bkt_website(env::AWSEnv, bkt::String)
    ro = RO(:DELETE, bkt, "")
    ro.sub_res=[("website", "")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function get_bkt(env::AWSEnv, bkt::String; options::GetBucketOptions=GetBucketOptions())
    ro = RO(:GET, bkt, "")
    ro.sub_res=get_subres(Tuple[], options)

    @req_n_process(env, ro, ListBucketResult)
end

function get_bkt_acl(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("acl", "")]

    @req_n_process(env, ro, AccessControlPolicy)
end

function get_bkt_cors(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("cors", "")]

    @req_n_process(env, ro, CORSConfiguration)
end

function get_bkt_lifecycle(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("lifecycle", "")]

    s3_resp = do_request(env, ro)
    s3_resp

end

function get_bkt_policy(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("policy", "")]

    s3_resp = do_request(env, ro)
    s3_resp

end

function get_bkt_location(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("location", "")]

    @req_n_process(env, ro, LocationConstraint)
end

function get_bkt_logging(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("logging", "")]

    @req_n_process(env, ro, BucketLoggingStatus)
end

function get_bkt_notification(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("notification", "")]

    @req_n_process(env, ro, NotificationConfiguration)
end

function get_bkt_tagging(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("tagging", "")]

    @req_n_process(env, ro, Tagging)

end

function get_bkt_object_versions(env::AWSEnv, bkt::String; options::GetBucketObjectVersionsOptions=GetBucketObjectVersionsOptions())
    ro = RO(:GET, bkt, "")
    ro.sub_res=get_subres([("versions", "")], options)

    @req_n_process(env, ro, ListVersionsResult)
end

function get_bkt_request_payment(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("requestPayment", "")]

    @req_n_process(env, ro, RequestPaymentConfiguration)
end

function get_bkt_versioning(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("versioning", "")]

    @req_n_process(env, ro, VersioningConfiguration)
end

function get_bkt_website(env::AWSEnv, bkt::String)
    ro = RO(:GET, bkt, "")
    ro.sub_res=[("website", "")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function test_bkt(env::AWSEnv, bkt::String)
    ro = RO(:HEAD, bkt, "")

    s3_resp = do_request(env, ro)
    s3_resp
end

"""
split a s3 path to bucket name and key
"""
function splits3(path::AbstractString)
    path = replace(path, r"^s3://", s"")
    bkt, key = split(path, "/", limit = 2)
    return String(bkt), String(key)
end

"""
transfer s3 file to local and return local file name
`Inputs:`
awsEnv: AWS awsEnviroment
s3name: String, s3 file path
lcname: String, local temporal folder path or local file name
`Outputs:`
lcname: String, local file name
"""
function Base.download(awsEnv::AWSEnv, s3fname::AbstractString, lcfname::AbstractString)
    # directly return if not s3 file
    if !iss3(s3fname)
        return s3fname
    end

    if isdir(lcfname)
        lcfname = joinpath(lcfname, basename(s3fname))
    end
    # remove existing file
    if isfile(lcfname)
        rm(lcfname)
    elseif !isdir(dirname(lcfname))
        # create local directory
        mkdir(dirname(lcfname))
    end
    # get bucket name and key
    bkt,key = splits3(s3fname)
    # download s3 file using awscli
    f = open(lcfname, "w")
    resp = S3.get_object(awsEnv, bkt, key)
    # check that the file exist
    @assert resp.http_code == 200
    write( f, resp.obj )
    close(f)
    # run(`aws s3 cp $(s3name) $(lcfname)`)
    return lcfname
end

function get_bkt_multipart_uploads(env::AWSEnv, bkt::String; options::GetBucketUploadsOptions=GetBucketUploadsOptions())
    ro = RO(:GET, bkt, "")
    ro.sub_res=get_subres([("uploads", "")], options)

    @req_n_process(env, ro, ListMultipartUploadsResult)
end

function create_bkt(env::AWSEnv, bkt::String; acl::Union{S3_ACL, Void}=nothing, config::Union{CreateBucketConfiguration, Void}=nothing)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    if isa(acl, AccessControlPolicy)
        ro.amz_hdrs = amz_headers(Tuple[], acl)
    end
    if (config != nothing)
        ro.body = xml(config)
    end

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_acl(env::AWSEnv, bkt::String, acl:: Union{AccessControlPolicy, S3_ACL})
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res=[("acl", "")]
    if isa(acl, AccessControlPolicy)
        ro.body = xml(acl)
    else
        ro.amz_hdrs = amz_headers(Tuple[], acl)
    end

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_cors(env::AWSEnv, bkt::String, cors::CORSConfiguration)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("cors", "")]
    ro.body = xml(cors)

    s3_resp = do_request(env, ro)
    s3_resp
end

#TODO : Create LifecycleCongig type
function set_bkt_lifecycle(env::AWSEnv, bkt::String, lifecycleconfig::String)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("lifecycle", "")]
    ro.body = lifecycleconfig

    s3_resp = do_request(env, ro)
    s3_resp

end

function set_bkt_policy(env::AWSEnv, bkt::String, policy_json::String)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("policy", "")]
    ro.body = policy_json

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_logging(env::AWSEnv, bkt::String, logging::BucketLoggingStatus)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("logging", "")]
    ro.body = xml(logging)

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_notification(env::AWSEnv, bkt::String, notif::NotificationConfiguration)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("notification", "")]
    ro.body = xml(notif)

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_tagging(env::AWSEnv, bkt::String, tagging::Tagging)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("tagging", "")]
    ro.body = xml(tagging)

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_request_payment(env::AWSEnv, bkt::String, pay::RequestPaymentConfiguration)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("requestPayment", "")]
    ro.body = xml(pay)

    s3_resp = do_request(env, ro)
    s3_resp
end

function set_bkt_versioning(env::AWSEnv, bkt::String, config::VersioningConfiguration; mfa::String="")
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("versioning", "")]
    ro.body = xml(config)
    if (mfa != "") ro.amz_hdrs = [("mfa", mfa)] end

    s3_resp = do_request(env, ro)
    s3_resp
end

#TODO: Define WebsiteConfiguration as a type and use that....
function set_bkt_website(env::AWSEnv, bkt::String, websiteconfig::String)
    ro = RO(:PUT, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("website", "")]
    ro.body = websiteconfig

    s3_resp = do_request(env, ro)
    s3_resp
end

######################################
### OPERATIONS ON OBJECTS
######################################

function del_object(env::AWSEnv, bkt::String, key::String; version_id::String="", mfa::String="")
    ro = RO(:DELETE, bkt, key)
    if (version_id != "") ro.sub_res = [("versionId", version_id)] end
    if (mfa != "") ro.amz_hdrs = [("mfa", mfa)] end

    s3_resp = do_request(env, ro)
    s3_resp
end

function del_object_multi(env::AWSEnv, bkt::String, delset::DeleteObjectsType; mfa::String="")
    ro = RO(:POST, bkt, "")
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("delete", "")]
    ro.body= xml(delset)
    if (mfa != "") ro.amz_hdrs = [("mfa", mfa)] end

    @req_n_process(env, ro, DeleteResult)
end

function get_object(env::AWSEnv, bkt::String, key::String; options::GetObjectOptions=GetObjectOptions(), out::Union{IO, String, Void}=nothing, version_id::String="")
    ro = RO(:GET, bkt, key)
    ro.sub_res = Tuple[]

    if (version_id != "") push!(ro.sub_res, ("versionId", version_id)) end
    ro.sub_res = query_params(ro.sub_res, options)

    ro.http_hdrs = http_headers(Tuple[], options)
    ro.ostream = out

    s3_resp = do_request(env, ro, conv_to_string=false)
    s3_resp
end

function get_object_acl(env::AWSEnv, bkt::String, key::String; version_id::String="")
    ro = RO(:GET, bkt, key)
    ro.sub_res=[("acl", "")]
    if (version_id != "") push!(ro.sub_res, ("versionId", version_id)) end

    @req_n_process(env, ro, AccessControlPolicy)
end

function get_object_torrent(env::AWSEnv, bkt::String, key::String, out::Union{IO, String}) # out is either an IOStream or a filename
    ro = RO(:GET, bkt, key)
    ro.sub_res=[("torrent","")]
    ro.ostream = out

    s3_resp = do_request(env, ro, conv_to_string=false)
    s3_resp
end

function describe_object(env::AWSEnv, bkt::String, key::String; options::GetObjectOptions=GetObjectOptions(), version_id::String="")
    ro = RO(:HEAD, bkt, key)
    ro.http_hdrs = http_headers(Tuple[], options)
    if (version_id != "") ro.sub_res=[("versionId", version_id)] end

    s3_resp = do_request(env, ro)
    s3_resp
end

function test_object(env::AWSEnv, bkt::String, key::String, origin::String, acc_ctrl_req_method::String; acc_ctrl_req_hdrs::Union{String, Void}=nothing)
    ro = RO(:OPTIONS, bkt, key)
    ro.http_hdrs = [("Origin", origin), ("Access-Control-Request-Method",acc_ctrl_req_method)]
    if headers != nothing
        push!(ro.http_hdrs, ("Access-Control-Request-Headers", acc_ctrl_req_hdrs))
    end

    s3_resp = do_request(env, ro)
    s3_resp
end

#post_object, a different version of put_object is not supported

function restore_object(env::AWSEnv, bkt::String, key::String, days::Int)
    # TODO qp = "restore"
    # convert days to RestoreRequest
    ro = RO(:POST, bkt, key)
	ro.cont_typ = "application/octet-stream"

    ro.sub_res=[("restore", "")]
    ro.body = "<RestoreRequest xmlns=\"http://$(amz_s3_endpoint(env))/doc/2006-3-01\"><Days>$(days)</Days></RestoreRequest>"

    s3_resp = do_request(env, ro)
    s3_resp
end

function put_object(env::AWSEnv, bkt::String, key::String, data:: Union{IO, String, Tuple}; content_type="", options::PutObjectOptions=PutObjectOptions(), version_id::String="")
    ro = RO(:PUT, bkt, key)
	ro.cont_typ = "application/octet-stream"

    ro.amz_hdrs = amz_headers(Tuple[], options)
    ro.http_hdrs = http_headers(Array{Tuple}(0), options)
    if (content_type != "") ro.cont_typ = content_type end

    if isa(data, String)
        ro.body = data
    elseif isa(data, IO)
        ro.istream = data
    elseif (isa(data, Tuple) && data[1] == :file)
        ro.istream = data[2]
    else
        error("Void to upload")
    end

    if (version_id != "") ro.sub_res=[("versionId", version_id)] end

    s3_resp = do_request(env, ro)
    s3_resp
end

function put_object_acl(env::AWSEnv, bkt::String, key::String, acl:: Union{AccessControlPolicy, S3_ACL})
    ro = RO(:PUT, bkt, key)
	ro.cont_typ = "application/octet-stream"
    ro.sub_res=[("acl", "")]
    if isa(acl, AccessControlPolicy)
        ro.body = xml(acl)
    else
        ro.amz_hdrs = amz_headers(Tuple[], acl)
    end

    s3_resp = do_request(env, ro)
    s3_resp
end

function copy_object(env::AWSEnv, dest_bkt::String, dest_key::String, options::CopyObjectOptions, version_id::String="")
    ro = RO(:PUT, dest_bkt, dest_key)
	ro.cont_typ = "application/octet-stream"
    ro.amz_hdrs = amz_headers(Tuple[], options)
    if (version_id != "")
        ro.sub_res=[("versionId", version_id)]
    end

    @req_n_process(env, ro, CopyObjectResult)
end

function initiate_multipart_upload(env::AWSEnv, bkt::String, key::String; content_type="", options::PutObjectOptions=PutObjectOptions())
    ro = RO(:POST, bkt, key)
	ro.cont_typ = "application/octet-stream"
    ro.amz_hdrs = amz_headers(Tuple[], options)
    ro.http_hdrs = http_headers(Array{Tuple}(0), options)
    if (content_type != "") ro.cont_typ = content_type end
    ro.sub_res=[("uploads", "")]

    @req_n_process(env, ro, InitiateMultipartUploadResult)
end

function upload_part(env::AWSEnv, bkt::String, key::String, part_number::String, upload_id::String, data:: Union{IO, String, Tuple})
    ro = RO(:PUT, bkt, key)
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("partNumber", "$(part_number)"), ("uploadId", "$(upload_id)")]

    if isa(data, String)
        ro.body = data
    elseif isa(data, IO)
        ro.istream = data
    elseif (isa(data, Tuple) && data[1] == :file)
        ro.istream = data[2]
    else
        error("Void to upload")
    end

    s3_resp = do_request(env, ro)
    s3_resp

end

function copy_upload_part(env::AWSEnv, bkt::String, key::String, part_number::String, upload_id::String, options::CopyUploadPartOptions)
    ro = RO(:PUT, bkt, key)
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("partNumber", "$(part_number)"), ("uploadId", "$(upload_id)")]
    ro.amz_hdrs = amz_headers(options)

    @req_n_process(env, ro, CopyPartResult)
end

function complete_multipart_upload(env::AWSEnv, bkt::String, key::String, upload_id::String, part_ids::Vector{S3PartTag})
    ro = RO(:POST, bkt, key)
	ro.cont_typ = "application/octet-stream"
    ro.sub_res = [("uploadId", "$(upload_id)")]

    ro.body = xml("CompleteMultipartUpload", part_ids)

    @req_n_process(env, ro, CompleteMultipartUploadResult)
end

function abort_multipart_upload(env::AWSEnv, bkt::String, key::String, upload_id::String)
    ro = RO(:DELETE, bkt, key)
    ro.sub_res = [("uploadId", "$(upload_id)")]

    s3_resp = do_request(env, ro)
    s3_resp
end

function list_upload_parts(env::AWSEnv, bkt::String, key::String,
    upload_id::String;
    max_parts:: Union{Int, Void}=nothing,
    part_number_marker:: Union{Int, Void}=nothing)

    ro = RO(:GET, bkt, key)
    ro.sub_res = [("uploadId", "$(upload_id)"), ("max-parts", "$max_parts"), ("part-number-marker", "$part_number_marker")]

    @req_n_process(env, ro, ListPartsResult)
end

function do_request(env::AWSEnv, ro::RO; conv_to_string=true)
    http_resp = do_http(env, ro)

    s3_resp = S3Response()

    s3_resp.http_code = http_resp.status
    cl_s = get(http_resp.headers, "Content-Length", [""])[1]
    if cl_s != ""
        s3_resp.content_length = Base.parse(Int, cl_s)
    end

    s3_resp.date = get(http_resp.headers, "Date", [""])
    s3_resp.server = get(http_resp.headers, "Server", [""])

    ## s3_resp.delete_marker = lowercase(get(http_resp.headers, "x-amz-delete-marker", ["false"])[1]) == "true" ? true : false
    ## s3_resp.id_2 = get(http_resp.headers, "x-amz-id-2", [""])[1]
    ## s3_resp.request_id = get(http_resp.headers, "x-amz-request-id", [""])[1]
    ## s3_resp.version_id = get(http_resp.headers, "x-amz-version-id", [""])[1]

    s3_resp.headers = http_resp.headers

    if (http_resp.status > 299)
        if  (search(Base.get(http_resp.headers, "Content-Type", [""]), "/xml") != 0:-1)
            s3_resp.obj = S3Error(LightXML.root(LightXML.parse_string(String(copy(http_resp.data)))))
        else
            s3_resp.pd = String(copy(http_resp.data))
        end
    else
        s3_resp.obj = conv_to_string ? String(copy(http_resp.data)) : http_resp.data
    end

    s3_resp
end

function do_http(env::AWSEnv, ro::RO)
    if isa(ro.body, String) && (ro.body != "")
#        digest = zeros(UInt8, 16)
#        MD5(body, length(body), digest)
        md5 = base64encode(Crypto.md5(ro.body))
    elseif isa(ro.istream, IO)
        # Read the entire istream to get the MD5 of the same.
        md5 = base64encode(Crypto.md5(ro.istream))
        seekstart(ro.istream)
    elseif isa(ro.istream, String)
        # The file will be read twice (once to get the MD5 and once while sending - no other way?
        md5 = base64encode(Crypto.md5_file(ro.istream))
    else
        md5 = ""
    end

    (new_amz_hdrs, full_path, s_b64) = canonicalize_and_sign(env, ro, md5)
	## @show new_amz_hdrs, full_path, s_b64

    all_hdrs = new_amz_hdrs
    (md5 != "") ? push!(all_hdrs, ("Content-MD5", md5)) : nothing
    (ro.cont_typ != "") ? push!(all_hdrs, ("Content-Type", ro.cont_typ)) : nothing

    # Remove Content-MD5 and Expect headers from the passed http_hdrs
    for (name, value) in ro.http_hdrs
        lname = lowercase(strip(name))
        if (lname != "content-md5") && (lname != "expect") && (lname != "content-type")
            push!(all_hdrs, (name, value))
        elseif (lname == "content-type") && (ro.cont_typ == "")
            # If content type has been set as part of RO, it will be used.
            push!(all_hdrs, (name, value))
        end
    end

    push!(all_hdrs, ("Authorization",  "AWS " * env.aws_id * ":" * s_b64))

    url = "https://$(amz_s3_endpoint(env))" * full_path

    http_options = RequestOptions(headers=all_hdrs, ostream=ro.ostream, request_timeout=env.timeout, auto_content_type=false)

    if env.dbg
        info("Verb : " * string(ro.verb))
        info("URL : " * url)
        info("Headers : ")
        for (k,v) in http_options.headers
            info("  $k : $v")
        end
        info("Body : " * string(ro.body))
    end

    if env.dry_run
        return Requests.Response()
    end

	headers = Dict(
        "User-Agent" => "Requests.jl/0.0.0",
        "Accept" => "text/html,application/xhtml+xml,application/xml;q=0.9,*/*;q=0.8"
    )
	for (k,v) in all_hdrs
        ## info("  $k : $v")
		setindex!(headers, v, k)
	end

	## @show headers
	## TODO To include ostream, auto_content_type and timeout values for all VERBS !

    if ro.verb == :GET
        http_resp = Requests.get(url; headers = headers)
    elseif (ro.verb == :PUT) || (ro.verb == :POST)
        senddata = isa(ro.istream, IO) ? ro.istream :
                   isa(ro.istream, String) ? (:file, ro.istream) :
                   (isa(ro.body, String) && (ro.body != "")) ? ro.body :
                   ""
#                   error("Must specify either a body or istream for PUT/POST")

        if (ro.verb == :PUT)
            http_resp = Requests.put(URIParser.parse_url(url), senddata; headers=headers)
        else
            http_resp = Requests.post(URIParser.parse_url(url), senddata; headers=headers)
        end

    elseif ro.verb == :DELETE
        http_resp = Requests.delete(URIParser.parse_url(url); headers=headers)

    elseif ro.verb == :HEAD
        http_resp = Requests.head(URIParser.parse_url(url); headers=headers)

    else
        error("Unknown HTTP verb $verb")
    end

     if env.dbg
         info("Response : ", http_resp)
     end

    return http_resp
end

function canonicalize_and_sign(env::AWSEnv, ro::RO, md5::String)
    if env.aws_token != ""
        push!(ro.amz_hdrs, ("x-amz-security-token", env.aws_token))
    end
    (new_amz_hdrs, amz_hdrs_str) = get_canon_amz_headers(ro.amz_hdrs)
    (full_path, sign_path) = get_canonicalized_resource(ro)

    str2sign = string(ro.verb) * "\n" *
    md5 * "\n" *
    ro.cont_typ * "\n" * "\n" * # Using x-amz-date instead of Date.
    amz_hdrs_str * sign_path

    if (env.dbg)
        info("String to sign:\n$(str2sign)")
    end

    s_b64 = base64encode(Crypto.hmacsha1_digest(str2sign, env.aws_seckey))

    return new_amz_hdrs, full_path, s_b64
end

const qstr_sign_list = Set([
    "acl",
    "lifecycle",
    "location",
    "logging",
    "notification",
    "partNumber",
    "policy",
    "requestPayment",
    "torrent",
    "uploadId",
    "uploads",
    "versionId",
    "versioning",
    "versions",
    "website",
    "response-content-type",
    "response-content-language",
    "response-expires",
    "response-cache-control",
    "response-content-disposition",
    "response-content-encoding",
    "delete"]
)

function get_canonicalized_resource(ro::RO)
    if length(ro.bkt) > 0
        part1 = (startswith(ro.bkt, "/") ? "" : "/") * ro.bkt
        if length(ro.key) > 0
            part1 = part1 * (startswith(ro.key, "/") ? "" : "/") * ro.key
        else
            part1 = part1 * "/"
        end
    else
        part1 = "/"
    end

    sorted = sort(ro.sub_res)
    signlist = filter(x -> begin (k,v) = x; in(k, qstr_sign_list) end, sorted)

    signparams = ""
    for x in signlist
        k,v = x;
        if (v != "")
            ep = string(k) * "=" * string(v)
        else
            ep = string(k)
        end
        if signparams == ""
            signparams = ep
        else
            signparams = signparams  * "&" * ep
        end
    end

    canon_res = part1 * "?" * Requests.format_query_str(sorted)

    canon_sign_res = part1
    if length(signlist) > 0
        canon_sign_res = canon_sign_res * "?" * signparams
    end

    return (canon_res, canon_sign_res)
end

function get_canon_amz_headers(headers::Vector{Tuple})
    lcase = [(lowercase(strip(k)), v) for (k,v) in headers]

    # remove any existing amz-date header
    lcase = filter(x -> begin k,v = x; k != "x-amz-date" end, lcase)

    push!(lcase, ("x-amz-date", rfc1123_date()))

    reduced = Dict{String, String}()
    for (k,v) in lcase
        if startswith(k,"x-amz-")
            new_v = strip(replace(v, "\n", ' '))
            if haskey(reduced, k)
                ev = reduced[k]
                reduced[k] = ev * "," * new_v
            else
                reduced[k] = new_v
            end
        end
    end

    # Use the sorted one in the final request too since the order of 'values'
    sorted = sort([(x[1], x[2]) for x in collect(reduced)])
    canon_hdr_str = ""
    for x in sorted
        (k,v) = x
        canon_hdr_str = canon_hdr_str * k * ":" * v * "\n"
    end

    return sorted, canon_hdr_str
end

amz_s3_endpoint(region::String) = ((isempty(region) || (region == "us-east-1")) ? "s3" : ("s3-"*region)) * ".amazonaws.com"
amz_s3_endpoint(env::AWSEnv) = isempty(env.ep_host) ? amz_s3_endpoint(env.region) : env.ep_host

rfc1123_date() = rfc1123_date(now(Dates.UTC))
function rfc1123_date(dt)
    y,m,d = Dates.yearmonthday(dt)
    h,mi,s = Dates.hour(dt),Dates.minute(dt),Dates.second(dt)

    dayofweek = Dates.dayabbr(dt)
    dd = lpad(d,2,"0")
    monthname = Dates.monthabbr(dt)
    yy = y < 0 ? @sprintf("%05i",y) : lpad(y,4,"0")
    hh = lpad(h,2,"0")
    mii = lpad(mi,2,"0")
    ss = lpad(s,2,"0")

    #format("EEE, dd MMM yyyy HH:mm:ss V")
    # "Tue, 24 Feb 2015 05:21:51 GMT"
    return "$dayofweek, $dd $monthname $yy $hh:$mii:$ss GMT"
end
rfc1123_date(d::Void) = nothing

"""
list objects of s3. no directory/folder in the list
`Inputs`:
awsEnv: AWS awsEnvironment
bkt: bucket name
path: path
`Outputs`:
ret: channel that produces object keys (String)
"""
function list_objects(awsEnv::AWSEnv, path::AbstractString; re::Regex = r"^\s*")
    bkt, prefix = splits3(path)
    return list_objects(awsEnv, bkt, prefix; re=re)
end
function list_objects(path::AbstractString; re::Regex = r"^\s*")
    bkt, prefix = splits3(path)
    return list_objects(bkt, prefix; re=re)
end
function list_objects(bkt::AbstractString, prefix::AbstractString; re::Regex = r"^\s*")
    awsEnv = AWS.AWSEnv()
    return list_objects(awsEnv, bkt, prefix; re=re)
end
function list_objects(awsEnv::AWSEnv, bkt::AbstractString, prefix::AbstractString; re::Regex = r"^\s*")
    (bkt, in_c) = list_objects_details(awsEnv, bkt, prefix)

    function filter(out_c::Channel)
        for content in in_c
            if startswith(content.key, prefix)
                fname = content.key[length(prefix)+1:end]
            else
                fname = content.key
            end

            if ismatch(re, fname)
                put!(out_c, content.key)
            end
        end
    end

    bkt, Channel(filter, ctype=String)
end

"""
list objects of s3. no directory/folder in the list
`Inputs`:
awsEnv: AWS awsEnvironment
bkt: bucket name
path: path
`Outputs`:
ret: channel with object metadata (AWS.S3.Contents)
"""
function list_objects_details(awsEnv::AWSEnv, bkt::AbstractString, prefix::AbstractString; max_keys = 1000)
    prefix = lstrip(prefix, '/')

    function producer(c::Channel)
        marker = nothing
        more = true

        while more
            bucket_options = AWS.S3.GetBucketOptions(delimiter="/", prefix=prefix, marker=marker, max_keys=max_keys)
            resp = AWS.S3.get_bkt(awsEnv, bkt; options=bucket_options)

            more = resp.obj.isTruncated
            marker = resp.obj.nextMarker

            put!.(c, resp.obj.contents)
        end
    end

    # csize=max_keys-1: start the next listing as soon as the caller starts
    # processing the previous listing
    return bkt, Channel(producer, ctype=AWS.S3.Contents, csize=max(0, max_keys-1))
end

end