text:
# Basics of the mathematica SIR Model

using DifferentialEquations
using Plots

# This is a test for the Lorenz ODE system
function lorenz!(du,u,p,t)
 du[1] = 10.0*(u[2]-u[1])
 du[2] = u[1]*(28.0-u[3]) - u[2]
 du[3] = u[1]*u[2] - (8/3)*u[3]
end

u0 = [1.0;0.0;0.0]
tspan = (0.0,100.0)
prob = ODEProblem(lorenz!,u0,tspan)
sol = solve(prob)

plot(sol,vars=(1,2,3))

# Now the basic SIR ODE system

beta = 0.25
gamma = 0.07
N = 7000
function SIR_basic!(du,u,p,t)
    du[1] = -beta* u[1]/N*u[2]
    du[2] = beta* u[1]/N*u[2] - gamma*u[2]
    du[3] = gamma*u[2]
end

   u0 = [N-15.0;15.0;0.0]
   tspan = (0.0,100.0)
   prob = ODEProblem(SIR_basic!,u0,tspan)
   sol = solve(prob)

   plot(sol,vars=(1,2,3))


text:
using LinearAlgebra
using Test
using SparseArrays
using StaticArrays

tol = 5e2*eps()
N = 10

# scalar test
A = sparse(randn(N,N))
ATr_list = (Matrix(transpose(A)),)
function flux1(UL,UR)
    return UL .+ UR
end
U = (randn(N),)
rhs_exact = (A*U[1] + U[1].*sum(A,dims=2),)
rhs1, rhs2 = zero.(U), zero.(U)
hadamard_sum_ATr!(rhs1, ATr_list,flux1,U)
hadamard_sum_ATr!(rhs2, sparse.(ATr_list),flux1,U)
@test sum(norm.(rhs1 .- rhs_exact)) < tol
@test sum(norm.(rhs_exact .- rhs2)) < tol
@test sum(norm.(rhs1 .- rhs2)) < tol

# 2 flux test
A,B = randn(N,N),randn(N,N)
ATr_list = map(A->Matrix(transpose(A)),(A,B))
function flux2(UL,UR)
    uL,vL = UL
    uR,vR = UR
    return SVector{2}(uL+uR,0),SVector{2}(0,2*(vL+vR))
end
u,v = (randn(N),randn(N))
U = (u,v)
# computes sum(A.*F1 + B.*F2,dims=2)
#rhs_exact = (sum(A,dims=2), sum(B,dims=2))
rhs_exact = (A*u + u.*sum(A,dims=2),2*B*v + 2*v.*sum(B,dims=2))
# rhs_exact = (A*u + u.*sum(A,dims=2), B*v + v.*sum(B,dims=2))
rhs1, rhs2 = zero.(U), zero.(U)
hadamard_sum_ATr!(rhs1,ATr_list,flux2,U)
hadamard_sum_ATr!(rhs2,sparse(ATr_list[1]),(uL,uR)->flux2(uL,uR)[1],U)
hadamard_sum_ATr!(rhs2,sparse(ATr_list[2]),(uL,uR)->flux2(uL,uR)[2],U)

@test sum(norm.(rhs1 .- rhs_exact)) < tol
@test sum(norm.(rhs_exact .- rhs2)) < tol
@test sum(norm.(rhs1 .- rhs2)) < tol

ext: jl


text:
module TestPool
using StructDatabaseMapping
const Pool = StructDatabaseMapping.QueuePool
struct TestConnection
    a::Integer
end
function test_pool()
    pool = Pool(x->TestConnection(5))

end
end

max_stars_repo_name: lucianolorenti/StructDatabaseMapping

max_issues_repo_name: lucianolorenti/StructDatabaseMapping

max_forks_repo_name: lucianolorenti/StructDatabaseMapping


text:
const REGISTER_SIZE = 16
const REGISTER_COUNT = 16
const FP256 = true
const CACHELINE_SIZE = 64
const CACHE_SIZE = ((1<<15,1<<20,0))
const NUM_CORES = (Sys.CPU_THREADS)
const FMA = false
const FMA3 = false
const AVX2 = false
const AVX512F = false
const AVX512ER = false
const AVX512PF = false
const AVX512VL = false
const AVX512BW = false
const AVX512DQ = false
const AVX512CD = false
const SIMD_NATIVE_INTEGERS = true


text:
function read_nodal_segy(fname::String,
  nn::String,
  s::TimeSpec,
  t::TimeSpec,
  chans::ChanSpec,
  memmap::Bool)

  # Preprocessing
  (d0, d1) = parsetimewin(s, t)
  t0 = DateTime(d0).instant.periods.value*1000 - dtconst
  f = memmap ? IOBuffer(Mmap.mmap(fname)) : open(fname, "r")
  trace_fh = Array{Int16, 1}(undef, 3)
  shorts  = getfield(BUF, :int16_buf)
  fhd = Dict{String,Any}()
  # ww = Array{String, 1}(undef, 0)

  # ==========================================================================
  # Read file header
  fhd["filehdr"]  = fastread(f, 3200)
  fhd["jobid"]    = bswap(fastread(f, Int32))
  fhd["lineid"]   = bswap(fastread(f, Int32))
  fhd["reelid"]   = bswap(fastread(f, Int32))
  fast_readbytes!(f, BUF.buf, 48)
  fillx_i16_be!(shorts, BUF.buf, 24, 0)
  fastskip(f, 240)
  for i in 25:27
    shorts[i] = read(f, Int16)
  end
  fastskip(f, 94)

  # check endianness; can be inconsistent; 0x0400 is a kludge
  # (stands for SEG Y rev 10.24, current is rev 2.0)
  # if (unsigned(shorts[25]) > 0x0400) || (shorts[26] < zero(Int16)) || (shorts[26] > one(Int16)) || (shorts[27] < zero(Int16))
  #   shorts[25:27] .= bswap.(shorts[25:27])
  #   push!(ww, "Inconsistent file header endianness")
  # end
  # fastskip(f, 94)

  # Process
  nh = max(zero(Int16), getindex(shorts, 27))
  fhd["exthdr"] = Array{String,1}(undef, nh)
  [fhd["exthdr"][i] = fastread(f, 3200) for i in 1:nh]
  for (j,i) in enumerate(String[ "ntr", "naux", "filedt", "origdt", "filenx",
                                 "orignx", "fmt", "cdpfold", "trasort", "vsum",
                                 "swst", "swen0", "swlen", "swtyp", "tapnum",
                                 "swtapst", "swtapen", "taptyp", "corrtra", "bgainrec",
                                 "amprec", "msys", "zupdn", "vibpol", "segyver",
                                 "isfixed", "n_exthdr" ])
    fhd[i] = shorts[j]
  end

  nt = getindex(shorts, 1)
  trace_fh[1] = getindex(shorts,3)
  trace_fh[2] = getindex(shorts,5)
  trace_fh[3] = getindex(shorts,7)

  # ==========================================================================
  # initialize NodalData container; set variables
  if isempty(chans)
    chans = 1:nt
  end
  fs = sμ / Float64(trace_fh[1])
  data = Array{Float32, 2}(undef, trace_fh[2], nt)
  S = NodalData(data, fhd, chans, t0)
  net = nn * "."
  cha = string("..O", getbandcode(fs), "0")

  # ==========================================================================
  # Read traces
  j = 0

  for i = 1:nt
    C = do_trace(f, false, true, 0x00, true, trace_fh)
    if i in chans
      j += 1
      S.id[j]       = string(net, lpad(i, 5, '0'), cha)
      S.name[j]     = string(C.misc["rec_no"], "_", i)
      S.fs[j]       = C.fs
      S.gain[j]     = C.gain
      S.misc[j]     = C.misc
      S.t[j]        = C.t
      S.data[:, j] .= C.x
    end
  end

  # TO DO: actually use s, t here

  # ==========================================================================
  # Cleanup
  close(f)
  resize!(BUF.buf, 65535)
  fill!(S.fs, fs)
  fill!(S.src, realpath(fname))
  fill!(S.units, "m/m")

  # Output warnings to STDOUT
  # if !isempty(ww)
  #   for i = 1:length(ww)
  #     @warn(ww[i])
  #   end
  # end
  # S.info["Warnings"] = ww
  return S
end


text:
using LabelledArrays, Test, InteractiveUtils

@testset "Basic interface" begin
    x = @LSliced [1.0 2.0; 3.0 4.0; 5.0 6.0] (:a,:b,:c), (:x, :y)

    syms1 = (:a,:b,:c)
    syms2 = (:x, :y)

    for (i,s1) in enumerate(syms1), (j,s2) in enumerate(syms2)
        @show i,s2,j,s2
        @test x[i,j] == x[s1,s2]
    end

    f(x) = x[1,1]
    g(x) = x.a.x
    @time f(x)
    @time f(x)
    @time g(x)
    @time g(x)

    #@code_warntype x.a
    #@inferred getindex(x,:a)
    @code_warntype g(x)
    @inferred g(x)

    x = @LSliced [1 2; 3 4; 5 6] (:a,:b,:c), (:x, :y)
    @test x .* x isa LArray
    @test x .+ 1 isa LArray
    @test x .+ 1. isa LArray
    z = x .+ ones(Int, 3)
    @test z isa LArray && eltype(z) === Int
    z = x .+ ones(Float64, 3)
    @test z isa LArray && eltype(z) === Float64
    @test eltype(x .+ 1.) === Float64

    z = @LSliced Float64 (2,2) (:a,:b), (:c,:d)
    w = rand(2,2)
    z .= w

    @test z[:a,:c] == w[1,1]
    @test z[:b,:c] == w[2,1]
    @test z[:a,:d] == w[1,2]
    @test z[:b,:d] == w[2,2]

    vals = [1 2; 3 4; 5 6]
    syms = (:a,:b,:c), (:x, :y)
    x = @LSliced vals syms
    x.a.x = 33
    x.b.y = 99
    x[:a,:y] = 44
    x[:b,:x] = 77
    @test x[1,1] == 33
    @test x[2,2] == 99
    @test x[1,2] == 44
    @test x[2,1] == 77
end

@testset "Alternate array backends" begin
    v = view([1 2 3; 4 5 6; 7 8 9], 2:3, 2:3)
    x = LArray{Tuple{(:a,:b),(:c,:d)}}(v)
    @test x.a.c == 5
    @test x[:b, :d] == 9

    s = similar(x)
    @test size(s) == size(x)
    @test typeof(s.__x) == Array{Int64,2}
    @test LabelledArrays.symnames(typeof(s)) == ((:a,:b), (:c,:d))
end

@testset "undef copy" begin
    z = @LSliced Float64 (4,2) (:a,:b,:c,:d), (:x, :y)
    t = similar(z, String) # t's elements are uninitialized
    @test_throws UndefRefError t[1]
    copy(t) # should be ok
    deepcopy(t) # should also be ok
end

@testset "accessing labels, i.e. symbols" begin
    A2 = @LSliced [1 2; 3 4; 5 6] (:a,:b,:c), (:x, :y)
    ret = symbols(A2)
    @test ret == ((:a, :b, :c), (:x, :y))
    #
    # higher dimensional still returns a tuple-Type object
    A3 = @LSliced cat(A2,10*A2;dims = 3) (:a,:b,:c), (:x, :y), (:u,:v)
    ret = symbols(A3)
    @test ret == Tuple{(:a, :b, :c),(:x, :y),(:u, :v)}
end

# @testset "accessing dimensions symbols of LSliced" begin
#     A = @LSliced [1 2; 3 4; 5 6] (:a,:b,:c), (:x, :y)
#     dim = 1
#     rows = dimSymbols(A,dim)
#     @test rows == (:a,:b,:c)
#     cols = dimSymbols(A,2)
#     @test cols == (:x,:y)
#     rows = rowSymbols(A)
#     @test rows == (:a,:b,:c)
#     cols = colSymbols(A)
#     @test cols == (:x,:y)
# end


text:
acceptAll(x) = true

"""
Fill volume with spheres of certain radius in hexagonal packing.

In brief, the model generates a box with spheres in hexagonal packaging, removes the ones outside the volume.
The positions are perturbed by the noise term.

# Parameters
 - **model** (Model) Compiled agent based model
 - **box** (Array{Float,2}) Maximum box where to fill the spheres.
 - **r** (Number) Radius of the spheres

# Optional keyword parameters
    - **fExtrude** (Function) Function that returns true if center of the sphere is inside the volume. Helpts to define non-cubic shapes
    - **N** (Int) Maximum number of particles inside the volume. If NaN (default), there is not upper bound.
    - **mediumN** (Array{Int,1) Grid dimensions of medium if Medium is declared.

![Figure](assets/FillVolumeSpheres.png)
Figure rendered with [Makie.jl](https://github.com/JuliaPlots/Makie.jl) using meshscatter function.
"""
function initialiseCommunityCompactHexagonal(model,box,r;fExtrude::Function=acceptAll,N=NaN,mediumN::Array{Int,1}=Array{Int,1}([]))

    if model.agent.dims != 3
        error("Initialiser only works with 3D models.")
    end

    X = compactHexagonal(box,r)

    #Extrude
    ext = fExtrude.([X[i,:] for i in 1:size(X)[1]])
    X = X[ext,:]

    #Remove exceding cells
    if N != NaN
        l = size(X)[1]
        while size(X)[1] > N
            rem = rand(1:length(volumeX),l-N)
            nrem = [i for i in 1:l if !(i in rem)]
            X = X[nrem,:];
            l = size(X)[1]
        end
    end

    com = Community(model,N=size(X)[1],mediumN=mediumN)
    com.x = X[:,1]
    com.y = X[:,2]
    com.z = X[:,3]

    return com
end

"""
Fill volume with spheres of certain radius in cubic packing.

In brief, the model generates a box with spheres in hexagonal packaging, removes the ones outside the volume.
The positions are perturbed by the noise term.

# Parameters
 - **model** (Model) Compiled agent based model
 - **box** (Array{Float,2}) Maximum box where to fill the spheres.
 - **r** (Number) Radius of the spheres

# Optional keyword parameters
    - **fExtrude** (Function) Function that returns true if center of the sphere is inside the volume. Helpts to define non-cubic shapes
    - **N** (Int) Maximum number of particles inside the volume. If NaN (default), there is not upper bound.
    - **mediumN** (Array{Int,1) Grid dimensions of medium if Medium is declared.

![Figure](assets/FillVolumeSpheres.png)
Figure rendered with [Makie.jl](https://github.com/JuliaPlots/Makie.jl) using meshscatter function.
"""
function initialiseCommunityCompactCubic(model,box,r;fExtrude=acceptAll,N=NaN,mediumN::Array{Int,1}=Array{Int,1}([]))

    if model.agent.dims != 3
        error("Initialiser only works with 3D models.")
    end

    X = cubic(box,r)

    #Extrude
    ext = fExtrude.([X[i,:] for i in 1:size(X)[1]])
    X = X[ext,:]

    #Remove exceding cells
    if N != NaN
        l = size(X)[1]
        while size(X)[1] > N
            rem = rand(1:length(volumeX),l-N)
            nrem = [i for i in 1:l if !(i in rem)]
            X = X[nrem,:];
            l = size(X)[1]
        end
    end

    com = Community(model,N=size(X)[1],mediumN=mediumN)
    com.x = X[:,1]
    com.y = X[:,2]
    com.z = X[:,3]

    return com
end


text:
"""
Type for storing data acquired with a steady state response (SSR) experimental paradigm.

In addition to the functions available for all EEG types,
the SSR type supports:

* `modulationrate()`

The following standard names are used when saving data to the processing dictionary.

* `Name`: The identifier for the participant
* `Side`: Side of stimulation
* `Carrier_Frequency`: Carrier frequency of the stimulus
* `Amplitude`: Amplitude of the stimulus
* `epochs`: The epochs extracted from the recording
* `sweeps`: The extracted sweeps from the recording

#### Example

Put an example here

```julia
s = read_SSR("filename")
s.modulationrate = 40.0391u"Hz"
s = rereference(s, "Cz")
```

"""
mutable struct SSR <: EEG
    data::Array
    sensors::Array{Sensor}
    triggers::Dict
    system_codes::Dict
    samplingrate::typeof(1.0u"Hz")
    modulationrate::typeof(1.0u"Hz")
    reference_channel::Array{AbstractString,1}
    file_path::AbstractString
    file_name::AbstractString
    processing::Dict
    header::Dict
end

#######################################
#
# SSR info
#
#######################################

"""
Return the modulation rate of a steady state type.
If no type is provided, the modulation rate is returned as a floating point.

#### Example

Return the modulation rate of a recording

```julia
s = read_SSR(filename)
modulationrate(s)
```
"""
modulationrate(t, s::SSR) = convert(t, ustrip(s.modulationrate))
modulationrate(s::SSR) = modulationrate(AbstractFloat, s)

#######################################
#
# Show
#
#######################################

import Base.show
function Base.show(io::IO, a::SSR)
    time_length = round.(size(a.data, 1) / samplingrate(a) / 60)
    println(
        io,
        "SSR measurement of $time_length mins with $(size(a.data,2)) channels sampled at $(a.samplingrate)",
    )
    println(io, "  Modulation frequency: $(a.modulationrate )")

    if haskey(a.processing, "Amplitude")
        println(io, "  Stimulation amplitude: $(a.processing["Amplitude"]) dB")
    end
    if haskey(a.processing, "Name")
        println(io, "  Participant name: $(a.processing["Name"] )")
    end
    if haskey(a.processing, "Side")
        println(io, "  Stimulation side: $(a.processing["Side"] )")
    end
    if haskey(a.processing, "Carrier_Frequency")
        println(io, "  Carrier frequency: $(a.processing["Carrier_Frequency"] ) Hz")
    end

end

#######################################
#
# Helper functions
#
#######################################

function assr_frequency(
    rounded_freq::Number;
    stimulation_samplingrate::Number = 32000,
    stimulation_frames_per_epoch::Number = 32768,
)

    round.(rounded_freq / (stimulation_samplingrate / stimulation_frames_per_epoch)) *
    stimulation_samplingrate / stimulation_frames_per_epoch
end

function assr_frequency(rounded_freq::AbstractVector)

    [assr_frequency(f) for f in rounded_freq]
end


text:
begin
    function format_new_call_expr(expr, orig_expr, call, orig_call)
        lparen_pos = line_pos(orig_call, first(children(orig_call)[2].span))
        # Heuristic: If the arugment indents (more spefically the first indented line), is further
        # left than the original lparen, don't unindent, it's unlikely they were aligned.
        if lparen_pos != 0
            indents = countindent_body(orig_expr)
            if !isempty(indents) && (indents[1] < lparen_pos)
                return (false, call)
            end
        end
        # A bit cheesy, but w/e
        buf = IOBuffer()
        print_replacement(buf, expr, false, false)
        x = String(take!(buf))
        nx = sizeof(x)-length(orig_expr.span)
        return true, format_addindent_body(call, nx, nothing, false)
    end

    # Special purpose formatter for `@compat` results
    function format_result(expr, orig_expr)
        CSTParser.defines_function(expr) || return expr
        call_idx = isexpr(expr, FunctionDef) ? 2 : 1
        callorwhere = children(expr)[call_idx]
        call = is_where_expr(callorwhere) ? children(callorwhere)[1] : callorwhere
        ret_expr = ChildReplacementNode(nothing, collect(children(expr)), expr)
        orig_call = children(children(orig_expr)[2])[call_idx]
        orig_call = is_where_expr(orig_call) ? children(orig_call)[1] : orig_call
        ok, call′ = format_new_call_expr(expr, orig_expr, call, orig_call)
        ok || return expr
        if is_where_expr(callorwhere)
            children(ret_expr)[call_idx] = ChildReplacementNode(ret_expr, collect(children(callorwhere)), callorwhere)
            children(children(ret_expr)[call_idx])[1] = call′
        else
            children(ret_expr)[call_idx] = call′
        end
        ret_expr
    end

    struct ObsoleteCompatMacro
        min_ver::VersionNumber
    end
    register(ObsoleteCompatMacro, Deprecation(
        "This compat macro is no longer required",
        "julia",
        v"0.5.0", v"0.5.0", typemax(VersionNumber)
    ))

    ObsoleteCompatMacro() = ObsoleteCompatMacro(v"0.6.0")
    function dep_for_vers(::Type{ObsoleteCompatMacro}, vers)
        ObsoleteCompatMacro(minimum(map(interval->interval.lower, vers["julia"].intervals)))
    end

    # This follows the structure of the _compat function from Compat.jl
    # to determine if the Compat macro would have done something
    istopsymbol(ex, mod, sym) = ex in (sym, Expr(:(.), mod, Expr(:quote, sym)))
    _compat(min_ver, ex) = false
    function new_style_typealias(ex::ANY)
        Base.Meta.isexpr(ex, :(=)) || return false
        ex = ex::Expr
        return length(children(ex)) == 2 && Base.Meta.isexpr(children(ex)[1], :curly)
    end
    is_index_style(ex::Expr) = ex == :(Compat.IndexStyle) || ex == :(Base.IndexStyle) ||
    (ex.head == :(.) && (children(ex)[1] == :Compat || children(ex)[1] == :Base) &&
        children(ex)[2] == Expr(:quote, :IndexStyle))
    is_index_style(arg) = false
    withincurly(ex) = Base.Meta.isexpr(ex, :curly) ? children(ex)[1] : ex
    function _compat(min_ver, ex::Expr)
        if ex.head === :call
            f = children(ex)[1]
            if min_ver < v"0.6.0-dev.826" && length(children(ex)) == 3 && # julia#18510
                    istopsymbol(withincurly(children(ex)[1]), :Base, :Nullable)
                return true
            end
        elseif ex.head === :curly
            f = children(ex)[1]
            if min_ver < v"0.6.0-dev.2575" && any(i->isa(children(ex)[i], Expr) && children(ex)[i].head == Symbol("<:"), 2:length(children(ex))) #20414
                return true
            end
        elseif ex.head === :quote && isa(children(ex)[1], Symbol)
            # Passthrough
            return false
        elseif min_ver < v"0.6.0-dev.2782" && new_style_typealias(ex)
            return true
        elseif min_ver < v"0.6.0-dev.2782" && ex.head === :const && length(children(ex)) == 1 && new_style_typealias(children(ex)[1])
            return true
        end
        if min_ver < v"0.6.0-dev.2840"
            if ex.head == :(=) && isa(children(ex)[1], Expr) && children(ex)[1].head == :call
                a = children(children(ex)[1])[1]
                if is_index_style(a)
                    return true
                elseif isa(a, Expr) && a.head == :curly
                    if is_index_style(children(a)[1])
                        return true
                    end
                end
            end
        end
        if min_ver < v"0.7.0-DEV.880"
            if ex.head == :curly && children(ex)[1] == :CartesianRange && length(children(ex)) >= 2
                a = children(ex)[2]
                if a != :CartesianIndex && !(isa(a, Expr) && a.head == :curly && children(a)[1] == :CartesianIndex)
                    return true
                end
            end
        end
        if min_ver < v"0.7.0-DEV.2562"
            if ex.head == :call && children(ex)[1] == :finalizer
                return true
            end
        end
        return any(x->_compat(min_ver, x), children(ex))
    end

    function is_at_compat_needed(dep, expr)
        if isexpr(expr, CSTParser.Primitive) || isexpr(expr, CSTParser.Abstract)
            return dep.min_ver < v"0.6.0-dev.2746"
        else
            _compat(dep.min_ver, Expr(expr.expr))
        end
    end

    match(ObsoleteCompatMacro, CSTParser.MacroCall) do x
        dep, expr, resolutions, context = x

        is_macroname(expr, "compat") || return
        args = filter(x->!isexpr(x, CSTParser.PUNCTUATION), children(expr)[2:end])
        length(args) == 1 || return

        if !is_at_compat_needed(dep, args[1])
            buf = IOBuffer()
            print_replacement(buf, format_result(args[1], expr), false, false)
            push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
        end
    end

    struct ObsoleteModuleRef
        min_ver::VersionNumber
    end
    register(ObsoleteModuleRef, Deprecation(
        "This reference to a module that used to be in Base (or Compat) can now be directly referenced.",
        "julia",
        v"0.5.0", v"0.5.0", typemax(VersionNumber)
    ))
    ObsoleteModuleRef() = ObsoleteModuleRef(v"0.6.0")
    function dep_for_vers(::Type{ObsoleteModuleRef}, vers)
        ObsoleteModuleRef(minimum(map(interval->interval.lower, vers["julia"].intervals)))
    end

    function should_replace(dep, refed)
        if (is_identifier(refed, "Test") ||
            is_identifier(refed, "Mmap") ||
            is_identifier(refed, "DelimitedFiles") ||
            is_identifier(refed, "SharedArrays")) && dep.min_ver > v"0.7.0-DEV.2005"
        elseif is_identifier(refed, "Dates") && dep.min_ver > v"0.7.0-DEV.2575"
        elseif is_identifier(refed, "Libdl") && dep.min_ver > v"0.7.0-DEV.3382"
        elseif is_identifier(refed, "Printf") && dep.min_ver > v"0.7.0-DEV.3052"
        elseif is_identifier(refed, "Markdown") && dep.min_ver > v"0.7.0-DEV.3589"
        elseif is_identifier(refed, "Random") && dep.min_ver > v"0.7.0-DEV.3406"
        else
            return false
        end
        return true
    end

    function process_module_getfield(x)
        dep, expr, resolutions, context = x

        isexpr(children(expr)[2], OPERATOR, Tokens.DOT) || return

        (is_identifier(children(expr)[1], "Compat") ||
         is_identifier(children(expr)[1], "Base")) || return

        refed = children(expr)[3]
        should_replace(dep, refed) || return

        buf = IOBuffer()
        print_replacement(buf, format_result(args[1], expr), false, false)
        push!(resolutions, TextReplacement(dep, expr.span, String(take!(buf))))
    end

    function process_module_using(x)
        dep, expr, resolutions, context = x
        i = 2
        while i < length(children(expr))
            id = children(expr)[i]
            i += 1
            if is_identifier(id, "Compat") || is_identifier(id, "Base")
                if !isexpr(children(expr)[i], PUNCTUATION, Tokens.DOT)
                    i += 1
                    continue
                end
                i += 1
                rid = children(expr)[i]
                if should_replace(dep, rid)
                    buf = IOBuffer()
                    print_replacement(buf, rid, false, false)
                    push!(resolutions, TextReplacement(dep,
                        first(id.span):last(rid.span), String(take!(buf))))
                end
            end
            while i < length(children(expr)) && !isexpr(children(expr)[i], PUNCTUATION, Tokens.COMMA)
                i += 1
            end
            i += 1
        end
    end

    match(ObsoleteModuleRef, CSTParser.BinarySyntaxOpCall) do x
        process_module_getfield(x)
    end

    match(ObsoleteModuleRef, CSTParser.Using) do x
        process_module_using(x)
    end

    match(ObsoleteModuleRef, CSTParser.Import) do x
        process_module_using(x)
    end
end


text:
"""
Ciphertext stealing

- Swaps the final two ciphertext blocks unconditionally [^1]
- Uses standard CBC interface [^2]
- In case the plaintext is one block long or less, the IV acts as the prior block and is modified

[^1]: https://en.wikipedia.org/wiki/Ciphertext_stealing#CS3
[^2]: https://en.wikipedia.org/wiki/Ciphertext_stealing#CBC_implementation_notes
"""
module CiphertextStealing

using Nettle: encrypt, decrypt, Encryptor, Decryptor
using PaddedViews: PaddedView

export aes_128_cbc_encode_cts!, aes_128_cbc_decode_cts!

#===== Helpers =====#

function swap!(a::AbstractVector{T}, b::AbstractVector{T}) where T
    for i in eachindex(a, b)
        @inbounds a[i], b[i] = b[i], a[i]
    end
end

function aes_128_ecb_encode(
    plaintext::AbstractVector{UInt8},
    key::AbstractVector{UInt8}
)::Vector{UInt8}
    @assert length(plaintext) % 16 == 0 && length(key) == 16
    enc = Encryptor("AES128", key)
    encrypt(enc, plaintext)
end

function aes_128_ecb_decode(
    ciphertext::AbstractVector{UInt8},
    key::AbstractVector{UInt8}
)::Vector{UInt8}
    @assert length(ciphertext) % 16 == 0 && length(key) == 16
    dec = Decryptor("AES128", key)
    decrypt(dec, ciphertext)
end

#===== The standard AES-128 CBC routine =====#

function aes_128_cbc_encode(
    plaintext::AbstractVector{UInt8},
    key::Vector{UInt8},
    iv::Vector{UInt8},
)::Vector{UInt8}
    @assert length(plaintext) % 16 == 0 && length(key) == 16 && length(iv) == 16

    ciphertext = UInt8[]
    sizehint!(ciphertext, length(plaintext))

    last_block = copy(iv)
    for b in Iterators.partition(plaintext, 16)
        last_block .⊻= b
        last_block = aes_128_ecb_encode(last_block, key)
        append!(ciphertext, last_block)
    end

    ciphertext
end

function aes_128_cbc_decode(
    ciphertext::AbstractVector{UInt8},
    key::Vector{UInt8},
    iv::Vector{UInt8},
)::Vector{UInt8}
    @assert length(ciphertext) % 16 == 0 && length(key) == 16 && length(iv) == 16

    plaintext = UInt8[]
    sizehint!(plaintext, length(ciphertext))

    last_block = iv
    for b in Iterators.partition(ciphertext, 16)
        dec_block = aes_128_ecb_decode(b, key) .⊻ last_block
        append!(plaintext, dec_block)
        last_block = b
    end

    plaintext
end

#===== Ciphertext Stealing =====#

# It will modify IV if the plaintext is one block or shorter.
function aes_128_cbc_encode_cts!(
    plaintext::AbstractVector{UInt8},
    key::Vector{UInt8},
    iv::Vector{UInt8},
)::Vector{UInt8}
    len = length(plaintext)
    @assert len > 0
    padded_len = len + mod(-len, 16)
    ciphertext = aes_128_cbc_encode(PaddedView(UInt8(0), plaintext, (padded_len,)), key, iv)

    ultimate = @view ciphertext[(end - 16 + 1):end]
    penultimate = @views len > 16 ? ciphertext[(end - 16 * 2 + 1):(end - 16)] : iv[:]
    swap!(ultimate, penultimate)
    resize!(ciphertext, len)

    ciphertext
end

# It will modify the penultimate block or (if the ciphertext is one block or shorter) IV.
function aes_128_cbc_decode_cts!(
    ciphertext::AbstractVector{UInt8},
    key::Vector{UInt8},
    iv::Vector{UInt8},
)::Vector{UInt8}
    len = length(ciphertext)
    @assert len > 0
    padding_len = mod(-len, 16)
    ultimate_len = 16 - padding_len
    padded_len = len + padding_len

    cipher_ultimate = @view ciphertext[(end - ultimate_len + 1):end]
    cipher_penultimate =
        if len > 16
            @view ciphertext[(end - ultimate_len - 16 + 1):(end - ultimate_len)]
        else
            @view iv[:]
        end
    decrypt_ultimate = aes_128_ecb_decode(cipher_penultimate, key)  # Performance penalty

    cipher_penultimate[1:ultimate_len] .= cipher_ultimate
    cipher_penultimate[(end - padding_len + 1):end] .= decrypt_ultimate[(end - padding_len + 1):end]

    plaintext = aes_128_cbc_decode(PaddedView(UInt8(0), ciphertext, (padded_len,)), key, iv)
    resize!(plaintext, len)
    plaintext[(end - ultimate_len + 1):end] .= cipher_ultimate .⊻ decrypt_ultimate[1:ultimate_len]

    plaintext
end

end # module

max_stars_repo_name: gyk/TrivialSolutions

max_issues_repo_name: gyk/TrivialSolutions

max_forks_repo_name: gyk/TrivialSolutions


text:
## Utils
instrument_reset(obj::Instrument)    = write(obj, "*RST")
instrument_clear(obj)    = write(obj, "*CLS")
instrument_get_id(obj)   = query(obj, "*IDN?")
instrument_beep_on(obj)  = write(obj,"SYST:BEEP:STAT on")
instrument_beep_off(obj) = write(obj,"SYST:BEEP:STAT off")
instrument_query_complete(obj) = query(obj,"*OPC?")
instrument_wait_complete(obj)  = query(obj,"*WAI")
instrument_query_error(obj)    = query(obj,"SYST:ERR?")
function instrument_empty_buffer(instr)
    dummy = read(instr)
    return nothing
end

## TRIGGER
instrument_trigger_set_source(obj;ch="1",sour="CHAN1")    = write(obj,"TRIG$ch:SOUR $sour")
instrument_trigger_set_delay(obj;ch="1",del=0)          = write(obj,"TRIG$ch:DEL $del")
instrument_trigger_set_slope(obj;ch="1",slp="positive") = write(obj,"TRIG$ch:SLOP $slp")
instrument_trigger_set_out_stat(obj;st="normal")        = write(obj,"OUTP:TRIG $st")
instrument_set_out_trig_slope(obj;slp="positive")       = write(obj,"OUTP:TRIG:SLOP $slp")
instrument_send_software_trigger(obj)                   = write(obj,"*TRG")

## OUTPUT FUNCTION
instrument_set_waveform(obj;ch="1",func="sin") = write(obj,"SOURce$ch:FUNC $func")
instrument_set_arbitrary_waveform(obj;ch="1",arb="mywfm") = write(obj,"SOURce$ch:FUNCtion:ARBitrary $arb")
instrument_set_amplitude(obj;ch="1",unit="vpp",volt=0.01) = write(obj,"SOURce$ch:VOLT:UNIT $unit\n;SOURce$ch:VOLT $volt")
instrument_set_offset(obj;ch="1",volt=0.01) = write(obj,"SOURce$ch:VOLT:OFFset $volt")
instrument_set_hilevel(obj;ch="1",volt=0.01) = write(obj,"SOURce$ch:VOLT:HIGH $volt")
instrument_set_lolevel(obj;ch="1",volt=0.01) = write(obj,"SOURce$ch:VOLT:LOW $volt")
instrument_set_frequency(obj;ch="1",fr=1000) = write(obj,"SOURce$ch:FREQ $fr")
instrument_set_period(obj;ch="1",func="sin",per=0.001) = write(obj,"SOURce$ch:FUNC:$func:PER $per")
instrument_set_phase(obj;ch="1",unit="deg",phase=90) = write(obj,"UNIT:ANGL $unit\n;SOURce$ch:PHAS $phase")
instrument_set_symmetry(obj;ch="1",func="ramp",sym=50) = write(obj,"SOURce$ch:FUNC:$func:SYMM $sym")
instrument_set_both_edge_times(obj;ch="1",func="pulse",dt=10e-9) = write(obj,"SOURce$ch:FUNC:$func:TRAN $dt")
#instrument_set_lead_edge_time(obj;ch="1",func="pulse",dt=10e-9) where {T<:F335x2} = write(obj,"SOURce$ch:FUNC:$func:TRAN:LEAD $dt")
#instrument_set_trail_edge_time(obj;ch="1",func="pulse",dt=10e-9) where {T<:F335x2} = write(obj,"SOURce$ch:FUNC:$func:TRAN:TRA $dt")
instrument_set_width_time(obj;ch="1",func="pulse",wd=1e-4) = write(obj,"SOURce$ch:FUNC:$func:WIDTh $wd")

## OUTPUT
instrument_set_load(obj;ch="1",load="INF") = write(obj,"OUTP$ch:LOAD $load")
instrument_set_polarity(obj;ch="1",pol="normal") = write(obj,"OUTP$ch:POL $pol")
instrument_set_output(obj;ch="1",st="off") = write(obj,"OUTP$ch $st")

## SYNC
instrument_sync_state_on(obj;ch="1") = write(obj,"OUTP$ch:SYNC on")
instrument_sync_state_off(obj;ch="1") = write(obj,"OUTP$ch:SYNC off")

## ARBITRARY WAVEFORMS
instrument_abritrary_waveform_sample_rate(obj;ch="1",srate=1e6) = write(obj,"SOURce$ch:FUNC:ARB:SRAT $srate")
instrument_arbitrary_waveforms_clear(obj;ch="1") = write(obj,"SOURce$ch:DATA:VOLatile:CLEar")


text:
struct DiffEqNLSolveTag end

struct DiffCache{T<:AbstractArray, S<:AbstractArray}
    du::T
    dual_du::S
end

Base.@pure function DiffCache(T, size, ::Type{Val{chunk_size}}) where chunk_size
  DiffCache(fill(zero(T), size...), fill(zero(Dual{typeof(ForwardDiff.Tag(DiffEqNLSolveTag(),T)),T,chunk_size}), size...))
end

Base.@pure DiffCache(u::AbstractArray) = DiffCache(eltype(u),size(u),Val{ForwardDiff.pickchunksize(length(u))})
Base.@pure DiffCache(u::AbstractArray,nlsolve) = DiffCache(eltype(u),size(u),Val{get_chunksize(nlsolve)})
Base.@pure DiffCache(u::AbstractArray,T::Type{Val{CS}}) where {CS} = DiffCache(eltype(u),size(u),T)

get_du(dc::DiffCache, ::Type{T}) where {T<:Dual} = dc.dual_du
get_du(dc::DiffCache, T) = dc.du

# Default nlsolve behavior, should move to DiffEqDiffTools.jl

Base.@pure determine_chunksize(u,alg::DiffEqBase.DEAlgorithm) = determine_chunksize(u,get_chunksize(alg))
Base.@pure function determine_chunksize(u,CS)
  if CS != 0
    return CS
  else
    return ForwardDiff.pickchunksize(length(u))
  end
end

struct NLSOLVEJL_SETUP{CS,AD} end
Base.@pure NLSOLVEJL_SETUP(;chunk_size=0,autodiff=true) = NLSOLVEJL_SETUP{chunk_size,autodiff}()
(::NLSOLVEJL_SETUP)(f,u0; kwargs...) = (res=NLsolve.nlsolve(f,u0; kwargs...); res.zero)
function (p::NLSOLVEJL_SETUP{CS,AD})(::Type{Val{:init}},f,u0_prototype) where {CS,AD}
  AD ? autodiff = :forward : autodiff = :central
  OnceDifferentiable(f, u0_prototype, u0_prototype, autodiff,
                     ForwardDiff.Chunk(determine_chunksize(u0_prototype,CS)))
end

get_chunksize(x) = 0
get_chunksize(x::NLSOLVEJL_SETUP{CS,AD}) where {CS,AD} = CS

"""
    calculate_residuals!(out, ũ, u₀, u₁, α, ρ)

Save element-wise residuals
```math
\\frac{ũ}{α+\\max{|u₀|,|u₁|}*ρ}
```
in `out`.
"""
@inline @muladd function calculate_residuals!(out, ũ, u₀, u₁, α, ρ, internalnorm)
    @. out = ũ / (α + max(internalnorm(u₀), internalnorm(u₁)) * ρ)
end

@inline @muladd function calculate_residuals!(out::Array{T}, ũ::Array{T}, u₀::Array{T},
                                              u₁::Array{T}, α::T2, ρ::Real,
                                              internalnorm) where {T<:Number,T2<:Number}
    @tight_loop_macros for i in eachindex(out)
        @inbounds out[i] = ũ[i] / (α + max(internalnorm(u₀[i]), internalnorm(u₁[i])) * ρ)
    end
end

"""
    calculate_residuals!(out, u₀, u₁, α, ρ)

Save element-wise residuals
```math
\\frac{ũ}{α+\\max{|u₀|,|u₁|}*ρ}
```
in `out`.
"""
@inline @muladd function calculate_residuals!(out, u₀, u₁, α, ρ, internalnorm)
    @. out = (u₁ - u₀) / (α + max(internalnorm(u₀), internalnorm(u₁)) * ρ)
end

@inline @muladd function calculate_residuals!(out::Array{T}, u₀::Array{T},
                                              u₁::Array{T}, α::T2, ρ::Real,
                                              internalnorm) where {T<:Number,T2<:Number}
    @tight_loop_macros for i in eachindex(out)
        @inbounds out[i] = (u₁[i] - u₀[i]) / (α + max(internalnorm(u₀[i]), internalnorm(u₁[i])) * ρ)
    end
end

"""
    calculate_residuals(ũ, u₀, u₁, α, ρ)

Calculate element-wise residuals
```math
\\frac{ũ}{α+\\max{|u₀|,|u₁|}*ρ}
```
"""
@inline @muladd function calculate_residuals(ũ, u₀, u₁, α, ρ, internalnorm)
    ũ / (α + max(internalnorm(u₀), internalnorm(u₁)) * ρ)
end

@inline @muladd function calculate_residuals(ũ::Array{T}, u₀::Array{T}, u₁::Array{T}, α::T2,
                                             ρ::Real, internalnorm) where
                                             {T<:Number,T2<:Number}
    out = similar(ũ)
    calculate_residuals!(out, ũ, u₀, u₁, α, ρ, internalnorm)
    out
end

"""
    calculate_residuals(u₀, u₁, α, ρ)

Calculate element-wise residuals
```math
\\frac{ũ}{α+\\max{|u₀|,|u₁|}*ρ}
```
"""
@inline @muladd function calculate_residuals(u₀, u₁, α, ρ, internalnorm)
    (u₁ - u₀) / (α + max(internalnorm(u₀), internalnorm(u₁)) * ρ)
end

@inline @muladd function calculate_residuals(u₀::Array{T}, u₁::Array{T}, α::T2,
                                             ρ::Real, internalnorm) where
                                             {T<:Number,T2<:Number}
    out = similar(u₀)
    calculate_residuals!(out, u₀, u₁, α, ρ, internalnorm)
    out
end

macro swap!(x,y)
  quote
    local tmp = $(esc(x))
    $(esc(x)) = $(esc(y))
    $(esc(y)) = tmp
  end
end

islinear(f) = f isa DiffEqBase.AbstractDiffEqLinearOperator && f.update_func === DEFAULT_UPDATE_FUNC


text:
#########################################################################
# Insolation
#########################################################################
#=
Code originally from:
https://github.com/climate-machine/RRTMGP.jl/blob/master/src/rte/SolarZenithAngle.jl

with examples at:
https://climate-machine.github.io/RRTMGP.jl/latest/RTE/SolarZenithAngle.html

The physics and details about the algorithm
can be found in Tapio Schneider’s book Physics of the Earth's Climate, section 3.4
http://climate-dynamics.org/wp-content/uploads/2017/04/Climate_Book.pdf
=#
export insolation, monthly_insolation
"""
    insolation(t, ϕ; kwargs...)
Calculate daily averaged insolation in W/m² at given time and latitude `t, φ`.
`φ` is given in **degrees**, and `t` in **days** (real number or date).

Keywords:
```
Ya = DAYS_IN_YEAR # = 365.26 # days
t_VE = 76.0 # days of vernal equinox
S_0 = 1362.0 # W/m^2
γ=23.44
ϖ=282.95
e=0.017 # eccentricity
```
"""
insolation(t::TimeType, φ; kwargs...) = insolation(float(dayofyear(t)), φ; kwargs...)
function insolation(t::Real, ϕ;
        Ya = DAYS_IN_YEAR, # = 365.26 # days
        t_VE = 76.0, # days since Jan 1
        S_0 = 1362.0, # W/m^2
        γ=23.44,
        ϖ=282.95,
        e=0.017,
    )

    # convert inputs from degrees to radians
    ϕ = ϕ * π / 180
    γ = γ * π / 180
    ϖ = ϖ * π / 180

    # step 1, calculate the mean anomaly at vernal equinox
    β = sqrt(1 - e^2)
    M_VE = -ϖ + (e + e^3 / 4) * (1 + β) * sin(ϖ)

    # step 2, calculate the mean anomaly
    M = (2 * π * (t - t_VE)) / (Ya) + M_VE

    # step 3, calculate the true anomaly
    A = M + (2 * e - e^3 / 4) * sin(M)

    # step 4, calculate the distance to the sun
    d = (1 - e^2) / (1 + e * cos(A))

    # step 5, calculate the solar longitude
    L_s = A + ϖ

    # step 6, calculate the declination angle
    δ = asin(sin(γ) * sin(L_s))

    # step 7, calculate the sunrise/sunset angle
    T = tan(ϕ) * tan(δ)
    if T >= 1
        η_d = π
    elseif T <= -1
        η_d = 0.0
    else
        η_d = acos(-1 * T)
    end

    # step 8, calculate the daily averaged cos(zenith angle)
    c1 = η_d * sin(ϕ) * sin(δ)
    c2 = cos(ϕ) * cos(δ) * sin(η_d)
    cosbar = (1 / π) * (c1 + c2)

    # step 9, calculate the flux
    F = S_0 * (1 / d)^2 * cosbar
    return F
end


text:
# module for playing monte carlo localization by particle filter
# particles are resampled by random sampling

module AnimeKldMcl
  include(joinpath(split(@__FILE__, "src")[1], "src/model/world/world.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/model/robot/differential_wheeled_robot/differential_wheeled_robot.jl"))
  include(joinpath(split(@__FILE__, "src")[1], "src/localization/particle_filter/kld_sampling/kld_mcl.jl"))

  function animate_per_time(time_interval, world, map, robot)
    draw(world)
    annotate!(-3.5, 4.5, "t = $(time_interval)", "black")
    draw!(map)
    draw!(robot)
  end

  function main(time_interval=0.1; is_test=false)
    w = World(-5.0, 5.0, -5.0, 5.0)

    m = Map()
    add_object(m, Object(2.0, -3.0, id=1))
    add_object(m, Object(3.0, 3.0, id=2))

    s = Sensor(m, dist_noise_rate=0.1, dir_noise=pi/90,
               dist_bias_rate_stddev=0.1, dir_bias_stddev=pi/90)

    init_pose = [0.0, 0.0, 0.0]
    e = KldMcl(init_pose, 1000, env_map=m)
    circling = Agent(0.2, 10.0/180*pi, estimator=e)
    r = DifferentialWheeledRobot(init_pose, 0.2, "black",
                                 circling, time_interval,
                                 noise_per_meter=5,
                                 noise_std=pi/30,
                                 bias_rate_stds=[0.1, 0.1],
                                 sensor=s)

    if is_test == false
      anime = @animate for t in 0:time_interval:30
        animate_per_time(t, w, m, r)
      end
      path = "src/localization/particle_filter/kld_sampling/anime_kld_mcl.gif"
      gif(anime, fps=15, joinpath(split(@__FILE__, "src")[1], path))
    else
      for t in 0:time_interval:10
        animate_per_time(t, w, m, r)
      end
    end
  end
end

max_stars_repo_name: ShisatoYano/JuliaAutonomy

max_issues_repo_name: ShisatoYano/JuliaAutonomy

max_forks_repo_name: ShisatoYano/JuliaAutonomy


text:
@with_kw mutable struct RMSpropParameters{T <: AbstractFloat} <: AbstractPolicyParameters
    ρ::T = 0.9
    ϵ::T = 1e-6
end

struct RMSprop <: AbstractSmoothing
    params::RMSpropParameters{Float64}
    grms::Vector{Float64}

    function (::Type{RMSprop})(; kwargs...)
        return new(RMSpropParameters(; kwargs...),Vector{Float64}())
    end
end

function initialize!(rmsprop::RMSprop,x₀::Vector{Float64})
    resize!(rmsprop.grms,length(x₀))
    rmsprop.grms .= zeros(length(x₀))
end

function smooth!(rmsprop::RMSprop,klocal::Integer,kglobal::Integer,x::AbstractVector,gprev::AbstractVector,gcurr::AbstractVector)
    @unpack ρ,ϵ = rmsprop.params
    rmsprop.grms .= ρ*rmsprop.grms + (1-ρ)*(gprev .* gprev)
    gcurr .= gprev ./ (.√rmsprop.grms .+ ϵ)
end


text:
using Conda

Conda.add_channel("conda-forge")
Conda.add("cdsapi")

max_stars_repo_name: henrikhodel/GlobalEnergyGIS

max_issues_repo_name: henrikhodel/GlobalEnergyGIS

max_forks_repo_name: henrikhodel/GlobalEnergyGIS


text:
function grav_potential_cartesian(u)
    x,y,z = u
    return -1.0/sqrt(x^2 + y^2 + z^2)
end

function grav_potential_cylindrical(u)
    r,θ,z = u
    return -1.0/sqrt(r^2 + z^2)
end

function grav_potential_spherical(u)
    r,θ,ϕ = u
    return -1.0/r
end

function grav_field(potential::Function, u::Coordinate)
    return -grad(potential, u)
end

# Scalar Field Functions
function scalar_field(u::Coordinate)
    return u[3]*u[1]^2 + sin(u[1])*u[2]*u[3]^2
end

function grad_scalar_field_cartesian(u::Coordinate)
    x,y,z = u
    return [2*x*z + cos(x)*y*z^2, sin(x)*z^2, x^2 + 2*sin(x)*y*z]
end

function laplacian_scalar_field_cartesian(u::Coordinate)
    x,y,z = u
    return (2*z - sin(x)*y*z^2) + (0.0) + (2*sin(x)*y)
end

function grad_scalar_field_cylindrical(u::Coordinate)
    r,θ,z = u
    return [2*r*z + cos(r)*θ*z^2, (sin(r)*z^2)/r, r^2 + 2*sin(r)*θ*z]
end

function laplacian_scalar_field_cylindrical(u::Coordinate)
    r,θ,z = u
    return (4*z + (θ*z^2)*(cos(r)/r - sin(r))) + 2*sin(r)*θ
end

function grad_scalar_field_spherical(u::Coordinate)
    r,θ,ϕ = u
    return [2*r*ϕ + cos(r)*θ*ϕ^2, (sin(r)*ϕ^2)/(r*sin(ϕ)), (r^2 + 2*sin(r)*θ*ϕ)/r]
end

function laplacian_scalar_field_spherical(u::Coordinate)
    r,θ,ϕ = u
    return (6*ϕ + (θ*ϕ^2)*(2*cos(r)/r - sin(r))) + (cot(ϕ) + (2*θ*sin(r)/(r^2))*(cot(ϕ)*ϕ + 1))
end

#Vector Field Functions
function vector_field(u)
    V = [u[2]*u[1]^2, u[1]*u[2]^2, u[1]*u[2]*u[3]]
    return ContravariantVector(V, u, unit_basis=true)
end

function div_vector_field_cartesian(u)
    x,y,z = u
    return 5*x*y
end

function curl_vector_field_cartesian(u)
    x,y,z = u
    return [x*z, -y*z, y^2 - x^2]
end

function vector_laplacian_cartesian(u)
    x,y,z = u
    return [2*y , 2*x, 0.0]
end

function div_vector_field_cylindrical(u)
    r,θ,z = u
    return 4*θ*r + 2*θ
end

function curl_vector_field_cylindrical(u)
    r,θ,z = u
    return [z, -θ*z, 2*θ^2 - r]
end

function vector_laplacian_cylindrical(u)
    r,θ,z = u
    return [ θ*(3 - 4/r), 2*(1 + 1/r), θ*z/r]
end

function div_vector_field_spherical(u)
    r,θ,ϕ = u
    return 4*θ*r + (2*θ)/sin(ϕ) + θ*(1 + ϕ*cot(ϕ))
end

function curl_vector_field_spherical(u)
    r,θ,ϕ = u
    return [cot(ϕ)*θ^2 - ϕ/sin(ϕ), 2*θ*ϕ, (r/sin(ϕ)) - 2*θ^2]
end

function vector_laplacian_spherical(u)
    r,θ,ϕ = u
    return [4*θ - (2*θ/r)*(1 + 2/sin(ϕ) + ϕ*cot(ϕ)),
            (2*θ^2)/r + (2-θ^2)/(r*sin(ϕ)^2) + 2*(1+ ϕ*cot(ϕ)/r)/sin(ϕ),
            (θ*ϕ/r)*(2 - 1/sin(ϕ)^2) + (θ*cot(ϕ)/r)*(1 - 4/sin(ϕ))]
end

function Φ1(c)
    tan(prod(c))
end

function Φ2(c)
    c[1]*sin(c[2]/c[3])^2
end

function differential_identities(c)

    grad_fg1 = grad(x->Φ1(x)*Φ2(x), c)
    grad_fg2 = Φ1(c)*grad(Φ2,c) + grad(Φ1,c)*Φ2(c)
    @test grad_fg1 ≈ grad_fg2

    div_fA1 = div(x->Φ1(x)*vector_field(x), c)
    div_fA2 = Φ1(c)*div(vector_field,c) + dot(vector_field(c),grad(Φ1,c))
    @test div_fA1 ≈ div_fA2

    curl_fA1 = curl(x->Φ1(x)*vector_field(x), c)
    curl_fA2 = Φ1(c)*curl(vector_field,c) + cross(grad(Φ1,c),vector_field(c))
    @test curl_fA1 ≈ curl_fA2

end

function random_cartesian()
    CartesianCoordinate(3*randn(3))
end

function random_cylindrical()
    x = [0.0, 0.0, -3.0] + [3.0,2pi,6.0].*rand(3)
    CylindricalCoordinate(x)
end

function random_spherical()
    x = [0.0, 0.0, 0.0] + [3.0,2pi,pi].*rand(3)
    SphericalCoordinate(x)
end

const ubc = unit_basis_components

@testset "Differential Operations" begin
    x = CartesianCoordinate([2.0,-1.0,3.0])
    cyl = CylindricalCoordinate([sqrt(x[1]^2 + x[2]^2), atan(x[2],x[1]), x[3]])
    sph = SphericalCoordinate([norm(x), atan(x[2],x[1]), acos(x[3]/norm(x))])

    f = norm(grav_field(grav_potential_cartesian, x))
    @test f ≈ norm(grav_field(grav_potential_cylindrical, cyl))
    @test f ≈ norm(grav_field(grav_potential_spherical, sph))

    x = random_cartesian()
    cyl = random_cylindrical()
    sph = random_spherical()
    println("Testing Gradient...")
    @testset "Gradient" begin
        @test ubc(grad(scalar_field, x)) ≈ grad_scalar_field_cartesian(x)
        @test ubc(grad(scalar_field, cyl)) ≈ grad_scalar_field_cylindrical(cyl)
        @test ubc(grad(scalar_field, sph)) ≈ grad_scalar_field_spherical(sph)
    end

    println("Testing Laplacian...")
    @testset "Laplacian" begin
        @test laplacian(scalar_field, x) ≈ laplacian_scalar_field_cartesian(x)
        @test laplacian(scalar_field, cyl) ≈ laplacian_scalar_field_cylindrical(cyl)
        @test laplacian(scalar_field, sph) ≈ laplacian_scalar_field_spherical(sph)
    end

    println("Testing Divergence...")
    @testset "Divergence" begin
        @test div(vector_field, x) ≈ div_vector_field_cartesian(x)
        @test div(vector_field, cyl) ≈ div_vector_field_cylindrical(cyl)
        @test div(vector_field, sph) ≈ div_vector_field_spherical(sph)
    end

    println("Testing Curl...")
    @testset "Curl" begin
        @test ubc(curl(vector_field, x)) ≈ curl_vector_field_cartesian(x)
        @test ubc(curl(vector_field, cyl)) ≈ curl_vector_field_cylindrical(cyl)
        @test ubc(curl(vector_field, sph)) ≈ curl_vector_field_spherical(sph)
    end

    println("Testing Vector Laplacian...")
    @testset "Vector Laplacian" begin
        @test ubc(vector_laplacian(vector_field, x)) ≈ vector_laplacian_cartesian(x)
        @test ubc(vector_laplacian(vector_field, cyl)) ≈ vector_laplacian_cylindrical(cyl)
        @test ubc(vector_laplacian(vector_field, sph)) ≈ vector_laplacian_spherical(sph)
    end

    println("Testing Differential Identities...")
    @testset "Differential Identities" begin
        @testset "Cartesian" begin
            differential_identities(x)
        end
        @testset "Cylindrical" begin
            differential_identities(cyl)
        end
        @testset "Spherical" begin
            differential_identities(sph)
        end
    end
end


text:
function _log_path(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame})
    joinpath(_obs_path_local(mission, obs_row; kind=:jaxtam), "JAXTAM/obs_log.jld2")
end

function _log_entry(; kwargs...)
    template = OrderedDict(
        :category     => Union{String,Symbol},
        :e_range      => Union{String,Missing},
        :bin_time     => Union{String,Missing},
        :bin_size_sec => Union{Real,Missing},
        :kind         => Symbol,
        :kind_order   => Union{Int,Missing},
        :value        => Union{Any,Missing},
        :group        => Union{Int,Missing},
        :file_name    => Union{String,Missing},
        :path         => Union{String,Missing},
        :notes        => Union{String,Missing},
    )

    kind_order = Dict{Symbol,Int}(
        :lcurve => 1,
        :pgram  => 2,
        :fspec  => 3,
        :sgram  => 4,
        :pulsec => 5,
        :espec  => 6,
        :DELETE => 0,
        :feather_cl => -1,
        :gtis   => -1
    )

    t = Array{Type,1}(collect(values(template)))
    n = collect(keys(template))

    log_entry = DataFrame(t, n, 1)
    for col in names(log_entry)
        if haskey(kwargs, col)
            if col == :e_range
                v = "$(kwargs[:e_range][1])_$(kwargs[:e_range][2])"
            elseif col == :bin_time
                v = "2pow$(Int(log2(kwargs[:bin_time])))"
            elseif col == :file_name
                if haskey(kwargs, :group)
                    v = string("groups/", kwargs[:group], "_", kwargs[:file_name])
                else
                    v = kwargs[:file_name]
                end
            else
                v = kwargs[col]
            end

            log_entry[col] = v
        elseif col == :kind_order
            log_entry[:kind_order] = kind_order[kwargs[:kind]]
        else
            log_entry[col] = missing
        end
    end

    if !haskey(kwargs, :path)
        path = "./"
        [path=joinpath(path, string(v)) for (k, v) in zip(names(log_entry), convert(Matrix, log_entry)) if !in(k, [:kind_order, :value, :group, :path, :notes]) && !ismissing(v)]
    end

    if !haskey(kwargs, :file_name)
        path = string(path, "/") # No file name? Make path a directory path
    end

    log_entry[:path] = path

    return log_entry
end

function _log_query_path(; kwargs...)
    return _log_entry(; kwargs...)[1, :path]
end

function _log_gen(path::String, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame})
    if !ispath(path)
        mkpath(dirname(path))
    end

    meta_info = Dict{Any,Any}(
        :obs_row    => obs_row,
        :downloaded => false
    )

    save(path, Dict{Any,Any}("meta"=>meta_info))
end

function _log_read(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame})
    path = _log_path(mission, obs_row)

    if !isfile(path)
        @info "Generating log for $(obs_row[:obsid])"
        _log_gen(path, obs_row)
    end

    Dict{Any,Any}(load(path))
end

function _merge(log::Dict, k::Union{Symbol,Tuple}, entry::Dict{T,DataFrames.DataFrame}) where T <: Any
    template_log = _log_entry(; category="DELETE", kind=:DELETE)
    deleterows!(template_log, 1)

    if haskey(log, k)
        append!(template_log, log[k])
    end

    append!(template_log, entry[k])
    unique!(template_log)

    log[k] = template_log

    return log
end

function _merge(log::Dict, k::Union{Symbol,Tuple}, entry::Dict)
    log   = convert(Dict{Any,Any}, log)
    entry = convert(Dict{Any,Any}, entry)

    merge!(log, entry)

    return log
end

function _log_add_recursive(log::Dict, entry::Dict)
    log = Dict{Any,Any}(log)
    for (k, v) in entry
        if typeof(v) <: Dict
            if haskey(log, k)
                log[k] = _log_add_recursive(log[k], entry[k])
            else
                log[k] = entry[k]
                return log
            end
        else
            return _merge(log, k, entry)
        end
    end

    return log
end

function _log_add(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}, entry::Dict{String,T}) where T <: Any
    path = _log_path(mission, obs_row)
    log  = _log_read(mission, obs_row)

    log_new = _log_add_recursive(log, entry)
    save(path, log_new)
    return log_new
end

function _log_query(mission::Mission, obs_row::DataFrames.DataFrameRow{DataFrames.DataFrame}, args...; surpress_warn=true)
    log = _log_read(mission, obs_row)

    reply = log
    for (i, key) in enumerate(args)
        if reply isa Dict
            if haskey(reply, key)
                reply = reply[key]
            else
                if !surpress_warn
                    @warn "Key '$key' not found in $(obs_row[:obsid]) log, available: $(keys(reply))"
                end
                return missing
            end
        else
            @warn "Entry not Dict, gone too deep? Accessing: $(args[1:i]) gives log at '$key' as $(typeof(reply)) not Dict"
            return missing
        end
    end

    return reply
end


text:
# Copyright (C) 2021 Ilya Baran.  This program is distributed under the terms of the MIT license.

using ModelBlocks
using Test

println("Variables Test starting");

variableArray = [
    Variable("foo", 3, (2, 4), "meters", "Basic foo variable"),
    Variable("bar", 4, (2, 5), "meters", "Basic bar variable"),
    Variable("m_3", 2, (1, 4), "moles", "This is m_3"),
];
variables = Variables(variableArray);

@test variables.foo == 3;
@test variables.bar == 4;
@test variables.m_3 == 2;
@test variables.values == [3, 4, 2];
@test variables.nameToIndex["foo"] == 1;
@test variables["foo"] == 3;

copyto!(variables.values, [4, 2, 1]);
@test variables.values == [4, 2, 1];

@test_throws ErrorException Variable("3foo", 3, (2, 4), "meters", "Invalid identifier")
@test_throws ErrorException Variable("foo", 3, (5, 4), "meters", "Invalid range")

resize!(variableArray, 4);
variableArray[4] = Variable("foo", 3, (2, 4), "meters", "Duplicate foo variable");
@test_throws ErrorException Variables(variableArray);

# Test combineVariables

v1 = Variables([
    Variable("a", 1, (1, 4), "", "A"),
    Variable("b", 2, (1, 5), "", "B"),
    Variable("c", 3, (1, 6), "", "C"),
]);

v2 = Variables([
    Variable("b", 5, (1, 16), "", "B"),
    Variable("d", 7, (1, 14), "", "D"),
    Variable("c", 6, (1, 15), "", "C"),
]);

vUnion = variablesUnion(v1, v2);
@test vUnion.values == [1, 5, 6, 7];

println("Variables Test finished");

max_stars_repo_name: ilya-baran-personal/modelBlocks

max_issues_repo_name: ilya-baran-personal/modelBlocks

max_forks_repo_name: ilya-baran-personal/modelBlocks


text:
# simple test that things get created
@testset "RSS gen" begin
    f = joinpath(p, "basic", "__site", "feed.xml")
    @test isfile(f)
    fc = prod(readlines(f, keep=true))

    @test occursin(raw"""<rss version="2.0" xmlns:atom="http://www.w3.org/2005/Atom">""", fc)
    @test occursin(raw"""<title><![CDATA[Franklin Template]]></title>""", fc)
    @test occursin(raw"""<description><![CDATA[Example website using Franklin
  ]]></description>""", fc)
    @test !occursin(raw"""<author>""", fc)
    @test occursin(raw"""<link>https://tlienart.github.io/FranklinTemplates.jl/menu1/</link>""", fc)
    @test occursin(raw"""<description><![CDATA[A short description of the page which would serve as blurb in a RSS feed;""", fc)
end

@testset "Sitemap gen" begin
    f = joinpath(p, "basic", "__site", "sitemap.xml")
    @test isfile(f)
    fc = prod(readlines(f, keep=true))

    @test occursin(raw"""
        <?xml version="1.0" encoding="utf-8" standalone="yes" ?>
        <urlset xmlns="http://www.sitemaps.org/schemas/sitemap/0.9">""", fc)
    # check pages
    for pg in ("", "menu1", "menu2", "menu3")
        pgg = joinpath(pg, "index.html")
        @test occursin("""
            <loc>https://tlienart.github.io/FranklinTemplates.jl/$pgg</loc>""", fc)
    end
end

@testset "Robots.txt gen" begin
    f = joinpath(p, "basic", "__site", "robots.txt")
    @test isfile(f)
    fc = prod(readlines(f, keep=true))

    @test occursin(raw"""
        Sitemap: https://tlienart.github.io/FranklinTemplates.jl/sitemap.xml""", fc)
end


text:
# "lines and linestyles"
begin
    s = Scene(resolution = (800, 800), camera = campixel!)
    scalar = 30
    points = Point2f[(1, 1), (1, 2), (2, 3), (2, 1)]
    linestyles = [
        :solid, :dash, :dot, :dashdot, :dashdotdot,
        [1, 2, 3], [1, 2, 4, 5]
    ]
    for linewidth in 1:10
        for (i, linestyle) in enumerate(linestyles)
            lines!(s,
                scalar .* (points .+ Point2f(linewidth*2, i * 3.25)),
                linewidth = linewidth,
                linestyle = linestyle,
            )
        end
    end
    s
end

# "lines with gaps"
begin
    s = Scene(resolution = (800, 800), camera = campixel!)
    points = [
        Point2f[(1, 0), (2, 0.5), (NaN, NaN), (4, 0.5), (5, 0)],
        Point2f[(NaN, NaN), (2, 0.5), (3, 0), (4, 0.5), (5, 0)],
        Point2f[(NaN, NaN), (2, 0.5), (3, 0), (4, 0.5), (NaN, NaN)],
        Point2f[(NaN, NaN), (2, 0.5), (NaN, NaN), (4, 0.5), (NaN, NaN)],
        Point2f[(NaN, NaN), (NaN, NaN), (3, 0), (NaN, NaN), (NaN, NaN)],
        Point2f[(NaN, NaN), (NaN, NaN), (NaN, NaN), (NaN, NaN), (NaN, NaN)],
    ]
    for (i, p) in enumerate(points)
        lines!(s, (p .+ Point2f(0, i)) .* 100, linewidth = 10)
    end
    s
end

# "scatters"
begin
    s = Scene(resolution = (800, 800), camera = campixel!)

    markersizes = 0:2:30
    markers = [:circle, :rect, :cross, :utriangle, :dtriangle,
        'a', 'x', 'h', 'g', 'Y', 'J', 'α', '↑'
    ]

    for (i, ms) in enumerate(markersizes)
        for (j, m) in enumerate(markers)
            scatter!(s,
                Point2f(i, j) .* 45,
                marker = m,
                markersize = ms,
            )
        end
    end
    s
end

# "scatter rotations"
begin
    s = Scene(resolution = (800, 800), camera = campixel!)

    rotations = range(0, 2pi, length = 15)
    markers = [:circle, :rect, :cross, :utriangle, :dtriangle,
        'a', 'x', 'h', 'g', 'Y', 'J', 'α', '↑'
    ]

    for (i, rot) in enumerate(rotations)
        for (j, m) in enumerate(markers)
            p = Point2f(i, j) .* 45
            scatter!(s,
                p,
                marker = m,
                markersize = 30,
                rotations = rot,
            )
            scatter!(s, p, color = :red, markersize = 6)
        end
    end
    s
end

# "poly"
begin
    s = Scene(resolution = (800, 800), camera = campixel!)

    scalefactor = 70
    Pol = Makie.GeometryBasics.Polygon

    translate(pol, p) = Pol(
        Ref(p) .+ decompose(Point2f, pol.exterior),
        map(x -> decompose(Point2f, x) .+ Ref(p), pol.interiors)
    )
    scale(pol, sca) = Pol(
        sca .* decompose(Point2f, pol.exterior),
        map(x -> decompose(Point2f, x) .* sca, pol.interiors)
    )

    polys = [
        # three points
        Pol(Point2f[(1, 1), (1, 2), (2, 1)]),
        # four points
        Pol(Point2f[(1, 1), (1, 2), (2, 2), (2, 1)]),
        # double point
        Pol(Point2f[(1, 1), (1, 2), (2, 2), (2, 2)]),
        # one hole
        Pol(
            Point2f[(1, 1), (1, 2), (2, 2), (2, 1)],
            [Point2f[(1.3, 1.3), (1.3, 1.7), (1.7, 1.7), (1.7, 1.3)]]
        ),
        # two holes
        Pol(
            Point2f[(1, 1), (1, 2), (2, 2), (2, 1)],
            [
                Point2f[(1.15, 1.15), (1.15, 1.85), (1.4, 1.85), (1.4, 1.15)],
                Point2f[(1.6, 1.15), (1.6, 1.85), (1.85, 1.85), (1.85, 1.15)],
            ]
        ),
        # hole half same as exterior
        Pol(
            Point2f[(1, 1), (1, 2), (2, 2), (2, 1)],
            [Point2f[(1, 1), (1, 2), (2, 2)]],
        ),
        # point self intersection
        Pol(
            Point2f[(1, 1), (2, 1), (2, 2), (1.5, 1), (1, 2)],
        ),
    ]

    linewidths = 0:2:9

    for (i, p) in enumerate(polys)
        for (j, lw) in enumerate(linewidths)
            poly!(
                s,
                translate(
                    scale(p, scalefactor),
                    scalefactor * Point2f(1.3 * (i-1), 1.3 * j)
                ),
                color = (:red, 0.5),
                strokewidth = lw,
            )
        end
    end

    s
end


text:
using StokesParticles, Test
import LinearAlgebra: norm

d = 1                                   # dimension of data
M = 10_000                              # number of targets
N = 10_000                              # number of sources
h = 0.2                                 # bandwidth
ε = 1e-2                                # max tolerance
x = rand(N*d)                           # source vectors, reshape'd to 1D array
y = collect(range(-2,stop=2,length=M))  # target vectors, reshape'd to 1D array
W = 1                                   # number of sources being evaluated
q = rand(N*W)                           # weights for sources
v1 = rand(M)                            # output values array
v2 = rand(M)                            # output 2
v3 = rand(M)                            # output 3

t1 = @elapsed StokesParticles.fgt!(v1,d,M,N,h,ε,x,y,q,W)
t2 = @elapsed StokesParticles.dgt!(v2,d,M,N,h,ε,x,y,q,W)
t3 = @elapsed StokesParticles.mydgt!(v3,d,M,N,h,ε,x,y,q,W)

@test t1 < t2 < t3
@test norm(v1 - v3)/norm(v3) < 1e-2     # approx alg, some error allowed
@test norm(v2 - v3)/norm(v3) < eps()    # same alg, so no error allowed


text:
using ShipstabNURBS
using Test

@testset "FindSpan" begin
    U = [0 0 0 0 2/5 3/5 3/5 1 1 1 1];
    t1,s1 = FindSpan(0,U)
    t2,s2 = FindSpan(2/5,U)
    t3,s3 = FindSpan(1,U)
  @test t1 == 4
  @test s1 == 3
  @test t2 == 7
  @test s2 == 2
  @test t3 == 7
  @test s3 == 0
end

@testset "BasisFuns" begin
    P = [0 1 3 5 6 8 9;0 2 4 2 0 0 3];
    U = [0 0 0 0 2/5 3/5 3/5 1 1 1 1];
    n = size(P,2) - 1;
    m = length(U) - 1;
    p = m-n-1;

    u = range(0,1,length=100);
    funs = zeros(length(u),n+1);

    for i = 1:length(u)
        funs[i,:] = BasisFuns(n,p,u[i],U);
    end
    @test sum(funs) == 100.0
end

max_stars_repo_name: barche/ShipstabNURBS

max_issues_repo_name: barche/ShipstabNURBS

max_forks_repo_name: barche/ShipstabNURBS


text:
# Test typing in a typed gripper domain
@testset "typing" begin

path = joinpath(dirname(pathof(PDDL)), "..", "test", "typing")

domain = load_domain(joinpath(path, "domain.pddl"))
@test domain.name == Symbol("gripper-typed")
@test convert(Term, domain.predicates[:free]) == pddl"(free ?g)"
@test domain.predicates[:carry].argtypes == (:ball, :gripper)
@test :gripper in PDDL.get_types(domain)

problem = load_problem(joinpath(path, "problem.pddl"))
@test problem.name == Symbol("gripper-problem")
@test problem.objects == @pddl("rooma", "roomb", "ball1", "ball2", "left", "right")
@test problem.objtypes[Const(:ball1)] == :ball

state = initstate(domain, problem)
implementations = [
    "concrete interpreter" => (domain, state),
    "abstract interpreter" => abstracted(domain, state),
    "ground interpreter" => (ground(domain, state), state),
    "concrete compiler" => compiled(domain, state),
    "abstract compiler" => compiled(abstracted(domain), state)
]

@testset "typing ($name)" for (name, (domain, _)) in implementations

    # Test forward execution of plans
    state = initstate(domain, problem)
    state = execute(domain, state, pddl"(pick ball1 rooma left)")
    @test domain[state => pddl"(carry ball1 left)"] ≃ true
    state = execute(domain, state, pddl"(move rooma roomb)")
    @test domain[state => pddl"(robbyat roomb)"] ≃ true
    state = execute(domain, state, pddl"(drop ball1 roomb left)")
    @test domain[state => pddl"(at ball1 roomb)"] ≃ true

    @test satisfy(domain, state, problem.goal) ≃ true

    # Test action availability
    state = initstate(domain, problem)
    @test Set{Term}(available(domain, state)) == Set{Term}(@pddl(
        "(pick ball1 rooma right)", "(pick ball1 rooma left)",
        "(pick ball2 rooma right)", "(pick ball2 rooma left)",
        "(move rooma roomb)", "(move rooma rooma)"
    ))

end

# Test backward regression of plans
state = goalstate(domain, problem)
state = regress(domain, state, pddl"(drop ball1 roomb left)")
@test domain[state => pddl"(carry ball1 left)"] == true
state = regress(domain, state, pddl"(move rooma roomb)")
@test domain[state => pddl"(robbyat rooma)"] == true
state = regress(domain, state, pddl"(pick ball1 rooma left)")
@test domain[state => pddl"(at ball1 rooma)"] == true
@test issubset(state, initstate(domain, problem))

# Test action relevance
state = goalstate(domain, problem)
@test Set{Term}(relevant(domain, state)) == Set{Term}(@pddl(
    "(drop ball1 roomb left)", "(drop ball1 roomb right)"
))

end # typing


text:
type GbGlb
    # REMEMBER -- all multi-dim arrays are column-indexing!  [i,j] => [j*N+i]
    particles::Array{Float64,1} # [Ndim x Ndens]      // means of selected particles
    variance::Array{Float64,1}  # [Ndim x Ndens]      //   variance of selected particles
    p::Array{Float64,1}         # [Np]                // probability of ith kernel
    ind::Array{Int64,1}                    # current indexes of MCMC step
    Malmost::Array{Float64,1}
    Calmost::Array{Float64,1}   #[Ndim x 1]   // Means & Cov. of all indices but j^th
    # Random number callback
    randU::Array{Float64,1}     # [Npoints * Ndens * (Niter+1)] uniformly distrib'd random variables
    randN::Array{Float64,1}     # [Ndim * Npoints] normally distrib'd random variables
    Ndim::Int64
    Ndens::Int64
    Nlevels::Int64
    dNp::Int64
    newPoints::Array{Float64,1}
    newWeights::Array{Float64,1}
    newIndices::Array{Int64,1}
    trees::Array{BallTreeDensity,1}
    levelList::Array{Int64,2}
    levelListNew::Array{Int64,2}
    dNpts::Array{Int64,1}
    ruptr::Int64
    rnptr::Int64
end

function makeEmptyGbGlb()
   return GbGlb(zeros(0),
                zeros(0),
                zeros(0),
                zeros(Int64,0),
                zeros(0),
                zeros(0),
                zeros(0),
                zeros(0),
                0,0,0,0,
                zeros(0),
                zeros(0),
                ones(Int64,0),
                Array(BallTreeDensity,1),
                ones(Int64,1,0),
                ones(Int64,1,0),
                zeros(Int64,0),
                0, 0)
end

type MSCompOpt
  pT::Float64
  tmpC::Float64
  tmpM::Float64
end

function levelInit!(glb::GbGlb)
  for j in 1:glb.Ndens
    glb.dNpts[j] = 1
    glb.levelList[j,1] = root(glb.trees[j])
  end
end

# this is a weird function -- not tested at all at this point
function initIndices!(glb::GbGlb)
  count=1
  for j in 1:glb.Ndens
    dNp = glb.dNpts[j]
    zz=glb.levelList[j,1]
    z = 1
    while z <= dNp # ++z, used to be a for loop
      glb.p[z] = weight(glb.trees[j], zz);  # init by sampling from weights
        z+=1
        if z<=dNp zz=glb.levelList[j,z]; end
    end

    for z in 2:dNp
        glb.p[z] += glb.p[z-1];
    end
    zz=glb.levelList[j,1]
    z = 1
    while z <= (dNp-1)
      if (glb.randU[glb.ruptr] <= glb.p[z]) break; end #count
        z+=1
        if z<dNp zz=glb.levelList[j,z]; end
    end
    glb.ind[j] = zz;
    count+=1
    glb.ruptr += 1
  end
  #glb.randU = glb.randU[count:end]
end

function calcIndices!(glb::GbGlb)
  #println("oooooooooooooooooooooooooooooooo")
  @fastmath @inbounds begin
    for j in 1:glb.Ndens
      for z in 1:glb.Ndim
        #@show j, z, glb.ind[j]
        #@show mean(glb.trees[j],glb.ind[j])
        glb.particles[z+glb.Ndim*(j-1)] = mean(glb.trees[j],glb.ind[j], z)#[z];
        glb.variance[z+glb.Ndim*(j-1)]  = bw(glb.trees[j],glb.ind[j], z);
      end
    end
  end
  #println("pppppppppppppppppppppppppppppppp")
end

function samplePoint!(X::Array{Float64,1}, glb::GbGlb, frm::Int64)
  #counter = 1
  for j in 1:glb.Ndim
    mn=0.0; vn=0.0;
    for z in 1:glb.Ndens                             # Compute mean and variances of
      vn += 1.0/glb.variance[j+glb.Ndim*(z-1)]       #   product of selected particles
      #@show j, z, j+glb.Ndim*(z-1)
      #@show round([mn, glb.particles[j+glb.Ndim*(z-1)], glb.variance[j+glb.Ndim*(z-1)]],3)
      mn += glb.particles[j+glb.Ndim*(z-1)]/glb.variance[j+glb.Ndim*(z-1)]
    end
    vn = 1.0/vn; mn *= vn;
    glb.rnptr += 1
    X[j+frm] = mn + sqrt(vn) * glb.randN[glb.rnptr] #counter       # then draw a sample from it
    #counter+=1
  end
  #glb.randN = glb.randN[counter:end]   # TODO ugh
  return Union{}
end

function levelDown!(glb::GbGlb)
  for j in 1:glb.Ndens
    z = 1
    for y in 1:(glb.dNpts[j])
      #@show "left", y, j, left(glb.trees[j], glb.levelList[j,y])
      if validIndex(glb.trees[j], left(glb.trees[j], glb.levelList[j,y]))
        #@show "a", y, j, z, left(glb.trees[j], glb.levelList[j,y])
        glb.levelListNew[j,z] = left(glb.trees[j], glb.levelList[j,y])
        z+=1
      end
      #@show "right", y, j, right(glb.trees[j], glb.levelList[j,y])
      if validIndex(glb.trees[j], right(glb.trees[j], glb.levelList[j,y]))
        #@show "b", y, j, z, right(glb.trees[j], glb.levelList[j,y])
        glb.levelListNew[j,z] = right(glb.trees[j], glb.levelList[j,y]);
        z+=1
      end
      #@show "last", glb.ind[j], glb.levelList[j,y]
      if (glb.ind[j] == glb.levelList[j,y])                      # make sure ind points to
        #@show "c", y, j, z, glb.levelListNew[j,z-1]
        glb.ind[j] = glb.levelListNew[j,z-1]                     #  a child of the old ind
      end
    end
    glb.dNpts[j] = z-1
    #println("levelDown! -- glb.dNpts[j=$j]=$(glb.dNpts[j])")
  end
  tmp = glb.levelList                            # make new list the current
  glb.levelList = glb.levelListNew              #   list and recycle the old
  glb.levelListNew=tmp

  #println("A, levelList $(glb.levelList[1,:]), $(glb.levelList[2,:])")
  #println("A, levelListNew $(glb.levelListNew[1,:]), $(glb.levelListNew[2,:])")
  #println("A, ind $(glb.ind)")
end

function sampleIndices!(X::Array{Float64,1}, cmoi::MSCompOpt, glb::GbGlb, frm::Int64)#pT::Array{Float64,1}
  counter=1
  zz=0
  #z = 1
  for j in 1:glb.Ndens
    dNp = glb.dNpts[j]    #trees[j].Npts();
    cmoi.pT = 0.0

    ##z = 1
    zz=glb.levelList[j,1]#z
    #println("start zz=glb.levelList[$(j),$(z)]=$(zz)")
    for z in 1:dNp ##while <=
      #println("starting first z=$(z), dNp=$(dNp), zz=$(zz)")
      glb.p[z] = 0.0
      for i in 1:glb.Ndim
        tmp = X[i+frm] - mean(glb.trees[j], zz, i)#[i]
        glb.p[z] += (tmp*tmp) / bw(glb.trees[j], zz, i)#[i]
        glb.p[z] += Base.Math.JuliaLibm.log(bw(glb.trees[j], zz, i))
      end
      glb.p[z] = exp( -0.5 * glb.p[z] ) * weight(glb.trees[j], zz)
      cmoi.pT += glb.p[z]
      #println("end ++z zz=glb.levelList[$(j),$(z)]=$(zz)")
      ##z+=1
      #println("1 ++z=$(z)")
      zz = z<dNp ? glb.levelList[j,z+1] : zz
    end
    #println("after first zz=$(zz)")
    @simd for z in 1:dNp
        glb.p[z] /= cmoi.pT
    end
    @simd for z in 2:dNp
        glb.p[z] += glb.p[z-1]              # construct CDF and sample a
    end

    z=1
    zz=glb.levelList[j,z]
    #println("mid zz=glb.levelList[$(j),$(z)]=$(zz)")
    while z<=(dNp-1)
      #println("second zz=glb.levelList[$(j),$(z)]=$(zz)")
      if (glb.randU[glb.ruptr] <= glb.p[z]) # counter
        break;                              # new kernel from jth density
      end
      z+=1
      #println("2 ++z=$(z)")
      if z<=dNp
        zz=glb.levelList[j,z]
      else
        error("This should never happend due to -1 MSGibbs01.jl")
      end
    end
    glb.ind[j] = zz                          #  using those weights
    #println("sampleIndices -- glb.ind[j=$(j)]=$(glb.ind[j])")
    counter+=1
    glb.ruptr += 1
  end
  #glb.randU = glb.randU[counter:end]
  calcIndices!(glb);                         # recompute particles, variance
end

## SLOWEST PIECE OF THE COMPUTATION -- TODO
# easy PARALLELs overhead here is much slower, already tried -- rather search for BLAS optimizations...
function makeFasterSampleIndex!(j::Int64, cmo::MSCompOpt, glb::GbGlb)
  #pT::Array{Float64,1}
  cmo.tmpC = 0.0
  cmo.tmpM = 0.0

  zz=glb.levelList[j,1]
  #z=1
  for z in 1:(glb.dNpts[j])#dNp
    glb.p[z] = 0.0
    for i in 1:glb.Ndim
      cmo.tmpC = bw(glb.trees[j], zz, i) + glb.Calmost[i]
      cmo.tmpM = mean(glb.trees[j], zz, i) - glb.Malmost[i]
      glb.p[z] += abs2(cmo.tmpM)/cmo.tmpC + log(cmo.tmpC) # This is the slowest piece
    end
    glb.p[z] = exp( -0.5 * glb.p[z] ) * weight(glb.trees[j].bt, zz) # slowest piece
    z < glb.dNpts[j] ? zz = glb.levelList[j,(z+1)] : nothing
  end

  # incorrect! remember zz
  #@simd for z in 1:glb.dNpts[j]
  #  glb.p[z] = exp( -0.5 * glb.p[z] ) * weight(glb.trees[j], zz)
  #end
  nothing
end

function sampleIndex(j::Int64, cmo::MSCompOpt, glb::GbGlb)
#pT::Array{Float64,1}
  #dNp = glb.dNpts[j];  #trees[j].Npts();
  cmo.pT = 0.0

  # determine product of selected particles from all but jth density
  for i in 1:glb.Ndim
    iCalmost = 0.0; iMalmost = 0.0;
    for k in 1:glb.Ndens
      if (k!=j) iCalmost += 1.0/glb.variance[i+glb.Ndim*(k-1)]; end
      if (k!=j) iMalmost += glb.particles[i+glb.Ndim*(k-1)]/glb.variance[i+glb.Ndim*(k-1)]; end
    end
    glb.Calmost[i] = 1/iCalmost;
    glb.Malmost[i] = iMalmost * glb.Calmost[i];
  end

  makeFasterSampleIndex!(j, cmo, glb)

  @inbounds @simd for k in 1:glb.dNpts[j]
    cmo.pT += glb.p[k]
  end

  @inbounds @simd for z in 1:glb.dNpts[j]#dNp
    glb.p[z] /= cmo.pT            # normalize weights# normalize weights
  end
  @inbounds @simd for z in 2:glb.dNpts[j]#dNp
    glb.p[z] += glb.p[z-1]    # construct CDF and sample
  end
  zz=glb.levelList[j,1]
  z=1
  while z<=(glb.dNpts[j]-1)#dNp
    if (glb.randU[glb.ruptr] <= glb.p[z]) break;  end   #1  #   a new kernel from the jth
    z+=1
    if z<=glb.dNpts[j]#dNp
      zz=glb.levelList[j,z]
      #println("sampleIndex 2++z, , zz=glb.levelList[$(j),$(z)]=$(zz)")
    end
  end
  glb.ind[j] = zz;                                          #   density using these weights
  #glb.randU = glb.randU[2:end]
  glb.ruptr += 1

  for i in 1:glb.Ndim
    #@show round(glb.trees[j].means, 3)
    #@show i, j, glb.ind[j], mean(glb.trees[j], glb.ind[j])
    glb.particles[i+glb.Ndim*(j-1)] = mean(glb.trees[j], glb.ind[j], i)#[i]
    glb.variance[i+glb.Ndim*(j-1)]  = bw(glb.trees[j], glb.ind[j], i)#[i];
  end
end

function printGlbs(g::GbGlb, tag=Union{})
    if tag==Union{}
        println("=========================================================================")
    else
        println(string(tag,"================================================================"))
    end
    println("Ndim=$(g.Ndim), Ndens=$(g.Ndens), Nlevels=$(g.Nlevels), dNp=$(g.dNp), dNpts=$(g.dNpts)")
    @show g.ind
    @show round(g.particles,2)
    @show round(g.variance,2)
    @show round(g.p,2)
    @show round(g.Malmost,2)
    @show round(g.Calmost,2)
    @show g.levelList
    @show g.levelListNew
    @show round(g.newPoints,4)
    @show g.newIndices
end

function gibbs1(Ndens::Int64, trees::Array{BallTreeDensity,1},
                Np::Int64, Niter::Int64,
                pts::Array{Float64,1}, ind::Array{Int64,1},
                randU::Array{Float64,1}, randN::Array{Float64,1})

    glbs = makeEmptyGbGlb()
    glbs.Ndens = Ndens
    glbs.trees = trees
    glbs.newPoints = pts
    glbs.newIndices = ind
    glbs.randU = randU
    glbs.randN = randN
    glbs.Ndim = trees[1].bt.dims

    maxNp = 0                         # largest # of particles we deal with
    for tree in trees
        if (maxNp < Npts(tree))
            maxNp = Npts(tree)
        end
    end

    glbs.ind = ones(Int64,Ndens)
    glbs.p = zeros(maxNp)
    glbs.Malmost = zeros(glbs.Ndim)
    glbs.Calmost = zeros(glbs.Ndim)
    glbs.Nlevels = floor(Int64,((log(maxNp)/log(2))+1))
    glbs.particles = zeros(glbs.Ndim*Ndens)
    glbs.variance  = zeros(glbs.Ndim*Ndens)
    glbs.dNpts = zeros(Int64,Ndens)
    glbs.levelList = ones(Int64,Ndens,maxNp)
    glbs.levelListNew = ones(Int64,Ndens,maxNp)
    cmo = MSCompOpt(0.0, 0.0, 0.0)
    cmoi = MSCompOpt(0.0, 0.0, 0.0)

    ##@show glbs.ruptr, size(glbs.randU)
    for s in 1:Np   #   (for each sample:)
          #println("+++++++++++++++++++++++++++++NEW POINT+++++++++++++++++++++++++++++++++++++++")
        frm = ((s-1)*glbs.Ndim)

        levelInit!(glbs)
        initIndices!(glbs)
        calcIndices!(glbs)

        #printGlbs(glbs, string("s=",s-1.5))

        for l in 1:glbs.Nlevels
          samplePoint!(glbs.newPoints, glbs, frm)
          levelDown!(glbs);
          sampleIndices!(glbs.newPoints, cmoi, glbs, frm);

          #printGlbs(glbs, string("s=",s-1.5+.1*l))

          for i in 1:Niter         #   perform Gibbs sampling
            for j in 1:glbs.Ndens
              @fastmath @inbounds sampleIndex(j, cmo, glbs);
            end
          end
        end

        for j in 1:glbs.Ndens                       #/ save and
          glbs.newIndices[(s-1)*glbs.Ndens+j] = getIndexOf(glbs.trees[j], glbs.ind[j])+1;  # return particle label
        end

        samplePoint!(glbs.newPoints, glbs, frm);                            # draw a sample from that label
               #glbs.newIndices = glbs.newIndices[glbs.Ndens:end]  # move pointers to next sample
               #glbs.newPoints  = glbs.newPoints[glbs.Ndim:end]
        #printGlbs(glbs, string("s=",s-1))
    end
    ##@show glbs.ruptr, size(glbs.randU), size(glbs.randU,1)+glbs.ruptr
    # deref the glbs data and let gc() remove it later, or force with gc()
    glbs = 0
    #@show @elapsed gc()
    #error("gibbs1 -- not implemented yet")
    nothing
end

# function remoteProdAppxMSGibbsS(npd0::BallTreeDensity,
#                           npds::Array{BallTreeDensity,1}, anFcns, anParams,
#                           Niter::Int64=5)
#
#   len = length(npds)
#   arr = Array{Array{Float64,2},1}(len)
#
#   d = npds[1].bt.dims
#   N = npds[1].bt.num_points
#
#
#   for i in 1:len
#     arr[i] = reshape(npds[i].bt.centers[(N*d+1):end],d,N)
#   end
#
#   pts, bw = remoteProd(arr)
#
#   return pts, -1
# end

function prodAppxMSGibbsS(npd0::BallTreeDensity,
                          npds::Array{BallTreeDensity,1}, anFcns, anParams,
                          Niter::Int64=5)
    # See  Ihler,Sudderth,Freeman,&Willsky, "Efficient multiscale sampling from products
    #         of Gaussian mixtures", in Proc. Neural Information Processing Systems 2003

    Ndens = length(npds)              # of densities
    Ndim  = npds[1].bt.dims           # of dimensions
    Np    = Npts(npd0)                # of points to sample

    # skipping analytic functions for now TODO ??

    UseAn = false
    #??pointsM = zeros(Ndim, Np)
    points = zeros(Ndim*Np)
    #??plhs[1] = mxCreateNumericMatrix(Ndens, Np, mxUINT32_CLASS, mxREAL);
    indices=ones(Int64,Ndens*Np)
    maxNp = Np                        # largest # of particles we deal with
    for tree in npds
        if (maxNp < Npts(tree))
            maxNp = Npts(tree)
        end
    end
    Nlevels = floor(Int64,(log(Float64(maxNp))/log(2.0))+1.0)  # how many levels to a balanced binary tree?

    # Generate enough random numbers to get us through the rest of this
    if true
      randU = rand(Int64(Np*Ndens*(Niter+2)*Nlevels))
      randN = randn(Int64(Ndim*Np*(Nlevels+1)))
    else
        randU = vec(readdlm("randU.csv"))
        randN = vec(readdlm("randN.csv"))
    end

    #@show size(randU), size(randN)
    gibbs1(Ndens, npds, Np, Niter, points, indices, randU, randN);

    return reshape(points, Ndim, Np), reshape(indices, Ndens, Np)
end

function *(p1::BallTreeDensity, p2::BallTreeDensity)
  numpts = round(Int,(Npts(p1)+Npts(p2))/2)
  d = Ndim(p1)
  d != Ndim(p2) ? error("kdes must have same dimension") : nothing
  dummy = kde!(rand(d,numpts),[1.0]);
  pGM, = prodAppxMSGibbsS(dummy, [p1;p2], Union{}, Union{}, 5)
  return kde!(pGM)
end


text:
include("../module/CustomStencil.jl")
include("../misc/misc.jl")
include("../misc/comp_count.jl")
include("../misc/cpu_stencils.jl")

## Star Stencil definition with radius = 1
coefs = [1;0.5;]
star_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1, 1, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i] +
            D[x-i,y,z] + D[x,y-i,z] + D[x,y,z-i]))
st_def = CreateStencilDefinition(star_stencil, coefs)
st_single = NewStencilInstance(st_def, m_step=false)
sym = st_single.stencil_sym
countComp(sym,1,7)
## Dense Stencil definition with radius = 1
coefs = [0.5;0.25]
dense_stencil = @def_stencil_expression(
            @sum(i,-1,1,
                @sum(j,-1,1,
                    @sum(k,-1,1, c[max(abs(i),abs(j),abs(k))]*D[x+i,y+j,z+k]))))

st_def = CreateStencilDefinition(dense_stencil, coefs)
st_single = NewStencilInstance(st_def, m_step=false)
sym = st_single.stencil_sym
countComp(sym,1,4)
## Flux summation Stencil definition with radius = 1
coefs = [1;0.2]
flux_stencil = @def_stencil_expression c[0]D[x,y,z] + @sum(i, 1,1, c[i]*(
            D[x+i,y,z] + D[x,y+i,z] + D[x,y,z+i]))
st_def = CreateStencilDefinition(flux_stencil, coefs)
st_single = NewStencilInstance(st_def, m_step=false)
sym = st_single.stencil_sym
a = countComp(sym,1,8)


text:
#= extract_pcd_points.jl
extract the radar-points in pcd file(ascii) =#

using XLSX;

prj_path = dirname(dirname(@__FILE__))
data_path = prj_path * raw"\data\data_no_anno\radar_ascii";

xls_file_path = prj_path * raw"\data\radar_process\radar_points.xlsx";
XLSX.openxlsx(xls_file_path, mode="w") do xf
    for file_id in 1:38
        file_path = data_path * "\\$file_id.pcd"
        f = open(file_path, "r");
        n = countlines(f)

        # skip the first 12 lines
        seekstart(f)
        for i = 1:12
            readline(f)
        end

        data = Array{Float64}(undef, n - 12, 18)
        for i = 13:n
            row_id = i - 12;
            str_list = split(readline(f), " ")

            x = parse(Float64, str_list[1]);
            data[row_id, 1] = x;

            y = parse(Float64, str_list[2]);
            data[row_id, 2] = y;

            z = parse(Float64, str_list[3]);
            data[row_id, 3] = z;

            dyn_prop = parse(Float64, str_list[4]);
            data[row_id, 4] = dyn_prop;

            id = parse(Float64, str_list[5]);
            data[row_id, 5] = id;

            rcs = parse(Float64, str_list[6]);
            data[row_id, 6] = rcs;

            vx = parse(Float64, str_list[7]);
            data[row_id, 7] = vx;

            vy = parse(Float64, str_list[8]);
            data[row_id, 8] = vy;

            vx_comp = parse(Float64, str_list[9]);
            data[row_id, 9] = vx_comp;

            vy_comp = parse(Float64, str_list[10]);
            data[row_id, 10] = vy_comp;

            is_quality_valid = parse(Float64, str_list[11]);
            data[row_id, 11] = is_quality_valid;

            ambig_state = parse(Float64, str_list[12]);
            data[row_id, 12] = ambig_state;

            x_rms = parse(Float64, str_list[13]);
            data[row_id, 13] = x_rms;

            y_rms = parse(Float64, str_list[14]);
            data[row_id, 14] = y_rms;

            invalid_state = parse(Float64, str_list[15]);
            data[row_id, 15] = invalid_state;

            pdh0 = parse(Float64, str_list[16]);
            data[row_id, 16] = pdh0;

            vx_rms = parse(Float64, str_list[17]);
            data[row_id, 17] = vx_rms;

            vy_rms = parse(Float64, str_list[18]);
            data[row_id, 18] = vy_rms;
        end

        if file_id == 1
            sheet = xf[file_id];
            XLSX.rename!(sheet, "$file_id.pcd");
        else
            XLSX.addsheet!(xf, "$file_id.pcd")
        end
        sheet = xf["$file_id.pcd"]

        sheet["A1", dim=2] = [
            "x";
            "y ";
            "z ";
            "dyn_prop ";
            "id ";
            "rcs ";
            "vx ";
            "vy ";
            "vx_comp";
            "vy_comp ";
            "is_quality_valid";
            "ambig_state";
            "x_rms";
            "y_rms ";
            "invalid_state";
            "pdh0";
            "vx_rms";
            "vy_rms";
        ];
        sheet["A2:R$(n - 11)"] = data;
    end
end