text:
using Base.Test
using CUDArt
using Knet

# import Knet: forw, back, ninputs, param, similar!, gpu, initforw, initback, push, pop, get1
# include("../src/net.jl")

include("isapprox.jl")

info("TEST 5")

isdefined(:MNIST) || include("mnist.jl")
net = [Drop(0.5), Conv(20,5), Bias(), Relu(), Pool(2),
       Drop(0.5), Conv(50,5), Bias(), Relu(), Pool(2),
       Drop(0.5), Mmul(500), Bias(), Relu(),
       Drop(0.5), Mmul(10), Bias(), XentLoss()]

nb = 128
x = KUdense(gpucopy(reshape(MNIST.xtrn[:,1:nb],28,28,1,nb)))
forw(net, copy(x))   # initializes the weights
rnn = Net(gpucopy(net)...)

setseed(42)
# @date y1 = forw(net, copy(x))
ybuf1 = Array(Any, length(net))
y = copy(x)
for n=1:length(net)
    y = forw(net[n], y)
    ybuf1[n] = cpucopy(y.arr)
end
y1 = y

setseed(42)
ybuf2 = Array(Any, length(net))
y = copy(x)
a = ()
r = rnn
inputs = Any[y]
trn = true
seq = false
using Knet: push, pop, dbg, forw, back, get1, initbatch

    initbatch(r, inputs...; trn=trn, seq=seq, a...)
    for i = 1:ninputs(r)
        n = i+nops(r)                           # input[i] goes into out[i+nops(r)]
        eltype(inputs[i]) == eltype(r.out0[n]) || error("Element type mismatch $i $n")
        trn && r.push[n] && push(r,n)         # t:140 save old input if necessary
        r.out[n] = copy!(r.out0[n], inputs[i]) 	# ; dbg(r,:out,n) # t:98 inputs can be any type of array, this will copy it to gpu or wherever
    end
    for n = 1:nops(r)
        trn && r.push[n] && push(r,n)         # t:327
        r.out[n] = forw(r.op[n], r.out[r.inputs[n]]...; y=r.out0[n], a...)     # ;dbg(r,:out,n) # t:2300
        ybuf2[n] = cpucopy(r.out[n].arr)
    end

y2 = r.out[nops(r)]

@show map(isequal, ybuf1, ybuf2)
@test @show to_host(y1.arr)==to_host(y2.arr)

y = KUdense(gpucopy(MNIST.ytrn[:,1:nb]))

# @date back(net, copy(y))
dy = copy(y)
dybuf1 = Array(Any, length(net))
for n=length(net):-1:1
    dy = back(net[n], dy)
    # @show (n, to_host(dy)[1])
    dybuf1[n] = cpucopy(dy.arr)
end
# @show length(dybuf1)

# @date back(rnn, copy(y))
r = rnn
dy = copy(y)
dx = nothing
seq = false
dybuf2 = Array(Any, nops(r))

    n = nops(r)
    if dy == nothing
        r.multi[n] || (r.dif[n] = nothing)
    elseif eltype(dy) != eltype(r.dif0[n])
        error("Element type mismatch $n")
    elseif r.multi[n]
        copy!(r.dif1[n], dy)
        r.dif[n] = axpy!(1,r.dif1[n],r.dif0[n])
    else
        r.dif[n] = copy!(r.dif0[n], dy)
    end										; dbg(r,:dif,n)
    for n = nops(r):-1:1
        if r.dif[n] == nothing
            for i in r.inputs[n]
                r.multi[i] || (r.dif[i] = nothing)
            end
        else
            dxn = Any[]
            for i in r.inputs[n]
                push!(dxn, r.multi[i] ? r.dif1[i] : r.dif0[i])
            end
            back(r.op[n], r.dif[n]; incr=seq, x=get1(r.out[r.inputs[n]]), y=r.out[n], dx=get1(dxn), a...) # t:2164

### back gives eq for layers 18..8, approxeq for layers 7..1: TODO: investigate why
        dybuf2[n] = cpucopy(get1(dxn).arr)
        dytest1 = isapprox(dybuf1[n], dybuf2[n])
        @test dytest1
        dytest2 = isequal(dybuf1[n], dybuf2[n])
        p1 = map(p->p.diff, params(rnn.op[n]))
        p2 = map(p->p.diff, params(net[n]))
        dwtest1 = reduce(&, map(isapprox, p1, p2))
        @test dwtest1
        dwtest2 = reduce(&, map(isequal, p1, p2))
        println((n, typeof(rnn.op[n]), dytest1, dytest2, dwtest1, dwtest2))
###

            for i in r.inputs[n]
                r.multi[i] && axpy!(1, r.dif1[i], r.dif0[i])            ; r.multi[i]&&dbg(r,:dif1,i)
                r.dif[i] = r.dif0[i]                                    ; dbg(r,:dif,i)
            end
            r.multi[n] && fill!(r.dif[n],0)                           ; r.multi[n]&&dbg(r,:dif,n) # t:157
        end
        r.push[n] && pop(r,n)                                    ; r.push[n]&&dbg(r,:out,n)
    end
    for i = ninputs(r):-1:1
        n = i+nops(r)
        r.push[n] && pop(r,n)                                    ; r.push[n] && dbg(r,:out,n)
        dx == nothing || copy!(dx[i], r.dif[n])
    end

# # initback(r)
# n = nops(r)

# r.dif[n] = copy!(r.dif0[n], dy)             # println("back:dy=$((idx1(getdy(r,nops(r))),))")
# for n = nops(r):-1:1
#     if r.dif[n] != nothing                  # 'nothing' represents 0 loss gradient
#         dx = map(r.inputs[n]) do i; r.inc[i]!=nothing ? r.inc[i] : (r.dif[i]=r.dif0[i]); end
#         dy = back(r.op[n], r.dif[n]; incr=true, x=get1(r.out[r.inputs[n]]), y=r.out[n], dx=get1(dx), a...)

# ### back gives eq for layers 18..8, approxeq for layers 7..1
#         dybuf2[n] = cpucopy(dy.arr)
#         dytest1 = isapprox(dybuf1[n], dybuf2[n])
#         @test dytest1
#         dytest2 = (dybuf1[n] == dybuf2[n])
#         p1 = map(p->p.diff, params(rnn.op[n]))
#         p2 = map(p->p.diff, params(net[n]))
#         @test dwtest1 = all(map(isapprox, p1, p2))
#         dwtest2 = all(map(isequal, p1, p2))
#         println((n, typeof(rnn.op[n]), dytest1, dytest2, dwtest1, dwtest2))
# ###

#         for i in r.inputs[n]; r.inc[i]!=nothing && axpy!(1, r.inc[i], r.dif[i]); end
#         r.inc[n]!=nothing && fill!(r.dif[n],0)
#     end                                     # println("op[$n]:$((typeof(r.op[n]),:x,map(idx1,getx(r,n))...,:y,idx1(gety(r,n)),:dy,idx1(getdy(r,n)),:dx,map(idx1,getdxbuf(r,n))...)) st=$(map(idx1,r.stack[1:r.sp]))")
#     r.push[n] && pop(r,n)                   # println("pop[$n]:y=$((idx1(gety(r,n)),)) st=$(map(idx1,r.stack[1:r.sp]))")
# end
# for i = ninputs(r):-1:1                     # println("in[$i]=$(map(idx1,(getinput(r,i),))) st=$(map(idx1,r.stack[1:r.sp]))")
#     n = i+nops(r)                           # println("in[$i]=$(map(idx1,(getinput(r,i),))) st=$(map(idx1,r.stack[1:r.sp]))")
#     r.push[n] && pop(r,n)
# end

if false

info("TEST 4")

isdefined(:MNIST) || include("mnist.jl")
setseed(42)
net = [Conv(20,5), Bias(), Relu(), Pool(2),
       Conv(50,5), Bias(), Relu(), Pool(2),
       Mmul(500), Bias(), Relu(),
       Mmul(10), Bias(), XentLoss()]

nb = 128
x = KUdense(gpucopy(reshape(MNIST.xtrn[:,1:nb],28,28,1,nb)))
y = KUdense(gpucopy(MNIST.ytrn[:,1:nb]))
@date y1 = forw(net, copy(x))   # also initializes the weights
rnn = Net(deepcopy(net)...)
init(rnn, copy(x); train=true)
@date y2 = forw(rnn, copy(x); y=similar(y1))
@test @show to_host(y1.arr)==to_host(y2.arr)

@date back(net, copy(y))
@date back(rnn, copy(y))
for i=nops(rnn):-1:1
    isempty(params(net[i])) && continue
    p1 = map(p->p.diff, params(rnn.op[i]))
    p2 = map(p->p.diff, params(net[i]))
    print("$i "); @test @show all(map(isapprox, p1, p2))
    print("$i "); @show all(map(isequal, p1, p2)) # 1,2,5 only approx
end

info("TEST 3")

x = KUdense(gpucopy(rand(784,10)))
net = Op[Mmul(10),QuadLoss()]
@date y1 = forw(net, copy(x))
rnn = Net(gpucopy(net)...)
init(rnn, copy(x))
@date y2 = forw(rnn, copy(x); y=similar(y1))
@test @show to_host(y1.arr)==to_host(y2.arr)
dy = rand!(copy(y1))
back(net, copy(dy))
back(rnn, copy(dy))
for i=1:nops(rnn)
    isempty(params(net[i])) && continue
    p1 = map(p->p.diff, params(rnn.op[i]))
    p2 = map(p->p.diff, params(net[i]))
    print("$i "); @test @show all(map(isequal, p1, p2)) # 1,2,5 only approx
end

info("TEST 2")

isdefined(:MNIST) || include("mnist.jl")
setseed(42)
nb = 100
x = KUdense(gpucopy(MNIST.xtrn[:,1:nb]))
net = [Mmul(64), Bias(), Relu(),
       Mmul(10), Bias(), XentLoss()]
@date y1 = forw(net, copy(x))

rnn = Net(gpucopy(net)...)
init(rnn, copy(x))
@date y2 = forw(rnn, copy(x); y=similar(y1))

# @show isapprox(y1,y2)
@test @show to_host(y1.arr)==to_host(y2.arr)

y = KUdense(gpucopy(MNIST.ytrn[:,1:nb]))
back(net, copy(y))
back(rnn, copy(y))
for i=1:nops(rnn)
    isempty(params(net[i])) && continue
    p1 = map(p->p.diff, params(rnn.op[i]))
    p2 = map(p->p.diff, params(net[i]))
    print("$i "); @test @show all(map(isequal, p1, p2)) # 1,2,5 only approx
end

info("TEST 1")
# irnn(h)=Net(Mmul(h), (Mmul(h),5), Add2(), Bias(), Relu())
# add1(h)=Net(Mmul(h), (Mmul(h),-1), Add2(), Bias(), Sigm())
# add2(h)=Net(Mmul(h), (Mmul(h),-1), Add2(), Bias(), Tanh())
# lstm(h)=Net((add1(h),0,9),      # 1. i
#             (add1(h),0,9),      # 2. f
#             (add1(h),0,9),      # 3. o
#             (add2(h),0,9),      # 4. cc
#             (Mul2(),1,4),       # 5. i*cc
#             (Mul2(),2,7),       # 6. f*c[t-1]
#             Add2(),             # 7. c
#             Tanh(),             # 8. tanh(c)
#             (Mul2(),3,8))       # 9. o * tanh(c)

function testops(a,ops)
    nops(a) == length(ops) || return false
    for i=1:nops(a)
        isa(a.op[i], ops[i]) || return false
    end
    return true
end

a = irnn(10)
aops = [Mmul, Mmul, Add2, Bias, Relu]
@test @show testops(a, aops)
@test @show a.inputs == Any[[6],[5],[1,2],[3],[4]]
@test @show a.ninputs == 1
@test @show find(a.push) == [5,6]
#@test @show a.y == [1,2,1,1,1,3]
#@test @show a.dy == [5,4,5,5,6,7]
#@test @show all(a.out .== nothing)
#@test @show all(a.out0 .== nothing)
@test @show a.stack == Any[]
@test @show a.sp == 0

b = Net(irnn(10),irnn(10))
bops = vcat(aops,aops)
@test @show testops(b, bops)
@test @show b.inputs == Any[[11],[5],[1,2],[3],[4],[5],[10],[6,7],[8],[9]]
@test @show b.ninputs == 1
@test @show find(b.push) == [5,10,11]
#@test @show b.y == [1,2,1,1,1,3,4,3,3,3,5]
#@test @show b.dy == [7,6,7,7,8,10,9,10,10,11,12]
#@test @show all(b.out .== nothing)
#@test @show b.out0 == b.out
@test @show b.stack == Any[]
@test @show b.sp == 0

c = lstm(10)
cops = [Mmul,Mmul,Add2,Bias,Sigm,Mmul,Mmul,Add2,Bias,Sigm,Mmul,Mmul,Add2,Bias,Sigm,Mmul,Mmul,Add2,Bias,Tanh,Mul2,Mul2,Add2,Tanh,Mul2]
@test @show testops(c, cops)
@test @show c.inputs == Any[[26],[25],[1,2],[3],[4],[26],[25],[6,7],[8],[9],[26],[25],[11,12],[13],[14],[26],[25],[16,17],[18],[19],[5,20],[10,23],[21,22],[23],[15,24]]
@test @show c.ninputs == 1
@test @show find(c.push) == [5,10,15,20,23,24,25,26]
#@test @show c.y == [1,2,1,1,1,3,4,3,3,3,5,6,5,5,5,7,8,7,7,7,9,10,9,11,12,13]
#@test @show c.dy == [15,14,15,15,15,17,16,17,17,17,19,18,19,19,19,21,20,21,21,21,22,23,24,25,26,27]
#@test @show all(c.out .== nothing)
#@test @show c.out0 == c.out
@test @show c.stack == Any[]
@test @show c.sp == 0

end # if false

:ok


text:
function fft(a)
    y1 = Any[]; y2 = Any[]
    n = length(a)
    if n ==1 return a end
    wn(n) = exp(-2*π*im/n)
    y_even = fft(a[1:2:end])
    y_odd = fft(a[2:2:end])
    w = 1
    for k in 1:Int(n/2)
        push!(y1, y_even[k] + w*y_odd[k])
        push!(y2, y_even[k] - w*y_odd[k])
        w = w*wn(n)
    end
    return vcat(y1,y2)
end

ext: jl

max_stars_repo_name: ethansaxenian/RosettaDecode

max_issues_repo_name: ethansaxenian/RosettaDecode

max_forks_repo_name: ethansaxenian/RosettaDecode


text:
################################################################################
# HighwayRoadwayOptions
################################################################################

@with_kw mutable struct HighwayRoadwayOptions{T} <: RoadwayOptions
    # Options
	"Length of the lanes."
	lane_length::T=4.0

	"Width of the lanes."
	lane_width::T=0.20
end

function build_roadway(opts::HighwayRoadwayOptions{T}) where {T}
	ll = opts.lane_length
	lw = opts.lane_width
	#     y
	#     ^          x3---------------------------------------------x2
	# 	  |							      w1
	# 	  |							      w2
	# (0,0) -->x     x4---------------------------------------------x1
	# 								      w3
	# 								      w4
	# 			     x5---------------------------------------------x6
	x1 = [ ll,  0.]
	x2 = [ ll,  lw]
	x3 = [ 0.,  lw]
	x4 = [ 0.,  0.]
	x5 = [ 0., -lw]
	x6 = [ ll, -lw]
	v = [0., 1.]
	w1 = Wall(x3, x2,  v)
	w2 = Wall(x4, x1, -v)
	w3 = Wall(x4, x1,  v)
	w4 = Wall(x5, x6, -v)

	# Centerlines
	f_left(s,v=0.)  = VehicleState(s,  lw/2, 0., v)
	f_right(s,v=0.) = VehicleState(s, -lw/2, 0., v)
	f_both(s,v=0.)  = VehicleState(s, -lw/2, 0., v)

	# Left lane
	left_lane = Lane(1,
					 :left_lane,
					 [w1,w2],
					 Vector{CircularWall}(),
					 StartingArea(
						 VehicleState(0.,  lw/2, 0.0,  0.00),
						 VehicleState(0.,  lw/2, 0.0, -0.05),
						 VehicleState(ll,  lw/2, 0.4,  0.05),
						 ),
					 f_left,
					 )

	# Right lane
	right_lane = Lane(2,
					  :right_lane,
					  [w3,w4],
					  Vector{CircularWall}(),
					  StartingArea(
						  VehicleState(0., -lw/2, 0.0,  0.00),
						  VehicleState(0., -lw/2, 0.0, -0.05),
					      VehicleState(ll, -lw/2, 0.4,  0.05),
					      ),
					  f_right,
					  )

	# Both lane
	both_lane = Lane(3,
					 :both_lane,
					 [w1,w4],
					 Vector{CircularWall}(),
					 StartingArea(
						 VehicleState(0., -lw/2, 0.0,  0.00),
						 VehicleState(0., -lw/2, 0.0, -0.05),
						 VehicleState(ll, -lw/2, 0.4,  0.05),
						 ),
					 f_both,
					 )

	lanes = [left_lane, right_lane, both_lane]
	roadway = Roadway(lanes, opts)
    return roadway
end


text:
function remove_interpretted_mode_overdub(primal_meth::Method)
    primal_sigp = parameter_typeof(primal_meth.sig)

    overdub_args_sigt = Tuple{MixModeDebugCtx, primal_sigp...}
    overdub_meth = only(methods(Cassette.overdub, overdub_args_sigt))
    return Base.delete_method(overdub_meth)
end

function remove_interpretted_mode_overdub(bp::BreakpointSignature)
    return make_interpretted_mode_overdub(bp.f)
end


text:
using ITensors,
      Test

@testset "MPO Basics" begin

  N = 10
  sites = SiteSet(N,2)
  O = MPO(sites)
  @test length(O) == N

  O[1] = ITensor(sites[1], prime(sites[1]))
  @test hasindex(O[1],sites[1])
  @test hasindex(O[1],prime(sites[1]))

  P = copy(O)
  @test hasindex(P[1],sites[1])
  @test hasindex(P[1],prime(sites[1]))

  @testset "inner" begin
    phi = randomMPS(sites)
    K = randomMPO(sites)
    psi = randomMPS(sites)
    phidag = dag(phi)
    prime!(phidag)
    phiKpsi = phidag[1]*K[1]*psi[1]
    for j = 2:N
      phiKpsi *= phidag[j]*K[j]*psi[j]
    end
    @test phiKpsi[] ≈ inner(phi,K,psi)
  end

end


text:
include("util.jl")

import Base: *

# No overloads for + and - methods
# Julia will automatically call the according per element op for addition/subtraction if dimensions are the same
# Otherwise a call to the dotted operation (.+) is needed which calls the broadcasted/per-element op

# To fix ambiguity with *(AbstractArray, AbstractArray) method in LinearAlgebra
# Define specialized function for every AbstractArray subtype (and Tensor struct)
for S in (Tensor, AbstractArray, AbstractVector, AbstractMatrix)
@eval begin
        function *(t::Tensor, a::$S)
            Tensor(data(t) * data(a), op=typeof(*), prev=(t, a))
        end
    end
end

max_stars_repo_name: leonvol/magrad

max_issues_repo_name: leonvol/magrad

max_forks_repo_name: leonvol/magrad


text:
using RLEVectors
using DataFrames
using CSV

results_file = "/Users/phaverty/.julia/dev/RLEVectors/benchmark/rle.timings.csv"

macro timeit(ex)
# like @time, but returning the timing rather than the computed value
  return quote
    #gc_disable()
    local val = $ex # compile
    local t0 = time()
    for i in 1:1e4 val = $ex end
    local t1 = time()
    #gc_enable()
    t1-t0
  end
end

foo = IntegerRle(Vector{Int32}(collect(1:1000)), Vector{Int32}(collect(5:5:5000)))
length(foo)
timings = DataFrame()
timings[:language] = "julia"
timings[:language_version] = VERSION
timings[:date] = chomp(read(`date "+%Y-%m-%d"`, String))
timings[:indexing] = @timeit foo[100]
timings[:range_indexing] = @timeit foo[801:900]
timings[:setting] = @timeit foo[800] = 5
timings[:range_setting] = @timeit foo[801:900] = 1:100
timings[:scalar_add] = @timeit broadcast(+,foo,4)
timings[:length] = @timeit length(foo)
timings[:nrun] = @timeit nrun(foo)
timings[:mean] = @timeit mean(foo)
timings[:max] = @timeit maximum(foo)
timings[:width] = @timeit rwidth(foo)
timings[:last] = @timeit rlast(foo)
timings[:first] = @timeit rfirst(foo)
timings[:add_two_rles] = @timeit broadcast(+,foo,foo)
timings[:disjoin] = @timeit disjoin(foo.runends,foo.runends)
timings[:which] = @timeit findall(in(foo),[800,200,357])
timings[:scalar_less] = @timeit broadcast(<,foo,3)
timings[:median] = @timeit median(foo)
timings[:which_max] = @timeit findmax(foo)

bdf = CSV.read(results_file,header=true);

for n in names(timings)
  if !(n in names(bdf))
    bdf[n] = NaN
  end
end

bdf = vcat(bdf,timings)

CSV.write(results_file, bdf)

jdf = timings
rdf = bdf[ bdf[:language] .== "R",:]
rdf = rdf[end,:]
rdf = rdf[4:end]
jdf = jdf[4:end]

for n in names(rdf)
       println(n)
       println( rdf[n] / jdf[n])
end
r_over_julia = zeros(ncol(rdf))

for i in 1:length(r_over_julia)
    r_over_julia[i] = log2(rdf[i][1] / jdf[i][1])
end

### Plotting
using VegaLite
using FileIO

## Performance Relative to R
date = bdf[end,:date]
relative_perf_file = "/Users/phaverty/.julia/dev/RLEVectors/benchmark/plots/benchmark_rle_vectors.$(date).svg"

df = DataFrame(:function => names(rdf), :ratio => r_over_julia)
df[:,:function] = string.(df[:,:function])
df |> @vlplot(
  :bar,
  x=:function,
  y=:ratio,
  width=500,
  title = "Relative Performance of R and Julia Rle Vectors") |> save(relative_perf_file)

  #Guide.ylabel("Elapsed Time: log2(R/julia)"),
  #Guide.xticks(orientation=:vertical),
  #Scale.color_continuous(minvalue=-15,maxvalue=15),color=r_over_julia,
  #Geom.hline(color="black"),yintercept=[0],Guide.xlabel("")

current_relative_perf_file = "/Users/phaverty/.julia/dev/RLEVectors/benchmark/plots/benchmark_rle_vectors.svg"
cp(relative_perf_file, current_relative_perf_file, force=true)

## Performance over time
timeline_file = "/Users/phaverty/.julia/dev/RLEVectors/benchmark/plots/benchmark_rle_vectors.$(date).timeline.svg"
jdf = bdf[ bdf[:,:language] .== "julia", 3:end ]
melted_bdf = melt(jdf, :date)
melted_bdf |>
  @vlplot(
    :line,
    x="date:t",
    y=:value,
    color=:variable,
    width=500
    ) |> save(timeline_file)
#Guide.xlabel("Date"), Geom.line,
#Scale.y_log10, Guide.ylabel("log2 elapsed seconds (1e4 runs)"))
current_timeline_file = "/Users/phaverty/.julia/dev/RLEVectors/benchmark/plots/benchmark_rle_vectors.timeline.svg";
cp(timeline_file, current_timeline_file, force=true)

## Profiling
using ProfileView
foo + foo; Profile.clear(); @profile for i in 1:1e4 foo + foo end; ProfileView.view()
foo .< 3; Profile.clear(); @profile for i in 1:1e4 foo .< 3 end; ProfileView.view()
foo .+ 3; Profile.clear(); @profile for i in 1:1e4 foo .+ 3 end; ProfileView.view()
sum(foo); Profile.clear(); @profile for i in 1:1e4 sum(foo) end; ProfileView.view()
findall(in([800,300,357]), foo); Profile.clear(); @profile for i in 1:1e5 findall(in([800,300,357]), foo) end; ProfileView.view()
median(foo); Profile.clear(); @profile for i in 1:1e4 median(foo) end; ProfileView.view()
rfirst(foo); Profile.clear(); @profile for i in 1:1e6 rfirst(foo) end; ProfileView.view()
rwidth(foo); Profile.clear(); @profile for i in 1:1e6 rwidth(foo) end; ProfileView.view()
foo[800] = 5; Profile.clear(); @profile for i in 1:1e6 foo[800] = 5 end; ProfileView.view()
foo[802] = 5; Profile.clear(); @profile for i in 1:1e6 foo[802] = 5 end; ProfileView.view()
foo[801:900] = 1:100; Profile.clear(); @profile for i in 1:1e5 foo[801:900] = 1:100 end; ProfileView.view()

x = collect(2:6:5000)
y = collect(6:4:5000)

foo = disjoin(x, y)
goo = disjoin2(x, y)
goo = goo[1:1667]
foo == goo

@time for i in 1:1e6 disjoin(x, y) end
@time for i in 1:1e6 disjoin2(x, y) end
Profile.clear(); @profile for i in 1:1e6 disjoin2(x, y) end; ProfileView.view()


text:
using VideoIO
opencamera()

max_stars_repo_name: terasakisatoshi/juliaExer

max_issues_repo_name: terasakisatoshi/juliaExer

max_forks_repo_name: terasakisatoshi/juliaExer


text:
push!(LOAD_PATH,"../../src/")
import UNSflow

alphadef1 = UNSflow.ConstDef(45. *pi/180)
alphadef2 = UNSflow.ConstDef(5. *pi/180)

hdef = UNSflow.ConstDef(0.)

udef = UNSflow.ConstDef(1.)

full_kinem1 = UNSflow.KinemDef(alphadef1, hdef, udef)
full_kinem2 = UNSflow.KinemDef(alphadef2, hdef, udef)

pvt = 0.0

geometry = "FlatPlate"

lespcrit = [0.11;]

surf1 = UNSflow.TwoDSurf(geometry, pvt, full_kinem1, lespcrit)
surf2 = UNSflow.TwoDSurf(geometry, pvt, full_kinem2, lespcrit, initpos = [15.5 ;-0.5 ])

surf = [surf1;surf2]

curfield = UNSflow.TwoDFlowField()

dtstar = UNSflow.find_tstep(alphadef)

t_tot = 5.

nsteps = Int(round(t_tot/dtstar))+1

startflag = 0

writeflag = 1

writeInterval = t_tot/18.

#delvort = delSpalart(500, 12, 1e-5)
delvort = UNSflow.delNone()

mat, surf, curfield = UNSflow.ldvm(surf, curfield, nsteps, dtstar,startflag, writeflag, writeInterval, delvort)

UNSflow.makeVortPlots2D()

UNSflow.makeForcePlots2D()

UNSflow.cleanWrite()


text:
@doc raw"""
    VehicleTracking(y::Matrix{Fl}, ρ::Fl, H::Matrix{Fl}, Q::Matrix{Fl}) where Fl

The vehicle tracking example illustrates a model where there are no hyperparameters,
the user defines the parameters ``\rho``, ``H`` and ``Q`` and the model gives the predicted and filtered speed and
position. In this case, $y_t$ is a $2 \times 1$ observation vector representing the corrupted measurements
of the vehicle's position on the two-dimensional plane in instant $t$.

The position and speed in each dimension compose the state of the vehicle. Let us refer to ``x_t^{(d)}`` as
the position on the axis $d$ and to ``\dot{x}^{(d)}_t`` as the speed on the axis $d$ in instant ``t``. Additionally,
let ``\eta^{(d)}_t`` be the input drive force on the axis ``d``, which acts as state noise. For a single dimension,
we can describe the vehicle dynamics as

```math
\begin{equation}
    \begin{aligned}
        & x_{t+1}^{(d)} = x_t^{(d)} + \Big( 1 - \frac{\rho \Delta_t}{2} \Big) \Delta_t \dot{x}^{(d)}_t + \frac{\Delta^2_t}{2} \eta_t^{(d)}, \\
        & \dot{x}^{(d)}_{t+1} = (1 - \rho) \dot{x}^{(d)}_{t} + \Delta_t \eta^{(d)}_t,
    \end{aligned}\label{eq_control}
\end{equation}
```

We can cast the dynamical system as a state-space model in the following manner:
```math
\begin{align*}
    y_t &= \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{bmatrix} \alpha_{t+1} + \varepsilon_t, \\
    \alpha_{t+1} &= \begin{bmatrix} 1 & (1 - \tfrac{\rho \Delta_t}{2}) \Delta_t & 0 & 0 \\ 0 & (1 - \rho) & 0 & 0 \\ 0 & 0 & 1 & (1 - \tfrac{\rho \Delta_t}{2}) \\ 0 & 0 & 0 & (1 - \rho) \end{bmatrix} \alpha_{t} + \begin{bmatrix} \tfrac{\Delta^2_t}{2} & 0 \\ \Delta_t & 0 \\ 0 & \tfrac{\Delta^2_t}{2} \\ 0 & \Delta_t \end{bmatrix} \eta_{t},
\end{align*}
```

See more on [Vehicle tracking](@ref)
"""
mutable struct VehicleTracking <: StateSpaceModel
    system::LinearMultivariateTimeInvariant

    function VehicleTracking(y::Matrix{Fl}, ρ::Fl, H::Matrix{Fl}, Q::Matrix{Fl}) where Fl
        p = 2
        Z = kron(Matrix{Fl}(I, p, p), [1.0 0.0])
        T = kron(Matrix{Fl}(I, p, p), [1 (1 - ρ / 2); 0 (1 - ρ)])
        R = kron(Matrix{Fl}(I, p, p), [0.5; 1])
        d = zeros(Fl, p)
        c = zeros(Fl, 4)
        H = H
        Q = Q

        system = LinearMultivariateTimeInvariant{Fl}(y, Z, T, R, d, c, H, Q)

        return new(system)
    end
end

function default_filter(model::VehicleTracking)
    Fl = typeof_model_elements(model)
    steadystate_tol = Fl(1e-5)
    a1 = zeros(Fl, num_states(model))
    skip_llk_instants = length(a1)
    P1 = Fl(1e6) .* Matrix{Fl}(I, num_states(model), num_states(model))
    return MultivariateKalmanFilter(2, a1, P1, skip_llk_instants, steadystate_tol)
end


text:
import Base: size

const NDIM = 2

import CartesianGrids: cellsize, origin
import RigidBodyTools: assign_velocity!
import ImmersedLayers: normals, areas

include("navierstokes/linesource.jl")

#=
The goal of this should be to change NavierStokes into
a more general type for any PDEs, similar to the structures
used for OrdinaryDiffEq.

Rule of thumb: keep the typing generic!

Other things to keep in mind:
- For MovingPoints problems, the immersion operators will get changed. These
   need to be easily accessible from the overall system. We can endow
   the term type with a parameter for types that must be updated if the points
   update. Each such type will need to be accompanied by an operator for
   setting it up (and possibly a separate one for updating it).
- Some operators, such as Rf and Ef, may get used by more than one term.
=#

"""
$(TYPEDEF)

A system type that utilizes a grid of `NX` x `NY` dual cells and `N` Lagrange forcing
points to solve the discrete Navier-Stokes equations in vorticity form. The
parameter `static_points` specifies whether the forcing points remain static in the
grid. It should be set to `false` if a supplied motion requires that the points move.

# Constructors:

`NavierStokes(Re,Δx,xlimits,ylimits,Δt
              [,freestream = (0.0, 0.0)]
              [,pulses=nothing])` specifies the Reynolds number `Re`, the grid
              spacing `Δx`, the dimensions of the domain in the tuples `xlimits`
              and `ylimits` (excluding the ghost cells), and the time step size `Δt`.
              The other arguments are optional. The `freestream` argument can be
              passed as either a tuple (a static freestream) or a `RigidBodyMotion`
              for a time-varying freestream. The `pulses` argument can be
              used to pass in one or more spatiotemporal pulses.

`NavierStokes(Re,Δx,xlimits,ylimits,Δt,bodies::Body/BodyList
              [,flow_side=ExternalInternalFlow]
              [,ddftype=CartesianGrids.Yang3])` passes the body
              information. This constructor
              sets the motions of the body/ies to be stationary.
              The same optional arguments used for the basic constructor
              also apply for this one. In addition, the `flow_side` can be set to
              `ExternalFlow` (default), `InternalFlow`, or `ExternalInternalFlow`.
              However, it is forced to `ExternalInternalFlow` for open Bodies
              (like `Plate` type).

`NavierStokes(Re,Δx,xlimits,ylimits,Δt,bodies::Body/BodyList,
              motions::Motion/MotionList
              [,static_points=false])`
              passes the body and associated motion information.
              The list of motions must be the same length as the list of bodies.
              The same optional arguments used for the other constructors
              also apply for this one. In addition, `static_points` can
              be set to `true` if the supplied motion should not cause the
              points to move.

"""
mutable struct NavierStokes{NX, NY, N, MT<:PointMotionType, FS<:FreestreamType, SD<:FlowSide, DDF<:CartesianGrids.DDFType, FT,
                            SP, WP, VP, SPP, GVP, VBP, SBP, FST,
                            RHST, CDT, COT, DLT, VCT, PCT, SCT, SVCT, STCT}
    # Physical Parameters
    "Reynolds number"
    Re::Float64
    "Free stream velocities"
    U∞::FST
    "Bodies"
    bodies::Union{BodyList,Nothing}
    "Body motions"
    motions::Union{RigidMotionList,DirectlySpecifiedMotionList,Nothing}
    "Pulses"
    pulses::Union{Vector{ModulatedField},Nothing}

    ## ADD τ_bc HERE (PRESCRIBED LINE SOURCE)
    τ_bc::Union{Vector{PrescribedLineSource},Nothing}

    # Discretization
    "Grid metadata"
    grid::CartesianGrids.PhysicalGrid{2}
    "Time step"
    Δt::Float64

    # Operators
    "Laplacian operator"
    L::CartesianGrids.Laplacian

    # Layers
    dlf::Union{DoubleLayer,Nothing} # used for viscous surface terms
    slc::Union{SingleLayer,Nothing} # used for scalar potential field in velocity
    sln::Union{SingleLayer,Nothing} # might not be used

    # Body coordinate data, if present
    # if a static problem, these coordinates are in inertial coordinates
    # if a non-static problem, in their own coordinate systems
    points::VectorData{N,Float64}

    # Pre-stored regularization and interpolation matrices (if present)
    Rf::Union{RegularizationMatrix,Nothing} # faces (edges)
    Ef::Union{InterpolationMatrix,Nothing}
    Cf::Union{AbstractMatrix,Nothing}

    # Operators
    f :: FT

    # prototypes
    state_prototype :: SP
    vorticity_prototype :: WP
    velocity_prototype :: VP
    scapot_prototype :: SPP
    gradv_prototype :: GVP
    surfacevel_prototype :: VBP
    surfacescalar_prototype :: SBP

    # Caches
    rhs_cache :: RHST
    convderiv_cache :: CDT
    constraintop_cache :: COT
    doublelayer_cache :: DLT
    velocity_cache :: VCT
    pressure_cache :: PCT
    streamfunction_cache :: SCT
    surfacevel_cache :: SVCT
    surfacetraction_cache :: STCT

end

function NavierStokes(Re::Real, Δx::Real, xlimits::Tuple{Real,Real},ylimits::Tuple{Real,Real}, Δt::Real;
                       freestream::FST = (0.0, 0.0),
                       bodies::Union{BodyList,Nothing} = nothing,
                       motions::Union{RigidMotionList,DirectlySpecifiedMotionList,Nothing} = nothing,
                       pulses::PT = nothing,
                       τ_bc = nothing,
                       static_points = true,
                       flow_side::Type{SD} = ExternalFlow,
                       ddftype=CartesianGrids.Yang3) where {FST,PT,SD<:FlowSide}

    g = PhysicalGrid(xlimits,ylimits,Δx)
    NX, NY = size(g)

    N = numpts(bodies)

    α = Δt/(Re*Δx^2)

    vorticity_prototype = Nodes{Dual,NX,NY,Float64}()
    velocity_prototype = Edges{Primal,NX,NY,Float64}()
    scapot_prototype = Nodes{Primal,NX,NY,Float64}()
    gradv_prototype = EdgeGradient{Primal,Dual,NX,NY,Float64}()
    surfacevel_prototype = VectorData(N)
    surfacescalar_prototype = ScalarData(N)

    L = plan_laplacian(vorticity_prototype,with_inverse=true)

    pulsefields = _process_pulses(pulses,vorticity_prototype,g)

    # for now, if there are any bodies that are Open,
    # then force flow_side to ExternalInternalFlow.
    # but should be more flexible here
    flow_side_internal = _any_open_bodies(bodies) ? ExternalInternalFlow : flow_side

    points, dlf, slc, sln, Rf, Ef, Cf =
              _immersion_operators(bodies,g,flow_side_internal,ddftype,
                        velocity_prototype,scapot_prototype,surfacevel_prototype)

    viscous_L = plan_laplacian(similar(vorticity_prototype),
                               factor=1/(Re*Δx^2))

    rhs_cache = RHSCache(similar(velocity_prototype))

    convderiv_cache = ConvectiveDerivativeCache(similar(velocity_prototype),
                                                similar(gradv_prototype),
                                                similar(gradv_prototype),
                                                similar(gradv_prototype))
    constraintop_cache = ConstraintOperatorCache(similar(velocity_prototype),
                                                 similar(surfacevel_prototype))

    doublelayer_cache = DoubleLayerCache(similar(velocity_prototype),
                                         similar(surfacevel_prototype))

    velocity_cache = VelocityCache(similar(vorticity_prototype),
                                   similar(velocity_prototype),
                                   similar(scapot_prototype),
                                   VectorData(N),
                                   ScalarData(N))

    pressure_cache = PressureCache(similar(velocity_prototype),
                                   similar(velocity_prototype),
                                   similar(scapot_prototype))

    streamfunction_cache = StreamfunctionCache(similar(vorticity_prototype))

    surfacevel_cache = SurfaceVelocityCache(similar(surfacevel_prototype))

    surfacetraction_cache = SurfaceTractionCache(similar(surfacevel_prototype),
                                                similar(surfacevel_prototype),
                                                similar(surfacescalar_prototype))

    if isnothing(bodies)
      state_prototype = solvector(state=vorticity_prototype)
      f = ConstrainedODEFunction(ns_rhs!,viscous_L,_func_cache=state_prototype)
    else
      if static_points
        state_prototype = solvector(state=vorticity_prototype,constraint=surfacevel_prototype)
        f = ConstrainedODEFunction(ns_rhs!,bc_constraint_rhs!,
                                      ns_op_constraint_force!,bc_constraint_op!,
                                      viscous_L,_func_cache=state_prototype)
      else
        state_prototype = solvector(state=vorticity_prototype,constraint=surfacevel_prototype,aux_state=zero_body_state(bodies))
        rhs! = ConstrainedSystems.r1vector(state_r1 = ns_rhs!,aux_r1 = rigid_body_rhs!)
        f = ConstrainedODEFunction(rhs!,bc_constraint_rhs!,
                                      ns_op_constraint_force!,bc_constraint_op!,
                                      viscous_L,_func_cache=state_prototype,
                                      param_update_func=update_immersion_operators!)
      end

    end

    NavierStokes{NX, NY, N, _motiontype(static_points), _fstype(FST), flow_side_internal,
                  ddftype, typeof(f), typeof(state_prototype),typeof(vorticity_prototype),
                  typeof(velocity_prototype),typeof(scapot_prototype),typeof(gradv_prototype),
                  typeof(surfacevel_prototype),typeof(surfacescalar_prototype),
                  FST,typeof(rhs_cache),
                  typeof(convderiv_cache), typeof(constraintop_cache), typeof(doublelayer_cache),
                  typeof(velocity_cache), typeof(pressure_cache),
                  typeof(streamfunction_cache),typeof(surfacevel_cache),typeof(surfacetraction_cache)}(
                          Re, freestream, bodies, motions, pulsefields,
                          _line_source(τ_bc,Edges(Primal,size(g)),g),
                          g, Δt, # rk,
                          L,
                          dlf,slc,sln,
                          points, Rf, Ef, Cf,
                          f,state_prototype,
                          vorticity_prototype, velocity_prototype, scapot_prototype,
                          gradv_prototype, surfacevel_prototype, surfacescalar_prototype,
                          rhs_cache, convderiv_cache, constraintop_cache,
                          doublelayer_cache,
                          velocity_cache, pressure_cache, streamfunction_cache,
                          surfacevel_cache,surfacetraction_cache)
end

NavierStokes(Re,Δx,xlim,ylim,Δt,bodies::BodyList;
        motions=RigidMotionList(map(x -> RigidBodyMotion(0.0,0.0),bodies)),kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt;bodies=bodies,motions=motions,kwargs...)

NavierStokes(Re,Δx,xlim,ylim,Δt,body::Body;kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt,BodyList([body]);kwargs...)

NavierStokes(Re,Δx,xlim,ylim,Δt,body::Body,τ_bc::PrescribedLineSource;kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt,BodyList([body]),τ_bc=τ_bc;kwargs...)

function NavierStokes(Re,Δx,xlim,ylim,Δt,bodies::BodyList,motions::Union{RigidMotionList,DirectlySpecifiedMotionList};static_points=false,kwargs...)
    length(bodies) == length(motions) || error("Inconsistent lengths of bodies and motions lists")
    NavierStokes(Re,Δx,xlim,ylim,Δt,bodies;motions=motions,static_points=static_points,kwargs...)
end

NavierStokes(Re,Δx,xlim,ylim,Δt,body::Body,motion::RigidBodyMotion,static_points=false;kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt,BodyList([body]),RigidMotionList([motion]),static_points=static_points;kwargs...)

NavierStokes(Re,Δx,xlim,ylim,Δt,body::Body,motion::DirectlySpecifiedMotion,static_points=false;kwargs...) =
        NavierStokes(Re,Δx,xlim,ylim,Δt,BodyList([body]),DirectlySpecifiedMotionList([motion]),static_points=static_points;kwargs...)

function Base.show(io::IO, sys::NavierStokes{NX,NY,N,MT,FS,SD}) where {NX,NY,N,MT,FS,SD}
    mtype = (MT == StaticPoints) ? "static" : "moving"
    fsmsg = (FS == StaticFreestream) ? "Static freestream = $(sys.U∞)" : "Variable freestream"
    sdmsg = (N == 0) ? "Unbounded" : ((SD == ExternalFlow) ? "External flow" : ((SD == InternalFlow) ? "Internal flow" : "External/internal"))
    println(io, "$sdmsg Navier-Stokes system on a grid of size $NX x $NY and $N $mtype immersed points")
    println(io, "   $fsmsg")
    if N > 0
      bdmsg = (length(sys.bodies) == 1) ? "1 body" : "$(length(sys.bodies)) bodies"
      println(io, "   $bdmsg")
    end
end

# Routines to set up the immersion operators

function _immersion_operators(bodies::BodyList,g::PhysicalGrid,flow_side::Type{SD},ddftype,
                              gridvector_prototype::GridData{NX,NY},gridscalar_prototype::GridData{NX,NY},
                              surfvector_prototype::PointData{N}) where {NX,NY,N,SD<:ViscousFlow.FlowSide}

  points = VectorData(collect(bodies))
  numpts(bodies) == N || error("Inconsistent size of bodies")

  body_areas = areas(bodies)
  body_normals = normals(bodies)

  if !(flow_side==ExternalInternalFlow)
    dlf = DoubleLayer(bodies,g,gridvector_prototype)
    slc = SingleLayer(bodies,g,gridscalar_prototype)
    #sln = SingleLayer(bodies,g,Sn)
    sln = nothing
  else
    dlf = nothing
    slc = nothing
    sln = nothing
  end

  #regop = Regularize(points,cellsize(g);I0=CartesianGrids.origin(g),weights=body_areas.data,ddftype=ddftype)
  regop = _regularization(points,g,bodies,ddftype)

  Rf = RegularizationMatrix(regop,surfvector_prototype,gridvector_prototype) # Used by B₁ᵀ
  Ef = InterpolationMatrix(regop,gridvector_prototype,surfvector_prototype) # Used by constraint_rhs! and B₂

  regopfilt = Regularize(points,cellsize(g);I0=CartesianGrids.origin(g),filter=true,weights=cellsize(g)^2,ddftype=ddftype)

  Ẽf = InterpolationMatrix(regopfilt,gridvector_prototype,surfvector_prototype)
  Cf = sparse(Ẽf*Rf)

  return points, dlf, slc, sln, Rf, Ef, Cf
end

_immersion_operators(::Nothing,a...) =
    VectorData(0), nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing, nothing

# For updating the system with body data

function update_immersion_operators!(sys::NavierStokes{NX,NY,N,MT,FS,SD,DDF},bodies::BodyList) where {NX,NY,N,MT,FS,SD,DDF<:CartesianGrids.DDFType}
    @unpack velocity_prototype, scapot_prototype, surfacevel_prototype = sys
    sys.bodies = deepcopy(bodies)
    sys.points, sys.dlf, sys.slc, sys.sln, sys.Rf, sys.Ef, sys.Cf =
      _immersion_operators(sys.bodies,sys.grid,SD,DDF,velocity_prototype,scapot_prototype,surfacevel_prototype)
    return sys
end

update_immersion_operators!(sys::NavierStokes,body::Body) = update_immersion_operators!(sys,BodyList([body]))

function update_immersion_operators!(sys::NavierStokes,x::AbstractVector)
    tl! = RigidTransformList(x)
    tl!(sys.bodies)
    update_immersion_operators!(sys,sys.bodies)
end

# The form passed to ConstrainedODEFunction
update_immersion_operators!(sys::NavierStokes,u,sys_old::NavierStokes,t) =
    update_immersion_operators!(sys,aux_state(u))

"""
    setstepsizes(Re[,gridRe=2][,cfl=0.5][,fourier=0.5]) -> Float64, Float64

Set the grid cell spacing and time step size based on the Reynolds number `Re`,
the grid Reynolds number `gridRe`, cfl number `cfl`, and grid Fourier number `fourier`.
The last three parameters all have default values.

# Example

Here is an example of setting parameters based on Reynolds number 100 (with
  default choices for grid Reynolds number, CFL number, and Fourier number):
```jldoctest
julia> Δx, Δt = setstepsizes(100)
(0.02, 0.01)
```
"""
function setstepsizes(Re::Real; gridRe = 2.0, cfl = 0.5, fourier = 0.5)
    Δx = gridRe/Re
    Δt = min(cfl*Δx,fourier*Δx^2*Re)
    return Δx, Δt
end

# some convenience functions
"""
    size(sys::NavierStokes,d::Int) -> Int

Return the number of indices of the grid used by `sys` along dimension `d`.
"""
size(sys::NavierStokes{NX,NY},d::Int) where {NX,NY} = d == 1 ? NX : NY

"""
    size(sys::NavierStokes) -> Tuple{Int,Int}

Return a tuple of the number of indices of the grid used by `sys`
"""
size(sys::NavierStokes{NX,NY}) where {NX,NY} = (size(sys,1),size(sys,2))

"""
    cellsize(sys::NavierStokes) -> Float64

Return the grid cell size of system `sys`
"""
cellsize(sys::NavierStokes) = cellsize(sys.grid)

"""
    timestep(sys::NavierStokes) -> Float64

Return the time step size of system `sys`
"""
timestep(sys::NavierStokes) = sys.Δt

"""
    origin(sys::NavierStokes) -> Tuple{Int,Int}

Return a tuple of the indices of the primal node that corresponds to the
physical origin of the coordinate system used by `sys`. Note that these
indices need not lie inside the range of indices occupied by the grid.
For example, if the range of physical coordinates occupied by the grid
is (1.0,3.0) x (2.0,4.0), then the origin is not inside the grid.
"""
origin(sys::NavierStokes) = origin(sys.grid)

"""
    timerange(tf,sys::NavierStokes)

Create a range of times, starting at the t = Δt (the time step of `sys`),
and ending at t = `tf`.
"""
timerange(tf,sys) = timestep(sys):timestep(sys):tf

# Wrap the output of the motion evaluation in VectorData
@inline assign_velocity!(u::VectorData,a...) = (assign_velocity!(u.u,u.v,a...); u)

@inline normals(sys::NavierStokes) = (!isnothing(sys.bodies)) ? normals(sys.bodies) : nothing

"""
    freestream(t,sys::NavierStokes) -> Tuple

Return the value of the freestream in `sys` at time `t` as a tuple.
"""
freestream(t::Real,sys::NavierStokes{NX,NY,N,MT,StaticFreestream}) where {NX,NY,N,MT} = sys.U∞

function freestream(t::Real,sys::NavierStokes{NX,NY,N,MT,VariableFreestream}) where {NX,NY,N,MT}
    _,ċ,_,_,_,_ = sys.U∞(t)
    return reim(ċ)
end

"""
    newstate(sys::NavierStokes)

Return a new (zero) instance of the state vector for `sys`.
"""
newstate(sys::NavierStokes) = zero(sys.state_prototype)

"""
    newstate(s::AbstractSpatialField,sys::NavierStokes)

Return an instance of the state vector for `sys`, assigned the
data in the spatial field `s`.
"""
function newstate(s::AbstractSpatialField,sys::NavierStokes)
  u = newstate(sys)
  gf = GeneratedField(state(u),s,sys.grid)
  state(u) .= cellsize(sys)*gf()
  return u
end

"""
    flowside(sys::NavierStokes) -> FlowSide

Returns the side of the body on which the flow is to be computed.
"""
flowside(::NavierStokes{NX,NY,N,MT,FS,SD}) where {NX,NY,N,MT,FS,SD} = SD

# Other functions
_hasfilter(sys::NavierStokes) = !isnothing(sys.Cf)
_motiontype(isstatic::Bool) = isstatic ? StaticPoints : MovingPoints
  _fstype(F) = F <: Union{RigidBodyMotion,Kinematics} ? VariableFreestream : StaticFreestream

_body_closure_type(b::T) where {T<:Body{N,C}} where {N,C} = C

_any_open_bodies(nothing) = false
_any_open_bodies(bodies::BodyList) =  any(b -> _body_closure_type(b) == RigidBodyTools.OpenBody,bodies)

_regularization(sys::NavierStokes{NX, NY, N, MT, FS, SD, DDF}) where {NX,NY,N,MT,FS,SD,DDF} =
        _regularization(sys.points,sys.grid,sys.bodies,DDF)

_regularization(points,g,bodies,ddftype) = Regularize(points,cellsize(g),
                                I0=CartesianGrids.origin(g),weights=areas(bodies).data,ddftype=ddftype)

function _line_source(lineparams::Vector{<:LineSourceParams},u::VectorGridData,g::PhysicalGrid)
  τ_bc = PrescribedLineSource[]
  for p in lineparams
    push!(τ_bc,PrescribedLineSource(p,u,g))
  end
  τ_bc
end

_line_source(lineparams::LineSourceParams,u::VectorGridData,g::PhysicalGrid) =
    _line_source([lineparams],u,g)

_line_source(::Nothing,u,g) = nothing

# function set_linesource_strength!(τcoor::body,qu::Scalar{T},qv::Scalar{T}) where {T<:Real}
    ### can only set a uniform strength

#     @unpack τ_bc = sys

#     TractionStrength_coor = VectorData(collect(τcoor))
#     fill!(TractionStrength_coor.u,qu)
#     fill!(TractionStrength_coor.v,qv)
#     set_linesource_strength!(sys.τ_bc[*],TractionStrength_coor)
### would need to know the order of τ to get *

# end

function set_linesource_strength!(sys::NavierStokes,q::Vector{T}) where {T<:Real}
    @unpack τ_bc = sys
    if isnothing(τ_bc)
      return sys
    else
      total_len = mapreduce(τl -> length(τl.τ),+,τ_bc)
      @assert total_len == length(q)
      qcpy = copy(q)
      first_index = 0
      for τl in τ_bc
        τl.τ .= view(q,first_index+1:first_index+length(τl.τ))
        first_index += length(τl.τ)
      end
    end
    return sys
end

include("navierstokes/surfacevelocities.jl")
include("navierstokes/fields.jl")
include("navierstokes/pointforce.jl")
include("navierstokes/basicoperators.jl")
include("navierstokes/rigidbodyoperators.jl")
include("navierstokes/movingbodyoperators.jl")
include("navierstokes/timemarching.jl")


text:
using Gadfly
using Colors
using MLLibrary

path=joinpath(Pkg.dir("MLLibrary"),"resources","datasets")
seed=123
threshold =.1
delimiter = ","
QRDecomposeMethod =MLLibrary.LinearAlgebra.HouseHolders.householdDecompose

#resources = ["3percent-miscategorization","10percent-miscatergorization","abalone.data","Seperable.csv"]
resources=["res"=>"abalone.data" ,"proxy"=>[Dict("column"=>1,"proxy"=>["M" "F" "I"])],"path"=>Pkg.dir("MLLibrary","resources","datasets")]
errors = Dict()

training, test = MLLibrary.DataSampling.ThresholdSampling.getSets(joinpath(resources["path"],resources["res"]),",",threshold,seed)

proxy = resources["proxy"]

training=[MLLibrary.DataManipulation.DataTransform.addProxyColumns(training[:,1:end-1],proxy) training[:,end]]

test=[MLLibrary.DataManipulation.DataTransform.addProxyColumns(test[:,1:end-1],proxy) test[:,end]]

# println(training)
# println(test)

unscaled_training = MLLibrary.DataManipulation.DataTransform.convertToFloat(training)
unscaled_test= MLLibrary.DataManipulation.DataTransform.convertToFloat(test)

training_y = unscaled_training[:,end]
test_y = unscaled_test[:,end]
training = unscaled_training[:,1:end-1]
test = unscaled_test[:,1:end-1]
println("Data is ready")

unscaled_training = copy(training)
unscaled_test = copy(test)
training =MLLibrary.DataManipulation.DataTransform.zScale(training)

test = MLLibrary.DataManipulation.DataTransform.zScale(test,
MLLibrary.Statistics.meanColumn(training),MLLibrary.Statistics.stdColumn(training))

outfile = open(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week4","Problem1centroids.txt"),"w")
outfile3 = open(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week4","Problem1wcss.txt"),"w")
outfile4 = open(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week4","Problem1stats.txt"),"w")

outfile2 = open(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week4","Problem2.txt"),"w")

wcss=Dict()
rmse=Dict()
for k in [1,2,4,8,16]
  centroids,ktrain_y = MLLibrary.LearningAlgorithms.KMeans.trainKmeans(training,k,seed)

  write(outfile, "K: ","$k","\n\n")
  write(outfile, "Centroids:","\n","$centroids","\n\n")

  #this section is purely for average/standarddeviation. Can actually be removed...
  for i=1:k

    newTraining = training[find(x->x==i,ktrain_y),:]
    average = MLLibrary.Statistics.meanColumn(newTraining)
    standarddeviation = MLLibrary.Statistics.stdColumn(newTraining)

    write(outfile4, "K: ","$k","\n\n")

    write(outfile4, "Cluster $i","\n")
    write(outfile4,"Mean:","$average" ,"\n\n","Std:", "$standarddeviation","\n\n"
    )

  end
  #caclulate the wcss
  wcss[k]=MLLibrary.Utils.ErrorAnalysis.calculateWCSS(centroids,training,ktrain_y)
  wcsstemp = wcss[k]
  write(outfile3, "K: ","$k","\n\n")
  write(outfile3, "WCSS:","\n","$wcsstemp","\n\n")

  #Problem 2, we now run on test set
  kmeans_y = MLLibrary.LearningAlgorithms.KMeans.runKMeans(centroids,test,k)

  predicted_y=zeros(Float64,size(test_y))

  #go through each cluster, clean cluster, QR train, and predict
  for j=1:k
    training_indices=find(x->x==j,ktrain_y)
    test_indices=find(x->x==j,kmeans_y)
    println(k,' ',j, ' ',size(training_indices))
    println(k,' ',j, ' ',size(test_indices))
    if !(0 in size(test_indices))
      #println(k,j)
      cluster_k= unscaled_training[training_indices,:];
      cluster_k_y = training_y[training_indices];
      test_cluster = unscaled_test[test_indices,:];
      # println(size(test_cluster))

      standarddeviation = MLLibrary.Statistics.stdColumn(cluster_k)
      latentvar = []
      for l=1:length(standarddeviation)
        if standarddeviation[l]>.001
          push!(latentvar,l)
        end
      end
      cluster_k=cluster_k[:,latentvar]
      test_cluster= test_cluster[:,latentvar]
      betacol = MLLibrary.LearningAlgorithms.LinearRegression.linearRegression(cluster_k,cluster_k_y,QRDecomposeMethod)
      predicted_y[test_indices,:] = test_cluster*betacol
    end
  end
  write(outfile2, "RMSES", "\n\n")
  write(outfile2, "K: ","$k","\n\n")
  rmse[k]= MLLibrary.Utils.ErrorAnalysis.calculateRootMeanStandardError(predicted_y,test_y);
  rmsetemp = rmse[k]
  write(outfile2, "RMSE:","\n","$rmsetemp","\n\n")

  # println(predicted_y,test_y)

end
close(outfile)
close(outfile2)
close(outfile3)
close(outfile4)

println(rmse);
p = plot(

  Guide.ylabel("WCSS"),
  Guide.xlabel("K"),
  Guide.title("Within Cluster Sum of Squares"),
  x=collect(keys(wcss)),
  y=collect(values(wcss)),
  Geom.line
  # Geom.label(position=:above),
  # label=map(string,(values(wcss)))
  #color="DataSet"
)

draw(PNG(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week4","WCSSErrors.png"), 800px, 600px), p)

p = plot(

  Guide.ylabel("RMSE"),
  Guide.xlabel("K"),
  Guide.title("Root Mean Square Error"),
  x=collect(keys(rmse)),
  y=collect(values(rmse)),
  Geom.line,
  Geom.label(position=:above),
  label=map(string,(values(rmse)))
  #color="DataSet"
)
draw(PNG(joinpath(Pkg.dir("MLLibrary"),"resources","paResources","Week4","RMSEErrors.png"), 800px, 600px), p)

max_stars_repo_name: ajermaky/CustomMLLibrary

max_issues_repo_name: ajermaky/CustomMLLibrary

max_forks_repo_name: ajermaky/CustomMLLibrary


text:
function vector_to_symmetric_matrix_converter()
    E1 = [
        1.0 0.0
        0.0 0.0
        0.0 1.0
    ]
    E2 = [
        0.0 0.0
        0.0 1.0
        1.0 0.0
    ]
    return [E1, E2]
end

function make_row_matrix(matrix, vals)
    return hcat([v * matrix for v in vals]...)
end

function interpolation_matrix(vals, ndofs)
    return make_row_matrix(diagm(ones(ndofs)), vals)
end

function transform_gradient(gradf, jacobian)
    return gradf / Diagonal(jacobian)
end

function displacement_bilinear_form(
    basis,
    quad,
    stiffness,
    jacobian,
    dim,
    vectosymmconverter,
)
    nf = number_of_basis_functions(basis)
    ndofs = dim * nf
    matrix = zeros(ndofs, ndofs)
    detjac = prod(jacobian)
    for (p, w) in quad
        grad = transform_gradient(gradient(basis, p), jacobian)
        NK = sum([
            make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim
        ])
        matrix .+= NK' * stiffness * NK * detjac * w
    end
    return matrix
end

function mass_matrix(basis, quad, ndofs, scale)
    nf = number_of_basis_functions(basis)
    totaldofs = nf * ndofs
    matrix = zeros(totaldofs, totaldofs)
    for (p, w) in quad
        vals = basis(p)

        NI = interpolation_matrix(vals, ndofs)
        matrix .+= NI' * NI * scale * w
    end
    return matrix
end

function mass_matrix(
    basis,
    quad1,
    quad2,
    ndofs,
    facescale;
    weight_tolerance = 1e-3,
)
    numqp = length(quad1)
    @assert length(quad2) == length(facescale) == numqp
    nf = number_of_basis_functions(basis)
    totaldofs = ndofs * nf
    matrix = zeros(totaldofs, totaldofs)
    for qpidx = 1:numqp
        p1, w1 = quad1[qpidx]
        p2, w2 = quad2[qpidx]
        @assert isapprox(w1, w2, atol = weight_tolerance)

        vals1 = basis(p1)
        vals2 = basis(p2)

        NI1 = interpolation_matrix(vals1, ndofs)
        NI2 = interpolation_matrix(vals2, ndofs)

        matrix .+= NI1' * NI2 * facescale[qpidx] * w1
    end
    return matrix
end

function component_mass_matrix(
    basis,
    quad1,
    quad2,
    components,
    ndofs,
    facescale;
    weight_tolerance = 1e-3,
)
    numqp = length(quad1)
    @assert length(quad2) == length(facescale) == numqp
    @assert size(components) == (ndofs, numqp)
    nf = number_of_basis_functions(basis)
    totaldofs = ndofs * nf
    matrix = zeros(totaldofs, totaldofs)
    for qpidx = 1:numqp
        p1, w1 = quad1[qpidx]
        p2, w2 = quad2[qpidx]
        @assert isapprox(w1, w2, atol = weight_tolerance)

        component = components[:, qpidx]
        projector = component * component'

        vals1 = basis(p1)
        vals2 = basis(p2)

        NI1 = interpolation_matrix(vals1, ndofs)
        NI2 = make_row_matrix(projector, vals2)

        matrix .+= NI1' * NI2 * facescale[qpidx] * w1
    end
    return matrix
end

function surface_traction_component_operator(
    basis,
    quad1,
    quad2,
    components,
    normals,
    stiffness,
    dim,
    scalearea,
    jac,
    vectosymmconverter;
    weight_tolerance = 1e-3,
)
    numqp = length(quad1)
    @assert length(quad2) ==
            size(normals)[2] ==
            length(scalearea) ==
            size(components)[2] ==
            numqp
    @assert size(normals)[1] == size(components)[1] == dim

    nf = number_of_basis_functions(basis)
    totalndofs = dim * nf
    matrix = zeros(totalndofs, totalndofs)

    for qpidx = 1:numqp
        p1, w1 = quad1[qpidx]
        p2, w2 = quad2[qpidx]
        @assert isapprox(w1, w2, atol = weight_tolerance)

        vals = basis(p1)
        grad = transform_gradient(gradient(basis, p2), jac)

        normal = normals[:, qpidx]
        component = components[:, qpidx]
        projector = component * component'

        NK = sum([
            make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim
        ])
        N = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])
        NI = interpolation_matrix(vals, dim)

        matrix .+= NI' * projector * N * stiffness * NK * scalearea[qpidx] * w1
    end
    return matrix
end

function surface_traction_operator(
    basis,
    quad1,
    quad2,
    normals,
    stiffness,
    dim,
    scalearea,
    jac,
    vectosymmconverter;
    weight_tolerance = 1e-3,
)

    numqp = length(quad1)
    @assert length(quad2) == size(normals)[2] == length(scalearea) == numqp
    @assert size(normals)[1] == dim

    nf = number_of_basis_functions(basis)
    totalndofs = dim * nf
    matrix = zeros(totalndofs, totalndofs)

    for qpidx = 1:numqp
        p1, w1 = quad1[qpidx]
        p2, w2 = quad2[qpidx]
        @assert isapprox(w1, w2, atol = weight_tolerance)

        vals = basis(p1)
        grad = transform_gradient(gradient(basis, p2), jac)
        normal = normals[:, qpidx]
        NK = sum([
            make_row_matrix(vectosymmconverter[k], grad[:, k]) for k = 1:dim
        ])
        N = sum([normal[k] * vectosymmconverter[k]' for k = 1:dim])
        NI = interpolation_matrix(vals, dim)

        matrix .+= NI' * N * stiffness * NK * scalearea[qpidx] * w1
    end
    return matrix
end

function linear_form(rhsfunc, basis, quad, cellmap, ndofs, detjac)
    nf = number_of_basis_functions(basis)
    rhs = zeros(ndofs * nf)
    for (p, w) in quad
        vals = rhsfunc(cellmap(p))
        @assert length(vals) == ndofs
        N = interpolation_matrix(basis(p), ndofs)
        rhs .+= N' * vals * detjac * w
    end
    return rhs
end

function component_linear_form(
    rhsfunc,
    basis,
    quad,
    component,
    cellmap,
    ndofs,
    detjac,
)
    @assert length(component) == ndofs
    nf = number_of_basis_functions(basis)
    rhs = zeros(ndofs * nf)
    for (p, w) in quad
        vals = rhsfunc(cellmap(p))
        N = interpolation_matrix(basis(p), ndofs)
        rhs .+= vals * N' * component * detjac * w
    end
    return rhs
end


text:
using Ipopt
using JuMP
using Juniper
using MINLPTests
using Test

const IPOPT =
    JuMP.optimizer_with_attributes(Ipopt.Optimizer, "print_level" => 0)

const JUNIPER = JuMP.optimizer_with_attributes() do
    return Juniper.Optimizer(nl_solver = IPOPT, log_levels = [])
end

const NLP_SOLVERS = [IPOPT]
const MINLP_SOLVERS = [JUNIPER]
const POLY_SOLVERS = [IPOPT]
const MIPOLY_SOLVERS = [JUNIPER]

@testset "JuMP Model Tests" begin
    @testset "$(solver): nlp" for solver in NLP_SOLVERS
        MINLPTests.test_nlp(
            solver,
            exclude = ["005_011"],  # Uses the function `\`
            debug = true,
        )
        MINLPTests.test_directory("nlp", solver, include = ["001_010"])
        MINLPTests.test_nlp_cvx(solver)
    end
    @testset "$(solver): nlp_mi" for solver in MINLP_SOLVERS
        MINLPTests.test_nlp_mi(solver, exclude = [
            "005_011",  # Uses the function `\`
            "006_010",  # Bug in Juniper - handling of user-defined functions.
        ])
        MINLPTests.test_nlp_mi_cvx(solver)
    end
    @testset "$(solver): poly" for solver in POLY_SOLVERS
        MINLPTests.test_poly(solver)
        MINLPTests.test_poly_cvx(solver)
    end
    @testset "$(solver): poly_mi" for solver in MIPOLY_SOLVERS
        MINLPTests.test_poly_mi(solver)
        MINLPTests.test_poly_mi_cvx(solver)
    end
end


text:
module integrands

constant(x, c=1, varargs...) = c

linear(x, a=1, b=0) = mean(a.*x .+ b)

sphere(x, varargs...) = vecdot(x, x) < 1.0 ? 1.0 : 0.0

sinus(x, varargs...) = sum(sin.(x)) / length(x)

sinusc(x, varargs...) = sum(sinc.(x)) / length(x)

oscillatory(x, c=1, r=1) = cos(2.0*pi*r + sum(c.*x))

productpeak(x, c=1, w=0) = 1 / prod(c.^(-2.0) .+ (x .- w).^2)

cornerpeak(x, c=1, r=1) = 1 / (1 + sum(abs.(c.*x)))^(length(x)+r)

gaussian(x, c=1, w=0) = exp(- sum(c.^2 .* (x.-w).^2))

continuous(x, c=1, w=0) = exp(- sum(c.* abs.(x.-w)))

discontinuous(x, c=1, w=0) = all(x.<=w) ? exp(sum(c.*x)) : 0.0

end

max_stars_repo_name: florianLepretre/ucbature

max_issues_repo_name: florianLepretre/ucbature

max_forks_repo_name: florianLepretre/ucbature


text:
using Base.Test

# test imwarp()
img = reshape(float(collect(1:121).%2), 11, 11) # 11x11 checkerboard
warp = img

tform = [1 0 0;
        0 1 0;
        0 0 1]
offset = [0, 0]
img_warped, warped_offset = imwarp(img, tform, offset)
@test_approx_eq warp img_warped
@test warped_offset == [0, 0]

tform = [1 0 0;
        0 1 0;
        0 0 1]
offset = [10, 20]
img_warped, warped_offset = imwarp(img, tform, offset)
@test_approx_eq warp img_warped
@test warped_offset == [10, 20]

tform = [1 0 0;
        0 1 0;
        10 10 1]
offset = [0, 0]
img_warped, warped_offset = imwarp(img, tform, offset)
@test_approx_eq warp img_warped
@test warped_offset == [10, 10]

tform = [cos(0.5) -sin(0.5) 0;
        sin(0.5) cos(0.5) 0;
        0 0 1]
offset = [0.0, 0.0]
img_warped, warped_offset = imwarp(img, tform, offset)
@test warped_offset == [0, -6]


text:
using MuJoCo

############################################################# mujoco startup

pm = mj_loadXML(modelfile, C_NULL)

@test pm != nothing

pd = mj_makeData(pm)

@test typeof(pd) == Ptr{mjData}

m = unsafe_load(pm)
d = unsafe_load(pd)

@test m.nq == 28
@test m.nv == 27
@test m.nbody == 14
@test m.opt.timestep == 0.002


text:
@reexport module SensitivityTests
    """ An abstract type indicating a sensitivity test. """
    abstract type SensitivityTest end

    export SensitivityTest
end

"""
    SensitivityTests

A module defining sensitivity tests.
"""
SensitivityTests


text:
"""Calculate the spectral range of a Hamiltonian `H` on the real axis.

```julia
E_min, E_max = specrange(H; method=:auto, kwargs...)
```

calculates the approximate lowest and highest eigenvalues of `H`. Any imaginary
part in the eigenvalues is ignored: the routine is intended for (although not
strictly limited to) a Hermitian `H`.

This delegates to

```julia
specrange(H, method; kwargs...)
```

for the different methods.

The default `method=:auto` chooses the best method for the given `H`. This is
`:diag` for small matrices, and `:arnoldi` otherwise. Keyword arguments not
relevant to the underlying implementation will be ignored.
"""
function specrange(H; method=Val(:auto), kwargs...)
    return specrange(H, method; kwargs...)
end

function specrange(H, method::Symbol; kwargs...)
    return specrange(H, Val(method); kwargs...)
end

function specrange(H, method::Val{:auto}; kwargs...)
    if size(H)[1] <= 32
        # TODO: benchmark what a good cross-over point from exact
        # diagonalization to Arnoldi is
        return specrange(H, Val(:diag); kwargs...)
    end
    return specrange(H, Val(:arnoldi); kwargs...)
end

"""
```julia
E_min, E_max = specrange(H, :arnoldi; state=random_state(H), m_min=20,
                         m_max=60, prec=1e-3, nom_min=1e-15, enlarge=true)
```

uses [Arnoldi iteration](https://en.wikipedia.org/wiki/Arnoldi_iteration) with
`state` as the starting vector. It approximates the eigenvalues of `H` with
between `m_min` and `m_max` Ritz values, until the lowest and highest
eigenvalue are stable to a releative precision of `prec`. The `norm_min`
parameter is passed to the underlying [`arnoldi!`](@ref).

If `enlarge=true` (default) the returned `E_min` and `E_max` will be enlarged
via a heuristic to slightly over-estimate the spectral radius instead of
under-estimating it.
"""
function specrange(H, method::Val{:arnoldi}; kwargs...)

    state = get(kwargs, :state, random_state(H))
    m_max = get(kwargs, :m_max, 60)
    m_min = max(5, min(get(kwargs, :m_min, 25), m_max-1))
    prec = get(kwargs, :prec, 1e-3)
    norm_min = get(kwargs, :norm_min, 1e-15)
    enlarge = get(kwargs, :enlarge, true)

    R = ritzvals(H, state, m_min, m_max; prec=prec, norm_min=norm_min)
    E_min = real(R[1])
    E_max = real(R[end])
    if enlarge
        # We want to overestimate the spec.rad., not underestimate it, so we
        # use the distance to the next eigenvalues as a buffer
        E_min = 2 * E_min - real(R[2])
        E_max = 2 * E_max - real(R[end-1])
    end
    return E_min,  E_max
end

"""
```julia
E_min, E_max = specrange(H, :diag)
```

uses exact diagonization via the standard `eigenvals` function to obtain the
smallest and largest eigenvalue. This should only be used for relatively small
matrices.
"""
function specrange(H, method::Val{:diag}; kwargs...)
    evals = eigenvals(Array(H))
    return evals[1], evals[end]
end

"""Random normalized quantum state.

```julia
    Ψ = random_state(H)
```

returns a random normalized state compatible with the Hamiltonian `H`. This is
intended to provide a starting vector for estimating the spectral radius of `H`
via an Arnoldi method.
"""
function random_state(H::AbstractMatrix)
    N = size(H)[2]
    Ψ = rand(N) .* exp.((2π * im) .* rand(N))
    Ψ ./= norm(Ψ)
    return Ψ
end

"""Calculate a vector for Ritz values converged to a given precision.

```julia
R = ritzvals(G, state, m_min, m_max=2*m_min; prec=1e-5, norm_min=1e-15)
```

calculates a complex vector `R` of at least `m_min` and at most `m_max` Ritz
values.
"""
function ritzvals(G, state, m_min, m_max=2*m_min; prec=1e-5, norm_min=1e-15)
    if m_max <= m_min
        throw(ArgumentError("m_max=$m_max must be smaller than m_min=$min"))
    end
    m = max(5, min(m_min, m_max-1))

    Hess = Matrix{ComplexF64}(undef, m_max, m_max)
    q = [similar(state) for _ in 1:(m_max+1)]

    # Get the Ritz eigenvalues for order m-1, so that we can decide whether
    # the values for order m are already precise enough
    m₀ = m - 1
    m₀ = arnoldi!(Hess, q, m₀, state, G; extended=false, norm_min=norm_min)
    @assert m₀ > 1
    eigenvals = diagonalize_hessenberg_matrix(Hess, m)
    v̲r₀ = minimum(real(eigenvals))
    v̄r₀ = maximum(real(eigenvals))
    v̄i₀ = maximum(abs.(imag(eigenvals)))
    if m₀ == m - 1
        # Starting from original order m, increase `m` until the Ritz values
        # reach the desired precision
        extend_arnoldi!(Hess, q, m, G; norm_min=norm_min)
        eigenvals = diagonalize_hessenberg_matrix(Hess, m)
        v̲r = minimum(real(eigenvals))
        v̄r = maximum(real(eigenvals))
        v̄i = maximum(abs.(imag(eigenvals)))
        e̲r = (v̲r₀ ≠ 0.0) ? abs(1.0 - v̲r/v̲r₀) : 0.0
        ēr = (v̄r₀ ≠ 0.0) ? abs(1.0 - v̄r/v̄r₀) : 0.0
        ēi = (v̄i₀ ≠ 0.0) ? abs(1.0 - v̄i/v̄i₀) : 0.0
        while (e̲r > prec) || (ēr > prec) || ((v̄i₀ > 1e-14) && ēi > prec)
            v̲r₀ = v̲r
            v̄r₀ = v̄r
            v̄i₀ = v̄i
            m₀ = m
            m = m + 1
            extend_arnoldi!(Hess, q, m, G; norm_min=norm_min)
            (m == m₀) && break  # dimensionality exhausted
            eigenvals = diagonalize_hessenberg_matrix(Hess, m)
            v̲r = minimum(real(eigenvals))
            v̄r = maximum(real(eigenvals))
            v̄i = maximum(abs.(imag(eigenvals)))
            e̲r = abs(1.0 - (v̲r/v̲r₀))
            ēr = abs(1.0 - (v̄r/v̄r₀))
            ēi = abs(1.0 - (v̄i/v̄i₀))
            if (m == m_max)
                @warn "Ritz values did not converge within m_max=$m_max"
                break
            end
        end
    end
    return eigenvals
end


text:
# Note that this script can accept some limited command-line arguments, run
# `julia build_tarballs.jl --help` to see a usage message.
using BinaryBuilder, Pkg

name = "Bullet_Physics_SDK"
version = v"3.22.0"

# Collection of sources required to complete build
sources = [
    ArchiveSource("https://github.com/bulletphysics/bullet3/archive/refs/tags/3.22b.tar.gz",
                  "c6cd89ecbc4bd73fee64723c831c1578daab056d88774755a6f56afc6f417b2b")
]

# Bash recipe for building across all platforms
script = raw"""
cd $WORKSPACE/srcdir/bullet*
mkdir build && cd build
cmake -DCMAKE_INSTALL_PREFIX=$prefix \
    -DCMAKE_TOOLCHAIN_FILE=${CMAKE_TARGET_TOOLCHAIN} \
    -DCMAKE_BUILD_TYPE=Release \
    -DCMAKE_CXX_STANDARD="11" \
    -DCMAKE_POSITION_INDEPENDENT_CODE=ON \
    -DBUILD_PYBULLET=OFF \
    -DBUILD_PYBULLET_NUMPY=OFF \
    -DBUILD_UNIT_TESTS=OFF \
    -DBUILD_CPU_DEMOS=OFF \
    -DBUILD_OPENGL3_DEMOS=OFF \
    -DBUILD_BULLET2_DEMOS=OFF \
    -DBT_USE_EGL=OFF \
    -DUSE_DOUBLE_PRECISION=OFF \
    -DBUILD_EXTRAS=OFF \
    -S .. \
    -B .
make -j${nproc}
make install
"""

# These are the platforms we will build for by default, unless further
# platforms are passed in on the command line
platforms = expand_cxxstring_abis(supported_platforms())

# The products that we will ensure are always built
products = [
        FileProduct("lib/libBullet3OpenCL_clew.a", :libBullet3OpenCL_clew),
        FileProduct("lib/libBullet2FileLoader.a", :libBullet2FileLoader),
        FileProduct("lib/libBullet3Dynamics.a", :libBullet3Dynamics),
        FileProduct("lib/libBullet3Collision.a", :libBullet3Collision),
        FileProduct("lib/libBullet3Geometry.a", :libBullet3Geometry),
        FileProduct("lib/libBulletInverseDynamics.a", :libBulletInverseDynamics),
        FileProduct("lib/libBulletSoftBody.a", :libBulletSoftBody),
        FileProduct("lib/libBulletCollision.a", :libBulletCollision),
        FileProduct("lib/libBulletDynamics.a", :libBulletDynamics),
        FileProduct("lib/libLinearMath.a", :libLinearMath),
        FileProduct("lib/libBullet3Common.a", :libBullet3Common),
]

# Dependencies that must be installed before this package can be built
dependencies = Dependency[
]

# Build the tarballs, and possibly a `build.jl` as well.
build_tarballs(ARGS, name, version, sources, script, platforms, products, dependencies; julia_compat="1.6", preferred_gcc_version=v"5")


text:
using BinaryBuilder

sources = [
    "https://github.com/staticfloat/small_bin/raw/d846f4a966883e7cc032a84acf4fa36695d05482/broken_symlink/broken_symlink.tar.gz" =>
    "470d47f1e6719df286dade223605e0c7e78e2740e9f0ecbfa608997d52a00445",
]

script = """
mkdir -p a/b/c/foobar
echo "Entering $(pwd)/a/b/c/foobar, and setting FOOBAR=1"
cd a/b/c/foobar
FOOBAR=1
exit 1
"""

products(prefix) = [LibraryProduct(prefix, "libfoo", :libfoo)]
build_args = ["--verbose", "--debug", "x86_64-linux-gnu"]
build_tarballs(build_args, "Broken", sources, script, [Linux(:x86_64)], products, [])


text:
module DispatchOps

export
    # Nexus/simplenexus
    Graph, Digraph, Arc, ArcIter,
    add_arc!, rem_arc!, add_node!, rem_node!,
    nn, na, nodes, arcs, fadj, badj, has_arc, has_node,
    src, tgt, key,

    # Nexus/metanexus
    MetaGraph, MetaDigraph, props, get_prop, has_prop,
    set_props!, set_prop!, rem_prop!, clear_props!,
    filter_arcs, filter_nodes,

    # Sim/interface
    Libraries, Params, States, Accumulators, Simulation,

    # Sim/Events/plan
    buildGraph, hard_holdover_model, soft_holdover_model, optimizeModel, plan!,

    # Sim/Events/fulfill
    fulfill!,

    # Sim/Events/transport
    transport!,

    # Sim/Utils/flow
    reset!, initiate!, schedule!, run!,

    # Sim/Utils/demand
    noisify_fixed, noisify_varied,

    # Utils/analysis
    lost_sales, total_cost, fixed_cost, variable_cost

include("Nexus/Nexus.jl")
using .Nexus

include("Sim/Sim.jl")
using .Sim

end


text:
using Documenter, Modia3D

makedocs(
  #modules  = [Modia3D],
  sitename = "Modia3D",
  authors  = "Andrea Neumayr, Martin Otter, Gerhard Hippmann (DLR-SR)",
  format   = Documenter.HTML(prettyurls = false),
  pages    = [
     "Home"   => "index.md",
     "Tutorial" => [
        "tutorial/Tutorial.md",
        "tutorial/GettingStarted.md",
        "tutorial/HierarchicalModels.md",
        "tutorial/RecursiveModels.md",
        "tutorial/CollisionHandling.md"
     ],
     "Components" => [
        "Components/Object3D.md"
        "Components/Shapes.md"
        "Components/Joints.md"
        "Components/Materials.md"
        "Components/GravityField.md"
     ],
     "Internal" => [
        "internal/Profiling.md"
        "internal/DynamicDispatch.md"
        "internal/ContactDetection.md"
        "internal/ContactForceLaw.md"
     ],
  ]
)


text:
using Distributions
using IndependentHypothesisWeighting
using Random
using MLDataPattern
import MLDataPattern: kfolds

using StatsBase
using CategoricalArrays

Random.seed!(1)
m= 10000
kf = kfolds(MLDataPattern.shuffleobs(1:m), 5)
Xs = CategoricalVector(sample(1:2, m))
Ps = rand(BetaUniformMixtureModel(0.7, 0.2), m) .* (Xs.==1) .+ rand(Uniform(), m) .* (Xs.==2)

typeof(Xs)
typeof(Ps)

α = 0.1
sum(adjust(Ps, BenjaminiHochberg()) .<= α ) # 580 discoveries

ihw_grenander = IHW(weight_learner = GrenanderLearner(), α = α)
 # 677 discoveries

ihw_grenander_fit1 = fit(ihw_grenander, Ps, Xs, kf)
sum(adjust(ihw_grenander_fit1) .<= α)

Random.seed!(1)
ihw_grenander_fit2 = fit(ihw_grenander, Ps, Xs)
sum(adjust(ihw_grenander_fit2) .<= α)


text:
using MultiResolutionIterators
using Test

struct Document{T}
    content::T
end

Base.iterate(doc::Document)=iterate(doc.content)
Base.iterate(doc::Document, state)=iterate(doc.content, state)

Base.length(doc::Document)=Base.length(doc.content)

function MultiResolutionIterators.apply(ff, aiter::Document)
    Document(MultiResolutionIterators.apply(ff, aiter.content))
end

animal_info = [
    Document([["Turtles", "are", "reptiles", "."],
     ["They", "have", "shells", "."],
     ["They", "live", "in", "the", "water", "."]]),
    Document([["Cats", "are", "mammals", "."],
     ["They", "live", "on", "the", "internet", "."]])
    ]

# Declare an indexer.
struct AnimalTextIndexer end;

# Overload `levelname_map` this so it knows the name mapping
MultiResolutionIterators.levelname_map(::AnimalTextIndexer) = [
    :documents=>1,
    :sentences=>2,
    :words=>3, :tokens=>3, # can have multiple aliases for same level
    :characters=>4 # As characters themselves are not iterable this name/level has little effect
]
indexer = AnimalTextIndexer();

@testset "Consolidate" begin
    @test consolidate(animal_info[1]) |> typeof <: Document
    @test consolidate_levels(animal_info[1],1) |> typeof <: Document
    @test full_consolidate(animal_info[1]) |> typeof <: Document
end

@testset "Demo tests" begin
    # Merge all sentences
    docs_of_words1 =  full_consolidate(flatten_levels(animal_info, lvls(indexer, :sentences)))
    @test typeof(docs_of_words1) == Vector{Document{Vector{String}}}

    # I.e keep documents, and words
    docs_of_words2 = full_consolidate(flatten_levels(animal_info, (!lvls)(indexer, :words, :documents)))
    @test typeof(docs_of_words2) == Vector{Document{Vector{String}}}
    @test collect(docs_of_words2[1]) == collect(docs_of_words1[1])
    @test collect(docs_of_words2[2]) == collect(docs_of_words1[2])
end


text:
using DataFrames, Dash, DashHtmlComponents, DashCoreComponents, PlotlyJS, UrlDownload

df1 = DataFrame(urldownload("https://raw.githubusercontent.com/plotly/datasets/master/gapminderDataFiveYear.csv"))

years = unique(df1[!, :year])

app = dash()

app.layout = html_div() do
    dcc_graph(id = "graph"),
    dcc_slider(
        id = "year-slider-1",
        min = minimum(years),
        max = maximum(years),
        marks = Dict([Symbol(v) => Symbol(v) for v in years]),
        value = minimum(years),
        step = nothing,
    )
end

callback!(
    app,
    Output("graph", "figure"),
    Input("year-slider-1", "value"),
) do selected_year
    return Plot(
        df1[df1.year .== selected_year, :],
        Layout(
            xaxis_type = "log",
            xaxis_title = "GDP Per Capita",
            yaxis_title = "Life Expectancy",
            legend_x = 0,
            legend_y = 1,
            hovermode = "closest",
            transition_duration = 500
        ),
        x = :gdpPercap,
        y = :lifeExp,
        text = :country,
        group = :continent,
        mode = "markers",
        marker_size = 15,
        marker_line_color = "white",
    )
end

run_server(app, "0.0.0.0", debug = true)

max_stars_repo_name: joelostblom/dash-docs

max_issues_repo_name: joelostblom/dash-docs

max_forks_repo_name: joelostblom/dash-docs


text:
using Base.Test

using Celeste: Model, Transform, SensitiveFloats, MCMC, Synthetic
using StatsBase

# helper to create synthetic data
function generate_single_star_data(; lnr=7.7251)
    # generate images (data)
    ea, vp, catalog = SampleData.gen_sample_star_dataset()

    # adjust fluxes to match input lnr
    old_fluxes = catalog[1].star_fluxes
    _, colors = Model.fluxes_to_colors(old_fluxes)
    new_fluxes = Model.colors_to_fluxes(lnr, colors)
    catalog[1].star_fluxes = new_fluxes

    # cache fluxes location, etc
    ce0 = deepcopy(catalog[1])
    dat_images = deepcopy(ea.images)

    # make mod/dat images
    Synthetic.gen_images!(dat_images, catalog; expectation=false);

    # make sure ea images data reflects dat images for vb fitting
    for b in 1:5
      ea.images[b].pixels[:] = dat_images[b].pixels[:]
    end

    # save truedf for python plot comparison
    truedf = MCMC.catalog_to_data_frame_row(ce0)
    return truedf, ce0, dat_images, ea, vp
end

###############
#### tests ####
###############

function test_mcmc_catalog_to_data_frame_row()
    # create slightly less bright log re
    truedf, ce0, dat_images, ea, vp = generate_single_star_data(; lnr=6.)
end

function test_star_inference_functions()
    # gen data and initial params
    truedf, ce0, images, ea, vp = generate_single_star_data(; lnr=6.)
    patches = Model.get_sky_patches(images, [ce0]; radius_override_pix=25)
    background_images = [zeros(size(img.pixels)) for img in images]
    star_loglike, star_logprior, star_logpost, sample_star_prior, ra_lim, dec_lim, uniform_to_deg, deg_to_uniform =
        MCMC.make_star_inference_functions(images, ce0;
                                           patches=patches[1,:], background_images=background_images,
                                           pos_delta = [2., 2.])
    th_rand = sample_star_prior()
    @test !isnan(star_loglike(th_rand))
    @test !isnan(star_logprior(th_rand))
    @test !isnan(star_logpost(th_rand))
end

function test_gal_inference_functions()
    # gen data and initial params
    truedf, ce0, images, ea, vp = generate_single_star_data(; lnr=6.)
    patches = Model.get_sky_patches(images, [ce0]; radius_override_pix=25)
    background_images = [zeros(size(img.pixels)) for img in images]
    gal_loglike, gal_logprior, gal_logpost, sample_gal_prior, ra_lim, dec_lim, uniform_to_deg, deg_to_uniform =
        MCMC.make_gal_inference_functions(images, ce0;
                                           patches=patches[1,:], background_images=background_images,
                                           pos_delta = [2., 2.])
    th_rand = sample_gal_prior()
    @test !isnan(gal_loglike(th_rand))
    @test !isnan(gal_logprior(th_rand))
    @test !isnan(gal_logpost(th_rand))
end

@testset "MCMC" begin
    test_star_inference_functions()
    test_gal_inference_functions()
end


text:
function IS_DQN_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        term = 𝒟[:done] .| 𝒟[:fail]
        term .* Base.log.(𝒟[:fail] .+ 1f-16) .+ (1.f0 .- term) .* logsumexp(𝒫[:xlogprobs] .+ value(π, 𝒟[:sp]), dims=1)
        # abs.(𝒟[:r]) .* (𝒟[:r] .< -10)  .+ (1.f0 .- 𝒟[:done]) .* γ .* sum(exp.(𝒫[:xlogprobs]) .* value(π, 𝒟[:sp]), dims=1)
        # 𝒟[:done] .* 𝒟[:fail] .+ (1.f0 .- 𝒟[:done]) .* sum(exp.(𝒫[:xlogprobs]) .* value(π, 𝒟[:sp]), dims=1)
end

function IS_estimate_log_pfail(π, px, s, Nsamples=10)
        xs_and_logpdfs = [exploration(π, s) for _ in 1:Nsamples]
        xs = [e[1] for e in xs_and_logpdfs]
        νlogpdfs = vcat([e[2] for e in xs_and_logpdfs]...)
        logpfails = vcat([value(π, s, x) for x in xs]...)
        logpxs = Float32.(vcat([logpdf.(px, x) for x in xs]...))

        -Float32(Base.log(Nsamples)) .+ logsumexp(νlogpdfs .+ logpfails .- logpxs, dims=1)
end

function estimate_log_pfail(π, px, s, Nsamples=10)
        xs = [reshape(rand(px, size(s,2)), :, size(s,2)) for _ in 1:Nsamples]
        logpfails = vcat([value(π, s, x) for x in xs]...)

        -Float32(Base.log(Nsamples)) .+ logsumexp(logpfails, dims=1)
end

function estimate_pfail(π, px, s, Nsamples=10)
        xs = [reshape(rand(px, size(s,2)), :, size(s,2)) for _ in 1:Nsamples]
        pfails = vcat([value(π, s, x) for x in xs]...)

        mean(pfails, dims=1)
end

function IS_Continuous_target(π, 𝒫, 𝒟, γ::Float32; kwargs...)
        term = 𝒟[:done] .| 𝒟[:fail]
        # term .* (𝒟[:fail] .== false) .* -100f0 .+ (1.f0 .- term) .* estimate_log_pfail(π, 𝒫[:px], 𝒟[:sp], 𝒫[:N_IS_Samples])
        term .* 𝒟[:fail] .+ (1.f0 .- term) .* estimate_pfail(π, 𝒫[:px], 𝒟[:sp], 𝒫[:N_IS_Samples])
end

function IS_L_KL_log(π, 𝒫, 𝒟; kwargs...)
        x, logνx = exploration(π, 𝒟[:s])
        logpfail_s = Zygote.ignore() do
                estimate_log_pfail(π, 𝒫[:px], 𝒟[:s], 𝒫[:N_IS_Samples])
        end
        logpx = logpdf(𝒫[:px], x)
        logpfail_sx = value(π, 𝒟[:s], x)

        -mean(exp.( logpx .+ logpfail_sx .- logνx .- logpfail_s) .* logνx)
end

function IS_L_KL(π, 𝒫, 𝒟; kwargs...)
        x, logνx = exploration(π, 𝒟[:s])
        νx = exp.(logνx)
        pfail_s = Zygote.ignore() do
                estimate_pfail(π, 𝒫[:px], 𝒟[:s], 𝒫[:N_IS_Samples])
        end
        px = exp.(logpdf(𝒫[:px], x))
        pfail_sx = value(π, 𝒟[:s], x)

        -mean(logνx .* px .* pfail_sx ./ (νx .* pfail_s))
end

function compute_IS_weight(𝒟, 𝒫; info=Dict())
        𝒟[:weight] .= exp.(sum(𝒫[:xlogprobs] .* 𝒟[:x], dims=1) .- 𝒟[:xlogprob])
        info[:mean_is_weight] = mean(𝒟[:weight])
        info[:std_is_weight] = std(𝒟[:weight])
        info[:min_is_weight] = minimum(𝒟[:weight])
        info[:max_is_weight] = maximum(𝒟[:weight])
end

function compute_IS_weight_continuous(𝒟, 𝒫; info=Dict())
        if 𝒫[:px] isa UnivariateDistribution
                𝒟[:weight] .= clamp.(exp.(reshape(logpdf.(𝒫[:px], 𝒟[:x]), 1, :) .- 𝒟[:xlogprob]), 0f0, 5f0)
        else
                𝒟[:weight] .= clamp.(exp.(reshape(logpdf(𝒫[:px], 𝒟[:x]), 1, :) .- 𝒟[:xlogprob]), 0f0, 5f0)
        end
        info[:mean_is_weight] = mean(𝒟[:weight])
        info[:std_is_weight] = std(𝒟[:weight])
        info[:min_is_weight] = minimum(𝒟[:weight])
        info[:max_is_weight] = maximum(𝒟[:weight])
end

ISARL_DQN(;kwargs...) = DQN(;c_loss=td_loss(name=:x_Qavg, a_key=:x), prefix="x_", target_fn=IS_DQN_target)
ISARL_DDPG(;kwargs...) = DQN(;a_loss=IS_L_KL, c_loss=td_loss(name=:x_Qavg, a_key=:x), prefix="x_", target_fn=IS_Continuous_target, )

function CERL(;𝒮_pro,
               𝒮_ant,
               px,
               log::NamedTuple=(;),
               train_pro_every::Int=1,
               train_ant_every::Int=1,
               buffer_size=1000, # Size of the buffer
               required_columns=[:fail, :x, :xlogprob, :weight],
               buffer::ExperienceBuffer=ExperienceBuffer(𝒮_pro.S, 𝒮_pro.agent.space, buffer_size, required_columns), # The replay buffer
               buffer_init::Int = max(max(𝒮_pro.c_opt.batch_size, 𝒮_ant.c_opt.batch_size), 200) # Number of observations to initialize the buffer with
               )

               # Set buffer parameters to be consistent between solvers (Since buffer is shared)
       𝒮_pro.required_columns = required_columns
       𝒮_pro.buffer = buffer
       𝒮_pro.buffer_init = buffer_init
       𝒮_ant.required_columns = required_columns
       𝒮_ant.buffer = buffer
       𝒮_ant.buffer_init = buffer_init
        AdversarialOffPolicySolver(;𝒮_pro=𝒮_pro,
                                    𝒮_ant=𝒮_ant,
                                    px=px,
                                    train_pro_every=train_pro_every,
                                    train_ant_every=train_ant_every,
                                    log=LoggerParams(;dir="log/isarl", log...),)
end


text:
# Use baremodule to shave off a few KB from the serialized `.ji` file
baremodule CImGuiPack_jll
using Base
using Base: UUID
import JLLWrappers

JLLWrappers.@generate_main_file_header("CImGuiPack")
JLLWrappers.@generate_main_file("CImGuiPack", UUID("333409e9-af72-5310-9767-d6ad21a76a05"))
end  # module CImGuiPack_jll


text:
"""
$(SIGNATURES)

Checks if upcoming integration is a new bit for generic BOC signal
"""
function is_upcoming_integration_new_bit(
    boc::BOCcos,
    prns,
    num_prns
)
    is_upcoming_integration_new_bit(boc.system, prns, num_prns)
end

function get_default_correlator(
    ::BOCcos,
    numAnts::NumAnts
)
    EarlyPromptLateCorrelator(numAnts)
end


text:
export chess960fen

function insertatempty!(array, item, index)
    for i ∈ 1:length(array)
        if array[i] ≠ '?' && i ≤ index
            index += 1
        end
    end
    array[index] = item
end

KNIGHT_TABLE_960 =
    [[1, 1], [1, 2], [1, 3], [1, 4], [2, 2], [2, 3], [2, 4], [3, 3], [3, 4], [4, 4]]

function backrank(n::Int)
    result = fill('?', 8)
    (n2, b1) = fldmod(n, 4)
    result[2b1+2] = 'b'
    (n3, b2) = fldmod(n2, 4)
    result[2b2+1] = 'b'
    (n4, q) = fldmod(n3, 6)
    insertatempty!(result, 'q', q + 1)
    knights = KNIGHT_TABLE_960[n4+1]
    insertatempty!(result, 'n', knights[1])
    insertatempty!(result, 'n', knights[2])
    insertatempty!(result, 'r', 1)
    insertatempty!(result, 'k', 1)
    insertatempty!(result, 'r', 1)
    reduce(*, result)
end

function castlestring960(setup, buffer)
    ch1 = '?'
    ch2 = '?'
    for i ∈ 8:-1:1
        if setup[i] == 'r'
            ch1 = Char(i - 1 + Int('a'))
            break
        end
    end
    for i ∈ 1:8
        if setup[i] == 'r'
            ch2 = Char(i - 1 + Int('a'))
            break
        end
    end
    print(buffer, uppercase(ch1))
    print(buffer, uppercase(ch2))
    print(buffer, lowercase(ch1))
    print(buffer, lowercase(ch2))
end

"""
    chess950fen(i)

Returns the FEN string of the Chess960 position with index `i`.

The parameter `i` must be an integer in the range 0-959.
"""
function chess960fen(i)
    setup = backrank(i)
    result = IOBuffer()
    print(result, setup)
    print(result, "/pppppppp/8/8/8/8/PPPPPPPP/")
    print(result, uppercase(setup))
    print(result, " w ")
    castlestring960(setup, result)
    print(result, " -")
    String(take!(result))
end


text:
function DumpInformation(sys::System,comm::MPI.Comm,root::Int64)
    MPI.Barrier(comm)
    if sys.trjdump!=NoDump()
        coords_all=SerializeArray(sys.atoms.coords_all)
        MPI.Gatherv!(SerializeArray(sys.atoms.coords),VBuffer(coords_all,sys.mpi.commsizes*3),root,comm)
        DeSerializeAllCoord!(sys,coords_all)
        if MPI.Comm_rank(comm)==root
            DumpTrajectory(sys,sys.trjdump)
        end
    end
    MPI.Barrier(comm)
end

function DumpTrajectory(sys::System,dump::LAMMPSTrj)
    if sys.first_step[1] || sys.current_step[1]%dump.freq==0
        string_out="ITEM: TIMESTEP"
        string_out*="\n"
        string_out*=string(sys.current_step[1])
        string_out*="\n"
        string_out*="ITEM: NUMBER OF ATOMS"
        string_out*="\n"
        string_out*=string(sys.atoms.total_N[1])
        string_out*="\n"
        string_out*="ITEM: BOX BOUNDS pp pp pp"
        string_out*="\n"
        string_out*="0 "*string(sys.pbc.Lx)
        string_out*="\n"
        string_out*="0 "*string(sys.pbc.Ly)
        string_out*="\n"
        string_out*="0 "*string(sys.pbc.Lz)
        string_out*="\n"

        string_out*="ITEM: ATOMS"
        string_out*=" id"
        string_out*=" type x y z"

        if dump.is_velocity==true
            string_out*=" vx vy vz"
        end
        if dump.is_force==true
            string_out*=" fx fy fz"
        end
        string_out*="\n"
        write(dump.io,string_out)
        for atomii=1:sys.atoms.total_N[1]
            @fastmath @inbounds begin
            string_out=""
            string_out*=string(atomii)*" "
            string_out*=string(sys.atoms.types_all[atomii])*" "*string(sys.atoms.coords_all[atomii][1])*" "*string(sys.atoms.coords_all[atomii][2])*" "*string(sys.atoms.coords_all[atomii][3])*" "

            if dump.is_velocity==true
                string_out*=string(sys.atoms.vels[atomii])*" "*string(sys.atoms.vels[atomii])*" "*string(sys.atoms.vels[atomii])*" "
            end
            if dump.is_force==true
                string_out*=string(sys.atoms.forces[atomii])*" "*string(sys.atoms.forces[atomii])*" "*string(sys.atoms.forces[atomii])*" "
            end
            string_out*="\n"
            write(dump.io,string_out)
            end
        end
    end
end


text:
println("Starting @async")
for i in 1:10
    @async begin
        sleep(rand())
        println(i)
    end
end
println("Done @async")

sleep(2)

println("\n\nStarting @sync-@async")
@sync for i in 1:10
    @async begin
        sleep(rand())
        println(i)
    end
end
println("Done @sync-@async")