text:
test_context("Data types and NA's")

test_group("NA's")
@assert length(NA) == 1
@assert size(NA) == ()
@assert isna(3 == NA)
@assert isna(NA == 3)
@assert isna(NA == NA)

test_group("DataVector creation")
dvint = DataVector[1, 2, NA, 4]
dvint2 = DataArray([5:8])
dvint3 = DataArray(5:8)
dvflt = DataVector[1.0, 2, NA, 4]
dvstr = DataVector["one", "two", NA, "four"]

@assert isa(dvint, DataVector{Int64})
@assert isa(dvint2, DataVector{Int64})
@assert isa(dvint3, DataVector{Int64})
@assert isa(dvflt, DataVector{Float64})
@assert isa(dvstr, DataVector{ASCIIString})
@test throws_exception(DataArray([5:8], falses(2)), Exception)

@assert isequal(DataArray(dvint), dvint)

test_group("PooledDataVector creation")
pdvstr = PooledDataVector["one", "one", "two", "two", NA, "one", "one"]
@assert isa(pdvstr, PooledDataVector{ASCIIString})
@test throws_exception(PooledDataVector["one", "one", 9], Exception)
@assert isequal(PooledDataArray(pdvstr), pdvstr)

test_group("PooledDataVector creation with predetermined pool")
pdvpp = PooledDataArray([1, 2, 2, 3], [1, 2, 3, 4])
@assert isequal(pdvpp.pool, [1, 2, 3, 4])
@assert string(pdvpp) == "[1, 2, 2, 3]"
@test throws_exception(PooledDataArray([1, 2, 3], [1, 2]), Exception)
pdvpp = PooledDataArray([1, 2, 2, 3, 2, 1], [1, 2, 3, 4])
@assert isequal(pdvpp.pool, [1, 2, 3, 4])
@assert string(pdvpp) == "[1, 2, 2, 3, 2, 1]"
pdvpp = PooledDataArray(["one", "two", "two"], ["one", "two", "three"])
@assert isequal(values(pdvpp), DataVector["one", "two", "two"])
@assert all(get_indices(pdvpp) .== uint16([1, 3, 3]))
@assert isequal(levels(pdvpp), DataVector["one", "three", "two"])
@assert isequal(pdvpp.pool, ["one", "three", "two"])
@assert string(pdvpp) == "[one, two, two]"
@test throws_exception(PooledDataArray(["one", "two", "four"], ["one", "two", "three"]), Exception)

test_group("PooledDataVector utf8 support")
pdvpp = PooledDataArray([utf8("hello")], [false])
@assert isa(pdvpp[1], UTF8String)
pdvpp = PooledDataArray([utf8("hello")])
@assert isa(pdvpp[1], UTF8String)

test_group("DataVector access")
@assert dvint[1] == 1
@assert isna(dvint[3])
@assert isequal(dvflt[3:4], DataVector[NA, 4.0])
@assert isequal(dvint[[true, false, true, false]], DataVector[1, NA])
@assert isequal(dvstr[[1, 2, 1, 4]], DataVector["one", "two", "one", "four"])
@assert isequal(dvstr[[1, 2, 1, 3]], DataVector["one", "two", "one", NA])

test_group("PooledDataVector access")
@assert pdvstr[1] == "one"
@assert isna(pdvstr[5])
@assert isequal(pdvstr[1:3], DataVector["one", "one", "two"])
@assert isequal(pdvstr[[true, false, true, false, true, false, true]], PooledDataVector["one", "two", NA, "one"])
@assert isequal(pdvstr[[1, 3, 1, 2]], DataVector["one", "two", "one", "one"])

test_group("DataVector methods")
@assert size(dvint) == (4,)
@assert length(dvint) == 4
@assert sum(isna(dvint)) == 1
@assert eltype(dvint) == Int64

test_group("PooledDataVector methods")
@assert size(pdvstr) == (7,)
@assert length(pdvstr) == 7
@assert sum(isna(pdvstr)) == 1
@assert eltype(pdvstr) == ASCIIString

test_group("DataVector operations")
@assert isequal(dvint + 1, DataArray([2, 3, 4, 5], [false, false, true, false]))
@assert isequal(dvint .* 2, DataVector[2, 4, NA, 8])
@assert isequal(dvint .== 2, DataVector[false, true, NA, false])
@assert isequal(dvint .> 1, DataVector[false, true, NA, true])

test_group("PooledDataVector operations")
# @assert isequal(pdvstr .== "two", PooledDataVector[false, false, true, true, NA, false, false])

test_group("DataVector to something else")
@assert all(removeNA(dvint) .== [1, 2, 4])
@assert all(replaceNA(dvint, 0) .== [1, 2, 0, 4])
@assert all(convert(Vector{Int}, dvint2) .== [5:8])
@assert all([i + 1 for i in dvint2] .== [6:9])
@assert all([length(x)::Int for x in dvstr] == [3, 3, 1, 4])
@assert repr(dvint) == "[1, 2, NA, 4]"

test_group("PooledDataVector to something else")
@assert all(removeNA(pdvstr) .== ["one", "one", "two", "two", "one", "one"])
@assert all(replaceNA(pdvstr, "nine") .== ["one", "one", "two", "two", "nine", "one", "one"])
@assert all([length(i)::Int for i in pdvstr] .== [3, 3, 3, 3, 1, 3, 3])
@assert string(pdvstr[1:3]) == "[one, one, two]"

test_group("DataVector Filter and Replace")
@assert isequal(removeNA(dvint), [1, 2, 4])
@assert isequal(replaceNA(dvint, 7), [1, 2, 7, 4])
@assert sum(removeNA(dvint)) == 7
@assert sum(replaceNA(dvint, 7)) == 14

test_group("PooledDataVector Filter and Replace")
@assert reduce(string, "", removeNA(pdvstr)) == "oneonetwotwooneone"
@assert reduce(string, "", replaceNA(pdvstr,"!")) == "oneonetwotwo!oneone"

test_group("DataVector assignment")
assigntest = DataVector[1, 2, NA, 4]
assigntest[1] = 8
@assert isequal(assigntest, DataVector[8, 2, NA, 4])
assigntest[1:2] = 9
@assert isequal(assigntest, DataVector[9, 9, NA, 4])
assigntest[[1,3]] = 10
@assert isequal(assigntest, DataVector[10, 9, 10, 4])
assigntest[[true, false, true, true]] = 11
@assert isequal(assigntest, DataVector[11, 9, 11, 11])
assigntest[1:2] = [12, 13]
@assert isequal(assigntest, DataVector[12, 13, 11, 11])
assigntest[[1, 4]] = [14, 15]
@assert isequal(assigntest, DataVector[14, 13, 11, 15])
assigntest[[true, false, true, false]] = [16, 17]
@assert isequal(assigntest, DataVector[16, 13, 17, 15])
assigntest[1] = NA
@assert isequal(assigntest, DataVector[NA, 13, 17, 15])
assigntest[[1, 2]] = NA
@assert isequal(assigntest, DataVector[NA, NA, 17, 15])
assigntest[[true, false, true, false]] = NA
@assert isequal(assigntest, DataVector[NA, NA, NA, 15])
assigntest[1] = 1
assigntest[2:4] = NA
@assert isequal(assigntest, DataVector[1, NA, NA, NA])

test_group("PooledDataVector assignment")
ret = (pdvstr[2] = "three")
@assert ret == "three"
@assert pdvstr[2] == "three"
ret = pdvstr[[1,2]] = "two"
@assert ret == "two"
@assert pdvstr[2] == "two"
pdvstr2 = PooledDataVector["one", "one", "two", "two"]
ret = (pdvstr2[[true, false, true, false]] = "three")
@assert ret == "three"
@assert pdvstr2[1] == "three"
ret = (pdvstr2[[false, true, false, true]] = ["four", "five"])
@assert isequal(ret, ["four", "five"])
@assert isequal(pdvstr2[3:4], DataVector["three", "five"])
pdvstr2 = PooledDataVector["one", "one", "two", "two"]
ret = (pdvstr2[2:3] = "three")
@assert ret == "three"
@assert isequal(pdvstr2[3:4], DataVector["three", "two"])
ret = (pdvstr2[2:3] = ["four", "five"])
@assert ret == ["four", "five"]
@assert isequal(pdvstr2[1:2], DataVector["one", "four"])
pdvstr2 = PooledDataVector["one", "one", "two", "two", "three"]
@assert isna(begin pdvstr2[1] = NA end)
@assert all(isna(begin pdvstr2[[1, 2]] = NA end))
@assert all(isna(begin pdvstr2[[false, false, true, false, false]] = NA end))
@assert all(isna(begin pdvstr2[4:5] = NA end))
@assert all(isna(pdvstr2))

test_group("PooledDataVector replace!")
pdvstr2 = PooledDataVector["one", "one", "two", "two", "three"]
@assert replace!(pdvstr2, "two", "four") == "four"
@assert replace!(pdvstr2, "three", "four") == "four"
@assert isna(replace!(pdvstr2, "one", NA))
@assert replace!(pdvstr2, NA, "five") == "five"
@assert isequal(pdvstr2, DataVector["five", "five", "four", "four", "four"])

test_context("DataFrames")

test_group("constructors")
df1 = DataFrame({dvint, dvstr}, ["Ints", "Strs"])
df2 = DataFrame({dvint, dvstr})
df3 = DataFrame({dvint})
df4 = DataFrame([1:4 1:4])
df5 = DataFrame({DataVector[1,2,3,4], dvstr})
df6 = DataFrame({dvint, dvint, dvstr}, ["A", "B", "C"])

test_group("description functions")
@assert nrow(df6) == 4
@assert ncol(df6) == 3
@assert all(colnames(df6) .== ["A", "B", "C"])
@assert all(colnames(df2) .== ["x1", "x2"])

test_group("ref")
@assert df6[2, 3] == "two"
@assert isna(df6[3, 3])
@assert df6[2, "C"] == "two"
@assert isequal(df6["B"], dvint)
@assert ncol(df6[[2,3]]) == 2
@assert nrow(df6[2,:]) == 1
@assert size(df6[[1, 3], [1, 3]]) == (2, 2)
@assert size(df6[1:2, 1:2]) == (2, 2)
@assert size(head(df6,2)) == (2, 3)
# lots more to do

test_group("combining")

dfc = cbind(df3, df4)
@assert ncol(dfc) == 3
@assert all(colnames(dfc) .== ["x1", "x1_1", "x2"])
@assert isequal(dfc["x1"], df3["x1"])

@assert isequal(dfc, [df3 df4])

dfr = rbind(df4, df4)
@assert nrow(dfr) == 8
@assert all(colnames(df4) .== colnames(dfr))
@assert isequal(dfr, [df4, df4])

test_group("show")
@assert repr(df1) == "4x2 DataFrame:\n        Ints   Strs\n[1,]       1  \"one\"\n[2,]       2  \"two\"\n[3,]      NA     NA\n[4,]       4 \"four\"\n"

test_group("assign")
df6[3] = DataVector["un", "deux", "troix", "quatre"]
@assert df6[1, 3] == "un"
df6["B"] = [4, 3, 2, 1]
@assert df6[1,2] == 4
df6["D"] = [true, false, true, false]
@assert df6[1,4] == true
delete!(df6, "D")
@assert all(colnames(df6) .== ["A", "B", "C"])
@assert ncol(df6) == 3

test_group("NA handling")
@assert nrow(df5[complete_cases(df5), :]) == 3

test_context("SubDataFrames")

test_group("constructors")
# single index is rows
sdf6a = sub(df6, 1)
sdf6b = sub(df6, 2:3)
sdf6c = sub(df6, [true, false, true, false])
@assert size(sdf6a) == (1,3)
sdf6d = sub(df6, [1,3], "B")
@assert size(sdf6d) == (2,1)

test_group("ref")
@assert sdf6a[1,2] == 4

test_context("Within")
test_group("Associative")

srand(1)
a1 = [:a => [1, 2], :b => [3, 4], :c => [5, 6]]
a2 = ["a" => [1, 2], "b" => [3, 4], "c" => [5, 6]]
a3 = {"a" => [1, 2], "b" => [3, 4], :c => [5, 6]}

@assert isequal(with(a1, :(c + 1)), a1[:c] + 1)
@assert isequal(with(a2, :(c + 1)), with(a1, :(c + 1)))
@assert isequal(with(a3, :(c + 1)), with(a1, :(c + 1)))
@assert isequal(with(a3, :(c + 1 + 0 * b)), with(a1, :(c + 1)))

a4 = within(a1, :( d = a + b ))
@assert isequal(a4[:d], a1[:a] + a1[:b])
@assert isequal(a4[:a], a1[:a])

a4c = @transform(copy(a1), d => a + b )
@assert isequal(a4[:d], a4c[:d])

a4 = within(a2, :( d = a + b ))
@assert isequal(a4["d"], a2["a"] + a2["b"])
@assert isequal(a4["a"], a2["a"])

a4c = @transform(copy(a2), d => a + b )
@assert isequal(a4c["d"], a4["d"])

a4 = within(a3, :( d = a + b ))
@assert isequal(a4[:d], a3["a"] + a3["b"])
@assert isequal(a4["a"], a3["a"])

# Note: The following won't work.
#       @transform will only find the keys that are symbols.
## a4c = @transform(copy(a3), d => a + b )
## @assert isequal(a4c[:d], a4[:d])

a4 = based_on(a1, :( d = a + b ))
@assert isequal(a4[:d], a1[:a] + a1[:b])

a4 = based_on(a2, :( d = a + b ))
@assert isequal(a4["d"], a2["a"] + a2["b"])

a4 = based_on(a3, :( d = a + b ))
@assert isequal(a4[:d], a3["a"] + a3["b"])

test_group("DataFrame")

srand(1)
N = 20
d1 = PooledDataArray(rand(1:2, N))
d2 = PooledDataVector["A", "B", NA][rand(1:3, N)]
d3 = DataArray(randn(N))
d4 = DataArray(randn(N))
df7 = DataFrame({d1, d2, d3}, ["d1", "d2", "d3"])

@assert isequal(with(df7, :(d3 + d3)), df7["d3"] + df7["d3"])
@assert isequal(with(df7, :(d3 + $d4)), df7["d3"] + d4)
x = df7 | with(:( d3 + d3 ))
@assert isequal(x, df7["d3"] + df7["d3"])

df8 = within(df7, :(d4 = d3 + d3 + 1))
@assert isequal(df7, df8[1:3])
@assert isequal(df8["d4"], df7["d3"] + df7["d3"] + 1)
within!(df8, :( d4 = d1 ))
@assert isequal(df8["d1"], df8["d4"])

df8 = @transform(copy(df7), d4 => d3 + 1)
@assert isequal(df7, df8[1:3])

df8 = based_on(df7, :( d1 = d3 ))
@assert isequal(df8["d1"], df7["d3"])
df8 = df7 | based_on(:( d1 = d3 ))
@assert isequal(df8["d1"], df7["d3"])
df8 = based_on(df7, :( sum_d3 = sum(d3) ))
@assert isequal(df8[1,1], sum(df7["d3"]))

#@assert all(df7[:( d2 .== "B" )]["d1"] .== PooledDataArray([1,2,1,1]))
# TODO: Remove all tests that depend upon srand(), which was just changed.
@assert all(df7[:( d2 .== "B" ), "d1"] .== PooledDataArray([2,1,1,1,1,2,2,1,2]))

test_group("groupby")

gd = groupby(df7, "d1")
@assert length(gd) == 2
# @assert isequal(gd[2]["d2"], PooledDataVector["A", "B", NA, "A", NA, NA, NA, NA])
@assert sum(gd[2]["d3"]) == sum(df7["d3"][removeNA(df7["d1"] .== 2)])

g1 = groupby(df7, ["d1", "d2"])
g2 = groupby(df7, ["d2", "d1"])
@assert sum(g1[1]["d3"]) == sum(g2[1]["d3"])

res = 0.0
for x in g1
    res += sum(x["d1"])
end
@assert res == sum(df7["d1"])

df8 = df7 | groupby(["d2"]) | :( d3sum = sum(d3); d3mean = mean(removeNA(d3)) )
@assert isequal(df8["d2"], PooledDataVector[NA, "A", "B"])

df9 = based_on(groupby(df7, "d2"),
               :( d3sum = sum(d3); d3mean = mean(removeNA(d3)) ))
@assert isequal(df9, df8)

df8 = within(groupby(df7, "d2"),
             :( d4 = d3 + 1; d1sum = sum(d1) ))
@assert all(df8[:( d2 .== "C" ), "d1sum"] .== 13)

## @assert isequal(with(g1, :( sum(d1) )), map(x -> sum(x["d1"]), g1))

df8 = colwise(df7[[1, 3]], :sum)
@assert df8[1, "d1_sum"] == sum(df7["d1"])

df8 = colwise(groupby(df7, "d2"), [:sum, :length])
@assert nrow(df8) == 3
@assert ncol(df8) == 5
#@assert df8[1, "d1_sum"] == 13
@assert df8[2, "d1_length"] == 8

df9 = df7 | groupby(["d2"]) | [:sum, :length]
@assert isequal(df9, df8)
df9 = by(df7, "d2", [:sum, :length])
@assert isequal(df9, df8)

test_group("reshape")

d1 = DataFrame(quote
    a = [1:3]
    b = [1:4]
    c = randn(12)
    d = randn(12)
end)

d1c = @DataFrame(a => [1:3],
                 b => [1:4],
                 c => randn(12),
                 d => randn(12))

@assert isequal(d1[1:2], d1c[1:2])

d1s = stack(d1, ["a", "b"])
d1s2 = stack(d1, ["c", "d"])
@assert isequal(d1s[1:12, "c"], d1["c"])
@assert isequal(d1s[13:24, "c"], d1["c"])
@assert all(colnames(d1s) .== ["variable", "value", "c", "d"])
d1s_df = stack_df(d1, ["a", "b"])
@assert isequal(d1s["variable"], d1s_df["variable"][:])
@assert isequal(d1s["value"], d1s_df["value"][:])
@assert isequal(d1s["c"], d1s_df["c"][:])
@assert isequal(d1s[1,:], d1s_df[1,:])

d1s["idx"] = [1:12, 1:12]
d1s2["idx"] = [1:12, 1:12]
d1us = unstack(d1s, "variable", "idx", "value")
d1us2 = unstack(d1s2, "variable", "idx", "value")
@assert isequal(d1us["a"], d1["a"])
@assert isequal(d1us2["d"], d1["d"])

d = DataFrame(quote
    a = letters[5:8]
    b = LETTERS[10:11]
    c = LETTERS[13 + [1, 1, 2, 2, 2, 1, 1, 2, 1, 2, 2, 1, 1, 2]]
    d = pi * [1:14]
end)

dpv = pivot_table(d, ["a", "b"], "c", "d")
@assert( dpv[1,"O"] == d[5,"d"])
@assert( nrow(dpv) == 4 )

dpv2 = pivot_table(d, ["a"], ["c", "b"], "d")
@assert( dpv2[1,"O_J"] == d[5,"d"])

dpv3 = pivot_table(d, ["a"], ["c", "b"], "d", length)
@assert( dpv3[1,"O_J"] == 1.0)

test_group("merge")

srand(1)
df1 = DataFrame(quote
    a = shuffle!([1:10])
    b = ["A","B"][rand(1:2, 10)]
    v1 = randn(10)
end)

df2 = DataFrame(quote
    a = shuffle!(reverse([1:5]))
    b2 = ["A","B","C"][rand(1:3, 5)]
    v2 = randn(3)    # test unequal lengths in the constructor
end)

m1 = merge(df1, df2, "a")
@assert isequal(m1["a"], DataVector[1, 2, 3, 4, 5])
m2 = merge(df1, df2, "a", "outer")
@assert isequal(m2["b2"], DataVector["B", "B", "B", "C", "B", NA, NA, NA, NA, NA])

df1 = DataFrame({"a" => [1, 2, 3],
                 "b" => ["America", "Europe", "Africa"]})
df2 = DataFrame({"a" => [1, 2, 4],
                 "c" => ["New World", "Old World", "New World"]})
m1 = merge(df1, df2, "a", "inner")
@assert isequal(m1["a"], DataVector[1, 2])
m2 = merge(df1, df2, "a", "left")
@assert isequal(m2["a"], DataVector[1, 2, 3])
m3 = merge(df1, df2, "a", "right")
@assert isequal(m3["a"], DataVector[1, 2, 4])
m4 = merge(df1, df2, "a", "outer")
@assert isequal(m4["a"], DataVector[1, 2, 3, 4])

test_group("extras")

srand(1)
x = rand(1:10,10) - 4.0
a1 = cut(x, 4)
a2 = cut(x, [-2, 3, 4.0])
#@assert a2[1] == "[-3.0,-2.0]"
#@assert a2[2] == "(-2.0,3.0]"
#@assert a2[4] == "(4.0,6.0]"

test_group("cut()")
@assert isequal(cut([2, 3, 5], [1, 3, 6]), PooledDataArray(["(1,3]", "(1,3]", "(3,6]"]))
@assert isequal(cut([2, 3, 5], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]"]))
@assert isequal(cut([2, 3, 5, 6], [3, 6]), PooledDataArray(["[2,3]", "[2,3]", "(3,6]", "(3,6]"]))
@assert isequal(cut([1, 2, 4], [1, 3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
@assert isequal(cut([1, 2, 4], [3, 6]), PooledDataArray(["[1,3]", "[1,3]", "(3,6]"]))
@assert isequal(cut([1, 2, 4], [3]), PooledDataArray(["[1,3]", "[1,3]", "(3,4]"]))
@assert isequal(cut([1, 5, 7], [3, 6]), PooledDataArray(["[1,3]", "(3,6]", "(6,7]"]))

ages = [20, 22, 25, 27, 21, 23, 37, 31, 61, 45, 41, 32]
bins = [18, 25, 35, 60, 100]
cats = cut(ages, bins)
pdv = PooledDataArray(["(18,25]", "(18,25]", "(18,25]",
                       "(25,35]", "(18,25]", "(18,25]",
                       "(35,60]", "(25,35]", "(60,100]",
                       "(35,60]", "(35,60]", "(25,35]"])
@assert isequal(cats, pdv)

test_group("New DataVector constructors")
dv = DataArray(Int64, 5)
@assert all(isna(dv))
dv = DataArray(Float64, 5)
@assert all(isna(dv))
dv = datazeros(5)
@assert all(dv .== 0.0)
dv = dataones(5)
@assert all(dv .== 1.0)

# No more NA corruption
dv = dataones(10_000)
@assert !any(isna(dv))

PooledDataArray(falses(2), falses(2))
PooledDataArray(falses(2), trues(2))

# Test vectorized comparisons work for DataVector's and PooledDataVector's
DataVector[1, 2, NA] .== 1
PooledDataVector[1, 2, NA] .== 1
DataVector["1", "2", NA] .== "1"
PooledDataVector["1", "2", NA] .== "1"

# Test unique()
test_group("unique()")
dv = DataArray(1:4)
dv[4] = NA
@assert contains(unique(dv), 1)
@assert contains(unique(dv), 2)
@assert contains(unique(dv), 3)
@assert contains(unique(dv), NA)

test_group("find()")
dv = DataArray([true, false, true])
@assert isequal(find(dv), [1, 3])

pdv = PooledDataArray([true, false, true])
@assert isequal(find(pdv), [1, 3])

dv[1] = NA
@assert isequal(find(dv), [3])

pdv[1] = NA
@assert isequal(find(pdv), [3])


text:
using OrdinaryDiffEq, DiffEqSensitivity, Flux, DiffEqGPU, CUDA, Test
CUDA.allowscalar(false)

function modelf(du,u,p,t)
    du[1] = 1.01 * u[1] * p[1] * p[2]
end

function model()
  prob = ODEProblem(modelf, u0, (0.0, 1.0), pa)

  function prob_func(prob, i, repeat)
    remake(prob, u0 = 0.5 .+ i/100 .* prob.u0)
  end

  ensemble_prob = EnsembleProblem(prob, prob_func = prob_func)
  solve(ensemble_prob, Tsit5(), EnsembleGPUArray(), saveat = 0.1, trajectories = 10)
end

# loss function
loss() = sum(abs2,1.0.-Array(model()))

data = Iterators.repeated((), 10)

cb = function () # callback function to observe training
  @show loss()
end

pa = [1.0,2.0]
u0 = [3.0]
opt = ADAM(0.1)
println("Starting to train")

l1 = loss()

Flux.@epochs 10 Flux.train!(loss, params([pa]), data, opt; cb = cb)

l2 = loss()
@test 3l2 < l1

ext: jl


text:
@testset "AWAP" begin
    using SimpleSDMDataSources: rastername, rasterurl, Values, Deciles,
          SoilMoisture, Lower, Upper, Evapotrans, Potential, Areal, Actual, RefCrop, Tall,
          DeepDrainage

    alwb_path = joinpath(ENV["ECODATASOURCES_PATH"], "ALWB")
    @test rasterpath(ALWB) == alwb_path
    @test rasterpath(ALWB{Values,Year}) == joinpath(alwb_path, "values/year")
    @test rastername(ALWB{Values,Year}, SoilMoisture{Lower}, Date(2001, 1)) == "ss_pct.nc"
    @test rastername(ALWB{Values,Month}, SoilMoisture{Lower}, Date(2001, 1)) == "ss_pct.nc"
    @test rastername(ALWB{Values,Day}, SoilMoisture{Lower}, Date(2001, 1)) == "ss_pct_2001.nc"
    @test rasterpath(ALWB{Values,Day}, SoilMoisture{Lower}, Date(2001, 1)) ==
        joinpath(alwb_path, "values/day/ss_pct_2001.nc")

    @test rasterurl(ALWB{Values,Year}, SoilMoisture{Lower}, Date(2001, 1)) ==
        "http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/year/ss_pct.nc"
    @test rasterurl(ALWB{Values,Day}, SoilMoisture{Lower}, Date(2001, 1)) ==
        "http://www.bom.gov.au/jsp/awra/thredds/fileServer/AWRACMS/values/day/ss_pct_2001.nc"

    dates = DateTime(2018, 01, 01), DateTime(2018, 03, 02)
    download_raster(ALWB{Values,Day}, SoilMoisture{Lower}; dates=dates)
    @test isfile(joinpath(alwb_path, "values/day/ss_pct_2018.nc"))
    download_raster(ALWB{Deciles,Month}, SoilMoisture{Upper}; dates=dates)
    @test isfile(joinpath(alwb_path, "deciles/month/s0_pct.nc"))

    download_raster(ALWB{Values,Day}, Evapotrans{Potential{Areal}}; dates=dates)
    @test isfile(joinpath(alwb_path, "values/day/ma_wet_2018.nc"))
    download_raster(ALWB{Values,Month}, Evapotrans{Actual}; dates=dates)
    @test isfile(joinpath(alwb_path, "values/month/etot.nc"))
    download_raster(ALWB{Values,Year}, Evapotrans{RefCrop{Tall}}; dates=dates)
    @test isfile(joinpath(alwb_path, "values/year/asce_pet.nc"))
    download_raster(ALWB{Values,Year}, DeepDrainage; dates=dates)
    @test isfile(joinpath(alwb_path, "values/year/dd.nc"))
end


text:
using BenchmarkTools

# ## Starting a Cluster
using Distributed
procs()

# change to 2 if started julia with -p2
addprocs(4)

procs()

# addprocs(["10.0.2.1", "10.0.2.2"]) ;

# ## Communication between Julia processes

a = remotecall(sqrt, 2, 4.0)

wait(a)

fetch(a)

remotecall_fetch(sqrt, 2, 4.0)

# # Programming parallel tasks

using Pkg
Pkg.add("Distributions")
using Distributions  # precompile

@everywhere using Distributions

@everywhere println(rand(Normal()))

# ## @spawn macro

a=@spawn randn(5,5)^2

fetch(a)

b=rand(5,5)

a=@spawn b^2

fetch(a)

@time begin
   A = rand(1000,1000)
   Bref = @spawn A^2
   fetch(Bref)
end;

@time begin
   Bref = @spawn rand(1000,1000)^2
   fetch(Bref)
end;

# ## @spawnat

r = remotecall(rand, 2, 2, 2)

s = @spawnat 3 1 .+ fetch(r)

fetch(s)

# ## @parallel for

function serial_add()
    s=0.0
    for i = 1:1000000
         s=s+randn()
    end
    return s
end

function parallel_add()
    return @distributed (+) for i=1:1000000
       randn()
    end
end

@btime serial_add()

@btime parallel_add()

# ## Parallel map

x=[rand(100,100) for i in 1:10];

@everywhere using LinearAlgebra

@btime map(svd, x);

@btime pmap(svd, x);

# ## Distributed Monte Carlo

@everywhere function darts_in_circle(N)
    n = 0
    for i in 1:N
        if rand()^2 + rand()^2 < 1
            n += 1
        end
    end
    return n
end

function pi_distributed(N, loops)
    n = sum(pmap((x)->darts_in_circle(N), 1:loops))
    4 * n / (loops * N)
end

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

@btime pi_distributed(1_000_000, 50)

@btime pi_serial(50_000_000)

# ## Distributed Arrays

using Pkg
Pkg.add("DistributedArrays")
using DistributedArrays
@everywhere using DistributedArrays
d=dzeros(12, 12)
x=rand(10,10);
dx = distribute(x)

@everywhere function par(I)
    d=(size(I[1], 1), size(I[2], 1))
    m = fill(myid(), d)
    return m
end

m = DArray(par, (800, 800))

d.indices

r = @spawnat 2 localpart(d)

fetch(r)

@distributed (+) for i in 1:nworkers()
   sum(localpart(m))
end

# ### Game of Life

function life_step(d::DArray)
   DArray(size(d),procs(d)) do I
      top = mod1(first(I[1])-1,size(d,1))  #outside edge
      bot = mod1( last(I[1])+1,size(d,1))
      left = mod1(first(I[2])-1,size(d,2))
      right = mod1( last(I[2])+1,size(d,2))
      old = Array{Bool}(undef, length(I[1])+2, length(I[2])+2) #temp array
      old[1 , 1 ] = d[top , left]  #get from remote
      old[2:end-1, 1 ] = d[I[1], left] # left
      old[end , 1 ] = d[bot , left]
      old[1 , end ] = d[top , right]
      old[2:end-1, end ] = d[I[1], right] # right
      old[end , end ] = d[bot , right]
      old[1 , 2:end-1] = d[top , I[2]] # top
      old[end , 2:end-1] = d[bot , I[2]] # bottom
      old[2:end-1, 2:end-1] = d[I[1], I[2]] # middle (local)

      life_rule(old) # Step!
   end
end

@everywhere function life_rule(old)
    m, n = size(old)
    new = similar(old, m-2, n-2)
    for j = 2:n-1
        @inbounds for i = 2:m-1
            nc = (+)(old[i-1,j-1], old[i-1,j], old[i-1,j+1],
                     old[i ,j-1], old[i ,j+1],
                     old[i+1,j-1], old[i+1,j], old[i+1,j+1])
            new[i-1,j-1] = (nc == 3 || nc == 2 && old[i,j])
        end
    end
    new
end

A = DArray(I->rand(Bool, length.(I)), (20,20))

using Pkg; Pkg.add("Colors")
using Colors
Gray.(A)

B = copy(A)
B = Gray.(life_step(B))

# # Shared Arrays
using SharedArrays
S=SharedArray{Float64}((100, 100, 5), pids=[2,3, 4, 5]);

pmap(x->S[x]=myid(), eachindex(S));

S

# ## Parallel Prefix sum

function prefix_shared!(y::SharedArray)
    l=length(y)
    k=ceil(Int, log2(l))
    for j=1:k
        @sync @distributed 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
        @sync @distributed 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


text:
@userplot GroupedBar

grouped_xy(x::AbstractVector, y::AbstractMatrix) = x, y
grouped_xy(y::AbstractMatrix) = 1:size(y,1), y

@recipe function f(g::GroupedBar)
    x, y = grouped_xy(g.args...)

    nr, nc = size(y)
    isstack = pop!(d, :bar_position, :dodge) == :stack

    # extract xnums and set default bar width.
    # might need to set xticks as well
    xnums = if eltype(x) <: Number
        bar_width --> (0.8 * mean(diff(x)))
        x
    else
        bar_width --> 0.8
        ux = unique(x)
        xnums = (1:length(ux)) - 0.5
        xticks --> (xnums, ux)
        xnums
    end
    @assert length(xnums) == nr

    # compute the x centers.  for dodge, make a matrix for each column
    x = if isstack
        x
    else
        bws = d[:bar_width] / nc
        bar_width := bws
        xmat = zeros(nr,nc)
        for r=1:nr
            bw = _cycle(bws, r)
            farleft = xnums[r] - 0.5 * (bw * nc)
            for c=1:nc
                xmat[r,c] = farleft + 0.5bw + (c-1)*bw
            end
        end
        xmat
    end

    # compute fillrange
    fillrange := if isstack
        # shift y/fillrange up
        y = copy(y)
        fr = zeros(nr, nc)
        for c=2:nc
            for r=1:nr
                fr[r,c] = y[r,c-1]
                y[r,c] += fr[r,c]
            end
        end
        fr
    else
        get(d, :fillrange, nothing)
    end

    seriestype := :bar
    x, y
end


text:
module TestMOIWrapper

using Test
using MathOptInterface
import Clp

const MOI = MathOptInterface

function runtests()
    for name in names(@__MODULE__; all = true)
        if startswith("$(name)", "test_")
            @testset "$(name)" begin
                getfield(@__MODULE__, name)()
            end
        end
    end
end

function test_SolverName()
    @test MOI.get(Clp.Optimizer(), MOI.SolverName()) == "Clp"
    return
end

function test_supports_default_copy_to()
    @test !MOI.supports_incremental_interface(Clp.Optimizer())
    return
end

function test_runtests()
    model = MOI.Bridges.full_bridge_optimizer(
        MOI.Utilities.CachingOptimizer(
            MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
            Clp.Optimizer(),
        ),
        Float64,
    )
    MOI.set(model, MOI.Silent(), true)
    MOI.Test.runtests(
        model,
        MOI.Test.Config(
            exclude = Any[MOI.DualObjectiveValue, MOI.ObjectiveBound],
        ),
        exclude = [
            # Unable to prove infeasibility
            "test_conic_NormInfinityCone_INFEASIBLE",
            "test_conic_NormOneCone_INFEASIBLE",
        ],
    )
    return
end

function test_Nonexistant_unbounded_ray()
    model = MOI.Utilities.CachingOptimizer(
        MOI.Utilities.UniversalFallback(MOI.Utilities.Model{Float64}()),
        Clp.Optimizer(),
    )
    MOI.set(model, MOI.Silent(), true)
    x = MOI.add_variables(model, 5)
    MOI.set(
        model,
        MOI.ObjectiveFunction{MOI.ScalarAffineFunction{Float64}}(),
        MOI.ScalarAffineFunction(MOI.ScalarAffineTerm.(1.0, x), 0.0),
    )
    MOI.set(model, MOI.ObjectiveSense(), MOI.MAX_SENSE)
    MOI.optimize!(model)
    status = MOI.get(model, MOI.TerminationStatus())
    @test status == MOI.DUAL_INFEASIBLE
    return
end

function test_RawOptimizerAttribute()
    model = Clp.Optimizer()
    MOI.set(model, MOI.RawOptimizerAttribute("LogLevel"), 1)
    @test MOI.get(model, MOI.RawOptimizerAttribute("LogLevel")) == 1
    MOI.set(model, MOI.RawOptimizerAttribute("LogLevel"), 2)
    @test MOI.get(model, MOI.RawOptimizerAttribute("LogLevel")) == 2

    MOI.set(model, MOI.RawOptimizerAttribute("SolveType"), 1)
    @test MOI.get(model, MOI.RawOptimizerAttribute("SolveType")) == 1
    MOI.set(model, MOI.RawOptimizerAttribute("SolveType"), 4)
    @test MOI.get(model, MOI.RawOptimizerAttribute("SolveType")) == 4

    MOI.set(model, MOI.RawOptimizerAttribute("PresolveType"), 1)
    @test MOI.get(model, MOI.RawOptimizerAttribute("PresolveType")) == 1
    MOI.set(model, MOI.RawOptimizerAttribute("PresolveType"), 0)
    @test MOI.get(model, MOI.RawOptimizerAttribute("PresolveType")) == 0
    return
end

function test_All_parameters()
    model = Clp.Optimizer()
    param = MOI.RawOptimizerAttribute("NotAnOption")
    @test !MOI.supports(model, param)
    @test_throws MOI.UnsupportedAttribute(param) MOI.get(model, param)
    @test_throws MOI.UnsupportedAttribute(param) MOI.set(model, param, false)
    for key in Clp.SUPPORTED_PARAMETERS
        @test MOI.supports(model, MOI.RawOptimizerAttribute(key))
        value = MOI.get(model, MOI.RawOptimizerAttribute(key))
        MOI.set(model, MOI.RawOptimizerAttribute(key), value)
        @test MOI.get(model, MOI.RawOptimizerAttribute(key)) == value
    end
    return
end

function test_copy_to_bug()
    model = MOI.Utilities.Model{Float64}()
    x = MOI.add_variable(model)
    con = [
        MOI.add_constraint(
            model,
            MOI.ScalarAffineFunction([MOI.ScalarAffineTerm(1.0, x)], 0.0),
            MOI.EqualTo(1.0),
        ) for i in 1:2
    ]
    clp = Clp.Optimizer()
    index_map = MOI.copy_to(clp, model)
    @test index_map[con[1]] != index_map[con[2]]
    return
end

function test_options_after_empty!()
    model = Clp.Optimizer()
    @test MOI.get(model, MOI.Silent()) == false
    MOI.set(model, MOI.Silent(), true)
    @test MOI.get(model, MOI.Silent()) == true
    MOI.empty!(model)
    @test MOI.get(model, MOI.Silent()) == true
    return
end

end  # module TestMOIWrapper

TestMOIWrapper.runtests()


text:
module SphericalLearning
include("SphericalTensor.jl")
include("SphericalPlotting.jl")
include("HelmholtzExpansion.jl")

end # module


text:
# This file was generated by the Julia Swagger Code Generator
# Do not modify this file directly. Modify the swagger specification instead.

@doc raw"""EnvVarSource represents a source for the value of an EnvVar.

    IoK8sApiCoreV1EnvVarSource(;
        configMapKeyRef=nothing,
        fieldRef=nothing,
        resourceFieldRef=nothing,
        secretKeyRef=nothing,
    )

    - configMapKeyRef::IoK8sApiCoreV1ConfigMapKeySelector : Selects a key of a ConfigMap.
    - fieldRef::IoK8sApiCoreV1ObjectFieldSelector : Selects a field of the pod: supports metadata.name, metadata.namespace, metadata.labels, metadata.annotations, spec.nodeName, spec.serviceAccountName, status.hostIP, status.podIP, status.podIPs.
    - resourceFieldRef::IoK8sApiCoreV1ResourceFieldSelector : Selects a resource of the container: only resources limits and requests (limits.cpu, limits.memory, limits.ephemeral-storage, requests.cpu, requests.memory and requests.ephemeral-storage) are currently supported.
    - secretKeyRef::IoK8sApiCoreV1SecretKeySelector : Selects a key of a secret in the pod's namespace
"""
mutable struct IoK8sApiCoreV1EnvVarSource <: SwaggerModel
    configMapKeyRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ConfigMapKeySelector } # spec name: configMapKeyRef
    fieldRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ObjectFieldSelector } # spec name: fieldRef
    resourceFieldRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1ResourceFieldSelector } # spec name: resourceFieldRef
    secretKeyRef::Any # spec type: Union{ Nothing, IoK8sApiCoreV1SecretKeySelector } # spec name: secretKeyRef

    function IoK8sApiCoreV1EnvVarSource(;configMapKeyRef=nothing, fieldRef=nothing, resourceFieldRef=nothing, secretKeyRef=nothing)
        o = new()
        validate_property(IoK8sApiCoreV1EnvVarSource, Symbol("configMapKeyRef"), configMapKeyRef)
        setfield!(o, Symbol("configMapKeyRef"), configMapKeyRef)
        validate_property(IoK8sApiCoreV1EnvVarSource, Symbol("fieldRef"), fieldRef)
        setfield!(o, Symbol("fieldRef"), fieldRef)
        validate_property(IoK8sApiCoreV1EnvVarSource, Symbol("resourceFieldRef"), resourceFieldRef)
        setfield!(o, Symbol("resourceFieldRef"), resourceFieldRef)
        validate_property(IoK8sApiCoreV1EnvVarSource, Symbol("secretKeyRef"), secretKeyRef)
        setfield!(o, Symbol("secretKeyRef"), secretKeyRef)
        o
    end
end # type IoK8sApiCoreV1EnvVarSource

const _property_map_IoK8sApiCoreV1EnvVarSource = Dict{Symbol,Symbol}(Symbol("configMapKeyRef")=>Symbol("configMapKeyRef"), Symbol("fieldRef")=>Symbol("fieldRef"), Symbol("resourceFieldRef")=>Symbol("resourceFieldRef"), Symbol("secretKeyRef")=>Symbol("secretKeyRef"))
const _property_types_IoK8sApiCoreV1EnvVarSource = Dict{Symbol,String}(Symbol("configMapKeyRef")=>"IoK8sApiCoreV1ConfigMapKeySelector", Symbol("fieldRef")=>"IoK8sApiCoreV1ObjectFieldSelector", Symbol("resourceFieldRef")=>"IoK8sApiCoreV1ResourceFieldSelector", Symbol("secretKeyRef")=>"IoK8sApiCoreV1SecretKeySelector")
Base.propertynames(::Type{ IoK8sApiCoreV1EnvVarSource }) = collect(keys(_property_map_IoK8sApiCoreV1EnvVarSource))
Swagger.property_type(::Type{ IoK8sApiCoreV1EnvVarSource }, name::Symbol) = Union{Nothing,eval(Base.Meta.parse(_property_types_IoK8sApiCoreV1EnvVarSource[name]))}
Swagger.field_name(::Type{ IoK8sApiCoreV1EnvVarSource }, property_name::Symbol) =  _property_map_IoK8sApiCoreV1EnvVarSource[property_name]

function check_required(o::IoK8sApiCoreV1EnvVarSource)
    true
end

function validate_property(::Type{ IoK8sApiCoreV1EnvVarSource }, name::Symbol, val)
end


text:
#
# Vector with values stored in a sorted dictionary.
#

immutable DictVector{T} <: AbstractVector{T}
    ptr::Vector{Int}
    dict::Vector{T}
end

function DictVector{T}(vals::AbstractVector{T})
    len = length(vals)
    v2i = Dict{T,Int}()
    ptr = Vector{Int}(len)
    for k = 1:len
        ptr[k] = get!(v2i, vals[k], k)
    end
    ulen = length(v2i)
    vs = collect(keys(v2i))
    is = collect(values(v2i))
    perm = sortperm(vs)
    recode = Vector{Int}(ulen)
    for k = 1:ulen
        recode[is[k]] = perm[k]
    end
    for k = 1:len
        ptr[k] = recode[ptr[k]]
    end
    dict = vs[perm]
    return DictVector{T}(ptr, dict)
end

keys(dv::DictVector) = dv.ptr
values(dv::DictVector) = dv.dict

size(dv::DictVector) = (length(dv.ptr),)
length(dv::DictVector) = length(dv.ptr)

getindex(dv::DictVector, i::Int) = dv.dict[dv.ptr[i]]

getindex(dv::DictVector, idxs::AbstractVector{Int}) =
    DictVector(dv.ptr[idxs], dv.dict)

max_stars_repo_name: rbt-lang/rbt-proto

max_issues_repo_name: rbt-lang/rbt-proto

max_forks_repo_name: rbt-lang/rbt-proto


text:
using SBMLBioModelsRepository
using Pkg, Test
Pkg.add(url="https://github.com/paulflang/SBML.jl/", rev="pl/mk-species-units")
using SBML
using ModelingToolkit, OrdinaryDiffEq, CSV, DataFrames
using Base.Threads, Glob

!isdir("logs/") && mkdir("logs/")

@testset "SBMLBioModelsRepository.jl" begin
    @testset "biomd" begin include("biomd.jl") end
    # @testset "test_suite" begin include("test_suite.jl") end
end


text:
using SimpleSparseArrays
using SimpleSparseArrays.Testing
using LinearAlgebra
using Test

@testset "basics" begin
    @test SimpleSparseArray(1,2,3) == SimpleSparseArray{Float64}(1,2,3)
end

@testset "no redundant zeros" begin
    arr = SimpleSparseArray{Float64}(10,10)
    @test isempty(eachstoredindex(arr))
    @test isempty(eachstoredindex(arr + arr))
    @test isempty(eachstoredindex(5 * arr))
    @test isempty(eachstoredindex(map(sin, arr)))
    @test isempty(eachstoredindex(map(+, arr, arr)))
    @test isempty(eachstoredindex(map(+, arr, arr, arr, arr)))
    @test isempty(eachstoredindex(sparsify(zeros(3))))

    @test_broken isempty(eachstoredindex(arr[:,1]))
    @test_broken isempty(eachstoredindex(arr[1:3, 2:4]))
    @test_broken isempty(eachstoredindex(sqrt.(arr)))
    @test_broken isempty(eachstoredindex(sin.(arr) .+ Ref(1) .* arr))
end

@testset "test against full" begin
    v1 = randn(3)
    v2 = randn(3)
    vs1  = sparsify(v1)
    vs2  = sparsify(v2)
    @test_broken vs1 == v1
    @test vs1 == vs1
    @test vs1 == deepcopy(vs1)
    @test vs1 ≈ v1
    @test !(vs1 ≈ vs2)
    @test !(v1 ≈ vs2)
    for i in eachindex(v1)
        @test v1[i] === vs1[i]
    end

    @test isapprox_full(+, vs1, vs2)
    @test isapprox_full(-, vs1, vs2)
    @test isapprox_full(norm, vs1)
    @test isapprox_full(sum, vs1)
    @test isapprox_full(map, +, vs1, vs2, vs1)
    @test isapprox_full(map, -, vs1, vs2)
    @test isapprox_full(map, sin, vs1)
    @test isapprox_full(map, cos, vs1)

end


text:
using Stheno, Random, Plots, Statistics
using Stheno: @model, EQ, Noise

###########################  Define and inspect our model  ###########################

rng = MersenneTwister(123456);

#=
In this example, `f` is an unknown real-valued function that we wish to infer. To achieve
this, we have access to two sensors. The first returns noisy estimates of `f`, where we have
been reliably informed by whoever designed the sensor that the mean of the noise is given by
`sin(x) - 5 + sqrt(abs(x))`, and that it's variance is low (1e-2). (how the designer
estimated this function, is why a sensor might possibly have such a weird mean error is
beyond the scope of this example) The second returns biased measurements of `f`, where the
bias is known to be 3.5. The model below specifies a model for this scenario.
=#
@model function model()

    # Define a smooth latent process that we wish to infer.
    f = GP(EQ())

    # Define the two noise processes described.
    noise1 = sqrt(1e-2) * GP(Noise()) + (x->sin.(x) .- 5.0 .+ sqrt.(abs.(x)))
    noise2 = sqrt(1e-1) * GP(3.5, Noise())

    # Define the processes that we get to observe.
    y1 = f + noise1
    y2 = f + noise2

    return f, noise1, noise2, y1, y2
end
f, noise₁, noise₂, y₁, y₂ = model();

# Generate some toy observations of `y1` and `y2`.
X₁, X₂ = sort(rand(rng, 3) * 10), sort(rand(rng, 10) * 10);
ŷ₁, ŷ₂ = rand(rng, [y₁(X₁), y₂(X₂)]);

# Compute the posterior processes.
(f′, y₁′, y₂′) = (f, y₁, y₂) | (y₁(X₁)←ŷ₁, y₂(X₂)←ŷ₂);

# Sample jointly from the posterior processes and compute posterior marginals.
Xp = range(-2.5, stop=12.5, length=500);
f′Xp, y₁′Xp, y₂′Xp = rand(rng, [f′(Xp, 1e-9), y₁′(Xp, 1e-9), y₂′(Xp, 1e-9)], 100);

ms1 = marginals(f′(Xp));
ms2 = marginals(y₁′(Xp));
ms3 = marginals(y₂′(Xp));

μf′, σf′ = mean.(ms1), std.(ms1);
μy₁′, σy₁′ = mean.(ms2), std.(ms2);
μy₂′, σy₂′ = mean.(ms3), std.(ms3);

###########################  Plot results  ###########################

plotly();

posterior_plot = plot();

# Plot posterior marginal std. dev.
plot!(posterior_plot, Xp, [μy₁′ μy₁′];
    linewidth=0.0,
    fillrange=[μy₁′ .- 3 .* σy₁′, μy₁′ .+ 3 * σy₁′],
    fillalpha=0.3,
    fillcolor=:red,
    label="");
plot!(posterior_plot, Xp, [μy₂′ μy₂′];
    linewidth=0.0,
    fillrange=[μy₂′ .- 3 .* σy₂′, μy₂′ .+ 3 * σy₂′],
    fillalpha=0.3,
    fillcolor=:green,
    label="");
plot!(posterior_plot, Xp, [μf′ μf′];
    linewidth=0.0,
    fillrange=[μf′.- 3  .* σf′ μf′ .+ 3 .* σf′],
    fillalpha=0.5,
    fillcolor=:blue,
    label="");

# Plot posterior marginal samples.
scatter!(posterior_plot, Xp, y₁′Xp,
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=0.5,
    markeralpha=0.3,
    label="");
scatter!(posterior_plot, Xp, y₂′Xp,
    markercolor=:green,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=0.5,
    markeralpha=0.3,
    label="");
plot!(posterior_plot, Xp, f′Xp;
    linecolor=:blue,
    linealpha=0.2,
    label="");

# Plot posterior means
plot!(posterior_plot, Xp, μy₁′;
    linecolor=:red,
    linewidth=2.0,
    label="");
plot!(posterior_plot, Xp, μy₂′;
    linecolor=:green,
    linewidth=2.0,
    label="");
plot!(posterior_plot, Xp, μf′;
    linecolor=:blue,
    linewidth=2.0,
    label="Latent Function");

# Plot posterior of first noise process.
scatter!(posterior_plot, X₁, ŷ₁;
    markercolor=:red,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Sensor 1");
scatter!(posterior_plot, X₂, ŷ₂;
    markercolor=:green,
    markershape=:circle,
    markerstrokewidth=0.0,
    markersize=4,
    markeralpha=0.8,
    label="Sensor 2");

display(posterior_plot);