text:
"""
initforw allocates register storage (out0) and sets register fields
and flags if these are uninitialized or there is a change in input
size, type, or keyword arguments.  It does not zero or reset the
registers, there is a separate reset!  function for that.  FNN's do
not need reset, and RNNs only do at the beginning of a sequence.  Here
RNN refers to any network with a read-before-write register regardless
of whether it is applied to a sequence.
"""
function initforw(f::Net, inputs...; o...)
    nextforw = (map(typeof, inputs), map(size, inputs), o)
    if isvoid(f,:lastforw) || f.lastforw!=nextforw
        #setp(f,:forwoverwrite,false); setp(f,:backoverwrite,false)
        if isvoid(f,:lastforw) || (o != f.lastforw[3]) # when conditions change, the comp graph changes, effects argv, fanout, save
            initcond(f; o...)       # o,name,cond -> :forw
            initargv(f)             # :forw,name,args -> argv
            initsave(f)             # :forw,op,argv -> :save
            initfanout(f)           # :forw,name,argv -> :fanout
        end                         # changing input size/type or argv may change size and type of registers
        initsize(f, inputs...)      # :forw,inputs,op -> :size
        inittype(f, inputs...)      # :forw,inputs,op -> :eltype,:outtype
        initout0(f)                 # :forw,argv,op,:outtype,:eltype,:size,:forwoverwrite,:save,:fanout -> out0
        f.lastforw = nextforw
        gpusync()
    end
    length(inputs) == ninputs(f) || throw(ArgumentError())
end

initout0(f::Net;o...)=(for p in forwregs(f); initout0(f,p;o...); end)

function initout0(f::Net, p::Reg; seq=false)
    at, et, sz = getp(p,:outtype), getp(p,:eltype), getp(p,:size)
    if checkarray(p, :out0, at, et, sz)
        # all done
    elseif !isvoid(p, :out0) && (isa(p.op, Par) || isa(p.op, Arr))
        error("Size or type change not allowed in parameters and constants")
    else
        p.out0 = nothing
        if canoverwrite(p.op) && getp(p, :forwoverwrite, true)
            for i in reverse(p.argv) # consider overwriting the last input first
                q = reg(f,i)
                if canshare(f,p,q)
                    p.out0 = q.out0
                    @dbg info("Overwrite $(findfirst(f.reg,q))->$(findfirst(f.reg,p))=$(Int(pointer(p.out0))%1000)")
                    break
                end
            end
        end
        if p.out0 == nothing
            p.out0 = newarray(at, et, sz)
            gpusync()
            @dbg info("Alloc $(findfirst(f.reg,p))=$(Int(pointer(p.out0))%1000)")
        end
    end
end

function checkarray(r::Reg, n::Symbol, atype::DataType, etype::DataType, dims::Dims)
    isdefined(r,n) &&
    isa(r.(n), atype) &&
    eltype(r.(n)) == etype &&
    size(r.(n)) == dims
end

function canshare(f::Net, p::Reg, q::Reg)
    # checkarray(q,:out0,at,et,sz) && !ispersistent(q) && !getp(q, :save) && getp(q, :fanout)==1
    ispersistent(q) && return false
    getp(q,:save) && return false
    getp(q,:fanout) > 1 && return false
    at, et, sz = getp(p,:outtype), getp(p,:eltype), getp(p,:size)
    checkarray(q,:out0,at,et,sz) || return false
    # a=b*c; c=sigm(a) should not share, dot cannot overwrite!
    # a=b*c; d=e*f; e=a+d; e should not share a, a overwrites it!
    p.out0 = q.out0
    ok = checkshare(f)
    p.out0 = nothing
    return ok
end

isvoid(x,n)=(!isdefined(x,n) || isa(x.(n),Void))

function checkshare(f::Net)
    d = ObjectIdDict()
    for i=1:length(f)
        ri=f.reg[i]
        getp(ri,:forw) || continue
        !isvoid(ri,:out0) || continue
        d[ri.out0]=i
    end
    for i=1:length(f)
        ri=f.reg[i]
        getp(ri,:forw) || continue
        for j in ri.argv
            rj=f.reg[j]
            !isvoid(rj,:out0) || continue
            k = d[rj.out0]
            k == j || return false # warn("$j was overwritten by $k before being read by $i.")
            !isvoid(ri,:out0) || continue
            !canoverwrite(ri.op) && ri.out0===rj.out0 && return false # warn("$i and $j share array in nonoverwriting op.")
        end
        !isvoid(ri,:out0) && (d[ri.out0] = i)
    end
    return true
end

newarray(::Type{Array}, t::DataType, d::Dims)=Array(t,d)
newarray(::Type{SparseMatrixCSC}, t::DataType, d::NTuple{2,Int})=SparseMatrixCSC(d[1], d[2], ones(Cint, d[2]+1), Array(Cint, 0), Array(t, 0))
@gpu newarray(::Type{CudaArray}, t::DataType, d::Dims)=CudaArray(t,d)
@gpu newarray(::Type{CudaSparseMatrixCSC}, t::DataType, d::NTuple{2,Int})=CudaSparseMatrixCSC(t, fillsync!(CudaArray(Cint, d[2]+1), 1), CudaArray(Cint, 0), CudaArray(t, 0), d)
@gpu newarray(::Type{CudaSparseMatrixCSR}, t::DataType, d::NTuple{2,Int})=CudaSparseMatrixCSR(t, fillsync!(CudaArray(Cint, d[1]+1), 1), CudaArray(Cint, 0), CudaArray(t, 0), d)
@gpu newarray(::Type{CudaSparseMatrixCSCU}, t::DataType, d::NTuple{2,Int})=CudaSparseMatrixCSCU(t, d...)
@gpu newarray(::Type{CudaSparseMatrixCSRU}, t::DataType, d::NTuple{2,Int})=CudaSparseMatrixCSRU(t, d...)

# We do not need to copy persistent registers, they are guaranteed not to change during forwback.
ispersistent(p::Reg)=(isa(p.op,Par) || isa(p.op,Arr))
isreturn(p::Reg)=(p.name==:return)

# initcond(f::Net; o...) sets the :forw property of each f.reg[n] which
# is determined by evaluating f.reg[n].cond in the context of the
# kwargs passed to forw.

function initcond(f::Net; o...)
    for p in regs(f)
        setp(p,:forw,false)
    end
    d = Dict(o)
    for p in regs(f)
        c = condeval(p.cond,d)
        setp(p,:forw,c)
        c && isreturn(p) && break
    end
end

function condeval(s,d::Dict)
    (isa(s, Bool) ? s :
     isa(s, Symbol) ? get(d,s,false) :
     !isa(s, Expr) ? error("Expected boolean expression got $s") :
     s.head == :&& ? mapreduce(x->condeval(x,d), &, true, s.args) :
     s.head == :|| ? mapreduce(x->condeval(x,d), |, false, s.args) :
     s.head == :call && s.args[1] == :! ? !condeval(s.args[2],d) :
     error("Expected boolean expression got $s"))
end

function forwregs(f::Net)
    a = Reg[]
    for r in regs(f)
        getp(r,:forw) && push!(a,r)
    end
    return a
end

function initargv(f::Net)
    # There could be repeated and read-before-write variables
    s2i = Dict{Symbol,Int}()
    for n=1:length(f)
        p = reg(f,n)
        getp(p,:forw) || continue
        s2i[p.name] = n
    end
    for n=1:length(f)
        p = reg(f,n)
        getp(p,:forw) || continue
        p.argv = convert(Vector{Int}, map(s->s2i[s], p.args))
        s2i[p.name] = n
    end
end

function initsize(f::Net, inputs...)
    for p in forwregs(f)
        setp(p,:size,nothing)
    end
    lastinput = 0
    for p in forwregs(f)
        isa(p.op, Input) || continue
        setp(p,:size,size(inputs[lastinput += 1]))
    end
    nodims(p)=(s=getp(p,:size); s==nothing || prod(s)==0)
    notfound = count(nodims, forwregs(f))
    while notfound > 0
        for p in forwregs(f)
            isa(p.op, Input) && continue
            dy = getp(p,:size)
            xx = inputregs(f,p)
            dx = map(x->getp(x,:size), xx)
            d = infersize(p.op, dx..., dy)
            d == nothing && continue
            setp(p,:size,d[end])
            for i=1:length(xx)
                setp(xx[i],:size,d[i])
            end
        end
        nf = count(nodims, forwregs(f))
        nf == notfound && error("Cannot infer sizes")
        notfound = nf
    end
end

# inittype: Determines element and array type of p.out0 for each register p.
function inittype(f::Net, inputs...)
    isempty(inputs) && error("Don't know how to infer eltype with inputless networks yet.")
    ftype = eltype(inputs[1])
    for i in inputs
        eltype(i) == ftype || error("Conflicting element type in input $i.")
    end
    lastinput = 0
    for p in forwregs(f)
        setp(p,:eltype,ftype)
        if isa(p.op,Input) && issparse(inputs[lastinput+=1])
            setp(p,:outtype, (gpu() ? CudaSparseMatrixCSC : SparseMatrixCSC))
        else
            setp(p,:outtype, (gpu() ? CudaArray : Array))
        end
    end
end

# initsave(f::Net) sets the :save property of each f.prog[n] which
# should be true if the result of f.prog[n] will be needed for back
# calculation.  We find this out using back_reads_x and back_reads_y
# on each op as well as checking if the output will be returned and
# thus will be needed for back loss and grad calculation.  Any
# register with a save flag should not be overwritten in forw.  We'll
# be conservative and also not overwrite in sforw, even though it
# saves registers on the stack.

function initsave(f::Net)
    for p in regs(f); setp(p,:save,false); end
    for p in forwregs(f)
        if isreturn(p) || back_reads_y(p.op)
            setp(p,:save,true)
        end
        if back_reads_x(p.op)
            for i in p.argv
                setp(reg(f,i),:save,true)
            end
        end
    end
end

# getp(p,:fanout)>1: these registers are read at multiple locations in
# the program.  Depends on :cond.  :incr flag set for both sback and
# back.  The return registers are read by the user and compared to
# gold, which should count as 1.  All Par registers become
# multi-output during sequences and need :incr in sback.  Multi-output
# registers can forwoverwrite but should not be overwritten.  More
# specifically they could be overwritten by the last reader, but that
# is probably not worth the complexity.

function initfanout(f::Net)
    regs = forwregs(f)
    for p in regs; setp(p,:fanout,0); end
    for p in regs
        isreturn(p) && incp(p,:fanout)
        for i in inputregs(f,p)
            @assert getp(i,:forw)
            incp(i,:fanout)
        end
    end
end

### DEAD CODE:

# """
# outputs(inputs) returns an array of output indices for each register.
# Note that the final register is the network output, so it could be
# read externally even if outputs[N] is empty.
# """
# function outputs(inputs)
#     N = length(inputs)
#     outputs = [ Int[] for n=1:N ]
#     for n=1:N
#         for i in inputs[n]
#             push!(outputs[i], n)
#         end
#     end
#     push!(outputs[N], 0)        # for network output
#     return outputs
# end

    # if ygold != nothing
    #     @assert issimilar2(ygold, r.out0[N])
    #     if (r.dif0[ N]!=nothing)
    #         @assert issimilar3(ygold, r.dif0[N])
    #     else
    #         r.dif0[N] = newarray(gpu(), stype(ygold), eltype(ygold), size(ygold))
    #     end
    # end

    # if !seq
    #     # @assert all((r.out .== nothing) | (r.out .== r.out0)) # t:110/244
    #     @assert !keepstate "meaningless keepstate in non-sequence run"
    #     fillsync!(r.out, nothing)
    # elseif keepstate
    #     copysync!(r.out, r.out0)
    # else
    #     # @assert all((r.out .== nothing) | (r.out .== r.out0))
    #     fillsync!(r.out, nothing)
    # end

# keepstate=false, seq=false

# Switching to user controlled register sharing instead of automatic optimization...
# Safer.  Also makes reg(f,:name) work as expected.

# function findout0(f::Net, n, sizes, nsparse)
#     qn = f.prog[n]
#     ispersistent(qn) && return
#     free = nothing               # search most recent written first to avoid copying in overwriting ops
#     for i = n-1:-1:1             # considering overwriting i with n
#         qi = f.prog[i]
#         ispersistent(qi) && continue
#         !overwrites(qn.op) && in(qi.output, qn.inputs) && continue
#         ri = f.reg[qi.output]
#         isdefined(ri,:out0) || continue
#         size(ri.out0) == sizes[n] || continue
#         stype(ri.out0) == nsparse || continue
#         willberead = false                              # is anybody going to read i before it is written again?
#         k = n
#         while true
#             k = mod1(k+1, length(f.prog))
#             qk = f.prog[k]
#             for j in qk.inputs
#                 rj = f.reg[j]
#                 isdefined(rj,:out0) && rj.out0 === ri.out0 && (willberead = true; break)
#             end
#             willberead && break
#             rk = f.reg[qk.output]
#             isdefined(rk,:out0) && rk.out0 === ri.out0 && break
#         end
#         !willberead && (free = ri.out0; break)
#     end
#     return free
# end

# function initout0(f::Net, inputs...; o...)
#     xtype = infertype(f, inputs...)
#     sizes = infersize(f, inputs...) # todo: fix infersize
#     lastinput = 0
#     for p in f.prog
#         r = f.reg[p.output]
#         # We may keep arrays pre-allocated due to parameter sharing or batch resizing
#         # Only safe for Par/Arr, others may have sharing that no longer is valid
#         if isdefined(r,:out0) && (isa(p.op,Par) || isa(p.op,Arr))
#             size(r.out0)==sizes[n] || error("Size incompatible with network parameters")
#         else
#             r.out0 = nothing
#         end
#     end
#     # todo: infersize should make sure same register has same size
#     # should we leave array sharing to the user? x=relu(x)
#     # what if they type x=dot(x,y)?
#     # how about array sharing going back?  different rules?
#     for n=1:nops(f)
#         p = f.prog[n]
#         # todo: this is all wrong, the register has out not  p.
#         p.out = nothing
#         if p.out0 != nothing && (isa(p.op,Par) || isa(p.op,Arr))
#             continue
#         end
#         ### ARRAY SHARING OPTIMIZATION:
#         st = (isa(p.op, Input) ? stype(inputs[lastinput += 1]) : nothing)
#         p.out0 = findout0(f, n, sizes, st)
#         if p.out0 == nothing
#             p.out0 = newarray(gpu(), st, xtype, sizes[n])
#         end
#     end
#     # todo: figure out tmp
# end

# it = nothing
# for i in inputs
#     t = eltype(i)
#     if it == nothing
#         it = t
#     elseif t != it
#         error("Conflicting input eltypes")
#     end
# end
# it == nothing && error("Cannot infer eltype")
# return it
# # todo: deal with inputless networks:

# ispersistent(p::Ins)=(isa(p.op,Par) || isa(p.op,Arr) || getp(p,:push))

# initsize infers the size of each register

        # This will happen with conditionals:
        # save && f.lastforw[1] && Base.warn_once("Net reset during RNN training")
        # initsave(f, save)       # set :push

ext: jl


text:
function wrap_cell(nc, ic, jc) #ic e jc são os indices de uma celula
    a = ic%nc[1]
    if a == 0
        iw = nc[1]
    elseif a > 0
        iw = a
    else
        iw = a+nc[1]
    end
    b = jc%nc[2]
    if b == 0
        jw = nc[2]
    elseif b > 0
        jw = b
    else
        jw = b+nc[2]
    end
    return iw, jw
end


text:
include("HelperFunctions.jl")

using LinearAlgebra
using SparseArrays

"""
    calculates initial conditions for parameter estimation.
    parameters:
        Y : DataFrame | Array{Number, 2}
            matrix of variables, size (n, z). Not including date column.
        dates : Array{Dates.Date, 1}
            array of corresponding dates for the rows of Y
        p : number of lags in transition equation
        blocks : DataFrame, size (n_variables, n_blocks)
            matrix of 1s or 0s for block loadings, i.e. included in block.
            ex for Y 2 x 3, 2 blocks:
                [1, 0]
                [1, 1]
                [0, 1]
        R_mat : Array{Int64, 2}
            quarterly-monthly aggregation scheme matrix
    returns: Dictionary
        :A => transition matrix
        :C => measurement matrix
        :Q => Covariance for transition equation residuals
        :R => Covariance for measurement equation residuals
        :Z0 => Initial value of state
        :V0 => Initial value of covariance matrix
"""
function initialize_conditions(Y; dates, p, blocks, R_mat)
    y_tmp = copy(Y) |> standardize
    monthly_quarterly_array = gen_monthly_quarterly(dates, y_tmp)
    # put all quarterly variables to the end
    y_tmp = [y_tmp[!, .!monthly_quarterly_array]  y_tmp[!, monthly_quarterly_array]]
    q = zeros(4, 1)
    nM = (monthly_quarterly_array .== false) |> sum # number of monthly variables
    nQ = monthly_quarterly_array |> sum # number of quarterly variables
    nBlocks = size(blocks)[2]

    tmp = fill_na(y_tmp)
    y_filled = tmp[:output]
    na_indices = tmp[:na_indices]
    n_rows = nrow(y_filled)
    n_cols = ncol(y_filled)

    res = copy(y_filled)
    res_tmp = copy(res)
    allowmissing!(res_tmp)
    for i in 1:ncol(res_tmp)
        res_tmp[na_indices[!, i], i] .= missing
    end

    pC = 5
    ppC = maximum([p, pC])
    C = missing
    A = missing
    Q = missing
    V0 = missing
    na_indices[1:(pC-1), :] .= true

    for i in 1:nBlocks
        C_i = zeros(n_cols, ppC) # Initialize state variable matrix helper
        idx_i = findall(x-> x==1.0, blocks[!, i]) # Indices of series loading in block i
        idx_iM = idx_i[idx_i .<= nM] # Monthly series indices that load in block i
        idx_iQ = idx_i[idx_i .> nM] # Quarterly series indices that loaded in block i
        cov_matrix = cov(res[!, idx_iM] |> Array) # Eigenvector of cov matrix of monthly data and largest eigenvalue
        d = eigvals(cov_matrix)[end] # Largest eigenvalue...
        v = eigvecs(cov_matrix)[:, end] |> x-> # ...and the associated eigenvector
            sum(x) < 0 ? -x : x # Flip sign for clearer output (not required, but facilitates reading)
        C_i[idx_iM, 1] = v
        f = Array(res[!, idx_iM]) * v # Data projection for eigenvector direction
        for j in 0:(maximum([p+1, pC]) - 1)
            global F = j == 0 ?  f[(pC-j):(length(f)-j)] : hcat(F, f[(pC-j):(length(f)-j)]) # Lag matrix
        end

        ff = F[:, 1:pC] # Projected data with lag structure, so pC-1 fewer entries
        for j in idx_iQ
            xx_j = res_tmp[pC:nrow(res_tmp), j] # For series j, values are dropped to accommodate lag structure
            if sum(ismissing.(xx_j)) < size(ff, 2) + 2
                xx_j = res[pC:nrow(res), j] # Replaces xx_j with spline if too many NaNs
            end
            ff_j = ff[.!ismissing.(xx_j),:]
            xx_j = xx_j[.!ismissing.(xx_j)]
            iff_j = inv(transpose(ff_j) * ff_j)
            Cc = iff_j * transpose(ff_j) * xx_j   # OLS
            Cc = Cc - iff_j * transpose(R_mat) * inv(R_mat * iff_j * transpose(R_mat)) * (R_mat * Cc - q) # Spline data monthly to quarterly conversion
            C_i[j, 1:pC] = transpose(Cc)   # Place in output matrix
        end

        ff = vcat(zeros(pC-1, pC), ff)   # Zeros in first pC-1 entries (replacing dropped from lag)
        res = res .- ff * transpose(C_i)   # Residual calculations
        res_tmp = copy(res)
        allowmissing!(res_tmp)
        for i in 1:ncol(res_tmp)
            res_tmp[na_indices[!, i], i] .= missing
        end
        C = ismissing(C) ? C_i : hcat(C, C_i) # Combine past loadings together

        # transition equation
        z = F[:, 1] # Projected data (no lag)
        Z = F[:, 2:(p+1)] # Data with lag 1
        A_i = zeros(ppC, ppC) # Initialize transition matrix
        A_tmp = inv(transpose(Z) * Z) * transpose(Z) * z # OLS: gives coefficient value AR(p) process
        A_i[1, 1:p] = transpose(A_tmp)
        A_i[2:size(A_i)[1], 1:(ppC-1)] = I(ppC-1)

        Q_i = zeros(ppC, ppC)
        e = z - Z * A_tmp # VAR residuals
        Q_i[1, 1] = cov(e) # VAR covariance matrix
        initV_i = reshape(inv(I(ppC ^ 2) - kron(A_i, A_i)) * reshape(Q_i, (prod(size(Q_i)),1)), (ppC, ppC))
        # Gives top left block for the transition matrix
        A = ismissing(A) ? A_i : blockdiag(sparse(A), sparse(A_i)) |> Array
        Q = ismissing(Q) ? Q_i : blockdiag(sparse(Q), sparse(Q_i)) |> Array
        V0 = ismissing(V0) ? initV_i : blockdiag(sparse(V0), sparse(initV_i)) |> Array
    end

    I_n_cols = I(n_cols) |> x-> # Used inside measurement matrix, only for monthly variables
        x[:, .!monthly_quarterly_array] |>
        Array

    C = hcat(C, I_n_cols)
    C = vcat(zeros(nM, 5 * nQ), transpose(kron(I(nQ), [1, 2, 3, 2, 1]))) |> x-> # Monthly-quarterly agreggation scheme
        hcat(C, x)
    R = res_tmp |> eachcol .|> skipmissing .|> var |> Diagonal # Initialize covariance matrix for transition matrix

    BM = zeros(nM, nM) # Initialize monthly transition matrix values
    SM = zeros(nM, nM) # Initialize monthly residual covariance matrix values

    for i in 1:nM # Loop for monthly variables
        R[i, i] = 1e-4 # Set measurement equation residual covariance matrix diagonal
        res_i = res_tmp[:, i] # Subsetting series residuals for series i
        # Returns number of leading/ending zeros (i.e. missings)
        leadZero = findall(!ismissing, res_i)[1]-1
        endZero = length(res_i) - (findall(!ismissing, res_i)[end])
        # Truncate leading and ending zeros
        res_i = res[!, i] |> x->
            x[(leadZero + 1):(length(res_i) - endZero)]
        # Linear regression: AR(1) process for monthly series residuals
        BM[i, i] = res_i[1:end-1] |> x->
            transpose(inv(transpose(x) * x) * x) * res_i[2:end]
        SM[i, i] = (res_i[2:end] .- res_i[1:end-1] * BM[i, i]) |> var # Residual covariance matrix
    end

    Rdiag = [R[i, i] for i in 1:size(R)[1]]
    σ_e = Rdiag[(nM + 1):n_cols] ./ 19
    Rdiag[(nM+1):n_cols] .= 1e-4
    R = Diagonal(Rdiag) |> collect # Covariance for obs matrix residuals

    # For BQ, SQ
    ρ0 = 0.1
    tmp = zeros(5,5)
    tmp[1, 1] = 1

    # Blocks for covariance matrices
    SQ = kron(Diagonal((1-ρ0^2)*σ_e), tmp)
    BQ = vcat(transpose([ρ0; zeros(4)]), [collect(I(4)) zeros(4)]) |> x->
        kron(collect(I(nQ)), x)
    initViQ = collect(I((5 * nQ)^2)) .- kron(BQ, BQ) |> x->
        inv(x) * reshape(SQ, (prod(size(SQ)),1)) |> x->
        reshape(x, 5*nQ, 5*nQ)
    initViM = collect(I(size(BM)[1])) .- BM^2 |> x->
        1 ./ [x[i, i] for i in 1:size(x)[1]] |>
        Diagonal |> collect |> x->
        x * SM

    # Output
    A = blockdiag(sparse(A), sparse(BM), sparse(BQ)) |> collect # Measurement matrix
    Q = blockdiag(sparse(Q), sparse(SM), sparse(SQ)) |> collect # Residual covariance matrix (transition)
    Z0 = zeros(size(A)[1]) # States
    V0 = blockdiag(sparse(V0), sparse(initViM), sparse(initViQ)) |> collect # Covariance of states

    return Dict(
        :A => A,
        :C => C,
        :Q => Q,
        :R => R,
        :Z0 => Z0,
        :V0 => V0
    )
end


text:
function gyration(df::DataFrame, id::Symbol, x::Symbol, y::Symbol)
    spt.add_dr!(df)
    N = maximum(df.TrackID)
    result = DataFrame(TrackID = Int[], value = Float64[])
    dr = mean(filter(!isnan, df.dr))
    @inbounds for n = 1:maximum(df[!, id])
        data = extract(df, n, id, x, y)
        T = size(data, 1)
        x̄ = mean(data[:, 1])
        ȳ = mean(data[:, 2])
        c = 0.0
        @inbounds for t = 1:T
            x = abs2(data[t, 1] - x̄)
            y = abs2(data[t, 2] - ȳ)
            c += (x + y)
        end
        c *= ((sqrt(π / 2) / dr) * 1 / T)
        push!(result, [n, c])
    end
    result
end

max_stars_repo_name: SosUts/spt

max_issues_repo_name: SosUts/spt

max_forks_repo_name: SosUts/spt


text:
using Gadfly, DataFrames

set_default_plot_size(6inch, 3inch)

points = DataFrame(
    index=[2,10,7,3,10,0,1,10,4,5,4,0,3,7,7,8,5,7,2,3,3,6,4,8,4,10,5,9,5,8],
    val=[6,10,4,2,4,9,1,1,7,4,8,2,2,2,7,9,1,3,7,10,2,3,10,7,8,9,2,6,6,10])
line = DataFrame(val=[4,3,10,6,10,2,3,2,4,7,4], index=collect(0:10))
pointLayer = layer(points, x="index", y="val", Geom.point,Theme(default_color="green"))
lineLayer = layer(line, x="index", y="val", Geom.line)
plot(pointLayer, lineLayer,
     Guide.manual_color_key("Legend", ["Points", "Line"], ["green", "deepskyblue"]))


text:
# SPD-License-Identifier: MIT
#
# This is a thin wrapper over the Knet package to have a more high-level
# API to construct the neural network.

using CUDA
using Knet

const DEFAULT_ARRAY_TYPE = Array{Float32}

trainmode!(model) = model
testmode!(model) = model
reallocate(obj, atype) = deepcopy(obj)

gpu(obj) = reallocate(obj, CuArray{Float32})
cpu(obj) = reallocate(obj, Array{Float32})

"""
    permutedims(perm)

apply a permutation. `perm` is an array of the indices of the
dimension, e.g. `perm = [2, 1, 3]` will permute dimension 2 and 1.
"""
struct PermuteDims
    perm
end
(p::PermuteDims)(X) = permutedims(X, p.perm)

"""
    AddExtraDim(perm)

Add an extra second dimension. This is usually to use 2d-convolution on
1d signals.
"""
struct AddExtraDim end
(::AddExtraDim)(X) = reshape(X, size(X,1), 1, size(X,2), size(X,3))

"""
    RemoveExtraDim(perm)

Remove the extra second dimension. See [`AddExtraDim`](@ref).
"""
struct RemoveExtraDim end
(::RemoveExtraDim)(X) = reshape(X, size(X,1), size(X,3), size(X,4))

"""
    Dropout(pdrop)

Dropout layer. `pdrop` is the probability of a value to be set to zero.
"""
mutable struct Dropout
    pdrop
    training
end
(d::Dropout)(X) = dropout(X, d.pdrop; drop = d.training)
Dropout(pdrop) = Dropout(pdrop, true)
trainmode!(d::Dropout) = d.training = true
testmode!(d::Dropout) = d.training = false

"""
    Dense(in, out, σ = identity)

Dense layer.
"""
struct Dense
    W
    b
    σ
end
function (d::Dense)(X)
    rX = reshape(X, size(X,1), :)
    Y = d.σ.( d.W * rX .+ d.b )
    reshape(Y, :, size(X)[2:end]...)
end

Dense(in::Int, out::Int, σ = identity) =
    Dense(param(out, in, atype = DEFAULT_ARRAY_TYPE),
          param0(out, atype = DEFAULT_ARRAY_TYPE), σ)

reallocate(d::Dense, atype) =
        Dense(Param(atype(value(d.W))), Param(atype(value(d.b))), d.σ)

"""
    Conv(kernelsize, inchannels => outchannels; stride = (1,1), padding = (0,0), dilation=(1,1))

2D-convolution layer.
"""
struct Conv
    W
    b
    stride
    padding
    dilation
end

function Conv(ksize, in_out; stride = 1, padding = 0, dilation = 1)
    W = param(ksize..., in_out.first, in_out.second, atype = DEFAULT_ARRAY_TYPE)
    b = param0(1, 1, in_out.second, 1, atype = DEFAULT_ARRAY_TYPE)
    Conv(W, b, stride, padding, dilation)
end

(c::Conv)(X) = conv4(c.W, X, padding = c.padding, stride = c.stride,
                     dilation = c.dilation) .+ c.b

function reallocate(c::Conv, atype)
    Conv(Param(atype(value(c.W))), Param(atype(value(c.b))), c.stride,
         c.padding, c.dilation)
end

"""
    BatchNorm(dim, σ = identity)

Batch normalization.
"""
struct BatchNorm
    moments
    params
    σ
end
BatchNorm(dim, σ = identity) = BatchNorm(bnmoments(),
                                         DEFAULT_ARRAY_TYPE(bnparams(dim)), σ)
(bn::BatchNorm)(X) = bn.σ.(batchnorm(X, bn.moments, bn.params))

function reallocate(bn::BatchNorm, atype)
    mean, var = bn.moments.mean, bn.moments.var
    BatchNorm(
        bnmoments(;
            momentum = bn.moments.momentum,
            mean = isnothing(mean) ? nothing : atype(mean),
            var = isnothing(var) ? nothing : atype(var),
            meaninit = bn.moments.meaninit,
            varinit = bn.moments.varinit
        ),
        atype(value(bn.params)),
        bn.σ
    )
end

"""
    chain(layer1, layer2, ...)

Build a sequential model.
"""
struct Chain
    layers
    Chain(layers...) = new(layers)
end

function (c::Chain)(X)
    for layer in c.layers
        X = layer(X)
    end
    X
end

reallocate(chain::Chain, atype) = Chain((reallocate.(chain.layers, atype))...)

function Base.show(io::IO, c::Chain)
    println(io, "Chain(")
    for layer in c.layers
        println(layer)
    end
    print(io, ")")
end

max_stars_repo_name: MartinKocour/SpeechLab

max_issues_repo_name: MartinKocour/SpeechLab

max_forks_repo_name: MartinKocour/SpeechLab


text:
#!/usr/bin/env julia

using DataFrames
using Measurements
using DataStructures

mutable struct Kernel{T<:AbstractFloat}
    name::String
    registers::Int
    runtime::T
end
runtime(kernel::Kernel) = kernel.runtime
registers(kernel::Kernel) = kernel.registers

mutable struct Iteration{T<:AbstractFloat}
    runtime::T

    # kernels
    kernels::Vector{Kernel{T}}

    # memory transfers
    transfer_size::Float64
    transfer_count::Int
    transfer_runtime::T

    # API calls
    api_count::Int
    api_runtime::T

    Iteration{T}() where {T} = new{T}(zero(T), Kernel{T}[], 0.0, 0, zero(T), 0, zero(T))
end

function Base.show(io::IO, it::Iteration)
    println(io, "Iteration taking $(round(it.runtime,2)) us:")
    println(io, " - $(length(it.kernels)) kernel launches: $(round(sum(runtime, it.kernels), 2)) us")
    println(io, " - $(it.transfer_count) memory transfers: $(round(it.transfer_size, 2)) MB in $(round(it.transfer_runtime, 2)) us")
    print(io, " - $(it.api_count) API calls: $(round(it.api_runtime, 2)) us")
end

function group_trace(table)
    its = Iteration[]
    it = nothing
    it_start = 0.
    for row in eachrow(table)
        # handle start/stop of a new iteration
        if contains(row[:Name], "[Range start]")
            @assert it == nothing
            it = Iteration{Float64}()
            it_start = row[:Start]
            continue
        end
        if contains(row[:Name], "[Marker]")
            continue
        end
        @assert it != nothing
        if contains(row[:Name], "[Range end]")
            it_stop = row[:Start]
            it.runtime = it_stop - it_start
            push!(its, it)
            it = nothing
            continue
        end

        # other entries are either API traces, memory transfers or executed kernels
        if startswith(row[:Name], "cu")
            it.api_count += 1
            # TODO: only count API calls that contribute to the total execution (ie. no
            # async calls, event queries, etc)
            if row[:Name] != "cuCtxSynchronize"
                it.api_runtime += row[:Duration]
            end
        elseif contains(row[:Name], "[CUDA memcpy")
            it.transfer_size += row[:Size]
            it.transfer_count += 1
            it.transfer_runtime += row[:Duration]
        else
            push!(it.kernels, Kernel{Float64}(row[:Name], row[:Registers_Per_Thread], row[:Duration]))
        end
    end

    return its
end

function average_trace(its)
    # NOTE: better just chuck everything in the DataFrame,
    #       and group+average when displaying

    popfirst!(its)

    avg_it = Iteration{Measurement{Float64}}()

    # kernels
    @assert all(it->length(its[1].kernels)==length(it.kernels), its)
    for i in 1:length(its[1].kernels)
        @assert all(it->its[1].kernels[i].name==it.kernels[i].name, its)
        @assert all(it->its[1].kernels[i].registers==it.kernels[i].registers, its)
        push!(avg_it.kernels,
            Kernel{Measurement{Float64}}(
                its[1].kernels[i].name,
                its[1].kernels[i].registers,
                mean([it.kernels[i].runtime for it in its]) ±
                 std([it.kernels[i].runtime for it in its]),
            )
        )
    end

    # other
    ## counters
    for field in (:transfer_size, :transfer_count, :api_count)
        @assert all(it->getfield(its[1], field) == getfield(it, field), its)
        setfield!(avg_it, field, getfield(its[1], field))
    end
    ## timings
    for field in (:runtime, :transfer_runtime, :api_runtime)
        val = mean([getfield(it, field) for it in its]) ± std([getfield(it, field) for it in its])
        setfield!(avg_it, field, val)
    end

    return avg_it
end

function read_trace(tag)
    df = open("$(tag)_trace.jls", "r") do io
        deserialize(io)
    end

    its = group_trace(df)
    return average_trace(its)
end

function read_metrics(tag)
    df = open("$(tag)_metrics.jls", "r") do io
        deserialize(io)
    end

    return df
end

function read(tag)
    return read_metrics(tag), read_trace(tag)
end

function get_metric(df, kernel, metric, col) # col: Min, Max or Avg
    mask = (df[:Kernel] .== kernel) .& (df[:Metric_Name] .== metric)
    any(mask) || return missing
    val = first(eachrow(df[mask, :]))[col]
    if isa(val, String)
        if endswith(val, "%")
            val = val[1:end-1]
            return parse(Float64, val) / 100
        else
            return parse(Float64, val)
        end
    else
        return val
    end
end

# basic timings

function add_row!(df, trace, metrics; language, fused=missing, control=missing, arity=missing, problem_size)
    push!(df, [language, fused, control, arity, problem_size,
               trace.runtime,
               length(trace.kernels), maximum(registers.(trace.kernels)),
               sum(runtime, trace.kernels),
               mean(kernel->get_metric(metrics, kernel.name, "achieved_occupancy",                  :Avg), trace.kernels),
               mean(kernel->get_metric(metrics, kernel.name, "branch_efficiency",                   :Avg), trace.kernels),
               mean(kernel->get_metric(metrics, kernel.name, "warp_execution_efficiency",           :Avg), trace.kernels),
               mean(kernel->get_metric(metrics, kernel.name, "warp_nonpred_execution_efficiency",   :Avg), trace.kernels),
               trace.transfer_size, trace.transfer_count, trace.transfer_runtime,
               trace.api_count, trace.api_runtime])
end

function process(dir)
    df = DataFrame(language = Symbol[],
                   # kernel type
                   fused=Union{Missing,Bool}[], control=Union{Missing,Symbol}[],
                   arity=Union{Missing,Int}[],
                   # parameters
                   problem_size=Int[],
                   # total
                   runtime = AbstractFloat[],
                   # only kernel
                   kernel_count = Int[], kernel_registers = Int[],
                   kernel_runtime = AbstractFloat[],
                   ## metrics
                   kernel_occupancy = Union{Missing,AbstractFloat}[],
                   kernel_branch_efficiency = Union{Missing,AbstractFloat}[],
                   kernel_warp_efficiency = Union{Missing,AbstractFloat}[],
                   kernel_warp_nonpred_efficiency = Union{Missing,AbstractFloat}[],
                   # only memory
                   transfer_size = Float64[], transfer_count = Int[], transfer_runtime = AbstractFloat[],
                   # only api
                   api_count = Int[], api_runtime = AbstractFloat[])

    cd(dir) do
        for problem_size in [512,1024,2048]
            try
                metrics, trace = read("python_$(problem_size)")
                add_row!(df, trace, metrics; problem_size=problem_size, language=:python, fused=false)
            catch exception
                isa(exception, SystemError) || rethrow()
            end

            for tfstyle in [true, false], uniform in [true, false]
                try
                    metrics, trace = read("julia_$(tfstyle ? "tfstyle" : "fused")_$(uniform ? "uniform" : "random")_$(problem_size)")
                    add_row!(df, trace, metrics; problem_size=problem_size, language=:julia,
                             fused=!tfstyle, control=(uniform ? :uniform : :random))
                catch exception
                    isa(exception, SystemError) || rethrow()
                end
            end

            for arity in 1:10
                try
                    metrics, trace = read("julia_arity$(arity)_$(problem_size)")
                    add_row!(df, trace, metrics; problem_size=problem_size, language=:julia, arity=arity)
                catch exception
                    isa(exception, SystemError) || rethrow()
                end
            end
        end
    end

    df
end

# main

"""Expand columns of type Measurement to separate value and error columns for CSV export."""
function expand_measurements!(df)
    for col in names(df)
        # Measurement <: AbstractFloat, so check actual values
        if any(val -> eltype(val) <: Measurement, df[col])
            df[Symbol("$(col)_err")] = Measurements.uncertainty.(df[col])
            df[col]                  = Measurements.value.(df[col])
        end
    end
end

function process_all()
    measurements = Dict(
        :titan   => joinpath(@__DIR__, "GeForce GTX TITAN"),
        :p100    => joinpath(@__DIR__, "Tesla P100-PCIE-16GB"),
        :v100    => joinpath(@__DIR__, "Tesla V100-PCIE-16GB"),
        # :gtx1080 => joinpath(@__DIR__, "GeForce GTX 1080 Ti"),
    )

    output = joinpath(dirname(@__DIR__), "img")

    df = nothing
    for (gpu, path) in measurements
        df_gpu = process(path)
        df_gpu[:gpu] = gpu
        df_gpu = df_gpu[[end, 1:end-1...]]
        if df === nothing
            df = df_gpu
        else
            append!(df, df_gpu)
        end
    end

    # basic HLMSTM timings

    let df = filter(row -> ismissing(row[:arity]) && row[:control] !== :uniform, df)
        delete!(df, [:arity, :control])

        expand_measurements!(df)
        writetable(joinpath(output, "hmlstm.csv"), df)
    end

    # arity scaling

    let df = filter(row -> !ismissing(row[:arity]), df)
        delete!(df, [:language, :control, :fused])

        expand_measurements!(df)
        writetable(joinpath(output, "arity.csv"), df)
    end

    # divergence analysis

    let df = filter(row -> !ismissing(row[:control]) && row[:fused], df)
        delete!(df, [:arity, :language, :fused])

        join_columns = [:gpu, :problem_size]
        other_columns = setdiff(names(df), [join_columns; :control])

        # split columns into separate ones based on the type of control

        control_dfs = DataFrame[]

        for control in (:random, :uniform)
            df_control = df[df[:control] .== control, :]
            delete!(df_control, :control)

            for col in other_columns
                rename!(df_control, col => Symbol("$(control)_$(col)"))
            end

            push!(control_dfs, df_control)
        end

        # join and divide control columns to get ratios

        df = join(control_dfs..., on=join_columns)

        for col in other_columns
            random = Symbol("random_$(col)")
            uniform = Symbol("uniform_$(col)")
            ratio = Symbol("$(col)_ratio")

            df[ratio] = df[random] ./ df[uniform]
            delete!(df, [random, uniform])
        end

        # drop all-one columns
        for col in names(df)
            if all(val -> val === 1.0, df[col])
                delete!(df, col)
            end
        end

        expand_measurements!(df)
        writetable(joinpath(output, "divergence.csv"), df)
    end
end


text:
"""
A package for point process modeling, simulation and inference.
"""
module PointProcesses

## Imports

using DataStructures
using GalacticOptim
using LinearAlgebra
using LogExpFunctions
using MeasureTheory
using NamedTupleTools
using OffsetArrays
using Optim
using Plots
using Quadrature
using Random
using SciMLBase
using StatsPlots
using TransformVariables

# Hidden names

using Random: GLOBAL_RNG

# Functions to extend

import Base: eltype, length, rand
import Distributions: fit, fit_mle, suffstats
import MeasureTheory: density, logdensity, sampletype, testvalue

## Includes and exports

# Reexports

export eltype, rand  # Base
export mean
export logpdf, fit, fit_mle, suffstats  # Distributions
export plot, scatter  # Plot

# History

include("history/abstract.jl")
export AbstractHistory

include("history/temporal.jl")
export TemporalHistory
export event_times, event_marks, min_time, max_time
export nb_events, has_events, duration, time_change

# Markov processes

include("markov/abstract.jl")
export AbstractMarkovChain
export nb_states
export AbstractMarkovChainPrior

include("markov/discrete_time.jl")
export DiscreteMarkovChain
export initial_distribution, transition_matrix, stationary_distribution
export DiscreteMarkovChainPrior, DiscreteMarkovChainStats

include("markov/continuous_time.jl")
export ContinuousMarkovChain
export rate_matrix, rate_diag, discretize_chain
export ContinuousMarkovChainPrior, ContinuousMarkovChainStats

# Point processes

include("point_processes/abstract.jl")
export AbstractPointProcess, build_transform

include("point_processes/temporal.jl")
export TemporalPointProcess, BoundedTemporalPointProcess
export intensity, mark_distribution, ground_intensity, ground_intensity_bound
export MultivariateTemporalPointProcess, all_marks, all_mark_probabilities
export integrated_ground_intensity, check_residuals

# Models

include("models/poisson.jl")
export PoissonProcess

include("models/poisson_multivariate_naive.jl")
export NaiveMultivariatePoissonProcess

# Hidden Markov models

include("hmm/hmm.jl")
export HiddenMarkovModel
export transitions, emissions, emission

include("hmm/baum_welch.jl")
export forward_nolog!, forward_log!
export backward_nolog!, backward_log!
export forward_backward_nolog!, forward_backward_log!
export update_obs_pdf!, update_obs_logpdf!
export baum_welch!, baum_welch

include("hmm/mmpp.jl")
export MarkovModulatedPoissonProcess

include("hmm/ryden.jl")
export forward_backward, ryden

# Utils

include("utils/overflow.jl")
export all_minus_inf, all_plus_inf, all_zeros, all_nan

include("utils/dists.jl")

include("utils/categorical.jl")

include("utils/plot.jl")
export plot_events, plot_intensity, qqplot_interevent_times

include("utils/randvals.jl")
export uniformprobvec, randprobvec
export uniformtransmat, randtransmat

end


text:
immutable NewtonSolver{T, TF, TJ, TL} <: AbstractNewtonSolver{T}
    @newton_solver_variables
    function NewtonSolver(x, Fparams, Jparams, linear_solver, nmax, atol, rtol, stol)
        new(x, Fparams, Jparams, linear_solver, NonlinearSolverParameters{T}(nmax, atol, rtol, stol), NonlinearSolverStatus{T}())
    end
end

function NewtonSolver(x::Vector, Fparams::NonlinearFunctionParameters; J=nothing, linear_solver=nothing, nmax=DEFAULT_nmax, atol=DEFAULT_atol, rtol=DEFAULT_rtol, stol=DEFAULT_stol, ϵ=DEFAULT_ϵ, autodiff=false)
    T = eltype(x)
    n = length(x)
    Jparams = getJacobianParameters(J, Fparams, ϵ, T, n, autodiff)
    linear_solver = getLinearSolver(T, n, linear_solver)
    NewtonSolver{T, typeof(Fparams), typeof(Jparams), typeof(linear_solver)}(x, Fparams, Jparams, linear_solver, nmax, atol, rtol, stol)
end

function solve!{T}(s::NewtonSolver{T})
    function_stages!(s.x, s.linear.b, s.Fparams)
    scale!(s.linear.b, -one(T))
    s.status.i  = 0
    s.status.rₐ = residual_absolute(s.linear.b)
    s.status.r₀ = s.status.rₐ

    if s.status.r₀ ≥ s.params.atol²
        for s.status.i = 1:s.params.nmax
            computeJacobian(s.x, s.linear.A, s.Jparams)
            factorize!(s.linear)
            solve!(s.linear)
            simd_xpy!(s.linear.b, s.x)
            s.status.rᵣ = residual_relative(s.linear.b, s.x)
            function_stages!(s.x, s.linear.b, s.Fparams)
            scale!(s.linear.b, -one(T))
            s.status.rₛ = s.status.rₐ
            s.status.rₐ = residual_absolute(s.linear.b)
            s.status.rₛ = abs(s.status.rₛ - s.status.rₐ)/s.status.r₀

            if s.status.rₐ < s.params.atol² || s.status.rᵣ < s.params.rtol || s.status.rₛ < s.params.stol²
                break
            end
        end
    end
    nothing
end


text:
using PyCall
# Method One
py"[1,2]"
py"""
import numpy as np

def sinpi(x):
    return np.sin(np.pi * x)
"""
py"sinpi"(1)

# Method Two: autoconverte the object and can easily import python modules.
np = pyimport("numpy")
pd = pyimport("pandas")

A = np.array([[1, 2], [3, 4]])
B = [1 2; 3 4]
B == A

torch = pyimport("torch")
torch.cuda.is_available()


text:
#=

  Euler #22 in Julia.

  Problem 22:
  """
  Using names.txt (right click and 'Save Link/Target As...'), a 46K
  text file containing over five-thousand first names, begin by sorting
  it into alphabetical order. Then working out the alphabetical value
  for each name, multiply this value by its alphabetical position in the
  list to obtain a name score.

  For example, when the list is sorted into alphabetical order, COLIN,
  which is worth 3 + 15 + 12 + 9 + 14 = 53, is the 938th name in
  the list. So, COLIN would obtain a score of 938 53 = 49714.

  What is the total of all the name scores in the file?")
  """

  This Julia program was created by Hakan Kjellerstrand, hakank@gmail.com
  See also my Julia page: http://www.hakank.org/julia/

=#

include("Euler.jl")

function to_code(s)
    return [Int(s[i])-64 for i in 1:length(s)]
end

# 0.00231866s
function euler22a()
    words = sort(split(readline("euler22_names.txt"),","))
    a = 1;
    s = 0;
    for word in words
        word = replace(word,"\""=>"")
        s += a*sum(to_code(word));
        a += 1;
    end
    return s;

end

# 0.02438064s
function euler22b()
    a = 1
    return sum(map(word->(word=replace(word,"\""=>"");a=a+1; (a-1)*sum(to_code(word))
                   ),
                  sort(split(readline("euler22_names.txt"),","))))

end

# Nicer chaining.
# Note the () before |>sum
# 0.02813539s
function euler22c()
    a = 1
    return (sort(split(readline("euler22_names.txt"),",")).|>word->(
                    word=replace(word,"\""=>"");
                    a=a+1; (a-1)*sum(to_code(word))
                    ))|>sum
end

run_euler(euler22a)
# run_euler(euler22b)
# run_euler(euler22c)

max_stars_repo_name: Wikunia/hakank

max_issues_repo_name: Wikunia/hakank

max_forks_repo_name: Wikunia/hakank


text:
@doc raw"""
    Ellipsoid(a, b, f, finv, e2, ϵ2, name)
    Ellipsoid(a, finv, name)

Earth ellipsoid model.

# Fields

- `a::Real`: semi-major axis (m).
- `b::Real`: semi-minor axis (m).
- `f::Real`: flattening (``f``).
- `finv::Real`: inverse of flattening (``f^{-1}``).
- `e2::Real`: eccentricity squared.
- `ϵ2::Real`: second eccentricity squared.
- `name::String`: ellipsoid name.

# Notes

``a`` is the semi-major axis, ``b`` is the semi-minor axis. Ellipsoids are normally defined
by ``a`` and ``f^{-1}``. The rest of parameters are derived.

- Flattening:
```math
f = 1 - \frac{b}{a}
```
- Eccentricity (first eccentricity):
```math
e^2 = \frac{a^2 - b^2}{a^2} = f (2 - f)
```
- Second eccentricity:
```math
ϵ^2 = \frac{a^2 - b^2}{b^2} = \frac{f (2 - f)}{(1 - f)^2}
```

# References

1. Stevens, B. L., Lewis, F. L., (1992). Aircraft control and simulation: dynamics, controls
   design, and autonomous systems. John Wiley & Sons. Section 1.6 (page 23)
2. Rogers, R. M. (2007). Applied mathematics in integrated navigation systems. American
   Institute of Aeronautics and Astronautics. Chapter 4 (Nomenclature differs from 1).
3. Bowring, B. R. (1976). Transformation from spatial to geographical coordinates. Survey
   review, 23(181), 323-327.
"""
struct Ellipsoid
    "Semi-major axis (m)."
    a::Real
    "Semi-minor axis (m)."
    b::Real
    "Flattening ``f = \frac{a-b}{a} = a - \frac{b}{a}``."
    f::Real
    "Inverse of flattening."
    finv::Real
    "Eccentricity. ``e^2 = \frac{a^2 - b^2}{a^2} = f (2 - f)``"
    e2::Real
    "Second eccentricity. ``ϵ^2 = \frac{a^2 - b^2}{b^2} = \frac{f (2 - f)}{(1 - f)^2}``."
    ϵ2::Real
    "Ellipsoid name."
    name::String
end

function Ellipsoid(a, finv, name)
    f = 1 / finv
    b = a * (1 - f)
    e2 = f * (2 - f)
    ϵ2 = e2 / (1 - f)^2
    return Ellipsoid(a, b, f, finv, e2, ϵ2, name)
end

# Rogers, R. M. (2007). Applied mathematics in integrated navigation systems.
# American Institute of Aeronautics and Astronautics. (Page 76, table 4.1)
const Clarke1866 = Ellipsoid(6378206.4, 294.9786982, "Clarke1866")
const Clarke1880 = Ellipsoid(6378249.145, 294.465, "Clarke1880")
const International = Ellipsoid(6378388.0, 297.0, "International")
const Bessel = Ellipsoid(6377397.155, 299.1528128, "Bessel")
const Everest = Ellipsoid(6377276.345, 300.8017, "Everest")
const ModifiedEverest = Ellipsoid(6377304.063, 300.8017, "ModifiedEverest")
const AustralianNational = Ellipsoid(6378160.0, 298.25, "AustralianNational")
const SouthAmerican1969 = Ellipsoid(6378160.0, 298.25, "SouthAmerican1969")
const Airy = Ellipsoid(6377564.396, 299.3249646, "Airy")
const ModifiedAiry = Ellipsoid(6377340.189, 299.3249646, "ModifiedAiry")
const Hough = Ellipsoid(6378270.0, 297.0, "Hough")
const Fischer1960SouthAsia = Ellipsoid(6378155.0, 298.3, "Fischer1960SouthAsia")
const Fischer1960Mercury = Ellipsoid(6378166.0, 298.3, "Fischer1960Mercury")
const Fischer1968 = Ellipsoid(6378150.0, 298.3, "Fischer1968")
const WGS60 = Ellipsoid(6378165.0, 298.3, "WGS60")
const WGS66 = Ellipsoid(6378145.0, 298.25, "WGS66")
const WGS72 = Ellipsoid(6378135.0, 298.26, "WGS72")
const WGS84 = Ellipsoid(6378137.0, 298.257223563, "WGS84")

"""
    llh2ecef(lat, lon, height; ellipsoid=WGS84)

Transform geodetic latitude, longitude (rad) and ellipsoidal height (m) to ECEF for the
given ellipsoid (default ellipsoid is WGS84).

# References

1. Rogers, R. M. (2007). Applied mathematics in integrated navigation systems. American
   Institute of Aeronautics and Astronautics. (Page 75, equations 4.20, 4.21, 4.22)
"""
function llh2ecef(lat, lon, height; ellipsoid=WGS84)
    f = ellipsoid.f
    a = ellipsoid.a
    e2 = ellipsoid.e2

    var = a / sqrt(1.0 - e2 * sin(lat)*sin(lat))

    px = (var + height) * cos(lat)*cos(lon)
    py = (var + height) * cos(lat)*sin(lon)
    pz = (var * (1.0 - e2) + height)* sin(lat)

    return [px, py, pz]
end

"""
    ecef2llh(xecef, yecef, zecef; ellipsoid=WGS84)

Transform ECEF coordinates to geodetic latitude, longitude (rad) and ellipsoidal height (m)
for the given ellipsoid (default ellipsoid is WGS84).

# Notes

* The transformation is direct without iterations as [1] introduced the need to iterate for
  near Earth positions.
* [2] is an update of increased accuracy of [1]. The former is used in this implementation
  although the latter implementation is commented in the code.
* Model becomes unstable if latitude is close to 90º. An alternative equation
  can be found in [2] equation (16) but has not been implemented.

 # References

1. Bowring, B. R. (1976). Transformation from spatial to geographical coordinates. Survey
   review, 23(181), 323-327.
2. Bowring, B. R. (1985). The accuracy of geodetic latitude and height equations. Survey
   Review, 28(218), 202-206.
"""
function ecef2llh(xecef, yecef, zecef; ellipsoid=WGS84)

    x, y, z = xecef, yecef, zecef
    e = sqrt(ellipsoid.e2);
    e2 = ellipsoid.e2
    ϵ2 = ellipsoid.ϵ2
    a = ellipsoid.a
    b = ellipsoid.b

    p = sqrt(x*x + y*y)
    R = sqrt(p*p + z*z)
    θ = atan(z, p)
    # [1] equation (1) does not change in [2]
    lon = atan(y, x)
    # u -> geographical latitude
    # [1] below equation (4) and [2] equation (6)
    # This lead to errors in latitud with maximum value 0.0018"
    #u = atan(a/b * z/x)

    # [2] equation (17) If the latitude is also required to be very accurate
    # for outer-space positions then the value of u for (6) should be obtained
    # from:
    u = atan(b*z / (a*p) * (1 + ϵ2 * b / R))
    # [1] equation (4)
    #lat = atan((z + ϵ2 * b * sin(u)^3) / (x - e2 * a * cos(u)^3))
    # [2] equation (18)
    lat = atan((z + ϵ2 * b * sin(u)^3) / (p - e2 * a * cos(u)^3))
    # [2] after equation (1)
    v = a / sqrt(1.0 - e2*sin(lat)^2)
    # [2] equation (7)
    # height = p*cos(lat) + z*sin(lat) - a*a / v
    # equivalent to [2] equation (8)
    height = R * cos(lat - θ) - a*a / v
    # which is insensitive to the error in latitude calculation (The influence
    # is of order 2 [2] equation (12))
    return [lat, lon, height]
end


text:
using ISA #References: Sandoval, Steven, and Phillip L. De Leon. "The Instantaneous Spectrum: A General Framework for Time-Frequency Analysis." IEEE Transactions on Signal Processing 66.21 (2018): 5679-5693.

using Plots

function main()
#DEFINE COMPONENT
a₀(t) = exp(-t^2)
ω₀(t) = 25
φ₀ = 0
𝐶₀ = Tuple([a₀,ω₀,φ₀])
ψ₀ = AMFMcomp(𝐶₀)

#DEFINE COMPONENT
a₁(t) = exp(-0.5t^2)
ω₁(t) = 10 + 3*sin(t)
φ₁ = 0
𝐶₁ = Tuple([a₁,ω₁,φ₁])
ψ₁ = AMFMcomp(𝐶₁)

t = -1.0:0.01:1.0

p1 = isaPlot3d(ψ₀,t)
p2 = isaPlot3d(AMFMmodel([ψ₀,ψ₁]),t)

#COMPONENT OBSERVATION
Ψ₀ = ψ₀(t)
Ψ₁ = ψ₁(t)

#ESTIMATE NUMERICAL COMPONENT
𝚿₀ = AMFMdemod(Ψ₀,t)
𝚿₁ = AMFMdemod(Ψ₁,t)

p3 = isaPlot3d(𝚿₀)
p4 = isaPlot3d([𝚿₀,𝚿₁])
plot(p1,p2,p3,p4)

end

main()

max_stars_repo_name: kricheso/ISA

max_issues_repo_name: kricheso/ISA

max_forks_repo_name: kricheso/ISA


text:
using SoapySDR

# Don't forget to add/import a device-specific plugin package!
# using xtrx_jll
# using SoapyLMS7_jll
# using SoapyRTLSDR_jll
# using SoapyPlutoSDR_jll
# using SoapyUHD_jll

# Get the Devices
devs = Devices()

# we can add arguments to the device constructor or filter multiple devices:
# devs[1]["refclk"] = "26000000"
# devs = filter(x -> haskey(x, "driver") && x["driver"] == "plutosdr", devs)

# Now open the first device and get the first RX and TX channel streams
dev = open(devs[1])
c_tx = dev.tx[1]
c_rx = dev.rx[1]

# Configure channel with appropriate parameters.  We choose to sneak into
# the very high end of the 2.4 GHz ISM band, above Wifi channel 14 (which
# itself is well outside of the typical US equipment range, which only goes
# up to channel 11, e.g. 2.473 GHz).
c_tx.bandwidth = 2u"MHz"
c_rx.bandwidth = 200u"kHz"
c_tx.frequency = 2.498u"GHz"
c_rx.frequency = 2.498u"GHz"
c_tx.gain = 52u"dB"
c_rx.gain = -2u"dB"
c_tx.sample_rate = 1u"MHz"
c_rx.sample_rate = 1u"MHz"

# Write out a sinusoid oscillating at 100KHz, lasting 10ms
t = (1:round(Int, 0.01*1e6))./1e6
data_tx = ComplexF32.(sin.(2π.*t.*100e3), 0.0f0)
data_tx_zeros = zeros(ComplexF32, length(data_tx))

# Open both RX and TX streams
s_tx = SoapySDR.Stream(ComplexF32, [c_tx])
s_rx = SoapySDR.Stream(ComplexF32, [c_rx])

# Then, do the actual writing and reading!
function loopback_test(s_tx, s_rx, num_buffers)
    # allocate some recieve buffers
    data_rx_buffs = Vector{ComplexF32}[]
    for _ in 1:num_buffers
        push!(data_rx_buffs, Vector(undef, length(data_tx)))
    end

    SoapySDR.activate!.((s_tx, s_rx))
    # Read/write `num_buffers`, writing zeros out except for one buffer.
    for idx in 1:num_buffers
        if idx == ceil(num_buffers/2)
            Base.write(s_tx, (data_tx,))
        else
            Base.write(s_tx, (data_tx_zeros,))
        end
        Base.read!(s_rx, (data_rx_buffs[idx], ))
    end
    SoapySDR.deactivate!.((s_tx, s_rx))
    return data_rx_buffs
end

# This should take about 100ms, since we're rx/tx'ing 10x buffers which should each be 10ms long.
loopback_test(s_tx, s_rx, 2)
data_rx_buffs = loopback_test(s_tx, s_rx, 3)

# Make sure we close the Streams and Devices where we are done with them
close.((s_tx, s_rx))
close(dev)

# Join all buffers together
data_rx = vcat(data_rx_buffs...)

# Should see a nice big spike of a sinusoid being transmitted in the middle
using Plots; plotly(size=(750,750)); plot(real.(data_rx))


text:
@testset "Inverter Components" begin
    converter = AverageConverter(690.0, 2750000.0) #S_rated goes in Watts
    @test converter isa PowerSystems.DynamicComponent
    dc_source = FixedDCSource(600.0) #Not in the original data, guessed.
    @test dc_source isa PowerSystems.DynamicComponent
    filter = LCLFilter(0.08, 0.003, 0.074, 0.2, 0.01)
    @test filter isa PowerSystems.DynamicComponent
    pll = KauraPLL(500.0, 0.084, 4.69)
    @test pll isa PowerSystems.DynamicComponent
    reduced_pll = ReducedOrderPLL(500.0, 0.084, 4.69)
    @test reduced_pll isa PowerSystems.DynamicComponent
    virtual_H = VirtualInertia(2.0, 400.0, 20.0, 2 * pi * 50.0)
    @test virtual_H isa PowerSystems.DeviceParameter
    P_control = ActivePowerDroop(0.2, 1000.0)
    @test P_control isa PowerSystems.DeviceParameter
    P_control_PI = ActivePowerPI(2.0, 20.0, 50.0)
    @test P_control_PI isa PowerSystems.DeviceParameter
    Q_control = ReactivePowerDroop(0.2, 1000.0)
    @test Q_control isa PowerSystems.DeviceParameter
    Q_control_PI = ReactivePowerPI(2.0, 20.0, 50.0)
    @test Q_control_PI isa PowerSystems.DeviceParameter
    outer_control = OuterControl(virtual_H, Q_control)
    @test outer_control isa PowerSystems.DynamicComponent
    test_accessors(outer_control)
    outer_control_droop = OuterControl(P_control, Q_control)
    @test outer_control_droop isa PowerSystems.DynamicComponent
    test_accessors(outer_control_droop)
    outer_control_PI = OuterControl(P_control_PI, Q_control_PI)
    @test outer_control_PI isa PowerSystems.DynamicComponent
    test_accessors(outer_control_PI)
    vsc = VoltageModeControl(0.59, 736.0, 0.0, 0.0, 0.2, 1.27, 14.3, 0.0, 50.0, 0.2)
    @test vsc isa PowerSystems.DynamicComponent
    vsc2 = VoltageModeControl(0.59, 736.0, 0.0, 0.0, 0.2, 1.27, 14.3, 0.0, 50.0, 0.2)
    @test vsc2 isa PowerSystems.DynamicComponent
    vsc3 = CurrentModeControl(1.27, 14.3, 0.0)
    @test vsc3 isa PowerSystems.DynamicComponent
    BESS_source = ZeroOrderBESS(
        (sqrt(8) / sqrt(3)) * 690.0,
        (sqrt(3) / sqrt(8)) * 2750000.0,
        370.0,
        0.001,
        4.63,
        3200.0,
        0.6,
        4.0,
        0.39,
        10.34,
        1.08,
    )
    @test BESS_source isa PowerSystems.DynamicComponent
end

@testset "Dynamic Inverter" begin
    sys = create_system_with_dynamic_inverter()
    inverters = collect(get_components(DynamicInverter, sys))
    @test length(inverters) == 1
    test_inverter = inverters[1]
    @test test_inverter isa PowerSystems.Component
    test_accessors(test_inverter)
end

@testset "Generic Renewable Models" begin
    converter_regca1 = RenewableEnergyConverterTypeA(
        T_g = 0.02,
        Rrpwr = 10.0,
        Brkpt = 0.9,
        Zerox = 0.4,
        Lvpl1 = 1.22,
        Vo_lim = 1.2,
        Lv_pnts = (0.5, 0.9),
        Io_lim = -1.3,
        T_fltr = 0.2,
        K_hv = 0.0,
        Iqr_lims = (-100.0, 100.0),
        Accel = 0.7,
        Lvpl_sw = 0,
    )
    @test converter_regca1 isa PowerSystems.DynamicComponent
    filt_current = RLFilter(rf = 0.0, lf = 0.1)
    @test filt_current isa PowerSystems.DynamicComponent
    inner_ctrl_typeB = RECurrentControlB(
        Q_Flag = 0,
        PQ_Flag = 0,
        Vdip_lim = (-99.0, 99.0),
        T_rv = 0.02,
        dbd_pnts = (-1.0, 1.0),
        K_qv = 0.0,
        Iqinj_lim = (-1.1, 1.1),
        V_ref0 = 0.0,
        K_vp = 10.0,
        K_vi = 60.0,
        T_iq = 0.02,
        I_max = 1.11,
    )
    @test inner_ctrl_typeB isa PowerSystems.DynamicComponent
    # Creates 2^5 = 32 combinations of flags for an outer control
    for (F_flag, VC_flag, R_flag, PF_flag, V_flag) in
        reverse.(Iterators.product(fill(0:1, 5)...))[:]
        P_control_typeAB = ActiveRenewableControllerAB(
            bus_control = 0,
            from_branch_control = 0,
            to_branch_control = 0,
            branch_id_control = "0",
            Freq_Flag = F_flag,
            K_pg = 1.0,
            K_ig = 0.05,
            T_p = 0.25,
            fdbd_pnts = (-1.0, 1.0),
            fe_lim = (-99.0, 99.0),
            P_lim = (0.0, 1.2),
            T_g = 0.1,
            D_dn = 0.0,
            D_up = 0.0,
            dP_lim = (-99.0, 99.0),
            P_lim_inner = (0.0, 1.2),
            T_pord = 0.02,
        )
        Q_control_typeAB = ReactiveRenewableControllerAB(
            bus_control = 0,
            from_branch_control = 0,
            to_branch_control = 0,
            branch_id_control = "0",
            VC_Flag = VC_flag,
            Ref_Flag = R_flag,
            PF_Flag = PF_flag,
            V_Flag = V_flag,
            T_fltr = 0.02,
            K_p = 18.0,
            K_i = 5.0,
            T_ft = 0.0,
            T_fv = 0.05,
            V_frz = 0.0,
            R_c = 0.0,
            X_c = 0.0,
            K_c = 0.0,
            e_lim = (-0.1, 0.1),
            dbd_pnts = (-1.0, 1.0),
            Q_lim = (-0.43, 0.43),
            T_p = 0.0,
            Q_lim_inner = (-0.44, 0.44),
            V_lim = (0.9, 1.05),
            K_qp = 0.0,
            K_qi = 0.01,
        )
        @test P_control_typeAB isa PowerSystems.DeviceParameter
        @test Q_control_typeAB isa PowerSystems.DeviceParameter
        outer_control_typeAB = OuterControl(P_control_typeAB, Q_control_typeAB)
        @test outer_control_typeAB isa PowerSystems.DynamicComponent
        test_accessors(outer_control_typeAB)
    end
end


text:
using LibCloud
using Base.Test

include("test_compute.jl")
include("test_storage.jl")
include("test_dns.jl")
include("test_containers.jl")
include("test_backup.jl")

using TestCompute
using TestStorage
using TestDNS
using TestContainers
using TestBackup

println("\nCOMPUTE\n===================")
TestCompute.dummy()
println("\nSTORAGE\n===================")
TestStorage.dummy()
println("\nDNS\n===================")
TestDNS.dummy()
println("\nContainers\n===================")
TestContainers.dummy()
println("\nBackup\n===================")
TestBackup.dummy()


text:
#using DataStructures, DataFrames, StatsFuns, GLM , JuMP,NLopt, StatStack, NLsolve, DBAPI, JDBC
#ENV["LD_LIBRARY_PATH"] = "/usr/lib/jvm/java-1.8.0-openjdk-1.8.0.71-2.b15.el7_2.x86_64/jre/lib/amd64/:/usr/lib/jvm/java-1.8.0-openjdk-1.8.0.71-2.b15.el7_2.x86_64/jre/lib/amd64/server/"
#JavaCall.addClassPath("/home/rmadmin/jdbc/ora/ojdbc7.jar")
#JavaCall.init(["-Duser.timezone=UTC","-Djava.class.path=/home/rmadmin/jdbc/ora/ojdbc7.jar"]);  #JDBC.init()

"""
ENV["LD_LIBRARY_PATH"] = "/usr/lib/jvm/java-1.7.0-openjdk-1.7.0.101.x86_64/jre/lib/amd64/:/usr/lib/jvm/java-1.7.0-openjdk-1.7.0.101.x86_64/jre/lib/amd64/server/"
JavaCall.addClassPath("/mapr/mapr04p/analytics0001/analytic_users/jpkg/v0.4/ojdbc7.jar")
JavaCall.init(["-Duser.timezone=UTC","-Djava.class.path=/mapr/mapr04p/analytics0001/analytic_users/jpkg/v0.4/ojdbc7.jar"]);

/opt/jdk1.7.0_79/jre/lib/amd64/server/libjvm.so

scp -r rmadmin@10.63.36.7:g/StatStack ./
0
"""

function saveDBora(sdf::DataFrame)
    #conn = DriverManager.getConnection("jdbc:oracle:thin:PRSPC/prspc@//ex02-scan.ch3.prod.i.com:1521/SV01DWHD")
    conn = DriverManager.getConnection("jdbc:oracle:thin:WH_MEDIA_SYND_P1/Passw0rd@//ex04-scan1.ch3.prod.i.com:1521/SV01DWHP")
    stmt = createStatement(conn)
    #rs = executeQuery(stmt, "select count(*) cnt from PRSPC.WH_PARAMS")

    rout=Int64[]
    rs = executeQuery(stmt, "select max(SCORE_MODEL_RUN_ID)+1 as cnt from WH_MEDIA_SYND_P1.LIFT_SCORE")
    for r in rs  push!(rout,getInt(r, 1))  end
    rid=rout[1]

    sql="insert into WH_MEDIA_SYND_P1.LIFT_SCORE_MODEL_RUN (SCORE_MODEL_RUN_ID,SCORE_MODEL_NAME,SCORE_MODEL_VER, RUN_DATE,RUN_STATUS, "
    sql = sql*"VIEW_DIM_TYPE, OPT_TYPE,OPT_MESR_TYPE, RELS_CD, BATCH_ID, REGISTRATION_REQUEST_ID, REGISTRATION_ID, "
    sql = sql*"TM_AGG_PERD,START_WEEK,END_WEEK) VALUES ( "
    sql = sql*string(rid)*", "
    sql = sql*"NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,NULL,188765,188765,99,9999,9999) "
    executeQuery(stmt, sql)

    pref = "insert into WH_MEDIA_SYND_P1.LIFT_SCORE "
    pref=pref*"    ( SCORE_MODEL_RUN_ID,MODEL_DESC,MODEL,DEPNDNT_VAR,CNT_HH_EXPSD,UNADJ_AVG_EXPSD_HH_PRE,UNADJ_AVG_CNTRL_HH_PRE,UNADJ_AVG_EXPSD_HH_PST,"
    pref=pref*"      UNADJ_AVG_CNTRL_HH_PST,UNADJ_DOD_EFFCT,UNADJ_DIFF_EFFCT,ADJ_MEAN_EXPSD_GRP,ADJ_MEAN_CNTRL_GRP,ADJ_DOD_EFFCT,TWOTAIL_PVAL,ONETAIL_PVAL, "
    pref=pref*"      ABS_DIFF,DOL_DIFF,ONETAIL_80_PCT_INTRVL_UB,ONETAIL_80_PCT_INTRVL_LB,ONETAIL_90_PCT_INTRVL_UB,ONETAIL_90_PCT_INTRVL_LB, "
    pref=pref*"  TWOTAIL_80_PCT_INTRVL_UB,TWOTAIL_80_PCT_INTRVL_LB,TWOTAIL_90_PCT_INTRVL_UB,TWOTAIL_90_PCT_INTRVL_LB,CNT_IMPRESNS,TWOTAIL_PVAL_TO_CAMPAIGN, "
    pref=pref*"      ONETAIL_PVAL_TO_CAMPAIGN,CNT_MODEL_HH,BATCH_ID "
    pref=pref*"    ) "
    pref=pref*"values "

    function isnax(v::Any) isna(v) ? true: isnan(v) ? true : false end

    for i in 1:length(sdf[1])
        println("Number : ",i)
        r=AbstractString[]
        push!(r,"'"*sdf[i,:MODEL_DESC]*"', " )
        isnax(sdf[1,:Model]) ?  push!(r,"NULL, ") : push!(r,"'"*sdf[1,:Model]*"', ")
        #isnax(sdf[i,:dependent_variable]) ? push!(r,"NULL, ") : push!(r,"'"*sdf[i,:dependent_variable]*"', " )
        push!(r,"'"*sdf[i,:dependent_variable]*"', " )
        isnax(sdf[i,:CNT_EXPSD_HH]) ? push!(r,"NULL, ") : push!(r, string(sdf[i,:CNT_EXPSD_HH])*", "   )
        isnax(sdf[i,:UNADJ_AVG_EXPSD_HH_PRE]) ? push!(r,"NULL, ") : push!(r,string(sdf[i,:UNADJ_AVG_EXPSD_HH_PRE])*", ")
        isnax(sdf[i,:UNADJ_AVG_CNTRL_HH_PRE]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:UNADJ_AVG_CNTRL_HH_PRE])*", ")
        isnax(sdf[i,:UNADJ_AVG_EXPSD_HH_PST]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:UNADJ_AVG_EXPSD_HH_PST])*", ")
        isnax(sdf[i,:UNADJ_AVG_CNTRL_HH_PST]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:UNADJ_AVG_CNTRL_HH_PST])*", ")
        isnax(sdf[i,:UNADJ_DOD_EFFCT]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:UNADJ_DOD_EFFCT])*", ")
        isnax(sdf[i,:UNADJ_DIFF_EFFCT]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:UNADJ_DIFF_EFFCT])*", ")
        isnax(sdf[i,:ADJ_MEAN_EXPSD_GRP]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ADJ_MEAN_EXPSD_GRP])*", ")
        isnax(sdf[i,:ADJ_MEAN_CNTRL_GRP]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ADJ_MEAN_CNTRL_GRP])*", ")
        isnax(sdf[i,:ADJ_DOD_EFFCT]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ADJ_DOD_EFFCT])*", ")
        isnax(sdf[i,:TWOTAIL_PVAL]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:TWOTAIL_PVAL])*", ")
        isnax(sdf[i,:ONETAIL_PVAL]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ONETAIL_PVAL])*", ")
        isnax(sdf[i,:ABS_DIFF]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ABS_DIFF])*", ")
        isnax(sdf[i,:DOL_DIFF]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:DOL_DIFF])*", ")
        isnax(sdf[i,:ONETAIL_80_PCT_INTRVL_UB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ONETAIL_80_PCT_INTRVL_UB])*", ")
        isnax(sdf[i,:ONETAIL_80_PCT_INTRVL_LB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ONETAIL_80_PCT_INTRVL_LB])*", ")
        isnax(sdf[i,:ONETAIL_90_PCT_INTRVL_UB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ONETAIL_90_PCT_INTRVL_UB])*", ")
        isnax(sdf[i,:ONETAIL_90_PCT_INTRVL_LB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ONETAIL_90_PCT_INTRVL_LB])*", ")
        isnax(sdf[i,:TWOTAIL_80_PCT_INTRVL_UB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:TWOTAIL_80_PCT_INTRVL_UB])*", ")
        isnax(sdf[i,:TWOTAIL_80_PCT_INTRVL_LB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:TWOTAIL_80_PCT_INTRVL_LB])*", ")
        isnax(sdf[i,:TWOTAIL_90_PCT_INTRVL_UB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:TWOTAIL_90_PCT_INTRVL_UB])*", ")
        isnax(sdf[i,:TWOTAIL_90_PCT_INTRVL_LB]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:TWOTAIL_90_PCT_INTRVL_LB])*", ")
        isnax(sdf[i,:CNT_IMPRESSIONS]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:CNT_IMPRESSIONS])*", ")
        isnax(sdf[i,:TWOTAIL_PVAL_to_Campaign]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:TWOTAIL_PVAL_to_Campaign])*", ")
        isnax(sdf[i,:ONETAIL_PVAL_to_Campaign]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:ONETAIL_PVAL_to_Campaign])*", ")
        isnax(sdf[i,:CNT_Model_HH]) ? push!(r,"NULL, "): push!(r,string(sdf[i,:CNT_Model_HH])*", ")

        sql=pref*" ("*string(rid)*", "
        for x in 1:length(r)  sql=sql*r[x] end
        sql=sql*"0.0)"
        #println(rout)
        executeQuery(stmt, sql)
    end

      """
        sql=pref*" ("*string(rid)*", "
        sql=sql*"'"*sdf[i,:MODEL_DESC]*"'"*", "
        sql=sql*"'"*sdf[i,:Model]*"'"*", "
        sql=sql*"'"*sdf[i,:DEPNDNT_VAR]*"'"*", "
        sql=sql+String(sdf[i,:CNT_HH_EXPSD])
        sql=sql+String(sdf[i,:UNADJ_AVG_EXPSD_HH_PRE])
        sql=sql+String(sdf[i,:UNADJ_AVG_CNTRL_HH_PRE])
        sql=sql+String(sdf[i,:UNADJ_AVG_EXPSD_HH_PST])
        sql=sql+String(sdf[i,:UNADJ_AVG_CNTRL_HH_PST])
        sql=sql+String(sdf[i,:UNADJ_DOD_EFFCT])
        sql=sql+String(sdf[i,:UNADJ_DIFF_EFFCT])
        sql=sql+String(sdf[i,:ADJ_MEAN_EXPSD_GRP])
        sql=sql+String(sdf[i,:ADJ_MEAN_CNTRL_GRP])
        sql=sql+String(sdf[i,:ADJ_DOD_EFFCT])
        sql=sql+String(sdf[i,:TWOTAIL_PVAL])
        sql=sql+String(sdf[i,:ONETAIL_PVAL])
        sql=sql+String(sdf[i,:ABS_DIFF])
        sql=sql+String(sdf[i,:DOL_DIFF])
        sql=sql+String(sdf[i,:ONETAIL_80_PCT_INTRVL_UB])
        sql=sql+String(sdf[i,:ONETAIL_80_PCT_INTRVL_LB])
        sql=sql+String(sdf[i,:ONETAIL_90_PCT_INTRVL_UB])
        sql=sql+String(sdf[i,:ONETAIL_90_PCT_INTRVL_LB])
        sql=sql+String(sdf[i,:TWOTAIL_80_PCT_INTRVL_UB])
        sql=sql+String(sdf[i,:TWOTAIL_80_PCT_INTRVL_LB])
        sql=sql+String(sdf[i,:TWOTAIL_90_PCT_INTRVL_UB])
        sql=sql+String(sdf[i,:TWOTAIL_90_PCT_INTRVL_LB])
        sql=sql+String(sdf[i,:CNT_IMPRESNS])
        sql=sql+String(sdf[i,:TWOTAIL_PVAL_TO_CAMPAIGN])
        sql=sql+String(sdf[i,:ONETAIL_PVAL_TO_CAMPAIGN])
        sql=sql+String(sdf[i,:CNT_MODEL_HH])
        sql=sql+"0.0"   #String(sdf[1,:BATCH_ID])  #BATCH_ID
    """

    #println(sql)
    #ji = JDBCRowIterator(rs); getInt(ji.rs, 1)
    #ji = JDBCRowIterator(rs); typeof(ji.get_methods[1](ji.rs, 1))
    #ji=JDBCRowIterator(rs); getInt(ji.get_methods[1](ji.rs, 1), 1)

    #rid=getInt(JDBCRowIterator(rs)[1], 1)
    #println(getInt(r, 1))
    #for r in rs  println(getInt(r, 1),  getString(r,"cnt"))  end

    #collect(or r in rs  getInt(r, 1)  end)[1]

    #execute!(csr, "insert into pi_table (pi_value) values (3.14);")
    #execute!(csr, "select * from my_table;")

    #csr = cursor(conn)
    #execute!(csr, "select count(*) cnt from PRSPC.WH_PARAMS")
    #rs = rows(csr)
    #for row in rs
    #    println(row[:cnt])
    #end
    close(conn)
end
#saveDBora(sdf)

max_stars_repo_name: werbac/jlib

max_issues_repo_name: werbac/jlib

max_forks_repo_name: werbac/jlib


text:
function add_entry!(names, values, entry; default)
    i = findfirst(==(entry), names)
    if isnothing(i)
        push!(names, entry)
        push!(values, default)
        i = lastindex(names)
    end
    return values[i]
end

struct LegendSection
    title::String
    names::Vector{String}
    plots::Vector{Vector{AbstractPlot}}
end
LegendSection(title::String="") = LegendSection(title, String[], Vector{AbstractPlot}[])

# Add an empty trace list with name `entry` to the legend section
function add_entry!(legendsection::LegendSection, entry::String)
    names, plots = legendsection.names, legendsection.plots
    return add_entry!(names, plots, entry; default=AbstractPlot[])
end

struct Legend
    names::Vector{String}
    sections::Vector{LegendSection}
end
Legend() = Legend(String[], LegendSection[])

# Add an empty section with name `entry` and title `title` to the legend
function add_entry!(legend::Legend, entry::String; title::String="")
    names, sections = legend.names, legend.sections
    return add_entry!(names, sections, entry; default=LegendSection(title))
end

function create_legend(scene, legend::Legend)
    sections = legend.sections
    MakieLayout.LLegend(
        scene,
        getproperty.(sections, :plots),
        getproperty.(sections, :names),
        getproperty.(sections, :title)
    )
end


text:
using Plots

# enable plotting by default
if !@isdefined do_visu; do_visu = true end

@views function diffusion_1D_2procs(; do_visu=true)
    # Physics
    Hl  = 10.0   # left  H
    Hr  = 1.0    # right H
    λ   = 1.0    # diffusion coeff
    nt  = 200    # number of time steps
    # Numerics
    nx  = 32     # number of local grid points
    dx  = 1.0    # cell size
    # Derived numerics
    dt  = dx^2/λ/2.1
    # Initial condition
    HL  = Hl*ones(nx)
    HR  = Hr*ones(nx)
    H   = [HL[1:end-1]; HR[2:end]]
    Hg  = copy(H)
    # Time loop
    for it = 1:nt
        # Compute physics locally
        HL[2:end-1] .= HL[2:end-1] + dt*λ*diff(diff(HL)/dx)/dx
        HR[2:end-1] .= HR[2:end-1] + dt*λ*diff(diff(HR)/dx)/dx
        # Update boundaries (MPI)
        HL[end] = HR[2]
        HR[1]   = HL[end-1]
        # Global picture
        H .= [HL[1:end-1]; HR[2:end]]
        # Compute physics globally (check)
        Hg[2:end-1] .= Hg[2:end-1] + dt*λ*diff(diff(Hg)/dx)/dx
        # Visualise
        if do_visu
            fontsize = 12
            plot(Hg, legend=false, linewidth=0, markershape=:circle, markersize=5, yaxis=font(fontsize, "Courier"), xaxis=font(fontsize, "Courier"), titlefontsize=fontsize, titlefont="Courier")
            display(plot!(H, legend=false, linewidth=3, framestyle=:box, xlabel="Lx", ylabel="H", title="diffusion (it=$(it))"))
        end
    end
    return
end

diffusion_1D_2procs(; do_visu=do_visu)


text:
using IterativeSolvers
using Test
using Random
using LinearAlgebra

@testset ("BiCGStab(l)") begin

Random.seed!(1234321)
n = 20

@testset "Matrix{$T}" for T in (Float32, Float64, ComplexF32, ComplexF64)
    A = rand(T, n, n) + 15I
    x = ones(T, n)
    b = A * x

    tol = √eps(real(T))

    @testset "BiCGStab($l)" for l = (2, 4)
        # Solve without preconditioner
        x1, his1 = bicgstabl(A, b, l, max_mv_products = 100, log = true, tol = tol)
        @test isa(his1, ConvergenceHistory)
        @test norm(A * x1 - b) / norm(b) ≤ tol

        # With an initial guess
        x_guess = rand(T, n)
        x2, his2 = bicgstabl!(x_guess, A, b, l, max_mv_products = 100, log = true, tol = tol)
        @test isa(his2, ConvergenceHistory)
        @test x2 == x_guess
        @test norm(A * x2 - b) / norm(b) ≤ tol

        # Do an exact LU decomp of a nearby matrix
        F = lu(A + rand(T, n, n))
        x3, his3 = bicgstabl(A, b, Pl = F, l, max_mv_products = 100, log = true, tol = tol)
        @test norm(A * x3 - b) / norm(b) ≤ tol
    end
end
end


text:
"""
Functionality for dealing with HTTP cookies.
"""
module Cookies

using HttpServer, HttpCommon, Genie, Encryption

"""
    get(res::Response, key::Union{String,Symbol}) :: Nullable{String}

Retrieves a value stored on the cookie as `key` from the `Respose` object.
"""
function get(res::Response, key::Union{String,Symbol}) :: Nullable{String}
  if haskey(res.cookies, string(key))
    return Nullable{String}(res.cookies[string(key)].value |> Encryption.decrypt)
  end

  Nullable{String}()
end

"""
    get(req::Request, key::Union{String,Symbol}) :: Nullable{String}

Retrieves a value stored on the cookie as `key` from the `Request` object.
"""
function get(req::Request, key::Union{String,Symbol}) :: Nullable{String}
  if haskey(req.headers, "Cookie")
    cookies = to_dict(req)
    if haskey(cookies, string(key))
      return Nullable{String}(cookies[string(key)] |> Encryption.decrypt)
    end
  end

  Nullable{String}()
end

"""
    set!(res::Response, key::Union{String,Symbol}, value::Any, attributes::Dict) :: Dict{String,HttpCommon.Cookie}
    set!(res::Response, key::Union{AbstractString,Symbol}, value::Any) :: Dict{String,HttpCommon.Cookie}

Sets `value` under the `key` label on the `Cookie`.
"""
function set!(res::Response, key::Union{String,Symbol}, value::Any, attributes::Dict) :: Dict{String,HttpCommon.Cookie}
  setcookie!(res, string(key), string(value) |> Encryption.encrypt, attributes)

  res.cookies
end
function set!(res::Response, key::Union{AbstractString,Symbol}, value::Any) :: Dict{String,HttpCommon.Cookie}
  set!(res, key, value, Dict())
end

"""
    to_dict(req::Request) :: Dict{String,String}

Extracts the `Cookie` data from the `Request` and converts it into a dict.
"""
function to_dict(req::Request) :: Dict{String,String}
  d = Dict{String,String}()
  for cookie in split(req.headers["Cookie"], ";")
    cookie_parts = split(cookie, "=")
    d[strip(cookie_parts[1])] = cookie_parts[2]
  end

  d
end

end


text:
export angular_distance, azimuth, angular_step, surface_distance, forward
export inverse, get_dist, get_azi, get_baz, get_loc

"""
angular_distance, azimuth, angular_step, surface_distance, forward and inverse
functions were written by Andy Nowacki in his (unregistered) Geodesics.jl
package https://github.com/anowacki/Geodesics.jl
"""

"Earth ellipsoid semi-axes in WGS84"
const EARTH_R_MAJOR_WGS84 = 6378137.0000
const EARTH_R_MINOR_WGS84 = 6356752.3142
"Flattening of the Earth in WGS84"
const F_WGS84 = (EARTH_R_MAJOR_WGS84 - EARTH_R_MINOR_WGS84)/EARTH_R_MAJOR_WGS84

"""
    angular_distance(lon0, lat0, lon1, lat1, degrees=true; f=F_WGS84) -> Δ
Return the angular distance between points (`lon0`,`lat0`) and (`lon1`,`lat1`)
on a flattened sphere.  The default flattening is $F_WGS84.  By default,
input and output are in degrees, but specify `degrees` as `false` to use radians.
Use `f=0` to perform computations on a sphere.
"""
function angular_distance(lon0, lat0, lon1, lat1, degrees::Bool=true; f=F_WGS84)
    if degrees
        lon0, lat0, lon1, lat1 = deg2rad(lon0), deg2rad(lat0), deg2rad(lon1), deg2rad(lat1)
    end
    gcarc, az, baz = inverse(lon0, lat0, lon1, lat1, 1.0, f)
    degrees ? rad2deg(gcarc) : gcarc
end

"""
    azimuth(lon0, lat0, lon1, lat1, degrees=true, f=F_WGS84)
Return the azimuth from point (`lon0`,`lat0`) to point (`lon1`,`lat1`) on a
flattened sphere.  The default flattening is $F_WGS84.  By default,
input and output are in degrees, but specify `degrees` as `false` to use radians.
Use `f=0` to perform computations on a sphere.
"""
function azimuth(lon0::T1, lat0::T2, lon1::T3, lat1::T4, degrees::Bool=true; f=F_WGS84) where {T1,T2,T3,T4}
    # Promote to Float64 at least to avoid error when one of the points is on the pole
    T = promote_type(Float64, T1, T2, T3, T4)
    lon0, lat0, lon1, lat1 = T(lon0), T(lat0), T(lon1), T(lat1)
    if degrees
        lon0, lat0, lon1, lat1 = deg2rad(lon0), deg2rad(lat0), deg2rad(lon1), deg2rad(lat1)
    end
    gcarc, az, baz = inverse(lon0, lat0, lon1, lat1, 1.0, f)
    degrees ? rad2deg(az) : az
end

"""
    step(lon, lat, azimuth, distance, degrees=false; f=F_WGS84) -> lon′, lat′, backazimuth
Return the longitude `lon′`, latitude `lat′` and backazimuth `baz` reached by
travelling an angular `distance` along `azimuth` from the starting point at
(`lon`,`lat`).  The default flattening is $F_WGS84.  By default,
input and output are in degrees, but specify `degrees` as `false` to use radians.
Use `f=0` to perform computations on a sphere.
"""
function angular_step(lon, lat, azimuth, distance, degrees::Bool=true; f=F_WGS84)
    if degrees
        lon, lat, azimuth, distance = deg2rad(lon), deg2rad(lat), deg2rad(azimuth), deg2rad(distance)
    end
    lon′, lat′, baz = forward(lon, lat, azimuth, distance, 1.0, f)
    if degrees
        rad2deg(lon′), rad2deg(lat′), rad2deg(baz)
    else
        lon′, lat′, baz
    end
end

"""
    surface_distance(lon0, lat0, lon1, lat1, a, degrees::Bool=false; f=F_WGS84)
Return the physical distance between points (`lon0`,`lat0`) and (`lon1`,`lat1`) on
the flattened sphere with flattening `f` and semimajor radius `a`.  Distance is given
in the same units as `a`.  By default, input angles are in degrees, but specify `degrees`
as `false` to use radians.  Use `f=0` to perform computations on a sphere.
"""
function surface_distance(lon0, lat0, lon1, lat1, a, degrees::Bool=true; f=F_WGS84)
    if degrees
        lon0, lat0, lon1, lat1 = deg2rad(lon0), deg2rad(lat0), deg2rad(lon1), deg2rad(lat1)
    end
    distance, az, baz = inverse(lon0, lat0, lon1, lat1, a, f)
    distance
end

"""
    forward(lon, lat, azimuth, distance, a, f) -> lon′, lat′, backazimuth
Return the longitude `lon′` and latitude `lat′` and `backazimuth` of a projected
point, reached by travelling along an `azimuth` for `distance` from an original
point at (`lon`, `lat`).  Specify the spheroid with the semimajor radius `a` and
flattening `f`.
Coordinates and azimuth are in radians.
Calculations use Vincenty's forward formula [1].
#### References
1. Vincenty, T. (1975). "Direct and Inverse Solutions of Geodesics on the Ellipsoid
   with application of nested equations" (PDF). Survey Review. XXIII (176): 88–93.
   doi:10.1179/sre.1975.23.176.88
"""
function forward(lon, lat, azimuth, distance, a, f)::Tuple{Float64,Float64,Float64}
    abs(lat <= π/2) || throw(ArgumentError("Latitude ($lat) must be in range [-π/2, π/2]"))
    a > 0 || throw(ArgumentError("Semimajor axis ($a) must be positive"))
    abs(f) < 1 || throw(ArgumentError("Magnitude of flattening ($f) must be less than 1"))
    # Calculations are done with Float64s internally as the tolerances are hard-wired
    lambda1, phi1, alpha12, s = Float64(lon), Float64(lat), Float64(azimuth), Float64(distance)
    a, f = Float64(a), Float64(f)
    alpha12 = mod(alpha12, 2π)
    b = a*(1 - f)

    TanU1 = (1 - f)*tan(phi1)
    U1 = atan(TanU1)
    sigma1 = atan( TanU1, cos(alpha12) )
    Sinalpha = cos(U1)*sin(alpha12)
    cosalpha_sq = 1.0 - Sinalpha*Sinalpha

    u2 = cosalpha_sq*(a*a - b*b )/(b*b)
    A = 1.0 + (u2/16384)*(4096 + u2*(-768 + u2*(320 - 175*u2)))
    B = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2)))

    # Starting with the approximation
    sigma = (s/(b*A))

    # Not moving anywhere. We can return the location that was passed in.
    if sigma == 0
        return phi1, lambda1, alpha12
    end

    last_sigma = 2*sigma + 2 # something impossible

    # Iterate the following three equations
    # until there is no significant change in sigma
    # two_sigma_m , delta_sigma
    while abs((last_sigma - sigma)/sigma) > 1.0e-9
        global two_sigma_m = 2*sigma1 + sigma
        delta_sigma = B*sin(sigma)*(cos(two_sigma_m) + (B/4)*(cos(sigma)*(-1 + 2*cos(two_sigma_m)^2 - (B/6)*cos(two_sigma_m)*(-3 + 4*sin(sigma)^2)*(-3 + 4*cos(two_sigma_m)^2 ))))
        last_sigma = sigma
        sigma = (s/(b*A)) + delta_sigma
    end

    phi2 = atan((sin(U1)*cos(sigma) + cos(U1)*sin(sigma)*cos(alpha12)),
        ((1-f)*sqrt(Sinalpha^2 + (sin(U1)*sin(sigma) - cos(U1)*cos(sigma)*cos(alpha12))^2)))

    lambda = atan((sin(sigma)*sin(alpha12)),
        (cos(U1)*cos(sigma) - sin(U1)*sin(sigma)*cos(alpha12)))

    C = (f/16)*cosalpha_sq*(4 + f*(4 - 3*cosalpha_sq))

    omega = lambda - (1-C)*f*Sinalpha*(sigma + C*sin(sigma)*(
        cos(two_sigma_m) + C*cos(sigma)*(-1 + 2*cos(two_sigma_m)^2)))

    lambda2 = lambda1 + omega

    alpha21 = atan(Sinalpha, (-sin(U1)*sin(sigma) + cos(U1)*cos(sigma)*cos(alpha12)))
    alpha21 = mod(alpha21 + π, 2π)

    return lambda2, phi2, alpha21
end

"""
    inverse(lon1, lat1, lon2, lat2, a, f) -> distance, azimuth, backazimuth
Return the `distance`, `azimuth` and `backazimuth` between two points with longitudes
`lon1` and `lon2`, and latitudes `lat1` and `lat2`.  Specify the spheroid with the
semimajor radius `a` and flattening `f`.
Coordinates and angles are in radians, whilst `distance` is in the same units as `a`.
Calculations use Vincenty's inverse formula [1].
#### References
1. Vincenty, T. (1975). "Direct and Inverse Solutions of Geodesics on the Ellipsoid
   with application of nested equations" (PDF). Survey Review. XXIII (176): 88–93.
   doi:10.1179/sre.1975.23.176.88
"""
function inverse(lon1, lat1, lon2, lat2, a, f)::Tuple{Float64,Float64,Float64}
    for lat in (lat1, lat2)
        abs(lat <= π/2) || throw(ArgumentError("Latitude ($lat) must be in range [-π/2, π/2]"))
    end
    a > 0 || throw(ArgumentError("Semimajor axis ($a) must be positive"))
    abs(f) < 1 || throw(ArgumentError("Magnitude of flattening ($f) must be less than 1"))
    lambda1, phi1, lambda2, phi2 = Float64(lon1), Float64(lat1), Float64(lon2), Float64(lat2)
    a, f = Float64(a), Float64(f)
    tol = 1.0e-8
    if (abs(phi2 - phi1) < tol) && (abs(lambda2 - lambda1) < tol)
        return 0.0, 0.0, 0.0
    end

    b = a*(1 - f)

    TanU1 = (1 - f)*tan(phi1)
    TanU2 = (1 - f)*tan(phi2)

    U1 = atan(TanU1)
    U2 = atan(TanU2)

    lambda = lambda2 - lambda1
    last_lambda = -4000000.0 # an impossibe value
    omega = lambda

    # Iterate the following equations until there is no significant change in lambda
    alpha, sigma, Sin_sigma, Cos2sigma_m, Cos_sigma, sqr_sin_sigma =
        -999999., -999999., -999999., -999999., -999999., -999999.
    while ((last_lambda < -3000000.0) || (lambda != 0)) &&
            (abs((last_lambda - lambda)/lambda) > 1.0e-9)
        sqr_sin_sigma = (cos(U2)*sin(lambda))^2 +
                         ((cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda)))^2
        Sin_sigma = sqrt(sqr_sin_sigma)
        Cos_sigma = sin(U1)*sin(U2) + cos(U1)*cos(U2)*cos(lambda)
        sigma = atan(Sin_sigma, Cos_sigma)

        Sin_alpha = cos(U1)*cos(U2)*sin(lambda)/sin(sigma)

        if Sin_alpha >= 1
            Sin_alpha = 1.0
        elseif Sin_alpha <= -1
            Sin_alpha = -1.0
        end

        alpha = asin(Sin_alpha)
        Cos2sigma_m = cos(sigma) - 2*sin(U1)*sin(U2)/cos(alpha)^2
        C = (f/16)*cos(alpha)^2*(4 + f*(4 - 3*cos(alpha)^2))
        last_lambda = lambda
        lambda = omega + (1 - C)*f*sin(alpha)*(sigma +
            C*sin(sigma)*(Cos2sigma_m + C*cos(sigma)*(-1 + 2*Cos2sigma_m^2)))
    end

    u2 = cos(alpha)^2*(a*a - b*b)/(b*b)
    A = 1 + (u2/16384)*(4096 + u2*(-768 + u2*(320 - 175*u2)))
    B = (u2/1024)*(256 + u2*(-128 + u2*(74 - 47*u2)))
    delta_sigma = B*Sin_sigma*(Cos2sigma_m + (B/4)*(
        Cos_sigma*(-1 + 2*Cos2sigma_m^2) -
        (B/6)*Cos2sigma_m*(-3 + 4*sqr_sin_sigma)*(-3 + 4*Cos2sigma_m^2)))
    s = b*A*(sigma - delta_sigma)

    alpha12 = atan((cos(U2)*sin(lambda)), ( cos(U1)*sin(U2) - sin(U1)*cos(U2)*cos(lambda)))
    alpha21 = atan((cos(U1)*sin(lambda)), (-sin(U1)*cos(U2) + cos(U1)*sin(U2)*cos(lambda)))

    alpha12 = mod(alpha12, 2π)
    alpha21 = mod(alpha21 + π, 2π)

    return s, alpha12, alpha21
end

"""
  get_dist(loc1,loc2)

Get distance bewteen source `loc1` and receiver `loc2` in km.
If either location is empty, return zero.
# Arguments
- `loc1::GeoLoc`: Location of source.
- `loc2::GeoLoc`: Location of receiver.
"""
function get_dist(loc1::GeoLoc, loc2::GeoLoc)
    if isempty(loc1) | isempty(loc2)
        dist =  zero(Float64)
    else
        dist =  surface_distance(loc1.lon,loc1.lat,loc2.lon,loc2.lat,
                               EARTH_R_MAJOR_WGS84/1000)
    end
    return dist
end

"""
  get_dist(FFT1,FFT2)

Get distance bewteen source `FFT1` and receiver `FFT2` in km.
If either location is empty, return zero.
# Arguments
- `FFT1::FFTData`: FFTData with source location.
- `FFT2::FFTData`: FFTData with receiver location.
"""
function get_dist(FFT1::FFTData, FFT2::FFTData)
    return get_dist(FFT1.loc,FFT2.loc)
end

"""
  get_azi(loc1,loc2)

Get azimuth from source `loc1` to receiver `loc2` in degrees.
If either location is empty, return zero.
# Arguments
- `loc1::GeoLoc`: Location of source.
- `loc2::GeoLoc`: Location of receiver.
"""
function get_azi(loc1::GeoLoc, loc2::GeoLoc)
    if isempty(loc1) | isempty(loc2)
        azi =  zero(Float64)
    else
        azi = azimuth(loc1.lon,loc1.lat,loc2.lon,loc2.lat)
    end
    return azi
end

"""
  get_azi(FFT1,FFT2)

Get azimuth from source `FFT1` to receiver `FFT2` in degrees.
If either location is empty, return zero.
# Arguments
- `FFT1::FFTData`: FFTData with source location.
- `FFT2::FFTData`: FFTData with receiver location.
"""
function get_azi(FFT1::FFTData, FFT2::FFTData)
    return get_azi(FFT1.loc,FFT2.loc)
end

"""
  get_baz(loc1,loc2)

Get back azimuth from source `loc1` to receiver `loc2` in degrees.
If either location is empty, return zero.
# Arguments
- `loc1::GeoLoc`: Location of source.
- `loc2::GeoLoc`: Location of receiver.
"""
function get_baz(loc1::GeoLoc, loc2::GeoLoc)
    if isempty(loc1) | isempty(loc2)
        baz =  zero(Float64)
    else
        baz = inverse(deg2rad.((loc1.lon,loc1.lat,loc2.lon,loc2.lat))...,
                          EARTH_R_MAJOR_WGS84, F_WGS84)[3]
        baz = rad2deg(baz)
    end
    return baz
end

"""
  get_baz(FFT1,FFT2)

Get back azimuth from source `FFT1` to receiver `FFT2` in degrees.
If either location is empty, return zero.
# Arguments
- `FFT1::FFTData`: FFTData with source location.
- `FFT2::FFTData`: FFTData with receiver location.
"""
function get_baz(FFT1::FFTData, FFT2::FFTData)
    return get_baz(FFT1.loc,FFT2.loc)
end

"""
  get_loc(loc1,azi, dist)

Get location of receiver given location of source and dist + azi to receiver.
If location is empty, return empty location.
# Arguments
- `loc1::GeoLoc`: Location of source.
- `azi::Float64`: Azimuth from source `loc1` to receiver.
- `dist::Float64`: Distance from source `loc1` to receiver in km.
"""
function get_loc(loc1::GeoLoc, azi::Float64, dist::Float64)
    if isempty(loc1)
        loc2 = GeoLoc()
    else
        lon2,lat2,baz2 = forward(deg2rad.((loc1.lon, loc1.lat,azi))...,
                         dist*1000, EARTH_R_MAJOR_WGS84, F_WGS84)
        lon2, lat2, baz2 = rad2deg.((lon2, lat2, baz2))
        loc2 = GeoLoc("",lat2,lon2, zero(Float64), zero(Float64), zero(Float64),
                      loc1.inc)
    end
    return loc2
end

"""
  get_loc(C)

Get location of receiver and source from CorrData `C`.
If either location is empty, return empty locations.
# Arguments
- `C::CorrData`: CorrData with source location and distance + azimuth to receiver.
"""
function get_loc(C::CorrData)
    if isempty(C.loc)
        loc2 = GeoLoc()
    else
        loc2 = get_loc(C.loc, C.azi, C.dist)
    end

    return C.loc, loc2
end


text:
function GBVector{T}(n) where {T}
    GBVector{T}(libgb.GrB_Vector_new(toGrB_Type(T),n))
end
GBVector{T}() where {T} = GBVector{T}(libgb.GxB_INDEX_MAX)

Base.unsafe_convert(::Type{libgb.GrB_Vector}, v::GBVector) = v.p

function Base.copy(v::GBVector{T}) where {T}
    return GBVector{T}(libgb.GrB_Vector_dup(v))
end

clear!(v::GBVector) = libgb.GrB_Vector_clear(v)

function Base.size(v::GBVector, dim = nothing)
    if dim === nothing || dim == 1
        return Int64(libgb.GrB_Vector_size(v))
    else
        return 1
    end
end

nnz(v::GBVector) = Int64(libgb.GrB_Vector_nvals(v))
Base.eltype(::Type{GBVector{T}}) where{T} = T

#Similar functions. This definition is kind of clunky, but there were overwrite errors if I did it in another way.
Base.similar(v::GBVector) = GBVector{eltype(v)}(size(v))
Base.similar(::GBVector{T}, dims::Integer) where {T} = GBVector{T}(dims)
Base.similar(v::GBVector, element_type::DataType) = GBVector{element_type}(size(v))
Base.similar(::GBVector, element_type::DataType, dims::Integer) = GBVector{element_type}(dims)

for T ∈ valid_vec
    if T ∈ GxB_vec
        prefix = :GxB
    else
        prefix = :GrB
    end
    func = Symbol(prefix, :_Vector_build_, suffix(T))
    @eval begin
        function build(v::GBVector{$T}, I::Vector, X::Vector{$T}; dup = BinaryOps.PLUS)
            nnz(v) == 0 || error("Cannot build vector with existing elements")
            length(X) == length(I) || DimensionMismatch("I and X must have the same length")
            libgb.$func(v, Vector{libgb.GrB_Index}(I), X, length(X), dup[$T])
        end
    end
    func = Symbol(prefix, :_Vector_setElement_, suffix(T))
    @eval begin
        function Base.setindex!(v::GBVector{$T}, x::$T, i::Integer)
            return libgb.$func(v, x, libgb.GrB_Index(i))
        end
    end
    func = Symbol(prefix, :_Vector_extractElement_, suffix(T))
    @eval begin
        function Base.getindex(v::GBVector{$T}, i::Integer)
            return libgb.$func(v, libgb.GrB_Index(i))
        end
    end
    func = Symbol(prefix, :_Vector_extractTuples_, suffix(T))
    @eval begin
        function SparseArrays.findnz(v::GBVector{$T})
            return libgb.$func(v)
        end
    end
end

#Determine the length of output vector for indexing operations.
function wlength(u, I)
    if I == ALL
        wlen = size(u)
    else
        wlen = length(I)
    end
    return wlen
end

#Determine the index vector depending on type and contents of I.

@kwargs function extract!(w::GBVector, u::GBVector, I)
    I, ni = idx(I)
    libgb.GrB_Vector_extract(w, mask, accum, u, I, ni, desc)
    return w
end

@kwargs function extract(u::GBVector, I)
    wlen = wlength(u, I)
    w = similar(u, wlen)
    return extract!(w, u, I; mask, accum, desc)
end

@kwargs Base.getindex(u::GBVector, ::Colon) = extract(u, ALL; mask, accum, desc)
@kwargs function Base.getindex(u::GBVector, i::Union{Vector, UnitRange, StepRange})
    return extract(u, i; mask, accum, desc)
end

@kwargs function subassign!(w::GBVector, u, I)
    I, ni = idx(I)
    u isa Vector && (u = GBVector(u))
    if u isa GBVector
        libgb.GxB_Vector_subassign(w, mask, accum, u, I, ni, desc)
    else
        libgb.scalarvecsubassign[eltype(u)](w, mask, accum, u, I, ni, desc)
    end
end

@kwargs function assign!(w::GBVector, u, I)
    I, ni = idx(I)
    u isa Vector && (u = GBVector(u))
    if u isa GBVector
        libgb.GrB_Vector_assign(w, mask, accum, u, I, ni, desc)
    else
        libgb.scalarvecassign[eltype(u)](w, mask, accum, u, I, ni, desc)
    end
end

#setindex! uses subassign rather than assign. This behavior may change in the future.
@kwargs Base.setindex!(u::GBVector, x, ::Colon) = subassign!(u, x, ALL; mask, accum, desc)
@kwargs function Base.setindex!(u::GBVector, x, I::Union{Vector, UnitRange, StepRange})
    subassign!(u, x, I; mask, accum, desc)
end

function GBVector(I::Vector, X::Vector{T}; dup = BinaryOps.PLUS) where {T}
    x = GBVector{T}(maximum(I))
    build(x, I, X, dup = dup)
    return x
end

#TEMPORARY: NEEDS IMPORT/EXPORT
function GBVector(X::Vector)
    return GBVector(collect(1:length(X)), X) #Collect... Ouch.
end

#TEMPORARY: NEEDS IMPORT/EXPORT
function GBVector(A::SparseVector)
    i, k = findnz(A)
    return GBVector(i, k)
end
function Base.deleteat!(v::GBVector, i)
    libgb.GrB_Vector_removeElement(v, libgb.GrB_Index(i))
    return v
end

function Base.resize!(v::GBVector, nrows_new)
    libgb.GrB_Vector_resize(v, nrows_new)
    return v
end

function diagonal(A::GBMatrix, k = 0, desc = C_NULL)
    return libgb.GxB_Vector_diag(A, k, desc)
end

function printtest(io::IO, M::GBVector)
    str = mktemp() do _, f
        cf = Libc.FILE(f)
        libgb.GxB_Vector_fprint(M, "Test", libgb.GxB_SHORT, cf)
        flush(f)
        seekstart(f)
        x = read(f, String)
        close(cf)
        x
    end
    print(str)
end

function Base.show(io::IO, ::MIME"text/plain", v::GBVector)
    printtest(io, v)
end

function Base.print(io::IO, v::GBVector)
    printtest(io, v)
end


text:
##### Numeric simplification

"""
simplify(x; rewriter=default_simplifier(),
            threaded=false,
            polynorm=true,
            thread_subtree_cutoff=100)

Simplify an expression (`x`) by applying `rewriter` until there are no changes.
`polynorm=true` applies `polynormalize` in the beginning of each fixpoint iteration.
"""
function simplify(x;
                  polynorm=false,
                  threaded=false,
                  thread_subtree_cutoff=100,
                  rewriter=nothing)
    f = if rewriter === nothing
        if threaded
            threaded_simplifier(thread_subtree_cutoff)
        elseif polynorm
            serial_polynormal_simplifier
        else
            serial_simplifier
        end
    else
        Fixpoint(rewriter)
    end

    PassThrough(f)(to_symbolic(x))
end

Base.@deprecate simplify(x, ctx; kwargs...)  simplify(x; rewriter=ctx, kwargs...)

"""
    substitute(expr, dict)

substitute any subexpression that matches a key in `dict` with
the corresponding value.
"""
function substitute(expr, dict; fold=true)
    rs = Prewalk(PassThrough(@rule ~x::(x->haskey(dict, x)) => dict[~x]))
    if fold
        rs(to_symbolic(expr)) |> SymbolicUtils.fold
    else
        rs(to_symbolic(expr))
    end
end


text:
## test_karhunen_loeve.jl : test Karhunen-Lo\`eve expansion for GRF generation

@testset "Karhunen-Lo\u00e8ve expansion " begin

## 1d Mat\'ern ##
cov = CovarianceFunction(1,Matern(0.3,1))
pts = 0:0.01:1
grf = GaussianRandomField(cov,KarhunenLoeve(500),pts)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Matern)
@test ndims(grf.cov) == 1
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{500}} where {C})
@test length(grf.pts[1]) == length(pts)
@test length(sample(grf)) == length(pts)

## test anisotropic
cov = CovarianceFunction(2,AnisotropicExponential(100*[1 0.2; 0.2 1]))
pts1 = range(0,stop = 1,length = 128)
pts2 = range(0,stop = 1,length = 128)
grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,pts2)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,AnisotropicExponential)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{1000}} where {C})
@test length(grf.pts) == 2
@test length(grf.pts[1]) == length(pts1)
@test length(grf.pts[2]) == length(pts2)
@test size(sample(grf),1) == length(pts1)
@test size(sample(grf),2) == length(pts2)

## variation of number of KL terms ##
@test_throws ArgumentError KarhunenLoeve(0)

cov = CovarianceFunction(2,Matern(0.3,1))
pts = 0:0.01:1
grf = GaussianRandomField(cov,KarhunenLoeve(1),pts,pts)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Matern)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{1}} where {C})
@test length(grf.pts[1]) == length(pts)
@test length(grf.pts[2]) == length(pts)
@test size(sample(grf),1) == length(pts)
@test size(sample(grf),2) == length(pts)

mat = Matern(0.1,1.0)
cov = CovarianceFunction(1,mat)
pts1 = 0:0.01:1
nterms = [1 2 5 10 20 50 100 200 500 1000]
err = Float64[]
for n in nterms
	grf = GaussianRandomField(cov,KarhunenLoeve(n),pts1)
	push!(err,rel_error(grf))
end
@test err[1]/err[end] > 1e7

## quadrature methods ##
mat = Matern(0.1,1.0)
cov = CovarianceFunction(1,mat)
pts1 = 0:0.025:1
grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=GaussLegendre())
@test rel_error(grf) < 1e-6
grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=EOLE())
@test rel_error(grf) < 1e-6
@suppress grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=Simpson())
@test rel_error(grf) < 1e-6
grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=Midpoint())
@test rel_error(grf) < 1e-6
grf = GaussianRandomField(cov,KarhunenLoeve(1000),pts1,quad=Trapezoidal())
@test rel_error(grf) < 1e-6

## test 3D Gaussian ##
cov = CovarianceFunction(3,Gaussian(1.0))
pts1 = 0:0.05:1
pts2 = 0:0.05:1
pts3 = 0:0.05:1
grf = GaussianRandomField(cov,KarhunenLoeve(100),pts1,pts2,pts3)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Gaussian)
@test ndims(grf.cov) == 3
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{100}} where {C})
@test length(grf.pts) == 3
@test length(grf.pts[1]) == length(pts1)
@test length(grf.pts[2]) == length(pts2)
@test length(grf.pts[3]) == length(pts3)
@test size(sample(grf),1) == length(pts1)
@test size(sample(grf),2) == length(pts2)
@test size(sample(grf),3) == length(pts3)

## variation of number of points ##
cov = CovarianceFunction(1,Matern(1,1))
pts = 0:0.1:0
@test_throws ArgumentError grf = GaussianRandomField(cov,KarhunenLoeve(300),pts)
pts = 0:0.1:0.1
grf = GaussianRandomField(cov,KarhunenLoeve(300),pts)
pts = 1:-0.1:0
grf = GaussianRandomField(cov,KarhunenLoeve(300),pts)
grf = GaussianRandomField(cov,KarhunenLoeve(300),pts,nq=1000)
@test_throws ArgumentError GaussianRandomField(cov,KarhunenLoeve(300),pts,nq=100)

## non-SPD covariance matrix ##
cov = CovarianceFunction(2,SquaredExponential(0.5))
@suppress grf = GaussianRandomField(cov,KarhunenLoeve(200),pts1,pts2)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,SquaredExponential)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,K} where {C,K<:KarhunenLoeve})
@test length(grf.pts[1]) == length(pts1)
@test length(grf.pts[2]) == length(pts2)
@test size(sample(grf),1) == length(pts1)
@test size(sample(grf),2) == length(pts2)

## test domain with negative sides ##
cov = CovarianceFunction(2,Matern(1.,2.5,σ=1.,p=2))
pts1 = collect(-2:0.1:0)
pts2 = collect(-1:0.1:1)
grf = GaussianRandomField(cov,KarhunenLoeve(458),pts1,pts2)
@test isa(grf,GaussianRandomField)
@test isa(grf.cov,CovarianceFunction)
@test isa(grf.cov.cov,Matern)
@test ndims(grf.cov) == 2
@test isa(grf,GaussianRandomField{C,KarhunenLoeve{458}} where {C})
@test length(grf.pts[1]) == length(pts1)
@test length(grf.pts[2]) == length(pts2)
@test size(sample(grf),1) == length(pts1)
@test size(sample(grf),2) == length(pts2)

end