text:
include("SortSparse.jl")
include("KoLesky.jl")
include("Cholesky.jl")
include("CovFuncs.jl")
include("Utils.jl")
using LinearAlgebra
using Distributions
using Distances
using Random
using Match
using JLD
using BlockArrays
using ArgParse
Random.seed!(123)

s = ArgParseSettings()
@add_arg_table s begin
  "maternOrder"
    help = "The order of the matern kernel to be used. Should be in (\"12\", \"32\", \"52\")"
    required = true
end
maternOrder = parse_args(s)["maternOrder"]

NTrain = 10000
NTest = 1
NSamples = 100
N = NTrain + NTest
d = 2
l = 0.5
# For now, only properly implemented for $σ = 0$.
σ = 0.0 * exp.(randn(N))
λ = 1.5
ρ = 3.0

# Colums
meanErrorPredFirst = zeros(NTest, NSamples)
stdErrorPredFirst = zeros(NTest, NSamples)
KLTrueApproxPredFirst = zeros(NSamples)
KLApproxTruePredFirst = zeros(NSamples)
meanErrorPredLast = zeros(NTest, NSamples)
stdErrorPredLast = zeros(NTest, NSamples)
KLTrueApproxPredLast = zeros(NSamples)
KLApproxTruePredLast = zeros(NSamples)
meanErrorNoPred = zeros(NTest, NSamples)
stdErrorNoPred = zeros(NTest, NSamples)
KLTrueApproxNoPred = zeros(NSamples)
KLApproxTrueNoPred = zeros(NSamples)

# Creating data on ellipsoid, testing data in the center
h = 1 / NTrain
xTrain = h : h : 1
xTrain = vcat(sin.(2 * π * xTrain)',
              2 * cos.(2 * π * xTrain)')
xTest = 0.1 * randn(2,NTest)
x = hcat(xTrain, xTest)

σ .= max.(σ, eps(Float64))
σ = PseudoBlockArray(σ, [NTest, NTrain])

covfunc12!(r) = matern12!(r,l)
covfunc32!(r) = matern32!(r,l)
covfunc52!(r) = matern52!(r,l)

# Make sure to avoid NANs
@match maternOrder begin
  "12" => (covfunc! = covfunc12!)
  "32" => (covfunc! = covfunc32!)
  "52" => (covfunc! = covfunc52!)
end

for k = 1 : NSamples
  # =====================================
  # Starting with interpolation:
  # =====================================

  # Creating reference solution
  Θ = covfunc!(pairwise(Euclidean(), x, dims=2)) + Diagonal(σ)
  L_Θ = cholesky(Matrix(Θ))
  y = L_Θ.L * randn(N)
  Θ = Θ + Diagonal(σ)

  # Creating data
  Θ = PseudoBlockArray(Θ, [NTrain, NTest], [NTrain, NTest])
  y = PseudoBlockArray(y, [NTrain, NTest])
  y[Block(1)] .+= Diagonal(sqrt.(σ[Block(2)])) * randn(NTrain)

  # Computing the true posterior mean
  μ = Θ[Block(2,1)] * (Θ[Block(1,1)] \ y[Block(1)])
  Σ = Θ[Block(2,2)] - Θ[Block(2,1)] * (Θ[Block(1,1)] \ Θ[Block(1,2)])

  # =====================================
  # Computing posterior mean for prediction variables last
  # =====================================

  # reordering points and forming skeletons
  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, xTest, ρ, 1)
  xOrd = x[:, P]
  xOrdTest = xOrd[:, 1 : NTest]
  xOrdTrain = xOrd[:, (NTest + 1) : end]
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  μOrd = μ[P[1 : NTest] .- NTrain]
  ΘOrd = Θ[P, P]
  ΣOrd = Σ[P[1 : NTest] .- NTrain, P[1 : NTest] .- NTrain]
  yOrd = PseudoBlockArray(y[P], [NTest, NTrain])
  σOrd = PseudoBlockArray(σ[P], [NTest, NTrain])

  # Creating vecchia approximation
  L = assembleL(skeletons, xOrd, covfunc!)
  L22 = L[(NTest + 1) : end, (NTest + 1) : end]
  U22 = squareSparse(L22) + spdiagm(0 => 1 ./ σOrd[Block(2)])
  icholU_high_level!(U22)
  noisecov = NoiseCov(L22, U22, σOrd[Block(2)])
  μPredLast = (L' \ ( L \ vcat(zero(y[Block(2)]), noisecov \ (yOrd[Block(2)]))))[1:NTest]

  ΣPredLast = inv(Matrix(L[1 : NTest, 1 : NTest] * L[1 : NTest, 1 : NTest]'))

  ΘPredLast = inv(Matrix(L * L'))

  meanErrorPredLast[:, k] = abs.(μOrd - μPredLast)
  stdErrorPredLast[:, k] = abs.(diag(ΣOrd) - diag(ΣPredLast))
  KLTrueApproxPredLast[k] = KL(ΘOrd, ΘPredLast)
  KLApproxTruePredLast[k] = KL(ΘPredLast, ΘOrd)

  # =====================================
  # Computing posterior mean for prediction variables first
  # =====================================

  # reordering points and forming skeletons
  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
  xOrdTrain = xTrain[:, P]
  xOrd = hcat(xOrdTrain, xTest)
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  μOrd = μ
  ΘOrd = Θ[vcat(P, (NTrain + 1) : N), vcat(P, (NTrain + 1) : N)]
  ΣOrd = Σ
  yOrd = y[Block(1)][P]
  σOrd = σ[P]

  μPredFirst, ΣPredFirst = predict(skeletons, (NTrain+1): N, yOrd, xOrd, covfunc!, σ)

  # Assembling the matrix implicit in predict, in order to compute KL divergence
  L = SparseMatrixCSC(NTrain, NTrain, colptr, rowval, ones(size(rowval, 1)))
  L = vcat(L, ones(NTest, NTrain))
  L = tril(hcat(L, ones(N, NTest)))
  colptr = L.colptr
  rowval = L.rowval
  distances = vcat(distances, fill(Inf, NTest))
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  L = assembleL(skeletons, xOrd, covfunc!)
  ΘPredFirst = inv(Matrix(L * L'))

  meanErrorPredFirst[:, k] = abs.(μOrd - μPredFirst)
  stdErrorPredFirst[:, k] = abs.(diag(ΣOrd) - ΣPredFirst)
  KLTrueApproxPredFirst[k] = KL(ΘOrd, ΘPredFirst)
  KLApproxTruePredFirst[k] = KL(ΘPredFirst, ΘOrd)

  # =====================================
  # Don't include prediction variables at all
  # =====================================

  # reordering points and forming skeletons
  colptr, rowval, P, revP, distances = sortSparseRev(xTrain, ρ, 1)
  xOrdTrain = xTrain[:, P]
  xOrd = hcat(xOrdTrain, xTest)
  skeletons = construct_skeletons(colptr, rowval, distances, λ)
  L = assembleL(skeletons, xOrdTrain, covfunc!, σ)
  μOrd = μ
  ΘOrd = Θ[vcat(P, (NTrain + 1) : N), vcat(P, (NTrain + 1) : N)]
  ΣOrd = Σ
  yOrd = y[Block(1)][P]
  σOrd = σ[P]

  TMat = pairwise(Euclidean(), xOrdTrain, xTest, dims=2)
  covfunc!(TMat)
  CMat = pairwise(Euclidean(), xTest, dims=2)
  covfunc!(CMat)

  μNoPred = vec(yOrd' * L * L' * TMat)
  ΣNoPred = CMat - TMat' * L * L' * TMat

  ΘNoPred = inv(Matrix(L * L'))
  ΘNoPred = vcat(ΘNoPred, TMat')
  ΘNoPred = hcat(ΘNoPred, vcat(TMat, CMat))

  meanErrorNoPred[:, k] = abs.(μOrd - μNoPred)
  stdErrorNoPred[:, k] = abs.(ΣOrd - ΣNoPred)
  KLTrueApproxNoPred[k] = KL(ΘOrd, ΘNoPred)
  KLApproxTrueNoPred[k] = KL(ΘNoPred, ΘOrd)
end

save("./out/jld/IncludePredCircle$maternOrder.jld",
     "meanErrorNoPred", meanErrorNoPred,
     "stdErrorNoPred", stdErrorNoPred,
     "KLTrueApproxNoPred", KLTrueApproxNoPred,
     "KLApproxTrueNoPred", KLApproxTrueNoPred,
     "meanErrorPredFirst", meanErrorPredFirst,
     "stdErrorPredFirst", stdErrorPredFirst,
     "KLTrueApproxPredFirst", KLTrueApproxPredFirst,
     "KLApproxTruePredFirst", KLApproxTruePredFirst,
     "meanErrorPredLast", meanErrorPredLast,
     "stdErrorPredLast", stdErrorPredLast,
     "KLTrueApproxPredLast", KLTrueApproxPredLast,
     "KLApproxTruePredLast", KLApproxTruePredLast)


text:
# Learn DataStructures.jl
# https://github.com/JuliaCollections/DataStructures.jl

# TODO: Combine learn_SortedDict.jl, learn_SortedSet.jl into here

using Test
#import DataStructures
using DataStructures
using Dates

"""
Example custom struct for testing with containers.

"""
struct Node
    id::Int64
end
# https://discourse.julialang.org/t/proper-way-to-overload-operators/19872
Base.:(<)(x::Node, y::Node) = x.id < y.id
Base.:(==)(x::Node,y::Node) = x.id == y.id

"""
Demonstrate `id` alone identify the uniqueness of a node and `value` is used
to compare magnitude.
"""
struct Node2
    id::Int64
    value::Float64
end
Base.:(<)(x::Node2, y::Node2) = x.value < y.value
Base.:(==)(x::Node2,y::Node2) = x.value == y.value
# isequal is used for ID
Base.isequal(x::Node2,y::Node2) = isequal(x.id, y.id) #error("Unsupported operation") #x.id == y.id
Base.hash(x::Node2, h::UInt64=UInt64(13)) = hash(x.id,h)

function learn_Node2()
    # N{id,value}
    N = Node2 # save some typing

    # (<)
    @test N(3, 3.3) < N(2, 4.4)
    @test N(3, 3.3) <= N(2, 3.3)
    @test N(2, 4.4) > N(3, 3.3)
    @test N(2, 4.4) >= N(3, 4.4)

    # (==)
    @test N(3, 3.3) != N(4, 4.4) # value !=
    @test N(3, 3.3) == N(4, 3.3) # value ==

    # hash(), isequal uses ID for uniqueness
    db = Dict{N,Int64}()
    @test isempty(db)
    db[N(3, 3.3)] = 3
    @test 1 == length(db)
    db[N(3, 4.4)] = 4     # same key as above
    @test 1 == length(db) # so still 1 element
    @test 4 == db[N(3, 5.5)]
    db[N(4, 6.6)] = 6     # add one element
    @test 2 == length(db)

    # (===)
    @test N(3, 3.3) === N(3, 3.3) # immutables are compared at bit level

    nothing
end

function learn_Dict()
    db = Dict{Node2,Int64}()
    # @show db

    db[Node2(1,2)] = 13
    db[Node2(1,3)] = 14
    @test 14 == db[Node2(1,2)]
    #@show getindex(db, Node2(1,2))

    @show db
end

function learn_OrderedSet()
    set = OrderedSet{Int64}()
    push!(set, 13)
    push!(set, 7)
    push!(set, 17)

    @test [13,7,17] == collect(set)
end

function learn_PriorityQueue()
    # ------------------ #
    # Value: Low -> High #
    # ------------------ #
    que = PriorityQueue{String,Int64}() # (key,value), value: low -> high
    @test isempty(que)

    enqueue!(que, "z", 100) # z,100
    @test !isempty(que)
    @test 1 == length(que)

    enqueue!(que, "w", 300) # w,300
    @test 2 == length(que)

    # Get key only
    x = dequeue!(que)
    @test 1 == length(que)
    @test "z" == x

    # Get (key,value)
    x = dequeue_pair!(que)
    @test isempty(que)
    @test "w" == x[1]
    @test 300 == x[2]

    # ------------------ #
    # Value: High -> Low #
    # ------------------ #
    que = PriorityQueue{String,Int64}(Base.Order.Reverse) # (key,value), value: high -> low
    @test isempty(que)
    enqueue!(que, "z", 100) # z,100
    @test !isempty(que)
    @test 1 == length(que)

    enqueue!(que, "p", 300) # p,300
    @test 2 == length(que)

    enqueue!(que, "w", 700) # w,700
    @test 3 == length(que)

    # Get key only
    x = dequeue!(que) # w,700
    @test 2 == length(que)
    @test "w" == x

    # Get (key,value)
    x = dequeue_pair!(que) # p,300
    @test 1 == length(que)
    @test "p" == x[1]
    @test 300 == x[2]

    # TODO
    # ---------
    # Use custom type as key::Node #
    # ---------------------------- #
    que = PriorityQueue{Node2,Float64}()
    enqueue!(que, Node2(1,2), 3.0)
    enqueue!(que, Node2(1,3), 4.0)
end

"""
    learn_Queue()

Use `methodswith(Queue)` to see methods for `Queue`.
"""
function learn_Queue()
    #ds = DataStructures

    que = Queue{Int64}()
    @test isempty(que)

    enqueue!(que, 7)
    @test !isempty(que)
    @test 1 == length(que)
    @test 7 == first(que)
    @test 7 == last(que)

    enqueue!(que, 13)
    @test 2 == length(que)
    @test 7 == first(que)
    @test 13 == last(que)

    x = dequeue!(que)
    @test 7 == x
    @test 1 == length(que)
    @test 13 == first(que)
    @test 13 == last(que)
end

"""
    learn_Set()

Learn how to use `Set`.  Common Set functions

```julia
julia> methodswith(Set)
[1] length(s::Set) in Base at set.jl:55
[2] sizehint!(s::Set, newsz) in Base at set.jl:74
[3] copy(s::Set) in Base at set.jl:68
[4] delete!(s::Set, x) in Base at set.jl:66
[5] empty!(s::Set) in Base at set.jl:75
[6] filter!(f, s::Set) in Base at set.jl:390
[7] in(x, s::Set) in Base at set.jl:56
[8] isempty(s::Set) in Base at set.jl:54
[9] iterate(s::Set, i...) in Base at set.jl:78
[10] pop!(s::Set) in Base at set.jl:61
[11] pop!(s::Set, x) in Base at set.jl:58
[12] pop!(s::Set, x, default) in Base at set.jl:59
[13] push!(s::Set, x) in Base at set.jl:57
[14] setdiff!(s::Set, t::Set) in Base at set.jl:87
[15] show(io::IO, s::Set) in Base at set.jl:39
```
"""
function learn_Set()
    # push!, pop!
    set = Set{Int64}()

    push!(set, 1)
    push!(set, 13)

    @test in(1, set)
    @test in(13, set)
    @test !in(2, set)

    x = pop!(set, 1)
    @test 1 == x
    @test !in(1, set)

    # empty!
    set = Set{Int64}()
    @test isempty(set)
    push!(set, 7)
    @test !isempty(set)
    x = empty!(set)
    @test x === set
    @test isempty(set)

    # Node
    set = Set{Node}()

    push!(set, Node(1))
    push!(set, Node(13))
    @test in(Node(1), set)
    @test in(Node(13), set)

    x = pop!(set, Node(13))
    @test Node(13) == x
    @test !in(Node(13), set)

end

function learn_SortedSet()
    #>> Sort by Date <<#
    d = SortedSet{Date}(Base.Order.Forward)
    push!(d, Date(2020, 12, 24))
    push!(d, Date(2020, 12, 22))
    push!(d, Date(2020, 12, 23))

    #@info d
    @test [Date("2020-12-22"), Date("2020-12-23"), Date("2020-12-24")] == collect(d)

    # Check for key
    @test in(Date(2020, 12, 24), d)
    @test !in(Date(2019, 1, 1), d)

    set = SortedSet{Int64}([8,2,3,1,7], Base.Order.Forward)
    @test [1,2,3,7,8] == collect(set)
    set = SortedSet{Int64}(Base.Order.Forward, [8,2,3,1,7])
    @test [1,2,3,7,8] == collect(set)
end

function learn_Stack()
    ds = DataStructures

    stack = ds.Stack{Int64}()
    push!(stack, 1)
    push!(stack, 2)

    # first in last out
    @test 2 == pop!(stack)
    @test 1 == pop!(stack)
end

@testset "Node2" begin
    learn_Node2()
end

@testset "Set" begin
    learn_OrderedSet()
    learn_SortedSet()
end
#learn_Dict()
#A = learn_PriorityQueue()
#A = learn_Queue()
#A = learn_Set()
#A = learn_Stack()
#learn_OrderedSet()
#A
nothing

ext: jl


text:
using DifferentialEquations, Plots; pyplot()

beta, delta, gamma = 0.25, 0.2, 0.1
initialInfect = 0.025
println("R0 = ", beta/gamma)

initX = [1-initialInfect, 0.0, initialInfect, 0.0]
tEnd = 100.0

RHS(x,parms,t) = [  -beta*x[1]*x[3],
                    beta*x[1]*x[3] - delta*x[2],
                    delta*x[2] - gamma*x[3],
                    gamma*x[3] ]

prob = ODEProblem(RHS, initX, (0,tEnd), 0)
sol = solve(prob)
println("Final infected proportion= ", sol.u[end][4])

plot(sol.t,((x)->x[1]).(sol.u),label = "Susceptible", c=:green)
plot!(sol.t,((x)->x[2]).(sol.u),label = "Exposed", c=:blue)
plot!(sol.t,((x)->x[3]).(sol.u),label = "Infected", c=:red)
plot!(sol.t,((x)->x[4]).(sol.u),label = "Removed", c=:yellow,
    xlabel = "Time", ylabel = "Proportion",legend = :top)

max_stars_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_issues_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook

max_forks_repo_name: Yoshinobu-Ishizaki/StatsWithJuliaBook


text:
using Documenter

makedocs(sitename="BTG Docs")

# Refer to https://juliadocs.github.io/Documenter.jl/stable/man/guide/#Package-Guide-1
# and https://juliadocs.github.io/Documenter.jl/stable/man/hosting/#Hosting-Documentation-1
# for details on hosting and building a doc site

max_stars_repo_name: dbindel/btg

max_issues_repo_name: dbindel/btg

max_forks_repo_name: dbindel/btg


text:
module Leaf

include("stomcond.jl")
include("ebal.jl")
# include("photosyn.jl")

# export submodule functions to the current level
using .StomCond
using .EBal

end  # module


text:
# Print "Hello, World".

println("Hello, World")


text:
# 26 in base 16 or 2
base(16, 26)
base(2, 26)

# Parse to integer
parse(Int, "1a", 16)
parse(Int, "101101", 2)

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
#largeeig.png
#Algorithm 1.4 of Random Eigenvalues by Alan Edelman

#Experiment:  Largest eigenvalue of random Hermitian matrices
#Plot:        Histogram of the normalized largest eigenvalues
#Theory:      Tracy-Widom as n->infinity

include("tracywidom.jl")

## Parameters
n=100        # matrix size
t=5000       # trials
dx=.2        # binsize
function largeeig_experiment(n,t,dx)
    ## Experiment
    v=Float64[] # eigenvalue samples
    for i=1:t
        a=randn(n,n)+im*randn(n,n) # random nxn complex matrix
        s=(a+a')/2                 # Hermitian matrix
        push!(v,max(eig(s)[1]))    # Largest Eigenvalue
    end
    v=n^(1/6)*(v-2*sqrt(n))        # normalized eigenvalues

    grid=-5:dx:2
    count=hist(v,grid)[2]/(t*dx)
    (grid, count)
end
(grid, count) = largeeig_experiment(n,t,dx)

## Theory
(t, f) = tracywidom(5., -8.)

## Plot
using Winston
p = FramedPlot()
h = Histogram(count, step(grid))
h.x0 = first(grid)
add(p, h)
add(p, Curve(t, f, "linewidth", 2, "color", "blue"))
if isinteractive()
    Winston.display(p)
else
    file(p, "largeeig.png")
end


text:
# Salsa20 implementation in Julia
# Copyright (C): Jaka Smrekar (vinctux) <vinctux@outlook.com>, 2016
#
# Reference:       https://cr.yp.to/snuffle/salsafamily-20071225.pdf
# x86 C reference: https://cr.yp.to/snuffle/salsa20/ref/salsa20.c
#                  https://github.com/keybase/python-salsa20/blob/master/libsodium-salsa20

using Krypto

global const σ = Array{UInt8, 1}("expand 32-byte k")

function LLE(x::Array{UInt8, 1}, l::Integer)
    a = 0
    for i in 0:3 a |= (Int32(x[l + i + 1]) << 8i) end
    return a
end
SLE(u::Int32) = [(u & [0x000000FF, 0x0000FF00, 0x00FF0000, 0xFF000000][i]) >> 8i for i in 1:4]

# Salsa20 cryptographic primitive
function SALSA20(I::Array{UInt8, 1}, k::Array{UInt8, 1}, c::Array{UInt8, 1}, R::Integer = 20)
    X = J = Int32[LLE(c, 0),  LLE(k, 0),  LLE(k, 4),  LLE(k, 8),
                  LLE(k, 12), LLE(c, 4),  LLE(I, 0),  LLE(I, 4),
                  LLE(I, 8),  LLE(I, 12), LLE(c, 8),  LLE(k, 16),
                  LLE(k, 20), LLE(k, 24), LLE(k, 28), LLE(c, 12)]
    for i in 1:2:R
        X[5]  ⊻= rotl( X[1] + X[13],  7)
        X[9]  ⊻= rotl( X[5] +  X[1],  9)
        X[13] ⊻= rotl( X[9] +  X[5], 13)
        X[1]  ⊻= rotl(X[13] +  X[9], 18)
        X[10] ⊻= rotl( X[6] +  X[2],  7)
        X[14] ⊻= rotl(X[10] +  X[6],  9)
        X[2]  ⊻= rotl(X[14] + X[10], 13)
        X[6]  ⊻= rotl( X[2] + X[14], 18)
        X[15] ⊻= rotl(X[11] +  X[7],  7)
        X[3]  ⊻= rotl(X[15] + X[11],  9)
        X[7]  ⊻= rotl( X[3] + X[15], 13)
        X[11] ⊻= rotl( X[7] +  X[3], 18)
        X[4]  ⊻= rotl(X[16] + X[12],  7)
        X[8]  ⊻= rotl( X[4] + X[16],  9)
        X[12] ⊻= rotl( X[8] +  X[4], 13)
        X[16] ⊻= rotl(X[12] +  X[8], 18)
        X[2]  ⊻= rotl( X[1] +  X[4],  7)
        X[3]  ⊻= rotl( X[2] +  X[1],  9)
        X[4]  ⊻= rotl( X[3] +  X[2], 13)
        X[1]  ⊻= rotl( X[4] +  X[3], 18)
        X[7]  ⊻= rotl( X[6] +  X[5],  7)
        X[8]  ⊻= rotl( X[7] +  X[6],  9)
        X[5]  ⊻= rotl( X[8] +  X[7], 13)
        X[6]  ⊻= rotl( X[5] +  X[8], 18)
        X[12] ⊻= rotl(X[11] + X[10],  7)
        X[9]  ⊻= rotl(X[12] + X[11],  9)
        X[10] ⊻= rotl( X[9] + X[12], 13)
        X[11] ⊻= rotl(X[10] +  X[9], 18)
        X[13] ⊻= rotl(X[16] + X[15],  7)
        X[14] ⊻= rotl(X[13] + X[16],  9)
        X[15] ⊻= rotl(X[14] + X[13], 13)
        X[16] ⊻= rotl(X[15] + X[14], 18)
    end
    for i in 1:16 X[i] += J[i] end
    return flat([SLE(x) for x in X])
end

# Salsa20 stream expander primitive
function Salsa20Stream(l::Integer, n::Array{UInt8, 1}, k::Array{UInt8, 1}, R::Integer = 20)
    assert(l > 0); assert(length(n) >= 8); assert(length(k) >= 32)
    I, C, q = vcat(n[1:8], zeros(UInt8, 8)), zeros(UInt8, l), l % 64
    for i in 0:div(l, 64)
        C[(64i + 1):64(i + 1)], u = SALSA20(I, k, σ, R), 1
        for j in 9:16
            I[j] = (u += UInt8(I[j]))
            u >>= 8
        end
    end
    if q > 0 C[1:q] = SALSA20(I, k, σ, R)[1:q] end
    return C
end

# Salsa20 xored stream expander primitive
function Salsa20StreamXOR(m::Array{UInt8, 1}, l::Integer, n::Array{UInt8, 1}, k::Array{UInt8, 1}, R::Integer = 20)
    assert(l > 0); assert(length(n) >= 8); assert(length(k) >= 32)
    I, C, q = vcat(n[1:8], zeros(UInt8, 8)), zeros(UInt8, l), l % 64
    for i in 0:(div(l, 64) - 1)
        B, u, t = SALSA20(I, k, σ, R), 1, (64i + 1):64(i + 1)
        C[t] = m[t] ⊻ B
        for j in 9:16
            I[j] = (u += UInt8(I[j]))
            u >>= 8
        end
    end
    if q > 0 C[1:q] = m[1:q] ⊻ SALSA20(I, k, σ, R)[1:q] end
    return C
end

encrypt(::Type{Krypto.Salsa20}, M::Array{UInt8, 1},
        IV::Array{UInt8, 1}, K::Array{UInt8, 1},
        R::Integer = 20) = Salsa20StreamXOR(M, length(M), IV, K, R)
decrypt(::Type{Krypto.Salsa20}, C::Array{UInt8, 1},
        IV::Array{UInt8, 1}, K::Array{UInt8, 1},
        R::Integer = 20) = encrypt(Krypto.Salsa20, C, IV, K, R)


text:
module Commands

using Discord

function echo(c::Client, msg::Message)
    content = lstrip(msg.content[6:end])  # Trim off the command part.
    reply(c, msg, content)
end

function main()
    c = Client(ENV["DISCORD_TOKEN"])
    add_command!(c, "!echo", echo)
    open(c)
    return c
end

end

if abspath(PROGRAM_FILE) == @__FILE__
    c = Commands.main()
    wait(c)
end


text:
# set up packages:
using Pkg
Pkg.activate(".")
Pkg.add("RCall")
#Pkg.add(PackageSpec(url="https://github.com/esa-esdl/ESDL.jl"))
#Pkg.add(PackageSpec(url="https://github.com/esa-esdl/ESDLPlots.jl"))
Pkg.add("ColorSchemes");
Pkg.add("Plots");
Pkg.add("Statistics")
Pkg.add("NetCDF")
Pkg.add("FFTW")
Pkg.add("DataFrames")
Pkg.add("Polynomials")

using RCall, Plots, ColorSchemes;
using Statistics;
using NetCDF;
using FFTW;
using DataFrames
using ESDL;
using Polynomials: fit

path = "Documents/Projects/DATA/ESA_ECV_ocean_salinity/"

Pkg.installed()["RCall"]

pwd()
cd("/Users/juanrocha/Documents/Projects/DATA/ESA_ECV_ocean_salinity/")

R"""
library(tidyverse)
fls <- fs::dir_ls(recurse = TRUE) |>
    str_subset( pattern = "nc$")
"""

@rget fls

## This formally reduces the computation to weekly observations:
fls = fls[1:7:length(fls)]

# this only gives me the files on the current directory, not recursively
# fls = readdir()
#
# fls = run(`ls -R `)
# read(fls, String)

## This is a lazy AbstractArray that does not load data on memory
c = NetCDF.open(fls[1])
ncinfo(fls[1])
## This is the data:
t = @timed dat = ncread(fls[1], "sss"); # 0.14s
t.time
## Extract lon and lat
lon = ncread(fls[1], "lon");
lat = ncread(fls[1], "lat");
size(dat)

m = map(!isnan, dat)

# it will be useful to set missing values for later
function set_missing(x)
    ## this one is for boolean, but I need one for nans.
    if x == 0
        x = missing
    else
        x
    end
    return(x)
end

function set_missing_nan(x)
    ## this one is for boolean, but I need one for nans.
    if isnan(x)
        x = missing
    else
        x
    end
    return(x)
end

m1 = map(set_missing, m)
m[250,250] == 1

## The bottle neck is reading over and over the files.  For each pixel time series, reading all 522 files takes over 11s. An alternative approach is to create the cube on memory:

# initialize one pixel
# pxl_ts = [0.0]
# # read the pixel in all the NetCDF files
# t = @timed for i in 1:length(fls)
#   pxl = ncread(fls[i], "sss", start = [250,250,1], count = [1,1,1])
#   pxl_ts = push!(pxl_ts, pxl[1,1,1])
# end
# @time deleteat!(pxl_ts, 1)
# @time fft(pxl_ts)

# initiallize the array
A = Array{Union{Float32,Missing}}(missing, 1388, 584, 522)

# 14secs, A is 1.9Gb! varinfo()
t = @timed for i in 1:length(fls)
    A[:,:,i] = ncread(fls[i], "sss")
end

A = map(set_missing_nan, A)

## testing: the problem with mapslices is that the result is a matrix
# mapslices(skipmissing ∘ )
T = A[20:30, 20:30, :]

@time l = size(A)[3]
@time outar = Matrix(undef, l, 4)

## plans for the Fourier transform
testar = zeros(Complex{Base.nonmissingtype(eltype(A[250,250,:]))}, l)
fftplan = plan_fft!(testar)
ifftplan = inv(fftplan)

#### fourier transform: inspiration from ESDL::filterTSFFT
## https://github.com/esa-esdl/ESDL.jl/blob/master/src/TSDecomposition.jl
function linreg(x,y)
  b = cov(x,y)/var(x)
  a = mean(y) - b*mean(x)
  a,b
end

function detrendTS!(outar::AbstractMatrix,xin::AbstractVector)
    x=collect(eachindex(xin))
    a,b=linreg(x,xin)
    for i in eachindex(xin)
        outar[i,1]=a+b*x[i]
    end
end

function tscale2ind(b::Float64,l::Int)
    i=round(Int,l/b+1)
    return i
end
mirror(i,l)=l-i+2

function filterFFT(outar, y, annfreq::Int=52, nharm::Int=3)
    detrendTS!(outar, y)
    fy = Complex{Base.nonmissingtype(eltype(y))}[y[i] - outar[i,1] for i=1:l]
    fftplan * fy # changes fy inplace
    fyout = similar(fy) # creates a new element array with same size
    czero = zero(eltype(fy))

    # remove annual cycle
    fill!(fyout, zero(eltype(fyout)))

    for jharm = 1:nharm
           iup=tscale2ind(annfreq*1.1/jharm,l)

           idown=tscale2ind(annfreq*0.9/jharm,l)

           for i=iup:idown
               fyout[i]  = fy[i]
               i2        = mirror(i,l)
               fyout[i2] = fy[i2]
               fy[i]     = czero
               fy[i2]    = czero
           end
       end

    ifftplan * fyout
    for i=1:l
        outar[i,3]=real(fyout[i])
    end

    iup   = tscale2ind(annfreq*1.1,l)
    idown = tscale2ind(annfreq*0.9,l)
    #Now split the remaining parts
    fill!(fyout,czero)
    for i=2:iup-1
        i2        = mirror(i,l)
        fyout[i]  = fy[i]
        fyout[i2] = fy[i2]
        fy[i]     = czero
        fy[i2]    = czero
    end
    ifftplan * fyout
    ifftplan * fy

    for i=1:l
        outar[i,2]=real(fyout[i])
        outar[i,4]=real(fy[i])
    end

end

@time filterFFT(outar, A[250,250,:])
#
# t = @timed df_all = DataFrame(lon = lon[250], lat = lat[250],
#    sss_slope = linreg(1:l, pxl_ts)[2],
#    sss_mean = mean(pxl_ts),
#    sss_var = var(pxl_ts),
#    sss_std_longterm = std(outar[:,2], mean=0),
#    sss_std_annual_cycle = std(outar[:,3], mean = 0),
#    sss_std_fast_oscillations = std(outar[:,4], mean = 0)
#    )
idx = 1 # counter
# precompute the computations worth doing, with few missing values
z = mapslices(sum, ismissing.(A), dims = [3])
# old conditional: ((!ismissing(T[i,j,:])) & (T[i,j,:] .|> ismissing |> sum < 100))
sum(z .< 20 )
response = Matrix(undef,  sum(z .< 20 ), 8)

t = @timed for i in 1:size(A)[1], j in 1:size(A)[2] # i = lons, j = lats
  if z[i,j] < 20
        outar = Matrix(undef, l, 4)
        x = gapfillpoly2(A[i,j,:])

        filterFFT(outar, x)

        response[idx,1] = lon[i] # longitud
        response[idx,2] = lat[j] # latitude
        response[idx,3] = linreg(1:l, x)[2] # sss_slope
        response[idx,4] = (mean ∘ skipmissing)(x) # sss_mean
        response[idx,5] = (var ∘ skipmissing)(x) # sss_var
        response[idx,6] = (std ∘ skipmissing)(outar[:,2]) # sss_std_longterm
        response[idx,7] = (std ∘ skipmissing)(outar[:,3]) # sss_std_annual_cycle
        response[idx,8] = (std ∘ skipmissing)(outar[:,4]) # sss_std_fast_oscillations
        idx += 1
  end
end
## It took ~24h to compllete almost 12k out of 505k. Cannot afford that
## with new approach

## save the results

R"""
df_all <- tibble(
    lon = unlist( $response[,1]),
    lat = unlist( $response[,2]),
    sss_slope = unlist( $response[,3]),
    sss_mean = unlist( $response[,4]),
    sss_var = unlist( $response[,5]),
    sss_std_longterm = unlist( $response[,6]),
    sss_std_annual_cycle = unlist( $response[,7]),
    sss_std_fast_oscillations = unlist( $response[,8])
)
"""

R"""
setwd("/Users/juanrocha/Documents/Projects/ESDL_earlyadopter/ESDL")
save(
  df_all,
  file = "Results/sampled_FFT_variables/sea_surface_salinity_FFT_.RData")
"""

x = Matrix(undef, size(A)[1], size(A)[2])
for i in 1:size(A)[1], j in 1:size(A)[2]
    if ismissing(A[i,j,:])
        x[i,j] = missing
    else
        x[i,j] = A[i,j,:] .|> ismissing |> sum
    end
end

## Additional unused code:

#### Some R manipulations
#
# R"tibble(ts = $pxl_ts, time = seq_along($pxl_ts)) %>% ggplot(aes(time, ts)) + geom_line() "
#
# R"""
# x <- as_tibble($dat)
# colnames(x) <- $lat
# x <- x %>% add_column(lon = $lon, .before = 1)
# x <- x |> pivot_longer(cols = 2:last_col(), names_to = "lat", values_to = "sss") |>
# mutate(lat = as.numeric(lat))
# """
#
# R"x"
#
# R"x |> ggplot(aes(lon,lat)) + geom_tile(aes(fill = sss))"
#
# # This reduces the dataset to 502k pixels
# R"x <- x |> filter(!is.nan(sss))"
#
# #### index of the lon and lat
# R"which(unique(x$lat) %in% $lat)"

function gapfillpoly2(x;max_gap=30,nbefore_after=10, polyorder = 2)
    x = replace(i->(!ismissing(i) && isfinite(i)) ? i : missing,x)
    a = copy(x)
    workx = Float64[]
    worky = Float64[]
    # find first nonmissing value
    idxstart = findfirst(!ismissing,a)
    idxstart === nothing && return a
    while true
        #find next missing value
        gapstart = findnext(ismissing,a,idxstart)
        gapstart === nothing && break
        gapstop = findnext(!ismissing,a,gapstart)
        gapstop === nothing && break
        if gapstop-gapstart < max_gap
            idxfirst = max(1,gapstart - nbefore_after)
            idxlast  = min(length(a), gapstop + nbefore_after - 1)
            idxr = idxfirst:idxlast
            for (ii,idx) in enumerate(idxr)
                if !ismissing(x[idx])
                    push!(workx,idxr[ii])
                    push!(worky,x[idx])
                end
            end
            if length(workx)>polyorder
                p = fit(workx,worky,polyorder)
                for idx in gapstart:(gapstop-1)
                    a[idx] = p(idx)
                end
            end
        end
        idxstart = gapstop
        empty!(workx)
        empty!(worky)
    end
    a
    # added this line to correct NAs on the ends of the vector
    y = replace(i->(!ismissing(i) && isfinite(i)) ? i : (mean ∘ skipmissing)(a),a)
    return y
end

@time x = gapfillpoly!(T[5,5,:])

mapslices(sum, ismissing.(T) ; dims = [3])

max_stars_repo_name: juanrocha/ESDL

max_issues_repo_name: juanrocha/ESDL

max_forks_repo_name: juanrocha/ESDL


text:
struct Node
    loc::Integer
    per::Integer
end

struct Mode
    Q::Integer
    con::Real
    dis::Real
    var::Real
end

struct Arc
    ori::Node #origin
    dst::Node #destination
    md::Mode #mode
end

#Accessors for AbstractArc and its sub
ori(a::Arc) = a.ori
dst(a::Arc) = a.dst
md(a::Arc) = a.md

struct Commodity
    src::Node #source
    snk::Node #sink
    ld::Integer #load
end

#Accessors for Commodities
src(k::Commodity) = k.src
snk(k::Commodity) = k.snk
ld(k::Commodity) = k.d

"""
    δ(n, k)
maps if n is 0 (intermediate node), 1 (source node), or -1 (sink node)
"""
function δ(n::Node, k::Commodity)
    a = (src(k) == n)
    b = -(snk(k) == n)

    return a + b
end

"""
    arc_fcnf(N,A,K)
builds the arc formulation of the fixed charge network flow (FCNF). The inputs are dictionaries and contains its iterators and all the information.
"""
function arc_fcnf(N::Vector{Node}, A::Vector{Arc}, K::Vector{Commodity})
    m = Model()

    @variable(m, 0 <= x[a = eachindex(A), k = eachindex(K)] <= 1)
    @variable(m, y[a = eachindex(A)], binary = true)

    @objective(m, Min,
        sum(f(A[a]) * y[a] for a in eachindex(A)) +
        sum(g(A[a]) * d(K[k]) * x[a, k] for a in eachindex(A), k in eachindex(K))
    )

    @constraint(m, flow[n = eachindex(N), k = eachindex(K)],
        sum(x[a, k] for a in eachindex(A) if ori(A[a]) == n) -
        sum(x[a, k] for a in eachindex(A) if dst(A[a]) == n) == δ(n, K[k])
    )

    @constraint(m, caps[a = eachindex(A)],
        sum(d(K[k]) * x[a, k] for k in eachindex(K)) <= u(A[a]) * y[a]
    )

    return m
end


text:
@json struct Branch
    name::String
    merged::Bool
    protected::Bool
    default::Bool
    developers_can_push::Bool
    developers_can_merge::Bool
    can_push::Bool
    commit::Commit
end

endpoint(::GitLabAPI, ::typeof(get_branch), owner::AStr, repo::AStr, branch::AStr) =
    Endpoint(:GET, "/projects/$(encode(owner, repo))/repository/branches/$branch")
into(::GitLabAPI, ::typeof(get_branch)) = Branch


text:
### Exports ###

export RewardModulator

### Types ###

mutable struct RewardModulator
  avgInput::Mean
end
RewardModulator() = RewardModulator(Mean(weight=ExponentialWeight()))

### Methods ###

early_modulate!(synapses, mod::Nothing, myconn, conns) = ()
late_modulate!(synapses, mod::Nothing, myconn, conns) = ()

function early_modulate!(synapses, mod::RewardModulator, myconn, conns)
  inputs = myconn[2][:]
  avgInput = mean(inputs)
  rewardModifier = avgInput - value(mod.avgInput)
  fit!(mod.avgInput, avgInput)

  # Update learnRate for each SynapticInput connection
  for conn in filter(conn->conn[1] isa SynapticInput, conns)
    conn[3][:learnRate] *= rewardModifier
  end
end
late_modulate!(synapses, mod::RewardModulator, myconn, conns) = ()

#=
@with_kw mutable struct FunctionalModulator{OF<:Function,IF}
  outer::OF
  inner::IF = nothing
end
modulate!(state, node, mod::FunctionalModulator) = mod.func(modulate!(state, node, mod.inner), node)
=#

#@with_kw mutable struct RewardModulator
#  avgReward::Mean = RewardModulator(Mean(weight=ExponentialWeight()))
#end
#=
function modulate!(state, node::GenericNeurons, mod::RewardModulator)
  rF = node.state.F
  rw = mean(rF)
  #fit!(learnRate, rw-value(mod.avgReward))
  #fit!(mod.avgReward, rw)
  # TODO
end
=#


text:
struct SamplePath{T,N,A,B} <: AbstractDiffEqArray{T, N, A}
    u::A # A <: AbstractVector{<: AbstractArray{T, N - 1}}
    t::B
end

Ensemble{T,N,A,B} = Vector{SamplePath{T,N,A,B}}

SamplePath(xs::AbstractVector{T}, ts, dims::NTuple{N}) where {T, N} = SamplePath{eltype(T), N, typeof(xs), typeof(ts)}(xs, ts)
# Assume that the first element is representative all all other elements
SamplePath(xs::AbstractVector, ts::AbstractVector) = SamplePath(xs, ts, (size(xs[1])..., length(xs)))

##### pretty printing
function Base.show(io::IO, xw::SamplePath)
    print(io,"t: "); show(io, xw.t)
    println(io);
    print(io,"x: "); show(io, xw.u)
    nothing
end

function Base.show(io::IO, m::MIME"text/plain", xw::SamplePath)
    print(io,"t: "); show(io,m,xw.t)
    println(io)
    print(io,"x: "); show(io,m,xw.u)
    nothing
end

##### update! logic
function update_samplepath!(f, xw::SamplePath, simulator, t, state, model, save_points)
    j = searchsortedlast(save_points, t)

    j == 0 && return xw

    i = j
    # backtrack in case the last event jumped over save points
    while (i ≥ 1) && !(save_points[i] in xw.t)
        i -= 1
    end

    # fill in the data for each save point
    for k in i+1:j
        push!(xw.u, f(simulator, state, model))
        push!(xw.t, save_points[k])
    end

    return xw
end

function update_samplepath!(f, xw::SamplePath, simulator, t, state, model, save_points::Nothing)
    push!(xw.u, f(simulator, state, model))
    push!(xw.t, t)

    return xw
end

# this just appends (t, state) to a SamplePath
function copy_to_samplepath!(xw::SamplePath, t, state, save_points)
    j = searchsortedlast(save_points, t)

    j == 0 && return xw

    i = j
    # backtrack in case the last event jumped over save points
    while (i ≥ 1) && !(save_points[i] in xw.t)
        i -= 1
    end

    # fill in the data for each save point
    for k in i+1:j
        push!(xw.u, state)
        push!(xw.t, save_points[k])
    end

    return xw
end

##### extracting save data
save_state(simulator, state::Vector{T}, model) where T <: Int = copy(state)
save_state(simulator, state::Lattice, model) = Configuration(state)
save_rates(simulator, state, model) = copy(jump_rates(simulator))

##### initializing output
function build_output(f, simulator, state, model, ::Nothing)

    xw = SamplePath([f(simulator, state, model)], [0.0])
    sizehint!(xw, 1_000)

    return xw
end

function build_output(f, simulator, state, model, ntrials::Integer)

    ensemble = [SamplePath([f(simulator, state, model)], [0.0]) for k in 1:ntrials]
    foreach(xw -> sizehint!(xw, 1000), ensemble)

    return ensemble
end

##### interpolation
function get_regular_path(xw::SamplePath, save_points)
    yw = SamplePath([xw.u[1]], [save_points[1]])
    t = xw.t
    u = xw.u

    # copy data from the sample path
    for j in 2:length(xw.u)
        copy_to_samplepath!(yw, t[j], u[j], save_points)
    end

    return yw
end

function get_regular_ensemble(x::Ensemble, save_points)
    return [get_regular_path(x[i], save_points) for i in eachindex(x)]
end

##### plotting recipes #####
__make_indices(vars::Integer) = vars
__make_indices(vars::NTuple{N,Integer}) where N = (N > 1 ? collect(vars) : vars[1])
__make_indices(vars::AbstractVector{T}) where T <: Integer = vars

@recipe function f(xw::SamplePath{T};
    summary = :trajectory,
    vars = nothing) where T <: Number

    if summary == :trajectory
        seriestype --> :steppre
        n = length(xw[1])
        if vars == nothing
            xw.t, xw'
        else
            idxs = __make_indices(vars)

            if idxs isa Integer
                xw.t, xw[idxs,:]
            else
                xw.t, xw[idxs,:]'
            end
        end
    elseif summary == :phase
        idxs = (vars == nothing) ? (1,2) : (vars[1], vars[2])
        xs = [xw[idxs[1],t] for t in eachindex(xw.t)]
        ys = [xw[idxs[2],t] for t in eachindex(xw.t)]

        legend --> nothing

        @series begin
            seriestype --> :scatter
            marker_z --> xw.t

            xs, ys
        end

        @series begin
            seriestype --> :path
            xs, ys
        end
    end
end

@recipe function f(ens::Ensemble;
    summary = :mean,
    save_points = nothing,
    error_style = :bars,
    timepoint = 1,
    vars = nothing)

    if save_points == nothing
        m = length(ens[1].t)
        if any(xw -> length(xw) != m, ens)
            t_values = 0:ens[1].t[end]
        else
            t_values = ens[1].t
        end
    else

    end

    reg = get_regular_ensemble(ens, t_values)

    if summary == :mean
        # extract the series data
        ts = reg[1].t

        tmp_xs = transpose(mean(reg))
        tmp_xstd = transpose(std(reg))

        if vars == nothing
            xs = tmp_xs
            xstd = tmp_xstd
        else
            idxs = __make_indices(vars)

            xs = tmp_xs[:, idxs]
            xstd = tmp_xstd[:, idxs]
        end

        # make the default a scatter plot
        # and add bars to represent standard deviation
        seriestype --> :scatter

        if error_style == :bars
            yerrorbar  --> xstd
        elseif error_style == :ribbon
            ribbon --> xstd
        else
            error("Unknown option for error_style: $(error_style)")
        end

        ts, xs
    elseif summary == :histogram
        t = max(1, searchsortedlast(reg[1].t, timepoint))
        default_idxs = 1:length(reg[1].u[1])
        idxs = (vars == nothing) ? default_idxs : vars

        seriestype --> :histogram

        xs = [reg[n][k,t] for n in eachindex(reg), k in idxs]
    elseif summary == :phase
        tmp = mean(reg)
        idxs = (vars == nothing) ? (1,2) : (vars[1], vars[2])
        xs = [tmp[idxs[1],t] for t in eachindex(reg[1].t)]
        ys = [tmp[idxs[2],t] for t in eachindex(reg[1].t)]

        legend --> nothing

        @series begin
            seriestype --> :scatter
            marker_z --> reg[1].t

            xs, ys
        end

        @series begin
            seriestype --> :path
            xs, ys
        end
    else
        error("Unknown option for summary_type: $(summary_type)")
    end
end


text:
############ Trades  ####################
"""
    stocks_trades(opts::PolyOpts, ticker::AbstractString, date::AbstractString; limit=10, reverse=true, kwargs...)

Get stock trades for a given ticker symbol on a specified date.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * ticker::AbstractString: The ticker symbol we want trades for.
 * date::AbstractString: The date/day of the trades to retrieve in the format YYYY-MM-DD.
 * limit::Int: Limit the size of the response, max 50000 and default 10.
 * reverse::Boolean: Reverse the order of the results.
 * kwargs::Any: A list of additional arguments to pass to the Polygon IO API.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_trades(opts, "AAPL", "2017-01-01")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array.
 * See https://polygon.io/docs/get_v2_ticks_stocks_trades__ticker___date__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_trades(opts::PolyOpts, ticker::AbstractString, date::AbstractString; limit=10, reverse=true, kwargs...)
    stocks_trades_url = "$stocks_trades_base_url/$ticker/$date"
    params = Dict(
        "apiKey" => opts.api_key,
        "limit" => limit,
        "reverse" => reverse
    )

    merge!(params, Dict(kwargs))

    return generate_output_from_url(YesSinkYesResults(), stocks_trades_url, params, opts.sink)
end

############ Quotes  ####################
"""
    stocks_quotes_nbbo(opts::PolyOpts, ticker::AbstractString, date::AbstractString; limit=10, reverse=true, kwargs...)

Get NBBO quotes for a given ticker symbol on a specified date.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * ticker::AbstractString: The ticker symbol we want quotes for.
 * date::AbstractString: The date/day of the quotes to retrieve in the format YYYY-MM-DD.
 * limit::Int: Limit the size of the response, max 50000 and default 10.
 * reverse::Bool: Reverse the order of the results.
 * kwargs::Any: A list of additional arguments to pass to the Polygon IO API.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_quotes_nbbo(opts, "AAPL", "2017-01-01")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_ticks_stocks_nbbo__ticker___date__anchor for documentation on response attributes and supported keyword arguments.

"""
function stocks_quotes_nbbo(opts::PolyOpts, ticker::AbstractString, date::AbstractString; limit=10, reverse=true, kwargs...)
    stocks_quotes_url = "$stocks_quotes_base_url/$ticker/$date"
    params = Dict(
        "apiKey" => opts.api_key,
        "limit" => limit,
        "reverse" => reverse
    )

    merge!(params, Dict(kwargs))

    return generate_output_from_url(YesSinkYesResults(), stocks_quotes_url, params, opts.sink)
end

# ############ Last Trade For A Symbol v2  ####################
"""
    stocks_last_trade_symbol(opts::PolyOpts, stocksTicker::AbstractString)

Get the most recent trade for a given stock.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString: The ticker symbol of the stock/equity.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_last_trade_symbol(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_last_trade__stocksTicker__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_last_trade_symbol(opts::PolyOpts, stocksTicker::AbstractString)
    stocks_last_trade_url = "$stocks_last_trade_base_url/$stocksTicker"
    params = Dict("apiKey" => opts.api_key)

    return generate_output_from_url(YesSinkYesResults(), stocks_last_trade_url, params, opts.sink)
end

# ############ Last Quote For A Symbol v2  ####################
"""
    stocks_last_quote_symbol(opts::PolyOpts, stocksTicker::AbstractString)

Get the most recent NBBO (Quote) tick for a given stock.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString: The ticker symbol of the stock/equity.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_last_quote_symbol(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_last_nbbo__stocksTicker__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_last_quote_symbol(opts::PolyOpts, stocksTicker::AbstractString)
    stocks_last_quote_url = "$stocks_last_quote_base_url/$stocksTicker"
    params = Dict("apiKey" => opts.api_key)

    return generate_output_from_url(YesSinkYesResults(), stocks_last_quote_url, params, opts.sink)
end

# ############ Daily Open/Close  ####################
"""
    stocks_daily_open_close(opts::PolyOpts, stocksTicker::AbstractString, date::AbstractString; adjusted=true)

Get the open, close and afterhours prices of a stock symbol on a certain date.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString: The ticker symbol of the stock/equity.
 * date::AbstractString: The date of the requested open/close in the format YYYY-MM-DD.
 * adjusted::Bool: Whether or not the results are adjusted for splits. By default, results are adjusted.
    Set this to false to get results that are NOT adjusted for splits.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_daily_open_close(opts, "AAPL", "2017-01-01")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v1_open-close__stocksTicker___date__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_daily_open_close(opts::PolyOpts, stocksTicker::AbstractString, date::AbstractString; adjusted=true)
    stocks_daily_open_close_url = "$stocks_daily_open_close_base_url/$stocksTicker/$date"
    params = Dict(
        "apiKey" => opts.api_key,
        "adjusted" => adjusted
    )

    return generate_output_from_url(YesSinkNoResults(), stocks_daily_open_close_url, params, opts.sink)
end

# ############ Grouped Daily (Bars)  ####################
"""
    stocks_grouped_daily_bars(opts::PolyOpts, date::AbstractString; adjusted=true)

Get the daily open, high, low, and close (OHLC) for the entire stocks/equities markets.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * date::AbstractString: The beginning date for the aggregate window.
 * adjusted::Bool: Whether or not the results are adjusted for splits. By default, results are adjusted.
    Set this to false to get results that are NOT adjusted for splits.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_grouped_daily_bars(opts, "2017-01-01")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_aggs_grouped_locale_us_market_stocks__date__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_grouped_daily_bars(opts::PolyOpts, date::AbstractString; adjusted=true)
    stocks_grouped_daily_bars_url = "$stocks_grouped_daily_bars_base_url/$date"
    params = Dict(
        "apiKey" => opts.api_key,
        "adjusted" => adjusted
    )

    return generate_output_from_url(YesSinkYesResults(), stocks_grouped_daily_bars_url, params, opts.sink)
end

# ############ Previous Close  ####################
"""
    stocks_previous_close(opts::PolyOpts, stocksTicker::AbstractString; adjusted=true)

Get the previous day's open, high, low, and close (OHLC) for the specified stock ticker.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString: The ticker symbol of the stock/equity.
 * adjusted::Bool: Whether or not the results are adjusted for splits.
    By default, results are adjusted. Set this to false to get results that are NOT adjusted for splits.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_previous_close(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_aggs_ticker__stocksTicker__prev_anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_previous_close(opts::PolyOpts, stocksTicker::AbstractString; adjusted=true)
    stocks_previous_close_url = "$stocks_previous_close_base_url/$stocksTicker/prev"
    params = Dict(
        "apiKey" => opts.api_key,
        "adjusted" => adjusted
    )

    return generate_output_from_url(YesSinkYesResults(), stocks_previous_close_url, params, opts.sink)
end

# ############ Aggregates (Bars)  ####################
"""
    stocks_aggregates_bars(opts::PolyOpts, stocksTicker::AbstractString;
                        multiplier=1, timespan="day", from="2020-10-14", to="2020-10-14",
                        adjusted=true, sort="asc", limit=120, kwargs...)

Get aggregate bars for a stock over a given date range in custom time window sizes.
For example, if timespan = ‘minute’ and multiplier = ‘5’ then 5-minute bars will be returned.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString: The ticker symbol of the stock/equity.
 * multiplier::Int: TThe size of the timespan multiplier.
 * timespan::AbstractString: The size of the time window.
 * from::AbstractString: The start of the aggregate time window.
 * to::AbstractString: The end of the aggregate time window.
 * adjusted::Bool: Whether or not the results are adjusted for splits.
    By default, results are adjusted. Set this to false to get results that are NOT adjusted for splits.
 * sort::AbstractString: Sort the results by timestamp. asc will return results in ascending order (oldest at the top),
    desc will return results in descending order (newest at the top).
 * limit::Int: Limits the number of base aggregates queried to create the aggregate results. Max 50000 and Default 120.
 * kwargs::AbstractString: Any additional arguments.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_aggregates_bars(opts, "AAPL", multiplier=5, timespan="minute", from="2017-01-01", to="2017-01-01")
```
# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_aggs_ticker__stocksTicker__range__multiplier___timespan___from___to__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_aggregates_bars(opts::PolyOpts, stocksTicker::AbstractString;
                            multiplier=1, timespan="day", from="2020-10-14", to="2020-10-14",
                            adjusted=true, sort="asc", limit=120, kwargs...)

    stocks_aggregates_bars_url = "$stocks_aggregates_bars_base_url/$stocksTicker/range/$multiplier/$timespan/$from/$to"

    params = Dict(
        "apiKey" => opts.api_key,
        "multiplier" => multiplier,
        "timespan" => timespan,
        "from" => from,
        "to" => to,
        "adjusted" => adjusted,
        "sort" => sort,
        "limit" => limit
    )

    merge!(params, Dict(kwargs))

    return generate_output_from_url(YesSinkYesResults(), stocks_aggregates_bars_url, params, opts.sink)

end

# ############ Snapshot - All Tickers  ####################
"""
    stocks_snapshot_all_tickers(opts::PolyOpts, tickers::AbstractString)

Get the current minute, day, and previous day’s aggregate, as well as the last trade and quote for all traded stock symbols.
Note: Snapshot data is cleared at 12am EST and gets populated as data is received from the exchanges. This can happen as early as 4am EST.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * tickers::AbstractString: A comma separated list of tickers to get snapshots for.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_snapshot(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_snapshot_locale_us_markets_stocks_tickers_anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_snapshot_all_tickers(opts::PolyOpts, tickers::AbstractString)
    stocks_snapshot_all_tickers_url = "$stocks_snapshot_all_tickers_base_url"

    params = Dict(
        "apiKey" => opts.api_key,
        "tickers" => tickers
    )

    return generate_output_from_url(NoSinkYesTickers(), stocks_snapshot_all_tickers_url, params, opts.sink)
end

# ############ Snapshot - Ticker  ####################
"""
    stocks_snapshot_ticker(opts::PolyOpts, stocksTicker::AbstractString)

Get the current minute, day, and previous day’s aggregate, as well as the last trade and quote for a single traded stock ticker.
Note: Snapshot data is cleared at 12am EST and gets populated as data is received from the exchanges. This can happen as early as 4am EST.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * stocksTicker::AbstractString: The ticker symbol of the stock/equity.

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_snapshot_ticker(opts, "AAPL")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_snapshot_locale_us_markets_stocks_tickers__stocksTicker__anchor for documentation on response attributes and supported keyword arguments.

"""
function stocks_snapshot_ticker(opts::PolyOpts, stocksTicker::AbstractString)
    stocks_snapshot_ticker_url = "$stocks_snapshot_ticker_base_url/$stocksTicker"
    params = Dict("apiKey" => opts.api_key)

    return generate_output_from_url(NoSinkYesTicker(), stocks_snapshot_ticker_url, params, opts.sink)
end

# ############ Snapshot - Gainers/Losers  ####################
"""
    stocks_snapshot_gainers_losers(opts::PolyOpts, direction="losers")

Get the current top 20 gainers or losers of the day in stocks/equities markets.
Top gainers are those tickers whose price has increased by the highest percentage since the previous day's close.
Top losers are those tickers whose price has decreased by the highest percentage since the previous day's close.
Note: Snapshot data is cleared at 12am EST and gets populated as data is received from the exchanges.

# Arguments
 * opts::PolyOpts: The PolyOpts object used to configure the request.
 * direction::AbstractString: The direction of the snapshot results to return. The direction can be "gainers" or "losers".

# Example
```julia-repl
julia> opts = PolyOpts(API_KEY, nothing)
julia> stocks_snapshot_gainers_losers(opts, "losers")
```

# Returns
 * A JSON3.Array or specified tabular representation of the JSON3.Array
 * See https://polygon.io/docs/get_v2_snapshot_locale_us_markets_stocks__direction__anchor for documentation on response attributes and supported keyword arguments.
"""
function stocks_snapshot_gainers_losers(opts::PolyOpts, direction="losers")
    stocks_snapshot_gainers_losers_url = "$stocks_snapshot_gainers_losers_base_url/$direction"
    params = Dict(
        "apiKey" => opts.api_key,
        "direction" => direction
    )

    return generate_output_from_url(NoSinkYesTickers(), stocks_snapshot_gainers_losers_url, params, opts.sink)
end


text:
# Intermittency simple Manneville map
using DrWatson
@quickactivate "NonlinearDynamicsTextbook"
include(srcdir("style.jl"))
using DynamicalSystems, PyPlot
lo = Systems.manneville_simple(0.454)

rs = [-0.1, 0.0001, 0.05]

fig = figure(figsize = (figx, 1.5figy))
ax1 = fig.add_subplot(3, 2, 1)
ax2 = fig.add_subplot(3, 2, 3; sharex = ax1)
ax3 = fig.add_subplot(3, 2, 5; sharex = ax1)
axs = (ax1, ax2, ax3)
axc = fig.add_subplot(1, 2, 2)
add_identifiers!(fig)

T = 180
xs = 0:0.0001:1

for i in 1:length(rs)
	set_parameter!(lo, 1, rs[i])
	x = trajectory(lo, T; Ttr = 100)
	chaotici = findall(a -> a>0.1, x)
	regulari = setdiff(1:length(x), chaotici)
	axs[i].plot(0:length(x)-1, x; lw = 0.5, color = "C0")
	axs[i].plot(regulari .- 1, x[regulari]; ls = "None", marker = "o", color = "C2")
	axs[i].plot(chaotici .- 1, x[chaotici]; ls = "None", marker = "o", color = "C0")
	axs[i].set_ylim(-0.1, 1.2)
	axs[i].set_yticks([0, 1])
	# axs[i].text(0.99, 0.90, "\$r=$(rs[i])\$"; bbox = bbox,
	# transform = axs[i].transAxes, va = :top, ha=:right, fontsize = 22)
end
# axs[1].set_xlim(0, T)
for ax in (ax1, ax2); ax.tick_params(labelbottom=false); end
axs[1].set_xlim(0, T)
axs[1].set_xticks(0:60:180)
axs[3].set_xlabel("\$n\$"; labelpad = -25)
axs[2].set_ylabel("\$x_n\$"; labelpad = -10)
axs[2].axvspan(21, 40; color = "C3", alpha = 0.5)

# laminar arrow
xc = (135 + 169)/2
xspan = (169 - 135)
nice_arrow!(ax2, xc, 1.0, xspan, 0; tex = "\$\\ell\$", yo = 0.0, xo = -xspan/2 - 5, color = "C2")

xc = (65 + 100)/2
xspan = abs(65 - 100)
nice_arrow!(ax2, xc, 1.0, xspan, 0; tex = "\$c\$", yo = 0.0, xo = xspan/2, color = "C0")

# Do the cobweb plot
function cobweb(t) # transform timeseries t into cobweb (point2D)
    cx = Float64[]; cy = Float64[]
    for i ∈ 1:length(t)-1
        push!(cx, t[i]); push!(cy, t[i])
		push!(cx, t[i]); push!(cy, t[i+1])
    end
    return cx, cy
end

axc.set_xlim(0.0, 1.0)
axc.set_ylim(0.0, 1.0)
axc.set_xticks([0,1])
axc.set_yticks([0,1])
axc.set_ylabel("\$x_{n+1}\$"; labelpad = -15)
axc.set_xlabel("\$x_{n}\$"; labelpad = -25)
axi = axc.inset_axes([0.05, 0.55, 0.3, 0.4])

for ax in (axc, axi)
	_r = rs[2]
	set_parameter!(lo, 1, _r)
	f = lo.f.(xs, Ref([_r]), 0)
	ax.plot(xs, f, color = "C1")
	ax.plot([0, 1], [0,1]; lw = 1.0, color = "k")
	x = trajectory(lo, 60; Ttr = 100)
	cx, cy = cobweb(x)
	ax.plot(cx, cy; lw = 1.0, alpha = 0.5, color = COLORS[1])
end

z = 0.01
w = 0.05
zbox = ((z, z), (z+w, z+w))
axis_zoomin!(axi, axc, zbox, zbox, "C3"; dir = :top)
axi.set_ylim(z, z+w)
axi.set_xlim(z, z+w)
axi.axis("off")

fig.tight_layout(pad = 0.3)
wsave(plotsdir("4", "intermittency_manneville"), fig)

max_stars_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_issues_repo_name: JuliaDynamics/NonlinearDynamicsTextbook

max_forks_repo_name: JuliaDynamics/NonlinearDynamicsTextbook


text:
using Base.Test
import ArchGDAL; const AG = ArchGDAL

@testset "Testing FeatureLayer Methods" begin
    AG.registerdrivers() do
        AG.read("data/point.geojson") do dataset
            AG.createcopy(dataset, "tmp/point.geojson") do tmpcopy
                @test AG.nlayer(tmpcopy) == 1
                # AG.deletelayer!(tmpcopy, 0)
                # @test AG.nlayer(tmpcopy) == 0
                tmplayer = AG.getlayer(tmpcopy,0)
                @test AG.isignored(AG.getgeomfielddefn(AG.getlayerdefn(tmplayer),0)) == false
                AG.setignoredfields!(tmplayer, ["OGR_GEOMETRY"])
                @test AG.isignored(AG.getgeomfielddefn(AG.getlayerdefn(tmplayer),0)) == true
                AG.synctodisk!(tmplayer)
            end
            rm("tmp/point.geojson")

            layer = AG.getlayer(dataset, 0)
            println(AG.getspatialref(layer))
            println("FID Column name: $(AG.getfidcolname(layer))")
            println("Geom Column name: $(AG.getgeomcolname(layer))")
            @test AG.nreference(layer) == 0
            AG.reference(layer)
            @test AG.nreference(layer) == 1
            AG.dereference(layer)
            @test AG.nreference(layer) == 0
            for feature in layer
                println(feature)
                @test AG.nreference(layer) == 0
                println("Features read: $(AG.getfeaturesread(layer))")
            end
            AG.getfeature(layer, 2) do feature
                print(feature)
            end
            AG.setspatialfilter!(layer,100,-1,100.1,1)
            println("Filtering all points in $(AG.getspatialfilter(layer))")
            AG.clone(AG.getspatialfilter(layer)) do poly
                println(poly)
                for feature in layer
                    println(feature)
                end
                AG.clearspatialfilter!(layer)
                println("Clearing Spatial Filter")
                println("New $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                AG.setspatialfilter!(layer, poly)
                println("Re-filter all points in $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                AG.clearspatialfilter!(layer)
                println("Clearing Spatial Filter")
                println("New $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                println("Test with setting to index of geomfield")
                AG.setspatialfilter!(layer, 0, poly)
                println("  Re-filter all points in $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                AG.clearspatialfilter!(layer, 0)
                println("  Clearing Spatial Filter")
                println("  New $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                AG.setspatialfilter!(layer, 0,100,-1,100.1,1)
                println("Re-filter all points in $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                AG.clearspatialfilter!(layer)
                println("Clearing Spatial Filter")
                println("New $(AG.getspatialfilter(layer))")
                for feature in layer
                    println(feature)
                end
                query = "FID = 2"
                println("Setting attribute filter: $query")
                AG.setattributefilter!(layer, query)
                for feature in layer
                    println(feature)
                end
                query = "FID = 3"
                println("Setting attribute filter: $query")
                AG.setattributefilter!(layer, query)
                for feature in layer
                    println(feature)
                end
                AG.clearattributefilter!(layer)
                println("After clearing attribute filter:")
                for feature in layer
                    println(feature)
                end
                println("Fast forward to index 2:")
                AG.setnextbyindex!(layer, 2)
                for feature in layer
                    println(feature)
                end
                print("Test capability for random access writing: ")
                println(AG.testcapability(layer,"OLCRandomWrite"))
            end
            println("FID exact index: $(AG.findfieldindex(layer,"FID", true))")
            println("FID index: $(AG.findfieldindex(layer,"FID", false))")
            println("pointname exact index: $(AG.findfieldindex(layer,"pointname", true))")
            println("pointname index: $(AG.findfieldindex(layer,"pointname", false))")
            println("geom exact findfieldindex: $(AG.findfieldindex(layer,"geom", true))")
            println("geom findfieldindex: $(AG.findfieldindex(layer,"geom", true))")
            println("rubbish exact index: $(AG.findfieldindex(layer,"rubbish", true))")
            println("rubbish index: $(AG.findfieldindex(layer,"rubbish", false))")

            println(AG.getextent(layer, 0, true))
            println(AG.getextent(layer, true))
        end
    end
end

# Untested:
# setfeature!(layer::FeatureLayer, feature::Feature)
# deletefeature!(layer::FeatureLayer, i::Integer)
# creategeomfield!(layer::FeatureLayer, field::GeomFieldDefn, approx::Bool = false)
# deletefield!(layer::FeatureLayer, i::Integer)
# reorderfields!(layer::FeatureLayer, indices::Vector{Cint})
# reorderfield!(layer::FeatureLayer, oldpos::Integer, newpos::Integer)
# alterfielddefn!(layer::FeatureLayer, i::Integer, newfielddefn::FieldDefn, flags::UInt8)
# setignoredfields!(layer::FeatureLayer, fieldnames)
# intersection(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.C.GDALDummyProgress, progressdata = C_NULL)
# union(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressdata = C_NULL, progressfunc::Function = GDAL.C.GDALDummyProgress)
# symdifference(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.C.GDALDummyProgress, progressdata = C_NULL)
# identity(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.C.GDALDummyProgress, progressdata = C_NULL)
# update(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.C.GDALDummyProgress, progressdata = C_NULL)
# clip(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.C.GDALDummyProgress, progressdata = C_NULL)
# erase(input::FeatureLayer, method::FeatureLayer, result::FeatureLayer; options = StringList(C_NULL), progressfunc::Function = GDAL.C.GDALDummyProgress, progressdata = C_NULL)


text:
"""
    NoOp <: AugmentorBase.AffineOperation

Identity transformation that does not do anything with the given
image, but instead passes it along unchanged (without copying).

Usually used in combination with [`Either`](@ref) to denote a
"branch" that does not perform any computation.
"""
struct NoOp <: AffineOperation end

@inline supports_eager(::Type{NoOp}) = false
@inline supports_stepview(::Type{NoOp}) = true
@inline supports_view(::Type{NoOp}) = true

# TODO: implement method for n-dim arrays
applyeager(::NoOp, input::AbstractArray, param) = maybe_copy(input)
applylazy(::NoOp, input::AbstractArray, param) = input

function applyview(::NoOp, input::AbstractArray, param)
    idx = map(i->1:length(i), axes(input))
    indirect_view(input, idx)
end

function applystepview(::NoOp, input::AbstractArray, param)
    idx = map(i->1:1:length(i), axes(input))
    indirect_view(input, idx)
end

function showconstruction(io::IO, op::NoOp)
    print(io, typeof(op).name.name, "()")
end

function Base.show(io::IO, op::NoOp)
    if get(io, :compact, false)
        print(io, "No operation")
    else
        print(io, "AugmentorBase.")
        showconstruction(io, op)
    end
end


text:
module Tree

importall ImmutableArrays
import ..Mesh.Face
import Base:cmp,search

type Leaf
	faces::Vector{Face}
end

type Node
	front::Union(Node,Leaf)
	back::Union(Node,Leaf)
	normal::Vector3{Float32}
	distance::Float32
end

cmp(pos::Vector3, node::Node) = sign(dot(pos,node.normal)-node.distance)

function search(tree::Node, pos::Vector3)
	if cmp(pos, tree) < 0
		return search(tree.back, pos)
	else
		return search(tree.front, pos)
	end
end
search(leaf::Leaf, ::Vector3) = leaf

end


text:
module CUDA

##### CUDA #####
if is_windows()
    const libcuda = Libdl.find_library(["nvcuda"])
else
    const libcuda = Libdl.find_library(["libcuda"])
end
isempty(libcuda) && error("CUDA library cannot be found.")

function check_curesult(status)
    status == CUDA_SUCCESS && return nothing
    warn("CUDA error triggered from:")
    Base.show_backtrace(STDOUT, backtrace())
    ref = Ptr{UInt8}[0]
    cuGetErrorString(status, ref)
    throw(unsafe_string(ref[1]))
end

ref = Cint[0]
ccall((:cuDriverGetVersion,libcuda),UInt32,(Ptr{Cint},),ref)
const driver_version = Int(ref[1])
const major = div(driver_version, 1000)
const minor = div(driver_version - major*1000, 10)

include("../lib/$(major).$(minor)/libcuda.jl")
include("../lib/$(major).$(minor)/libcuda_types.jl")

include("context.jl")
include("device.jl")
include("module.jl")
include("function.jl")
include("headers.jl")

info("CUDA driver version: $driver_version")
initctx()
infodevices()

ctype(::Type{Int64}) = :int
ctype(::Type{Float32}) = :float
ctype(::Type{Float64}) = :double

include("pointer.jl")
include("abstractarray.jl")
include("array.jl")
include("arraymath.jl")
include("reducedim.jl")
include("subarray.jl")
include("Interop.jl")

##### NVRTC #####
include("NVRTC.jl")

##### CUBLAS #####
include("CUBLAS.jl")
using .CUBLAS

##### CUDNN #####
include("cudnn/CUDNN.jl")

end


text:
################################################################################
#
# Generate synthetic data for the BG Compass model
#
# Author: Mathias Loubout and Gabbrio Rizzuti
# Date: September 2020
################################################################################
### Module loading

using JUDI, JLD2

### Load true model
base_path = dirname(pathof(JUDI))*"/../examples/twri/"
@load string(base_path*"data/GaussLens.jld") n d o m

# Model

model_true = Model(n, d, o, m; nb=80)

### Acquisition geometrye
dt = 1f0
T = 2000f0 # total recording time [ms]

# Source wavelet
freq_peak = 0.006f0
wavelet = ricker_wavelet(T, dt, freq_peak)

# Sources
ix_src = 3:14:199
nsrc = length(ix_src)
iz_src = 3
x_src = convertToCell((ix_src.-1)*d[1])
y_src = convertToCell(range(0f0, stop = 0f0, length = nsrc))
z_src = convertToCell(range((iz_src-1)*d[2], stop = (iz_src-1)*d[2], length = nsrc))

# Receivers
ix_rcv = 3:199
iz_rcv = 199
nrcv = length(ix_rcv)
x_rcv = (ix_rcv.-1)*d[1]
y_rcv = 0f0
z_rcv = range((iz_rcv-1)*d[2], stop = (iz_rcv-1)*d[2], length = nrcv)

# Geometry structures
src_geom = Geometry(x_src, y_src, z_src; dt = dt, t = T)
rcv_geom = Geometry(x_rcv, y_rcv, z_rcv; dt = dt, t = T, nsrc = nsrc)

# Source function
fsrc = judiVector(src_geom, wavelet)

# Set up info structure for linear operators
ntComp = get_computational_nt(src_geom, rcv_geom, model_true)
info = Info(prod(n), nsrc, ntComp)

# Setup operators
F = judiModeling(info,model_true, src_geom,rcv_geom)
dat = F*fsrc

### Saving data
base_path = dirname(pathof(JUDI))*"/../examples/twri/"
@save string(base_path*"data/GaussLens_data_acou.jld") model_true fsrc dat

max_stars_repo_name: nogueirapeterson/JUDI

max_issues_repo_name: nogueirapeterson/JUDI

max_forks_repo_name: nogueirapeterson/JUDI


text:
module EOQ

using DecFP
export @d128_str, @d32_str, @d64_str, @d_str, Dec128, Dec32, Dec64, DecFP, exponent10, ldexp10, sigexp
"""
Simple EOQ model with required arguments:
    - demand (constant per unit time), order (set up) cost, holding (carrying) cost
    - holding cost can be calculated as (unit cost)*(holding rate)
Optional arguments:
    - time units ---> Tuple ---> default is ('m' for months, 12) => 1 year
    - lead time (default is 0)
"""
function eoq_execute(;demand::Signed, c::Union{Signed,AbstractFloat}, ct::Union{Signed,AbstractFloat},
    ce::Union{Signed,AbstractFloat}, time_units::Tuple{String,Signed}, lead_time::Signed)
    # Basic Output
    Q = sqrt(2*ct*demand/ce)
    T_raw = Q/demand
    T = Dec64(T_raw*time_units[2])
    N = 1/T_raw

    # Inventory Position
    # ROP_IP = demand*lead_time/time_units[2]

    # Costs (not accounting for pipeline inventory cost, i.e. for lead_time > 0)
    OrderingCost = ct*N
    HoldingCost = ce*(Q/2)
    TRC() = OrderingCost + HoldingCost
    TC() = TRC() + c*demand

    # Inventory Policy
    policy() = Dict("Q"=>Q, "T"=>T, "N"=>N, "TRC"=>TRC(), "TC"=>TC())

    function policy_for(s::String, value::Union{Signed, AbstractFloat})
        if s == "Q"
            Q_for = value
            T_raw_for = value/demand
            T_for = T_raw_for*time_units[2]
            N_for = 1/T_raw_for
        elseif s == "T"
            T_for = value
            N_for = time_units[2]/value
            Q_for = demand/N_for
        else
            return "Invalid input"
        end
        OrderingCost_for = ct*N_for
        HoldingCost_for = ce*(Q_for/2)
        TRC_for() = OrderingCost_for + HoldingCost_for
        TC_for() = TRC_for() + c*demand
        return Dict("Q*"=>Q_for, "T"=>T_for, "N"=>N_for, "TRC"=>TRC_for(), "TC" => TC_for(), "ΔTC"=> TC_for()/TC())
    end

    # [Required] Returning a function with "public" elements
    ()->(Q, T_raw, T, N,
    OrderingCost, HoldingCost, TRC, TC,
    policy, policy_for)
end

"""
Intermediary function to simplify input for eoq_execute()
"""
function eoq(;demand, c=0, ct, ce, time_units=("m",12), lead_time=0)
    d = Dict(:demand=>demand, :c=>c, :ct=>ct, :ce=>ce, :time_units=>time_units, :lead_time=>lead_time)
    return eoq_execute(;d...)
end

end # Module EOQ


text:
function evalcompound(u_rad, A::SHType{T}, B::FIDModelType{T,SST})::Complex{T} where {T <: Real, SST}

    #u_rad = 2*π*u

    out_sys = evalcLcompoundviapartitions(u_rad, A.αs, A.Ωs, B.ss_params,
    B.λ0, A.Δc_bar, A.part_inds_compound)

    out_singlets = evalsinglets(u_rad, B.d_singlets, A.αs_singlets, A.Ωs_singlets,
    B.β_singlets, B.λ0, B.κs_λ_singlets)

    return out_sys + out_singlets
end

function evalmixture(u_rad, As::Vector{SHType{T}}, Bs::Vector{FIDModelType{T,SST}};
    w::Vector{T} = ones(T, length(As)))::Complex{T} where {T <: Real, SST}

    #u_rad = 2*π*u

    out = zero(Complex{T})
    for n = 1:length(As)

        out += w[n]*evalcompound(u_rad, As[n], Bs[n])
    end

    return out
end

"""
Evaluates an element from the partition of a spin group.
r := 2*π*u-d, d is chem shift of this element in radians.
Does not evaluate the complex phase β.
"""
function evalcLpartitionelement(r,
    α::Vector{T}, Ω::Vector{T}, λ::T)::Complex{T} where T <: Real

    out = sum( α[l]/(λ+im*(r-Ω[l])) for l = 1:length(α) )

    return out
end

function evalcLcompoundviapartitions(u_rad,
    αs::Vector{Vector{T}}, Ωs::Vector{Vector{T}},
    x::SpinSysParamsType1{T}, λ0::T,
    c, part_inds_compound)::Complex{T} where T <: Real

    #u_rad = 2*π*u

    out = zero(Complex{T})
    for i = 1:length(αs)
        r = u_rad - x.d[i]

        for k = 1:length(part_inds_compound[i])
            inds = part_inds_compound[i][k]

            out += evalcLpartitionelement(r, αs[i][inds],
                Ωs[i][inds], x.κs_λ[i]*λ0)*exp(im*dot(x.κs_β[i], c[i][k]))
        end
    end

    return out
end

function evalcLcompoundviapartitions(u_rad,
    αs::Vector{Vector{T}}, Ωs::Vector{Vector{T}},
    x::SpinSysParamsType2{T}, λ0::T,
    c, part_inds_compound)::Complex{T} where T <: Real

    #u_rad = 2*π*u

    out = zero(Complex{T})
    for i = 1:length(αs)

        for k = 1:length(part_inds_compound[i])
            r = u_rad - x.d[i][k]
            inds = part_inds_compound[i][k]

            out += evalcLpartitionelement(r, αs[i][inds],
                Ωs[i][inds], x.κs_λ[i][k]*λ0)*exp(im*dot(x.κs_β[i], c[i][k]))
        end
    end

    return out
end