text:
# dict
h = Dict()
for i=1:10000
    h[i] = i+1
end
for i=1:10000
    @test (h[i] == i+1)
end
for i=1:2:10000
    delete!(h, i)
end
for i=1:2:10000
    h[i] = i+1
end
for i=1:10000
    @test (h[i] == i+1)
end
for i=1:10000
    delete!(h, i)
end
@test isempty(h)
h[77] = 100
@test h[77]==100
for i=1:10000
    h[i] = i+1
end
for i=1:2:10000
    delete!(h, i)
end
for i=10001:20000
    h[i] = i+1
end
for i=2:2:10000
    @test h[i]==i+1
end
for i=10000:20000
    @test h[i]==i+1
end
h = {"a" => 3}
@test h["a"] == 3
h["a","b"] = 4
@test h["a","b"] == h[("a","b")] == 4
h["a","b","c"] = 4
@test h["a","b","c"] == h[("a","b","c")] == 4

let
    z = Dict()
    get_KeyError = false
    try
        z["a"]
    catch _e123_
        get_KeyError = isa(_e123_,KeyError)
    end
    @test get_KeyError
end

_d = {"a"=>0}
@test isa([k for k in filter(x->length(x)==1, collect(keys(_d)))], Vector{Any})

# issue #1821
let
    d = Dict{UTF8String, Vector{Int}}()
    d["a"] = [1, 2]
    @test_throws MethodError d["b"] = 1
    @test isa(repr(d), String)  # check that printable without error
end

# issue #2344
let
    local bar
    bestkey(d, key) = key
    bestkey{K<:String,V}(d::Associative{K,V}, key) = string(key)
    bar(x) = bestkey(x, :y)
    @test bar([:x => [1,2,5]]) == :y
    @test bar(["x" => [1,2,5]]) == "y"
end

# issue #1438
type I1438T
    id
end
import Base.hash
hash(x::I1438T) = x.id

begin
    local seq, xs, s
    seq = [26,28,29,30,31,32,33,34,35,36,-32,-35,-34,-28,37,38,39,40,-30,
           -31,41,42,43,44,-33,-36,45,46,47,48,-37,-38,49,50,51,52,-46,-50,53]
    xs = [ I1438T(id) for id=1:53 ]
    s = Set()
    for id in seq
        if id > 0
            x = xs[id]
            push!(s, x)
            @test in(x, s)                 # check that x can be found
        else
            delete!(s, xs[-id])
        end
    end
end

@test  isequal(Dict(), Dict())
@test  isequal({1 => 1}, {1 => 1})
@test !isequal({1 => 1}, {})
@test !isequal({1 => 1}, {1 => 2})
@test !isequal({1 => 1}, {2 => 1})

# Generate some data to populate dicts to be compared
data_in = [ (rand(1:1000), randstring(2)) for _ in 1:1001 ]

# Populate the first dict
d1 = Dict{Int, String}()
for (k,v) in data_in
    d1[k] = v
end
data_in = collect(d1)
# shuffle the data
for i in 1:length(data_in)
    j = rand(1:length(data_in))
    data_in[i], data_in[j] = data_in[j], data_in[i]
end
# Inserting data in different (shuffled) order should result in
# equivalent dict.
d2 = Dict{Int, String}()
for (k,v) in data_in
    d2[k] = v
end

@test  isequal(d1, d2)
d3 = copy(d2)
d4 = copy(d2)
# Removing an item gives different dict
delete!(d1, data_in[rand(1:length(data_in))][1])
@test !isequal(d1, d2)
# Changing a value gives different dict
d3[data_in[rand(1:length(data_in))][1]] = randstring(3)
!isequal(d1, d3)
# Adding a pair gives different dict
d4[1001] = randstring(3)
@test !isequal(d1, d4)

@test isequal(Dict(), sizehint(Dict(),96))

# Here is what currently happens when dictionaries of different types
# are compared. This is not necessarily desirable. These tests are
# descriptive rather than proscriptive.
@test !isequal({1 => 2}, {"dog" => "bone"})
@test isequal(Dict{Int, Int}(), Dict{String, String}())

# get! (get with default values assigned to the given location)

let f(x) = x^2,
    d = {8=>19},
    def = {}

    @test get!(d, 8, 5) == 19
    @test get!(d, 19, 2) == 2

    @test get!(d, 42) do  # d is updated with f(2)
        f(2)
    end == 4

    @test get!(d, 42) do  # d is not updated
        f(200)
    end == 4

    @test get(d, 13) do   # d is not updated
        f(4)
    end == 16

    @test d == {8=>19, 19=>2, 42=>4}
end

# issue #2540
d = {x => 1
    for x in ['a', 'b', 'c']}
@test d == {'a'=>1, 'b'=>1, 'c'=> 1}

# issue #2629
d = (String => String)[ a => "foo" for a in ["a","b","c"]]
@test d == ["a"=>"foo","b"=>"foo","c"=>"foo"]

# issue #5886
d5886 = Dict()
for k5886 in 1:11
   d5886[k5886] = 1
end
for k5886 in keys(d5886)
   # undefined ref if not fixed
   d5886[k5886] += 1
end

# ############# end of dict tests #############

# #################### set ####################

# show

# isempty
@test  isempty(Set())
@test !isempty(Set([1]))
@test !isempty(Set(["banana", "apple"]))
@test !isempty(Set({1, 1:10, "pear"}))

# ordering
@test Set() < Set([1])
@test Set([1]) < Set([1,2])
@test !(Set([3]) < Set([1,2]))
@test !(Set([3]) > Set([1,2]))
@test Set([1,2,3]) > Set([1,2])
@test !(Set([3]) <= Set([1,2]))
@test !(Set([3]) >= Set([1,2]))
@test Set([1]) <= Set([1,2])
@test Set([1,2]) <= Set([1,2])
@test Set([1,2]) >= Set([1,2])
@test Set([1,2,3]) >= Set([1,2])
@test !(Set([1,2,3]) >= Set([1,2,4]))
@test !(Set([1,2,3]) <= Set([1,2,4]))

# add, length
s = Set()
@test isempty(s)
for i in 1:1000
    push!(s, i)
    @test (length(s) == i)
end

# delete!, has, in
for i in 1:2:1000
    delete!(s, i)
end
for i in 1:2:1000
    @test !in(i  , s)
    @test  in(i+1, s)
end

# elements
data_in = (1,"banana", ())
s = Set(data_in)
data_out = collect(s)
@test is(typeof(data_out), Array{Any,1})
@test all(map(d->in(d,data_out), data_in))
@test all(map(data_in) do d
              in(d,data_out)
          end)
@test length(data_out) == length(data_in)

# homogeneous sets
@test is(typeof(Set([1,2,3])), Set{Int})
@test is(typeof(Set{Int}([3])), Set{Int})

# eltype
@test is(eltype(Set({1,"hello"})), Any)
@test is(eltype(Set{String}()), String)

# no duplicates
s = Set([1,2,3])
@test length(s) == 3
push!(s,2)
@test length(s) == 3
delete!(s,2)
@test length(s) == 2
push!(s,2)
@test length(s) == 3
pop!(s,2)
@test length(s) == 2
pop!(s)
@test length(s) == 1

# union
s = union(Set(1,2), Set(3,4))
@test isequal(s, Set(1,2,3,4))
s = union(Set(5,6,7,8), Set(7,8,9))
@test isequal(s, Set(5,6,7,8,9))

# intersect
s = intersect(Set(1,2), Set(3,4))
@test isequal(s, Set())
s = intersect(Set(5,6,7,8), Set(7,8,9))
@test isequal(s, Set(7,8))
@test isequal(intersect(Set(2,3,1), Set(4,2,3), Set(5,4,3,2)), Set(2,3))

# setdiff
@test isequal(setdiff(Set(1,2,3), Set()), Set(1,2,3))
@test isequal(setdiff(Set(1,2,3), Set(1)),  Set(2,3))
@test isequal(setdiff(Set(1,2,3), Set(1,2)),  Set(3))
@test isequal(setdiff(Set(1,2,3), Set(1,2,3)), Set())
@test isequal(setdiff(Set(1,2,3), Set(4)),  Set(1,2,3))
@test isequal(setdiff(Set(1,2,3), Set(4,1)),  Set(2,3))

for (l,r) in ((Set(1,2),     Set(3,4)),
              (Set(5,6,7,8), Set(7,8,9)),
              (Set(1,2),     Set(3,4)),
              (Set(5,6,7,8), Set(7,8,9)),
              (Set(1,2,3),   Set()),
              (Set(1,2,3),   Set(1)),
              (Set(1,2,3),   Set(1,2)),
              (Set(1,2,3),   Set(1,2,3)),
              (Set(1,2,3),   Set(4)),
              (Set(1,2,3),   Set(4,1)))
    @test issubset(intersect(l,r), l)
    @test issubset(intersect(l,r), r)
    @test issubset(l, union(l,r))
    @test issubset(r, union(l,r))
    @test isequal(union(intersect(l,r),symdiff(l,r)), union(l,r))
end

@test setdiff(IntSet(1, 2, 3, 4), IntSet(2, 4, 5, 6)) == IntSet(1, 3)
@test setdiff(Set(1, 2, 3, 4), Set(2, 4, 5, 6)) == Set(1, 3)

@test symdiff(IntSet(1, 2, 3, 4), IntSet(2, 4, 5, 6)) == IntSet(1, 3, 5, 6)
@test symdiff(Set(1, 2, 3, 4), Set(2, 4, 5, 6)) == Set(1, 3, 5, 6)

s1 = Set(1, 2, 3, 4)
setdiff!(s1, Set(2, 4, 5, 6))

@test s1 == Set(1, 3)

s2 = IntSet(1, 2, 3, 4)
setdiff!(s2, IntSet(2, 4, 5, 6))

@test s2 == IntSet(1, 3)

# union!
s = Set(1,3,5,7)
union!(s,(2,3,4,5))
@test isequal(s,Set(1,2,3,4,5,7))

# setdiff!
s = Set(1,3,5,7)
setdiff!(s,(3,5))
@test isequal(s,Set(1,7))

# similar
s = similar(Set([1,"Banana"]))
@test length(s) == 0
@test typeof(s) == Set{Any}
s = similar(Set{Float32}([2.0f0,3.0f0,4.0f0]))
@test length(s) == 0
@test typeof(s) == Set{Float32}

# copy
data_in = (1,2,9,8,4)
s = Set(data_in)
c = copy(s)
@test isequal(s,c)
push!(s,100)
push!(c,200)
@test !in(100, c)
@test !in(200, s)

# start, done, next
for data_in in ((7,8,4,5),
                ("hello", 23, 2.7, (), [], (1,8)))
    s = Set(data_in)

    s_new = Set()
    for el in s
        push!(s_new, el)
    end
    @test isequal(s, s_new)

    t = tuple(s...)
    @test length(t) == length(s)
    for e in t
        @test in(e,s)
    end
end

# zip
let i = 0
x = 1:2:8
y = 2:2:8
xy = 1:8
for (thisx, thisy) in zip(x, y)
    @test thisx == xy[i+=1]
    @test thisy == xy[i+=1]
end
end

# pop!
origs = Set([1,2,3,"apple"])
s = copy(origs)
for i in 1:length(origs)
    el = pop!(s)
    @test !in(el, s)
    @test in(el, origs)
end
@test isempty(s)

# hash
s1 = Set{ASCIIString}(["bar", "foo"])
s2 = Set{ASCIIString}(["foo", "bar"])
s3 = Set{ASCIIString}(["baz"])
@test hash(s1) == hash(s2)
@test hash(s1) != hash(s3)

# isequal
@test  isequal(Set(), Set())
@test !isequal(Set(), Set(1))
@test  isequal(Set{Any}({1,2}), Set{Int}([1,2]))
@test !isequal(Set{Any}({1,2}), Set{Int}([1,2,3]))

# Comparison of unrelated types seems rather inconsistent

@test  isequal(Set{Int}(), Set{String}())
@test !isequal(Set{Int}(), Set{String}([""]))
@test !isequal(Set{String}(), Set{Int}([0]))
@test !isequal(Set{Int}([1]), Set{String}())

@test  isequal(Set{Any}([1,2,3]), Set{Int}([1,2,3]))
@test  isequal(Set{Int}([1,2,3]), Set{Any}([1,2,3]))

@test !isequal(Set{Any}([1,2,3]), Set{Int}([1,2,3,4]))
@test !isequal(Set{Int}([1,2,3]), Set{Any}([1,2,3,4]))

@test !isequal(Set{Any}([1,2,3,4]), Set{Int}([1,2,3]))
@test !isequal(Set{Int}([1,2,3,4]), Set{Any}([1,2,3]))

# ########## end of set tests ##########

## IntSet

s = IntSet(0,1,10,20,200,300,1000,10000,10002)
@test last(s) == 10002
@test first(s) == 0
@test length(s) == 9
@test pop!(s) == 10002
@test length(s) == 8
@test shift!(s) == 0
@test length(s) == 7
@test !in(0,s)
@test !in(10002,s)
@test in(10000,s)
@test_throws ErrorException first(IntSet())
@test_throws ErrorException last(IntSet())

# Ranges

@test isempty((1:4)[5:4])
@test_throws BoundsError (1:10)[8:-1:-2]

max_stars_repo_name: tobig/julia

max_issues_repo_name: tobig/julia

max_forks_repo_name: tobig/julia


text:
#gets the index of the first high value for each pulse in a sequence
#if pad_first is true then the first index counts if it is high
function find_pulse_starts(pulses::AbstractVector{Bool}; pad_first = true)
    output = Int[]
    if pad_first && pulses[1]
        push!(output, 1)
    end
    is_high = pulses[1]
    for i = 2:length(pulses)
        if !is_high
            if pulses[i]
                push!(output, i)
                is_high = true
            end
        elseif !pulses[i]
                is_high = false
        end
    end
    return output
end

#gets the index of the last high value for each pulse in a sequence
#if pad_last is true then the last index counts if it is high
function find_pulse_stops(pulses::AbstractVector{Bool}; pad_last = true)
    output = Int[]
    is_low = !pulses[1]
    for i = 2:length(pulses)
        if !is_low
            if !pulses[i]
                push!(output, i-1)
                is_low = true
            end
        elseif pulses[i]
                is_low = false
        end
    end
    if pad_last && pulses[end]
        push!(output, length(pulses))
    end
    return output
end

find_pulse_starts(pulses::AbstractVector{Bool}, thresh) = find_pulse_starts(pulses)
find_pulse_stops(pulses::AbstractVector{Bool}, thresh) = find_pulse_stops(pulses)
find_pulse_starts(pulses::AbstractVector{Bool}, thresh::Bool) = find_pulse_starts(pulses)
find_pulse_stops(pulses::AbstractVector{Bool}, thresh::Bool) = find_pulse_stops(pulses)

find_pulse_starts(pulses::AbstractVector{T}, thresh::T) where T = find_pulse_starts(pulses.>=thresh)
find_pulse_stops(pulses::AbstractVector{T}, thresh::T) where T = find_pulse_stops(pulses.>=thresh)
find_pulse_starts(com::ImagineSignal; thresh = 1.15 * Unitful.V, sampmap = :world) = find_pulse_starts(get_samples(com; sampmap=sampmap).data, thresh)
find_pulse_stops(com::ImagineSignal; thresh = 1.15 * Unitful.V, sampmap = :world) = find_pulse_stops(get_samples(com; sampmap=sampmap).data, thresh)

count_pulses(pulses::AbstractVector{T}, thresh::T) where T = length(find_pulse_starts(pulses, thresh))

#default threshold is half of 3.3V TTL pulse -> 1.15V
count_pulses(com::ImagineSignal; thresh = isdigital(com) ? true : 1.15 * Unitful.V, sampmap = :world) = count_pulses(get_samples(com; sampmap = sampmap).data, thresh)

ext: jl

max_forks_repo_name: HolyLab/ImagineInterface


text:
"""
    filter_exp(N, Nc, s, V, invV)

Construct exponential filter for nodal solution

- @arg N: number of coefficients
- @arg Nc: cutoff location
- @arg s: order of filter (must be even)
- @arg V: Vandermonde matrix
"""
function filter_exp(N, s, V, Nc = 0, invV = inv(V))
    nv = size(V, 1)
    if nv == N + 1
        filterdiag = KitBase.filter_exp1d(N, s, Nc)
    elseif nv == (N + 1) * (N + 2) ÷ 2
        filterdiag = filter_exp2d(N, s, Nc)
    end
    F = V * diagm(filterdiag) * invV

    return F
end

"""
    filter_exp2d(N, Nc, s)

Construct exponential filter for modal solution

- @arg N: degree of polynomials
- @arg s: order of filter (must be even)
- @arg Nc: cutoff location
"""
function filter_exp2d(Norder, sp, Nc = 0)
    alpha = -log(eps())

    filterdiag = ones((Norder + 1) * (Norder + 2) ÷ 2)
    sk = 1
    for i = 0:Norder
        for j = 0:Norder-i
            if i + j >= Nc
                filterdiag[sk] = exp(-alpha * ((i + j - Nc) / (Norder - Nc))^sp)
            end
            sk += 1
        end
    end

    return filterdiag
end

"""
    basis_norm(deg)

Calculate norm of polynomial basis
"""
function basis_norm(deg)
    NxHat = 100
    xHat = range(-1, stop = 1, length = NxHat) |> collect
    dxHat = xHat[2] - xHat[1]

    nLocal = deg + 1
    PhiL1 = zeros(nLocal)
    for i = 1:nLocal
        PhiL1[i] = dxHat * sum(abs.(JacobiP(xHat, 0, 0, i - 1)))
    end

    return PhiL1
end


text:
using ArgParse
using PyPlot, NNlib, StatsFuns
using Serialization, DelimitedFiles

"""
    parse_commandline()

This function parses the arguments of the call `julia banana.jl`.
"""
function parse_commandline()
    s = ArgParseSettings()

    @add_arg_table! s begin

        "datapath"
            help = "dataset path"
            required = true

        "resultspath"
            help = "results path"
            required = true

        "uid"
            help = "uid"
            required = true

        "kernelstring"
            help = "kernel string"
            required = true
    end

    return parse_args(s)
end

parsed_args = parse_commandline()

datapath = parsed_args["datapath"]
datadir = parsed_args["resultspath"]
uid = parsed_args["uid"]
kernel = parsed_args["kernelstring"]

include("functions.jl")
include("helper_functions.jl")
include("utils.jl")

samples = deserialize(joinpath(datadir, "banana_samples_$(uid)_$(kernel).jd"))
stats = deserialize(joinpath(datadir, "banana_stats_$(uid)_$(kernel).jd"))
K, D, n_samples, n_adapts, ad, kernel, seed, pν = deserialize(joinpath(datadir, "banana_setup_$(uid)_$(kernel).jd"))

X = readdlm(joinpath(datapath, "banana_datapoints_subset.csv"), ',')
y = vec(Int.(readdlm(joinpath(datapath, "banana_classes_subset.csv"), ',')))

gridlength = 3.75
xmin, ymin = -gridlength, -gridlength
xmax, ymax = gridlength, gridlength

grid = mapreduce(yi -> mapreduce(xi -> [xi, yi], hcat, range(xmin, xmax, length=300)), hcat, range(ymin, ymax, length=300))'
yhat, σhat = predict_banana(grid, samples, kernel, K)

fig1, ax1 = plt.subplots()

ax1.imshow(reshape(1 .- yhat, 300, 300)', cmap="RdBu", origin="lower", vmin = 0, vmax = 1, extent=[-gridlength, gridlength, -gridlength, gridlength])
ax1.plot(X[y .== false,1], X[y .== false,2], "o", markersize=6.5, markerfacecolor="white", markeredgewidth=1.5, markeredgecolor="blue", alpha=1.0)
ax1.plot(X[y .== true,1], X[y .== true,2], "v", markersize=6.5, markerfacecolor="white", markeredgewidth=1.5, markeredgecolor="red", alpha=1.0)

ax1.axis("off")

fig1.tight_layout()
savefig("banana_$(kernel)_expectation.pdf", bbox_inches = "tight", pad_inches = 0)

display(gcf())

max_stars_repo_name: AaltoML/PeriodicBNN

max_issues_repo_name: AaltoML/PeriodicBNN

max_forks_repo_name: AaltoML/PeriodicBNN


text:
#Run the Chronic simulations in parallel
#change the number of cores used as needed
using Distributed
addprocs(40)
#The @everywhere macro runs the command/defines the thing on all the processors
@everywhere include("ProblemGenerator.jl")

#Given a new parameter set, determine average IFN at end of simulation
@everywhere function ChronicIFN(kcat8,τ7,prob)
    #Make a copy of the parameter container
    θ = deepcopy(prob.p)

    #Need to check if each cell has a unique parameter set
    if isa(prob.p.par[1],Array)
      #Set the parameters to desired input
      θ.par[13] .= kcat8
      θ.par[27] .= τ7
    else
      #Set the parameters to desired input
      θ.par[13] = kcat8
      θ.par[27] = τ7
    end

    @show (kcat8,τ7)
    #Generate the new problem (run for longer to ensure steady state)
    probNew = remake(prob; p=θ,tspan=(0.0,168.0))
    sol = solve(probNew,CVODE_BDF(linear_solver=:GMRES),saveat=1.0)

    return mean(sol[:,:,7,end])
end

#Generate a problem
@everywhere prob = ModelSetup(:ISD,:notStochastic,:Hetero)

#We need to know if every cell has been given a unique parameter set
if isa(prob.p.par[1],Array) #Is the first parameter a single number or array?
  kcat8True = prob.p.par[13][1]
  τ7True = prob.p.par[27][1]
else
  kcat8True,τ7True = prob.p.par[[13,27]]
end

#Vary STAT production (kcat8) and IFN Degradation (τ7)
parRange = 20
kcat8Vals = range(0.5*kcat8True,2.0*kcat8True,length=parRange)
τ7Vals = range(0.5*τ7True,2.0*τ7True,length=parRange)

iter = Iterators.product(kcat8Vals,τ7Vals)
meanIFNβ = pmap(x -> ChronicIFN(x[1],x[2],prob),iter)

#plt = heatmap(kcat8Vals,τ7Vals,meanIFNβ,color=:viridis)

kvals = zeros(parRange^2,2)
for ((kcat8,τ7),i) in zip(iter,1:parRange^2)
    kvals[i,:] = [kcat8,τ7]
end

chronicIFN = DataFrame(kcat8 = kvals[:,1],tau7 = kvals[:,2], IFN = reverse(meanIFNβ[:]))

#Save the simulation to plot
CSV.write("./ServerSimulations/VirusFigure6_Det_Hetero.csv",chronicIFN)


text:
#================================================

AbstractNode API definition, and typedefs for
child types

=================================================#

using Base.Threads

#===========================
Type definitions
===========================#

struct Node{T,M} <: AbstractNode{T,M}
    vals::Vector{T}
    next::Vector{Node{T,M}}
    capacity::Int64
    flags::UInt8
end

#===========================
Shared AbstractNode constructors
===========================#

for dtype in (:Node,)
    @eval begin
        $dtype{M}(val; kws...) where {T,M} =
            $dtype{eltype(val),M}(val; kws...)

        $dtype{M}(val, height; kws...) where {T,M} =
            $dtype{eltype(val),M}(val, height; kws...)

        $dtype{T,M}(val; p = DEFAULT_P, max_height = DEFAULT_MAX_HEIGHT, kws...) where {T,M} =
            $dtype{T,M}(val, random_height(p; max_height=max_height); kws...)
    end
end

#===========================
Shared AbstractNode API
===========================#

height(node::AbstractNode) = length(node.next)

Base.show(io::IO, node::AbstractNode) = write(io, string(node))
Base.display(node::AbstractNode) = println(string(node))

# Node links

function link_nodes!(src, dst, level)
    src.next[level] = dst
end

next(src::AbstractNode, level) = src.next[level]

# Flags

has_flag(node, flag) = (node.flags & flag) != 0
is_left_sentinel(node) = has_flag(node, FLAG_IS_LEFT_SENTINEL)
is_right_sentinel(node) = has_flag(node, FLAG_IS_RIGHT_SENTINEL)
is_sentinel(node) = has_flag(node, IS_SENTINEL)

# Node comparison

Base.:(<)(node::AbstractNode, val) = !(val ≤ node)
Base.:(<)(val, node::AbstractNode) = !(node ≤ val)
Base.:(<)(node_1::AbstractNode, node_2::AbstractNode) = !(node_2 ≤ node_1)

Base.:(<=)(node::AbstractNode, val) =
    is_sentinel(node) ? is_left_sentinel(node) : (key(node) ≤ val)

Base.:(<=)(val, node::AbstractNode) =
    is_sentinel(node) ? is_right_sentinel(node) : (val ≤ key(node))

function Base.:(<=)(node_1::AbstractNode, node_2::AbstractNode)
    if is_sentinel(node_1)
        is_left_sentinel(node_1) || is_right_sentinel(node_2)
    elseif is_sentinel(node_2)
        is_right_sentinel(node_2)
    else
        key(node_1) ≤ key(node_2)
    end
end

Base.:(==)(node::AbstractNode, val) = is_sentinel(node) ? false : key(node) == val
Base.:(==)(val, node::AbstractNode) = (node == val)
Base.:(==)(node_1::AbstractNode, node_2::AbstractNode) =
    (is_sentinel(node_1) || is_sentinel(node_2)) ?
    false :
    key(node_1) == key(node_2)


text:
using Serialization
using Parameters
using Setfield
using Dates
using Suppressor
using LinearAlgebra

""" Configuration struct for running experiments """
@with_kw struct TrainConfig
    # metadata
    name::String
    result_path::String = "results"
    data_path::String
    log_path::String = ""
    dtype::Type = Float64
    # assume by default that the system has two threads per core
    blas_threads::Int = Integer(0.5 * length(Sys.cpu_info()))
    # train configuration
    num_epochs::Int
    p_init_fn = initial_params
    patience::Int
    verbose::Bool = true
    logevery::Int = 1
    # model configuration
    ode::O where {O<:Ode}
    link::L where {L<:Link}
    mod::M where {M<:Modifier} = Nominal()
    feed::F where {F<:Feedthrough} = NoFeedthrough()
    # DynamicsTrainer arguments
    solver::Any
    adaptive::Bool = true
    sensealg::Any
    abstol = 1e-6
    reltol = 1e-3
    dt = 1.0
    initial_condition = solve_ss
    loss = Flux.Losses.mse
    regularization = no_regularization
    regularization_coef = 0.0
    optim = ADAMW(dtype(1e-3))
    wdecay = WeightDecay(dtype(0.0))
    sched = Step(λ = dtype(1e-3), γ = dtype(1.0), 1)
    validation_split::Float64 = 0.0
    test_split::Float64 = 0.0
    train_split::Float64 = 1.0
    batch_size::Int
    time_dropout::Float64
    time_horizon::Int
    time_horizon_increment::Int
    time_horizon_max::Int
    callback_tasks = []
end

""" Unpack and run model train configuration """
function run_config(c::TrainConfig; runs = 1)
    @unpack_TrainConfig c

    # set number of threads to be used by BLAS
    LinearAlgebra.BLAS.set_num_threads(blas_threads)

    # display multithreading information
    println(
        "Number of BLAS threads: $(ccall((:openblas_get_num_threads64_, Base.libblas_name), Cint, ())
), Number of Julia threads: $(Threads.nthreads())",
    )

    for _ = 1:runs
        # create experiment result path
        savepath = joinpath(result_path, name, Dates.format(now(), "HH_MM_SS_MS"))
        if ~ispath(savepath)
            mkpath(savepath)
        end

        # write config
        serialize(joinpath(savepath, "config.jls"), c)

        # write human readable config as txt file
        config_string = @capture_out show(c)
        open(joinpath(savepath, "config.txt"), "w") do io
            write(io, config_string)
        end

        # read data
        dataset = deserialize(data_path)
        t, u, y = dataset["t"], dataset["u"], dataset["y"]

        # model training
        trainer = DynamicsTrainer(
            ode,
            link,
            mod,
            feed,
            solver,
            sensealg,
            Array{dtype,1}(t),
            Array{dtype,3}(u),
            Array{dtype,3}(y);
            savepath = savepath,
            dtype = dtype,
            validation_split = validation_split,
            test_split = test_split,
            abstol = abstol,
            reltol = reltol,
            dt = dt,
            adaptive = adaptive,
            initial_condition = initial_condition,
            loss = loss,
            regularization = regularization,
            regularization_coef = regularization_coef,
            optim = optim,
            wdecay = wdecay,
            sched = sched,
            batch_size = batch_size,
            callback_tasks = callback_tasks,
        )

        train(
            trainer,
            num_epochs;
            p_init = p_init_fn(trainer; log_path = log_path),
            time_dropout = time_dropout,
            time_horizon = time_horizon,
            time_horizon_increment = time_horizon_increment,
            time_horizon_max = time_horizon_max,
            patience = patience,
            logevery = logevery,
            verbose = verbose,
        )
    end
end

max_stars_repo_name: syanga/ISSNeuralODE

max_issues_repo_name: syanga/ISSNeuralODE

max_forks_repo_name: syanga/ISSNeuralODE


text:
using Potential, Preprocessing, Optimization

function VelocityVerlet(x, v, q, t, a, b, c, pbc, app, potentials)

    dim = Int32(size(x,1))
    inum = Int32(length(t))

    ensemblemode = app.ensemblemode;
    time = app.time;             # initial time
    dt = app.dt
    ntimesteps =  app.ntimesteps;       # # time steps
    globalfreq =  app.globalfreq;
    skin = app.neighskin;
    every = app.neighevery;      # perform neighbor rebuild list for "every" iterations
    delay = app.neighdelay;      # delay neighbor rebuild
    distcheck = app.neighcheck;  # 0 -> neighbor rebuild by delay and every settings,
    outputflag = app.outputflag
    eta = app.eta[:]
    kappa = app.kappa[:]
    atommass = [0; app.atommasses[:]]

    if (lowercase(ensemblemode) == "nve")
        ensemblemodenum = 0;
    elseif (lowercase(ensemblemode) == "nvelimit")
        ensemblemodenum = 1;
    elseif (lowercase(ensemblemode) == "nvt")
        ensemblemodenum = 2;
    else
        ensemblemodenum = -1;
    end
    ensemblemodenum = Int32(ensemblemodenum)

    units = Preprocessing.setunits(app.unitstyle)
    boltz = Int32(units.boltz)
    mvv2e = Int32(units.mvv2e)
    nktv2p = Int32(units.nktv2p)

    dtarray, tarray, vlimitsq,  eta_mass, eta_thermostat, eta_dot, eta_dotdot,
        eta_mass_flag, biasflag, mtchain, nc_tchain = ensemble(app, dim, inum)

    dom = setdomainstruct(pbc, a, b, c)
    boxlo = dom.boxlo
    boxhi = dom.boxhi
    hi_lambda = dom.boxhi_lamda
    lo_lambda = dom.boxlo_lamda
    h = dom.h
    h_inv = dom.h_inv
    h_rate = dom.h_rate
    vdeform = dom.vdeform
    triclinic = dom.triclinic
    inv_volume = 1.0 / (dom.h[1] * dom.h[2] * dom.h[3])

    rcutmax = skin + Potential.maximumcutoff(potentials)
    x, alist, neighlist, neighnum = Potential.fullneighborlist(x, a, b, c, pbc, rcutmax);
    etot, eatom, f, vatom, vtot = Potential.emlpotential(x, q, t, alist, neighlist, neighnum, eta, kappa, potentials)
    pe, temp, pres, ke = thermooutput(v, t, alist, atommass, vtot, etot, mvv2e, boltz, nktv2p, inv_volume, outputflag)
    print("Total number of atoms = $inum\n");
    print("Time step     Temperature        Potential energy        Total energy            Pressure \n");
    printoutput(0, ntimesteps, temp, pe, ke, pres)

    xhold = 1.0*x
    for istep=1:ntimesteps
        time = time + dt
        gnum = Int32(length(alist))
        dtarray[6] = istep
        xstart = 1.0*x

        InitialIntegrate(x, v, f, atommass, dtarray, tarray, eta_mass, eta_thermostat,
            eta_dot, eta_dotdot, vlimitsq, Int32.(t[:]), Int32.(alist.-1), eta_mass_flag, biasflag, mtchain,
            nc_tchain, ensemblemodenum, dim, inum);

        PostInitialIntegration(x, v, f, t, eatom, vatom, atommass, dtarray, app.setvelocity, app.wallreflect)

        build = checkbuild(x, xhold, skin, delay, every, distcheck, istep, dim, inum)
        if build==1
            print("Rebuild neighbor list at timestep $istep...\n");
            if sum(pbc) > 0
                ApplyPBC(x, v, image, boxhi, boxlo, hi_lambda, lo_lambda, h, h_inv,
                    h_rate, pbc, vdeform, triclinic, dim, inum)
            end
            x, alist, neighlist, neighnum = Potential.fullneighborlist(x[:,1:inum], a, b, c, pbc, rcutmax);
            xhold = 1.0*x
        else
            if sum(pbc) > 0
                # calculate the distance between x and xstart: xstart = x - xstart
                xstart = x - xstart
                # shift periodic images by the above distance
                xtm = x[:,(inum+1):end]
                atm = alist[(inum+1):end]
                ArrayPlusAtColumnIndex(xtm, xstart, Int32.(atm.-1), dim, gnum);
                x[:,(inum+1):end] = xtm
            end
        end

        etot, eatom, f, vatom, vtot = Potential.emlpotential(x, q, t, alist, neighlist, neighnum, eta, kappa, potentials)

        PostForceComputation(x, v, f, eatom, vatom, app.setforce, app.lineforce, app.planeforce,
                    app.wallharmonic, app.walllj126, app.wallmorse)

        FinalIntegrate(x, v, f, atommass, dtarray, tarray, eta_mass, eta_thermostat, eta_dot,
            eta_dotdot, vlimitsq, Int32.(t[:]), Int32.(alist.-1), eta_mass_flag, biasflag,
            mtchain, nc_tchain,  ensemblemodenum, dim, inum);

        PostFinalIntegration(x, v, f, t, atommass, dtarray, tarray,
                biasflag, app.setvelocity, app.velocityrescaling)

        if (istep % globalfreq == 0)
            pe, temp, pres, ke = thermooutput(v, t, alist, atommass, vtot, etot,
                mvv2e, boltz, nktv2p, inv_volume, outputflag)
            printoutput(istep, ntimesteps, temp, pe, ke, pres)
        end
    end
end

function VelocityVerlet(config, app, potentials, ci=1)

    natom = [0; cumsum(config.natom[:])];

    x = config.x[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
    v = config.v[:,(natom[ci]+1):natom[ci+1]] # atom positions of configuration ci
    t = config.t[:,(natom[ci]+1):natom[ci+1]] # atom types of configuration ci
    if length(config.q)>0
        q = config.q[:,(natom[ci]+1):natom[ci+1]]
    else
        q = []
    end
    if size(config.a,2) > 1
        a = config.a[:,ci]
        b = config.b[:,ci]
        c = config.c[:,ci]
    else
        a = config.a[:,1]
        b = config.b[:,1]
        c = config.c[:,1]
    end
    if size(config.pbc,2) > 1
        pbc = config.pbc[:,ci]
    else
        pbc = config.pbc[:,1]
    end

    VelocityVerlet(x, v, q, t, a, b, c, pbc, app, potentials)
end

function printoutput(istep, ntimesteps, temp, pe, ke, pres)
    n = length(digits(ntimesteps))
    s0 = repeat(' ', n)
    s1 = replacestring(s0, string(istep)) * "        "
    s2 = num2string(temp, 16) * "      "
    s3 = num2string(pe, 16) * "      "
    s4 = num2string(pe+ke, 16) * "      "
    s5 = num2string(pres, 16) * "      "
    print(s1 * s2 * s3 * s4 * s5 * "\n");
end

function ensemble(app, dim, inum)

    ensemblemode = app.ensemblemode
    dtarray = zeros(10)
    tarray = zeros(10)
    eta_mass = zeros(10)
    eta_thermostat = zeros(10)
    eta_dot = zeros(10)
    eta_dotdot = zeros(10)

    units = Preprocessing.setunits(app.unitstyle)
    boltz = units.boltz
    mvv2e = units.mvv2e
    ftm2v = units.ftm2v
    tdof = (inum-1)*dim;

    dt = app.dt
    dtarray[1] = dt;     # dt
    dtarray[2] = 0.5*dt*ftm2v; # dtf
    dtarray[3] = dt;     # dtv
    dtarray[4] = 0;             # beginstep
    dtarray[5] = app.ntimesteps; # endstep
    vlimitsq = 0.0;
    eta_mass_flag = 1;
    biasflag = 0;
    mtchain = 1;
    nc_tchain = 1;
    if ensemblemode != "nve"
        if (ensemblemode == "nvelimit")
            nveparam = app.nveparam
            vlimitsq = (nveparam[1]/dt) * (nveparam[1]/dt);
        end
        if (ensemblemode == "nvt")  # NVT ensemble
            nvtparam = app.nvtparam
            tarray[1] = nvtparam[0]; # temp start
            tarray[2] = nvtparam[1]; # temp stop
            tarray[3] = 1/nvtparam[2]; # temp frequency
            #delta = (ntimestep - beginstep)/(endstep - beginstep)
            #tarray[4] = t_start + delta * (t_stop - t_start); # t_target
            #t_current = cpuComputeTempScalar(v, mass, tfactor, type, ilist, dim, inum);
            #tarray[5] = t_current;
            tarray[6] = tdof;
            tarray[7] = boltz;
            tarray[8] = (length(nvtparam)>3) ? nvtparam[3] : 0.0; # drag factor added to barostat/thermostat
            tarray[9] = mvv2e;
            mtchain   = (length(nvtparam)>4) ? Int32(nvtparam[4]) : 3;
            nc_tchain = (length(nvtparam)>5) ? Int32(nvtparam[5]) : 1;
        end
    end

    eta_mass_flag = Int32(eta_mass_flag)
    biasflag = Int32(biasflag)
    mtchain = Int32(mtchain)
    nc_tchain = Int32(nc_tchain)
    return dtarray, tarray, vlimitsq,  eta_mass, eta_thermostat, eta_dot, eta_dotdot, eta_mass_flag, biasflag, mtchain, nc_tchain
end

function thermooutput(v, t, ilist, atommass, vtot, etot, mvv2e, boltz, nktv2p, inv_volume, flag=0)

    dim, inum = size(v)
    tdof = (inum-1)*dim;
    tfactor = mvv2e/(tdof * boltz)

    # function ComputeKEAtom(ke::Vector{Float64}, mass::Vector{Float64}, v::Vector{Float64},
    #     mvv2e::Float64, type::Vector{Int32}, ilist::Vector{Int32}, dim::Int32, inum::Int32)

    keatom = zeros(inum)
    ComputeKEAtom(keatom, atommass, v[:], 2.0, Int32.(t[:]), Int32.(ilist .- 1), Int32(dim), Int32(inum));
    kesum = sum(keatom)

    pe = etot
    temp = tfactor*kesum
    ke = 0.5 * tdof * boltz * temp
    pres = (tdof * boltz * temp + vtot[1] + vtot[2] + vtot[3]) / 3.0 * inv_volume * nktv2p

    if (flag == 0)
        pe = pe/inum
        ke = ke/inum
    end

    return pe, temp, pres, ke
end

function PostInitialIntegration(x, v, f, t, eatom, vatom, mass,
            dtarray, setvelocity=nothing, wallreflect=nothing)

    dim, anum = size(v)
    alist = Int32.(0:(anum-1))
    eflag_atom = 0;
    vflag_atom = 0;
    dtf = dtarray[2]
    dtv = dtarray[3]

    x = x[:]
    v = v[:]
    f = f[:]
    if (setvelocity !== nothing)
        nsetvelocity = length(setvelocity)
        for  i = 1:nsetvelocity
            fparam = setvelocity[i].fparam
            iparam = setvelocity[i].iparam
            ind =  findall(t == setvelocity[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            SetVelocityInitialIntegrate(x, v, f, mass, fparam, dtf,
                    dtv, t, ilist, iparam, dim, inum);
        end
    end

    if (wallreflect !== nothing)
        nwallreflect = length(wallreflect)
        for i = 1:nwallreflect
            fparam = wallreflect[i].fparam
            iparam = wallreflect[i].iparam
            ind =  findall(t == wallreflect[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixWallReflect(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end
end

function checkbuild(x, xhold, skin, delay, every, distcheck, istep, dim, inum)
    build = 0;
    ago = istep + 1;

    if ((ago >= delay) && (ago % every == 0))
        if (distcheck == 0)
            build = 1;
        else
            ds = zeros(inum)
            ArrayDistSquareSum(ds, x[:], xhold[:], dim, inum);
            maxdistsquare = maximum(ds)
            if (maxdistsquare > 0.25*skin*skin)
                build = 1;
            end
        end
    end

    return build
end

function PostForceComputation(x, v, f, t, eatom, vatom, setforce=nothing, lineforce=nothing,
    planeforce=nothing, wallharmonic=nothing, walllj126=nothing, wallmorse=nothing)

    dim, anum = size(v)
    alist = Int32.(0:(anum-1))

    eflag_atom = 0;
    vflag_atom = 0;

    if (setforce !== nothing)
        nsetforce = length(setforce)
        for i = 1:nsetforce
            fparam = setforce[i].fparam
            iparam = setforce[i].iparam
            ind =  findall(t == setforce[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixSetForce(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end

    if (lineforce !== nothing)
        nlineforce = length(lineforce)
        for i = 1:nlineforce
            fparam = lineforce[i].fparam
            iparam = lineforce[i].iparam
            ind =  findall(t == lineforce[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixLineForce(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end

    if (planeforce !== nothing)
        nplaneforce = length(planeforce)
        for i = 1:nplaneforce
            fparam = planeforce[i].fparam
            iparam = planeforce[i].iparam
            ind =  findall(t == planeforce[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixPlaneForce(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end

    if (wallharmonic !== nothing)
        nwallharmonic = length(wallharmonic)
        for i = 1:nwallharmonic
            fparam = wallharmonic[i].fparam
            iparam = wallharmonic[i].iparam
            ind =  findall(t == wallharmonic[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixWallHarmonic(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end

    if (walllj126 !== nothing)
        nwalllj126 = length(walllj126)
        for i = 1:nwalllj126
            fparam = walllj126[i].fparam
            iparam = walllj126[i].iparam
            ind =  findall(t == walllj126[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixWallLJ126(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end

    if (wallmorse !== nothing)
        nwallmorse = length(wallmorse)
        for i = 1:nwallmorse
            fparam = wallmorse[i].fparam
            iparam = wallmorse[i].iparam
            ind =  findall(t == wallmorse[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            FixWallMorse(x, v, f, eatom, vatom, fparam,
                iparam, ilist, eflag_atom, vflag_atom, dim, inum);
        end
    end
end

function PostFinalIntegration(x, v, f, t, mass, dtarray, tarray,
        biasflag, setvelocity=nothing, velocityrescaling=nothing)

    dim, anum = size(v)
    alist = Int32.(0:(anum-1))
    energy = 0.0
    dtf = dtarray[2]
    #dtv = dtarray[3]

    if (setvelocity !== nothing)
        nsetvelocity = length(setvelocity)
        for  i = 1:nsetvelocity
            #fparam = setvelocity[i].fparam
            iparam = setvelocity[i].iparam
            ind =  findall(t == setvelocity[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            SetVelocityFinalIntegrate(x, v, f, mass, dtf, t, ilist, iparam, dim, inum);
        end
    end

    if (velocityrescaling !== nothing)
        nvelocityrescaling = length(velocityrescaling)
        second = zeros(1)
        seed = Int32.(ones(1))
        save = Int32.(zeros(1))
        for  i = 1:nvelocityrescaling
            fparam = velocityrescaling[i].fparam
            iparam = velocityrescaling[i].iparam
            energy = fparam[1];
            mode = iparam[1];
            ind =  findall(t == velocityrescaling[i].groupnum)
            ilist = alist[ind]
            inum = length(ilist)
            energy = VelocityRescalingThermostat(v, mass, dtarray, tarray, second, energy,
                    t, ilist, seed, save, biasflag, mode, dim, inum);
        end
    end

    return energy
end

max_stars_repo_name: cesmix-mit/MDP

max_issues_repo_name: cesmix-mit/MDP

max_forks_repo_name: cesmix-mit/MDP


text:
# This file is a part of AstroLib.jl. License is MIT "Expat".
# Copyright (C) 2016 Mosè Giordano.

function aitoff(l::T, b::T) where {T<:AbstractFloat}
    l = rem(l, T(360), RoundNearest)
    alpha2 = deg2rad(l/2)
    delta = deg2rad(b)
    r2 = sqrt(T(2))
    f = 2*r2/pi
    sin_alpha2, cos_alpha2 = sincos(alpha2)
    sin_delta, cos_delta = sincos(delta)
    denom = sqrt(1 + cos_delta * cos_alpha2)
    return rad2deg(cos_delta * sin_alpha2 * 2 * r2 / denom / f),
           rad2deg(sin_delta * r2 / denom / f)
end

"""
    aitoff(l, b) -> x, y

### Purpose ###

Convert longitude `l` and latitude `b` to `(x, y)` using an Aitoff projection.

### Explanation ###

This function can be used to create an all-sky map in Galactic coordinates with
an equal-area Aitoff projection.  Output map coordinates are zero longitude
centered.

### Arguments ###

* `l`: longitude, scalar or vector, in degrees.
* `b`: latitude, number of elements as `l`, in degrees.

Coordinates can be given also as a 2-tuple `(l, b)`.

### Output ###

2-tuple `(x, y)`.

* `x`: x coordinate, same number of elements as `l`.  `x` is normalized to be in
  \$[-180, 180]\$.
* `y`: y coordinate, same number of elements as `l`.  `y` is normalized to be in
  \$[-90, 90]\$.

### Example ###

Get \$(x ,y)\$ Aitoff coordinates of Sirius, whose Galactic coordinates are
\$(227.23, -8.890)\$.

```jldoctest
julia> using AstroLib

julia> x, y = aitoff(227.23, -8.890)
(-137.92196683723276, -11.772527357473054)
```

### Notes ###

See AIPS memo No. 46
(ftp://ftp.aoc.nrao.edu/pub/software/aips/TEXT/PUBL/AIPSMEMO46.PS), page 4, for
details of the algorithm.  This version of `aitoff` assumes the projection is
centered at b=0 degrees.

Code of this function is based on IDL Astronomy User's Library.
"""
aitoff(l::Real, b::Real) = aitoff(promote(float(l), float(b))...)

aitoff(lb::Tuple{Real, Real}) = aitoff(lb...)

function aitoff(l::AbstractArray{L}, b::AbstractArray{B}) where {L<:Real,B<:Real}
    @assert length(l) == length(b)
    typel = float(L)
    x = similar(l, typel)
    y = similar(b, typel)
    for i in eachindex(l)
        x[i], y[i] = aitoff(l[i], b[i])
    end
    return x, y
end


text:
using Genie
Genie.newapp_webservice("MyGenieApp")

# julia --project=@. bootstrap.jl -h
# julia --project=@. bootstrap.jl s
#
# Genie.loadapp()
# up(8000, "127.0.0.1", aync=false)
# up(8000, aync=false)
# startup(8000, "127.0.0.1", aync=false)


text:
# kuramoto_8-7--19.jl

# Plots takes about 8 seconds, and DifferentialEquations
# about 25, to load in.

@time using BenchmarkTools
@time using Plots
@time using DifferentialEquations

plot_dir = "./diary/images/"

# First, let's grab our previous Kuramoto code, and modify
# it so it doesn't wrap around S_1 again.

function KuramotoModel(N, K, ω_nat, wrapS1=false)
    @assert(N==size(ω_nat)[1], "number of oscillators != N")

    function Kuramoto(dθ, θ, p, t)
        for i in 1:N
            if wrapS1
                θ[i] = mod2pi(θ[i])
            end

            stabilizers = map(θ_j -> sin(θ_j - θ[i]), θ)
            stabilizing_term = (K / N) * sum(stabilizers)
            dθ[i] = ω_nat[i] + stabilizing_term

            if wrapS1
                θ[i] = mod2pi(θ[i])
            end
        end
    end

    return Kuramoto
end

# Scheiss, I completely forget how this code worked again.
# Glad I have notes.

tspan = (0.0, 10.0) # Time span.

ω_n = [1.0          # Natural frequencies of the oscs.
       2.0
       5.0
       1.0
       0.5]

K5_no_wrap = KuramotoModel(5, 1.0, ω_n)
θ_init = [0.0
          0.0
          0.0
          0.0
          0.0]

prob = ODEProblem(K5_no_wrap, θ_init, tspan)

@time solution = solve(prob)

plot!(title = "5 Kuramoto modelled oscillators")

plot(solution.t, solution.u)

# Together, these two take about half a minute together.

png(string(plot_dir, "8-7-19_k5-no-wrap_01.png"))

#----------------------------------------------------------#

# Next question: How do I take the solution and wrap it
# around S_1 correctly?

# If we naively modulate it by 2pi, we'll probably lose the
# correct curving on the plot that makes it clear to see
# what's going on. On the other hand, if we modulate inside
# the ODE itself, we're stuck with fixed timesteps.


text:
#`=@BDH("GCNY10Yr Index", "TRADE, BID, ASK", "1/1/2000 00:00:00", "7/7/2020 05:30:00",
#"IntrRw", "True", "Size", "S", "Type", "S","TimeZone=New York",)

PARAM[:geoff] = Dict{Symbol, Any}(
  :tickers => ["GCNY30YR Index", "GCNY10YR Index", "GCNY5YR Index", "GCNY2YR Index", "GCNY1YR Index"],
  :pullfields => ["TRADE", "BID", "ASK"],
  :file => "geoffresults")

function geoffpull(S;
  livepull::Bool="live pull is required",
  tickers::Vector{String}=error("tickers kw are required"),
  fields::Vector{String}=error("fields kw are required"),
  startdatetime = DateTime(2000,1,1,0,0,0),
  enddatetime = DateTime(2020,7,8,0,0,0), #WARNING- time zone not adjusted
  pullname::String=error("pullname is required"),
  outbinstream=OUT_BIN_STREAM,
  binext = BIN_EXTENSION)

  local results
  pullpath::String = "$(PARAM[:inputpath])\\$pullname.$binext"

  #set up a parameter dict

  #this will hold each of the result dictionaries
  #allunique(tickers) || error("Tickers are not unique!")

  #batch the tickers for performance
  Ntickers::Int = length(tickers)

  resultparts::Vector{DataFrame} = Vector{DataFrame}(undef, Ntickers)
  if livepull
    options::Dict{String, Any} = Dict(
      #"Size"=>"S",
      #"Type"=>"S",
      #"TimeZone" => "New York",
      #"DateFormat"=>DateFormat,
      )

    #query BB for each partition
    @sync for j ∈ 1:Ntickers
      @async resultparts[j] = BLPData.bdh_intraday_ticks(S, tickers[j], fields,
        startdatetime, enddatetime,
        options=options) |> DataFrame!
    end

    #label each df
    badtickers::Vector{String} = Vector{String}()
    for j ∈ 1:Ntickers
      if ncol(resultparts[j]) > 1#length(setdiff(fields, names(resultparts[i][j]))) == 0
        resultparts[j][!,:ticker] .= tickers[j]
      else
        push!(badtickers, "$(tickers[j])")
        resultparts[j] = DataFrame()
      end
    end

    @info "Individual dataframes formed. $(length(badtickers)) bad tickers:\n$(badtickers)"

    #combine the results
    @assert length(resultparts) == Ntickers
    vcatunion(df1,df2) = vcat(df1,df2, cols=:union)
    results = reduce(vcatunion, resultparts)

    #write out the results
    outbinstream(pullpath, results)
  else
    results = inbinstream(pullpath)
  end

  #println(results)
  #println(results)
  results
end

function writepull(df::AbstractDataFrame, outname::String;
  compress::Bool = false,#size(df,1) > 150_000, #simple heuristic for compressing the CSV
  outputpath::String = PARAM[:outputpath],
  futuresfile::String = PARAM[:geoff][:file],
  outcsvstream = compress ? OUT_CSV_STREAM : OUT_SMALL_CSV_STREAM,
  csvext = compress ? CSV_EXTENSION : SMALL_CSV_EXTENSION
  )

  p::String = "$outputpath\\$outname.$csvext"
  outcsvstream(p, df)
end

#main entry point
function bbscript(S, ::Val{:geoffpull};
  livepull=error("Live pull is required"))

  fields=PARAM[:geoff][:pullfields]
  tickers=PARAM[:geoff][:tickers]

  results = geoffpull(S,
    livepull=livepull,
    tickers=tickers,
    fields=fields,
    pullname="geoff")

  writepull(results, "$(PARAM[:geoff][:file])")

end

max_stars_repo_name: clintonTE/CCA

max_issues_repo_name: clintonTE/CCA

max_forks_repo_name: clintonTE/CCA


text:
import Pkg; Pkg.add(Pkg.PackageSpec(url="https://github.com/JuliaComputing/JuliaAcademyData.jl"))

using JuliaAcademyData; activate("Deep learning with Flux")

# # Building a single neural network layer using `Flux.jl`
#
# In this notebook, we'll move beyond binary classification. We'll try to distinguish between three fruits now, instead of two. We'll do this using **multiple** neurons arranged in a **single layer**.

#-

# ## Read in and process data

#-

# We can start by loading the necessary packages and getting our data into working order with similar code we used at the beginning of the previous notebooks, except that now we will combine three different apple data sets, and will add in some grapes to the fruit salad!

using CSV, DataFrames, Flux, Plots

#-

## Load apple data in CSV.read for each file
apples1 = DataFrame(CSV.File(datapath("data/Apple_Golden_1.dat"), delim='\t', allowmissing=:none, normalizenames=true))
apples2 = DataFrame(CSV.File(datapath("data/Apple_Golden_2.dat"), delim='\t', allowmissing=:none, normalizenames=true))
apples3 = DataFrame(CSV.File(datapath("data/Apple_Golden_3.dat"), delim='\t', allowmissing=:none, normalizenames=true))
## And then concatenate them all together
apples = vcat(apples1, apples2, apples3)

# And now let's build an array called `x_apples` that stores data from the `red` and `blue` columns of `apples`. From `applecolnames1`, we can see that these are the 3rd and 5th columns of `apples`:

size(apples)

#-

apples[1:1, :]

#-

x_apples  = [ [apples[i, :red], apples[i, :blue]] for i in 1:size(apples, 1) ]

# Similarly, let's create arrays called `x_bananas` and `x_grapes`:

## Load data from *.dat files
bananas = DataFrame(CSV.File(datapath("data/Banana.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes1 = DataFrame(CSV.File(datapath("data/Grape_White.dat"), delim='\t', allowmissing=:none, normalizenames=true))
grapes2 = DataFrame(CSV.File(datapath("data/Grape_White_2.dat"), delim='\t', allowmissing=:none, normalizenames=true))

## Concatenate data from two grape files together
grapes = vcat(grapes1, grapes2)

## Build x_bananas and x_grapes from bananas and grapes
x_bananas  = [ [bananas[i, :red], bananas[i, :blue]] for i in 1:size(bananas, 1) ]
x_grapes = [ [grapes[i, :red], grapes[i, :blue]] for i in 1:size(grapes, 1) ]

#-

xs = vcat(x_apples, x_bananas, x_grapes)

# ## One-hot vectors

#-

# Now we wish to classify *three* different types of fruit. It is not clear how to encode these three types using a single output variable; indeed, in general this is not possible.
#
# Instead, we have the idea of encoding $n$ output types from the classification into *vectors of length $n$*, called "one-hot vectors":
#
# $$
# \textrm{apple} = \begin{pmatrix} 1 \\ 0 \\ 0 \end{pmatrix};
# \quad
# \textrm{banana} = \begin{pmatrix} 0 \\ 1 \\ 0 \end{pmatrix};
# \quad
# \textrm{grape} = \begin{pmatrix} 0 \\ 0 \\ 1 \end{pmatrix}.
# $$
#
# The term "one-hot" refers to the fact that each vector has a single $1$, and is $0$ otherwise.
#
# Effectively, the first neuron will learn whether or not (1 or 0) the data corresponds to an apple, the second whether or not (1 or 0) it corresponds to a banana, etc.

#-

# `Flux.jl` provides an efficient representation for one-hot vectors, using advanced features of Julia so that it does not actually store these vectors, which would be a waste of memory; instead `Flux` just records in which position the non-zero element is. To us, however, it looks like all the information is being stored:

using Flux: onehot

onehot(1, 1:3)

# #### Exercise 1
#
# Make an array `labels` that gives the labels (1, 2 or 3) of each data point. Then use `onehot` to encode the information about the labels as a vector of `OneHotVector`s.

#-

# **Solution**:

labels = [ones(length(x_apples)); 2*ones(length(x_bananas)); 3*ones(length(x_grapes))]

#-

ys = [onehot(label, 1:3) for label in labels]  # onehotbatch(labels, 1:3)

# The input data is in `xs` and the one-hot vectors are in `ys`.

#-

# ## Single layer in Flux

#-

# Let's suppose that there are two pieces of input data, as in the previous single neuron notebook. Then the network has 2 inputs and 3 outputs:

include(datapath("scripts/draw_neural_net.jl"))
draw_network([2, 3])

# `Flux` allows us to express this again in a simple way:

using Flux

#-

model = Dense(2, 3, σ)

# #### Exercise 2
#
# Now what do the weights inside `model` look like? How does this compare to the diagram of the network layer above?

#-

# #### Solution

model.W

# Each of the 6 lines in the figure denotes a weight of the neuron on the right, taking as input the output of the neuron on the left. These weights are collected in the **matrix** `W`. Note that it seems to be "backwards", since it is designed to multiply vectors of length 2 (the input size):

x = rand(2)
model.W * x

# The whole `model` object represents a set of three sigmoidal neurons. Again, we can use `model` as a function and apply it to input data to get the model's prediction from that data:

model(x)

# This is equivalent to the following calculation:

σ.( (model.W) * x .+ (model.b) )

# Note that here we have again used Julia's **broadcasting** capability, in which the function $\sigma$ is applied to each element of the vector `W * x` in turn. This elementwise application of the function is implicit (i.e. not written out) in most of the literature on machine learning, but it is very much clearer to make this explicit, as Julia forces us to do.
#
# [Recall that we also use `.` to access the fields `W` and `b` in the `model` object; this is a different use of `.`.]

#-

# ## Training the model

#-

# Despite the fact that the model is now more complicated than the single neuron from the previous notebook, the beauty of `Flux.jl` is that the rest of the training process **looks exactly the same**!

#-

# #### Exercise 3
#
# Implement training for this model.

#-

# **Solution**:

loss(x, y) = Flux.mse(model(x), y)

#-

data = zip(xs, ys);

#-

opt = SGD(params(model), 0.01)
## give a list of the parameters that will be modified

#-

for i in 1:100
    Flux.train!(loss, data, opt)
end

#-

model.W

#-

model.b

# #### Exercise 4
#
# Visualize the result of the learning for each neuron. Since each neuron is sigmoidal, we can get a good idea of the function by just plotting a single contour level where the function takes the value 0.5, using the `contour` function with keyword argument `levels=[0.5, 0.501]`.

using Plots
plot()

contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[1], levels=[0.5, 0.501], color = cgrad([:blue, :blue]))
contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[2], levels=[0.5,0.501], color = cgrad([:green, :green]))
contour!(0:0.01:1, 0:0.01:1, (x,y)->model([x,y]).data[3], levels=[0.5,0.501], color = cgrad([:red, :red]))

scatter!(first.(x_apples), last.(x_apples), m=:cross, label="apples")
scatter!(first.(x_bananas), last.(x_bananas), m=:circle, label="bananas")
scatter!(first.(x_grapes), last.(x_grapes), m=:square, label="grapes")

# #### Exercise 5
#
# Interpret the results by checking which fruit each neuron was supposed to learn and what it managed to achieve.

#-

# #### Solution

#-

# We see that two of the hyperplanes have been learnt correctly, the one that separates bananas from the rest, and the one that separates grapes from the rest. However, the hyperplane for the neuron that was supposed to learn to distinguish apples is not at all correct.
#
# It is hopefully clear why this is so: there *is no way* to separate apples from non-apples with a *single* hyperplane, given this data.
#
# The problem is that this network is so simple that it is *not capable of learning a good representation of the data*: the class of functions that are modelled by this single layer, with only 9 parameters, is not complex enough.
#
# In order to learn how to classify three different types of fruit, we must thus move to more complex neural networks, and towards **deep learning**!

max_stars_repo_name: ricardoprins/manbearpig

max_issues_repo_name: ricardoprins/manbearpig

max_forks_repo_name: ricardoprins/manbearpig


text:
module Stingray

#  Packages/core modules, Ex - using Dates

end


text:
import CurveFit
using Printf
using stoosc
rc("font", family="arial")

MODEL = "model02b"  # model02 or model02b

# Maximal considered input intensity
I_max = 3.0

## Deterministic Arnold tongue
filename = joinpath("arnold_tongues", MODEL, "ode.csv")
arnold_tongue = load_data(filename)
entrainment_area_ode = estimate_entrainment_area(arnold_tongue, I_max)
println("ode entrainment area = $(entrainment_area_ode)")
fig2, ax2 = subplots(figsize=(2.5, 2.25))
plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
ax2.set_title("Deterministic")
fig2.tight_layout()
filename = joinpath("arnold_tongues", MODEL, "arnold_deterministic.svg")
save_figure(fig2, filename)

## Arnold tongue for a stochastic model with many cells
arnold_tongue = load_data("arnold_tongues", MODEL, "sde_n=1000_σ=0.2.csv")
entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
println("entrainment area = $(entrainment_area)")
fig2, ax2 = subplots(figsize=(2.5, 2.25))
plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
ax2.set_title("n = 1000, σ = 0.2")
fig2.tight_layout()
filename = joinpath("arnold_tongues", MODEL, "arnold_many_cells.svg")
save_figure(fig2, filename)

## Arnold tongue for a stochastic model with few cells
arnold_tongue = load_data("arnold_tongues", MODEL, "sde_n=10_σ=0.2.csv")
entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
println("entrainment area = $(entrainment_area)")
fig2, ax2 = subplots(figsize=(2.5, 2.25))
plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
ax2.set_title("n = 10, σ = 0.2")
fig2.tight_layout()
filename = joinpath("arnold_tongues", MODEL, "arnold_few_cells.svg")
save_figure(fig2, filename)

## Increasing noise for constant population size
fig, ax = subplots(figsize=(7, 3))
ax.plot(log10.([0.01, 1.0]), [entrainment_area_ode, entrainment_area_ode], "--",
    color="black", label="Deterministic")

# Save the values for the maximal entrainment area
maximum_entrainment_area = Dict(
    "n" => Float64[],
    "σ" => Float64[],
    "entrainment_area" => Float64[]
)

x = [0.01, 0.02, 0.05, 0.1, 0.2, 0.5, 1.0]  # noise intensitiess
for n = [1, 10, 100, 1000]  # number of cells
    # for each population size

    y = Float64[]  # entrainment area

    for σ = x
        # for each noise value

        # Load the arnold tongue for given n and σ
        filename = joinpath("arnold_tongues", MODEL, "sde_n=$(n)_σ=$(σ).csv")
        arnold_tongue = load_data(filename)

        # Save arnold tongues for the individual points in the figure
        fig2, ax2 = subplots(figsize=(2.5, 2.25))
        plot_arnold_tongue(arnold_tongue, I_max, ax=ax2)
        ax2.set_title("n = $(n); σ = $(σ)")
        fig2.tight_layout()
        filename = joinpath("arnold_tongues", MODEL, "sde_n_$(n)_σ_$(σ).svg")
        save_figure(fig2, filename)

        # Calculate and save the entrainment area
        entrainment_area = estimate_entrainment_area(arnold_tongue, I_max)
        push!(y, entrainment_area)

    end

    # Plot entrainment area versus noise intensity
    x_log = log10.(x)
    if MODEL == "model02"
        f_fit = CurveFit.curve_fit(CurveFit.Polynomial, x_log, y, 2)
    elseif MODEL == "model02b"
        f_fit = CurveFit.curve_fit(CurveFit.Polynomial, x_log, y, 3)
    end
    h = ax.plot(x_log, y, "o", label="n = $n")

    # Fit a polynomial for a smooth curve interpolation and plot it
    x_new = range(x_log[1], x_log[end], length=100)
    y_new = f_fit.(x_new)
    ax.plot(x_new, y_new, "-", color=h[1].get_color())

    # Find point of maximal entrainment area and save it into a dictionary
    index = argmax(y_new)
    push!(maximum_entrainment_area["σ"], x_new[index])
    push!(maximum_entrainment_area["n"], n)
    push!(maximum_entrainment_area["entrainment_area"], y_new[index])
    # ax.plot(x_new[index], y_new[index], "x", markersize=10, color=h[1].get_color())
    println("Maximum for n = $n is at σ ≈ $(10^x_new[index])")

end

# Set labels for the figure
ax.set_xlabel("log\$_{10}\$(Noise intensity)")
ax.set_ylabel("Entrainment area")
ax.set_title("Limit cycle oscillator")
ax.legend(edgecolor="black", framealpha=1.0, loc="best", ncol=2)
fig.tight_layout()

# Save the main figure
filename = joinpath("arnold_tongues", MODEL, "overview_n.svg")
save_figure(fig, filename)

# ## Plot maximal entrainment area and corresponding noise intenisty for each n
# fig, ax_arr = subplots(1, 2, figsize=(7, 3))
# n = log10.(maximum_entrainment_area["n"])
# σ = maximum_entrainment_area["σ"]
# entrainment_area = maximum_entrainment_area["entrainment_area"]

# ax_arr[1].plot(n, σ, "o", color="black")
# ax_arr[1].set_xlabel("log\$_{10}\$(Number of cells, n)")
# ax_arr[1].set_ylabel("log\$_{10}\$(Noise intensity, σ)")
# ax_arr[1].set_title("Noise intensity for maximal entrainment area")

# ax_arr[2].plot(n, entrainment_area, "o", color="black")
# ax_arr[2].set_xlabel("log\$_{10}\$(Number of cells, n)")
# ax_arr[2].set_ylabel("Entrainment area")
# ax_arr[2].set_title("Maximal entrainment area")

# fig.tight_layout()
# save_figure(fig, "arnold_tongues", MODEL, "maximal_entrainment_n.svg")

max_stars_repo_name: vkumpost/stoosc

max_issues_repo_name: vkumpost/stoosc

max_forks_repo_name: vkumpost/stoosc


text:
module ElectronDisplay

using Electron, Base64

struct ElectronDisplayType <: Base.AbstractDisplay end

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("image/png")}, x)
    img = stringmime(MIME("image/png"), x)

    payload = string("")

    w = Electron.Window(payload)
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("image/png")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("image/svg+xml")}, x)
    payload = stringmime(MIME("image/svg+xml"), x)

    w = Electron.Window(payload)
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("image/svg+xml")}) = true

asset(url...) = replace(normpath(joinpath(@__DIR__, "..", "assets", url...)), "\\" => "/")

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}, x)
    payload = stringmime(MIME("application/vnd.vegalite.v2+json"), x)

    html_page = """

    """

    w = Electron.Window(html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vegalite.v2+json")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v3+json")}, x)
    payload = stringmime(MIME("application/vnd.vega.v3+json"), x)

    html_page = """

    """

    w = Electron.Window(html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.vega.v3+json")}) = true

function Base.display(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.plotly.v1+json")}, x)
    payload = stringmime(MIME("application/vnd.plotly.v1+json"), x)

    html_page = """

    """

    w = Electron.Window(html_page, options=Dict("webPreferences" => Dict("webSecurity" => false)))
end

Base.displayable(d::ElectronDisplayType, ::MIME{Symbol("application/vnd.plotly.v1+json")}) = true

function Base.display(d::ElectronDisplayType, x)
    if showable("application/vnd.vegalite.v2+json", x)
        display(d,MIME("application/vnd.vegalite.v2+json"), x)
    elseif showable("application/vnd.vega.v3+json", x)
            display(d,MIME("application/vnd.vega.v3+json"), x)
    elseif showable("application/vnd.plotly.v1+json", x)
            display(d,MIME("application/vnd.plotly.v1+json"), x)
    elseif showable("image/svg+xml", x)
        display(d,"image/svg+xml", x)
    elseif showable("image/png", x)
        display(d,"image/png", x)
    else
        throw(MethodError(Base.display,(d,x)))
    end
end

function __init__()
    Base.Multimedia.pushdisplay(ElectronDisplayType())
end

end # module


text:
##### Beginning of file

function _print_welcome_message()::Nothing
    organizationsnapshots_version::VersionNumber = version()
    organizationsnapshots_pkgdir::String = package_directory()
    @info(string("This is Snapshots, version ",organizationsnapshots_version,),)
    @info(string("For help, please visit https://github.com/UnofficialJuliaMirrorSnapshots/Snapshots.jl",),)
    @debug(string("Snapshots package directory: ",organizationsnapshots_pkgdir,),)
    return nothing
end

##### End of file


text:
@symbol_func function K_VT(cur_reactor::AbstractReactor)
  cur_nu_n = cur_reactor.nu_n

  cur_nu_T = cur_reactor.nu_T

  cur_K = K_n(cur_reactor)

  cur_K *= cur_reactor.epsilon ^ 2

  cur_K *= g_p(cur_reactor)

  cur_K *= ( 1 + cur_reactor.f_D )

  cur_K *= ( 1 + cur_nu_n )

  cur_K *= ( 1 + cur_nu_T )

  cur_K /= 1 + cur_nu_n + cur_nu_T

  cur_K
end


text:
# using Distributed
# addprocs(5)

using Caesar, RoME, Profile
using RoMEPlotting, Gadfly, KernelDensityEstimatePlotting
using KernelDensityEstimate
using IncrementalInference
using DocStringExtensions
using DelimitedFiles
# using TransformUtils

include(joinpath(@__DIR__,"slamUtils.jl"))

## Default parameters
N = 100
pμ = [0.0,0,1,0]
pσ = Matrix(Diagonal([0.1;0.1;0.1;0.1].^2))

symbolstart = 1;
windowstart = 70;
saswindow = 5;

poses = [Symbol("x$i") for i in symbolstart:(symbolstart+saswindow-1)]
sasframes = collect(windowstart:1:(windowstart+saswindow-1))

# dataDir = joinpath(ENV["HOME"],"data","sas","06_20_sample");
dataDir = joinpath(ENV["HOME"],"data","kayaks","20_gps_pos");
posData = importdata_nav(sasframes, datadir=dataDir);
navchecked, errorind = sanitycheck_nav(posData)

cfgFile = joinpath(ENV["HOME"],"data","sas","SAS2D.yaml");
cfgd=loadConfigFile(cfgFile)
chirpFile = joinpath(ENV["HOME"],"data","sas","chirp250.txt");

## Build the factor graph
fg = initfg();

beacon = :l1
addVariable!(fg, beacon, Point2 )

dposData = deepcopy(posData)
cumulativeDrift!(dposData,[0.0;0],[0.1,0.1])
waveformData = importdata_waveforms(sasframes,2, datadir=dataDir);
tcurrent = windowstart*1_000_000

for sym in poses
  global tcurrent
  addVariable!(fg, sym, DynPoint2(ut=tcurrent))
  tcurrent += 1_000_000
end

priors = [1]
#Priors
for i in priors
    xdotp = posData[i+1,1] - posData[i,1];
    ydotp = posData[i+1,2] - posData[i,2];
    dpμ = [posData[i,1];posData[i,2];xdotp;ydotp];
    dpσ = Matrix(Diagonal([0.1;0.1;0.1;0.1].^2))
    pp = DynPoint2VelocityPrior(MvNormal(dpμ,dpσ))
    addFactor!(fg, [poses[i];], pp, autoinit=false)
end

vps = [1,2,3,4]
for i in vps
    xdotp = dposData[i+1,1] - dposData[i,1];
    ydotp = dposData[i+1,2] - dposData[i,2];
    dpμ = [xdotp;ydotp;0;0];
    dpσ = Matrix(Diagonal([0.1;0.1;0.1;0.1].^2))
    vp = VelPoint2VelPoint2(MvNormal(dpμ,dpσ))
    addFactor!(fg, [poses[i];poses[i+1]], vp, autoinit=false)
    # IncrementalInference.doautoinit!(fg,[getVariable(fg,poses[i])])
end

# IncrementalInference.doautoinit!(fg,[getVariable(fg,poses[5])])
sas2d = prepareSAS2DFactor(saswindow, waveformData, rangemodel=:Correlator,
                           cfgd=cfgd, chirpFile=chirpFile)
addFactor!(fg, [beacon;poses], sas2d, autoinit=false)

writeGraphPdf(fg, engine="dot")

getSolverParams(fg).drawtree = true
getSolverParams(fg).showtree = true
# getSolverParams(fg).downsolve

## solve the factor graph
tree, smt, hist = solveTree!(fg, recordcliqs=ls(fg))

printCliqHistorySummary(tree,:x5)
drawTree(tree)

Profile.init(n = 10^7, delay = 0.01)
@profile stuff = sandboxCliqResolveStep(tree,:x3,6)
ProfileView.view()
Profile.print(format=:flat)

## Plot the SAS factor pairs
fsym = :l1x1x2x3x4x5f1
pl = plotSASPair(fg, fsym, show=true, filepath="/tmp/testDP_2pp_3vpvp_init.pdf");

plk=[]
X1 = getKDEMean(getBelief(fg,:x1))
push!(plk, layer(x=X1[1,:],y=X1[2,:], Geom.point))
tplt = Gadfly.plot(plk...); tplt |> SVG("/tmp/test.svg") ; tplt |> PDF("/tmp/test.pdf")

plotKDE(fg, ls(fg,r"x"), dims=[1;2])

## More debugging below

stuff = IIF.localProduct(fg, :x1)
pl = plotKDE(stuff[1], dims=[1;2], levels=3, c=["blue"])


text:
using Sockets

function matgen(n, seed)
    tmp = seed / n / n
    [tmp * (i - j) * (i + j - 2) for i = 1:n, j = 1:n]
end

function mul(a, b)
    m, n = size(a)
    q, p = size(b)
    @assert n == q

    # transpose a for cache-friendliness
    aT = zeros(n, m)
    @simd for i = 1:m
        for j = 1:n
            @inbounds aT[j, i] = a[i, j]
        end
    end

    out = zeros(m, p)
    @simd for i = 1:m
        for j = 1:p
            z = 0.0
            for k = 1:n
                @inbounds z += aT[k, i] * b[k, j]
            end
            @inbounds out[i, j] = z
        end
    end
    out
end

function calc(n)
    n = round(Int, n / 2) * 2
    a = matgen(n, 1.0)
    b = matgen(n, 2.0)
    c = mul(a, b)
    c[Int(n / 2)+1, Int(n / 2)+1]
end

function notify(msg)
    try
        socket = connect("localhost", 9001)
        write(socket, msg)
        close(socket)
    catch
        # standalone usage
    end
end

if abspath(PROGRAM_FILE) == @__FILE__
    n = length(ARGS) > 0 ? parse(Int, ARGS[1]) : 100

    left = calc(101)
    right = -18.67
    if abs(left - right) > 0.1
        println(stderr, "$(left) != $(right)")
        exit(1)
    end

    notify("Julia (no BLAS)\t$(getpid())")
    t = time()
    results = calc(n)
    elapsed = time() - t
    notify("stop")

    println(results)
    println("time: $(elapsed) s")
end


text:
#Run sample code
#-------------------------------------------------------------------------------

#==Show results
===============================================================================#

for i in 1:14
	file = "./demo$i.jl"
	sepline = "---------------------------------------------------------------------"
	println("\nExecuting $file...")
	println(sepline)
	evalfile(file)
end

:SampleCode_Executed


text:
function get_best_bound(m::Model)
    objval_p = Array(Cdouble, 1)
    stat = @cpx_ccall(getbestobjval, Cint, (Ptr{Void}, Ptr{Void}, Ptr{Cdouble}), m.env.ptr, m.lp, objval_p)
    if stat != 0
        throw(CplexError(m.env.ptr, stat))
    end
    return objval_p[1]
end

type CallbackData
    cbdata::Ptr{Void}
    model::Model
end

export CallbackData

function setcallbackcut(cbdata::CallbackData, where::Cint, ind::Vector{Cint}, val::Vector{Cdouble}, sense::Char, rhs::Cdouble, purgeable::Cint)
    len = length(ind)
    @assert length(val) == len
    sns = convert(Cint, sense)
    stat = @cpx_ccall(cutcallbackadd, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Cint,
                      Cdouble,
                      Cint,
                      Ptr{Cint},
                      Ptr{Cdouble},
                      Cint
                      ),
                      cbdata.model.env.ptr, cbdata.cbdata, where, len, rhs, sns, ind.-1, val, purgeable)
    if stat != 0
        throw(CplexError(cbdata.model.env.ptr, stat))
    end
end

function setcallbackcutlocal(cbdata::CallbackData, where::Cint, ind::Vector{Cint}, val::Vector{Cdouble}, sense::Char, rhs::Cdouble, purgeable::Cint)
    len = length(ind)
    @assert length(val) == len
    sns = convert(Cint, sense)
    stat = @cpx_ccall(cutcallbackaddlocal, Cint, (
                      Ptr{Void},
                      Ptr{Void},
                      Cint,
                      Cint,
                      Cdouble,
                      Cint,
                      Ptr{Cint},
                      Ptr{Cdouble},
                      Cint
                      ),
                      cbdata.model.env.ptr, cbdata.cbdata, where, len, rhs, sns, ind.-1, val, purgeable)
    if stat != 0
        throw(CplexError(cbdata.model.env.ptr, stat))
    end
end

cbcut(cbdata::CallbackData, where::Cint, ind::Vector{Cint}, val::Vector{Cdouble}, sense::Char, rhs::Cdouble) =
        setcallbackcut(cbdata, where, ind, val, sense, rhs, convert(Cint,CPX_USECUT_PURGE))

cbcutlocal(cbdata::CallbackData, where::Cint, ind::Vector{Cint}, val::Vector{Cdouble}, sense::Char, rhs::Cdouble) =
        setcallbackcutlocal(cbdata, where, ind, val, sense, rhs, convert(Cint,CPX_USECUT_PURGE))

cblazy(cbdata::CallbackData, where::Cint, ind::Vector{Cint}, val::Vector{Cdouble}, sense::Char, rhs::Cdouble) =
        setcallbackcut(cbdata, where, ind, val, sense, rhs, convert(Cint,CPX_USECUT_FORCE))

cblazylocal(cbdata::CallbackData, where::Cint, ind::Vector{Cint}, val::Vector{Cdouble}, sense::Char, rhs::Cdouble) =
        setcallbackcutlocal(cbdata, where, ind, val, sense, rhs, convert(Cint,CPX_USECUT_FORCE))

function cbbranch(cbdata::CallbackData, where::Cint, idx::Cint, LU::Cchar, bd::Cdouble, nodeest::Cdouble)
    seqnum = Array(Cint,1)
    stat = @cpx_ccall(branchcallbackbranchbds, Cint, (Ptr{Void},Ptr{Void},Cint,Cint,Ptr{Cint},Ptr{Cchar},Ptr{Cdouble},Cdouble,Ptr{Void},Ptr{Cint}),
                      cbdata.model.env.ptr,cbdata.cbdata,where,convert(Cint,1),[idx],[LU],[bd],nodeest,C_NULL,seqnum)
    if stat != 0
        throw(CplexError(cbdata.model.env.ptr, stat))
    end
    return seqnum[1]
end

function cbbranchconstr(cbdata::CallbackData, where::Cint, indices::Vector{Cint}, coeffs::Vector{Cdouble}, rhs::Cdouble, sense::Cchar, nodeest::Cdouble)
    seqnum = Array(Cint,1)
    stat = @cpx_ccall(branchcallbackbranchconstraints, Cint,
                      (Ptr{Void},
                       Ptr{Void},
                       Cint,
                       Cint,
                       Cint,
                       Ptr{Cdouble},
                       Ptr{Cchar},
                       Ptr{Cint},
                       Ptr{Cint},
                       Ptr{Cdouble},
                       Cdouble,
                       Ptr{Void},
                       Ptr{Cint}),
                      cbdata.model.env.ptr,
                      cbdata.cbdata,
                      where,
                      convert(Cint,1),
                      convert(Cint,length(indices)),
                      [rhs],
                      [sense],
                      Cint[0],
                      indices,
                      coeffs,
                      nodeest,
                      C_NULL,
                      seqnum)
    if stat != 0
        throw(CplexError(cbdata.model.env.ptr, stat))
    end
    return seqnum[1]
end


text:
function metadata(ds::Dataset)
    # dimensions

    dim = OrderedDict()
    unlimited_dims = unlimited(ds.dim)

    for (dimname,dimlen) in ds.dim
        dim[dimname] = Dict(
            "name" => dimname,
            "length" => dimlen,
            "unlimited" => dimname in unlimited_dims
        )
    end

    # variables
    vars = OrderedDict()

    for (varname,ncvar) in ds
        storage,chunksizes = chunking(ncvar.var)
        isshuffled,isdeflated,deflatelevel = deflate(ncvar.var)

        vars[varname] = OrderedDict(
            "name" => varname,
            "size" => size(ncvar),
            "eltype" => eltype(ncvar.var),
            "attrib" => OrderedDict(ncvar.attrib),
            "dimensions" => dimnames(ncvar),
            "chunksize" => chunksizes,
            "fillvalue" => fillvalue(ncvar.var),
            "shuffle" => isshuffled,
            "deflatelevel" => deflatelevel
        )
    end

    group = OrderedDict()
    for (groupname,ncgroup) in ds.group
        group[groupname] = metadata(ncgroup)
    end

    return OrderedDict(
        "dim" => dim,
        "var" => vars,
        "attrib" => OrderedDict(ds.attrib),
        "group" => group
        )
end

function DeferDataset(filename,mode,info)
    Dataset(filename,mode) do ds
        r = Resource(filename,mode,info)
        groupname = "/"
        da = DeferAttributes(r,"/",r.metadata["attrib"])
        dd = DeferDimensions(r,r.metadata["dim"])
        dg = DeferGroups(r,r.metadata["group"])
        return DeferDataset(r,groupname,da,dd,dg,info)
    end
end

function DeferDataset(filename,mode = "r")
    Dataset(filename,mode) do ds
        info = metadata(ds)
        r = Resource(filename,mode,info)
        groupname = "/"
        da = DeferAttributes(r,"/",r.metadata["attrib"])
        dd = DeferDimensions(r,r.metadata["dim"])
        dg = DeferGroups(r,r.metadata["group"])
        return DeferDataset(r,groupname,da,dd,dg,info)
    end
end

close(dds::DeferDataset) = nothing
groupname(dds::DeferDataset) = dds.groupname
path(dds::DeferDataset) = dds.r.filename
Base.keys(dds::DeferDataset) = collect(keys(dds.data["var"]))

function Dataset(f::Function, r::Resource)
    Dataset(r.filename,r.mode) do ds
        f(ds)
    end
end

function Dataset(f::Function, dds::DeferDataset)
    Dataset(dds.r.filename,dds.r.mode) do ds
        f(ds)
    end
end

function Variable(f::Function, dv::DeferVariable)
    Dataset(dv.r.filename,dv.r.mode) do ds
        f(variable(ds,dv.varname))
    end
end

function variable(dds::DeferDataset,varname::AbstractString)
    data = dds.data["var"][varname]
    T = data["eltype"]
    N = length(data["dimensions"])
    da = DeferAttributes(dds.r,varname,data["attrib"])

    return DeferVariable{T,N}(dds.r,varname,da,data)
end

function Base.getindex(dv::DeferVariable,indexes::Union{Int,Colon,UnitRange{Int},StepRange{Int,Int}}...)
    Variable(dv) do v
        return v[indexes...]
    end
end

Base.size(dv::DeferVariable) = dv.data["size"]
dimnames(dv::DeferVariable) = dv.data["dimensions"]
name(dv::DeferVariable) = dv.varname

#----------------------------------------------
Base.keys(dd::DeferDimensions) = collect(keys(dd.data))
Base.getindex(dd::DeferDimensions,name::AbstractString) = dd.data[name]["length"]

Base.keys(da::DeferAttributes) = collect(keys(da.data))
Base.getindex(da::DeferAttributes,name::AbstractString) = da.data[name]

#------------------------------------------------

Base.keys(dg::DeferGroups) = collect(keys(dg.data))
function Base.getindex(dg::DeferGroups,name::AbstractString)
    data = dg.data[name]
    da = DeferAttributes(dg.r,"/",data["attrib"])
    dd = DeferDimensions(dg.r,data["dim"])
    dg = DeferGroups(dg.r,data["group"])
    return DeferDataset(dg.r,name,da,dd,dg,data)
end


text:
@testset "4.4 Cotangent" begin
    include("4.4.0 (a trg)^m (b cot)^n.jl")
    include("4.4.1.2 (d csc)^m (a+b cot)^n.jl")
    include("4.4.1.3 (d cos)^m (a+b cot)^n.jl")
    include("4.4.10 (c+d x)^m (a+b cot)^n.jl")
    include("4.4.2.1 (a+b cot)^m (c+d cot)^n.jl")
    include("4.4.7 (d trig)^m (a+b (c cot)^n)^p.jl")
    include("4.4.9 trig^m (a+b cot^n+c cot^(2 n))^p.jl")
end

max_stars_repo_path: test/4 Trig functions/4.4 Cotangent/.jl

max_issues_repo_path: test/4 Trig functions/4.4 Cotangent/.jl

max_forks_repo_path: test/4 Trig functions/4.4 Cotangent/.jl


text:
module DataAssim

export timecorr_noise
export resample
export enkf

include("utils_pfilter.jl")
include("utils_enkf.jl")

end


text:
export bicgstabl, bicgstabl!, bicgstabl_iterator, bicgstabl_iterator!, BiCGStabIterable

import Base: start, next, done

mutable struct BiCGStabIterable{precT, matT, solT, vecT <: AbstractVector, smallMatT <: AbstractMatrix, realT <: Real, scalarT <: Number}
    A::matT
    l::Int

    x::solT
    r_shadow::vecT
    rs::smallMatT
    us::smallMatT

    max_mv_products::Int
    mv_products::Int
    reltol::realT
    residual::realT

    Pl::precT

    γ::vecT
    ω::scalarT
    σ::scalarT
    M::smallMatT
end

function bicgstabl_iterator!(x, A, b, l::Int = 2;
    Pl = Identity(),
    max_mv_products = size(A, 2),
    initial_zero = false,
    tol = sqrt(eps(real(eltype(b))))
)
    T = eltype(x)
    n = size(A, 1)
    mv_products = 0

    # Large vectors.
    r_shadow = rand(T, n)
    rs = Matrix{T}(n, l + 1)
    us = zeros(T, n, l + 1)

    residual = view(rs, :, 1)

    # Compute the initial residual rs[:, 1] = b - A * x
    # Avoid computing A * 0.
    if initial_zero
        copy!(residual, b)
    else
        A_mul_B!(residual, A, x)
        residual .= b .- residual
        mv_products += 1
    end

    # Apply the left preconditioner
    A_ldiv_B!(Pl, residual)

    γ = zeros(T, l)
    ω = σ = one(T)

    nrm = norm(residual)

    # For the least-squares problem
    M = zeros(T, l + 1, l + 1)

    # Stopping condition based on relative tolerance.
    reltol = nrm * tol

    BiCGStabIterable(A, l, x, r_shadow, rs, us,
        max_mv_products, mv_products, reltol, nrm,
        Pl,
        γ, ω, σ, M
    )
end

@inline converged(it::BiCGStabIterable) = it.residual ≤ it.reltol
@inline start(::BiCGStabIterable) = 0
@inline done(it::BiCGStabIterable, iteration::Int) = it.mv_products ≥ it.max_mv_products || converged(it)

function next(it::BiCGStabIterable, iteration::Int)
    T = eltype(it.x)
    L = 2 : it.l + 1

    it.σ = -it.ω * it.σ

    ## BiCG part
    for j = 1 : it.l
        ρ = dot(it.r_shadow, view(it.rs, :, j))
        β = ρ / it.σ

        # us[:, 1 : j] .= rs[:, 1 : j] - β * us[:, 1 : j]
        it.us[:, 1 : j] .= it.rs[:, 1 : j] .- β .* it.us[:, 1 : j]

        # us[:, j + 1] = Pl \ (A * us[:, j])
        next_u = view(it.us, :, j + 1)
        A_mul_B!(next_u, it.A, view(it.us, :, j))
        A_ldiv_B!(it.Pl, next_u)

        it.σ = dot(it.r_shadow, next_u)
        α = ρ / it.σ

        it.rs[:, 1 : j] .-= α .* it.us[:, 2 : j + 1]

        # rs[:, j + 1] = Pl \ (A * rs[:, j])
        next_r = view(it.rs, :, j + 1)
        A_mul_B!(next_r, it.A , view(it.rs, :, j))
        A_ldiv_B!(it.Pl, next_r)

        # x = x + α * us[:, 1]
        it.x .+= α .* view(it.us, :, 1)
    end

    # Bookkeeping
    it.mv_products += 2 * it.l

    ## MR part

    # M = rs' * rs
    Ac_mul_B!(it.M, it.rs, it.rs)

    # γ = M[L, L] \ M[L, 1]
    F = lufact!(view(it.M, L, L))
    A_ldiv_B!(it.γ, F, view(it.M, L, 1))

    # This could even be BLAS 3 when combined.
    BLAS.gemv!('N', -one(T), view(it.us, :, L), it.γ, one(T), view(it.us, :, 1))
    BLAS.gemv!('N', one(T), view(it.rs, :, 1 : it.l), it.γ, one(T), it.x)
    BLAS.gemv!('N', -one(T), view(it.rs, :, L), it.γ, one(T), view(it.rs, :, 1))

    it.ω = it.γ[it.l]
    it.residual = norm(view(it.rs, :, 1))

    it.residual, iteration + 1
end

# Classical API

"""
    bicgstabl(A, b, l; kwargs...) -> x, [history]

Same as [`bicgstabl!`](@ref), but allocates a solution vector `x` initialized with zeros.
"""
bicgstabl(A, b, l = 2; kwargs...) = bicgstabl!(zerox(A, b), A, b, l; initial_zero = true, kwargs...)

"""
    bicgstabl!(x, A, b, l; kwargs...) -> x, [history]

# Arguments

- `A`: linear operator;
- `b`: right hand side (vector);
- `l::Int = 2`: Number of GMRES steps.

## Keywords

- `max_mv_products::Int = size(A, 2)`: maximum number of matrix vector products.
For BiCGStab(l) this is a less dubious term than "number of iterations";
- `Pl = Identity()`: left preconditioner of the method;
- `tol::Real = sqrt(eps(real(eltype(b))))`: tolerance for stopping condition `|r_k| / |r_0| ≤ tol`.
   Note that (1) the true residual norm is never computed during the iterations,
   only an approximation; and (2) if a preconditioner is given, the stopping condition is based on the
   *preconditioned residual*.

# Return values

**if `log` is `false`**

- `x`: approximate solution.

**if `log` is `true`**

- `x`: approximate solution;
- `history`: convergence history.
"""
function bicgstabl!(x, A, b, l = 2;
    tol = sqrt(eps(real(eltype(b)))),
    max_mv_products::Int = size(A, 2),
    log::Bool = false,
    verbose::Bool = false,
    Pl = Identity(),
    kwargs...
)
    history = ConvergenceHistory(partial = !log)
    history[:tol] = tol

    # This doesn't yet make sense: the number of iters is smaller.
    log && reserve!(history, :resnorm, max_mv_products)

    # Actually perform CG
    iterable = bicgstabl_iterator!(x, A, b, l; Pl = Pl, tol = tol, max_mv_products = max_mv_products, kwargs...)

    if log
        history.mvps = iterable.mv_products
    end

    for (iteration, item) = enumerate(iterable)
        if log
            nextiter!(history)
            history.mvps = iterable.mv_products
            push!(history, :resnorm, iterable.residual)
        end
        verbose && @printf("%3d\t%1.2e\n", iteration, iterable.residual)
    end

    verbose && println()
    log && setconv(history, converged(iterable))
    log && shrink!(history)

    log ? (iterable.x, history) : iterable.x
end


text:
struct Foo

end

function bar(self::Foo)::String
    return "a"
end

function main()
    f = Foo()
    b = bar(f)
    println(join([b], " "))
end

main()


text:
using OrdinaryDiffEq, DiffEqProblemLibrary, Test, DiffEqBase, DiffEqCallbacks

using DiffEqProblemLibrary.ODEProblemLibrary: importodeproblems; importodeproblems()
import DiffEqProblemLibrary.ODEProblemLibrary: prob_ode_linear

prob = prob_ode_linear

time_choice(integrator) = rand()+integrator.t
affect!(integrator) =  integrator.u*=2
cb = IterativeCallback(time_choice,affect!)

sol = solve(prob,Tsit5(),callback=cb)


text:
using ExactConvolution
using Documenter

DocMeta.setdocmeta!(ExactConvolution, :DocTestSetup, :(using ExactConvolution); recursive=true)

makedocs(;
    modules=[ExactConvolution],
    authors="Ido Kessler",
    repo="https://github.com/kessido/ExactConvolution.jl/blob/{commit}{path}#{line}",
    sitename="ExactConvolution.jl",
    format=Documenter.HTML(;
        prettyurls=get(ENV, "CI", "false") == "true",
        canonical="https://kessido.github.io/ExactConvolution.jl",
        assets=String[],
    ),
    pages=[
        "Home" => "index.md",
    ],
)

deploydocs(;
    repo="github.com/kessido/ExactConvolution.jl",
    devbranch="main",
)


text:
_apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
    LAPACK.ormqr!('L','T',A,τ,B)
_apply_qr!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
    LAPACK.ormqr!('L','C',A,τ,B)
_apply_qr!(_, _, _, A::AbstractMatrix, τ::AbstractVector, B::AbstractVecOrMat) = lmul!(MatrixFactorizations.QRPackedQ(A, τ)', B)
apply_qr!(A, τ, B) = _apply_qr!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)

qlf!(A,τ) = _qlf!(MemoryLayout(A),MemoryLayout(τ),A,τ)
_qlf!(::AbstractColumnMajor,::AbstractStridedLayout,A::AbstractMatrix{T},τ::AbstractVector{T}) where T<:BlasFloat =
    LAPACK.geqlf!(A,τ)
_apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasReal =
    LAPACK.ormql!('L','T',A,τ,B)
_apply_ql!(::AbstractColumnMajor, ::AbstractStridedLayout, ::AbstractStridedLayout, A::AbstractMatrix{T}, τ::AbstractVector{T}, B::AbstractVecOrMat{T}) where T<:BlasComplex =
    LAPACK.ormql!('L','C',A,τ,B)
apply_ql!(A, τ, B) = _apply_ql!(MemoryLayout(A), MemoryLayout(τ), MemoryLayout(B), A, τ, B)

function _blockbanded_qr!(A::AbstractMatrix, τ::AbstractVector)
    M,N = Block.(blocksize(A))
    (M < N ? axes(A,1) : axes(A,2)) == axes(τ,1) || throw(DimensionMismatch(""))
    _blockbanded_qr!(A, τ, min(N,M))
    QR(A,τ.blocks)
end

function _blockbanded_qr!(A::AbstractMatrix, τ::AbstractVector, NCOLS::Block{1})
    l,u = blockbandwidths(A)
    M,N = Block.(blocksize(A))
    for K = Block(1):NCOLS
        KR = K:min(K+l,M)
        V = view(A,KR,K)
        t = view(τ,K)
        MatrixFactorizations.qrfactUnblocked!(V,t)
        for J = K+1:min(K+u,N)
            apply_qr!(V, t, view(A,KR,J))
        end
    end
    A,τ
end

function qr!(A::BlockBandedMatrix{T}) where T
    M,N = blocksize(A)
    ax1 = M < N ? axes(A,1) : axes(A,2)
    _blockbanded_qr!(A, PseudoBlockVector(zeros(T,length(ax1)), (ax1,)))
end

function ql!(A::BlockBandedMatrix{T}) where T
    l,u = blockbandwidths(A)
    M,N = blocksize(A)

    ax2 = if M < N
        throw(ArgumentError("Wide block-QL not implented"))
    else
        axes(A,2)
    end
    τ = PseudoBlockVector{T}(undef, (ax2,))

    for K = N:-1:max(N - M + 1,1)
        μ = M+K-N
        KR = Block.(max(K-u,1):μ)
        V = view(A,KR,Block(K))
        t = view(τ,Block(K-N+min(M,N)))
        qlf!(V,t)
        for J = K-1:-1:max(K-l,1)
            apply_ql!(V, t, view(A,KR,Block(J)))
        end
    end
    QL(A,τ.blocks)
end

_qr(::AbstractBlockBandedLayout, ::Tuple{Integer,Integer}, A) = qr!(BlockBandedMatrix(A, (blockbandwidth(A,1), blockbandwidth(A,1)+blockbandwidth(A,2))))
_qr(lay::AbstractBlockBandedLayout, ax::Tuple{AbstractUnitRange{Int},AbstractUnitRange{Int}}, A) = _qr(lay, map(length, ax), A)
_ql(::AbstractBlockBandedLayout, _, A) = ql!(BlockBandedMatrix(A, (blockbandwidth(A,1)+blockbandwidth(A,2),blockbandwidth(A,2))))
_factorize(::AbstractBlockBandedLayout, _, A) = qr(A)

function materialize!(Mul::MatLmulVec{<:AdjQRPackedQLayout{<:AbstractBlockBandedLayout}})
    adjQ,Bin = Mul.A,Mul.B
    Q = parent(adjQ)
    A = Q.factors
    l,u = blockbandwidths(A)
    N,M = blocksize(A)
    # impose block structure
    ax1,ax2 = axes(A)
    τ  = PseudoBlockArray(Q.τ, (length(ax1) ≤ length(ax2) ? ax1 : ax2,))
    B = PseudoBlockArray(Bin, (ax1,))
    for K = 1:min(N,M)
        KR = Block.(K:min(K+l,N))
        V = view(A,KR,Block(K))
        t = view(τ,Block(K))
        apply_qr!(V, t, view(B,KR))
    end
    Bin
end
function materialize!(Mul::MatLmulVec{<:AdjQLPackedQLayout{<:AbstractBlockBandedLayout}})
    Q = parent(Mul.A)
    A = Q.factors
    l,u = blockbandwidths(A)
    N,M = blocksize(A)
    # impose block structure
    ax1,ax2 = axes(A)
    τ  = PseudoBlockArray(Q.τ, (length(ax1) ≤ length(ax2) ? ax1 : ax2,))
    B = PseudoBlockArray(Mul.B, (ax1,))
    for K = N:-1:1
        KR = Block.(max(1,K-u):K)
        V = view(A,KR,Block(K))
        t = view(τ,Block(K))
        apply_ql!(V, t, view(B,KR))
    end
    Mul.B
end

function materialize!(Mul::MatLmulMat{<:AdjQRPackedQLayout{<:AbstractBlockBandedLayout}})
    adjQ,Bin = Mul.A,Mul.B
    Q = parent(adjQ)
    A = Q.factors
    l,u = blockbandwidths(A)
    N,M = blocksize(A)
    # impose block structure
    ax1,ax2 = axes(A)
    τ  = PseudoBlockArray(Q.τ, (length(ax1) ≤ length(ax2) ? ax1 : ax2,))
    B = PseudoBlockArray(Bin, (ax1,axes(Bin,2)))
    for K = 1:min(N,M), J = 1:blocksize(Bin,2)
        KR = Block.(K:min(K+l,N))
        V = view(A,KR,Block(K))
        t = view(τ,Block(K))
        apply_qr!(V, t, view(B,KR,Block(J)))
    end
    Bin
end

# avoid LinearALgebra Strided obsession

for Typ in (:StridedVector, :StridedMatrix, :AbstractVector, :AbstractMatrix, :LayoutMatrix, :LayoutVector)
    @eval function ldiv!(A::QR{<:Any,<:BlockSkylineMatrix}, B::$Typ)
        lmul!(adjoint(A.Q), B)
        M,N = blocksize(A.factors)
        if M == N
            materialize!(Ldiv(UpperTriangular(A.factors), view(B,1:size(A.factors,1),:)))
        else
            MN = min(M,N)
            V = view(A.factors,Block.(1:MN), Block.(1:MN))
            materialize!(Ldiv(UpperTriangular(V), view(B,1:size(V,1),:)))
        end
        B
    end
end

function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractVector)
    lmul!(adjoint(A.Q), B)
    materialize!(Ldiv(LowerTriangular(A.factors), B))
end

function ldiv!(A::QL{<:Any,<:BlockSkylineMatrix}, B::AbstractMatrix)
    lmul!(adjoint(A.Q), B)
    materialize!(Ldiv(LowerTriangular(A.factors), B))
end


text:
export LinearFilter
export fcoef, length, offset, range

import Base: length,range,isapprox

#+LinearFilter L:LinearFilter
#
# Abstract type defining a linear filter. A linear filter is defined by its coefficients and by its domain
#
abstract type LinearFilter{T<:Number} end

#+LinearFilter
# Returns filter coefficients
fcoef(c::LinearFilter) = c._fcoef

#+LinearFilter
# Returns filter length
length(c::LinearFilter)::Int = length(fcoef(c))
#+LinearFilter
# Returns filter offset
#
# *Caveat:* the first position is *0* (and not *1*)
#
# See: [[range_filter][]]
offset(c::LinearFilter)::Int = c._offset

#+LinearFilter,Convolution,Internal
#
# Computes [[range_filter][]] using primitive types.
# This allows reuse by =directConv!= for instance.
#
# *Caveat:* do not overload Base.range !!!
filter_range(size::Int,offset::Int)::UnitRange = UnitRange(-offset,size-offset-1)

#+LinearFilter L:range_filter
# Returns filter range $\Omega$
#
# Filter support is defined by
# $$
# \Omega_\alpha = [ -\text{offset}(\alpha) , \text{size}(\alpha) -\text{offset}(\alpha)  - 1 ]
# $$
range(c::LinearFilter)::UnitRange = filter_range(length(c),offset(c))

#+LinearFilter,Internal
# For convenience only, used in utests
function isapprox(f::LinearFilter{T},v::AbstractArray{T,1}) where {T<:Number}
    return isapprox(fcoef(f),v)
end

#+LinearFilter,Internal
#
# Default linear filter
#
struct LinearFilter_Default{T<:Number,N} <: LinearFilter{T}
    _fcoef::SVector{N,T}
    _offset::Int
end

#+LinearFilter,Internal
# Creates a linear filter from a coefficient vector and its associated offset
#
# *Example:*
#!linear_filter=LinearFilter(rand(3),5)
#!offset(linear_filter)
#!range(linear_filter)
#
function LinearFilter_Default(c::AbstractArray{T,1},offset::Int)  where {T<:Number}
    v=SVector{length(c),T}(c)
    return LinearFilter_Default{T,length(c)}(v,offset)
end

#+LinearFilter,Internal
#
# Default *centered* linear filter
#
# Array length has to be odd, 2n+1. Filter offset is n by construction.
#
struct LinearFilter_DefaultCentered{T<:Number,N} <: LinearFilter{T}
    _fcoef::SVector{N,T}
end

#+LinearFilter,Internal
function LinearFilter_DefaultCentered(c::AbstractArray{T,1}) where {T<:Number}
    N = length(c)
    @assert isodd(length(c)) "Centered filters must have an odd number of coefficients, $N is even"
    return LinearFilter_DefaultCentered{T,N}(SVector{N,T}(c))
end

#+LinearFilter,Internal
offset(f::LinearFilter_DefaultCentered{T,N}) where {T<:Number,N} = (N-1)>>1

# Once that we have defined LinearFilter_Default and
# LinearFilter_DefaultCentered we can unify construction. Switch to
# the right type is decided according to arguments

#+LinearFilter
#
# Creates a linear filter from its coefficients and an offset
#
# The *offset* is the position of the filter coefficient to be aligned with zero, see [[range_filter][]].
#
# *Example:*
#!f=LinearFilter([0:5;],4);
#!hcat([range(f);],fcoef(f))
#
function LinearFilter(c::AbstractArray{T,1},offset::Int)::LinearFilter  where {T}
    return LinearFilter_Default(c,offset)
end

#+LinearFilter
#
# Creates a centered linear filter, it must have an odd number of
# coefficients, $2n+1$ and is centered by construction (offset=n)
#
# *Example:*
#!f=LinearFilter([0:4;]);
#!hcat([range(f);],fcoef(f))
#
function LinearFilter(c::AbstractArray{T,1})::LinearFilter  where {T}
    return LinearFilter_DefaultCentered(c)
end


text:
export getData

function getData(m,pFor::EikonalInvParam,doClear::Bool=false)

    # extract pointers
    Mesh   	= pFor.Mesh
    Q     	= pFor.Sources
    P     	= pFor.Receivers
	n_nodes = pFor.Mesh.n+1;
    nrec  	= size(P,2)
    nsrc  	= size(Q,2)

    # allocate space for data and fields
    D  = zeros(nrec,nsrc)
    pEik = Array{EikonalParam}(nsrc);

	if pFor.useFilesForFields
		tfilename = getFieldsFileName();
		tfile = matopen(tfilename, "w");
	end

	pMem = getEikonalTempMemory(n_nodes);
	ntup = tuple(n_nodes...);
	T = zeros(Float64,ntup) # V is a temporary array of Float64.
    T1_temp = zeros(Float32,ntup);
	ordering_temp = zeros(FactoredEikonalFastMarching.EikIdxType,prod(n_nodes));
	OP_temp 	  = zeros(Int8,ntup)
	for k=1:nsrc
		src_k_loc = Q.rowval[k];
		if Mesh.dim==2
			src = zeros(Int64,2);
			cs2loc(src,src_k_loc,n_nodes);
		else
			src = zeros(Int64,3);
			cs2loc3D(src,src_k_loc,n_nodes);
		end
		m = reshape(m,ntup);
		pEik[k] = getEikonalParam(Mesh,m,src,pFor.HO);
		pEik[k].T1 = T; # Here, we set T (of Float64) for the calculation, but reuse the memory.
		pEik[k].ordering = ordering_temp;
		pEik[k].OP		 = OP_temp;
		solveFastMarchingUpwindGrad(pEik[k],pMem);
		T1_temp[:] = T;
		if Mesh.dim==2
			selfMultiplyWithAnalyticSolution2D(n_nodes,Mesh.h,src,T);
		else
			selfMultiplyWithAnalyticSolution3D(n_nodes,Mesh.h,src,T)
		end
		D[:,k] = P'*T[:];
		if doClear
			FactoredEikonalFastMarching.clear!(pEik[k]);
		elseif pFor.useFilesForFields
			write(tfile,string("T1_",k),T1_temp);
			write(tfile,string("ordering_",k),pEik[k].ordering);
			write(tfile,string("OP_",k),pEik[k].OP);
			FactoredEikonalFastMarching.clear!(pEik[k]);
		else
			# Only here we actually save the data.
			pEik[k].T1 = copy(T1_temp);
			pEik[k].ordering = copy(ordering_temp);
			pEik[k].OP = copy(OP_temp);
		end
	end
	if pFor.useFilesForFields
		close(tfile);
		gc()
	end
	pFor.eikonalParams = pEik;
    return D,pFor
end


text:
using Test
using LoopVectorization

const var"@_avx" = LoopVectorization.var"@_turbo"

using LinearAlgebra
function clenshaw(x, coeff)
  len_c = length(coeff)
  tmp = zero(x)
  ret = zero(x)
  for i = len_c:-1:2
    ret = muladd(x, 2tmp, coeff[i] - ret)
    ret, tmp = tmp, ret
  end
  ret = muladd(x, tmp, coeff[1] - ret)
  return ret
end

"""
Causes `check_args` to fail.
"""
struct FallbackArrayWrapper{T,N} <: AbstractArray{T,N}
  data::Array{T,N}
end
Base.size(A::FallbackArrayWrapper) = size(A.data)
Base.@propagate_inbounds Base.getindex(
  A::FallbackArrayWrapper,
  i::Vararg{Int,N},
) where {N} = getindex(A.data, i...)
Base.@propagate_inbounds Base.setindex!(
  A::FallbackArrayWrapper,
  v,
  i::Vararg{Int,N},
) where {N} = setindex!(A.data, v, i...)
Base.IndexStyle(::Type{<:FallbackArrayWrapper}) = IndexLinear()

const RUN_SLOW_TESTS =
  LoopVectorization.register_count() ≤ 16 ||
  !parse(Bool, get(ENV, "GITHUB_ACTIONS", "false"))